source: branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate.f90 @ 6319

Last change on this file since 6319 was 6319, checked in by josefine.ghattas, 5 years ago

Improvment for the ESM CO2 configuration:

  • Separate variable fco2_lu into 3 parts: fco2_lu, fco2_wh and fco2_ha
  • Move calculation of co2_flux from dt_sechiba time-step to daily time-step (in the part for do_slow)
  • Removed co2_flux and fco2_lu from stomate_intialize argument list. These variables were never used in the intialization phase.
  • Add co2_flux, and fco2_wh, fco2_ha to restart file
  • Corrected output unit for nee to be consistent with LMDZ and stomate output variables. It is now in kgC/m2/s.
  • Corrected output for znetco2
  • Added fCO2_fWoodharvest and fCO2_fHarvest as new possible tracers in LMDZ (intersurf).
  • Added diagnostic output for fCO2_fWoodharvest and fCO2_fHarvest
  1. Cadule
  • Property svn:keywords set to HeadURL Date Author Revision
File size: 241.1 KB
Line 
1! =================================================================================================================================
2! MODULE       : stomate
3!
4! CONTACT      : orchidee-help _at_ listes.ipsl.fr
5!
6! LICENCE      : IPSL (2006)
7! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8!
9!>\BRIEF       Groups the subroutines that: (1) initialize all variables in
10!! stomate, (2) read and write forcing files of stomate and the soil component,
11!! (3) aggregates and convert variables to handle the different time steps
12!! between sechiba and stomate, (4) call subroutines that govern major stomate
13!! processes (litter, soil, and vegetation dynamics) and (5) structures these tasks
14!! in stomate_main
15!!
16!!\n DESCRIPTION : None
17!!
18!! RECENT CHANGE(S) : None
19!!
20!! REFERENCE(S) : None
21!!
22!! SVN :
23!! $HeadURL$
24!! $Date$
25!! $Revision$
26!! \n
27!_ ================================================================================================================================
28
29MODULE stomate
30
31  ! Modules used:
32  USE netcdf
33  USE defprec
34  USE grid
35  USE time, ONLY : one_day, one_year, dt_sechiba, dt_stomate, LastTsYear, LastTsMonth
36  USE time, ONLY : year_end, month_end, day_end, sec_end
37  USE constantes
38  USE constantes_soil
39  USE pft_parameters
40  USE stomate_io
41  USE stomate_data
42  USE stomate_season
43  USE stomate_lpj
44  USE stomate_litter
45  USE stomate_vmax
46  USE stomate_soilcarbon
47  USE stomate_resp
48  USE mod_orchidee_para
49  USE ioipsl_para 
50  USE xios_orchidee
51
52  USE matrix_resolution
53 
54  IMPLICIT NONE
55
56  ! Private & public routines
57
58  PRIVATE
59  PUBLIC stomate_main,stomate_clear,init_forcing, stomate_forcing_read, stomate_initialize, stomate_finalize
60
61  INTERFACE stomate_accu
62     MODULE PROCEDURE stomate_accu_r1d, stomate_accu_r2d, stomate_accu_r3d
63  END INTERFACE
64
65  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:):: biomass              !! Biomass per ground area @tex $(gC m^{-2})$ @endtex
66!$OMP THREADPRIVATE(biomass)
67  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: veget_cov_max        !! Maximal fractional coverage: maximum share of a pixel
68                                                                         !! taken by a PFT
69!$OMP THREADPRIVATE(veget_cov_max)
70  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: ind                  !! Vegetation density, number of individuals per unit
71                                                                         !! ground area @tex $(m^{-2})$ @endtex
72!$OMP THREADPRIVATE(ind)
73  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: age                  !! Age of PFT it normalized by biomass - can increase and
74                                                                         !! decrease - (years)
75!$OMP THREADPRIVATE(age)
76  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: adapted              !! Winter too cold for PFT to survive (0-1, unitless)
77!$OMP THREADPRIVATE(adapted)
78  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: regenerate           !! Winter sufficiently cold to produce viable seeds
79                                                                         !! (0-1, unitless)
80!$OMP THREADPRIVATE(regenerate)
81  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: everywhere           !! Is the PFT everywhere in the grid box or very localized
82                                                                         !! (after its intoduction)
83!$OMP THREADPRIVATE(everywhere)
84  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: fireindex            !! Probability of fire (unitless)
85!$OMP THREADPRIVATE(fireindex)
86  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: veget_lastlight      !! Vegetation fractions (on ground) after last light
87                                                                         !! competition (unitless)
88!$OMP THREADPRIVATE(veget_lastlight)
89  REAL(r_std), ALLOCATABLE,SAVE,DIMENSION(:,:)   :: fpc_max              !! "maximal" coverage fraction of a grid box (LAI ->
90                                                                         !! infinity) on ground. [??CHECK??] It's set to zero here,
91                                                                         !! and then is used once in lpj_light.f90 to test if
92                                                                         !! fpc_nat is greater than it. Something seems missing
93!$OMP THREADPRIVATE(fpc_max)
94  LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:,:)        :: PFTpresent           !! PFT exists (equivalent to veget > 0 for natural PFTs)
95!$OMP THREADPRIVATE(PFTpresent)
96  LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:,:)        :: senescence           !! The PFT is senescent
97!$OMP THREADPRIVATE(senescence)
98  LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:,:)        :: begin_leaves         !! Signal to start putting leaves on (true/false)
99!$OMP THREADPRIVATE(begin_leaves)
100  LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:,:)        :: need_adjacent        !! This PFT needs to be in present in an adjacent gridbox
101                                                                         !! if it is to be introduced in a new gridbox
102!$OMP THREADPRIVATE(need_adjacent)
103!--
104  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: humrel_daily         !! Daily plant available water -root profile weighted
105                                                                         !! (0-1, unitless)
106!$OMP THREADPRIVATE(humrel_daily)
107  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: humrel_week          !! "Weekly" plant available water -root profile weighted
108                                                                         !! (0-1, unitless)
109!$OMP THREADPRIVATE(humrel_week)
110  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: humrel_month         !! "Monthly" plant available water -root profile weighted
111                                                                         !! (0-1, unitless)
112!$OMP THREADPRIVATE(humrel_month)
113  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: maxhumrel_lastyear   !! Last year's max plant available water -root profile
114                                                                         !! weighted (0-1, unitless)
115!$OMP THREADPRIVATE(maxhumrel_lastyear)
116  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: maxhumrel_thisyear   !! This year's max plant available water -root profile
117                                                                         !! weighted (0-1, unitless)
118!$OMP THREADPRIVATE(maxhumrel_thisyear)
119  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: minhumrel_lastyear   !! Last year's min plant available water -root profile
120                                                                         !! weighted (0-1, unitless) 
121!$OMP THREADPRIVATE(minhumrel_lastyear)
122  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: minhumrel_thisyear   !! This year's minimum plant available water -root profile
123                                                                         !! weighted (0-1, unitless)
124!$OMP THREADPRIVATE(minhumrel_thisyear)
125!--- 
126  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: t2m_daily            !! Daily air temperature at 2 meter (K)
127!$OMP THREADPRIVATE(t2m_daily)
128
129  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: Tseason              !! "seasonal" 2 meter temperatures (K)
130!$OMP THREADPRIVATE(Tseason)
131  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: Tseason_length       !! temporary variable to calculate Tseason
132!$OMP THREADPRIVATE(Tseason_length)
133  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: Tseason_tmp          !! temporary variable to calculate Tseason
134!$OMP THREADPRIVATE(Tseason_tmp)
135  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: Tmin_spring_time     !! Number of days after begin_leaves (leaf onset)
136!$OMP THREADPRIVATE(Tmin_spring_time)
137  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: onset_date           !! Date in the year at when the leaves started to grow(begin_leaves), only for diagnostics.
138!$OMP THREADPRIVATE(onset_date)
139  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: t2m_week             !! Mean "weekly" (default 7 days) air temperature at 2
140                                                                         !! meter (K) 
141!$OMP THREADPRIVATE(t2m_week)
142  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: t2m_month            !! Mean "monthly" (default 20 days) air temperature at 2
143                                                                         !! meter (K)
144!$OMP THREADPRIVATE(t2m_month)
145  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: t2m_longterm         !! Mean "Long term" (default 3 years) air temperature at
146                                                                         !! 2 meter (K)
147!$OMP THREADPRIVATE(t2m_longterm)
148  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: t2m_min_daily        !! Daily minimum air temperature at 2 meter (K)
149!$OMP THREADPRIVATE(t2m_min_daily)
150  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: tsurf_daily          !! Daily surface temperatures (K)
151!$OMP THREADPRIVATE(tsurf_daily)
152!---
153  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: precip_daily         !! Daily precipitations sum @tex $(mm day^{-1})$ @endtex
154!$OMP THREADPRIVATE(precip_daily)
155  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: precip_lastyear      !! Last year's annual precipitation sum
156                                                                         !! @tex $??(mm year^{-1})$ @endtex
157!$OMP THREADPRIVATE(precip_lastyear)
158  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: precip_thisyear      !! This year's annual precipitation sum
159                                                                         !! @tex $??(mm year^{-1})$ @endtex
160!$OMP THREADPRIVATE(precip_thisyear)
161!---
162  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: soilhum_daily        !! Daily soil humidity (0-1, unitless)
163!$OMP THREADPRIVATE(soilhum_daily)
164  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: soilhum_month        !! Soil humidity - integrated over a month (0-1, unitless)
165!$OMP THREADPRIVATE(soilhum_month)
166  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: tsoil_daily          !! Daily soil temperatures (K)
167!$OMP THREADPRIVATE(tsoil_daily)
168  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: tsoil_month          !! Soil temperatures at each soil layer integrated over a
169                                                                         !! month (K)
170!$OMP THREADPRIVATE(tsoil_month)
171!---
172  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: litterhum_daily      !! Daily litter humidity (0-1, unitless)
173!$OMP THREADPRIVATE(litterhum_daily)
174!---
175  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: control_moist        !! Moisture control of heterotrophic respiration
176                                                                         !! (0-1, unitless)
177!$OMP THREADPRIVATE(control_moist)
178  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: control_temp         !! Temperature control of heterotrophic respiration at the
179                                                                         !! different soil levels (0-1, unitless)
180!$OMP THREADPRIVATE(control_temp)
181  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: control_moist_daily  !! Moisture control of heterotrophic respiration daily
182                                                                         !! (0-1, unitless)
183!$OMP THREADPRIVATE(control_moist_daily)
184  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: control_temp_daily   !! Temperature control of heterotrophic respiration, above
185                                                                         !! and below daily (0-1, unitless)
186!$OMP THREADPRIVATE(control_temp_daily)
187!---
188  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: gdd_init_date        !! inital date for gdd count
189!$OMP THREADPRIVATE(gdd_init_date)
190
191  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: gdd_from_growthinit  !! gdd from beginning of season (C)
192!$OMP THREADPRIVATE(gdd_from_growthinit)
193  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: gdd0_lastyear        !! Last year's annual Growing Degree Days,
194                                                                         !! threshold 0 deg C (K)
195!$OMP THREADPRIVATE(gdd0_lastyear)
196  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: gdd0_thisyear        !! This year's annual Growing Degree Days,
197                                                                         !! threshold 0 deg C (K)
198!$OMP THREADPRIVATE(gdd0_thisyear)
199  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: gdd_m5_dormance      !! Growing degree days for onset of growing season,
200                                                                         !! threshold -5 deg C (K)
201!$OMP THREADPRIVATE(gdd_m5_dormance)
202  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: gdd_midwinter        !! Growing degree days for onset of growing season,
203                                                                         !! since midwinter (K)
204!$OMP THREADPRIVATE(gdd_midwinter)
205  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: ncd_dormance         !! Number of chilling days since leaves were lost (days)
206!$OMP THREADPRIVATE(ncd_dormance)
207  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: ngd_minus5           !! Number of growing days, threshold -5 deg C (days)
208!$OMP THREADPRIVATE(ngd_minus5)
209  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: hum_min_dormance     !! Minimum moisture during dormance (0-1, unitless)
210!$OMP THREADPRIVATE(hum_min_dormance)
211!---
212  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: gpp_daily            !! Daily gross primary productivity per ground area
213                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
214!$OMP THREADPRIVATE(gpp_daily)
215  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: gpp_week             !! Mean "weekly" (default 7 days) GPP 
216                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
217!$OMP THREADPRIVATE(gpp_week)
218  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: maxgppweek_lastyear  !! Last year's maximum "weekly" GPP 
219                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
220!$OMP THREADPRIVATE(maxgppweek_lastyear)
221  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: maxgppweek_thisyear  !! This year's maximum "weekly" GPP 
222                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex 
223!$OMP THREADPRIVATE(maxgppweek_thisyear)
224!---
225  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: npp_daily            !! Daily net primary productivity per ground area
226                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
227!$OMP THREADPRIVATE(npp_daily)
228  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: npp_longterm         !! "Long term" (default 3 years) net primary productivity
229                                                                         !! per ground area 
230                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex   
231!$OMP THREADPRIVATE(npp_longterm)
232  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: npp_equil            !! Equilibrium NPP written to forcesoil
233                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
234!$OMP THREADPRIVATE(npp_equil)
235  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: npp_tot              !! Total NPP written to forcesoil
236                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
237!$OMP THREADPRIVATE(npp_tot)
238!---
239  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: resp_maint_part_radia!! Maintenance respiration of different plant parts per
240                                                                         !! total ground area at Sechiba time step 
241                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
242!$OMP THREADPRIVATE(resp_maint_part_radia)
243  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: resp_maint_part      !! Maintenance respiration of different plant parts per
244                                                                         !! total ground area at Stomate time step
245                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
246!$OMP THREADPRIVATE(resp_maint_part)
247  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: resp_maint_radia     !! Maintenance respiration per ground area at Sechiba time
248                                                                         !! step   
249                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
250!$OMP THREADPRIVATE(resp_maint_radia)
251  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: resp_maint_d         !! Maintenance respiration per ground area at Stomate time
252                                                                         !! step 
253                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
254!$OMP THREADPRIVATE(resp_maint_d)
255  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: resp_growth_d        !! Growth respiration per ground area
256                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
257!$OMP THREADPRIVATE(resp_growth_d)
258  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: resp_hetero_d        !! Heterotrophic respiration per ground area
259                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
260!$OMP THREADPRIVATE(resp_hetero_d)
261  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: resp_hetero_radia    !! Heterothrophic respiration per ground area at Sechiba
262                                                                         !! time step
263                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
264!$OMP THREADPRIVATE(resp_hetero_radia)
265!---
266  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)     :: turnover_time       !! Turnover time of grasses
267                                                                         !! @tex $(dt_stomate^{-1})$ @endtex
268!$OMP THREADPRIVATE(turnover_time)
269  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:) :: turnover_daily      !! Senescence-driven turnover (better: mortality) of
270                                                                         !! leaves and roots 
271                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
272!$OMP THREADPRIVATE(turnover_daily)
273  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:) :: turnover_littercalc !! Senescence-driven turnover (better: mortality) of
274                                                                         !! leaves and roots at Sechiba time step
275                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
276!$OMP THREADPRIVATE(turnover_littercalc)
277  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:) :: turnover_longterm   !! "Long term" (default 3 years) senescence-driven
278                                                                         !! turnover (better: mortality) of leaves and roots
279                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
280!$OMP THREADPRIVATE(turnover_longterm)
281  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:) :: bm_to_litter        !! Background (not senescence-driven) mortality of biomass
282                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
283!$OMP THREADPRIVATE(bm_to_litter)
284  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:) :: bm_to_littercalc    !! conversion of biomass to litter per ground area at
285                                                                         !! Sechiba time step
286                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
287!$OMP THREADPRIVATE(bm_to_littercalc)
288  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: dead_leaves          !! Metabolic and structural pools of dead leaves on ground
289                                                                         !! per PFT @tex $(gC m^{-2})$ @endtex
290!$OMP THREADPRIVATE(dead_leaves)
291  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:,:):: litter             !! Above and below ground metabolic and structural litter
292                                                                         !! per ground area
293                                                                         !! @tex $(gC m^{-2})$ @endtex
294!$OMP THREADPRIVATE(litter)
295  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: litterpart           !! Fraction of litter above the ground belonging to
296                                                                         !! different litter pools (unitless)
297!$OMP THREADPRIVATE(litterpart)
298  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: firelitter           !! Total litter above the ground that could potentially
299                                                                         !! burn @tex $(gC m^{-2})$ @endtex
300!$OMP THREADPRIVATE(firelitter)
301  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:):: soilcarbon_input     !! Quantity of carbon going into carbon pools from litter
302                                                                         !! decomposition per ground area  at Sechiba time step
303                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex
304!$OMP THREADPRIVATE(soilcarbon_input)
305  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: soilcarbon_input_daily !! Daily quantity of carbon going into carbon pools from
306                                                                           !! litter decomposition per ground area
307                                                                           !! @tex $(gC m^{-2} day^{-1})$ @endtex
308!$OMP THREADPRIVATE(soilcarbon_input_daily)
309  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: carbon               !! Soil carbon pools per ground area: active, slow, or
310                                                                         !! passive, @tex $(gC m^{-2})$ @endtex
311!$OMP THREADPRIVATE(carbon)
312  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: lignin_struc         !! Ratio Lignine/Carbon in structural litter for above and
313                                                                         !! below ground compartments (unitless)
314!$OMP THREADPRIVATE(lignin_struc)
315  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: lm_lastyearmax       !! Last year's maximum leaf mass per ground area for each
316                                                                         !! PFT @tex $(gC m^{-2})$ @endtex 
317!$OMP THREADPRIVATE(lm_lastyearmax)
318  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: lm_thisyearmax       !! This year's maximum leaf mass per ground area for each
319                                                                         !! PFT @tex $(gC m^{-2})$ @endtex 
320!$OMP THREADPRIVATE(lm_thisyearmax)
321  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: maxfpc_lastyear      !! Last year's maximum fpc for each natural PFT, on ground
322                                                                         !! [??CHECK] fpc but this ones look ok (computed in
323                                                                         !! season, used in light)??
324!$OMP THREADPRIVATE(maxfpc_lastyear)
325  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: maxfpc_thisyear      !! This year's maximum fpc for each PFT, on ground (see
326                                                                         !! stomate_season), [??CHECK] fpc but this ones look ok
327                                                                         !! (computed in season, used in light)??
328!$OMP THREADPRIVATE(maxfpc_thisyear)
329!---
330  REAL(r_std), ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: leaf_age             !! Age of different leaf classes (days)
331!$OMP THREADPRIVATE(leaf_age)
332  REAL(r_std), ALLOCATABLE,SAVE,DIMENSION(:,:,:)  :: leaf_frac            !! PFT fraction of leaf mass in leaf age class (0-1,
333                                                                         !! unitless)
334!$OMP THREADPRIVATE(leaf_frac)
335  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: when_growthinit      !! Days since beginning of growing season (days)
336!$OMP THREADPRIVATE(when_growthinit)
337  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: herbivores           !! Time constant of probability of a leaf to be eaten by a
338                                                                         !! herbivore (days)
339!$OMP THREADPRIVATE(herbivores)
340  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: RIP_time             !! How much time ago was the PFT eliminated for the last
341                                                                         !! time (year)
342!$OMP THREADPRIVATE(RIP_time)
343  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: time_hum_min         !! Time elapsed since strongest moisture limitation (days)
344!$OMP THREADPRIVATE(time_hum_min)
345!---
346  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: clay_fm              !! Soil clay content (0-1, unitless), parallel computing
347  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: clay_fm_g            !! Soil clay content (0-1, unitless), parallel computing
348  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: precip_fm            !! Daily precipitations sum @tex $(mm day^{-1})$ @endtex,
349                                                                         !! parallel computing
350  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: precip_fm_g          !! Daily precipitations sum @tex $(mm day^{-1})$ @endtex,
351                                                                         !! parallel computing
352  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: litterhum_daily_fm   !! Daily relative humidity of litter (0-1, unitless),
353                                                                         !! parallel computing
354  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: litterhum_daily_fm_g !! Daily relative humidity of litter (0-1, unitless),
355                                                                         !! parallel computing
356  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: t2m_daily_fm         !! Daily air temperature at 2 meter (K), parallel
357                                                                         !! computing
358  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: t2m_daily_fm_g       !! Daily air temperature at 2 meter (K), parallel
359                                                                         !! computing
360  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: t2m_min_daily_fm     !! Daily minimum air temperature at 2 meter (K),
361                                                                         !! parallel computing
362  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: t2m_min_daily_fm_g   !! Daily minimum air temperature at 2 meter (K),
363                                                                         !! parallel computing
364  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: tsurf_daily_fm       !! Daily surface temperatures (K), parallel
365                                                                         !! computing
366  REAL(r_std),ALLOCATABLE,DIMENSION(:,:)         :: tsurf_daily_fm_g     !! Daily surface temperatures (K), parallel
367                                                                         !! computing
368  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: tsoil_daily_fm       !! Daily soil temperatures (K), parallel computing
369  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: tsoil_daily_fm_g     !! Daily soil temperatures (K), parallel computing
370  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: soilhum_daily_fm     !! Daily soil humidity (0-1, unitless), parallel computing
371  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: soilhum_daily_fm_g   !! Daily soil humidity (0-1, unitless), parallel computing
372  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: humrel_daily_fm      !! Daily relative humidity of atmosphere (0-1, unitless),
373                                                                         !! parallel computing
374  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: humrel_daily_fm_g    !! Daily relative humidity of atmosphere (0-1, unitless),
375                                                                         !! parallel computing
376  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: gpp_daily_fm         !! Daily gross primary productivity per ground area
377                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex,
378                                                                         !! parallel computing
379  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: gpp_daily_fm_g       !! Daily gross primary productivity per ground area
380                                                                         !! @tex $(gC m^{-2} day^{-1})$ @endtex,
381                                                                         !! parallel computing
382  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: veget_fm             !! Vegetation coverage taking into account non-biological
383                                                                         !! coverage (unitless), parallel computing
384  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: veget_fm_g           !! Vegetation coverage taking into account non-biological
385                                                                         !! coverage (unitless), parallel computing
386  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: veget_max_fm         !! Maximum vegetation coverage taking into account
387                                                                         !! non-biological coverage (unitless), parallel computing
388  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: veget_max_fm_g       !! Maximum vegetation coverage taking into account none
389                                                                         !! biological coverage (unitless), parallel computing
390  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: lai_fm               !! Leaf area index @tex $@tex $(m^2 m^{-2})$ @endtex$ @endtex,
391                                                                         !! parallel computing
392  REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)       :: lai_fm_g             !! Leaf area index @tex $@tex $(m^2 m^{-2})$ @endtex$ @endtex,
393                                                                         !! parallel computing
394!---
395  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: co2_fire             !! Carbon emitted to the atmosphere by burning living
396                                                                         !! and dead biomass
397                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
398!$OMP THREADPRIVATE(co2_fire)
399  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: co2_to_bm_dgvm       !! Psuedo-photosynthesis,C used to provide seedlings with
400                                                                         !! an initial biomass, arbitrarily removed from the
401                                                                         !! atmosphere 
402                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
403!$OMP THREADPRIVATE(co2_to_bm_dgvm)
404  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: nep_daily            !! Daily net CO2 flux (positive from atmosphere to land)
405                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
406!$OMP THREADPRIVATE(nep_daily)
407  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: nep_monthly          !! Monthly net CO2 flux (positive from atmosphere to land)
408                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
409!$OMP THREADPRIVATE(nep_monthly)
410  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: prod10               !! Wood products remaining in the 10 year-turnover pool
411                                                                         !! after the annual release for each compartment
412                                                                         !! @tex $(gC m^{-2})$ @endtex   
413                                                                         !! (0:10 input from year of land cover change),
414                                                                         !! dimension(#pixels,0:10 years
415!$OMP THREADPRIVATE(prod10)
416  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: prod100              !! Wood products remaining in the 100 year-turnover pool
417                                                                         !! after the annual release for each compartment
418                                                                         !! @tex $(gC m^{-2})$ @endtex 
419                                                                         !! (0:100 input from year of land cover change),
420                                                                         !! dimension(#pixels,0:100 years)
421!$OMP THREADPRIVATE(prod100)
422  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: flux10               !! Wood decomposition from the 10 year-turnover pool
423                                                                         !! compartments
424                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
425                                                                         !! dimension(#pixels,0:10) 
426!$OMP THREADPRIVATE(flux10)
427  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: flux100              !! Wood decomposition from the 100 year-turnover pool
428                                                                         !! compartments
429                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
430                                                                         !! dimension(#pixels,0:100)
431!$OMP THREADPRIVATE(flux100)
432  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: co2_flux             !! CO2 flux between atmosphere and biosphere
433                                                                         !! @tex $(gC m^{-2} one_day^{-1})$ @endtex
434!$OMP THREADPRIVATE(co2_flux)
435  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: fco2_lu              !! CO2 flux between atmosphere and biosphere from land-use
436                                                                         !! (without forest management)
437                                                                         !! @tex $(gC m^{-2} one_day^{-1})$ @endtex
438!$OMP THREADPRIVATE(fco2_lu)
439  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: fco2_wh              !! CO2 Flux to Atmosphere from Wood Harvesting (positive from atm to land)
440                                                                         !! @tex $(gC m^{-2} one_day^{-1})$ @endtex
441!$OMP THREADPRIVATE(fco2_wh)
442  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: fco2_ha              !! CO2 Flux to Atmosphere from Crop Harvesting (positive from atm to land)
443                                                                         !! @tex $(gC m^{-2} one_day^{-1})$ @endtex
444!$OMP THREADPRIVATE(fco2_ha)
445
446  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: convflux             !! Release during first year following land cover change
447                                                                         !! (paper, burned, etc...)
448                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex 
449!$OMP THREADPRIVATE(convflux)
450  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: cflux_prod10         !! Total annual release from the 10 year-turnover pool
451                                                                         !! sum of flux10 
452                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
453!$OMP THREADPRIVATE(cflux_prod10)
454  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: cflux_prod100        !! Total annual release from the 100 year-turnover pool
455                                                                         !! sum of flux100
456                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
457!$OMP THREADPRIVATE(cflux_prod100)
458  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: prod10_harvest       !! Wood products remaining in the 10 year-turnover pool
459                                                                         !! after the annual release for each compartment
460                                                                         !! @tex $(gC m^{-2})$ @endtex   
461                                                                         !! (0:10 input from year of wood harvest),
462                                                                         !! dimension(#pixels,0:10 years
463!$OMP THREADPRIVATE(prod10_harvest)
464  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: prod100_harvest      !! Wood products remaining in the 100 year-turnover pool
465                                                                         !! after the annual release for each compartment
466                                                                         !! @tex $(gC m^{-2})$ @endtex 
467                                                                         !! (0:100 input from year of wood harvest),
468                                                                         !! dimension(#pixels,0:100 years)
469!$OMP THREADPRIVATE(prod100_harvest)
470  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: flux10_harvest       !! Wood decomposition from the 10 year-turnover pool
471                                                                         !! compartments
472                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
473                                                                         !! dimension(#pixels,0:10) 
474!$OMP THREADPRIVATE(flux10_harvest)
475  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: flux100_harvest      !! Wood decomposition from the 100 year-turnover pool
476                                                                         !! compartments
477                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
478                                                                         !! dimension(#pixels,0:100)
479!$OMP THREADPRIVATE(flux100_harvest)
480  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: convflux_harvest     !! Release during first year following wood harvest
481                                                                         !! (paper, burned, etc...)
482                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex 
483!$OMP THREADPRIVATE(convflux_harvest)
484  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: cflux_prod10_harvest !! Total annual release from the 10 year-turnover pool
485                                                                         !! sum of flux10 
486                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
487!$OMP THREADPRIVATE(cflux_prod10_harvest)
488  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: cflux_prod100_harvest!! Total annual release from the 100 year-turnover pool
489                                                                         !! sum of flux100
490                                                                         !! @tex $(gC m^{-2} year^{-1})$ @endtex
491!$OMP THREADPRIVATE(cflux_prod100_harvest)
492  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: convfluxpft          !! Convflux per PFT                     
493!$OMP THREADPRIVATE(convfluxpft)
494  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: fDeforestToProduct   !! Deforested biomass into product pool due to anthropogenic                                                                                                           
495                                                                         !! land use change                   
496!$OMP THREADPRIVATE(fDeforestToProduct)   
497  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: fLulccResidue        !! Carbon mass flux into soil and litter due to anthropogenic land use or land cover change                                                                         
498!$OMP THREADPRIVATE(fLulccResidue)
499  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)    :: fHarvestToProduct    !! Deforested biomass into product pool due to anthropogenic                                                                                       
500                                                                         !! land use
501!$OMP THREADPRIVATE(fHarvestToProduct)
502  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:):: woodharvestpft       !! New year wood harvest per  PFT
503!$OMP THREADPRIVATE(woodharvestpft)
504  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: harvest_above        !! Harvest of above ground biomass for agriculture -not
505                                                                         !! just from land use change
506                                                                         !! @tex $(??gC m^{-2} dt_stomate^{-1})$ @endtex
507!$OMP THREADPRIVATE(harvest_above)
508  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: carb_mass_total      !! Total on-site and off-site C pool
509                                                                         !! @tex $(??gC m^{-2})$ @endtex                       
510!$OMP THREADPRIVATE(carb_mass_total)
511!---
512  REAL(r_std), SAVE                              :: tau_longterm
513!$OMP THREADPRIVATE(tau_longterm)
514  REAL(r_std),SAVE                               :: dt_days=zero         !! Time step of STOMATE (days)
515!$OMP THREADPRIVATE(dt_days)
516  INTEGER(i_std),SAVE                            :: days_since_beg=0     !! Number of full days done since the start of the simulation
517!$OMP THREADPRIVATE(days_since_beg)
518  INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:)   :: nforce               !! Number of states calculated for the soil forcing
519                                                                         !! variables (unitless), dimension(::nparan*::nbyear) both
520                                                                         !! given in the run definition file   
521!$OMP THREADPRIVATE(nforce)
522  INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:)   :: isf                  !! Index for number of time steps that can be stored in
523                                                                         !! memory (unitless), dimension (#nsfm)
524!$OMP THREADPRIVATE(isf)
525  INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:)   :: nf_cumul             !! Number of years over which the average is calculated in
526                                                                         !! forcesoil when cumul flag is set, dimension (#nsft)
527                                                                         !! [??CHECK] definition the dimension is number of
528                                                                         !! timesteps in a year?
529!$OMP THREADPRIVATE(nf_cumul)
530  INTEGER(i_std), SAVE                           :: spinup_period        !! Period of years used to calculate the resolution of the system for spinup analytic.
531                                                                         !! This period correspond in most cases to the period of years of forcing data used
532!$OMP THREADPRIVATE(spinup_period)
533  INTEGER,PARAMETER                              :: r_typ = nf90_real4   !! Specify data format (server dependent)
534  LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:)          :: nf_written           !! Flag indicating whether the forcing data have been
535                                                                         !! written
536!$OMP THREADPRIVATE(nf_written)
537!---
538  LOGICAL, SAVE                                  :: do_slow=.FALSE.      !! Flag that determines whether stomate_accu calculates
539                                                                         !! the sum(do_slow=.FALSE.) or the mean
540                                                                         !! (do_slow=.TRUE.)
541!$OMP THREADPRIVATE(do_slow)
542  LOGICAL, SAVE                                  :: l_first_stomate = .TRUE.!! Is this the first call of stomate?
543!$OMP THREADPRIVATE(l_first_stomate)
544  LOGICAL, SAVE                                  :: cumul_forcing=.FALSE.!! flag for cumul of forcing if teststomate
545!$OMP THREADPRIVATE(cumul_forcing)
546  LOGICAL, SAVE                                  :: cumul_Cforcing=.FALSE.  !! Flag, if internal parameter cumul_Cforcing is
547                                                                            !! TRUE then ::nbyear (defined in run definition
548                                                                            !! file will be forced to 1 later in this module. If
549                                                                            !! FALSE the mean over ::nbyear is written in forcesoil
550!$OMP THREADPRIVATE(cumul_Cforcing)
551!---   
552  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: harvest_above_monthly   !! [??CHECK] post-processing - should be removed?
553!$OMP THREADPRIVATE(harvest_above_monthly)
554  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:)      :: cflux_prod_monthly      !! [??CHECK] post-processing - should be removed?
555!$OMP THREADPRIVATE(cflux_prod_monthly)
556!---
557  INTEGER(i_std), SAVE                               :: global_years        !! Global counter of years (year)
558!$OMP THREADPRIVATE(global_years)
559  LOGICAL, ALLOCATABLE, SAVE, DIMENSION(:)           :: ok_equilibrium      !! Logical array marking the points where the resolution is ok
560                                                                            !! (true/false)
561!$OMP THREADPRIVATE(ok_equilibrium)
562  LOGICAL, ALLOCATABLE, SAVE, DIMENSION(:)           :: carbon_eq           !! Logical array to mark the carbon pools at equilibrium ?
563                                                                            !! If true, the job stops. (true/false)
564!$OMP THREADPRIVATE(carbon_eq)
565  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)       :: nbp_accu            !! Accumulated Net Biospheric Production over the year (gC.m^2 )
566!$OMP THREADPRIVATE(nbp_accu)
567  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)       :: nbp_flux            !! Net Biospheric Production (gC.m^2.day^{-1})
568!$OMP THREADPRIVATE(nbp_flux)
569  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:)       :: matrixA             !! Matrix containing the fluxes between the carbon pools
570                                                                            !! per sechiba time step
571                                                                            !! @tex $(gC.m^2.day^{-1})$ @endtex
572!$OMP THREADPRIVATE(matrixA)
573  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)         :: vectorB             !! Vector containing the litter increase per sechiba time step
574                                                                            !! @tex $(gC m^{-2})$ @endtex
575!$OMP THREADPRIVATE(vectorB)
576  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: MatrixV             !! Matrix containing the accumulated values of matrixA
577!$OMP THREADPRIVATE(MatrixV)
578  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: VectorU             !! Matrix containing the accumulated values of VectorB
579!$OMP THREADPRIVATE(VectorU)
580  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: MatrixW             !! Matrix containing the opposite of matrixA
581!$OMP THREADPRIVATE(MatrixW)
582  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: previous_stock      !! Array containing the carbon stock calculated by the analytical
583                                                                            !! method in the previous resolution
584!$OMP THREADPRIVATE(previous_stock)
585  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: current_stock       !! Array containing the carbon stock calculated by the analytical
586                                                                            !! method in the current resolution
587!$OMP THREADPRIVATE(current_stock)
588  REAL(r_std), SAVE                                  :: eps_carbon          !! Stopping criterion for carbon pools (unitless,0-1)
589!$OMP THREADPRIVATE(eps_carbon)
590  REAL(r_std),SAVE                                   :: dt_forcesoil        !! Time step of soil forcing file (days)
591!$OMP THREADPRIVATE(dt_forcesoil)
592  INTEGER(i_std),PARAMETER                           :: nparanmax=366       !! Maximum number of time steps per year for forcesoil
593  INTEGER(i_std),SAVE                                :: nparan              !! Number of time steps per year for forcesoil read from run definition (unitless)
594!$OMP THREADPRIVATE(nparan)
595  INTEGER(i_std),SAVE                                :: nbyear=1            !! Number of years saved for forcesoil (unitless)
596!$OMP THREADPRIVATE(nbyear)
597  INTEGER(i_std),SAVE                                :: iatt                !! Time step of forcing of soil processes (iatt = 1 to ::nparan*::nbyear)
598!$OMP THREADPRIVATE(iatt)
599  INTEGER(i_std),SAVE                                :: iatt_old=1          !! Previous ::iatt
600!$OMP THREADPRIVATE(iatt_old)
601  INTEGER(i_std),SAVE                                :: nsfm                !! Number of time steps that can be stored in memory (unitless)
602!$OMP THREADPRIVATE(nsfm)
603  INTEGER(i_std),SAVE                                :: nsft                !! Number of time steps in a year (unitless)
604!$OMP THREADPRIVATE(nsft)
605  INTEGER(i_std),SAVE                                :: iisf                !! Current pointer for teststomate (unitless)
606!$OMP THREADPRIVATE(iisf)
607  CHARACTER(LEN=100), SAVE                           :: forcing_name        !! Name of forcing file 1
608!$OMP THREADPRIVATE(forcing_name)
609  CHARACTER(LEN=100), SAVE                           :: Cforcing_name       !! Name of forcing file 2
610!$OMP THREADPRIVATE(Cforcing_name)
611  INTEGER(i_std),SAVE                                :: Cforcing_id         !! File identifer of file 2
612!$OMP THREADPRIVATE(Cforcing_id)   
613  INTEGER(i_std),PARAMETER                           :: ndm = 10            !! Maximum number of dimensions (unitless)
614
615 
616PUBLIC clay_fm, humrel_daily_fm, litterhum_daily_fm, t2m_daily_fm, &
617   & t2m_min_daily_fm, tsurf_daily_fm, tsoil_daily_fm, soilhum_daily_fm, &
618   & precip_fm, gpp_daily_fm, veget_fm, veget_max_fm, lai_fm
619PUBLIC  dt_days, days_since_beg, do_slow
620PUBLIC isf, nf_written
621
622CONTAINS
623 
624
625!! ================================================================================================================================
626!! SUBROUTINE   : stomate_initialize
627!!
628!>\BRIEF        Initialization routine for stomate module.
629!!
630!! DESCRIPTION  : Initialization routine for stomate module. Read options from parameter file, allocate variables, read variables
631!!                from restart file and initialize variables if necessary.
632!!               
633!! \n
634!_ ================================================================================================================================
635
636SUBROUTINE stomate_initialize &
637        (kjit,           kjpij,             kjpindex,                        &
638         rest_id_stom,   hist_id_stom,      hist_id_stom_IPCC,               &
639         index,          lalo,              neighbours,   resolution,        &
640         contfrac,       totfrac_nobio,     clay,         temp_air,          &
641         lai,            veget,             veget_max,                       &
642         co2_to_bm_radia,deadleaf_cover,    assim_param,  temp_growth )
643
644
645    IMPLICIT NONE
646    !! 0. Variable and parameter declaration
647    !! 0.1 Input variables
648    INTEGER(i_std),INTENT(in)                       :: kjit              !! Time step number (unitless)
649    INTEGER(i_std),INTENT(in)                       :: kjpij             !! Total size of the un-compressed grid (unitless)
650    INTEGER(i_std),INTENT(in)                       :: kjpindex          !! Domain size - terrestrial pixels only (unitless)
651    INTEGER(i_std),INTENT(in)                       :: rest_id_stom      !! STOMATE's _Restart_ file identifier (unitless)
652    INTEGER(i_std),INTENT(in)                       :: hist_id_stom      !! STOMATE's _history_ file identifier (unitless)
653    INTEGER(i_std),INTENT(in)                       :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file identifier(unitless)
654    INTEGER(i_std),DIMENSION(kjpindex),INTENT(in)   :: index             !! The indices of the terrestrial pixels only (unitless)
655    REAL(r_std),DIMENSION(kjpindex,2),INTENT(in)    :: lalo              !! Geographical coordinates (latitude,longitude) for pixels (degrees)
656    INTEGER(i_std),DIMENSION(kjpindex,NbNeighb),INTENT(in) :: neighbours !! Neighoring grid points if land for the DGVM (unitless)
657    REAL(r_std),DIMENSION(kjpindex,2),INTENT(in)    :: resolution        !! Size in x an y of the grid (m) - surface area of the gridbox
658    REAL(r_std),DIMENSION (kjpindex), INTENT (in)   :: contfrac          !! Fraction of continent in the grid cell (unitless)
659    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: totfrac_nobio     !! Fraction of grid cell covered by lakes, land ice, cities, ... (unitless)
660    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: clay              !! Clay fraction of soil (0-1, unitless)
661    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: temp_air          !! Air temperature at first atmospheric model layer (K)
662    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: lai               !! Leaf area inex @tex $(m^2 m^{-2})$ @endtex
663    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: veget             !! Fraction of vegetation type including
664                                                                         !! non-biological fraction (unitless)
665    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: veget_max         !! Maximum fraction of vegetation type including
666                                                                         !! non-biological fraction (unitless)
667
668    !! 0.2 Output variables
669
670    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: co2_to_bm_radia   !! virtual gpp flux between atmosphere and biosphere
671    REAL(r_std),DIMENSION(kjpindex),INTENT(out)     :: deadleaf_cover    !! Fraction of soil covered by dead leaves (unitless)
672    REAL(r_std),DIMENSION(kjpindex,nvm,npco2),INTENT(out) :: assim_param !! min+max+opt temperatures (K) & vmax for photosynthesis 
673                                                                         !! @tex $(\mu mol m^{-2}s^{-1})$ @endtex 
674    REAL(r_std),DIMENSION(kjpindex),INTENT(out)     :: temp_growth       !! Growth temperature (°C) 
675                                                                         !! Is equal to t2m_month
676    !! 0.3 Local variables
677    REAL(r_std)                                   :: dt_days_read             !! STOMATE time step read in restart file (days)
678    INTEGER(i_std)                                :: l,k,ji, jv, i, j, m      !! indices   
679    REAL(r_std),PARAMETER                         :: max_dt_days = 5.         !! Maximum STOMATE time step (days)
680    REAL(r_std),DIMENSION(kjpindex,nvm)           :: rprof                    !! Coefficient of the exponential functions that
681                                                                              !! relates root density to soil depth (unitless)
682    REAL(r_std),DIMENSION(kjpindex,nvm)           :: gpp_daily_x              !! "Daily" gpp for teststomate 
683                                                                              !! @tex $(??gC m^{-2} dt_stomate^{-1})$ @endtex
684    REAL(r_std),DIMENSION(kjpindex,nvm)           :: veget_cov                !! Fractional coverage: actually share of the pixel
685                                                                              !! covered by a PFT (fraction of ground area),
686                                                                              !! taking into account LAI ??(= grid scale fpc)??
687    INTEGER(i_std)                                :: ier                      !! Check errors in netcdf call (unitless)
688
689    INTEGER(i_std)                                :: max_totsize              !! Memory management - maximum memory size (Mb)
690    INTEGER(i_std)                                :: totsize_1step            !! Memory management - memory required to store one
691                                                                              !! time step on one processor (Mb)
692    INTEGER(i_std)                                :: totsize_tmp              !! Memory management - memory required to store one
693                                                                              !! time step on all processors(Mb)
694    INTEGER(i_std)                                :: vid                      !! Variable identifer of netCDF (unitless)
695    INTEGER(i_std)                                :: nneigh                   !! Number of neighbouring pixels
696    INTEGER(i_std)                                :: direct                   !!
697    INTEGER(i_std),DIMENSION(ndm)                 :: d_id                     !!
698
699
700!_ ================================================================================================================================
701   
702    !! 1. Initialize variable
703    !! Update flag
704    l_first_stomate = .FALSE.
705   
706    !! 1.1 Store current time step in a common variable
707    itime = kjit
708   
709   
710    !! 1.3 PFT rooting depth across pixels, humescte is pre-defined
711    ! (constantes_veg.f90). It is defined as the coefficient of an exponential
712    ! function relating root density to depth
713    DO j=1,nvm
714       rprof(:,j) = 1./humcste(j)
715    ENDDO
716   
717    !! 1.4.0 Parameters for spinup
718    !
719    eps_carbon = 0.01
720    !Config Key   = EPS_CARBON
721    !Config Desc  = Allowed error on carbon stock
722    !Config If    = SPINUP_ANALYTIC
723    !Config Def   = 0.01
724    !Config Help  =
725    !Config Units = [%]   
726    CALL getin_p('EPS_CARBON',eps_carbon)       
727   
728   
729    !Config Key   = SPINUP_PERIOD
730    !Config Desc  = Period to calulcate equilibrium during spinup analytic
731    !Config If    = SPINUP_ANALYTIC
732    !Config Def   = -1
733    !Config Help  = Period corresponds in most cases to the number of years of forcing data used in the spinup.
734    !Config Units = [years]   
735    spinup_period = -1
736    CALL getin_p('SPINUP_PERIOD',spinup_period)       
737   
738    ! Check spinup_period values.
739    ! For periods uptil 6 years, to obtain equilibrium, a bigger period have to be used
740    ! and therefore spinup_period is adjusted to 10 years.
741    IF (spinup_analytic) THEN
742       IF (spinup_period <= 0) THEN
743          WRITE(numout,*) 'Error in parameter spinup_period. This parameter must be > 0 : spinup_period=',spinup_period
744          CALL ipslerr_p (3,'stomate_initialize', &
745               'Parameter spinup_period must be set to a positive integer.', &
746               'Set this parameter to the number of years of forcing data used for the spinup.', &
747               '')
748       ELSE IF (spinup_period <= 6) THEN
749          ! Adjust to bigger period. The period must be a multiple of the original period.
750          WRITE(numout,*) 'Initial spinup_period =',spinup_period,' will be adjusted.'
751          spinup_period = spinup_period*(INT(9/spinup_period)+1)
752       END IF
753       IF (printlev >=1) WRITE(numout,*) 'Spinup analytic is activated using eps_carbon=',&
754            eps_carbon, ' and spinup_period=',spinup_period
755    END IF
756   
757
758    !! 1.4.1 Allocate memory for all variables in stomate
759    ! Allocate memory for all variables in stomate, build new index
760    ! tables accounting for the PFTs, read and check flags and set file
761    ! identifier for restart and history files.
762    CALL stomate_init (kjpij, kjpindex, index, lalo, &
763         rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
764   
765    !! 1.4.2 Initialization of PFT specific parameters
766    ! Initialization of PFT specific parameters i.e. sla from leaf life,
767    ! sapling characteristics (biomass), migration speed, critical diameter,
768    ! coldest tolerable temperature, critical values for phenology, maximum
769    ! life time of leaves, respiration coefficients and photosynthesis.
770    ! The subroutine also communicates settings read by stomate_constant_init.
771    CALL data (kjpindex, lalo)
772   
773    !! 1.4.3 Initial conditions
774   
775    !! 1.4.3.1 Read initial values for STOMATE's variables from the _restart_ file
776
777    ! Get values from _restart_ file. Note that only ::kjpindex, ::index, ::lalo
778    ! and ::resolution are input variables, all others are output variables.
779    CALL readstart &
780         (kjpindex, index, lalo, resolution, temp_air, &
781         dt_days_read, days_since_beg, &
782         ind, adapted, regenerate, &
783         humrel_daily, gdd_init_date, litterhum_daily, &
784         t2m_daily, t2m_min_daily, tsurf_daily, tsoil_daily, &
785         soilhum_daily, precip_daily, &
786         gpp_daily, npp_daily, turnover_daily, &
787         humrel_month, humrel_week, &
788         t2m_longterm, tau_longterm, t2m_month, t2m_week, &
789         tsoil_month, soilhum_month, fireindex, firelitter, &
790         maxhumrel_lastyear, maxhumrel_thisyear, &
791         minhumrel_lastyear, minhumrel_thisyear, &
792         maxgppweek_lastyear, maxgppweek_thisyear, &
793         gdd0_lastyear, gdd0_thisyear, &
794         precip_lastyear, precip_thisyear, &
795         gdd_m5_dormance,  gdd_from_growthinit, gdd_midwinter, ncd_dormance, ngd_minus5, &
796         PFTpresent, npp_longterm, lm_lastyearmax, lm_thisyearmax, &
797         maxfpc_lastyear, maxfpc_thisyear, &
798         turnover_longterm, gpp_week, biomass, resp_maint_part, &
799         leaf_age, leaf_frac, &
800         senescence, when_growthinit, age, &
801         resp_hetero_d, resp_maint_d, resp_growth_d, co2_fire, co2_to_bm_dgvm, co2_to_bm_radia, &
802         veget_lastlight, everywhere, need_adjacent, RIP_time, &
803         time_hum_min, hum_min_dormance, &
804         litterpart, litter, dead_leaves, &
805         carbon, lignin_struc,turnover_time,&
806         co2_flux, fco2_lu, fco2_wh, fco2_ha, &
807         prod10,prod100,flux10, flux100, &
808         convflux, cflux_prod10, cflux_prod100, &
809         prod10_harvest,prod100_harvest,flux10_harvest, flux100_harvest, &
810         convflux_harvest, cflux_prod10_harvest, cflux_prod100_harvest, &
811         convfluxpft, fDeforestToProduct, fLulccResidue,fHarvestToProduct, &
812         woodharvestpft, bm_to_litter, carb_mass_total, &
813         Tseason, Tseason_length, Tseason_tmp, &
814         Tmin_spring_time, begin_leaves, onset_date, &
815         global_years, ok_equilibrium, nbp_accu, nbp_flux, &
816         MatrixV, VectorU, previous_stock, current_stock, assim_param)
817   
818    !! 1.4.5 Check time step
819       
820    !! 1.4.5.1 Allow STOMATE's time step to change although this is dangerous
821    IF (dt_days /= dt_days_read) THEN
822       WRITE(numout,*) 'slow_processes: STOMATE time step changes:', &
823            & dt_days_read,' -> ',dt_days
824    ENDIF
825   
826    !! 1.4.5.2 Time step has to be a multiple of a full day
827    IF ( ( dt_days-REAL(NINT(dt_days),r_std) ) > min_stomate ) THEN
828       WRITE(numout,*) 'slow_processes: STOMATE time step is not a mutiple of a full day:', &
829            & dt_days,' days.'
830       STOP
831    ENDIF
832   
833    !! 1.4.5.3 upper limit to STOMATE's time step
834    IF ( dt_days > max_dt_days ) THEN
835       WRITE(numout,*) 'slow_processes: STOMATE time step exceeds the maximum value:', &
836            & dt_days,' days > ', max_dt_days, ' days.' 
837       STOP
838    ENDIF
839   
840    !! 1.4.5.4 STOMATE time step must not be less than the forcing time step
841    IF ( dt_sechiba > dt_days*one_day ) THEN
842       WRITE(numout,*) &
843            & 'slow_processes: STOMATE time step ::dt_days smaller than forcing time step ::dt_sechiba'
844       STOP
845    ENDIF
846   
847    !! 1.4.5.6 Final message on time step
848    IF (printlev >=2) WRITE(numout,*) 'Slow_processes, STOMATE time step (days): ', dt_days
849   
850    !! 1.4.6 Write forcing file for teststomate
851    IF (allow_forcing_write) THEN
852       
853       !Config Key   = STOMATE_FORCING_NAME
854       !Config Desc  = Name of STOMATE's forcing file
855       !Config If    = OK_STOMATE
856       !Config Def   = NONE
857       !Config Help  = Name that will be given
858       !Config         to STOMATE's offline forcing file
859       !Config         Compatible with Nicolas Viovy's driver
860       !Config Units = [FILE]
861       forcing_name = stomate_forcing_name
862       CALL getin_p('STOMATE_FORCING_NAME',forcing_name)
863       
864       IF ( TRIM(forcing_name) /= 'NONE' ) THEN
865         
866          !! 1.4.6.1 Calculate steps that can be stored in memory
867          ! Action for the root processor only (parallel computing) 
868          IF (is_root_prc) CALL SYSTEM ('rm -f '//TRIM(forcing_name))
869          IF (printlev>=2) WRITE(numout,*) 'writing a forcing file for STOMATE.'
870         
871          !Config Key   = STOMATE_FORCING_MEMSIZE
872          !Config Desc  = Size of STOMATE forcing data in memory
873          !Config If    = OK_STOMATE
874          !Config Def   = 50
875          !Config Help  = This variable determines how many
876          !Config         forcing states will be kept in memory.
877          !Config         Must be a compromise between memory
878          !Config         use and frequeny of disk access.
879          !Config Units = [MegaBytes]
880          max_totsize = 50
881          CALL getin_p('STOMATE_FORCING_MEMSIZE', max_totsize)     
882          max_totsize = max_totsize*1000000
883         
884          totsize_1step = &
885               SIZE(clay)*KIND(clay) &
886               +SIZE(humrel_daily)*KIND(humrel_daily) &
887               +SIZE(litterhum_daily)*KIND(litterhum_daily) &
888               +SIZE(t2m_daily)*KIND(t2m_daily) &
889               +SIZE(t2m_min_daily)*KIND(t2m_min_daily) &
890               +SIZE(tsurf_daily)*KIND(tsurf_daily) &
891               +SIZE(tsoil_daily)*KIND(tsoil_daily) &
892               +SIZE(soilhum_daily)*KIND(soilhum_daily) &
893               +SIZE(precip_daily)*KIND(precip_daily) &
894               +SIZE(gpp_daily_x)*KIND(gpp_daily_x) &
895               +SIZE(veget)*KIND(veget) &
896               +SIZE(veget_max)*KIND(veget_max) &
897               +SIZE(lai)*KIND(lai)
898         
899          ! Totsize_1step is the size on a single processor, sum
900          ! all processors and send to all processors
901          CALL reduce_sum(totsize_1step,totsize_tmp)
902          CALL bcast(totsize_tmp)
903          totsize_1step=totsize_tmp
904         
905          ! Total number of forcing steps
906          nsft = INT(one_year/(dt_stomate/one_day))
907         
908          ! Number of forcing steps in memory
909          nsfm = MIN(nsft, &
910               MAX(1,NINT( REAL(max_totsize,r_std) &
911               /REAL(totsize_1step,r_std))))
912           
913             
914          !! 1.6.4.2 Allocate memory for variables containing forcing data 
915          ! and initialize variables (set to zero).
916          CALL init_forcing (kjpindex,nsfm,nsft)
917         
918          ! Indexing for writing forcing file
919          isf(:) = (/ (i,i=1,nsfm) /)
920          nf_written(:) = .FALSE.
921          nf_cumul(:) = 0
922          iisf = 0
923         
924          !! 1.6.4.3 Create netcdf file
925          ! Create, define and populate a netcdf file containing the forcing data.
926          ! For the root processor only (parallel computing). NF90_ are functions
927          ! from and external library. 
928          IF (is_root_prc) THEN
929             
930             ! Create new netCDF dataset
931             ier = NF90_CREATE (TRIM(forcing_name),NF90_SHARE,forcing_id)
932             
933             ! Add variable attribute
934             ! Note ::iim_g and ::jjm_g are dimensions of the global field and
935             ! ::nbp_glo is the number of global continental points
936             ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL,'dt_sechiba',dt_sechiba)
937             ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL,'dt_stomate',dt_stomate)
938             ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL, &
939                  'nsft',REAL(nsft,r_std))
940             ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL, &
941                  'kjpij',REAL(iim_g*jjm_g,r_std))
942             ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL, &
943                  'kjpindex',REAL(nbp_glo,r_std))
944             
945             ! Add new dimension
946             ier = NF90_DEF_DIM (forcing_id,'points',nbp_glo,d_id(1))
947             ier = NF90_DEF_DIM (forcing_id,'layers',nslm,d_id(2))
948             ier = NF90_DEF_DIM (forcing_id,'pft',nvm,d_id(3))
949             direct=2
950             ier = NF90_DEF_DIM (forcing_id,'direction',direct,d_id(4))
951             nneigh=8
952             ier = NF90_DEF_DIM (forcing_id,'nneigh',nneigh,d_id(5))
953             ier = NF90_DEF_DIM (forcing_id,'time',NF90_UNLIMITED,d_id(6))
954             ier = NF90_DEF_DIM (forcing_id,'nbparts',nparts,d_id(7))
955             
956             ! Add new variable
957             ier = NF90_DEF_VAR (forcing_id,'points',    r_typ,d_id(1),vid)
958             ier = NF90_DEF_VAR (forcing_id,'layers',    r_typ,d_id(2),vid)
959             ier = NF90_DEF_VAR (forcing_id,'pft',       r_typ,d_id(3),vid)
960             ier = NF90_DEF_VAR (forcing_id,'direction', r_typ,d_id(4),vid)
961             ier = NF90_DEF_VAR (forcing_id,'nneigh',    r_typ,d_id(5),vid)
962             ier = NF90_DEF_VAR (forcing_id,'time',      r_typ,d_id(6),vid)
963             ier = NF90_DEF_VAR (forcing_id,'nbparts',   r_typ,d_id(7),vid)
964             ier = NF90_DEF_VAR (forcing_id,'index',     r_typ,d_id(1),vid)
965             ier = NF90_DEF_VAR (forcing_id,'contfrac',  r_typ,d_id(1),vid) 
966             ier = NF90_DEF_VAR (forcing_id,'lalo', &
967                  r_typ,(/ d_id(1),d_id(4) /),vid)
968             ier = NF90_DEF_VAR (forcing_id,'neighbours', &
969                  r_typ,(/ d_id(1),d_id(5) /),vid)
970             ier = NF90_DEF_VAR (forcing_id,'resolution', &
971                  r_typ,(/ d_id(1),d_id(4) /),vid)
972             ier = NF90_DEF_VAR (forcing_id,'clay', &
973                  r_typ,(/ d_id(1),d_id(6) /),vid)
974             ier = NF90_DEF_VAR (forcing_id,'humrel', &
975                  r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
976             ier = NF90_DEF_VAR (forcing_id,'litterhum', &
977                  r_typ,(/ d_id(1),d_id(6) /),vid)
978             ier = NF90_DEF_VAR (forcing_id,'t2m', &
979                  r_typ,(/ d_id(1),d_id(6) /),vid)
980             ier = NF90_DEF_VAR (forcing_id,'t2m_min', &
981                  r_typ,(/ d_id(1),d_id(6) /),vid)
982             ier = NF90_DEF_VAR (forcing_id,'tsurf', &
983                  r_typ,(/ d_id(1),d_id(6) /),vid)
984             ier = NF90_DEF_VAR (forcing_id,'tsoil', &
985                  r_typ,(/ d_id(1),d_id(2),d_id(6) /),vid)
986             ier = NF90_DEF_VAR (forcing_id,'soilhum', &
987                  r_typ,(/ d_id(1),d_id(2),d_id(6) /),vid)
988             ier = NF90_DEF_VAR (forcing_id,'precip', &
989                  r_typ,(/ d_id(1),d_id(6) /),vid)
990             ier = NF90_DEF_VAR (forcing_id,'gpp', &
991                  r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
992             ier = NF90_DEF_VAR (forcing_id,'veget', &
993                  r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
994             ier = NF90_DEF_VAR (forcing_id,'veget_max', &
995                  r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
996             ier = NF90_DEF_VAR (forcing_id,'lai', &
997                  r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
998             ier = NF90_ENDDEF (forcing_id)
999             
1000             ! Given the name of a varaible, nf90_inq_varid finds the variable
1001             ! ID (::vid). Put data value(s) into variable ::vid
1002             ier = NF90_INQ_VARID (forcing_id,'points',vid)
1003             ier = NF90_PUT_VAR (forcing_id,vid, &
1004                  (/(REAL(i,r_std),i=1,nbp_glo) /))
1005             ier = NF90_INQ_VARID (forcing_id,'layers',vid)
1006             ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,nslm)/))
1007             ier = NF90_INQ_VARID (forcing_id,'pft',vid)
1008             ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,nvm)/))
1009             ier = NF90_INQ_VARID (forcing_id,'direction',vid)
1010             ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,2)/))
1011             ier = NF90_INQ_VARID (forcing_id,'nneigh',vid)
1012             ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,8)/))
1013             ier = NF90_INQ_VARID (forcing_id,'time',vid)
1014             ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,nsft)/))
1015             ier = NF90_INQ_VARID (forcing_id,'nbparts',vid)
1016             ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,nparts)/))
1017             ier = NF90_INQ_VARID (forcing_id,'index',vid) 
1018             ier = NF90_PUT_VAR (forcing_id,vid,REAL(index_g,r_std))
1019             ier = NF90_INQ_VARID (forcing_id,'contfrac',vid)
1020             ier = NF90_PUT_VAR (forcing_id,vid,REAL(contfrac_g,r_std))
1021             ier = NF90_INQ_VARID (forcing_id,'lalo',vid)
1022             ier = NF90_PUT_VAR (forcing_id,vid,lalo_g)
1023             !ym attention a neighbours, a modifier plus tard     
1024             ier = NF90_INQ_VARID (forcing_id,'neighbours',vid)
1025             ier = NF90_PUT_VAR (forcing_id,vid,REAL(neighbours_g,r_std))
1026             ier = NF90_INQ_VARID (forcing_id,'resolution',vid)
1027             ier = NF90_PUT_VAR (forcing_id,vid,resolution_g)
1028          ENDIF ! is_root_prc
1029       ENDIF ! (forcing_name) /= 'NONE'
1030    ENDIF 
1031   
1032    !! 1.4.7 write forcing file for forcesoil
1033    !! 1.4.7.1 Initialize
1034    !Config Key   = STOMATE_CFORCING_NAME
1035    !Config Desc  = Name of STOMATE's carbon forcing file
1036    !Config If    = OK_STOMATE
1037    !Config Def   = NONE
1038    !Config Help  = Name that will be given to STOMATE's carbon
1039    !Config         offline forcing file
1040    !Config         Compatible with Nicolas Viovy's driver
1041    !Config Units = [FILE]
1042    Cforcing_name = stomate_Cforcing_name
1043    CALL getin_p('STOMATE_CFORCING_NAME',Cforcing_name)
1044   
1045    IF ( TRIM(Cforcing_name) /= 'NONE' ) THEN
1046       
1047       ! Time step of forcesoil
1048       !Config Key   = FORCESOIL_STEP_PER_YEAR
1049       !Config Desc  = Number of time steps per year for carbon spinup.
1050       !Config If    = OK_STOMATE
1051       !Config Def   = 365
1052       !Config Help  = Number of time steps per year for carbon spinup.
1053       !Config Units = [days, months, year]
1054       nparan = 365
1055       CALL getin_p('FORCESOIL_STEP_PER_YEAR', nparan)
1056       
1057       ! Correct if setting is out of bounds
1058       IF ( nparan < 1 ) nparan = 1
1059       
1060       !Config Key   = FORCESOIL_NB_YEAR
1061       !Config Desc  = Number of years saved for carbon spinup.
1062       !Config If    = OK_STOMATE
1063       !Config Def   = 1
1064       !Config Help  = Number of years saved for carbon spinup. If internal parameter cumul_Cforcing is TRUE in stomate.f90
1065       !Config         Then this parameter is forced to one.
1066       !Config Units = [years]
1067       CALL getin_p('FORCESOIL_NB_YEAR', nbyear)
1068       
1069       ! Set ::nbyear to 1. if ::cumul_Cforcing=.TRUE.
1070       IF ( cumul_Cforcing ) THEN
1071          CALL ipslerr_p (1,'stomate', &
1072               'Internal parameter cumul_Cforcing is TRUE in stomate.f90', &
1073               'Parameter FORCESOIL_NB_YEAR is therefore forced to 1.', &
1074               '::nbyear is thus set to 1.')
1075          nbyear=1
1076       ENDIF
1077       
1078       ! Make use of ::nparan to calculate ::dt_forcesoil
1079       dt_forcesoil = zero
1080       nparan = nparan+1
1081       DO WHILE ( dt_forcesoil < dt_stomate/one_day )
1082          nparan = nparan-1
1083          IF ( nparan < 1 ) THEN
1084             STOP 'Problem with number of soil forcing time steps ::nparan < 1.'
1085          ENDIF
1086          dt_forcesoil = one_year/REAL(nparan,r_std)
1087       ENDDO
1088       IF ( nparan > nparanmax ) THEN
1089          STOP 'Problem with number of soil forcing time steps ::nparan > ::nparanmax'
1090       ENDIF
1091       IF (printlev>=2) WRITE(numout,*) 'Time step of soil forcing (d): ',dt_forcesoil
1092       
1093       ! Allocate memory for the forcing variables of soil dynamics
1094       ALLOCATE( nforce(nparan*nbyear))
1095       nforce(:) = 0
1096       ALLOCATE(control_moist(kjpindex,nlevs,nparan*nbyear))
1097       ALLOCATE(npp_equil(kjpindex,nparan*nbyear))
1098       ALLOCATE(npp_tot(kjpindex))
1099       ALLOCATE(control_temp(kjpindex,nlevs,nparan*nbyear))
1100       ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan*nbyear)) 
1101       
1102       ! Initialize variables, set to zero
1103       control_moist(:,:,:) = zero
1104       npp_equil(:,:) = zero
1105       npp_tot(:) = zero
1106       control_temp(:,:,:) = zero
1107       soilcarbon_input(:,:,:,:) = zero
1108       
1109    ENDIF ! Cforcing_name) /= 'NONE'
1110   
1111    !! 1.4.8 Calculate STOMATE's vegetation fractions from veget, veget_max
1112    DO j=1,nvm
1113       WHERE ((1.-totfrac_nobio(:)) > min_sechiba)       
1114          ! Pixels with vegetation
1115          veget_cov(:,j) = veget(:,j)/( 1.-totfrac_nobio(:) )
1116          veget_cov_max(:,j) = veget_max(:,j)/( 1.-totfrac_nobio(:) )
1117       ELSEWHERE
1118          ! Pixels without vegetation
1119          veget_cov(:,j) = zero
1120          veget_cov_max(:,j) = zero
1121       ENDWHERE
1122    ENDDO ! Loop over PFTs
1123
1124    !! 1.4.9 Initialize non-zero variables
1125    CALL stomate_var_init &
1126         (kjpindex, veget_cov_max, leaf_age, leaf_frac, &
1127         dead_leaves, &
1128         veget, lai, deadleaf_cover, assim_param)
1129   
1130    ! Initialize land cover change variable
1131    ! ??Should be integrated in the subroutine??
1132    harvest_above(:) = zero
1133   
1134    ! Initialize temp_growth
1135    temp_growth(:)=t2m_month(:)-tp_00 
1136
1137     
1138  END SUBROUTINE stomate_initialize
1139 
1140
1141!! ================================================================================================================================
1142!! SUBROUTINE   : stomate_main
1143!!
1144!>\BRIEF        Manages variable initialisation, reading and writing forcing
1145!! files, aggregating data at stomate's time step (dt_stomate), aggregating data
1146!! at longer time scale (i.e. for phenology) and uses these forcing to calculate
1147!! CO2 fluxes (NPP and respirations) and C-pools (litter, soil, biomass, ...)
1148!!
1149!! DESCRIPTION  : The subroutine manages
1150!! divers tasks:
1151!! (1) Initializing all variables of stomate (first call)
1152!! (2) Reading and writing forcing data (last call)
1153!! (3) Adding CO2 fluxes to the IPCC history files
1154!! (4) Converting the time steps of variables to maintain consistency between
1155!! sechiba and stomate
1156!! (5) Use these variables to call stomate_lpj, maint_respiration, littercalc,
1157!! soilcarbon. The called subroutines handle: climate constraints
1158!! for PFTs, PFT dynamics, Phenology, Allocation, NPP (based on GPP and
1159!! authothropic respiration), fire, mortality, vmax, assimilation temperatures,
1160!! all turnover processes, light competition, sapling establishment, lai, 
1161!! land cover change and litter and soil dynamics.
1162!! (6) Use the spin-up method developed by Lardy (2011)(only if SPINUP_ANALYTIC
1163!! is set to TRUE).
1164!!
1165!! RECENT CHANGE(S) : None
1166!!
1167!! MAIN OUTPUT VARIABLE(S): deadleaf_cover, assim_param, lai, height, veget,
1168!! veget_max, resp_maint, resp_hetero, resp_growth,
1169!! co2_flux_out, fco2_lu_out, fco2_wh_out, fco2_ha_out.
1170!!
1171!! REFERENCES   :
1172!! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium,
1173!! Environmental Modelling & Software (2011), doi:10.1016|j.envsoft.2011.05.016
1174!!
1175!! FLOWCHART    :
1176!! \latexonly
1177!! \includegraphics[scale=0.5]{stomatemainflow.png}
1178!! \endlatexonly
1179!! \n
1180!_ ================================================================================================================================
1181 
1182SUBROUTINE stomate_main &
1183       & (kjit, kjpij, kjpindex, &
1184       &  index, lalo, neighbours, resolution, contfrac, totfrac_nobio, clay, &
1185       &  temp_air, temp_sol, stempdiag, &
1186       &  humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, &
1187       &  gpp, deadleaf_cover, assim_param, &
1188       &  lai, frac_age, height, veget, veget_max, &
1189       &  veget_max_new, woodharvest, totfrac_nobio_new, fraclut, &
1190       &  rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
1191       &  co2_flux_out, fco2_lu_out, fco2_wh_out, fco2_ha_out, &
1192       &  resp_maint, resp_hetero, resp_growth, co2_to_bm_radia, temp_growth)
1193   
1194    IMPLICIT NONE
1195
1196   
1197  !! 0. Variable and parameter declaration
1198
1199    !! 0.1 Input variables
1200
1201    INTEGER(i_std),INTENT(in)                       :: kjit              !! Time step number (unitless)
1202    INTEGER(i_std),INTENT(in)                       :: kjpindex          !! Domain size - terrestrial pixels only (unitless)
1203    INTEGER(i_std),INTENT(in)                       :: kjpij             !! Total size of the un-compressed grid (unitless)
1204    INTEGER(i_std),INTENT(in)                       :: rest_id_stom      !! STOMATE's _Restart_ file identifier (unitless)
1205    INTEGER(i_std),INTENT(in)                       :: hist_id_stom      !! STOMATE's _history_ file identifier (unitless)
1206    INTEGER(i_std),INTENT(in)                       :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file identifier
1207                                                                         !! (unitless)
1208    INTEGER(i_std),DIMENSION(kjpindex),INTENT(in)   :: index             !! Indices of the pixels on the map. Stomate uses a
1209                                                                         !! reduced grid excluding oceans. ::index contains
1210                                                                         !! the indices of the terrestrial pixels only
1211                                                                         !! (unitless)
1212    INTEGER(i_std),DIMENSION(kjpindex,NbNeighb),INTENT(in) :: neighbours !! Neighoring grid points if land for the DGVM
1213                                                                         !! (unitless)
1214    REAL(r_std),DIMENSION(kjpindex,2),INTENT(in)    :: lalo              !! Geographical coordinates (latitude,longitude)
1215                                                                         !! for pixels (degrees)
1216    REAL(r_std),DIMENSION(kjpindex,2),INTENT(in)    :: resolution        !! Size in x an y of the grid (m) - surface area of
1217                                                                         !! the gridbox
1218    REAL(r_std),DIMENSION (kjpindex), INTENT (in)   :: contfrac          !! Fraction of continent in the grid cell (unitless)
1219    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: totfrac_nobio     !! Fraction of grid cell covered by lakes, land
1220                                                                         !! ice, cities, ... (unitless)
1221    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: clay              !! Clay fraction of soil (0-1, unitless)
1222    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: humrel            !! Relative humidity ("moisture availability")
1223                                                                         !! (0-1, unitless)
1224    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: temp_air          !! Air temperature at first atmosperic model layer (K)
1225    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: temp_sol          !! Surface temperature (K)
1226    REAL(r_std),DIMENSION(kjpindex,nslm),INTENT(in) :: stempdiag         !! Soil temperature (K)
1227    REAL(r_std),DIMENSION(kjpindex,nslm),INTENT(in) :: shumdiag          !! Relative soil moisture (0-1, unitless)
1228    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: litterhumdiag     !! Litter humidity (0-1, unitless)
1229    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: precip_rain       !! Rain precipitation 
1230                                                                         !! @tex $(mm dt_stomate^{-1})$ @endtex
1231    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: precip_snow       !! Snow precipitation 
1232                                                                         !! @tex $(mm dt_stomate^{-1})$ @endtex
1233    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: gpp               !! GPP of total ground area 
1234                                                                         !! @tex $(gC m^{-2} time step^{-1})$ @endtex
1235                                                                         !! Calculated in sechiba, account for vegetation
1236                                                                         !! cover and effective time step to obtain ::gpp_d
1237    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: veget_max_new     !! New "maximal" coverage fraction of a PFT: only if
1238                                                                         !! vegetation is updated in slowproc
1239    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: woodharvest       !! Harvested wood biomass (gC m-2 yr-1)
1240    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: totfrac_nobio_new !! New fraction of nobio per gridcell
1241    REAL(r_std),DIMENSION(kjpindex, nlut),INTENT(in):: fraclut           !! Fraction of landuse tiles
1242
1243    !! 0.2 Output variables
1244
1245    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: co2_flux_out      !! CO2 flux between atmosphere and biosphere per
1246                                                                         !! average ground area
1247                                                                         !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex 
1248    REAL(r_std),DIMENSION(kjpindex),INTENT(out)     :: fco2_lu_out       !! CO2 flux between atmosphere and biosphere from
1249                                                                         !! land-use (without forest management) (gC/m2/dt_stomate)
1250    REAL(r_std),DIMENSION(kjpindex),INTENT(out)     :: fco2_wh_out       !! CO2 Flux to Atmosphere from Wood Harvesting (gC/m2/dt_stomate)
1251    REAL(r_std),DIMENSION(kjpindex),INTENT(out)     :: fco2_ha_out       !! CO2 Flux to Atmosphere from Crop Harvesting (gC/m2/dt_stomate)
1252    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: resp_maint        !! Maitenance component of autotrophic respiration in
1253                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex
1254    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: resp_growth       !! Growth component of autotrophic respiration in
1255                                                                         !! @tex ($gC m^{-2} dt_stomate^{-1}$) @endtex
1256    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: resp_hetero       !! Heterotrophic respiration in 
1257                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex 
1258    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: co2_to_bm_radia   !! Virtual gpp created for equilibrium of carbon mass 
1259                                                                         !! @tex $(gC m^{-2} dt_stomate^{-1})$ @endtex 
1260    REAL(r_std),DIMENSION(kjpindex),INTENT(out)     :: temp_growth       !! Growth temperature (°C) 
1261                                                                         !! Is equal to t2m_month
1262
1263    !! 0.3 Modified
1264   
1265    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(inout)       :: lai            !! Leaf area inex @tex $(m^2 m^{-2})$ @endtex
1266    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)          :: veget          !! Fraction of vegetation type including
1267                                                                              !! non-biological fraction (unitless)
1268    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(inout)       :: veget_max      !! Maximum fraction of vegetation type including
1269                                                                              !! non-biological fraction (unitless)
1270    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(inout)       :: height         !! Height of vegetation (m)
1271    REAL(r_std),DIMENSION(kjpindex,nvm,npco2),INTENT(inout) :: assim_param    !! min+max+opt temperatures (K) & vmax for
1272                                                                              !! photosynthesis 
1273                                                                              !! @tex $(\mu mol m^{-2}s^{-1})$ @endtex 
1274    REAL(r_std),DIMENSION(kjpindex),INTENT(inout)           :: deadleaf_cover !! Fraction of soil covered by dead leaves
1275                                                                              !! (unitless)
1276    REAL(r_std),DIMENSION(kjpindex,nvm,nleafages),INTENT(inout):: frac_age    !! Age efficacity from STOMATE
1277
1278    !! 0.4 local variables
1279   
1280    REAL(r_std)                                   :: dt_days_read             !! STOMATE time step read in restart file (days)
1281    INTEGER(i_std)                                :: l,k,ji, jv, i, j, m      !! indices   
1282    REAL(r_std),PARAMETER                         :: max_dt_days = 5.         !! Maximum STOMATE time step (days)
1283    REAL(r_std)                                   :: hist_days                !! Writing frequency for history file (days)
1284    REAL(r_std),DIMENSION(0:nslm)                 :: z_soil                   !! Variable to store depth of the different soil
1285                                                                              !! layers (m)
1286    REAL(r_std),DIMENSION(kjpindex,nvm)           :: rprof                    !! Coefficient of the exponential functions that
1287                                                                              !! relates root density to soil depth (unitless)
1288    REAL(r_std),DIMENSION(kjpindex)               :: cvegtot                  !! Total "vegetation" cover (unitless)
1289    REAL(r_std),DIMENSION(kjpindex)               :: precip                   !! Total liquid and solid precipitation 
1290                                                                              !! @tex $(??mm dt_stomate^{-1})$ @endtex
1291    REAL(r_std),DIMENSION(kjpindex,nvm)           :: gpp_d                    !! Gross primary productivity per ground area
1292                                                                              !! @tex $(??gC m^{-2} dt_stomate^{-1})$ @endtex 
1293    REAL(r_std),DIMENSION(kjpindex,nvm)           :: gpp_daily_x              !! "Daily" gpp for teststomate 
1294                                                                              !! @tex $(??gC m^{-2} dt_stomate^{-1})$ @endtex
1295    REAL(r_std),DIMENSION(kjpindex,nvm)           :: resp_hetero_litter       !! Litter heterotrophic respiration per ground area
1296                                                                              !! @tex $(gC m^{-2} day^{-1})$ @endtex 
1297                                                                              !! ??Same variable is also used to
1298                                                                              !! store heterotrophic respiration per ground area
1299                                                                              !! over ::dt_sechiba??
1300    REAL(r_std),DIMENSION(kjpindex,nvm)           :: resp_hetero_soil         !! soil heterotrophic respiration 
1301                                                                              !! @tex $(gC m^{-2} day^{-1})$ @endtex
1302    REAL(r_std),DIMENSION(kjpindex,nvm)           :: veget_cov                !! Fractional coverage: actually share of the pixel
1303                                                                              !! covered by a PFT (fraction of ground area),
1304                                                                              !! taking into account LAI ??(= grid scale fpc)??
1305    REAL(r_std),DIMENSION(kjpindex,nvm)           :: veget_cov_max_new        !! New value for maximal fractional coverage (unitless)
1306    REAL(r_std),DIMENSION(kjpindex,nvm)           :: vcmax                    !! Maximum rate of carboxylation
1307                                                                              !! @tex $(\mumol m^{-2} s^{-1})$ @endtex
1308    REAL(r_std),DIMENSION(kjpindex,nlevs)         :: control_moist_inst       !! Moisture control of heterotrophic respiration
1309                                                                              !! (0-1, unitless)
1310    REAL(r_std),DIMENSION(kjpindex,nlevs)         :: control_temp_inst        !! Temperature control of heterotrophic
1311                                                                              !! respiration, above and below (0-1, unitless)
1312    REAL(r_std),DIMENSION(kjpindex,ncarb,nvm)     :: soilcarbon_input_inst    !! Quantity of carbon going into carbon pools from
1313                                                                              !! litter decomposition
1314                                                                              !! @tex $(gC m^{-2} day^{-1})$ @endtex
1315   
1316    INTEGER(i_std)                                :: ier                      !! Check errors in netcdf call (unitless)
1317    REAL(r_std)                                   :: sf_time                  !! Intermediate variable to calculate current time
1318                                                                              !! step
1319    INTEGER(i_std)                                :: max_totsize              !! Memory management - maximum memory size (Mb)
1320    INTEGER(i_std)                                :: totsize_1step            !! Memory management - memory required to store one
1321                                                                              !! time step on one processor (Mb)
1322    INTEGER(i_std)                                :: totsize_tmp              !! Memory management - memory required to store one
1323                                                                              !! time step on all processors(Mb)
1324    REAL(r_std)                                   :: xn                       !! How many times have we treated in this forcing
1325    REAL(r_std), DIMENSION(kjpindex)              :: vartmp                   !! Temporary variable
1326    INTEGER(i_std)                                :: vid                      !! Variable identifer of netCDF (unitless)
1327    INTEGER(i_std)                                :: nneigh                   !! Number of neighbouring pixels
1328    INTEGER(i_std)                                :: direct                   !! ??
1329    INTEGER(i_std),DIMENSION(ndm)                 :: d_id                     !! ??
1330    REAL(r_std)                                   :: net_nep_monthly          !! Integrated nep_monthly over all grid-cells on local domain
1331    REAL(r_std)                                   :: net_nep_monthly_sum      !! Integrated nep_monthly over all grid-cells on total domain(global)
1332    REAL(r_std),DIMENSION(nbp_glo)                :: clay_g                   !! Clay fraction of soil (0-1, unitless), parallel
1333                                                                              !! computing
1334    REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:,:)    :: soilcarbon_input_g       !! Quantity of carbon going into carbon pools from
1335                                                                              !! litter decomposition 
1336                                                                              !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex, parallel
1337                                                                              !! computing
1338    REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)      :: control_moist_g          !! Moisture control of heterotrophic respiration
1339                                                                              !! (0-1, unitless), parallel computing
1340    REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)      :: control_temp_g           !! Temperature control of heterotrophic respiration
1341                                                                              !! (0-1, unitless), parallel computing
1342    REAL(r_std),ALLOCATABLE,DIMENSION(:,:)        :: npp_equil_g              !! Equilibrium NPP written to forcesoil
1343                                                                              !! @tex $(gC m^{-2} year^{-1})$ @endtex, parallel
1344                                                                              !! computing
1345
1346    REAL(r_std)                                   :: net_cflux_prod_monthly_sum    !! AR5 output?? gC m2 month-1 (one variable for
1347                                                                                   !! reduce_sum and one for bcast??), parallel
1348                                                                                   !! computing
1349    REAL(r_std)                                   :: net_cflux_prod_monthly_tot    !! AR5 output?? gC m2 month-1 (one variable for
1350                                                                                   !! reduce_sum and one for bcast??), parallel
1351                                                                                   !! computing
1352    REAL(r_std)                                   :: net_harvest_above_monthly_sum !! AR5 output?? gC m2 month-1 (one variable for
1353                                                                                   !! reduce_sum and one for bcast??), parallel
1354                                                                                   !! computing
1355    REAL(r_std)                                   :: net_harvest_above_monthly_tot !! AR5 output?? gC m2 month-1 (one variable for
1356                                                                                   !! reduce_sum and one for bcast??), parallel
1357                                                                                   !! computing
1358    REAL(r_std)                                   :: net_biosp_prod_monthly_sum    !! AR5 output?? gC m2 month-1 (one variable for
1359                                                                                   !! reduce_sum and one for bcast??), parallel
1360                                                                                   !! computing
1361    REAL(r_std)                                   :: net_biosp_prod_monthly_tot    !! AR5 output?? gC m2 month-1 (one variable for
1362                                                                                   !! reduce_sum and one for bcast??), parallel
1363                                                                                   !! computing
1364    REAL(r_std), DIMENSION(kjpindex,nvm,nbpools)  :: carbon_stock                  !! Array containing the carbon stock for each pool
1365                                                                                   !! used by ORCHIDEE
1366
1367!_ ================================================================================================================================
1368   
1369  !! 1. Initialize variables
1370
1371    !! 1.1 Store current time step in a common variable
1372    itime = kjit
1373   
1374    !! 1.3 PFT rooting depth across pixels, humescte is pre-defined
1375    ! (constantes_veg.f90). It is defined as the coefficient of an exponential
1376    ! function relating root density to depth
1377    DO j=1,nvm
1378       rprof(:,j) = 1./humcste(j)
1379    ENDDO
1380   
1381    !! 1.4 Initialize first call
1382    ! Set growth respiration to zero
1383    resp_growth=zero
1384
1385    ! Check that initialization is done
1386    IF (l_first_stomate) CALL ipslerr_p(3,'stomate_main','Initialization not yet done.','','')
1387   
1388    IF (printlev >= 4) THEN
1389       WRITE(numout,*) 'stomate_main: date=',days_since_beg,' ymds=', year_end, month_end, day_end, sec_end, &
1390            ' itime=', itime, ' do_slow=',do_slow
1391    ENDIF
1392
1393!! 3. Special treatment for some input arrays.
1394   
1395    !! 3.1 Sum of liquid and solid precipitation
1396    precip(:) = ( precip_rain(:) + precip_snow(:) )*one_day/dt_sechiba
1397   
1398    !! 3.2 Calculate STOMATE's vegetation fractions from veget and veget_max
1399    DO j=1,nvm 
1400       WHERE ((1.-totfrac_nobio(:)) > min_sechiba)
1401          ! Pixels with vegetation
1402          veget_cov(:,j) = veget(:,j)/( 1.-totfrac_nobio(:) )
1403          veget_cov_max(:,j) = veget_max(:,j)/( 1.-totfrac_nobio(:) )
1404       ELSEWHERE
1405          ! Pixels without vegetation
1406          veget_cov(:,j) = zero
1407          veget_cov_max(:,j) = zero
1408       ENDWHERE
1409    ENDDO
1410
1411    IF ( do_now_stomate_lcchange ) THEN
1412       DO j=1,nvm
1413          WHERE ((1.-totfrac_nobio_new(:)) > min_sechiba)
1414             ! Pixels with vegetation
1415             veget_cov_max_new(:,j) = veget_max_new(:,j)/( 1.-totfrac_nobio_new(:) )
1416          ELSEWHERE
1417             ! Pixels without vegetation
1418             veget_cov_max_new(:,j) = zero
1419          ENDWHERE
1420       ENDDO
1421    ENDIF
1422
1423    !! 3.3 Adjust time step of GPP
1424    ! No GPP for bare soil
1425    gpp_d(:,1) = zero
1426    ! GPP per PFT
1427    DO j = 2,nvm   
1428       WHERE (veget_cov_max(:,j) > min_stomate)
1429          ! The PFT is available on the pixel
1430          gpp_d(:,j) =  gpp(:,j)/ veget_cov_max(:,j)* one_day/dt_sechiba 
1431       ELSEWHERE
1432          ! The PFT is absent on the pixel
1433          gpp_d(:,j) = zero
1434       ENDWHERE
1435    ENDDO
1436
1437  !! 4. Calculate variables for dt_stomate (i.e. "daily")
1438
1439    ! Note: If dt_days /= 1, then variables 'xx_daily' (eg. half-daily or bi-daily) are by definition
1440    ! not expressed on a daily basis. This is not a problem but could be
1441    ! confusing
1442
1443    !! 4.1 Accumulate instantaneous variables (do_slow=.FALSE.)
1444    ! Accumulate instantaneous variables (do_slow=.FALSE.) and eventually
1445    ! calculate daily mean value (do_slow=.TRUE.)
1446    CALL stomate_accu (do_slow, humrel,        humrel_daily)
1447    CALL stomate_accu (do_slow, litterhumdiag, litterhum_daily)
1448    CALL stomate_accu (do_slow, temp_air,      t2m_daily)
1449    CALL stomate_accu (do_slow, temp_sol,      tsurf_daily)
1450    CALL stomate_accu (do_slow, stempdiag,     tsoil_daily)
1451    CALL stomate_accu (do_slow, shumdiag,      soilhum_daily)
1452    CALL stomate_accu (do_slow, precip,        precip_daily)
1453    CALL stomate_accu (do_slow, gpp_d,         gpp_daily)
1454   
1455    !! 4.2 Daily minimum temperature
1456    t2m_min_daily(:) = MIN( temp_air(:), t2m_min_daily(:) )
1457
1458    !! 4.3 Calculate maintenance respiration
1459    ! Note: lai is passed as output argument to overcome previous problems with
1460    ! natural and agricultural vegetation types.
1461    CALL maint_respiration &
1462         & (kjpindex,lai,temp_air,t2m_longterm,stempdiag,height,veget_cov_max, &
1463         & rprof,biomass,resp_maint_part_radia)
1464   
1465    ! Aggregate maintenance respiration across the different plant parts
1466    resp_maint_radia(:,:) = zero
1467    DO j=2,nvm
1468       DO k= 1, nparts
1469          resp_maint_radia(:,j) = resp_maint_radia(:,j) &
1470               & + resp_maint_part_radia(:,j,k)
1471       ENDDO
1472    ENDDO
1473   
1474    ! Maintenance respiration separated by plant parts
1475    resp_maint_part(:,:,:) = resp_maint_part(:,:,:) &
1476         & + resp_maint_part_radia(:,:,:)
1477   
1478    !! 4.4 Litter dynamics and litter heterothropic respiration
1479    ! Including: litter update, lignin content, PFT parts, litter decay,
1480    ! litter heterotrophic respiration, dead leaf soil cover.
1481    ! Note: there is no vertical discretisation in the soil for litter decay.
1482    turnover_littercalc(:,:,:,:) = turnover_daily(:,:,:,:) * dt_sechiba/one_day
1483    bm_to_littercalc(:,:,:,:)    = bm_to_litter(:,:,:,:) * dt_sechiba/one_day       
1484    CALL littercalc (kjpindex, &
1485         turnover_littercalc, bm_to_littercalc, &
1486         veget_cov_max, temp_sol, stempdiag, shumdiag, litterhumdiag, &
1487         litterpart, litter, dead_leaves, lignin_struc, &
1488         deadleaf_cover, resp_hetero_litter, &
1489         soilcarbon_input_inst, control_temp_inst, control_moist_inst, &
1490         matrixA, vectorB)
1491   
1492    ! Heterothropic litter respiration during time step ::dt_sechiba @tex $(gC m^{-2})$ @endtex
1493    resp_hetero_litter(:,:) = resp_hetero_litter(:,:) * dt_sechiba/one_day
1494   
1495    !! 4.5 Soil carbon dynamics and soil heterotrophic respiration
1496    ! Note: there is no vertical discretisation in the soil for litter decay.
1497    CALL soilcarbon (kjpindex, clay, &
1498         soilcarbon_input_inst, control_temp_inst, control_moist_inst, veget_cov_max, &
1499         carbon, resp_hetero_soil, matrixA)
1500   
1501    ! Heterothropic soil respiration during time step ::dt_sechiba @tex $(gC m^{-2})$ @endtex
1502    resp_hetero_soil(:,:) = resp_hetero_soil(:,:) * dt_sechiba/one_day
1503
1504    ! Total heterothrophic respiration during time step ::dt_sechiba @tex $(gC m^{-2})$ @endtex
1505    resp_hetero_radia(:,:) = resp_hetero_litter(:,:) + resp_hetero_soil(:,:)
1506    resp_hetero_d(:,:) = resp_hetero_d(:,:) + resp_hetero_radia(:,:)
1507   
1508    !! 4.6 Accumulate instantaneous variables (do_slow=.FALSE.)
1509    ! Accumulate instantaneous variables (do_slow=.FALSE.) and eventually
1510    ! calculate daily mean value (do_slow=.TRUE.)
1511    CALL stomate_accu (do_slow, control_moist_inst, control_moist_daily)
1512    CALL stomate_accu (do_slow, control_temp_inst,  control_temp_daily)
1513    CALL stomate_accu (do_slow, soilcarbon_input_inst, soilcarbon_input_daily)
1514   
1515!! 5. Daily processes - performed at the end of the day
1516   
1517    IF (do_slow) THEN
1518
1519       !! 5.1 Update lai
1520       ! Use lai from stomate
1521       ! ?? check if this is the only time ok_pheno is used??
1522       ! ?? Looks like it is the only time. But this variables probably is defined
1523       ! in stomate_constants or something, in which case, it is difficult to track.
1524       IF (ok_pheno) THEN
1525          !! 5.1.1 Update LAI
1526          ! Set lai of bare soil to zero
1527          lai(:,ibare_sechiba) = zero
1528          ! lai for all PFTs
1529          DO j = 2, nvm
1530             lai(:,j) = biomass(:,j,ileaf,icarbon)*sla(j)
1531          ENDDO
1532          frac_age(:,:,:) = leaf_frac(:,:,:)
1533       ELSE 
1534          ! 5.1.2 Use a prescribed lai
1535          ! WARNING: code in setlai is identical to the lines above
1536          ! Update subroutine if LAI has to be forced
1537          CALL  setlai(kjpindex,lai) 
1538          frac_age(:,:,:) = zero
1539       ENDIF
1540
1541       !! 5.2 Calculate long-term "meteorological" and biological parameters
1542       ! mainly in support of calculating phenology. If LastTsYear=.TRUE.
1543       ! annual values are update (i.e. xx_lastyear).
1544       CALL season &
1545            &          (kjpindex, dt_days, &
1546            &           veget_cov, veget_cov_max, &
1547            &           humrel_daily, t2m_daily, tsoil_daily, soilhum_daily, lalo, &
1548            &           precip_daily, npp_daily, biomass, &
1549            &           turnover_daily, gpp_daily, when_growthinit, &
1550            &           maxhumrel_lastyear, maxhumrel_thisyear, &
1551            &           minhumrel_lastyear, minhumrel_thisyear, &
1552            &           maxgppweek_lastyear, maxgppweek_thisyear, &
1553            &           gdd0_lastyear, gdd0_thisyear, &
1554            &           precip_lastyear, precip_thisyear, &
1555            &           lm_lastyearmax, lm_thisyearmax, &
1556            &           maxfpc_lastyear, maxfpc_thisyear, &
1557            &           humrel_month, humrel_week, t2m_longterm, tau_longterm, &
1558            &           t2m_month, t2m_week, tsoil_month, soilhum_month, &
1559            &           npp_longterm, turnover_longterm, gpp_week, &
1560            &           gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
1561            &           time_hum_min, hum_min_dormance, gdd_init_date, &
1562            &           gdd_from_growthinit, herbivores, &
1563            &           Tseason, Tseason_length, Tseason_tmp, &
1564            &           Tmin_spring_time, t2m_min_daily, begin_leaves, onset_date)
1565       
1566       !! 5.3 Use all processes included in stomate
1567
1568       !! 5.3.1  Activate stomate processes
1569       ! Activate stomate processes (the complete list of processes depends
1570       ! on whether the DGVM is used or not). Processes include: climate constraints
1571       ! for PFTs, PFT dynamics, Phenology, Allocation, NPP (based on GPP and
1572       ! authothropic respiration), fire, mortality, vmax, assimilation temperatures,
1573       ! all turnover processes, light competition, sapling establishment, lai and
1574       ! land cover change.
1575       CALL StomateLpj &
1576            &            (kjpindex, dt_days, &
1577            &             neighbours, resolution, &
1578            &             clay, herbivores, &
1579            &             tsurf_daily, tsoil_daily, t2m_daily, t2m_min_daily, &
1580            &             litterhum_daily, soilhum_daily, &
1581            &             maxhumrel_lastyear, minhumrel_lastyear, &
1582            &             gdd0_lastyear, precip_lastyear, &
1583            &             humrel_month, humrel_week, t2m_longterm, t2m_month, t2m_week, &
1584            &             tsoil_month, soilhum_month, &
1585            &             gdd_m5_dormance,  gdd_from_growthinit, gdd_midwinter, ncd_dormance, ngd_minus5, &
1586            &             turnover_longterm, gpp_daily, &
1587            &             time_hum_min, maxfpc_lastyear, resp_maint_part,&
1588            &             PFTpresent, age, fireindex, firelitter, &
1589            &             leaf_age, leaf_frac, biomass, ind, adapted, regenerate, &
1590            &             senescence, when_growthinit, litterpart, litter, &
1591            &             dead_leaves, carbon, lignin_struc, &
1592            &             veget_cov_max, veget_cov_max_new, woodharvest, fraclut, npp_longterm, lm_lastyearmax, &
1593            &             veget_lastlight, everywhere, need_adjacent, RIP_time, &
1594            &             lai, rprof,npp_daily, turnover_daily, turnover_time,&
1595            &             control_moist_inst, control_temp_inst, soilcarbon_input_inst, &
1596            &             co2_to_bm_dgvm, co2_fire, &
1597            &             resp_hetero_d, resp_maint_d, resp_growth_d, &
1598            &             height, deadleaf_cover, vcmax, &
1599            &             bm_to_litter,&
1600            &             prod10, prod100, flux10, flux100, &
1601            &             convflux, cflux_prod10, cflux_prod100, &
1602            &             prod10_harvest, prod100_harvest, flux10_harvest, flux100_harvest, &
1603            &             convflux_harvest, cflux_prod10_harvest, cflux_prod100_harvest, woodharvestpft, & 
1604            &             convfluxpft, fDeforestToProduct, fLulccResidue,fHarvestToProduct, &
1605            &             harvest_above, carb_mass_total, &
1606            &             fpc_max, matrixA, &
1607            &             Tseason, Tmin_spring_time, begin_leaves, onset_date)
1608       
1609
1610       !! 5.3.2 Calculate the total CO2 flux from land use change
1611
1612       ! CO2 from land-use change
1613       fco2_lu(:) = convflux(:) + cflux_prod10(:) + cflux_prod100(:) 
1614
1615       ! CO2 from wood harvest
1616       fco2_wh(:) = convflux_harvest(:) + cflux_prod10_harvest(:) + cflux_prod100_harvest(:)
1617       
1618       ! CO2 from harvest
1619       fco2_ha(:) = harvest_above(:) 
1620             
1621       !! 5.4 Calculate veget and veget_max
1622       veget_max(:,:) = zero 
1623       DO j = 1, nvm
1624          veget_max(:,j) = veget_max(:,j) + &
1625               & veget_cov_max(:,j) * ( 1.-totfrac_nobio(:) )
1626       ENDDO
1627       
1628       !! 5.5 Photosynthesis parameters
1629       assim_param(:,:,ivcmax) = zero
1630       DO j = 2,nvm
1631          assim_param(:,j,ivcmax) = vcmax(:,j)
1632       ENDDO
1633       
1634       !! 5.6 Update forcing variables for soil carbon
1635       IF (TRIM(Cforcing_name) /= 'NONE') THEN
1636          npp_tot(:) = 0
1637          DO j=2,nvm
1638             npp_tot(:) = npp_tot(:) + npp_daily(:,j)
1639          ENDDO
1640          ! ::nbyear Number of years saved for carbon spinup
1641          sf_time = MODULO(REAL(days_since_beg,r_std)-1,one_year*REAL(nbyear,r_std))
1642          iatt=FLOOR(sf_time/dt_forcesoil) + 1
1643          IF (iatt == 0) iatt = iatt_old + 1
1644          IF ((iatt<iatt_old) .and. (.not. cumul_Cforcing)) THEN
1645             nforce(:)=0
1646             soilcarbon_input(:,:,:,:) = zero
1647             control_moist(:,:,:) = zero
1648             control_temp(:,:,:) = zero
1649             npp_equil(:,:) = zero
1650          ENDIF
1651          iatt_old = iatt
1652          ! Update forcing
1653          nforce(iatt) = nforce(iatt)+1
1654          soilcarbon_input(:,:,:,iatt) = soilcarbon_input(:,:,:,iatt) + soilcarbon_input_daily(:,:,:)
1655          control_moist(:,:,iatt) = control_moist(:,:,iatt) + control_moist_daily(:,:)
1656          control_temp(:,:,iatt) = control_temp(:,:,iatt) + control_temp_daily(:,:)
1657          npp_equil(:,iatt) = npp_equil(:,iatt) + npp_tot(:)
1658       ENDIF
1659       
1660       !! 5.8 Write forcing file
1661       ! Note: if STOMATE is run in coupled mode the forcing file is written
1662       ! If run in stand-alone mode, the forcing file is read!
1663       IF ( TRIM(forcing_name) /= 'NONE' ) THEN
1664         
1665          !! 5.8.1 Convert GPP to sechiba time steps
1666          ! GPP is multiplied by coverage to obtain forcing @tex $(gC m^{-2} dt_stomate^{-1})$\f \end@tex $(m^2 m^{-2})$ @endtexonly
1667          ! @tex$ m^{-2}$ @endtex remains in the units because ::veget_cov_max is a fraction, not a
1668          ! surface area. In sechiba values are ponderated by surface and frac_no_bio.
1669          ! At the beginning of stomate, the units are converted.
1670          ! When we use forcesoil we call sechiba_main and so we need the have the same units as in sechiba.
1671          gpp_daily_x(:,:) = zero
1672          DO j = 2, nvm             
1673             gpp_daily_x(:,j) = gpp_daily_x(:,j) + &
1674              & gpp_daily(:,j) * dt_stomate / one_day * veget_cov_max(:,j)
1675          ENDDO
1676         
1677          ! Bare soil moisture availability has not been treated
1678          ! in STOMATE, update it here
1679          humrel_daily(:,ibare_sechiba) = humrel(:,ibare_sechiba)   
1680
1681          ! Update index to store the next forcing step in memory
1682          iisf = iisf+1
1683
1684          ! How many times have we treated this forcing state
1685          xn = REAL(nf_cumul(isf(iisf)),r_std)
1686         
1687          !! 5.8.2 Cumulate forcing variables
1688          ! Cumulate forcing variables (calculate average)
1689          ! Note: precipitation is multiplied by dt_stomate/one_day to be consistent with
1690          ! the units in sechiba
1691          IF (cumul_forcing) THEN
1692             clay_fm(:,iisf) = (xn*clay_fm(:,iisf)+clay(:))/(xn+1.)
1693             humrel_daily_fm(:,:,iisf) = &
1694                  & (xn*humrel_daily_fm(:,:,iisf) + humrel_daily(:,:))/(xn+1.)
1695             litterhum_daily_fm(:,iisf) = &
1696                  & (xn*litterhum_daily_fm(:,iisf)+litterhum_daily(:))/(xn+1.)
1697             t2m_daily_fm(:,iisf) = &
1698                  & (xn*t2m_daily_fm(:,iisf)+t2m_daily(:))/(xn+1.)
1699             t2m_min_daily_fm(:,iisf) = &
1700                  & (xn*t2m_min_daily_fm(:,iisf)+t2m_min_daily(:))/(xn+1.)
1701             tsurf_daily_fm(:,iisf) = &
1702                  & (xn*tsurf_daily_fm(:,iisf)+tsurf_daily(:))/(xn+1.)
1703             tsoil_daily_fm(:,:,iisf) = &
1704                  & (xn*tsoil_daily_fm(:,:,iisf)+tsoil_daily(:,:))/(xn+1.)
1705             soilhum_daily_fm(:,:,iisf) = &
1706                  & (xn*soilhum_daily_fm(:,:,iisf)+soilhum_daily(:,:))/(xn+1.)
1707             precip_fm(:,iisf) = &
1708                  & (xn*precip_fm(:,iisf)+precip_daily(:)*dt_stomate/one_day)/(xn+1.)
1709             gpp_daily_fm(:,:,iisf) = &
1710                  & (xn*gpp_daily_fm(:,:,iisf) + gpp_daily_x(:,:))/(xn+1.)
1711             veget_fm(:,:,iisf) = &
1712                  & (xn*veget_fm(:,:,iisf) + veget(:,:) )/(xn+1.)
1713             veget_max_fm(:,:,iisf) = &
1714                  & (xn*veget_max_fm(:,:,iisf) + veget_max(:,:) )/(xn+1.)
1715             lai_fm(:,:,iisf) = &
1716                  & (xn*lai_fm(:,:,iisf) + lai(:,:) )/(xn+1.)
1717          ELSE
1718             ! Here we just calculate the values
1719             clay_fm(:,iisf) = clay(:)
1720             humrel_daily_fm(:,:,iisf) = humrel_daily(:,:)
1721             litterhum_daily_fm(:,iisf) = litterhum_daily(:)
1722             t2m_daily_fm(:,iisf) = t2m_daily(:)
1723             t2m_min_daily_fm(:,iisf) =t2m_min_daily(:)
1724             tsurf_daily_fm(:,iisf) = tsurf_daily(:)
1725             tsoil_daily_fm(:,:,iisf) =tsoil_daily(:,:)
1726             soilhum_daily_fm(:,:,iisf) =soilhum_daily(:,:)
1727             precip_fm(:,iisf) = precip_daily(:)
1728             gpp_daily_fm(:,:,iisf) =gpp_daily_x(:,:)
1729             veget_fm(:,:,iisf) = veget(:,:)
1730             veget_max_fm(:,:,iisf) =veget_max(:,:)
1731             lai_fm(:,:,iisf) =lai(:,:)
1732          ENDIF
1733          nf_cumul(isf(iisf)) = nf_cumul(isf(iisf))+1
1734
1735          ! 5.8.3 Do we have to write the forcing states?
1736          IF (iisf == nsfm) THEN
1737
1738             !! 5.8.3.1 Write these forcing states
1739             CALL forcing_write(forcing_id,1,nsfm)
1740             ! determine which forcing states must be read
1741             isf(1) = isf(nsfm)+1
1742             IF ( isf(1) > nsft ) isf(1) = 1
1743             DO iisf = 2, nsfm
1744                isf(iisf) = isf(iisf-1)+1
1745                IF (isf(iisf) > nsft)  isf(iisf) = 1
1746             ENDDO
1747
1748             ! Read forcing variables - for debug use only
1749             ! CALL forcing_read(forcing_id,nsfm)
1750             iisf = 0
1751
1752          ENDIF
1753
1754       ENDIF
1755
1756
1757       !! 5.9 Compute daily CO2 flux diagnostics
1758       ! CO2 flux in @tex $gC m^{-2} s^{-1}$ @endtex (positive from atmosphere to land) is sum of:
1759       !   (1) co2 taken up by photosyntyhesis + (2) co2 taken up in the DGVM to establish saplings
1760       ! - (3) plants maintenance respiration  - (4) plants growth respiration
1761       ! - (5) heterotrophic respiration from ground
1762       ! - (6) co2 emission from fire
1763       ! co2_to_bm is not added as it is already encounted in gpp
1764       nep_daily(:,:)= gpp_daily(:,:)       &
1765                     - resp_maint_d(:,:)  - resp_growth_d(:,:)   &
1766                     - resp_hetero_d(:,:) - co2_fire(:,:) 
1767
1768       CALL xios_orchidee_send_field("nep",SUM(nep_daily*veget_cov_max,dim=2)/1e3/one_day)
1769       CALL xios_orchidee_send_field("rhSoil",SUM(resp_hetero_soil*veget_cov_max,dim=2)/1e3)
1770       CALL xios_orchidee_send_field("rhLitter",SUM(resp_hetero_litter*veget_cov_max,dim=2)/1e3)
1771
1772       ! Calculate co2_flux as (-1)*nep_daily*veget_cov_max.
1773       ! This variable will be used for the coupling to LMDZ for ESM configuration.
1774       co2_flux(:,:) = (resp_hetero_d(:,:) + resp_maint_d(:,:) + resp_growth_d(:,:) &
1775            + co2_fire(:,:) -  gpp_daily(:,:))*veget_cov_max
1776     
1777       IF ( hist_id_stom_IPCC > 0 ) THEN
1778          vartmp(:) = SUM(nep_daily*veget_cov_max,dim=2)/1e3/one_day*contfrac
1779          CALL histwrite_p (hist_id_stom_IPCC, "nep", itime, &
1780               vartmp, kjpindex, hori_index)
1781       ENDIF
1782
1783       ! Cumulate nep, harvest and land use change fluxes
1784       nep_monthly(:,:) = nep_monthly(:,:) + nep_daily(:,:)
1785       harvest_above_monthly(:) = harvest_above_monthly(:) + harvest_above(:)
1786       cflux_prod_monthly(:) = cflux_prod_monthly(:) + convflux(:) + & 
1787        & cflux_prod10(:) + cflux_prod100(:) + convflux_harvest(:) + & 
1788        & cflux_prod10_harvest(:) + cflux_prod100_harvest(:)
1789     
1790       !! 5.10 Compute monthly CO2 fluxes
1791       IF ( LastTsMonth ) THEN
1792          !! 5.10.1 Write history file for monthly fluxes
1793          CALL histwrite_p (hist_id_stomate, 'CO2FLUX', itime, &
1794               nep_monthly, kjpindex*nvm, horipft_index)
1795         
1796          ! Integrate nep_monthly over all grid-cells on local domain
1797          net_nep_monthly = zero
1798          DO ji=1,kjpindex
1799             DO j=2,nvm
1800                net_nep_monthly = net_nep_monthly + &
1801                     nep_monthly(ji,j)*resolution(ji,1)*resolution(ji,2)*contfrac(ji)*veget_cov_max(ji,j)
1802             ENDDO
1803          ENDDO
1804          ! Change unit from gC/m2 grid-cell into PgC/m2 grid-cell
1805          net_nep_monthly = net_nep_monthly*1e-15
1806
1807     
1808          !! 5.10.2 Cumulative fluxes of land use cover change, harvest and net biosphere production
1809          ! Parallel processing, gather the information from different processors. first argument is the
1810          ! local variable, the second argument is the global variable. bcast send it to all processors.
1811          net_cflux_prod_monthly_sum = &
1812              &  SUM(cflux_prod_monthly(:)*resolution(:,1)*resolution(:,2)*contfrac(:))*1e-15
1813          CALL reduce_sum(net_cflux_prod_monthly_sum,net_cflux_prod_monthly_tot)
1814          CALL bcast(net_cflux_prod_monthly_tot)
1815          net_harvest_above_monthly_sum = &
1816             &   SUM(harvest_above_monthly(:)*resolution(:,1)*resolution(:,2)*contfrac(:))*1e-15
1817          CALL reduce_sum(net_harvest_above_monthly_sum,net_harvest_above_monthly_tot)
1818          CALL bcast(net_harvest_above_monthly_tot)
1819          CALL reduce_sum(net_nep_monthly,net_nep_monthly_sum)
1820          CALL bcast(net_nep_monthly_sum)
1821          net_biosp_prod_monthly_tot =  net_cflux_prod_monthly_tot + net_harvest_above_monthly_tot - net_nep_monthly_sum
1822         
1823          WRITE(numout,9010) 'GLOBAL net_cflux_prod_monthly    (Peta gC/month)  = ',net_cflux_prod_monthly_tot
1824          WRITE(numout,9010) 'GLOBAL net_harvest_above_monthly (Peta gC/month)  = ',net_harvest_above_monthly_tot
1825          WRITE(numout,9010) 'GLOBAL net_nep_monthly           (Peta gC/month)  = ',net_nep_monthly_sum
1826          WRITE(numout,9010) 'GLOBAL net_biosp_prod_monthly    (Peta gC/month)  = ',net_biosp_prod_monthly_tot
1827
18289010  FORMAT(A52,F17.14)
1829
1830          ! Reset Monthly values
1831          nep_monthly(:,:) = zero
1832          harvest_above_monthly(:) = zero
1833          cflux_prod_monthly(:)    = zero
1834
1835       ENDIF ! Monthly processes - at the end of the month
1836       
1837       IF (spinup_analytic) THEN
1838          nbp_accu(:) = nbp_accu(:) + (SUM(nep_daily(:,:) * veget_max(:,:),dim=2) - (convflux(:) + cflux_prod10(:) + &
1839                    cflux_prod100(:)) - (convflux_harvest(:) + cflux_prod10_harvest(:) + &
1840                    cflux_prod100_harvest(:))  - harvest_above(:))/1e3 
1841       ENDIF
1842
1843       !! 5.11 Reset daily variables
1844       humrel_daily(:,:) = zero
1845       litterhum_daily(:) = zero
1846       t2m_daily(:) = zero
1847       t2m_min_daily(:) = large_value
1848       tsurf_daily(:) = zero
1849       tsoil_daily(:,:) = zero
1850       soilhum_daily(:,:) = zero
1851       precip_daily(:) = zero
1852       gpp_daily(:,:) = zero
1853       resp_maint_part(:,:,:)=zero
1854       resp_hetero_d=zero
1855       IF (printlev >= 3) THEN
1856          WRITE(numout,*) 'stomate_main: daily processes done'
1857       ENDIF
1858
1859    ENDIF  ! Daily processes - at the end of the day
1860   
1861  !! 6. Outputs from Stomate
1862
1863    !! 6.1 Respiration and fluxes
1864    resp_maint(:,:) = resp_maint_radia(:,:)*veget_cov_max(:,:)
1865    resp_maint(:,ibare_sechiba) = zero
1866    resp_growth(:,:)= resp_growth_d(:,:)*veget_cov_max(:,:)*dt_sechiba/one_day
1867    co2_to_bm_radia(:,:)=co2_to_bm_dgvm(:,:)*veget_cov_max(:,:)*dt_sechiba/one_day
1868    resp_hetero(:,:) = resp_hetero_radia(:,:)*veget_cov_max(:,:)
1869   
1870    temp_growth(:)=t2m_month(:)-tp_00 
1871
1872
1873    ! Copy module variables into local variables to allow them to be in the argument output list of the subroutine
1874    co2_flux_out(:,:)=co2_flux(:,:)
1875    fco2_lu_out(:)=fco2_lu(:)
1876    fco2_wh_out(:)=fco2_wh(:)
1877    fco2_ha_out(:)=fco2_ha(:)
1878
1879
1880  !! 7. Analytical spinup
1881
1882    IF (spinup_analytic) THEN
1883
1884       !! 7.1. Update V and U at sechiba time step
1885       DO m = 2,nvm
1886          DO j = 1,kjpindex 
1887             ! V <- A * V
1888             MatrixV(j,m,:,:) = MATMUL(matrixA(j,m,:,:),MatrixV(j,m,:,:))
1889             ! U <- A*U + B
1890             VectorU(j,m,:) = MATMUL(matrixA(j,m,:,:),VectorU(j,m,:)) + vectorB(j,m,:)
1891          ENDDO ! loop pixels
1892       ENDDO ! loop PFTS
1893
1894
1895       !! 7.2. What happened at the end of the year ?
1896       IF (LastTsYear) THEN
1897
1898          !
1899          ! 7.2.1 Increase the years counter every LastTsYear which is the last sechiba time step of each year
1900          !
1901          global_years = global_years + 1 
1902
1903
1904          !
1905          ! 7.2.3 Is global_years is a multiple of the period time ?
1906          !
1907
1908          !
1909          ! 3.2.1 When global_years is a multiple of the spinup_period, we calculate :
1910          !       1) the mean nbp flux over the period. This value is restarted
1911          !       2) we solve the matrix system by Gauss Jordan method
1912          !       3) We test if a point is at equilibrium : if yes, we mark the point (ok_equilibrium array)
1913          !       4) Then we reset the matrix
1914          !       5) We erase the carbon_stock calculated by ORCHIDEE by the one found by the method
1915          IF( MOD(global_years, spinup_period) == 0 ) THEN
1916             WRITE(numout,*) 'Spinup analytic : Calculate if system is in equlibrium. global_years=',global_years
1917             ! The number total of days during the forcing period is given by :
1918             !    spinup_period*365 (we consider only the noleap calendar)
1919             nbp_flux(:) = nbp_accu(:) / ( spinup_period * 365.)
1920             ! Reset the values
1921             nbp_accu(:) = zero
1922
1923             carbon_stock(:,ibare_sechiba,:) = zero
1924             ! Prepare the matrix for the resolution
1925             ! Add a temporary matrix W which contains I-MatrixV
1926             ! we should take the opposite of matrixV and add the identitiy : we solve (I-MatrixV)*C = VectorU
1927             MatrixW(:,:,:,:) = moins_un * MatrixV(:,:,:,:)
1928             DO jv = 1,nbpools
1929                MatrixW(:,:,jv,jv) =  MatrixW(:,:,jv,jv) + un
1930             ENDDO
1931             carbon_stock(:,:,:) = VectorU(:,:,:)
1932
1933             !
1934             !  Solve the linear system
1935             !
1936             DO m = 2,nvm
1937                DO j = 1,kjpindex
1938                   ! the solution will be stored in VectorU : so it should be restarted before
1939                   ! loop over npts and nvm, so we solved npts*(nvm-1) (7,7) linear systems
1940                   CALL gauss_jordan_method(nbpools,MatrixW(j,m,:,:),carbon_stock(j,m,:))
1941                ENDDO ! loop pixels
1942             ENDDO ! loop PFTS
1943
1944             ! Reset temporary matrixW
1945             MatrixW(:,:,:,:) = zero 
1946
1947
1948             previous_stock(:,:,:) = current_stock(:,:,:)
1949             current_stock(:,:,:) = carbon_stock(:,:,:) 
1950             ! The relative error is calculated over the passive carbon pool (sum over the pfts) over the pixel.
1951             CALL error_L1_passive(kjpindex,nvm, nbpools, current_stock, previous_stock, veget_max, &
1952                  &                eps_carbon, carbon_eq)   
1953
1954             !! ok_equilibrium is saved,
1955             WHERE( carbon_eq(:) .AND. .NOT.(ok_equilibrium(:)) )
1956                ok_equilibrium(:) = .TRUE. 
1957             ENDWHERE
1958
1959             ! Reset matrixV for the pixel to the identity matrix and vectorU to zero
1960             MatrixV(:,:,:,:) = zero
1961             VectorU(:,:,:) = zero
1962             DO jv = 1,nbpools
1963                MatrixV(:,:,jv,jv) = un
1964             END DO
1965             IF (printlev >= 2) WRITE(numout,*) 'Reset for matrixV and VectorU done'   
1966
1967             !! Write the values found in the standard outputs of ORCHIDEE
1968             litter(:,istructural,:,iabove,icarbon) = carbon_stock(:,:,istructural_above)
1969             litter(:,istructural,:,ibelow,icarbon) = carbon_stock(:,:,istructural_below)
1970             litter(:,imetabolic,:,iabove,icarbon) = carbon_stock(:,:,imetabolic_above)
1971             litter(:,imetabolic,:,ibelow,icarbon) = carbon_stock(:,:,imetabolic_below)
1972             carbon(:,iactive,:) = carbon_stock(:,:,iactive_pool)
1973             carbon(:,islow,:) = carbon_stock(:,:,islow_pool)
1974             carbon(:,ipassive,:) = carbon_stock(:,:,ipassive_pool) 
1975
1976             ! Final step, test if all points at the local domain are at equilibrium
1977             ! The simulation can be stopped when all local domains have reached the equilibrium
1978             IF (printlev >=1) THEN
1979                IF (ALL(ok_equilibrium)) THEN
1980                   WRITE(numout,*) 'Spinup analytic : Equilibrium for carbon pools is reached for current local domain'
1981                ELSE
1982                   WRITE(numout,*) 'Spinup analytic : Equilibrium for carbon pools is not yet reached for current local domain'
1983                END IF
1984             END IF
1985          ENDIF ! ( MOD(global_years,spinup_period) == 0)
1986       ENDIF ! (LastTsYear)
1987
1988    ENDIF !(spinup_analytic)
1989   
1990    IF (printlev >= 4) WRITE(numout,*) 'Leaving stomate_main'
1991
1992  END SUBROUTINE stomate_main
1993
1994!! ================================================================================================================================
1995!! SUBROUTINE   : stomate_finalize
1996!!
1997!>\BRIEF        Write variables to restart file
1998!!
1999!! DESCRIPTION  : Write variables to restart file
2000!! RECENT CHANGE(S) : None
2001!!
2002!! MAIN OUTPUT VARIABLE(S):
2003!!
2004!! REFERENCES   :
2005!!
2006!! \n
2007!_ ================================================================================================================================
2008
2009  SUBROUTINE stomate_finalize (kjit, kjpindex, index, clay, co2_to_bm_radia, assim_param) 
2010   
2011    IMPLICIT NONE
2012   
2013    !! 0. Variable and parameter declaration
2014    !! 0.1 Input variables
2015    INTEGER(i_std),INTENT(in)                       :: kjit              !! Time step number (unitless)
2016    INTEGER(i_std),INTENT(in)                       :: kjpindex          !! Domain size - terrestrial pixels only (unitless)
2017    INTEGER(i_std),DIMENSION(kjpindex),INTENT(in)   :: index             !! Indices of the terrestrial pixels only (unitless)
2018    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: clay              !! Clay fraction of soil (0-1, unitless)
2019    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: co2_to_bm_radia   !! virtual gpp flux between atmosphere and biosphere
2020    REAL(r_std),DIMENSION(kjpindex,nvm,npco2),INTENT(in) :: assim_param  !! min+max+opt temperatures (K) & vmax for photosynthesis 
2021
2022    !! 0.4 Local variables
2023    REAL(r_std)                                   :: dt_days_read             !! STOMATE time step read in restart file (days)
2024    INTEGER(i_std)                                :: l,k,ji, jv, i, j, m      !! indices   
2025    REAL(r_std),PARAMETER                         :: max_dt_days = 5.         !! Maximum STOMATE time step (days)
2026    REAL(r_std)                                   :: hist_days                !! Writing frequency for history file (days)
2027    REAL(r_std),DIMENSION(0:nslm)                 :: z_soil                   !! Variable to store depth of the different soil layers (m)
2028    REAL(r_std),DIMENSION(kjpindex)               :: cvegtot                  !! Total "vegetation" cover (unitless)
2029    REAL(r_std),DIMENSION(kjpindex)               :: precip                   !! Total liquid and solid precipitation 
2030                                                                              !! @tex $(??mm dt_stomate^{-1})$ @endtex
2031    REAL(r_std),DIMENSION(kjpindex,nvm)           :: gpp_d                    !! Gross primary productivity per ground area
2032                                                                              !! @tex $(??gC m^{-2} dt_stomate^{-1})$ @endtex 
2033    REAL(r_std),DIMENSION(kjpindex,nvm)           :: gpp_daily_x              !! "Daily" gpp for teststomate 
2034                                                                              !! @tex $(??gC m^{-2} dt_stomate^{-1})$ @endtex
2035    REAL(r_std),DIMENSION(kjpindex,nvm)           :: veget_cov                !! Fractional coverage: actually share of the pixel
2036                                                                              !! covered by a PFT (fraction of ground area),
2037                                                                              !! taking into account LAI ??(= grid scale fpc)??
2038    REAL(r_std),DIMENSION(kjpindex,nvm)           :: vcmax                    !! Maximum rate of carboxylation
2039                                                                              !! @tex $(\mumol m^{-2} s^{-1})$ @endtex
2040    REAL(r_std),DIMENSION(kjpindex,nlevs)         :: control_moist_inst       !! Moisture control of heterotrophic respiration
2041                                                                              !! (0-1, unitless)
2042    REAL(r_std),DIMENSION(kjpindex,nlevs)         :: control_temp_inst        !! Temperature control of heterotrophic
2043                                                                              !! respiration, above and below (0-1, unitless)
2044    REAL(r_std),DIMENSION(kjpindex,ncarb,nvm)     :: soilcarbon_input_inst    !! Quantity of carbon going into carbon pools from
2045                                                                              !! litter decomposition
2046                                                                              !! @tex $(gC m^{-2} day^{-1})$ @endtex
2047   
2048    INTEGER(i_std)                                :: ier                      !! Check errors in netcdf call (unitless)
2049    REAL(r_std)                                   :: sf_time                  !! Intermediate variable to calculate current time
2050                                                                              !! step
2051    INTEGER(i_std)                                :: max_totsize              !! Memory management - maximum memory size (Mb)
2052    INTEGER(i_std)                                :: totsize_1step            !! Memory management - memory required to store one
2053                                                                              !! time step on one processor (Mb)
2054    INTEGER(i_std)                                :: totsize_tmp              !! Memory management - memory required to store one
2055                                                                              !! time step on all processors(Mb)
2056    REAL(r_std)                                   :: xn                       !! How many times have we treated in this forcing
2057    REAL(r_std), DIMENSION(kjpindex)              :: vartmp                   !! Temporary variable
2058    INTEGER(i_std)                                :: vid                      !! Variable identifer of netCDF (unitless)
2059    INTEGER(i_std)                                :: nneigh                   !! Number of neighbouring pixels
2060    INTEGER(i_std)                                :: direct                   !! ??
2061    INTEGER(i_std),DIMENSION(ndm)                 :: d_id                     !! ??
2062    REAL(r_std),DIMENSION(nbp_glo)                :: clay_g                   !! Clay fraction of soil (0-1, unitless), parallel
2063                                                                              !! computing
2064    REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:,:)    :: soilcarbon_input_g       !! Quantity of carbon going into carbon pools from
2065                                                                              !! litter decomposition 
2066                                                                              !! @tex $(gC m^{-2} dt_sechiba^{-1})$ @endtex, parallel
2067                                                                              !! computing
2068    REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)      :: control_moist_g          !! Moisture control of heterotrophic respiration
2069                                                                              !! (0-1, unitless), parallel computing
2070    REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:)      :: control_temp_g           !! Temperature control of heterotrophic respiration
2071                                                                              !! (0-1, unitless), parallel computing
2072    REAL(r_std),ALLOCATABLE,DIMENSION(:,:)        :: npp_equil_g              !! Equilibrium NPP written to forcesoil
2073                                                                              !! @tex $(gC m^{-2} year^{-1})$ @endtex, parallel
2074                                                                              !! computing
2075
2076    REAL(r_std)                                   :: net_cflux_prod_monthly_sum    !! AR5 output?? gC m2 month-1 (one variable for
2077                                                                                   !! reduce_sum and one for bcast??), parallel
2078                                                                                   !! computing
2079    REAL(r_std)                                   :: net_cflux_prod_monthly_tot    !! AR5 output?? gC m2 month-1 (one variable for
2080                                                                                   !! reduce_sum and one for bcast??), parallel
2081                                                                                   !! computing
2082    REAL(r_std)                                   :: net_harvest_above_monthly_sum !! AR5 output?? gC m2 month-1 (one variable for
2083                                                                                   !! reduce_sum and one for bcast??), parallel
2084                                                                                   !! computing
2085    REAL(r_std)                                   :: net_harvest_above_monthly_tot !! AR5 output?? gC m2 month-1 (one variable for
2086                                                                                   !! reduce_sum and one for bcast??), parallel
2087                                                                                   !! computing
2088    REAL(r_std)                                   :: net_biosp_prod_monthly_sum    !! AR5 output?? gC m2 month-1 (one variable for
2089                                                                                   !! reduce_sum and one for bcast??), parallel
2090                                                                                   !! computing
2091    REAL(r_std)                                   :: net_biosp_prod_monthly_tot    !! AR5 output?? gC m2 month-1 (one variable for
2092                                                                                   !! reduce_sum and one for bcast??), parallel
2093                                                                                   !! computing
2094    REAL(r_std), DIMENSION(kjpindex,nvm,nbpools)  :: carbon_stock                  !! Array containing the carbon stock for each pool
2095                                                                                   !! used by ORCHIDEE
2096
2097!_ ================================================================================================================================
2098   
2099    !! 1. Write restart file for stomate
2100    IF (printlev>=3) WRITE (numout,*) 'Write restart file for STOMATE'
2101       
2102    CALL writerestart &
2103         (kjpindex, index, &
2104         dt_days, days_since_beg, &
2105         ind, adapted, regenerate, &
2106         humrel_daily, gdd_init_date, litterhum_daily, &
2107         t2m_daily, t2m_min_daily, tsurf_daily, tsoil_daily, &
2108         soilhum_daily, precip_daily, &
2109         gpp_daily, npp_daily, turnover_daily, &
2110         humrel_month, humrel_week, &
2111         t2m_longterm, tau_longterm, t2m_month, t2m_week, &
2112         tsoil_month, soilhum_month, fireindex, firelitter, &
2113         maxhumrel_lastyear, maxhumrel_thisyear, &
2114         minhumrel_lastyear, minhumrel_thisyear, &
2115         maxgppweek_lastyear, maxgppweek_thisyear, &
2116         gdd0_lastyear, gdd0_thisyear, &
2117         precip_lastyear, precip_thisyear, &
2118         gdd_m5_dormance, gdd_from_growthinit, gdd_midwinter, ncd_dormance, ngd_minus5, &
2119         PFTpresent, npp_longterm, lm_lastyearmax, lm_thisyearmax, &
2120         maxfpc_lastyear, maxfpc_thisyear, &
2121         turnover_longterm, gpp_week, biomass, resp_maint_part, &
2122         leaf_age, leaf_frac, &
2123         senescence, when_growthinit, age, &
2124         resp_hetero_d, resp_maint_d, resp_growth_d, co2_fire, co2_to_bm_dgvm, co2_to_bm_radia, &
2125         veget_lastlight, everywhere, need_adjacent, &
2126         RIP_time, &
2127         time_hum_min, hum_min_dormance, &
2128         litterpart, litter, dead_leaves, &
2129         carbon, lignin_struc,turnover_time,&
2130         co2_flux, fco2_lu, fco2_wh, fco2_ha, &
2131         prod10,prod100,flux10, flux100, &
2132         convflux, cflux_prod10, cflux_prod100, &
2133         prod10_harvest,prod100_harvest,flux10_harvest, flux100_harvest, &
2134         convflux_harvest, cflux_prod10_harvest, cflux_prod100_harvest, &
2135         convfluxpft, fDeforestToProduct, fLulccResidue,fHarvestToProduct, &
2136         woodharvestpft, bm_to_litter, carb_mass_total, &
2137         Tseason, Tseason_length, Tseason_tmp, &
2138         Tmin_spring_time, begin_leaves, onset_date, &
2139         global_years, ok_equilibrium, nbp_accu, nbp_flux, &
2140         MatrixV, VectorU, previous_stock, current_stock, assim_param)
2141   
2142    !! 2. Write file with variables that force general processes in stomate
2143    IF ( allow_forcing_write ) THEN
2144       IF ( TRIM(forcing_name) /= 'NONE' ) THEN 
2145          CALL forcing_write(forcing_id,1,iisf)
2146          ! Close forcing file
2147          IF (is_root_prc) ier = NF90_CLOSE (forcing_id)
2148          forcing_id=-1
2149       END IF
2150    END IF
2151   
2152    !! 3. Collect variables that force the soil processes in stomate
2153    IF (TRIM(Cforcing_name) /= 'NONE' ) THEN 
2154       
2155       !! Collet variables
2156       IF (printlev >= 1) WRITE(numout,*) 'stomate: writing the forcing file for carbon spinup'
2157       DO iatt = 1, nparan*nbyear
2158          IF ( nforce(iatt) > 0 ) THEN
2159             soilcarbon_input(:,:,:,iatt) = &
2160                  & soilcarbon_input(:,:,:,iatt)/REAL(nforce(iatt),r_std)
2161             control_moist(:,:,iatt) = &
2162                  & control_moist(:,:,iatt)/REAL(nforce(iatt),r_std)
2163             control_temp(:,:,iatt) = &
2164                  & control_temp(:,:,iatt)/REAL(nforce(iatt),r_std)
2165             npp_equil(:,iatt) = &
2166                  & npp_equil(:,iatt)/REAL(nforce(iatt),r_std)
2167          ELSE
2168             IF (printlev >= 1) THEN
2169                WRITE(numout,*) 'We have no soil carbon forcing data for this time step:', iatt
2170                WRITE(numout,*) ' -> we set them to zero'
2171             END IF
2172             soilcarbon_input(:,:,:,iatt) = zero
2173             control_moist(:,:,iatt) = zero
2174             control_temp(:,:,iatt) = zero
2175             npp_equil(:,iatt) = zero
2176          ENDIF
2177       ENDDO
2178       
2179       ! Allocate memory for parallel computing
2180       IF (is_root_prc) THEN
2181          ALLOCATE(soilcarbon_input_g(nbp_glo,ncarb,nvm,nparan*nbyear))
2182          ALLOCATE(control_moist_g(nbp_glo,nlevs,nparan*nbyear))
2183          ALLOCATE(control_temp_g(nbp_glo,nlevs,nparan*nbyear))
2184          ALLOCATE(npp_equil_g(nbp_glo,nparan*nbyear))
2185       ENDIF
2186       
2187       ! Gather distributed variables
2188       CALL gather(clay,clay_g)
2189       CALL gather(soilcarbon_input,soilcarbon_input_g)
2190       CALL gather(control_moist,control_moist_g)
2191       CALL gather(control_temp,control_temp_g)
2192       CALL gather(npp_equil,npp_equil_g)
2193       
2194       !! Create netcdf
2195       ! Create, define and populate a netcdf file containing the forcing data.
2196       ! For the root processor only (parallel computing). NF90_ are functions
2197       ! from and external library. 
2198       IF (is_root_prc) THEN
2199          IF (printlev>=2) WRITE (numout,*) 'Create Cforcing file : ',TRIM(Cforcing_name)
2200          ! Create new netCDF dataset
2201          ier = NF90_CREATE (TRIM(Cforcing_name),NF90_64BIT_OFFSET ,Cforcing_id)
2202          IF (ier /= NF90_NOERR) THEN
2203             WRITE (numout,*) 'Error in creating Cforcing file : ',TRIM(Cforcing_name)
2204             CALL ipslerr_p (3,'stomate_finalize', &
2205                  &        'PROBLEM creating Cforcing file', &
2206                  &        NF90_STRERROR(ier),'')
2207          END IF
2208         
2209          ! Add variable attribute
2210          ! Note ::nbp_glo is the number of global continental points
2211          ier = NF90_PUT_ATT (Cforcing_id,NF90_GLOBAL, &
2212               &                        'kjpindex',REAL(nbp_glo,r_std))
2213          ier = NF90_PUT_ATT (Cforcing_id,NF90_GLOBAL, &
2214               &                        'nparan',REAL(nparan,r_std))
2215          ier = NF90_PUT_ATT (Cforcing_id,NF90_GLOBAL, &
2216               &                        'nbyear',REAL(nbyear,r_std))
2217         
2218          ! Add new dimension
2219          ier = NF90_DEF_DIM (Cforcing_id,'points',nbp_glo,d_id(1))
2220          ier = NF90_DEF_DIM (Cforcing_id,'carbtype',ncarb,d_id(2))
2221          ier = NF90_DEF_DIM (Cforcing_id,'vegtype',nvm,d_id(3))
2222          ier = NF90_DEF_DIM (Cforcing_id,'level',nlevs,d_id(4))
2223          ier = NF90_DEF_DIM (Cforcing_id,'time_step',NF90_UNLIMITED,d_id(5))
2224         
2225          ! Add new variable
2226          ier = NF90_DEF_VAR (Cforcing_id,'points',    r_typ,d_id(1),vid)
2227          ier = NF90_DEF_VAR (Cforcing_id,'carbtype',  r_typ,d_id(2),vid)
2228          ier = NF90_DEF_VAR (Cforcing_id,'vegtype',   r_typ,d_id(3),vid)
2229          ier = NF90_DEF_VAR (Cforcing_id,'level',     r_typ,d_id(4),vid)
2230          ier = NF90_DEF_VAR (Cforcing_id,'time_step', r_typ,d_id(5),vid)
2231          ier = NF90_DEF_VAR (Cforcing_id,'index',     r_typ,d_id(1),vid)
2232          ier = NF90_DEF_VAR (Cforcing_id,'clay',      r_typ,d_id(1),vid)
2233          ier = NF90_DEF_VAR (Cforcing_id,'soilcarbon_input',r_typ, &
2234               &                        (/ d_id(1),d_id(2),d_id(3),d_id(5) /),vid)
2235          ier = NF90_DEF_VAR (Cforcing_id,'control_moist',r_typ, &
2236               &                        (/ d_id(1),d_id(4),d_id(5) /),vid)
2237          ier = NF90_DEF_VAR (Cforcing_id,'control_temp',r_typ, &
2238               &                        (/ d_id(1),d_id(4),d_id(5) /),vid)
2239          ier = NF90_DEF_VAR (Cforcing_id,'npp_equil',r_typ, &
2240               &                        (/ d_id(1),d_id(5) /),vid)
2241          ier = NF90_ENDDEF (Cforcing_id)
2242         
2243          ! Given the name of a varaible, nf90_inq_varid finds the variable
2244          ! ID (::vid). Put data value(s) into variable ::vid
2245          ier = NF90_INQ_VARID (Cforcing_id,'points',vid)
2246          ier = NF90_PUT_VAR (Cforcing_id,vid, &
2247               &                          (/(REAL(i,r_std),i=1,nbp_glo)/))
2248          ier = NF90_INQ_VARID (Cforcing_id,'carbtype',vid)
2249          ier = NF90_PUT_VAR (Cforcing_id,vid, &
2250               &                        (/(REAL(i,r_std),i=1,ncarb)/))
2251          ier = NF90_INQ_VARID (Cforcing_id,'vegtype',vid)
2252          ier = NF90_PUT_VAR (Cforcing_id,vid, &
2253               &                            (/(REAL(i,r_std),i=1,nvm)/))
2254          ier = NF90_INQ_VARID (Cforcing_id,'level',vid)
2255          ier = NF90_PUT_VAR (Cforcing_id,vid, &
2256               &                          (/(REAL(i,r_std),i=1,nlevs)/))
2257          ier = NF90_INQ_VARID (Cforcing_id,'time_step',vid)
2258          ier = NF90_PUT_VAR (Cforcing_id,vid, &
2259               &                          (/(REAL(i,r_std),i=1,nparan*nbyear)/))
2260          ier = NF90_INQ_VARID (Cforcing_id,'index',vid)
2261          ier = NF90_PUT_VAR (Cforcing_id,vid, REAL(index_g,r_std) )
2262          ier = NF90_INQ_VARID (Cforcing_id,'clay',vid)
2263          ier = NF90_PUT_VAR (Cforcing_id,vid, clay_g )
2264          ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',vid)
2265          ier = NF90_PUT_VAR (Cforcing_id,vid, soilcarbon_input_g )
2266          ier = NF90_INQ_VARID (Cforcing_id,'control_moist',vid)
2267          ier = NF90_PUT_VAR (Cforcing_id,vid, control_moist_g )
2268          ier = NF90_INQ_VARID (Cforcing_id,'control_temp',vid)
2269          ier = NF90_PUT_VAR (Cforcing_id,vid, control_temp_g )
2270          ier = NF90_INQ_VARID (Cforcing_id,'npp_equil',vid)
2271          ier = NF90_PUT_VAR (Cforcing_id,vid, npp_equil_g )
2272         
2273          ! Close netCDF
2274          ier = NF90_CLOSE (Cforcing_id)
2275          IF (ier /= NF90_NOERR) THEN
2276             CALL ipslerr_p (3,'stomate_finalize', &
2277                  &        'PROBLEM in closing Cforcing file', &
2278                  &        NF90_STRERROR(ier),'')
2279          END IF
2280         
2281          Cforcing_id = -1
2282       ENDIF
2283
2284       ! Clear memory
2285       IF (is_root_prc) THEN
2286          DEALLOCATE(soilcarbon_input_g)
2287          DEALLOCATE(control_moist_g)
2288          DEALLOCATE(control_temp_g)
2289          DEALLOCATE(npp_equil_g)
2290       ENDIF
2291       
2292    ENDIF
2293 
2294  END SUBROUTINE stomate_finalize
2295
2296
2297!! ================================================================================================================================
2298!! SUBROUTINE   : stomate_init
2299!!
2300!>\BRIEF        The routine is called only at the first simulation. At that
2301!! time settings and flags are read and checked for internal consistency and
2302!! memory is allocated for the variables in stomate.
2303!!
2304!! DESCRIPTION  : The routine reads the
2305!! following flags from the run definition file:
2306!! -ipd (index of grid point for online diagnostics)\n
2307!! -ok_herbivores (flag to activate herbivores)\n
2308!! -treat_expansion (flag to activate PFT expansion across a pixel\n
2309!! -harvest_agri (flag to harvest aboveground biomass from agricultural PFTs)\n
2310!! \n
2311!! Check for inconsistent setting between the following flags:
2312!! -ok_stomate\n
2313!! -ok_dgvm\n
2314!! \n
2315!! Memory is allocated for all the variables of stomate and new indexing tables
2316!! are build. New indexing tables are needed because a single pixel can conatin
2317!! several PFTs. The new indexing tables have separate indices for the different
2318!! PFTs. Similar index tables are build for land use cover change.\n
2319!! \n
2320!! Several global variables and land cover change variables are initialized to
2321!! zero.\n
2322!!
2323!! RECENT CHANGE(S) : None
2324!!
2325!! MAIN OUTPUT VARIABLE(S): Strictly speaking the subroutine has no output
2326!! variables. However, the routine allocates memory and builds new indexing
2327!! variables for later use.\n
2328!!
2329!! REFERENCE(S) : None
2330!!
2331!! FLOWCHART    : None
2332!! \n
2333!_ ================================================================================================================================
2334
2335  SUBROUTINE stomate_init &
2336       &  (kjpij, kjpindex, index, lalo, &
2337       &   rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
2338
2339  !! 0. Variable and parameter declaration
2340
2341    !! 0.1 Input variables
2342
2343    INTEGER(i_std),INTENT(in)                    :: kjpij             !! Total size of the un-compressed grid, including
2344                                                                      !! oceans (unitless)
2345    INTEGER(i_std),INTENT(in)                    :: kjpindex          !! Domain size - number of terrestrial pixels
2346                                                                      !! (unitless)
2347    INTEGER(i_std),INTENT(in)                    :: rest_id_stom      !! STOMATE's _Restart_ file identifier
2348    INTEGER(i_std),INTENT(in)                    :: hist_id_stom      !! STOMATE's _history_ file identifier
2349    INTEGER(i_std),INTENT(in)                    :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file identifier
2350    INTEGER(i_std),DIMENSION(kjpindex),INTENT(in):: index             !! Indices of the terrestrial pixels on the global
2351                                                                      !! map
2352    REAL(r_std),DIMENSION(kjpindex,2),INTENT(in) :: lalo              !! Geogr. coordinates (latitude,longitude) (degrees)
2353   
2354    !! 0.2 Output variables
2355
2356    !! 0.3 Modified variables
2357
2358    !! 0.4 Local variables
2359
2360    LOGICAL                                      :: l_error           !! Check errors in netcdf call
2361    INTEGER(i_std)                               :: ier               !! Check errors in netcdf call
2362    INTEGER(i_std)                               :: ji,j,ipd,l        !! Indices
2363!_ ================================================================================================================================
2364   
2365  !! 1. Online diagnostics
2366
2367    IF ( kjpindex > 0 ) THEN
2368       !Config  Key  = STOMATE_DIAGPT
2369       !Config  Desc = Index of grid point for online diagnostics
2370       !Config If    = OK_STOMATE
2371       !Config  Def  = 1
2372       !Config  Help = This is the index of the grid point which
2373       !               will be used for online diagnostics.
2374       !Config Units = [-]
2375       ! By default ::ipd is set to 1
2376       ipd = 1
2377       ! Get ::ipd from run definition file
2378       CALL getin_p('STOMATE_DIAGPT',ipd)
2379       ipd = MIN( ipd, kjpindex )
2380       IF ( printlev >=3 ) THEN
2381          WRITE(numout,*) 'Stomate: '
2382          WRITE(numout,*) '  Index of grid point for online diagnostics: ',ipd
2383          WRITE(numout,*) '  Lon, lat:',lalo(ipd,2),lalo(ipd,1)
2384          WRITE(numout,*) '  Index of this point on GCM grid: ',index(ipd)
2385       END IF
2386    ENDIF
2387   
2388  !! 2. Check consistency of flags
2389
2390    IF ( ( .NOT. ok_stomate ) .AND. ok_dgvm ) THEN
2391       WRITE(numout,*) 'Cannot do dynamical vegetation without STOMATE.'
2392       WRITE(numout,*) 'Inconsistency between ::ok_stomate and ::ok_dgvm'
2393       WRITE(numout,*) 'Stop: fatal error'
2394       STOP
2395    ENDIF
2396
2397  !! 3. Communicate settings
2398   
2399    IF (printlev >=2) THEN
2400       WRITE(numout,*) 'stomate first call - overview of the activated flags:'
2401       WRITE(numout,*) '  STOMATE: ', ok_stomate
2402       WRITE(numout,*) '  LPJ: ', ok_dgvm
2403    END IF
2404  !! 4. Allocate memory for STOMATE's variables
2405
2406    l_error = .FALSE.
2407
2408    ALLOCATE(veget_cov_max(kjpindex,nvm),stat=ier)
2409    l_error = l_error .OR. (ier /= 0)
2410    IF (l_error) THEN
2411       WRITE(numout,*) 'Memory allocation error for veget_cov_max. We stop. We need kjpindex*nvm words',kjpindex,nvm
2412       STOP 'stomate_init'
2413    ENDIF
2414
2415    ALLOCATE(ind(kjpindex,nvm),stat=ier)
2416    l_error = l_error .OR. (ier /= 0)
2417    IF (l_error) THEN
2418       WRITE(numout,*) 'Memory allocation error for ind. We stop. We need kjpindex*nvm words',kjpindex,nvm
2419       STOP 'stomate_init'
2420    ENDIF
2421
2422    ALLOCATE(adapted(kjpindex,nvm),stat=ier)
2423    l_error = l_error .OR. (ier /= 0)
2424    IF (l_error) THEN
2425       WRITE(numout,*) 'Memory allocation error for adapted. We stop. We need kjpindex*nvm words',kjpindex,nvm
2426       STOP 'stomate_init'
2427    ENDIF
2428
2429    ALLOCATE(regenerate(kjpindex,nvm),stat=ier)
2430    l_error = l_error .OR. (ier /= 0)
2431    IF (l_error) THEN
2432       WRITE(numout,*) 'Memory allocation error for regenerate. We stop. We need kjpindex*nvm words',kjpindex,nvm
2433       STOP 'stomate_init'
2434    ENDIF
2435
2436    ALLOCATE(humrel_daily(kjpindex,nvm),stat=ier)
2437    l_error = l_error .OR. (ier /= 0)
2438    IF (l_error) THEN
2439       WRITE(numout,*) 'Memory allocation error for humrel_daily. We stop. We need kjpindex*nvm words',kjpindex,nvm
2440       STOP 'stomate_init'
2441    ENDIF
2442
2443    ALLOCATE(litterhum_daily(kjpindex),stat=ier)
2444    l_error = l_error .OR. (ier /= 0)
2445    IF (l_error) THEN
2446       WRITE(numout,*) 'Memory allocation error for litterhum_daily. We stop. We need kjpindex words',kjpindex
2447       STOP 'stomate_init'
2448    ENDIF
2449
2450    ALLOCATE(t2m_daily(kjpindex),stat=ier)
2451    l_error = l_error .OR. (ier /= 0)
2452    IF (l_error) THEN
2453       WRITE(numout,*) 'Memory allocation error for t2m_daily. We stop. We need kjpindex words',kjpindex
2454       STOP 'stomate_init'
2455    ENDIF
2456
2457    ALLOCATE(t2m_min_daily(kjpindex),stat=ier)
2458    l_error = l_error .OR. (ier /= 0)
2459    IF (l_error) THEN
2460       WRITE(numout,*) 'Memory allocation error for t2m_min_daily. We stop. We need kjpindex words',kjpindex
2461       STOP 'stomate_init'
2462    ENDIF
2463
2464    ALLOCATE(tsurf_daily(kjpindex),stat=ier)
2465    l_error = l_error .OR. (ier /= 0)
2466    IF (l_error) THEN
2467       WRITE(numout,*) 'Memory allocation error for tsurf_daily. We stop. We need kjpindex words',kjpindex
2468       STOP 'stomate_init'
2469    ENDIF
2470
2471    ALLOCATE(tsoil_daily(kjpindex,nslm),stat=ier)
2472    l_error = l_error .OR. (ier /= 0)
2473    IF (l_error) THEN
2474       WRITE(numout,*) 'Memory allocation error for tsoil_daily. We stop. We need kjpindex*nslm words',kjpindex,nslm
2475       STOP 'stomate_init'
2476    ENDIF
2477
2478    ALLOCATE(soilhum_daily(kjpindex,nslm),stat=ier)
2479    l_error = l_error .OR. (ier /= 0)
2480    IF (l_error) THEN
2481       WRITE(numout,*) 'Memory allocation error for soilhum_daily. We stop. We need kjpindex*nslm words',kjpindex,nslm
2482       STOP 'stomate_init'
2483    ENDIF
2484
2485    ALLOCATE(precip_daily(kjpindex),stat=ier)
2486    l_error = l_error .OR. (ier /= 0)
2487    IF (l_error) THEN
2488       WRITE(numout,*) 'Memory allocation error for precip_daily. We stop. We need kjpindex words',kjpindex,nvm
2489       STOP 'stomate_init'
2490    ENDIF
2491
2492    ALLOCATE(gpp_daily(kjpindex,nvm),stat=ier)
2493    l_error = l_error .OR. (ier /= 0)
2494    IF (l_error) THEN
2495       WRITE(numout,*) 'Memory allocation error for gpp_daily. We stop. We need kjpindex*nvm words',kjpindex,nvm
2496       STOP 'stomate_init'
2497    ENDIF
2498
2499    ALLOCATE(npp_daily(kjpindex,nvm),stat=ier)
2500    l_error = l_error .OR. (ier /= 0)
2501    IF (l_error) THEN
2502       WRITE(numout,*) 'Memory allocation error for npp_daily. We stop. We need kjpindex*nvm words',kjpindex,nvm
2503       STOP 'stomate_init'
2504    ENDIF
2505
2506    ALLOCATE(turnover_daily(kjpindex,nvm,nparts,nelements),stat=ier)
2507    l_error = l_error .OR. (ier /= 0)
2508    IF (l_error) THEN
2509       WRITE(numout,*) 'Memory allocation error for turnover_daily. We stop. We need kjpindex*nvm*nparts*nelements words', &
2510       &   kjpindex,nvm,nparts,nelements
2511       STOP 'stomate_init'
2512    ENDIF
2513
2514    ALLOCATE(turnover_littercalc(kjpindex,nvm,nparts,nelements),stat=ier)
2515    l_error = l_error .OR. (ier /= 0)
2516    IF (l_error) THEN
2517       WRITE(numout,*) 'Memory allocation error for turnover_littercalc. We stop. We need kjpindex*nvm*nparts*nelements words', & 
2518        &  kjpindex,nvm,nparts,nelements
2519       STOP 'stomate_init'
2520    ENDIF
2521
2522    ALLOCATE(humrel_month(kjpindex,nvm),stat=ier)
2523    l_error = l_error .OR. (ier /= 0)
2524    IF (l_error) THEN
2525       WRITE(numout,*) 'Memory allocation error for humrel_month. We stop. We need kjpindex*nvm words',kjpindex,nvm
2526       STOP 'stomate_init'
2527    ENDIF
2528
2529    ALLOCATE(humrel_week(kjpindex,nvm),stat=ier)
2530    l_error = l_error .OR. (ier /= 0)
2531    IF (l_error) THEN
2532       WRITE(numout,*) 'Memory allocation error for humrel_week. We stop. We need kjpindex*nvm words',kjpindex,nvm
2533       STOP 'stomate_init'
2534    ENDIF
2535
2536    ALLOCATE(t2m_longterm(kjpindex),stat=ier)
2537    l_error = l_error .OR. (ier /= 0)
2538    IF (l_error) THEN
2539       WRITE(numout,*) 'Memory allocation error for t2m_longterm. We stop. We need kjpindex*nvm words',kjpindex,nvm
2540       STOP 'stomate_init'
2541    ENDIF
2542
2543    ALLOCATE(t2m_month(kjpindex),stat=ier)
2544    l_error = l_error .OR. (ier /= 0)
2545    IF (l_error) THEN
2546       WRITE(numout,*) 'Memory allocation error for t2m_month. We stop. We need kjpindex words',kjpindex
2547       STOP 'stomate_init'
2548    ENDIF
2549
2550    ALLOCATE(Tseason(kjpindex),stat=ier)
2551    l_error = l_error .OR. (ier /= 0)
2552    IF (l_error) THEN
2553       WRITE(numout,*) 'Memory allocation error for Tseason. We stop. We need kjpindex words',kjpindex
2554       STOP 'stomate_init'
2555    ENDIF
2556
2557    ALLOCATE(Tseason_length(kjpindex),stat=ier)
2558    l_error = l_error .OR. (ier /= 0)
2559    IF (l_error) THEN
2560       WRITE(numout,*) 'Memory allocation error for Tseason_length. We stop. We need kjpindex words',kjpindex
2561       STOP 'stomate_init'
2562    ENDIF
2563
2564    ALLOCATE(Tseason_tmp(kjpindex),stat=ier)
2565    l_error = l_error .OR. (ier /= 0)
2566    IF (l_error) THEN
2567       WRITE(numout,*) 'Memory allocation error for Tseason_tmp. We stop. We need kjpindex words',kjpindex
2568       STOP 'stomate_init'
2569    ENDIF
2570
2571    ALLOCATE(Tmin_spring_time(kjpindex,nvm),stat=ier)
2572    l_error = l_error .OR. (ier /= 0)
2573    IF (l_error) THEN
2574       WRITE(numout,*) 'Memory allocation error for Tmin_spring_time. We stop. We need kjpindex*nvm words',kjpindex,nvm
2575       STOP 'stomate_init'
2576    ENDIF
2577
2578    ALLOCATE(onset_date(kjpindex,nvm),stat=ier)
2579    l_error = l_error .OR. (ier /= 0)
2580    IF (l_error) THEN
2581       WRITE(numout,*) 'Memory allocation error for onset_date. We stop. We need kjpindex*nvm*nparts words',kjpindex,nvm,2
2582       STOP 'stomate_init'
2583    ENDIF
2584
2585    ALLOCATE(t2m_week(kjpindex),stat=ier)
2586    l_error = l_error .OR. (ier /= 0)
2587    IF (l_error) THEN
2588       WRITE(numout,*) 'Memory allocation error for t2m_week. We stop. We need kjpindex words',kjpindex
2589       STOP 'stomate_init'
2590    ENDIF
2591
2592    ALLOCATE(tsoil_month(kjpindex,nslm),stat=ier)
2593    l_error = l_error .OR. (ier /= 0)
2594    IF (l_error) THEN
2595       WRITE(numout,*) 'Memory allocation error for tsoil_month. We stop. We need kjpindex*nslm words',kjpindex,nslm
2596       STOP 'stomate_init'
2597    ENDIF
2598
2599    ALLOCATE(soilhum_month(kjpindex,nslm),stat=ier)
2600    l_error = l_error .OR. (ier /= 0)
2601    IF (l_error) THEN
2602       WRITE(numout,*) 'Memory allocation error for soilhum_month. We stop. We need kjpindex*nslm words',kjpindex,nslm
2603       STOP 'stomate_init'
2604    ENDIF
2605
2606    ALLOCATE(fireindex(kjpindex,nvm),stat=ier) 
2607    l_error = l_error .OR. (ier /= 0)
2608    IF (l_error) THEN
2609       WRITE(numout,*) 'Memory allocation error for fireindex. We stop. We need kjpindex*nvm words',kjpindex,nvm
2610       STOP 'stomate_init'
2611    ENDIF
2612
2613    ALLOCATE(firelitter(kjpindex,nvm),stat=ier)
2614    l_error = l_error .OR. (ier /= 0)
2615    IF (l_error) THEN
2616       WRITE(numout,*) 'Memory allocation error for firelitter. We stop. We need kjpindex*nvm words',kjpindex,nvm
2617       STOP 'stomate_init'
2618    ENDIF
2619
2620    ALLOCATE(maxhumrel_lastyear(kjpindex,nvm),stat=ier)
2621    l_error = l_error .OR. (ier /= 0)
2622    IF (l_error) THEN
2623       WRITE(numout,*) 'Memory allocation error for maxhumrel_lastyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2624       STOP 'stomate_init'
2625    ENDIF
2626
2627    ALLOCATE(maxhumrel_thisyear(kjpindex,nvm),stat=ier)
2628    l_error = l_error .OR. (ier /= 0)
2629    IF (l_error) THEN
2630       WRITE(numout,*) 'Memory allocation error for maxhumrel_thisyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2631       STOP 'stomate_init'
2632    ENDIF
2633
2634    ALLOCATE(minhumrel_lastyear(kjpindex,nvm),stat=ier)
2635    l_error = l_error .OR. (ier /= 0)
2636    IF (l_error) THEN
2637       WRITE(numout,*) 'Memory allocation error for minhumrel_lastyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2638       STOP 'stomate_init'
2639    ENDIF
2640
2641    ALLOCATE(minhumrel_thisyear(kjpindex,nvm),stat=ier)
2642    l_error = l_error .OR. (ier /= 0)
2643    IF (l_error) THEN
2644       WRITE(numout,*) 'Memory allocation error for minhumrel_thisyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2645       STOP 'stomate_init'
2646    ENDIF
2647
2648    ALLOCATE(maxgppweek_lastyear(kjpindex,nvm),stat=ier)
2649    l_error = l_error .OR. (ier /= 0)
2650    IF (l_error) THEN
2651       WRITE(numout,*) 'Memory allocation error for maxgppweek_lastyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2652       STOP 'stomate_init'
2653    ENDIF
2654
2655    ALLOCATE(maxgppweek_thisyear(kjpindex,nvm),stat=ier)
2656    l_error = l_error .OR. (ier /= 0)
2657    IF (l_error) THEN
2658       WRITE(numout,*) 'Memory allocation error for maxgppweek_thisyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2659       STOP 'stomate_init'
2660    ENDIF
2661
2662    ALLOCATE(gdd0_lastyear(kjpindex),stat=ier)
2663    l_error = l_error .OR. (ier /= 0)
2664    IF (l_error) THEN
2665       WRITE(numout,*) 'Memory allocation error for gdd0_lastyear. We stop. We need kjpindex words',kjpindex
2666       STOP 'stomate_init'
2667    ENDIF
2668
2669    ALLOCATE(gdd0_thisyear(kjpindex),stat=ier)
2670    l_error = l_error .OR. (ier /= 0)
2671    IF (l_error) THEN
2672       WRITE(numout,*) 'Memory allocation error for gdd0_thisyear. We stop. We need kjpindex words',kjpindex
2673       STOP 'stomate_init'
2674    ENDIF
2675
2676    ALLOCATE(gdd_init_date(kjpindex,2),stat=ier)
2677    l_error = l_error .OR. (ier /= 0)
2678    IF (l_error) THEN
2679       WRITE(numout,*) 'Memory allocation error for gdd_init_date. We stop. We need kjpindex*2 words',kjpindex,2
2680       STOP 'stomate_init'
2681    ENDIF
2682
2683    ALLOCATE(gdd_from_growthinit(kjpindex,nvm),stat=ier)
2684    l_error = l_error .OR. (ier /= 0)
2685    IF (l_error) THEN
2686       WRITE(numout,*) 'Memory allocation error for gdd_from_growthinit. We stop. We need kjpindex*nvm words',kjpindex,nvm
2687       STOP 'stomate_init'
2688    ENDIF
2689
2690    ALLOCATE(precip_lastyear(kjpindex),stat=ier)
2691    l_error = l_error .OR. (ier /= 0)
2692    IF (l_error) THEN
2693       WRITE(numout,*) 'Memory allocation error for precip_lastyear. We stop. We need kjpindex*nvm words',kjpindex
2694       STOP 'stomate_init'
2695    ENDIF
2696
2697    ALLOCATE(precip_thisyear(kjpindex),stat=ier)
2698    l_error = l_error .OR. (ier /= 0)
2699    IF (l_error) THEN
2700       WRITE(numout,*) 'Memory allocation error for precip_thisyear. We stop. We need kjpindex words',kjpindex
2701       STOP 'stomate_init'
2702    ENDIF
2703
2704    ALLOCATE(gdd_m5_dormance(kjpindex,nvm),stat=ier)
2705    l_error = l_error .OR. (ier /= 0)
2706    IF (l_error) THEN
2707       WRITE(numout,*) 'Memory allocation error for gdd_m5_dormance. We stop. We need kjpindex*nvm words',kjpindex,nvm
2708       STOP 'stomate_init'
2709    ENDIF
2710
2711    ALLOCATE(gdd_midwinter(kjpindex,nvm),stat=ier)
2712    l_error = l_error .OR. (ier /= 0)
2713    IF (l_error) THEN
2714       WRITE(numout,*) 'Memory allocation error for gdd_midwinter. We stop. We need kjpindex*nvm words',kjpindex,nvm
2715       STOP 'stomate_init'
2716    ENDIF
2717
2718    ALLOCATE(ncd_dormance(kjpindex,nvm),stat=ier)
2719    l_error = l_error .OR. (ier /= 0)
2720    IF (l_error) THEN
2721       WRITE(numout,*) 'Memory allocation error for ncd_dormance. We stop. We need kjpindex*nvm words',kjpindex,nvm
2722       STOP 'stomate_init'
2723    ENDIF
2724
2725    ALLOCATE(ngd_minus5(kjpindex,nvm),stat=ier)
2726    l_error = l_error .OR. (ier /= 0)
2727    IF (l_error) THEN
2728       WRITE(numout,*) 'Memory allocation error for ngd_minus5. We stop. We need kjpindex*nvm words',kjpindex,nvm
2729       STOP 'stomate_init'
2730    ENDIF
2731
2732    ALLOCATE(PFTpresent(kjpindex,nvm),stat=ier)
2733    l_error = l_error .OR. (ier /= 0)
2734    IF (l_error) THEN
2735       WRITE(numout,*) 'Memory allocation error for PFTpresent. We stop. We need kjpindex*nvm words',kjpindex,nvm
2736       STOP 'stomate_init'
2737    ENDIF
2738
2739    ALLOCATE(npp_longterm(kjpindex,nvm),stat=ier)
2740    l_error = l_error .OR. (ier /= 0)
2741    IF (l_error) THEN
2742       WRITE(numout,*) 'Memory allocation error for npp_longterm. We stop. We need kjpindex*nvm words',kjpindex,nvm
2743       STOP 'stomate_init'
2744    ENDIF
2745
2746    ALLOCATE(lm_lastyearmax(kjpindex,nvm),stat=ier)
2747    l_error = l_error .OR. (ier /= 0)
2748    IF (l_error) THEN
2749       WRITE(numout,*) 'Memory allocation error for lm_lastyearmax. We stop. We need kjpindex*nvm words',kjpindex,nvm
2750       STOP 'stomate_init'
2751    ENDIF
2752
2753    ALLOCATE(lm_thisyearmax(kjpindex,nvm),stat=ier)
2754    l_error = l_error .OR. (ier /= 0)
2755    IF (l_error) THEN
2756       WRITE(numout,*) 'Memory allocation error for lm_thisyearmax. We stop. We need kjpindex*nvm words',kjpindex,nvm
2757       STOP 'stomate_init'
2758    ENDIF
2759
2760    ALLOCATE(maxfpc_lastyear(kjpindex,nvm),stat=ier)
2761    l_error = l_error .OR. (ier /= 0)
2762    IF (l_error) THEN
2763       WRITE(numout,*) 'Memory allocation error for maxfpc_lastyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2764       STOP 'stomate_init'
2765    ENDIF
2766
2767    ALLOCATE(maxfpc_thisyear(kjpindex,nvm),stat=ier)
2768    l_error = l_error .OR. (ier /= 0)
2769    IF (l_error) THEN
2770       WRITE(numout,*) 'Memory allocation error for maxfpc_thisyear. We stop. We need kjpindex*nvm words',kjpindex,nvm
2771       STOP 'stomate_init'
2772    ENDIF
2773
2774    ALLOCATE(turnover_longterm(kjpindex,nvm,nparts,nelements),stat=ier)
2775    l_error = l_error .OR. (ier /= 0)
2776    IF (l_error) THEN
2777       WRITE(numout,*) 'Memory allocation error for turnover_longterm. We stop. We need kjpindex*nvm*nparts*nelements words', & 
2778       &    kjpindex,nvm,nparts,nelements
2779       STOP 'stomate_init'
2780    ENDIF
2781
2782    ALLOCATE(gpp_week(kjpindex,nvm),stat=ier)
2783    l_error = l_error .OR. (ier /= 0)
2784    IF (l_error) THEN
2785       WRITE(numout,*) 'Memory allocation error for gpp_week. We stop. We need kjpindex*nvm words',kjpindex,nvm
2786       STOP 'stomate_init'
2787    ENDIF
2788
2789    ALLOCATE(biomass(kjpindex,nvm,nparts,nelements),stat=ier)
2790    l_error = l_error .OR. (ier /= 0)
2791    IF (l_error) THEN
2792       WRITE(numout,*) 'Memory allocation error for biomass. We stop. We need kjpindex*nvm*nparts*nelements words', &
2793       &    kjpindex,nvm,nparts,nelements
2794       STOP 'stomate_init'
2795    ENDIF
2796
2797    ALLOCATE(senescence(kjpindex,nvm),stat=ier)
2798    l_error = l_error .OR. (ier /= 0)
2799    IF (l_error) THEN
2800       WRITE(numout,*) 'Memory allocation error for senescence. We stop. We need kjpindex*nvm words',kjpindex,nvm
2801       STOP 'stomate_init'
2802    ENDIF
2803
2804    ALLOCATE(begin_leaves(kjpindex,nvm),stat=ier)
2805    l_error = l_error .OR. (ier /= 0)
2806    IF (l_error) THEN
2807       WRITE(numout,*) 'Memory allocation error for begin_leaves. We stop. We need kjpindex*nvm words',kjpindex,nvm
2808       STOP 'stomate_init'
2809    ENDIF
2810
2811    ALLOCATE(when_growthinit(kjpindex,nvm),stat=ier)
2812    l_error = l_error .OR. (ier /= 0)
2813    IF (l_error) THEN
2814       WRITE(numout,*) 'Memory allocation error for when_growthinit. We stop. We need kjpindex*nvm words',kjpindex,nvm
2815       STOP 'stomate_init'
2816    ENDIF
2817
2818    ALLOCATE(age(kjpindex,nvm),stat=ier)
2819    l_error = l_error .OR. (ier /= 0)
2820    IF (l_error) THEN
2821       WRITE(numout,*) 'Memory allocation error for age. We stop. We need kjpindex*nvm words',kjpindex,nvm
2822       STOP 'stomate_init'
2823    ENDIF
2824
2825    ALLOCATE(resp_hetero_d(kjpindex,nvm),stat=ier)
2826    l_error = l_error .OR. (ier /= 0)
2827    IF (l_error) THEN
2828       WRITE(numout,*) 'Memory allocation error for resp_hetero_d. We stop. We need kjpindex*nvm words',kjpindex,nvm
2829       STOP 'stomate_init'
2830    ENDIF
2831
2832    ALLOCATE(resp_hetero_radia(kjpindex,nvm),stat=ier)
2833    l_error = l_error .OR. (ier /= 0)
2834    IF (l_error) THEN
2835       WRITE(numout,*) 'Memory allocation error for resp_hetero_radia. We stop. We need kjpindex*nvm words',kjpindex,nvm
2836       STOP 'stomate_init'
2837    ENDIF
2838
2839    ALLOCATE(resp_maint_d(kjpindex,nvm),stat=ier)
2840    l_error = l_error .OR. (ier /= 0)
2841    IF (l_error) THEN
2842       WRITE(numout,*) 'Memory allocation error for resp_maint_d. We stop. We need kjpindex*nvm words',kjpindex,nvm
2843       STOP 'stomate_init'
2844    ENDIF
2845
2846    ALLOCATE(resp_growth_d(kjpindex,nvm),stat=ier)
2847    l_error = l_error .OR. (ier /= 0)
2848    IF (l_error) THEN
2849       WRITE(numout,*) 'Memory allocation error for resp_growth_d. We stop. We need kjpindex*nvm words',kjpindex,nvm
2850       STOP 'stomate_init'
2851    ENDIF
2852
2853    ALLOCATE(co2_fire(kjpindex,nvm),stat=ier)
2854    l_error = l_error .OR. (ier /= 0)
2855    IF (l_error) THEN
2856       WRITE(numout,*) 'Memory allocation error for co2_fire. We stop. We need kjpindex*nvm words',kjpindex,nvm
2857       STOP 'stomate_init'
2858    ENDIF
2859
2860    ALLOCATE(co2_to_bm_dgvm(kjpindex,nvm),stat=ier)
2861    l_error = l_error .OR. (ier /= 0)
2862    IF (l_error) THEN
2863       WRITE(numout,*) 'Memory allocation error for co2_to_bm_dgvm. We stop. We need kjpindex*nvm words',kjpindex,nvm
2864       STOP 'stomate_init'
2865    ENDIF
2866
2867    ALLOCATE(veget_lastlight(kjpindex,nvm),stat=ier)
2868    l_error = l_error .OR. (ier /= 0)
2869    IF (l_error) THEN
2870       WRITE(numout,*) 'Memory allocation error for veget_lastlight. We stop. We need kjpindex*nvm words',kjpindex,nvm
2871       STOP 'stomate_init'
2872    ENDIF
2873
2874    ALLOCATE(everywhere(kjpindex,nvm),stat=ier)
2875    l_error = l_error .OR. (ier /= 0)
2876    IF (l_error) THEN
2877       WRITE(numout,*) 'Memory allocation error for everywhere. We stop. We need kjpindex*nvm words',kjpindex,nvm
2878       STOP 'stomate_init'
2879    ENDIF
2880
2881    ALLOCATE(need_adjacent(kjpindex,nvm),stat=ier)
2882    l_error = l_error .OR. (ier /= 0)
2883    IF (l_error) THEN
2884       WRITE(numout,*) 'Memory allocation error for need_adjacent. We stop. We need kjpindex*nvm words',kjpindex,nvm
2885       STOP 'stomate_init'
2886    ENDIF
2887
2888    ALLOCATE(leaf_age(kjpindex,nvm,nleafages),stat=ier)
2889    l_error = l_error .OR. (ier /= 0)
2890    IF (l_error) THEN
2891       WRITE(numout,*) 'Memory allocation error for leaf_age. We stop. We need kjpindex*nvm*nleafages words', & 
2892       &      kjpindex,nvm,nleafages
2893       STOP 'stomate_init'
2894    ENDIF
2895
2896    ALLOCATE(leaf_frac(kjpindex,nvm,nleafages),stat=ier)
2897    l_error = l_error .OR. (ier /= 0)
2898    IF (l_error) THEN
2899       WRITE(numout,*) 'Memory allocation error for leaf_frac. We stop. We need kjpindex*nvm*nleafages words', & 
2900       &      kjpindex,nvm,nleafages
2901       STOP 'stomate_init'
2902    ENDIF
2903
2904    ALLOCATE(RIP_time(kjpindex,nvm),stat=ier)
2905    l_error = l_error .OR. (ier /= 0)
2906    IF (l_error) THEN
2907       WRITE(numout,*) 'Memory allocation error for RIP_time. We stop. We need kjpindex*nvm words',kjpindex,nvm
2908       STOP 'stomate_init'
2909    ENDIF
2910
2911    ALLOCATE(time_hum_min(kjpindex,nvm),stat=ier)
2912    l_error = l_error .OR. (ier /= 0)
2913    IF (l_error) THEN
2914       WRITE(numout,*) 'Memory allocation error for time_hum_min. We stop. We need kjpindex*nvm words',kjpindex,nvm
2915       STOP 'stomate_init'
2916    ENDIF
2917
2918    ALLOCATE(hum_min_dormance(kjpindex,nvm),stat=ier)
2919    l_error = l_error .OR. (ier /= 0)
2920    IF (l_error) THEN
2921       WRITE(numout,*) 'Memory allocation error for hum_min_dormance. We stop. We need kjpindex*nvm words',kjpindex,nvm
2922       STOP 'stomate_init'
2923    ENDIF
2924
2925    ALLOCATE(litterpart(kjpindex,nvm,nlitt),stat=ier)
2926    l_error = l_error .OR. (ier /= 0)
2927    IF (l_error) THEN
2928       WRITE(numout,*) 'Memory allocation error for litterpart. We stop. We need kjpindex*nvm*nlitt words',  &
2929       &  kjpindex,nvm,nlitt
2930       STOP 'stomate_init'
2931    ENDIF
2932
2933    ALLOCATE(litter(kjpindex,nlitt,nvm,nlevs,nelements),stat=ier)
2934    l_error = l_error .OR. (ier /= 0)
2935    IF (l_error) THEN
2936       WRITE(numout,*) 'Memory allocation error for litter. We stop. We need kjpindex*nlitt*nvm*nlevs*nelements words', & 
2937       &    kjpindex,nlitt,nvm,nlevs,nelements
2938       STOP 'stomate_init'
2939    ENDIF
2940
2941    ALLOCATE(dead_leaves(kjpindex,nvm,nlitt),stat=ier)
2942    l_error = l_error .OR. (ier /= 0)
2943    IF (l_error) THEN
2944       WRITE(numout,*) 'Memory allocation error for dead_leaves. We stop. We need kjpindex*nvm*nlitt words', & 
2945       &   kjpindex,nvm,nlitt
2946       STOP 'stomate_init'
2947    ENDIF
2948
2949    ALLOCATE(carbon(kjpindex,ncarb,nvm),stat=ier)
2950    l_error = l_error .OR. (ier /= 0)
2951    IF (l_error) THEN
2952       WRITE(numout,*) 'Memory allocation error for carbon. We stop. We need kjpindex*ncarb*nvm words',kjpindex,ncarb,nvm
2953       STOP 'stomate_init'
2954    ENDIF
2955
2956    ALLOCATE(lignin_struc(kjpindex,nvm,nlevs),stat=ier)
2957    l_error = l_error .OR. (ier /= 0)
2958    IF (l_error) THEN
2959       WRITE(numout,*) 'Memory allocation error for lignin_struc. We stop. We need kjpindex*nvm*nlevs words',kjpindex,nvm,nlevs
2960       STOP 'stomate_init'
2961    ENDIF
2962
2963    ALLOCATE(turnover_time(kjpindex,nvm),stat=ier)
2964    l_error = l_error .OR. (ier /= 0)
2965    IF (l_error) THEN
2966       WRITE(numout,*) 'Memory allocation error for turnover_time. We stop. We need kjpindex*nvm words',kjpindex,nvm
2967       STOP 'stomate_init'
2968    ENDIF
2969
2970    ALLOCATE(nep_daily(kjpindex,nvm),stat=ier)
2971    l_error = l_error .OR. (ier /= 0)
2972    IF (l_error) THEN
2973       WRITE(numout,*) 'Memory allocation error for nep_daily. We stop. We need kjpindex*nvm words',kjpindex,nvm
2974       STOP 'stomate_init'
2975    ENDIF
2976
2977    ALLOCATE(nep_monthly(kjpindex,nvm),stat=ier)
2978    l_error = l_error .OR. (ier /= 0)
2979    IF (l_error) THEN
2980       WRITE(numout,*) 'Memory allocation error for nep_monthly. We stop. We need kjpindex*nvm words',kjpindex,nvm
2981       STOP 'stomate_init'
2982    ENDIF
2983
2984    ALLOCATE (cflux_prod_monthly(kjpindex), stat=ier)
2985    l_error = l_error .OR. (ier /= 0)
2986    IF (l_error) THEN
2987       WRITE(numout,*) 'Memory allocation error for cflux_prod_monthly. We stop. We need kjpindex words',kjpindex
2988       STOP 'stomate_init'
2989    ENDIF
2990 
2991    ALLOCATE (harvest_above_monthly(kjpindex), stat=ier)
2992    l_error = l_error .OR. (ier /= 0)
2993    IF (l_error) THEN
2994       WRITE(numout,*) 'Memory allocation error for harvest_above_monthly. We stop. We need kjpindex words',kjpindex
2995       STOP 'stomate_init'
2996    ENDIF
2997
2998    ALLOCATE(bm_to_litter(kjpindex,nvm,nparts,nelements),stat=ier)
2999    l_error = l_error .OR. (ier /= 0)
3000    IF (l_error) THEN
3001       WRITE(numout,*) 'Memory allocation error for bm_to_litter. We stop. We need kjpindex*nvm*nparts*nelements words', & 
3002       &    kjpindex,nvm,nparts,nelements
3003       STOP 'stomate_init'
3004    ENDIF
3005
3006    ALLOCATE(bm_to_littercalc(kjpindex,nvm,nparts,nelements),stat=ier)
3007    l_error = l_error .OR. (ier /= 0)
3008    IF (l_error) THEN
3009       WRITE(numout,*) 'Memory allocation error for bm_to_littercalc. We stop. We need kjpindex*nvm*nparts*nelements words', &
3010       &   kjpindex,nvm,nparts,nelements
3011       STOP 'stomate_init'
3012    ENDIF
3013
3014    ALLOCATE(herbivores(kjpindex,nvm),stat=ier)
3015    l_error = l_error .OR. (ier /= 0)
3016    IF (l_error) THEN
3017       WRITE(numout,*) 'Memory allocation error for herbivores. We stop. We need kjpindex*nvm words',kjpindex,nvm
3018       STOP 'stomate_init'
3019    ENDIF
3020
3021    ALLOCATE(hori_index(kjpindex),stat=ier)
3022    l_error = l_error .OR. (ier /= 0)
3023    IF (l_error) THEN
3024       WRITE(numout,*) 'Memory allocation error for hori_index. We stop. We need kjpindex words',kjpindex
3025       STOP 'stomate_init'
3026    ENDIF
3027
3028    ALLOCATE(horipft_index(kjpindex*nvm),stat=ier)
3029    l_error = l_error .OR. (ier /= 0)
3030    IF (l_error) THEN
3031       WRITE(numout,*) 'Memory allocation error for horipft_index. We stop. We need kjpindex*nvm words',kjpindex*nvm
3032       STOP 'stomate_init'
3033    ENDIF
3034
3035    ALLOCATE(resp_maint_part_radia(kjpindex,nvm,nparts),stat=ier)
3036    l_error = l_error .OR. (ier /= 0)
3037    IF (l_error) THEN
3038       WRITE(numout,*) 'Memory allocation error for resp_maint_part_radia. We stop. We need kjpindex*nvm*nparts words', &
3039       &  kjpindex,nvm,nparts
3040       STOP 'stomate_init'
3041    ENDIF
3042
3043    ALLOCATE(resp_maint_radia(kjpindex,nvm),stat=ier)
3044    l_error = l_error .OR. (ier /= 0)
3045    IF (l_error) THEN
3046       WRITE(numout,*) 'Memory allocation error for resp_maint_radia. We stop. We need kjpindex*nvm words',kjpindex,nvm
3047       STOP 'stomate_init'
3048    ENDIF
3049
3050    ALLOCATE(resp_maint_part(kjpindex,nvm,nparts),stat=ier)
3051    l_error = l_error .OR. (ier /= 0)
3052    IF (l_error) THEN
3053       WRITE(numout,*) 'Memory allocation error for resp_maint_part. We stop. We need kjpindex*nvm*nparts words', &
3054       &    kjpindex,nvm,nparts
3055       STOP 'stomate_init'
3056    ENDIF
3057    resp_maint_part(:,:,:) = zero
3058
3059    ALLOCATE (horip10_index(kjpindex*10), stat=ier)
3060    l_error = l_error .OR. (ier /= 0)
3061    IF (l_error) THEN
3062       WRITE(numout,*) 'Memory allocation error for horip10_index. We stop. We need kjpindex*10 words',kjpindex,10
3063       STOP 'stomate_init'
3064    ENDIF
3065
3066    ALLOCATE (horip100_index(kjpindex*100), stat=ier)
3067    l_error = l_error .OR. (ier /= 0)
3068    IF (l_error) THEN
3069       WRITE(numout,*) 'Memory allocation error for horip100_index. We stop. We need kjpindex*100 words',kjpindex,100
3070       STOP 'stomate_init'
3071    ENDIF
3072
3073    ALLOCATE (horip11_index(kjpindex*11), stat=ier)
3074    l_error = l_error .OR. (ier /= 0)
3075    IF (l_error) THEN
3076       WRITE(numout,*) 'Memory allocation error for horip11_index. We stop. We need kjpindex*11 words',kjpindex,11
3077       STOP 'stomate_init'
3078    ENDIF
3079
3080    ALLOCATE (horip101_index(kjpindex*101), stat=ier)
3081    l_error = l_error .OR. (ier /= 0)
3082    IF (l_error) THEN
3083       WRITE(numout,*) 'Memory allocation error for horip101_index. We stop. We need kjpindex*101 words',kjpindex,101
3084       STOP 'stomate_init'
3085    ENDIF
3086
3087    ALLOCATE (co2_flux(kjpindex,nvm), stat=ier)
3088    l_error = l_error .OR. (ier /= 0)
3089    IF (l_error) THEN
3090       WRITE(numout,*) 'Memory allocation error for co2_flux. We stop. We need kjpindex words',kjpindex,nvm
3091       STOP 'stomate_init'
3092    ENDIF
3093
3094    ALLOCATE (fco2_lu(kjpindex), stat=ier)
3095    l_error = l_error .OR. (ier /= 0)
3096    IF (l_error) THEN
3097       WRITE(numout,*) 'Memory allocation error for fco2_lu. We stop. We need kjpindex words',kjpindex
3098       STOP 'stomate_init'
3099    ENDIF
3100
3101    ALLOCATE (fco2_wh(kjpindex), stat=ier)
3102    l_error = l_error .OR. (ier /= 0)
3103    IF (l_error) THEN
3104       WRITE(numout,*) 'Memory allocation error for fco2_wh. We stop. We need kjpindex words',kjpindex
3105       STOP 'stomate_init'
3106    ENDIF
3107
3108    ALLOCATE (fco2_ha(kjpindex), stat=ier)
3109    l_error = l_error .OR. (ier /= 0)
3110    IF (l_error) THEN
3111       WRITE(numout,*) 'Memory allocation error for fco2_ha. We stop. We need kjpindex words',kjpindex
3112       STOP 'stomate_init'
3113    ENDIF
3114
3115    ALLOCATE (prod10(kjpindex,0:10), stat=ier)
3116    l_error = l_error .OR. (ier /= 0)
3117    IF (l_error) THEN
3118       WRITE(numout,*) 'Memory allocation error for prod10. We stop. We need kjpindex*11 words',kjpindex,11
3119       STOP 'stomate_init'
3120    ENDIF
3121
3122    ALLOCATE (prod100(kjpindex,0:100), stat=ier)
3123    l_error = l_error .OR. (ier /= 0)
3124    IF (l_error) THEN
3125       WRITE(numout,*) 'Memory allocation error for prod100. We stop. We need kjpindex*101 words',kjpindex,101
3126       STOP 'stomate_init'
3127    ENDIF
3128
3129    ALLOCATE (flux10(kjpindex,10), stat=ier)
3130    l_error = l_error .OR. (ier /= 0)
3131    IF (l_error) THEN
3132       WRITE(numout,*) 'Memory allocation error for flux10. We stop. We need kjpindex*10 words',kjpindex,10
3133       STOP 'stomate_init'
3134    ENDIF
3135
3136    ALLOCATE (flux100(kjpindex,100), stat=ier)
3137    l_error = l_error .OR. (ier /= 0)
3138    IF (l_error) THEN
3139       WRITE(numout,*) 'Memory allocation error for flux100. We stop. We need kjpindex*100 words',kjpindex,100
3140       STOP 'stomate_init'
3141    ENDIF
3142
3143    ALLOCATE (convflux(kjpindex), stat=ier)
3144    l_error = l_error .OR. (ier /= 0)
3145    IF (l_error) THEN
3146       WRITE(numout,*) 'Memory allocation error for convflux. We stop. We need kjpindex words',kjpindex
3147       STOP 'stomate_init'
3148    ENDIF
3149
3150    ALLOCATE (cflux_prod10(kjpindex), stat=ier)
3151    l_error = l_error .OR. (ier /= 0)
3152    IF (l_error) THEN
3153       WRITE(numout,*) 'Memory allocation error for cflux_prod10. We stop. We need kjpindex words',kjpindex
3154       STOP 'stomate_init'
3155    ENDIF
3156
3157    ALLOCATE (cflux_prod100(kjpindex), stat=ier)
3158    l_error = l_error .OR. (ier /= 0)
3159    IF (l_error) THEN
3160       WRITE(numout,*) 'Memory allocation error for cflux_prod100. We stop. We need kjpindex words',kjpindex
3161       STOP 'stomate_init'
3162    ENDIF
3163
3164    ALLOCATE (prod10_harvest(kjpindex,0:10), stat=ier)
3165    l_error = l_error .OR. (ier /= 0)
3166    IF (l_error) THEN
3167       WRITE(numout,*) 'Memory allocation error for prod10_harvest. We stop. We need kjpindex*11 words',kjpindex,11
3168       STOP 'stomate_init'
3169    ENDIF
3170
3171    ALLOCATE (prod100_harvest(kjpindex,0:100), stat=ier)
3172    l_error = l_error .OR. (ier /= 0)
3173    IF (l_error) THEN
3174       WRITE(numout,*) 'Memory allocation error for prod100_harvest. We stop. We need kjpindex*101 words',kjpindex,101
3175       STOP 'stomate_init'
3176    ENDIF
3177
3178    ALLOCATE (flux10_harvest(kjpindex,10), stat=ier)
3179    l_error = l_error .OR. (ier /= 0)
3180    IF (l_error) THEN
3181       WRITE(numout,*) 'Memory allocation error for flux10_harvest. We stop. We need kjpindex*10 words',kjpindex,10
3182       STOP 'stomate_init'
3183    ENDIF
3184
3185    ALLOCATE (flux100_harvest(kjpindex,100), stat=ier)
3186    l_error = l_error .OR. (ier /= 0)
3187    IF (l_error) THEN
3188       WRITE(numout,*) 'Memory allocation error for flux100_harvest. We stop. We need kjpindex*100 words',kjpindex,100
3189       STOP 'stomate_init'
3190    ENDIF
3191
3192    ALLOCATE (convflux_harvest(kjpindex), stat=ier)
3193    l_error = l_error .OR. (ier /= 0)
3194    IF (l_error) THEN
3195       WRITE(numout,*) 'Memory allocation error for convflux_harvest. We stop. We need kjpindex words',kjpindex
3196       STOP 'stomate_init'
3197    ENDIF
3198
3199    ALLOCATE (cflux_prod10_harvest(kjpindex), stat=ier)
3200    l_error = l_error .OR. (ier /= 0)
3201    IF (l_error) THEN
3202       WRITE(numout,*) 'Memory allocation error for cflux_prod10_harvest. We stop. We need kjpindex words',kjpindex
3203       STOP 'stomate_init'
3204    ENDIF
3205
3206    ALLOCATE (cflux_prod100_harvest(kjpindex), stat=ier)
3207    l_error = l_error .OR. (ier /= 0)
3208    IF (l_error) THEN
3209       WRITE(numout,*) 'Memory allocation error for cflux_prod100_harvest. We stop. We need kjpindex words',kjpindex
3210       STOP 'stomate_init'
3211    ENDIF
3212
3213    ALLOCATE (woodharvestpft(kjpindex,nvm), stat=ier)
3214    l_error = l_error .OR. (ier /= 0)
3215    IF (l_error) THEN
3216       WRITE(numout,*) 'Memory allocation error for woodharvestpft. We stop. We need kjpindex*nvm words',kjpindex*nvm
3217       STOP 'stomate_init'
3218    ENDIF
3219
3220    ALLOCATE (convfluxpft(kjpindex,nvm), stat=ier)
3221    l_error = l_error .OR. (ier /= 0)
3222    IF (l_error) THEN
3223       WRITE(numout,*) 'Memory allocation error for convfluxpft. We stop. We need kjpindex*nvm words',kjpindex*nvm
3224       STOP 'stomate_init'
3225    ENDIF
3226
3227    ALLOCATE (fDeforestToProduct(kjpindex,nvm), stat=ier)
3228    l_error = l_error .OR. (ier /= 0)
3229    IF (l_error) THEN
3230       WRITE(numout,*) 'Memory allocation error for fDeforestToProduct. We stop. We need kjpindex*nvm words',kjpindex*nvm
3231       STOP 'stomate_init'
3232    ENDIF
3233
3234    ALLOCATE (fLulccResidue(kjpindex,nvm), stat=ier)
3235    l_error = l_error .OR. (ier /= 0)
3236    IF (l_error) THEN
3237       WRITE(numout,*) 'Memory allocation error for fLulccResidue. We stop. We need kjpindex*nvm words',kjpindex*nvm
3238       STOP 'stomate_init'
3239    ENDIF
3240
3241    ALLOCATE (fHarvestToProduct(kjpindex,nvm), stat=ier)
3242    l_error = l_error .OR. (ier /= 0)
3243    IF (l_error) THEN
3244       WRITE(numout,*) 'Memory allocation error for fHarvestToProduct. We stop. We need kjpindex*nvm words',kjpindex*nvm
3245       STOP 'stomate_init'
3246    ENDIF
3247
3248    ALLOCATE (harvest_above(kjpindex), stat=ier)
3249    l_error = l_error .OR. (ier /= 0)
3250    IF (l_error) THEN
3251       WRITE(numout,*) 'Memory allocation error for harvest_above. We stop. We need kjpindex words',kjpindex
3252       STOP 'stomate_init'
3253    ENDIF
3254
3255    ALLOCATE (carb_mass_total(kjpindex), stat=ier)
3256    l_error = l_error .OR. (ier /= 0)
3257    IF (l_error) THEN
3258       WRITE(numout,*) 'Memory allocation error for carb_mass_total. We stop. We need kjpindex words',kjpindex
3259       STOP 'stomate_init'
3260    ENDIF
3261
3262    ALLOCATE (soilcarbon_input_daily(kjpindex,ncarb,nvm), stat=ier)
3263    l_error = l_error .OR. (ier /= 0)
3264    IF (l_error) THEN
3265       WRITE(numout,*) 'Memory allocation error for soilcarbon_input_daily. We stop. We need kjpindex*ncarb*nvm words', & 
3266       &    kjpindex,ncarb,nvm
3267       STOP 'stomate_init'
3268    ENDIF
3269
3270    ALLOCATE (control_temp_daily(kjpindex,nlevs), stat=ier)
3271    l_error = l_error .OR. (ier /= 0)
3272    IF (l_error) THEN
3273       WRITE(numout,*) 'Memory allocation error for control_temp_daily. We stop. We need kjpindex*nlevs words',kjpindex,nlevs
3274       STOP 'stomate_init'
3275    ENDIF
3276
3277    ALLOCATE (control_moist_daily(kjpindex,nlevs), stat=ier)
3278    l_error = l_error .OR. (ier /= 0)
3279    IF (l_error) THEN
3280       WRITE(numout,*) 'Memory allocation error for control_moist_daily. We stop. We need kjpindex*nlevs words',kjpindex,nlevs
3281       STOP 'stomate_init'
3282    ENDIF
3283
3284    ALLOCATE (fpc_max(kjpindex,nvm), stat=ier)
3285    l_error = l_error .OR. (ier /= 0)
3286    IF (l_error) THEN
3287       WRITE(numout,*) 'Memory allocation error for fpc_max. We stop. We need kjpindex*nvm words',kjpindex,nvm
3288       STOP 'stomate_init'
3289    ENDIF
3290
3291    ALLOCATE(ok_equilibrium(kjpindex),stat=ier)
3292    l_error = l_error .OR. (ier /= 0) 
3293    IF (l_error) THEN
3294       WRITE(numout,*) 'Memory allocation error for ok_equilibrium. We stop. We need kjpindex words',kjpindex
3295       STOP 'stomate_init'
3296    ENDIF
3297
3298    ALLOCATE(carbon_eq(kjpindex),stat=ier)
3299    l_error = l_error .OR. (ier /= 0)
3300    IF (l_error) THEN
3301       WRITE(numout,*) 'Memory allocation error for carbon_eq. We stop. We need kjpindex words',kjpindex
3302       STOP 'stomate_init'
3303    ENDIF
3304
3305    ALLOCATE(nbp_accu(kjpindex),stat=ier)
3306    l_error = l_error .OR. (ier /= 0)
3307    IF (l_error) THEN
3308       WRITE(numout,*) 'Memory allocation error for nbp_accu. We stop. We need kjpindex words',kjpindex
3309       STOP 'stomate_init'
3310    ENDIF
3311
3312    ALLOCATE(nbp_flux(kjpindex),stat=ier)
3313    l_error = l_error .OR. (ier /= 0)
3314    IF (l_error) THEN
3315       WRITE(numout,*) 'Memory allocation error for nbp_flux. We stop. We need kjpindex words',kjpindex
3316       STOP 'stomate_init'
3317    ENDIF
3318
3319    ALLOCATE(matrixA(kjpindex,nvm,nbpools,nbpools),stat=ier)
3320    l_error = l_error .OR. (ier /= 0)
3321    IF (l_error) THEN
3322       WRITE(numout,*) 'Memory allocation error for matrixA. We stop. We need kjpindex*nvm*nbpools*nbpools words',  & 
3323       &     kjpindex, nvm, nbpools, nbpools
3324       STOP 'stomate_init'
3325    ENDIF
3326
3327    ALLOCATE(vectorB(kjpindex,nvm,nbpools),stat=ier)
3328    l_error = l_error .OR. (ier /= 0)
3329    IF (l_error) THEN
3330       WRITE(numout,*) 'Memory allocation error for vectorB. We stop. We need kjpindex*nvm*nbpools words',  & 
3331       &     kjpindex, nvm, nbpools
3332       STOP 'stomate_init'
3333    ENDIF
3334
3335    ALLOCATE(VectorU(kjpindex,nvm,nbpools),stat=ier)
3336    l_error = l_error .OR. (ier /= 0)
3337    IF (l_error) THEN
3338       WRITE(numout,*) 'Memory allocation error for VectorU. We stop. We need kjpindex*nvm*nbpools words',  & 
3339       &     kjpindex, nvm, nbpools
3340       STOP 'stomate_init'
3341    ENDIF
3342
3343    ALLOCATE(MatrixV(kjpindex,nvm,nbpools,nbpools),stat=ier)
3344    l_error = l_error .OR. (ier /= 0)
3345    IF (l_error) THEN
3346       WRITE(numout,*) 'Memory allocation error for MatrixV. We stop. We need kjpindex*nvm*nbpools*nbpools words',  & 
3347       &     kjpindex, nvm, nbpools, nbpools
3348       STOP 'stomate_init'
3349    ENDIF
3350
3351    ALLOCATE(MatrixW(kjpindex,nvm,nbpools,nbpools),stat=ier)
3352    l_error = l_error .OR. (ier /= 0)
3353    IF (l_error) THEN
3354       WRITE(numout,*) 'Memory allocation error for MatrixW. We stop. We need kjpindex*nvm*nbpools*nbpools words',  & 
3355       &     kjpindex, nvm, nbpools, nbpools
3356       STOP 'stomate_init'
3357    ENDIF
3358
3359    ALLOCATE(previous_stock(kjpindex,nvm,nbpools),stat=ier)
3360    l_error = l_error .OR. (ier /= 0)
3361    IF (l_error) THEN
3362       WRITE(numout,*) 'Memory allocation error for previous_stock. We stop. We need kjpindex*nvm*nbpools words',  & 
3363       &     kjpindex, nvm, nbpools
3364       STOP 'stomate_init'
3365    ENDIF
3366
3367    ALLOCATE(current_stock(kjpindex,nvm,nbpools),stat=ier)
3368    l_error = l_error .OR. (ier /= 0)
3369    IF (l_error) THEN
3370       WRITE(numout,*) 'Memory allocation error for current_stock. We stop. We need kjpindex*nvm*nbpools words',  & 
3371       &     kjpindex, nvm, nbpools
3372       STOP 'stomate_init'
3373    ENDIF
3374   
3375  !! 5. File definitions
3376
3377    ! Store history and restart files in common variables
3378    hist_id_stomate = hist_id_stom
3379    hist_id_stomate_IPCC = hist_id_stom_IPCC
3380    rest_id_stomate = rest_id_stom
3381   
3382    ! In STOMATE reduced grids are used containing only terrestrial pixels.
3383    ! Build a new indexing table for the vegetation fields separating
3384    ! between the different PFTs. Note that ::index has dimension (kjpindex)
3385    ! wheras ::indexpft has dimension (kjpindex*nvm).
3386
3387    hori_index(:) = index(:)
3388
3389    DO j = 1, nvm
3390       DO ji = 1, kjpindex
3391          horipft_index((j-1)*kjpindex+ji) = index(ji)+(j-1)*kjpij + offset_omp - offset_mpi
3392       ENDDO
3393    ENDDO
3394
3395    ! Similar index tables are build for the land cover change variables
3396    DO j = 1, 10
3397       DO ji = 1, kjpindex
3398          horip10_index((j-1)*kjpindex+ji) = index(ji)+(j-1)*kjpij + offset_omp - offset_mpi
3399       ENDDO
3400    ENDDO
3401
3402    DO j = 1, 100
3403       DO ji = 1, kjpindex
3404          horip100_index((j-1)*kjpindex+ji) = index(ji)+(j-1)*kjpij + offset_omp - offset_mpi
3405       ENDDO
3406    ENDDO
3407
3408    DO j = 1, 11
3409       DO ji = 1, kjpindex
3410          horip11_index((j-1)*kjpindex+ji) = index(ji)+(j-1)*kjpij + offset_omp - offset_mpi
3411       ENDDO
3412    ENDDO
3413
3414    DO j = 1, 101
3415       DO ji = 1, kjpindex
3416          horip101_index((j-1)*kjpindex+ji) = index(ji)+(j-1)*kjpij + offset_omp - offset_mpi
3417       ENDDO
3418    ENDDO
3419
3420  !! 6. Initialization of global and land cover change variables.
3421
3422    ! All variables are cumulative variables. bm_to_litter is not and is therefore
3423    ! excluded
3424    !   bm_to_litter(:,:,:) = zero
3425    turnover_daily(:,:,:,:) = zero
3426    resp_hetero_d(:,:) = zero
3427    nep_daily(:,:) = zero
3428    nep_monthly(:,:) = zero
3429    cflux_prod_monthly(:) = zero
3430    harvest_above_monthly(:) = zero
3431    control_moist_daily(:,:) = zero
3432    control_temp_daily(:,:) = zero
3433    soilcarbon_input_daily(:,:,:) = zero
3434    ! Land cover change variables
3435    prod10(:,:)  = zero
3436    prod100(:,:) = zero
3437    flux10(:,:)  = zero
3438    flux100(:,:) = zero
3439    convflux(:)  = zero
3440    cflux_prod10(:) = zero
3441    cflux_prod100(:) = zero
3442    prod10_harvest(:,:)  = zero
3443    prod100_harvest(:,:) = zero
3444    flux10_harvest(:,:)  = zero
3445    flux100_harvest(:,:) = zero
3446    convflux_harvest(:)  = zero
3447    cflux_prod10_harvest(:) = zero
3448    cflux_prod100_harvest(:) = zero
3449    woodharvestpft(:,:) = zero
3450    fpc_max(:,:)=zero
3451   
3452    convfluxpft(:,:)=zero
3453    fDeforestToProduct(:,:)=zero
3454    fLulccResidue(:,:)=zero
3455    fHarvestToProduct(:,:)=zero
3456  END SUBROUTINE stomate_init
3457
3458
3459!! ================================================================================================================================
3460!! SUBROUTINE   : stomate_clear
3461!!
3462!>\BRIEF        Deallocate memory of the stomate variables.
3463!!
3464!! DESCRIPTION  : None
3465!!
3466!! RECENT CHANGE(S) : None
3467!!
3468!! MAIN OUTPUT VARIABLE(S): None
3469!!
3470!! REFERENCES   : None
3471!!
3472!! FLOWCHART    : None
3473!! \n
3474!_ ================================================================================================================================
3475 
3476  SUBROUTINE stomate_clear
3477
3478  !! 1. Deallocate all dynamics variables
3479
3480    IF (ALLOCATED(veget_cov_max)) DEALLOCATE(veget_cov_max)
3481    IF (ALLOCATED(ind)) DEALLOCATE(ind)
3482    IF (ALLOCATED(adapted)) DEALLOCATE(adapted)
3483    IF (ALLOCATED(regenerate)) DEALLOCATE(regenerate)
3484    IF (ALLOCATED(humrel_daily)) DEALLOCATE(humrel_daily)
3485    IF (ALLOCATED(gdd_init_date)) DEALLOCATE(gdd_init_date)
3486    IF (ALLOCATED(litterhum_daily)) DEALLOCATE(litterhum_daily)
3487    IF (ALLOCATED(t2m_daily))  DEALLOCATE(t2m_daily)
3488    IF (ALLOCATED(t2m_min_daily))  DEALLOCATE(t2m_min_daily)
3489    IF (ALLOCATED(tsurf_daily))  DEALLOCATE(tsurf_daily)
3490    IF (ALLOCATED(tsoil_daily)) DEALLOCATE(tsoil_daily)
3491    IF (ALLOCATED(soilhum_daily)) DEALLOCATE(soilhum_daily)
3492    IF (ALLOCATED(precip_daily)) DEALLOCATE(precip_daily)
3493    IF (ALLOCATED(gpp_daily)) DEALLOCATE(gpp_daily)
3494    IF (ALLOCATED(npp_daily)) DEALLOCATE(npp_daily)
3495    IF (ALLOCATED(turnover_daily)) DEALLOCATE(turnover_daily)
3496    IF (ALLOCATED(turnover_littercalc)) DEALLOCATE(turnover_littercalc)
3497    IF (ALLOCATED(humrel_month)) DEALLOCATE(humrel_month)
3498    IF (ALLOCATED(humrel_week)) DEALLOCATE(humrel_week)
3499    IF (ALLOCATED(t2m_longterm)) DEALLOCATE(t2m_longterm)
3500    IF (ALLOCATED(t2m_month)) DEALLOCATE(t2m_month)
3501    IF (ALLOCATED(Tseason)) DEALLOCATE(Tseason)
3502    IF (ALLOCATED(Tseason_length)) DEALLOCATE(Tseason_length)
3503    IF (ALLOCATED(Tseason_tmp)) DEALLOCATE(Tseason_tmp)
3504    IF (ALLOCATED(Tmin_spring_time)) DEALLOCATE(Tmin_spring_time)
3505    IF (ALLOCATED(onset_date)) DEALLOCATE(onset_date)
3506    IF (ALLOCATED(begin_leaves)) DEALLOCATE(begin_leaves)
3507    IF (ALLOCATED(t2m_week)) DEALLOCATE(t2m_week)
3508    IF (ALLOCATED(tsoil_month)) DEALLOCATE(tsoil_month)
3509    IF (ALLOCATED(soilhum_month)) DEALLOCATE(soilhum_month)
3510    IF (ALLOCATED(fireindex)) DEALLOCATE(fireindex)
3511    IF (ALLOCATED(firelitter)) DEALLOCATE(firelitter)
3512    IF (ALLOCATED(maxhumrel_lastyear)) DEALLOCATE(maxhumrel_lastyear)
3513    IF (ALLOCATED(maxhumrel_thisyear)) DEALLOCATE(maxhumrel_thisyear)
3514    IF (ALLOCATED(minhumrel_lastyear)) DEALLOCATE(minhumrel_lastyear)
3515    IF (ALLOCATED(minhumrel_thisyear)) DEALLOCATE(minhumrel_thisyear)
3516    IF (ALLOCATED(maxgppweek_lastyear)) DEALLOCATE(maxgppweek_lastyear)
3517    IF (ALLOCATED(maxgppweek_thisyear)) DEALLOCATE(maxgppweek_thisyear)
3518    IF (ALLOCATED(gdd0_lastyear)) DEALLOCATE(gdd0_lastyear)
3519    IF (ALLOCATED(gdd0_thisyear)) DEALLOCATE(gdd0_thisyear)
3520    IF (ALLOCATED(precip_lastyear)) DEALLOCATE(precip_lastyear)
3521    IF (ALLOCATED(precip_thisyear)) DEALLOCATE(precip_thisyear)
3522    IF (ALLOCATED(gdd_m5_dormance)) DEALLOCATE(gdd_m5_dormance)
3523    IF (ALLOCATED(gdd_from_growthinit)) DEALLOCATE(gdd_from_growthinit)
3524    IF (ALLOCATED(gdd_midwinter)) DEALLOCATE(gdd_midwinter)
3525    IF (ALLOCATED(ncd_dormance)) DEALLOCATE(ncd_dormance)
3526    IF (ALLOCATED(ngd_minus5))  DEALLOCATE(ngd_minus5)
3527    IF (ALLOCATED(PFTpresent)) DEALLOCATE(PFTpresent)
3528    IF (ALLOCATED(npp_longterm)) DEALLOCATE(npp_longterm)
3529    IF (ALLOCATED(lm_lastyearmax)) DEALLOCATE(lm_lastyearmax)
3530    IF (ALLOCATED(lm_thisyearmax)) DEALLOCATE(lm_thisyearmax)
3531    IF (ALLOCATED(maxfpc_lastyear)) DEALLOCATE(maxfpc_lastyear)
3532    IF (ALLOCATED(maxfpc_thisyear)) DEALLOCATE(maxfpc_thisyear)
3533    IF (ALLOCATED(turnover_longterm)) DEALLOCATE(turnover_longterm)
3534    IF (ALLOCATED(gpp_week)) DEALLOCATE(gpp_week)
3535    IF (ALLOCATED(biomass)) DEALLOCATE(biomass)
3536    IF (ALLOCATED(senescence)) DEALLOCATE(senescence)
3537    IF (ALLOCATED(when_growthinit)) DEALLOCATE(when_growthinit)
3538    IF (ALLOCATED(age))  DEALLOCATE(age)
3539    IF (ALLOCATED(resp_hetero_d)) DEALLOCATE(resp_hetero_d)
3540    IF (ALLOCATED(resp_hetero_radia)) DEALLOCATE(resp_hetero_radia)
3541    IF (ALLOCATED(resp_maint_d)) DEALLOCATE(resp_maint_d)
3542    IF (ALLOCATED(resp_growth_d)) DEALLOCATE(resp_growth_d)
3543    IF (ALLOCATED(co2_fire)) DEALLOCATE(co2_fire)
3544    IF (ALLOCATED(co2_to_bm_dgvm)) DEALLOCATE(co2_to_bm_dgvm)
3545    IF (ALLOCATED(veget_lastlight)) DEALLOCATE(veget_lastlight)
3546    IF (ALLOCATED(everywhere)) DEALLOCATE(everywhere)
3547    IF (ALLOCATED(need_adjacent)) DEALLOCATE(need_adjacent)
3548    IF (ALLOCATED(leaf_age)) DEALLOCATE(leaf_age)
3549    IF (ALLOCATED(leaf_frac)) DEALLOCATE(leaf_frac)
3550    IF (ALLOCATED(RIP_time)) DEALLOCATE(RIP_time)
3551    IF (ALLOCATED(time_hum_min)) DEALLOCATE(time_hum_min)
3552    IF (ALLOCATED(hum_min_dormance)) DEALLOCATE(hum_min_dormance)
3553    IF (ALLOCATED(litterpart)) DEALLOCATE(litterpart)
3554    IF (ALLOCATED(litter)) DEALLOCATE(litter)
3555    IF (ALLOCATED(dead_leaves)) DEALLOCATE(dead_leaves)
3556    IF (ALLOCATED(carbon)) DEALLOCATE(carbon)
3557    IF (ALLOCATED(lignin_struc)) DEALLOCATE(lignin_struc)
3558    IF (ALLOCATED(turnover_time)) DEALLOCATE(turnover_time)
3559    IF (ALLOCATED(nep_daily)) DEALLOCATE(nep_daily)
3560    IF (ALLOCATED(nep_monthly)) DEALLOCATE(nep_monthly)
3561    IF (ALLOCATED(harvest_above_monthly)) DEALLOCATE (harvest_above_monthly)
3562    IF (ALLOCATED(cflux_prod_monthly)) DEALLOCATE (cflux_prod_monthly)
3563    IF (ALLOCATED(bm_to_litter)) DEALLOCATE(bm_to_litter)
3564    IF (ALLOCATED(bm_to_littercalc)) DEALLOCATE(bm_to_littercalc)
3565    IF (ALLOCATED(herbivores)) DEALLOCATE(herbivores)
3566    IF (ALLOCATED(resp_maint_part_radia)) DEALLOCATE(resp_maint_part_radia)
3567    IF (ALLOCATED(resp_maint_radia)) DEALLOCATE(resp_maint_radia)
3568    IF (ALLOCATED(resp_maint_part)) DEALLOCATE(resp_maint_part)
3569    IF (ALLOCATED(hori_index)) DEALLOCATE(hori_index)
3570    IF (ALLOCATED(horipft_index)) DEALLOCATE(horipft_index)
3571    IF (ALLOCATED(clay_fm)) DEALLOCATE(clay_fm)
3572    IF (ALLOCATED(humrel_daily_fm)) DEALLOCATE(humrel_daily_fm)
3573    IF (ALLOCATED(litterhum_daily_fm))  DEALLOCATE(litterhum_daily_fm)
3574    IF (ALLOCATED(t2m_daily_fm))  DEALLOCATE(t2m_daily_fm)
3575    IF (ALLOCATED(t2m_min_daily_fm))  DEALLOCATE(t2m_min_daily_fm)
3576    IF (ALLOCATED(tsurf_daily_fm)) DEALLOCATE(tsurf_daily_fm)
3577    IF (ALLOCATED(tsoil_daily_fm)) DEALLOCATE(tsoil_daily_fm)
3578    IF (ALLOCATED(soilhum_daily_fm))  DEALLOCATE(soilhum_daily_fm)
3579    IF (ALLOCATED(precip_fm)) DEALLOCATE(precip_fm)
3580    IF (ALLOCATED(gpp_daily_fm))  DEALLOCATE(gpp_daily_fm)
3581    IF (ALLOCATED(veget_fm)) DEALLOCATE(veget_fm)
3582    IF (ALLOCATED(veget_max_fm)) DEALLOCATE(veget_max_fm)
3583    IF (ALLOCATED(lai_fm))  DEALLOCATE(lai_fm)
3584    !
3585    IF (ALLOCATED(ok_equilibrium)) DEALLOCATE(ok_equilibrium)
3586    IF (ALLOCATED(carbon_eq)) DEALLOCATE(carbon_eq)
3587    IF (ALLOCATED(matrixA)) DEALLOCATE(matrixA)
3588    IF (ALLOCATED(vectorB)) DEALLOCATE(vectorB)
3589    IF (ALLOCATED(MatrixV)) DEALLOCATE(MatrixV)
3590    IF (ALLOCATED(VectorU)) DEALLOCATE(VectorU)
3591    IF (ALLOCATED(MatrixW)) DEALLOCATE(MatrixW)
3592    IF (ALLOCATED(previous_stock)) DEALLOCATE(previous_stock)
3593    IF (ALLOCATED(current_stock)) DEALLOCATE(current_stock) 
3594    IF (ALLOCATED(nbp_accu)) DEALLOCATE(nbp_accu)
3595    IF (ALLOCATED(nbp_flux)) DEALLOCATE(nbp_flux)
3596
3597    IF (ALLOCATED(clay_fm_g)) DEALLOCATE(clay_fm_g)
3598    IF (ALLOCATED(humrel_daily_fm_g)) DEALLOCATE(humrel_daily_fm_g)
3599    IF (ALLOCATED(litterhum_daily_fm_g))  DEALLOCATE(litterhum_daily_fm_g)
3600    IF (ALLOCATED(t2m_daily_fm_g))  DEALLOCATE(t2m_daily_fm_g)
3601    IF (ALLOCATED(t2m_min_daily_fm_g))  DEALLOCATE(t2m_min_daily_fm_g)
3602    IF (ALLOCATED(tsurf_daily_fm_g)) DEALLOCATE(tsurf_daily_fm_g)
3603    IF (ALLOCATED(tsoil_daily_fm_g)) DEALLOCATE(tsoil_daily_fm_g)
3604    IF (ALLOCATED(soilhum_daily_fm_g))  DEALLOCATE(soilhum_daily_fm_g)
3605    IF (ALLOCATED(precip_fm_g)) DEALLOCATE(precip_fm_g)
3606    IF (ALLOCATED(gpp_daily_fm_g))  DEALLOCATE(gpp_daily_fm_g)
3607    IF (ALLOCATED(veget_fm_g)) DEALLOCATE(veget_fm_g)
3608    IF (ALLOCATED(veget_max_fm_g)) DEALLOCATE(veget_max_fm_g)
3609    IF (ALLOCATED(lai_fm_g))  DEALLOCATE(lai_fm_g)
3610   
3611    IF (ALLOCATED(isf)) DEALLOCATE(isf)
3612    IF (ALLOCATED(nf_written)) DEALLOCATE(nf_written)
3613    IF (ALLOCATED(nf_cumul)) DEALLOCATE(nf_cumul)
3614    IF (ALLOCATED(nforce)) DEALLOCATE(nforce)
3615    IF (ALLOCATED(control_moist)) DEALLOCATE(control_moist)
3616    IF (ALLOCATED(control_temp)) DEALLOCATE(control_temp)
3617    IF (ALLOCATED(soilcarbon_input)) DEALLOCATE(soilcarbon_input)
3618    IF ( ALLOCATED (horip10_index)) DEALLOCATE (horip10_index)
3619    IF ( ALLOCATED (horip100_index)) DEALLOCATE (horip100_index)
3620    IF ( ALLOCATED (horip11_index)) DEALLOCATE (horip11_index)
3621    IF ( ALLOCATED (horip101_index)) DEALLOCATE (horip101_index)
3622    IF ( ALLOCATED (co2_flux)) DEALLOCATE (co2_flux)
3623    IF ( ALLOCATED (fco2_lu)) DEALLOCATE (fco2_lu)
3624    IF ( ALLOCATED (fco2_wh)) DEALLOCATE (fco2_wh)
3625    IF ( ALLOCATED (fco2_ha)) DEALLOCATE (fco2_ha)
3626    IF ( ALLOCATED (prod10)) DEALLOCATE (prod10)
3627    IF ( ALLOCATED (prod100)) DEALLOCATE (prod100)
3628    IF ( ALLOCATED (flux10)) DEALLOCATE (flux10)
3629    IF ( ALLOCATED (flux100)) DEALLOCATE (flux100)
3630    IF ( ALLOCATED (convflux)) DEALLOCATE (convflux)
3631    IF ( ALLOCATED (cflux_prod10)) DEALLOCATE (cflux_prod10)
3632    IF ( ALLOCATED (cflux_prod100)) DEALLOCATE (cflux_prod100)
3633    IF ( ALLOCATED (prod10_harvest)) DEALLOCATE (prod10_harvest)
3634    IF ( ALLOCATED (prod100_harvest)) DEALLOCATE (prod100_harvest)
3635    IF ( ALLOCATED (flux10_harvest)) DEALLOCATE (flux10_harvest)
3636    IF ( ALLOCATED (flux100_harvest)) DEALLOCATE (flux100_harvest)
3637    IF ( ALLOCATED (convflux_harvest)) DEALLOCATE (convflux_harvest)
3638    IF ( ALLOCATED (cflux_prod10_harvest)) DEALLOCATE (cflux_prod10_harvest)
3639    IF ( ALLOCATED (cflux_prod100_harvest)) DEALLOCATE (cflux_prod100_harvest)
3640    IF ( ALLOCATED (woodharvestpft)) DEALLOCATE (woodharvestpft)
3641    IF ( ALLOCATED (convfluxpft)) DEALLOCATE (convfluxpft)
3642    IF ( ALLOCATED (fDeforestToProduct)) DEALLOCATE (fDeforestToProduct)
3643    IF ( ALLOCATED (fLulccResidue)) DEALLOCATE (fLulccResidue)
3644    IF ( ALLOCATED (fHarvestToProduct)) DEALLOCATE (fHarvestToProduct)
3645    IF ( ALLOCATED (harvest_above)) DEALLOCATE (harvest_above)
3646    IF ( ALLOCATED (soilcarbon_input_daily)) DEALLOCATE (soilcarbon_input_daily)
3647    IF ( ALLOCATED (control_temp_daily)) DEALLOCATE (control_temp_daily)
3648    IF ( ALLOCATED (control_moist_daily)) DEALLOCATE (control_moist_daily)
3649
3650    IF ( ALLOCATED (fpc_max)) DEALLOCATE (fpc_max)
3651
3652 !! 2. reset l_first
3653
3654    l_first_stomate=.TRUE.
3655
3656 !! 3. call to clear functions
3657
3658    CALL season_clear
3659    CALL stomatelpj_clear
3660    CALL littercalc_clear
3661    CALL vmax_clear
3662 
3663  END SUBROUTINE stomate_clear
3664
3665
3666!! ================================================================================================================================
3667!! SUBROUTINE   : stomate_var_init
3668!!
3669!>\BRIEF        Initialize variables of stomate with a none-zero initial value.
3670!! Subroutine is called only if ::ok_stomate = .TRUE. STOMATE diagnoses some
3671!! variables for SECHIBA : assim_param, deadleaf_cover, etc. These variables can
3672!! be recalculated from STOMATE's prognostic variables. Note that height is
3673!! saved in SECHIBA.
3674!!
3675!! DESCRIPTION  : None
3676!!
3677!! RECENT CHANGE(S) : None
3678!!
3679!! MAIN OUTPUT VARIABLE(S): leaf age (::leaf_age) and fraction of leaves in leaf
3680!! age class (::leaf_frac). The maximum water on vegetation available for
3681!! interception, fraction of soil covered by dead leaves
3682!! (::deadleaf_cover) and assimilation parameters (:: assim_param).
3683!!
3684!! REFERENCE(S) : None
3685!!
3686!! FLOWCHART    : None
3687!! \n
3688!_ ================================================================================================================================
3689 
3690  SUBROUTINE stomate_var_init &
3691       &  (kjpindex, veget_cov_max, leaf_age, leaf_frac, &
3692       &   dead_leaves, &
3693       &   veget, lai, deadleaf_cover, assim_param)
3694
3695
3696  !! 0. Variable and parameter declaration
3697
3698    !! 0.1 Input variables
3699
3700    INTEGER(i_std),INTENT(in)                             :: kjpindex        !! Domain size - terrestrial pixels only
3701    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)        :: veget           !! Fraction of pixel covered by PFT. Fraction
3702                                                                             !! accounts for none-biological land covers
3703                                                                             !! (unitless)
3704    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)        :: veget_cov_max   !! Fractional coverage: maximum share of the pixel
3705                                                                             !! covered by a PFT (unitless)
3706    REAL(r_std),DIMENSION(kjpindex,nvm,nlitt),INTENT(in)  :: dead_leaves     !! Metabolic and structural fraction of dead leaves
3707                                                                             !! per ground area
3708                                                                             !! @tex $(gC m^{-2})$ @endtex
3709    REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)        :: lai             !! Leaf area index
3710                                                                             !! @tex $(m^2 m{-2})$ @endtex
3711    REAL(r_std),DIMENSION(kjpindex,nvm,nleafages),INTENT(in) :: leaf_age     !! Age of different leaf classes per PFT (days)
3712    REAL(r_std),DIMENSION(kjpindex,nvm,nleafages),INTENT(in) :: leaf_frac    !! Fraction of leaves in leaf age class per PFT
3713                                                                             !! (unitless; 1)     
3714
3715    !! 0.2 Modified variables
3716    REAL(r_std),DIMENSION(kjpindex,nvm,npco2),INTENT(inout) :: assim_param   !! min+max+opt temperatures (K) & vmax for
3717                                                                             !! photosynthesis 
3718   
3719    !! 0.3 Output variables
3720
3721    REAL(r_std),DIMENSION(kjpindex), INTENT (out)         :: deadleaf_cover  !! Fraction of soil covered by dead leaves
3722                                                                             !! (unitless)
3723
3724
3725    ! 0.4 Local variables
3726   
3727    REAL(r_std),PARAMETER                                 :: dt_0 = zero     !! Dummy time step, must be zero
3728    REAL(r_std),DIMENSION(kjpindex,nvm)                   :: vcmax           !! Dummy vcmax
3729                                                                             !! @tex $(\mu mol m^{-2} s^{-1})$ @endtex
3730    REAL(r_std),DIMENSION(kjpindex,nvm,nleafages)         :: leaf_age_tmp    !! Temporary variable
3731    REAL(r_std),DIMENSION(kjpindex,nvm,nleafages)         :: leaf_frac_tmp   !! Temporary variable
3732                                                                             !! (unitless; 1)     
3733    INTEGER(i_std)                                        :: j               !! Index (untiless)
3734   
3735!_ ================================================================================================================================   
3736
3737
3738    ! Calculate assim_param if it was not found in the restart file
3739    IF (ALL(assim_param(:,:,:)==val_exp)) THEN
3740       ! Use temporary leaf_age_tmp and leaf_frac_tmp to preserve the input variables from being modified by the subroutine vmax.
3741       leaf_age_tmp(:,:,:)=leaf_age(:,:,:)
3742       leaf_frac_tmp(:,:,:)=leaf_frac(:,:,:)
3743
3744       !! 1.1 Calculate a temporary vcmax (stomate_vmax.f90)
3745       CALL vmax (kjpindex, dt_0, leaf_age_tmp, leaf_frac_tmp, vcmax )
3746
3747       !! 1.2 transform into nvm vegetation types
3748       assim_param(:,:,ivcmax) = zero
3749       DO j = 2, nvm
3750          assim_param(:,j,ivcmax)=vcmax(:,j)
3751       ENDDO
3752    END IF
3753   
3754    !! 2. Dead leaf cover (stomate_litter.f90)
3755    CALL deadleaf (kjpindex, veget_cov_max, dead_leaves, deadleaf_cover)     
3756   
3757  END SUBROUTINE stomate_var_init
3758
3759
3760!! ================================================================================================================================
3761!! INTERFACE    : stomate_accu
3762!!
3763!>\BRIEF        Accumulate a variable for the time period specified by
3764!! dt_sechiba or calculate the mean value over the period of dt_stomate
3765!!
3766!! DESCRIPTION : Accumulate a variable for the time period specified by
3767!! dt_sechiba or calculate the mean value over the period of dt_stomate.
3768!! stomate_accu interface can be used for variables having 1, 2 or 3 dimensions.
3769!! The corresponding subruoutine stomate_accu_r1d, stomate_accu_r2d or
3770!! stomate_accu_r3d will be selected through the interface depending on the number of dimensions.
3771!!
3772!! RECENT CHANGE(S) : None
3773!!
3774!! MAIN OUTPUT VARIABLE(S): accumulated or mean variable ::field_out::
3775!!
3776!! REFERENCE(S) : None
3777!!
3778!! FLOWCHART    : None
3779!! \n
3780!_ ================================================================================================================================
3781 
3782  SUBROUTINE stomate_accu_r1d (ldmean, field_in, field_out)
3783   
3784  !! 0. Variable and parameter declaration
3785
3786    !! 0.1 Input variables
3787    LOGICAL,INTENT(in)                     :: ldmean    !! Flag to calculate the mean over
3788    REAL(r_std),DIMENSION(:),INTENT(in)    :: field_in  !! Field that needs to be accumulated
3789   
3790    !! 0.2 Modified variables
3791    REAL(r_std),DIMENSION(:),INTENT(inout) :: field_out !! Accumulated or mean field
3792
3793!_ ================================================================================================================================
3794
3795  !! 1. Accumulate field
3796
3797    field_out(:) = field_out(:)+field_in(:)*dt_sechiba
3798   
3799  !! 2. Mean fields
3800
3801    IF (ldmean) THEN
3802       field_out(:) = field_out(:)/dt_stomate
3803    ENDIF
3804
3805  END SUBROUTINE stomate_accu_r1d
3806
3807  SUBROUTINE stomate_accu_r2d (ldmean, field_in, field_out)
3808   
3809  !! 0. Variable and parameter declaration
3810
3811    !! 0.1 Input variables
3812    LOGICAL,INTENT(in)                       :: ldmean    !! Flag to calculate the mean over
3813    REAL(r_std),DIMENSION(:,:),INTENT(in)    :: field_in  !! Field that needs to be accumulated
3814   
3815    !! 0.2 Modified variables
3816    REAL(r_std),DIMENSION(:,:),INTENT(inout) :: field_out !! Accumulated or mean field
3817
3818!_ ================================================================================================================================
3819
3820  !! 1. Accumulate field
3821
3822    field_out(:,:) = field_out(:,:)+field_in(:,:)*dt_sechiba
3823   
3824  !! 2. Mean fields
3825
3826    IF (ldmean) THEN
3827       field_out(:,:) = field_out(:,:)/dt_stomate
3828    ENDIF
3829
3830  END SUBROUTINE stomate_accu_r2d
3831
3832  SUBROUTINE stomate_accu_r3d (ldmean, field_in, field_out)
3833   
3834  !! 0. Variable and parameter declaration
3835
3836    !! 0.1 Input variables
3837    LOGICAL,INTENT(in)                         :: ldmean    !! Flag to calculate the mean over
3838    REAL(r_std),DIMENSION(:,:,:),INTENT(in)    :: field_in  !! Field that needs to be accumulated
3839   
3840    !! 0.2 Modified variables
3841    REAL(r_std),DIMENSION(:,:,:),INTENT(inout) :: field_out !! Accumulated or mean field
3842
3843!_ ================================================================================================================================
3844
3845  !! 1. Accumulate field
3846
3847    field_out(:,:,:) = field_out(:,:,:)+field_in(:,:,:)*dt_sechiba
3848   
3849  !! 2. Mean fields
3850
3851    IF (ldmean) THEN
3852       field_out(:,:,:) = field_out(:,:,:)/dt_stomate
3853    ENDIF
3854
3855  END SUBROUTINE stomate_accu_r3d
3856
3857!! ================================================================================================================================
3858!! SUBROUTINE   : init_forcing
3859!!
3860!>\BRIEF        Allocate memory for the variables containing the forcing data.
3861!! The maximum size of the allocated memory is specified in run definition file
3862!! (::max_totsize) and needs to be a compromise between charging the memory and
3863!! accessing disks to get the forcing data.
3864!!
3865!! DESCRIPTION : None
3866!!
3867!! RECENT CHANGE(S) : None
3868!!
3869!! MAIN OUTPUT VARIABLE(S): Strictly speaking the subroutine has no output
3870!! variables. However, the routine allocates memory for later use.
3871!!
3872!! REFERENCE(S) : None
3873!!
3874!! FLOWCHART    : None
3875!! \n
3876!_ ================================================================================================================================
3877 
3878  SUBROUTINE init_forcing (kjpindex,nsfm,nsft_loc)
3879   
3880  !! 0. Variable and parameter declaration
3881
3882    !! 0.1 Input variables
3883    INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
3884    INTEGER(i_std),INTENT(in) :: nsfm     !! Number of time steps that can be stored in memory (unitless)
3885    INTEGER(i_std),INTENT(in) :: nsft_loc !! Number of time steps in a year (unitless)
3886
3887   !! 0.2 Output variables
3888
3889   !! 0.3 Modified variables
3890
3891   !! 0.4 Local variables
3892
3893    LOGICAL                   :: l_error  !! Check errors in netcdf call
3894    INTEGER(i_std)            :: ier      !! Check errors in netcdf call
3895!_ ================================================================================================================================
3896   
3897  !! 1. Allocate memory
3898
3899    ! Note ::nvm is number of PFTs and ::nslm is number of soil layers
3900    l_error = .FALSE.
3901    ALLOCATE(clay_fm(kjpindex,nsfm),stat=ier)
3902    l_error = l_error .OR. (ier /= 0)
3903    IF (l_error) THEN
3904       WRITE(numout,*) 'Problem with memory allocation: forcing variables clay_fm ',kjpindex,nsfm
3905       STOP 'init_forcing'
3906    ENDIF
3907    ALLOCATE(humrel_daily_fm(kjpindex,nvm,nsfm),stat=ier)
3908    l_error = l_error .OR. (ier /= 0)
3909    IF (l_error) THEN
3910       WRITE(numout,*) 'Problem with memory allocation: forcing variables humrel_daily_fm ',kjpindex,nvm,nsfm
3911       STOP 'init_forcing'
3912    ENDIF
3913    ALLOCATE(litterhum_daily_fm(kjpindex,nsfm),stat=ier)
3914    l_error = l_error .OR. (ier /= 0)
3915    IF (l_error) THEN
3916       WRITE(numout,*) 'Problem with memory allocation: forcing variables litterhum_daily_fm ',kjpindex,nsfm
3917       STOP 'init_forcing'
3918    ENDIF
3919    ALLOCATE(t2m_daily_fm(kjpindex,nsfm),stat=ier)
3920    l_error = l_error .OR. (ier /= 0)
3921    IF (l_error) THEN
3922       WRITE(numout,*) 'Problem with memory allocation: forcing variables t2m_daily_fm ',kjpindex,nsfm
3923       STOP 'init_forcing'
3924    ENDIF
3925    ALLOCATE(t2m_min_daily_fm(kjpindex,nsfm),stat=ier)
3926    l_error = l_error .OR. (ier /= 0)
3927    IF (l_error) THEN
3928       WRITE(numout,*) 'Problem with memory allocation: forcing variables t2m_min_daily_fm ',kjpindex,nsfm
3929       STOP 'init_forcing'
3930    ENDIF
3931    ALLOCATE(tsurf_daily_fm(kjpindex,nsfm),stat=ier)
3932    l_error = l_error .OR. (ier /= 0)
3933    IF (l_error) THEN
3934       WRITE(numout,*) 'Problem with memory allocation: forcing variables tsurf_daily_fm ',kjpindex,nsfm
3935       STOP 'init_forcing'
3936    ENDIF
3937    ALLOCATE(tsoil_daily_fm(kjpindex,nslm,nsfm),stat=ier)
3938    l_error = l_error .OR. (ier /= 0)
3939    IF (l_error) THEN
3940       WRITE(numout,*) 'Problem with memory allocation: forcing variables tsoil_daily_fm ',kjpindex,nslm,nsfm
3941       STOP 'init_forcing'
3942    ENDIF
3943    ALLOCATE(soilhum_daily_fm(kjpindex,nslm,nsfm),stat=ier)
3944    l_error = l_error .OR. (ier /= 0)
3945    IF (l_error) THEN
3946       WRITE(numout,*) 'Problem with memory allocation: forcing variables soilhum_daily_fm ',kjpindex,nslm,nsfm
3947       STOP 'init_forcing'
3948    ENDIF
3949    ALLOCATE(precip_fm(kjpindex,nsfm),stat=ier)
3950    l_error = l_error .OR. (ier /= 0)
3951    IF (l_error) THEN
3952       WRITE(numout,*) 'Problem with memory allocation: forcing variables precip_fm ',kjpindex,nsfm
3953       STOP 'init_forcing'
3954    ENDIF
3955    ALLOCATE(gpp_daily_fm(kjpindex,nvm,nsfm),stat=ier)
3956    l_error = l_error .OR. (ier /= 0)
3957    IF (l_error) THEN
3958       WRITE(numout,*) 'Problem with memory allocation: forcing variables gpp_daily_fm ',kjpindex,nvm,nsfm
3959       STOP 'init_forcing'
3960    ENDIF
3961    ALLOCATE(veget_fm(kjpindex,nvm,nsfm),stat=ier)
3962    l_error = l_error .OR. (ier /= 0)
3963    IF (l_error) THEN
3964       WRITE(numout,*) 'Problem with memory allocation: forcing variables veget_fm ',kjpindex,nvm,nsfm
3965       STOP 'init_forcing'
3966    ENDIF
3967    ALLOCATE(veget_max_fm(kjpindex,nvm,nsfm),stat=ier)
3968    l_error = l_error .OR. (ier /= 0)
3969    IF (l_error) THEN
3970       WRITE(numout,*) 'Problem with memory allocation: forcing variables veget_max_fm ',kjpindex,nvm,nsfm
3971       STOP 'init_forcing'
3972    ENDIF
3973    ALLOCATE(lai_fm(kjpindex,nvm,nsfm),stat=ier)
3974    l_error = l_error .OR. (ier /= 0)
3975    IF (l_error) THEN
3976       WRITE(numout,*) 'Problem with memory allocation: forcing variables lai_fm ',kjpindex,nvm,nsfm
3977       STOP 'init_forcing'
3978    ENDIF
3979    ALLOCATE(isf(nsfm),stat=ier)
3980    l_error = l_error .OR. (ier /= 0)
3981    IF (l_error) THEN
3982       WRITE(numout,*) 'Problem with memory allocation: forcing variables isf ',nsfm
3983       STOP 'init_forcing'
3984    ENDIF
3985    ALLOCATE(nf_written(nsft_loc),stat=ier)
3986    l_error = l_error .OR. (ier /= 0)
3987    IF (l_error) THEN
3988       WRITE(numout,*) 'Problem with memory allocation: forcing variables nf_written ',nsft_loc
3989       STOP 'init_forcing'
3990    ENDIF
3991    ALLOCATE(nf_cumul(nsft_loc),stat=ier)
3992    l_error = l_error .OR. (ier /= 0)
3993    IF (l_error) THEN
3994       WRITE(numout,*) 'Problem with memory allocation: forcing variables nf_cumul ',nsft_loc
3995       STOP 'init_forcing'
3996    ENDIF
3997   
3998  !! 2. Allocate memory for the root processor only (parallel computing)
3999
4000    ! Where, ::nbp_glo is the number of global continental points
4001    IF (is_root_prc) THEN
4002       ALLOCATE(clay_fm_g(nbp_glo,nsfm),stat=ier)
4003       l_error = l_error .OR. (ier /= 0)
4004       IF (l_error) THEN
4005          WRITE(numout,*) 'Problem with memory allocation: forcing variables clay_fm_g ',nbp_glo,nsfm
4006          STOP 'init_forcing'
4007       ENDIF
4008       ALLOCATE(humrel_daily_fm_g(nbp_glo,nvm,nsfm),stat=ier)
4009       l_error = l_error .OR. (ier /= 0)
4010       IF (l_error) THEN
4011          WRITE(numout,*) 'Problem with memory allocation: forcing variables humrel_daily_fm_g ',nbp_glo,nvm,nsfm
4012          STOP 'init_forcing'
4013       ENDIF
4014       ALLOCATE(litterhum_daily_fm_g(nbp_glo,nsfm),stat=ier)
4015       l_error = l_error .OR. (ier /= 0)
4016       IF (l_error) THEN
4017          WRITE(numout,*) 'Problem with memory allocation: forcing variables litterhum_daily_fm_g ',nbp_glo,nsfm
4018          STOP 'init_forcing'
4019       ENDIF
4020       ALLOCATE(t2m_daily_fm_g(nbp_glo,nsfm),stat=ier)
4021       l_error = l_error .OR. (ier /= 0)
4022       IF (l_error) THEN
4023          WRITE(numout,*) 'Problem with memory allocation: forcing variables t2m_daily_fm_g ',nbp_glo,nsfm
4024          STOP 'init_forcing'
4025       ENDIF
4026       ALLOCATE(t2m_min_daily_fm_g(nbp_glo,nsfm),stat=ier)
4027       l_error = l_error .OR. (ier /= 0)
4028       IF (l_error) THEN
4029          WRITE(numout,*) 'Problem with memory allocation: forcing variables t2m_min_daily_fm_g ',nbp_glo,nsfm
4030          STOP 'init_forcing'
4031       ENDIF
4032       ALLOCATE(tsurf_daily_fm_g(nbp_glo,nsfm),stat=ier)
4033       l_error = l_error .OR. (ier /= 0)
4034       IF (l_error) THEN
4035          WRITE(numout,*) 'Problem with memory allocation: forcing variables tsurf_daily_fm_g ',nbp_glo,nsfm
4036          STOP 'init_forcing'
4037       ENDIF
4038       ALLOCATE(tsoil_daily_fm_g(nbp_glo,nslm,nsfm),stat=ier)
4039       l_error = l_error .OR. (ier /= 0)
4040       IF (l_error) THEN
4041          WRITE(numout,*) 'Problem with memory allocation: forcing variables tsoil_daily_fm_g ',nbp_glo,nslm,nsfm
4042          STOP 'init_forcing'
4043       ENDIF
4044       ALLOCATE(soilhum_daily_fm_g(nbp_glo,nslm,nsfm),stat=ier)
4045       l_error = l_error .OR. (ier /= 0)
4046       IF (l_error) THEN
4047          WRITE(numout,*) 'Problem with memory allocation: forcing variables soilhum_daily_fm_g ',nbp_glo,nslm,nsfm
4048          STOP 'init_forcing'
4049       ENDIF
4050       ALLOCATE(precip_fm_g(nbp_glo,nsfm),stat=ier)
4051       l_error = l_error .OR. (ier /= 0)
4052       IF (l_error) THEN
4053          WRITE(numout,*) 'Problem with memory allocation: forcing variables precip_fm_g ',nbp_glo,nsfm
4054          STOP 'init_forcing'
4055       ENDIF
4056       ALLOCATE(gpp_daily_fm_g(nbp_glo,nvm,nsfm),stat=ier)
4057       l_error = l_error .OR. (ier /= 0)
4058       IF (l_error) THEN
4059          WRITE(numout,*) 'Problem with memory allocation: forcing variables gpp_daily_fm_g ',nbp_glo,nvm,nsfm
4060          STOP 'init_forcing'
4061       ENDIF
4062       ALLOCATE(veget_fm_g(nbp_glo,nvm,nsfm),stat=ier)
4063       l_error = l_error .OR. (ier /= 0)