mo_mpr_smhorizons.f90

Go to the documentation of this file.

!> \file mo_mpr_smhorizons.f90
!> \brief \copybrief mo_mpr_smhorizons
!> \details \copydetails mo_mpr_smhorizons
 
!> \brief setting up the soil moisture horizons
!> \details This module sets up the soil moisture horizons
!> \authors Stephan Thober, 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_smhorizons
 
  use mo_kind, only : i4, dp
  use mo_common_constants, only : nodata_dp
 
  implicit none
 
  public :: mpr_smhorizons
 
  private
 
contains
 
  ! ----------------------------------------------------------------------------
 
  !    NAME
  !        mpr_SMhorizons
 
  !    PURPOSE
  !>       \brief upscale soil moisture horizons
 
  !>       \details calculate soil properties at the level 1.
  !>       Global parameters needed (see mhm_parameter.nml):
  !>       - param(1) = rootFractionCoefficient_forest
  !>       - param(2) = rootFractionCoefficient_impervious
  !>       - param(3) = rootFractionCoefficient_pervious
  !>       - param(4) = infiltrationShapeFactor
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(:) :: param"               parameters
  !>       \param[in] "integer(i4), dimension(:, :) :: processMatrix" - matrix specifying user defined processes
  !>       \param[in] "integer(i4) :: iFlag_soil"                     - flags for handling multiple soil databases
  !>       \param[in] "integer(i4) :: nHorizons_mHM"                  - number of horizons to model
  !>       \param[in] "real(dp), dimension(:) :: HorizonDepth"        [10^-3 m] horizon depth from
  !>       surface, postive downwards
  !>       \param[in] "integer(i4), dimension(:) :: LCOVER0"          Land cover at level 0
  !>       \param[in] "integer(i4), dimension(:, :) :: soilID0"       soil ID at level 0
  !>       \param[in] "integer(i4), dimension(:) :: nHorizons"        horizons per soil type
  !>       \param[in] "integer(i4), dimension(:) :: nTillHorizons"    Number of Tillage horizons
  !>       \param[in] "real(dp), dimension(:, :, :) :: thetaS_till"   saturated water content of soil
  !>       horizons upto tillage depth,
  !>       f(OM, management)
  !>       \param[in] "real(dp), dimension(:, :, :) :: thetaFC_till"  Field capacity of tillage
  !>       layers; LUC dependent,
  !>       f(OM, management)
  !>       \param[in] "real(dp), dimension(:, :, :) :: thetaPW_till"  Permament wilting point of
  !>       tillage layers; LUC dependent,
  !>       f(OM, management)
  !>       \param[in] "real(dp), dimension(:, :) :: thetaS"           saturated water content of soil
  !>       horizons after tillage depth
  !>       \param[in] "real(dp), dimension(:, :) :: thetaFC"          Field capacity of deeper layers
  !>       \param[in] "real(dp), dimension(:, :) :: thetaPW"          Permanent wilting point of
  !>       deeper layers
  !>       \param[in] "real(dp), dimension(:, :, :) :: Wd"            weights of mHM Horizons
  !>       according to horizons provided
  !>       in soil database
  !>       \param[in] "real(dp), dimension(:, :, :) :: Db"            Bulk density
  !>       \param[in] "real(dp), dimension(:, :) :: DbM"              mineral Bulk density
  !>       \param[in] "real(dp), dimension(:) :: RZdepth"             [mm]       Total soil depth
  !>       \param[in] "logical, dimension(:, :) :: mask0"             mask at L0
  !>       \param[in] "integer(i4), dimension(:) :: cell_id0"         Cell ids of hi res field
  !>       \param[in] "integer(i4), dimension(:) :: upp_row_L1"       Upper row of hi res block
  !>       \param[in] "integer(i4), dimension(:) :: low_row_L1"       Lower row of hi res block
  !>       \param[in] "integer(i4), dimension(:) :: lef_col_L1"       Left column of hi res block
  !>       \param[in] "integer(i4), dimension(:) :: rig_col_L1"       Right column of hi res block
  !>       \param[in] "integer(i4), dimension(:) :: nL0_in_L1"        Number of L0 cells within a L1 cel
 
  !    INTENT(INOUT)
  !>       \param[inout] "real(dp), dimension(:, :) :: L1_beta"   Parameter that determines the
  !>       relative contribution to SM, upscaled
  !>       Bulk density
  !>       \param[inout] "real(dp), dimension(:, :) :: L1_SMs"    [10^-3 m] depth of saturated SM cont
  !>       \param[inout] "real(dp), dimension(:, :) :: L1_FC"     [10^-3 m] field capacity
  !>       \param[inout] "real(dp), dimension(:, :) :: L1_PW"     [10^-3 m] permanent wilting point
  !>       \param[inout] "real(dp), dimension(:, :) :: L1_fRoots" fraction of roots in soil horizons
 
  !    HISTORY
  !>       \authors Luis Samaniego, Rohini Kumar, Stephan Thober
 
  !>       \date Dec 2012
 
  ! Modifications:
  ! Stephan Thober                  Jan 2013 - updated calling sequence for upscaling operators
  ! Juliane Mai                     Oct 2013 - OLD parametrization
  !                                                --> param(1) = rootFractionCoefficient_forest
  !                                                --> param(2) = rootFractionCoefficient_impervious
  !                                                --> param(3) = rootFractionCoefficient_pervious
  !                                                --> param(4) = infiltrationShapeFactor
  !                                             -------------------------------
  !                                             rootFractionCoeff_perv   = rootFractionCoeff_forest - delta_1
  !                                             -------------------------------
  !                                             NEW parametrization
  !                                                --> param(1) = rootFractionCoefficient_forest
  !                                                --> param(2) = rootFractionCoefficient_impervious
  !                                                --> param(3) = delta_1
  !                                                --> param(4) = infiltrationShapeFactor
  !                                                ! if processMatrix(3,1) = 3 and 4 additionally
  !                                                --> param(5) = rootFractionCoefficient_sand
  !                                                --> param(6) = rootFractionCoefficient_clay
  !                                                --> param(7) = FCmin_glob
  !                                                --> param(8) = FCdelta_glob
  ! Stephan Thober                  Mar 2014 - added omp parallelization
  ! Rohini Kumar                    Mar 2016 - changes for handling multiple soil database options
  ! M. Cuneyd Demirel, Simon Stisen Apr 2017 - added FC dependency on root fraction coefficient
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
  ! M. Cuneyd Demirel, Simon Stisen Jun 2020 - added Feddes and FC dependency on root fraction coefficient processCase(3) = 4
  ! M. Cuneyd Demirel, Simon Stisen Feb 2021 - Bug fix normalization of FCnorm
 
  subroutine mpr_smhorizons(param, processMatrix, iFlag_soil, nHorizons_mHM, HorizonDepth, LCOVER0, soilID0, nHorizons, &
                           nTillHorizons, thetaS_till, thetaFC_till, thetaPW_till, thetaS, thetaFC, thetaPW, Wd, Db, &
                           DbM, RZdepth, mask0, cell_id0, upp_row_L1, low_row_L1, lef_col_L1, rig_col_L1, nL0_in_L1, &
                           L1_beta, L1_SMs, L1_FC, L1_PW, L1_fRoots, &
                           ! neutron count
                           latWat_till   , & ! lattice water upto tillage depth
                           COSMIC_L3_till, & ! COSMIC parameter L3 upto tillage depth
                           latWat        , & ! lattice water
                           COSMIC_L3     , & ! COSMIC paramter L3
                           L1_bulkDens   , & ! L bulk density
                           L1_latticeWater,& ! L1 lattice water content
                           L1_COSMICL3   )   ! L1 COSMIC L3 parameter from neutron module
 
    use mo_message, only : message, error_message
    use mo_string_utils, only : num2str
    use mo_upscaling_operators, only : upscale_harmonic_mean
    !$ use omp_lib
 
    implicit none
 
    ! parameters
    real(dp), dimension(:), intent(in) :: param
 
    ! - matrix specifying user defined processes
    integer(i4), dimension(:, :), intent(in) :: processmatrix
 
    ! - flags for handling multiple soil databases
    integer(i4), intent(in) :: iflag_soil
 
    ! - number of horizons to model
    integer(i4), intent(in) :: nhorizons_mhm
 
    ! [10^-3 m] horizon depth from
    ! surface, postive downwards
    real(dp), dimension(:), intent(in) :: horizondepth
 
    ! Land cover at level 0
    integer(i4), dimension(:), intent(in) :: lcover0
 
    ! soil ID at level 0
    integer(i4), dimension(:, :), intent(in) :: soilid0
 
    ! horizons per soil type
    integer(i4), dimension(:), intent(in) :: nhorizons
 
    ! Number of Tillage horizons
    integer(i4), dimension(:), intent(in) :: ntillhorizons
 
    ! saturated water content of soil
    ! horizons upto tillage depth,
    ! f(OM, management)
    real(dp), dimension(:, :, :), intent(in) :: thetas_till
 
    ! Field capacity of tillage
    ! layers; LUC dependent,
    ! f(OM, management)
    real(dp), dimension(:, :, :), intent(in) :: thetafc_till
 
    ! Permament wilting point of
    ! tillage layers; LUC dependent,
    ! f(OM, management)
    real(dp), dimension(:, :, :), intent(in) :: thetapw_till
 
    ! saturated water content of soil
    ! horizons after tillage depth
    real(dp), dimension(:, :), intent(in) :: thetas
 
    ! Field capacity of deeper layers
    real(dp), dimension(:, :), intent(in) :: thetafc
 
    ! Permanent wilting point of
    ! deeper layers
    real(dp), dimension(:, :), intent(in) :: thetapw
 
    ! weights of mHM Horizons
    ! according to horizons provided
    ! in soil database
    real(dp), dimension(:, :, :), intent(in) :: wd
 
    ! Bulk density
    real(dp), dimension(:, :, :), intent(in) :: db
 
    ! mineral Bulk density
    real(dp), dimension(:, :), intent(in) :: dbm
 
    ! [mm]       Total soil depth
    real(dp), dimension(:), intent(in) :: rzdepth
 
    ! mask at L0
    logical, dimension(:, :), intent(in) :: mask0
 
    ! Cell ids of hi res field
    integer(i4), dimension(:), intent(in) :: cell_id0
 
    ! Upper row of hi res block
    integer(i4), dimension(:), intent(in) :: upp_row_l1
 
    ! Lower row of hi res block
    integer(i4), dimension(:), intent(in) :: low_row_l1
 
    ! Left column of hi res block
    integer(i4), dimension(:), intent(in) :: lef_col_l1
 
    ! Right column of hi res block
    integer(i4), dimension(:), intent(in) :: rig_col_l1
 
    ! Number of L0 cells within a L1 cel
    integer(i4), dimension(:), intent(in) :: nl0_in_l1
 
    ! Parameter that determines the
    ! relative contribution to SM, upscaled
    ! Bulk density
    real(dp), dimension(:, :), intent(inout) :: l1_beta
 
    ! [10^-3 m] depth of saturated SM cont
    real(dp), dimension(:, :), intent(inout) :: l1_sms
 
    ! [10^-3 m] field capacity
    real(dp), dimension(:, :), intent(inout) :: l1_fc
 
    ! [10^-3 m] permanent wilting point
    real(dp), dimension(:, :), intent(inout) :: l1_pw
 
    ! fraction of roots in soil horizons
    real(dp), dimension(:, :), intent(inout) :: l1_froots
 
    ! neutron count
    real(dp), dimension(:,:,:), intent(in) :: latwat_till   ! lattice water
    real(dp), dimension(:,:,:), intent(in) :: cosmic_l3_till! COSMIC parameter L3
    real(dp), dimension(:,:),   intent(in) :: latwat        ! lattice water
    real(dp), dimension(:,:),   intent(in) :: cosmic_l3     ! COSMIC paramter L3
    ! out
    real(dp), dimension(:,:), intent(inout) :: l1_bulkdens
    real(dp), dimension(:,:), intent(inout) :: l1_latticewater
    real(dp), dimension(:,:), intent(inout) :: l1_cosmicl3
 
 
    ! loop index
    integer(i4) :: h
 
    ! loop index
    integer(i4) :: k
 
    ! loop index
    integer(i4) :: l
 
    ! loop index
    integer(i4) :: s
 
    real(dp) :: dpth_f
 
    real(dp) :: dpth_t
 
    real(dp) :: ftotroots
 
    ! beta 0
    real(dp), dimension(size(LCOVER0, 1)) :: beta0
 
    ! [10^3 kg/m3] Bulk density
    real(dp), dimension(size(LCOVER0, 1)) :: bd0
 
    ! [10^-3 m] depth of saturated SM cont
    real(dp), dimension(size(LCOVER0, 1)) :: sms0
 
    ! [10^-3 m] field capacity
    real(dp), dimension(size(LCOVER0, 1)) :: fc0
 
    ! [10^-3 m] permanent wilting point
    real(dp), dimension(size(LCOVER0, 1)) :: pw0
 
    ! neutron count
    real(dp), dimension(size(LCOVER0,1))    :: lw0     ! lattice water
    real(dp), dimension(size(LCOVER0,1))    :: l30     ! COSMIC parameter L3
 
 
    ! fraction of roots in soil horizons
    real(dp), dimension(size(LCOVER0, 1)) :: froots0
 
    real(dp) :: tmp_rootfractioncoefficient_forest
 
    real(dp) :: tmp_rootfractioncoefficient_impervious
 
    real(dp) :: tmp_rootfractioncoefficient_pervious
 
    ! Field capacity dependent
    ! root frac coeffiecient
    real(dp) :: tmp_rootfractioncoefficient_perviousfc
 
    ! Model parameter describing the threshold for
    ! actual ET reduction for sand
    real(dp) :: tmp_rootfractioncoefficient_sand
 
    ! Model parameter describing the threshold for actual
    ! ET reduction for clay
    real(dp) :: tmp_rootfractioncoefficient_clay
 
    real(dp) :: fcmin_glob
 
    ! real(dp) :: FCdelta_glob
 
    real(dp) :: fcmax_glob
 
    real(dp) :: fcnorm
    ! the minimum number of till horizons
    integer(i4) :: min_nth
 
 
    min_nth = minval(ntillhorizons(:))
    tmp_rootfractioncoefficient_forest     = param(1)     ! min(1.0_dp, param(2) + param(3) + param(1))
    tmp_rootfractioncoefficient_impervious = param(2)
 
    ! decide which parameterization should be used for route fraction:
    select case (processmatrix(3, 1))
    case(1,2)
       ! 1 and 2 - dependent on land cover
       tmp_rootfractioncoefficient_pervious = param(1) - param(3) ! min(1.0_dp, param(2) + param(3))
       !write(*,*) 'tmp_rootFractionCoefficient_forest     = ', tmp_rootFractionCoefficient_forest
       !write(*,*) 'tmp_rootFractionCoefficient_impervious = ', tmp_rootFractionCoefficient_impervious
       !write(*,*) 'tmp_rootFractionCoefficient_pervious   = ', tmp_rootFractionCoefficient_pervious
    case(3,4)
       ! 3 and 4 - dependent on land cover and additionally soil texture
       tmp_rootfractioncoefficient_pervious = param(3) ! min(1.0_dp, param(2) + param(3))
       !delta approach is used as in tmp_rootFractionCoefficient_pervious
       tmp_rootfractioncoefficient_sand = param(6) - param(5)
       !the value in parameter namelist is before substraction i.e. param(5)
       tmp_rootfractioncoefficient_clay = param(6)
       fcmin_glob=param(7)
       fcmax_glob=param(7)+param(8)
       !write(*,*) 'tmp_rootFractionCoefficient_forest     = ', tmp_rootFractionCoefficient_forest
       !write(*,*) 'tmp_rootFractionCoefficient_impervious = ', tmp_rootFractionCoefficient_impervious
       !write(*,*) 'tmp_rootFractionCoefficient_pervious   = ', tmp_rootFractionCoefficient_pervious
       !write(*,*) 'tmp_rootFractionCoefficient_sand       = ', tmp_rootFractionCoefficient_sand
       !write(*,*) 'tmp_rootFractionCoefficient_clay       = ', tmp_rootFractionCoefficient_clay
       !write(*,*) 'FCmin_glob = ', FCmin_glob
       !write(*,*) 'FCmax_glob = ', FCmax_glob
    end select
 
 
    ! select case according to a given soil database flag
    SELECT CASE(iflag_soil)
      ! classical mHM soil database format
    CASE(0)
       do h = 1, nhorizons_mhm
 
          bd0  = nodata_dp
          sms0 = nodata_dp
          fc0  = nodata_dp
          pw0  = nodata_dp
          froots0 = nodata_dp
          tmp_rootfractioncoefficient_perviousfc = nodata_dp
 
          ! neutron count
          lw0 = nodata_dp
          l30 = nodata_dp
 
          ! Initalise mHM horizon depth
          ! Last layer depth is soil type dependent, and hence it assigned within the inner loop
          ! by default for the first soil layer
          dpth_f = 0.0_dp
          dpth_t = horizondepth(h)
          ! check for the layer (2, ... n-1 layers) update depth
          if( (h .gt. 1) .and. (h .lt. nhorizons_mhm) ) then
             dpth_f = horizondepth(h-1)
             dpth_t = horizondepth(h)
          end if
 
          !$OMP PARALLEL
          !$OMP DO PRIVATE( l, s ) SCHEDULE( STATIC )
          cellloop0 : do k = 1, size(lcover0, 1)
             l = lcover0(k)
             s = soilid0(k, 1)  ! >> in this case the second dimension of soilId0 = 1
             ! depth weightage bulk density
             bd0(k) = sum(db(s, : ntillhorizons(s), l) * wd(s, h, 1 : ntillhorizons(s)), &
                  wd(s, h, 1 : ntillhorizons(s)) > 0.0_dp) &
                  + sum(dbm(s, ntillhorizons(s) + 1 : nhorizons(s)) &
                  * wd(s, h, ntillhorizons(s) + 1 : nhorizons(s)), &
                  wd(s, h, ntillhorizons(s) + 1 : nhorizons(s)) >= 0.0_dp)
             ! depth weightage thetaS
             sms0(k) = sum(thetas_till(s, : ntillhorizons(s), l) &
                  * wd(s, h, 1 : ntillhorizons(s)), &
                  wd(s, h, 1 : ntillhorizons(s)) > 0.0_dp) &
                  + sum(thetas(s, ntillhorizons(s) + 1 - min_nth : nhorizons(s) - min_nth) &
                  * wd(s, h, ntillhorizons(s) + 1 : nhorizons(s)), &
                  wd(s, h, ntillhorizons(s) + 1 : nhorizons(s)) > 0.0_dp)
             ! depth weightage FC
             fc0(k) = sum(thetafc_till(s, : ntillhorizons(s), l) &
                  * wd(s, h, 1 : ntillhorizons(s)), &
                  wd(s, h, 1 : ntillhorizons(s)) > 0.0_dp) &
                  + sum(thetafc(s, ntillhorizons(s) + 1 - min_nth : nhorizons(s) - min_nth) &
                  * wd(s, h, ntillhorizons(s) + 1 : nhorizons(s)), &
                  wd(s, h, ntillhorizons(s) + 1 : nhorizons(s)) > 0.0_dp)
             ! depth weightage PWP
             pw0(k) = sum(thetapw_till(s, : ntillhorizons(s), l) &
                  * wd(s, h, 1 : ntillhorizons(s)), &
                  wd(s, h, 1 : ntillhorizons(s)) > 0.0_dp) &
                  + sum(thetapw(s, ntillhorizons(s) + 1 - min_nth : nhorizons(s) - min_nth) &
                  * wd(s, h, ntillhorizons(s) + 1 : nhorizons(s)), &
                  wd(s, h, ntillhorizons(s) + 1 : nhorizons(s)) > 0.0_dp)
 
             ! neutron count --> depth weightage LW and L30
             lw0(k) = sum( latwat_till(s, : ntillhorizons(s), l) &
                  * wd(s, h, 1 : ntillhorizons(s) ), &
                  wd(s, h, 1 : ntillhorizons(s)) > 0.0_dp ) &
                  + sum( latwat(s,ntillhorizons(s) + 1 - min_nth : nhorizons(s) - min_nth) &
                  * wd(s, h, ntillhorizons(s) + 1 : nhorizons(s)), &
                  wd(s, h, ntillhorizons(s) + 1 : nhorizons(s)) > 0.0_dp )
 
             l30(k) = sum( cosmic_l3_till(s, : ntillhorizons(s), l) &
                  * wd(s, h, 1 : ntillhorizons(s) ), &
                  wd(s, h, 1 : ntillhorizons(s)) > 0.0_dp ) &
                  + sum( cosmic_l3(s,ntillhorizons(s) + 1 - min_nth : nhorizons(s) - min_nth) &
                  * wd(s, h, ntillhorizons(s) + 1 : nhorizons(s)), &
                  wd(s, h, ntillhorizons(s) + 1 : nhorizons(s)) > 0.0_dp )
 
 
             ! Horizon depths: last soil horizon is varying, and thus the depth
             ! of the horizon too...
             if(h .eq. nhorizons_mhm) then
                dpth_f = horizondepth(nhorizons_mhm - 1)
                dpth_t = rzdepth(s)
             end if
             ! other soil properties [SMs, FC, PWP in mm]
             sms0(k) = sms0(k) * (dpth_t - dpth_f)
             fc0(k)  = fc0(k)  * (dpth_t - dpth_f)
             pw0(k)  = pw0(k)  * (dpth_t - dpth_f)
             lw0(k)  = lw0(k)  * (dpth_t - dpth_f)
          end do cellloop0
          !$OMP END DO
          !$OMP END PARALLEL
 
 
          !$OMP PARALLEL
          !$OMP DO PRIVATE( l, tmp_rootFractionCoefficient_perviousFC, FCnorm ) SCHEDULE( STATIC )
          celllloop0 : do k = 1, size(lcover0, 1)
             l = lcover0(k)
             !---------------------------------------------------------------------
             ! Effective root fractions in soil horizon...
             !  as weightage sum (according to LC fraction)
             !---------------------------------------------------------------------
             ! vertical root distribution = f(LC), following asymptotic equation
             ! [for refrence see, Jackson et. al., Oecologia, 1996. 108(389-411)]
 
             ! Roots(H) = 1 - beta^d
             !  where,
             !   Roots(H) = cumulative root fraction [-], range: 0-1
             !   beta     = fitted extinction cofficient parameter [-], as a f(LC)
             !   d        = soil surface to depth [cm]
 
             ! NOTES **
             !  sum(fRoots) for soil horions = 1.0
 
             !  if [sum(fRoots) for soil horions < 1.0], then
             !  normalise fRoot distribution such that all roots end up
             !  in soil horizon and thus satisfying the constrain that
             !  sum(fRoots) = 1
 
             !  The above constrains means that there are not roots below the soil horizon.
             !  This may or may not be realistic but it has been coded here to satisfy the
             !  conditions of the EVT vales, otherwise which the EVT values would be lesser
             !  than the acutal EVT from whole soil layers.
 
             !  Code could be modified in a way that a portion of EVT comes from the soil layer
             !  which is in between unsaturated and saturated zone or if necessary the saturated
             !  layer (i.e. Groundwater layer) can also contribute to EVT. Note that the above
             !  modification should be done only if and only if [sum(fRoots) for soil horions < 1.0].
             !  In such cases, you have to judiciously decide which layers (either soil layer between
             !  unsaturated and saturated zone or saturated zone) will contribute to EVT and in which
             !  proportions. Also note that there are no obervations on the depth avialable ata a
             !  moment on these layers.
             !------------------------------------------------------------------------
             select case(l)
             case(1)
                ! forest
                froots0(k) = (1.0_dp - tmp_rootfractioncoefficient_forest**(dpth_t * 0.1_dp)) &
                     - (1.0_dp - tmp_rootfractioncoefficient_forest**(dpth_f * 0.1_dp))
             case(2)
                ! impervious
                froots0(k) = (1.0_dp - tmp_rootfractioncoefficient_impervious**(dpth_t * 0.1_dp)) &
                     - (1.0_dp - tmp_rootfractioncoefficient_impervious**(dpth_f * 0.1_dp))
             case(3)
                ! permeable
                select case (processmatrix(3, 1))
                case(1,2)
                   ! permeable
                   froots0(k) = (1.0_dp - tmp_rootfractioncoefficient_pervious**(dpth_t * 0.1_dp)) &
                        - (1.0_dp - tmp_rootfractioncoefficient_pervious**(dpth_f * 0.1_dp))
                case(3,4)
                   ! permeable
                   ! introducing global FC dependency on root frac. coef. by Simon Stisen and M. Cuneyd Demirel from GEUS.dk
                   ! The normalization is based on Demirel et al 2018 (doi: 10.5194/hess-22-1299-2018)
                   ! Case 3 is based on Jarvis (doi: 10.1016/0022-1694(89)90050-4)
                   ! Case 4 is based on Feddes (doi: 10.1016/0022-1694(76)90017-2)
 
                   fcnorm = (((fc0(k) / (dpth_t - dpth_f)) - fcmin_glob) / (fcmax_glob - fcmin_glob))
 
                   if(fcnorm .lt. 0.0_dp) then
                      ! print*, "FCnorm is below 0, will become 0", FCnorm
                      fcnorm=0.0_dp
                   else if(fcnorm .gt. 1.0_dp) then
                      ! print*, "FCnorm is above 1, will become 1", FCnorm
                      fcnorm=1.0_dp
                   end if
 
                   tmp_rootfractioncoefficient_perviousfc = (fcnorm * tmp_rootfractioncoefficient_clay) &
                        + ((1 - fcnorm) * tmp_rootfractioncoefficient_sand)
 
                   froots0(k) = (1.0_dp - tmp_rootfractioncoefficient_perviousfc**(dpth_t * 0.1_dp)) &
                        - (1.0_dp - tmp_rootfractioncoefficient_perviousfc**(dpth_f * 0.1_dp))
 
                end select
 
                if((froots0(k) .lt. 0.0_dp) .OR. (froots0(k) .gt. 1.0_dp)) then
                  ! why is this not stopping here?
                  call message('***ERROR: Fraction of roots out of range [0,1]. Cell', &
                        num2str(k), ' has value ', num2str(froots0(k)))
                  ! stop
                end if
             end select
 
          end do celllloop0
          !$OMP END DO
          !$OMP END PARALLEL
 
          beta0 = bd0 * param(4)
 
          !---------------------------------------------
          ! Upscale the soil related parameters
          !---------------------------------------------
          l1_sms(:, h) = upscale_harmonic_mean(nl0_in_l1, upp_row_l1, low_row_l1, &
               lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, sms0)
          l1_beta(:, h) = upscale_harmonic_mean(nl0_in_l1, upp_row_l1, low_row_l1, &
               lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, beta0)
          l1_pw(:, h) = upscale_harmonic_mean(nl0_in_l1, upp_row_l1, low_row_l1, &
               lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, pw0)
          l1_fc(:, h) = upscale_harmonic_mean(nl0_in_l1, upp_row_l1, low_row_l1, &
               lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, fc0)
          l1_froots(:, h) = upscale_harmonic_mean(nl0_in_l1, upp_row_l1, low_row_l1, &
               lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, froots0)
 
          ! !> neutron count
          l1_bulkdens(:,h) = upscale_harmonic_mean( nl0_in_l1, upp_row_l1, low_row_l1, &
               lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, bd0 )
          l1_latticewater(:,h) = upscale_harmonic_mean( nl0_in_l1, upp_row_l1, low_row_l1, &
               lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, lw0 )
          l1_cosmicl3(:,h) = upscale_harmonic_mean( nl0_in_l1, upp_row_l1, low_row_l1, &
               lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, l30 )
 
       end do
 
      ! to handle multiple soil horizons with unique soil class
    CASE(1)
      ! horizon wise calculation
      do h = 1, nhorizons_mhm
        bd0 = nodata_dp
        sms0 = nodata_dp
        fc0 = nodata_dp
        pw0 = nodata_dp
        froots0 = nodata_dp
        tmp_rootfractioncoefficient_perviousfc = nodata_dp
 
        ! neutron count
        lw0     = nodata_dp
        l30     = nodata_dp
 
        ! initalise mHM horizon depth
        if (h .eq. 1) then
          dpth_f = 0.0_dp
          dpth_t = horizondepth(h)
          ! check for the layer (2, ... n-1 layers) update depth
        else
          dpth_f = horizondepth(h - 1)
          dpth_t = horizondepth(h)
        end if
        ! need to be done for every layer to get fRoots
        !$OMP PARALLEL
        !$OMP DO PRIVATE( l, s ) SCHEDULE( STATIC )
        cellloop1 : do k = 1, size(lcover0, 1)
          l = lcover0(k)
          s = soilid0(k, h)
          if (h .le. ntillhorizons(1)) then
            bd0(k) = db(s, 1, l)
            sms0(k)= thetas_till(s, 1, l) * (dpth_t - dpth_f)  ! in mm
            fc0(k) = thetafc_till(s, 1, l) * (dpth_t - dpth_f)  ! in mm
            pw0(k) = thetapw_till(s, 1, l) * (dpth_t - dpth_f)  ! in mm
            lw0(k) = latwat_till(s, 1, l)  * (dpth_t - dpth_f)  ! in mm  ! >> neutron count
          else
            bd0(k) = dbm(s, 1)
            sms0(k)= thetas(s, 1) * (dpth_t - dpth_f)  ! in mm
            fc0(k) = thetafc(s, 1) * (dpth_t - dpth_f)  ! in mm
            pw0(k) = thetapw(s, 1) * (dpth_t - dpth_f)  ! in mm
            lw0(k) = latwat(s, 1)  * (dpth_t - dpth_f)  ! in mm  ! >> neutron count
          end if
        end do cellloop1
        !$OMP END DO
        !$OMP END PARALLEL
 
 
        !$OMP PARALLEL
        !$OMP DO PRIVATE( l, tmp_rootFractionCoefficient_perviousFC, FCnorm ) SCHEDULE( STATIC )
 
        celllloop1 : do k = 1, size(lcover0, 1)
          l = lcover0(k)
          !================================================================================
          ! fRoots = f[LC] --> (fRoots(H) = 1 - beta^d)
          ! see below for comments and references for the use of this simple equation
          ! NOTE that in this equation the unit of soil depth is in cm
          !================================================================================
 
          select case(l)
          case(1)
            ! forest
            froots0(k) = (1.0_dp - tmp_rootfractioncoefficient_forest**(dpth_t * 0.1_dp)) &
                    - (1.0_dp - tmp_rootfractioncoefficient_forest**(dpth_f * 0.1_dp))
          case(2)
            ! impervious
            froots0(k) = (1.0_dp - tmp_rootfractioncoefficient_impervious**(dpth_t * 0.1_dp)) &
                    - (1.0_dp - tmp_rootfractioncoefficient_impervious**(dpth_f * 0.1_dp))
          case(3)
 
            select case (processmatrix(3, 1))
            case(1,2)
              ! permeable
              froots0(k) = (1.0_dp - tmp_rootfractioncoefficient_pervious**(dpth_t * 0.1_dp)) &
                      - (1.0_dp - tmp_rootfractioncoefficient_pervious**(dpth_f * 0.1_dp))
            case(3,4)
              ! permeable
              ! introducing global FC dependency on root frac. coef. by Simon Stisen and M. Cuneyd Demirel from GEUS.dk
              ! The normalization is based on Demirel et al 2018 (doi: 10.5194/hess-22-1299-2018)
              ! Case 3 is based on Jarvis (doi: 10.1016/0022-1694(89)90050-4)
              ! Case 4 is based on Feddes (doi: 10.1016/0022-1694(76)90017-2)
              fcnorm = (((fc0(k) / (dpth_t - dpth_f)) - fcmin_glob) / (fcmax_glob - fcmin_glob))
 
              if(fcnorm .lt. 0.0_dp) then
                ! print*, "FCnorm is below 0, will become 0", FCnorm
                fcnorm=0.0_dp
              else if(fcnorm .gt. 1.0_dp) then
                ! print*, "FCnorm is above 1, will become 1", FCnorm
                fcnorm=1.0_dp
              end if
 
              tmp_rootfractioncoefficient_perviousfc = (fcnorm * tmp_rootfractioncoefficient_clay) &
                      + ((1 - fcnorm) * tmp_rootfractioncoefficient_sand)
 
 
              froots0(k) = (1.0_dp - tmp_rootfractioncoefficient_perviousfc**(dpth_t * 0.1_dp)) &
                      - (1.0_dp - tmp_rootfractioncoefficient_perviousfc**(dpth_f * 0.1_dp))
 
            end select
 
            if((froots0(k) .lt. 0.0_dp) .OR. (froots0(k) .gt. 1.0_dp)) then
              ! why is this not stopping here?
              call message('***ERROR: Fraction of roots out of range [0,1]. Cell', &
                    num2str(k), ' has value ', num2str(froots0(k)))
              ! stop
            end if
          end select
 
        end do celllloop1
        !$OMP END DO
        !$OMP END PARALLEL
 
        ! beta parameter
        beta0 = bd0 * param(4)
 
        ! Upscale the soil related parameters
        l1_sms(:, h) = upscale_harmonic_mean(nl0_in_l1, upp_row_l1, low_row_l1, &
                lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, sms0)
        l1_beta(:, h) = upscale_harmonic_mean(nl0_in_l1, upp_row_l1, low_row_l1, &
                lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, beta0)
        l1_pw(:, h) = upscale_harmonic_mean(nl0_in_l1, upp_row_l1, low_row_l1, &
                lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, pw0)
        l1_fc(:, h) = upscale_harmonic_mean(nl0_in_l1, upp_row_l1, low_row_l1, &
                lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, fc0)
        l1_froots(:, h) = upscale_harmonic_mean(nl0_in_l1, upp_row_l1, low_row_l1, &
             lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, froots0)
 
        ! neutron count
        l1_bulkdens(:,h) = upscale_harmonic_mean( nl0_in_l1, upp_row_l1, low_row_l1, &
             lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, bd0 )
       l1_latticewater(:,h)   = upscale_harmonic_mean( nl0_in_l1, upp_row_l1, low_row_l1, &
             lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, lw0  )
        l1_cosmicl3(:,h)   = upscale_harmonic_mean( nl0_in_l1, upp_row_l1, low_row_l1, &
             lef_col_l1, rig_col_l1, cell_id0, mask0, nodata_dp, l30  )
 
      end do
      ! anything else
    CASE DEFAULT
      call error_message('***ERROR: iFlag_soilDB option given does not exist. Only 0 and 1 is taken at the moment.')
    END SELECT
 
 
    ! below operations are common to all soil databases flags
    !$OMP PARALLEL
    !------------------------------------------------------------------------
    ! CHECK LIMITS OF PARAMETERS
    !   [PW <= FC <= ThetaS]
    !  units of all variables are now in [mm]
    ! If by any means voilation of this rule appears (e.g. numerical errors)
    ! than correct it --> threshold limit = 1% of the upper ones
    !------------------------------------------------------------------------
    l1_fc = merge(l1_sms - 0.01_dp * l1_sms, l1_fc, l1_fc .gt. l1_sms)
    l1_pw = merge(l1_fc - 0.01_dp  * l1_fc,  l1_pw, l1_pw .gt. l1_fc )
    ! check the physical limit
    l1_sms = merge(0.0001_dp, l1_sms, l1_sms .lt. 0.0_dp)
    l1_fc  = merge(0.0001_dp, l1_fc,  l1_fc  .lt. 0.0_dp)
    l1_pw  = merge(0.0001_dp, l1_pw,  l1_pw  .lt. 0.0_dp)
    ! Normalise the vertical root distribution profile such that [sum(fRoots) = 1.0]
    !$OMP DO PRIVATE( fTotRoots ) SCHEDULE( STATIC )
    do k = 1, size(l1_froots, 1)
      ftotroots = sum(l1_froots(k, :), l1_froots(k, :) .gt. 0.0_dp)
      ! This if clause is necessary for test program but may be redundant in actual program
      if (ftotroots .gt. 0.0_dp) then
        l1_froots(k, :) = l1_froots(k, :) / ftotroots
      else
        l1_froots(k, :) = 0.0_dp
      end If
    end do
 
    !$OMP END DO
    !$OMP END PARALLEL
!close(1)
  end subroutine mpr_smhorizons
 
end module mo_mpr_smhorizons