Context Navigation

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

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

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

Download in other formats: