Context Navigation

stomate.f90 @ 6319

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

Improvment for the ESM CO2 configuration:

Separate variable fco2_lu into 3 parts: fco2_lu, fco2_wh and fco2_ha
Move calculation of co2_flux from dt_sechiba time-step to daily time-step (in the part for do_slow)
Removed co2_flux and fco2_lu from stomate_intialize argument list. These variables were never used in the intialization phase.
Add co2_flux, and fco2_wh, fco2_ha to restart file
Corrected output unit for nee to be consistent with LMDZ and stomate output variables. It is now in kgC/m2/s.
Corrected output for znetco2
Added fCO2_fWoodharvest and fCO2_fHarvest as new possible tracers in LMDZ (intersurf).
Added diagnostic output for fCO2_fWoodharvest and fCO2_fHarvest

Cadule

Property svn:keywords set to HeadURL Date Author Revision

File size: 241.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: