Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-182/src_stomate/stomate.f90 @ 7442

Last change on this file since 7442 was 1303, checked in by josefine.ghattas, 11 years ago
Change to have unlimited time axis for stomate_forcing.nc and stomate_Cforcing.nc files. This simplifies for concatenation afterwards.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 224.6 KB

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

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

Download in other formats: