LCOV - code coverage report
Current view: top level - mRM - mo_mrm_routing.f90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 53 58 91.4 %
Date: 2024-04-30 08:53:32 Functions: 2 2 100.0 %

          Line data    Source code
       1             : !> \file mo_mrm_routing.f90
       2             : !> \brief   \copybrief mo_mrm_routing
       3             : !> \details \copydetails mo_mrm_routing
       4             : 
       5             : !> \brief Performs runoff routing for mHM at level L11.
       6             : !> \details This module performs flood routing at a given time step through the stream network at level L11 to the sink cell.
       7             : !! The Muskingum flood routing algorithm is used.
       8             : !> \changelog
       9             : !! - Stephan Thober Aug 2015
      10             : !!   - adapted to mRM
      11             : !! - Sebastian Mueller Jun 2020
      12             : !!   - outsourcing helper functions
      13             : !> \authors Luis Samaniego
      14             : !> \date Dec 2012
      15             : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
      16             : !! mHM is released under the LGPLv3+ license \license_note
      17             : !> \ingroup f_mrm
      18             : MODULE mo_mrm_routing
      19             : 
      20             :   ! This module performs runoff flood routing for mHM.
      21             : 
      22             :   ! Written Luis Samaniego, Dec 2012
      23             : 
      24             :   USE mo_kind, ONLY : i4, dp
      25             : 
      26             :   IMPLICIT NONE
      27             : 
      28             :   PRIVATE
      29             : 
      30             :   PUBLIC :: mRM_routing
      31             : 
      32             :   ! ------------------------------------------------------------------
      33             : 
      34             : CONTAINS
      35             : 
      36             :   ! ------------------------------------------------------------------
      37             : 
      38             :   !    NAME
      39             :   !        mRM_routing
      40             : 
      41             :   !    PURPOSE
      42             :   !>       \brief route water given runoff
      43             : 
      44             :   !>       \details This routine first performs mpr for the routing variables
      45             :   !>       if required, then accumulates the runoff to the routing resolution
      46             :   !>       and eventually routes the water in a third step. The last step is
      47             :   !>       repeated multiple times if the routing timestep is smaller than
      48             :   !>       the timestep of the hydrological timestep
      49             : 
      50             :   !    INTENT(IN)
      51             :   !>       \param[in] "logical :: read_states"                            whether states are derived from restart file
      52             :   !>       \param[in] "integer(i4) :: processCase"                        Process switch for routing
      53             :   !>       \param[in] "real(dp), dimension(:) :: global_routing_param"    routing parameters
      54             :   !>       \param[in] "real(dp), dimension(:) :: L1_total_runoff"         total runoff from L1 grid cells
      55             :   !>       \param[in] "real(dp), dimension(:) :: L1_areaCell"             L1 cell area
      56             :   !>       \param[in] "integer(i4), dimension(:) :: L1_L11_Id"            L1 cell ids on L11
      57             :   !>       \param[in] "real(dp), dimension(:) :: L11_areaCell"            L11 cell area
      58             :   !>       \param[in] "integer(i4), dimension(:) :: L11_L1_Id"            L11 cell ids on L1
      59             :   !>       \param[in] "integer(i4), dimension(:) :: L11_netPerm"          L11 routing order
      60             :   !>       \param[in] "integer(i4), dimension(:) :: L11_fromN"            L11 source grid cell order
      61             :   !>       \param[in] "integer(i4), dimension(:) :: L11_toN"              L11 target grid cell order
      62             :   !>       \param[in] "integer(i4) :: L11_nOutlets"                       L11 number of outlets/sinks
      63             :   !>       \param[in] "integer(i4) :: timestep"                           simulation timestep in [h]
      64             :   !>       \param[in] "real(dp) :: tsRoutFactor"                          factor between routing timestep and
      65             :   !>       hydrological timestep
      66             :   !>       \param[in] "integer(i4) :: nNodes"                             number of nodes
      67             :   !>       \param[in] "integer(i4) :: nInflowGauges"                      number of inflow gauges
      68             :   !>       \param[in] "integer(i4), dimension(:) :: InflowGaugeIndexList" index list of inflow gauges
      69             :   !>       \param[in] "logical, dimension(:) :: InflowGaugeHeadwater"     flag for headwater cell of inflow gauge
      70             :   !>       \param[in] "integer(i4), dimension(:) :: InflowGaugeNodeList"  gauge node list at L11
      71             :   !>       \param[in] "real(dp), dimension(:) :: InflowDischarge"         inflowing discharge at discharge gauge at
      72             :   !>       current day
      73             :   !>       \param[in] "integer(i4) :: nGauges"                            number of recording gauges
      74             :   !>       \param[in] "integer(i4), dimension(:) :: gaugeIndexList"       index list for outflow gauges
      75             :   !>       \param[in] "integer(i4), dimension(:) :: gaugeNodeList"        gauge node list at L11
      76             :   !>       \param[in] "logical :: map_flag"                               flag indicating whether routing resolution
      77             :   !>       iscoarser than hydrologic resolution
      78             :   !>       \param[in] "real(dp), dimension(:) :: L11_length"              L11 link length
      79             :   !>       \param[in] "real(dp), dimension(:) :: L11_slope"               L11 slope
      80             :   !>       \param[in] "real(dp), dimension(:) :: L11_FracFPimp"           L11 fraction of flood plain with impervios
      81             :   !>       cover
      82             : 
      83             :   !    INTENT(INOUT)
      84             :   !>       \param[inout] "real(dp), dimension(:) :: L11_C1"         L11 muskingum parameter 1
      85             :   !>       \param[inout] "real(dp), dimension(:) :: L11_C2"         L11 muskingum parameter 2
      86             :   !>       \param[inout] "real(dp), dimension(:) :: L11_qOut"       total runoff from L11 grid cells
      87             :   !>       \param[inout] "real(dp), dimension(:, :) :: L11_qTIN"    L11 inflow to the reach
      88             :   !>       \param[inout] "real(dp), dimension(:, :) :: L11_qTR"     L11 routed outflow
      89             :   !>       \param[inout] "real(dp), dimension(:) :: L11_qMod"       modelled discharge at each grid cell
      90             :   !>       \param[inout] "real(dp), dimension(:) :: GaugeDischarge" modelled discharge at each gauge
      91             : 
      92             :   !    HISTORY
      93             :   !>       \authors Stephan Thober
      94             : 
      95             :   !>       \date Aug 2015
      96             : 
      97             :   ! Modifications:
      98             :   ! Stephan Thober Sep 2015 - using arguments instead of global variables
      99             :   ! Stephan Thober Sep 2015 - added variables for routing resolution higher than hydrologic resolution
     100             :   ! Stephan Thober May 2016 - added check whether gauge is actually inside modelling domain before copying simulated runoff
     101             :   ! Stephan Thober Nov 2016 - implemented second routing process i.e. adaptive timestep
     102             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     103             : 
     104      634236 :   subroutine mRM_routing( &
     105      634236 :     read_states, processCase, global_routing_param, L1_total_runoff, L1_areaCell, L1_L11_Id, &
     106     1268472 :     L11_areaCell, L11_L1_Id, L11_netPerm, L11_fromN, L11_toN, L11_nOutlets, timestep, tsRoutFactor, &
     107      634236 :     nNodes, nInflowGauges, InflowGaugeIndexList, InflowGaugeHeadwater, InflowGaugeNodeList, &
     108     1902708 :     InflowDischarge, nGauges, gaugeIndexList, gaugeNodeList, map_flag, L11_length, L11_slope, &
     109     2536944 :     L11_FracFPimp, L11_C1, L11_C2, L11_qOut, L11_qTIN, L11_qTR, L11_qMod, GaugeDischarge &
     110             :   )
     111             : 
     112             :     use mo_constants, only : T0_dp
     113             :     use mo_mrm_global_variables, only : riv_temp_pcs, is_start
     114             :     use mo_mrm_mpr, only : reg_rout
     115             :     use mo_mrm_pre_routing, only : L11_runoff_acc, add_inflow
     116             :     ! use mo_mrm_riv_temp_class, only : riv_temp_type
     117             : 
     118             :     implicit none
     119             : 
     120             :     ! whether states are derived from restart file
     121             :     logical, intent(in) :: read_states
     122             :     ! Process switch for routing
     123             :     integer(i4), intent(in) :: processCase
     124             :     ! routing parameters
     125             :     real(dp), dimension(:), intent(in) :: global_routing_param
     126             :     ! total runoff from L1 grid cells
     127             :     real(dp), dimension(:), intent(in) :: L1_total_runoff
     128             :     ! L1 cell area
     129             :     real(dp), dimension(:), intent(in) :: L1_areaCell
     130             :     ! L1 cell ids on L11
     131             :     integer(i4), dimension(:), intent(in) :: L1_L11_Id
     132             :     ! L11 cell area
     133             :     real(dp), dimension(:), intent(in) :: L11_areaCell
     134             :     ! L11 cell ids on L1
     135             :     integer(i4), dimension(:), intent(in) :: L11_L1_Id
     136             :     ! L11 routing order
     137             :     integer(i4), dimension(:), intent(in) :: L11_netPerm
     138             :     ! L11 source grid cell order
     139             :     integer(i4), dimension(:), intent(in) :: L11_fromN
     140             :     ! L11 target grid cell order
     141             :     integer(i4), dimension(:), intent(in) :: L11_toN
     142             :     ! L11 number of outlets/sinks
     143             :     integer(i4), intent(in) :: L11_nOutlets
     144             :     ! simulation timestep in [h]
     145             :     integer(i4), intent(in) :: timestep
     146             :     ! factor between routing timestep and hydrological timestep
     147             :     real(dp), intent(in) :: tsRoutFactor
     148             :     ! number of nodes
     149             :     integer(i4), intent(in) :: nNodes
     150             :     ! number of inflow gauges
     151             :     integer(i4), intent(in) :: nInflowGauges
     152             :     ! index list of inflow gauges
     153             :     integer(i4), dimension(:), intent(in) :: InflowGaugeIndexList
     154             :     ! flag for headwater cell of inflow gauge
     155             :     logical, dimension(:), intent(in) :: InflowGaugeHeadwater
     156             :     ! gauge node list at L11
     157             :     integer(i4), dimension(:), intent(in) :: InflowGaugeNodeList
     158             :     ! inflowing discharge at discharge gauge at current day
     159             :     real(dp), dimension(:), intent(in) :: InflowDischarge
     160             :     ! number of recording gauges
     161             :     integer(i4), intent(in) :: nGauges
     162             :     ! index list for outflow gauges
     163             :     integer(i4), dimension(:), intent(in) :: gaugeIndexList
     164             :     ! gauge node list at L11
     165             :     integer(i4), dimension(:), intent(in) :: gaugeNodeList
     166             :     ! flag indicating whether routing resolution iscoarser than hydrologic resolution
     167             :     logical, intent(in) :: map_flag
     168             :     ! L11 link length
     169             :     real(dp), dimension(:), intent(in) :: L11_length
     170             :     ! L11 slope
     171             :     real(dp), dimension(:), intent(in) :: L11_slope
     172             :     ! L11 fraction of flood plain with impervios cover
     173             :     real(dp), dimension(:), intent(in) :: L11_FracFPimp
     174             :     ! L11 muskingum parameter 1
     175             :     real(dp), dimension(:), intent(inout) :: L11_C1
     176             :     ! L11 muskingum parameter 2
     177             :     real(dp), dimension(:), intent(inout) :: L11_C2
     178             :     ! total runoff from L11 grid cells
     179             :     real(dp), dimension(:), intent(inout) :: L11_qOut
     180             :     ! L11 inflow to the reach
     181             :     real(dp), dimension(:, :), intent(inout) :: L11_qTIN
     182             :     ! L11 routed outflow
     183             :     real(dp), dimension(:, :), intent(inout) :: L11_qTR
     184             :     ! modelled discharge at each grid cell
     185             :     real(dp), dimension(:), intent(inout) :: L11_qMod
     186             :     ! modelled discharge at each gauge
     187             :     real(dp), dimension(:), intent(inout) :: GaugeDischarge
     188             : 
     189             :     integer(i4) :: s11, e11 ! only for riv temp routing
     190             :     integer(i4) :: gg
     191             :     integer(i4) :: tt
     192             :     ! number of routing loops
     193             :     integer(i4) :: rout_loop
     194             :     ! variable for accumulation
     195    28516260 :     real(dp), dimension(size(L11_qMod, dim = 1)) :: L11_qAcc
     196      634236 :     real(dp), dimension(:), allocatable :: L11_E_Acc
     197             : 
     198      634236 :     if ( riv_temp_pcs%active ) then
     199             :       ! allocate accumulated temperature energy
     200       26208 :       allocate(L11_E_Acc(size(L11_qMod, dim = 1)))
     201      229320 :       L11_E_Acc = 0._dp ! init to zero
     202             :       ! get shortcuts for start end ending of current L11-domain
     203        6552 :       s11 = riv_temp_pcs%s11
     204        6552 :       e11 = riv_temp_pcs%e11
     205             :     end if
     206             : 
     207      634236 :     if (is_start) is_start = .false.
     208             : 
     209             :     ! this is using the sealed fraction for determining the routing parameters
     210             :     ! MPR has already been done
     211      634236 :     if (processCase .eq. 1_i4 .AND. (.not. read_states)) then
     212             :       ! for a single node model run
     213      495096 :       if (nNodes .GT. 1) then
     214           0 :         call reg_rout(global_routing_param, &
     215      495096 :                 L11_length, L11_slope, L11_FracFPimp(: nNodes - L11_nOutlets), &
     216      990192 :                 real(timeStep, dp), L11_C1(: nNodes - L11_nOutlets), L11_C2(: nNodes - L11_nOutlets))
     217             :       end if
     218             :     end if
     219             : 
     220             :     ! =====================================================================
     221             :     ! NOW, EXECUTE ROUTING
     222             :     ! ====================================================================
     223             :     ! calculate number of routing loops
     224      634236 :     rout_loop = max(1_i4, nint(1._dp / tsRoutFactor))
     225             : 
     226             :     ! runoff accumulation from L1 to L11 level
     227             :     call L11_runoff_acc( &
     228             :       L1_total_runoff, L1_areaCell, L1_L11_Id, &
     229             :       L11_areaCell, L11_L1_Id, timeStep, & ! Intent IN
     230             :       map_flag, & ! Intent IN
     231             :       L11_qOut & ! Intent OUT
     232      634236 :     )
     233             :     ! add inflow
     234             :     call add_inflow( &
     235             :       nInflowGauges, &
     236             :       InflowGaugeIndexList, &
     237             :       InflowGaugeHeadwater, &
     238             :       InflowGaugeNodeList, &
     239             :       InflowDischarge, & ! Intent IN
     240             :       L11_qOUT & ! Intent INOUT
     241      634236 :     )
     242             :     ! for a single node model run
     243      634236 :     if(nNodes .GT. 1) then
     244    28516260 :       L11_qAcc = 0._dp
     245             :       ! routing multiple times if timestep is smaller than 1
     246     1268472 :       do tt = 1, rout_loop
     247             :         ! routing of water within river reaches
     248             :         call L11_routing( &
     249             :           nNodes, &
     250             :           nNodes - L11_nOutlets, &
     251             :           L11_netPerm, &
     252             :           L11_fromN, & ! Intent IN
     253             :           L11_toN, & ! Intent IN
     254             :           L11_C1, & ! Intent IN
     255             :           L11_C2, & ! Intent IN
     256             :           L11_qOut, & ! Intent IN
     257             :           nInflowGauges, & ! Intent IN
     258             :           InflowGaugeHeadwater, & ! Intent IN
     259             :           InflowGaugeNodeList, & ! Intent IN
     260             :           L11_qTIN, & ! Intent INOUT
     261             :           L11_qTR, & ! Intent INOUT
     262           0 :           L11_Qmod & ! Intent OUT
     263      634236 :         )
     264             :         ! accumulate values of individual subtimesteps
     265    28516260 :         L11_qAcc = L11_qAcc + L11_qMod
     266             :         ! do the temperature routing
     267     1268472 :         if ( riv_temp_pcs%active ) then
     268             :           call riv_temp_pcs%L11_routing_E( &
     269             :             nNodes - L11_nOutlets, &
     270             :             L11_netPerm, &
     271             :             L11_fromN, & ! Intent IN
     272             :             L11_toN, & ! Intent IN
     273             :             L11_C1, & ! Intent IN
     274             :             L11_C2, & ! Intent IN
     275             :             nInflowGauges, & ! Intent IN
     276             :             InflowGaugeHeadwater, & ! Intent IN
     277             :             InflowGaugeNodeList, & ! Intent IN
     278             :             L11_qTR(:, 1), & ! Intent IN
     279             :             L11_Qmod & ! Intent IN
     280        6552 :           )
     281             :         end if
     282             :       end do
     283             :       ! calculate mean over routing period (timestep)
     284    28516260 :       L11_qMod = L11_qAcc / real(rout_loop, dp)
     285             :     else
     286           0 :       L11_Qmod = L11_qOUT
     287           0 :       if ( riv_temp_pcs%active ) riv_temp_pcs%river_temp(s11 : e11) = &
     288           0 :         max(riv_temp_pcs%delta_T, riv_temp_pcs%netNode_E_out(s11 : e11) / L11_qOUT - T0_dp)
     289             :     end if
     290             : 
     291             :     !----------------------------------------------------------------------
     292             :     ! FOR STORING the optional arguments
     293             :     !
     294             :     ! FOR RUNOFF
     295             :     ! NOTE:: Node ID for a given gauging station is stored at gaugeindex's
     296             :     !        index in runoff. In consequence the gauges in runoff are
     297             :     !        ordered corresponing to gauge%Q(:,:)
     298             :     !----------------------------------------------------------------------
     299     1294704 :     do gg = 1, nGauges
     300     1294704 :       GaugeDischarge(gaugeIndexList(gg)) = L11_Qmod(gaugeNodeList(gg))
     301             :     end do
     302             : 
     303      634236 :   end subroutine mRM_routing
     304             : 
     305             :   ! ------------------------------------------------------------------
     306             : 
     307             :   !    NAME
     308             :   !        L11_routing
     309             : 
     310             :   !    PURPOSE
     311             :   !>       \brief Performs runoff routing for mHM at L11 upscaled network
     312             :   !>       (\ref fig_routing "Routing Network").
     313             :   !>       \details
     314             :   !>       Hydrograph routing is carried out with the Muskingum algorithm
     315             :   !>       \cite CMM1988.  This simplification of the St. Venant
     316             :   !>       equations is justified in mHM because the potential areas of
     317             :   !>       application of this model would hardly exhibit abruptly
     318             :   !>       changing hydrographs with supercritical flows.  The discharge
     319             :   !>       leaving the river reach located on cell \f$ i \f$ \f$
     320             :   !>       Q_{i}^{1}(t) \f$ at time step \f$ t \f$ can be determined by
     321             :   !>       \f[ Q_{i}^{1}(t) =  Q_{i}^{1}(t-1)
     322             :   !>       + c_{1} \left( Q_{i}^{0}(t-1) - Q_{i}^{1}(t-1) \right)
     323             :   !>       + c_{2} \left( Q_{i}^{0}(t)   - Q_{i}^{0}(t-1) \right) \f]
     324             :   !>       with
     325             :   !>       \f[  Q_{i}^{0}(t) = Q_{i'}(t) + Q_{i'}^{1}(t) \f]
     326             :   !>       \f[ c_{1}= \frac{\Delta t} { \kappa (1- \xi ) + \frac{\Delta t}{2} } \f]
     327             :   !>       \f[ c_{2}= \frac{ \frac{\Delta t}{2} - \kappa \xi} { \kappa (1- \xi)
     328             :   !>       + \frac{\Delta t}{2} } \f]
     329             :   !>       where
     330             :   !>       \f$ Q_{i}^{0} \f$ and \f$ Q_{i}^{1} \f$ denote the discharge
     331             :   !>       entering and leaving the river reach located on cell \f$ i \f$
     332             :   !>       respectively.
     333             :   !>       \f$ Q_{i'} \f$ is the contribution from the upstream cell \f$
     334             :   !>       i'\f$.
     335             :   !>       \f$ \kappa \f$ Muskingum travel time parameter.
     336             :   !>       \f$ \xi \f$ Muskingum attenuation parameter.
     337             :   !>       \f$ \Delta t \f$ time interval in hours.
     338             :   !>       \f$ t \f$ Time index for each \f$ \Delta t \f$ interval.
     339             :   !>       To improve performance, a routing sequence "netPerm" is
     340             :   !>       required. This permutation is determined in the mo_init_mrm
     341             :   !>       routine.
     342             : 
     343             :   !>       \details TODO: add description
     344             : 
     345             :   !    INTENT(IN)
     346             :   !>       \param[in] "integer(i4) :: nNodes"                       number of network nodes = nCells1
     347             :   !>       \param[in] "integer(i4) :: nLinks"                       number of stream segment (reaches)
     348             :   !>       \param[in] "integer(i4), dimension(:) :: netPerm"        routing order of a given domain (permutation)
     349             :   !>       \param[in] "integer(i4), dimension(:) :: netLink_fromN"  from node
     350             :   !>       \param[in] "integer(i4), dimension(:) :: netLink_toN"    to node
     351             :   !>       \param[in] "real(dp), dimension(:) :: netLink_C1"        routing parameter  C1 (\cite CMM1988 p. 25-41)
     352             :   !>       \param[in] "real(dp), dimension(:) :: netLink_C2"        routing parameters C2 (id)
     353             :   !>       \param[in] "real(dp), dimension(:) :: netNode_qOUT"      Total outflow from cells (given domain) L11 at time
     354             :   !>       tt in [m3 s-1]
     355             :   !>       \param[in] "integer(i4) :: nInflowGauges"                [-]      number of inflow points
     356             :   !>       \param[in] "logical, dimension(:) :: InflowHeadwater"    [-]      if to consider headwater cells of inflow
     357             :   !>       gauge
     358             :   !>       \param[in] "integer(i4), dimension(:) :: InflowNodeList" [-]      L11 ID of inflow points
     359             : 
     360             :   !    INTENT(INOUT)
     361             :   !>       \param[inout] "real(dp), dimension(:, :) :: netNode_qTIN" [m3 s-1] Total inputs at t-1 and t
     362             :   !>       \param[inout] "real(dp), dimension(:, :) :: netNode_qTR"  [m3 s-1] Transformed outflow leaving
     363             :   !>       node I (Muskingum)
     364             : 
     365             :   !    INTENT(OUT)
     366             :   !>       \param[out] "real(dp), dimension(nNodes) :: netNode_Qmod" [m3 s-1] Simulated routed discharge
     367             : 
     368             :   !    HISTORY
     369             :   !>       \authors Luis Samaniego
     370             : 
     371             :   !>       \date Dec 2005
     372             : 
     373             :   ! Modifications:
     374             :   ! Luis Samaniego Feb 2008 - routing module (cells)
     375             :   ! Rohini Kumar   Aug 2011 - vector version of mHM-UFZ
     376             :   !                Nov 2011 - parallel version
     377             :   ! Luis Samaniego Jan 2013 - modularization, documentation
     378             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     379             : 
     380     1902708 :   subroutine L11_routing(nNodes, nLinks, netPerm, netLink_fromN, netLink_toN, netLink_C1, netLink_C2, netNode_qOUT, &
     381     1902708 :                         nInflowGauges, InflowHeadwater, InflowNodeList, netNode_qTIN, netNode_qTR, netNode_Qmod)
     382             :     implicit none
     383             : 
     384             :     ! number of network nodes = nCells1
     385             :     integer(i4), intent(in) :: nNodes
     386             :     ! number of stream segment (reaches)
     387             :     integer(i4), intent(in) :: nLinks
     388             :     ! routing order of a given domain (permutation)
     389             :     integer(i4), dimension(:), intent(in) :: netPerm
     390             :     ! from node
     391             :     integer(i4), dimension(:), intent(in) :: netLink_fromN
     392             :     ! to node
     393             :     integer(i4), dimension(:), intent(in) :: netLink_toN
     394             :     ! routing parameter  C1 (\cite CMM1988 p. 25-41)
     395             :     real(dp), dimension(:), intent(in) :: netLink_C1
     396             :     ! routing parameters C2 (id)
     397             :     real(dp), dimension(:), intent(in) :: netLink_C2
     398             :     ! Total outflow from cells (given domain) L11 at time tt in [m3 s-1]
     399             :     real(dp), dimension(:), intent(in) :: netNode_qOUT
     400             :     ! [-]      number of inflow points
     401             :     integer(i4), intent(in) :: nInflowGauges
     402             :     ! [-]      if to consider headwater cells of inflow gauge
     403             :     logical, dimension(:), intent(in) :: InflowHeadwater
     404             :     ! [-]      L11 ID of inflow points
     405             :     integer(i4), dimension(:), intent(in) :: InflowNodeList
     406             :     ! [m3 s-1] Total inputs at t-1 and t
     407             :     real(dp), dimension(:, :), intent(inout) :: netNode_qTIN
     408             :     ! [m3 s-1] Transformed outflow leaving node I (Muskingum)
     409             :     real(dp), dimension(:, :), intent(inout) :: netNode_qTR
     410             :     ! [m3 s-1] Simulated routed discharge
     411             :     real(dp), dimension(nNodes), intent(out) :: netNode_Qmod
     412             : 
     413             :     integer(i4) :: g, i, k, iNode, tNode
     414             :     ! current routing state (2)
     415             :     integer(i4), parameter :: IT = 2
     416             :     ! past routing state (1)
     417             :     integer(i4), parameter :: IT1 = 1
     418             : 
     419             :     ! Entry value for the auxiliary vectors
     420             :     !   netNode_qTIN(iNode,:)
     421             :     !   netNode_qTR(iNode,:)
     422             :     ! which store current and past states of
     423             :     ! incoming and outgoing of discharge at iNode
     424             :     !--------------------------------------------------------------------------
     425             :     !                             Muskingum Flood Routing
     426             :     !--------------------------------------------------------------------------
     427             :     ! initialize total input at point time IT in all nodes
     428    28516260 :     netNode_qTIN(:, IT) = 0.0_dp
     429             :     !--------------------------------------------------------------------------
     430             :     ! Links in sequential mode .... with single node
     431             :     !--------------------------------------------------------------------------
     432             :     ! ST - decent parallelization has to be done!!!
     433             :     !!$OMP parallel
     434             :     !!$OMP do private(g, i, inode, tnode)
     435    27882024 :     do k = 1, nLinks
     436             :       ! get LINK routing order -> i
     437    27247788 :       i = netPerm(k)
     438    27247788 :       iNode = netLink_fromN(i)
     439    27247788 :       tNode = netLink_toN(i)
     440             : 
     441             :       ! accumulate all inputs in iNode
     442    27247788 :       netNode_qTIN(iNode, IT) = netNode_qTIN(iNode, IT) + netNode_qOUT(iNode)
     443             : 
     444             :       ! routing iNode
     445    27247788 :       netNode_qTR(iNode, IT) = netNode_qTR(iNode, IT1)                               &
     446    27247788 :               + netLink_C1(i) * (netNode_qTIN(iNode, IT1) - netNode_qTR (iNode, IT1)) &
     447    27247788 :               + netLink_C2(i) * (netNode_qTIN(iNode, IT) - netNode_qTIN(iNode, IT1))
     448             : 
     449             :       ! check if the inflow from upstream cells should be deactivated
     450    27247788 :       if (nInflowGauges .GT. 0) then
     451     1729728 :         do g = 1, nInflowGauges
     452             :           ! check if downstream Node (tNode) is inflow gauge and headwaters should be ignored
     453     1729728 :           if ((tNode == InflowNodeList(g)) .AND. (.NOT. InflowHeadwater(g))) netNode_qTR(iNode, IT) = 0.0_dp
     454             :         end do
     455             :       end if
     456             : 
     457             :       ! add routed water to downstream node
     458    27882024 :       netNode_qTIN(tNode, IT) = netNode_qTIN(tNode, IT) + netNode_qTR(iNode, IT)
     459             :     end do
     460             :     !!$OMP end do
     461             :     !!$OMP end parallel
     462             : 
     463             :     !--------------------------------------------------------------------------
     464             :     ! Accumulate all inputs in tNode (netNode_qOUT) ONLY for last link
     465             :     !--------------------------------------------------------------------------
     466      634236 :     tNode = netLink_toN(netPerm(nLinks))
     467      634236 :     netNode_qTIN(tNode, IT) = netNode_qTIN(tNode, IT) + netNode_qOUT(tNode)
     468             : 
     469             :     !--------------------------------------------------------------------------
     470             :     ! save modeled discharge at time step tt then shift flow storages
     471             :     ! (NOTE aggregation to daily values to be done outside)
     472             :     !--------------------------------------------------------------------------
     473             :     ! !!$OMP parallel
     474             :     ! store generated discharge
     475    29150496 :     netNode_Qmod(1 : nNodes) = netNode_qTIN(1 : nNodes, IT)
     476             :     ! backflow t-> t-1
     477    29150496 :     netNode_qTR(1 : nNodes, IT1) = netNode_qTR(1 : nNodes, IT)
     478    29150496 :     netNode_qTIN(1 : nNodes, IT1) = netNode_qTIN(1 : nNodes, IT)
     479             :     ! !!$OMP end parallel
     480             : 
     481      634236 :   end subroutine L11_routing
     482             : 
     483             : END MODULE mo_mrm_routing

Generated by: LCOV version 1.16