ORCHIDEE-CN-CAN

ORCHIDEE-CN-CAN is a combination of three code bases: the trunk as it is in 2019, the ORCHIDEE-CN branch in 2016 and some key functionality of ORCHIDEE-CAN (aka ORCHIDEE-DOFOCO on the svn server). The committed version can be considered a light version of ORCHIDEE because some of the available functionality has not yet committed (for example, the DGVM). Nevertheless, this light version is believed to contain all functionality that affects the parameterization of the core of the model. Consider the model as two parts: 1) core processes that effect fluxes and pools on the per square meter level, and 2) perirpheric processes that involve everything else. In this picture, the functionality that is not yet included depends on core processes like vegetation growth, soil hydrology or the energy budget, but it does not affect these processes at the pixel level.

• The available functionality is described in more detail below
• Missing functionalities compared to the trunk: DGVM, IPCC woodharvest, prescribed LAI, and fire disturbances
• Missing functionalities compared to ORCHIDEE-CN: none
• Missing functionalities compared to ORCHIDEE-CAN: none
• Better functionalities than those available in the trunk, ORCHIDEE-CN or ORCHIDEE-CN-CAN: other branches have used an improved fire module (Spitfire) for fire disturbances and gross land cover changes, and we could consider adding those to Tag 4.0
• Coupling: LMDZ coupled with ORCHIDEE-CAN has been used only with nudging and zoomed over Europe. Given that several key parameters for the coupling have changed: albedo, roughness is now dynamic, more landscape heterogeneity when using age classes, transpiration, mass balance and quality control functions … it is expected that lots of testing will be required to come to a satisfying result with LMDzOR-CN-CAN.

What is the future of the version?

Now that the AR6v1 is stable, a robust well tested version of ORCHIDEE has been tagged. This gives the ORCHIDEE-team the opportunity to prepare for the future by making a series of far reaching changes to the code. The first ambition is to add CN into AR6v0 to obtain Tag 3.0 before the end of 2019. In parallel, ORCHIDEE-CN-CAN has been merged with the TRUNK during the year of 2019 to prepare for Tag 4.0. We will continue testing and discussing Tag 4.0 before making this the trunk for both coupled and off-line simulations. The basis of Tag 4.0 is the current ORCHIDEE trunk (the version to become Tag 3.0) into which CAN functionality is merged. At this stage key functionality from other branches may be merged as well (to be discussed in tandem with the fair use policy). In addition to the recent improvements to the trunk, the new trunk (Tag 4.0) should have an increased functionality (see below), increased internal consistency (see below, for example albedo and photosynthesis) and some long awaited fixes (see below, for example tot_bare_soil).

How can I contribute to this version?

Given the large number of code changes, there is no guarantee that the current version of ORCHIDEE-CN-CAN (future Tag 4.0 of the TRUNK) will run flawless for all set-ups and possible combinations of flags currently being used by ORCHIDEE users and developers.

• All users and developers are invited to install future Tag 4.0 and set-up their favourite simulation. Use the run.def in config/ORCHIDEE_OL/OOL_SEC_STO_FG2/ as a starting point, or play with the Python script in config/ORCHIDEE_OL/MAKE_RUN_DEF directory for your favourite set-up (see below). Report on the outcome of your tests during the ORCHIDEE meetings.
• All developers are invited to contribute bug fixes, and code for improved functionality to Tag 4.0.
• Adjust the interface such that Tag 4.0 can be coupled to the latest version of LMDz and can be used for regional simulations with WRF.

How do I run this version?

The instructions to download and extract ORCHIDEE-CN-CAN are found on the wiki (see ​https://forge.ipsl.jussieu.fr/orchidee/wiki/Documentation/UserGuide/ORCHIDEEDOFOCOInstall). When you check-out the model, you also check-out OCRHIDEE_OL, where OOL_SEC_STO_FG1trans, OOL_SEC_STO_FG2 and SPINUP_ANALYTIC_FG1 have been updated to contain the configurations needed to run ORCHIDEE-CN-CAN. The configurations are based on the ORCHIDEE_OL from the trunk as of December 2018. In ORCHIDEE_OL there is a MAKE_RUN_DEF folder. This folder contains the parameter values and a script that will make all run.defs CN-CAN users have requested this far (i.e. 13, 15, 28 and 37 PFTs both for the FG2 and ENSEMBLE configuration). If you would like to use another run.def let us know and we can help you to adjust the script.

When running the model, the default setting should be XIOS, thus compile the model with gmake.

If you chose to use your own configurations to run ORHICDEE-CN-CAN, be aware that new output files have been added to ORCHIDEE-CN-CAN, and this has to be specified in the cards and drivers.

Below is given an example for stomate_history_4dim.
In stomate.driver, add the following within the shown if statement:

    # Modify file_def_orchidee.def if XIOS is activated
    if [ X${orchidee_ol_UserChoices_XIOS} = Xy ] ; then
        # Modify file_def_orchidee.xml
        IGCM_comp_modifyXmlFile nonblocker file_def_orchidee.xml stomate3 enabled ${stomate_enabled}
        IGCM_comp_modifyXmlFile nonblocker file_def_orchidee.xml stomate3 output_freq ${stomate_freq}
    fi

In stomate.card add the new file to list of output files:

[OutputFiles]
List=   (stomate_history.nc, ${R_OUT_SBG_O_M}/${PREFIX}_1M_stomate_history.nc, Post_1M_stomate_history) \
        (stomate_ipcc_history.nc, ${R_OUT_SBG_O_M}/${PREFIX}_1M_stomate_ipcc_history.nc, Post_1M_stomate_ipcc_history),\
        (stomate_history_4dim.nc, ${R_OUT_SBG_O_M}$ {PREFIX}_1M_stomate_history_4dim.nc, NONE)

For sechiba, the following changes are needed in sechiba.driver:

    if [ X${orchidee_ol_UserChoices_XIOS} = Xy ] ; then
      # Modify file_def_orchidee.xml file
      IGCM_comp_modifyXmlFile nonblocker file_def_orchidee.xml sechiba1 enabled ${sechiba1_enabled}
      IGCM_comp_modifyXmlFile nonblocker file_def_orchidee.xml sechiba1 output_freq ${sechiba1_freq}
      IGCM_comp_modifyXmlFile nonblocker file_def_orchidee.xml sechiba3 enabled ${sechiba1_enabled}
      IGCM_comp_modifyXmlFile nonblocker file_def_orchidee.xml sechiba3 output_freq ${sechiba1_freq}
      IGCM_comp_modifyXmlFile nonblocker file_def_orchidee.xml sechiba2 enabled ${sechiba2_enabled}
      IGCM_comp_modifyXmlFile nonblocker file_def_orchidee.xml sechiba2 output_freq ${sechiba1_freq}
    fi

If you are seeing the following error when you run:

Error [StdIStream& operator>>(StdIStream& in , CDuration& duration)] : In file '/home/users/mmcgrath/ORCHIDEE-CN-CAN/modeles/XIOS/src/duration.cpp', line 74 -> Bad duration format: impossible to read a pair (value, unit).

then it is possible you have not added the required lines in stomate.card, stomate.driver, sechiba.card, and sechiba.driver. This error arises because the .xml files, when they have been copied to the working directory, cannot have values like "_AUTO_" in them. libIGCM is supposed to replace such values, which exist in the template .xml files in modeles/ORCHIDEE/src_xml/, when it copies files from the run directory to the working directory. If there is something missing in the .driver files, this replacement doesn't happen and XIOS throws the parsing error above.

If the models keeps recompiling XIOS (which has been known to happen on obelix) the following may help. To avoid this and reduce compiling time, change the following line in the Makefile from

with_xios : xios ioipsl driver_xios verif

to

with_xios : ioipsl driver_xios verif

Functionalities (alphabetical order)

Age classes

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 when an older age class is killed.

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-CAN, 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 run.def. A Python-script was written to create this kind of run.def. Increasing the number of PFTs has important consequences for the speed of the model and the memory use, although ORCHIDEE-CAN 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. 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), 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 attached run.def for a functional example. 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, C3 grass, C4 grass, C3 crop and C4 crop and 4 times 8 for the 8 forest PFTs. Thus NVM = 37
• 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=13
• 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 (general)

ORCHIDEE-CN-CAN 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 avaialble in forcing data. In ORCHIDEE-CN-CAN 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 expontential 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-CN-CAN 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 Albedo (background)).

Albedo (background)

If covered by snow, the background albedo is calculated by the snow module and accounts for snow age and snow density (needs to be checked – last time snow did not account for NIR). If not covered by snow, the background albedo is not simulated but prescribed by the parameters bgrd_ref_vis and bgrd_ref_nir. 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, reflecting 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-CN-CAN has two options to prescribe the background albedo:

• The background albedo is prescribed per PFT but is constant throughout the year. This is the option that has been used in ORCHIDEE-CAN and is the option that has been validated 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 validated yet. Set alb_bg_modis = y.

Albedo (snow)

The snow albedo could be either prescribed (in condveg_init.f90) or calculated following Chalita and Treut (1994) do_new_snow_albedo = n or calculated following CLM3 do_new_snow_albedo = y. The difference between the latter two methods has not been tested yet. The CLM method was added to CN-CAN, the Chalita and Treut method was added in parallel to the runk. When merging both versions we ended up with two options.

Allocation

ORCHIDEE-CN-CAN uses the allometric allocation as developed in OCN. In ORCHIDEE-CAN 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, the allocation scheme is identical between all 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. 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 too much. This dampening was observed in ORCHIDEE-CN and ORCHIDEE-CN-CAN. IN ORCHIDEE-CNP this temperature relationship was removed because the interannual variability became unrealistic.

ORCHIDEE-CN-CAN 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.

There are no options to revert to the allocation based on resource limitation. All references and parameters for allocation based on resource limitation have been removed from the code (those that were overlloked can be removed). Allometric allocation makes use of the following PFT-specific parameters: sla, tau_root, tau_leaf, tau_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_root, and k_leaf (see also plant water stress). As such there is a functional link between C and N-allocation and the hydraulic architecture of a plant. Details on the parameters can be found in the SI of Naudts et al 2015 in GMD or in src_parameters/constantes_mtc.f90.

Anthropogenic species change

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-CAN. 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. 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

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 as a function of veget. When a deciduous PFT sheds its leaves, the gaps in the forest will contribute to 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 deserts as the background albedo of forest gaps. 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. Moving gaps to the bare soil is no longer needed. So, if ok_bare_soil_new is set to TRUE, canopy gaps no longer contribute to the bare soil. It needs to be tested what will happen with the evaporation in the single-layer model. The multi-layer energy budget should be able to correctly deal with the gaps in the canopy because the diffusivity will increase when the canopy is becoming sparser.

Bark beetles

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 wich 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.

C13

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.

Configuration

The model comes with the following configuration set-ups (config/ORCHIDEE_OL):

• SPINUP: under development.
• ENSEMBLE: under development.
• SPINUP_ANALYTIC_FG1: 2x2 degrees CRU-NCEP forcing cycles over 1901-1910 between 1901 and 2240. Start from scratch. 15 PFTs, no land cover changes, nitrogen input for 1850, 1860 or 1900 depending on the species, and fixed CO2 concentrations.
• OOL_SEC_STO_FG1trans: 2x2 degrees CRU-NCEP forcing cycles over 1901-1910 between 1861 and 1900. Restart from SPINUP_ANALYTIC_FG1. 15 PFTs, no land cover changes, nitrogen input for 1900, and annual CO2 concentrations.
• OOL_SEC_STO_FG2: 2x2 degrees annual CRU-NCEP forcing between 1901 and 2010. Restart from FG1trans. 15 PFTs, annual land cover changes, annual nitrogen input and annual CO2 concentrations.
• OOL_SEC_STO_FG3: 0.5x0.5 degrees annual WFDEI_GPCC forcing between 1979-2009. Start from scratch. 15 PFTs, annual land cover changes, annual nitrogen input and annual CO2 concentrations.
• OOL_SEC_STO_FG3nd: 0.5x0.5 degrees annual WFDEI_GPCC forcing between 1979-2013. Start from scratch. 15 PFTs, annual land cover changes, annual nitrogen input and annual CO2 concentrations. This configuration has been set-up and committed but has not yet been tested. Waiting for the new driver to be tested in CN-CAN.
• OOL_SEC_STO_FG4: ~0.5x0.5 degrees annual CRU-NCEP forcing between 1901 and 2010. 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: 1x1 degrees annual IPSL RCP 4.5 forcing between 1911 and 2100. Start from OOL_SEC_STO_FG4. XX 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. This configuration is under development. Waiting for the boundary files to be copied to orchideeshare. The species and management change functionality needs to be tested within CN-CAN.
• OOL_SEC: 2x2 degrees annual CRU-NCEP forcing between 1901 and 2010. 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 disactivated. This configuration is under development, waiting for a global LAI map (in case of CN-CAN this has become biomass map) to be tested.
• TESTSTOMATE: Not maintained for ORCHIDEE-CN-CAN. Will be added back once teststomate.exe has been restored.
• FORCESOIL: Not maintained for ORCHIDEE-CN-CAN. Will be added back once forcesoil.exe has been restored.

Consistency checks

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-CN-CAN 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.

CWRR vs Choinel

ORCHIDEE-CN-CAN was developed and tested with CWRR. The Choinel code is still available but was never used in ORCHIDEE-CN-CAN and can thus be removed. The hydrological schemes were not touched during the development of ORCHIDEE-CN-CAN. The hydraulic architecture in ORCHIDEE-CN-CAN relies on CWRR.

Diameter classes

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-CAN and ORCHIDEE-CN-CAN are coded and tested for NCIRC = 3 but the model has been run with one diameter class as well +++REDO. Until further testing, three diameter classes are considered sufficient.

Given earlier choices in ORCHIDEE, 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.

Forced clear cut

OK_CLEARCUT is a flag used to force ORCHIDEE-CN-CAN to clearcut a forest after one year of simulation. This flag is set to TRUE in order to restart a new stand at the begining 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).

Forest management and management changes

70% of the global forests are managed, which contradicts the assumption in previous versions of ORCHIDEE that forests are long-lived natural vegetation. 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, only the principles of ORCHIDEE-FM, i.e., the allocation rule of Deleuze and Dhote that allocates carbon to different diameter classes based on the basal area of the tree, and relative-density index (RDI)-based management which enforces thinning and harvest operations based on the current tree density and the self-thinning density, were retained.

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. 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. Therefore, the biomass pools in ORCHIDEE-CN-CAN no longer depend on a prescribed longevity.

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-CN-CAN 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. Unless, one wants to run forest management over Europe the user will have to create his/her forest management maps first. Set read_fm_map to y and specify the location of the forest management map in COMP/stomate.card. Check the existing forest management maps for Europe for an example of how the map should be defined.

+++UPDATE+++ Dynamic RDI was implemented and prescribe has been changed When prescribing a forest stand (independent of forest management) the Initial density nmaxtrees, and the range of the initial tree height of the seedlings needs to be specified height_init_min and height_init_max. Irrespective of the management strategy the maximum carrying capacity needs to be described. Carrying capacity was formalized through the self-thinning relationship which makes use of two parameters alpha_self_thinning 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 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 : first thinning height h_first, stand replacing density ntrees_dia_profit, harvest diameter max_harvest_dia, coppice diameter coppice_diameter, rotation length src_rot_length, number of rotations src_nrots, fuelwood diameter fuelwood_diameter and the minimum and maximum alpha and beta (thus 4 parameters) specifying the RDI range alpha_rdi_upper, alpha_rdi_lower, beta_rdi_upper and beta_rdi_lower.

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-CAN and ORCHIDEE-CN-CAN. 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.

While developing the code some conflicts were encountered between RDI and self-thinning. As a first solution an additional threshold was introduced rdi_limit_upper. When debugging progressed this threshold was set to 0.99 (if set to 1.00 there is no correction any more). The initial problem was resolved but the initial fix has not been removed yet. For the time being set rdi_limit_upper to 0.99. Persons interested in simplifying the code could start with this removing this parameter.

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-CAN. 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.

Grasslands

In the DOFOCO branch, grasslands were treated as evergreens as dead grasslands in winter at the time had too high of albedo compared to real grasslands. After the work of Nicolas Vuichard, we have moved CAN back to deciduous phenology, trying to follow what is done in the TRUNK as closely as possible. For the moment we are 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.

Harvest

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!). 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 should be replaced by 2, 17 and 50 years pools in Europe.

Leaf area index map

Four flags have been identified that control the model behavior in terms of lai: ok_stomate, ok_pheno, impose_veg 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 the ORCHIDEE and ORCHIDEE-CN, canopy structure is prescribed by a single variable lai and the assumption of a turbid medium (Lambert-Beer). Consequently reading an lai value suffice to prescribe the entire canopy. In ORCHIDEE-CN-CAN, 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-CN-CAN would require a substantial number of assumptions to turn an aggregated 1D lai value into a disaggregated 3D canopy at the PFT level.

Land cover change (with age classes)

Land cover change now accounts for age classes. It is controlled by the flags land_cover_change and veget_update. Set land_cover_change = n and 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 burried 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.

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 farctions 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 critrion 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.

Litter decomposition

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

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-CN-CAN 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.

Mortality

ORCHIDEE-CN-CAN distinguished 3 types of natural mortality. The first two options are similar to those in previous version of ORCHIDEE and are set by the flag constant_mortality. If constant_mortality = y, the background mortality of a forests is calculated as a constant, prescribed fraction. In ORCHIDEE-CN-CAN, this fraction is given by residence_time (see also forest management). If constant_mortality = n, the 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, tested and can be used in ORCHIDEE-CN-CAN. However, because of the introduction of self-thinning (the third type of natural mortality) in ORCHIDEE-CN-CAN, constant_mortality = y soon became the default setting. In ORCHIDEE-CN-CAN, 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. 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. It needs to be tested how the interplay between background mortality and self-thinning will work out when constant_mortality = n is used.

Notice 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

ORCHIDEE-CN-CAN strictly follows ORCHIDEE-CN 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

Parameter files

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 toghether 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 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 expecation 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.

Phenology (forced)

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 particulary welcome in coupled simulations but likley 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 forced phenology is used, 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, CN-CAN is all about vegetation dynamics and thus instabilities in the vegetation cover. In CN-CAN there are two processes that can deal with dying PFts including evergreens PFTs. First, ok_recruitment could used. If ok_recruitment = .TRUE. a decreas 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. Generation can take over from each other without loosing the canopy cover entirely. Second, if there are insufficient reserves to grow no leaves, there will be no or insufficient gpp, the crabon 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. CN-CAN was developed to work with always_init = .FALSE. so this has become the default value, contrary to the trunk where always_init = .TRUE. is the default.

Photosynthesis

Photosynthesis is calculated as in ORCHIDEE-CAN, and ORCHIDEE-CN 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-CAN used nlevels_photo = 10. This values was retained in ORCHIDEE-CN-CAN. Contrary to ORCHIDEE-CN and previous versions of ORCHIDEE where each canopy layer contained 0.5 units of LAI, in ORCHIDEE-CN-CAN, 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

With ORCHIDEE-CN-CAN there is two option to calculate waters stress.

• Same as in the trunk, 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.

Moreover, as a further refinement to the hydraulic architecture, a more detailed description of the soil to root resistance has been included (now the default setting) and is activated with the flag ok_soil_root. This overcomes the need of psi_soil_tune in ORCHIDEE-CN-CAN as needed in ORCIDEE-CAN. When ok_soil_root = y, an adjustment of the soil water potential (psi_soil_root) due to soil to root resistance is allowed. 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.

If you wish to make simulations with the hydraulic architecture the CWRR hydrology scheme needs to be activated (hydrol_cwrr=y). Moreover, both mleb and ok_hydrol_arch needs to be true. These can be controlled by the overall flag multi_layer_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.

Prescribe initial vegetation

+++UPDATE+++ 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 is 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 height was chosen because it is easy to (mentally) visualize the prescribed vegetation. In the run.def height_init_min and height_init_max need to be prescribed. Typical values are 2 to 5 meter. If more than one diameter class is used, height_init_max needs to larger than height_init_min. The larger the difference between the min and max, the more vegetation layers the canopy will be composed from.

In addition, the initial number of seedlings needs to be prescribed as well. For this the parameter nmaxtrees needs to be set. nmaxtrees 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 (2-3 m) the calculated number of trees could in the millions and so the GPP had to be distributed over too many individuals.

The way the initial vegetation is prescribed can also be used to prescribe a mature vegetation right at the start of the simulation. One could set height_init_min and height_init_max to for example 15 and 20. If nmaxtrees is not adjusted, massive self-thinning will take place on the first day. It is better to set nmaxtrees to a value just below the exact self-thinning value. Prescribing mature vegetation has been tested and works but it is very sensitive. If the combination of nmaxtrees, height_init_min and height_init_max are far from realistic, the model may crash (the change in KF following the change in gap probability could result in problems in allocation) and N-feedbacks will become apparent in for example the leaf mass (likely because when mature vegetation is prescribed at the start of a run, there is not enough soil N to maintain a mature vegetation). Although this was a convenient feature in ORCHIDEE-CAN, it lost much of it power after merging it into ORCHIDEE-CN-CAN. ++++++++++++

Pseudo sugar loading

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 (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. See 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. Pseudo sugar loading is now the default setting but it can be switched of by setting the flag ok_sugar_loading=n in the run.def.

Recruitment

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-CN-CAN 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. When developing and testing the code it was considered too time consuming to change the settings at the PFT level so a flag ignoring the setting for recruitment_pft was introduced. If the ok_recruitment flag is set to true, the settings of recruitment_pft will be used. If the ok_recruitment flag is set to false, the settings of recruitment_pft will be ignored and the model runs without recruitment.

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. 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.

River routing

The river routing code was not touched during the development of ORCHIDEE-CN-CAN. This functionality has not been evaluated yet for ORCHIDEE-CN-CAN, although it technically runs. Unless rivers are your main interest, set river_routing to n.

run.def

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 some driver settings (written in run.def), a bunch of 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 will generate

Soil maps

ORCHIDEE-CN-CAN runs for two soil maps: zobler (3 soil types) and usda (12 soil types). The usda map will become the default setting. The code had to be adjusted in sechiba_hydrol_architect.f90 and stomate_windthrow.f90. The soil map can be selected in the run.def by setting the value for soiltype_classif:

SOILTYPE_CLASSIF= usda
SOILCLASS_FILE=soils_param_usda.nc
SOILCLASS_BULK_FILE=soils_param.nc 

or

SOILTYPE_CLASSIF= zobler
SOILCLASS_FILE=never used zobler since the files were split
SOILCLASS_BULK_FILE=never used zobler since the files were split

Single layer vs multi layer energy budget

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. This can, however, be done by two different methods that gives 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 collapses the multi-layer energy scheme to the original enerbil code, but if you wish to take hydraulic water stress with feedback on stomatal conductance into account in your simulation, you need to choose this method.

Several flags exist to control these choices. For simplicity, a flag has been made to automatically control several of the flags related to the multi-layering, called multi_layer_control. Multi_layer_control has 4 possibilities

• 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)

The default setting of multi_layer_control in ORHIDEE-CN-CAN is 1. This is well tested and considered save to use. More details of the flags controlling the multi-layer can be found below and within the model code. The flags controlled by the multi_layer_control are:

• ok_hydrol_arch: Flag that activates the hydraulic architecture routine (true/false). The trunk version of ORCHIDEE (false) uses soil water as a proxy for water stress and applies the stress to Vcmax. When set to true the hydraulic architecture of the vegetation is accounted for to calculate the amount of water that can be transported through the plant given the soil and leaf potential and the conductivities of the roots, wood and leaves. Water supply through the plant is compared against the atmospheric demand for water. If the supply is smaller then the demand, the plant experiences water stress and the stomata will be closed (water stress is now on gs rather than Vcmax). Note that whether stomatal regulation is used or not is controled by a separate flag: ok_gs_feedback.
• 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. Be aware that if you wish to run with hydraulic architechture ok_mleb needs to be se to true as well. Furthermore, if you wish to run with the original energy scheme (enerbil), set the layers for mleb to 1.
• 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.

The multi_layer_control set these flags in the following manner:

• single layer, using the multi layer energy scheme

         ok_hydrol_arch = .TRUE.;  
         ok_gs_feedback = .TRUE.; 
         ok_mleb = .TRUE.; 
         ok_impose_can_structure = .TRUE.; 
         ok_mleb_history_file = .FALSE.
  

• multi-layer energy budget

         ok_hydrol_arch = .TRUE.; 
         ok_gs_feedback = .TRUE.; 
         ok_mleb = .TRUE.; 
         ok_impose_can_structure = .FALSE.; 
         ok_mleb_history_file = .FALSE.
  

• user specific as set in the run.def

         All read from the run.def
  

• single layer, using the original enerbil scheme

         ok_hydrol_arch = .FALSE.; 
         ok_gs_feedback = .FALSE.; 
         ok_mleb = .FALSE.; 
         ok_impose_can_structure = .FALSE.; 
         ok_mleb_history_file = .FALSE.
  

Specific leaf area

Similar to ORCHIDEE-CN, ORCHIDEE-CN-CAN 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

The analytical spin-up works with ORCHIDEE-CN-CAN. 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 systems 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.

Start and restart (under development)

The table below shows typical configurations to start and restart ORCHIDEE-CN-CAN

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 true, ok_stomate should be false. 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.

Sechiba vs stomate

ORCHIDEE-CN-CAN strengthen the links between sechiba and stomate. As in previous versions, stomate makes use of variables calculated in sechiba but in ORCHIDEE-CAN and ORCHIDEE-CN-CAN, sechiba requires information from stomate to work properly. It is possible to prescribe the canopy structure and thus only run sechiba. Set lai_map = y, the data for a canopy structure map need to come from an ORCHIDEE-CN-CAN run with stomate. All the required information is stored in the sechiba restart file to enable restarting the model without stomate.

Wind throw

The calculation of wind storm damage can be activated by setting ok_windthrow to y in the run.def. This module calculated 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 of the stands are moved to the first age class (if more than 1 age class is used).

Parameterization and Evaluation

Workflow

The following outlines the strategy for parameterizing and evaluating the performance of ORCHIDEE-CN-CAN in simulating several key forest ecosystem processes including tree growth dynamics, energy exchange, C-N cycling and plant hydrology. The work will likely require several iteration containing all or just a coupled of the following:

• running the model;
• tuning key parameters;
• re-running the model;
• spin-up and;
• evaluating to model, to discern parameter values that allow us to reproduce as closely as possible the global patterns and trends for several key processes.

The tests will move across scales starting at pixel scale moving to longitudinal band, European scale and ending at global scale. Although the longitudinal bands are of little use in the evaluation itself they can be considered pre-tests before global tests are run. It is hoped that longitudinal bands would speed up the tests and reduce the computational cost because they could show the incapability of the model to simulate major spatial gradients at a much lower cost compared to global runs.

Table 1. Workflow of the parameterization and evaluation of ORCHIDEE-CN-CAN.

Phase	Work to be done
Prepare and check	Check whether the model runs, the key functionality, and the spin-up are bug free (Table 2). Check whether the scripts, tools and data are available and still working (Table 3).
Initial evaluation	Run longitudinal spin-up and check order of magnitude in soil carbon pools. Use the same runs to check response gradients to temperature, N-deposition, precipitation and management (see Settings for spin-up).
Parameterization	Use FLUXNET and some additional data sources to parameterize the model at the site level (Table 4).
Final evaluation	Run a spin-up + transient simulation over Europe compare the simulation against the spatially explicit data (Table 5).

Apparent bug-free

Before starting the work proposed in Table 1 it needs to be confirmed that the model is technically capable of the tasks presented in Table 2.

Table 2. Essential technical capabilities before evaluating ORCHIDEE-CN-CAN.

Description of the task	Status
Is the model stable? Can it be run for 500 - 1000 years over a few individual pixels? Can it be run for all PFTs?	Recent problems with grasslands
Is the N-cycle working? Can we run simulations with IMPOSE_CN=n and IMPOSE_CN=y?	OK
Is land cover change working? Can we run simulations with (changing) land cover maps (impose_veg=no)?	Needs to be checked
Is forest management working? Can we read forest management maps, are the temporal trends in biomass what is expected for unmanaged, high-stand and coppice management?	OK
Is litter raking working? Can we read litter raking maps and is the litter pool in forest decreasing and the litter pool in croplands increasing?	OK
Does the model run in parallel?	OK
Is the analytical spin-up working?	OK
Do we have all the driver maps for the years 1600 to 2000 for European simulations? Climate, land cover change, forest management, litter raking, and N-deposition?	OK
Do we have all the driver maps for the years 1600 to 2000 for global simulations?	Still need litter raking and FM
Read N-deposition maps through COMP	OK

Data availability

Observational data products for model-data comparison can be found on OBELIX at: /home/satellites5/maignan/ORCHIDEE/EVALUATION.
1 - Phenology

• GIMMS (LAI and FPAR3g)

2 - Forest structure

• Biomass of EU forest from JRC (Europe)
• Global 1-degree Maps of Forest Area, Carbon Stocks, and Biomass, 1950-2010 (ORNL DAAC)
• Forest basal area (Europe)
• Canopy height (Global)

3 - NPP

• Site-level NPP database Luyssaert et al 2007

4 - NEP

• FLUXNET site-level data

5 - TER

• FLUXNET site-level data

6 - GPP

• FLUXNET site-level data
• EC-based upscaled GPP, i.e., Jung

7 - NPP/GPP

• site-level data and regional patterns, i.e., Campioli et al 2015

8 - Soil hydrology

• ESA CCI ECV
• measurements from Brazil (ABRACOS product)
• River discharge records from selected gauges (Global coverage)

9 - Albedo

• MODIS or GlobAlbedo? for albedo

10 - Energy (sensible and latent heat)

• GLEAM for evapotranspiration
• EC-based upscaled evapotranspiration, i.e., Jung

11 - Tree ring width (not on /home/satellites5/)

• ITRDB

12 – LCC

• Luyssaert et al 2014 – FLUXNET changes in Rn, LE, H, G, albedo
• Duveiller et al 2018 – Remote sensing changes in Rn, LE, H+G, albedo

13 – NFI

• France, Spain, Germany and Sweden (/home/satellites5/)
• EU-wide data through the VERIFY project?

There are several data tools (ATLAS, Mapper and Jerome’s) to help compare model outputs with observation, which we might be able to use. If we will not be using the available tools for comparison, we need to pre-process the observational data products to produce global means, time series, decadal averages of spatial patterns etc. The analyses presented in Naudts et al 2015 are a good starting point for the evaluation of ORCHIDEE-CN-CAN. The following scripts are available and have been tested and re-activated

Table 3. Model evaluations scripts available per November 2018

Description	Status	Path
Extract species-level productivity from French NFI data and compare against simulated productivity – looking to replace the French data by EU-wide data through the VERIFY project	OK	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/validation_species_level
Extract Jung GPP and compare against spatially explicit simulations for 8 regions in Europe	OK	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/validation_spatial_explicit
Extract Jung evapotranspiration and compare against spatially explicit simulations for 8 regions in Europe	OK	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/validation_spatial_explicit
Extract MODIS albedo (NIR & VIS) and compare against spatially explicit simulations for 8 regions in Europe	OK	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/validation_spatial_explicit
Extract height and compare against spatially explicit simulations for 8 regions in Europe	OK	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/validation_spatial_explicit
Extract effective LAI and compare against spatially explicit simulations for 8 regions in Europe	OK	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/validation_spatial_explicit
Extract BA and compare against spatially explicit simulations for 8 regions in Europe	OK	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/validation_spatial_explicit
Extract half-hourly EC observations and compare with half-hourly site level simulations for NEP, GPP, TER, and evapotranspiration	???	ENSEMBLE/FLUXNET
Plot the various components of the water budget for sites	Being developed	???

Parameterization

As we will compare the simulations to observations, the simulations need to be our best shot to resemble reality. Hence, all site-level simulations made for the parameterization will have the following configuration:

• Start from a spin-up
• impose_cn = no (i.e., accounting for N-deposition )
• impose_veg = yes (i.e., prescribe the vegetation)
• impose_soilt = yes (i.e., prescribe the soil)
• Forest management = 2 (i.e., accounting for a thin and fell type of forest management)
• litter raking = no

European simulations will have the following configuration:

• Start from a spin-up
• impose_cn = no (i.e., accounting for N-deposition )
• impose_veg = no (i.e., accounting for PFT distribution)
• impose_soilt = no (i.e., use soil map)
• Forest management from map
• litter raking = no (except for the transient simulations to construct the control)

The following parameterization approach – making use of parameters that were already shown to be sensitive to tuning – is proposed:

Table 4. Proposed order, parameters and tools to parameterize ORCHIDEE-CN-CAN. Once all steps has been performed, the NPP to GPP ratio should be re-evaluated, possibly resulting in another tuning cycle for some of the PFTs.

Process	Parameter(s)	Tool
Onset of growing season and start of senescence	thresholds for phenology	FLUXNET
NPP/GPP ratio	coeff for maintenance respiration	FLUXNET + Campioli
Magnitude of LAI	k_latosa_min and k_latosa_max	NFI + Luyssaert et al 2007
Magnitude of GPP	LL_alpha, Vcmax, J_max	FLUXNET
Magnitude of NPP	Implicit through NPP/GPP and GPP	Luyssaert et al 2007
Evapotranspiration	LAI_top, water stress	FLUXNET
Water stress	To be discovered	FLUXNET
Magnitude of Rh	To be checked - Rh	Luyssaert et al 2007
Magnitude of TER	To be checked - Rh	FLUXNET
Diameter, height, and biomass	form factor	NFI + Luyssaert et al 2007
Density and biomass	self-thinning parameters	NFI + Luyssaert et al 2007
Harvest and biomass	Self-thinning and RDI parameters	NFI + ???
Albedo	Implicit through LAI and forest structure	FLUXNET?
Tree ring width	Self-thinning, recruitment	ITRDB or VERIFY Fig 4

Settings for the FLUXNET comparison

The parameterization starts with several 1-pixel test cases coinciding with long-term flux-net sites to test whether the model captures the growth dynamics such as phenology, max LAI, GPP, etc. These tests require a spin-up. The 1-pixel test cases will allow for both parameter tuning and changes in the code to improve the model behaviour. The majority of the data represent mature forests, hence, the modelled forests should be mature as well. The model will be run for 80 years, before any output will be compared to the FLUXNET measurements. An iterative process is be planned:

• 80 years to reach mature forest → parameterize
• Re-run the 80 years to reach mature forest with the new parameters → parameterize
• Re-run spin-up and 80 year simulation to reach mature forest with the new parameters → parameterize
• Continue until satisfied

The current scripts for FLUXNET evaluation are broken down into three parts. 1) An initial looping over the driver file (1-10 years, depending on the file), 2) 500 years of spinup (regardless of the length of driver file), 3) one final loop over the driver file, for production. Tests for the grassland and cropland sites can easily use the existing setup. Given the temporal evolution of forest structure and their fluxes, testing now needs 80 years from planting after spinup as a production run over the forcing file to avoid trees dying and biasing our results. The FLUXNET script will need to be adjusted to do this. Possible issues with age classes (i.e., changes in the PFT of interest) will be avoided by using a run.def without age classes.

Only if we experience too many difficulties with manual tuning (if there are too many non-linearities in the model), we will use the multi-site optimization tool developed by Vlad . When the simulated growth dynamics are satisfying, 140 years long tests will be performed to check cumulative variables such as basal area, tree height, tree diameter, stand density, standing biomass, and harvest. To evaluate net ecosystem exchange of carbon and soil carbon and nitrogen pools a spin-up is required. Note that the spin-up depends on the parameters used in ORCHIDEE and that the sensitivity of parameters in ORCHIDEE depends on the spin-up. There is no easy way to break this dependency. We should avoid to ‘over-tune’ the 1-pixel FLUXNET comparisons. Instead, we will continue evaluating the model over longitudinal bands.

Settings for longitudinal bands

Test with longitudinal bands will be used to ensure that we have the expected gradients across regions (response to: climate, precipitation, management, and nitrogen). Longitudinal bands will not be (formally) compared against data, they should be considered as intermediate tests. If the expected responses are not present in the longitudinal bands, there is no reason to expect that they will be present in the European/global runs. Therefore, the issue should be investigated/fixed before launching large scale simulations. Three longitudinal tests will be conducted covering 2 pixels wide bands from N-Canada to Southern Chili, N-Norway to S-Africa and one band N-E Russia to Australia (including W-Australia → only place where something grows in Australia). If the evaluation is limited to checking whether response gradients are indeed present in the simulations, the configuration could be limited to running the analytical spin-up until equilibrium. Once we moved on to the longitudinal bands the focus will be on tuning of parameters. Any problems with the model and its functionality should have been caught during the 1-pixel and longitudinal tests.

Longitudinal bands can be run largely independent from the 1-pixel tests. After concluding the 1-pixels tests it probably makes sense to run the longitudinal bands to check whether the latest parameters still produce the expected response gradients.

Settings for European tests

Once happy with the pixel and longitudinal tests we will move on to the European scale. Ultimately, the aim is to produce the control run for future simulation experiments. The model should be spin-up for 1600 by making use of the 1901-1920 climatology (because those are the coldest decades in the climate forcing). During the spin-up the forest management map and N-deposition maps for 1600 should be used. As a consequence, equilibrium soil carbon should be tested at the regional scale rather than the individual pixels.

Subsequently, a transient spin-up will be run starting from the 1600 spin-up. During this spin-up, changes in forest management, CO2, N-deposition, litter-raking, and land cover change will be accounted for. The climate data will still have to cycle over 1901-1920. The transient simulation will run until 1750 and, pending on successful evaluation, extended until 1950 and 2015. 1750, 1950 and 2015 are common starting points for simulation experiments.

Whether we take 20 or 30 years climate-cycles depends on the exact length of the simulation; the target is that by the time the simulation reaches the year 1901, a cycle has been completed and so the cycle in the climate forcing is synchronized with the simulation years.

The European control should be run from the year 1601 to 2015 includes:

• A spin-up as initial condition to make a European 1° x 1° CONTROL simulation
• 64 PFTs (some with 4 age classes)
• With forest management
• Litter raking
• Increasing atmospheric CO2 concentrations read from file
• LCC
• CRU-NCEP meteorological forcing, cyclic meteorology (1901-1920) until 1900, then use the corresponding year
• A dynamic N-cycle (impose_cn=n)
• N-deposition, we have N-deposition files from approximately 1860 until present.

Table 5. Observational data products to evaluate European control simulation

Observations	Tool
EU-wide NFI data (Fig. 3) through the VERIFY project	???
EU-wide NFI data (Fig. 4) through the VERIFY project	In progress
Jung GPP – relies on FLUXNET thus partly circular	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/
Jung evapotranspiration – relies on FLUXNET thus partly circular	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/
MODIS albedo NIR & VIS	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/
Remote sensing tree height	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/
Extract effective LAI from MODIS	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/
BA from JRC – relies on NFI thus partly circular	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/
Biomass from JRC – relies on NFI thus partly circular	/ccc/work/cont003/dofoco/dofoco/SCRIPTS/
Harvest from Schelhaas - data through the VERIFY project	???
LCC – Luyssaert et al 2014	XXX
LCC – Duveiller et al 2018	???

Settings for spin-up and re-parameterization

As a spin-up is costly in both time and computer resources, we need a strategy to avoid wasting these resources. Thus, the spin-up will like the parameterization, be done across scale moving from pixel to global scale. The spin-up will be done in parallel with the parameterization. Often the problems with the spin-up have a technical characters and show up for the pixels with extreme climate conditions. Before launching a longitudinal, regional or global spin-up, we should agree on the model version to use, because structural changes to the code will necessitate re-running the spin-up. The model version to use, will most likely be the version ready once the parameterization at the 1-pixel level is satisfying, and no more changes to the code need be added.

• Identify the variables that are targeted by the spin-up (such as NEP, heterotrophic respiration, decomposition etc). The spin-up will reveal whether parameters affecting these variables need to be tuned. The spin-up itself is also an interesting test case that could be loosely compared against data.
• We could compare the spin-ups to maps of soil carbon stocks to check the order of magnitude. Soil carbon maps should only be formally compared with the control run (spin-up + transient simulation) for the year 2000 because that run includes the simulated effects of N-deposition, management, litter raking and land cover change. The more simple configurations of the spin-up do not account for these processes or do not account for the right sequence of processes.

Testing guidelines before committing new code to ORCHIDEE-CN-CAN on svn

In some rare cases after bugfixes or implementation of new code, problems with reproducibility or 1+1=2 might be introduced unintentionally. Often these are related to incorrect variable dimensions in different sub-routines, memory issues or the lack of variables in the restart files. Such issues are easier to catch sooner than later. Thus, to minimize the time spent on debugging reproducibility and 1+1=2 issues, the following simple test are suggested/required before each commit of substantial code changes*:

1+1=2

If you do not run these test globally, make sure to use impose_veg=y. The standard F2 run.def settings have been tested and 1+1=2 from revision r6272. Thus, please always make the test for the standard settings. In case of other run.def settings during your developments, make same tests for your settings also. More recent tests have shown that 1+1=2 for LCC with r6279 at the global scale.

The standard test

1) 1Y vs. 12*1M, i.e. do two simulations for a full year; one with period length of 1 year; the other with period length of 1 month. Afterwards compare their final restart files both from stomate and sechiba.

Most issues should be caught with (1). In case of problems, it will make the debugging easier, if you can track down the onset of difference between the restart files (i.e. start of year, onset of growing season, end of year etc.) Thus, continue with test like

2) 1D+1D=2D (compare the final restart files)

3) 1M+1M=2M (compare the final restart files)

How to compare netcdf files

The comparison is easiest if the same variables are contained in the two netcdf files and the variables are in the same order. The differ100.sh script by Josefine Ghattas, nicely does this. Moreover, it uses cdo diffv to compare the files. Howeve,r 5dim variables are ignored by the cdo diffv command, thus not all variables in the restart files can be compared by the differ100.sh

Have to check for differences between to netcdf files that have variables with dimensions higher than 4

The matlab function nccmp are able to compare all variables contained within two netcdf files. The original version can be found here: ​https://fr.mathworks.com/matlabcentral/fileexchange/47857-comparing-two-netcdf-files. I have made some small modifications such that the information produced by the script are put into a file instead of being printed to the screen. The update version can be found here on IRENE:/ccc/work/cont003/dofoco/dofoco/SCRIPTS/debug/nccmp.m or here on obelix:/home/data03/dofoco/SCRIPTS_obelix/debug/nccmp.m.

Sadly, matlab is not on obelix, but on IRENE. To open matlab on IRENE type Matlab or if you wish to run from the terminal type matlab -nodesktop.

Next run the function by typing:

NCCMP(ncfile1,ncfile2,tolerance,forceCompare)

Tolerance is if you allow some variation in the variables between the two files. We want identical files thus put [] here.

forceCompare can be set to true or false.

• True - write all occurrences of differences in a variable (specifically gives all the indices) to the file: all_diff.txt.

• False - only write if there is differences in a variable and its first occurrence of such differences to the file: first_diff.txt.

For global simulation the True option can produce a large file and the information might be hard to process, if there are many differences between the compared restart files. In addition, the True option makes the much script slower. However, for small simulation the true option is very useful.

I recommend that you use the re-ordered files from the difffer100.sh script as inputs to nccmp.

Debugging suggestions:

• If possible limit the spatial scale (to maximize speed).
• Track down the onset of the deviation between the restart files.
• Track down the problem. Hopefully, the differences in the restarts files will give you a clue on which variable to start the investigation from. The best approach depend on the source of the problem (memory issue or lack variable in the restart file etc.). For memory issue a debugger could be the best choice. For lack of variables in restart file it is best to run two identical runs with different period lenghts – either manually or by Totalview while tracking down which variables are causing the differences.
• Once you have fixed the problem, verify that it is also valid at the global scale (i.e. run the global tests again, if you chose to zoom in on a smaller region)

Merging CN-CAN into the trunk

In 2016 the CAN functionality was added to the CN branch with revision 3238. The end of CN branch is revision 5638. The first update of the CN-CAN branch incorporates the changes what were done to the CN branch between revision 3238 and 5638. Note that this already includes most of the changes done on the trunk because revision 5638 of the CN branch is mostly up to date with the trunk.

The path of the working version, this is the merge prepared by svn, is located at: /ccc/work/cont003/gen6328/p86ghatt/ORCHIDEE_OL/CN-CAN/modipsl-cn-can/modeles/ORCHIDEE. The working version was created with the following command: svn merge --dry-run -r 3238:5638 ​svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-CN --accept postpone

svn identified 67 files that were in conflict and 22 files that were present in one but absent in the other version. All files where svn identified conflicts as well as some files where we expected some conflicts or code that required some attention where merged manually. Whenever possible the version of CN was copied and the changes done in CN-CAN were integrated. For several of the stomate routine, the opposite approach was used.

Tickets

Not all problems could be solved on the fly. The TICKETS that still need to be taken care of were registered on the ticket server.

Notes

Not all problems could be solved on the fly. So the following NOTES may come in handy during debugging:

• sechiba.f90. The variable ksave which is present in CN-CAN is not needed in sechiba. It could be local to hydrol_arch. It was declared as an intent(in) but never calculated!

• diffuco.f90, diffuco_trans_co2. co2_to_bm is added to this routine, see ticket #419. Because the CN-CAN code closes it mass balance for both C and N, SL hopes this correction is not needed in the merge.

• diffuco.f90, diffuco_trans_co2. CN trunk uses temp_air/qair but CAN uses temp_sol/qsurf. Note before merge calculation of T_Sco and T_gm used t2m in CAN. This inconsistencies were corrected.

• slowproc.f90. cn_leaf_init_2D was not allocated or initialized. readcnleaf is not called... JG removed frac_bare because it was not used. Frac_bare is output in hydrol.f90.

• condveg.f90. JG removed use slowproc. It is not clear why it was there.

• stomate_turnover.f90 was copied from CAN. Only some small changes from CN were integrated such as using min_stomate instead of zero. Note that some arguments in IF-statements we change to account for the dimensions used in CAN.

• somate_growth_fun_all.f90. The If loop for impose_cn was taken from the CN branch therefore the structure of the CAN loops have been changed.

• stomate_lpf.f90. Bug correction: wstress_month, wstress_season were swapped in the call to stomate_lpj. Corrected in the merge.

• stomate_lpj.f90. The subroutine harvest was removed from stomate_lpj.f90 because in the merge its functionality is included in the subroutine called crop_harvest which is included in a separate module called sapiens_agriculture.f90.

• stomate_lpj.f90. The AR5 and AR6 output was added into CN-CAN. Variables calculated from prod_x and/or harvest_pool could be in the wrong units because in CN-CAN those variables were stored in gC per pixels. If the output values are 10^{10 times too big. This should be the first issue to look into.}

• stomate and stomate_lpj. Variabale names that refer explicitly to carbon, i.e., cflux were replaced by flux because the code now uses a dimension for icarbon and initrogen. The new variable names better match the content of the variables.

• stomate_soilcarbon.f90. the merge of som_dynamics started from the version of CN. The merge of som_dynamics started from the version of CN. calculation of n_minerilazation as done in CAN. The merge of nitrogen_dynamics started from the CAN version and integrated changes from CN. fluxtot has dimension also on nvm in CAN. This is not necessary as it is only used local in the nvm loop. Change in calculation of FixNH4 and f_drain as done in CN. nitrification changes, i.e. *(1.0 - anvf(:,m)) is no longer used, follows CN. kn_conv changed as in CN, mu_no3, mu_no2, mu_n2o changed as in CN, denitrification changed as in CN, frac_nh3 changed as in CN, and emm_fac is now a parameter read with getin. Its default value is 0.2 but was set to 1 in nitrogen_dynamics previously.

• constantes_mtc.f90. The following variables changed values: slainit, Vcmax25_mtc, stress_gs_mtc(from 0 to 1), stress_gm_mtc(from 0 to 1), nue_opt_mtc, height_init_max_mtc, tau_root_mtc, leafagecrit_mtc, k_root_mtc, k_sap_mtc, cn_leaf_init_mtc, alb_bg_modis(F->T), rough_dyn(F->T), sla_dyn(F->T), sandfraction_default, siltfraction_default amd, snowa_aged_mtc etc. were splitted into separate variables for vis/nir.

• readdim2.f90. The code from the latest CN version was used and we integrated changeset 5657. But one IF (dt_force >=3600 ) was removed from the trunk. This block of code was marked with a CHECK in the merged version. This was discussed with Agnes. The if that was removed in CN but not in CAN was in a section for daily_interpol and therfore dt_force was always greater than 1 hour.

• constantes.f90. FW_ change name to FWNIT_ as done in rev 4988 of the CN branch. Kept getin LAIMAX for now; it was not originally in CN. Kept IF (impveg) around getin IMPOSE_SOILT as in CN. It was removed from CAN.

• Reading of PFTmap.nc with XIOS must be checked and adapted to indice nvmap : DONE

• Added pgflux in xml files. This was missing already in the version CN-CAN before current update.