hydro_soildif.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #########
      SUBROUTINE hydro_soildif(IO, KK, PK, PEK, PTSTEP, PPG, PLETR, PLEG, PEVAPCOR, &
                               PF2WGHT, PPS, PQSAT, PQSATI, PDRAIN, PHORTON, KMAX_LAYER, PQSB)
!     ##########################################################################
!
!
!!****  *HYDRO_SOILDIF*  
!
!!    PURPOSE
!!    -------
!     This subroutine solves the 1-D (z) mass conservation equation
!     (mixed-form Richard's equation: tendency in terms of volumetric
!     water content, gradient in terms of matric potential)
!     for liquid water (using Darcy's Law for the vertical flux)
!     together with the Clapp and Hornberger (1978) simplification to
!     the Brooks and Corey (1966) empirical model for relating matric
!     potential and hydraulic conductivity to soil water content.
!     Any set of parameters can be used (eg. Clapp and Hornberger 1978;
!     Cosby et al. 1984; etc.) Modifications to the equations is also
!     made for the Van Genucten model/relationships. The equations
!     also incorporate vapor transfer for dry soils. The soil porosity
!     is modified in the presence of soil ice. Soil ice content
!     is also updated herein due to changes resulting from sublimation.
!     The layer averaged set of equations are time differenced
!     using an implicit time scheme. 
!     The equations/model *assume* a heterogeneous soil texture profile,
!     however, if the soil properties are homogeneous, the equations
!     collapse into the standard homogeneous approach (i.e. give the
!     same results as). The eqs are solved rapidly by taking advantage of the
!     fact that the matrix is tridiagonal. 
!     Note that the exponential profile of hydraulic conductivity with soil
!     depth is applied to interfacial conductivity (or interblock)
!
!     
!!**  METHOD
!!    ------
!
!     Direct calculation
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    USE MODD_CST
!!    USE MODI_TRIDIAG_GROUND
!!      
!!    REFERENCE
!!    ---------
!!
!!    Boone (2000)
!!    Boone et al. (2000)
!!      
!!    AUTHOR
!!    ------
!!      A. Boone          * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    16/02/00 Boone
!!      Modif       04/2010  B.Decharme: geometric mean for interfacial conductivity
!!                                       Brook and Corey 
!!      Modif       08/2011  B.Decharme: Optimization using global loops
!!                  10/12    B.Decharme: EVAPCOR snow correction in DIF
!!                  04/13    B.Decharme: Subsurface runoff if SGH (DIF option only)
!!                  07/2013  B.Decharme: Surface / Water table depth coupling
!!                  01/2016  B.Decharme: Bug : if no surface runoff (HRUNOFF=WSAT) then no Horton
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_k_t, isba_p_t, isba_pe_t
!
USE modd_surf_par, ONLY : xundef, nundef
USE modd_csts,     ONLY : xrholw
USE modd_isba_par, ONLY : xwgmin
!
USE mode_hydro_dif
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_k_t), INTENT(INOUT) :: KK
TYPE(isba_p_t), INTENT(INOUT) :: PK
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
!
REAL, INTENT(IN)                    :: PTSTEP ! Model time step (s)
!
REAL, DIMENSION(:), INTENT(IN)      :: PPS, PPG, PLETR, PLEG, PEVAPCOR
!                                      PPS    = surface pressure (Pa)
!                                      PPG    = throughfall rate: 
!                                               rate at which water reaches the surface
!                                               of the soil (from snowmelt, rain, canopy
!                                               drip, etc...) (m/s)
!                                      PLETR  = transpiration rate (m/s)
!                                      PLEG   = bare-soil evaporation rate (m/s)
!                                      PEVAPCOR = correction for any excess evaporation 
!                                                from snow as it completely ablates (m/s)
!
REAL, DIMENSION(:,:), INTENT(IN)    :: PQSAT,PQSATI
!                                      specific humidity at saturation
!
REAL, DIMENSION(:,:), INTENT(IN)    :: PF2WGHT
!                                      PF2WGHT      = root-zone transpiration weights (-)
!
INTEGER,               INTENT(IN)   :: KMAX_LAYER  
!                                      KMAX_LAYER = Max number of soil moisture layers (DIF option)
!
REAL, DIMENSION(:), INTENT(OUT)     :: PDRAIN, PHORTON
!                                      PDRAIN   = drainage (flux out of model base) (kg m-2 s-1)
!
REAL, DIMENSION(:),   INTENT(OUT)   :: PQSB     !Lateral subsurface flow [kg/m²/s]
!
!
!*      0.2    declarations of local variables
!
INTEGER                             :: JJ, JL    ! loop control
!
INTEGER                             :: INJ, INL, IDEPTH ! Number of point and grid layers
!
REAL, DIMENSION(SIZE(PK%XDZG,1))       :: ZINFILTC, ZEXCESS, ZDGN, ZWGTOT, ZPSIWTD, ZWTD
!                                      ZEXCESS    = working variable: excess soil water
!                                                   which is used as a constraint 
!                                      ZDGN       = Depth of the last node (m)
!                                                   and the water table (m s-1)
!                                      ZWGTOT    = total soil moisture for ZINFNEG computation
!                                      ZPSIWTD   = matric potential at saturation for water table depth coupling
!                                      ZWTD      = water table depth positive below soil surface (m)
!
REAL, DIMENSION(SIZE(PK%XDZG,1),SIZE(PK%XDZG,2)) :: ZWFLUX, ZDFLUXDT1, ZDFLUXDT2, ZWFLUXN
!                                      ZWFLUX    = vertical soil water flux (+ up) (m s-1)
!                                      ZDFLUXDT  = total vertical flux derrivative
!                                      ZDFLUXDT1 = vertical flux derrivative: dF_j/dw_j
!                                      ZDFLUXDT2 = vertical flux derrivative: dF_j/dw_j+1
!                                      ZWFLUXN   = vertical soil water flux at end of time 

REAL, DIMENSION(SIZE(PK%XDZG,1),SIZE(PK%XDZG,2)) :: ZPSI, ZK, ZNU, ZWSAT, ZHEAD, &
                                              ZVAPCOND, ZFRZ, ZKI, ZCAPACITY, ZINFNEG
!                                      ZNU       = interfacial total conductivity (m s-1)
!                                                  at level z_j
!                                      ZK        = hydraulic conductivity (m s-1)
!                                      ZHEAD     = matric potential gradient (-)
!                                      ZWSAT     = ice modified soil porosity (m3 m-3)
!                                      ZFRZ      = diffusion coefficient for freezing (-)
!                                      ZVAPCOND  = vapor conductivity (m s-1) 
!                                      ZKI       = interfacial hydraulic conductivity (m s-1) at level z_j
!                                      ZCAPACITY = simple volumetric water holding capacity estimate for
!                                                  wetting front penetration (-) 
!                                      ZINFNEG   = Negative infiltration (m s-1)
!
REAL, DIMENSION(SIZE(PK%XDZG,1),SIZE(PK%XDZG,2)) :: ZAMTRX, ZBMTRX, ZCMTRX, ZFRC, ZSOL, &
                                              ZTOPQS
!                                      ZAMTRX    = leftmost diagonal element of tri-diagonal
!                                                  coefficient matrix 
!                                      ZBMTRX    = center diagonal element (vector)
!                                      ZCMTRX    = rightmost diagonal element (vector)
!                                      ZFRC      = forcing function (vector)
!                                      ZSOL      = solution vector
!
REAL, DIMENSION(SIZE(PK%XDZG,1),SIZE(PK%XDZG,2))  :: ZINFLAYER
!
REAL, PARAMETER                     :: ZWGHT = 0.5  ! time scheme weight for calculating flux.
!                                                     varies from 0 (explicit time scheme)
!                                                     to 1 (backward difference implicit)
!                                                     Default is 1/2 (Crank-Nicholson)
!
REAL, PARAMETER                     :: ZEICE = 6.0  ! Ice vertical diffusion impedence factor 
!
REAL                                :: ZLOG10, ZS, ZLOG, ZWDRAIN
!
REAL                                :: ZDKDT1, ZDKDT2, ZDHEADDT1, ZDHEADDT2
!                                      ZDKDT1    = hydraulic conductivity derrivative w/r/t upper layer water content
!                                      ZDKDT2    = "" lower layer water content
!                                      ZDHEADDT1 = matric potential gradient derrivative w/r/t upper layer water content
!                                      ZDHEADDT2 = "" lower layer water content
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
! 0. Initialization:
!    ---------------
!
IF (lhook) CALL dr_hook('HYDRO_SOILDIF',0,zhook_handle)
!
inj = SIZE(pk%XDZG,1)
inl = kmax_layer
!
zlog10 = log(10.0)
!
pdrain(:) = 0.0
pqsb(:) = 0.0
phorton(:) = 0.0
zinfiltc(:) = 0.0
zexcess(:) = 0.0
!
zdgn(:) = xundef
zpsiwtd(:) = xundef
zwtd(:) = xundef
!
zinfneg(:,:) = 0.0
zinflayer(:,:) = 0.0
zfrc(:,:) = 0.0
zfrz(:,:) = 0.0
!
zwsat(:,:) = xundef
zcapacity(:,:) = xundef
zpsi(:,:) = xundef
zk(:,:) = xundef
zvapcond(:,:) = xundef
zki(:,:) = xundef
zsol(:,:) = xundef
znu(:,:) = xundef
zhead(:,:) = xundef
zwflux(:,:) = xundef
zwfluxn(:,:) = xundef
zdfluxdt1(:,:) = xundef
zdfluxdt2(:,:) = xundef
zamtrx(:,:) = xundef
zbmtrx(:,:) = xundef
zcmtrx(:,:) = xundef
!
! Modification/addition of frozen soil parameters
! -----------------------------------------------
!
DO jl=1,inl
  DO jj=1,inj
!
    idepth=pk%NWG_LAYER(jj)
    IF(jl<=idepth)THEN
!
!     Modify soil porosity as ice assumed to become part
!     of solid soil matrix (with respect to liquid flow):
      zwsat(jj,jl) = max(xwgmin, kk%XWSAT(jj,jl)-pek%XWGI(jj,jl))   
!
!     Factor from (Johnsson and Lundin 1991), except here it is normalized so that it
!     goes to zero in the limit as all available pore space is filled up with ice.
!     For now, a simple constant is used for all soils. Further modifications
!     will be made as research warrents.
!     Old : 10.**(-ZEICE*(PEK%XWGI(JJ,JL)/(PEK%XWGI(JJ,JL)+PEK%XWG(JJ,JL))))
      zfrz(jj,jl) = exp(zlog10*(-zeice*(pek%XWGI(jj,jl)/(pek%XWGI(jj,jl)+pek%XWG(jj,jl)))))
!
    ENDIF
!
  ENDDO
ENDDO
!
! Lateral sub-surface flow (m s-1) if Topmodel
! --------------------------------------------
!
IF(io%CRUNOFF=='SGH')THEN
  ztopqs(:,:)=zfrz(:,:)*pk%XTOPQS(:,:)
ELSE
  ztopqs(:,:)=0.0
ENDIF
!
! 1. Infiltration at "t"
!    -------------------
!
! Surface flux term (m s-1): Infiltration (and surface runoff)
! Surface fluxes are limited a Green-Ampt approximation from Abramopoulos et al
! (1988) and Entekhabi and Eagleson (1989).
! Note : when Horton option is used, infiltration already calculated in hydro_sgh
!
!Surface cumulative infiltration  (m)
zinfiltc(:) = max(0.0,ppg(:))*ptstep

!
! 2. Initialise soil moisture profile according to infiltration terms at "t"
!    ----------------------------------------------------------------------
!
DO jl=1,inl
  DO jj=1,inj
    idepth=pk%NWG_LAYER(jj)
    IF(jl<=idepth)THEN
!     Simple volumetric water holding capacity estimate for wetting front penetration
      zcapacity(jj,jl) = max(0.0,zwsat(jj,jl)-pek%XWG(jj,jl))*pk%XDZG(jj,jl)
!     Infiltration terms (m) :
      zinflayer(jj,jl) = min(zinfiltc(jj),zcapacity(jj,jl))
!     Soil moisture (m3/m3) :
      pek%XWG(jj,jl) = pek%XWG(jj,jl)+zinflayer(jj,jl)/pk%XDZG(jj,jl)
!     Put remainding infiltration into the next layer (m)
      zinfiltc(jj) = zinfiltc(jj) - zinflayer(jj,jl)
    ENDIF
  ENDDO
ENDDO
!
! 3. Compute Fast(temporal)-response runoff and Possible negative infiltration
!    -------------------------------------------------------------------------
!
!Possible negative infiltration  (m s-1)
zwgtot(:)=0.0 
DO jl=1,inl
  DO jj=1,inj
    idepth=pk%NWG_LAYER(jj)
    IF(jl<idepth)THEN
      zinfneg(jj,jl) = (min(0.0,ppg(jj))-pevapcor(jj))*pk%XDZG(jj,jl)*pek%XWG(jj,jl)
      zwgtot(jj   ) = zwgtot(jj)+pk%XDZG(jj,jl)*pek%XWG(jj,jl)
    ENDIF
  ENDDO
ENDDO
DO jl=1,inl
  DO jj=1,inj
    zinfneg(jj,jl) = zinfneg(jj,jl)/zwgtot(jj)
  ENDDO
ENDDO 
!
!Fast(temporal)-response runoff (surface excess) (kg m2 s-1):
!special case : if infiltration > total soil capacity
phorton(:)=(phorton(:)+zinfiltc(:)/ptstep)*xrholw
!
!
! 4. Initialise matric potential and hydraulic conductivity at "t"
!    -------------------------------------------------------------
!
DO jl=1,inl
  DO jj=1,inj    
    idepth=pk%NWG_LAYER(jj)
    IF(jl<=idepth)THEN
!     Matric potential (m) :
!     psi=mpotsat*(w/wsat)**(-bcoef)
      zs          = min(1.0,pek%XWG(jj,jl)/zwsat(jj,jl))
      zlog        = kk%XBCOEF(jj,jl)*log(zs)
      zpsi(jj,jl) = kk%XMPOTSAT(jj,jl)*exp(-zlog)
!     Hydraulic conductivity from matric potential (m s-1):
!     k=frz*condsat*(psi/mpotsat)**(-2-3/bcoef)
      zlog      = -zlog*(2.0+3.0/kk%XBCOEF(jj,jl))
      zk(jj,jl) = zfrz(jj,jl)*pk%XCONDSAT(jj,jl)*exp(-zlog)
    ENDIF
  ENDDO
ENDDO    
!
! Prepare water table depth coupling
! ----------------------------------
!
DO jj=1,inj
  idepth=pk%NWG_LAYER(jj)   
! Depth of the last node
  zdgn(jj) = 0.5*(pk%XDG(jj,idepth)+pk%XDG(jj,idepth-1))
  zpsiwtd(jj) = zpsi(jj,idepth)
  IF(kk%XWTD(jj)/=xundef)THEN  
!   Water table depth
    zwtd(jj)    = max(pk%XDG(jj,idepth),-kk%XWTD(jj))
!   Modify matric potential at saturation for water table coupling
    zs          = min(1.0,zwsat(jj,idepth)/kk%XWSAT(jj,idepth))
    zlog        = kk%XBCOEF(jj,idepth)*log(zs)
    zpsiwtd(jj) = kk%XMPOTSAT(jj,idepth)*exp(-zlog)
  ENDIF
ENDDO
!
! 5. Vapor diffusion conductivity (m s-1)
!    ------------------------------------
!
zvapcond(:,:) = vapcondcf(pek%XTG, pps, pek%XWG, pek%XWGI, zpsi,&
                          kk%XWSAT, kk%XWFC, pqsat, pqsati, pk%NWG_LAYER, inl)
zvapcond(:,:) = zfrz(:,:)*zvapcond(:,:)
!
! 6. Linearized water flux: values at "t"
!    ------------------------------------
!    calculate flux at the beginning of the time step:
!
DO jl=1,inl
  DO jj=1,inj
!
    idepth=pk%NWG_LAYER(jj)
    IF(jl<idepth)THEN
!
!     Total interfacial conductivity (m s-1) And Potential gradient (dimensionless):
      zki(jj,jl) = sqrt(zk(jj,jl)*zk(jj,jl+1))
      znu(jj,jl) = zki(jj,jl) + sqrt(zvapcond(jj,jl)*zvapcond(jj,jl+1))
      zhead(jj,jl) = (zpsi(jj,jl)-zpsi(jj,jl+1))/pk%XDZDIF(jj,jl)
!
!     Total Sub-surface soil water fluxes (m s-1): (+ up, - down) using Darcy's
!     Law with an added linear drainage term:
      zwflux(jj,jl) = -znu(jj,jl) * zhead(jj,jl) - zki(jj,jl)
!
    ELSEIF(jl==idepth)THEN !Last layers   
!        
!     Total interfacial conductivity (m s-1) And Potential gradient (dimensionless):
      zki(jj,idepth) = zk(jj,idepth)
      znu(jj,idepth) = zk(jj,idepth) * kk%XFWTD(jj)
      zhead(jj,idepth) = (zpsi(jj,idepth)-zpsiwtd(jj))/(zwtd(jj)-zdgn(jj))
!
!     Total Sub-surface soil water fluxes (m s-1): (+ up, - down) using Darcy's
!     Law with an added linear drainage term:
      zwflux(jj,idepth) = -znu(jj,idepth) * zhead(jj,idepth) - zki(jj,idepth)
!
    ENDIF
!
  ENDDO
ENDDO
!
!
! 7. Linearized water flux: values at "t+dt"
!    ---------------------------------------
! Flux Derrivative terms, see A. Boone thesis (Annexe E).
!
DO jl=1,inl
  DO jj=1,inj
    idepth=pk%NWG_LAYER(jj)        
    IF(jl<idepth)THEN                
      zdheaddt1 = -kk%XBCOEF(jj,jl  )*zpsi(jj,jl  )/(pek%XWG(jj,jl  )*pk%XDZDIF(jj,jl))
      zdheaddt2 = -kk%XBCOEF(jj,jl+1)*zpsi(jj,jl+1)/(pek%XWG(jj,jl+1)*pk%XDZDIF(jj,jl))
      zdkdt1    = (2.*kk%XBCOEF(jj,jl  )+3.)*zki(jj,jl)/(2.0*pek%XWG(jj,jl  ))
      zdkdt2    = (2.*kk%XBCOEF(jj,jl+1)+3.)*zki(jj,jl)/(2.0*pek%XWG(jj,jl+1))
!     Total Flux derrivative terms:
      zdfluxdt1(jj,jl) = -zdkdt1*zhead(jj,jl) - znu(jj,jl)*zdheaddt1 - zdkdt1
      zdfluxdt2(jj,jl) = -zdkdt2*zhead(jj,jl) + znu(jj,jl)*zdheaddt2 - zdkdt2  
    ELSEIF(jl==idepth)THEN !Last layers
      zdheaddt1 = -kk%XBCOEF(jj,idepth)*zpsi(jj,idepth)/(pek%XWG  (jj,idepth)*(zwtd(jj)-zdgn(jj))) &
                  +kk%XBCOEF(jj,idepth)*zpsiwtd(jj       )/(zwsat(jj,idepth)    *(zwtd(jj)-zdgn(jj)))
      zdheaddt2 = 0.0
      zdkdt1    = (2.*kk%XBCOEF(jj,idepth)+3.)*zk(jj,idepth)/pek%XWG(jj,idepth)
      zdkdt2    = 0.0                
!     Total Flux derrivative terms:
      zdfluxdt1(jj,idepth) = -zdkdt1*zhead(jj,idepth)*kk%XFWTD(jj) - znu(jj,idepth)*zdheaddt1 - zdkdt1
      zdfluxdt2(jj,idepth) = -zdkdt2*zhead(jj,idepth)*kk%XFWTD(jj) + znu(jj,idepth)*zdheaddt2 - zdkdt2  
    ENDIF
  ENDDO
ENDDO
!
! 8. Jacobian Matrix coefficients and Forcing function
!    -------------------------------------------------
!     
!surface layer:
zfrc(:,1) = zwflux(:,1) - pleg(:) - pf2wght(:,1)*pletr(:) + zinfneg(:,1) - ztopqs(:,1)
zamtrx(:,1) = 0.0
zbmtrx(:,1) = (pk%XDZG(:,1)/ptstep) - zwght*zdfluxdt1(:,1)
zcmtrx(:,1) = -zwght*zdfluxdt2(:,1)
!
!Other sub-surface layers:       
DO jl=2,inl
  DO jj=1,inj   
    idepth=pk%NWG_LAYER(jj)
    IF(jl<=idepth)THEN
      zfrc(jj,jl) = zwflux(jj,jl) - zwflux(jj,jl-1) - pf2wght(jj,jl)*pletr(jj) + zinfneg(jj,jl) - ztopqs(jj,jl)
      zamtrx(jj,jl) = zwght*zdfluxdt1(jj,jl-1)
      zbmtrx(jj,jl) = (pk%XDZG(jj,jl)/ptstep) - zwght*(zdfluxdt1(jj,jl)-zdfluxdt2(jj,jl-1))       
      zcmtrx(jj,jl) = -zwght*zdfluxdt2(jj,jl)
    ENDIF
  ENDDO
ENDDO
!
!Solve Matrix Equation: tridiagonal system: solve for soil
!water (volumetric water content) tendencies:
!
 CALL tridiag_dif(zamtrx,zbmtrx,zcmtrx,zfrc,pk%NWG_LAYER(:),inl,zsol)
!
! 9. Final calculations and diagnostics:
!    -----------------------------------
!
!
DO jl=1,inl
  DO jj=1,inj
!   
    idepth=pk%NWG_LAYER(jj)
    IF(jl<idepth)THEN
! 
!     Update liquid water content (m3 m-3):
      pek%XWG(jj,jl)   = pek%XWG(jj,jl) + zsol(jj,jl)    
!
!     Supersaturated drainage (kg m-2 s-1):
      zexcess(jj)  = max(0.0, pek%XWG(jj,jl) - zwsat(jj,jl))
      pek%XWG(jj,jl  ) = min(pek%XWG(jj,jl),zwsat(jj,jl))
      pek%XWG(jj,jl+1) = pek%XWG(jj,jl+1) + zexcess(jj)*(pk%XDZG(jj,jl)/pk%XDZG(jj,jl+1))
!
!     final fluxes (at end of time step) (m s-1):
      zwfluxn(jj,jl) = zwflux(jj,jl) + zdfluxdt1(jj,jl)*zsol(jj,jl) + zdfluxdt2(jj,jl)*zsol(jj,jl+1)
!
!     Total topmodel subsurface flow
      pqsb(jj) = pqsb(jj) + ztopqs(jj,jl)*xrholw
!
    ELSEIF(jl==idepth)THEN
! 
!     Update liquid water content (m3 m-3):
      pek%XWG(jj,idepth)   = pek%XWG(jj,idepth) + zsol(jj,idepth)    
!        
!     Supersaturated drainage (kg m-2 s-1):
      zexcess(jj)    = max(0.0, pek%XWG(jj,idepth) - zwsat(jj,idepth))
      pek%XWG(jj,idepth) = min(pek%XWG(jj,idepth),zwsat(jj,idepth))   
      pdrain(jj)    = pdrain(jj) + zexcess(jj)*pk%XDZG(jj,idepth)*xrholw/ptstep
!   
!     final fluxes (at end of time step) (m s-1):
      zwfluxn(jj,idepth) = zwflux(jj,idepth) + zdfluxdt1(jj,idepth)*zsol(jj,idepth)
!   
!     Drainage or baseflow out of bottom of model (slow time response) (kg m-2 s-1):
!     Final fluxes (if needed) over the time step (kg m-2 s-1)
!     would be calculated as (for water budget checks) as F = [ wgt*F(t+dt) + (1.-wgt)*F(t) ]*XRHOLW
      pdrain(jj) = pdrain(jj)-(zwght*zwfluxn(jj,idepth)+(1.-zwght)*zwflux(jj,idepth))*xrholw
!
!     Total topmodel subsurface flow
      pqsb(jj) = pqsb(jj) + ztopqs(jj,idepth)*xrholw
!
    ENDIF
!
  ENDDO
ENDDO
!
! Possible correction/Constraint application: 
!
DO jl=1,inl
  DO jj=1,inj   
    idepth=pk%NWG_LAYER(jj)
    IF(jl<idepth)THEN
!     if the soil water happens to fall below the minimum, then
!     extract needed water from the layer below: this should
!     generally be non-existant: but added to ensure conservation
!     even for the most extreme events.              
      zexcess(jj)  = max(0., xwgmin  - pek%XWG(jj,jl))
      pek%XWG(jj,jl)   = pek%XWG(jj,jl)   + zexcess(jj) 
      pek%XWG(jj,jl+1) = pek%XWG(jj,jl+1) - zexcess(jj)*pk%XDZG(jj,jl)/pk%XDZG(jj,jl+1)
    ELSEIF(jl==idepth.AND.pek%XWG(jj,idepth)<xwgmin)THEN
!     NOTE, negative moisture can arise for *completely* dry/dessicated soils 
!     owing to the above check because vertical fluxes
!     can be *very* small but nonzero. Here correct owing to
!     small numerical drainage.
      pdrain(jj)     = pdrain(jj) + min(0.0,pek%XWG(jj,idepth)-xwgmin)*pk%XDZG(jj,idepth)*xrholw/ptstep
      pek%XWG(jj,idepth) = xwgmin
    ENDIF
  ENDDO
ENDDO
!
IF (lhook) CALL dr_hook('HYDRO_SOILDIF',1,zhook_handle)
!
END SUBROUTINE hydro_soildif