mo_common_restart.f90

Go to the documentation of this file.

!> \file mo_common_restart.f90
!> \brief   \copybrief mo_common_read_data
!> \details \copydetails mo_common_read_data
 
!> \brief common restart tools
!> \details Routines to deal with grid infos for restart files
!> \authors Robert Schweppe
!> \date Jun 2018
!> \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_common
module mo_common_restart
 
  use mo_kind, only : i4, dp
  use mo_message, only: message, error_message
 
  IMPLICIT NONE
 
  PRIVATE
 
  PUBLIC :: write_grid_info
  PUBLIC :: read_grid_info     ! read restart files for configuration from a given path
  PUBLIC :: read_nlai_and_check_dims
 
 
  !> \brief check consistency of two given items
  INTERFACE check_consistency_element
    MODULE PROCEDURE check_consistency_element_i4, check_consistency_element_dp
  end interface check_consistency_element
 
 
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     May 2019
  !!   - where statement for gnu73 to translate level0 mask
  !> \authors Stephan Thober
  !> \date Jun 2014
  subroutine write_grid_info(grid_in, level_name, nc)
 
    use mo_common_constants, only : nodata_dp, nodata_i4
    use mo_common_types, only: grid
    use mo_netcdf, only : ncdataset, ncdimension, ncvariable
 
    implicit none
 
    !> level to be written
    type(grid), intent(in) :: grid_in
 
    !> level_id
    character(*), intent(in) :: level_name
 
    !> NcDataset to write information to
    type(ncdataset), intent(inout) :: nc
 
    ! dummy for gnu73
    integer(i4), allocatable :: dummy(:, :)
 
    type(ncdimension) :: rows, cols
 
    type(ncvariable) :: var
 
 
    rows = nc%setDimension("nrows" // trim(level_name), grid_in%nrows)
    cols = nc%setDimension("ncols" // trim(level_name), grid_in%ncols)
 
    ! now set everything related to the grid
    var = nc%setVariable("L" // trim(level_name) // "_domain_mask", "i32", (/rows, cols/))
    call var%setFillValue(nodata_i4)
    ! transform from logical to i32
    ! ST: where statement is used because gnu73 does not properly translate with merge
    allocate(dummy(size(grid_in%mask, 1), size(grid_in%mask, 2)))
    dummy = 0_i4
    where(grid_in%mask) dummy = 1_i4
    call var%setData(dummy)
    deallocate(dummy)
    call var%setAttribute("long_name", "Mask at level " // trim(level_name))
 
    var = nc%setVariable("L" // trim(level_name) // "_domain_lat", "f64", (/rows, cols/))
    call var%setFillValue(nodata_dp)
    call var%setData(grid_in%y)
    call var%setAttribute("long_name", "Latitude at level " // trim(level_name))
 
    var = nc%setVariable("L" // trim(level_name) // "_domain_lon", "f64", (/rows, cols/))
    call var%setFillValue(nodata_dp)
    call var%setData(grid_in%x)
    call var%setAttribute("long_name", "Longitude at level " // trim(level_name))
 
    var = nc%setVariable("L" // trim(level_name) // "_domain_cellarea", "f64", (/rows, cols/))
    call var%setFillValue(nodata_dp)
    call var%setData(unpack(grid_in%CellArea * 1.0e-6_dp, grid_in%mask, nodata_dp))
    call var%setAttribute("long_name", "Cell area at level " // trim(level_name))
 
    call nc%setAttribute("xllcorner_L" // trim(level_name), grid_in%xllcorner)
    call nc%setAttribute("yllcorner_L" // trim(level_name), grid_in%yllcorner)
    call nc%setAttribute("cellsize_L" // trim(level_name), grid_in%cellsize)
    call nc%setAttribute("nrows_L" // trim(level_name), grid_in%nrows)
    call nc%setAttribute("ncols_L" // trim(level_name), grid_in%ncols)
    call nc%setAttribute("nCells_L" // trim(level_name), grid_in%nCells)
 
  end subroutine write_grid_info
 
 
  !> \brief reads configuration apart from Level 11 configuration from a restart directory
  !> \details read configuration variables from a given restart
  !> directory and initializes all configuration variables,
  !> that are initialized in the subroutine initialise,
  !> contained in module mo_startup.
  !> \changelog
  !! - David Schaefer     Nov 2015
  !!   - mo_netcdf
  !! - 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     May 2019
  !!   - added allocation check for mask and cellArea because cellArea needs to be read by mRM, but mask is created before by mHM
  !> \authors Stephan Thober
  !> \date Apr 2013
  subroutine read_grid_info(InFile, level_name, new_grid)
 
    use mo_common_types, only: grid
    use mo_netcdf, only : ncdataset, ncvariable
 
    implicit none
 
    !> Input Path including trailing slash
    character(256), intent(in) :: infile
 
    !> level_name (id)
    character(*), intent(in) :: level_name
 
    !> grid to save information to
    type(grid), intent(inout) :: new_grid
 
    ! dummy, 2 dimension I4
    integer(i4), dimension(:, :), allocatable :: dummyi2
 
    ! dummy, 2 dimension DP
    real(dp), dimension(:, :), allocatable :: dummyd2
 
    character(256) :: fname
 
    type(ncdataset) :: nc
 
    type(ncvariable) :: var
 
    integer(i4) :: k
 
 
    ! read config
    fname = trim(infile)
    call message('    Reading config from     ', trim(adjustl(fname)), ' ...')
 
    nc = ncdataset(fname, "r")
 
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! Read L1 variables <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! read the grid properties
    call nc%getAttribute("xllcorner_L" // trim(level_name), new_grid%xllcorner)
    call nc%getAttribute("yllcorner_L" // trim(level_name), new_grid%yllcorner)
    call nc%getAttribute("nrows_L" // trim(level_name), new_grid%nrows)
    call nc%getAttribute("ncols_L" // trim(level_name), new_grid%ncols)
    call nc%getAttribute("cellsize_L" // trim(level_name), new_grid%cellsize)
    call nc%getAttribute("nCells_L" // trim(level_name), new_grid%nCells)
 
    if (.not. allocated(new_grid%mask)) allocate(new_grid%mask(new_grid%nrows, new_grid%ncols))
    if (.not. allocated(new_grid%x)) allocate(new_grid%x(new_grid%nrows, new_grid%ncols))
    if (.not. allocated(new_grid%y)) allocate(new_grid%y(new_grid%nrows, new_grid%ncols))
    ! read L1 mask
    var = nc%getVariable("L" // trim(level_name) // "_domain_mask")
    ! read integer
    call var%getData(dummyi2)
    ! transform to logical
    new_grid%mask = (dummyi2 .eq. 1_i4)
 
    var = nc%getVariable("L" // trim(level_name) // "_domain_lat")
    call var%getData(new_grid%y)
 
    var = nc%getVariable("L" // trim(level_name) // "_domain_lon")
    call var%getData(new_grid%x)
 
    var = nc%getVariable("L" // trim(level_name) // "_domain_cellarea")
    call var%getData(dummyd2)
    if (.not. allocated(new_grid%CellArea)) new_grid%CellArea = pack(dummyd2 / 1.0e-6_dp, new_grid%mask)
    ! new_grid%CellArea = pack(dummyD2 / 1.0E-6_dp, new_grid%mask)
 
    call nc%close()
 
    new_grid%Id = (/ (k, k = 1, new_grid%nCells) /)
 
  end subroutine read_grid_info
 
 
  !> \brief read nubmer of LAI time steps and check dimension configurations read from restart file
  !> \author Robert Schweppe
  !> \date Aug 2019
  !> \author Sebastian Mueller
  !> \date Feb 2023
  subroutine read_nlai_and_check_dims(iDomain, InFile)
 
    use mo_mpr_global_variables, only: nlai, laiboundaries ! may read from restart
    use mo_netcdf, only : ncdataset, ncvariable, ncdimension
    use mo_common_constants, only : soilhorizonsvarname, landcoverperiodsvarname, laivarname
    use mo_common_mhm_mrm_variables, only: read_old_style_restart_bounds
 
    implicit none
 
    !> domain counter (not ID)
    integer(i4), intent(in) :: idomain
    !> Input Path including trailing slash
    character(256), intent(in) :: infile
 
    character(256) :: fname
    type(ncdataset) :: nc
    type(ncvariable) :: var
    type(ncdimension) :: nc_dim
 
    integer(i4) :: nsoilhorizons_temp, nlais_temp, nlandcoverperiods_temp
    real(dp), dimension(:), allocatable :: landcoverperiodboundaries_temp, soilhorizonboundaries_temp, &
            laiboundaries_temp
 
    ! dummy, 2 dimension
    real(dp), dimension(:, :), allocatable :: dummyd2, dummyd2_tmp
 
    integer(i4) :: ii
 
    ! read config
    fname = trim(infile)
    call message('    Reading and checking LAI, land-cover and soil-horizons from     ', trim(adjustl(fname)), ' ...')
 
    nc = ncdataset(fname, "r")
 
    ! get the dimensions
    var = nc%getVariable(trim(soilhorizonsvarname)//'_bnds')
    call var%getData(dummyd2_tmp)
    if (allocated(dummyd2)) deallocate(dummyd2)
    if ( read_old_style_restart_bounds ) then
      allocate(dummyd2(size(dummyd2_tmp,2), size(dummyd2_tmp,1)))
      dummyd2 = transpose(dummyd2_tmp)
    else
      allocate(dummyd2(size(dummyd2_tmp,1), size(dummyd2_tmp,2)))
      dummyd2 = dummyd2_tmp
    end if
    deallocate(dummyd2_tmp)
    nsoilhorizons_temp = size(dummyd2, 2)
    allocate(soilhorizonboundaries_temp(nsoilhorizons_temp+1))
    soilhorizonboundaries_temp(1:nsoilhorizons_temp) = dummyd2(1, :)
    soilhorizonboundaries_temp(nsoilhorizons_temp+1) = dummyd2(2, nsoilhorizons_temp)
 
    ! get the landcover dimension
    var = nc%getVariable(trim(landcoverperiodsvarname)//'_bnds')
    call var%getData(dummyd2_tmp)
    if (allocated(dummyd2)) deallocate(dummyd2)
    if ( read_old_style_restart_bounds ) then
      allocate(dummyd2(size(dummyd2_tmp,2), size(dummyd2_tmp,1)))
      dummyd2 = transpose(dummyd2_tmp)
    else
      allocate(dummyd2(size(dummyd2_tmp,1), size(dummyd2_tmp,2)))
      dummyd2 = dummyd2_tmp
    end if
    deallocate(dummyd2_tmp)
    nlandcoverperiods_temp = size(dummyd2, 2)
    allocate(landcoverperiodboundaries_temp(nlandcoverperiods_temp+1))
    landcoverperiodboundaries_temp(1:nlandcoverperiods_temp) = dummyd2(1, :)
    landcoverperiodboundaries_temp(nlandcoverperiods_temp+1) = dummyd2(2, nlandcoverperiods_temp)
 
    ! get the LAI dimension
    if (nc%hasVariable(trim(laivarname)//'_bnds')) then
      var = nc%getVariable(trim(laivarname)//'_bnds')
      call var%getData(dummyd2_tmp)
      if (allocated(dummyd2)) deallocate(dummyd2)
      if ( read_old_style_restart_bounds ) then
        allocate(dummyd2(size(dummyd2_tmp,2), size(dummyd2_tmp,1)))
        dummyd2 = transpose(dummyd2_tmp)
      else
        allocate(dummyd2(size(dummyd2_tmp,1), size(dummyd2_tmp,2)))
        dummyd2 = dummyd2_tmp
      end if
      deallocate(dummyd2_tmp)
      nlais_temp = size(dummyd2, 2)
      allocate(laiboundaries_temp(nlais_temp+1))
      laiboundaries_temp(1:nlais_temp) = dummyd2(1, :)
      laiboundaries_temp(nlais_temp+1) = dummyd2(2, nlais_temp)
    else if (nc%hasDimension('L1_LAITimesteps')) then
      nc_dim = nc%getDimension('L1_LAITimesteps')
      nlais_temp = nc_dim%getLength()
      allocate(laiboundaries_temp(nlais_temp+1))
      laiboundaries_temp = [(ii, ii=1, nlais_temp+1)]
    else
      call error_message('***ERROR: no LAI information in restart file for reading')
    end if
 
    ! check LAI for consistency on all domains (-1 indicates first reading)
    if (nlai == -1_i4) then
      nlai = nlais_temp
      allocate(laiboundaries(nlai + 1_i4))
      laiboundaries = laiboundaries_temp
    end if
 
    call check_dimension_consistency(idomain, nsoilhorizons_temp, soilhorizonboundaries_temp, &
          nlais_temp, laiboundaries_temp, nlandcoverperiods_temp, landcoverperiodboundaries_temp)
 
    call nc%close()
 
  end subroutine read_nlai_and_check_dims
 
  !> \brief checks dimension configurations read from restart file
  !> \authors Robert Schweppe
  !> \date Aug 2019
  subroutine check_dimension_consistency(iDomain, nSoilHorizons_temp, soilHorizonBoundaries_temp, &
      nLAIs_temp, LAIBoundaries_temp, nLandCoverPeriods_temp, landCoverPeriodBoundaries_temp)
 
    use mo_mpr_global_variables, only: nsoilhorizons_mhm, horizondepth_mhm, nlai, laiboundaries ! may read from restart
    use mo_common_variables, only: nlcoverscene, lc_year_start, lc_year_end ! read from nml
    use mo_string_utils, only: compress, num2str
 
    integer(i4), intent(in) :: iDomain
 
    integer(i4), intent(in) :: nSoilHorizons_temp, nLAIs_temp, nLandCoverPeriods_temp
    real(dp), dimension(:), intent(inout) :: landCoverPeriodBoundaries_temp, soilHorizonBoundaries_temp, &
          LAIBoundaries_temp
    character(256) :: errorString
 
    integer(i4) :: k
 
    ! compare local to global
    call check_consistency_element(nlcoverscene, nlandcoverperiods_temp, 'number of land cover periods', idomain)
    call check_consistency_element(nsoilhorizons_mhm, nsoilhorizons_temp, 'number of soil horizons', idomain)
    call check_consistency_element(nlai, nlais_temp, 'number of LAI timesteps', idomain)
 
    ! now check the boundaries
    do k=1, nlcoverscene
      errorstring = compress(trim(num2str(k)))//'th land cover boundary'
      call check_consistency_element(real(lc_year_start(k), dp), landcoverperiodboundaries_temp(k), errorstring, idomain)
    end do
    errorstring = 'last land cover boundary (with 1 year added due to real/int conversion) '
    call check_consistency_element(real(lc_year_end(nlcoverscene) + 1_i4, dp), &
          landcoverperiodboundaries_temp(nlcoverscene+1), errorstring, idomain)
 
    ! last soil horizon is spatially variable, so this is not checked yet
    ! first soil horizon 0 and not contained in HorizonDepth_mHM, so skip that, too
    do k=2, nsoilhorizons_mhm
      errorstring = compress(trim(num2str(k)))//'th soil horizon boundary'
      call check_consistency_element(horizondepth_mhm(k-1), soilhorizonboundaries_temp(k), errorstring, idomain)
    end do
 
    do k=1, nlai+1
      errorstring = compress(trim(num2str(k)))//'th LAI period boundary'
      call check_consistency_element(laiboundaries(k), laiboundaries_temp(k), errorstring, idomain)
    end do
 
  end subroutine check_dimension_consistency
 
  !> \copydoc check_consistency_element
  subroutine check_consistency_element_dp(item1, item2, name, iDomain)
    use mo_utils, only: ne
    use mo_string_utils, only: compress, num2str
 
    real(dp), intent(in) :: item1, item2
    character(*), intent(in) :: name
    integer(i4), intent(in) :: iDomain
 
    if (ne(item1, item2)) then
    call error_message( &
      'The ', trim(name),&
      ' as set in the configuration file (', &
      compress(trim(num2str(item1))), &
      ') does not conform with domain ', &
      compress(trim(num2str(idomain))), ' (', compress(trim(num2str(item2))), ').')
    end if
  end subroutine check_consistency_element_dp
 
  !> \copydoc check_consistency_element
  subroutine check_consistency_element_i4(item1, item2, name, iDomain)
    use mo_string_utils, only: compress, num2str
 
    integer(i4), intent(in) :: item1, item2
    character(*), intent(in) :: name
    integer(i4), intent(in) :: iDomain
 
    if (item1 /= item2) then
    call error_message( &
      'The ', trim(name),&
      ' as set in the configuration file (', &
      compress(trim(num2str(item1))), &
      ') does not conform with domain ', &
      compress(trim(num2str(idomain))), ' (', compress(trim(num2str(item2))), ').')
    end if
  end subroutine check_consistency_element_i4
 
end module mo_common_restart