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(IO, KK, PK, PEK, AG, DEK, DMK, OMEB, PTSTEP, PVEG, &
                       PWRMAX, PSNOW_THRUFAL, PEVAPCOR, PSUBVCOR, PSOILHCAPZ, &
                       PF2WGHT, PF2, PPS, PIRRIG_GR, PDELHEATG, PDELHEATG_SFC,&
                       PDELPHASEG, PDELPHASEG_SFC               )
!     #####################################################################
!
!!****  *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_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_k_t, isba_p_t, isba_pe_t
USE modd_agri_n, ONLY : agri_t
USE modd_diag_evap_isba_n, ONLY : diag_evap_isba_t
USE modd_diag_misc_isba_n, ONLY : diag_misc_isba_t
!
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
!
!
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
TYPE(agri_t), INTENT(INOUT) :: AG
TYPE(diag_evap_isba_t), INTENT(INOUT) :: DEK
TYPE(diag_misc_isba_t), INTENT(INOUT) :: DMK
!
LOGICAL, INTENT(IN)                :: OMEB   ! True  = patch with multi-energy balance 
!                                            ! False = patch with classical (composite) ISBA
REAL, INTENT(IN)                    :: PTSTEP
!                                      timestep of the integration
!
REAL, DIMENSION(:), INTENT(IN)    :: PVEG, PWRMAX
!                                      PVEG = fraction of vegetation 
!                                      PWRMAX = maximum equivalent water content
!                                               in the vegetation canopy
!
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)  :: PF2WGHT
!                                    PF2WGHT   = water stress factor (profile) (-)
!
REAL, DIMENSION(:,:), INTENT(IN) :: PSOILHCAPZ
!                                   PSOILHCAPZ = ISBA-DF Soil heat capacity profile [J/(m3 K)]
!
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(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_GR ! ground irrigation rate (kg/m2/s)
!
!*      0.2    declarations of local variables
!
!
INTEGER                         :: JJ, JL      ! loop control                                       
INTEGER                         :: INDT, JDT   ! Time splitting indicies
INTEGER                         :: INJ, INL, IDEKTH ! (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(PEK%XWG,1),SIZE(PEK%XWG,2)) :: ZQSAT, ZQSATI, ZTI, ZPS
!                                           For specific humidity at saturation computation (ISBA-DIF)
!
REAL, DIMENSION(SIZE(PEK%XWG,1),SIZE(PEK%XWG,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(:)           = dmk%XRRSFC(:)
!
zdrain(:)        = 0.
zhorton(:)       = 0.
zrunoff(:)       = 0.
zwgi_excess(:)   = 0.
zevapcor(:)      = 0.
zqsb(:)      = 0.
!
dek%XDRAIN(:)        = 0.
dek%XRUNOFF(:)       = 0.
dek%XHORT(:)       = 0.
dek%XQSB   (:)       = 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(:)          = dek%XLEV(:)
   zletr(:)         = dek%XLETR(:)
   zleg(:)          = dek%XLEG(:)
   zlegi(:)         = dek%XLEGI(:)
   zpsnv(:)         = 0.0
!
   zevapcor(:)      = pevapcor(:) + psubvcor(:)
!
ELSE
!
! Initialize evaporation components: variable definitions
! depend on snow scheme:
!
   IF(pek%TSNOW%SCHEME == '3-L' .OR. pek%TSNOW%SCHEME == 'CRO' .OR. io%CISBA == 'DIF')THEN
      zlev(:)          = (1.0-pek%XPSNV(:)-kk%XFFV(:)) * dek%XLEV(:)
      zletr(:)         = (1.0-pek%XPSNV(:)-kk%XFFV(:)) * dek%XLETR(:)
      zleg(:)          = (1.0-pek%XPSNG(:)-kk%XFFG(:)) * dek%XLEG(:)
      zlegi(:)         = (1.0-pek%XPSNG(:)-kk%XFFG(:)) * dek%XLEGI(:)
      zpsnv(:)         = 0.0
   ELSE
      zlev(:)          = dek%XLEV(:)
      zletr(:)         = dek%XLETR(:)
      zleg(:)          = dek%XLEG(:)
      zlegi(:)         = dek%XLEGI(:)
      zpsnv(:)         = pek%XPSNV(:)+kk%XFFV(:)
   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(io%CISBA == '2-L' .OR. io%CISBA == '3-L')THEN  
   zwsat_avg(:)     = kk%XWSAT(:,1)
   zwwilt_avg(:)    = kk%XWWILT(:,1)
   zwfc_avg(:)      = kk%XWFC(:,1)
ENDIF
!
IF (io%CISBA == '3-L') THEN                                   
   zdg3(:) = pk%XDG(:,3)
   zwg3(:) = pek%XWG(:,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.
!
   dek%XIRRIG_FLUX(:)=0.0
!
!* add irrigation over vegetation to liquid precipitation (rr)
!
!
   IF (ASSOCIATED(ag%LIRRIGATE)) THEN
     IF (SIZE(ag%LIRRIGATE)>0) THEN
       WHERE (ag%LIRRIGATE(:) .AND. pek%XIRRIG(:)>0. .AND. pek%XIRRIG(:) /= xundef .AND. (pf2(:)<ag%XTHRESHOLDSPT(:)) )
         dek%XIRRIG_FLUX(:) = pek%XWATSUP(:) / xday           
         zrr(:) = zrr(:) + pek%XWATSUP(:) / xday
         ag%LIRRIDAY(:)    = .true.           
       END WHERE
     ENDIF
   ENDIF
!
!* interception reservoir and dripping computation
!
   CALL hydro_veg(io%CRAIN, ptstep, kk%XMUF, zrr, zlev, zletr, pveg, &
                  zpsnv,  pek%XWR(:), pwrmax, zpg, dek%XDRIP, dek%XRRVEG, pk%XLVTT  ) 
!
ELSE
!
! For MEB case, interception interactions already computed and DMK%XRRSFC 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(:) = dmk%XRRSFC(:)
!
ENDIF
!
!* add irrigation over ground to potential soil infiltration (pg)
!
dek%XIRRIG_FLUX(:) = dek%XIRRIG_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(pek%TSNOW%SCHEME == '3-L' .OR. pek%TSNOW%SCHEME == 'CRO' .OR. io%CISBA == '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.
!
  dek%XMELT(:) = dek%XMELT(:) + psnow_thrufal(:)          ! [kg/(m2 s)]
!
ELSE
  !
  CALL hydro_snow(io%LGLACIER, ptstep, pk%XVEGTYPE_PATCH(:,:), dmk%XSRSFC, &
                  dek%XLES, dek%XMELT, pek%TSNOW, 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(io, kk, pk, pek, dek, dmk, ptstep, zpg, zpg_melt, zdunne )         
!
!----------------------------------------------------------------------------
!
!*       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(pek%XLAI(:)/=xundef .AND. pveg(:)/=0.)
    zksfc_iveg(:) = (1.0-zinsolfrz_veg*pveg(:)) * min(max(1.0-(pek%XLAI(:)/zinsolfrz_lai),0.0),1.0)    
ELSEWHERE
    zksfc_iveg(:) = 1.0 ! No vegetation
ENDWHERE
!
!
zwgi0(:,:) = pek%XWGI(:,:) ! save initial ice content before phase changes and sublimation
!
IF (io%CISBA=='DIF') THEN                
!
  inj = SIZE(pk%XDG(:,:),1)
  inl = maxval(pk%NWG_LAYER(:))
!
! Initialize some field
! ---------------------
!
  zps(:,:)=xundef
  zti(:,:)=xundef
  DO jl=1,inl
     DO jj=1,inj
        idekth=pk%NWG_LAYER(jj)
        IF(jl<=idekth)THEN
          zps(jj,jl) = pps(jj)
          zti(jj,jl) = min(xtt,pek%XTG(jj,jl))
        ENDIF
     ENDDO
  ENDDO
!
! Compute specific humidity at saturation for the vapor conductivity
! ------------------------------------------------------------------
!
  zqsat(:,:) = qsat(pek%XTG,zps,pk%NWG_LAYER,inl)
  zqsati(:,:) = qsati(zti,zps,pk%NWG_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*pk%XLVTT(:))
  zletr(:) = (zletr(:)/zf2(:))/(xrholw*pk%XLVTT(:))
  zlegi(:) = zlegi(:)        /(xrholw*pk%XLSTT(:))
!
  DO jdt = 1,indt
!                      
    CALL hydro_soildif(io, kk, pk, pek, ztstep, zpg, zletr, zleg, zevapcor,  &
                       pf2wght, pps, zqsat, zqsati, zdrain, zhorton, inl, zqsb )
!

    CALL ice_soildif(kk, pk, pek, ztstep, zksfc_iveg, zlegi, psoilhcapz, zwgi_excess )
!
    dek%XDRAIN(:) = dek%XDRAIN(:) + (zdrain(:)+zqsb(:)+zwgi_excess(:))/REAL(indt)
    dek%XQSB  (:) = dek%XQSB  (:) + zqsb(:)/REAL(indt)
    dek%XHORT (:) = dek%XHORT (:) + 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(:) = (pek%XWGI(:,1)-zwgi0(:,1))*(xlmtt*xrholw/ptstep)*pk%XDZG(:,1) + &
                         zlegi(:)*(xrholw*xlstt)
    pdelphaseg(:)     = pdelphaseg_sfc(:)
    DO jl=2,inl
      DO jj=1,inj
        pdelphaseg(jj) = pdelphaseg(jj) + (pek%XWGI(jj,jl)-zwgi0(jj,jl))*&
                           (xlmtt*xrholw/ptstep)*pk%XDZG(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(io, kk, pk, pek, dmk, ztstep, zksfc_iveg, zdwgi1, zdwgi2)
!
    CALL hydro_soil(io, kk, pk, pek, dmk, ztstep, zletr, zleg, zpg, zevapcor, zdg3,    &
                    zwsat_avg, zwfc_avg, zdwgi1, zdwgi2, zlegi, zwg3, zrunoff, zdrain, zwwilt_avg )
!
    dek%XDRAIN (:)  = dek%XDRAIN (:) + zdrain(:)/REAL(indt)
    dek%XRUNOFF(:)  = dek%XRUNOFF(:) + 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(:)    = (pek%XWGI(:,1)-zwgi0(:,1))*(xlmtt*xrholw/ptstep)*pk%XDG(:,1) + zlegi(:)
  pdelphaseg(:)        = (pek%XWGI(:,2)-zwgi0(:,2))*(xlmtt*xrholw/ptstep)*pk%XDG(:,2)
  pdelheatg_sfc(:)     = pdelheatg_sfc(:) + pdelphaseg_sfc(:)
  pdelheatg(:)         = pdelheatg(:)     + pdelphaseg(:)
!
  IF (io%CISBA == '3-L') pek%XWG(:,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)) + &
                                          (dek%XRUNOFF(jj)+ dek%XHORT(jj))*xatop(nmaskt_patch(jj))*ptstep
          IF (io%CRUNOFF=='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:
!
dek%XRUNOFF(:) = dek%XRUNOFF(:) + zdunne(:) + dek%XHORT(:)
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('HYDRO',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE hydro