Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

stomate.f90 @ 4114

Last change on this file since 4114 was 119, checked in by martial.mancip, 13 years ago
PC, CE, MM : corrections for IPSLCM5A_C couple carbon model.
File size: 129.9 KB

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

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

Context Navigation

source: tags/ORCHIDEE_1_9_5_1/ORCHIDEE/src_stomate/stomate.f90 @ 4114

Download in other formats: