ORCHIDEE trunk

Description of usage of the main functionalities added in the trunk since 2018.

Good to know

The trunk since tag 4.0 is a combination of three code bases: the previous trunk (branched as ORCHIDEE 3) as it was in 2020, the ORCHIDEE-CN branch in 2016 and the key functionality of ORCHIDEE-CAN.

• The available functionality is described in more detail below
• Missing functionalities compared to version ORCHIDEE 2 (trunk until March 2018): DGVM, IPCC wood harvest, prescribed LAI, FORCESOIL, and fire disturbances
• Missing functionalities compared to ORCHIDEE-CN/ORCHIDEE 3 (trunk between March 2018- March 2020): none
• Missing functionalities compared to ORCHIDEE-CAN(merged in the trunk March 2020): multi-layer energy budget
• Coupling: ORCHIDEE trunk is coupled to LMDZ

Fair use policy

Before downloading and using the trunk, please, read and adhere to the ​''FAIR USE POLICY''.

Warranty

ORCHIDEE (milestone 4) comes with a 12 month warranty (starting on March 10th 2020)! This warranty is valid for the head of the trunk if no code changes were made and for the default configurations (analytical spinup, FG1, FG2, FG4 and FG5). The warranty applies to model crashes, ipsl errors and obvious bugs (e.g., spikes, sudden changes on January 1st, crazy values, …). Following the establishment of the problem, you need to extract a 1 pixel or 3x3 pixel test case that reproduces the crash on Obelix. A script to extract the restart files for this kind of restricted test cases can be downloaded from svn ( ​download here). Once the TRUNK FORCE has received the test case, the first response time is 5 working days (outside the holiday period). If your problem falls outside the warranty we can still give assistance and advice.

How do I run the trunk?

The instructions to download and extract ORCHIDEE are found here: wiki:Documentation/UserGuide.

Functionalities (alphabetical order)

Age classes

Describes r6614. Age classes were introduced to better handle heterogeneity at the landscape level. The feature allows us to distinguish between different successional stages of the same PFT (e.g., a newly grown forest vs. a mature forest). Age classes are independent of the number of diameter classes. Using age classes adds a lot of details to both the biophysics and the biogeochemistry following natural disturbances, forest management and land cover change. If half of a grassland is afforested with a PFT that already exists in the pixel, previous versions of ORCHIDEE will combine this newly forest land and the existing forest in a single PFT. This will result in, for example, a low albedo, a high roughness, and other properties. When age classes are used, the newly afforested and the existing forest will end up in separate PFTs. One will have a high albedo, the other a low albedo, and other properties may differ significantly as well. In CAN with age classes, PFTs are only merged if the youngest age class for a PFT already has biomass.

Age classes are defined as separate PFTs. Different age classes of the same PFT could therefore be, in principle, run with different parameters. This option has not been tested yet because it is expected to result in discontinuities when the biomass is moved from one age class to another. The number of age classes is fixed for the whole simulation, but for each PFT it can be decided whether age classes are used or not. In other words, if the user chooses four age classes for the simulation, each PFT can have either 1 or 4 age classes. This adds a lot of flexibility to the model. ORCHIDEE trunk, for example, has been run with 64 PFTs, using age classes for European forest and using no age classes for all forests outside the domain of interest. Setting-up a simulation with age classes will require some thinking when creating the orchidee_pft.def. A Python-script was written to create this kind of run.def and is stored in config/ORCHIDEE_OL/MAKE_RUN_DEF. Increasing the number of PFTs has important consequences for the speed of the model and the memory use, although ORCHIDEE trunk does make extensive use of "CYCLE" statements to avoid calculations where no biomass is present. Because a single run can contain PFTs with and PFTs without age classes, processing of the simulation output needs to account for the relationship between PFTs of the same species but a different age class.

It does not make sense to use age classes for runs ignoring land cover change and forest management. Age classes can be used for site level simulations that involve land cover change. Only use age classes if you really need them (e.g., land cover change, forest management, other disturbances), as not using age classes will make post processing of the simulation results considerably easier.

The number of age classes is defined by the parameter NAGEC. Setting this parameter to 1 is a good start unless you have a special interest in using age classes. When NAGEC is set to more than 1, PFT_TO_MTC', AGEC_GROUP and PFT_NAME will all need to be carefully defined. See the different orchidee_pft.defs that come with the standard configurations in config/ORCHIDEE_OL for functional examples. See below for some principles:

• NAGEC = 4
• Assume we want to use four age classes for all forests. We will end-up with 37 PFTs: one each for bare soil, three C3 grasses, C4 grass, C3 crop and C4 crop and 4 times 8 for the 8 forest PFTs. Thus NVM = 39
• Because we still use the 13 default MTCs we can use the default maps. Let the model know how many MTCs it should find on the maps: NVMAP=15
• If you want to use IMPOSE_VEG=y then vegetation should only be placed in the youngest age class. ORCHIDEE will update the vegetation fractions during the simulations

  SECHIBA_VEG_01=0.0769230769231
  SECHIBA_VEG_02=0.0769230769231
  SECHIBA_VEG_03=0.0
  SECHIBA_VEG_04=0.0
  SECHIBA_VEG_05=0.0
  ...
  SECHIBA_VEGMAX_01=0.0769230769231
  SECHIBA_VEGMAX_02=0.0769230769231
  SECHIBA_VEGMAX_03=0.0
  SECHIBA_VEGMAX_04=0.0
  SECHIBA_VEGMAX_05=0.0
  ...
  

• Link PFTs to MTCs

  PFT_TO_MTC_01=1
  PFT_TO_MTC_02=2
  PFT_TO_MTC_03=2
  PFT_TO_MTC_04=2
  PFT_TO_MTC_05=2
  ...
  

• Tell ORCHIDEE which PFTs have a successional relationship

  AGEC_GROUP_01=1
  AGEC_GROUP_02=2
  AGEC_GROUP_03=2
  AGEC_GROUP_04=2
  AGEC_GROUP_05=2
  ...
  

• Give all PFTs a name for clarity

  PFT_NAME__01=Soilbare
  PFT_NAME__02=Broadleavedevergreentropical_age01
  PFT_NAME__03=Broadleavedevergreentropical_age02
  PFT_NAME__04=Broadleavedevergreentropical_age03
  PFT_NAME__05=Broadleavedevergreentropical_age04
  ...
  

Albedo

Describes r6614. ORCHIDEE trunk makes use of a two stream radiative transfer scheme through the canopy, extended to multiple canopy levels (​https://doi.org/10.5194/gmd-2016-280). The scheme is based on Pinty et al 2006. This approach accounts not only for the leaf mass but also for the vertical and horizontal distribution of the leaf mass (=canopy structure), calculating an effective LAI based on the solar angle. Light from collimated (black sky) and diffuse (white sky) sources are used, and both are weighted equally as information about this partitioning is not yet available in forcing data. In ORCHIDEE trunk the same scheme is used to simulate the reflected, transmitted and absorbed light, of which the absorbed light as a function of canopy level is passed to the photosynthesis routines and used in place of the exponential LAI layering found in older versions of the TRUNK (see the section Photosynthesis). This implies that albedo and photosynthesis are now fully consistent as well as the light reaching the forest floor (the latter is used in for example recruitment). ORCHIDEE trunk cannot revert to previous approaches for calculating albedo. The radiative transfer through the canopy is controlled by 3 parameters for each wavelength/band: single leaf scattering leaf_ssa_xxx, forward scattering leaf_psd_xxx and background reflectance bgrd_ref_xxx. At present, both the visible (VIS) and near-infrared (NIR) spectra have been parameterized. Parameterization is based on running an inverse radiation scheme on the MODIS albedo product while accounting for the different land cover types. The inverted parameters are provided by the JRC as the JRC TIP product. Seasonal variation in the background albedo was observed but small and therefore not accounted for.

When snow is present in a pixel, all snow is assumed to reach the ground and the background albedo and the snow albedo (calculated as a function of snow age) are weighted according to their cover fractions (see Background albedo).

Albedo (background)

Describes r7253. If covered by snow, the background albedo is calculated in the VIS and NIR wavelengths by the snow module and accounts for snow age and snow density. If not covered by snow, the background albedo is not simulated but prescribed by the parameters bgd_reflectance(:,ivis) and bgd_reflectance(:,inir). If the background is partly covered by snow, the snow albedo and the background albedo are merged, which allows snow to settle under the canopy. In deciduous forest, grasslands and croplands, the background albedo is known to be strongly affected by the phenology and senescence of the understory vegetation.

ORCHIDEE trunk has three options to calculate the background albedo for the bare soil (PFT1):

• It can read the soil color (set alb_bg_model = y and alb_bg_modis = n) and account for soil moisture
• It can read the soil color (set alb_bg_model = n and alb_bg_modis = n) and ignore soil moisture
• It can read the MODIS background albedo (alb_bg_modis = y , value of alb_bg_model doesn't matter). This option also ignores the soil moisture.

ORCHIDEE trunk has two options to calculate the background albedo for vegetated pixels:

• The background albedo is prescribed per PFT but is constant throughout the year. This option has been evaluated over Europe. Set alb_bg_modis = n.
• The background albedo is constant across PFTs. This option reads background maps. Given that those maps are based on the JRC TIP product, they should be compatible with the new albedo scheme. This option, however, has not been evaluated yet. Set alb_bg_modis = y.

When alb_bg_modis = y this is the option used to calculate the background albedo for PFT1 and vegetated PFTs.

Albedo (snow)

Describes r6614. The snow albedo could be either prescribed (the value of fixed_snow_albedo should then be set in the run.def and will be read in src_parameters/constantes_var.f90) or calculated following Chalita and Treut (1994) ok_snow_albedo_clm3 = n or calculated following CLM3 ok_snow_albedo_clm3 = y (recommended). The difference between the latter two methods has not been tested yet. The CLM method was implemented in ORCHIDEE-CN-CAN and has been tested in combination with the other changes to the albedo calculations. The Chalita and Treut method was added in parallel to the ORCHIDEE trunk 3.0. When merging both versions, we ended up with two option.

Allocation

Describes r7066. ORCHIDEE trunk uses the allometric allocation as developed in O-CN. In ORCHIDEE trunk the approach was adjusted to work for more than one diameter class. Since it was developed this allocation has been used in ORCHIDEE-CN and ORCHIDEE-CNP. In those branches only a single diameter class was used. Except for the way the reserves and labile pools are calculated (incl. the pseudo sugar loading), the allocation scheme remained rather similar between the aforementioned versions. The model is, however, very sensitive to the way the reserves and labile pools are calculated. The allocation makes use of a labile pool for which the activity is calculated based on the temperature. This sensitivity is important at the start and the end of the growing seasons when temperatures may be low. As such the model addresses the sink/source discussion initiated by Körner. Whereas this approach resulted in an acceptable interannual variability in for example NPP in ORCHIDEE-CAN, adding N seems to have dampen the interannual variability a lot/too much. This dampening was observed in ORCHIDEE-CN and ORCHIDEE-CN-CAN. IN ORCHIDEE-CNP the temperature relationship was removed (hence NPP and GPP are strictly coupled) because the interannual variability became unrealistic.

One of the big challenges for all deciduous PFTs is to grow first (in spring), restore the reserves towards the end of the growing season (late summer) and still have enough reserves to grow leaves the next year. To do so a phenological stage called pre-senescence was introduced. This is justified by the observation that tree rings stop growing weeks before senescence starts. Basically allocation occurs when the plant_status is icanopy. When the ratio between Rm and GPP (weekly mean value) drops below a PFT-specific trhreshold, plant_status is changed to ipresenescence. This phenological phase is used to restore the labile and reserve pools. The relationship between allocation and plant status is detailed in the table below.

Table. Allocation is a function of the plant status.

Process	Plant_status	Snipets of the code and comments
Fresh GPP destination	ibudbreak & icanopy	GPP goes to labile pool
	ipresenescence	GPP goes to reserve pool. Divert fresh GPP to reserve pool so that it’s not readily available for allocation. But still be available after rebalance between reserve and labile pool.
	isenescent	GPP goes to reserve pool
	idormant	GPP is zero
Calculating bm_alloc_tot	ibudbreak & icanopy	(1) gtemp > 0; (2) IF use_reserve, shift 5% of reserve to labile pool; and (3) bm_alloc_tot > 0. bm_alloc_tot is taken from labile carbon. bm_alloc_tot(ipts,j) = gtemp(ipts,j)*tmp_bm(ipts,j,ilabile,icarbon)
	ipresenescence	(1) gtemp > 0; (2) suppress the shift of 5% of reserve to labile pool; and (3) bm_alloc_tot > 0. bm_alloc_tot is taken from labile carbon. bm_alloc_tot(ipts,j) = gtemp(ipts,j)*tmp_bm(ipts,j,ilabile,icarbon)
	isenescent	gtemp = 0 => bm_alloc_tot = 0. bm_alloc_tot is taken from labile carbon. bm_alloc_tot(ipts,j) = gtemp(ipts,j)*tmp_bm(ipts,j,ilabile,icarbon)
	idormant	gtemp = 0 => bm_alloc_tot = 0. bm_alloc_tot is taken from labile carbon. bm_alloc_tot(ipts,j) = gtemp(ipts,j)*tmp_bm(ipts,j,ilabile,icarbon)
Maintenance Respiration (Rm)	Irrespective of plant_status	Rm deficit will be paid by reserve if Rm>bm_alloc_tot, otherwise limit Rm to “bm_alloc_tot + reserve”. Rm could deplete all reserves or reduce it to a too low value, giving trouble to leaf onset in next growing season.
	ibudbreak & icanopy	Rm deficit will be paid by reserve if Rm>bm_alloc_tot, otherwise limit Rm to “bm_alloc_tot + reserve”.
	ipresenescence	Rm deficit will be paid by reserve if Rm>bm_alloc_tot, otherwise limit Rm to “bm_alloc_tot + reserve”. Stop tissue allocation but do not truncate Rm. If Rm>bm_alloc_tot: it will be paid by labile; Otherwise, if Rm<bm_alloc_tot, the extra bm_alloc_tot will not go to tissue but go back to reserve. Divert the gpp not consumed by Rm to reserve pool, in order to increase pool.
	isenescent	Rm deficit will be paid by reserve if Rm>bm_alloc_tot, otherwise limit Rm to “bm_alloc_tot + reserve”. Rm = 0. Suppress Rm to not to lose the accumulated reserve pool.
	idormant & ibudavail	Rm deficit will be paid by reserve if Rm>bm_alloc_tot, otherwise limit Rm to “bm_alloc_tot + reserve”. Rm = 0. Suppress Rm to not to lose the accumulated reserve pool.
Rebalance between reserve and labile pool	Independent of plant_status	Allow exchange between reserve and labile pool to reach their optimal share

There are no options to revert to the allocation based on resource limitation (Friedlingstein et al. 1999). All references and parameters for allocation based on resource limitation have been removed from the code (those that were overlooked can be removed). Allometric allocation makes use of the following PFT-specific parameters: sla, longevity_root, longevity_leaf, longevity_sap, pipe_density, tree_ff, pipe_tune_x, k_latosa_max, and k_latosa_min. In addition to this set of parameters that mainly describe the allometric relationships and the longevity of the different tissues, the calculation of the allocation coefficients makes use PFT-specific tissue conductivities, i.e., k_sap, k_belowground, and k_leaf (see also plant water stress). Details on the parameters can be found in the SI of Naudts et al 2015 in GMD or in src_parameters/constantes_mtc.f90.

Previously there was a functional link between C and N-allocation and the hydraulic architecture of a plant because both approaches used the same parameter k_root. In ORCHIDEE-CAN k_root described the conductivity of the fine roots and the soil. In ORCHIDEE trunk this joined conductivity has been split in a fine root conductivity and a soil to root conductivity. Allocation should make use of both conductivities but soil to root conductivity cannot be easily calculated when needed in the allocation. This is subject to future developments. Accounting for the soil to root conductivity in the allocation would imply an adaptation of plant growth to its environment.

Crown area is calculated from tree diameter by making use of pipe_tune1 and pipe_tune_exp_coeff. Since r6933, pipe_tune1 and pipe_tune_exp_coeff can be calculated as a function of respectively pipe_tune2 and pipe_tune3. This has became the default setting. For that reason pipe_tune1 and pipe_tune_exp_coeff are undefined but the parameters crown_to_height and crown_to_height_exp need to be set. If one wants the crown diameter to be 66% of the tree height, crown_to_height needs to be set to 0.66 and crown_to_height_exp to 2. As long as crown_to_height_exp equals to 2, crown_to_height should range between 0 and 0.99. If a different values for crown_to_height_exp it should be checked that the crown diameter never exceeds the tree height as this will make the model crash with a clear error message.

Anthropogenic species change

Describes r6614. Following a disturbance (which could be a clear cut), tree species changes and forest management change can be prescribed or read from a map in ORCHIDEE trunk. Set ok_change_species = y, read_species_change_map = y, and read_desired_fm_map = y and specify the paths of those maps in the COMP/stomate.card. A example of such a configuration can be found in config/ORCHIDEE_OL/OOL_SEC_STO_FG5. This functionality replaces the DGVM in areas where humans rather than nature govern species distribution, for example, Europe. Note that there are some constraints on the possible species changes. If the forest is unmanaged (fm=1), the code assumes that nature will determine the species rather than humans. Anthropogenic species change has not been developed to work together with land cover change. For the moment it is one or the other. When testing this functionality read_species_change_map and/or read_desired_fm_map could be set to n. The new forest management strategy can then be simply prescribed by setting the parameter fm_change_force to one of the four fm strategies. Likewise the new species can be prescribed by setting the parameter species_change_force to the desired PFT number.

Bare soil

Describes r6783. The flag ok_bare_soil_new controls how the bare soil is perceived and calculated. If set to FALSE the total bare soil is still calculated veget_max_1 + sum(veget_max_i - veget_i) with i from 2 to the number of PFTs. When a deciduous PFT sheds its leaves, the gaps in the forest will thus contribute to the bare soil fraction in the grid. Although this approach was introduced a long time ago to get acceptable evaporation estimates from forest, the approach also resulted in using the albedo of PFT1 deserts as the background albedo of forest gaps. In ORCHIDEE v2.1 the background albedo has been reparameterized and this issue may have been largely resolved now if alb_bg_modis = y (​more details). From many points of view a dynamic bare soil fraction is strange, e.g., bare soil has its own water column so when moving the forest gaps from to forest to PFT1 the soil water content, soil carbon content, litter layer, etc all changes temporary. If ok_bare_soil_new is set to TRUE, canopy gaps no longer contribute to the bare soil. The new albedo scheme (see Albedo and Background albedo) considers a specific background albedo for each PFT and calculates the albedo of the PFT including the canopy gaps but the calculation of bare soil evaporation underneath a canopy would be problematic. For that reason ok_bare_soil_new is recommended only to be used with the multi-layer energy budget (when run with more than 1 layer). The multi-layer energy budget accounts for within canopy turbulence and can therefore deal with evaporation from beneath a canopy. At present the settings for ok_bare_soil_new are included in the energy_control flag (​more details). The default settings combine the new albedo scheme with the single layer energy budget (enerbil) and ok_bare_soil_new = n.

Bare soil evaporation

The calculation of vbeat4 in r7278 was (numerically) confirmed to be identical to the calculation of vbeta4 in ORC2 for a pixel with 100% bare soil coverage. At sun set the bare soil evaporation starts to deviate between both models because differences in the albedo calculation. When prescribing the ORC2 albedo in r7278 the exact same bare soil evapotranspiration is obtained for the whole 10 years of the test run. When running global simulations there is a 40% higher bare soil evaporation in r7278 compared to ORC2. Approximative factorial simulations with fixed albedo and fixed veget_max suggest that the difference in veget for vegetated PFTs contributes most to the difference in bare soil evapotranspiration. Given that veget is calculated through Pgap, the nature of Pgap seems to be the cause.

ORC2 calculates veget with Lambert Beer and vbeta4 was most likely tuned to obtain reasonable bare soil evaporation (is it reasonable to start with? CAn it be tuned?). ORC4 uses Pgap which is from an ecological point of view superior to Lambert Beer because it accounts for clumping and could account for stand structure. When combing Pgap with ORC2 vbeta4, bare soil evapotranspiration becomes too high. In the long term, Pgap should be kept and the whole definition of "bare soil"should be revised. This may happen as part of the multi-layer energy budget calculations.

Bark beetles

Describes r6614. The bark beetle module was developed such that it interacts with the windthrow module. If you want to activate it use the flag OK_PEST=y in the run.def. To specify which PFT may be affected by bark beetles use BEETLE_PFT_xxx= TRUE where XXX is the pft you interested in. Note that for the moment bark beetles in ORCHIDEE is parameterized only to work with Picea abies over Europe.

C13

Describes r6614. The concentration of C13 in the leaves can be calculated by setting ok_c13 to y in the run.def. This calculation follows Farquhar's approach and is currently limited to the fractionation in the leaves. Following changes to the recalculation of GPP under plant water stress the calculation of Ci is no longer accurate. This may have broken the functionality to calculate C13. Needs to be tested.

Comparing to forest inventory observations

Describes r7773. ORCHIDEE calculates the HEIGHT, DIAMETER, BA, and IND for the whole stand including saplings. Forest inventories use a diameter threshold below which trees are not sampled. The HEIGHT, DIAMETER, BA, and IND from a forest inventory cannot simply be compared to the ORCHIDEE simulations. This cut off diameter has been implemented in ORCHIDEE (DIA_THRESH_INV) and enables the model to calculate the simulated HEIGHT, DIAMETER, BA, and IND for the whole stand as well as these variables as if they would be measured in a forest inventory. These variables are written to the history file as HEIGHT, HEIGHT_DOM, HEIGHT_INV, DIAMETER, DIA_DOM, DIA_INV, BA, BA_INV, IND, and IND_INV.

Configurations

The model ORCHIDEE can run in offline mode. Different predefined experiment set ups can be found in folder ORCHIDEE_OL, see section Experiments further below.

The model can also be coupled to an atomspheric model and be run in coupled mode. For example in the coupled model IPSL-CM, ORCHIDEE is coupled to the atmospheric model LMDZ and the ocean model NEMO.

Consistency checks

Describes r6947. The code distinguishes between three options to check for mass and surface conservation. These options are controlled by the parameter err_act. Always use err_act = 3 when developing and testing the code. Note that in addition to checking for mass balance closure ORCHIDEE trunk will also check for the conservation of veget_max and frac_nobio. This is useful to make sure no surface area is lost when moving biomass from one PFT to another following natural disturbances, forest management, land cover changes and when using age classes. In some parts of the code, for example, modules that deal with disturbances, it is assumed that the tallest trees are stored in the last diameter class. When the difference in diameter between diameter classes becomes very small, this assumption could be violated. Therefore, the diameter classes are sorted to enforce the assumed order and where needed the order is checked.

• err_act = 1 is recommended when running global long-term simulations. Under this option, mass balance closure is checked for all biogeochemical processes but only at the highest level thus stomate.f90 and stomate_lpj.f90. Although the mass balance checks are not very expensive in terms of computer time, skipping the numerous lower level checks is expected to save some time. Under this option the total mass balance error is only written to the history file. No information is provided in which subroutine the problem occurred.
• err_act = 2 is recommended when developing and testing the model. Now the mass balance is explicitly checked in stomate.f90, stomate_lpj.f90 and all its subroutines. Under this option the mass balance error is written to the history file and if the mass balance is not closed, the warning message will indicate in which subroutine the problem likely originated.
• arr_act = 3 is recommended when having a problem with mass balance closure. The mass balance is explicitly checked in stomate.f90, stomate_lpj.f90 and all its subroutines. If a mass balance error occurs, the model is stopped. It is expected that within the time step of the subroutine a precision of 1e-8 is achieved. If not the model will crash. Additionally consistency between flux-based and pool-based NBP is checked. This check represent an integrated check which give confidence for the technical implementation of the C and N cycles. Comparing the flux and pool based nbp estimates is an instantaneous check that covers just a single day in the model (thus 48 time steps for the soil and litter calcuations which are called from stomate.f90). The cross-consistency check for the flux and pool based NBP requires mass balance closure and area preservation from all all subroutines being called from stomate.f90 including stomate_lpj.f90. This check should always be satisfied, if not there is a bug somewhere! The mass balance of individual subroutines teaches us about the calculations performed within that subroutine. The flux and pool based NBP consistency check teaches us of how these fluxes and pools should be accounted for, passing this check requires additional understanding of the flow of the code (i.e., bm_to_litter) and how the different time steps come together.

As an additional consistency check the time integral of the NBP (from 0 -> t) should equal the C and N stock at time step t. Comparing the time integral of NEP vs the all pools includes all the accumulated precision errors. If those errors would be random there would be no problem but they are more likely biases. Hence, when the fluxes and pools increase several order of magnitudes (e.g toward 1e4) and a small bias (1e-12) is integrated over 10,000 days (1e4) the divergence between the time integral of NEP and the pools will approach/exceed 1-8 and woukd stop the model. This is a fundamental problem related to machine precision and the way we do the calculations (think of the so called catastrophic rounding errors that may occur with floating numbers). For this reason this check is written to the history files but does not make the model crash. The user will have to decide whether (s)he is facing a precision error or a real bug. Fixing precision errors will require to improve the precision of the calculations in one or several subroutines.

The variable MBC_NBP3_c was added to enable monitoring whether the NBP send to LMDZ in intersurf.f90 is identical to the NBP calculated in stomate_lpj.f90. If ERR_ACT=3 is used, the model will crash if the NBP send in intersurf differs from the NBP calculated in stomate_lpj.f90. Irrespective of the value of ERR_ACT, the mismatch between both NBPs will be send to MBC_NBP3_c.

Croplands

Describes r7767. Makes use of sapiens_planting and sapiens_harvest. Planting and harvest is determined by temperature and soil moisture using the moigdd phenology type. At senescence, 50% of the aboveground biomass is placed into the harvest_pool (just like wood) and its decomposition is accounted for in the fluxes from the product use pool. The remaining 50% of aboveground biomass enters the litter pool, along with all of the roots.

Crops: ipresenescence

With the flag, ok_crops_ipresenescence, it is now possible for crops to enter presenescence. Entering presenescence enables growth to be directed towards reserves and labile pools and reducing LAI. Two new PFT parameters are introduced in the check for ipresenescence: pre_gdd and pre_harv. If the flag is set to false, crops move directly from icanopy to isenescent. This means that leaf growth continues until senescence, resulting in high LAI.

Crop planting: nitrogen from soil

At the time of crop planting, the required nitrogen can be taken from one of the soil nitrogen pools. To use this feature, set ok_crop_planting_n_soil to TRUE. Nitrogen is preferentially taken from the mineralised nitrogen pool (soil_n_min) and if there is not enough (there usually isn’t), the nitrogen is taken from the active, slow, or passive som pools, depending on availability. In the case that there is not enough nitrogen available in the soil_n_min or som pools, nitrogen is taken from the atmosphere (atm_to_bm). The instances where nitrogen is taken from the atmosphere at the time of planting are recorded using the N_ESTAB_ATM output variable. If ok_crop_planting_n_soil is set to false, nitrogen is taken from the atmosphere at the time of planting (atm_to_bm).

CWRR vs Choinel

Describes r6614. Following ORCHIDEE trunk 3, the Choinel code is no longer available in ORCHIDEE trunk. The hydrological schemes were not touched during the development of ORCHIDEE trunk. The hydraulic architecture in ORCHIDEE trunk relies on CWRR.

Diameter classes

Describes r6614. Diameter classes were introduced to better simulate the canopy structure, they are a tool to simulate heterogeneity within a single PFT. Given that the canopy is the interface between the land and the atmosphere this feature has effects well beyond forest management. Stand structure was observed to affect albedo, transpiration, photosynthesis, soil temperature, roughness length, and recruitment. Using diameter classes adds very little complexity to setting-up the simulations as well as to the output files. The complexity is mostly within the code.

The computational cost of using diameter classes is negligible and when a reasonable low number of diameter classes is used, the memory cost remains very small as the dimensions of only two variables are increased. The number of diameter classes is the same for all PFTs and is set by the parameter NCIRC. ORCHIDEE-CN, and ORCHIDEE-CNP are coded and used for NCIRC = 1. ORCHIDEE trunk was coded and tested for NCIRC = 3. Until further testing, three diameter classes are considered sufficient although for some specific application, i.e., simulating tree ring width NCIRC = 5 resulted in slightly better model performance. The branch ORCHIDEE-CN-CAN that preceded this trunk version has been run with one diameter class as well (this should be confirmed for ORCHIDEE trunk).

Given earlier choices in the ORCHIDEE model, we either need to define the boundaries of each diameter class or the diameter distribution. While developing the code, we considered the second approach the most flexible. To allow maximal flexibility, each diameter class needs to be defined separately by the variable CIRC_CLASS_DIST_0000X, where X is the number of the diameter class. The smallest number presents the smallest diameter class. From literature it is known that a truncated exponential distribution is a good first guess:

CIRC_CLASS_DIST_00001=9
CIRC_CLASS_DIST_00002=5
CIRC_CLASS_DIST_00003=1 

The above declaration implies that 9/15th of the trees will always be in the smallest diameter class, 5/15th will be in the medium class and 1 tree out of 15 will be in the largest diameter class. These ratios are kept throughout the simulations and the boundaries of the diameter classes are adjusted to respect this constraint. Consequently, an even-aged stand will be simulated with three diameter classes where the diameter of the first class may be, for example, 20.3 cm, the diameter of the second class 20.4 cm and the diameter of the third class 20.5 cm. The same code and set-up allows to simulate, in the same simulation, an uneven-aged stand for the same PFT but in a different pixel with, for example, the smallest diameter 7 cm, the medium diameter 25 cm and the largest diameter 45 cm.

Dynamic self thinning coefficient

Prior to rxxxx the parameter describing the site quality in the self thinning relationship (i.e., alpha_self_thinning) was PFT-specific. This parameter is used in the calculation of the stand density and thus stand biomass. In combination with upper and lower values for RDI, it contributes to the standing biomass and thus the C-sink. From rxxxx the option TEST_DYNAMIC_SELF_THINNING has been added (default = FALSE). The user can now chose between a fixed PFT-dependent value for ALPHA_SELF_THINNING or to calculate a pixel and PFT-dependent value for ALPHA_SELF_THINNING. The pixel and PFT dependent value is calculated as a function of a reference GPP.

At the first time step the reference GPP is prescribed from the parameter file. Subsequently the reference GPP is updated during the spinup (but not during a transient or historical run). At the first time step the initial value for alpha_self_thinning is prescribed from the parameter file. Each update of the reference gpp will also result in an update of alpha_self_thinning which in turn affects gpp. The spinup is now needed to reach an equilibrium between alpha_self_thinning and gpp. The initial values of the reference GPP and alpha_self_thinning will determine the outcome of the spinup.

Dynamic tree height

Prior to r8316 the parameter describing the tree height for a tree with a diameter of 1 m was PFT-specific. This parameter is used in the calculation of tree height from tree diameter. From r8316 the option TEST_DYNAMIC_HEIGHT has been added (default=FALSE). The user can now chose between a fixed PFT-dependent value for PIPE_TUNE2 or to calculate a pixel-dependent value for PIPE_TUNE2. The pixel-dependent value is calculated as a function of precipitation. The underlying approach and parameters are documented in src_parameters/dynamic_parameters.f90. The functionality can be activated by setting TEST_DYNAMIC_HEIGHT=y. Note that this functionality has been coded as a hack because we need to decide which approach we want to use. Only one of the options should be maintained. When running the code with or without TEST_DYNAMIC_HEIGHT the temporal dynamics of tree height vary quite a lot. As a consequence PIPE_TUNE3 is very different. constantes_mtc now defines a pipe_tune3a and pipe_tune3b for the two options. Parameters for the height-precipitation relationship has been fitted on GEDI by making use of CRUJRA. Although these parameters should not be changed/tuned unless the fitting is redone, these parameters (ALPHA_HEIGHT_PRECIP and BETA_HEIGHT_PRECIP) have been externalized.

In r8316 the initial value of precip_longterm was set to 0 and the initial pipe_tune2 value came from the parameter min_pipe_tune2 (7m). In r8318 the initial value of precip_longterm was set to 400/365 (in mm day-1) which results in a pipe_tune2 of ~22m. min_pipe_tune2 is now redundant. Following more testing it needs to be decided which approach should be used to initialize this functionality.

Experiments

Describes r8017. The model comes with the following well-tested offline experiments found in folder config/ORCHIDEE_OL. Well-tested means that these set up of experiment have been tested on different servers and some of the experiments are part of the trusting:

• FG_CRUJRA_SPIN: 2x2 degrees CRU-JRA forcing cycles over 1901-1910 between 1901 and 2240 (​about the forcing). Start from scratch. 15 PFTs, no land cover changes, nitrogen input for 1850, 1860 or 1900 depending on the species, and fixed CO2 concentrations (​details and context).
• FG_CRUJRA_TRANS: 2x2 degrees CRU-JRA forcing cycles over 1901-1910 between 1861 and 1900 (​about the forcing). Restart from FG_CRUJRA_SPIN. 15 PFTs, no land cover changes, nitrogen input for 1900, and annual CO2 concentrations.
• FG_CRUJRA_HIST: 2x2 degrees annual CRU-JRA forcing between 1901 and 2010 (​about the forcing). Restart from FG_CRUJRA_TRANS. 15 PFTs, annual land cover changes, annual nitrogen input and annual CO2 concentrations.
• FR_SAFRAN_SPIN, FR_SAFRAN_TRANS, FR_SAFRAN_HIST: experiement set up for more detailed studies over France.
• FG_CRUJRA_SECHIBA_OD_HIST: 2x2 degrees annual CRU-JRA forcing between 1901 and 2010 (​about the forcing). Start from scratch. 15 PFTs prescribed by reading in a biomass map. Land cover changes, input deposition, CO2 concentrations, and all other biogeochemical and ecological processes are deactivated. This configuration is under development, waiting for a global LAI map (in case of CN-CAN this has become biomass map) to be tested.

Experiments that are used but only on specific servers:

• SPINUP: used for single pixel FLUXNET runs.
• ENSEMBLE: used for single pixel FLUXNET runs (more details).

Experiments that are no longer used or removed:

• OOL_SEC_STO_FG3 (removed in r8012): 0.5x0.5 degrees annual WFDEI_GPCC forcing between 1979-2009 (​about the forcing). Start from scratch. 15 PFTs, annual land cover changes, annual nitrogen input and annual CO2 concentrations.
• OOL_SEC_STO_FG3nd (removed in r8012): 0.5x0.5 degrees annual WFDEI_GPCC forcing between 1979-2013 (​about the forcing). Start from scratch. 15 PFTs, annual land cover changes, annual nitrogen input and annual CO2 concentrations.
• FORCESOIL: was restored for ORCHIDEE trunk 3 but was not maintained for ORCHIDEE trunk (see ticket #684).
• SPINUP_ANALYTIC_FG0(nd): removed in r7671 because better results were obtained by using recruitment for unmanaged forests in the spinup (​Read more about the new challenge to spinup a model with stand structure and abrupt mortality.). SPINUP_ANALYTIC_FG0(nd): 2x2 degrees CRU-JRA forcing cycles over 1901-1910 between 1801 and 1900. Start from scratch. 15 PFTs, no land cover changes, nitrogen input for 1850, 1860 or 1900 depending on the species, and fixed CO2 concentrations. During this pre-spinup each forested PFT is clear cut once. This makes the age structure more heterogeneous. Clear cutting is prescribed according to a map such that the process is reproducible. When FG0 was used, the final files of FG0 were to be used as restart files for FG1.
• OOL_SEC_STO_FG4 (can be found in svn before tag 4.1 or r7686): ~0.5x0.5 degrees annual CRU-JRA forcing between 1901 and 2010 (​about the forcing). Start from scratch. 64 PFTs, European forest are defined at the species level with 4 age classes, forests outside of Europe are defined at the MTC level with 1 age class, annual land cover and tree species changes, annual input deposition, annual CO2 concentrations, annual forest management, and annual litter raking.
• OOL_SEC_STO_FG5 (can be found in svn before tag 4.1 or r7686): 1x1 degrees annual IPSL RCP 4.5 forcing between 1911 and 2100. Start from OOL_SEC_STO_FG4. 64 PFTs, no land cover and changes, annual input deposition, annual CO2 concentrations, prescribed species and management changes following annual a stand replacing disturbance, litter raking for 2010.

Forced clear cut

Describes r6614. FORCED_CLEAR_CUT is a flag used to force ORCHIDEE trunk to clearcut a forest after one year of simulation. This flag is set to TRUE in order to start a new stand at the beginning of the FIN step in ENSEMBLE runs. It helps us to control the stand age at the end of the HIST step. If you want to use this flag for other purposes, do not forget that a clearcut means that the majority of the living biomass is removed (circ_class_biomass for sapwood and heartwood), but the other pools are transferred in the litter pool (leaf, branches, fruit, fine root). Note that if FORCED_CLEAR_CUT is used, the model will clearcut at the end of every year in that run. The typical set-up should be: 300 years of spin-up with FORCED_CLEAR_CUT set to FALSE, 1 year with FORCED_CLEAR_CUT set to TRUE, a simulation with the length similar to the age of the forest with FORCED_CLEAR_CUT set to FALSE.

Note that the clear cut happens on the last day of the year. If monthly or annual history files are used it is not possible to see whether the clear cut actually happened because XIOS outputs 30 or 364 normal values and only 1 value after the clear cut. The clear cut should be clearly visible in the first year of the follow-up simulation.

Forest management and management changes

Describes r7451. 70% of the global forests are managed, which contradicts the assumption in previous versions of ORCHIDEE that forests are long-lived natural vegetation types. Forest management, inspired by ORCHIDEE-FM, was implemented in ORCHIDEE-CAN. Owing to the allometric allocation scheme, the introduction of diameter classes and a canopy structure, one of the key principles of ORCHIDEE-FM were retained. This principle is reflected in the rule of Deleuze and Dhote that allocates carbon to different diameter classes based on the basal area of the tree. Tree growth is followed by relative-density index (RDI)-based management which enforces thinning and harvest operations based on the current tree density, RDI and the maximum density which in turn is based on the observed (large-scale) self-thinning density. Good to know: the forest management approach is embedded in the carrying capacity of a stand which has been formalized through the self-thinning relationship which makes use of two parameters alpha_self_thinnin and beta_self_thinning. As a fail-safe option the longevity of a stand is still defined but will only be used when all other criteria fail to kill the stand (not observed). Longevity is defined by the parameter residence_time.

The forest management strategy can either be forced as a single value to all PFTs and grid cells, or be read from an input map to allow for spatially and temporally varying strategies. If the forest management strategy is not specified the default value "unmanaged" (FM = 1) is used. This implies that the stand is never thinned or harvested. Mortality is governed by an RDI which in turn is a function of quadratic mean stand diameter. Once the stand density drops below the threshold or the tree diameter exceeds a different threshold, a stand replacing disturbance occurs and a new stand is prescribed in the next time step. When developing and testing the model, a single forest management strategy can be applied for all pixels and PFTs. Set read_fm_map to n and specify the desired management strategy (1-4) through forest_managed_forced. ORCHIDEE trunk distinguishes 4 different strategies:

• FM=1 unmanaged
• FM=2 high stand management: with RDI based thinnings and density/diameter based final harvest
• FM=3 coppice
• FM=4 short rotation coppice with willow or poplar

For applications that focus on forestry or require landscape heterogeneity, a PFT-specific management strategy can be read from a spatially explicit map. Thus, the same PFT in different pixels can be assigned a different management strategy. However, within a pixel a single PFT can only have one management strategy. ORCHIDEE now has a rather detailed management reconstruction for Europe and a rather coarser global management reconstruction. Set read_fm_map to y and specify the location of the forest management map in COMP/stomate.card. Check OOL_SEC_STO_FG4 for an example of how the existing map should be defined and used.

The details of each of the 4 management strategies can be refined through a set of PFT-specific parameters. Note that not every management strategy makes use of all parameters. For more details see the SI of Naudts et al 2015 (last table). The different management strategies require parameter values for : harvest diameter max_harvest_dia, coppice diameter coppice_diameter, rotation length src_rot_length, number of rotations src_nrots, fuelwood diameter fuelwood_diameter and 4 parameters specifying the upper RDI and 4 parameters specifying the lower RDI: a_rdi_xxx_yyy (intercept), b_rdi_lower_yyy (slope), c_rdi_xxx_yyy (upper limit) and d_rdi_xxx_yyy (lower limit) where xxx stands for lower or upper. This set of 8 parameters is defined for yyy = man (ifm = 2, thin and fell) and yyy = unman (ifm =1, unmanaged). According to economic theory, high-stand forest are harvested when the actual growth drops below the long-term growth. This has been implemented in ORCHIDEE trunk. This feature was found to be very sensitive to the time frame for which actual increment was calculated. This option can be by-passed by setting this period unrealistically high, for example, n_pai =1000. Persons interested in further testing/developing this feature should set this parameter (unit: years) to 5 or 10 and go from there. Following a disturbance (which could be a clear cut), tree species changes and forest management change can be prescribed or read from a map in ORCHIDEE trunk. Set lchange_species = y, read_species_change_map = y, and read_desired_fm_map = y and specify the paths of those maps in the COMP/stomate.card. This functionality replaces the DGVM in areas where humans rather than nature govern species distribution, for example, Europe. Note that there are some constraints on the possible management changes. If the species is a conifer, for example, the new management strategy cannot be coppicing (fm=3). Anthropogenic species change has not been developed to work together with land cover change (this would require some good bookkeeping for veget_max). Forest management change has not been tested together with land cover change but because they affect different variables, i.e., forest_managed and veget_max combining both functionalities seems not overly complex.

Forcing files by ORCHIDEE

Describes r7912. ORCHIDEE can now write an output file that can be used as as a forcing file. This functionality is intended to be used in the coupled mode where ORCHIDEE would then output the simulated climate variables in a format that the ORCHIDEE offline driver can read. This functionality is disabled by default. To activate it go to src_xml/file_orchidee.def and search for Forcing_by_ORCHIDEE, set "enabled" to ".TRUE.". Add this file to one of the cards in the COMP folder (it should be orchidee.card when using the coupled model). Note that this functionality uses several implicit assumptions on the data format expected by the ORCHIDEE driver. Also, correct reading of the forcing files relies on the netcdf attributes coordinates and cell_method, the values of these attributes are set by XIOS. So, changes in XIOS could break this functionality.

Grasslands

Describes r6614. In terms of grassland phenology ORCHIDEE trunk follows branch 3 as closely as possible. For the moment ORCHIDEE trunk is also using Nicolas's patch that enforces senescence if too many of the leaves are in the oldest leaf age class. This patch is controlled by the LNVGRASSPATCH flag, and is set to T in the standard run.def.

Part of the turnover goes into the harvest pool (harvest_fraction is now hard coded to 0.5 for sapwood and leaves and 0 for roots in stomate_turnover.f90 - see ticket #687). This simple parameter should be replaced by a grassland management module but for the moment 50% of the turnover goes into the harvest_pool (just like wood) and its decomposition is accounted for in the fluxes from the product use pool.

Harvest

Describes r6614. All biomass harvest is recorded in a harvest variable harvest_pool, this variable also stores the mass and dimensions of the harvest/mortality (absolute numbers thus accounting for the pixel area!). Related variables were introduced: harvest_type, harvest_cut, and harvest_area. Wood product pools and fluxes from wood and biomass decomposition are calculated in a separate module dedicated to wood use. The dimension of the wood use pools were externalized and can be changed to better address regional differences when aiming for regional simulations. The default 1, 10 and 100 year pools were replaced by 2, 17 and 50 years which is closer to the reality for Europe. For most parts of the world a 100 year wood pool is very optimistic.

Height (of vegetation)

Describes r7451. Vegetation height is dynamic; as the PFT (in particular forests) grow, the height will grow through the allocation. When multiple circumferences classes are used, height has an impact on where the carbon goes. Due to intra-stand competition, more carbon gets allocated to taller circumference classes. Disturbances can change the average stand height. Self-thinning is a natural process in which small trees die because they are out-competed for resources. Environmental mortality is another natural process in which larger trees are killed as they are assumed to be more sensitive to environmental stress. Windthrow is an (optional) process which will disproportionally impact taller trees. In runs using forest management (not currently the default), forests can be thinned from above or below. Forests thinned from above lose tall trees to harvest, while those thinned from below lose small trees. Both actions change the average height of the stand. The parameters pipe_tune3 and pipe_tune2 link tree diameter to tree height. Where pipe_tune2 is the tree height (m) for a tree with a diameter of 1 meter. Those parameters determine the diameter-height relationship but the height growth is mostly driven by stand mortality. Many trees need to die early on during stand development to reach realistic height growth rates. That is the reason why the RDI of young stands is rather low (0.3 to 0.4 -> determined by d_rdi_xxx_yyy) Vegetation height has an impact on roughness height (aerodynamic properties of the surface). Two variables are output for each PFT: "HEIGHT", which is the average height of the vegetation (averaged across all circumference classes), and "HEIGHT_DOM", which is the dominant height (the height of the largest circumference class). A switch was introduced in r7083 which allows one or the other to be used in the calculation of the roughness height.

Hydraulic architecture

Describes r7829. A new numerical solution for hydraulic architecture has been implemented in the hydrol.f90 module of ORCHIDEE trunk. This hydraulic architecture is controlled thanks to a set of flags. The flag ok_hydrol_arch imposes the model to solve a hydraulic architecture (note that energy_control should be set to 1 to control this). By default, the hydraulic architecture developed in ORCHIDEE-CAN is solved (in hydraulic_arch.f90 module). To use the hydraulic architecture developed following Tuzet et al. (2017) model, another flag needs to be set : is_tuzet_hydrol_arch. When this hydrualic architecture is set up, several options can be added.

Two different ways of solving root absorption are implemented. One with low calculation time using a dynamic resistance or another one solving radially Richards' diffusion equation around the roots. The flag ok_tuzet_hydrol_arch_muff permits to activate the radial resolution around the roots (by default, it is set to False and the dynamic resistance is used).

This hydraulic architecture also permits to represent water storage in vegetation. This process is controlled thanks to the flag ok_tuzet_hydrol_arch_storage. By default, the flag is set to False and no storage pools are modelled.

Following discussions, it has been decided to also implement dynamic hydraulic resistances in the vegetation following Yao et al. (2022). Those dynamic resistances (which for now lead to too strong instabilities) are controlled thanks to the flag ok_tuzet_hydrol_arch_dyn_resist. By default, the flag is set to False and fixed parameterised resistances are set.

Finally, for soil resolution purposes, two additional flags have been set up. The first one, split_soil_properties along with the parameters mcr_sup, mcr_inf, mcs_sup, mcs_inf permits to impose different soil properties in the two soil horizons modelled in the hydraulic architecture (this is only used in the hydraulic architecture and does not change anything in hydrol vertical resolution). By default, the flag is set to False and the soil properties are the same in both horizons and are equal to the ones used in hydrol. The other one, is_vg permits to change the way soil water potential and soil hydraulic conductivity are calculated. By default, it is set to False and the hydraulic architecture follows the calculation in Campbell (1985). If it is set to True, the hydraulic architecture will use Van Genuchten (1980) calculation.

Land cover change (with age classes)

Describes r6614. Land cover change now accounts for age classes. It is controlled by veget_update. Set veget_update = 0Y if land cover change should be disabled. The wood pool and its subsequent fluxes were moved from the land cover change routine to a separate routine. Furthermore, land cover change also deals with the change of biological land uses to non biological land uses (of which the most important change is probably urbanization). If urbanization happens, all the carbon an nitrogen are stored in a series of variables burried_xxx where xxx stands for a different pool, e.g., litter, soil, .... Burried_xxx are cumulative variables thus increasing over time . There is a place holder in sapiens_lcchange.90 to also develop the release of the buried carbon and nitrogen following de-urbanization (see ticket #616). The series of the burried_xxx variables are not yet written to an output file but this could be easily added (they are already defined in the xml files).

An interesting parameter is min_vegfrac. When reading in a land cover map, PFTs with a fraction below min_vegfrac are removed. Likewise the fraction cover of a PFT after a land cover change should not be less than min_vegfrac either. This requirement seems to have been solely established to avoid ending up with too many PFTs with very small fractions. Because the the non-biological and biological fraction covers of each pixel should sum up to one, removing even these very small fractions implies that these fractions need to be added to one of the remaining PFTs. First it is tried to add the fraction to the bare soil (this will only be accepted if the new fraction of the bare soil exceeds min_vegfrac), then the code tries to allocate the residual fraction to the largest vegetated fraction. If age classes are used this should be the largest vegetated fraction in the first age class of a PFT. If all of the above failed, the residual fraction is added to frac_nobio irrespective of whether frac_nobio exceeds min_vegfrac. Everytime this happens, the failure to meet the min_vegfrac criterion is registered in the variable failed_vegfrac. This variable is not yet added to an output file.

Note that the min_vegfrac criterion could be the reason of why very small land cover changes occur. Another consequence is that the land cover fractions in the model are not exactly the same as those read in from the maps. Deviations should remain small and should not accumulate over time. Assume that in y0 the fraction of PFT2 = 0. In y1 the map tells us the fraction is half of min_vegfrac. The model will keep the PFT fraction to zero. The model and the map will no longer be in line with each other. In y2 the map tells us the fraction is twice min_vegfrac. The model will now accept the change. The model and the map will be in line with each other.

HACK_VEGET_MAX_NEW can be used to prescribe veget_max from run.def rather than a PFT map. This is a useful feature to debug simplified test cases with land cover changes. One possible approach to set up a simple test case is to run ORCHIDEE for 1 year (but could be any length depending on your needs) with IMPOSE_VEG = y, VEGET_UPDATE = 0Y. Set the SECHIBA_VEGMAX according to your test case (for example SECHIBA_VEGMAX06=0.5 SECHIBA_VEGMAX10=0.5 all others set to 0). Next set-up the model to start from the restart of the previous set-up. This new set-up should not be longer than 1y (the model should run for longer than 1 y but the hack will be implemented every year). Set IMPOSE_VEG = n (this is needed to make veget_update work), VEGET_UPDATE = 1Y (this is needed to make HACK_VEGET_MAX_NEW work), HACK_VEGET_MAX_NEW = y and now set FRAC_NOBIO_NEW and VEGET_MAX_NEW to your needs. For example VEGET_MAX_NEW06 = 1 and VEGET_MAX_NEW10 = 0.0 (note that VEGET_MAX_NEW is defined at the PFT level). For more complex test also FRAC_NOBIO_NEW can be set (defined at the pixel level). In the example above the 0.5 fraction of grasslands (PFT10) will be converted into summergreen deciduous forest (PFT6).

Leaf area

Describes rXXXX. There is no longer a parameter for maximum LAI in the trunk. Given that LAI is a key variable that links the biogeochemistry to the biophysics in ORCHIDEE it deserves some background information. LAI is now calculated as a prognostic variable and is the net result of growth, allocation and turnover. There are few parameters in stomate that do not affect the LAI but three parameters stand out:

• Growth is largely determined by GPP and Ra. (a) GPP is mainly driven by the LAI so this results in an interesting feedback. If the LAI is very low, GPP will be low, growth will be low, LAI will remain low, etc. When GPP is too low to make the LAI, the vegetation will die. This could be rather quick, i.e., within a year, or it could be rather slow, i.e., several years. Given this is a feedback loop there is no single parameter that can be used to control this feedback. (b) The ratio of NPP/GPP or the fraction of GPP that is used in autothropic respiration and especially maintenance respiration will to a large extent determine how much C is availble for growth. If the LAI is very high or low, always check the NPP/GPP ratio. Ra and thus the NPP/GPP ratio can be tuned through the MAINT_RESP_SLOPE_xxx variables and the COEFF_MAINT_INIT.
• Allocation distributes the growth across the different plant components. The parameter that controls the fraction that ends up in the leaves is K_LATOSA_xxx. This variable can be used to tune LAI but the effect of increasing K_LATOSA_xxx decreases for larger values (possibly due to N-limitation). The range of observed K_LATOSA_xxx is quite large. The model shows a nice sensitivity within the observed range. Note that allocation is a "zero sum game", hence, allocating less C to other components will result in more C going to the leaves. Note that the leaves (GPP) and roots (N-uptake) are necessary to grow in the first place (see section Leaf and roots).
• Turnover will decrease the leaf area. The parameter that controls the turnover is LONGEVITY_LEAF. For evergreen PFTs the relationship between longevity_leaf and turnover is straightforward. For deciduous PFTs the relationship is more confusing (see section Leaf longevity and turnover). Turnover also affects the leaf age and young leaves have a higher VCMAX than old leaves. The age-effect on vcmax is strong enough to make that a small young canopy may have a higher GPP than a large old canopy.

Leaf area index map

Describes r6614. Four flags have been identified that control the model behavior in terms of lai: ok_stomate, ok_pheno, impveg and read_lai. There is a 5th implicit flag which is whether restart files are used or not. If a restart file is used, the lai values will come from the sechiba restart file which is read at t=48. Given that each flag can take two values, we have 32 configurations in total. Out of these 32 configurations 10 are defined of which about 5 to 7 seem to be intended (for more details see Start and restart - Table 1). Many of the remaining 22 settings are inconsistent (i.e. running stomate to calculate a lai and reading an lai_map to prescribe lai), duplicate other settings, or would require further developments to work properly. Furthermore, the current code does not stop or warn when inconsistent settings are selected. The table (see Start and restart) proposes a scheme with 2 flags which can run with our without restart files, thus resulting in 8 different ways to control the lai in sechiba or the initial lai in stomate. The remaining 2 combinations are inconsistent and will stop the model.

In previous ORCHIDEE trunks, canopy structure is prescribed by a single variable called lai and the assumption of a turbid medium (Lambert-Beer). Consequently reading an lai value suffice to prescribe the entire canopy. In ORCHIDEE trunk, however, canopy structure has become a 3D property that can be calculated from the leaf biomass, stem biomass, the number of individuals and the assumptions that the trees follow a Poisson distribution in the horizontal plain, that the crowns are spherical and that the leaf biomass is uniformly distributed within the crowns.

The functionality to simply prescribe an lai value (either through impose or by reading a map) will need to be replaced by functionality that prescribes or reads the biomass, leaf age, and number of individuals. Given that the current lai map is based on a previous ORCHIDEE run, the same approach could be used to generate a spatially explicit canopy map that contains biomass, individuals and leaf age distribution. Nevertheless, reading an observed lai, i.e., from MODIS, and using it to force ORCHIDEE trunk would require a substantial number of assumptions to turn an aggregated 1D lai value into a disaggregated 3D canopy at the PFT level. For the moment the code (slowproc.f90) contains a simplified place holder that only works with NCIRC=3 and it doesn’t read leaf_age.

See also the item on "Stomate bypass"

Leaf longevity and turnover

Describes r7767. Longevity_leaf is a parameter for the longevity of a typical leaf/needle at the average temperature for that PFT. For PFTs where a large range in longevity has been observed, the large range is accounted for in the calculation of leaf_age_crit. Leaf_age_crit is thus a location-specific critical age of a leaf. For example, longevity_leaf = 220, leaf_age_crit = 165 at 25C, 220 at 15C and 330 at 5C. If the age of a leaf/needle exceeds half the value for leaf_age_crit, leaf, root and sapwood turnover will start (thus at 110 days, in example for the temperate zone, which is within the observations). Also the critical age drives the efficiency on NUE and thus the photosynthesis per unit of leaf area. Moreover, leaf growth may continue at the same time as leaf turnover(depending on plant_status). The mean leaf age is the outcome of the interplay of the leaf turnover and the leaf growth.

For evergreen PFTs turnover from ageing is easy to understand as it is the only source of turnover. A longevity_leaf of 730 days will result in a two year canopy if the PFT is growing correctly. For evergreen PFTs, the control of leaf_age_crit on the simulated longevity (reflected in the mean leaf age) is thus straightforward. Because of the way leaf_age_crit is calculated the simulated longevity differs from the prescribed longevity but the spatial patterns make sense and follow the temperature gradient that was imposed through the temperature dependency of leaf_age_crit.

For deciduous PFTs, the control of leaf_age_crit on the simulated longevity (reflected in the mean leaf age) is not at all straightforward. longevity_leaf affects the calculation of leaf_age_crit and once half the value of leaf_age_crit is reached it will determine the magnitude of the turnover but the age of the leaves is mostly driven by the phenology (senescence). For deciduous species turnover from ageing can thus be understood as the leaf turnover between plant_status = icanopy and plant_status = isenescent. The contribution leaf turnover from ageing can be monitored by two variables: LEAF_TURN_AGEING_c and LEAF_M_MAX_c. When expressed as a percentage, the within season turnover varies between 0 and 30% for trees and reaches 100 to 150% for grasses. This seems to be acceptable. At large spatial scales the relationship seems to be the inverse of what is expected, i.e., leaves live longer in warm compared to cold regions. It is not clear yet what is driving this observation. The impact of NUE through leaf_age_crit seems small, it could be that in warm regions rm/gpp is higher resulting in an early start of presenescence which would prevents the growth of new leaves. In colder regions rm/gpp could be lower resulting in a later start of presenescence and thus a longer period in which new leaves are grown. Adding new leaves decreases the mean age of the leaf pool.

The current implementation of leaf phenology and senescence would benefit from some conceptual cleaning. Although all processes make sense when considered in isolation, their combination is likely overly complex and may implement the same controls in different places in the code. Longevity describes the time a leaf can grow (but is a single parameter for the whole PFT), the rm/gpp accounts for longevity from a plant perspective, leaf_age_crit accounts for longevity from a climate perspective, and also phenology and senescence account for leaf longevity from a climatic perspective. Moreover, There are alternative approaches of simulating phenology and especially senescence (Rm/GPP, climatic, leaf age, mixed). This abundance of approaches doesn't help to understand the ecology and physiology implemented in ORCHIDEE. See ticket #876 Turnover will result in an imbalance between the vegetation components and this balance will be restored in the allocation routines. Leaf turnover will thus replace old leaves (with a low NUE) by young leaves (with a high NUE). Despite the fact that leaves are being killed, the GPP may increase and eventually the LAI may even increase.

Litter decomposition

Describes r6614. After large-scale dieback events (with a closed n-cycle, i.e., impose_cn = n), so much soil mineral N becomes immobilized to decompose litter that too little N is left for plant regrowth. To address this, we implicitly represent the action of fungivores, which eat the decomposing fungi and release N for the plants and increase N turnover rates. We set aside a fraction of qd (stomate_litter.f90) which becomes available for plant uptake in nitrogen_dynamics. This fraction is calculated and is at its maximum when the litter pool is large compared to the biomass pool. The fraction is at its lowest when the living biomass is high compared to the litter biomass. The implemented principle mimics a Lokta-Volta dynamic where the predator are the fungivores and the prey the fungi. The share of the N contained in the decomposing fungi that is released as an excrement from the fungivores ranges between 0 and 1 and is calculated.

After a die-back (with a closed n-cycle, i.e., impose_cn = n) too little N is available to reach the target C:N ratio of the receiving soil pool, meaning that carbon litter cannot decompose. If the C:N ratio is not controlled it spins out of control reaching values of more than 600. The target C:N ratio, which is a variable in stomate_litter.f90, can now vary between 1 and max_cn where the latter can be set in the run.def. Good results were obtained by setting max_cn to 250. max_cn is bound between 1 and 400, the C:N ratio of wood which is the pool with the highest C:N ratio (and biomass that contributes to the litter). During the simulation the target C:N ratio varies and a couple of decades after the dieback reaches the expected value between 5 and 20. The ecological justification behind a max_cn above the C:N ratio of the fungi is that mycorrhizae receive their C from the plant photosynthesis products. Thus, litter can still decompose even though C is not incorporated into the plant biomass. max_cn can be set in the run.def:

MAX_CN = 250

ORCHIDEE-CN allows negative n_mineralisation values to be passed to stomate_soilcarbon where they are treated as immobilization. If we have persistent negative n_mineralisation, however, we can reach a situation where we completely deplete our soil_n_min pools to satisfy the demand for immobilization. The model cannot recover from such a condition. This tends to happen when there is no plant recruitment, in which case a generation of trees will all reach their maximum diameter and die at the same time. This creates an enormous amount of litter, meaning our som_input values become relatively large and when we subtract som_input from n_mineralisation, the values become negative. To remedy this, we'll truncate som_input based on the current mineralisation pool, but we'll also limit the occurrence of negative values. We'll set a minimum min_n for n_mineralisation so we can keep some N in the soil (and so the simulation can recover). This problem was discovered and fixed before adding the fungivores and max_cn. It is not clear whether this fix is still needed but it looks like a good problem-catcher. min_n can be set in the run.def:

MIN_N = .0001

Litter raking

Describes r6614. Tree litter was collected from the forest and used in the winter in stables instead of straw. In spring the litter and manure was spread on the croplands. This lateral flow of C and N between PFTS in the same pixel can be accounted for in ORCHIDEE trunk by setting use_litter_raking = y. If litter raking is to be used, the model will search for annual maps. The path of these maps needs to be specified in COMP/stomate.card. Litter raking maps were prepared for Europe. Unless liter raking is your research topic set use_litter_raking = n. An example of to set up the model to make use of the historical litter raking maps can be found in config/ORCHIDEE_OL/OOL_SEC_STO_FG4

Litter resistance to evaporation

Describes r6783. The default for the flag do_rsoil has changed to TRUE. When set to TRUE the flag reduces the bare soil evaporation because it accounts for the moisture content of the litter layer (Sellers et al., 1992). The presence of litter always acts as an extra barrier (resistance) to bare soil evaporation. If the litter layer dries out the resistance will even further increase. As in previous version, ORCHIDEE calculates the bare soil fraction as PFT1 + SUM(veget_max_i - veget_i) where i goes from 2 to the number of PFTs. veget is now calculated by using Pgap which resulted in more bare soil, the extra resistance from litter helps to keep to evapotranspiration in check.

Mass balance closure

Mass balance closure is checked for carbon and nitrogen for all subroutines called from stomate and below. Depending on the setting of ERR_ACT the model will be stop runnin, a warning will be generated or the mismatch in the mass balance will only be written to the history files.If the model can be run without mass balance errors, we can be confident in the calculations but this does not imply that we can have the same confidence in the history files because some apparent and real mass balance problems may be introduced in between the calculation of a variable and the variable being written to the history files. Apparent problems come from : (1) temporal averaging by XIOS, (2) unit convertions (especially when the area of the pixel has to be taken into account). Apparent problems are detailed in ticket #798. Real problems can come from sending variables to the history files that are not included in the actual mass balance checks.

Mortality

Describes r6614. ORCHIDEE trunk distinguishes 2 types of natural mortality: (1) explicitly considering mortality from disturbances and self-thinning, and (2) implicitly considering background mortality. Ideally, approach (1) should be further developed such that all underlying agents driving background mortality are represented in the model (i.e., gap-scale mortality, pests, disease, windthrow, etc.) such that it can replace approach (2). Two options of background mortality may be chosen: constant background mortality and dynamic background mortality. To use the first option, set the flag constant_mortality = y. The background mortality of a forests is calculated as a constant, prescribed fraction. In trunk, this fraction is given by residence_time (see also forest management). Otherwise set constant_mortality = n, the dynamic background mortality of a forest is a function of its net primary production (npp). If npp decreases, mortality will increase Both options have been developed but only constant_mortality = y has been tested in ORCHIDEE trunk.

However, because of the introduction of self-thinning in ORCHIDEE trunk, constant_mortality = y became the default setting. In ORCHIDEE trunk, the total mortality is the maximum of the background mortality and the mortality from self-thinning. Only if self-thinning is absent or too low, background mortality will play a role. This approach implies that when constant_mortality = y is used in combination with self-thinning, background mortality will only play a role in the first years to decade before self-thinning starts (the latest calculations of RDI - see Prescribe - the role of the background mortality has further decreased). Despite its limited use, it represents an essential process: owing to background mortality, the number of individuals decreases, the remaining individuals grow faster and thus manage to reach self-thinning in a reasonable amount of time.

ORCHIDEE trunk calculates the number of individuals and uses this as a criterion to initiate a stand replacing disturbance. This approach, guided by the self-thinning relationship, avoids the need for a stand-level turnover time. ORCHIDEE-CN, and ORCHIDEE-CNP still make use of stand-level turnover. Note that the meaning of residence_time is very different between the CAN branch and the trunk. In the trunk biomass has no age and thus the residence time accounts for all forest dynamics including self-thinning, pests, diseases and windthrow. In the CAN branch, biomass does have an age and self-thinning is explicitly accounted for, hence, the residence time should be much higher as it only accounts for pest, diseases and windthrow. Even the latter is not exact because as long as those disturbances are small scale they are probably accounted for in the parametrization of self-thinning.

Nitrogen cycle

Describes r6614. ORCHIDEE trunk strictly follows ORCHIDEE trunk 3 where it concerns the implementation of the N-cycle. Following mass balance problems caused by negative N mineralization and followed by negative immobilization, the code has been slightly adjusted to ensure mass balance closure. First the flag stomate_ok_ncycle needs to be set to y, to run the model with a N-cycle. Subsequently the parameter impose_cn is used to control the N-cycle calculations. If set to y, C/N ratios are calculated but whenever N appears to be limiting, it is taken from the atmosphere to satisfy this need. This is the preferred setting when testing/developing the code without a proper spin-up. N-limitation will only be accounted for when setting impose_cn = n. With this setting the N-cycle is closed (checked when checking for mass balance closure) it requires a spin-up to produce reasonable results.

The paths of the N-inputs from atmospheric N-deposition, fertilization, and biological nitrogen fixation (Nammonium_FILE, Nnitrate_FILE, Nfert_FILE , NManure_FILE and Nbnf_FILE) are set in the stomate.card. Moreover, the N inputs you wish to include in your simulation can be specified in the run.def. If set to NONE, the given N input files will not be read.

Nammonium_FILE= ndep_nhx.nc
Nnitrate_FILE= ndep_noy.nc
Nnitrate_VAR=noy
Nammonium_VAR=nhx

Nfert_cropland_FILE = nfert_cropland.nc
Nfert_cropland_VAR = nfer
Nmanure_cropland_FILE = nmanure_cropland.nc
Nmanure_cropland_VAR = Nmanure

Nfert_pasture_FILE = nfert_pasture.nc
Nfert_pasture_VAR = Nfer
Nmanure_pasture_FILE = nmanure_pasture.nc
Nmanure_pasture_VAR = Nmanure

Nbnf_FILE= bnf.nc
Nbnf_VAR= BNF_MGN_PERM2_PERYR

NINPUT_UPDATE=_AUTO_
INPUT_SUFFIX_YEAR = n

The following summarizes the soil dynamics of nitrogen in both ORC3 and TAG4.1 (r7761), in part based on ​https://doi.org/10.5194/gmd-12-4751-2019 and in part based on the source code. Both nitrification (conversion of NH4+ to NO3-) and denitrification (conversion of N03- to N2) are present. Only four species are considered: nitrate (NO3-), nitrogen oxides (NOx), nitrous oxide (N2O), and nitrogen gas (N2). The pathways are impacted by pH, temperature, soil moisture, and soil carbon (which dictates denitrifier bacteria biomass). N2O is produced during both processes, although it is also consumed during denitrification. If oxygen availability is high, the fraction of anaerobic microsites is low and nitrification is reduced; oxygen availability is limited by oxygen diffusion from the atmosphere as well as oxygen lost to respiration of soil carbon. Oxygen does not appear to impact denitrification.

In TAG4.1, ammonium in the soil can be lost through leaching and adsorption onto clays, while nitrate can be lost through leaching. Leaching is mentioned as an output of nitrogen in ​https://doi.org/10.5194/gmd-12-4751-2019, but not which species. Clays are not mentioned, though both loss mechanisms are present in r7104 of ORC3.

Nitrogen reserve pools

Describes r6904. As for carbon, the model simulates two non-structural nitrogen pools: labile and reserve. Nitrogen uptake is stored in the labile pool, and after the growth, nitrogen in excess of the calculated target labile pool is moved to reserve nitrogen. Before r6825, nitrogen limitation during plant growth was rather common but was solely based on the labile nitrogen pool. It was rather common that the labile pool was empty but that at the same time the reserve pool was very well filled. Plant growth experienced n-stress despite having a good amount of N in the reserves. To deal with this issue, r6825 also considers the reserve nitrogen for estimating nitrogen limitation. The parameter use_reserve_n controls the functionality along with nitrogen uptake feedback which was added r6869 (See Root nitrogen uptake).

Parameter files

Describes r6614. A model run now requires three different parameter files: run.def, orchidee.def and orchidee_pft.def. LibIGCM reads these files in COMP/orchidee_ol.card.

• The run.def contains all the information for the driver. This is also the file where the default global domain can be reduced (LIMIT_N, LIMIT_S, LIMIT_E and LIMIT_W). When the same experiment is to run in off-line and coupled mode, only the run.def should be changed.
• The orchidee.def contains the details of the experiment together with the settings of the orchidee_ol, sechiba and stomate cards in COMP. Each of the predefined configurations in config/ORCHIDEE_OL will have different settings. These files need to set up manually because each experiment will differ. Copy the best matching predefined experiment from config/ORCHIDEE_OL and use it as a template to set-up your own simulation experiment.
• The orchidee_pft.def is available in different configurations. The configuration depends on the interest of the user in vegetation and its dynamics. The default setting for the moment is the same as the trunk, i.e., 15 pfts with a single age class. Nevertheless other useful orchidee_pft.def files are available. Users can make their own files by making use of the script in the MAKE_RUN_DEF folder. Note that different numbers of PFTs will require matching land cover maps. Also the order of the PFTs in the run.def should follow the order of the PFTs in the map. There is no quality control to check whether definition files and input maps following the same PFT order. If the PFT maps are available, users can run each of the predefined simulation experiments (config/ORCHIDEE_OL) with all available orchidee_pft.def files. The number of PFTs largely determines the speed of the model.

How to use a different orchidee_pft.def file?

• Go to COMP/orchidee_ol.card and change the value of DefSuffix. DefSuffix = 15pft.1ac will make use of orchidee_pft.def_15pft.1ac, change DefSuffix = 3pft.1ac to have the model use orchidee_pft.def_3pft.1ac. You can also use your own customized orchidee_pft.def, just make sure the follow the naming guideline.
• make sure the PFT map has the correct number of units nvm vs. nvmap. These numbers are set in orchidee_pft.def but affect the expectation of the model when reading the land cover map. The correct land cover map is set in COMP/sechiba.card
• orchidee_pft.def with 3 PFTs make for a powerful test-case. Most of the output variables in stomate have their units in something per square meter of vegetation. So two identical PFTs with different vegetation fractions should still result in the same stomate outputs. If this is not the case, something went wrong with the units (and likely the way veget and/or veget_max are used in the model. To use this test case remember to set impose_veg = y (PARAM/orchidee.def) and veget_update = 0Y (COMP/sechiba.card). The test case can be run over large spatial domains which may also help to find bugs.
• orchidee_pft.def with 64 PFTs demonstrates how MTCs, PFTs, species and age classes can all be combined in a single set-up. The species parameters are limited to the European domain, the MTCs are used outside of Europe. In Europe 4 age classes are distinguished, outside of Europe a single age is used. This is the set-up that was used in Naudts et al 2016 and Luyssaert et al 2018.

Parameters - good to know

k_latosa K_LATOSA regulates the leaf mass as a function of the sapwood mass. A higher K_LATOSA will give a higher leaf mass (and thus LAI) for the same amount of sapwood. Higher leaf mass typically results in higher GPP. The literature reports a relatively large range of K_LATOSA values (500-17,000) which gives us some room to tune this parameter.

While parameterizing CAN (a former carbon only branch), the response to k_latosa appeared almost linear. More recent tests with r6818 and a focus on PFT4 showed a clear nonlinear relationship between k_latosa and annual mean LAI. LAI saturated when k_latosa increased: 2500 to 4500 in steps of 500. Around r7124 we think we understand this saturation. GPP is the input of carbon. We can only allocate to tissue if we have enough N. If we don't, it goes to reserve pool. If we have too much in the reserves, then we reduce GPP with sugar loading. Higher LAIs will initially result in higher GPP but if the system lacks N to allocate the feedback through sugar loading may reduce the GPP on the long term. The model needs a way to "burn off" the GPP that could not be allocated. Leaching of labile carbon to the soil (thought to be responsible for up to 20-30% of the NPP - textbook by Chapin 2012) comes to mind. As leaching has not been implemented yet, we could add leaching artificially by turning over the roots. So now we have high GPP, same amount of C goes to tissue, but rather than the excess going to the reserves (and decreasing GPP through the sugar loading) it turns over the roots. This could also increase the soil carbon, which is still desirable around r7113.

As a functional constraint the ratio of labile+reserves over leaf mass could be used. As long as carbon is allocated to the leaf mass rather than the labile+reserves tuning K_LATOSA makes sense. If the labile and reserves are increasing faster than the leaf mass, tuning K_LATOSA reaches its limit (i.e. its effect saturates). The described responses are more complex than what could be captured in a simple ratio.

The obvious data to constraint the tuning of K_LATOSA are LAI and GPP. The correlation between LAI and GPP is so strong that one of these data streams is probably enough. Note that ORCHIDEE 2 can downregulate GPP due to sugar loading. Plotting the variable sugar loading and the sum of labile and reserve carbon gives an idea of whether the pixel x PFT experiences N-limitations.

coeff_maint_init COEFF_MAINT_INIT regulates the autotropic respiration as a function of temperature. Increasing the parameter will result in more respiration and thus less allocatable biomass. Given that GPP-Ra = NPP, tuning COEFF_MAINT_INIT gives a direct control over the NPP/GPP ratio. A high NPP/GPP implies that the model needs to find N to convert a large share of the GPP into tissue. If this is not possible due to N-limitation, the allocatable carbon that cannot be put away into tissue will end up into the labile and reserve pools. Likewise aiming for a lower NPP/GPP implies that more of the GPP will be required, leaving less allocatable biomass and thus reducing the N-demand to build tissue. As a consequence lower NPP/GPP ratios are likely to results in lower reserve and labile pools but it does not guarantee acceptable labile+reserves (i.e., in cases with extreme N-limitation, even low NPP/GPP ratios may still accumulate C in the labile+reserves).

If NPP becomes so low that it can no longer compensate for the biomass turnover, root mass may start to decrease resulting in reduced n-uptake, likewise LAI should decrease and thus GPP. It is not clear yet whether this parameter can be tuned outside a functionally desirable range (implying unwanted feedbacks). A single pft-dependent temperature relationship makes it difficult to have reasonable NPP/GPP ratios throughout a PFT covering a large temperature gradient.
In ORCHIDEE Ra also accounts also for leaching, BVOCs, C-subsidies to mycorrhizae, etc. Hence Ra should be higher than observed Ra. The data constraint for tuning COEFF_MAINT_INIT is thus the observed NPP/GPP.

vmax_uptake VMAX_UPTAKE has a linear impact on the n-uptake, unless its value is too large so that excessive N in the plant begins to have a negative feedback on N uptake (at least this is what embedded in the model). Higher VMAX_UPTAKE values result in more n-uptake for a given amount of root mass. More plant available N implies that a larger share of the allocatable biomass can be converted into tissue likely resulting in less excess carbon which in ORCHIDEE 4.0 will be reflected in less labile+reserve carbon. This relationship was demonstrated to hold over large areas and many PFTs but there is suspicion that there is an unwanted feedback: high plant available N will also result in higher LAI which will result in higher GPP which results in more allocatable carbon and thus finally results in higher biomass. It is conceivable that the C-availability (through GPP) increases faster than the N-availability (through root uptake). Excessive N uptake might also deplete soil N pool to have a negative feedback. These points need to be confirmed/rejected.

As a functional constraint one could look at the ratio between labile+reserves over n-uptake. Higher n-uptake without relative increase in the labile+reserve carbon is desirable. A more direct way of assessing this issue is to monitor the ratio between leaf mass and root mass. If the ratio between root and leaf mass decreases considerably, tuning VMAX_UPTAKE may have reached its limits.

Parameter compensations

• When setting up factorial experiments to quantify model sensitivity, we often assume additive effect and ignore the interactions between parameters. Tests with r6818 and a focus on PFT4 demonstrated clear non-additive effects. Five sets of parameters for PFT4 were stepwise replaced with the parameters for PFT7.
  
  

set1: k_latosa, longevity_leaf, sla
set2: arjv, vcmax25, coeff_maint_init
set3: pipe_density, pipe_tune_*, tree_ff, alloc_min
set4: psi_50, psi_leaf, k_sap, r_froot
set5: leaf_psd_*, leaf_ssa_*, bgrd_ref_*, nir, vis

Changes of set1 and set1 combined with set2 to set 5 (set1+set2, set1+set3, set1+set4 and set1+set5) increased LAI much (~2.5 -> ~7). However, change of all sets (set1+set2+set3+set4+set5) showed a moderate increase in LAI (~2.5 -> ~4) which implies interactions between parameters.

Allocation and plant hydraulic architecture
When tuning allocation to match observed wood or root masses productivity or their stocks, it may be necessary to change longevity_sap and longevity_root. Note that this will affect the ratio between leaf and root mass because Mr = Ml / LF where LF = c0_alloc * KF and c0_alloc = sqrt(k_belowground(pft)/k_sap(pft)*longevity_sap/longevity_root*sapwood_density). Likewise when tuning the hydraulic architecture to match observed transpiration, c0_alloc will change and so will allocation. Note that in c0_alloc k_belowground is used whereas the calculation of hydraulic architecture uses k_root. The difference is explained in the code.

Towards formal parameter tuning This far ORCHIDEE 4.1 has been tuned manually. The following steps were used: (1) Tune K_LATOSA to obtain a reasonable LAI (2) Tune COEFF_MAINT_INIT to obtain a reasonable NPP/GPP. If needed iterate with (1) (3) Limit SUGAR_LOAD_MIN to 0.75 to make sure the model is sensitive to accumulating excess C in its labile and reserve pools (4) Tune COEFF_MAINT_INIT, MAIN_RESP_SLOPE_1, and MAIN_RESP_SLOPE_2 to obtain a good NPP/GPP ratio, a good LAI and the lowest possible number of pixels with excessive C in the labile and reserve pools. (5) Set SUGAR_LOAD_MIN to 0.05 and thus give back the feedback between excessive C and GPP. Note that this has an effect on the NPP/GPP ratio. For the majority of the pixels sugar loading should remain within the range of 1-0.75. Sugar loading should prevent the model from accumulating crazy C-reserves in few pixels (which were so crazy that they affected global numbers). (6) Fine tune k_LATOSA

In a formal optimization, the objective function should maximize the LAI while minimizing the labile+reserves while maintaining a leaf_m_c/root_m_c ratio between 0.5 and 2.0. From a plant perspective K_LATOSA should be as high as possible, COEFF_MAINT_INIT should be as low as possible and VMAX_UPTAKE should be as high as possible. Nevertheless, GPP that has to be stored as excessive labile+reserves is a sign of N-limitation and/or a disequilibrium between these parameters. A plant that is well adapted to its environment should only take up the N it can used, develop a canopy that results in a GPP that can be allocated and has a high NPP/GPP ratio as that gives a clear competitive advantage.

The combined impact of tuning all three variables may be seen in the accumulated biomass. Biomass, diameter, and tree height remain valuable data sources to evaluate the parameterization of the allocation scheme. It remains difficult to disentangle the impact of single parameters (Previously, Jina has observed some non-linear interactions between parameters. Attempts to refine the parameterization of r7334 resulted in the following insights. Understanding the changes in the absolute values is mostly straightforward but rather subtle changes of several parameters at the same time were required to change the functioning (e.g. higher gpp for the same lai).

Table. Responses and interactions between variables (columns) and parameters (row). The description assumes that the parameters values were INCREASED.

PARAMETER	LEAF_M_c (LAI)	ROOT_M_c	SUGAR_LOAD	PLANT_N_UPTAKE
K_LATOSA	K_LATOSA is used in the calculation of the leaf mass, tuning K_LATOSA will thus directly affect LAI. Increasing K_LATOSA will increase LAI up to the point where sugar_load becomes limiting	Root mass will follow the leaf mass as a fixed function of KF and c0_alloc	Increasing LAI requires more N to sustain the higher NPP, sugar_load kicks in when N-limitations occurs and will decrease GPP and thus also NPP	Uptake increases proportional to the root mass but the N-demand will increase proportional to the leaf mass
K_BELOWGROUND	Lai increases because less C is now allocated to the roots so there is more C for the leaves. GPP will increase but not as fast as LAI because there might be not enough plant available N	K_BELOWGROUND is used in the calculation of c0_alloc, tuning K_BELOWGROUND alters the ratio between leaf and root allocation. Increasing K_BELOWGROUND decreases the root mass relative to the leaf mass but the absolute root mass may increase due to the increase in leaf mass.	N-stress will most likely increase (the value of sugar_load will therefore decrease) because the root mass and decrease relative to the leaf mass. Hence, the supply of nitrogen decreases whereas the demand is increasing.	PLANT_N_UPTAKE decreases because of a decrease in root mass.
VMAX_UPTAKE	Increasing the N supply will allow growing larger canopies even if plant growth before parameter tuning was not N-limited. Interestingly these large canopies may generate lots of gpp resulting in N-limitation which will be reflected in a decreasing value of sugar loading.	Root mass will follow LAI. The absolute value of root mass will increase jointly with the increase in LAI.	N-stress will decrease (so sugar_load will increase). Nevertheless, the effect of increasing VMAX_UPTAKE was observed to level off.	PLANT_N_UPTAKE will increase up to the point that soil_n_min becomes the limiting factor. If this point is reached there will be consequences for litter and soil carbon decomposition and thus heterotrophic respiration.
COEFF_MAINT_RESP	More respiration leaves less carbon to be allocated to leaves and roots for a given LAI. Note that the absolute value of LAI will decrease.	More respiration leaves less carbon to be allocated to leaves and roots for a given LAI. Note that the absolute value of LAI will decrease.	For a given LAI the N-demand will decrease (as a smaller fraction of the GPP will be allocated towards biomass growth which requires N) resulting in less N-stress (thus a higher value for sugar_load)	Following the decrease in the allocatable C, root mass will decrease resulting in less nitrogen uptake

Phenology

Describes r6918. The pft-specific parameter always_init controls whether the phenology depends on the reserves (set to .FALSE.) or is forced (set to .TRUE.). Note that a forced phenology (thus always_init = .TRUE.) has no ecophysiological basis, it is a numerical approach to stabilize the vegetation cover. A stable vegetation cover is particularly welcome in coupled simulations but likely hides real vegetation dynamics (especially under future climate conditions) or problems in other routines or parameter settings. If a PFT keeps dying in an area where it is currently present, this would hint at a problem with the current model/parameters. If a PFT keeps dying under future conditions, it may be a real response (depending on the PFT). If phenology is always initialized, plants will develop an initial canopy in phenology irrespective of whether the plant had sufficient carbon and nitrogen reserves and for evergreen species irrespective of whether the canopy was viable at all. This setting basically overcomes a mortality event at the expense of taking up carbon and nitrogen from the atmosphere. When used in combination with impose_cn = n, an inconsistency is introduced: impose_cn = n reflect the desire to close the nitrogen cycle, always_init = y opens a backdoor in the nitrogen cycle.

From a conceptual point of view, ORCHIDEE trunk is all about vegetation dynamics and thus instabilities in the vegetation cover. In ORCHIDEE trunk there are two processes that can deal with dying PFTs including evergreens PFTs. First, ok_recruitment could used. If ok_recruitment = .TRUE. a decrease in the canopy cover will result in more light reaching the forest floor which in turn should trigger recruitment of -for the moment- the same PFT. Generations (cohorts) can take over from each other without loosing the canopy cover entirely. Second, if there are insufficient reserves to grow leaves, there will be no or insufficient gpp, the carbon reserves will be consumed by respiration processes, the plants will be killed, the biomass transferred to the litter pools and the same or another PFT (see section on species change) will be replanted. ORCHIDEE trunk was developed to work with always_init = .FALSE. so this has become the default value, contrary to previous versions of the trunk where always_init = .TRUE. is the default.

Temperature is a very consistent climatological variable in the sense that each year the summer is warmer than the winter (even in so-called cold years). This implies that temperature driven phenology is very robust in the sense that the PFT will develop a canopy each year. Soil moisture, however, is much less predictable and in very dry or very yet years the seasonal patterns may differ a lot. Especially when soil moisture criteria are combined with other criteria (pheno model: moigdd) there are years that the criteria are never met and that the PFT never develops leaves. This is not at all observed. Even farmers would plant crops when the conditions are far from optimal (typically because they have to plant the crops before they know that the conditions will be far from optimal). If set the yes, the flagok_force_pheno will start the growth of a canopy a predefined number of days (force_pheno_mtc) after the average budbreak data for that PFT at that location. Test simulations over the Americas showed that over a 30 year period forced phenology was used 5 to 10% of the years between 1901 and 1930 for PFT10, 12 and 13. This flag did little or nothing to the phenology of the other PFTs with the pheno model moigdd. The variable qc_pheno_event registers if budbreak happened during the year. The variable qc_forced_pheno_event registers if the phenology event had to be forced rather then being simulated as a function of the climatological variables. The ratio of these output variables can be used to plot in which regions and years phenology had to be forced. Not surprisingly it is most used in dry regions of the PFT distribution.

Pgap

The Pgap calculation in r7278 was checked line by line by three persons. No obvious problems were found. Why was it checked in the first place? Pgap based veget is higher than Lambert-Beer based veget. Given the differences in definitions/approach this is not really a surprise but it seemed worth to check. This line by line check revealed that the Pgap that goes into the calculation of veget assumes zenith solar angle. We have two options but both might be problematic: (a)assume zenith angle equals to 0 degrees when calculating Pgap. This woul be OK for interception and turbulence but problematic for radiation-processes, and (b) use the zenith angle at noon. This is an OK for radiative processes (but will always overestimate the gap fraction but it is underestimates the gap fraction for interception. A third solution could be to introduce two or more vegets based on different solar angles? Easy to calculate but then the problem might be shifted to sechiba: which veget to use when?

Since r789X the calculation of Pgap was changed such that the total crown volume cannot exceed the entire canopy space. In earlier version this happened quite often but the impact on porosity was rather small. The same version calculates Pgap for the local zenith angle at noon. Canopy shape was made to vary between the tropics, temperate and boreal PFTs. Also the parameter hack_pgap was added. Its default value is FALSE meaning that Pgap (is iused to calculate veget). If set to TRUE, veget will be calculated by making use of Lambert Beer.

Photosynthesis

Describes r6614. Photosynthesis is calculated as in ORCHIDEE trunk, and ORCHIDEE trunk 3 but the way the canopy levels are defined has changed. The reason of this change is in the albedo and energy budget for which physical layers are required. For this reason the space between the bottom and the top of the canopy is divided into x layers. x is set by the parameter nlevels_photo. Applications with ORCHIDEE trunk used nlevels_photo = 10. Contrary to previous trunk versions of ORCHIDEE where each canopy layer contained 0.5 units of LAI, in ORCHIDEE trunk, each canopy layer will have the same depth (in m) but will contain a different amount of LAI because tree canopies are assumed spherical.

Plant water stress

Describes r6614. With ORCHIDEE trunk there are two option to calculate waters stress:

• Same as in branch 3 and previous ORCHIDEE trunks, where a soil moisture stress functions limits C assimilation through constraints on the carboxylation capacity.
• The second possibility takes hydraulic architecture of plants into account, when calculating the plant water supply. This scheme, based on Hickler et al (2006), calculates the water supply as the ratio of the pressure difference between the soil and leaves. The plant water supply is the amount of water the plant can transport from the soil to the stomate accounting for resistance of water transport in the roots, sapwood and leaves. The resistances are inversely proportional to the conductivities in these different tissues, with the sapwood conductivity decreasing when cavitation occurs. If transpiration rates exceeds plant water supply, stomatal conductance is reduced. If you wish to make simulations with the hydraulic architecture choose the correct setting for energy_control (see section on Single layer vs. multi-layer energy budget for more elaboration on these flags).

The following PFT dependent parameters are needed for the calculations accounting for plant hydraulic architecture; minimal leaf water potential psi_leaf, sapwood leaf water potential that causes 50 % loss of xylem psi_50, maximum sapwood conductivity k_sap, root conductivity k_root, leaf conductivity k_leaf, specific root lenght srl, fine root radius r_froot, minimum root water potential psi_root. Moreover, as a further refinement to the hydraulic architecture, a more detailed description of the soil to root resistance has been included. This overcomes the need of psi_soil_tune in ORCHIDEE trunk as needed in ORCIDEE-CAN. Psi_soil is weighted by the fraction of evapotranspiration supplied by each soil layer, which depends on the soil to root resistance per layer that accounts for different root properties and hydraulic conductivity. In ORCHIDEE-CAN k_root described the conductivity of the fine roots and the soil. In ORCHIDEE trunk this joined conductivity has been split in a fine root conductivity and a soil to root conductivity. k_root is solely used in the calculation of the hydraulic architecture (because k soil to root is calculated separately in that subroutine), k_belowground is solely used in allocation. Hence, allocation and hydraulic architecture are no longer intimately linked.

The soil hydrology model distinguishes three water columns, the so called soil tiles. The first soil tile simulates the water contained under PFT1, the second the soil water contained under forest PFTs and the third the water contained under short vegetation PFTs (grasses and crops). PFTs are linked to soil tiles by means of the parameter pref_soil_veg (in src_parameters/constantes_mtc.f90). Note that both transpiration (veget) and evaporation (frac_bare_ns) is taken from the same soil tile.

Potential evaporation

Describes TAG4.1 (r7761).

Uses the method proposed by Milly (1992) to predict the potential evaporation at a specified temperature from a surface well-supplied with water based on the potential evaporation at a temperature which has been calculated using the actual surface evaporation. This latter method is what results from solving the energy budget, but the former allows for calculating the evaporation with a simple empirical function for the moisture availability that is well-constrained by observations. By using this correction, the model can use this moisture availiability function without recalculating the energy budget. The function in enerbil_flux has been updated with a few small bug fixes, but the correction does not appear to have been touched since ORC2 (or before).

Prescribe initial vegetation

Describes r6614. At the start of the model run or after a die-back or clear-cut new vegetation needs to be planted as ORCHIDEE does not grow vegetation from seeds. The initial dimensions of the vegetation are thus prescribed. Given that the allocation follows allometric relationships, any of the tree dimensions or any mass of any component could have been used to prescribe. The variable diameter was chosen because it is easy to (mentally) visualize the prescribed vegetation. In the run.def dia_init_min and dia_init_max need to be prescribed. Typical values are 1 to 3 cm. If more than one diameter class is used, dia_init_max needs to larger than dia_init_min. The larger the difference between the min and max, the more vegetation layers the canopy will be composed from. Grasses are prescribed by the parameters height_init_min and height_init_max because we have no diameter-height relationship for grasses. Note that height_init_min and height_init_max are the same. The fact that we have two values is a left over from previous versions of stomate_prescribe.f90 (see ticket #688).

In addition, the initial number of seedlings needs to be prescribed as well. For this a parameter needs to be set. This is a critical parameter to obtain acceptable model behaviour. If it is too high, lai saturates but the stand-level GPP will be distributed over too many individuals, each individual will grow very little and so it will take very long before self-thinning is reached. If it is set too low, LAI will be too low resulting in a too low GPP and thus very slow growth. A good starting values is a bit below self-thinning. That way the vegetation starts growing, individual are killed thanks to the background mortality and within 10 to 20 years self-thinning is reached. Why not starting at self-thinning? During code development it was tried to have the model start at the exact number of trees at which self-thinning will start given the diameter of the tree. One issues was that when prescribing small individuals (1 to 3 cm) the calculated number of trees could in the millions and so the GPP had to be distributed over too many individuals.

Depending on whether a stand is managed or not, the initial number of trees is calculated in a different way (see stomate_prescribe.f90). For managed stands (FM = 2, 3 and 4) alpha_rdi_upper and beta_rdi_upper are used to calculate a dynamically changing RDI. For unmanaged forest (FM = 1), grad_thin, alpha_rdi_upper, beta_rdi_upper are used to calculate RDI as long as the stand has not reached self-thinning. For better model performance the parameter grad_thin should be made PFT-dependent (see ticket #686)

Pseudo sugar loading

Describes r6614. The model code does not control the C/N ratio of the labile pool hence, even if there is a strong N-limitation, the model can accumulate lots of carbon in the labile pool. The first CN-version was indeed doing this as the plant could easily store several 1000 gC m-2. As this was considered unrealistic, the excess C in the labile pool was burned-off by some excess respiration. Although this luxury/wastage respiration has been suggested in the literature (see Amthor et al 2000 and Chamber et al 2004) it is not confirmed by many observations. It was therefore decided to control the size of the labile pool. The model already had an estimate of the optimal pool size of the labile and carbres pools. If the plant has more labile carbon than the optimal, GPP is downregulated because too much sugars in the leaves will increase the viscosity and hamper the sapflow in the phloem. The viscosity can be decreased again by closing the stomata and transpiring less of the sapflow in the xylem. By closing the stomata, GPP will be downregulated (Holtta et al 2017 doi:10.1093/treephys/tpx011). Because ORCHIDEE has no sapflow, turgor and viscosity yet, we used a simple ratio to downregulate NUE. Hence, the feature was called pseudo sugar loading to make clear this is not the real thing yet.

The regulation is smoothened by setting boundaries to avoid sudden decreases in GPP (which are not apparent in the data). Smoothing is taken care of in stomate_vmax.f90. If the plant has less carbon in its labile and carbres pools than wanted, the NUE is upregulated. Up regulation is also capped to avoid crazy NUE values and high frequency changes between up and downregulation. Up and downregulation are done in stomate_vmax.f90. The parameters than control the smoothing are sugar_load_min and sugar_load_max and can be set in the run.def.

Sugar loading can have a huge effect on the GPP even up to the point that realist LAIs come with way too low GPPs. This issues was eventually solved by setting the roots, leaves and sapwood longivity to observed values (for sapwood turnover this is not a very well constrainted parameter). Most important is to aim for a good NPP/GPP ratio. Interestingly when approaching the values reported in Campiolli et al 2015 and Vicca et al 2012, the reserves drop to more realistic values. At this point sugar loading is becoming much less dominant. From r7135 onwards, sugar loading cannot reduce the GPP by more than 25%. This implies that reserves could become crazy again and should be monitored regularly.

Recruitment

Describes r7784. When stands grow old the tree density decreases and under certain conditions the LAI can no longer cover the ground area. When this happens productivity will start to decrease. In nature the decrease in LAI comes with an increase in the amount of light reaching the forest floor which enables seedlings to grow and to restore the LAI. This process is known as recruitment. Note that in managed forest and forest with a lot of stand replacing disturbances (for example, fire) the forest may never reach the stage where the canopy sufficiently opens up to enable recruitment.

ORCHIDEE trunk can simulate recruitment for each PFT separately by setting recruitment_pft to true or false. When using age classes it makes sense to have the same setting for all age classes of the same species. The light-sensitivity of recruitment is determined by two parameters recruitment_alpha and recruitment_beta where recruitment_beta determines the shape of the relationship. When tuning these parameters one should take care that the number of recruits for 100% of the light reaching the forest floor should be similar to the number of recruits that is being prescribed for stand establishment. For the same recruitment_beta the model is very sensitive to the value of recruitment_alpha. When tuning this parameter change by steps of 0.1.

Recruitment has been developed, tested and validated for tropical forests. There is no reason why it shouldn't work for other forests. Initial test for temperate regions show that it works there as well. Also forest productivity at higher ages seems relatively sensitive to recruitment. If in unmanaged forests recruitment is high it will cause cyclic behaviour. The stand will reach the upper RDI, self-thinning occurs, the gap allows light to the forest floor, recruitment happens, trees grow, gaps are closed, the upper RDI is reached and the cycle starts all over again. The cycle itself is not a problem as it could be considered a simplified simulation of real gap dynamics. It becomes very regular when cycling over the same climate data and running without any other disturbances than self-thinning. Also recruitment is only a function of the light environment (not of the climate) which adds to the regularity. If, however, a too low recruitment is used, the frequency of the gap dynamics decreases and the biomass is getting very high.

At present recruitment was introduced at the PFT level. It probably makes more sense to link it to the management strategy than to the PFT. This needs to be checked (see ticket #737).

River routing

Describes r6614. This functionality works has not been evaluated yet for ORCHIDEE trunk, although it technically runs. Unless rivers are your main interest, set river_routing to n. The functionality still has problem when combined with the zoomed grids because pixels far away from the zoom become so big that some rivers can no longer reach the sea. Stand alone applications of ORCHIDEE trunk with river routing being activated require that the whole watershed in include in the simulation domain. If not the model will crash right at the start.

Root nitrogen uptake

Describes r7208. In previous versions of trunk 4, the calculation of root nitrogen uptake in nitrogen_dynamice had two parts. First, the maximum amount of root nitrogen uptake is calculated based on soil chemical properties and long-term root carbon mass. And then, nitrogen uptake is adjusted by the feedback from nitrogen reserve; if the long-term ratio of actual to potential nitrogen reserve pool is over one, then root nitrogen uptake decreases. The uptake does not decrease further than p_n_uptake times maximum root nitrogen uptake. The adjustment was added due to avoid carbon-limited growth (evident from too much nitrogen in the plants) which became more common after the code was changed such that the reserve nitrogen was used as much as possible. This setting is described in Nitrogen reserve pools), and it is controlled by the parameter use_reserve_n to y. A better use of the reserve nitrogen in the code was implemented after observing that the labile nitrogen pool was empty resulting in n-limited growth whereas the reserve pool was still well filled.

Later versions of trunk 4 suffered from too little n uptake rather than too much. The code limiting the n-uptake was removed because it was realized that the variable f_NCUptake had a very similar function. Enhanced plant_n_uptake was achieved by increasing the parameter vmax_uptake.

Root profile

Describes r7149. There are two different ways of looking at a root profile in the code. It could reflect "structure" or "function". Root structure is probably how most of us think about roots (i.e. digging a whole and observing where the roots are). When thinking at root structure the profile should be relatively constant over time. A logic time integrator to determine this constancy would be longevity_root as the profile cannot grow faster then the roots grow and die. In ORCHIDEE the structure profile is fixed it over time because the parameter that describes it (humcste) and the rooting depth are constant over time. Note that when ok_soil_carbon_discretization is used, the depth of the roots may depend on the thickness of the active layer. Likewise it would not be unlogic to make rooting depth a function of plant biomass as it could be expected that small plants need some time (and biomass) to develop a deep root system. When doing so the literature on this topic should be checked as in semi-arid regions plants probably propagate through vegetative reproduction because that would give them the opportunity to develop deep roots before becoming independent from the parent plant.

In r7149 a new subroutine calles hydrol_root_profile was added to hydrol.f90. This subroutine caldculates the rooting depth based on the PFT-dependent variable max_root_depth and the thickness of the active layer as calculated in stomate_soil_carbon_discretization. The nodes and interfaces of the soil layers follow the scheme proposed by De Rosnay (PhD thesis). This is calculated in vertical_soil.f90. The subroutine continues with calculating root structure as an exponentially decreasing function with depth which is the default in ORCHIDEE. This approach is based on root mass observations and, therefore, reflects a structural approach. The structural approach is used in the calculation of kfact_root which is water infiltration along roots (accounted for in hydrol.f90) and the input of soil carbon and nitrogen at depth due to the turnover of roots which is accounted for stomate_soil_carbon_discretization.f90. Furthermore, it is used to calculate the root temperature in stomate_resp.f90.

When thinking about root function it is not so important where the roots are located but it is more important at which depth the roots will be active. This seems a fair approach. Plants probably have way more roots than they need exactly to be able to quickly adjust to changes in the soil moisture environment. The function approach is used in the functional root profile where the plants can take most of the soil water from the layers where the soil water happens to be (feature added by D.Zhu in r4363. I'm not sure which versions this is in but we had it in trunk 4.0). The nodes and interfaces of the soil layers follow the scheme proposed by De Rosnay (PhD thesis). This is calculated in vertical_soil.f90. This way of looking at the roots is similar to how we look at the canopy where we have a lot of leaves at places in the canopy where little light can penetrate and where a large part of the photosynthesis is taken care of by the leaves in the top layers of the canopy. This functional approach is used in hydrol.f90 to calculate the water uptake from the soil for transpiration. It is also used in hydraulic architecture to calculate the water supply to the stomata and also the uptake of soil water from each soil layer.

The implementation of the functional root profile would benefit from a more careful literature study. Here are some references on the topic: doi.org/10.1002/eco.1523, doi.org/10.1007/BF00582230, doi.org/10.1029/2006WR005541, doi.org/10.1093/treephys/16.1-2.263 and doi.org/10.1073/pnas.1712381114. An interesting paper related to N and P uptake: doi.org/10.1126/science.aba0196

In ok_soil_carbon_discretization many more soil depths are introduced. When testing this functionality with trunk 4.0 some of the options could be removed and the use of different soil depths could be checked for consistency, i.e., are the same nodes and interfaces used as in the rest of the model? Is it really a new depth or just another name of the same depth?

Table The root profile interacts with the use of hydraulic architecture and soil discretization. Each of these combinations should be checked and tested.

ok_hydraul_arch	ok_soil_discretization	Comments
-	-	Default (r7135). Nodes and interfaces of the discretization have been checked. Extensively tested at the global scale.
-	+	Nodes and interfaces of the discretization have been checked. Tested for a limited test case.
+	-	Nodes and interfaces of the discretization have been checked. Tested for a limited test case.
+	+	Nodes and interfaces of the discretization have been checked. Tested for a limited test case.

Roughness

Describes r6790. The approached used to calculate z0m and z0h are identical to version 2.2 and 3.0. However, because r6790 simulates structured canopies and thus accounts for height growth, the differences might (needs to be confirmed) differ. In ORCHIDEE 2.2 and 3.0 the height is prescribed and fixed. So, at day 1 of the simulation a tropical forest with maybe 100 gC m-2 will already 27 m tall as there is no relationship between biomass and height. Given the PFT-specific values height values that were prescribed, a more correct name would have been "dominant height" rather than "height". In a forest the height could be 3 times lower than the dominant height because 80% of the trees could be in the understory and could be much smaller than the 20% dominant trees.

In r6790 height is a function of accumulated biomass and so it may take several decades before the simulated height of the dominant trees in a tropical forest would reach 27m. Furthermore, in r6790 the distinction is made between mean height (variable name: height) and dominant height (variable name height_dom). The use can choose whether they want to use "height" or "dominant height" in the calculation of roughness by setting the flag USE_HEIGHT_DOM. In r6790 the default value is .FALSE. which means that the model uses the exact same equations as ORCHIDEE2.2 and 3.0 but uses a poorly-defined height with regard to roughness. When set to .TRUE. the same equations are used but instead of height, the dominant height is used. Although the equations have changed this setting is conceptually closer to what was done in ORCHIDEE2.2.

run.def

Describes r6614. The format of the run.def has been modified to use the same structure for both coupled and off-line runs. The run.def is split between its driver settings (written in run.def), general items which control the ORCHIDEE model (in the orchidee.def file), and PFT-dependent parameters (in the orchidee_pft.def file). A script is now included with the config/ORCHIDEE_OL directory called MAKE_RUN_DEF. This scripts generates different orchidee_pfts.def files. Default model configurations (and thus different orchidee.def files) can be found in config/ORCHIDEE_OL. If one is interested in running two identical simulations with the land-only and the land-atmosphere model, only the run.def needs to be replaced. If one is interested in running the model in the same mode (land-only or land-atmosphere) but change the model functionality, only orchidee.def needs to be changed. If one is interested in running the same mode and same model functionality but only wants to change the PFTs tho orchidee_pft.def file will have to be changed. This modularity is expected to enhance the transparency of the model configuration.

Snags

Describes r8446. Snags (= standing dead trees that decompose slowly because there is no ground contact) were introduced as a separate component of the litter pool. This functionality can be activated by setting OK_SNAGS to .TRUE. (the default is .FALSE.). This functionality still needs to be parameterized. Once parameterized the flag could be removed and the model should always account for the snags. r8433 still uses hard coded values of 0.5 in stomate_litter.f90. These values should be externalized and parameterized as a function of the disturbance, for example, drought mortality and insect outbreaks result in 100% snags whereas uprooting following a wind storm results in 0% snags.

Soil maps

Describes r7352. ORCHIDEE trunk runs for two soil maps: zobler (3 soil types) and usda (12 soil types). Although the usda map has been the default for a while, ORCHIDEE trunk is back at using Zobler. The code has been adjusted in hydraulic_architect.f90 and stomate_windthrow.f90 so it works with Zobler and usda. The soil map can be selected in the orchidee.def by setting the value for soiltype_classif:

SOILTYPE_CLASSIF= zobler

The file names of the maps that need to be read can be set in the COMP/sechiba.card:

 (${R_IN}/SRF/SOIL/soils_param.nc, soils_param.nc)

If you want the background albedo to be calculated from the soil color, the model will find the information it needs it the ${R_IN}/SRF/SOIL/soils_param.nc file. If you want to use the background albedo based on MODIS set ALB_BG_MODIS to Y in the orchidee.def and specify the following in COMP/sechiba.card

 (${R_IN}/SRF/albedo/alb_bg_modisopt_2D_ESA_v2.nc, alb_bg.nc)

The usda map does not contain the soil colors so when using the albedo background map (thus alb_bg_modis = y) set:

SOILTYPE_CLASSIF= usda

in the orchidee.def and change the COMP/sechiba.card:

 (${R_IN}/SRF/SOIL/soils_param_usdatop4.nc, soils_param.nc),\
 (${R_IN}/SRF/albedo/alb_bg_modisopt_2D_ESA_v2.nc, alb_bg.nc)

when using usda without the albedo background map (thus alb_bg_modis = n) set:

SOILTYPE_CLASSIF= usda
SOILALB_FILE=soilcolor.nc

in the orchidee.def and change the COMP/sechiba.card:

 (${R_IN}/SRF/SOIL/soils_param_usdatop4.nc, soils_param.nc),\
 (${R_IN}/SRF/SOIL/soils_param.nc, soilcolor.nc)

The model will now read the soil color from the zobler map.

Single vs multi layer energy budget

Describes r6614. There are still some issues with the multi-layer energy budget, and it is currently only possible to run for one PFT. Thus, we suggest you use the single layer energy budget. In this can, however, be done by two different methods that in theory should give the exact same results: A - Use the energy scheme as found in the original enerbil.f90 B - Use the multi-layer energy scheme for a single canopy layer. This should reduce the multi-layer energy scheme to the original enerbil code but this is no longer the case in r6614. This needs to be restored (see ticket #276, #379, #457)

Several flags control which code is used to calculate the energy budget a functionality that is intimately linked to the way the plant water stress has to be calculated. For simplicity, a flag has been made to automatically control several of the flags related to the energy budget and the plant water stress, called energy_control.

Energy_control has 5 possible settings:

• Single layer, using the multi layer energy scheme (i.e. choice B above)
• Multi-layer energy budget
• User specific as set in the run.def
• Single layer, using the original enerbil scheme (i.e. choice A above)

1 - DEFAULT uses the enerbil module in combination with the hydraulic architecture (ok_hydrol_arch and ok_gs_feedback true, while ok_mleb, ok_impose_canopy_structure and ok_bare_soil_new are set to false). 2 - option to use enerbil module and original water stress (not hydraulic architecture). 3 - The energy budget is calculated using the multi-layer energy scheme with a 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 is 1,jnlvls_under is 0,jnlvls_canopy is 1,jnlvls_over is 0). No mleb output, ok_mleb_history_file and ok_bare_soil_new are set to FALSE. 4 - multi-layer energy budget: ok_hydrol_arch, ok_gs_feedback and ok_mleb all TRUE. ok_impose_canopy_structure is False, and the energy budget is calculated for multiple layers (jnlvls is 29,jnlvls_under is 10,jnlvls_canopy is 10,jnlvls_over is 9). No mleb output, ok_mleb_history_file is set to FALSE. ok_bare_soil_new is set to TRUE 5 - user specific: user specific settings for these controls and layers as defined in the run.def by the user.

The default setting of energy_control in ORCHIDEE trunk is 1. This is well tested and considered safe to use. More details of the flags controlling the energy budget and water stress calculations can be found in src_parameters/control.f90. The flags controlled by the energy_control are:

• ok_hydrol_arch: Flag that activates the hydraulic architecture routine (true/false). Previous trunk versions of ORCHIDEE used soil water as a proxy for water stress and applied 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 controlled by a separate flag: ok_gs_feedback.
• 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.
• 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 9 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.
• 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.
• 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 slows down the code.
• ok_bare_soil_new: This flag is described in a separate ​section.

Specific leaf area

Describes r6614. Similar to branch 3.0, ORCHIDEE trunk users can choose to use a constant specific leaf area (SLA) or a dynamic (SLA) by setting the flag dyn_sla. SLA is a fundamental parameter in the allocation scheme and it is well established that SLA is dynamic in nature especially during leaf onset. Nevertheless, the effect of this flag is much more limited than one could expect based on the perceived importance of sla.

Spin-up

Describes r7178. The analytical spin-up works with ORCHIDEE trunk. To ensure growth at the onset of the analytic spin-up for all PFTs initial values are needed for the SOM pools, and that is irrespective of whether IMPOSE_CN is y or n. The initial pools are set in stomate_io.f90, and they can be specified in the run.def by:

SOM_INIT_ACTIVE = 1000
SOM_INIT_SLOW = 3000
SOM_INIT_PASSIVE = 3000
SOM_INIT_SURFACE = 1000

The initial values of carbon for the four SOM pools are used globally. Sensitivity test have shown that the analytic spin-up is not sensitive to the actual size of the initial values. The different ecosystems will after a while settle into their own states (N limited or saturated) in spite of having the same initial state. Convergence appears to be faster if initial values are set.

Whether the analytical spinup will be used or not is set in COMP/stomate.card by the parameter SPINUP_ANALYTIC. By default the SPINUP_PERIOD is calculated as the length of the time series used in cyclic_begin and cyclic_end as set in the config.card. This calculation is done in the COMP/stomate.driver

### Begin Code in COMP/stomate.driver                                           
    if [ X${stomate_UserChoices_SPINUP_ANALYTIC} = Xy ] ; then                  
       IGCM_comp_modifyDefFile nonblocker orchidee.def SPINUP_ANALYTIC y        
                                                                                
       # Test if CyclicBegin and CyclicEnd is set in config.card                
       if ( [ X${config_UserChoices_CyclicBegin} = X ] || [ X${config_UserChoices_CyclicEnd} = X ] ) ; then
           IGCM_debug_Exit "CyclicBegin and CyclicEnd must be set in config.card to run option spinup_analytic."
       fi                                                                       
       # Calculate and set number of years of forcing data                      
       CycleNb=$(( ${config_UserChoices_CyclicEnd} - ${config_UserChoices_CyclicBegin} + 1 ))
       IGCM_comp_modifyDefFile nonblocker orchidee.def SPINUP_PERIOD ${CycleNb}
    else                                                                        
        IGCM_comp_modifyDefFile nonblocker orchidee.def SPINUP_ANALYTIC n       
        IGCM_comp_modifyDefFile nonblocker orchidee.def SPINUP_PERIOD -1        
    fi                                                                                                                                       
### End Code in COMP/stomate.driver  

The calculated value can be overwritten by setting SPINUP_PERIOD in PARAM/orchidee.def. In rXXXX a pre-spinup was introduced as SPINUP_ANALYTIC_FG0 and SPINUP_ANALYTIC_FG0nd. The purpose of the pre-spinup is to break the synchronicity of neighboring pixels. With FG0, all pixels will be zero years old at t=0 of the FG1 spinup. Because neighboring pixels often experience the same climate and have the same soil properties they will all mature and die within the same decade. This results in some cyclic NBP for some PFTs. FG0 is run for 200 years and within this 200 years each pixel x PFT is cut once at random (the random cutting is prescribed through a map and is therefore reproducible). This random cutting breaks the synchronicity of neighboring pixels and is thought to result quicker in a more stable steady-state (to be confimed).

Start and restart

Describes r6614. The table below shows typical configurations to start and restart ORCHIDEE trunk

≈≈ -

Sechiba				stomate
Restart	impveg	impsoil	laimap	ok_stomate	restart	Comment
+	-	1. Typical for a spin-up from scratch
-	-	-	+	-	n.a.	2. Typical for a spin-up with sechiba only
-	+	+	-	+	-	3. Typical for a site-level spin-up
-	+	-	-	+	-	4a. Typical for a test run with prescribed PFTs
-	-	+	-	+	-	4b. Typical for a test run with prescribed soil properties
+	-	-	-	+	+	Typical for a simulation following the spinup described in 1
+	-	-	+	-	n.a.	Typical for a simulation following the spin-up described in 2
+	+	+	-	+	+	Typical for a simulation that will be compared against site-level observations and that was spun-up with 3
+	+	-	-	+	+	Typical for restarting a test run following a test run described in 4a
+	-	+	-	+	+	Typical for restarting a test run following a test run described in 4a

The code allows for many more combinations than described above. Some of those configuration could be useful in specific cases but be aware of the following:

• If the restart file for sechiba does not come from the same simulation/year as the restart file for stomate, the values for the variables that are stored in both sechiba and stomate will be overwritten by the values stored in the stomate restart file because that file is read later in the initialization phase. Some variables are duplicated because that adds the flexibility to the configuration that can be used to restart the model.
• When read_lai is not -1, ok_stomate should be false (but this is not yet enforced by the model - see ticket #685). It may work but it is flawed from a conceptual point of view. read_lai was developed for cases where the model cannot simulate its own canopy. ok_stomate was developed so the model would simulate its own canopy.
• Combining a sechiba restart from with impveg, impsoil or laimap may result in values from the restart file being overwritten because impveg, impsoil and laimap are checked after reading the sechiba restart file.
• NCO commands could be used to hack into the restart file. This could be powerful when stomate is used and one wants to restart the model after a spin-up but also want that the forest starts growing the first year of the simulation. In that case circ_class_biomass and circ_class_n should be set to zero. This approach is straightforward when no age classes are used. When age classes are used many more variables need to be moved to the correct PFT to avoid conflicts.

Stomate bypass

Describes r7830. For debug or test purposes at site level, a new subroutine has been implemented to bypass the stomate calculation. The integer read_lai controls this routine. If read_lai is set to 1, the model skip the calculation of stomate module and enters the subroutine slowproc_impose_site_level_lai. In this subroutine, the vegetation variables used in sechiba are either imposed via the run.def or calculated using the same equations as in stomate (direct calculation from the imposed variables instead of running the whole stomate module). The routine has several ways of functioning according to the parameters set in the flag. The user can either impose a seasonal dynamic of the LAI by prescribing the parameters maximum_lai, minimum_lai, phenology_doy_start, senescence_doy_end, phenology_delta and senescence_delta (from this, the biomasses and the canopy structure are directly calculated) or impose the full canopy structure via the variable lvllai_prescribe which is directly given to lai_per_level. Finally, the user can impose the heights of the canopy thanks to the variables min_canopy_height_prescribed and max_canopy_height_prescribed.

If read_lai is set to 2, the user can impose a LAI map. See also the item on "Leaf area index map"

Sechiba vs stomate

Describes r6614. In ORCHIDEE trunk the links between sechiba and stomate have been strengthened compared to previous ORCHIDEE trunks. As in previous versions, stomate makes use of variables calculated in sechiba but in ORCHIDEE trunk, sechiba requires information from stomate to work properly. It is possible to prescribe the canopy structure and thus only run sechiba. Set read_lai to 1 or 2, then, according to the way of functioning, the data for a canopy structure is either given thanks to a map from an ORCHIDEE trunk run with stomate (all the required information is stored in the sechiba restart file to enable restarting the model without stomate) or imposed through the run.def for site level tests.

Tree ring width

Describes r6721. Tree ring width is always calculated and is not controlled by any flags or additional parameters. The basal area increment is central in the allocation scheme when making use of the Deleuze and Dhote rule to account for within stand competition between trees with different dimensions. It is saved to the output variable CCDELTABA (mm2 day-1). Basal area increment can be used to calculate the tree ring width when assuming that the cross-section of a tree trunk is circular. The output variable CCTRW represents the radial increase of individual trees for each circumference class (mm day-1). Since CCTRW is simply sqrt(ba_later/pi) - sqrt(ba_beginning/pi), it reconstructs trees from zero at year 1 removing the effect of BA_init. An implication of this approach is that the value at the first time step does not only include the growth calculated in stomate_growth_fun_all but also the basal area that was prescribed in stomate_prescribe.f90. The contribution of the initial tree dimensions prescribed in stomate_prescribe becomes more important if bigger trees are prescribed. Both CCDELTABA and CCTRW are written to stomate_4dim.nc.

Trusting

Describes r8058. The flag ANNUAL_PROC_END_OF_YEAR with default True is added. If it is set to false in orchidee.def, then the annual processes are done in the end of June. As long as the 1Y+1Y=2Y test cannot be performed in the trusting, we have no means to check whether the processes done at the end of year are restartable. One work around is to add the functionality to do the annual processes in the end of June instead of the end of the year. This function should only be used in the trusting it has no operational applications. This option could later be added to the default trusting tests.

Describes rXXXX. When set to y the flag TRUSTING_HACK_AGE_CLASS skips the calculation of redistributing biomass between age classes. In age_class_distr their is a calculation that redistributes the biomass between age classes. Although mass is conserved during this calculation it results in precision errors that propagate through the model calculations. In the trusting we could have one test case that helps to check the technical integrity of the age class distribution code: if the model is run without LCC and disturbances, A simulation without age classes should generate the exact same results as one with age classes (but the dimensions of these simulations differ). If biomass is recalculated this test fails. With this hack, the biomass is explicitly moved from one age class to another without any calculations and thus without any precision errors. With this hack set to .TRUE. the model can be trusted for the age class code. This hack should not be used in simulation experiments as biomass from different age classes will NOT be merged but simply overwritten. If all goes well this should result in a mass balance error. This option could later be added to the default trusting tests.

Vertical discretization soil carbon

Describes r6721. To be completed.

Wind throw

Describes r6614. The calculation of wind storm damage can be activated by setting ok_windthrow to y in the orchidee.def. This module recognizes storms and calculates the critical wind speed of a stand taking stand and soil properties into account. If the stand is managed, the damaged trees are salvaged logged. If the stand is unmanaged the damaged trees are left on-site to decompose. The default setting for ok_windthrow is n. The damaged fractions (if this fraction is more than the replace_threshold of the total biomass per PFT) of the stands are replanted and moved to the first age class (if more than 1 age class is used).

Winthrow damage calculation can be divided into two parts: storm recognition and damage calculation.

Storm recognition is sensitive to 3 parameters

wind _ratio_threshold	The ratio of the current wind speed to the 5-year-mean wind speed will be summed if the ratio is above this threshold. A lower threshold causes more frequent storms. Note that excessively low values can cause continuous storms in the coastal areas.
wind_sum_threshold	The first condition to identify a storm. If the total sum of that ratio exceeds this threshold, then it is considered a storm. A lower threshold will result in the storm affecting a larger area.
wind_speed_storm_thr	The second condition to identify a storm. Maximum wind speed of 3 days (legacy_years_wind) should be over this threshold. Lower threshold causes more frequent storm.

Calculation for windthrow damage is sensitive to 2 parameters

daily _max_tune

This represents wind gustiness. ORCHIDEE uses interpolated 6-hourly wind speed data, which means that it cannot capture wind gusts that occur at finer spatial and temporal scales. A higher value leads to greater damage.

factor_closer and sfactor_further (PFT-level parameters)

This parameter alters the shape of the sigmoid function which represents the damage probability within a pixel. A low value makes the function stepwise, which causes near-zero damage when the wind speed is lower than the critical wind speed, and near the maximum damage when the wind speed is higher than the critical wind speed. A high value smoothes the function, which allows for small gaps from wind damage and damage in a broader area.

One setting that needs caution is use_fluxnet. This flag was created for the case in which the model utilizes fluxnet forcing but it can also be applied to coupled simulations. If use_fluxnet is set to y, 1) the wind speed is not retrieved (Check stomate.f90). Because loaded wind speeds from 6-hourly forcing become lower before interpolcation to translocate the wind (10m -> 2m above zhd) it was retrieved before storm recognition. 2) Gustiness is not considered.

XIOS

Describes r7994. When using output_level_stomate_history=10 the model may crash. This is because we are trying to write too many data to XIOS. There is a limitation related to HDF5 version that causes an XIOS error (terminate called after throwing an instance of 'xios::CNetCdfException' what(): Error when calling function nc_enddef(ncId), NetCDF: HDF error). This kind of error is due to a combination of the number of output variables, the number of metadata variables AND the parallel writing. There is no solution except to reduce the number of variables or writing in sequential mode (that means with only 1 XIOS server or in "multiple file" mode rather than "one file" mode or by activating 2nd level server). We hope that will be solved in next releases of HDF5. The solution is to reduce the output_level_stomate_history or to manually select the variables you really need (src_xml/file_def_orchidee.xml - no need to recompile).