Context Navigation

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

Last change on this file since 277 was 277, checked in by didier.solyga, 13 years ago
Update the externalized version with the last commit of the trunk (revision 275)
File size: 130.8 KB

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

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

Download in other formats: