cotwores.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 cotwores(PTSTEP, IO, OSHADE, PK, PEK, PDMAX, PPOI, PCSP, &
                    PTG, PF2, PSW_RAD, PQA, PQSAT, PPSNV, PDELTA, PRHOA, &
                    PZENITH, PFFV, PIACAN_SUNLIT, PIACAN_SHADE, PFRAC_SUN, &
                    PIACAN, PABC, PRS, PGPP, PRESP_LEAF     ) 
!   #########################################################################
!
!!****  *COTWORES*  
!!
!!    PURPOSE
!!    -------
!!
!!    Calculates net assimilation of CO2 and leaf conductance.
!!              
!!**  METHOD
!!    ------
!!    Calvet et al. 1998 Forr. Agri. Met. [from model of Jacobs(1994)]
!!
!!    EXTERNAL
!!    --------
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!      
!!    USE MODD_CST
!!    USE MODD_CO2V_PAR
!!    USE MODI_COTWO
!!    USE MODI_CCETR
!!    USE MODE_THERMOS
!!
!!    REFERENCE
!!    ---------
!!
!!    Calvet et al. 1998 Forr. Agri. Met. 
!!      
!!    AUTHOR
!!    ------
!!
!!      A. Boone           * Meteo-France *
!!      (following Belair)
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    27/10/97 
!!      V. Masson and V. Rivailland 12/2003 modificatino of ISBA routines order
!!      L. Jarlan   27/10/04 : add of T2 to calculate soil respiration and use
!!                              of CRESPSL key to manage the calculation of soil
!!                              respiration
!!                             PEK%XAN et PPST in kgCO2 m-2 s-1 to be fully
!!                              compatible with vegetation growth module (lailoss.f90)
!!      P Le Moigne 09/2005 AGS modifs of L. Jarlan
!!      S. Lafont      03/09 : change units of EPSO GAMM ANDAY
!!      A.L. Gibelin   06/09 : suppress evolution of [CO2] in canopy
!!      A.L. Gibelin   06/09 : move calculations of some CO2 fluxes
!!      A.L. Gibelin   06/09 : add RESP_LEAF
!!      A.L. Gibelin   07/09 : ensure coherence between cotwores and cotworestress
!!      A.L. Gibelin   07/09 : Suppress PPST and PPSTF as outputs, and diagnose GPP
!!        S. Lafont    03/11 : Correct a bug fopr grassland below wilting point
!!      D. Carrer      04/11 : new radiative transfert 
!!      A. Boone       11/11 : add rsmax to MODD_VEG_PAR
!!      B. Decharme    05/12 : Bug : flood fraction in COTWORES
!!                                   Optimization
!!      R. Alkama      04/12 : add 6 new tree vegtype (9 instead 3)
!!      C. Delire      01/14 : vertical profile of dark respiration for tropical forest 
!!                             (GTROP)   with Carrer radiative transfer (IO%LTR_ML = T)               
!!Seferian & Delire  06/2015 : generalization of (i) linear water-stress reponse
!                              and (ii) exponential decrease of autothrophic respiration to all woody PFTs
!!      B. Decharme    07/15 : Suppress some numerical adjustement for F2 
!!
!-------------------------------------------------------------------------------
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_p_t, isba_pe_t
!
USE modd_csts,           ONLY : xmd, xtt, xlvtt
USE modd_isba_par,       ONLY : xrs_max, xdenom_min
USE modd_co2v_par,       ONLY : xparcf, xmco2,       &
                                xdmaxx, xdmaxn, xaw, xbw, xasw                              
USE modd_data_cover_par, ONLY : nvt_tebd, nvt_trbe, nvt_bone,   &
                                nvt_trbd, nvt_tebe, nvt_tene,   &
                                nvt_bobd, nvt_bond, nvt_shrb
USE modd_surf_par,       ONLY : xundef
!
USE modi_ccetr
USE modi_cotwo
!
!*       0.     DECLARATIONS
!               ------------
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
REAL,                INTENT(IN)  :: PTSTEP      ! time step
!
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_p_t), INTENT(INOUT) :: PK
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
!
LOGICAL, DIMENSION(:),INTENT(IN) :: OSHADE
!
REAL, DIMENSION(:),  INTENT(IN)  :: PPOI     ! Gaussian weights (as above)
!
REAL,DIMENSION(:),   INTENT(IN)  :: PDMAX  
!                                    PDMAX     = maximum saturation deficit of 
!                                                  atmosphere tolerate by vegetation
!
REAL,DIMENSION(:),   INTENT(IN)  :: PCSP, PTG, PF2, PSW_RAD 
!                                    PCSP  = atmospheric concentration of CO2
!                                    PTG   = updated leaf temperature
!                                    PF2   = normalized soil water stress factor
!                                    PSW_RAD = incident solar radiation
!
REAL,DIMENSION(:),   INTENT(IN)  :: PQA, PQSAT, PPSNV, PDELTA, PRHOA
!                                    PQA   = atmospheric mixing ratio
!                                    PQSAT = surface saturation mixing ratio
!                                    PPSNV = snow cover fraction
!                                    PDELTA= fraction of the foliage covered
!                                        by intercepted water
!                                    PEK%XRESA = air density
!
REAL,DIMENSION(:),    INTENT(IN)  :: PZENITH
!                                    PZENITH = solar zenith angle needed 
!                                    for computation of diffusuion of solar
!                                    radiation: for CO2 model.
!
REAL, DIMENSION(:,:), INTENT(IN)    :: PIACAN_SUNLIT, PIACAN_SHADE, PFRAC_SUN
!
REAL, DIMENSION(:), INTENT(IN)      :: PFFV ! Floodplain fraction over vegetation
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: PIACAN ! PAR in the canopy at different gauss level
!
REAL,DIMENSION(:),  INTENT(INOUT) :: PABC, PRS, PGPP
!                                    PABC  = Carrer radiative transfer: normalized heigh of considered layer (bottom=0, top=1)
!                                            Calvet radiative transfer: abcissa of the 3-points Gaussian quadrature 
!                                                (Goudriaan, Agric&For.Meteor, 38,1986)    
!                                    PRS   = stomatal resistance
!                                    PGPP  = Gross Primary Production (kg_CO2/kg_air * m/s)
!
REAL,DIMENSION(:),    INTENT(OUT) :: PRESP_LEAF
!                                    PRESP_LEAF = dark respiration over canopy
!
!*      0.2    declarations of local variables
!
REAL, PARAMETER                :: ZRS_MIN     = 1.e-4  ! minimum canopy resistance (s m-1)
!
INTEGER                     :: JINT, JJ ! index for loops
!
REAL, DIMENSION(SIZE(PEK%XLAI,1)) :: ZANF ! ZANF  = total assimilation over canopy
REAL, DIMENSION(SIZE(PEK%XLAI,1)) :: ZCONVE1, ZTSPC, ZIA
!                                 ZTSPC = temperature conversion (K to C) 
!                                 ZIA   = absorbed PAR
REAL, DIMENSION(SIZE(PEK%XLAI,1)) :: ZLAI, ZGMEST, ZFZERO, ZDMAX
!                                 ZLAI = LAI 
!                                 ZFZERO  = ideal value of F, no photorespiration or 
!                                            saturation deficit
!                                 ZDMAX   = maximum saturation deficit of atmosphere
!                                           tolerate by vegetation
!
REAL, DIMENSION(SIZE(PEK%XLAI,1)) :: ZGAMMT, ZDSP, ZANMAX
!                                 ZGAMMT  = compensation point 
!                                 ZDSP    = saturation deficit of atmosphere 
!                                           verses the leaf surface (with correction)
!
REAL, DIMENSION(SIZE(PEK%XLAI,1)) :: ZXMUS, ZTAN, ZTGS, ZXIA, ZAN0, ZGS0, ZXTGS, ZRDK,ZLAITOP,ZTRDK,ZZLAI  
!                                           ZXMUS = cosine of solar zenith angle
!                                           ZTAN  = canopy integrated net assimilation 
!                                           ZTGS  = canopy integrated  leaf conductance
!                                           ZXIA  = incident radiation after diffusion
!                                           ZAN0  = net asimilation at each interval
!                                                   in the canopy
!                                           ZGS0  = leaf conductance at each interval
!                                                   in the canopy        
!                                           ZXTGS = total canopy conductance
!                                           ZRDK  = dark respiration
!                                           ZLAITOP = LAI (thickness of canopy) above considered layer 
!                                           ZTRDK = canopy integrated dark respiration
!                                           ZZLAI = LAI, used for dark respiration profile
!
REAL, DIMENSION(SIZE(PEK%XLAI,1)) :: ZAN0_,ZGS0_,ZRDK_ ! parameters for shaded leaves
!
REAL, DIMENSION(SIZE(PEK%XLAI,1)) :: ZEPSO
!                                           ZEPSO conversion of PEPSO in kgCO2/kgair m/s
!
REAL, DIMENSION(SIZE(PEK%XLAI,1)) :: ZDMAXSTAR, ZFZEROSTAR, ZFZERON, ZGMESTN  
!                                 ZDMAXSTAR  = maximum saturation deficit of atmosphere
!                                              tolerate by vegetation without soil water stress
!                                 ZFZEROSTAR = initial optimal ratio Ci/Cs for woody vegetation
!                                 ZFZERON    = minimum value for "fzero" in defensive woody strategy
!                                 ZGMESTN    = gmest value at zf2=zf2i in offensive woody strategy
!
!
REAL :: ZABC, ZWEIGHT
!                                           ZABC    = abscissa needed for integration
!                                                     of net assimilation and stomatal 
!                                                     conductance over canopy depth 
!                                                     (working scalar)
!
REAL, DIMENSION(SIZE(PEK%XLAI,1))    :: ZWORK !Work array
!
LOGICAL, DIMENSION(SIZE(PEK%XLAI,1)) :: GHERB, GWOOD, GF2_INF_F2I, GTROP
!
INTEGER, DIMENSION(1)          :: IDMAX
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
! STOMATAL RESISTANCE: ENTRY VARIABLES TO CO2 ROUTINE:
!   CS        = CO2 concentration (kgCO2 kgair-1) cs
!   DSP       = specific humidity deficit (kgH2O kgair-1) ds
!   TSM       = surface temperature (C) ts
!   RG        = global radiation (W m-2) rg
!
! initialisation: convert from ppm to mg/m-3
!
IF (lhook) CALL dr_hook('COTWORES',0,zhook_handle)
!
zconve1(:) = xmco2*prhoa/xmd
!
! initialisation: convert from K to C
!
ztspc(:)  = ptg(:) - xtt               
!
zlai(:)   = pek%XLAI(:)
zgmest(:) = pek%XGMES(:)
zfzero(:) = pk%XFZERO(:)
!
!GTROP = linear stress in case of tropical evergreen forest 
!        (with fixed f0=0.74 for Carrer rad. transf., f0=0.7 with Calvet rad. transf.)
gtrop(:) = (pk%XVEGTYPE_PATCH(:,nvt_trbe) > 0.8) 
!
gherb(:) = (pk%XVEGTYPE_PATCH(:,nvt_tebd) + pk%XVEGTYPE_PATCH(:,nvt_trbe) + pk%XVEGTYPE_PATCH(:,nvt_bone)   &
           +pk%XVEGTYPE_PATCH(:,nvt_trbd) + pk%XVEGTYPE_PATCH(:,nvt_tebe) + pk%XVEGTYPE_PATCH(:,nvt_tene)   & 
           +pk%XVEGTYPE_PATCH(:,nvt_bobd) + pk%XVEGTYPE_PATCH(:,nvt_bond) + pk%XVEGTYPE_PATCH(:,nvt_shrb)<0.5)
gwood(:) = (.NOT.gherb(:))
!
!
WHERE (pek%XLAI(:)==xundef) zlai(:)=0.0
!
!    See (Varlet-Granchet C., M. Chartier, G. Gosse,  and R. Bonhomme, 1981: 
!    Rayonnement utilise pour la photosynthese des vegetaux en
!    conditions naturelles: caracterisation et variations. 
!    Oecol. Plant. 2(16), 189-202.)
!
!-------------------------------------
! Add soil moisture stress effect to leaf conductance:
! OFFENSIVE and DEFENSIVE water stress response 
!  
zdmax(:)  = pdmax(:)
!
gf2_inf_f2i(:) = (pf2(:)<pek%XF2I(:))
!
! -HERBACEOUS-
!
WHERE(gherb(:).AND.pek%LSTRESS(:))
  zdmax(:) = xdmaxn
ENDWHERE
WHERE(gherb(:).AND..NOT.pek%LSTRESS(:))
  zdmax(:) = xdmaxx
ENDWHERE
!
! PAH and PBH are original coefficients of Calvet 2000
WHERE(gherb(:).AND.(.NOT.gf2_inf_f2i(:)))
  zdmaxstar(:) = exp((log(zgmest(:)*1000.)-pk%XAH(:))/pk%XBH(:))/1000.
  zdmax(:) = zdmaxstar(:) - (zdmaxstar(:)-zdmax(:))*(1.-pf2(:))/(1.-pek%XF2I(:))
ENDWHERE
!
WHERE(gherb(:))
  zgmest(:) = exp(pk%XAH(:)+pk%XBH(:)*log(zdmax(:)*1000.))/1000.
ENDWHERE
!
WHERE (gherb(:).AND.gf2_inf_f2i(:).AND.pek%LSTRESS(:))
  zgmest(:) = zgmest(:) * pf2(:)/pek%XF2I(:)
ENDWHERE
WHERE(gherb(:).AND.gf2_inf_f2i(:).AND.(.NOT.pek%LSTRESS(:)))
  zdmax(:) = zdmax(:) * pf2(:)/pek%XF2I(:)
ENDWHERE
!
! to limit photosynthesis under wilting point
WHERE (gherb(:).AND.(.NOT.pek%LSTRESS(:)).AND.zdmax(:)<=xdmaxn)
  zdmax(:)  = xdmaxn
  zgmest(:) = (exp(pk%XAH(:)+pk%XBH(:)*log(xdmaxn*1000.))/1000.)*pf2(:)/pek%XF2I(:)
ENDWHERE
!
! -WOODY but not tropical forest-
!
WHERE(gwood(:))
  zfzerostar(:) = ( xaw  - log(zgmest(:)*1000.) )/xbw
ENDWHERE
!
WHERE (gwood(:).AND.pek%LSTRESS(:))
  zgmestn(:) = zgmest(:)
ENDWHERE
WHERE(gwood(:).AND.(.NOT.pek%LSTRESS(:)))
  zgmestn(:) = exp(xasw - xbw*zfzerostar(:))/1000.
ENDWHERE
!
WHERE (gwood(:).AND.gf2_inf_f2i(:)) 
  zgmestn(:) = zgmestn(:)*pf2(:)/pek%XF2I(:)
ENDWHERE
!
WHERE(gwood(:))
  zwork(:) = max( xdenom_min, zgmestn(:) )
  zfzeron(:) = (xasw - log(zwork(:)*1000.))/xbw
ENDWHERE
!
WHERE(gwood(:).AND.(.NOT.gf2_inf_f2i(:)).AND.pek%LSTRESS(:))
  zfzero(:) = zfzerostar(:)
  zfzero(:) = zfzero(:) - (zfzero(:)-zfzeron(:))*(1.-pf2(:))/(1.-pek%XF2I(:))  
ENDWHERE    
WHERE(gwood(:).AND.(.NOT.gf2_inf_f2i(:)).AND.(.NOT.pek%LSTRESS(:)))
  zfzero(:) = zfzerostar(:)
  zgmest(:) = zgmest(:) - (zgmest(:)-zgmestn(:))*(1.-pf2(:))/(1.-pek%XF2I(:))  
ENDWHERE    
!
WHERE(gwood(:).AND.gf2_inf_f2i(:))
  zfzero(:) = min(.95, zfzeron(:))
  zgmest(:) = zgmestn(:)
ENDWHERE  
!
! -Tropical Forest-
!
WHERE(gtrop(:))
 zfzero(:) = pk%XFZERO(:)
 zgmest(:) = pek%XGMES(:)*pf2(:)
ENDWHERE
!
!-------------------------
!
! compensation point (ppm): temperature response
!
!before optimization (with non log PQDGAMM) : 
!ZGAMMT(:) = PGAMM(:)*PQDGAMM(:)**(0.1*(ZTSPC(:)-25.0))
zwork(:) = (0.1*(ztspc(:)-25.0)) * pk%XQDGAMM(:)
zgammt(:) = pk%XGAMM(:) * exp(zwork(:))
!
! specific humidity deficit (kg kg-1)
!
zdsp(:)   = max( 0.0, pqsat(:) - pqa(:) - pek%XLE(:)*pek%XRESA(:)/(prhoa*xlvtt) )
!
! cosine of solar zenith angle 
!
zxmus(:) = max(cos(pzenith(:)),0.01)
!
!
! Compute temperature response functions:
!
! kg/m2/s
!before optimization (with non log PQDAMAX) : 
!ZANMAX(:) = ( PAMAX(:)*PQDAMAX(:)**(0.1*(ZTSPC(:)-25.0)) ) / ...
zwork(:) = (0.1*(ztspc(:)-25.0)) * pk%XQDAMAX(:)
zanmax(:) = ( pk%XAMAX(:) * exp(zwork(:))  ) &
          / ( (1.0+exp(0.3*(pk%XT1AMAX(:)-ztspc(:))))* (1.0+exp(0.3*(ztspc(:)-pk%XT2AMAX(:)))) )
!
! m/s
!before optimization (with non log PQDGMES) : 
!ZGMEST(:) = ( ZGMEST(:)*PQDGMES(:)**(0.1*(ZTSPC(:)-25.0)) ) / ...
zwork(:) = (0.1*(ztspc(:)-25.0)) * pk%XQDGMES(:)
zgmest(:) = ( zgmest(:) * exp(zwork(:)) ) &
          / ( (1.0+exp(0.3*(pk%XT1GMES(:)-ztspc(:))))*  (1.0+exp(0.3*(ztspc(:)-pk%XT2GMES(:)))) )  
!
!
! Integration over the canopy: SIZE(PABC) increments
! are used to approximate the integral.
!
ztan(:) = 0.0
ztgs(:) = 0.0
ztrdk(:)= 0.0
!
! Unit conversion
! ZANMAX and ZEPSO from kgCO2/m2/s to kgCO2/kgair m/s by dividing by RHOA (kgair/m3)
! ZGAMMT from ppm to kgCO2/kgair
zgammt(:)  = zgammt(:) * xmco2 / xmd * 1e-6
zanmax(:) = zanmax(:) / prhoa
zepso(:)  = pk%XEPSO(:)  / prhoa
!
zia(:)     = psw_rad(:)*xparcf
!
DO jint = 1, SIZE(pabc)
  !
  !  Diffusion of incident radiation:
  !
  IF (io%LTR_ML) THEN
    !
    zabc = 1.
    IF (jint.LT.SIZE(pabc)) zabc = pabc(jint+1)
    zweight = zabc - pabc(jint)
    zxia(:) = piacan_sunlit(:,jint)
    !
  ELSE
    !
    zabc = pabc(jint)
    zweight = ppoi(jint)
    !
    CALL ccetr(zxia,zia,zxmus,zabc,zlai)
    !
    ! PAR at different Gauss  level in micmolphot/m2/s
    !
    piacan(:,jint)= zxia(:)
    !
  ENDIF
  !
  ! Compute conductance and assimilation of CO2: 
  !
  !Extinction of respiration depends on LAI above only for tropical evergreen forest
  zlaitop(:) = 0.
  zzlai(:) = 1.
  IF (io%LTR_ML) THEN         
    WHERE(gwood(:))  
      zlaitop(:) = (1.-(pabc(jint)+zabc)/2.)*zlai(:)
      zzlai(:) = zlai(:)
    ENDWHERE
  ENDIF
  CALL cotwo(pcsp, pf2, zxia, zdsp, zgammt,             &
             zfzero, zepso, zanmax, zgmest, pek%XGC(:), zdmax, &  
             zan0, zgs0, zrdk, zlaitop, zzlai           )
  !
  IF (io%LTR_ML) THEN
    !
    zxia(:) = piacan_shade(:,jint)
    CALL cotwo(pcsp, pf2, zxia, zdsp, zgammt,             &
               zfzero, zepso, zanmax, zgmest, pek%XGC(:), zdmax, &  
               zan0_, zgs0_, zrdk_, zlaitop, zzlai        )
    !
    WHERE (oshade(:))
      !ponderate sum.
      zan0(:)=pfrac_sun(:,jint)*zan0(:)+(1.-pfrac_sun(:,jint))*zan0_(:)
      zrdk(:)=pfrac_sun(:,jint)*zrdk(:)+(1.-pfrac_sun(:,jint))*zrdk_(:)
      zgs0(:)=pfrac_sun(:,jint)*zgs0(:)+(1.-pfrac_sun(:,jint))*zgs0_(:)
    ENDWHERE
    !
  ENDIF
  !
  ! kgCO2/kgair m/s
  ztan(:) = ztan(:) + zan0(:)*zweight
  ztgs(:) = ztgs(:) + zgs0(:)*zweight
  ztrdk(:) = ztrdk(:) + zrdk(:)*zweight
  !
END DO
!
!
! Total assimilation
!
zanf(:)= ztan(:)
!
! Net assimilation over canopy
!
pek%XAN(:) = (1.0-pdelta(:))*(1.0-ppsnv(:)-pffv(:))*zanf(:)*zlai(:)
!
! Dark respiration over canopy (does not depend on radiation, 
! no need to integrate over vertical dimension)
!
presp_leaf(:) = (1.0-pdelta(:))*(1.0-ppsnv(:)-pffv(:))*ztrdk(:)*zlai(:)
!
! Gross primary production over canopy
!
pgpp(:) = pek%XAN(:) + presp_leaf(:)
!
! Cumulated daily net assimilation over canopy (kgCO2/m2/day)
!
pek%XANDAY(:) = pek%XANDAY(:) + pek%XAN(:) * ptstep * prhoa
!
! Adjust maximum leaf assimilation:
!
pek%XANFM(:) = max( zanf(:), pek%XANFM(:) )
!
! Total conductance over canopy 
!
zxtgs(:) = ztgs(:)*zlai(:)
!
! Canopy resistance from Ags:
!
prs(:) = min( 1.0/(zxtgs(:)+xdenom_min), xrs_max)
!
prs(:) = max( prs(:), zrs_min)
!
IF (lhook) CALL dr_hook('COTWORES',1,zhook_handle)
!
END SUBROUTINE cotwores