Context Navigation

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

Last change on this file since 7346 was 5691, checked in by sebastiaan.luyssaert, 6 years ago
DEV: tested with 13, 37 and 64 PFTs with LCC on different pixels. Some configuration run for 20 years on a given pixel, other crash on another pixel. There is a mass balance problem in sapiens_lcc (ticket #482). This commit fixes a problem with PFT1 in littercalc. This PFT is now fully integrated in LCC and subsequent litter and soil dynamics. veget_max was changed in veget_cov_max where appropriate, a typo in enerbil was corrected.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 157.4 KB

Line
1
2	! ================================================================================================================================
3	!
4	! MODULE : stomate_lpj
5	!
6	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
7	!
8	! LICENCE : IPSL (2006)
9	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
10	!
11	!>\BRIEF Main entry point for daily processes in STOMATE and LPJ (phenology,
12	!! allocation, kill, turn, light, establish, crown, cover, lcchange)
13	!!
14	!!\n DESCRIPTION: None
15	!!
16	!! RECENT CHANGE(S): None
17	!!
18	!! REFERENCE(S) : None
19	!!
20	!! SVN :
21	!! $HeadURL$
22	!! $Date$
23	!! $Revision$
24	!! \n
25	!_ ================================================================================================================================
26
27	MODULE stomate_lpj
28
29	! modules used:
30
31	USE ioipsl_para
32	USE xios_orchidee
33	USE grid
34	USE stomate_data
35	USE constantes
36	USE constantes_soil
37	USE pft_parameters
38	USE structures
39	USE lpj_constraints
40	USE lpj_pftinout
41	USE lpj_kill
42	USE lpj_crown
43	USE lpj_fire
44	USE lpj_gap
45	USE lpj_light
46	USE lpj_establish
47	USE lpj_cover
48	USE sapiens_agriculture
49	USE sapiens_lcchange
50	USE sapiens_kill, ONLY : anthropogenic_mortality
51	USE sapiens_forestry
52	USE sapiens_product_use
53	USE stomate_prescribe
54	USE stomate_phenology
55	USE stomate_mark_kill
56	USE stomate_kill
57	USE stomate_growth_fun_all
58	USE stomate_stand_structure
59	USE stomate_turnover
60	USE stomate_litter
61	USE stomate_laieff, ONLY : find_lai_per_level, calculate_z_level_photo,&
62	combine_lai_levels, fitting_laieff
63	USE stomate_som_dynamics
64	USE stomate_vmax
65	USE stomate_windthrow, ONLY : wind_damage
66
67	USE function_library, ONLY: wood_to_qmdia, cc_to_lai,&
68	wood_to_qmheight, check_biomass_sync, &
69	check_surface_area, check_mass_balance, &
70	wood_to_volume
71
72	IMPLICIT NONE
73
74	! private & public routines
75
76	PRIVATE
77	PUBLIC stomate_lpj_vegetation,stomate_lpj_clear
78
79	CONTAINS
80
81
82	!! ================================================================================================================================
83	!! SUBROUTINE : stomate_lpj_clear
84	!!
85	!>\BRIEF Re-initialisation of variable
86	!!
87	!! DESCRIPTION : This subroutine reinitializes variables. To be used if we want to relaunch
88	!! ORCHIDEE but the routine is not used in current version.
89	!!
90	!! RECENT CHANGE(S) : None
91	!!
92	!! MAIN OUTPUT VARIABLE(S): None
93	!!
94	!! REFERENCE(S) : None
95	!!
96	!! FLOWCHART : None
97	!! \n
98	!_ ================================================================================================================================
99
100	SUBROUTINE stomate_lpj_clear
101
102	CALL prescribe_clear
103	CALL phenology_clear
104	CALL turnover_clear
105	CALL som_dynamics_clear
106	CALL constraints_clear
107	CALL establish_clear
108	CALL fire_clear
109	CALL gap_clear
110	CALL light_clear
111	CALL pftinout_clear
112
113	END SUBROUTINE stomate_lpj_clear
114
115
116	!! ================================================================================================================================
117	!! SUBROUTINE : stomate_lpj
118	!!
119	!>\BRIEF Main entry point for daily processes in STOMATE and LPJ, structures the call sequence
120	!! to the different processes such as dispersion, establishment, competition and mortality of PFT's.
121	!!
122	!! DESCRIPTION : This routine is the main entry point to all processes calculated on a
123	!! daily time step. Is mainly devoted to call the different STOMATE and LPJ routines
124	!! depending of the ok_dgvm (is dynamic veg used) and lpj_constant_mortality (is background mortality used).
125	!! It also prepares the cumulative
126	!! fluxes or pools (e.g TOTAL_M TOTAL_BM_LITTER etc...)
127	!!
128	!! This routine makes frequent use of "weekly", "monthly" and "long term" variables. Quotion is used because
129	!! by default "weekly" denotes 7 days, by default "monthly" denotes 20 days and by default "Long term" denotes
130	!! 3 years. dtslow refers to 24 hours (1 day).
131	!!
132	!!
133	!! RECENT CHANGE(S) : None
134	!!
135	!! MAIN OUTPUT VARIABLE(S): All variables related to stomate and required for LPJ dynamic vegetation mode.
136	!!
137	!! REFERENCE(S) :
138	!! - Krinner, G., N. Viovy, N. de Noblet-DucoudrÃ©, J. OgeÃ©, J. Polcher, P. Friedlingstein, P. Ciais, S. Sitch,
139	!! and I. C. Prentice. 2005. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere
140	!! system. Global Biogeochemical Cycles 19:GB1015, doi:1010.1029/2003GB002199.
141	!! - Sitch, S., B. Smith, I. C. Prentice, A. Arneth, A. Bondeau, W. Cramer, J. O. Kaplan, S. Levis, W. Lucht,
142	!! M. T. Sykes, K. Thonicke, and S. Venevsky. 2003. Evaluation of ecosystem dynamics, plant geography and
143	!! terrestrial carbon cycling in the LPJ dynamic global vegetation model. Global Change Biology 9:161-185.
144	!!
145	!! FLOWCHART : Update with existing flowchart from N Viovy (Jan 19, 2012)
146	!! \n
147	!_ ================================================================================================================================
148
149	SUBROUTINE stomate_lpj_vegetation (npts, dt_days, EndOfYear, &
150	neighbours, resolution, &
151	herbivores, &
152	tsurf_daily, tsoil_daily, t2m_daily, t2m_min_daily, &
153	litterhum_daily, soilhum_daily, &
154	maxmoiavail_lastyear, minmoiavail_lastyear, &
155	gdd0_lastyear, precip_lastyear, &
156	moiavail_month, moiavail_week, &
157	t2m_longterm, t2m_month, t2m_week, &
158	tsoil_month, soilhum_month, &
159	gdd_m5_dormance, gdd_from_growthinit, gdd_midwinter, ncd_dormance, ngd_minus5, &
160	turnover_longterm, gpp_daily, &
161	time_hum_min, hum_min_dormance, maxfpc_lastyear, resp_maint_part, &
162	PFTpresent, age, fireindex, firelitter, &
163	leaf_age, leaf_frac, adapted, regenerate, &
164	plant_status, when_growthinit, litter, &
165	dead_leaves, som, lignin_struc, lignin_wood, &
166	veget_cov_max, veget_cov, veget_cov_max_new, npp_longterm, lm_lastyearmax, &
167	veget_lastlight, everywhere, need_adjacent, RIP_time, &
168	rprof, npp_daily, turnover_daily, turnover_time,&
169	control_moist, control_temp, som_input, &
170	atm_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, &
171	height, deadleaf_cover, vcmax, &
172	nue,bm_to_litter, &
173	prod_s, prod_m, prod_l, flux_s, flux_m, flux_l, &
174	flux_prod_s, flux_prod_m, flux_prod_l, carb_mass_total, &
175	fpc_max, MatrixA, MatrixV, VectorB, VectorU, &
176	Tseason, Tmin_spring_time, onset_date, KF, k_latosa_adapt, &
177	cn_leaf_min_season, nstress_season, moiavail_growingseason, soil_n_min, &
178	rue_longterm, plant_n_uptake_daily, &
179	circ_class_n, circ_class_biomass, forest_managed, forest_managed_lastyear, &
180	tau_eff_leaf, tau_eff_sap, tau_eff_root, &
181	species_change_map, fm_change_map, lpft_replant, &
182	age_stand, rotation_n, last_cut, mai, pai, previous_wood_volume, &
183	mai_count, coppice_dens, &
184	lab_fac, circ_class_dist, store_sum_delta_ba, &
185	harvest_pool_bound, harvest_pool, harvest_type, harvest_cut, harvest_area, &
186	lai_per_level, laieff_fit, h_array_out, z_array_out, max_height_store, &
187	wstress_season, wstress_month, st_dist, litter_demand, &
188	mass_balance_closure, light_tran_to_level_season, p_O2, bact, &
189	CN_som_litter_longterm,nbp_pool_start, &
190	wind_speed_daily, soil_temp_daily, harvest_5y_area)
191
192	!! 0. Variable and parameter declaration
193
194	!! 0.1 input
195
196	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
197	REAL(r_std), INTENT(in) :: dt_days !! Time step of stomate (days)
198	LOGICAL, INTENT(in) :: EndOfYear !! Flag set on the last day of the year used
199	!! to update "yearly variables". This
200	!! variable must be .TRUE. once a year
201	INTEGER(i_std), DIMENSION(npts,8), INTENT(in) :: neighbours !! Indices of the 8 neighbours of each grid
202	!! point [1=N, 2=NE, 3=E, 4=SE,
203	!! 5=S, 6=SW, 7=W, 8=NW]
204	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: resolution !! Resolution at each grid point (m)
205	!! [1=E-W, 2=N-S]
206	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: herbivores !! Time constant of probability of a leaf to
207	!! be eaten by a herbivore (days)
208	REAL(r_std), DIMENSION(npts), INTENT(in) :: tsurf_daily !! Daily surface temperatures (K)
209	REAL(r_std), DIMENSION(npts,nbdl), INTENT(in) :: tsoil_daily !! Daily soil temperatures (K)
210	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_daily !! Daily 2 meter temperatures (K)
211	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_min_daily !! Daily minimum 2 meter temperatures (K)
212	REAL(r_std), DIMENSION(npts), INTENT(in) :: litterhum_daily !! Daily litter humidity (0 to 1, unitless)
213	REAL(r_std), DIMENSION(npts,nbdl), INTENT(in) :: soilhum_daily !! Daily soil humidity (0 to 1, unitless)
214	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: maxmoiavail_lastyear !! Last year's maximum moisture availability
215	!! (0 to 1, unitless)
216	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: minmoiavail_lastyear !! Last year's minimum moisture availability
217	!! (0 to 1, unitless)
218	REAL(r_std), DIMENSION(npts), INTENT(in) :: gdd0_lastyear !! Last year's GDD0 (K)
219	REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_lastyear !! Lastyear's precipitation
220	!! @tex $(mm year^{-1})$ @endtex
221	!! to determine if establishment possible
222	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
223	!! unitless)
224	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
225	!! (0 to 1, unitless)
226	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_longterm !! "Long term" 2 meter reference
227	!! temperatures (K)
228	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_month !! "Monthly" 2-meter temperatures (K)
229	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_week !! "Weekly" 2-meter temperatures (K)
230	REAL(r_std), DIMENSION(npts), INTENT(in) :: Tseason !! "seasonal" 2-meter temperatures (K)
231	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: Tmin_spring_time
232	REAL(r_std), DIMENSION(npts,nvm,2), INTENT(in) :: onset_date
233	REAL(r_std), DIMENSION(npts,nbdl), INTENT(in) :: tsoil_month !! "Monthly" soil temperatures (K)
234	REAL(r_std), DIMENSION(npts,nbdl), INTENT(in) :: soilhum_month !! "Monthly" soil humidity
235	!! (0 to 1, unitless)
236	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
237	!! C (for phenology)
238	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit for crops
239	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
240	!! (for phenology) - this is written to the history files
241	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
242	!! leaves were lost (for phenology)
243	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: turnover_longterm !! "Long term" turnover rate
244	!! @tex $(gC m^{-2} year^{-1})$ @endtex
245	REAL(r_std), DIMENSION(npts), INTENT(inout) :: wind_speed_daily !! Actual wind speed calculated from the actual half hourly
246	REAL(r_std), DIMENSION(npts), INTENT(inout) :: soil_temp_daily !! Soil temperature at 80 cm below ground
247
248	!! @tex $(gC m^{-2} year^{-1})$ @endtex
249	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
250	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
251
252	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: maxfpc_lastyear !! Last year's maximum foliage projected
253	!! coverage for each natural PFT,
254	!! @tex $(m^2 m^{-2})$ @endtex
255	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: resp_maint_part !! Maintenance respiration of different
256	!! plant parts
257	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
258	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: fpc_max !! "Maximal" coverage fraction of a PFT (LAI
259	!! -> infinity) on ground
260	!! @tex $(m^2 m^{-2})$ @endtex
261	!!$ REAL(r_std), DIMENSION(npts,nvm),INTENT(in) :: veget_cov_max_new !! Old timestep "maximal" coverage fraction of a PFT
262	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: cn_leaf_min_season !! Seasonal min CN ratio of leaves
263	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: nstress_season !! N-related seasonal stress (used for allocation)
264	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: moiavail_growingseason !! mean growing season moisture availability
265	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_root !! Effective root turnover time that accounts
266	!! waterstress (days)
267	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_sap !! Effective sapwood turnover time that accounts
268	!! waterstress (days)
269	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_leaf !! Effective leaf turnover time that accounts
270	!! waterstress (days)
271	REAL(r_std), DIMENSION(0:), INTENT(in) :: harvest_pool_bound !! The boundaries of the diameter classes
272	!! in the wood harvest pools
273	!! @tex $(m)$ @endtex
274	REAL(r_std), DIMENSION(:), INTENT(in) :: st_dist !! The probability distribution of trees
275	!! removed from a circ class in case of
276	!! self thinning (unitless).
277	INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: species_change_map !! A map which gives the PFT number that each
278	!! PFT will be replanted as in case of a clearcut.
279	!! (1-nvm,unitless)
280	INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: fm_change_map !! A map which gives the desired FM strategy when
281	!! the PFT will be replanted after a clearcut.
282	!! (1-nvm,unitless)
283	REAL(r_std), DIMENSION(:,:,:),INTENT(in) :: light_tran_to_level_season !! Mean seasonal fraction of light transmitted
284	!! to canopy levels
285
286	!! 0.2 Output variables
287
288	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: npp_daily !! Net primary productivity
289	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
290	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out) :: turnover_daily !! Turnover rates
291	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: co2_fire !! Carbon emitted into the atmosphere by
292	!! fire (living and dead biomass)
293	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
294	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: resp_hetero !! Heterotrophic respiration
295	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
296	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_maint !! Maintenance respiration
297	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
298	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_growth !! Growth respiration
299	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
300
301	REAL(r_std), DIMENSION(npts), INTENT(inout) :: deadleaf_cover !! Fraction of soil covered by dead leaves
302	!! (0 to 1, unitless)
303	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: vcmax !! Maximum rate of carboxylation
304	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out):: bm_to_litter !! Conversion of biomass to litter
305	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
306	REAL(r_std),DIMENSION(npts,nvm), INTENT(out) :: nue !! Nitrogen use Efficiency with impact of leaf
307	!! age (umol CO2 (gN)-1 s-1)
308
309	REAL(r_std), DIMENSION(:,:), INTENT(out) :: lab_fac !! Activity of labile pool factor (??units??)
310	REAL(r_std), DIMENSION(:,:,:), INTENT(out) :: lai_per_level !! This is the LAI per vertical level
311	!! @tex $(m^{2} m^{-2})$
312	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_tot), INTENT(out) &
313	:: h_array_out !! An output of h_array, to use in sechiba
314	REAL(r_std),DIMENSION(npts,nvm,ncirc,nlevels_tot), INTENT(out) &
315	:: z_array_out !! Height above soil of the Pgap points.
316	!! @tex $(m)$ @endtex
317	REAL(r_std), DIMENSION(npts, nvm), INTENT(out) :: max_height_store !! ???
318	TYPE(laieff_type),DIMENSION (:,:,:),INTENT(out) :: laieff_fit !! Fitted parameters for the effective LAI
319
320
321	!! 0.3 Modified variables
322	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
323	!! availability (days)
324	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: hum_min_dormance !! minimum moisture during dormance
325	!! (0-1, unitless)
326	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
327	!! -5 deg C (for phenology)
328	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: height !! Height of vegetation (m)
329	REAL(r_std), DIMENSION(npts,nlevs), INTENT(inout) :: control_moist !! Moisture control of heterotrophic
330	!! respiration (0 to 1, unitless)
331	REAL(r_std), DIMENSION(npts,nlevs), INTENT(inout) :: control_temp !! Temperature control of heterotrophic
332	!! respiration, above and below
333	!! (0 to 1, unitless)
334	REAL(r_std), DIMENSION(npts,ncarb,nvm,nelements), INTENT(inout) :: som_input !! Quantity of carbon going into carbon
335	!! pools from litter decomposition
336	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
337	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: rprof !! Prescribed root depth (m)
338	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: PFTpresent !! Tab indicating which PFTs are present in
339	!! each pixel
340	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age !! Age (years)
341	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: fireindex !! Probability of fire (0 to 1, unitless)
342	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: firelitter !! Longer term litter above the ground that
343	!! can be burned, @tex $(gC m^{-2})$ @endtex
344	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_age !! Leaf age (days)
345	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac !! Fraction of leaves in leaf age class,
346	!! (0 to 1, unitless)
347	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: adapted !! Adaptation of PFT (killed if too cold)
348	!! (0 to 1, unitless)
349	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: regenerate !! "Fitness": Winter sufficiently cold for
350	!! PFT regeneration ? (0 to 1, unitless)
351	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: plant_status !! Growth and phenological status of the plant
352	!! istatus = Phases defined in constantes
353	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
354	!! the growing season (days)
355	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout):: litter !! Metabolic and structural litter, above
356	!! and below ground
357	!! @tex $(gC m^{-2})$ @endtex
358	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: dead_leaves !! Dead leaves on ground, per PFT, metabolic
359	!! and structural,
360	!! @tex $(gC m^{-2})$ @endtex
361	REAL(r_std), DIMENSION(npts,ncarb,nvm,nelements), INTENT(inout) :: som !! Carbon pool: active, slow, or passive,
362	!! @tex $(gC m^{-2})$ @endtex
363	REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout) :: lignin_struc !! Ratio of Lignin/Carbon in structural
364	!! litter, above and below ground,
365	!! @tex $(gC m^{-2})$ @endtex
366	REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout) :: lignin_wood !! Ratio of Lignin/Carbon in woody
367	!! litter, above and below ground,
368	!! @tex $(gC m^{-2})$ @endtex
369	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_cov_max !! "Maximal" coverage fraction of a PFT (LAI
370	!! -> infinity) on ground
371	REAL(r_std),DIMENSION(npts,nvm), INTENT(in) :: veget_cov !! Fractional coverage: actually share of the pixel
372	!! covered by a PFT (fraction of ground area),
373	!! taking into account LAI ??(= grid scale fpc)??
374
375
376	!! 0.3 Modified variables
377
378	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: npp_longterm !! "Long term" mean yearly primary
379	!! productivity
380	!! @tex $(gC m^{-2} year^{-1})$ @endtex
381	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: lm_lastyearmax !! Last year's maximum leaf mass, for each
382	!! PFT @tex $(gC m^{-2})$ @endtex
383	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_lastlight !! Vegetation fractions (on ground) after
384	!! last light competition
385	!! @tex $(m^2 m^{-2})$ @endtex
386	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: everywhere !! Is the PFT everywhere in the grid box or
387	!! very localized (after its introduction)
388	!! (unitless)
389	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: need_adjacent !! In order for this PFT to be introduced,
390	!! does it have to be present in an
391	!! adjacent grid box?
392	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: RIP_time !! How much time ago was the PFT eliminated
393	!! for the last time (y)
394	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout) :: turnover_time !! Turnover_time of leaves for grasses
395	!! (days)
396	REAL(r_std), DIMENSION(:,:),INTENT(inout) :: veget_cov_max_new !! New "maximal" coverage fraction of a PFT
397	!! (LAI -> infinity) Includes
398	!! the nobio fraction, so may sum to
399	!! less than unity. (unitless)
400	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: flux_prod_s !! C-released during first years (short term)
401	!! following land cover change
402	!! @tex ($gC year^{-1}$) @endtex
403	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: flux_prod_m !! Total annual release from decomposition of
404	!! the medium-lived product pool
405	!! @tex ($gC year^{-1}$) @endtex
406	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: flux_prod_l !! Total annual release from decomposition of
407	!! the long-lived product pool
408	!! @tex ($gC year^{-1}$) @endtex
409	REAL(r_std), DIMENSION(:,0:,:), INTENT(inout) :: prod_s !! Short-lived product pool after the annual
410	!! release of each compartment (short + 1 :
411	!! input from year of land cover change)
412	!! @tex ($gC$) @endtex
413	REAL(r_std), DIMENSION(:,0:,:), INTENT(inout) :: prod_m !! Medium-lived product pool after the annual
414	!! release of each compartment (medium + 1 :
415	!! input from year of land cover change)
416	!! @tex ($gC$) @endtex
417	REAL(r_std), DIMENSION(:,0:,:), INTENT(inout) :: prod_l !! long-lived product pool after the annual
418	!! release of each compartment (long + 1 :
419	!! input from year of land cover change)
420	!! @tex ($gC$) @endtex
421	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: flux_s !! Annual release from the short-lived product
422	!! pool @tex ($gC) @endtex
423	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: flux_m !! Annual release from the medium-lived product
424	!! pool compartments @tex ($gC) @endtex
425	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: flux_l !! Annual release from the long-lived product
426	!! pool @tex ($gC) @endtex
427
428
429	REAL(r_std), DIMENSION(npts), INTENT(inout) :: carb_mass_total !! Carbon Mass total (soil, litter, veg)
430	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: KF !! Scaling factor to convert sapwood mass
431	!! into leaf mass (m)
432	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: k_latosa_adapt !! Leaf to sapwood area adapted for long
433	!! term water stress (m)
434	REAL(r_std), DIMENSION(npts,nvm,nionspec), INTENT(inout) :: plant_n_uptake_daily !! Uptake of soil N by plants
435	!! (gN/m**2/day)
436	REAL(r_std), DIMENSION(npts,nvm,nnspec), INTENT(inout) :: soil_n_min !! mineral nitrogen in the soil (gN/m**2)
437	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: rue_longterm !! Longterm radiation use efficiency
438	!! (??units??)
439	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: circ_class_biomass !! Biomass of the componets of the model
440	!! tree within a circumference
441	!! class @tex $(gC ind^{-1})$ @endtex
442	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: circ_class_n !! Number of individuals in each circ class
443	!! @tex $(m^{-2})$ @endtex
444	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: store_sum_delta_ba !! Sum of increase of ba (m^2) which doesn't pass
445	!! through the gap clean
446	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: MatrixA !! Matrix containing the fluxes
447	!! between the carbon pools
448	!! per sechiba time step
449	!! @tex $(gC.m^2.day^{-1})$ @endtex
450
451	REAL(r_std), DIMENSION(npts,nvm,nbpools), &
452	INTENT(inout) :: VectorB !! Vector containing the litter increase per
453	!! sechiba time step
454	!! @tex $(gC m^{-2})$ @endtex
455	REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), &
456	INTENT(inout) :: MatrixV !! Matrix containing the accumulated values of matrixA
457	REAL(r_std), DIMENSION(npts,nvm,nbpools), &
458	INTENT(inout) :: VectorU !! Matrix containing the accumulated values of VectorB
459
460	INTEGER(i_std), DIMENSION (:,:), INTENT(inout) :: forest_managed !! forest management flag (is the forest
461	!! being managed?)
462	INTEGER(i_std), DIMENSION (:,:), INTENT(inout) :: forest_managed_lastyear !! forest management flag for the
463	!! previous year
464	INTEGER(i_std), DIMENSION(:,:), INTENT(inout) :: age_stand !! Age of stand (years)
465	INTEGER(i_std), DIMENSION(:,:), INTENT(inout) :: rotation_n !! Rotation number (number of rotation
466	!! since pft is managed)
467	INTEGER(i_std), DIMENSION(:,:), INTENT(inout) :: last_cut !! Years since last thinning (years)
468
469	LOGICAL, DIMENSION(:,:), INTENT(inout) :: lpft_replant !! Set to true if a PFT has been clearcut
470	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: atm_to_bm !! C and N taken from atmosphere when
471	!! introducing a new PFT (introduced for
472	!! carbon balance closure)
473	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
474	REAL(r_std), DIMENSION(:), INTENT(inout) :: circ_class_dist !! The probability distribution of trees
475	!! in a circ class in case of a
476	!! redistribution (unitless).
477	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: harvest_pool !! The wood and biomass that have been
478	!! havested by humans @tex $(gC)$ @endtex
479	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_type !! Type of management that resulted
480	!! in the harvest (unitless)
481	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_cut !! Type of cutting that was used for the harvest
482	!! (unitless)
483	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_area !! Harvested area (m^{2})
484	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: harvest_5y_area !! Harvested area accumulated in previous five year (m^{2})
485	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: mai !! The mean annual increment
486	!! @tex $(m3 / m2 / year)$ @endtex
487	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: pai !! The period annual increment
488	!! @tex $(m3 / m2 / year)$ @endtex
489	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: previous_wood_volume !! The volume of the tree trunks
490	!! in a stand for the previous year.
491	!! @tex $(m3 / m2 )$ @endtex
492	INTEGER(i_std), DIMENSION(:,:),INTENT(inout) :: mai_count !! The number of times we've
493	!! calculated the volume increment
494	!! for a stand
495	REAL(r_std), DIMENSION(:,:),INTENT(inout) :: coppice_dens !! The density of a coppice at the first
496	!! cutting @tex $( 1 / m**2 )$ @endtex
497	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: wstress_month !! Water stress factor, based on hum_rel_daily
498	!! (unitless, 0-1)
499	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: wstress_season !! Water stress factor, based on hum_rel_daily
500	!! (unitless, 0-1)
501	REAL(r_std), DIMENSION(:,:),INTENT(inout) :: mass_balance_closure !! Missing mass in the carbon balance
502	!! @tex $(gC(N) pixel^{-1} dt^{-1})$@endtex
503	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: p_O2 !! partial pressure of oxigen in the soil (hPa)
504	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: bact !! denitrifier biomass (gC/m**2)
505	REAL(r_std), DIMENSION(:,:,:),INTENT(inout) :: CN_som_litter_longterm !! Longterm CN ratio of litter and som pools (gC/gN)
506	REAL(r_std), DIMENSION(:,:),INTENT(inout) :: nbp_pool_start !! Variable used to calculate the pool based NBP
507	!! @tex $(gC day^{-1} m^{-2})$@endtex
508
509	!! 0.4 Local variables
510	REAL(r_std), DIMENSION(npts,nvm) :: var_real !! temporary variable used to convert an integer
511	!! when writing to a history file
512	REAL(r_std), DIMENSION(npts,nelements) :: prod_s_total !! Total products remaining in the short-lived
513	!! pool after the annual decomposition
514	!! @tex $(gC)$ @endtex
515	REAL(r_std), DIMENSION(npts,nelements) :: prod_m_total !! Total products remaining in the medium-lived
516	!! pool after the annual decomposition
517	!! @tex $(gC)$ @endtex
518	REAL(r_std), DIMENSION(npts,nelements) :: prod_l_total !! Total products remaining in the long-lived
519	!! pool after the annual release
520	!! @tex $(gC)$ @endtex
521	REAL(r_std), DIMENSION(npts,nelements) :: flux_prod_total !! Total flux from decomposition of the short,
522	!! medium and long lived product pools
523	!! @tex $(gC year^{-1})$ @endtex
524
525
526
527	REAL(r_std), DIMENSION(npts,nvm,nelements) :: tot_bm_to_litter !! Total conversion of biomass to litter
528	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
529	REAL(r_std), DIMENSION(npts,nvm,nelements) :: tot_live_biomass !! Total living biomass
530	!! @tex $(gC m{-2})$ @endtex
531	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements):: bm_alloc !! Biomass increase, i.e. NPP per plant part
532	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
533	REAL(r_std), DIMENSION(npts,nvm,nelements) :: tot_turnover !! Total turnover rate
534	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
535	REAL(r_std), DIMENSION(npts,nvm) :: tot_litter_soil_carb !! Total soil and litter carbon
536	!! @tex $(gC m^{-2})$ @endtex
537	REAL(r_std), DIMENSION(npts,nvm) :: tot_litter_carb !! Total litter carbon
538	!! @tex $(gC m^{-2})$ @endtex
539	REAL(r_std), DIMENSION(npts,nvm) :: tot_soil_carb !! Total soil carbon
540	!! @tex $(gC m^{-2})$ @endtex
541	REAL(r_std), DIMENSION(npts) :: carb_mass_variation !! Carbon Mass variation
542	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
543	REAL(r_std), DIMENSION(npts,nvm) :: cn_ind !! Crown area of individuals
544	!! @tex $(m^{2})$ @endtex
545	REAL(r_std), DIMENSION(npts,nvm) :: lai !! Leaf area index OF AN INDIVIDUAL PLANT,
546	!! where a PFT contains n indentical plants
547	!! i.e., using the mean individual approach
548	!! @tex $(m^2 m^{-2})$ @endtex
549	REAL(r_std), DIMENSION(npts,nvm,ncirc,nfm_types,ncut_times) &
550	:: circ_class_kill !! Number of trees within a circ class that needs
551	!! to be killed @tex $(ind m^{-2})$ @endtex
552	!! IMPORTANT: See the note in constantes.f90
553	!! regarding the indicies for this variable.
554	REAL(r_std), DIMENSION(npts,nvm) :: woodmass_ind !! Woodmass of individuals (gC)
555	REAL(r_std), DIMENSION(npts,nvm,nparts) :: f_alloc !! Fraction that goes into plant part
556	!! (0 to 1, unitless)
557	REAL(r_std), DIMENSION(npts) :: avail_tree !! Space availability for trees
558	!! (0 to 1, unitless)
559	REAL(r_std), DIMENSION(npts) :: avail_grass !! Space availability for grasses
560	!! (0 to 1, unitless)
561	INTEGER :: j,k,l,ipts,ivm,icut,ifm,ispec
562	INTEGER :: ilev, ilit,icarb,iele, ipar, imbc, icir, iyear
563	!!$ REAL(r_std),DIMENSION(npts,nvm) :: veget_cov_max_tmp !! "Maximal" coverage fraction of a PFT
564	!!$ !! (LAI-> infinity) on ground (unitless)
565	REAL(r_std), DIMENSION(npts,nvm) :: mortality !! Fraction of individual dying this time
566	!! step (0 to 1, unitless)
567	REAL(r_std), DIMENSION(npts) :: vartmp !! Temporary variable used to add history
568	REAL(r_std), DIMENSION(npts,nvm) :: histvar !! History variables
569	REAL(r_std), DIMENSION(npts,nvm) :: use_reserve !! Mass taken from carbohydrate reserve
570	!! @tex $(gC m^{-2})$ @endtex
571	CHARACTER(LEN=2), DIMENSION(nelements) :: element_str !! string suffix indicating element
572	REAL(r_std), DIMENSION(npts,nvm) :: vcmax_new
573	REAL(r_std), DIMENSION(npts,nvm) :: gammas !! Slope for individual tree growth (m)
574	REAL(r_std), DIMENSION(npts,nvm) :: sigma !! Threshold for indivudal tree growth (m)
575	REAL(r_std), DIMENSION(npts,nvm) :: lrake_frac !! Relative amount of litter that raked (-)
576	REAL(r_std), DIMENSION(npts,nvm) :: qm_height !! Quadratic mean height of vegetation (m)
577	REAL(r_std), DIMENSION(npts,nvm) :: qm_dia !! Quadratic mean diameter of vegetation (m)
578	REAL(r_std), DIMENSION(npts,nvm) :: rdi !! relative density index (ntrees/Nmax)
579	REAL(r_std), DIMENSION(npts,nvm) :: rdi_target_upper !! relative density index (ntrees/Nmax)
580	REAL(r_std), DIMENSION(npts,nvm) :: rdi_target_lower !! relative density index (ntrees/Nmax)
581	REAL(r_std), DIMENSION(npts,nvm,nlevels) :: lai_per_level_temp !! This is the LAI per vertical level
582	!! @tex $(m^{2} m^{-2})$
583	REAL(r_std), DIMENSION(npts,nvm,nlevels_tot) :: z_level_photo !! The height of the levels that we will
584	!! use to calculate the effective LAI for
585	!! the albedo routines and photosynthesis.
586	!! @tex $(m)$ @endtex
587	REAL(r_std), DIMENSION(npts) :: litter_demand !! The amount of litter which will
588	!! be moved from forest to crop pools
589	!! at the end of the year.
590	!! @tex $( gC year^{-1} )$ @endtex
591	CHARACTER(30) :: var_name !! Temporary variable name for histwrite
592	REAL(r_std), DIMENSION(npts) :: pixel_area !! surface area of the pixel @tex ($m^{2}$) @endtex
593	REAL(r_std), DIMENSION(npts,nvm) :: N_support !! Nitrogen which is added to the ecosystem to
594	!! support vegetation growth (only used when
595	!! IMPOSE_CN .EQ. true) (gC/m2/day)
596	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
597	!! mass balance chech for this routine
598	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
599	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
600	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
601	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start pool of this routine
602	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
603	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! End pool of this routine
604	REAL(r_std), DIMENSION(npts,nvm) :: temp
605	REAL(r_std), DIMENSION(npts,nvm,nlevels_tot) :: dummy !! dummy variable for prescribe
606
607	! Debug
608	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements) :: temp_bio !! Biomass @tex $(gC m^{-2})$ @endtex
609	CHARACTER(30) :: losses !! Loss distribution across age classes
610	REAL(r_std), DIMENSION(npts,nvm,nupdates,nelements) :: atm_to_bm_hist !! History of atm_to_bm for C and N for the
611	!! different subroutines where it is used recruitment
612	REAL(r_std), DIMENSION(npts,nvm,nupdates) :: veget_cov_max_hist !! History of veget_cov_max for C and N for the
613	integer, dimension(3) :: pft_t
614
615
616	! Windthrow
617	REAL(r_std), DIMENSION(npts,nvm,wind_years) :: tmp_5y_area !! Temporary array for harvested area in previous 5years (m^{2})
618	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements,nfm_types,ncut_times) &
619	:: biomass_cut !! Biomass change due to ncut
620	REAL(r_std), DIMENSION(npts,ncut_times) :: wood_volume_pix_cut !! Wood volume change due to ncut
621
622	!_ ================================================================================================================================
623
624	IF (printlev>=2) WRITE(numout,*) 'Entering stomate_lpj'
625
626	!!$ !! 0. Calculate the LAI from the biomass.
627	!!$ !! This is done to prevent passing an incorrect
628	!!$ !! LAI around, since LAI is diagnostic.
629	!!$ !! We can't do this if we use an LAI map! So this breaks DO_STOMATE=N
630	!!$ lai(:,ibare_sechiba) = zero
631	!!$ DO ipts = 1, npts
632	!!$ DO ivm = 1, nvm
633	!!$ lai(ipts,ivm) = cc_to_lai(circ_class_biomass(ipts,ivm,:,ileaf,icarbon),&
634	!!$ circ_class_n(ipts,ivm,:),ivm)
635	!!$ ENDDO
636	!!$ ENDDO
637
638	!! 1. Initializations
639
640	!! 1.1 Initialize local printlev
641	!printlev_loc=get_printlev('stomate_lpj')
642
643	!! 1.2 Initialize variables to zero
644	co2_fire(:,:) = zero
645	npp_daily(:,:) = zero
646	resp_maint(:,:) = zero
647	resp_growth(:,:) = zero
648	bm_to_litter(:,:,:,:) = zero
649	cn_ind(:,:) = zero
650	woodmass_ind(:,:) = zero
651	turnover_daily(:,:,:,:) = zero
652	circ_class_kill(:,:,:,:,:) = zero
653	bm_alloc(:,:,:,:) = zero
654	rdi(:,:) = zero
655	rdi_target_upper(:,:) = zero
656	rdi_target_lower(:,:) = zero
657	lrake_frac(:,:) = zero
658	use_reserve(:,:) = zero
659	atm_to_bm_hist(:,:,:,:) = zero
660	veget_cov_max_hist(:,:,:) = zero
661	harvest_pool(:,:,:,:) = zero
662	biomass_cut(:,:,:,:,:,:) = zero
663
664	!! 1.3 Surface area (m2) of the pixel
665	pixel_area(:) = resolution(:,1) * resolution(:,2)
666
667	!! 1.4 Initialize check for mass balance closure
668	! This test is always performed. If err_act.EQ.1 then
669	! the value of the mass balance error -if any- is
670	! written to the history file.
671
672	! Based on the dimensions of the variables it looks like
673	! we can distinguish different PFTs but that is not true.
674	! We need to account for the wood products and those are
675	! only defined at the pixel level. The highest resolution
676	! we can thus check for mass balance closure is the pixel.
677	check_intern(:,:,:,:) = zero
678	pool_start(:,:,:) = zero
679	DO iele = 1,nelements
680
681	! atm_to_bm has as intent inout, the variable
682	! accumulates carbon over the course of a day.
683	! Use the difference between start and the end of
684	! this routine
685	DO ipar = 1,nparts
686	DO icir=1, ncirc
687	pool_start(:,:,iele) = pool_start(:,:,iele) + &
688	(circ_class_biomass(:,:,icir,ipar,iele) * &
689	circ_class_n(:,:,icir) * veget_cov_max(:,:))
690	ENDDO
691	ENDDO
692
693	! The biomass harvest pool shouldn't be multiplied by veget_cov_max
694	! its units are already in gC or gN pixel-1. The total amount
695	! for the pixel was stored in harvest_pool. Account for C and N
696	! stored in the wood product pools. There are no longer PFTs
697	pool_start(:,1,iele) = pool_start(:,1,iele) + &
698	(SUM(SUM(harvest_pool(:,:,:,iele),3),2) + &
699	SUM(prod_l(:,:,iele),2) + SUM(prod_m(:,:,iele),2) + &
700	SUM(prod_s(:,:,iele),2) ) / pixel_area(:)
701
702	ENDDO
703
704	! Account for the N-uptake calculated in stomate.f90
705	pool_start(:,:,initrogen) = pool_start(:,:,initrogen) + &
706	(plant_n_uptake_daily(:,:,iammonium) + &
707	plant_n_uptake_daily(:,:,initrate)) * veget_cov_max(:,:)
708
709	!! DEBUG !!
710	! Debug...only use if you know what you are doing, not for general use
711	IF (.FALSE.) THEN
712	pft_t=(/ 13, 16, 1/)
713	DO j=1, size(pft_t)
714	ivm=pft_t(j)
715	DO iele = 1,nelements
716	WRITE(numout,*) 'veget_cov_max: ', veget_cov_max(test_grid,ivm), 'ipts:',test_grid, 'ivm:',ivm
717	WRITE(numout,) 'atm_to_bm: ', atm_to_bm_hist(test_grid,ivm,:,iele)veget_cov_max(test_grid,ivm)
718	WRITE(numout,*) 'element ',iele
719	WRITE(numout,) 'biomass',SUM(SUM(circ_class_biomass(test_grid,ivm,:,:,iele),1)circ_class_n(test_grid,ivm,:)) * &
720	veget_cov_max(test_grid,ivm)
721	WRITE(numout,) 'in bm_to_litter:',SUM(bm_to_litter(test_grid,ivm,:,iele))veget_cov_max(test_grid,ivm)
722	WRITE(numout,)'in turnover_daily:',SUM(turnover_daily(test_grid,ivm,:,iele))veget_cov_max(test_grid,ivm)
723	WRITE(numout,) 'in litter:',SUM(SUM(litter(test_grid,:,ivm,:,iele),2))veget_cov_max(test_grid,ivm)
724	WRITE(numout,) 'in som:',SUM(som(test_grid,:,ivm,iele))veget_cov_max(test_grid,ivm)
725	WRITE(numout,*) 'in harvest_pool:',SUM(harvest_pool(test_grid,ivm,:,iele),1)/pixel_area(test_grid)
726	!wrong dimensions WRITE(numout,*) 'in prod_s:',prod_s(test_grid,ivm,iele)/pixel_area(test_grid)
727	!wrong dimensions WRITE(numout,*) 'in prod_m:',prod_m(test_grid,ivm,iele)/pixel_area(test_grid)
728	!wrong dimensions WRITE(numout,*) 'in prod_l:',prod_l(test_grid,ivm,iele)/pixel_area(test_grid)
729	WRITE(numout,) 'in gpp_daily:',gpp_daily(test_grid,ivm)veget_cov_max(test_grid,ivm)
730	WRITE(numout,) 'inresp_maint:',resp_maint(test_grid,ivm)veget_cov_max(test_grid,ivm)
731	WRITE(numout,) 'in resp_growth:',resp_growth(test_grid,ivm)veget_cov_max(test_grid,ivm)
732	!wrong dimensions WRITE(numout,*) 'in flux_s:',flux_s(test_grid,ivm,iele)/pixel_area(test_grid)
733	!wrong dimensions WRITE(numout,*) 'in flux_m:',flux_m(test_grid,ivm,iele)/pixel_area(test_grid)
734	!wrong dimensions WRITE(numout,*) 'in flux_l:',flux_l(test_grid,ivm,iele)/pixel_area(test_grid)
735	WRITE(numout,*) 'in plant_n_uptake_daily:',SUM(plant_n_uptake_daily(test_grid,ivm,:))
736	ENDDO
737	ENDDO
738	ENDIF
739	!!! END DEBUG !!!
740
741	! Initialize check for area conservation
742	veget_cov_max_hist(:,:,ibeg) = veget_cov_max(:,:)
743	atm_to_bm_hist(:,:,ibeg,:) = atm_to_bm(:,:,:)
744
745
746	! +++CHECK+++
747	! Check and document why we do this
748	IF (EndOfYear) THEN
749	forest_managed_lastyear(:,:) = forest_managed(:,:)
750	ENDIF
751	! +++++++++++
752
753	!! 1.6 Update stand age
754	! Here we increment the age of the stand, if it's the end of the
755	! year. Notice that we need to do this before the forestry routines,
756	! since the forest has had a chance to grow since the last time
757	! we called the routine. Doing it this way may cause a small
758	! problem since forestry is only done every year. For example,
759	! let us say that a forest die-off occurs due to natural reasons
760	! at the end of December. The age of that stand will be set to zero.
761	! However, when the end of year happens the age will be incremented
762	! by one here, even though the forest is only a couple weeks old.
763	! There doesn't seem to be a better way to do this, though, so a
764	! forest age of x should be taken to mean that the forest is
765	! between x-1 and x years old. I see this only being an issue
766	! in short rotation coppices, but those forests should never have
767	! die off because they are harvested every couple years.
768	IF (EndOfYear) THEN
769	DO ipts=1,npts
770	DO ivm=1,nvm
771	IF(veget_cov_max(ipts,ivm) .GT. min_stomate)THEN
772	age_stand(ipts,ivm) = age_stand(ipts,ivm) + 1
773
774	! last_cut is similar to stand age, but it measures the
775	! time since any human intervention, either thinning or
776	! clearcutting.
777	last_cut(ipts,ivm) = last_cut(ipts,ivm) + 1
778
779	ELSE
780
781	age_stand(ipts,ivm) = 0
782	last_cut(ipts,ivm) = 0
783
784	ENDIF
785	ENDDO
786	ENDDO
787	ENDIF
788
789	!! 1.7 Calculate some vegetation characteristics
790
791	! +++CHECK++++
792	! Seems useless except for the veget_cov_max issue which I don't
793	! understand. Crown and height are now prognostic and can
794	! be calculated from biomass whenever needed. No need to
795	! call crown here.
796	!! 1.7.1 Calculate some vegetation characteristics
797	! Calculate cn_ind (individual crown mass) and individual height from
798	! state variables if running DGVM or dynamic mortality in static cover mode
799	!?? Explain (maybe in the header once) why you mulitply with veget_cov_max in the DGVM
800	!?? and why you don't multiply with veget_cov_max in stomate.
801	!!$ IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
802	!!$ IF(ok_dgvm) THEN
803	!!$ WHERE (ind(:,:).GT.min_stomate)
804	!!$ woodmass_ind(:,:) = &
805	!!$ ((biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
806	!!$ +biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon)) &
807	!!$ *veget_cov_max(:,:))/ind(:,:)
808	!!$ ENDWHERE
809	!!$ ELSE
810	!!$ WHERE (ind(:,:).GT.min_stomate)
811	!!$ woodmass_ind(:,:) = &
812	!!$ (biomass(:,:,isapabove,icarbon)+biomass(:,:,isapbelow,icarbon) &
813	!!$ +biomass(:,:,iheartabove,icarbon)+biomass(:,:,iheartbelow,icarbon))/ind(:,:)
814	!!$ ENDWHERE
815	!!$ ENDIF
816	!!$
817	!!$ CALL crown (npts, PFTpresent, &
818	!!$ ind, biomass, woodmass_ind, &
819	!!$ veget_cov_max, cn_ind, height)
820	!!$ ENDIF
821
822	!! 2. Prescribe vegetation characteristics if the vegetation is not dynamic
823	! At the first call this routine takes atmospheric CO2 and converts it
824	! into carbon pools to basically avoid modelling seed germination.
825	! When lpft_replant is used we ware not sure whether we want to replant
826	! with the same or another species. Replanting with another species is
827	! a lcchange and therefore when lpf_replant is used it will be dealt with
828	! at the end of the year in sapiens_lcchange. In this case replanting
829	! will only take place at the end of the year. If a PFT dies during the
830	! year it will be left fallow until the last day of the year.
831	! +++CHECK+++
832	! Old issue, check whether the code can deal with it - I [SL] think this
833	! should not be an issue. The model now always has density and individuals.
834	! IF the DGVM is not activated, the density of individuals and their crown
835	! areas don't matter, but they should be defined for the case we switch on
836	! the DGVM afterwards. At the first call, if the DGVM is not activated,
837	! impose a minimum biomass for prescribed PFTs and declare them present.
838	! +++++++++++
839
840	! Debug
841	!printlev_loc=get_printlev('stomate_lpj')
842	IF (printlev_loc>=4) CALL debug_write(npts, 'before prescribe', &
843	circ_class_biomass, circ_class_n, circ_class_kill, &
844	plant_status, veget_cov_max)
845	!-
846
847	dummy(:,:,:) = un
848	CALL prescribe (npts, veget_cov_max, dt_days, PFTpresent, &
849	everywhere, when_growthinit, leaf_frac, circ_class_dist, &
850	circ_class_n, circ_class_biomass, atm_to_bm, &
851	forest_managed, KF, plant_status, age, npp_longterm, &
852	lm_lastyearmax, tau_eff_leaf, tau_eff_sap, tau_eff_root, &
853	k_latosa_adapt, dummy, EndOfYear, lpft_replant, species_change_map)
854
855	!We store intermediate veget_cov_max and atm_to_bm in order to calculate
856	!the mass balance closure at the end odf the subroutine
857	atm_to_bm_hist(:,:,ipre,:) = atm_to_bm(:,:,:)
858	veget_cov_max_hist(:,:,ipre) = veget_cov_max(:,:)
859	! atm_to_bm is always reset to zero begore going to another subroutine
860	atm_to_bm(:,:,:) = zero
861
862
863	!! 2. Climatic constraints for PFT presence and regenerativeness
864
865	! Call this even when DGVM is not activated so that "adapted"
866	! and "regenerate"
867	! are kept up to date for the moment when the DGVM is activated.
868	! +++CHECK+++
869	! Uncomment when working on the DGVM
870	!!$ CALL constraints (npts, dt_days, &
871	!!$ t2m_month, t2m_min_daily,when_growthinit, Tseason, &
872	!!$ adapted, regenerate)
873	! +++++++++++
874
875	!! 3. Determine introduction and elimination of PTS based on climate criteria
876	! +++CHECK+++
877	! Uncomment when working on the DGVM
878	!!$ IF ( ok_dgvm ) THEN
879	!!$
880	!!$ !! 3.1 Calculate introduction and elimination
881	!!$ CALL pftinout (npts, dt_days, adapted, regenerate, &
882	!!$ neighbours, veget_cov_max, &
883	!!$ biomass, ind, cn_ind, age, leaf_frac, npp_longterm, lm_lastyearmax, senescence, &
884	!!$ PFTpresent, everywhere, when_growthinit, need_adjacent, RIP_time, &
885	!!$ atm_to_bm, avail_tree, avail_grass)
886	!!$
887	!!$ !! 3.2 Reset attributes for eliminated PFTs.
888	!!$ ! This also kills PFTs that had 0 leafmass during the last year. The message
889	!!$ ! "... after pftinout" is misleading in this case.
890	!!$ CALL kill (npts, 'pftinout ', lm_lastyearmax, &
891	!!$ ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
892	!!$ lai, age, leaf_age, leaf_frac, npp_longterm, &
893	!!$ when_growthinit, everywhere, veget_cov_max, bm_to_litter)
894	!!$
895	!!$ !! 3.3 Calculate new crown area and diameter
896	!!$ ! Calculate new crown area, diameter and maximum vegetation cover**[No longer used in the subroutine]
897	!!$ ! unsure whether this is really required
898	!!$ ! - in theory this could ONLY be done at the END of stomate_lpj
899	!!$ ! calculate woodmass of individual tree
900	!!$ WHERE ((ind(:,:).GT.min_stomate))
901	!!$ WHERE ( veget_cov_max(:,:) .GT. min_stomate)
902	!!$ woodmass_ind(:,:) = &
903	!!$ ((biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
904	!!$ + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))*veget_cov_max(:,:))/ind(:,:)
905	!!$ ELSEWHERE
906	!!$ woodmass_ind(:,:) =(biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
907	!!$ + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))/ind(:,:)
908	!!$ ENDWHERE
909	!!$
910	!!$ ENDWHERE
911	!!$
912	!!$ ! Calculate crown area and diameter for all PFTs (including the newly established)
913	!!$ CALL crown (npts, PFTpresent, &
914	!!$ ind, biomass, woodmass_ind, &
915	!!$ veget_cov_max, cn_ind, height)
916	!!$
917	!!$ ENDIF
918	! +++++++++++
919
920	!! 4. Phenology
921
922	! Forests and grasses or the so called natural vegetation in ORCHIDEE
923	! undergo a real phenology based on climatology. This means that the
924	! plants are planted the first day of the year and that they live from
925	! their reserves until the day of bud burst. Crops are planted
926	! (and harvested) every year. We plant them the day of bud burst. The
927	! biggest difference is that crops do not need reserves to live
928	! between January 1st and the day of bud burst.
929
930	! Debug
931	!printlev_loc=get_printlev('stomate_lpj')
932	IF (printlev_loc>=4) CALL debug_write(npts, 'before phenology', &
933	circ_class_biomass, circ_class_n, circ_class_kill, &
934	plant_status, veget_cov_max)
935
936	!
937	CALL phenology (npts, dt_days, PFTpresent, &
938	veget_cov_max, &
939	t2m_longterm, t2m_month, t2m_week, gpp_daily, &
940	maxmoiavail_lastyear, minmoiavail_lastyear, &
941	moiavail_month, moiavail_week, &
942	gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
943	plant_status, time_hum_min, &
944	leaf_frac, leaf_age, &
945	when_growthinit, atm_to_bm, circ_class_n, &
946	circ_class_biomass, KF, &
947	tau_eff_leaf, tau_eff_sap, tau_eff_root, age, &
948	everywhere, npp_longterm, lm_lastyearmax, k_latosa_adapt, &
949	cn_leaf_min_season, plant_n_uptake_daily)
950
951	!We store intermediate veget_cov_max and atm_to_bm in order to calculate
952	!the mass balance closure at the end odf the subroutine
953	atm_to_bm_hist(:,:,iphe,:) = atm_to_bm(:,:,:)
954	veget_cov_max_hist(:,:,iphe) = veget_cov_max(:,:)
955	atm_to_bm(:,:,:) = zero
956
957	!! 5. Allocate C to different plant parts
958
959	! Allometry based allocation and intra-stand competition (based on
960	! Sitch et al 2003, Zaehle et al 2010 and Deleuze 2004)
961
962	! Debug
963	!printlev_loc=get_printlev('stomate_lpj')
964	IF (printlev_loc>=4) CALL debug_write(npts, 'before allocation', &
965	circ_class_biomass, circ_class_n, circ_class_kill, &
966	plant_status, veget_cov_max)
967
968	!-
969	CALL growth_fun_all (npts, dt_days, veget_cov_max, veget_cov, PFTpresent, &
970	plant_status, when_growthinit, t2m_week, &
971	nstress_season, moiavail_growingseason, &
972	gpp_daily, resp_maint_part, resp_maint, &
973	resp_growth, npp_daily, bm_alloc, age, &
974	leaf_age, leaf_frac, use_reserve, &
975	lab_fac, rue_longterm, circ_class_n, &
976	circ_class_biomass, KF, sigma, &
977	gammas, tau_eff_leaf, tau_eff_sap, tau_eff_root, &
978	k_latosa_adapt, forest_managed, circ_class_dist, &
979	cn_leaf_min_season, atm_to_bm, store_sum_delta_ba)
980
981	check_intern(:,:,iatm2land,icarbon) = gpp_daily(:,:) * veget_cov_max(:,:) * dt_days
982	check_intern(:,:,iland2atm,icarbon) = &
983	-un(resp_maint(:,:) + resp_growth(:,:)) veget_cov_max(:,:) * dt_days
984
985
986	!We store intermediate veget_cov_max and atm_to_bm in order to calculate
987	!the mass balance closure at the end odf the subroutine
988	veget_cov_max_hist(:,:,igro) = veget_cov_max(:,:)
989	atm_to_bm_hist(:,:,igro,:) = atm_to_bm(:,:,:)
990	atm_to_bm(:,:,:)=zero
991
992	!! 6. Land cover change and land management
993
994	! We only want to do management at the end of the year. Whether
995	! LCC takes place before land management appears to be a matter
996	! of taste but given that LCC results in the (partial) destruction
997	! of the PFT it appears somehow logic to deal with LCC before FM.
998	! After all, it seems unlikely that a land owner will first
999	! carefully thin a forest and then cut it to turn it into a
1000	! grassland or cropland.
1001	IF ( EndOfYear) THEN
1002
1003	! Debug
1004	!printlev_loc=get_printlev('stomate_lpj')
1005	IF (printlev_loc>=4) CALL debug_write(npts,'before age class distr', &
1006	circ_class_biomass, circ_class_n, circ_class_kill, &
1007	plant_status, veget_cov_max)
1008	!-
1009
1010	! Following allocation, the trees have grown and it needs to be
1011	! checked whether they are still in the correct age class. If not,
1012	! redistribute the age classes in line with a prescribe diameter
1013	! threshold for each age class. This function could be applied on
1014	! daily time step but once per year is probably enough. This
1015	! routines distributes biomass within a species group so it can
1016	! be used in combination with species changes.
1017	CALL age_class_distr(npts, circ_class_n, circ_class_biomass, &
1018	veget_cov_max, circ_class_dist, wstress_season, &
1019	lm_lastyearmax, age, leaf_frac, atm_to_bm, &
1020	everywhere, litter, som, lignin_struc, &
1021	lignin_wood, bm_to_litter,turnover_daily, PFTpresent, when_growthinit,&
1022	forest_managed, KF, plant_status, &
1023	npp_longterm, gpp_daily, leaf_age, lab_fac, &
1024	gdd_from_growthinit, gdd_midwinter, time_hum_min, &
1025	hum_min_dormance, gdd_m5_dormance, &
1026	ncd_dormance, moiavail_month, moiavail_week, ngd_minus5, &
1027	resp_maint, resp_growth, npp_daily, &
1028	use_reserve, rue_longterm, mai, pai, &
1029	mai_count, previous_wood_volume, moiavail_growingseason,&
1030	matrixA, matrixV, VectorB, VectorU, age_stand, last_cut, &
1031	k_latosa_adapt, fm_change_map, lpft_replant, cn_leaf_min_season,&
1032	nstress_season, soil_n_min, p_O2, bact, store_sum_delta_ba, &
1033	CN_som_litter_longterm,nbp_pool_start)
1034
1035	! We store intermediate veget_cov_max and atm_to_bm in order to calculate
1036	! the mass balance closure at the end odf the subroutine
1037	veget_cov_max_hist(:,:,iage) = veget_cov_max(:,:)
1038	atm_to_bm_hist(:,:,iage,:) = atm_to_bm(:,:,:)
1039	atm_to_bm(:,:,:)=zero
1040
1041
1042	!! 5.2 Land cover change
1043	IF (do_now_stomate_lcchange) THEN
1044
1045	! When there is a loss in veget_cov_max (i.e., oaks are changed
1046	! to grasslands the) standard code will distribute this loss
1047	! proportionaly (i.e., all age classes with oak will be partly
1048	! converted) to the populated age classes of the species group.
1049	! In case of a species change we only want to change the surface
1050	! area of the youngest age class (because we will only plant
1051	! trees in the youngest age class). set losses = proportional
1052	losses = 'proportional'
1053	CALL land_cover_change_main(npts, dt_days, veget_cov_max, veget_cov, &
1054	harvest_pool, harvest_type, harvest_cut, harvest_area, &
1055	litter, som, lignin_struc, lignin_wood, &
1056	PFTpresent, everywhere, when_growthinit, &
1057	leaf_frac, circ_class_n, circ_class_biomass, &
1058	atm_to_bm, forest_managed, &
1059	KF, wstress_season, wstress_month, plant_status, npp_longterm, &
1060	age, lm_lastyearmax, harvest_pool_bound, resolution, &
1061	bm_to_litter, turnover_daily,leaf_age, circ_class_dist, &
1062	tau_eff_leaf, tau_eff_sap, &
1063	tau_eff_root, veget_cov_max_new, age_stand, last_cut, &
1064	k_latosa_adapt, losses, light_tran_to_level_season, &
1065	EndOfYear, lpft_replant, soil_n_min, bact, species_change_map)
1066
1067	! Set the flag done_stomate_lcchange to be used in the end of sechiba_main
1068	! to update the fractions.
1069	do_now_stomate_lcchange=.FALSE.
1070	done_stomate_lcchange=.TRUE.
1071
1072	!We store intermediate veget_cov_max and atm_to_bm in order to
1073	!calculate
1074	!the mass balance closure at the end odf the subroutine
1075	atm_to_bm_hist(:,:,ilcc,:) = atm_to_bm(:,:,:)
1076	veget_cov_max_hist(:,:,ilcc) = veget_cov_max(:,:)
1077	atm_to_bm(:,:,:) = zero
1078
1079	ENDIF
1080
1081	!! 5.3 Forest management
1082	! Thin and harvest species. Flag opportunities for species
1083	! changes.
1084
1085	! Debug
1086	!printlev_loc=get_printlev('stomate_lpj')
1087	IF (printlev_loc>=4) CALL debug_write(npts,'before forestry', &
1088	circ_class_biomass, circ_class_n, circ_class_kill, &
1089	plant_status, veget_cov_max)
1090	!-
1091	CALL forestry (npts, age_stand, last_cut, &
1092	circ_class_n, circ_class_kill, forest_managed, rdi,&
1093	circ_class_biomass, mai, pai, previous_wood_volume, &
1094	mai_count, coppice_dens, veget_cov_max, rdi_target_upper,&
1095	rdi_target_lower, fm_change_map, species_change_map)
1096
1097	!! 5.4 Crop harvest
1098	! According to the logic of this module crop harvest should
1099	! take place here but it depends on :: turnover_daily and that
1100	! variable is not calculated until the call to turnover. Also
1101	! the implemented approach harvests biomass daily rather once
1102	! per year (as it should be).
1103
1104	!! 5.5 Grazing
1105	! This is an excellent place to account for grazing.
1106
1107	!! 5.6 Litter raking
1108	! If you are NOT doing historical runs, litter raking almost certainly
1109	! does not apply to you.
1110	!
1111	! In order to feed animals through the winter, as well as absorb
1112	! their manure to spread on the fields during the spring, farmers
1113	! in historical Europe would collect litter from forests and place
1114	! it in the stables during winter. We will simulate this by
1115	! looking at the litter pools at the end of every year and moving
1116	! around a certain amount of litter from the forest to the crop PFTs.
1117	! We do this after the forest management under the assumption that
1118	! branches left on site after thinning operations would, historically,
1119	! have been collected for this purpose. Of course, litter raking fell
1120	! out of fashion in the mid-19th century, around the same time that
1121	! scientific sylviculture became popular, so perhaps it doesn't
1122	! matter.
1123	IF(ok_litter_raking)THEN
1124
1125	! Calculate fluxes from litter raking
1126	CALL litter_raking(npts, veget_cov_max, resolution, litter_demand, &
1127	litter, lrake_frac)
1128
1129	ENDIF
1130
1131	ENDIF ! EndOfYear
1132
1133	!! 6. Fire, Wind and Pests
1134
1135	!! 6.1 Call wind_damage module
1136	IF (ok_windthrow) THEN
1137
1138	! Calculate the critical wind speed for uprooting and breakage
1139	! critical windspeeds are calculated every half-hour. A
1140	! cummulated critical windspeed is then passed to stomate_lpj
1141	! to calculate the actual stand damage and mortality
1142	! IF(ld_windfall)THEN
1143	WRITE(numout,*) 'Call wind_damage module for calculating critical wind speed'
1144	WRITE(numout,*) 'Daily maximum windspeed at pixel level:', &
1145	wind_speed_daily
1146	WRITE(numout,*) 'Daily maximum soil temperature near 80 cm below ground at pixel level:',&
1147	soil_temp_daily
1148
1149	! ENDIF
1150
1151	CALL wind_damage(npts, nlevels_tot, circ_class_biomass, &
1152	veget_cov_max, circ_class_n, wind_speed_daily, plant_status, &
1153	resolution, circ_class_kill, harvest_5y_area, forest_managed, &
1154	soil_temp_daily)
1155
1156	ENDIF
1157	!! 6.2 Call the fire module
1158
1159
1160	!! 7. Kill PFT's
1161
1162	! Kill slow growing PFTs in DGVM or STOMATE with dynamic or
1163	! constant mortality. Mark trees that died from self-thinning
1164	! or environmental mortality. The biomass pools are not touched.
1165	! Vegetation is only marked for killing
1166
1167	! Debug
1168	!printlev_loc=get_printlev('stomate_lpj')
1169	IF (printlev_loc>=4) CALL debug_write(npts, 'before mark_to_kill', &
1170	circ_class_biomass, circ_class_n, circ_class_kill, &
1171	plant_status, veget_cov_max)
1172	!-
1173	CALL mark_to_kill (npts, 'growth ', lm_lastyearmax, &
1174	PFTpresent, npp_longterm, circ_class_biomass, &
1175	circ_class_n, circ_class_kill, rue_longterm, turnover_longterm, &
1176	dt_days, veget_cov_max, st_dist, forest_managed)
1177
1178
1179	! +++CHECK+++
1180	! Uncomment when working on the DGVM
1181	!!$ IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN
1182	!!$ CALL kill (npts, 'npp ', lm_lastyearmax, &
1183	!!$ ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
1184	!!$ lai, age, leaf_age, leaf_frac, npp_longterm, &
1185	!!$ when_growthinit, everywhere, veget_cov_max, bm_to_litter)
1186	!!$
1187	!!$ !! 6.2.1 Update wood biomass
1188	!!$ ! For the DGVM
1189	!!$ IF(ok_dgvm) THEN
1190	!!$ WHERE (ind(:,:).GT.min_stomate)
1191	!!$ woodmass_ind(:,:) = &
1192	!!$ ((biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
1193	!!$ + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon)) &
1194	!!$ *veget_cov_max(:,:))/ind(:,:)
1195	!!$ ENDWHERE
1196	!!$
1197	!!$ ! For all pixels with individuals
1198	!!$ ELSE
1199	!!$ WHERE (ind(:,:).GT.min_stomate)
1200	!!$ woodmass_ind(:,:) = &
1201	!!$ (biomass(:,:,isapabove,icarbon) + biomass(:,:,isapbelow,icarbon) &
1202	!!$ + biomass(:,:,iheartabove,icarbon) + biomass(:,:,iheartbelow,icarbon))/ind(:,:)
1203	!!$ ENDWHERE
1204	!!$ ENDIF ! ok_dgvm
1205	!!$
1206	!!$ !! 6.2.2 New crown area and maximum vegetation cover after growth
1207	!!$ CALL crown (npts, PFTpresent, &
1208	!!$ ind, biomass, woodmass_ind,&
1209	!!$ veget_cov_max, cn_ind, height)
1210	!!$
1211	!!$ ENDIF ! ok_dgvm
1212	! +++++++++++
1213
1214	!! 7. fire
1215
1216	!! 7.1. Burn PFTs
1217	! +++CHECK+++
1218	! Uncomment when working on the DGVM
1219	!!$ CALL fire (npts, dt_days, &
1220	!!$ litterhum_daily, t2m_daily, lignin_struc, lignin_wood, veget_cov_max, &
1221	!!$ fireindex, firelitter, biomass, ind, &
1222	!!$ litter, dead_leaves, bm_to_litter, &
1223	!!$ co2_fire, matrixA)
1224	! +++++++++++
1225
1226	!! 7.2 Kill PFTs in DGVM
1227	! +++CHECK+++
1228	! Uncomment when working on the DGVM
1229	!!$ IF ( ok_dgvm ) THEN
1230	!!$
1231	!!$ ! reset attributes for eliminated PFTs
1232	!!$ CALL kill (npts, 'fire ', lm_lastyearmax, &
1233	!!$ ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
1234	!!$ lai, age, leaf_age, leaf_frac, npp_longterm, &
1235	!!$ when_growthinit, everywhere, veget_cov_max, bm_to_litter)
1236	!!$
1237	!!$ ENDIF ! ok_dgvm
1238	! +++++++++++
1239
1240	!! 8. Tree mortality
1241
1242	! Here is where the trees actually die and biomass is moved around.
1243	! First, we kill all the trees that died because of human intervention.
1244	! When we mark trees for killing the reason of their dead is also
1245	! stored in circ_class_kill. This is where we make use of that
1246	! information.
1247
1248	!+++CHECK+++
1249	! This is called always. It is not clear whether this relates to
1250	! the DGVM or not. The temperature dependencies suggest it does.
1251	! check when working on the DGVM. It is not needed for when the
1252	! PFTs are prescribed.
1253	!!$ ! Does not depend on age, therefore does not change crown area.
1254	!!$ CALL gap (npts, dt_days, &
1255	!!$ npp_longterm, turnover_longterm, lm_lastyearmax, &
1256	!!$ PFTpresent, t2m_min_daily, Tmin_spring_time, &
1257	!!$ biomass, ind, bm_to_litter, mortality)
1258	!+++++++++++
1259
1260	! Debug
1261	!printlev_loc=get_printlev('stomate_lpj')
1262	IF (printlev_loc>=4) CALL debug_write(npts,'before anthrop. mortality', &
1263	circ_class_biomass, circ_class_n, circ_class_kill, &
1264	plant_status, veget_cov_max)
1265	!-
1266	CALL anthropogenic_mortality (npts, &
1267	bm_to_litter, circ_class_biomass, circ_class_kill, &
1268	circ_class_n, harvest_pool_bound, harvest_pool, harvest_type, &
1269	harvest_cut, harvest_area, resolution, &
1270	veget_cov_max, age_stand, last_cut, mai_count, circ_class_dist, &
1271	coppice_dens, plant_status, biomass_cut)
1272
1273	!+++CHECK+++
1274	! Code for the DGVM - in addition to the killing
1275	!!$ IF ( ok_dgvm ) THEN
1276	!!$ ! reset attributes for eliminated PFTs
1277	!!$ CALL kill (npts, 'gap ', lm_lastyearmax, &
1278	!!$ ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
1279	!!$ lai, age, leaf_age, leaf_frac, npp_longterm, &
1280	!!$ when_growthinit, everywhere, veget_cov_max, bm_to_litter)
1281	!!$ ENDIF
1282	!+++++++++++
1283
1284	! Now we kill all the plants which died due to natural causes. This comes
1285	! after the plant died due to human causes, since in theory a forest
1286	! that is thinned will not suffer as high of natural mortality.
1287	! When we mark trees for killing the reason of their dead is also
1288	! stored in circ_class_kill. This is where we make use of that information.
1289
1290	! Debug
1291	!printlev_loc=get_printlev('stomate_lpj')
1292	IF (printlev_loc>=4) CALL debug_write(npts, 'before natural mortality', &
1293	circ_class_biomass, circ_class_n, circ_class_kill, &
1294	plant_status, veget_cov_max)
1295	!-
1296
1297	CALL natural_mortality (npts, bm_to_litter, &
1298	circ_class_biomass, circ_class_kill, circ_class_n, &
1299	veget_cov_max, biomass_cut)
1300
1301	!! 9. Change forestry
1302
1303	! PFTs have been managed and killed based on the outcome of
1304	! the operations and processes it is now decided whether we
1305	! will replant the PFT possible with a different species.
1306	IF(ok_change_species)THEN
1307	CALL flag_species_change(npts, veget_cov_max, circ_class_biomass,&
1308	circ_class_n,lpft_replant, forest_managed)
1309	ENDIF
1310
1311	!! 10. Clean the stands after mortality
1312
1313	! Now we need to do some checking. If there is no biomass on the site,
1314	! some variables have to be reset. If some of the circumference classes
1315	! are empty, we need to redistribute the biomass among them. Notice this
1316	! must be done after both the natural and human killing is done, because
1317	! circ_class_kill for the natural killing is calculated based on
1318	! circ_class_biomass, and if we redistribute we change circ_class_biomass.
1319	! If a PFT is found empty (= killed) it is getting replanted unless
1320	! species change is used. In that case we wait for the end of the year
1321	! to replant with a different species.
1322
1323	! Debug
1324	!printlev_loc=get_printlev('stomate_lpj')
1325	IF (printlev_loc>=4) CALL debug_write(npts, 'before mortality_clean', &
1326	circ_class_biomass, circ_class_n, circ_class_kill, &
1327	plant_status, veget_cov_max)
1328
1329	!-
1330
1331	CALL mortality_clean (npts, circ_class_biomass, &
1332	circ_class_n, age, everywhere, leaf_age,&
1333	leaf_frac, plant_status, when_growthinit, circ_class_dist,&
1334	age_stand, PFTpresent, last_cut, mai_count, &
1335	veget_cov_max, npp_longterm, KF, atm_to_bm, &
1336	wstress_season, lm_lastyearmax, lignin_struc, lignin_wood, &
1337	som, litter, dt_days, forest_managed, bm_to_litter, lab_fac, &
1338	gpp_daily, resp_maint,&
1339	resp_growth, npp_daily, use_reserve, mai, pai,&
1340	rue_longterm, previous_wood_volume,species_change_map,&
1341	tau_eff_leaf, tau_eff_sap, tau_eff_root, moiavail_growingseason,&
1342	k_latosa_adapt, lpft_replant, cn_leaf_min_season, &
1343	nstress_season, soil_n_min, plant_n_uptake_daily, p_O2, bact, &
1344	CN_som_litter_longterm, matrixA, matrixV, vectorB, vectorU)
1345
1346	atm_to_bm_hist(:,:,imor,:) = atm_to_bm(:,:,:)
1347	veget_cov_max_hist(:,:,imor) = veget_cov_max(:,:)
1348	atm_to_bm(:,:,:) = zero
1349
1350	!! 11. Leaf senescence and other turnover processes
1351
1352	! Calculate all turnover processes which are responsible
1353	! for moving carbon/nitrogen to bifferent living biomass pools in
1354	! the case of sapwood to heartwood turnover or moving carbon/
1355	! nitrogen from the living biomass to the litter pools. The
1356	! turnover module also deals with leaf fall in autumn.
1357
1358	! Debug
1359	!printlev_loc=get_printlev('stomate_lpj')
1360	IF (printlev_loc>=4) CALL debug_write(npts,'before turnover', &
1361	circ_class_biomass, circ_class_n, circ_class_kill, &
1362	plant_status, veget_cov_max)
1363	!-
1364
1365	CALL turn_over (npts, dt_days, PFTpresent, herbivores, &
1366	maxmoiavail_lastyear, minmoiavail_lastyear, moiavail_week, &
1367	moiavail_month, t2m_longterm, t2m_month, t2m_week, &
1368	veget_cov_max, gdd_from_growthinit, leaf_age, leaf_frac, age, &
1369	turnover_daily, plant_status, turnover_time, &
1370	circ_class_biomass, circ_class_n, &
1371	when_growthinit, tau_eff_leaf, tau_eff_sap, tau_eff_root, &
1372	harvest_pool, harvest_type, harvest_cut, harvest_area, &
1373	resolution, wstress_month)
1374
1375	!! 12. Product use
1376
1377	! +++CHECK+++
1378	! Write the harvest pools to the history files here because they will
1379	! be reset when the subroutine basic_product_use is called. Note that
1380	! due to some peculiarities in the coding the units of the harvest
1381	! pool differ for grasses, crops and forests. For forest, the
1382	! harvest_pool is filled and emptied the last day of the year. Hence
1383	! all daily values are 0, except for the last day. When writing,
1384	! for example, monthly output files the mean monthly value will be
1385	! written. To obtain the absolute harvest, the value in the history
1386	! file should be multiplied by 30 or 31 depending on whether a 360 or
1387	! 365 day calendre is used. For grasslands the harvest_pool is filled
1388	! daily with the turnover. So the mean value written to the history file
1389	! is the correct value. Crops are typically harvested before December
1390	! 31st but the harvest_pool is only emptied on December 31st (just as
1391	! for grasses and forests). So between the harvest date and December
1392	! 31st the same biomass remains in the harvest pools. If the harvest
1393	! takes place before December 1st the same average value is written for
1394	! every day of December and therefore the last month contains the correct
1395	! value for harvest.
1396	! All of this would become straightforward if we would only write
1397	! harvest_pool the last day of the year. We tried using the operator
1398	! once in histwrite but were unhappy with the outcome. The IPSL team
1399	! no longer maintains histwrite as it should get replaced by the
1400	! ioxserver.
1401	CALL histwrite_p (hist_id_stomate, 'HARVEST_POOL_C', itime, &
1402	SUM(harvest_pool(:,:,:,icarbon),3), npts*nvm, horipft_index)
1403	CALL histwrite_p (hist_id_stomate, 'HARVEST_POOL_N', itime, &
1404	SUM(harvest_pool(:,:,:,initrogen),3), npts*nvm, horipft_index)
1405	CALL histwrite_p (hist_id_stomate, 'HARVEST_TYPE', itime, &
1406	harvest_type(:,:), npts*nvm, horipft_index)
1407	CALL histwrite_p (hist_id_stomate, 'HARVEST_CUT', itime, &
1408	harvest_cut(:,:), npts*nvm, horipft_index)
1409	CALL histwrite_p (hist_id_stomate, 'HARVEST_AREA', itime, &
1410	harvest_area(:,:), npts*nvm, horipft_index)
1411
1412	CALL xios_orchidee_send_field("HARVEST_POOL_C",SUM(harvest_pool(:,:,:,icarbon),3))
1413	CALL xios_orchidee_send_field("HARVEST_POOL_N",SUM(harvest_pool(:,:,:,initrogen),3))
1414	CALL xios_orchidee_send_field("HARVEST_TYPE",harvest_type(:,:))
1415	CALL xios_orchidee_send_field("HARVEST_CUT",harvest_cut(:,:))
1416	CALL xios_orchidee_send_field("HARVEST_AREA",harvest_area(:,:))
1417	!+++++++++++
1418
1419	! Agricultural harvest now goes into the short-lived product pool.
1420	! Hence, product use can only be estimated after crop_harvest which
1421	! in turn depends on the calculation of turnover. In previous versions
1422	! product use was an integral part of LCC. Maintaining this approach
1423	! would result in inconsistencies with harvest from forestry and
1424	! agriculture. Note that there is still no real harvest from grasslands.
1425	! Grasslands are thus wrongly considered natural PFTs without any human
1426	! use. To partly overcome this issue the turnover in grasslands is moved
1427	! into the harvest pool but this approach should be considered as a
1428	! patch rather than a permanent solution.
1429	IF(EndofYear)THEN
1430
1431	! Debug
1432	!printlev_loc=get_printlev('stomate_lpj')
1433	IF (printlev_loc>=4) CALL debug_write(npts,'before wood use', &
1434	circ_class_biomass, circ_class_n, circ_class_kill,&
1435	plant_status, veget_cov_max)
1436	!-
1437
1438	IF(ok_dimensional_product_use)THEN
1439
1440	! Wood harvest with small dimensions is being used as fire wood
1441	! the remaining wood goes into the short, medium and long product
1442	! pools. This implies that when the dimensions of the harvest
1443	! change, the product pools will change as well
1444	CALL dim_product_use(npts, dt_days, harvest_pool_bound, harvest_pool, &
1445	harvest_type, harvest_cut, harvest_area, &
1446	flux_prod_s, flux_prod_m, flux_prod_l, prod_s, &
1447	prod_m, prod_l, prod_s_total, prod_m_total, &
1448	prod_l_total, flux_prod_total, flux_s, flux_m, &
1449	flux_l, resolution, veget_cov_max)
1450
1451	ELSE
1452
1453	! The scheme does not account for the harvest dimensions. All
1454	! wood is distributed across the short, medium and long product
1455	! pools. This basically means that wood use in 1600, 1700, 1800,
1456	! 1900 and 2000 was identical.
1457	CALL basic_product_use(npts, dt_days, harvest_pool_bound, harvest_pool, &
1458	harvest_type, harvest_cut, harvest_area, &
1459	flux_prod_s, flux_prod_m, flux_prod_l, prod_s, &
1460	prod_m, prod_l, prod_s_total, prod_m_total, &
1461	prod_l_total, flux_prod_total, flux_s, flux_m, &
1462	flux_l, resolution, veget_cov_max)
1463
1464	ENDIF ! dimensional product use
1465
1466	ELSE ! EndOfYear
1467
1468	! These variables are all intent(out) so they should be initialized
1469	! if it's not the end of the year
1470	CALL product_init(flux_prod_s, flux_prod_m, flux_prod_l, &
1471	prod_s_total, prod_m_total, prod_l_total, flux_s, &
1472	flux_m, flux_l, flux_prod_total)
1473
1474	ENDIF
1475
1476	!! 13. Establishment of saplings/recruitment
1477
1478	IF(ok_recruitment .AND. EndOfYear)THEN
1479
1480	! Debug
1481	IF(printlev_loc.GT.4) WRITE(numout,*) 'We will call recruitment from stomate_lpj.f90 '
1482	!-
1483
1484	!! 12.1 Call recruitment process in
1485	CALL prescribe (npts, veget_cov_max, dt_days, PFTpresent, &
1486	everywhere, when_growthinit, leaf_frac, circ_class_dist, &
1487	circ_class_n, circ_class_biomass, atm_to_bm, &
1488	forest_managed, KF, plant_status, age, npp_longterm, &
1489	lm_lastyearmax, tau_eff_leaf, tau_eff_sap, tau_eff_root, &
1490	k_latosa_adapt, light_tran_to_level_season, &
1491	EndOfYear, lpft_replant, species_change_map)
1492
1493	!We store intermediate veget_cov_max and atm_to_bm in order to calculate
1494	!the mass balance closure at the end odf the subroutine
1495	atm_to_bm_hist(:,:,irec,:) = atm_to_bm(:,:,:)
1496	veget_cov_max_hist(:,:,irec) = veget_cov_max(:,:)
1497	atm_to_bm(:,:,:) = zero
1498
1499	ENDIF ! establishment/recruitment
1500
1501	!! 13. Establishment of saplings
1502
1503	!! IF (.NOT. ok_constant_mortality) THEN
1504
1505	! Establish new plants based on the functional allocation
1506	! +++CHECK+++
1507	! The comments below were written during the first round of
1508	! ORCHIDEE-CAN code development. In the mean time a recruitment
1509	! module was added to prescribe. Some the functionality
1510	! that was aimed for in the initial establish routine can now
1511	! be accounted for in prescribe.
1512
1513	! Old comments:
1514	! How does this intefer with self-thinning it seems that this
1515	! code tries to maintain a +/- high biomass by adding individuals
1516	! when the fpc (foliage) projected cover decreases). This is
1517	! exactely what we want to avoid as we want to introduce gaps
1518	! in the static approach. We could try to adjust the
1519	! proposed aproach for unmanaged systems but I doubt whether
1520	! the rule of Deleuze and Dhote could deal with such an extreme
1521	! range of diameter sizes. Most unmanged forest are, however,
1522	! under a disturbance regime and may be aproximated
1523	! by an 'even-age' evolution.
1524
1525	! NOTE: We have decided not to create this routine. The
1526	! reasoning is that we really want gaps in the system and don't
1527	! want to continually establish. All of this should be taken
1528	! care of in prescribe at the end
1529	! of the year. If someone wants to run with DGVM and the new
1530	! allocation, this will require some effort here.
1531
1532	!CALL establish_prognostic (npts, dt_days, PFTpresent, regenerate, &
1533	! neighbours, resolution, need_adjacent, herbivores, &
1534	! precip_lastyear, gdd0_lastyear, lm_lastyearmax, avail_tree, &
1535	! avail_grass, npp_longterm, leaf_age, leaf_frac, &
1536	! ind, biomass, age, everywhere, &
1537	! co2_to_bm, veget_cov_max, circ_class_biomass, circ_class_n)
1538	! ++++++++++++
1539
1540	!! ELSE
1541
1542	! Plant establishment does not occur when constant mortality
1543	! is used outside the DGVM. Write the values of the flags to
1544	! make sure this condition is satisfied
1545	!! IF (firstcall) THEN
1546
1547	!! WRITE(numout,*) 'Plant establishment does not occur because '
1548	!! WRITE(numout,*) 'constant mortality is used outside the DGVM'
1549	!! WRITE(numout,*) 'control%ok_dgvm, ', ok_dgvm
1550	!! WRITE(numout,*) 'control%ok_constant_mortality, ', &
1551	!! ok_constant_mortality
1552	!!
1553	!! ENDIF
1554
1555	!! ENDIF ! constant mortality
1556
1557	!! 14. Replant PFTs with a different species.
1558
1559	IF(EndOfYear)THEN
1560
1561	IF(ok_change_species)THEN
1562
1563	IF(printlev_loc>4) WRITE(numout,*) 'species change lpft_replant'
1564
1565	! This is an ad-hoc method to change the species after they die.
1566	! It is only meant to change species through human intervention.
1567	! In theory, this is a function that should be done in a DGVM, which
1568	! is why this code is considered only temporary. Because we passed
1569	! the mortality routines the PFTs may have changed. For example, in
1570	! sapiens_forestry the oldest age class was marked for harvest. In
1571	! kill the biomass is removed and the surface area taken by the
1572	! oldest age class was moved to the youngest age class of that group
1573	! Here we will move the youngest age class of the intial group to
1574	! the younest age class of the new species group. First create
1575	! the vectors that are then used to calculate the losses and gains
1576	! in surface area (veget_cov_max)
1577	CALL species_change(npts, lpft_replant, veget_cov_max_new, &
1578	forest_managed, species_change_map, veget_cov_max, fm_change_map)
1579
1580	! Use the land cover change code to move the biomass to the
1581	! correct PFT. When veget_cov_max is lost (i.e., oaks are changed
1582	! to grasslandsthe) standard code will distribute this loss
1583	! proportionaly (i.e., all age classes with oak will be partly
1584	! coverted) to the populated age classes of the species group.
1585	! In case of species change we only want to change the surface
1586	! area of the youngest age class (because we will only plant
1587	! trees in the youngest age class). set losses = youngest
1588	losses = 'youngest'
1589
1590	CALL land_cover_change_main (npts, dt_days, veget_cov_max, &
1591	veget_cov, harvest_pool, &
1592	harvest_type, harvest_cut, harvest_area, &
1593	litter, som, lignin_struc, lignin_wood, &
1594	PFTpresent, everywhere, when_growthinit, &
1595	leaf_frac, circ_class_n, circ_class_biomass, &
1596	atm_to_bm, forest_managed, &
1597	KF, wstress_season, wstress_month, plant_status, npp_longterm, &
1598	age, lm_lastyearmax, harvest_pool_bound, resolution, &
1599	bm_to_litter, turnover_daily,leaf_age, circ_class_dist, &
1600	tau_eff_leaf, tau_eff_sap, &
1601	tau_eff_root, veget_cov_max_new, age_stand, last_cut, &
1602	k_latosa_adapt, losses, light_tran_to_level_season, &
1603	EndOfYear, lpft_replant, soil_n_min, bact, species_change_map)
1604
1605	veget_cov_max_hist(:,:,ispc) = veget_cov_max(:,:)
1606	atm_to_bm_hist(:,:,ispc,:) = atm_to_bm(:,:,:)
1607	atm_to_bm(:,:,:) = zero
1608
1609	ENDIF ! checking for species change
1610
1611	ENDIF ! checking for the end of the year
1612
1613	!! 11. Light competition
1614	! +++CHECK+++
1615	! Uncomment when working on the DGVM
1616
1617	!!$ !! If not using constant mortality then kill with light competition
1618	!!$! IF ( ok_dgvm .OR. .NOT.(lpj_gap_const_mort) ) THEN
1619	!!$ IF ( ok_dgvm ) THEN
1620	!!$
1621	!!$ !! 11.1 Light competition
1622	!!$ CALL light (npts, dt_days, &
1623	!!$ veget_cov_max, fpc_max, PFTpresent, cn_ind, maxfpc_lastyear, &
1624	!!$ lm_lastyearmax, ind, biomass, veget_lastlight, bm_to_litter, mortality)
1625	!!$
1626	!!$ !! 11.2 Reset attributes for eliminated PFTs
1627	!!$ CALL kill (npts, 'light ', lm_lastyearmax, &
1628	!!$ ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
1629	!!$ lai, age, leaf_age, leaf_frac, npp_longterm, &
1630	!!$ when_growthinit, everywhere, veget_cov_max, bm_to_litter)
1631	!!$
1632	!!$ ENDIF
1633	!!$
1634	!!$
1635	!!$
1636	!!$ !! 12. Establishment of saplings
1637	!!$
1638	!!$ IF ( ok_dgvm .OR. .NOT.lpj_gap_const_mort ) THEN
1639	!!$
1640	!!$ !! 12.1 Establish new plants
1641	!!$ CALL establish (npts, dt_days, PFTpresent, regenerate, &
1642	!!$ neighbours, resolution, need_adjacent, herbivores, &
1643	!!$ precip_lastyear, gdd0_lastyear, lm_lastyearmax, &
1644	!!$ cn_ind, lai, avail_tree, avail_grass, npp_longterm, &
1645	!!$ leaf_age, leaf_frac, &
1646	!!$ ind, biomass, age, everywhere, atm_to_bm, &
1647	!!$ soil_n_min, plant_n_uptake_daily, nstress_season, veget_cov_max, woodmass_ind, &
1648	!!$ mortality, bm_to_litter)
1649	!!$
1650	!!$ IF (printlev_loc>=3) WRITE (numout,*) 'after establish soil_n_min(test_grid,test_pft,:):',soil_n_min(test_grid,test_pft,:)
1651	!!$ !! 12.2 Calculate new crown area (and maximum vegetation cover)
1652	!!$ CALL crown (npts, PFTpresent, &
1653	!!$ ind, biomass, woodmass_ind, &
1654	!!$ veget_cov_max, cn_ind, height)
1655	!!$
1656	!!$ ENDIF
1657	! +++++++++++
1658
1659
1660	!! 13. Calculate final LAI and vegetation cover
1661	! +++CHECK+++
1662	! Move to DGVM - not needed when using prescribed veget_cov_max
1663	!!$ CALL cover (npts, cn_ind, ind, biomass, &
1664	!!$ veget_cov_max, veget_cov_max_tmp, lai, &
1665	!!$ litter, som, turnover_daily, bm_to_litter, &
1666	!!$ atm_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, gpp_daily, &
1667	!!$ lignin_struc, lignin_wood, soil_n_min)
1668	! +++++++++++
1669
1670
1671
1672	!! 16. Calculate vcmax
1673
1674	CALL vmax (npts, dt_days, leaf_age, leaf_frac, vcmax, nue)
1675
1676	!! 17. Compute the effective leaf area index.
1677
1678	! Need to put in a function here to compute the height of the
1679	! photosynthesis levels given the height of the energy levels and the
1680	! vegetation on the grid square. The hightest levels will be a
1681	! function of the height of the vegetation so that we don't waste
1682	! computational time on empty levels.
1683
1684	! Debug
1685	! The ipslerr_p issue is most likely causing the crash here.
1686	! printlev_loc=get_printlev('stomate_lpj')
1687	! IF (printlev_loc>=4)
1688	! CALL debug_write(npts,'before lai effective', &
1689	! circ_class_biomass, circ_class_n, circ_class_kill, &
1690	! plant_status, veget_cov_max)
1691
1692	!-
1693	CALL calculate_z_level_photo(npts, circ_class_biomass, circ_class_n, &
1694	z_level_photo)
1695
1696	! Finding the true LAI per level is different from finding the
1697	! effective LAI per level, so we'll do that here. This only
1698	! changes once every day so hopefully it is not too expensive.
1699	! It will eventually be needed by the energy budget. It is
1700	! also needed by effective_lai for grasses and crops.
1701	CALL find_lai_per_level(npts, z_level_photo, &
1702	circ_class_biomass, circ_class_n, lai_per_level, &
1703	max_height_store)
1704
1705	! This is a function of the solar angle, and needed for the albedo
1706	! and the photosynthesis. It was very expensive to compute, so instead
1707	! of computing it at every timestep we compute it at a couple points
1708	! throughout the day and fit a function to it. We do this now because
1709	! the canopy as it appears now will be used for the next day.
1710	! Now we actually find the effective LAI and fit the function
1711	! we'll use later on
1712	CALL fitting_laieff(npts, z_level_photo, circ_class_biomass, &
1713	circ_class_n, veget_cov_max, lai_per_level, laieff_fit, &
1714	h_array_out, z_array_out)
1715
1716	!! 18. Check numerical consistency of this routine
1717
1718	! This test is always performed. If err_act.EQ.1 then
1719	! the value of the mass balance error -if any- is written
1720	! to the history file.
1721	!!$ printlev_loc=get_printlev('stomate_lpj')
1722
1723	! 14.2 Check surface area
1724	CALL check_surface_area("stomate_lpj", npts, veget_cov_max_hist(:,:,ibeg), &
1725	veget_cov_max,'pixel')
1726
1727	IF (printlev_loc .GE. 4) THEN
1728	DO ivm=1,nvm
1729	WRITE(numout,*) "step, veget_cov_max, atm_to_bm for ivm:", ivm
1730	WRITE(numout,*)"IBEG",veget_cov_max_hist(test_grid,ivm,ibeg),&
1731	atm_to_bm_hist(test_grid,ivm,ibeg,1)
1732	WRITE(numout,*)"IPRE",veget_cov_max_hist(test_grid,ivm,ipre),&
1733	atm_to_bm_hist(test_grid,ivm,ipre,1)
1734	WRITE(numout,*)"IPHE",veget_cov_max_hist(test_grid,ivm,iphe),&
1735	atm_to_bm_hist(test_grid,ivm,iphe,1)
1736	WRITE(numout,*)"IGRO",veget_cov_max_hist(test_grid,ivm,igro),&
1737	atm_to_bm_hist(test_grid,ivm,igro,1)
1738	WRITE(numout,*)"IAGE",veget_cov_max_hist(test_grid,ivm,iage),&
1739	atm_to_bm_hist(test_grid,ivm,iage,1)
1740	WRITE(numout,*)"ILCC",veget_cov_max_hist(test_grid,ivm,ilcc),&
1741	atm_to_bm_hist(test_grid,ivm,ilcc,1)
1742	WRITE(numout,*)"IMOR",veget_cov_max_hist(test_grid,ivm,imor),&
1743	atm_to_bm_hist(test_grid,ivm,imor,1)
1744	WRITE(numout,*)"IREC",veget_cov_max_hist(test_grid,ivm,irec),&
1745	atm_to_bm_hist(test_grid,ivm,irec,1)
1746	WRITE(numout,*)"ISPC",veget_cov_max_hist(test_grid,ivm,ispc),&
1747	atm_to_bm_hist(test_grid,ivm,ispc,1)
1748	ENDDO
1749	ENDIF
1750
1751	! 14.3 Mass balance closure
1752	! 14.3.1 Calculate final biomass
1753	pool_end(:,:,:) = zero
1754	DO iele = 1,nelements
1755	DO ipar = 1,nparts
1756	DO icir=1,ncirc
1757	pool_end(:,:,iele) = pool_end(:,:,iele) + &
1758	circ_class_biomass(:,:,icir,ipar,iele)* &
1759	circ_class_n(:,:,icir) * veget_cov_max(:,:)
1760	ENDDO
1761	pool_end(:,:,iele) = pool_end(:,:,iele) + &
1762	(bm_to_litter(:,:,ipar,iele) + &
1763	turnover_daily(:,:,ipar,iele)) * veget_cov_max(:,:)
1764	ENDDO
1765	pool_end(:,1,iele) = pool_end(:,1,iele) + &
1766	( SUM(SUM(harvest_pool(:,:,:,iele),3),2) + &
1767	SUM(prod_l(:,:,iele),2) + SUM(prod_m(:,:,iele),2) + &
1768	SUM(prod_s(:,:,iele),2) ) / pixel_area(:)
1769	ENDDO
1770
1771	! Common processes. Atm_to_bm should be empty when the model enters
1772	! stomate_lpj but its value can change at several call in stomate_lpj.
1773	! for this reason the history of atm_to_bm needs to be stored and
1774	! it has be multiplied by the history of veget_cov_max within stomate_lpj
1775	! to track all the carbon down.
1776	DO iele = 1,nelements
1777	check_intern(:,:,iatm2land,iele) = &
1778	check_intern(:,:,iatm2land,iele) + &
1779	SUM(atm_to_bm_hist(:,:,:,iele) * veget_cov_max_hist(:,:,:) * dt_days,3)
1780	check_intern(:,1,iland2atm,iele) = check_intern(:,1,iland2atm,iele) &
1781	-un * ( SUM(flux_s(:,:,iele),2) + SUM(flux_m(:,:,iele),2) + &
1782	SUM(flux_l(:,:,iele),2) ) / pixel_area(:)
1783	check_intern(:,:,ilat2out,iele) = zero
1784	check_intern(:,:,ilat2in,iele) = -un * zero
1785	check_intern(:,:,ipoolchange,iele) = -un * (pool_end(:,:,iele) - &
1786	pool_start(:,:,iele))
1787	ENDDO
1788
1789	closure_intern = zero
1790	DO ivm=1, nvm
1791	DO imbc = 1,nmbcomp
1792	DO iele=1,nelements
1793	! Debug
1794	IF (printlev_loc>=4) THEN
1795	IF(iele .EQ. 1) WRITE(numout,"(9999(G12.5,:,','))") &
1796	'imbc, ivm, check_intern', test_grid , imbc, &
1797	ivm, check_intern(test_grid,ivm,imbc,iele)
1798	ENDIF
1799	closure_intern(:,ivm,iele) = closure_intern(:,ivm,iele) + &
1800	check_intern(:,ivm,imbc,iele)
1801	ENDDO
1802	ENDDO
1803	ENDDO
1804	CALL check_mass_balance("stomate_lpj", closure_intern, npts, &
1805	pool_end, pool_start, veget_cov_max, 'pixel')
1806
1807	! Add to the mass balance errors of stomate_dt_sechiba
1808	DO iele = 1,nelements
1809	mass_balance_closure(:,iele) = mass_balance_closure(:,iele) + &
1810	closure_intern(:,1,iele)
1811	ENDDO
1812
1813
1814	!+++CHECK+++
1815	! This does not look essential. Just making daily totals.
1816	! Shouldn we do this where we calculate heterothropic respiration?
1817	!!$ !! 16. Total heterotrophic respiration
1818	!!$
1819	!!$ tot_soil_carb(:,:) = zero
1820	!!$ tot_litter_carb(:,:) = zero
1821	!!$ DO j=2,nvm
1822	!!$
1823	!!$ tot_litter_carb(:,j) = tot_litter_carb(:,j) + (litter(:,istructural,j,iabove,icarbon) + &
1824	!!$ & litter(:,imetabolic,j,iabove,icarbon) + litter(:,iwoody,j,iabove,icarbon) + &
1825	!!$ & litter(:,istructural,j,ibelow,icarbon) + litter(:,imetabolic,j,ibelow,icarbon)+ &
1826	!!$ & litter(:,iwoody,j,ibelow,icarbon))
1827	!!$
1828	!!$ tot_soil_carb(:,j) = tot_soil_carb(:,j) + (som(:,iactive,j,icarbon) + &
1829	!!$ & som(:,islow,j,icarbon)+ som(:,ipassive,j,icarbon)+ som(:,isurface,j,icarbon))
1830	!!$
1831	!!$ ENDDO
1832	!!$ tot_litter_soil_carb(:,:) = tot_litter_carb(:,:) + tot_soil_carb(:,:)
1833	!!$
1834	!!$! DO k = 1, nelements ! Loop over # elements
1835	!!$! tot_live_biomass(:,:,k) = biomass(:,:,ileaf,k) + biomass(:,:,isapabove,k) + biomass(:,:,isapbelow,k) +&
1836	!!$! & biomass(:,:,iheartabove,k) + biomass(:,:,iheartbelow,k) + &
1837	!!$! & biomass(:,:,iroot,k) + biomass(:,:,ifruit,k) + biomass(:,:,icarbres,k)
1838	!!$! END DO ! Loop over # elements
1839	!!$
1840	!!$ tot_live_biomass(:,:,:) = biomass(:,:,ileaf,:) + biomass(:,:,isapabove,:) + biomass(:,:,isapbelow,:) +&
1841	!!$ & biomass(:,:,iheartabove,:) + biomass(:,:,iheartbelow,:) + &
1842	!!$ & biomass(:,:,iroot,:) + biomass(:,:,ifruit,:) + biomass(:,:,icarbres,:) + biomass(:,:,ilabile,:)
1843	!!$
1844	!!$
1845	!!$ tot_turnover(:,:,:) = turnover_daily(:,:,ileaf,:) + turnover_daily(:,:,isapabove,:) + &
1846	!!$ & turnover_daily(:,:,isapbelow,:) + turnover_daily(:,:,iheartabove,:) + &
1847	!!$ & turnover_daily(:,:,iheartbelow,:) + turnover_daily(:,:,iroot,:) + &
1848	!!$ & turnover_daily(:,:,ifruit,:) + turnover_daily(:,:,icarbres,:) + turnover_daily(:,:,ilabile,:)
1849	!!$
1850	!!$ tot_bm_to_litter(:,:,:) = bm_to_litter(:,:,ileaf,:) + bm_to_litter(:,:,isapabove,:) +&
1851	!!$ & bm_to_litter(:,:,isapbelow,:) + bm_to_litter(:,:,iheartbelow,:) +&
1852	!!$ & bm_to_litter(:,:,iheartabove,:) + bm_to_litter(:,:,iroot,:) + &
1853	!!$ & bm_to_litter(:,:,ifruit,:) + bm_to_litter(:,:,icarbres,:) + bm_to_litter(:,:,ilabile,:)
1854	!!$
1855	!!$ carb_mass_variation(:)=-carb_mass_total(:)
1856	!!$ carb_mass_total(:)=SUM((tot_live_biomass(:,:,icarbon)+tot_litter_carb+tot_soil_carb)*veget_cov_max,dim=2) + &
1857	!!$ & (prod10_total + prod100_total)
1858	!!$ carb_mass_variation(:)=carb_mass_total(:)+carb_mass_variation(:)
1859	! +++++++++++
1860
1861
1862	!! 19. Write history
1863
1864	!! 19.1 Calculate the latest values before writing
1865	lai(:,1) = zero
1866	qm_dia(:,1) = zero
1867	qm_height(:,1) = zero
1868
1869	vcmax_new=zero
1870
1871	DO ipts = 1,npts
1872	DO ivm=2,nvm
1873
1874	lai(ipts,ivm) = cc_to_lai(circ_class_biomass(ipts,ivm,:,ileaf,icarbon),&
1875	circ_class_n(ipts,ivm,:),ivm)
1876	qm_dia(ipts,ivm)=wood_to_qmdia(circ_class_biomass(ipts,ivm,:,:,icarbon), &
1877	circ_class_n(ipts,ivm,:), ivm)
1878	qm_height(ipts,ivm)=wood_to_qmheight(circ_class_biomass(ipts,ivm,:,:,icarbon), &
1879	circ_class_n(ipts,ivm,:), ivm)
1880
1881	IF(SUM(circ_class_biomass(ipts,ivm,:,ileaf,icarbon)*&
1882	circ_class_n(ipts,ivm,:)) .GT. min_stomate)THEN
1883
1884	vcmax_new(ipts,ivm)=nue(ipts,ivm)SUM(circ_class_biomass(ipts,ivm,:,ileaf,initrogen)&
1885	circ_class_n(ipts,ivm,:))* ext_coeff_N(ivm) / &
1886	( 1.-exp(-ext_coeff_N(ivm) * lai(ipts,ivm)))
1887
1888	ENDIF
1889	ENDDO
1890	ENDDO
1891
1892	! Soil, litter, biomass and llc
1893	! Total living biomass (gC m-2)i
1894	DO l=1,nelements
1895	tot_live_biomass(:,:,l) = SUM((circ_class_biomass(:,:,:,ileaf,l) +&
1896	circ_class_biomass(:,:,:,isapabove,l) + &
1897	circ_class_biomass(:,:,:,isapbelow,l) + &
1898	circ_class_biomass(:,:,:,iheartabove,l) + &
1899	circ_class_biomass(:,:,:,iheartbelow,l) + &
1900	circ_class_biomass(:,:,:,iroot,l) + &
1901	circ_class_biomass(:,:,:,ifruit,l) + &
1902	circ_class_biomass(:,:,:,icarbres,l))*circ_class_n(:,:,:),3)
1903	ENDDO
1904
1905	!! Calculate the wood volume change at pixel level by ncuts
1906	DO ivm = 1,nvm
1907	IF(is_tree(ivm))THEN
1908	DO icut= 1,ncut_times
1909
1910	! This variable is used to monitor the change in biomass due
1911	! to different reasons for cutting trees. For this variable we
1912	! don't care about the exact forest management but we do care
1913	! whether the forest was managed or not. If the forest is
1914	! managed the cuttings will be exported. If the forest is not
1915	! managed the cuttings are left on-site.
1916	wood_volume_pix_cut(:,icut) = wood_volume_pix_cut(:,icut) + &
1917	((biomass_cut(:,ivm,isapabove,icarbon,ifm_none,icut) + &
1918	biomass_cut(:,ivm,iheartabove,icarbon,ifm_none,icut)) * &
1919	(1-branch_ratio(ivm)) / pipe_density(ivm) + &
1920	(biomass_cut(:,ivm,isapabove,icarbon,ifm_thin,icut) + &
1921	biomass_cut(:,ivm,iheartabove,icarbon,ifm_thin,icut)) * &
1922	(1-branch_ratio(ivm)) / pipe_density(ivm) + &
1923	(biomass_cut(:,ivm,isapabove,icarbon,ifm_cop,icut) + &
1924	biomass_cut(:,ivm,iheartabove,icarbon,ifm_cop,icut)) * &
1925	(1-branch_ratio(ivm)) / pipe_density(ivm) + &
1926	(biomass_cut(:,ivm,isapabove,icarbon,ifm_src,icut) + &
1927	biomass_cut(:,ivm,iheartabove,icarbon,ifm_src,icut)) * &
1928	(1-branch_ratio(ivm)) / pipe_density(ivm)) * &
1929	veget_cov_max(:,ivm)
1930
1931	ENDDO !icut
1932	ENDIF ! is_tree
1933	ENDDO !ivm
1934
1935	CALL xios_orchidee_send_field("WOOD_VOL_PIX_CUT",wood_volume_pix_cut)
1936
1937	!! 19.2 Write using XIOS
1938	!! 19.2.1 stomate history
1939	CALL xios_orchidee_send_field("RESOLUTION_X",resolution(:,1))
1940	CALL xios_orchidee_send_field("RESOLUTION_Y",resolution(:,2))
1941	CALL xios_orchidee_send_field("VEGET_COV_MAX",veget_cov_max)
1942	CALL xios_orchidee_send_field("CONTFRAC_STOMATE",contfrac(:))
1943
1944	! climatic conditions
1945	CALL xios_orchidee_send_field("T2M_MONTH",t2m_month)
1946	CALL xios_orchidee_send_field("T2M_WEEK",t2m_week)
1947	CALL xios_orchidee_send_field("TSEASON",Tseason)
1948	CALL xios_orchidee_send_field("TMIN_SPRING_TIME",Tmin_spring_time)
1949	CALL xios_orchidee_send_field("FPC_MAX",fpc_max)
1950	CALL xios_orchidee_send_field("MAXFPC_LASTYEAR",maxfpc_lastyear)
1951	!CALL xios_orchidee_send_field("ONSET_DATE",onset_date(:,:,2))
1952	CALL xios_orchidee_send_field("LITTERHUM",litterhum_daily)
1953	CALL xios_orchidee_send_field("FIREINDEX",fireindex(:,:))
1954
1955	! Fluxes
1956	CALL xios_orchidee_send_field("NPP_STOMATE",npp_daily)
1957	CALL xios_orchidee_send_field("GPP",gpp_daily)
1958	CALL xios_orchidee_send_field("HET_RESP",resp_hetero(:,:))
1959	CALL xios_orchidee_send_field("CO2_TAKEN",SUM(atm_to_bm_hist(:,:,:,initrogen),3))
1960	IF(ok_ncycle) CALL xios_orchidee_send_field("N_TAKEN",SUM(atm_to_bm_hist(:,:,:,initrogen),3))
1961	CALL xios_orchidee_send_field("RDI",rdi(:,:))
1962	CALL xios_orchidee_send_field("RDI_TARGET_UPPER",rdi_target_upper(:,:))
1963	CALL xios_orchidee_send_field("RDI_TARGET_LOWER",rdi_target_lower(:,:))
1964	CALL xios_orchidee_send_field("MAINT_RESP",resp_maint)
1965	CALL xios_orchidee_send_field("GROWTH_RESP",resp_growth)
1966	CALL xios_orchidee_send_field("PLANT_STATUS",plant_status)
1967
1968	! Fire results
1969	CALL xios_orchidee_send_field("CO2_FIRE",co2_fire)
1970
1971	! Soil, litter, biomass and llc
1972	! Total living biomass (gC m-2)
1973	DO l=1,nelements
1974	tot_live_biomass(:,:,l) = SUM((circ_class_biomass(:,:,:,ileaf,l) +&
1975	circ_class_biomass(:,:,:,isapabove,l) + &
1976	circ_class_biomass(:,:,:,isapbelow,l) + &
1977	circ_class_biomass(:,:,:,iheartabove,l) + &
1978	circ_class_biomass(:,:,:,iheartbelow,l) + &
1979	circ_class_biomass(:,:,:,iroot,l) + &
1980	circ_class_biomass(:,:,:,ifruit,l) + &
1981	circ_class_biomass(:,:,:,icarbres,l))*circ_class_n(:,:,:),3)
1982	ENDDO
1983
1984	DO l=1,nelements
1985	IF (l == icarbon) THEN
1986	element_str(l) = '_c'
1987	ELSEIF (l == initrogen) THEN
1988	element_str(l) = '_n'
1989	ELSE
1990	STOP 'Define element_str'
1991	ENDIF
1992
1993	!! The writting could be more compact, but writting all the components to
1994	!! seperate variables might make it more userfriendly i.e. the user need not to
1995	!! remember the order of the indexes of nelement and nparts.(asl)
1996	! Biomass
1997	CALL xios_orchidee_send_field("TOTAL_M"//TRIM(element_str(l)),tot_live_biomass(:,:,l))
1998	CALL xios_orchidee_send_field("LEAF_M"//TRIM(element_str(l)),&
1999	SUM(circ_class_biomass(:,:,:,ileaf,l)*circ_class_n(:,:,:),3))
2000	CALL xios_orchidee_send_field("SAP_M_AB"//TRIM(element_str(l)),&
2001	SUM(circ_class_biomass(:,:,:,isapabove,l)*circ_class_n(:,:,:),3))
2002	CALL xios_orchidee_send_field("SAP_M_BE"//TRIM(element_str(l)),&
2003	SUM(circ_class_biomass(:,:,:,isapbelow,l)*circ_class_n(:,:,:),3))
2004	CALL xios_orchidee_send_field("HEART_M_AB"//TRIM(element_str(l)),&
2005	SUM(circ_class_biomass(:,:,:,iheartabove,l)*circ_class_n(:,:,:),3))
2006	CALL xios_orchidee_send_field("HEART_M_BE"//TRIM(element_str(l)),&
2007	SUM(circ_class_biomass(:,:,:,iheartbelow,l)*circ_class_n(:,:,:),3))
2008	CALL xios_orchidee_send_field("ROOT_M"//TRIM(element_str(l)),&
2009	SUM(circ_class_biomass(:,:,:,iroot,l)*circ_class_n(:,:,:),3))
2010	CALL xios_orchidee_send_field("FRUIT_M"//TRIM(element_str(l)),&
2011	SUM(circ_class_biomass(:,:,:,ifruit,l)*circ_class_n(:,:,:),3))
2012	CALL xios_orchidee_send_field("LABILE_M"//TRIM(element_str(l)),&
2013	SUM(circ_class_biomass(:,:,:,ilabile,l)*circ_class_n(:,:,:),3))
2014	CALL xios_orchidee_send_field("RESERVE_M"//TRIM(element_str(l)),&
2015	SUM(circ_class_biomass(:,:,:,icarbres,l)*circ_class_n(:,:,:),3))
2016
2017	!CALL xios_orchidee_send_field("TOTAL_TURN"//TRIM(element_str(l)),tot_turnover(:,:,l))
2018	CALL xios_orchidee_send_field("LEAF_TURN"//TRIM(element_str(l)),turnover_daily(:,:,ileaf,l))
2019	CALL xios_orchidee_send_field("SAP_AB_TURN"//TRIM(element_str(l)),turnover_daily(:,:,isapabove,l))
2020	CALL xios_orchidee_send_field("ROOT_TURN"//TRIM(element_str(l)),turnover_daily(:,:,iroot,l))
2021	CALL xios_orchidee_send_field("FRUIT_TURN"//TRIM(element_str(l)),turnover_daily(:,:,ifruit,l))
2022
2023	!CALL xios_orchidee_send_field("TOTAL_BM_LITTER"//TRIM(element_str(l)),tot_bm_to_litter(:,:,l))
2024	CALL xios_orchidee_send_field("LEAF_BM_LITTER"//TRIM(element_str(l)),bm_to_litter(:,:,ileaf,l))
2025	CALL xios_orchidee_send_field("SAP_AB_BM_LITTER"//TRIM(element_str(l)),bm_to_litter(:,:,isapabove,l))
2026	CALL xios_orchidee_send_field("SAP_BE_BM_LITTER"//TRIM(element_str(l)),bm_to_litter(:,:,isapbelow,l))
2027	CALL xios_orchidee_send_field("HEART_AB_BM_LITTER"//TRIM(element_str(l)),bm_to_litter(:,:,iheartabove,l))
2028	CALL xios_orchidee_send_field("HEART_BE_BM_LITTER"//TRIM(element_str(l)),bm_to_litter(:,:,iheartbelow,l))
2029	CALL xios_orchidee_send_field("ROOT_BM_LITTER"//TRIM(element_str(l)),bm_to_litter(:,:,iroot,l))
2030	CALL xios_orchidee_send_field("FRUIT_BM_LITTER"//TRIM(element_str(l)),bm_to_litter(:,:,ifruit,l))
2031	CALL xios_orchidee_send_field("LABILE_BM_LITTER"//TRIM(element_str(l)),bm_to_litter(:,:,ilabile,l))
2032	CALL xios_orchidee_send_field("RESERVE_BM_LITTER"//TRIM(element_str(l)),bm_to_litter(:,:,icarbres,l))
2033
2034	! Litter
2035	CALL xios_orchidee_send_field("LITTER_STR_AB"//TRIM(element_str(l)),litter(:,istructural,:,iabove,l))
2036	CALL xios_orchidee_send_field("LITTER_MET_AB"//TRIM(element_str(l)),litter(:,imetabolic,:,iabove,l))
2037	CALL xios_orchidee_send_field("LITTER_WOD_AB"//TRIM(element_str(l)),litter(:,iwoody,:,iabove,l))
2038	CALL xios_orchidee_send_field("LITTER_STR_BE"//TRIM(element_str(l)),litter(:,istructural,:,ibelow,l))
2039	CALL xios_orchidee_send_field("LITTER_MET_BE"//TRIM(element_str(l)),litter(:,imetabolic,:,ibelow,l))
2040	CALL xios_orchidee_send_field("LITTER_WOD_BE"//TRIM(element_str(l)),litter(:,iwoody,:,ibelow,l))
2041
2042	! Soil
2043	CALL xios_orchidee_send_field("SOIL_ACTIVE"//TRIM(element_str(l)),som(:,iactive,:,l))
2044	CALL xios_orchidee_send_field("SOIL_SLOW"//TRIM(element_str(l)),som(:,islow,:,l))
2045	CALL xios_orchidee_send_field("SOIL_PASSIVE"//TRIM(element_str(l)),som(:,ipassive,:,l))
2046	CALL xios_orchidee_send_field("SOIL_SURF"//TRIM(element_str(l)),som(:,isurface,:,l))
2047
2048	! land cover change
2049	CALL xios_orchidee_send_field('FLUX_PROD_S'//TRIM(element_str(l)),flux_prod_s(:,l))
2050	CALL xios_orchidee_send_field('FLUX_PROD_M'//TRIM(element_str(l)),flux_prod_m(:,l))
2051	CALL xios_orchidee_send_field('FLUX_PROD_L'//TRIM(element_str(l)),flux_prod_l(:,l))
2052	CALL xios_orchidee_send_field('PROD_S'//TRIM(element_str(l)),prod_s(:,:,l))
2053	CALL xios_orchidee_send_field('PROD_M'//TRIM(element_str(l)),prod_m(:,:,l))
2054	CALL xios_orchidee_send_field('PROD_L'//TRIM(element_str(l)),prod_l(:,:,l))
2055	CALL xios_orchidee_send_field('FLUX_S'//TRIM(element_str(l)),flux_s(:,:,l))
2056	CALL xios_orchidee_send_field('FLUX_M'//TRIM(element_str(l)),flux_m(:,:,l))
2057	CALL xios_orchidee_send_field('FLUX_L'//TRIM(element_str(l)),flux_l(:,:,l))
2058
2059	ENDDO
2060
2061	!CN longterm ratios
2062	IF (spinup_analytic) THEN
2063	CALL xios_orchidee_send_field("CN_LONGTERM_ACTIVE",CN_som_litter_longterm(:,:,iactive_pool))
2064	CALL xios_orchidee_send_field("CN_LONGTERM_SLOW",CN_som_litter_longterm(:,:,islow_pool))
2065	CALL xios_orchidee_send_field("CN_LONGTERM_PASSIVE",CN_som_litter_longterm(:,:,ipassive_pool))
2066	CALL xios_orchidee_send_field("CN_LONGTERM_SURF",CN_som_litter_longterm(:,:,isurface_pool))
2067	CALL xios_orchidee_send_field("CN_LONGTERM_LITTER_STR_AB",CN_som_litter_longterm(:,:,istructural_above))
2068	CALL xios_orchidee_send_field("CN_LONGTERM_LITTER_STR_BE",CN_som_litter_longterm(:,:,istructural_below))
2069	CALL xios_orchidee_send_field("CN_LONGTERM_LITTER_MET_AB",CN_som_litter_longterm(:,:,imetabolic_above))
2070	CALL xios_orchidee_send_field("CN_LONGTERM_LITTER_MET_BE",CN_som_litter_longterm(:,:,imetabolic_below))
2071	CALL xios_orchidee_send_field("CN_LONGTERM_LITTER_WOD_AB",CN_som_litter_longterm(:,:,iwoody_above))
2072	CALL xios_orchidee_send_field("CN_LONGTERM_LITTER_WOD_BE",CN_som_litter_longterm(:,:,iwoody_below))
2073	ENDIF
2074
2075	! Stand structure
2076	CALL xios_orchidee_send_field("DIAMETER",qm_dia)
2077	CALL xios_orchidee_send_field("LAI_PER_LEVEL",lai_per_level)
2078	CALL xios_orchidee_send_field("CLEVEL_HEIGHT",z_level_photo)
2079	CALL xios_orchidee_send_field("lai_ipcc",SUM(laiveget_cov_max,dim=2)contfrac)
2080	CALL xios_orchidee_send_field("LAI",lai) ! this is written twice!
2081	CALL xios_orchidee_send_field("HEIGHT",qm_height)
2082	CALL xios_orchidee_send_field("IND",SUM(circ_class_n(:,:,:),3))
2083	CALL xios_orchidee_send_field("WOODMASS_IND",woodmass_ind)
2084	CALL xios_orchidee_send_field("AGE",age)
2085	CALL xios_orchidee_send_field("DEADLEAF_COVER",deadleaf_cover)
2086	CALL xios_orchidee_send_field("TOTAL_SOIL_CARB",tot_litter_soil_carb)
2087	CALL xios_orchidee_send_field("CN_IND",cn_ind)
2088
2089	! Growth stress and factors
2090	CALL xios_orchidee_send_field("MOISTRESS",moiavail_week)
2091	CALL xios_orchidee_send_field("VCMAX",vcmax)
2092	CALL xios_orchidee_send_field("VCMAX_NEW",vcmax_new)
2093	CALL xios_orchidee_send_field("TURNOVER_TIME_LEAF",turnover_time(:,:,ileaf))
2094	CALL xios_orchidee_send_field("TURNOVER_TIME_ROOT",turnover_time(:,:,iroot))
2095	CALL xios_orchidee_send_field("TURNOVER_TIME_SAP_AB",turnover_time(:,:,isapabove))
2096	CALL xios_orchidee_send_field("TURNOVER_TIME_FRUIT",turnover_time(:,:,ifruit))
2097
2098
2099	! for functional allocation, we have all of this information
2100	! at every step so we can write it out.
2101	var_real=REAL(age_stand)
2102	CALL xios_orchidee_send_field("AGE_STAND",var_real)
2103	var_real=REAL(forest_managed)
2104	CALL xios_orchidee_send_field("FOREST_MANAGED",var_real)
2105	var_real=REAL(last_cut)
2106	CALL xios_orchidee_send_field("LAST_CUT",var_real)
2107	var_real=REAL(rotation_n)
2108	CALL xios_orchidee_send_field("ROTATION_N",var_real)
2109	CALL xios_orchidee_send_field("LITTER_RAKE_FRAC",lrake_frac)
2110	CALL xios_orchidee_send_field("MAI",mai)
2111	CALL xios_orchidee_send_field("PAI",pai)
2112	CALL xios_orchidee_send_field("SIGMA",sigma)
2113	CALL xios_orchidee_send_field("GAMMA",gammas)
2114
2115
2116	!! 19.2.2 ipcc history
2117	CALL xios_orchidee_send_field('RESOLUTION_X',resolution(:,1))
2118	CALL xios_orchidee_send_field('RESOLUTION_Y',resolution(:,2))
2119	CALL xios_orchidee_send_field('CONTFRAC_STOMATE',contfrac(:))
2120
2121	!CALL xios_orchidee_send_field("cVeg",SUM(tot_live_biomass(:,:,icarbon)*veget_cov_max,dim=2)/1e3)
2122	!CALL xios_orchidee_send_field("cLitter",SUM(tot_litter_carb*veget_cov_max,dim=2)/1e3)
2123	!CALL xios_orchidee_send_field("cSoil",SUM(tot_soil_carb*veget_cov_max,dim=2)/1e3)
2124	!!$ CALL xios_orchidee_send_field("cProduct",(prod10_total + prod100_total)/1e3)
2125	! CALL xios_orchidee_send_field("cMassVariation",carb_mass_variation/1e3/one_day)
2126
2127	CALL xios_orchidee_send_field("lai_ipcc",SUM(lai*veget_cov_max,dim=2)) ! m2/m2
2128
2129	! Carbon fluxes transformed from gC/m2/d into kgC/m2/s
2130	CALL xios_orchidee_send_field("gpp_ipcc",SUM(gpp_daily*veget_cov_max,dim=2)/1e3/one_day)
2131	CALL xios_orchidee_send_field("ra",SUM((resp_maint+resp_growth)*veget_cov_max,dim=2)/1e3/one_day)
2132	CALL xios_orchidee_send_field("npp_ipcc",SUM(npp_daily*veget_cov_max,dim=2)/1e3/one_day)
2133	CALL xios_orchidee_send_field("rh",SUM(resp_hetero*veget_cov_max,dim=2)/1e3/one_day)
2134	CALL xios_orchidee_send_field("fFire",SUM(co2_fire*veget_cov_max,dim=2)/1e3/one_day)
2135	! As LLC is not acitivated yet, flux_prod_total is not calculated.
2136	! Likewise, harvest_above is no longer used. Possibly, it can be replaced
2137	! by harvest_pool? If yes, be careful with the units, and should be done further
2138	! up the code, as harvest_pool is cleared in basic_product_use.
2139	! CALL xios_orchidee_send_field("fHarvest",harvest_above/1e3/one_day)
2140	! CALL xios_orchidee_send_field("fLuc",flux_prod_total/1e3/one_day)
2141	! CALL xios_orchidee_send_field("nbp",(SUM((gpp_daily + atm_to_bm(:,:,icarbon) - &
2142	! (resp_maint+resp_growth+resp_hetero)-co2_fire) *veget_cov_max,dim=2)- &
2143	! flux_prod_total-harvest_above)/1e3/one_day)
2144	!CALL xios_orchidee_send_field("fVegLitter",SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))*&
2145	! veget_cov_max,dim=2)/1e3/one_day)
2146	CALL xios_orchidee_send_field("fLitterSoil",SUM(SUM(som_input(:,:,:,icarbon),dim=2)*veget_cov_max,dim=2)/1e3/one_day)
2147	CALL xios_orchidee_send_field("cLeaf",SUM(SUM(circ_class_biomass(:,:,:,ileaf,icarbon)*&
2148	circ_class_n(:,:,:),3)*veget_cov_max,dim=2)/1e3)
2149	CALL xios_orchidee_send_field("cWood",SUM(SUM((circ_class_biomass(:,:,:,isapabove,icarbon)+&
2150	circ_class_biomass(:,:,:,iheartabove,icarbon))circ_class_n(:,:,:),3)veget_cov_max,dim=2)/1e3)
2151	CALL xios_orchidee_send_field("cRoot",SUM(SUM((circ_class_biomass(:,:,:,iroot,icarbon) + &
2152	circ_class_biomass(:,:,:,isapbelow,icarbon) + &
2153	circ_class_biomass(:,:,:,iheartbelow,icarbon) )circ_class_n(:,:,:),3)veget_cov_max,dim=2)/1e3)
2154	CALL xios_orchidee_send_field("cMisc",SUM(SUM((circ_class_biomass(:,:,:,icarbres,icarbon) + &
2155	circ_class_biomass(:,:,:,ifruit,icarbon))circ_class_n(:,:,:),3)veget_cov_max,dim=2)/1e3)
2156	CALL xios_orchidee_send_field("cLitterAbove",SUM((litter(:,istructural,:,iabove,icarbon)+&
2157	litter(:,imetabolic,:,iabove,icarbon)+ litter(:,iwoody,:,iabove,icarbon))*veget_cov_max,dim=2)/1e3)
2158	CALL xios_orchidee_send_field("cLitterBelow",SUM((litter(:,istructural,:,ibelow,icarbon)+&
2159	litter(:,imetabolic,:,ibelow,icarbon)+ litter(:,iwoody,:,ibelow,icarbon))*veget_cov_max,dim=2)/1e3)
2160	CALL xios_orchidee_send_field("cSoilFast",SUM((som(:,iactive,:,icarbon)+som(:,isurface,:,icarbon)) * &
2161	veget_cov_max,dim=2)/1e3)
2162	CALL xios_orchidee_send_field("cSoilMedium",SUM(som(:,islow,:,icarbon)*veget_cov_max,dim=2)/1e3)
2163	CALL xios_orchidee_send_field("cSoilSlow",SUM(som(:,ipassive,:,icarbon)*veget_cov_max,dim=2)/1e3)
2164
2165	! Vegetation fractions [0,100]
2166	CALL xios_orchidee_send_field("landCoverFrac",veget_cov_max*100)
2167	vartmp(:)=zero
2168	DO j = 2,nvm
2169	IF (is_deciduous(j)) THEN
2170	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2171	ENDIF
2172	ENDDO
2173
2174	CALL xios_orchidee_send_field("treeFracPrimDec",vartmp)
2175	vartmp(:)=zero
2176	DO j = 2,nvm
2177	IF (is_evergreen(j)) THEN
2178	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2179	ENDIF
2180	ENDDO
2181
2182	CALL xios_orchidee_send_field("treeFracPrimEver",vartmp)
2183	vartmp(:)=zero
2184	DO j = 2,nvm
2185	IF ( .NOT.(is_c4(j)) ) THEN
2186	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2187	ENDIF
2188	ENDDO
2189
2190	CALL xios_orchidee_send_field("c3PftFrac",vartmp)
2191	vartmp(:)=zero
2192	DO j = 2,nvm
2193	IF ( is_c4(j) ) THEN
2194	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2195	ENDIF
2196	ENDDO
2197
2198	CALL xios_orchidee_send_field("c4PftFrac",vartmp)
2199
2200	! Carbon fluxes transformed from gC/m2/d into kgC/m2/s
2201	CALL xios_orchidee_send_field("rGrowth",SUM(resp_growth*veget_cov_max,dim=2)/1e3/one_day)
2202	CALL xios_orchidee_send_field("rMaint",SUM(resp_maint*veget_cov_max,dim=2)/1e3/one_day)
2203	CALL xios_orchidee_send_field("nppLeaf",SUM(bm_alloc(:,:,ileaf,icarbon)*&
2204	veget_cov_max,dim=2)/1e3/one_day)
2205	CALL xios_orchidee_send_field("nppWood",SUM(bm_alloc(:,:,isapabove,icarbon)*&
2206	veget_cov_max,dim=2)/1e3/one_day)
2207	CALL xios_orchidee_send_field("nppRoot",SUM(( bm_alloc(:,:,isapbelow,icarbon) + &
2208	bm_alloc(:,:,iroot,icarbon) )*veget_cov_max,dim=2)/1e3/one_day)
2209
2210	! 19.3 Write to file using histwrite
2211	CALL histwrite_p (hist_id_stomate, 'RESOLUTION_X', itime, &
2212	resolution(:,1), npts, hori_index)
2213	CALL histwrite_p (hist_id_stomate, 'RESOLUTION_Y', itime, &
2214	resolution(:,2), npts, hori_index)
2215	CALL histwrite_p (hist_id_stomate, 'VEGET_COV_MAX', itime, &
2216	veget_cov_max, npts*nvm, horipft_index)
2217	CALL histwrite_p (hist_id_stomate, 'CONTFRAC', itime, &
2218	contfrac(:), npts, hori_index)
2219
2220	! Climatic conditions
2221	CALL histwrite_p (hist_id_stomate, 'T2M_MONTH', itime, &
2222	t2m_month, npts, hori_index)
2223	CALL histwrite_p (hist_id_stomate, 'T2M_WEEK', itime, &
2224	t2m_week, npts, hori_index)
2225	CALL histwrite_p (hist_id_stomate, 'TSEASON', itime, &
2226	Tseason, npts, hori_index)
2227	CALL histwrite_p (hist_id_stomate, 'TMIN_SPRING_TIME', itime, &
2228	Tmin_spring_time, npts*nvm, horipft_index)
2229	! CALL histwrite_p (hist_id_stomate, 'ONSET_DATE', itime, &
2230	! onset_date(:,:,2), npts*nvm, horipft_index)
2231	CALL histwrite_p (hist_id_stomate, 'FPC_MAX', itime, &
2232	fpc_max, npts*nvm, horipft_index)
2233	CALL histwrite_p (hist_id_stomate, 'MAXFPC_LASTYEAR', itime, &
2234	maxfpc_lastyear, npts*nvm, horipft_index)
2235	CALL histwrite_p (hist_id_stomate, 'FIREINDEX', itime, &
2236	fireindex(:,:), npts*nvm, horipft_index)
2237	CALL histwrite_p (hist_id_stomate, 'LITTERHUM', itime, &
2238	litterhum_daily, npts, hori_index)
2239
2240	! Fire results
2241	CALL histwrite_p (hist_id_stomate, 'CO2_FIRE', itime, &
2242	co2_fire, npts*nvm, horipft_index)
2243
2244	! Stand structure
2245	CALL histwrite_p (hist_id_stomate, 'LAI', itime, &
2246	lai, npts*nvm, horipft_index)
2247	CALL histwrite_p (hist_id_stomate, 'DIAMETER', itime, &
2248	qm_dia, npts*nvm, horipft_index)
2249	CALL histwrite_p (hist_id_stomate, 'HEIGHT', itime, &
2250	qm_height, npts*nvm, horipft_index)
2251	CALL histwrite_p (hist_id_stomate, 'IND', itime, &
2252	SUM(circ_class_n(:,:,:),3), npts*nvm, horipft_index)
2253
2254	! Calculate lai per level. Also print out the real height of
2255	! the canopy level.
2256	DO ivm=1,nvm
2257
2258	WRITE(var_name,'(A,I3.3)') 'LAI_PER_LEVEL_',ivm
2259	CALL histwrite_p (hist_id_stomate, var_name, itime, &
2260	lai_per_level(:,ivm,:), npts*nlevels_tot, horican_index)
2261
2262	WRITE(var_name,'(A,I3.3)') 'CLEVEL_HEIGHT_',ivm
2263	CALL histwrite_p (hist_id_stomate, var_name, itime, &
2264	z_level_photo(:,ivm,:), npts*nlevels_tot, horican_index)
2265
2266	ENDDO
2267
2268	!!$ CALL histwrite_p (hist_id_stomate, 'CN_IND', itime, &
2269	!!$ cn_ind, npts*nvm, horipft_index)
2270	CALL histwrite_p (hist_id_stomate, 'WOODMASS_IND', itime, &
2271	woodmass_ind, npts*nvm, horipft_index)
2272	CALL histwrite_p (hist_id_stomate, 'AGE', itime, &
2273	age, npts*nvm, horipft_index)
2274
2275	! Growth stress and factors
2276	CALL histwrite_p (hist_id_stomate, 'MOISTRESS', itime, &
2277	moiavail_week, npts*nvm, horipft_index)
2278	CALL histwrite_p (hist_id_stomate, 'VCMAX', itime, &
2279	vcmax, npts*nvm, horipft_index)
2280	CALL histwrite_p (hist_id_stomate, 'VCMAX_NEW', itime, &
2281	vcmax_new, npts*nvm, horipft_index)
2282
2283	!+++CHECK+++
2284	! New dimensions
2285	!! CALL histwrite_p (hist_id_stomate, 'TURNOVER_TIME_LEAF', itime, &
2286	!! turnover_time(:,:,ileaf), npts*nvm, horipft_index)
2287	!! CALL histwrite_p (hist_id_stomate, 'TURNOVER_TIME_ROOT', itime, &
2288	!! turnover_time(:,:,iroot), npts*nvm, horipft_index)
2289	!! CALL histwrite_p (hist_id_stomate, 'TURNOVER_TIME_SAP', itime, &
2290	!! turnover_time(:,:,isapabove), npts*nvm, horipft_index)
2291	!! CALL histwrite_p (hist_id_stomate, 'TURNOVER_TIME_FRUIT', itime, &
2292	!! turnover_time(:,:,ifruit), npts*nvm, horipft_index)
2293	!+++++++++++
2294
2295	! for functional allocation, we have all of this information
2296	! at every step so we can write it out.
2297	var_real=REAL(age_stand)
2298	CALL histwrite_p (hist_id_stomate, 'AGE_STAND', itime, &
2299	var_real, npts*nvm, horipft_index)
2300	CALL histwrite_p (hist_id_stomate, 'MAI', itime, &
2301	mai, npts*nvm, horipft_index)
2302	CALL histwrite_p (hist_id_stomate, 'PAI', itime, &
2303	pai, npts*nvm, horipft_index)
2304	var_real=REAL(rotation_n)
2305	CALL histwrite_p (hist_id_stomate, 'ROTATION_N', itime, &
2306	var_real, npts*nvm, horipft_index)
2307	var_real=REAL(last_cut)
2308	CALL histwrite_p (hist_id_stomate, 'LAST_CUT', itime, &
2309	var_real, npts*nvm, horipft_index)
2310	var_real=REAL(forest_managed)
2311	CALL histwrite_p (hist_id_stomate, 'FOREST_MANAGED', itime, &
2312	var_real, npts*nvm, horipft_index)
2313	CALL histwrite_p (hist_id_stomate, 'LITTER_RAKE_FRAC', itime, &
2314	lrake_frac, npts*nvm, horipft_index)
2315	CALL histwrite_p (hist_id_stomate, 'SIGMA', itime, &
2316	sigma, npts*nvm, horipft_index)
2317	CALL histwrite_p (hist_id_stomate, 'GAMMA', itime, &
2318	gammas, npts*nvm, horipft_index)
2319
2320
2321	! Fluxes
2322	CALL histwrite_p (hist_id_stomate, 'NPP', itime, &
2323	npp_daily, npts*nvm, horipft_index)
2324	CALL histwrite_p (hist_id_stomate, 'GPP', itime, &
2325	gpp_daily, npts*nvm, horipft_index)
2326	CALL histwrite_p (hist_id_stomate, 'HET_RESP', itime, &
2327	resp_hetero(:,:), npts*nvm, horipft_index)
2328	CALL histwrite_p (hist_id_stomate, 'CO2_TAKEN', itime, &
2329	SUM(atm_to_bm_hist(:,:,:,icarbon),3), npts*nvm, horipft_index)
2330	IF(ok_ncycle) CALL histwrite_p (hist_id_stomate, 'N_TAKEN', itime, &
2331	SUM(atm_to_bm_hist(:,:,:,initrogen),3), npts*nvm, horipft_index)
2332	CALL histwrite_p (hist_id_stomate, 'RDI', itime, &
2333	rdi, npts*nvm, horipft_index)
2334	CALL histwrite_p (hist_id_stomate, 'RDI_TARGET_UPPER', itime, &
2335	rdi_target_upper, npts*nvm, horipft_index)
2336	CALL histwrite_p (hist_id_stomate, 'RDI_TARGET_LOWER', itime, &
2337	rdi_target_lower, npts*nvm, horipft_index)
2338	CALL histwrite_p (hist_id_stomate, 'MAINT_RESP', itime, &
2339	resp_maint, npts*nvm, horipft_index)
2340	CALL histwrite_p (hist_id_stomate, 'GROWTH_RESP', itime, &
2341	resp_growth, npts*nvm, horipft_index)
2342	CALL histwrite_p (hist_id_stomate, 'PLANT_STATUS', itime, &
2343	plant_status, npts*nvm, horipft_index)
2344
2345	! Soil, litter and biomass
2346	CALL histwrite_p (hist_id_stomate, 'DEADLEAF_COVER', itime, &
2347	deadleaf_cover, npts, hori_index)
2348	! CALL histwrite_p (hist_id_stomate, 'TOTAL_SOIL_CARB', itime, &
2349	! tot_litter_soil_carb, npts*nvm, horipft_index)
2350	! +++CHECK+++
2351	! If needed declare in ioipscontrol first
2352	!!$ CALL histwrite_p (hist_id_stomate, 'LIGNIN_STRUC_AB', itime, &
2353	!!$ lignin_struc(:,:,iabove), npts*nvm, horipft_index)
2354	!!$ CALL histwrite_p (hist_id_stomate, 'LIGNIN_STRUC_BE', itime, &
2355	!!$ lignin_struc(:,:,ibelow), npts*nvm, horipft_index)
2356	!!$ CALL histwrite_p (hist_id_stomate, 'LIGNIN_WOOD_AB', itime, &
2357	!!$ lignin_wood(:,:,iabove), npts*nvm, horipft_index)
2358	!!$ CALL histwrite_p (hist_id_stomate, 'LIGNIN_WOOD_BE', itime, &
2359	!!$ lignin_wood(:,:,ibelow), npts*nvm, horipft_index)
2360	! +++++++++++
2361
2362	!!$ WRITE(numout,*) 'Back in stomate_lpj - 4'
2363
2364	DO l=1,nelements
2365	IF (l == icarbon) THEN
2366	element_str(l) = '_c'
2367	ELSEIF (l == initrogen) THEN
2368	element_str(l) = '_n'
2369	ELSE
2370	STOP 'Define element_str'
2371	ENDIF
2372
2373	! Soil and litter per element
2374	CALL histwrite_p (hist_id_stomate, 'LITTER_STR_AB'//TRIM(element_str(l)), itime, &
2375	litter(:,istructural,:,iabove,l), npts*nvm, horipft_index)
2376	CALL histwrite_p (hist_id_stomate, 'LITTER_MET_AB'//TRIM(element_str(l)), itime, &
2377	litter(:,imetabolic,:,iabove,l), npts*nvm, horipft_index)
2378	CALL histwrite_p (hist_id_stomate, 'LITTER_STR_BE'//TRIM(element_str(l)), itime, &
2379	litter(:,istructural,:,ibelow,l), npts*nvm, horipft_index)
2380	CALL histwrite_p (hist_id_stomate, 'LITTER_MET_BE'//TRIM(element_str(l)), itime, &
2381	litter(:,imetabolic,:,ibelow,l), npts*nvm, horipft_index)
2382	CALL histwrite_p (hist_id_stomate, 'LITTER_WOD_AB'//TRIM(element_str(l)), itime, &
2383	litter(:,iwoody,:,iabove,l), npts*nvm, horipft_index)
2384	CALL histwrite_p (hist_id_stomate, 'LITTER_WOD_BE'//TRIM(element_str(l)), itime, &
2385	litter(:,iwoody,:,ibelow,l), npts*nvm, horipft_index)
2386	CALL histwrite_p (hist_id_stomate, 'SOIL_ACTIVE'//TRIM(element_str(l)), itime, &
2387	som(:,iactive,:,l), npts*nvm, horipft_index)
2388	CALL histwrite_p (hist_id_stomate, 'SOIL_SLOW'//TRIM(element_str(l)), itime, &
2389	som(:,islow,:,l), npts*nvm, horipft_index)
2390	CALL histwrite_p (hist_id_stomate, 'SOIL_PASSIVE'//TRIM(element_str(l)), itime, &
2391	som(:,ipassive,:,l), npts*nvm, horipft_index)
2392	CALL histwrite_p (hist_id_stomate, 'SOIL_SURF'//TRIM(element_str(l)), itime, &
2393	som(:,isurface,:,l), npts*nvm, horipft_index)
2394
2395	! Biomass per element (g C m-2)
2396	CALL histwrite_p (hist_id_stomate, 'TOTAL_M'//TRIM(element_str(l)), itime, &
2397	tot_live_biomass(:,:,l), npts*nvm, horipft_index)
2398	CALL histwrite_p (hist_id_stomate, 'LEAF_M'//TRIM(element_str(l)), itime, &
2399	SUM(circ_class_biomass(:,:,:,ileaf,l)*circ_class_n(:,:,:),3),&
2400	npts*nvm, horipft_index)
2401	CALL histwrite_p (hist_id_stomate, 'SAP_M_AB'//TRIM(element_str(l)), itime, &
2402	SUM(circ_class_biomass(:,:,:,isapabove,l)circ_class_n(:,:,:),3), nptsnvm, horipft_index)
2403	CALL histwrite_p (hist_id_stomate, 'SAP_M_BE'//TRIM(element_str(l)), itime, &
2404	SUM(circ_class_biomass(:,:,:,isapbelow,l)circ_class_n(:,:,:),3), nptsnvm, horipft_index)
2405	CALL histwrite_p (hist_id_stomate, 'HEART_M_AB'//TRIM(element_str(l)), itime, &
2406	SUM(circ_class_biomass(:,:,:,iheartabove,l)circ_class_n(:,:,:),3), nptsnvm, horipft_index)
2407	CALL histwrite_p (hist_id_stomate, 'HEART_M_BE'//TRIM(element_str(l)), itime, &
2408	SUM(circ_class_biomass(:,:,:,iheartbelow,l)circ_class_n(:,:,:),3), nptsnvm, horipft_index)
2409	CALL histwrite_p (hist_id_stomate, 'ROOT_M'//TRIM(element_str(l)), itime, &
2410	SUM(circ_class_biomass(:,:,:,iroot,l)circ_class_n(:,:,:),3), nptsnvm, horipft_index)
2411	CALL histwrite_p (hist_id_stomate, 'FRUIT_M'//TRIM(element_str(l)), itime, &
2412	SUM(circ_class_biomass(:,:,:,ifruit,l)circ_class_n(:,:,:),3), nptsnvm, horipft_index)
2413	CALL histwrite_p (hist_id_stomate, 'RESERVE_M'//TRIM(element_str(l)), itime, &
2414	SUM(circ_class_biomass(:,:,:,icarbres,l)circ_class_n(:,:,:),3), nptsnvm, horipft_index)
2415	CALL histwrite_p (hist_id_stomate, 'LABILE_M'//TRIM(element_str(l)), itime, &
2416	SUM(circ_class_biomass(:,:,:,ilabile,l)circ_class_n(:,:,:),3), nptsnvm, horipft_index)
2417	!CALL histwrite_p (hist_id_stomate, 'TOTAL_TURN'//TRIM(element_str(l)), itime, &
2418	! tot_turnover(:,:,l), npts*nvm, horipft_index)
2419	CALL histwrite_p (hist_id_stomate, 'LEAF_TURN'//TRIM(element_str(l)), itime, &
2420	turnover_daily(:,:,ileaf,l), npts*nvm, horipft_index)
2421	CALL histwrite_p (hist_id_stomate, 'SAP_AB_TURN'//TRIM(element_str(l)), itime, &
2422	turnover_daily(:,:,isapabove,l), npts*nvm, horipft_index)
2423	CALL histwrite_p (hist_id_stomate, 'ROOT_TURN'//TRIM(element_str(l)), itime, &
2424	turnover_daily(:,:,iroot,l), npts*nvm, horipft_index)
2425	CALL histwrite_p (hist_id_stomate, 'FRUIT_TURN'//TRIM(element_str(l)), itime, &
2426	turnover_daily(:,:,ifruit,l), npts*nvm, horipft_index)
2427	!CALL histwrite_p (hist_id_stomate, 'TOTAL_BM_LITTER'//TRIM(element_str(l)), itime, &
2428	! tot_bm_to_litter(:,:,l), npts*nvm, horipft_index)
2429	CALL histwrite_p (hist_id_stomate, 'LEAF_BM_LITTER'//TRIM(element_str(l)), itime, &
2430	bm_to_litter(:,:,ileaf,l), npts*nvm, horipft_index)
2431	CALL histwrite_p (hist_id_stomate, 'SAP_AB_BM_LITTER'//TRIM(element_str(l)), itime, &
2432	bm_to_litter(:,:,isapabove,l), npts*nvm, horipft_index)
2433	CALL histwrite_p (hist_id_stomate, 'SAP_BE_BM_LITTER'//TRIM(element_str(l)), itime, &
2434	bm_to_litter(:,:,isapbelow,l), npts*nvm, horipft_index)
2435	CALL histwrite_p (hist_id_stomate, 'HEART_AB_BM_LITTER'//TRIM(element_str(l)), itime, &
2436	bm_to_litter(:,:,iheartabove,l), npts*nvm, horipft_index)
2437	CALL histwrite_p (hist_id_stomate, 'HEART_BE_BM_LITTER'//TRIM(element_str(l)), itime, &
2438	bm_to_litter(:,:,iheartbelow,l), npts*nvm, horipft_index)
2439	CALL histwrite_p (hist_id_stomate, 'ROOT_BM_LITTER'//TRIM(element_str(l)), itime, &
2440	bm_to_litter(:,:,iroot,l), npts*nvm, horipft_index)
2441	CALL histwrite_p (hist_id_stomate, 'FRUIT_BM_LITTER'//TRIM(element_str(l)), itime, &
2442	bm_to_litter(:,:,ifruit,l), npts*nvm, horipft_index)
2443	CALL histwrite_p (hist_id_stomate, 'RESERVE_BM_LITTER'//TRIM(element_str(l)), itime, &
2444	bm_to_litter(:,:,icarbres,l), npts*nvm, horipft_index)
2445	CALL histwrite_p (hist_id_stomate, 'LABILE_BM_LITTER'//TRIM(element_str(l)), itime, &
2446	bm_to_litter(:,:,ilabile,l), npts*nvm, horipft_index)
2447
2448	! land cover change
2449	CALL histwrite_p (hist_id_stomate, 'FLUX_PROD_S'//TRIM(element_str(l)), itime, &
2450	flux_prod_s(:,l), npts, hori_index)
2451	CALL histwrite_p (hist_id_stomate, 'FLUX_PROD_M'//TRIM(element_str(l)), itime, &
2452	flux_prod_m(:,l), npts, hori_index)
2453	CALL histwrite_p (hist_id_stomate, 'FLUX_PROD_L'//TRIM(element_str(l)), itime, &
2454	flux_prod_l(:,l), npts, hori_index)
2455	CALL histwrite_p (hist_id_stomate, 'PROD_S'//TRIM(element_str(l)), itime, &
2456	prod_s(:,:,l), npts*(nshort+1), horip_ss_index)
2457	CALL histwrite_p (hist_id_stomate, 'PROD_M'//TRIM(element_str(l)), itime, &
2458	prod_m(:,:,l), npts*(nmedium+1), horip_mm_index)
2459	CALL histwrite_p (hist_id_stomate, 'PROD_L'//TRIM(element_str(l)), itime, &
2460	prod_l(:,:,l), npts*(nlong+1), horip_ll_index)
2461	CALL histwrite_p (hist_id_stomate, 'FLUX_S'//TRIM(element_str(l)), itime, &
2462	flux_s(:,:,l), npts*nshort, horip_s_index)
2463	CALL histwrite_p (hist_id_stomate, 'FLUX_M'//TRIM(element_str(l)), itime, &
2464	flux_m(:,:,l), npts*nmedium, horip_m_index)
2465	CALL histwrite_p (hist_id_stomate, 'FLUX_L'//TRIM(element_str(l)), itime, &
2466	flux_l(:,:,l), npts*nlong, horip_l_index)
2467
2468	ENDDO
2469
2470	! Wind throw
2471	CALL histwrite_p (hist_id_stomate, 'WOOD_VPCUT', itime, &
2472	wood_volume_pix_cut, npts*ncut_times, horicut_index)
2473
2474	! IPCC
2475	IF ( hist_id_stomate_IPCC > 0 ) THEN
2476	vartmp(:)=SUM(tot_live_biomass(:,:,icarbon)*veget_cov_max,dim=2)/1e3
2477	CALL histwrite_p (hist_id_stomate_IPCC, "cVeg", itime, &
2478	vartmp, npts, hori_index)
2479	!vartmp(:)=SUM(tot_litter_carb*veget_cov_max,dim=2)/1e3
2480	!CALL histwrite_p (hist_id_stomate_IPCC, "cLitter", itime, &
2481	! vartmp, npts, hori_index)
2482	!vartmp(:)=SUM(tot_soil_carb*veget_cov_max,dim=2)/1e3
2483	!CALL histwrite_p (hist_id_stomate_IPCC, "cSoil", itime, &
2484	! vartmp, npts, hori_index)
2485	!!$ vartmp(:)=(prod10_total + prod100_total)/1e3
2486	!!$ CALL histwrite_p (hist_id_stomate_IPCC, "cProduct", itime, &
2487	!!$ vartmp, npts, hori_index)
2488	!vartmp(:)=carb_mass_variation/1e3/one_day
2489	!CALL histwrite_p (hist_id_stomate_IPCC, "cMassVariation", itime, &
2490	! vartmp, npts, hori_index)
2491	vartmp(:)=SUM(lai*veget_cov_max,dim=2)
2492	CALL histwrite_p (hist_id_stomate_IPCC, "lai", itime, &
2493	vartmp, npts, hori_index)
2494	vartmp(:)=SUM(gpp_daily*veget_cov_max,dim=2)/1e3/one_day
2495	CALL histwrite_p (hist_id_stomate_IPCC, "gpp", itime, &
2496	vartmp, npts, hori_index)
2497	vartmp(:)=SUM((resp_maint+resp_growth)*veget_cov_max,dim=2)/1e3/one_day
2498	CALL histwrite_p (hist_id_stomate_IPCC, "ra", itime, &
2499	vartmp, npts, hori_index)
2500	vartmp(:)=SUM(npp_daily*veget_cov_max,dim=2)/1e3/one_day
2501	CALL histwrite_p (hist_id_stomate_IPCC, "npp", itime, &
2502	vartmp, npts, hori_index)
2503	vartmp(:)=SUM(resp_hetero*veget_cov_max,dim=2)/1e3/one_day
2504	CALL histwrite_p (hist_id_stomate_IPCC, "rh", itime, &
2505	vartmp, npts, hori_index)
2506	vartmp(:)=SUM(co2_fire*veget_cov_max,dim=2)/1e3/one_day
2507	CALL histwrite_p (hist_id_stomate_IPCC, "fFire", itime, &
2508	vartmp, npts, hori_index)
2509	!!$ vartmp(:)=harvest_above/1e3/one_day
2510	!!$ CALL histwrite_p (hist_id_stomate_IPCC, "fHarvest", itime, &
2511	!!$ vartmp, npts, hori_index)
2512	!!$ vartmp(:)=cflux_prod_total/1e3/one_day
2513	!!$ CALL histwrite_p (hist_id_stomate_IPCC, "fLuc", itime, &
2514	!!$ vartmp, npts, hori_index)
2515	!!$ vartmp(:)=(SUM((gpp_daily + atm_to_bm(:,:,icarbon) -(resp_maint+resp_growth+resp_hetero)-co2_fire) &
2516	!!$ & * veget_cov_max,dim=2)-cflux_prod_total-harvest_above)/1e3/one_day
2517	!!$ CALL histwrite_p (hist_id_stomate_IPCC, "nbp", itime, &
2518	!!$ vartmp, npts, hori_index)
2519	!vartmp(:)=SUM((tot_bm_to_litter(:,:,icarbon) + tot_turnover(:,:,icarbon))*&
2520	! veget_cov_max,dim=2)/1e3/one_day
2521	!CALL histwrite_p (hist_id_stomate_IPCC, "fVegLitter", itime, &
2522	! vartmp, npts, hori_index)
2523	vartmp(:)=SUM(SUM(som_input(:,:,:,icarbon),dim=2)*veget_cov_max,dim=2)/1e3/one_day
2524	CALL histwrite_p (hist_id_stomate_IPCC, "fLitterSoil", itime, &
2525	vartmp, npts, hori_index)
2526	vartmp(:)=SUM(SUM(circ_class_biomass(:,:,:,ileaf,icarbon)*&
2527	circ_class_n(:,:,:),3)*veget_cov_max,dim=2)/1e3
2528	CALL histwrite_p (hist_id_stomate_IPCC, "cLeaf", itime, &
2529	vartmp, npts, hori_index)
2530	vartmp(:)=SUM(SUM((circ_class_biomass(:,:,:,isapabove,icarbon)+&
2531	circ_class_biomass(:,:,:,iheartabove,icarbon))circ_class_n(:,:,:),3) &
2532	veget_cov_max,dim=2)/1e3
2533	CALL histwrite_p (hist_id_stomate_IPCC, "cWood", itime, &
2534	vartmp, npts, hori_index)
2535	vartmp(:)=SUM(SUM(( circ_class_biomass(:,:,:,iroot,icarbon) + &
2536	circ_class_biomass(:,:,:,isapbelow,icarbon) + &
2537	circ_class_biomass(:,:,:,iheartbelow,icarbon) )circ_class_n(:,:,:),3)&
2538	veget_cov_max,dim=2)/1e3
2539	CALL histwrite_p (hist_id_stomate_IPCC, "cRoot", itime, &
2540	vartmp, npts, hori_index)
2541	vartmp(:)=SUM(SUM((circ_class_biomass(:,:,:,icarbres,icarbon) + &
2542	circ_class_biomass(:,:,:,ilabile,icarbon) + &
2543	circ_class_biomass(:,:,:,ifruit,icarbon))circ_class_n(:,:,:),3)&
2544	veget_cov_max,dim=2)/1e3
2545	CALL histwrite_p (hist_id_stomate_IPCC, "cMisc", itime, &
2546	vartmp, npts, hori_index)
2547	vartmp(:)=SUM((litter(:,istructural,:,iabove,icarbon)+ &
2548	litter(:,imetabolic,:,iabove,icarbon)+litter(:,iwoody,:,iabove,icarbon))*&
2549	veget_cov_max,dim=2)/1e3
2550	CALL histwrite_p (hist_id_stomate_IPCC, "cLitterAbove", itime, &
2551	vartmp, npts, hori_index)
2552	vartmp(:)=SUM((litter(:,istructural,:,ibelow,icarbon)+ &
2553	litter(:,imetabolic,:,ibelow,icarbon)+litter(:,iwoody,:,ibelow,icarbon))*&
2554	veget_cov_max,dim=2)/1e3
2555	CALL histwrite_p (hist_id_stomate_IPCC, "cLitterBelow", itime, &
2556	vartmp, npts, hori_index)
2557	vartmp(:)=SUM((som(:,iactive,:,icarbon)+som(:,isurface,:,icarbon))* &
2558	veget_cov_max,dim=2)/1e3
2559	CALL histwrite_p (hist_id_stomate_IPCC, "cSoilFast", itime, &
2560	vartmp, npts, hori_index)
2561	vartmp(:)=SUM(som(:,islow,:,icarbon)*veget_cov_max,dim=2)/1e3
2562	CALL histwrite_p (hist_id_stomate_IPCC, "cSoilMedium", itime, &
2563	vartmp, npts, hori_index)
2564	vartmp(:)=SUM(som(:,ipassive,:,icarbon)*veget_cov_max,dim=2)/1e3
2565	CALL histwrite_p (hist_id_stomate_IPCC, "cSoilSlow", itime, &
2566	vartmp, npts, hori_index)
2567	DO j=1,nvm
2568	histvar(:,j)=veget_cov_max(:,j)*100
2569	ENDDO
2570	CALL histwrite_p (hist_id_stomate_IPCC, "landCoverFrac", itime, &
2571	histvar, npts*nvm, horipft_index)
2572	!-
2573	vartmp(:)=zero
2574	DO j = 2,nvm
2575	IF (is_deciduous(j)) THEN
2576	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2577	ENDIF
2578	ENDDO
2579	CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimDec", itime, &
2580	vartmp, npts, hori_index)
2581	!-
2582	vartmp(:)=zero
2583	DO j = 2,nvm
2584	IF (is_evergreen(j)) THEN
2585	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2586	ENDIF
2587	ENDDO
2588	CALL histwrite_p (hist_id_stomate_IPCC, "treeFracPrimEver", itime, &
2589	vartmp, npts, hori_index)
2590	!-
2591	vartmp(:)=zero
2592	DO j = 2,nvm
2593	IF ( .NOT.(is_c4(j)) ) THEN
2594	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2595	ENDIF
2596	ENDDO
2597	CALL histwrite_p (hist_id_stomate_IPCC, "c3PftFrac", itime, &
2598	vartmp, npts, hori_index)
2599	!-
2600	vartmp(:)=zero
2601	DO j = 2,nvm
2602	IF ( is_c4(j) ) THEN
2603	vartmp(:) = vartmp(:) + veget_cov_max(:,j)*100
2604	ENDIF
2605	ENDDO
2606	CALL histwrite_p (hist_id_stomate_IPCC, "c4PftFrac", itime, &
2607	vartmp, npts, hori_index)
2608	!-
2609	vartmp(:)=SUM(resp_growth*veget_cov_max,dim=2)/1e3/one_day
2610	CALL histwrite_p (hist_id_stomate_IPCC, "rGrowth", itime, &
2611	vartmp, npts, hori_index)
2612	vartmp(:)=SUM(resp_maint*veget_cov_max,dim=2)/1e3/one_day
2613	CALL histwrite_p (hist_id_stomate_IPCC, "rMaint", itime, &
2614	vartmp, npts, hori_index)
2615	vartmp(:)=SUM(bm_alloc(:,:,ileaf,icarbon)*veget_cov_max,dim=2)/1e3/one_day
2616	CALL histwrite_p (hist_id_stomate_IPCC, "nppLeaf", itime, &
2617	vartmp, npts, hori_index)
2618	vartmp(:)=SUM(bm_alloc(:,:,isapabove,icarbon)*veget_cov_max,dim=2)/1e3/one_day
2619	CALL histwrite_p (hist_id_stomate_IPCC, "nppWood", itime, &
2620	vartmp, npts, hori_index)
2621	vartmp(:)=SUM(( bm_alloc(:,:,isapbelow,icarbon) + bm_alloc(:,:,iroot,icarbon) )*&
2622	veget_cov_max,dim=2)/1e3/one_day
2623	CALL histwrite_p (hist_id_stomate_IPCC, "nppRoot", itime, &
2624	vartmp, npts, hori_index)
2625
2626	CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_X', itime, &
2627	resolution(:,1), npts, hori_index)
2628	CALL histwrite_p (hist_id_stomate_IPCC, 'RESOLUTION_Y', itime, &
2629	resolution(:,2), npts, hori_index)
2630	CALL histwrite_p (hist_id_stomate_IPCC, 'CONTFRAC', itime, &
2631	contfrac(:), npts, hori_index)
2632
2633	ENDIF
2634
2635	IF (printlev_loc>=3) WRITE (numout,*) 'soil_n_min(test_grid,test_pft,:):', &
2636	soil_n_min(test_grid,test_pft,:)
2637
2638	IF (printlev>=3) WRITE(numout,*) 'Leaving stomate_lpj'
2639
2640	END SUBROUTINE stomate_lpj_vegetation
2641
2642
2643	!! ================================================================================================================================
2644	!! SUBROUTINE : debug_write
2645	!!
2646	!>\BRIEF Write a set of variables to the out_orchidee file
2647	!!
2648	!! DESCRIPTION : Write a set of variables to the out_orchidee file after a routine
2649	!! a routine was called.
2650	!!
2651	!! RECENT CHANGE(S) : None
2652	!!
2653	!! MAIN OUTPUT VARIABLE(S): None
2654	!!
2655	!! REFERENCE(S) : None
2656	!!
2657	!! FLOWCHART : None
2658	!! \n
2659	!_ ================================================================================================================================
2660
2661	SUBROUTINE debug_write(npts, check_point, circ_class_biomass, &
2662	circ_class_n, circ_class_kill, plant_status, veget_cov_max)
2663
2664	!! 0. Variable and parameter description
2665
2666	!! 0.1 Input variables
2667	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
2668	REAL(r_std), DIMENSION(:,:), INTENT(in) :: plant_status !! Growth and phenological status of the plant
2669	!! istatus = Phases defined in constantes
2670	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Biomass of the componets of the model
2671	!! tree within a circumference
2672	!! class @tex $(gC ind^{-1})$ @endtex
2673	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_kill !! Number of trees within a circ class that needs
2674	!! to be killed @tex $(ind m^{-2})$ @endtex
2675	!! IMPORTANT: See the note in constantes.f90
2676	!! regarding the indicies for this variable.
2677	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Number of individuals in each circ class
2678	!! @tex $(m^{-2})$ @endtex
2679	CHARACTER(*),INTENT(in) :: check_point !! A flag to indicate at which
2680	!! point in the code we're doing
2681	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_cov_max !! "maximal" coverage fraction of a PFT
2682	!! (LAI -> infinity) on ground. May sum to
2683	!! less than unity if the pixel has
2684	!! nobio area. (unitless; 0-1)
2685
2686	!! 0.4 Local variables
2687	INTEGER :: icut, ifm, ipts, ivm
2688	!_ ================================================================================================================================
2689
2690	WRITE(numout,*) check_point
2691	WRITE(numout,*) 'test_pft, ', test_pft
2692	WRITE(numout,*) 'circ_class_n, ', SUM(circ_class_n(test_grid,test_pft,:))
2693	WRITE(numout,*) 'circ_class_biomass - labile N, ',&
2694	SUM(circ_class_biomass(test_grid,test_pft,:,ilabile,initrogen)*&
2695	circ_class_n(test_grid,test_pft,:),1)
2696	WRITE(numout,*) 'circ_class_biomass - C, ',&
2697	SUM(SUM(circ_class_biomass(test_grid,test_pft,:,:,icarbon),2)),&
2698	SUM(SUM(circ_class_biomass(test_grid,test_pft,:,:,icarbon),2)*&
2699	circ_class_n(test_grid,test_pft,:),1)
2700	WRITE(numout,*) 'circ_class_biomass - N, ',&
2701	SUM(SUM(circ_class_biomass(test_grid,test_pft,:,:,initrogen),2)),&
2702	SUM(SUM(circ_class_biomass(test_grid,test_pft,:,:,initrogen),2)*&
2703	circ_class_n(test_grid,test_pft,:),1)
2704	WRITE(numout,*) 'plant_status, ', plant_status(:,test_pft)
2705	! DO ifm = 1, nfm_types
2706	! DO icut = 1, ncut_times
2707	! WRITE(numout,*) 'circ_class_kill, ifm, icut, ', ifm, icut, &
2708	! circ_class_kill(test_grid,test_pft,:,ifm,icut)
2709	! ENDDO
2710	! ENDDO
2711	WRITE(numout,*) 'veget_cov_max, ', veget_cov_max(:,test_pft), SUM(veget_cov_max(:,:),2)
2712
2713	DO ipts = 1,npts
2714	DO ivm = 1,nvm
2715	IF (veget_cov_max(ipts,ivm) .GT. min_stomate .AND. &
2716	SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1)) .LT. min_stomate .AND.&
2717	SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1)) .NE. zero) THEN
2718	WRITE(numout,*) check_point
2719	WRITE(numout,*) 'veget_cov_max, ',ivm, veget_cov_max(ipts,ivm)
2720	WRITE(numout,*) 'biomass -C, ', &
2721	SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1))
2722	WRITE(numout,*) 'biomass -N, ', &
2723	SUM(SUM(circ_class_biomass(ipts,ivm,:,:,initrogen),1))
2724	WRITE(numout,*) 'WARNING - some very small amount of C left'
2725	CALL ipslerr_p(2,'stomate_lpj','Very small amount of C left',&
2726	'Could cause problems elsewhere','')
2727	ENDIF
2728	ENDDO
2729	ENDDO
2730
2731
2732	END SUBROUTINE debug_write
2733
2734	END MODULE stomate_lpj

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

Download in other formats: