soil.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 soil( IO, KK, PK, PEK, DMI, PVEG, PCS, PFROZEN1, PFFG_NOSNOW, PFFV_NOSNOW               )  
!     ##########################################################################
!
!!****  *SOIL*  
!!
!!    PURPOSE
!!    -------
!
!     Calculates the coefficients related to the soil (i.e., CG, CT,
!     C1, C2, WGEQ) and to the snow canopy (i.e., Cs, ps, psng, psnv, and psnz0)
!         
!     
!!**  METHOD
!!    ------
!
!     Direct calculation
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Noilhan and Planton (1989)
!!    Belair (1995)
!!      
!!    AUTHOR
!!    ------
!!      S. Belair           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    13/03/95 
!!                  20/03/96      (Masson)  error in the threshold for PCG
!!                  04/09/98      (Masson)  error in C1 normalization
!!                  16/09/98      (Masson)  frozen water in the soil
!!                  07/10/98      (Masson)  new C1 formulation
!!                  26/11/98      (Boone)   C1 option (old vs new formulations)
!!                  15/03/99      (Boone)   Soil ice modifiactions: scale C1sat,
!!                                          use surface ice-weighted CG, GB method
!!                                          for dry conditions uses ZWWILT for MAX
!!                  25/03/99      (Boone)   Added Johansen (1975)/Peters-Lidard 
!!                                          option to explicitly compute CG
!!                  25/05/08     (Decharme) Added flood properties 
!!                  22/06/10     (Chauvin)  XWGMIN added as a limit value of ZWG2              
!!                                          Modification of the formula for DMI%XWGEQ                                      
!!                                          to solve numerical problem
!!                     10/10     (Decharme) The previous computation of WGEQ as ( 1.-ZX(JJ)**(IP%XPCOEF(JJ)*8.) )
!!                                          can introduced some model explosions for heavy clay soil
!!                     12/14     (LeMoigne) EBA scheme update
!-------------------------------------------------------------------------------
!
!*       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 : xpi, xci, xrholi, xday, xcl, xrholw, xcondi
USE modd_isba_par,   ONLY : xcondwtr, xwgmin
USE modd_surf_par,   ONLY : xundef
USE modd_deepsoil,   ONLY : lphysdomc
!
!
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
!
REAL, DIMENSION(:), INTENT(OUT)   :: PCS, PFROZEN1
!                                      soil and snow coefficients
!                                      PCS = heat capacity of the snow
!                                      PFROZEN1   = fraction of ice in superficial
!                                               soil
!
REAL, DIMENSION(:), INTENT(IN)   :: PFFG_NOSNOW, PFFV_NOSNOW
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PVEG))   :: ZLAMS,                         &
!                                              conductivity of snow
!
                                  zcw1max, zx2, zy1, zy2,     &
                                  zlymy1, zza, zzb, zdelta,   &
                                  za, zb,                          &
!                                              temporary variables for the 
!                                              calculation of DMI%XC1 in the case
!                                              where PWG < PWWILT (i.e., dry soils)
!
                                  zx,                              &
!                                              temporary variable for the 
!                                              calculation of DMI%XWGEQ
                                  zwsat,                           &
!                                              Wsat when ice is present in ground 
                                  zwsat1,                          &
!                                              Wsat when ice is present in surface
!                                              ground layer
                                  zwwilt,                         &
!                                              Wwilt when ice is present in ground
                                  zc1sat  
!                                              C1sat scaled due to soil ice
!
!
!                               Thermal conductivity option:                   
!                               Johansen (1975) parameters (as presented by Peters-
!                               Lidard, 1998, JAS). Used to compute CG.
!
REAL, DIMENSION(SIZE(PVEG)) ::  ZFROZEN2, ZUNFROZEN2, ZCONDSAT, ZSATDEG, ZKERSTEN, &
                                  ZCOND, ZHCAP  
!                               ZFROZEN2   = fraction of total soil layer frozen
!                               ZUNFROZEN2 = unfrozen fraction available to liquid
!                               ZCONDSAT   = saturated conductivity (water)
!                               ZSATDEG    = degree of saturation
!                               ZKERSTEN   = Kersten number
!                               ZCOND      = soil thermal conductivity (explicitly
!                                            includes soil, water and ice)
!                               ZHCAP      = Soil heat capacity
!
REAL, DIMENSION(SIZE(PVEG)) :: ZWG2
!                              ZWG2 = adjusted root-zone soil water content
!
REAL, DIMENSION(SIZE(PVEG)) :: ZCF !heat capacity of the flood
REAL, DIMENSION(SIZE(PVEG)) :: ZFF !Fraction of floodplain at the surface without snow
!
INTEGER                 :: JJ
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('SOIL',0,zhook_handle)
zwwilt(:)   = 0.
!
zfrozen2(:)  = 0.
zunfrozen2(:)= 0.
zcondsat(:)  = 0.
zsatdeg(:)   = 0.
zkersten(:)  = 0.
zcond(:)     = 0.
zhcap(:)     = 0.
!
zlams(:)     = 0.
zx(:)        = 0.
zcw1max(:)   = 0.
zy1(:)       = 0.
zx2(:)       = 0.
zy2(:)       = 0.
zlymy1(:)    = 0.
zza(:)    = 0.
zzb(:)    = 0.
zdelta(:)    = 0.
za(:)    = 0.
zb(:)    = 0.
!
zcf(:)       = xundef
!
pcs(:)       = xundef
!
!-------------------------------------------------------------------------------
!
!*       1.     FROZEN WATER FRACTION IN THE SOIL
!               ---------------------------------
!
pfrozen1(:) = 0.
WHERE (pek%XWGI(:,1) + pek%XWG(:,1) .NE. 0.) 
  pfrozen1(:) = pek%XWGI(:,1) / (pek%XWGI(:,1) + pek%XWG(:,1))
END WHERE
!
DO jj=1,SIZE(kk%XWSAT,1)
!
  zwsat(jj)    = max(kk%XWSAT(jj,1) - pek%XWGI(jj,2),xwgmin)
!
  zwsat1(jj)   = max(kk%XWSAT(jj,1) - pek%XWGI(jj,1),xwgmin)
!
  zwwilt(jj)   = kk%XWWILT(jj,1) * (zwsat1(jj) / kk%XWSAT(jj,1))
!
ENDDO
!-------------------------------------------------------------------------------
!
!*       2.     THE HEAT CAPACITY OF BARE-GROUND
!               --------------------------------
!
IF(io%CSCOND == 'NP89')THEN  
!
!                                      Actually, all the 'C' coefficients in
!                                      ISBA do not represent heat capacities,
!                                      but rather the inverse.  So in the
!                                      following formulation, CG is large
!                                      when W2 is small, thus leading to small
!                                      values for the heat capacity.  In other
!                                      words, a small amount of energy will
!                                      result in great temperature variations
!                                      (for dry soils).
!
!
! Now calculate the thermal inertia of the soil weighted
! by soil ice content (including the soil ice thermal inertia):
!
  dmi%XCG(:) = (1.-pek%XWGI(:,2)) * kk%XCGSAT(:) * ( zwsat(:)/pek%XWG(:,2) ) **( 0.5*kk%XBCOEF(:,1)/log(10.) )  &
              +  pek%XWGI(:,2)  * 2. * sqrt(xpi/(xcondi*xci*xrholi*xday))  
!
!
ELSE
!
!                                      Method of Johansen (1975) as presented by 
!                                      Peters-Lidard (JAS: 1998) for thermal
!                                      Conductivity of soil. Explicit calculation for
!                                      now (as opposed to implicit method of
!                                      Noilhan and Planton 1989). NP89 uses the
!                                      method of McCumber and Pielke (1981)
!                                      with parameters of Clapp and Hornberger (1978).
!
  DO jj=1,SIZE(pek%XWG,1)
!
! Total fraction of soil frozen:
!
    zfrozen2(jj)   = pek%XWGI(jj,2)/(pek%XWGI(jj,2) + pek%XWG(jj,2))
!
! Unfrozen fraction:
!
    zunfrozen2(jj) = (1.0-zfrozen2(jj))*kk%XWSAT(jj,1)
!
! Saturated thermal conductivity:
!
    zcondsat(jj)   = (kk%XCONDSLD(jj,1)**(1.0-kk%XWSAT(jj,1)))* (xcondi**(kk%XWSAT(jj,1)-zunfrozen2(jj)))*   &
                      (xcondwtr**zunfrozen2(jj))  
!
! Degree of saturation of soil:
!
    zsatdeg(jj)   = max(0.1, (pek%XWGI(jj,2)+pek%XWG(jj,2))/kk%XWSAT(jj,1))
!
! Kersten number:
!
    zkersten(jj)  = log10(zsatdeg(jj)) + 1.0
!
! Put in a smooth transition from thawed to frozen soils:
! simply linearly weight Kersten number by frozen fraction 
! in soil:
!
    zkersten(jj)  = (1.0-zfrozen2(jj))*zkersten(jj) +           &
                         zfrozen2(jj) *zsatdeg(jj)  
!
! Thermal conductivity of soil:
!
    zcond(jj)     = zkersten(jj)*(zcondsat(jj)-kk%XCONDDRY(jj,1)) + kk%XCONDDRY(jj,1)
!
! Heat capacity of soil:
!
    zhcap(jj)     = (1.0-kk%XWSAT(jj,1)) * kk%XHCAPSOIL(jj,1) +     &
                         pek%XWG (jj,2) * xcl * xrholw   +     &
                         pek%XWGI(jj,2) * xci * xrholi       
!
! Explicit CG calculation:
!
    dmi%XCG(jj)       = 2.*sqrt(xpi/zcond(jj)/zhcap(jj)/xday)
!
  ENDDO
!
ENDIF
!
!                                              Cg must be smaller than 2.E-5
!
dmi%XCG(:) = min( dmi%XCG(:), io%XCGMAX )
!
!-------------------------------------------------------------------------------
!
!*       4.     THE HEAT CAPACITY OF THE SNOW AND FLOOD
!               ---------------------------------------
!
WHERE (kk%XFF(:) > 0.)                                                 
       zcf(:) = 2.0 * sqrt( xpi/(xcondwtr*xrholw*xcl*xday) )
END WHERE  
!
IF(pek%TSNOW%SCHEME == 'D95' .OR. (pek%TSNOW%SCHEME == 'EBA' .AND. io%LGLACIER) )THEN
!
   WHERE (pek%XPSN(:) > 0.)
      zlams(:) = xcondi * (pek%TSNOW%RHO(:,1)/xrholw)**1.885              ! first calculate the
!                                                                   ! conductivity of snow
      pcs(:)   = 2.0 * sqrt( xpi/(zlams(:)*pek%TSNOW%RHO(:,1)*xci*xday) )
   END WHERE
!
!-------------------------------------------------------------------------------
!
!*      5.      GRID-AVERAGED HEAT CAPACITY
!               ---------------------------
!
! With contribution from the ground, vegetation, flood and snow areas
! for composite (Force-Restore) snow scheme option:
!
   dmi%XCT(:) = 1. / ( (1.-pveg(:))*(1.-pek%XPSNG(:)-kk%XFFG(:)) / dmi%XCG(:)     &
                      +  pveg(:)   *(1.-pek%XPSNV(:)-kk%XFFV(:)) / pek%XCV(:)     &
                      +                kk%XFF(:)                 / zcf(:)     &
                      +                pek%XPSN(:)               / pcs(:)     )  

!
ELSE
!
  DO jj=1,SIZE(pveg)
!
    zff(jj) = pveg(jj)*pffv_nosnow(jj) + (1.-pveg(jj))*pffg_nosnow(jj)
!
! With contribution from the ground and vegetation for explicit
! (ISBA-ES) snow scheme option:
!
     dmi%XCT(jj) = 1. / ( (1.-pveg(jj))*(1.-pffg_nosnow(jj)) / dmi%XCG(jj)     &
                      +  pveg(jj)      *(1.-pffv_nosnow(jj)) / pek%XCV(jj)     &
                      +  zff(jj)                            / zcf(jj)     )  
  ENDDO
!
ENDIF
!
!
!-------------------------------------------------------------------------------
!
!*      6.      COEFFICIENT C1
!               --------------
!                                      Scale the C1SAT coefficient as a function 
!                                      of the soil ice content
!
zc1sat(:) = pk%XC1SAT(:)*sqrt(zwsat1(:)/kk%XWSAT(:,1))
!
!
!                                      The coefficient C1 is calculated two
!                                      different ways depending on the humidity
!                                      of the soil
!
WHERE (pek%XWG(:,1) > zwwilt(:))
!                                    ! First situation:  humid soil
!                                      Then the calculation follows eq. (19)
!                                      of Noilhan and Planton(1989)
!
   dmi%XC1(:)    = zc1sat(:) * ( zwsat1(:)/pek%XWG(:,1) )**( 0.5*kk%XBCOEF(:,1) + 1 )
!
END WHERE
!
!
!                                     Calculate C1 coefficient for dry soil.
!                                     The default option is the continuous
!                                     formulation of Giard and Bazile. The
!                                     alternate approach is a discontinuous
!                                     formulation by Giordani (1993) and
!                                     Braud et al. (1993). This method
!                                     is perhaps more accurate from a physical
!                                     standpoint, as it is an explicit function
!                                     of temperature, whereas the continuous method
!                                     assumes a constant temperature. 
!
IF(io%CC1DRY=='GB93')THEN
!
  DO jj=1,SIZE(pek%XWG,1)
!  
    IF (pek%XWG(jj,1) <= zwwilt(jj)) THEN
!
!                                   ! Second situation: dry soil
!                                      We use the Gaussian formulation of
!                                      Giordanni (1993) and Braud et al. (1993)
!
!* maximum of C1 curve (computed with true Wwilt)
!
       zcw1max(jj)    = ( 1.19*zwwilt(jj)-5.09 )*pek%XTG(jj,1)  + (-146.4*zwwilt(jj)+1786.)
!
!* Giordanni (1993) and Braud et al. (1993)
!
       za(jj)         =   (-1.815e-2*pek%XTG(jj,1)+6.41)*zwwilt(jj) + (6.5e-3*pek%XTG(jj,1)-1.4)  
       zb(jj)         = za(jj)*zwwilt(jj)
       zdelta(jj)     = ( zb(jj)*zb(jj) )  / ( 2.*log( zcw1max(jj) ) )  
!
       dmi%XC1(jj) = zcw1max(jj)*(1. - 2.*pveg(jj)*( 1.-pveg(jj) ))             &
                  *exp( -(pek%XWG(jj,1)-zb(jj))*(pek%XWG(jj,1)-zb(jj)) / (2.*zdelta(jj)) )  
!
    ENDIF
!
  ENDDO
!
ELSE
!
  DO jj=1,SIZE(pek%XWG,1)
!
    IF (pek%XWG(jj,1) <= zwwilt(jj)) THEN
!
!* maximum of C1 curve (computed with true Wwilt)
!
       zcw1max(jj)    = ( 1.19*zwwilt(jj)-5.09 )*pek%XTG(jj,1)  + (-146.4*zwwilt(jj)+1786.)
!
!* C1 value at Wg = zero
!
       zy1(jj) = 10.
!
!* C1 value at Wg = wwilt
!
       zx2(jj) = zwwilt(jj)
       zy2(jj) = zc1sat(jj)*(zwsat1(jj)/zwwilt(jj))**( 0.5*kk%XBCOEF(jj,1) + 1)
!
!* correction of maximum of C1 curve for frozen soils
!
       zcw1max(jj)    = max(max(zcw1max(jj),zy2(jj)),zy1(jj))
!
!* Giard-Bazile formulation (resolution of a second order equation)
!
       zlymy1(jj) =   log( zcw1max(jj)/zy1(jj))
       zza(jj) = - log( zy2(jj)/zy1(jj))
       zzb(jj) = 2. * zx2(jj)    * zlymy1(jj)
       zdelta(jj) = 4. * (zlymy1(jj)+zza(jj)) * zlymy1(jj) * zx2(jj)**2
!
       za(jj) = (-zzb(jj)+sqrt(zdelta(jj))) / (2.*zza(jj))
!
       zb(jj) = za(jj)**2 / zlymy1(jj)
!
!
       dmi%XC1(jj) = zcw1max(jj) * exp( - (pek%XWG(jj,1)-za(jj))**2 / zb(jj) )
!
    ENDIF
!
  ENDDO
!
ENDIF
!-------------------------------------------------------------------------------
!
!*      6.      COEFFICIENT C2
!               --------------
! Including vertical diffusion limiting factor for surface soil ice:
!
IF(io%CKSAT=='SGH' .OR. io%CKSAT=='EXP')THEN
!
! Adjusted root-zone soil water content
!
  DO jj=1,SIZE(pek%XWG,1)
     zwg2(jj)=pek%XWG(jj,2)*(pk%XCONDSAT(jj,2)/pk%XCONDSAT(jj,1))**(1./(2.*kk%XBCOEF(jj,1)+3))
  ENDDO
  zwg2(:)=max(zwg2(:),xwgmin)
!
ELSE
!
   zwg2(:)=pek%XWG(:,2)
!
ENDIF
!
DO jj=1,SIZE(zwsat)
!
!Including vertical diffusion limiting factor for surface soil ice:
!
  dmi%XC2(jj) = (pk%XC2REF(jj)*zwg2(jj) / ( zwsat(jj)-zwg2(jj) + 0.01 ))           &
              *(1.0-(pek%XWGI(jj,1)/(kk%XWSAT(jj,1)-xwgmin)))  
!
!-------------------------------------------------------------------------------
!
!*      7.      EQUILIBRIUM VOLUMETRIC WATER CONTENT WGEQ
!               -----------------------------------------
!
  zx(jj) = zwg2(jj)/zwsat(jj)
!
  dmi%XWGEQ(jj) = zwg2(jj) - zwsat(jj)*kk%XACOEF(jj) *  zx(jj)**kk%XPCOEF(jj)           &
                               *( 1.-exp(kk%XPCOEF(jj)*8.*log(zx(jj))))  
!
ENDDO
!-------------------------------------------------------------------------------
!
!*      8.      SPECIAL CASE OF POLAR REGIONS
!               -----------------------------
!
IF (lphysdomc) THEN
   dmi%XCT(:) = 9.427757e-6   ! corresponds to a density of 350kg/m3 for snow
ENDIF        
IF (lhook) CALL dr_hook('SOIL',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE soil