vegetation_evol.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 vegetation_evol(IO, DTI, PK, PEK, OAGRIP, PTSTEP, KMONTH, KDAY, PTIME, &
                               PLAT, PRHOA, P_CO2, ISSK, PRESP_BIOMASS_INST, PSWDIR)  
!   ###############################################################
!!****  *VEGETATION EVOL*
!!
!!    PURPOSE
!!    -------
!
!     performs the time evolution of vegetation parameters
!     at solar midnight in the case of interactive vegetation (ISBA-Ags)
!              
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!      
!!    none
!!
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      V. Masson          * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/03/03 
!!      P. Le Moigne 12/2004 : NIT version 
!!      P Le Moigne  09/2005 : AGS modifs of L. Jarlan
!!      A.L. Gibelin 04/2009 : BIOMASS and RESP_BIOMASS arrays
!!      A.L. Gibelin 04/2009 : Add NCB option 
!!      D. Carrer    01/2012 : representation of nitrogen dilution fct of CO2 (from Calvet et al. 2008)
!!      B. Decharme  05/2012 : Optimization and ISBA-DIF coupling
!!      C. Delire    01/2014 : IBIS respiration for tropical evergreen
!!      R. Seferian  05/2015 : expanding of Nitrogen dilution option to the complete formulation proposed by Yin et al. GCB 2002 
!!Seferian & Delire  06/2015 : accouting for living woody biomass respiration (expanding work of E Joetzjer to all woody PFTs) 
!!      B. Decharme    01/16 : Bug when vegetation veg, z0 and emis are imposed whith interactive vegetation
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_data_isba_n, ONLY : data_isba_t
USE modd_isba_n, ONLY : isba_p_t, isba_pe_t
!
USE modd_sso_n, ONLY : sso_t
!
USE modd_co2v_par,       ONLY : xmc, xmco2, xpcco2, xrespfactor_nit,       &
                                xcoeff_maint_resp_zero, xslope_maint_resp, &
                                xparcf, xdiludec, itransfert_esg
USE modd_csts,           ONLY : xday, xtt, xmd
!
USE modi_albedo
USE modi_laigain
USE modi_lailoss
USE modi_nitro_decline
USE modi_emis_from_veg
USE modi_veg_from_lai
USE modi_z0v_from_lai
USE modi_subscale_z0eff
USE modd_type_date_surf
USE modd_data_cover_par, ONLY : nvegtype_ecosg, nvegtype, &
                                nvt_tebd, nvt_trbe, nvt_bone,   &
                                nvt_trbd, nvt_tebe, nvt_tene,   &
                                nvt_bobd, nvt_bond, nvt_shrb,   &
                                nvt_trbe
USE modd_surf_par
!
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(data_isba_t), INTENT(INOUT) :: DTI
TYPE(isba_p_t), INTENT(INOUT) :: PK
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
!
LOGICAL,              INTENT(IN)    :: OAGRIP  ! agricultural practices
!
REAL,                 INTENT(IN)    :: PTSTEP  ! time step
INTEGER,              INTENT(IN)    :: KMONTH  ! current month
INTEGER,              INTENT(IN)    :: KDAY    ! current day
REAL,                 INTENT(IN)    :: PTIME   ! current time since midnight
REAL,   DIMENSION(:), INTENT(IN)    :: PLAT    ! latitude of each grid point
REAL,   DIMENSION(:), INTENT(IN)    :: PRHOA   ! air density
!
REAL,   DIMENSION(:), INTENT(IN)    :: P_CO2 ! CO2 concentration [ppmm]
!
TYPE(sso_t), INTENT(INOUT) :: ISSK
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: PRESP_BIOMASS_INST ! instantaneous respiration of biomass (kgCO2/kgair m/s)
!
REAL, DIMENSION(:),   INTENT(IN),   OPTIONAL :: PSWDIR    ! Global incoming shortwave radiation (W m-2)
!
!*      0.2    declarations of local parameter
!
REAL, PARAMETER                   :: ZCOEF1 = 10.0
REAL, PARAMETER                   :: ZCOEF2 = 25.0
REAL, PARAMETER                   :: ZDEPTH = 1.0   !Temp depth m
!
REAL, PARAMETER                   :: ZWOOD_IBIS=0.0125
REAL, PARAMETER                   :: ZROOT_IBIS=1.25 
REAL, PARAMETER                   :: ZCIBIS1   =3500.
REAL, PARAMETER                   :: ZCIBIS2   =1./288.
REAL, PARAMETER                   :: ZNDAY     =365.
!
REAL, PARAMETER                   :: ZCDILU1 = -0.048
REAL, PARAMETER                   :: ZCDILU2 = 6.3
REAL, PARAMETER                   :: ZCDILU3 = 371.
! Required for Yin et al., nitrogen dilu param
REAL, PARAMETER                   :: ZPHOTON    = 2.010402e-3 ! conversion coef for W m-2 in photon m-2
REAL, PARAMETER                   :: ZDEPTH_VEG = 0.40        !Depth in meters for daily temperature
REAL, PARAMETER                   :: ZTEMP_VEG  = 23.         !Average temperature of the vegetation
REAL, PARAMETER                   :: ZDECIDUS   = 0.75        !Coef for decidus trees
!
!*      0.3    declarations of local variables
!
REAL, DIMENSION(SIZE(PEK%XRESP_BIOMASS,1),SIZE(PEK%XRESP_BIOMASS,2)) :: ZRESP_BIOMASS_LAST ! biomass at t-1 (kg_DM/m2/day)
REAL,    DIMENSION(SIZE(PEK%XLAI,1))    :: ZBIOMASS_LEAF   ! temporary leaf biomass 
REAL,    DIMENSION(SIZE(PEK%XLAI,1))    :: ZBSLAI_NITRO    ! (Calvet et al. 2008) ratio of biomass to LAI
                                                     ! with representation of nitrogen dilution
REAL,    DIMENSION(SIZE(PEK%XLAI,1)) :: ZCO2, ZCNA_NITRO   ! fct of CO2        
REAL,    DIMENSION(SIZE(PEK%XLAI,1)) :: ZPARAM
REAL,    DIMENSION(SIZE(PEK%XLAI,1)) :: ZHTREE, ZSAPFRAC   ! tree height & sap fraction used for estimation of 
                                                     ! sapwood fraction
!
REAL                              :: ZLOG2, ZWORK
!
REAL, DIMENSION(SIZE(PEK%XTG,1))      :: ZTG_VEG      ! surface temperature   (C)
REAL, DIMENSION(SIZE(PEK%XTG,1))      :: ZTG_SOIL     ! soil temperature   (C)
REAL, DIMENSION(SIZE(PEK%XTG,1))      :: ZDG_SOIL     ! soil depth for DIF (m)
REAL                              :: ZWGHT_SOIL   ! Weight for DIF (m)
!
LOGICAL, DIMENSION(SIZE(PEK%XLAI,1))    :: GWOOD,GHERB
LOGICAL, DIMENSION(SIZE(PEK%XLAI,1))    :: GMASK_AGRI
LOGICAL                           :: GMASK
INTEGER                           :: INI, INL, JI, JL, IDEPTH, JTYPE
!
REAL,    DIMENSION(SIZE(PK%XVEGTYPE_PATCH,1),SIZE(PK%XVEGTYPE_PATCH,2)) :: ZPARAM_TYPE
!
! * Azote
REAL,    DIMENSION(SIZE(PEK%XLAI,1)) :: ZFERT
!
REAL :: ZDILUDEC
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-----------------------------------------------------------------
!
!*      1.     Preliminaries
!              -------------
!
IF (lhook) CALL dr_hook('VEGETATION_EVOL',0,zhook_handle)
!
ini=SIZE(pek%XTG,1)
inl=SIZE(pek%XTG,2)
!
zlog2 = log(2.0)
!
ztg_soil(:) = 0.0
ztg_veg(:) = 0.0
!
! Define herbaceous and woody patches
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(:))
!
! Mask where vegetation evolution is performed (just before solar midnight)
gmask = ( ptime - ptstep < 0. ) .AND. ( ptime >= 0. )
!
! Save RESP_BIOMASS at t-1
IF (gmask) THEN
  pek%XRESP_BIOMASS(:,1) = 0.0
  zresp_biomass_last(:,:) = 0.0
ELSE
  pek%XRESP_BIOMASS(:,1) = pek%XRESP_BIOMASS(:,1) + presp_biomass_inst(:,1) * (ptstep*prhoa(:)*xmc)/(xpcco2*xmco2)
  zresp_biomass_last(:,:) = pek%XRESP_BIOMASS(:,:)
ENDIF
!
!*      2.     Interactive vegetation
!              ----------------------
!
!  LAI daily mortality and assimilation
!
zbiomass_leaf(:) = pek%XBIOMASS(:,1)
!
IF (gmask) THEN
  !        
  pk%XINCREASE (:,:) = 0.0
  pk%XTURNOVER(:,:) = 0.0
  zbslai_nitro(:  ) = pk%XBSLAI_NITRO(:) 
  !
  IF(io%LNITRO_DILU)THEN
!
!     * Compute Vegetation temperature
!       We use the temperature of the second layer of the soil (<40cm)
!       since the parametrization employs a daily temperature
!
    IF(io%CISBA/='DIF')THEN        
      ztg_veg(:) = pek%XTG(:,2)
    ELSE 
      DO ji=1,ini
         idepth=pk%NWG_LAYER(ji)
         zdg_soil(ji)=min(zdepth_veg,pk%XDG(ji,idepth))
      ENDDO  
      DO jl=1,inl
        DO ji=1,ini     
          zwght_soil=min(pk%XDZG(ji,jl),max(0.0,zdg_soil(ji)-pk%XDG(ji,jl)+pk%XDZG(ji,jl)))        
          ztg_veg(ji)=ztg_veg(ji)+pek%XTG(ji,jl)*zwght_soil/zdg_soil(ji)
        ENDDO
      ENDDO 
    ENDIF
!
    zparam(:) = 0.0
    zfert(:) = 0.0
    DO jtype=1,SIZE(pk%XVEGTYPE_PATCH,2)
      IF (nvegtype==nvegtype_ecosg) THEN
        zdiludec = xdiludec(itransfert_esg(jtype))
      ELSE
        zdiludec = xdiludec(jtype)
      ENDIF
      DO ji = 1,ini
        zparam_type(ji,jtype) = zdiludec * (zdecidus + 1.1 * zphoton * xparcf * pswdir(ji)       &
                              + (ztg_veg(ji)-xtt)/ztemp_veg - 0.33 * zfert(ji))                         &
                              + (1 - zdiludec) * (1.1 * zphoton * xparcf * pswdir(ji) &
                              + (ztg_veg(ji)-xtt)/ztemp_veg - 0.33 * zfert(ji))
        zparam(ji) = zparam(ji) + zparam_type(ji,jtype) * pk%XVEGTYPE_PATCH(ji,jtype)
      ENDDO 
    ENDDO  

    WHERE((pek%XCE_NITRO(:)*pek%XCNA_NITRO(:)+pek%XCF_NITRO(:))/=0.0.AND.pek%XCNA_NITRO(:)/=0.0)
      zco2(:) = p_co2(:)*(xmd/(1.e-6*xmco2))  ! (ppmm ->  ppm)
      zcna_nitro(:) = pek%XCNA_NITRO(:) * &
                   exp(zcdilu1*exp(zparam(:)-pek%XCNA_NITRO(:)/zcdilu2) * alog(max(1.,zco2(:)/zcdilu3)))
      zbslai_nitro(:) = 1. / (pek%XCE_NITRO(:)*zcna_nitro(:)+pek%XCF_NITRO(:))
    ENDWHERE
!
  ENDIF
!    
  IF(any(pek%XLAI(:)/=xundef))THEN
    CALL nitro_decline(io, pk, pek, gwood, zbslai_nitro, plat, zbiomass_leaf)
    CALL laigain(zbslai_nitro, pek, zbiomass_leaf)
  ENDIF
  !  
  ! reinitialise  PEK%XANDAY(:) and PEK%XANFM(:) 
  pek%XANDAY(:)=0.0
  pek%XANFM(:) =0.0
  !
ENDIF
!
!
! * soil temperature in K (over 1m depth for DIF)
!
ztg_veg(:) = pek%XTG(:,1)
!
IF(io%CISBA/='DIF')THEN        
  ztg_soil(:) = pek%XTG(:,2)
ELSE       
  DO ji=1,ini
    idepth=pk%NWG_LAYER(ji)
    zdg_soil(ji)=min(zdepth,pk%XDG(ji,idepth))
  ENDDO  
  DO jl=1,inl
    DO ji=1,ini     
      zwght_soil=min(pk%XDZG(ji,jl),max(0.0,zdg_soil(ji)-pk%XDG(ji,jl)+pk%XDZG(ji,jl)))        
      ztg_soil(ji)=ztg_soil(ji)+pek%XTG(ji,jl)*zwght_soil/zdg_soil(ji)
    ENDDO
  ENDDO 
ENDIF
!
!
! * Respiration of structural biomass pools
!
WHERE(gwood(:))
  ! IBIS respiration with either respiration factor rwood=0.0125 - otherwise rroot=1.25 
  ! (Kucharik et al, 2000, eq 6-8) Soil temp in K         
  pek%XRESP_BIOMASS(:,2) = pek%XRESP_BIOMASS(:,2) + pek%XBIOMASS(:,2) * ptstep &
                              * max(0.,zroot_ibis*exp(zcibis1*(zcibis2-1./ztg_veg(:)))/(znday*xday)) 
ELSEWHERE 
  pek%XRESP_BIOMASS(:,2) = pek%XRESP_BIOMASS(:,2) + pek%XBIOMASS(:,2) * xrespfactor_nit    &
                              * exp((zlog2/zcoef1)*(ztg_veg(:)-xtt-zcoef2)) * ptstep  
  ! before optimization                   * 2.0**((PEK%XTG(:,2)-XTT-ZCOEF2)/ZCOEF1) * PTSTEP               
ENDWHERE
!
IF (io%CPHOTO == 'NIT') THEN
  !
  pek%XRESP_BIOMASS(:,3) = pek%XRESP_BIOMASS(:,3) + pek%XBIOMASS(:,3) * xrespfactor_nit &
                            * exp((zlog2/zcoef1)*(ztg_soil(:)-xtt-zcoef2)) * ptstep  
  ! before optimization                   * 2.0**((PEK%XTG(:,2)-XTT-ZCOEF2)/ZCOEF1) * PTSTEP               
  !
ELSEIF (io%CPHOTO == 'NCB') THEN
  !
  pek%XRESP_BIOMASS(:,2) = min(pek%XRESP_BIOMASS(:,2), pek%XBIOMASS(:,2))
  ! 
  pek%XRESP_BIOMASS(:,3) = pek%XRESP_BIOMASS(:,3) + pek%XBIOMASS(:,3) * &
            max( 0., xcoeff_maint_resp_zero * (1. + xslope_maint_resp*(ztg_veg(:)-xtt))) * ptstep  
  pek%XRESP_BIOMASS(:,3) = min(pek%XRESP_BIOMASS(:,3), pek%XBIOMASS(:,3))
  ! 
  WHERE(gwood(:))
    ! Resp IBIS (Soil temp in K)
    pek%XRESP_BIOMASS(:,4) = pek%XRESP_BIOMASS(:,4) + pek%XBIOMASS(:,4) * ptstep &
                        * max(0.,zroot_ibis * exp(zcibis1*(zcibis2-1./ztg_soil(:)))/(znday*xday))
  ELSEWHERE 
    pek%XRESP_BIOMASS(:,4) = pek%XRESP_BIOMASS(:,4) + pek%XBIOMASS(:,4) * &
             max( 0., xcoeff_maint_resp_zero * (1. + xslope_maint_resp*(ztg_soil(:)-xtt))) * ptstep  
  ENDWHERE
  !
  pek%XRESP_BIOMASS(:,4) = min(pek%XRESP_BIOMASS(:,4), pek%XBIOMASS(:,4))
  !
  WHERE( (gwood(:)).AND.(pek%XBIOMASS(:,5)>0.) )
    ! IBIS estimation of sapwood fraction based on the height of tree, sapspeed and 
    ! max transpiration rates. Conversion from DM to C. To be changed with DGVM.  (Soil temp in K)        
    zhtree(:) = 2.5*0.75*(pek%XBIOMASS(:,1)+pek%XBIOMASS(:,2)+pek%XBIOMASS(:,3)+&
                          pek%XBIOMASS(:,4)+pek%XBIOMASS(:,5)+pek%XBIOMASS(:,6))*0.4
    zsapfrac(:) = min(0.5, max(0.05,0.0025/25.*zhtree(:)*0.75*400/(pek%XBIOMASS(:,5)*0.4)))
    !ZSAPFRAC(:) = 0.5
    
    pek%XRESP_BIOMASS(:,5) = pek%XRESP_BIOMASS(:,5) + pek%XBIOMASS(:,5) * zsapfrac(:) * ptstep &
                               * max(0.,zwood_ibis*exp(zcibis1*(zcibis2-1./ztg_veg(:)))/(znday*xday))
    pek%XRESP_BIOMASS(:,5) = min(pek%XRESP_BIOMASS(:,5), pek%XBIOMASS(:,5))
  ELSEWHERE
    pek%XRESP_BIOMASS(:,5) = 0.0
  ENDWHERE
  !
ENDIF
!
! * Instantaneous respiration (kgCO2/kgair m/s)
!
DO jl=2,SIZE(pek%XRESP_BIOMASS(:,:),2)
   presp_biomass_inst(:,jl) = (pek%XRESP_BIOMASS(:,jl) - zresp_biomass_last(:,jl)) &
                                 * xpcco2*xmco2/(ptstep*prhoa(:)*xmc)                              
ENDDO
!
!*      3.     Agricultural practices
!              ----------------------
!
IF (oagrip) THEN
  !
  gmask_agri(:) = .false.
  WHERE ( pek%TSEED(:)%TDATE%MONTH /= nundef .AND. ( kmonth < pek%TSEED(:)%TDATE%MONTH .OR. &
         (kmonth == pek%TSEED(:)%TDATE%MONTH .AND. kday < pek%TSEED(:)%TDATE%DAY) ) )  gmask_agri(:) = .true.
  WHERE ( pek%TREAP(:)%TDATE%MONTH /= nundef .AND. ( kmonth > pek%TREAP(:)%TDATE%MONTH .OR. &
         (kmonth == pek%TREAP(:)%TDATE%MONTH .AND. kday >= pek%TREAP(:)%TDATE%DAY) ) ) gmask_agri(:) = .true. 
  !
  WHERE (gmask_agri(:))
    pek%XLAI(:)             = pek%XLAIMIN(:)
    zbiomass_leaf(:)    = pek%XLAI(:) * zbslai_nitro(:)
  END WHERE

  WHERE (gmask_agri(:))
    pek%XBIOMASS(:,1)       = 0.0
    pek%XBIOMASS(:,2)       = 0.0
    pek%XBIOMASS(:,3)       = 0.0
    pek%XRESP_BIOMASS(:,2)  = 0.0
    pek%XRESP_BIOMASS(:,3)  = 0.0
  END WHERE
  !
  IF (io%CPHOTO == 'NCB') THEN
    !
    WHERE (gmask_agri(:)) 
      pek%XBIOMASS(:,4)       = 0.0
      pek%XBIOMASS(:,5)       = 0.0
      pek%XBIOMASS(:,6)       = 0.0
      pek%XRESP_BIOMASS(:,4)  = 0.0
    END WHERE
    !
  ENDIF
  !
ENDIF
!
!*      4.     Physical parameters depending on vegetation
!              -------------------------------------------
!
IF (gmask) THEN
  !
  ! Evolution of vegetation fraction and roughness length due to LAI change
  IF(.NOT.dti%LIMP_Z0) THEN
    WHERE( pek%XVEG(:) > 0. ) &
      pek%XZ0 (:) = z0v_from_lai(pek%XLAI(:),pk%XH_TREE(:),pk%XVEGTYPE_PATCH(:,:),io%LAGRI_TO_GRASS) 
  ENDIF
  IF(.NOT.dti%LIMP_VEG) THEN
    WHERE( pek%XVEG(:) > 0. ) &
      pek%XVEG(:) = veg_from_lai(pek%XLAI(:),pk%XVEGTYPE_PATCH(:,:),io%LAGRI_TO_GRASS)
  ENDIF
  !
  ! Evolution of radiative parameters due to vegetation fraction change
  IF(.NOT.dti%LIMP_EMIS) THEN
    WHERE( pek%XVEG(:) > 0. ) pek%XEMIS(:)= emis_from_veg(pek%XVEG(:),pk%XVEGTYPE_PATCH(:,:))
  ENDIF
  !
  CALL albedo(io%CALBEDO, pek )  
  !
  ! Evolution of effective roughness length due to new surface roughness length
  !
  IF (ASSOCIATED(issk%XAOSIP)) THEN
    IF (SIZE(issk%XAOSIP)>0) THEN
      CALL subscale_z0eff(issk,pek%XZ0(:),.false. )
    ENDIF
  ENDIF
  !
ENDIF
!
IF (lhook) CALL dr_hook('VEGETATION_EVOL',1,zhook_handle)
!-----------------------------------------------------------------
!
END SUBROUTINE vegetation_evol