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 (PMESH_SIZE, IO, S, K, NK, NP, NPE, PRAIN )
!     #######################################################################
!
!!****  *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_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_s_t, isba_k_t, isba_nk_t, isba_np_t, isba_npe_t, &
                        isba_p_t, isba_pe_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
!
!
REAL, DIMENSION(:), INTENT(IN) :: PMESH_SIZE
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_k_t), INTENT(INOUT) :: K
TYPE(isba_nk_t), INTENT(INOUT) :: NK
TYPE(isba_np_t), INTENT(INOUT) :: NP
TYPE(isba_npe_t), INTENT(INOUT) :: NPE

!
REAL, DIMENSION(:), INTENT(IN)   :: PRAIN
!                                   PRAIN   = rain rate (kg/m2/s)
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PRAIN))          :: ZDIST, ZBETA, ZFSAT    ! IO%CRAIN = SGH
!                                        ZDIST  = the cell scale (in km)
!                                        ZBETA  = cell scale dependency parameter
!
REAL, DIMENSION(SIZE(PRAIN))          :: ZD_TOP, ZW_TOP, ZQTOP  ! IO%CRUNOFF = 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  ! IO%CRUNOFF = 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
!
TYPE(isba_p_t), POINTER :: PK
TYPE(isba_k_t), POINTER :: KK
TYPE(isba_pe_t), POINTER :: PEK
!
REAL                                  :: ZW_UP, ZW_DOWN
REAL                                  :: ZF_UP, ZF_DOWN, ZSLOPEF
REAL                                  :: ZQ_UP, ZQ_DOWN, ZSLOPEQ
!
INTEGER                               :: INJ, JJ, JI, JP, JTAB, ICOUNT, &
                                         JL, IMASK
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('ISBA_SGH_UPDATE',0,zhook_handle)
!
inj=SIZE(prain,1)
!
zfsat(:) = 0.0
!
!*   1.0 Spatial distribution of precipitation
!    ---------------------------------------------
!
IF(io%CRAIN=='SGH')THEN
!
  WHERE(prain(:)>0.0)
    k%XMUF (:) =1.0
  ELSEWHERE
    k%XMUF (:) =0.0
  ENDWHERE

!        
! calculate the cell scale (in km)
!
  zdist(:) = sqrt(pmesh_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
!   
    k%XMUF (:) = 1.0 - exp(-zbeta(:)*(prain(:)*xstohr))
!
  ENDWHERE
!
    DO jp = 1,io%NPATCH
      pk => np%AL(jp)
      kk => nk%AL(jp)
      IF (pk%NSIZE_P>0 )THEN
        DO ji = 1,pk%NSIZE_P
          imask = pk%NR_P(ji)
          kk%XMUF(ji) = k%XMUF(imask)
        ENDDO
      ENDIF
    ENDDO

ENDIF
!
!*   2.0 Computation of the saturated fraction given by TOPMODEL 
!    -----------------------------------------------------------
!
IF(io%CRUNOFF=='SGH')THEN
!
! Calculation of the ative TOPMODEL-soil moisture at 't' (m)
! ---------------------------------------------------------------
!
  zd_top(:) = 0.0
  zwsat_avg(:) = 0.0
  zwwilt_avg(:) = 0.0
  zw_top(:) = 0.0
!
  IF(io%CISBA=='DIF')THEN        
!
    DO jp = 1,io%NPATCH
      kk => nk%AL(jp)
      pk => np%AL(jp)
      pek => npe%AL(jp)

      IF (pk%NSIZE_P>0 )THEN

        DO jl = 1,io%NLAYER_DUN
          DO ji = 1,pk%NSIZE_P
            imask = pk%NR_P(ji)
            !
            zd_top(imask) = zd_top(imask) + pk%XPATCH(ji)*pk%XSOILWGHT(ji,jl)
            zwsat_avg(imask) = zwsat_avg(imask) + pk%XPATCH(ji)*pk%XSOILWGHT(ji,jl)*k%XWSAT(imask,jl)
            zwwilt_avg(imask) = zwwilt_avg(imask) + pk%XPATCH(ji)*pk%XSOILWGHT(ji,jl)*k%XWD0 (imask,jl)
            zw_top(imask) = zw_top(imask) + pk%XPATCH(ji)*pk%XSOILWGHT(ji,jl)*pek%XWG (ji,jl)
            !
          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 jp = 1,io%NPATCH
      pk => np%AL(jp)
      pek => npe%AL(jp)

      IF (pk%NSIZE_P>0 )THEN
        DO ji = 1,pk%NSIZE_P
          imask = pk%NR_P(ji)
          !          
          zd_top(imask) = zd_top(imask)+pk%XRUNOFFD(ji)*pk%XPATCH(ji)
          zw_top(imask) = zw_top(imask)+pk%XRUNOFFD(ji)*pk%XPATCH(ji)*pek%XWG(ji,2)
          !
        ENDDO
      ENDIF
    ENDDO
!  
    WHERE(zd_top(:)>0.0)
      zw_top(:) = zw_top(:) / zd_top(:)
    ENDWHERE
!      
    zwsat_avg(:) = k%XWSAT(:,1)
    zwwilt_avg(:) = k%XWD0 (:,1)
!
  ENDIF
!
! Find the boundary
! -----------------
!
  nmask(:)=0
  icount=0
  DO ji=1,inj  
     IF((s%XTI_MEAN(ji)/=xundef.AND.zw_top(ji)<zwsat_avg(ji).AND.zw_top(ji)>zwwilt_avg(ji))) THEN
       icount=icount+1
       nmask(icount)=ji       
     ENDIF
     IF(zw_top(ji)>=zwsat_avg(ji))THEN
        zfsat(ji) = 1.0
     ENDIF
  ENDDO
!     
! compare wt_array and WT
! -----------------------
!
  DO jtab=1,ndimtab
     DO jj=1,icount
        ji = nmask(jj)    
        IF(s%XTAB_WTOP(ji,jtab)>zw_top(ji))THEN
          iup(jj)=jtab
          idown(jj)=jtab+1
        ELSEIF(s%XTAB_WTOP(ji,jtab)==zw_top(ji))THEN
          iup(jj)=jtab
          idown(jj)=jtab
        ENDIF
     ENDDO    
  ENDDO 
!    
! calculate fsat
! --------------
!   
  zqtop(:) = 0.0
!
  DO jj=1,icount
!  
     ji = nmask(jj)
!     
!    new range
     zf_up   = s%XTAB_FSAT(ji,iup(jj))
     zf_down = s%XTAB_FSAT(ji,idown(jj))
     zq_up   = s%XTAB_QTOP(ji,iup(jj))
     zq_down = s%XTAB_QTOP(ji,idown(jj))     
     zw_up   = s%XTAB_WTOP(ji,iup(jj))
     zw_down = s%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
!     
     zfsat(ji) = zf_down+(zw_top(ji)-zw_down)*zslopef
     zqtop(ji) = zq_down+(zw_top(ji)-zw_down)*zslopeq
!     
  ENDDO
!
  DO jp=1,io%NPATCH
    kk => nk%AL(jp)
    pk => np%AL(jp)
    IF(pk%NSIZE_P>0)THEN
      DO ji=1,pk%NSIZE_P
        imask = pk%NR_P(ji)
        kk%XFSAT(ji) = zfsat(imask)
      ENDDO
    ENDIF
  ENDDO
!
! Subsurface flow by layer (m/s)
! ------------------------------
!
  IF(io%CISBA=='DIF')THEN        
!
    DO jp=1,io%NPATCH
      kk => nk%AL(jp)
      pk => np%AL(jp)
      IF(pk%NSIZE_P>0)THEN
        DO jl=1,io%NLAYER_DUN
          DO ji=1,pk%NSIZE_P
            imask = pk%NR_P(ji)
            pk%XTOPQS(ji,jl) = k%XKANISO(imask,jl) * pk%XCONDSAT(ji,1) * zqtop(imask) * &
                                       pk%XSOILWGHT(ji,jl) / pk%XRUNOFFD(ji)
           ENDDO
        ENDDO
      ENDIF
    ENDDO
!
  ENDIF
!
ENDIF
!
IF (lhook) CALL dr_hook('ISBA_SGH_UPDATE',1,zhook_handle)
!
!-------------------------------------------------------------
!
END SUBROUTINE isba_sgh_update