module SoilHydrologyMod !----------------------------------------------------------------------- !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, & 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 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 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 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. 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, & 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) 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 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, & 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 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 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 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 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 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 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 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 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 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 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 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 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, & 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_irrig(:) !irrigation flux (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: ! integer :: c,j,fc,i !indices 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)) (-) real(r8) :: available_h2osoi_liq !available soil liquid water in a layer !----------------------------------------------------------------------- ! 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_irrig => clm3%g%l%c%cwf%qflx_irrig 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 do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) dzmm(c,j) = dz(c,j)*1.e3_r8 end do end do ! Initial set 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 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 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 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 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 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 do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) if (h2osoi_liq(c,j) < watmin) then 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 ! Get water for bottom layer from layers above if possible j = nlevsoi 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 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) do fc = 1, num_urbanc c = filter_urbanc(fc) if (ctype(c) /= icol_road_perv) then qflx_drain(c) = 0._r8 qflx_irrig(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