init_top.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.
!     ######spl
      SUBROUTINE init_top (IO, S, K, NK, NP, KLUOUT, PM    )
!
!     #####################################################################
!
!!****  *INIT_TOP*  
!!
!!    PURPOSE
!!    =======
!
!     Calculates the new array of the Datin-Saulnier TOPMODEL framework fonction for xsat and compute each 
!     satured fraction for each xsat value of the grids cells but also the active TOPMODEL-layer array, the 
!     driest fraction array and the normalized mean deficit array.
!     For calculate new array, we use the incomplete gamma function. (see gamma_inc.f for more detail)
! 
!     Note that over land point where topographic index do not exist, a VIC
!     distribution is used with Bcoef at least equal to 0.1. This value can be
!     change in namelist
!
!!
!!    AUTHOR
!!    ------
!!     B. Decharme
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2008
!       B. decharme 04/2013 : DIF lateral drainage
!                             CTI linear regression done in PGD (HTOPREG deleted)
!
!-------------------------------------------------------------------------------
!
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY: isba_s_t, isba_np_t, isba_nk_t, isba_p_t, isba_k_t
!
USE modd_surf_par,ONLY : xundef
!
USE modd_sgh_par, ONLY : x2, x4, xregp, xrega
!
USE modi_dgam
USE modi_gammas
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
!*       0.     DECLARATIONS
!        ===================
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
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
!
INTEGER, INTENT(IN)                  :: KLUOUT
!
REAL,    DIMENSION(:), INTENT(INOUT) :: PM
!                                       PM = exponential decay factor of the local deficit
!
!*      0.2    declarations of local variables
!
TYPE(isba_k_t), POINTER :: KK
TYPE(isba_p_t), POINTER :: PK
!
REAL, DIMENSION(SIZE(PM)) :: ZD_TOP, ZWSAT_AVG, ZWD0_AVG
!                            ZD_TOP  = Topmodel active layer
!
REAL                  :: ZXI, ZPHI, ZNU, ZTI_MEAN, ZTI_MIN, ZTI_MAX, ZTI_STD, ZTI_SKEW
!                        ZTI_MEAN= ti mean after regression
!                        ZXI     = ti pdf parameter
!                        ZPHI    = ti pdf parameter
!                        ZNU     = ti pdf parameter
!
REAL                  :: ZFTOT, ZYMAX, ZYMIN
REAL                  :: ZGYMAX, ZGYMIN
!                        ZFTOT   = total fraction of a grid cell
!                        ZYMAX   = yi maximum variable
!                        ZYMIN   = yi minimum variable
!                        ZGYMAX  = incomplete gamma function for ymax (GAMSTAR result)
!                        ZGYMIN  = incomplete gamma function for ymin (GAMSTAR result)
!
REAL                  :: ZXSAT_IND, ZYSAT, ZY0, ZDMOY, ZXMOY, ZFMED, ZF0
!                        ZXSAT_IND = Satured index for all index 
!                        ZYSAT     = changing variable of satured index
!                        ZY0       = changing variable of dry index
!                        ZDMOY     = grid cell average deficit (= Dbar/M)
!                        ZXMOY     = ti mean value on wet (1-fsat-f0) fraction
!                        ZF0       = dry fraction
!                        ZFMED     = wet (1-fsat-f0) fraction
!
REAL                  :: ZG, ZGYSAT, ZGY0
!                        ZG     = GAM result 
!                        ZGYSAT = the incomplete gamma function for ysat (GAMSTAR result)
!                        ZGY0   = the incomplete gamma function for y0 (GAMSTAR result)
!
REAL                  :: ZD0, ZPAS, ZCOR
!                        ZD0  = Normalized TOPMODEL maximum deficit D0/M coefficient
!                        ZPAS = pas for calculate the new xsat values array
!
INTEGER               :: IFLG, IFLGST
!                        IFLG   = incomplete gamma function error flag (GAM result)
!                        IFLGST = incomplete gamma function error flag (GAMSTAR result)
!                                 (see gamma_inc.f for more detail)
!
REAL                  :: ZNO, ZAR, ZTOT
!
REAL                  :: ZFUP, ZFDOWN, ZQUP, ZQDOWN, ZSLOPEQ, ZWUP, ZWDOWN, ZSLOPEW
!
INTEGER, DIMENSION (1):: ID
!
INTEGER               :: INI, JI, IND, JSI_MIN, JSI_MAX, IPAS, &
                         JL, INL, JP, IMASK
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
!       1 TOPMODEL SURFACE RUNOFF SCHEME
!       ================================
!
!  1.1     initialisation of local variable
!  ----------------------------------------
!
! Grid cells number
!
IF (lhook) CALL dr_hook('INIT_TOP',0,zhook_handle)
!
ini = SIZE(pm(:))
inl = io%NGROUND_LAYER
ipas = SIZE(s%XTAB_FSAT,2)
!
! GAM result (not use here !)
!
zg  = 0.0
!
zd_top(:) = 0.0
zwsat_avg(:) = 0.0
zwd0_avg(:) = 0.0
!
! soil properties for runoff (m)
!
IF (io%CISBA == 'DIF') THEN                                   
!
  DO jp=1,io%NPATCH

    pk => np%AL(jp)
    kk => nk%AL(jp)

    IF (pk%NSIZE_P == 0 ) cycle

    DO jl=1,inl
      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)
        zwd0_avg(imask) = zwd0_avg(imask) + pk%XPATCH(ji) * pk%XSOILWGHT(ji,jl) * k%XWD0 (imask,jl)
        !
      ENDDO
    ENDDO
  ENDDO
!
  WHERE(zd_top(:)>0.0)
    zwsat_avg(:)=zwsat_avg(:)/zd_top(:)
    zwd0_avg(:)=zwd0_avg(:)/zd_top(:)
  ENDWHERE
!
ELSE
!     
  DO jp=1,io%NPATCH
    !
    pk => np%AL(jp)
    kk => nk%AL(jp)    
    !
     IF (pk%NSIZE_P == 0 ) cycle
     !
     DO ji=1,pk%NSIZE_P
       !
       imask = pk%NR_P(ji)
       !
       zd_top(imask)=zd_top(imask) + pk%XRUNOFFD(ji) * pk%XPATCH(ji)
       !
     ENDDO
     !
  ENDDO
  !
  zwsat_avg(:) = k%XWSAT(:,1)
  zwd0_avg(:) = k%XWD0 (:,1)
  !      
ENDIF
!
!
!  1.2     Algorithme
!  ------------------
!
!grid cells loops 
!
zar = 0.0
ztot= 0.0
zno = 0.0
!
DO ji=1,ini
!   
   IF(s%XTI_MEAN(ji)==xundef)THEN
!
!    *Case where the Topographics index are not defined.
!    --------------------------------------------------------         
     zno=zno+1.0
     s%XTAB_FSAT(ji,:)=0.0     
     s%XTAB_WTOP(ji,:)=xundef
     s%XTAB_QTOP(ji,:)=0.0
!     
     pm(ji) =xundef
!
   ELSE 
!
     ztot = ztot + 1.0
!
!    1.2.0 first initialisation
!    --------------------------
!
     zxi       = 0.0        
     zphi      = 0.0
     znu       = 0.0
!      
!    New version : Regressions directly in the pgd
!    1000 meter DEM to 2m DEM (PAN AND KING 2012)
!             
     zti_mean=s%XTI_MEAN(ji)
     zti_min =s%XTI_MIN (ji)
     zti_max =s%XTI_MAX (ji)
     zti_std =s%XTI_STD (ji)
     zti_skew=s%XTI_SKEW(ji)
!
!    Calculate topographic index pdf parameters 
!
!    Numerical problem especialy over Greenland
     IF(zti_skew<=0.2)THEN
!     
       zti_skew=0.2
!       
       WRITE(kluout,*)'TI_SKEW is too low or negatif (=',zti_skew,'),' 
       WRITE(kluout,*)'then PHI is too big for the grid-cell',ji,'So,GAMMA(PHI) -> +inf.'
       WRITE(kluout,*)'The applied solution is to put TI_SKEW = 0.2'
       IF(zti_std<1.0)THEN
         WRITE(kluout,*)'In addition TI_STD is too low (=',zti_std,'),' 
         WRITE(kluout,*)'The applied solution is to put TI_STD = 1.0'
         zti_std=1.0
       ENDIF               
!
       zar  = zar +1.0
!     
       zxi  = zti_skew*zti_std/x2 
       zphi = (zti_std/zxi)**x2
!
     ELSE
!
       zxi  = zti_skew*zti_std/x2 
       zphi = (zti_std/zxi)**x2
!
     ENDIF
!
     znu  = zti_mean-zphi*zxi
!
!    Exponential decay factor of the local deficit
!
     pm(ji) =(zwsat_avg(ji)-zwd0_avg(ji))*zd_top(ji)/x4 
!
!    1.2.1 Calculate grid cell pdf total density FTOT = F(ymin --> ymax)
!    -------------------------------------------------------------------
!
!    Normalized TOPMODEL maximum deficit D0/M coefficient
!
     zd0 = (zwsat_avg(ji)-zwd0_avg(ji))*zd_top(ji)/pm(ji)
!
!    Initialise
!
     zgymax = 0.0
     zgymin = 0.0
     zftot  = 0.0
     zymin  = 0.0
     zymax  = 0.0
!
!    variable changing yi ---> (ti-nu)/xi
!
     zymin = (zti_min-znu)/zxi
     zymax = (zti_max-znu)/zxi
!  
!    Supress numerical artifact
!
     zcor  = abs(min(0.0,zymin))
!
     zymin = max(0.0,zymin+zcor)
     zymax = zymax+zcor
!
!    Errors flags indicating a number of error condition in G and GSTAR
!    (see gamma_inc.f for more detail)
!
     iflg        =0
     iflgst      =0
!
!    Computation of F(0 --> ymin)
!
     CALL dgam(zphi,zymin,10.,zg,zgymin,iflg,iflgst)
!      
!    if the incomplete gamma function don't work, print why
!
     IF (iflgst/=0)THEN
        WRITE(kluout,*)'GRID-CELL =',ji,'FLGST= ',iflgst,'PHI= ',zphi,'YMIN= ',zymin 
        CALL abor1_sfx('INIT_TOP: (1) PROBLEM WITH DGAM FUNCTION')
     ENDIF      
!
!    Computation of F(0 --> ymax)
!
     CALL dgam(zphi,zymax,10.,zg,zgymax,iflg,iflgst)
!      
!    if the incomplete gamma function don't work, print why
!
     IF (iflgst/=0)THEN
        WRITE(kluout,*)'GRID-CELL =',ji,'FLGST= ',iflgst,'PHI= ',zphi,'YMAX= ',zymax
        CALL abor1_sfx('INIT_TOP: (2) PROBLEM WITH DGAM FUNCTION')
     ENDIF      
!
!    FTOT = F(0 --> ymax) - F(0 --> ymin)
!
     zftot=zgymax-zgymin
!
!    initialization water content and fraction
!
     s%XTAB_WTOP(ji,1) = zwsat_avg(ji)
     s%XTAB_FSAT(ji,1) = 1.0
     s%XTAB_QTOP(ji,1) = 0.0
!     
     s%XTAB_WTOP(ji,ipas) = zwd0_avg(ji)
     s%XTAB_FSAT(ji,ipas) = 0.0
     s%XTAB_QTOP(ji,ipas) = 0.0
!
!    Define the new limits for the satured index loop
!
     jsi_min = 2
     jsi_max = ipas-1
     zpas    = (zti_max-zti_min)/REAL(ipas-1)
!
!    1.2.2 Calculate all topmodel arrays
!    -----------------------------------
!
!    Satured index loop
!
     DO ind=jsi_min,jsi_max
!
!       initialize of loops variables
!
        zxsat_ind = 0.0 
        zysat     = 0.0
        zy0       = 0.0
        zdmoy     = 0.0
        zxmoy     = 0.0
        zfmed     = 0.0
!
!       Initialize of incomplete gamma function flags and variables
!
        iflg   = 0
        iflgst = 0
        zgysat = 0.0
        zgy0   = 0.0
!
!       calculate xsat for all new index
!
        zxsat_ind=zti_min+REAL(ind-1)*ZPAS
!
!       Changing variable to compute incomplete gamma function 
!
        zysat=(zxsat_ind-znu)/zxi 
        zy0  =((zxsat_ind-zd0)-znu)/zxi
!      
!       Calculate Y0 and ysat and assume ymin < y0 < ymax !

        zysat=max(zymin,min(zysat+zcor,zymax))
        zy0  =max(zymin,min(zy0  +zcor,zymax))
!
!       call incomplete gamma function for xsat
!
        CALL dgam(zphi,zysat,10.,zg,zgysat,iflg,iflgst)  
!
!       if the incomplete gamma function don't works, print why
!
        IF (iflgst/=0)THEN
           WRITE(kluout,*)'GRID-CELL= ',ji,'FLGST= ',iflgst,'PHI= ',zphi,'YSAT= ',zysat
           CALL abor1_sfx('INIT_TOP: (3) PROBLEM WITH DGAM FUNCTION')
        ENDIF 
!
!       call incomplete gamma function for xsat
!
        CALL dgam(zphi,zy0,10.,zg,zgy0,iflg,iflgst)
!
!       if the incomplete gamma function don't works, print why
!
        IF (iflgst/=0)THEN
           WRITE(kluout,*)'GRID-CELL= ',ji,'FLGST= ',iflgst,'PHI= ',zphi,'Y0= ',zy0
           CALL abor1_sfx('INIT_TOP: (4) PROBLEM WITH DGAM FUNCTION')
        ENDIF 
!
!       compute satured fraction as FSAT = F(0 --> ymax) - F(0 --> ysat)
!       
        s%XTAB_FSAT(ji,ind)=max(0.0,(zgymax-zgysat)/zftot)
!
!       Compute driest fraction
!
        zf0=max(0.0,(zgy0-zgymin)/zftot)
!
!       Calculate FMED
!        
        zfmed=(1.0-s%XTAB_FSAT(ji,ind)-zf0)
!
        IF (zfmed/=0.0) THEN
!
!          Compute the new x mean, xmoy', over the wet fraction Fwet
!
           zxmoy = znu+zxi*(zphi-zcor+(exp(-zy0)*(zy0**(zphi/x2))*(zy0**(zphi/x2))     &
                 - exp(-zysat)*(zysat**(zphi/x2))*(zysat**(zphi/x2)))/(zfmed*gammas(zphi)))
!
!          supress numerical artifacs
!
           zxmoy =max((zxsat_ind-zd0),min(zxsat_ind,zxmoy))
!
!          Calculate the mean normalysed deficit as Dbar/M = (1-fsat-f0)*(xsat-xmoy')+f0*D0/M
!
           zdmoy = zfmed*(zxsat_ind-zxmoy)+zf0*zd0
!
        ENDIF
!
!       supress numerical artifacs
!
        zdmoy = max(0.0,min(zdmoy,zd0))
!
!       Solves Dbar = (Wsat-WT)*d_top with Dbar/M (=ZDMOY) = (Wsat-WT)*d_top/M
!
        s%XTAB_WTOP(ji,ind) = zwsat_avg(ji)-(pm(ji)*zdmoy/zd_top(ji))

!       Solves Qs = FMED * M * Ks * exp(-Xsat) / Ks (dimentionless)
!
        s%XTAB_QTOP(ji,ind) = zfmed*pm(ji)*exp(-zxsat_ind)
!        
      ENDDO
!
    ENDIF
!
ENDDO
!
! supress numerical artifacs for boundaries conditions
!
DO ji=1,ini
!
!  Upper boundary
!
   IF(s%XTAB_WTOP(ji,2)==zwsat_avg(ji))THEN
!
     zfup=s%XTAB_FSAT(ji,1)
     zwup=s%XTAB_WTOP(ji,1)
     zqup=s%XTAB_QTOP(ji,1)
!   
     id(:)=maxloc(s%XTAB_WTOP(ji,:),s%XTAB_WTOP(ji,:)<zwsat_avg(ji))
!   
     zfdown=s%XTAB_FSAT(ji,id(1))
     zwdown=s%XTAB_WTOP(ji,id(1))
     zqdown=s%XTAB_QTOP(ji,id(1))
!     
     zslopew=(zwup-zwdown)/(zfup-zfdown)   
     zslopeq=(zqup-zqdown)/(zfup-zfdown)   
!
     DO ind=2,id(1)-1
        s%XTAB_WTOP(ji,ind)=zwdown+(s%XTAB_FSAT(ji,ind)-zfdown)*zslopew
        s%XTAB_QTOP(ji,ind)=zqdown+(s%XTAB_FSAT(ji,ind)-zfdown)*zslopeq
     ENDDO
!   
   ENDIF
!
!  Lower boundary
!
   WHERE(s%XTAB_FSAT(ji,:)<=0.0      )
         s%XTAB_WTOP(ji,:)=zwd0_avg(ji)
         s%XTAB_QTOP(ji,:)=0.0
   ENDWHERE
   WHERE(s%XTAB_WTOP(ji,:)<=zwd0_avg(ji))
         s%XTAB_FSAT(ji,:)=0.0
         s%XTAB_QTOP(ji,:)=0.0
   ENDWHERE
!   
ENDDO
!
WRITE(kluout,*)'-------------------TOPMODEL SUM-UP-------------------------'
WRITE(kluout,*)'Number of grid-cells ',ini,'Number of Topmodel points',int(ztot)
IF(ini/=0) THEN
  WRITE(kluout,*)'Percentage of non TOPMODEL grid-cells',(100.*zno/float(ini))
ENDIF
IF(ztot>0.0)THEN
  WRITE(kluout,*)'Percentage of arranged (TI-SKE=0.2) grid-cells',(100.*zar/ztot)
ENDIF
WRITE(kluout,*)'-----------------------------------------------------------'
!
IF (lhook) CALL dr_hook('INIT_TOP',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE init_top