Merge branch 'development' of https://github.com/CIMM-ORG/SWMM5plus i…

…nto development

definitions/define_indexes.f90

@@ -50,8 +50,9 @@ module define_indexes
         enumerator :: li_P_image             ! image number assigned from BIPquick
         enumerator :: li_parent_link         ! A map to the corresponding SWMM link after a BIPquick link-split
         enumerator :: li_weir_EndContractions ! (0,1) to indicate contraction
-        enumerator :: li_RoadSurface           ! roadsurface type for roadway weir
+        enumerator :: li_RoadSurface         ! roadsurface type for roadway weir
         enumerator :: li_curve_id            ! curve id if the link is associated with any curve
+        enumerator :: li_lateralInflowNode   ! downstream node from which the lateral inflow is coming from
         enumerator :: li_first_elem_idx
         enumerator :: li_last_elem_idx
         enumerator :: li_transect_idx         ! transect index if the link is associated with an irregular geometry transect
@@ -100,6 +101,7 @@ module define_indexes
         enumerator :: lr_Setting               !% the 0 to 1 open/close setting of EPA-SWMM
         enumerator :: lr_TargetSetting         !% target setting of a control action
         enumerator :: lr_TimeLastSet           !% the time (in seconds) the link setting was last changed
+        enumerator :: lr_VolumeFractionMetric  !% a metric to distribute the lateral inflows across a link
         enumerator :: lr_ZbottomUp             ! Z bottom of upstream node
         enumerator :: lr_ZbottomDn             ! Z bottom of downstream node
         enumerator :: lr_lastplusone !% must be last enum item
@@ -195,6 +197,7 @@ module define_indexes
         enumerator :: nr_head
         enumerator :: nr_Volume
         enumerator :: nr_Flooding
+        enumerator :: nr_UpLinksFullVolume
         enumerator :: nr_JunctionBranch_Kfactor
         enumerator :: nr_lastplusone !% must be last enum item
     end enum
@@ -285,6 +288,7 @@ module define_indexes
          enumerator :: ei_geometryType              !% KEY cross-sectional geometry type  (static)
          enumerator :: ei_barrels                   !% Integer number of barrels
          enumerator :: ei_HeqType                   !% KEY type of head equation (static)
+         enumerator :: ei_lateralInflowNode         !% SWMM node from which the element will get lateral inflow
          enumerator :: ei_link_Gidx_SWMM            !% link index from global SWMM network  (static)
          enumerator :: ei_link_Gidx_BIPquick        !% link index from global BIPquick network  (static)
          enumerator :: ei_link_pos                  !% position (elem from upstream = 1 to downstream = n) in link
@@ -403,6 +407,7 @@ module define_indexes
         enumerator :: er_VolumePonded               !% volume ponding this time step
         enumerator :: er_VolumePondedTotal          !% Volume ponded total
         enumerator :: er_VolumeStore                !% temporary storage used for adjacent AC and ETM elements
+        enumerator :: er_VolumeFractionMetric       !% metric to distribute lateral inflow across channel/conduit elements
         enumerator :: er_WaveSpeed                  !% wave speed in element
         enumerator :: er_YoverYfull                 !% ratio of depth to full depth for tabular geometry
         enumerator :: er_Zbottom                    !% bottom elevation of element (static)

definitions/define_settings.f90

@@ -405,6 +405,7 @@ module define_settings
         ! real(8) :: MinElemLengthFactor = 0.5d0           !% define the minimum allowable fraction of an element size to help with the cfl
         !integer :: MinElemLengthMethod = ElemLengthAdjust
         logical :: UseNominalElemLength = .true.
+        logical :: DistributeOpenChannelInflowsTF = .false.
         real(8) :: NominalElemLength   = 10.0d0
         integer :: MinElementPerLink   = 3               !% force a minimum number of elements per link
         logical :: UseEquivalentOrifice = .false.        !% replace small conduits with equivalent orifice
@@ -1110,6 +1111,11 @@ subroutine define_settings_load()
         if (found) setting%Discretization%UseEquivalentOrifice = logical_value
         if ((.not. found) .and. (jsoncheck)) stop "Error - json file - setting " //'Discretization.UseEquivalentOrifice not found'
 
+        !%                       Discretization.DistributeOpenChannelInflowsTF
+        call json%get('Discretization.DistributeOpenChannelInflowsTF', logical_value, found)
+        if (found) setting%Discretization%DistributeOpenChannelInflowsTF = logical_value
+        if ((.not. found) .and. (jsoncheck)) stop "Error - json file - setting " //'Discretization.DistributeOpenChannelInflowsTF not found'
+
         !%                       Discretization.EquivalentOrificeDischargeCoeff
         call json%get('Discretization.EquivalentOrificeDischargeCoeff', real_value, found)
         if (found) setting%Discretization%EquivalentOrificeDischargeCoeff = real_value
@@ -1497,6 +1503,21 @@ subroutine define_settings_load()
         if (found) setting%Output%DataOut%isVolumeOut = logical_value
         if ((.not. found) .and. (jsoncheck)) stop "Error - json file - setting " // 'Output.DataOut.isVolumeOut not found'
 
+        !%                       Dataout.isVolumeConsOut
+        call json%get('Output.DataOut.isVolumeConsOut', logical_value, found)
+        if (found) setting%Output%DataOut%isVolumeConsOut = logical_value
+        if ((.not. found) .and. (jsoncheck)) stop "Error - json file - setting " // 'Output.DataOut.isVolumeConsOut not found'
+
+        !%                       Dataout.isVolumeOverflowOut
+        call json%get('Output.DataOut.isVolumeOverflowOut', logical_value, found)
+        if (found) setting%Output%DataOut%isVolumeOverflowOut = logical_value
+        if ((.not. found) .and. (jsoncheck)) stop "Error - json file - setting " // 'Output.DataOut.isVolumeOverflowOut not found'
+
+        !%                       Dataout.isVolumePondedOut
+        call json%get('Output.DataOut.isVolumePondedOut', logical_value, found)
+        if (found) setting%Output%DataOut%isVolumePondedOut = logical_value
+        if ((.not. found) .and. (jsoncheck)) stop "Error - json file - setting " // 'Output.DataOut.isVolumePondedOut not found'
+
         !%                       Dataout.isWaveSpeedOut
         call json%get('Output.DataOut.isWaveSpeedOut', logical_value, found)
         if (found) setting%Output%DataOut%isWaveSpeedOut = logical_value

initialization/initial_condition.f90

@@ -908,14 +908,17 @@ subroutine init_IC_get_flow_and_roughness_from_linkdata (thisLink)
         !%------------------------------------------------------------------  
         !% --- handle all the initial conditions that don't depend on geometry type
         where (elemI(:,ei_link_Gidx_BIPquick) == thisLink)
-            elemR(:,er_Flowrate)           = link%R(thisLink,lr_FlowrateInitial) / link%I(thisLink,li_barrels)
-            elemR(:,er_Flowrate_N0)        = link%R(thisLink,lr_FlowrateInitial) / link%I(thisLink,li_barrels)
-            elemR(:,er_Flowrate_N1)        = link%R(thisLink,lr_FlowrateInitial) / link%I(thisLink,li_barrels)
-            elemR(:,er_ManningsN)          = link%R(thisLink,lr_Roughness)
+            elemR(:,er_Flowrate)             = link%R(thisLink,lr_FlowrateInitial) / link%I(thisLink,li_barrels)
+            elemR(:,er_Flowrate_N0)          = link%R(thisLink,lr_FlowrateInitial) / link%I(thisLink,li_barrels)
+            elemR(:,er_Flowrate_N1)          = link%R(thisLink,lr_FlowrateInitial) / link%I(thisLink,li_barrels)
+            elemR(:,er_ManningsN)            = link%R(thisLink,lr_Roughness)
             !% --- distribute minor losses uniformly over all the elements in thi link
-            elemR(:,er_Kconduit_MinorLoss) = link%R(thisLink,lr_Kconduit_MinorLoss) / (real(lastelem - firstelem + oneI,8))
-            elemR(:,er_FlowrateLimit)      = link%R(thisLink,lr_FlowrateLimit)
-            elemR(:,er_SeepRate)           = link%R(thisLink,lr_SeepRate)
+            elemR(:,er_Kconduit_MinorLoss)   = link%R(thisLink,lr_Kconduit_MinorLoss) / (real(lastelem - firstelem + oneI,8))
+            !% --- distribute volume fraction for lateral inflow across elements
+            elemR(:,er_VolumeFractionMetric) = link%R(thisLink,lr_VolumeFractionMetric) * elemR(:,er_Length) / link%R(thisLink,lr_Length)
+            elemI(:,ei_lateralInflowNode)    = link%I(thisLink,li_lateralInflowNode)
+            elemR(:,er_FlowrateLimit)        = link%R(thisLink,lr_FlowrateLimit)
+            elemR(:,er_SeepRate)             = link%R(thisLink,lr_SeepRate)
         endwhere
 
         !% --- assign minor losses

initialization/initialization.f90

@@ -475,7 +475,7 @@ subroutine init_linknode_arrays()
         !%   EPA-SWMM functionalities.
         !%------------------------------------------------------------------
         !% Declarations   
-            integer          :: ii, total_n_links
+            integer          :: ii, jj, total_n_links, link_idx
             integer, pointer :: linkUp, linkDn
             logical          :: noerrorfound
             character(64)    :: subroutine_name = 'init_linknode_arrays'
@@ -573,7 +573,7 @@ subroutine init_linknode_arrays()
             node%I(link%I(ii,li_Mnode_d), ni_idx_base1 + node%I(link%I(ii,li_Mnode_d), ni_N_link_u)) = ii
 
             !% --- increment the connection counter for the node upstream
-            node%I(link%I(ii,li_Mnode_u), ni_N_link_d) = node%I(link%I(ii,li_Mnode_u), ni_N_link_d) + 1
+            node%I(link%I(ii,li_Mnode_u), ni_N_link_d) = node%I(link%I(ii,li_Mnode_u), ni_N_link_d) + oneI
 
             !% --- set the maps for the upstream node to the downstream link (ni_Mlink_d#)
             !%     NOTE: this makes use of the ordering of the ni_Mlink_d# in define_indexes
@@ -737,7 +737,7 @@ subroutine init_linknode_arrays()
                 call util_crashpoint(387666)
             end if
 
-            if (node%I(ii, ni_N_link_u) > max_dn_branch_per_node) then
+            if (node%I(ii, ni_N_link_d) > max_dn_branch_per_node) then
                 if (this_image() == 1) then
                     write(*,*) 'USER CONFIGURATION ERROR'
                     write(*,"(A,i4,A)") 'One or more nodes have more than ',max_dn_branch_per_node,' downstream connections'
@@ -1135,10 +1135,43 @@ subroutine init_linknode_arrays()
                 if (node%I(ii, ni_node_type) == nJ1) node%I(ii, ni_node_type) = nBCup
             end if
 
+
+            !% --- calculate the full volumes of the links upstream of a node
+            !%     this volume will be used later the distribute lateral inflows across links
+            node%R(ii,nr_UpLinksFullVolume) = zeroR
+
+            do jj = 1,node%I(ii,ni_N_link_u)
+                !% pointer to the upstream links of that node
+                link_idx = node%I(ii,ni_idx_base1 + jj)
+
+                !% only calculate the volumes of open channels 
+                if (link%I(link_idx,li_link_type) == lchannel) then
+                    !% add the the volume of the upstream links
+                    node%R(ii,nr_UpLinksFullVolume) = node%R(ii,nr_UpLinksFullVolume) &
+                                    + link%R(link_idx,lr_FullArea) * link%R(link_idx,lr_Length)    
+                end if
+
+            end do
+
             !% --- note pattern initialization MUST be called after inflows are set
             node%I(ii,ni_pattern_resolution) = interface_get_BC_resolution(ii)
         end do
 
+        !% --- find the volume fraction metric for channel links to distribute lateral inflows
+        do ii = 1, setting%SWMMinput%N_link
+            if (link%I(ii,li_link_type) == lchannel) then
+                link%R(ii,lr_VolumeFractionMetric) = link%R(ii,lr_FullArea) * link%R(ii,lr_Length) &
+                                                / node%R(link%I(ii,li_Mnode_d),nr_UpLinksFullVolume)
+            else
+                !% --- lateral inflows will not be distributed for non channel links.
+                link%R(ii,lr_VolumeFractionMetric) = zeroR
+            end if
+
+            !% save the node from which lateral inflows will be distributed
+            link%I(ii,li_lateralInflowNode)    = link%I(ii,li_Mnode_d)
+
+        end do
+
         !% --- error checking for disconnected nodes
         noerrorfound = .true.
         do ii = 1,N_node