#include <misc.h>
#include <preproc.h>
module SoilHydrologyMod 1
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: SoilHydrologyMod
!
! !DESCRIPTION:
! Calculate soil hydrology
!
! !PUBLIC TYPES:
implicit none
save
!
! !PUBLIC MEMBER FUNCTIONS:
public :: SurfaceRunoff ! Calculate surface runoff
public :: Infiltration ! Calculate infiltration into surface soil layer
public :: SoilWater ! Calculate soil hydrology
public :: Drainage ! Calculate subsurface drainage
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
! 04/25/07 Keith Oleson: CLM3.5 hydrology
!
!EOP
!-----------------------------------------------------------------------
contains
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: SurfaceRunoff
!
! !INTERFACE:
subroutine SurfaceRunoff (lbc, ubc, lbp, ubp, num_hydrologyc, filter_hydrologyc, & 1,6
num_urbanc, filter_urbanc, vol_liq, icefrac)
!
! !DESCRIPTION:
! Calculate surface runoff
!
! !USES:
use shr_kind_mod , only : r8 => shr_kind_r8
use clmtype
use clm_varcon , only : denice, denh2o, wimp, pondmx_urban, &
icol_roof, icol_sunwall, icol_shadewall, &
icol_road_imperv, icol_road_perv
use clm_varpar , only : nlevsoi, maxpatch_pft
use clm_time_manager, only : get_step_size
!
! !ARGUMENTS:
implicit none
integer , intent(in) :: lbc, ubc ! column bounds
integer , intent(in) :: lbp, ubp ! pft bounds
integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter
integer , intent(in) :: filter_hydrologyc(ubc-lbc+1) ! column filter for soil points
integer , intent(in) :: num_urbanc ! number of column urban points in column filter
integer , intent(in) :: filter_urbanc(ubc-lbc+1) ! column filter for urban points
real(r8), intent(out) :: vol_liq(lbc:ubc,1:nlevsoi) ! partial volume of liquid water in layer
real(r8), intent(out) :: icefrac(lbc:ubc,1:nlevsoi) ! fraction of ice in layer (-)
!
! !CALLED FROM:
! subroutine Hydrology2 in module Hydrology2Mod
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 12 November 1999: Z.-L. Yang and G.-Y. Niu
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
! 2/26/02, Peter Thornton: Migrated to new data structures.
! 4/26/05, David Lawrence: Made surface runoff for dry soils a function
! of rooting fraction in top three soil layers.
! 04/25/07 Keith Oleson: Completely new routine for CLM3.5 hydrology
!
! !LOCAL VARIABLES:
!
! local pointers to original implicit in arguments
!
integer , pointer :: cgridcell(:) ! gridcell index for each column
integer , pointer :: ctype(:) ! column type index
real(r8), pointer :: qflx_top_soil(:) !net water input into soil from top (mm/s)
real(r8), pointer :: watsat(:,:) !volumetric soil water at saturation (porosity)
real(r8), pointer :: hkdepth(:) !decay factor (m)
real(r8), pointer :: zwt(:) !water table depth (m)
real(r8), pointer :: fcov(:) !fractional impermeable area
real(r8), pointer :: fsat(:) !fractional area with water table at surface
real(r8), pointer :: dz(:,:) !layer depth (m)
real(r8), pointer :: h2osoi_ice(:,:) !ice lens (kg/m2)
real(r8), pointer :: h2osoi_liq(:,:) !liquid water (kg/m2)
real(r8), pointer :: wtfact(:) !maximum saturated fraction for a gridcell
real(r8), pointer :: hksat(:,:) ! hydraulic conductivity at saturation (mm H2O /s)
real(r8), pointer :: bsw(:,:) ! Clapp and Hornberger "b"
real(r8), pointer :: sucsat(:,:) ! minimum soil suction (mm)
integer , pointer :: snl(:) ! minus number of snow layers
real(r8), pointer :: qflx_evap_grnd(:) ! ground surface evaporation rate (mm H2O/s) [+]
real(r8), pointer :: zi(:,:) ! interface level below a "z" level (m)
!
! local pointers to original implicit out arguments
!
real(r8), pointer :: qflx_surf(:) ! surface runoff (mm H2O /s)
real(r8), pointer :: eff_porosity(:,:) ! effective porosity = porosity - vol_ice
real(r8), pointer :: fracice(:,:) !fractional impermeability (-)
!
!EOP
!
! !OTHER LOCAL VARIABLES:
!
integer :: c,j,fc,g !indices
real(r8) :: dtime ! land model time step (sec)
real(r8) :: xs(lbc:ubc) ! excess soil water above urban ponding limit
real(r8) :: vol_ice(lbc:ubc,1:nlevsoi) !partial volume of ice lens in layer
real(r8) :: fff(lbc:ubc) !decay factor (m-1)
real(r8) :: s1 !variable to calculate qinmax
real(r8) :: su !variable to calculate qinmax
real(r8) :: v !variable to calculate qinmax
real(r8) :: qinmax !maximum infiltration capacity (mm/s)
!-----------------------------------------------------------------------
! Assign local pointers to derived subtype components (column-level)
ctype => clm3%g%l%c%itype
qflx_top_soil => clm3%g%l%c%cwf%qflx_top_soil
qflx_surf => clm3%g%l%c%cwf%qflx_surf
watsat => clm3%g%l%c%cps%watsat
hkdepth => clm3%g%l%c%cps%hkdepth
dz => clm3%g%l%c%cps%dz
h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice
h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq
fcov => clm3%g%l%c%cws%fcov
fsat => clm3%g%l%c%cws%fsat
eff_porosity => clm3%g%l%c%cps%eff_porosity
wtfact => clm3%g%l%c%cps%wtfact
zwt => clm3%g%l%c%cws%zwt
fracice => clm3%g%l%c%cps%fracice
hksat => clm3%g%l%c%cps%hksat
bsw => clm3%g%l%c%cps%bsw
sucsat => clm3%g%l%c%cps%sucsat
snl => clm3%g%l%c%cps%snl
qflx_evap_grnd => clm3%g%l%c%cwf%pwf_a%qflx_evap_grnd
zi => clm3%g%l%c%cps%zi
! Get time step
dtime = get_step_size()
do j = 1,nlevsoi
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
! Porosity of soil, partial volume of ice and liquid, fraction of ice in each layer,
! fractional impermeability
vol_ice(c,j) = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice))
eff_porosity(c,j) = max(0.01_r8,watsat(c,j)-vol_ice(c,j))
vol_liq(c,j) = min(eff_porosity(c,j), h2osoi_liq(c,j)/(dz(c,j)*denh2o))
icefrac(c,j) = min(1._r8,h2osoi_ice(c,j)/(h2osoi_ice(c,j)+h2osoi_liq(c,j)))
fracice(c,j) = max(0._r8,exp(-3._r8*(1._r8-icefrac(c,j)))- exp(-3._r8))/(1.0_r8-exp(-3._r8))
end do
end do
! Saturated fraction
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
fff(c) = 0.5_r8
fsat(c) = wtfact(c) * exp(-0.5_r8*fff(c)*zwt(c))
fcov(c) = (1._r8 - fracice(c,1)) * fsat(c) + fracice(c,1)
end do
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
! Maximum infiltration capacity
s1 = max(0.01_r8,vol_liq(c,1)/max(wimp,eff_porosity(c,1)))
su = max(0._r8,(s1-fcov(c)) / (max(0.01_r8,1._r8-fcov(c))))
v = -bsw(c,1)*sucsat(c,1)/(0.5_r8*dz(c,1)*1000._r8)
qinmax = (1._r8+v*(su-1._r8))*hksat(c,1)
! Surface runoff
qflx_surf(c) = fcov(c) * qflx_top_soil(c) + &
(1._r8-fcov(c)) * max(0._r8, qflx_top_soil(c)-qinmax)
end do
! Determine water in excess of ponding limit for urban roof and impervious road.
! Excess goes to surface runoff. No surface runoff for sunwall and shadewall.
!dir$ concurrent
!cdir nodep
do fc = 1, num_urbanc
c = filter_urbanc(fc)
if (ctype(c) == icol_roof .or. ctype(c) == icol_road_imperv) then
! If there are snow layers then all qflx_top_soil goes to surface runoff
if (snl(c) < 0) then
qflx_surf(c) = max(0._r8,qflx_top_soil(c))
else
xs(c) = max(0._r8, &
h2osoi_liq(c,1)/dtime + qflx_top_soil(c) - qflx_evap_grnd(c) - &
pondmx_urban/dtime)
if (xs(c) > 0.) then
h2osoi_liq(c,1) = pondmx_urban
else
h2osoi_liq(c,1) = max(0._r8,h2osoi_liq(c,1)+ &
(qflx_top_soil(c)-qflx_evap_grnd(c))*dtime)
end if
qflx_surf(c) = xs(c)
end if
else if (ctype(c) == icol_sunwall .or. ctype(c) == icol_shadewall) then
qflx_surf(c) = 0._r8
end if
end do
end subroutine SurfaceRunoff
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Infiltration
!
! !INTERFACE:
subroutine Infiltration(lbc, ubc, num_hydrologyc, filter_hydrologyc, & 1,3
num_urbanc, filter_urbanc)
!
! !DESCRIPTION:
! Calculate infiltration into surface soil layer (minus the evaporation)
!
! !USES:
use shr_kind_mod, only : r8 => shr_kind_r8
use clm_varcon , only : icol_roof, icol_road_imperv, icol_sunwall, icol_shadewall, &
icol_road_perv
use clmtype
!
! !ARGUMENTS:
implicit none
integer, intent(in) :: lbc, ubc ! column bounds
integer, intent(in) :: num_hydrologyc ! number of column soil points in column filter
integer, intent(in) :: filter_hydrologyc(ubc-lbc+1) ! column filter for soil points
integer, intent(in) :: num_urbanc ! number of column urban points in column filter
integer, intent(in) :: filter_urbanc(ubc-lbc+1) ! column filter for urban points
!
! !CALLED FROM:
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 12 November 1999: Z.-L. Yang and G.-Y. Niu
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
! 2/27/02, Peter Thornton: Migrated to new data structures.
!
! !LOCAL VARIABLES:
!
! local pointers to original implicit in arguments
!
integer , pointer :: ctype(:) ! column type index
integer , pointer :: snl(:) ! minus number of snow layers
real(r8), pointer :: qflx_top_soil(:) ! net water input into soil from top (mm/s)
real(r8), pointer :: qflx_surf(:) ! surface runoff (mm H2O /s)
real(r8), pointer :: qflx_evap_grnd(:)! ground surface evaporation rate (mm H2O/s) [+]
!
! local pointers to original implicit out arguments
!
real(r8), pointer :: qflx_infl(:) !infiltration (mm H2O /s)
!
!EOP
!
! !OTHER LOCAL VARIABLES:
!
integer :: c, fc !indices
!-----------------------------------------------------------------------
! Assign local pointers to derived type members (column-level)
ctype => clm3%g%l%c%itype
snl => clm3%g%l%c%cps%snl
qflx_top_soil => clm3%g%l%c%cwf%qflx_top_soil
qflx_surf => clm3%g%l%c%cwf%qflx_surf
qflx_infl => clm3%g%l%c%cwf%qflx_infl
qflx_evap_grnd => clm3%g%l%c%cwf%pwf_a%qflx_evap_grnd
! Infiltration into surface soil layer (minus the evaporation)
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
if (snl(c) >= 0) then
qflx_infl(c) = qflx_top_soil(c) - qflx_surf(c) - qflx_evap_grnd(c)
else
qflx_infl(c) = qflx_top_soil(c) - qflx_surf(c)
end if
end do
! No infiltration for impervious urban surfaces
!dir$ concurrent
!cdir nodep
do fc = 1, num_urbanc
c = filter_urbanc(fc)
if (ctype(c) /= icol_road_perv) then
qflx_infl(c) = 0._r8
end if
end do
end subroutine Infiltration
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: SoilWater
!
! !INTERFACE:
subroutine SoilWater(lbc, ubc, num_hydrologyc, filter_hydrologyc, & 1,10
num_urbanc, filter_urbanc, &
vol_liq, dwat, hk, dhkdw)
!
! !DESCRIPTION:
! Soil hydrology
! Soil moisture is predicted from a 10-layer model (as with soil
! temperature), in which the vertical soil moisture transport is governed
! by infiltration, runoff, gradient diffusion, gravity, and root
! extraction through canopy transpiration. The net water applied to the
! surface layer is the snowmelt plus precipitation plus the throughfall
! of canopy dew minus surface runoff and evaporation.
! CLM3.5 uses a zero-flow bottom boundary condition.
!
! The vertical water flow in an unsaturated porous media is described by
! Darcy's law, and the hydraulic conductivity and the soil negative
! potential vary with soil water content and soil texture based on the work
! of Clapp and Hornberger (1978) and Cosby et al. (1984). The equation is
! integrated over the layer thickness, in which the time rate of change in
! water mass must equal the net flow across the bounding interface, plus the
! rate of internal source or sink. The terms of water flow across the layer
! interfaces are linearly expanded by using first-order Taylor expansion.
! The equations result in a tridiagonal system equation.
!
! Note: length units here are all millimeter
! (in temperature subroutine uses same soil layer
! structure required but lengths are m)
!
! Richards equation:
!
! d wat d d wat d psi
! ----- = - -- [ k(----- ----- - 1) ] + S
! dt dz dz d wat
!
! where: wat = volume of water per volume of soil (mm**3/mm**3)
! psi = soil matrix potential (mm)
! dt = time step (s)
! z = depth (mm)
! dz = thickness (mm)
! qin = inflow at top (mm h2o /s)
! qout= outflow at bottom (mm h2o /s)
! s = source/sink flux (mm h2o /s)
! k = hydraulic conductivity (mm h2o /s)
!
! d qin d qin
! qin[n+1] = qin[n] + -------- d wat(j-1) + --------- d wat(j)
! d wat(j-1) d wat(j)
! ==================|=================
! < qin
!
! d wat(j)/dt * dz = qin[n+1] - qout[n+1] + S(j)
!
! > qout
! ==================|=================
! d qout d qout
! qout[n+1] = qout[n] + --------- d wat(j) + --------- d wat(j+1)
! d wat(j) d wat(j+1)
!
!
! Solution: linearize k and psi about d wat and use tridiagonal
! system of equations to solve for d wat,
! where for layer j
!
!
! r_j = a_j [d wat_j-1] + b_j [d wat_j] + c_j [d wat_j+1]
!
! !USES:
use shr_kind_mod, only: r8 => shr_kind_r8
use clmtype
use clm_varcon , only : wimp, icol_roof, icol_road_imperv
use clm_varpar , only : nlevsoi, max_pft_per_col
use clm_varctl , only : iulog
use shr_const_mod , only : SHR_CONST_TKFRZ, SHR_CONST_LATICE, SHR_CONST_G
use TridiagonalMod, only : Tridiagonal
use clm_time_manager , only : get_step_size
!
! !ARGUMENTS:
implicit none
integer , intent(in) :: lbc, ubc ! column bounds
integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter
integer , intent(in) :: filter_hydrologyc(ubc-lbc+1) ! column filter for soil points
integer , intent(in) :: num_urbanc ! number of column urban points in column filter
integer , intent(in) :: filter_urbanc(ubc-lbc+1) ! column filter for urban points
real(r8), intent(in) :: vol_liq(lbc:ubc,1:nlevsoi) ! soil water per unit volume [mm/mm]
real(r8), intent(out) :: dwat(lbc:ubc,1:nlevsoi) ! change of soil water [m3/m3]
real(r8), intent(out) :: hk(lbc:ubc,1:nlevsoi) ! hydraulic conductivity [mm h2o/s]
real(r8), intent(out) :: dhkdw(lbc:ubc,1:nlevsoi) ! d(hk)/d(vol_liq)
!
! !CALLED FROM:
! subroutine Hydrology2 in module Hydrology2Mod
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
! 2/27/02, Peter Thornton: Migrated to new data structures. Includes
! treatment of multiple PFTs on a single soil column.
! 04/25/07 Keith Oleson: CLM3.5 hydrology
!
! !LOCAL VARIABLES:
!
! local pointers to original implicit in arguments
!
integer , pointer :: ctype(:) ! column type index
integer , pointer :: npfts(:) ! column's number of pfts - ADD
real(r8), pointer :: pwtcol(:) ! weight relative to column for each pft
real(r8), pointer :: pwtgcell(:) ! weight relative to gridcell for each pft
real(r8), pointer :: z(:,:) ! layer depth (m)
real(r8), pointer :: dz(:,:) ! layer thickness (m)
real(r8), pointer :: smpmin(:) ! restriction for min of soil potential (mm)
real(r8), pointer :: qflx_infl(:) ! infiltration (mm H2O /s)
real(r8), pointer :: qflx_tran_veg_pft(:) ! vegetation transpiration (mm H2O/s) (+ = to atm)
real(r8), pointer :: qflx_tran_veg_col(:) ! vegetation transpiration (mm H2O/s) (+ = to atm)
real(r8), pointer :: eff_porosity(:,:) ! effective porosity = porosity - vol_ice
real(r8), pointer :: watsat(:,:) ! volumetric soil water at saturation (porosity)
real(r8), pointer :: hksat(:,:) ! hydraulic conductivity at saturation (mm H2O /s)
real(r8), pointer :: bsw(:,:) ! Clapp and Hornberger "b"
real(r8), pointer :: sucsat(:,:) ! minimum soil suction (mm)
real(r8), pointer :: t_soisno(:,:) ! soil temperature (Kelvin)
real(r8), pointer :: rootr_pft(:,:) ! effective fraction of roots in each soil layer
integer , pointer :: pfti(:) ! beginning pft index for each column
real(r8), pointer :: fracice(:,:) ! fractional impermeability (-)
real(r8), pointer :: h2osoi_vol(:,:) ! volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3]
real(r8), pointer :: qcharge(:) ! aquifer recharge rate (mm/s)
real(r8), pointer :: hkdepth(:) ! decay factor (m)
real(r8), pointer :: zwt(:) ! water table depth (m)
real(r8), pointer :: zi(:,:) ! interface level below a "z" level (m)
!
! local pointers to original implicit inout arguments
!
real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2)
!
! local pointer s to original implicit out arguments
!
real(r8), pointer :: rootr_col(:,:) ! effective fraction of roots in each soil layer
real(r8), pointer :: smp_l(:,:) ! soil matrix potential [mm]
real(r8), pointer :: hk_l(:,:) ! hydraulic conductivity (mm/s)
!
!EOP
!
! !OTHER LOCAL VARIABLES:
!
integer :: p,c,fc,j ! do loop indices
integer :: jtop(lbc:ubc) ! top level at each column
real(r8) :: dtime ! land model time step (sec)
real(r8) :: amx(lbc:ubc,1:nlevsoi+1) ! "a" left off diagonal of tridiagonal matrix
real(r8) :: bmx(lbc:ubc,1:nlevsoi+1) ! "b" diagonal column for tridiagonal matrix
real(r8) :: cmx(lbc:ubc,1:nlevsoi+1) ! "c" right off diagonal tridiagonal matrix
real(r8) :: rmx(lbc:ubc,1:nlevsoi+1) ! "r" forcing term of tridiagonal matrix
real(r8) :: zmm(lbc:ubc,1:nlevsoi+1) ! layer depth [mm]
real(r8) :: dzmm(lbc:ubc,1:nlevsoi+1) ! layer thickness [mm]
real(r8) :: den ! used in calculating qin, qout
real(r8) :: dqidw0(lbc:ubc,1:nlevsoi+1) ! d(qin)/d(vol_liq(i-1))
real(r8) :: dqidw1(lbc:ubc,1:nlevsoi+1) ! d(qin)/d(vol_liq(i))
real(r8) :: dqodw1(lbc:ubc,1:nlevsoi+1) ! d(qout)/d(vol_liq(i))
real(r8) :: dqodw2(lbc:ubc,1:nlevsoi+1) ! d(qout)/d(vol_liq(i+1))
real(r8) :: dsmpdw(lbc:ubc,1:nlevsoi+1) ! d(smp)/d(vol_liq)
real(r8) :: num ! used in calculating qin, qout
real(r8) :: qin(lbc:ubc,1:nlevsoi+1) ! flux of water into soil layer [mm h2o/s]
real(r8) :: qout(lbc:ubc,1:nlevsoi+1) ! flux of water out of soil layer [mm h2o/s]
real(r8) :: s_node ! soil wetness
real(r8) :: s1 ! "s" at interface of layer
real(r8) :: s2 ! k*s**(2b+2)
real(r8) :: smp(lbc:ubc,1:nlevsoi) ! soil matrix potential [mm]
real(r8) :: sdamp ! extrapolates soiwat dependence of evaporation
integer :: pi ! pft index
real(r8) :: temp(lbc:ubc) ! accumulator for rootr weighting
integer :: jwt(lbc:ubc) ! index of the soil layer right above the water table (-)
real(r8) :: smp1,dsmpdw1,wh,wh_zwt,ka
real(r8) :: dwat2(lbc:ubc,1:nlevsoi+1)
real(r8) :: dzq ! used in calculating qin, qout (difference in equilbirium matric potential)
real(r8) :: zimm(lbc:ubc,0:nlevsoi) ! layer interface depth [mm]
real(r8) :: zq(lbc:ubc,1:nlevsoi+1) ! equilibrium matric potential for each layer [mm]
real(r8) :: vol_eq(lbc:ubc,1:nlevsoi+1) ! equilibrium volumetric water content
real(r8) :: tempi ! temp variable for calculating vol_eq
real(r8) :: temp0 ! temp variable for calculating vol_eq
real(r8) :: voleq1 ! temp variable for calculating vol_eq
real(r8) :: zwtmm(lbc:ubc) ! water table depth [mm]
!-----------------------------------------------------------------------
! Assign local pointers to derived type members (column-level)
qcharge => clm3%g%l%c%cws%qcharge
hkdepth => clm3%g%l%c%cps%hkdepth
zi => clm3%g%l%c%cps%zi
zwt => clm3%g%l%c%cws%zwt
ctype => clm3%g%l%c%itype
npfts => clm3%g%l%c%npfts
z => clm3%g%l%c%cps%z
dz => clm3%g%l%c%cps%dz
smpmin => clm3%g%l%c%cps%smpmin
watsat => clm3%g%l%c%cps%watsat
hksat => clm3%g%l%c%cps%hksat
bsw => clm3%g%l%c%cps%bsw
sucsat => clm3%g%l%c%cps%sucsat
eff_porosity => clm3%g%l%c%cps%eff_porosity
rootr_col => clm3%g%l%c%cps%rootr_column
t_soisno => clm3%g%l%c%ces%t_soisno
h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq
h2osoi_vol => clm3%g%l%c%cws%h2osoi_vol
qflx_infl => clm3%g%l%c%cwf%qflx_infl
fracice => clm3%g%l%c%cps%fracice
qflx_tran_veg_col => clm3%g%l%c%cwf%pwf_a%qflx_tran_veg
pfti => clm3%g%l%c%pfti
smp_l => clm3%g%l%c%cws%smp_l
hk_l => clm3%g%l%c%cws%hk_l
! Assign local pointers to derived type members (pft-level)
qflx_tran_veg_pft => clm3%g%l%c%p%pwf%qflx_tran_veg
rootr_pft => clm3%g%l%c%p%pps%rootr
pwtcol => clm3%g%l%c%p%wtcol
pwtgcell => clm3%g%l%c%p%wtgcell
! Get time step
dtime = get_step_size()
! Because the depths in this routine are in mm, use local
! variable arrays instead of pointers
do j = 1, nlevsoi
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
zmm(c,j) = z(c,j)*1.e3_r8
dzmm(c,j) = dz(c,j)*1.e3_r8
zimm(c,j) = zi(c,j)*1.e3_r8
end do
end do
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
zimm(c,0) = 0.0_r8
zwtmm(c) = zwt(c)*1.e3_r8
end do
! First step is to calculate the column-level effective rooting
! fraction in each soil layer. This is done outside the usual
! PFT-to-column averaging routines because it is not a simple
! weighted average of the PFT level rootr arrays. Instead, the
! weighting depends on both the per-unit-area transpiration
! of the PFT and the PFTs area relative to all PFTs.
temp(:) = 0._r8
do j = 1, nlevsoi
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
rootr_col(c,j) = 0._r8
end do
end do
do pi = 1,max_pft_per_col
do j = 1,nlevsoi
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
if (pi <= npfts(c)) then
p = pfti(c) + pi - 1
if (pwtgcell(p)>0._r8) then
rootr_col(c,j) = rootr_col(c,j) + rootr_pft(p,j) * qflx_tran_veg_pft(p) * pwtcol(p)
end if
end if
end do
end do
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
if (pi <= npfts(c)) then
p = pfti(c) + pi - 1
if (pwtgcell(p)>0._r8) then
temp(c) = temp(c) + qflx_tran_veg_pft(p) * pwtcol(p)
end if
end if
end do
end do
do j = 1, nlevsoi
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
if (temp(c) /= 0._r8) then
rootr_col(c,j) = rootr_col(c,j)/temp(c)
end if
end do
end do
!compute jwt index
! The layer index of the first unsaturated layer, i.e., the layer right above
! the water table
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
jwt(c) = nlevsoi
do j = 2,nlevsoi
if(zwt(c) <= zi(c,j)) then
jwt(c) = j-1
exit
end if
enddo
end do
! calculate the equilibrium water content based on the water table depth
do j=1,nlevsoi
do fc=1, num_hydrologyc
c = filter_hydrologyc(fc)
if ((zwtmm(c) .lt. zimm(c,j-1))) then !fully saturated when wtd is less than the layer top
vol_eq(c,j) = watsat(c,j)
! use the weighted average from the saturated part (depth > wtd) and the equilibrium solution for the
! rest of the layer
else if ((zwtmm(c) .lt. zimm(c,j)) .and. (zwtmm(c) .gt. zimm(c,j-1))) then
tempi = 1.0_r8
temp0 = (((sucsat(c,j)+zwtmm(c)-zimm(c,j-1))/sucsat(c,j)))**(1._r8-1._r8/bsw(c,j))
voleq1 = -sucsat(c,j)*watsat(c,j)/(1._r8-1._r8/bsw(c,j))/(zwtmm(c)-zimm(c,j-1))*(tempi-temp0)
vol_eq(c,j) = (voleq1*(zwtmm(c)-zimm(c,j-1)) + watsat(c,j)*(zimm(c,j)-zwtmm(c)))/(zimm(c,j)-zimm(c,j-1))
vol_eq(c,j) = min(watsat(c,j),vol_eq(c,j))
vol_eq(c,j) = max(vol_eq(c,j),0.0_r8)
else
tempi = (((sucsat(c,j)+zwtmm(c)-zimm(c,j))/sucsat(c,j)))**(1._r8-1._r8/bsw(c,j))
temp0 = (((sucsat(c,j)+zwtmm(c)-zimm(c,j-1))/sucsat(c,j)))**(1._r8-1._r8/bsw(c,j))
vol_eq(c,j) = -sucsat(c,j)*watsat(c,j)/(1._r8-1._r8/bsw(c,j))/(zimm(c,j)-zimm(c,j-1))*(tempi-temp0)
vol_eq(c,j) = max(vol_eq(c,j),0.0_r8)
vol_eq(c,j) = min(watsat(c,j),vol_eq(c,j))
endif
zq(c,j) = -sucsat(c,j)*(max(vol_eq(c,j)/watsat(c,j),0.01_r8))**(-bsw(c,j))
zq(c,j) = max(smpmin(c), zq(c,j))
end do
end do
! If water table is below soil column calculate zq for the 11th layer
j = nlevsoi
do fc=1, num_hydrologyc
c = filter_hydrologyc(fc)
if(jwt(c) == nlevsoi) then
tempi = 1._r8
temp0 = (((sucsat(c,j)+zwtmm(c)-zimm(c,j))/sucsat(c,j)))**(1._r8-1._r8/bsw(c,j))
vol_eq(c,j+1) = -sucsat(c,j)*watsat(c,j)/(1._r8-1._r8/bsw(c,j))/(zwtmm(c)-zimm(c,j))*(tempi-temp0)
vol_eq(c,j+1) = max(vol_eq(c,j+1),0.0_r8)
vol_eq(c,j+1) = min(watsat(c,j),vol_eq(c,j+1))
zq(c,j+1) = -sucsat(c,j)*(max(vol_eq(c,j+1)/watsat(c,j),0.01_r8))**(-bsw(c,j))
zq(c,j+1) = max(smpmin(c), zq(c,j+1))
end if
end do
! Hydraulic conductivity and soil matric potential and their derivatives
sdamp = 0._r8
do j = 1, nlevsoi
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
s1 = 0.5_r8*(h2osoi_vol(c,j) + h2osoi_vol(c,min(nlevsoi, j+1))) / &
(0.5_r8*(watsat(c,j)+watsat(c,min(nlevsoi, j+1))))
s1 = min(1._r8, s1)
s2 = hksat(c,j)*s1**(2._r8*bsw(c,j)+2._r8)
hk(c,j) = (1._r8-0.5_r8*(fracice(c,j)+fracice(c,min(nlevsoi, j+1))))*s1*s2
dhkdw(c,j) = (1._r8-0.5_r8*(fracice(c,j)+fracice(c,min(nlevsoi, j+1))))* &
(2._r8*bsw(c,j)+3._r8)*s2*0.5_r8/watsat(c,j)
s_node = max(h2osoi_vol(c,j)/watsat(c,j), 0.01_r8)
s_node = min(1.0_r8, s_node)
smp(c,j) = -sucsat(c,j)*s_node**(-bsw(c,j))
smp(c,j) = max(smpmin(c), smp(c,j))
dsmpdw(c,j) = -bsw(c,j)*smp(c,j)/(s_node*watsat(c,j))
smp_l(c,j) = smp(c,j)
hk_l(c,j) = hk(c,j)
end do
end do
! aquifer (11th) layer
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
zmm(c,nlevsoi+1) = 0.5*(1.e3_r8*zwt(c) + zmm(c,nlevsoi))
if(jwt(c) < nlevsoi) then
dzmm(c,nlevsoi+1) = dzmm(c,nlevsoi)
else
dzmm(c,nlevsoi+1) = (1.e3_r8*zwt(c) - zmm(c,nlevsoi))
end if
end do
! Set up r, a, b, and c vectors for tridiagonal solution
! Node j=1 (top)
j = 1
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
qin(c,j) = qflx_infl(c)
den = (zmm(c,j+1)-zmm(c,j))
dzq = (zq(c,j+1)-zq(c,j))
num = (smp(c,j+1)-smp(c,j)) - dzq
qout(c,j) = -hk(c,j)*num/den
dqodw1(c,j) = -(-hk(c,j)*dsmpdw(c,j) + num*dhkdw(c,j))/den
dqodw2(c,j) = -( hk(c,j)*dsmpdw(c,j+1) + num*dhkdw(c,j))/den
rmx(c,j) = qin(c,j) - qout(c,j) - qflx_tran_veg_col(c) * rootr_col(c,j)
amx(c,j) = 0._r8
bmx(c,j) = dzmm(c,j)*(sdamp+1._r8/dtime) + dqodw1(c,j)
cmx(c,j) = dqodw2(c,j)
end do
! Nodes j=2 to j=nlevsoi-1
do j = 2, nlevsoi - 1
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
den = (zmm(c,j) - zmm(c,j-1))
dzq = (zq(c,j)-zq(c,j-1))
num = (smp(c,j)-smp(c,j-1)) - dzq
qin(c,j) = -hk(c,j-1)*num/den
dqidw0(c,j) = -(-hk(c,j-1)*dsmpdw(c,j-1) + num*dhkdw(c,j-1))/den
dqidw1(c,j) = -( hk(c,j-1)*dsmpdw(c,j) + num*dhkdw(c,j-1))/den
den = (zmm(c,j+1)-zmm(c,j))
dzq = (zq(c,j+1)-zq(c,j))
num = (smp(c,j+1)-smp(c,j)) - dzq
qout(c,j) = -hk(c,j)*num/den
dqodw1(c,j) = -(-hk(c,j)*dsmpdw(c,j) + num*dhkdw(c,j))/den
dqodw2(c,j) = -( hk(c,j)*dsmpdw(c,j+1) + num*dhkdw(c,j))/den
rmx(c,j) = qin(c,j) - qout(c,j) - qflx_tran_veg_col(c)*rootr_col(c,j)
amx(c,j) = -dqidw0(c,j)
bmx(c,j) = dzmm(c,j)/dtime - dqidw1(c,j) + dqodw1(c,j)
cmx(c,j) = dqodw2(c,j)
end do
end do
! Node j=nlevsoi (bottom)
j = nlevsoi
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
if(j > jwt(c)) then !water table is in soil column
den = (zmm(c,j) - zmm(c,j-1))
dzq = (zq(c,j)-zq(c,j-1))
num = (smp(c,j)-smp(c,j-1)) - dzq
qin(c,j) = -hk(c,j-1)*num/den
dqidw0(c,j) = -(-hk(c,j-1)*dsmpdw(c,j-1) + num*dhkdw(c,j-1))/den
dqidw1(c,j) = -( hk(c,j-1)*dsmpdw(c,j) + num*dhkdw(c,j-1))/den
qout(c,j) = 0._r8
dqodw1(c,j) = 0._r8
rmx(c,j) = qin(c,j) - qout(c,j) - qflx_tran_veg_col(c)*rootr_col(c,j)
amx(c,j) = -dqidw0(c,j)
bmx(c,j) = dzmm(c,j)/dtime - dqidw1(c,j) + dqodw1(c,j)
cmx(c,j) = 0._r8
!scs: next set up aquifer layer; hydrologically inactive
rmx(c,j+1) = 0._r8
amx(c,j+1) = 0._r8
bmx(c,j+1) = dzmm(c,j+1)/dtime
cmx(c,j+1) = 0._r8
else ! water table is below soil column
!scs: compute aquifer soil moisture as average of layer 10 and saturation
s_node = max(0.5*(1.0_r8+h2osoi_vol(c,j)/watsat(c,j)), 0.01_r8)
s_node = min(1.0_r8, s_node)
!scs: compute smp for aquifer layer
smp1 = -sucsat(c,j)*s_node**(-bsw(c,j))
smp1 = max(smpmin(c), smp1)
!scs: compute dsmpdw for aquifer layer
dsmpdw1 = -bsw(c,j)*smp1/(s_node*watsat(c,j))
!scs: first set up bottom layer of soil column
den = (zmm(c,j) - zmm(c,j-1))
dzq = (zq(c,j)-zq(c,j-1))
num = (smp(c,j)-smp(c,j-1)) - dzq
qin(c,j) = -hk(c,j-1)*num/den
dqidw0(c,j) = -(-hk(c,j-1)*dsmpdw(c,j-1) + num*dhkdw(c,j-1))/den
dqidw1(c,j) = -( hk(c,j-1)*dsmpdw(c,j) + num*dhkdw(c,j-1))/den
den = (zmm(c,j+1)-zmm(c,j))
dzq = (zq(c,j+1)-zq(c,j))
num = (smp1-smp(c,j)) - dzq
qout(c,j) = -hk(c,j)*num/den
dqodw1(c,j) = -(-hk(c,j)*dsmpdw(c,j) + num*dhkdw(c,j))/den
dqodw2(c,j) = -( hk(c,j)*dsmpdw1 + num*dhkdw(c,j))/den
rmx(c,j) = qin(c,j) - qout(c,j) - qflx_tran_veg_col(c)*rootr_col(c,j)
amx(c,j) = -dqidw0(c,j)
bmx(c,j) = dzmm(c,j)/dtime - dqidw1(c,j) + dqodw1(c,j)
cmx(c,j) = dqodw2(c,j)
!scs: next set up aquifer layer; den/num unchanged, qin=qout
qin(c,j+1) = qout(c,j)
dqidw0(c,j+1) = -(-hk(c,j)*dsmpdw(c,j) + num*dhkdw(c,j))/den
dqidw1(c,j+1) = -( hk(c,j)*dsmpdw1 + num*dhkdw(c,j))/den
qout(c,j+1) = 0._r8 ! zero-flow bottom boundary condition
dqodw1(c,j+1) = 0._r8 ! zero-flow bottom boundary condition
rmx(c,j+1) = qin(c,j+1) - qout(c,j+1)
amx(c,j+1) = -dqidw0(c,j+1)
bmx(c,j+1) = dzmm(c,j+1)/dtime - dqidw1(c,j+1) + dqodw1(c,j+1)
cmx(c,j+1) = 0._r8
endif
end do
! Solve for dwat
jtop(:) = 1
call Tridiagonal(lbc, ubc, 1, nlevsoi+1, jtop, num_hydrologyc, filter_hydrologyc, &
amx, bmx, cmx, rmx, dwat2 )
!scs: set dwat
do fc = 1,num_hydrologyc
c = filter_hydrologyc(fc)
do j = 1, nlevsoi
dwat(c,j)=dwat2(c,j)
end do
end do
! Renew the mass of liquid water
!scs: also compute qcharge from dwat in aquifer layer
!scs: update in drainage for case jwt < nlevsoi
!dir$ concurrent
!cdir nodep
do fc = 1,num_hydrologyc
c = filter_hydrologyc(fc)
do j = 1, nlevsoi
h2osoi_liq(c,j) = h2osoi_liq(c,j) + dwat2(c,j)*dzmm(c,j)
end do
!scs: calculate qcharge for case jwt < nlevsoi
if(jwt(c) < nlevsoi) then
wh_zwt = 0._r8 !since wh_zwt = -sucsat - zq_zwt, where zq_zwt = -sucsat
s_node = max(h2osoi_vol(c,jwt(c))/watsat(c,jwt(c)), 0.01_r8)
s_node = min(1.0_r8, s_node)
!scs: use average moisture between water table and layer jwt
s1 = 0.5_r8*(1.0+s_node)
s1 = min(1._r8, s1)
!scs: this is the expression for unsaturated hk
ka = hksat(c,jwt(c))*s1**(2._r8*bsw(c,jwt(c))+3._r8)
! Recharge rate qcharge to groundwater (positive to aquifer)
smp1 = -sucsat(c,jwt(c))*s_node**(-bsw(c,jwt(c)))
smp1 = max(smpmin(c), smp(c,jwt(c)))
wh = smp1 - zq(c,jwt(c))
qcharge(c) = -ka * (wh_zwt-wh) /((zwt(c)-z(c,jwt(c)))*1000._r8)
! To limit qcharge (for the first several timesteps)
qcharge(c) = max(-10.0_r8/dtime,qcharge(c))
qcharge(c) = min( 10.0_r8/dtime,qcharge(c))
else
!scs: if water table is below soil column, compute qcharge from dwat2(11)
qcharge(c) = dwat2(c,nlevsoi+1)*dzmm(c,nlevsoi+1)/dtime
endif
end do
end subroutine SoilWater
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Drainage
!
! !INTERFACE:
subroutine Drainage(lbc, ubc, num_hydrologyc, filter_hydrologyc, & 1,6
num_urbanc, filter_urbanc, vol_liq, hk, &
icefrac)
!
! !DESCRIPTION:
! Calculate subsurface drainage
!
! !USES:
use shr_kind_mod, only : r8 => shr_kind_r8
use clmtype
use clm_time_manager, only : get_step_size
use clm_varcon , only : pondmx, tfrz, icol_roof, icol_road_imperv, icol_road_perv, watmin
use clm_varpar , only : nlevsoi
!
! !ARGUMENTS:
implicit none
integer , intent(in) :: lbc, ubc ! column bounds
integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter
integer , intent(in) :: num_urbanc ! number of column urban points in column filter
integer , intent(in) :: filter_urbanc(ubc-lbc+1) ! column filter for urban points
integer , intent(in) :: filter_hydrologyc(ubc-lbc+1) ! column filter for soil points
real(r8), intent(in) :: vol_liq(lbc:ubc,1:nlevsoi) ! partial volume of liquid water in layer
real(r8), intent(in) :: hk(lbc:ubc,1:nlevsoi) ! hydraulic conductivity (mm h2o/s)
real(r8), intent(in) :: icefrac(lbc:ubc,1:nlevsoi) ! fraction of ice in layer
!
! !CALLED FROM:
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 12 November 1999: Z.-L. Yang and G.-Y. Niu
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
! 4/26/05, Peter Thornton and David Lawrence: Turned off drainage from
! middle soil layers for both wet and dry fractions.
! 04/25/07 Keith Oleson: Completely new routine for CLM3.5 hydrology
! 27 February 2008: Keith Oleson; Saturation excess modification
!
! !LOCAL VARIABLES:
!
! local pointers to original implicit in arguments
!
integer , pointer :: ctype(:) !column type index
integer , pointer :: snl(:) !number of snow layers
real(r8), pointer :: qflx_snwcp_liq(:) !excess rainfall due to snow capping (mm H2O /s) [+]
real(r8), pointer :: qflx_dew_grnd(:) !ground surface dew formation (mm H2O /s) [+]
real(r8), pointer :: qflx_dew_snow(:) !surface dew added to snow pack (mm H2O /s) [+]
real(r8), pointer :: qflx_sub_snow(:) !sublimation rate from snow pack (mm H2O /s) [+]
real(r8), pointer :: dz(:,:) !layer depth (m)
real(r8), pointer :: bsw(:,:) !Clapp and Hornberger "b"
real(r8), pointer :: eff_porosity(:,:) !effective porosity = porosity - vol_ice
real(r8), pointer :: t_soisno(:,:) !soil temperature (Kelvin)
real(r8), pointer :: hksat(:,:) !hydraulic conductivity at saturation (mm H2O /s)
real(r8), pointer :: sucsat(:,:) !minimum soil suction (mm)
real(r8), pointer :: z(:,:) !layer depth (m)
real(r8), pointer :: zi(:,:) !interface level below a "z" level (m)
real(r8), pointer :: watsat(:,:) !volumetric soil water at saturation (porosity)
real(r8), pointer :: hkdepth(:) !decay factor (m)
real(r8), pointer :: zwt(:) !water table depth (m)
real(r8), pointer :: wa(:) !water in the unconfined aquifer (mm)
real(r8), pointer :: wt(:) !total water storage (unsaturated soil water + groundwater) (mm)
real(r8), pointer :: qcharge(:) !aquifer recharge rate (mm/s)
!
! local pointers to original implicit inout arguments
!
real(r8), pointer :: h2osoi_ice(:,:) !ice lens (kg/m2)
real(r8), pointer :: h2osoi_liq(:,:) !liquid water (kg/m2)
!
! local pointers to original implicit out arguments
!
real(r8), pointer :: qflx_drain(:) !sub-surface runoff (mm H2O /s)
real(r8), pointer :: qflx_qrgwl(:) !qflx_surf at glaciers, wetlands, lakes (mm H2O /s)
real(r8), pointer :: eflx_impsoil(:) !implicit evaporation for soil temperature equation
real(r8), pointer :: qflx_rsub_sat(:) !soil saturation excess [mm h2o/s]
!
!EOP
!
! !OTHER LOCAL VARIABLES:
!
!KO integer :: c,j,fc !indices
!KO
integer :: c,j,fc,i !indices
!KO
real(r8) :: dtime !land model time step (sec)
real(r8) :: xs(lbc:ubc) !water needed to bring soil moisture to watmin (mm)
real(r8) :: dzmm(lbc:ubc,1:nlevsoi) !layer thickness (mm)
integer :: jwt(lbc:ubc) !index of the soil layer right above the water table (-)
real(r8) :: rsub_bot(lbc:ubc) !subsurface runoff - bottom drainage (mm/s)
real(r8) :: rsub_top(lbc:ubc) !subsurface runoff - topographic control (mm/s)
real(r8) :: fff(lbc:ubc) !decay factor (m-1)
real(r8) :: xsi(lbc:ubc) !excess soil water above saturation at layer i (mm)
real(r8) :: xsia(lbc:ubc) !available pore space at layer i (mm)
real(r8) :: xs1(lbc:ubc) !excess soil water above saturation at layer 1 (mm)
real(r8) :: smpfz(1:nlevsoi) !matric potential of layer right above water table (mm)
real(r8) :: wtsub !summation of hk*dzmm for layers below water table (mm**2/s)
real(r8) :: rous !aquifer yield (-)
real(r8) :: wh !smpfz(jwt)-z(jwt) (mm)
real(r8) :: wh_zwt !water head at the water table depth (mm)
real(r8) :: ws !summation of pore space of layers below water table (mm)
real(r8) :: s_node !soil wetness (-)
real(r8) :: dzsum !summation of dzmm of layers below water table (mm)
real(r8) :: icefracsum !summation of icefrac*dzmm of layers below water table (-)
real(r8) :: fracice_rsub(lbc:ubc) !fractional impermeability of soil layers (-)
real(r8) :: ka !hydraulic conductivity of the aquifer (mm/s)
real(r8) :: dza !fff*(zwt-z(jwt)) (-)
!KO
real(r8) :: available_h2osoi_liq !available soil liquid water in a layer
!KO
!-----------------------------------------------------------------------
! Assign local pointers to derived subtypes components (column-level)
ctype => clm3%g%l%c%itype
! cgridcell => clm3%g%l%c%gridcell
snl => clm3%g%l%c%cps%snl
dz => clm3%g%l%c%cps%dz
bsw => clm3%g%l%c%cps%bsw
t_soisno => clm3%g%l%c%ces%t_soisno
hksat => clm3%g%l%c%cps%hksat
sucsat => clm3%g%l%c%cps%sucsat
z => clm3%g%l%c%cps%z
zi => clm3%g%l%c%cps%zi
watsat => clm3%g%l%c%cps%watsat
hkdepth => clm3%g%l%c%cps%hkdepth
zwt => clm3%g%l%c%cws%zwt
wa => clm3%g%l%c%cws%wa
wt => clm3%g%l%c%cws%wt
qcharge => clm3%g%l%c%cws%qcharge
eff_porosity => clm3%g%l%c%cps%eff_porosity
qflx_snwcp_liq => clm3%g%l%c%cwf%pwf_a%qflx_snwcp_liq
qflx_dew_grnd => clm3%g%l%c%cwf%pwf_a%qflx_dew_grnd
qflx_dew_snow => clm3%g%l%c%cwf%pwf_a%qflx_dew_snow
qflx_sub_snow => clm3%g%l%c%cwf%pwf_a%qflx_sub_snow
qflx_drain => clm3%g%l%c%cwf%qflx_drain
qflx_qrgwl => clm3%g%l%c%cwf%qflx_qrgwl
qflx_rsub_sat => clm3%g%l%c%cwf%qflx_rsub_sat
eflx_impsoil => clm3%g%l%c%cef%eflx_impsoil
h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq
h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice
! Get time step
dtime = get_step_size()
! Convert layer thicknesses from m to mm
do j = 1,nlevsoi
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
dzmm(c,j) = dz(c,j)*1.e3_r8
end do
end do
! Initial set
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
qflx_drain(c) = 0._r8
rsub_bot(c) = 0._r8
qflx_rsub_sat(c) = 0._r8
rsub_top(c) = 0._r8
fracice_rsub(c) = 0._r8
end do
! The layer index of the first unsaturated layer, i.e., the layer right above
! the water table
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
jwt(c) = nlevsoi
do j = 2,nlevsoi
if(zwt(c) <= zi(c,j)) then
jwt(c) = j-1
exit
end if
enddo
end do
! Topographic runoff
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
fff(c) = 1._r8/ hkdepth(c)
dzsum = 0._r8
icefracsum = 0._r8
do j = jwt(c), nlevsoi
dzsum = dzsum + dzmm(c,j)
icefracsum = icefracsum + icefrac(c,j) * dzmm(c,j)
end do
fracice_rsub(c) = max(0._r8,exp(-3._r8*(1._r8-(icefracsum/dzsum)))- exp(-3._r8))/(1.0_r8-exp(-3._r8))
rsub_top(c) = (1._r8 - fracice_rsub(c)) * 5.5e-3_r8 * exp(-fff(c)*zwt(c))
end do
rous = 0.2_r8
! Water table calculation
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
! Water storage in aquifer + soil
wt(c) = wt(c) + (qcharge(c) - rsub_top(c)) * dtime
if(jwt(c) == nlevsoi) then ! water table is below the soil column
wa(c) = wa(c) + (qcharge(c) -rsub_top(c)) * dtime
wt(c) = wa(c)
zwt(c) = (zi(c,nlevsoi) + 25._r8) - wa(c)/1000._r8/rous
h2osoi_liq(c,nlevsoi) = h2osoi_liq(c,nlevsoi) + max(0._r8,(wa(c)-5000._r8))
wa(c) = min(wa(c), 5000._r8)
else ! water table within soil layers
if (jwt(c) == nlevsoi-1) then ! water table within bottom soil layer
zwt(c) = zi(c,nlevsoi)- (wt(c)-rous*1000._r8*25._r8) /eff_porosity(c,nlevsoi)/1000._r8
else ! water table within soil layers 1-9
ws = 0._r8 ! water used to fill soil air pores regardless of water content
do j = jwt(c)+2,nlevsoi
ws = ws + eff_porosity(c,j) * 1000._r8 * dz(c,j)
enddo
zwt(c) = zi(c,jwt(c)+1)-(wt(c)-rous*1000_r8*25._r8-ws) /eff_porosity(c,jwt(c)+1)/1000._r8
endif
wtsub = 0._r8
do j = jwt(c)+1, nlevsoi
wtsub = wtsub + hk(c,j)*dzmm(c,j)
end do
! Remove subsurface runoff
do j = jwt(c)+1, nlevsoi
h2osoi_liq(c,j) = h2osoi_liq(c,j) - rsub_top(c)*dtime*hk(c,j)*dzmm(c,j)/wtsub
end do
end if
zwt(c) = max(0.05_r8,zwt(c))
zwt(c) = min(80._r8,zwt(c))
end do
! excessive water above saturation added to the above unsaturated layer like a bucket
! if column fully saturated, excess water goes to runoff
do j = nlevsoi,2,-1
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
xsi(c) = max(h2osoi_liq(c,j)-eff_porosity(c,j)*dzmm(c,j),0._r8)
h2osoi_liq(c,j) = min(eff_porosity(c,j)*dzmm(c,j), h2osoi_liq(c,j))
h2osoi_liq(c,j-1) = h2osoi_liq(c,j-1) + xsi(c)
end do
end do
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
xs1(c) = max(max(h2osoi_liq(c,1),0._r8)-max(0._r8,(pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_ice(c,1))),0._r8)
h2osoi_liq(c,1) = min(max(0._r8,pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_ice(c,1)), h2osoi_liq(c,1))
qflx_rsub_sat(c) = xs1(c) / dtime
end do
! Limit h2osoi_liq to be greater than or equal to watmin.
! Get water needed to bring h2osoi_liq equal watmin from lower layer.
! If insufficient water in soil layers, get from aquifer water
do j = 1, nlevsoi-1
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
!KO if (h2osoi_liq(c,j) < 0._r8) then
!KO
if (h2osoi_liq(c,j) < watmin) then
!KO
xs(c) = watmin - h2osoi_liq(c,j)
else
xs(c) = 0._r8
end if
h2osoi_liq(c,j ) = h2osoi_liq(c,j ) + xs(c)
h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) - xs(c)
end do
end do
!KO j = nlevsoi
!KO!dir$ concurrent
!KO!cdir nodep
!KO do fc = 1, num_hydrologyc
!KO c = filter_hydrologyc(fc)
!KO if (h2osoi_liq(c,j) < watmin) then
!KO xs(c) = watmin-h2osoi_liq(c,j)
!KO else
!KO xs(c) = 0._r8
!KO end if
!KO h2osoi_liq(c,j) = h2osoi_liq(c,j) + xs(c)
!KO wa(c) = wa(c) - xs(c)
!KO wt(c) = wt(c) - xs(c)
!KO end do
!KO
! Get water for bottom layer from layers above if possible
j = nlevsoi
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
if (h2osoi_liq(c,j) < watmin) then
xs(c) = watmin-h2osoi_liq(c,j)
searchforwater: do i = nlevsoi-1, 1, -1
available_h2osoi_liq = max(h2osoi_liq(c,i)-watmin-xs(c),0._r8)
if (available_h2osoi_liq .ge. xs(c)) then
h2osoi_liq(c,j) = h2osoi_liq(c,j) + xs(c)
h2osoi_liq(c,i) = h2osoi_liq(c,i) - xs(c)
xs(c) = 0._r8
exit searchforwater
else
h2osoi_liq(c,j) = h2osoi_liq(c,j) + available_h2osoi_liq
h2osoi_liq(c,i) = h2osoi_liq(c,i) - available_h2osoi_liq
xs(c) = xs(c) - available_h2osoi_liq
end if
end do searchforwater
else
xs(c) = 0._r8
end if
! Needed in case there is no water to be found
h2osoi_liq(c,j) = h2osoi_liq(c,j) + xs(c)
wt(c) = wt(c) - xs(c)
! Instead of removing water from aquifer where it eventually
! shows up as excess drainage to the ocean, take it back out of
! drainage
rsub_top(c) = rsub_top(c) - xs(c)/dtime
end do
!KO
!dir$ concurrent
!cdir nodep
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
! Sub-surface runoff and drainage
qflx_drain(c) = qflx_rsub_sat(c) + rsub_top(c)
! Set imbalance for snow capping
qflx_qrgwl(c) = qflx_snwcp_liq(c)
! Implicit evaporation term is now zero
eflx_impsoil(c) = 0._r8
! Renew the ice and liquid mass due to condensation
if (snl(c)+1 >= 1) then
h2osoi_liq(c,1) = h2osoi_liq(c,1) + qflx_dew_grnd(c) * dtime
h2osoi_ice(c,1) = h2osoi_ice(c,1) + (qflx_dew_snow(c) * dtime)
if (qflx_sub_snow(c)*dtime > h2osoi_ice(c,1)) then
qflx_sub_snow(c) = h2osoi_ice(c,1)/dtime
h2osoi_ice(c,1) = 0._r8
else
h2osoi_ice(c,1) = h2osoi_ice(c,1) - (qflx_sub_snow(c) * dtime)
end if
end if
end do
! No drainage for urban columns (except for pervious road as computed above)
!dir$ concurrent
!cdir nodep
do fc = 1, num_urbanc
c = filter_urbanc(fc)
if (ctype(c) /= icol_road_perv) then
qflx_drain(c) = 0._r8
! This must be done for roofs and impervious road (walls will be zero)
qflx_qrgwl(c) = qflx_snwcp_liq(c)
eflx_impsoil(c) = 0._r8
end if
! Renew the ice and liquid mass due to condensation for urban roof and impervious road
if (ctype(c) == icol_roof .or. ctype(c) == icol_road_imperv) then
if (snl(c)+1 >= 1) then
h2osoi_liq(c,1) = h2osoi_liq(c,1) + qflx_dew_grnd(c) * dtime
h2osoi_ice(c,1) = h2osoi_ice(c,1) + (qflx_dew_snow(c) * dtime)
if (qflx_sub_snow(c)*dtime > h2osoi_ice(c,1)) then
qflx_sub_snow(c) = h2osoi_ice(c,1)/dtime
h2osoi_ice(c,1) = 0._r8
else
h2osoi_ice(c,1) = h2osoi_ice(c,1) - (qflx_sub_snow(c) * dtime)
end if
end if
end if
end do
end subroutine Drainage
end module SoilHydrologyMod