soilstress.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 soilstress( HISBA, PF2,                                   &
                  prootfrac, pwsat, pwfc, pwwilt,                    &
                  pwg, pwgi, kwg_layer, pf2wght, pf5                 )  
!     ####################################################################
!
!!****  *SOILSTRESS*  
!!
!!    PURPOSE
!!    -------
!
!     Calculates the vegetation stress due to soil water
!         
!     
!!**  METHOD
!!    ------
!
!     Calculates the F2 coefficient.
!     
!
!!    EXTERNAL
!!    --------
!!
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!    none
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Noilhan and Planton (1989)
!!    Belair (1995)
!!      
!!    AUTHOR
!!    ------
!!
!!      S. Belair           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original      13/03/95 
!!     (P.Jabouille)  13/11/96    mininum value for ZF1
!!     (V. Masson)    28/08/98    add PF2 for Calvet (1998) CO2 computations
!!     (B. Decharme)     07/15    Suppress numerical adjustement for PF2 
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!
USE modd_isba_par,  ONLY : xdenom_min
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
 CHARACTER(LEN=*),     INTENT(IN)   :: hisba   ! type of soil (Force-Restore OR Diffusion)
!                                             ! '2-L'
!                                             ! '3-L'
!                                             ! 'DIF'   ISBA-DF
!
REAL, DIMENSION(:,:), INTENT(IN) :: prootfrac, pwsat, pwfc, pwwilt,       &
                                      pwg, pwgi  
!                                     PROOTFRAC = cumulative root fraction (-)
!                                     PWFC      = field capacity profile (m3/m3)
!                                     PWWILT    = wilting point profile (m3/m3)
!                                     PWSAT     = porosity profile (m3/m3)
!                                     PWG       = soil liquid volumetric water content (m3/m3)
!                                     PWGI      = soil frozen volumetric water content (m3/m3)
!
INTEGER, DIMENSION(:), INTENT(IN) :: kwg_layer  
!                                    KWG_LAYER = Number of soil moisture layers (DIF option)
!
REAL, DIMENSION(:), INTENT(OUT)  :: pf2      ! water stress coefficient
!
REAL, DIMENSION(:), INTENT(OUT)  :: pf5      ! water stress coefficient for Hv (based on F2):
!                                            ! Verify that Etv=>0 as F2=>0
!
REAL, DIMENSION(:,:), INTENT(OUT):: pf2wght  ! water stress coefficient profile (ISBA-DF)
!
!
!*      0.2    declarations of local variables
!
!
REAL, DIMENSION(SIZE(PWFC,1)) ::  zwfc_avgz, zwsat_avgz, zwwilt_avgz
!                                  ZWFC_AVGZ   = field capacity averaged over entire soil column
!                                  ZWSAT_AVGZ  = porosity averaged over entire soil column
!                                  ZWWILT_AVGZ = wilting point averaged over entire soil column
!
! ISBA-DF:
!
REAL, DIMENSION(SIZE(PWG,1)) :: zwsat, zwfc, zwwilt,zf2wght
!                               ZWSAT     = ice-adjusted porosity profile (m3/m3)
!                               ZWFC      = ice-adjusted field capacity profile (m3/m3)
!                               ZWWILT    = ice-adjusted wilting point profile (m3/m3)
!
REAL    :: zrootfracn
!          ZROOTFRACN = Normalized root fraction weights
!
INTEGER :: ini, inl, jj, jl, idepth
!
!
!*      0.3    declarations of local parameters:
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
!*       0.     Initialization of variables:
!               ----------------------------
!
IF (lhook) CALL dr_hook('SOILSTRESS',0,zhook_handle)
!
ini=SIZE(pwg,1)
inl=maxval(kwg_layer(:))
!
pf2(:)      = 0.0
pf2wght(:,:)    = 0.0
!
zwfc_avgz(:)    = 0.
zwsat_avgz(:)   = 0.
zwwilt_avgz(:)  = 0.
!
zf2wght(:)      = 0.
!
!-------------------------------------------------------------------------------
!
!*       2.     THE 'PF2' FACTOR
!               ----------------
!               This factor takes into account the effect
!               of the water stress on the surface
!               resistance
!
! - For humid soils (> WFC), this factor does not
!   increase the stomatal resistance                                   
! - The stomatal resistance should be large
!   when the soil is very dry (< WILT)
!
IF(hisba =='DIF')THEN      
!
! If using the diffusion option, then calculate transpiration weights
! and the mean root-zone soil water stress factor F2:
!
!---------------------------------------------------------
! First layer
!---------------------------------------------------------
!
! Due to the presence of ice, modify soil parameters:
  zwsat(:) = pwsat(:,1) - pwgi(:,1)
  zwfc(:) = pwfc(:,1) * zwsat(:)/pwsat(:,1)
  zwwilt(:) = pwwilt(:,1) * zwsat(:)/pwsat(:,1)
!
! Calculate the soil water stress factor for each layer:
  pf2wght(:,1) = max(xdenom_min,min(1.0,(pwg(:,1)-zwwilt(:))/(zwfc(:)-zwwilt(:))))
  zf2wght(:)   = max(       0.0,min(1.0,(pwg(:,1)-zwwilt(:))/(zwfc(:)-zwwilt(:))))
!
! Normalize the transpiration weights by root fraction:
  pf2wght(:,1) = prootfrac(:,1)*pf2wght(:,1)
  zf2wght(:)   = prootfrac(:,1)*zf2wght(:)
!
! Net soil water stress for entire root zone:
  pf2(:)       = zf2wght(:)
!
!---------------------------------------------------------
! Other layers
!---------------------------------------------------------
!
  DO jl=2,inl
     DO jj=1,ini
!
      idepth=kwg_layer(jj)
      IF(jl<=idepth)THEN
!
!       Due to the presence of ice, modify soil parameters:
        zwsat(jj) = pwsat(jj,jl) - pwgi(jj,jl)
        zwfc(jj) = pwfc(jj,jl) * zwsat(jj)/pwsat(jj,jl)
        zwwilt(jj) = pwwilt(jj,jl) * zwsat(jj)/pwsat(jj,jl)
!
!       Calculate normalized root fraction weights:
        zrootfracn = prootfrac(jj,jl) - prootfrac(jj,jl-1)
!
!       Calculate the soil water stress factor for each layer:
!!        PF2WGHT(JJ,JL) = MAX(0.0,MIN(1.0,(PWG(JJ,JL)-ZWWILT(JJ))/(ZWFC(JJ)-ZWWILT(JJ))))

        pf2wght(jj,jl) = max(xdenom_min,min(1.0,(pwg(jj,jl)-zwwilt(jj))/(zwfc(jj)-zwwilt(jj))))
        zf2wght(jj)    = max(       0.0,min(1.0,(pwg(jj,jl)-zwwilt(jj))/(zwfc(jj)-zwwilt(jj))))
!
!       Normalize the transpiration weights by root fraction:                                                
        pf2wght(jj,jl) = zrootfracn*pf2wght(jj,jl)
        zf2wght(jj)    = zrootfracn*zf2wght(jj)
!
!       Net soil water stress for entire root zone:
        pf2(jj)        = pf2(jj) + zf2wght(jj)
!        
      ENDIF
!
     ENDDO
  ENDDO
!
ELSE
!
! When using the Force-Restore (FR) soil option, the soil properties
! are assumed to be homogeneous in the vertical. Therefore, simply
! use the properties for the uppermost soil layer when defining
! soil properties for local computations.
!
! Due to the presence of ice, modify soil parameters:
!
   zwsat_avgz(:)  = pwsat(:,1) - pwgi(:,2)
   zwfc_avgz(:)   = pwfc(:,1)*zwsat_avgz(:)/pwsat(:,1)
   zwwilt_avgz(:) = pwwilt(:,1)*zwsat_avgz(:)/pwsat(:,1)
!
! Compute the water stress factor:
!
   pf2(:) = (pwg(:,2)-zwwilt_avgz(:))/(zwfc_avgz(:)-zwwilt_avgz(:))
   pf2(:) = max(0.0,min(1.0, pf2(:)))
!
!
ENDIF
!
! Function to cause Etv to approach 0 as F2 goes to 0:
!
pf5(:) = pf2(:)
!
IF (lhook) CALL dr_hook('SOILSTRESS',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE soilstress