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

          Line data    Source code
       1             : !> \file mo_restart.f90
       2             : !> \brief   \copybrief mo_restart
       3             : !> \details \copydetails mo_restart
       4             : 
       5             : !> \brief reading and writing states, fluxes and configuration for restart of mHM.
       6             : !> \details routines are seperated for reading and writing variables for:
       7             : !!       - states and fluxes, and
       8             : !!       - configuration.
       9             : !!
      10             : !! Reading of L11 configuration is also seperated from the rest,
      11             : !! since it is only required when routing is activated.
      12             : !> \authors Stephan Thober
      13             : !> \date Jul 2013
      14             : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
      15             : !! mHM is released under the LGPLv3+ license \license_note
      16             : !> \ingroup f_mhm
      17             : MODULE mo_restart
      18             : 
      19             :   use mo_common_constants, only : soilHorizonsVarName, landCoverPeriodsVarName, LAIVarName
      20             : 
      21             :   IMPLICIT NONE
      22             : 
      23             :   PRIVATE
      24             : 
      25             :   PUBLIC :: read_restart_states     ! read restart files for state variables from a given path
      26             :   PUBLIC :: write_restart_files     ! write restart files for configuration to a given path
      27             : 
      28             : CONTAINS
      29             : 
      30             :   !> \brief write restart files for each domain
      31             :   !> \details write restart files for each domain. For each domain
      32             :   !! three restart files are written. These are xxx_states.nc,
      33             :   !! xxx_L11_config.nc, and xxx_config.nc (xxx being the three digit
      34             :   !! domain index). If a variable is added here, it should also be added
      35             :   !! in the read restart routines below.
      36             :   !> \changelog
      37             :   !! - Stephan Thober     Aug  2015
      38             :   !!   - moved write of routing states to mRM
      39             :   !! - David Schaefer     Nov  2015
      40             :   !!   - mo_netcdf
      41             :   !! - Stephan Thober     Nov  2016
      42             :   !!   - moved processMatrix to common variables
      43             :   !! - Zink M. Demirel C. Mar 2017
      44             :   !!   - Added Jarvis soil water stress function at SM process(3)
      45             :   !! - Robert Schweppe    Feb 2018
      46             :   !!   - Removed all L0 references
      47             :   !! - Robert Schweppe    Jun 2018
      48             :   !!   - refactoring and reformatting
      49             :   !! - Stephan Thober     Dec 2022
      50             :   !!   - added grid info for level0
      51             :   !> \authors Stephan Thober
      52             :   !> \date Jun 2014
      53          19 :   subroutine write_restart_files(OutFile)
      54             : 
      55             :     use mo_common_constants, only : nodata_dp
      56             :     use mo_common_restart, only : write_grid_info
      57             :     use mo_common_variables, only : level0, level1, nLCoverScene, domainMeta, LC_year_start, LC_year_end
      58             :     use mo_global_variables, only : L1_Inter, L1_Throughfall, L1_aETCanopy, L1_aETSealed, L1_aETSoil, L1_baseflow, &
      59             :                                     L1_fastRunoff, L1_infilSoil, L1_melt, L1_percol, L1_preEffect, L1_rain, &
      60             :                                     L1_runoffSeal, L1_satSTW, L1_sealSTW, L1_slowRunoff, L1_snow, L1_snowPack, &
      61             :                                     L1_soilMoist, L1_total_runoff, L1_unsatSTW
      62             : 
      63             :     use mo_kind, only : dp, i4
      64             :     use mo_message, only : message
      65             :     use mo_mpr_global_variables, only : nLAI, nSoilHorizons_mHM, HorizonDepth_mHM, LAIBoundaries
      66             :     use mo_mpr_restart, only : write_eff_params
      67             :     use mo_netcdf, only : NcDataset, NcDimension, NcVariable
      68             :     use mo_string_utils, only : num2str
      69             : 
      70             :     implicit none
      71             : 
      72             :     character(256) :: Fname
      73             : 
      74             :     ! Output Path for each domain
      75             :     character(256), dimension(:), intent(in) :: OutFile
      76             : 
      77             :     integer(i4) :: iDomain, domainID
      78             : 
      79             :     integer(i4) :: ii
      80             : 
      81             :     ! start index at level 1
      82             :     integer(i4) :: s1
      83             : 
      84             :     ! end index at level 1
      85             :     integer(i4) :: e1
      86             : 
      87             :     ! mask at level 1
      88           9 :     logical, dimension(:, :), allocatable :: mask1
      89             : 
      90             :     ! dummy variable
      91           9 :     real(dp), dimension(:, :, :), allocatable :: dummy_3D
      92           9 :     real(dp), dimension(:), allocatable :: dummy_1D
      93             : 
      94             :     integer(i4) :: max_extent
      95             : 
      96             :     type(NcDataset) :: nc
      97             : 
      98             :     type(NcDimension) :: rows1, cols1, soil1, lcscenes, lais
      99             : 
     100             :     type(NcVariable) :: var
     101             : 
     102             : 
     103             :     ! get maximum extent of one dimension 2 or 3
     104           9 :     max_extent = max(nSoilHorizons_mHM, nLCoverScene, nLAI)
     105             : 
     106          25 :     domain_loop : do iDomain = 1, domainMeta%nDomains
     107          16 :       domainID = domainMeta%indices(iDomain)
     108             : 
     109             :       ! write restart file for iDomain
     110          16 :       Fname = trim(OutFile(iDomain))
     111             :       ! print a message
     112          16 :       call message("    Writing Restart-file: ", trim(adjustl(Fname)), " ...")
     113             : 
     114          16 :       nc = NcDataset(fname, "w")
     115             : 
     116          16 :       call write_grid_info(level1(iDomain), "1", nc)
     117          16 :       call write_grid_info(level0(domainMeta%L0DataFrom(iDomain)), "0", nc)
     118             : 
     119          16 :       rows1 = nc%getDimension("nrows1")
     120          16 :       cols1 = nc%getDimension("ncols1")
     121             : 
     122             :       ! write the dimension to the file and also save bounds
     123          48 :       allocate(dummy_1D(nSoilHorizons_mHM+1))
     124          16 :       dummy_1D(1) = 0.0_dp
     125          48 :       dummy_1D(2:nSoilHorizons_mHM+1) = HorizonDepth_mHM(:)
     126          16 :       soil1 = nc%setCoordinate(trim(soilHorizonsVarName), nSoilHorizons_mHM, dummy_1D, 2_i4)
     127          16 :       deallocate(dummy_1D)
     128          48 :       allocate(dummy_1D(nLCoverScene+1))
     129          48 :       dummy_1D(1:nLCoverScene) = LC_year_start(:)
     130             :       ! this is done because bounds are always stored as real so e.g.
     131             :       ! 1981-1990,1991-2000 is thus saved as 1981.0-1991.0,1991.0-2001.0
     132             :       ! it is translated back into ints correctly during reading
     133          16 :       dummy_1D(nLCoverScene+1) = LC_year_end(nLCoverScene) + 1
     134          16 :       lcscenes = nc%setCoordinate(trim(landCoverPeriodsVarName), nLCoverScene, dummy_1D, 0_i4)
     135          16 :       deallocate(dummy_1D)
     136             :       ! write the dimension to the file
     137          16 :       lais = nc%setCoordinate(trim(LAIVarName), nLAI, LAIBoundaries, 0_i4)
     138             : 
     139             :       ! for appending and intialization
     140          80 :       allocate(dummy_3D(rows1%getLength(), cols1%getLength(), max_extent))
     141          64 :       allocate(mask1(rows1%getLength(), cols1%getLength()))
     142          16 :       s1 = level1(iDomain)%iStart
     143          16 :       e1 = level1(iDomain)%iEnd
     144        1316 :       mask1 = level1(iDomain)%mask
     145             : 
     146          48 :       var = nc%setVariable("L1_Inter", "f64", (/rows1, cols1/))
     147          16 :       call var%setFillValue(nodata_dp)
     148          16 :       call var%setData(unpack(L1_inter(s1 : e1), mask1, nodata_dp))
     149          16 :       call var%setAttribute("long_name", "Interception storage at level 1")
     150             : 
     151          48 :       var = nc%setVariable("L1_snowPack", "f64", (/rows1, cols1/))
     152          16 :       call var%setFillValue(nodata_dp)
     153          16 :       call var%setData(unpack(L1_snowPack(s1 : e1), mask1, nodata_dp))
     154          16 :       call var%setAttribute("long_name", "Snowpack at level 1")
     155             : 
     156          48 :       var = nc%setVariable("L1_sealSTW", "f64", (/rows1, cols1/))
     157          16 :       call var%setFillValue(nodata_dp)
     158          16 :       call var%setData(unpack(L1_sealSTW(s1 : e1), mask1, nodata_dp))
     159          16 :       call var%setAttribute("long_name", "Retention storage of impervious areas at level 1")
     160             : 
     161          48 :       do ii = 1, nSoilHorizons_mHM
     162          48 :         dummy_3D(:, :, ii) = unpack(L1_soilMoist(s1 : e1, ii), mask1, nodata_dp)
     163             :       end do
     164             : 
     165          64 :       var = nc%setVariable("L1_soilMoist", "f64", (/rows1, cols1, soil1/))
     166          16 :       call var%setFillValue(nodata_dp)
     167          16 :       call var%setData(dummy_3D(:, :, 1 : nSoilHorizons_mHM))
     168          16 :       call var%setAttribute("long_name", "soil moisture at level 1")
     169             : 
     170          48 :       var = nc%setVariable("L1_unsatSTW", "f64", (/rows1, cols1/))
     171          16 :       call var%setFillValue(nodata_dp)
     172          16 :       call var%setData(unpack(L1_unsatSTW(s1 : e1), mask1, nodata_dp))
     173          16 :       call var%setAttribute("long_name", "upper soil storage at level 1")
     174             : 
     175          48 :       var = nc%setVariable("L1_satSTW", "f64", (/rows1, cols1/))
     176          16 :       call var%setFillValue(nodata_dp)
     177          16 :       call var%setData(unpack(L1_satSTW(s1 : e1), mask1, nodata_dp))
     178          16 :       call var%setAttribute("long_name", "groundwater storage at level 1")
     179             : 
     180          48 :       do ii = 1, nSoilHorizons_mHM
     181          48 :         dummy_3D(:, :, ii) = unpack(L1_aETSoil(s1 : e1, ii), mask1, nodata_dp)
     182             :       end do
     183             : 
     184          64 :       var = nc%setVariable("L1_aETSoil", "f64", (/rows1, cols1, soil1/))
     185          16 :       call var%setFillValue(nodata_dp)
     186          16 :       call var%setData(dummy_3D(:, :, 1 : nSoilHorizons_mHM))
     187          16 :       call var%setAttribute("long_name", "soil actual ET at level 1")
     188             : 
     189          48 :       var = nc%setVariable("L1_aETCanopy", "f64", (/rows1, cols1/))
     190          16 :       call var%setFillValue(nodata_dp)
     191          16 :       call var%setData(unpack(L1_aETCanopy(s1 : e1), mask1, nodata_dp))
     192          16 :       call var%setAttribute("long_name", "canopy actual ET at level 1")
     193             : 
     194          48 :       var = nc%setVariable("L1_aETSealed", "f64", (/rows1, cols1/))
     195          16 :       call var%setFillValue(nodata_dp)
     196          16 :       call var%setData(unpack(L1_aETSealed(s1 : e1), mask1, nodata_dp))
     197          16 :       call var%setAttribute("long_name", "sealed actual ET at level 1")
     198             : 
     199          48 :       var = nc%setVariable("L1_baseflow", "f64", (/rows1, cols1/))
     200          16 :       call var%setFillValue(nodata_dp)
     201          16 :       call var%setData(unpack(L1_baseflow(s1 : e1), mask1, nodata_dp))
     202          16 :       call var%setAttribute("long_name", "baseflow at level 1")
     203             : 
     204          48 :       do ii = 1, nSoilHorizons_mHM
     205          48 :         dummy_3D(:, :, ii) = unpack(L1_infilSoil(s1 : e1, ii), mask1, nodata_dp)
     206             :       end do
     207             : 
     208          64 :       var = nc%setVariable("L1_infilSoil", "f64", (/rows1, cols1, soil1/))
     209          16 :       call var%setFillValue(nodata_dp)
     210          16 :       call var%setData(dummy_3D(:, :, 1 : nSoilHorizons_mHM))
     211          16 :       call var%setAttribute("long_name", "soil in-exfiltration at level 1")
     212             : 
     213          48 :       var = nc%setVariable("L1_fastRunoff", "f64", (/rows1, cols1/))
     214          16 :       call var%setFillValue(nodata_dp)
     215          16 :       call var%setData(unpack(L1_fastRunoff(s1 : e1), mask1, nodata_dp))
     216          16 :       call var%setAttribute("long_name", "fast runoff")
     217             : 
     218          48 :       var = nc%setVariable("L1_percol", "f64", (/rows1, cols1/))
     219          16 :       call var%setFillValue(nodata_dp)
     220          16 :       call var%setData(unpack(L1_percol(s1 : e1), mask1, nodata_dp))
     221          16 :       call var%setAttribute("long_name", "percolation at level 1")
     222             : 
     223          48 :       var = nc%setVariable("L1_melt", "f64", (/rows1, cols1/))
     224          16 :       call var%setFillValue(nodata_dp)
     225          16 :       call var%setData(unpack(L1_melt(s1 : e1), mask1, nodata_dp))
     226          16 :       call var%setAttribute("long_name", "snow melt at level 1")
     227             : 
     228          48 :       var = nc%setVariable("L1_preEffect", "f64", (/rows1, cols1/))
     229          16 :       call var%setFillValue(nodata_dp)
     230          16 :       call var%setData(unpack(L1_preEffect(s1 : e1), mask1, nodata_dp))
     231          16 :       call var%setAttribute("long_name", "effective precip. depth (snow melt + rain) at level 1")
     232             : 
     233          48 :       var = nc%setVariable("L1_rain", "f64", (/rows1, cols1/))
     234          16 :       call var%setFillValue(nodata_dp)
     235          16 :       call var%setData(unpack(L1_rain(s1 : e1), mask1, nodata_dp))
     236          16 :       call var%setAttribute("long_name", "rain (liquid water) at level 1")
     237             : 
     238          48 :       var = nc%setVariable("L1_runoffSeal", "f64", (/rows1, cols1/))
     239          16 :       call var%setFillValue(nodata_dp)
     240          16 :       call var%setData(unpack(L1_runoffSeal(s1 : e1), mask1, nodata_dp))
     241          16 :       call var%setAttribute("long_name", "runoff from impervious area at level 1")
     242             : 
     243          48 :       var = nc%setVariable("L1_slowRunoff", "f64", (/rows1, cols1/))
     244          16 :       call var%setFillValue(nodata_dp)
     245          16 :       call var%setData(unpack(L1_slowRunoff(s1 : e1), mask1, nodata_dp))
     246          16 :       call var%setAttribute("long_name", "slow runoff at level 1")
     247             : 
     248          48 :       var = nc%setVariable("L1_snow", "f64", (/rows1, cols1/))
     249          16 :       call var%setFillValue(nodata_dp)
     250          16 :       call var%setData(unpack(L1_snow(s1 : e1), mask1, nodata_dp))
     251          16 :       call var%setAttribute("long_name", "snow (solid water) at level 1")
     252             : 
     253          48 :       var = nc%setVariable("L1_Throughfall", "f64", (/rows1, cols1/))
     254          16 :       call var%setFillValue(nodata_dp)
     255          16 :       call var%setData(unpack(L1_Throughfall(s1 : e1), mask1, nodata_dp))
     256          16 :       call var%setAttribute("long_name", "throughfall at level 1")
     257             : 
     258          48 :       var = nc%setVariable("L1_total_runoff", "f64", (/rows1, cols1/))
     259          16 :       call var%setFillValue(nodata_dp)
     260          16 :       call var%setData(unpack(L1_total_runoff(s1 : e1), mask1, nodata_dp))
     261          16 :       call var%setAttribute("long_name", "total runoff at level 1")
     262             : 
     263          16 :       call write_eff_params(mask1, s1, e1, rows1, cols1, soil1, lcscenes, lais, nc)
     264             : 
     265          16 :       call nc%close()
     266             : 
     267          57 :       deallocate(dummy_3D, mask1)
     268             :     end do domain_loop
     269             : 
     270           9 :   end subroutine write_restart_files
     271             : 
     272             : 
     273             :   !> \brief reads fluxes and state variables from file
     274             :   !> \details read fluxes and state variables from given
     275             :   !! restart directory and initialises all state variables
     276             :   !! that are initialized in the subroutine initialise,
     277             :   !! contained in module mo_startup.
     278             :   !> \changelog
     279             :   !! - Stephan Thober Aug  2015
     280             :   !!   - moved read of routing states to mRM
     281             :   !! - David Schaefer Nov  2015
     282             :   !!   - mo_netcdf
     283             :   !! - Stephan Thober Nov  2016
     284             :   !!   - moved processMatrix to common variables
     285             :   !! - Robert Schweppe Jun 2018
     286             :   !!   - refactoring and reformatting
     287             :   !! - Sebastian Müller Mar 2023
     288             :   !!   - compatibility layer for 2D/3D data
     289             :   !!   - move reading of nLAI to mo_startup (needed beforehand)
     290             :   !> \authors Stephan Thober
     291             :   !> \date Apr 2013
     292           1 :   subroutine read_restart_states(iDomain, domainID, InFile)
     293             : 
     294           9 :     use mo_common_variables, only : LC_year_end, LC_year_start, level1, nLCoverScene, processMatrix
     295             :     use mo_global_variables, only : L1_Inter, L1_Throughfall, L1_aETCanopy, &
     296             :                                     L1_aETSealed, L1_aETSoil, L1_baseflow, L1_fastRunoff, L1_infilSoil, L1_melt, &
     297             :                                     L1_percol, L1_preEffect, L1_rain, L1_runoffSeal, L1_satSTW, L1_sealSTW, &
     298             :                                     L1_slowRunoff, L1_snow, L1_snowPack, L1_soilMoist, L1_total_runoff, L1_unsatSTW
     299             :     use mo_kind, only : dp, i4
     300             :     use mo_mpr_global_variables, only : L1_HarSamCoeff, &
     301             :                                         L1_PrieTayAlpha, L1_aeroResist, L1_alpha, L1_degDay, L1_degDayInc, L1_degDayMax, &
     302             :                                         L1_degDayNoPre, L1_fAsp, L1_fRoots, L1_fSealed, L1_jarvis_thresh_c1, &
     303             :                                         L1_kBaseFlow, L1_kPerco, L1_kSlowFlow, L1_karstLoss, L1_kfastFlow, L1_maxInter, &
     304             :                                         L1_petLAIcorFactor, L1_sealedThresh, L1_soilMoistExp, L1_soilMoistFC, &
     305             :                                         L1_soilMoistSat, L1_surfResist, L1_tempThresh, L1_unsatThresh, L1_wiltingPoint, &
     306             :                                         nLAI, nSoilHorizons_mHM, &
     307             :                                         ! neutron count
     308             :                                         L1_No_Count, L1_bulkDens, L1_latticeWater, L1_COSMICL3
     309             : 
     310             :     use mo_netcdf, only : NcDataset, NcDimension, NcVariable
     311             :     use mo_string_utils, only : num2str
     312             :     use mo_message, only: message, error_message
     313             : 
     314             :     implicit none
     315             : 
     316             :     !> number of domain
     317             :     integer(i4), intent(in) :: iDomain
     318             :     !> ID of domain
     319             :     integer(i4), intent(in) :: domainID
     320             :     !> Input Path including trailing slash
     321             :     character(256), intent(in) :: InFile
     322             : 
     323             :     character(256) :: Fname
     324             :     ! variable rank
     325             :     integer(i4) :: var_rank
     326             :     ! loop index
     327             :     integer(i4) :: ii, jj
     328             :     ! start index at level 1
     329             :     integer(i4) :: s1
     330             :     ! end index at level 1
     331             :     integer(i4) :: e1
     332             :     ! mask at level 1
     333           1 :     logical, dimension(:, :), allocatable :: mask1
     334             :     ! dummy, 2 dimension
     335           1 :     real(dp), dimension(:, :), allocatable :: dummyD2
     336             :     ! dummy, 3 dimension
     337           1 :     real(dp), dimension(:, :, :), allocatable :: dummyD3
     338             :     ! dummy, 4 dimension
     339           1 :     real(dp), dimension(:, :, :, :), allocatable :: dummyD4
     340             : 
     341             :     type(NcDataset) :: nc
     342             :     type(NcVariable) :: var
     343             : 
     344             : 
     345           1 :     Fname = trim(InFile)
     346           1 :     call message('    Reading states from ', trim(adjustl(Fname)),' ...')
     347             : 
     348             :     ! get domain information at level 1
     349           4 :     allocate(mask1 (level1(iDomain)%nrows, level1(iDomain)%ncols))
     350          65 :     mask1 = level1(iDomain)%mask
     351           1 :     s1 = level1(iDomain)%iStart
     352           1 :     e1 = level1(iDomain)%iEnd
     353             : 
     354           1 :     nc = NcDataset(fname, "r")
     355             : 
     356           1 :     if (nc%hasVariable('L1_Inter')) then
     357             :       !-------------------------------------------
     358             :       ! STATE VARIABLES (optionally)
     359             :       !-------------------------------------------
     360             : 
     361             :       ! Interception
     362           1 :       var = nc%getVariable("L1_Inter")
     363           1 :       call var%getData(dummyD2)
     364           1 :       L1_inter(s1 : e1) = pack(dummyD2, mask1)
     365             : 
     366             :       ! Snowpack
     367           1 :       var = nc%getVariable("L1_snowPack")
     368           1 :       call var%getData(dummyD2)
     369           1 :       L1_snowPack(s1 : e1) = pack(dummyD2, mask1)
     370             : 
     371             :       ! Retention storage of impervious areas
     372           1 :       var = nc%getVariable("L1_sealSTW")
     373           1 :       call var%getData(dummyD2)
     374           1 :       L1_sealSTW(s1 : e1) = pack(dummyD2, mask1)
     375             : 
     376             :       ! upper soil storage
     377           1 :       var = nc%getVariable("L1_unsatSTW")
     378           1 :       call var%getData(dummyD2)
     379           1 :       L1_unsatSTW(s1 : e1) = pack(dummyD2, mask1)
     380             : 
     381             :       ! groundwater storage
     382           1 :       var = nc%getVariable("L1_satSTW")
     383           1 :       call var%getData(dummyD2)
     384           1 :       L1_satSTW(s1 : e1) = pack(dummyD2, mask1)
     385             : 
     386             :       ! Soil moisture of each horizon
     387           1 :       var = nc%getVariable("L1_soilMoist")
     388           1 :       call var%getData(dummyD3)
     389           3 :       do ii = 1, nSoilHorizons_mHM
     390           3 :         L1_soilMoist(s1 : e1, ii) = pack(dummyD3(:, :, ii), mask1)
     391             :       end do
     392             : 
     393             :       !-------------------------------------------
     394             :       ! FLUXES (optionally)
     395             :       !-------------------------------------------
     396             : 
     397             :       !  soil actual ET
     398           1 :       var = nc%getVariable("L1_aETSoil")
     399           1 :       call var%getData(dummyD3)
     400           3 :       do ii = 1, nSoilHorizons_mHM
     401           3 :         L1_aETSoil(s1 : e1, ii) = pack(dummyD3(:, :, ii), mask1)
     402             :       end do
     403             : 
     404             :       ! canopy actual ET
     405           1 :       var = nc%getVariable("L1_aETCanopy")
     406           1 :       call var%getData(dummyD2)
     407           1 :       L1_aETCanopy(s1 : e1) = pack(dummyD2, mask1)
     408             : 
     409             :       ! sealed area actual ET
     410           1 :       var = nc%getVariable("L1_aETSealed")
     411           1 :       call var%getData(dummyD2)
     412           1 :       L1_aETSealed(s1 : e1) = pack(dummyD2, mask1)
     413             : 
     414             :       ! baseflow
     415           1 :       var = nc%getVariable("L1_baseflow")
     416           1 :       call var%getData(dummyD2)
     417           1 :       L1_baseflow(s1 : e1) = pack(dummyD2, mask1)
     418             : 
     419             :       ! soil in-exfiltration
     420           1 :       var = nc%getVariable("L1_infilSoil")
     421           1 :       call var%getData(dummyD3)
     422           3 :       do ii = 1, nSoilHorizons_mHM
     423           3 :         L1_infilSoil(s1 : e1, ii) = pack(dummyD3(:, :, ii), mask1)
     424             :       end do
     425             : 
     426             :       ! fast runoff
     427           1 :       var = nc%getVariable("L1_fastRunoff")
     428           1 :       call var%getData(dummyD2)
     429           1 :       L1_fastRunoff(s1 : e1) = pack(dummyD2, mask1)
     430             : 
     431             :       ! snow melt
     432           1 :       var = nc%getVariable("L1_melt")
     433           1 :       call var%getData(dummyD2)
     434           1 :       L1_melt(s1 : e1) = pack(dummyD2, mask1)
     435             : 
     436             :       ! percolation
     437           1 :       var = nc%getVariable("L1_percol")
     438           1 :       call var%getData(dummyD2)
     439           1 :       L1_percol(s1 : e1) = pack(dummyD2, mask1)
     440             : 
     441             :       ! effective precip. depth (snow melt + rain)
     442           1 :       var = nc%getVariable("L1_preEffect")
     443           1 :       call var%getData(dummyD2)
     444           1 :       L1_preEffect(s1 : e1) = pack(dummyD2, mask1)
     445             : 
     446             :       ! rain (liquid water)
     447           1 :       var = nc%getVariable("L1_rain")
     448           1 :       call var%getData(dummyD2)
     449           1 :       L1_rain(s1 : e1) = pack(dummyD2, mask1)
     450             : 
     451             :       ! runoff from impervious area
     452           1 :       var = nc%getVariable("L1_runoffSeal")
     453           1 :       call var%getData(dummyD2)
     454           1 :       L1_runoffSeal(s1 : e1) = pack(dummyD2, mask1)
     455             : 
     456             :       ! slow runoff
     457           1 :       var = nc%getVariable("L1_slowRunoff")
     458           1 :       call var%getData(dummyD2)
     459           1 :       L1_slowRunoff(s1 : e1) = pack(dummyD2, mask1)
     460             : 
     461             :       ! snow (solid water)
     462           1 :       var = nc%getVariable("L1_snow")
     463           1 :       call var%getData(dummyD2)
     464           1 :       L1_snow(s1 : e1) = pack(dummyD2, mask1)
     465             : 
     466             :       ! throughfall
     467           1 :       var = nc%getVariable("L1_Throughfall")
     468           1 :       call var%getData(dummyD2)
     469           1 :       L1_Throughfall(s1 : e1) = pack(dummyD2, mask1)
     470             : 
     471             :       ! total runoff
     472           1 :       var = nc%getVariable("L1_total_runoff")
     473           1 :       call var%getData(dummyD2)
     474           1 :       L1_total_runoff(s1 : e1) = pack(dummyD2, mask1)
     475             :     end if
     476             : 
     477             :     !-------------------------------------------
     478             :     ! EFFECTIVE PARAMETERS
     479             :     !-------------------------------------------
     480           1 :     var = nc%getVariable("L1_fSealed")
     481           1 :     call var%getData(dummyD3)
     482           3 :     do ii = 1, nLCoverScene
     483           3 :       L1_fSealed(s1 : e1, 1, ii) = pack(dummyD3(:, :, ii), mask1)
     484             :     end do
     485             : 
     486             :     ! exponent for the upper reservoir
     487           1 :     var = nc%getVariable("L1_alpha")
     488           1 :     var_rank = var%getrank()
     489           0 :     select case(var_rank)
     490             :       case(2)
     491             :         ! distribute over all land cover scenes
     492           0 :         call var%getData(dummyD2)
     493           0 :         do ii = 1, nLCoverScene
     494           0 :           L1_alpha(s1 : e1, 1, ii) = pack(dummyD2, mask1)
     495             :         end do
     496             :       case(3)
     497           1 :         call var%getData(dummyD3)
     498           3 :         do ii = 1, nLCoverScene
     499           3 :           L1_alpha(s1 : e1, 1, ii) = pack(dummyD3(:, :, ii), mask1)
     500             :         end do
     501             :       case default
     502           1 :         call error_message("Restart: L1_alpha rank needs to be 2 or 3")
     503             :     end select
     504             : 
     505             :     ! increase of the Degree-day factor per mm of increase in precipitation
     506           1 :     var = nc%getVariable("L1_degDayInc")
     507           1 :     call var%getData(dummyD3)
     508           3 :     do ii = 1, nLCoverScene
     509           3 :       L1_degDayInc(s1 : e1, 1, ii) = pack(dummyD3(:, :, ii), mask1)
     510             :     end do
     511             : 
     512             :     ! maximum degree-day factor
     513           1 :     var = nc%getVariable("L1_degDayMax")
     514           1 :     call var%getData(dummyD3)
     515           3 :     do ii = 1, nLCoverScene
     516           3 :       L1_degDayMax(s1 : e1, 1, ii) = pack(dummyD3(:, :, ii), mask1)
     517             :     end do
     518             : 
     519             :     ! degree-day factor with no precipitation
     520           1 :     var = nc%getVariable("L1_degDayNoPre")
     521           1 :     call var%getData(dummyD3)
     522           3 :     do ii = 1, nLCoverScene
     523           3 :       L1_degDayNoPre(s1 : e1, 1, ii) = pack(dummyD3(:, :, ii), mask1)
     524             :     end do
     525             : 
     526             :     ! degree-day factor
     527           1 :     var = nc%getVariable("L1_degDay")
     528           1 :     var_rank = var%getrank()
     529           0 :     select case(var_rank)
     530             :       case(2)
     531             :         ! distribute over all land cover scenes
     532           0 :         call var%getData(dummyD2)
     533           0 :         do ii = 1, nLCoverScene
     534           0 :           L1_degDay(s1 : e1, 1, ii) = pack(dummyD2, mask1)
     535             :         end do
     536             :       case(3)
     537           1 :         call var%getData(dummyD3)
     538           3 :         do ii = 1, nLCoverScene
     539           3 :           L1_degDay(s1 : e1, 1, ii) = pack(dummyD3(:, :, ii), mask1)
     540             :         end do
     541             :       case default
     542           1 :         call error_message("Restart: L1_degDay rank needs to be 2 or 3")
     543             :     end select
     544             : 
     545             :     ! Karstic percolation loss
     546           1 :     var = nc%getVariable("L1_karstLoss")
     547           1 :     call var%getData(dummyD2)
     548           1 :     L1_karstLoss(s1 : e1, 1, 1) = pack(dummyD2, mask1)
     549             : 
     550             :     ! Fraction of roots in soil horizons
     551           1 :     var = nc%getVariable("L1_fRoots")
     552           1 :     call var%getData(dummyD4)
     553           3 :     do jj = 1, nLCoverScene
     554           7 :       do ii = 1, nSoilHorizons_mHM
     555           6 :         L1_fRoots(s1 : e1, ii, jj) = pack(dummyD4(:, :, ii, jj), mask1)
     556             :       end do
     557             :     end do
     558             : 
     559             :     ! Maximum interception
     560           1 :     var = nc%getVariable("L1_maxInter")
     561           1 :     call var%getData(dummyD3)
     562          13 :     do ii = 1, nLAI
     563          13 :       L1_maxInter(s1 : e1, ii, 1) = pack(dummyD3(:, :, ii), mask1)
     564             :     end do
     565             : 
     566             :     ! fast interflow recession coefficient ("L1_kFastFlow" or "L1_kfastFlow")
     567           1 :     if (nc%hasvariable("L1_kFastFlow")) then
     568           0 :       var = nc%getVariable("L1_kFastFlow")
     569             :     else
     570           1 :       var = nc%getVariable("L1_kfastFlow")
     571             :     end if
     572           1 :     call var%getData(dummyD3)
     573           3 :     do ii = 1, nLCoverScene
     574           3 :       L1_kfastFlow(s1 : e1, 1, ii) = pack(dummyD3(:, :, ii), mask1)
     575             :     end do
     576             : 
     577             :     ! slow interflow recession coefficient
     578           1 :     var = nc%getVariable("L1_kSlowFlow")
     579           1 :     var_rank = var%getrank()
     580           0 :     select case(var_rank)
     581             :       case(2)
     582             :         ! distribute over all land cover scenes
     583           0 :         call var%getData(dummyD2)
     584           0 :         do ii = 1, nLCoverScene
     585           0 :           L1_kSlowFlow(s1 : e1, 1, ii) = pack(dummyD2, mask1)
     586             :         end do
     587             :       case(3)
     588           1 :         call var%getData(dummyD3)
     589           3 :         do ii = 1, nLCoverScene
     590           3 :           L1_kSlowFlow(s1 : e1, 1, ii) = pack(dummyD3(:, :, ii), mask1)
     591             :         end do
     592             :       case default
     593           1 :         call error_message("Restart: L1_kSlowFlow rank needs to be 2 or 3")
     594             :     end select
     595             : 
     596             :     ! baseflow recession coefficient
     597           1 :     var = nc%getVariable("L1_kBaseFlow")
     598           1 :     var_rank = var%getrank()
     599           0 :     select case(var_rank)
     600             :       case(2)
     601             :         ! distribute over all land cover scenes
     602           0 :         call var%getData(dummyD2)
     603           0 :         do ii = 1, nLCoverScene
     604           0 :           L1_kBaseFlow(s1 : e1, 1, ii) = pack(dummyD2, mask1)
     605             :         end do
     606             :       case(3)
     607           1 :         call var%getData(dummyD3)
     608           3 :         do ii = 1, nLCoverScene
     609           3 :           L1_kBaseFlow(s1 : e1, 1, ii) = pack(dummyD3(:, :, ii), mask1)
     610             :         end do
     611             :       case default
     612           1 :         call error_message("Restart: L1_kBaseFlow rank needs to be 2 or 3")
     613             :     end select
     614             : 
     615             :     ! percolation coefficient
     616           1 :     var = nc%getVariable("L1_kPerco")
     617           1 :     var_rank = var%getrank()
     618           0 :     select case(var_rank)
     619             :       case(2)
     620             :         ! distribute over all land cover scenes
     621           0 :         call var%getData(dummyD2)
     622           0 :         do ii = 1, nLCoverScene
     623           0 :           L1_kPerco(s1 : e1, 1, ii) = pack(dummyD2, mask1)
     624             :         end do
     625             :       case(3)
     626           1 :         call var%getData(dummyD3)
     627           3 :         do ii = 1, nLCoverScene
     628           3 :           L1_kPerco(s1 : e1, 1, ii) = pack(dummyD3(:, :, ii), mask1)
     629             :         end do
     630             :       case default
     631           1 :         call error_message("Restart: L1_kPerco rank needs to be 2 or 3")
     632             :     end select
     633             : 
     634             :     ! Soil moisture below which actual ET is reduced linearly till PWP
     635             :     ! for processCase(3) = 1
     636           1 :     var = nc%getVariable("L1_soilMoistFC")
     637           1 :     call var%getData(dummyD4)
     638           3 :     do jj = 1, nLCoverScene
     639           7 :       do ii = 1, nSoilHorizons_mHM
     640           6 :         L1_soilMoistFC(s1 : e1, ii, jj) = pack(dummyD4(:, :, ii, jj), mask1)
     641             :       end do
     642             :     end do
     643             : 
     644             :     ! Saturation soil moisture for each horizon [mm]
     645           1 :     var = nc%getVariable("L1_soilMoistSat")
     646           1 :     call var%getData(dummyD4)
     647           3 :     do jj = 1, nLCoverScene
     648           7 :       do ii = 1, nSoilHorizons_mHM
     649           6 :         L1_soilMoistSat(s1 : e1, ii, jj) = pack(dummyD4(:, :, ii, jj), mask1)
     650             :       end do
     651             :     end do
     652             : 
     653             :     ! Exponential parameter to how non-linear is the soil water retention
     654           1 :     var = nc%getVariable("L1_soilMoistExp")
     655           1 :     call var%getData(dummyD4)
     656           3 :     do jj = 1, nLCoverScene
     657           7 :       do ii = 1, nSoilHorizons_mHM
     658           6 :         L1_soilMoistExp(s1 : e1, ii, jj) = pack(dummyD4(:, :, ii, jj), mask1)
     659             :       end do
     660             :     end do
     661             : 
     662           3 :     if (any(processMatrix(3, 1) == (/2, 3/))) then
     663             :       ! jarvis critical value for normalized soil water content
     664           0 :       var = nc%getVariable("L1_jarvis_thresh_c1")
     665           0 :       call var%getData(dummyD2)
     666           0 :       L1_jarvis_thresh_c1(s1 : e1, 1, 1) = pack(dummyD2, mask1)
     667             :     end if
     668             : 
     669             :     ! Threshold temperature for snow/rain
     670           1 :     var = nc%getVariable("L1_tempThresh")
     671           1 :     call var%getData(dummyD3)
     672           3 :     do ii = 1, nLCoverScene
     673           3 :       L1_tempThresh(s1 : e1, 1, ii) = pack(dummyD3(:, :, ii), mask1)
     674             :     end do
     675             : 
     676             :     ! Threshold water depth controlling fast interflow
     677           1 :     var = nc%getVariable("L1_unsatThresh")
     678           1 :     call var%getData(dummyD2)
     679           1 :     L1_unsatThresh(s1 : e1, 1, 1) = pack(dummyD2, mask1)
     680             : 
     681             :     ! Threshold water depth for surface runoff in sealed surfaces
     682           1 :     var = nc%getVariable("L1_sealedThresh")
     683           1 :     call var%getData(dummyD2)
     684           1 :     L1_sealedThresh(s1 : e1, 1, 1) = pack(dummyD2, mask1)
     685             : 
     686             :     ! Permanent wilting point
     687           1 :     var = nc%getVariable("L1_wiltingPoint")
     688           1 :     call var%getData(dummyD4)
     689           3 :     do jj = 1, nLCoverScene
     690           7 :       do ii = 1, nSoilHorizons_mHM
     691           6 :         L1_wiltingPoint(s1 : e1, ii, jj) = pack(dummyD4(:, :, ii, jj), mask1)
     692             :       end do
     693             :     end do
     694             : 
     695             :     ! different parameters dependent on PET formulation
     696           1 :     select case (processMatrix(5, 1))
     697             :     case(-1) ! PET is input
     698             : 
     699             :       ! PET correction factor due to LAI
     700           0 :       var = nc%getVariable("L1_petLAIcorFactor")
     701           0 :       call var%getData(dummyD4)
     702           0 :       do jj = 1, nLCoverScene
     703           0 :         do ii = 1, nLAI
     704           0 :           L1_petLAIcorFactor(s1 : e1, ii, jj) = pack(dummyD4(:, :, ii, jj), mask1)
     705             :         end do
     706             :       end do
     707             : 
     708             :     case(0) ! PET is input
     709             : 
     710             :       ! PET correction factor due to terrain aspect
     711           0 :       var = nc%getVariable("L1_fAsp")
     712           0 :       call var%getData(dummyD2)
     713           0 :       L1_fAsp(s1 : e1, 1, 1) = pack(dummyD2, mask1)
     714             : 
     715             :     case(1) ! Hargreaves-Samani
     716             : 
     717             :       ! PET correction factor due to terrain aspect
     718           1 :       var = nc%getVariable("L1_fAsp")
     719           1 :       call var%getData(dummyD2)
     720           1 :       L1_fAsp(s1 : e1, 1, 1) = pack(dummyD2, mask1)
     721             : 
     722             :       ! Hargreaves Samani coeffiecient
     723           1 :       var = nc%getVariable("L1_HarSamCoeff")
     724           1 :       call var%getData(dummyD2)
     725           1 :       L1_HarSamCoeff(s1 : e1, 1, 1) = pack(dummyD2, mask1)
     726             : 
     727             :     case(2) ! Priestely-Taylor
     728             : 
     729             :       ! Priestley Taylor coeffiecient (alpha)
     730           0 :       var = nc%getVariable("L1_PrieTayAlpha")
     731           0 :       call var%getData(dummyD3)
     732           0 :       do ii = 1, nLAI
     733           0 :         L1_PrieTayAlpha(s1 : e1, ii, 1) = pack(dummyD3(:, :, ii), mask1)
     734             :       end do
     735             : 
     736             :     case(3) ! Penman-Monteith
     737             : 
     738             :       ! aerodynamical resitance
     739           0 :       var = nc%getVariable("L1_aeroResist")
     740           0 :       call var%getData(dummyD4)
     741           0 :       do jj = 1, nLCoverScene
     742           0 :         do ii = 1, nLAI
     743           0 :           L1_aeroResist(s1 : e1, ii, jj) = pack(dummyD4(:, :, ii, jj), mask1)
     744             :         end do
     745             :       end do
     746             : 
     747             :       ! bulk surface resitance
     748           0 :       var = nc%getVariable("L1_surfResist")
     749           0 :       call var%getData(dummyD3)
     750           1 :       do ii = 1, nLAI
     751           0 :         L1_surfResist(s1 : e1, ii, 1) = pack(dummyD3(:, :, ii), mask1)
     752             :       end do
     753             : 
     754             :     end select
     755             : 
     756             : 
     757             :    ! neutron count
     758           1 :    select case (processMatrix(10, 1))
     759             :    case(1) ! deslet
     760             :       ! N0 count
     761           0 :       var = nc%getVariable("L1_No_Count")
     762           0 :       call var%getData(dummyD2)
     763           0 :       L1_No_Count(s1:e1, 1, 1) = pack(dummyD2, mask1)
     764             : 
     765             :       ! Bulk density
     766           0 :       var = nc%getVariable("L1_bulkDens")
     767           0 :       call var%getData(dummyD4)
     768           0 :       do jj = 1, nLCoverScene
     769           0 :          do ii = 1, nSoilHorizons_mHM
     770           0 :             L1_bulkDens(s1:e1, ii, jj) = pack(dummyD4(:, :, ii, jj), mask1)
     771             :          end do
     772             :       end do
     773             : 
     774             :       ! Lattice water
     775           0 :       var = nc%getVariable("L1_latticeWater")
     776           0 :       call var%getData(dummyD4)
     777           0 :       do jj = 1, nLCoverScene
     778           0 :          do ii = 1, nSoilHorizons_mHM
     779           0 :             L1_latticeWater(s1:e1, ii, jj) = pack(dummyD4(:, :, ii, jj), mask1)
     780             :          end do
     781             :       end do
     782             : 
     783             :    case(2) ! COSMIC
     784             :       ! N0 count
     785           0 :       var = nc%getVariable("L1_No_Count")
     786           0 :       call var%getData(dummyD2)
     787           0 :       L1_No_Count(s1 : e1, 1, 1) = pack(dummyD2, mask1)
     788             : 
     789             :       ! Bulk density
     790           0 :       var = nc%getVariable("L1_bulkDens")
     791           0 :       call var%getData(dummyD4)
     792           0 :       do jj = 1, nLCoverScene
     793           0 :          do ii = 1, nSoilHorizons_mHM
     794           0 :             L1_bulkDens(s1:e1, ii, jj) = pack(dummyD4(:, :, ii, jj), mask1)
     795             :          end do
     796             :       end do
     797             : 
     798             :       ! Lattice water
     799           0 :       var = nc%getVariable("L1_latticeWater")
     800           0 :       call var%getData(dummyD4)
     801           0 :       do jj = 1, nLCoverScene
     802           0 :          do ii = 1, nSoilHorizons_mHM
     803           0 :             L1_latticeWater(s1: e1, ii, jj) = pack(dummyD4(:, :, ii, jj), mask1)
     804             :          end do
     805             :       end do
     806             : 
     807             :       ! COSMIC L3 parameter
     808           0 :       var = nc%getVariable("L1_COSMICL3")
     809           0 :       call var%getData(dummyD4)
     810           1 :       do jj = 1, nLCoverScene
     811           0 :          do ii = 1, nSoilHorizons_mHM
     812           0 :             L1_COSMICL3(s1:e1, ii, jj) = pack(dummyD4(:, :, ii, jj), mask1)
     813             :          end do
     814             :       end do
     815             : 
     816             :    end select
     817             : 
     818             :    ! close file
     819           1 :    call nc%close()
     820             : 
     821           1 :   end subroutine read_restart_states
     822             : 
     823             : END MODULE mo_restart

Generated by: LCOV version 1.16