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(PTSTEP, PTAUICE, PKSFC_IVEG, PLEGI,                 &
                            psoilhcapz, pwsatz, pmpotsatz, pbcoefz,              &
                            ptg, pwgi, pwg, kwg_layer,                           &
                            pdzg, 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   Boone
!!      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_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
!
REAL, INTENT(IN)                   :: ptstep  ! Model time step (s)
!
REAL, DIMENSION(:), INTENT(IN)      :: ptauice, pksfc_iveg, plegi
!                                      PKSFC_IVEG = effect of surface layer insolation on phase changes
!                                                    Giard and Bazile (2000): non-dimensional
!                                      PTAUICE    = soil phase change characteristic time scale (s)
!                                      PLEGI      = ice sublimation (m s-1)
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psoilhcapz, pwsatz
!                                      PSOILHCAPZ = soil heat capacity [J/(m3 K)]
!                                      PWSATZ     = soil porosity (m3/m3)
!
REAL, DIMENSION(:,:), INTENT(IN)    :: pdzg
!                                      PDZG   = Layer thickness (DIF option)
!
REAL, DIMENSION(:,:), INTENT(IN)    :: pmpotsatz, pbcoefz
!                                      PMPOTSATZ  = matric potential at saturation (m)
!                                      PBCOEFZ    = slope of the water retention curve (-)
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: ptg, pwgi, pwg
!                                      PTG    = soil temperature (K)
!                                      PWGI   = soil volumetric ice content (m3/m3)
!                                      PWGI   = soil volumetric liquid water content (m3/m3)
!
INTEGER, DIMENSION(:), INTENT(IN)   :: kwg_layer  
!                                      KWG_LAYER = Number of soil moisture layers (DIF option)
!
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(PTG,1),SIZE(PTG,2)) :: zk, zexcessfc
!
REAL, DIMENSION(SIZE(PTG,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(ptg(:,:),1)
inl = maxval(kwg_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=kwg_layer(jj)
    IF(jl<=idepth)THEN
!
      zwgim = pwgi(jj,jl)
      zwgm  = pwg(jj,jl)
      ztgm  = ptg(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(pmpotsatz(jj,jl),xlmtt*(ptg(jj,jl)-xtt)/(xg*ptg(jj,jl)))
!        
      zwork  = zpsimax/pmpotsatz(jj,jl)
      zlog   = log(zwork)/pbcoefz(jj,jl)
      zwgmax = pwsatz(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 = pwg(jj,jl)/pwsatz(jj,jl)
      zlog  = pbcoefz(jj,jl)*log(zwork)
      zpsi  = pmpotsatz(jj,jl)*exp(-zlog)
!
      zdeltat = ptg(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*pbcoefz(jj,jl)*ztgm*ztgm)
      ENDIF
!
!     *Melt* ice if energy and ice available:
      zphasem  = (ptstep/ptauice(jj))*min(zk(jj,jl)*xci*xrholi*max(0.0,zdeltat),zwgim*xlmtt*xrholw)
!
!     *Freeze* liquid water if energy and water available:
      zphasef  = (ptstep/ptauice(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
      ptg(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)*(ptg(jj,jl)-ztgm)
!
!     Adjust ice and liquid water conents (m3/m3) accordingly :
      pwgi(jj,jl) = zwgim + zphase/(xlmtt*xrholw)     
      pwg(jj,jl) = zwgm  - zphase/(xlmtt*xrholw) 
!
    ENDIF
  ENDDO
ENDDO
!
! 3. Adjust surface soil ice content for sublimation
!    -----------------------------------------------
!
pwgi(:,1) = pwgi(:,1) - plegi(:)*ptstep/pdzg(:,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,  - pwgi(:,1))
pwg(:,1)   = pwg(:,1) - zexcess(:)
pwgi(:,1)   = pwgi(:,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=kwg_layer(jj)
    IF(jl<=idepth)THEN
      zexcess(jj)       = max(0.0, pwgi(jj,jl) - (pwsatz(jj,jl)-xwgmin) )
      pwgi(jj,jl)       = pwgi(jj,jl)   - zexcess(jj)
      zexcessfc(jj,jl)  = zexcessfc(jj,jl) + zexcess(jj)
      IF(jl<idepth)THEN
        pwgi(jj,jl+1)     = pwgi(jj,jl+1) + zexcess(jj)*(pdzg(jj,jl)/pdzg(jj,jl+1))
        zexcessfc(jj,jl+1)= zexcessfc(jj,jl+1) - zexcess(jj)*(pdzg(jj,jl)/pdzg(jj,jl+1))
      ELSE
        pwgi_excess(jj)      = zexcess(jj)*pdzg(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=kwg_layer(jj)  
    IF(jl<=idepth.AND.pwgi(jj,jl)>0.0.AND.pwgi(jj,jl)<1.0e-6)THEN
      pwg(jj,jl) = pwg(jj,jl) + pwgi(jj,jl)
      zexcessfc(jj,jl) = zexcessfc(jj,jl) + pwgi(jj,jl)
      pwgi(jj,jl) = 0.0
    ENDIF
    ptg(jj,jl) = ptg(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