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 (I, &
                            hisba, kluout, ppatch, prunoffd,          &
                           pwd0, pwsat, pti_min,                     &
                           pti_max, pti_mean, pti_std, pti_skew,     &
                           psoilwght, ptab_fsat, ptab_wtop,          &
                           ptab_qtop, 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_n, ONLY : isba_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_t), INTENT(INOUT) :: i
!
 CHARACTER(LEN=*), INTENT(IN)         :: hisba   ! type of ISBA version:
!                                               ! '2-L' (default)
!                                               ! '3-L'
!                                               ! 'DIF'
!
INTEGER, INTENT(IN)                  :: kluout
!
REAL,    DIMENSION(:,:), INTENT(IN)  :: ppatch
!
REAL, DIMENSION(:,:,:),INTENT(IN)    :: psoilwght  ! ISBA-DIF: weights for vertical
!                                                  ! integration of soil water and properties
!
REAL,    DIMENSION(:,:), INTENT(IN)  :: prunoffd ! depth over which sub-grid runoff is computed
!
REAL,    DIMENSION(:,:), INTENT(IN)  :: pwd0, pwsat
!                                       PWD0   = water content equivalent to 
!                                                D0 maximum deficit (m3 m-3)
!                                       PWSAT  = saturation volumetric water content
!                                                of the soil (m3 m-3)
!
REAL,    DIMENSION(:), INTENT(IN)    :: pti_min, pti_max, pti_std, &
                                        pti_skew, pti_mean
!                                       PTI_MEAN  = ti mean
!                                       PTI_MIN  = ti min
!                                       PTI_MAX  = ti max
!                                       PTI_STD  = ti standard deviation
!                                       PTI_SKEW = ti skewness
!
REAL, DIMENSION(:,:), INTENT(INOUT)  :: ptab_fsat, ptab_wtop, ptab_qtop
!                                       PTAB_FSAT = Satured fraction array
!                                       PTAB_WTOP = Active TOPMODEL-layer array
!                                       PTAB_QTOP = Subsurface flow TOPMODEL array
!
REAL,    DIMENSION(:), INTENT(INOUT) :: pm
!                                       PM = exponential decay factor of the local deficit
!
!*      0.2    declarations of local variables
!
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, jpatch
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 = SIZE(pwsat,2)
ipas = SIZE(ptab_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 (hisba == 'DIF') THEN                                   
!
  DO jpatch=1,i%NPATCH
    IF (i%NSIZE_NATURE_P(jpatch) == 0 ) cycle
    DO jl=1,inl
       DO ji=1,ini
          zd_top(ji) = zd_top(ji) + ppatch(ji,jpatch)*psoilwght(ji,jl,jpatch)
          zwsat_avg(ji) = zwsat_avg(ji) + ppatch(ji,jpatch)*psoilwght(ji,jl,jpatch)*pwsat(ji,jl)
          zwd0_avg(ji) = zwd0_avg(ji) + ppatch(ji,jpatch)*psoilwght(ji,jl,jpatch)*pwd0(ji,jl)
       ENDDO
    ENDDO
  ENDDO
!
  WHERE(zd_top(:)>0.0)
        zwsat_avg(:)=zwsat_avg(:)/zd_top(:)
        zwd0_avg(:)=zwd0_avg(:)/zd_top(:)
  ENDWHERE
!
ELSE
!     
  DO jpatch=1,i%NPATCH
     IF (i%NSIZE_NATURE_P(jpatch) == 0 ) cycle
     DO ji=1,ini
        zd_top(ji)=zd_top(ji)+prunoffd(ji,jpatch)*ppatch(ji,jpatch)
     ENDDO
  ENDDO
!     
  zwsat_avg(:) = pwsat(:,1)
  zwd0_avg(:) = pwd0(:,1)
!      
ENDIF
!
!
!  1.2     Algorithme
!  ------------------
!
!grid cells loops 
!
zar = 0.0
ztot= 0.0
zno = 0.0
!
DO ji=1,ini
!   
   IF(pti_mean(ji)==xundef)THEN
!
!    *Case where the Topographics index are not defined.
!    --------------------------------------------------------         
     zno=zno+1.0
     ptab_fsat(ji,:)=0.0     
     ptab_wtop(ji,:)=xundef
     ptab_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=pti_mean(ji)
     zti_min =pti_min(ji)
     zti_max =pti_max(ji)
     zti_std =pti_std(ji)
     zti_skew=pti_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
!
     ptab_wtop(ji,1) = zwsat_avg(ji)
     ptab_fsat(ji,1) = 1.0
     ptab_qtop(ji,1) = 0.0
!     
     ptab_wtop(ji,ipas) = zwd0_avg(ji)
     ptab_fsat(ji,ipas) = 0.0
     ptab_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)
!       
        ptab_fsat(ji,ind)=max(0.0,(zgymax-zgysat)/zftot)
!
!       Compute driest fraction
!
        zf0=max(0.0,(zgy0-zgymin)/zftot)
!
!       Calculate FMED
!        
        zfmed=(1.0-ptab_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
!
        ptab_wtop(ji,ind) = zwsat_avg(ji)-(pm(ji)*zdmoy/zd_top(ji))
!
!       Solves Qs = FMED * M * Ks * exp(-Xsat) / Ks (dimentionless)
!
        ptab_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(ptab_wtop(ji,2)==zwsat_avg(ji))THEN
!
     zfup=ptab_fsat(ji,1)
     zwup=ptab_wtop(ji,1)
     zqup=ptab_qtop(ji,1)
!  
     
     id(:)=maxloc(ptab_wtop(ji,:),ptab_wtop(ji,:)<zwsat_avg(ji))
     IF (id(1)==0) id(1) = 2
!   
     zfdown=ptab_fsat(ji,id(1))
     zwdown=ptab_wtop(ji,id(1))
     zqdown=ptab_qtop(ji,id(1))
!     
     zslopew=(zwup-zwdown)/(zfup-zfdown)   
     zslopeq=(zqup-zqdown)/(zfup-zfdown)   
!
     DO ind=2,id(1)-1
        ptab_wtop(ji,ind)=zwdown+(ptab_fsat(ji,ind)-zfdown)*zslopew
        ptab_qtop(ji,ind)=zqdown+(ptab_fsat(ji,ind)-zfdown)*zslopeq
     ENDDO
!   
   ENDIF
!
!  Lower boundary
!
   WHERE(ptab_fsat(ji,:)<=0.0      )
         ptab_wtop(ji,:)=zwd0_avg(ji)
         ptab_qtop(ji,:)=0.0
   ENDWHERE
   WHERE(ptab_wtop(ji,:)<=zwd0_avg(ji))
         ptab_fsat(ji,:)=0.0
         ptab_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