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(HRUNOFF,HHORT,PTSTEP,                         &
                               pbcoef,pwsat,pcondsat,pmpotsat,pwfc,pdg,pdzg, &
                               pdzdif,ppg,pletr,pleg,pevapcor,pf2wght,       &
                               pwg,pwgi,ptg,pps,pqsat,pqsati,                &
                               pdrain,phorton,pfsat,kwg_layer,               &
                               kmax_layer,klayer_hort,ptopqs,pqsb,pfwtd,pwtd )
!     ##########################################################################
!
!
!!****  *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_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
!
 CHARACTER(LEN=*),     INTENT(IN)    :: hrunoff ! surface runoff formulation
!                                              ! 'WSAT'
!                                              ! 'DT92'
!                                              ! 'SGH ' Topmodel
!
 CHARACTER(LEN=*),     INTENT(IN)    :: hhort    ! Hortonian runoff
!
REAL, INTENT(IN)                    :: ptstep ! Model time step (s)
!
REAL, DIMENSION(:), INTENT(IN)      :: pps, ppg, pletr, pleg, pevapcor, pfsat, &
                                       pfwtd, pwtd
!                                      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)
!                                      PFSAT  = Saturated fraction
!                                      PFWTD  = grid-cell fraction of water table to rise
!                                      PWTD   = water table depth negative below soil surface (m)
!
REAL, DIMENSION(:,:), INTENT(IN)    :: pqsat,pqsati
!                                      specific humidity at saturation
!
REAL, DIMENSION(:,:), INTENT(IN)    :: ptg, pdg, pdzg, pdzdif
!                                      PTG = layer-average soil temperatures (K)
!                                      PDG = soil layer depth       (m)
!                                      PDZG= soil layer thicknesses (m)
!                                      PDZDIF = distance between the midpoints of consecutive layers (m)
!
REAL, DIMENSION(:,:), INTENT(IN)    :: pwsat, pcondsat, pf2wght, pwfc
!                                      PWSAT        = porosity profile (m3 m-3)
!                                      PCONDSAT     = hydraulic conductivity at saturation (m s-1)
!                                      PF2WGHT      = root-zone transpiration weights (-)
!                                      PWFC         = field capacity water content (m3 m-3)
!
REAL, DIMENSION(:,:), INTENT(IN)    :: pbcoef,pmpotsat
!                                      PMPOTSAT = matric potential at saturation (m) (BC parameters)
!                                      PBCOEF   = slope of the retention curve (-) (BC parameters)
!
INTEGER, DIMENSION(:), INTENT(IN)   :: kwg_layer  
!                                      KWG_LAYER = Number of soil moisture layers
!
INTEGER,               INTENT(IN)   :: kmax_layer  
!                                      KMAX_LAYER = Max number of soil moisture layers (DIF option)
!
INTEGER,               INTENT(IN)   :: klayer_hort! DIF optimization
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: pwg, pwgi
!                                      PWG  = volumetric liquid water content (m3 m-3) 
!                                      PWGI = volumetric ice content (m3 m-3)
!
REAL, DIMENSION(:), INTENT(OUT)     :: pdrain, phorton
!                                      PDRAIN   = drainage (flux out of model base) (kg m-2 s-1)
!                                      PHORTON  = runoff (due to saturation (lateral) (kg m-2 s-1)

REAL, DIMENSION(:,:), INTENT(IN)    :: ptopqs
!                                      PTOPQS = Topmodel subsurface flow by layer (m/s)
!
REAL, DIMENSION(:),   INTENT(OUT)   :: pqsb     !Lateral subsurface flow [kg/m²/s]
!
!
!*      0.2    declarations of local variables
!
INTEGER                             :: jj, jl    ! loop control
!
INTEGER                             :: ini, inl, idepth ! Number of point and grid layers
!
REAL, DIMENSION(SIZE(PDZG,1))       :: zinfiltmax, zinfiltc, zexcess, zdgn, zwgtot, zpsiwtd, zwtd
!                                      ZINFILTMAX = maximum allowable infiltration rate
!                                                   (from Darcy's Law) (m s-1)
!                                      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(PDZG,1),SIZE(PDZG,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(PDZG,1),SIZE(PDZG,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(PDZG,1),SIZE(PDZG,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(PDZG,1),SIZE(PDZG,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)
!
ini = SIZE(pdzg(:,:),1)
inl = kmax_layer
!
zlog10 = log(10.0)
!
pdrain(:) = 0.0
pqsb(:) = 0.0
phorton(:) = 0.0
zinfiltc(:) = 0.0
zexcess(:) = 0.0
zinfiltmax(:) = 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,ini
!
      idepth=kwg_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, pwsat(jj,jl)-pwgi(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*(PWGI(JJ,JL)/(PWGI(JJ,JL)+PWG(JJ,JL))))
        zfrz(jj,jl) = exp(zlog10*(-zeice*(pwgi(jj,jl)/(pwgi(jj,jl)+pwg(jj,jl)))))
!
     ENDIF
!
   ENDDO
ENDDO
!
! Lateral sub-surface flow (m s-1) if Topmodel
! --------------------------------------------
!
  IF(hrunoff=='SGH')THEN
    ztopqs(:,:)=zfrz(:,:)*ptopqs(:,:)
  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
IF(hrunoff/='WSAT'.AND.hhort/='SGH')THEN
!  Green-Ampt approximation for maximum infiltration (derived form)
   zinfiltmax(:) = infmax_func(pwg,zwsat,zfrz,pcondsat,pmpotsat,pbcoef,pdzg,pdg,klayer_hort)
!  Fast(temporal)-response runoff (surface excess) (m s-1):
   phorton(:) = (1.-pfsat(:)) * max(0.0,ppg(:)-zinfiltmax(:))
ENDIF
!
!
! 2. Initialise soil moisture profile according to infiltration terms at "t"
!    ----------------------------------------------------------------------
!
!Surface cumulative infiltration  (m)
zinfiltc(:) = max(0.0,ppg(:)-phorton(:))*ptstep
!
DO jl=1,inl
   DO jj=1,ini
      idepth=kwg_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)-pwg(jj,jl))*pdzg(jj,jl)
!       Infiltration terms (m) :
        zinflayer(jj,jl) = min(zinfiltc(jj),zcapacity(jj,jl))
!       Soil moisture (m3/m3) :
        pwg(jj,jl) = pwg(jj,jl)+zinflayer(jj,jl)/pdzg(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,ini
      idepth=kwg_layer(jj)
      IF(jl<idepth)THEN
        zinfneg(jj,jl) = (min(0.0,ppg(jj))-pevapcor(jj))*pdzg(jj,jl)*pwg(jj,jl)
        zwgtot(jj   ) = zwgtot(jj)+pdzg(jj,jl)*pwg(jj,jl)
      ENDIF
   ENDDO
ENDDO
DO jl=1,inl
   DO jj=1,ini
      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,ini    
     idepth=kwg_layer(jj)
     IF(jl<=idepth)THEN
!       Matric potential (m) :
!       psi=mpotsat*(w/wsat)**(-bcoef)
        zs          = min(1.0,pwg(jj,jl)/zwsat(jj,jl))
        zlog        = pbcoef(jj,jl)*log(zs)
        zpsi(jj,jl) = pmpotsat(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/pbcoef(jj,jl))
        zk(jj,jl) = zfrz(jj,jl)*pcondsat(jj,jl)*exp(-zlog)
     ENDIF
   ENDDO
ENDDO    
!
! Prepare water table depth coupling
! ----------------------------------
!
DO jj=1,ini
   idepth=kwg_layer(jj)   
!  Depth of the last node
   zdgn(jj) = 0.5*(pdg(jj,idepth)+pdg(jj,idepth-1))
   zpsiwtd(jj) = zpsi(jj,idepth)
   IF(pwtd(jj)/=xundef)THEN  
!    Water table depth
     zwtd(jj)    = max(pdg(jj,idepth),-pwtd(jj))
!    Modify matric potential at saturation for water table coupling
     zs          = min(1.0,zwsat(jj,idepth)/pwsat(jj,idepth))
     zlog        = pbcoef(jj,idepth)*log(zs)
     zpsiwtd(jj) = pmpotsat(jj,idepth)*exp(-zlog)
   ENDIF
ENDDO
!
! 5. Vapor diffusion conductivity (m s-1)
!    ------------------------------------
!
zvapcond(:,:) = vapcondcf(ptg,pps,pwg,pwgi,zpsi,pwsat,pwfc,pqsat,pqsati,kwg_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,ini
!
      idepth=kwg_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))/pdzdif(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) * pfwtd(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,ini
      idepth=kwg_layer(jj)        
      IF(jl<idepth)THEN                
         zdheaddt1 = -pbcoef(jj,jl  )*zpsi(jj,jl  )/(pwg(jj,jl  )*pdzdif(jj,jl))
         zdheaddt2 = -pbcoef(jj,jl+1)*zpsi(jj,jl+1)/(pwg(jj,jl+1)*pdzdif(jj,jl))
         zdkdt1    = (2.*pbcoef(jj,jl  )+3.)*zki(jj,jl)/(2.0*pwg(jj,jl  ))
         zdkdt2    = (2.*pbcoef(jj,jl+1)+3.)*zki(jj,jl)/(2.0*pwg(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 = -pbcoef(jj,idepth)*zpsi(jj,idepth)/(pwg(jj,idepth)*(zwtd(jj)-zdgn(jj))) &
                     +pbcoef(jj,idepth)*zpsiwtd(jj       )/(zwsat(jj,idepth)*(zwtd(jj)-zdgn(jj)))
         zdheaddt2 = 0.0
         zdkdt1    = (2.*pbcoef(jj,idepth)+3.)*zk(jj,idepth)/pwg(jj,idepth)
         zdkdt2    = 0.0                
!        Total Flux derrivative terms:
         zdfluxdt1(jj,idepth) = -zdkdt1*zhead(jj,idepth)*pfwtd(jj) - znu(jj,idepth)*zdheaddt1 - zdkdt1
         zdfluxdt2(jj,idepth) = -zdkdt2*zhead(jj,idepth)*pfwtd(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) = (pdzg(:,1)/ptstep) - zwght*zdfluxdt1(:,1)
zcmtrx(:,1) = -zwght*zdfluxdt2(:,1)
!
!Other sub-surface layers:       
DO jl=2,inl
   DO jj=1,ini   
      idepth=kwg_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) = (pdzg(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,kwg_layer,inl,zsol)
!
! 9. Final calculations and diagnostics:
!    -----------------------------------
!
!
DO jl=1,inl
   DO jj=1,ini
!   
      idepth=kwg_layer(jj)
      IF(jl<idepth)THEN
! 
!       Update liquid water content (m3 m-3):
        pwg(jj,jl)   = pwg(jj,jl) + zsol(jj,jl)    
!
!       Supersaturated drainage (kg m-2 s-1):
        zexcess(jj)  = max(0.0, pwg(jj,jl) - zwsat(jj,jl))
        pwg(jj,jl  ) = min(pwg(jj,jl),zwsat(jj,jl))
        pwg(jj,jl+1) = pwg(jj,jl+1) + zexcess(jj)*(pdzg(jj,jl)/pdzg(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):
        pwg(jj,idepth)   = pwg(jj,idepth) + zsol(jj,idepth)    
!        
!       Supersaturated drainage (kg m-2 s-1):
        zexcess(jj)    = max(0.0, pwg(jj,idepth) - zwsat(jj,idepth))
        pwg(jj,idepth) = min(pwg(jj,idepth),zwsat(jj,idepth))   
        pdrain(jj)    = pdrain(jj) + zexcess(jj)*pdzg(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,ini   
      idepth=kwg_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  - pwg(jj,jl))
        pwg(jj,jl)   = pwg(jj,jl)   + zexcess(jj) 
        pwg(jj,jl+1) = pwg(jj,jl+1) - zexcess(jj)*pdzg(jj,jl)/pdzg(jj,jl+1)
      ELSEIF(jl==idepth.AND.pwg(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,pwg(jj,idepth)-xwgmin)*pdzg(jj,idepth)*xrholw/ptstep
        pwg(jj,idepth) = xwgmin
      ENDIF
   ENDDO
ENDDO
!
IF (lhook) CALL dr_hook('HYDRO_SOILDIF',1,zhook_handle)
!
END SUBROUTINE hydro_soildif