Context Navigation

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

Last change on this file since 257 was 257, checked in by didier.solyga, 13 years ago
Externalized version merged with the trunk
File size: 130.5 KB

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

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

Download in other formats: