Context Navigation

stomate.f90 @ 6368

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

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

Context Navigation

source: tags/ORCHIDEE_2_0/ORCHIDEE/src_stomate/stomate.f90 @ 6368

Download in other formats: