mo_mpr_soilmoist.f90

Go to the documentation of this file.

!> \file mo_mpr_soilmoist.f90
!> \brief \copybrief mo_mpr_soilmoist
!> \details \copydetails mo_mpr_soilmoist
 
!> \brief Multiscale parameter regionalization (MPR) for soil moisture
!> \details This module contains all routines required for parametrizing soil moisture processes.
!> \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_soilmoist
 
  use mo_kind, only : i4, dp
  use mo_message, only : message, error_message
 
  implicit none
 
  public :: mpr_sm
 
  private
 
contains
  ! ----------------------------------------------------------------------------
 
  !    NAME
  !        mpr_sm
 
  !    PURPOSE
  !>       \brief multiscale parameter regionalization for soil moisture
 
  !>       \details This subroutine is a wrapper around all soil moisture
  !>       parameter routines. This subroutine requires 13 parameters. These
  !>       parameters have to correspond to the parameters in the original
  !>       parameter array at the following locations: 10-12, 13-18, 27-30.
  !>       Global parameters needed (see mhm_parameter.nml):
  !>       - param( 1) = orgMatterContent_forest
  !>       - param( 2) = orgMatterContent_impervious
  !>       - param( 3) = orgMatterContent_pervious
  !>       - param( 4) = PTF_lower66_5_constant
  !>       - param( 5) = PTF_lower66_5_clay
  !>       - param( 6) = PTF_lower66_5_Db
  !>       - param( 7) = PTF_higher66_5_constant
  !>       - param( 8) = PTF_higher66_5_clay
  !>       - param( 9) = PTF_higher66_5_Db
  !>       - param(10) = PTF_Ks_constant
  !>       - param(11) = PTF_Ks_sand
  !>       - param(12) = PTF_Ks_clay
  !>       - param(13) = PTF_Ks_curveSlope
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(13) :: param"           global parameters
  !>       \param[in] "integer(i4), dimension(:) :: is_present"    indicates whether soiltype is present
  !>       \param[in] "integer(i4), dimension(:) :: nHorizons"     Number of Horizons per soiltype
  !>       \param[in] "integer(i4), dimension(:) :: nTillHorizons" Number of Tillage Horizons
  !>       \param[in] "real(dp), dimension(:, :) :: sand"          sand content
  !>       \param[in] "real(dp), dimension(:, :) :: clay"          clay content
  !>       \param[in] "real(dp), dimension(:, :) :: DbM"           mineral Bulk density
  !>       \param[in] "integer(i4), dimension(:) :: ID0"           cell ids at level 0
  !>       \param[in] "integer(i4), dimension(:, :) :: soilId0"    soil ids at level 0
  !>       \param[in] "integer(i4), dimension(:) :: LCOVER0"       land cover ids at level 0
 
  !    INTENT(OUT)
  !>       \param[out] "real(dp), dimension(:, :, :) :: thetaS_till"  saturated soil moisture tillage layer
  !>       \param[out] "real(dp), dimension(:, :, :) :: thetaFC_till" field capacity tillage layer
  !>       \param[out] "real(dp), dimension(:, :, :) :: thetaPW_till" permanent wilting point tillage layer
  !>       \param[out] "real(dp), dimension(:, :) :: thetaS"          saturated soil moisture
  !>       \param[out] "real(dp), dimension(:, :) :: thetaFC"         field capacity
  !>       \param[out] "real(dp), dimension(:, :) :: thetaPW"         permanent wilting point
  !>       \param[out] "real(dp), dimension(:, :, :) :: Ks"           saturated hydraulic conductivity
  !>       \param[out] "real(dp), dimension(:, :, :) :: Db"           Bulk density
  !>       \param[out] "real(dp), dimension(:) :: KsVar_H0"           rel. var. of Ks for horizontal flow
  !>       \param[out] "real(dp), dimension(:) :: KsVar_V0"           rel. var. of Ks for vertical flow
  !>       \param[out] "real(dp), dimension(:) :: SMs_FC0"            soil mositure deficit from
  !>       field cap. w.r.t to saturation
 
  !    HISTORY
  !>       \authors Stephan Thober, Rohini Kumar
 
  !>       \date Dec 2012
 
  ! Modifications:
  ! Juliane Mai    Oct 2013 - OLD parametrization --> param(1) = orgMatterContent_forest
  !                                               --> param(2) = orgMatterContent_impervious
  !                                               --> param(3) = orgMatterContent_pervious
  !                                               --> param(4:13) = ... orgMatterContent_forest = orgMatterContent_perv
  !                                                                 + delta_1 NEW parametrization
  !                                               --> param(1) = delta_1 --> param(2) = orgMatterContent_impervious
  !                                               --> param(3) = orgMatterContent_pervious --> param(4:13) = ...
  ! Matthias Zink  Nov 2013 - documentation, inouts --> out moved constants to mhm_constants
  ! Stephan Thober Mar 2014 - separated cell loop from soil loop for better scaling in parallelisation
  ! David Schaefer Mar 2015 - Added dummy variable to avoid redundant computations
  !                           -> Total number of instruction is reduced by ~25%
  !                           (tested on packaged example/gnu48/{release,debug})
  ! Rohini Kumar   Mar 2016 - changes for handling multiple soil database options
  ! 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
  ! Rohini Kumar                    Oct 2021 - Neutron count module to mHM integrate into develop branch (5.11.2)
 
  subroutine mpr_sm(param, processMatrix, is_present, nHorizons, nTillHorizons, sand, clay, DbM, ID0, soilId0, LCover0, &
                   thetaS_till, thetaFC_till, thetaPW_till, thetaS, thetaFC, thetaPW, Ks, Db, KsVar_H0, KsVar_V0, SMs_FC0)
 
    use mo_common_constants, only : nodata_dp, nodata_i4
    use mo_mpr_constants, only : bulkdens_orgmatter
    use mo_mpr_global_variables, only : iflag_soildb
    !$ use omp_lib
 
    implicit none
 
    ! global parameters
    real(dp), dimension(13), intent(in) :: param
 
    ! - matrix specifying user defined processes
    integer(i4), dimension(:, :), intent(in) :: processmatrix
 
 
    ! indicates whether soiltype is present
    integer(i4), dimension(:), intent(in) :: is_present
 
    ! Number of Horizons per soiltype
    integer(i4), dimension(:), intent(in) :: nhorizons
 
    ! Number of Tillage Horizons
    integer(i4), dimension(:), intent(in) :: ntillhorizons
 
    ! sand content
    real(dp), dimension(:, :), intent(in) :: sand
 
    ! clay content
    real(dp), dimension(:, :), intent(in) :: clay
 
    ! mineral Bulk density
    real(dp), dimension(:, :), intent(in) :: dbm
 
    ! cell ids at level 0
    integer(i4), dimension(:), intent(in) :: id0
 
    ! soil ids at level 0
    integer(i4), dimension(:, :), intent(in) :: soilid0
 
    ! land cover ids at level 0
    integer(i4), dimension(:), intent(in) :: lcover0
 
    ! saturated soil moisture tillage layer
    real(dp), dimension(:, :, :), intent(out) :: thetas_till
 
    ! field capacity tillage layer
    real(dp), dimension(:, :, :), intent(out) :: thetafc_till
 
    ! permanent wilting point tillage layer
    real(dp), dimension(:, :, :), intent(out) :: thetapw_till
 
    ! saturated soil moisture
    real(dp), dimension(:, :), intent(out) :: thetas
 
    ! field capacity
    real(dp), dimension(:, :), intent(out) :: thetafc
 
    ! permanent wilting point
    real(dp), dimension(:, :), intent(out) :: thetapw
 
    ! saturated hydraulic conductivity
    real(dp), dimension(:, :, :), intent(out) :: ks
 
    ! Bulk density
    real(dp), dimension(:, :, :), intent(out) :: db
 
    ! rel. var. of Ks for horizontal flow
    real(dp), dimension(:), intent(out) :: ksvar_h0
 
    ! rel. var. of Ks for vertical flow
    real(dp), dimension(:), intent(out) :: ksvar_v0
 
    ! soil mositure deficit from
    ! field cap. w.r.t to saturation
    real(dp), dimension(:), intent(out) :: sms_fc0
 
    ! loop index
    integer(i4) :: i
 
    ! loop index
    integer(i4) :: j
 
    ! loop index
    integer(i4) :: l
 
    ! dummy variable for storing soil class
    integer(i4) :: s
 
    integer(i4) :: tmp_minsoilhorizon
 
    real(dp) :: pm
 
    real(dp) :: pom
 
    ! temporal saturated hydr. cond
    real(dp) :: ks_tmp
 
    ! van Genuchten shape param
    real(dp) :: genu_mual_n
 
    ! van Genuchten shape param
    real(dp) :: genu_mual_alpha
 
    real(dp) :: tmp_orgmattercontent_forest
 
    real(dp) :: tmp_orgmattercontent_pervious
 
    real(dp) :: tmp_orgmattercontent_impervious
 
    ! total saturated soil mositure content
    real(dp), dimension(:), allocatable :: sms_tot0
 
    ! maximum LCover class in L0
    integer(i4) :: max_lcover
 
    ! saturated hydraulic conductivity
    real(dp), dimension(:, :), allocatable :: ks_non_till
 
    ! Case 1 and 4 is based on Jarvis https://doi.org/10.1016/0022-1694(89)90050-4
    ! Case 2 and 3 is based on Feddes https://doi.org/10.1016/0022-1694(76)90017-2
    select case (processmatrix(3, 1))
    case(1,2)
       tmp_orgmattercontent_forest = param(3) + param(1)
    case(3,4)
       tmp_orgmattercontent_forest = param(1)
    end select
 
    tmp_orgmattercontent_impervious = param(2)
    tmp_orgmattercontent_pervious   = param(3)
    !write(*,*) 'tmp_orgMatterContent_forest = ', tmp_orgMatterContent_forest
    !write(*,*) 'tmp_orgMatterContent_impervious = ', tmp_orgMatterContent_impervious
    !write(*,*) 'tmp_orgMatterContent_pervious = ', tmp_orgMatterContent_pervious
 
    tmp_minsoilhorizon = minval(ntillhorizons(:))
 
    ! allocatable local variables
    ! total saturated soil moisture content
    allocate(sms_tot0(size(id0, 1)))
 
    ! some additional variables for iFlag_soil = 1
    if(iflag_soildb .eq. 1) then
      s = size(is_present, 1)
      ! note: although second dimension is not required
      ! we assign value 1 to be comparable to other assigned variables
      allocate(ks_non_till(s, 1))
    end if
 
    ! initializing soil hydraulic properties related params
    ksvar_h0 = merge(0.0_dp, nodata_dp, id0 .ne. nodata_i4)
    ksvar_v0 = merge(0.0_dp, nodata_dp, id0 .ne. nodata_i4)
    sms_tot0 = merge(0.0_dp, nodata_dp, id0 .ne. nodata_i4)
    sms_fc0 = merge(0.0_dp, nodata_dp, id0 .ne. nodata_i4)
 
    ! initialization
    thetas_till = 0.0_dp
    thetafc_till = 0.0_dp
    thetapw_till = 0.0_dp
    thetas = 0.0_dp
    thetafc = 0.0_dp
    thetapw = 0.0_dp
    ks = 0.0_dp
    db = 0.0_dp
    if(allocated(ks_non_till)) ks_non_till = 0.0_dp
    max_lcover = maxval(lcover0)
    ! select case according to a given soil database flag
    SELECT CASE(iflag_soildb)
      ! classical mHM soil database format
    CASE(0)
      !$OMP PARALLEL default(shared)
      !$OMP DO &
      !$OMP PRIVATE( i, j, L, pOM, pM, Ks_tmp, Genu_Mual_alpha, Genu_Mual_n ) &
      !$OMP SCHEDULE( STATIC )
      do i = 1, size(is_present)
        if (is_present(i) .lt. 1) cycle
        horizon : do j = 1, nhorizons(i)
          ! calculating vertical hydraulic conductivity
          call hydro_cond(ks_tmp, param(10 : 13), sand(i, j), clay(i, j))
          ks(i, j, :) = ks_tmp
          ! calculating other soil hydraulic properties
          ! tillage horizons
          if (j .le. ntillhorizons(i)) then
            ! LC class
            do l = 1, max_lcover
              select case (l)
              case(1)               ! forest
                pom = tmp_orgmattercontent_forest
              case(2)               ! impervious
                pom = tmp_orgmattercontent_impervious !param(2)
              case(3)               ! permeable
                pom = tmp_orgmattercontent_pervious
              case default
                 call error_message('Error mpr_sm: pOM used uninitialized.')
              end select
              pm = 100.0_dp - pom
              ! bulk density acording to Rawl's (1982) paper
              db(i, j, l) = 100.0_dp / ((pom / bulkdens_orgmatter) + (pm / dbm(i, j)))
              ! Effect of organic matter content
              ! This is taken into account in a simplified form by using
              ! the ratio of(Bd / BdOM)
              ks(i, j, l) = ks(i, j, l) * (dbm(i, j) / db(i, j, l))
              ! estimated SMs_till & van Genuchten's shape parameter (n)
              call genuchten(thetas_till(i, j, l), genu_mual_n, genu_mual_alpha, &
                      param(4 : 9), sand(i, j), clay(i, j), db(i, j, l))
              ! estimating field capacity
              call field_cap(thetafc_till(i, j, l), ks(i, j, l), thetas_till(i, j, l), genu_mual_n)
              ! estimating permanent wilting point
              call pwp(genu_mual_n, genu_mual_alpha, thetas_till(i, j, l), thetapw_till(i, j, l))
            end do
            ! deeper layers
          else
            ! estimate SMs & van Genuchten's shape parameter (n)
            call genuchten(thetas(i, j - tmp_minsoilhorizon), genu_mual_n, genu_mual_alpha, &
                    param(4 : 9), sand(i, j), clay(i, j), dbm(i, j))
            ! estimate field capacity
            call field_cap(thetafc(i, j - tmp_minsoilhorizon), &
                    ks(i, j, 1), thetas(i, j - tmp_minsoilhorizon), genu_mual_n)
            ! estimate permanent wilting point
            call pwp(genu_mual_n, genu_mual_alpha, thetas(i, j - tmp_minsoilhorizon), &
                    thetapw(i, j - tmp_minsoilhorizon))
          end if
        end do horizon
      end do
      !$OMP END DO
 
      ! calculate other soil properties at each location [L0] for regionalising model parameters
      !$OMP DO PRIVATE( s, j ) SCHEDULE( STATIC )
      cellloop : do i = 1, size(soilid0, 1) ! >> here = ncells0
        s = soilid0(i, 1)                    ! >> in this case the second dimension of soilId0 = 1
        do j = 1, nhorizons(s)
          if (j .le. ntillhorizons(s)) then
            ! Soil properties over the whole soil coloum depth
            ksvar_h0(i) = ksvar_h0(i) + thetas_till(s, j, lcover0(i)) * ks(s, j, lcover0(i))
            ksvar_v0(i) = ksvar_v0(i) + thetas_till(s, j, lcover0(i)) / ks(s, j, lcover0(i))
            sms_fc0(i) = sms_fc0(i) + thetafc_till(s, j, lcover0(i))
            sms_tot0(i) = sms_tot0(i) + thetas_till(s, j, lcover0(i))
          else
            ! soil_properties over the whole soil column
            ksvar_h0(i) = ksvar_h0(i) + thetas(s, j - tmp_minsoilhorizon) * ks(s, j, 1)
            ksvar_v0(i) = ksvar_v0(i) + thetas(s, j - tmp_minsoilhorizon) / ks(s, j, 1)
            sms_fc0(i) = sms_fc0(i) + thetafc(s, j - tmp_minsoilhorizon)
            sms_tot0(i) = sms_tot0(i) + thetas(s, j - tmp_minsoilhorizon)
          end if
        end do
        ! ------------------------------------------------------------------
        ! DETERMINE RELATIVE VARIABILITIES OF
        !   Ks FOR HORIZONTAL FLOW (KsVar_H)
        !               &
        !   Ks FOR VERTICAL FLOW (KsVar_V)
        ! ------------------------------------------------------------------
        ! soil moisture saturation deficit relative to the field capacity soil moisture
        sms_fc0(i) = (sms_tot0(i) - sms_fc0(i)) / sms_tot0(i)
        ! Ks variability over the whole soil coloum depth for
        ! both horizontal and vertical flows including relative variabilities
        ksvar_h0(i) = ksvar_h0(i) / sms_tot0(i) / param(13)
        ksvar_v0(i) = sms_tot0(i) / ksvar_v0(i) / param(13)
      end do cellloop
      !$OMP END DO
      !$OMP END PARALLEL
 
      ! to handle multiple soil horizons with unique soil class
    CASE(1)
      do i = 1, size(is_present)
        if (is_present(i) .lt. 1) cycle
        ! **** FOR THE TILLAGE TYPE OF SOIL *****
        ! there is actually no soil horizons/soil type in this case
        ! but we assign of j = 1 to use variables as defined in the classical option (iFlag_soil = 0)
        do j = 1, 1
          ! calculating vertical hydraulic conductivity
          call hydro_cond(ks_tmp, param(10 : 13), sand(i, j), clay(i, j))
          ks_non_till(i, j) = ks_tmp  ! >> non-till
          ks(i, j, :) = ks_tmp  ! >> till layers
          ! calculating other soil hydraulic properties
          ! tillage horizons properties depending on the LC class
          do l = 1, max_lcover
            select case (l)
            case(1)               ! forest
              pom = tmp_orgmattercontent_forest
            case(2)               ! impervious
              pom = tmp_orgmattercontent_impervious !param(2)
            case(3)               ! permeable
              pom = tmp_orgmattercontent_pervious
            case default
               call error_message('Error mpr_sm: pOM used is not initialized.')
            end select
            pm = 100.0_dp - pom
            ! bulk density acording to Rawl's (1982) paper
            db(i, j, l) = 100.0_dp / ((pom / bulkdens_orgmatter) + (pm / dbm(i, j)))
            ! Effect of organic matter content on Ks estimates
            ! This is taken into account in a simplified form by using
            ! the ratio of (Bd/BdOM)
            ks_tmp = ks_tmp * (dbm(i, j) / db(i, j, l))
            ks(i, j, l) = ks_tmp
            ! estimated SMs_till & van Genuchten's shape parameter (n)
            call genuchten(thetas_till(i, j, l), genu_mual_n, genu_mual_alpha, &
                    param(4 : 9), sand(i, j), clay(i, j), db(i, j, l))
            ! estimating field capacity
            call field_cap(thetafc_till(i, j, l), ks_tmp, thetas_till(i, j, l), genu_mual_n)
            ! estimating permanent wilting point
            call pwp(genu_mual_n, genu_mual_alpha, thetas_till(i, j, l), thetapw_till(i, j, l))
          end do
 
          ! *** FOR NON-TILLAGE TYPE OF SOILS ***
          ! note j = 1
          ! since Ks_tmp has changed earlier ... get the original Ks once again
          ks_tmp = ks_non_till(i, j)
          ! estimate SMs & van Genuchten's shape parameter (n)
          call genuchten(thetas(i, j), genu_mual_n, genu_mual_alpha, param(4 : 9), sand(i, j), clay(i, j), dbm(i, j))
          ! estimate field capacity
          call field_cap(thetafc(i, j), ks_tmp, thetas(i, j), genu_mual_n)
          ! estimate permanent wilting point
          call pwp(genu_mual_n, genu_mual_alpha, thetas(i, j), thetapw(i, j))
 
        end do  ! >> HORIZON
      end do   ! >> SOIL TYPE
 
      ! calculate other soil properties at each location [L0] for regionalising model parameters
      do i = 1, size(soilid0, 1)     !! over all cells
        do j = 1, size(soilid0, 2)  !! over horizons
          s = soilid0(i, j)
          if (j .le. ntillhorizons(1)) then
            ! soil properties over the whole soil coloum depth
            ksvar_h0(i) = ksvar_h0(i) + thetas_till(s, 1, lcover0(i)) * ks(s, 1, lcover0(i))
            ksvar_v0(i) = ksvar_v0(i) + thetas_till(s, 1, lcover0(i)) / ks(s, 1, lcover0(i))
            sms_fc0(i) = sms_fc0(i) + thetafc_till(s, 1, lcover0(i))
            sms_tot0(i) = sms_tot0(i) + thetas_till(s, 1, lcover0(i))
          else
            ! soil_properties over the whole soil column
            ksvar_h0(i) = ksvar_h0(i) + thetas(s, 1) * ks_non_till(s, 1)
            ksvar_v0(i) = ksvar_v0(i) + thetas(s, 1) / ks_non_till(s, 1)
            sms_fc0(i) = sms_fc0(i) + thetafc(s, 1)
            sms_tot0(i) = sms_tot0(i) + thetas(s, 1)
          end if
        end do
        ! ------------------------------------------------------------------
        ! DETERMINE RELATIVE VARIABILITIES OF
        ! Ks FOR HORIZONTAL FLOW (KsVar_H) & Ks FOR VERTICAL FLOW (KsVar_V)
        ! ------------------------------------------------------------------
        ! soil moisture saturation deficit relative to the field capacity soil moisture
        sms_fc0(i) = (sms_tot0(i) - sms_fc0(i)) / sms_tot0(i)
        ! Ks variability over the whole soil coloum depth for
        ! both horizontal and vertical flows including relative variabilities
        ksvar_h0(i) = ksvar_h0(i) / sms_tot0(i) / param(13)
        ksvar_v0(i) = sms_tot0(i) / ksvar_v0(i) / param(13)
      end do
 
    CASE DEFAULT
      call error_message('***ERROR: iFlag_soilDB option given does not exist. Only 0 and 1 is taken at the moment.')
    END SELECT
 
    ! free space **
    deallocate(sms_tot0)
    if(allocated(ks_non_till)) deallocate(ks_non_till)
 
  end subroutine mpr_sm
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        PWP
 
  !    PURPOSE
  !>       \brief Permanent Wilting point
 
  !>       \details This subroutine calculates the permanent wilting
  !>       point according to Zacharias et al. (2007, Soil Phy.) and
  !>       using van Genuchten 1980's equation. For the water retention curve at
  !>       a matrix potential of -1500 kPa, it is assumed that thetaR = 0.
  !>       ADDITIONAL INFORMATION
  !>       Zacharias et al. 2007, Soil Phy.
 
  !    INTENT(IN)
  !>       \param[in] "real(dp) :: Genu_Mual_n"     - Genuchten shape parameter
  !>       \param[in] "real(dp) :: Genu_Mual_alpha" - Genuchten shape parameter
  !>       \param[in] "real(dp) :: thetaS"          - saturated water content
 
  !    INTENT(OUT)
  !>       \param[out] "real(dp) :: thetaPWP" - Permanent Wilting point
 
  !    HISTORY
  !>       \authors Stephan Thober, Rohini Kumar
 
  !>       \date Dec, 2012
 
  ! Modifications:
  ! Matthias Zink Nov 2013 - documentation, moved constants to mhm_constants
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  elemental pure subroutine pwp(Genu_Mual_n, Genu_Mual_alpha, thetaS, thetaPWP)
 
    use mo_mpr_constants, only : pwp_c, pwp_matpot_thetar
 
    implicit none
 
    ! - Genuchten shape parameter
    real(dp), intent(in) :: genu_mual_n
 
    ! - Genuchten shape parameter
    real(dp), intent(in) :: genu_mual_alpha
 
    ! - saturated water content
    real(dp), intent(in) :: thetas
 
    ! - Permanent Wilting point
    real(dp), intent(out) :: thetapwp
 
    real(dp) :: x
 
    ! Genuchten shape parameter
    real(dp) :: genu_mual_m
 
 
    genu_mual_m = pwp_c - (pwp_c / genu_mual_n)
    x = pwp_c + exp(genu_mual_n * log(genu_mual_alpha * pwp_matpot_thetar))
    x = exp(genu_mual_m * log(x))
    ! constrain
    if (x < 1.0_dp) x = 1.0_dp
    thetapwp = thetas / x
 
  end subroutine pwp
 
  ! ----------------------------------------------------------------------------
 
  !    NAME
  !        field_cap
 
  !    PURPOSE
  !>       \brief calculates the field capacity
 
  !>       \details estimate Field capacity; FC -- Flux based
  !>       approach (Twarakavi, et. al. 2009, WRR)
  !>       According to the
  !>       above reference FC is defined as the soil water content at
  !>       which the drainage from a profile ceases under natural
  !>       conditions. Since drainage from a soil profile in a simulation
  !>       never becomes zero, we assume that drainage ceases when the
  !>       bottom flux from the soil reaches a value that is equivalent to
  !>       the minimum amount of precipitation that could be recorded
  !>       (i.e. 0.01 cm/d == 1 mm/d). It is assumed that ThetaR = 0.0_dp
  !>       ADDITIONAL INFORMATION
  !>       Twarakavi, et. al. 2009, WRR
 
  !    INTENT(OUT)
  !>       \param[out] "real(dp) :: thetaFC" - Field capacity
 
  !    INTENT(IN)
  !>       \param[in] "real(dp) :: Ks"          - saturated hydraulic conductivity
  !>       \param[in] "real(dp) :: thetaS"      - saturated water content
  !>       \param[in] "real(dp) :: Genu_Mual_n" - Genuchten shape parameter
 
  !    HISTORY
  !>       \authors Stephan Thober, Rohini Kumar
 
  !>       \date Dec 2012
 
  ! Modifications:
  ! Matthias Zink Nov 2013 - documentation, moved constants to mhm_constants
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  elemental pure subroutine field_cap(thetaFC, Ks, thetaS, Genu_Mual_n)
 
    use mo_mpr_constants, only : field_cap_c1, field_cap_c2
 
    implicit none
 
    ! - saturated hydraulic conductivity
    real(dp), intent(in) :: ks
 
    ! - saturated water content
    real(dp), intent(in) :: thetas
 
    ! - Genuchten shape parameter
    real(dp), intent(in) :: genu_mual_n
 
    ! - Field capacity
    real(dp), intent(out) :: thetafc
 
    real(dp) :: x
 
 
    x = (field_cap_c1) * (field_cap_c2 + log10(ks))
    thetafc = thetas * exp(x * log(genu_mual_n))
 
  end subroutine field_cap
 
  ! ----------------------------------------------------------------------------
 
  !    NAME
  !        Genuchten
 
  !    PURPOSE
  !>       \brief calculates the Genuchten shape parameter
 
  !>       \details estimate SMs_till & van Genuchten's shape parameter (n)
  !>       (Zacharias et al, 2007, soil Phy.)
  !>       Global parameters needed (see mhm_parameter.nml):
  !>       - param( 1) = PTF_lower66_5_constant
  !>       - param( 2) = PTF_lower66_5_clay
  !>       - param( 3) = PTF_lower66_5_Db
  !>       - param( 4) = PTF_higher66_5_constant
  !>       - param( 5) = PTF_higher66_5_clay
  !>       - param( 6) = PTF_higher66_5_Db
  !>       ADDITIONAL INFORMATION
  !>       Zacharias et al, 2007, soil Phy.
 
  !    INTENT(OUT)
  !>       \param[out] "real(dp) :: thetaS"          - saturated water content
  !>       \param[out] "real(dp) :: Genu_Mual_n"     - van Genuchten shape parameter
  !>       \param[out] "real(dp) :: Genu_Mual_alpha" - van Genuchten shape parameter
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(6) :: param" parameters
  !>       \param[in] "real(dp) :: sand"                - [%] sand content
  !>       \param[in] "real(dp) :: clay"                - [%] clay content
  !>       \param[in] "real(dp) :: Db"                  - [10^3 kg/m3] bulk density
 
  !    HISTORY
  !>       \authors Stephan Thober, Rohini Kumar
 
  !>       \date Dec 2012
 
  ! Modifications:
  ! Rohini Kumar Mar 2014 - ThetaS limit changed from 0 to 0.001
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine genuchten(thetaS, Genu_Mual_n, Genu_Mual_alpha, param, sand, clay, Db)
 
    use mo_mpr_constants, only : vgenuchtenn_c1, vgenuchtenn_c10, vgenuchtenn_c11, vgenuchtenn_c12, vgenuchtenn_c13, &
                                 vgenuchtenn_c14, vgenuchtenn_c15, vgenuchtenn_c16, vgenuchtenn_c17, vgenuchtenn_c18, &
                                 vgenuchtenn_c2, vgenuchtenn_c3, vgenuchtenn_c4, vgenuchtenn_c5, vgenuchtenn_c6, &
                                 vgenuchtenn_c7, vgenuchtenn_c8, vgenuchtenn_c9, vgenuchten_sandtresh
 
    implicit none
 
    ! parameters
    real(dp), dimension(6), intent(in) :: param
 
    ! - [%] sand content
    real(dp), intent(in) :: sand
 
    ! - [%] clay content
    real(dp), intent(in) :: clay
 
    ! - [10^3 kg/m3] bulk density
    real(dp), intent(in) :: Db
 
    ! - saturated water content
    real(dp), intent(out) :: thetaS
 
    ! - van Genuchten shape parameter
    real(dp), intent(out) :: Genu_Mual_n
 
    ! - van Genuchten shape parameter
    real(dp), intent(out) :: Genu_Mual_alpha
 
    ! temporal variable
    real(dp) :: x
 
 
    ! estimate SMs_till & van Genuchten's parameters (alpha and n)
    if (sand < vgenuchten_sandtresh) then
      thetas = param(1) + param(2) * clay + param(3) * db
      genu_mual_n = vgenuchtenn_c1 - vgenuchtenn_c2 * (sand**(vgenuchtenn_c3)) + &
              vgenuchtenn_c4 * (clay**(vgenuchtenn_c5))
      x = vgenuchtenn_c6 + vgenuchtenn_c7 * sand + vgenuchtenn_c8 * clay - &
              vgenuchtenn_c9 * db
    else
      thetas = param(4) + param(5) * clay + param(6) * db
      genu_mual_n = vgenuchtenn_c10 + vgenuchtenn_c11 * (sand**(vgenuchtenn_c12)) + &
              vgenuchtenn_c13 * (clay**(vgenuchtenn_c14))
      x = vgenuchtenn_c15 + vgenuchtenn_c16 * sand + vgenuchtenn_c17 * clay - &
              vgenuchtenn_c18 * db
    end if
 
    ! Mualem alpha
    genu_mual_alpha = exp(x)
 
    ! hard coded limits, according to (Zacharias et al, 2007, soil Phy.)
    if (thetas < 0.01_dp) then
      call message('thetaS below threshold limit 1e-2, reset.')
      ! Put constrains on theta_S
      thetas = 0.01_dp
    end if
    if (thetas > 1.0_dp) then
      call message('thetaS above 1, reset.')
      ! Put constrains on theta_S
      thetas = 1.0_dp
    end if
    if (genu_mual_n < 1.01000_dp) then
      call message('Genu_Mual_n below threshold limit 1.01, reset.')
      genu_mual_n = 1.01000_dp
    end if
    if (genu_mual_alpha < 0.00001_dp) then
      call message('Genu_Mual_alpha below threshold limit 1e-5, reset.')
      genu_mual_alpha = 0.00001_dp
    end if
 
  end subroutine genuchten
 
  ! ----------------------------------------------------------------------------
 
  !    NAME
  !        hydro_cond
 
  !    PURPOSE
  !>       \brief calculates the hydraulic conductivity Ks
 
  !>       \details By default save this value of Ks, particularly for the
  !>       deeper layers where OM content plays relatively low or no role
  !>       Global parameters needed (see mhm_parameter.nml):
  !>       - param(1) = PTF_Ks_constant
  !>       - param(2) = PTF_Ks_sand
  !>       - param(3) = PTF_Ks_clay
  !>       - param(4) = PTF_Ks_curveSlope
  !>       ADDITIONAL INFORMATION
  !>       Written,  Stephan Thober, Dec 2012
 
  !    INTENT(OUT)
  !>       \param[out] "real(dp) :: KS"
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(4) :: param"
  !>       \param[in] "real(dp) :: sand"                - [%] sand content
  !>       \param[in] "real(dp) :: clay"                - [%] clay content
 
  !    HISTORY
  !>       \authors Stephan Thober, Rohini Kumar
 
  !>       \date Dec 2012
 
  ! Modifications:
  ! Matthias Zink  Nov 2013 - documentation, moved constants to mhm_constants
  ! Matthias Cuntz Jun 2014 - suggested to fix param(4)
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine hydro_cond(KS, param, sand, clay)
 
    use mo_mpr_constants, only : ks_c
 
    implicit none
 
    real(dp), dimension(4), intent(in) :: param
 
    ! - [%] sand content
    real(dp), intent(in) :: sand
 
    ! - [%] clay content
    real(dp), intent(in) :: clay
 
    real(dp), intent(out) :: KS
 
    ! temporal variable
    real(dp) :: x
 
 
    ! saturated vertical hydraulic conductivity, Ks (cm/d)
    !   from Cosby et. al. (WRR 1984) Table 4
    ! param(4) is the unit conversion from inch/h to cm/d and should be a constant.
    ! Fix it in the namelist, i.e. in
    ! mhm_parameter.nml set the 4th value (=FLAG) to 0 and the third value to 60.96
    !   PTF_Ks_curveSlope = 60.96, 60.96, 60.96, 0, 1
    x = param(1) + param(2) * sand - param(3) * clay
    ks = param(4) * exp(x * log(ks_c))
 
    if (ks < 1.10_dp) then
      call message('JMJMJM-Ks-BAD')
    end if
 
    ! minimum value of Ks = 1.1cm/d
    if (ks < 1.10_dp) ks = 1.10_dp
 
  end subroutine hydro_cond
 
end module mo_mpr_soilmoist