= 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. [[BR]] 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.[[BR]] 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.[[BR]] 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 [attachment: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.[[BR]] 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 :[[BR]] - Full Orchidee (sechiba + stomate, no forcesoil)[[BR]] - 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.[[BR]] 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) : [[Image(comp_LITTER_STR_AB_LITTER_STR_AB_star.gif)]] [[Image(comp_LITTER_STR_BE_LITTER_STR_BE_star.gif)]] [[Image(comp_LITTER_MET_AB_LITTER_MET_AB_star.gif)]] [[Image(comp_LITTER_MET_BE_LITTER_MET_BE_star.gif)]] [[Image(comp_CARBON_ACTIVE_CARBON_ACTIVE_star.gif)]] [[Image(comp_CARBON_SLOW_CARBON_SLOW_star.gif)]] [[Image(comp_CARBON_PASSIVE_CARBON_PASSIVE_star.gif)]] [[Image(comparaison_evolution_stock.gif)]] 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) : [[Image(error_LITTER_MET_AB_inferior_to_0.5%.gif)]] 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_LITTER_STR_AB_inferior_to_0.2%.gif)]] [[Image(error_LITTER_STR_BE_inferior_to_0.2%.gif)]] [[Image(error_LITTER_MET_AB_inferior_to_0.2%.gif)]] [[Image(error_LITTER_MET_BE_inferior_to_0.5%.gif)]] [[Image(error_CARBON_ACTIVE_inferior_to_0.5%.gif)]] [[[Image(error_CARBON_SLOW_inferior_to_0.5%.gif)]] [[Image(error_CARBON_PASSIVE_inferior_to_2%.gif)]] [[Image(error_stock.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!). [[BR]] 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.[[BR]] ORCHIDEE+forcesoil runs only 500Y to find a error of 0.05%.[[BR]] 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. [[BR]] Our analytic method is 3 to 4 more faster than ORCHIDEE alone but less performing than forcesoil. [[BR]] The analytic method found the same results after 500Y than ORCHIDEE+forcesoil in 51Y.[[BR]] (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 .[[BR]] 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 [[BR]] - Implementation of a stopping condition based on the maximum of all the relative errors calculated for each pool[[BR]] - 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[[BR]] - Tests made on the temperate regions for pfts 10 and 6 [[BR]] - 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[[BR]] - 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)! [[BR]] 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).[[BR]] For testing this point, I make a 250Y run without the method. Then I solve for 5Y after restarting this job. [[BR]] [[BR]] For the PFT 10 and the PFT 6, I obtain good results (0.05% relative error to the referenced values - 7000Y ORCHIDEE run). [[BR]] We divide the simulation lenght by 6 for the grasses and by 12 to have the same relative error (see above)! [[BR]] 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.[[BR]] 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.[[BR]] 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). [[BR]] 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.[[BR]] 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). [[BR]] 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. [[BR]] We choose to set the threshold for the carbon to 0.01%.[[BR]] I make a run whose the simulation lenght was 500Y max. [[BR]] If the biomass equilibrium is not reached, the resolution start after 495Y for the last five years.[[BR]] 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! [[BR]] Because of the grasses and of the intrisic instabilities of the model, we decide to make a third change. [[BR]] For the biomass criterions, we focus now only on the fluxes : turnover_daily_accu and bm_to_litter_accu. [[BR]] 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.[[BR]] 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. [[BR]] 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. [[BR]] 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 || PASSIVE THRESHOLD || SPINUP LENGTH || SOLUTION || || 0.1% || 0.01% || 302Y || 2983 | 1154 | 208.2| 119.1 | 388 | 8208 | 12849 || Results for PFT 10 (7E-50N): || BIOMASS THRESHOLD || PASSIVE THRESHOLD || SPINUP LENGTH || SOLUTION || || 0.1% || 0.01% || ??Y || 2428| 559.7 | 112.1 | 48.93 | 254.8 | 5821 | 9149 || Results for PFT 6 (50%) + PFT 10 (50%) (7E-50N): || BIOMASS THRESHOLD || PASSIVE THRESHOLD || SPINUP LENGTH || SOLUTION || || 0.1% || 0.01% || 352Y || 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 || PASSIVE THRESHOLD || SPINUP LENGTH || SOLUTION || || 0.1% || 0.01% || 252Y || 1781 | 443.8 | 164.1 | 31.23 | 218.8 | 3499 | 5518 || Results for PFT11 - C4 grass (8W-10N) : || BIOMASS THRESHOLD || PASSIVE THRESHOLD || SPINUP LENGTH || SOLUTION || || 0.1% || 0.01% || 47Y || 2534 | 307.6 | 24.23 | 17.6 | 160.9 | 4816 | 7519 || Results for PFT7 - Boreal needleleaf evergreen (7E-60N) : || BIOMASS THRESHOLD || PASSIVE THRESHOLD || SPINUP LENGTH || SOLUTION || || 0.1% || 0.01% || 542Y|| 2534 | 307.6 | 24.23 | 17.6 | 160.9 | 4816 | 7519 || ''' Conclusion :''' The reset works for the boreal region (1 resolution less) but not on the PFT 2 (2 resolutions more).[[BR]] I implemented also a diagnostic variable nbp_flux to see the evolution of the nbp during the forcing period. === 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