mo_mrm_restart.f90

Go to the documentation of this file.

!> \file mo_mrm_restart.f90
!> \brief \copybrief mo_mrm_restart
!> \details \copydetails mo_mrm_restart
 
!> \brief Restart routines
!> \details This module contains the subroutines for reading and writing routing related variables to file.
!> \authors Stephan Thober
!> \date Aug 2015
!> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
!! mHM is released under the LGPLv3+ license \license_note
!> \ingroup f_mrm
module mo_mrm_restart
  use mo_kind, only : i4, dp
  use mo_message, only : message, error_message
 
  implicit none
 
  public :: mrm_write_restart
  public :: mrm_read_restart_states
  public :: mrm_read_restart_config
contains
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        mrm_write_restart
 
  !    PURPOSE
  !>       \brief write routing states and configuration
 
  !>       \details write configuration and state variables to a given restart
  !>       directory.
 
  !    INTENT(IN)
  !>       \param[in] "integer(i4) :: iDomain"                   number of domain
  !>       \param[in] "character(256), dimension(:) :: OutFile" list of Output paths per Domain
 
  !    HISTORY
  !>       \authors Stephan Thober
 
  !>       \date Aug 2015
 
  ! Modifications:
  ! Stephan Thober    Sep 2015 - added all write restart commands in this subroutine
  ! Stephan Thober    Sep 2015 - added L11_areaCell L1_ID and L1_L11_Id for routing
  !                              resolution higher than hydrology resolution
  ! David Schaefer    Nov 2015 - mo_netcdf
  ! Stephan Thober    May 2016 - split L0_OutletCoord into L0_rowOutlet & L0_colOutlet
  !                              because multiple outlets could exist
  ! Stephan Thober    Nov 2016 - added L11_TSrout, ProcessMatrix
  ! Matthias Kelbling Aug 2017 - added L11_fAcc, L0_slope, L0_celerity, L11_celerity, L11_meandering
  ! Robert Schweppe   Jun 2018 - refactoring and reformatting
  ! Stephan Thober    Jun 2018 - including varying celerity functionality
  ! Stephan Thober    May 2019 - added L0 info required for Process 3
 
  subroutine mrm_write_restart(iDomain, domainID, OutFile)
 
    use mo_common_constants, only : nodata_dp, nodata_i4
    use mo_common_restart, only : write_grid_info
    use mo_common_variables, only : level0, level1, nlcoverscene, processmatrix, domainmeta, &
            lc_year_start, lc_year_end
    use mo_common_constants, only : landcoverperiodsvarname
    use mo_mrm_constants, only : nroutingstates
    use mo_mpr_global_variables, only : l0_slope
    use mo_mrm_global_variables, only : l0_fdir, l0_facc, l0_streamnet, &
                                        l1_l11_id, &
                                        l11_c1, l11_c2, l11_k, l11_l1_id, l11_qmod, &
                                        l11_tsrout, l11_afloodplain, l11_colout, l11_colout, l11_fcol, l11_fdir, &
                                        l11_facc, l11_frow, l11_fromn, l11_label, l11_length, l11_nlinkfracfpimp, &
                                        l11_netperm, l11_qout, l11_qtin, l11_qtr, l11_rorder, l11_rowout, l11_rowout, &
                                        l11_sink, l11_slope, l11_tcol, l11_trow, l11_ton, l11_xi, l11_celerity, &
                                        level11, domain_mrm
    use mo_netcdf, only : ncdataset, ncdimension, ncvariable
    use mo_string_utils, only : num2str
 
    implicit none
 
    ! number of domain
    integer(i4), intent(in) :: idomain
 
    ! domain
    integer(i4), intent(in) :: domainid
 
    ! list of Output paths per Domain
    character(256), dimension(:), intent(in) :: outfile
 
    character(256) :: fname
 
    integer(i4) :: ii
 
    ! number of outlets at Level 0
    integer(i4) :: noutlet
 
    ! start index at level 0
    integer(i4) :: s0
 
    ! end index at level 0
    integer(i4) :: e0
 
    ! mask at level 0
    logical, dimension(:, :), allocatable :: mask0
 
    ! start index at level 1
    integer(i4) :: s1
 
    ! end index at level 1
    integer(i4) :: e1
 
    ! mask at level 1
    logical, dimension(:, :), allocatable :: mask1
 
    ! start index at level 11
    integer(i4) :: s11
 
    ! end index at level 11
    integer(i4) :: e11
 
    ! number of colums at level 11
    integer(i4) :: ncols11
 
    ! number of rows at level 11
    integer(i4) :: nrows11
 
    ! dummy variable for writing L11_sink
    integer(i4), allocatable :: dummy(:)
 
    ! mask at level 11
    logical, dimension(:, :), allocatable :: mask11
 
    ! dummy variable
    real(dp), dimension(:, :, :), allocatable :: dummy_d3
    real(dp), dimension(:), allocatable :: dummy_d1
 
    type(ncdataset) :: nc
 
    type(ncdimension) :: rows0, cols0, rows1, cols1, rows11, cols11, it11, lcscenes, nout
 
    type(ncdimension) :: links, nts, nproc
 
    type(ncvariable) :: var
 
 
    ! get Level1 and Level11 information about the Domain
    noutlet = domain_mrm(domainmeta%L0DataFrom(idomain))%L0_noutlet
    s0 = level0(domainmeta%L0DataFrom(idomain))%iStart
    e0 = level0(domainmeta%L0DataFrom(idomain))%iEnd
    mask0 = level0(domainmeta%L0DataFrom(idomain))%mask
    s1 = level1(idomain)%iStart
    e1 = level1(idomain)%iEnd
    mask1 = level1(idomain)%mask
    s11 = level11(idomain)%iStart
    e11 = level11(idomain)%iEnd
    mask11 = level11(idomain)%mask
    ncols11 = level11(idomain)%ncols
    nrows11 = level11(idomain)%nrows
    allocate(dummy_d3(nrows11, ncols11, nroutingstates))
 
    ! set restart file name
    fname = trim(outfile(idomain))
 
    call message('    Writing mRM restart file to ' // trim(fname) // ' ...')
 
    nc = ncdataset(fname, "w")
 
    call write_grid_info(level0(domainmeta%L0DataFrom(idomain)), "0", nc)
    call write_grid_info(level1(idomain), "1", nc)
    call write_grid_info(level11(idomain), "11", nc)
 
    nout = nc%setDimension("Noutlet", noutlet)
    rows0 = nc%getDimension("nrows0")
    cols0 = nc%getDimension("ncols0")
    rows1 = nc%getDimension("nrows1")
    cols1 = nc%getDimension("ncols1")
    rows11 = nc%getDimension("nrows11")
    cols11 = nc%getDimension("ncols11")
 
    it11 = nc%setDimension("nIT", nroutingstates)
    links = nc%setDimension("nLinks", size(l11_length(s11 : e11)))
    nts = nc%setDimension("TS", 1)
    nproc = nc%setDimension("Nprocesses", size(processmatrix, dim = 1))
    allocate(dummy_d1(nlcoverscene+1))
    dummy_d1(1:nlcoverscene) = lc_year_start(:)
    ! this is done because bounds are always stored as real so e.g.
    ! 1981-1990,1991-2000 is thus saved as 1981.0-1991.0,1991.0-2001.0
    ! it is translated back into ints correctly during reading
    dummy_d1(nlcoverscene+1) = lc_year_end(nlcoverscene) + 1
    lcscenes = nc%setCoordinate(trim(landcoverperiodsvarname), nlcoverscene, dummy_d1, 0_i4)
    deallocate(dummy_d1)
 
 
    ! add processMatrix
    var = nc%setVariable("ProcessMatrix", "i32", (/nproc/))
    call var%setFillValue(nodata_i4)
    call var%setData(processmatrix(:, 1))
    call var%setAttribute("long_name", "Process Matrix")
 
    ! add L0 variables if processmatrix is equal to 3
    if (processmatrix(8, 1) .eq. 3_i4) then
      ! add L0_fdir, L0_fAcc, L0_slope, L0_streamnet
      var = nc%setVariable("L0_fDir", "i32", (/rows0, cols0/))
      call var%setFillValue(nodata_i4)
      call var%setData(unpack(l0_fdir(s0:e0), mask0, nodata_i4))
      call var%setAttribute("long_name", "flow direction at level 0")
 
      var = nc%setVariable("L0_fAcc", "i32", (/rows0, cols0/))
      call var%setFillValue(nodata_i4)
      call var%setData(unpack(l0_facc(s0:e0), mask0, nodata_i4))
      call var%setAttribute("long_name", "flow accumulation at level 0")
 
      var = nc%setVariable("L0_slope", "f64", (/rows0, cols0/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l0_slope(s0:e0), mask0, nodata_dp))
      call var%setAttribute("long_name", "slope at level 0")
 
      var = nc%setVariable("L0_streamnet", "i32", (/rows0, cols0/))
      call var%setFillValue(nodata_i4)
      call var%setData(unpack(l0_streamnet(s0:e0), mask0, nodata_i4))
      call var%setAttribute("long_name", "streamnet at level 0")
    end if
 
    var = nc%setVariable("L1_Id", "i32", (/rows1, cols1/))
    call var%setFillValue(nodata_i4)
    call var%setData(unpack(level1(idomain)%Id(1:e1-s1+1), mask1, nodata_i4))
    call var%setAttribute("long_name", "cell IDs at level 1")
 
    var = nc%setVariable("L1_L11_Id", "i32", (/rows1, cols1/))
    call var%setFillValue(nodata_i4)
    call var%setData(unpack(l1_l11_id(s1 : e1), mask1, nodata_i4))
    call var%setAttribute("long_name", "Mapping of L1 Id on L11")
 
    var = nc%setVariable("L11_Qmod", "f64", (/rows11, cols11/))
    call var%setFillValue(nodata_dp)
    call var%setData(unpack(l11_qmod(s11 : e11), mask11, nodata_dp))
    call var%setAttribute("long_name", "simulated discharge at each node at level 11")
 
    var = nc%setVariable("L11_qOUT", "f64", (/rows11, cols11/))
    call var%setFillValue(nodata_dp)
    call var%setData(unpack(l11_qout(s11 : e11), mask11, nodata_dp))
    call var%setAttribute("long_name", "Total outflow from cells L11 at time tt at level 11")
 
    do ii = 1, size(dummy_d3, 3)
      dummy_d3(:, :, ii) = unpack(l11_qtin(s11 : e11, ii), mask11, nodata_dp)
    end do
    var = nc%setVariable("L11_qTIN", "f64", (/rows11, cols11, it11/))
    call var%setFillValue(nodata_dp)
    call var%setData(dummy_d3)
    call var%setAttribute("long_name", "Total discharge inputs at t-1 and t at level 11")
 
    do ii = 1, size(dummy_d3, 3)
      dummy_d3(:, :, ii) = unpack(l11_qtr(s11 : e11, ii), mask11, nodata_dp)
    end do
    var = nc%setVariable("L11_qTR", "f64", (/rows11, cols11, it11/))
    call var%setFillValue(nodata_dp)
    call var%setData(dummy_d3)
    call var%setAttribute("long_name", "Routed outflow leaving a node at level 11")
 
    var = nc%setVariable("L11_K", "f64", (/rows11, cols11/))
    call var%setFillValue(nodata_dp)
    call var%setData(unpack(l11_k(s11 : e11), mask11, nodata_dp))
    call var%setAttribute("long_name", "kappa: Muskingum travel time parameter at level 11")
 
    var = nc%setVariable("L11_xi", "f64", (/rows11, cols11/))
    call var%setFillValue(nodata_dp)
    call var%setData(unpack(l11_xi(s11 : e11), mask11, nodata_dp))
    call var%setAttribute("long_name", "xi: Muskingum diffusion parameter at level 11")
 
    var = nc%setVariable("L11_C1", "f64", (/rows11, cols11/))
    call var%setFillValue(nodata_dp)
    call var%setData(unpack(l11_c1(s11 : e11), mask11, nodata_dp))
    call var%setAttribute("long_name", "Routing parameter C1=f(K,xi, DT) (Chow, 25-41) at level 11")
 
    var = nc%setVariable("L11_C2", "f64", (/rows11, cols11/))
    call var%setFillValue(nodata_dp)
    call var%setData(unpack(l11_c2(s11 : e11), mask11, nodata_dp))
    call var%setAttribute("long_name", "Routing parameter C2=f(K,xi, DT) (Chow, 25-41) at level 11")
 
    deallocate(dummy_d3)
    allocate(dummy_d3(nrows11, ncols11, nlcoverscene))
    do ii = 1, size(dummy_d3, 3)
      dummy_d3(:, :, ii) = unpack(l11_nlinkfracfpimp(s11 : e11, ii), mask11, nodata_dp)
    end do
    var = nc%setVariable("L11_nLinkFracFPimp", "f64", (/rows11, cols11, lcscenes/))
    call var%setFillValue(nodata_dp)
    call var%setData(dummy_d3)
    call var%setAttribute("long_name", "Fraction of the flood plain with impervious cover at level 11")
 
    ! ----------------------------------------------------------
    ! L11 config set
    ! ----------------------------------------------------------
    var = nc%setVariable("L11_domain_Mask", "i32", (/rows11, cols11/))
    call var%setFillValue(nodata_i4)
    call var%setData(merge(1_i4, 0_i4,  mask11))
    call var%setAttribute("long_name", "Mask at Level 11")
 
    var = nc%setVariable("L11_TSrout", "i32", (/nts/))
    call var%setFillValue(nodata_i4)
    call var%setData(l11_tsrout(idomain))
    call var%setAttribute("long_name", "routing resolution at Level 11")
    call var%setAttribute("units", "s")
 
    var = nc%setVariable("L11_Id", "i32", (/rows11, cols11/))
    call var%setFillValue(nodata_i4)
    call var%setData(unpack(level11(idomain)%Id(1:e11-s11+1), mask11, nodata_i4))
    call var%setAttribute("long_name", "cell Ids at Level 11")
 
    var = nc%setVariable("L11_fDir", "i32", (/rows11, cols11/))
    call var%setFillValue(nodata_i4)
    call var%setData(unpack(l11_fdir(s11 : e11), mask11, nodata_i4))
    call var%setAttribute("long_name", "flow Direction at Level 11")
 
    var = nc%setVariable("L11_fAcc", "f64", (/rows11, cols11/))
    call var%setFillValue(nodata_dp)
    call var%setData(unpack(l11_facc(s11:e11), mask11, nodata_dp))
    call var%setAttribute("long_name", "flow accumulation at Level 11")
 
    var = nc%setVariable("L11_rowOut", "i32", (/rows11, cols11/))
    call var%setFillValue(nodata_i4)
    call var%setData(unpack(l11_rowout(s11 : e11), mask11, nodata_i4))
    call var%setAttribute("long_name", "Grid vertical location of the Outlet at Level 11")
 
    var = nc%setVariable("L11_colOut", "i32", (/rows11, cols11/))
    call var%setFillValue(nodata_i4)
    call var%setData(unpack(l11_colout(s11 : e11), mask11, nodata_i4))
    call var%setAttribute("long_name", "Grid horizontal location of the Outlet at Level 11")
 
    var = nc%setVariable("L11_fromN", "i32", (/links/))
    call var%setFillValue(nodata_i4)
    call var%setData(l11_fromn(s11 : e11))
    call var%setAttribute("long_name", "From Node")
 
    var = nc%setVariable("L11_toN", "i32", (/links/))
    call var%setFillValue(nodata_i4)
    call var%setData(l11_ton(s11 : e11))
    call var%setAttribute("long_name", "To Node")
 
    var = nc%setVariable("L11_rOrder", "i32", (/links/))
    call var%setFillValue(nodata_i4)
    call var%setData(l11_rorder(s11 : e11))
    call var%setAttribute("long_name", "Network routing order at Level 11")
 
    var = nc%setVariable("L11_label", "i32", (/links/))
    call var%setFillValue(nodata_i4)
    call var%setData(l11_label(s11 : e11))
    call var%setAttribute("long_name", "Label Id [0='', 1=HeadWater, 2=Sink] at Level 11")
 
    var = nc%setVariable("L11_sink", "i32", (/links/))
    call var%setFillValue(nodata_i4)
    allocate(dummy(e11 - s11 + 1))
    dummy = 0_i4
    where(l11_sink(s11 : e11))
      dummy = 1_i4
    end where
    ! call var%setData(merge(1_i4, 0_i4, L11_sink(s11 : e11)))
    call var%setData(dummy)
    deallocate(dummy)
    call var%setAttribute("long_name", ".true. if sink node reached at Level 11")
 
    var = nc%setVariable("L11_netPerm", "i32", (/links/))
    call var%setFillValue(nodata_i4)
    call var%setData(l11_netperm(s11 : e11))
    call var%setAttribute("long_name", "Routing sequence (permutation of L11_rOrder) at Level 11")
 
    var = nc%setVariable("L11_fRow", "i32", (/links/))
    call var%setFillValue(nodata_i4)
    call var%setData(l11_frow(s11 : e11))
    call var%setAttribute("long_name", "From row in L0 grid at Level 11")
 
    var = nc%setVariable("L11_fCol", "i32", (/links/))
    call var%setFillValue(nodata_i4)
    call var%setData(l11_fcol(s11 : e11))
    call var%setAttribute("long_name", "From col in L0 grid at Level 11")
 
    var = nc%setVariable("L11_tRow", "i32", (/links/))
    call var%setFillValue(nodata_i4)
    call var%setData(l11_trow(s11 : e11))
    call var%setAttribute("long_name", "To row in L0 grid at Level 11")
 
    var = nc%setVariable("L11_tCol", "i32", (/links/))
    call var%setFillValue(nodata_i4)
    call var%setData(l11_tcol(s11 : e11))
    call var%setAttribute("long_name", "To Col in L0 grid at Level 11")
 
    var = nc%setVariable("L11_length", "f64", (/links/))
    call var%setFillValue(nodata_dp)
    call var%setData(l11_length(s11 : e11))
    call var%setAttribute("long_name", "Total length of river link [m]")
 
    var = nc%setVariable("L11_aFloodPlain", "f64", (/links/))
    call var%setFillValue(nodata_dp)
    call var%setData(l11_afloodplain(s11 : e11))
    call var%setAttribute("long_name", "Area of the flood plain [m2]")
 
    var = nc%setVariable("L11_slope", "f64", (/links/))
    call var%setFillValue(nodata_dp)
    call var%setData(l11_slope(s11 : e11))
    call var%setAttribute("long_name", "Average slope of river link")
 
    var = nc%setVariable("L11_L1_Id", "i32", (/rows11, cols11/))
    call var%setFillValue(nodata_i4)
    call var%setData(unpack(l11_l1_id(s11 : e11), mask11, nodata_i4))
    call var%setAttribute("long_name", "Mapping of L1 Id on L11")
 
    var = nc%setVariable("gaugeNodeList", "i32", &
            (/nc%setDimension("Ngauges", size(domain_mrm(idomain)%gaugeNodeList(:)))/) &
            )
    call var%setFillValue(nodata_i4)
    call var%setData(domain_mrm(idomain)%gaugeNodeList(:))
    call var%setAttribute("long_name", "cell ID of gauges")
 
    var = nc%setVariable("InflowGaugeNodeList", "i32", &
            (/nc%setDimension("nInflowGauges", size(domain_mrm(idomain)%InflowGaugeNodeList(:)))/) &
            )
    call var%setFillValue(nodata_i4)
    call var%setData(domain_mrm(idomain)%InflowGaugeNodeList(:))
    call var%setAttribute("long_name", "cell ID of gauges")
 
    if (processmatrix(8, 1) .eq. 3) then
      var = nc%setVariable("L11_celerity", "f64", (/links/))   ! celerity
      call var%setFillValue(nodata_dp)
      call var%setData(l11_celerity(s11:e11))
      call var%setAttribute("long_name", "celerity at Level 11")
    end if
 
    ! free dummy variables
    deallocate(dummy_d3)
    call nc%close()
 
  end subroutine mrm_write_restart
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        mrm_read_restart_states
 
  !    PURPOSE
  !>       \brief read routing states
 
  !>       \details This subroutine reads the routing states from
  !>       mRM_states_<domain_id>.nc that has to be in the given
  !>       path directory. This subroutine has to be called directly
  !>       each time the mHM_eval or mRM_eval is called such that the
  !>       the states are always the same at the first simulation time
  !>       step, crucial for optimization.
 
  !    INTENT(IN)
  !>       \param[in] "integer(i4) :: iDomain"    number of domains
  !>       \param[in] "character(256) :: InFile" Input Path including trailing slash
 
  !    HISTORY
  !>       \authors Stephan Thober
 
  !>       \date Sep 2015
 
  ! Modifications:
  ! David Schaefer Mar 2016 - mo_netcdf
  ! Stephan Thober May 2016 - split L0_OutletCoord into L0_rowOutlet & L0_colOutlet because multiple outlets could exist
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine mrm_read_restart_states(iDomain, domainID, InFile)
 
    use mo_common_variables, only : nlcoverscene
    use mo_mrm_constants, only : nroutingstates
    use mo_mrm_global_variables, only : l11_c1, l11_c2, l11_k, l11_qmod, l11_qout, l11_nlinkfracfpimp, l11_qtin, l11_qtr, &
                                        l11_xi, level11
    use mo_netcdf, only : ncdataset, ncvariable
    use mo_string_utils, only : num2str
 
    implicit none
 
    ! number of Domain
    integer(i4), intent(in) :: idomain
 
    integer(i4), intent(in) :: domainid
 
    ! Input Path including trailing slash
    character(256), intent(in) :: infile
 
    integer(i4) :: ii
 
    ! start index at level 11
    integer(i4) :: s11
 
    ! end index at level 11
    integer(i4) :: e11
 
    ! mask at level 11
    logical, dimension(:, :), allocatable :: mask11
 
    ! dummy, 2 dimension
    real(dp), dimension(:, :), allocatable :: dummyd2
 
    ! dummy, 3 dimension
    real(dp), dimension(:, :, :), allocatable :: dummyd3
 
    character(256) :: fname
 
    type(ncdataset) :: nc
 
    type(ncvariable) :: var
 
    !TODO-RIV-TEMP: read/write restart for riv-temp process
 
    ! set file name
    fname = trim(infile)
 
    ! get Domain info at L11 including ncells, start, end, and mask
    s11 = level11(idomain)%iStart
    e11 = level11(idomain)%iEnd
    mask11 = level11(idomain)%mask
 
    nc = ncdataset(fname, "r")
 
    ! simulated discharge at each node
    var = nc%getVariable("L11_Qmod")
    call var%getData(dummyd2)
    l11_qmod(s11 : e11) = pack(dummyd2, mask11)
 
    ! Total outflow from cells L11 at time tt
    var = nc%getVariable("L11_qOUT")
    call var%getData(dummyd2)
    l11_qout(s11 : e11) = pack(dummyd2, mask11)
 
    ! Total discharge inputs at t-1 and t
    var = nc%getVariable("L11_qTIN")
    call var%getData(dummyd3)
    do ii = 1, nroutingstates
      l11_qtin(s11 : e11, ii) = pack(dummyd3(:, :, ii), mask11)
    end do
 
    !  Routed outflow leaving a node
    var = nc%getVariable("L11_qTR")
    call var%getData(dummyd3)
    do ii = 1, nroutingstates
      l11_qtr(s11 : e11, ii) = pack(dummyd3(:, :, ii), mask11)
    end do
 
    ! kappa: Muskingum travel time parameter.
    var = nc%getVariable("L11_K")
    call var%getData(dummyd2)
    l11_k(s11 : e11) = pack(dummyd2, mask11)
 
    !  xi:    Muskingum diffusion parameter
    var = nc%getVariable("L11_xi")
    call var%getData(dummyd2)
    l11_xi(s11 : e11) = pack(dummyd2, mask11)
 
    ! Routing parameter C1=f(K,xi, DT) (Chow, 25-41)
    var = nc%getVariable("L11_C1")
    call var%getData(dummyd2)
    l11_c1(s11 : e11) = pack(dummyd2, mask11)
 
    ! Routing parameter C2 =f(K,xi, DT) (Chow, 25-41)
    var = nc%getVariable("L11_C2")
    call var%getData(dummyd2)
    l11_c2(s11 : e11) = pack(dummyd2, mask11)
 
    ! Fraction of the flood plain with impervious cover
    var = nc%getVariable("L11_nLinkFracFPimp")
    deallocate(dummyd3)
    call var%getData(dummyd3)
    do ii = 1, nlcoverscene
      l11_nlinkfracfpimp(s11 : e11, ii) = pack(dummyd3(:, :, ii), mask11)
    end do
 
    ! free memory
    deallocate(dummyd2, dummyd3)
 
  end subroutine mrm_read_restart_states
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        mrm_read_restart_config
 
  !    PURPOSE
  !>       \brief reads Level 11 configuration from a restart directory
 
  !>       \details read Level 11 configuration variables from a given restart
  !>       directory and initializes all Level 11 configuration variables,
  !>       that are initialized in L11_variable_init,
  !>       contained in module mo_startup.
 
  !    INTENT(IN)
  !>       \param[in] "integer(i4) :: iDomain"    number of Domain
  !>       \param[in] "character(256) :: InFile" Input Path including trailing slash
 
  !    HISTORY
  !>       \authors Stephan Thober
 
  !>       \date Apr 2013
 
  ! Modifications:
  ! Matthias Zink  Apr 2014 - added inflow gauge
  ! Stephan Thober Aug 2015 - adapted for mRM
  ! Stephan Thober Sep 2015 - added all read restart commands in this subroutine
  ! Stephan Thober Sep 2015 - added L11_areaCell, L1_ID and L1_L11_Id for routing resolution higher than hydrology resolution
  ! David Schaefer Mar 2016 - mo_netcdf
  ! Stephan Thober Nov 2016 - added L11_TSrout, ProcessMatrix
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
  ! Stephan Thober May 2019 - added L0 info required for Process 3
 
  subroutine mrm_read_restart_config(iDomain, domainID, InFile)
 
    use mo_append, only : append
    use mo_common_constants, only : nodata_dp
    use mo_common_variables, only : level0, level1, domainmeta, processmatrix, domainmeta
    use mo_kind, only : dp, i4
    use mo_mpr_global_variables, only : l0_slope
    use mo_mrm_global_variables, only : l0_fdir, l0_facc, l0_streamnet, &
                                        l11_l1_id, l11_tsrout, l11_afloodplain, l11_colout, l11_fcol, &
                                        l11_fdir, l11_facc, l11_frow, l11_fromn, l11_label, l11_length, l11_noutlets, l11_netperm, &
                                        l11_rorder, l11_rowout, l11_sink, l11_slope, l11_tcol, l11_trow, l11_ton, &
                                        l1_l11_id, domain_mrm, level11
    use mo_netcdf, only : ncdataset, ncvariable
    use mo_string_utils, only : num2str
 
    implicit none
 
    ! number of Domain
    integer(i4), intent(in) :: idomain
 
    ! domain
    integer(i4), intent(in) :: domainid
 
    ! Input Path including trailing slash
    character(256), intent(in) :: infile
 
    character(256) :: fname
 
    ! Mask at Level 0
    logical, allocatable, dimension(:, :) :: mask0
 
    ! Mask at Level 1
    logical, allocatable, dimension(:, :) :: mask1
 
    ! Mask at Level 11
    logical, allocatable, dimension(:, :) :: mask11
 
    ! dummy, 1 dimension I4
    integer(i4), allocatable, dimension(:) :: dummyi1
 
    ! dummy, 2 dimension I4
    integer(i4), allocatable, dimension(:, :) :: dummyi2
 
    ! dummy, DP
    real(dp), allocatable, dimension(:) :: dummyd1
    real(dp), allocatable, dimension(:, :) :: dummyd2
    real(dp), allocatable :: dummyd0
 
    type(ncdataset) :: nc
 
    type(ncvariable) :: var
 
 
    ! set file name
    fname = trim(infile)
    call message('        Reading mRM restart file:  ', trim(adjustl(fname)), ' ...')
 
    ! get Domain info at L11 mask
    mask0 = level0(domainmeta%L0DataFrom(idomain))%mask
    mask1 = level1(idomain)%mask
    mask11 = level11(idomain)%mask
 
    nc = ncdataset(fname, "r")
 
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! Read Process Matrix for check <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    var = nc%getVariable("ProcessMatrix")
    allocate(dummyi1(size(processmatrix, dim = 1)))
    call var%getData(dummyi1)
    if (dummyi1(8) .ne. processmatrix(8, 1)) then
      call error_message('***ERROR: process description for routing', raise=.false.)
      call error_message('***ERROR: given in restart file does not match', raise=.false.)
      call error_message('***ERROR: that in namelist', raise=.false.)
      call error_message('***ERROR: restart file value:. ' // num2str(dummyi1(8), '(i2)'), raise=.false.)
      call error_message('***ERROR: namelist value:..... ' // num2str(processmatrix(8, 1), '(i2)'), raise=.false.)
      call error_message('ERROR: mrm_read_restart_config')
    end if
    deallocate(dummyi1)
 
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! Read L0 variables <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    if (processmatrix(8, 1) .eq. 3_i4) then
      ! add L0_fdir, L0_fAcc, L0_slope, L0_streamnet
      var = nc%getVariable("L0_fDir")
      call var%getData(dummyi2)
      call append(l0_fdir, pack(dummyi2, mask0))
 
      var = nc%getVariable("L0_fAcc")
      call var%getData(dummyi2)
      call append(l0_facc, pack(dummyi2, mask0))
 
      var = nc%getVariable("L0_slope")
      call var%getData(dummyd2)
      call append(l0_slope, pack(dummyd2, mask0))
 
      var = nc%getVariable("L0_streamnet")
      call var%getData(dummyi2)
      call append(l0_streamnet, pack(dummyi2, mask0))
    end if
 
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! Read L1 variables <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    var = nc%getVariable("L1_L11_Id")
    call var%getData(dummyi2)
    call append(l1_l11_id, pack(dummyi2, mask1))
 
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! Read L11 variables <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! read L11_TSrout
    if (idomain .eq. 1) then
      allocate(l11_tsrout(domainmeta%nDomains))
      l11_tsrout = nodata_dp
    end if
    var = nc%getVariable("L11_TSrout")
    call var%getData(dummyd0)
    l11_tsrout(idomain) = dummyd0
 
    ! L11 data sets
    ! Mapping of L1 Ids on L1
    var = nc%getVariable("L11_L1_Id")
    call var%getData(dummyi2)
    call append(l11_l1_id, pack(dummyi2, mask11))
 
    ! Flow direction (standard notation)
    var = nc%getVariable("L11_fDir")
    call var%getData(dummyi2)
    call append(l11_fdir, pack(dummyi2, mask11))
    ! append Number of Outlets at Level 11 (where facc == 0 )
    call append(l11_noutlets, count((dummyi2 .eq. 0_i4)))
 
    ! Flow accumulation
    var = nc%getVariable("L11_fAcc")
    call var%getData(dummyd2)
    call append(l11_facc, pack(dummyd2, mask11))
 
    ! Grid vertical location of the Outlet
    var = nc%getVariable("L11_rowOut")
    call var%getData(dummyi2)
    call append(l11_rowout, pack(dummyi2, mask11))
 
    ! Grid horizontal location  of the Outlet
    var = nc%getVariable("L11_colOut")
    call var%getData(dummyi2)
    call append(l11_colout, pack(dummyi2, mask11))
 
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! From Node
    var = nc%getVariable("L11_fromN")
    call var%getData(dummyi1)
    call append(l11_fromn, dummyi1)
 
    ! To Node
    var = nc%getVariable("L11_toN")
    call var%getData(dummyi1)
    call append(l11_ton, dummyi1)
 
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! Network routing order
    var = nc%getVariable("L11_rOrder")
    call var%getData(dummyi1)
    call append(l11_rorder, dummyi1)
 
    ! Label Id [0='', 1=HeadWater, 2=Sink]
    var = nc%getVariable("L11_label")
    call var%getData(dummyi1)
    call append(l11_label, dummyi1)
 
    ! .true. if sink node reached
    var = nc%getVariable("L11_sink")
    call var%getData(dummyi1)
    call append(l11_sink, (dummyi1 .eq. 1_i4))
 
    ! Routing sequence (permutation of L11_rOrder)
    var = nc%getVariable("L11_netPerm")
    call var%getData(dummyi1)
    call append(l11_netperm, dummyi1)
 
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! From row in L0 grid
    var = nc%getVariable("L11_fRow")
    call var%getData(dummyi1)
    call append(l11_frow, dummyi1)
 
    ! From col in L0 grid
    var = nc%getVariable("L11_fCol")
    call var%getData(dummyi1)
    call append(l11_fcol, dummyi1)
 
    ! To row in L0 grid
    var = nc%getVariable("L11_tRow")
    call var%getData(dummyi1)
    call append(l11_trow, dummyi1)
 
    ! To col in L0 grid
    var = nc%getVariable("L11_tCol")
    call var%getData(dummyi1)
    call append(l11_tcol, dummyi1)
 
 
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! read gaugenodelist
    var = nc%getVariable("gaugeNodeList")
    call var%getData(dummyi1)
    domain_mrm(idomain)%gaugeNodeList(:) = dummyi1
 
    ! read InflowGaugeNodelist
    if (domain_mrm(idomain)%nInflowGauges > 0) then
      var = nc%getVariable("InflowGaugeNodeList")
      call var%getData(dummyi1)
      domain_mrm(idomain)%InflowgaugeNodeList(:) = dummyi1
    end if
 
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! L11 network data sets
    ! [m]     Total length of river link
    var = nc%getVariable("L11_length")
    call var%getData(dummyd1)
    call append(l11_length, dummyd1)
 
    ! [m2]    Area of the flood plain
    var = nc%getVariable("L11_aFloodPlain")
    call var%getData(dummyd1)
    call append(l11_afloodplain, dummyd1)
 
    ! Average slope of river link
    var = nc%getVariable("L11_slope")
    call var%getData(dummyd1)
    call append(l11_slope, dummyd1)
 
  end subroutine mrm_read_restart_config
 
end module mo_mrm_restart