LCOV - code coverage report
Current view: top level - mRM - mo_mrm_restart.f90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 366 377 97.1 %
Date: 2024-04-30 08:53:32 Functions: 3 3 100.0 %

          Line data    Source code
       1             : !> \file mo_mrm_restart.f90
       2             : !> \brief \copybrief mo_mrm_restart
       3             : !> \details \copydetails mo_mrm_restart
       4             : 
       5             : !> \brief Restart routines
       6             : !> \details This module contains the subroutines for reading and writing routing related variables to file.
       7             : !> \authors Stephan Thober
       8             : !> \date Aug 2015
       9             : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
      10             : !! mHM is released under the LGPLv3+ license \license_note
      11             : !> \ingroup f_mrm
      12             : module mo_mrm_restart
      13             :   use mo_kind, only : i4, dp
      14             :   use mo_message, only : message, error_message
      15             : 
      16             :   implicit none
      17             : 
      18             :   public :: mrm_write_restart
      19             :   public :: mrm_read_restart_states
      20             :   public :: mrm_read_restart_config
      21             : contains
      22             : 
      23             :   ! ------------------------------------------------------------------
      24             : 
      25             :   !    NAME
      26             :   !        mrm_write_restart
      27             : 
      28             :   !    PURPOSE
      29             :   !>       \brief write routing states and configuration
      30             : 
      31             :   !>       \details write configuration and state variables to a given restart
      32             :   !>       directory.
      33             : 
      34             :   !    INTENT(IN)
      35             :   !>       \param[in] "integer(i4) :: iDomain"                   number of domain
      36             :   !>       \param[in] "character(256), dimension(:) :: OutFile" list of Output paths per Domain
      37             : 
      38             :   !    HISTORY
      39             :   !>       \authors Stephan Thober
      40             : 
      41             :   !>       \date Aug 2015
      42             : 
      43             :   ! Modifications:
      44             :   ! Stephan Thober    Sep 2015 - added all write restart commands in this subroutine
      45             :   ! Stephan Thober    Sep 2015 - added L11_areaCell L1_ID and L1_L11_Id for routing
      46             :   !                              resolution higher than hydrology resolution
      47             :   ! David Schaefer    Nov 2015 - mo_netcdf
      48             :   ! Stephan Thober    May 2016 - split L0_OutletCoord into L0_rowOutlet & L0_colOutlet
      49             :   !                              because multiple outlets could exist
      50             :   ! Stephan Thober    Nov 2016 - added L11_TSrout, ProcessMatrix
      51             :   ! Matthias Kelbling Aug 2017 - added L11_fAcc, L0_slope, L0_celerity, L11_celerity, L11_meandering
      52             :   ! Robert Schweppe   Jun 2018 - refactoring and reformatting
      53             :   ! Stephan Thober    Jun 2018 - including varying celerity functionality
      54             :   ! Stephan Thober    May 2019 - added L0 info required for Process 3
      55             : 
      56          41 :   subroutine mrm_write_restart(iDomain, domainID, OutFile)
      57             : 
      58             :     use mo_common_constants, only : nodata_dp, nodata_i4
      59             :     use mo_common_restart, only : write_grid_info
      60             :     use mo_common_variables, only : level0, level1, nLCoverScene, processMatrix, domainMeta, &
      61             :             LC_year_start, LC_year_end
      62             :     use mo_common_constants, only : landCoverPeriodsVarName
      63             :     use mo_mrm_constants, only : nRoutingStates
      64             :     use mo_mpr_global_variables, only : L0_slope
      65             :     use mo_mrm_global_variables, only : L0_fdir, L0_fAcc, L0_streamnet, &
      66             :                                         L1_L11_Id, &
      67             :                                         L11_C1, L11_C2, L11_K, L11_L1_Id, L11_Qmod, &
      68             :                                         L11_TSrout, L11_aFloodPlain, L11_colOut, L11_colOut, L11_fCol, L11_fDir, &
      69             :                                         L11_fAcc, L11_fRow, L11_fromN, L11_label, L11_length, L11_nLinkFracFPimp, &
      70             :                                         L11_netPerm, L11_qOUT, L11_qTIN, L11_qTR, L11_rOrder, L11_rowOut, L11_rowOut, &
      71             :                                         L11_sink, L11_slope, L11_tCol, L11_tRow, L11_toN, L11_xi, L11_celerity, &
      72             :                                         level11, domain_mrm
      73             :     use mo_netcdf, only : NcDataset, NcDimension, NcVariable
      74             :     use mo_string_utils, only : num2str
      75             : 
      76             :     implicit none
      77             : 
      78             :     ! number of domain
      79             :     integer(i4), intent(in) :: iDomain
      80             : 
      81             :     ! domain
      82             :     integer(i4), intent(in) :: domainID
      83             : 
      84             :     ! list of Output paths per Domain
      85             :     character(256), dimension(:), intent(in) :: OutFile
      86             : 
      87             :     character(256) :: Fname
      88             : 
      89             :     integer(i4) :: ii
      90             : 
      91             :     ! number of outlets at Level 0
      92             :     integer(i4) :: Noutlet
      93             : 
      94             :     ! start index at level 0
      95             :     integer(i4) :: s0
      96             : 
      97             :     ! end index at level 0
      98             :     integer(i4) :: e0
      99             : 
     100             :     ! mask at level 0
     101          20 :     logical, dimension(:, :), allocatable :: mask0
     102             : 
     103             :     ! start index at level 1
     104             :     integer(i4) :: s1
     105             : 
     106             :     ! end index at level 1
     107             :     integer(i4) :: e1
     108             : 
     109             :     ! mask at level 1
     110          20 :     logical, dimension(:, :), allocatable :: mask1
     111             : 
     112             :     ! start index at level 11
     113             :     integer(i4) :: s11
     114             : 
     115             :     ! end index at level 11
     116             :     integer(i4) :: e11
     117             : 
     118             :     ! number of colums at level 11
     119             :     integer(i4) :: ncols11
     120             : 
     121             :     ! number of rows at level 11
     122             :     integer(i4) :: nrows11
     123             : 
     124             :     ! dummy variable for writing L11_sink
     125          20 :     integer(i4), allocatable :: dummy(:)
     126             : 
     127             :     ! mask at level 11
     128          20 :     logical, dimension(:, :), allocatable :: mask11
     129             : 
     130             :     ! dummy variable
     131          20 :     real(dp), dimension(:, :, :), allocatable :: dummy_d3
     132          20 :     real(dp), dimension(:), allocatable :: dummy_d1
     133             : 
     134             :     type(NcDataset) :: nc
     135             : 
     136             :     type(NcDimension) :: rows0, cols0, rows1, cols1, rows11, cols11, it11, lcscenes, nout
     137             : 
     138             :     type(NcDimension) :: links, nts, nproc
     139             : 
     140             :     type(NcVariable) :: var
     141             : 
     142             : 
     143             :     ! get Level1 and Level11 information about the Domain
     144          20 :     noutlet = domain_mrm(domainMeta%L0DataFrom(iDomain))%L0_noutlet
     145          20 :     s0 = level0(domainMeta%L0DataFrom(iDomain))%iStart
     146          20 :     e0 = level0(domainMeta%L0DataFrom(iDomain))%iEnd
     147     2363004 :     mask0 = level0(domainMeta%L0DataFrom(iDomain))%mask
     148          20 :     s1 = level1(iDomain)%iStart
     149          20 :     e1 = level1(iDomain)%iEnd
     150        1787 :     mask1 = level1(iDomain)%mask
     151          20 :     s11 = level11(iDomain)%iStart
     152          20 :     e11 = level11(iDomain)%iEnd
     153        1616 :     mask11 = level11(iDomain)%mask
     154          20 :     ncols11 = level11(iDomain)%ncols
     155          20 :     nrows11 = level11(iDomain)%nrows
     156         100 :     allocate(dummy_d3(nrows11, ncols11, nRoutingStates))
     157             : 
     158             :     ! set restart file name
     159          20 :     Fname = trim(OutFile(iDomain))
     160             : 
     161          20 :     call message('    Writing mRM restart file to ' // trim(Fname) // ' ...')
     162             : 
     163          20 :     nc = NcDataset(fname, "w")
     164             : 
     165          20 :     call write_grid_info(level0(domainMeta%L0DataFrom(iDomain)), "0", nc)
     166          20 :     call write_grid_info(level1(iDomain), "1", nc)
     167          20 :     call write_grid_info(level11(iDomain), "11", nc)
     168             : 
     169          20 :     nout = nc%setDimension("Noutlet", Noutlet)
     170          20 :     rows0 = nc%getDimension("nrows0")
     171          20 :     cols0 = nc%getDimension("ncols0")
     172          20 :     rows1 = nc%getDimension("nrows1")
     173          20 :     cols1 = nc%getDimension("ncols1")
     174          20 :     rows11 = nc%getDimension("nrows11")
     175          20 :     cols11 = nc%getDimension("ncols11")
     176             : 
     177          20 :     it11 = nc%setDimension("nIT", nRoutingStates)
     178          20 :     links = nc%setDimension("nLinks", size(L11_length(s11 : e11)))
     179          20 :     nts = nc%setDimension("TS", 1)
     180          20 :     nproc = nc%setDimension("Nprocesses", size(processMatrix, dim = 1))
     181          60 :     allocate(dummy_d1(nLCoverScene+1))
     182          60 :     dummy_d1(1:nLCoverScene) = LC_year_start(:)
     183             :     ! this is done because bounds are always stored as real so e.g.
     184             :     ! 1981-1990,1991-2000 is thus saved as 1981.0-1991.0,1991.0-2001.0
     185             :     ! it is translated back into ints correctly during reading
     186          20 :     dummy_d1(nLCoverScene+1) = LC_year_end(nLCoverScene) + 1
     187          20 :     lcscenes = nc%setCoordinate(trim(landCoverPeriodsVarName), nLCoverScene, dummy_d1, 0_i4)
     188          20 :     deallocate(dummy_d1)
     189             : 
     190             : 
     191             :     ! add processMatrix
     192          40 :     var = nc%setVariable("ProcessMatrix", "i32", (/nproc/))
     193          20 :     call var%setFillValue(nodata_i4)
     194          20 :     call var%setData(processMatrix(:, 1))
     195          20 :     call var%setAttribute("long_name", "Process Matrix")
     196             : 
     197             :     ! add L0 variables if processmatrix is equal to 3
     198          20 :     if (processMatrix(8, 1) .eq. 3_i4) then
     199             :       ! add L0_fdir, L0_fAcc, L0_slope, L0_streamnet
     200           9 :       var = nc%setVariable("L0_fDir", "i32", (/rows0, cols0/))
     201           3 :       call var%setFillValue(nodata_i4)
     202           3 :       call var%setData(unpack(L0_fdir(s0:e0), mask0, nodata_i4))
     203           3 :       call var%setAttribute("long_name", "flow direction at level 0")
     204             : 
     205           9 :       var = nc%setVariable("L0_fAcc", "i32", (/rows0, cols0/))
     206           3 :       call var%setFillValue(nodata_i4)
     207           3 :       call var%setData(unpack(L0_fAcc(s0:e0), mask0, nodata_i4))
     208           3 :       call var%setAttribute("long_name", "flow accumulation at level 0")
     209             : 
     210           9 :       var = nc%setVariable("L0_slope", "f64", (/rows0, cols0/))
     211           3 :       call var%setFillValue(nodata_dp)
     212           3 :       call var%setData(unpack(L0_slope(s0:e0), mask0, nodata_dp))
     213           3 :       call var%setAttribute("long_name", "slope at level 0")
     214             : 
     215           9 :       var = nc%setVariable("L0_streamnet", "i32", (/rows0, cols0/))
     216           3 :       call var%setFillValue(nodata_i4)
     217           3 :       call var%setData(unpack(L0_streamnet(s0:e0), mask0, nodata_i4))
     218           3 :       call var%setAttribute("long_name", "streamnet at level 0")
     219             :     end if
     220             : 
     221          60 :     var = nc%setVariable("L1_Id", "i32", (/rows1, cols1/))
     222          20 :     call var%setFillValue(nodata_i4)
     223          20 :     call var%setData(unpack(level1(iDomain)%Id(1:e1-s1+1), mask1, nodata_i4))
     224          20 :     call var%setAttribute("long_name", "cell IDs at level 1")
     225             : 
     226          60 :     var = nc%setVariable("L1_L11_Id", "i32", (/rows1, cols1/))
     227          20 :     call var%setFillValue(nodata_i4)
     228          20 :     call var%setData(unpack(L1_L11_Id(s1 : e1), mask1, nodata_i4))
     229          20 :     call var%setAttribute("long_name", "Mapping of L1 Id on L11")
     230             : 
     231          60 :     var = nc%setVariable("L11_Qmod", "f64", (/rows11, cols11/))
     232          20 :     call var%setFillValue(nodata_dp)
     233          20 :     call var%setData(unpack(L11_Qmod(s11 : e11), mask11, nodata_dp))
     234          20 :     call var%setAttribute("long_name", "simulated discharge at each node at level 11")
     235             : 
     236          60 :     var = nc%setVariable("L11_qOUT", "f64", (/rows11, cols11/))
     237          20 :     call var%setFillValue(nodata_dp)
     238          20 :     call var%setData(unpack(L11_qOUT(s11 : e11), mask11, nodata_dp))
     239          20 :     call var%setAttribute("long_name", "Total outflow from cells L11 at time tt at level 11")
     240             : 
     241          60 :     do ii = 1, size(dummy_d3, 3)
     242          60 :       dummy_d3(:, :, ii) = unpack(L11_qTIN(s11 : e11, ii), mask11, nodata_dp)
     243             :     end do
     244          80 :     var = nc%setVariable("L11_qTIN", "f64", (/rows11, cols11, it11/))
     245          20 :     call var%setFillValue(nodata_dp)
     246          20 :     call var%setData(dummy_d3)
     247          20 :     call var%setAttribute("long_name", "Total discharge inputs at t-1 and t at level 11")
     248             : 
     249          60 :     do ii = 1, size(dummy_d3, 3)
     250          60 :       dummy_d3(:, :, ii) = unpack(L11_qTR(s11 : e11, ii), mask11, nodata_dp)
     251             :     end do
     252          80 :     var = nc%setVariable("L11_qTR", "f64", (/rows11, cols11, it11/))
     253          20 :     call var%setFillValue(nodata_dp)
     254          20 :     call var%setData(dummy_d3)
     255          20 :     call var%setAttribute("long_name", "Routed outflow leaving a node at level 11")
     256             : 
     257          60 :     var = nc%setVariable("L11_K", "f64", (/rows11, cols11/))
     258          20 :     call var%setFillValue(nodata_dp)
     259          20 :     call var%setData(unpack(L11_K(s11 : e11), mask11, nodata_dp))
     260          20 :     call var%setAttribute("long_name", "kappa: Muskingum travel time parameter at level 11")
     261             : 
     262          60 :     var = nc%setVariable("L11_xi", "f64", (/rows11, cols11/))
     263          20 :     call var%setFillValue(nodata_dp)
     264          20 :     call var%setData(unpack(L11_xi(s11 : e11), mask11, nodata_dp))
     265          20 :     call var%setAttribute("long_name", "xi: Muskingum diffusion parameter at level 11")
     266             : 
     267          60 :     var = nc%setVariable("L11_C1", "f64", (/rows11, cols11/))
     268          20 :     call var%setFillValue(nodata_dp)
     269          20 :     call var%setData(unpack(L11_C1(s11 : e11), mask11, nodata_dp))
     270          20 :     call var%setAttribute("long_name", "Routing parameter C1=f(K,xi, DT) (Chow, 25-41) at level 11")
     271             : 
     272          60 :     var = nc%setVariable("L11_C2", "f64", (/rows11, cols11/))
     273          20 :     call var%setFillValue(nodata_dp)
     274          20 :     call var%setData(unpack(L11_C2(s11 : e11), mask11, nodata_dp))
     275          20 :     call var%setAttribute("long_name", "Routing parameter C2=f(K,xi, DT) (Chow, 25-41) at level 11")
     276             : 
     277          20 :     deallocate(dummy_d3)
     278         100 :     allocate(dummy_d3(nrows11, ncols11, nLCoverScene))
     279          80 :     do ii = 1, size(dummy_d3, 3)
     280          60 :       dummy_d3(:, :, ii) = unpack(L11_nLinkFracFPimp(s11 : e11, ii), mask11, nodata_dp)
     281             :     end do
     282          80 :     var = nc%setVariable("L11_nLinkFracFPimp", "f64", (/rows11, cols11, lcscenes/))
     283          20 :     call var%setFillValue(nodata_dp)
     284          20 :     call var%setData(dummy_d3)
     285          20 :     call var%setAttribute("long_name", "Fraction of the flood plain with impervious cover at level 11")
     286             : 
     287             :     ! ----------------------------------------------------------
     288             :     ! L11 config set
     289             :     ! ----------------------------------------------------------
     290          60 :     var = nc%setVariable("L11_domain_Mask", "i32", (/rows11, cols11/))
     291          20 :     call var%setFillValue(nodata_i4)
     292        1556 :     call var%setData(merge(1_i4, 0_i4,  mask11))
     293          20 :     call var%setAttribute("long_name", "Mask at Level 11")
     294             : 
     295          40 :     var = nc%setVariable("L11_TSrout", "i32", (/nts/))
     296          20 :     call var%setFillValue(nodata_i4)
     297          20 :     call var%setData(L11_TSrout(iDomain))
     298          20 :     call var%setAttribute("long_name", "routing resolution at Level 11")
     299          20 :     call var%setAttribute("units", "s")
     300             : 
     301          60 :     var = nc%setVariable("L11_Id", "i32", (/rows11, cols11/))
     302          20 :     call var%setFillValue(nodata_i4)
     303          20 :     call var%setData(unpack(level11(iDomain)%Id(1:e11-s11+1), mask11, nodata_i4))
     304          20 :     call var%setAttribute("long_name", "cell Ids at Level 11")
     305             : 
     306          60 :     var = nc%setVariable("L11_fDir", "i32", (/rows11, cols11/))
     307          20 :     call var%setFillValue(nodata_i4)
     308          20 :     call var%setData(unpack(L11_fDir(s11 : e11), mask11, nodata_i4))
     309          20 :     call var%setAttribute("long_name", "flow Direction at Level 11")
     310             : 
     311          60 :     var = nc%setVariable("L11_fAcc", "f64", (/rows11, cols11/))
     312          20 :     call var%setFillValue(nodata_dp)
     313          20 :     call var%setData(unpack(L11_fAcc(s11:e11), mask11, nodata_dp))
     314          20 :     call var%setAttribute("long_name", "flow accumulation at Level 11")
     315             : 
     316          60 :     var = nc%setVariable("L11_rowOut", "i32", (/rows11, cols11/))
     317          20 :     call var%setFillValue(nodata_i4)
     318          20 :     call var%setData(unpack(L11_rowOut(s11 : e11), mask11, nodata_i4))
     319          20 :     call var%setAttribute("long_name", "Grid vertical location of the Outlet at Level 11")
     320             : 
     321          60 :     var = nc%setVariable("L11_colOut", "i32", (/rows11, cols11/))
     322          20 :     call var%setFillValue(nodata_i4)
     323          20 :     call var%setData(unpack(L11_colOut(s11 : e11), mask11, nodata_i4))
     324          20 :     call var%setAttribute("long_name", "Grid horizontal location of the Outlet at Level 11")
     325             : 
     326          40 :     var = nc%setVariable("L11_fromN", "i32", (/links/))
     327          20 :     call var%setFillValue(nodata_i4)
     328          20 :     call var%setData(L11_fromN(s11 : e11))
     329          20 :     call var%setAttribute("long_name", "From Node")
     330             : 
     331          40 :     var = nc%setVariable("L11_toN", "i32", (/links/))
     332          20 :     call var%setFillValue(nodata_i4)
     333          20 :     call var%setData(L11_toN(s11 : e11))
     334          20 :     call var%setAttribute("long_name", "To Node")
     335             : 
     336          40 :     var = nc%setVariable("L11_rOrder", "i32", (/links/))
     337          20 :     call var%setFillValue(nodata_i4)
     338          20 :     call var%setData(L11_rOrder(s11 : e11))
     339          20 :     call var%setAttribute("long_name", "Network routing order at Level 11")
     340             : 
     341          40 :     var = nc%setVariable("L11_label", "i32", (/links/))
     342          20 :     call var%setFillValue(nodata_i4)
     343          20 :     call var%setData(L11_label(s11 : e11))
     344          20 :     call var%setAttribute("long_name", "Label Id [0='', 1=HeadWater, 2=Sink] at Level 11")
     345             : 
     346          40 :     var = nc%setVariable("L11_sink", "i32", (/links/))
     347          20 :     call var%setFillValue(nodata_i4)
     348          80 :     allocate(dummy(e11 - s11 + 1))
     349         840 :     dummy = 0_i4
     350         820 :     where(L11_sink(s11 : e11))
     351             :       dummy = 1_i4
     352             :     end where
     353             :     ! call var%setData(merge(1_i4, 0_i4, L11_sink(s11 : e11)))
     354          20 :     call var%setData(dummy)
     355          20 :     deallocate(dummy)
     356          20 :     call var%setAttribute("long_name", ".true. if sink node reached at Level 11")
     357             : 
     358          40 :     var = nc%setVariable("L11_netPerm", "i32", (/links/))
     359          20 :     call var%setFillValue(nodata_i4)
     360          20 :     call var%setData(L11_netPerm(s11 : e11))
     361          20 :     call var%setAttribute("long_name", "Routing sequence (permutation of L11_rOrder) at Level 11")
     362             : 
     363          40 :     var = nc%setVariable("L11_fRow", "i32", (/links/))
     364          20 :     call var%setFillValue(nodata_i4)
     365          20 :     call var%setData(L11_fRow(s11 : e11))
     366          20 :     call var%setAttribute("long_name", "From row in L0 grid at Level 11")
     367             : 
     368          40 :     var = nc%setVariable("L11_fCol", "i32", (/links/))
     369          20 :     call var%setFillValue(nodata_i4)
     370          20 :     call var%setData(L11_fCol(s11 : e11))
     371          20 :     call var%setAttribute("long_name", "From col in L0 grid at Level 11")
     372             : 
     373          40 :     var = nc%setVariable("L11_tRow", "i32", (/links/))
     374          20 :     call var%setFillValue(nodata_i4)
     375          20 :     call var%setData(L11_tRow(s11 : e11))
     376          20 :     call var%setAttribute("long_name", "To row in L0 grid at Level 11")
     377             : 
     378          40 :     var = nc%setVariable("L11_tCol", "i32", (/links/))
     379          20 :     call var%setFillValue(nodata_i4)
     380          20 :     call var%setData(L11_tCol(s11 : e11))
     381          20 :     call var%setAttribute("long_name", "To Col in L0 grid at Level 11")
     382             : 
     383          40 :     var = nc%setVariable("L11_length", "f64", (/links/))
     384          20 :     call var%setFillValue(nodata_dp)
     385          20 :     call var%setData(L11_length(s11 : e11))
     386          20 :     call var%setAttribute("long_name", "Total length of river link [m]")
     387             : 
     388          40 :     var = nc%setVariable("L11_aFloodPlain", "f64", (/links/))
     389          20 :     call var%setFillValue(nodata_dp)
     390          20 :     call var%setData(L11_aFloodPlain(s11 : e11))
     391          20 :     call var%setAttribute("long_name", "Area of the flood plain [m2]")
     392             : 
     393          40 :     var = nc%setVariable("L11_slope", "f64", (/links/))
     394          20 :     call var%setFillValue(nodata_dp)
     395          20 :     call var%setData(L11_slope(s11 : e11))
     396          20 :     call var%setAttribute("long_name", "Average slope of river link")
     397             : 
     398          60 :     var = nc%setVariable("L11_L1_Id", "i32", (/rows11, cols11/))
     399          20 :     call var%setFillValue(nodata_i4)
     400          20 :     call var%setData(unpack(L11_L1_Id(s11 : e11), mask11, nodata_i4))
     401          20 :     call var%setAttribute("long_name", "Mapping of L1 Id on L11")
     402             : 
     403             :     var = nc%setVariable("gaugeNodeList", "i32", &
     404           0 :             (/nc%setDimension("Ngauges", size(domain_mrm(iDomain)%gaugeNodeList(:)))/) &
     405          40 :             )
     406          20 :     call var%setFillValue(nodata_i4)
     407          20 :     call var%setData(domain_mrm(iDomain)%gaugeNodeList(:))
     408          20 :     call var%setAttribute("long_name", "cell ID of gauges")
     409             : 
     410             :     var = nc%setVariable("InflowGaugeNodeList", "i32", &
     411           0 :             (/nc%setDimension("nInflowGauges", size(domain_mrm(iDomain)%InflowGaugeNodeList(:)))/) &
     412          40 :             )
     413          20 :     call var%setFillValue(nodata_i4)
     414          20 :     call var%setData(domain_mrm(iDomain)%InflowGaugeNodeList(:))
     415          20 :     call var%setAttribute("long_name", "cell ID of gauges")
     416             : 
     417          20 :     if (processMatrix(8, 1) .eq. 3) then
     418           6 :       var = nc%setVariable("L11_celerity", "f64", (/links/))   ! celerity
     419           3 :       call var%setFillValue(nodata_dp)
     420           3 :       call var%setData(L11_celerity(s11:e11))
     421           3 :       call var%setAttribute("long_name", "celerity at Level 11")
     422             :     end if
     423             : 
     424             :     ! free dummy variables
     425          20 :     deallocate(dummy_d3)
     426          20 :     call nc%close()
     427             : 
     428          80 :   end subroutine mrm_write_restart
     429             : 
     430             :   ! ------------------------------------------------------------------
     431             : 
     432             :   !    NAME
     433             :   !        mrm_read_restart_states
     434             : 
     435             :   !    PURPOSE
     436             :   !>       \brief read routing states
     437             : 
     438             :   !>       \details This subroutine reads the routing states from
     439             :   !>       mRM_states_<domain_id>.nc that has to be in the given
     440             :   !>       path directory. This subroutine has to be called directly
     441             :   !>       each time the mHM_eval or mRM_eval is called such that the
     442             :   !>       the states are always the same at the first simulation time
     443             :   !>       step, crucial for optimization.
     444             : 
     445             :   !    INTENT(IN)
     446             :   !>       \param[in] "integer(i4) :: iDomain"    number of domains
     447             :   !>       \param[in] "character(256) :: InFile" Input Path including trailing slash
     448             : 
     449             :   !    HISTORY
     450             :   !>       \authors Stephan Thober
     451             : 
     452             :   !>       \date Sep 2015
     453             : 
     454             :   ! Modifications:
     455             :   ! David Schaefer Mar 2016 - mo_netcdf
     456             :   ! Stephan Thober May 2016 - split L0_OutletCoord into L0_rowOutlet & L0_colOutlet because multiple outlets could exist
     457             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     458             : 
     459           1 :   subroutine mrm_read_restart_states(iDomain, domainID, InFile)
     460             : 
     461          20 :     use mo_common_variables, only : nLCoverScene
     462             :     use mo_mrm_constants, only : nRoutingStates
     463             :     use mo_mrm_global_variables, only : L11_C1, L11_C2, L11_K, L11_Qmod, L11_Qout, L11_nLinkFracFPimp, L11_qTIN, L11_qTR, &
     464             :                                         L11_xi, level11
     465             :     use mo_netcdf, only : NcDataset, NcVariable
     466             :     use mo_string_utils, only : num2str
     467             : 
     468             :     implicit none
     469             : 
     470             :     ! number of Domain
     471             :     integer(i4), intent(in) :: iDomain
     472             : 
     473             :     integer(i4), intent(in) :: domainID
     474             : 
     475             :     ! Input Path including trailing slash
     476             :     character(256), intent(in) :: InFile
     477             : 
     478             :     integer(i4) :: ii
     479             : 
     480             :     ! start index at level 11
     481             :     integer(i4) :: s11
     482             : 
     483             :     ! end index at level 11
     484             :     integer(i4) :: e11
     485             : 
     486             :     ! mask at level 11
     487           1 :     logical, dimension(:, :), allocatable :: mask11
     488             : 
     489             :     ! dummy, 2 dimension
     490           1 :     real(dp), dimension(:, :), allocatable :: dummyD2
     491             : 
     492             :     ! dummy, 3 dimension
     493           1 :     real(dp), dimension(:, :, :), allocatable :: dummyD3
     494             : 
     495             :     character(256) :: fname
     496             : 
     497             :     type(NcDataset) :: nc
     498             : 
     499             :     type(NcVariable) :: var
     500             : 
     501             :     !TODO-RIV-TEMP: read/write restart for riv-temp process
     502             : 
     503             :     ! set file name
     504           1 :     fname = trim(InFile)
     505             : 
     506             :     ! get Domain info at L11 including ncells, start, end, and mask
     507           1 :     s11 = level11(iDomain)%iStart
     508           1 :     e11 = level11(iDomain)%iEnd
     509          67 :     mask11 = level11(iDomain)%mask
     510             : 
     511           1 :     nc = NcDataset(fname, "r")
     512             : 
     513             :     ! simulated discharge at each node
     514           1 :     var = nc%getVariable("L11_Qmod")
     515           1 :     call var%getData(dummyD2)
     516           1 :     L11_Qmod(s11 : e11) = pack(dummyD2, mask11)
     517             : 
     518             :     ! Total outflow from cells L11 at time tt
     519           1 :     var = nc%getVariable("L11_qOUT")
     520           1 :     call var%getData(dummyD2)
     521           1 :     L11_qOUT(s11 : e11) = pack(dummyD2, mask11)
     522             : 
     523             :     ! Total discharge inputs at t-1 and t
     524           1 :     var = nc%getVariable("L11_qTIN")
     525           1 :     call var%getData(dummyD3)
     526           3 :     do ii = 1, nRoutingStates
     527           3 :       L11_qTIN(s11 : e11, ii) = pack(dummyD3(:, :, ii), mask11)
     528             :     end do
     529             : 
     530             :     !  Routed outflow leaving a node
     531           1 :     var = nc%getVariable("L11_qTR")
     532           1 :     call var%getData(dummyD3)
     533           3 :     do ii = 1, nRoutingStates
     534           3 :       L11_qTR(s11 : e11, ii) = pack(dummyD3(:, :, ii), mask11)
     535             :     end do
     536             : 
     537             :     ! kappa: Muskingum travel time parameter.
     538           1 :     var = nc%getVariable("L11_K")
     539           1 :     call var%getData(dummyD2)
     540           1 :     L11_K(s11 : e11) = pack(dummyD2, mask11)
     541             : 
     542             :     !  xi:    Muskingum diffusion parameter
     543           1 :     var = nc%getVariable("L11_xi")
     544           1 :     call var%getData(dummyD2)
     545           1 :     L11_xi(s11 : e11) = pack(dummyD2, mask11)
     546             : 
     547             :     ! Routing parameter C1=f(K,xi, DT) (Chow, 25-41)
     548           1 :     var = nc%getVariable("L11_C1")
     549           1 :     call var%getData(dummyD2)
     550           1 :     L11_C1(s11 : e11) = pack(dummyD2, mask11)
     551             : 
     552             :     ! Routing parameter C2 =f(K,xi, DT) (Chow, 25-41)
     553           1 :     var = nc%getVariable("L11_C2")
     554           1 :     call var%getData(dummyD2)
     555           1 :     L11_C2(s11 : e11) = pack(dummyD2, mask11)
     556             : 
     557             :     ! Fraction of the flood plain with impervious cover
     558           1 :     var = nc%getVariable("L11_nLinkFracFPimp")
     559           1 :     deallocate(dummyD3)
     560           1 :     call var%getData(dummyD3)
     561           3 :     do ii = 1, nLCoverScene
     562           3 :       L11_nLinkFracFPimp(s11 : e11, ii) = pack(dummyD3(:, :, ii), mask11)
     563             :     end do
     564             : 
     565             :     ! free memory
     566           1 :     deallocate(dummyD2, dummyD3)
     567             : 
     568           1 :   end subroutine mrm_read_restart_states
     569             : 
     570             :   ! ------------------------------------------------------------------
     571             : 
     572             :   !    NAME
     573             :   !        mrm_read_restart_config
     574             : 
     575             :   !    PURPOSE
     576             :   !>       \brief reads Level 11 configuration from a restart directory
     577             : 
     578             :   !>       \details read Level 11 configuration variables from a given restart
     579             :   !>       directory and initializes all Level 11 configuration variables,
     580             :   !>       that are initialized in L11_variable_init,
     581             :   !>       contained in module mo_startup.
     582             : 
     583             :   !    INTENT(IN)
     584             :   !>       \param[in] "integer(i4) :: iDomain"    number of Domain
     585             :   !>       \param[in] "character(256) :: InFile" Input Path including trailing slash
     586             : 
     587             :   !    HISTORY
     588             :   !>       \authors Stephan Thober
     589             : 
     590             :   !>       \date Apr 2013
     591             : 
     592             :   ! Modifications:
     593             :   ! Matthias Zink  Apr 2014 - added inflow gauge
     594             :   ! Stephan Thober Aug 2015 - adapted for mRM
     595             :   ! Stephan Thober Sep 2015 - added all read restart commands in this subroutine
     596             :   ! Stephan Thober Sep 2015 - added L11_areaCell, L1_ID and L1_L11_Id for routing resolution higher than hydrology resolution
     597             :   ! David Schaefer Mar 2016 - mo_netcdf
     598             :   ! Stephan Thober Nov 2016 - added L11_TSrout, ProcessMatrix
     599             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     600             :   ! Stephan Thober May 2019 - added L0 info required for Process 3
     601             : 
     602           1 :   subroutine mrm_read_restart_config(iDomain, domainID, InFile)
     603             : 
     604           1 :     use mo_append, only : append
     605             :     use mo_common_constants, only : nodata_dp
     606             :     use mo_common_variables, only : level0, level1, domainMeta, processMatrix, domainMeta
     607             :     use mo_kind, only : dp, i4
     608             :     use mo_mpr_global_variables, only : L0_slope
     609             :     use mo_mrm_global_variables, only : L0_fdir, L0_fAcc, L0_streamnet, &
     610             :                                         L11_L1_Id, L11_TSrout, L11_aFloodPlain, L11_colOut, L11_fCol, &
     611             :                                         L11_fDir, L11_fAcc, L11_fRow, L11_fromN, L11_label, L11_length, L11_nOutlets, L11_netPerm, &
     612             :                                         L11_rOrder, L11_rowOut, L11_sink, L11_slope, L11_tCol, L11_tRow, L11_toN, &
     613             :                                         L1_L11_Id, domain_mrm, level11
     614             :     use mo_netcdf, only : NcDataset, NcVariable
     615             :     use mo_string_utils, only : num2str
     616             : 
     617             :     implicit none
     618             : 
     619             :     ! number of Domain
     620             :     integer(i4), intent(in) :: iDomain
     621             : 
     622             :     ! domain
     623             :     integer(i4), intent(in) :: domainID
     624             : 
     625             :     ! Input Path including trailing slash
     626             :     character(256), intent(in) :: InFile
     627             : 
     628             :     character(256) :: fname
     629             : 
     630             :     ! Mask at Level 0
     631           1 :     logical, allocatable, dimension(:, :) :: mask0
     632             : 
     633             :     ! Mask at Level 1
     634           1 :     logical, allocatable, dimension(:, :) :: mask1
     635             : 
     636             :     ! Mask at Level 11
     637           1 :     logical, allocatable, dimension(:, :) :: mask11
     638             : 
     639             :     ! dummy, 1 dimension I4
     640           1 :     integer(i4), allocatable, dimension(:) :: dummyI1
     641             : 
     642             :     ! dummy, 2 dimension I4
     643           1 :     integer(i4), allocatable, dimension(:, :) :: dummyI2
     644             : 
     645             :     ! dummy, DP
     646           1 :     real(dp), allocatable, dimension(:) :: dummyD1
     647           1 :     real(dp), allocatable, dimension(:, :) :: dummyD2
     648           1 :     real(dp), allocatable :: dummyD0
     649             : 
     650             :     type(NcDataset) :: nc
     651             : 
     652             :     type(NcVariable) :: var
     653             : 
     654             : 
     655             :     ! set file name
     656           1 :     fname = trim(InFile)
     657           1 :     call message('        Reading mRM restart file:  ', trim(adjustl(Fname)), ' ...')
     658             : 
     659             :     ! get Domain info at L11 mask
     660      124852 :     mask0 = level0(domainMeta%L0DataFrom(iDomain))%mask
     661          67 :     mask1 = level1(iDomain)%mask
     662          67 :     mask11 = level11(iDomain)%mask
     663             : 
     664           1 :     nc = NcDataset(fname, "r")
     665             : 
     666             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     667             :     ! Read Process Matrix for check <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     668             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     669           1 :     var = nc%getVariable("ProcessMatrix")
     670           1 :     allocate(dummyI1(size(processMatrix, dim = 1)))
     671           1 :     call var%getData(dummyI1)
     672           1 :     if (dummyI1(8) .ne. processMatrix(8, 1)) then
     673           0 :       call error_message('***ERROR: process description for routing', raise=.false.)
     674           0 :       call error_message('***ERROR: given in restart file does not match', raise=.false.)
     675           0 :       call error_message('***ERROR: that in namelist', raise=.false.)
     676           0 :       call error_message('***ERROR: restart file value:. ' // num2str(dummyI1(8), '(i2)'), raise=.false.)
     677           0 :       call error_message('***ERROR: namelist value:..... ' // num2str(processMatrix(8, 1), '(i2)'), raise=.false.)
     678           0 :       call error_message('ERROR: mrm_read_restart_config')
     679             :     end if
     680           1 :     deallocate(dummyI1)
     681             : 
     682             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     683             :     ! Read L0 variables <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     684             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     685           1 :     if (processMatrix(8, 1) .eq. 3_i4) then
     686             :       ! add L0_fdir, L0_fAcc, L0_slope, L0_streamnet
     687           1 :       var = nc%getVariable("L0_fDir")
     688           1 :       call var%getData(dummyI2)
     689           1 :       call append(L0_fdir, pack(dummyI2, mask0))
     690             : 
     691           1 :       var = nc%getVariable("L0_fAcc")
     692           1 :       call var%getData(dummyI2)
     693           1 :       call append(L0_fAcc, pack(dummyI2, mask0))
     694             : 
     695           1 :       var = nc%getVariable("L0_slope")
     696           1 :       call var%getData(dummyD2)
     697           1 :       call append(L0_Slope, pack(dummyD2, mask0))
     698             : 
     699           1 :       var = nc%getVariable("L0_streamnet")
     700           1 :       call var%getData(dummyI2)
     701           1 :       call append(L0_streamnet, pack(dummyI2, mask0))
     702             :     end if
     703             : 
     704             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     705             :     ! Read L1 variables <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     706             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     707           1 :     var = nc%getVariable("L1_L11_Id")
     708           1 :     call var%getData(dummyI2)
     709           1 :     call append(L1_L11_Id, pack(dummyI2, mask1))
     710             : 
     711             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     712             :     ! Read L11 variables <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     713             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     714             :     ! read L11_TSrout
     715           1 :     if (iDomain .eq. 1) then
     716           3 :       allocate(L11_TSrout(domainMeta%nDomains))
     717           3 :       L11_tsRout = nodata_dp
     718             :     end if
     719           1 :     var = nc%getVariable("L11_TSrout")
     720           1 :     call var%getData(dummyD0)
     721           1 :     L11_TSrout(iDomain) = dummyD0
     722             : 
     723             :     ! L11 data sets
     724             :     ! Mapping of L1 Ids on L1
     725           1 :     var = nc%getVariable("L11_L1_Id")
     726           1 :     call var%getData(dummyI2)
     727           1 :     call append(L11_L1_Id, pack(dummyI2, mask11))
     728             : 
     729             :     ! Flow direction (standard notation)
     730           1 :     var = nc%getVariable("L11_fDir")
     731           1 :     call var%getData(dummyI2)
     732           1 :     call append(L11_fDir, pack(dummyI2, mask11))
     733             :     ! append Number of Outlets at Level 11 (where facc == 0 )
     734          64 :     call append(L11_nOutlets, count((dummyI2 .eq. 0_i4)))
     735             : 
     736             :     ! Flow accumulation
     737           1 :     var = nc%getVariable("L11_fAcc")
     738           1 :     call var%getData(dummyD2)
     739           1 :     call append(L11_fAcc, pack(dummyD2, mask11))
     740             : 
     741             :     ! Grid vertical location of the Outlet
     742           1 :     var = nc%getVariable("L11_rowOut")
     743           1 :     call var%getData(dummyI2)
     744           1 :     call append(L11_rowOut, pack(dummyI2, mask11))
     745             : 
     746             :     ! Grid horizontal location  of the Outlet
     747           1 :     var = nc%getVariable("L11_colOut")
     748           1 :     call var%getData(dummyI2)
     749           1 :     call append(L11_colOut, pack(dummyI2, mask11))
     750             : 
     751             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     752             :     ! From Node
     753           1 :     var = nc%getVariable("L11_fromN")
     754           1 :     call var%getData(dummyI1)
     755           1 :     call append(L11_fromN, dummyI1)
     756             : 
     757             :     ! To Node
     758           1 :     var = nc%getVariable("L11_toN")
     759           1 :     call var%getData(dummyI1)
     760           1 :     call append(L11_toN, dummyI1)
     761             : 
     762             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     763             :     ! Network routing order
     764           1 :     var = nc%getVariable("L11_rOrder")
     765           1 :     call var%getData(dummyI1)
     766           1 :     call append(L11_rOrder, dummyI1)
     767             : 
     768             :     ! Label Id [0='', 1=HeadWater, 2=Sink]
     769           1 :     var = nc%getVariable("L11_label")
     770           1 :     call var%getData(dummyI1)
     771           1 :     call append(L11_label, dummyI1)
     772             : 
     773             :     ! .true. if sink node reached
     774           1 :     var = nc%getVariable("L11_sink")
     775           1 :     call var%getData(dummyI1)
     776          35 :     call append(L11_sink, (dummyI1 .eq. 1_i4))
     777             : 
     778             :     ! Routing sequence (permutation of L11_rOrder)
     779           1 :     var = nc%getVariable("L11_netPerm")
     780           1 :     call var%getData(dummyI1)
     781           1 :     call append(L11_netPerm, dummyI1)
     782             : 
     783             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     784             :     ! From row in L0 grid
     785           1 :     var = nc%getVariable("L11_fRow")
     786           1 :     call var%getData(dummyI1)
     787           1 :     call append(L11_fRow, dummyI1)
     788             : 
     789             :     ! From col in L0 grid
     790           1 :     var = nc%getVariable("L11_fCol")
     791           1 :     call var%getData(dummyI1)
     792           1 :     call append(L11_fCol, dummyI1)
     793             : 
     794             :     ! To row in L0 grid
     795           1 :     var = nc%getVariable("L11_tRow")
     796           1 :     call var%getData(dummyI1)
     797           1 :     call append(L11_tRow, dummyI1)
     798             : 
     799             :     ! To col in L0 grid
     800           1 :     var = nc%getVariable("L11_tCol")
     801           1 :     call var%getData(dummyI1)
     802           1 :     call append(L11_tCol, dummyI1)
     803             : 
     804             : 
     805             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     806             :     ! read gaugenodelist
     807           1 :     var = nc%getVariable("gaugeNodeList")
     808           1 :     call var%getData(dummyI1)
     809           2 :     domain_mrm(iDomain)%gaugeNodeList(:) = dummyI1
     810             : 
     811             :     ! read InflowGaugeNodelist
     812           1 :     if (domain_mrm(iDomain)%nInflowGauges > 0) then
     813           0 :       var = nc%getVariable("InflowGaugeNodeList")
     814           0 :       call var%getData(dummyI1)
     815           0 :       domain_mrm(iDomain)%InflowgaugeNodeList(:) = dummyI1
     816             :     end if
     817             : 
     818             :     ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
     819             :     ! L11 network data sets
     820             :     ! [m]     Total length of river link
     821           1 :     var = nc%getVariable("L11_length")
     822           1 :     call var%getData(dummyD1)
     823           1 :     call append(L11_length, dummyD1)
     824             : 
     825             :     ! [m2]    Area of the flood plain
     826           1 :     var = nc%getVariable("L11_aFloodPlain")
     827           1 :     call var%getData(dummyD1)
     828           1 :     call append(L11_aFloodPlain, dummyD1)
     829             : 
     830             :     ! Average slope of river link
     831           1 :     var = nc%getVariable("L11_slope")
     832           1 :     call var%getData(dummyD1)
     833           1 :     call append(L11_slope, dummyD1)
     834             : 
     835           1 :   end subroutine mrm_read_restart_config
     836             : 
     837             : end module mo_mrm_restart

Generated by: LCOV version 1.16