Context Navigation

source: tags/ORCHIDEE_1_9_5_2/ORCHIDEE/src_stomate/stomate.f90 @ 2061

Last change on this file since 2061 was 371, checked in by martial.mancip, 13 years ago
Thank's to Didier, forget to delete a deallocation of non used resp_maint_part_fm (not used in teststomate).
File size: 131.2 KB

Line
1	!$Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/stomate.f90,v 1.41 2010/05/27 08:30:35 ssipsl Exp $
2	!IPSL (2006)
3	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
4	!
5	MODULE stomate
6	!---------------------------------------------------------------------
7	! Daily update of leaf area index etc.
8	!---------------------------------------------------------------------
9	USE netcdf
10	!-
11	USE ioipsl
12	USE defprec
13	USE grid
14	USE constantes
15	USE constantes_veg
16	USE constantes_co2
17	USE stomate_constants
18	USE stomate_io
19	USE stomate_data
20	USE stomate_season
21	USE stomate_lpj
22	USE stomate_assimtemp
23	USE stomate_litter
24	USE stomate_vmax
25	USE stomate_soilcarbon
26	USE stomate_resp
27	USE parallel
28	! USE Write_field_p
29	!-
30	IMPLICIT NONE
31	PRIVATE
32	PUBLIC stomate_main,stomate_clear,init_forcing,forcing_read
33	!
34	INTEGER,PARAMETER :: r_typ =nf90_real4
35	!-
36	! variables used inside stomate module : declaration and initialisation
37	!-
38	! total natural space (fraction of total space)
39	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: veget_cov_max
40	! new year total natural space (fraction of total space)
41	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: veget_cov_max_new
42	! carbon pool: active, slow, or passive, (gC/(m**2 of ground))
43	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: carbon
44	! density of individuals (1/(m**2 of ground))
45	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: ind
46	! daily moisture availability
47	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: humrel_daily
48	! daily litter humidity
49	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: litterhum_daily
50	! daily 2 meter temperatures (K)
51	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: t2m_daily
52	! daily minimum 2 meter temperatures (K)
53	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: t2m_min_daily
54	! daily surface temperatures (K)
55	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: tsurf_daily
56	! daily soil temperatures (K)
57	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: tsoil_daily
58	! daily soil humidity
59	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: soilhum_daily
60	! daily precipitations (mm/day) (for phenology)
61	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: precip_daily
62	! daily gross primary productivity (gC/(m**2 of ground)/day)
63	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: gpp_daily
64	! daily net primary productivity (gC/(m**2 of ground)/day)
65	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: npp_daily
66	! Turnover rates (gC/(m**2 of ground)/day)
67	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: turnover_daily
68	! Turnover rates (gC/(m*2 of ground)dtradia)
69	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: turnover_littercalc
70	! Probability of fire
71	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: fireindex
72	! Longer term total litter above the ground (gC/m**2 of ground)
73	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: firelitter
74	! "monthly" moisture availability
75	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: humrel_month
76	! "weekly" moisture availability
77	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: humrel_week
78	! "long term" 2 meter temperatures (K)
79	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: t2m_longterm
80	! "long term" 2 meter reference temperatures (K)
81	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: tlong_ref
82	! "monthly" 2 meter temperatures (K)
83	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: t2m_month
84	! "weekly" 2 meter temperatures (K)
85	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: t2m_week
86	! "monthly" soil temperatures (K)
87	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: tsoil_month
88	! "monthly" soil humidity
89	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: soilhum_month
90	! growing degree days, threshold -5 deg C (for phenology)
91	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: gdd_m5_dormance
92	! growing degree days, since midwinter (for phenology)
93	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: gdd_midwinter
94	! number of chilling days since leaves were lost (for phenology)
95	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: ncd_dormance
96	! number of growing days, threshold -5 deg C (for phenology)
97	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: ngd_minus5
98	! last year's maximum moisture availability
99	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: maxhumrel_lastyear
100	! this year's maximum moisture availability
101	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: maxhumrel_thisyear
102	! last year's minimum moisture availability
103	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: minhumrel_lastyear
104	! this year's minimum moisture availability
105	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: minhumrel_thisyear
106	! last year's maximum "weekly" GPP (gC/m**2 covered/day)
107	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: maxgppweek_lastyear
108	! this year's maximum "weekly" GPP (gC/m**2 covered/day)
109	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: maxgppweek_thisyear
110	! last year's annual GDD0
111	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: gdd0_lastyear
112	! this year's annual GDD0
113	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: gdd0_thisyear
114	! last year's annual precipitation (mm/year)
115	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: precip_lastyear
116	! this year's annual precipitation (mm/year)
117	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: precip_thisyear
118	! PFT exists (equivalent to veget > 0 for natural PFTs)
119	LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: PFTpresent
120	! "long term" net primary productivity (gC/(m**2 of ground)/year)
121	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: npp_longterm
122	! last year's maximum leaf mass, for each PFT (gC/(m**2 of ground))
123	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: lm_lastyearmax
124	! this year's maximum leaf mass, for each PFT (gC/(m**2 of ground))
125	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: lm_thisyearmax
126	! last year's maximum fpc for each natural PFT, on ground
127	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: maxfpc_lastyear
128	! this year's maximum fpc for each PFT, on ground (see stomate_season)
129	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: maxfpc_thisyear
130	! "long term" turnover rate (gC/(m**2 of ground)/year)
131	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: turnover_longterm
132	! "weekly" GPP (gC/day/(m**2 covered)
133	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: gpp_week
134	! biomass (gC/(m**2 of ground))
135	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: biomass
136	! is the plant senescent?
137	! (only for deciduous trees - carbohydrate reserve)
138	LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: senescence
139	! how many days ago was the beginning of the growing season
140	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: when_growthinit
141	! age (years). For trees, mean stand age.
142	! For grasses, ears since introduction of PFT
143	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: age
144	! Winter too cold? between 0 and 1
145	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: adapted
146	! Winter sufficiently cold? between 0 and 1
147	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: regenerate
148	! fraction of litter above the ground belonging to different PFTs
149	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: litterpart
150	! metabolic and structural litter, above and below ground
151	! (gC/(m**2 of ground))
152	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:):: litter
153	! dead leaves on ground, per PFT, metabolic and structural,
154	! in gC/(m**2 of ground)
155	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: dead_leaves
156	! black carbon on the ground (gC/(m**2 of total ground))
157	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: black_carbon
158	! ratio Lignine/Carbon in structural litter, above and below ground,
159	! (gC/(m**2 of ground))
160	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: lignin_struc
161	! carbon emitted into the atmosphere by fire (living and dead biomass)
162	! (in gC/m**2 of average ground/day)
163	!NV devient 2D
164	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: co2_fire
165	! co2 taken up (gC/(m**2 of total ground)/day)
166	!NV devient 2D
167	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: co2_to_bm_dgvm
168	! heterotrophic respiration (gC/day/m**2 of total ground)
169	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: resp_hetero_d
170	! maintenance respiration (gC/day/(m**2 of total ground))
171	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: resp_hetero_radia
172	! maintenance respiration (gC/day/(m**2 of total ground))
173	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: resp_maint_d
174	! growth respiration (gC/day/(m**2 of total ground))
175	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: resp_growth_d
176	! vegetation fractions (on ground) after last light competition
177	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: veget_lastlight
178	! is the PFT everywhere in the grid box or very localized (after its intoduction)
179	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: everywhere
180	! in order for this PFT to be introduced,
181	! does it have to be present in an adjacent grid box?
182	LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: need_adjacent
183	! leaf age (d)
184	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: leaf_age
185	! fraction of leaves in leaf age class
186	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: leaf_frac
187	! How much time ago was the PFT eliminated for the last time (y)
188	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: RIP_time
189	! duration of dormance (d)
190	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: time_lowgpp
191	! time elapsed since strongest moisture availability (d)
192	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: time_hum_min
193	! minimum moisture during dormance
194	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: hum_min_dormance
195	! time constant of probability of a leaf to be eaten by a herbivore (days)
196	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: herbivores
197	! npp total written for forcesoil...
198	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: npp_equil
199	! npp total ...
200	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: npp_tot
201	! moisture control of heterotrophic respiration
202	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: control_moist
203	! temperature control of heterotrophic respiration, above and below
204	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: control_temp
205	! quantity of carbon going into carbon pools from litter decomposition
206	! (gC/(m**2 of ground)/day)
207	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:):: soilcarbon_input
208	! quantity of carbon going into carbon pools from litter decomposition daily
209	! (gC/(m**2 of ground)/day)
210	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:):: soilcarbon_input_daily
211	! moisture control of heterotrophic respiration daily
212	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: control_moist_daily
213	! temperature control of heterotrophic respiration, above and below daily
214	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: control_temp_daily
215	! how many states were calculated for a given soil forcing time step
216	! turnover time of leaves
217	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: turnover_time
218	! daily total CO2 flux (gC/m**2/day)
219	!NV champs 2D
220	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: co2_flux_daily
221	! monthly total CO2 flux (gC/m**2)
222	!NV champs 2D
223	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: co2_flux_monthly
224	! conversion of biomass to litter (gC/(m**2 of ground))/day
225	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: bm_to_litter
226	! conversion of biomass to litter (gC/(m*2 of ground))dtradia
227	REAL(r_std), ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: bm_to_littercalc
228
229	INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: nforce
230
231	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: harvest_above_monthly, cflux_prod_monthly
232
233	! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground
234	REAL(r_std), ALLOCATABLE,SAVE,DIMENSION(:,:) :: fpc_max
235
236	! Date and EndOfYear, intialize and update in slowproc
237	! (Now managed in slowproc for land_use)
238	! time step of STOMATE in days
239	REAL(r_std),SAVE :: dt_days=zero ! Time step in days for stomate
240	! to check
241	REAL(r_std),SAVE :: day_counter=zero ! count each sechiba (dtradia) time step each day
242	! date (d)
243	INTEGER(i_std),SAVE :: date=0
244	! is it time for Stomate or update of LAI etc. ?
245	LOGICAL, SAVE :: do_slow=.FALSE.
246	! Do update of yearly variables ?
247	! This variable must be .TRUE. once a year
248	LOGICAL, SAVE :: EndOfYear=.FALSE.
249	! Land cover change flag
250	LOGICAL,SAVE :: lcchange=.FALSE.
251	! Do update of monthly variables ?
252	! This variable must be .TRUE. once a month
253	LOGICAL, SAVE :: EndOfMonth=.FALSE.
254	PUBLIC dt_days, day_counter, date, do_slow, EndOfYear, lcchange
255
256
257	! forcing data in memory
258	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: clay_fm
259	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: humrel_daily_fm
260	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: litterhum_daily_fm
261	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: t2m_daily_fm
262	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: t2m_min_daily_fm
263	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: tsurf_daily_fm
264	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: tsoil_daily_fm
265	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: soilhum_daily_fm
266	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: precip_fm
267	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: gpp_daily_fm
268	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: veget_fm
269	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: veget_max_fm
270	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: lai_fm
271	PUBLIC clay_fm, humrel_daily_fm, litterhum_daily_fm, t2m_daily_fm, t2m_min_daily_fm, tsurf_daily_fm, tsoil_daily_fm, &
272	soilhum_daily_fm, precip_fm, gpp_daily_fm, veget_fm, veget_max_fm, lai_fm
273
274	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: clay_fm_g
275	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: humrel_daily_fm_g
276	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: litterhum_daily_fm_g
277	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: t2m_daily_fm_g
278	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: t2m_min_daily_fm_g
279	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: tsurf_daily_fm_g
280	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: tsoil_daily_fm_g
281	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: soilhum_daily_fm_g
282	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: precip_fm_g
283	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: gpp_daily_fm_g
284	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: veget_fm_g
285	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: veget_max_fm_g
286	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: lai_fm_g
287
288	INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: isf
289	LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:) :: nf_written
290	INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: nf_cumul
291	PUBLIC isf, nf_written
292
293	! first call
294	LOGICAL,SAVE :: l_first_stomate = .TRUE.
295	! flag for cumul of forcing if teststomate
296	LOGICAL,SAVE :: cumul_forcing=.FALSE.
297	! flag for cumul of forcing if forcesoil
298	LOGICAL,SAVE :: cumul_Cforcing=.FALSE.
299	!
300	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: resp_maint_part_radia
301	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: resp_maint_part
302	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: resp_maint_radia
303
304	! Land cover change variables
305	! products remaining in the 10/100 year-turnover pool after the annual release for each compartment
306	! (10 or 100 + 1 : input from year of land cover change)
307	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: prod10
308	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: prod100
309	! annual release from the 10/100 year-turnover pool compartments
310	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: flux10
311	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: flux100
312	! release during first year following land cover change
313	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: convflux
314	! total annual release from the 10/100 year-turnover pool
315	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: cflux_prod10, cflux_prod100
316
317	! harvest above ground biomass for agriculture
318	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: harvest_above
319
320	! Carbon Mass total
321	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: carb_mass_total
322
323	CONTAINS
324	!
325	!=
326	!
327	SUBROUTINE stomate_main &
328	& (kjit, kjpij, kjpindex, dtradia, dt_slow, &
329	& ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, &
330	& index, lalo, neighbours, resolution, contfrac, totfrac_nobio, clay, &
331	& t2m, t2m_min, temp_sol, stempdiag, &
332	& humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, &
333	& gpp, deadleaf_cover, assim_param, &
334	& lai, height, veget, veget_max, qsintmax, &
335	& veget_max_new, totfrac_nobio_new, &
336	& hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
337	& co2_flux, fco2_lu, resp_maint,resp_hetero,resp_growth)
338	!---------------------------------------------------------------------
339	!
340	IMPLICIT NONE
341	! 0 interface description
342	!
343	! 0.1 input
344	!
345	! 0.1.1 input scalar
346	!
347	! Time step number
348	INTEGER(i_std),INTENT(in) :: kjit
349	! Domain size
350	INTEGER(i_std),INTENT(in) :: kjpindex
351	! Total size of the un-compressed grid
352	INTEGER(i_std),INTENT(in) :: kjpij
353	! Time step of SECHIBA
354	REAL(r_std),INTENT(in) :: dtradia
355	! Time step of STOMATE
356	REAL(r_std),INTENT(in) :: dt_slow
357	! Logical for _restart_ file to read
358	LOGICAL,INTENT(in) :: ldrestart_read
359	! Logical for _restart_ file to write
360	LOGICAL,INTENT(in) :: ldrestart_write
361	! Logical for _forcing_ file to write
362	LOGICAL,INTENT(in) :: ldforcing_write
363	! Logical for _carbon_forcing_ file to write
364	LOGICAL,INTENT(in) :: ldcarbon_write
365	! SECHIBA's _history_ file identifier
366	INTEGER(i_std),INTENT(in) :: hist_id
367	! SECHIBA's _history_ file 2 identifier
368	INTEGER(i_std),INTENT(in) :: hist2_id
369	! STOMATE's _Restart_ file file identifier
370	INTEGER(i_std),INTENT(in) :: rest_id_stom
371	! STOMATE's _history_ file file identifier
372	INTEGER(i_std),INTENT(in) :: hist_id_stom
373	! STOMATE's IPCC _history_ file file identifier
374	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC
375	!
376	! 0.1.2 input fields
377	!
378	! Indeces of the points on the map
379	INTEGER(i_std),DIMENSION(kjpindex),INTENT(in) :: index
380	! Geogr. coordinates (latitude,longitude) (degrees)
381	REAL(r_std),DIMENSION(kjpindex,2),INTENT(in) :: lalo
382	! neighoring grid points if land
383	INTEGER(i_std),DIMENSION(kjpindex,8),INTENT(in) :: neighbours
384	! size in x an y of the grid (m)
385	REAL(r_std),DIMENSION(kjpindex,2),INTENT(in) :: resolution
386	! fraction of continent in the grid
387	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac
388	! fraction of grid cell covered by lakes, land ice, cities, ...
389	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: totfrac_nobio
390	! clay fraction
391	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: clay
392	! Relative humidity ("moisture availability")
393	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in) :: humrel
394	! 2 m air temperature (K)
395	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: t2m
396	! min. 2 m air temp. during forcing time step (K)
397	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: t2m_min
398	! Surface temperature
399	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_sol
400	! Soil temperature
401	REAL(r_std),DIMENSION(kjpindex,nbdl),INTENT(in) :: stempdiag
402	! Relative soil moisture
403	REAL(r_std),DIMENSION(kjpindex,nbdl),INTENT(in) :: shumdiag
404	! Litter humidity
405	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: litterhumdiag
406	! Rain precipitation
407	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: precip_rain
408	! Snow precipitation
409	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: precip_snow
410	! GPP (gC/(m**2 of total ground)/time step)
411	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in) :: gpp
412	! new "maximal" coverage fraction of a PFT (LAI -> infinity) on ground
413	! only if EndOfYear is activated
414	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in) :: veget_max_new
415	! new fraction of grid cell covered by lakes, land ice, cities, ...
416	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: totfrac_nobio_new
417	!
418	! 0.2 output
419	!
420	! 0.2.1 output scalar
421	!
422	! 0.2.2 output fields
423	!
424	! CO2 flux in gC/m**2 of average ground/dt
425	!NV champs 2D
426	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: co2_flux
427	REAL(r_std),DIMENSION(kjpindex),INTENT(out) :: fco2_lu
428	! autotrophic respiration in gC/m**2 of surface/dt
429	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: resp_maint
430	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: resp_growth
431	!NV champs2D
432	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: resp_hetero
433	!
434	! 0.3 modified
435	!
436	! 0.3.1 modified scalar
437	! 0.3.2 modified fields
438	!
439	! Surface foliere
440	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(inout) :: lai
441	! Fraction of vegetation type from hydrological module. Takes into account ice etc.
442	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(inout) :: veget
443	! Maximum fraction of vegetation type from hydrological module. Takes into account ice etc.
444	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(inout) :: veget_max
445	! height (m)
446	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(inout) :: height
447	! min+max+opt temps & vmax for photosynthesis
448	REAL(r_std),DIMENSION(kjpindex,nvm,npco2),INTENT(inout) :: assim_param
449	! fraction of soil covered by dead leaves
450	REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: deadleaf_cover
451	! Maximum water on vegetation for interception
452	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(inout) :: qsintmax
453	!
454	! 0.4 local declaration
455	!
456	! 0.4.1 variables defined on nvm in SECHIBA, nvm in STOMATE/LPJ
457
458	! gross primary productivity (gC/(m**2 of total ground)/day)
459	REAL(r_std),DIMENSION(kjpindex,nvm) :: gpp_d
460	! fractional coverage: actually covered space, taking into account LAI (= grid scale fpc).
461	! Fraction of ground.
462	REAL(r_std),DIMENSION(kjpindex,nvm) :: veget_cov
463
464	! 0.4.2 other
465	!
466	! soil level used for LAI
467	INTEGER(i_std),SAVE :: lcanop
468	! counts time until next STOMATE time step
469	REAL(r_std) :: day_counter_read
470	! STOMATE time step read in restart file
471	REAL(r_std) :: dt_days_read
472	! STOMATE date
473	INTEGER(i_std) :: date_read
474	! Maximum STOMATE time step (days)
475	REAL(r_std),PARAMETER :: max_dt_days = 5.
476	! Writing frequency for history file (d)
477	REAL(r_std) :: hist_days
478	! precipitation (mm/day)
479	REAL(r_std),DIMENSION(kjpindex) :: precip
480	! Maximum rate of carboxylation
481	REAL(r_std),DIMENSION(kjpindex,nvm) :: vcmax
482	! Maximum rate of RUbp regeneration
483	REAL(r_std),DIMENSION(kjpindex,nvm) :: vjmax
484	! Min temperature for photosynthesis (deg C)
485	REAL(r_std),DIMENSION(kjpindex,nvm) :: t_photo_min
486	! Opt temperature for photosynthesis (deg C)
487	REAL(r_std),DIMENSION(kjpindex,nvm) :: t_photo_opt
488	! Max temperature for photosynthesis (deg C)
489	REAL(r_std),DIMENSION(kjpindex,nvm) :: t_photo_max
490	!NV supprimé
491	! total autotrophic respiration (gC/day/(m**2 of total ground))
492	! REAL(r_std),DIMENSION(kjpindex) :: resp_auto_tot
493	!NV gpp_tot supprimé
494	! total photosynthesis (gC/day/(m**2 of total ground))
495	! REAL(r_std),DIMENSION(kjpindex) :: gpp_tot
496	! -- LOOP
497	REAL(r_std) :: net_co2_flux_monthly, net_co2_flux_monthly_sum
498	INTEGER :: ios
499	! for forcing file: "daily" gpp
500	REAL(r_std),DIMENSION(kjpindex,nvm) :: gpp_daily_x
501	! total "vegetation" cover
502	REAL(r_std),DIMENSION(kjpindex) :: cvegtot
503	INTEGER(i_std) :: ji, jv, i, j
504	REAL(r_std) :: trans_veg
505	REAL(r_std) :: tmp_day(1)
506	!
507
508	INTEGER(i_std) :: ier
509	! moisture control of heterotrophic respiration
510	REAL(r_std),DIMENSION(kjpindex,nlevs) :: control_moist_inst
511	! temperature control of heterotrophic respiration, above and below
512	REAL(r_std),DIMENSION(kjpindex,nlevs) :: control_temp_inst
513	! quantity of carbon going into carbon pools from litter decomposition
514	! (gC/(m**2 of ground)/day)
515	REAL(r_std),DIMENSION(kjpindex,ncarb,nvm) :: soilcarbon_input_inst
516	! time step of soil forcing file (in days)
517	REAL(r_std),SAVE :: dt_forcesoil
518	INTEGER(i_std),SAVE :: nparan ! Number of time steps per year for carbon spinup
519	INTEGER(i_std),SAVE :: nbyear
520	INTEGER(i_std),PARAMETER :: nparanmax=366 ! Number max of time steps per year for carbon spinup
521	REAL(r_std) :: sf_time
522	INTEGER(i_std),SAVE :: iatt
523	INTEGER(i_std),SAVE :: iatt_old=1
524	INTEGER(i_std) :: max_totsize, totsize_1step,totsize_tmp
525	REAL(r_std) :: xn
526	INTEGER(i_std),SAVE :: nsfm, nsft
527	INTEGER(i_std),SAVE :: iisf
528	!-
529	CHARACTER(LEN=100), SAVE :: forcing_name,Cforcing_name
530	INTEGER(i_std),SAVE :: Cforcing_id
531	INTEGER(i_std),PARAMETER :: ndm = 10
532	INTEGER(i_std) :: vid
533	INTEGER(i_std) :: nneigh,direct
534	INTEGER(i_std),DIMENSION(ndm) :: d_id
535
536
537	! root temperature (convolution of root and soil temperature profiles)
538	REAL(r_std),DIMENSION(kjpindex,nvm) :: t_root
539	REAL(r_std),DIMENSION(kjpindex,nvm,nparts) :: coeff_maint
540	! temperature which is pertinent for maintenance respiration (K)
541	REAL(r_std),DIMENSION(kjpindex,nparts) :: t_maint_radia
542	! integration constant for root profile
543	REAL(r_std),DIMENSION(kjpindex) :: rpc
544	! long term annual mean temperature, C
545	REAL(r_std),DIMENSION(kjpindex) :: tl
546	! slope of maintenance respiration coefficient (1/K)
547	REAL(r_std),DIMENSION(kjpindex) :: slope
548	! soil levels (m)
549	REAL(r_std),DIMENSION(0:nbdl) :: z_soil
550	! root depth. This will, one day, be a prognostic variable.
551	! It will be calculated by
552	! STOMATE (save in restart file & give to hydrology module!),
553	! probably somewhere
554	! in the allocation routine. For the moment, it is prescribed.
555	REAL(r_std),DIMENSION(kjpindex,nvm) :: rprof
556	INTEGER(i_std) :: l,k
557	! litter heterotrophic respiration (gC/day/m**2 by PFT)
558	REAL(r_std),DIMENSION(kjpindex,nvm) :: resp_hetero_litter
559	! soil heterotrophic respiration (gC/day/m**2 by PFT)
560	REAL(r_std),DIMENSION(kjpindex,nvm) :: resp_hetero_soil
561	INTEGER(i_std) :: iyear
562
563	REAL(r_std),DIMENSION(nbp_glo) :: clay_g
564	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:,:) :: soilcarbon_input_g
565	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: control_moist_g
566	REAL(r_std),ALLOCATABLE,DIMENSION(:,:,:) :: control_temp_g
567	REAL(r_std),ALLOCATABLE,DIMENSION(:,:) :: npp_equil_g
568
569	REAL(r_std), DIMENSION(kjpindex) :: vartmp
570	REAL(r_std) :: net_cflux_prod_monthly_sum , net_cflux_prod_monthly_tot
571	REAL(r_std) :: net_harvest_above_monthly_sum, net_harvest_above_monthly_tot
572	REAL(r_std) :: net_biosp_prod_monthly_sum , net_biosp_prod_monthly_tot
573	!---------------------------------------------------------------------
574	! first of all: store time step in common value
575	itime = kjit
576
577	z_soil(0) = zero
578	z_soil(1:nbdl) = diaglev(1:nbdl)
579	DO j=1,nvm
580	rprof(:,j) = 1./humcste(j)
581	ENDDO
582	!-
583	! 1 do initialisation
584	!-
585	resp_growth=zero
586	IF (l_first_stomate) THEN
587	IF (long_print) THEN
588	WRITE (numout,*) ' l_first_stomate : call stomate_init'
589	ENDIF
590	!
591	! 1.0 Initialization of constant tabs without bare_soil
592	!
593	CALL stomate_constants_init ()
594	!
595	! 1.1 allocation, file definitions. Set flags.
596	!
597	CALL stomate_init (kjpij, kjpindex, index, ldforcing_write, lalo, &
598	rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
599
600	!
601	! 1.2 read PFT data
602	!
603	CALL data (kjpindex, lalo)
604	!
605	! 1.3 Initial conditions
606	!
607	! 1.3.1 read STOMATE's start file
608	!
609	co2_flux(:,:) = zero
610	fco2_lu(:) = zero
611	!
612	CALL readstart &
613	& (kjpindex, index, lalo, resolution, &
614	& day_counter_read, dt_days_read, date_read, &
615	& ind, adapted, regenerate, &
616	& humrel_daily, litterhum_daily, &
617	& t2m_daily, t2m_min_daily, tsurf_daily, tsoil_daily, &
618	& soilhum_daily, precip_daily, &
619	& gpp_daily, npp_daily, turnover_daily, &
620	& humrel_month, humrel_week, &
621	& t2m_longterm, tlong_ref, t2m_month, t2m_week, &
622	& tsoil_month, soilhum_month, fireindex, firelitter, &
623	& maxhumrel_lastyear, maxhumrel_thisyear, &
624	& minhumrel_lastyear, minhumrel_thisyear, &
625	& maxgppweek_lastyear, maxgppweek_thisyear, &
626	& gdd0_lastyear, gdd0_thisyear, &
627	& precip_lastyear, precip_thisyear, &
628	& gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
629	& PFTpresent, npp_longterm, lm_lastyearmax, lm_thisyearmax, &
630	& maxfpc_lastyear, maxfpc_thisyear, &
631	& turnover_longterm, gpp_week, biomass, resp_maint_part, &
632	& leaf_age, leaf_frac, &
633	& senescence, when_growthinit, age, &
634	& resp_hetero_d, resp_maint_d, resp_growth_d, co2_fire, co2_to_bm_dgvm, &
635	& veget_lastlight, everywhere, need_adjacent, RIP_time, time_lowgpp, &
636	& time_hum_min, hum_min_dormance, &
637	& litterpart, litter, dead_leaves, &
638	& carbon, black_carbon, lignin_struc,turnover_time,&
639	& prod10,prod100,flux10, flux100, &
640	& convflux, cflux_prod10, cflux_prod100, bm_to_litter, carb_mass_total)
641
642	! 1.4 read the boundary conditions
643	!
644	CALL readbc (kjpindex, lalo, resolution, tlong_ref)
645	!
646	! 1.5 check time step
647	!
648	! 1.5.1 allow STOMATE's time step to change
649	! although this is dangerous
650	!
651	IF (dt_days /= dt_days_read) THEN
652	WRITE(numout,*) 'slow_processes: STOMATE time step changes:', &
653	dt_days_read,' -> ',dt_days
654	ENDIF
655
656	! 1.5.2 time step has to be a multiple of a full day
657
658	IF ( ( dt_days-REAL(NINT(dt_days),r_std) ) > min_stomate ) THEN
659	WRITE(numout,*) 'slow_processes: STOMATE time step is wrong:', &
660	& dt_days,' days.'
661	STOP
662	ENDIF
663
664	! 1.5.3 upper limit to STOMATE's time step
665
666	IF ( dt_days > max_dt_days ) THEN
667	WRITE(numout,*) 'slow_processes: STOMATE time step too large:', &
668	& dt_days,' days.'
669	STOP
670	ENDIF
671
672	! 1.5.4 STOMATE time step must not be less than the forcing time step
673
674	IF ( dtradia > dt_days*one_day ) THEN
675	WRITE(numout,*) &
676	& 'slow_processes: STOMATE time step smaller than forcing time step.'
677	STOP
678	ENDIF
679
680	! 1.5.5 For teststomate : test day_counter
681	IF ( abs(day_counter - day_counter_read) > min_stomate ) THEN
682	WRITE(numout,*) 'slow_processes: STOMATE day counter changes:', &
683	day_counter_read,' -> ',day_counter
684	ENDIF
685
686	! 1.5.6 date check
687	IF (date /= date_read) THEN
688	WRITE(numout,*) 'slow_processes: STOMATE date changes:', &
689	date_read,' -> ',date
690	ENDIF
691
692	! 1.5.6 some more messages
693
694	WRITE(numout,*) 'slow_processes, STOMATE time step (d): ', dt_days
695
696	!
697	! 1.6 write forcing file for stomate?
698	!
699	IF (ldforcing_write) THEN
700
701	!Config Key = STOMATE_FORCING_NAME
702	!Config Desc = Name of STOMATE's forcing file
703	!Config Def = NONE
704	!Config Help = Name that will be given
705	!Config to STOMATE's offline forcing file
706	!-
707	forcing_name = stomate_forcing_name ! compatibilité avec driver Nicolas
708	CALL getin_p('STOMATE_FORCING_NAME',forcing_name)
709
710	IF ( TRIM(forcing_name) /= 'NONE' ) THEN
711	IF (is_root_prc) CALL SYSTEM ('rm -f '//TRIM(forcing_name))
712	WRITE(numout,*) 'writing a forcing file for STOMATE.'
713
714	!Config Key = STOMATE_FORCING_MEMSIZE
715	!Config Desc = Size of STOMATE forcing data in memory (MB)
716	!Config Def = 50
717	!Config Help = This variable determines how many
718	!Config forcing states will be kept in memory.
719	!Config Must be a compromise between memory
720	!Config use and frequeny of disk access.
721
722	max_totsize = 50
723	CALL getin_p('STOMATE_FORCING_MEMSIZE', max_totsize)
724	max_totsize = max_totsize*1000000
725
726	totsize_1step = &
727	& SIZE(clay)*KIND(clay) &
728	& +SIZE(humrel_daily)*KIND(humrel_daily) &
729	& +SIZE(litterhum_daily)*KIND(litterhum_daily) &
730	& +SIZE(t2m_daily)*KIND(t2m_daily) &
731	& +SIZE(t2m_min_daily)*KIND(t2m_min_daily) &
732	& +SIZE(tsurf_daily)*KIND(tsurf_daily) &
733	& +SIZE(tsoil_daily)*KIND(tsoil_daily) &
734	& +SIZE(soilhum_daily)*KIND(soilhum_daily) &
735	& +SIZE(precip_daily)*KIND(precip_daily) &
736	& +SIZE(gpp_daily_x)*KIND(gpp_daily_x) &
737	& +SIZE(veget)*KIND(veget) &
738	& +SIZE(veget_max)*KIND(veget_max) &
739	& +SIZE(lai)*KIND(lai)
740
741	CALL reduce_sum(totsize_1step,totsize_tmp)
742	CALL bcast(totsize_tmp)
743	totsize_1step=totsize_tmp
744	! total number of forcing steps
745	nsft = INT(one_year/(dt_slow/one_day))
746
747	! number of forcing steps in memory
748	nsfm = MIN(nsft, &
749	& MAX(1,NINT( REAL(max_totsize,r_std) &
750	& /REAL(totsize_1step,r_std))))
751
752	CALL init_forcing (kjpindex,nsfm,nsft)
753
754	isf(:) = (/ (i,i=1,nsfm) /)
755	nf_written(:) = .FALSE.
756	nf_cumul(:) = 0
757
758	iisf = 0
759
760	!-
761	IF (is_root_prc) THEN
762	ier = NF90_CREATE (TRIM(forcing_name),NF90_SHARE,forcing_id)
763	ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL,'dtradia',dtradia)
764	ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL,'dt_slow',dt_slow)
765	ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL, &
766	& 'nsft',REAL(nsft,r_std))
767	ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL, &
768	& 'kjpij',REAL(iim_g*jjm_g,r_std))
769	ier = NF90_PUT_ATT (forcing_id,NF90_GLOBAL, &
770	& 'kjpindex',REAL(nbp_glo,r_std))
771	!-
772	ier = NF90_DEF_DIM (forcing_id,'points',nbp_glo,d_id(1))
773	ier = NF90_DEF_DIM (forcing_id,'layers',nbdl,d_id(2))
774	ier = NF90_DEF_DIM (forcing_id,'pft',nvm,d_id(3))
775	direct=2
776	ier = NF90_DEF_DIM (forcing_id,'direction',direct,d_id(4))
777	nneigh=8
778	ier = NF90_DEF_DIM (forcing_id,'nneigh',nneigh,d_id(5))
779	ier = NF90_DEF_DIM (forcing_id,'time',nsft,d_id(6))
780	ier = NF90_DEF_DIM (forcing_id,'nbparts',nparts,d_id(7))
781	!-
782	ier = NF90_DEF_VAR (forcing_id,'points', r_typ,d_id(1),vid)
783	ier = NF90_DEF_VAR (forcing_id,'layers', r_typ,d_id(2),vid)
784	ier = NF90_DEF_VAR (forcing_id,'pft', r_typ,d_id(3),vid)
785	ier = NF90_DEF_VAR (forcing_id,'direction', r_typ,d_id(4),vid)
786	ier = NF90_DEF_VAR (forcing_id,'nneigh', r_typ,d_id(5),vid)
787	ier = NF90_DEF_VAR (forcing_id,'time', r_typ,d_id(6),vid)
788	ier = NF90_DEF_VAR (forcing_id,'nbparts', r_typ,d_id(7),vid)
789	ier = NF90_DEF_VAR (forcing_id,'index', r_typ,d_id(1),vid)
790	ier = NF90_DEF_VAR (forcing_id,'contfrac', r_typ,d_id(1),vid)
791	ier = NF90_DEF_VAR (forcing_id,'lalo', &
792	& r_typ,(/ d_id(1),d_id(4) /),vid)
793	ier = NF90_DEF_VAR (forcing_id,'neighbours', &
794	& r_typ,(/ d_id(1),d_id(5) /),vid)
795	ier = NF90_DEF_VAR (forcing_id,'resolution', &
796	& r_typ,(/ d_id(1),d_id(4) /),vid)
797	ier = NF90_DEF_VAR (forcing_id,'clay', &
798	& r_typ,(/ d_id(1),d_id(6) /),vid)
799	ier = NF90_DEF_VAR (forcing_id,'humrel', &
800	& r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
801	ier = NF90_DEF_VAR (forcing_id,'litterhum', &
802	& r_typ,(/ d_id(1),d_id(6) /),vid)
803	ier = NF90_DEF_VAR (forcing_id,'t2m', &
804	& r_typ,(/ d_id(1),d_id(6) /),vid)
805	ier = NF90_DEF_VAR (forcing_id,'t2m_min', &
806	& r_typ,(/ d_id(1),d_id(6) /),vid)
807	ier = NF90_DEF_VAR (forcing_id,'tsurf', &
808	& r_typ,(/ d_id(1),d_id(6) /),vid)
809	ier = NF90_DEF_VAR (forcing_id,'tsoil', &
810	& r_typ,(/ d_id(1),d_id(2),d_id(6) /),vid)
811	ier = NF90_DEF_VAR (forcing_id,'soilhum', &
812	& r_typ,(/ d_id(1),d_id(2),d_id(6) /),vid)
813	ier = NF90_DEF_VAR (forcing_id,'precip', &
814	& r_typ,(/ d_id(1),d_id(6) /),vid)
815	ier = NF90_DEF_VAR (forcing_id,'gpp', &
816	& r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
817	ier = NF90_DEF_VAR (forcing_id,'veget', &
818	& r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
819	ier = NF90_DEF_VAR (forcing_id,'veget_max', &
820	& r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
821	ier = NF90_DEF_VAR (forcing_id,'lai', &
822	& r_typ,(/ d_id(1),d_id(3),d_id(6) /),vid)
823	ier = NF90_ENDDEF (forcing_id)
824	!-
825	ier = NF90_INQ_VARID (forcing_id,'points',vid)
826	ier = NF90_PUT_VAR (forcing_id,vid, &
827	& (/(REAL(i,r_std),i=1,nbp_glo) /))
828	ier = NF90_INQ_VARID (forcing_id,'layers',vid)
829	ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,nbdl)/))
830	ier = NF90_INQ_VARID (forcing_id,'pft',vid)
831	ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,nvm)/))
832	ier = NF90_INQ_VARID (forcing_id,'direction',vid)
833	ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,2)/))
834	ier = NF90_INQ_VARID (forcing_id,'nneigh',vid)
835	ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,8)/))
836	ier = NF90_INQ_VARID (forcing_id,'time',vid)
837	ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,nsft)/))
838	ier = NF90_INQ_VARID (forcing_id,'nbparts',vid)
839	ier = NF90_PUT_VAR (forcing_id,vid,(/(REAL(i,r_std),i=1,nparts)/))
840	ier = NF90_INQ_VARID (forcing_id,'index',vid)
841	ier = NF90_PUT_VAR (forcing_id,vid,REAL(index_g,r_std))
842	ier = NF90_INQ_VARID (forcing_id,'contfrac',vid)
843	ier = NF90_PUT_VAR (forcing_id,vid,REAL(contfrac_g,r_std))
844	ier = NF90_INQ_VARID (forcing_id,'lalo',vid)
845	ier = NF90_PUT_VAR (forcing_id,vid,lalo_g)
846	!ym attention a neighbours, a modifier plus tard
847	ier = NF90_INQ_VARID (forcing_id,'neighbours',vid)
848	ier = NF90_PUT_VAR (forcing_id,vid,REAL(neighbours_g,r_std))
849	ier = NF90_INQ_VARID (forcing_id,'resolution',vid)
850	ier = NF90_PUT_VAR (forcing_id,vid,resolution_g)
851	ENDIF
852	ENDIF
853	ENDIF
854	!
855	! 1.7 write forcing file for the soil?
856	!
857	IF (ldcarbon_write) THEN
858	!
859	!Config Key = STOMATE_CFORCING_NAME
860	!Config Desc = Name of STOMATE's carbon forcing file
861	!Config Def = NONE
862	!Config Help = Name that will be given to STOMATE's carbon
863	!Config offline forcing file
864	!-
865	Cforcing_name = stomate_Cforcing_name ! compatibilité avec driver Nicolas
866	CALL getin_p('STOMATE_CFORCING_NAME',Cforcing_name)
867
868	IF ( TRIM(Cforcing_name) /= 'NONE' ) THEN
869	IF (is_root_prc) CALL SYSTEM ('rm -f '//TRIM(Cforcing_name))
870	!
871	! time step of forcesoil
872	!
873	!Config Key = FORCESOIL_STEP_PER_YEAR
874	!Config Desc = Number of time steps per year for carbon spinup.
875	!Config Def = 365
876	!Config Help = Number of time steps per year for carbon spinup.
877	nparan = 365
878	CALL getin_p('FORCESOIL_STEP_PER_YEAR', nparan)
879
880	IF ( nparan < 1 ) nparan = 1
881
882	!Config Key = FORCESOIL_NB_YEAR
883	!Config Desc = Number of years saved for carbon spinup.
884	!Config If internal parameter cumul_Cforcing is TRUE in stomate.f90
885	!config Then this parameter is forced to one.
886	!Config Def = 1
887	!Config Help = Number of years saved for carbon spinup.
888	nbyear=1
889	CALL getin_p('FORCESOIL_NB_YEAR', nbyear)
890	IF ( cumul_Cforcing ) THEN
891	CALL ipslerr (1,'stomate', &
892	& 'Internal parameter cumul_Cforcing is TRUE in stomate.f90', &
893	& 'Parameter FORCESOIL_NB_YEAR is forced to 1.', &
894	& 'Then variable nbyear is forced to 1.')
895	nbyear=1
896	ENDIF
897
898	dt_forcesoil = zero
899	nparan = nparan+1
900	DO WHILE (dt_forcesoil < dt_slow/one_day)
901	nparan = nparan-1
902	IF (nparan < 1) THEN
903	STOP 'Problem 1 with number of soil forcing time steps.'
904	ENDIF
905	dt_forcesoil = one_year/REAL(nparan,r_std)
906	ENDDO
907
908	IF ( nparan > nparanmax ) THEN
909	STOP 'Problem 2 with number of soil forcing time steps.'
910	ENDIF
911
912	WRITE(numout,*) 'time step of soil forcing (d): ',dt_forcesoil
913
914	ALLOCATE( nforce(nparan*nbyear))
915	nforce(:) = 0
916
917	ALLOCATE(control_moist(kjpindex,nlevs,nparan*nbyear))
918	ALLOCATE(npp_equil(kjpindex,nparan*nbyear))
919	ALLOCATE(npp_tot(kjpindex))
920	ALLOCATE(control_temp(kjpindex,nlevs,nparan*nbyear))
921	ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan*nbyear))
922
923	npp_equil(:,:) = zero
924	npp_tot(:) = zero
925	control_moist(:,:,:) = zero
926	control_temp(:,:,:) = zero
927	soilcarbon_input(:,:,:,:) = zero
928	ENDIF
929	ENDIF
930	!
931	! 1.8 calculate STOMATE's vegetation fractions
932	! from veget, veget_max
933	!
934	DO j=1,nvm
935	WHERE ((1.-totfrac_nobio(:)) > min_sechiba)
936	veget_cov(:,j)=veget(:,j)/( 1.-totfrac_nobio(:) )
937	veget_cov_max(:,j)=veget_max(:,j)/( 1.-totfrac_nobio(:) )
938	ELSEWHERE
939	veget_cov(:,j) = zero
940	veget_cov_max(:,j) = zero
941	ENDWHERE
942	ENDDO
943	IF ( EndOfYear ) THEN
944	DO j=1,nvm
945	WHERE ((1.-totfrac_nobio_new(:)) > min_sechiba)
946	veget_cov_max_new(:,j)=veget_max_new(:,j)/( 1.-totfrac_nobio_new(:) )
947	ELSEWHERE
948	veget_cov_max_new(:,j) = zero
949	ENDWHERE
950	ENDDO
951	ENDIF
952	!
953	! 1.9 initialize some variables
954	! STOMATE diagnoses some variables for SECHIBA :
955	! assim_param, deadleaf_cover, etc.
956	! These variables can be recalculated easily
957	! from STOMATE's prognostic variables.
958	! height is saved in Sechiba.
959	!
960	IF (control%ok_stomate) THEN
961	CALL stomate_var_init &
962	& (kjpindex, veget_cov, veget_cov_max, leaf_age, leaf_frac, &
963	& tlong_ref, t2m_month, dead_leaves, &
964	& veget, lai, qsintmax, deadleaf_cover, assim_param)
965	harvest_above = zero
966	ENDIF
967	!
968	! 1.10 reset flag
969	!
970	l_first_stomate = .FALSE.
971	!
972	! 1.11 return
973	!
974	RETURN
975	ENDIF ! first call
976	IF (bavard >= 4) THEN
977	WRITE(numout,*) 'DATE ',date,' ymds', year, month, day, sec, '-- stp --', itime, do_slow
978	ENDIF
979	!-
980	! 2 prepares restart file for the next simulation
981	!-
982	IF (ldrestart_write) THEN
983	IF (long_print) THEN
984	WRITE (numout,*) &
985	& ' we have to complete restart file with STOMATE variables'
986	ENDIF
987	CALL writerestart &
988	& (kjpindex, index, &
989	& day_counter, dt_days, date, &
990	& ind, adapted, regenerate, &
991	& humrel_daily, litterhum_daily, &
992	& t2m_daily, t2m_min_daily, tsurf_daily, tsoil_daily, &
993	& soilhum_daily, precip_daily, &
994	& gpp_daily, npp_daily, turnover_daily, &
995	& humrel_month, humrel_week, &
996	& t2m_longterm, tlong_ref, t2m_month, t2m_week, &
997	& tsoil_month, soilhum_month, fireindex, firelitter, &
998	& maxhumrel_lastyear, maxhumrel_thisyear, &
999	& minhumrel_lastyear, minhumrel_thisyear, &
1000	& maxgppweek_lastyear, maxgppweek_thisyear, &
1001	& gdd0_lastyear, gdd0_thisyear, &
1002	& precip_lastyear, precip_thisyear, &
1003	& gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
1004	& PFTpresent, npp_longterm, lm_lastyearmax, lm_thisyearmax, &
1005	& maxfpc_lastyear, maxfpc_thisyear, &
1006	& turnover_longterm, gpp_week, biomass, resp_maint_part, &
1007	& leaf_age, leaf_frac, &
1008	& senescence, when_growthinit, age, &
1009	& resp_hetero_d, resp_maint_d, resp_growth_d, co2_fire, co2_to_bm_dgvm, &
1010	& veget_lastlight, everywhere, need_adjacent, &
1011	& RIP_time, time_lowgpp, &
1012	& time_hum_min, hum_min_dormance, &
1013	& litterpart, litter, dead_leaves, &
1014	& carbon, black_carbon, lignin_struc,turnover_time,&
1015	& prod10,prod100,flux10, flux100, &
1016	& convflux, cflux_prod10, cflux_prod100, bm_to_litter,carb_mass_total)
1017
1018	IF (ldforcing_write .AND. TRIM(forcing_name) /= 'NONE' ) THEN
1019	CALL forcing_write(forcing_id,1,iisf)
1020	!
1021	IF (is_root_prc) ier = NF90_CLOSE (forcing_id)
1022	forcing_id=-1
1023	ENDIF
1024
1025	IF (ldcarbon_write .AND. TRIM(Cforcing_name) /= 'NONE' ) THEN
1026	WRITE(numout,*) &
1027	& 'stomate: writing the forcing file for carbon spinup'
1028	!
1029	DO iatt = 1, nparan*nbyear
1030	IF ( nforce(iatt) > 0 ) THEN
1031	soilcarbon_input(:,:,:,iatt) = &
1032	& soilcarbon_input(:,:,:,iatt)/REAL(nforce(iatt),r_std)
1033	control_moist(:,:,iatt) = &
1034	& control_moist(:,:,iatt)/REAL(nforce(iatt),r_std)
1035	control_temp(:,:,iatt) = &
1036	& control_temp(:,:,iatt)/REAL(nforce(iatt),r_std)
1037	npp_equil(:,iatt) = &
1038	& npp_equil(:,iatt)/REAL(nforce(iatt),r_std)
1039	ELSE
1040	WRITE(numout,*) &
1041	& 'We have no soil carbon forcing data for this time step:', &
1042	& iatt
1043	WRITE(numout,*) ' -> we set them to zero'
1044	soilcarbon_input(:,:,:,iatt) = zero
1045	control_moist(:,:,iatt) = zero
1046	control_temp(:,:,iatt) = zero
1047	npp_equil(:,iatt) = zero
1048	ENDIF
1049	ENDDO
1050
1051	IF (is_root_prc) THEN
1052	ALLOCATE(soilcarbon_input_g(nbp_glo,ncarb,nvm,nparan*nbyear))
1053	ALLOCATE(control_moist_g(nbp_glo,nlevs,nparan*nbyear))
1054	ALLOCATE(control_temp_g(nbp_glo,nlevs,nparan*nbyear))
1055	ALLOCATE(npp_equil_g(nbp_glo,nparan*nbyear))
1056	ENDIF
1057
1058	CALL gather(clay,clay_g)
1059	CALL gather(soilcarbon_input,soilcarbon_input_g)
1060	CALL gather(control_moist,control_moist_g)
1061	CALL gather(control_temp,control_temp_g)
1062	CALL gather(npp_equil,npp_equil_g)
1063
1064	!-
1065	IF (is_root_prc) THEN
1066	ier = NF90_CREATE (TRIM(Cforcing_name),NF90_WRITE,Cforcing_id)
1067	ier = NF90_PUT_ATT (Cforcing_id,NF90_GLOBAL, &
1068	& 'kjpindex',REAL(nbp_glo,r_std))
1069	ier = NF90_PUT_ATT (Cforcing_id,NF90_GLOBAL, &
1070	& 'nparan',REAL(nparan,r_std))
1071	ier = NF90_PUT_ATT (Cforcing_id,NF90_GLOBAL, &
1072	& 'nbyear',REAL(nbyear,r_std))
1073	ier = NF90_DEF_DIM (Cforcing_id,'points',nbp_glo,d_id(1))
1074	ier = NF90_DEF_DIM (Cforcing_id,'carbtype',ncarb,d_id(2))
1075	ier = NF90_DEF_DIM (Cforcing_id,'vegtype',nvm,d_id(3))
1076	ier = NF90_DEF_DIM (Cforcing_id,'level',nlevs,d_id(4))
1077	ier = NF90_DEF_DIM (Cforcing_id,'time_step',nparan*nbyear,d_id(5))
1078	!-
1079	ier = NF90_DEF_VAR (Cforcing_id,'points', r_typ,d_id(1),vid)
1080	ier = NF90_DEF_VAR (Cforcing_id,'carbtype', r_typ,d_id(2),vid)
1081	ier = NF90_DEF_VAR (Cforcing_id,'vegtype', r_typ,d_id(3),vid)
1082	ier = NF90_DEF_VAR (Cforcing_id,'level', r_typ,d_id(4),vid)
1083	ier = NF90_DEF_VAR (Cforcing_id,'time_step', r_typ,d_id(5),vid)
1084	ier = NF90_DEF_VAR (Cforcing_id,'index', r_typ,d_id(1),vid)
1085	ier = NF90_DEF_VAR (Cforcing_id,'clay', r_typ,d_id(1),vid)
1086	ier = NF90_DEF_VAR (Cforcing_id,'soilcarbon_input',r_typ, &
1087	& (/ d_id(1),d_id(2),d_id(3),d_id(5) /),vid)
1088	ier = NF90_DEF_VAR (Cforcing_id,'control_moist',r_typ, &
1089	& (/ d_id(1),d_id(4),d_id(5) /),vid)
1090	ier = NF90_DEF_VAR (Cforcing_id,'control_temp',r_typ, &
1091	& (/ d_id(1),d_id(4),d_id(5) /),vid)
1092	ier = NF90_DEF_VAR (Cforcing_id,'npp_equil',r_typ, &
1093	& (/ d_id(1),d_id(5) /),vid)
1094	ier = NF90_ENDDEF (Cforcing_id)
1095	!-
1096	ier = NF90_INQ_VARID (Cforcing_id,'points',vid)
1097	ier = NF90_PUT_VAR (Cforcing_id,vid, &
1098	& (/(REAL(i,r_std),i=1,nbp_glo)/))
1099	ier = NF90_INQ_VARID (Cforcing_id,'carbtype',vid)
1100	ier = NF90_PUT_VAR (Cforcing_id,vid, &
1101	& (/(REAL(i,r_std),i=1,ncarb)/))
1102	ier = NF90_INQ_VARID (Cforcing_id,'vegtype',vid)
1103	ier = NF90_PUT_VAR (Cforcing_id,vid, &
1104	& (/(REAL(i,r_std),i=1,nvm)/))
1105	ier = NF90_INQ_VARID (Cforcing_id,'level',vid)
1106	ier = NF90_PUT_VAR (Cforcing_id,vid, &
1107	& (/(REAL(i,r_std),i=1,nlevs)/))
1108	ier = NF90_INQ_VARID (Cforcing_id,'time_step',vid)
1109	ier = NF90_PUT_VAR (Cforcing_id,vid, &
1110	& (/(REAL(i,r_std),i=1,nparan*nbyear)/))
1111	ier = NF90_INQ_VARID (Cforcing_id,'index',vid)
1112	ier = NF90_PUT_VAR (Cforcing_id,vid, REAL(index_g,r_std) )
1113	ier = NF90_INQ_VARID (Cforcing_id,'clay',vid)
1114	ier = NF90_PUT_VAR (Cforcing_id,vid, clay_g )
1115	ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',vid)
1116	ier = NF90_PUT_VAR (Cforcing_id,vid, soilcarbon_input_g )
1117	ier = NF90_INQ_VARID (Cforcing_id,'control_moist',vid)
1118	ier = NF90_PUT_VAR (Cforcing_id,vid, control_moist_g )
1119	ier = NF90_INQ_VARID (Cforcing_id,'control_temp',vid)
1120	ier = NF90_PUT_VAR (Cforcing_id,vid, control_temp_g )
1121	ier = NF90_INQ_VARID (Cforcing_id,'npp_equil',vid)
1122	ier = NF90_PUT_VAR (Cforcing_id,vid, npp_equil_g )
1123	!-
1124	ier = NF90_CLOSE (Cforcing_id)
1125	Cforcing_id = -1
1126	ENDIF
1127
1128	IF (is_root_prc) THEN
1129	DEALLOCATE(soilcarbon_input_g)
1130	DEALLOCATE(control_moist_g)
1131	DEALLOCATE(control_temp_g)
1132	DEALLOCATE(npp_equil_g)
1133	ENDIF
1134	ENDIF
1135	RETURN
1136	ENDIF ! write restart-files
1137	!
1138	! 4 Special treatment for some input arrays.
1139	!
1140	!
1141	! 4.1 Sum of liquid and solid precipitation
1142	!
1143	precip = ( precip_rain+precip_snow )*one_day/dtradia
1144	!
1145	! 4.2 Copy.
1146	!
1147	! 4.2.1 calculate STOMATE's vegetation fractions
1148	! from veget and veget_max
1149	!
1150	DO j=1,nvm
1151	WHERE ((1.-totfrac_nobio(:)) > min_sechiba)
1152	veget_cov(:,j)=veget(:,j)/( 1.-totfrac_nobio(:) )
1153	veget_cov_max(:,j)=veget_max(:,j)/( 1.-totfrac_nobio(:) )
1154	ELSEWHERE
1155	veget_cov(:,j) = zero
1156	veget_cov_max(:,j) = zero
1157	ENDWHERE
1158	ENDDO
1159	IF ( EndOfYear ) THEN
1160	DO j=1,nvm
1161	WHERE ((1.-totfrac_nobio(:)) > min_sechiba)
1162	veget_cov_max_new(:,j)=veget_max_new(:,j)/( 1.-totfrac_nobio(:) )
1163	ELSEWHERE
1164	veget_cov_max_new(:,j) = zero
1165	ENDWHERE
1166	ENDDO
1167	ENDIF
1168	gpp_d(:,1) = zero
1169	DO j = 2,nvm
1170	! 4.2.2 GPP
1171	! gpp in gC/m**2 of total ground/day
1172	WHERE (veget_cov_max(:,j) > min_stomate)
1173	gpp_d(:,j) = gpp(:,j)/ veget_cov_max(:,j)* one_day/dtradia
1174	ELSEWHERE
1175	gpp_d(:,j) = zero
1176	ENDWHERE
1177	ENDDO
1178
1179	IF ( day == 1 .AND. sec .LT. dtradia ) THEN
1180	EndOfMonth=.TRUE.
1181	ELSE
1182	EndOfMonth=.FALSE.
1183	ENDIF
1184	!
1185	! 5 "daily" variables
1186	! Note: If dt_days /= 1, then xx_daily are not daily variables,
1187	! but that is not really a problem.
1188	!
1189	!
1190	! 5.1 accumulate instantaneous variables
1191	! and eventually calculate daily mean value
1192	!
1193
1194	CALL stomate_accu (kjpindex, nvm, dt_slow, dtradia, &
1195	& do_slow, humrel, humrel_daily)
1196	CALL stomate_accu (kjpindex, 1, dt_slow, dtradia, &
1197	& do_slow, litterhumdiag, litterhum_daily)
1198	CALL stomate_accu (kjpindex, 1, dt_slow, dtradia, &
1199	& do_slow, t2m, t2m_daily)
1200	CALL stomate_accu (kjpindex, 1, dt_slow, dtradia, &
1201	& do_slow, temp_sol, tsurf_daily)
1202	CALL stomate_accu (kjpindex, nbdl, dt_slow, dtradia, &
1203	& do_slow, stempdiag, tsoil_daily)
1204	CALL stomate_accu (kjpindex, nbdl, dt_slow, dtradia, &
1205	& do_slow, shumdiag, soilhum_daily)
1206	CALL stomate_accu (kjpindex, 1, dt_slow, dtradia, &
1207	& do_slow, precip, precip_daily)
1208	CALL stomate_accu (kjpindex, nvm, dt_slow, dtradia, &
1209	& do_slow, gpp_d, gpp_daily)
1210
1211	!
1212	! 5.2 daily minimum temperature
1213	!
1214	t2m_min_daily(:) = MIN( t2m_min(:), t2m_min_daily(:) )
1215	!
1216	! 5.3 calculate respiration (NV 14/5/2002)
1217	!
1218	! 5.3.1 calculate maintenance respiration
1219	!
1220	!NV maintenant lai est passé comme argument car on n'a plus les problemes de nat/agri!
1221	!donc plus à etre calculé ici....
1222	CALL maint_respiration &
1223	& (kjpindex,dtradia,lai,t2m,tlong_ref,stempdiag,height,veget_cov_max, &
1224	& rprof,biomass,resp_maint_part_radia)
1225
1226	resp_maint_radia(:,:) = zero
1227	DO j=2,nvm
1228	DO k= 1, nparts
1229	resp_maint_radia(:,j) = resp_maint_radia(:,j) &
1230	& +resp_maint_part_radia(:,j,k)
1231	ENDDO
1232	ENDDO
1233	resp_maint_part(:,:,:)= resp_maint_part(:,:,:) &
1234	& +resp_maint_part_radia(:,:,:)
1235	!
1236	! 5.3.2 the whole litter stuff:
1237	! litter update, lignin content, PFT parts, litter decay,
1238	! litter heterotrophic respiration, dead leaf soil cover.
1239	! No vertical discretisation in the soil for litter decay.
1240	!
1241	turnover_littercalc = turnover_daily * dtradia/one_day
1242	bm_to_littercalc = bm_to_litter*dtradia/one_day
1243
1244	CALL littercalc (kjpindex, dtradia/one_day, &
1245	turnover_littercalc, bm_to_littercalc, &
1246	veget_cov_max, temp_sol, stempdiag, shumdiag, litterhumdiag, &
1247	litterpart, litter, dead_leaves, lignin_struc, &
1248	deadleaf_cover, resp_hetero_litter, &
1249	soilcarbon_input_inst, control_temp_inst, control_moist_inst)
1250	resp_hetero_litter=resp_hetero_litter*dtradia/one_day
1251	!
1252	! 5.3.3 soil carbon dynamics: heterotrophic respiration from the soil.
1253	! For the moment, no vertical discretisation.
1254	! We might later introduce a vertical discretisation.
1255	! However, in that case, we would have to treat the vertical
1256	! exchanges of carbon between the different levels.
1257	!
1258	CALL soilcarbon (kjpindex, dtradia/one_day, clay, &
1259	soilcarbon_input_inst, control_temp_inst, control_moist_inst, &
1260	carbon, resp_hetero_soil)
1261	resp_hetero_soil=resp_hetero_soil*dtradia/one_day
1262	resp_hetero_radia = resp_hetero_litter+resp_hetero_soil
1263	resp_hetero_d = resp_hetero_d + resp_hetero_radia
1264	CALL stomate_accu (kjpindex, nlevs, dt_slow, dtradia, &
1265	& do_slow, control_moist_inst, control_moist_daily)
1266	CALL stomate_accu (kjpindex, nlevs, dt_slow, dtradia, &
1267	& do_slow, control_temp_inst, control_temp_daily)
1268	DO i=1,ncarb
1269	CALL stomate_accu (kjpindex, nvm, dt_slow, dtradia, &
1270	& do_slow, soilcarbon_input_inst(:,i,:), soilcarbon_input_daily(:,i,:))
1271	ENDDO
1272	!
1273	! 6 Daily processes
1274	!
1275	IF (do_slow) THEN
1276
1277	! 6.0 update lai
1278	IF (control%ok_pheno) THEN ! lai from stomate
1279	lai(:,ibare_sechiba) = zero
1280	DO j = 2, nvm
1281	lai(:,j) = biomass(:,j,ileaf)*sla(j)
1282	ENDDO
1283	ELSE
1284	CALL setlai(kjpindex,lai) ! lai prescribed
1285	ENDIF
1286
1287	!
1288	! 6.1 total natural space
1289	!
1290	!
1291	! 6.2 Calculate longer-term "meteorological" and biological parameters
1292	!
1293	CALL season &
1294	& (kjpindex, dt_days, EndOfYear, &
1295	& veget_cov, veget_cov_max, &
1296	& humrel_daily, t2m_daily, tsoil_daily, soilhum_daily, &
1297	& precip_daily, npp_daily, biomass, &
1298	& turnover_daily, gpp_daily, when_growthinit, &
1299	& maxhumrel_lastyear, maxhumrel_thisyear, &
1300	& minhumrel_lastyear, minhumrel_thisyear, &
1301	& maxgppweek_lastyear, maxgppweek_thisyear, &
1302	& gdd0_lastyear, gdd0_thisyear, &
1303	& precip_lastyear, precip_thisyear, &
1304	& lm_lastyearmax, lm_thisyearmax, &
1305	& maxfpc_lastyear, maxfpc_thisyear, &
1306	& humrel_month, humrel_week, t2m_longterm, &
1307	& tlong_ref, t2m_month, t2m_week, tsoil_month, soilhum_month, &
1308	& npp_longterm, turnover_longterm, gpp_week, &
1309	& gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
1310	& time_lowgpp, time_hum_min, hum_min_dormance, herbivores)
1311	!
1312	! 6.3 STOMATE: allocation, phenology, etc.
1313	!
1314	IF (control%ok_stomate) THEN
1315
1316	! 6.3.1 call stomate
1317
1318	CALL StomateLpj &
1319	& (kjpindex, dt_days, EndOfYear, EndOfMonth, &
1320	& neighbours, resolution, &
1321	& clay, herbivores, &
1322	& tsurf_daily, tsoil_daily, t2m_daily, t2m_min_daily, &
1323	& litterhum_daily, soilhum_daily, &
1324	& maxhumrel_lastyear, minhumrel_lastyear, &
1325	& gdd0_lastyear, precip_lastyear, &
1326	& humrel_month, humrel_week, tlong_ref, t2m_month, t2m_week, &
1327	& tsoil_month, soilhum_month, &
1328	& gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
1329	& turnover_longterm, gpp_daily, time_lowgpp, &
1330	& time_hum_min, maxfpc_lastyear, resp_maint_part,&
1331	& PFTpresent, age, fireindex, firelitter, &
1332	& leaf_age, leaf_frac, biomass, ind, adapted, regenerate, &
1333	& senescence, when_growthinit, litterpart, litter, &
1334	& dead_leaves, carbon, black_carbon, lignin_struc, &
1335	& veget_cov_max, veget_cov, npp_longterm, lm_lastyearmax, &
1336	& veget_lastlight, everywhere, need_adjacent, RIP_time, &
1337	& lai, rprof,npp_daily, turnover_daily, turnover_time,&
1338	& control_moist_inst, control_temp_inst, soilcarbon_input_inst, &
1339	& co2_to_bm_dgvm, co2_fire, &
1340	& resp_hetero_d, resp_maint_d, resp_growth_d, &
1341	& height, deadleaf_cover, vcmax, vjmax, &
1342	& t_photo_min, t_photo_opt, t_photo_max,bm_to_litter,&
1343	& prod10, prod100, flux10, flux100, veget_cov_max_new,&
1344	& convflux, cflux_prod10, cflux_prod100, harvest_above, carb_mass_total, lcchange,&
1345	& fpc_max)
1346
1347	!
1348	! fco2_lu --> luccarb
1349	fco2_lu(:) = convflux(:) &
1350	& + cflux_prod10(:) &
1351	& + cflux_prod100(:) &
1352	& + harvest_above(:)
1353	!
1354	! 6.4 output: transform from dimension nvm to dimension nvm
1355	! Several Stomate-PFTs may correspond
1356	! to a single Sechiba-PFT (see 4.2).
1357	! We sum up the vegetation cover over these Stomate-PFTs.
1358	! Mean LAI, height, and Vmax is calculated
1359	! by ponderating with (maximum) vegetation cover.
1360	!
1361
1362	! 6.4.1 calculate veget, veget_max,
1363	! from veget_cov and veget_cov_max
1364	veget(:,:) = zero
1365	veget_max(:,:) = zero
1366	!
1367	DO j = 1, nvm
1368	veget(:,j) = veget(:,j) + &
1369	& veget_cov(:,j) * ( 1.-totfrac_nobio(:) )
1370	veget_max(:,j) = veget_max(:,j) + &
1371	& veget_cov_max(:,j) * ( 1.-totfrac_nobio(:) )
1372	ENDDO
1373	!
1374	! 6.4.2 photosynthesis parameters
1375
1376	assim_param(:,:,ivcmax) = zero
1377	assim_param(:,:,ivjmax) = zero
1378	assim_param(:,:,itmin) = zero
1379	assim_param(:,:,itopt) = zero
1380	assim_param(:,:,itmax) = zero
1381
1382	DO j = 2,nvm
1383	assim_param(:,j,ivcmax)=vcmax(:,j)
1384	ENDDO
1385
1386	DO j = 2, nvm
1387	assim_param(:,j,ivjmax)=vjmax(:,j)
1388	ENDDO
1389
1390	DO j = 2, nvm
1391	assim_param(:,j,itmin)=t_photo_min(:,j)
1392	ENDDO
1393
1394	DO j = 2, nvm
1395	assim_param(:,j,itopt)=t_photo_opt(:,j)
1396	ENDDO
1397
1398	DO j = 2, nvm
1399	assim_param(:,j,itmax)=t_photo_max(:,j)
1400	ENDDO
1401
1402	!
1403	! 6.5 update forcing variables for soil carbon
1404	!
1405	IF (ldcarbon_write .AND. TRIM(Cforcing_name) /= 'NONE') THEN
1406	npp_tot(:)=0
1407
1408	DO j=2,nvm
1409	npp_tot(:)=npp_tot(:) + npp_daily(:,j)
1410	ENDDO
1411	sf_time = MODULO(REAL(date,r_std)-1,one_year*REAL(nbyear,r_std))
1412	iatt=FLOOR(sf_time/dt_forcesoil)
1413	IF (iatt == 0) iatt = iatt_old + 1
1414
1415	IF ((iatt<iatt_old) .and. (.not. cumul_Cforcing)) THEN
1416	nforce(:)=0
1417	soilcarbon_input(:,:,:,:) = zero
1418	control_moist(:,:,:) = zero
1419	control_temp(:,:,:) = zero
1420	npp_equil(:,:) = zero
1421	ENDIF
1422	iatt_old=iatt
1423
1424
1425	nforce(iatt) = nforce(iatt)+1
1426	soilcarbon_input(:,:,:,iatt) = soilcarbon_input(:,:,:,iatt)+ soilcarbon_input_daily(:,:,:)
1427	control_moist(:,:,iatt) = control_moist(:,:,iatt)+ control_moist_daily(:,:)
1428	control_temp(:,:,iatt) = control_temp(:,:,iatt)+ control_temp_daily(:,:)
1429	npp_equil(:,iatt) = npp_equil(:,iatt)+npp_tot(:)
1430	ENDIF
1431	!
1432	! 6.6 updates qsintmax
1433	!
1434	DO j=2,nvm
1435	qsintmax(:,j) = qsintcst * veget(:,j) * lai(:,j)
1436	ENDDO
1437	ENDIF
1438	!
1439	! 6.7 write forcing file?
1440	! ldforcing_write should only be .TRUE.
1441	! if STOMATE is run in coupled mode.
1442	! In stand-alone mode, the forcing file is read!
1443	!
1444	IF ( ldforcing_write .AND. TRIM(forcing_name) /= 'NONE' ) THEN
1445
1446	gpp_daily_x(:,:) = zero
1447	!gpp needs to be multiplied by coverage for forcing (see above)
1448	DO j = 2, nvm
1449	gpp_daily_x(:,j) = gpp_daily_x(:,j) + &
1450	& gpp_daily(:,j) * dt_slow / one_day * veget_cov_max(:,j)
1451	ENDDO
1452	!
1453	! bare soil moisture availability has not been treated
1454	! in STOMATE (doesn't matter)
1455	!
1456	humrel_daily(:,ibare_sechiba) = humrel(:,ibare_sechiba)
1457
1458	! next forcing step in memory
1459	iisf = iisf+1
1460
1461	! how many times have we treated this forcing state
1462	xn = REAL(nf_cumul(isf(iisf)),r_std)
1463
1464	! cumulate. be careful :
1465	! precipitation is multiplied by dt_slow/one_day
1466	IF (cumul_forcing) THEN
1467	clay_fm(:,iisf) = (xn*clay_fm(:,iisf)+clay(:))/(xn+1.)
1468	humrel_daily_fm(:,:,iisf) = &
1469	& (xn*humrel_daily_fm(:,:,iisf) + humrel_daily(:,:))/(xn+1.)
1470	litterhum_daily_fm(:,iisf) = &
1471	& (xn*litterhum_daily_fm(:,iisf)+litterhum_daily(:))/(xn+1.)
1472	t2m_daily_fm(:,iisf) = &
1473	& (xn*t2m_daily_fm(:,iisf)+t2m_daily(:))/(xn+1.)
1474	t2m_min_daily_fm(:,iisf) = &
1475	& (xn*t2m_min_daily_fm(:,iisf)+t2m_min_daily(:))/(xn+1.)
1476	tsurf_daily_fm(:,iisf) = &
1477	& (xn*tsurf_daily_fm(:,iisf)+tsurf_daily(:))/(xn+1.)
1478	tsoil_daily_fm(:,:,iisf) = &
1479	& (xn*tsoil_daily_fm(:,:,iisf)+tsoil_daily(:,:))/(xn+1.)
1480	soilhum_daily_fm(:,:,iisf) = &
1481	& (xn*soilhum_daily_fm(:,:,iisf)+soilhum_daily(:,:))/(xn+1.)
1482	precip_fm(:,iisf) = &
1483	& (xnprecip_fm(:,iisf)+precip_daily(:)dt_slow/one_day)/(xn+1.)
1484	gpp_daily_fm(:,:,iisf) = &
1485	& (xn*gpp_daily_fm(:,:,iisf) + gpp_daily_x(:,:))/(xn+1.)
1486	veget_fm(:,:,iisf) = &
1487	& (xn*veget_fm(:,:,iisf) + veget(:,:) )/(xn+1.)
1488	veget_max_fm(:,:,iisf) = &
1489	& (xn*veget_max_fm(:,:,iisf) + veget_max(:,:) )/(xn+1.)
1490	lai_fm(:,:,iisf) = &
1491	& (xn*lai_fm(:,:,iisf) + lai(:,:) )/(xn+1.)
1492	ELSE
1493	clay_fm(:,iisf) = clay(:)
1494	humrel_daily_fm(:,:,iisf) = humrel_daily(:,:)
1495	litterhum_daily_fm(:,iisf) = litterhum_daily(:)
1496	t2m_daily_fm(:,iisf) = t2m_daily(:)
1497	t2m_min_daily_fm(:,iisf) =t2m_min_daily(:)
1498	tsurf_daily_fm(:,iisf) = tsurf_daily(:)
1499	tsoil_daily_fm(:,:,iisf) =tsoil_daily(:,:)
1500	soilhum_daily_fm(:,:,iisf) =soilhum_daily(:,:)
1501	precip_fm(:,iisf) = precip_daily(:)
1502	gpp_daily_fm(:,:,iisf) =gpp_daily_x(:,:)
1503	veget_fm(:,:,iisf) = veget(:,:)
1504	veget_max_fm(:,:,iisf) =veget_max(:,:)
1505	lai_fm(:,:,iisf) =lai(:,:)
1506	ENDIF
1507	nf_cumul(isf(iisf)) = nf_cumul(isf(iisf))+1
1508
1509	! do we have to write the forcing states?
1510	IF (iisf == nsfm) THEN
1511
1512	! write these forcing states
1513	CALL forcing_write(forcing_id,1,nsfm)
1514	! determine which forcing states must be read
1515	isf(1) = isf(nsfm)+1
1516	IF ( isf(1) > nsft ) isf(1) = 1
1517	DO iisf = 2, nsfm
1518	isf(iisf) = isf(iisf-1)+1
1519	IF (isf(iisf) > nsft) isf(iisf) = 1
1520	ENDDO
1521
1522	! read them
1523	!for debug only! CALL forcing_read(forcing_id,nsfm)
1524
1525	iisf = 0
1526
1527	ENDIF
1528
1529	ENDIF
1530
1531
1532	! 6.8 compute daily co2_flux
1533
1534	! CO2 flux in gC/m**2/sec
1535	! (positive towards the atmosphere) is sum of:
1536	! 1/ heterotrophic respiration from ground
1537	! 2/ maintenance respiration from the plants
1538	! 3/ growth respiration from the plants
1539	! 4/ co2 created by fire
1540	! 5/ - co2 taken up in the DGVM to establish saplings.
1541	! 6/ - co2 taken up by photosyntyhesis
1542
1543	!NV co2_flux devient un champs dépendant du PFT
1544	co2_flux_daily(:,:)= &
1545	& resp_maint_d(:,:) + resp_growth_d(:,:) + resp_hetero_d(:,:) + &
1546	& co2_fire(:,:) - co2_to_bm_dgvm(:,:) - gpp_daily(:,:)
1547
1548	IF ( hist_id_stom_IPCC > 0 ) THEN
1549	vartmp(:)=SUM(co2_flux_dailyveget_cov_max,dim=2)/1e3/one_daycontfrac
1550	CALL histwrite (hist_id_stom_IPCC, "nep", itime, &
1551	vartmp, kjpindex, hori_index)
1552	ENDIF
1553	!
1554	co2_flux_monthly(:,:) = co2_flux_monthly(:,:) + co2_flux_daily(:,:)
1555	! Monthly Cumulative fluxes of fluc and harvest
1556	harvest_above_monthly(:) = harvest_above_monthly(:) + harvest_above(:)
1557	cflux_prod_monthly(:) = cflux_prod_monthly(:) + convflux(:) + cflux_prod10(:) + cflux_prod100(:)
1558	IF ( EndOfMonth ) THEN
1559	IF ( control%ok_stomate ) THEN
1560	CALL histwrite (hist_id_stomate, 'CO2FLUX', itime, &
1561	co2_flux_monthly, kjpindex*nvm, horipft_index)
1562	ENDIF
1563	!MM
1564	! Si on supprimer le cumul par mois,
1565	! il ne faut pas oublier cette modif resolution(:,1)resolution(:,2)contfrac(:)
1566	DO j=2, nvm
1567	co2_flux_monthly(:,j) = co2_flux_monthly(:,j)* &
1568	resolution(:,1)resolution(:,2)contfrac(:)
1569	ENDDO
1570	!NV modif calcul du net_co2_flux
1571	net_co2_flux_monthly = zero
1572	DO ji=1,kjpindex
1573	DO j=2,nvm
1574	net_co2_flux_monthly = net_co2_flux_monthly + &
1575	& co2_flux_monthly(ji,j)*veget_cov_max(ji,j)
1576	ENDDO
1577	ENDDO
1578	! Total ( land) Cumulative fluxes of fluc and harvest
1579	net_cflux_prod_monthly_sum=&
1580	& SUM(cflux_prod_monthly(:)resolution(:,1)resolution(:,2)contfrac(:))1e-15
1581	CALL reduce_sum(net_cflux_prod_monthly_sum,net_cflux_prod_monthly_tot)
1582	CALL bcast(net_cflux_prod_monthly_tot)
1583
1584	net_harvest_above_monthly_sum=&
1585	& SUM(harvest_above_monthly(:)resolution(:,1)resolution(:,2)contfrac(:))1e-15
1586	CALL reduce_sum(net_harvest_above_monthly_sum,net_harvest_above_monthly_tot)
1587	CALL bcast(net_harvest_above_monthly_tot)
1588
1589	net_co2_flux_monthly = net_co2_flux_monthly*1e-15
1590	CALL reduce_sum(net_co2_flux_monthly,net_co2_flux_monthly_sum)
1591	CALL bcast(net_co2_flux_monthly_sum)
1592
1593	WRITE(numout,9010) 'GLOBAL net_cflux_prod_monthly (Peta gC/month) = ',net_cflux_prod_monthly_tot
1594	WRITE(numout,9010) 'GLOBAL net_harvest_above_monthly (Peta gC/month) = ',net_harvest_above_monthly_tot
1595	WRITE(numout,9010) 'GLOBAL net_co2_flux_monthly (Peta gC/month) = ',net_co2_flux_monthly_sum
1596
1597	! Calculation of net biospheric production
1598	net_biosp_prod_monthly_tot = &
1599	& ( net_co2_flux_monthly_sum + net_cflux_prod_monthly_tot + net_harvest_above_monthly_tot )
1600	WRITE(numout,9010) 'GLOBAL net_biosp_prod_monthly (Peta gC/month) = ',net_biosp_prod_monthly_tot
1601
1602	9010 FORMAT(A52,F17.14)
1603	!!$ IF ( control%ok_stomate ) THEN
1604	!!$ vartmp(:)=net_co2_flux_monthly_sum
1605	!!$ CALL histwrite (hist_id_stomate, 'CO2FLUX_MONTHLY_SUM', itime, &
1606	!!$ vartmp, kjpindex, hori_index )
1607	!!$ ENDIF
1608	!!$ IF (is_root_prc) THEN
1609	!!$ OPEN( unit=39, &
1610	!!$ file="stomate_co2flux.data", &
1611	!!$ action="write", &
1612	!!$ position="append", &
1613	!!$ iostat=ios )
1614	!!$ IF ( ios /= 0 ) THEN
1615	!!$ STOP "Erreur lors de la lecture/ecriture du fichier stomate_co2flux.data"
1616	!!$ ELSE
1617	!!$ WRITE(numout,*)
1618	!!$ WRITE(numout,*) "Ecriture du fichier stomate_co2flux.data"
1619	!!$ WRITE(numout,*)
1620	!!$ END IF
1621	!!$ WRITE(39,*) net_co2_flux_monthly_sum
1622	!!$ CLOSE( unit=39 )
1623	!!$ ENDIF
1624	co2_flux_monthly(:,:) = zero
1625	harvest_above_monthly(:) = zero
1626	cflux_prod_monthly(:) = zero
1627	ENDIF
1628	!
1629	! 6.9 reset daily variables
1630	!
1631	humrel_daily(:,:) = zero
1632	litterhum_daily(:) = zero
1633	t2m_daily(:) = zero
1634	t2m_min_daily(:) = large_value
1635	tsurf_daily(:) = zero
1636	tsoil_daily(:,:) = zero
1637	soilhum_daily(:,:) = zero
1638	precip_daily(:) = zero
1639	gpp_daily(:,:) = zero
1640	resp_maint_part(:,:,:)=zero
1641	resp_hetero_d=zero
1642	IF (bavard >= 3) THEN
1643	WRITE(numout,*) 'stomate_main: daily processes done'
1644	ENDIF
1645
1646	ENDIF ! daily processes?
1647	!
1648	! 7 Outputs from Stomate
1649	! co2_flux is assigned a value only if STOMATE is activated.
1650	! Otherwise, the calling hydrological module must do this itself.
1651	!
1652	IF ( control%ok_stomate ) THEN
1653
1654	resp_maint(:,:) = resp_maint_radia(:,:)*veget_cov_max(:,:)
1655	resp_maint(:,ibare_sechiba) = zero
1656	resp_growth(:,:)= resp_growth_d(:,:)veget_cov_max(:,:)dtradia/one_day
1657	!
1658	resp_hetero(:,:) = resp_hetero_radia(:,:)*veget_cov_max(:,:)
1659
1660	! CO2 flux in gC/m**2/sec (positive towards the atmosphere) is sum of:
1661	! 1/ heterotrophic respiration from ground
1662	! 2/ maintenance respiration from the plants
1663	! 3/ growth respiration from the plants
1664	! 4/ co2 created by fire
1665	! 5/ - co2 taken up in the DGVM to establish saplings.
1666	! 6/ - co2 taken up by photosyntyhesis
1667
1668	co2_flux(:,:) = resp_hetero(:,:) + resp_maint(:,:) + resp_growth(:,:) &
1669	& + (co2_fire(:,:)-co2_to_bm_dgvm(:,:))*veget_cov_max(:,:)/one_day &
1670	& - gpp(:,:)
1671
1672	ENDIF
1673	!
1674	! 8 message
1675	!
1676	IF ( (bavard >= 2).AND.EndOfYear.AND.do_slow) THEN
1677	WRITE(numout,*) 'stomate: EndOfYear'
1678	ENDIF
1679	IF (bavard >= 4) WRITE(numout,*) 'Leaving stomate_main'
1680	IF (long_print) WRITE (numout,*) ' stomate_main done '
1681	!--------------------------
1682	END SUBROUTINE stomate_main
1683	!
1684	!=
1685	!
1686	SUBROUTINE stomate_init &
1687	& (kjpij, kjpindex, index, ldforcing_write, lalo, &
1688	& rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
1689	!---------------------------------------------------------------------
1690	! interface description
1691	! input scalar
1692	! Total size of the un-compressed grid
1693	INTEGER(i_std),INTENT(in) :: kjpij
1694	! Domain size
1695	INTEGER(i_std),INTENT(in) :: kjpindex
1696	! Logical for _forcing_ file to write
1697	LOGICAL,INTENT(in) :: ldforcing_write
1698	! Geogr. coordinates (latitude,longitude) (degrees)
1699	REAL(r_std),DIMENSION(kjpindex,2),INTENT(in) :: lalo
1700	! STOMATE's _Restart_ file file identifier
1701	INTEGER(i_std),INTENT(in) :: rest_id_stom
1702	! STOMATE's _history_ file file identifier
1703	INTEGER(i_std),INTENT(in) :: hist_id_stom
1704	! STOMATE's IPCC _history_ file file identifier
1705	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC
1706	! input fields
1707	! Indeces of the points on the map
1708	INTEGER(i_std),DIMENSION(kjpindex),INTENT(in) :: index
1709	! local declaration
1710	REAL(r_std) :: tmp_day(1)
1711	! soil depth taken for canopy
1712	REAL(r_std) :: zcanop
1713	! soil depths at diagnostic levels
1714	REAL(r_std),DIMENSION(nbdl) :: zsoil
1715	! Index
1716	INTEGER(i_std) :: l
1717	! allocation error
1718	LOGICAL :: l_error
1719	INTEGER(i_std) :: ier
1720	! indices
1721	INTEGER(i_std) :: ji,j,ipd
1722	!---------------------------------------------------------------------
1723	!
1724	! 1 online diagnostics
1725	! (by default, "bavard" is set to 1 in stomate_constants)
1726	!
1727	!Config Key = BAVARD
1728	!Config Desc = level of online diagnostics in STOMATE (0-4)
1729	!Config Def = 1
1730	!Config Help = With this variable, you can determine
1731	!Config how much online information STOMATE
1732	!Config gives during the run. 0 means
1733	!Config virtually no info.
1734	!
1735	bavard = 1
1736	CALL getin_p('BAVARD', bavard)
1737
1738	IF ( kjpindex > 0 ) THEN
1739	!
1740	!Config Key = STOMATE_DIAGPT
1741	!Config Desc = Index of grid point for online diagnostics
1742	!Config Def = 1
1743	!Config Help = This is the index of the grid point which
1744	! will be used for online diagnostics.
1745	ipd = 1
1746	CALL getin_p('STOMATE_DIAGPT',ipd)
1747	ipd = MIN( ipd, kjpindex )
1748	WRITE(numout,*) 'Stomate: '
1749	WRITE(numout,*) ' Index of grid point for online diagnostics: ',ipd
1750	WRITE(numout,*) ' Lon, lat:',lalo(ipd,2),lalo(ipd,1)
1751	WRITE(numout,*) ' Index of this point on GCM grid: ',index(ipd)
1752	!
1753	ENDIF
1754	!
1755	! 2 check consistency of flags
1756	!
1757	IF ( ( .NOT. control%ok_stomate ) .AND. control%ok_dgvm ) THEN
1758	WRITE(numout,*) 'Cannot do dynamical vegetation without STOMATE.'
1759	WRITE(numout,*) 'We stop.'
1760	STOP
1761	ENDIF
1762
1763	IF ((.NOT.control%ok_co2).AND.control%ok_stomate) THEN
1764	WRITE(numout,*) 'Cannot call STOMATE without GPP.'
1765	WRITE(numout,*) 'We stop.'
1766	STOP
1767	ENDIF
1768
1769	IF ( ( .NOT. control%ok_co2 ) .AND. ldforcing_write ) THEN
1770	WRITE(numout,*) &
1771	& 'Cannot write forcing file if photosynthesis is not activated'
1772	WRITE(numout,*) 'We stop.'
1773	STOP
1774	ENDIF
1775	!
1776	! 3 messages
1777	!
1778	WRITE(numout,*) 'stomate: first call'
1779	WRITE(numout,*) ' Photosynthesis: ', control%ok_co2
1780	WRITE(numout,*) ' STOMATE: ', control%ok_stomate
1781	WRITE(numout,*) ' LPJ: ', control%ok_dgvm
1782	!
1783	! 4 allocation of STOMATE's variables
1784	!
1785	l_error = .FALSE.
1786	ALLOCATE(veget_cov_max(kjpindex,nvm),stat=ier)
1787	l_error = l_error .OR. (ier /= 0)
1788	ALLOCATE(veget_cov_max_new(kjpindex,nvm),stat=ier)
1789	l_error = l_error .OR. (ier /= 0)
1790	ALLOCATE(ind(kjpindex,nvm),stat=ier)
1791	l_error = l_error .OR. (ier /= 0)
1792	ALLOCATE(adapted(kjpindex,nvm),stat=ier)
1793	l_error = l_error .OR. (ier /= 0)
1794	ALLOCATE(regenerate(kjpindex,nvm),stat=ier)
1795	l_error = l_error .OR. (ier /= 0)
1796	ALLOCATE(humrel_daily(kjpindex,nvm),stat=ier)
1797	l_error = l_error .OR. (ier /= 0)
1798	ALLOCATE(litterhum_daily(kjpindex),stat=ier)
1799	l_error = l_error .OR. (ier /= 0)
1800	ALLOCATE(t2m_daily(kjpindex),stat=ier)
1801	l_error = l_error .OR. (ier /= 0)
1802	ALLOCATE(t2m_min_daily(kjpindex),stat=ier)
1803	l_error = l_error .OR. (ier /= 0)
1804	ALLOCATE(tsurf_daily(kjpindex),stat=ier)
1805	l_error = l_error .OR. (ier /= 0)
1806	ALLOCATE(tsoil_daily(kjpindex,nbdl),stat=ier)
1807	l_error = l_error .OR. (ier /= 0)
1808	ALLOCATE(soilhum_daily(kjpindex,nbdl),stat=ier)
1809	l_error = l_error .OR. (ier /= 0)
1810	ALLOCATE(precip_daily(kjpindex),stat=ier)
1811	l_error = l_error .OR. (ier /= 0)
1812	ALLOCATE(gpp_daily(kjpindex,nvm),stat=ier)
1813	l_error = l_error .OR. (ier /= 0)
1814	ALLOCATE(npp_daily(kjpindex,nvm),stat=ier)
1815	l_error = l_error .OR. (ier /= 0)
1816	ALLOCATE(turnover_daily(kjpindex,nvm,nparts),stat=ier)
1817	l_error = l_error .OR. (ier /= 0)
1818	ALLOCATE(turnover_littercalc(kjpindex,nvm,nparts),stat=ier)
1819	l_error = l_error .OR. (ier /= 0)
1820	ALLOCATE(humrel_month(kjpindex,nvm),stat=ier)
1821	l_error = l_error .OR. (ier /= 0)
1822	ALLOCATE(humrel_week(kjpindex,nvm),stat=ier)
1823	l_error = l_error .OR. (ier /= 0)
1824	ALLOCATE(t2m_longterm(kjpindex),stat=ier)
1825	l_error = l_error .OR. (ier /= 0)
1826	ALLOCATE(tlong_ref(kjpindex),stat=ier)
1827	l_error = l_error .OR. (ier /= 0)
1828	ALLOCATE(t2m_month(kjpindex),stat=ier)
1829	l_error = l_error .OR. (ier /= 0)
1830	ALLOCATE(t2m_week(kjpindex),stat=ier)
1831	l_error = l_error .OR. (ier /= 0)
1832	ALLOCATE(tsoil_month(kjpindex,nbdl),stat=ier)
1833	l_error = l_error .OR. (ier /= 0)
1834	ALLOCATE(soilhum_month(kjpindex,nbdl),stat=ier)
1835	l_error = l_error .OR. (ier /= 0)
1836	ALLOCATE(fireindex(kjpindex,nvm),stat=ier)
1837	l_error = l_error .OR. (ier /= 0)
1838	ALLOCATE(firelitter(kjpindex,nvm),stat=ier)
1839	l_error = l_error .OR. (ier /= 0)
1840	ALLOCATE(maxhumrel_lastyear(kjpindex,nvm),stat=ier)
1841	l_error = l_error .OR. (ier /= 0)
1842	ALLOCATE(maxhumrel_thisyear(kjpindex,nvm),stat=ier)
1843	l_error = l_error .OR. (ier /= 0)
1844	ALLOCATE(minhumrel_lastyear(kjpindex,nvm),stat=ier)
1845	l_error = l_error .OR. (ier /= 0)
1846	ALLOCATE(minhumrel_thisyear(kjpindex,nvm),stat=ier)
1847	l_error = l_error .OR. (ier /= 0)
1848	ALLOCATE(maxgppweek_lastyear(kjpindex,nvm),stat=ier)
1849	l_error = l_error .OR. (ier /= 0)
1850	ALLOCATE(maxgppweek_thisyear(kjpindex,nvm),stat=ier)
1851	l_error = l_error .OR. (ier /= 0)
1852	ALLOCATE(gdd0_lastyear(kjpindex),stat=ier)
1853	l_error = l_error .OR. (ier /= 0)
1854	ALLOCATE(gdd0_thisyear(kjpindex),stat=ier)
1855	l_error = l_error .OR. (ier /= 0)
1856	ALLOCATE(precip_lastyear(kjpindex),stat=ier)
1857	l_error = l_error .OR. (ier /= 0)
1858	ALLOCATE(precip_thisyear(kjpindex),stat=ier)
1859	l_error = l_error .OR. (ier /= 0)
1860	ALLOCATE(gdd_m5_dormance(kjpindex,nvm),stat=ier)
1861	l_error = l_error .OR. (ier /= 0)
1862	ALLOCATE(gdd_midwinter(kjpindex,nvm),stat=ier)
1863	l_error = l_error .OR. (ier /= 0)
1864	ALLOCATE(ncd_dormance(kjpindex,nvm),stat=ier)
1865	l_error = l_error .OR. (ier /= 0)
1866	ALLOCATE(ngd_minus5(kjpindex,nvm),stat=ier)
1867	l_error = l_error .OR. (ier /= 0)
1868	ALLOCATE(PFTpresent(kjpindex,nvm),stat=ier)
1869	l_error = l_error .OR. (ier /= 0)
1870	ALLOCATE(npp_longterm(kjpindex,nvm),stat=ier)
1871	l_error = l_error .OR. (ier /= 0)
1872	ALLOCATE(lm_lastyearmax(kjpindex,nvm),stat=ier)
1873	l_error = l_error .OR. (ier /= 0)
1874	ALLOCATE(lm_thisyearmax(kjpindex,nvm),stat=ier)
1875	l_error = l_error .OR. (ier /= 0)
1876	ALLOCATE(maxfpc_lastyear(kjpindex,nvm),stat=ier)
1877	l_error = l_error .OR. (ier /= 0)
1878	ALLOCATE(maxfpc_thisyear(kjpindex,nvm),stat=ier)
1879	l_error = l_error .OR. (ier /= 0)
1880	ALLOCATE(turnover_longterm(kjpindex,nvm,nparts),stat=ier)
1881	l_error = l_error .OR. (ier /= 0)
1882	ALLOCATE(gpp_week(kjpindex,nvm),stat=ier)
1883	l_error = l_error .OR. (ier /= 0)
1884	ALLOCATE(biomass(kjpindex,nvm,nparts),stat=ier)
1885	l_error = l_error .OR. (ier /= 0)
1886	ALLOCATE(senescence(kjpindex,nvm),stat=ier)
1887	l_error = l_error .OR. (ier /= 0)
1888	ALLOCATE(when_growthinit(kjpindex,nvm),stat=ier)
1889	l_error = l_error .OR. (ier /= 0)
1890	ALLOCATE(age(kjpindex,nvm),stat=ier)
1891	l_error = l_error .OR. (ier /= 0)
1892	ALLOCATE(resp_hetero_d(kjpindex,nvm),stat=ier)
1893	l_error = l_error .OR. (ier /= 0)
1894	ALLOCATE(resp_hetero_radia(kjpindex,nvm),stat=ier)
1895	l_error = l_error .OR. (ier /= 0)
1896	ALLOCATE(resp_maint_d(kjpindex,nvm),stat=ier)
1897	l_error = l_error .OR. (ier /= 0)
1898	ALLOCATE(resp_growth_d(kjpindex,nvm),stat=ier)
1899	l_error = l_error .OR. (ier /= 0)
1900	!NV devient 2D
1901	ALLOCATE(co2_fire(kjpindex,nvm),stat=ier)
1902	l_error = l_error .OR. (ier /= 0)
1903	!NV devient 2D
1904	ALLOCATE(co2_to_bm_dgvm(kjpindex,nvm),stat=ier)
1905	l_error = l_error .OR. (ier /= 0)
1906	ALLOCATE(veget_lastlight(kjpindex,nvm),stat=ier)
1907	l_error = l_error .OR. (ier /= 0)
1908	ALLOCATE(everywhere(kjpindex,nvm),stat=ier)
1909	l_error = l_error .OR. (ier /= 0)
1910	ALLOCATE(need_adjacent(kjpindex,nvm),stat=ier)
1911	l_error = l_error .OR. (ier /= 0)
1912	ALLOCATE(leaf_age(kjpindex,nvm,nleafages),stat=ier)
1913	l_error = l_error .OR. (ier /= 0)
1914	ALLOCATE(leaf_frac(kjpindex,nvm,nleafages),stat=ier)
1915	l_error = l_error .OR. (ier /= 0)
1916	ALLOCATE(RIP_time(kjpindex,nvm),stat=ier)
1917	l_error = l_error .OR. (ier /= 0)
1918	ALLOCATE(time_lowgpp(kjpindex,nvm),stat=ier)
1919	l_error = l_error .OR. (ier /= 0)
1920	ALLOCATE(time_hum_min(kjpindex,nvm),stat=ier)
1921	l_error = l_error .OR. (ier /= 0)
1922	ALLOCATE(hum_min_dormance(kjpindex,nvm),stat=ier)
1923	l_error = l_error .OR. (ier /= 0)
1924	ALLOCATE(litterpart(kjpindex,nvm,nlitt),stat=ier)
1925	l_error = l_error .OR. (ier /= 0)
1926	ALLOCATE(litter(kjpindex,nlitt,nvm,nlevs),stat=ier)
1927	l_error = l_error .OR. (ier /= 0)
1928	ALLOCATE(dead_leaves(kjpindex,nvm,nlitt),stat=ier)
1929	l_error = l_error .OR. (ier /= 0)
1930	ALLOCATE(carbon(kjpindex,ncarb,nvm),stat=ier)
1931	l_error = l_error .OR. (ier /= 0)
1932	ALLOCATE(black_carbon(kjpindex),stat=ier)
1933	l_error = l_error .OR. (ier /= 0)
1934	ALLOCATE(lignin_struc(kjpindex,nvm,nlevs),stat=ier)
1935	l_error = l_error .OR. (ier /= 0)
1936	ALLOCATE(turnover_time(kjpindex,nvm),stat=ier)
1937	l_error = l_error .OR. (ier /= 0)
1938	ALLOCATE(co2_flux_daily(kjpindex,nvm),stat=ier)
1939	l_error = l_error .OR. (ier /= 0)
1940	ALLOCATE(co2_flux_monthly(kjpindex,nvm),stat=ier)
1941	l_error = l_error .OR. (ier /= 0)
1942	ALLOCATE (cflux_prod_monthly(kjpindex), stat=ier)
1943	l_error = l_error .OR. (ier.NE.0)
1944	ALLOCATE (harvest_above_monthly(kjpindex), stat=ier)
1945	l_error = l_error .OR. (ier.NE.0)
1946	ALLOCATE(bm_to_litter(kjpindex,nvm,nparts),stat=ier)
1947	l_error = l_error .OR. (ier /= 0)
1948	ALLOCATE(bm_to_littercalc(kjpindex,nvm,nparts),stat=ier)
1949	l_error = l_error .OR. (ier /= 0)
1950	ALLOCATE(herbivores(kjpindex,nvm),stat=ier)
1951	l_error = l_error .OR. (ier /= 0)
1952	ALLOCATE(hori_index(kjpindex),stat=ier)
1953	l_error = l_error .OR. (ier /= 0)
1954	ALLOCATE(horipft_index(kjpindex*nvm),stat=ier)
1955	l_error = l_error .OR. (ier /= 0)
1956	ALLOCATE(resp_maint_part_radia(kjpindex,nvm,nparts),stat=ier)
1957	l_error = l_error .OR. (ier /= 0)
1958	ALLOCATE(resp_maint_radia(kjpindex,nvm),stat=ier)
1959	l_error = l_error .OR. (ier /= 0)
1960	ALLOCATE(resp_maint_part(kjpindex,nvm,nparts),stat=ier)
1961	l_error = l_error .OR. (ier /= 0)
1962	resp_maint_part(:,:,:)=zero
1963	ALLOCATE (horip10_index(kjpindex*10), stat=ier)
1964	l_error = l_error .OR. (ier.NE.0)
1965	ALLOCATE (horip100_index(kjpindex*100), stat=ier)
1966	l_error = l_error .OR. (ier.NE.0)
1967	ALLOCATE (horip11_index(kjpindex*11), stat=ier)
1968	l_error = l_error .OR. (ier.NE.0)
1969	ALLOCATE (horip101_index(kjpindex*101), stat=ier)
1970	l_error = l_error .OR. (ier.NE.0)
1971	ALLOCATE (prod10(kjpindex,0:10), stat=ier)
1972	l_error = l_error .OR. (ier.NE.0)
1973	ALLOCATE (prod100(kjpindex,0:100), stat=ier)
1974	l_error = l_error .OR. (ier.NE.0)
1975	ALLOCATE (flux10(kjpindex,10), stat=ier)
1976	l_error = l_error .OR. (ier.NE.0)
1977	ALLOCATE (flux100(kjpindex,100), stat=ier)
1978	l_error = l_error .OR. (ier.NE.0)
1979	ALLOCATE (convflux(kjpindex), stat=ier)
1980	l_error = l_error .OR. (ier.NE.0)
1981	ALLOCATE (cflux_prod10(kjpindex), stat=ier)
1982	l_error = l_error .OR. (ier.NE.0)
1983	ALLOCATE (cflux_prod100(kjpindex), stat=ier)
1984	l_error = l_error .OR. (ier.NE.0)
1985	ALLOCATE (harvest_above(kjpindex), stat=ier)
1986	l_error = l_error .OR. (ier.NE.0)
1987	ALLOCATE (carb_mass_total(kjpindex), stat=ier)
1988	l_error = l_error .OR. (ier.NE.0)
1989	ALLOCATE (soilcarbon_input_daily(kjpindex,ncarb,nvm), stat=ier)
1990	l_error = l_error .OR. (ier.NE.0)
1991	ALLOCATE (control_temp_daily(kjpindex,nlevs), stat=ier)
1992	l_error = l_error .OR. (ier.NE.0)
1993	ALLOCATE (control_moist_daily(kjpindex,nlevs), stat=ier)
1994	l_error = l_error .OR. (ier.NE.0)
1995	!
1996	ALLOCATE (fpc_max(kjpindex,nvm), stat=ier)
1997	l_error = l_error .OR. (ier.NE.0)
1998	!
1999	IF (l_error) THEN
2000	STOP 'stomate_init: error in memory allocation'
2001	ENDIF
2002	!
2003	! 5 file definitions: stored in common variables
2004	!
2005	hist_id_stomate = hist_id_stom
2006	hist_id_stomate_IPCC = hist_id_stom_IPCC
2007	rest_id_stomate = rest_id_stom
2008	hori_index(:) = index(:)
2009
2010	! Get the indexing table for the vegetation fields.
2011	! In STOMATE we work on
2012	! reduced grids but to store in the full 3D filed vegetation variable
2013	! we need another index table : indexpft
2014
2015	DO j = 1, nvm
2016	DO ji = 1, kjpindex
2017	horipft_index((j-1)kjpindex+ji) = index(ji)+(j-1)kjpij
2018	ENDDO
2019	ENDDO
2020
2021	! indexing tables added for land cover change fields
2022	DO j = 1, 10
2023	DO ji = 1, kjpindex
2024	horip10_index((j-1)kjpindex+ji) = index(ji)+(j-1)kjpij
2025	ENDDO
2026	ENDDO
2027
2028	DO j = 1, 100
2029	DO ji = 1, kjpindex
2030	horip100_index((j-1)kjpindex+ji) = index(ji)+(j-1)kjpij
2031	ENDDO
2032	ENDDO
2033
2034	DO j = 1, 11
2035	DO ji = 1, kjpindex
2036	horip11_index((j-1)kjpindex+ji) = index(ji)+(j-1)kjpij
2037	ENDDO
2038	ENDDO
2039
2040	DO j = 1, 101
2041	DO ji = 1, kjpindex
2042	horip101_index((j-1)kjpindex+ji) = index(ji)+(j-1)kjpij
2043	ENDDO
2044	ENDDO
2045	!
2046	! 6 some flags
2047	!
2048	!
2049	!Config Key = HERBIVORES
2050	!Config Desc = herbivores allowed?
2051	!Config Def = n
2052	!Config Help = With this variable, you can determine
2053	!Config if herbivores are activated
2054	!
2055	ok_herbivores = .FALSE.
2056	CALL getin_p('HERBIVORES', ok_herbivores)
2057	!
2058	WRITE(numout,*) 'herbivores activated: ',ok_herbivores
2059	!
2060	!Config Key = TREAT_EXPANSION
2061	!Config Desc = treat expansion of PFTs across a grid cell?
2062	!Config Def = n
2063	!Config Help = With this variable, you can determine
2064	!Config whether we treat expansion of PFTs across a
2065	!Config grid cell.
2066	!
2067	treat_expansion = .FALSE.
2068	CALL getin_p('TREAT_EXPANSION', treat_expansion)
2069	!
2070	WRITE(numout,*) &
2071	& 'expansion across a grid cell is treated: ',treat_expansion
2072
2073	!Config Key = LPJ_GAP_CONST_MORT
2074	!Config Desc = prescribe mortality if not using DGVM?
2075	!Config Def = y
2076	!Config Help = set to TRUE if constant mortality is to be activated
2077	! ignored if DGVM=true!
2078	!
2079	lpj_gap_const_mort=.TRUE.
2080	CALL getin('LPJ_GAP_CONST_MORT', lpj_gap_const_mort)
2081	WRITE(numout,*) 'LPJ GAP: constant mortality:', lpj_gap_const_mort
2082
2083	!Config Key = HARVEST_AGRI
2084	!Config Desc = Harvert model for agricol PFTs.
2085	!Config Def = y
2086	!Config Help = Compute harvest above ground biomass for agriculture.
2087	!Config Change daily turnover.
2088	harvest_agri = .TRUE.
2089	CALL getin_p('HARVEST_AGRI', harvest_agri)
2090
2091	!
2092	! 7 some global initializations
2093	!
2094	! edit shilong
2095	! bm_to_litter(:,:,:) = zero
2096	turnover_daily(:,:,:) = zero
2097	resp_hetero_d(:,:) = zero
2098	co2_flux_daily(:,:) = zero
2099	co2_flux_monthly(:,:) = zero
2100	cflux_prod_monthly(:) = zero
2101	harvest_above_monthly(:) = zero
2102	control_moist_daily(:,:) = zero
2103	control_temp_daily(:,:) = zero
2104	soilcarbon_input_daily(:,:,:) = zero
2105
2106	! initialisation of land cover change variables
2107	prod10(:,:) = zero
2108	prod100(:,:) = zero
2109	flux10(:,:) = zero
2110	flux100(:,:) = zero
2111	convflux(:) = zero
2112	cflux_prod10(:) = zero
2113	cflux_prod100(:)= zero
2114
2115	fpc_max(:,:)=zero
2116	!--------------------------
2117	END SUBROUTINE stomate_init
2118	!
2119	!=
2120	!
2121	SUBROUTINE stomate_clear
2122	!---------------------------------------------------------------------
2123	! 1. Deallocate all dynamics variables
2124	IF (ALLOCATED(veget_cov_max)) DEALLOCATE(veget_cov_max)
2125	IF (ALLOCATED(veget_cov_max_new)) DEALLOCATE(veget_cov_max_new)
2126	IF (ALLOCATED(ind)) DEALLOCATE(ind)
2127	IF (ALLOCATED(adapted)) DEALLOCATE(adapted)
2128	IF (ALLOCATED(regenerate)) DEALLOCATE(regenerate)
2129	IF (ALLOCATED(humrel_daily)) DEALLOCATE(humrel_daily)
2130	IF (ALLOCATED(litterhum_daily)) DEALLOCATE(litterhum_daily)
2131	IF (ALLOCATED(t2m_daily)) DEALLOCATE(t2m_daily)
2132	IF (ALLOCATED(t2m_min_daily)) DEALLOCATE(t2m_min_daily)
2133	IF (ALLOCATED(tsurf_daily)) DEALLOCATE(tsurf_daily)
2134	IF (ALLOCATED(tsoil_daily)) DEALLOCATE(tsoil_daily)
2135	IF (ALLOCATED(soilhum_daily)) DEALLOCATE(soilhum_daily)
2136	IF (ALLOCATED(precip_daily)) DEALLOCATE(precip_daily)
2137	IF (ALLOCATED(gpp_daily)) DEALLOCATE(gpp_daily)
2138	IF (ALLOCATED(npp_daily)) DEALLOCATE(npp_daily)
2139	IF (ALLOCATED(turnover_daily)) DEALLOCATE(turnover_daily)
2140	IF (ALLOCATED(turnover_littercalc)) DEALLOCATE(turnover_littercalc)
2141	IF (ALLOCATED(humrel_month)) DEALLOCATE(humrel_month)
2142	IF (ALLOCATED(humrel_week)) DEALLOCATE(humrel_week)
2143	IF (ALLOCATED(t2m_longterm)) DEALLOCATE(t2m_longterm)
2144	IF (ALLOCATED(tlong_ref)) DEALLOCATE(tlong_ref)
2145	IF (ALLOCATED(t2m_month)) DEALLOCATE(t2m_month)
2146	IF (ALLOCATED(t2m_week)) DEALLOCATE(t2m_week)
2147	IF (ALLOCATED(tsoil_month)) DEALLOCATE(tsoil_month)
2148	IF (ALLOCATED(soilhum_month)) DEALLOCATE(soilhum_month)
2149	IF (ALLOCATED(fireindex)) DEALLOCATE(fireindex)
2150	IF (ALLOCATED(firelitter)) DEALLOCATE(firelitter)
2151	IF (ALLOCATED(maxhumrel_lastyear)) DEALLOCATE(maxhumrel_lastyear)
2152	IF (ALLOCATED(maxhumrel_thisyear)) DEALLOCATE(maxhumrel_thisyear)
2153	IF (ALLOCATED(minhumrel_lastyear)) DEALLOCATE(minhumrel_lastyear)
2154	IF (ALLOCATED(minhumrel_thisyear)) DEALLOCATE(minhumrel_thisyear)
2155	IF (ALLOCATED(maxgppweek_lastyear)) DEALLOCATE(maxgppweek_lastyear)
2156	IF (ALLOCATED(maxgppweek_thisyear)) DEALLOCATE(maxgppweek_thisyear)
2157	IF (ALLOCATED(gdd0_lastyear)) DEALLOCATE(gdd0_lastyear)
2158	IF (ALLOCATED(gdd0_thisyear)) DEALLOCATE(gdd0_thisyear)
2159	IF (ALLOCATED(precip_lastyear)) DEALLOCATE(precip_lastyear)
2160	IF (ALLOCATED(precip_thisyear)) DEALLOCATE(precip_thisyear)
2161	IF (ALLOCATED(gdd_m5_dormance)) DEALLOCATE(gdd_m5_dormance)
2162	IF (ALLOCATED(gdd_midwinter)) DEALLOCATE(gdd_midwinter)
2163	IF (ALLOCATED(ncd_dormance)) DEALLOCATE(ncd_dormance)
2164	IF (ALLOCATED(ngd_minus5)) DEALLOCATE(ngd_minus5)
2165	IF (ALLOCATED(PFTpresent)) DEALLOCATE(PFTpresent)
2166	IF (ALLOCATED(npp_longterm)) DEALLOCATE(npp_longterm)
2167	IF (ALLOCATED(lm_lastyearmax)) DEALLOCATE(lm_lastyearmax)
2168	IF (ALLOCATED(lm_thisyearmax)) DEALLOCATE(lm_thisyearmax)
2169	IF (ALLOCATED(maxfpc_lastyear)) DEALLOCATE(maxfpc_lastyear)
2170	IF (ALLOCATED(maxfpc_thisyear)) DEALLOCATE(maxfpc_thisyear)
2171	IF (ALLOCATED(turnover_longterm)) DEALLOCATE(turnover_longterm)
2172	IF (ALLOCATED(gpp_week)) DEALLOCATE(gpp_week)
2173	IF (ALLOCATED(biomass)) DEALLOCATE(biomass)
2174	IF (ALLOCATED(senescence)) DEALLOCATE(senescence)
2175	IF (ALLOCATED(when_growthinit)) DEALLOCATE(when_growthinit)
2176	IF (ALLOCATED(age)) DEALLOCATE(age)
2177	IF (ALLOCATED(resp_hetero_d)) DEALLOCATE(resp_hetero_d)
2178	IF (ALLOCATED(resp_hetero_radia)) DEALLOCATE(resp_hetero_radia)
2179	IF (ALLOCATED(resp_maint_d)) DEALLOCATE(resp_maint_d)
2180	IF (ALLOCATED(resp_growth_d)) DEALLOCATE(resp_growth_d)
2181	IF (ALLOCATED(co2_fire)) DEALLOCATE(co2_fire)
2182	IF (ALLOCATED(co2_to_bm_dgvm)) DEALLOCATE(co2_to_bm_dgvm)
2183	IF (ALLOCATED(veget_lastlight)) DEALLOCATE(veget_lastlight)
2184	IF (ALLOCATED(everywhere)) DEALLOCATE(everywhere)
2185	IF (ALLOCATED(need_adjacent)) DEALLOCATE(need_adjacent)
2186	IF (ALLOCATED(leaf_age)) DEALLOCATE(leaf_age)
2187	IF (ALLOCATED(leaf_frac)) DEALLOCATE(leaf_frac)
2188	IF (ALLOCATED(RIP_time)) DEALLOCATE(RIP_time)
2189	IF (ALLOCATED(time_lowgpp)) DEALLOCATE(time_lowgpp)
2190	IF (ALLOCATED(time_hum_min)) DEALLOCATE(time_hum_min)
2191	IF (ALLOCATED(hum_min_dormance)) DEALLOCATE(hum_min_dormance)
2192	IF (ALLOCATED(litterpart)) DEALLOCATE(litterpart)
2193	IF (ALLOCATED(litter)) DEALLOCATE(litter)
2194	IF (ALLOCATED(dead_leaves)) DEALLOCATE(dead_leaves)
2195	IF (ALLOCATED(carbon)) DEALLOCATE(carbon)
2196	IF (ALLOCATED(black_carbon)) DEALLOCATE(black_carbon)
2197	IF (ALLOCATED(lignin_struc)) DEALLOCATE(lignin_struc)
2198	IF (ALLOCATED(turnover_time)) DEALLOCATE(turnover_time)
2199	IF (ALLOCATED(co2_flux_daily)) DEALLOCATE(co2_flux_daily)
2200	IF (ALLOCATED(co2_flux_monthly)) DEALLOCATE(co2_flux_monthly)
2201	IF (ALLOCATED(harvest_above_monthly)) DEALLOCATE (harvest_above_monthly)
2202	IF (ALLOCATED(cflux_prod_monthly)) DEALLOCATE (cflux_prod_monthly)
2203	IF (ALLOCATED(bm_to_litter)) DEALLOCATE(bm_to_litter)
2204	IF (ALLOCATED(bm_to_littercalc)) DEALLOCATE(bm_to_littercalc)
2205	IF (ALLOCATED(herbivores)) DEALLOCATE(herbivores)
2206	IF (ALLOCATED(resp_maint_part_radia)) DEALLOCATE(resp_maint_part_radia)
2207	IF (ALLOCATED(resp_maint_radia)) DEALLOCATE(resp_maint_radia)
2208	IF (ALLOCATED(resp_maint_part)) DEALLOCATE(resp_maint_part)
2209	IF (ALLOCATED(hori_index)) DEALLOCATE(hori_index)
2210	IF (ALLOCATED(horipft_index)) DEALLOCATE(horipft_index)
2211	IF (ALLOCATED(clay_fm)) DEALLOCATE(clay_fm)
2212	IF (ALLOCATED(humrel_daily_fm)) DEALLOCATE(humrel_daily_fm)
2213	IF (ALLOCATED(litterhum_daily_fm)) DEALLOCATE(litterhum_daily_fm)
2214	IF (ALLOCATED(t2m_daily_fm)) DEALLOCATE(t2m_daily_fm)
2215	IF (ALLOCATED(t2m_min_daily_fm)) DEALLOCATE(t2m_min_daily_fm)
2216	IF (ALLOCATED(tsurf_daily_fm)) DEALLOCATE(tsurf_daily_fm)
2217	IF (ALLOCATED(tsoil_daily_fm)) DEALLOCATE(tsoil_daily_fm)
2218	IF (ALLOCATED(soilhum_daily_fm)) DEALLOCATE(soilhum_daily_fm)
2219	IF (ALLOCATED(precip_fm)) DEALLOCATE(precip_fm)
2220	IF (ALLOCATED(gpp_daily_fm)) DEALLOCATE(gpp_daily_fm)
2221	IF (ALLOCATED(veget_fm)) DEALLOCATE(veget_fm)
2222	IF (ALLOCATED(veget_max_fm)) DEALLOCATE(veget_max_fm)
2223	IF (ALLOCATED(lai_fm)) DEALLOCATE(lai_fm)
2224
2225	IF (is_root_prc) THEN
2226	IF (ALLOCATED(clay_fm_g)) DEALLOCATE(clay_fm_g)
2227	IF (ALLOCATED(humrel_daily_fm_g)) DEALLOCATE(humrel_daily_fm_g)
2228	IF (ALLOCATED(litterhum_daily_fm_g)) DEALLOCATE(litterhum_daily_fm_g)
2229	IF (ALLOCATED(t2m_daily_fm_g)) DEALLOCATE(t2m_daily_fm_g)
2230	IF (ALLOCATED(t2m_min_daily_fm_g)) DEALLOCATE(t2m_min_daily_fm_g)
2231	IF (ALLOCATED(tsurf_daily_fm_g)) DEALLOCATE(tsurf_daily_fm_g)
2232	IF (ALLOCATED(tsoil_daily_fm_g)) DEALLOCATE(tsoil_daily_fm_g)
2233	IF (ALLOCATED(soilhum_daily_fm_g)) DEALLOCATE(soilhum_daily_fm_g)
2234	IF (ALLOCATED(precip_fm_g)) DEALLOCATE(precip_fm_g)
2235	IF (ALLOCATED(gpp_daily_fm_g)) DEALLOCATE(gpp_daily_fm_g)
2236	IF (ALLOCATED(veget_fm_g)) DEALLOCATE(veget_fm_g)
2237	IF (ALLOCATED(veget_max_fm_g)) DEALLOCATE(veget_max_fm_g)
2238	IF (ALLOCATED(lai_fm_g)) DEALLOCATE(lai_fm_g)
2239	ENDIF
2240
2241	IF (ALLOCATED(isf)) DEALLOCATE(isf)
2242	IF (ALLOCATED(nf_written)) DEALLOCATE(nf_written)
2243	IF (ALLOCATED(nf_cumul)) DEALLOCATE(nf_cumul)
2244	IF (ALLOCATED(nforce)) DEALLOCATE(nforce)
2245	IF (ALLOCATED(control_moist)) DEALLOCATE(control_moist)
2246	IF (ALLOCATED(control_temp)) DEALLOCATE(control_temp)
2247	IF (ALLOCATED(soilcarbon_input)) DEALLOCATE(soilcarbon_input)
2248	IF ( ALLOCATED (horip10_index)) DEALLOCATE (horip10_index)
2249	IF ( ALLOCATED (horip100_index)) DEALLOCATE (horip100_index)
2250	IF ( ALLOCATED (horip11_index)) DEALLOCATE (horip11_index)
2251	IF ( ALLOCATED (horip101_index)) DEALLOCATE (horip101_index)
2252	IF ( ALLOCATED (prod10)) DEALLOCATE (prod10)
2253	IF ( ALLOCATED (prod100)) DEALLOCATE (prod100)
2254	IF ( ALLOCATED (flux10)) DEALLOCATE (flux10)
2255	IF ( ALLOCATED (flux100)) DEALLOCATE (flux100)
2256	IF ( ALLOCATED (convflux)) DEALLOCATE (convflux)
2257	IF ( ALLOCATED (cflux_prod10)) DEALLOCATE (cflux_prod10)
2258	IF ( ALLOCATED (cflux_prod100)) DEALLOCATE (cflux_prod100)
2259	IF ( ALLOCATED (harvest_above)) DEALLOCATE (harvest_above)
2260	IF ( ALLOCATED (soilcarbon_input_daily)) DEALLOCATE (soilcarbon_input_daily)
2261	IF ( ALLOCATED (control_temp_daily)) DEALLOCATE (control_temp_daily)
2262	IF ( ALLOCATED (control_moist_daily)) DEALLOCATE (control_moist_daily)
2263
2264	IF ( ALLOCATED (fpc_max)) DEALLOCATE (fpc_max)
2265
2266	! 2. reset l_first
2267	l_first_stomate=.TRUE.
2268	! 3. call to clear functions
2269	CALL get_reftemp_clear
2270	CALL season_clear
2271	CALL stomatelpj_clear
2272	CALL littercalc_clear
2273	CALL vmax_clear
2274	!---------------------------
2275	END SUBROUTINE stomate_clear
2276	!
2277	!=
2278	!
2279	SUBROUTINE stomate_var_init &
2280	& (kjpindex, veget_cov, veget_cov_max, leaf_age, leaf_frac, &
2281	& tlong_ref, t2m_month, dead_leaves, &
2282	& veget, lai, qsintmax, deadleaf_cover, assim_param)
2283
2284	!---------------------------------------------------------------------
2285	! this subroutine outputs values of assim_param etc.
2286	! only if ok_stomate = .TRUE.
2287	! otherwise,the calling procedure must do it itself.
2288	!
2289	! interface description
2290	! input scalar
2291	! Domain size
2292	INTEGER(i_std),INTENT(in) :: kjpindex
2293	! input fields
2294	! fractional coverage: actually covered space
2295	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in) :: veget_cov
2296	! fractional coverage: maximum covered space
2297	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in) :: veget_cov_max
2298	! "long term" 2 meter reference temperatures (K)
2299	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: tlong_ref
2300	! "monthly" 2 meter temperatures (K)
2301	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: t2m_month
2302	! dead leaves on ground, per PFT, metabolic and structural,
2303	! in gC/(m**2 of ground)
2304	REAL(r_std),DIMENSION(kjpindex,nvm,nlitt),INTENT(in) :: dead_leaves
2305	! Fraction of vegetation type
2306	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in) :: veget
2307	! Surface foliere
2308	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in) :: lai
2309	! modified fields (actually NOT modified)
2310	! leaf age (d)
2311	REAL(r_std),DIMENSION(kjpindex,nvm,nleafages),INTENT(inout) :: leaf_age
2312	! fraction of leaves in leaf age class
2313	REAL(r_std),DIMENSION(kjpindex,nvm,nleafages),INTENT(inout) :: leaf_frac
2314	! output scalar
2315	! output fields
2316	! Maximum water on vegetation for interception
2317	REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(out) :: qsintmax
2318	! fraction of soil covered by dead leaves
2319	REAL(r_std),DIMENSION(kjpindex), INTENT (out) :: deadleaf_cover
2320	! min+max+opt temps & vmax for photosynthesis
2321	REAL(r_std),DIMENSION(kjpindex,nvm,npco2),INTENT(out) :: assim_param
2322
2323	!-
2324	! local declaration
2325	!-
2326	! dummy time step, must be zero
2327	REAL(r_std),PARAMETER :: dt_0 = zero
2328	REAL(r_std),DIMENSION(kjpindex,nvm) :: vcmax
2329	REAL(r_std),DIMENSION(kjpindex,nvm) :: vjmax
2330	! Min temperature for photosynthesis (deg C)
2331	REAL(r_std),DIMENSION(kjpindex,nvm) :: t_photo_min
2332	! Opt temperature for photosynthesis (deg C)
2333	REAL(r_std),DIMENSION(kjpindex,nvm) :: t_photo_opt
2334	! Max temperature for photosynthesis (deg C)
2335	REAL(r_std),DIMENSION(kjpindex,nvm) :: t_photo_max
2336	! Index
2337	INTEGER(i_std) :: j
2338	!---------------------------------------------------------------------
2339
2340	IF (control%ok_stomate) THEN
2341	!
2342	! 1 photosynthesis parameters
2343	!
2344	!
2345	! 1.1 vcmax
2346	! only if STOMATE is activated
2347	!
2348	CALL vmax (kjpindex, dt_0, leaf_age, leaf_frac, vcmax, vjmax)
2349	!
2350	! 1.2 assimilation temperatures
2351	!
2352	CALL assim_temp(kjpindex, tlong_ref, t2m_month, &
2353	& t_photo_min, t_photo_opt, t_photo_max)
2354	!
2355	! 1.3 transform into nvm vegetation types
2356	!
2357	assim_param(:,:,ivcmax) = zero
2358
2359	DO j = 2, nvm
2360	assim_param(:,j,ivcmax)=vcmax(:,j)
2361	ENDDO
2362
2363	assim_param(:,:,ivjmax) = zero
2364
2365	DO j = 2, nvm
2366	assim_param(:,j,ivjmax)=vjmax(:,j)
2367	ENDDO
2368
2369	assim_param(:,:,itmin) = zero
2370
2371	DO j = 2, nvm
2372	assim_param(:,j,itmin)=t_photo_min(:,j)
2373	ENDDO
2374
2375	assim_param(:,:,itopt) = zero
2376
2377	DO j = 2, nvm
2378	assim_param(:,j,itopt)=t_photo_opt(:,j)
2379	ENDDO
2380
2381	assim_param(:,:,itmax) = zero
2382
2383	DO j = 2, nvm
2384	assim_param(:,j,itmax)=t_photo_max(:,j)
2385	ENDDO
2386
2387	!
2388	! 2 dead leaf cover
2389	!
2390	! edit shilong
2391	! CALL deadleaf (kjpindex, dead_leaves, deadleaf_cover)
2392	CALL deadleaf (kjpindex, veget_cov_max, dead_leaves, deadleaf_cover)
2393	!
2394	! 3 qsintmax
2395	!
2396	qsintmax(:,ibare_sechiba) = zero
2397	DO j=2,nvm
2398	qsintmax(:,j) = qsintcstveget(:,j)lai(:,j)
2399	ENDDO
2400
2401	ENDIF ! ok_stomate = .TRUE.
2402	!--------------------------------
2403	END SUBROUTINE stomate_var_init
2404	!
2405	!=
2406	!
2407	SUBROUTINE stomate_accu &
2408	& (npts, n_dim2, dt_tot, dt, ldmean, field_in, field_out)
2409	!---------------------------------------------------------------------
2410	!
2411	! 0 declarations
2412	!
2413	! 0.1 input
2414	!
2415	! Domain size
2416	INTEGER(i_std),INTENT(in) :: npts
2417	! 2nd dimension (1 or nvm)
2418	INTEGER(i_std),INTENT(in) :: n_dim2
2419	! Time step of STOMATE (days)
2420	REAL(r_std),INTENT(in) :: dt_tot
2421	! Time step in days
2422	REAL(r_std),INTENT(in) :: dt
2423	! Calculate mean ?
2424	LOGICAL,INTENT(in) :: ldmean
2425	! Daily field
2426	REAL(r_std),DIMENSION(npts,n_dim2),INTENT(in) :: field_in
2427	!
2428	! 0.2 modified field
2429	!
2430	! Annual field
2431	REAL(r_std),DIMENSION(npts,n_dim2),INTENT(inout) :: field_out
2432	!---------------------------------------------------------------------
2433	!
2434	! 1 accumulation
2435	!
2436	field_out(:,:) = field_out(:,:)+field_in(:,:)*dt
2437	!
2438	! 2 mean fields
2439	!
2440	IF (ldmean) THEN
2441	field_out(:,:) = field_out(:,:)/dt_tot
2442	ENDIF
2443	!---------------------------------------------------------------------
2444	END SUBROUTINE stomate_accu
2445
2446	SUBROUTINE init_forcing (kjpindex,nsfm,nsft)
2447	!---------------------------------------------------------------------
2448	INTEGER(i_std),INTENT(in) :: kjpindex
2449	INTEGER(i_std),INTENT(in) :: nsfm
2450	INTEGER(i_std),INTENT(in) :: nsft
2451	!
2452	LOGICAL :: l_error
2453	INTEGER(i_std) :: ier
2454	!---------------------------------------------------------------------
2455	l_error = .FALSE.
2456	!
2457	ALLOCATE(clay_fm(kjpindex,nsfm),stat=ier)
2458	l_error = l_error .OR. (ier /= 0)
2459	ALLOCATE(humrel_daily_fm(kjpindex,nvm,nsfm),stat=ier)
2460	l_error = l_error .OR. (ier /= 0)
2461	ALLOCATE(litterhum_daily_fm(kjpindex,nsfm),stat=ier)
2462	l_error = l_error .OR. (ier /= 0)
2463	ALLOCATE(t2m_daily_fm(kjpindex,nsfm),stat=ier)
2464	l_error = l_error .OR. (ier /= 0)
2465	ALLOCATE(t2m_min_daily_fm(kjpindex,nsfm),stat=ier)
2466	l_error = l_error .OR. (ier /= 0)
2467	ALLOCATE(tsurf_daily_fm(kjpindex,nsfm),stat=ier)
2468	l_error = l_error .OR. (ier /= 0)
2469	ALLOCATE(tsoil_daily_fm(kjpindex,nbdl,nsfm),stat=ier)
2470	l_error = l_error .OR. (ier /= 0)
2471	ALLOCATE(soilhum_daily_fm(kjpindex,nbdl,nsfm),stat=ier)
2472	l_error = l_error .OR. (ier /= 0)
2473	ALLOCATE(precip_fm(kjpindex,nsfm),stat=ier)
2474	l_error = l_error .OR. (ier /= 0)
2475	ALLOCATE(gpp_daily_fm(kjpindex,nvm,nsfm),stat=ier)
2476	l_error = l_error .OR. (ier /= 0)
2477	ALLOCATE(veget_fm(kjpindex,nvm,nsfm),stat=ier)
2478	l_error = l_error .OR. (ier /= 0)
2479	ALLOCATE(veget_max_fm(kjpindex,nvm,nsfm),stat=ier)
2480	l_error = l_error .OR. (ier /= 0)
2481	ALLOCATE(lai_fm(kjpindex,nvm,nsfm),stat=ier)
2482	l_error = l_error .OR. (ier /= 0)
2483	ALLOCATE(isf(nsfm),stat=ier)
2484	l_error = l_error .OR. (ier /= 0)
2485	ALLOCATE(nf_written(nsft),stat=ier)
2486	l_error = l_error .OR. (ier /= 0)
2487	ALLOCATE(nf_cumul(nsft),stat=ier)
2488	l_error = l_error .OR. (ier /= 0)
2489	IF (l_error) THEN
2490	WRITE(numout,*) 'Problem with memory allocation: forcing variables'
2491	STOP 'init_forcing'
2492	ENDIF
2493
2494	IF (is_root_prc) THEN
2495	ALLOCATE(clay_fm_g(nbp_glo,nsfm),stat=ier)
2496	l_error = l_error .OR. (ier /= 0)
2497	ALLOCATE(humrel_daily_fm_g(nbp_glo,nvm,nsfm),stat=ier)
2498	l_error = l_error .OR. (ier /= 0)
2499	ALLOCATE(litterhum_daily_fm_g(nbp_glo,nsfm),stat=ier)
2500	l_error = l_error .OR. (ier /= 0)
2501	ALLOCATE(t2m_daily_fm_g(nbp_glo,nsfm),stat=ier)
2502	l_error = l_error .OR. (ier /= 0)
2503	ALLOCATE(t2m_min_daily_fm_g(nbp_glo,nsfm),stat=ier)
2504	l_error = l_error .OR. (ier /= 0)
2505	ALLOCATE(tsurf_daily_fm_g(nbp_glo,nsfm),stat=ier)
2506	l_error = l_error .OR. (ier /= 0)
2507	ALLOCATE(tsoil_daily_fm_g(nbp_glo,nbdl,nsfm),stat=ier)
2508	l_error = l_error .OR. (ier /= 0)
2509	ALLOCATE(soilhum_daily_fm_g(nbp_glo,nbdl,nsfm),stat=ier)
2510	l_error = l_error .OR. (ier /= 0)
2511	ALLOCATE(precip_fm_g(nbp_glo,nsfm),stat=ier)
2512	l_error = l_error .OR. (ier /= 0)
2513	ALLOCATE(gpp_daily_fm_g(nbp_glo,nvm,nsfm),stat=ier)
2514	l_error = l_error .OR. (ier /= 0)
2515	ALLOCATE(veget_fm_g(nbp_glo,nvm,nsfm),stat=ier)
2516	l_error = l_error .OR. (ier /= 0)
2517	ALLOCATE(veget_max_fm_g(nbp_glo,nvm,nsfm),stat=ier)
2518	l_error = l_error .OR. (ier /= 0)
2519	ALLOCATE(lai_fm_g(nbp_glo,nvm,nsfm),stat=ier)
2520	l_error = l_error .OR. (ier /= 0)
2521	IF (l_error) THEN
2522	WRITE(numout,*) 'Problem with memory allocation: forcing variables'
2523	STOP 'init_forcing'
2524	ENDIF
2525	ELSE
2526	ALLOCATE(clay_fm_g(0,nsfm),stat=ier)
2527	ALLOCATE(humrel_daily_fm_g(0,nvm,nsfm),stat=ier)
2528	ALLOCATE(litterhum_daily_fm_g(0,nsfm),stat=ier)
2529	ALLOCATE(t2m_daily_fm_g(0,nsfm),stat=ier)
2530	ALLOCATE(t2m_min_daily_fm_g(0,nsfm),stat=ier)
2531	ALLOCATE(tsurf_daily_fm_g(0,nsfm),stat=ier)
2532	ALLOCATE(tsoil_daily_fm_g(0,nbdl,nsfm),stat=ier)
2533	ALLOCATE(soilhum_daily_fm_g(0,nbdl,nsfm),stat=ier)
2534	ALLOCATE(precip_fm_g(0,nsfm),stat=ier)
2535	ALLOCATE(gpp_daily_fm_g(0,nvm,nsfm),stat=ier)
2536	ALLOCATE(veget_fm_g(0,nvm,nsfm),stat=ier)
2537	ALLOCATE(veget_max_fm_g(0,nvm,nsfm),stat=ier)
2538	ALLOCATE(lai_fm_g(0,nvm,nsfm),stat=ier)
2539	ENDIF
2540	!
2541	IF (l_error) THEN
2542	WRITE(numout,*) 'Problem with memory allocation: forcing variables'
2543	STOP 'init_forcing'
2544	ENDIF
2545	!
2546	CALL forcing_zero
2547	!--------------------------
2548	END SUBROUTINE init_forcing
2549	!
2550	!=
2551	!
2552	SUBROUTINE forcing_zero
2553	!---------------------------------------------------------------------
2554	clay_fm(:,:) = zero
2555	humrel_daily_fm(:,:,:) = zero
2556	litterhum_daily_fm(:,:) = zero
2557	t2m_daily_fm(:,:) = zero
2558	t2m_min_daily_fm(:,:) = zero
2559	tsurf_daily_fm(:,:) = zero
2560	tsoil_daily_fm(:,:,:) = zero
2561	soilhum_daily_fm(:,:,:) = zero
2562	precip_fm(:,:) = zero
2563	gpp_daily_fm(:,:,:) = zero
2564	veget_fm(:,:,:) = zero
2565	veget_max_fm(:,:,:) = zero
2566	lai_fm(:,:,:) = zero
2567	!--------------------------------
2568	END SUBROUTINE forcing_zero
2569	!
2570	!=
2571	!
2572	SUBROUTINE forcing_write(forcing_id,ibeg,iend)
2573	!---------------------------------------------------------------------
2574	INTEGER(i_std),INTENT(in) :: forcing_id
2575	INTEGER(i_std),INTENT(in) :: ibeg, iend
2576	!
2577	INTEGER(i_std) :: iisf, iblocks, nblocks
2578	INTEGER(i_std) :: ier
2579	INTEGER(i_std),DIMENSION(0:2) :: ifirst, ilast
2580	INTEGER(i_std),PARAMETER :: ndm = 10
2581	INTEGER(i_std),DIMENSION(ndm) :: start, count_force
2582	INTEGER(i_std) :: ndim, vid
2583	!---------------------------------------------------------------------
2584	!
2585	! determine blocks of forcing states that are contiguous in memory
2586	!
2587	nblocks = 0
2588	ifirst(:) = 1
2589	ilast(:) = 1
2590	!
2591	DO iisf = ibeg, iend
2592	IF ( (nblocks /= 0) &
2593	& .AND.(isf(iisf) == isf(ilast(nblocks))+1)) THEN
2594	! element is contiguous with last element found
2595	ilast(nblocks) = iisf
2596	ELSE
2597	! found first element of new block
2598	nblocks = nblocks+1
2599	IF (nblocks > 2) STOP 'Problem in forcing_write'
2600	ifirst(nblocks) = iisf
2601	ilast(nblocks) = iisf
2602	ENDIF
2603	ENDDO
2604	!
2605	CALL gather(clay_fm,clay_fm_g)
2606	CALL gather(humrel_daily_fm,humrel_daily_fm_g)
2607	CALL gather(litterhum_daily_fm,litterhum_daily_fm_g)
2608	CALL gather(t2m_daily_fm,t2m_daily_fm_g)
2609	CALL gather(t2m_min_daily_fm,t2m_min_daily_fm_g)
2610	CALL gather(tsurf_daily_fm,tsurf_daily_fm_g)
2611	CALL gather(tsoil_daily_fm,tsoil_daily_fm_g)
2612	CALL gather(soilhum_daily_fm,soilhum_daily_fm_g)
2613	CALL gather(precip_fm,precip_fm_g)
2614	CALL gather(gpp_daily_fm,gpp_daily_fm_g)
2615	CALL gather(veget_fm,veget_fm_g)
2616	CALL gather(veget_max_fm,veget_max_fm_g)
2617	CALL gather(lai_fm,lai_fm_g)
2618	IF (is_root_prc) THEN
2619	DO iblocks = 1, nblocks
2620	IF (ifirst(iblocks) /= ilast(iblocks)) THEN
2621	ndim = 2
2622	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2623	count_force(1:ndim) = SHAPE(clay_fm_g)
2624	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2625	ier = NF90_INQ_VARID (forcing_id,'clay',vid)
2626	ier = NF90_PUT_VAR (forcing_id,vid, &
2627	& clay_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
2628	& start=start(1:ndim), count=count_force(1:ndim))
2629	ndim = 3;
2630	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2631	count_force(1:ndim) = SHAPE(humrel_daily_fm_g)
2632	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2633	ier = NF90_INQ_VARID (forcing_id,'humrel',vid)
2634	ier = NF90_PUT_VAR (forcing_id, vid, &
2635	& humrel_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2636	& start=start(1:ndim), count=count_force(1:ndim))
2637	ndim = 2;
2638	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2639	count_force(1:ndim) = SHAPE(litterhum_daily_fm_g)
2640	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2641	ier = NF90_INQ_VARID (forcing_id,'litterhum',vid)
2642	ier = NF90_PUT_VAR (forcing_id, vid, &
2643	& litterhum_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
2644	& start=start(1:ndim), count=count_force(1:ndim))
2645	ndim = 2;
2646	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2647	count_force(1:ndim) = SHAPE(t2m_daily_fm_g)
2648	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2649	ier = NF90_INQ_VARID (forcing_id,'t2m',vid)
2650	ier = NF90_PUT_VAR (forcing_id, vid, &
2651	& t2m_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
2652	& start=start(1:ndim), count=count_force(1:ndim))
2653	ndim = 2;
2654	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2655	count_force(1:ndim) = SHAPE(t2m_min_daily_fm_g)
2656	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2657	ier = NF90_INQ_VARID (forcing_id,'t2m_min',vid)
2658	ier = NF90_PUT_VAR (forcing_id, vid, &
2659	& t2m_min_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
2660	& start=start(1:ndim), count=count_force(1:ndim))
2661	ndim = 2;
2662	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2663	count_force(1:ndim) = SHAPE(tsurf_daily_fm_g)
2664	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2665	ier = NF90_INQ_VARID (forcing_id,'tsurf',vid)
2666	ier = NF90_PUT_VAR (forcing_id, vid, &
2667	& tsurf_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
2668	& start=start(1:ndim), count=count_force(1:ndim))
2669	ndim = 3;
2670	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2671	count_force(1:ndim) = SHAPE(tsoil_daily_fm_g)
2672	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2673	ier = NF90_INQ_VARID (forcing_id,'tsoil',vid)
2674	ier = NF90_PUT_VAR (forcing_id, vid, &
2675	& tsoil_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2676	& start=start(1:ndim), count=count_force(1:ndim))
2677	ndim = 3;
2678	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2679	count_force(1:ndim) = SHAPE(soilhum_daily_fm_g)
2680	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2681	ier = NF90_INQ_VARID (forcing_id,'soilhum',vid)
2682	ier = NF90_PUT_VAR (forcing_id, vid, &
2683	& soilhum_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2684	& start=start(1:ndim), count=count_force(1:ndim))
2685	ndim = 2;
2686	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2687	count_force(1:ndim) = SHAPE(precip_fm_g)
2688	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2689	ier = NF90_INQ_VARID (forcing_id,'precip',vid)
2690	ier = NF90_PUT_VAR (forcing_id, vid, &
2691	& precip_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
2692	& start=start(1:ndim), count=count_force(1:ndim))
2693	ndim = 3;
2694	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2695	count_force(1:ndim) = SHAPE(gpp_daily_fm_g)
2696	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2697	ier = NF90_INQ_VARID (forcing_id,'gpp',vid)
2698	ier = NF90_PUT_VAR (forcing_id, vid, &
2699	& gpp_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2700	& start=start(1:ndim), count=count_force(1:ndim))
2701	ndim = 3;
2702	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2703	count_force(1:ndim) = SHAPE(veget_fm_g)
2704	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2705	ier = NF90_INQ_VARID (forcing_id,'veget',vid)
2706	ier = NF90_PUT_VAR (forcing_id, vid, &
2707	& veget_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2708	& start=start(1:ndim), count=count_force(1:ndim))
2709	ndim = 3;
2710	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2711	count_force(1:ndim) = SHAPE(veget_max_fm_g)
2712	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2713	ier = NF90_INQ_VARID (forcing_id,'veget_max',vid)
2714	ier = NF90_PUT_VAR (forcing_id, vid, &
2715	& veget_max_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2716	& start=start(1:ndim), count=count_force(1:ndim))
2717	ndim = 3;
2718	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2719	count_force(1:ndim) = SHAPE(lai_fm_g)
2720	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2721	ier = NF90_INQ_VARID (forcing_id,'lai',vid)
2722	ier = NF90_PUT_VAR (forcing_id, vid, &
2723	& lai_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2724	& start=start(1:ndim), count=count_force(1:ndim))
2725	ENDIF
2726	ENDDO
2727	ENDIF
2728	nf_written(isf(:)) = .TRUE.
2729	!---------------------------
2730	END SUBROUTINE forcing_write
2731	!
2732	!=
2733	!
2734	SUBROUTINE forcing_read(forcing_id,nsfm)
2735	!---------------------------------------------------------------------
2736	INTEGER(i_std),INTENT(in) :: forcing_id
2737	INTEGER(i_std),INTENT(in) :: nsfm
2738	!
2739	INTEGER(i_std) :: iisf, iblocks, nblocks
2740	INTEGER(i_std) :: ier
2741	LOGICAL :: a_er
2742	INTEGER(i_std),DIMENSION(0:2) :: ifirst, ilast
2743	INTEGER(i_std),PARAMETER :: ndm = 10
2744	INTEGER(i_std),DIMENSION(ndm) :: start, count_force
2745	INTEGER(i_std) :: ndim, vid
2746
2747	LOGICAL, PARAMETER :: check=.FALSE.
2748
2749	IF (check) WRITE(numout,*) "forcing_read "
2750	!---------------------------------------------------------------------
2751	!
2752	! set to zero if the corresponding forcing state
2753	! has not yet been written into the file
2754	!
2755	DO iisf = 1, nsfm
2756	IF (.NOT.nf_written(isf(iisf))) THEN
2757	clay_fm(:,iisf) = zero
2758	humrel_daily_fm(:,:,iisf) = zero
2759	litterhum_daily_fm(:,iisf) = zero
2760	t2m_daily_fm(:,iisf) = zero
2761	t2m_min_daily_fm(:,iisf) = zero
2762	tsurf_daily_fm(:,iisf) = zero
2763	tsoil_daily_fm(:,:,iisf) = zero
2764	soilhum_daily_fm(:,:,iisf) = zero
2765	precip_fm(:,iisf) = zero
2766	gpp_daily_fm(:,:,iisf) = zero
2767	veget_fm(:,:,iisf) = zero
2768	veget_max_fm(:,:,iisf) = zero
2769	lai_fm(:,:,iisf) = zero
2770	ENDIF
2771	ENDDO
2772	!
2773	! determine blocks of forcing states that are contiguous in memory
2774	!
2775	nblocks = 0
2776	ifirst(:) = 1
2777	ilast(:) = 1
2778	!
2779	DO iisf = 1, nsfm
2780	IF (nf_written(isf(iisf))) THEN
2781	IF ( (nblocks /= 0) &
2782	& .AND.(isf(iisf) == isf(ilast(nblocks))+1)) THEN
2783	! element is contiguous with last element found
2784	ilast(nblocks) = iisf
2785	ELSE
2786	! found first element of new block
2787	nblocks = nblocks+1
2788	IF (nblocks > 2) STOP 'Problem in forcing_read'
2789	!
2790	ifirst(nblocks) = iisf
2791	ilast(nblocks) = iisf
2792	ENDIF
2793	ENDIF
2794	ENDDO
2795	IF (check) WRITE(numout,*) "forcing_read nblocks, ifirst, ilast",nblocks, ifirst, ilast
2796	!
2797	IF (is_root_prc) THEN
2798	DO iblocks = 1, nblocks
2799	IF (check) WRITE(numout,*) "forcing_read iblocks, ifirst(iblocks), ilast(iblocks)",iblocks, &
2800	ifirst(iblocks), ilast(iblocks)
2801	IF (ifirst(iblocks) /= ilast(iblocks)) THEN
2802	a_er=.FALSE.
2803	ndim = 2;
2804	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2805	count_force(1:ndim) = SHAPE(clay_fm_g)
2806	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2807	ier = NF90_INQ_VARID (forcing_id,'clay',vid)
2808	a_er = a_er.OR.(ier.NE.0)
2809	ier = NF90_GET_VAR (forcing_id, vid, &
2810	& clay_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
2811	& start=start(1:ndim), count=count_force(1:ndim))
2812	a_er = a_er.OR.(ier.NE.0)
2813	!---------
2814	ndim = 3;
2815	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2816	count_force(1:ndim) = SHAPE(humrel_daily_fm_g)
2817	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2818	ier = NF90_INQ_VARID (forcing_id,'humrel',vid)
2819	a_er = a_er.OR.(ier.NE.0)
2820	ier = NF90_GET_VAR (forcing_id, vid, &
2821	& humrel_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2822	& start=start(1:ndim), count=count_force(1:ndim))
2823	a_er = a_er.OR.(ier.NE.0)
2824	!---------
2825	ndim = 2;
2826	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2827	count_force(1:ndim) = SHAPE(litterhum_daily_fm_g)
2828	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2829	ier = NF90_INQ_VARID (forcing_id,'litterhum',vid)
2830	a_er = a_er.OR.(ier.NE.0)
2831	ier = NF90_GET_VAR (forcing_id, vid, &
2832	& litterhum_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
2833	& start=start(1:ndim), count=count_force(1:ndim))
2834	a_er = a_er.OR.(ier.NE.0)
2835	!---------
2836	ndim = 2;
2837	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2838	count_force(1:ndim) = SHAPE(t2m_daily_fm_g)
2839	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2840	ier = NF90_INQ_VARID (forcing_id,'t2m',vid)
2841	a_er = a_er.OR.(ier.NE.0)
2842	ier = NF90_GET_VAR (forcing_id, vid, &
2843	& t2m_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
2844	& start=start(1:ndim), count=count_force(1:ndim))
2845	a_er = a_er.OR.(ier.NE.0)
2846	!---------
2847	ndim = 2;
2848	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2849	count_force(1:ndim) = SHAPE(t2m_min_daily_fm_g)
2850	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2851	ier = NF90_INQ_VARID (forcing_id,'t2m_min',vid)
2852	a_er = a_er.OR.(ier.NE.0)
2853	ier = NF90_GET_VAR (forcing_id, vid, &
2854	& t2m_min_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
2855	& start=start(1:ndim), count=count_force(1:ndim))
2856	a_er = a_er.OR.(ier.NE.0)
2857	!---------
2858	ndim = 2;
2859	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2860	count_force(1:ndim) = SHAPE(tsurf_daily_fm_g)
2861	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2862	ier = NF90_INQ_VARID (forcing_id,'tsurf',vid)
2863	a_er = a_er.OR.(ier.NE.0)
2864	ier = NF90_GET_VAR (forcing_id, vid, &
2865	& tsurf_daily_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
2866	& start=start(1:ndim), count=count_force(1:ndim))
2867	a_er = a_er.OR.(ier.NE.0)
2868	!---------
2869	ndim = 3;
2870	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2871	count_force(1:ndim) = SHAPE(tsoil_daily_fm_g)
2872	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2873	ier = NF90_INQ_VARID (forcing_id,'tsoil',vid)
2874	a_er = a_er.OR.(ier.NE.0)
2875	ier = NF90_GET_VAR (forcing_id, vid, &
2876	& tsoil_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2877	& start=start(1:ndim), count=count_force(1:ndim))
2878	a_er = a_er.OR.(ier.NE.0)
2879	!---------
2880	ndim = 3;
2881	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2882	count_force(1:ndim) = SHAPE(soilhum_daily_fm_g)
2883	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2884	ier = NF90_INQ_VARID (forcing_id,'soilhum',vid)
2885	a_er = a_er.OR.(ier.NE.0)
2886	ier = NF90_GET_VAR (forcing_id, vid, &
2887	& soilhum_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2888	& start=start(1:ndim), count=count_force(1:ndim))
2889	a_er = a_er.OR.(ier.NE.0)
2890	!---------
2891	ndim = 2;
2892	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2893	count_force(1:ndim) = SHAPE(precip_fm_g)
2894	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2895	ier = NF90_INQ_VARID (forcing_id,'precip',vid)
2896	a_er = a_er.OR.(ier.NE.0)
2897	ier = NF90_GET_VAR (forcing_id, vid, &
2898	& precip_fm_g(:,ifirst(iblocks):ilast(iblocks)), &
2899	& start=start(1:ndim), count=count_force(1:ndim))
2900	a_er = a_er.OR.(ier.NE.0)
2901	!---------
2902	ndim = 3;
2903	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2904	count_force(1:ndim) = SHAPE(gpp_daily_fm_g)
2905	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2906	ier = NF90_INQ_VARID (forcing_id,'gpp',vid)
2907	a_er = a_er.OR.(ier.NE.0)
2908	ier = NF90_GET_VAR (forcing_id, vid, &
2909	& gpp_daily_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2910	& start=start(1:ndim), count=count_force(1:ndim))
2911	a_er = a_er.OR.(ier.NE.0)
2912	!---------
2913	ndim = 3;
2914	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2915	count_force(1:ndim) = SHAPE(veget_fm_g)
2916	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2917	ier = NF90_INQ_VARID (forcing_id,'veget',vid)
2918	a_er = a_er.OR.(ier.NE.0)
2919	ier = NF90_GET_VAR (forcing_id, vid, &
2920	& veget_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2921	& start=start(1:ndim), count=count_force(1:ndim))
2922	a_er = a_er.OR.(ier.NE.0)
2923	!---------
2924	ndim = 3;
2925	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2926	count_force(1:ndim) = SHAPE(veget_max_fm_g)
2927	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2928	ier = NF90_INQ_VARID (forcing_id,'veget_max',vid)
2929	a_er = a_er.OR.(ier.NE.0)
2930	ier = NF90_GET_VAR (forcing_id, vid, &
2931	& veget_max_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2932	& start=start(1:ndim), count=count_force(1:ndim))
2933	a_er = a_er.OR.(ier.NE.0)
2934	!---------
2935	ndim = 3;
2936	start(1:ndim) = 1; start(ndim) = isf(ifirst(iblocks));
2937	count_force(1:ndim) = SHAPE(lai_fm_g)
2938	count_force(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
2939	ier = NF90_INQ_VARID (forcing_id,'lai',vid)
2940	a_er = a_er.OR.(ier.NE.0)
2941	ier = NF90_GET_VAR (forcing_id, vid, &
2942	& lai_fm_g(:,:,ifirst(iblocks):ilast(iblocks)), &
2943	& start=start(1:ndim), count=count_force(1:ndim))
2944	a_er = a_er.OR.(ier.NE.0)
2945	IF (a_er) THEN
2946	CALL ipslerr (3,'forcing_read', &
2947	& 'PROBLEM when read forcing file', &
2948	& '','')
2949	ENDIF
2950	ENDIF
2951	ENDDO
2952	ENDIF
2953	CALL scatter(clay_fm_g,clay_fm)
2954	CALL scatter(humrel_daily_fm_g,humrel_daily_fm)
2955	CALL scatter(litterhum_daily_fm_g,litterhum_daily_fm)
2956	CALL scatter(t2m_daily_fm_g,t2m_daily_fm)
2957	CALL scatter(t2m_min_daily_fm_g,t2m_min_daily_fm)
2958	CALL scatter(tsurf_daily_fm_g,tsurf_daily_fm)
2959	CALL scatter(tsoil_daily_fm_g,tsoil_daily_fm)
2960	CALL scatter(soilhum_daily_fm_g,soilhum_daily_fm)
2961	CALL scatter(precip_fm_g,precip_fm)
2962	CALL scatter(gpp_daily_fm_g,gpp_daily_fm)
2963	CALL scatter(veget_fm_g,veget_fm)
2964	CALL scatter(veget_max_fm_g,veget_max_fm)
2965	CALL scatter(lai_fm_g,lai_fm)
2966	!--------------------------
2967	END SUBROUTINE forcing_read
2968	!
2969	!=
2970	!
2971	SUBROUTINE setlai(npts,lai)
2972	!---------------------------------------------------------------------
2973	! routine to force the lai in STOMATE (for assimilation procedures)
2974	! for the moment setlai only gives the lai from stomate,
2975	! this routine must be written in the future
2976	!
2977	! 0 declarations
2978	!
2979	! 0.1 input
2980	!
2981	! Domain size
2982	INTEGER(i_std),INTENT(in) :: npts
2983	! 0.3 output
2984	REAL(r_std),DIMENSION(npts,nvm),INTENT(out) :: lai
2985	! 0.4 local definitions
2986	INTEGER(i_std) :: j
2987	!---------------------------------------------------------------------
2988	lai(:,ibare_sechiba) = zero
2989	DO j=2,nvm
2990	lai(:,j) = biomass(:,j,ileaf)*sla(j)
2991	ENDDO
2992	!--------------------
2993	END SUBROUTINE setlai
2994	!
2995	!
2996	!-----------------
2997	END MODULE stomate

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

Download in other formats: