! =================================================================================================================================
! MODULE 	: stomate_data
!
! CONTACT      : orchidee-help _at_ listes.ipsl.fr
!
! LICENCE      : IPSL (2006)
! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF         "stomate_data" module defines the values about the PFT parameters. It will print
!! the values of the parameters for STOMATE in the standard outputs. 
!!
!!\n DESCRIPTION: None 
!!
!! RECENT CHANGE(S): Sonke Zaehle: Reich et al, 1992 find no statistically significant differences 
!!                  between broadleaved and coniferous forests, specifically, the assumption that grasses grow 
!!                  needles is not justified. Replacing the function with the one based on Reich et al. 1997. 
!!                  Given that sla=100cm2/gDW at 9 months, sla is:
!!                  sla=exp(5.615-0.46*ln(leaflon in months))
!!
!! REFERENCE(S)	: None
!!
!! SVN          :
!! $HeadURL: svn://forge.ipsl.fr/orchidee/trunk/ORCHIDEE/src_stomate/stomate_data.f90 $
!! $Date: 2023-11-22 20:28:54 +0100 (mer. 22 nov. 2023) $
!! $Revision: 8316 $
!! \n
!_ ================================================================================================================================

MODULE stomate_data

  ! modules used:

  USE constantes
  USE time, ONLY : one_day, dt_sechiba, one_year
  USE pft_parameters
  USE dynamic_parameters
  USE defprec
  USE function_library, ONLY: get_printlev  

  IMPLICIT NONE

  INTEGER(i_std), SAVE, PRIVATE                :: printlev_loc   !! local printlev for this module
!$OMP THREADPRIVATE(printlev_loc)

  INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: hori_index     !! Move to Horizontal indices
!$OMP THREADPRIVATE(hori_index)

  INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: horipft_index  !! Horizontal + PFT indices
!$OMP THREADPRIVATE(horipft_index)

  INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: horican_index  !! Horizontal + canopy levels indices
!$OMP THREADPRIVATE(horican_index)
  
  INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: horicut_index  !! Horizontal + cut times indices
!$OMP THREADPRIVATE(horicut_index)

  !
  ! Land cover change
  !
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION(:) :: horip_s_index   !! Horizontal + P short indices
!$OMP THREADPRIVATE(horip_s_index)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION(:) :: horip_m_index   !! Horizontal + P medium indices
!$OMP THREADPRIVATE(horip_m_index)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION(:) :: horip_l_index   !! Horizontal + P long indice
!$OMP THREADPRIVATE(horip_l_index)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION(:) :: horip_ss_index  !! Horizontal + P short indices
!$OMP THREADPRIVATE(horip_ss_index)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION(:) :: horip_mm_index  !! Horizontal + P medium indices
!$OMP THREADPRIVATE(horip_mm_index)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION(:) :: horip_ll_index  !! Horizontal + P long indices
!$OMP THREADPRIVATE(horip_ll_index)

  INTEGER(i_std),SAVE :: itime                 !! time step
!$OMP THREADPRIVATE(itime)
  INTEGER(i_std),SAVE :: hist_id_stomate       !! STOMATE history file ID
!$OMP THREADPRIVATE(hist_id_stomate)
  INTEGER(i_std),SAVE :: hist_id_stomate_IPCC  !! STOMATE history file ID for IPCC output
!$OMP THREADPRIVATE(hist_id_stomate_IPCC)
  INTEGER(i_std),SAVE :: rest_id_stomate       !! STOMATE restart file ID
!$OMP THREADPRIVATE(rest_id_stomate)

  REAL(r_std),PARAMETER :: adapted_crit = 1. - ( 1. / euler ) !! critical value for being adapted (1-1/e) (unitless)
  REAL(r_std),PARAMETER :: regenerate_crit = 1. / euler       !! critical value for being regenerative (1/e) (unitless)


  ! private & public routines

  PUBLIC data

CONTAINS

!! ================================================================================================================================
!! SUBROUTINE 	: data
!!
!>\BRIEF         This routine defines the values of the PFT parameters. It will print the values of the parameters for STOMATE
!!               in the standard outputs of ORCHIDEE. 
!!
!! DESCRIPTION : This routine defines PFT parameters. It initializes the pheno_crit structure by tabulated parameters.\n
!!               Some initializations are done for parameters. The SLA is calculated according *to* Reich et al (1992).\n
!!               Another formulation by Reich et al(1997) could be used for the computation of the SLA.
!!               The geographical coordinates might be used for defining some additional parameters
!!               (e.g. frequency of anthropogenic fires, irrigation of agricultural surfaces, etc.). \n
!!               For the moment, this possibility is not used. \n
!!               The specifc leaf area (SLA) is calculated according Reich et al, 1992 by :
!!               \latexonly
!!               \input{stomate_data_SLA.tex}
!!               \endlatexonly
!!               The sapling (young) biomass for trees and for each compartment of biomass is calculated by :
!!               \latexonly
!!               \input{stomate_data_sapl_tree.tex}
!!               \endlatexonly
!!               The sapling biomass for grasses and for each compartment of biomass is calculated by :
!!               \latexonly
!!               \input{stomate_data_sapl_grass.tex}
!!               \endlatexonly
!!               The critical stem diameter is given by the following formula :
!!               \latexonly
!!               \input{stomate_data_stem_diameter.tex}
!!               \endlatexonly
!!
!! RECENT CHANGE(S): Sonke Zaehle: Reich et al, 1992 find no statistically significant differences 
!!                  between broadleaved and coniferous forests, specifically, the assumption that grasses grow 
!!                  needles is not justified. Replacing the function with the one based on Reich et al. 1997. 
!!                  Given that sla=100cm2/gDW at 9 months, sla is:
!!                  sla=exp(5.615-0.46*ln(leaflon in months)) 
!!                   \latexonly
!!                   \input{stomate_data_SLA_Reich_97.tex}
!!                   \endlatexonly
!!
!! MAIN OUTPUT VARIABLE(S): 
!!
!! REFERENCE(S) :
!! - Reich PB, Walters MB, Ellsworth DS, (1992), Leaf life-span in relation to leaf, plant and 
!! stand characteristics among diverse ecosystems. Ecological Monographs, Vol 62, pp 365-392.
!! - Reich PB, Walters MB, Ellsworth DS (1997) From tropics to tundra: global convergence in plant 
!!  functioning. Proc Natl Acad Sci USA, 94:13730 13734
!!
!! FLOWCHART    :
!! \n
!_ ================================================================================================================================

  SUBROUTINE data


    INTEGER(i_std)                               :: i,j     !! Index (unitless)

!_ ================================================================================================================================

    ! Initialize local printlev
    printlev_loc=get_printlev('data')
    

    !- pheno_gdd_crit
    pheno_gdd_crit(:,1) = pheno_gdd_crit_c(:)
    pheno_gdd_crit(:,2) = pheno_gdd_crit_b(:)         
    pheno_gdd_crit(:,3) = pheno_gdd_crit_a(:) 
    !
    !- senescence_temp
    senescence_temp(:,1) = senescence_temp_c(:)
    senescence_temp(:,2) = senescence_temp_b(:)
    senescence_temp(:,3) = senescence_temp_a(:)
    ! 
    !-LC 
    LC(:,ileaf) = LC_leaf(:) 
    LC(:,isapabove) = LC_sapabove(:) 
    LC(:,isapbelow) = LC_sapbelow(:) 
    LC(:,iheartabove) = LC_heartabove(:) 
    LC(:,iheartbelow) = LC_heartbelow(:) 
    LC(:,iroot) = LC_root(:) 
    LC(:,ifruit) = LC_fruit(:) 
    LC(:,icarbres) = LC_carbres(:) 
    LC(:,ilabile) = LC_labile(:) 

    IF ( printlev_loc >= 2 ) WRITE(numout,*) 'data: PFT characteristics'

    DO j = 2,nvm ! Loop over # PFTS 

       IF ( printlev_loc >= 2 ) WRITE(numout,'(a,i3,a,a)') '    > PFT#',j,': ', PFT_name(j)

       !
       ! 1 tree? (true/false)
       !
       IF ( printlev_loc >= 2 ) WRITE(numout,*) '       tree: (::is_tree) ', is_tree(j)

       !
       ! 2 flamability (0-1, unitless)
       !

       IF ( printlev_loc >= 2 ) WRITE(numout,*) '       litter flamability (::flam) :', flam(j)

       !
       ! 3 fire resistance (unitless)
       !

       IF ( printlev_loc >= 2 ) WRITE(numout,*) '       fire resistance (::resist) :', resist(j)

       !
       ! 4 specific leaf area per mass carbon = 2 * sla / dry mass (m^2.gC^{-1})
       !

       ! S. Zaehle: Reich et al, 1992 find no statistically significant differences between broadleaved and coniferous
       ! forests, specifically, the assumption that grasses grow needles is not justified. Replacing the function
       ! with the one based on Reich et al. 1997. Given that sla=100cm2/gDW at 9 months, sla is:
       ! sla=exp(5.615-0.46*ln(leaflon in months))

       ! Oct 2010 : sla values are prescribed by values given by N.Viovy 

       ! includes conversion from 
       !!       sla(j) = 2. * 1e-4 * EXP(5.615 - 0.46 * log(12./(longevity_eff_leaf(j)/365)))
       !!\latexonly
       !!\input{stomate_data_SLA.tex}
       !!\endlatexonly
!       IF ( leaf_tab(j) .EQ. 2 ) THEN
!
!          ! needle leaved tree
!          sla(j) = 2. * ( 10. ** ( 2.29 - 0.4 * LOG10(12./(longevity_eff_leaf(j)/365)) ) ) *1e-4
!
!       ELSE
!
!          ! broad leaved tree or grass (Reich et al 1992)
!          sla(j) = 2. * ( 10. ** ( 2.41 - 0.38 * LOG10(12./(longevity_eff_leaf(j)/365)) ) ) *1e-4
!
!       ENDIF

!!!$      IF ( leaf_tab(j) .EQ. 1 ) THEN
!!!$
!!!$        ! broad leaved tree
!!!$
!!!$        sla(j) = 2. * ( 10. ** ( 2.41 - 0.38 * LOG10(12./(longevity_eff_leaf(j)/365)) ) ) *1e-4
!!!$
!!!$      ELSE
!!!$
!!!$        ! needle leaved or grass (Reich et al 1992)
!!!$
!!!$        sla(j) = 2. * ( 10. ** ( 2.29 - 0.4 * LOG10(12./(longevity_eff_leaf(j)/365)) ) ) *1e-4
!!!$
!!!$      ENDIF
!!!$
!!!$      IF ( ( leaf_tab(j) .EQ. 2 ) .AND. ( pheno_type_tab(j) .EQ. 2 ) ) THEN
!!!$
!!!$        ! summergreen needle leaf
!!!$
!!!$        sla(j) = 1.25 * sla(j)
!!!$
!!!$      ENDIF

!!$       IF ( printlev_loc >= 2 ) WRITE(numout,*) '       specific leaf area (m**2/gC) (::sla):', sla(j), 12./(longevity_eff_leaf(j)/365)


       !+++CHECK+++
       ! May be needed for the DGVM but no longer needed for stomate_prescribe. There
       ! is a separate subroutine make_sapling that can be called whenever the initial
       ! biomass of a sapling needs to be calculated. Estimating sapling biomass 
       ! on-the-fly rather than once at the start of a simulation has the advantage that 
       ! allocation factors can be adapted to the site conditions.
       !+++++++++++

       !
       ! 6 migration speed (m/year)
       !

       IF ( is_tree(j) ) THEN

          migrate(j) = migrate_tree

       ELSE

          ! can be any value as grasses are, per *definition*, everywhere (big leaf).
          migrate(j) = migrate_grass

       ENDIF !( is_tree(j) )

       IF ( printlev_loc >= 2 ) WRITE(numout,*) '       migration speed (m/year): (::migrate(j))', migrate(j)

       !
       ! 7 critical stem diameter: beyond this diameter, the crown area no longer
       !     increases (m)
       !

       IF ( is_tree(j) ) THEN

          !!\latexonly
          !!\input{stomate_data_stem_diameter.tex}
          !!\endlatexonly          
          cn_sapl(j) = cn_sapl_init !crown of individual tree, first year
         
       ELSE

          cn_sapl(j)=1

       ENDIF !( is_tree(j) )
   
       
       !
       ! 8 Coldest tolerable temperature (K)
       !

       IF ( ABS( tmin_crit(j) - undef ) .GT. min_stomate ) THEN
          tmin_crit(j) = tmin_crit(j) + ZeroCelsius
       ELSE
          tmin_crit(j) = undef
       ENDIF 

       IF ( printlev_loc >= 2 ) &
            WRITE(numout,*) '       coldest tolerable temperature (K): (::tmin_crit(j))', tmin_crit(j)

       !
       ! 9 Maximum temperature of the coldest month: need to be below this temperature
       !      for a certain time to regrow leaves next spring *(vernalization)* (K)
       !

       IF ( ABS ( tcm_crit(j) - undef ) .GT. min_stomate ) THEN
          tcm_crit(j) = tcm_crit(j) + ZeroCelsius
       ELSE
          tcm_crit(j) = undef
       ENDIF

       IF ( printlev_loc >= 2 ) &
            WRITE(numout,*) '       vernalization temperature (K): (::tcm_crit(j))', tcm_crit(j)

       !
       ! 10 critical values for phenology
       !

       ! 10.1 model used

       IF ( printlev_loc >= 2 ) &
            WRITE(numout,*) '       phenology model used: (::pheno_model(j)) ',pheno_model(j)

       ! 10.2 growing degree days. What kind of gdd is meant (i.e. threshold 0 or -5 deg C
       !        or whatever), depends on how this is used in stomate_phenology.


       IF ( ( printlev_loc >= 2 ) .AND. ( ALL(pheno_gdd_crit(j,:) .NE. undef) ) ) THEN
          WRITE(numout,*) '         critical GDD is a function of long term T (C): (::gdd)'
          WRITE(numout,*) '          ',pheno_gdd_crit(j,1), &
               ' + T *',pheno_gdd_crit(j,2), &
               ' + T^2 *',pheno_gdd_crit(j,3)
       ENDIF

       ! consistency check

       IF ( ( ( pheno_model(j) .EQ. 'moigdd' ) .OR. &
            ( pheno_model(j) .EQ. 'humgdd' )       ) .AND. &
            ( ANY(pheno_gdd_crit(j,:) .EQ. undef) )                      ) THEN
          CALL ipslerr_p(3,'stomate_data','problem with phenology parameters, critical GDD. (::pheno_model)','','')
       ENDIF

       ! 10.3 number of growing days

       IF ( ( printlev_loc >= 2 ) .AND. ( ngd_crit(j) .NE. undef ) ) &
            WRITE(numout,*) '         critical NGD: (::ngd_crit(j))', ngd_crit(j)

       ! 10.4 critical temperature for ncd vs. gdd function in phenology (C)

       IF ( ( printlev_loc >= 2 ) .AND. ( ncdgdd_temp(j) .NE. undef ) ) &
            WRITE(numout,*) '         critical temperature for NCD vs. GDD (C): (::ncdgdd_temp(j))', &
            ncdgdd_temp(j)

       ! 10.5 humidity fractions (0-1, unitless)

       IF ( ( printlev_loc >= 2 ) .AND. ( hum_frac(j) .NE. undef ) ) &
            WRITE(numout,*) '         critical humidity fraction: (::hum_frac(j))', &
            &  hum_frac(j)


       ! 10.6 minimum time elapsed since moisture minimum (days)

       IF ( ( printlev_loc >= 2 ) .AND. ( hum_min_time(j) .NE. undef ) ) &
            WRITE(numout,*) '         time to wait after moisture min (d): (::hum_min_time(j))', &
        &    hum_min_time(j)

       !
       ! 11 critical values for senescence
       !

       ! 11.1 type of senescence

       IF ( printlev_loc >= 2 ) WRITE(numout,*) '       type of senescence: (::senescence_type(j))',&
            senescence_type(j)

       ! 11.2 critical temperature for senescence (C)

       IF ( ( printlev_loc >= 2 ) .AND. ( ALL(senescence_temp(j,:) .NE. undef) ) ) THEN
          WRITE(numout,*) '         critical temperature for senescence (C) is'
          WRITE(numout,*) '          a function of long term T (C): (::senescence_temp)'
          WRITE(numout,*) '          ',senescence_temp(j,1), &
               ' + T *',senescence_temp(j,2), &
               ' + T^2 *',senescence_temp(j,3)
       ENDIF

       ! consistency check

       IF ( ( ( senescence_type(j) .EQ. 'cold' ) .OR. &
            ( senescence_type(j) .EQ. 'mixed' )      ) .AND. &
            ( ANY(senescence_temp(j,:) .EQ. undef ) )           ) THEN
          CALL ipslerr_p(3,'stomate_data','Problem with senescence parameters, temperature. (::senescence_type)','','')
       ENDIF

       ! 11.3 critical relative moisture availability for senescence

       IF ( ( printlev_loc >= 2 ) .AND. ( senescence_hum(j) .NE. undef ) ) THEN
          WRITE(numout,*)  ' max. critical relative moisture availability for' 
          WRITE(numout,*)  ' senescence: (::senescence_hum(j))',  &
               & senescence_hum(j)
       ENDIF

       ! consistency check

       IF ( ( ( senescence_type(j) .EQ. 'dry' ) .OR. &
            ( senescence_type(j) .EQ. 'mixed' )     ) .AND. &
            ( senescence_hum(j) .EQ. undef )                   ) THEN
          CALL ipslerr_p(3,'stomate_data','Problem with senescence parameters, humidity.(::senescence_type)','','')
       ENDIF

       ! 14.3 relative moisture availability above which there is no moisture-related
       !      senescence (0-1, unitless)

       IF ( ( printlev_loc >= 2 ) .AND. ( nosenescence_hum(j) .NE. undef ) ) THEN
          WRITE(numout,*) '         relative moisture availability above which there is' 
          WRITE(numout,*) '             no moisture-related senescence: (::nosenescence_hum(j))', &
               &  nosenescence_hum(j)
       ENDIF

       ! consistency check

       IF ( ( ( senescence_type(j) .EQ. 'dry' ) .OR. &
            ( senescence_type(j) .EQ. 'mixed' )     ) .AND. &
            ( nosenescence_hum(j) .EQ. undef )                   ) THEN
          CALL ipslerr_p(3,'stomate_data','Problem with senescence parameters, humidity. (::senescence_type)','','')
       ENDIF

       !
       ! 12 sapwood -> heartwood conversion time (days)
       !

       IF ( printlev_loc >= 2 ) &
            WRITE(numout,*) '       sapwood -> heartwood conversion time (d): (::longevity_sap(j))', longevity_sap(j)

       !
       ! 13 fruit lifetime (days)
       !

       IF ( printlev_loc >= 2 ) WRITE(numout,*) '       fruit lifetime (d): (::longevity_fruit(j))', longevity_fruit(j)

       !
       ! 14 length of leaf death (days)
       !      For evergreen trees, this variable determines the lifetime of the leaves.
       !      Note that it is different from the value given in (longevity_eff_leaf/365).
       !

       IF ( printlev_loc >= 2 ) &
            WRITE(numout,*) '       length of leaf death (d): (::leaffall(j))', leaffall(j)


       !
       ! 17 minimum lai, initial (m^2.m^{-2})
       !

       IF ( is_tree(j) ) THEN
          lai_initmin(j) = lai_initmin_tree
       ELSE
          lai_initmin(j) = lai_initmin_grass
       ENDIF !( is_tree(j) )

       IF ( printlev_loc >= 2 ) &
            WRITE(numout,*) '       initial LAI: (::lai_initmin(j))', lai_initmin(j)

       !
       ! 19 maximum LAI (m^2.m^{-2})
       !

       IF ( printlev_loc >= 2 ) &
            WRITE(numout,*) '       critical LAI above which no leaf allocation: (::lai_max(j))', lai_max(j)

       !
       ! 20 fraction of primary leaf and root allocation put into reserve (0-1, unitless)
       !

       IF ( printlev_loc >= 2 ) &
            WRITE(numout,*) '       reserve allocation factor: (::ecureuil(j))', ecureuil(j)


       ! 23 natural ?
       !

       IF ( printlev_loc >= 2 ) &
            WRITE(numout,*) '       Natural: (::natural(j))', natural(j)

       !
       ! 24 Vcmax et Vjmax (umol.m^{-2}.s^{-1}) 
       !

       IF ( printlev_loc >= 2 ) &
            WRITE(numout,*) '       Maximum rate of carboxylation: (::Vcmax_25(j))', vcmax25(j)
       !
       ! 25 constants for photosynthesis temperatures
       !

       IF ( printlev_loc >= 2 ) THEN


          !
          ! 26 Properties
          !
          WRITE(numout,*) '       C4 photosynthesis: (::is_c4(j))', is_c4(j)

       ENDIF

       !
       ! 27 extinction coefficient of the Monsi and Saeki (1953) relationship 
       !
       IF ( printlev_loc >= 2 ) THEN
          WRITE(numout,*) '       extinction coefficient: (::ext_coeff(j))', ext_coeff(j)
       ENDIF

       !
       ! 30 fraction of allocatable biomass which is lost as growth respiration (0-1, unitless)
       !
       IF ( printlev_loc >= 2 ) &
            WRITE(numout,*) '       growth respiration fraction: (::frac_growthresp(j))', frac_growthresp(j)

    ENDDO ! Loop over # PFTS 

    !
    ! 29 time scales for phenology and other processes (in days)
    !

    tau_longterm_max = coeff_tau_longterm * one_year

    IF ( printlev_loc >= 2 ) THEN

       WRITE(numout,*) '   > time scale for ''monthly'' moisture availability (d): (::tau_hum_month)', &
            tau_hum_month
       WRITE(numout,*) '   > time scale for ''weekly'' moisture availability (d): (::tau_hum_week)', &
           tau_hum_week
       WRITE(numout,*) '   > time scale for ''monthly'' 2 meter temperature (d): (::tau_t2m_month)', &
            tau_t2m_month
       WRITE(numout,*) '   > time scale for ''weekly'' 2 meter temperature (d): (::tau_t2m_week)', &
            tau_t2m_week
       WRITE(numout,*) '   > time scale for ''weekly'' GPP (d): (::tau_gpp_week)', &
            tau_gpp_week
       WRITE(numout,*) '   > time scale for ''monthly'' soil temperature (d): (::tau_tsoil_month)', &
            tau_tsoil_month
       WRITE(numout,*) '   > time scale for vigour calculations (y): (::tau_longterm_max / one_year)', &
            tau_longterm_max / one_year

    ENDIF

    !
    ! 30 Maintenance respiration
    !

    maint_resp_slope(:,1) = maint_resp_slope_c(:)               
    maint_resp_slope(:,2) = maint_resp_slope_b(:) 
    maint_resp_slope(:,3) = maint_resp_slope_a(:)
 
    IF (printlev_loc >= 4) WRITE(numout,*) 'Leaving stomate_data'

  END SUBROUTINE data

END MODULE stomate_data
