Context Navigation

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

Last change on this file since 6367 was 6367, checked in by josefine.ghattas, 5 years ago
Integrated correction of diagnostic variables rhSoil and rhLitter as done in the trunk rev [6362]. See ticket #630
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 240.7 KB

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

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

Download in other formats: