isba_sgh_update.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 isba_sgh_update (IG, I, &
                                  hisba,hrunoff,hrain,prain,pmuf,pfsat,ptopqs)
!     #######################################################################
!
!!****  *SGH_UPDATE*  
!!
!!    PURPOSE
!!    -------
!
!     Calculates the evolution of the fraction, mu, of the grid cell
!     reached by the rain, the Topmodel saturated fraction and the diagnostic
!     wetland fraction.
!         
!     
!!**  METHOD
!!    ------
!
!!    EXTERNAL
!!    --------
!!
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      B. Decharme           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      07/2011 (B. Decharme) : Add fsat diag for dt92
!!      (B. Decharme) 04/2013 : DIF lateral subsurface drainage
!!
!-------------------------------------------------------------------------------
!
!
!*       0.     DECLARATIONS
!               ------------
!
!
!
!
USE modd_isba_grid_n, ONLY : isba_grid_t
USE modd_isba_n, ONLY : isba_t
!
USE modd_sgh_par,     ONLY : ndimtab, xmtokm, xstohr, x001,      &
                             xmuregp, xmurega
!
USE modd_surf_par,   ONLY : xundef
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(isba_grid_t), INTENT(INOUT) :: ig
TYPE(isba_t), INTENT(INOUT) :: i
!
 CHARACTER(LEN=*), INTENT(IN)     :: hisba  ! type of ISBA version:
!                                          ! '2-L' (default)
!                                          ! '3-L'
!                                          ! 'DIF'
!
 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
!
REAL, DIMENSION(:), INTENT(IN)   :: prain
!                                   PRAIN   = rain rate (kg/m2/s)
!
REAL, DIMENSION(:), INTENT(OUT)  :: pmuf
!                                   PMUF = fraction of the grid cell reached by the precipitation
!
REAL, DIMENSION(:), INTENT(OUT)  :: pfsat
!                                   PFSAT   = Topmodel satured fraction
!
REAL, DIMENSION(:,:,:), INTENT(OUT):: ptopqs
!                                     PTOPQS   = Topmodel subsurface flow by layer (m/s)
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PRAIN))          :: zdist, zbeta    ! HRAIN = SGH
!                                        ZDIST  = the cell scale (in km)
!                                        ZBETA  = cell scale dependency parameter
!
REAL, DIMENSION(SIZE(PRAIN))          :: zd_top, zw_top, zqtop  ! HRUNOFF = SGH
!                                        ZW_TOP = ative TOPMODEL-soil moisture at 't' (m3 m-3)
!                                        ZD_TOP = Topmodel active layer
!                                        ZQTOP  = Topmodel lateral sub-surface flow (-)
!
INTEGER, DIMENSION(SIZE(PRAIN))       :: iup,idown  ! HRUNOFF = SGH
!                                        change in xsat (or fsat) index
!
INTEGER, DIMENSION(SIZE(PRAIN))       :: nmask      ! indices correspondance between arrays
!
REAL, DIMENSION(SIZE(PRAIN))          :: zwsat_avg, zwwilt_avg
!                                        Average soil properties content
!
REAL                                  :: zw_up, zw_down
REAL                                  :: zf_up, zf_down, zslopef
REAL                                  :: zq_up, zq_down, zslopeq
!
INTEGER                               :: ini, jj, ji, jpatch, jtab, icount, &
                                         jl
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('ISBA_SGH_UPDATE',0,zhook_handle)
!
ini=SIZE(prain,1)
pfsat(:) = 0.0
!
!*   1.0 Spatial distribution of precipitation
!    ---------------------------------------------
!
IF(hrain=='SGH')THEN
!
  WHERE(prain(:)>0.0)
        pmuf(:) =1.0
  ELSEWHERE
        pmuf(:) =0.0
  ENDWHERE

!        
! calculate the cell scale (in km)
!
  zdist(:) = sqrt(ig%XMESH_SIZE(:))/xmtokm
!
  WHERE(zdist(:)>=15.0)
!
!       calculate beta for the mu calculation
!         
        zbeta(:) = xmurega + xmuregp * exp(-x001*zdist(:))
!
!       calculate mu, precip is in mm/hr
!         
        pmuf(:) = 1.0 - exp(-zbeta(:)*(prain(:)*xstohr))
!
  ENDWHERE
!
ENDIF
!
!*   2.0 Computation of the saturated fraction given by TOPMODEL 
!    -----------------------------------------------------------
!
IF(hrunoff=='SGH')THEN
!
! Calculation of the ative TOPMODEL-soil moisture at 't' (m)
! ---------------------------------------------------------------
!
  zqtop(:) = 0.0
  zw_top(:) = 0.0
  zd_top(:) = 0.0
  zwsat_avg(:) = 0.0
  zwwilt_avg(:) = 0.0
!
  IF(hisba=='DIF')THEN        
!
    DO jpatch=1,i%NPATCH
      IF (i%NSIZE_NATURE_P(jpatch)>0 )THEN
      DO jl=1,i%NLAYER_DUN
        DO jj=1,ini
          zd_top(jj) = zd_top(jj) + i%XPATCH(jj,jpatch)*i%XSOILWGHT(jj,jl,jpatch)
          zwsat_avg(jj) = zwsat_avg(jj) + i%XPATCH(jj,jpatch)*i%XSOILWGHT(jj,jl,jpatch)*i%XWSAT(jj,jl)
          zwwilt_avg(jj) = zwwilt_avg(jj) + i%XPATCH(jj,jpatch)*i%XSOILWGHT(jj,jl,jpatch)*i%XWD0 (jj,jl)
          zw_top(jj) = zw_top(jj) + i%XPATCH(jj,jpatch)*i%XSOILWGHT(jj,jl,jpatch)*i%XWG(jj,jl,jpatch)
        ENDDO
      ENDDO
      ENDIF
    ENDDO
!
    WHERE(zd_top(:)>0.0)
         zwsat_avg(:) = zwsat_avg(:)/zd_top(:)
         zwwilt_avg(:) = zwwilt_avg(:)/zd_top(:)
         zw_top(:) = zw_top(:)/zd_top(:)
    ENDWHERE
!
  ELSE
!    
    DO jpatch=1,i%NPATCH
      IF (i%NSIZE_NATURE_P(jpatch)>0 )THEN
        DO jj=1,ini
          zd_top(jj) = zd_top(jj)+i%XRUNOFFD(jj,jpatch)*i%XPATCH(jj,jpatch)
          zw_top(jj) = zw_top(jj)+i%XRUNOFFD(jj,jpatch)*i%XPATCH(jj,jpatch)*i%XWG(jj,2,jpatch)
        ENDDO
      ENDIF
    ENDDO
!  
    WHERE(zd_top(:)>0.0)
          zw_top(:) = zw_top(:) / zd_top(:)
    ENDWHERE
!      
    zwsat_avg(:) = i%XWSAT(:,1)
    zwwilt_avg(:) = i%XWD0 (:,1)
!
  ENDIF
!
! Find the boundary
! -----------------
!
  nmask(:)=0
  icount=0
  DO jj=1,ini  
     IF((i%XTI_MEAN(jj)/=xundef.AND.zw_top(jj)<zwsat_avg(jj).AND.zw_top(jj)>zwwilt_avg(jj)))THEN     
       icount=icount+1
       nmask(icount)=jj       
     ENDIF
     IF(zw_top(jj)>=zwsat_avg(jj))THEN
        pfsat(jj) = 1.0
     ENDIF
  ENDDO
!     
! compare wt_array and WT
! -----------------------
!
  DO jtab=1,ndimtab
     DO jj=1,icount
        ji = nmask(jj)    
        IF(i%XTAB_WTOP(ji,jtab)>zw_top(ji))THEN
          iup(jj)=jtab
          idown(jj)=jtab+1
        ELSEIF(i%XTAB_WTOP(ji,jtab)==zw_top(ji))THEN
          iup(jj)=jtab
          idown(jj)=jtab
        ENDIF
     ENDDO    
  ENDDO 
!    
! calculate fsat
! --------------
!     
  DO jj=1,icount
!  
     ji = nmask(jj)
!     
!    new range
     zf_up   = i%XTAB_FSAT(ji,iup(jj))
     zf_down = i%XTAB_FSAT(ji,idown(jj))
     zq_up   = i%XTAB_QTOP(ji,iup(jj))
     zq_down = i%XTAB_QTOP(ji,idown(jj))     
     zw_up   = i%XTAB_WTOP(ji,iup(jj))
     zw_down = i%XTAB_WTOP(ji,idown(jj))
!     
!    Calculate new FSAT
     zslopef = 0.0
     zslopeq = 0.0
     IF(iup(jj)/=idown(jj))THEN
       zslopef = (zf_up-zf_down)/(zw_up-zw_down)
       zslopeq = (zq_up-zq_down)/(zw_up-zw_down)
     ENDIF
!     
     pfsat(ji) = zf_down+(zw_top(ji)-zw_down)*zslopef
     zqtop(ji) = zq_down+(zw_top(ji)-zw_down)*zslopeq
!     
  ENDDO
!
! Subsurface flow by layer (m/s)
! ------------------------------
!
  IF(hisba=='DIF')THEN        
!
    DO jpatch=1,i%NPATCH
      IF(i%NSIZE_NATURE_P(jpatch)>0)THEN
        DO jl=1,i%NLAYER_DUN
           DO jj=1,ini
            ptopqs(jj,jl,jpatch)=i%XKANISO(jj,jl)*i%XCONDSAT(jj,1,jpatch)*zqtop(jj)*i%XSOILWGHT(jj,jl,jpatch)/i%XRUNOFFD(jj,jpatch)
           ENDDO
        ENDDO
      ENDIF
    ENDDO
!
  ENDIF
!
ENDIF
!
IF (lhook) CALL dr_hook('ISBA_SGH_UPDATE',1,zhook_handle)
!
!-------------------------------------------------------------
!
END SUBROUTINE isba_sgh_update