source: branches/ORCHIDEE_3_CMIP6/ORCHIDEE/src_stomate/stomate.f90 @ 8367

Last change on this file since 8367 was 7657, checked in by nicolas.vuichard, 2 years ago

add one missing OMP directive

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