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( HC1DRY, HSCOND, HSNOW_ISBA, OGLACIER,                     &
                         psnowrhom, pveg,                                        &
                         pcgsat, pcgmax,                                         &
                         pc1sat, pc2ref, pacoef, ppcoef, pcv,                    &
                         ppsn, ppsng, ppsnv, pffg, pffv, pff,                    &
                         pcg, pc1, pc2, pwgeq, pct, pcs, pfrozen1,               &
                         ptg, pwg, pwgi,                                         &
                         phcapsoilz, pconddryz, pcondsldz,                       &
                         pbcoef, pwsat, pwwilt,                                  &
                         hksat, pcondsat, 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 PWGEQ                                      
!!                                          to solve numerical problem
!!                     10/10     (Decharme) The previous computation of WGEQ as ( 1.-ZX(JJ)**(PPCOEF(JJ)*8.) )
!!                                          can introduced some model explosions for heavy clay soil
!!                     12/14     (LeMoigne) EBA scheme update
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
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
!
!
 CHARACTER(LEN=*),     INTENT(IN)   :: hc1dry  ! C1 for dry soil formulation
!                                             ! 'DEF' Default: Giard and Bazile
!                                             ! 'GB93' Giordani 1993, Braud 1993
!                                             ! (discontinuous at WILT)
!
 CHARACTER(LEN=*),     INTENT(IN)   :: hscond  ! thermal conductivity formulation
!                                             ! 'NP89':  Noilhan and Planton 
!                                             !  (1989: McCumber-Pielke (1981) and
!                                             !  Clapp and Hornberger (1978))
!                                             ! 'PL98' Method of Johansen (1975) as
!                                             ! presented by Peters-Lidard (JAS: 1998)
!
 CHARACTER(LEN=*),     INTENT(IN)  :: hsnow_isba ! 'DEF' = Default F-R snow scheme
!                                               !         (Douville et al. 1995)
!                                               ! '3-L' = 3-L snow scheme (option)
!                                               !         (Boone and Etchevers 2000)
!
LOGICAL, INTENT(IN)               :: oglacier   ! T = Over permanent snow and ice, 
!                                               !     initialise WGI=WSAT, Hsnow>=10m 
!                                               !     and allow 0.8<SNOWALB<0.85
!                                               ! F = No specific treatment
!                                                  
REAL, DIMENSION(:), INTENT(IN)    :: psnowrhom
!                                      Prognostic variables of ISBA at 't-dt'
!                                      PSNOWRHOM = density of snow
!
REAL,               INTENT(IN)    :: pcgmax
!                                      Maximum soil heat capacity
!
REAL, DIMENSION(:), INTENT(IN)    :: pveg, pcgsat, pc1sat
REAL, DIMENSION(:), INTENT(IN)    :: pc2ref, pacoef, ppcoef, pcv
!                                      Soil and vegetation parameters
!                                      PVEG = fraction of vegetation
!                                      PCGSAT = soil thermal coefficient at saturation
!                                      PC1SAT = value of C1 at saturation
!                                      PC2REF = reference value of C2
!                                      PACOEF, PPCOEF = a and p coefficients for
!                                                       the wgeq formulation
!                                      PCV = the heat capacity of the vegetation
!
REAL, DIMENSION(:), INTENT(IN)    :: ppsn, ppsng, ppsnv
!                                      PPSN = fraction of the grid covered by snow
!                                      PPSNG = fraction of the bare soil covered
!                                              by snow
!                                      PPSNV = fraction of the vegetation covered
!                                              by snow
!
REAL, DIMENSION(:), INTENT(IN)    :: phcapsoilz, pconddryz, pcondsldz 
!                                    PHCAPSOILZ = soil heat capacity [J/(K m3)]
!                                    PCONDDRYZ  = soil dry thermal conductivity 
!                                                 [W/(m K)] 
!                                    PCONDSLDZ  = soil solids thermal conductivity 
!                                                 [W/(m K)]
!
REAL, DIMENSION(:), INTENT(IN)    :: pbcoef, pwsat, pwwilt
!                                    PBCOEF   = b-parameter (-)
!                                    PWSAT    = porosity (m3/m3)
!                                    PWWILT   = wilting point volumetric water content(m3/m3)
!
REAL, DIMENSION(:), INTENT(IN)    :: ptg
!                                    PTG    = surface temperature (K)
!
REAL, DIMENSION(:,:), INTENT(IN)  :: pwg, pwgi 
!                                    PWG    = soil liquid water content (m3/m3)
!                                    PWGI   = soil frozen water content (m3/m3)
!
REAL, DIMENSION(:), INTENT(OUT)   :: pcg, pc1, pc2, pwgeq, pct
REAL, DIMENSION(:), INTENT(OUT)   :: pcs, pfrozen1
!                                      soil and snow coefficients
!                                      PCG = heat capacity of the bare soil
!                                      C1, C2 = coefficients for the moisture calculations
!                                      PWGEQ = equilibrium surface volumetric moisture
!                                      PCT = averaged heat capacity of the grid
!                                      PCS = heat capacity of the snow
!                                      PFROZEN1   = fraction of ice in superficial
!                                               soil
!
 CHARACTER(LEN=*),     INTENT(IN)  :: hksat      ! soil hydraulic profil option
!                                               ! 'DEF'  = ISBA homogenous soil
!                                               ! 'SGH'  = ksat exponential decay
!
REAL, DIMENSION(:,:), INTENT(IN)  :: pcondsat   ! hydraulic conductivity at saturation (m/s)
                                                ! use in the force restore scheme for
                                                ! calculate the exponential decay factor coefficients
!
REAL, DIMENSION(:), INTENT(IN)   :: pffv, pffg, pff, pffg_nosnow, pffv_nosnow
!                                   PFFG = Floodplain fraction over the ground
!                                   PFFV = Floodplain fraction over vegetation
!                                   PFF  = Floodplain fraction at the surface
!
!
!*      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 PC1 in the case
!                                              where PWG < PWWILT (i.e., dry soils)
!
                                  zx,                              &
!                                              temporary variable for the 
!                                              calculation of PWGEQ
                                  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 (pwgi(:,1) + pwg(:,1) .NE. 0.) 
  pfrozen1(:) = pwgi(:,1) / (pwgi(:,1) + pwg(:,1))
END WHERE
!
DO jj=1,SIZE(pwsat)
!
  zwsat(jj)    = max(pwsat(jj) - pwgi(jj,2),xwgmin)
!
  zwsat1(jj)   = max(pwsat(jj) - pwgi(jj,1),xwgmin)
!
  zwwilt(jj)   = pwwilt(jj) * (zwsat1(jj) / pwsat(jj))
!
ENDDO
!-------------------------------------------------------------------------------
!
!*       2.     THE HEAT CAPACITY OF BARE-GROUND
!               --------------------------------
!
IF(hscond == '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):
!
  pcg(:) = (1.-pwgi(:,2)) * pcgsat(:) * ( zwsat(:)/pwg(:,2) )                &
                                            **( 0.5*pbcoef(:)/log(10.) )       &
              +  pwgi(:,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(pwg,1)
!
! Total fraction of soil frozen:
!
    zfrozen2(jj)   = pwgi(jj,2)/(pwgi(jj,2) + pwg(jj,2))
!
! Unfrozen fraction:
!
    zunfrozen2(jj) = (1.0-zfrozen2(jj))*pwsat(jj)
!
! Saturated thermal conductivity:
!
    zcondsat(jj)   = (pcondsldz(jj)**(1.0-pwsat(jj)))*             &
                      (xcondi**(pwsat(jj)-zunfrozen2(jj)))*         &
                      (xcondwtr**zunfrozen2(jj))  
!
! Degree of saturation of soil:
!
    zsatdeg(jj)   = max(0.1, (pwgi(jj,2)+pwg(jj,2))/pwsat(jj))
!
! 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)-pconddryz(jj)) + pconddryz(jj)
!
! Heat capacity of soil:
!
    zhcap(jj)     = (1.0-pwsat(jj)) * phcapsoilz(jj) +     &
                         pwg(jj,2) * xcl * xrholw   +     &
                         pwgi(jj,2) * xci * xrholi       
!
! Explicit CG calculation:
!
    pcg(jj)       = 2.*sqrt(xpi/zcond(jj)/zhcap(jj)/xday)
!
  ENDDO
!
ENDIF
!
!                                              Cg must be smaller than 2.E-5
!
pcg(:) = min( pcg(:), pcgmax )
!
!-------------------------------------------------------------------------------
!
!*       4.     THE HEAT CAPACITY OF THE SNOW AND FLOOD
!               ---------------------------------------
!
WHERE (pff(:) > 0.)                                                 
       zcf(:) = 2.0 * sqrt( xpi/(xcondwtr*xrholw*xcl*xday) )
END WHERE  
!
IF(hsnow_isba == 'D95' .OR. (hsnow_isba == 'EBA' .AND. oglacier) )THEN
!
   WHERE (ppsn > 0.)
      zlams(:) = xcondi * (psnowrhom(:)/xrholw)**1.885              ! first calculate the
!                                                                   ! conductivity of snow
      pcs(:)   = 2.0 * sqrt( xpi/(zlams(:)*psnowrhom(:)*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:
!
   pct(:) = 1. / ( (1.-pveg(:))*(1.-ppsng(:)-pffg(:)) / pcg(:)     &
                      +  pveg(:) *(1.-ppsnv(:)-pffv(:)) / pcv(:)     &
                      +                pff(:)          / zcf(:)     &
                      +                ppsn(:)          / 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:
!
     pct(jj) = 1. / ( (1.-pveg(jj))*(1.-pffg_nosnow(jj)) / pcg(jj)     &
                      +  pveg(jj) *(1.-pffv_nosnow(jj)) / pcv(jj)     &
                      +  zff(jj)                      / zcf(jj)     )  
  ENDDO
!
ENDIF
!
!
!-------------------------------------------------------------------------------
!
!*      6.      COEFFICIENT C1
!               --------------
!                                      Scale the C1SAT coefficient as a function 
!                                      of the soil ice content
!
zc1sat(:) = pc1sat(:)*sqrt(zwsat1(:)/pwsat(:))
!
!
!                                      The coefficient C1 is calculated two
!                                      different ways depending on the humidity
!                                      of the soil
!
WHERE (pwg(:,1) > zwwilt(:))
!                                    ! First situation:  humid soil
!                                      Then the calculation follows eq. (19)
!                                      of Noilhan and Planton(1989)
!
   pc1(:)    = zc1sat(:) * ( zwsat1(:)/pwg(:,1) )**( 0.5*pbcoef(:) + 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(hc1dry=='GB93')THEN
!
  DO jj=1,SIZE(pwg,1)
!  
    IF (pwg(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 )*ptg(jj)  + (-146.4*zwwilt(jj)+1786.)
!
!* Giordanni (1993) and Braud et al. (1993)
!
       za(jj)         =   (-1.815e-2*ptg(jj)+6.41)*zwwilt(jj)            &
                          + (6.5e-3*ptg(jj)-1.4)  
       zb(jj)         = za(jj)*zwwilt(jj)
       zdelta(jj)     = ( zb(jj)*zb(jj) )  /                               &
                         ( 2.*log( zcw1max(jj) ) )  
!
       pc1(jj) = zcw1max(jj)*(1. - 2.*pveg(jj)*( 1.-pveg(jj) ))             &
                  *exp( -(pwg(jj,1)-zb(jj))*(pwg(jj,1)-zb(jj)) /              &
                  (2.*zdelta(jj)) )  
!
    ENDIF
!
  ENDDO
!
ELSE
!
  DO jj=1,SIZE(pwg,1)
!
    IF (pwg(jj,1) <= zwwilt(jj)) THEN
!
!* maximum of C1 curve (computed with true Wwilt)
!
       zcw1max(jj)    = ( 1.19*zwwilt(jj)-5.09 )*ptg(jj)  + (-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*pbcoef(jj) + 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)
!
!
       pc1(jj) = zcw1max(jj) * exp( - (pwg(jj,1)-za(jj))**2 / zb(jj) )
!
    ENDIF
!
  ENDDO
!
ENDIF
!-------------------------------------------------------------------------------
!
!*      6.      COEFFICIENT C2
!               --------------
! Including vertical diffusion limiting factor for surface soil ice:
!
IF(hksat=='SGH' .OR. hksat=='EXP')THEN
!
! Adjusted root-zone soil water content
!
  DO jj=1,SIZE(pwg,1)
     zwg2(jj)=pwg(jj,2)*(pcondsat(jj,2)/pcondsat(jj,1))**(1./(2.*pbcoef(jj)+3))
  ENDDO
  zwg2(:)=max(zwg2(:),xwgmin)
!
ELSE
!
   zwg2(:)=pwg(:,2)
!
ENDIF
!
DO jj=1,SIZE(zwsat)
!
!Including vertical diffusion limiting factor for surface soil ice:
!
  pc2(jj) = (pc2ref(jj)*zwg2(jj) / ( zwsat(jj)-zwg2(jj) + 0.01 ))           &
              *(1.0-(pwgi(jj,1)/(pwsat(jj)-xwgmin)))  
!
!-------------------------------------------------------------------------------
!
!*      7.      EQUILIBRIUM VOLUMETRIC WATER CONTENT WGEQ
!               -----------------------------------------
!
  zx(jj) = zwg2(jj)/zwsat(jj)
!
  pwgeq(jj) = zwg2(jj) - zwsat(jj)*pacoef(jj)                      & 
                               *     zx(jj)** ppcoef(jj)           &
                               *( 1.-exp(ppcoef(jj)*8.*log(zx(jj))))  
!
ENDDO
!-------------------------------------------------------------------------------
!
!*      8.      SPECIAL CASE OF POLAR REGIONS
!               -----------------------------
!
IF (lphysdomc) THEN
   pct(:) = 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