hydro.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(HISBA, HSNOW_ISBA, HRUNOFF, HSOILFRZ, OMEB, OGLACIER,&
                         oflood, ptstep, pvegtype,                          &
                         prr, psr, plev, pletr, pleg, ples,                 &
                         prunoffb, pwdrain,                                 &
                         pc1, pc2, pc3, pc4b, pc4ref, pwgeq, pcg, pct,      &
                         pveg, plai, pwrmax, pmelt, ptauice, plegi,         &
                         prunoffd, psoilwght, klayer_hort, klayer_dun,      &
                         ppsnv, ppsng,                                      &
                         psnow_thrufal, pevapcor, psubvcor,                 &
                         pwr, psoilhcapz,                                   &
                         psnowswe, psnowalb, psnowrho,                      &
                         pbcoef, pwsat, pcondsat, pmpotsat, pwfc,           &
                         pwwilt, pf2wght, pf2, pd_g, pdzg, pdzdif, pps,     &
                         pwg, pwgi, ptg, kwg_layer,                         &
                         pdrain, prunoff, ptopqs,                           &
                         pirrig, pwatsup, pthreshold, lirriday, lirrigate,  &
                         hksat, hrain, hhort, pmuf, pfsat, pksat_ice,       &
                         pd_ice, phorton, pdrip, pffg, pffv , pfflood,      &
                         ppiflood, piflood, ppflood, prrveg, pirrig_flux,   &
                         pirrig_gr, pqsb, pfwtd, pwtd,                      &
                         pdelheatg, pdelheatg_sfc,                          &
                         pdelphaseg, pdelphaseg_sfc, plvtt, plstt           )
!     #####################################################################
!
!!****  *HYDRO*  
!!
!!    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., SNOWALB and SNOWRHO).  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)
!!                  19/07/05 (P. LeMoigne) bug in runoff computation if isba-2L
!!                  10/10/05 (P. LeMoigne) bug in hydro-soil calling sequence
!!                  25/05/08 (B. Decharme) Add floodplains
!!                  27/11/09 (A. Boone)    Add possibility to do time-splitting when
!!                                         calling hydro_soildif (DIF option only)
!!                                         for *very* large time steps (30min to 1h+).
!!                                         For *usual* sized time steps, time step
!!                                         NOT split.
!!                     08/11 (B. Decharme) DIF optimization
!!                     09/12 (B. Decharme) Bug in wg2 ice energy budget
!!                     10/12 (B. Decharme) EVAPCOR snow correction in DIF
!!                                         Add diag IRRIG_FLUX
!!                     04/13 (B. Decharme) Pass soil phase changes routines here
!!                                         Apply physical limits on wg in hydro_soil.F90
!!                                         Subsurface runoff if SGH (DIF option only)
!!                                         water table / surface coupling
!!                  02/2013  (C. de Munck) specified irrigation rate of ground added
!!                  10/2014  (A. Boone)    MEB added
!!                  07/15    (B. Decharme) Numerical adjustement for F2 soilstress function
!!                  03/16    (B. Decharme) Limit flood infiltration
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_csts,      ONLY : xrholw, xday, xtt, xlstt, xlmtt
USE modd_isba_par,  ONLY : xwgmin, xdenom_min
USE modd_surf_par,  ONLY : xundef, nundef
!
#ifdef TOPD
USE modd_coupling_topd, ONLY : lcoupl_topd, xas_nature, xatop, xrunoff_top, nmaskt_patch
#endif
!
USE modi_hydro_veg
USE modi_hydro_snow
USE modi_hydro_soil
USE modi_hydro_soildif                                          
USE modi_hydro_sgh
USE modi_ice_soildif              
USE modi_ice_soilfr              
!
USE mode_thermos
!
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'
!                                             ! 'DIF'   ISBA-DF
 CHARACTER(LEN=*),     INTENT(IN)  :: hsnow_isba ! 'DEF' = Default F-R snow scheme
!                                               !         (Douville et al. 1995)
!                                               ! '3-L' = 3-L snow scheme (option)
!                                               !         (Boone and Etchevers 2001)
 CHARACTER(LEN=*),     INTENT(IN)   :: hrunoff ! surface runoff formulation
!                                             ! 'WSAT'
!                                             ! 'DT92'
!                                             ! 'SGH ' Topmodel
 CHARACTER(LEN=*),   INTENT(IN)      :: hsoilfrz   ! soil freezing-physics option
!                                                 ! 'DEF'   Default (Boone et al. 2000; Giard and Bazile 2000)
!                                                 ! 'LWT'   phase changes as above, but relation between unfrozen 
!                                                         water and temperature considered
LOGICAL, INTENT(IN)                :: oglacier ! True = Over permanent snow and ice, 
!                                                initialise WGI=WSAT,
!                                                Hsnow>=10m and allow 0.8<SNOALB<0.85
!                                                False = No specific treatment
!
LOGICAL, INTENT(IN)                :: omeb   ! True  = patch with multi-energy balance 
!                                            ! False = patch with classical (composite) ISBA 
!
LOGICAL, INTENT(IN)                :: oflood ! Flood scheme 
!
REAL, INTENT(IN)                    :: ptstep
!                                      timestep of the integration
!
REAL, DIMENSION(:,:), INTENT(IN)  :: pvegtype ! fraction of each vegetation
!
REAL, DIMENSION(:), INTENT(IN)    :: prr, psr, plev, pletr, pleg, ples
!                                      PRR = rain rate
!                                      PSR = snow rate
!                                      PLEV = latent heat of evaporation over vegetation
!                                      PLETR = evapotranspiration of the vegetation
!                                      PLEG = latent heat of evaporation over the ground
!                                      PLES = latent heat of sublimation over the snow
!
REAL, DIMENSION(:), INTENT(IN)    :: prunoffb      ! slope of the runoff curve
REAL, DIMENSION(:), INTENT(IN)    :: plvtt, plstt  ! latent heat of vaporization and sublimation (J/kg)  
REAL, DIMENSION(:), INTENT(IN)    :: pwdrain       ! minimum Wg for drainage (m3/m3)
!
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 WG2 calculation
!                                      PWGEQ = equilibrium surface volumetric moisture
!                                      PCG = soil heat capacity
!                                      PCT = grid-averaged heat capacity
!
REAL, DIMENSION(:), INTENT(IN)    :: pveg, plai, prunoffd, pwrmax
!                                      PVEG = fraction of vegetation 
!                                      PLAI = Leaf Area Index 
!                                      PRUNOFFD = depth over which sub-grid runoff calculated (m)
!                                      PWRMAX = maximum equivalent water content
!                                               in the vegetation canopy
!
REAL, DIMENSION(:), INTENT(IN)    :: pc4b, pc4ref
!                                      PC4REF, PC4B = fiiting soil paramter for vertical diffusion (C4)
!
REAL, DIMENSION(:), INTENT(IN)    :: ppsnv, ppsng
!                                      PPSNV = vegetation covered by snow
!                                      PPSNG = baresoil covered by snow
!
REAL, DIMENSION(:), INTENT(IN)    :: ptauice, plegi
!                                    PTAUICE = characteristic time scale for soil water
!                                                phase changes (s) 
!                                    PLEGI   = surface soil ice sublimation 
!
REAL, DIMENSION(:), INTENT(IN)    :: psnow_thrufal, pevapcor, psubvcor
!                                    PSNOW_THRUFAL = rate that liquid water leaves snow pack: 
!                                               *ISBA-ES* [kg/(m2 s)]
!                                    PEVAPCOR = correction if evaporation from snow exceeds
!                                               actual amount on the surface [kg/(m2 s)]
!                                    PSUBVCOR = correction if sublimation from snow intercepted 
!                                               on the MEB canopy exceeds snow available as it 
!                                               disappears [kg/(m2 s)]
!
REAL, DIMENSION(:), INTENT(IN)    :: pps, pf2                                       
!                                    PPS  = surface pressure (Pa)
!                                    PF2  = total water stress factor (-)
!
REAL, DIMENSION(:,:), INTENT(IN)  :: pwsat, pcondsat, pwfc, pd_g, pf2wght, pwwilt
!                                    PD_G   = Depth of bottom of Soil layers (m)
!                                    PWFC     = field capacity profile (m3/m3)
!                                    PWWILT   = wilting point profile (m3/m3)
!                                    PWSAT    = porosity profile (m3/m3)
!                                    PCONDSAT = hydraulic conductivity at saturation (m/s)
!                                    PF2WGHT   = water stress factor (profile) (-)

REAL, DIMENSION(:,:), INTENT(IN)  :: pdzdif, pdzg
!                                    PDZDIF = distance between consecuative layer mid-points
!                                    PDZG   = soil layers thicknesses
!
REAL, DIMENSION(:,:), INTENT(IN) :: psoilhcapz
!                                   PSOILHCAPZ = ISBA-DF Soil heat capacity profile [J/(m3 K)]
!
REAL, DIMENSION(:,:), INTENT(IN)  :: psoilwght  ! ISBA-DIF: weights for vertical
!                                               ! integration of soil water and properties
INTEGER,             INTENT(IN)   :: klayer_hort! DIF optimization
INTEGER,             INTENT(IN)   :: klayer_dun ! DIF optimization
!
REAL, DIMENSION(:,:), INTENT(IN)  :: pmpotsat,pbcoef
!                                    PMPOTSAT = matric potential at saturation (m)
!                                    PBCOEF   = slope of the retention curve (-)
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: ptg
!                                    PTG   = soil layer average temperatures (K)
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: pwgi, pwg
!                                      PWGI  = soil frozen volumetric water content (m3/m3)
!                                      PWG  = soil liquid volumetric water content (m3/m3)
!                                      Prognostic variables of ISBA at 't-dt'
!                                      PWGI(:,1) = surface-soil volumetric ice content
!                                      PWGI(:,2) = deep-soil volumetric ice content
!
INTEGER, DIMENSION(:), INTENT(IN) :: kwg_layer  
!                                    KWG_LAYER = Number of soil moisture layers (DIF option)
!
REAL, DIMENSION(:), INTENT(INOUT) :: pdelheatg, pdelheatg_sfc
!                                      PDELHEATG_SFC = change in heat storage of the surface soil layer over the current time step (W m-2)
!                                      PDELHEATG     = change in heat storage of the entire soil column over the current time step (W m-2)
!
REAL, DIMENSION(:), INTENT(INOUT) :: pwr, psnowswe, psnowalb, psnowrho
REAL, DIMENSION(:), INTENT(INOUT) :: pmelt
REAL, DIMENSION(:), INTENT(OUT)   :: pdrain, prunoff
!                                      PWR = liquid water retained on the foliage
!                                             of the vegetation at time 't+dt'
!                                      PSNOWSWE = equivalent water content of the
!                                             snow reservoir at time 't+dt'
!                                      PSNOWALB = albedo of the snow at 't+dt'
!                                      PSNOWRHO = density of the snow at 't+dt' 
!                                      PMELT = melting rate of the snow
!
REAL, DIMENSION(:), INTENT(OUT)   :: pdelphaseg, pdelphaseg_sfc
!                                     PDELPHASEG     = latent heating due to soil freeze-thaw in the entire soil column  (W m-2)
!                                     PDELPHASEG_SFC = latent heating due to soil freeze-thaw in the surface soil layer  (W m-2)
!
!
REAL   ,DIMENSION(:),INTENT(IN)    :: pirrig
REAL   ,DIMENSION(:),INTENT(IN)    :: pwatsup
REAL   ,DIMENSION(:),INTENT(IN)    :: pthreshold
LOGICAL,DIMENSION(:),INTENT(INOUT) :: lirriday
LOGICAL,DIMENSION(:),INTENT(IN)    :: lirrigate
REAL   ,DIMENSION(:),INTENT(INOUT) :: pirrig_flux ! irrigation rate (kg/m2/s)
REAL   ,DIMENSION(:),INTENT(IN)    :: pirrig_gr ! ground irrigation rate (kg/m2/s)
!
!
!
 CHARACTER(LEN=*),     INTENT(IN)   :: hksat   ! soil hydraulic profil option
!                                             ! 'DEF'  = ISBA homogenous soil
!                                             ! 'SGH'  = ksat exponential decay
!
 CHARACTER(LEN=*), INTENT(IN)       :: hrain   ! Rainfall spatial distribution
                                              ! 'DEF' = No rainfall spatial distribution
                                              ! 'SGH' = Rainfall exponential spatial distribution
                                              ! 
!
 CHARACTER(LEN=*), INTENT(IN)       :: hhort   ! Horton runoff
                                              ! 'DEF' = no Horton runoff
                                              ! 'SGH' = Horton runoff
!                                        
REAL, DIMENSION(:),  INTENT(IN)   :: pd_ice   !depth of the soil column for the calculation
!                                              of the frozen soil fraction (m)
REAL, DIMENSION(:),  INTENT(IN)   :: pksat_ice!hydraulic conductivity at saturation (m/s)
!                                            
REAL, DIMENSION(:),  INTENT(IN)   :: pmuf     !fraction of the grid cell reached by the rainfall
REAL, DIMENSION(:),  INTENT(INOUT):: pfsat    !Topmodel/dt92 saturated fraction
!
REAL, DIMENSION(:),  INTENT(OUT)  :: phorton   !Horton runoff (kg/m2/s)
!
REAL, DIMENSION(:),  INTENT(INOUT) :: pdrip    !Dripping from the vegetation (kg/m2/s)
REAL, DIMENSION(:),  INTENT(INOUT) :: prrveg   !Precip. intercepted by vegetation (kg/m2/s)
!
REAL, DIMENSION(:),  INTENT(IN)    :: pffg,pffv
REAL, DIMENSION(:),  INTENT(IN)    :: pfflood  !Floodplain effective fraction
REAL, DIMENSION(:),  INTENT(IN)    :: ppiflood !Floodplain potential infiltration [kg/m²/s]
!
REAL, DIMENSION(:), INTENT(INOUT)  :: piflood  !Floodplain real infiltration      [kg/m²/s]
REAL, DIMENSION(:), INTENT(INOUT)  :: ppflood  !Floodplain interception           [kg/m²/s]
!
REAL, DIMENSION(:,:), INTENT(IN)   :: ptopqs   !Topmodel (HRUNOFF=SGH) subsurface flow by layer (m/s)
REAL, DIMENSION(:),   INTENT(OUT)  :: pqsb     !Topmodel (HRUNOFF=SGH) lateral subsurface flow [kg/m²/s]
!
REAL, DIMENSION(:), INTENT(IN)     :: pfwtd    !grid-cell fraction of water table to rise
REAL, DIMENSION(:), INTENT(IN)     :: pwtd     !water table depth (m)
!
!*      0.2    declarations of local variables
!
!
INTEGER                         :: jj, jl      ! loop control                                       
INTEGER                         :: indt, jdt   ! Time splitting indicies
INTEGER                         :: ini, inl, idepth ! (ISBA-DF option)
!
REAL                            :: ztstep      ! maximum time split time step (<= PTSTEP)
!                                              ! ONLY used for DIF option.
!
REAL, DIMENSION(SIZE(PVEG))     :: zpg, zpg_melt, zdunne,                            &
                                   zlev, zleg, zlegi, zletr, zpsnv,                  &
                                   zrr, zdg3, zwg3, zwsat_avg, zwwilt_avg, zwfc_avg, &
                                   zrunoff, zdrain, zhorton, zevapcor, zqsb 
!                                      Prognostic variables of ISBA at 't-dt'
!                                      ZPG = total water reaching the ground
!                                      ZPG_MELT = snowmelt reaching the ground 
!                                      ZDUNNE  = Dunne runoff
!                                 ZLEV, ZLEG, ZLEGI, ZLETR = Evapotranspiration amounts
!                                      from the non-explicit snow area *ISBA-ES*
!                                 ZPSNV = used to calculate interception of liquid
!                                      water by the vegetation in FR snow method:
!                                      For ES snow method, precipitation already modified
!                                      so set this to zero here for this option.
!                                 ZWSAT_AVG, ZWWILT_AVG, ZWFC_AVG = Average water and ice content
!                                      values over the soil depth D2 (for calculating surface runoff)
!                                 ZDG3, ZWG3, ZRUNOFF, ZDRAIN, ZQSB and ZHORTON are working variables only used for DIF option
!                                 ZEVAPCOR = correction if evaporation from snow exceeds
!                                               actual amount on the surface [m/s]
!
REAL, DIMENSION(SIZE(PVEG))     :: zdwgi1, zdwgi2, zksfc_iveg
!                                      ZDWGI1 = surface layer liquid water equivalent 
!                                               volumetric ice content time tendency
!                                      ZDWGI2 = deep-soil layer liquid water equivalent 
!                                               volumetric ice content time tendency
!                                      ZKSFC_IVEG = non-dimensional vegetation insolation coefficient
!
REAL, DIMENSION(SIZE(PVEG))    :: zwgi_excess, zf2
!                                 ZWGI_EXCESS = Soil ice excess water content
!                                 ZF2         = Soilstress function for transpiration
!
REAL, DIMENSION(SIZE(PWG,1),SIZE(PWG,2)) :: zqsat, zqsati, zti, zps
!                                           For specific humidity at saturation computation (ISBA-DIF)
!
REAL, DIMENSION(SIZE(PWG,1),SIZE(PWG,2)) :: zwgi0
!                                      ZWGI0 = initial soil ice content (m3 m-3) before update
!                                              for budget diagnostics     
!
!*      0.3    declarations of local parameters
!
REAL, PARAMETER             :: zinsolfrz_veg = 0.20  ! (-)       Vegetation insolation coefficient
!
REAL, PARAMETER             :: zinsolfrz_lai = 30.0  ! (m2 m-2)  Vegetation insolation coefficient

REAL, PARAMETER             :: ztimemax      = 900.  ! s  Maximum timescale without time spliting
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
!*       0.     Initialization:
!               ---------------
!
IF (lhook) CALL dr_hook('HYDRO',0,zhook_handle)
!
jdt    = 0
indt   = 0
ztstep = 0.0
!
zpg(:)           = 0.0
zpg_melt(:)      = 0.0
zdunne(:)        = 0.0
!
zwsat_avg(:)     = 0.0
zwwilt_avg(:)    = 0.0
zwfc_avg(:)      = 0.0
!
zrr(:)           = prr(:)
!
zdrain(:)        = 0.
zhorton(:)       = 0.
zrunoff(:)       = 0.
zwgi_excess(:)   = 0.
zevapcor(:)      = 0.
zqsb(:)      = 0.
!
pdrain(:)        = 0.
prunoff(:)       = 0.
phorton(:)       = 0.
pqsb(:)       = 0.
!
pdelphaseg(:)    = 0.0
pdelphaseg_sfc(:)= 0.0
zwgi0(:,:)       = 0.0
!
zf2(:) = max(xdenom_min,pf2(:))
!
! Initialize evaporation components: variable definitions
! depend on snow or explicit canopy scheme:
!
IF(omeb)THEN
!
! MEB uses explicit snow scheme by default, but fluxes already aggregated
! for snow and floods so no need to multiply by fractions here. 
!
   zlev(:)          = plev(:)
   zletr(:)         = pletr(:)
   zleg(:)          = pleg(:)
   zlegi(:)         = plegi(:)
   zpsnv(:)         = 0.0
!
   zevapcor(:)      = pevapcor(:) + psubvcor(:)
!
ELSE
!
! Initialize evaporation components: variable definitions
! depend on snow scheme:
!
   IF(hsnow_isba == '3-L' .OR. hsnow_isba == 'CRO' .OR. hisba == 'DIF')THEN
      zlev(:)          = (1.0-ppsnv(:)-pffv(:)) * plev(:)
      zletr(:)         = (1.0-ppsnv(:)-pffv(:)) * pletr(:)
      zleg(:)          = (1.0-ppsng(:)-pffg(:)) * pleg(:)
      zlegi(:)         = (1.0-ppsng(:)-pffg(:)) * plegi(:)
      zpsnv(:)         = 0.0
   ELSE
      zlev(:)          = plev(:)
      zletr(:)         = pletr(:)
      zleg(:)          = pleg(:)
      zlegi(:)         = plegi(:)
      zpsnv(:)         = ppsnv(:)+pffv(:)
   ENDIF
!
   zevapcor(:)         = pevapcor(:) 

ENDIF
!
! Initialize average soil hydrological parameters
! over the entire soil column: if Isba Force-Restore
! is in use, then parameter profile is constant
! so simply use first element of this array: if
! the Diffusion option is in force, the relevant
! calculation is done later within this routine.
!
IF(hisba == '2-L' .OR. hisba == '3-L')THEN  
   zwsat_avg(:)     = pwsat(:,1)
   zwwilt_avg(:)    = pwwilt(:,1)
   zwfc_avg(:)      = pwfc(:,1)
ENDIF
!
IF (hisba == '3-L') THEN                                   
   zdg3(:) = pd_g(:,3)
   zwg3(:) = pwg(:,3)
ELSE
   zdg3(:) = xundef
   zwg3(:) = xundef
END IF
!
!-------------------------------------------------------------------------------
!
!*       1.     EVOLUTION OF THE EQUIVALENT WATER CONTENT Wr
!               --------------------------------------------
!
!
!
IF(.NOT.omeb)THEN ! Canopy Int & Irrig Already accounted for if MEB in use.
!
   pirrig_flux(:)=0.0
!
!* add irrigation over vegetation to liquid precipitation (rr)
!
!
   IF (SIZE(lirrigate)>0) THEN
      WHERE (lirrigate(:) .AND. pirrig(:)>0. .AND. pirrig(:) /= xundef .AND. (pf2(:)<pthreshold(:)) )
         pirrig_flux(:) = pwatsup(:) / xday           
         zrr(:) = zrr(:) + pwatsup(:) / xday
         lirriday(:) = .true.           
      END WHERE
   ENDIF
!
!* interception reservoir and dripping computation
!
   CALL hydro_veg(hrain, ptstep, pmuf,                    &
                   zrr, zlev, zletr, pveg, zpsnv,         &
                   pwr, pwrmax, zpg, pdrip, prrveg, plvtt ) 
!
!
!
ELSE
!
! For MEB case, interception interactions already computed and PRR represents
! water falling (drip and not intercepted by vegetation) outside of snow covered
! areas. Part for snow covered areas (net outflow at base of snowpack) accounted
! for in PSNOW_THRUFAL.
!
   zpg(:) = prr(:)
!
ENDIF
!
!* add irrigation over ground to potential soil infiltration (pg)
!
pirrig_flux(:) = pirrig_flux(:) + pirrig_gr(:)
!
zpg(:) = zpg(:) + pirrig_gr(:)
!
!-------------------------------------------------------------------------------
!
!*       2.     EVOLUTION OF THE EQUIVALENT WATER CONTENT snowSWE 
!               -------------------------------------------------
!
!*       3.     EVOLUTION OF SNOW ALBEDO 
!               ------------------------
!
!*       4.     EVOLUTION OF SNOW DENSITY 
!               -------------------------
!
! Boone and Etchevers '3-L' snow option
IF(hsnow_isba == '3-L' .OR. hsnow_isba == 'CRO' .OR. hisba == 'DIF')THEN
!
  zpg_melt(:)   = zpg_melt(:)   + psnow_thrufal(:)          ! [kg/(m2 s)]
!
! Note that 'melt' is referred to as rain and meltwater
! running off from the snowpack in a timestep for ISBA-ES,
! not the actual amount of ice converted to liquid.
!
  pmelt(:) = pmelt(:) + psnow_thrufal(:)          ! [kg/(m2 s)]
!
ELSE
  !
  CALL hydro_snow(oglacier, ptstep, pvegtype,          &
                  psr, ples, pmelt,                    &
                  psnowswe, psnowalb, psnowrho,zpg_melt)
  !
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       5.     Sub Grid Hydrology
!               ------------------
!
! - Dunne runoff  : Dumenil et Todini (1992) or Topmodel
! - Horton runoff : Direct or exponential precipitation distribution
! - Floodplains interception and infiltration
!
 CALL hydro_sgh(hisba,hrunoff,hrain,hhort,        &
                 ptstep,pd_g,pdzg,                &
                 pwsat,pwfc,pwwilt,               &
                 pwg, pwgi, kwg_layer,            &
                 zpg, zpg_melt, pmuf,             &
                 pcondsat, pbcoef,                &
                 pmpotsat, pksat_ice, pd_ice,     &
                 pfsat, phorton, zdunne, pfflood, &
                 ppiflood, piflood, ppflood,      &
                 prunoffb, prunoffd, pcg,         &
                 psoilwght, oflood, klayer_hort,  &
                 klayer_dun                       )         
!
!-------------------------------------------------------------------------------
!
!*       6.     EVOLUTION OF THE SOIL WATER CONTENT
!               -----------------------------------
!
!*       7.     EFFECT OF MELTING/FREEZING ON SOIL ICE AND LIQUID WATER CONTENTS
!               ----------------------------------------------------------------
!
!*       8.     DRAINAGE FROM THE DEEP SOIL
!               ---------------------------
!
!*      9.     RUN-OFF 
!               -------
!                                     when the soil water exceeds saturation, 
!                                     there is fast-time-response runoff
!
!
! -----------------------------------------------------------------
! Time splitting parameter for *very large time steps* since Richard
! and/or soil freezing equations are very non-linear 
! NOTE for NWP/GCM type applications, the time step is generally not split
! (usually just for offline applications with a time step on order of 
! 15 minutes to an hour for example)
! ------------------------------------------------------------------
!
indt = 1
IF(ptstep>=ztimemax)THEN
  indt = max(2,nint(ptstep/ztimemax))
ENDIF
!
ztstep  = ptstep/REAL(indt)
!
! ------------------------------------------------------------------
! The values for the two coefficients (multiplied by VEG and LAI) 
! in the expression below are from 
! Giard and Bazile (2000), Mon. Wea. Rev.: they model the effect of insolation due to
! vegetation cover. This used by both 'DEF' (code blocks 3.-4.) and 'DIF' options.
! ------------------------------------------------------------------
!
WHERE(plai(:)/=xundef .AND. pveg(:)/=0.)
    zksfc_iveg(:) = (1.0-zinsolfrz_veg*pveg(:)) * (1.0-(plai(:)/zinsolfrz_lai))
ELSEWHERE
    zksfc_iveg(:) = 1.0 ! No vegetation
ENDWHERE
!
!
zwgi0(:,:) = pwgi(:,:) ! save initial ice content before phase changes and sublimation
!
IF (hisba=='DIF') THEN                
!
  ini = SIZE(pd_g(:,:),1)
  inl = maxval(kwg_layer(:))
!
! Initialize some field
! ---------------------
!
  zps(:,:)=xundef
  zti(:,:)=xundef
  DO jl=1,inl
     DO jj=1,ini
        idepth=kwg_layer(jj)
        IF(jl<=idepth)THEN
          zps(jj,jl) = pps(jj)
          zti(jj,jl) = min(xtt,ptg(jj,jl))
        ENDIF
     ENDDO
  ENDDO
!
! Compute specific humidity at saturation for the vapor conductivity
! ------------------------------------------------------------------
!
  zqsat(:,:) = qsat(ptg(:,:),zps(:,:),kwg_layer(:),inl)
  zqsati(:,:) = qsati(zti(:,:),zps(:,:),kwg_layer(:),inl)
!
! Soil water sink terms: convert from (W m-2) and (kg m-2 s-1) to (m s-1)
! ------------------------------------------------------------------
!
  zpg(:) =  zpg(:)        / xrholw
  zevapcor(:) = zevapcor(:)        / xrholw
  zleg(:) =  zleg(:)        /(xrholw*plvtt(:))
  zletr(:) = (zletr(:)/zf2(:))/(xrholw*plvtt(:))
  zlegi(:) = zlegi(:)        /(xrholw*plstt(:))
!
  DO jdt = 1,indt
!                      
     CALL hydro_soildif(hrunoff, hhort, ztstep,                     &
                pbcoef, pwsat, pcondsat, pmpotsat, pwfc,            &
                pd_g, pdzg, pdzdif, zpg, zletr, zleg, zevapcor,     &
                pf2wght, pwg, pwgi, ptg, pps, zqsat, zqsati,        &
                zdrain, zhorton, pfsat, kwg_layer, inl,             &
                klayer_hort, ptopqs, zqsb, pfwtd, pwtd              )
!
     CALL ice_soildif(ztstep, ptauice, zksfc_iveg, zlegi,    &
                     psoilhcapz, pwsat, pmpotsat, pbcoef,    &
                     ptg, pwgi, pwg, kwg_layer,              &
                     pdzg,  zwgi_excess                      )
!
     pdrain(:)  = pdrain(:) + (zdrain(:)+zqsb(:)+zwgi_excess(:))/REAL(indt)
     pqsb(:)  = pqsb(:) + zqsb(:)/REAL(indt)
     phorton(:)  = phorton(:) + zhorton(:)/REAL(indt)
!
! Output diagnostics:
! Compute latent heating from phase change only in surface layer and total soil column,
! then adjust surface and total soil heat content to maintain balance.
!
     pdelphaseg_sfc(:)    = (pwgi(:,1)-zwgi0(:,1))*(xlmtt*xrholw/ptstep)*pdzg(:,1) + zlegi(:)*(xrholw*xlstt)
     pdelphaseg(:)        = pdelphaseg_sfc(:)
     DO jl=2,inl
        DO jj=1,ini
           pdelphaseg(jj) = pdelphaseg(jj) + (pwgi(jj,jl)-zwgi0(jj,jl))*(xlmtt*xrholw/ptstep)*pdzg(jj,jl)
        ENDDO
     ENDDO
     pdelheatg_sfc(:)     = pdelheatg_sfc(:) + pdelphaseg_sfc(:)
     pdelheatg(:)         = pdelheatg(:)     + pdelphaseg(:)

  ENDDO
!
ELSE
!
  DO jdt = 1,indt
!
!   Only layer 1 and 2 are used for soil freezing (ZWG3 not used)
    CALL ice_soilfr(hsnow_isba, hsoilfrz, ztstep, zksfc_iveg, pcg, pct,  &
                    ppsng, pffg, ptauice, zdwgi1, zdwgi2, pwsat,         &
                    pmpotsat, pbcoef, pd_g, ptg, pwgi, pwg               )
!
    CALL hydro_soil(hisba,                                           &
                    ztstep,                                          &
                    zletr, zleg, zpg, zevapcor,                      &
                    pwdrain,                                         &
                    pc1, pc2, pc3, pc4b, pc4ref, pwgeq,              &
                    pd_g(:,2), zdg3, zwsat_avg, zwfc_avg,            &
                    zdwgi1, zdwgi2, zlegi, pd_g(:,1), pcg, pct,      &
                    ptg(:,1), ptg(:,2),                              &
                    pwg(:,1), pwg(:,2), zwg3(:),                     &
                    pwgi(:,1), pwgi(:,2),                            &
                    zrunoff,zdrain,hksat,zwwilt_avg                  )
!
    pdrain(:)  = pdrain(:) + zdrain(:)/REAL(indt)
    prunoff(:)  = prunoff(:) + zrunoff(:)/REAL(indt)
!    
  ENDDO
!
! Output diagnostics:
! Compute latent heating from phase change only in surface layer and total soil column,
! then adjust surface and total soil heat content to maintain balance.
!
  pdelphaseg_sfc(:)    = (pwgi(:,1)-zwgi0(:,1))*(xlmtt*xrholw/ptstep)*pd_g(:,1) + zlegi(:)
  pdelphaseg(:)        = (pwgi(:,2)-zwgi0(:,2))*(xlmtt*xrholw/ptstep)*pd_g(:,2)
  pdelheatg_sfc(:)     = pdelheatg_sfc(:) + pdelphaseg_sfc(:)
  pdelheatg(:)         = pdelheatg(:)     + pdelphaseg(:)
!
  IF (hisba == '3-L') pwg(:,3) = zwg3(:)
!
#ifdef TOPD
  IF (lcoupl_topd) THEN
    !runoff topo cumule (kg/m²)
    DO jj=1,SIZE(nmaskt_patch)
      IF  (nmaskt_patch(jj)/=0) THEN
        IF ( xatop(nmaskt_patch(jj))/=xundef) THEN
          xrunoff_top(nmaskt_patch(jj)) = xrunoff_top(nmaskt_patch(jj)) + &
                                          (prunoff(jj)+ phorton(jj))*xatop(nmaskt_patch(jj))*ptstep
          IF (hrunoff=='TOPD') THEN
            xrunoff_top(nmaskt_patch(jj)) = xrunoff_top(nmaskt_patch(jj)) + zdunne(jj)*ptstep
          ELSE
            ! ZDUNNE contains only saturated pixels on mesh so only catchment
            xrunoff_top(nmaskt_patch(jj)) = xrunoff_top(nmaskt_patch(jj)) + zdunne(jj)*xatop(nmaskt_patch(jj))*ptstep
          ENDIF  
        ENDIF
      ENDIF
      ! ZDUNNE concerns all the mesh so not only catchment =>*XATOP
    ENDDO
  ENDIF
#endif
  !
ENDIF
!
!-------------------------------------------------------------------------------
!
! Add sub-grid surface and subsurface runoff to saturation excess:
!
prunoff(:) = prunoff(:) + zdunne(:) + phorton(:)
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('HYDRO',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE hydro