mo_mpr_startup.f90

Go to the documentation of this file.

!> \file mo_mpr_startup.f90
!> \brief \copybrief mo_mpr_startup
!> \details \copydetails mo_mpr_startup
 
!> \brief Startup procedures for mHM.
!> \details This module initializes all variables required to run mHM. This
!! module needs to be run only one time at the beginning of a simulation if re-starting files do not exist.
!> \authors Luis Samaniego, Rohini Kumar
!> \date Dec 2012
!> \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_mpr
MODULE mo_mpr_startup
 
  USE mo_kind, ONLY : i4, dp
  use mo_common_constants, only : nodata_i4, nodata_dp   ! global nodata values (i4, dp)
 
  IMPLICIT NONE
 
  PRIVATE
 
  PUBLIC :: mpr_initialize, init_eff_params                      ! initialization sequence
 
CONTAINS
 
  !> \brief Initialize main mHM variables
  !> \details Initialize main mHM variables for a given domain.
  !! Calls the following procedures in this order:
  !! - Constant initialization.
  !! - Generate soil database.
  !! - Checking inconsistencies input fields.
  !! - Variable initialization at level-0.
  !! - Variable initialization at level-1.
  !! - Variable initialization at level-11.
  !! - Space allocation of remaining variable/parameters.
  !! Global variables will be used at this stage.
  !> \changelog
  !! - Luis Samaniego Mar 2008
  !!   - fully distributed multilayer
  !! - Rohini Kumar   Oct 2010
  !!   - matrix to vector version
  !!   - openmp parallelization
  !!   - routing level 11
  !! - Luis Samaniego Jul 2012
  !!   - removal of IMSL dependencies
  !! - Luis Samaniego Dec 2012
  !!   - modular version
  !! - Rohini Kumar   May 2013
  !!   - code cleaned and error checks
  !! - Rohini Kumar   Nov 2013
  !!   - updated documentation
  !! - Stephan Thober Jun 2014
  !!   - copied L2 initialization from mo_meteo_forcings
  !! - Stephan Thober Jun 2014
  !!   - updated flag for read_restart
  !! - Stephan Thober Aug 2015
  !!   - removed initialisation of routing
  !! - Rohini Kumar   Mar 2016
  !!   - changes for handling multiple soil database options
  !! - Robert Schweppe Jun 2018
  !!   - refactoring and reformatting
  !> \authors Luis Samaniego, Rohini Kumar
  !> \date Dec 2012
  subroutine mpr_initialize
 
    use mo_common_variables, only : l0_l1_remap, level0, level1, domainmeta, resolutionhydrology
    use mo_grid, only : init_lowres_level, set_domain_indices
    use mo_kind, only : i4
    use mo_read_latlon, only : read_latlon
    use mo_soil_database, only : generate_soil_database
 
    implicit none
 
    integer(i4) :: idomain
 
 
    ! soilDB common for all domains
    call generate_soil_database()
 
    allocate(level1(domainmeta%nDomains))
    allocate(l0_l1_remap(domainmeta%nDomains))
 
    ! L0 and L1 initialization
    do idomain = 1, domainmeta%nDomains
      if (idomain .eq. 1) then
        call l0_check_input(domainmeta%L0DataFrom(idomain))
        call l0_variable_init(domainmeta%L0DataFrom(idomain))
      ! ToDo: adopt to parallel
      ! ToDo: check change
      else if (domainmeta%L0DataFrom(idomain) == idomain) then
        ! this needs only be done if there is new input
        call l0_check_input(domainmeta%L0DataFrom(idomain))
        call l0_variable_init(domainmeta%L0DataFrom(idomain))
      end if
 
      call init_lowres_level(level0(domainmeta%L0DataFrom(idomain)), resolutionhydrology(idomain), &
      level1(idomain), l0_l1_remap(idomain))
      ! read lat lon coordinates for level 1
      call read_latlon(idomain, "lon", "lat", "level1", level1(idomain))
 
      ! Parameter fields have to be allocated in any case
      call init_eff_params(level1(idomain)%nCells)
 
    end do
 
    call set_domain_indices(level1)
 
  end subroutine mpr_initialize
 
 
  !> \brief Check for errors in L0 input data
  !> \details Check for possible errors in input data (morphological and land cover) at level-0
  !> \changelog
  !! - Rohini Kumar   Aug  2013
  !!   - added iFlag_LAI_data_format to handle LAI options, and changed within the code made accordingly
  !! - Rohini  Kumar  Sep 2013
  !!   - read input data for routing processes according & Stephan Thober, to process_matrix flag
  !! - Stephan Thober Aug 2015
  !!   - moved check of L0 routing variables to mRM
  !! - Rohini Kumar   Mar 2016
  !!   - changes for handling multiple soil database options
  !! - Robert Schweppe Jun 2018
  !!   - refactoring and reformatting
  !> \authors Rohini Kumar
  !> \date Jan 2013
  subroutine l0_check_input(iDomain)
 
    use mo_common_constants, only : eps_dp
    use mo_common_variables, only : l0_lcover, l0_elev, level0, nlcoverscene
    use mo_message, only : error_message
    use mo_mpr_global_variables, only : l0_asp, l0_geounit, l0_gridded_lai, &
                                        l0_slope, l0_soilid, iflag_soildb, nsoilhorizons_mhm, timestep_lai_input
    use mo_string_utils, only : num2str
    use mo_utils, only : eq
 
    implicit none
 
    !> domain id
    integer(i4), intent(in) :: iDomain
 
    integer(i4) :: k, n, nH
 
    CHARACTER(len=1024) :: message_text = ''
 
    ! START CHECKING VARIABLES
    do k = level0(idomain)%iStart, level0(idomain)%iEnd
 
      ! elevation [m]
      if (abs(l0_elev(k) - nodata_dp) .lt. eps_dp) then
        message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(idomain, '(I5)'))
        call error_message(' Error: elevation has missing value within the valid masked area at cell in domain ', &
                trim(message_text))
      end if
 
      ! slope [%]
      if (abs(l0_slope(k) - nodata_dp) .lt. eps_dp) then
        message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(idomain, '(I5)'))
        call error_message(' Error: slope has missing value within the valid masked area at cell in domain ', &
                trim(message_text))
      end if
 
      ! aspect [degree]
      if (abs(l0_asp(k) - nodata_dp) .lt. eps_dp) then
        message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(idomain, '(I5)'))
        call error_message(' Error: aspect has missing values within the valid masked area at cell in domain ', &
                trim(message_text))
      end if
 
      ! soil-Id [-]
      nh = 1 !> by default; when iFlag_soilDB = 0
      if (iflag_soildb .eq. 1) nh = nsoilhorizons_mhm
      ! another option to handle multiple soil horizons properties
      do n = 1, nh
        if (l0_soilid(k, n) .eq. nodata_i4) then
          message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(idomain, '(I5)')) // ',' // trim(num2str(n, '(I5)'))
          call error_message(' Error: soil id has missing values within the valid masked area at cell in domain and horizon ', &
                  trim(message_text))
        end if
      end do
 
      ! geological-Id [-]
      if (l0_geounit(k) .eq. nodata_i4) then
        message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(idomain, '(I5)'))
        call error_message(' Error: geological formation id has missing values within the valid masked area at cell in domain ', &
                trim(message_text))
      end if
 
      ! landcover scenes
      do  n = 1, nlcoverscene
        if (l0_lcover(k, n) .eq. nodata_i4) then
          message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(idomain, '(I5)')) // ',' // trim(num2str(n, '(I5)'))
          call error_message(' Error: land cover id has missing values within the valid masked area at cell in domain and scene ', &
                  trim(message_text))
        end if
      end do
 
      ! land cover scenes related to LAI
      if(timestep_lai_input .EQ. 0) then
        if (eq(l0_gridded_lai(k, 1), nodata_dp)) then
          message_text = trim(num2str(k, '(G5.3)')) // ',' // trim(num2str(idomain, '(I5)'))
          call error_message(' Error: gridded LAI has missing values within the valid masked area at cell in domain ', &
                  trim(message_text))
        end if
      end if
 
    end do
 
  end subroutine l0_check_input
 
 
  !> \brief level 0 variable initialization
  !> \details following tasks are performed for L0 data sets
  !! -  cell id & numbering
  !! -  storage of cell cordinates (row and coloum id)
  !! -  empirical dist. of terrain slope
  !! -  flag to determine the presence of a particular soil id
  !! in this configuration of the model run
  !! If a variable is added or removed here, then it also has to
  !! be added or removed in the subroutine config_variables_set in
  !! module mo_restart and in the subroutine set_config in module
  !! mo_set_netcdf_restart
  !> \changelog
  !! - Rohini Kumar & Matthias Cuntz  May 2014
  !!   - cell area calulation based on a regular lat-lon grid or on a regular X-Y coordinate system
  !! - Matthias Cuntz           May 2014
  !!   - changed empirical distribution function so that doubles get the same value
  !! - Matthias Zink & Matthias Cuntz Feb 2016
  !!   - code speed up due to reformulation of CDF calculation
  !! - Rohini Kumar             Mar 2016
  !!   - changes for handling multiple soil database options
  !! - Maren Kaluza             Feb 2018
  !!   - removed slope_val, temp, only sort the index to speed up finding the empirical distribution slope_emp
  !! - Robert Schweppe          Jun 2018
  !!   - refactoring and reformatting
  !> \authors Rohini Kumar
  !> \date Jan 2013
  subroutine l0_variable_init(iDomain)
 
    use mo_append, only : append
    use mo_common_variables, only : level0
    use mo_grid, only : l0_grid_setup
    use mo_mpr_global_variables, only : l0_slope, l0_slope_emp, l0_soilid, iflag_soildb, nsoilhorizons_mhm, &
            nsoiltypes, soildb
    use mo_orderpack, only : sort_index
    use mo_utils, only : eq
 
    implicit none
 
    !> domain id
    integer(i4), intent(in) :: iDomain
 
    real(dp), dimension(:), allocatable :: slope_emp
 
    integer(i4), dimension(:), allocatable :: slope_sorted_index
 
    integer(i4) :: i, j, k, nH, i_sort, i_sortpost
 
    ! STEPS ::
 
 
    !--------------------------------------------------------
    ! 1) Estimate each variable locally for a given domain
    ! 2) Pad each variable to its corresponding global one
    !--------------------------------------------------------
    !------------------------------------------------------
    ! Assign whether a given soil type is present or not
    !------------------------------------------------------
    if (idomain .eq. 1) then
      allocate(soildb%is_present(nsoiltypes))
      soildb%is_present(:) = 0_i4
    end if
 
    call l0_grid_setup(level0(idomain))
 
    !---------------------------------------------------
    ! Estimate empirical distribution of slope
    !---------------------------------------------------
    allocate(slope_emp(level0(idomain)%nCells), slope_sorted_index(level0(idomain)%nCells))
 
    ! get sorted data and sorted indexes to remap later
    slope_sorted_index = sort_index(l0_slope(level0(idomain)%iStart : level0(idomain)%iEnd))
 
    ! empirical distribution of slopes = cumulated number points with slopes that are <= the slope at this point
    !
    !       sorted data                     emp. CDF
    ! 9 |             x x       7/8 |             x x
    !   |                           |
    ! 8 |           x           5/8 |           x
    !   |                           |
    ! 5 |     x x x             4/8 |     x x x
    !   |                           |
    ! 2 |  x                    1/8 |  x
    !   |__________________         |__________________
    !
    ! highest slope value = highest rank or No. of data points / (data points + 1)
    slope_emp(slope_sorted_index(level0(idomain)%nCells)) = real(level0(idomain)%nCells, dp) / &
            real(level0(iDomain)%nCells + 1_i4, dp)
 
    ! backward loop to check if the preceding data point has the same slope value
    do i = level0(idomain)%nCells - 1, 1, -1
      i_sort=slope_sorted_index(i)
      i_sortpost=slope_sorted_index(i+1)
      if (eq(l0_slope(level0(idomain)%iStart-1_i4+i_sort), l0_slope(level0(idomain)%iStart-1_i4+i_sortpost))) then
        ! if yes: assign the same probabitity
        slope_emp(i_sort) = slope_emp(i_sortpost)
      else
        ! if not: assign rank / (data points + 1)
        slope_emp(i_sort) = real(i, dp) / real(level0(idomain)%nCells + 1_i4, dp)
      end if
    end do
 
    ! EXAMPLE
    ! in      = [  7, 20, 31, 31, 12, 31, 42 ]
    ! sorted  = [  7, 12, 20, 31, 31, 31, 42 ]
    ! index   = [  1,  5,  2,  3,  4,  6,  7 ]
    ! temp    = [  1,  2,  3,  6,  6,  6,  7 ]
    ! out     = [  1,  3,  6,  6,  2,  6,  7 ] / (len(out) + 1 )
 
    !--------------------------------------------------------
    ! Start padding up local variables to global variables
    !--------------------------------------------------------
    call append(l0_slope_emp, slope_emp)
 
    nh = 1_i4 !> by default; when iFlag_soilDB = 0
    if (iflag_soildb .eq. 1) nh = nsoilhorizons_mhm
    do i = 1, nh
      do k = level0(idomain)%iStart, level0(idomain)%iEnd
        j = l0_soilid(k, i)
        soildb%is_present(j) = 1_i4
      end do
    end do
 
    ! free space
    deallocate(slope_emp, slope_sorted_index)
 
 
  end subroutine l0_variable_init
 
 
  !> \brief Allocation of space for mHM related L1 and L11 variables.
  !> \details Allocation of space for mHM related L1 and L11 variables (e.g., states,
  !! fluxes, and parameters) for a given domain. Variables allocated here is
  !! defined in them mo_global_variables.f90 file. After allocating any variable
  !! in this routine, initalize them in the following variables_default_init
  !! subroutine.
  !> \changelog
  !! - R. Kumar           Sep 2013
  !!   - documentation added according to the template
  !! - S. Thober          Aug 2015
  !!   - removed routing related variables
  !! - Zink M. Demirel C. Mar 2017
  !!   - Init Jarvis soil water stress variable at SM process(3)
  !! - Robert Schweppe    Dec 2017
  !!   - restructured allocation in variables_alloc, expanded dimensions of effective parameters
  !! - Robert Schweppe    Jun 2018
  !!   - refactoring and reformatting
  !! - Rohini Kumar       Oct 2021
  !!   - Added Neutron count module to mHM integrate into develop branch (5.11.2)
  !! - Sebastian Müller Mar 2023
  !!   - made L1_alpha, L1_kSlowFlow, L1_kBaseFlow and L1_kPerco land cover dependent
  !> \authors Rohini Kumar
  !> \date Jan 2013
  subroutine init_eff_params(ncells1)
 
    use mo_append, only : append
    use mo_constants, only : yearmonths
    use mo_common_constants, only : p1_initstatefluxes
    use mo_common_variables, only : nlcoverscene
    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, &
                                        l1_no_count, l1_bulkdens, l1_latticewater, l1_cosmicl3
 
    implicit none
 
    !> number of L1 cells
    integer(i4), intent(in) :: ncells1
 
    real(dp), dimension(:, :, :), allocatable :: dummy_3d
 
    integer(i4) :: max_extent
 
 
    ! get maximum extent of one dimension 2 or 3
    max_extent = max(nsoilhorizons_mhm, int(yearmonths, i4), nlcoverscene, nlai)
 
    ! for appending and intialization
    allocate(dummy_3d(ncells1, max_extent, nlcoverscene))
 
    dummy_3d = p1_initstatefluxes
 
    !-------------------------------------------
    ! EFFECTIVE PARAMETERS
    !-------------------------------------------
    call append(l1_fsealed, dummy_3d(:, 1 : 1, 1 : nlcoverscene))
    ! exponent for the upper reservoir
    call append(l1_alpha, dummy_3d(:, 1 : 1, 1 : nlcoverscene))
    ! increase of the Degree-day factor per mm of increase in precipitation
    call append(l1_degdayinc, dummy_3d(:, 1 : 1, 1 : nlcoverscene))
    ! maximum degree-day factor
    call append(l1_degdaymax, dummy_3d(:, 1 : 1, 1 : nlcoverscene))
    ! degree-day factor with no precipitation
    call append(l1_degdaynopre, dummy_3d(:, 1 : 1, 1 : nlcoverscene))
    ! degree-day factor
    call append(l1_degday, dummy_3d(:, 1 : 1, 1 : nlcoverscene))
    ! Karstic percolation loss
    call append(l1_karstloss, dummy_3d(:, 1 : 1, 1 : 1))
    ! PET correction factor due to terrain aspect
    call append(l1_fasp, dummy_3d(:, 1 : 1, 1 : 1))
    ! PET correction factor due to LAI
    call append(l1_petlaicorfactor, dummy_3d(:, 1 : nlai, 1 : nlcoverscene))
    ! PET Hargreaves Samani coefficient
    call append(l1_harsamcoeff, dummy_3d(:, 1 : 1, 1 : 1))
    ! PET Prietley Taylor coefficient
    call append(l1_prietayalpha, dummy_3d(:, 1 : nlai, 1 : 1))
    ! PET aerodynamical resistance
    call append(l1_aeroresist, dummy_3d(:, 1 : nlai, 1 : nlcoverscene))
    ! PET bulk surface resistance
    call append(l1_surfresist, dummy_3d(:, 1 : nlai, 1 : 1))
    ! Fraction of roots in soil horizons
    call append(l1_froots, dummy_3d(:, 1 : nsoilhorizons_mhm, 1 : nlcoverscene))
    ! Maximum interception
    call append(l1_maxinter, dummy_3d(:, 1 : nlai, 1 : 1))
    ! fast interflow recession coefficient
    call append(l1_kfastflow, dummy_3d(:, 1 : 1, 1 : nlcoverscene))
    ! slow interflow recession coefficient
    call append(l1_kslowflow, dummy_3d(:, 1 : 1, 1 : nlcoverscene))
    ! baseflow recession coefficient
    call append(l1_kbaseflow, dummy_3d(:, 1 : 1, 1 : nlcoverscene))
    ! percolation coefficient
    call append(l1_kperco, dummy_3d(:, 1 : 1, 1 : nlcoverscene))
    ! Soil moisture below which actual ET is reduced linearly till PWP
    call append(l1_soilmoistfc, dummy_3d(:, 1 : nsoilhorizons_mhm, 1 : nlcoverscene))
    ! Saturation soil moisture for each horizon [mm]
    call append(l1_soilmoistsat, dummy_3d(:, 1 : nsoilhorizons_mhm, 1 : nlcoverscene))
    ! jarvis critical value for normalized soil water content
    call append(l1_jarvis_thresh_c1, dummy_3d(:, 1 : 1, 1 : 1))
    ! Exponential parameter to how non-linear is the soil water retention
    call append(l1_soilmoistexp, dummy_3d(:, 1 : nsoilhorizons_mhm, 1 : nlcoverscene))
    ! Threshold temperature for snow/rain
    call append(l1_tempthresh, dummy_3d(:, 1 : 1, 1 : nlcoverscene))
    ! Threshold water depth controlling fast interflow
    call append(l1_unsatthresh, dummy_3d(:, 1 : 1, 1 : 1))
    ! Threshold water depth for surface runoff in sealed surfaces
    call append(l1_sealedthresh, dummy_3d(:, 1 : 1, 1 : 1))
    ! Permanent wilting point
    call append(l1_wiltingpoint, dummy_3d(:, 1 : nsoilhorizons_mhm, 1 : nlcoverscene))
    ! neutron count related parameters
    call append(l1_no_count,     dummy_3d(:, 1:1,                 1:1))            ! N0 count
    call append(l1_bulkdens,     dummy_3d(:, 1:nsoilhorizons_mhm, 1:nlcoverscene)) ! bulk density
    call append(l1_latticewater, dummy_3d(:, 1:nsoilhorizons_mhm, 1:nlcoverscene)) ! lattice water
    call append(l1_cosmicl3,     dummy_3d(:, 1:nsoilhorizons_mhm, 1:nlcoverscene)) ! cosmic L3 parameter
 
    ! free space
    if (allocated(dummy_3d)) deallocate(dummy_3d)
 
  end subroutine init_eff_params
 
END MODULE mo_mpr_startup