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(IO, KK, PK, PEK, DMI, PVEG, PFROZEN1, PFFG, PFFV, PSOILCONDZ, PSOILHCAPZ   )
!     ##########################################################################
!
!!****  *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_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_k_t, isba_p_t, isba_pe_t
USE modd_diag_misc_isba_n, ONLY : diag_misc_isba_t
!
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
!
!
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_k_t), INTENT(INOUT) :: KK
TYPE(isba_p_t), INTENT(INOUT) :: PK
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
TYPE(diag_misc_isba_t), INTENT(INOUT) :: DMI
!
REAL, DIMENSION(:), INTENT(IN)    :: PVEG
!                                      Soil and vegetation parameters
!                                      PVEG = fraction of vegetation

!
REAL, DIMENSION(:), INTENT(OUT)   :: PFROZEN1
!                                      PFROZEN1 = fraction of ice in superficial soil
!
REAL, DIMENSION(:), INTENT(IN)   :: PFFV, PFFG
!                                   PFFG = Floodplain fraction over the ground
!                                   without snow (ES)
!                                   PFFV = Floodplain fraction over vegetation
!                                   without snow (ES)
!
REAL, DIMENSION(:,:), INTENT(OUT) :: PSOILCONDZ, PSOILHCAPZ
!                                    PSOILHCAP = soil heat capacity        (J m-3 K-1)
!                                    PSOILCOND = soil thermal conductivity (W m-1 K-1)
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PEK%XTG,1),SIZE(PEK%XTG,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(pek%XWG,1)
inl=SIZE(pek%XWG,2)
!
zff(:) = 0.0
zcf(:) = xundef
!
!-------------------------------------------------------------------------------
!
!*       1.     WATER TABLE DETH ADJUSTMENT FOR ISBA (m)
!               -----------------------------------------
!
WHERE(kk%XWTD(:)==xundef)
! no water table / surface coupling over some regions        
  zwtd(:) = xwtd_maxdepth 
ELSEWHERE
  zwtd(:) = kk%XFWTD(:)/max(-kk%XWTD(:),0.001) + (1.0-kk%XFWTD(:))/max(-kk%XWTD(:),xwtd_maxdepth)
  zwtd(:) = 1.0/zwtd(:)
ENDWHERE
!
!-------------------------------------------------------------------------------
!
!*       2.     MATRIC POTENTIAL AND MOISTURE EXTRAPOLATION
!               -------------------------------------------
!
DO jl=1,inl
   DO jj=1,ini
!   
      idepth=pk%NWG_LAYER(jj)
      IF(jl>idepth)THEN
!                           
!       total matric potential
        zwork1  = min(1.0,(pek%XWG(jj,idepth)+pek%XWGI(jj,idepth))/kk%XWSAT(jj,idepth))
        zlog    = kk%XBCOEF(jj,idepth)*log(zwork1)
        zmatpot(jj,idepth) = kk%XMPOTSAT(jj,idepth)*exp(-zlog)

!       extrapolation of total matric potential
        zwork1         = 0.5*(pk%XDG(jj,idepth) + pk%XDG(jj,idepth-1))
        zwork2         = 0.5*(pk%XDG(jj,jl    ) + pk%XDG(jj,jl-1))
        zwork3         = max(0.0,(zwtd(jj)-zwork2)/(zwork2-zwork1))
        zmatpot(jj,jl) = (kk%XMPOTSAT(jj,jl)+zwork3*zmatpot(jj,idepth))/(1.0+zwork3)
!
!       total soil water content computation
        zwork1      = max(1.0,zmatpot(jj,jl)/kk%XMPOTSAT(jj,jl))
        zlog        = log(zwork1)/kk%XBCOEF(jj,jl)
        zwtot       = kk%XWSAT(jj,jl)*exp(-zlog)
        zwtot       = max(xwgmin,zwtot)
!
!       soil liquid water content computation
        zmatpot(jj,jl) = min(kk%XMPOTSAT(jj,jl),xlmtt*(pek%XTG(jj,jl)-xtt)/(xg*pek%XTG(jj,jl)))
!        
        zwork1      = max(1.0,zmatpot(jj,jl)/kk%XMPOTSAT(jj,jl))
        zlog        = log(zwork1)/kk%XBCOEF(jj,jl)
        zwl         = kk%XWSAT(jj,jl)*exp(-zlog)
        zwl         = max(zwl,xwgmin)
        pek%XWG (jj,jl) = min(zwl,zwtot )
!        
!       soil ice computation        
        pek%XWGI(jj,jl) = max(0.0,zwtot-pek%XWG(jj,jl))
!        
      ENDIF
   ENDDO
ENDDO
!
!-------------------------------------------------------------------------------
!
!*       3.     SURFACE FROZEN FRACTION
!               -----------------------
!
!
! Surface soil water reservoir frozen fraction:
!
pfrozen1(:) = pek%XWGI(:,1)/(pek%XWGI(:,1) + max(pek%XWG(:,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(:,:) = kk%XCONDDRY (:,:)
zcondsldz(:,:) = kk%XCONDSLD (:,:)
!
IF(io%CDIFSFCOND == 'MLCH') THEN
!  
  DO jl=1,inl
     DO jj=1,ini  
!
        zvegmulch(jj,jl) = pveg(jj)*min(pk%XDZG(jj,jl),max(0.0,zthickm-pk%XDG(jj,jl)+pk%XDZG(jj,jl)))/pk%XDZG(jj,jl)     
!
        IF(zvegmulch(jj,jl)>0.0)THEN
           zconddryz(jj,jl) = 1.0/((1.0-zvegmulch(jj,jl))/kk%XCONDDRY(jj,jl)+zvegmulch(jj,jl)/xomconddry)
           zcondsldz(jj,jl) = 1.0/((1.0-zvegmulch(jj,jl))/kk%XCONDSLD(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   = pek%XWGI(jj,jl)/(pek%XWGI(jj,jl) + max(pek%XWG(jj,jl),xwgmin))
      zunfrozen2df = (1.0-zfrozen2df)*kk%XWSAT(jj,jl)
!
!Old: CONDSATDF=(CONDSLDZ**(1.0-WSAT))*(CONDI**(WSAT-UNFROZEN2DF))*(CONDWTR**UNFROZEN2DF)  
      zwork1      = log(zcondsldz(jj,jl))*(1.0-kk%XWSAT(jj,jl))
      zwork2      = zlog_condi*(kk%XWSAT(jj,jl)-zunfrozen2df)
      zwork3      = zlog_condwtr*zunfrozen2df
      zcondsatdf  = exp(zwork1+zwork2+zwork3)
!
      zsatdegdf   = max(0.1, (pek%XWGI(jj,jl)+pek%XWG(jj,jl))/kk%XWSAT(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-kk%XWSAT(jj,jl))*kk%XHCAPSOIL(jj,jl) +       &
                           pek%XWG  (jj,jl) *xcl*xrholw + pek%XWGI (jj,jl) *xci*xrholi    
   ENDDO
ENDDO
!
! Surface soil thermal inertia [(m2 K)/J]
!
dmi%XCG(:) = 1.0 / ( pk%XDG(:,1) * psoilhcapz(:,1) )
!
dmi%XCG(:) = min(zctmax,dmi%XCG(:))
!
!-------------------------------------------------------------------------------
!
!*       7.     THE HEAT CAPACITY OF VEGETATION
!               --------------------------------
!
! Vegetation thermal inertia [(m2 K)/J]
!
zcv(:) = 1.0 / ( xcvheatf/pek%XCV(:) +  xcl * pek%XWR(:) )
!
zcv(:) = min(zctmax,zcv(:))
!
!-------------------------------------------------------------------------------
!
!*       8.     THE HEAT CAPACITY OF FLOOD
!               --------------------------------
!
IF(io%LFLOOD)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):
!
dmi%XCT(:) = 1. / ( (1.-pveg(:))*(1.-pffg(:)) / dmi%XCG(:)     &
                 +  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=pk%NWG_LAYER(jj)
      IF(jl>idepth)THEN
        pek%XWG (jj,jl) = xundef
        pek%XWGI(jj,jl) = xundef
      ENDIF
   ENDDO
ENDDO
!
IF (lhook) CALL dr_hook('SOILDIF',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE soildif