ice_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 ice_soildif(KK, PK, PEK, PTSTEP, PKSFC_IVEG, PLEGI, PSOILHCAPZ, PWGI_EXCESS)  
!     ##########################################################################
!
!!****  *ICE_SOILDIF*  
!
!!    PURPOSE
!!    -------
!     This subroutine calculates soil water phase changes using the
!     available/excess energy approach. Soil temperature and volumetric
!     ice content are adjusted due to phase changes. See the references
!     Boone et al., 39, JAM, 2000 and Giard and Bazile, 128, MWR, 2000.
!     NOTE that more recently a modification was made: freeze/thaw follows
!     a relationship between liquid water and temperature derriving from
!     the Clausius Clapeyron Eq. This results in little to no freezing for
!     sufficiently dry but cold (below freezing) soils. Scatter about this
!     curve results due to 'phase change efficiencies' and the surface insolation
!     coefficient.
!     
!!**  METHOD
!!    ------
!
!     Direct calculation
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Boone (2000)
!!    Boone et al., (2000)
!!      
!!    AUTHOR
!!    ------
!!      A. Boone          * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    28/02/00   Boonei
!!      Modified    24/11/09   Boone
!!                             Limit energy available for phase change by
!                              local amount owing to diffusion. Has almost
!                              no impact except under rare circumstances
!                              (avoids rare but possible oscillatory behavior)
!                              Also, add minimum (numerical) melt/freeze efficieny to prevent
!                              prolonged periods of small ice amounts approaching zero.
!
!!      Modified    01/06/11   Boone
!                              Use apparent heat capacity linearization for freezing
!                              (when temperature depenence on ice change is direct)
!                              Do away with efficiency coefficients as they acted to provide
!                              numerical stability: not needed as apparent heat capacity increases
!                              the effective heat capacity thus increasing stability.
!                              NOTE: for now considers Brooks & Corey type water retention.
!
!      Modified    08/2011     Decharme
!                              Optimization
!      Modified    10/2013     Boone
!                              Slight edit to phase computation to improve enthalpy conservation
!                              
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_isba_n, ONLY : isba_k_t, isba_p_t, isba_pe_t
!
USE modd_csts,     ONLY : xlmtt, xtt, xg, xci, xrholi, xrholw
USE modd_isba_par, ONLY : xwgmin
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
TYPE(isba_k_t), INTENT(INOUT) :: KK
TYPE(isba_p_t), INTENT(INOUT) :: PK
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
!
REAL, INTENT(IN)                   :: PTSTEP  ! Model time step (s)
!
REAL, DIMENSION(:), INTENT(IN)      :: PKSFC_IVEG, PLEGI
!                                      PKSFC_IVEG = effect of surface layer insolation on phase changes
!                                                    Giard and Bazile (2000): non-dimensional
!                                      PLEGI      = ice sublimation (m s-1)
!
REAL, DIMENSION(:,:), INTENT(IN)    :: PSOILHCAPZ
!                                      PSOILHCAPZ = soil heat capacity [J/(m3 K)]
!
REAL, DIMENSION(:), INTENT(OUT)     :: PWGI_EXCESS
!                                      PWGI_EXCESS = Soil ice excess water content
!
!*      0.2    declarations of local variables
!
INTEGER                             :: JJ, JL   ! loop control
!
INTEGER                             :: INI    ! Number of point
INTEGER                             :: INL    ! Number of explicit soil layers
INTEGER                             :: IDEPTH ! Total moisture soil depth
!
REAL, DIMENSION(SIZE(PEK%XTG,1),SIZE(PEK%XTG,2)) :: ZK, ZEXCESSFC
!
REAL, DIMENSION(SIZE(PEK%XTG,1))             :: ZEXCESS
!
REAL                                     :: ZWGMAX, ZPSIMAX, ZPSI, ZDELTAT,  &
                                            ZPHASE, ZTGM, ZWGM, ZWGIM, ZLOG, &
                                            ZEFFIC, ZPHASEM, ZPHASEF, ZWORK, &
                                            ZAPPHEATCAP
!                                            
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
! Initialization:
! ---------------
!
IF (lhook) CALL dr_hook('ICE_SOILDIF',0,zhook_handle)
!
ini = SIZE(pek%XTG,1)
inl = maxval(pk%NWG_LAYER)
!
zexcessfc(:,:)=0.0
zexcess(:  )=0.0
pwgi_excess(:  )=0.0
!
!-------------------------------------------------------------------------------
!
! 1. Surface layer vegetation insulation coefficient (-)
!    ---------------------------------------------------
!
zk(:,:) = 1.0
zk(:,1) = pksfc_iveg(:)
!
! 2. Soil ice evolution computation:
!    -------------------------------
!
DO jl=1,inl
  DO jj=1,ini                 
    idepth=pk%NWG_LAYER(jj)
    IF(jl<=idepth)THEN
!
      zwgim = pek%XWGI(jj,jl)
      zwgm  = pek%XWG (jj,jl)
      ztgm  = pek%XTG (jj,jl)

!     The maximum liquid water content as
!     as function of temperature (sub-freezing)
!     based on Gibbs free energy (Fuchs et al., 1978):
!
      zpsimax = min(kk%XMPOTSAT(jj,jl),xlmtt*(pek%XTG(jj,jl)-xtt)/(xg*pek%XTG(jj,jl)))
!        
      zwork  = zpsimax/kk%XMPOTSAT(jj,jl)
      zlog   = log(zwork)/kk%XBCOEF(jj,jl)
      zwgmax = kk%XWSAT(jj,jl)*exp(-zlog)
!
!     Calculate maximum temperature for ice based on Gibbs free energy: first
!     compute soil water potential using Brook and Corey (1966) model:
!     psi=mpotsat*(w/wsat)**(-bcoef)
!
      zwork = pek%XWG(jj,jl)/kk%XWSAT(jj,jl)
      zlog  = kk%XBCOEF(jj,jl)*log(zwork)
      zpsi  = kk%XMPOTSAT(jj,jl)*exp(-zlog)
!
      zdeltat = pek%XTG(jj,jl) - xlmtt*xtt/(xlmtt-xg*zpsi)
!
!     Compute apparent heat capacity. This is considered
!     only when there is available energy (cold) and liquid water
!     available...freezing front.
!     This also has the secondary effect of increasing numerical stability
!     during freezing (as there is a strong temperature dependence) by
!     i) potentially significantly increasing the "apparent" heat capacity and
!     ii) this part is also treated implicitly herein.
!
      zwork = (xci*xrholi/(xlmtt*xrholw))*zk(jj,jl)*max(0.0,-zdeltat)
!        
      zappheatcap=0.0
      IF(zdeltat<0.0.AND.zwgm>=zwgmax.AND.zwork>=max(0.0,zwgm-zwgmax))THEN
        zappheatcap = -(xtt*xrholw*xlmtt*xlmtt/xg)*zwgmax/(zpsimax*kk%XBCOEF(jj,jl)*ztgm*ztgm)
      ENDIF
!
!     *Melt* ice if energy and ice available:
      zphasem  = (ptstep/pk%XTAUICE(jj))*min(zk(jj,jl)*xci*xrholi*max(0.0,zdeltat),zwgim*xlmtt*xrholw)
!
!     *Freeze* liquid water if energy and water available:
      zphasef  = (ptstep/pk%XTAUICE(jj))*min(zk(jj,jl)*xci*xrholi*max(0.0,-zdeltat),max(0.0,zwgm-zwgmax)*xlmtt*xrholw)
!
!     Update heat content if melting or freezing
      pek%XTG(jj,jl) = ztgm + (zphasef - zphasem)/(psoilhcapz(jj,jl)+zappheatcap)
!
!     Get estimate of actual total phase change (J/m3) for equivalent soil water changes:
      zphase = (psoilhcapz(jj,jl)+zappheatcap)*(pek%XTG(jj,jl)-ztgm)
!
!     Adjust ice and liquid water conents (m3/m3) accordingly :
      pek%XWGI(jj,jl) = zwgim + zphase/(xlmtt*xrholw)     
      pek%XWG(jj,jl) = zwgm  - zphase/(xlmtt*xrholw) 
!
    ENDIF
  ENDDO
ENDDO
!
! 3. Adjust surface soil ice content for sublimation
!    -----------------------------------------------
!
pek%XWGI(:,1) = pek%XWGI(:,1) - plegi(:)*ptstep/pk%XDZG(:,1)
!
! The remaining code in this block are merely constraints to ensure a highly
! accurate water budget: most of the time this code will not have any
! effect on the soil water profile.
! If sublimation causes all of the remaining ice to be extracted, remove
! some of the liquid (a correction): NOTE that latent heating already accounted
! for in sublimation term, so no need to alter soil temperature.
!
zexcess(:)  = max(0.0,  - pek%XWGI(:,1))
pek%XWG(:,1)   = pek%XWG  (:,1) - zexcess(:)
pek%XWGI(:,1)  = pek%XWGI (:,1) + zexcess(:)
zexcessfc(:,1) = zexcessfc(:,1) - zexcess(:)
!
! 4. Prevent some possible problems
!    ------------------------------
!
! If sublimation is negative (condensation), make sure ice does not
! exceed maximum possible. If it does, then put excess ice into layer below:
! This correction should rarely if ever cause any ice accumulation in the
! sub-surface layer: this is especially true of deeper layers but it is
! accounted for none-the-less.
!
DO jl=1,inl
  DO jj=1,ini
    idepth=pk%NWG_LAYER(jj)
    IF(jl<=idepth)THEN
      zexcess(jj)       = max(0.0, pek%XWGI(jj,jl) - (kk%XWSAT(jj,jl)-xwgmin) )
      pek%XWGI(jj,jl)   = pek%XWGI(jj,jl)  - zexcess(jj)
      zexcessfc(jj,jl)  = zexcessfc(jj,jl) + zexcess(jj)
      IF(jl<idepth)THEN
        pek%XWGI(jj,jl+1)  = pek%XWGI(jj,jl+1)  + zexcess(jj)*(pk%XDZG(jj,jl)/pk%XDZG(jj,jl+1))
        zexcessfc(jj,jl+1) = zexcessfc(jj,jl+1) - zexcess(jj)*(pk%XDZG(jj,jl)/pk%XDZG(jj,jl+1))
      ELSE
        pwgi_excess(jj)    = zexcess(jj)*pk%XDZG(jj,idepth)*xrholw/ptstep
      ENDIF
    ENDIF
  ENDDO
ENDDO
!   
! Prevent keeping track of very small numbers for ice content: (melt it)
! and conserve energy:
!
DO jl=1,inl
  DO jj=1,ini 
    idepth=pk%NWG_LAYER(jj)  
    IF(jl<=idepth.AND.pek%XWGI(jj,jl)>0.0.AND.pek%XWGI(jj,jl)<1.0e-6)THEN
      pek%XWG   (jj,jl)  = pek%XWG(jj,jl) + pek%XWGI(jj,jl)
      zexcessfc(jj,jl) = zexcessfc(jj,jl) + pek%XWGI(jj,jl)
      pek%XWGI(jj,jl) = 0.0
    ENDIF
    pek%XTG(jj,jl) = pek%XTG(jj,jl) - zexcessfc(jj,jl)*xlmtt*xrholw/psoilhcapz(jj,jl)           
  ENDDO
ENDDO
!
!
IF (lhook) CALL dr_hook('ICE_SOILDIF',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE ice_soildif