mo_restart.f90

Go to the documentation of this file.

!> \file mo_restart.f90
!> \brief   \copybrief mo_restart
!> \details \copydetails mo_restart
 
!> \brief reading and writing states, fluxes and configuration for restart of mHM.
!> \details routines are seperated for reading and writing variables for:
!!       - states and fluxes, and
!!       - configuration.
!!
!! Reading of L11 configuration is also seperated from the rest,
!! since it is only required when routing is activated.
!> \authors Stephan Thober
!> \date Jul 2013
!> \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_mhm
MODULE mo_restart
 
  use mo_common_constants, only : soilhorizonsvarname, landcoverperiodsvarname, laivarname
 
  IMPLICIT NONE
 
  PRIVATE
 
  PUBLIC :: read_restart_states     ! read restart files for state variables from a given path
  PUBLIC :: write_restart_files     ! write restart files for configuration to a given path
 
CONTAINS
 
  !> \brief write restart files for each domain
  !> \details write restart files for each domain. For each domain
  !! three restart files are written. These are xxx_states.nc,
  !! xxx_L11_config.nc, and xxx_config.nc (xxx being the three digit
  !! domain index). If a variable is added here, it should also be added
  !! in the read restart routines below.
  !> \changelog
  !! - Stephan Thober     Aug  2015
  !!   - moved write of routing states to mRM
  !! - David Schaefer     Nov  2015
  !!   - mo_netcdf
  !! - Stephan Thober     Nov  2016
  !!   - moved processMatrix to common variables
  !! - Zink M. Demirel C. Mar 2017
  !!   - Added Jarvis soil water stress function at SM process(3)
  !! - Robert Schweppe    Feb 2018
  !!   - Removed all L0 references
  !! - Robert Schweppe    Jun 2018
  !!   - refactoring and reformatting
  !! - Stephan Thober     Dec 2022
  !!   - added grid info for level0
  !> \authors Stephan Thober
  !> \date Jun 2014
  subroutine write_restart_files(OutFile)
 
    use mo_common_constants, only : nodata_dp
    use mo_common_restart, only : write_grid_info
    use mo_common_variables, only : level0, level1, nlcoverscene, domainmeta, lc_year_start, lc_year_end
    use mo_global_variables, only : l1_inter, l1_throughfall, l1_aetcanopy, l1_aetsealed, l1_aetsoil, l1_baseflow, &
                                    l1_fastrunoff, l1_infilsoil, l1_melt, l1_percol, l1_preeffect, l1_rain, &
                                    l1_runoffseal, l1_satstw, l1_sealstw, l1_slowrunoff, l1_snow, l1_snowpack, &
                                    l1_soilmoist, l1_total_runoff, l1_unsatstw
 
    use mo_kind, only : dp, i4
    use mo_message, only : message
    use mo_mpr_global_variables, only : nlai, nsoilhorizons_mhm, horizondepth_mhm, laiboundaries
    use mo_mpr_restart, only : write_eff_params
    use mo_netcdf, only : ncdataset, ncdimension, ncvariable
    use mo_string_utils, only : num2str
 
    implicit none
 
    character(256) :: fname
 
    ! Output Path for each domain
    character(256), dimension(:), intent(in) :: outfile
 
    integer(i4) :: idomain, domainid
 
    integer(i4) :: ii
 
    ! start index at level 1
    integer(i4) :: s1
 
    ! end index at level 1
    integer(i4) :: e1
 
    ! mask at level 1
    logical, dimension(:, :), allocatable :: mask1
 
    ! dummy variable
    real(dp), dimension(:, :, :), allocatable :: dummy_3d
    real(dp), dimension(:), allocatable :: dummy_1d
 
    integer(i4) :: max_extent
 
    type(ncdataset) :: nc
 
    type(ncdimension) :: rows1, cols1, soil1, lcscenes, lais
 
    type(ncvariable) :: var
 
 
    ! get maximum extent of one dimension 2 or 3
    max_extent = max(nsoilhorizons_mhm, nlcoverscene, nlai)
 
    domain_loop : do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
 
      ! write restart file for iDomain
      fname = trim(outfile(idomain))
      ! print a message
      call message("    Writing Restart-file: ", trim(adjustl(fname)), " ...")
 
      nc = ncdataset(fname, "w")
 
      call write_grid_info(level1(idomain), "1", nc)
      call write_grid_info(level0(domainmeta%L0DataFrom(idomain)), "0", nc)
 
      rows1 = nc%getDimension("nrows1")
      cols1 = nc%getDimension("ncols1")
 
      ! write the dimension to the file and also save bounds
      allocate(dummy_1d(nsoilhorizons_mhm+1))
      dummy_1d(1) = 0.0_dp
      dummy_1d(2:nsoilhorizons_mhm+1) = horizondepth_mhm(:)
      soil1 = nc%setCoordinate(trim(soilhorizonsvarname), nsoilhorizons_mhm, dummy_1d, 2_i4)
      deallocate(dummy_1d)
      allocate(dummy_1d(nlcoverscene+1))
      dummy_1d(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_1d(nlcoverscene+1) = lc_year_end(nlcoverscene) + 1
      lcscenes = nc%setCoordinate(trim(landcoverperiodsvarname), nlcoverscene, dummy_1d, 0_i4)
      deallocate(dummy_1d)
      ! write the dimension to the file
      lais = nc%setCoordinate(trim(laivarname), nlai, laiboundaries, 0_i4)
 
      ! for appending and intialization
      allocate(dummy_3d(rows1%getLength(), cols1%getLength(), max_extent))
      allocate(mask1(rows1%getLength(), cols1%getLength()))
      s1 = level1(idomain)%iStart
      e1 = level1(idomain)%iEnd
      mask1 = level1(idomain)%mask
 
      var = nc%setVariable("L1_Inter", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_inter(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "Interception storage at level 1")
 
      var = nc%setVariable("L1_snowPack", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_snowpack(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "Snowpack at level 1")
 
      var = nc%setVariable("L1_sealSTW", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_sealstw(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "Retention storage of impervious areas at level 1")
 
      do ii = 1, nsoilhorizons_mhm
        dummy_3d(:, :, ii) = unpack(l1_soilmoist(s1 : e1, ii), mask1, nodata_dp)
      end do
 
      var = nc%setVariable("L1_soilMoist", "f64", (/rows1, cols1, soil1/))
      call var%setFillValue(nodata_dp)
      call var%setData(dummy_3d(:, :, 1 : nsoilhorizons_mhm))
      call var%setAttribute("long_name", "soil moisture at level 1")
 
      var = nc%setVariable("L1_unsatSTW", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_unsatstw(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "upper soil storage at level 1")
 
      var = nc%setVariable("L1_satSTW", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_satstw(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "groundwater storage at level 1")
 
      do ii = 1, nsoilhorizons_mhm
        dummy_3d(:, :, ii) = unpack(l1_aetsoil(s1 : e1, ii), mask1, nodata_dp)
      end do
 
      var = nc%setVariable("L1_aETSoil", "f64", (/rows1, cols1, soil1/))
      call var%setFillValue(nodata_dp)
      call var%setData(dummy_3d(:, :, 1 : nsoilhorizons_mhm))
      call var%setAttribute("long_name", "soil actual ET at level 1")
 
      var = nc%setVariable("L1_aETCanopy", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_aetcanopy(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "canopy actual ET at level 1")
 
      var = nc%setVariable("L1_aETSealed", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_aetsealed(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "sealed actual ET at level 1")
 
      var = nc%setVariable("L1_baseflow", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_baseflow(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "baseflow at level 1")
 
      do ii = 1, nsoilhorizons_mhm
        dummy_3d(:, :, ii) = unpack(l1_infilsoil(s1 : e1, ii), mask1, nodata_dp)
      end do
 
      var = nc%setVariable("L1_infilSoil", "f64", (/rows1, cols1, soil1/))
      call var%setFillValue(nodata_dp)
      call var%setData(dummy_3d(:, :, 1 : nsoilhorizons_mhm))
      call var%setAttribute("long_name", "soil in-exfiltration at level 1")
 
      var = nc%setVariable("L1_fastRunoff", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_fastrunoff(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "fast runoff")
 
      var = nc%setVariable("L1_percol", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_percol(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "percolation at level 1")
 
      var = nc%setVariable("L1_melt", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_melt(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "snow melt at level 1")
 
      var = nc%setVariable("L1_preEffect", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_preeffect(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "effective precip. depth (snow melt + rain) at level 1")
 
      var = nc%setVariable("L1_rain", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_rain(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "rain (liquid water) at level 1")
 
      var = nc%setVariable("L1_runoffSeal", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_runoffseal(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "runoff from impervious area at level 1")
 
      var = nc%setVariable("L1_slowRunoff", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_slowrunoff(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "slow runoff at level 1")
 
      var = nc%setVariable("L1_snow", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_snow(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "snow (solid water) at level 1")
 
      var = nc%setVariable("L1_Throughfall", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_throughfall(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "throughfall at level 1")
 
      var = nc%setVariable("L1_total_runoff", "f64", (/rows1, cols1/))
      call var%setFillValue(nodata_dp)
      call var%setData(unpack(l1_total_runoff(s1 : e1), mask1, nodata_dp))
      call var%setAttribute("long_name", "total runoff at level 1")
 
      call write_eff_params(mask1, s1, e1, rows1, cols1, soil1, lcscenes, lais, nc)
 
      call nc%close()
 
      deallocate(dummy_3d, mask1)
    end do domain_loop
 
  end subroutine write_restart_files
 
 
  !> \brief reads fluxes and state variables from file
  !> \details read fluxes and state variables from given
  !! restart directory and initialises all state variables
  !! that are initialized in the subroutine initialise,
  !! contained in module mo_startup.
  !> \changelog
  !! - Stephan Thober Aug  2015
  !!   - moved read of routing states to mRM
  !! - David Schaefer Nov  2015
  !!   - mo_netcdf
  !! - Stephan Thober Nov  2016
  !!   - moved processMatrix to common variables
  !! - Robert Schweppe Jun 2018
  !!   - refactoring and reformatting
  !! - Sebastian Müller Mar 2023
  !!   - compatibility layer for 2D/3D data
  !!   - move reading of nLAI to mo_startup (needed beforehand)
  !> \authors Stephan Thober
  !> \date Apr 2013
  subroutine read_restart_states(iDomain, domainID, InFile)
 
    use mo_common_variables, only : lc_year_end, lc_year_start, level1, nlcoverscene, processmatrix
    use mo_global_variables, only : l1_inter, l1_throughfall, l1_aetcanopy, &
                                    l1_aetsealed, l1_aetsoil, l1_baseflow, l1_fastrunoff, l1_infilsoil, l1_melt, &
                                    l1_percol, l1_preeffect, l1_rain, l1_runoffseal, l1_satstw, l1_sealstw, &
                                    l1_slowrunoff, l1_snow, l1_snowpack, l1_soilmoist, l1_total_runoff, l1_unsatstw
    use mo_kind, only : dp, i4
    use mo_mpr_global_variables, only : l1_harsamcoeff, &
                                        l1_prietayalpha, l1_aeroresist, l1_alpha, l1_degday, l1_degdayinc, l1_degdaymax, &
                                        l1_degdaynopre, l1_fasp, l1_froots, l1_fsealed, l1_jarvis_thresh_c1, &
                                        l1_kbaseflow, l1_kperco, l1_kslowflow, l1_karstloss, l1_kfastflow, l1_maxinter, &
                                        l1_petlaicorfactor, l1_sealedthresh, l1_soilmoistexp, l1_soilmoistfc, &
                                        l1_soilmoistsat, l1_surfresist, l1_tempthresh, l1_unsatthresh, l1_wiltingpoint, &
                                        nlai, nsoilhorizons_mhm, &
                                        ! neutron count
                                        l1_no_count, l1_bulkdens, l1_latticewater, l1_cosmicl3
 
    use mo_netcdf, only : ncdataset, ncdimension, ncvariable
    use mo_string_utils, only : num2str
    use mo_message, only: message, error_message
 
    implicit none
 
    !> number of domain
    integer(i4), intent(in) :: idomain
    !> ID of domain
    integer(i4), intent(in) :: domainid
    !> Input Path including trailing slash
    character(256), intent(in) :: infile
 
    character(256) :: fname
    ! variable rank
    integer(i4) :: var_rank
    ! loop index
    integer(i4) :: ii, jj
    ! start index at level 1
    integer(i4) :: s1
    ! end index at level 1
    integer(i4) :: e1
    ! mask at level 1
    logical, dimension(:, :), allocatable :: mask1
    ! dummy, 2 dimension
    real(dp), dimension(:, :), allocatable :: dummyd2
    ! dummy, 3 dimension
    real(dp), dimension(:, :, :), allocatable :: dummyd3
    ! dummy, 4 dimension
    real(dp), dimension(:, :, :, :), allocatable :: dummyd4
 
    type(ncdataset) :: nc
    type(ncvariable) :: var
 
 
    fname = trim(infile)
    call message('    Reading states from ', trim(adjustl(fname)),' ...')
 
    ! get domain information at level 1
    allocate(mask1(level1(idomain)%nrows, level1(idomain)%ncols))
    mask1 = level1(idomain)%mask
    s1 = level1(idomain)%iStart
    e1 = level1(idomain)%iEnd
 
    nc = ncdataset(fname, "r")
 
    if (nc%hasVariable('L1_Inter')) then
      !-------------------------------------------
      ! STATE VARIABLES (optionally)
      !-------------------------------------------
 
      ! Interception
      var = nc%getVariable("L1_Inter")
      call var%getData(dummyd2)
      l1_inter(s1 : e1) = pack(dummyd2, mask1)
 
      ! Snowpack
      var = nc%getVariable("L1_snowPack")
      call var%getData(dummyd2)
      l1_snowpack(s1 : e1) = pack(dummyd2, mask1)
 
      ! Retention storage of impervious areas
      var = nc%getVariable("L1_sealSTW")
      call var%getData(dummyd2)
      l1_sealstw(s1 : e1) = pack(dummyd2, mask1)
 
      ! upper soil storage
      var = nc%getVariable("L1_unsatSTW")
      call var%getData(dummyd2)
      l1_unsatstw(s1 : e1) = pack(dummyd2, mask1)
 
      ! groundwater storage
      var = nc%getVariable("L1_satSTW")
      call var%getData(dummyd2)
      l1_satstw(s1 : e1) = pack(dummyd2, mask1)
 
      ! Soil moisture of each horizon
      var = nc%getVariable("L1_soilMoist")
      call var%getData(dummyd3)
      do ii = 1, nsoilhorizons_mhm
        l1_soilmoist(s1 : e1, ii) = pack(dummyd3(:, :, ii), mask1)
      end do
 
      !-------------------------------------------
      ! FLUXES (optionally)
      !-------------------------------------------
 
      !  soil actual ET
      var = nc%getVariable("L1_aETSoil")
      call var%getData(dummyd3)
      do ii = 1, nsoilhorizons_mhm
        l1_aetsoil(s1 : e1, ii) = pack(dummyd3(:, :, ii), mask1)
      end do
 
      ! canopy actual ET
      var = nc%getVariable("L1_aETCanopy")
      call var%getData(dummyd2)
      l1_aetcanopy(s1 : e1) = pack(dummyd2, mask1)
 
      ! sealed area actual ET
      var = nc%getVariable("L1_aETSealed")
      call var%getData(dummyd2)
      l1_aetsealed(s1 : e1) = pack(dummyd2, mask1)
 
      ! baseflow
      var = nc%getVariable("L1_baseflow")
      call var%getData(dummyd2)
      l1_baseflow(s1 : e1) = pack(dummyd2, mask1)
 
      ! soil in-exfiltration
      var = nc%getVariable("L1_infilSoil")
      call var%getData(dummyd3)
      do ii = 1, nsoilhorizons_mhm
        l1_infilsoil(s1 : e1, ii) = pack(dummyd3(:, :, ii), mask1)
      end do
 
      ! fast runoff
      var = nc%getVariable("L1_fastRunoff")
      call var%getData(dummyd2)
      l1_fastrunoff(s1 : e1) = pack(dummyd2, mask1)
 
      ! snow melt
      var = nc%getVariable("L1_melt")
      call var%getData(dummyd2)
      l1_melt(s1 : e1) = pack(dummyd2, mask1)
 
      ! percolation
      var = nc%getVariable("L1_percol")
      call var%getData(dummyd2)
      l1_percol(s1 : e1) = pack(dummyd2, mask1)
 
      ! effective precip. depth (snow melt + rain)
      var = nc%getVariable("L1_preEffect")
      call var%getData(dummyd2)
      l1_preeffect(s1 : e1) = pack(dummyd2, mask1)
 
      ! rain (liquid water)
      var = nc%getVariable("L1_rain")
      call var%getData(dummyd2)
      l1_rain(s1 : e1) = pack(dummyd2, mask1)
 
      ! runoff from impervious area
      var = nc%getVariable("L1_runoffSeal")
      call var%getData(dummyd2)
      l1_runoffseal(s1 : e1) = pack(dummyd2, mask1)
 
      ! slow runoff
      var = nc%getVariable("L1_slowRunoff")
      call var%getData(dummyd2)
      l1_slowrunoff(s1 : e1) = pack(dummyd2, mask1)
 
      ! snow (solid water)
      var = nc%getVariable("L1_snow")
      call var%getData(dummyd2)
      l1_snow(s1 : e1) = pack(dummyd2, mask1)
 
      ! throughfall
      var = nc%getVariable("L1_Throughfall")
      call var%getData(dummyd2)
      l1_throughfall(s1 : e1) = pack(dummyd2, mask1)
 
      ! total runoff
      var = nc%getVariable("L1_total_runoff")
      call var%getData(dummyd2)
      l1_total_runoff(s1 : e1) = pack(dummyd2, mask1)
    end if
 
    !-------------------------------------------
    ! EFFECTIVE PARAMETERS
    !-------------------------------------------
    var = nc%getVariable("L1_fSealed")
    call var%getData(dummyd3)
    do ii = 1, nlcoverscene
      l1_fsealed(s1 : e1, 1, ii) = pack(dummyd3(:, :, ii), mask1)
    end do
 
    ! exponent for the upper reservoir
    var = nc%getVariable("L1_alpha")
    var_rank = var%getrank()
    select case(var_rank)
      case(2)
        ! distribute over all land cover scenes
        call var%getData(dummyd2)
        do ii = 1, nlcoverscene
          l1_alpha(s1 : e1, 1, ii) = pack(dummyd2, mask1)
        end do
      case(3)
        call var%getData(dummyd3)
        do ii = 1, nlcoverscene
          l1_alpha(s1 : e1, 1, ii) = pack(dummyd3(:, :, ii), mask1)
        end do
      case default
        call error_message("Restart: L1_alpha rank needs to be 2 or 3")
    end select
 
    ! increase of the Degree-day factor per mm of increase in precipitation
    var = nc%getVariable("L1_degDayInc")
    call var%getData(dummyd3)
    do ii = 1, nlcoverscene
      l1_degdayinc(s1 : e1, 1, ii) = pack(dummyd3(:, :, ii), mask1)
    end do
 
    ! maximum degree-day factor
    var = nc%getVariable("L1_degDayMax")
    call var%getData(dummyd3)
    do ii = 1, nlcoverscene
      l1_degdaymax(s1 : e1, 1, ii) = pack(dummyd3(:, :, ii), mask1)
    end do
 
    ! degree-day factor with no precipitation
    var = nc%getVariable("L1_degDayNoPre")
    call var%getData(dummyd3)
    do ii = 1, nlcoverscene
      l1_degdaynopre(s1 : e1, 1, ii) = pack(dummyd3(:, :, ii), mask1)
    end do
 
    ! degree-day factor
    var = nc%getVariable("L1_degDay")
    var_rank = var%getrank()
    select case(var_rank)
      case(2)
        ! distribute over all land cover scenes
        call var%getData(dummyd2)
        do ii = 1, nlcoverscene
          l1_degday(s1 : e1, 1, ii) = pack(dummyd2, mask1)
        end do
      case(3)
        call var%getData(dummyd3)
        do ii = 1, nlcoverscene
          l1_degday(s1 : e1, 1, ii) = pack(dummyd3(:, :, ii), mask1)
        end do
      case default
        call error_message("Restart: L1_degDay rank needs to be 2 or 3")
    end select
 
    ! Karstic percolation loss
    var = nc%getVariable("L1_karstLoss")
    call var%getData(dummyd2)
    l1_karstloss(s1 : e1, 1, 1) = pack(dummyd2, mask1)
 
    ! Fraction of roots in soil horizons
    var = nc%getVariable("L1_fRoots")
    call var%getData(dummyd4)
    do jj = 1, nlcoverscene
      do ii = 1, nsoilhorizons_mhm
        l1_froots(s1 : e1, ii, jj) = pack(dummyd4(:, :, ii, jj), mask1)
      end do
    end do
 
    ! Maximum interception
    var = nc%getVariable("L1_maxInter")
    call var%getData(dummyd3)
    do ii = 1, nlai
      l1_maxinter(s1 : e1, ii, 1) = pack(dummyd3(:, :, ii), mask1)
    end do
 
    ! fast interflow recession coefficient ("L1_kFastFlow" or "L1_kfastFlow")
    if (nc%hasvariable("L1_kFastFlow")) then
      var = nc%getVariable("L1_kFastFlow")
    else
      var = nc%getVariable("L1_kfastFlow")
    end if
    call var%getData(dummyd3)
    do ii = 1, nlcoverscene
      l1_kfastflow(s1 : e1, 1, ii) = pack(dummyd3(:, :, ii), mask1)
    end do
 
    ! slow interflow recession coefficient
    var = nc%getVariable("L1_kSlowFlow")
    var_rank = var%getrank()
    select case(var_rank)
      case(2)
        ! distribute over all land cover scenes
        call var%getData(dummyd2)
        do ii = 1, nlcoverscene
          l1_kslowflow(s1 : e1, 1, ii) = pack(dummyd2, mask1)
        end do
      case(3)
        call var%getData(dummyd3)
        do ii = 1, nlcoverscene
          l1_kslowflow(s1 : e1, 1, ii) = pack(dummyd3(:, :, ii), mask1)
        end do
      case default
        call error_message("Restart: L1_kSlowFlow rank needs to be 2 or 3")
    end select
 
    ! baseflow recession coefficient
    var = nc%getVariable("L1_kBaseFlow")
    var_rank = var%getrank()
    select case(var_rank)
      case(2)
        ! distribute over all land cover scenes
        call var%getData(dummyd2)
        do ii = 1, nlcoverscene
          l1_kbaseflow(s1 : e1, 1, ii) = pack(dummyd2, mask1)
        end do
      case(3)
        call var%getData(dummyd3)
        do ii = 1, nlcoverscene
          l1_kbaseflow(s1 : e1, 1, ii) = pack(dummyd3(:, :, ii), mask1)
        end do
      case default
        call error_message("Restart: L1_kBaseFlow rank needs to be 2 or 3")
    end select
 
    ! percolation coefficient
    var = nc%getVariable("L1_kPerco")
    var_rank = var%getrank()
    select case(var_rank)
      case(2)
        ! distribute over all land cover scenes
        call var%getData(dummyd2)
        do ii = 1, nlcoverscene
          l1_kperco(s1 : e1, 1, ii) = pack(dummyd2, mask1)
        end do
      case(3)
        call var%getData(dummyd3)
        do ii = 1, nlcoverscene
          l1_kperco(s1 : e1, 1, ii) = pack(dummyd3(:, :, ii), mask1)
        end do
      case default
        call error_message("Restart: L1_kPerco rank needs to be 2 or 3")
    end select
 
    ! Soil moisture below which actual ET is reduced linearly till PWP
    ! for processCase(3) = 1
    var = nc%getVariable("L1_soilMoistFC")
    call var%getData(dummyd4)
    do jj = 1, nlcoverscene
      do ii = 1, nsoilhorizons_mhm
        l1_soilmoistfc(s1 : e1, ii, jj) = pack(dummyd4(:, :, ii, jj), mask1)
      end do
    end do
 
    ! Saturation soil moisture for each horizon [mm]
    var = nc%getVariable("L1_soilMoistSat")
    call var%getData(dummyd4)
    do jj = 1, nlcoverscene
      do ii = 1, nsoilhorizons_mhm
        l1_soilmoistsat(s1 : e1, ii, jj) = pack(dummyd4(:, :, ii, jj), mask1)
      end do
    end do
 
    ! Exponential parameter to how non-linear is the soil water retention
    var = nc%getVariable("L1_soilMoistExp")
    call var%getData(dummyd4)
    do jj = 1, nlcoverscene
      do ii = 1, nsoilhorizons_mhm
        l1_soilmoistexp(s1 : e1, ii, jj) = pack(dummyd4(:, :, ii, jj), mask1)
      end do
    end do
 
    if (any(processmatrix(3, 1) == (/2, 3/))) then
      ! jarvis critical value for normalized soil water content
      var = nc%getVariable("L1_jarvis_thresh_c1")
      call var%getData(dummyd2)
      l1_jarvis_thresh_c1(s1 : e1, 1, 1) = pack(dummyd2, mask1)
    end if
 
    ! Threshold temperature for snow/rain
    var = nc%getVariable("L1_tempThresh")
    call var%getData(dummyd3)
    do ii = 1, nlcoverscene
      l1_tempthresh(s1 : e1, 1, ii) = pack(dummyd3(:, :, ii), mask1)
    end do
 
    ! Threshold water depth controlling fast interflow
    var = nc%getVariable("L1_unsatThresh")
    call var%getData(dummyd2)
    l1_unsatthresh(s1 : e1, 1, 1) = pack(dummyd2, mask1)
 
    ! Threshold water depth for surface runoff in sealed surfaces
    var = nc%getVariable("L1_sealedThresh")
    call var%getData(dummyd2)
    l1_sealedthresh(s1 : e1, 1, 1) = pack(dummyd2, mask1)
 
    ! Permanent wilting point
    var = nc%getVariable("L1_wiltingPoint")
    call var%getData(dummyd4)
    do jj = 1, nlcoverscene
      do ii = 1, nsoilhorizons_mhm
        l1_wiltingpoint(s1 : e1, ii, jj) = pack(dummyd4(:, :, ii, jj), mask1)
      end do
    end do
 
    ! different parameters dependent on PET formulation
    select case (processmatrix(5, 1))
    case(-1) ! PET is input
 
      ! PET correction factor due to LAI
      var = nc%getVariable("L1_petLAIcorFactor")
      call var%getData(dummyd4)
      do jj = 1, nlcoverscene
        do ii = 1, nlai
          l1_petlaicorfactor(s1 : e1, ii, jj) = pack(dummyd4(:, :, ii, jj), mask1)
        end do
      end do
 
    case(0) ! PET is input
 
      ! PET correction factor due to terrain aspect
      var = nc%getVariable("L1_fAsp")
      call var%getData(dummyd2)
      l1_fasp(s1 : e1, 1, 1) = pack(dummyd2, mask1)
 
    case(1) ! Hargreaves-Samani
 
      ! PET correction factor due to terrain aspect
      var = nc%getVariable("L1_fAsp")
      call var%getData(dummyd2)
      l1_fasp(s1 : e1, 1, 1) = pack(dummyd2, mask1)
 
      ! Hargreaves Samani coeffiecient
      var = nc%getVariable("L1_HarSamCoeff")
      call var%getData(dummyd2)
      l1_harsamcoeff(s1 : e1, 1, 1) = pack(dummyd2, mask1)
 
    case(2) ! Priestely-Taylor
 
      ! Priestley Taylor coeffiecient (alpha)
      var = nc%getVariable("L1_PrieTayAlpha")
      call var%getData(dummyd3)
      do ii = 1, nlai
        l1_prietayalpha(s1 : e1, ii, 1) = pack(dummyd3(:, :, ii), mask1)
      end do
 
    case(3) ! Penman-Monteith
 
      ! aerodynamical resitance
      var = nc%getVariable("L1_aeroResist")
      call var%getData(dummyd4)
      do jj = 1, nlcoverscene
        do ii = 1, nlai
          l1_aeroresist(s1 : e1, ii, jj) = pack(dummyd4(:, :, ii, jj), mask1)
        end do
      end do
 
      ! bulk surface resitance
      var = nc%getVariable("L1_surfResist")
      call var%getData(dummyd3)
      do ii = 1, nlai
        l1_surfresist(s1 : e1, ii, 1) = pack(dummyd3(:, :, ii), mask1)
      end do
 
    end select
 
 
   ! neutron count
   select case (processmatrix(10, 1))
   case(1) ! deslet
      ! N0 count
      var = nc%getVariable("L1_No_Count")
      call var%getData(dummyd2)
      l1_no_count(s1:e1, 1, 1) = pack(dummyd2, mask1)
 
      ! Bulk density
      var = nc%getVariable("L1_bulkDens")
      call var%getData(dummyd4)
      do jj = 1, nlcoverscene
         do ii = 1, nsoilhorizons_mhm
            l1_bulkdens(s1:e1, ii, jj) = pack(dummyd4(:, :, ii, jj), mask1)
         end do
      end do
 
      ! Lattice water
      var = nc%getVariable("L1_latticeWater")
      call var%getData(dummyd4)
      do jj = 1, nlcoverscene
         do ii = 1, nsoilhorizons_mhm
            l1_latticewater(s1:e1, ii, jj) = pack(dummyd4(:, :, ii, jj), mask1)
         end do
      end do
 
   case(2) ! COSMIC
      ! N0 count
      var = nc%getVariable("L1_No_Count")
      call var%getData(dummyd2)
      l1_no_count(s1 : e1, 1, 1) = pack(dummyd2, mask1)
 
      ! Bulk density
      var = nc%getVariable("L1_bulkDens")
      call var%getData(dummyd4)
      do jj = 1, nlcoverscene
         do ii = 1, nsoilhorizons_mhm
            l1_bulkdens(s1:e1, ii, jj) = pack(dummyd4(:, :, ii, jj), mask1)
         end do
      end do
 
      ! Lattice water
      var = nc%getVariable("L1_latticeWater")
      call var%getData(dummyd4)
      do jj = 1, nlcoverscene
         do ii = 1, nsoilhorizons_mhm
            l1_latticewater(s1: e1, ii, jj) = pack(dummyd4(:, :, ii, jj), mask1)
         end do
      end do
 
      ! COSMIC L3 parameter
      var = nc%getVariable("L1_COSMICL3")
      call var%getData(dummyd4)
      do jj = 1, nlcoverscene
         do ii = 1, nsoilhorizons_mhm
            l1_cosmicl3(s1:e1, ii, jj) = pack(dummyd4(:, :, ii, jj), mask1)
         end do
      end do
 
   end select
 
   ! close file
   call nc%close()
 
  end subroutine read_restart_states
 
END MODULE mo_restart