hydro_sgh.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_sgh(IO, KK, PK, PEK, DEK, DMK, PTSTEP, PPG, PPG_MELT, PDUNNE )  
!
!     #####################################################################
!
!!****  *HYDRO_SGH*  
!!
!!    PURPOSE
!!    =======
!
!     1. Determine the Horton runoff that take account of a spatial subgri
!        exponential distribution of the precipitation and of the surface ksat.
!     1. Determine the surface saturated fraction (dt92 or Topmodel).
!     3. Determine the Dunne runoff (dt92 or Topmodel).
!     4. Determine the infiltration rate.
!     5. Determine the flooplains interception and infiltration rate.
!
!!    MODIFICATIONS
!!    -------------
!!
!!         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_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, xcl, xci, xrholi
USE modd_isba_par,  ONLY : xwgmin, xsphsoil, xdrywght
USE modd_surf_par,  ONLY : xundef
USE modd_sgh_par,   ONLY : xhort_depth
!
#ifdef TOPD
USE modd_coupling_topd, ONLY : lcoupl_topd, xas_nature, xatop, nmaskt_patch
#endif
!
USE modi_hydro_dt92
!
USE mode_hydro_dif
!
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(diag_evap_isba_t), INTENT(INOUT) :: DEK
TYPE(diag_misc_isba_t), INTENT(INOUT) :: DMK
!
REAL, INTENT(IN)                 :: PTSTEP
!                                   timestep of the integration
!
REAL, DIMENSION(:), INTENT(INOUT):: PPG
REAL, DIMENSION(:), INTENT(IN)   :: PPG_MELT
!                                   PPG      = water reaching the ground
!                                   PPG_MELT = snowmelt reaching the ground
!
REAL, DIMENSION(:), INTENT(OUT)  :: PDUNNE
!                                   PDUNNE  = Dunne runoff
!
!*      0.2    declarations of local variables
!
REAL, PARAMETER                            :: ZEICE = 6.0  ! Ice vertical diffusion impedence factor 
!
REAL, DIMENSION(SIZE(PPG))                 :: ZPG_INI, ZFROZEN, ZIMAX_ICE, ZIMAX, &
                                              ZHORT_R, ZHORT_M, ZSOILMAX, ZIF_MAX,&
                                              ZPIFLDMAX
!                                             ZFROZEN  = frozen soil fraction for runoff
!                                             ZIMAX_ICE    = maximum infiltration rate for frozen soil
!                                             ZIMAX     = maximum infiltration rate for unfrozen soil
!                                             ZPIFLDMAX    = maximum floodplains infiltration during 1 day (kg/m2/s)
REAL, DIMENSION(SIZE(PPG))                 :: ZWG2_AVG, ZWGI2_AVG, ZWSAT_AVG, ZWWILT_AVG
!                                             Average water and ice content
!                                             values over the soil depth D2 (for calculating surface runoff)
!
REAL, DIMENSION(SIZE(PK%XDG,1),SIZE(PK%XDG,2)) :: ZWSAT, ZWFC, ZFRZ
!
REAL, DIMENSION(SIZE(PPG))                 :: ZPG_WORK, ZRUISDT, ZNL_HORT, ZDEPTH
!
REAL, DIMENSION(SIZE(PPG))                 :: ZRUNOFF_TOPD
!
REAL                                       :: ZEFFICE, ZLOG10, ZLOG, ZS, ZD_H
!
REAL                                       :: ZTDIURN, ZSOILHEATCAP
!                                             ZTDIURN      = thermal penetration depth for restore (m)
!                                             ZSOILHEATCAP = Total soil volumetric heat capacity [J/(m3 K)]
!
INTEGER                                    :: INJ, INL, JJ, JL, IDEPTH
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
!allocate
!
IF (lhook) CALL dr_hook('HYDRO_SGH',0,zhook_handle)
!
!initialize
!
zfrozen(:)  = 0.0
zimax_ice(:)  = 0.0
zimax(:)  = 0.0
!
zwsat(:,:)  = 0.0
zwfc(:,:)  = 0.0
!
zlog10 = log(10.0)
!
!IO%CRUNOFF = DT92 ZFSAT calculation
zwg2_avg(:)   = 0.0
zwgi2_avg(:)  = 0.0
zwsat_avg(:)  = 0.0
zwwilt_avg(:) = 0.0
!
!IO%CHORT=SGH
zhort_r(:) = 0.0
zhort_m(:) = 0.0
!
!DEK%XIFLOOD calculation
zsoilmax(:)  = 0.0
zif_max(:)   = 0.0
zpifldmax(:) = 0.0
!
!IO%CRUNOFF = DT92 DUNNE calculation
zpg_work(:)   = 0.0
zruisdt(:)    = 0.0
!
!IO%CRUNOFF = TOPD DUNNE calculation
zrunoff_topd(:) = 0.0
!
!to limit numerical artifacts
zpg_ini(:) = ppg(:) + ppg_melt(:)
!
!
inj=SIZE(pk%XDG,1)
inl=maxval(pk%NWG_LAYER)
!
!-------------------------------------------------------------------------------
!
!*           1. Surface saturated fraction
!            -----------------------------
!
IF( io%CRUNOFF=='DT92' .OR. io%CRUNOFF == 'TOPD' )THEN
!
! Calculate the layer average water content for the sub-grid
! surface runoff computation: use PK%XRUNOFFD(:,1) as the depth over which
! runoff is calculated.
!
! First, determine a weight for each layer's contribution
! to thickness averaged water content and soil properties for runoff.
!
  IF (io%CISBA == 'DIF') THEN
!
! Vertically averaged soil properties and moisture for surface runoff computation:
!
    DO jl=1,io%NLAYER_DUN
      DO jj=1,inj
        idepth=pk%NWG_LAYER(jj)
        IF(jl<=idepth)THEN
          zwg2_avg(jj) = zwg2_avg(jj) + pk%XSOILWGHT(jj,jl)*pek%XWG (jj,jl) /max(1.e-6,pk%XRUNOFFD(jj))
          zwgi2_avg(jj) = zwgi2_avg(jj) + pk%XSOILWGHT(jj,jl)*pek%XWGI(jj,jl) /max(1.e-6,pk%XRUNOFFD(jj))
          zwsat_avg(jj) = zwsat_avg(jj) + pk%XSOILWGHT(jj,jl)*kk%XWSAT (jj,jl)/max(1.e-6,pk%XRUNOFFD(jj))
          zwwilt_avg(jj) = zwwilt_avg(jj) + pk%XSOILWGHT(jj,jl)*kk%XWWILT(jj,jl)/max(1.e-6,pk%XRUNOFFD(jj))
        ENDIF
      ENDDO
    ENDDO
!
  ELSE
!           
    zwg2_avg(:)   = pek%XWG (:,2)
    zwgi2_avg(:)  = pek%XWGI(:,2)
    zwsat_avg(:)  = kk%XWSAT (:,1)
    zwwilt_avg(:) = kk%XWWILT(:,1)
!      
  ENDIF
!
  IF(io%CHORT=='SGH')THEN
    !runoff over frozen soil explicitly calculated
    zwgi2_avg(:)=0.0
  ENDIF
!
   DO jj=1,inj
     zs = min(1.0,(zwg2_avg(jj)+zwgi2_avg(jj)-zwwilt_avg(jj))/(zwsat_avg(jj)-zwwilt_avg(jj)))
     kk%XFSAT(jj) = 1.0-(1.0-max(0.0,zs))**(kk%XRUNOFFB(jj)/(kk%XRUNOFFB(jj)+1.))
   ENDDO        
!
ENDIF
!
!*           2. Horton runoff
!            ----------------
!
IF(io%CHORT=='SGH'.OR.io%LFLOOD)THEN  
!
  IF(io%CISBA == 'DIF')THEN
!
!   no subgrid frozen soil fraction of the grid cells
    zfrozen(:) = 0.0
!    
    DO jl=1,io%NLAYER_HORT
      DO jj=1,inj   
!              
!       Modify soil porosity as ice assumed to become part
!       of solid soil matrix (with respect to liquid flow):                
        zwsat(jj,jl) = max(xwgmin, kk%XWSAT(jj,jl)-pek%XWGI(jj,jl))
        zwfc(jj,jl) = kk%XWFC(jj,jl)*zwsat(jj,jl)/kk%XWSAT(jj,jl) 
!        
!       Impedance Factor from (Johnsson and Lundin 1991).
        zfrz(jj,jl) = exp(zlog10*(-zeice*(pek%XWGI(jj,jl)/(pek%XWGI(jj,jl)+pek%XWG(jj,jl)))))
!
      ENDDO
    ENDDO    
!
!   Calculate infiltration MAX using green-ampt approximation (derived form)
    zimax(:) = infmax_func(pek%XWG, zwsat, zfrz, pk%XCONDSAT, kk%XMPOTSAT, kk%XBCOEF, &
                           pk%XDZG, pk%XDG, io%NLAYER_HORT)
!  
  ELSE
!
    DO jj=1,inj
!
!    Total soil volumetric heat capacity [J/(m3 K)]:
!
      zsoilheatcap = xcl*xrholw*pek%XWG (jj,2) +                           &
                     xci*xrholi*pek%XWGI(jj,2) +                           &
                     xsphsoil*xdrywght*(1.0-kk%XWSAT(jj,1))
!                     
!     Soil thickness which corresponds to the diurnal surface temperature
!     wave penetration depth as T2 is the average temperature for this layer:
!
      ztdiurn   = min(pk%XDG(jj,2), 4./(zsoilheatcap*dmk%XCG(jj)))
!    
!     Effective frozen depth penetration 
!
      zeffice = pk%XDG(jj,2)*pek%XWGI(jj,2)/(pek%XWGI(jj,2)+pek%XWG(jj,2))
!
!     Modify soil porosity as ice assumed to become part
!     of solid soil matrix (with respect to liquid flow):
!
      zwsat(jj,1) = max(xwgmin, kk%XWSAT(jj,1)-pek%XWGI(jj,2)) 
!
!     calculate the subgrid frozen soil fraction of the grid cells
!
      zfrozen(jj) = min(1.,zeffice/max(pk%XD_ICE(jj),ztdiurn))
!
!     Impedance Factor from (Johnsson and Lundin 1991).
!
      zfrz(jj,1) = exp(zlog10*(-zeice*min(1.,zeffice/ztdiurn)))
!
!     Calculate infiltration MAX on frozen soil as Johnsson and Lundin (1991).
!     The max infiltration is equal ,1to the unsaturated conductivity function at a
!     water content c,1orresponding to the total porosity less the ice-filled volume.
!
      zs =min(1.,zwsat(jj,1)/kk%XWSAT(jj,1))
      zimax_ice(jj)=zfrz(jj,1)*pk%XKSAT_ICE(jj)*(zs**(2*kk%XBCOEF(jj,1)+3.))
!
!     Calculate infiltration MAX on unfrozen soil using green-ampt approximation
!    
      zs   =min(1.,pek%XWG(jj,2)/zwsat(jj,1))
      zd_h =min(0.10,pk%XDG(jj,2))
      zimax(jj)=pk%XCONDSAT(jj,1)*(kk%XBCOEF(jj,1)*kk%XMPOTSAT(jj,1)*(zs-1.0)/zd_h+1.0)
!
    ENDDO
!
  ENDIF
!
ENDIF
!
IF(io%CHORT=='SGH')THEN
!
! calculate the Horton runoff generated by the rainfall rate
!
  IF(io%CRAIN=='SGH')THEN
!
    WHERE(ppg(:)>0.)
       zhort_r(:) = (1.- zfrozen(:))* ppg(:)/((zimax(:)*xrholw*kk%XMUF(:)/ppg(:)) + 1.) & !unfrozen soil
                  +      zfrozen(:) * ppg(:)/((zimax_ice(:)*xrholw*kk%XMUF(:)/ppg(:)) + 1.)   !frozen soil
    END WHERE        
!
  ELSE
!
    zhort_r(:) = (1.- zfrozen(:))* max(0.,ppg(:)-zimax(:)*xrholw)          & !unfrozen soil
               +      zfrozen(:) * max(0.,ppg(:)-zimax_ice(:)*xrholw)            !frozen soil
!
  ENDIF
!
! calculate the Horton runoff generated by the snow melt
!        
  zhort_m(:) = (1.- zfrozen(:))* max(0.,ppg_melt(:)-zimax(:)*xrholw)          & !unfrozen soil
             +      zfrozen(:) * max(0.,ppg_melt(:)-zimax_ice(:)*xrholw)            !frozen soil
!
! calculate the  total Horton runoff 
!
  WHERE(kk%XFFLOOD(:)<=kk%XFSAT(:))
        dek%XHORT(:) = (1. - kk%XFSAT(:))   * (zhort_r(:) + zhort_m(:))
  ELSEWHERE
        dek%XHORT(:) = (1. - kk%XFFLOOD(:)) * (zhort_r(:) + zhort_m(:))
  ENDWHERE
!
ELSE
!
  dek%XHORT(:) = 0.0
!
ENDIF
!
! calculate all water reaching the ground
!
ppg(:) = ppg(:) + ppg_melt(:)        
!
!
!*           3. Dunne runoff and flood interception
!            --------------------------------------
!
! Interception by the flooplains
!
IF(io%LFLOOD)THEN
  dek%XPFLOOD(:)=kk%XFFLOOD(:)*max(0.0,ppg(:))
ELSE
  dek%XPFLOOD(:)=0.0
ENDIF
!
IF(io%CRUNOFF=='SGH ')THEN
!        
! calculate the Dunne runoff with TOPMODEL
!
  pdunne(:) = max(ppg(:),0.0) * max(kk%XFSAT(:)-kk%XFFLOOD(:),0.0)
!
ELSEIF (io%CRUNOFF=='DT92' .OR. io%CRUNOFF=='TOPD')THEN
!
!*       Dumenil et Todini (1992)  RUNOFF SCHEME
!        ---------------------------------------         
!
! surface runoff done only on the Fsat-Fflood fraction
!
  zpg_work(:) = ppg(:) - dek%XHORT(:) - dek%XPFLOOD(:)
!
#ifdef TOPD
  IF ( lcoupl_topd.AND.io%CRUNOFF == 'TOPD' )THEN
    !
    DO jj=1,SIZE(nmaskt_patch)
      IF (nmaskt_patch(jj)/=0) THEN
        IF ( xatop(nmaskt_patch(jj))/=0. .AND. xas_nature(nmaskt_patch(jj))/=xundef ) THEN
          zrunoff_topd(jj) = max(ppg(jj),0.0) * max(xas_nature(nmaskt_patch(jj)),0.0)
        ENDIF
      ENDIF 
    ENDDO
    !
  ENDIF
#endif
  !
  CALL hydro_dt92(ptstep, kk%XRUNOFFB, zwwilt_avg,  pk%XRUNOFFD, zwsat_avg, &
                  zwg2_avg, zwgi2_avg, zpg_work, zruisdt           )
!
  pdunne(:) = zruisdt(:)*pk%XRUNOFFD(:)*xrholw/ptstep
  !
#ifdef TOPD
  IF (lcoupl_topd.AND.io%CRUNOFF == 'TOPD') THEN
    pdunne(:) = zrunoff_topd(:) +  pdunne(:)*(1-xatop(nmaskt_patch(:)))
  ENDIF
#endif
  !
  IF(io%LFLOOD)THEN
    WHERE(kk%XFFLOOD(:)>=kk%XFSAT(:).AND.kk%XFFLOOD(:)>0.0)pdunne(:) = 0.0
  ENDIF   
  !
ELSE
! 
! Default case (no subgrid runoff)
!
  kk%XFSAT (:) = 0.0
  pdunne(:) = 0.0
!
ENDIF
!
! calculate the infiltration rate after runoff
!
ppg(:) = ppg(:) - pdunne(:) - dek%XHORT(:) - dek%XPFLOOD(:)
!
! Supress numerical artifacts:
!
WHERE (zpg_ini(:)<0.0)
  ppg(:)     = zpg_ini(:)
  dek%XHORT(:) = 0.0
  pdunne(:) = 0.0
  dek%XPFLOOD(:) = 0.0
ENDWHERE
!
!*           4. infiltration rate from floodplains (à revoir pour DF !!!)
!            -------------------------------------
!
IF(io%LFLOOD)THEN
!
! calculate the maximum flood infiltration
!
  zpifldmax(:) = min(kk%XPIFLOOD(:),xrholw/xday) ! no more than 1 meter of water per days
!
  zif_max(:) = max(0.,(1.- zfrozen(:))) * zimax(:)*xrholw &   !unfrozen soil
             +             zfrozen(:)   * zimax_ice(:)*xrholw     !frozen soil
!
  IF(io%CISBA == 'DIF')THEN
    zdepth(:)=0.0
    DO jl=1,io%NLAYER_HORT
      DO jj=1,inj
        IF(zdepth(jj)<xhort_depth)THEN
          zsoilmax(jj) = zsoilmax(jj)+max(0.0,zwfc(jj,jl)-pek%XWG(jj,jl))*pk%XDZG(jj,jl)*xrholw/ptstep
          zdepth(jj) = pk%XDG(jj,jl)
        ENDIF
      ENDDO
    ENDDO
  ELSE
    DO jj=1,inj
      zwsat(jj,1)  = max(xwgmin, kk%XWSAT(jj,1)-pek%XWGI(jj,2)) 
      zsoilmax(jj) = max(0.0,zwsat(jj,1)-pek%XWG(jj,2))*pk%XDG(jj,2)*xrholw/ptstep
    ENDDO
  ENDIF
!
  zsoilmax(:) = min(zsoilmax(:),zif_max(:))
!
  dek%XIFLOOD(:) = max(0.0,(kk%XFFLOOD(:)-kk%XFSAT(:))) * min(zpifldmax(:),zsoilmax(:))
!
ELSE
!
  dek%XIFLOOD(:)=0.0
!
ENDIF
!
!calculate the infiltration rate
!
ppg(:) = ppg(:) + dek%XIFLOOD(:)
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('HYDRO_SGH',1,zhook_handle)
!
END SUBROUTINE hydro_sgh