ice_soilfr.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 ice_soilfr(HSNOW_ISBA, HSOILFRZ, PTSTEP, PKSFC_IVEG, PCG, PCT, &
                            ppsng, pffg, ptauice, pdwgi1, pdwgi2, pwsatz,       &
                            pmpotsatz, pbcoefz, pd_g, ptg, pwgi, pwg            )   
!!     ##########################################################################
!
!!****  *ICE_SOILFR*  
!!
!!    PURPOSE
!!    -------
!
!     In ISBA-FR: calculates the 
!     1.) evolution of the surface and deep-soil temperature(s)
!         (i.e., Ts and T2 if Force-Restore, TN if DIFfusion) due to soil water 
!         phase changes
!     
!!**  METHOD
!!    ------
!
!     - latent heating from soil ice phase changes
!
!!    EXTERNAL
!!    --------
!!
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Noilhan and Planton (1989)
!!    Belair (1995)
!!    Boone et al. (2000)
!!      
!!    AUTHOR
!!    ------
!!
!!      A. Boone           * Meteo-France *
!!      B. Decharme        * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original      14/03/95 
!!      (A.Boone)     08/11/00 soil ice phase changes herein
!!      (A.Boone)     06/05/02 Updates, ordering. Addition of 'HSOILFRZ' option
!!      (B. Decharme) 03/2009  BUG : effect of insolation due to vegetation cover
!!                                  at 1 for bare soil
!!      (B. Decharme) 07/2012  Time spliting for soil ice
!!      (A.Boone)     02/2013  Split from isba_fluxes.F90
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_csts,       ONLY : xcl, xtt, xpi, xday, xci, xrholi,     &
                            xlmtt, xrholw, xg, xcondi
USE modd_surf_par,   ONLY : xundef
USE modd_isba_par,   ONLY : xwgmin, xsphsoil, xdrywght
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
 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 2001)
!
 CHARACTER(LEN=*),   INTENT(IN)      :: hsoilfrz   ! soil freezing-physics option
!                                                 ! 'DEF'   Default (Boone et al. 2000; Giard and Bazile 2000)
!                                                 ! 'LWT'   phase changes as above, but relation between unfrozen 
!                                                         water and temperature considered
!
REAL, INTENT (IN)                   :: ptstep     ! model time step (s)
!
!
REAL, DIMENSION(:), INTENT(IN)      :: pksfc_iveg
!                                      PKSFC_IVEG= non-dimensional vegetation insolation coefficient
!
REAL, DIMENSION(:), INTENT(IN)      :: ptauice
!                                      PTAUICE = characteristic time scale for soil water phase changes (s)
!
REAL, DIMENSION(:),  INTENT(IN)     :: ppsng
!                                      PPSNG = snow fractions over ground
!
REAL, DIMENSION(:),  INTENT(IN)     :: pffg
!                                      PPSNG = flood fractions over ground
!
REAL, DIMENSION(:), INTENT (IN)     :: pcg, pct
!                                      PCT    = area-averaged heat capacity (K m2 J-1)
!                                      PCG    = heat capacity of the soil (K m2 J-1)
!
REAL, DIMENSION(:,:), INTENT(IN)    :: pd_g, pwsatz
!                                      PD_G   = Depth of bottom of Soil layers (m)
!                                      PWSATZ    = porosity (m3/m3)
!
REAL, DIMENSION(:,:), INTENT(IN)    :: pmpotsatz, pbcoefz
!                                      PMPOTSATZ = matric potential at saturation (m)
!                                      PBCOEFZ   = slope of the water retention curve (-)
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: pwg, pwgi, ptg 
!                                      PWGI   = soil frozen volumetric water content (m3/m3)
!                                      PWG    = soil liquid volumetric water content (m3/m3)
!                                      PTG    = soil temperature profile (K)
!
REAL, DIMENSION(:), INTENT(OUT)     :: pdwgi1, pdwgi2
!                                      PDWGI1   = near-surface liquid water equivalent
!                                                 volumetric ice content tendency
!                                      PDWGI2   = deep ground liquid water equivalent
!                                                 volumetric ice content tendency
!
!
!*      0.2    declarations of local variables
!
REAL                        ::   zksoil     ! coefficient for soil freeze/thaw
!
REAL, DIMENSION(SIZE(PCG)) ::   zksfc_frz, zfreezing, zice_melt, zwim,       &
                                 zwit, zwgi1, zwgi2, zwm, zsoilheatcap,       &
                                 ziceeff, zeffic, ztauice,                    &
                                 zwgmin, ztgmax, zmatpot, zdeltat
!                                ZKSFC_FRZ    = surface insolation coefficient (kg m-2 K-1)
!                                ZFREEZING    = rate for freezing soil water (kg m-2)
!                                ZICE_MELT    = rate for melting soil ice (kg m-2)
!                                ZWIM,ZWIT    = available ice content (m3 m-3)
!                                ZWGI1,ZWGI2  = volumetric ice contents (m3 m-3)
!                                ZWM          = available liquid water content (m3 m-3)
!                                ZSOILHEATCAP = soil heat capacity (J  m-3 K-1)
!                                ZICEEFF      = effective soil ice penetration depth (m)
!                                ZEFFIC       = phase change efficiency
!                                ZMATPOT      = soil matric potential (m)
!                                ZWGMIN       = volumetric water content above which soil water can
!                                               be unfrozen (if energy and mass available)(m3 m-3)
!                                ZTGMAX       = temperature below which liquid water 
!                                               can be frozen (if energy and mass available)(K)
!                                ZDELTAT      = Freezing or melting temperature depression (K) after 
!                                               possible flux correction
!
REAL, DIMENSION(SIZE(PCG)) ::  zwsat_avgz
!                               ZWSAT_AVGZ = soil column average porosity (m3 m-3)
!
REAL, DIMENSION(SIZE(PCG)) :: zpsng
!                               ZPSNG = snow fractions corresponding to
!                                       dummy argument PPSNG
!                                       if HSNOW_ISBA = 'DEF' (composite
!                                       or Force-Restore snow scheme), else
!                                       they are zero for explicit snow case
!                                       as snow fluxes calculated outside of
!                                       this routine using the 
!                                       HSNOW_ISBA = '3-L' or 'CRO' option.
!
!*      0.3    declarations of local parameters
!
REAL, PARAMETER             :: zinsolfrz_veg = 0.20  ! (-)       Vegetation insolation coefficient
!
REAL, PARAMETER             :: zinsolfrz_lai = 30.0  ! (m2 m-2)  Vegetation insolation coefficient
!
REAL, PARAMETER             :: zeffic_min    = 0.01  ! (-)   (0 <= ZEFFIC_MIN << 1)
!                                                                This parameter ensures
!                                                                a small minimum melt or freeze efficiency...
!                                                                It is numerical. When it is small, it has
!                                                                a only small impact on results, except
!                                                                that it keeps very small values of ice from persisting
!                                                                over long periods of time as they approach zero.
!                                                                If it is zero, then this effect off.
!
!
INTEGER         :: ini, jj
!
REAL, DIMENSION(SIZE(PCG))          :: zwork1, zwork2, ztdiurn
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
!*       0.     Initialization
!               --------------
IF (lhook) CALL dr_hook('ICE_SOILFR',0,zhook_handle)
!
zfreezing(:)    = 0.0
zksfc_frz(:)    = 0.0
zeffic(:)       = 0.0
zice_melt(:)    = 0.0
zwgi1(:)        = 0.0
zwim(:)         = 0.0
zsoilheatcap(:) = 0.0
zwit(:)         = 0.0
zwgi2(:)        = 0.0
ztgmax(:)       = 0.0
zwgmin(:)       = 0.0
zmatpot(:)      = 0.0
zwsat_avgz(:)   = xundef
zdeltat(:)      = 0.0
ztdiurn(:)      = 0.0
!
ini = SIZE(ptg,1)
!
!-------------------------------------------------------------------------------
!
! If ISBA-ES option in use, then snow covered surface
! fluxes calculated outside of this routine, so set
! the local snow fractions here to zero:
! 
IF(hsnow_isba == '3-L' .OR. hsnow_isba == 'CRO')THEN
   zpsng(:)     = 0.0
ELSE
   zpsng(:)     = ppsng(:)+pffg(:)
ENDIF
!
!*       1.    Melting/freezing normalized coefficient
!               ---------------------------------------
!
zksoil       = (0.5 * sqrt(xcondi*xci*xrholi*xday/xpi))/xlmtt
!
ztauice(:) = max(ptstep,ptauice(:))
!
DO jj=1,ini
!-------------------------------------------------------------------------------
!*       2.     EFFECT OF THE MELTING/FREEZING 
!               ON THE SURFACE-SOIL HEAT AND ICE CONTENTS
!               ('DEF' or Force-Restore soil option)
!               -----------------------------------------
!
!        2.0    Average soil-column porosity
!               ----------------------------
!               if Force-Restore option in use, then vertical
!               profiles of soil hydrological parameters are constant,
!               so use the values in uppermost element (arbitrary)
!
    zwsat_avgz(jj) = pwsatz(jj,1)
!
!               Influence of vegetation insolation on surface:
!
    zksfc_frz(jj) = zksoil * pksfc_iveg(jj)
!
ENDDO
!*       2.2    Water freezing
!               --------------
!
IF(hsoilfrz == 'LWT')THEN
!
! use option to control phase changes based on a relationship
! between the unfrozen liquid water content and temperature.
! Uses the Clapp and Hornberger model for water potential.
! The energy-limit method used by Boone et al. 2000 and
! Giard and Bazile (2000) is the default. 
!
  DO jj=1,ini
      zmatpot(jj)   = min(pmpotsatz(jj,1), xlmtt*(ptg(jj,1)-xtt)/(xg*ptg(jj,1)) )
      zwgmin(jj)    = zwsat_avgz(jj)*( (zmatpot(jj)/pmpotsatz(jj,1))**(-1./pbcoefz(jj,1)) )

      zmatpot(jj)   = pmpotsatz(jj,1)*( (pwg(jj,1)/zwsat_avgz(jj))**(-pbcoefz(jj,1)) )
      ztgmax(jj)    = xlmtt*xtt/(xlmtt - xg* zmatpot(jj))
  ENDDO
ELSE
    zwgmin(:)    = xwgmin
    ztgmax(:)    = xtt
ENDIF
!
DO jj=1,ini
! 
    zdeltat(jj)  = ptg(jj,1) - ztgmax(jj) ! initial temperature depression
!
    zwork2(jj) = xrholw*pd_g(jj,1)
    zeffic(jj)    = max(zeffic_min,(pwg(jj,1)-xwgmin)/zwsat_avgz(jj))
    zfreezing(jj) = min( max(0.0,pwg(jj,1)-zwgmin(jj))*zwork2(jj),    &  
                  zksfc_frz(jj)*zeffic(jj)*max( -zdeltat(jj), 0.) )
!
!*       2.3    Ground Ice melt
!               ---------------
!
    zeffic(jj)    =  max(zeffic_min,pwgi(jj,1)/(zwsat_avgz(jj)-xwgmin))
    zice_melt(jj) = min( pwgi(jj,1)*zwork2(jj),                      &
                  zksfc_frz(jj)*zeffic(jj)*max( zdeltat(jj), 0. ) )
!
!*       2.4    Ice reservoir evolution
!               -----------------------
!
! Melting of ice/freezing of water:
!
    zwgi1(jj) = pwgi(jj,1) + (ptstep/ztauice(jj))*(1.0-zpsng(jj))*        &
              (zfreezing(jj) - zice_melt(jj))/zwork2(jj) 
!
!
    zwgi1(jj)  = max( zwgi1(jj) , 0.             )
    zwgi1(jj)  = min( zwgi1(jj) , zwsat_avgz(jj)-xwgmin)
!
! Time tendency:
!
    pdwgi1(jj) = zwgi1(jj) - pwgi(jj,1)
!
!
!*       2.5    Effect on temperature
!               ---------------------
!
    ptg(jj,1)   = ptg(jj,1) + pdwgi1(jj)*xlmtt*pct(jj)*zwork2(jj)
!
!-------------------------------------------------------------------------------
!
!*       3.     EFFECT OF THE MELTING/FREEZING 
!               ON THE DEEP-SOIL HEAT AND ICE CONTENTS
!               ('DEF' or Force-Restore soil option)
!               --------------------------------------
!
    zwork1(jj) = pd_g(jj,1)/pd_g(jj,2)
!*       3.1  Available Deep ice content
!             --------------------------
!
    zwim(jj) = ( pwgi(jj,2) - zwork1(jj) * pwgi(jj,1) )  / ( 1. - zwork1(jj) )
!
    zwim(jj) = max(0.,zwim(jj))  ! Just in case of round-off errors
!
!*       3.2  Deep liquid water content
!             -------------------------
!
    zwm(jj)  = ( pwg(jj,2) - zwork1(jj) * pwg(jj,1) )  / ( 1. - zwork1(jj) )
!
!*       3.3    Water freezing
!               --------------
!
! Total soil volumetric heat capacity [J/(m3 K)]:
!
    zsoilheatcap(jj) = xcl*xrholw*pwg(jj,2)  +                           &
                     xci*xrholi*pwgi(jj,2) +                           &
                     xsphsoil*xdrywght*(1.0-zwsat_avgz(jj))*(1.0-zwsat_avgz(jj))
!
! Soil thickness which corresponds to T2 (m): 2 times the diurnal
! surface temperature wave penetration depth as T2 is the average
! temperature for this layer:
!
    ztdiurn(jj)   = min(pd_g(jj,2), 4./(zsoilheatcap(jj)*pcg(jj)))
!
! Effective soil ice penetration depth (m):
!
    ziceeff(jj)   = (pwgi(jj,2)/(pwgi(jj,2)+pwg(jj,2)))*pd_g(jj,2)
!
ENDDO
!
IF(hsoilfrz == 'LWT')THEN
!
! as for the surface layer (above)JJ 
! Note also that if the 'DIF'
! soil option is not in force, then the soil parameters are assumed
! to be homogeneous (in the verticalJJ thus we use 1st element of 2nd dimension
! of the 2D-soil parameter arrays).
!
  DO jj=1,ini

       zmatpot(jj)   = min(pmpotsatz(jj,1), xlmtt*(ptg(jj,2)-xtt)/(xg*ptg(jj,2)) )
       zwgmin(jj)    = zwsat_avgz(jj)*( (zmatpot(jj)/pmpotsatz(jj,1))**(-1./pbcoefz(jj,1)) )

       zmatpot(jj)   = pmpotsatz(jj,1)*( (pwg(jj,2)/zwsat_avgz(jj))**(-pbcoefz(jj,1)) )
       ztgmax(jj)    = xlmtt*xtt/(xlmtt - xg* zmatpot(jj))
  ENDDO
ELSE
    zwgmin(:)    = xwgmin
    ztgmax(:)    = xtt
ENDIF
!
! Allow freezing by T2 up to a certain depth so that
! T2 energy can not be used to freeze soil water
! at levels sufficiently deep in the soil.
!
DO jj=1,ini
!
    zdeltat(jj)  = ptg(jj,2) - ztgmax(jj) ! initial temperature depression 
!  
    zwork1(jj) = pd_g(jj,1)/pd_g(jj,2)
    zwork2(jj) = xrholw*(pd_g(jj,2)-pd_g(jj,1))

    zfreezing(jj) = 0.0
    IF (ziceeff(jj) <= ztdiurn(jj)) THEN
!
      zeffic(jj)    = max(zeffic_min, max(0.0,zwm(jj) - xwgmin)/zwsat_avgz(jj))
      zfreezing(jj) = min( max(0.0, zwm(jj) - zwgmin(jj))* zwork2(jj),            &
                     zksoil*zeffic(jj)*max( -zdeltat(jj) , 0. ) )
    ENDIF
!
!
!*       3.4    Ground Ice melt
!               ---------------
!
    zeffic(jj)    = max(zeffic_min, zwim(jj)/(zwsat_avgz(jj)-xwgmin))
    zice_melt(jj) = min( zwim(jj)*zwork2(jj),             &
                  zksoil*zeffic(jj)*max( zdeltat(jj) , 0. ) )
!
!
!*       3.5    Deep-part of deep-soil Ice reservoir evolution
!               ----------------------------------------------
!
    zwit(jj)   = zwim(jj) + (ptstep/ztauice(jj))*(1.0-zpsng(jj))*       &
               ((zfreezing(jj) - zice_melt(jj))/ zwork2(jj))
!
    zwit(jj)   = max( zwit(jj) , 0.             )
    zwit(jj)   = min( zwit(jj) , zwsat_avgz(jj)-xwgmin)
!
!
!*       3.6    Add reservoir evolution from surface freezing (WI2 contains WI1)
!               ----------------------------------------------------------------
!
    zwgi2(jj)  = (1.-zwork1(jj))*zwit(jj) +  zwork1(jj)*zwgi1(jj)
!
    pdwgi2(jj) = zwgi2(jj) - pwgi(jj,2)
!
!
!*       3.7    Effect on temperature
!               ---------------------
!
    ptg(jj,2) = ptg(jj,2) + pdwgi2(jj)*xlmtt*pcg(jj)*xrholw*pd_g(jj,2)
!
ENDDO
!
IF (lhook) CALL dr_hook('ICE_SOILFR',1,zhook_handle)
!-------------------------------------------------------------------------------
END SUBROUTINE ice_soilfr