!!=================================================================================================================================
! MODULE       : constantes_var
!
! CONTACT      : orchidee-help _at_ listes.ipsl.fr
!
! LICENCE      : IPSL (2006)
! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF        constantes_var module contains most constantes like pi, Earth radius, etc...
!!              and all externalized parameters except pft-dependent constants.
!!
!!\n DESCRIPTION: This module contains most constantes and the externalized parameters of ORCHIDEE which 
!!                are not pft-dependent.\n
!!                In this module, you can set the flag diag_qsat in order to detect the pixel where the
!!                temperature is out of range (look qsatcalc and dev_qsatcalc in qsat_moisture.f90).\n
!!                The Earth radius is approximated by the Equatorial radius.The Earth's equatorial radius a,
!!                or semi-major axis, is the distance from its center to the equator and equals 6,378.1370 km.
!!                The equatorial radius is often used to compare Earth with other planets.\n
!!                The meridional mean is well approximated by the semicubic mean of the two axe yielding 
!!                6367.4491 km or less accurately by the quadratic mean of the two axes about 6,367.454 km
!!                or even just the mean of the two axes about 6,367.445 km.\n
!!                This module is already USE in module constantes. Therefor no need to USE it seperatly except
!!                if the subroutines in module constantes are not needed.\n
!!                
!! RECENT CHANGE(S):
!!
!! REFERENCE(S)	: 
!! - Louis, Jean-Francois (1979), A parametric model of vertical eddy fluxes in the atmosphere. 
!! Boundary Layer Meteorology, 187-202.\n
!!
!! SVN          :
!! $HeadURL: svn://forge.ipsl.fr/orchidee/trunk/ORCHIDEE/src_parameters/constantes_var.f90 $
!! $Date: 2024-05-14 12:29:33 +0200 (mar. 14 mai 2024) $
!! $Revision: 8556 $
!! \n
!_ ================================================================================================================================

MODULE constantes_var

  USE defprec

  IMPLICIT NONE
!-

                         !-----------------------!
                         !  ORCHIDEE CONSTANTS   !
                         !-----------------------!

  !
  ! FLAGS 
  !

  INTEGER, SAVE       :: energy_control                       !! Flag that automatically controls several other flags related to 
                                                              !! multi-layering (1/2/3).
                                                              !! 1 - DEFAULT used the enerbil module in combination with the
                                                              !! hydraulic architecture ((ok_hydrol_arch and ok_gs_feedback
                                                              !! true, while ok_mleb and ok_impose_canopy_structure are set
                                                              !! to false
                                                              !! 2 - option to use enerbil module and original water stress
                                                              !! (not hydraulic architecture)
                                                              !! 3 - single layer: ok_hydrol_arch, ok_gs_feedback, ok_impose_canopy_structure 
                                                              !! and ok_mleb all true, but the energy budget is only calculated for a single layer 
                                                              !! (jnlvls=1,jnlvls_under=0,jnlvls_canopy=1,jnlvls_over=0).
                                                              !! 4 - multi-layer: ok_hydrol_arch, ok_gs_feedback,  ok_impose_canopy_structure,
                                                              !! ok_mleb all true, and the energy budget is calculated for multiple layers 
                                                              !! (jnlvls=29,jnlvls_under=10,jnlvls_canopy=10,jnlvls_over=9).    
                                                              !! 5 - user specific: user specific settings for these controls 
                                                              !! and layers as defined in the run.def by the user.
!$OMP THREADPRIVATE(energy_control)
  LOGICAL, SAVE       :: ok_hydrol_arch                       !! Flag that activates the hydraulic architecture routine (true/false) 
                                                              !! The trunk version of ORCHIDEE (false) uses soil water as a 
                                                              !! proxy for water stress and applies the stress to Vcmax.
                                                              !! When set to true the hydraulic architecture of the vegetation
                                                              !! is accounted for to calculate the amount of water that 
                                                              !! can be transported through the plant given the soil and leaf
                                                              !! potential and the conductivities of the roots, wood and 
                                                              !! leaves. Water supply through the plant is compared against
                                                              !! the atmospheric demand for water. If the supply is smaller
                                                              !! then the demand, the plant experiences water stress and the
                                                              !! stomata will be closed (water stress is now on gs rather 
                                                              !! than Vcmax). Note that whether stomatal regulation is used or 
                                                              !! not is controled by a separate flag: ok_gs_feedback.
!$OMP THREADPRIVATE(ok_hydrol_arch)


  LOGICAL       :: is_tuzet_hydrol_arch = .FALSE.             !! Flag that activates the new version of hydraulic architecture
                                                              !! based on the model developed by Tuzet et al. (2017). If the flag
                                                              !! is set to FALSE, the previous version of the hydraulic architecture
                                                              !! will be running.
!$OMP THREADPRIVATE(is_tuzet_hydrol_arch)

  LOGICAL       :: ok_tuzet_hydrol_arch_muff = .FALSE.        !! Flag that activates the radial resolution of Richard's equation
                                                              !! around the roots. This resolution is only made when the
                                                              !! hydraulic architecture is running. If set to .FALSE., the root 
                                                              !! absorption is modelled by a classical resistance term.
!$OMP THREADPRIVATE(ok_tuzet_hydrol_arch_muff)

  LOGICAL       :: ok_tuzet_hydrol_arch_dyn_resist = .FALSE.  !! Flag that activates the dynamic resistance equations used by
                                                              !! Yao et al. (2021) (Under review at the time this is developed)
                                                              !! instead of having parameterised hydraulic resistances as in Tuzet 
                                                              !! et al. (2017), the resistances now depend on the water potential of the stage.
!$OMP THREADPRIVATE(ok_tuzet_hydrol_arch_dyn_resist)

  LOGICAL       :: split_soil_properties = .FALSE.            !! Flag that permits to distinguish properties between the superficial
                                                              !! soil horizon and the inferior one in the new hydraulic architecture.
!$OMP THREADPRIVATE(split_soil_properties)

  LOGICAL       :: ok_tuzet_hydrol_arch_storage = .FALSE.     !! Flag that activates the resolution of water storage inside the
                                                              !! vegetation as in Tuzet et al. (2017).
!$OMP THREADPRIVATE(ok_tuzet_hydrol_arch_storage)


  LOGICAL, SAVE       :: ok_gs_feedback                       !! Flag that activates water stress on stomata (true/false)
                                                              !! This flag is for debugging only! It allows developers
                                                              !! to calculate GPP without any water stress. If the model is 
                                                              !! used in production mode and ok_hydrol_arch is true this 
                                                              !! flag should be true as well.
!$OMP THREADPRIVATE(ok_gs_feedback)   
  LOGICAL, SAVE       :: ok_mleb                              !! Flag that activates the multilayer energy budget (true/false)
                                                              !! The model uses 10 (default) canopy layers to calculate
                                                              !! the albedo, transmittance, absorbance and GPP. These canopy
                                                              !! layers can be combined with 10 (default) layers below and 
                                                              !! 10 layers above the canopy to calculate the energy budget
                                                              !! (ok_mleb=y). If set to no, this flag will make the model
                                                              !! use 10 layers for the canopy albedo, transmittance, 
                                                              !! absorbance and GPP and just a single layer for the energy 
                                                              !! budget. Be aware that if you wish to run with hydraulic 
                                                              !! architechture ok_mleb needs to be se to true as well. Furthermore
                                                              !! if you  wish to run with the original energy scheme (enerbil),
                                                              !! set the layers for mleb to 1.
!$OMP THREADPRIVATE(ok_mleb)
  LOGICAL, SAVE       :: ok_impose_can_structure              !! This flag is for debugging only! It allows developers
                                                              !! to use a prescribed canopy structure rather then the 
                                                              !! structure calculate by ORCHIDEE. The flag activates the 
                                                              !! sections of code which directly link the energy budget 
                                                              !! scheme to the the size and LAI profile of the canopy for the 
                                                              !! respective PFT and age class that is calculated in stomate, 
                                                              !! for the albedo. If set to TRUE and the multi-layer budget 
                                                              !! is activated the model takes LAI profile information and 
                                                              !! canopy level heights from the run.def. If set to FALSE, and 
                                                              !! the multi-layer energy budget is used the profile 
                                                              !! information and canopy levels heights comes from the
                                                              !! PGap-based processes for calculation of stand profile 
                                                              !! information in stomate.
!$OMP THREADPRIVATE(ok_impose_can_structure)
  LOGICAL, SAVE       :: ok_mleb_history_file                 !! Flag that controls the writing of an output file with the 
                                                              !! multi-layer energy simulations (true/false). Note that this
                                                              !! is a large file and writing it slows down the code.
!$OMP THREADPRIVATE(ok_mleb_history_file)
  LOGICAL, SAVE :: ok_read_fm_map = .FALSE.                   !! A logical flag determining if we read
                                                              !! in the forest management strategy from a map.
                                                              !! This should be .TRUE. for all applications
                                                              !! except for debugging and pixel-level simulations.
!$OMP THREADPRIVATE(ok_read_fm_map)
  LOGICAL, SAVE :: ok_read_sp_clearcut_map = .FALSE.          !! A logical flag determining if we read
                                                              !! in a map determining for each pixel and PFT,
                                                              !! whether a clearcut during spinup will happen.
!$OMP THREADPRIVATE(ok_read_sp_clearcut_map)
  LOGICAL, SAVE :: ok_change_species = .FALSE.                !! A logical flag determining if we change
                                                              !! species after a clearcut
!$OMP THREADPRIVATE(ok_change_species)
  LOGICAL, SAVE :: ok_read_species_change_map = .FALSE.       !! A logical flag determining if we read
                                                              !! in a map which changes species after a clearcut
                                                              !! To be used with ok_change_species = .TRUE. or
                                                              !! set species_change_force (see below)
!$OMP THREADPRIVATE(ok_read_species_change_map)
  LOGICAL, SAVE :: ok_read_desired_fm_map = .FALSE.           !! A logical flag determining if we read
                                                              !! in the desired forest management strategy from a map.
                                                              !! To be used with ok_change_species = .TRUE. or set
                                                              !! fm_change_force (see below)
!$OMP THREADPRIVATE(ok_read_desired_fm_map)
  LOGICAL, SAVE :: ok_litter_raking = .FALSE.                 !! If TRUE, this flag will simulate litter raking in
	                                                      !! in grid squares.  This has the effect of moving litter
	                                                      !! once a year from forest PFTs to agricultural PFTs, if they
	                                                      !! are present on this pixel.  If TRUE, you must also provide
	                                                      !! a map with the litter demand so we know how much litter
	                                                      !! to remove for each pixel. Litter raking is a historical
                                                              !! land use so you reconstrauctions to use this option.
!$OMP THREADPRIVATE(ok_litter_raking)
 
  LOGICAL       :: ok_dimensional_product_use = .TRUE.        !! Once the wood is harvested it ends up in wood product pools
                                                              !! Two options were coded: (1) the product use and thus 
                                                              !! its on longevity depends on the dimensions of the harvest
                                                              !! (2) the dimensions are ignored and the wood is used according
                                                              !! to fixed ratios.
 
!$OMP THREADPRIVATE(ok_dimensional_product_use)

  LOGICAL       :: ok_constant_mortality                      !! Use constant mortality or calculate mortality 
                                                              !! as a function of last years´s NPP
!$OMP THREADPRIVATE(ok_constant_mortality)

  LOGICAL, SAVE :: ok_c13                                     !! Activate carbon isotope concetration of biomass
!$OMP THREADPRIVATE(ok_c13)

  LOGICAL, SAVE :: ok_windthrow = .FALSE.                     !! Activate the wind throw module. Trees will be killed
                                                              !! if the wind speed exceeds the critical wind speed of
                                                              !! of the forest (PFT).                             
!$OMP THREADPRIVATE(ok_windthrow)

  LOGICAL, SAVE :: ok_pest = .FALSE.                          !! Activate pest  module. Trees will be killed by bark beetle
 !$OMP THREADPRIVATE(ok_pest)

  LOGICAL, SAVE :: ok_bare_soil_new                           !! Choose between the two options to calculate the bare soil.
                                                              !! False = classic view: gaps within a canopy should be treated
                                                              !! as bare soil. True = ecological view: gaps within a canopy are
                                                              !! part of the ecosystem and should be treated as such.
                                                              !! The parameter is set in control.f90 as part of the ENERGY_CONTROL 
                                                              !! flag.
!$OMP THREADPRIVATE(ok_bare_soil_new)

  LOGICAL, SAVE :: ok_snow_albedo_clm3 = .FALSE.              !! Choose between two options to calculate the snow albedo. If 
                                                              !! set to .TRUE. the snow albedo is calculated with CLM3 scheme having the diffuse
                                                              !! and direct light separately and decays with age. If set to .FALSE.
                                                              !! a scheme from Chalita and Le Treut is used.
!$OMP THREADPRIVATE(ok_snow_albedo_clm3)

  CHARACTER(LEN=30), SAVE :: maint_resp_control               !! Choose between two options to calculate the maint respiration. If
                                                              !! set to 'nitrogen' the plant nitrogen pools and temperature will 
                                                              !! drive maint_resp. If set to 'c/n', maint_resp will be adjusted for 
                                                              !! the expected c/n ratio of the plant biomass. Also the temperature 
                                                              !! control itself is calculated according to Krinner et al 2005 instead
                                                              !! of Sitch et al 2003.

!$OMP THREADPRIVATE(maint_resp_control) 

  LOGICAL, SAVE :: ok_force_pheno = .TRUE.                    !! Temperature phenology is very predictable and happens every year but
                                                              !! moisture-driven phenology is less stable because the conditions
                                                              !! may not be satisfied during 12 months resulting in a year without
                                                              !! a canopy. If the model has passed the average budbreak day by a
                                                              !! prescribed PFT-specific offset (::force_pheno), budbreak will be forced
                                                              !! (given the condition that buds are avaialble)
!$OMP THREADPRIVATE(ok_force_pheno)

  LOGICAL, SAVE :: nc_restart_compression !! activate netcdf restart compression
!$OMP THREADPRIVATE(nc_restart_compression)
  CHARACTER(LEN=80), SAVE :: which_rt        !! Flag for radiative transfer: 'Iterative', 'Matrixial'
!$OMP THREADPRIVATE(which_rt)
  LOGICAL :: river_routing                   !! activate river routing
!$OMP THREADPRIVATE(river_routing)
  CHARACTER(LEN=80), SAVE :: routing_method  !! 'standard', 'highres', 'simple' or 'none': Character string used to switch between routing modules
!$OMP THREADPRIVATE(routing_method) 
  LOGICAL, SAVE :: ok_nudge_mc  !! Activate nudging of soil moisture 
!$OMP THREADPRIVATE(ok_nudge_mc)
  LOGICAL, SAVE :: ok_nudge_snow!! Activate nudging of snow variables
!$OMP THREADPRIVATE(ok_nudge_snow)
  LOGICAL, SAVE :: nudge_interpol_with_xios  !! Activate reading and interpolation with XIOS for nudging fields
!$OMP THREADPRIVATE(nudge_interpol_with_xios)
  LOGICAL :: do_floodplains     !! activate flood plains
!$OMP THREADPRIVATE(do_floodplains)
  LOGICAL :: do_irrigation      !! activate computation of irrigation flux
!$OMP THREADPRIVATE(do_irrigation)
  LOGICAL :: ok_sechiba         !! activate physic of the model
!$OMP THREADPRIVATE(ok_sechiba)
  LOGICAL :: ok_stomate         !! activate carbon cycle
!$OMP THREADPRIVATE(ok_stomate)
  LOGICAL :: ok_ncycle          !! activate nitrogen cycle
!$OMP THREADPRIVATE(ok_ncycle)
  LOGICAL :: impose_cn          !! impose the CN ratio of leaves
!$OMP THREADPRIVATE(impose_cn)
  LOGICAL :: reset_impose_cn    !! reset the CN ratio of leaves
!$OMP THREADPRIVATE(reset_impose_cn)
  LOGICAL :: read_cn          !! read the CN ratio of leaves
!$OMP THREADPRIVATE(read_cn)
  LOGICAL :: ok_dgvm            !! activate dynamic vegetation
!$OMP THREADPRIVATE(ok_dgvm)
  LOGICAL :: do_wood_harvest    !! activate wood harvest
!$OMP THREADPRIVATE(do_wood_harvest)
  LOGICAL :: ok_pheno           !! activate the calculation of lai using stomate rather than a prescription
!$OMP THREADPRIVATE(ok_pheno)
  LOGICAL :: ok_bvoc            !! activate biogenic volatile organic coumpounds
!$OMP THREADPRIVATE(ok_bvoc)
  LOGICAL :: ok_leafage         !! activate leafage
!$OMP THREADPRIVATE(ok_leafage)
  LOGICAL :: ok_snags           !! create snag litter pool. When set to yes, snag
                                !! pool will be accounted for, with snag decay and snag fall into woody pool
!$OMP THREADPRIVATE(ok_snags)
  LOGICAL :: ok_radcanopy       !! use canopy radiative transfer model
!$OMP THREADPRIVATE(ok_radcanopy)
  LOGICAL :: ok_multilayer      !! use canopy radiative transfer model with multi-layers
!$OMP THREADPRIVATE(ok_multilayer)
  LOGICAL :: ok_pulse_NOx       !! calculate NOx emissions with pulse
!$OMP THREADPRIVATE(ok_pulse_NOx)
  LOGICAL :: ok_bbgfertil_NOx   !! calculate NOx emissions with bbg fertilizing effect
!$OMP THREADPRIVATE(ok_bbgfertil_NOx)
  LOGICAL :: ok_cropsfertil_NOx !! calculate NOx emissions with fertilizers use
!$OMP THREADPRIVATE(ok_cropsfertil_NOx)

  LOGICAL :: ok_co2bvoc_poss    !! CO2 inhibition on isoprene activated following Possell et al. (2005) model
!$OMP THREADPRIVATE(ok_co2bvoc_poss)
  LOGICAL :: ok_co2bvoc_wilk    !! CO2 inhibition on isoprene activated following Wilkinson et al. (2006) model
!$OMP THREADPRIVATE(ok_co2bvoc_wilk)
  LOGICAL :: use_fluxnet
!$OMP THREADPRIVATE(use_fluxnet)
  LOGICAL, SAVE :: OFF_LINE_MODE = .FALSE.  !! ORCHIDEE detects if it is coupled with a GCM or 
                                            !! just use with one driver in OFF-LINE. (true/false)
!$OMP THREADPRIVATE(OFF_LINE_MODE)  
  CHARACTER(LEN=80), SAVE     :: restname_in       = 'NONE'                 !! Input Restart files name for Sechiba component  
!$OMP THREADPRIVATE(restname_in)
  CHARACTER(LEN=80), SAVE :: restname_out = 'sechiba_rest_out.nc'  !! Output Restart files name for Sechiba component
!$OMP THREADPRIVATE(restname_out)
  CHARACTER(LEN=80), SAVE :: stom_restname_in = 'NONE'        !! Input Restart files name for Stomate component
!$OMP THREADPRIVATE(stom_restname_in)
  CHARACTER(LEN=80), SAVE :: stom_restname_out = 'stomate_rest_out.nc'  !! Output Restart files name for Stomate component
!$OMP THREADPRIVATE(stom_restname_out)
  INTEGER, SAVE :: printlev=2       !! Standard level for text output [0, 1, 2, 3]
!$OMP THREADPRIVATE(printlev)

  !
  ! HACKS
  !
  LOGICAL, SAVE :: hack_enerbil_hydrol = .FALSE.!! For debugging only! Flag to skip a particular block of code in enerbil.f90 which results in
                                                !! incorrect results for large scale simulations.
!$OMP THREADPRIVATE(hack_enerbil_hydrol)
  LOGICAL, SAVE :: hack_pgap = .FALSE.          !! For debugging only! If =.TRUE. the model uses Lambert Beer to calculate veget instead of Pgap.
!$OMP THREADPRIVATE(hack_pgap)
  REAL(r_std), SAVE :: hack_vessel_loss = -9999 !! Impose a constant vessel_loss in hydraulic_rachitecture when ok_vessel_mortality is TRUE. When set outside the range from 0 to 1 it will not be used.
!$OMP THREADPRIVATE(hack_vessel_loss)
  LOGICAL, SAVE :: hack_e_frac = .FALSE.        !! By pass the use of root length in the calculation of psi_soilroot.
!$OMP THREADPRIVATE(hack_e_frac)
  LOGICAL, SAVE :: hack_veget_max_new = .FALSE. !! Prescribe veget_max from run.def rather than a PFT map. Useful feature to debug simplified LCCs. 
                                                !! See slowproc.f90 for details on this feature.
!$OMP THREADPRIVATE(hack_veget_max_new) 
  LOGICAL, SAVE :: hack_kb_m1 = .FALSE.         !! Use higher threshold for small lai in the calculation of kB_m1 in condveg
!$OMP THREADPRIVATE(hack_kb_m1)
  LOGICAL, SAVE :: trusting_hack_age_class = .FALSE.!! Skip the calculation of the biomass redistribution of the age classes
!$OMP THREADPRIVATE(trusting_hack_age_class)
  LOGICAL, SAVE :: test_dynamic_max_height = .FALSE. !! Use a dynamic parameter for pipe_tune2.
!$OMP THREADPRIVATE(test_dynamic_max_height)  
  LOGICAL, SAVE ::  test_Q10_on_plant_uptake = .FALSE. !! Use a Q10 function to control the plant uptake.
!$OMP THREADPRIVATE(test_Q10_on_plant_uptake) 
  !
  ! PARAMETERS
  !
  REAL(r_std), SAVE :: min_n = 0.0001          !! Minimum allowable n_mineralisation when truncating som_input_total(:,initrogen) in stomate_litter. 
!$OMP THREADPRIVATE(min_n)
  REAL(r_std), SAVE :: max_cn = 250            !! Maximum allowable ratio of som_input_total(:,icarbon) to som_input_total(:,initrogen).
!$OMP THREADPRIVATE(max_cn)
  REAL(r_std), SAVE :: snc = 0.004             !! Structural nitrogen  [gN g-1C] based on C:N of dead wood (White et al., 2000)
                                               !! assuming carbon dry matter ratio of 0.5  
!$OMP THREADPRIVATE(snc) 
  REAL(r_std), SAVE :: slope_ra = 40           !! Reduction factor to make resp_maint less temperature sensitive
!$OMP THREADPRIVATE(slope_ra)

  !
  ! TIME
  !
  INTEGER(i_std), PARAMETER  :: spring_days_max = 40  !! Maximum number of days during which we watch for possible spring frost damage

  !
  ! SPECIAL VALUES 
  !
  INTEGER(i_std), PARAMETER :: undef_int = 999999999     !! undef integer for integer arrays (unitless)
  !-
  REAL(r_std), SAVE :: val_exp = 999999.                 !! Specific value if no restart value  (unitless)
!$OMP THREADPRIVATE(val_exp)
  REAL(r_std), PARAMETER :: undef = -9999.               !! Special value for stomate (unitless)
  
  REAL(r_std), PARAMETER :: min_sechiba = 1.E-8_r_std    !! Epsilon to detect a near zero floating point (unitless)
  REAL(r_std), PARAMETER :: undef_sechiba = 1.E+20_r_std !! The undef value used in SECHIBA (unitless)
  
  REAL(r_std), PARAMETER :: min_stomate = 1.E-8_r_std    !! Epsilon to detect a near zero floating point (unitless)
  REAL(r_std), PARAMETER :: large_value = 1.E33_r_std    !! some large value (for stomate) (unitless)


  REAL(r_std), SAVE :: alpha_nudge_mc                    !! Nudging constant for soil moisture 
!$OMP THREADPRIVATE(alpha_nudge_mc)
  REAL(r_std), SAVE :: alpha_nudge_snow                  !! Nudging constant for snow variables
!$OMP THREADPRIVATE(alpha_nudge_snow)

  !
  !  DIMENSIONING AND INDICES PARAMETERS  
  ! 
  INTEGER(i_std), PARAMETER :: istart = 1        !! Index to store values at the start
  INTEGER(i_std), PARAMETER :: iend = 2          !! Index to store values at the end
  INTEGER(i_std), PARAMETER :: ibare_sechiba = 1 !! Index for bare soil in Sechiba (unitless)
  INTEGER(i_std), PARAMETER :: ivis = 1          !! index for albedo in visible range (unitless)
  INTEGER(i_std), PARAMETER :: inir = 2          !! index for albeod i near-infrared range (unitless) 
  INTEGER(i_std), PARAMETER :: n_spectralbands=2 !! number of spectral bands
  INTEGER(i_std), PARAMETER :: nnobio = 1        !! Number of other surface types: land ice (lakes,cities, ...) (unitless)
  INTEGER(i_std), PARAMETER :: iice = 1          !! Index for land ice (see nnobio) (unitless)
  !-
  !! Soil
  INTEGER(i_std), PARAMETER :: classnb = 9       !! Levels of soil colour classification (unitless)
  !-
  INTEGER(i_std), PARAMETER :: nleafages = 4     !! leaf age discretisation ( 1 = no discretisation )(unitless)
  !-
  !! litter fractions: indices (unitless)
  INTEGER(i_std), PARAMETER :: ileaf = 1         !! Index for leaf compartment (unitless)
  INTEGER(i_std), PARAMETER :: isapabove = 2     !! Index for sapwood above compartment (unitless)
  INTEGER(i_std), PARAMETER :: isapbelow = 3     !! Index for sapwood below compartment (unitless)
  INTEGER(i_std), PARAMETER :: iheartabove = 4   !! Index for heartwood above compartment (unitless)
  INTEGER(i_std), PARAMETER :: iheartbelow = 5   !! Index for heartwood below compartment (unitless)
  INTEGER(i_std), PARAMETER :: iroot = 6         !! Index for roots compartment (unitless)
  INTEGER(i_std), PARAMETER :: ifruit = 7        !! Index for fruits compartment (unitless)
  INTEGER(i_std), PARAMETER :: icarbres = 8      !! Index for reserve compartment (unitless)
  INTEGER(i_std), PARAMETER :: ilabile = 9       !! Index for reserve compartment (unitless) 
  INTEGER(i_std), PARAMETER :: nparts = 9        !! Number of biomass compartments (unitless)
  !-
  !! indices for assimilation parameters 
  INTEGER(i_std), PARAMETER :: ivcmax = 1        !! Index for vcmax (assimilation parameters) (unitless)
  INTEGER(i_std), PARAMETER :: inue = 2          !! Index for nue (assimilation parameters) (unitless)
  INTEGER(i_std), PARAMETER :: ileafN = 3        !! Index for leaf N (assimilation parameters) (unitless)
  INTEGER(i_std), PARAMETER :: npco2 = 3         !! Number of assimilation parameters (unitless)

  !-
  !! trees and litter: indices for the parts of heart-
  !! and sapwood above and below the ground 
  INTEGER(i_std), PARAMETER :: iabove = 1       !! Index for above part (unitless)
  INTEGER(i_std), PARAMETER :: ibelow = 2       !! Index for below part (unitless)
  INTEGER(i_std), PARAMETER :: nlevs = 2        !! Number of levels for trees and litter (unitless)
  !-
  !! litter: indices for metabolic and structural part
  INTEGER(i_std), PARAMETER :: imetabolic = 1   !! Index for metabolic litter (unitless)
  INTEGER(i_std), PARAMETER :: istructural = 2  !! Index for structural litter (unitless)
  INTEGER(i_std), PARAMETER :: iwoody = 3       !! Index for woody litter (unitless)
  INTEGER(i_std), PARAMETER :: isnag = 4        !! Index for snag litter (unitless)
  INTEGER(i_std), PARAMETER :: nlitt = 4        !! Number of levels for litter compartments (unitless)
  !-
  !! carbon pools: indices
  INTEGER(i_std), PARAMETER :: iactive = 1      !! Index for active carbon pool (unitless)
  INTEGER(i_std), PARAMETER :: islow = 2        !! Index for slow carbon pool (unitless)
  INTEGER(i_std), PARAMETER :: ipassive = 3     !! Index for passive carbon pool (unitless)
  INTEGER(i_std), PARAMETER :: isurface = 4     !! Index for passive carbon pool (unitless)
  INTEGER(i_std), PARAMETER :: ncarb = 4        !! Number of soil carbon pools (unitless)
  !-
  !! For isotopes and nitrogen
  INTEGER(i_std), PARAMETER :: nelements = 2    !! Number of elements considered
  INTEGER(i_std), PARAMETER :: icarbon = 1      !! Index for carbon 
  INTEGER(i_std), PARAMETER :: initrogen = 2    !! Index for nitrogen 
  !! N-cycle : indices
  INTEGER(i_std), PARAMETER :: iammonium = 1    !! Index for Ammonium 
  INTEGER(i_std), PARAMETER :: initrate  = 2    !! Index for Nitrate
  INTEGER(i_std), PARAMETER :: inox      = 3    !! Index for NOX
  INTEGER(i_std), PARAMETER :: initrous  = 4    !! Index for N2O
  INTEGER(i_std), PARAMETER :: idinitro  = 5    !! Index for N2
  INTEGER(i_std), PARAMETER :: nnspec    = 5    !! Number of N-species considered
  INTEGER(i_std), PARAMETER :: nionspec  = 2    !! Number of ions form considered (ammonium, nitrate)
  !! For N-deposition
  INTEGER(i_std), PARAMETER :: iatm_ammo = 1    !! Index for N input from Ammonium N atmospheric deposition
  INTEGER(i_std), PARAMETER :: iatm_nitr = 2    !! Index for N input from Nitrate N atmospheric deposition
  INTEGER(i_std), PARAMETER :: ibnf      = 3    !! Index for N input from BNF
  INTEGER(i_std), PARAMETER :: imanure   = 4    !! Index for N input from Manure
  INTEGER(i_std), PARAMETER :: ifert_ammo = 5   !! Index for N input from Fertilisation
  INTEGER(i_std), PARAMETER :: ifert_nitr = 6   !! Index for N input from Fertilisation
  INTEGER(i_std), PARAMETER :: ninput    = 6    !! Number of N-input considered  
 

  INTEGER(i_std), PARAMETER :: i_nh4_to_no3 = 1 !! Index for NO3 production
  INTEGER(i_std), PARAMETER :: i_nh4_to_no  = 2 !! Index for NO production
  INTEGER(i_std), PARAMETER :: i_nh4_to_n2o = 3 !! Index for N2O production
  INTEGER(i_std), PARAMETER :: n_nh4_to_x = 3   !! Number of NH4 pathways

  INTEGER(i_std), PARAMETER :: i_no3_to_nox = 1 !! Index for NO3 consumption
  INTEGER(i_std), PARAMETER :: i_nox_to_n2o  = 2!! Index for NO/Nox consumption
  INTEGER(i_std), PARAMETER :: i_n2o_to_n2 = 3  !! Index for N2O consumption
  INTEGER(i_std), PARAMETER :: n_n_to_x = 3     !! Number of N pathways

  INTEGER(i_std), PARAMETER :: nmonth = 12      !! Months in a year; used for input .nc files with monthly arrays
  
  !! Updates for the mass balance closure in stomate_lpj
  INTEGER(i_std), PARAMETER :: ibeg      = 1    !! At the begining of the routine
  INTEGER(i_std), PARAMETER :: ipre      = 2    !! After precribe
  INTEGER(i_std), PARAMETER :: iphe      = 3    !! After phenology
  INTEGER(i_std), PARAMETER :: igro      = 4    !! After growth functional allocation
  INTEGER(i_std), PARAMETER :: iage      = 5    !! After Age class distribution
  INTEGER(i_std), PARAMETER :: iluc      = 6    !! After land cover change
  INTEGER(i_std), PARAMETER :: idis      = 7    !! After land cover change after disturbances 
  INTEGER(i_std), PARAMETER :: icle      = 8    !! After mortality_clean
  INTEGER(i_std), PARAMETER :: irec      = 9    !! After recruitment
  INTEGER(i_std), PARAMETER :: ispc      = 10   !! After Species Change
  INTEGER(i_std), PARAMETER :: nupdate1  = 10   !! Number of step in stomate_lpj where veget_max and atm_to_bm is update 
  INTEGER(i_std), PARAMETER :: ihar      = 11   !! After wood harvest
  INTEGER(i_std), PARAMETER :: imor      = 12   !! After mortality
  INTEGER(i_std), PARAMETER :: itur      = 13   !! After wood harvest
  INTEGER(i_std), PARAMETER :: ilpj      = 14   !! After all subroutines in lpj
  INTEGER(i_std), PARAMETER :: nupdate2  = 14   !! Number of step in stomate_lpj where wood volume is update

  ! These next sets of parameters are now used for both circ_class_kill and
  ! for the harvest_pool.  One source of confusion is what to do with trees that
  ! die from self-thinning or forest dieoffs.  These happen in all forests, regardless
  ! of management strategy.  I decided to put death of this kind into ifm_none, since
  ! it is the only type of mortality found in an unmanaged forest.  If the mortality
  ! does not kill the whole forest (e.g. self thinning), it goes into icut_thin.  If it
  ! does (forest dieoff), it goes into icut_clear.  The biomass is killed in lpj_gap.

  !! Indices used for forest management strategies
  INTEGER(i_std), PARAMETER :: nfm_types = 7             !! The total number of forest management 
                                                         !! strategies we can use
  INTEGER(i_std), PARAMETER :: ifm_none = 1              !! No human intervention in the forests.
  INTEGER(i_std), PARAMETER :: ifm_thin = 2              !! Regular thinning and harvesting of 
                                                         !! wood based on RDI.
  INTEGER(i_std), PARAMETER :: ifm_cop = 3               !! Coppicing for fuelwood.
  INTEGER(i_std), PARAMETER :: ifm_src = 4               !! Short rotation coppices for biomass 
                                                         !! production.
  INTEGER(i_std), PARAMETER :: ifm_uneven = 5            !! Uneven-aged management
  INTEGER(i_std), PARAMETER :: ifm_grass = 6             !! Grazing or cutting
  INTEGER(i_std), PARAMETER :: ifm_crop = 7              !! Crop harvest

  !! Parameter for thinning
  REAL(r_std), PARAMETER    :: r_max = 0.66              !! At r_max*diameter_max the thinning regime switches
                                                         !! from thinning from below to thinning from above.
  
  !! flag to trigger clearcut during spinup 
  INTEGER(i_std), PARAMETER :: flag_spinup_clearcut = 1  !! 

  !! Indices used for harvest pools
  INTEGER(i_std), PARAMETER :: ncut_times = 12           !! The total number of times when trees 
                                                         !! are cut and wood harvested.
  INTEGER(i_std), PARAMETER :: icut_clear = 1            !! A clearcut where all biomass is removed.
  INTEGER(i_std), PARAMETER :: icut_thin = 2             !! Thinning of biomass to reduce the 
                                                         !! number of trees.
  INTEGER(i_std), PARAMETER :: icut_lcc_wood = 3         !! Wood harvest following land cover 
                                                         !! change (LCC)
  INTEGER(i_std), PARAMETER :: icut_lcc_res = 4          !! Site clearing, removal of the stumps 
                                                         !! and branches following LCC
  INTEGER(i_std), PARAMETER :: icut_crop = 5             !! Crop harvest
  INTEGER(i_std), PARAMETER :: icut_grass = 6            !! Grazing or cutting
  INTEGER(i_std), PARAMETER :: icut_cop1 = 7             !! The first coppice cut
  INTEGER(i_std), PARAMETER :: icut_cop2 = 8             !! The second (and subsequent) coppice cut
  INTEGER(i_std), PARAMETER :: icut_cop3 = 9             !! The last coppice cut (only for SRC)
  INTEGER(i_std), PARAMETER :: icut_storm_break = 10     !! Stem breakage due to storm
  INTEGER(i_std), PARAMETER :: icut_storm_uproot = 11    !! Tee uprooting due to storm
  INTEGER(i_std), PARAMETER :: icut_beetle = 12          !! Tree killed due to bark beetle attacks

  !! Indices used to define the product pools
  !! Numbers based on Eggers 2008 - EFI report
  INTEGER(i_std), PARAMETER :: nshort = 1                !! Length in years of the short-lived product pool (GE 1)
  INTEGER(i_std), PARAMETER :: nmedium =17               !! Length in years of the medium-lived product pool (GT 4)
  INTEGER(i_std), PARAMETER :: nlong = 50                !! Length in years of the long-lived product pool (GT 4)

  !! Indices used for land use output variables
  INTEGER(i_std), PARAMETER :: nlanduse = 2              !! Total number of different land uses
  INTEGER(i_std), PARAMETER :: iharvest = 1              !! Biomass removals due to all sorts of harvest
  INTEGER(i_std), PARAMETER :: ilcc = 2                  !! Biomass removals due to all sorts of land cover changes

  !! Indices used to check the mass balance closure
  INTEGER(i_std), PARAMETER :: nmbcomp = 5               !! The total number of components in 
                                                         !! our mass balance check
  INTEGER(i_std), PARAMETER :: iatm2land = 1             !! atmosphere to land fluxes such as GPP 
                                                         !! and co2_2_bm
  INTEGER(i_std), PARAMETER :: iland2atm = 2             !! land to atmosphere fluxes such as Rh, 
                                                         !! Ra and product decomposition
  INTEGER(i_std), PARAMETER :: ilat2out = 3              !! outgoing lateral flux i.e. DOC leaching 
                                                         !! for the litter routine
  INTEGER(i_std), PARAMETER :: ilat2in = 4               !! incoming lateral flux i.e. N deposition 
                                                         !! for the land
  INTEGER(i_std), PARAMETER :: ipoolchange = 5           !! change in pool size i.e. change in biomass

  !! Indices used for warning tracking
  INTEGER(i_std), PARAMETER :: nwarns = 1                !! The total number of warnings we track
  INTEGER(i_std), PARAMETER :: iwphoto = 1               !! A warning about division by zero in photosynthesis

  !! Indices used for wind damage  
  INTEGER(i_std), PARAMETER :: ibreakage = 1             !! The index for stem breakage dur to wind damage
  INTEGER(i_std), PARAMETER :: ioverturning = 2          !! The index for the tree overtuning due to wind damage

  !! Indices for orphan fluxes
  INTEGER(i_std), PARAMETER :: norphans = 8              !! Total number of orphan fluxes (unitless)
  INTEGER(i_std), PARAMETER :: ivegold = 1               !! Index for veget_max before LCC
  INTEGER(i_std), PARAMETER :: ivegnew = 2               !! Index for veget_max before LCC (includes veget_max of orphan fluxes)
  INTEGER(i_std), PARAMETER :: igpp = 3                  !! Index for gpp_daily
  INTEGER(i_std), PARAMETER :: ico2bm = 4                !! Index for co2_to_bm
  INTEGER(i_std), PARAMETER :: irmain = 5                !! Index for maintenance respiration
  INTEGER(i_std), PARAMETER :: irgrow = 6                !! Index for growth respiration 
  INTEGER(i_std), PARAMETER :: inpp = 7                  !! Index for npp_daily
  INTEGER(i_std), PARAMETER :: irhet = 8                 !! Index for total heterotrophic respiration
  !

  !! Indices for phenology
  ! The variable ::plant_status replaces several variables (senescence, 
  ! begin_leaves and allow_phenoinit) that describe the phenological status of the plant
  ! by storing these different aspects of phenology in a single variable, inconsistencies
  ! become impossible or at least easier to check.
  ! When the model starts from scratch the status is set to iprescribe this allows us
  ! to grow leaves from the first year onwards. The plant should then go through the 
  ! different growth phases: ibudsavail, ibudbreak, icanopy, isenescent, idormant and 
  ! finally idead. Following idormant the status should return ibudsavail to initiate 
  ! another cycle in the subsequent growing season. Following idead new vegetation 
  ! should be prescribed.
  INTEGER(i_std), PARAMETER :: inone = 0                 !! No plants thus no status
  INTEGER(i_std), PARAMETER :: iprescribe = 1            !! Prescribe a PFT
  INTEGER(i_std), PARAMETER :: ibudsavail = 2            !! Buds are present
  INTEGER(i_std), PARAMETER :: ibudbreak = 3             !! Day that the buds break and leaf 
                                                         !! on-set begins
  INTEGER(i_std), PARAMETER :: icanopy = 4               !! Canopy is present
  INTEGER(i_std), PARAMETER :: ipresenescence = 5        !! Wood growth cessation
  INTEGER(i_std), PARAMETER :: isenescent = 6            !! The plant is senescent
  INTEGER(i_std), PARAMETER :: idormant = 7              !! The plant is dormant
  INTEGER(i_std), PARAMETER :: idead = 8                 !! The plant was killed

  !
  !! Indices used for analytical spin-up
  INTEGER(i_std), PARAMETER :: nbpools = 12             !! Total number of carbon pools (unitless)
  INTEGER(i_std), PARAMETER :: istructural_above = 1    !! Index for structural litter above (unitless)
  INTEGER(i_std), PARAMETER :: istructural_below = 2    !! Index for structural litter below (unitless)
  INTEGER(i_std), PARAMETER :: imetabolic_above = 3     !! Index for metabolic litter above (unitless)
  INTEGER(i_std), PARAMETER :: imetabolic_below = 4     !! Index for metabolic litter below (unitless)
  INTEGER(i_std), PARAMETER :: iwoody_above = 5         !! Index for woody litter above (unitless)
  INTEGER(i_std), PARAMETER :: iwoody_below = 6         !! Index for woody litter below (unitless)
  INTEGER(i_std), PARAMETER :: isnag_above = 7          !! Index for snag litter above (unitless)
  INTEGER(i_std), PARAMETER :: isnag_below = 8          !! Index for snag litter below (unitless)
  INTEGER(i_std), PARAMETER :: iactive_pool = 9         !! Index for active carbon pool (unitless)
  INTEGER(i_std), PARAMETER :: islow_pool   = 10        !! Index for slow carbon pool (unitless)
  INTEGER(i_std), PARAMETER :: ipassive_pool = 11       !! Index for passive carbon pool (unitless)
  INTEGER(i_std), PARAMETER :: isurface_pool = 12       !! Index for surface carbon pool (unitless)

  !
  !! Indices for the burried carbon and nitrogen following LUC
  INTEGER(i_std), PARAMETER :: nbury=6                  !! Total number of buried pools
  INTEGER(i_std), PARAMETER :: ilitter=1                 !! Index for litter
  INTEGER(i_std), PARAMETER :: ifreshlitter=2           !! Index for fresh litter
  INTEGER(i_std), PARAMETER :: ifreshsom=3              !! Index for fresh som
  INTEGER(i_std), PARAMETER :: ibact=4                  !! Index for bacteria
  INTEGER(i_std), PARAMETER :: isom=5                   !! Index for som
  INTEGER(i_std), PARAMETER :: iminnitrogen=6           !! Index for mineral soil nitrogen

  !
  !! Indicies used for output variables on Landuse tiles defined according to LUMIP project
  !! Note that ORCHIDEE do not represent pasture and urban land. Therefor the variables will have 
  !! val_exp as missing value for these tiles. 
  INTEGER(i_std), PARAMETER :: nlut=4                   !! Total number of landuse tiles according to LUMIP
  INTEGER(i_std), PARAMETER :: id_psl=1                 !! Index for primary and secondary land
  INTEGER(i_std), PARAMETER :: id_crp=2                 !! Index for crop land
  INTEGER(i_std), PARAMETER :: id_pst=3                 !! Index for pasture land
  INTEGER(i_std), PARAMETER :: id_urb=4                 !! Index for urban land

  !! Indices used for canopy structure (Pgap & eff lai) 
  INTEGER(i_std),PARAMETER   :: ndist_types=6           !! the number of distributions we need in the LAI effective routines
  INTEGER(i_std),PARAMETER   :: iheight=1               !! the tree height distribution
  INTEGER(i_std),PARAMETER   :: idiameter=2             !! the trunk diameter distribution
  INTEGER(i_std),PARAMETER   :: icnvol=3                !! the crown volume distribution
  INTEGER(i_std),PARAMETER   :: icnarea=4               !! the crown area distribution
  INTEGER(i_std),PARAMETER   :: icndiaver=5             !! the verticle crown diameter distribution
  INTEGER(i_std),PARAMETER   :: icndiahor=6             !! the horizontal crown diameter distribution

  !! indices used in AR5 output
  INTEGER(i_std),PARAMETER   :: nlctypes=3              !! Number of land cover types
  INTEGER(i_std),PARAMETER   :: iforest=1               !! forest
  INTEGER(i_std),PARAMETER   :: igrass=2                !! grass
  INTEGER(i_std),PARAMETER   :: icrop=3                 !! crops

  !! Indices used in the root profile
  INTEGER(i_std),PARAMETER   :: nroot_prof=2            !! Number of different root profiles
  INTEGER(i_std),PARAMETER   :: istruc=1                !! Structural root profile
  INTEGER(i_std),PARAMETER   :: ifunc=2                 !! Functional root profile for water uptake

  !! Indices used to defile the vertical root rofile
  INTEGER(i_std),PARAMETER   :: ndepths=2               !! Number of components in the definition of a vertical soil profile
  INTEGER(i_std),PARAMETER   :: inode=1                 !! Node
  INTEGER(i_std),PARAMETER   :: iinterface=2            !! Interface


  !
  ! NUMERICAL AND PHYSICS CONSTANTS
  !
  !-
  ! 1. Mathematical and numerical constants
  !-
  REAL(r_std), PARAMETER :: pi = 3.141592653589793238   !! pi souce : http://mathworld.wolfram.com/Pi.html (unitless)
  REAL(r_std), PARAMETER :: euler = 2.71828182845904523 !! e source : http://mathworld.wolfram.com/e.html (unitless)
  REAL(r_std), PARAMETER :: zero = 0._r_std             !! Numerical constant set to 0 (unitless)
  REAL(r_std), PARAMETER :: undemi = 0.5_r_std          !! Numerical constant set to 1/2 (unitless)
  REAL(r_std), PARAMETER :: un = 1._r_std               !! Numerical constant set to 1 (unitless)
  REAL(r_std), PARAMETER :: moins_un = -1._r_std        !! Numerical constant set to -1 (unitless)
  REAL(r_std), PARAMETER :: deux = 2._r_std             !! Numerical constant set to 2 (unitless)
  REAL(r_std), PARAMETER :: trois = 3._r_std            !! Numerical constant set to 3 (unitless)
  REAL(r_std), PARAMETER :: quatre = 4._r_std           !! Numerical constant set to 4 (unitless)
  REAL(r_std), PARAMETER :: cinq = 5._r_std             !![DISPENSABLE] Numerical constant set to 5 (unitless)
  REAL(r_std), PARAMETER :: six = 6._r_std              !![DISPENSABLE] Numerical constant set to 6 (unitless)
  REAL(r_std), PARAMETER :: huit = 8._r_std             !! Numerical constant set to 8 (unitless)
  REAL(r_std), PARAMETER :: mille = 1000._r_std         !! Numerical constant set to 1000 (unitless)

  !-
  ! 2 . Physics
  !-
  REAL(r_std), PARAMETER :: R_Earth = 6378000.              !! radius of the Earth : Earth radius ~= Equatorial radius (m)
  REAL(r_std), PARAMETER :: mincos  = 0.0001                !! Minimum cosine value used for interpolation (unitless) 
  REAL(r_std), PARAMETER :: pb_std = 1013.                  !! standard pressure (hPa)
  REAL(r_std), PARAMETER :: ZeroCelsius = 273.15            !! 0 degre Celsius in degre Kelvin (K)
  REAL(r_std), PARAMETER :: tp_00 = 273.15                  !! 0 degre Celsius in degre Kelvin (K)
  REAL(r_std), PARAMETER :: chalsu0 = 2.8345E06             !! Latent heat of sublimation (J.kg^{-1})
  REAL(r_std), PARAMETER :: chalev0 = 2.5008E06             !! Latent heat of evaporation (J.kg^{-1}) 
  REAL(r_std), PARAMETER :: chalfu0 = chalsu0-chalev0       !! Latent heat of fusion (J.kg^{-1}) 
  REAL(r_std), PARAMETER :: c_stefan = 5.6697E-8            !! Stefan-Boltzman constant (W.m^{-2}.K^{-4})
  REAL(r_std), PARAMETER :: cp_air = 1004.675               !! Specific heat of dry air (J.kg^{-1}.K^{-1}) 
  REAL(r_std), PARAMETER :: cte_molr = 287.05               !! Specific constant of dry air (kg.mol^{-1}) 
  REAL(r_std), PARAMETER :: kappa = cte_molr/cp_air         !! Kappa : ratio between specific constant and specific heat 
                                                            !! of dry air (unitless)
  REAL(r_std), PARAMETER :: msmlr_air = 28.964E-03          !! Molecular weight of dry air (kg.mol^{-1})
  REAL(r_std), PARAMETER :: msmlr_h2o = 18.02E-03           !! Molecular weight of water vapor (kg.mol^{-1}) 
  REAL(r_std), PARAMETER :: cp_h2o = &                      !! Specific heat of water vapor (J.kg^{-1}.K^{-1}) 
       & cp_air*(quatre*msmlr_air)/( 3.5_r_std*msmlr_h2o) 
  REAL(r_std), PARAMETER :: cte_molr_h2o = cte_molr/quatre  !! Specific constant of water vapor (J.kg^{-1}.K^{-1}) 
  REAL(r_std), PARAMETER :: retv = msmlr_air/msmlr_h2o-un   !! Ratio between molecular weight of dry air and water 
                                                            !! vapor minus 1(unitless)  
  REAL(r_std), PARAMETER :: rvtmp2 = cp_h2o/cp_air-un       !! Ratio between specific heat of water vapor and dry air
  
  REAL(r_std), PARAMETER :: rho_h2o= 0.9991_r_std           !! Density of water at 15°C (g cm-3)                                                          !! minus 1 (unitless)
  REAL(r_std), PARAMETER :: cepdu2 = (0.1_r_std)**2         !! Squared wind shear (m^2.s^{-2}) 
  REAL(r_std), PARAMETER :: ct_karman = 0.40_r_std          !! Van Karmann Constant (unitless)
  REAL(r_std), PARAMETER :: cte_grav = 9.80665_r_std        !! Acceleration of the gravity (m.s^{-2})
  REAL(r_std), PARAMETER :: pa_par_hpa = 100._r_std         !! Transform pascal into hectopascal (unitless)
  REAL(r_std), PARAMETER :: RR = 8.314                      !! Ideal gas constant (J.mol^{-1}.K^{-1})
  REAL(r_std), PARAMETER :: Sct = 1370.                     !! Solar constant (W.m^{-2}) 

  REAL(r_std), PARAMETER :: mm_m = 1000._r_std              !! conversion from milimeters to meters

  INTEGER(i_std), SAVE :: testpft = 6
!$OMP THREADPRIVATE(testpft)

  !-
  ! 3. Climatic constants
  !-
  !! Constantes of the Louis scheme 
  REAL(r_std), SAVE :: cb = 5._r_std              !! Constant of the Louis scheme (unitless);
                                                  !! reference to Louis (1979)
!$OMP THREADPRIVATE(cb)
  REAL(r_std), SAVE :: cc = 5._r_std              !! Constant of the Louis scheme (unitless);
                                                  !! reference to Louis (1979)
!$OMP THREADPRIVATE(cc)
  REAL(r_std), SAVE :: cd = 5._r_std              !! Constant of the Louis scheme (unitless);
                                                  !! reference to Louis (1979)
!$OMP THREADPRIVATE(cd)
  REAL(r_std), SAVE :: rayt_cste = 125.           !! Constant in the computation of surface resistance (W.m^{-2})
!$OMP THREADPRIVATE(rayt_cste)
  REAL(r_std), SAVE :: defc_plus = 23.E-3         !! Constant in the computation of surface resistance (K.W^{-1})
!$OMP THREADPRIVATE(defc_plus)
  REAL(r_std), SAVE :: defc_mult = 1.5            !! Constant in the computation of surface resistance (K.W^{-1})
!$OMP THREADPRIVATE(defc_mult)

  !-
  ! 4. Soil thermodynamics constants
  !-
  ! Look at constantes_soil.f90


  !
  ! OPTIONAL PARTS OF THE MODEL
  !
  LOGICAL,PARAMETER :: diag_qsat = .TRUE.         !! One of the most frequent problems is a temperature out of range
                                                  !! we provide here a way to catch that in the calling procedure. 
                                                  !! (from Jan Polcher)(true/false) 
  LOGICAL, SAVE     :: almaoutput =.FALSE.        !! Selects the type of output for the model.(true/false)
                                                  !! Value is read from run.def in intersurf_history
!$OMP THREADPRIVATE(almaoutput)

  !
  ! DIVERSE
  !
  CHARACTER(LEN=100), SAVE :: stomate_forcing_name='NONE'  !! NV080800 Name of STOMATE forcing file (unitless)
                                                           ! Compatibility with Nicolas Viovy driver.
!$OMP THREADPRIVATE(stomate_forcing_name)
  CHARACTER(LEN=100), SAVE :: stomate_Cforcing_name='NONE' !! NV080800 Name of soil forcing file (unitless)
                                                           ! Compatibility with Nicolas Viovy driver.
!$OMP THREADPRIVATE(stomate_Cforcing_name)
  CHARACTER(LEN=100), SAVE :: stomate_Cforcing_discretization_name='NONE' !! Name of soil carbon discretization forcing file (unitless)
!$OMP THREADPRIVATE(stomate_Cforcing_discretization_name)
  INTEGER(i_std), SAVE :: forcing_id                 !! Index of the forcing file (unitless)
!$OMP THREADPRIVATE(forcing_id)
  LOGICAL, SAVE :: allow_forcing_write=.TRUE.        !! Allow writing of stomate_forcing file. 
                                                     !! This variable will be set to false for teststomate.
!$OMP THREADPRIVATE(allow_forcing_write)
  LOGICAL, SAVE :: ok_sort_lon_lat=.FALSE.           !! Activate sort in increasing order of latitudes and longitudes read from the forcing file
!$OMP THREADPRIVATE(ok_sort_lon_lat)



                         !------------------------!
                         !  SECHIBA PARAMETERS    !
                         !------------------------!
 

  !
  ! GLOBAL PARAMETERS   
  !
  REAL(r_std), SAVE :: min_wind = 0.1      !! The minimum wind (m.s^{-1})
!$OMP THREADPRIVATE(min_wind)
  REAL(r_std), PARAMETER :: min_qc = 1.e-4 !! The minimum value for qc (qc=drag*wind) used in coupled(enerbil) and forced mode (enerbil and diffuco) 
  REAL(r_std), SAVE :: snowcri = 1.5       !! Sets the amount above which only sublimation occurs (kg.m^{-2})
!$OMP THREADPRIVATE(snowcri)


  !
  ! FLAGS ACTIVATING SUB-MODELS
  !
  LOGICAL, SAVE :: treat_expansion = .FALSE.           !! Do we treat PFT expansion across a grid point after introduction? (true/false)
!$OMP THREADPRIVATE(treat_expansion)
  LOGICAL, SAVE :: ok_herbivores = .FALSE.             !! flag to activate herbivores (true/false)
!$OMP THREADPRIVATE(ok_herbivores)
  LOGICAL, SAVE :: harvest_agri = .TRUE.               !! flag to harvest aboveground biomass from agricultural PFTs)(true/false)
!$OMP THREADPRIVATE(harvest_agri)
  LOGICAL, SAVE :: lpj_gap_const_mort                  !! constant moratlity (true/false). Default value depend on OK_DGVM.
!$OMP THREADPRIVATE(lpj_gap_const_mort)
  LOGICAL, SAVE :: disable_fire = .TRUE.               !! flag that disable fire (true/false)
!$OMP THREADPRIVATE(disable_fire)
  LOGICAL, SAVE :: spinup_analytic = .FALSE.           !! Flag to activate analytical resolution for spinup (true/false)
!$OMP THREADPRIVATE(spinup_analytic)
  LOGICAL, SAVE :: ok_soil_carbon_discretization       !! Flag to activate soil carbon discretization (vertical carbon and soil carbon thermal insulation)
!$OMP THREADPRIVATE(ok_soil_carbon_discretization)
  LOGICAL, SAVE :: new_carbinput_intdepzlit = .FALSE.  !! Flag to activate a different method to redistribute the litter input into the different soil layers.
!$OMP THREADPRIVATE(new_carbinput_intdepzlit)
LOGICAL, SAVE :: ok_forcesoil_write                    !! Flag to activate soil carbon discretization output file to be used for forcesoil program
!$OMP THREADPRIVATE(ok_forcesoil_write)
  LOGICAL, SAVE :: ok_vessel_mortality                 !! Flag to activate death and recovery of vegetation following hydraulic failure.
!$OMP THREADPRIVATE(ok_vessel_mortality)
  LOGICAL, SAVE :: ok_ice_sheet                !! Flag to activate ice scheme on ice-sheet area
!$OMP THREADPRIVATE(ok_ice_sheet)

  !
  ! CONFIGURATION VEGETATION
  !
  LOGICAL, SAVE :: agriculture = .TRUE.    !! allow agricultural PFTs (true/false)
!$OMP THREADPRIVATE(agriculture)
  LOGICAL, SAVE :: impveg = .FALSE.        !! Impose vegetation ? (true/false)
!$OMP THREADPRIVATE(impveg)
  LOGICAL, SAVE :: impsoilt = .FALSE.      !! Impose soil ? (true/false)
!$OMP THREADPRIVATE(impsoilt)
  LOGICAL, SAVE :: impslope = .FALSE.      !! Impose reinf_slope ? (true/false)
!$OMP THREADPRIVATE(impslope)
  LOGICAL, SAVE :: impose_ninput_dep = .FALSE. !! Impose N input values from the atmosphere ? (true/false)
!$OMP THREADPRIVATE(impose_ninput_dep)
  LOGICAL, SAVE :: impose_ninput_fert = .FALSE. !! Impose N input values from fertilizer ? (true/false)        
!$OMP THREADPRIVATE(impose_ninput_fert)    
 LOGICAL, SAVE :: impose_ninput_manure = .FALSE. !! Impose N input values from manure? (true/false)        
!$OMP THREADPRIVATE(impose_ninput_manure)                     
  LOGICAL, SAVE :: impose_ninput_bnf = .FALSE. !! Impose N input values from biological nitrogen fixation (BNF) ? (true/false)        
!$OMP THREADPRIVATE(impose_ninput_bnf) 

  LOGICAL, SAVE :: do_now_stomate_lcchange = .FALSE.  !! Time to call lcchange in stomate_lpj
!$OMP THREADPRIVATE(do_now_stomate_lcchange)
  LOGICAL, SAVE :: do_now_stomate_woodharvest = .FALSE.  !! Time to call woodharvest in stomate_lpj
!$OMP THREADPRIVATE(do_now_stomate_woodharvest)
  LOGICAL, SAVE :: do_now_stomate_product_use = .FALSE.  !! Time to call wood product use and decomposition
!$OMP THREADPRIVATE(do_now_stomate_product_use)
  LOGICAL, SAVE :: done_stomate_lcchange = .FALSE.    !! If true, call lcchange in stomate_lpj has just been done. 
!$OMP THREADPRIVATE(done_stomate_lcchange)
  LOGICAL, SAVE :: do_now_recruit = .FALSE.  !! Time to call prescribe to calculate recruits in stomate_lpj
!$OMP THREADPRIVATE(do_now_recruit)
  INTEGER, SAVE :: read_lai                  !! Option for LAI if STOMATE is not activated (0: no LAI reading, 1: site-based LAI reading, 2: read LAI map)
!$OMP THREADPRIVATE(read_lai)
  LOGICAL, SAVE :: vegetmap_reset = .FALSE.  !! Reset the vegetation map and reset carbon related variables
!$OMP THREADPRIVATE(vegetmap_reset)
  INTEGER(i_std) , SAVE :: veget_update      !! Update frequency in years for landuse (nb of years)
!$OMP THREADPRIVATE(veget_update)
  LOGICAL, SAVE :: ninput_reinit = .TRUE.    !! To change N INPUT file in a run. (true/false)
!$OMP THREADPRIVATE(ninput_reinit)


  !
  ! PARAMETERS USED BY BOTH HYDROLOGY MODELS
  !
  REAL(r_std), SAVE :: max_snow_age = 50._r_std !! Maximum period of snow aging (days)
!$OMP THREADPRIVATE(max_snow_age)
  REAL(r_std), SAVE :: snow_trans = 0.2_r_std   !! Transformation time constant for snow (m), reduced from the value 0.3 (04/07/2016)
!$OMP THREADPRIVATE(snow_trans)
  REAL(r_std), SAVE :: snow_trans_nobio = 1.0_r_std   !! Transformation time constant for snow on nobio (m), (18/02/2019)
!$OMP THREADPRIVATE(snow_trans_nobio)
  REAL(r_std), SAVE :: snowa_aged_nobio_vis = 0.58_r_std   !! Minimum snow albedo value for nobio type after aging (dirty old snow), visible albedo (09/12/2022)
!$OMP THREADPRIVATE(snowa_aged_nobio_vis)
  REAL(r_std), SAVE :: snowa_aged_nobio_nir = 0.35_r_std   !! Minimum snow albedo value for nobio type after aging (dirty old snow), near infrared albedo (09/12/2022)
!$OMP THREADPRIVATE(snowa_aged_nobio_nir)
  REAL(r_std), SAVE :: snowa_dec_nobio_vis = 0.28_r_std   !! Decay rate of snow albedo value for nobio type, visible albedo (09/12/2022)
!$OMP THREADPRIVATE(snowa_dec_nobio_vis)
  REAL(r_std), SAVE :: snowa_dec_nobio_nir = 0.45_r_std   !! Decay rate of snow albedo value for nobio type, near infrared albedo (09/12/2022)
!$OMP THREADPRIVATE(snowa_dec_nobio_nir)
  REAL(r_std), SAVE :: sneige                   !! Lower limit of snow amount (kg.m^{-2})
!$OMP THREADPRIVATE(sneige)
  REAL(r_std), SAVE :: maxmass_snow = 3000.     !! The maximum mass of snow (kg.m^{-2})
!$OMP THREADPRIVATE(maxmass_snow)

  !! Heat capacity
  REAL(r_std), PARAMETER :: rho_water = 1000.           !! Density of water (kg/m3)
  REAL(r_std), PARAMETER :: rho_ice = 920.              !! Density of ice (kg/m3)
  REAL(r_std), PARAMETER :: rho_soil = 2700.            !! Density of soil particles (kg/m3), value from Peters-Lidard et al. 1998

  !! Thermal conductivities
  REAL(r_std), PARAMETER :: cond_water = 0.6            !! Thermal conductivity of liquid water (W/m/K)
  REAL(r_std), PARAMETER :: cond_ice = 2.2              !! Thermal conductivity of ice (W/m/K)
  REAL(r_std), PARAMETER :: cond_solid = 2.32           !! Thermal conductivity of mineral soil particles (W/m/K)

  !! Time constant of long-term soil humidity (s) 
  REAL(r_std), PARAMETER :: lhf = 0.3336*1.E6           !! Latent heat of fusion (J/kg)

  INTEGER(i_std), SAVE      :: nsnow                    !! Number of levels in the snow for explicit snow scheme
!$OMP THREADPRIVATE(nsnow)
  INTEGER(i_std), SAVE      :: nice                     !! Number of levels in ice (for ice-sheet area)
!$OMP THREADPRIVATE(nice)
  REAL(r_std), SAVE         :: frac_min_ice             !! Limit of fraction of ice in the grid cell where the ice mask is considered as 1
!$OMP THREADPRIVATE(frac_min_ice)
  LOGICAL, SAVE             :: read_ice_sheet_mask      !! If true, read the ice sheet mask from file
!$OMP THREADPRIVATE(read_ice_sheet_mask)

  REAL(r_std), PARAMETER    :: XMD    = 28.9644E-3 
  REAL(r_std), PARAMETER    :: XBOLTZ      = 1.380658E-23 
  REAL(r_std), PARAMETER    :: XAVOGADRO   = 6.0221367E+23 
  REAL(r_std), PARAMETER    :: XRD    = XAVOGADRO * XBOLTZ / XMD 
  REAL(r_std), PARAMETER    :: XCPD   = 7.* XRD /2. 
  REAL(r_std), PARAMETER    :: phigeoth = 0.057 ! 0. DKtest 
  REAL(r_std), PARAMETER    :: thick_min_snow = .01 

  !! The maximum snow density and water holding characterisicts 
  REAL(r_std), SAVE         :: xrhosmax = 750.  ! (kg m-3) 
!$OMP THREADPRIVATE(xrhosmax)
  REAL(r_std), SAVE         :: xwsnowholdmax1   = 0.03  ! (-) 
!$OMP THREADPRIVATE(xwsnowholdmax1)
  REAL(r_std), SAVE         :: xwsnowholdmax2   = 0.10  ! (-) 
!$OMP THREADPRIVATE(xwsnowholdmax2)
  REAL(r_std), SAVE         :: xsnowrhohold     = 200.0 ! (kg/m3) 
!$OMP THREADPRIVATE(xsnowrhohold)
  REAL(r_std), SAVE         :: xrhosmin = 50. 
!$OMP THREADPRIVATE(xrhosmin)
  REAL(r_std), PARAMETER    :: xci = 2.106e+3 
  REAL(r_std), PARAMETER    :: xrv = 6.0221367e+23 * 1.380658e-23 /18.0153e-3 

  !! ISBA-ES Critical snow depth at which snow grid thicknesses constant 
  REAL(r_std), PARAMETER    :: xsnowcritd = 0.03  ! (m) 

  !! The threshold of snow depth used for preventing numerical problem in thermal calculations
  REAL(r_std), PARAMETER    :: snowcritd_thermal = 0.01  ! (m)  
  
  !! ISBA-ES CROCUS (Pahaut 1976): snowfall density coefficients: 
  REAL(r_std), PARAMETER       :: snowfall_a_sn = 109.0  !! (kg/m3) 
  REAL(r_std), PARAMETER       :: snowfall_b_sn =   6.0  !! (kg/m3/K) 
  REAL(r_std), PARAMETER       :: snowfall_c_sn =  26.0  !! [kg/(m7/2 s1/2)] 

  REAL(r_std), PARAMETER       :: dgrain_new_max=  2.0e-4!! (m) : Maximum grain size of new snowfall 
  
  !! Used in explicitsnow to prevent numerical problems as snow becomes vanishingly thin. 
  REAL(r_std), PARAMETER                :: psnowdzmin = .0001   ! m 
  REAL(r_std), PARAMETER                :: xsnowdmin = .000001  ! m 

  REAL(r_std), PARAMETER                :: ph2o = 1000.         !! Water density [kg/m3] 
  
  ! ISBA-ES Thermal conductivity coefficients from Anderson (1976): 
  ! see Boone, Meteo-France/CNRM Note de Centre No. 70 (2002) 
  REAL(r_std), SAVE                     :: ZSNOWTHRMCOND1 = 0.02    ! [W/m/K] 
!$OMP THREADPRIVATE(ZSNOWTHRMCOND1)
  REAL(r_std), SAVE                     :: ZSNOWTHRMCOND2 = 2.5E-6  ! [W m5/(kg2 K)] 
!$OMP THREADPRIVATE(ZSNOWTHRMCOND2)
  
  ! Thermal ice conductivity coefficients Yin CHao Yen 1981 (Temp > 195K): 
  REAL(r_std), SAVE                     :: ZICETHRMCOND1 = 6.727   ! [W/m/K] 
  REAL(r_std), SAVE                     :: ZICETHRMCOND2 = -0.0041 !
  ! Thermal ice heat capacity coefficients : 
  REAL(r_std), SAVE                     :: ZICETHRMHEAT1 = 2115.3  !
  REAL(r_std), SAVE                     :: ZICETHRMHEAT2 = 7.79293 !
  
  ! ISBA-ES Thermal conductivity: Implicit vapor diffn effects 
  ! (sig only for new snow OR high altitudes) 
  ! from Sun et al. (1999): based on data from Jordan (1991) 
  ! see Boone, Meteo-France/CNRM Note de Centre No. 70 (2002) 
  ! 
  REAL(r_std), SAVE                       :: ZSNOWTHRMCOND_AVAP  = -0.06023 ! (W/m/K) 
!$OMP THREADPRIVATE(ZSNOWTHRMCOND_AVAP)
  REAL(r_std), SAVE                       :: ZSNOWTHRMCOND_BVAP  = -2.5425  ! (W/m) 
!$OMP THREADPRIVATE(ZSNOWTHRMCOND_BVAP)
  REAL(r_std), SAVE                       :: ZSNOWTHRMCOND_CVAP  = -289.99  ! (K) 
!$OMP THREADPRIVATE(ZSNOWTHRMCOND_CVAP)
  REAL(r_std),SAVE :: xansmax = 0.85      !! Maxmimum snow albedo
!$OMP THREADPRIVATE(xansmax)
  REAL(r_std),SAVE :: xansmin = 0.50      !! Miniumum snow albedo
!$OMP THREADPRIVATE(xansmin)
  REAL(r_std),SAVE :: xans_todry = 0.008  !! Albedo decay rate for dry snow
!$OMP THREADPRIVATE(xans_todry)
  REAL(r_std),SAVE :: xans_t = 0.240      !! Albedo decay rate for wet snow
!$OMP THREADPRIVATE(xans_t)

  ! ISBA-ES Thermal conductivity coefficients from Anderson (1976):
  ! see Boone, Meteo-France/CNRM Note de Centre No. 70 (2002)
  REAL(r_std), PARAMETER                  :: XP00 = 1.E5

  ! ISBA-ES Thermal conductivity: Implicit vapor diffn effects
  ! (sig only for new snow OR high altitudes)
  ! from Sun et al. (1999): based on data from Jordan (1991)
  ! see Boone, Meteo-France/CNRM Note de Centre No. 70 (2002)
  !
  REAL(r_std), SAVE          :: ZSNOWCMPCT_RHOD  = 150.0        !! (kg/m3)
!$OMP THREADPRIVATE(ZSNOWCMPCT_RHOD)
  REAL(r_std), SAVE          :: ZSNOWCMPCT_ACM   = 2.8e-6       !! (1/s)
!$OMP THREADPRIVATE(ZSNOWCMPCT_ACM)
  REAL(r_std), SAVE          :: ZSNOWCMPCT_BCM   = 0.04         !! (1/K) 
!$OMP THREADPRIVATE(ZSNOWCMPCT_BCM)
  REAL(r_std), SAVE          :: ZSNOWCMPCT_CCM   = 460.         !! (m3/kg)
!$OMP THREADPRIVATE(ZSNOWCMPCT_CCM)
  REAL(r_std), SAVE          :: ZSNOWCMPCT_V0    = 3.7e7        !! (Pa/s) 
!$OMP THREADPRIVATE(ZSNOWCMPCT_V0)
  REAL(r_std), SAVE          :: ZSNOWCMPCT_VT    = 0.081        !! (1/K)
!$OMP THREADPRIVATE(ZSNOWCMPCT_VT)
  REAL(r_std), SAVE          :: ZSNOWCMPCT_VR    = 0.018        !! (m3/kg)
!$OMP THREADPRIVATE(ZSNOWCMPCT_VR)


  ! Parameters for viscosity if T > 264 K
  REAL(r_std), SAVE          :: ZSNOWCMPCT_W0    = 2.0e7        !! (Pa/s) 
  REAL(r_std), SAVE          :: ZSNOWCMPCT_WT    = 0.35         !! (1/K)
  REAL(r_std), SAVE          :: ZSNOWCMPCT_WR    = 0.01         !! (m3/kg)

  ! Constants for use of the compaction scheme from Decharme et al.(2016)
  REAL(r_std), SAVE          :: xrhosref         = 250.          !!  (kg/m3)
  REAL(r_std), SAVE          :: xrhosmob         = 295.          !!  (kg/m3)
  REAL(r_std), SAVE          :: xrhoswind        = 350.          !!  (kg/m3)

  REAL(r_std), SAVE          :: ZSNOWCMPCT_X0    = 7.62237e6     ! 7.82237e6     !! (Pa/s)

  REAL(r_std), SAVE          :: ZSNOWCMPCT_XT    = 0.1           !! (1/K)
  REAL(r_std), SAVE          :: ZSNOWCMPCT_XR    = 0.023         !! (m3/kg)
  REAL(r_std), SAVE          :: DTref            = 5.            !! (K)

  REAL(r_std), SAVE          :: ACMPCT           = 2.868
  REAL(r_std), SAVE          :: BCMPCT           = 0.085         !! (s/m)
  REAL(r_std), SAVE          :: KCMPCT            = 1.25         !! (unitless)

  REAL(r_std), SAVE          :: a_tau            = 10            !! 
  REAL(r_std), SAVE          :: b_tau            = 3.25          !!

  REAL(r_std), SAVE          :: Amob             = 1.25          !! (unitless) 

  !
  ! PARAMETERS USED FOR CANOPY LAYERS (Albedo, photosynthesis, energy budget)
  !

  INTEGER(i_std), PARAMETER   :: nlevels = 1                  !! Originally the number of levels in the canopy used in  
                                                              !! calculation of the energy budget. After the
                                                              !! mleb calculations have been implemented in enerbil, 
                                                              !! the jnlvls levels are determining the levels used 
                                                              !! within the multi-layer energy budget calculations. However, 
                                                              !! nlevels are still used in calculate_z_level_photo. Cannot 
                                                              !! be deleted before any decision regarding the vertical 
                                                              !! layering has been made.
  INTEGER(i_std), SAVE        :: nlai = 10                    !! Number of levels in the canopy used in the photosynthesis
                                                              !! routine per level dictacted by nlevels. For example, if
                                                              !! if nlevels = 2 and nlai = 3, the photosynthesis
                                                              !! will be calculated for nlevels_tot=2*3=6 total levels. 
!$OMP THREADPRIVATE(nlai)
  INTEGER(i_std), SAVE        :: nlevels_tot                  !! Total number of levels, nlevels*nlai. Currently, 
                                                              !! nlevels=1 is used, thus nlevels_tot=nlai.
                                                              !! Note that when using the multi-layer budget nlai 
                                                              !! needs to be nlai=jnlvls_canopy+1

!$OMP THREADPRIVATE(nlevels_tot)
  INTEGER(i_std), SAVE        :: jnlvls=29                    !! Number of levels in the multilayer energy budget scheme
!$OMP THREADPRIVATE(jnlvls)
  INTEGER(i_std), SAVE        :: jnlvls_under=10              !! Number of levels in the understorey of the multilayer energy budget scheme
!$OMP THREADPRIVATE(jnlvls_under)
  INTEGER(i_std), SAVE        :: jnlvls_canopy=9              !! Number of levels in the canopy of the multilayer energy budget scheme
!$OMP THREADPRIVATE(jnlvls_canopy)
  INTEGER(i_std), SAVE        :: jnlvls_over=10               !! Number of levels in the overstorey of the multilayer energy budget scheme
!$OMP THREADPRIVATE(jnlvls_over)


  INTEGER(i_std), SAVE        :: nlev_top                     !! Maximum number of canopy levels that are used to construct the "top"
                                                              !! layer of the canopy. The top layer is used in the calculation
                                                              !! transpiration.
!$OMP THREADPRIVATE(nlev_top) 
  REAL(r_std), PARAMETER, &
         DIMENSION (nlevels) :: z_level = (/ 0.0 /)           !! The height of the bottom of each canopy layer                                              
                                                              !! @tex $(m)$ @endtex
!$OMP THREADPRIVATE(z_level)



  !
  ! Parameters for determining the effective LAI for use in Pinty's albedo scheme
  !
  REAL(r_std), SAVE          ::  laieff_solar_angle           !! the zenith angle of the sun which determines our effective LAI
                                                              !! Pinty et al recommend a value of 60 degrees for this regadless of the true
                                                              !! solar zenith angle
!$OMP THREADPRIVATE(laieff_solar_angle)
  REAL(r_std), SAVE          ::  laieff_zero_cutoff           !! an arbitrary cutoff to prevent too low of values from being passed to
                                                              !! routines in the calculation of the effective LAI
!$OMP THREADPRIVATE(laieff_zero_cutoff)
  REAL(r_std), SAVE          ::  direct_light_weight          !! an arbitrary weight (between 0 and 1) used when calculating the albedo and
                                                              !! absorbed light in the canopy.  The amounts are first calculated using one unit of
                                                              !! incoming light from either source, and then combined afterwards with this fraction.
                                                              !! This fraction represents the light coming from
                                                              !! the sun directory, while (1.0-this fraction) gives the amount coming
                                                              !! from light first reflected off clouds and albedo.
!$OMP THREADPRIVATE(direct_light_weight)

  !
  ! PARAMETERS FOR HYDRAULIC ARCHITECTURE
  !

  REAL(r_std), SAVE, DIMENSION(2)         :: a_viscosity = (/0.556,0.022/) !! Empirical parameters to adjust the resistance of fine
                                                                           !! root and sapwood to the temperature dependency of the
                                                                           !! viscosity of water Cochard et al 2000
!$OMP THREADPRIVATE(a_viscosity)

  !
  ! BVOC : Biogenic activity  for each age class
  !
  REAL(r_std), SAVE, DIMENSION(nleafages) :: iso_activity = (/0.5, 1.5, 1.5, 0.5/)     !! Biogenic activity for each 
                                                                                       !! age class : isoprene (unitless)
!$OMP THREADPRIVATE(iso_activity)
  REAL(r_std), SAVE, DIMENSION(nleafages) :: methanol_activity = (/1., 1., 0.5, 0.5/)  !! Biogenic activity for each
                                                                                       !! age class : methanol (unnitless)
!$OMP THREADPRIVATE(methanol_activity)


 !
 ! Parameters for irrigation scheme
 !
  REAL(r_std), SAVE :: irrig_dosmax = 3.              !! The maximum irrigation water injected per hour (kg.m^{-2}/hour)
!$OMP THREADPRIVATE(irrig_dosmax)
  REAL(r_std), SAVE :: cum_dh_thr = 0.64              !! Cumulated nroot threshoold to define root zone, and calculate water deficit for irrigation (-)
!$OMP THREADPRIVATE(cum_dh_thr)
  LOGICAL, SAVE :: irrigated_soiltile = .FALSE.       !! Do we introduce a new soil tile for irrigated croplands? (true/false)
!$OMP THREADPRIVATE(irrigated_soiltile)
  LOGICAL, SAVE :: old_irrig_scheme = .FALSE.         !! Do we run with the old irrigation scheme? (true/false)  , add to compatiblity
!$OMP THREADPRIVATE(old_irrig_scheme)
  INTEGER, SAVE :: irrig_st = 3                       !! Which is the soil tile with irrigation flux
!$OMP THREADPRIVATE(irrig_st)
  REAL(r_std), SAVE, DIMENSION(3) :: avail_reserve = (/0.9,0.9,0.9/)     !! Available water from routing reservoirs, to withdraw for irrigation
                                                      !! IMPORTANT: As the routing model uses 3 reservoirs, dimension is set to 3
                                                      !! IMPORTANT: Order of available water must be in this order: streamflow, fast, and slow reservoir
!$OMP THREADPRIVATE(avail_reserve)
  REAL(r_std), SAVE :: beta_irrig = 0.9               !! Threshold multiplier of Target SM to calculate root deficit(unitless)
!$OMP THREADPRIVATE(beta_irrig)
  REAL(r_std), SAVE :: lai_irrig_min = 0.1            !! Minimum LAI to trigger irrigation (kg.m^{-2}/hour)
!$OMP THREADPRIVATE(lai_irrig_min)
  LOGICAL, SAVE :: irrig_map_dynamic_flag = .FALSE.   !! Do we use a dynamic irrig map?
!$OMP THREADPRIVATE(irrig_map_dynamic_flag)
  LOGICAL, SAVE :: select_source_irrig = .FALSE.      !! Do we use the new priorization scheme, based on maps of equipped area with surface water?
!$OMP THREADPRIVATE(select_source_irrig)
  LOGICAL, SAVE :: Reinfiltr_IrrigField = .FALSE.     !! Do we reinfiltrate all runoff from crop soil tile?O
!$OMP THREADPRIVATE(Reinfiltr_IrrigField)
  REAL, SAVE :: reinf_slope_cropParam = 0.8           !! Externalized for irrigated cropland, when Reinfiltr_IrrigField=.TRUE.
                                                      !! Max value of reinf_slope in irrig_st  
!$OMP THREADPRIVATE(reinf_slope_cropParam)
  REAL, SAVE :: a_stream_adduction = 0.05             !! Externalized for available volume to adduction
!$OMP THREADPRIVATE(a_stream_adduction)



  !
  ! condveg.f90
  !

  ! 1. Scalar

  ! 1.1 Flags used inside the module

  LOGICAL, SAVE :: alb_bare_model = .FALSE. !! Switch for choosing values of bare soil 
                                            !! albedo (see header of subroutine)
                                            !! (true/false)
!$OMP THREADPRIVATE(alb_bare_model)
  LOGICAL, SAVE :: alb_bg_modis = .FALSE.   !! Switch for choosing values of bare soil 
                                            !! albedo read from file
                                            !! (true/false)
!$OMP THREADPRIVATE(alb_bg_modis)
  LOGICAL, SAVE :: impaze = .FALSE.         !! Switch for choosing surface parameters
                                            !! (see header of subroutine).  
                                            !! (true/false)
!$OMP THREADPRIVATE(impaze)
  LOGICAL, SAVE :: rough_dyn = .TRUE.       !! Chooses between two methods to calculate the 
                                            !! the roughness height : static or dynamic (varying with LAI)
                                            !! (true/false)
!$OMP THREADPRIVATE(rough_dyn)
  LOGICAL, SAVE :: use_ratio_z0m_z0h = .TRUE. !! To impose a constant ratio as done in ROUGH_DYN=F
!$OMP THREADPRIVATE(use_ratio_z0m_z0h)

  LOGICAL, SAVE :: sla_dyn = .TRUE.        !! Chooses between two methods to calculate the 
                                            !! specific leaf area: static or dynamic (varying with LAI or biomass)
                                            !! (true/false)
!$OMP THREADPRIVATE(sla_dyn)

  LOGICAL, SAVE :: new_watstress = .FALSE.
!$OMP THREADPRIVATE(new_watstress)

  REAL(r_std), SAVE :: alpha_watstress = 1.
!$OMP THREADPRIVATE(alpha_watstress)
  LOGICAL, SAVE :: use_height_dom = .FALSE.       !! A logical flag determining if we use
                                                  !! the dominant height of the vegetation (.TRUE.) 
                                                  !! or the quadratic mean height (.FALSE.) when
                                                  !! the roughness height is calculated
!$OMP THREADPRIVATE(use_height_dom)

  ! 1.2 Others 

  REAL(r_std), SAVE :: height_displacement = 0.66        !! Factor to calculate the zero-plane displacement
                                                         !! height from vegetation height (m)
!$OMP THREADPRIVATE(height_displacement)
  REAL(r_std), SAVE :: crown_packing = 1.0               !! Parameter to describe how much of the canopy space could 
                                                         !! be filled with crowns (think: close packing of equal spheres)
!$OMP THREADPRIVATE(crown_packing)
  REAL(r_std), SAVE :: alpha_height_precip = 3.22        !! alpha parameter in height-precipitation relationship
                                                         !! used when test_dynamic_max_height=y
  !$OMP THREADPRIVATE(alpha_height_precip)
  REAL(r_std), SAVE :: beta_height_precip = 0.33877      !! beta parameter in height-precipitation relationship
                                                         !! used when test_dynamic_max_height=y
!$OMP THREADPRIVATE(beta_height_precip)  
  REAL(r_std), SAVE :: z0_bare = 0.01                    !! bare soil roughness length (m)
!$OMP THREADPRIVATE(z0_bare)
  REAL(r_std), SAVE :: z0_ice = 0.001                    !! ice roughness length (m) --> Standard value
!$OMP THREADPRIVATE(z0_ice)
  REAL(r_std), SAVE :: tcst_snowa = 10.0                 !! Time constant of the albedo decay of snow (days)
!$OMP THREADPRIVATE(tcst_snowa)
  REAL(r_std), SAVE :: tcst_snowa_nobio = 2.0            !! Time constant of the albedo decay of snow on nobio (days) (18/02/2018)
!$OMP THREADPRIVATE(tcst_snowa_nobio)
REAL(r_std), SAVE :: omg1 = 7.0                          !! First tuning constant for snow ageing on nobio areas (14/10/2020)
!$OMP THREADPRIVATE(omg1)
REAL(r_std), SAVE :: omg2 = 4.0                          !! Second tuning constant for snow ageing on nobio areas (14/10/2020)
!$OMP THREADPRIVATE(omg2)
  REAL(r_std), SAVE :: snowcri_alb = 10.                 !! Critical value for computation of snow albedo (cm) 
!$OMP THREADPRIVATE(snowcri_alb)
  REAL(r_std), SAVE :: fixed_snow_albedo = undef_sechiba !! To choose a fixed snow albedo value (unitless)
!$OMP THREADPRIVATE(fixed_snow_albedo)
  REAL(r_std), SAVE :: z0_scal = 0.15                    !! Surface roughness height imposed (m)
!$OMP THREADPRIVATE(z0_scal)
  REAL(r_std), SAVE :: roughheight_scal = zero           !! Effective roughness Height depending on zero-plane 
                                                         !! displacement height (m) (imposed)
!$OMP THREADPRIVATE(roughheight_scal)
  REAL(r_std), SAVE :: emis_scal = 1.0                   !! Surface emissivity imposed (unitless)
!$OMP THREADPRIVATE(emis_scal)

  REAL(r_std), SAVE :: c1 = 0.32                         !! Constant used in the formulation of the ratio of 
!$OMP THREADPRIVATE(c1)                                  !! friction velocity to the wind speed at the canopy top
                                                         !! see Ershadi et al. (2015) for more info
  REAL(r_std), SAVE :: c2 = 0.264                        !! Constant used in the formulation of the ratio of 
!$OMP THREADPRIVATE(c2)                                  !! friction velocity to the wind speed at the canopy top
                                                         !! see Ershadi et al. (2015) for more info
  REAL(r_std), SAVE :: c3 = 15.1                         !! Constant used in the formulation of the ratio of 
!$OMP THREADPRIVATE(c3)                                  !! friction velocity to the wind speed at the canopy top
                                                         !! see Ershadi et al. (2015) for more info
  REAL(r_std), SAVE :: Cdrag_foliage = 0.2               !! Drag coefficient of the foliage
!$OMP THREADPRIVATE(Cdrag_foliage)                       !! See Ershadi et al. (2015) and Su et. al (2001) for more info
  REAL(r_std), SAVE :: Ct = 0.01                         !! Heat transfer coefficient of the leaf
!$OMP THREADPRIVATE(Ct)                                  !! See Ershadi et al. (2015) and Su et. al (2001) for more info
  REAL(r_std), SAVE :: Prandtl = 0.71                    !! Prandtl number used in the calculation of Ct_star
!$OMP THREADPRIVATE(Prandtl)                             !! See Su et. al (2001) for more info



  ! 2. Arrays

  
! 2. Arrays
 
  REAL(r_std), SAVE, DIMENSION(n_spectralbands) :: alb_deadleaf = (/ .12, .35/)    !! albedo of dead leaves, VIS+NIR (unitless)
!$OMP THREADPRIVATE(alb_deadleaf)
  REAL(r_std), SAVE, DIMENSION(n_spectralbands) :: alb_ice = (/ .60, .20/)         !! albedo of ice, VIS+NIR (unitless)
!$OMP THREADPRIVATE(alb_ice)
  REAL(r_std), SAVE, DIMENSION(n_spectralbands) :: albedo_scal = (/ 0.25, 0.25 /)                !! Albedo values for visible and near-infrared 
                                                                                   !! used imposed (unitless) 
!$OMP THREADPRIVATE(albedo_scal)
  REAL(r_std), DIMENSION(n_spectralbands), SAVE :: alb_snow_0 = (/.95, .65/)       !! VIS and NIR albedo of fresh snow according to CLM3
!$OMP THREADPRIVATE(alb_snow_0)
  REAL(r_std), DIMENSION(n_spectralbands), SAVE :: c_albedo = (/0.2, 0.5/)         !! VIS and NIR albedo constant according to CLM3
!$OMP THREADPRIVATE(c_albedo)
  REAL(r_std) , SAVE, DIMENSION(classnb) :: vis_dry = (/0.24,&
       &0.22, 0.20, 0.18, 0.16, 0.14, 0.12, 0.10, 0.27/)  !! Soil albedo values to soil colour classification:
                                                          !! dry soil albedo values in visible range
!$OMP THREADPRIVATE(vis_dry)
  REAL(r_std), SAVE, DIMENSION(classnb) :: nir_dry = (/0.48,&
       &0.44, 0.40, 0.36, 0.32, 0.28, 0.24, 0.20, 0.55/)  !! Soil albedo values to soil colour classification:
                                                          !! dry soil albedo values in near-infrared range 
!$OMP THREADPRIVATE(nir_dry)
  REAL(r_std), SAVE, DIMENSION(classnb) :: vis_wet = (/0.12,&
       &0.11, 0.10, 0.09, 0.08, 0.07, 0.06, 0.05, 0.15/)  !! Soil albedo values to soil colour classification:
                                                          !! wet soil albedo values in visible range 
!$OMP THREADPRIVATE(vis_wet)
  REAL(r_std), SAVE, DIMENSION(classnb) :: nir_wet = (/0.24,&
       &0.22, 0.20, 0.18, 0.16, 0.14, 0.12, 0.10, 0.31/)  !! Soil albedo values to soil colour classification:
                                                          !! wet soil albedo values in near-infrared range
!$OMP THREADPRIVATE(nir_wet)
  REAL(r_std), SAVE, DIMENSION(classnb) :: albsoil_vis = (/ &
       &0.18, 0.16, 0.16, 0.15, 0.12, 0.105, 0.09, 0.075, 0.25/)   !! Soil albedo values to soil colour classification:
                                                                   !! Averaged of wet and dry soil albedo values
                                                                   !! in visible and near-infrared range
!$OMP THREADPRIVATE(albsoil_vis) 
  REAL(r_std), SAVE, DIMENSION(classnb) :: albsoil_nir = (/ &
       &0.36, 0.34, 0.34, 0.33, 0.30, 0.25, 0.20, 0.15, 0.45/)  !! Soil albedo values to soil colour classification:
                                                                !! Averaged of wet and dry soil albedo values
                                                                !! in visible and near-infrared range
!$OMP THREADPRIVATE(albsoil_nir)
  REAL(r_std) :: alb_threshold = 0.0000000001_r_std       !! A threshold for the iteration of the
                                                          !! multilevel albedo.  Could be externalised.
                                                          !! Fairly arbitrary, although if a level has
                                                          !! no LAI the absorption often ends up being
                                                          !! equal to this value, so it should not
                                                          !! be high.
!$OMP THREADPRIVATE(alb_threshold)
  REAL(r_std), SAVE, DIMENSION(8) :: ZSICOEF1  = (/0.01, 0.05, 0.15, 0.5, 1., 5., 10., 50./)    !! Ice grid parameters 8 layers : thickness of ice levels
!$OMP THREADPRIVATE(zsicoef1)

  !
  ! diffuco.f90
  !

  ! 0. Constants

  REAL(r_std), PARAMETER :: Tetens_1 = 0.622         !! Ratio between molecular weight of water vapor and molecular weight  
                                                     !! of dry air (unitless)
  REAL(r_std), PARAMETER :: Tetens_2 = 0.378         !!
  REAL(r_std), PARAMETER :: ratio_H2O_to_CO2 = 1.6   !! Ratio of water vapor diffusivity to the CO2 diffusivity (unitless)
  REAL(r_std), PARAMETER :: mol_to_m_1 = 0.0244      !!
  REAL(r_std), PARAMETER :: RG_to_PAR = 0.5          !!
  REAL(r_std), PARAMETER :: W_to_mol = 4.6          !! W_to_mmol * RG_to_PAR = 2.3

  ! 1. Scalar

  LOGICAL, SAVE :: ldq_cdrag_from_gcm = .FALSE. !! Set to .TRUE. if you want q_cdrag coming from GCM
!$OMP THREADPRIVATE(ldq_cdrag_from_gcm)
  REAL(r_std), SAVE :: laimax = 12.             !! Maximal LAI used for splitting LAI into N layers (m^2.m^{-2})
!$OMP THREADPRIVATE(laimax)
  LOGICAL, SAVE :: downregulation_co2 = .TRUE.             !! Set to .TRUE. if you want CO2 downregulation.
!$OMP THREADPRIVATE(downregulation_co2)
  REAL(r_std), SAVE :: downregulation_co2_baselevel = 380. !! CO2 base level (ppm)
!$OMP THREADPRIVATE(downregulation_co2_baselevel)
  REAL(r_std), SAVE :: gb_ref = 1./25.                     !! Leaf bulk boundary layer resistance (s m-1)
!$OMP THREADPRIVATE(gb_ref)

  ! 3. Coefficients of equations

  REAL(r_std), SAVE :: lai_level_depth = 0.15  !!
!$OMP THREADPRIVATE(lai_level_depth)
!
  REAL(r_std), SAVE :: x1_coef =  0.177        !! Multiplicative factor for calculating the pseudo first order rate constant
                                               !! of assimilation response to
                                               !co2 kt (unitless)
!$OMP THREADPRIVATE(x1_coef)
  REAL(r_std), SAVE :: x1_Q10 =  0.069         !! Exponential factor in the equation defining kt (unitless)
!$OMP THREADPRIVATE(x1_Q10)
  REAL(r_std), SAVE :: quantum_yield =  0.092  !!
!$OMP THREADPRIVATE(quantum_yield)
  REAL(r_std), SAVE :: kt_coef = 0.7           !! Multiplicative factor in the equation defining kt (unitless)
!$OMP THREADPRIVATE(kt_coef)
  REAL(r_std), SAVE :: kc_coef = 39.09         !! Multiplicative factor for calculating the Michaelis-Menten
                                               !! coefficient Kc (unitless)
!$OMP THREADPRIVATE(kc_coef)
  REAL(r_std), SAVE :: Ko_Q10 = 0.085          !! Exponential factor for calculating the Michaelis-Menten coefficients
                                               !! Kc and Ko (unitless)
!$OMP THREADPRIVATE(Ko_Q10)
  REAL(r_std), SAVE :: Oa = 210000.            !! Intercellular concentration of O2 (ppm)
!$OMP THREADPRIVATE(Oa)
  REAL(r_std), SAVE :: Ko_coef =  2.412        !! Multiplicative factor for calculating the Michaelis-Menten 
                                               !! coefficient Ko (unitless)
!$OMP THREADPRIVATE(Ko_coef)
  REAL(r_std), SAVE :: CP_0 = 42.              !! Multiplicative factor for calculating the CO2 compensation 
                                               !! point CP (unitless)
!$OMP THREADPRIVATE(CP_0)
  REAL(r_std), SAVE :: CP_temp_coef = 9.46     !! Exponential factor for calculating the CO2 compensation point CP (unitless)
!$OMP THREADPRIVATE(CP_temp_coef)
  REAL(r_std), SAVE :: CP_temp_ref = 25.       !! Reference temperature for the CO2 compensation point CP (C)
!$OMP THREADPRIVATE(CP_temp_ref)
  !
  REAL(r_std), SAVE, DIMENSION(2) :: rt_coef = (/ 0.8, 1.3 /)    !!
!$OMP THREADPRIVATE(rt_coef)
  REAL(r_std), SAVE, DIMENSION(2) :: vc_coef = (/ 0.39, 0.3 /)   !!
!$OMP THREADPRIVATE(vc_coef)
  REAL(r_std), SAVE               :: c13_a = 4.4      !! fractionation against during diffusion
!$OMP THREADPRIVATE(c13_a)
  REAL(r_std), SAVE               :: c13_b = 27.      !! fractionation against during carboxylation
!$OMP THREADPRIVATE(c13_b)
  REAL(r_std), SAVE               :: threshold_c13_assim = 0.01 !! If assimilation falls below this threshold
                                                                !! the delta_c13
                                                                !is set to zero
!$OMP THREADPRIVATE(threshold_c13_assim)
  !
  REAL(r_std), SAVE, DIMENSION(6) :: dew_veg_poly_coeff = &            !! coefficients of the 5 degree polynomomial used
  & (/ 0.887773, 0.205673, 0.110112, 0.014843,  0.000824,  0.000017 /) !! in the equation of coeff_dew_veg
!$OMP THREADPRIVATE(dew_veg_poly_coeff)
!
  REAL(r_std), SAVE               :: Oi=210000.    !! Intercellular oxygen partial pressure (ubar)
!$OMP THREADPRIVATE(Oi)
  !
  ! slowproc.f90 
  !

  ! 1. Scalar

  INTEGER(i_std), SAVE :: ninput_year_orig = 0       !!  first year for N inputs (number)
!$OMP THREADPRIVATE(ninput_year_orig)
  LOGICAL, SAVE :: ninput_suffix_year = .FALSE.      !! Do the Ninput datasets have a 'year' suffix ? (y/n)  
!$OMP THREADPRIVATE(ninput_suffix_year)
  REAL(r_std), SAVE :: bulk_default = 1000.0           !! Default value for bulk density of soil (kg/m3)
!$OMP THREADPRIVATE(bulk_default)
  REAL(r_std), SAVE :: ph_default = 5.5              !! Default value for pH of soil (-)
!$OMP THREADPRIVATE(ph_default)

  REAL(r_std), SAVE :: min_vegfrac = 0.001           !! Minimal fraction of mesh a vegetation type can occupy (0-1, unitless)
!$OMP THREADPRIVATE(min_vegfrac)
  REAL(r_std), SAVE :: frac_nobio_fixed_test_1 = 0.0 !! Value for frac_nobio for tests in 0-dim simulations (0-1, unitless)
!$OMP THREADPRIVATE(frac_nobio_fixed_test_1)
  
  REAL(r_std), SAVE :: stempdiag_bid = 280.          !! only needed for an initial LAI if there is no restart file
!$OMP THREADPRIVATE(stempdiag_bid)


                           !-----------------------------!
                           !  STOMATE AND LPJ PARAMETERS !
                           !-----------------------------!


  ! Allometric relationships
  REAL(r_std), SAVE :: exp_kf = 3.0              !! Tuning parameter for the relationship between tree height and k_latosas
!$OMP THREADPRIVATE(exp_kf)
  !
  ! lpj_constraints.f90
  !
  
  ! 1. Scalar

  REAL(r_std), SAVE  :: too_long = 5.      !! longest sustainable time without 
                                           !! regeneration (vernalization) (years)
!$OMP THREADPRIVATE(too_long)


  !
  ! lpj_establish.f90
  !

  ! 1. Scalar

  REAL(r_std), SAVE :: estab_max_tree = 0.12   !! Maximum tree establishment rate (ind/m2/dt_stomate)
!$OMP THREADPRIVATE(estab_max_tree)
  REAL(r_std), SAVE :: estab_max_grass = 0.12  !! Maximum grass establishment rate (ind/m2/dt_stomate)
!$OMP THREADPRIVATE(estab_max_grass)
  
  ! 3. Coefficients of equations

  REAL(r_std), SAVE :: establish_scal_fact = 5.  !!
!$OMP THREADPRIVATE(establish_scal_fact)
  REAL(r_std), SAVE :: max_tree_coverage = 0.98  !! (0-1, unitless)
!$OMP THREADPRIVATE(max_tree_coverage)
  REAL(r_std), SAVE :: ind_0_estab = 0.2         !! = ind_0 * 10.
!$OMP THREADPRIVATE(ind_0_estab)


  !
  ! lpj_fire.f90
  !

  ! 1. Scalar

  REAL(r_std), SAVE :: tau_fire = 30.           !! Time scale for memory of the fire index (days).
!$OMP THREADPRIVATE(tau_fire)
  REAL(r_std), SAVE :: litter_crit = 200.       !! Critical litter quantity for fire
                                                !! below which iginitions extinguish 
                                                !! @tex $(gC m^{-2})$ @endtex
!$OMP THREADPRIVATE(litter_crit)
  REAL(r_std), SAVE :: fire_resist_lignin = 0.5 !!
!$OMP THREADPRIVATE(fire_resist_lignin)
  ! 2. Arrays

  REAL(r_std), SAVE, DIMENSION(nparts) :: co2frac = &    !! The fraction of the different biomass 
       & (/ .95, .95, 0., 0.3, 0., 0., .95, .95, .95 /)       !! compartments emitted to the atmosphere 
!$OMP THREADPRIVATE(co2frac)                                                         !! when burned (unitless, 0-1)  

  ! 3. Coefficients of equations

  REAL(r_std), SAVE, DIMENSION(3) :: bcfrac_coeff = (/ .3,  1.3,  88.2 /)         !! (unitless)
!$OMP THREADPRIVATE(bcfrac_coeff)
  REAL(r_std), SAVE, DIMENSION(4) :: firefrac_coeff = (/ 0.45, 0.8, 0.6, 0.13 /)  !! (unitless)
!$OMP THREADPRIVATE(firefrac_coeff)

  !
  ! lpj_gap.f90
  !

  ! 1. Scalar

  REAL(r_std), SAVE :: ref_greff = 0.035         !! Asymptotic maximum mortality rate
                                                 !! @tex $(year^{-1})$ @endtex
!$OMP THREADPRIVATE(ref_greff)

  !               
  ! lpj_light.f90 
  !              

  ! 1. Scalar
  
  LOGICAL, SAVE :: annual_increase = .TRUE. !! for diagnosis of fpc increase, compare today's fpc to last year's maximum (T) or
                                            !! to fpc of last time step (F)? (true/false)
!$OMP THREADPRIVATE(annual_increase)
  REAL(r_std), SAVE :: min_cover = 0.05     !! For trees, minimum fraction of crown area occupied
                                            !! (due to its branches etc.) (0-1, unitless)
                                            !! This means that only a small fraction of its crown area
                                            !! can be invaded by other trees.
!$OMP THREADPRIVATE(min_cover)
  !
  ! lpj_pftinout.f90 
  !

  ! 1. Scalar

  REAL(r_std), SAVE :: min_avail = 0.01         !! minimum availability
!$OMP THREADPRIVATE(min_avail)
  REAL(r_std), SAVE :: ind_0 = 0.02             !! initial density of individuals
!$OMP THREADPRIVATE(ind_0)
  ! 3. Coefficients of equations
  
  REAL(r_std), SAVE :: RIP_time_min = 1.25      !! test whether the PFT has been eliminated lately (years)
!$OMP THREADPRIVATE(RIP_time_min)
  REAL(r_std), SAVE :: npp_longterm_init = 10.  !! Initialisation value for npp_longterm (gC.m^{-2}.year^{-1})
!$OMP THREADPRIVATE(npp_longterm_init)
  REAL(r_std), SAVE :: everywhere_init = 0.05   !!
!$OMP THREADPRIVATE(everywhere_init)




  !
  ! stomate_data.f90 
  !

  ! 1. Scalar 

  ! 1.2 climatic parameters 

  REAL(r_std), SAVE :: precip_crit = 100.        !! minimum precip, in (mm/year)
!$OMP THREADPRIVATE(precip_crit)
  REAL(r_std), SAVE :: gdd_crit_estab = 150.     !! minimum gdd for establishment of saplings
!$OMP THREADPRIVATE(gdd_crit_estab)
  REAL(r_std), SAVE :: fpc_crit = 0.95           !! critical fpc, needed for light competition and establishment (0-1, unitless)
!$OMP THREADPRIVATE(fpc_crit)

  ! 1.3 sapling characteristics

  REAL(r_std), SAVE :: alpha_grass = 0.5         !! alpha coefficient for grasses (unitless)
!$OMP THREADPRIVATE(alpha_grass)
  REAL(r_std), SAVE :: alpha_tree = 1.           !! alpha coefficient for trees (unitless)
!$OMP THREADPRIVATE(alpha_tree)
  REAL(r_std), SAVE :: struct_to_leaves = 0.05  !! Fraction of structural carbon in grass and crops as a share of the leaf
                                                !! carbon pool. Only used for grasses and crops (thus NOT for trees)
	                                        !! (unitless)
!$OMP THREADPRIVATE(struct_to_leaves)

  REAL(r_std), SAVE :: labile_to_total = 0.01   !! Fraction of the labile pool in trees, grasses and crops as a share of the
                                                !! total carbon pool (accounting for the N-content of the different tissues).
                                                !! (unitless)
!$OMP THREADPRIVATE(labile_to_total)



  ! 1.4  time scales for phenology and other processes (in days)
REAL(r_std), SAVE :: tau_month = 20.              !! (days)       
!$OMP THREADPRIVATE(tau_month)
  REAL(r_std), SAVE :: tau_week = 7.              !! (days)  
!$OMP THREADPRIVATE(tau_week)
  REAL(r_std), SAVE :: tau_hum_month = 20.        !! (days)       
!$OMP THREADPRIVATE(tau_hum_month)
  REAL(r_std), SAVE :: tau_hum_week = 7.          !! (days)  
!$OMP THREADPRIVATE(tau_hum_week)
  REAL(r_std), SAVE :: tau_t2m_month = 20.        !! (days)      
!$OMP THREADPRIVATE(tau_t2m_month)
  REAL(r_std), SAVE :: tau_t2m_week = 7.          !! (days)  
!$OMP THREADPRIVATE(tau_t2m_week)
  REAL(r_std), SAVE :: tau_tsoil_month = 20.      !! (days)     
!$OMP THREADPRIVATE(tau_tsoil_month)
  REAL(r_std), SAVE :: tau_gpp_week = 7.          !! (days)  
!$OMP THREADPRIVATE(tau_gpp_week)
  REAL(r_std), SAVE :: tau_sugarload_week = 7.    !! (days)
!$OMP THREADPRIVATE(tau_sugarload_week)
  REAL(r_std), SAVE :: tau_gdd = 40.              !! (days)  
!$OMP THREADPRIVATE(tau_gdd)
  REAL(r_std), SAVE :: tau_ngd = 43.              !! (days)  
!$OMP THREADPRIVATE(tau_ngd)
  REAL(r_std), SAVE :: coeff_tau_longterm = 3.    !! (unitless)
!$OMP THREADPRIVATE(coeff_tau_longterm)
  REAL(r_std), SAVE :: tau_longterm_max           !! (days)  
!$OMP THREADPRIVATE(tau_longterm_max)

  ! 3. Coefficients of equations

  REAL(r_std), SAVE :: bm_sapl_carbres = 5.             !!
!$OMP THREADPRIVATE(bm_sapl_carbres)
  REAL(r_std), SAVE :: bm_sapl_labile = 5.             !!
!$OMP THREADPRIVATE(bm_sapl_labile)
  REAL(r_std), SAVE :: bm_sapl_sapabove = 0.5           !!
!$OMP THREADPRIVATE(bm_sapl_sapabove)
  REAL(r_std), SAVE :: bm_sapl_heartabove = 2.          !!
!$OMP THREADPRIVATE(bm_sapl_heartabove)
  REAL(r_std), SAVE :: bm_sapl_heartbelow = 2.          !!
!$OMP THREADPRIVATE(bm_sapl_heartbelow)
  REAL(r_std), SAVE :: init_sapl_mass_leaf_nat = 0.1    !!
!$OMP THREADPRIVATE(init_sapl_mass_leaf_nat)
  REAL(r_std), SAVE :: init_sapl_mass_leaf_agri = 1.    !!
!$OMP THREADPRIVATE(init_sapl_mass_leaf_agri)
  REAL(r_std), SAVE :: init_sapl_mass_carbres = 5.      !!
!$OMP THREADPRIVATE(init_sapl_mass_carbres)
  REAL(r_std), SAVE :: init_sapl_mass_labile = 5.      !!
!$OMP THREADPRIVATE(init_sapl_mass_labile)
  REAL(r_std), SAVE :: init_sapl_mass_root = 0.1        !!
!$OMP THREADPRIVATE(init_sapl_mass_root)
  REAL(r_std), SAVE :: init_sapl_mass_fruit = 0.3       !!  
!$OMP THREADPRIVATE(init_sapl_mass_fruit)
  REAL(r_std), SAVE :: cn_sapl_init = 0.5               !!
!$OMP THREADPRIVATE(cn_sapl_init)
  REAL(r_std), SAVE :: migrate_tree = 10.*1.E3          !!
!$OMP THREADPRIVATE(migrate_tree)
  REAL(r_std), SAVE :: migrate_grass = 10.*1.E3         !!
!$OMP THREADPRIVATE(migrate_grass)
  REAL(r_std), SAVE :: lai_initmin_tree = 0.3           !!
!$OMP THREADPRIVATE(lai_initmin_tree)
  REAL(r_std), SAVE :: lai_initmin_grass = 0.1          !!
!$OMP THREADPRIVATE(lai_initmin_grass)
  REAL(r_std), SAVE, DIMENSION(2) :: dia_coeff = (/ 4., 0.5 /)            !!
!$OMP THREADPRIVATE(dia_coeff)
  REAL(r_std), SAVE, DIMENSION(4) :: bm_sapl_leaf = (/ 4., 4., 0.8, 5./)  !!
!$OMP THREADPRIVATE(bm_sapl_leaf)


  !
  ! stomate_litter.f90 
  !

  ! 0. Constants

  REAL(r_std), PARAMETER :: Q10 = 10.               !!

  ! 1. Scalar

  REAL(r_std), SAVE :: z_decomp = 0.2               !!  Maximum depth for soil decomposer's activity (m)
!$OMP THREADPRIVATE(z_decomp)

  ! 2. Arrays

  REAL(r_std), SAVE :: frac_soil_struct_sua = 0.4    !! corresponding to frac_soil(istructural,isurface,iabove) 
!$OMP THREADPRIVATE(frac_soil_struct_sua)
  REAL(r_std), SAVE :: frac_soil_struct_aa = 0.4     !! corresponding to frac_soil(istructural,iactive,iabove) 
!$OMP THREADPRIVATE(frac_soil_struct_aa)
 REAL(r_std), SAVE :: frac_soil_struct_ab = 0.45   !! corresponding to frac_soil(istructural,iactive,ibelow)
!$OMP THREADPRIVATE(frac_soil_struct_ab)
  REAL(r_std), SAVE :: frac_soil_struct_sa = 0.7    !! corresponding to frac_soil(istructural,islow,iabove)
!$OMP THREADPRIVATE(frac_soil_struct_sa)
  REAL(r_std), SAVE :: frac_soil_struct_sb = 0.7    !! corresponding to frac_soil(istructural,islow,ibelow)
!$OMP THREADPRIVATE(frac_soil_struct_sb)
  REAL(r_std), SAVE :: frac_soil_metab_sua = 0.4    !! corresponding to frac_soil(imetabolic,isurface,iabove)
!$OMP THREADPRIVATE(frac_soil_metab_sua)
  REAL(r_std), SAVE :: frac_soil_metab_aa = 0.4    !! corresponding to frac_soil(imetabolic,iactive,iabove)
!$OMP THREADPRIVATE(frac_soil_metab_aa)
  REAL(r_std), SAVE :: frac_soil_metab_ab = 0.45    !! corresponding to frac_soil(imetabolic,iactive,ibelow)
!$OMP THREADPRIVATE(frac_soil_metab_ab)
  REAL(r_std), SAVE :: frac_woody = 0.65   !! Coefficient for determining the lignin fraction of woody litter
!$OMP THREADPRIVATE(frac_woody)
  REAL(r_std), SAVE :: frac_snag = 0.65   !! Coefficient for determining the lignin fraction of snag litter (unitless)
!$OMP THREADPRIVATE(frac_snag)
  REAL(r_std), SAVE, DIMENSION(nparts) :: CN_fix = & !! C/N ratio of each plant pool (0-100, unitless)
       & (/ 40., 40., 40., 40., 40., 40., 40., 40., 40. /) 
!$OMP THREADPRIVATE(CN_fix)

  ! 3. Coefficients of equations

  REAL(r_std), SAVE :: metabolic_ref_frac = 0.85    !! used by litter and soilcarbon (0-1, unitless)
!$OMP THREADPRIVATE(metabolic_ref_frac)
  REAL(r_std), SAVE :: metabolic_LN_ratio = 0.018   !! (0-1, unitless)   
!$OMP THREADPRIVATE(metabolic_LN_ratio)
  ! Turnover rate (yr-1) - From Parton et al., 1993
  REAL(r_std), SAVE :: turn_metabolic = 15           !!
!$OMP THREADPRIVATE(turn_metabolic)
  REAL(r_std), SAVE :: turn_struct = 4                !!
!$OMP THREADPRIVATE(turn_struct)
  REAL(r_std), SAVE :: turn_woody = 1.33              !! from DOFOCO
!$OMP THREADPRIVATE(turn_woody)
  REAL(r_std), SAVE :: turn_snag = 0.00133             !! Decomposition / Turnover of snag pool to soil organic matter
                                                       !! [years]
!$OMP THREADPRIVATE(turn_snag)
  REAL(r_std), SAVE :: turn_snag_fall = 0.133          !! Snag fall / Turnover of snag pool to woody litter pool
                                                       !! [years]
!$OMP THREADPRIVATE(turn_snag_fall)
  REAL(r_std), SAVE :: soil_Q10 = 0.69              !!= ln 2
!$OMP THREADPRIVATE(soil_Q10)
  REAL(r_std), SAVE :: soil_Q10_uptake = 0.45         
!$OMP THREADPRIVATE(soil_Q10_uptake)
  REAL(r_std), SAVE :: tsoil_ref = 30.              !!
!$OMP THREADPRIVATE(tsoil_ref)
  REAL(r_std), SAVE :: litter_struct_coef = 3.      !! 
!$OMP THREADPRIVATE(litter_struct_coef)
  REAL(r_std), SAVE, DIMENSION(3) :: moist_coeff = (/ 1.1,  2.4,  0.29 /) !!
!$OMP THREADPRIVATE(moist_coeff)
  REAL(r_std), SAVE :: moistcont_min = 0.25  !! minimum soil wetness to limit the heterotrophic respiration
!$OMP THREADPRIVATE(moistcont_min)
  REAL(r_std), SAVE :: finerootdepthratio = 0.5  !! the ratio of fine root to overall root e-folding depth (for C inputs)
!$OMP THREADPRIVATE(finerootdepthratio)
  REAL(r_std), SAVE :: altrootratio = 0.5        !! the maximum ratio of fine root depth to active layer thickness (for C inputs)
!$OMP THREADPRIVATE(altrootratio)

  !
  ! stomate_lpj.f90
  !

  ! 1. Scalar

  REAL(r_std), SAVE :: frac_turnover_daily = 0.55  !! (0-1, unitless)
!$OMP THREADPRIVATE(frac_turnover_daily)


  !
  ! stomate_npp.f90 
  !

  ! 1. Scalar

  REAL(r_std), SAVE :: tax_max = 0.8 !! Maximum fraction of allocatable biomass used 
                                     !! for maintenance respiration (0-1, unitless)
!$OMP THREADPRIVATE(tax_max)


  !
  ! stomate_phenology.f90
  !

  ! 1. Scalar

  REAL(r_std), SAVE :: min_growthinit_time = 300.  !! minimum time since last beginning of a growing season (days)
!$OMP THREADPRIVATE(min_growthinit_time)

  ! 3. Coefficients of equations
  
  REAL(r_std), SAVE :: gddncd_ref = 964.           !! Empirical parameter to determine budbreak for ncdgdd phenology type. Reprensents gdd.
!$OMP THREADPRIVATE(gddncd_ref)
  REAL(r_std), SAVE :: gddncd_curve = 0.0058       !! Empirical parameter to determine budbreak for ncdgdd phenology type
!$OMP THREADPRIVATE(gddncd_curve)
  REAL(r_std), SAVE :: gddncd_offset = 12.8        !! Empirical parameter to determine budbreak for ncdgdd phenology type
!$OMP THREADPRIVATE(gddncd_offset)


  !
  ! stomate_resp.f90
  !

  ! 3. Coefficients of equations

  REAL(r_std), SAVE :: maint_resp_min_vmax = 0.3   !!
!$OMP THREADPRIVATE(maint_resp_min_vmax)
  REAL(r_std), SAVE :: maint_resp_coeff = 1.4      !!
!$OMP THREADPRIVATE(maint_resp_coeff)


  !
  ! stomate_som_dynamics.f90 (in stomate_soilcarbon.f90 or
  ! stomate_soil_carbon_discretization.f90)   
  !

  ! 2. Arrays 

  ! 2.1 Fixed fraction from one pool to another (or to CO2 emission)

  REAL(r_std), SAVE :: active_to_pass_ref_frac = 0.003  !! from active pool: depends on clay content  (0-1, unitless)
                                                        !! corresponding to frac_carb(:,iactive,ipassive)
!$OMP THREADPRIVATE(active_to_pass_ref_frac)
  REAL(r_std), SAVE :: surf_to_slow_ref_frac = 0.4      !! from surface pool
                                                        !! corresponding to frac_carb(:,isurf,islow)
!$OMP THREADPRIVATE(surf_to_slow_ref_frac)

  REAL(r_std), SAVE :: active_to_CO2_ref_frac  = 0.85   !! from active pool: depends on clay content  (0-1, unitless)
                                                        !! corresponding to frac_resp(:,iactive)
!$OMP THREADPRIVATE(active_to_CO2_ref_frac)
  REAL(r_std), SAVE :: slow_to_pass_ref_frac   = 0.03  !! from slow pool: depends on clay content  (0-1, unitless) 
                                                        !! corresponding to frac_carb(:,islow,ipassive)
!$OMP THREADPRIVATE(slow_to_pass_ref_frac)
  REAL(r_std), SAVE :: slow_to_CO2_ref_frac    = 0.55   !! from slow pool (0-1, unitless) 
                                                        !! corresponding to frac_resp(:,islow)
!$OMP THREADPRIVATE(slow_to_CO2_ref_frac)
  REAL(r_std), SAVE :: pass_to_active_ref_frac = 0.45   !! from passive pool (0-1, unitless)
                                                        !! corresponding to frac_carb(:,ipassive,iactive)
!$OMP THREADPRIVATE(pass_to_active_ref_frac)
  REAL(r_std), SAVE :: pass_to_slow_ref_frac   = 0.0    !! from passive pool (0-1, unitless)
                                                        !! corresponding to frac_carb(:,ipassive,islow)
!$OMP THREADPRIVATE(pass_to_slow_ref_frac)

  ! 2.2 som carbon pools

  REAL(r_std), SAVE :: som_init_active = 1000           !! Initial active SOM carbon (g m-2) 
!$OMP THREADPRIVATE(som_init_active)
  REAL(r_std), SAVE :: som_init_slow = 3000             !! Initial slow SOM carbon (g m-2)
!$OMP THREADPRIVATE(som_init_slow)
  REAL(r_std), SAVE :: som_init_passive = 5000          !! Initial passive SOM carbon (g m-2)
!$OMP THREADPRIVATE(som_init_passive)
  REAL(r_std), SAVE :: som_init_surface = 1000          !! Initial surface SOM carbon (g m-2)
!$OMP THREADPRIVATE(som_init_surface)

  ! 3. Define Variable fraction from one pool to another (function of silt and clay fraction)
  REAL(r_std), SAVE :: active_to_pass_clay_frac     = 0.032
!$OMP THREADPRIVATE(active_to_pass_clay_frac)
  REAL(r_std), SAVE :: active_to_CO2_clay_silt_frac = 0.68
!$OMP THREADPRIVATE(active_to_CO2_clay_silt_frac)
  REAL(r_std), SAVE :: slow_to_pass_clay_frac   = 0 
!$OMP THREADPRIVATE(slow_to_pass_clay_frac)

  ! C to N target ratios of differnt pools
  REAL(r_std), SAVE ::  CN_target_iactive_ref  = 15. !! CN target ratio of active pool for soil min N = 0
!$OMP THREADPRIVATE(CN_target_iactive_ref)
  REAL(r_std), SAVE ::  CN_target_islow_ref    = 20. !! CN target ratio of slow pool for soil min N = 0
!$OMP THREADPRIVATE(CN_target_islow_ref)
  REAL(r_std), SAVE ::  CN_target_ipassive_ref = 10. !! CN target ratio of passive pool for soil min N = 0
!$OMP THREADPRIVATE(CN_target_ipassive_ref)
  REAL(r_std), SAVE ::  CN_target_isurface_ref = 20. !! CN target ratio of surface pool for litter nitrogen content = 0
!$OMP THREADPRIVATE(CN_target_isurface_ref)
  REAL(r_std), SAVE ::  CN_target_iactive_Nmin  = -6.  !! CN target ratio change per mineral N unit (g m-2) for active pool
!$OMP THREADPRIVATE(CN_target_iactive_Nmin) 
  REAL(r_std), SAVE ::  CN_target_islow_Nmin    = -4.  !! CN target ratio change per mineral N unit (g m-2) for slow pool 
!$OMP THREADPRIVATE(CN_target_islow_Nmin)
  REAL(r_std), SAVE ::  CN_target_ipassive_Nmin = -1.5 !! CN target ratio change per mineral N unit (g m-2) for passive pool 
!$OMP THREADPRIVATE(CN_target_ipassive_Nmin)
  REAL(r_std), SAVE ::  CN_target_isurface_pnc  = -5.  !! CN target ratio change per plant nitrogen content unit (%) for surface pool
  !! Turnover in SOM pools (year-1)
!$OMP THREADPRIVATE(CN_target_isurface_pnc)
  REAL(r_std), SAVE :: som_turn_isurface = 6.0           !! turnover of surface pool (year-1)
!$OMP THREADPRIVATE(som_turn_isurface)
  REAL(r_std), SAVE :: som_turn_iactive     = 7.3        !! turnover of active pool (year-1)
!$OMP THREADPRIVATE(som_turn_iactive)
  REAL(r_std), SAVE :: som_turn_islow       = 0.1        !! turnover of slow pool (year-1)
!$OMP THREADPRIVATE(som_turn_islow)
  REAL(r_std), SAVE :: som_turn_ipassive    = 0.0025     !! turnover of passive pool (year-1)
!$OMP THREADPRIVATE(som_turn_ipassive)
  REAL(r_std), SAVE :: fslow = 37.0                      !! convertiing factor to go from active pool turnover to slow pool turnover from Guimberteau et al 2018 GMD
!$OMP THREADPRIVATE(fslow)
  REAL(r_std), SAVE :: fpassive = 1617.45                !! convertiing factor to go from active pool turnover to passive pool turnover from Guimberteau et al 2018 GMD
!$OMP THREADPRIVATE(fpassive)
  REAL(r_std), SAVE :: stomate_tau = 4.699E6 
!$OMP THREADPRIVATE(stomate_tau)
  REAL(r_std), SAVE :: depth_modifier = 1.E6             !! e-folding depth of turnover rates,following Koven et al.,2013,Biogeosciences. A very large value means no depth modification
!$OMP THREADPRIVATE(depth_modifier)
  REAL(r_std), SAVE :: som_turn_iactive_clay_frac = 0.75 !! clay-dependant parameter impacting on turnover rate of active pool 
                                                         !! Tm parameter of Parton et al. 1993 (-)
!$OMP THREADPRIVATE(som_turn_iactive_clay_frac)

  !
  ! stomate_turnover.f90
  !

  ! 3. Coefficients of equations

  REAL(r_std), SAVE :: new_turnover_time_ref = 20. !!(days)
!$OMP THREADPRIVATE(new_turnover_time_ref)


  !
  ! stomate_vmax.f90
  !
 
  ! 1. Scalar

  REAL(r_std), SAVE :: vmax_offset = 0.3        !! minimum leaf efficiency (unitless)
!$OMP THREADPRIVATE(vmax_offset)
  REAL(r_std), SAVE :: leafage_firstmax = 0.03  !! relative leaf age at which efficiency
                                                !! reaches 1 (unitless)
!$OMP THREADPRIVATE(leafage_firstmax)
  REAL(r_std), SAVE :: leafage_lastmax = 0.5    !! relative leaf age at which efficiency
                                                !! falls below 1 (unitless)
!$OMP THREADPRIVATE(leafage_lastmax)
  REAL(r_std), SAVE :: leafage_old = 1.         !! relative leaf age at which efficiency
                                                !! reaches its minimum (vmax_offset) 
                                                !! (unitless)
!$OMP THREADPRIVATE(leafage_old)
  REAL(r_std), SAVE :: sugar_load_min = 0.05    !! Lower bound for sugar loading when used to regulate NUE. 
                                                !! Sugar loafding results in a strong reduction of GPP. By
                                                !! taking a rather high lower limit the impact of sugar loading
                                                !! is limited. This will result in high C-reserves unless 
                                                !! leaching is implemented. 
!$OMP THREADPRIVATE(sugar_load_min)
  REAL(r_std), SAVE :: sugar_load_max = 1.0     !! Upper bound for sugar loading when used to regulate NUE
!$OMP THREADPRIVATE(sugar_load_max)

  !
  ! nitrogen_dynamics (in stomate_soilcarbon.f90) 
  !

  ! 0. Constants
  REAL(r_std), PARAMETER :: D_air = 1.73664     !! Oxygen diffusion rate in the air = 0.07236 m2/h
                                               !! from Table 2 of Li et al, 2000
                                               !! (m**2/day)

  REAL(r_std), PARAMETER :: C_molar_mass = 12  !! Carbon Molar mass (gC mol-1) 

  REAL(r_std), PARAMETER :: Pa_to_hPa    = 0.01      !! Conversion factor from Pa to hPa (-)
  REAL(r_std), PARAMETER :: V_O2         = 0.209476  !! Volumetric fraction of O2 in air (-)

  REAL(r_std), PARAMETER :: pk_NH4 = 9.25      !! The negative logarithm of the acid dissociation constant K_NH4     
                                               !! See Table 4 of Li et al. 1992 and Appendix A of Zhang et al. 2002    

                                     
  ! 1. Scalar

  ! Coefficients for defining maximum porosity
  ! From Saxton, K.E., Rawls, W.J., Romberger, J.S., Papendick, R.I., 1986
  ! Estimationg generalized soil-water characteristics from texture. 
  ! Soil Sci. Soc. Am. J. 50, 1031-1036
  ! Cited in Table 5 (page 444) of
  ! Y. Pachepsky, W.J. Rawls
  ! Development of Pedotransfer Functions in Soil Hydrology
  ! Elsevier, 23 nov. 2004 - 542 pages
  ! http://books.google.fr/books?id=ar_lPXaJ8QkC&printsec=frontcover&hl=fr#v=onepage&q&f=false
  REAL(r_std), SAVE :: h_saxton = 0.332          !! h coefficient 
!$OMP THREADPRIVATE(h_saxton)
  REAL(r_std), SAVE :: j_saxton = -7.251*1e-4    !! j coefficient 
!$OMP THREADPRIVATE(j_saxton)
  REAL(r_std), SAVE :: k_saxton = 0.1276         !! k coefficient
!$OMP THREADPRIVATE(k_saxton)

  ! Values of the power used in the equation defining the diffusion of oxygen in soil
  ! from Table 2 of Li et al, 2000
  REAL(r_std), SAVE :: diffusionO2_power_1 = 3.33 !! (unitless)
!$OMP THREADPRIVATE(diffusionO2_power_1)
  REAL(r_std), SAVE :: diffusionO2_power_2 = 2.0  !! (unitless) 
!$OMP THREADPRIVATE(diffusionO2_power_2)

  ! Temperature-related Factors impacting on Oxygen diffusion rate
  ! From eq. 2 of Table 2 (Li et al, 2000)
  REAL(r_std), SAVE ::   F_nofrost = 1.2          !! (unitless)
!$OMP THREADPRIVATE(F_nofrost)
  REAL(r_std), SAVE ::   F_frost   = 0.8          !! (unitless)
!$OMP THREADPRIVATE(F_frost)

  ! Coefficients used in the calculation of Volumetric fraction of anaerobic microsites 
  ! a and b constants are not specified in Li et al., 2000
  ! S. Zaehle used a=0.85 and b=1 without mention to any publication
  REAL(r_std), SAVE ::   a_anvf    = 0.85 !! (-)
!$OMP THREADPRIVATE(a_anvf)
  REAL(r_std), SAVE ::   b_anvf    = 1.   !! (-)
!$OMP THREADPRIVATE(b_anvf)

  ! Coefficients used in the calculation of the Fraction of adsorbed NH4+
  ! Li et al. 1992, JGR, Table 4
  REAL(r_std), SAVE ::   a_FixNH4 = 0.41  !! (-)
!$OMP THREADPRIVATE(a_FixNH4)
  REAL(r_std), SAVE ::   b_FixNH4 = -0.47 !! (-)
!$OMP THREADPRIVATE(b_FixNH4)
  REAL(r_std), SAVE ::   clay_max = 0.63  !! (-)
!$OMP THREADPRIVATE(clay_max)

  ! Coefficients used in the calculation of the Response of Nitrification
  ! to soil moisture
  ! Zhang et al. 2002, Ecological Modelling, appendix A, page 101
  REAL(r_std), SAVE ::   fw_nit_0 =  -0.0243  !! (-)
!$OMP THREADPRIVATE(fw_nit_0)
  REAL(r_std), SAVE ::   fw_nit_1 =   0.9975  !! (-)
!$OMP THREADPRIVATE(fw_nit_1)
  REAL(r_std), SAVE ::   fw_nit_2 =  -5.5368  !! (-)
!$OMP THREADPRIVATE(fw_nit_2)
  REAL(r_std), SAVE ::   fw_nit_3 =  17.651   !! (-)
!$OMP THREADPRIVATE(fw_nit_3)
  REAL(r_std), SAVE ::   fw_nit_4 = -12.904   !! (-)
!$OMP THREADPRIVATE(fw_nit_4)

  ! Coefficients used in the calculation of the Response of Nitrification
  ! to Temperature
  ! Zhang et al. 2002, Ecological Modelling, appendix A, page 101
  REAL(r_std), SAVE ::   ft_nit_0 =  -0.0233 !! (-)
!$OMP THREADPRIVATE(ft_nit_0)
  REAL(r_std), SAVE ::   ft_nit_1 =   0.3094 !! (-)
!$OMP THREADPRIVATE(ft_nit_1)
  REAL(r_std), SAVE ::   ft_nit_2 =  -0.2234 !! (-)
!$OMP THREADPRIVATE(ft_nit_2)
  REAL(r_std), SAVE ::   ft_nit_3 =   0.1566 !! (-)
!$OMP THREADPRIVATE(ft_nit_3)
  REAL(r_std), SAVE ::   ft_nit_4 =  -0.0272 !! (-)
!$OMP THREADPRIVATE(ft_nit_4)

  ! Coefficients used in the calculation of the Response of Nitrification
  ! to pH
  ! Zhang et al. 2002, Ecological Modelling, appendix A, page 101
  REAL(r_std), SAVE ::   fph_0 = -1.2314  !! (-)
!$OMP THREADPRIVATE(fph_0)
  REAL(r_std), SAVE ::   fph_1 = 0.7347   !! (-)
!$OMP THREADPRIVATE(fph_1)
  REAL(r_std), SAVE ::   fph_2 = -0.0604  !! (-)
!$OMP THREADPRIVATE(fph_2)

  ! Coefficients used in the calculation of the response of NO2 or NO 
  ! production during nitrificationof to Temperature
  ! Zhang et al. 2002, Ecological Modelling, appendix A, page 102
  REAL(r_std), SAVE ::   ftv_0 = 2.72   !! (-)
!$OMP THREADPRIVATE(ftv_0)
  REAL(r_std), SAVE ::   ftv_1 = 34.6   !! (-)
!$OMP THREADPRIVATE(ftv_1)
  REAL(r_std), SAVE ::   ftv_2 = 9615.  !! (-) 
!$OMP THREADPRIVATE(ftv_2)

  REAL(r_std), SAVE ::   k_nitrif = 2.0         !! Nitrification rate at 20 ◦C and field capacity (day-1)
                                                !! Schmid et al., 2001 give a value of 0.2 per day.
                                                !! (https://doi.org/10.1023/A:1012694218748)
                                                !! OCN used a value of 2.0.  We keep the OCN value.
!$OMP THREADPRIVATE(k_nitrif)

  REAL(r_std), SAVE ::   n2o_nitrif_p = 0.0006  !! Reference n2o production per N-NO3 produced g N-N2O  (g N-NO3)-1
                                                !! From Zhang et al., 2002 - Appendix A p. 102
!$OMP THREADPRIVATE(n2o_nitrif_p)
  REAL(r_std), SAVE ::   no_nitrif_p = 0.0025   !! Reference NO production per N-NO3 produced g N-NO  (g N-NO3)-1
                                                !! From Zhang et al., 2002 - Appendix A p. 102
!$OMP THREADPRIVATE(no_nitrif_p)

  ! NO production from chemodenitrification
  ! based on Kesik et al., 2005, Biogeosciences
  ! Coefficients used in the calculation of the Response to Temperature
  REAL(r_std), SAVE ::   chemo_t0  = -31494. !! (-)
!$OMP THREADPRIVATE(chemo_t0)
  ! Coefficients use in the calculation of the Response to pH
  REAL(r_std), SAVE ::   chemo_ph0 = -1.62   !! (-)
!$OMP THREADPRIVATE(chemo_ph0)
  ! Coefficients used in the calculation of NO production from chemodenitrification
  REAL(r_std), SAVE ::   chemo_0   = 30.     !! (-)
!$OMP THREADPRIVATE(chemo_0)
  REAL(r_std), SAVE ::   chemo_1   = 16565.  !! (-)
!$OMP THREADPRIVATE(chemo_1)

  ! Denitrification processes
  ! Li et al, 2000, JGR Table 4 eq 1, 2 and 4
  !
  ! Coefficients used in the Temperature response of 
  ! relative growth rate of total denitrifiers - Eq. 2 Table 4 of Li et al., 2000
  REAL(r_std), SAVE ::   ft_denit_0 = 2.     !! (-)
!$OMP THREADPRIVATE(ft_denit_0)
  REAL(r_std), SAVE ::   ft_denit_1 = 22.5   !! (-)
!$OMP THREADPRIVATE(ft_denit_1)
  REAL(r_std), SAVE ::   ft_denit_2 = 10.    !! (-)
!$OMP THREADPRIVATE(ft_denit_2)
  !
  ! Coefficients used in the pH response of 
  ! relative growth rate of total denitrifiers - Eq. 2 Table 4 of Li et al., 2000
  REAL(r_std), SAVE ::   fph_no3_0  = 4.25    !! (-) 
!$OMP THREADPRIVATE(fph_no3_0)
  REAL(r_std), SAVE ::   fph_no3_1  = 0.5     !! (-)
!$OMP THREADPRIVATE(fph_no3_1)
  REAL(r_std), SAVE ::   fph_no_0  = 5.25     !! (-)
!$OMP THREADPRIVATE(fph_no_0)
  REAL(r_std), SAVE ::   fph_no_1  = 1.       !! (-)
!$OMP THREADPRIVATE(fph_no_1)
  REAL(r_std), SAVE ::   fph_n2o_0  = 6.25    !! (-)
!$OMP THREADPRIVATE(fph_n2o_0)
  REAL(r_std), SAVE ::   fph_n2o_1  = 1.5     !! (-)
!$OMP THREADPRIVATE(fph_n2o_1)

  REAL(r_std), SAVE ::   Kn = 0.083           !! Half Saturation of N oxydes (kgN/m3)
                                              !! Table 4 of Li et al., 2000
!$OMP THREADPRIVATE(Kn)

  REAL(r_std), SAVE ::   cte_bact = 0.00015
!$OMP THREADPRIVATE(cte_bact)

  ! Maximum Relative growth rate of Nox denitrifiers
  ! Eq.1 Table 4 Li et al., 2000 
  REAL(r_std), SAVE ::   mu_no3_max = 0.67   !! (hour-1)
!$OMP THREADPRIVATE(mu_no3_max)
  REAL(r_std), SAVE ::   mu_no_max  = 0.67   !! (hour-1)
!$OMP THREADPRIVATE(mu_no_max)
  REAL(r_std), SAVE ::   mu_n2o_max = 0.67   !! (hour-1)
!$OMP THREADPRIVATE(mu_n2o_max)

  ! Maximum growth yield of NOx denitrifiers on N oxydes
  ! Table 4 Li et al., 2000
  REAL(r_std), SAVE ::   Y_no3 = 0.401 !! (kgC / kgN)
!$OMP THREADPRIVATE(Y_no3)
  REAL(r_std), SAVE ::   Y_no  = 0.428 !! (kgC / kgN)
!$OMP THREADPRIVATE(Y_no)
  REAL(r_std), SAVE ::   Y_n2o = 0.151 !! (kgC / kgN)
!$OMP THREADPRIVATE(Y_n2o)

  ! Maintenance coefficient on N oxyde
  ! Table 4 Li et al., 2000
  REAL(r_std), SAVE ::   M_no3 = 0.09   !! (kgN / kgC / hour)
!$OMP THREADPRIVATE(M_no3)
  REAL(r_std), SAVE ::   M_no  = 0.035  !! (kgN / kgC / hour)
!$OMP THREADPRIVATE(M_no)
  REAL(r_std), SAVE ::   M_n2o = 0.079  !! (kgN / kgC / hour)
!$OMP THREADPRIVATE(M_n2o)

        
  REAL(r_std), SAVE ::   Maint_c = 0.0076    !! Maintenance coefficient of carbon (kgC/kgC/h)
                                        !! Table 4 Li et al., 2000
!$OMP THREADPRIVATE(Maint_c)
  REAL(r_std), SAVE ::   Yc = 0.503     !! Maximum growth yield on soluble carbon (kgC/kgC)
                                        !! Table 4 Li et al., 2000
!$OMP THREADPRIVATE(Yc)

  !! Coefficients used in the eq. defining the response of N-emission to clay fraction (-)
  !! from  Table 4, Li et al. 2000
  REAL(r_std), SAVE ::   F_clay_0 = 0.13   
!$OMP THREADPRIVATE(F_clay_0) 
  REAL(r_std), SAVE ::   F_clay_1 = -0.079
!$OMP THREADPRIVATE(F_clay_1)


  REAL(r_std), SAVE ::   ratio_nh4_fert = 0.5  !! Proportion of ammonium in the fertilizers (ammo-nitrate) 
                                           
!$OMP THREADPRIVATE(ratio_nh4_fert)

  REAL(r_std), SAVE ::   cn_ratio_manure = 13.7  !! C:N ratio of organic fertilizer (average over table 3-1-1 from Fuchs et al. 2014) 
 
!$OMP THREADPRIVATE(cn_ratio_manure)

  ! 2. Arrays
  REAL(r_std), SAVE, DIMENSION(2)  :: K_N_min = (/ 30., 30. /)           !! [NH4+] (resp. [NO3-]) for which the Nuptake 
                                                                         !! equals vmax/2.   (umol per litter)
                                                                         !! from Kronzucker, 1995
!$OMP THREADPRIVATE(K_N_min)

  REAL(r_std), SAVE, DIMENSION(2)  :: low_K_N_min = (/ 0.0002, 0.0002 /) !! Rate of N uptake not associated with 
                                                                         !! Michaelis- Menten Kinetics for Ammonium 
                                                                         !! (ind.1) and Nitrate (ind.2)
                                                                         !! from Kronzucker, 1995 ((umol)-1)
!$OMP THREADPRIVATE(low_K_N_min)

  REAL(r_std), SAVE      :: emm_fac = 0.2         !! Factor for reducing NH3 emission  
!$OMP THREADPRIVATE(emm_fac)
  REAL(r_std), SAVE      :: fact_kn_no = 0.012    !! Factor for adusting kn constant for NOx production
!$OMP THREADPRIVATE(fact_kn_no)
  REAL(r_std), SAVE      :: fact_kn_n2o = 0.04    !! Factor for adusting kn constant for N2O production
!$OMP THREADPRIVATE(fact_kn_n2o)
  REAL(r_std), SAVE      :: kfwdenit = -5.        !! Factor for adjusting sensitivity of denitrification to water content                     
!$OMP THREADPRIVATE(kfwdenit)  
 REAL(r_std), SAVE      :: fwdenitfc = 0.05       !! Value at field capacity of the sensitivity function of denitrification to water content                     
!$OMP THREADPRIVATE(fwdenitfc)          
  REAL(r_std), SAVE      :: fracn_drainage = 1.0  !! Fraction of NH3/NO3 loss by drainage 
!$OMP THREADPRIVATE(fracn_drainage)
  REAL(r_std), SAVE      :: fracn_runoff = 1.0    !! Fraction of NH3/NO3 loss by runoff
!$OMP THREADPRIVATE(fracn_runoff)


  !! Other N-related parameters
  REAL(r_std), SAVE      :: Dmax = 0.25           !! Maximal elasticity of foliage N concentrations (Zaehle et al 2010, SI2, Table 1)
!$OMP THREADPRIVATE(Dmax)

  REAL(r_std), SAVE :: reserve_time_tree = 30.     !! Maximum number of days during which
                                                   !! carbohydrate reserve may be used for 
                                                   !! trees (days)
!$OMP THREADPRIVATE(reserve_time_tree)
  
  REAL(r_std), SAVE :: reserve_time_grass = 20.    !! Maximum number of days during which
                                                   !! carbohydrate reserve may be used for 
                                                   !! grasses (days)
!$OMP THREADPRIVATE(reserve_time_grass)
  REAL(r_std), SAVE :: p_n_uptake = 0.6            !! The minimum correction factor for nitrogen uptake in case of 
                                                   !! enough nitrogen in the reserve 
!$OMP THREADPRIVATE(p_n_uptake)


  !! Scalars for slowproc_impose_site_level_lai (read_lai == 1)
  REAL(r_std), SAVE      :: max_canopy_height_prescribed = 50.0         !! Maximum height of the canopy
!$OMP THREADPRIVATE(max_canopy_height_prescribed)
  REAL(r_std), SAVE      :: min_canopy_height_prescribed = 10.0         !! Starting height of the canopy
!$OMP THREADPRIVATE(min_canopy_height_prescribed)
  REAL(r_std), SAVE      :: canopy_density_prescribed = 0.15            !! Prescribed minimum LAI 
!$OMP THREADPRIVATE(canopy_density_prescribed)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: lvllai_prescribed     !! Prescribed LAI per level
!$OMP THREADPRIVATE(lvllai_prescribed)
  REAL(r_std), SAVE      :: maximum_lai = 5.0                           !! Prescribed maximum LAI 
!$OMP THREADPRIVATE(maximum_lai)
  REAL(r_std), SAVE      :: minimum_lai = 0.0                           !! Prescribed minimum LAI 
!$OMP THREADPRIVATE(minimum_lai)
  INTEGER(i_std), SAVE   :: phenology_doy_start = 1                     !! Day of Year at which phenology starts (leaves appear)
!$OMP THREADPRIVATE(phenology_doy_start)
  INTEGER(i_std), SAVE   :: senescence_doy_end = 365                    !! Day of Year at which senescence ends (leaves fall)
                                                                        !!(between 1 and 365, higher than phenology_doy_start)
!$OMP THREADPRIVATE(senescence_doy_end)
  INTEGER(i_std), SAVE   :: phenology_delta = 30                        !! Delta of Doy between start and end of phenology
!$OMP THREADPRIVATE(phenology_delta)
  INTEGER(i_std), SAVE   :: senescence_delta = 30                       !! Delta of Doy between start and end of senescence
!$OMP THREADPRIVATE(senescence_delta)


  !
  ! stomate_windthrow.f90 / stomtate_pest.f90
  !

  ! 0. Constants

  REAL(r_std), SAVE :: one_third = 0.333              !! This value is used on multiple occasions in
                                                      !! stomate_windthrow.f90
                                                      !!(unitless)
!$OMP THREADPRIVATE(one_third)
  REAL(r_std), SAVE :: dbh_height_standard = 1.3      !! The height where the diameter of the tree stem is
                                                      !! measured by default.
                                                      !@tex $(m)$ @endtex
!$OMP THREADPRIVATE(dbh_height_standard)
  REAL(r_std), SAVE :: dbh_height_stump = zero        !! The height where the diameter of the tree stem is
                                                      !! measured if the middle of the canopy is below 1.3 m.
                                                      !! @tex $(m)$ @endtex
!$OMP THREADPRIVATE(dbh_height_stump)
  REAL(r_std), SAVE :: snow_density = 150.0           !! Density of snow (kg/m3). It should be considered
                                                      !! to calculate this value for simulations during
                                                      !  future development.
!$OMP THREADPRIVATE(snow_density)
  REAL(r_std), SAVE :: clear_cut_max = 20000.0        !! The maximum contiguous area allowed to be clearfelled
                                                      !! @tex $(m^{2})$ @endtex
!$OMP THREADPRIVATE(clear_cut_max)
  REAL(r_std), SAVE :: c_surface = 0.003              !! Surface Drag Coefficient (Raupach 1994) (unitless)
!$OMP THREADPRIVATE(c_surface)
  REAL(r_std), SAVE :: c_drag = 0.3                   !! Element Drag Coefficient (Raupach 1994) (unitless)
!$OMP THREADPRIVATE(c_drag)
  REAL(r_std), SAVE :: c_displacement = 7.5           !! Used by Raupach to calculate the zero-plane displacement (Raupach 1994) (unitless)
!$OMP THREADPRIVATE(c_displacement)
  REAL(r_std), SAVE :: c_roughness = 2.0              !! Used by Raupach to calculate the surface roughness length (Raupach 1994) (unitless)
!$OMP THREADPRIVATE(c_roughness)
  REAL(r_std), SAVE :: air_density = 1.2226           !! The value of air density (kg*m-3). If needed, this can be derived dynamically from
                                                      !! other modules of ORCHIDEE, but considering the range of values it can hold, it is probably
                                                      !! not worth additional calculations for being used in WINDTHROW.
!$OMP THREADPRIVATE(air_density)
  REAL(r_std), SAVE :: f_crown_weight = 1.136         !! This factor represents the weight of the overhanging crown when the tree stem is bent.
                                                      !! The origin of 1.136 is described in the supplementary material of Hale et al. 2015.
!$OMP THREADPRIVATE(f_crown_weight)

  INTEGER(i_std), SAVE :: legacy_years_wind = 5       !! The years used to calculate the maximum/mean/sum related to wind speed
!$OMP THREADPRIVATE(legacy_years_wind)
  INTEGER(i_std), SAVE :: wind_days = 3               !! The number of timesteps to caluclate mean wind speed to determine calcualte wind damage or not
!$OMP THREADPRIVATE(wind_days)
  REAL(r_std), SAVE :: wind_sum_threshold = 200       !! The threshold for the sum of wind_days wind speed for all time step to calculate wind damage
!$OMP THREADPRIVATE(wind_sum_threshold)
  REAL(r_std), SAVE :: wind_ratio_threshold = 3       !! The threshold for the ratio of wind to longterm wind for each time step to sum ratios
!$OMP THREADPRIVATE(wind_ratio_threshold)
  REAL(r_std), SAVE :: gap_threshold = 0.3            !! The threshold for the basal area loss to be considered to create the gap with edges
!$OMP THREADPRIVATE(gap_threshold)
  INTEGER(i_std), SAVE :: legacy_years = 3            !! The years used to calculate the total damage from disturbances with in legacy_years and the default is 3 years.
!$OMP THREADPRIVATE(legacy_years)
  INTEGER(i_std), SAVE :: legacy_years_wood = 1       !! The years used to calculate the wood leftover legacy.
!$OMP THREADPRIVATE(legacy_years_wood)
  INTEGER(i_std), SAVE :: beetle_legacy = 6           !! The years used to calculate wood leftover.
!$OMP THREADPRIVATE(beetle_legacy)

! 1. Scalar
  REAL(r_std), SAVE :: wind_speed_storm_thr=20.0      !! the wind speed threshold above which is_storm flag is set to TRUE 
!$OMP THREADPRIVATE(wind_speed_storm_thr)
  REAL(r_std), SAVE :: daily_max_tune=0.1155          !! This is a linear tunning factor to adjust the calculated daily maximum wind speed from forcing dataset.
!$OMP THREADPRIVATE(daily_max_tune)
  REAL(r_std), SAVE :: elevate_wind=10.0              !! Height conversion of critical wind speeds above zero-plane displacement
!$OMP THREADPRIVATE(elevate_wind)
  REAL(r_std), SAVE :: spacing_adj=2.75
!$OMP THREADPRIVATE(spacing_adj)
  REAL(r_std), SAVE :: h_threshold_wind = 0.3         !! Ratio of height threshold to pipe_tune2 for truncating wind damage
!$OMP THREADPRIVATE(h_threshold_wind)
  INTEGER(i_std), SAVE :: nb_days_storm=5             !! the number of days at which the max wind speed is less than wind_speed_storm_thr
!$OMP THREADPRIVATE(nb_days_storm)
  INTEGER(i_std), SAVE :: nb_years_bgi=3              !! the number of years needed to average the bark beetle generation index 
!$OMP THREADPRIVATE(nb_years_bgi)
  REAL(r_std), SAVE :: replace_threshold=0.3          !! the threshold for the ratio of number of individuals revmoved from the naturl 
                                                      !! disturbances to replace stand to the youngest age class
!$OMP THREADPRIVATE(replace_threshold)


  ! stomate_season.f90 
  !

  ! 1. Scalar

  REAL(r_std), SAVE :: gppfrac_dormance = 0.2  !! report maximal GPP/GGP_max for dormance (0-1, unitless)
!$OMP THREADPRIVATE(gppfrac_dormance)
  REAL(r_std), SAVE :: tau_climatology = 20.   !! tau for "climatologic variables (years)
!$OMP THREADPRIVATE(tau_climatology)
  REAL(r_std), SAVE :: hvc1 = 0.019            !! parameters for herbivore activity (unitless)
!$OMP THREADPRIVATE(hvc1)
  REAL(r_std), SAVE :: hvc2 = 1.38             !! parameters for herbivore activity (unitless)
!$OMP THREADPRIVATE(hvc2)
  REAL(r_std), SAVE :: leaf_frac_hvc = 0.33    !! leaf fraction (0-1, unitless)
!$OMP THREADPRIVATE(leaf_frac_hvc)
  REAL(r_std), SAVE :: tlong_ref_max = 303.1   !! maximum reference long term temperature (K)
!$OMP THREADPRIVATE(tlong_ref_max)
  REAL(r_std), SAVE :: tlong_ref_min = 253.1   !! minimum reference long term temperature (K)
!$OMP THREADPRIVATE(tlong_ref_min)

  ! 3. Coefficients of equations

  REAL(r_std), SAVE :: ncd_max_year = 3.
!$OMP THREADPRIVATE(ncd_max_year)
  REAL(r_std), SAVE :: ngd_threshold = -5.
!$OMP THREADPRIVATE(ngd_threshold)
  REAL(r_std), SAVE :: green_age_ever = 2.
!$OMP THREADPRIVATE(green_age_ever)
  REAL(r_std), SAVE :: green_age_dec = 0.5
!$OMP THREADPRIVATE(green_age_dec)
  
  REAL(r_std), SAVE :: ngd_min_dormance = 55.
!$OMP THREADPRIVATE(ngd_min_dormance)

  !
  ! sapiens_forestry.f90
  !

  INTEGER(i_std), SAVE      :: ncirc = 1                  !! Number of circumference classes used to calculate C allocation. This creates within stand heterogeneity (i.e., structured canopies).
!$OMP THREADPRIVATE(ncirc)
  INTEGER(i_std), SAVE      :: nagec = 1                  !! Number of age classes. This creates landscape-level heterogeneity (i.e. the effects of deforestation and harvesting are longer seen in the energy budget)
!$OMP THREADPRIVATE(nagec)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: age_class_bound !! Boundaries of the age classes defined as threshold diameters (m)
!$OMP THREADPRIVATE(age_class_bound)
  LOGICAL, SAVE             :: forced_clear_cut = .FALSE. !! flag at wich the stand is force to clear-cut
!$OMP THREADPRIVATE(forced_clear_cut)
  INTEGER(i_std), SAVE      :: noutdiaclass = 16          !! Number of diameter classes in the output files. ncirc is the number of circumference classes
                                                          !! ORCHIDEE. For tree demography it could be more convenient to get the results in diameter
                                                          !! classes with fixed boundaries.
!$OMP THREADPRIVATE(noutdiaclass)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:) :: out_dia_class !! Boundaries of the diameter classes (no longer circumference as in
                                                          !! ORCHIEE!) (m). Given how the code is setup, this list should start with zero and end with
                                                          !! a very high and unlikely diameter. The length is thus noutdiaclass+1.
!$OMP THREADPRIVATE(out_dia_class) 

 !+++++++++ DEBUG +++++++
  LOGICAL, SAVE             :: ld_fake_height=.FALSE.     !! a flag to turn on the statements
!$OMP THREADPRIVATE(ld_fake_height)
  INTEGER(i_std), SAVE      :: test_pft = 2               !! Number of PFT for which detailed output 
!$OMP THREADPRIVATE(test_pft)
  INTEGER(i_std), SAVE      :: test_grid = 1              !! Number of the grid square for which detailed output 
                                                          !! If the default value is not one, this can cause crashes in debugging for small regions.
!$OMP THREADPRIVATE(test_grid)
  LOGICAL, SAVE             :: lnvgrasspatch=.FALSE.    !! a patch for grasslands
!$OMP THREADPRIVATE(lnvgrasspatch)
!++++++++++++++++++++++++
  
  LOGICAL,PARAMETER                       :: jr_nextstep = .FALSE.!! set this to TRUE to activate the 
 
  INTEGER(i_std), SAVE      :: ndia_harvest             !! The number of diameter classes used for 
                                                        !! the wood harvest pools.
!$OMP THREADPRIVATE(ndia_harvest)

  REAL(r_std), SAVE         :: max_harvest_dia = 0.7    !! The largest diameter for the harvest pools to
                                                        !! keep track of harvested wood from forests.
!$OMP THREADPRIVATE(max_harvest_dia)

  INTEGER(i_std), SAVE      :: n_pai = 1000             !! The number of years used for the cumulative
                                                        !! averages of the periodic annual increment. By
                                                        !! setting this value to 1000, this functionality is not used.
                                                        !! This way of management is extremly rational. Owners are
                                                        !! to follow yield tables than to monitor stand growth every 5 years.                     
!$OMP THREADPRIVATE(n_pai)

  INTEGER(i_std), SAVE      :: ntrees_profit            !! The number of trees over which the average
                                                        !! height is calculated to determine if the
                                                        !! stand will be profitable to thin.
!$OMP THREADPRIVATE(ntrees_profit)

INTEGER, SAVE             :: species_change_force      !! This is the PFT number which is replanted after a
  	                                               !! clearcut, if such a thing is being done.
                                                       !! To be used with lchange_species = .TRUE. and
	                                               !! lread_species_change_map = .FALSE. The
	                                               !! forced value is mainly useful for debugging
!$OMP THREADPRIVATE(species_change_force)

INTEGER, SAVE             :: fm_change_force           !! This is the FM strategy which is used for the replant
	                                               !! after a clearcut, if such a thing is being done.
 	                                               !! To be used with lchange_species = .TRUE. and
       	                                               !! lread_desired_fm_map = .FALSE. The forced value is
	                                               !! mainly useful for debugging
!$OMP THREADPRIVATE(fm_change_force)

 REAL(r_std), SAVE        :: min_water_stress = 0.1    !! Minimal value for wstress_fac (unitless, 0-1)
!$OMP THREADPRIVATE(min_water_stress)

REAL(r_std), SAVE         :: max_delta_KF = 0.1        !! Maximum change in KF from one time step to another (m)
                                                       !! This is a bit arbitrary. 
!$OMP THREADPRIVATE(max_delta_KF)

REAL(r_std), SAVE         :: maint_from_gpp = 0.8      !! Some carbon needs to remain to support the growth, hence, 
                                                       !! respiration will be limited. In this case resp_maint 
                                                       !! (gC m-2 dt-1) should not be more than 80% (::maint_from_gpp) 
                                                       !! of the GPP (gC m-2 s-1)
!$OMP THREADPRIVATE(maint_from_gpp)
 
  REAL(r_std), PARAMETER :: m2_to_ha = 10000.          !! Conversion from m2 to hectares
  REAL(r_std), PARAMETER :: ha_to_m2 = 0.0001          !! Conversion from hectares (forestry) to m2 (rest of the code)
  REAL(r_std), PARAMETER :: m_to_cm = 100.             !! Conversion from m to cm
  REAL(r_std), PARAMETER :: cm_to_m = 0.01             !! Conversion from cm to m
  REAL(r_std), PARAMETER :: peta_to_unit = 1.0E15      !! Convert Peta to unit
  REAL(r_std), PARAMETER :: tera_to_unit = 1.0E12      !! Convert Tera to unit 
  REAL(r_std), PARAMETER :: giga_to_unit = 1.0E09      !! Convert Giga to unit
  REAL(r_std), PARAMETER :: mega_to_unit = 1.0E06      !! Convert Mega to unit
  REAL(r_std), PARAMETER :: kilo_to_unit = 1.0E03      !! Convert Kilo to unit
  REAL(r_std), PARAMETER :: centi_to_unit = 1.0E02     !! Convert centi to unit
  REAL(r_std), PARAMETER :: milli_to_unit = 1.0E-03    !! Convert milli to unit
  REAL(r_std), PARAMETER :: carbon_to_kilo = 2.0E-03   !! Convert g carbon to kilo biomass

  !
  ! Debugging
  !
  INTEGER(i_std), SAVE :: err_act = 1                  !! There are three levels of error checking 
                                                       !! see constantes.f90 for more details
!$OMP THREADPRIVATE(err_act)
  INTEGER(i_std), SAVE :: plev = 0                     !! print level of the subroutine ipslerr_p 
                                                       !! (1:note, 2: warn and 3:stop) 
!$OMP THREADPRIVATE(plev) 
  REAL(r_std), SAVE    :: sync_threshold = 0.0001      !! The threshold above which a warning is generated when the
!$OMP THREADPRIVATE(sync_threshold) 

  ! Soil carbon related variables
  LOGICAL, SAVE                   :: soilc_isspinup = .FALSE.            !! Activate soil carbon spinup
!$OMP THREADPRIVATE(soilc_isspinup)
  REAL(r_std), PARAMETER          :: O2_init_conc = 298.3                !! gO2/m**3 mean for Cherskii
  REAL(r_std), PARAMETER          :: CH4_init_conc = 0.001267            !! gCH4/m**3 mean for Cherskii
  REAL(r_std), PARAMETER          :: z_root_max = 0.5                    !! Depth at which litter carbon input decays e-fold (root depth); 0.5 for compar w/WH
  REAL(r_std), PARAMETER          :: diffO2_air = 1.596E-5               !! oxygen diffusivity in air (m**2/s)
  REAL(r_std), PARAMETER          :: diffO2_w = 1.596E-9                 !! oxygen diffusivity in water (m**2/s)
  REAL(r_std), PARAMETER          :: O2_surf = 0.209                     !! oxygen concentration in surface air (molar fraction)
  REAL(r_std), PARAMETER          :: diffCH4_air = 1.702E-5              !! methane diffusivity in air (m**2/s)
  REAL(r_std), PARAMETER          :: diffCH4_w = 2.0E-9                  !! methane diffusivity in water (m**2/s)
  REAL(r_std), PARAMETER          :: CH4_surf = 1700.E-9                 !! methane concentration in surface air (molar fraction)
  REAL(r_std), SAVE               :: tetasat =  .5                       !! volumetric water content at saturation (porosity)
!$OMP THREADPRIVATE(tetasat)
  REAL(r_std), SAVE               :: tetamoss =  0.92                    !! porosity of moss
!$OMP THREADPRIVATE(tetamoss)
  REAL(r_std), SAVE               :: rho_moss                            !! density of moss
!$OMP THREADPRIVATE(rho_moss)
  REAL(r_std), PARAMETER          :: zmoss = 0.2                         !! thickness of moss layer (in permafrost regions,m) 0. ! 0.001 DKtest for compar w/WH
  REAL(r_std), PARAMETER          :: h_snowmoss = 0.2                    !! snow height above which we consider the moss layer to be to compressed to be effective (m)
  REAL(r_std), PARAMETER          :: BunsenO2 = 0.038                    !! Bunsen coefficient for O2 (10C, 1bar)
  REAL(r_std), PARAMETER          :: BunsenCH4 = 0.043                   !! Bunsen coefficient for CH4 (10C, Wiesenburg et Guinasso, Jr., 1979)
  REAL(r_std), PARAMETER          :: ebuthr = 0.9                        !! Soil humidity threshold for ebullition
  REAL(r_std), PARAMETER          :: wCH4 = 16.                          !! molar weight of CH4 (g/mol)
  REAL(r_std), PARAMETER          :: wO2 = 32.                           !! molar weight of O2 (g/mol)
  REAL(r_std), PARAMETER          :: wC = 12.                            !! molar weight of C (g/mol)
  REAL(r_std), PARAMETER          :: avm = .01                           !! minimum air volume (m**3 air/m**3 soil)
  REAL(r_std), PARAMETER          :: hmin_tcalc = .001                   !! minimum total snow layer thickness below which we ignore diffusion across the snow layer

  

END MODULE constantes_var
