Context Navigation

stomate_lpj.f90 @ 7346

Last change on this file since 7346 was 6249, checked in by chunjing.qiu, 5 years ago
till litter and soil
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 169.8 KB

Line
1	! ================================================================================================================================
2	! MODULE : stomate_lpj
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.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_glcchange_fh
56	USE stomate_glcchange_SinAgeC_fh
57	USE stomate_lai
58	!gmjc
59	USE grassland_management
60	!end gmjc
61	! USE Orch_Write_field_p
62
63	!pss:+
64	! USE stomate_cste_wetlands
65	USE stomate_wet_ch4_pt_ter_0
66	USE stomate_wet_ch4_pt_ter_wet1
67	USE stomate_wet_ch4_pt_ter_wet2
68	USE stomate_wet_ch4_pt_ter_wet3
69	USE stomate_wet_ch4_pt_ter_wet4
70	!pss:-
71
72
73
74	IMPLICIT NONE
75
76	! private & public routines
77
78	PRIVATE
79	PUBLIC StomateLpj,StomateLpj_clear
80
81	LOGICAL, SAVE :: firstcall = .TRUE. !! first call
82	!$OMP THREADPRIVATE(firstcall)
83	!gmjc
84	! flag that enable grazing
85	LOGICAL, SAVE :: enable_grazing
86	!end gmjc
87
88	CONTAINS
89
90
91	!! ================================================================================================================================
92	!! SUBROUTINE : StomateLpj_clear
93	!!
94	!>\BRIEF Re-initialisation of variable
95	!!
96	!! DESCRIPTION : This subroutine reinitializes variables. To be used if we want to relaunch
97	!! ORCHIDEE but the routine is not used in current version.
98	!!
99	!! RECENT CHANGE(S) : None
100	!!
101	!! MAIN OUTPUT VARIABLE(S): None
102	!!
103	!! REFERENCE(S) : None
104	!!
105	!! FLOWCHART : None
106	!! \n
107	!_ ================================================================================================================================
108
109	SUBROUTINE StomateLpj_clear
110
111	CALL prescribe_clear
112	CALL phenology_clear
113	CALL npp_calc_clear
114	CALL turn_clear
115	CALL soilcarbon_clear
116	CALL constraints_clear
117	CALL establish_clear
118	CALL fire_clear
119	CALL spitfire_clear
120	CALL gap_clear
121	CALL light_clear
122	CALL pftinout_clear
123	CALL alloc_clear
124
125	CALL grassmanag_clear
126	!pss:+
127	CALL ch4_wet_flux_density_clear_0
128	CALL ch4_wet_flux_density_clear_wet1
129	CALL ch4_wet_flux_density_clear_wet2
130	CALL ch4_wet_flux_density_clear_wet3
131	CALL ch4_wet_flux_density_clear_wet4
132	!pss:-
133
134	END SUBROUTINE StomateLpj_clear
135
136
137	!! ================================================================================================================================
138	!! SUBROUTINE : StomateLPJ
139	!!
140	!>\BRIEF Main entry point for daily processes in STOMATE and LPJ, structures the call sequence
141	!! to the different processes such as dispersion, establishment, competition and mortality of PFT's.
142	!!
143	!! DESCRIPTION : This routine is the main entry point to all processes calculated on a
144	!! daily time step. Is mainly devoted to call the different STOMATE and LPJ routines
145	!! depending of the ok_dgvm (is dynamic veg used) and lpj_constant_mortality (is background mortality used).
146	!! It also prepares the cumulative
147	!! fluxes or pools (e.g TOTAL_M TOTAL_BM_LITTER etc...)
148	!!
149	!! This routine makes frequent use of "weekly", "monthly" and "long term" variables. Quotion is used because
150	!! by default "weekly" denotes 7 days, by default "monthly" denotes 20 days and by default "Long term" denotes
151	!! 3 years. dtslow refers to 24 hours (1 day).
152	!!
153	!!
154	!! RECENT CHANGE(S) : None
155	!!
156	!! MAIN OUTPUT VARIABLE(S): All variables related to stomate and required for LPJ dynamic vegetation mode.
157	!!
158	!! REFERENCE(S) :
159	!! - Krinner, G., N. Viovy, N. de Noblet-DucoudrÃ©, J. OgeÃ©, J. Polcher, P. Friedlingstein, P. Ciais, S. Sitch,
160	!! and I. C. Prentice. 2005. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere
161	!! system. Global Biogeochemical Cycles 19:GB1015, doi:1010.1029/2003GB002199.
162	!! - Sitch, S., B. Smith, I. C. Prentice, A. Arneth, A. Bondeau, W. Cramer, J. O. Kaplan, S. Levis, W. Lucht,
163	!! M. T. Sykes, K. Thonicke, and S. Venevsky. 2003. Evaluation of ecosystem dynamics, plant geography and
164	!! terrestrial carbon cycling in the LPJ dynamic global vegetation model. Global Change Biology 9:161-185.
165	!!
166	!! FLOWCHART : Update with existing flowchart from N Viovy (Jan 19, 2012)
167	!! \n
168	!_ ================================================================================================================================
169
170	SUBROUTINE StomateLpj (npts, dt_days, &
171	lalo, neighbours, resolution, contfrac, &
172	clay, herbivores, &
173	tsurf_daily, tsoil_daily, t2m_daily, t2m_min_daily, &
174	!spitfire
175	t2m_max_daily, precip_daily, wspeed_daily, lightn, popd, a_nd, &
176	read_observed_ba, observed_ba, &
177	read_cf_fine,cf_fine,read_cf_coarse,cf_coarse,read_ratio_flag,ratio_flag,read_ratio,ratio,date,&
178	!endspit
179	litterhum_daily, soilhum_daily, &
180	maxmoiavail_lastyear, minmoiavail_lastyear, &
181	gdd0_lastyear, precip_lastyear, &
182	moiavail_month, moiavail_week, t2m_longterm, t2m_month, t2m_week, &
183	tsoil_month, soilhum_month, &
184	gdd_m5_dormance, gdd_from_growthinit, gdd_midwinter, ncd_dormance, ngd_minus5, &
185	turnover_longterm, gpp_daily, gpp_week, &
186	time_hum_min, maxfpc_lastyear, resp_maint_part, &
187	PFTpresent, age, fireindex, firelitter, &
188	leaf_age, leaf_frac, biomass, ind, adapted, regenerate, &
189	senescence, when_growthinit, &
190	litterpart, litter, litter_avail, litter_not_avail, litter_avail_frac, &
191	dead_leaves, carbon,carbon_surf, lignin_struc, &
192	!spitfire
193	ni_acc,fire_numday,fuel_1hr,fuel_10hr,fuel_100hr,fuel_1000hr, &
194	lcc,bafrac_deforest_accu,emideforest_litter_accu,emideforest_biomass_accu,&
195	deforest_litter_remain,deforest_biomass_remain,&
196	def_fuel_1hr_remain,def_fuel_10hr_remain,&
197	def_fuel_100hr_remain,def_fuel_1000hr_remain,&
198	!endspit
199	veget_max, veget_max_new, npp_longterm, lm_lastyearmax, &
200	veget_lastlight, everywhere, need_adjacent, RIP_time, &
201	lai, rprof,npp_daily, turnover_daily, turnover_time,&
202	control_moist, control_temp, soilcarbon_input, &
203	co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, &
204	height, deadleaf_cover, vcmax, &
205	bm_to_litter, &
206	prod10,prod100,flux10, flux100, &
207	vegetnew_firstday, glccNetLCC, &
208	glccSecondShift,glccPrimaryShift, &
209	harvest_matrix, bound_spa, glcc_pft, &
210	convflux,cflux_prod10,cflux_prod100, harvest_above, carb_mass_total, &
211	fpc_max, Tseason, Tseason_length, Tseason_tmp, &
212	Tmin_spring_time, begin_leaves, onset_date, &
213	MatrixA,&
214	!!!! crop variables
215	pdlai, slai, &
216	! for crop allocation
217	in_cycle, deltai, dltaisen, ssla, pgrain, deltgrain, reprac, &
218	nger, nlev, ndrp, nlax, &
219	c_reserve, c_leafb, nmat, nrec, N_limfert, tday_counter, &
220	!!!! end crop variables, xuhui
221	zz_coef_deep, deepC_a, deepC_s, deepC_p, & !pss:+
222	ch4_flux_density_tot_0, ch4_flux_density_dif_0, ch4_flux_density_bub_0, ch4_flux_density_pla_0,&
223	ch4_flux_density_tot_wet1,ch4_flux_density_dif_wet1,ch4_flux_density_bub_wet1,ch4_flux_density_pla_wet1,&
224	ch4_flux_density_tot_wet2,ch4_flux_density_dif_wet2,ch4_flux_density_bub_wet2,ch4_flux_density_pla_wet2,&
225	ch4_flux_density_tot_wet3,ch4_flux_density_dif_wet3,ch4_flux_density_bub_wet3,ch4_flux_density_pla_wet3,&
226	ch4_flux_density_tot_wet4,ch4_flux_density_dif_wet4,ch4_flux_density_bub_wet4,ch4_flux_density_pla_wet4,&
227	tsurf_year, &!) !pss:-
228	!gmjc
229	wshtotsum, sr_ugb, compt_ugb, nb_ani, grazed_frac, &
230	import_yield, sla_age1, t2m_14, sla_calc, snowfall_daily, day_of_year, &
231	when_growthinit_cut, nb_grazingdays, &
232	EndOfYear, &
233	moiavail_daily,tmc_topgrass_daily,fc_grazing,snowmass_daily,&
234	after_snow, after_wet, wet1day, wet2day, &
235	!end gmjc
236	!!!qcj++ peatland
237	tcarbon_acro,tcarbon_cato, carbon_acro, carbon_cato,height_acro,height_cato, &
238	resp_acro_oxic_d, resp_acro_anoxic_d,resp_cato_d,acro_to_cato_d,litter_to_acro_d,&
239	carbon_save,deepC_a_save,deepC_s_save,deepC_p_save,delta_fsave,liqwt_max_lastyear,veget_max_adjusted)
240
241	!! 0. Variable and parameter declaration
242
243	!! 0.1 input
244
245	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
246	REAL(r_std), INTENT(in) :: dt_days !! Time step of Stomate (days)
247	INTEGER(i_std), DIMENSION(npts,NbNeighb), INTENT(in) :: neighbours !! Indices of the 8 neighbours of each grid
248	!! point [1=North and then clockwise]
249	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: resolution !! Resolution at each grid point (m)
250	!! [1=E-W, 2=N-S]
251	REAL(r_std),DIMENSION(npts,2),INTENT(in) :: lalo !! Geographical coordinates (latitude,longitude)
252	!! for pixels (degrees)
253	REAL(r_std),DIMENSION (npts), INTENT (in) :: contfrac !! Fraction of continent in the grid cell (unitless)
254	REAL(r_std), DIMENSION(npts), INTENT(in) :: clay !! Clay fraction (0 to 1, unitless)
255	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: herbivores !! Time constant of probability of a leaf to
256	!! be eaten by a herbivore (days)
257	REAL(r_std), DIMENSION(npts), INTENT(in) :: tsurf_daily !! Daily surface temperatures (K)
258	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: tsoil_daily !! Daily soil temperatures (K)
259	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_daily !! Daily 2 meter temperatures (K)
260	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_min_daily !! Daily minimum 2 meter temperatures (K)
261	!spitfire
262	INTEGER(i_std),INTENT(in) :: date !! Date (days)
263	! daily maximum 2 meter temperatures (K)
264	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_max_daily
265	! daily precip (mm/day)
266	REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_daily
267	! Wind speed
268	REAL(r_std), DIMENSION(npts), INTENT(in) :: wspeed_daily
269	! Lightning flash rate
270	REAL(r_std), DIMENSION(npts), INTENT(in) :: lightn
271	! Population density rate
272	REAL(r_std), DIMENSION(npts), INTENT(inout) :: popd !popd declared and allocated and input in slowproc.f90
273	! Flag for read in observed burned area
274	LOGICAL, INTENT (in) :: read_observed_ba
275	! Observed burned area
276	REAL(r_std),DIMENSION (npts), INTENT (in) :: observed_ba
277
278	! Flag for read in observed burned area
279	LOGICAL, INTENT (in) :: read_cf_coarse
280	! Observed burned area
281	REAL(r_std),DIMENSION (npts), INTENT (in) :: cf_coarse
282	! Flag for read in observed burned area
283	LOGICAL, INTENT (in) :: read_cf_fine
284	! Observed burned area
285	REAL(r_std),DIMENSION (npts), INTENT (in) :: cf_fine
286	! Flag for read in observed burned area
287	LOGICAL, INTENT (in) :: read_ratio
288	! Observed burned area
289	REAL(r_std),DIMENSION (npts), INTENT (in) :: ratio
290	! Flag for read in observed burned area
291	LOGICAL, INTENT (in) :: read_ratio_flag
292	! Observed burned area
293	REAL(r_std),DIMENSION (npts), INTENT (in) :: ratio_flag
294
295	!endspit
296	!!!qcj++ peatland
297	REAL(r_std),DIMENSION(npts,nvm), INTENT(in) :: carbon_acro
298	REAL(r_std),DIMENSION(npts,nvm), INTENT(in) :: carbon_cato
299	REAL(r_std),DIMENSION(npts), INTENT(in) :: height_acro
300	REAL(r_std),DIMENSION(npts), INTENT(in) :: height_cato
301	REAL(r_std),DIMENSION(npts), INTENT(in) :: tcarbon_acro
302	REAL(r_std),DIMENSION(npts), INTENT(in) :: tcarbon_cato
303	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) ::resp_acro_oxic_d !!respiration of acrotelm( oxic ),daily
304	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) ::resp_acro_anoxic_d !!respiration of acrotelm( anoxic ),daily
305	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) ::resp_cato_d !!respiration of catotelm,daily
306	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) ::litter_to_acro_d
307	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) ::acro_to_cato_d
308	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) ::veget_max_adjusted
309	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: carbon_save
310	REAL(r_std), DIMENSION(npts,ndeep), INTENT(inout) :: deepC_a_save
311	REAL(r_std), DIMENSION(npts,ndeep), INTENT(inout) :: deepC_s_save
312	REAL(r_std), DIMENSION(npts,ndeep), INTENT(inout) :: deepC_p_save
313	REAL(r_std), DIMENSION(npts), INTENT(inout) :: delta_fsave
314	REAL(r_std),DIMENSION(npts), INTENT(in) :: liqwt_max_lastyear
315	REAL(r_std),DIMENSION(npts,nvm,nparts,nelements) :: biomass_remove
316
317
318	REAL(r_std), DIMENSION(npts), INTENT(in) :: litterhum_daily !! Daily litter humidity (0 to 1, unitless)
319	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: soilhum_daily !! Daily soil humidity (0 to 1, unitless)
320	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: maxmoiavail_lastyear !! Last year's maximum moisture availability
321	!! (0 to 1, unitless)
322	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: minmoiavail_lastyear !! Last year's minimum moisture availability
323	!! (0 to 1, unitless)
324	REAL(r_std), DIMENSION(npts), INTENT(in) :: gdd0_lastyear !! Last year's GDD0 (K)
325	REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_lastyear !! Lastyear's precipitation
326	!! @tex $(mm year^{-1})$ @endtex
327	!! to determine if establishment possible
328	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
329	!! unitless)
330	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
331	!! (0 to 1, unitless)
332	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_longterm !! "Long term" 2 meter reference
333	!! temperatures (K)
334	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_month !! "Monthly" 2-meter temperatures (K)
335	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_week !! "Weekly" 2-meter temperatures (K)
336	! "seasonal" 2-meter temperatures (K)
337	REAL(r_std), DIMENSION(npts), INTENT(in) :: Tseason
338	! temporary variable to calculate Tseason
339	REAL(r_std), DIMENSION(npts), INTENT(in) :: Tseason_length
340	! temporary variable to calculate Tseason
341	REAL(r_std), DIMENSION(npts), INTENT(in) :: Tseason_tmp
342
343	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: Tmin_spring_time !! Number of days after begin_leaves (leaf onset)
344	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: onset_date !! Date in the year at when the leaves started to grow(begin_leaves)
345
346	!pss:+
347	REAL(r_std), DIMENSION(npts), INTENT(in) :: tsurf_year ! annual surface temperatures (K)
348	!pss:-
349	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: tsoil_month !! "Monthly" soil temperatures (K)
350	REAL(r_std), DIMENSION(npts,nslm), INTENT(in) :: soilhum_month !! "Monthly" soil humidity
351	!! (0 to 1, unitless)
352	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
353	!! C (for phenology)
354	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit for crops
355	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
356	!! (for phenology) - this is written to the history files
357	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
358	!! leaves were lost (for phenology)
359	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
360	!! -5 deg C (for phenology)
361	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: turnover_longterm !! "Long term" turnover rate
362	!! @tex $(gC m^{-2} year^{-1})$ @endtex
363	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
364	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
365	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
366	!! @tex $(gC m^{-2} day^{-1})$ @endtex
367	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
368	!! availability (days)
369	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: maxfpc_lastyear !! Last year's maximum foliage projected
370	!! coverage for each natural PFT,
371	!! @tex $(m^2 m^{-2})$ @endtex
372	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: resp_maint_part !! Maintenance respiration of different
373	!! plant parts
374	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
375	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: fpc_max !! "Maximal" coverage fraction of a PFT (LAI
376	!! -> infinity) on ground
377	!! @tex $(m^2 m^{-2})$ @endtex
378	REAL(r_std), DIMENSION(ndeep), INTENT (in) :: zz_coef_deep !! deep vertical profile
379	!!!!! crop variables
380
381	! FOR CROP---STICS
382	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: pdlai
383	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: slai
384
385	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: in_cycle
386	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: deltai
387	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: dltaisen
388	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ssla
389	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: pgrain
390	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: deltgrain
391	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: reprac
392	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nger
393	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nlev
394	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: ndrp
395	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nmat
396	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nlax
397
398	! REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: N_limfert !!
399	! defined already in GM
400	INTEGER(i_std), INTENT(in) :: tday_counter
401
402	!!!!! end crop variables, xuhui
403
404	!gmjc
405	LOGICAL, INTENT(in) :: EndOfYear
406	!end gmjc
407	!! 0.2 Output variables
408
409	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: npp_daily !! Net primary productivity
410	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
411	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out) :: turnover_daily !! Turnover rates
412	!! @tex $(gC m^{-2} year^{-1})$ @endtex
413	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: co2_to_bm !! CO2 taken up from atmosphere when
414	!! introducing a new PFT (introduced for
415	!! carbon balance closure)
416	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
417	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: co2_fire !! Carbon emitted into the atmosphere by
418	!! fire (living and dead biomass)
419	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
420	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
421	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
422	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_maint !! Maintenance respiration
423	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
424	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_growth !! Growth respiration
425	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
426
427	REAL(r_std), DIMENSION(npts), INTENT(inout) :: deadleaf_cover !! Fraction of soil covered by dead leaves
428	!! (0 to 1, unitless)
429	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: vcmax !! Maximum rate of carboxylation
430	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out):: bm_to_litter !! Conversion of biomass to litter
431	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
432	LOGICAL, DIMENSION(npts,nvm), INTENT(out) :: begin_leaves !! signal to start putting leaves on (true/false)
433
434	! Wetland CH4 methane density
435	!pss:+
436	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_tot_0
437	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_dif_0
438	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_bub_0
439	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_pla_0
440
441	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_tot_wet1
442	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_dif_wet1
443	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_bub_wet1
444	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_pla_wet1
445
446	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_tot_wet2
447	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_dif_wet2
448	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_bub_wet2
449	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_pla_wet2
450
451	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_tot_wet3
452	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_dif_wet3
453	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_bub_wet3
454	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_pla_wet3
455
456	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_tot_wet4
457	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_dif_wet4
458	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_bub_wet4
459	REAL(r_std), DIMENSION(npts), INTENT(in) :: ch4_flux_density_pla_wet4
460	!pss:-
461
462	!!!!! crop variables
463	! for crop c pools
464	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: c_reserve
465	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: c_leafb
466
467	! for crop turnover
468	INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in) :: nrec !harvest date
469	!!!!! end crop variables, xuhui
470
471	!! 0.3 Modified variables
472
473	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: height !! Height of vegetation (m)
474	REAL(r_std), DIMENSION(npts,nlevs), INTENT(inout) :: control_moist !! Moisture control of heterotrophic
475	!! respiration (0 to 1, unitless)
476	REAL(r_std), DIMENSION(npts,nlevs), INTENT(inout) :: control_temp !! Temperature control of heterotrophic
477	!! respiration, above and below
478	!! (0 to 1, unitless)
479	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: soilcarbon_input !! Quantity of carbon going into carbon
480	!! pools from litter decomposition
481	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
482	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lai !! Leaf area index OF AN INDIVIDUAL PLANT,
483	!! where a PFT contains n indentical plants
484	!! i.e., using the mean individual approach
485	!! @tex $(m^2 m^{-2})$ @endtex
486	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: rprof !! Prescribed root depth (m)
487	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
488	!! each pixel
489	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age !! Age (years)
490	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: fireindex !! Probability of fire (0 to 1, unitless)
491	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: firelitter !! Longer term litter above the ground that
492	!! can be burned, @tex $(gC m^{-2})$ @endtex
493	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_age !! Leaf age (days)
494	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac !! Fraction of leaves in leaf age class,
495	!! (0 to 1, unitless)
496	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass !! Biomass @tex $(gC m^{-2})$ @endtex
497	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind !! Density of individuals
498	!! @tex $(m^{-2})$ @endtex
499	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: adapted !! Adaptation of PFT (killed if too cold)
500	!! (0 to 1, unitless)
501	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: regenerate !! "Fitness": Winter sufficiently cold for
502	!! PFT regeneration ? (0 to 1, unitless)
503	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
504	!! for deciduous trees)
505	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
506	!! the growing season (days)
507	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: litterpart !! Fraction of litter above the ground
508	!! belonging to different PFTs
509	!! (0 to 1, unitless)
510	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout):: litter !! Metabolic and structural litter, above
511	!! and below ground
512	!! @tex $(gC m^{-2})$ @endtex
513	!gmjc for grazing litter
514	REAL(r_std), DIMENSION(npts,nlitt,nvm), INTENT(out):: litter_avail
515	REAL(r_std), DIMENSION(npts,nlitt,nvm) , INTENT(out):: litter_not_avail
516	REAL(r_std), DIMENSION(npts,nlitt,nvm), INTENT(in):: litter_avail_frac
517	!end gmjc
518	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: dead_leaves !! Dead leaves on ground, per PFT, metabolic
519	!! and structural,
520	!! @tex $(gC m^{-2})$ @endtex
521	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: carbon !! Carbon pool: active, slow, or passive,
522	!! @tex $(gC m^{-2})$ @endtex
523	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: carbon_surf
524	!! @tex $(gC m^{-2})$ @endtex
525	REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout) :: lignin_struc !! Ratio of Lignin/Carbon in structural
526	!! litter, above and below ground,
527	!! @tex $(gC m^{-2})$ @endtex
528	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max !! "Maximal" coverage fraction of a PFT (LAI
529	!! -> infinity) on ground
530	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max_new !! "Maximal" coverage fraction of a PFT
531	!! (LAI-> infinity) on ground (unitless)
532	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary
533	!! productivity
534	!! @tex $(gC m^{-2} year^{-1})$ @endtex
535	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lm_lastyearmax !! Last year's maximum leaf mass, for each
536	!! PFT @tex $(gC m^{-2})$ @endtex
537	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_lastlight !! Vegetation fractions (on ground) after
538	!! last light competition
539	!! @tex $(m^2 m^{-2})$ @endtex
540	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: everywhere !! Is the PFT everywhere in the grid box or
541	!! very localized (after its introduction)
542	!! (unitless)
543	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: need_adjacent !! In order for this PFT to be introduced,
544	!! does it have to be present in an
545	!! adjacent grid box?
546	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: RIP_time !! How much time ago was the PFT eliminated
547	!! for the last time (y)
548	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: turnover_time !! Turnover_time of leaves for grasses
549	!! (days)
550	REAL(r_std), DIMENSION(npts,nvm),INTENT(in) :: vegetnew_firstday !! New "maximal" coverage fraction of a PFT
551	!! (LAI -> infinity) (unitless)
552	REAL(r_std),DIMENSION(npts,12),INTENT(inout) :: glccNetLCC
553	REAL(r_std),DIMENSION(npts,12),INTENT(inout) :: glccSecondShift
554	REAL(r_std),DIMENSION(npts,12),INTENT(inout) :: glccPrimaryShift
555	REAL(r_std), DIMENSION(npts,12),INTENT(inout) :: harvest_matrix !! The gross land use change matrix in case
556	!! of gross land cover change is simulated.
557	REAL(r_std), DIMENSION(npts,nvm),INTENT(inout) :: bound_spa !! The gross land use change matrix in case
558	!! of gross land cover change is simulated.
559	REAL(r_std),DIMENSION(npts,0:10,nwp), INTENT(inout) :: prod10 !! Products remaining in the 10
560	!! year-turnover pool after the annual
561	!! release for each compartment (10
562	!! + 1 : input from year of land cover
563	!! change) @tex $(gC m^{-2})$ @endtex
564	REAL(r_std),DIMENSION(npts,0:100,nwp), INTENT(inout) :: prod100 !! Products remaining in the 100
565	!! year-turnover pool after the annual
566	!! release for each compartment (100
567	!! + 1 : input from year of land cover
568	!! change) @tex $(gC m^{-2})$ @endtex
569	REAL(r_std),DIMENSION(npts,10,nwp), INTENT(inout) :: flux10 !! Annual release from the 10
570	!! year-turnover pool compartments
571	!! @tex $(gC m^{-2})$ @endtex
572	REAL(r_std),DIMENSION(npts,100,nwp), INTENT(inout) :: flux100 !! Annual release from the 100
573	!! year-turnover pool compartments
574	!! @tex $(gC m^{-2})$ @endtex
575	REAL(r_std),DIMENSION(npts,nwp), INTENT(inout) :: convflux !! Release during first year following land
576	!! cover change @tex $(gC m^{-2})$ @endtex
577	REAL(r_std),DIMENSION(npts,nwp), INTENT(inout) :: cflux_prod10 !! Total annual release from the 10
578	!! year-turnover pool
579	!! @tex $(gC m^{-2})$ @endtex
580	REAL(r_std),DIMENSION(npts,nwp), INTENT(inout) :: cflux_prod100 !! Total annual release from the 100
581	!! year-turnover pool
582	!! @tex $(gC m^{-2})$ @endtex
583	REAL(r_std), DIMENSION(npts), INTENT(inout) :: harvest_above !! Harvest above ground biomass for
584	!! agriculture @tex $(gC m^{-2})$ @endtex
585	REAL(r_std), DIMENSION(npts), INTENT(inout) :: carb_mass_total !! Carbon Mass total (soil, litter, veg)
586	!! @tex $(gC m^{-2})$ @endtex
587	REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(inout) :: MatrixA !! Matrix containing the fluxes
588	!! between the carbon pools
589	!! per sechiba time step
590	!! @tex $(gC.m^2.day^{-1})$ @endtex
591	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_a !! permafrost soil carbon (g/m**3) active
592	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_s !! permafrost soil carbon (g/m**3) slow
593	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_p !! permafrost soil carbon (g/m**3) passive
594	!gmjc
595
596	!glcc
597
598	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: glcc_pft !! a temporary variable to hold the fraction of ipft->ivma, i.e., from
599	!! PFT_{ipft} to the youngest age class of MTC_{ivma}
600	!spitfire
601	! Nesterov index accumulated
602	REAL(r_std), DIMENSION(npts), INTENT(inout) :: ni_acc
603	REAL(r_std), DIMENSION(npts), INTENT(inout) :: fire_numday
604	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: bafrac_deforest_accu
605	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: emideforest_litter_accu
606	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: emideforest_biomass_accu
607	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) :: deforest_litter_remain !! Vegetmax-weighted remaining litter on the ground for
608	!! deforestation region.
609	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: deforest_biomass_remain !! Vegetmax-weighted remaining biomass on the ground for
610	!! deforestation region.
611	! fuel classes (1, 10, 100, 1000 hours)
612	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: fuel_1hr
613	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: fuel_10hr
614	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: fuel_100hr
615	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: fuel_1000hr
616	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: def_fuel_1hr_remain
617	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: def_fuel_10hr_remain
618	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: def_fuel_100hr_remain
619	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: def_fuel_1000hr_remain
620	REAL(r_std), DIMENSION(npts) :: d_area_burnt
621	REAL(r_std), DIMENSION(npts) :: d_numfire
622	REAL(r_std), DIMENSION(npts) :: fc_crown
623	! parameter for potential human-caused ignitions, ignitions ind^{-1}day{-1}, used in lpj_spitfire.f90
624	REAL(r_std), DIMENSION(npts), INTENT(in) :: a_nd
625	!endspit
626	!gmjc
627	! snowfall_daily (mm/d?)
628	REAL(r_std), DIMENSION(npts), INTENT(in) :: snowfall_daily
629	! snowmass_daily (kg/m2)
630	REAL(r_std), DIMENSION(npts), INTENT(in) :: snowmass_daily
631	! "14days" 2-meter temperatures (K)
632	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_14
633	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: sla_calc
634	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: wshtotsum
635	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: sr_ugb
636	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: compt_ugb
637	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: nb_ani
638	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: grazed_frac
639	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: import_yield
640	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: sla_age1
641	INTEGER(i_std), INTENT(in) :: day_of_year
642	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: N_limfert
643	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: when_growthinit_cut
644	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: nb_grazingdays
645	! REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: resp_hetero_litter_d
646	! REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: resp_hetero_soil_d
647	! top 5 layer grassland soil moisture for grazing
648	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: moiavail_daily
649	!! Daily moisture availability (0-1, unitless)
650	REAL(r_std),DIMENSION (npts), INTENT(in) :: tmc_topgrass_daily
651	REAL(r_std),DIMENSION (npts), INTENT(in) :: fc_grazing
652	REAL(r_std),DIMENSION (npts), INTENT(inout) :: after_snow
653	REAL(r_std),DIMENSION (npts), INTENT(inout) :: after_wet
654	REAL(r_std),DIMENSION (npts), INTENT(inout) :: wet1day
655	REAL(r_std),DIMENSION (npts), INTENT(inout) :: wet2day
656	!end gmjc
657	!! 0.4 Local variables
658
659	REAL(r_std), DIMENSION(npts,nvm,nelements) :: tot_bm_to_litter !! Total conversion of biomass to litter
660	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
661	REAL(r_std), DIMENSION(npts,nvm,nelements) :: tot_live_biomass !! Total living biomass
662	!! @tex $(gC m{-2})$ @endtex
663	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements) :: bm_alloc !! Biomass increase, i.e. NPP per plant part
664	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
665	REAL(r_std), DIMENSION(npts,nvm,nelements) :: tot_turnover !! Total turnover rate
666	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
667	REAL(r_std), DIMENSION(npts,nvm) :: tot_litter_soil_carb!! Total soil and litter carbon
668	!! @tex $(gC m^{-2})$ @endtex
669	REAL(r_std), DIMENSION(npts,nvm) :: tot_litter_carb !! Total litter carbon
670	!! @tex $(gC m^{-2})$ @endtex
671	REAL(r_std), DIMENSION(npts,nvm) :: tot_soil_carb !! Total soil carbon
672	!! @tex $(gC m^{-2})$ @endtex
673	REAL(r_std), DIMENSION(npts) :: carb_mass_variation !! Carbon Mass variation
674	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
675	REAL(r_std), DIMENSION(npts,nvm) :: cn_ind !! Crown area of individuals
676	!! @tex $(m^{2})$ @endtex
677	REAL(r_std), DIMENSION(npts,nvm) :: woodmass_ind !! Woodmass of individuals (gC)
678	REAL(r_std), DIMENSION(npts,nvm,nparts) :: f_alloc !! Fraction that goes into plant part
679	!! (0 to 1, unitless)
680	REAL(r_std), DIMENSION(npts) :: avail_tree !! Space availability for trees
681	!! (0 to 1, unitless)
682	REAL(r_std), DIMENSION(npts) :: avail_grass !! Space availability for grasses
683	!! (0 to 1, unitless)
684	INTEGER :: j,ivm,ivma,ipts,ilev,ilitt
685	REAL(r_std),DIMENSION(npts) :: prod10_total !! Total products remaining in the pool
686	!! after the annual release
687	!! @tex $(gC m^{-2})$ @endtex
688	REAL(r_std),DIMENSION(npts) :: prod100_total !! Total products remaining in the pool
689	!! after the annual release
690	!! @tex $(gC m^{-2})$ @endtex
691	REAL(r_std),DIMENSION(npts) :: cflux_prod_total !! Total flux from conflux and the 10/100
692	!! year-turnover pool
693	!! @tex $(gC m^{-2} year^{-1})$ @endtex
694	REAL(r_std),DIMENSION(npts,nvm) :: veget_max_tmp !! "Maximal" coverage fraction of a PFT
695	!! (LAI-> infinity) on ground (unitless)
696	!!!!! crop
697	REAL(r_std), DIMENSION(npts,nvm) :: crop_export !! Cropland export (harvest & straw)
698	!!!!! end crop, xuhui
699	REAL(r_std), DIMENSION(npts,nvm) :: mortality !! Fraction of individual dying this time
700	!! step (0 to 1, unitless)
701	REAL(r_std), DIMENSION(npts) :: vartmp !! Temporary variable used to add history
702	REAL(r_std), DIMENSION(npts,nvm) :: histvar !! History variables
703	REAL(r_std), DIMENSION(npts,nvm) :: lcc !! land cover change, i.e., loss of each PFT,
704	!! positive value indicates loss.
705	REAL(r_std),DIMENSION(npts,nvm) :: deflitsup_total
706	REAL(r_std),DIMENSION(npts,nvm) :: defbiosup_total
707	!glcc
708	REAL(r_std), DIMENSION(npts,nvm,nvmap) :: glcc_pftmtc !! a temporary variable to hold the fractions each PFT is going to lose
709	REAL(r_std), DIMENSION(npts,12) :: glccReal !! The "real" glcc matrix that we apply in the model
710	!! after considering the consistency between presribed
711	!! glcc matrix and existing vegetation fractions.
712	REAL(r_std), DIMENSION(npts,12) :: IncreDeficit !! "Increment" deficits, negative values mean that
713	!! there are not enough fractions in the source PFTs
714	!! /vegetations to target PFTs/vegetations. I.e., these
715	!! fraction transfers are presribed in LCC matrix but
716	!! not realized.
717	REAL(r_std), DIMENSION(npts) :: Deficit_pf2yf_final !!
718	REAL(r_std), DIMENSION(npts) :: Deficit_sf2yf_final !!
719	REAL(r_std), DIMENSION(npts) :: pf2yf_compen_sf2yf !!
720	REAL(r_std), DIMENSION(npts) :: sf2yf_compen_pf2yf !!
721	INTEGER :: f2g=1, f2p=2, f2c=3
722	INTEGER :: g2f=4, g2p=5, g2c=6, p2f=7, p2g=8, p2c=9, c2f=10, c2g=11, c2p=12
723	REAL(r_std) :: valtmp
724
725	!gmjc lcchange of managed grassland
726	! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground
727	INTEGER(i_std) :: ier
728	LOGICAL :: l_error =.FALSE.
729	! variables to get closure carbon cycle for nbp
730	REAL(r_std), DIMENSION(npts,nvm) :: harvest_gm
731	REAL(r_std), DIMENSION(npts,nvm) :: ranimal_gm
732	REAL(r_std), DIMENSION(npts,nvm) :: ch4_pft_gm
733	REAL(r_std), DIMENSION(npts) :: ch4_gm
734	REAL(r_std), DIMENSION(npts,nvm) :: cinput_gm
735	REAL(r_std), DIMENSION(npts) :: co2_gm
736	REAL(r_std),DIMENSION(npts,nvm) :: veget_max_gm
737	REAL(r_std), DIMENSION(npts) :: veget_exist_gm
738	REAL(r_std),DIMENSION(npts,nvm) :: n2o_pft_gm
739	REAL(r_std), DIMENSION(npts) :: n2o_gm
740	!end gmjc
741
742
743
744	REAL(r_std), DIMENSION(npts,7) :: outflux_sta,outflux_end
745	REAL(r_std), DIMENSION(npts,3) :: influx_sta,influx_end
746	REAL(r_std), DIMENSION(npts,4) :: pool_sta,pool_end
747	INTEGER(i_std) :: ind_biomass,ind_litter,ind_soil,ind_prod,ind_co2tobm,ind_gpp,ind_npp,&
748	ind_bm2lit,ind_resph,ind_respm,ind_respg,ind_convf,ind_cflux,ind_fire
749
750
751	REAL(r_std), DIMENSION(npts) :: biomass_check
752	INTEGER ::i,k,m,l
753	!_ ================================================================================================================================
754
755	IF (printlev>=3) WRITE(numout,*) 'Entering stomate_lpj'
756
757	!! 1. Initializations
758
759	lcc(:,:) = zero
760	glccReal(:,:) = zero
761	glcc_pftmtc(:,:,:) = zero
762
763	!gmjc
764	IF (firstcall) THEN
765
766	firstcall = .FALSE.
767
768	!Config Key = GRM_ENABLE_GRAZING
769	!Config Desc = grazing allowed
770	!Config Def = n
771	!Config Help = flag for choose if you want animals or not.
772	!
773	enable_grazing = .FALSE.
774	CALL getin_p('GRM_ENABLE_GRAZING',enable_grazing)
775	WRITE (numout,*) 'enable_grazing',enable_grazing
776	WRITE (numout,*) 'manage',is_grassland_manag
777	WRITE (numout,*) 'cut',is_grassland_cut
778	WRITE (numout,*) 'grazed',is_grassland_grazed
779
780	emideforest_biomass_accu(:,:,:,:) = zero
781	emideforest_litter_accu(:,:,:,:) = zero
782	bafrac_deforest_accu(:,:) = zero
783
784	IF (do_now_stomate_lcchange) THEN
785	IF (use_age_class) THEN
786	IF (SingleAgeClass) THEN
787	CALL gross_glcc_firstday_SinAgeC_fh(npts,veget_max,harvest_matrix, &
788	glccSecondShift,glccPrimaryShift,glccNetLCC,&
789	glccReal,glcc_pft,glcc_pftmtc,IncreDeficit, &
790	Deficit_pf2yf_final, Deficit_sf2yf_final, &
791	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
792
793	ELSE
794	CALL gross_glcc_firstday_fh(npts,veget_max,harvest_matrix, &
795	glccSecondShift,glccPrimaryShift,glccNetLCC,&
796	glccReal,glcc_pft,glcc_pftmtc,IncreDeficit, &
797	Deficit_pf2yf_final, Deficit_sf2yf_final, &
798	pf2yf_compen_sf2yf, sf2yf_compen_pf2yf)
799	ENDIF
800	! We put only the conversion of tree->Notree as deforestation
801	DO ipts = 1,npts
802	DO ivm = 1,nvm
803	DO ivma = 1,nvmap
804	IF (is_tree(ivm) .AND. .NOT. is_tree(start_index(ivma))) THEN
805	lcc(ipts,ivm) = lcc(ipts,ivm) + glcc_pftmtc(ipts,ivm,ivma)
806	ENDIF
807	ENDDO
808	ENDDO
809	ENDDO
810	ELSE ! (.NOT. use_age_class), i.e., net land cover change.
811	! note here veget_max is the last-year veget_max; vegetnew_firstday
812	! is the veget_max of next year, we need lcc to be >0 where forest
813	! area decreased, i.e., when vegetnew_firstday < veget_max
814	lcc(:,:) = veget_max(:,:) - vegetnew_firstday(:,:)
815	ENDIF
816
817	IF (.NOT. allow_deforest_fire) lcc(:,:) = zero
818
819	! we creat the proxy that's needed for deforestation fire simulation.
820	DO ipts = 1,npts
821	DO ivm = 1,nvm
822	deforest_litter_remain(ipts,:,ivm,:,:) = litter(ipts,:,ivm,:,:)*lcc(ipts,ivm)
823	deforest_biomass_remain(ipts,ivm,:,:) = biomass(ipts,ivm,:,:)*lcc(ipts,ivm)
824	def_fuel_1hr_remain(ipts,ivm,:,:) = fuel_1hr(ipts,ivm,:,:)*lcc(ipts,ivm)
825	def_fuel_10hr_remain(ipts,ivm,:,:) = fuel_10hr(ipts,ivm,:,:)*lcc(ipts,ivm)
826	def_fuel_100hr_remain(ipts,ivm,:,:) = fuel_100hr(ipts,ivm,:,:)*lcc(ipts,ivm)
827	def_fuel_1000hr_remain(ipts,ivm,:,:) = fuel_1000hr(ipts,ivm,:,:)*lcc(ipts,ivm)
828	ENDDO
829	ENDDO
830	ENDIF
831
832
833	END IF !firstcall
834
835
836	!! 1.1 Initialize variables to zero
837	co2_to_bm(:,:) = zero
838	co2_fire(:,:) = zero
839	npp_daily(:,:) = zero
840	resp_maint(:,:) = zero
841	resp_growth(:,:) = zero
842	harvest_above(:) = zero
843	bm_to_litter(:,:,:,:) = zero
844	cn_ind(:,:) = zero
845	woodmass_ind(:,:) = zero
846	turnover_daily(:,:,:,:) = zero
847	crop_export(:,:) = zero
848	deflitsup_total(:,:) = zero
849	defbiosup_total(:,:) = zero
850	!gmjc
851	!! Initialize gm variables for nbp to zero
852	harvest_gm(:,:) = zero
853	ranimal_gm(:,:) = zero
854	ch4_pft_gm(:,:) = zero
855	cinput_gm(:,:) = zero
856	co2_gm(:) = zero
857	ch4_gm(:) = zero
858	n2o_gm(:) = zero
859	n2o_pft_gm(:,:) = zero
860	veget_max_gm(:,:) = zero
861	veget_exist_gm(:) = zero
862	!end gmjc
863
864	!! 1.2 Initialize variables to veget_max
865	veget_max_tmp(:,:) = veget_max(:,:)
866
867	!! 1.3 Calculate some vegetation characteristics
868	!! 1.3.1 Calculate some vegetation characteristics
869	! Calculate cn_ind (individual crown mass) and individual height from
870	! state variables if running DGVM or dynamic mortality in static cover mode
871	!?? Explain (maybe in the header once) why you mulitply with veget_max in the DGVM
872	!?? and why you don't multiply with veget_max in stomate.
873	IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
874	IF(ok_dgvm) THEN
875	WHERE (ind(:,:).GT.min_stomate)
876	woodmass_ind(:,:) = &
877	((biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
878	+biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon)) &
879	*veget_max(:,:))/ind(:,:)
880	ENDWHERE
881	ELSE
882	WHERE (ind(:,:).GT.min_stomate)
883	woodmass_ind(:,:) = &
884	(biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
885	+biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon))/ind(:,:)
886	ENDWHERE
887	ENDIF
888
889	CALL crown (npts, PFTpresent, &
890	ind, biomass, woodmass_ind, &
891	veget_max, cn_ind, height)
892	ENDIF
893
894	!! 1.3.2 Prescribe characteristics if the vegetation is not dynamic
895	! IF the DGVM is not activated, the density of individuals and their crown
896	! areas don't matter, but they should be defined for the case we switch on
897	! the DGVM afterwards. At the first call, if the DGVM is not activated,
898	! impose a minimum biomass for prescribed PFTs and declare them present.
899	CALL prescribe (npts, &
900	veget_max, dt_days, PFTpresent, everywhere, when_growthinit, &
901	biomass, leaf_frac, ind, cn_ind, co2_to_bm)
902
903
904	!! 2. Climatic constraints for PFT presence and regenerativeness
905
906	! Call this even when DGVM is not activated so that "adapted" and "regenerate"
907	! are kept up to date for the moment when the DGVM is activated.
908	CALL constraints (npts, dt_days, &
909	t2m_month, t2m_min_daily,when_growthinit, &
910	adapted, regenerate, Tseason)
911
912
913	!! 3. Determine introduction and elimination of PTS based on climate criteria
914
915	IF ( ok_dgvm ) THEN
916
917	!! 3.1 Calculate introduction and elimination
918	CALL pftinout (npts, dt_days, adapted, regenerate, &
919	neighbours, veget_max, &
920	biomass, ind, cn_ind, age, leaf_frac, npp_longterm, lm_lastyearmax, senescence, &
921	PFTpresent, everywhere, when_growthinit, need_adjacent, RIP_time, &
922	co2_to_bm, &
923	avail_tree, avail_grass, &
924	!gmjc
925	sla_calc)
926	!end gmjc
927	!! 3.2 Reset attributes for eliminated PFTs.
928	! This also kills PFTs that had 0 leafmass during the last year. The message
929	! "... after pftinout" is misleading in this case.
930	CALL kill (npts, 'pftinout ', lm_lastyearmax, &
931	ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
932	lai, age, leaf_age, leaf_frac, npp_longterm, &
933	when_growthinit, everywhere, veget_max, bm_to_litter)
934
935
936	!! 3.3 Calculate new crown area and diameter
937	! Calculate new crown area, diameter and maximum vegetation cover**[No longer used in the subroutine]
938	! unsure whether this is really required
939	! - in theory this could ONLY be done at the END of stomate_lpj
940	! calculate woodmass of individual tree
941	IF(ok_dgvm) THEN
942	WHERE (ind(:,:).GT.min_stomate)
943	woodmass_ind(:,:) = &
944	((biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
945	+biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon)) &
946	*veget_max(:,:))/ind(:,:)
947	ENDWHERE
948	ELSE
949	WHERE (ind(:,:).GT.min_stomate)
950	woodmass_ind(:,:) = &
951	(biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
952	+biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon))/ind(:,:)
953	ENDWHERE
954	ENDIF
955
956	! Calculate crown area and diameter for all PFTs (including the newly established)
957	CALL crown (npts, PFTpresent, &
958	ind, biomass, woodmass_ind, &
959	veget_max, cn_ind, height)
960
961	ENDIF
962
963	!! 4. Phenology
964
965	!! 4.1 Write values to history file
966	! Current values for ::when_growthinit
967	CALL xios_orchidee_send_field("WHEN_GROWTHINIT",when_growthinit)
968
969	CALL histwrite_p (hist_id_stomate, 'WHEN_GROWTHINIT', itime, when_growthinit, npts*nvm, horipft_index)
970
971	! Set and write values for ::PFTpresent
972	WHERE(PFTpresent)
973	histvar=un
974	ELSEWHERE
975	histvar=zero
976	ENDWHERE
977
978	CALL xios_orchidee_send_field("PFTPRESENT",histvar)
979
980	CALL histwrite_p (hist_id_stomate, 'PFTPRESENT', itime, histvar, npts*nvm, horipft_index)
981
982	! Set and write values for gdd_midwinter
983	WHERE(gdd_midwinter.EQ.undef)
984	histvar=val_exp
985	ELSEWHERE
986	histvar=gdd_midwinter
987	ENDWHERE
988
989	CALL xios_orchidee_send_field("GDD_MIDWINTER",histvar)
990
991	CALL histwrite_p (hist_id_stomate, 'GDD_MIDWINTER', itime, histvar, npts*nvm, horipft_index)
992
993	! Set and write values for gdd_m5_dormance
994	WHERE(gdd_m5_dormance.EQ.undef)
995	histvar=val_exp
996	ELSEWHERE
997	histvar=gdd_m5_dormance
998	ENDWHERE
999
1000	CALL xios_orchidee_send_field('GDD_M5_DORMANCE',histvar)
1001	CALL histwrite_p (hist_id_stomate, 'GDD_M5_DORMANCE', itime, histvar, npts*nvm, horipft_index)
1002
1003	! Set and write values for ncd_dormance
1004	WHERE(ncd_dormance.EQ.undef)
1005	histvar=val_exp
1006	ELSEWHERE
1007	histvar=ncd_dormance
1008	ENDWHERE
1009
1010	CALL xios_orchidee_send_field("NCD_DORMANCE",histvar)
1011
1012	CALL histwrite_p (hist_id_stomate, 'NCD_DORMANCE', itime, histvar, npts*nvm, horipft_index)
1013
1014	!! 4.2 Calculate phenology
1015	! WRITE(numout,*) 'slai before lpj_phenology: ',slai(1,12:14)
1016	CALL phenology (npts, dt_days, PFTpresent, &
1017	veget_max, &
1018	t2m_longterm, t2m_month, t2m_week, gpp_daily, &
1019	maxmoiavail_lastyear, minmoiavail_lastyear, &
1020	moiavail_month, moiavail_week, &
1021	gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
1022	senescence, time_hum_min, &
1023	biomass, leaf_frac, leaf_age, &
1024	when_growthinit, co2_to_bm, &
1025	pdlai, slai, deltai, ssla, &
1026	begin_leaves, &!)
1027	!gmjc
1028	sla_calc)
1029	!end gmjc
1030	! IF (printlev>=4) THEN
1031	! WRITE(,) 'biomass reserve after phenology is', biomass(:,12:14,icarbres, icarbon)
1032	! ENDIF
1033
1034	! IF ( ANY(biomass(1,14,:,icarbon) .lt. 0. ) ) THEN
1035	! WRITE(numout,*) 'biomass low0 after phenology'
1036	! WRITE(numout,*) 'biomass(1,14,:,icarbon)',biomass(1,14,:,icarbon)
1037	! ENDIF
1038
1039	!! 5. Allocate C to different plant parts
1040	!WRITE(numout,*) 'slai before lpj_alloc: ',slai(1,12:14)
1041	CALL alloc (npts, dt_days, &
1042	lai, veget_max, senescence, when_growthinit, &
1043	moiavail_week, tsoil_month, soilhum_month, &
1044	biomass, age, leaf_age, leaf_frac, rprof, f_alloc, &!)
1045	deltai, ssla, & !added for crop by xuhui
1046	!gmjc
1047	sla_calc, when_growthinit_cut)
1048	!end gmjc
1049	! IF (printlev>=4) THEN
1050	! WRITE(,) 'biomass reserve after alloc is', biomass(:,12:14, icarbres,icarbon)
1051	! ENDIF
1052
1053	! IF ( ANY(biomass(1,14,:,icarbon) .lt. 0. ) ) THEN
1054	! WRITE(numout,*) 'biomass low0 after alloc'
1055	! WRITE(numout,*) 'biomass(1,14,:,icarbon)',biomass(1,14,:,icarbon)
1056	! ENDIF
1057	!! 5.1. Recalculate lai
1058	!! This should be done whenever biomass is modified
1059	CALL setlai(biomass, sla_calc, slai)
1060
1061	!! 6. NPP, maintenance and growth respiration
1062	! WRITE(numout,*) 'slai before lpj_npp_calc: ',slai(1,12:14)
1063	!! 6.1 Calculate NPP and respiration terms
1064	CALL npp_calc (npts, dt_days, &
1065	PFTpresent, &
1066	t2m_daily, tsoil_daily, lai, rprof, &
1067	gpp_daily, f_alloc, bm_alloc, resp_maint_part,&
1068	biomass, leaf_age, leaf_frac, age, &
1069	resp_maint, resp_growth, npp_daily, &!)
1070	! for crop bm_alloc
1071	!!! crop variables
1072	in_cycle, deltai, dltaisen, ssla, pgrain, deltgrain, reprac, &
1073	nger, nlev, ndrp, nlax, nmat, nrec, &
1074	c_reserve, c_leafb, slai, tday_counter, veget_max, &
1075	!!! end crop, xuhui
1076	!gmjc
1077	sla_calc, sla_age1,N_limfert)
1078	!end gmjc
1079	IF (printlev>=4) THEN
1080	WRITE(,) 'biomass reserve after npp_calc is', biomass(1,:,icarbres, icarbon)
1081	ENDIF
1082
1083	! IF ( ANY(biomass(1,14,:,icarbon) .lt. 0. ) ) THEN
1084	! WRITE(numout,*) 'biomass low0 after npp_calc'
1085	! WRITE(numout,*) 'biomass(1,14,:,icarbon)',biomass(1,14,:,icarbon)
1086	! ENDIF
1087
1088
1089	!! 6.2 Kill slow growing PFTs in DGVM or STOMATE with constant mortality
1090	IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
1091	CALL kill (npts, 'npp ', lm_lastyearmax, &
1092	ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
1093	lai, age, leaf_age, leaf_frac, npp_longterm, &
1094	when_growthinit, everywhere, veget_max, bm_to_litter)
1095
1096	!! 6.2.1 Update wood biomass
1097	! For the DGVM
1098	IF(ok_dgvm) THEN
1099	WHERE (ind(:,:).GT.min_stomate)
1100	woodmass_ind(:,:) = &
1101	((biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
1102	+ biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon)) &
1103	*veget_max(:,:))/ind(:,:)
1104	ENDWHERE
1105
1106	! For all pixels with individuals
1107	ELSE
1108	WHERE (ind(:,:).GT.min_stomate)
1109	woodmass_ind(:,:) = &
1110	(biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
1111	+ biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))/ind(:,:)
1112	ENDWHERE
1113	ENDIF ! ok_dgvm
1114
1115	!! 6.2.2 New crown area and maximum vegetation cover after growth
1116	CALL crown (npts, PFTpresent, &
1117	ind, biomass, woodmass_ind,&
1118	veget_max, cn_ind, height)
1119
1120	ENDIF ! ok_dgvm
1121
1122
1123	!! 7. fire
1124	!WRITE(numout,*) 'slai before lpj_fire: ',slai(1,12:14)
1125	!! 7.1. Burn PFTs
1126	!CALL fire (npts, dt_days, litterpart, &
1127	! litterhum_daily, t2m_daily, lignin_struc, veget_max, &
1128	! fireindex, firelitter, biomass, ind, &
1129	! litter, dead_leaves, bm_to_litter, &
1130	! co2_fire, MatrixA)
1131
1132	!gmjc update available and not available litter for grazing litter
1133	!spitfire
1134
1135	!disable_fire and allow_deforest_fire are defined as constants in src_parameters/constantes_var.f90
1136	!disable_fire to DISABLE fire when being TRUE
1137	!allow_deforest_fire to activate deforestation fire module when being TRUE
1138	IF(.NOT.disable_fire) THEN
1139	CALL spitfire(npts, dt_days, veget_max,resolution,contfrac, &
1140	PFTpresent,t2m_min_daily,t2m_max_daily, &
1141	precip_daily,wspeed_daily,soilhum_daily(:,1), &
1142	lightn,litter(:,:,:,:,icarbon),ni_acc,fire_numday, &
1143	fuel_1hr(:,:,:,icarbon),fuel_10hr(:,:,:,icarbon),fuel_100hr(:,:,:,icarbon), &
1144	fuel_1000hr(:,:,:,icarbon),ind,biomass(:,:,:,icarbon),popd,a_nd,height, &
1145	read_observed_ba, observed_ba, read_cf_fine,cf_fine, &
1146	read_cf_coarse,cf_coarse,read_ratio_flag, &
1147	ratio_flag,read_ratio,ratio,date, &
1148	bm_to_litter(:,:,:,icarbon),co2_fire, &
1149	lcc,bafrac_deforest_accu,emideforest_litter_accu(:,:,:,icarbon),emideforest_biomass_accu(:,:,:,icarbon),&
1150	deforest_litter_remain(:,:,:,:,icarbon),deforest_biomass_remain(:,:,:,icarbon),&
1151	def_fuel_1hr_remain(:,:,:,icarbon),def_fuel_10hr_remain(:,:,:,icarbon),&
1152	def_fuel_100hr_remain(:,:,:,icarbon),def_fuel_1000hr_remain(:,:,:,icarbon))
1153	ELSE
1154	ni_acc=0.
1155	fire_numday=0.
1156
1157	ENDIF
1158	!endspit
1159	! after fire burning
1160	litter_avail(:,:,:) = litter(:,:,:,iabove,icarbon) * &
1161	litter_avail_frac(:,:,:)
1162	litter_not_avail(:,:,:) = litter(:,:,:,iabove,icarbon) * &
1163	(1.0 - litter_avail_frac(:,:,:))
1164	!end gmjc
1165	!! 7.2 Kill PFTs in DGVM
1166	IF ( ok_dgvm ) THEN
1167
1168	! reset attributes for eliminated PFTs
1169	CALL kill (npts, 'fire ', lm_lastyearmax, &
1170	ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
1171	lai, age, leaf_age, leaf_frac, npp_longterm, &
1172	when_growthinit, everywhere, veget_max, bm_to_litter)
1173
1174	ENDIF ! ok_dgvm
1175
1176	!! 8. Tree mortality
1177
1178	! Does not depend on age, therefore does not change crown area.
1179	CALL gap (npts, dt_days, &
1180	npp_longterm, turnover_longterm, lm_lastyearmax, &
1181	PFTpresent, biomass, ind, bm_to_litter, mortality, t2m_min_daily, Tmin_spring_time, &!)
1182	!gmjc
1183	sla_calc)
1184	!end gmjc
1185
1186	IF ( ok_dgvm ) THEN
1187
1188	! reset attributes for eliminated PFTs
1189	CALL kill (npts, 'gap ', lm_lastyearmax, &
1190	ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
1191	lai, age, leaf_age, leaf_frac, npp_longterm, &
1192	when_growthinit, everywhere, veget_max, bm_to_litter)
1193
1194	ENDIF
1195
1196	!! 10. Leaf senescence, new lai and other turnover processes
1197	! IF ( ANY(biomass(1,14,:,icarbon) .lt. 0. ) ) THEN
1198	! WRITE(numout,*) 'biomass low0 before turn'
1199	! WRITE(numout,*) 'biomass(1,14,:,icarbon)',biomass(1,14,:,icarbon)
1200	! ENDIF
1201
1202	CALL turn (npts, dt_days, PFTpresent, &
1203	herbivores, &
1204	maxmoiavail_lastyear, minmoiavail_lastyear, &
1205	moiavail_week, moiavail_month,t2m_longterm, t2m_month, t2m_week, veget_max, &
1206	gdd_from_growthinit, leaf_age, leaf_frac, age, lai, biomass, &
1207	turnover_daily, senescence,turnover_time, &!)
1208	nrec,crop_export, &
1209	!gmjc
1210	sla_calc)
1211	!end gmjc
1212	! IF ( ANY(biomass(1,14,:,icarbon) .lt. 0. ) ) THEN
1213	! WRITE(numout,*) 'biomass low0 after turn'
1214	! WRITE(numout,*) 'biomass(1,14,:,icarbon)',biomass(1,14,:,icarbon)
1215	! ENDIF
1216	!! 11. Light competition
1217
1218	!! If not using constant mortality then kill with light competition
1219	! IF ( ok_dgvm .OR. .NOT.(lpj_gap_const_mort) ) THEN
1220	IF ( ok_dgvm ) THEN
1221
1222	!! 11.1 Light competition
1223	CALL light (npts, dt_days, &
1224	veget_max, fpc_max, PFTpresent, cn_ind, lai, maxfpc_lastyear, &
1225	lm_lastyearmax, ind, biomass, veget_lastlight, bm_to_litter, mortality, &!)
1226	!gmjc
1227	sla_calc)
1228	!end gmjc
1229
1230	!! 11.2 Reset attributes for eliminated PFTs
1231	CALL kill (npts, 'light ', lm_lastyearmax, &
1232	ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
1233	lai, age, leaf_age, leaf_frac, npp_longterm, &
1234	when_growthinit, everywhere, veget_max, bm_to_litter)
1235
1236	ENDIF
1237
1238
1239	!! 12. Establishment of saplings
1240
1241	IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort ) THEN
1242
1243	!! 12.1 Establish new plants
1244	CALL establish (npts, lalo, dt_days, PFTpresent, regenerate, &
1245	neighbours, resolution, need_adjacent, herbivores, &
1246	precip_lastyear, gdd0_lastyear, lm_lastyearmax, &
1247	cn_ind, lai, avail_tree, avail_grass, npp_longterm, &
1248	leaf_age, leaf_frac, &
1249	ind, biomass, age, everywhere, co2_to_bm, veget_max, woodmass_ind, &
1250	mortality, bm_to_litter, &
1251	!gmjc
1252	sla_calc)
1253	!end gmjc
1254
1255	!! 12.2 Calculate new crown area (and maximum vegetation cover)
1256	CALL crown (npts, PFTpresent, &
1257	ind, biomass, woodmass_ind, &
1258	veget_max, cn_ind, height)
1259
1260	ENDIF
1261	!gmjc Grassland_management
1262	!
1263	! 13 calculate grazing by animals or cutting for forage
1264	!
1265	IF (enable_grazing) THEN
1266	CALL main_grassland_management(&
1267	npts, lalo, neighbours, resolution, contfrac, &
1268	dt_days , &
1269	day_of_year , &
1270	t2m_daily , &
1271	t2m_min_daily , &
1272	t2m_14 , &
1273	tsurf_daily , &
1274	snowfall_daily , &
1275	biomass , &
1276	bm_to_litter , &
1277	litter , &
1278	litter_avail , &
1279	litter_not_avail , &
1280	!spitfire
1281	fuel_1hr(:,:,:,icarbon), &
1282	fuel_10hr(:,:,:,icarbon), &
1283	fuel_100hr(:,:,:,icarbon), &
1284	fuel_1000hr(:,:,:,icarbon), &
1285	!end spitfire
1286	.TRUE. , &
1287	EndOfYear , &
1288	when_growthinit_cut, nb_grazingdays, &
1289	lai,sla_calc,leaf_age,leaf_frac, &
1290	wshtotsum,sr_ugb,compt_ugb, &
1291	nb_ani,grazed_frac,import_yield,N_limfert, &
1292	moiavail_daily,tmc_topgrass_daily,fc_grazing, snowmass_daily, &
1293	after_snow, after_wet, wet1day, wet2day, &
1294	harvest_gm, ranimal_gm, ch4_pft_gm, cinput_gm, n2o_pft_gm)
1295	ENDIF
1296	!end gmjc
1297
1298	!! 13. Calculate final LAI and vegetation cover
1299
1300	!!!qcj++ peatland
1301	IF (update_peatfrac) THEN
1302	CALL lpj_cover_peat(npts,lalo, cn_ind, ind, biomass, veget_max_new, veget_max, &
1303	veget_max_tmp, litter, litter_avail, litter_not_avail, carbon, &
1304	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1305	turnover_daily, bm_to_litter, &
1306	co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, gpp_daily, &
1307	deepC_a, deepC_s, deepC_p, &
1308	dt_days, age, PFTpresent, senescence, when_growthinit,&
1309	everywhere, leaf_frac, lm_lastyearmax, npp_longterm,&
1310	carbon_save,deepC_a_save,deepC_s_save,deepC_p_save,delta_fsave,liqwt_max_lastyear)
1311	update_peatfrac = .FALSE.
1312	done_update_peatfrac = .TRUE.
1313	ELSE
1314	![chaoyue] veget_max_tmp is used as veget_max_old in cover SUBROUTINE
1315	CALL cover (npts, cn_ind, ind, biomass, &
1316	veget_max, veget_max_tmp, lai, &
1317	litter, litter_avail, litter_not_avail, carbon, &
1318	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1319	turnover_daily, bm_to_litter, &
1320	co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, gpp_daily, &
1321	deepC_a, deepC_s,deepC_p)
1322	ENDIF
1323
1324	!! 14. Update litter pools to account for harvest
1325
1326	! the whole litter stuff:
1327	! litter update, lignin content, PFT parts, litter decay,
1328	! litter heterotrophic respiration, dead leaf soil cover.
1329	! No vertical discretisation in the soil for litter decay.\n
1330	! added by shilong for harvest
1331	IF(harvest_agri) THEN !!! DO NOT activate harves_agri if ok_LAIdev
1332	CALL harvest(npts, dt_days, veget_max, &
1333	bm_to_litter, turnover_daily, &
1334	harvest_above)
1335	ENDIF
1336
1337	!! 15. Land cover change
1338
1339	IF(EndOfYear) THEN
1340	IF (use_age_class) THEN
1341	CALL age_class_distr(npts, lalo, resolution, bound_spa, &
1342	biomass, veget_max, ind, &
1343	lm_lastyearmax, leaf_frac, co2_to_bm, &
1344	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
1345	everywhere, litter, carbon, lignin_struc, &
1346	deepC_a, deepC_s, deepC_p, &
1347	bm_to_litter, PFTpresent, when_growthinit,&
1348	senescence, npp_longterm, gpp_daily, leaf_age, &
1349	gdd_from_growthinit, gdd_midwinter, time_hum_min, gdd_m5_dormance, &
1350	ncd_dormance, moiavail_month, moiavail_week, ngd_minus5, &
1351	gpp_week, resp_maint, resp_growth, npp_daily)
1352	ENDIF
1353	ENDIF
1354
1355	IF (do_now_stomate_lcchange) THEN
1356	! Initialize LUC specific variables
1357	convflux(:,:) = zero
1358	cflux_prod10(:,:) = zero
1359	cflux_prod100(:,:) = zero
1360
1361	prod10(:,0,:) = zero
1362	prod100(:,0,:) = zero
1363
1364	IF (use_age_class) THEN
1365	IF (SingleAgeClass) THEN
1366	CALL gross_glcchange_SinAgeC_fh (npts, dt_days, harvest_matrix, &
1367	glccSecondShift,glccPrimaryShift,glccNetLCC,&
1368	def_fuel_1hr_remain, def_fuel_10hr_remain, &
1369	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
1370	deforest_litter_remain, deforest_biomass_remain, &
1371	convflux, cflux_prod10, cflux_prod100, &
1372	glccReal, IncreDeficit, glcc_pft, glcc_pftmtc, &
1373	veget_max, prod10, prod100, flux10, flux100, &
1374	PFTpresent, senescence, moiavail_month, moiavail_week, &
1375	gpp_week, ngd_minus5, resp_maint, resp_growth, &
1376	resp_hetero, npp_daily, when_growthinit, npp_longterm, &
1377	ind, lm_lastyearmax, everywhere, age, &
1378	co2_to_bm, gpp_daily, co2_fire, &
1379	time_hum_min, gdd_midwinter, gdd_from_growthinit, &
1380	gdd_m5_dormance, ncd_dormance, &
1381	lignin_struc, carbon, leaf_frac, &
1382	deepC_a, deepC_s, deepC_p, &
1383	leaf_age, bm_to_litter, biomass, litter, &
1384	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
1385
1386	! Multiple age class + wood harvest
1387	ELSE
1388	CALL gross_glcchange_fh (npts, dt_days, harvest_matrix, &
1389	glccSecondShift,glccPrimaryShift,glccNetLCC,&
1390	def_fuel_1hr_remain, def_fuel_10hr_remain, &
1391	def_fuel_100hr_remain, def_fuel_1000hr_remain, &
1392	deforest_litter_remain, deforest_biomass_remain, &
1393	convflux, cflux_prod10, cflux_prod100, &
1394	glccReal, IncreDeficit, glcc_pft, glcc_pftmtc, &
1395	veget_max, prod10, prod100, flux10, flux100, &
1396	PFTpresent, senescence, moiavail_month, moiavail_week, &
1397	gpp_week, ngd_minus5, resp_maint, resp_growth, &
1398	resp_hetero, npp_daily, when_growthinit, npp_longterm, &
1399	ind, lm_lastyearmax, everywhere, age, &
1400	co2_to_bm, gpp_daily, co2_fire, &
1401	time_hum_min, gdd_midwinter, gdd_from_growthinit, &
1402	gdd_m5_dormance, ncd_dormance, &
1403	lignin_struc, carbon, leaf_frac, &
1404	deepC_a, deepC_s, deepC_p, &
1405	leaf_age, bm_to_litter, biomass, litter, &
1406	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr)
1407	ENDIF !(SingleAgeClass)
1408
1409	ELSE ! .NOT. use_age_class
1410	IF (allow_deforest_fire) THEN
1411	![chaoyue] veget_max is used as the old veget_max in lcchange_deffire
1412	! veget_max_new is used as the new veget_max in lcchange_deffire
1413	CALL lcchange_deffire (npts, dt_days, veget_max, veget_max_new, &
1414	biomass, ind, age, PFTpresent, senescence, when_growthinit, everywhere, &
1415	co2_to_bm, bm_to_litter, turnover_daily, bm_sapl, cn_ind, &
1416	flux10(:,:,iwplcc),flux100(:,:,iwplcc), &
1417	prod10(:,:,iwplcc),prod100(:,:,iwplcc),&
1418	convflux(:,iwplcc),&
1419	cflux_prod10(:,iwplcc),cflux_prod100(:,iwplcc), leaf_frac,&
1420	npp_longterm, lm_lastyearmax, litter, litter_avail, litter_not_avail, &
1421	carbon,&
1422	deepC_a, deepC_s, deepC_p,&
1423	fuel_1hr,fuel_10hr,fuel_100hr,fuel_1000hr,&
1424	lcc,bafrac_deforest_accu,emideforest_litter_accu,emideforest_biomass_accu,&
1425	deflitsup_total,defbiosup_total)
1426	ELSE
1427	![chaoyue] veget_max is used as veget_max_old in lcchange_main
1428	! veget_max_new is used as the new veget_max in lcchange_main
1429	!!!qcj++ peatland
1430	IF (agri_peat) THEN
1431	CALL lcchange_main_agripeat (npts, dt_days, veget_max_new, veget_max, &
1432	biomass, ind, age, PFTpresent, senescence, when_growthinit, everywhere, &
1433	co2_to_bm, bm_to_litter, turnover_daily, bm_sapl, cn_ind,&
1434	flux10(:,:,iwplcc),flux100(:,:,iwplcc), &
1435	prod10(:,:,iwplcc),prod100(:,:,iwplcc),convflux(:,iwplcc), &
1436	cflux_prod10(:,iwplcc),cflux_prod100(:,iwplcc),leaf_frac,&
1437	npp_longterm, lm_lastyearmax, litter, litter_avail, litter_not_avail, &
1438	carbon, &
1439	deepC_a, deepC_s, deepC_p,&
1440	fuel_1hr,fuel_10hr,fuel_100hr,fuel_1000hr, &
1441	veget_max_adjusted,lalo,carbon_save,deepC_a_save,deepC_s_save,deepC_p_save,delta_fsave,biomass_remove)
1442	ELSE
1443	IF (.NOT. dyn_peat) THEN
1444	CALL lcchange_main (npts, dt_days, veget_max_new, veget_max, &
1445	biomass, ind, age, PFTpresent, senescence, when_growthinit, everywhere, &
1446	co2_to_bm, bm_to_litter, turnover_daily, bm_sapl, cn_ind,&
1447	flux10(:,:,iwplcc),flux100(:,:,iwplcc), &
1448	prod10(:,:,iwplcc),prod100(:,:,iwplcc),convflux(:,iwplcc), &
1449	cflux_prod10(:,iwplcc),cflux_prod100(:,iwplcc),leaf_frac,&
1450	npp_longterm, lm_lastyearmax, litter, litter_avail, litter_not_avail, &
1451	carbon, &
1452	deepC_a, deepC_s, deepC_p,&
1453	fuel_1hr,fuel_10hr,fuel_100hr,fuel_1000hr)
1454	ENDIF
1455	ENDIF
1456	ENDIF
1457	ENDIF ! (use_age_class)
1458
1459	do_now_stomate_lcchange=.FALSE.
1460
1461	! Set the flag done_stomate_lcchange to be used in the end of sechiba_main to update the fractions.
1462	done_stomate_lcchange=.TRUE.
1463
1464	! [Added by chaoyue]
1465	! As it's possbile people may forget to check if their land use change
1466	! data (either gross or net land cover change allows the conservation of
1467	! sum of veget_cov_max (i.e., equal to 1 or
1468	! zero), so will adjust them here in case this happens.
1469
1470	DO ipts = 1,npts
1471	valtmp = SUM(veget_max(ipts,:))
1472	IF (valtmp .GT. min_stomate) THEN
1473	DO j = 1,nvm
1474	veget_max(ipts,j) = veget_max(ipts,j)/valtmp
1475	ENDDO
1476	ELSE
1477	veget_max(ipts,:) = zero
1478	ENDIF
1479	ENDDO
1480
1481	ENDIF ! do_now_stomate_lcchange
1482
1483	!MM dÃ©placement pour initialisation correcte des grandeurs cumulÃ©es :
1484	cflux_prod_total(:) = SUM(convflux(:,:) + cflux_prod10(:,:) + cflux_prod100(:,:),DIM=2)
1485	prod10_total(:)=SUM(SUM(prod10,dim=2),DIM=2)
1486	prod100_total(:)=SUM(SUM(prod100,dim=2),DIM=2)
1487
1488	!! 16. Recalculate lai
1489	!! This should be done whenever biomass is modified
1490	CALL setlai(biomass, sla_calc, slai)
1491	CALL setlai(biomass, sla_calc, lai)
1492
1493	!! 17. Calculate vcmax
1494	CALL vmax (npts, dt_days, &
1495	leaf_age, leaf_frac, &
1496	vcmax, &!)
1497	N_limfert)
1498
1499	!! 18. Total heterotrophic respiration
1500
1501	tot_soil_carb(:,:) = zero
1502	tot_litter_carb(:,:) = zero
1503	DO j=1,nvm
1504
1505	tot_litter_carb(:,j) = tot_litter_carb(:,j) + (litter(:,istructural,j,iabove,icarbon) + &
1506	& litter(:,imetabolic,j,iabove,icarbon) + &
1507	& litter(:,istructural,j,ibelow,icarbon) + litter(:,imetabolic,j,ibelow,icarbon))
1508
1509	tot_soil_carb(:,j) = tot_soil_carb(:,j) + (carbon(:,iactive,j) + &
1510	& carbon(:,islow,j)+ carbon(:,ipassive,j))
1511
1512	ENDDO
1513	tot_litter_soil_carb(:,:) = tot_litter_carb(:,:) + tot_soil_carb(:,:)
1514
1515	!!$ DO k = 1, nelements ! Loop over # elements
1516	!!$ tot_live_biomass(:,:,k) = biomass(:,:,ileaf,k) + biomass(:,:,isapabove,k) + biomass(:,:,isapbelow,k) +&
1517	!!$ & biomass(:,:,iheartabove,k) + biomass(:,:,iheartbelow,k) + &
1518	!!$ & biomass(:,:,iroot,k) + biomass(:,:,ifruit,k) + biomass(:,:,icarbres,k)
1519	!!$ END DO ! Loop over # elements
1520
1521	tot_live_biomass(:,:,:) = biomass(:,:,ileaf,:) + biomass(:,:,isapabove,:) + biomass(:,:,isapbelow,:) +&
1522	& biomass(:,:,iheartabove,:) + biomass(:,:,iheartbelow,:) + &
1523	& biomass(:,:,iroot,:) + biomass(:,:,ifruit,:) + biomass(:,:,icarbres,:)
1524
1525
1526	tot_turnover(:,:,:) = turnover_daily(:,:,ileaf,:) + turnover_daily(:,:,isapabove,:) + &
1527	& turnover_daily(:,:,isapbelow,:) + turnover_daily(:,:,iheartabove,:) + &
1528	& turnover_daily(:,:,iheartbelow,:) + turnover_daily(:,:,iroot,:) + &
1529	& turnover_daily(:,:,ifruit,:) + turnover_daily(:,:,icarbres,:)
1530
1531	tot_bm_to_litter(:,:,:) = bm_to_litter(:,:,ileaf,:) + bm_to_litter(:,:,isapabove,:) +&
1532	& bm_to_litter(:,:,isapbelow,:) + bm_to_litter(:,:,iheartbelow,:) +&
1533	& bm_to_litter(:,:,iheartabove,:) + bm_to_litter(:,:,iroot,:) + &
1534	& bm_to_litter(:,:,ifruit,:) + bm_to_litter(:,:,icarbres,:)
1535
1536	carb_mass_variation(:)=-carb_mass_total(:)
1537	carb_mass_total(:)=SUM((tot_live_biomass(:,:,icarbon)+tot_litter_carb+tot_soil_carb)*veget_max,dim=2) + &
1538	& (prod10_total + prod100_total)
1539	carb_mass_variation(:)=carb_mass_total(:)+carb_mass_variation(:)
1540
1541	!! 17. Write history
1542
1543	CALL xios_orchidee_send_field("RESOLUTION_X",resolution(:,1))
1544	CALL xios_orchidee_send_field("RESOLUTION_Y",resolution(:,2))
1545	CALL xios_orchidee_send_field("CONTFRAC_STOMATE",contfrac(:))
1546	CALL xios_orchidee_send_field("T2M_MONTH",t2m_month)
1547	CALL xios_orchidee_send_field("T2M_WEEK",t2m_week)
1548	CALL xios_orchidee_send_field("HET_RESP",resp_hetero(:,:))
1549	CALL xios_orchidee_send_field("CO2_FIRE",co2_fire)
1550	CALL xios_orchidee_send_field("CO2_TAKEN",co2_to_bm)
1551	CALL xios_orchidee_send_field("LAI",lai)
1552	! DO ipts = 1,npts
1553	! WRITE (numout,*) 'chunjing end of stomate_lpj',veget_max(ipts,:)
1554	! ENDDO
1555	CALL xios_orchidee_send_field("VEGET_MAX",veget_max)
1556	CALL xios_orchidee_send_field("NPP_STOMATE",npp_daily)
1557	CALL xios_orchidee_send_field("GPP",gpp_daily)
1558	CALL xios_orchidee_send_field("IND",ind)
1559	CALL xios_orchidee_send_field("CN_IND",cn_ind)
1560	CALL xios_orchidee_send_field("WOODMASS_IND",woodmass_ind)
1561	! CALL xios_orchidee_send_field("TOTAL_M",tot_live_biomass(:,:,icarbon))
1562	CALL xios_orchidee_send_field("MOISTRESS",moiavail_week)
1563	CALL xios_orchidee_send_field("LEAF_M",biomass(:,:,ileaf,icarbon))
1564	CALL xios_orchidee_send_field("SAP_M_AB",biomass(:,:,isapabove,icarbon))
1565	CALL xios_orchidee_send_field("SAP_M_BE",biomass(:,:,isapbelow,icarbon))
1566	CALL xios_orchidee_send_field("HEART_M_AB",biomass(:,:,iheartabove,icarbon))
1567	CALL xios_orchidee_send_field("HEART_M_BE",biomass(:,:,iheartbelow,icarbon))
1568	CALL xios_orchidee_send_field("ROOT_M",biomass(:,:,iroot,icarbon))
1569	CALL xios_orchidee_send_field("FRUIT_M",biomass(:,:,ifruit,icarbon))
1570	! HERE WE ADD A OUTPUT VARIABLE, NAMED CROPYIELD
1571	CALL xios_orchidee_send_field("CROPYIELD",biomass(:,:,ifruit,icarbon))
1572	CALL xios_orchidee_send_field("BIOMYIELD",tot_live_biomass(:,:,icarbon))
1573	! END DEFINE
1574	CALL xios_orchidee_send_field("RESERVE_M",biomass(:,:,icarbres,icarbon))
1575	! CALL xios_orchidee_send_field("TOTAL_TURN",tot_turnover)
1576	CALL xios_orchidee_send_field("LEAF_TURN",turnover_daily(:,:,ileaf,icarbon))
1577	CALL xios_orchidee_send_field("MAINT_RESP",resp_maint)
1578	CALL xios_orchidee_send_field("GROWTH_RESP",resp_growth)
1579	CALL xios_orchidee_send_field("SAP_AB_TURN",turnover_daily(:,:,isapabove,icarbon))
1580	CALL xios_orchidee_send_field("ROOT_TURN",turnover_daily(:,:,iroot,icarbon))
1581	CALL xios_orchidee_send_field("FRUIT_TURN",turnover_daily(:,:,ifruit,icarbon))
1582	CALL xios_orchidee_send_field("TOTAL_BM_LITTER",tot_bm_to_litter(:,:,icarbon))
1583	CALL xios_orchidee_send_field("LEAF_BM_LITTER",bm_to_litter(:,:,ileaf,icarbon))
1584	CALL xios_orchidee_send_field("SAP_AB_BM_LITTER",bm_to_litter(:,:,isapabove,icarbon))
1585	CALL xios_orchidee_send_field("SAP_BE_BM_LITTER",bm_to_litter(:,:,isapbelow,icarbon))
1586	CALL xios_orchidee_send_field("HEART_AB_BM_LITTER",bm_to_litter(:,:,iheartabove,icarbon))
1587	CALL xios_orchidee_send_field("HEART_BE_BM_LITTER",bm_to_litter(:,:,iheartbelow,icarbon))
1588	CALL xios_orchidee_send_field("ROOT_BM_LITTER",bm_to_litter(:,:,iroot,icarbon))
1589	CALL xios_orchidee_send_field("FRUIT_BM_LITTER",bm_to_litter(:,:,ifruit,icarbon))
1590	CALL xios_orchidee_send_field("RESERVE_BM_LITTER",bm_to_litter(:,:,icarbres,icarbon))
1591	CALL xios_orchidee_send_field("LITTER_STR_AB",litter(:,istructural,:,iabove,icarbon))
1592	CALL xios_orchidee_send_field("LITTER_MET_AB",litter(:,imetabolic,:,iabove,icarbon))
1593	CALL xios_orchidee_send_field("LITTER_STR_BE",litter(:,istructural,:,ibelow,icarbon))
1594	CALL xios_orchidee_send_field("LITTER_MET_BE",litter(:,imetabolic,:,ibelow,icarbon))
1595	CALL xios_orchidee_send_field("DEADLEAF_COVER",deadleaf_cover)
1596	CALL xios_orchidee_send_field("TOTAL_SOIL_CARB",tot_litter_soil_carb)
1597	CALL xios_orchidee_send_field("CARBON_ACTIVE",carbon(:,iactive,:))
1598	CALL xios_orchidee_send_field("CARBON_SLOW",carbon(:,islow,:))
1599	CALL xios_orchidee_send_field("CARBON_PASSIVE",carbon(:,ipassive,:))
1600	CALL xios_orchidee_send_field("LITTERHUM",litterhum_daily)
1601	CALL xios_orchidee_send_field("TURNOVER_TIME",turnover_time)
1602	! CALL xios_orchidee_send_field("PROD10",prod10)
1603	! CALL xios_orchidee_send_field("FLUX10",flux10)
1604	! CALL xios_orchidee_send_field("PROD100",prod100)
1605	! CALL xios_orchidee_send_field("FLUX100",flux100)
1606	CALL xios_orchidee_send_field("CONVFLUX",convflux)
1607	CALL xios_orchidee_send_field("CFLUX_PROD10",cflux_prod10)
1608	CALL xios_orchidee_send_field("CFLUX_PROD100",cflux_prod100)
1609	CALL xios_orchidee_send_field("HARVEST_ABOVE",harvest_above)
1610	CALL xios_orchidee_send_field("VCMAX",vcmax)
1611	CALL xios_orchidee_send_field("AGE",age)
1612	CALL xios_orchidee_send_field("HEIGHT",height)
1613	CALL xios_orchidee_send_field("FIREINDEX",fireindex(:,:))
1614	!gmjc
1615	CALL xios_orchidee_send_field("LITTER_STR_AVAIL",litter_avail(:,istructural,:))
1616	CALL xios_orchidee_send_field("LITTER_MET_AVAIL",litter_avail(:,imetabolic,:))
1617	CALL xios_orchidee_send_field("LITTER_STR_NAVAIL",litter_not_avail(:,istructural,:))
1618	CALL xios_orchidee_send_field("LITTER_MET_NAVAIL",litter_not_avail(:,imetabolic,:))
1619	CALL xios_orchidee_send_field("LITTER_STR_AVAILF",litter_avail_frac(:,istructural,:))
1620	CALL xios_orchidee_send_field("LITTER_MET_AVAILF",litter_avail_frac(:,imetabolic,:))
1621	CALL xios_orchidee_send_field("N_LIMFERT",N_limfert)
1622	!!!qcj++ peatland
1623	IF (ok_peat) THEN
1624	CALL xios_orchidee_send_field("TCARBON_ACRO", tcarbon_acro)
1625	CALL xios_orchidee_send_field("TCARBON_CATO", tcarbon_cato)
1626	CALL xios_orchidee_send_field("CARBON_ACRO",carbon_acro(:,:))
1627	CALL xios_orchidee_send_field("CARBON_CATO",carbon_cato(:,:))
1628	CALL xios_orchidee_send_field("HEIGHT_ACRO",height_acro)
1629	CALL xios_orchidee_send_field("HEIGHT_CATO",height_cato)
1630	CALL xios_orchidee_send_field("RESP_ACRO_oxic",resp_acro_oxic_d(:,:))
1631	CALL xios_orchidee_send_field("RESP_ACRO_anoxic",resp_acro_anoxic_d(:,:))
1632	CALL xios_orchidee_send_field("RESP_CATO",resp_cato_d(:,:))
1633	CALL xios_orchidee_send_field("LITTER_to_ACRO",litter_to_acro_d(:,:))
1634	CALL xios_orchidee_send_field("ACRO_to_CATO",acro_to_cato_d)
1635	ENDIF
1636	CALL xios_orchidee_send_field("deepC_a_save",deepC_a_save)
1637	CALL xios_orchidee_send_field("deepC_s_save",deepC_s_save)
1638	CALL xios_orchidee_send_field("deepC_p_save",deepC_p_save)
1639	CALL xios_orchidee_send_field("delta_fsave",delta_fsave)
1640	CALL xios_orchidee_send_field("carbon_a_save",carbon_save(:,iactive,:))
1641	CALL xios_orchidee_send_field("carbon_s_save",carbon_save(:,islow,:))
1642	CALL xios_orchidee_send_field("carbon_p_save",carbon_save(:,ipassive,:))
1643	CALL xios_orchidee_send_field("biomass_remove_leaf",biomass_remove(:,:,ileaf,icarbon))
1644	CALL xios_orchidee_send_field("biomass_remove_sapabove",biomass_remove(:,:,isapabove,icarbon))
1645	CALL xios_orchidee_send_field("biomass_remove_sapbelow",biomass_remove(:,:,isapbelow,icarbon))
1646	CALL xios_orchidee_send_field("biomass_remove_heartabove",biomass_remove(:,:,iheartabove,icarbon))
1647	CALL xios_orchidee_send_field("biomass_remove_heartbelow",biomass_remove(:,:,iheartbelow,icarbon))
1648	CALL xios_orchidee_send_field("biomass_remove_root",biomass_remove(:,:,iroot,icarbon))
1649	CALL xios_orchidee_send_field("biomass_remove_fruit",biomass_remove(:,:,ifruit,icarbon))
1650	CALL xios_orchidee_send_field("biomass_remove_carbres",biomass_remove(:,:,icarbres,icarbon))
1651
1652	!calculate grassland co2 fluxes
1653	DO j=2,nvm
1654	IF ((.NOT. is_tree(j)) .AND. natural(j)) THEN
1655	veget_max_gm(:,j) = veget_max(:,j)
1656	ENDIF
1657	END DO ! nvm
1658	veget_exist_gm(:) = SUM(veget_max_gm,dim=2)
1659	WHERE (veget_exist_gm(:) .GT. 0.0)
1660	co2_gm(:) = SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire &
1661	-harvest_gm-ranimal_gm-ch4_pft_gm+cinput_gm) &
1662	*veget_max_gm,dim=2)/veget_exist_gm
1663	ch4_gm(:) = SUM(ch4_pft_gm*veget_max_gm,dim=2)/veget_exist_gm
1664	n2o_gm(:) = SUM(n2o_pft_gm*veget_max_gm,dim=2)/veget_exist_gm
1665	ELSEWHERE
1666	co2_gm(:) = zero
1667	ch4_gm(:) = zero
1668	n2o_gm(:) = zero
1669	ENDWHERE
1670	CALL xios_orchidee_send_field("CO2_GM",co2_gm)
1671	CALL xios_orchidee_send_field("CH4_GM",ch4_gm)
1672	CALL xios_orchidee_send_field("N2O_GM",n2o_gm)
1673	CALL xios_orchidee_send_field("N2O_PFT_GM",n2o_pft_gm)
1674	!end gmjc
1675
1676	! ipcc history
1677	CALL xios_orchidee_send_field("cVeg",SUM(tot_live_biomass(:,:,icarbon)veget_max,dim=2)/1e3contfrac)
1678	CALL xios_orchidee_send_field("cLitter",SUM(tot_litter_carbveget_max,dim=2)/1e3contfrac)
1679	CALL xios_orchidee_send_field("cSoil",SUM(tot_soil_carbveget_max,dim=2)/1e3contfrac)
1680	CALL xios_orchidee_send_field("cProduct",(prod10_total + prod100_total)/1e3)
1681	CALL xios_orchidee_send_field("cMassVariation",carb_mass_variation/1e3/one_day*contfrac)
1682	CALL xios_orchidee_send_field("lai_ipcc",SUM(laiveget_max,dim=2)contfrac)
1683	CALL xios_orchidee_send_field("gpp_ipcc",SUM(gpp_dailyveget_max,dim=2)/1e3/one_daycontfrac)
1684	CALL xios_orchidee_send_field("ra",SUM((resp_maint+resp_growth)veget_max,dim=2)/1e3/one_daycontfrac)
1685	CALL xios_orchidee_send_field("npp_ipcc",SUM(npp_dailyveget_max,dim=2)/1e3/one_daycontfrac)
1686	CALL xios_orchidee_send_field("rh",SUM(resp_heteroveget_max,dim=2)/1e3/one_daycontfrac)
1687	CALL xios_orchidee_send_field("fFire",SUM(co2_fireveget_max,dim=2)/1e3/one_daycontfrac)
1688	CALL xios_orchidee_send_field("fHarvest",harvest_above/1e3/one_day*contfrac)
1689	CALL xios_orchidee_send_field("fLuc",cflux_prod_total/1e3/one_day*contfrac)
1690	!gmjc
1691	IF (enable_grazing) THEN
1692	CALL xios_orchidee_send_field("nbp",&
1693	(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire &
1694	-harvest_gm-ranimal_gm-ch4_pft_gm+cinput_gm) & ! specific items for gmjc
1695	& veget_max,dim=2)-cflux_prod_total-harvest_above)/1e3/one_daycontfrac)
1696	ELSE
1697	CALL xios_orchidee_send_field("nbp",&
1698	(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire) &
1699	& veget_max,dim=2)-cflux_prod_total-harvest_above)/1e3/one_daycontfrac)
1700	ENDIF
1701	!end gmjc
1702	! CALL xios_orchidee_send_field("nbp",(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire) &
1703	! & veget_max,dim=2)-cflux_prod_total-harvest_above)/1e3/one_daycontfrac)
1704
1705	CALL xios_orchidee_send_field("fVegLitter",SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))*&
1706	veget_max,dim=2)/1e3/one_day*contfrac)
1707	CALL xios_orchidee_send_field("fLitterSoil",SUM(SUM(soilcarbon_input,dim=2)veget_max,dim=2)/1e3/one_daycontfrac)
1708	CALL xios_orchidee_send_field("cLeaf",SUM(biomass(:,:,ileaf,icarbon)veget_max,dim=2)/1e3contfrac)
1709	CALL xios_orchidee_send_field("cWood",SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon))*&
1710	veget_max,dim=2)/1e3*contfrac)
1711	CALL xios_orchidee_send_field("cRoot",SUM(( biomass(:,:,iroot,icarbon) + biomass(:,:,isapbelow,icarbon) + &
1712	biomass(:,:,iheartbelow,icarbon) )veget_max,dim=2)/1e3contfrac)
1713	CALL xios_orchidee_send_field("cMisc",SUM(( biomass(:,:,icarbres,icarbon) + biomass(:,:,ifruit,icarbon))*&
1714	veget_max,dim=2)/1e3*contfrac)
1715	CALL xios_orchidee_send_field("cLitterAbove",SUM((litter(:,istructural,:,iabove,icarbon)+&
1716	litter(:,imetabolic,:,iabove,icarbon))veget_max,dim=2)/1e3contfrac)
1717	CALL xios_orchidee_send_field("cLitterBelow",SUM((litter(:,istructural,:,ibelow,icarbon)+&
1718	litter(:,imetabolic,:,ibelow,icarbon))veget_max,dim=2)/1e3contfrac)
1719	CALL xios_orchidee_send_field("cSoilFast",SUM(carbon(:,iactive,:)veget_max,dim=2)/1e3contfrac)
1720	CALL xios_orchidee_send_field("cSoilMedium",SUM(carbon(:,islow,:)veget_max,dim=2)/1e3contfrac)
1721	CALL xios_orchidee_send_field("cSoilSlow",SUM(carbon(:,ipassive,:)veget_max,dim=2)/1e3contfrac)
1722	DO j=1,nvm
1723	histvar(:,j)=veget_max(:,j)contfrac(:)100
1724	ENDDO
1725	CALL xios_orchidee_send_field("landCoverFrac",histvar)
1726	vartmp(:)=zero
1727	DO j = 2,nvm
1728	IF (is_deciduous(j)) THEN
1729	vartmp(:) = vartmp(:) + veget_max(:,j)contfrac100
1730	ENDIF
1731	ENDDO
1732	CALL xios_orchidee_send_field("treeFracPrimDec",vartmp)
1733	vartmp(:)=zero
1734	DO j = 2,nvm
1735	IF (is_evergreen(j)) THEN
1736	vartmp(:) = vartmp(:) + veget_max(:,j)contfrac100
1737	ENDIF
1738	ENDDO
1739	CALL xios_orchidee_send_field("treeFracPrimEver",vartmp)
1740	vartmp(:)=zero
1741	DO j = 2,nvm
1742	IF ( .NOT.(is_c4(j)) ) THEN
1743	vartmp(:) = vartmp(:) + veget_max(:,j)contfrac100
1744	ENDIF
1745	ENDDO
1746	CALL xios_orchidee_send_field("c3PftFrac",vartmp)
1747	vartmp(:)=zero
1748	DO j = 2,nvm
1749	IF ( is_c4(j) ) THEN
1750	vartmp(:) = vartmp(:) + veget_max(:,j)contfrac100
1751	ENDIF
1752	ENDDO
1753	CALL xios_orchidee_send_field("c4PftFrac",vartmp)
1754	CALL xios_orchidee_send_field("rGrowth",SUM(resp_growthveget_max,dim=2)/1e3/one_daycontfrac)
1755	CALL xios_orchidee_send_field("rMaint",SUM(resp_maintveget_max,dim=2)/1e3/one_daycontfrac)
1756	CALL xios_orchidee_send_field("nppLeaf",SUM(bm_alloc(:,:,ileaf,icarbon)veget_max,dim=2)/1e3/one_daycontfrac)
1757	CALL xios_orchidee_send_field("nppWood",SUM(bm_alloc(:,:,isapabove,icarbon)veget_max,dim=2)/1e3/one_daycontfrac)
1758	CALL xios_orchidee_send_field("nppRoot",SUM(( bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iroot,icarbon) )*&
1759	veget_max,dim=2)/1e3/one_day*contfrac)
1760
1761
1762	CALL histwrite_p (hist_id_stomate, 'RESOLUTION_X', itime, &
1763	resolution(:,1), npts, hori_index)
1764	CALL histwrite_p (hist_id_stomate, 'RESOLUTION_Y', itime, &
1765	resolution(:,2), npts, hori_index)
1766	CALL histwrite_p (hist_id_stomate, 'CONTFRAC', itime, &
1767	contfrac(:), npts, hori_index)
1768
1769	CALL histwrite_p (hist_id_stomate, 'LITTER_STR_AB', itime, &
1770	litter(:,istructural,:,iabove,icarbon), npts*nvm, horipft_index)
1771	CALL histwrite_p (hist_id_stomate, 'LITTER_MET_AB', itime, &
1772	litter(:,imetabolic,:,iabove,icarbon), npts*nvm, horipft_index)
1773	CALL histwrite_p (hist_id_stomate, 'LITTER_STR_BE', itime, &
1774	litter(:,istructural,:,ibelow,icarbon), npts*nvm, horipft_index)
1775	CALL histwrite_p (hist_id_stomate, 'LITTER_MET_BE', itime, &
1776	litter(:,imetabolic,:,ibelow,icarbon), npts*nvm, horipft_index)
1777	!spitfiretest
1778	!CALL histwrite (hist_id_stomate, 'fuel_1hr_met_b', itime, &
1779	! fuel_1hr(:,:,imetabolic,icarbon), npts*nvm, horipft_index)
1780	!CALL histwrite (hist_id_stomate, 'fuel_1hr_str_b', itime, &
1781	! fuel_1hr(:,:,istructural,icarbon), npts*nvm, horipft_index)
1782	!CALL histwrite (hist_id_stomate, 'fuel_10hr_met_b', itime, &
1783	! fuel_10hr(:,:,imetabolic,icarbon), npts*nvm, horipft_index)
1784	!CALL histwrite (hist_id_stomate, 'fuel_10hr_str_b', itime, &
1785	! fuel_10hr(:,:,istructural,icarbon), npts*nvm, horipft_index)
1786	!CALL histwrite (hist_id_stomate, 'fuel_100hr_met_b', itime, &
1787	! fuel_100hr(:,:,imetabolic,icarbon), npts*nvm, horipft_index)
1788	!CALL histwrite (hist_id_stomate, 'fuel_100hr_str_b', itime, &
1789	! fuel_100hr(:,:,istructural,icarbon), npts*nvm, horipft_index)
1790	!CALL histwrite (hist_id_stomate, 'fuel_1000hr_met_b', itime, &
1791	! fuel_1000hr(:,:,imetabolic,icarbon), npts*nvm, horipft_index)
1792	!CALL histwrite (hist_id_stomate, 'fuel_1000hr_str_b', itime, &
1793	! fuel_1000hr(:,:,istructural,icarbon), npts*nvm, horipft_index)
1794	!endspittest
1795
1796	CALL histwrite_p (hist_id_stomate, 'DEADLEAF_COVER', itime, &
1797	deadleaf_cover, npts, hori_index)
1798
1799	CALL histwrite_p (hist_id_stomate, 'TOTAL_SOIL_CARB', itime, &
1800	tot_litter_soil_carb, npts*nvm, horipft_index)
1801	CALL histwrite_p (hist_id_stomate, 'CARBON_ACTIVE', itime, &
1802	carbon(:,iactive,:), npts*nvm, horipft_index)
1803	CALL histwrite_p (hist_id_stomate, 'CARBON_SLOW', itime, &
1804	carbon(:,islow,:), npts*nvm, horipft_index)
1805	CALL histwrite_p (hist_id_stomate, 'CARBON_PASSIVE', itime, &
1806	carbon(:,ipassive,:), npts*nvm, horipft_index)
1807
1808	CALL histwrite_p (hist_id_stomate, 'CARBON_ACTIVE_SURF', itime, &
1809	carbon_surf(:,iactive,:), npts*nvm, horipft_index)
1810	CALL histwrite_p (hist_id_stomate, 'CARBON_SLOW_SURF', itime, &
1811	carbon_surf(:,islow,:), npts*nvm, horipft_index)
1812	CALL histwrite_p (hist_id_stomate, 'CARBON_PASSIVE_SURF', itime, &
1813	carbon_surf(:,ipassive,:), npts*nvm, horipft_index)
1814
1815	!!!! Wetland CH4 methane
1816	!pss:+
1817	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_TOT_0', itime, &
1818	ch4_flux_density_tot_0, npts, hori_index)
1819	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_DIF_0', itime, &
1820	ch4_flux_density_dif_0, npts, hori_index)
1821	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_BUB_0', itime, &
1822	ch4_flux_density_bub_0, npts, hori_index)
1823	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_PLA_0', itime, &
1824	ch4_flux_density_pla_0, npts, hori_index)
1825	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_TOT_wet1', itime, &
1826	ch4_flux_density_tot_wet1, npts, hori_index)
1827	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_DIF_wet1', itime, &
1828	ch4_flux_density_dif_wet1, npts, hori_index)
1829	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_BUB_wet1', itime, &
1830	ch4_flux_density_bub_wet1, npts, hori_index)
1831	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_PLA_wet1', itime, &
1832	ch4_flux_density_pla_wet1, npts, hori_index)
1833	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_TOT_wet2', itime, &
1834	ch4_flux_density_tot_wet2, npts, hori_index)
1835	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_DIF_wet2', itime, &
1836	ch4_flux_density_dif_wet2, npts, hori_index)
1837	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_BUB_wet2', itime, &
1838	ch4_flux_density_bub_wet2, npts, hori_index)
1839	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_PLA_wet2', itime, &
1840	ch4_flux_density_pla_wet2, npts, hori_index)
1841	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_TOT_wet3', itime, &
1842	ch4_flux_density_tot_wet3, npts, hori_index)
1843	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_DIF_wet3', itime, &
1844	ch4_flux_density_dif_wet3, npts, hori_index)
1845	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_BUB_wet3', itime, &
1846	ch4_flux_density_bub_wet3, npts, hori_index)
1847	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_PLA_wet3', itime, &
1848	ch4_flux_density_pla_wet3, npts, hori_index)
1849	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_TOT_wet4', itime, &
1850	ch4_flux_density_tot_wet4, npts, hori_index)
1851	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_DIF_wet4', itime, &
1852	ch4_flux_density_dif_wet4, npts, hori_index)
1853	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_BUB_wet4', itime, &
1854	ch4_flux_density_bub_wet4, npts, hori_index)
1855	CALL histwrite_p (hist_id_stomate, 'CH4_FLUX_PLA_wet4', itime, &
1856	ch4_flux_density_pla_wet4, npts, hori_index)
1857
1858	CALL histwrite_p (hist_id_stomate, 'TSURF_YEAR', itime, &
1859	tsurf_year, npts, hori_index)
1860	!pss:-
1861
1862
1863	CALL histwrite_p (hist_id_stomate, 'T2M_MONTH', itime, &
1864	t2m_month, npts, hori_index)
1865	CALL histwrite_p (hist_id_stomate, 'T2M_WEEK', itime, &
1866	t2m_week, npts, hori_index)
1867	CALL histwrite_p (hist_id_stomate, 'TSEASON', itime, &
1868	Tseason, npts, hori_index)
1869	CALL histwrite_p (hist_id_stomate, 'TMIN_SPRING_TIME', itime, &
1870	Tmin_spring_time, npts*nvm, horipft_index)
1871	CALL histwrite_p (hist_id_stomate, 'ONSET_DATE', itime, &
1872	onset_date(:,:), npts*nvm, horipft_index)
1873
1874	CALL histwrite_p (hist_id_stomate, 'HET_RESP', itime, &
1875	resp_hetero(:,:), npts*nvm, horipft_index)
1876
1877	! gmjc
1878	CALL histwrite_p(hist_id_stomate ,'T2M_14' ,itime, &
1879	t2m_14, npts, hori_index)
1880	! CALL histwrite (hist_id_stomate, 'LITTER_RESP', itime, &
1881	! resp_hetero_litter_d(:,:), npts*nvm, horipft_index)
1882	! CALL histwrite (hist_id_stomate, 'ACTIVE_RESP', itime, &
1883	! resp_hetero_soil_d(:,iactive,:), npts*nvm, horipft_index)
1884	! CALL histwrite (hist_id_stomate, 'SLOW_RESP', itime, &
1885	! resp_hetero_soil_d(:,islow,:), npts*nvm, horipft_index)
1886	! CALL histwrite (hist_id_stomate, 'PASSIVE_RESP', itime, &
1887	! resp_hetero_soil_d(:,ipassive,:), npts*nvm, horipft_index)
1888	CALL histwrite_p (hist_id_stomate, 'LITTER_STR_AVAIL', itime, &
1889	litter_avail(:,istructural,:), npts*nvm, horipft_index)
1890	CALL histwrite_p (hist_id_stomate, 'LITTER_MET_AVAIL', itime, &
1891	litter_avail(:,imetabolic,:), npts*nvm, horipft_index)
1892	CALL histwrite_p (hist_id_stomate, 'LITTER_STR_NAVAIL', itime, &
1893	litter_not_avail(:,istructural,:), npts*nvm, horipft_index)
1894	CALL histwrite_p (hist_id_stomate, 'LITTER_MET_NAVAIL', itime, &
1895	litter_not_avail(:,imetabolic,:), npts*nvm, horipft_index)
1896	CALL histwrite_p (hist_id_stomate, 'LITTER_STR_AVAILF', itime, &
1897	litter_avail_frac(:,istructural,:), npts*nvm, horipft_index)
1898	CALL histwrite_p (hist_id_stomate, 'LITTER_MET_AVAILF', itime, &
1899	litter_avail_frac(:,imetabolic,:), npts*nvm, horipft_index)
1900
1901	! Crop is enabled
1902	IF (ANY(ok_LAIdev)) THEN
1903	CALL histwrite_p (hist_id_stomate, 'N_LIMFERT', itime, &
1904	N_limfert, npts*nvm, horipft_index)
1905	ENDIF
1906	! end gmjc
1907	CALL histwrite_p (hist_id_stomate, 'FIREINDEX', itime, &
1908	fireindex(:,:), npts*nvm, horipft_index)
1909	CALL histwrite_p (hist_id_stomate, 'LITTERHUM', itime, &
1910	litterhum_daily, npts, hori_index)
1911	CALL histwrite_p (hist_id_stomate, 'CO2_FIRE', itime, &
1912	co2_fire, npts*nvm, horipft_index)
1913	CALL histwrite_p (hist_id_stomate, 'CO2_TAKEN', itime, &
1914	co2_to_bm, npts*nvm, horipft_index)
1915	! land cover change
1916	CALL histwrite_p (hist_id_stomate, 'CONVFLUX_LCC', itime, &
1917	convflux(:,iwplcc), npts, hori_index)
1918	CALL histwrite_p (hist_id_stomate, 'CONVFLUX_HAR', itime, &
1919	convflux(:,iwphar), npts, hori_index)
1920	CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD10_LCC', itime, &
1921	cflux_prod10(:,iwplcc), npts, hori_index)
1922	CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD10_HAR', itime, &
1923	cflux_prod10(:,iwphar), npts, hori_index)
1924	CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD100_LCC', itime, &
1925	cflux_prod100(:,iwplcc), npts, hori_index)
1926	CALL histwrite_p (hist_id_stomate, 'CFLUX_PROD100_HAR', itime, &
1927	cflux_prod100(:,iwphar), npts, hori_index)
1928	CALL histwrite_p (hist_id_stomate, 'HARVEST_ABOVE', itime, &
1929	harvest_above, npts, hori_index)
1930	CALL histwrite_p (hist_id_stomate, 'PROD10_LCC', itime, &
1931	prod10(:,:,iwplcc), npts*11, horip11_index)
1932	CALL histwrite_p (hist_id_stomate, 'PROD10_HAR', itime, &
1933	prod10(:,:,iwphar), npts*11, horip11_index)
1934	CALL histwrite_p (hist_id_stomate, 'PROD100_LCC', itime, &
1935	prod100(:,:,iwplcc), npts*101, horip101_index)
1936	CALL histwrite_p (hist_id_stomate, 'PROD100_HAR', itime, &
1937	prod100(:,:,iwphar), npts*101, horip101_index)
1938	CALL histwrite_p (hist_id_stomate, 'FLUX10_LCC', itime, &
1939	flux10(:,:,iwplcc), npts*10, horip10_index)
1940	CALL histwrite_p (hist_id_stomate, 'FLUX10_HAR', itime, &
1941	flux10(:,:,iwphar), npts*10, horip10_index)
1942	CALL histwrite_p (hist_id_stomate, 'FLUX100_LCC', itime, &
1943	flux100(:,:,iwplcc), npts*100, horip100_index)
1944	CALL histwrite_p (hist_id_stomate, 'FLUX100_HAR', itime, &
1945	flux100(:,:,iwphar), npts*100, horip100_index)
1946	CALL histwrite_p (hist_id_stomate, 'DefLitSurplus', itime, &
1947	deflitsup_total, npts*nvm, horipft_index)
1948	CALL histwrite_p (hist_id_stomate, 'DefBioSurplus', itime, &
1949	defbiosup_total, npts*nvm, horipft_index)
1950
1951	IF (use_bound_spa) THEN
1952	CALL histwrite_p (hist_id_stomate, 'bound_spa', itime, &
1953	bound_spa, npts*nvm, horipft_index)
1954	ENDIF
1955
1956	IF (do_now_stomate_lcchange) THEN
1957	CALL histwrite_p (hist_id_stomate, 'LCC', itime, &
1958	lcc, npts*nvm, horipft_index)
1959	ENDIF
1960
1961	CALL histwrite_p (hist_id_stomate, 'LAI', itime, &
1962	lai, npts*nvm, horipft_index)
1963	CALL histwrite_p (hist_id_stomate, 'FPC_MAX', itime, &
1964	fpc_max, npts*nvm, horipft_index)
1965	CALL histwrite_p (hist_id_stomate, 'MAXFPC_LASTYEAR', itime, &
1966	maxfpc_lastyear, npts*nvm, horipft_index)
1967	CALL histwrite_p (hist_id_stomate, 'VEGET_MAX', itime, &
1968	veget_max, npts*nvm, horipft_index)
1969	CALL histwrite_p (hist_id_stomate, 'NPP', itime, &
1970	npp_daily, npts*nvm, horipft_index)
1971	CALL histwrite_p (hist_id_stomate, 'GPP', itime, &
1972	gpp_daily, npts*nvm, horipft_index)
1973	CALL histwrite_p (hist_id_stomate, 'IND', itime, &
1974	ind, npts*nvm, horipft_index)
1975	CALL histwrite_p (hist_id_stomate, 'CN_IND', itime, &
1976	cn_ind, npts*nvm, horipft_index)
1977	CALL histwrite_p (hist_id_stomate, 'WOODMASS_IND', itime, &
1978	woodmass_ind, npts*nvm, horipft_index)
1979	CALL histwrite_p (hist_id_stomate, 'TOTAL_M', itime, &
1980	tot_live_biomass(:,:,icarbon), npts*nvm, horipft_index)
1981	CALL histwrite_p (hist_id_stomate, 'LEAF_M', itime, &
1982	biomass(:,:,ileaf,icarbon), npts*nvm, horipft_index)
1983	CALL histwrite_p (hist_id_stomate, 'SAP_M_AB', itime, &
1984	biomass(:,:,isapabove,icarbon), npts*nvm, horipft_index)
1985	CALL histwrite_p (hist_id_stomate, 'SAP_M_BE', itime, &
1986	biomass(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
1987	CALL histwrite_p (hist_id_stomate, 'HEART_M_AB', itime, &
1988	biomass(:,:,iheartabove,icarbon), npts*nvm, horipft_index)
1989	CALL histwrite_p (hist_id_stomate, 'HEART_M_BE', itime, &
1990	biomass(:,:,iheartbelow,icarbon), npts*nvm, horipft_index)
1991	CALL histwrite_p (hist_id_stomate, 'ROOT_M', itime, &
1992	biomass(:,:,iroot,icarbon), npts*nvm, horipft_index)
1993	CALL histwrite_p (hist_id_stomate, 'FRUIT_M', itime, &
1994	biomass(:,:,ifruit,icarbon), npts*nvm, horipft_index)
1995	!!!!! crop variables
1996	CALL histwrite_p (hist_id_stomate, 'CROPYIELD', itime, &
1997	biomass(:,:,ifruit,icarbon), npts*nvm, horipft_index)
1998
1999	CALL histwrite_p (hist_id_stomate, 'BIOMYIELD', itime, &
2000	biomass(:,:,ileaf,icarbon)+biomass(:,:,isapabove,icarbon) &
2001	+biomass(:,:,ifruit,icarbon)+biomass(:,:,icarbres,icarbon), npts*nvm,horipft_index)
2002
2003	CALL histwrite_p (hist_id_stomate, 'CROP_EXPORT', itime, &
2004	crop_export, npts*nvm, horipft_index)
2005	!!!!! end crop variables, xuhui
2006	CALL histwrite_p (hist_id_stomate, 'RESERVE_M', itime, &
2007	biomass(:,:,icarbres,icarbon), npts*nvm, horipft_index)
2008	CALL histwrite_p (hist_id_stomate, 'TOTAL_TURN', itime, &
2009	tot_turnover(:,:,icarbon), npts*nvm, horipft_index)
2010	CALL histwrite_p (hist_id_stomate, 'LEAF_TURN', itime, &
2011	turnover_daily(:,:,ileaf,icarbon), npts*nvm, horipft_index)
2012	CALL histwrite_p (hist_id_stomate, 'SAP_AB_TURN', itime, &
2013	turnover_daily(:,:,isapabove,icarbon), npts*nvm, horipft_index)
2014	CALL histwrite_p (hist_id_stomate, 'ROOT_TURN', itime, &
2015	turnover_daily(:,:,iroot,icarbon), npts*nvm, horipft_index)
2016	CALL histwrite_p (hist_id_stomate, 'FRUIT_TURN', itime, &
2017	turnover_daily(:,:,ifruit,icarbon), npts*nvm, horipft_index)
2018	CALL histwrite_p (hist_id_stomate, 'TOTAL_BM_LITTER', itime, &
2019	tot_bm_to_litter(:,:,icarbon), npts*nvm, horipft_index)
2020	CALL histwrite_p (hist_id_stomate, 'LEAF_BM_LITTER', itime, &
2021	bm_to_litter(:,:,ileaf,icarbon), npts*nvm, horipft_index)
2022	CALL histwrite_p (hist_id_stomate, 'SAP_AB_BM_LITTER', itime, &
2023	bm_to_litter(:,:,isapabove,icarbon), npts*nvm, horipft_index)
2024	CALL histwrite_p (hist_id_stomate, 'SAP_BE_BM_LITTER', itime, &
2025	bm_to_litter(:,:,isapbelow,icarbon), npts*nvm, horipft_index)
2026	CALL histwrite_p (hist_id_stomate, 'HEART_AB_BM_LITTER', itime, &
2027	bm_to_litter(:,:,iheartabove,icarbon), npts*nvm, horipft_index)
2028	CALL histwrite_p (hist_id_stomate, 'HEART_BE_BM_LITTER', itime, &
2029	bm_to_litter(:,:,iheartbelow,icarbon), npts*nvm, horipft_index)
2030	CALL histwrite_p (hist_id_stomate, 'ROOT_BM_LITTER', itime, &
2031	bm_to_litter(:,:,iroot,icarbon), npts*nvm, horipft_index)
2032	CALL histwrite_p (hist_id_stomate, 'FRUIT_BM_LITTER', itime, &
2033	bm_to_litter(:,:,ifruit,icarbon), npts*nvm, horipft_index)
2034	CALL histwrite_p (hist_id_stomate, 'RESERVE_BM_LITTER', itime, &
2035	bm_to_litter(:,:,icarbres,icarbon), npts*nvm, horipft_index)
2036	CALL histwrite_p (hist_id_stomate, 'MAINT_RESP', itime, &
2037	resp_maint, npts*nvm, horipft_index)
2038	CALL histwrite_p (hist_id_stomate, 'GROWTH_RESP', itime, &
2039	resp_growth, npts*nvm, horipft_index)
2040	CALL histwrite_p (hist_id_stomate, 'AGE', itime, &
2041	age, npts*nvm, horipft_index)
2042	CALL histwrite_p (hist_id_stomate, 'HEIGHT', itime, &
2043	height, npts*nvm, horipft_index)
2044	CALL histwrite_p (hist_id_stomate, 'MOISTRESS', itime, &
2045	moiavail_week, npts*nvm, horipft_index)
2046	CALL histwrite_p (hist_id_stomate, 'VCMAX', itime, &
2047	vcmax, npts*nvm, horipft_index)
2048	CALL histwrite_p (hist_id_stomate, 'TURNOVER_TIME', itime, &
2049	turnover_time, npts*nvm, horipft_index)
2050	!!DZADD
2051	! CALL histwrite_p (hist_id_stomate, 'LEAF_FRAC1', itime, leaf_frac(:,:,1), npts*nvm, horipft_index)
2052	! CALL histwrite_p (hist_id_stomate, 'LEAF_FRAC2', itime, leaf_frac(:,:,2), npts*nvm, horipft_index)
2053	! CALL histwrite_p (hist_id_stomate, 'LEAF_FRAC3', itime, leaf_frac(:,:,3), npts*nvm, horipft_index)
2054	! CALL histwrite_p (hist_id_stomate, 'LEAF_FRAC4', itime, leaf_frac(:,:,4), npts*nvm, horipft_index)
2055	!!ENDDZADD
2056
2057	IF ( hist_id_stomate_IPCC > 0 ) THEN
2058	vartmp(:)=SUM(tot_live_biomass(:,:,icarbon)veget_max,dim=2)/1e3contfrac
2059	CALL histwrite_p (hist_id_stomate_IPCC, "cVeg", itime, &
2060	vartmp, npts, hori_index)
2061	vartmp(:)=SUM(tot_litter_carbveget_max,dim=2)/1e3contfrac
2062	CALL histwrite_p (hist_id_stomate_IPCC, "cLitter", itime, &
2063	vartmp, npts, hori_index)
2064	vartmp(:)=SUM(tot_soil_carbveget_max,dim=2)/1e3contfrac
2065	CALL histwrite_p (hist_id_stomate_IPCC, "cSoil", itime, &
2066	vartmp, npts, hori_index)
2067	vartmp(:)=(prod10_total + prod100_total)/1e3
2068	CALL histwrite_p (hist_id_stomate_IPCC, "cProduct", itime, &
2069	vartmp, npts, hori_index)
2070	vartmp(:)=carb_mass_variation/1e3/one_day*contfrac
2071	CALL histwrite_p (hist_id_stomate_IPCC, "cMassVariation", itime, &
2072	vartmp, npts, hori_index)
2073	vartmp(:)=SUM(laiveget_max,dim=2)contfrac
2074	CALL histwrite_p (hist_id_stomate_IPCC, "lai", itime, &
2075	vartmp, npts, hori_index)
2076	vartmp(:)=SUM(gpp_dailyveget_max,dim=2)/1e3/one_daycontfrac
2077	CALL histwrite_p (hist_id_stomate_IPCC, "gpp", itime, &
2078	vartmp, npts, hori_index)
2079	vartmp(:)=SUM((resp_maint+resp_growth)veget_max,dim=2)/1e3/one_daycontfrac
2080	CALL histwrite_p (hist_id_stomate_IPCC, "ra", itime, &
2081	vartmp, npts, hori_index)
2082	vartmp(:)=SUM(npp_dailyveget_max,dim=2)/1e3/one_daycontfrac
2083	CALL histwrite_p (hist_id_stomate_IPCC, "npp", itime, &
2084	vartmp, npts, hori_index)
2085	vartmp(:)=SUM(resp_heteroveget_max,dim=2)/1e3/one_daycontfrac
2086	CALL histwrite_p (hist_id_stomate_IPCC, "rh", itime, &
2087	vartmp, npts, hori_index)
2088	vartmp(:)=SUM(co2_fireveget_max,dim=2)/1e3/one_daycontfrac
2089	CALL histwrite_p (hist_id_stomate_IPCC, "fFire", itime, &
2090	vartmp, npts, hori_index)
2091	vartmp(:)=harvest_above/1e3/one_day*contfrac
2092	CALL histwrite_p (hist_id_stomate_IPCC, "fHarvest", itime, &
2093	vartmp, npts, hori_index)
2094	vartmp(:)=cflux_prod_total/1e3/one_day*contfrac
2095	CALL histwrite_p (hist_id_stomate_IPCC, "fLuc", itime, &
2096	vartmp, npts, hori_index)
2097	!gmjc
2098	IF (enable_grazing) THEN
2099	vartmp(:)=(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire &
2100	-harvest_gm-ranimal_gm-ch4_pft_gm+cinput_gm) & ! specific
2101	& veget_max,dim=2)-cflux_prod_total-harvest_above)/1e3/one_daycontfrac
2102	ELSE
2103	vartmp(:)=(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire) &
2104	& veget_max,dim=2)-cflux_prod_total-harvest_above)/1e3/one_daycontfrac
2105	ENDIF
2106	!end gmjc
2107	! vartmp(:)=(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire) &
2108	! & veget_max,dim=2)-cflux_prod_total-harvest_above)/1e3/one_daycontfrac
2109	CALL histwrite_p (hist_id_stomate_IPCC, "nbp", itime, &
2110	vartmp, npts, hori_index)
2111	vartmp(:)=SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))veget_max,dim=2)/1e3/one_daycontfrac
2112	CALL histwrite_p (hist_id_stomate_IPCC, "fVegLitter", itime, &
2113	vartmp, npts, hori_index)
2114	vartmp(:)=SUM(SUM(soilcarbon_input,dim=2)veget_max,dim=2)/1e3/one_daycontfrac
2115	CALL histwrite_p (hist_id_stomate_IPCC, "fLitterSoil", itime, &
2116	vartmp, npts, hori_index)
2117	vartmp(:)=SUM(biomass(:,:,ileaf,icarbon)veget_max,dim=2)/1e3contfrac
2118	CALL histwrite_p (hist_id_stomate_IPCC, "cLeaf", itime, &
2119	vartmp, npts, hori_index)
2120	vartmp(:)=SUM((biomass(:,:,isapabove,icarbon)+biomass(:,:,iheartabove,icarbon))veget_max,dim=2)/1e3contfrac
2121	CALL histwrite_p (hist_id_stomate_IPCC, "cWood", itime, &
2122	vartmp, npts, hori_index)
2123	vartmp(:)=SUM(( biomass(:,:,iroot,icarbon) + biomass(:,:,isapbelow,icarbon) + biomass(:,:,iheartbelow,icarbon) ) &
2124	& veget_max,dim=2)/1e3contfrac
2125	CALL histwrite_p (hist_id_stomate_IPCC, "cRoot", itime, &
2126	vartmp, npts, hori_index)
2127	vartmp(:)=SUM(( biomass(:,:,icarbres,icarbon) + biomass(:,:,ifruit,icarbon))veget_max,dim=2)/1e3contfrac
2128	CALL histwrite_p (hist_id_stomate_IPCC, "cMisc", itime, &
2129	vartmp, npts, hori_index)
2130	vartmp(:)=SUM((litter(:,istructural,:,iabove,icarbon)+litter(:,imetabolic,:,iabove,icarbon))*&
2131	veget_max,dim=2)/1e3*contfrac
2132	CALL histwrite_p (hist_id_stomate_IPCC, "cLitterAbove", itime, &
2133	vartmp, npts, hori_index)
2134	vartmp(:)=SUM((litter(:,istructural,:,ibelow,icarbon)+litter(:,imetabolic,:,ibelow,icarbon))*&
2135	veget_max,dim=2)/1e3*contfrac
2136	CALL histwrite_p (hist_id_stomate_IPCC, "cLitterBelow", itime, &
2137	vartmp, npts, hori_index)
2138	vartmp(:)=SUM(carbon(:,iactive,:)veget_max,dim=2)/1e3contfrac
2139	CALL histwrite_p (hist_id_stomate_IPCC, "cSoilFast", itime, &
2140	vartmp, npts, hori_index)
2141	vartmp(:)=SUM(carbon(:,islow,:)veget_max,dim=2)/1e3contfrac
2142	CALL histwrite_p (hist_id_stomate_IPCC, "cSoilMedium", itime, &
2143	vartmp, npts, hori_index)
2144	vartmp(:)=SUM(carbon(:,ipassive,:)veget_max,dim=2)/1e3contfrac
2145	CALL histwrite_p (hist_id_stomate_IPCC, "cSoilSlow", itime, &
2146	vartmp, npts, hori_index)
2147	DO j=1,nvm
2148	histvar(:,j)=veget_max(:,j)contfrac(:)100
2149	ENDDO
2150	CALL histwrite_p (hist_id_stomate_IPCC, "landCoverFrac", itime, &
2151	histvar, npts*nvm, horipft_index)
2152	!-
2153	vartmp(:)=zero
2154	DO j = 2,nvm
2155	IF (is_deciduous(j)) THEN
2156	vartmp(:) = vartmp(:) + veget_max(:,j)contfrac100
2157	ENDIF
2158	ENDDO
2159	CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimDec", itime, &
2160	vartmp, npts, hori_index)
2161	!-
2162	vartmp(:)=zero
2163	DO j = 2,nvm
2164	IF (is_evergreen(j)) THEN
2165	vartmp(:) = vartmp(:) + veget_max(:,j)contfrac100
2166	ENDIF
2167	ENDDO
2168	CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimEver", itime, &
2169	vartmp, npts, hori_index)
2170	!-
2171	vartmp(:)=zero
2172	DO j = 2,nvm
2173	IF ( .NOT.(is_c4(j)) ) THEN
2174	vartmp(:) = vartmp(:) + veget_max(:,j)contfrac100
2175	ENDIF
2176	ENDDO
2177	CALL histwrite_p (hist_id_stomate_IPCC, "c3PftFrac", itime, &
2178	vartmp, npts, hori_index)
2179	!-
2180	vartmp(:)=zero
2181	DO j = 2,nvm
2182	IF ( is_c4(j) ) THEN
2183	vartmp(:) = vartmp(:) + veget_max(:,j)contfrac100
2184	ENDIF
2185	ENDDO
2186	CALL histwrite_p (hist_id_stomate_IPCC, "c4PftFrac", itime, &
2187	vartmp, npts, hori_index)
2188	!-
2189	vartmp(:)=SUM(resp_growthveget_max,dim=2)/1e3/one_daycontfrac
2190	CALL histwrite_p (hist_id_stomate_IPCC, "rGrowth", itime, &
2191	vartmp, npts, hori_index)
2192	vartmp(:)=SUM(resp_maintveget_max,dim=2)/1e3/one_daycontfrac
2193	CALL histwrite_p (hist_id_stomate_IPCC, "rMaint", itime, &
2194	vartmp, npts, hori_index)
2195	vartmp(:)=SUM(bm_alloc(:,:,ileaf,icarbon)veget_max,dim=2)/1e3/one_daycontfrac
2196	CALL histwrite_p (hist_id_stomate_IPCC, "nppLeaf", itime, &
2197	vartmp, npts, hori_index)
2198	vartmp(:)=SUM(bm_alloc(:,:,isapabove,icarbon)veget_max,dim=2)/1e3/one_daycontfrac
2199	CALL histwrite_p (hist_id_stomate_IPCC, "nppWood", itime, &
2200	vartmp, npts, hori_index)
2201	vartmp(:)=SUM(( bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iroot,icarbon) )veget_max,dim=2)/1e3/one_daycontfrac
2202	CALL histwrite_p (hist_id_stomate_IPCC, "nppRoot", itime, &
2203	vartmp, npts, hori_index)
2204
2205	CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_X', itime, &
2206	resolution(:,1), npts, hori_index)
2207	CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_Y', itime, &
2208	resolution(:,2), npts, hori_index)
2209	CALL histwrite_p (hist_id_stomate_IPCC, 'CONTFRAC', itime, &
2210	contfrac(:), npts, hori_index)
2211
2212	ENDIF
2213
2214	IF (printlev>=4) WRITE(numout,*) 'Leaving stomate_lpj'
2215
2216	END SUBROUTINE StomateLpj
2217
2218
2219	!! ================================================================================================================================
2220	!! SUBROUTINE : harvest
2221	!!
2222	!>\BRIEF Harvest of croplands
2223	!!
2224	!! DESCRIPTION : To take into account biomass harvest from crop (mainly to take
2225	!! into account for the reduced litter input and then decreased soil carbon. it is a
2226	!! constant (40\%) fraction of above ground biomass.
2227	!!
2228	!! RECENT CHANGE(S) : None
2229	!!
2230	!! MAIN OUTPUT VARIABLE(S): ::harvest_above the harvested biomass
2231	!!
2232	!! REFERENCE(S) :
2233	!! - Piao, S., P. Ciais, P. Friedlingstein, N. de Noblet-Ducoudre, P. Cadule, N. Viovy, and T. Wang. 2009.
2234	!! Spatiotemporal patterns of terrestrial carbon cycle during the 20th century. Global Biogeochemical
2235	!! Cycles 23:doi:10.1029/2008GB003339.
2236	!!
2237	!! FLOWCHART : None
2238	!! \n
2239	!_ ================================================================================================================================
2240
2241	SUBROUTINE harvest(npts, dt_days, veget_max, &
2242	bm_to_litter, turnover_daily, &
2243	harvest_above)
2244
2245	!! 0. Variable and parameter declaration
2246
2247	!! 0.1 Input variables
2248
2249	INTEGER, INTENT(in) :: npts !! Domain size (unitless)
2250	REAL(r_std), INTENT(in) :: dt_days !! Time step (days)
2251	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max !! new "maximal" coverage fraction of a PFT (LAI ->
2252	!! infinity) on ground @tex $(m^2 m^{-2})$ @endtex
2253
2254	!! 0.2 Output variables
2255
2256	!! 0.3 Modified variables
2257
2258	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: bm_to_litter !! [DISPENSABLE] conversion of biomass to litter
2259	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
2260	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: turnover_daily !! Turnover rates
2261	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
2262	REAL(r_std), DIMENSION(npts), INTENT(inout) :: harvest_above !! harvest above ground biomass for agriculture
2263	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
2264	!! 0.4 Local variables
2265
2266	INTEGER(i_std) :: i, j, k, l, m !! indices
2267	REAL(r_std) :: above_old !! biomass of previous time step
2268	!! @tex $(gC m^{-2})$ @endtex
2269	!_ ================================================================================================================================
2270
2271	!! 1. Yearly initialisation
2272
2273	above_old = zero
2274	harvest_above = zero
2275
2276	DO i = 1, npts
2277	DO j = 1,nvm
2278	IF ((.NOT. natural(j)) .AND. (.NOT. is_peat(j))) THEN
2279	above_old = turnover_daily(i,j,ileaf,icarbon) + turnover_daily(i,j,isapabove,icarbon) + &
2280	& turnover_daily(i,j,iheartabove,icarbon) + turnover_daily(i,j,ifruit,icarbon) + &
2281	& turnover_daily(i,j,icarbres,icarbon) + turnover_daily(i,j,isapbelow,icarbon) + &
2282	& turnover_daily(i,j,iheartbelow,icarbon) + turnover_daily(i,j,iroot,icarbon)
2283
2284	turnover_daily(i,j,ileaf,icarbon) = turnover_daily(i,j,ileaf,icarbon)*frac_turnover_daily
2285	turnover_daily(i,j,isapabove,icarbon) = turnover_daily(i,j,isapabove,icarbon)*frac_turnover_daily
2286	turnover_daily(i,j,isapbelow,icarbon) = turnover_daily(i,j,isapbelow,icarbon)*frac_turnover_daily
2287	turnover_daily(i,j,iheartabove,icarbon) = turnover_daily(i,j,iheartabove,icarbon)*frac_turnover_daily
2288	turnover_daily(i,j,iheartbelow,icarbon) = turnover_daily(i,j,iheartbelow,icarbon)*frac_turnover_daily
2289	turnover_daily(i,j,iroot,icarbon) = turnover_daily(i,j,iroot,icarbon)*frac_turnover_daily
2290	turnover_daily(i,j,ifruit,icarbon) = turnover_daily(i,j,ifruit,icarbon)*frac_turnover_daily
2291	turnover_daily(i,j,icarbres,icarbon) = turnover_daily(i,j,icarbres,icarbon)*frac_turnover_daily
2292	harvest_above(i) = harvest_above(i) + veget_max(i,j) * above_old *(un - frac_turnover_daily)
2293	ENDIF
2294	ENDDO
2295	ENDDO
2296
2297	!!$ harvest_above = harvest_above
2298	END SUBROUTINE harvest
2299
2300	!! ================================================================================================================================
2301	!! SUBROUTINE : lpj_cover_peat
2302	!!
2303	!>\BRIEF Peatland area fraction is calculated by TOPMODEL
2304	! 1. IF ok_dgvm, then non-peat vegetations' area is calculated by DGVM (bioclimatic conditions)
2305	! 2. IF .NOT. ok_dgvm, then non-peat vegetations' area is read from vegetation maps
2306	! We adjust non-peatland vegetations coverage according to peatland cover.
2307	! a. Peatland initiates/expands ---> all non-peatland vegetations shrink, proportionally
2308	! b. Peatland shrink/disappear ----> all non-peatland vegetations expand, proportionally
2309	! When peatland shrinks, change of area and C density will be saved (*_save). Then when pealtand expands, peatland will take
2310	! carbon from the *_save first, if fsave is not enough, extra carbon will be taken from natural PFTs
2311	!!
2312	!! \n
2313	!_ ================================================================================================================================
2314
2315	SUBROUTINE lpj_cover_peat(npts, lalo,cn_ind, ind, biomass, veget_max_new, veget_max, &
2316	veget_max_old, litter, litter_avail, litter_not_avail, carbon, &
2317	fuel_1hr, fuel_10hr, fuel_100hr, fuel_1000hr, &
2318	turnover_daily, bm_to_litter, &
2319	co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, gpp_daily, &
2320	deepC_a, deepC_s, deepC_p, &
2321	dt_days,age, PFTpresent, senescence, when_growthinit,&
2322	everywhere, leaf_frac, lm_lastyearmax, npp_longterm,&
2323	carbon_save,deepC_a_save,deepC_s_save,deepC_p_save,delta_fsave,liqwt_max_lastyear)
2324
2325	!! 0. Variable and parameter declaration
2326
2327	!! 0.1 Input variables
2328
2329	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
2330	REAL(r_std),DIMENSION(npts,2),INTENT(in) :: lalo
2331	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: cn_ind !! Crown area
2332	!! @tex $(m^2)$ @endtex per individual
2333	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind !! Number of individuals
2334	!! @tex $(m^{-2})$ @endtex
2335	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max_new !! saved fpeat
2336	!! @tex ($gC individual^{-1}$) @endtex
2337	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate
2338	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max_old
2339	REAL(r_std), DIMENSION(npts), INTENT(in) :: liqwt_max_lastyear
2340	!! 0.2 Output variables
2341
2342	!! 0.3 Modified variables
2343
2344	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) :: litter !! Metabolic and structural litter, above and
2345	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: fuel_1hr
2346	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: fuel_10hr
2347	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: fuel_100hr
2348	REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements), INTENT(inout) :: fuel_1000hr
2349	!! below ground @tex $(gC m^{-2})$ @endtex
2350	REAL(r_std), DIMENSION(npts,nlitt,nvm), INTENT(inout):: litter_avail
2351	REAL(r_std), DIMENSION(npts,nlitt,nvm) , INTENT(inout):: litter_not_avail
2352	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: carbon !! Carbon pool: active, slow, or passive
2353	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass !! Biomass @tex $(gC m^{-2})$ @endtex
2354	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max !! "Maximal" coverage fraction of a PFT (LAI->
2355	!! infinity) on ground (unitless)
2356	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: turnover_daily !! Turnover rates (gC m^{-2} day^{-1})
2357	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: bm_to_litter !! Conversion of biomass to litter
2358	!! @tex $(gC m^{-2} day^{-1})$ @endtex
2359	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: co2_to_bm !! biomass up take for establishment
2360	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: co2_fire
2361	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: resp_hetero
2362	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: resp_maint
2363	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: resp_growth
2364	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gpp_daily
2365
2366	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_a !! Permafrost soil carbon (g/m**3) active
2367	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_s !! Permafrost soil carbon (g/m**3) slow
2368	REAL(r_std), DIMENSION(npts,ndeep,nvm), INTENT(inout) :: deepC_p !! Permafrost soil carbon (g/m**3) passive
2369
2370	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age !! mean age (years)
2371	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: senescence !! plant senescent (only for deciduous trees) Set
2372	!! to .FALSE. if PFT is introduced or killed
2373	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: PFTpresent !! Is pft there (unitless)
2374	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or very
2375	!! localized (unitless)
2376	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: when_growthinit !! how many days ago was the beginning of the
2377	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class
2378	!! (unitless)
2379	!! growing season (days)
2380	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
2381	!! @tex ($gC m^{-2}$) @endtex
2382	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: npp_longterm !! "long term" net primary productivity
2383	!! @tex ($gC m^{-2} year^{-1}$) @endtex
2384	REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: carbon_save
2385	REAL(r_std), DIMENSION(npts,ndeep), INTENT(inout) :: deepC_a_save
2386	REAL(r_std), DIMENSION(npts,ndeep), INTENT(inout) :: deepC_s_save
2387	REAL(r_std), DIMENSION(npts,ndeep), INTENT(inout) :: deepC_p_save
2388	REAL(r_std), DIMENSION(npts), INTENT(inout) :: delta_fsave
2389
2390	!! 0.4 Local variables
2391	REAL(r_std), DIMENSION(npts,ncarb,nvm) :: carbon_save_tmp
2392	INTEGER(i_std) :: i,j,k,m,l !! Index (unitless)
2393	REAL(r_std), DIMENSION(npts,nlitt,nlevs,nelements) :: dilu_lit !! Litter dilution @tex $(gC m^{-2})$ @endtex
2394	REAL(r_std), DIMENSION(npts,nparts,nelements) :: dilu_bio !! Biomass dilution
2395	REAL(r_std), DIMENSION(npts,nparts,nelements) :: dilu_turnover_daily
2396	REAL(r_std), DIMENSION(npts,nparts,nelements) :: dilu_bm_to_litter
2397	REAL(r_std), DIMENSION(npts) :: dilu_co2flux_new
2398	REAL(r_std), DIMENSION(npts) :: dilu_gpp_daily
2399	REAL(r_std), DIMENSION(npts) :: dilu_resp_growth
2400	REAL(r_std), DIMENSION(npts) :: dilu_resp_maint
2401	REAL(r_std), DIMENSION(npts) :: dilu_resp_hetero
2402	REAL(r_std), DIMENSION(npts) :: dilu_co2_to_bm
2403	REAL(r_std), DIMENSION(npts) :: dilu_co2_fire
2404	REAL(r_std), DIMENSION(npts,nvm) :: co2flux_new
2405	REAL(r_std), DIMENSION(npts,nvm) :: co2flux_old
2406	REAL(r_std), DIMENSION(npts,ncarb,nvm) :: carbon_old
2407	REAL(r_std), DIMENSION(npts,ncarb,nvm) :: carbon_tmp
2408
2409	REAL(r_std), DIMENSION(nvm) :: delta_veg !! Conversion factors (unitless)
2410	REAL(r_std) :: delta_veg_sum !! Conversion factors (unitless)
2411	REAL(r_std) :: exp_nat_sum
2412	REAL(r_std) :: nat_sum
2413	REAL(r_std), DIMENSION(npts,nlitt,nelements) :: dilu_f1hr !! Litter dilution @tex $(gC m^{-2})$ @endtex
2414	REAL(r_std), DIMENSION(npts,nlitt,nelements) :: dilu_f10hr !! Litter dilution @tex $(gC m^{-2})$ @endtex
2415	REAL(r_std), DIMENSION(npts,nlitt,nelements) :: dilu_f100hr !! Litter dilution @tex $(gC m^{-2})$ @endtex
2416	REAL(r_std), DIMENSION(npts,nlitt,nelements) :: dilu_f1000hr !! Litter dilution @tex $(gC m^{-2})$ @endtex
2417	REAL(r_std),DIMENSION(npts,nvm) :: delta_ind !! change in number of individuals
2418	REAL(r_std), DIMENSION(npts) :: frac_nat
2419	REAL(r_std), DIMENSION(npts) :: sum_veget_natveg !! Conversion factors (unitless)
2420
2421	REAL(r_std) :: sum_veg
2422	REAL(r_std) :: exp_nat_sum2
2423	REAL(r_std) :: sumvpeat_old, sumvpeat ! last an new sum of peatland vegetation
2424	REAL(r_std) :: rapport ! (S-B) / (S-A)
2425	REAL(r_std) :: SUMveg
2426	REAL(r_std), DIMENSION(npts,nvm) :: veget_tmp
2427	REAL(r_std), DIMENSION(npts) :: delta_fpeat
2428	REAL(r_std), DIMENSION(npts) :: diff_fpeat
2429	REAL(r_std), DIMENSION(npts) :: extra_fra
2430	REAL(r_std), DIMENSION(npts) :: sum_nat
2431	REAL(r_std), DIMENSION(npts,ncarb,nvm) :: dilu_carbon
2432	REAL(r_std), DIMENSION(npts,ndeep,nvm) :: dilu_a
2433	REAL(r_std), DIMENSION(npts,ndeep,nvm) :: dilu_s
2434	REAL(r_std), DIMENSION(npts,ndeep,nvm) :: dilu_p
2435	REAL(r_std),DIMENSION(npts,nparts,nelements) :: bm_new
2436	REAL(r_std),DIMENSION(npts,nparts,nelements) :: bm_gain
2437	REAL(r_std), DIMENSION(npts,nvm) :: co2_take
2438	REAL(r_std), DIMENSION(npts,ncarb) :: carbon_obtain
2439	REAL(r_std), DIMENSION(npts) :: biomass_before
2440	REAL(r_std), DIMENSION(npts) :: biomass_after
2441	REAL(r_std), DIMENSION(npts) :: litter_before
2442	REAL(r_std), DIMENSION(npts) :: litter_after
2443	REAL(r_std), DIMENSION(npts) :: soilc_before
2444	REAL(r_std), DIMENSION(npts) :: soilc_after
2445	REAL(r_std), DIMENSION(npts) :: delta_Cpool
2446	REAL(r_std), DIMENSION(npts,ncarb) :: excess_C
2447	REAL(r_std) :: excess_tmp
2448	REAL(r_std), DIMENSION(npts) :: flux_before
2449	REAL(r_std), DIMENSION(npts) :: flux_after
2450	!_
2451	!================================================================================================================================
2452
2453	IF (printlev>=3) WRITE(numout,*) 'Entering lpj_cover_peat'
2454
2455	! WRITE(numout,*) 'chunjing Entering lpj_cover_peat'
2456
2457	biomass_before(:)=zero
2458	biomass_after(:)=zero
2459	litter_before(:)=zero
2460	litter_after(:)=zero
2461	soilc_before(:)=zero
2462	soilc_after(:)=zero
2463	carbon_save_tmp(:,:,:)=zero
2464	excess_C(:,:)=zero
2465	! flux_before(:)= zero
2466	! flux_after(:)=zero
2467
2468	DO i=1, npts
2469	! WRITE (numout,*) 'qcj check veget_old',veget_max_old(i,:)
2470	! WRITE (numout,*) 'qcj check veget_new',veget_max_new(i,:)
2471	DO l=1,ncarb
2472	DO j=1,nvm
2473	soilc_before(i)=soilc_before(i)+carbon(i,l,j)*veget_max_old(i,j)
2474	ENDDO
2475	ENDDO
2476	ENDDO
2477
2478	!! 1.1 Calculate initial values of vegetation cover
2479	frac_nat(:) = un
2480	sum_veget_natveg(:) = zero
2481	veget_max(:,ibare_sechiba) = un
2482	co2flux_new = undef
2483	co2flux_old = undef
2484
2485	carbon_old(:,:,:)=carbon(:,:,:)
2486
2487	bm_new(:,:,:)=zero
2488	co2_take(:,:)=zero
2489	bm_gain(:,:,:)=zero
2490	carbon_obtain(:,:)=zero
2491
2492	IF (ok_dgvm) THEN
2493	!!! 1. Fraction of non-peat vegetations are results of bioclimatic conditions
2494	DO j = 2,nvm ! loop over PFTs
2495	IF ( natural(j) .AND. .NOT. pasture(j) ) THEN
2496
2497	! Summation of individual tree crown area to get total foliar projected coverage
2498	veget_max(:,j) = ind(:,j) * cn_ind(:,j)
2499	sum_veget_natveg(:) = sum_veget_natveg(:) + veget_max(:,j)
2500	ELSE
2501
2502	!fraction occupied by agriculture needs to be substracted for the DGVM
2503	!this is used below to constrain veget for natural vegetation, see below
2504	frac_nat(:) = frac_nat(:) - veget_max(:,j)
2505
2506	ENDIF
2507
2508	ENDDO ! loop over PFTs
2509
2510	DO i = 1, npts ! loop over grid points
2511
2512	! Recalculation of vegetation projected coverage when ::frac_nat was below ::sum_veget_natveg
2513	! It means that non-natural vegetation will recover ::veget_max as natural vegetation
2514	IF (sum_veget_natveg(i) .GT. frac_nat(i) .AND. frac_nat(i) .GT. min_stomate) THEN
2515	DO j = 2,nvm ! loop over PFTs
2516	IF( natural(j) .AND. .NOT. pasture(j)) THEN
2517	veget_max(i,j) = veget_max(i,j) * frac_nat(i) / sum_veget_natveg(i)
2518	ENDIF
2519	ENDDO ! loop over PFTs
2520	ENDIF
2521	ENDDO ! loop over grid points
2522
2523	! Renew veget_max of bare soil as 0 to difference of veget_max (ibare_sechiba)
2524	! to current veget_max
2525	DO j = 2,nvm ! loop over PFTs
2526	veget_max(:,ibare_sechiba) = veget_max(:,ibare_sechiba) - veget_max(:,j)
2527	ENDDO ! loop over PFTs
2528	veget_max(:,ibare_sechiba) = MAX( veget_max(:,ibare_sechiba), zero )
2529
2530	veget_tmp(:,:)=veget_max(:,:)
2531
2532	!!! 2. Fraction of peatland vegetation are calculated by TOPMODEL
2533	!!! First establishment of pealtand vegetation, companied by increase in biomass
2534	DO i = 1, npts ! Loop over # pixels - domain size
2535	delta_veg(:) = veget_max_new(i,:)-veget_tmp(i,:)
2536	DO j=1, nvm ! Loop over # PFTs
2537	IF (is_peat(j)) THEN
2538	!!!change veget_max to veget_max_new
2539	IF (delta_veg(j) .LT. -min_stomate) THEN
2540	!!!peatland contracts, no limitation
2541	IF (veget_max_new(i,j) > min_stomate) THEN
2542	veget_max(i,j)=veget_max_new(i,j)
2543	ELSE
2544	veget_max(i,j)=zero
2545	ENDIF
2546	ELSEIF(delta_veg(j) .GT. min_stomate) THEN
2547	!!!peatland expands
2548	IF (veget_tmp(i,j) .LE. zero) THEN !!!veget_tmp(i,j) .LE. zero
2549	!!!first initiation of peatland
2550	veget_max(i,j)=veget_max_new(i,j)
2551	ELSE
2552	!!peatland expands from non-zero, the soil need to be wet enough to support growing of peat vegetations
2553	IF (liqwt_max_lastyear(i) > 0.6) THEN
2554	veget_max(i,j)=veget_max_new(i,j)
2555	ENDIF
2556	ENDIF
2557	ENDIF
2558	IF ( veget_max(i,j)-veget_tmp(i,j) > min_stomate) THEN
2559	!! Initial setting of peatland vegetation new establishment
2560	IF (veget_tmp(i,j) .LE. zero) THEN
2561	cn_ind(i,j)=un
2562	ind(i,j)= veget_max(i,j) / cn_ind(i,j)
2563	PFTpresent(i,j) = .TRUE.
2564	everywhere(i,j) = un
2565	senescence(i,j) = .FALSE.
2566	age(i,j) = zero
2567	when_growthinit(i,j) = large_value
2568	leaf_frac(i,j,1) = 1.0
2569	npp_longterm(i,j) = npp_longterm_init
2570	lm_lastyearmax(i,j) = bm_sapl(j,ileaf,icarbon) * ind(i,j)
2571	! bm_new(i,:,icarbon) = bm_sapl(j,:,icarbon) * ind(i,j) /veget_max(i,j)
2572	! co2_take(i,j) = SUM(bm_new(i,:,icarbon))/ dt_days
2573	ENDIF
2574	ENDIF
2575	ENDIF
2576	ENDDO
2577	ENDDO
2578
2579	!!! 3. Adjust fractions of Non-peatland vegetations accordingly
2580
2581	DO i = 1, npts ! Loop over # pixels - domain size
2582	sum_veg=SUM(veget_tmp(i,:))
2583	sumvpeat=zero
2584	sumvpeat_old=zero
2585	DO j = 1,nvm
2586	IF (is_peat(j)) THEN
2587	sumvpeat=sumvpeat+veget_max(i,j)
2588	sumvpeat_old=sumvpeat_old+veget_tmp(i,j)
2589	ENDIF
2590	ENDDO
2591
2592	!!!!3.1 Peat vegetation area increase:
2593	IF (sumvpeat .GT. sumvpeat_old) THEN
2594	! All non-pealtand PFTs (natural vegetation +baresoil) decreases, keep the the proportion of natural PFTs
2595	rapport = ( sum_veg - sumvpeat ) / ( sum_veg - sumvpeat_old)
2596	DO j = 1, nvm
2597	IF ( natural(j) ) THEN
2598	veget_max(i,j) = veget_tmp(i,j) * rapport
2599	ENDIF
2600	ENDDO
2601	ELSE
2602	! Peat vegetation decrease: natural vegetation fractions will not change, the decrease of peat is replaced by bare soil.
2603	! The DGVM will re-introduce natural PFT's.
2604	DO j = 1, nvm
2605	IF (j==1) THEN
2606	veget_max(i,j)=veget_tmp(i,j)+sumvpeat_old-sumvpeat
2607	ENDIF
2608	ENDDO
2609	ENDIF
2610	ENDDO
2611	ELSE !(ok_dgvm=False)
2612	!!!non-peat vegetations and crops' area are read from PFT maps, peatland fraction is calculated by TOPMODEL
2613	!!!substract peatland fraction from natural vegetations (area of crops will not change and no crops growing in peatland)
2614	!!! 1. Fraction of peatland vegetation are calculated by TOPMODEL
2615	!!! First establishment of pealtand vegetation, companied by increase in biomass
2616	DO i = 1, npts ! Loop over # pixels - domain size
2617	delta_veg(:) = veget_max_new(i,:)-veget_max_old(i,:)
2618	DO j=1, nvm ! Loop over # PFTs
2619	IF (is_peat(j)) THEN
2620	!!!change veget_max to veget_max_new
2621	IF (delta_veg(j) .LT. -min_stomate) THEN
2622	!!!peatland contracts, no limitation
2623	IF (veget_max_new(i,j) > min_stomate) THEN
2624	veget_max(i,j)=veget_max_new(i,j)
2625	ELSE
2626	veget_max(i,j)=zero
2627	ENDIF
2628	ELSEIF(delta_veg(j) .GT. min_stomate) THEN
2629	!!!peatland expands
2630	IF (veget_max_old(i,j) .LE. zero) THEN
2631	!!!first initiation of peatland
2632	veget_max(i,j)=veget_max_new(i,j)
2633	ELSE
2634	!!peatland expands from non-zero, the soil need to be wet enough to support growing of peat vegetations
2635	IF (liqwt_max_lastyear(i) > 0.6) THEN
2636	veget_max(i,j)=veget_max_new(i,j)
2637	ENDIF
2638	ENDIF
2639	ENDIF
2640	ENDIF
2641	ENDDO
2642	ENDDO
2643
2644	!!! 2. adjust fractions of Non-peatland natural vegetations accordingly
2645	!!!!!!! no crops growing in peatland
2646	DO j = 1,nvm ! loop over PFTs
2647	IF ( natural(j) .AND. .NOT. pasture(j) ) THEN
2648	sum_veget_natveg(:) = sum_veget_natveg(:) + veget_max_new(:,j)
2649	ENDIF
2650	ENDDO
2651	DO i=1, npts
2652	!!!area of crops (PFT12, PFT13) are from PFT maps and will not occupied by peatland
2653	veget_max(i,12)= veget_max_new(i,12)
2654	veget_max(i,13)= veget_max_new(i,13)
2655	IF (veget_max(i,14) .LE. sum_veget_natveg(i)) THEN
2656	IF (sum_veget_natveg(i) .GT. min_stomate) THEN
2657	DO j=1,nvm
2658	IF ( natural(j) .AND. .NOT. pasture(j) ) THEN
2659	veget_max(i,j)=veget_max_new(i,j)-veget_max_new(i,j)/sum_veget_natveg(i)*veget_max(i,14)
2660	ENDIF
2661	ENDDO
2662	ELSE
2663	veget_max(i,14)=zero
2664	veget_max(i,1:13)=veget_max_new(i,1:13)
2665	veget_max(i,15:16)=veget_max_new(i,15:16) !!!=0
2666	ENDIF
2667	ELSE
2668	veget_max(i,14)=sum_veget_natveg(i)
2669	DO j=1,nvm
2670	IF ( natural(j) .AND. .NOT. pasture(j) ) THEN
2671	veget_max(i,j)=zero
2672	ENDIF
2673	ENDDO
2674	veget_max(i,15:16)=veget_max_new(i,15:16) !!!=0
2675	ENDIF
2676	ENDDO
2677
2678	DO i=1, npts
2679	DO j = 1,nvm
2680	IF (is_peat(j)) THEN
2681	IF ( veget_max(i,j)-veget_max_old(i,j) > min_stomate) THEN
2682	IF (veget_max_old(i,j) .LE. zero) THEN
2683	cn_ind(i,j)=un
2684	ind(i,j)= veget_max(i,j) / cn_ind(i,j)
2685	PFTpresent(i,j) = .TRUE.
2686	everywhere(i,j) = un
2687	senescence(i,j) = .FALSE.
2688	age(i,j) = zero
2689	when_growthinit(i,j) = large_value
2690	leaf_frac(i,j,1) = 1.0
2691	npp_longterm(i,j) = npp_longterm_init
2692	lm_lastyearmax(i,j) = bm_sapl(j,ileaf,icarbon) * ind(i,j)
2693	! bm_new(i,:,icarbon) = bm_sapl(j,:,icarbon) * ind(i,j)/veget_max(i,j)
2694	! co2_take(i,j) = SUM(bm_new(i,:,icarbon))/ dt_days
2695	ENDIF
2696	ENDIF
2697	ENDIF
2698	ENDDO
2699	ENDDO
2700	ENDIF
2701
2702	!!! Correct vegetation fraction, normalize fractions of veget_max smaller than min_vegfrac
2703	!!! to avoid numerical error
2704	DO i = 1, npts ! loop over grid points
2705	DO j=1,nvm
2706	IF ( veget_max(i,j) .LT. min_vegfrac ) THEN
2707	veget_max(i,j) = zero
2708	ENDIF
2709	ENDDO
2710	SUMveg =SUM(veget_max(i,:))
2711	veget_max(i,:) = veget_max(i,:)/SUMveg
2712	ENDDO
2713
2714	!!!Adjust C fluxes and C pools
2715
2716	DO i = 1, npts ! loop over grid points
2717	IF ( ABS( SUM(veget_max(i,:)) - SUM(veget_max_new(i,1:13))) > 1000.*min_stomate ) THEN
2718	WRITE(numout,*) 'qcj check lpj_cover_peat,veget',lalo(i,:)
2719	WRITE(numout,*) 'qcj check veget,sum_new',SUM(veget_max(i,:))
2720	ENDIF
2721
2722	! Calculate the change in veget_max between previous time step and current time step
2723	delta_veg(:) = veget_max(i,:)-veget_max_old(i,:)
2724	delta_veg_sum = SUM(delta_veg,MASK=delta_veg.LT.zero)
2725	delta_fpeat(i)=delta_veg(14)
2726
2727	exp_nat_sum=zero
2728	DO j=1,nvm
2729	IF ((.NOT. is_peat(j)) .AND. (delta_veg(j) .GT. min_stomate)) THEN
2730	exp_nat_sum=exp_nat_sum+delta_veg(j)
2731	ENDIF
2732	ENDDO
2733
2734	dilu_lit(i,:,:,:) = zero
2735	dilu_f1hr(i,:,:) = zero
2736	dilu_f10hr(i,:,:) = zero
2737	dilu_f100hr(i,:,:) = zero
2738	dilu_f1000hr(i,:,:) = zero
2739
2740	dilu_turnover_daily(i,:,:)=zero
2741	dilu_bm_to_litter(i,:,:)=zero
2742	dilu_co2flux_new(i)=zero
2743	dilu_gpp_daily(i)=zero
2744	dilu_resp_growth(i)=zero
2745	dilu_resp_maint(i)=zero
2746	dilu_resp_hetero(i)=zero
2747	dilu_co2_to_bm(i)=zero
2748	dilu_co2_fire(i)=zero
2749
2750	dilu_bio(i,:,:) = zero
2751
2752	dilu_carbon(i,:,:)=zero
2753	dilu_a(i,:,:)=zero
2754	dilu_s(i,:,:)=zero
2755	dilu_p(i,:,:)=zero
2756
2757	!!!!!!!!!!!!!!!!!!!!!!!!!!!!Deal with Biomass C
2758	DO j=1, nvm
2759	IF (is_peat(j)) THEN
2760	IF (delta_fpeat(i) .GT. min_stomate) THEN
2761	bm_new(i,:,:)=delta_fpeat(i)* bm_sapl(j,:,:)
2762	ENDIF
2763
2764	IF (veget_max(i,j).GT.min_stomate) THEN
2765	co2_to_bm(i,j) = co2_to_bm(i,j) + (SUM(bm_new(i,:,icarbon))/ (dt_days*veget_max(i,j)))
2766	biomass(i,j,:,:)= (biomass(i,j,:,:)* veget_max_old(i,j)+ bm_new(i,:,:))/ veget_max(i,j)
2767	ENDIF
2768	ELSE
2769	IF(veget_max(i,j).GT.min_stomate) THEN
2770	biomass(i,j,:,:) = biomass(i,j,:,:) * veget_max_old(i,j) / veget_max(i,j)
2771	ENDIF
2772	ENDIF
2773	ENDDO
2774
2775
2776	!!!!!!!!!!!!!!!!!!!!!!!!!!!!Deal with fluxes
2777	DO j=1, nvm
2778	co2flux_old(i,j)=resp_maint(i,j)+resp_growth(i,j)+resp_hetero(i,j)+co2_fire(i,j)-co2_to_bm(i,j)-gpp_daily(i,j)
2779	co2flux_new(i,j)=resp_maint(i,j)+resp_growth(i,j)+resp_hetero(i,j)+co2_fire(i,j)-co2_to_bm(i,j)-gpp_daily(i,j)
2780	ENDDO
2781
2782	!!!if vegetation coverage decreases, compute dilution of litter, biomass, turnover...
2783	DO j=1, nvm ! loop over PFTs
2784	IF ( delta_veg(j) < -min_stomate ) THEN
2785	dilu_lit(i,:,:,:) = dilu_lit(i,:,:,:) + delta_veg(j) * litter(i,:,j,:,:) / delta_veg_sum
2786	dilu_f1hr(i,:,:) = dilu_f1hr(i,:,:) + delta_veg(j) * fuel_1hr(i,j,:,:) / delta_veg_sum
2787	dilu_f10hr(i,:,:) = dilu_f10hr(i,:,:) + delta_veg(j) * fuel_10hr(i,j,:,:) / delta_veg_sum
2788	dilu_f100hr(i,:,:) = dilu_f100hr(i,:,:) + delta_veg(j) * fuel_100hr(i,j,:,:) / delta_veg_sum
2789	dilu_f1000hr(i,:,:) = dilu_f1000hr(i,:,:) + delta_veg(j) * fuel_1000hr(i,j,:,:) / delta_veg_sum
2790	dilu_turnover_daily(i,:,:)=dilu_turnover_daily(i,:,:)+delta_veg(j)*turnover_daily(i,j,:,:)/delta_veg_sum
2791	dilu_bm_to_litter(i,:,:)=dilu_bm_to_litter(i,:,:)+delta_veg(j)*bm_to_litter(i,j,:,:)/delta_veg_sum
2792	dilu_co2flux_new(i)=dilu_co2flux_new(i)+delta_veg(j)*co2flux_old(i,j)/delta_veg_sum
2793	dilu_gpp_daily(i)=dilu_gpp_daily(i)+delta_veg(j)*gpp_daily(i,j)/delta_veg_sum
2794	dilu_resp_growth(i)=dilu_resp_growth(i)+delta_veg(j)*resp_growth(i,j)/delta_veg_sum
2795	dilu_resp_maint(i)=dilu_resp_maint(i)+delta_veg(j)*resp_maint(i,j)/delta_veg_sum
2796	dilu_resp_hetero(i)=dilu_resp_hetero(i)+delta_veg(j)*resp_hetero(i,j)/delta_veg_sum
2797	dilu_co2_to_bm(i)=dilu_co2_to_bm(i)+delta_veg(j)*co2_to_bm(i,j)/delta_veg_sum
2798	dilu_co2_fire(i)=dilu_co2_fire(i)+delta_veg(j)*co2_fire(i,j)/delta_veg_sum
2799	dilu_bio(i,:,:) = dilu_bio(i,:,:)-delta_veg(j)*biomass(i,j,:,:)
2800	ENDIF
2801	ENDDO ! loop over PFTs
2802
2803	!!!PFTs that increases, recalculate the litter, biomass... with taking into accout the change in veget_max
2804	DO j=1, nvm
2805	IF ( delta_veg(j) > min_stomate) THEN
2806	litter(i,:,j,:,:)=(litter(i,:,j,:,:) * veget_max_old(i,j) + dilu_lit(i,:,:,:) * delta_veg(j)) / veget_max(i,j)
2807	fuel_1hr(i,j,:,:)=(fuel_1hr(i,j,:,:) * veget_max_old(i,j) + dilu_f1hr(i,:,:) * delta_veg(j)) / veget_max(i,j)
2808	fuel_10hr(i,j,:,:)=(fuel_10hr(i,j,:,:) * veget_max_old(i,j) + dilu_f10hr(i,:,:) * delta_veg(j)) / veget_max(i,j)
2809	fuel_100hr(i,j,:,:)=(fuel_100hr(i,j,:,:) * veget_max_old(i,j) + dilu_f100hr(i,:,:) * delta_veg(j)) / veget_max(i,j)
2810	fuel_1000hr(i,j,:,:)=(fuel_1000hr(i,j,:,:) * veget_max_old(i,j) + dilu_f1000hr(i,:,:) * delta_veg(j)) / veget_max(i,j)
2811	IF (is_grassland_manag(j) .AND. is_grassland_grazed(j)) THEN
2812	litter_avail(i,:,j) = litter_avail(i,:,j) * veget_max_old(i,j) / veget_max(i,j)
2813	litter_not_avail(i,:,j) = litter(i,:,j,iabove,icarbon) - litter_avail(i,:,j)
2814	ENDIF
2815	turnover_daily(i,j,:,:)=(turnover_daily(i,j,:,:)veget_max_old(i,j)+dilu_turnover_daily(i,:,:)delta_veg(j))/veget_max(i,j)
2816	bm_to_litter(i,j,:,:)=(bm_to_litter(i,j,:,:)veget_max_old(i,j)+dilu_bm_to_litter(i,:,:)delta_veg(j))/veget_max(i,j)
2817	co2flux_new(i,j)=(co2flux_old(i,j)veget_max_old(i,j)+dilu_co2flux_new(i)delta_veg(j))/veget_max(i,j)
2818	gpp_daily(i,j)=(gpp_daily(i,j)veget_max_old(i,j)+dilu_gpp_daily(i)delta_veg(j))/veget_max(i,j)
2819	resp_growth(i,j)=(resp_growth(i,j)veget_max_old(i,j)+dilu_resp_growth(i)delta_veg(j))/veget_max(i,j)
2820	resp_maint(i,j)=(resp_maint(i,j)veget_max_old(i,j)+dilu_resp_maint(i)delta_veg(j))/veget_max(i,j)
2821	resp_hetero(i,j)=(resp_hetero(i,j)veget_max_old(i,j)+dilu_resp_hetero(i)delta_veg(j))/veget_max(i,j)
2822	co2_fire(i,j)=(co2_fire(i,j)veget_max_old(i,j)+dilu_co2_fire(i)delta_veg(j))/veget_max(i,j)
2823	co2_to_bm(i,j)=(co2_to_bm(i,j)veget_max_old(i,j)+dilu_co2_to_bm(i)delta_veg(j))/veget_max(i,j)
2824	ENDIF
2825	ENDDO
2826
2827	!!!!!!!!!!!!!!!!!!!!!!!!!!!!Deal with Soil carbon pools
2828	DO j=1, nvm
2829	!!!peatland shrink, *_save record old_peat area and C (old_peat)
2830	IF ( delta_veg(j) .LE. zero) THEN
2831	IF (is_peat(j)) THEN
2832	! WRITE (numout,*) 'CQcheck ,peatland shrink'
2833	! WRITE (numout,*) 'CQcheck ,old carbon_save',carbon_save(i,2,j),'fsave',delta_fsave(i)
2834	carbon_save_tmp(i,:,j)=carbon_old(i,:,j)* (-delta_veg(j))
2835	carbon_save(i,:,j) = carbon_save(i,:,j) + carbon_save_tmp(i,:,j)
2836	delta_fsave(i)=delta_fsave(i)-delta_veg(j)
2837	IF ( ok_pc ) THEN
2838	deepC_a_save(i,:)= deepC_a_save(i,:) + deepC_a(i,:,j)* (-delta_veg(j))
2839	deepC_s_save(i,:)= deepC_s_save(i,:) + deepC_s(i,:,j)* (-delta_veg(j))
2840	deepC_p_save(i,:)= deepC_p_save(i,:) + deepC_p(i,:,j)* (-delta_veg(j))
2841	ENDIF
2842	! WRITE (numout,*) 'CQcheck ,new carbon_save',carbon_save(i,2,j),'fsave',delta_fsave(i)
2843	ELSE
2844	IF (delta_veg(j) < - min_stomate) THEN
2845	!!!C from shrinked non-peat PFTs
2846	dilu_carbon(i,:,j) = dilu_carbon(i,:,j) - delta_veg(j) * carbon_old(i,:,j)
2847	ENDIF
2848	ENDIF
2849	ENDIF
2850	ENDDO
2851
2852	IF (delta_fpeat(i) .LE. zero) THEN
2853	!!! peatland shrink
2854	DO j=1,nvm
2855	IF ((delta_veg(j) > min_stomate) .AND. (.NOT. is_peat(j)) .AND. (exp_nat_sum > min_stomate)) THEN
2856	!!! expanding non-peat PFTs get C from shrinking non-peat PFTs and C of old_peat
2857	carbon(i,:,j)=(carbon_old(i,:,j) * veget_max_old(i,j)+SUM(dilu_carbon(i,:,:),dim=2)*delta_veg(j)/exp_nat_sum+&
2858	carbon_save_tmp(i,:,14)*delta_veg(j)/exp_nat_sum)/veget_max(i,j)
2859
2860	IF ( ok_pc ) THEN
2861	IF (carbon_old(i,iactive,j).GT. min_stomate) THEN
2862	deepC_a(i,:,j)=deepC_a(i,:,j)*carbon(i,iactive,j)/carbon_old(i,iactive,j)
2863	! WRITE (numout,*), 'CQ check deepC_a0,pft',j, deepC_a(i,:,j)
2864	ENDIF
2865	IF (carbon_old(i,islow,j).GT. min_stomate) THEN
2866	deepC_s(i,:,j)=deepC_s(i,:,j)*carbon(i,islow,j)/carbon_old(i,islow,j)
2867	ENDIF
2868	IF (carbon_old(i,ipassive,j).GT. min_stomate) THEN
2869	deepC_p(i,:,j)=deepC_p(i,:,j)*carbon(i,ipassive,j)/carbon_old(i,ipassive,j)
2870	ENDIF
2871	ENDIF
2872	ENDIF
2873	ENDDO
2874	ENDIF
2875
2876
2877	!!! peatland expand, occupy old_peat first
2878	IF (delta_fpeat(i) .GT. zero) THEN
2879	IF (delta_fpeat(i) .LE. delta_fsave(i)) THEN
2880	!!!if there is old_peat and old_peat is large enough for peatland expansion, peatland obtain C only from old_peat
2881	carbon_obtain(i,:)= MIN(un,delta_fpeat(i)/delta_fsave(i))*carbon_save(i,:,14)
2882	carbon_save(i,:,14)=MAX((un-delta_fpeat(i)/delta_fsave(i)),zero)*carbon_save(i,:,14)
2883	delta_fsave(i)=MAX(delta_fsave(i)-delta_fpeat(i),zero)
2884	IF ( ok_pc ) THEN
2885	deepC_a_save(i,:)=MAX((un-delta_fpeat(i)/delta_fsave(i)),zero)*deepC_a_save(i,:)
2886	deepC_s_save(i,:)=MAX((un-delta_fpeat(i)/delta_fsave(i)),zero)*deepC_s_save(i,:)
2887	deepC_p_save(i,:)=MAX((un-delta_fpeat(i)/delta_fsave(i)),zero)*deepC_p_save(i,:)
2888	ENDIF
2889	ELSE
2890	!!!if there is no old_peat or old_peat is not large enough, get old_peat first
2891	!then get some c from mineral soil
2892	diff_fpeat(i)=delta_fpeat(i)-delta_fsave(i)
2893	nat_sum=zero
2894	DO j=1,nvm
2895	IF (natural(j) .AND. (veget_max_old(i,j) .GT. min_stomate)) THEN
2896	nat_sum=nat_sum+veget_max_old(i,j)
2897	ENDIF
2898	ENDDO
2899	DO j=1,nvm
2900	IF (natural(j) .AND. (veget_max_old(i,j) .GT. min_stomate)) THEN
2901	carbon_obtain(i,:)=carbon_obtain(i,:)+diff_fpeat(i)(veget_max_old(i,j)/nat_sum)carbon_old(i,:,j)
2902	ENDIF
2903	ENDDO
2904	carbon_obtain(i,:)= carbon_obtain(i,:)+carbon_save(i,:,14)
2905	carbon_save(i,:,14)=zero
2906	delta_fsave(i)= zero
2907	IF (ok_pc) THEN
2908	deepC_a_save(i,:)=zero
2909	deepC_s_save(i,:)=zero
2910	deepC_p_save(i,:)=zero
2911	ENDIF
2912	ENDIF
2913
2914	!!!add *_obtain to peatland
2915	DO j=1,nvm
2916	IF (is_peat(j)) THEN
2917	carbon(i,:,j)=(carbon_old(i,:,j)*veget_max_old(i,j)+ carbon_obtain(i,:))/veget_max(i,j)
2918	IF ( ok_pc ) THEN
2919	IF (carbon_old(i,iactive,j).GT. min_stomate) THEN
2920	deepC_a(i,:,j)=deepC_a(i,:,j)*carbon(i,iactive,j)/carbon_old(i,iactive,j)
2921	ENDIF
2922	IF (carbon_old(i,islow,j).GT. min_stomate) THEN
2923	deepC_s(i,:,j)=deepC_s(i,:,j)*carbon(i,islow,j)/carbon_old(i,islow,j)
2924	ENDIF
2925	IF (carbon_old(i,ipassive,j).GT. min_stomate) THEN
2926	deepC_p(i,:,j)=deepC_p(i,:,j)*carbon(i,ipassive,j)/carbon_old(i,ipassive,j)
2927	ENDIF
2928	ENDIF
2929	ENDIF
2930	ENDDO
2931
2932	carbon_tmp(i,:,:) = carbon(i,:,:)
2933	!!!substract *_obtain from non-peat PFTs
2934	DO l=1, ncarb
2935	!! if contracting non-peat PFTs can provide the *_obtain
2936	IF (SUM(dilu_carbon(i,l,:)) .GE. carbon_obtain(i,l)) THEN
2937	!! subtract *_obtain (has been given to peatland) from dilu_carbon,and then give the remaining dilu_carbon to expanding non-peat PFTs
2938	IF (ANY(delta_veg(1:13) .GT. min_stomate)) THEN
2939	DO j=1,nvm
2940	IF ((delta_veg(j) > min_stomate) .AND. (natural(j)) .AND. (exp_nat_sum > min_stomate)) THEN
2941	carbon(i,l,j)=(carbon_old(i,l,j) * veget_max_old(i,j)+ &
2942	(SUM(dilu_carbon(i,l,:))-carbon_obtain(i,l))*delta_veg(j)/exp_nat_sum)/veget_max(i,j)
2943	! WRITE (numout,*) 'CQ check, C substracted from nonpeat1',SUM(dilu_carbon(i,l,:))-carbon_obtain(i,l)
2944	IF ( ok_pc ) THEN
2945	IF ( (l .EQ. iactive) .AND. (carbon_old(i,iactive,j).GT. min_stomate)) THEN
2946	deepC_a(i,:,j)=deepC_a(i,:,j)*carbon(i,iactive,j)/carbon_old(i,iactive,j)
2947	! WRITE (numout,*), 'CQ check deepC_a2,pft',j,deepC_a(i,:,j)
2948	ENDIF
2949	IF ( (l .EQ. islow) .AND. (carbon_old(i,islow,j).GT. min_stomate)) THEN
2950	deepC_s(i,:,j)=deepC_s(i,:,j)*carbon(i,islow,j)/carbon_old(i,islow,j)
2951	ENDIF
2952	IF ( (l .EQ. ipassive) .AND. (carbon_old(i,ipassive,j).GT. min_stomate)) THEN
2953	deepC_p(i,:,j)=deepC_p(i,:,j)*carbon(i,ipassive,j)/carbon_old(i,ipassive,j)
2954	ENDIF
2955	ENDIF
2956	ENDIF
2957	ENDDO
2958	ELSE
2959	!! peatland is the only PFT that is expanding, add all dilu_carbon to peatland
2960	DO j=1,nvm
2961	IF (is_peat(j)) THEN
2962	carbon(i,l,j)= carbon(i,l,j)+(SUM(dilu_carbon(i,l,:))-carbon_obtain(i,l))/veget_max(i,j)
2963	IF ( ok_pc ) THEN
2964	IF ((carbon_tmp(i,iactive,j).GT. min_stomate).AND. (l .EQ. iactive)) THEN
2965	deepC_a(i,:,j)=deepC_a(i,:,j)*carbon(i,iactive,j)/carbon_tmp(i,iactive,j)
2966	ENDIF
2967	IF ((carbon_tmp(i,islow,j).GT. min_stomate).AND. (l .EQ. islow)) THEN
2968	deepC_s(i,:,j)=deepC_s(i,:,j)*carbon(i,islow,j)/carbon_tmp(i,islow,j)
2969	ENDIF
2970	IF ((carbon_tmp(i,ipassive,j).GT. min_stomate) .AND. (l .EQ. ipassive)) THEN
2971	deepC_p(i,:,j)=deepC_p(i,:,j)*carbon(i,ipassive,j)/carbon_tmp(i,ipassive,j)
2972	ENDIF
2973	ENDIF
2974	ENDIF
2975	ENDDO
2976	ENDIF
2977	ELSE
2978	!! If contracting non-peat PFTs can't provide enough C to *_obtain, all dilu_carbon will be added to peat
2979	!! expanding non-peat PFTs can't receive any C, adjust their C density
2980	DO j=1,nvm
2981	IF ((delta_veg(j) > min_stomate) .AND. (natural(j))) THEN
2982	carbon(i,l,j)=(carbon_old(i,l,j) * veget_max_old(i,j))/veget_max(i,j)
2983	ENDIF
2984	ENDDO
2985	!! Then the difference between dilu_carbon and *_obtain need to be substracted from non-peat PFTs
2986	excess_C(i,l)=SUM(dilu_carbon(i,l,:))-carbon_obtain(i,l)
2987	nat_sum=zero
2988	DO j=1,nvm
2989	IF (natural(j) .AND. (veget_max(i,j) .GT. min_stomate)) THEN
2990	nat_sum=nat_sum+veget_max(i,j)*carbon(i,l,j)
2991	ENDIF
2992	ENDDO
2993	DO j=1,nvm
2994	IF (natural(j) .AND. (veget_max(i,j) .GT. min_stomate) .AND. (excess_C(i,l) .LT. zero) .AND. (nat_sum .GT. zero)) THEN
2995	excess_tmp=MAX(-carbon(i,l,j) veget_max(i,j),excess_C(i,l)veget_max(i,j)*carbon(i,l,j)/nat_sum)
2996	nat_sum=nat_sum-veget_max(i,j)*carbon(i,l,j)
2997	carbon(i,l,j)=(carbon(i,l,j)*veget_max(i,j)+ excess_tmp)/veget_max(i,j)
2998	excess_C(i,l)=excess_C(i,l)-excess_tmp
2999
3000	IF (ok_pc) THEN
3001	IF ( (l .EQ. iactive) .AND. (carbon_old(i,iactive,j).GT. min_stomate)) THEN
3002	deepC_a(i,:,j)=deepC_a(i,:,j)*carbon(i,iactive,j)/carbon_old(i,iactive,j)
3003	! WRITE (numout,*), 'CQ check deepC_a3,pft',j, deepC_a(i,:,j)
3004	ENDIF
3005	IF ( (l .EQ. islow) .AND. (carbon_old(i,islow,j).GT. min_stomate)) THEN
3006	deepC_s(i,:,j)=deepC_s(i,:,j)*carbon(i,islow,j)/carbon_old(i,islow,j)
3007	ENDIF
3008	IF ( (l .EQ. ipassive) .AND. (carbon_old(i,ipassive,j).GT. min_stomate)) THEN
3009	deepC_p(i,:,j)=deepC_p(i,:,j)*carbon(i,ipassive,j)/carbon_old(i,ipassive,j)
3010	ENDIF
3011	ENDIF
3012	ENDIF
3013	ENDDO
3014	IF (excess_C(i,l) .LT. zero) THEN !!!!Sum of C from all mineral soils is smaller than old_peat C
3015	DO j=1,nvm
3016	IF (is_peat(j)) THEN
3017	carbon(i,l,j)=(carbon(i,l,j)*veget_max(i,j)+ excess_C(i,l))/veget_max(i,j)
3018	IF ( ok_pc ) THEN
3019	IF ((carbon_tmp(i,iactive,j).GT. min_stomate) .AND. (l .EQ. iactive)) THEN
3020	deepC_a(i,:,j)=deepC_a(i,:,j)*carbon(i,iactive,j)/carbon_tmp(i,iactive,j)
3021	ENDIF
3022	IF ((carbon_tmp(i,islow,j).GT. min_stomate) .AND. (l .EQ. islow)) THEN
3023	deepC_s(i,:,j)=deepC_s(i,:,j)*carbon(i,islow,j)/carbon_tmp(i,islow,j)
3024	ENDIF
3025	IF ((carbon_tmp(i,ipassive,j).GT. min_stomate) .AND. (l .EQ. ipassive)) THEN
3026	deepC_p(i,:,j)=deepC_p(i,:,j)*carbon(i,ipassive,j)/carbon_tmp(i,ipassive,j)
3027	ENDIF
3028	ENDIF
3029	ENDIF
3030	ENDDO
3031	ENDIF
3032	ENDIF
3033	ENDDO !(ncarb)
3034	ENDIF !(delta_fpeat(i) > min_stomate)
3035	ENDDO !(npts)
3036
3037
3038	DO i=1, npts
3039	DO l=1,ncarb
3040	DO j=1,nvm
3041	soilc_after(i)=soilc_after(i)+carbon(i,l,j)*veget_max(i,j)
3042	ENDDO
3043	ENDDO
3044	ENDDO
3045
3046
3047	! WRITE (numout,*) 'qcj check veget_max last lpj_peat',veget_max(:,:)
3048	DO i=1, npts
3049	! IF (ABS(biomass_after(i)-biomass_before(i)) .GT. 1000.*min_stomate) THEN
3050	! WRITE(numout,*) ' qcj check lpj_cover_peat,biomass',lalo(i,:)
3051	! WRITE(numout,*) ' qcj check biomass,after-before',biomass_after(i),biomass_before(i)
3052	! ENDIF
3053	IF ( ABS(soilc_after(i)-soilc_before(i)) .GT. 1000.*min_stomate) THEN
3054	! WRITE (numout,*) 'qcj check veget_old',veget_max_old(i,:)
3055	! WRITE (numout,*) 'qcj check veget',veget_max(i,:)
3056	WRITE(numout,*) ' qcj check lpj_cover_peat,C',lalo(i,:)
3057	WRITE(numout,*) ' qcj check C,after-before',soilc_after(i),soilc_before(i)
3058	ENDIF
3059	ENDDO
3060
3061	END SUBROUTINE lpj_cover_peat
3062	!_
3063	!================================================================================================================================
3064	END MODULE stomate_lpj

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd_mict_peat_ch4/src_stomate/stomate_lpj.f90 @ 7346

Download in other formats: