LCOV - code coverage report
Current view: top level - mRM - mo_mrm_net_startup.f90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 784 883 88.8 %
Date: 2024-04-30 08:53:32 Functions: 14 15 93.3 %

          Line data    Source code
       1             : !> \file mo_mrm_net_startup.f90
       2             : !> \brief \copybrief mo_mrm_net_startup
       3             : !> \details \copydetails mo_mrm_net_startup
       4             : 
       5             : !> \brief Startup drainage network for mHM.
       6             : !> \details This module initializes the drainage network at L11 in mHM.
       7             : !! - Delineation of drainage network at level 11.
       8             : !! - Setting network topology (i.e. nodes and link).
       9             : !! - Determining routing order.
      10             : !! - Determining cell locations for network links.
      11             : !! - Find drainage outlet.
      12             : !! - Determine stream (links) features.
      13             : !!
      14             : !> \changelog
      15             : !! - Rohini Kumar May 2014
      16             : !!   - cell area calulation based on a regular lat-lon grid or on a regular X-Y coordinate system
      17             : !! - Robert Schweppe Jun 2018
      18             : !!   - refactoring and reformatting
      19             : !> \authors Luis Samaniego
      20             : !> \date Dec 2012
      21             : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
      22             : !! mHM is released under the LGPLv3+ license \license_note
      23             : !> \ingroup f_mrm
      24             : module mo_mrm_net_startup
      25             :   use mo_kind, only : i4, dp
      26             :   use mo_message, only : message, error_message
      27             : 
      28             :   implicit none
      29             : 
      30             :   PUBLIC :: L11_L1_mapping
      31             :   PUBLIC :: L11_flow_direction
      32             :   PUBLIC :: L11_set_network_topology
      33             :   PUBLIC :: L11_routing_order
      34             :   PUBLIC :: L11_link_location
      35             :   PUBLIC :: L11_set_drain_outlet_gauges
      36             :   PUBLIC :: L11_stream_features
      37             :   PUBLIC :: L11_fraction_sealed_floodplain
      38             :   PUBLIC :: get_distance_two_lat_lon_points
      39             :   PUBLIC :: L11_flow_accumulation
      40             :   PUBLIC :: L11_calc_celerity
      41             : contains
      42             : 
      43             :   !    NAME
      44             :   !        L11_L1_mapping
      45             : 
      46             :   !    PURPOSE
      47             :   !>       \brief TODO: add description
      48             : 
      49             :   !>       \details TODO: add description
      50             : 
      51             :   !    INTENT(IN)
      52             :   !>       \param[in] "integer(i4) :: iDomain" domain
      53             : 
      54             :   !    HISTORY
      55             :   !>       \authors Robert Schweppe
      56             : 
      57             :   !>       \date Jun 2018
      58             : 
      59             :   ! Modifications:
      60             : 
      61          24 :   subroutine L11_L1_mapping(iDomain)
      62             : 
      63             :     use mo_append, only : append
      64             :     use mo_common_constants, only : nodata_i4
      65             :     use mo_common_variables, only : level1
      66             :     use mo_mrm_global_variables, only : L11_L1_ID, L1_L11_ID, level11
      67             : 
      68             :     implicit none
      69             : 
      70             :     ! domain
      71             :     integer(i4), intent(in) :: iDomain
      72             : 
      73             :     integer(i4) :: nrows1, ncols1
      74             : 
      75             :     integer(i4) :: nrows11, ncols11
      76             : 
      77             :     integer(i4) :: kk
      78             : 
      79             :     integer(i4) :: icc, jcc
      80             : 
      81             :     integer(i4) :: iu, id, jl, jr
      82             : 
      83             :     ! mapping of L11 Id on L1
      84          24 :     integer(i4), dimension(:, :), allocatable :: L11Id_on_L1
      85             : 
      86             :     ! mapping of L1 Id on L11
      87          24 :     integer(i4), dimension(:, :), allocatable :: L1Id_on_L11
      88             : 
      89             :     ! dummy ID
      90          24 :     integer(i4), dimension(:, :), allocatable :: dummy_2d_id
      91             : 
      92          24 :     real(dp) :: cellFactorRbyH
      93             : 
      94             :     integer(i4) :: cellFactorRbyH_inv
      95             : 
      96             : 
      97          24 :     nrows1 = level1(iDomain)%nrows
      98          24 :     nrows11 = level11(iDomain)%nrows
      99          24 :     ncols1 = level1(iDomain)%ncols
     100          24 :     ncols11 = level11(iDomain)%ncols
     101             : 
     102             :     ! allocate variables for mapping L11 Ids and L1 Ids
     103          96 :     allocate (L11Id_on_L1  (nrows1, ncols1))
     104          96 :     allocate (L1Id_on_L11  (nrows11, ncols11))
     105        1983 :     L11Id_on_L1(:, :) = nodata_i4
     106        1812 :     L1Id_on_L11(:, :) = nodata_i4
     107             : 
     108             :     ! set cell factor for routing
     109          24 :     cellFactorRbyH = level11(iDomain)%cellsize / level1(iDomain)%cellsize
     110             : 
     111             :     ! set mapping
     112             :     ! create mapping between L11 and L1 for L11 resolution higher than L1 resolution
     113          24 :     if (cellFactorRbyH .lt. 1._dp) then
     114           0 :       allocate (dummy_2d_id  (nrows1, ncols1))
     115           0 :       dummy_2d_id = unpack(level1(iDomain)%Id, level1(iDomain)%mask, nodata_i4)
     116           0 :       cellFactorRbyH_inv = int(1. / cellFactorRbyH, i4)
     117           0 :       kk = 0
     118           0 :       do jcc = 1, ncols1
     119           0 :         do icc = 1, nrows1
     120           0 :           if(.not. level1(iDomain)%mask(icc, jcc)) cycle
     121           0 :           kk = kk + 1
     122             :           !
     123           0 :           iu = (icc - 1) * cellFactorRbyH_inv + 1
     124           0 :           id = min(icc * cellFactorRbyH_inv, nrows11)
     125           0 :           jl = (jcc - 1) * cellFactorRbyH_inv + 1
     126           0 :           jr = min(jcc * cellFactorRbyH_inv, ncols11)
     127             : 
     128           0 :           L1Id_on_L11(iu : id, jl : jr) = merge(dummy_2d_id(icc, jcc), nodata_i4, level11(iDomain)%mask(iu : id, jl : jr))
     129             :         end do
     130             :       end do
     131             :     else
     132          72 :       allocate (dummy_2d_id  (nrows11, ncols11))
     133          24 :       dummy_2d_id = unpack(level11(iDomain)%Id, level11(iDomain)%mask, nodata_i4)
     134             : 
     135          24 :       kk = 0
     136         250 :       do jcc = 1, ncols11
     137        1812 :         do icc = 1, nrows11
     138        1562 :           if(.not. level11(iDomain)%mask(icc, jcc)) cycle
     139         936 :           kk = kk + 1
     140             : 
     141             :           ! coord. of all corners L11 -> of finer scale level-1
     142         936 :           iu = (icc - 1) * nint(cellFactorRbyH, i4) + 1
     143         936 :           id = icc * nint(cellFactorRbyH, i4)
     144         936 :           jl = (jcc - 1) * nint(cellFactorRbyH, i4) + 1
     145         936 :           jr = jcc * nint(cellFactorRbyH, i4)
     146             : 
     147             :           ! constrain the range of up, down, left, and right boundaries
     148         936 :           if(iu < 1) iu = 1
     149         936 :           if(id > nrows1) id = nrows1
     150         936 :           if(jl < 1) jl = 1
     151         936 :           if(jr > ncols1) jr = ncols1
     152             : 
     153             :           ! Delimitation of level-11 cells on level-1 for L11 resolution lower than L1 resolution
     154        3170 :           L11Id_on_L1(iu : id, jl : jr) = merge(dummy_2d_id(icc, jcc), nodata_i4, level1(iDomain)%mask(iu : id, jl : jr))
     155             :         end do
     156             :       end do
     157             :     end if
     158             : 
     159             :     ! L1 data sets
     160          24 :     call append(L1_L11_Id, pack (L11Id_on_L1(:, :), level1(iDomain)%mask))
     161             :     ! L11 data sets
     162          24 :     call append(L11_L1_Id, pack (L1Id_on_L11(:, :), level11(iDomain)%mask))
     163             :     ! free space
     164          24 :     deallocate(L11Id_on_L1, L1Id_on_L11, dummy_2d_id)
     165             : 
     166          24 :   end subroutine L11_L1_mapping
     167             :   ! --------------------------------------------------------------------------
     168             : 
     169             :   !    NAME
     170             :   !        L11_flow_direction
     171             : 
     172             :   !    PURPOSE
     173             :   !>       \brief Determine the flow direction of the upscaled river
     174             :   !>       network at level L11.
     175             : 
     176             :   !>       \details The hydrographs generated at each cell are routed
     177             :   !>       through the drainage network at level-11 towards their
     178             :   !>       outlets. The drainage network at level-11 is conceptualized as a
     179             :   !>       graph whose nodes are hypothetically located at the center of
     180             :   !>       each grid cell connected by links that represent the river
     181             :   !>       reaches. The flow direction of a link correspond to the
     182             :   !>       direction towards a neighboring cell in which the net flow
     183             :   !>       accumulation (outflows minus inflows) attains its maximum
     184             :   !>       value. The net flow accumulation across a cell's boundary at
     185             :   !>       level-11 is estimated based on flow direction and flow
     186             :   !>       accumulation obtained at level-0 (\ref fig_routing "Routing
     187             :   !>       Network"). Note: level-1 denotes the modeling level, whereas
     188             :   !>       level-L11 is at least as coarse as level-1. Experience has
     189             :   !>       shown that routing can be done at a coarser resolution as
     190             :   !>       level-1, hence the level-11 was introduced.
     191             :   !>       \image html  routing.png "Upscaling routing network from L0 to L1 (or L11)"
     192             :   !>       \anchor fig_routing \image latex routing.pdf "Upscaling routing network from L0 to L1 (or L11)" width=14cm
     193             :   !>       The left panel depicts a schematic derivation of a drainage
     194             :   !>       network at the level-11 based on level-0 flow direction and
     195             :   !>       flow accumulation. The dotted line circle denotes the point
     196             :   !>       with the highest flow accumulation within a grid cell. The
     197             :   !>       topology of a tipical drainage routing network at level-11 is
     198             :   !>       shown in the right panel. Gray color areas denote the flood
     199             :   !>       plains estimated in mo_init_mrm, where the network
     200             :   !>       upscaling is also carried out.
     201             :   !>       For the sake of simplicity, it is assumed that all runoff leaving
     202             :   !>       a given cell would exit through a major direction.
     203             :   !>       Note that multiple outlets can exist within the modelling domain.
     204             :   !>       If a variable is added or removed here, then it also has to
     205             :   !>       be added or removed in the subroutine L11_config_set in
     206             :   !>       module mo_restart and in the subroutine set_L11_config in module
     207             :   !>       mo_set_netcdf_restart
     208             :   !>       ADDITIONAL INFORMATION
     209             :   !>       L11_flow_direction
     210             : 
     211             :   !    INTENT(IN)
     212             :   !>       \param[in] "integer(i4) :: iDomain" Domain Id
     213             : 
     214             :   !    HISTORY
     215             :   !>       \authors Luis Samaniego
     216             : 
     217             :   !>       \date Dec 2005
     218             : 
     219             :   ! Modifications:
     220             :   ! Luis Samaniego Jan 2013 - modular version
     221             :   ! Rohini Kumar   Apr 2014 - Case of L0 is same as L11 implemented
     222             :   ! Stephan Thober Aug 2015 - ported to mRM
     223             :   ! Stephan Thober Sep 2015 - create mapping between L11 and L1 if L11 resolution is higher than L1 resolution
     224             :   ! Stephan Thober May 2016 - introducing multiple outlets
     225             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     226             : 
     227          24 :   subroutine L11_flow_direction(iDomain)
     228             : 
     229          24 :     use mo_append, only : append
     230             :     use mo_common_constants, only : nodata_i4
     231             :     use mo_common_types, only: Grid
     232             :     use mo_common_variables, only : domainMeta, level0
     233             :     use mo_mrm_global_variables, only : L0_draSC, L0_fAcc, L0_fDir, L0_l11_remap, L11_colOut, L11_fDir, &
     234             :                                         L11_nOutlets, L11_rowOut, domain_mrm, level11
     235             :     use mo_string_utils, only : num2str
     236             : 
     237             :     implicit none
     238             : 
     239             :     ! Domain Id
     240             :     integer(i4), intent(in) :: iDomain
     241             : 
     242             :     integer(i4) :: nCells0
     243             : 
     244             :     integer(i4) :: nrows0, ncols0
     245             : 
     246             :     integer(i4) :: s0, e0
     247             : 
     248             :     integer(i4) :: nrows11, ncols11
     249             : 
     250             :     ! =  ncells11
     251             :     integer(i4) :: nNodes
     252             : 
     253             :     integer(i4) :: ii, jj, kk, ic, jc
     254             : 
     255             :     integer(i4) :: iu, id
     256             : 
     257             :     integer(i4) :: jl, jr
     258             : 
     259             :     integer(i4) :: iRow, jCol
     260             : 
     261          24 :     integer(i4), dimension(:, :), allocatable :: iD0
     262             : 
     263          24 :     integer(i4), dimension(:, :), allocatable :: fDir0
     264             : 
     265          24 :     integer(i4), dimension(:, :), allocatable :: fAcc0
     266             : 
     267          24 :     integer(i4), dimension(:, :), allocatable :: fDir11
     268             : 
     269             :     ! northing cell loc. of the Outlet
     270          24 :     integer(i4), dimension(:), allocatable :: rowOut
     271             : 
     272             :     ! easting cell loc. of the Outlet
     273          24 :     integer(i4), dimension(:), allocatable :: colOut
     274             : 
     275          24 :     integer(i4), dimension(:, :), allocatable :: draSC0
     276             : 
     277             :     ! output location in L0
     278          24 :     integer(i4), dimension(:, :), allocatable :: oLoc
     279             : 
     280             :     integer(i4) :: side
     281             : 
     282             :     integer(i4) :: fAccMax, idMax
     283             : 
     284             :     ! Number of outlet found
     285             :     integer(i4) :: Noutlet
     286             : 
     287             :     ! Number of outlets before this Domain
     288             :     integer(i4) :: old_Noutlet
     289             : 
     290             :     ! flag whether outlet is found
     291             :     logical :: is_outlet
     292             : 
     293             :     type(Grid), pointer :: level0_iDomain => null()
     294             : 
     295             : 
     296             :     !--------------------------------------------------------
     297             :     ! STEPS:
     298             :     ! 1) Estimate each variable locally for a given Domain
     299             :     ! 2) Pad each variable to its corresponding global one
     300             :     !--------------------------------------------------------
     301             :     ! initialize
     302          24 :     Noutlet = 0_i4
     303          24 :     level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
     304             :     !------------------------------------------------------------------
     305             :     !                Set Flow Direction at Level 11
     306             :     !                       Searching order
     307             :     !                             jl    jr
     308             :     !                          iu  +  4  +
     309             :     !                              3     1
     310             :     !                          id  +  2  +
     311             :     !------------------------------------------------------------------
     312             : 
     313             :     ! flow direction at level-11
     314             : 
     315             :     ! allocate
     316          24 :     nrows0 = level0_iDomain%nrows
     317          24 :     ncols0 = level0_iDomain%ncols
     318          24 :     nCells0 = level0_iDomain%ncells
     319          24 :     nrows11 = level11(iDomain)%nrows
     320          24 :     ncols11 = level11(iDomain)%ncols
     321          24 :     nNodes = level11(iDomain)%ncells
     322          24 :     s0 = level0_iDomain%iStart
     323          24 :     e0 = level0_iDomain%iEnd
     324             : 
     325          96 :     allocate (iD0         (nrows0, ncols0))
     326          72 :     allocate (fAcc0       (nrows0, ncols0))
     327          72 :     allocate (fDir0       (nrows0, ncols0))
     328          72 :     allocate (draSC0      (nrows0, ncols0))
     329          96 :     allocate (fDir11      (nrows11, ncols11))
     330          72 :     allocate (rowOut      (nNodes))
     331          48 :     allocate (colOut      (nNodes))
     332          24 :     allocate (oLoc        (1, 2))
     333             : 
     334             :     ! initialize
     335     2862360 :     iD0(:, :) = nodata_i4
     336     2862360 :     fAcc0(:, :) = nodata_i4
     337     2862360 :     fDir0(:, :) = nodata_i4
     338     2862360 :     draSC0(:, :) = nodata_i4
     339        1812 :     fDir11(:, :) = nodata_i4
     340         960 :     rowOut(:) = nodata_i4
     341         960 :     colOut(:) = nodata_i4
     342         120 :     oLoc(:, :) = nodata_i4
     343             : 
     344             : 
     345             :     ! get iD, fAcc, fDir at L0
     346          24 :     iD0(:, :) = UNPACK(level0_iDomain%Id, level0_iDomain%mask, nodata_i4)
     347          24 :     fAcc0(:, :) = UNPACK(L0_fAcc (s0 : e0), level0_iDomain%mask, nodata_i4)
     348          24 :     fDir0(:, :) = UNPACK(L0_fDir (s0 : e0), level0_iDomain%mask, nodata_i4)
     349             : 
     350             :     ! case where routing and input data scale is similar
     351          24 :     IF(nCells0 .EQ. nNodes) THEN
     352           0 :       oLoc(1, :) = maxloc(fAcc0, level0_iDomain%mask)
     353           0 :       kk = L0_L11_remap(iDomain)%lowres_id_on_highres(oLoc(1, 1), oLoc(1, 2))
     354             :       ! for a single node model run
     355           0 :       if(nCells0 .EQ. 1) then
     356           0 :         fDir11(1, 1) = fDir0(oLoc(1, 1), oLoc(1, 2))
     357             :       else
     358           0 :         fDir11(:, :) = fDir0(:, :)
     359             :       end if
     360           0 :       fDir11 (level11(iDomain)%CellCoor(kk, 1), level11(iDomain)%CellCoor(kk, 2)) = 0
     361             :       ! set location of main outlet in L11
     362           0 :       do kk = 1, nNodes
     363           0 :         ii = level11(iDomain)%CellCoor(kk, 1)
     364           0 :         jj = level11(iDomain)%CellCoor(kk, 2)
     365           0 :         rowOut(kk) = ii
     366           0 :         colOut(kk) = jj
     367             :       end do
     368           0 :       do kk = 1, ncells0
     369           0 :         ii = level0_iDomain%CellCoor(kk, 1)
     370           0 :         jj = level0_iDomain%CellCoor(kk, 2)
     371           0 :         draSC0(ii, jj) = kk
     372             :       end do
     373             : 
     374             :       ! case where routing and input data scale differs
     375             :     ELSE
     376             :       ! =======================================================================
     377             :       ! ST: find all cells whose downstream cells are outside the domain
     378             :       ! =======================================================================
     379     1066064 :       do ii = 1, nCells0
     380     1066040 :         iRow = level0_iDomain%CellCoor(ii, 1)
     381     1066040 :         jCol = level0_iDomain%CellCoor(ii, 2)
     382     1066040 :         call moveDownOneCell(fDir0(iRow, jCol), iRow, jCol)
     383             :         ! check whether new location is inside bound
     384     1066040 :         is_outlet = .False.
     385             :         if ((iRow .le. 0_i4) .or. (iRow .gt. nrows0) .or. &
     386     1066040 :             (jCol .le. 0_i4) .or. (jCol .gt. ncols0)) then
     387             :           is_outlet = .True.
     388             :         else
     389     1066040 :           if (fdir0(iRow, jCol) .le. 0) is_outlet = .True.
     390             :         end if
     391             :         !
     392     1066064 :         if (is_outlet) then
     393          24 :           Noutlet = Noutlet + 1_i4
     394             :           ! cell is an outlet
     395          24 :           if (Noutlet .eq. 1) then
     396          72 :             oLoc(1, :) = level0_iDomain%CellCoor(ii, :)
     397             :           else
     398           0 :             call append(oLoc, level0_iDomain%CellCoor(ii : ii, :))
     399             :           end if
     400             :           ! drain this cell into corresponding L11 cell
     401          24 :           kk = L0_L11_remap(iDomain)%lowres_id_on_highres(oLoc(Noutlet, 1), oLoc(Noutlet, 2))
     402          24 :           draSC0(oLoc(Noutlet, 1), oLoc(Noutlet, 2)) = kk
     403             :           ! check whether cell has maximum flow accumulation
     404             :           ! coord. of all corners
     405          24 :           iu = l0_l11_remap(iDomain)%upper_bound   (kk)
     406          24 :           id = l0_l11_remap(iDomain)%lower_bound (kk)
     407          24 :           jl = l0_l11_remap(iDomain)%left_bound (kk)
     408          24 :           jr = l0_l11_remap(iDomain)%right_bound(kk)
     409       52992 :           if (maxval(facc0(iu : id, jl : jr)) .eq. facc0(oLoc(Noutlet, 1), oLoc(Noutlet, 2)))  then
     410             :             ! set location of outlet at L11
     411          24 :             rowOut(kk) = oLoc(Noutlet, 1)
     412          24 :             colOut(kk) = oLoc(Noutlet, 2)
     413          24 :             fdir11(level11(iDomain)%CellCoor(kk, 1), level11(iDomain)%CellCoor(kk, 2)) = 0
     414             :           end if
     415             :         end if
     416             :       end do
     417             : 
     418             :       ! finding cell L11 outlets -  using L0_fAcc
     419             : 
     420         960 :       do kk = 1, nNodes
     421             : 
     422             :         ! exclude outlet L11
     423         936 :         if (rowOut(kk) > 0) cycle
     424             : 
     425         912 :         ic = level11(iDomain)%CellCoor(kk, 1)
     426         912 :         jc = level11(iDomain)%CellCoor(kk, 2)
     427             : 
     428             :         ! coord. of all corners
     429         912 :         iu = l0_l11_remap(iDomain)%upper_bound (kk)
     430         912 :         id = l0_l11_remap(iDomain)%lower_bound (kk)
     431         912 :         jl = l0_l11_remap(iDomain)%left_bound (kk)
     432         912 :         jr = l0_l11_remap(iDomain)%right_bound(kk)
     433             : 
     434         912 :         fAccMax = -9
     435         912 :         idMax = 0
     436         912 :         side = -1
     437             :         ! searching on side 4
     438       39504 :         do jj = jl, jr
     439       38592 :           if ((fAcc0(iu, jj) > fAccMax)  .and. &
     440             :               (fDir0(iu, jj) ==  32 .or.  &
     441             :               fDir0(iu, jj) ==  64 .or.  &
     442         912 :               fDir0(iu, jj) == 128)) then
     443        1858 :             fAccMax = fAcc0(iu, jj)
     444        1858 :             idMax = id0(iu, jj)
     445        1858 :             side = 4
     446             :           end if
     447             :         end do
     448             : 
     449             :         ! searching on side 1
     450       39504 :         do ii = iu, id
     451       38592 :           if ((fAcc0(ii, jr) > fAccMax)  .and. &
     452             :               (fDir0(ii, jr) ==   1 .or.  &
     453             :               fDir0(ii, jr) ==   2 .or.  &
     454         912 :               fDir0(ii, jr) == 128)) then
     455         418 :             fAccMax = fAcc0(ii, jr)
     456         418 :             idMax = id0(ii, jr)
     457         418 :             side = 1
     458             :           end if
     459             :         end do
     460             : 
     461             :         ! searching on side 2
     462       39504 :         do jj = jl, jr
     463       38592 :           if ((fAcc0(id, jj) > fAccMax)  .and. &
     464             :               (fDir0(id, jj) ==   2 .or.  &
     465             :               fDir0(id, jj) ==   4 .or.  &
     466         912 :               fDir0(id, jj) ==   8)) then
     467         546 :             fAccMax = fAcc0(id, jj)
     468         546 :             idMax = id0(id, jj)
     469         546 :             side = 2
     470             :           end if
     471             :         end do
     472             : 
     473             :         ! searching on side 3
     474       39504 :         do ii = iu, id
     475       38592 :           if ((fAcc0(ii, jl) > fAccMax)  .and. &
     476             :               (fDir0(ii, jl) ==   8 .or.  &
     477             :               fDir0(ii, jl) ==  16 .or.  &
     478         912 :               fDir0(ii, jl) ==  32)) then
     479         628 :             fAccMax = fAcc0(ii, jl)
     480         628 :             idMax = id0(ii, jl)
     481         628 :             side = 3
     482             :           end if
     483             :         end do
     484             : 
     485             :         ! set location of the cell-outlet (row, col) in L0
     486         912 :         ii = level0_iDomain%CellCoor(idMax, 1)
     487         912 :         jj = level0_iDomain%CellCoor(idMax, 2)
     488         912 :         rowOut(kk) = ii
     489         912 :         colOut(kk) = jj
     490         912 :         draSC0(ii, jj) = kk
     491             : 
     492             :         ! set fDir at L11
     493         936 :         if     (ii == iu .and.  jj == jl) then
     494           4 :           select case (fDir0(ii, jj))
     495             :           case (8, 16)
     496           0 :             fDir11(ic, jc) = 16
     497             :           case (32)
     498           0 :             fDir11(ic, jc) = 32
     499             :           case (64, 128)
     500           4 :             fDir11(ic, jc) = 64
     501             :           end select
     502         908 :         elseif (ii == iu .and.  jj == jr) then
     503           4 :           select case (fDir0(ii, jj))
     504             :           case (32, 64)
     505           2 :             fDir11(ic, jc) = 64
     506             :           case (128)
     507           0 :             fDir11(ic, jc) = 128
     508             :           case (1, 2)
     509           2 :             fDir11(ic, jc) = 1
     510             :           end select
     511         906 :         elseif (ii == id .and.  jj == jl) then
     512          62 :           select case (fDir0(ii, jj))
     513             :           case (2, 4)
     514           6 :             fDir11(ic, jc) = 4
     515             :           case (8)
     516           2 :             fDir11(ic, jc) = 8
     517             :           case (16, 32)
     518          56 :             fDir11(ic, jc) = 16
     519             :           end select
     520         850 :         elseif (ii == id .and.  jj == jr) then
     521           2 :           select case (fDir0(ii, jj))
     522             :           case (128, 1)
     523           0 :             fDir11(ic, jc) = 1
     524             :           case (2)
     525           0 :             fDir11(ic, jc) = 2
     526             :           case (4, 8)
     527           2 :             fDir11(ic, jc) = 4
     528             :           end select
     529             :         else
     530             :           ! cell on one side
     531          50 :           select case (side)
     532             :           case (1)
     533          50 :             fDir11(ic, jc) = 1
     534             :           case (2)
     535         130 :             fDir11(ic, jc) = 4
     536             :           case (3)
     537         354 :             fDir11(ic, jc) = 16
     538             :           case (4)
     539         314 :             fDir11(ic, jc) = 64
     540             :           case default
     541         872 :              call error_message('Error L11_flow_direction: side = -1')
     542             :           end select
     543             :         end if
     544             : 
     545             :       end do
     546             : 
     547             :     END IF
     548             :     !--------------------------------------------------------
     549             :     ! Start padding up local variables to global variables
     550             :     !--------------------------------------------------------
     551             : 
     552             :     ! allocate space for row and col Outlet
     553          24 :     allocate(domain_mrm(iDomain)%L0_rowOutlet(1))
     554          24 :     allocate(domain_mrm(iDomain)%L0_colOutlet(1))
     555          24 :     domain_mrm(iDomain)%L0_Noutlet = nodata_i4
     556          72 :     domain_mrm(iDomain)%L0_rowOutlet = nodata_i4
     557          72 :     domain_mrm(iDomain)%L0_colOutlet = nodata_i4
     558             : 
     559             :     ! L0 data sets
     560             :     ! check whether L0 data is shared
     561             :     !ToDo: check if change is correct
     562          24 :     if (iDomain .eq. 1) then
     563          12 :       call append(L0_draSC, PACK (draSC0(:, :), level0_iDomain%mask))
     564          12 :     else if (domainMeta%L0DataFrom(iDomain) == iDomain) then
     565           9 :       call append(L0_draSC, PACK (draSC0(:, :), level0_iDomain%mask))
     566             :     end if
     567             : 
     568          24 :     domain_mrm(iDomain)%L0_Noutlet = Noutlet
     569             :     ! set L0 outlet coordinates
     570          24 :     old_Noutlet = size(domain_mrm(iDomain)%L0_rowOutlet, dim = 1)
     571          24 :     if (Noutlet .le. old_Noutlet) then
     572          48 :       domain_mrm(iDomain)%L0_rowOutlet(: Noutlet) = oLoc(:, 1)
     573          48 :       domain_mrm(iDomain)%L0_colOutlet(: Noutlet) = oLoc(:, 2)
     574             :     else
     575             :       ! store up to size of old_Noutlet
     576           0 :       domain_mrm(iDomain)%L0_rowOutlet(: old_Noutlet) = oLoc(: old_Noutlet, 1)
     577           0 :       domain_mrm(iDomain)%L0_colOutlet(: old_Noutlet) = oLoc(: old_Noutlet, 2)
     578             :       ! enlarge rowOutlet and colOutlet in domain_mrm structure
     579             :       !TODO: do other domains also need to be enlarged accordingly???
     580           0 :       call append(domain_mrm(iDomain)%L0_rowOutlet, oLoc(old_Noutlet + 1 :, 1))
     581           0 :       call append(domain_mrm(iDomain)%L0_colOutlet, oLoc(old_Noutlet + 1 :, 2))
     582             :     end if
     583             : 
     584             :     ! L11 data sets
     585        1812 :     call append(L11_nOutlets, count(fdir11 .eq. 0_i4))
     586          24 :     call append(L11_fDir, PACK (fDir11(:, :), level11(iDomain)%mask))
     587          24 :     call append(L11_rowOut, rowOut(:))
     588          24 :     call append(L11_colOut, colOut(:))
     589             : 
     590             :     ! communicate
     591          24 :     call message('')
     592          24 :     call message('    Domain: ' // num2str(iDomain, '(i3)'))
     593          24 :     call message('      Number of outlets found at Level 0:.. ' // num2str(Noutlet, '(i7)'))
     594        1812 :     call message('      Number of outlets found at Level 11:. ' // num2str(count(fdir11 .eq. 0_i4), '(i7)'))
     595             : 
     596             :     ! free space
     597          24 :     deallocate(fDir0, fAcc0, fDir11, rowOut, colOut, draSC0)
     598             : 
     599          24 :   end subroutine L11_flow_direction
     600             : 
     601             :   ! ------------------------------------------------------------------
     602             : 
     603             :   !    NAME
     604             :   !        L11_set_network_topology
     605             : 
     606             :   !    PURPOSE
     607             :   !>       \brief Set network topology
     608             : 
     609             :   !>       \details Set network topology from and to node for all links
     610             :   !>       at level-11 (\ref fig_routing "Routing Network").
     611             :   !>       If a variable is added or removed here, then it also has to
     612             :   !>       be added or removed in the subroutine L11_config_set in
     613             :   !>       module mo_restart and in the subroutine set_L11_config in module
     614             :   !>       mo_set_netcdf_restart.
     615             : 
     616             :   !    INTENT(IN)
     617             :   !>       \param[in] "integer(i4) :: iDomain" Domain Id
     618             : 
     619             :   !    HISTORY
     620             :   !>       \authors Luis Samaniego
     621             : 
     622             :   !>       \date Dec 2005
     623             : 
     624             :   ! Modifications:
     625             :   ! Luis Samaniego Jan 2013 - modular version
     626             :   ! Stephan Thober Aug 2015 - ported to mRM
     627             :   ! Stephan Thober May 2016 - moved calculation of sink here
     628             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     629             : 
     630          24 :   subroutine L11_set_network_topology(iDomain)
     631             : 
     632          24 :     use mo_append, only : append
     633             :     use mo_common_constants, only : nodata_i4
     634             :     use mo_mrm_global_variables, only : L11_fDir, L11_fromN, L11_toN, level11
     635             : 
     636             :     implicit none
     637             : 
     638             :     ! Domain Id
     639             :     integer(i4), intent(in) :: iDomain
     640             : 
     641          24 :     integer(i4), dimension(:, :), allocatable :: fDir11
     642             : 
     643          24 :     integer(i4), dimension(:, :), allocatable :: dummy_2d_id
     644             : 
     645             :     integer(i4) :: jj, kk, ic, jc
     646             : 
     647             :     integer(i4) :: fn, tn
     648             : 
     649          24 :     integer(i4), dimension(:), allocatable :: nLinkFromN, nLinkToN
     650             : 
     651             : 
     652             :     !     Routing network vectors have nNodes size instead of nLinks to
     653             :     !     avoid the need of having two extra indices to identify a Domain.
     654             :     ! allocate
     655          72 :     allocate (nLinkFromN (level11(iDomain)%nCells))  ! valid from (1 : nLinks)
     656          48 :     allocate (nLinkToN   (level11(iDomain)%nCells))  ! "
     657          96 :     allocate (fDir11     (level11(iDomain)%nrows, level11(iDomain)%ncols))
     658          72 :     allocate (dummy_2d_id(level11(iDomain)%nrows, level11(iDomain)%ncols))
     659          24 :     dummy_2d_id = unpack(level11(iDomain)%Id, level11(iDomain)%mask, nodata_i4)
     660             : 
     661             : 
     662             :     ! initialize
     663         960 :     nLinkFromN(:) = nodata_i4
     664         960 :     nLinkToN(:) = nodata_i4
     665        1812 :     fDir11(:, :) = nodata_i4
     666             : 
     667             :     ! get grids of L11
     668          24 :     fDir11(:, :) = UNPACK(L11_fDir (level11(iDomain)%iStart : level11(iDomain)%iEnd), level11(iDomain)%mask, nodata_i4)
     669             : 
     670             :     ! ------------------------------------------------------------------
     671             :     !  network topology
     672             :     ! ------------------------------------------------------------------
     673             : 
     674          24 :     jj = 0
     675         960 :     do kk = 1, level11(iDomain)%nCells
     676         936 :       ic = level11(iDomain)%CellCoor(kk, 1)
     677         936 :       jc = level11(iDomain)%CellCoor(kk, 2)
     678         936 :       fn = kk
     679         936 :       call moveDownOneCell(fDir11(ic, jc), ic, jc)
     680         936 :       tn = dummy_2d_id(ic, jc)
     681         936 :       if (fn == tn) cycle
     682         912 :       jj = jj + 1
     683         912 :       nLinkFromN(jj) = fn
     684        1872 :       nLinkToN(jj) = tn
     685             :     end do
     686             : 
     687             :     !--------------------------------------------------------
     688             :     ! Start padding up local variables to global variables
     689             :     !--------------------------------------------------------
     690             : 
     691             :     ! L11 data sets
     692          24 :     call append(L11_fromN, nLinkFromN(:)) ! sinks are at the end
     693          24 :     call append(L11_toN, nLinkToN(:))
     694             : 
     695             :     ! free space
     696          24 :     deallocate (fDir11, nLinkFromN, nLinkToN)
     697             : 
     698          24 :   end subroutine L11_set_network_topology
     699             : 
     700             :   ! ------------------------------------------------------------------
     701             : 
     702             :   !    NAME
     703             :   !        L11_routing_order
     704             : 
     705             :   !    PURPOSE
     706             :   !>       \brief Find routing order, headwater cells and sink
     707             : 
     708             :   !>       \details Find routing order, headwater cells and sink.
     709             :   !>       If a variable is added or removed here, then it also has to
     710             :   !>       be added or removed in the subroutine L11_config_set in
     711             :   !>       module mo_restart and in the subroutine set_L11_config in module
     712             :   !>       mo_set_netcdf_restart
     713             : 
     714             :   !    INTENT(IN)
     715             :   !>       \param[in] "integer(i4) :: iDomain" Domain Id
     716             : 
     717             :   !    HISTORY
     718             :   !>       \authors Luis Samaniego
     719             : 
     720             :   !>       \date Dec 2005
     721             : 
     722             :   ! Modifications:
     723             :   ! Luis Samaniego Jan 2013 - modular version
     724             :   ! Sa. Ku.        Jan 2015 - corrected initialization of nLinkSink
     725             :   ! Stephan Thober Aug 2015 - ported to mRM
     726             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     727             : 
     728          24 :   subroutine L11_routing_order(iDomain)
     729             : 
     730          24 :     use mo_append, only : append
     731             :     use mo_common_constants, only : nodata_i4
     732             :     use mo_mrm_global_variables, only : L11_fDir, L11_fromN, L11_label, L11_nOutlets, L11_netPerm, L11_rOrder, L11_sink, L11_toN, &
     733             :                                         level11
     734             : 
     735             :     implicit none
     736             : 
     737             :     ! Domain Id
     738             :     integer(i4), intent(in) :: iDomain
     739             : 
     740             :     integer(i4) :: nLinks
     741             : 
     742             :     ! from node
     743          24 :     integer(i4), dimension(:), allocatable :: nLinkFromN
     744             : 
     745             :     ! to node
     746          24 :     integer(i4), dimension(:), allocatable :: nLinkToN
     747             : 
     748             :     ! network routing order
     749          24 :     integer(i4), dimension(:), allocatable :: nLinkROrder
     750             : 
     751             :     ! label Id [0='', 1=HeadWater, 2=Sink]
     752          24 :     integer(i4), dimension(:), allocatable :: nLinkLabel
     753             : 
     754             :     ! == .true. if sink node reached
     755          24 :     logical, dimension(:), allocatable :: nLinkSink
     756             : 
     757             :     ! routing order (permutation)
     758          24 :     integer(i4), dimension(:), allocatable :: netPerm
     759             : 
     760             :     integer(i4) :: ii, jj, kk
     761             : 
     762             :     logical :: flag
     763             : 
     764             : 
     765          24 :     nLinks = level11(iDomain)%nCells - L11_nOutlets(iDomain)
     766             :     !  Routing network vectors have nNodes size instead of nLinks to
     767             :     !  avoid the need of having two extra indices to identify a Domain.
     768             : 
     769             :     ! allocate
     770          72 :     allocate (nLinkFromN  (level11(iDomain)%nCells))  ! all vectors valid from (1 : nLinks)
     771          48 :     allocate (nLinkToN    (level11(iDomain)%nCells))
     772          48 :     allocate (nLinkROrder (level11(iDomain)%nCells))
     773          48 :     allocate (nLinkLabel  (level11(iDomain)%nCells))
     774          48 :     allocate (nLinkSink   (level11(iDomain)%nCells))
     775          48 :     allocate (netPerm     (level11(iDomain)%nCells))
     776             :     ! initialize
     777         960 :     nLinkFromN(:) = nodata_i4
     778         960 :     nLinkToN(:) = nodata_i4
     779         960 :     nLinkROrder(1 : nLinks) = 1
     780          24 :     nLinkROrder(level11(iDomain)%nCells) = nodata_i4
     781         936 :     nLinkLabel(1 : nLinks) = 0
     782          24 :     nLinkLabel(level11(iDomain)%nCells) = nodata_i4
     783         960 :     nLinkSink(:) = .FALSE.
     784         960 :     netPerm(:) = nodata_i4
     785             : 
     786             :     ! for a single node model run
     787          24 :     if(level11(iDomain)%nCells .GT. 1) then
     788             :       ! get network vectors of L11
     789         960 :       nLinkFromN(:) = L11_fromN (level11(iDomain)%iStart : level11(iDomain)%iEnd)
     790         960 :       nLinkToN(:) = L11_toN   (level11(iDomain)%iStart : level11(iDomain)%iEnd)
     791             : 
     792         936 :       loop1 : do ii = 1, nLinks
     793       31812 :         loop2 : do jj = 1, nLinks
     794       31474 :           if (jj == ii) cycle loop2
     795       30714 :           if (nLinkFromN(ii) == nLinkToN(jj)) then
     796         598 :             nLinkROrder(ii) = -9
     797             :           end if
     798       31028 :           if (nLinkROrder(ii) == -9) cycle loop1
     799             :         end do loop2
     800             :       end do loop1
     801             :       ! counting headwaters
     802             :       kk = 0
     803         936 :       do ii = 1, nLinks
     804         936 :         if (nLinkROrder(ii) == 1) then
     805         314 :           kk = kk + 1
     806         314 :           nLinkROrder(ii) = kk
     807         314 :           nLinkLabel(ii) = 1  ! 'Head Water'
     808             :         end if
     809             :       end do
     810             :       ! counting downstream
     811       11748 :       do while (minval(nLinkROrder(1 : nLinks)) < 0)
     812             :         !!  print *, count(nLinkROrder .lt. 0), minval(nLinkROrder)
     813       10836 :         loop3 : do ii = 1, nLinks
     814       10584 :           if (.NOT. nLinkROrder(ii) == -9) cycle loop3
     815       76524 :           flag = .TRUE.
     816       76524 :           loop4 : do jj = 1, nLinks
     817       76524 :             if (jj == ii .OR. nLinkFromN(ii)  /=  nLinkToN(jj)) then
     818             :               cycle loop4
     819        3716 :             else if (.NOT. (nLinkFromN(ii)  == nLinkToN(jj)  .AND. nLinkROrder(jj) > 0)) then
     820             :               flag = .FALSE.
     821             :               exit loop4
     822             :             else
     823             :             end if
     824             :           end do loop4
     825             : 
     826        2508 :           if (flag) then
     827         598 :             kk = kk + 1
     828         598 :             nLinkROrder(ii) = kk
     829             :           end if
     830             :         end do loop3
     831             :       end do
     832             : 
     833             :       ! identify sink cells
     834         936 :       do ii = 1, nLinks
     835         936 :         if (L11_fdir(level11(iDomain)%iStart + nLinkToN(ii) - 1_i4) .eq. 0_i4) nlinksink(ii) = .True.
     836             :       end do
     837         960 :       where(nlinksink) nLinkLabel = 2 !  'Sink'
     838             : 
     839             :       ! keep routing order
     840         936 :       do ii = 1, nLinks
     841         936 :         netPerm(nLinkROrder(ii)) = ii
     842             :       end do
     843             : 
     844             :       ! end of multi-node network design loop
     845             :     end if
     846             : 
     847             :     !--------------------------------------------------------
     848             :     ! Start padding up local variables to global variables
     849             :     !--------------------------------------------------------
     850             :     ! L11 network data sets
     851          24 :     call append(L11_rOrder, nLinkROrder(:))
     852          24 :     call append(L11_label, nLinkLabel(:))
     853          24 :     call append(L11_sink, nLinkSink(:))
     854          24 :     call append(L11_netPerm, netPerm(:))
     855             : 
     856             :     ! free space
     857          24 :     deallocate (nLinkFromN, nLinkToN, nLinkROrder, nLinkLabel, nLinkSink, netPerm)
     858             : 
     859          24 :   end subroutine L11_routing_order
     860             : 
     861             :   ! ------------------------------------------------------------------
     862             : 
     863             :   !    NAME
     864             :   !        L11_link_location
     865             : 
     866             :   !    PURPOSE
     867             :   !>       \brief Estimate the LO (row,col) location for each routing link at level L11
     868             : 
     869             :   !>       \details If a variable is added or removed here, then it also has to
     870             :   !>       be added or removed in the subroutine L11_config_set in
     871             :   !>       module mo_restart and in the subroutine set_L11_config in module
     872             :   !>       mo_set_netcdf_restart
     873             : 
     874             :   !    INTENT(IN)
     875             :   !>       \param[in] "integer(i4) :: iDomain" Domain Id
     876             : 
     877             :   !    HISTORY
     878             :   !>       \authors Luis Samaniego
     879             : 
     880             :   !>       \date Dec 2005
     881             : 
     882             :   ! Modifications:
     883             :   ! Luis Samaniego Jan 2013 - modular version
     884             :   ! Stephan Thober Aug 2015 - ported to mRM
     885             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     886             : 
     887          24 :   subroutine L11_link_location(iDomain)
     888             : 
     889          24 :     use mo_append, only : append
     890             :     use mo_common_constants, only : nodata_i4
     891             :     use mo_common_types, only: Grid
     892             :     use mo_common_variables, only : domainMeta, level0
     893             :     use mo_mrm_global_variables, only : L0_draSC, L0_fDir, L11_colOut, L11_fCol, L11_fRow, L11_fromN, &
     894             :                                         L11_nOutlets, L11_netPerm, L11_rowOut, L11_tCol, L11_tRow, domain_mrm, level11
     895             :     use mo_string_utils, only : num2str
     896             : 
     897             :     implicit none
     898             : 
     899             :     ! Domain Id
     900             :     integer(i4), intent(in) :: iDomain
     901             : 
     902             :     integer(i4) :: nLinks
     903             : 
     904             :     ! northing cell loc. of the Outlet
     905          24 :     integer(i4), dimension(:), allocatable :: rowOut
     906             : 
     907             :     ! easting cell loc. of the Outlet
     908          24 :     integer(i4), dimension(:), allocatable :: colOut
     909             : 
     910          24 :     integer(i4), dimension(:), allocatable :: nLinkFromN
     911             : 
     912          24 :     integer(i4), dimension(:), allocatable :: netPerm
     913             : 
     914          24 :     integer(i4), dimension(:), allocatable :: nLinkFromRow
     915             : 
     916          24 :     integer(i4), dimension(:), allocatable :: nLinkFromCol
     917             : 
     918          24 :     integer(i4), dimension(:), allocatable :: nLinkToRow
     919             : 
     920          24 :     integer(i4), dimension(:), allocatable :: nLinkToCol
     921             : 
     922          24 :     integer(i4), dimension(:, :), allocatable :: fDir0
     923             : 
     924          24 :     integer(i4), dimension(:, :), allocatable :: draSC0
     925             : 
     926             :     integer(i4) :: ii, rr, kk, s0, e0
     927             : 
     928             :     integer(i4) :: iNode, iRow, jCol, prevRow, prevCol
     929             : 
     930             :     ! output location in L0
     931          24 :     integer(i4), dimension(:, :), allocatable :: oLoc
     932             : 
     933             :     ! number of outlets in Domain
     934             :     integer(i4) :: nOutlets
     935             : 
     936             :     ! flag for finding outlet
     937             :     logical :: is_outlet
     938             : 
     939             :     type(Grid), pointer :: level0_iDomain => null()
     940             : 
     941             : 
     942          24 :     level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
     943          24 :     s0 = level0_iDomain%iStart
     944          24 :     e0 = level0_iDomain%iEnd
     945             : 
     946          24 :     nOutlets = L11_nOutlets(iDomain)
     947             : 
     948          24 :     nLinks = level11(iDomain)%nCells - nOutlets
     949             : 
     950             :     !  Routing network vectors have level11(iDomain)%nCells size instead of nLinks to
     951             :     !  avoid the need of having two extra indices to identify a Domain.
     952             :     ! allocate
     953          72 :     allocate (rowOut        (level11(iDomain)%nCells))
     954          48 :     allocate (colOut        (level11(iDomain)%nCells))
     955          48 :     allocate (nLinkFromN    (level11(iDomain)%nCells))  ! all network vectors valid from (1 : nLinks)
     956          48 :     allocate (netPerm       (level11(iDomain)%nCells))
     957          48 :     allocate (nLinkFromRow  (level11(iDomain)%nCells))
     958          48 :     allocate (nLinkFromCol  (level11(iDomain)%nCells))
     959          48 :     allocate (nLinkToRow    (level11(iDomain)%nCells))
     960          48 :     allocate (nLinkToCol    (level11(iDomain)%nCells))
     961          96 :     allocate (fDir0         (level0_iDomain%nrows, level0_iDomain%ncols))
     962          72 :     allocate (draSC0        (level0_iDomain%nrows, level0_iDomain%ncols))
     963             : 
     964             :     ! initialize
     965         984 :     rowOut = nodata_i4
     966         984 :     colOut = nodata_i4
     967         984 :     nLinkFromN = nodata_i4
     968         984 :     netPerm = nodata_i4
     969         984 :     nLinkFromRow = nodata_i4
     970         984 :     nLinkFromCol = nodata_i4
     971         984 :     nLinkToRow = nodata_i4
     972         984 :     nLinkToCol = nodata_i4
     973     2862384 :     fDir0 = nodata_i4
     974     2862384 :     draSC0 = nodata_i4
     975             : 
     976             :     ! for a single node model run
     977          24 :     if(level11(iDomain)%nCells .GT. 1) then
     978             :       ! get fDir at L0
     979          24 :       fDir0(:, :) = UNPACK(L0_fDir  (s0 : e0), level0_iDomain%mask, nodata_i4)
     980          24 :       draSC0(:, :) = UNPACK(L0_draSC (s0 : e0), level0_iDomain%mask, nodata_i4)
     981             : 
     982             :       ! get network vectors of L11
     983         960 :       nLinkFromN(:) = L11_fromN   (level11(iDomain)%iStart : level11(iDomain)%iEnd)
     984         960 :       netPerm(:) = L11_netPerm (level11(iDomain)%iStart : level11(iDomain)%iEnd)
     985         960 :       rowOut(:) = L11_rowOut  (level11(iDomain)%iStart : level11(iDomain)%iEnd)
     986         960 :       colOut(:) = L11_colOut  (level11(iDomain)%iStart : level11(iDomain)%iEnd)
     987             : 
     988             :       ! finding main outlet (row, col) in L0
     989          72 :       allocate(oLoc(Noutlets, 2))
     990          48 :       oLoc(:, 1) = domain_mrm(iDomain)%L0_rowOutlet(: Noutlets)
     991          48 :       oLoc(:, 2) = domain_mrm(iDomain)%L0_colOutlet(: Noutlets)
     992             : 
     993             :       ! Location of the stream-joint cells  (row, col)
     994         936 :       do rr = 1, nLinks
     995             : 
     996         912 :         ii = netPerm(rr)
     997         912 :         iNode = nLinkFromN(ii)
     998         912 :         iRow = rowOut(iNode)
     999         912 :         jCol = colOut(iNode)
    1000         912 :         call moveDownOneCell(fDir0(iRow, jcol), iRow, jcol)
    1001             :         ! set "from" cell
    1002         912 :         nLinkFromRow(ii) = iRow
    1003         912 :         nLinkFromCol(ii) = jCol
    1004             : 
    1005             :         ! check whether this location is an outlet
    1006         912 :         is_outlet = .False.
    1007        1824 :         do kk = 1, Noutlets
    1008        1824 :           if (iRow .eq. oLoc(kk, 1) .and. jCol .eq. oLoc(kk, 2)) is_outlet = .True.
    1009             :         end do
    1010             : 
    1011        1848 :         if (is_outlet) then
    1012             : 
    1013           0 :           nLinkToRow(ii) = iRow
    1014           0 :           nLinkToCol(ii) = jCol
    1015             : 
    1016             :         else
    1017             : 
    1018       35952 :           do while (.not. (draSC0(iRow, jCol) > 0))
    1019       35040 :             prevRow = iRow
    1020       35040 :             prevCol = jCol
    1021             :             call moveDownOneCell(fDir0(iRow, jcol), iRow, jCol)
    1022             :             ! check whether this location is an outlet and exit
    1023       70056 :             do kk = 1, Noutlets
    1024       70056 :               if (iRow .eq. oLoc(kk, 1) .and. jCol .eq. oLoc(kk, 2)) exit
    1025             :             end do
    1026       35952 :             if (prevRow .eq. iRow .and. prevCol .eq. jCol) then
    1027             :               call error_message('Something went wrong during L11_link_location, ', &
    1028             :                   'movedownonecell got stuck in infinite loop at cell (', num2str(iRow), ' ', &
    1029           0 :                   num2str(jCol))
    1030             :             end if
    1031             :           end do
    1032             :           ! set "to" cell (when an outlet is reached)
    1033         912 :           nLinkToRow(ii) = iRow
    1034         912 :           nLinkToCol(ii) = jCol
    1035             : 
    1036             :         end if
    1037             :       end do
    1038             : 
    1039             :       ! end of multi-node network design loop
    1040             :     end if
    1041             : 
    1042             :     !--------------------------------------------------------
    1043             :     ! Start padding up local variables to global variables
    1044             :     !--------------------------------------------------------
    1045             :     ! L11 network data sets
    1046          24 :     call append(L11_fRow, nLinkFromRow(:))
    1047          24 :     call append(L11_fCol, nLinkFromCol(:))
    1048          24 :     call append(L11_tRow, nLinkToRow(:))
    1049          24 :     call append(L11_tCol, nLinkToCol(:))
    1050             : 
    1051             :     ! free space
    1052             :     deallocate (rowOut, colOut, nLinkFromN, netPerm, nLinkFromRow, &
    1053          24 :             nLinkFromCol, nLinkToRow, nLinkToCol, fDir0, draSC0)
    1054             : 
    1055          24 :   end subroutine L11_link_location
    1056             : 
    1057             :   ! ------------------------------------------------------------------
    1058             : 
    1059             :   !    NAME
    1060             :   !        L11_set_drain_outlet_gauges
    1061             : 
    1062             :   !    PURPOSE
    1063             :   !>       \brief Draining cell identification and Set gauging node
    1064             : 
    1065             :   !>       \details Perform the following tasks:
    1066             :   !>       - Draining cell identification (cell at L0 to draining cell outlet at L11).
    1067             :   !>       - Set gauging nodes
    1068             :   !>       If a variable is added or removed here, then it also has to
    1069             :   !>       be added or removed in the subroutine L11_config_set in
    1070             :   !>       module mo_restart and in the subroutine set_L11_config in module
    1071             :   !>       mo_set_netcdf_restart
    1072             : 
    1073             :   !    INTENT(IN)
    1074             :   !>       \param[in] "integer(i4) :: iDomain" Domain Id
    1075             : 
    1076             :   !    HISTORY
    1077             :   !>       \authors Luis Samaniego
    1078             : 
    1079             :   !>       \date Dec 2005
    1080             : 
    1081             :   ! Modifications:
    1082             :   ! Luis Samaniego Jan 2013 - modular version
    1083             :   ! Matthias Zink  Mar 2014 - bugfix, added inflow gauge
    1084             :   ! Rohini Kumar   Apr 2014 - variable index is changed to index_gauge
    1085             :   ! Stephan Thober Aug 2015 - ported to mRM
    1086             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
    1087             : 
    1088          24 :   subroutine L11_set_drain_outlet_gauges(iDomain)
    1089             : 
    1090          24 :     use mo_append, only : append
    1091             :     use mo_common_constants, only : nodata_i4
    1092             :     use mo_common_types, only: Grid
    1093             :     use mo_common_variables, only : domainMeta, level0
    1094             :     use mo_mrm_global_variables, only : L0_InflowgaugeLoc, L0_draCell, L0_draSC, L0_fDir, L0_gaugeLoc, domain_mrm, &
    1095             :                                         l0_l11_remap
    1096             : 
    1097             :     implicit none
    1098             : 
    1099             :     ! Domain Id
    1100             :     integer(i4), intent(in) :: iDomain
    1101             : 
    1102          24 :     integer(i4), dimension(:, :), allocatable :: draSC0
    1103             : 
    1104          24 :     integer(i4), dimension(:, :), allocatable :: fDir0
    1105             : 
    1106          24 :     integer(i4), dimension(:, :), allocatable :: gaugeLoc0
    1107             : 
    1108          24 :     integer(i4), dimension(:, :), allocatable :: InflowGaugeLoc0
    1109             : 
    1110          24 :     integer(i4), dimension(:, :), allocatable :: draCell0
    1111             : 
    1112             :     integer(i4) :: ii, jj, kk, ll, s0, e0
    1113             : 
    1114             :     integer(i4) :: iSc
    1115             : 
    1116             :     integer(i4) :: iRow, jCol
    1117             : 
    1118             :     type(Grid), pointer :: level0_iDomain => null()
    1119             : 
    1120             : 
    1121          24 :     level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
    1122          24 :     s0 = level0_iDomain%iStart
    1123          24 :     e0 = level0_iDomain%iEnd
    1124             : 
    1125             : 
    1126             :     ! allocate
    1127          96 :     allocate (draSC0          (level0_iDomain%nrows, level0_iDomain%ncols))
    1128          72 :     allocate (fDir0           (level0_iDomain%nrows, level0_iDomain%ncols))
    1129          72 :     allocate (gaugeLoc0       (level0_iDomain%nrows, level0_iDomain%ncols))
    1130          72 :     allocate (InflowGaugeLoc0 (level0_iDomain%nrows, level0_iDomain%ncols))
    1131          72 :     allocate (draCell0        (level0_iDomain%nrows, level0_iDomain%ncols))
    1132             : 
    1133             :     ! initialize
    1134     2862360 :     draSC0(:, :) = nodata_i4
    1135     2862360 :     fDir0(:, :) = nodata_i4
    1136     2862360 :     gaugeLoc0(:, :) = nodata_i4
    1137     2862360 :     InflowGaugeLoc0(:, :) = nodata_i4
    1138     2862360 :     draCell0(:, :) = nodata_i4
    1139             : 
    1140           0 :     draSC0(:, :) = UNPACK(L0_draSC          (s0 : e0), &
    1141          24 :             level0_iDomain%mask, nodata_i4)
    1142           0 :     fDir0(:, :) = UNPACK(L0_fDir           (s0 : e0), &
    1143          24 :             level0_iDomain%mask, nodata_i4)
    1144           0 :     gaugeLoc0(:, :) = UNPACK(L0_gaugeLoc       (s0 : e0), &
    1145          24 :             level0_iDomain%mask, nodata_i4)
    1146           0 :     InflowGaugeLoc0(:, :) = UNPACK(L0_InflowgaugeLoc (s0 : e0), &
    1147          24 :             level0_iDomain%mask, nodata_i4)
    1148             : 
    1149     1066064 :     do kk = 1, level0_iDomain%nCells
    1150     1066040 :       ii = level0_iDomain%CellCoor(kk, 1)
    1151     1066040 :       jj = level0_iDomain%CellCoor(kk, 2)
    1152     1066040 :       iSc = draSC0(ii, jj)
    1153             :       ! find drainage path
    1154     1066040 :       iRow = ii
    1155     1066040 :       jCol = jj
    1156    43644186 :       do while (.NOT. iSC > 0)
    1157             :         ! move downstream
    1158    42578146 :         call moveDownOneCell(fDir0(iRow, jCol), iRow, jCol)
    1159    42578146 :         iSC = draSC0(iRow, jCol)
    1160             :       end do
    1161     1066040 :       draCell0(ii, jj) = iSC
    1162             : 
    1163             :       ! find cell at L11 corresponding to gauges in Domain at L0 ! >> L11_on_L0 is Id of
    1164             :       ! the routing cell at level-11
    1165     1066040 :       if (gaugeLoc0(ii, jj) .NE. nodata_i4) then
    1166             :         ! evaluation gauges
    1167          76 :         do ll = 1, domain_mrm(iDomain)%nGauges
    1168             :           ! search for gaugeID in L0 grid and save ID on L11
    1169          76 :           if (domain_mrm(iDomain)%gaugeIdList(ll) .EQ. gaugeLoc0(ii, jj)) then
    1170          21 :             domain_mrm(iDomain)%gaugeNodeList(ll) = L0_L11_remap(iDomain)%lowres_id_on_highres(ii, jj)
    1171             :           end if
    1172             :         end do
    1173             :       end if
    1174             : 
    1175     1066064 :       if (InflowGaugeLoc0(ii, jj) .NE. nodata_i4) then
    1176             :         ! inflow gauges
    1177          39 :         do ll = 1, domain_mrm(iDomain)%nInflowGauges
    1178             :           ! search for gaugeID in L0 grid and save ID on L11
    1179           4 :           if (domain_mrm(iDomain)%InflowGaugeIdList(ll) .EQ. InflowGaugeLoc0(ii, jj)) &
    1180          37 :               domain_mrm(iDomain)%InflowGaugeNodeList(ll) = L0_L11_remap(iDomain)%lowres_id_on_highres(ii, jj)
    1181             :         end do
    1182             :       end if
    1183             :     end do
    1184             : 
    1185             :     !--------------------------------------------------------
    1186             :     ! Start padding up local variables to global variables
    1187             :     !--------------------------------------------------------
    1188             :     ! L0 data sets
    1189             :     ! check whether L0 data is shared
    1190             :     ! ToDo: check if change is correct
    1191          24 :     if (iDomain .eq. 1) then
    1192          12 :       call append(L0_draCell, PACK(draCell0(:, :), level0_iDomain%mask))
    1193          12 :     else if (domainMeta%L0DataFrom(iDomain) == iDomain) then
    1194           9 :       call append(L0_draCell, PACK(draCell0(:, :), level0_iDomain%mask))
    1195             :     end if
    1196             : 
    1197             :     ! free space
    1198          24 :     deallocate (draSC0, fDir0, gaugeLoc0, draCell0)
    1199             : 
    1200          24 :   end subroutine  L11_set_drain_outlet_gauges
    1201             : 
    1202             :   ! ------------------------------------------------------------------
    1203             : 
    1204             :   !    NAME
    1205             :   !        L11_stream_features
    1206             : 
    1207             :   !    PURPOSE
    1208             :   !>       \brief Stream features (stream network and floodplain)
    1209             : 
    1210             :   !>       \details Stream features (stream network and floodplain)
    1211             :   !>       If a variable is added or removed here, then it also has to
    1212             :   !>       be added or removed in the subroutine L11_config_set in
    1213             :   !>       module mo_restart and in the subroutine set_L11_config in module
    1214             :   !>       mo_set_netcdf_restart
    1215             : 
    1216             :   !    INTENT(IN)
    1217             :   !>       \param[in] "integer(i4) :: iDomain" Domain Id
    1218             : 
    1219             :   !    HISTORY
    1220             :   !>       \authors Luis Samaniego
    1221             : 
    1222             :   !>       \date Dec 2005
    1223             : 
    1224             :   ! Modifications:
    1225             :   ! Luis Samaniego Jan 2013 - modular version
    1226             :   ! R. Kumar       Oct 2013 - stack size increased from nNodes to 100
    1227             :   ! Stephan Thober Aug 2015 - ported to mRM
    1228             :   ! Stephan Thober Nov 2016 - only read flood plain area if processMatrix for routing equals 1
    1229             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
    1230             :   ! Stephan Thober Jul 2018 - introduced cut off Length at 40 percentile to neglect short paths in headwaters for adaptive timesteps
    1231             :   ! Stephan Thober, Pallav Kumar Shrestha, Sep 2020 - bug fix in cut off Length at 40 percentile, neglecting links with -9999. that occur if multiple outlets are present
    1232             : 
    1233          24 :   subroutine L11_stream_features(iDomain)
    1234             : 
    1235          24 :     use mo_append, only : append
    1236             :     use mo_common_constants, only : nodata_dp, nodata_i4
    1237             :     use mo_common_types, only: Grid
    1238             :     use mo_common_variables, only : domainMeta, L0_elev, iFlag_cordinate_sys, level0, processMatrix
    1239             :     use mo_mrm_global_variables, only : L0_fDir, &
    1240             :         L0_floodPlain, L0_streamNet, L11_aFloodPlain, L11_fCol, L11_fRow, L11_length, &
    1241             :         L11_nOutlets, L11_netPerm, L11_slope, L11_tCol, L11_tRow, level11
    1242             :     use mo_percentile, only: percentile
    1243             : 
    1244             :     implicit none
    1245             : 
    1246             :     ! Domain Id
    1247             :     integer(i4), intent(in) :: iDomain
    1248             : 
    1249             :     integer(i4) :: nLinks
    1250             : 
    1251          24 :     integer(i4), dimension(:, :), allocatable :: iD0
    1252             : 
    1253          24 :     integer(i4), dimension(:, :), allocatable :: fDir0
    1254             : 
    1255          24 :     real(dp), dimension(:, :), allocatable :: elev0
    1256             : 
    1257          24 :     real(dp), dimension(:, :), allocatable :: cellarea0
    1258             : 
    1259          24 :     integer(i4), dimension(:, :), allocatable :: streamNet0
    1260             : 
    1261          24 :     integer(i4), dimension(:, :), allocatable :: floodPlain0
    1262             : 
    1263             :     ! routing order (permutation)
    1264          24 :     integer(i4), dimension(:), allocatable :: netPerm
    1265             : 
    1266          24 :     integer(i4), dimension(:), allocatable :: nLinkFromRow
    1267             : 
    1268          24 :     integer(i4), dimension(:), allocatable :: nLinkFromCol
    1269             : 
    1270          24 :     integer(i4), dimension(:), allocatable :: nLinkToRow
    1271             : 
    1272          24 :     integer(i4), dimension(:), allocatable :: nLinkToCol
    1273             : 
    1274          24 :     real(dp), dimension(:), allocatable :: nLinkLength
    1275             : 
    1276          24 :     real(dp), dimension(:), allocatable :: nLinkAFloodPlain
    1277             : 
    1278          24 :     real(dp), dimension(:), allocatable :: nLinkSlope
    1279             : 
    1280             :     integer(i4) :: ii, rr, ns, s0, e0
    1281             : 
    1282             :     integer(i4) :: frow, fcol
    1283             : 
    1284             :     integer(i4) :: fId, tId
    1285             : 
    1286          24 :     integer(i4), dimension(:, :), allocatable :: stack, append_chunk
    1287             : 
    1288          24 :     integer(i4), dimension(:), allocatable :: dummy_1d
    1289             : 
    1290          24 :     real(dp) :: length
    1291             : 
    1292          24 :     integer(i4), dimension(:, :), allocatable :: nodata_i4_tmp
    1293             : 
    1294          24 :     real(dp), dimension(:, :), allocatable :: nodata_dp_tmp
    1295             : 
    1296             :     type(Grid), pointer :: level0_iDomain => null()
    1297             : 
    1298             : 
    1299          24 :     level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
    1300          24 :     s0 = level0_iDomain%iStart
    1301          24 :     e0 = level0_iDomain%iEnd
    1302          24 :     nLinks = level11(iDomain)%nCells - L11_nOutlets(iDomain)
    1303             : 
    1304             : 
    1305             :     ! allocate
    1306          96 :     allocate (iD0           (level0_iDomain%nrows, level0_iDomain%ncols))
    1307          96 :     allocate (elev0         (level0_iDomain%nrows, level0_iDomain%ncols))
    1308          96 :     allocate (fDir0         (level0_iDomain%nrows, level0_iDomain%ncols))
    1309          96 :     allocate (cellarea0     (level0_iDomain%nrows, level0_iDomain%ncols))
    1310          96 :     allocate (streamNet0    (level0_iDomain%nrows, level0_iDomain%ncols))
    1311          72 :     allocate (floodPlain0   (level0_iDomain%nrows, level0_iDomain%ncols))
    1312             : 
    1313             :     !  Routing network vectors have nNodes size instead of nLinks to
    1314             :     !  avoid the need of having two extra indices to identify a Domain.
    1315          96 :     allocate (stack             (level11(iDomain)%nCells, 2)) ! >> stack(nNodes, 2)
    1316          24 :     allocate (dummy_1d          (2))
    1317          24 :     allocate (append_chunk      (8, 2))
    1318          72 :     allocate (netPerm           (level11(iDomain)%nCells))
    1319          48 :     allocate (nLinkFromRow      (level11(iDomain)%nCells))
    1320          48 :     allocate (nLinkFromCol      (level11(iDomain)%nCells))
    1321          48 :     allocate (nLinkToRow        (level11(iDomain)%nCells))
    1322          48 :     allocate (nLinkToCol        (level11(iDomain)%nCells))
    1323          72 :     allocate (nLinkLength       (level11(iDomain)%nCells))
    1324          48 :     allocate (nLinkAFloodPlain  (level11(iDomain)%nCells))
    1325          48 :     allocate (nLinkSlope        (level11(iDomain)%nCells))
    1326             : 
    1327          96 :     allocate (nodata_i4_tmp      (level0_iDomain%nrows, level0_iDomain%ncols))
    1328          96 :     allocate (nodata_dp_tmp      (level0_iDomain%nrows, level0_iDomain%ncols))
    1329             : 
    1330             :     ! initialize
    1331     2862360 :     iD0(:, :) = nodata_i4
    1332     2862360 :     elev0(:, :) = nodata_dp
    1333     2862360 :     fDir0(:, :) = nodata_i4
    1334     2862360 :     cellarea0(:, :) = nodata_dp
    1335     2862360 :     streamNet0(:, :) = nodata_i4
    1336     2862360 :     floodPlain0(:, :) = nodata_i4
    1337             : 
    1338        1944 :     stack(:, :) = nodata_i4
    1339         456 :     append_chunk(:, :) = nodata_i4
    1340         960 :     netPerm(:) = nodata_i4
    1341         960 :     nLinkFromRow(:) = nodata_i4
    1342         960 :     nLinkFromCol(:) = nodata_i4
    1343         960 :     nLinkToRow(:) = nodata_i4
    1344         960 :     nLinkToCol(:) = nodata_i4
    1345         960 :     nLinkLength(:) = nodata_dp
    1346         960 :     nLinkAFloodPlain(:) = nodata_dp
    1347         960 :     nLinkSlope(:) = nodata_dp
    1348             : 
    1349     2862360 :     nodata_i4_tmp(:, :) = nodata_i4
    1350     2862360 :     nodata_dp_tmp(:, :) = nodata_dp
    1351             : 
    1352             :     ! for a single node model run
    1353          24 :     if(level11(iDomain)%nCells .GT. 1) then
    1354             :       ! get L0 fields
    1355          24 :       iD0(:, :) = UNPACK(level0_iDomain%Id, level0_iDomain%mask, nodata_i4_tmp)
    1356           0 :       elev0(:, :) = UNPACK(L0_elev (s0 : e0), &
    1357          24 :           level0_iDomain%mask, nodata_dp_tmp)
    1358           0 :       fDir0(:, :) = UNPACK(L0_fDir (s0 : e0), &
    1359          24 :           level0_iDomain%mask, nodata_i4_tmp)
    1360          24 :       cellarea0(:, :) = UNPACK(level0_iDomain%CellArea, level0_iDomain%mask, nodata_dp_tmp)
    1361             : 
    1362             :       ! get network vectors of L11
    1363         960 :       netPerm(:) = L11_netPerm (level11(iDomain)%iStart : level11(iDomain)%iEnd)
    1364         960 :       nLinkFromRow(:) = L11_fRow    (level11(iDomain)%iStart : level11(iDomain)%iEnd)
    1365         960 :       nLinkFromCol(:) = L11_fCol    (level11(iDomain)%iStart : level11(iDomain)%iEnd)
    1366         960 :       nLinkToRow(:) = L11_tRow    (level11(iDomain)%iStart : level11(iDomain)%iEnd)
    1367         960 :       nLinkToCol(:) = L11_tCol    (level11(iDomain)%iStart : level11(iDomain)%iEnd)
    1368             : 
    1369             :       ! Flood plains:  stream network delineation
    1370     2862360 :       streamNet0(:, :) = nodata_i4
    1371     2862360 :       floodPlain0(:, :) = nodata_i4
    1372             : 
    1373         936 :       do rr = 1, nLinks
    1374             : 
    1375         912 :         ii = netPerm(rr)
    1376         912 :         frow = nLinkFromRow(ii)
    1377         912 :         fcol = nLinkFromCol(ii)
    1378             : 
    1379             :         ! Init
    1380         912 :         streamNet0(frow, fcol) = ii
    1381         912 :         floodPlain0(frow, fcol) = ii
    1382      144024 :         stack = 0
    1383       18240 :         append_chunk = 0
    1384         912 :         ns = 1
    1385         912 :         stack(ns, 1) = frow
    1386         912 :         stack(ns, 2) = fcol
    1387             : 
    1388         912 :         call cellLength(iDomain, fDir0(frow, fcol), fRow, fCol, iFlag_cordinate_sys, nLinkLength(ii))
    1389         912 :         nLinkSlope(ii) = elev0(frow, fcol)
    1390             : 
    1391         912 :         fId = iD0(frow, fcol)
    1392         912 :         tId = iD0(nLinkToRow(ii), nLinkToCol(ii))
    1393             : 
    1394       35952 :         do while (.NOT. (fId == tId))
    1395             :           ! Search flood plain from point(frow,fcol) upwards, keep co-ordinates in STACK
    1396     6659002 :           do while (ns > 0)
    1397     6623962 :             if (ns + 8 .gt. size(stack, 1)) then
    1398         118 :               call append(stack, append_chunk)
    1399             :             end if
    1400     6623962 :             call moveUp(elev0, fDir0, frow, fcol, stack, ns)
    1401     6623962 :             stack(1, 1) = 0
    1402     6623962 :             stack(1, 2) = 0
    1403             :             ! stack = cshift(stack, SHIFT = 1, DIM = 1)
    1404             :             ! substitute cshift <<<
    1405    26495848 :             dummy_1d = stack(1, :)
    1406  2009564780 :             stack(: size(stack, dim = 1) - 1, :) = stack(2 :, :)
    1407    19871886 :             stack(size(stack, dim = 1), :) = dummy_1d
    1408             :             ! substitute cshift >>>
    1409     6623962 :             if (stack(1, 1) > 0 .and. stack(1, 2) > 0) floodPlain0(stack(1, 1), stack(1, 2)) = ii
    1410  1018065354 :             ns = count(stack > 0) / 2
    1411             :           end do
    1412             : 
    1413             :           ! move downstream
    1414       35040 :           call moveDownOneCell(fDir0(frow, fcol), frow, fcol)
    1415       35040 :           streamNet0(frow, fcol) = ii
    1416       35040 :           floodPlain0(frow, fcol) = ii
    1417       35040 :           fId = iD0(frow, fcol)
    1418     5303016 :           stack = 0
    1419       35040 :           ns = 1
    1420       35040 :           stack(ns, 1) = frow
    1421       35040 :           stack(ns, 2) = fcol
    1422       35040 :           call cellLength(iDomain, fDir0(fRow, fCol), fRow, fCol, iFlag_cordinate_sys, length)
    1423       35040 :           nLinkLength(ii) = nLinkLength(ii) + length
    1424             : 
    1425             :         end do
    1426             : 
    1427             :         ! stream bed slope
    1428         912 :         nLinkSlope(ii) = (nLinkSlope(ii) - elev0(frow, fcol)) / nLinkLength(ii)
    1429             : 
    1430         912 :         if (nLinkSlope(ii) < 0.0001_dp) nLinkSlope(ii) = 0.0001_dp
    1431             : 
    1432             :         ! calculate area of floodplains (avoid overwriting)
    1433   111450936 :         nLinkAFloodPlain(ii) = sum (cellarea0(:, :), mask = (floodPlain0(:, :) == ii))
    1434             :         !  old > real( count( floodPlain0(:,:,) == i), dp ) * cellarea0
    1435             : 
    1436             :       end do
    1437             : 
    1438             :       ! end of multi-node network design loop
    1439             :     end if
    1440             : 
    1441             :     ! cut off Length at 40 percentile to neglect short paths in headwaters
    1442          24 :     if ((processMatrix(8, 1) .eq. 2) .or. (processMatrix(8, 1) .eq. 3)) then
    1443         280 :       if (count(nLinkLength(:) .ge. 0._dp) .gt. 2) then
    1444         280 :         length = percentile(pack(nLinkLength(:), nLinkLength(:) .ge. 0._dp), 40._dp)
    1445         280 :         nLinkLength(:) = merge(nLinkLength(:), length, (nLinkLength(:) .gt. length))
    1446             :       end if
    1447             :     end if
    1448             : 
    1449             :     !--------------------------------------------------------
    1450             :     ! Start padding up local variables to global variables
    1451             :     !--------------------------------------------------------
    1452             : 
    1453             :     ! L0 data sets
    1454             :     ! check whether L0 data is shared
    1455             :     ! ToDo: check if change is correct
    1456          24 :     if (iDomain .eq. 1) then
    1457          12 :       call append(L0_streamNet, PACK (streamNet0(:, :), level0_iDomain%mask))
    1458          12 :       call append(L0_floodPlain, PACK (floodPlain0(:, :), level0_iDomain%mask))
    1459          12 :     else if (domainMeta%L0DataFrom(iDomain) == iDomain) then
    1460           9 :       call append(L0_streamNet, PACK (streamNet0(:, :), level0_iDomain%mask))
    1461           9 :       call append(L0_floodPlain, PACK (floodPlain0(:, :), level0_iDomain%mask))
    1462             :     end if
    1463             : 
    1464             : 
    1465             :     ! L11 network data sets
    1466          24 :     call append(L11_length, nLinkLength(:))
    1467          24 :     call append(L11_aFloodPlain, nLinkAFloodPlain(:))
    1468          24 :     call append(L11_slope, nLinkSlope(:))
    1469             : 
    1470             :     ! free space
    1471           0 :     deallocate (&
    1472           0 :         iD0, elev0, fDir0, streamNet0, floodPlain0, &
    1473           0 :         cellarea0, stack, netPerm, nLinkFromRow, nLinkFromCol, nLinkToRow, nLinkToCol, &
    1474          24 :         nLinkLength, nLinkAFloodPlain, nLinkSlope, dummy_1d)
    1475          24 :     deallocate(nodata_i4_tmp, nodata_dp_tmp)
    1476             : 
    1477          24 :   end subroutine L11_stream_features
    1478             : 
    1479             :   ! ------------------------------------------------------------------
    1480             : 
    1481             :   !    NAME
    1482             :   !        L11_fraction_sealed_floodplain
    1483             : 
    1484             :   !    PURPOSE
    1485             :   !>       \brief Fraction of the flood plain with impervious cover
    1486             : 
    1487             :   !>       \details Fraction of the flood plain with impervious cover (\ref fig_routing "Routing
    1488             :   !>       Network"). This proportion is used to regionalize the Muskingum parameters.
    1489             :   !>       Samaniego et al. \cite SB05 found out that this fraction is one of the statistically
    1490             :   !>       significant predictor variables of peak discharge in mesoscale Domains.
    1491             :   !>       If a variable is added or removed here, then it also has to
    1492             :   !>       be added or removed in the subroutine L11_config_set in
    1493             :   !>       module mo_restart and in the subroutine set_L11_config in module
    1494             :   !>       mo_set_netcdf_restart
    1495             : 
    1496             :   !    INTENT(IN)
    1497             :   !>       \param[in] "integer(i4) :: LCClassImp" Impervious land cover class Id, e.g. = 2 (old code)
    1498             :   !>       \param[in] "logical :: do_init"
    1499             : 
    1500             :   !    HISTORY
    1501             :   !>       \authors Luis Samaniego
    1502             : 
    1503             :   !>       \date Dec 2005
    1504             : 
    1505             :   ! Modifications:
    1506             :   ! Luis Samaniego Jan 2013 - modular version
    1507             :   ! Stephan Thober Aug 2015 - ported to mRM
    1508             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
    1509             : 
    1510          13 :   subroutine L11_fraction_sealed_floodplain(LCClassImp, do_init)
    1511             : 
    1512          24 :     use mo_append, only : append
    1513             :     use mo_common_constants, only : nodata_dp
    1514             :     use mo_common_types, only: Grid
    1515             :     use mo_common_variables, only : domainMeta, L0_LCover, level0, domainMeta, nLCoverScene
    1516             :     use mo_mrm_global_variables, only : L0_floodPlain, L11_aFloodPlain, &
    1517             :                                         L11_nLinkFracFPimp, L11_nOutlets, level11
    1518             : 
    1519             :     implicit none
    1520             : 
    1521             :     ! Impervious land cover class Id, e.g. = 2 (old code)
    1522             :     integer(i4), intent(in) :: LCClassImp
    1523             : 
    1524             :     logical, intent(in) :: do_init
    1525             : 
    1526             :     integer(i4) :: nLinks
    1527             : 
    1528             :     real(dp), dimension(:), pointer :: nLinkAFloodPlain => null()
    1529             : 
    1530          13 :     real(dp), dimension(:,:), allocatable :: temp_array
    1531             : 
    1532             :     integer(i4) :: ii, iDomain, iiLC, s0, e0
    1533             : 
    1534             :     type(Grid), pointer :: level0_iDomain => null()
    1535             : 
    1536             : 
    1537             :     ! initialization
    1538          38 :     do iDomain = 1, domainMeta%nDomains
    1539         100 :       allocate(temp_array(level11(iDomain)%nCells, nLCoverScene))
    1540        2040 :       temp_array = nodata_dp
    1541          25 :       if (do_init) then
    1542          16 :         level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
    1543             : 
    1544          16 :         s0 = level0_iDomain%iStart
    1545          16 :         e0 = level0_iDomain%iEnd
    1546          16 :         nLinks = level11(iDomain)%nCells + 1 - L11_nOutlets(iDomain)
    1547          16 :         nLinkAFloodPlain => L11_aFloodPlain(level11(iDomain)%iStart : level11(iDomain)%iEnd)
    1548             : 
    1549          48 :         do iiLC = 1, nLCoverScene
    1550             :           ! for a single node model run
    1551          48 :           if(nLinks .GT. 0) then
    1552        1360 :             do ii = 1, nLinks
    1553        2656 :               temp_array(ii, iiLC) = sum(level0_iDomain%CellArea(:), &
    1554           0 :                   mask = (L0_floodPlain(s0 : e0) == ii .and. L0_LCover(s0 : e0, iiLC) == LCClassImp)) &
    1555    59876256 :                   / nLinkAFloodPlain(ii)
    1556             :             end do
    1557             :           end if
    1558             :         end do
    1559             :       end if
    1560          25 :       call append(L11_nLinkFracFPimp, temp_array(:,:))
    1561          38 :       deallocate(temp_array)
    1562             :     end do
    1563             : 
    1564          13 :   end subroutine L11_fraction_sealed_floodplain
    1565             : 
    1566             :   ! ------------------------------------------------------------------
    1567             :   !  MOVE UPSTREAM FROM-TO
    1568             :   ! ------------------------------------------------------------------
    1569             :   !    NAME
    1570             :   !        moveUp
    1571             : 
    1572             :   !    PURPOSE
    1573             :   !>       \brief TODO: add description
    1574             : 
    1575             :   !>       \details TODO: add description
    1576             : 
    1577             :   !    INTENT(IN)
    1578             :   !>       \param[in] "real(dp), dimension(:, :) :: elev0"
    1579             :   !>       \param[in] "integer(i4), dimension(:, :) :: fDir0"
    1580             :   !>       \param[in] "integer(i4) :: fi, fj"                 co-ordinate of the stream bed
    1581             :   !>       \param[in] "integer(i4) :: fi, fj"                 co-ordinate of the stream bed
    1582             : 
    1583             :   !    INTENT(INOUT)
    1584             :   !>       \param[inout] "integer(i4), dimension(:, :) :: ss"
    1585             :   !>       \param[inout] "integer(i4) :: nn"
    1586             : 
    1587             :   !    HISTORY
    1588             :   !>       \authors Robert Schweppe
    1589             : 
    1590             :   !>       \date Jun 2018
    1591             : 
    1592             :   ! Modifications:
    1593             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
    1594             : 
    1595     6623962 :   subroutine moveUp(elev0, fDir0, fi, fj, ss, nn)
    1596             : 
    1597          13 :     use mo_mrm_constants, only : deltaH
    1598             :     use mo_utils, only : ge, le
    1599             : 
    1600             :     implicit none
    1601             : 
    1602             :     real(dp), dimension(:, :), allocatable, intent(IN) :: elev0
    1603             : 
    1604             :     integer(i4), dimension(:, :), allocatable, intent(IN) :: fDir0
    1605             : 
    1606             :     ! co-ordinate of the stream bed
    1607             :     integer(i4), intent(IN) :: fi, fj
    1608             : 
    1609             :     integer(i4), dimension(:, :), intent(INOUT) :: ss
    1610             : 
    1611             :     integer(i4), intent(INOUT) :: nn
    1612             : 
    1613             :     integer(i4) :: ii, jj, ip, im, jp, jm
    1614             : 
    1615             :     integer(i4) :: nrows, ncols
    1616             : 
    1617             : 
    1618     6623962 :     ii = ss(1, 1)
    1619     6623962 :     jj = ss(1, 2)
    1620     6623962 :     ip = ii + 1
    1621     6623962 :     im = ii - 1
    1622     6623962 :     jp = jj + 1
    1623     6623962 :     jm = jj - 1
    1624             : 
    1625     6623962 :     nrows = size(fDir0, 1)
    1626     6623962 :     ncols = size(fDir0, 2)
    1627             : 
    1628             :     !E
    1629     6623962 :     if   (jp                <= ncols) then
    1630    13247924 :       if ((fdir0(ii, jp) == 16)                   .and. &
    1631    26495848 :           (le((elev0(ii, jp) - elev0(fi, fj)), deltaH))  .and. &
    1632     6623962 :           (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp))        &
    1633             :           ) then
    1634     1415524 :         nn = nn + 1
    1635     1415524 :         ss(nn, 1) = ii
    1636     1415524 :         ss(nn, 2) = jp
    1637             :         !print *, i,jp
    1638             :       end if
    1639             :     end if
    1640             : 
    1641             :     !SE
    1642     6623962 :     if ((ip                <= nrows) .and. &
    1643             :         (jp                <= ncols)) then
    1644           0 :       if ((fdir0(ip, jp) == 32)                   .and. &
    1645    13247924 :           (le((elev0(ip, jp) - elev0(fi, fj)), deltaH))  .and. &
    1646     6623962 :           (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp))        &
    1647             :           ) then
    1648      470784 :         nn = nn + 1
    1649      470784 :         ss(nn, 1) = ip
    1650      470784 :         ss(nn, 2) = jp
    1651             :         !print *, ip,jp
    1652             :       end if
    1653             :     end if
    1654             : 
    1655             :     !S
    1656     6623962 :     if ((ip               <= nrows)  .and. &
    1657             :         (jp               <= ncols)) then
    1658           0 :       if ((fdir0(ip, jj) == 64)                 .and. &
    1659    13247924 :           (le((elev0(ip, jj) - elev0(fi, fj)), deltaH))  .and. &
    1660     6623962 :           (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp))        &
    1661             :           ) then
    1662     1125660 :         nn = nn + 1
    1663     1125660 :         ss(nn, 1) = ip
    1664     1125660 :         ss(nn, 2) = jj
    1665             :         !print *, ip,j
    1666             :       end if
    1667             :     end if
    1668             : 
    1669             :     !SW
    1670             :     if ((ip                <= nrows) .and. &
    1671     6623962 :         (jp                <= ncols) .and. &
    1672             :         (jm                >= 1)) then
    1673           0 :       if ((fdir0(ip, jm) == 128)                 .and. &
    1674    13247924 :           (le((elev0(ip, jm) - elev0(fi, fj)), deltaH))  .and. &
    1675     6623962 :           (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp))        &
    1676             :           ) then
    1677      335846 :         nn = nn + 1
    1678      335846 :         ss(nn, 1) = ip
    1679      335846 :         ss(nn, 2) = jm
    1680             :         !print *, ip,jm
    1681             :       end if
    1682             :     end if
    1683             : 
    1684             :     !W
    1685     6623962 :     if ((jm                 >= 1) .and. &
    1686             :         (jp                 <= ncols)) then
    1687           0 :       if ((fdir0(ii, jm)  == 1)                 .and. &
    1688    13247924 :           (le((elev0(ii, jm) - elev0(fi, fj)), deltaH))  .and. &
    1689     6623962 :           (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp))        &
    1690             :           ) then
    1691      912382 :         nn = nn + 1
    1692      912382 :         ss(nn, 1) = ii
    1693      912382 :         ss(nn, 2) = jm
    1694             :         !print *, i,jm
    1695             :       end if
    1696             :     end if
    1697             : 
    1698             :     !NW
    1699             :     if ((im                >= 1) .and. &
    1700     6623962 :         (jp                <= ncols) .and. &
    1701             :         (jm                >= 1))  then
    1702           0 :       if ((fdir0(im, jm) == 2)                 .and. &
    1703    13247924 :           (le((elev0(im, jm) - elev0(fi, fj)), deltaH))  .and. &
    1704     6623962 :           (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp))        &
    1705             :           ) then
    1706      492132 :         nn = nn + 1
    1707      492132 :         ss(nn, 1) = im
    1708      492132 :         ss(nn, 2) = jm
    1709             :         !print *, im,jm
    1710             :       end if
    1711             :     end if
    1712             : 
    1713             :     !N
    1714     6623962 :     if ((im                >= 1) .and. &
    1715             :         (jp                 <= ncols)) then
    1716           0 :       if ((fdir0(im, jj)  == 4)                 .and. &
    1717    13247924 :           (le((elev0(im, jj) - elev0(fi, fj)), deltaH))  .and. &
    1718     6623962 :           (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp))        &
    1719             :           ) then
    1720     1179748 :         nn = nn + 1
    1721     1179748 :         ss(nn, 1) = im
    1722     1179748 :         ss(nn, 2) = jj
    1723             :         !print *, im,j
    1724             :       end if
    1725             :     end if
    1726             : 
    1727             :     !NE
    1728     6623962 :     if ((im                >= 1) .and. &
    1729             :         (jp                <= ncols))  then
    1730           0 :       if ((fdir0(im, jp) == 8)               .and. &
    1731    13247924 :           (le((elev0(im, jp) - elev0(fi, fj)), deltaH))  .and. &
    1732     6623962 :           (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp))        &
    1733             :           ) then
    1734      656846 :         nn = nn + 1
    1735      656846 :         ss(nn, 1) = im
    1736      656846 :         ss(nn, 2) = jp
    1737             :         !print *, im,jp
    1738             :       end if
    1739             :     end if
    1740             : 
    1741     6623962 :   end subroutine moveUp
    1742             : 
    1743             :   ! ------------------------------------------------------------------
    1744             :   !  MOVE DOWNSTREAM
    1745             :   ! ------------------------------------------------------------------
    1746             :   !    NAME
    1747             :   !        moveDownOneCell
    1748             : 
    1749             :   !    PURPOSE
    1750             :   !>       \brief TODO: add description
    1751             : 
    1752             :   !>       \details TODO: add description
    1753             : 
    1754             :   !    INTENT(IN)
    1755             :   !>       \param[in] "integer(i4) :: fDir"
    1756             : 
    1757             :   !    INTENT(INOUT)
    1758             :   !>       \param[inout] "integer(i4) :: iRow, jCol"
    1759             :   !>       \param[inout] "integer(i4) :: iRow, jCol"
    1760             : 
    1761             :   !    HISTORY
    1762             :   !>       \authors Robert Schweppe
    1763             : 
    1764             :   !>       \date Jun 2018
    1765             : 
    1766             :   ! Modifications:
    1767             : 
    1768    43740628 :   subroutine moveDownOneCell(fDir, iRow, jCol)
    1769             :     implicit none
    1770             : 
    1771             :     integer(i4), intent(IN) :: fDir
    1772             : 
    1773             :     integer(i4), intent(INOUT) :: iRow, jCol
    1774             : 
    1775             : 
    1776    45179838 :     select case (fDir)
    1777             :     case(1)   !E
    1778     1439210 :       jCol = jCol + 1
    1779             :     case(2)   !SE
    1780     1682387 :       iRow = iRow + 1
    1781     1682387 :       jCol = jCol + 1
    1782             :     case(4)   !S
    1783     3804397 :       iRow = iRow + 1
    1784             :     case(8)   !SW
    1785     5753487 :       iRow = iRow + 1
    1786     5753487 :       jCol = jCol - 1
    1787             :     case(16)  !W
    1788    11095532 :       jCol = jCol - 1
    1789             :     case(32)  !NW
    1790    10002909 :       iRow = iRow - 1
    1791    10002909 :       jCol = jCol - 1
    1792             :     case(64)  !N
    1793     7261091 :       iRow = iRow - 1
    1794             :     case(128) !NE
    1795     2701591 :       iRow = iRow - 1
    1796     2701591 :       jCol = jCol + 1
    1797             :     case default !sink
    1798             :       ! do nothing
    1799             :     end select
    1800             : 
    1801     6623962 :   end subroutine moveDownOneCell
    1802             : 
    1803             :   ! ------------------------------------------------------------------
    1804             :   !  CELL LENGTH
    1805             :   ! ------------------------------------------------------------------
    1806             :   !    NAME
    1807             :   !        cellLength
    1808             : 
    1809             :   !    PURPOSE
    1810             :   !>       \brief TODO: add description
    1811             : 
    1812             :   !>       \details TODO: add description
    1813             : 
    1814             :   !    INTENT(IN)
    1815             :   !>       \param[in] "integer(i4) :: iDomain"
    1816             :   !>       \param[in] "integer(i4) :: fDir"
    1817             :   !>       \param[in] "integer(i4) :: iRow"
    1818             :   !>       \param[in] "integer(i4) :: jCol"
    1819             :   !>       \param[in] "integer(i4) :: iCoorSystem"
    1820             : 
    1821             :   !    INTENT(OUT)
    1822             :   !>       \param[out] "real(dp) :: length"
    1823             : 
    1824             :   !    HISTORY
    1825             :   !>       \authors Robert Schweppe
    1826             : 
    1827             :   !>       \date Jun 2018
    1828             : 
    1829             :   ! Modifications:
    1830             : 
    1831       35952 :   subroutine cellLength(iDomain, fDir, iRow, jCol, iCoorSystem, length)
    1832             : 
    1833    43740628 :     use mo_common_types, only: Grid
    1834             :     use mo_common_variables, only : domainMeta, level0
    1835             :     use mo_constants, only : SQRT2_dp
    1836             : 
    1837             :     implicit none
    1838             : 
    1839             :     integer(i4), intent(IN) :: iDomain
    1840             : 
    1841             :     integer(i4), intent(IN) :: fDir
    1842             : 
    1843             :     integer(i4), intent(IN) :: iRow
    1844             : 
    1845             :     integer(i4), intent(IN) :: jCol
    1846             : 
    1847             :     integer(i4), intent(IN) :: iCoorSystem
    1848             : 
    1849             :     real(dp), intent(OUT) :: length
    1850             : 
    1851             :     integer(i4) :: iRow_to, jCol_to
    1852             : 
    1853       35952 :     real(dp) :: lat_1, long_1, lat_2, long_2
    1854             : 
    1855             :     type(Grid), pointer :: level0_iDomain => null()
    1856             : 
    1857             : 
    1858       35952 :     level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
    1859             : 
    1860             :     ! regular X-Y cordinate system
    1861       35952 :     IF(iCoorSystem .EQ. 0) THEN
    1862             : 
    1863       54786 :       select case (fDir)
    1864             :       case(1, 4, 16, 64)       ! E, S, W, N
    1865       18834 :         length = 1.0_dp
    1866             :       case(2, 8, 32, 128)      ! SE, SW, NW, NE
    1867       35952 :         length = SQRT2_dp
    1868             :       end select
    1869       35952 :       length = length * level0_iDomain%cellsize
    1870             : 
    1871             :       ! regular lat-lon cordinate system
    1872           0 :     ELSE IF(iCoorSystem .EQ. 1) THEN
    1873           0 :       iRow_to = iRow
    1874           0 :       jCol_to = jCol
    1875             : 
    1876             :       ! move in the direction of flow
    1877           0 :       call moveDownOneCell(fDir, iRow_to, jCol_to)
    1878             : 
    1879             :       ! estimate lat-lon points
    1880             :       lat_1 = level0_iDomain%yllcorner + real((level0_iDomain%ncols - jCol), dp) * level0_iDomain%cellsize + &
    1881           0 :               0.5_dp * level0_iDomain%cellsize
    1882             :       long_1 = level0_iDomain%xllcorner + real((iRow - 1), dp) * level0_iDomain%cellsize + &
    1883           0 :               0.5_dp * level0_iDomain%cellsize
    1884             : 
    1885             :       lat_2 = level0_iDomain%yllcorner + real((level0_iDomain%ncols - jCol_to), dp) * level0_iDomain%cellsize + &
    1886           0 :               0.5_dp * level0_iDomain%cellsize
    1887             :       long_2 = level0_iDomain%xllcorner + real((iRow_to - 1), dp) * level0_iDomain%cellsize + &
    1888           0 :               0.5_dp * level0_iDomain%cellsize
    1889             :       ! get distance between two points
    1890           0 :       call get_distance_two_lat_lon_points(lat_1, long_1, lat_2, long_2, length)
    1891             : 
    1892             :     END IF
    1893             :     !
    1894       35952 :   end subroutine cellLength
    1895             : 
    1896             : 
    1897             :   ! --------------------------------------------------------------------------
    1898             : 
    1899             :   !    NAME
    1900             :   !        get_distance_two_lat_lon_points
    1901             : 
    1902             :   !    PURPOSE
    1903             :   !>       \brief estimate distance in [m] between two points in a lat-lon
    1904             : 
    1905             :   !>       \details estimate distance in [m] between two points in a lat-lon
    1906             :   !>       Code is based on one that is implemented in the VIC-3L model
    1907             : 
    1908             :   !    INTENT(IN)
    1909             :   !>       \param[in] "real(dp) :: lat1, long1, lat2, long2" latitude  of point-1
    1910             :   !>       \param[in] "real(dp) :: lat1, long1, lat2, long2" longitude of point-1
    1911             :   !>       \param[in] "real(dp) :: lat1, long1, lat2, long2" latitude  of point-2
    1912             :   !>       \param[in] "real(dp) :: lat1, long1, lat2, long2" longitude of point-2
    1913             : 
    1914             :   !    INTENT(OUT)
    1915             :   !>       \param[out] "real(dp) :: distance_out" distance between two points [m]
    1916             : 
    1917             :   !    HISTORY
    1918             :   !>       \authors Rohini Kumar
    1919             : 
    1920             :   !>       \date May 2014
    1921             : 
    1922             :   ! Modifications:
    1923             :   ! Stephan Thober Aug 2015 - ported to mRM
    1924             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
    1925             : 
    1926           0 :   subroutine get_distance_two_lat_lon_points(lat1, long1, lat2, long2, distance_out)
    1927             : 
    1928       35952 :     use mo_constants, only : RadiusEarth_dp, TWOPI_dp
    1929             : 
    1930             :     implicit none
    1931             : 
    1932             :     ! longitude of point-2
    1933             :     real(dp), intent(in) :: lat1, long1, lat2, long2
    1934             : 
    1935             :     ! distance between two points [m]
    1936             :     real(dp), intent(out) :: distance_out
    1937             : 
    1938           0 :     real(dp) :: theta1
    1939             : 
    1940           0 :     real(dp) :: phi1
    1941             : 
    1942           0 :     real(dp) :: theta2
    1943             : 
    1944           0 :     real(dp) :: phi2
    1945             : 
    1946             :     real(dp) :: dtor
    1947             : 
    1948           0 :     real(dp) :: term1
    1949             : 
    1950           0 :     real(dp) :: term2
    1951             : 
    1952           0 :     real(dp) :: term3
    1953             : 
    1954           0 :     real(dp) :: temp
    1955             : 
    1956             : 
    1957           0 :     dtor = TWOPI_dp / 360.0_dp
    1958           0 :     theta1 = dtor * long1
    1959           0 :     phi1 = dtor * lat1
    1960           0 :     theta2 = dtor * long2
    1961           0 :     phi2 = dtor * lat2
    1962             : 
    1963           0 :     term1 = cos(phi1) * cos(theta1) * cos(phi2) * cos(theta2)
    1964           0 :     term2 = cos(phi1) * sin(theta1) * cos(phi2) * sin(theta2)
    1965           0 :     term3 = sin(phi1) * sin(phi2)
    1966           0 :     temp = term1 + term2 + term3
    1967           0 :     if(temp .GT. 1.0_dp) temp = 1.0_dp
    1968             : 
    1969           0 :     distance_out = RadiusEarth_dp * acos(temp);
    1970             : 
    1971           0 :   end subroutine get_distance_two_lat_lon_points
    1972             : 
    1973             :   ! --------------------------------------------------------------------------
    1974             : 
    1975             :   !     NAME
    1976             :   !         L11_flow_accumulation
    1977             : 
    1978             :   !     PURPOSE
    1979             : 
    1980             :   !>       \brief Calculates L11 flow accumulation per grid cell
    1981             :   !>       \details Calculates L11 flow accumulation per grid cell using L11_fDir
    1982             :   !>                and L11_cellarea. L11_flow_accumulation contains the recursiv subroutine
    1983             :   !>                calculate_L11_flow_accumulation
    1984             : 
    1985             :   !     INTENT(IN)
    1986             :   !>        iDomain
    1987             : 
    1988             :   !     INTENT(INOUT)
    1989             :   !         None
    1990             : 
    1991             :   !     INTENT(OUT)
    1992             :   !         L11_fAcc, L11_LinkIn_fAcc
    1993             : 
    1994             :   !     INTENT(IN), OPTIONAL
    1995             :   !         None
    1996             : 
    1997             :   !     INTENT(INOUT), OPTIONAL
    1998             :   !         None
    1999             : 
    2000             :   !     INTENT(OUT), OPTIONAL
    2001             :   !         None
    2002             : 
    2003             :   !     RETURN
    2004             :   !         None
    2005             : 
    2006             :   !     RESTRICTIONS
    2007             :   !         None
    2008             : 
    2009             :   !     EXAMPLE
    2010             :   !         None
    2011             : 
    2012             :   !     LITERATURE
    2013             :   !         None
    2014             : 
    2015             :   !     HISTORY
    2016             :   !>        \author Matthias Kelbling
    2017             :   !>        \date   Aug 2017
    2018             :   !     Modified
    2019             :   !         Stephan Thober, Jun 2018 - refactored to fit MPR_extract
    2020             : 
    2021             :   ! --------------------------------------------------------------------------
    2022          24 :   subroutine L11_flow_accumulation(iDomain)
    2023             : 
    2024             :     use mo_mrm_global_variables, only: &
    2025           0 :         L11_fDir,          & !  IN: flow direction at L11 (standard notation)
    2026             :         L11_fAcc             ! OUT: flow accumulation at L11 [km^2]
    2027             :     use mo_mrm_global_variables, only : level11
    2028             :     use mo_common_constants, only : nodata_i4, nodata_dp
    2029             :     use mo_append, only : append
    2030             : 
    2031             :     implicit none
    2032             : 
    2033             :     ! local variables
    2034             :     integer(i4), intent(in)                   :: iDomain
    2035          24 :     real(dp),    dimension(:,:), allocatable  :: fAcc11             ! L11_fAcc array
    2036             :     integer(i4)                               :: ii, jj             ! row and col index
    2037             :     integer(i4)                               :: s11, e11           ! Vec. position iDomain - fAcc
    2038             :     integer(i4)                               :: nrows11, ncols11   ! array size Domain
    2039          24 :     integer(i4), dimension(:,:), allocatable  :: fDir11             ! L11_fDir array
    2040          24 :     logical,     dimension(:,:), allocatable  :: mask11             ! Domain mask
    2041             : 
    2042             :     ! initialize Domain info
    2043          24 :     nrows11 = level11(iDomain)%nrows
    2044          24 :     ncols11 = level11(iDomain)%ncols
    2045          24 :     s11 = level11(iDomain)%iStart
    2046          24 :     e11 = level11(iDomain)%iEnd
    2047        1860 :     mask11 = level11(iDomain)%mask
    2048             : 
    2049             :     ! allocate arrays
    2050          96 :     allocate(fDir11      (nrows11, ncols11))
    2051          96 :     allocate(fAcc11      (nrows11, ncols11))
    2052             : 
    2053             :     ! initialize
    2054        1812 :     fDir11(:,:)     = nodata_i4
    2055             : 
    2056             :     ! get data
    2057          24 :     fDir11(:,:)  = UNPACK( L11_fDir(s11:e11),  mask11, nodata_i4 )
    2058             : 
    2059             :     ! initialize fAcc11 with cell area
    2060         960 :     fAcc11 = UNPACK( level11(iDomain)%cellarea * 1.e-6,  mask11, nodata_dp )
    2061             : 
    2062             :     ! For each sink call "calculate_L11_flow_accumulation"
    2063         250 :     do jj=1, ncols11
    2064        1812 :       do ii=1, nrows11
    2065        1788 :         if (fDir11(ii,jj) .eq. 0) then
    2066             :           call calculate_L11_flow_accumulation(fDir = fDir11, &
    2067             :               fAcc = fAcc11, &
    2068             :               ii = ii, &
    2069             :               jj = jj, &
    2070             :               nrow = nrows11, &
    2071          24 :               ncol = ncols11)
    2072             :         end if
    2073             :       end do
    2074             :     end do
    2075             : 
    2076             :     ! Append
    2077          24 :     call append( L11_fAcc, pack(fAcc11(:,:),mask11))
    2078             : 
    2079             :     ! free space
    2080          24 :     deallocate(fDir11, fAcc11, mask11)
    2081             : 
    2082             :   contains
    2083             : 
    2084         936 :     recursive subroutine calculate_L11_flow_accumulation(fDir, fAcc, ii, jj, nrow, ncol)
    2085             : 
    2086             :       implicit none
    2087             : 
    2088             :       integer(i4), intent(in)            :: fDir(:,:)      ! flow Direction
    2089             :       real(dp), intent (inout)           :: fAcc(:,:)      ! flow accumulation
    2090             :       integer(i4), intent(in)            :: ii, jj         ! row and col index
    2091             :       integer(i4), intent(in)            :: nrow, ncol     ! number of rows,cols in array
    2092             : 
    2093             :       ! Scan order:
    2094             :       !
    2095             :       !    6 7 8
    2096             :       !    5 x 1
    2097             :       !    4 3 2
    2098             :       !
    2099             :       ! Each:
    2100             :       ! 1. Check if cell is inflow cell to current grid
    2101             :       ! 2. If yes: Call calculate_subroutine and add result
    2102             : 
    2103         936 :       if (jj+1 .le. ncol) then
    2104         912 :         if (fDir(ii,jj+1) .eq. 16_i4) then
    2105         402 :           call calculate_L11_flow_accumulation(fDir, fAcc, ii, jj+1, nrow, ncol)
    2106         402 :           fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii,jj+1)
    2107             :         end if
    2108             :       end if
    2109             : 
    2110         936 :       if ((ii+1 .le. nrow) .and. (jj+1 .le. ncol)) then
    2111         798 :         if (fDir(ii+1,jj+1) .eq. 32_i4) then
    2112           0 :           call calculate_L11_flow_accumulation(fDir, fAcc, ii+1, jj+1, nrow, ncol)
    2113           0 :           fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii+1,jj+1)
    2114             :         end if
    2115             :       end if
    2116             : 
    2117         936 :       if (ii+1 .le. nrow) then
    2118         822 :         if (fDir(ii+1,jj) .eq. 64_i4) then
    2119         320 :           call calculate_L11_flow_accumulation(fDir, fAcc, ii+1, jj, nrow, ncol)
    2120         320 :           fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii+1,jj)
    2121             :         end if
    2122             :       end if
    2123             : 
    2124         936 :       if ((ii+1 .le. nrow) .and. (jj-1 .ge. 1)) then
    2125         770 :         if (fDir(ii+1,jj-1) .eq. 128_i4) then
    2126           0 :           call calculate_L11_flow_accumulation(fDir, fAcc, ii+1, jj-1, nrow, ncol)
    2127           0 :           fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii+1,jj-1)
    2128             :         end if
    2129             :       end if
    2130             : 
    2131         936 :       if (jj-1 .ge. 1) then
    2132         884 :         if (fDir(ii,jj-1) .eq. 1_i4) then
    2133          50 :           call calculate_L11_flow_accumulation(fDir, fAcc, ii, jj-1, nrow, ncol)
    2134          50 :           fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii,jj-1)
    2135             :         end if
    2136             :       end if
    2137             : 
    2138         936 :       if ((ii-1 .ge. 1) .and. (jj-1 .ge. 1)) then
    2139         842 :         if (fDir(ii-1,jj-1) .eq. 2_i4) then
    2140           0 :           call calculate_L11_flow_accumulation(fDir, fAcc, ii-1, jj-1, nrow, ncol)
    2141           0 :           fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii-1,jj-1)
    2142             :         end if
    2143             :       end if
    2144             : 
    2145         936 :       if (ii-1 .ge. 1) then
    2146         892 :         if (fDir(ii-1,jj) .eq. 4_i4) then
    2147         138 :           call calculate_L11_flow_accumulation(fDir, fAcc, ii-1, jj, nrow, ncol)
    2148         138 :           fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii-1,jj)
    2149             :         end if
    2150             :       end if
    2151             : 
    2152         936 :       if ((ii-1 .ge. 1) .and. (jj+1 .le. ncol)) then
    2153         868 :         if (fDir11(ii-1,jj+1) .eq. 8_i4) then
    2154           2 :           call calculate_L11_flow_accumulation(fDir, fAcc, ii-1, jj+1, nrow, ncol)
    2155           2 :           fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii-1,jj+1)
    2156             :         end if
    2157             :       end if
    2158             : 
    2159             :       ! print *, ii, jj, nrow, ncol, fAcc(ii, jj)
    2160             : 
    2161          24 :     end subroutine calculate_L11_flow_accumulation
    2162             : 
    2163             :   end subroutine L11_flow_accumulation
    2164             : 
    2165             : 
    2166             :   ! ------------------------------------------------------------------
    2167             : 
    2168             :   !     NAME
    2169             :   !         L11_calc_celerity
    2170             : 
    2171             :   !     PURPOSE
    2172             :   !>        \brief L11 celerity based on L0 elevation and L0 fAcc
    2173             : 
    2174             :   !>        \details L11 celerity based on L0 elevation and L0 fAcc
    2175             : 
    2176             :   !     INTENT(IN)
    2177             :   !>        \param[in]
    2178             : 
    2179             :   !     INTENT(INOUT)
    2180             :   !         None
    2181             : 
    2182             :   !     INTENT(OUT)
    2183             :   !         None
    2184             : 
    2185             :   !     INTENT(IN), OPTIONAL
    2186             :   !         None
    2187             : 
    2188             :   !     INTENT(INOUT), OPTIONAL
    2189             :   !         None
    2190             : 
    2191             :   !     INTENT(OUT), OPTIONAL
    2192             :   !         None
    2193             : 
    2194             :   !     RETURN
    2195             :   !         None
    2196             : 
    2197             :   !     RESTRICTIONS
    2198             :   !         None
    2199             : 
    2200             :   !     EXAMPLE
    2201             :   !         None
    2202             : 
    2203             :   !     LITERATURE
    2204             :   !         None
    2205             : 
    2206             :   !     HISTORY
    2207             :   !>        \author  Matthias Kelbling
    2208             :   !>        \date    Oct 2017
    2209             : 
    2210             :   ! ------------------------------------------------------------------
    2211             : 
    2212          17 :   subroutine L11_calc_celerity(iDomain, param)
    2213             :     use mo_common_constants, only: nodata_i4, nodata_dp
    2214             :     use mo_mad, only: mad
    2215             :     use mo_append, only: append
    2216             :     use mo_mpr_global_variables, only: &
    2217             :         L0_slope               ! IN:    slope [%]
    2218             :     use mo_common_types, only: Grid
    2219             :     use mo_common_variables, only: &
    2220             :         domainMeta,          & ! IN:    for L0 Domain indexer
    2221             :         level0                 ! IN:    level 0 grid
    2222             :     use mo_mrm_global_variables, only: &
    2223             :         L0_fDir,             & ! IN:    flow direction (standard notation) L0
    2224             :         L0_fAcc,             & ! IN:    flow accumulation (number of cells)?
    2225             :         L0_streamNet,        & ! IN:    stream Network at Level 0
    2226             :         level11,             & ! IN:    level 11 grid
    2227             :         L11_fRow,            & ! IN:    from row in L0 grid
    2228             :         L11_fCol,            & ! IN:    from col in L0 grid
    2229             :         L11_tRow,            & ! IN:    to row in L0 grid
    2230             :         L11_tCol,            & ! IN:    to col in L0 grid
    2231             :         L11_netPerm,         & ! IN:    routing order (permutation)
    2232             :         L11_nOutlets,        & ! IN:    Number of Outlets/Sinks
    2233             :         L11_celerity,        & ! INOUT: averaged celerity
    2234             :         L0_celerity            ! INOUT
    2235             : 
    2236             :     implicit none
    2237             : 
    2238             :     integer(i4), intent(in)                  :: iDomain         ! Domain
    2239             :     real(dp), dimension(:), intent(in)       :: param
    2240             : 
    2241             :     ! local
    2242             :     integer(i4)                              :: nCells0
    2243             :     integer(i4)                              :: nNodes
    2244             :     integer(i4)                              :: nLinks
    2245             :     integer(i4)                              :: nrows0, ncols0
    2246             :     integer(i4)                              :: iStart0, iEnd0
    2247             :     integer(i4)                              :: nrows11, ncols11
    2248             :     integer(i4)                              :: iStart11, iEnd11
    2249          17 :     logical,     dimension(:,:), allocatable :: mask0
    2250          17 :     integer(i4), dimension(:,:), allocatable :: iD0
    2251          17 :     integer(i4), dimension(:,:), allocatable :: fDir0
    2252          17 :     integer(i4), dimension(:,:), allocatable :: fAcc0
    2253          17 :     real(dp),    dimension(:,:), allocatable :: slope0
    2254          17 :     real(dp),    dimension(:), allocatable :: slope_tmp
    2255          17 :     real(dp),    dimension(:,:), allocatable :: cellarea0
    2256          17 :     integer(i4), dimension(:),   allocatable :: netPerm         ! routing order (permutation)
    2257          17 :     integer(i4), dimension(:),   allocatable :: nLinkFromRow
    2258          17 :     integer(i4), dimension(:),   allocatable :: nLinkFromCol
    2259          17 :     integer(i4), dimension(:),   allocatable :: nLinkToRow
    2260          17 :     integer(i4), dimension(:),   allocatable :: nLinkToCol
    2261             :     integer(i4)                              :: ii, rr, ns
    2262             :     integer(i4)                              :: frow, fcol
    2263             :     integer(i4)                              :: fId,  tId
    2264          17 :     real(dp),    dimension(:),   allocatable :: stack, append_chunk ! Stacks celerity along the L0 river-path
    2265          17 :     integer(i4), dimension(:),   allocatable :: dummy_1d
    2266             : 
    2267          17 :     real(dp)                                 :: L0_link_slope
    2268          17 :     real(dp),    dimension(:),   allocatable :: celerity11
    2269          17 :     real(dp),    dimension(:,:), allocatable :: celerity0
    2270             : 
    2271          17 :     integer(i4), dimension(:,:), allocatable :: nodata_i4_tmp
    2272          17 :     real(dp),    dimension(:,:), allocatable :: nodata_dp_tmp
    2273          17 :     logical,     dimension(:),   allocatable :: slopemask0
    2274             : 
    2275             :     type(Grid), pointer :: level0_iDomain
    2276             : 
    2277             :     ! level-0 information
    2278          17 :     level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
    2279          17 :     nrows0 = level0_iDomain%nrows
    2280          17 :     ncols0 = level0_iDomain%ncols
    2281          17 :     nCells0 = level0_iDomain%ncells
    2282          17 :     iStart0 = level0_iDomain%iStart
    2283          17 :     iEnd0 = level0_iDomain%iEnd
    2284     2122467 :     mask0 = level0_iDomain%mask
    2285             : 
    2286             :     ! level-11 information
    2287          17 :     iStart11 = level11(iDomain)%iStart
    2288          17 :     iEnd11 = level11(iDomain)%iEnd
    2289          17 :     nrows11 = level11(iDomain)%nrows
    2290          17 :     ncols11 = level11(iDomain)%ncols
    2291          17 :     nNodes = level11(iDomain)%ncells
    2292             : 
    2293          17 :     nLinks  = nNodes - L11_nOutlets(iDomain)
    2294             : 
    2295             :     ! allocate
    2296          68 :     allocate ( iD0         ( nrows0, ncols0 ) )
    2297          68 :     allocate ( slope0      ( nrows0, ncols0 ) )
    2298          51 :     allocate ( fDir0       ( nrows0, ncols0 ) )
    2299          51 :     allocate ( fAcc0       ( nrows0, ncols0 ) )
    2300          51 :     allocate ( cellarea0   ( nrows0, ncols0 ) )
    2301          51 :     allocate ( celerity0   ( nrows0, ncols0 ) )
    2302          51 :     allocate ( slopemask0  ( ncells0 ) )
    2303             : 
    2304             :     !  Routing network vectors have nNodes size instead of nLinks to
    2305             :     !  avoid the need of having two extra indices to identify a Domain.
    2306          17 :     allocate ( stack       ( 1 ) )
    2307          17 :     allocate ( append_chunk( 1 ) )
    2308          17 :     allocate ( dummy_1d    ( 2 ))
    2309          51 :     allocate ( netPerm     ( nNodes ) )
    2310          34 :     allocate ( nLinkFromRow( nNodes ) )
    2311          34 :     allocate ( nLinkFromCol( nNodes ) )
    2312          34 :     allocate ( nLinkToRow  ( nNodes ) )
    2313          34 :     allocate ( nLinkToCol  ( nNodes ) )
    2314          51 :     allocate ( celerity11  ( nNodes ) )
    2315          34 :     allocate ( slope_tmp   ( nNodes ) )
    2316             : 
    2317          51 :     allocate (nodata_i4_tmp( nrows0, ncols0 ) )
    2318          51 :     allocate (nodata_dp_tmp( nrows0, ncols0 ) )
    2319             : 
    2320             :     ! initialize
    2321     2122433 :     iD0(:,:)           = nodata_i4
    2322     2122433 :     fDir0(:,:)         = nodata_i4
    2323     2122433 :     fAcc0(:,:)         = nodata_i4
    2324     2122433 :     cellarea0(:,:)     = nodata_dp
    2325     2122433 :     slope0(:,:)        = nodata_dp
    2326             : 
    2327          34 :     stack(:)           = nodata_dp
    2328          34 :     append_chunk(:)    = nodata_dp
    2329         595 :     netPerm(:)         = nodata_i4
    2330         595 :     nLinkFromRow(:)    = nodata_i4
    2331         595 :     nLinkFromCol(:)    = nodata_i4
    2332         595 :     nLinkToRow(:)      = nodata_i4
    2333         595 :     nLinkToCol(:)      = nodata_i4
    2334         595 :     celerity11(:)      = nodata_dp
    2335     2122433 :     celerity0(:,:)     = nodata_dp
    2336      791282 :     slopemask0(:)      = .False.
    2337             : 
    2338     2122433 :     nodata_i4_tmp(:,:) = nodata_i4
    2339     2122433 :     nodata_dp_tmp(:,:) = nodata_dp
    2340             : 
    2341             :     ! for a single node model run
    2342             : 
    2343          17 :     if(nNodes .GT. 1) then
    2344             :       ! get L0 fields
    2345          17 :       iD0(:,:) = UNPACK(level0_iDomain%Id(1:nCells0), mask0, nodata_i4_tmp)
    2346          17 :       fDir0(:,:) = UNPACK(L0_fDir(iStart0:iEnd0), mask0, nodata_i4_tmp)
    2347          17 :       fAcc0(:,:) = UNPACK(L0_fAcc(iStart0:iEnd0), mask0, nodata_i4_tmp)
    2348          17 :       cellarea0(:,:) = UNPACK(level0_iDomain%cellarea(1:nCells0), mask0, nodata_dp_tmp)
    2349             : 
    2350             :       ! smoothing river slope
    2351      791299 :       slope_tmp = L0_slope(iStart0:iEnd0)
    2352      791282 :       where ( slope_tmp .lt. 0.1_dp ) slope_tmp = 0.1_dp
    2353             : 
    2354      791282 :       slopemask0(:) = (L0_streamNet(iStart0:iEnd0) .ne. nodata_i4)
    2355             : 
    2356             :       ! smooth river cells if there is more than one cell
    2357      791282 :       if( count(slopemask0) .GT. 1) then
    2358     1582564 :         slope_tmp = mad(arr = slope_tmp, z = 2.25_dp, mask = slopemask0, tout='u', mval=0.1_dp)
    2359             :       end if
    2360          17 :       slope0(:,:) = UNPACK(slope_tmp,  mask0, nodata_dp_tmp )
    2361             : 
    2362             :       ! get network vectors of L11
    2363         595 :       netPerm(:)      = L11_netPerm ( iStart11 : iEnd11 )
    2364         595 :       nLinkFromRow(:) = L11_fRow    ( iStart11 : iEnd11 )
    2365         595 :       nLinkFromCol(:) = L11_fCol    ( iStart11 : iEnd11 )
    2366         595 :       nLinkToRow(:)   = L11_tRow    ( iStart11 : iEnd11 )
    2367         595 :       nLinkToCol(:)   = L11_tCol    ( iStart11 : iEnd11 )
    2368             : 
    2369         578 :       do rr = 1, nLinks
    2370             : 
    2371         561 :         ii   = netPerm(rr)
    2372         561 :         frow = nLinkFromRow(ii)
    2373         561 :         fcol = nLinkFromCol(ii)
    2374             : 
    2375             :         ! Init
    2376        1122 :         stack(:) = 0_dp
    2377         561 :         ns = 1
    2378             : 
    2379         561 :         fId = iD0( frow, fcol )
    2380         561 :         tId = iD0( nLinkToRow(ii) , nLinkToCol(ii) )
    2381       24514 :         do
    2382       25075 :           L0_link_slope = slope0(frow, fcol) / 100._dp
    2383             : 
    2384             :           ! celerity parametrization
    2385       25075 :           stack(ns) = param(1) * sqrt(L0_link_slope)
    2386             : 
    2387       25075 :           celerity0(frow, fcol) = stack(ns)
    2388       25075 :           ns = ns + 1
    2389       25075 :           fId = iD0(frow, fcol)
    2390       25075 :           if( .NOT. (fID == tID)) then
    2391       24514 :             call append(stack, append_chunk)
    2392             :           else
    2393             :             exit
    2394             :           end if
    2395             :           ! move downstream
    2396       24514 :           call moveDownOneCell( fDir0(frow,fcol), frow, fcol )
    2397             :         end do
    2398             : 
    2399       25636 :         celerity11(ii) = size(stack) / sum(1/stack(:))
    2400         561 :         deallocate(stack)
    2401         578 :         allocate(stack(1))
    2402             : 
    2403             :       end do
    2404             : 
    2405             :     else
    2406             : 
    2407             :       ! There is only one cell, so no routing is taking place
    2408             :       ! set dummy value of 1 m / s
    2409           0 :       celerity11(:) =  1._dp
    2410             : 
    2411             :     end if
    2412             : 
    2413             :     ! Write celerity
    2414         595 :     L11_celerity(iStart11:iEnd11) = celerity11(:)
    2415          17 :     L0_celerity(iStart0:iEnd0) = PACK(celerity0(:,:), mask0)
    2416             : 
    2417             :     ! free space
    2418           0 :     deallocate (&
    2419           0 :         mask0, iD0, slope_tmp, slopemask0, &
    2420           0 :         slope0, fDir0, cellarea0,   &
    2421          17 :         stack, netPerm, nLinkFromRow, nLinkFromCol, nLinkToRow, nLinkToCol)
    2422             : 
    2423          17 :   end subroutine L11_calc_celerity
    2424             : 
    2425             : end module mo_mrm_net_startup

Generated by: LCOV version 1.16