Context Navigation

stomate_lpj.f90 @ 6737

Last change on this file since 6737 was 6538, checked in by chao.yue, 4 years ago
Rollback to r5301; previous commit does not achieve the right purpose
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 138.2 KB

Line
1	! ================================================================================================================================
2	! MODULE : stomate_lpj
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Main entry point for daily processes in STOMATE and LPJ (phenology,
10	!! allocation, npp_calc, kill, turn, light, establish, crown, cover, lcchange)
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S): None
15	!!
16	!! REFERENCE(S) : None
17	!!
18	!! SVN :
19	!! $HeadURL$
20	!! $Date$
21	!! $Revision$
22	!! \n
23	!_ ================================================================================================================================
24
25	MODULE stomate_lpj
26
27	! modules used:
28
29	USE ioipsl_para
30	USE xios_orchidee
31	USE grid
32	USE stomate_data
33	USE constantes
34	USE constantes_soil
35	USE pft_parameters
36	USE lpj_constraints
37	USE lpj_pftinout
38	USE lpj_kill
39	USE lpj_crown
40	USE lpj_fire
41	USE lpj_spitfire
42	USE lpj_gap
43	USE lpj_light
44	USE lpj_establish
45	USE lpj_cover
46	USE stomate_prescribe
47	USE stomate_phenology
48	USE stomate_alloc
49	USE stomate_npp
50	USE stomate_turnover
51	USE stomate_litter
52	USE stomate_soilcarbon
53	USE stomate_vmax
54	USE stomate_lcchange
55	USE stomate_gluc_common
56	USE stomate_glcchange_SinAgeC
57	USE stomate_glcchange_MulAgeC
58	USE stomate_fharvest_SinAgeC
59	USE stomate_fharvest_MulAgeC
60	USE stomate_glcc_bioe1
61	USE stomate_lai
62	USE stomate_gluc_constants
63	!gmjc
64	USE grassland_management
65	!end gmjc
66	USE stomate_check
67
68	IMPLICIT NONE
69
70	! private & public routines
71
72	PRIVATE
73	PUBLIC StomateLpj,StomateLpj_clear
74
75	LOGICAL, SAVE :: firstcall = .TRUE. !! first call
76	!$OMP THREADPRIVATE(firstcall)
77	!gmjc
78	! flag that enable grazing
79	LOGICAL, SAVE :: enable_grazing
80	!$OMP THREADPRIVATE(enable_grazing)
81	!end gmjc
82
83	! check mass balance for stomate
84	LOGICAL, SAVE :: STO_CHECK_MASSBALANCE = .TRUE.
85	!$OMP THREADPRIVATE(STO_CHECK_MASSBALANCE)
86
87	CONTAINS
88
89
90	!! ================================================================================================================================
91	!! SUBROUTINE : StomateLpj_clear
92	!!
93	!>\BRIEF Re-initialisation of variable
94	!!
95	!! DESCRIPTION : This subroutine reinitializes variables. To be used if we want to relaunch
96	!! ORCHIDEE but the routine is not used in current version.
97	!!
98	!! RECENT CHANGE(S) : None
99	!!
100	!! MAIN OUTPUT VARIABLE(S): None
101	!!
102	!! REFERENCE(S) : None
103	!!
104	!! FLOWCHART : None
105	!! \n
106	!_ ================================================================================================================================
107
108	SUBROUTINE StomateLpj_clear
109
110	CALL prescribe_clear
111	CALL phenology_clear
112	CALL npp_calc_clear
113	CALL turn_clear
114	CALL soilcarbon_clear
115	CALL constraints_clear
116	CALL establish_clear
117	CALL fire_clear
118	CALL spitfire_clear
119	CALL gap_clear
120	CALL light_clear
121	CALL pftinout_clear
122	CALL alloc_clear
123
124	CALL grassmanag_clear
125
126	END SUBROUTINE StomateLpj_clear
127
128
129	!! ================================================================================================================================
130	!! SUBROUTINE : StomateLPJ
131	!!
132	!>\BRIEF Main entry point for daily processes in STOMATE and LPJ, structures the call sequence
133	!! to the different processes such as dispersion, establishment, competition and mortality of PFT's.
134	!!
135	!! DESCRIPTION : This routine is the main entry point to all processes calculated on a
136	!! daily time step. Is mainly devoted to call the different STOMATE and LPJ routines
137	!! depending of the ok_dgvm (is dynamic veg used) and lpj_constant_mortality (is background mortality used).
138	!! It also prepares the cumulative
139	!! fluxes or pools (e.g TOTAL_M TOTAL_BM_LITTER etc...)
140	!!
141	!! This routine makes frequent use of "weekly", "monthly" and "long term" variables. Quotion is used because
142	!! by default "weekly" denotes 7 days, by default "monthly" denotes 20 days and by default "Long term" denotes
143	!! 3 years. dtslow refers to 24 hours (1 day).
144	!!
145	!!
146	!! RECENT CHANGE(S) : None
147	!!
148	!! MAIN OUTPUT VARIABLE(S): All variables related to stomate and required for LPJ dynamic vegetation mode.
149	!!
150	!! REFERENCE(S) :
151	!! - Krinner, G., N. Viovy, N. de Noblet-DucoudrÃ©, J. OgeÃ©, J. Polcher, P. Friedlingstein, P. Ciais, S. Sitch,
152	!! and I. C. Prentice. 2005. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere
153	!! system. Global Biogeochemical Cycles 19:GB1015, doi:1010.1029/2003GB002199.
154	!! - Sitch, S., B. Smith, I. C. Prentice, A. Arneth, A. Bondeau, W. Cramer, J. O. Kaplan, S. Levis, W. Lucht,
155	!! M. T. Sykes, K. Thonicke, and S. Venevsky. 2003. Evaluation of ecosystem dynamics, plant geography and
156	!! terrestrial carbon cycling in the LPJ dynamic global vegetation model. Global Change Biology 9:161-185.
157	!!
158	!! FLOWCHART : Update with existing flowchart from N Viovy (Jan 19, 2012)
159	!! \n
160	!_ ================================================================================================================================
161
162	SUBROUTINE StomateLpj (npts, dt_days, &
163	lalo, neighbours, resolution, contfrac, &
164	clay, herbivores, &
165	tsurf_daily, tsoil_daily, t2m_daily, t2m_min_daily, &
166	!spitfire
167	t2m_max_daily, precip_daily, wspeed_daily, lightn, popd, a_nd, &
168	read_observed_ba, observed_ba, &
169	read_cf_fine,cf_fine,read_cf_coarse,cf_coarse,read_ratio_flag,ratio_flag,read_ratio,ratio,date,&
170	!endspit
171	litterhum_daily, soilhum_daily, &
172	maxmoiavail_lastyear, minmoiavail_lastyear, &
173	gdd0_lastyear, precip_lastyear, &
174	moiavail_month, moiavail_week, t2m_longterm, t2m_month, t2m_week, &
175	tsoil_month, soilhum_month, &
176	gdd_m5_dormance, gdd_from_growthinit, gdd_midwinter, ncd_dormance, ngd_minus5, &
177	turnover_longterm, gpp_daily, gpp_week, &
178	time_hum_min, hum_min_dormance, maxfpc_lastyear, resp_maint_part, &
179	PFTpresent, age, fireindex, firelitter, &
180	leaf_age, leaf_frac, biomass, ind, adapted, regenerate, &
181	senescence, when_growthinit, &
182	litterpart, litter, litter_avail, litter_not_avail, litter_avail_frac, &
183	dead_leaves, carbon,carbon_surf, lignin_struc, &
184	!spitfire
185	ni_acc,fire_numday,fuel_1hr,fuel_10hr,fuel_100hr,fuel_1000hr, &
186	lcc,bafrac_deforest_accu,emideforest_litter_accu,emideforest_biomass_accu,&
187	deforest_litter_remain,deforest_biomass_remain,&
188	def_fuel_1hr_remain,def_fuel_10hr_remain,&
189	def_fuel_100hr_remain,def_fuel_1000hr_remain,&
190	!endspit
191	veget_cov_max, veget_cov_max_new, npp_longterm, lm_lastyearmax, &
192	veget_lastlight, everywhere, need_adjacent, RIP_time, &
193	lai, rprof,npp_daily, turnover_daily, turnover_time,&
194	control_moist, control_temp, soilcarbon_input, &
195	co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, &
196	height, deadleaf_cover, vcmax, &
197	bm_to_litter, &
198	prod10,prod100,flux10, flux100, &
199	vegetnew_firstday, glccNetLCC, &
200	glccSecondShift,glccPrimaryShift, &
201	harvest_matrix, harvest_biomass, bound_spa, newvegfrac, glcc_pft, &
202	convflux,cflux_prod10,cflux_prod100, &
203	harvest_above, carb_mass_total, &
204	fpc_max, Tseason, Tseason_length, Tseason_tmp, &
205	Tmin_spring_time, begin_leaves, onset_date, &
206	MatrixA,&
207	!!!! crop variables
208	pdlai, slai, &
209	! for crop allocation
210	in_cycle, deltai, dltaisen, ssla, pgrain, deltgrain, reprac, &
211	nger, nlev, ndrp, nlax, &
212	c_reserve, c_leafb, nmat, nrec, N_limfert, tday_counter, &
213	!!!! end crop variables, xuhui
214	zz_coef_deep, deepC_a, deepC_s, deepC_p, &
215	tsurf_year, & !pss:-
216	!gmjc
217	wshtotsum, sr_ugb, compt_ugb, nb_ani, grazed_frac, &
218	import_yield, sla_age1, t2m_14, sla_calc, snowfall_daily, day_of_year, &
219	when_growthinit_cut, nb_grazingdays, &
220	EndOfYear, &
221	moiavail_daily,tmc_topgrass_daily,fc_grazing,snowmass_daily, &
222	after_snow, after_wet, wet1day, wet2day)
223	!end gmjc
224
225	!! 0. Variable and parameter declaration
226
227	!! 0.1 input
228
229	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
230	REAL(r_std), INTENT(in) :: dt_days !! Time step of Stomate (days)
231	INTEGER(i_std), DIMENSION(npts,NbNeighb), INTENT(in) :: neighbours !! Indices of the 8 neighbours of each grid
232	!! point [1=North and then clockwise]
233	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: resolution !! Resolution at each grid point (m)
234	!! [1=E-W, 2=N-S]
235	REAL(r_std),DIMENSION(npts,2),INTENT(in) :: lalo !! Geographical coordinates (latitude,longitude)
236	!! for pixels (degrees)
237	REAL(r_std),DIMENSION (npts), INTENT (in) :: contfrac !! Fraction of continent in the grid cell (unitless)
238	REAL(r_std), DIMENSION(npts), INTENT(in) :: clay !! Clay fraction (0 to 1, unitless)
239	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: herbivores !! Time constant of probability of a leaf to
240	!! be eaten by a herbivore (days)
241	REAL(r_std), DIMENSION(npts), INTENT(in) :: tsurf_daily !! Daily surface temperatures (K)
242	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: tsoil_daily !! Daily soil temperatures (K)
243	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_daily !! Daily 2 meter temperatures (K)
244	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_min_daily !! Daily minimum 2 meter temperatures (K)
245	!spitfire
246	INTEGER(i_std),INTENT(in) :: date !! Date (days)
247	! daily maximum 2 meter temperatures (K)
248	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_max_daily
249	! daily precip (mm/day)
250	REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_daily
251	! Wind speed
252	REAL(r_std), DIMENSION(npts), INTENT(in) :: wspeed_daily
253	! Lightning flash rate
254	REAL(r_std), DIMENSION(npts), INTENT(in) :: lightn
255	! Population density rate
256	REAL(r_std), DIMENSION(npts), INTENT(inout) :: popd !popd declared and allocated and input in slowproc.f90
257	! Flag for read in observed burned area
258	LOGICAL, INTENT (in) :: read_observed_ba
259	! Observed burned area
260	REAL(r_std),DIMENSION (npts), INTENT (in) :: observed_ba
261
262	! Flag for read in observed burned area
263	LOGICAL, INTENT (in) :: read_cf_coarse
264	! Observed burned area
265	REAL(r_std),DIMENSION (npts), INTENT (in) :: cf_coarse
266	! Flag for read in observed burned area
267	LOGICAL, INTENT (in) :: read_cf_fine
268	! Observed burned area
269	REAL(r_std),DIMENSION (npts), INTENT (in) :: cf_fine
270	! Flag for read in observed burned area
271	LOGICAL, INTENT (in) :: read_ratio
272	! Observed burned area
273	REAL(r_std),DIMENSION (npts), INTENT (in) :: ratio
274	! Flag for read in observed burned area
275	LOGICAL, INTENT (in) :: read_ratio_flag
276	! Observed burned area
277	REAL(r_std),DIMENSION (npts), INTENT (in) :: ratio_flag
278
279	!endspit
280
281	REAL(r_std), DIMENSION(npts), INTENT(in) :: litterhum_daily !! Daily litter humidity (0 to 1, unitless)
282	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: soilhum_daily !! Daily soil humidity (0 to 1, unitless)
283	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: maxmoiavail_lastyear !! Last year's maximum moisture availability
284	!! (0 to 1, unitless)
285	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: minmoiavail_lastyear !! Last year's minimum moisture availability
286	!! (0 to 1, unitless)
287	REAL(r_std), DIMENSION(npts), INTENT(in) :: gdd0_lastyear !! Last year's GDD0 (K)
288	REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_lastyear !! Lastyear's precipitation
289	!! @tex $(mm year^{-1})$ @endtex
290	!! to determine if establishment possible
291	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
292	!! unitless)
293	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
294	!! (0 to 1, unitless)
295	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_longterm !! "Long term" 2 meter reference
296	!! temperatures (K)
297	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_month !! "Monthly" 2-meter temperatures (K)
298	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_week !! "Weekly" 2-meter temperatures (K)
299	! "seasonal" 2-meter temperatures (K)
300	REAL(r_std), DIMENSION(npts), INTENT(in) :: Tseason
301	! temporary variable to calculate Tseason
302	REAL(r_std), DIMENSION(npts), INTENT(in) :: Tseason_length
303	! temporary variable to calculate Tseason
304	REAL(r_std), DIMENSION(npts), INTENT(in) :: Tseason_tmp
305
306	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: Tmin_spring_time !! Number of days after begin_leaves (leaf onset)
307	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: onset_date !! Date in the year at when the leaves started to grow(begin_leaves)
308
309	!pss:+
310	REAL(r_std), DIMENSION(npts), INTENT(in) :: tsurf_year ! annual surface temperatures (K)
311	!pss:-
312	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: tsoil_month !! "Monthly" soil temperatures (K)
313	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: soilhum_month !! "Monthly" soil humidity
314	!! (0 to 1, unitless)
315	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
316	!! C (for phenology)
317	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit for crops
318	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
319	!! (for phenology) - this is written to the history files
320	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
321	!! leaves were lost (for phenology)
322	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
323	!! -5 deg C (for phenology)
324	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: turnover_longterm !! "Long term" turnover rate
325	!! @tex $(gC m^{-2} year^{-1})$ @endtex
326	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
327	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
328	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
329	!! @tex $(gC m^{-2} day^{-1})$ @endtex
330	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
331	!! availability (days)
332	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: hum_min_dormance !! minimum moisture during dormance
333	!! (0-1, unitless)
334	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: maxfpc_lastyear !! Last year's maximum foliage projected
335	!! coverage for each natural PFT,
336	!! @tex $(m^2 m^{-2})$ @endtex
337	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: resp_maint_part !! Maintenance respiration of different
338	!! plant parts
339	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
340	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: fpc_max !! "Maximal" coverage fraction of a PFT (LAI
341	!! -> infinity) on ground
342	!! @tex $(m^2 m^{-2})$ @endtex
343	REAL(r_std), DIMENSION(ndeep), INTENT (in) :: zz_coef_deep !! deep vertical profile
344	!!!!! crop variables
345	! FOR CROP---STICS
346	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: pdlai
347	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: slai
348
349	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: in_cycle
350	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: deltai
351	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: dltaisen
352	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ssla
353	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: pgrain
354	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: deltgrain
355	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: reprac
356	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nger
357	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nlev
358	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: ndrp
359	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nmat
360	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nlax
361
362	! REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: N_limfert !!
363	! defined already in GM
364	INTEGER(i_std), INTENT(in) :: tday_counter
365
366	!!!!! end crop variables, xuhui
367
368	!gmjc
369	LOGICAL, INTENT(in) :: EndOfYear
370	!end gmjc
371	!! 0.2 Output variables
372
373	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: npp_daily !! Net primary productivity
374	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
375	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out) :: turnover_daily !! Turnover rates
376	!! @tex $(gC m^{-2} year^{-1})$ @endtex
377	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: co2_to_bm !! CO2 taken up from atmosphere when
378	!! introducing a new PFT (introduced for
379	!! carbon balance closure)
380	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
381	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: co2_fire !! Carbon emitted into the atmosphere by
382	!! fire (living and dead biomass)
383	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
384	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
385	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
386	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_maint !! Maintenance respiration
387	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
388	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_growth !! Growth respiration
389	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
390
391	REAL(r_std), DIMENSION(npts), INTENT(inout) :: deadleaf_cover !! Fraction of soil covered by dead leaves
392	!! (0 to 1, unitless)
393	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: vcmax !! Maximum rate of carboxylation
394	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out):: bm_to_litter !! Conversion of biomass to litter
395	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
396	LOGICAL, DIMENSION(npts,nvm), INTENT(out) :: begin_leaves !! signal to start putting leaves on (true/false)
397
398	!!!!! crop variables
399	! for crop c pools
400	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: c_reserve
401	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: c_leafb
402
403	! for crop turnover
404	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nrec !harvest date
405	!!!!! end crop variables, xuhui
406
407	!! 0.3 Modified variables
408
409	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: height !! Height of vegetation (m)
410	REAL(r_std), DIMENSION(npts,nlevs), INTENT(inout) :: control_moist !! Moisture control of heterotrophic
411	!! respiration (0 to 1, unitless)
412	REAL(r_std), DIMENSION(npts,nlevs), INTENT(inout) :: control_temp !! Temperature control of heterotrophic
413	!! respiration, above and below
414	!! (0 to 1, unitless)
415	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: soilcarbon_input !! Quantity of carbon going into carbon
416	!! pools from litter decomposition
417	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
418	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lai !! Leaf area index OF AN INDIVIDUAL PLANT,
419	!! where a PFT contains n indentical plants
420	!! i.e., using the mean individual approach
421	!! @tex $(m^2 m^{-2})$ @endtex
422	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: rprof !! Prescribed root depth (m)
423	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
424	!! each pixel
425	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age !! Age (years)
426	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: fireindex !! Probability of fire (0 to 1, unitless)
427	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: firelitter !! Longer term litter above the ground that
428	!! can be burned, @tex $(gC m^{-2})$ @endtex
429	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_age !! Leaf age (days)
430	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac !! Fraction of leaves in leaf age class,
431	!! (0 to 1, unitless)
432	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass !! Biomass @tex $(gC m^{-2})$ @endtex
433	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind !! Density of individuals
434	!! @tex $(m^{-2})$ @endtex
435	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: adapted !! Adaptation of PFT (killed if too cold)
436	!! (0 to 1, unitless)
437	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: regenerate !! "Fitness": Winter sufficiently cold for
438	!! PFT regeneration ? (0 to 1, unitless)
439	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
440	!! for deciduous trees)
441	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
442	!! the growing season (days)
443	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: litterpart !! Fraction of litter above the ground
444	!! belonging to different PFTs
445	!! (0 to 1, unitless)
446	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout):: litter !! Metabolic and structural litter, above
447	!! and below ground
448	!! @tex $(gC m^{-2})$ @endtex
449	!gmjc for grazing litter
450	REAL(r_std), DIMENSION(npts,nlitt,nvm), INTENT(out):: litter_avail
451	REAL(r_std), DIMENSION(npts,nlitt,nvm) , INTENT(out):: litter_not_avail
452	REAL(r_std), DIMENSION(npts,nlitt,nvm), INTENT(in):: litter_avail_frac
453	!end gmjc
454	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: dead_leaves !! Dead leaves on ground, per PFT, metabolic
455	!! and structural,
456	!! @tex $(gC m^{-2})$ @endtex
457	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: carbon !! Carbon pool: active, slow, or passive,
458	!! @tex $(gC m^{-2})$ @endtex
459	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: carbon_surf
460	!! @tex $(gC m^{-2})$ @endtex
461	REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout) :: lignin_struc !! Ratio of Lignin/Carbon in structural
462	!! litter, above and below ground,
463	!! @tex $(gC m^{-2})$ @endtex
464	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_cov_max !! "Maximal" coverage fraction of a PFT (LAI
465	!! -> infinity) on ground
466	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_cov_max_new !! "Maximal" coverage fraction of a PFT
467	!! (LAI-> infinity) on ground (unitless)
468	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary
469	!! productivity
470	!! @tex $(gC m^{-2} year^{-1})$ @endtex
471	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lm_lastyearmax !! Last year's maximum leaf mass, for each
472	!! PFT @tex $(gC m^{-2})$ @endtex
473	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_lastlight !! Vegetation fractions (on ground) after
474	!! last light competition
475	!! @tex $(m^2 m^{-2})$ @endtex
476	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: everywhere !! Is the PFT everywhere in the grid box or
477	!! very localized (after its introduction)
478	!! (unitless)
479	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: need_adjacent !! In order for this PFT to be introduced,
480	!! does it have to be present in an
481	!! adjacent grid box?
482	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: RIP_time !! How much time ago was the PFT eliminated
483	!! for the last time (y)
484	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: turnover_time !! Turnover_time of leaves for grasses
485	!! (days)
486	REAL(r_std), DIMENSION(npts,nvm),INTENT(in) :: vegetnew_firstday !! New "maximal" coverage fraction of a PFT
487	!! (LAI -> infinity) (unitless)
488	REAL(r_std),DIMENSION(npts,12),INTENT(inout) :: glccNetLCC
489	REAL(r_std),DIMENSION(npts,12),INTENT(inout) :: glccSecondShift
490	REAL(r_std),DIMENSION(npts,12),INTENT(inout) :: glccPrimaryShift
491	REAL(r_std), DIMENSION(npts,12),INTENT(inout) :: harvest_matrix !! The gross land use change matrix in case
492	!! of gross land cover change is simulated.
493	REAL(r_std), DIMENSION(npts,12),INTENT(in) :: harvest_biomass !! The gross land use change matrix in case
494	!! of gross land cover change is simulated.
495	REAL(r_std), DIMENSION(npts,nvm),INTENT(inout) :: bound_spa !! The gross land use change matrix in case
496	!! of gross land cover change is simulated.
497	REAL(r_std), DIMENSION(npts,nvmap),INTENT(in) :: newvegfrac !! The gross land use change matrix in case
498	!! of gross land cover change is simulated.
499	REAL(r_std),DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! Products remaining in the 10
500	!! year-turnover pool after the annual
501	!! release for each compartment (10
502	!! + 1 : input from year of land cover
503	!! change) @tex $(gC m^{-2})$ @endtex
504	REAL(r_std),DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! Products remaining in the 100
505	!! year-turnover pool after the annual
506	!! release for each compartment (100
507	!! + 1 : input from year of land cover
508	!! change) @tex $(gC m^{-2})$ @endtex
509	REAL(r_std),DIMENSION(npts,10,nwp), INTENT(inout) :: flux10 !! Annual release from the 10
510	!! year-turnover pool compartments
511	!! @tex $(gC m^{-2})$ @endtex
512	REAL(r_std),DIMENSION(npts,100,nwp), INTENT(inout) :: flux100 !! Annual release from the 100
513	!! year-turnover pool compartments
514	!! @tex $(gC m^{-2})$ @endtex
515	REAL(r_std),DIMENSION(npts,nwp), INTENT(inout) :: convflux !! Release during first year following land
516	!! cover change @tex $(gC m^{-2})$ @endtex
517	REAL(r_std),DIMENSION(npts,nwp), INTENT(inout) :: cflux_prod10 !! Total annual release from the 10
518	!! year-turnover pool
519	!! @tex $(gC m^{-2})$ @endtex
520	REAL(r_std),DIMENSION(npts,nwp), INTENT(inout) :: cflux_prod100 !! Total annual release from the 100
521	!! year-turnover pool
522	!! @tex $(gC m^{-2})$ @endtex
523	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: harvest_above !! Harvest above ground biomass for
524	!! agriculture @tex $(gC m^{-2})$ @endtex
525	REAL(r_std), DIMENSION(npts), INTENT(inout) :: carb_mass_total !! Carbon Mass total (soil, litter, veg)
526	!! @tex $(gC m^{-2})$ @endtex
527	REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(inout) :: MatrixA !! Matrix containing the fluxes
528	!! between the carbon pools
529	!! per sechiba time step
530	!! @tex $(gC.m^2.day^{-1})$ @endtex
531	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_a !! permafrost soil carbon (g/m**3) active
532	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_s !! permafrost soil carbon (g/m**3) slow
533	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_p !! permafrost soil carbon (g/m**3) passive
534	!gmjc
535
536	!glcc
537
538	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! a temporary variable to hold the fraction of ipft->ivma, i.e., from
539	!! PFT_{ipft} to the youngest age class of MTC_{ivma}
540	!spitfire
541	! Nesterov index accumulated
542	REAL(r_std), DIMENSION(npts), INTENT(inout) :: ni_acc
543	REAL(r_std), DIMENSION(npts), INTENT(inout) :: fire_numday
544	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: bafrac_deforest_accu
545	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: emideforest_litter_accu
546	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: emideforest_biomass_accu
547	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) :: deforest_litter_remain !! Vegetmax-weighted remaining litter on the ground for
548	!! deforestation region.
549	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
550	!! deforestation region.
551	! fuel classes (1, 10, 100, 1000 hours)
552	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: fuel_1hr
553	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: fuel_10hr
554	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: fuel_100hr
555	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: fuel_1000hr
556	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: def_fuel_1hr_remain
557	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: def_fuel_10hr_remain
558	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: def_fuel_100hr_remain
559	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: def_fuel_1000hr_remain
560	REAL(r_std), DIMENSION(npts) :: d_area_burnt
561	REAL(r_std), DIMENSION(npts) :: d_numfire
562	REAL(r_std), DIMENSION(npts) :: fc_crown
563	! parameter for potential human-caused ignitions, ignitions ind^{-1}day{-1}, used in lpj_spitfire.f90
564	REAL(r_std), DIMENSION(npts), INTENT(in) :: a_nd
565	!endspit
566	!gmjc
567	! snowfall_daily (mm/d?)
568	REAL(r_std), DIMENSION(npts), INTENT(in) :: snowfall_daily
569	! snowmass_daily (kg/m2)
570	REAL(r_std), DIMENSION(npts), INTENT(in) :: snowmass_daily
571	! "14days" 2-meter temperatures (K)
572	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_14
573	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: sla_calc
574	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: wshtotsum
575	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: sr_ugb
576	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: compt_ugb
577	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: nb_ani
578	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: grazed_frac
579	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: import_yield
580	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: sla_age1
581	INTEGER(i_std), INTENT(in) :: day_of_year
582	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: N_limfert
583	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: when_growthinit_cut
584	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: nb_grazingdays
585	! REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: resp_hetero_litter_d
586	! REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: resp_hetero_soil_d
587	! top 5 layer grassland soil moisture for grazing
588	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: moiavail_daily
589	!! Daily moisture availability (0-1, unitless)
590	REAL(r_std),DIMENSION (npts), INTENT(in) :: tmc_topgrass_daily
591	REAL(r_std),DIMENSION (npts), INTENT(in) :: fc_grazing
592	REAL(r_std),DIMENSION (npts), INTENT(inout) :: after_snow
593	REAL(r_std),DIMENSION (npts), INTENT(inout) :: after_wet
594	REAL(r_std),DIMENSION (npts), INTENT(inout) :: wet1day
595	REAL(r_std),DIMENSION (npts), INTENT(inout) :: wet2day
596	!end gmjc
597	!! 0.4 Local variables
598
599	REAL(r_std), DIMENSION(npts,nvm,nelements) :: tot_bm_to_litter !! Total conversion of biomass to litter
600	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
601	REAL(r_std), DIMENSION(npts,nvm,nelements) :: tot_live_biomass !! Total living biomass
602	!! @tex $(gC m{-2})$ @endtex
603	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements) :: bm_alloc !! Biomass increase, i.e. NPP per plant part
604	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
605	REAL(r_std), DIMENSION(npts,nvm,nelements) :: tot_turnover !! Total turnover rate
606	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
607	REAL(r_std), DIMENSION(npts,nvm) :: tot_litter_soil_carb!! Total soil and litter carbon
608	!! @tex $(gC m^{-2})$ @endtex
609	REAL(r_std), DIMENSION(npts,nvm) :: tot_litter_carb !! Total litter carbon
610	!! @tex $(gC m^{-2})$ @endtex
611	REAL(r_std), DIMENSION(npts,nvm) :: tot_soil_carb !! Total soil carbon
612	!! @tex $(gC m^{-2})$ @endtex
613	REAL(r_std), DIMENSION(npts) :: carb_mass_variation !! Carbon Mass variation
614	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
615	REAL(r_std), DIMENSION(npts,nvm) :: cn_ind !! Crown area of individuals
616	!! @tex $(m^{2})$ @endtex
617	REAL(r_std), DIMENSION(npts,nvm) :: woodmass_ind !! Woodmass of individuals (gC)
618	REAL(r_std), DIMENSION(npts,nvm,nparts) :: f_alloc !! Fraction that goes into plant part
619	!! (0 to 1, unitless)
620	REAL(r_std), DIMENSION(npts) :: avail_tree !! Space availability for trees
621	!! (0 to 1, unitless)
622	REAL(r_std), DIMENSION(npts) :: avail_grass !! Space availability for grasses
623	!! (0 to 1, unitless)
624	INTEGER :: j,ivm,ivma,ipts,ilev,ilitt
625	REAL(r_std),DIMENSION(npts) :: prod10_total !! Total products remaining in the pool
626	!! after the annual release
627	!! @tex $(gC m^{-2})$ @endtex
628	REAL(r_std),DIMENSION(npts) :: prod100_total !! Total products remaining in the pool
629	!! after the annual release
630	!! @tex $(gC m^{-2})$ @endtex
631	REAL(r_std),DIMENSION(npts) :: cflux_prod_total !! Total flux from conflux and the 10/100
632	!! year-turnover pool
633	!! @tex $(gC m^{-2} year^{-1})$ @endtex
634	REAL(r_std),DIMENSION(npts) :: prod10_harvest_total!! Total products remaining in the pool
635	!! after the annual release
636	!! @tex $(gC m^{-2})$ @endtex
637	REAL(r_std),DIMENSION(npts) :: prod100_harvest_total!! Total products remaining in the pool
638	!! after the annual release
639	!! @tex $(gC m^{-2})$ @endtex
640	REAL(r_std),DIMENSION(npts) :: cflux_prod_harvest_total!! Total flux from conflux and the 10/100
641	!! year-turnover pool
642	!! @tex $(gC m^{-2} year^{-1})$ @endtex
643	REAL(r_std),DIMENSION(npts,nvm) :: veget_cov_max_tmp !! "Maximal" coverage fraction of a PFT
644	REAL(r_std), DIMENSION(npts) :: area_land_m2 !! Land area within gridcel excluding water body [m2]
645	!! (LAI-> infinity) on ground (unitless)
646	!!!!! crop
647	REAL(r_std), DIMENSION(npts,nvm) :: crop_export !! Cropland export (harvest & straw)
648	!!!!! end crop, xuhui
649	REAL(r_std), DIMENSION(npts,nvm) :: mortality !! Fraction of individual dying this time
650	!! step (0 to 1, unitless)
651	REAL(r_std), DIMENSION(npts) :: vartmp !! Temporary variable used to add history
652	REAL(r_std), DIMENSION(npts,nvm) :: histvar !! History variables
653	REAL(r_std), DIMENSION(npts,nvm) :: lcc !! land cover change, i.e., loss of each PFT,
654	!! positive value indicates loss.
655	REAL(r_std),DIMENSION(npts,nvm) :: deflitsup_total
656	REAL(r_std),DIMENSION(npts,nvm) :: defbiosup_total
657	!glcc
658	REAL(r_std), DIMENSION(npts,nvm,nvmap) :: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
659	REAL(r_std), DIMENSION(npts,12) :: glccReal !! The "real" glcc matrix that we apply in the model
660	!! after considering the consistency between presribed
661	!! glcc matrix and existing vegetation fractions.
662	REAL(r_std), DIMENSION(npts,12) :: IncreDeficit !! "Increment" deficits, negative values mean that
663	!! there are not enough fractions in the source PFTs
664	!! /vegetations to target PFTs/vegetations. I.e., these
665	!! fraction transfers are presribed in LCC matrix but
666	!! not realized.
667	REAL(r_std), DIMENSION(npts,12) :: glccZero !! "Increment" deficits, negative values mean that
668	REAL(r_std), DIMENSION(npts) :: Deficit_pf2yf_final !!
669	REAL(r_std), DIMENSION(npts) :: Deficit_sf2yf_final !!
670	REAL(r_std), DIMENSION(npts) :: pf2yf_compen_sf2yf !!
671	REAL(r_std), DIMENSION(npts) :: sf2yf_compen_pf2yf !!
672	REAL(r_std), DIMENSION(npts) :: fuelfrac !!
673	REAL(r_std) :: valtmp
674
675	!gmjc lcchange of managed grassland
676	! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground
677	INTEGER(i_std) :: ier
678	LOGICAL :: l_error =.FALSE.
679	! variables to get closure carbon cycle for nbp
680	REAL(r_std), DIMENSION(npts,nvm) :: harvest_gm
681	REAL(r_std), DIMENSION(npts,nvm) :: ranimal_gm
682	REAL(r_std), DIMENSION(npts,nvm) :: ch4_pft_gm
683	REAL(r_std), DIMENSION(npts) :: ch4_gm
684	REAL(r_std), DIMENSION(npts,nvm) :: cinput_gm
685	REAL(r_std), DIMENSION(npts) :: co2_gm
686	REAL(r_std),DIMENSION(npts,nvm) :: veget_cov_max_gm
687	REAL(r_std), DIMENSION(npts) :: veget_exist_gm
688	REAL(r_std),DIMENSION(npts,nvm) :: n2o_pft_gm
689	REAL(r_std), DIMENSION(npts) :: n2o_gm
690	!end gmjc
691
692
693
694	REAL(r_std), DIMENSION(npts,7) :: outflux_sta,outflux_end
695	REAL(r_std), DIMENSION(npts,3) :: influx_sta,influx_end
696	REAL(r_std), DIMENSION(npts,4) :: pool_sta,pool_end
697	INTEGER(i_std) :: ind_biomass,ind_litter,ind_soil,ind_prod,ind_co2tobm,ind_gpp,ind_npp,&
698	ind_bm2lit,ind_resph,ind_respm,ind_respg,ind_convf,ind_cflux,ind_fire,&
699	l,m
700
701	REAL(r_std), DIMENSION(npts,nvm,nelements) :: mass_before
702	REAL(r_std), DIMENSION(npts,nvm,nelements) :: mass_change !! positive
703	! sign as inrease in mass
704	REAL(r_std), DIMENSION(npts,nvm,nelements) :: mass_after !! Temporary variable
705
706	REAL(r_std), DIMENSION(npts,nelements) :: mass_before_2rd
707	REAL(r_std), DIMENSION(npts,nelements) :: mass_change_2rd !! positive
708	! sign as inrease in mass
709	REAL(r_std), DIMENSION(npts,nelements) :: mass_after_2rd !! Temporary variable
710	REAL(r_std), DIMENSION(npts,nelements) :: mass_balance_2rd !! Temporary variable
711	INTEGER :: f2g=1, f2p=2, f2c=3
712	INTEGER :: g2f=4, g2p=5, g2c=6, p2f=7, p2g=8, p2c=9, c2f=10, c2g=11, c2p=12
713	!_ ================================================================================================================================
714
715	IF (printlev>=3) WRITE(numout,*) 'Entering stomate_lpj'
716
717	!! 1. Initializations
718
719	lcc(:,:) = zero
720	glccZero = zero
721	area_land_m2(:) = (resolution(:,1)resolution(:,2)contfrac(:)) !unit:m2
722	glccReal(:,:) = zero
723	glcc_pftmtc(:,:,:) = zero
724
725	mass_before = zero
726	mass_change = zero
727	mass_after = zero
728
729	!gmjc
730	IF (firstcall) THEN
731
732	firstcall = .FALSE.
733
734	!Config Key = GRM_ENABLE_GRAZING
735	!Config Desc = grazing allowed
736	!Config Def = n
737	!Config Help = flag for choose if you want animals or not.
738	!
739	enable_grazing = .FALSE.
740	CALL getin_p('GRM_ENABLE_GRAZING',enable_grazing)
741	WRITE (numout,*) 'enable_grazing',enable_grazing
742	WRITE (numout,*) 'manage',is_grassland_manag
743	WRITE (numout,*) 'cut',is_grassland_cut
744	WRITE (numout,*) 'grazed',is_grassland_grazed
745
746	!Config Key = STO_CHECK_MASSBALANCE
747	!Config Desc = Check for mass balance in stomate_lpj
748	!Config Def = n
749	!Config Help = flag to enable mass balance in stomate_lpj
750	!
751	STO_CHECK_MASSBALANCE = .TRUE.
752	CALL getin_p('STO_CHECK_MASSBALANCE',STO_CHECK_MASSBALANCE)
753
754	emideforest_biomass_accu(:,:,:,:) = zero
755	emideforest_litter_accu(:,:,:,:) = zero
756	bafrac_deforest_accu(:,:) = zero
757
758
759	IF (do_now_stomate_lcchange) THEN
760	IF (use_age_class) THEN
761	IF (SingleAgeClass) THEN
762	! CALL glcc_SinAgeC_firstday(npts,veget_cov_max,newvegfrac,harvest_matrix, &
763	! glccSecondShift,glccPrimaryShift,glccNetLCC,&
764	! glccReal,glcc_pft,glcc_pftmtc,IncreDeficit, &
765	! Deficit_pf2yf_final, Deficit_sf2yf_final, &
766	! pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
767
768	ELSE
769	! CALL glcc_MulAgeC_firstday(npts,veget_cov_max,newvegfrac,harvest_matrix, &
770	! glccSecondShift,glccPrimaryShift,glccNetLCC,&
771	! glccReal,glcc_pft,glcc_pftmtc,IncreDeficit, &
772	! Deficit_pf2yf_final, Deficit_sf2yf_final, &
773	! pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
774	ENDIF
775	! We put only the conversion of tree->Notree as deforestation
776	DO ipts = 1,npts
777	DO ivm = 1,nvm
778	DO ivma = 1,nvmap
779	IF (is_tree(ivm) .AND. .NOT. is_tree(start_index(ivma))) THEN
780	lcc(ipts,ivm) = lcc(ipts,ivm) + glcc_pftmtc(ipts,ivm,ivma)
781	ENDIF
782	ENDDO
783	ENDDO
784	ENDDO
785	ELSE ! (.NOT. use_age_class), i.e., net land cover change.
786	! note here veget_cov_max is the last-year veget_cov_max; vegetnew_firstday
787	! is the veget_cov_max of next year, we need lcc to be >0 where forest
788	! area decreased, i.e., when vegetnew_firstday < veget_cov_max
789	lcc(:,:) = veget_cov_max(:,:) - vegetnew_firstday(:,:)
790	ENDIF
791
792	IF (.NOT. allow_deforest_fire) lcc(:,:) = zero
793
794	! we creat the proxy that's needed for deforestation fire simulation.
795	DO ipts = 1,npts
796	DO ivm = 1,nvm
797	deforest_litter_remain(ipts,:,ivm,:,:) = litter(ipts,:,ivm,:,:)*lcc(ipts,ivm)
798	deforest_biomass_remain(ipts,ivm,:,:) = biomass(ipts,ivm,:,:)*lcc(ipts,ivm)
799	def_fuel_1hr_remain(ipts,ivm,:,:) = fuel_1hr(ipts,ivm,:,:)*lcc(ipts,ivm)
800	def_fuel_10hr_remain(ipts,ivm,:,:) = fuel_10hr(ipts,ivm,:,:)*lcc(ipts,ivm)
801	def_fuel_100hr_remain(ipts,ivm,:,:) = fuel_100hr(ipts,ivm,:,:)*lcc(ipts,ivm)
802	def_fuel_1000hr_remain(ipts,ivm,:,:) = fuel_1000hr(ipts,ivm,:,:)*lcc(ipts,ivm)
803	ENDDO
804	ENDDO
805	ENDIF
806
807
808	END IF !firstcall
809
810	!! 1.1 Initialize variables to zero
811	co2_to_bm(:,:) = zero
812	co2_fire(:,:) = zero
813	npp_daily(:,:) = zero
814	resp_maint(:,:) = zero
815	resp_growth(:,:) = zero
816	harvest_above(:,:) = zero
817	bm_to_litter(:,:,:,:) = zero
818	cn_ind(:,:) = zero
819	woodmass_ind(:,:) = zero
820	turnover_daily(:,:,:,:) = zero
821	crop_export(:,:) = zero
822	deflitsup_total(:,:) = zero
823	defbiosup_total(:,:) = zero
824	!gmjc
825	!! Initialize gm variables for nbp to zero
826	harvest_gm(:,:) = zero
827	ranimal_gm(:,:) = zero
828	ch4_pft_gm(:,:) = zero
829	cinput_gm(:,:) = zero
830	co2_gm(:) = zero
831	ch4_gm(:) = zero
832	n2o_gm(:) = zero
833	n2o_pft_gm(:,:) = zero
834	veget_cov_max_gm(:,:) = zero
835	veget_exist_gm(:) = zero
836	!end gmjc
837
838	! GLUC
839	IncreDeficit = zero
840
841	!! 1.2 Initialize variables to veget_cov_max
842	veget_cov_max_tmp(:,:) = veget_cov_max(:,:)
843
844	!! 1.3 Calculate some vegetation characteristics
845	!! 1.3.1 Calculate some vegetation characteristics
846	! Calculate cn_ind (individual crown mass) and individual height from
847	! state variables if running DGVM or dynamic mortality in static cover mode
848	!?? Explain (maybe in the header once) why you mulitply with veget_cov_max in the DGVM
849	!?? and why you don't multiply with veget_cov_max in stomate.
850	IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
851	IF(ok_dgvm) THEN
852	WHERE (ind(:,:).GT.min_stomate)
853	woodmass_ind(:,:) = &
854	((biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
855	+biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon)) &
856	*veget_cov_max(:,:))/ind(:,:)
857	ENDWHERE
858	ELSE
859	WHERE (ind(:,:).GT.min_stomate)
860	woodmass_ind(:,:) = &
861	(biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
862	+biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon))/ind(:,:)
863	ENDWHERE
864	ENDIF
865
866	CALL crown (npts, PFTpresent, &
867	ind, biomass, woodmass_ind, &
868	veget_cov_max, cn_ind, height)
869	ENDIF
870
871	!! 1.3.2 Prescribe characteristics if the vegetation is not dynamic
872	! IF the DGVM is not activated, the density of individuals and their crown
873	! areas don't matter, but they should be defined for the case we switch on
874	! the DGVM afterwards. At the first call, if the DGVM is not activated,
875	! impose a minimum biomass for prescribed PFTs and declare them present.
876
877
878	!DSG mass conservation ========================================
879	mass_before(:,:,:) = SUM(biomass(:,:,:,:),DIM=3)
880
881	CALL prescribe (npts, &
882	veget_cov_max, dt_days, PFTpresent, everywhere, when_growthinit, &
883	biomass, leaf_frac, ind, cn_ind, co2_to_bm)
884
885
886	IF(STO_CHECK_MASSBALANCE) THEN
887	!DSG mass conservation ============================================
888	CALL stomate_check_mass_values(npts,biomass(:,:,:,:),'lpj: after prescribe')
889	WRITE (numout,*) 'No mass cons test; prescribe violates mass per design' !:DSG right?
890	!DSG mass_change(:,:,icarbon) = zero
891	!DSG CALL stomate_check_cons_mass(npts, nvm, nelements, & ! dimensions
892	!DSG mass_before(:,:,:), & ! mass before
893	!DSG SUM(biomass(:,:,:,:),DIM=3), & ! mass after
894	!DSG mass_change(:,:,:) & ! net of fluxes
895	!DSG )
896	ENDIF
897
898	!! 2. Climatic constraints for PFT presence and regenerativeness
899
900	! Call this even when DGVM is not activated so that "adapted" and "regenerate"
901	! are kept up to date for the moment when the DGVM is activated.
902	CALL constraints (npts, dt_days, &
903	t2m_month, t2m_min_daily,when_growthinit, &
904	adapted, regenerate, Tseason)
905
906
907	!! 3. Determine introduction and elimination of PTS based on climate criteria
908
909	IF ( ok_dgvm ) THEN
910
911	!DSG mass conservation ========================================
912	mass_before(:,:,:) = SUM(biomass(:,:,:,:),DIM=3)
913
914	!! 3.1 Calculate introduction and elimination
915	CALL pftinout (npts, dt_days, adapted, regenerate, &
916	neighbours, veget_cov_max, &
917	biomass, ind, cn_ind, age, leaf_frac, npp_longterm, lm_lastyearmax, senescence, &
918	PFTpresent, everywhere, when_growthinit, need_adjacent, RIP_time, &
919	co2_to_bm, &
920	avail_tree, avail_grass, &
921	!gmjc
922	sla_calc)
923	!end gmjc
924	IF (STO_CHECK_MASSBALANCE) THEN
925	CALL stomate_check_mass_values(npts,biomass(:,:,:,:), 'lpj: after pftinout')
926	!DSG mass conservation ============================================
927	WRITE (numout,*) 'MASS CONSERVATON: pftinout is designed to violate against mass conservation'
928	WRITE (numout,*) 'but it should be not detectable according to comments'
929	mass_change(:,:,icarbon) = zero
930	mass_after = SUM(biomass(:,:,:,:),DIM=3)
931
932	CALL stomate_check_cons_mass(lalo, mass_before(:,:,:), & ! mass before
933	mass_after, & ! mass after
934	mass_change(:,:,:), & ! net of fluxes
935	'lpj: after pftinout' )
936	ENDIF
937
938	!! 3.2 Reset attributes for eliminated PFTs.
939	! This also kills PFTs that had 0 leafmass during the last year. The message
940	! "... after pftinout" is misleading in this case.
941
942
943	!DSG mass conservation ========================================
944	mass_before(:,:,:) = SUM(biomass(:,:,:,:),DIM=3)
945
946	CALL kill (npts, 'pftinout ', lm_lastyearmax, &
947	ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
948	lai, age, leaf_age, leaf_frac, npp_longterm, &
949	when_growthinit, everywhere, veget_cov_max, bm_to_litter)
950
951	IF(STO_CHECK_MASSBALANCE) THEN
952	CALL stomate_check_mass_values(npts,biomass(:,:,:,:), 'lpj: after kill')
953	!DSG mass conservation ============================================
954	mass_change(:,:,icarbon) = SUM(bm_to_litter(:,:,:,icarbon),DIM=3)
955	mass_after = SUM(biomass(:,:,:,:),DIM=3)
956	CALL stomate_check_cons_mass(lalo, mass_before(:,:,:), & ! mass before
957	mass_after, & ! mass after
958	mass_change(:,:,:), & ! net of fluxes
959	'lpj: after kill' )
960	ENDIF
961
962	!! 3.3 Calculate woodmass of individual tree
963	IF(ok_dgvm) THEN
964	WHERE ((ind(:,:).GT.min_stomate))
965	woodmass_ind(:,:) = &
966	((biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
967	+biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon)) &
968	*veget_cov_max(:,:))/ind(:,:)
969	ENDWHERE
970	ELSE
971	WHERE (ind(:,:).GT.min_stomate)
972	woodmass_ind(:,:) = &
973	(biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
974	+biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon))/ind(:,:)
975	ENDWHERE
976	ENDIF
977
978	! Calculate crown area and diameter for all PFTs (including the newly established)
979	CALL crown (npts, PFTpresent, &
980	ind, biomass, woodmass_ind, &
981	veget_cov_max, cn_ind, height)
982
983	ENDIF
984
985	!! 4. Phenology
986
987	!! 4.1 Write values to history file
988	! Current values for ::when_growthinit
989	CALL xios_orchidee_send_field("WHEN_GROWTHINIT",when_growthinit)
990
991	CALL histwrite_p (hist_id_stomate, 'WHEN_GROWTHINIT', itime, when_growthinit, npts*nvm, horipft_index)
992
993	! Set and write values for ::PFTpresent
994	WHERE(PFTpresent)
995	histvar=un
996	ELSEWHERE
997	histvar=zero
998	ENDWHERE
999
1000	CALL xios_orchidee_send_field("PFTPRESENT",histvar)
1001
1002	CALL histwrite_p (hist_id_stomate, 'PFTPRESENT', itime, histvar, npts*nvm, horipft_index)
1003
1004	! Set and write values for gdd_midwinter
1005	WHERE(gdd_midwinter.EQ.undef)
1006	histvar=val_exp
1007	ELSEWHERE
1008	histvar=gdd_midwinter
1009	ENDWHERE
1010
1011	CALL xios_orchidee_send_field("GDD_MIDWINTER",histvar)
1012
1013	CALL histwrite_p (hist_id_stomate, 'GDD_MIDWINTER', itime, histvar, npts*nvm, horipft_index)
1014
1015	! Set and write values for gdd_m5_dormance
1016	WHERE(gdd_m5_dormance.EQ.undef)
1017	histvar=val_exp
1018	ELSEWHERE
1019	histvar=gdd_m5_dormance
1020	ENDWHERE
1021
1022	CALL xios_orchidee_send_field('GDD_M5_DORMANCE',histvar)
1023	CALL histwrite_p (hist_id_stomate, 'GDD_M5_DORMANCE', itime, histvar, npts*nvm, horipft_index)
1024
1025	! Set and write values for ncd_dormance
1026	WHERE(ncd_dormance.EQ.undef)
1027	histvar=val_exp
1028	ELSEWHERE
1029	histvar=ncd_dormance
1030	ENDWHERE
1031
1032	CALL xios_orchidee_send_field("NCD_DORMANCE",histvar)
1033
1034	CALL histwrite_p (hist_id_stomate, 'NCD_DORMANCE', itime, histvar, npts*nvm, horipft_index)
1035
1036	!DSG mass conservation ========================================
1037	mass_before(:,:,:) = SUM(biomass(:,:,:,:),DIM=3)
1038
1039	!! 4.2 Calculate phenology
1040	CALL phenology (npts, dt_days, PFTpresent, &
1041	veget_cov_max, &
1042	t2m_longterm, t2m_month, t2m_week, gpp_daily, &
1043	maxmoiavail_lastyear, minmoiavail_lastyear, &
1044	moiavail_month, moiavail_week, &
1045	gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
1046	senescence, time_hum_min, &
1047	biomass, leaf_frac, leaf_age, &
1048	when_growthinit, co2_to_bm, &
1049	pdlai, slai, deltai, ssla, &
1050	begin_leaves, &!)
1051	!gmjc
1052	sla_calc)
1053	!end gmjc
1054
1055	IF (STO_CHECK_MASSBALANCE) THEN
1056
1057	CALL stomate_check_mass_values(npts,biomass(:,:,:,:), 'lpj: after phenology_prognostic')
1058	WRITE (numout,*)'no mass conservation check after phenology: as the routine is designed to violate mass'
1059	!DSG as we sometimes have the case 'There is leaf or root carbon that
1060	!should not be here, something could have gone' see stomate_phenology.f90
1061
1062	!DSG mass conservation ========================================
1063	mass_before(:,:,:) = SUM(biomass(:,:,:,:),DIM=3)
1064	ENDIF
1065
1066	!! 5. Allocate C to different plant parts
1067	!WRITE(numout,*) 'slai before lpj_alloc: ',slai(1,12:14)
1068	CALL alloc (npts, dt_days, &
1069	lai, veget_cov_max, senescence, when_growthinit, &
1070	moiavail_week, tsoil_month, soilhum_month, &
1071	biomass, age, leaf_age, leaf_frac, rprof, f_alloc, &!)
1072	deltai, ssla, & !added for crop by xuhui
1073	!gmjc
1074	sla_calc, when_growthinit_cut)
1075	!end gmjc
1076
1077	!! 5.1. Recalculate lai
1078	!! This should be done whenever biomass is modified
1079	CALL setlai(biomass, sla_calc, slai)
1080
1081	!! 6. NPP, maintenance and growth respiration
1082	!! 6.1 Calculate NPP and respiration terms
1083	CALL npp_calc (npts, dt_days, &
1084	PFTpresent, &
1085	t2m_daily, tsoil_daily, lai, rprof, &
1086	gpp_daily, f_alloc, bm_alloc, resp_maint_part,&
1087	biomass, leaf_age, leaf_frac, age, &
1088	resp_maint, resp_growth, npp_daily, &!)
1089	! for crop bm_alloc
1090	!!! crop variables
1091	in_cycle, deltai, dltaisen, ssla, pgrain, deltgrain, reprac, &
1092	nger, nlev, ndrp, nlax, nmat, nrec, &
1093	c_reserve, c_leafb, slai, tday_counter, veget_cov_max, &
1094	!!! end crop, xuhui
1095	!gmjc
1096	sla_calc, sla_age1,N_limfert)
1097	!end gmjc
1098
1099
1100	!! 6.2 Kill slow growing PFTs in DGVM or STOMATE with constant mortality
1101	IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
1102	CALL kill (npts, 'npp ', lm_lastyearmax, &
1103	ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
1104	lai, age, leaf_age, leaf_frac, npp_longterm, &
1105	when_growthinit, everywhere, veget_cov_max, bm_to_litter)
1106
1107	!! 6.2.1 Update wood biomass
1108	! For the DGVM
1109	IF(ok_dgvm) THEN
1110	WHERE (ind(:,:).GT.min_stomate)
1111	woodmass_ind(:,:) = &
1112	((biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
1113	+ biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon)) &
1114	*veget_cov_max(:,:))/ind(:,:)
1115	ENDWHERE
1116
1117	! For all pixels with individuals
1118	ELSE
1119	WHERE (ind(:,:).GT.min_stomate)
1120	woodmass_ind(:,:) = &
1121	(biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
1122	+ biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))/ind(:,:)
1123	ENDWHERE
1124	ENDIF ! ok_dgvm
1125
1126	!! 6.2.2 New crown area and maximum vegetation cover after growth
1127	CALL crown (npts, PFTpresent, &
1128	ind, biomass, woodmass_ind,&
1129	veget_cov_max, cn_ind, height)
1130
1131	ENDIF ! ok_dgvm
1132
1133	!! 7. fire
1134	!! 7.1. Burn PFTs
1135	!CALL fire (npts, dt_days, litterpart, &
1136	! litterhum_daily, t2m_daily, lignin_struc, veget_cov_max, &
1137	! fireindex, firelitter, biomass, ind, &
1138	! litter, dead_leaves, bm_to_litter, &
1139	! co2_fire, MatrixA)
1140
1141	!gmjc update available and not available litter for grazing litter
1142	!spitfire
1143
1144	!disable_fire and allow_deforest_fire are defined as constants in src_parameters/constantes_var.f90
1145	!disable_fire to DISABLE fire when being TRUE
1146	!allow_deforest_fire to activate deforestation fire module when being TRUE
1147	IF(.NOT.disable_fire) THEN
1148	CALL spitfire(npts, dt_days, veget_cov_max,resolution,contfrac, &
1149	PFTpresent,t2m_min_daily,t2m_max_daily, &
1150	precip_daily,wspeed_daily,soilhum_daily(:,1), &
1151	lightn,litter(:,:,:,:,icarbon),ni_acc,fire_numday, &
1152	fuel_1hr(:,:,:,icarbon),fuel_10hr(:,:,:,icarbon),fuel_100hr(:,:,:,icarbon), &
1153	fuel_1000hr(:,:,:,icarbon),ind,biomass(:,:,:,icarbon),popd,a_nd,height, &
1154	read_observed_ba, observed_ba, read_cf_fine,cf_fine, &
1155	read_cf_coarse,cf_coarse,read_ratio_flag, &
1156	ratio_flag,read_ratio,ratio,date, &
1157	bm_to_litter(:,:,:,icarbon),co2_fire, &
1158	lcc,bafrac_deforest_accu,emideforest_litter_accu(:,:,:,icarbon),emideforest_biomass_accu(:,:,:,icarbon),&
1159	deforest_litter_remain(:,:,:,:,icarbon),deforest_biomass_remain(:,:,:,icarbon),&
1160	def_fuel_1hr_remain(:,:,:,icarbon),def_fuel_10hr_remain(:,:,:,icarbon),&
1161	def_fuel_100hr_remain(:,:,:,icarbon),def_fuel_1000hr_remain(:,:,:,icarbon))
1162	ELSE
1163	ni_acc=0.
1164	fire_numday=0.
1165
1166	ENDIF
1167	!endspit
1168	! after fire burning
1169	litter_avail(:,:,:) = litter(:,:,:,iabove,icarbon) * &
1170	litter_avail_frac(:,:,:)
1171	litter_not_avail(:,:,:) = litter(:,:,:,iabove,icarbon) * &
1172	(1.0 - litter_avail_frac(:,:,:))
1173	!end gmjc
1174	!! 7.2 Kill PFTs in DGVM
1175	IF ( ok_dgvm ) THEN
1176
1177	! reset attributes for eliminated PFTs
1178	CALL kill (npts, 'fire ', lm_lastyearmax, &
1179	ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
1180	lai, age, leaf_age, leaf_frac, npp_longterm, &
1181	when_growthinit, everywhere, veget_cov_max, bm_to_litter)
1182
1183	ENDIF ! ok_dgvm
1184
1185	IF (STO_CHECK_MASSBALANCE) THEN
1186	CALL stomate_check_mass_values(npts,biomass(:,:,:,:), 'lpj: after fire')
1187	ENDIF
1188	!! 8. Tree mortality
1189	!DSG mass conservation ========================================
1190	mass_before(:,:,:) = SUM(biomass(:,:,:,:),DIM=3) + SUM(bm_to_litter(:,:,:,:),DIM=3)
1191
1192	! Does not depend on age, therefore does not change crown area.
1193	CALL gap (npts, dt_days, &
1194	npp_longterm, turnover_longterm, lm_lastyearmax, &
1195	PFTpresent, biomass, ind, bm_to_litter, mortality, t2m_min_daily, Tmin_spring_time, &!)
1196	!gmjc
1197	sla_calc)
1198	!end gmjc
1199
1200	IF ( ok_dgvm ) THEN
1201
1202	! reset attributes for eliminated PFTs
1203	CALL kill (npts, 'gap ', lm_lastyearmax, &
1204	ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
1205	lai, age, leaf_age, leaf_frac, npp_longterm, &
1206	when_growthinit, everywhere, veget_cov_max, bm_to_litter)
1207
1208	ENDIF
1209
1210	IF(STO_CHECK_MASSBALANCE) THEN
1211	!DSG debug: ===================================================
1212	CALL stomate_check_mass_values(npts,biomass(:,:,:,:), 'lpj: after gap n kill')
1213
1214	!DSG mass conservation ============================================
1215	mass_change(:,:,icarbon) = -SUM(bm_to_litter(:,:,:,icarbon),DIM=3)
1216
1217	mass_after = SUM(biomass(:,:,:,:),DIM=3)
1218	CALL stomate_check_cons_mass(lalo, mass_before(:,:,:), & ! mass before
1219	mass_after, & ! mass after
1220	mass_change(:,:,:), & ! net of fluxes
1221	'lpj: after gap n kill')
1222	ENDIF
1223
1224
1225	!! =====================================
1226	!! DSG: mass conservation issue: START
1227
1228	!! 10. Leaf senescence, new lai and other turnover processes
1229	!DSG mass conservation ========================================
1230	mass_before(:,:,:) = SUM(biomass(:,:,:,:),DIM=3) + SUM(bm_to_litter(:,:,:,:),DIM=3)
1231
1232	CALL turn (npts, dt_days, PFTpresent, &
1233	herbivores, &
1234	maxmoiavail_lastyear, minmoiavail_lastyear, &
1235	moiavail_week, moiavail_month,t2m_longterm, t2m_month, t2m_week, veget_cov_max, &
1236	gdd_from_growthinit, leaf_age, leaf_frac, age, lai, biomass, &
1237	turnover_daily, senescence,turnover_time, &!)
1238	nrec,crop_export, &
1239	!gmjc
1240	sla_calc)
1241	!end gmjc
1242	!! 10. Light competition
1243
1244	!! If not using constant mortality then kill with light competition
1245	! IF ( ok_dgvm .OR. .NOT.(lpj_gap_const_mort) ) THEN
1246	IF ( ok_dgvm ) THEN
1247
1248	!! 10.1 Light competition
1249	CALL light (npts, dt_days, &
1250	veget_cov_max, fpc_max, PFTpresent, cn_ind, lai, maxfpc_lastyear, &
1251	lm_lastyearmax, ind, biomass, veget_lastlight, bm_to_litter, mortality, &!)
1252	!gmjc
1253	sla_calc)
1254	!end gmjc
1255
1256	!! 10.2 Reset attributes for eliminated PFTs
1257	CALL kill (npts, 'light ', lm_lastyearmax, &
1258	ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
1259	lai, age, leaf_age, leaf_frac, npp_longterm, &
1260	when_growthinit, everywhere, veget_cov_max, bm_to_litter)
1261
1262	ENDIF
1263	IF(STO_CHECK_MASSBALANCE) THEN
1264	!DSG debug: ===================================================
1265	CALL stomate_check_mass_values(npts,biomass(:,:,:,:), 'lpj: after turn')
1266
1267	!DSG mass conservation ============================================
1268	mass_change(:,:,icarbon) = -SUM(turnover_daily(:,:,:,icarbon),DIM=3) &
1269	-SUM(bm_to_litter(:,:,:,icarbon),DIM=3) ! in case ok_dvgm=true
1270
1271	mass_after = SUM(biomass(:,:,:,:),DIM=3)
1272	CALL stomate_check_cons_mass(lalo, mass_before(:,:,:), & ! mass before
1273	mass_after, & ! mass afte
1274	mass_change(:,:,:), & ! net of fluxes
1275	'lpj: after turn')
1276	ENDIF
1277
1278	!! DSG: mass conservation issue: END
1279	!! =====================================
1280
1281	!! 11. Establishment of saplings
1282
1283	IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort ) THEN
1284	!DSG mass conservation ========================================
1285	mass_before(:,:,:) = SUM(biomass(:,:,:,:),DIM=3) + SUM(bm_to_litter(:,:,:,:),DIM=3)
1286
1287	!! 11.1 Establish new plants
1288	CALL establish (npts, lalo, dt_days, PFTpresent, regenerate, &
1289	neighbours, resolution, need_adjacent, herbivores, &
1290	precip_lastyear, gdd0_lastyear, lm_lastyearmax, &
1291	cn_ind, lai, avail_tree, avail_grass, npp_longterm, &
1292	leaf_age, leaf_frac, &
1293	ind, biomass, age, everywhere, co2_to_bm, veget_cov_max, woodmass_ind, &
1294	mortality, bm_to_litter, &
1295	!gmjc
1296	sla_calc)
1297	!end gmjc
1298
1299	!! 11.2 Calculate new crown area (and maximum vegetation cover)
1300	CALL crown (npts, PFTpresent, &
1301	ind, biomass, woodmass_ind, &
1302	veget_cov_max, cn_ind, height)
1303
1304	IF(STO_CHECK_MASSBALANCE) THEN
1305	!DSG debug: ===================================================
1306	CALL stomate_check_mass_values(npts,biomass(:,:,:,:), 'lpj: after establish')
1307
1308	!DSG mass conservation ============================================
1309	mass_change(:,:,icarbon) = -SUM(bm_to_litter(:,:,:,icarbon),DIM=3)
1310
1311	mass_after = SUM(biomass,DIM=3)
1312	CALL stomate_check_cons_mass(lalo, mass_before, & ! mass before
1313	mass_after, & ! mass after
1314	mass_change, & ! net of fluxes
1315	'lpj: after establish')
1316	ENDIF
1317
1318	ENDIF
1319	!gmjc Grassland_management
1320	!
1321	! 13 calculate grazing by animals or cutting for forage
1322	!
1323	IF (enable_grazing) THEN
1324	CALL main_grassland_management(&
1325	npts, lalo, neighbours, resolution, contfrac, &
1326	dt_days , &
1327	day_of_year , &
1328	t2m_daily , &
1329	t2m_min_daily , &
1330	t2m_14 , &
1331	tsurf_daily , &
1332	snowfall_daily , &
1333	biomass , &
1334	bm_to_litter , &
1335	litter , &
1336	litter_avail , &
1337	litter_not_avail , &
1338	!spitfire
1339	fuel_1hr(:,:,:,icarbon), &
1340	fuel_10hr(:,:,:,icarbon), &
1341	fuel_100hr(:,:,:,icarbon), &
1342	fuel_1000hr(:,:,:,icarbon), &
1343	!end spitfire
1344	.TRUE. , &
1345	EndOfYear , &
1346	when_growthinit_cut, nb_grazingdays, &
1347	lai,sla_calc,leaf_age,leaf_frac, &
1348	wshtotsum,sr_ugb,compt_ugb, &
1349	nb_ani,grazed_frac,import_yield,N_limfert, &
1350	moiavail_daily,tmc_topgrass_daily,fc_grazing, snowmass_daily, &
1351	after_snow, after_wet, wet1day, wet2day, &
1352	harvest_gm, ranimal_gm, ch4_pft_gm, cinput_gm, n2o_pft_gm)
1353	ENDIF
1354	!end gmjc
1355	!! 12. Calculate final LAI and vegetation cover
1356
1357	![chaoyue] veget_cov_max_tmp is used as veget_cov_max_old in cover SUBROUTINE
1358	CALL cover (npts, cn_ind, ind, biomass, &
1359	veget_cov_max, veget_cov_max_tmp, lai, &
1360	litter, litter_avail, litter_not_avail, carbon, &
1361	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1362	turnover_daily, bm_to_litter, &
1363	co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, gpp_daily, &
1364	deepC_a, deepC_s,deepC_p)
1365
1366	!! 13. Update litter pools to account for harvest
1367
1368	! the whole litter stuff:
1369	! litter update, lignin content, PFT parts, litter decay,
1370	! litter heterotrophic respiration, dead leaf soil cover.
1371	! No vertical discretisation in the soil for litter decay.\n
1372	! added by shilong for harvest
1373	IF(harvest_agri) THEN !!! DO NOT activate harves_agri if ok_LAIdev
1374	CALL harvest(npts, dt_days, veget_cov_max, &
1375	bm_to_litter, turnover_daily, &
1376	harvest_above)
1377	ENDIF
1378
1379
1380	! We initialize these variables to zero. If LUC modules are called, their
1381	! values will be changed. Otherwise they stay as zero as is expected.
1382	convflux(:,:) = zero
1383	cflux_prod10(:,:) = zero
1384	cflux_prod100(:,:) = zero
1385	IF (do_now_stomate_lcchange) THEN
1386
1387	! Initialize LUC specific variables
1388	prod10(:,0,:) = zero
1389	prod100(:,0,:) = zero
1390
1391	IF (use_age_class) THEN
1392
1393	! Build the constants specific for gross land use change
1394	CALL stomate_gluc_constants_init
1395
1396	WRITE(numout,*) 'Buiding indecies for age classes:'
1397	WRITE(numout,*) 'Buiding indecies for tress:'
1398	WRITE(numout,*) indall_tree
1399	WRITE(numout,*) indagec_tree
1400	WRITE(numout,*) indold_tree
1401	WRITE(numout,*) 'Buiding indecies for natural grass:'
1402	WRITE(numout,*) indall_grass
1403	WRITE(numout,*) indagec_grass
1404	WRITE(numout,*) indold_grass
1405	WRITE(numout,*) 'Buiding indecies for pasture:'
1406	WRITE(numout,*) indall_pasture
1407	WRITE(numout,*) indagec_pasture
1408	WRITE(numout,*) indold_pasture
1409	WRITE(numout,*) 'Buiding indecies for crop:'
1410	WRITE(numout,*) indall_crop
1411	WRITE(numout,*) indagec_crop
1412	WRITE(numout,*) indold_crop
1413	WRITE(numout,*) 'Buiding indecies for bioe1:'
1414	WRITE(numout,*) indall_bioe1
1415	WRITE(numout,*) indagec_bioe1
1416	WRITE(numout,*) indold_bioe1
1417
1418	IF (gluc_use_harvest_biomass) THEN
1419	! calculate the fuel fraction when input is the harvest biomass
1420	fuelfrac = harvest_biomass(:,3)
1421	ELSE
1422	fuelfrac = zero
1423	ENDIF
1424
1425	IF (SingleAgeClass) THEN
1426	CALL fharvest_SinAgeC (npts, dt_days, harvest_matrix,newvegfrac, &
1427	fuelfrac, &
1428	glccZero,glccZero,glccZero,&
1429	def_fuel_1hr_remain, def_fuel_10hr_remain, &
1430	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
1431	deforest_litter_remain, deforest_biomass_remain, &
1432	convflux, cflux_prod10, cflux_prod100, &
1433	glccReal, IncreDeficit, glcc_pft, glcc_pftmtc, &
1434	veget_cov_max, prod10, prod100, flux10, flux100, &
1435	PFTpresent, senescence, moiavail_month, moiavail_week, &
1436	gpp_week, ngd_minus5, resp_maint, resp_growth, &
1437	resp_hetero, npp_daily, when_growthinit, npp_longterm, &
1438	ind, lm_lastyearmax, everywhere, age, &
1439	co2_to_bm, gpp_daily, co2_fire, &
1440	time_hum_min, gdd_midwinter, gdd_from_growthinit, &
1441	gdd_m5_dormance, ncd_dormance, &
1442	lignin_struc, carbon, leaf_frac, &
1443	deepC_a, deepC_s, deepC_p, &
1444	leaf_age, bm_to_litter, biomass, litter, &
1445	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
1446
1447	CALL glcc_SinAgeC (npts, dt_days, glccZero,newvegfrac, &
1448	glccSecondShift,glccPrimaryShift,glccNetLCC,&
1449	def_fuel_1hr_remain, def_fuel_10hr_remain, &
1450	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
1451	deforest_litter_remain, deforest_biomass_remain, &
1452	convflux, cflux_prod10, cflux_prod100, &
1453	glccReal, IncreDeficit, glcc_pft, glcc_pftmtc, &
1454	veget_cov_max, prod10, prod100, flux10, flux100, &
1455	PFTpresent, senescence, moiavail_month, moiavail_week, &
1456	gpp_week, ngd_minus5, resp_maint, resp_growth, &
1457	resp_hetero, npp_daily, when_growthinit, npp_longterm, &
1458	ind, lm_lastyearmax, everywhere, age, &
1459	co2_to_bm, gpp_daily, co2_fire, &
1460	time_hum_min, gdd_midwinter, gdd_from_growthinit, &
1461	gdd_m5_dormance, ncd_dormance, &
1462	lignin_struc, carbon, leaf_frac, &
1463	deepC_a, deepC_s, deepC_p, &
1464	leaf_age, bm_to_litter, biomass, litter, &
1465	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
1466
1467	! Multiple age class + wood harvest
1468	ELSE
1469
1470	IF (allow_forestry_harvest) THEN
1471
1472	IF (gluc_use_harvest_biomass) THEN
1473	CALL Get_harvest_matrix (npts,veget_cov_max,newvegfrac, &
1474	harvest_matrix, &
1475	biomass, harvest_biomass, area_land_m2, &
1476	glcc_pft,glcc_pftmtc, &
1477	Deficit_pf2yf_final, Deficit_sf2yf_final, &
1478	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
1479	ENDIF
1480
1481	CALL prepare_balance_check(outflux_sta,influx_sta,pool_sta, &
1482	veget_cov_max, &
1483	co2_to_bm,gpp_daily,npp_daily, &
1484	biomass,litter,carbon,prod10,prod100, &
1485	bm_to_litter,resp_maint,resp_growth,resp_hetero, &
1486	convflux,cflux_prod10,cflux_prod100,co2_fire)
1487
1488	CALL fharvest_MulAgeC (npts, dt_days, harvest_matrix,newvegfrac, &
1489	fuelfrac,&
1490	def_fuel_1hr_remain, def_fuel_10hr_remain, &
1491	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
1492	deforest_litter_remain, deforest_biomass_remain, &
1493	convflux, &
1494	glcc_pft, glcc_pftmtc, &
1495	veget_cov_max, prod10, prod100, &
1496	PFTpresent, senescence, moiavail_month, moiavail_week, &
1497	gpp_week, ngd_minus5, resp_maint, resp_growth, &
1498	resp_hetero, npp_daily, when_growthinit, npp_longterm, &
1499	ind, lm_lastyearmax, everywhere, age, &
1500	co2_to_bm, gpp_daily, co2_fire, &
1501	time_hum_min, gdd_midwinter, gdd_from_growthinit, &
1502	gdd_m5_dormance, ncd_dormance, &
1503	lignin_struc, carbon, leaf_frac, &
1504	deepC_a, deepC_s, deepC_p, &
1505	leaf_age, bm_to_litter, biomass, litter, &
1506	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
1507
1508	CALL prepare_balance_check(outflux_end,influx_end,pool_end, &
1509	veget_cov_max, &
1510	co2_to_bm,gpp_daily,npp_daily, &
1511	biomass,litter,carbon,prod10,prod100, &
1512	bm_to_litter,resp_maint,resp_growth,resp_hetero, &
1513	convflux,cflux_prod10,cflux_prod100,co2_fire)
1514
1515	CALL luc_balance_check(outflux_sta,influx_sta,pool_sta, &
1516	outflux_end,influx_end,pool_end, &
1517	npts,lalo,'wood harvest')
1518
1519	ENDIF
1520
1521	CALL prepare_balance_check(outflux_sta,influx_sta,pool_sta, &
1522	veget_cov_max, &
1523	co2_to_bm,gpp_daily,npp_daily, &
1524	biomass,litter,carbon,prod10,prod100, &
1525	bm_to_litter,resp_maint,resp_growth,resp_hetero, &
1526	convflux,cflux_prod10,cflux_prod100,co2_fire)
1527
1528	CALL glcc_MulAgeC (npts, dt_days, newvegfrac, &
1529	glccSecondShift,glccPrimaryShift,glccNetLCC, &
1530	def_fuel_1hr_remain, def_fuel_10hr_remain, &
1531	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
1532	deforest_litter_remain, deforest_biomass_remain, &
1533	convflux, &
1534	glccReal, IncreDeficit, glcc_pft, glcc_pftmtc, &
1535	veget_cov_max, prod10, prod100, &
1536	PFTpresent, senescence, moiavail_month, moiavail_week, &
1537	gpp_week, ngd_minus5, resp_maint, resp_growth, &
1538	resp_hetero, npp_daily, when_growthinit, npp_longterm, &
1539	ind, lm_lastyearmax, everywhere, age, &
1540	co2_to_bm, gpp_daily, co2_fire, &
1541	time_hum_min, gdd_midwinter, gdd_from_growthinit, &
1542	gdd_m5_dormance, ncd_dormance, &
1543	lignin_struc, carbon, leaf_frac, &
1544	deepC_a, deepC_s, deepC_p, &
1545	leaf_age, bm_to_litter, biomass, litter, &
1546	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
1547
1548	CALL prepare_balance_check(outflux_end,influx_end,pool_end, &
1549	veget_cov_max, &
1550	co2_to_bm,gpp_daily,npp_daily, &
1551	biomass,litter,carbon,prod10,prod100, &
1552	bm_to_litter,resp_maint,resp_growth,resp_hetero, &
1553	convflux,cflux_prod10,cflux_prod100,co2_fire)
1554
1555	CALL luc_balance_check(outflux_sta,influx_sta,pool_sta, &
1556	outflux_end,influx_end,pool_end, &
1557	npts,lalo,'gross land cover change')
1558
1559	IF (gluc_allow_trans_bioe) THEN
1560
1561	CALL prepare_balance_check(outflux_sta,influx_sta,pool_sta, &
1562	veget_cov_max, &
1563	co2_to_bm,gpp_daily,npp_daily, &
1564	biomass,litter,carbon,prod10,prod100, &
1565	bm_to_litter,resp_maint,resp_growth,resp_hetero, &
1566	convflux,cflux_prod10,cflux_prod100,co2_fire)
1567
1568	CALL glcc_bioe1 (npts, dt_days, newvegfrac, &
1569	glccSecondShift*0, &
1570	def_fuel_1hr_remain, def_fuel_10hr_remain, &
1571	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
1572	deforest_litter_remain, deforest_biomass_remain, &
1573	convflux, cflux_prod10, cflux_prod100, &
1574	glcc_pft, glcc_pftmtc, &
1575	veget_cov_max, prod10, prod100, flux10, flux100, &
1576	PFTpresent, senescence, moiavail_month, moiavail_week, &
1577	gpp_week, ngd_minus5, resp_maint, resp_growth, &
1578	resp_hetero, npp_daily, when_growthinit, npp_longterm, &
1579	ind, lm_lastyearmax, everywhere, age, &
1580	co2_to_bm, gpp_daily, co2_fire, &
1581	time_hum_min, gdd_midwinter, gdd_from_growthinit, &
1582	gdd_m5_dormance, ncd_dormance, &
1583	lignin_struc, carbon, leaf_frac, &
1584	deepC_a, deepC_s, deepC_p, &
1585	leaf_age, bm_to_litter, biomass, litter, &
1586	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
1587
1588	CALL prepare_balance_check(outflux_end,influx_end,pool_end, &
1589	veget_cov_max, &
1590	co2_to_bm,gpp_daily,npp_daily, &
1591	biomass,litter,carbon,prod10,prod100, &
1592	bm_to_litter,resp_maint,resp_growth,resp_hetero, &
1593	convflux,cflux_prod10,cflux_prod100,co2_fire)
1594
1595	CALL luc_balance_check(outflux_sta,influx_sta,pool_sta, &
1596	outflux_end,influx_end,pool_end, &
1597	npts,lalo,'glcc with bioenergy PFT')
1598
1599	ENDIF
1600
1601	ENDIF !(SingleAgeClass)
1602
1603	! clear the constants specific for gross land use change
1604	CALL stomate_gluc_constants_init_clear
1605
1606	ELSE ! .NOT. use_age_class
1607	IF (allow_deforest_fire) THEN
1608	![chaoyue] veget_cov_max is used as the old veget_cov_max in lcchange_deffire
1609	! veget_cov_max_new is used as the new veget_cov_max in lcchange_deffire
1610	CALL lcchange_deffire (npts, dt_days, veget_cov_max, veget_cov_max_new, &
1611	biomass, ind, age, PFTpresent, senescence, when_growthinit, everywhere, &
1612	co2_to_bm, bm_to_litter, turnover_daily, bm_sapl, cn_ind, &
1613	flux10,flux100, &
1614	prod10,prod100,&
1615	convflux,&
1616	cflux_prod10,cflux_prod100, leaf_frac,&
1617	npp_longterm, lm_lastyearmax, litter, litter_avail, litter_not_avail, &
1618	carbon,&
1619	deepC_a, deepC_s, deepC_p,&
1620	fuel_1hr,fuel_10hr,fuel_100hr,fuel_1000hr,&
1621	lcc,bafrac_deforest_accu,emideforest_litter_accu,emideforest_biomass_accu,&
1622	deflitsup_total,defbiosup_total)
1623	ELSE
1624	![chaoyue] veget_cov_max is used as veget_cov_max_old in lcchange_main
1625	! veget_cov_max_new is used as the new veget_cov_max in lcchange_main
1626	CALL lcchange_main (npts, dt_days, veget_cov_max_new, veget_cov_max, &
1627	biomass, ind, age, PFTpresent, senescence, when_growthinit, everywhere, &
1628	co2_to_bm, bm_to_litter, turnover_daily, bm_sapl, cn_ind,&
1629	flux10,flux100, &
1630	prod10,prod100,convflux, &
1631	cflux_prod10,cflux_prod100,leaf_frac,&
1632	npp_longterm, lm_lastyearmax, litter, litter_avail, litter_not_avail, &
1633	carbon, &
1634	deepC_a, deepC_s, deepC_p,&
1635	fuel_1hr,fuel_10hr,fuel_100hr,fuel_1000hr )
1636	ENDIF
1637	ENDIF ! (use_age_class)
1638
1639	!! Update the 10-year-turnover pool content following flux emission
1640	!! (linear decay (10%) of the initial carbon input)
1641	DO l = 0, 8
1642	m = 10 - l
1643	cflux_prod10(:,:) = cflux_prod10(:,:) + flux10(:,m,:)
1644	prod10(:,m,:) = prod10(:,m-1,:) - flux10(:,m-1,:)
1645	flux10(:,m,:) = flux10(:,m-1,:)
1646	ENDDO
1647
1648	cflux_prod10(:,:) = cflux_prod10(:,:) + flux10(:,1,:)
1649	flux10(:,1,:) = 0.1 * prod10(:,0,:)
1650	prod10(:,1,:) = prod10(:,0,:)
1651
1652	!! Update the 100-year-turnover pool content following flux emission\n
1653	DO l = 0, 98
1654	m = 100 - l
1655	cflux_prod100(:,:) = cflux_prod100(:,:) + flux100(:,m,:)
1656	prod100(:,m,:) = prod100(:,m-1,:) - flux100(:,m-1,:)
1657	flux100(:,m,:) = flux100(:,m-1,:)
1658	ENDDO
1659
1660	cflux_prod100(:,:) = cflux_prod100(:,:) + flux100(:,1,:)
1661	flux100(:,1,:) = 0.01 * prod100(:,0,:)
1662	prod100(:,1,:) = prod100(:,0,:)
1663	prod10(:,0,:) = zero
1664	prod100(:,0,:) = zero
1665
1666	! set the flag
1667	do_now_stomate_lcchange=.FALSE.
1668
1669	! Set the flag done_stomate_lcchange to be used in the end of sechiba_main to update the fractions.
1670	done_stomate_lcchange=.TRUE.
1671
1672	CALL stomate_check_mass_values(npts,biomass(:,:,:,:), 'lpj: after land use change')
1673
1674	ENDIF ! do_now_stomate_lcchange
1675
1676	! Update the status of age classes
1677	IF(EndOfYear) THEN
1678
1679	! start the mass balance check
1680	mass_before_2rd = zero
1681	mass_change_2rd = zero
1682	mass_after_2rd = zero
1683	mass_balance_2rd = zero
1684	pool_sta = zero
1685	influx_sta = zero
1686	outflux_sta = zero
1687
1688	ind_biomass = 1
1689	ind_litter = 2
1690	ind_soil = 3
1691	pool_sta(:,ind_biomass) = SUM(SUM(biomass(:,:,:,icarbon),DIM=3) * veget_cov_max,DIM=2)
1692	pool_sta(:,ind_litter) = SUM(SUM(SUM(litter(:,:,:,:,icarbon),DIM=4),DIM=2) * veget_cov_max,DIM=2)
1693	pool_sta(:,ind_soil) = SUM(SUM(carbon(:,:,:),DIM=2) * veget_cov_max,DIM=2)
1694	mass_before_2rd(:,icarbon) = SUM(pool_sta(:,:),DIM=2)
1695
1696	ind_co2tobm = 1
1697	ind_gpp = 2
1698	ind_npp = 3
1699	influx_sta(:,ind_co2tobm) = SUM(co2_to_bm*veget_cov_max,DIM=2)
1700	influx_sta(:,ind_gpp) = SUM(gpp_daily*veget_cov_max,DIM=2)
1701	influx_sta(:,ind_npp) = SUM(npp_daily*veget_cov_max,DIM=2)
1702
1703	ind_bm2lit = 1
1704	ind_respm = 2
1705	ind_respg = 3
1706	ind_resph = 4
1707	outflux_sta(:,ind_bm2lit) = SUM(SUM(bm_to_litter(:,:,:,icarbon),DIM=3) * veget_cov_max,DIM=2)
1708	outflux_sta(:,ind_respm) = SUM(resp_maint*veget_cov_max,DIM=2)
1709	outflux_sta(:,ind_respg) = SUM(resp_growth*veget_cov_max,DIM=2)
1710	outflux_sta(:,ind_resph) = SUM(resp_hetero*veget_cov_max,DIM=2)
1711
1712
1713	IF (use_age_class) THEN
1714	CALL age_class_distr(npts, lalo, resolution, bound_spa, &
1715	biomass, veget_cov_max, ind, &
1716	lm_lastyearmax, leaf_frac, co2_to_bm, &
1717	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1718	everywhere, litter, carbon, lignin_struc, &
1719	deepC_a, deepC_s, deepC_p, &
1720	bm_to_litter, PFTpresent, when_growthinit,&
1721	senescence, npp_longterm, gpp_daily, leaf_age, age, &
1722	gdd_from_growthinit, gdd_midwinter, time_hum_min, hum_min_dormance, &
1723	gdd_m5_dormance, &
1724	ncd_dormance, moiavail_month, moiavail_week, ngd_minus5, &
1725	gpp_week, resp_maint, resp_growth, npp_daily, resp_hetero)
1726	ENDIF
1727
1728	pool_end = zero
1729	influx_end = zero
1730	outflux_end = zero
1731
1732	ind_biomass = 1
1733	ind_litter = 2
1734	ind_soil = 3
1735	pool_end(:,ind_biomass) = SUM(SUM(biomass(:,:,:,icarbon),DIM=3) * veget_cov_max,DIM=2)
1736	pool_end(:,ind_litter) = SUM(SUM(SUM(litter(:,:,:,:,icarbon),DIM=4),DIM=2) * veget_cov_max,DIM=2)
1737	pool_end(:,ind_soil) = SUM(SUM(carbon(:,:,:),DIM=2) * veget_cov_max,DIM=2)
1738
1739	ind_co2tobm = 1
1740	ind_gpp = 2
1741	ind_npp = 3
1742	influx_end(:,ind_co2tobm) = SUM(co2_to_bm*veget_cov_max,DIM=2)
1743	influx_end(:,ind_gpp) = SUM(gpp_daily*veget_cov_max,DIM=2)
1744	influx_end(:,ind_npp) = SUM(npp_daily*veget_cov_max,DIM=2)
1745
1746
1747	ind_bm2lit = 1
1748	ind_respm = 2
1749	ind_respg = 3
1750	ind_resph = 4
1751	outflux_end(:,ind_bm2lit) = SUM(SUM(bm_to_litter(:,:,:,icarbon),DIM=3) * veget_cov_max,DIM=2)
1752	outflux_end(:,ind_respm) = SUM(resp_maint*veget_cov_max,DIM=2)
1753	outflux_end(:,ind_respg) = SUM(resp_growth*veget_cov_max,DIM=2)
1754	outflux_end(:,ind_resph) = SUM(resp_hetero*veget_cov_max,DIM=2)
1755
1756	CALL stomate_check_mass_values(npts,biomass(:,:,:,:), 'lpj: after age_class_distr')
1757
1758	! mass_change_2rd is the mass increase
1759	mass_change_2rd(:,icarbon) = SUM(influx_end(:,:),DIM=2) - SUM(influx_sta(:,:),DIM=2) + &
1760	+ SUM(outflux_sta(:,:),DIM=2) - SUM(outflux_end(:,:),DIM=2)
1761
1762	mass_after_2rd(:,icarbon) = SUM(pool_end(:,:),dim=2)
1763	mass_balance_2rd(:,icarbon) = mass_before_2rd(:,icarbon) - mass_after_2rd(:,icarbon) &
1764	+ mass_change_2rd(:,icarbon)
1765
1766	DO ipts = 1,npts
1767	IF (ABS(mass_balance_2rd(ipts,icarbon)) .GE. 1e-3) THEN
1768	WRITE (numout,*) 'FATAL Error'
1769	WRITE (numout,*) 'Mass conservation failed after age_class_distr'
1770	WRITE (numout,*) ' limit: ' , 1e-3
1771	WRITE (numout,*) ' Mismatch :' , mass_balance_2rd(ipts,icarbon)
1772	WRITE (numout,*) ' Coordinates :' , lalo(ipts,:)
1773	WRITE (numout,*) ' gridpoint: ' , ipts , ' of ngrids: ',npts
1774	STOP
1775	ENDIF
1776	ENDDO
1777	WRITE (numout,*) 'mass balance check successful after age_class_distr'
1778
1779	ENDIF
1780
1781
1782	!! 16. Recalculate lai
1783	!! This should be done whenever biomass is modified
1784	CALL setlai(biomass, sla_calc, slai)
1785	CALL setlai(biomass, sla_calc, lai)
1786
1787	!! 17. Calculate vcmax
1788	CALL vmax (npts, dt_days, &
1789	leaf_age, leaf_frac, &
1790	vcmax, &!)
1791	N_limfert)
1792
1793	!MM deplacement pour initialisation correcte des grandeurs cumulees :
1794	cflux_prod_total(:) = convflux(:,iwplcc) + cflux_prod10(:,iwplcc) + cflux_prod100(:,iwplcc)
1795	prod10_total(:)=SUM(prod10(:,:,iwplcc),dim=2)
1796	prod100_total(:)=SUM(prod100(:,:,iwplcc),dim=2)
1797
1798	cflux_prod_harvest_total = zero
1799	prod10_harvest_total = zero
1800	prod100_harvest_total = zero
1801
1802	!! 18. Total heterotrophic respiration
1803
1804	tot_soil_carb(:,:) = zero
1805	tot_litter_carb(:,:) = zero
1806	DO j=1,nvm
1807
1808	tot_litter_carb(:,j) = tot_litter_carb(:,j) + (litter(:,istructural,j,iabove,icarbon) + &
1809	& litter(:,imetabolic,j,iabove,icarbon) + &
1810	& litter(:,istructural,j,ibelow,icarbon) + litter(:,imetabolic,j,ibelow,icarbon))
1811
1812	tot_soil_carb(:,j) = tot_soil_carb(:,j) + (carbon(:,iactive,j) + &
1813	& carbon(:,islow,j)+ carbon(:,ipassive,j))
1814
1815	ENDDO
1816	tot_litter_soil_carb(:,:) = tot_litter_carb(:,:) + tot_soil_carb(:,:)
1817
1818	!!$ DO k = 1, nelements ! Loop over # elements
1819	!!$ tot_live_biomass(:,:,k) = biomass(:,:,ileaf,k) + biomass(:,:,isapabove,k) + biomass(:,:,isapbelow,k) +&
1820	!!$ & biomass(:,:,iheartabove,k) + biomass(:,:,iheartbelow,k) + &
1821	!!$ & biomass(:,:,iroot,k) + biomass(:,:,ifruit,k) + biomass(:,:,icarbres,k)
1822	!!$ END DO ! Loop over # elements
1823
1824	tot_live_biomass(:,:,:) = biomass(:,:,ileaf,:) + biomass(:,:,isapabove,:) + biomass(:,:,isapbelow,:) +&
1825	& biomass(:,:,iheartabove,:) + biomass(:,:,iheartbelow,:) + &
1826	& biomass(:,:,iroot,:) + biomass(:,:,ifruit,:) + biomass(:,:,icarbres,:)
1827
1828
1829	tot_turnover(:,:,:) = turnover_daily(:,:,ileaf,:) + turnover_daily(:,:,isapabove,:) + &
1830	& turnover_daily(:,:,isapbelow,:) + turnover_daily(:,:,iheartabove,:) + &
1831	& turnover_daily(:,:,iheartbelow,:) + turnover_daily(:,:,iroot,:) + &
1832	& turnover_daily(:,:,ifruit,:) + turnover_daily(:,:,icarbres,:)
1833
1834	tot_bm_to_litter(:,:,:) = bm_to_litter(:,:,ileaf,:) + bm_to_litter(:,:,isapabove,:) +&
1835	& bm_to_litter(:,:,isapbelow,:) + bm_to_litter(:,:,iheartbelow,:) +&
1836	& bm_to_litter(:,:,iheartabove,:) + bm_to_litter(:,:,iroot,:) + &
1837	& bm_to_litter(:,:,ifruit,:) + bm_to_litter(:,:,icarbres,:)
1838
1839	carb_mass_variation(:)=-carb_mass_total(:)
1840	carb_mass_total(:)=SUM((tot_live_biomass(:,:,icarbon)+tot_litter_carb+tot_soil_carb)*veget_cov_max,dim=2) + &
1841	& (prod10_total + prod100_total) + (prod10_harvest_total + prod100_harvest_total)
1842	carb_mass_variation(:)=carb_mass_total(:)+carb_mass_variation(:)
1843
1844	!! 17. Write history
1845
1846	CALL xios_orchidee_send_field("RESOLUTION_X",resolution(:,1))
1847	CALL xios_orchidee_send_field("RESOLUTION_Y",resolution(:,2))
1848	CALL xios_orchidee_send_field("CONTFRAC_STOMATE",contfrac(:))
1849	CALL xios_orchidee_send_field("T2M_MONTH",t2m_month)
1850	CALL xios_orchidee_send_field("T2M_WEEK",t2m_week)
1851	CALL xios_orchidee_send_field("HET_RESP",resp_hetero(:,:))
1852	CALL xios_orchidee_send_field("CO2_FIRE",co2_fire)
1853	CALL xios_orchidee_send_field("CO2_TAKEN",co2_to_bm)
1854	CALL xios_orchidee_send_field("LAI",lai)
1855	CALL xios_orchidee_send_field("VEGET_COV_MAX",veget_cov_max)
1856	CALL xios_orchidee_send_field("NPP_STOMATE",npp_daily)
1857	CALL xios_orchidee_send_field("GPP",gpp_daily)
1858	CALL xios_orchidee_send_field("IND",ind)
1859	CALL xios_orchidee_send_field("CN_IND",cn_ind)
1860	CALL xios_orchidee_send_field("WOODMASS_IND",woodmass_ind)
1861	CALL xios_orchidee_send_field("TOTAL_M",tot_live_biomass(:,:,icarbon))
1862	CALL xios_orchidee_send_field("MOISTRESS",moiavail_week)
1863	CALL xios_orchidee_send_field("LEAF_M",biomass(:,:,ileaf,icarbon))
1864	CALL xios_orchidee_send_field("SAP_M_AB",biomass(:,:,isapabove,icarbon))
1865	CALL xios_orchidee_send_field("SAP_M_BE",biomass(:,:,isapbelow,icarbon))
1866	CALL xios_orchidee_send_field("HEART_M_AB",biomass(:,:,iheartabove,icarbon))
1867	CALL xios_orchidee_send_field("HEART_M_BE",biomass(:,:,iheartbelow,icarbon))
1868	CALL xios_orchidee_send_field("ROOT_M",biomass(:,:,iroot,icarbon))
1869	CALL xios_orchidee_send_field("FRUIT_M",biomass(:,:,ifruit,icarbon))
1870	! HERE WE ADD A OUTPUT VARIABLE, NAMED CROPYIELD
1871	CALL xios_orchidee_send_field("CROPYIELD",biomass(:,:,ifruit,icarbon))
1872	CALL xios_orchidee_send_field("BIOMYIELD",tot_live_biomass(:,:,icarbon))
1873	! END DEFINE
1874	CALL xios_orchidee_send_field("RESERVE_M",biomass(:,:,icarbres,icarbon))
1875	! CALL xios_orchidee_send_field("TOTAL_TURN",tot_turnover)
1876	CALL xios_orchidee_send_field("LEAF_TURN",turnover_daily(:,:,ileaf,icarbon))
1877	CALL xios_orchidee_send_field("MAINT_RESP",resp_maint)
1878	CALL xios_orchidee_send_field("GROWTH_RESP",resp_growth)
1879	CALL xios_orchidee_send_field("SAP_AB_TURN",turnover_daily(:,:,isapabove,icarbon))
1880	CALL xios_orchidee_send_field("ROOT_TURN",turnover_daily(:,:,iroot,icarbon))
1881	CALL xios_orchidee_send_field("FRUIT_TURN",turnover_daily(:,:,ifruit,icarbon))
1882	CALL xios_orchidee_send_field("TOTAL_BM_LITTER",tot_bm_to_litter(:,:,icarbon))
1883	CALL xios_orchidee_send_field("LEAF_BM_LITTER",bm_to_litter(:,:,ileaf,icarbon))
1884	CALL xios_orchidee_send_field("SAP_AB_BM_LITTER",bm_to_litter(:,:,isapabove,icarbon))
1885	CALL xios_orchidee_send_field("SAP_BE_BM_LITTER",bm_to_litter(:,:,isapbelow,icarbon))
1886	CALL xios_orchidee_send_field("HEART_AB_BM_LITTER",bm_to_litter(:,:,iheartabove,icarbon))
1887	CALL xios_orchidee_send_field("HEART_BE_BM_LITTER",bm_to_litter(:,:,iheartbelow,icarbon))
1888	CALL xios_orchidee_send_field("ROOT_BM_LITTER",bm_to_litter(:,:,iroot,icarbon))
1889	CALL xios_orchidee_send_field("FRUIT_BM_LITTER",bm_to_litter(:,:,ifruit,icarbon))
1890	CALL xios_orchidee_send_field("RESERVE_BM_LITTER",bm_to_litter(:,:,icarbres,icarbon))
1891	CALL xios_orchidee_send_field("LITTER_STR_AB",litter(:,istructural,:,iabove,icarbon))
1892	CALL xios_orchidee_send_field("LITTER_MET_AB",litter(:,imetabolic,:,iabove,icarbon))
1893	CALL xios_orchidee_send_field("LITTER_STR_BE",litter(:,istructural,:,ibelow,icarbon))
1894	CALL xios_orchidee_send_field("LITTER_MET_BE",litter(:,imetabolic,:,ibelow,icarbon))
1895	CALL xios_orchidee_send_field("DEADLEAF_COVER",deadleaf_cover)
1896	CALL xios_orchidee_send_field("TOTAL_SOIL_CARB",tot_litter_soil_carb)
1897	CALL xios_orchidee_send_field("CARBON_ACTIVE",carbon(:,iactive,:))
1898	CALL xios_orchidee_send_field("CARBON_SLOW",carbon(:,islow,:))
1899	CALL xios_orchidee_send_field("CARBON_PASSIVE",carbon(:,ipassive,:))
1900	CALL xios_orchidee_send_field("LITTERHUM",litterhum_daily)
1901	CALL xios_orchidee_send_field("TURNOVER_TIME",turnover_time)
1902	! CALL xios_orchidee_send_field("PROD10",prod10)
1903	! CALL xios_orchidee_send_field("FLUX10",flux10)
1904	! CALL xios_orchidee_send_field("PROD100",prod100)
1905	! CALL xios_orchidee_send_field("FLUX100",flux100)
1906	CALL xios_orchidee_send_field("CONVFLUX",convflux)
1907	CALL xios_orchidee_send_field("CFLUX_PROD10",cflux_prod10)
1908	CALL xios_orchidee_send_field("CFLUX_PROD100",cflux_prod100)
1909	CALL xios_orchidee_send_field("HARVEST_ABOVE",SUM(harvest_above,DIM=2))
1910	CALL xios_orchidee_send_field("VCMAX",vcmax)
1911	CALL xios_orchidee_send_field("AGE",age)
1912	CALL xios_orchidee_send_field("HEIGHT",height)
1913	CALL xios_orchidee_send_field("FIREINDEX",fireindex(:,:))
1914	!gmjc
1915	CALL xios_orchidee_send_field("LITTER_STR_AVAIL",litter_avail(:,istructural,:))
1916	CALL xios_orchidee_send_field("LITTER_MET_AVAIL",litter_avail(:,imetabolic,:))
1917	CALL xios_orchidee_send_field("LITTER_STR_NAVAIL",litter_not_avail(:,istructural,:))
1918	CALL xios_orchidee_send_field("LITTER_MET_NAVAIL",litter_not_avail(:,imetabolic,:))
1919	CALL xios_orchidee_send_field("LITTER_STR_AVAILF",litter_avail_frac(:,istructural,:))
1920	CALL xios_orchidee_send_field("LITTER_MET_AVAILF",litter_avail_frac(:,imetabolic,:))
1921	IF (ANY(ok_LAIdev)) CALL xios_orchidee_send_field("N_LIMFERT",N_limfert)
1922	!calculate grassland co2 fluxes
1923	DO j=2,nvm
1924	IF ((.NOT. is_tree(j)) .AND. natural(j)) THEN
1925	veget_cov_max_gm(:,j) = veget_cov_max(:,j)
1926	ENDIF
1927	END DO ! nvm
1928	veget_exist_gm(:) = SUM(veget_cov_max_gm,dim=2)
1929	WHERE (veget_exist_gm(:) .GT. 0.0)
1930	co2_gm(:) = SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire &
1931	-harvest_gm-ranimal_gm-ch4_pft_gm+cinput_gm) &
1932	*veget_cov_max_gm,dim=2)/veget_exist_gm
1933	ch4_gm(:) = SUM(ch4_pft_gm*veget_cov_max_gm,dim=2)/veget_exist_gm
1934	n2o_gm(:) = SUM(n2o_pft_gm*veget_cov_max_gm,dim=2)/veget_exist_gm
1935	ELSEWHERE
1936	co2_gm(:) = zero
1937	ch4_gm(:) = zero
1938	n2o_gm(:) = zero
1939	ENDWHERE
1940	CALL xios_orchidee_send_field("CO2_GM",co2_gm)
1941	CALL xios_orchidee_send_field("CH4_GM",ch4_gm)
1942	CALL xios_orchidee_send_field("N2O_GM",n2o_gm)
1943	CALL xios_orchidee_send_field("N2O_PFT_GM",n2o_pft_gm)
1944	!end gmjc
1945
1946
1947	! ipcc history
1948	! Carbon stock transformed from gC/m2 into kgC/m2
1949	CALL xios_orchidee_send_field("cVeg",SUM(tot_live_biomass(:,:,icarbon)*veget_cov_max,dim=2)/1e3)
1950	CALL xios_orchidee_send_field("cLitter",SUM(tot_litter_carb*veget_cov_max,dim=2)/1e3)
1951	CALL xios_orchidee_send_field("cSoil",SUM(tot_soil_carb*veget_cov_max,dim=2)/1e3)
1952	CALL xios_orchidee_send_field("cProduct",(prod10_total + prod100_total + prod10_harvest_total + prod100_harvest_total)/1e3)
1953	CALL xios_orchidee_send_field("cMassVariation",carb_mass_variation/1e3/one_day)
1954
1955	CALL xios_orchidee_send_field("lai_ipcc",SUM(lai*veget_cov_max,dim=2)) ! m2/m2
1956
1957	! Carbon fluxes transformed from gC/m2/d into kgC/m2/s
1958	CALL xios_orchidee_send_field("gpp_ipcc",SUM(gpp_daily*veget_cov_max,dim=2)/1e3/one_day)
1959	CALL xios_orchidee_send_field("ra",SUM((resp_maint+resp_growth)*veget_cov_max,dim=2)/1e3/one_day)
1960	CALL xios_orchidee_send_field("npp_ipcc",SUM(npp_daily*veget_cov_max,dim=2)/1e3/one_day)
1961	CALL xios_orchidee_send_field("rh",SUM(resp_hetero*veget_cov_max,dim=2)/1e3/one_day)
1962	CALL xios_orchidee_send_field("fFire",SUM(co2_fire*veget_cov_max,dim=2)/1e3/one_day)
1963	CALL xios_orchidee_send_field("fHarvest",SUM(harvest_above,DIM=2)/1e3/one_day)
1964	CALL xios_orchidee_send_field("fLuc",cflux_prod_total/1e3/one_day)
1965	CALL xios_orchidee_send_field("fWoodharvest",cflux_prod_harvest_total/1e3/one_day)
1966	CALL xios_orchidee_send_field("nbp",(SUM((gpp_daily+co2_to_bm-(resp_maint+resp_growth+resp_hetero)-co2_fire) * &
1967	veget_cov_max,dim=2)-cflux_prod_total-cflux_prod_harvest_total-SUM(harvest_above,DIM=2))/1e3/one_day)
1968	CALL xios_orchidee_send_field("fVegLitter",SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))*&
1969	veget_cov_max,dim=2)/1e3/one_day)
1970	CALL xios_orchidee_send_field("fLitterSoil",SUM(SUM(soilcarbon_input,dim=2)*veget_cov_max,dim=2)/1e3/one_day)
1971
1972	! Carbon stock transformed from gC/m2 into kgC/m2
1973	CALL xios_orchidee_send_field("cLeaf",SUM(biomass(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3)
1974	CALL xios_orchidee_send_field("cWood",SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon))*&
1975	veget_cov_max,dim=2)/1e3)
1976	CALL xios_orchidee_send_field("cRoot",SUM(( biomass(:,:,iroot,icarbon) + biomass(:,:,isapbelow,icarbon) + &
1977	biomass(:,:,iheartbelow,icarbon) )*veget_cov_max,dim=2)/1e3)
1978	CALL xios_orchidee_send_field("cMisc",SUM(( biomass(:,:,icarbres,icarbon) + biomass(:,:,ifruit,icarbon))*&
1979	veget_cov_max,dim=2)/1e3)
1980	CALL xios_orchidee_send_field("cLitterAbove",SUM((litter(:,istructural,:,iabove,icarbon)+&
1981	litter(:,imetabolic,:,iabove,icarbon))*veget_cov_max,dim=2)/1e3)
1982	CALL xios_orchidee_send_field("cLitterBelow",SUM((litter(:,istructural,:,ibelow,icarbon)+&
1983	litter(:,imetabolic,:,ibelow,icarbon))*veget_cov_max,dim=2)/1e3)
1984	CALL xios_orchidee_send_field("cSoilFast",SUM(carbon(:,iactive,:)*veget_cov_max,dim=2)/1e3)
1985	CALL xios_orchidee_send_field("cSoilMedium",SUM(carbon(:,islow,:)*veget_cov_max,dim=2)/1e3)
1986	CALL xios_orchidee_send_field("cSoilSlow",SUM(carbon(:,ipassive,:)*veget_cov_max,dim=2)/1e3)
1987
1988	! Vegetation fractions [0,100]
1989	CALL xios_orchidee_send_field("landCoverFrac",veget_cov_max*100)
1990	vartmp(:)=zero
1991	DO j = 2,nvm
1992	IF (is_deciduous(j)) THEN
1993	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
1994	ENDIF
1995	ENDDO
1996	CALL xios_orchidee_send_field("treeFracPrimDec",vartmp)
1997	vartmp(:)=zero
1998	DO j = 2,nvm
1999	IF (is_evergreen(j)) THEN
2000	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2001	ENDIF
2002	ENDDO
2003	CALL xios_orchidee_send_field("treeFracPrimEver",vartmp)
2004	vartmp(:)=zero
2005	DO j = 2,nvm
2006	IF ( .NOT.(is_c4(j)) ) THEN
2007	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2008	ENDIF
2009	ENDDO
2010	CALL xios_orchidee_send_field("c3PftFrac",vartmp)
2011	vartmp(:)=zero
2012	DO j = 2,nvm
2013	IF ( is_c4(j) ) THEN
2014	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2015	ENDIF
2016	ENDDO
2017	CALL xios_orchidee_send_field("c4PftFrac",vartmp)
2018
2019	! Carbon fluxes transformed from gC/m2/d into kgC/m2/s
2020	CALL xios_orchidee_send_field("rGrowth",SUM(resp_growth*veget_cov_max,dim=2)/1e3/one_day)
2021	CALL xios_orchidee_send_field("rMaint",SUM(resp_maint*veget_cov_max,dim=2)/1e3/one_day)
2022	CALL xios_orchidee_send_field("nppLeaf",SUM(bm_alloc(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3/one_day)
2023	CALL xios_orchidee_send_field("nppWood",SUM(bm_alloc(:,:,isapabove,icarbon)*veget_cov_max,dim=2)/1e3/one_day)
2024	CALL xios_orchidee_send_field("nppRoot",SUM(( bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iroot,icarbon) ) * &
2025	veget_cov_max,dim=2)/1e3/one_day)
2026
2027
2028	CALL histwrite_p (hist_id_stomate, 'RESOLUTION_X', itime, &
2029	resolution(:,1), npts, hori_index)
2030	CALL histwrite_p (hist_id_stomate, 'RESOLUTION_Y', itime, &
2031	resolution(:,2), npts, hori_index)
2032	CALL histwrite_p (hist_id_stomate, 'CONTFRAC', itime, &
2033	contfrac(:), npts, hori_index)
2034
2035	CALL histwrite_p (hist_id_stomate, 'LITTER_STR_AB', itime, &
2036	litter(:,istructural,:,iabove,icarbon), npts*nvm, horipft_index)
2037	CALL histwrite_p (hist_id_stomate, 'LITTER_MET_AB', itime, &
2038	litter(:,imetabolic,:,iabove,icarbon), npts*nvm, horipft_index)
2039	CALL histwrite_p (hist_id_stomate, 'LITTER_STR_BE', itime, &
2040	litter(:,istructural,:,ibelow,icarbon), npts*nvm, horipft_index)
2041	CALL histwrite_p (hist_id_stomate, 'LITTER_MET_BE', itime, &
2042	litter(:,imetabolic,:,ibelow,icarbon), npts*nvm, horipft_index)
2043	!spitfiretest
2044	!CALL histwrite (hist_id_stomate, 'fuel_1hr_met_b', itime, &
2045	! fuel_1hr(:,:,imetabolic,icarbon), npts*nvm, horipft_index)
2046	!CALL histwrite (hist_id_stomate, 'fuel_1hr_str_b', itime, &
2047	! fuel_1hr(:,:,istructural,icarbon), npts*nvm, horipft_index)
2048	!CALL histwrite (hist_id_stomate, 'fuel_10hr_met_b', itime, &
2049	! fuel_10hr(:,:,imetabolic,icarbon), npts*nvm, horipft_index)
2050	!CALL histwrite (hist_id_stomate, 'fuel_10hr_str_b', itime, &
2051	! fuel_10hr(:,:,istructural,icarbon), npts*nvm, horipft_index)
2052	!CALL histwrite (hist_id_stomate, 'fuel_100hr_met_b', itime, &
2053	! fuel_100hr(:,:,imetabolic,icarbon), npts*nvm, horipft_index)
2054	!CALL histwrite (hist_id_stomate, 'fuel_100hr_str_b', itime, &
2055	! fuel_100hr(:,:,istructural,icarbon), npts*nvm, horipft_index)
2056	!CALL histwrite (hist_id_stomate, 'fuel_1000hr_met_b', itime, &
2057	! fuel_1000hr(:,:,imetabolic,icarbon), npts*nvm, horipft_index)
2058	!CALL histwrite (hist_id_stomate, 'fuel_1000hr_str_b', itime, &
2059	! fuel_1000hr(:,:,istructural,icarbon), npts*nvm, horipft_index)
2060	!endspittest
2061
2062	CALL histwrite_p (hist_id_stomate, 'DEADLEAF_COVER', itime, &
2063	deadleaf_cover, npts, hori_index)
2064
2065	CALL histwrite_p (hist_id_stomate, 'TOTAL_SOIL_CARB', itime, &
2066	tot_litter_soil_carb, npts*nvm, horipft_index)
2067	CALL histwrite_p (hist_id_stomate, 'CARBON_ACTIVE', itime, &
2068	carbon(:,iactive,:), npts*nvm, horipft_index)
2069	CALL histwrite_p (hist_id_stomate, 'CARBON_SLOW', itime, &
2070	carbon(:,islow,:), npts*nvm, horipft_index)
2071	CALL histwrite_p (hist_id_stomate, 'CARBON_PASSIVE', itime, &
2072	carbon(:,ipassive,:), npts*nvm, horipft_index)
2073
2074	CALL histwrite_p (hist_id_stomate, 'CARBON_ACTIVE_SURF', itime, &
2075	carbon_surf(:,iactive,:), npts*nvm, horipft_index)
2076	CALL histwrite_p (hist_id_stomate, 'CARBON_SLOW_SURF', itime, &
2077	carbon_surf(:,islow,:), npts*nvm, horipft_index)
2078	CALL histwrite_p (hist_id_stomate, 'CARBON_PASSIVE_SURF', itime, &
2079	carbon_surf(:,ipassive,:), npts*nvm, horipft_index)
2080
2081
2082	CALL histwrite_p (hist_id_stomate, 'TSURF_YEAR', itime, &
2083	tsurf_year, npts, hori_index)
2084	!pss:-
2085
2086
2087	CALL histwrite_p (hist_id_stomate, 'T2M_MONTH', itime, &
2088	t2m_month, npts, hori_index)
2089	CALL histwrite_p (hist_id_stomate, 'T2M_WEEK', itime, &
2090	t2m_week, npts, hori_index)
2091	CALL histwrite_p (hist_id_stomate, 'TSEASON', itime, &
2092	Tseason, npts, hori_index)
2093	CALL histwrite_p (hist_id_stomate, 'TMIN_SPRING_TIME', itime, &
2094	Tmin_spring_time, npts*nvm, horipft_index)
2095	CALL histwrite_p (hist_id_stomate, 'ONSET_DATE', itime, &
2096	onset_date(:,:), npts*nvm, horipft_index)
2097
2098	CALL histwrite_p (hist_id_stomate, 'HET_RESP', itime, &
2099	resp_hetero(:,:), npts*nvm, horipft_index)
2100	! gmjc
2101	CALL histwrite_p(hist_id_stomate ,'T2M_14' ,itime, &
2102	t2m_14, npts, hori_index)
2103	! CALL histwrite (hist_id_stomate, 'LITTER_RESP', itime, &
2104	! resp_hetero_litter_d(:,:), npts*nvm, horipft_index)
2105	! CALL histwrite (hist_id_stomate, 'ACTIVE_RESP', itime, &
2106	! resp_hetero_soil_d(:,iactive,:), npts*nvm, horipft_index)
2107	! CALL histwrite (hist_id_stomate, 'SLOW_RESP', itime, &
2108	! resp_hetero_soil_d(:,islow,:), npts*nvm, horipft_index)
2109	! CALL histwrite (hist_id_stomate, 'PASSIVE_RESP', itime, &
2110	! resp_hetero_soil_d(:,ipassive,:), npts*nvm, horipft_index)
2111	CALL histwrite_p (hist_id_stomate, 'LITTER_STR_AVAIL', itime, &
2112	litter_avail(:,istructural,:), npts*nvm, horipft_index)
2113	CALL histwrite_p (hist_id_stomate, 'LITTER_MET_AVAIL', itime, &
2114	litter_avail(:,imetabolic,:), npts*nvm, horipft_index)
2115	CALL histwrite_p (hist_id_stomate, 'LITTER_STR_NAVAIL', itime, &
2116	litter_not_avail(:,istructural,:), npts*nvm, horipft_index)
2117	CALL histwrite_p (hist_id_stomate, 'LITTER_MET_NAVAIL', itime, &
2118	litter_not_avail(:,imetabolic,:), npts*nvm, horipft_index)
2119	CALL histwrite_p (hist_id_stomate, 'LITTER_STR_AVAILF', itime, &
2120	litter_avail_frac(:,istructural,:), npts*nvm, horipft_index)
2121	CALL histwrite_p (hist_id_stomate, 'LITTER_MET_AVAILF', itime, &
2122	litter_avail_frac(:,imetabolic,:), npts*nvm, horipft_index)
2123
2124	! Crop is enabled
2125	IF (ANY(ok_LAIdev)) THEN
2126	CALL histwrite_p (hist_id_stomate, 'N_LIMFERT', itime, &
2127	N_limfert, npts*nvm, horipft_index)
2128	ENDIF
2129	! end gmjc
2130	CALL histwrite_p (hist_id_stomate, 'FIREINDEX', itime, &
2131	fireindex(:,:), npts*nvm, horipft_index)
2132	CALL histwrite_p (hist_id_stomate, 'LITTERHUM', itime, &
2133	litterhum_daily, npts, hori_index)
2134	CALL histwrite_p (hist_id_stomate, 'CO2_FIRE', itime, &
2135	co2_fire, npts*nvm, horipft_index)
2136	CALL histwrite_p (hist_id_stomate, 'CO2_TAKEN', itime, &
2137	co2_to_bm, npts*nvm, horipft_index)
2138	! land cover change
2139	CALL histwrite_p (hist_id_stomate, 'CONVFLUX_HAR', itime, &
2140	convflux(:,iwphar), npts, hori_index)
2141	CALL histwrite_p (hist_id_stomate, 'CONVFLUX_LCC', itime, &
2142	convflux(:,iwplcc), npts, hori_index)
2143	CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD10_HAR', itime, &
2144	cflux_prod10(:,iwphar), npts, hori_index)
2145	CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD10_LCC', itime, &
2146	cflux_prod10(:,iwplcc), npts, hori_index)
2147	CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD100_HAR', itime, &
2148	cflux_prod100(:,iwphar), npts, hori_index)
2149	CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD100_LCC', itime, &
2150	cflux_prod100(:,iwplcc), npts, hori_index)
2151
2152	CALL histwrite_p (hist_id_stomate, 'HARVEST_ABOVE', itime, &
2153	harvest_above, npts, horipft_index)
2154	CALL histwrite_p (hist_id_stomate, 'PROD10_HAR', itime, &
2155	prod10(:,:,iwphar), npts*11, horip11_index)
2156	CALL histwrite_p (hist_id_stomate, 'PROD10_LCC', itime, &
2157	prod10(:,:,iwplcc), npts*11, horip11_index)
2158	CALL histwrite_p (hist_id_stomate, 'PROD100_HAR', itime, &
2159	prod100(:,:,iwphar), npts*101, horip101_index)
2160	CALL histwrite_p (hist_id_stomate, 'PROD100_LCC', itime, &
2161	prod100(:,:,iwplcc), npts*101, horip101_index)
2162	CALL histwrite_p (hist_id_stomate, 'FLUX10_HAR', itime, &
2163	flux10(:,:,iwphar), npts*10, horip10_index)
2164	CALL histwrite_p (hist_id_stomate, 'FLUX10_LCC', itime, &
2165	flux10(:,:,iwplcc), npts*10, horip10_index)
2166	CALL histwrite_p (hist_id_stomate, 'FLUX100_HAR', itime, &
2167	flux100(:,:,iwphar), npts*100, horip100_index)
2168	CALL histwrite_p (hist_id_stomate, 'FLUX100_LCC', itime, &
2169	flux100(:,:,iwplcc), npts*100, horip100_index)
2170	CALL histwrite_p (hist_id_stomate, 'DefLitSurplus', itime, &
2171	deflitsup_total, npts*nvm, horipft_index)
2172	CALL histwrite_p (hist_id_stomate, 'DefBioSurplus', itime, &
2173	defbiosup_total, npts*nvm, horipft_index)
2174
2175	IF (use_bound_spa) THEN
2176	CALL histwrite_p (hist_id_stomate, 'bound_spa', itime, &
2177	bound_spa, npts*nvm, horipft_index)
2178	ENDIF
2179
2180	IF (do_now_stomate_lcchange) THEN
2181	CALL histwrite_p (hist_id_stomate, 'LCC', itime, &
2182	lcc, npts*nvm, horipft_index)
2183	ENDIF
2184
2185	CALL histwrite_p (hist_id_stomate, 'LAI', itime, &
2186	lai, npts*nvm, horipft_index)
2187	CALL histwrite_p (hist_id_stomate, 'FPC_MAX', itime, &
2188	fpc_max, npts*nvm, horipft_index)
2189	CALL histwrite_p (hist_id_stomate, 'MAXFPC_LASTYEAR', itime, &
2190	maxfpc_lastyear, npts*nvm, horipft_index)
2191	CALL histwrite_p (hist_id_stomate, 'VEGET_COV_MAX', itime, &
2192	veget_cov_max, npts*nvm, horipft_index)
2193	CALL histwrite_p (hist_id_stomate, 'NPP', itime, &
2194	npp_daily, npts*nvm, horipft_index)
2195	CALL histwrite_p (hist_id_stomate, 'GPP', itime, &
2196	gpp_daily, npts*nvm, horipft_index)
2197	CALL histwrite_p (hist_id_stomate, 'IND', itime, &
2198	ind, npts*nvm, horipft_index)
2199	CALL histwrite_p (hist_id_stomate, 'CN_IND', itime, &
2200	cn_ind, npts*nvm, horipft_index)
2201	CALL histwrite_p (hist_id_stomate, 'WOODMASS_IND', itime, &
2202	woodmass_ind, npts*nvm, horipft_index)
2203	CALL histwrite_p (hist_id_stomate, 'TOTAL_M', itime, &
2204	tot_live_biomass(:,:,icarbon), npts*nvm, horipft_index)
2205	CALL histwrite_p (hist_id_stomate, 'LEAF_M', itime, &
2206	biomass(:,:,ileaf,icarbon), npts*nvm, horipft_index)
2207	CALL histwrite_p (hist_id_stomate, 'SAP_M_AB', itime, &
2208	biomass(:,:,isapabove,icarbon), npts*nvm, horipft_index)
2209	CALL histwrite_p (hist_id_stomate, 'SAP_M_BE', itime, &
2210	biomass(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
2211	CALL histwrite_p (hist_id_stomate, 'HEART_M_AB', itime, &
2212	biomass(:,:,iheartabove,icarbon), npts*nvm, horipft_index)
2213	CALL histwrite_p (hist_id_stomate, 'HEART_M_BE', itime, &
2214	biomass(:,:,iheartbelow,icarbon), npts*nvm, horipft_index)
2215	CALL histwrite_p (hist_id_stomate, 'ROOT_M', itime, &
2216	biomass(:,:,iroot,icarbon), npts*nvm, horipft_index)
2217	CALL histwrite_p (hist_id_stomate, 'FRUIT_M', itime, &
2218	biomass(:,:,ifruit,icarbon), npts*nvm, horipft_index)
2219	!!!!! crop variables
2220	CALL histwrite_p (hist_id_stomate, 'CROPYIELD', itime, &
2221	biomass(:,:,ifruit,icarbon), npts*nvm, horipft_index)
2222
2223	CALL histwrite_p (hist_id_stomate, 'BIOMYIELD', itime, &
2224	biomass(:,:,ileaf,icarbon)+biomass(:,:,isapabove,icarbon) &
2225	+biomass(:,:,ifruit,icarbon)+biomass(:,:,icarbres,icarbon), npts*nvm,horipft_index)
2226
2227	CALL histwrite_p (hist_id_stomate, 'CROP_EXPORT', itime, &
2228	crop_export, npts*nvm, horipft_index)
2229	!!!!! end crop variables, xuhui
2230	CALL histwrite_p (hist_id_stomate, 'RESERVE_M', itime, &
2231	biomass(:,:,icarbres,icarbon), npts*nvm, horipft_index)
2232	CALL histwrite_p (hist_id_stomate, 'TOTAL_TURN', itime, &
2233	tot_turnover(:,:,icarbon), npts*nvm, horipft_index)
2234	CALL histwrite_p (hist_id_stomate, 'LEAF_TURN', itime, &
2235	turnover_daily(:,:,ileaf,icarbon), npts*nvm, horipft_index)
2236	CALL histwrite_p (hist_id_stomate, 'SAP_AB_TURN', itime, &
2237	turnover_daily(:,:,isapabove,icarbon), npts*nvm, horipft_index)
2238	CALL histwrite_p (hist_id_stomate, 'ROOT_TURN', itime, &
2239	turnover_daily(:,:,iroot,icarbon), npts*nvm, horipft_index)
2240	CALL histwrite_p (hist_id_stomate, 'FRUIT_TURN', itime, &
2241	turnover_daily(:,:,ifruit,icarbon), npts*nvm, horipft_index)
2242	CALL histwrite_p (hist_id_stomate, 'TOTAL_BM_LITTER', itime, &
2243	tot_bm_to_litter(:,:,icarbon), npts*nvm, horipft_index)
2244	CALL histwrite_p (hist_id_stomate, 'LEAF_BM_LITTER', itime, &
2245	bm_to_litter(:,:,ileaf,icarbon), npts*nvm, horipft_index)
2246	CALL histwrite_p (hist_id_stomate, 'SAP_AB_BM_LITTER', itime, &
2247	bm_to_litter(:,:,isapabove,icarbon), npts*nvm, horipft_index)
2248	CALL histwrite_p (hist_id_stomate, 'SAP_BE_BM_LITTER', itime, &
2249	bm_to_litter(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
2250	CALL histwrite_p (hist_id_stomate, 'HEART_AB_BM_LITTER', itime, &
2251	bm_to_litter(:,:,iheartabove,icarbon), npts*nvm, horipft_index)
2252	CALL histwrite_p (hist_id_stomate, 'HEART_BE_BM_LITTER', itime, &
2253	bm_to_litter(:,:,iheartbelow,icarbon), npts*nvm, horipft_index)
2254	CALL histwrite_p (hist_id_stomate, 'ROOT_BM_LITTER', itime, &
2255	bm_to_litter(:,:,iroot,icarbon), npts*nvm, horipft_index)
2256	CALL histwrite_p (hist_id_stomate, 'FRUIT_BM_LITTER', itime, &
2257	bm_to_litter(:,:,ifruit,icarbon), npts*nvm, horipft_index)
2258	CALL histwrite_p (hist_id_stomate, 'RESERVE_BM_LITTER', itime, &
2259	bm_to_litter(:,:,icarbres,icarbon), npts*nvm, horipft_index)
2260	CALL histwrite_p (hist_id_stomate, 'MAINT_RESP', itime, &
2261	resp_maint, npts*nvm, horipft_index)
2262	CALL histwrite_p (hist_id_stomate, 'GROWTH_RESP', itime, &
2263	resp_growth, npts*nvm, horipft_index)
2264	CALL histwrite_p (hist_id_stomate, 'AGE', itime, &
2265	age, npts*nvm, horipft_index)
2266	CALL histwrite_p (hist_id_stomate, 'HEIGHT', itime, &
2267	height, npts*nvm, horipft_index)
2268	CALL histwrite_p (hist_id_stomate, 'MOISTRESS', itime, &
2269	moiavail_week, npts*nvm, horipft_index)
2270	CALL histwrite_p (hist_id_stomate, 'VCMAX', itime, &
2271	vcmax, npts*nvm, horipft_index)
2272	CALL histwrite_p (hist_id_stomate, 'TURNOVER_TIME', itime, &
2273	turnover_time, npts*nvm, horipft_index)
2274	!!DZADD
2275	! CALL histwrite_p (hist_id_stomate, 'LEAF_FRAC1', itime, leaf_frac(:,:,1), npts*nvm, horipft_index)
2276	! CALL histwrite_p (hist_id_stomate, 'LEAF_FRAC2', itime, leaf_frac(:,:,2), npts*nvm, horipft_index)
2277	! CALL histwrite_p (hist_id_stomate, 'LEAF_FRAC3', itime, leaf_frac(:,:,3), npts*nvm, horipft_index)
2278	! CALL histwrite_p (hist_id_stomate, 'LEAF_FRAC4', itime, leaf_frac(:,:,4), npts*nvm, horipft_index)
2279	!!ENDDZADD
2280
2281	IF ( hist_id_stomate_IPCC > 0 ) THEN
2282	vartmp(:)=SUM(tot_live_biomass(:,:,icarbon)*veget_cov_max,dim=2)/1e3
2283	CALL histwrite_p (hist_id_stomate_IPCC, "cVeg", itime, &
2284	vartmp, npts, hori_index)
2285	vartmp(:)=SUM(tot_litter_carb*veget_cov_max,dim=2)/1e3
2286	CALL histwrite_p (hist_id_stomate_IPCC, "cLitter", itime, &
2287	vartmp, npts, hori_index)
2288	vartmp(:)=SUM(tot_soil_carb*veget_cov_max,dim=2)/1e3
2289	CALL histwrite_p (hist_id_stomate_IPCC, "cSoil", itime, &
2290	vartmp, npts, hori_index)
2291	vartmp(:)=(prod10_total + prod100_total + prod10_harvest_total + prod100_harvest_total)/1e3
2292	CALL histwrite_p (hist_id_stomate_IPCC, "cProduct", itime, &
2293	vartmp, npts, hori_index)
2294	vartmp(:)=carb_mass_variation/1e3/one_day
2295	CALL histwrite_p (hist_id_stomate_IPCC, "cMassVariation", itime, &
2296	vartmp, npts, hori_index)
2297	vartmp(:)=SUM(lai*veget_cov_max,dim=2)
2298	CALL histwrite_p (hist_id_stomate_IPCC, "lai", itime, &
2299	vartmp, npts, hori_index)
2300	vartmp(:)=SUM(gpp_daily*veget_cov_max,dim=2)/1e3/one_day
2301	CALL histwrite_p (hist_id_stomate_IPCC, "gpp", itime, &
2302	vartmp, npts, hori_index)
2303	vartmp(:)=SUM((resp_maint+resp_growth)*veget_cov_max,dim=2)/1e3/one_day
2304	CALL histwrite_p (hist_id_stomate_IPCC, "ra", itime, &
2305	vartmp, npts, hori_index)
2306	vartmp(:)=SUM(npp_daily*veget_cov_max,dim=2)/1e3/one_day
2307	CALL histwrite_p (hist_id_stomate_IPCC, "npp", itime, &
2308	vartmp, npts, hori_index)
2309	vartmp(:)=SUM(resp_hetero*veget_cov_max,dim=2)/1e3/one_day
2310	CALL histwrite_p (hist_id_stomate_IPCC, "rh", itime, &
2311	vartmp, npts, hori_index)
2312	vartmp(:)=SUM(co2_fire*veget_cov_max,dim=2)/1e3/one_day
2313	CALL histwrite_p (hist_id_stomate_IPCC, "fFire", itime, &
2314	vartmp, npts, hori_index)
2315	vartmp(:)=SUM(harvest_above,DIM=2)/1e3/one_day
2316	CALL histwrite_p (hist_id_stomate_IPCC, "fHarvest", itime, &
2317	vartmp, npts, hori_index)
2318	vartmp(:)=cflux_prod_total/1e3/one_day
2319	CALL histwrite_p (hist_id_stomate_IPCC, "fLuc", itime, &
2320	vartmp, npts, hori_index)
2321	!gmjc
2322	IF (enable_grazing) THEN
2323	vartmp(:)=(SUM((gpp_daily+co2_to_bm-(resp_maint+resp_growth+resp_hetero)-co2_fire &
2324	-harvest_gm-ranimal_gm-ch4_pft_gm+cinput_gm) & ! specific
2325	& *veget_cov_max,dim=2)-cflux_prod_total-cflux_prod_harvest_total-SUM(harvest_above,DIM=2))/1e3/one_day
2326	ELSE
2327	vartmp(:)=(SUM((gpp_daily+co2_to_bm-(resp_maint+resp_growth+resp_hetero)-co2_fire) &
2328	& *veget_cov_max,dim=2)-cflux_prod_total-cflux_prod_harvest_total-SUM(harvest_above,DIM=2))/1e3/one_day
2329	ENDIF
2330
2331	CALL histwrite_p (hist_id_stomate_IPCC, "fWoodharvest", itime, &
2332	vartmp, npts, hori_index)
2333	!end gmjc
2334	! vartmp(:)=(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire) &
2335	! & *veget_cov_max,dim=2)-cflux_prod_total-harvest_above)/1e3/one_day
2336	CALL histwrite_p (hist_id_stomate_IPCC, "nbp", itime, &
2337	vartmp, npts, hori_index)
2338	vartmp(:)=SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))*veget_cov_max,dim=2)/1e3/one_day
2339	CALL histwrite_p (hist_id_stomate_IPCC, "fVegLitter", itime, &
2340	vartmp, npts, hori_index)
2341	vartmp(:)=SUM(SUM(soilcarbon_input,dim=2)*veget_cov_max,dim=2)/1e3/one_day
2342	CALL histwrite_p (hist_id_stomate_IPCC, "fLitterSoil", itime, &
2343	vartmp, npts, hori_index)
2344	vartmp(:)=SUM(biomass(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3
2345	CALL histwrite_p (hist_id_stomate_IPCC, "cLeaf", itime, &
2346	vartmp, npts, hori_index)
2347	vartmp(:)=SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon))*veget_cov_max,dim=2)/1e3
2348	CALL histwrite_p (hist_id_stomate_IPCC, "cWood", itime, &
2349	vartmp, npts, hori_index)
2350	vartmp(:)=SUM(( biomass(:,:,iroot,icarbon) + biomass(:,:,isapbelow,icarbon) + biomass(:,:,iheartbelow,icarbon) ) &
2351	& *veget_cov_max,dim=2)/1e3
2352	CALL histwrite_p (hist_id_stomate_IPCC, "cRoot", itime, &
2353	vartmp, npts, hori_index)
2354	vartmp(:)=SUM(( biomass(:,:,icarbres,icarbon) + biomass(:,:,ifruit,icarbon))*veget_cov_max,dim=2)/1e3
2355	CALL histwrite_p (hist_id_stomate_IPCC, "cMisc", itime, &
2356	vartmp, npts, hori_index)
2357	vartmp(:)=SUM((litter(:,istructural,:,iabove,icarbon)+litter(:,imetabolic,:,iabove,icarbon))*&
2358	veget_cov_max,dim=2)/1e3
2359	CALL histwrite_p (hist_id_stomate_IPCC, "cLitterAbove", itime, &
2360	vartmp, npts, hori_index)
2361	vartmp(:)=SUM((litter(:,istructural,:,ibelow,icarbon)+litter(:,imetabolic,:,ibelow,icarbon))*&
2362	veget_cov_max,dim=2)/1e3
2363	CALL histwrite_p (hist_id_stomate_IPCC, "cLitterBelow", itime, &
2364	vartmp, npts, hori_index)
2365	vartmp(:)=SUM(carbon(:,iactive,:)*veget_cov_max,dim=2)/1e3
2366	CALL histwrite_p (hist_id_stomate_IPCC, "cSoilFast", itime, &
2367	vartmp, npts, hori_index)
2368	vartmp(:)=SUM(carbon(:,islow,:)*veget_cov_max,dim=2)/1e3
2369	CALL histwrite_p (hist_id_stomate_IPCC, "cSoilMedium", itime, &
2370	vartmp, npts, hori_index)
2371	vartmp(:)=SUM(carbon(:,ipassive,:)*veget_cov_max,dim=2)/1e3
2372	CALL histwrite_p (hist_id_stomate_IPCC, "cSoilSlow", itime, &
2373	vartmp, npts, hori_index)
2374	DO j=1,nvm
2375	histvar(:,j)=veget_cov_max(:,j)*100
2376	ENDDO
2377	CALL histwrite_p (hist_id_stomate_IPCC, "landCoverFrac", itime, &
2378	histvar, npts*nvm, horipft_index)
2379	!-
2380	vartmp(:)=zero
2381	DO j = 2,nvm
2382	IF (is_deciduous(j)) THEN
2383	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2384	ENDIF
2385	ENDDO
2386	CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimDec", itime, &
2387	vartmp, npts, hori_index)
2388	!-
2389	vartmp(:)=zero
2390	DO j = 2,nvm
2391	IF (is_evergreen(j)) THEN
2392	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2393	ENDIF
2394	ENDDO
2395	CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimEver", itime, &
2396	vartmp, npts, hori_index)
2397	!-
2398	vartmp(:)=zero
2399	DO j = 2,nvm
2400	IF ( .NOT.(is_c4(j)) ) THEN
2401	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2402	ENDIF
2403	ENDDO
2404	CALL histwrite_p (hist_id_stomate_IPCC, "c3PftFrac", itime, &
2405	vartmp, npts, hori_index)
2406	!-
2407	vartmp(:)=zero
2408	DO j = 2,nvm
2409	IF ( is_c4(j) ) THEN
2410	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2411	ENDIF
2412	ENDDO
2413	CALL histwrite_p (hist_id_stomate_IPCC, "c4PftFrac", itime, &
2414	vartmp, npts, hori_index)
2415	!-
2416	vartmp(:)=SUM(resp_growth*veget_cov_max,dim=2)/1e3/one_day
2417	CALL histwrite_p (hist_id_stomate_IPCC, "rGrowth", itime, &
2418	vartmp, npts, hori_index)
2419	vartmp(:)=SUM(resp_maint*veget_cov_max,dim=2)/1e3/one_day
2420	CALL histwrite_p (hist_id_stomate_IPCC, "rMaint", itime, &
2421	vartmp, npts, hori_index)
2422	vartmp(:)=SUM(bm_alloc(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3/one_day
2423	CALL histwrite_p (hist_id_stomate_IPCC, "nppLeaf", itime, &
2424	vartmp, npts, hori_index)
2425	vartmp(:)=SUM(bm_alloc(:,:,isapabove,icarbon)*veget_cov_max,dim=2)/1e3/one_day
2426	CALL histwrite_p (hist_id_stomate_IPCC, "nppWood", itime, &
2427	vartmp, npts, hori_index)
2428	vartmp(:)=SUM(( bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iroot,icarbon) )*veget_cov_max,dim=2)/1e3/one_day
2429	CALL histwrite_p (hist_id_stomate_IPCC, "nppRoot", itime, &
2430	vartmp, npts, hori_index)
2431
2432	CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_X', itime, &
2433	resolution(:,1), npts, hori_index)
2434	CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_Y', itime, &
2435	resolution(:,2), npts, hori_index)
2436	CALL histwrite_p (hist_id_stomate_IPCC, 'CONTFRAC', itime, &
2437	contfrac(:), npts, hori_index)
2438
2439	ENDIF
2440
2441	IF (printlev>=4) WRITE(numout,*) 'Leaving stomate_lpj'
2442
2443	END SUBROUTINE StomateLpj
2444
2445
2446	!! ================================================================================================================================
2447	!! SUBROUTINE : harvest
2448	!!
2449	!>\BRIEF Harvest of croplands
2450	!!
2451	!! DESCRIPTION : To take into account biomass harvest from crop (mainly to take
2452	!! into account for the reduced litter input and then decreased soil carbon. it is a
2453	!! constant (40\%) fraction of above ground biomass.
2454	!!
2455	!! RECENT CHANGE(S) : None
2456	!!
2457	!! MAIN OUTPUT VARIABLE(S): ::harvest_above the harvested biomass
2458	!!
2459	!! REFERENCE(S) :
2460	!! - Piao, S., P. Ciais, P. Friedlingstein, N. de Noblet-Ducoudre, P. Cadule, N. Viovy, and T. Wang. 2009.
2461	!! Spatiotemporal patterns of terrestrial carbon cycle during the 20th century. Global Biogeochemical
2462	!! Cycles 23:doi:10.1029/2008GB003339.
2463	!!
2464	!! FLOWCHART : None
2465	!! \n
2466	!_ ================================================================================================================================
2467
2468	SUBROUTINE harvest(npts, dt_days, veget_cov_max, &
2469	bm_to_litter, turnover_daily, &
2470	harvest_above)
2471
2472	!! 0. Variable and parameter declaration
2473
2474	!! 0.1 Input variables
2475
2476	INTEGER, INTENT(in) :: npts !! Domain size (unitless)
2477	REAL(r_std), INTENT(in) :: dt_days !! Time step (days)
2478	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_cov_max !! new "maximal" coverage fraction of a PFT (LAI ->
2479	!! infinity) on ground @tex $(m^2 m^{-2})$ @endtex
2480
2481	!! 0.2 Output variables
2482
2483	!! 0.3 Modified variables
2484
2485	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: bm_to_litter !! [DISPENSABLE] conversion of biomass to litter
2486	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
2487	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: turnover_daily !! Turnover rates
2488	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
2489	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: harvest_above !! harvest above ground biomass for agriculture
2490	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
2491	!! 0.4 Local variables
2492
2493	INTEGER(i_std) :: i, j, k, l, m !! indices
2494	REAL(r_std) :: above_old !! biomass of previous time step
2495	!! @tex $(gC m^{-2})$ @endtex
2496	!_ ================================================================================================================================
2497
2498	!! 1. Yearly initialisation
2499
2500	above_old = zero
2501	harvest_above = zero
2502
2503	DO i = 1, npts
2504	DO j = 1,nvm
2505	IF (.NOT. natural(j)) THEN
2506	above_old = turnover_daily(i,j,ileaf,icarbon) + turnover_daily(i,j,isapabove,icarbon) + &
2507	& turnover_daily(i,j,iheartabove,icarbon) + turnover_daily(i,j,ifruit,icarbon) + &
2508	& turnover_daily(i,j,icarbres,icarbon) + turnover_daily(i,j,isapbelow,icarbon) + &
2509	& turnover_daily(i,j,iheartbelow,icarbon) + turnover_daily(i,j,iroot,icarbon)
2510
2511	turnover_daily(i,j,ileaf,icarbon) = turnover_daily(i,j,ileaf,icarbon)*frac_turnover_daily
2512	turnover_daily(i,j,isapabove,icarbon) = turnover_daily(i,j,isapabove,icarbon)*frac_turnover_daily
2513	turnover_daily(i,j,isapbelow,icarbon) = turnover_daily(i,j,isapbelow,icarbon)*frac_turnover_daily
2514	turnover_daily(i,j,iheartabove,icarbon) = turnover_daily(i,j,iheartabove,icarbon)*frac_turnover_daily
2515	turnover_daily(i,j,iheartbelow,icarbon) = turnover_daily(i,j,iheartbelow,icarbon)*frac_turnover_daily
2516	turnover_daily(i,j,iroot,icarbon) = turnover_daily(i,j,iroot,icarbon)*frac_turnover_daily
2517	turnover_daily(i,j,ifruit,icarbon) = turnover_daily(i,j,ifruit,icarbon)*frac_turnover_daily
2518	turnover_daily(i,j,icarbres,icarbon) = turnover_daily(i,j,icarbres,icarbon)*frac_turnover_daily
2519	harvest_above(i,j) = veget_cov_max(i,j) * above_old *(un - frac_turnover_daily)
2520	ENDIF
2521	ENDDO
2522	ENDDO
2523
2524	!!$ harvest_above = harvest_above
2525	END SUBROUTINE harvest
2526	END MODULE stomate_lpj

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

Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_stomate/stomate_lpj.f90 @ 6737

Download in other formats: