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  (HISBA,HRUNOFF,HRAIN,HHORT,PTSTEP,           &
                               pd_g,pdzg,pwsat,pwfc,pwwilt,pwg,pwgi,     &
                               kwg_layer,ppg,ppg_melt,pmuf,              &
                               pcondsat,pbcoef,pmpotsat,                 &
                               pksat_ice,pd_ice,pfsat,phorton,pdunne,    &
                               pfflood,ppiflood,piflood,ppflood,         &
                               prunoffb,prunoffd,pcg,psoilwght,          &
                               oflood,klayer_hort,klayer_dun             )  
!
!     #####################################################################
!
!!****  *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_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
!
!
!
 CHARACTER(LEN=*),INTENT(IN)      :: hisba   ! hydrology/soil:
!                                           ! '2-L'  = single column
!                                           ! '3-L'  = root zone/baseflow layer
!                                           ! 'DIF'  = N-layer diffusion: Richard's Eq.
!                                           !         (Boone and Etchevers 2001)
!
 CHARACTER(LEN=*),     INTENT(IN) :: hrunoff ! surface runoff formulation
!                                           ! 'WSAT'
!                                           ! 'DT92'
!                                           ! 'SGH ' Topmodel
!
!
 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
!
LOGICAL, INTENT(IN)              :: oflood ! Flood scheme 
!
REAL, INTENT(IN)                 :: ptstep
!                                   timestep of the integration
!
REAL, DIMENSION(:,:), INTENT(IN) :: pwg,pwgi
!                                   PWG   = layer average liquid volumetric water content (m3 m-3)
!                                   PWGI  = layer average frozen volumetric water content (m3 m-3)
!
INTEGER, DIMENSION(:),INTENT(IN) :: kwg_layer  
!                                   KWG_LAYER = Number of soil moisture layers (DIF option)
!
REAL, DIMENSION(:,:), INTENT(IN) :: pd_g,pdzg,pwsat,pwfc,pwwilt
REAL, DIMENSION(:,:), INTENT(IN) :: pcondsat
!                                   PD_G  = layer depth (m)
!                                   PDZG= layer thickness (m)
!                                   PCONDSAT = surface saturated hydraulic conductivity
!                                   PWSAT = soil porosity (m3 m-3)
!
REAL, DIMENSION(:,:), INTENT(IN) :: pbcoef,pmpotsat
!                                   PMPOTSAT = matric potential at saturation (m) (BC parameters)
!                                   PBCOEF   = slope of the retention curve (-) (BC parameters)
!
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(INOUT):: pfsat
!                                   PFSAT = satured fraction
!
REAL, DIMENSION(:), INTENT(INOUT):: ppg
REAL, DIMENSION(:), INTENT(IN)   :: ppg_melt, pmuf
!                                   PPG      = water reaching the ground
!                                   PPG_MELT = snowmelt reaching the ground
!                                   PMUF      = wet fraction reached by rain
!
REAL, DIMENSION(:), INTENT(IN)   :: pksat_ice, pd_ice
!                                   PKSAT_ICE = hydraulic conductivity at saturation (m s-1)
!                                               on frozen soil depth (Horton calculation)
!                                   PD_ICE    = depth of the soil column for
!                                               fraction of frozen soil calculation (m)
!
REAL, DIMENSION(:), INTENT(OUT)  :: pdunne, phorton
!                                   PDUNNE  = Dunne runoff
!                                   PHORTON = Horton runoff
!
REAL, DIMENSION(:), INTENT(IN   ) :: pfflood
REAL, DIMENSION(:), INTENT(IN   ) :: ppiflood
!                                    PPIFLOOD = Floodplain potential infiltration [kg/m²/s]
!                                             = Floodplain mass
REAL, DIMENSION(:), INTENT(INOUT) :: piflood, ppflood
!                                    PIFLOOD = Floodplain infiltration     [kg/m²/s]
!                                    PPFLOOD = Floodplain interception     [kg/m²/s]
!
REAL, DIMENSION(:), INTENT(IN)    :: prunoffb ! slope of the runoff curve
REAL, DIMENSION(:), INTENT(IN)    :: prunoffd
!                                    PRUNOFFD = depth over which sub-grid runoff calculated (m)
!
REAL, DIMENSION(:), INTENT(IN)    :: pcg
!                                   PCG = soil heat capacity to compute thermal penetration depth
!
!*      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(PD_G,1),SIZE(PD_G,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                                    :: ini, 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)
!
!HRUNOFF = DT92 ZFSAT calculation
zwg2_avg(:)   = 0.0
zwgi2_avg(:)  = 0.0
zwsat_avg(:)  = 0.0
zwwilt_avg(:) = 0.0
!
!HHORT=SGH
zhort_r(:) = 0.0
zhort_m(:) = 0.0
!
!PIFLOOD calculation
zsoilmax(:) = 0.0
zif_max(:) = 0.0
zpifldmax(:) = 0.0
!
!HRUNOFF = DT92 DUNNE calculation
zpg_work(:)   = 0.0
zruisdt(:)    = 0.0
!
!HRUNOFF = TOPD DUNNE calculation
zrunoff_topd(:) = 0.0
!
!to limit numerical artifacts
zpg_ini(:) = ppg(:) + ppg_melt(:)
!
!
ini=SIZE(pd_g,1)
inl=maxval(kwg_layer(:))
!
!-------------------------------------------------------------------------------
!
!*           1. Surface saturated fraction
!            -----------------------------
!
IF( hrunoff=='DT92' .OR. hrunoff == 'TOPD' )THEN
!
! Calculate the layer average water content for the sub-grid
! surface runoff computation: use PRUNOFFD 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 (hisba == 'DIF') THEN
!
! Vertically averaged soil properties and moisture for surface runoff computation:
!
      DO jl=1,klayer_dun
         DO jj=1,ini
            idepth=kwg_layer(jj)
            IF(jl<=idepth)THEN
              zwg2_avg(jj) = zwg2_avg(jj) + psoilwght(jj,jl)*pwg(jj,jl)/max(1.e-6,prunoffd(jj))
              zwgi2_avg(jj) = zwgi2_avg(jj) + psoilwght(jj,jl)*pwgi(jj,jl)/max(1.e-6,prunoffd(jj))
              zwsat_avg(jj) = zwsat_avg(jj) + psoilwght(jj,jl)*pwsat(jj,jl)/max(1.e-6,prunoffd(jj))
              zwwilt_avg(jj) = zwwilt_avg(jj) + psoilwght(jj,jl)*pwwilt(jj,jl)/max(1.e-6,prunoffd(jj))
            ENDIF
         ENDDO
      ENDDO
!
   ELSE
!           
      zwg2_avg(:)   = pwg(:, 2)
      zwgi2_avg(:)  = pwgi(:,2)
      zwsat_avg(:)  = pwsat(:,1)
      zwwilt_avg(:) = pwwilt(:,1)
!      
   ENDIF
!
   IF(hhort=='SGH')THEN
     !runoff over frozen soil explicitly calculated
     zwgi2_avg(:)=0.0
   ENDIF
!
   DO jj=1,ini
      zs=min(1.0,(zwg2_avg(jj)+zwgi2_avg(jj)-zwwilt_avg(jj))/(zwsat_avg(jj)-zwwilt_avg(jj)))
      pfsat(jj) = 1.0-(1.0-max(0.0,zs))**(prunoffb(jj)/(prunoffb(jj)+1.))
   ENDDO        
!
ENDIF
!
!*           2. Horton runoff
!            ----------------
!
IF(hhort=='SGH'.OR.oflood)THEN  
!
  IF(hisba == 'DIF')THEN
!
!   no subgrid frozen soil fraction of the grid cells
    zfrozen(:) = 0.0
!    
    DO jl=1,klayer_hort
      DO jj=1,ini   
!              
!       Modify soil porosity as ice assumed to become part
!       of solid soil matrix (with respect to liquid flow):                
        zwsat(jj,jl) = max(xwgmin, pwsat(jj,jl)-pwgi(jj,jl)) 
        zwfc(jj,jl) = pwfc(jj,jl)*zwsat(jj,jl)/pwsat(jj,jl)
!        
!       Impedance Factor from (Johnsson and Lundin 1991).
        zfrz(jj,jl) = exp(zlog10*(-zeice*(pwgi(jj,jl)/(pwgi(jj,jl)+pwg(jj,jl)))))
!
      ENDDO
    ENDDO    
!
!   Calculate infiltration MAX using green-ampt approximation (derived form)
    zimax(:) = infmax_func(pwg,zwsat,zfrz,pcondsat,pmpotsat,pbcoef,pdzg,pd_g,klayer_hort)
!  
  ELSE
!
    DO jj=1,ini
!
!    Total soil volumetric heat capacity [J/(m3 K)]:
!
      zsoilheatcap = xcl*xrholw*pwg(jj,2) +                           &
                     xci*xrholi*pwgi(jj,2) +                           &
                     xsphsoil*xdrywght*(1.0-pwsat(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(pd_g(jj,2), 4./(zsoilheatcap*pcg(jj)))
!    
!     Effective frozen depth penetration 
!
      zeffice=pd_g(jj,2)*pwgi(jj,2)/(pwgi(jj,2)+pwg(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, pwsat(jj,1)-pwgi(jj,2)) 
!
!     calculate the subgrid frozen soil fraction of the grid cells
!
      zfrozen(jj) = min(1.,zeffice/max(pd_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 to the unsaturated conductivity function at a
!     water content corresponding to the total porosity less the ice-filled volume.
!
      zs =min(1.,zwsat(jj,1)/pwsat(jj,1))
      zimax_ice(jj)=zfrz(jj,1)*pksat_ice(jj)*(zs**(2*pbcoef(jj,1)+3.))
!
!     Calculate infiltration MAX on unfrozen soil using green-ampt approximation
!    
      zs   =min(1.,pwg(jj,2)/zwsat(jj,1))
      zd_h =min(0.10,pd_g(jj,2))
      zimax(jj)=pcondsat(jj,1)*(pbcoef(jj,1)*pmpotsat(jj,1)*(zs-1.0)/zd_h+1.0)
!
    ENDDO
!
  ENDIF
!
ENDIF
!
IF(hhort=='SGH')THEN
!
! calculate the Horton runoff generated by the rainfall rate
!
  IF(hrain=='SGH')THEN
!
    WHERE(ppg(:)>0.)
       zhort_r(:) = (1.- zfrozen(:))* ppg(:)/((zimax(:)*xrholw*pmuf(:)/ppg(:)) + 1.) & !unfrozen soil
                  +      zfrozen(:) * ppg(:)/((zimax_ice(:)*xrholw*pmuf(:)/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(pfflood(:)<=pfsat(:))
        phorton(:) = (1. -   pfsat(:)) * (zhort_r(:) + zhort_m(:))
  ELSEWHERE
        phorton(:) = (1. - pfflood(:)) * (zhort_r(:) + zhort_m(:))
  ENDWHERE
!
ELSE
!
  phorton(:) = 0.0
!
ENDIF
!
! calculate all water reaching the ground
!
ppg(:) = ppg(:) + ppg_melt(:)        
!
!
!*           3. Dunne runoff and flood interception
!            --------------------------------------
!
! Interception by the flooplains
!
IF(oflood)THEN
  ppflood(:)=pfflood(:)*max(0.0,ppg(:))
ELSE
  ppflood(:)=0.0
ENDIF
!
IF(hrunoff=='SGH ')THEN
!        
! calculate the Dunne runoff with TOPMODEL
!
  pdunne(:) = max(ppg(:),0.0) * max(pfsat(:)-pfflood(:),0.0)
!
ELSEIF (hrunoff=='DT92' .OR. hrunoff=='TOPD')THEN
!
!*       Dumenil et Todini (1992)  RUNOFF SCHEME
!        ---------------------------------------         
!
! surface runoff done only on the Fsat-Fflood fraction
!
  zpg_work(:) = ppg(:) - phorton(:) - ppflood(:)
!
#ifdef TOPD
  IF ( lcoupl_topd.AND.hrunoff == '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,                                &
                  prunoffb, zwwilt_avg,                  &
                  prunoffd, zwsat_avg,                   &
                  zwg2_avg, zwgi2_avg,                   &
                  zpg_work, zruisdt                      )
!
  pdunne(:) = zruisdt(:)*prunoffd(:)*xrholw/ptstep
  !
#ifdef TOPD
  IF (lcoupl_topd.AND.hrunoff == 'TOPD') THEN
    pdunne(:) = zrunoff_topd(:) +  pdunne(:)*(1-xatop(nmaskt_patch(:)))
  ENDIF
#endif
  !
  IF(oflood)THEN
    WHERE(pfflood(:)>=pfsat(:).AND.pfflood(:)>0.0)pdunne(:) = 0.0
  ENDIF   
  !
ELSE
! 
! Default case (no subgrid runoff)
!
  pfsat(:) = 0.0
  pdunne(:) = 0.0
!
ENDIF
!
! calculate the infiltration rate after runoff
!
ppg(:) = ppg(:) - pdunne(:) - phorton(:) - ppflood(:)
!
! Supress numerical artifacts:
!
WHERE (zpg_ini(:)<0.0)
       ppg(:)     = zpg_ini(:)
       phorton(:) = 0.0
       pdunne(:) = 0.0
       ppflood(:) = 0.0
ENDWHERE
!
!*           4. infiltration rate from floodplains (à revoir pour DF !!!)
!            -------------------------------------
!
IF(oflood)THEN
!
! calculate the maximum flood infiltration
!
  zpifldmax(:) = min(ppiflood(:),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(hisba == 'DIF')THEN
    zdepth(:)=0.0
    DO jl=1,klayer_hort
       DO jj=1,ini
          IF(zdepth(jj)<xhort_depth)THEN
            zsoilmax(jj) = zsoilmax(jj)+max(0.0,zwfc(jj,jl)-pwg(jj,jl))*pdzg(jj,jl)*xrholw/ptstep
            zdepth(jj) = pd_g(jj,jl)
          ENDIF
       ENDDO
    ENDDO
  ELSE
    DO jj=1,ini
       zwsat(jj,1)  = max(xwgmin, pwsat(jj,1)-pwgi(jj,2)) 
       zsoilmax(jj) = max(0.0,zwsat(jj,1)-pwg(jj,2))*pd_g(jj,2)*xrholw/ptstep
    ENDDO
  ENDIF
!
  zsoilmax(:) = min(zsoilmax(:),zif_max(:))
!
  piflood(:) = max(0.0,(pfflood(:)-pfsat(:))) * min(zpifldmax(:),zsoilmax(:))
!
ELSE
!
  piflood(:)=0.0
!
ENDIF
!
!calculate the infiltration rate
!
ppg(:) = ppg(:) + piflood(:)
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('HYDRO_SGH',1,zhook_handle)
!
END SUBROUTINE hydro_sgh