soildif.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 soildif(HDIFSFCOND, OFLOOD,                                     &
                         pveg, pcv, pffg, pffv,                                  &
                         pcg, pct, pfrozen1,                                     &
                         pd_g, pdzg, ptg, pwg, pwgi, kwg_layer,                  &
                         phcapsoilz, pconddryz, pcondsldz,                       &
                         pbcoef, pwsat, pmpotsat, psoilcondz, psoilhcapz,        &
                         pfwtd, pwtd, pwr                                        )
!     ##########################################################################
!
!!****  *SOIL*  
!!
!!    PURPOSE
!!    -------
!
!     Calculates the coefficients related to the soil (i.e., surface heat capacities, CG, CT,
!     and thermal conductivity and heat capacity profiles)
!         
!     
!!**  METHOD
!!    ------
!
!     Direct calculation
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    USE MODD_CST
!!    USE MODD_PARAMETERS
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Noilhan and Planton (1989)
!!    Belair (1995)
!!    Boone (2000)
!!      
!!    AUTHOR
!!    ------
!!      A. Boone           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    
!!                  25/03/99      (Boone)   Added Johansen (1975)/Peters-Lidard 
!!                                          option to explicitly compute CG
!!                  08/25/02      (Boone)   DIF option code only
!!                  25/05/08     (Decharme) Add Flood properties
!!                  03/08/11     (Decharme) Optimization
!!                     04/13     (Decharme) good soil moisture extrapolation computation
!!                  23/07/13     (Decharme) Surface / Water table depth coupling
!!                  23/10/14     (Decharme) revise all thermo properties
!!                                          delete NP89 option for thermal cond
!!                                          because not physical with explicit soil.
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_csts,       ONLY : xcl, xci, xrholw, xrholi, xpi, xday, xcondi, xtt, xlmtt, xg
USE modd_isba_par,   ONLY : xcondwtr, xwgmin, xwtd_maxdepth, & 
                            xomrho, xomsph, xomconddry,      &
                            xomcondsld, xcvheatf
!
USE modd_surf_par,   ONLY : xundef
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
 CHARACTER(LEN=*),     INTENT(IN)  :: hdifsfcond ! NOTE: Only used when HISBA = DIF
!                                               ! MLCH' = include the insulating effect of leaf
!                                               !         litter/mulch on the surface thermal cond.
!                                               ! 'DEF' = no mulch effect
!
LOGICAL, INTENT(IN)               :: oflood ! Flood scheme 
!
REAL, DIMENSION(:), INTENT(IN)    :: pveg, pfwtd, pwtd, pcv, pwr
!                                      Soil and vegetation parameters
!                                      PVEG = fraction of vegetation
!                                      PFWTD= grid-cell fraction of water table to rise
!                                      PWTD = water table depth (negative below soil surface)
!                                      PCV  = the heat capacity of the vegetation
!                                      PWR  = canopy intercepted water
!
REAL, DIMENSION(:,:), INTENT(IN)  :: phcapsoilz, pconddryz, pcondsldz, pd_g, pdzg
!                                    PHCAPSOILZ = soil heat capacity [J/(K m3)]
!                                    PCONDDRYZ  = soil dry thermal conductivity 
!                                                 [W/(m K)] 
!                                    PCONDSLDZ  = soil solids thermal conductivity 
!                                                 [W/(m K)]
!                                    PD_G       = soil layer depth [m]
!                                    PDZG       = soil layers thicknesses [m]
!
REAL, DIMENSION(:,:), INTENT(IN)  :: pbcoef, pwsat, pmpotsat, ptg
!                                    PBCOEF   = profile of b-parameter (-)
!                                    PWSAT    = profile of porosity (m3/m3)
!                                    PMPOTSAT = profile of matric potential at saturation (m)
!
REAL, DIMENSION(:,:),INTENT(INOUT):: pwg, pwgi
!                                    PWG    = soil liquid water content (m3/m3)
!                                    PWGI   = soil frozen water content (m3/m3)
!
INTEGER, DIMENSION(:), INTENT(IN) :: kwg_layer  
!                                    KWG_LAYER = Number of soil moisture layers (DIF option)
!
REAL, DIMENSION(:), INTENT(OUT)   :: pfrozen1, pcg, pct
!                                      PFROZEN1 = fraction of ice in superficial soil
!                                      PCT      = averaged surface heat capacity of the grid (m2 K J-1)
!                                      PCG      = averaged surface soil heat capacity (m2 K J-1)
!
REAL, DIMENSION(:,:), INTENT(OUT) :: psoilcondz, psoilhcapz
!                                    PSOILHCAP = soil heat capacity        (J m-3 K-1)
!                                    PSOILCOND = soil thermal conductivity (W m-1 K-1)
!
!
REAL, DIMENSION(:), INTENT(IN)   :: pffv, pffg
!                                   PFFG = Floodplain fraction over the ground
!                                   without snow (ES)
!                                   PFFV = Floodplain fraction over vegetation
!                                   without snow (ES)
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PTG,1),SIZE(PTG,2)) :: zmatpot, zconddryz, zcondsldz, zvegmulch
!                                           ZMATPOT    = soil matric potential (m)
!
REAL                         :: zfrozen2df, zunfrozen2df, zcondsatdf, zlog_condi, zlog_condwtr,  &
                                zsatdegdf, zkerstendf, zwork1, zwork2, zwork3, zlog, zwtot, zwl
!    
REAL, PARAMETER              :: zctmax = 1.e-4 ! Maximum thermal inertia
!
REAL, PARAMETER              :: zthickm = 0.04 ! Mulch thickness (m)
!
REAL, DIMENSION(SIZE(PVEG)) :: zff, zcf, zcv !Thermal inertia of the flood or vegetation
!
REAL, DIMENSION(SIZE(PVEG)) :: zwtd ! Water table depth if no coupling (m)  
!
INTEGER :: ini, inl, jj, jl, idepth
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
!*       0.     Initialization:
!               ---------------
!
IF (lhook) CALL dr_hook('SOILDIF',0,zhook_handle)
!
ini=SIZE(pwg,1)
inl=SIZE(pwg,2)
!
zff(:) = 0.0
zcf(:) = xundef
!
!-------------------------------------------------------------------------------
!
!*       1.     WATER TABLE DETH ADJUSTMENT FOR ISBA (m)
!               -----------------------------------------
!
WHERE(pwtd(:)==xundef)
! no water table / surface coupling over some regions        
  zwtd(:) = xwtd_maxdepth 
ELSEWHERE
  zwtd(:) = pfwtd(:)/max(-pwtd(:),0.001) + (1.0-pfwtd(:))/max(-pwtd(:),xwtd_maxdepth)
  zwtd(:) = 1.0/zwtd(:)
ENDWHERE
!
!-------------------------------------------------------------------------------
!
!*       2.     MATRIC POTENTIAL AND MOISTURE EXTRAPOLATION
!               -------------------------------------------
!
DO jl=1,inl
   DO jj=1,ini
!   
      idepth=kwg_layer(jj)
      IF(jl>idepth)THEN
!                           
!       total matric potential
        zwork1  = min(1.0,(pwg(jj,idepth)+pwgi(jj,idepth))/pwsat(jj,idepth))
        zlog    = pbcoef(jj,idepth)*log(zwork1)
        zmatpot(jj,idepth) = pmpotsat(jj,idepth)*exp(-zlog)

!       extrapolation of total matric potential
        zwork1         = 0.5*(pd_g(jj,idepth)+pd_g(jj,idepth-1))
        zwork2         = 0.5*(pd_g(jj,jl)+pd_g(jj,jl-1))
        zwork3         = max(0.0,(zwtd(jj)-zwork2)/(zwork2-zwork1))
        zmatpot(jj,jl) = (pmpotsat(jj,jl)+zwork3*zmatpot(jj,idepth))/(1.0+zwork3)
!
!       total soil water content computation
        zwork1      = max(1.0,zmatpot(jj,jl)/pmpotsat(jj,jl))
        zlog        = log(zwork1)/pbcoef(jj,jl)
        zwtot       = pwsat(jj,jl)*exp(-zlog)
        zwtot       = max(xwgmin,zwtot)
!
!       soil liquid water content computation
        zmatpot(jj,jl) = min(pmpotsat(jj,jl),xlmtt*(ptg(jj,jl)-xtt)/(xg*ptg(jj,jl)))
!        
        zwork1      = max(1.0,zmatpot(jj,jl)/pmpotsat(jj,jl))
        zlog        = log(zwork1)/pbcoef(jj,jl)
        zwl         = pwsat(jj,jl)*exp(-zlog)
        zwl         = max(zwl,xwgmin)
        pwg(jj,jl) = min(zwl,zwtot )
!        
!       soil ice computation        
        pwgi(jj,jl) = max(0.0,zwtot-pwg(jj,jl))
!        
      ENDIF
   ENDDO
ENDDO
!
!-------------------------------------------------------------------------------
!
!*       3.     SURFACE FROZEN FRACTION
!               -----------------------
!
!
! Surface soil water reservoir frozen fraction:
!
pfrozen1(:) = pwgi(:,1)/(pwgi(:,1) + max(pwg(:,1),xwgmin))
!
!-------------------------------------------------------------------------------
!
!*       4.     SIMPLE LITTER/MULCH EFFECT
!               --------------------------
!
! This takes into account the insulating effect of dead vegetation/leaf litter/mulch on
! uppermost soil layers thermal properties. Use organic matter thermal properties.
!
!
zconddryz(:,:) = pconddryz(:,:)
zcondsldz(:,:) = pcondsldz(:,:)
!
IF(hdifsfcond == 'MLCH') THEN
!  
  DO jl=1,inl
     DO jj=1,ini  
!
        zvegmulch(jj,jl) = pveg(jj)*min(pdzg(jj,jl),max(0.0,zthickm-pd_g(jj,jl)+pdzg(jj,jl)))/pdzg(jj,jl)     
!
        IF(zvegmulch(jj,jl)>0.0)THEN
           zconddryz(jj,jl) = 1.0/((1.0-zvegmulch(jj,jl))/pconddryz(jj,jl)+zvegmulch(jj,jl)/xomconddry)
           zcondsldz(jj,jl) = 1.0/((1.0-zvegmulch(jj,jl))/pcondsldz(jj,jl)+zvegmulch(jj,jl)/xomcondsld)
        ENDIF
!
     ENDDO
  ENDDO
!
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       5.     THE THERMAL CONDUCTIVITY OF BARE-GROUND
!               ---------------------------------------
!
! Calculate thermal conductivity using PL98 :
!
zlog_condi   = log(xcondi)
zlog_condwtr = log(xcondwtr)
!
DO jl=1,inl
   DO jj=1,ini
!     
      zfrozen2df   = pwgi(jj,jl)/(pwgi(jj,jl) + max(pwg(jj,jl),xwgmin))
      zunfrozen2df = (1.0-zfrozen2df)*pwsat(jj,jl)
!
!Old: CONDSATDF=(CONDSLDZ**(1.0-WSAT))*(CONDI**(WSAT-UNFROZEN2DF))*(CONDWTR**UNFROZEN2DF)  
      zwork1      = log(zcondsldz(jj,jl))*(1.0-pwsat(jj,jl))
      zwork2      = zlog_condi*(pwsat(jj,jl)-zunfrozen2df)
      zwork3      = zlog_condwtr*zunfrozen2df
      zcondsatdf  = exp(zwork1+zwork2+zwork3)
!
      zsatdegdf   = max(0.1, (pwgi(jj,jl)+pwg(jj,jl))/pwsat(jj,jl))
      zsatdegdf   = min(1.0,zsatdegdf)
      zkerstendf  = log10(zsatdegdf) + 1.0
      zkerstendf  = (1.0-zfrozen2df)*zkerstendf + zfrozen2df *zsatdegdf  
!
! Thermal conductivity of soil:
!
      psoilcondz(jj,jl) = zkerstendf*(zcondsatdf-zconddryz(jj,jl)) + zconddryz(jj,jl)  
!
   ENDDO
ENDDO
!
!-------------------------------------------------------------------------------
!
!*       6.     THE HEAT CAPACITY OF BARE-GROUND
!               --------------------------------
!
! Soil Heat capacity [J/(m3 K)]
!
DO jl=1,inl
   DO jj=1,ini
      psoilhcapz(jj,jl) = (1.0-pwsat(jj,jl))*phcapsoilz(jj,jl) +         &
                               pwg(jj,jl) *xcl*xrholw        +         &
                               pwgi(jj,jl) *xci*xrholi    
   ENDDO
ENDDO
!
! Surface soil thermal inertia [(m2 K)/J]
!
pcg(:) = 1.0 / ( pd_g(:,1) * psoilhcapz(:,1) )
!
pcg(:) = min(zctmax,pcg(:))
!
!-------------------------------------------------------------------------------
!
!*       7.     THE HEAT CAPACITY OF VEGETATION
!               --------------------------------
!
! Vegetation thermal inertia [(m2 K)/J]
!
zcv(:) = 1.0 / ( xcvheatf/pcv(:) +  xcl * pwr(:) )
!
zcv(:) = min(zctmax,zcv(:))
!
!-------------------------------------------------------------------------------
!
!*       8.     THE HEAT CAPACITY OF FLOOD
!               --------------------------------
!
IF(oflood)THEN
!
  zff(:) = pveg(:)*pffv(:) + (1.-pveg(:))*pffg(:)
!
  WHERE(zff(:)>0.0)
    zcf(:) = 2.0 * sqrt( xpi/(xcondwtr*xrholw*xcl*xday) )
  ENDWHERE
!
ENDIF
!
!-------------------------------------------------------------------------------
!
!*      9.      GRID-AVERAGED HEAT CAPACITY
!               ---------------------------
!
! With contribution from the ground, flood and vegetation for explicit
! (ISBA-ES) snow scheme option (i.e. no snow effects included here):
!
pct(:) = 1. / ( (1.-pveg(:))*(1.-pffg(:)) / pcg(:)     &
                 +  pveg(:) *(1.-pffv(:)) / zcv(:)     &
                 +  zff(:)               / zcf(:)     )  
!
!-------------------------------------------------------------------------------
!
!*      10.     RESTORE DEFAULT VALUES
!               ----------------------
!
! restore default moisture and ice values under moisture soil depth
!
DO jl=1,inl
   DO jj=1,ini
      idepth=kwg_layer(jj)
      IF(jl>idepth)THEN
        pwg(jj,jl) = xundef
        pwgi(jj,jl) = xundef
      ENDIF
   ENDDO
ENDDO
!
IF (lhook) CALL dr_hook('SOILDIF',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE soildif