cotwoinitn.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 cotwoinit_n (IO, S, PK, PEK, PCO2  )  
!     #######################################################################
!
!!****  *COTWOINIT*  
!!
!!    PURPOSE
!!    -------
!
!     Initialize model to calculate net assimilation of 
!     CO2 and leaf conductance.
!              
!!**  METHOD
!!    ------
!     Calvet at al (1998) [from model of Jacobs(1994)]
!!
!!    EXTERNAL
!!    --------
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!      
!!    USE MODD_CO2V_PAR
!!    USE MODI_COTWO  
!!
!!    REFERENCE
!!    ---------
!!
!!    Calvet et al. (1998)
!!      
!!    AUTHOR
!!    ------
!!
!!      A. Boone           * Meteo-France *
!!      (following Belair)
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    27/10/97
!!      (V. Rivalland) 10/04/02  Add: PK%XAH and PK%XBH coefficients for
!!                               herbaceous water stress response
!!      (P. LeMoigne) 03/2004:   computation of zgmest in SI units
!!      (P. LeMoigne) 10/2004:   possibility of 2 different FZERO
!!      (L. Jarlan)   10/2004:   initialization of DMAX
!!      P Le Moigne   09/2005    AGS modifs of L. Jarlan
!!      S. Lafont     03/2009    change unit of AMAX
!!      A.L. Gibelin  04/2009    TAU_WOOD for NCB option 
!!      A.L. Gibelin  04/2009    Suppress useless GPP and RDK arguments 
!!      A.L. Gibelin  07/2009    Suppress PPST and PPSTF as outputs
!!      B. Decharme   05/2012    Optimization
!!      R. Alkama     05/2012    add 7 new vegtype (19  instead 12)
!!      C. Delire     01/2014    Define a dummy LAI from top and total lai for Dark respiration 
!!
!-------------------------------------------------------------------------------
!
USE modd_data_cover_par, ONLY : 
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_p_t, isba_s_t, isba_pe_t
!
USE modd_data_cover_par, ONLY : nvegtype_ecosg, nvegtype, nvt_c3, nvt_c3w, nvt_c3s, nvt_c4, nvt_irr, &
                                nvt_trog, nvt_tebd, nvt_bone, nvt_trbe, nvt_trbd,    &
                                nvt_tebe, nvt_tene, nvt_bobd, nvt_bond, nvt_shrb, nvt_gras
USE modd_csts,           ONLY : xmd
USE modd_co2v_par,       ONLY : xtopt, xfzero1, xfzero2, xfzerotrop, xepso, xgamm, xqdgamm, &
                                  xqdgmes, xt1gmes, xt2gmes, xamax, itransfert_esg,         &
                                  xqdamax, xt1amax, xt2amax, xah, xbh,            &
                                  xdspopt, xiaopt, xaw, xbw, xmco2, xmc, xtau_wood  
! 
USE modi_cotwo  
!
!*       0.     DECLARATIONS
!               ------------
!
!
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_p_t), INTENT(INOUT) :: PK
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
TYPE(isba_s_t), INTENT(INOUT) :: S
!
REAL,DIMENSION(:),INTENT(IN)  :: PCO2
!                                     PK%XGMES     = mesophyll conductance (m s-1)
!                                     PCO2      = atmospheric CO2 concentration
!
!*      0.2    declaration of local variables
!
INTEGER                           :: JCLASS    ! indexes for loops
INTEGER                           :: ICLASS    ! indexes for loops
INTEGER                           :: ICO2TYPE  ! type of CO2 vegetation
INTEGER                           :: IRAD      ! with or without new radiative transfer
!
REAL, DIMENSION(SIZE(PK%XANMAX))     :: ZGS, ZGAMMT, ZTOPT, ZANMAX, ZGMEST, ZGPP, ZRDK, ZEPSO
!                                    ZTOPT     = optimum  temperature for compensation 
!                                                point
!                                    ZANMAX    = maximum photosynthesis rate
!                                    ZGS       = leaf conductance
!                                    ZGAMMT    = temperature compensation point
!                                    ZGPP      = gross primary production
!                                    ZRDK      = dark respiration
!
!
REAL, DIMENSION(SIZE(PK%XANMAX))     :: ZCO2INIT3, ZCO2INIT4, ZCO2INIT5, ZCO2INIT2,ZCO2INIT1
!                                    working arrays for initializing surface 
!                                    temperature, saturation deficit, global radiation,
!                                    optimum temperature for determining maximum 
!                                    photosynthesis rate, and soil water stress (none)
REAL, DIMENSION(SIZE(PK%XDMAX))      :: ZDMAX
REAL, DIMENSION(SIZE(PK%XDMAX))      :: ZWORK
!                                    Local variable in order to initialise DMAX
!                                    following Calvet, 2000 (AST or LST cases)
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('COTWOINIT_N',0,zhook_handle)
!
ztopt(:) = 0.
pk%XFZERO (:) = 0.
pk%XEPSO  (:) = 0.
pk%XGAMM  (:) = 0.
pk%XQDGAMM(:) = 0.
pk%XQDGMES(:) = 0.
pk%XT1GMES(:) = 0.
pk%XT2GMES(:) = 0.
pk%XAMAX  (:) = 0.
pk%XQDAMAX(:) = 0.
pk%XT1AMAX(:) = 0.
pk%XT2AMAX(:) = 0.
pk%XTAU_WOOD(:) = 0.
!
pk%XAH    (:) = 0.
pk%XBH    (:) = 0.
!
zepso(:) = 0.
zgpp(:) = 0.
zrdk(:) = 0.
zgammt(:) = 0.
zanmax(:) = 0.
zgmest(:) = 0.
zco2init3(:) = 0.
zco2init4(:) = 0.
zco2init5(:) = 0.
!
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
!
!
! INITIALIZE VARIOUS PARAMETERS FOR CO2 MODEL:
! --------------------------------------------
! as a function of CO2 vegetation class, C3=>1, C4=>2
!
DO jclass=1,nvegtype
  !
  IF (jclass==nvt_c4 .OR. jclass==nvt_irr .OR. jclass==nvt_trog) THEN
    ico2type = 2   ! C4 type
  ELSE
    ico2type = 1   ! C3 type
  END IF
  IF(io%LAGRI_TO_GRASS.AND.(jclass==nvt_c4 .OR. jclass==nvt_irr)) ico2type = 1
  IF (io%LTR_ML) THEN
    irad = 1   ! running with new radiative transfer
  ELSE
    irad = 2
  ENDIF
  !
  ztopt(:) = ztopt(:) + xtopt(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  IF((jclass==nvt_tebd) .OR. (jclass==nvt_bone) .OR.                         &
    (jclass==nvt_trbd) .OR. (jclass==nvt_tebe) .OR. (jclass==nvt_tene) .OR. &
    (jclass==nvt_bobd) .OR. (jclass==nvt_bond) .OR. (jclass==nvt_shrb)) THEN
    pk%XFZERO (:) = pk%XFZERO (:) + ((xaw - log(pek%XGMES(:)*1000.0))/xbw)*pk%XVEGTYPE_PATCH(:,jclass)
  ELSE IF (jclass==nvt_trbe) THEN
    pk%XFZERO (:) = pk%XFZERO (:) + xfzerotrop(irad) * pk%XVEGTYPE_PATCH(:,jclass)
  ELSE
    pk%XFZERO (:) = pk%XFZERO (:) + xfzero2(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  ENDIF
  !
  pk%XEPSO  (:) = pk%XEPSO  (:) + xepso(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  pk%XGAMM  (:) = pk%XGAMM  (:) + xgamm(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  pk%XQDGAMM(:) = pk%XQDGAMM(:) + xqdgamm(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  pk%XQDGMES(:) = pk%XQDGMES(:) + xqdgmes(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  pk%XT1GMES(:) = pk%XT1GMES(:) + xt1gmes(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  pk%XT2GMES(:) = pk%XT2GMES(:) + xt2gmes(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  pk%XQDAMAX(:) = pk%XQDAMAX(:) + xqdamax(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  pk%XT1AMAX(:) = pk%XT1AMAX(:) + xt1amax(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  pk%XT2AMAX(:) = pk%XT2AMAX(:) + xt2amax(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  pk%XAH    (:) = pk%XAH    (:) + xah(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  pk%XBH    (:) = pk%XBH    (:) + xbh(ico2type) * pk%XVEGTYPE_PATCH(:,jclass)
  !
  IF(io%LAGRI_TO_GRASS.AND.(jclass==nvt_c3 .OR. jclass==nvt_c3w .OR. jclass==nvt_c3s .OR. &
                            jclass==nvt_c4 .OR. jclass==nvt_irr))THEN
    iclass=nvt_gras
  ELSE
    iclass=jclass
  ENDIF    
  !
  IF (nvegtype==nvegtype_ecosg) THEN
    pk%XTAU_WOOD(:) = pk%XTAU_WOOD(:) + xtau_wood(itransfert_esg(iclass)) * pk%XVEGTYPE_PATCH(:,jclass)
    pk%XAMAX    (:) = pk%XAMAX    (:) + xamax(itransfert_esg(iclass)) * pk%XVEGTYPE_PATCH(:,jclass)          
  ELSE
    pk%XTAU_WOOD(:) = pk%XTAU_WOOD(:) + xtau_wood(iclass) * pk%XVEGTYPE_PATCH(:,jclass)
    pk%XAMAX    (:) = pk%XAMAX    (:) + xamax(iclass) * pk%XVEGTYPE_PATCH(:,jclass)
  ENDIF
  !
END DO
!
pk%XQDGAMM(:)=log(pk%XQDGAMM(:))
pk%XQDGMES(:)=log(pk%XQDGMES(:))
pk%XQDAMAX(:)=log(pk%XQDAMAX(:))
!
!
! INITIALIZE VARIOUS VARIABLES FOR CO2 MODEL:
! -------------------------------------------
!
!
! compute temperature responses:
!
!before optimization (with non log PK%XQDGAMM) : 
!ZGAMMT(:) = PK%XGAMM(:)*(PK%XQDGAMM(:)**(0.1*(ZTOPT(:)-25.0)))
zwork(:) = (0.1*(ztopt(:)-25.0)) * pk%XQDGAMM(:)
zgammt(:) = pk%XGAMM(:)*exp(zwork(:))
!
!before optimization (with non log PK%XQDAMAX) :
!ZANMAX(:) = ( PK%XAMAX(:)*PK%XQDAMAX(:)**(0.1*(ZTOPT(:)-25.0)) ) / ...
zwork(:) = (0.1*(ztopt(:)-25.0)) * pk%XQDAMAX(:)
zanmax(:) = ( pk%XAMAX(:)*exp(zwork(:)) )                   &
               /( (1.0+exp(0.3*(pk%XT1AMAX(:)-ztopt(:))))*  &
                  (1.0+exp(0.3*(ztopt(:)-pk%XT2AMAX(:)))) )  
!
!before optimization (with non log PK%XQDGMES) :
!ZGMEST(:) = ( PEK%XGMES(:)*PK%XQDGMES(:)**(0.1*(ZTOPT(:)-25.0)) )    &
zwork(:) = (0.1*(ztopt(:)-25.0)) * pk%XQDGMES(:)
zgmest(:) = ( pek%XGMES(:)*exp(zwork(:)) )                   &
               /( (1.0+exp(0.3*(pk%XT1GMES(:)-ztopt(:))))*  &
                  (1.0+exp(0.3*(ztopt(:)-pk%XT2GMES(:)))) )  
!
!
! initialize other variables: (using optimum values for some variables)
!
zco2init3(:) = xdspopt
zco2init4(:) = xiaopt
zco2init5(:) = 1.0
!
! Define a dummy LAI from top (zco2init2=0.1) and total lai (zco2init=1) for Dark respiration extinction parameterization 
!
zco2init2(:) = 0.1
zco2init1(:) = 1.0
!
! Add soil moisture stress effect to leaf conductance:
!
zgmest(:) = zgmest(:)*zco2init5(:)
!
! Initialise DMAX following Calvet (2000) in the case of 'AST' or 'LST' photosynthesis option
!
IF(io%CPHOTO/='NON') THEN
   zdmax(:) = exp((log(zgmest(:)*1000.)-pk%XAH(:))/pk%XBH(:))/1000.
ELSE
   zdmax(:) = pk%XDMAX(:)
ENDIF
!
! Compute maximum/initial/optimum net assimilation of CO2:
!
! Unit conversion with a constant value of 1.2 for PRHOA as it is not known here
! ZANMAX and ZEPSO from kgCO2/m2/s to kgCO2/kgair m/s by dividing by RHOA (kgair/m3)
! ZGAMMT from ppm to kgCO2/kgair
zanmax(:)=zanmax(:)/1.2
zepso(:)=pk%XEPSO(:)/1.2
zgammt(:)=zgammt(:)*xmco2/xmd*1e-6
!
 CALL cotwo(pco2, zco2init5, zco2init4, zco2init3, zgammt, &
           pk%XFZERO(:), zepso, zanmax, zgmest, pek%XGC, zdmax,     &
           pk%XANMAX(:), zgs, zrdk, zco2init2, zco2init1        )                     
! change by sebastien PK%XEPSO change into ZEPSO for units consistency
!
!
!
IF (lhook) CALL dr_hook('COTWOINIT_N',1,zhook_handle)
!
END SUBROUTINE cotwoinit_n