INTERFACE:
subroutine SoilWater(lbc, ubc, num_hydrologyc, filter_hydrologyc, & 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_sizeARGUMENTS:
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 Hydrology2ModREVISION 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 hydrologyLOCAL 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)