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(IO, KK, PK, PEK, DMK, PTSTEP, PKSFC_IVEG, PDWGI1, PDWGI2 )   
!!     ##########################################################################
!
!!****  *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 'IO%CSOILFRZ' 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_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, 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
!
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) :: DMK
!
REAL, INTENT (IN)                   :: PTSTEP     ! model time step (s)
!
REAL, DIMENSION(:), INTENT(IN)      :: PKSFC_IVEG
!                                      PKSFC_IVEG= non-dimensional vegetation insolation coefficient
!
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(DMK%XCG)) ::   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(DMK%XCG)) ::  ZWSAT_AVGZ
!                               ZWSAT_AVGZ = soil column average porosity (m3 m-3)
!
REAL, DIMENSION(SIZE(DMK%XCG)) :: ZPSNG
!                               ZPSNG = snow fractions corresponding to
!                                       dummy argument PEK%XPSNG(:)
!                                       if PEK%TSNOW%SCHEME = '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 
!                                       PEK%TSNOW%SCHEME = '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         :: INJ, JJ
!
REAL, DIMENSION(SIZE(DMK%XCG))          :: 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
!
inj = SIZE(pek%XTG,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(pek%TSNOW%SCHEME == '3-L' .OR. pek%TSNOW%SCHEME == 'CRO')THEN
  zpsng(:)     = 0.0
ELSE
  zpsng(:)     = pek%XPSNG(:)+kk%XFFG(:)
ENDIF
!
!*       1.    Melting/freezing normalized coefficient
!               ---------------------------------------
!
zksoil      = (0.5 * sqrt(xcondi*xci*xrholi*xday/xpi))/xlmtt
!
ztauice(:) = max(ptstep,pk%XTAUICE(:))
!
DO jj=1,inj
!-------------------------------------------------------------------------------
!*       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) = kk%XWSAT(jj,1)
!
!               Influence of vegetation insolation on surface:
!
  zksfc_frz(jj) = zksoil * pksfc_iveg(jj)
!
ENDDO
!*       2.2    Water freezing
!               --------------
!
IF(io%CSOILFRZ == '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,inj
    zmatpot(jj)   = min(kk%XMPOTSAT(jj,1), xlmtt*(pek%XTG(jj,1)-xtt)/(xg*pek%XTG(jj,1)) )
    zwgmin(jj)    = zwsat_avgz(jj)*( (zmatpot(jj)/kk%XMPOTSAT(jj,1))**(-1./kk%XBCOEF(jj,1)) )

    zmatpot(jj)   = kk%XMPOTSAT(jj,1)*( (pek%XWG(jj,1)/zwsat_avgz(jj))**(-kk%XBCOEF(jj,1)) )
    ztgmax(jj)    = xlmtt*xtt/(xlmtt - xg* zmatpot(jj))
  ENDDO
ELSE
  zwgmin(:)    = xwgmin
  ztgmax(:)    = xtt
ENDIF
!
DO jj=1,inj
! 
  zdeltat(jj)  = pek%XTG(jj,1) - ztgmax(jj) ! initial temperature depression
!
  zwork2(jj) = xrholw*pk%XDG(jj,1)
  zeffic(jj)    = max(zeffic_min,(pek%XWG(jj,1)-xwgmin)/zwsat_avgz(jj))
  zfreezing(jj) = min( max(0.0,pek%XWG(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,pek%XWGI(jj,1)/(zwsat_avgz(jj)-xwgmin))
  zice_melt(jj) = min( pek%XWGI(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) = pek%XWGI(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) - pek%XWGI(jj,1)
!
!
!*       2.5    Effect on temperature
!               ---------------------
!
  pek%XTG(jj,1)   = pek%XTG(jj,1) + pdwgi1(jj)*xlmtt*dmk%XCT(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) = pk%XDG(jj,1)/pk%XDG(jj,2)
!*       3.1  Available Deep ice content
!             --------------------------
!
  zwim(jj) = ( pek%XWGI(jj,2) - zwork1(jj) * pek%XWGI(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)  = ( pek%XWG(jj,2) - zwork1(jj) * pek%XWG(jj,1) )  / ( 1. - zwork1(jj) )
!
!*       3.3    Water freezing
!               --------------
!
! Total soil volumetric heat capacity [J/(m3 K)]:
!
  zsoilheatcap(jj) = xcl*xrholw*pek%XWG(jj,2)  +                           &
                     xci*xrholi*pek%XWGI(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(pk%XDG(jj,2), 4./(zsoilheatcap(jj)*dmk%XCG(jj)))
!
! Effective soil ice penetration depth (m):
!
  ziceeff(jj)   = (pek%XWGI(jj,2)/(pek%XWGI(jj,2)+pek%XWG(jj,2)))*pk%XDG(jj,2)
!
ENDDO
!
IF(io%CSOILFRZ == '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,inj

    zmatpot(jj)   = min(kk%XMPOTSAT(jj,1), xlmtt*(pek%XTG(jj,2)-xtt)/(xg*pek%XTG(jj,2)) )
    zwgmin(jj)    = zwsat_avgz(jj)*( (zmatpot(jj)/kk%XMPOTSAT(jj,1))**(-1./kk%XBCOEF(jj,1)) )

    zmatpot(jj)   = kk%XMPOTSAT(jj,1)*( (pek%XWG(jj,2)/zwsat_avgz(jj))**(-kk%XBCOEF(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,inj
!
  zdeltat(jj)  = pek%XTG(jj,2) - ztgmax(jj) ! initial temperature depression 
!  
  zwork1(jj) = pk%XDG(jj,1)/pk%XDG(jj,2)
  zwork2(jj) = xrholw*(pk%XDG(jj,2)-pk%XDG(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) - pek%XWGI(jj,2)
!
!
!*       3.7    Effect on temperature
!               ---------------------
!
  pek%XTG(jj,2) = pek%XTG(jj,2) + pdwgi2(jj)*xlmtt*dmk%XCG(jj)*xrholw*pk%XDG(jj,2)
!
ENDDO
!
IF (lhook) CALL dr_hook('ICE_SOILFR',1,zhook_handle)
!-------------------------------------------------------------------------------
END SUBROUTINE ice_soilfr