Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

stomate.f90 @ 8

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

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

Context Navigation

source: tags/ORCHIDEE_1_9_5/ORCHIDEE/src_stomate/stomate.f90 @ 8

Download in other formats: