LCOV - code coverage report
Current view: top level - MPR - mo_mpr_soilmoist.f90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 116 166 69.9 %
Date: 2024-04-30 08:53:32 Functions: 5 5 100.0 %

          Line data    Source code
       1             : !> \file mo_mpr_soilmoist.f90
       2             : !> \brief \copybrief mo_mpr_soilmoist
       3             : !> \details \copydetails mo_mpr_soilmoist
       4             : 
       5             : !> \brief Multiscale parameter regionalization (MPR) for soil moisture
       6             : !> \details This module contains all routines required for parametrizing soil moisture processes.
       7             : !> \authors Stephan Thober, Rohini Kumar
       8             : !> \date Dec 2012
       9             : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
      10             : !! mHM is released under the LGPLv3+ license \license_note
      11             : !> \ingroup f_mpr
      12             : module mo_mpr_soilmoist
      13             : 
      14             :   use mo_kind, only : i4, dp
      15             :   use mo_message, only : message, error_message
      16             : 
      17             :   implicit none
      18             : 
      19             :   public :: mpr_sm
      20             : 
      21             :   private
      22             : 
      23             : contains
      24             :   ! ----------------------------------------------------------------------------
      25             : 
      26             :   !    NAME
      27             :   !        mpr_sm
      28             : 
      29             :   !    PURPOSE
      30             :   !>       \brief multiscale parameter regionalization for soil moisture
      31             : 
      32             :   !>       \details This subroutine is a wrapper around all soil moisture
      33             :   !>       parameter routines. This subroutine requires 13 parameters. These
      34             :   !>       parameters have to correspond to the parameters in the original
      35             :   !>       parameter array at the following locations: 10-12, 13-18, 27-30.
      36             :   !>       Global parameters needed (see mhm_parameter.nml):
      37             :   !>       - param( 1) = orgMatterContent_forest
      38             :   !>       - param( 2) = orgMatterContent_impervious
      39             :   !>       - param( 3) = orgMatterContent_pervious
      40             :   !>       - param( 4) = PTF_lower66_5_constant
      41             :   !>       - param( 5) = PTF_lower66_5_clay
      42             :   !>       - param( 6) = PTF_lower66_5_Db
      43             :   !>       - param( 7) = PTF_higher66_5_constant
      44             :   !>       - param( 8) = PTF_higher66_5_clay
      45             :   !>       - param( 9) = PTF_higher66_5_Db
      46             :   !>       - param(10) = PTF_Ks_constant
      47             :   !>       - param(11) = PTF_Ks_sand
      48             :   !>       - param(12) = PTF_Ks_clay
      49             :   !>       - param(13) = PTF_Ks_curveSlope
      50             : 
      51             :   !    INTENT(IN)
      52             :   !>       \param[in] "real(dp), dimension(13) :: param"           global parameters
      53             :   !>       \param[in] "integer(i4), dimension(:) :: is_present"    indicates whether soiltype is present
      54             :   !>       \param[in] "integer(i4), dimension(:) :: nHorizons"     Number of Horizons per soiltype
      55             :   !>       \param[in] "integer(i4), dimension(:) :: nTillHorizons" Number of Tillage Horizons
      56             :   !>       \param[in] "real(dp), dimension(:, :) :: sand"          sand content
      57             :   !>       \param[in] "real(dp), dimension(:, :) :: clay"          clay content
      58             :   !>       \param[in] "real(dp), dimension(:, :) :: DbM"           mineral Bulk density
      59             :   !>       \param[in] "integer(i4), dimension(:) :: ID0"           cell ids at level 0
      60             :   !>       \param[in] "integer(i4), dimension(:, :) :: soilId0"    soil ids at level 0
      61             :   !>       \param[in] "integer(i4), dimension(:) :: LCOVER0"       land cover ids at level 0
      62             : 
      63             :   !    INTENT(OUT)
      64             :   !>       \param[out] "real(dp), dimension(:, :, :) :: thetaS_till"  saturated soil moisture tillage layer
      65             :   !>       \param[out] "real(dp), dimension(:, :, :) :: thetaFC_till" field capacity tillage layer
      66             :   !>       \param[out] "real(dp), dimension(:, :, :) :: thetaPW_till" permanent wilting point tillage layer
      67             :   !>       \param[out] "real(dp), dimension(:, :) :: thetaS"          saturated soil moisture
      68             :   !>       \param[out] "real(dp), dimension(:, :) :: thetaFC"         field capacity
      69             :   !>       \param[out] "real(dp), dimension(:, :) :: thetaPW"         permanent wilting point
      70             :   !>       \param[out] "real(dp), dimension(:, :, :) :: Ks"           saturated hydraulic conductivity
      71             :   !>       \param[out] "real(dp), dimension(:, :, :) :: Db"           Bulk density
      72             :   !>       \param[out] "real(dp), dimension(:) :: KsVar_H0"           rel. var. of Ks for horizontal flow
      73             :   !>       \param[out] "real(dp), dimension(:) :: KsVar_V0"           rel. var. of Ks for vertical flow
      74             :   !>       \param[out] "real(dp), dimension(:) :: SMs_FC0"            soil mositure deficit from
      75             :   !>       field cap. w.r.t to saturation
      76             : 
      77             :   !    HISTORY
      78             :   !>       \authors Stephan Thober, Rohini Kumar
      79             : 
      80             :   !>       \date Dec 2012
      81             : 
      82             :   ! Modifications:
      83             :   ! Juliane Mai    Oct 2013 - OLD parametrization --> param(1) = orgMatterContent_forest
      84             :   !                                               --> param(2) = orgMatterContent_impervious
      85             :   !                                               --> param(3) = orgMatterContent_pervious
      86             :   !                                               --> param(4:13) = ... orgMatterContent_forest = orgMatterContent_perv
      87             :   !                                                                 + delta_1 NEW parametrization
      88             :   !                                               --> param(1) = delta_1 --> param(2) = orgMatterContent_impervious
      89             :   !                                               --> param(3) = orgMatterContent_pervious --> param(4:13) = ...
      90             :   ! Matthias Zink  Nov 2013 - documentation, inouts --> out moved constants to mhm_constants
      91             :   ! Stephan Thober Mar 2014 - separated cell loop from soil loop for better scaling in parallelisation
      92             :   ! David Schaefer Mar 2015 - Added dummy variable to avoid redundant computations
      93             :   !                           -> Total number of instruction is reduced by ~25%
      94             :   !                           (tested on packaged example/gnu48/{release,debug})
      95             :   ! Rohini Kumar   Mar 2016 - changes for handling multiple soil database options
      96             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
      97             :   ! M. Cuneyd Demirel, Simon Stisen Jun 2020 - added Feddes and FC dependency on root fraction coefficient processCase(3) = 4
      98             :   ! Rohini Kumar                    Oct 2021 - Neutron count module to mHM integrate into develop branch (5.11.2)
      99             : 
     100        1242 :   subroutine mpr_sm(param, processMatrix, is_present, nHorizons, nTillHorizons, sand, clay, DbM, ID0, soilId0, LCover0, &
     101         414 :                    thetaS_till, thetaFC_till, thetaPW_till, thetaS, thetaFC, thetaPW, Ks, Db, KsVar_H0, KsVar_V0, SMs_FC0)
     102             : 
     103             :     use mo_common_constants, only : nodata_dp, nodata_i4
     104             :     use mo_mpr_constants, only : BulkDens_OrgMatter
     105             :     use mo_mpr_global_variables, only : iFlag_soilDB
     106             :     !$ use omp_lib
     107             : 
     108             :     implicit none
     109             : 
     110             :     ! global parameters
     111             :     real(dp), dimension(13), intent(in) :: param
     112             : 
     113             :     ! - matrix specifying user defined processes
     114             :     integer(i4), dimension(:, :), intent(in) :: processMatrix
     115             : 
     116             : 
     117             :     ! indicates whether soiltype is present
     118             :     integer(i4), dimension(:), intent(in) :: is_present
     119             : 
     120             :     ! Number of Horizons per soiltype
     121             :     integer(i4), dimension(:), intent(in) :: nHorizons
     122             : 
     123             :     ! Number of Tillage Horizons
     124             :     integer(i4), dimension(:), intent(in) :: nTillHorizons
     125             : 
     126             :     ! sand content
     127             :     real(dp), dimension(:, :), intent(in) :: sand
     128             : 
     129             :     ! clay content
     130             :     real(dp), dimension(:, :), intent(in) :: clay
     131             : 
     132             :     ! mineral Bulk density
     133             :     real(dp), dimension(:, :), intent(in) :: DbM
     134             : 
     135             :     ! cell ids at level 0
     136             :     integer(i4), dimension(:), intent(in) :: ID0
     137             : 
     138             :     ! soil ids at level 0
     139             :     integer(i4), dimension(:, :), intent(in) :: soilId0
     140             : 
     141             :     ! land cover ids at level 0
     142             :     integer(i4), dimension(:), intent(in) :: LCOVER0
     143             : 
     144             :     ! saturated soil moisture tillage layer
     145             :     real(dp), dimension(:, :, :), intent(out) :: thetaS_till
     146             : 
     147             :     ! field capacity tillage layer
     148             :     real(dp), dimension(:, :, :), intent(out) :: thetaFC_till
     149             : 
     150             :     ! permanent wilting point tillage layer
     151             :     real(dp), dimension(:, :, :), intent(out) :: thetaPW_till
     152             : 
     153             :     ! saturated soil moisture
     154             :     real(dp), dimension(:, :), intent(out) :: thetaS
     155             : 
     156             :     ! field capacity
     157             :     real(dp), dimension(:, :), intent(out) :: thetaFC
     158             : 
     159             :     ! permanent wilting point
     160             :     real(dp), dimension(:, :), intent(out) :: thetaPW
     161             : 
     162             :     ! saturated hydraulic conductivity
     163             :     real(dp), dimension(:, :, :), intent(out) :: Ks
     164             : 
     165             :     ! Bulk density
     166             :     real(dp), dimension(:, :, :), intent(out) :: Db
     167             : 
     168             :     ! rel. var. of Ks for horizontal flow
     169             :     real(dp), dimension(:), intent(out) :: KsVar_H0
     170             : 
     171             :     ! rel. var. of Ks for vertical flow
     172             :     real(dp), dimension(:), intent(out) :: KsVar_V0
     173             : 
     174             :     ! soil mositure deficit from
     175             :     ! field cap. w.r.t to saturation
     176             :     real(dp), dimension(:), intent(out) :: SMs_FC0
     177             : 
     178             :     ! loop index
     179             :     integer(i4) :: i
     180             : 
     181             :     ! loop index
     182             :     integer(i4) :: j
     183             : 
     184             :     ! loop index
     185             :     integer(i4) :: l
     186             : 
     187             :     ! dummy variable for storing soil class
     188             :     integer(i4) :: s
     189             : 
     190             :     integer(i4) :: tmp_minSoilHorizon
     191             : 
     192         414 :     real(dp) :: pM
     193             : 
     194         414 :     real(dp) :: pOM
     195             : 
     196             :     ! temporal saturated hydr. cond
     197         414 :     real(dp) :: Ks_tmp
     198             : 
     199             :     ! van Genuchten shape param
     200         414 :     real(dp) :: Genu_Mual_n
     201             : 
     202             :     ! van Genuchten shape param
     203         414 :     real(dp) :: Genu_Mual_alpha
     204             : 
     205         414 :     real(dp) :: tmp_orgMatterContent_forest
     206             : 
     207         414 :     real(dp) :: tmp_orgMatterContent_pervious
     208             : 
     209         414 :     real(dp) :: tmp_orgMatterContent_impervious
     210             : 
     211             :     ! total saturated soil mositure content
     212         414 :     real(dp), dimension(:), allocatable :: SMs_tot0
     213             : 
     214             :     ! maximum LCover class in L0
     215             :     integer(i4) :: max_LCover
     216             : 
     217             :     ! saturated hydraulic conductivity
     218         414 :     real(dp), dimension(:, :), allocatable :: Ks_non_till
     219             : 
     220             :     ! Case 1 and 4 is based on Jarvis https://doi.org/10.1016/0022-1694(89)90050-4
     221             :     ! Case 2 and 3 is based on Feddes https://doi.org/10.1016/0022-1694(76)90017-2
     222         414 :     select case (processMatrix(3, 1))
     223             :     case(1,2)
     224         410 :        tmp_orgMatterContent_forest = param(3) + param(1)
     225             :     case(3,4)
     226         414 :        tmp_orgMatterContent_forest = param(1)
     227             :     end select
     228             : 
     229         414 :     tmp_orgMatterContent_impervious = param(2)
     230         414 :     tmp_orgMatterContent_pervious   = param(3)
     231             :     !write(*,*) 'tmp_orgMatterContent_forest = ', tmp_orgMatterContent_forest
     232             :     !write(*,*) 'tmp_orgMatterContent_impervious = ', tmp_orgMatterContent_impervious
     233             :     !write(*,*) 'tmp_orgMatterContent_pervious = ', tmp_orgMatterContent_pervious
     234             : 
     235      611064 :     tmp_minSoilHorizon = minval(nTillHorizons(:))
     236             : 
     237             :     ! allocatable local variables
     238             :     ! total saturated soil moisture content
     239        1242 :     allocate(SMs_tot0(size(ID0, 1)))
     240             : 
     241             :     ! some additional variables for iFlag_soil = 1
     242         414 :     if(iFlag_soilDB .eq. 1) then
     243           0 :       s = size(is_present, 1)
     244             :       ! note: although second dimension is not required
     245             :       ! we assign value 1 to be comparable to other assigned variables
     246           0 :       allocate(Ks_non_till(s, 1))
     247             :     end if
     248             : 
     249             :     ! initializing soil hydraulic properties related params
     250    19168378 :     KsVar_H0 = merge(0.0_dp, nodata_dp, ID0 .ne. nodata_i4)
     251    19168378 :     KsVar_V0 = merge(0.0_dp, nodata_dp, ID0 .ne. nodata_i4)
     252    19168378 :     SMs_tot0 = merge(0.0_dp, nodata_dp, ID0 .ne. nodata_i4)
     253    19168378 :     SMs_FC0 = merge(0.0_dp, nodata_dp, ID0 .ne. nodata_i4)
     254             : 
     255             :     ! initialization
     256     1834848 :     thetaS_till = 0.0_dp
     257     1834848 :     thetaFC_till = 0.0_dp
     258     1834848 :     thetaPW_till = 0.0_dp
     259     1833606 :     thetaS = 0.0_dp
     260     1833606 :     thetaFC = 0.0_dp
     261     1833606 :     thetaPW = 0.0_dp
     262     5501232 :     Ks = 0.0_dp
     263     5501232 :     Db = 0.0_dp
     264         414 :     if(allocated(Ks_non_till)) Ks_non_till = 0.0_dp
     265    19167964 :     max_LCover = maxval(LCOVER0)
     266             :     ! select case according to a given soil database flag
     267         414 :     SELECT CASE(iFlag_soilDB)
     268             :       ! classical mHM soil database format
     269             :     CASE(0)
     270             :       !$OMP PARALLEL default(shared)
     271             :       !$OMP DO &
     272             :       !$OMP PRIVATE( i, j, L, pOM, pM, Ks_tmp, Genu_Mual_alpha, Genu_Mual_n ) &
     273             :       !$OMP SCHEDULE( STATIC )
     274      611064 :       do i = 1, size(is_present)
     275      610650 :         if (is_present(i) .lt. 1) cycle
     276       57318 :         horizon : do j = 1, nHorizons(i)
     277             :           ! calculating vertical hydraulic conductivity
     278       42678 :           call hydro_cond(Ks_tmp, param(10 : 13), sand(i, j), clay(i, j))
     279      170712 :           Ks(i, j, :) = Ks_tmp
     280             :           ! calculating other soil hydraulic properties
     281             :           ! tillage horizons
     282      653328 :           if (j .le. nTillHorizons(i)) then
     283             :             ! LC class
     284       56904 :             do L = 1, max_LCover
     285             :               select case (L)
     286             :               case(1)               ! forest
     287       14226 :                 pOM = tmp_orgMatterContent_forest
     288             :               case(2)               ! impervious
     289       14226 :                 pOM = tmp_orgMatterContent_impervious !param(2)
     290             :               case(3)               ! permeable
     291       14226 :                 pOM = tmp_orgMatterContent_pervious
     292             :               case default
     293       42678 :                  call error_message('Error mpr_sm: pOM used uninitialized.')
     294             :               end select
     295       42678 :               pM = 100.0_dp - pOM
     296             :               ! bulk density acording to Rawl's (1982) paper
     297       42678 :               Db(i, j, L) = 100.0_dp / ((pOM / BulkDens_OrgMatter) + (pM / DbM(i, j)))
     298             :               ! Effect of organic matter content
     299             :               ! This is taken into account in a simplified form by using
     300             :               ! the ratio of(Bd / BdOM)
     301       42678 :               Ks(i, j, L) = Ks(i, j, L) * (DbM(i, j) / Db(i, j, L))
     302             :               ! estimated SMs_till & van Genuchten's shape parameter (n)
     303       42678 :               call Genuchten(thetaS_till(i, j, L), Genu_Mual_n, Genu_Mual_alpha, &
     304       42678 :                       param(4 : 9), sand(i, j), clay(i, j), Db(i, j, L))
     305             :               ! estimating field capacity
     306       42678 :               call field_cap(thetaFC_till(i, j, L), Ks(i, j, L), thetaS_till(i, j, L), Genu_Mual_n)
     307             :               ! estimating permanent wilting point
     308       56904 :               call PWP(Genu_Mual_n, Genu_Mual_alpha, thetaS_till(i, j, L), thetaPW_till(i, j, L))
     309             :             end do
     310             :             ! deeper layers
     311             :           else
     312             :             ! estimate SMs & van Genuchten's shape parameter (n)
     313       28452 :             call Genuchten(thetaS(i, j - tmp_minSoilHorizon), Genu_Mual_n, Genu_Mual_alpha, &
     314       56904 :                     param(4 : 9), sand(i, j), clay(i, j), DbM(i, j))
     315             :             ! estimate field capacity
     316       28452 :             call field_cap(thetaFC(i, j - tmp_minSoilHorizon), &
     317       28452 :                     Ks(i, j, 1), thetaS(i, j - tmp_minSoilHorizon), Genu_Mual_n)
     318             :             ! estimate permanent wilting point
     319             :             call PWP(Genu_Mual_n, Genu_Mual_alpha, thetaS(i, j - tmp_minSoilHorizon), &
     320       28452 :                     thetaPW(i, j - tmp_minSoilHorizon))
     321             :           end if
     322             :         end do horizon
     323             :       end do
     324             :       !$OMP END DO
     325             : 
     326             :       ! calculate other soil properties at each location [L0] for regionalising model parameters
     327             :       !$OMP DO PRIVATE( s, j ) SCHEDULE( STATIC )
     328    19168378 :       cellloop : do i = 1, size(soilId0, 1) ! >> here = ncells0
     329    19167550 :         s = soilId0(i, 1)                    ! >> in this case the second dimension of soilId0 = 1
     330    76670200 :         do j = 1, nHorizons(s)
     331    76670200 :           if (j .le. nTillHorizons(s)) then
     332             :             ! Soil properties over the whole soil coloum depth
     333    19167550 :             KsVar_H0(i) = KsVar_H0(i) + thetaS_till(s, j, LCover0(i)) * Ks(s, j, LCover0(i))
     334    19167550 :             KsVar_V0(i) = KsVar_V0(i) + thetaS_till(s, j, LCover0(i)) / Ks(s, j, LCover0(i))
     335    19167550 :             SMs_FC0(i) = SMs_FC0(i) + thetaFC_till(s, j, LCover0(i))
     336    19167550 :             SMs_tot0(i) = SMs_tot0(i) + thetaS_till (s, j, LCover0(i))
     337             :           else
     338             :             ! soil_properties over the whole soil column
     339    38335100 :             KsVar_H0(i) = KsVar_H0(i) + thetaS(s, j - tmp_minSoilHorizon) * Ks(s, j, 1)
     340    38335100 :             KsVar_V0(i) = KsVar_V0(i) + thetaS(s, j - tmp_minSoilHorizon) / Ks(s, j, 1)
     341    38335100 :             SMs_FC0(i) = SMs_FC0(i) + thetaFC(s, j - tmp_minSoilHorizon)
     342    38335100 :             SMs_tot0(i) = SMs_tot0(i) + thetaS (s, j - tmp_minSoilHorizon)
     343             :           end if
     344             :         end do
     345             :         ! ------------------------------------------------------------------
     346             :         ! DETERMINE RELATIVE VARIABILITIES OF
     347             :         !   Ks FOR HORIZONTAL FLOW (KsVar_H)
     348             :         !               &
     349             :         !   Ks FOR VERTICAL FLOW (KsVar_V)
     350             :         ! ------------------------------------------------------------------
     351             :         ! soil moisture saturation deficit relative to the field capacity soil moisture
     352    19167550 :         SMs_FC0(i) = (SMs_tot0(i) - SMs_FC0(i)) / SMs_tot0(i)
     353             :         ! Ks variability over the whole soil coloum depth for
     354             :         ! both horizontal and vertical flows including relative variabilities
     355    19167550 :         KsVar_H0(i) = KsVar_H0(i) / SMs_tot0(i) / param(13)
     356    19167964 :         KsVar_V0(i) = SMs_tot0(i) / KsVar_V0(i) / param(13)
     357             :       end do cellloop
     358             :       !$OMP END DO
     359             :       !$OMP END PARALLEL
     360             : 
     361             :       ! to handle multiple soil horizons with unique soil class
     362             :     CASE(1)
     363           0 :       do i = 1, size(is_present)
     364           0 :         if (is_present(i) .lt. 1) cycle
     365             :         ! **** FOR THE TILLAGE TYPE OF SOIL *****
     366             :         ! there is actually no soil horizons/soil type in this case
     367             :         ! but we assign of j = 1 to use variables as defined in the classical option (iFlag_soil = 0)
     368           0 :         do j = 1, 1
     369             :           ! calculating vertical hydraulic conductivity
     370           0 :           call hydro_cond(Ks_tmp, param(10 : 13), sand(i, j), clay(i, j))
     371           0 :           Ks_non_till(i, j) = Ks_tmp  ! >> non-till
     372           0 :           Ks(i, j, :) = Ks_tmp  ! >> till layers
     373             :           ! calculating other soil hydraulic properties
     374             :           ! tillage horizons properties depending on the LC class
     375           0 :           do L = 1, max_LCover
     376             :             select case (L)
     377             :             case(1)               ! forest
     378           0 :               pOM = tmp_orgMatterContent_forest
     379             :             case(2)               ! impervious
     380           0 :               pOM = tmp_orgMatterContent_impervious !param(2)
     381             :             case(3)               ! permeable
     382           0 :               pOM = tmp_orgMatterContent_pervious
     383             :             case default
     384           0 :                call error_message('Error mpr_sm: pOM used is not initialized.')
     385             :             end select
     386           0 :             pM = 100.0_dp - pOM
     387             :             ! bulk density acording to Rawl's (1982) paper
     388           0 :             Db(i, j, L) = 100.0_dp / ((pOM / BulkDens_OrgMatter) + (pM / DbM(i, j)))
     389             :             ! Effect of organic matter content on Ks estimates
     390             :             ! This is taken into account in a simplified form by using
     391             :             ! the ratio of (Bd/BdOM)
     392           0 :             Ks_tmp = Ks_tmp * (DbM(i, j) / Db(i, j, L))
     393           0 :             Ks(i, j, L) = Ks_tmp
     394             :             ! estimated SMs_till & van Genuchten's shape parameter (n)
     395           0 :             call Genuchten(thetaS_till(i, j, L), Genu_Mual_n, Genu_Mual_alpha, &
     396           0 :                     param(4 : 9), sand(i, j), clay(i, j), Db(i, j, L))
     397             :             ! estimating field capacity
     398           0 :             call field_cap(thetaFC_till(i, j, L), Ks_tmp, thetaS_till(i, j, L), Genu_Mual_n)
     399             :             ! estimating permanent wilting point
     400           0 :             call PWP(Genu_Mual_n, Genu_Mual_alpha, thetaS_till(i, j, L), thetaPW_till(i, j, L))
     401             :           end do
     402             : 
     403             :           ! *** FOR NON-TILLAGE TYPE OF SOILS ***
     404             :           ! note j = 1
     405             :           ! since Ks_tmp has changed earlier ... get the original Ks once again
     406           0 :           Ks_tmp = Ks_non_till(i, j)
     407             :           ! estimate SMs & van Genuchten's shape parameter (n)
     408           0 :           call Genuchten(thetaS(i, j), Genu_Mual_n, Genu_Mual_alpha, param(4 : 9), sand(i, j), clay(i, j), DbM(i, j))
     409             :           ! estimate field capacity
     410           0 :           call field_cap(thetaFC(i, j), Ks_tmp, thetaS(i, j), Genu_Mual_n)
     411             :           ! estimate permanent wilting point
     412           0 :           call PWP(Genu_Mual_n, Genu_Mual_alpha, thetaS(i, j), thetaPW(i, j))
     413             : 
     414             :         end do  ! >> HORIZON
     415             :       end do   ! >> SOIL TYPE
     416             : 
     417             :       ! calculate other soil properties at each location [L0] for regionalising model parameters
     418           0 :       do i = 1, size(soilId0, 1)     !! over all cells
     419           0 :         do j = 1, size(soilId0, 2)  !! over horizons
     420           0 :           s = soilId0(i, j)
     421           0 :           if (j .le. nTillHorizons(1)) then
     422             :             ! soil properties over the whole soil coloum depth
     423           0 :             KsVar_H0(i) = KsVar_H0(i) + thetaS_till (s, 1, LCover0(i)) * Ks(s, 1, LCover0(i))
     424           0 :             KsVar_V0(i) = KsVar_V0(i) + thetaS_till (s, 1, LCover0(i)) / Ks(s, 1, LCover0(i))
     425           0 :             SMs_FC0(i) = SMs_FC0 (i) + thetaFC_till(s, 1, LCover0(i))
     426           0 :             SMs_tot0(i) = SMs_tot0(i) + thetaS_till (s, 1, LCover0(i))
     427             :           else
     428             :             ! soil_properties over the whole soil column
     429           0 :             KsVar_H0(i) = KsVar_H0(i) + thetaS (s, 1) * Ks_non_till(s, 1)
     430           0 :             KsVar_V0(i) = KsVar_V0(i) + thetaS (s, 1) / Ks_non_till(s, 1)
     431           0 :             SMs_FC0(i) = SMs_FC0 (i) + thetaFC(s, 1)
     432           0 :             SMs_tot0(i) = SMs_tot0(i) + thetaS (s, 1)
     433             :           end if
     434             :         end do
     435             :         ! ------------------------------------------------------------------
     436             :         ! DETERMINE RELATIVE VARIABILITIES OF
     437             :         ! Ks FOR HORIZONTAL FLOW (KsVar_H) & Ks FOR VERTICAL FLOW (KsVar_V)
     438             :         ! ------------------------------------------------------------------
     439             :         ! soil moisture saturation deficit relative to the field capacity soil moisture
     440           0 :         SMs_FC0(i) = (SMs_tot0(i) - SMs_FC0(i)) / SMs_tot0(i)
     441             :         ! Ks variability over the whole soil coloum depth for
     442             :         ! both horizontal and vertical flows including relative variabilities
     443           0 :         KsVar_H0(i) = KsVar_H0(i) / SMs_tot0(i) / param(13)
     444           0 :         KsVar_V0(i) = SMs_tot0(i) / KsVar_V0(i) / param(13)
     445             :       end do
     446             : 
     447             :     CASE DEFAULT
     448         414 :       call error_message('***ERROR: iFlag_soilDB option given does not exist. Only 0 and 1 is taken at the moment.')
     449             :     END SELECT
     450             : 
     451             :     ! free space **
     452         414 :     deallocate(SMs_tot0)
     453         414 :     if(allocated(Ks_non_till)) deallocate(Ks_non_till)
     454             : 
     455         414 :   end subroutine mpr_sm
     456             : 
     457             :   ! ------------------------------------------------------------------
     458             : 
     459             :   !    NAME
     460             :   !        PWP
     461             : 
     462             :   !    PURPOSE
     463             :   !>       \brief Permanent Wilting point
     464             : 
     465             :   !>       \details This subroutine calculates the permanent wilting
     466             :   !>       point according to Zacharias et al. (2007, Soil Phy.) and
     467             :   !>       using van Genuchten 1980's equation. For the water retention curve at
     468             :   !>       a matrix potential of -1500 kPa, it is assumed that thetaR = 0.
     469             :   !>       ADDITIONAL INFORMATION
     470             :   !>       Zacharias et al. 2007, Soil Phy.
     471             : 
     472             :   !    INTENT(IN)
     473             :   !>       \param[in] "real(dp) :: Genu_Mual_n"     - Genuchten shape parameter
     474             :   !>       \param[in] "real(dp) :: Genu_Mual_alpha" - Genuchten shape parameter
     475             :   !>       \param[in] "real(dp) :: thetaS"          - saturated water content
     476             : 
     477             :   !    INTENT(OUT)
     478             :   !>       \param[out] "real(dp) :: thetaPWP" - Permanent Wilting point
     479             : 
     480             :   !    HISTORY
     481             :   !>       \authors Stephan Thober, Rohini Kumar
     482             : 
     483             :   !>       \date Dec, 2012
     484             : 
     485             :   ! Modifications:
     486             :   ! Matthias Zink Nov 2013 - documentation, moved constants to mhm_constants
     487             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     488             : 
     489       71130 :   elemental pure subroutine PWP(Genu_Mual_n, Genu_Mual_alpha, thetaS, thetaPWP)
     490             : 
     491         414 :     use mo_mpr_constants, only : PWP_c, PWP_matPot_ThetaR
     492             : 
     493             :     implicit none
     494             : 
     495             :     ! - Genuchten shape parameter
     496             :     real(dp), intent(in) :: Genu_Mual_n
     497             : 
     498             :     ! - Genuchten shape parameter
     499             :     real(dp), intent(in) :: Genu_Mual_alpha
     500             : 
     501             :     ! - saturated water content
     502             :     real(dp), intent(in) :: thetaS
     503             : 
     504             :     ! - Permanent Wilting point
     505             :     real(dp), intent(out) :: thetaPWP
     506             : 
     507       71130 :     real(dp) :: x
     508             : 
     509             :     ! Genuchten shape parameter
     510       71130 :     real(dp) :: Genu_Mual_m
     511             : 
     512             : 
     513       71130 :     Genu_Mual_m = PWP_c - (PWP_c / Genu_Mual_n)
     514       71130 :     x = PWP_c + exp(Genu_Mual_n * log(Genu_Mual_alpha * PWP_matPot_ThetaR))
     515       71130 :     x = exp(Genu_Mual_m * log(x))
     516             :     ! constrain
     517           0 :     if (x < 1.0_dp) x = 1.0_dp
     518       71130 :     thetaPWP = thetaS / x
     519             : 
     520       71130 :   end subroutine PWP
     521             : 
     522             :   ! ----------------------------------------------------------------------------
     523             : 
     524             :   !    NAME
     525             :   !        field_cap
     526             : 
     527             :   !    PURPOSE
     528             :   !>       \brief calculates the field capacity
     529             : 
     530             :   !>       \details estimate Field capacity; FC -- Flux based
     531             :   !>       approach (Twarakavi, et. al. 2009, WRR)
     532             :   !>       According to the
     533             :   !>       above reference FC is defined as the soil water content at
     534             :   !>       which the drainage from a profile ceases under natural
     535             :   !>       conditions. Since drainage from a soil profile in a simulation
     536             :   !>       never becomes zero, we assume that drainage ceases when the
     537             :   !>       bottom flux from the soil reaches a value that is equivalent to
     538             :   !>       the minimum amount of precipitation that could be recorded
     539             :   !>       (i.e. 0.01 cm/d == 1 mm/d). It is assumed that ThetaR = 0.0_dp
     540             :   !>       ADDITIONAL INFORMATION
     541             :   !>       Twarakavi, et. al. 2009, WRR
     542             : 
     543             :   !    INTENT(OUT)
     544             :   !>       \param[out] "real(dp) :: thetaFC" - Field capacity
     545             : 
     546             :   !    INTENT(IN)
     547             :   !>       \param[in] "real(dp) :: Ks"          - saturated hydraulic conductivity
     548             :   !>       \param[in] "real(dp) :: thetaS"      - saturated water content
     549             :   !>       \param[in] "real(dp) :: Genu_Mual_n" - Genuchten shape parameter
     550             : 
     551             :   !    HISTORY
     552             :   !>       \authors Stephan Thober, Rohini Kumar
     553             : 
     554             :   !>       \date Dec 2012
     555             : 
     556             :   ! Modifications:
     557             :   ! Matthias Zink Nov 2013 - documentation, moved constants to mhm_constants
     558             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     559             : 
     560       71130 :   elemental pure subroutine field_cap(thetaFC, Ks, thetaS, Genu_Mual_n)
     561             : 
     562       71130 :     use mo_mpr_constants, only : field_cap_c1, field_cap_c2
     563             : 
     564             :     implicit none
     565             : 
     566             :     ! - saturated hydraulic conductivity
     567             :     real(dp), intent(in) :: Ks
     568             : 
     569             :     ! - saturated water content
     570             :     real(dp), intent(in) :: thetaS
     571             : 
     572             :     ! - Genuchten shape parameter
     573             :     real(dp), intent(in) :: Genu_Mual_n
     574             : 
     575             :     ! - Field capacity
     576             :     real(dp), intent(out) :: thetaFC
     577             : 
     578       71130 :     real(dp) :: x
     579             : 
     580             : 
     581       71130 :     x = (field_cap_c1) * (field_cap_c2 + log10(Ks))
     582       71130 :     thetaFC = thetaS * exp(x * log(Genu_Mual_n))
     583             : 
     584       71130 :   end subroutine field_cap
     585             : 
     586             :   ! ----------------------------------------------------------------------------
     587             : 
     588             :   !    NAME
     589             :   !        Genuchten
     590             : 
     591             :   !    PURPOSE
     592             :   !>       \brief calculates the Genuchten shape parameter
     593             : 
     594             :   !>       \details estimate SMs_till & van Genuchten's shape parameter (n)
     595             :   !>       (Zacharias et al, 2007, soil Phy.)
     596             :   !>       Global parameters needed (see mhm_parameter.nml):
     597             :   !>       - param( 1) = PTF_lower66_5_constant
     598             :   !>       - param( 2) = PTF_lower66_5_clay
     599             :   !>       - param( 3) = PTF_lower66_5_Db
     600             :   !>       - param( 4) = PTF_higher66_5_constant
     601             :   !>       - param( 5) = PTF_higher66_5_clay
     602             :   !>       - param( 6) = PTF_higher66_5_Db
     603             :   !>       ADDITIONAL INFORMATION
     604             :   !>       Zacharias et al, 2007, soil Phy.
     605             : 
     606             :   !    INTENT(OUT)
     607             :   !>       \param[out] "real(dp) :: thetaS"          - saturated water content
     608             :   !>       \param[out] "real(dp) :: Genu_Mual_n"     - van Genuchten shape parameter
     609             :   !>       \param[out] "real(dp) :: Genu_Mual_alpha" - van Genuchten shape parameter
     610             : 
     611             :   !    INTENT(IN)
     612             :   !>       \param[in] "real(dp), dimension(6) :: param" parameters
     613             :   !>       \param[in] "real(dp) :: sand"                - [%] sand content
     614             :   !>       \param[in] "real(dp) :: clay"                - [%] clay content
     615             :   !>       \param[in] "real(dp) :: Db"                  - [10^3 kg/m3] bulk density
     616             : 
     617             :   !    HISTORY
     618             :   !>       \authors Stephan Thober, Rohini Kumar
     619             : 
     620             :   !>       \date Dec 2012
     621             : 
     622             :   ! Modifications:
     623             :   ! Rohini Kumar Mar 2014 - ThetaS limit changed from 0 to 0.001
     624             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     625             : 
     626       71130 :   subroutine Genuchten(thetaS, Genu_Mual_n, Genu_Mual_alpha, param, sand, clay, Db)
     627             : 
     628       71130 :     use mo_mpr_constants, only : vGenuchtenN_c1, vGenuchtenN_c10, vGenuchtenN_c11, vGenuchtenN_c12, vGenuchtenN_c13, &
     629             :                                  vGenuchtenN_c14, vGenuchtenN_c15, vGenuchtenN_c16, vGenuchtenN_c17, vGenuchtenN_c18, &
     630             :                                  vGenuchtenN_c2, vGenuchtenN_c3, vGenuchtenN_c4, vGenuchtenN_c5, vGenuchtenN_c6, &
     631             :                                  vGenuchtenN_c7, vGenuchtenN_c8, vGenuchtenN_c9, vGenuchten_sandtresh
     632             : 
     633             :     implicit none
     634             : 
     635             :     ! parameters
     636             :     real(dp), dimension(6), intent(in) :: param
     637             : 
     638             :     ! - [%] sand content
     639             :     real(dp), intent(in) :: sand
     640             : 
     641             :     ! - [%] clay content
     642             :     real(dp), intent(in) :: clay
     643             : 
     644             :     ! - [10^3 kg/m3] bulk density
     645             :     real(dp), intent(in) :: Db
     646             : 
     647             :     ! - saturated water content
     648             :     real(dp), intent(out) :: thetaS
     649             : 
     650             :     ! - van Genuchten shape parameter
     651             :     real(dp), intent(out) :: Genu_Mual_n
     652             : 
     653             :     ! - van Genuchten shape parameter
     654             :     real(dp), intent(out) :: Genu_Mual_alpha
     655             : 
     656             :     ! temporal variable
     657       71130 :     real(dp) :: x
     658             : 
     659             : 
     660             :     ! estimate SMs_till & van Genuchten's parameters (alpha and n)
     661       71130 :     if (sand < vGenuchten_sandtresh) then
     662       48360 :       thetaS = param(1) + param(2) * clay + param(3) * Db
     663             :       Genu_Mual_n = vGenuchtenN_c1 - vGenuchtenN_c2 * (sand**(vGenuchtenN_c3)) + &
     664       48360 :               vGenuchtenN_c4 * (clay**(vGenuchtenN_c5))
     665             :       x = vGenuchtenN_c6 + vGenuchtenN_c7 * sand + vGenuchtenN_c8 * clay - &
     666       48360 :               vGenuchtenN_c9 * Db
     667             :     else
     668       22770 :       thetaS = param(4) + param(5) * clay + param(6) * Db
     669             :       Genu_Mual_n = vGenuchtenN_c10 + vGenuchtenN_c11 * (sand**(vGenuchtenN_c12)) + &
     670       22770 :               vGenuchtenN_c13 * (clay**(vGenuchtenN_c14))
     671             :       x = vGenuchtenN_c15 + vGenuchtenN_c16 * sand + vGenuchtenN_c17 * clay - &
     672       22770 :               vGenuchtenN_c18 * Db
     673             :     end if
     674             : 
     675             :     ! Mualem alpha
     676       71130 :     Genu_Mual_alpha = exp(x)
     677             : 
     678             :     ! hard coded limits, according to (Zacharias et al, 2007, soil Phy.)
     679       71130 :     if (thetaS < 0.01_dp) then
     680           0 :       call message('thetaS below threshold limit 1e-2, reset.')
     681             :       ! Put constrains on theta_S
     682           0 :       thetaS = 0.01_dp
     683             :     end if
     684       71130 :     if (thetaS > 1.0_dp) then
     685           0 :       call message('thetaS above 1, reset.')
     686             :       ! Put constrains on theta_S
     687           0 :       thetaS = 1.0_dp
     688             :     end if
     689       71130 :     if (Genu_Mual_n < 1.01000_dp) then
     690           0 :       call message('Genu_Mual_n below threshold limit 1.01, reset.')
     691           0 :       Genu_Mual_n = 1.01000_dp
     692             :     end if
     693       71130 :     if (Genu_Mual_alpha < 0.00001_dp) then
     694           0 :       call message('Genu_Mual_alpha below threshold limit 1e-5, reset.')
     695           0 :       Genu_Mual_alpha = 0.00001_dp
     696             :     end if
     697             : 
     698       71130 :   end subroutine Genuchten
     699             : 
     700             :   ! ----------------------------------------------------------------------------
     701             : 
     702             :   !    NAME
     703             :   !        hydro_cond
     704             : 
     705             :   !    PURPOSE
     706             :   !>       \brief calculates the hydraulic conductivity Ks
     707             : 
     708             :   !>       \details By default save this value of Ks, particularly for the
     709             :   !>       deeper layers where OM content plays relatively low or no role
     710             :   !>       Global parameters needed (see mhm_parameter.nml):
     711             :   !>       - param(1) = PTF_Ks_constant
     712             :   !>       - param(2) = PTF_Ks_sand
     713             :   !>       - param(3) = PTF_Ks_clay
     714             :   !>       - param(4) = PTF_Ks_curveSlope
     715             :   !>       ADDITIONAL INFORMATION
     716             :   !>       Written,  Stephan Thober, Dec 2012
     717             : 
     718             :   !    INTENT(OUT)
     719             :   !>       \param[out] "real(dp) :: KS"
     720             : 
     721             :   !    INTENT(IN)
     722             :   !>       \param[in] "real(dp), dimension(4) :: param"
     723             :   !>       \param[in] "real(dp) :: sand"                - [%] sand content
     724             :   !>       \param[in] "real(dp) :: clay"                - [%] clay content
     725             : 
     726             :   !    HISTORY
     727             :   !>       \authors Stephan Thober, Rohini Kumar
     728             : 
     729             :   !>       \date Dec 2012
     730             : 
     731             :   ! Modifications:
     732             :   ! Matthias Zink  Nov 2013 - documentation, moved constants to mhm_constants
     733             :   ! Matthias Cuntz Jun 2014 - suggested to fix param(4)
     734             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     735             : 
     736       42678 :   subroutine hydro_cond(KS, param, sand, clay)
     737             : 
     738       71130 :     use mo_mpr_constants, only : Ks_c
     739             : 
     740             :     implicit none
     741             : 
     742             :     real(dp), dimension(4), intent(in) :: param
     743             : 
     744             :     ! - [%] sand content
     745             :     real(dp), intent(in) :: sand
     746             : 
     747             :     ! - [%] clay content
     748             :     real(dp), intent(in) :: clay
     749             : 
     750             :     real(dp), intent(out) :: KS
     751             : 
     752             :     ! temporal variable
     753       42678 :     real(dp) :: x
     754             : 
     755             : 
     756             :     ! saturated vertical hydraulic conductivity, Ks (cm/d)
     757             :     !   from Cosby et. al. (WRR 1984) Table 4
     758             :     ! param(4) is the unit conversion from inch/h to cm/d and should be a constant.
     759             :     ! Fix it in the namelist, i.e. in
     760             :     ! mhm_parameter.nml set the 4th value (=FLAG) to 0 and the third value to 60.96
     761             :     !   PTF_Ks_curveSlope = 60.96, 60.96, 60.96, 0, 1
     762       42678 :     x = param(1) + param(2) * sand - param(3) * clay
     763       42678 :     Ks = param(4) * exp(X * log(Ks_c))
     764             : 
     765       42678 :     if (Ks < 1.10_dp) then
     766           0 :       call message('JMJMJM-Ks-BAD')
     767             :     end if
     768             : 
     769             :     ! minimum value of Ks = 1.1cm/d
     770       42678 :     if (Ks < 1.10_dp) Ks = 1.10_dp
     771             : 
     772       42678 :   end subroutine hydro_cond
     773             : 
     774             : end module mo_mpr_soilmoist

Generated by: LCOV version 1.16