Accelerate Spinup jobs in ORCHIDEE

This page presents a work which consists to implement a numerical method that can accelerate the convergence of the SPINUP.

ORCHIDEE usually uses a iterative and specific method based on forcesoil, consisting by running the module stomate_soilcarbon.f90 stand alone between 5000 and 10000 times.

Before forcesoil, you have to run ORCHIDEE a given length just to have the litter pools at equilibrium. The problem with the FORCESOIL job is that you got some biaises and takes longer at the global scale.

That's why we implement a new method for accelerating the spinup. This method is based used by the model PASIM(send a email to get the original papers or look at the presentation available see Presentation_comité_suivi_Spinup.odp​ at the end of the page).

Description of the method used by PASIM

PASIM is a model of prairial management developped by the INRA institute.

The soil module of PASIM is based on the CENTURY model like ORCHIDEE. It consists on 5 carbon which are linked by fluxes : litter structural, litter metabolic, carbon active, carbon slow and carbon passive. we can represent symbolically the fluxes between the pools by a matrix of size (5,5).

Each line of the matrix represents a pool and each row the fluxes received (or leaving) by the pool. All these pools are set by the same differential equation.(see presentation at the end of the page).

Using the approximation of the derivative in time and two sequences, we obtain a equation with the carbon stock and these two sequences.

Mathematically, equilibrium means that carbon stock is a fixed point of the previous equation. So we obtain a linear system which can be solved by a direct algebraic method like gauss-jordan. The solution is the carbon stock at equilibrium. In PASIM, we got one linear system for each point grid. The method is stopped when a criterion is reached (relative error on the current total carbon stock calculated and the previous one).

PS : PASIM takes into account the nitrogen cycle and a similar method is used for calculating the nitrogen stock. This method can be easily extended.

Adaptation to the method for ORCHIDEE

As it is written above, the soil module of ORCHIDEE is based on the CENTURY model like PASIM. There are three main differences between PASIM and ORCHIDEE :

1. ORCHIDEE takes into account 7 pools. The metabolic litter and the structural litter are

both divided in 2 sub-pools : above and below. So the method for ORCHIDEE consists to search the following solution (litter structural above, litter structural below, litter metabolic above, litter metabolic below, carbon active, carbon slow, carbon passive)

2. There are no PFTs in PASIM. So we have to have NVM systems per each grid cell.

3. ORCHIDEE considers fire (module lpj_fire.f90). Some part of the above litter (structural + metabolic) is burned.

So we have to take into account 2 more fluxes in the matrix.

We decided not to have the same stopping condition used by PASIM. We want for the moment to evaluate the relative error for each carbon pool and to consider the maximum of all these errors.

Implementation of the method in ORCHIDEE

Here you find a text file which explains the implementation in ORCHIDEE (compare with the original papers) :

Validation of the method

The most difficult part is the validation of the method. We have to solve many problems to find the good setup. Moreover, we want ORCHIDEE to stop when the stopping criterion is reached even if the job is not finished.

Finally, we consider the following reference job :

• Full Orchidee (sechiba + stomate, no forcesoil)
  

• Calendar type : no_leap (we modify a forcing file for the year 1982)

• One grid cell ((6°-7°lon,49°-50°lat + IMPOSE_VEG=y (20% bare soil, 80% C3 grass)) + no fire

• Length of the simulation : 7000Y

• loop on the same forcing file (option norestart=y in COMP/driver.card + modification in readdim2.f90)

• Annual outputs

The results are accessible on /dmnfs12/cont003/p529sol/IGCM_OUT/OL2/OOL_SEC_STO7000Y. The results obtained are (gC/m²) :

Pools
litter structural above	2339
litter structural below	559.2
litter metabolic above	142
litter metabolic below	55
carbon active	253
carbon slow	5873
carbon passive	9226

To be precise, these values are taken from the last restart file of the run. The writing of the restart file happens at the same instant than the computation of the solution by the matrix method.

In order to compare the results, I have to take the values for the carbon pools found in the restart files (variables LITTER_xx_AB, LITTER_xx_BE ans CARBON_xx).

Results obtained by the matrix method

Here we list the results for each pool at different lenght of simulation for the matrix method :

Pools	50Y	100Y	200Y	300Y	400Y	500Y	1000Y	1500Y	2000Y	7000Y
litter structural above	2333	2339	2339	2339	2339	2339	2339	2339	2339	2339
litter structural below	558.9	559.3	559.2	559.2	559.3	559.1	559.2	559.3	559.2	559.2
litter metabolic above	142.2	142.3	142.2	142	142.3	142	142.2	142.3	142.4	142.1
litter metabolic below	55.03	55.09	55.03	54.96	55.09	55.09	55.03	55.09	55.1	55
carbon active	250.4	252.9	252.8	253	253.2	252.8	252.9	253.2	253.2	253
carbon slow	5706	5856	5873	5873	5873	5873	5873	5873	5873	5873
carbon passive	8728	8980	9111	9155	9176	9190	9215	9223	9226	9229

We notice that the convergence of the method is controlled by the carbon passive pool. We focus on the evolution of the relative error of the carbon passive obtained by the matrix method to the reference one :

100Y	200Y	300Y	400Y	500Y	1000Y	1500Y	2000Y	7000Y
2.66%	1.24%	0.76%	0.54%	0.39%	0.14%	0.05%	0.0%	0.05%

We obtained the following time series (in red the results obtained by the matrix method) :

The biaises in the previous figures can be easily explained. The values given by the algebraic method are instantaneous, because they are calculated one time at the end of the year. With the algebraic method, we are unable to represent the inter-annual variability of a carbon pool. This phenomenom is particularly striking for the carbon pools with a strong variability like the metabolic litter.

Relative error

The following figures show the evolution of the relative error for the algebraic. We assume that the value found in the last restart file of the job represents the real "equilibrium".

In order to reduce the oscillations due to the decimal part of the values, we "smooth" the curves. Without this treatment, we obtain in the case of the metabolic litter (below) :

So we use the following instruction in our ferret script :

let error_$1=100*abs($2 - $1_star[d=2,k=10,l=1:7000@sbx:40])/$2

This means we smooth our curves except the last one with the global stock. Finally, we got :

[Image(error_CARBON_SLOW_inferior_to_0.5%.gif)?

The vertical axis correspond to the error in percentage. Notice that the carbon stock could be considered at equilibrium before the passive pool is at equilibrium (we got more than 1000 years difference!).

For this reason, we have chosen in a first time step to implement a stopping condition based on the maximum of the relative errors for each pool. It is still not satisfactory. The next step will be to implement a stopping condition based on the carbon flux.

Comparison with Forcesoil

We want to compare now with the forcesoil job. In ORCHIDEE, it consists to run soilcarbon a certain number of times after and before running full ORCHIDEE. We consider the same point, we loop over the same forcing file and we use a noleap calendar. This is the spinup configuration we use :

    # SPINUP configuration : 
    # ----------------------
# Initialisation for spin-up :
# sechiba alone (!!! only if ok_stomate == n !!!)
duree_nostomate=0
# sechiba and stomate
duree_inistomate=250
# teststomate (only if duree_nostomate or duree_inistomate > 0)
duree_offlineini=0

# Loop configuration for spin-up :
# The whole job is restarted n_iter times
n_iter=1
# orchidee with sechiba (and stomate if ok_stomate=y below)
duree_sechiba=10
# teststomate
duree_stomate=0
# forcesoil
duree_carbonsol=10000

# Finalization for spin-up :
# all orchidee
duree_final=250
# This last parameter must be non-zero.

We obtain the following results at the end of the simulation:

Pools
litter structural above	2339
litter structural below	559.2
litter metabolic above	142
litter metabolic below	54.9
carbon active	253
carbon slow	5873
carbon passive	9223

The output files can be found at : /dmnfs12/cont003/p529sol/IGCM_OUT/OL2/SpinUp_classic_setup_nofire

Forcesoil is very performing. Let's take the example of the passive pool.

ORCHIDEE+forcesoil runs only 500Y to find a error of 0.05%.

The algebraic method have to run 1500Y to have a similar error (and 500Y more to find the exact result) and ORCHIDEE alone got the same value as forcesoil after 6200Y.

Our analytic method is 3 to 4 more faster than ORCHIDEE alone but less performing than forcesoil.

The analytic method found the same results after 500Y than ORCHIDEE+forcesoil in 51Y.
(see /dmnfs12/cont003/p529sol/IGCM_OUT/OL2/50Y_ORC_FORCESOIL_2/)

Another test

A another test was made using the same point and the pft 6. The results obtained are globally the same : they are not detailed but could be found here : /dmnfs12/cont003/p529sol/IGCM_OUT/OL2/spinup_7000Y_pft6.tar

Summary

Work done

• Implementation of the method of PASIM in ORCHIDEE : adaptation to 7 pools + fire flux
  

• Implementation of a stopping condition based on the maximum of all the relative errors calculated for each pool
  

• Modification of the script and job of OOL_SEC_STO : if the stopping condition chosen by the user is reached before the end of the run, the job will stop
  

• Tests made on the temperate regions for pfts 10 and 6
  

• Creation of a branche on ths svn server called Spinup_analytic based on the externalized version

To be done

• New tests on a boreal region, tropical region and global
  

• Implementation of a new stopping condition based on the evaluation of the flux and not on the stock

MEETING (20/09/2011)

A meeting was organized in order to discuss the following problem : why is the algebraic method not more faster?

Presents : Nicolas Viovy, Nicolas Vuichard, Philippe Peylin Bertrand Guenet and Martial Mancip

After presenting the method in details (see presentation), we conclude that :

• the matrix method should not start at the beginning of the run but after the litter and the biomass are at equilibrium
• that the stopping criterion should be based on the fluxes ( maybe two criterions in percentage and gC)
• should test others points (boreal and tropical)

Results :

We launched ORCHIDEE for 250Y on the same point described above then we restart and solve the system for 5 years. We obtained the same results as ORCHIDEE in 7000Y (we tested the PFT 6 also)!

In order to optimize the length of the run, we have planned to add :

• a starting condition for the analytic solving of the system based on the equilibrium reached by the biomass and the litter (the length should not be more than 250Y)
• when the starting condition is OK, we solve the system until the stopping condition is reached (that will stop the job)

PS

• If you want to know more about this work, look at the presentation below and/or send a email to Didier.Solyga.at.lsce.ipsl.fr