Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE/src_stomate/stomate.f90 @ 64

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

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

Download in other formats: