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 .
Here are the results (7000Y run) :

Pools
litter structural above	2606
litter structural below	1013
litter metabolic above	196
litter metabolic below	107
carbon active	326
carbon slow	7251
carbon passive	11372

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 :

I launched ORCHIDEE for 250Y on the same point described above then I 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)

WORK DECEMBER 2011/JANUARY 2012

After the meeting of September 2011, we agree that is unuseful to start the algebraic method without the equilibrium of the biomass (litter inputs).

For testing this point, I make a 250Y run without the method. Then I solve for 5Y after restarting this job.

For the PFT 10 and the PFT 6, I obtain good results (0.05% relative error to the referenced values - 7000Y ORCHIDEE run).
We divide the simulation lenght by 6 for the grasses and by 12 to have the same relative error (see above)!

It proves well that the analytic spinup is very sensitive to the biomass and its equilibrium.

1) First implementation :

In order to optimize the time for the equilibrium of the biomass, we look at 3 criterions (2 fluxes and one stock):

• turnover_daily_accu : the cumulated sum of the turnover daily over the forcing period
• bm_to_litter_accu : the cumulated sum of bm_to_litter over the forcing period
• biomass : the stock of biomass

If the relatives errors for the three previous criterions between two forcing period is less than a given threshold (in %) simultaneously, we start the algebraic method.

The error is calculated by taking the maximum relative error over all the biomass compartment. To simlify the survey, we set the same

threshold for the three criterions.

If this condition cannot be reached, we force after the resolution after a prescribed length defined by the user (typically the estimated time for the biomass to be at equilibrium).

Similarly, the carbon equilibrium is not reached until a threshold defined by the user is reached. We consider the maximum relative error over the carbon pools.

The simulation finishes when the equilibrium is reached or normally.

Conclusion :

I test this method for the PFT 6 and PFT 10.

I conclude that we cannot set a strict criterion over the carbon (0.01%) and over the biomass (0.1% didn't work) due to the reactif pools (metabolic litter : a little variation implies a "big" relative error).
So we decide to focus only to the passive pool for computing the error.

2) Second implementation :

The same as above except the relative error for the carbon is based on the passive pool only.
We choose to set the threshold for the carbon to 0.01%.

I make a run whose the simulation lenght was 500Y max.
If the biomass equilibrium is not reached, the resolution start after 495Y for the last five years.

Here are the results for the PFT 6 :

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	SOLUTION
0.5%	85Y	0.01%	246Y	2617.5 | 1020.1 | 195 | 106.8 | 328.2 | 7269.6 | 11268
0.3%	103Y	0.01%	222Y	2611.2 | 1020.6 | 194.9 | 106.8 | 327.8 | 7252.2 | 11295.8
0.1%	495Y (forced resolution)	0.01%	498Y	2625.1 | 1019.4 | 195.2 | 106.8 | 328.8 | 7295.9 | 11448

Conclusion :

The threshold 0.1 % for the biomass stock cannot be reached so I have to wait 200Y for the equilibrium of the biomass for the grasses!
Because of the grasses and of the intrisic instabilities of the model, we decide to make a third change.
For the biomass criterions, we focus now only on the fluxes : turnover_daily_accu and bm_to_litter_accu.
We look at the relative errors every 10 years (decennal mean of the fluxes).

3) Third implementation :

In order to reduce intrisic instabilities of the model which affect the method, we look at the relative error on the decennal mean fo the biomass fluxes.

For the same reason, we solve the spinup system every 10 years after the starting of the method. I used the same protocol for the PFT 6 and 10 described above :

Results for PFT 6 (7E-50N):

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.5%	90Y	0.01%	330Y	15	2624.3 | 1019.5 | 195.3 | 106.7 | 328.8 | 7293.2 | 11433.3
0.25%	210Y	0.01%	311Y	11	2623.7 | 1019.6 | 195.3 | 106.7 | 328.7 | 7291.8 | 11436.2
0.15%	230Y	0.01%	241Y	2	2617.9 | 1020.1 | 195.1 | 106.8 | 328.4 | 7284.6 | 11430.6
0.1%	250Y	0.01%	260Y	2	2620.8 | 1019.8 | 195.2 | 106.8 | 328.6 | 7289.2 | 11437.7
0.01%	340Y	0.01%	361Y	3	2624.9 | 1019.5 | 195.3 | 106.7 | 328.8 | 7295.6 | 11447.5

Results for PFT 10 (maximum length : 200Y ; forcing resolution : 195Y) (only one test made)(7E-50N):

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	40Y	0.01%	60Y	3	2339.4 | 559.2 | 142.3 | 55.1 | 253.1 | 5873.7 | 9229.9

Results for PFT 6 (50%) + PFT 10 (50%) (7E-50N):

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	250	0.01%	281	4	2612.53 | 1013.9 | 194.4 | 107.7 | 324.8 | 7204.2 | 11306.8
					2691.13 | 516.3 | 154.4 | 47 |273.7 | 6411 | 10053.6

Results for PFT2 - tropical evergreen (65W-5S) :

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	210Y	0.01%	231Y	3	1597.3 | 373.1 | 58.8 | 13.2 | 137.9 | 3256.2 | 5059.2

Results for PFT11 - C4 grass (8W-10N) :

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	30Y	0.01%	41Y	2	3332.2 | 325.8 | 359.8 | 44.8 | 220.9 | 5501.6 | 8628

Results for PFT7 - Boreal needleleaf evergreen (7E-60N) :

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	410Y	0.01%	511Y	11	2872.2 | 698 | 122.2 | 32.3 | 270.6 | 6080.7 | 9527

Conclusion :

We decide now to reduce the number of system resolutions.
For that, I will implement a "reset" button : if the relative error on the passive pool is too high, I reset the matrix to zero. I hope to reduce the length simulation.
That means also we may not need biomass criterion anymore.

4) Fourth implementation :

Implementation of the reset for the matrix. Result for PFT 2 - tropical evergreen (65W-5S) :

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	210Y	0.01%	251Y	5	1597.7 | 372.9 | 58.8 | 13.2 | 138 | 3257 | 5060.6

Result for PFT 7 - Boreal needleleaf evergreen (7E-60N) :

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	410Y	0.01%	501Y	10	2872.5 | 698 | 122.3 | 32.3 | 270.6 | 6083.5 | 9534.1

5) Results for multi forcing years jobs :

After modifying the spinup job, I launched the same simulations for 5 years forcing years (NCC modified, 1980-1984, noleap). I obtain the following results :

Results for PFT 6 (7E-50N):

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	251Y	0.01%	302Y	6	2983 | 1154 | 208.2| 119.1 | 388 | 8208 | 12849

Results for PFT 10 (7E-50N):

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	31Y	0.01%	222Y	20 (?)	2428| 559.7 | 112.1 | 48.93 | 254.8 | 5821 | 9149

Results for PFT 6 (50%) + PFT 10 (50%) (7E-50N):

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	251Y	0.01%	352Y	11	2947 | 1130| 206.8 | 122 | 369.1 | 7940 | 12445
					2854 | 528 | 101.1 | 36.9 | 309 | 6684 | 10467

Results for PFT2 - tropical evergreen (65W-5S) :

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	211Y	0.01%	252Y	5	1781 | 443.8 | 164.1 | 31.23 | 218.8 | 3499 | 5518

Results for PFT11 - C4 grass (8W-10N) :

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	31Y	0.01%	47Y	4	2534 | 307.6 | 24.23 | 17.6 | 160.9 | 4816 | 7519

Results for PFT7 - Boreal needleleaf evergreen (7E-60N) :

BIOMASS THRESHOLD	BIOMASS EQUILIBRIUM TIME	PASSIVE THRESHOLD	SPINUP LENGTH	NUMBER OF SYSTEM RESOLUTIONS	SOLUTION
0.1%	411Y	0.01%	542Y	14	2636 | 703.5 | 90.66 | 27.8 | 263.7 | 6127 | 9590

Conclusion :

The reset works for the boreal region (1 resolution less) but not on the PFT 2 (2 resolutions more).
I implemented also a diagnostic variable nbp_flux to see the evolution of the nbp during the forcing period.

TEST OVER EUROPE :

Two tests over Europe were done with this version. I used 5 forcing years (1980-1984) for a 300Y simulation (prescribed resolution at 280Y).
In the first test, we use the complete algorithm of the spinup with strict criterions over biomass(0,1%) and passive pool(0,01%).
For the second test, we set a huge error on the biomass so the method can start immediately. I obtained in this case 3 snapshots at 100Y, 200Y and 300Y (see figures below).
I plot for each time the four following variables : total_carbon_soil, nbp, nbp/gpp et nbp/total_carbon_stock_prev (which corresponds to the variation of carbon between two resolutions).
Results after 100Y :

Results after 200Y :

Final results :

MEETING 13/04/2012

Presents : Ncolas Vuichard, Nicolas Viovy, Philippe Peylin, Fabienne Maignan, Patricia Cadule see Spinup_13042012.pdf​.

Actions :

• Simplification of the spinup algorithm : deletion of the biomass tests; start the resolution from scratch; reset the method at every period; still solve the system even if the point is at equiibrium
• Diagnostic : keep relative error on passive pool and nbp value diagnostic; the relative error is calculated over the sum of the passive pools over the pixel and not by pft anymore.
• TO DO :Test in parallel version + coupling with teststomate + improve scripts.

FINAL PRESENTATION ABOUT THIS WORK (05/07/2012)

The work was finalized between April and June 2012 and was presented during an ORCHIDEE-DEV meeting in July. The presentation is available here : ​SPINUP_final_presentation

SYNTHESIS DOCUMENT (SEPTEMBER 2012)

A technical document summarizing the work is available here : ​SPINUP_technical_documentation
It is essentially based on the presentation given in July 2012 and contains information about the mathematical principles, the modifications done to the code, etc. Extension to forest management : see here Branches/AccelerationSpinup/ExtensionMethod

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