hydro_soil.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_soil(HISBA,                                        &
                         ptstep,                                            &
                         pletr, pleg, ppg, pevapcor,                        &
                         pwdrain,                                           &
                         pc1, pc2, pc3, pc4b, pc4ref, pwgeq,                &
                         pd_g2, pd_g3, pwsat, pwfc,                         &
                         pdwgi1, pdwgi2, plegi, pd_g1, pcg, pct,            &
                         ptg, ptg2,                                         &
                         pwg1, pwg2, pwg3, pwgi1, pwgi2,                    &
                         prunoff, pdrain, hksat, pwwilt                     )  
!     #####################################################################
!
!!****  *HYDRO_SOIL*  
!!
!!    PURPOSE
!!    -------
!
!     Calculates the evolution of the water variables, i.e., the superficial
!     and deep-soil volumetric water content (wg and w2), the equivalent
!     liquid water retained in the vegetation canopy (Wr), the equivalent
!     water of the snow canopy (Ws), and also of the albedo and density of
!     the snow (i.e., ALBS and RHOS).  Also determine the runoff and drainage
!     into the soil.
!         
!     
!!**  METHOD
!!    ------
!
!!    EXTERNAL
!!    --------
!!
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Noilhan and Planton (1989)
!!    Belair (1995)
!!      
!!    AUTHOR
!!    ------
!!
!!      S. Belair           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!
!!      Original    14/03/95 
!!                  31/08/98 (V. Masson and F. Habets) add Dumenil et Todini
!!                           runoff scheme
!!                  31/08/98 (V. Masson and A. Boone) add the third soil-water
!!                           reservoir (WG3,D3)
!!                  15/03/99 A. Boone soil ice modification: advance
!!                           soil ice and liquid water, do budget checks
!!                  25/01/00 A. Boone : fully implicit method for WG2, WG3
!!                  05/03/07 A. Boone : changed drainage diagnostic computation for
!!                                      single bulk-soil option...i.e. for
!!                                      cases when HISBA=2-L or d2>=d3 (HISBA=3-L)
!!                                      for tighter water budget closure 
!!                  07/08/12 B. Decharme : Soil ice energy conservation
!!                     04/13 B. Decharme : Apply physical limits on wg here instead of in hydro.F90
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_csts,     ONLY : xlvtt, xrholw, xlmtt, xlstt, xday
USE modd_isba_par, ONLY : xwgmin
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
 CHARACTER(LEN=*),     INTENT(IN)   :: hisba   ! type of ISBA version:
!                                             ! '2-L' (default)
!                                             ! '3-L'
!
REAL, INTENT(IN)                  :: ptstep
!                                      timestep of the integration (s)
!
REAL, DIMENSION(:), INTENT(IN)    :: pd_g1
!                                      depth of surface ice reservoir (m)
!
REAL, DIMENSION(:), INTENT(IN)    :: pletr, pleg, ppg, pevapcor
!                                      PLETR    = evapotranspiration of the vegetation (W m-2)
!                                      PLEG     = latent heat of evaporation over the ground (W m-2)
!                                      PPG      = total water reaching the ground (kg m-2 s-1)
!                                      PEVAPCOR = correction for any excess evaporation 
!                                                from snow as it completely ablates (kg m-2 s-1)
!
REAL, DIMENSION(:), INTENT(IN)    :: pwdrain  ! minimum Wg for drainage (m3 m-3)
!
REAL, DIMENSION(:), INTENT(IN)    :: pc1, pc2, pwgeq, pcg, pct
REAL, DIMENSION(:,:), INTENT(IN)  :: pc3
!                                      soil coefficients
!                                      C1, C2 = coefficients for the moisture calculations (-)
!                                      C3     = coefficient for drainage calculation (m)
!                                      PWGEQ  = equilibrium surface volumetric moisture (m3 m-3)
!                                      PCG    = soil heat capacity
!                                      PCT    = grid-averaged heat capacity (K m2 J-1)
!
REAL, DIMENSION(:), INTENT(IN)    :: pd_g2, pd_g3, pwsat, pwfc
!                                      PD_G2 = root depth (m)
!                                      PD_G3 = depth of the soil column (m)
!                                      PWSAT = saturation volumetric water content
!                                              of the soil (m3 m-3)
!                                      PWFC  = field capacity volumetric water
!                                              content (m3 m-3)
!
REAL, DIMENSION(:), INTENT(IN)    :: pc4b, pc4ref
!                                      PC4REF, PC4B = fiiting soil paramter for vertical 
!                                      diffusion (C4) (-)
!
REAL, DIMENSION(:), INTENT(IN)    :: pdwgi1, pdwgi2, plegi
!                                      PDWGI1 = surface layer liquid water equivalent 
!                                               volumetric ice content time tendency (m3 m-3)
!                                      PDWGI2 = deep-soil layer liquid water equivalent  
!                                               volumetric ice content time tendency  (m3 m-3)
!                                      PLEGI  = surface soil ice sublimation (W m-2)
!
!
REAL, DIMENSION(:), INTENT(INOUT) :: ptg, ptg2
!                                    PTG  = surface temperature at 't' (K)
!                                    PTG2 = soil temperature at 't' (K)
!
REAL, DIMENSION(:), INTENT(INOUT) :: pwg1, pwg2, pwg3, pwgi1, pwgi2
!                                      PWG1   = near-surface soil moisture at 't+dt' (m3 m-3)
!                                      PWG2   = bulk root-soil moisture at 't+dt' (m3 m-3)
!                                      PWG3   = bulk deep-soil moisture at 't+dt' (m3 m-3)
!                                      PWGI1  = bulk surface-soil ice at 't+dt' (m3 m-3)
!                                      PWGI2  = bulk deep-soil ice at 't+dt' (m3 m-3)
!
REAL, DIMENSION(:), INTENT(OUT)   :: prunoff, pdrain
!                                      PRUNOFF = runoff (kg m-2 s-1)
!                                      PDRAIN  = drainage (kg m-2 s-1)
!
 CHARACTER(LEN=*),     INTENT(IN)  :: hksat      ! soil hydraulic profil option
!                                               ! 'DEF'  = ISBA homogenous soil
!                                               ! 'SGH'  = ksat exponential decay
!
REAL, DIMENSION(:), INTENT(IN)    :: pwwilt
!                                    PWWILT = wilting point volumetric water
!                                             content (m3 m-3)
!
!*      0.2    declarations of local variables
!
!
REAL, DIMENSION(SIZE(PTG))   :: zwgi1m, zwg2m, zwg3m, zwgi2m
!                                      Prognostic variables of ISBA at 't-dt'
!                                      ZWG2M = root-soil volumetric water content
!                                      ZWG3M = deep-soil volumetric water content
!                                      ZWGI1M = surface-soil volumetric ice content
!                                      ZWGI2M = deep-soil volumetric ice content
!
REAL, DIMENSION(SIZE(PTG))   :: zetr, zeg
!                                             ZETR = evapotranspiration rate
!                                             ZEG = evaporation rate from the ground
!
! 
REAL, DIMENSION(SIZE(PTG))   :: zwsat, zwfc
!                                             ZWSAT = Wsat  when ice is present
!                                             ZWFC  = Wfc   when ice is present
!
REAL, DIMENSION(SIZE(PWG3))  :: zc4
REAL, DIMENSION(SIZE(PWG3))  :: zwavg, zsink2,                                    &
                                  zfactor, zdraincf2, zdraincf3, zdrain2,           &
                                  zdelta2, zdelta3, zdelta22, zdelta33,             &
                                  zwdrain2, zwdrain3  
!
REAL, DIMENSION(SIZE(PWG3))  :: zexcessf, za2, zb2, zc2, za3, zb3, zc3, zwdrain, zexcessfc
!
REAL, DIMENSION(SIZE(PWG3))  :: zwlim2, zwlim3, zwwilt
!
INTEGER                 :: jj
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
!*       0.     Initialization
!               --------------
!
!
IF (lhook) CALL dr_hook('HYDRO_SOIL',0,zhook_handle)
zwsat(:)     = 0.
zwfc(:)      = 0.
zwwilt(:)    = 0.
!
pdrain(:)   = 0.
prunoff(:)   = 0.
!
zdrain2(:)   = 0.
zdraincf2(:) = 0.
zdraincf3(:) = 0.
zdelta2(:)   = 0.0
zdelta3(:)   = 0.0
zdelta22(:)  = 0.0
zdelta33(:)  = 0.0
zsink2(:)    = 0.
zwdrain(:)   = 0.0
zwdrain2(:)  = 0.0
zwdrain3(:)  = 0.0
za2(:)       = 0.0
zb2(:)       = 0.0
zc2(:)       = 0.0
za3(:)       = 0.0
zb3(:)       = 0.0
zc3(:)       = 0.0
!
zexcessf(:)  = 0.0
!
! Fields at time t-dt
!
zwg2m(:)     = pwg2(:)
zwg3m(:)     = pwg3(:)
zwgi1m(:)    = pwgi1(:)
zwgi2m(:)    = pwgi2(:)

!-------------------------------------------------------------------------------
!
DO jj=1,SIZE(ptg)
!
!*       1.     New Wsat
!               --------
!
  zwsat(jj) = pwsat(jj) - zwgi2m(jj)      ! per definition ZWGI2M<PWSAT-0.001
!                                             ! then ZWSAT>0.001
!
  zwfc(jj)   = pwfc(jj)   * zwsat(jj) / pwsat(jj)
!
  zwwilt(jj) = pwwilt(jj) * zwsat(jj) / pwsat(jj)
!
!
!                                           evaporation rates
!
  zetr(jj) = pletr(jj) / xlvtt
!
! Remove sublimation from total soil evaporation: this is done
! because ZEG is a liquid water sink term. Sublimation
! is removed from the surface ice tendency equation (below):
! PLEG represents the total of sublimation and evaporation.
! Also, add an additional possible correction which is used if
! evaporation from the snow surface exceeds the last amount of
! snow on the surface as the snow completely ablates (this
! is a very small term and is generally negligible HOWEVER
! it is retained here for a high order water budget):
!
  zeg(jj)  = pleg(jj) / xlvtt + pevapcor(jj)
!
!-------------------------------------------------------------------------------
!
!*       4.     EVOLUTION OF THE SUPERFICIAL WATER CONTENT WG1
!               ----------------------------------------------
!
!                                           updated values for wg 
!
  pwg1(jj) = (pwg1(jj) - ptstep * (pc1(jj) * (zeg(jj) - ppg(jj)) / xrholw    &
                                -  pc2(jj) * pwgeq(jj) / xday) )            &
              / (1. + ptstep * pc2(jj) / xday)  
!
!
ENDDO
!
IF(hksat=='SGH' .OR. hksat=='EXP') THEN
  zwlim2(:)=zwwilt(:)
  zwlim3(:)=pwwilt(:)
ELSE
  zwlim2(:)=xwgmin
  zwlim3(:)=xwgmin
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       5.     EVOLUTION OF THE DEEP WATER CONTENT WG2 and WG3
!               -----------------------------------------------
!
!*       5.1    2-L ISBA version
!               ----------------
!
IF (hisba=='2-L') THEN
! 
  DO jj=1,SIZE(ptg)

!
    pwg2(jj) = zwg2m(jj) - ptstep*(zeg(jj) + zetr(jj) - ppg(jj))   &
                  / (pd_g2(jj) * xrholw)  
!
!*       6.     DRAINAGE FROM THE DEEP SOIL
!               ------------------
!
    zwdrain(jj)   = pwdrain(jj) * max(0.0, min(zwfc(jj),pwg2(jj))-zwlim2(jj))/(zwfc(jj)-zwlim2(jj))

    zdrain2(jj)   =  max( min(zwdrain(jj),pwg2(jj)) , pwg2(jj)-zwfc(jj) )*pc3(jj,1)      &
                      / (pd_g2(jj)*xday) * ptstep  
!
!                                      the deep-soil volumetric water content w2
!                                      is modified consequently
!
    pwg2(jj)    = pwg2(jj) - zdrain2(jj)

    pdrain(jj)  = zdrain2(jj)*pd_g2(jj)*xrholw/ptstep  ! Final output units: kg m-2 s-1
!
  ENDDO
!
ELSE
!
  DO jj=1,SIZE(ptg)

!*       5.2    3-L ISBA version (with only 2 active layers)
!               ----------------
!
    IF (pd_g2(jj) >= pd_g3(jj)) THEN

      pwg2(jj) = zwg2m(jj) - ptstep*(zeg(jj) + zetr(jj) - ppg(jj))   &
                    / (pd_g2(jj) * xrholw)  

!*       6.     DRAINAGE FROM THE DEEP SOIL
!               ------------------
!                                      when w2 > wfc, there is drainage
!
      zwdrain(jj) = pwdrain(jj) * max(0.0, min(zwfc(jj),pwg2(jj))-zwlim2(jj))/(zwfc(jj)-zwlim2(jj))

      zdrain2(jj) = max( min(zwdrain(jj),pwg2(jj)) , pwg2(jj)-zwfc(jj) )*pc3(jj,1)        &
                   / (pd_g2(jj)*xday) * ptstep  

      pwg2(jj)    = pwg2(jj) -  zdrain2(jj)
      pwg3(jj)    = pwg2(jj)

      pdrain(jj)  = zdrain2(jj)*pd_g2(jj)*xrholw/ptstep  ! Final output units: kg m-2 s-1
!
    ELSE
!
!*       5.3    3-L ISBA version (with 3 active layers)
!               ----------------
!
! Linear drainage term check (m3 m-3):
! Use a linear scaling to prevent it from completely drying out the soil. This
! term will, under most conditions, not result in such a drying. Until a more
! rhobust (but more complicated) method is developed for maintaining a minimum river
! flow under dry conditions, this method will be used.
!
      zwdrain2(jj)   = pwdrain(jj)* max(0.0, min(zwfc(jj),zwg2m(jj))-zwlim2(jj))/(zwfc(jj)-zwlim2(jj))
      zwdrain3(jj)   = pwdrain(jj)* max(0.0, min(pwfc(jj),zwg3m(jj))-zwlim3(jj))/(pwfc(jj)-zwlim3(jj))
!
! Delta functions:
!
      zdelta2(jj)    = 0.0
      IF ( zwg2m(jj) - zwfc(jj) > zwdrain2(jj) ) zdelta2(jj) = 1.0
  !
      zdelta3(jj)    = 0.0
      IF ( zwg3m(jj) - pwfc(jj) > zwdrain3(jj) ) zdelta3(jj) = 1.0
  !
  !

! evaluate inter-facial water content, grid factor, and diffusion coefficient:

      zwavg(jj)     = ( ( (zwg2m(jj)**6)* pd_g2(jj)          +                       &
                           (zwg3m(jj)**6)*(pd_g3(jj)-pd_g2(jj)) )/pd_g3(jj) )**(1./6.)  

      zfactor(jj)   = pd_g2(jj)/(pd_g3(jj)-pd_g2(jj))

      zc4(jj)   = pc4ref(jj)*(zwavg(jj)**pc4b(jj))                &
                     *(10.**(-pc4b(jj)*pwgi2(jj)/(pwsat(jj)-xwgmin)))  
!
! calculate sources/sinks
!
      zsink2(jj)   = -(zeg(jj) + zetr(jj) - ppg(jj) )/(pd_g2(jj)*xrholw)

! Compute evolution of water content using linearized equations
! (see Boone 2000, Appendix F.2 for details)
!
! sink terms are treated explicitly, other terms are implicit
!
      zdraincf2(jj) = pc3(jj,1) / (pd_g2(jj) * xday)
      zdelta22(jj)  = zdelta2(jj)*zwfc(jj) - (1.0-zdelta2(jj))*zwdrain2(jj)
      zc2(jj)       = 1.0 + ptstep*(zdelta2(jj)*zdraincf2(jj) + (zc4(jj)/xday)) 
      zb2(jj)       = ptstep*zc4(jj)/(xday*zc2(jj))
      za2(jj)       = ( zwg2m(jj) + ptstep*(zsink2(jj) + zdraincf2(jj)*zdelta22(jj)) )/zc2(jj)
!
      zdraincf3(jj) = pc3(jj,2) / ( (pd_g3(jj)-pd_g2(jj)) * xday)
      zdelta33(jj)  = zdelta3(jj)*pwfc(jj) - (1.0-zdelta3(jj))*zwdrain3(jj)
      zc3(jj)       = 1.0 + ptstep*(zdelta3(jj)*zdraincf3(jj) + zfactor(jj)*(zc4(jj)/xday)) 
      zb3(jj)       = ptstep*zfactor(jj)*(zdelta2(jj)*zdraincf2(jj) + (zc4(jj)/xday) )/zc3(jj)
      za3(jj)       = ( zwg3m(jj) + ptstep*(                                                 &
                         - zfactor(jj)*zdraincf2(jj)*zdelta22(jj)                               &
                         +            zdraincf3(jj)*zdelta33(jj)) )/zc3(jj)  
!
! Advance volumetric water content values in time:
! system of 2 linear equations:
!
      pwg2(jj)      = ( za2(jj)+zb2(jj)*za3(jj) )/(1.0 - zb2(jj)*zb3(jj))
      pwg3(jj)      = za3(jj) + zb3(jj)*pwg2(jj)
!
! Drainage (kg m-2 s-1): this term is implicit and is extracted directly from
!                        the drainage computation in the above equations.
!
      zwdrain(jj)   = (xrholw*pc3(jj,2)/xday)*                                                 &
                       ( zdelta3(jj)*(pwg3(jj)-pwfc(jj)) + (1.0-zdelta3(jj))*zwdrain3(jj) )  
!
! As drainage is implicit, perform a check to prevent any negative drainage
! (can arise rarely and is generally negligible, but to ensure a high order conservation):
!
      pdrain(jj)    = max(0.0, zwdrain(jj))
      pwg3(jj)      = pwg3(jj) + (pdrain(jj) - zwdrain(jj))*ptstep/((pd_g3(jj)-pd_g2(jj))*xrholw)
!
    ENDIF
  ENDDO
END IF
!
!-------------------------------------------------------------------------------
!
DO jj=1,SIZE(ptg)
!*       7.     EFFECT OF THE MELTING/FREEZING ON THE SOIL WATER CONTENT
!               --------------------------------------------------------
!
!*       7.1    Effect on surface water liquid and ice reservoirs
!               -------------------------------------------------
!
! First, advance the surface layer ice content using the time tendency
! and sublimation (*heat effect* of sublimation already accounted
! for in latent heat flux calculation):
!           
  pwgi1(jj) = zwgi1m(jj) + pdwgi1(jj) - plegi(jj)*ptstep/(xlstt*pd_g1(jj)*xrholw)
!
! Next, update the liquid water content:
!
  pwg1(jj)  = pwg1(jj)  - pdwgi1(jj) 
!
! Make sure that ice has not dropped below
! zero due to sublimation (the ONLY way
! it can drop below 0). If it has, then
! freeze the needed liquid water (thus warming the
! layer) to ensure sublimation is accomodated,
! then extract this from liquid water store.
! This is a correction used ONLY when last traces
! of soil ice sublimate within a model time step.
!
  zexcessfc(jj)= 0.0
!
  zexcessf(jj) = max(0.0, - pwgi1(jj))
  pwg1(jj)     = pwg1(jj)  - zexcessf(jj)
  pwgi1(jj)    = pwgi1(jj) + zexcessf(jj)
  zexcessfc(jj)= zexcessfc(jj) - zexcessf(jj)
!
! Modif H.Douville 26/08/03 (global scale)
! Make sure that ice has not raised above 
! saturation (minus XWGMIN) due to sublimation.
! If it has, then melt the needed frozen water 
! (thus cooling the layer) to ensure sublimation
! is accomodated, then extract this from frozen water store.
!
  zexcessf(jj) = min(0.0, pwsat(jj) - xwgmin - pwgi1(jj))
  pwg1(jj)     = pwg1(jj)  - zexcessf(jj)
  pwgi1(jj)    = pwgi1(jj) + zexcessf(jj)
  zexcessfc(jj)= zexcessfc(jj) - zexcessf(jj)
!
! Make sure that liquid has not dropped below
! minimum threshold: this could arise due to
! evaporation and ice formation. Melt the needed water.
! Normally simply a budget check, i.e. usually small but accounted
! for none-the-less to assure high accuracy.
!
  zexcessf(jj) = max(0.0, xwgmin - pwg1(jj))
  pwgi1(jj)    = pwgi1(jj)  - zexcessf(jj)
  pwg1(jj)     = pwg1(jj)   + zexcessf(jj)
  zexcessfc(jj)= zexcessfc(jj) + zexcessf(jj)
!
! removes very small values due to computation precision
!
  IF(pwgi1(jj) < 1.0e-10) THEN
    zexcessf(jj) = pwgi1(jj)
    pwg1(jj)     = pwg1(jj) + zexcessf(jj)
    pwgi1(jj)    = 0.0
    zexcessfc(jj)= zexcessfc(jj) + zexcessf(jj)
  ENDIF
!
! Cummulative phase change for the ice/liquid budget corrections:
!
  ptg(jj)        = ptg(jj) - zexcessfc(jj)*xlmtt*pct(jj)*xrholw*pd_g1(jj)
!
!
!*       7.2    Effect on deep-soil liquid and ice reservoirs
!               ---------------------------------------------
!
! Update the ice content using the ice tendency:
! Since this reservoir includes surface reservoir, add
! any changes in ice content due to sublimation:
!
  pwgi2(jj) = zwgi2m(jj) + pdwgi2(jj) - plegi(jj)*ptstep/(xlstt*pd_g2(jj)*xrholw)
!
! Update the liquid water content:
!
  pwg2(jj)   = pwg2(jj)   - pdwgi2(jj)
!
! Make sure that ice has not dropped below
! zero due to sublimation (as above).
!
  zexcessfc(jj)= 0.0
!
  zexcessf(jj) = max(0.0, -pwgi2(jj))
  pwg2(jj)     = pwg2(jj)  - zexcessf(jj)
  pwgi2(jj)    = pwgi2(jj) + zexcessf(jj)
  zexcessfc(jj)= zexcessfc(jj) - zexcessf(jj)
!
! Budget check of minimum threshold for liquid
! water as for surface: MUCH LESS likely
! to be utilized, but retained for accuracy
! in energy and water balance (as above).
!
  zexcessf(jj) = max(0.0, xwgmin - pwg2(jj))
  pwgi2(jj)    = pwgi2(jj)  - zexcessf(jj)
  pwg2(jj)     = pwg2(jj)   + zexcessf(jj)
  zexcessfc(jj)= zexcessfc(jj) + zexcessf(jj)
!
! removes very small values due to computation precision
!
  IF (pwgi2(jj) < 1.0e-10 * ptstep) THEN
      zexcessf(jj) = pwgi2(jj)
      pwg2(jj)    = pwg2(jj) + zexcessf(jj)
      pwgi2(jj)    = 0.
      zexcessfc(jj)= zexcessfc(jj) + zexcessf(jj)
  ENDIF
!
! Cummulative phase change for the ice/liquid budget corrections:
!
  ptg2(jj) = ptg2(jj) - zexcessfc(jj)*xlmtt*pcg(jj)*xrholw*pd_g2(jj)
!
!
ENDDO
!
!-------------------------------------------------------------------------------
!
!*       8.     PHYSICAL LIMITS AND RUNOFF
!               --------------------------
!
! runoff of second layer
!
prunoff(:) = max( 0., pwg2(:)+pwgi2(:)-pwsat(:) )*pd_g2(:) * xrholw / ptstep
!
! now apply limits:
!
pwg1(:) = min( pwg1(:), pwsat(:) - pwgi1(:) )
pwg1(:) = max( pwg1(:), xwgmin              )
!
pwg2(:) = min( pwg2(:), pwsat(:) - pwgi2(:) )
pwg2(:) = max( pwg2(:), xwgmin              )
!
!runoff of third layer added to drainage
!
IF (hisba=='3-L') THEN
   pdrain(:) = pdrain(:) + max( 0., pwg3(:)-pwsat(:) )* (pd_g3(:)-pd_g2(:)) * xrholw / ptstep  
   pwg3(:) = min( pwg3(:), pwsat(:)         )
   pwg3(:) = max( pwg3(:), xwgmin           )
END IF
!
IF (lhook) CALL dr_hook('HYDRO_SOIL',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE hydro_soil