Context Navigation

source: branches/publications/ORCHIDEE_CAN_r3069/src_stomate/stomate_prescribe.f90

Last change on this file was 3069, checked in by sebastiaan.luyssaert, 9 years ago
PROD: tested gloabally for the zoomed European grid for 20 years. Do not use the previous commit as there were problems with svn it does not contain the described changes. This commit contains the bug fixes to species change and added functionality to change forest management following mortality or a harvest.
File size: 78.7 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_prescribe
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Initialize and update density, height, basal area and crown area. For the resource
10	! limited allocation and the allometric-based allocation
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S): None
15	!!
16	!! REFERENCE(S) :
17	!!
18	!! SVN :
19	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_stomate/stomate_prescribe.f90 $
20	!! $Date: 2013-01-04 16:50:56 +0100 (Fri, 04 Jan 2013) $
21	!! $Revision: 1126 $
22	!! \n
23	!_ ================================================================================================================================
24
25	MODULE stomate_prescribe
26
27	! modules used:
28
29	USE ioipsl_para
30	USE stomate_data
31	USE pft_parameters
32	USE function_library, ONLY: calculate_c0_alloc
33	USE constantes
34
35	IMPLICIT NONE
36
37	! private & public routines
38
39	PRIVATE
40	PUBLIC prescribe_diagnostic, prescribe_prognostic, prescribe_clear
41
42	! first call
43	LOGICAL, SAVE :: firstcall = .TRUE.
44	!$OMP THREADPRIVATE(firstcall)
45
46	CONTAINS
47
48
49	! =================================================================================================================================
50	!! SUBROUTINE : prescribe_clear
51	!!
52	!>\BRIEF : Set the firstcall flag back to .TRUE. to prepare for the next simulation.
53	!_=================================================================================================================================
54
55	SUBROUTINE prescribe_clear
56	firstcall=.TRUE.
57	END SUBROUTINE prescribe_clear
58
59
60
61	!! ================================================================================================================================
62	!! SUBROUTINE : prescribe_diagnostic
63	!!
64	!>\BRIEF Works only with static vegetation and agricultural PFT. Initialize biomass,
65	!! density, crown area in the first call and update them in the following calls.
66	!!
67	!! DESCRIPTION (functional, design, flags): \n
68	!! This module works only with static vegetation and agricultural PFT. In the first call, initialize
69	!! density of individuals, biomass, crown area, and leaf age distribution to some reasonable value.
70	!! In the following calls, these variables are updated.
71	!!
72	!! To fulfill these purposes, relationships from the pipe model are used:
73	!! \latexonly
74	!! \input{prescribe1.tex}
75	!! \input{prescribe2.tex}
76	!! \endlatexonly
77	!!
78	!! RECENT CHANGE(S): None
79	!!
80	!! MAIN OUTPUT VARIABLES(S): ::ind, ::cn_ind, ::leaf_frac
81	!!
82	!! REFERENCES :
83	!! - Krinner, G., N. Viovy, et al. (2005). "A dynamic global vegetation model
84	!! for studies of the coupled atmosphere-biosphere system." Global
85	!! Biogeochemical Cycles 19: GB1015, doi:1010.1029/2003GB002199.
86	!! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics,
87	!! plant geography and terrestrial carbon cycling in the LPJ dynamic
88	!! global vegetation model, Global Change Biology, 9, 161-185.
89	!!
90	!! FLOWCHART : None
91	!! \n
92	!_ ================================================================================================================================
93
94	SUBROUTINE prescribe_diagnostic (npts, veget_max, dt, PFTpresent, everywhere, &
95	when_growthinit, biomass, leaf_frac, ind, cn_ind, &
96	co2_to_bm)
97
98
99	!! 0. Parameters and variables declaration
100
101	!! 0.1 Input variables
102
103	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
104	REAL(r_std), INTENT(in) :: dt !! time step (dt_days)
105	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT
106	!! (LAI -> infinity) on ground.
107	!! (unitless;0-1)
108	!! 0.2 Output variables
109
110	!! 0.3 Modified variables
111
112	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: PFTpresent !! PFT present (0 or 1)
113	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
114	!! very localized (after its introduction)
115	!! (unitless;0-1)
116	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: when_growthinit !! how many days ago was the beginning of
117	!! the growing season (days)
118	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
119	INTENT(inout) :: biomass !! biomass @tex $(gC.m^{-2})$ @endtex
120	REAL(r_std), DIMENSION(npts,nvm,nleafages), &
121	INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age
122	!! class (unitless;0-1)
123	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind !! density of individuals
124	!! @tex $(m^{-2})$ @endtex
125	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: cn_ind !! crown area per individual
126	!! @tex $(m^{2})$ @endtex
127	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: co2_to_bm !! co2 taken up by carbohydrate
128	!! reserve at the beginning of the
129	!! growing season @tex (gC.m^{-2}dt^{-1})$ @endtex
130
131	!! 0.4 Local variables
132
133	REAL(r_std), DIMENSION(npts) :: dia !! stem diameter (m)
134	REAL(r_std), DIMENSION(npts,nelements) :: woodmass !! woodmass @tex $(gC.m^{-2})$ @endtex
135	REAL(r_std), DIMENSION(npts,nelements) :: woodmass_ind !! woodmass of an individual
136	!! @tex $(gC.tree^{-1})$ @endtex
137	INTEGER(i_std) :: i,j !! index (unitless)
138
139	!_ ================================================================================================================================
140
141	!! 1. Calculations at every call of the routine
142
143	DO j = 2,nvm
144
145	! only when the DGVM is not activated or agricultural PFT.
146	IF ( ( .NOT. control%ok_dgvm .AND. control%ok_constant_mortality ) .OR. ( .NOT. natural(j) ) ) THEN
147
148	!! 1.1 Update crown area
149	cn_ind(:,j) = zero
150
151	IF ( is_tree(j) ) THEN
152
153	!! 1.1.1 Trees
154	dia(:) = zero
155	DO i = 1, npts ! loop over grid points
156
157	IF ( veget_max(i,j) .GT. zero ) THEN
158
159	!! 1.1.1.1 calculate wood mass on an area basis, which include sapwood and heartwood aboveground and belowground.
160	woodmass(i,icarbon) = (biomass(i,j,isapabove,icarbon) + biomass(i,j,isapbelow,icarbon) + &
161	biomass(i,j,iheartabove,icarbon) + biomass(i,j,iheartbelow,icarbon)) * veget_max(i,j)
162
163	IF ( woodmass(i,icarbon) .GT. min_stomate ) THEN
164
165	!! 1.1.1.2 calculate critical individual density
166	!?? the logic for 1.1.3 and 1.1.2 is strange, it should be the case that first to calculate critical woodmass
167	!??per individual, then calculate critical density.
168	! how to derive the following equation:
169	! first, TreeHeight=pipe_tune2 * Diameter^{pipe_tune3}
170	! we assume the tree is an ideal cylinder, so it volume is:
171	! Volume = pi(Dia/2)^2Height = pi/4 * Dia * pipe_tune2*Dia^{pipe_tune3}
172	! = pi/4 * pipe_tune2 * Dia^{2+pipe_tune3}
173	! last, the woodmass_per_individual = pipe_density * Volume = pipe_densitypi/4.pipe_tune2 * Dia^{2+pipe_tune3}
174	ind(i,j) = woodmass(i,icarbon) / &
175	( pipe_density(j)pi/4.pipe_tune2(j) * maxdia(j)**(2.+pipe_tune3(j)) )
176
177
178	!! 1.1.1.3 individual biomass corresponding to this critical density of individuals
179	woodmass_ind(i,icarbon) = woodmass(i,icarbon) / ind(i,j)
180
181	!! 1.1.1.4 calculate stem diameter per individual tree
182	dia(i) = ( woodmass_ind(i,icarbon) / ( pipe_density(j) * pi/4. * pipe_tune2(j) ) ) ** &
183	( un / ( 2. + pipe_tune3(j) ) )
184
185	!! 1.1.1.5 calculate provisional tree crown area for per individual tree
186	! equation: CrownArea=pipe_tune1 * Diameter^{1.6}
187	cn_ind(i,j) = pipe_tune1(j) * MIN( maxdia(j), dia(i) ) ** pipe_tune_exp_coeff(j)
188
189	!! 1.1.1.6 If total tree crown area for this tree PFT exceeds its veget_max, tree density is recalculated.
190	IF ( cn_ind(i,j) * ind(i,j) .GT. 1.002* veget_max(i,j) ) THEN
191
192	ind(i,j) = veget_max(i,j) / cn_ind(i,j)
193
194	ELSE
195
196	ind(i,j) = ( veget_max(i,j) / ( pipe_tune1(j) * (woodmass(i,icarbon)/(pipe_density(j)pi/4.pipe_tune2(j))) &
197	& **(pipe_tune_exp_coeff(j)/(2.+pipe_tune3(j))) ) ) &
198	& ** (1./(1.-(pipe_tune_exp_coeff(j)/(2.+pipe_tune3(j)))))
199	woodmass_ind(i,icarbon) = woodmass(i,icarbon) / ind(i,j)
200	dia(i) = ( woodmass_ind(i,icarbon) / ( pipe_density(j) * pi/4. * pipe_tune2(j) ) ) ** &
201	& ( un / ( 2. + pipe_tune3(j) ) )
202	cn_ind(i,j) = pipe_tune1(j) * MIN( maxdia(j), dia(i) ) ** pipe_tune_exp_coeff(j)
203
204	ENDIF
205
206	ELSE ! woodmass=0 => impose some value
207
208	dia(:) = maxdia(j)
209	cn_ind(i,j) = pipe_tune1(j) * MIN( maxdia(j), dia(i) ) ** pipe_tune_exp_coeff(j)
210
211	ENDIF ! IF ( woodmass .GT. min_stomate )
212
213	ENDIF ! veget_max .GT. 0.
214
215	ENDDO ! loop over grid points
216
217	!! 1.2 grasses: crown area always set to 1m**2
218	ELSE
219
220	WHERE ( veget_max(:,j) .GT. zero )
221
222	cn_ind(:,j) = un
223
224	ENDWHERE
225
226	ENDIF !IF ( is_tree(j) )
227
228	!! 1.3 density of individuals
229	WHERE ( veget_max(:,j) .GT. zero )
230
231	ind(:,j) = veget_max(:,j) / cn_ind(:,j)
232
233	ELSEWHERE
234
235	ind(:,j) = zero
236
237	ENDWHERE
238
239	ENDIF ! IF ( ( .NOT. control%ok_dgvm .AND. control%ok_constant_mortality ) .OR. ( .NOT. natural(j) ) )
240
241	ENDDO ! loop over PFTs
242
243
244	!! 2. If it's the first call for this module,
245
246	IF (( firstcall ) .AND. (TRIM(stom_restname_in) == 'NONE')) THEN
247
248	! impose some biomass if zero and PFT prescribed
249	WRITE(numout,*) 'prescribe:'
250	WRITE(numout,*) ' > Imposing initial biomass for prescribed trees, '// &
251	'initial reserve mass for prescribed grasses.'
252	WRITE(numout,*) ' > Declaring prescribed PFTs present.'
253
254	DO j = 2,nvm ! loop over PFTs
255
256	DO i = 1, npts ! loop over grid points
257
258	! is vegetation static or PFT agricultural?
259	! Static vegetation or agricultural PFT
260	IF ( ( .NOT. control%ok_dgvm ) .OR. &
261	( ( .NOT. natural(j) ) .AND. ( veget_max(i,j) .GT. min_stomate ) ) ) THEN
262
263
264	!! 2.1 if tree biomass is extremely small, prescribe the biomass by assuming they have sapling biomass,
265	!! which is a constant in the model then set all the leaf age as 1.
266	!! if tree PFT and biomass too small, prescribe the biomass to a value.
267	IF ( is_tree(j) .AND. ( veget_max(i,j) .GT. min_stomate ) .AND. &
268	( SUM( biomass(i,j,:,icarbon) ) .LE. min_stomate ) ) THEN
269
270	!!? here the code is redundant, as "veget_max(i,j) .GT. min_stomate" is already met in the above if
271	!!? condition.
272	IF (veget_max(i,j) .GT. min_stomate) THEN
273
274	biomass(i,j,:,:) = (bm_sapl_rescale * bm_sapl_old(j,:,:) * ind(i,j)) / veget_max(i,j)
275
276	ELSE
277
278	biomass(i,j,:,:) = zero
279
280	ENDIF
281
282	! set leaf age classes
283	leaf_frac(i,j,:) = zero
284	leaf_frac(i,j,1) = un
285
286	! set time since last beginning of growing season
287	when_growthinit(i,j) = large_value
288
289	! seasonal trees: no leaves at beginning
290	IF ( pheno_model(j) .NE. 'none' ) THEN
291
292	biomass(i,j,ileaf,icarbon) = zero
293	leaf_frac(i,j,1) = zero
294
295	ENDIF
296
297	co2_to_bm(i,j) = co2_to_bm(i,j) + ( SUM(biomass(i,j,:,icarbon)) / dt )
298
299	ENDIF
300
301	!! 2.2 for grasses, set only the carbon reserve pool to "sapling" carbon reserve pool.
302	!! and set all leaf age to 1.
303	IF ( ( .NOT. is_tree(j) ) .AND. ( veget_max(i,j) .GT. min_stomate ) .AND. &
304	( SUM( biomass(i,j,:,icarbon) ) .LE. min_stomate ) ) THEN
305
306	!+++++ CHECK MJM ++++++
307	! bm_sapl will also be a function of the circumference class...just using the first class right now
308	biomass(i,j,icarbres,:) = bm_sapl_old(j,icarbres,:) * ind(i,j) / veget_max(i,j)
309
310	! set leaf age classes
311	leaf_frac(i,j,:) = zero
312	leaf_frac(i,j,1) = un
313
314	! set time since last beginning of growing season
315	when_growthinit(i,j) = large_value
316
317	co2_to_bm(i,j) = co2_to_bm(i,j) + ( biomass(i,j,icarbres,icarbon) / dt )
318
319	ENDIF
320
321	!! 2.3 declare all PFTs with positive veget_max as present everywhere in that grid box
322	IF ( veget_max(i,j) .GT. min_stomate ) THEN
323
324	PFTpresent(i,j) = .TRUE.
325	everywhere(i,j) = un
326
327	ENDIF
328
329	ENDIF ! not control%ok_dgvm or agricultural
330
331	ENDDO ! loop over grid points
332
333	ENDDO ! loop over PFTs
334
335	firstcall = .FALSE.
336
337	ENDIF
338
339	END SUBROUTINE prescribe_diagnostic
340
341
342	!! ================================================================================================================================
343	!! SUBROUTINE : prescribe_prognostic
344	!!
345	!>\BRIEF Works only with static vegetation and agricultural PFT. Initialize biomass,
346	!! density, presence in the first call and update them in the following.
347	!!
348	!! DESCRIPTION (functional, design, flags): \n
349	!! This module works only with static vegetation and agricultural PFT.
350	!! In the first call, initialize density of individuals, biomass,
351	!! and leaf age distribution to some reasonable value. In the following calls,
352	!! these variables are updated.
353	!!
354	!! To fulfill these purposes, pipe model are used:
355	!! \latexonly
356	!! \input{prescribe1.tex}
357	!! \input{prescribe2.tex}
358	!! \endlatexonly
359	!!
360	!! RECENT CHANGE(S): None
361	!!
362	!! MAIN OUTPUT VARIABLES(S): ::ind, ::cn_ind, ::leaf_frac
363	!!
364	!! REFERENCES :
365	!! - Krinner, G., N. Viovy, et al. (2005). "A dynamic global vegetation model
366	!! for studies of the coupled atmosphere-biosphere system." Global
367	!! Biogeochemical Cycles 19: GB1015, doi:1010.1029/2003GB002199.
368	!! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics,
369	!! plant geography and terrestrial carbon cycling in the LPJ dynamic
370	!! global vegetation model, Global Change Biology, 9, 161-185.
371	!! - McDowell, N., Barnard, H., Bond, B.J., Hinckley, T., Hubbard, R.M., Ishii,
372	!! H., KÃ¶stner, B., Magnani, F. Marshall, J.D., Meinzer, F.C., Phillips, N.,
373	!! Ryan, M.G., Whitehead D. 2002. The relationship between tree height and leaf
374	!! area: sapwood area ratio. Oecologia, 132:12â20.
375	!! - Novick, K., Oren, R., Stoy, P., Juang, F.-Y., Siqueira, M., Katul, G. 2009.
376	!! The relationship between reference canopy conductance and simplified hydraulic
377	!! architecture. Advances in water resources 32, 809-819.
378	!!
379	!! FLOWCHART : None
380	!! \n
381	!_ ================================================================================================================================
382
383	SUBROUTINE prescribe_prognostic (npts, veget_max, veget, dt, PFTpresent, &
384	everywhere, when_growthinit, biomass, leaf_frac, &
385	ind, circ_class_n, circ_class_biomass, co2_to_bm, &
386	forest_managed, KF, &
387	senescence, age, npp_longterm, lm_lastyearmax, &
388	tau_eff_leaf, tau_eff_sap, tau_eff_root, k_latosa_adapt,&
389	lpft_replant)
390
391	!! 0. Parameters and variables declaration
392
393	!! 0.1 Input variables
394
395	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
396	INTEGER(i_std), DIMENSION (:,:), INTENT(in) :: forest_managed !! forest management flag
397	REAL(r_std), INTENT(in) :: dt !! time step (dt_days)
398	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT
399	!! (LAI -> infinity) on ground. May sum to
400	!! less than unity if the pixel has
401	!! nobio area. (unitless; 0-1)
402	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget !! Fraction of vegetation type including
403	!! non-biological fraction (unitless)
404	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_root !! Effective root turnover time that accounts
405	!! waterstress (days)
406	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_sap !! Effective sapwood turnover time that accounts
407	!! waterstress (days)
408	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_leaf !! Effective leaf turnover time that accounts
409	!! waterstress (days)
410	LOGICAL, DIMENSION(npts,nvm), OPTIONAL, INTENT(in):: lpft_replant !! Set to true if a PFT has been clearcut
411	!! and needs to be replaced by another species
412
413
414	!! 0.2 Output variables
415
416
417	!! 0.3 Modified variables
418
419	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! PFT present (0 or 1)
420	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
421	!! for deciduous trees)
422	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
423	!! very localized (after its introduction) (?)
424	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! how many days ago was the beginning of
425	!! the growing season (days)
426	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Density of individuals at the stand level
427	!! @tex $(m^{-2})$ @endtex
428	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "long term" net primary productivity
429	!! @tex ($gC m^{-2} year^{-1}$) @endtex
430	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
431	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
432	!! @tex ($gC m^{-2}$) @endtex
433	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: circ_class_n !! Number of individuals in each circ class
434	!! @tex $(ind m^{-2})$ @endtex
435	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass
436	!! tex $(gC.m^{-2})$ @endtex
437	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: circ_class_biomass !! Biomass components of the model tree
438	!! within a circumference class
439	!! class @tex $(g C ind^{-1})$ @endtex
440	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C
441	!! to create the seedlings
442	!! @tex (gC.m^{-2}dt^{-1})$ @endtex
443	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age
444	!! class (unitless;0-1)
445	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: KF !! Scaling factor to convert sapwood mass
446	!! into leaf mass (m) - this variable is
447	!! passed to other routines
448	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: k_latosa_adapt !! Leaf to sapwood area adapted for long
449	!! term water stress (m)
450
451
452	!! 0.4 Local variables
453
454	REAL(r_std), DIMENSION(npts,nvm) :: c0_alloc !! Root to sapwood tradeoff parameter
455	INTEGER(i_std) :: ipts, ivm, ipar !! index (unitless)
456	INTEGER(i_std) :: icir, iele, imbc !! index (unitless)
457	INTEGER(i_std) :: deb,fin, imaxt !! index (unitless)
458	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements) :: sync_biomass !! Temporary stand level biomass
459	!! @tex $(gC.m^{-2})$ @endtex
460	REAL(r_std), DIMENSION(npts,nvm) :: sync_ind !! Temporary density of individuals at the
461	!! stand level @tex $(m^{-2})$ @endtex
462	REAL(r_std), DIMENSION(npts,nvm) :: k_latosa !! Height dependent base value to calculate
463	!! KF (-)
464	REAL(r_std), DIMENSION(nvm,ncirc,nparts,nelements) :: bm_sapl !! Sapling biomass for the functional
465	!! allocation with a dimension for the
466	!! circumference classes
467	!! @tex $(gC.ind^{-1})$ @endtex
468	REAL(r_std), DIMENSION(npts,nvm) :: LF !! Scaling factor to convert sapwood mass
469	!! into root mass (unitless)
470	REAL(r_std), DIMENSION(npts,nvm) :: lstress_fac !! Light stress factor, based on total
471	!! transmitted light (unitless, 0-1)
472	REAL(r_std), DIMENSION(npts,nvm,ncirc) :: circ_class !! circumference of individual trees
473	REAL(r_std) :: lambda !! lambda of the truncated exponential
474	!! distribution of initial diameters
475	!! (1/cm**2). Parameterized based on Dhote
476	!! (2003): lambda = sqrt(2)/Dg and
477	!! Dginit = 1 cm.
478	REAL(r_std) :: p_max !! Cumulated probability level at which the
479	!! exponential distribution truncated
480	!! (dimensionless, 0-1)
481	REAL(r_std),DIMENSION(ncirc) :: circ_ij0 !! Circumference of all individual trees
482	REAL(r_std),DIMENSION(ncirc) :: height_ij0 !! Height of each Circumference of all
483	!! individual trees
484	!! for 1 ha at gridpoint i and for PFT j
485	!! (m). As above, 0 is before thinning
486	REAL(r_std), DIMENSION(nparts) :: bm_init !! Biomass needed to initiate the next
487	!! planting @tex $(gC m^{-2})$ @endtex
488	REAL(r_std), DIMENSION(ncirc) :: nb_trees_i !! Number of trees in each twentith
489	!! circumference quantile of the
490	!! distribution (ind)
491	REAL(r_std) :: excedent !! Number of trees after truncation to be
492	!! reallocated to smaller quantiles of the
493	!! distribution (ind)
494	REAL(r_std) :: max_circ_init !! Maximum initial circumferences of the
495	!! truncated exponential distribution (cm)
496	REAL(r_std), DIMENSION(ncirc) :: ave_tree_height !! The height of the ideal tree in each
497	!! circumference class of the
498	!! distribution (ind)...not saved since
499	!! it should only be used for prescribing
500	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
501	!! mass balance chech for this routine
502	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
503	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
504	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
505	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start and end pool of this routine
506	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
507	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! Start and end pool of this routine
508	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
509	REAL(r_std) :: delta_KF !! Difference between new and old estimate
510	!! of KF while iterating
511	REAL(r_std) :: min_height_init
512	REAL(r_std) :: max_height_init
513	REAL(r_std) :: sapwood_density
514	REAL(r_std) :: dia_init
515	REAL(r_std) :: wood_init
516	REAL(r_std), DIMENSION(ncirc) :: height_dist
517	INTEGER(i_std) :: iloop
518
519	!_ ================================================================================================================================
520
521	!! 1. Initialize biomass at first call
522
523	!! 1.1.1 Initialize check for mass balance closure
524	! The mass balance is calculated at the end of this routine
525	! in section 4
526	! Initial biomass pool
527	pool_start(:,:,:) = zero
528	DO ipar = 1,nparts
529	DO iele = 1,nelements
530	pool_start(:,:,iele) = pool_start(:,:,iele) + &
531	(biomass(:,:,ipar,iele) * veget_max(:,:))
532	ENDDO
533	ENDDO
534
535	! co2_to_bm is has as intent inout, the variable accumulates
536	! carbon over the course of a day. Use the difference between
537	! start and the end of this routine
538	check_intern(:,:,iatm2land,icarbon) = - un * co2_to_bm(:,:) * veget_max(:,:) * dt
539
540	! Prescribe was taken out of the first call loop. Because of the firstcall here,
541	! nothing was done when the vegetation was killed in lpj_gap and the vegetation
542	! thus stayed dead forever. We want it to regrow. So instead, let's loop over
543	! every point and PFT. If there is supposed to be vegetation here (according to
544	! veget_max) but there isn't (according to ind), then we grow it.
545	DO ivm = 2,nvm ! Loop over PFTs
546
547	DO ipts = 1, npts ! Loop over pixels
548
549	! If we are regrowing different species after managed forests
550	! are killed/die, we sometimes need to prevent them regrowing.
551	! This is only true if they die in the middle of the year due to natural
552	! causes, in which case they will be replanted at the end of the year.
553	! This loop checks to see if we regrow this PFT now or later. We can
554	! regrow a PFT after lpft_replant is set back to false so before that
555	! we need to process all the information such that the correct species
556	! will be regrown. lpft_replant is an optional variable. First test
557	! whether it is present.
558	IF(lchange_species)THEN
559	IF(PRESENT(lpft_replant))THEN
560	IF(lpft_replant(ipts,ivm))THEN
561	CYCLE
562	ENDIF
563	ENDIF
564	ENDIF
565
566	! We are supposed to have vegetation, but we don't have any: prescribe some.
567	IF((veget_max(ipts,ivm) .GT. min_stomate) .AND. (ind(ipts,ivm) .LE. min_stomate))THEN
568
569	! Initilaize tree level variables
570	circ_class_n(ipts,ivm,:) = zero
571	biomass(ipts,ivm,:,:) = zero
572	circ_class_biomass(ipts,ivm,:,:,:) = zero
573
574	! Assume that there is no memory for k_latosa between
575	! different generations (this is probably true when trees
576	! are planted but seems less likely for natural regeneration
577	! Allowing for memory has caused problems with the LAI from
578	! the second generation onwards
579	k_latosa_adapt(ipts,ivm) = k_latosa_min(ivm)
580
581	! Write output
582	IF (ld_presc) THEN
583	WRITE(numout,*) 'Prescribe (stomate_prescribe.f90):'
584	WRITE(numout,*) ' > Imposing initial biomass for prescribed trees, '// &
585	'initial reserve mass for prescribed grasses.'
586	WRITE(numout,*) ' > Declaring prescribed PFTs present.'
587	WRITE(numout,*) ' > Grid point: ',ipts,' PFT Type: ',ivm
588	ENDIF
589
590	!! 2. Calculate the vegetation characteristics of a newly established vegetation
591
592	!! 2.1 Stress factors
593
594	! Note that light and water stress have an opposite effect on KF. Waterstress
595	! will decrease KF because more C should be allocated to the roots. Light
596	! stress will increase KF because more C should be allocated to the leaves.
597
598	! We will need the c0_alloc, it only depends on PFT and the effective
599	! longiveties
600	c0_alloc(ipts,ivm) = calculate_c0_alloc(ipts, ivm, tau_eff_root(ipts,ivm), &
601	tau_eff_sap(ipts,ivm))
602
603	! Lightstress varies from 0 to 1 and is calculated from the canopy structure (veget)
604	! Given that there is no vegetation at this point, the lightstress cannot be calculated
605	! and is set to 1. However as soon as we grow a canopy it will experience light stress
606	! and so KF should be adjusted. If this is not done, we can't prescribe tall vegetation
607	! which is a useful feature to speed up optimisation and testing. We will have iterate
608	! over light stress and KF to prescribe vegetation that is in balance with its
609	! light environment.
610	lstress_fac(ipts,ivm) = un
611
612	! Initial vegetation has to be prescribed only when the vegetation is static
613	IF ( ( .NOT. control%ok_dgvm ) .AND. &
614	( veget_max(ipts,ivm) .GT. min_stomate ) ) THEN
615
616	!! 2.2 Initialize woody PFT's
617	! Use veget_max to check whether the PFT is present. If the PFT is present but it has no
618	! biomass, prescribe its biomass (gC m-{2}).
619	IF ( is_tree(ivm) .AND. &
620	( veget_max(ipts,ivm) .GT. min_stomate ) .AND. &
621	( SUM( biomass(ipts,ivm,:,icarbon) ) .LE. min_stomate ) ) THEN
622
623	! PFT is present so it needs to be initialized
624	IF (veget_max(ipts,ivm) .GT. min_stomate) THEN
625
626	! The forest stand starts with the number of individuals as prescribed in constantes_mtc.f90
627	! This number is defined per hectare whereas these calculations are per m2
628	ind(ipts,ivm) = nmaxtrees(ivm) * ha_to_m2
629
630	!! 2.3 Initialize height distribution
631	! Initializing height distribution...I've taken this from the circumference initilization of Valentin
632	! we want ncirc bins between height_init_min and height_init_max
633	lambda = SQRT(2.0_r_std)/pi
634	min_height_init = height_init_min(ivm)
635	max_height_init = height_init_max(ivm)
636	p_max = 100./nmaxtrees(ivm)
637	height_dist(:) = zero
638
639	!! 2.4 Calculate height distribution
640
641	!! 2.4.1 Height distribution of a high stand
642	! Deterministic initial distribution following a truncated exponential law
643	! if they are coppices, we handle this differently
644
645	!++++ CHECK+++++
646	! I am treating SRC the same way as everything else for now.
647	!IF (forest_managed(ipts,ivm) .NE. ifm_src) THEN
648	nb_trees_i(:) = 0
649
650	! The distribution is truncated when P(X>max_circ_init)<p_max
651	! The circumference distribution is divided into ncirc quantiles of equal height range.
652	! An integer number of trees is
653	! allocated to each quantile according to the truncated exponential distribution.
654	DO icir = 1,ncirc ! loop over circumference quantiles
655
656	nb_trees_i(icir) = NINT(nmaxtrees(ivm)*&
657	(EXP(-lambdamax_height_init(icir-1)/ncirc)-EXP(-lambdamax_height_initicir/ncirc)))
658	! this array is going to be used to create the sapling for this class so that we can create
659	! the total amount of biomass
660	! in this class...after this we will always calculate the height/circ/diam/whatever
661	! from the biomass in the class
662	ave_tree_height(icir)=min_height_init+(REAL(icir,r_std)-0.5_r_std)*&
663	(max_height_init-min_height_init)/REAL(ncirc,r_std)
664
665	! I am changing this so that our lowest class is min_height_init and our biggest class
666	! is max_height_init. With this change, we can use the same code the redistribute
667	! the trees if the biomass in one of our height classes is equal to zero later on due to mortality or
668	! thinning.
669	!+++ This caused some strange behavior, so change it back for now.
670	! ave_tree_height(icir)=min_height_init+REAL(icir-1,r_std)/REAL(ncirc-1,r_std)*&
671	! (max_height_init-min_height_init)
672
673	ENDDO
674
675	excedent = REAL(nmaxtrees(ivm) - SUM(nb_trees_i))
676	nb_trees_i = nb_trees_i + FLOOR(excedent/ncirc)
677	nb_trees_i(1) = nb_trees_i(1) + NINT(excedent)-ncirc*FLOOR(excedent/ncirc)
678
679	IF (ld_presc) THEN
680	WRITE(numout,*)"min initial height ",min_height_init,'max initial height ',max_height_init
681	WRITE(numout,*)'nb_trees_i',nb_trees_i
682	ENDIF
683
684	! ENDIF ! check FM type
685
686	! we need to store the number of trees that we have per square meter for each circumference class
687	circ_class_n(ipts,ivm,:) = REAL(nb_trees_i(:),r_std) * ha_to_m2
688
689	!! 2.4.2 Circumference distribution of a coppice
690	! I'm not see why this should be any different than a standard prescribe.
691	! Ideally, the plantings should be stems 20-25 cm long stuck into the soil,
692	! but the allocation scheme will have the same issues with this that it does
693	! with planting a small sapling, i.e. it won't like it. The stems do not
694	! follow standard allocation rules. Then again, neither do coppices.
695
696	!!$ IF (forest_managed(ipts,ivm) == ifm_src) THEN
697	!!$ CALL ipslerr(3,'stomate_prescribe.f90','Problem with SRC','','')
698	!!$ ENDIF
699
700	!! 2.5 Allocation factors
701	! Sapwood to root ratio
702	! Following Magnani et al. 2000 "In order to decreases hydraulic resistance,
703	! the investment of carbon in fine roots or sapwood yields to the plant very
704	! different returns, both because of different hydraulic conductivities and
705	! because of the strong impact of plant height on shoot resistance. On the
706	! other hand, fine roots and sapwood have markedly different longevities and
707	! the cost of production, discounted for turnover, will differ accordingly.
708	! Optimal growth under hydraulic constraints requires that the ratio of marginal
709	! hydraulic returns to marginal annual cost for carbon investment in either roots
710	! or sapwood be the same (Bloom, Chapin & Mooney 1985; Case & Fair 1989). This
711	! is formalized in equation (13) and further derived to obtain equation (17)
712	! in Magnani et al 2000. The latter is implemented here.
713	! Pipe_density is given in gC/m-3, convert to kg/m-3. And apply equation (17)
714	! in Magnani et al 2000. Note that c0_alloc was calculated at the start of this
715	! routine. The calculation itself is done in function_library
716
717	! Calculate leaf area to sapwood area
718	! To be consistent with the hydraulic limitations and pipe theory,
719	! k_latosa is calculated from equation (18) in Magnani et al.
720	! To do so, total hydraulic resistance and tree height need to known. This
721	! poses a problem as the resistance depends on the leaf area and the leaf
722	! area on the resistance. There is no independent equation and equations 12
723	! and 18 depend on each other and substitution would be circular. Hence
724	! prescribed k_latosa values were obtained from observational records
725	! and are given in mtc_parameters.f90.
726
727	! The relationship between height and k_latosa as reported in McDowell
728	! et al 2002 and Novick et al 2009 is implemented to adjust k_latosa for
729	! the height of the stand. This did NOT result in a realistic model behavior
730	!!$ k_latosa(ipts,ivm) = wstress_fac(ipts,ivm) * &
731	!!$ (k_latosa_max(ivm) - (latosa_height(ivm) * &
732	!!$ (SUM( nb_trees_i(:) * ave_tree_height(:) ) / SUM( nb_trees_i(:) ))))
733
734	! Alternatively, k_latosa is also reported to be a function of diameter
735	! (i.e. stand thinning, Simonin et al 2006, Tree Physiology, 26:493-503).
736	! Here the relationship with thinning was interpreted as a realtionship with
737	! light stress. Note that light stress cannot be calculated at this time in
738	! the model because there is no canopy (that's why we are in prescribe!) and
739	! so there is no lightstress. lstress was therefore set to one (see above).
740	! We prefered this redundancy in the code because it makes it clear that
741	! k_latosa is calculated in the same way in prescribe.f90 and growth_fun_all.f90
742	! +++CHECK+++
743	! How do we want to deal with waterstress
744	!!$ k_latosa(ipts,ivm) = k_latosa_min(ivm) + &
745	!!$ (wstress_fac(ipts,ivm) * lstress_fac(ipts,ivm) * &
746	!!$ (k_latosa_max(ivm)-k_latosa_min(ivm)))
747	!!$ k_latosa(ipts,ivm) = wstress_fac(ipts,ivm) * (k_latosa_min(ivm) + &
748	!!$ (lstress_fac(ipts,ivm) * &
749	!!$ (k_latosa_max(ivm)-k_latosa_min(ivm))))
750	k_latosa(ipts,ivm) = (k_latosa_adapt(ipts,ivm) + &
751	(lstress_fac(ipts,ivm) * &
752	(k_latosa_max(ivm)-k_latosa_min(ivm))))
753	! ++++++++++++
754
755
756	! Also k_latosa has been reported to be a function of CO2 concentration
757	! (Atwell et al. 2003, Tree Physiology, 23:13-21 and Pakati et al. 2000,
758	! Global Change Biology, 6:889-897). This effect is not accounted for in
759	! the current code
760
761	! Calculate conversion coefficient for sapwood area to leaf area
762	! (1) The scaling parameter between leaf and sapwood mass is derived from
763	! LA_ind = k_latosa * SA_ind, where LA_ind = leaf area of an individual, SA_ind is the
764	! sapwood area of an individual and k_latosa a pipe-model parameter
765	! (2) LA_ind = Cl * sla
766	! (3) Cs = SA_ind * height * wooddensity * tree_ff
767	! Substitute (2) and (3) in (1)
768	! Cl = Cs * k1 / (wooddensity * sla * tree_ff * height)
769	! Cl = Cs*KF/height, where KF is in (m)
770	! KF is passed to the allocation routine and it is saved in the restart file.
771	KF(ipts,ivm) = k_latosa(ipts,ivm) / ( sla(ivm) * pipe_density(ivm) * tree_ff(ivm))
772
773	! Initialize delta_KF to get the DO WHILE started
774	delta_KF = un
775
776	iloop=0
777	DO WHILE (delta_KF .GT. max_delta_KF)
778
779	! If there is a WHILE loop, there always needs to be a check on the number
780	! of loops to make sure we don't get stuck in an infinite loop. This
781	! number is completely arbitrary.
782	iloop=iloop+1
783	IF(iloop > 1000)THEN
784	WRITE(numout,*) 'Taking too long to converge the delta_KF loop in prescribe!'
785	WRITE(numout,*) 'iloop,delta_KF,max_delta_KF,ivm,ipts: ',&
786	iloop,delta_KF,max_delta_KF,ivm,ipts
787	CALL ipslerr_p (3,'stomate_prescribe',&
788	'Taking too long to converge the delta_KF','','')
789
790	ENDIF
791
792	!! 2.6 Create saplings for each height class
793	! Now we create the saplings for each class, based on the height
794	DO icir = 1,ncirc
795
796	! The assumption we make is that we plant trees of 2 to 3 years old rather than
797	! growing trees from seeds. The allometric relationship between height and
798	! diameter is derived from mature tree and likely unrealistic for saplings.
799	! The height of the saplings is prescribed and determines the reserves which are
800	! especially important for deciduous species which need to survive on their
801	! reserves for the first year (new phenology scheme requires annual mean values
802	! to get started)
803	dia_init = ( ave_tree_height(icir) / pipe_tune2(ivm) ) ** ( 1. / pipe_tune3(ivm) )
804	wood_init = ( ave_tree_height(icir) * pi / 4. * (dia_init) ** 2. ) * &
805	pipe_density(ivm) * tree_ff(ivm)
806
807	! The woody biomass is contained in four components. Thus, wood_init = isapabove +
808	! isapbelow + iheartabove + iheartbelow. Given that isapbelow = isapbelow and
809	! iheartabove = iheartbelow = bm_sapl_heartabove*isapabove. If bm_sapl_heartbelow =
810	! bm_sapl_heartabove = 0.2, then isapabove = wood_init/2.4
811	bm_sapl(ivm,icir,isapabove,icarbon) = wood_init / (2. + bm_sapl_heartabove + bm_sapl_heartbelow)
812	bm_sapl(ivm,icir,isapbelow,icarbon) = bm_sapl(ivm,icir,isapabove,icarbon)
813	bm_sapl(ivm,icir,iheartabove,icarbon) = bm_sapl_heartabove * bm_sapl(ivm,icir,isapabove,icarbon)
814	bm_sapl(ivm,icir,iheartbelow,icarbon) = bm_sapl_heartbelow * bm_sapl(ivm,icir,isapbelow,icarbon)
815
816	! Use the allometric relationships to calculate initial leaf and root mass
817	bm_sapl(ivm,icir,ileaf,icarbon) = ( bm_sapl(ivm,icir,isapabove,icarbon) + &
818	bm_sapl(ivm,icir,isapbelow,icarbon) ) * KF(ipts,ivm) / ave_tree_height(icir)
819	!+++CHECK+++
820	!How do we want to deal with water stress? wstress is accounted for through c0
821	!!$ bm_sapl(ivm,icir,iroot,icarbon) = bm_sapl(ivm,icir,ileaf,icarbon) / ( KF(ipts,ivm) * c0_alloc(ivm) )
822	bm_sapl(ivm,icir,iroot,icarbon) = bm_sapl(ivm,icir,ileaf,icarbon) / &
823	( KF(ipts,ivm) * c0_alloc(ipts,ivm) )
824	!+++++++++++
825
826	! Write initial values
827	IF (ld_presc) THEN
828	WRITE(numout,*) ' Circumference class: ',icir,' PFT type: ',ivm
829	WRITE(numout,*) ' root to sapwood tradeoff p :', c0_alloc(ipts,ivm)
830	WRITE(numout,*) 'height_init, dia_init, wood_init, ', &
831	ave_tree_height(icir) , dia_init, wood_init
832	WRITE(numout,*) 'pipe_density, ',pipe_density(ivm)
833	ENDIF
834
835	!++++++ CHECK ++++++
836	! The carbohydrate reserves do not seem to be set before this line. This
837	! is a problem since it then uses an unitililized value. Therefore, I
838	! will initilize it.
839	! Should this really be zero? If so the code below is nonesensical and icarbres
840	! could be omitted. nonesensical code = 2 * (bm_sapl(ivm,icir,icarbres,icarbon) + &
841	! bm_sapl(ivm,icir,ileaf,icarbon) + bm_sapl(ivm,icir,iroot,icarbon))
842	bm_sapl(ivm,icir,icarbres,icarbon)=zero
843	!++++++++++++++++++
844
845	! Pools that are defined in the same way for trees and grasses
846	bm_sapl(ivm,icir,ifruit,icarbon) = zero
847
848	!+++CHECK+++
849	! There is an inconsistency in the calculation - most pools are in gN
850	! but leaves is in gC. The correction is proposed, that implies that
851	! the parameter labile_reserve will need to be tuned
852	!!$ bm_sapl(ivm,icir,ilabile,icarbon) = labile_to_total * &
853	!!$ (bm_sapl(ivm,icir,ileaf,icarbon) / cn_leaf_prescribed(ivm) + &
854	!!$ fcn_root(ivm) * bm_sapl(ivm,icir,iroot,icarbon) + fcn_wood(ivm) * &
855	!!$ (bm_sapl(ivm,icir,isapabove,icarbon) + bm_sapl(ivm,icir,isapbelow,icarbon) + &
856	!!$ bm_sapl(ivm,icir,icarbres,icarbon)))
857
858	bm_sapl(ivm,icir,ilabile,icarbon) = labile_to_total * (bm_sapl(ivm,icir,ileaf,icarbon) + &
859	fcn_root(ivm) * bm_sapl(ivm,icir,iroot,icarbon) + fcn_wood(ivm) * &
860	(bm_sapl(ivm,icir,isapabove,icarbon) + bm_sapl(ivm,icir,isapbelow,icarbon) + &
861	bm_sapl(ivm,icir,icarbres,icarbon)))
862	!+++++++++++
863
864	! Avoid deciduous PFTs to have leaves out at establishment
865	! Whether the saplings have leaves or don't have leaves the first year doesn't really matter
866	! Either the approach is correct in the northern hemisphere or in the southern hemisphere. Note
867	! that the resource limitation approach starts with the sapling having leaves.
868	! Anyhow, a spin-up is needed to avoid issues with the initial conditions
869	IF ( pheno_type(ivm) .NE. 1 ) THEN
870
871	! Not evergreen. Deciduous PFTs now need to survive an extra year before bud burst. To
872	! ensure survival there are several options: (a) either the height of the initial
873	! vegetation is increased (this results in more reserves) or (b) the reserves could be
874	! increased. The second option may result in numerical issues further down
875	! the code as the optimal reserve level is calculated from the other biomass pools.
876	! Also some initial tests showed that higher results simply resulted in more respiration.
877	! Use taller trees to start with.
878	bm_sapl(ivm,icir,icarbres,icarbon) = (bm_sapl(ivm,icir,icarbres,icarbon) + &
879	bm_sapl(ivm,icir,ileaf,icarbon) + bm_sapl(ivm,icir,iroot,icarbon))
880	bm_sapl(ivm,icir,ileaf,icarbon) = zero
881	bm_sapl(ivm,icir,iroot,icarbon) = zero
882
883	! When deciduous trees have no leaves they are senescent
884	senescence(ipts,ivm) = .TRUE.
885
886	ELSE
887
888	! Initilize carbohydrate reserves for evergreen PFTs
889	bm_sapl(ivm,icir,icarbres,icarbon) = zero
890
891	! Evergreen trees never go into senescence
892	senescence(ipts,ivm) = .FALSE.
893
894	ENDIF
895
896	IF (ld_presc) THEN
897	WRITE(numout,*) ' sapling biomass (gC):',icir,ivm,ipts
898	WRITE(numout,*) ' leaves: (::bm_sapl(ivm,icir,ileaf,icarbon))',&
899	bm_sapl(ivm,icir,ileaf,icarbon)
900	WRITE(numout,*) ' sap above ground: (::bm_sapl(ivm,icir,ispabove,icarbon)):',&
901	bm_sapl(ivm,icir,isapabove,icarbon)
902	WRITE(numout,*) ' sap below ground: (::bm_sapl(ivm,icir,isapbelow,icarbon))',&
903	bm_sapl(ivm,icir,isapbelow,icarbon)
904	WRITE(numout,*) ' heartwood above ground: (::bm_sapl(ivm,icir,iheartabove,icarbon))',&
905	bm_sapl(ivm,icir,iheartabove,icarbon)
906	WRITE(numout,*) ' heartwood below ground: (::bm_sapl(ivm,icir,iheartbelow,icarbon))',&
907	bm_sapl(ivm,icir,iheartbelow,icarbon)
908	WRITE(numout,*) ' roots: (::bm_sapl(ivm,icir,iroot,icarbon))',&
909	bm_sapl(ivm,icir,iroot,icarbon)
910	WRITE(numout,*) ' fruits: (::bm_sapl(ivm,icir,ifruit,icarbon))',&
911	bm_sapl(ivm,icir,ifruit,icarbon)
912	WRITE(numout,*) ' carbohydrate reserve: (::bm_sapl(ivm,icir,icarbres,icarbon))',&
913	bm_sapl(ivm,icir,icarbres,icarbon)
914	WRITE(numout,*) ' labile reserve: (::bm_sapl(ivm,icir,ilabile,icarbon))',&
915	bm_sapl(ivm,icir,ilabile,icarbon)
916	ENDIF
917
918	ENDDO
919
920	IF (ld_presc) THEN
921	WRITE(numout,*)'Initial distribution, method 2',ivm
922	WRITE(numout,*)'Average trees height (m): ',ave_tree_height(:)
923	WRITE(numout,*) 'lambda',lambda,'nmaxtrees(ivm)',nmaxtrees(ivm)
924	ENDIF
925
926	!! 2.7 Determine the biomass in each circumference class.
927	! I do this based on the biomass in each sapling and the number of trees in each
928	! circumference class...we need the biomass in an average tree
929	DO icir=1,ncirc
930
931	circ_class_biomass(ipts,ivm,icir,:,icarbon) = bm_sapl(ivm,icir,:,icarbon)
932
933	ENDDO
934
935	! Now the total biomass is just a sum over all of these
936	DO ipar = 1,nparts
937
938	biomass(ipts,ivm,ipar,icarbon)= &
939	SUM(circ_class_biomass(ipts,ivm,:,ipar,icarbon)*circ_class_n(ipts,ivm,:))
940
941	ENDDO
942
943	! The light stress should be calculated making use of Pgap so it accounts for LAI
944	! crown dimensions and tree distribution. However, this would be computationally
945	! expensive so we just use a first order estimate based on light attenuation model
946	! by Lambert-Beer. When LAI is low, a lot of light reaches the forest floor and so
947	! KF should increase to make use of the available light by growing leaves
948	lstress_fac = exp(-biomass(ipts,ivm,ileaf,icarbon) * sla(ivm) * 0.5)
949	! Causing large differences between first and second prescribe
950	delta_KF = ABS (KF(ipts,ivm) - ((k_latosa_adapt(ipts,ivm) + &
951	(lstress_fac(ipts,ivm) * (k_latosa_max(ivm)-k_latosa_min(ivm))))) / &
952	( sla(ivm) * tree_ff(ivm) * pipe_density(ivm) ))
953	KF(ipts,ivm) = ((k_latosa_adapt(ipts,ivm) + &
954	(lstress_fac(ipts,ivm) * &
955	(k_latosa_max(ivm)-k_latosa_min(ivm)))) / &
956	( sla(ivm) * tree_ff(ivm) * pipe_density(ivm) ))
957
958	IF(ld_presc)THEN
959	WRITE(numout,*) 'prescribe delta_KF, ', delta_KF
960	WRITE(numout,*) 'prescribe lstress, ', lstress_fac
961	ENDIF
962	END DO
963
964	IF (ld_presc) THEN
965	WRITE(numout,*)'Initial biomass distribution'
966	DO icir=1,ncirc
967	WRITE(numout,*) circ_class_biomass(ipts,ivm,icir,:,icarbon),circ_class_n(ipts,ivm,icir)
968	ENDDO
969	WRITE(numout,*)'End initial biomass distribution',SUM(circ_class_n(ipts,ivm,:))
970	WRITE(numout,*)"biomass(ipts,ivm,:,icarbon)",biomass(ipts,ivm,:,icarbon)
971	END IF
972
973	IF (ivm .EQ. test_pft .AND. ld_presc) THEN
974	WRITE(numout,*) 'Check prescribe'
975	WRITE(numout,*) 'stomate_prescribe::init circ_class_n ',&
976	ipts,ivm,circ_class_n(ipts,ivm,:)
977	DO ipar = 1,nparts
978	WRITE(numout,*) 'biomass vs circ_class_biomass (gC m-2), ',&
979	ipts,ivm,ipar,biomass(ipts,ivm,ipar,icarbon), &
980	SUM(circ_class_biomass(ipts,ivm,:,ipar,icarbon)*circ_class_n(ipts,ivm,:))
981	WRITE(numout,*) 'stomate_prescribe::init biomass ',&
982	ipts,ivm,circ_class_biomass(ipts,ivm,:,ipar,icarbon)
983	ENDDO
984	ENDIF
985
986	! PFT is not present
987	ELSE
988
989	! At the tree level
990	circ_class_n(ipts,ivm,:) = zero
991	circ_class_biomass(ipts,ivm,:,:,icarbon) = zero
992
993	! At the stand level
994	biomass(ipts,ivm,:,icarbon) = zero
995	ind(ipts,ivm) = zero
996
997	ENDIF
998
999	! Set leaf age classes, all leaves are current year leaves
1000	leaf_frac(ipts,ivm,:) = zero
1001	leaf_frac(ipts,ivm,1) = un
1002
1003	!+++CHECK+++
1004	! Set time since last beginning of growing season but only
1005	! for the first day of the whole simulation. When the model
1006	! is initialized when_growthinit is set to undef. In subsequent
1007	! time steps it should have a value. For trees without phenology
1008	! the growing season starts at the moment the PFT is prescribed
1009	IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
1010	when_growthinit(ipts,ivm) = 200
1011	ENDIF
1012	!+++++++++++
1013
1014	! Seasonal trees have no leaves at beginning
1015	! Saplings of evergreen trees have a leaf mass on day 1 and mass in the other components
1016	! saplings of deciduous trees have no leaves on day 1 but mass in the other components.
1017	IF ( pheno_model(ivm) .NE. 'none' ) THEN
1018
1019	! Add the carbon from the leaves to the reserve pool
1020	biomass(ipts,ivm,icarbres,icarbon) = biomass(ipts,ivm,icarbres,icarbon) + biomass(ipts,ivm,ileaf,icarbon)
1021	biomass(ipts,ivm,ileaf,icarbon) = zero
1022	leaf_frac(ipts,ivm,1) = zero
1023
1024	!+++CHECK+++
1025	! Set time since last beginning of growing season but only
1026	! for the first day of the whole simulation. When the model
1027	! is initialized when_growthinit is set to undef. In subsequent
1028	! time steps it should have a value. The phenology module
1029	! prevents leaf onset soon after senescence, by setting
1030	! ::when_growthinit to a value, leaf offset
1031	! will occur at the first opportunity
1032	!!$ when_growthinit(ipts,ivm) = large_value
1033	IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
1034	when_growthinit(ipts,ivm) = 200
1035	ENDIF
1036	!++++++++++++
1037
1038	! Redundant, flag has already been set. When there are no leaves, the tree is in senescence
1039	senescence(ipts,ivm) = .TRUE.
1040
1041	ENDIF ! pheno_model(ivm)
1042
1043	! The biomass to build the saplings is taken from the atmosphere, keep track of
1044	! amount to calculate the C-balance closure
1045	co2_to_bm(ipts,ivm) = co2_to_bm(ipts,ivm) + ( SUM(biomass(ipts,ivm,:,icarbon)) / dt )
1046
1047	ENDIF ! tree(ivm)
1048
1049	!! 2.8 Initialize grassy PFTs
1050	!! Use veget_max to check whether the PFT is present. If the PFT is present but it
1051	!! has no biomass, prescribe its biomass (gC m-{2}). It is assumed that at day 1
1052	!! grasses have all their biomass in the reserve pool. The criteria exclude crops.
1053	!! Crops are no longer prescribed but planted the day that begin_leaves is true.
1054
1055	!+++TEMP+++
1056	IF(ld_presc .AND. test_pft == ivm)THEN
1057	WRITE(numout,*) 'prescribe - total biomass, ',SUM(biomass(ipts,ivm,:,icarbon))
1058	ENDIF
1059	!++++++++++
1060
1061	IF ( ( .NOT. is_tree(ivm) ) .AND. &
1062	( natural(ivm) ) .AND. &
1063	( veget_max(ipts,ivm) .GT. min_stomate ) .AND. &
1064	( SUM( biomass(ipts,ivm,:,icarbon) ) .LE. min_stomate ) ) THEN
1065
1066	!+++TEMP+++
1067	IF(ld_presc .AND. test_pft == ivm)THEN
1068	WRITE(numout,*) 'We will prescribe a new vegetation, the old one died'
1069	ENDIF
1070	!++++++++++
1071
1072	! For grasses we assume that an individual grass is 1 m2 of grass. This is set
1073	! in nmaxtrees in pft_parameters.f90. It could be set to another value but
1074	! with the current code this should not have any meaning. The grassland
1075	! does not necessarily covers the whole 1 m2 so adjust for the canopy cover
1076	ind(ipts,ivm) = nmaxtrees(ivm) * ha_to_m2 * canopy_cover(ivm)
1077
1078	! now we generate the size of a single grass sapling...this was all taken from
1079	! stomate_data.f90...we do not deal with circumference classes for grasses
1080	! and crops, but we want to keep the arrays the same as for the trees so
1081	! we put the sapling information into the first circumference class
1082
1083	!+++CHECK+++
1084	! Calculate the sapwood to leaf mass in a similar way as has been done for trees.
1085	! For trees this approach had been justified by observations. For grasses such
1086	! justification is not supported by observations but we didn't try to find it.
1087	! Needs more work by someone interested in grasses. There might be a more elegant
1088	! solution making use of a well observed parameter.
1089	!!$ k_latosa(ipts,ivm) = k_latosa_min(ivm) + &
1090	!!$ (wstress_fac(ipts,ivm) * lstress_fac(ipts,ivm) * &
1091	!!$ (k_latosa_max(ivm)-k_latosa_min(ivm)))
1092	!!$ k_latosa(ipts,ivm) = wstress_fac(ipts,ivm) * (k_latosa_min(ivm) + &
1093	!!$ (lstress_fac(ipts,ivm) * &
1094	!!$ (k_latosa_max(ivm)-k_latosa_min(ivm))))
1095	k_latosa(ipts,ivm) = (k_latosa_adapt(ipts,ivm) + &
1096	(lstress_fac(ipts,ivm) * &
1097	(k_latosa_max(ivm)-k_latosa_min(ivm))))
1098
1099	! The mass of the structural carbon relates to the mass of the leaves through
1100	! a prescribed parameter ::k_latosa
1101	KF(ipts,ivm) = k_latosa(ipts,ivm)
1102	!+++++++++++
1103
1104	! Calculate leaf to root area
1105	LF(ipts,ivm) = c0_alloc(ipts,ivm) * KF(ipts,ivm)
1106
1107	!---TEMP---
1108	IF(ld_presc .AND. test_pft == ivm)THEN
1109	WRITE(numout,*) 'KF, ', k_latosa(ipts,ivm), KF(ipts,ivm)
1110	WRITE(numout,) 'LF, c0_alloc, ', c0_alloc(ipts,ivm) KF(ipts,ivm), &
1111	c0_alloc(ipts,ivm)
1112	ENDIF
1113	!----------
1114
1115	! initialize everything to make sure there are not random values floating around
1116	! for ncirc != 1
1117	bm_sapl(ivm,:,:,icarbon) = val_exp
1118
1119	! Similar as for trees, the initial height of the vegetation was defined
1120	bm_sapl(ivm,1,ileaf,icarbon) = height_init_min(ivm) / lai_to_height(ivm) / sla(ivm)
1121
1122	! Use allometric relationships to define the root mass based on leaf mass. Some
1123	! sapwood mass is needed to store the reserves. An arbitrairy fraction of 5% was
1124	! used
1125	bm_sapl(ivm,1,iroot,icarbon) = bm_sapl(ivm,1,ileaf,icarbon) / LF(ipts,ivm)
1126	bm_sapl(ivm,1,isapabove,icarbon) = bm_sapl(ivm,1,ileaf,icarbon) / KF(ipts,ivm)
1127
1128	! Some of the biomass components that exist for trees are undefined for grasses
1129	bm_sapl(ivm,1,isapbelow,icarbon) = zero
1130	bm_sapl(ivm,1,iheartabove,icarbon) = zero
1131	bm_sapl(ivm,1,iheartbelow,icarbon) = zero
1132	bm_sapl(ivm,1,ifruit,icarbon) = zero
1133	bm_sapl(ivm,1,icarbres,icarbon) = zero
1134
1135	! Pools that are defined in the same way for trees and grasses
1136	bm_sapl(ivm,1,ilabile,icarbon) = labile_to_total * (bm_sapl(ivm,1,ileaf,icarbon) + &
1137	fcn_root(ivm) * bm_sapl(ivm,1,iroot,icarbon) + fcn_wood(ivm) * &
1138	(bm_sapl(ivm,1,isapabove,icarbon) + bm_sapl(ivm,1,isapbelow,icarbon) + &
1139	bm_sapl(ivm,1,icarbres,icarbon)))
1140
1141	! Avoid deciduous PFTs to have leaves out at establishment
1142	! Whether the saplings have leaves or don't have leaves the first year doesn't really matter
1143	! Either the approach is correct in the northern hemisphere or in the southern hemisphere. Note
1144	! that the resource limitation approach starts with the sapling having leaves.
1145	! Anyhow, a spin-up is needed to avoid issues with the initial conditions
1146	IF ( pheno_type(ivm) .NE. 1 ) THEN
1147
1148	! Not evergreen. Deciduous PFTs now need to survive an extra year before bud burst. To
1149	! ensure survival there are several options: (a) either the height of the initial
1150	! vegetation is increased (this results in more reserves) or (b) the reserves could be
1151	! increased. The second option may result in result in numerical issues further down
1152	! the code as the optimal reserve level is calculated from the other biomass pools
1153	bm_sapl(ivm,1,icarbres,icarbon) = bm_sapl(ivm,1,icarbres,icarbon) + bm_sapl(ivm,1,ileaf,icarbon) + &
1154	bm_sapl(ivm,1,iroot,icarbon)
1155	bm_sapl(ivm,1,ileaf,icarbon) = zero
1156	bm_sapl(ivm,1,iroot,icarbon) = zero
1157
1158	! When there are no leaves, the crop/grass is in senescence
1159	senescence(ipts,ivm) = .TRUE.
1160
1161	ELSE
1162
1163	! Initilize carbohydrate reserves for evergreen PFTs
1164	bm_sapl(ivm,1,icarbres,icarbon) = zero
1165
1166	! Evergreen plants never go into senescence
1167	senescence(ipts,ivm) = .FALSE.
1168
1169	ENDIF
1170
1171	IF (ld_presc) THEN
1172	WRITE(numout,*) ' sapling biomass (gC):',1,ivm,ipts
1173	WRITE(numout,*) ' leaves: (::bm_sapl(ivm,1,ileaf,icarbon))',&
1174	bm_sapl(ivm,1,ileaf,icarbon)
1175	WRITE(numout,*) ' sap above ground: (::bm_sapl(ivm,1,ispabove,icarbon)):',&
1176	bm_sapl(ivm,1,isapabove,icarbon)
1177	WRITE(numout,*) ' sap below ground: (::bm_sapl(ivm,1,isapbelow,icarbon))',&
1178	bm_sapl(ivm,1,isapbelow,icarbon)
1179	WRITE(numout,*) ' heartwood above ground: (::bm_sapl(ivm,1,iheartabove,icarbon))',&
1180	bm_sapl(ivm,1,iheartabove,icarbon)
1181	WRITE(numout,*) ' heartwood below ground: (::bm_sapl(ivm,1,iheartbelow,icarbon))',&
1182	bm_sapl(ivm,1,iheartbelow,icarbon)
1183	WRITE(numout,*) ' roots: (::bm_sapl(ivm,1,iroot,icarbon))',&
1184	bm_sapl(ivm,1,iroot,icarbon)
1185	WRITE(numout,*) ' fruits: (::bm_sapl(ivm,1,ifruit,icarbon))',&
1186	bm_sapl(ivm,1,ifruit,icarbon)
1187	WRITE(numout,*) ' carbohydrate reserve: (::bm_sapl(ivm,1,icarbres,icarbon))',&
1188	bm_sapl(ivm,1,icarbres,icarbon)
1189	WRITE(numout,*) ' labile reserve: (::bm_sapl(ivm,1,ilabile,icarbon))',&
1190	bm_sapl(ivm,1,ilabile,icarbon)
1191	ENDIF
1192
1193	! Write initial values
1194	IF (ld_presc) THEN
1195	WRITE(numout,*) ' root to sapwood tradeoff (LF) : ', c0_alloc(ipts,ivm)
1196	WRITE(numout,*) ' grass sapling biomass: ',bm_sapl(ivm,1,:,icarbon)
1197	ENDIF
1198
1199	! Initial biomass (g C m-2)
1200	biomass(ipts,ivm,:,icarbon) = bm_sapl(ivm,1,:,icarbon) * ind(ipts,ivm)
1201	IF (ld_presc) THEN
1202	WRITE(numout,*) 'bm_grass, prescribe, ', biomass(ipts,ivm,:,icarbon)
1203	ENDIF
1204
1205	! Synchronize biomass and circ_class_biomass (gC tree-1). This
1206	! allows to more easily check the mass balance closure
1207	circ_class_biomass(ipts,ivm,1,:,icarbon) = bm_sapl(ivm,1,:,icarbon)
1208	circ_class_n(ipts,ivm,1) = ind(ipts,ivm)
1209
1210	! Set leaf age classes -> all leaves will be current year leaves
1211	leaf_frac(ipts,ivm,:) = zero
1212	leaf_frac(ipts,ivm,1) = un
1213
1214	! Set time since last beginning of growing season but only
1215	! for the first day of the whole simulation. When the model
1216	! is initialized when_growthinit is set to undef. In subsequent
1217	! time steps it should have a value.
1218	!!$ when_growthinit(ipts,ivm) = large_value
1219	IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
1220	when_growthinit(ipts,ivm) = 200
1221	ENDIF
1222
1223	! The biomass to build the saplings is taken from the atmosphere, keep track of
1224	! amount to calculate the C-balance closure
1225	co2_to_bm(ipts,ivm) = co2_to_bm(ipts,ivm) + ( SUM(biomass(ipts,ivm,:,icarbon)) / dt )
1226
1227	ELSEIF (.NOT. natural(ivm)) THEN
1228
1229	! Initialize croplands - else leaves_begin will never become true
1230	! Set time since last beginning of growing season but only
1231	! for the first day of the whole simulation. When the model
1232	! is initialized when_growthinit is set to undef. In subsequent
1233	! time steps it should have a value.
1234	IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
1235	when_growthinit(ipts,ivm) = 200
1236	ENDIF
1237
1238	ENDIF ! .NOT. tree(ivm)
1239
1240	!! 2.3 Declare PFT present
1241	!! Now that the PFT has biomass it should be declared 'present'
1242	!! everywhere in that grid box. Assign some additional properties
1243	PFTpresent(ipts,ivm) = .TRUE.
1244	everywhere(ipts,ivm) = un
1245	age(ipts,ivm) = zero
1246	npp_longterm(ipts,ivm) = npp_longterm_init
1247	lm_lastyearmax(ipts,ivm) = zero
1248
1249	ENDIF ! not control%ok_dgvm or agricultural
1250
1251	ENDIF ! IF (veget_max .GT. zero .AND. ind .EQ. zero)
1252
1253	ENDDO ! loop over pixels
1254
1255	ENDDO ! loop over PFTs
1256
1257	!! 3. Scynchronize biomass and circ_class_biomass
1258
1259	! For trees both biomass and circ_class_biomass are recorded. Where biomass
1260	! contains the stand level biomass of the different components (gC m-2),
1261	! circ_class_biomass contains the biomass of a model tree in each circumference
1262	! class (gC tree-1). When multiplied with the number of trees in each circumference
1263	! class (::circ_class_n), biomass and circ_class_biomass should be identical.
1264	! This has been checked in all routines where biomass is calculated and at the
1265	! daily level consistency between both variables is maintained. However, due
1266	! to precision issues and the fact tha biomass and circ_class_biomass are
1267	! cumulative variables, the difference between both variables accumulates to
1268	! 0.00001 at the annual level. To avoid propagation of the precision error
1269	! both variables are synchronized daily.
1270
1271	! Initialze the variables for the synchronisation
1272	sync_ind(:,:) = zero
1273	sync_biomass(:,:,:,:) = zero
1274
1275	! Synchronize
1276	DO ivm = 1,nvm
1277
1278	IF (is_tree(ivm)) THEN
1279
1280	DO icir = 1,ncirc
1281
1282	! Temporary stand density
1283	sync_ind(:,ivm) = sync_ind(:,ivm) + circ_class_n(:,ivm,icir)
1284
1285	DO iele = 1,nelements
1286
1287	! Temporary biomass
1288	DO ipar = 1,nparts
1289
1290	sync_biomass(:,ivm,ipar,iele) = sync_biomass(:,ivm,ipar,iele) + &
1291	circ_class_n(:,ivm,icir) * circ_class_biomass(:,ivm,icir,ipar,iele)
1292
1293	ENDDO ! nparts
1294
1295	ENDDO ! nelements
1296
1297	ENDDO ! ncirc
1298
1299	! Synchronize stand density
1300	DO ipts = 1,npts
1301
1302	DO iele = 1,nelements
1303
1304	IF (ivm .EQ. test_pft .AND. ABS(sync_ind(ipts,ivm) - ind(ipts,ivm)) .GT. min_stomate ) THEN
1305
1306	IF (ld_presc) THEN
1307	WRITE(numout,*) 'Divergence of density between ::ind and ::circ_class_n'
1308	WRITE(numout,*) 'This problem has been identified in stomate_prescribe but could'
1309	WRITE(numout,*) 'be caused in any routine calculating stand density'
1310	ENDIF
1311
1312	ELSE
1313
1314	ind(ipts,ivm) = sync_ind(ipts,ivm)
1315
1316	ENDIF
1317
1318	! Synchronize biomass
1319	IF (ld_presc .AND. ivm .EQ. test_pft)THEN
1320	IF( SUM(biomass(ipts,ivm,:,iele)) .NE. zero)THEN
1321	IF( SUM( ABS(sync_biomass(ipts,ivm,:,iele) - biomass(ipts,ivm,:,iele) ) ) / &
1322	SUM(biomass(ipts,ivm,:,iele)) .GT. 0.0001 ) THEN
1323
1324	WRITE(numout,*) 'Divergence of biomass between ::biomass and ::circ_class_biomass'
1325	WRITE(numout,*) 'This problem has been identified in stomate_prescribe but could'
1326	WRITE(numout,*) 'caused in any routine calculating biomass components'
1327	WRITE(numout,*) '% divergence', ivm, &
1328	SUM( ABS(sync_biomass(ipts,ivm,:,iele) - biomass(ipts,ivm,:,iele) ) ) / &
1329	SUM(biomass(ipts,ivm,:,iele)), &
1330	SUM(sync_biomass(ipts,ivm,:,iele)), SUM(biomass(ipts,ivm,:,iele))
1331	ENDIF
1332	ENDIF
1333	ENDIF
1334
1335	biomass(ipts,ivm,:,iele) = sync_biomass(ipts,ivm,:,iele)
1336
1337
1338	ENDDO ! nelements
1339
1340	ENDDO ! npts
1341
1342	ENDIF ! is_tree
1343
1344	ENDDO ! # PFT's
1345
1346	!! 4. Calculate components of the mass balance
1347
1348	!! 4.1 Calculate final biomass
1349	pool_end(:,:,:) = zero
1350	DO ipar = 1,nparts
1351	DO iele = 1,nelements
1352	pool_end(:,:,iele) = pool_end(:,:,iele) + &
1353	(biomass(:,:,ipar,iele) * veget_max(:,:))
1354	ENDDO
1355	ENDDO
1356
1357	!! 4.2 Calculate mass balance
1358	check_intern(:,:,iatm2land,icarbon) = check_intern(:,:,iatm2land,icarbon) + &
1359	co2_to_bm(:,:) * veget_max(:,:) * dt
1360	check_intern(:,:,iland2atm,icarbon) = -un * zero
1361	check_intern(:,:,ilat2out,icarbon) = zero
1362	check_intern(:,:,ilat2in,icarbon) = -un * zero
1363	check_intern(:,:,ipoolchange,icarbon) = -un * (pool_end(:,:,icarbon) - pool_start(:,:,icarbon))
1364	closure_intern = zero
1365	DO imbc = 1,nmbcomp
1366	closure_intern(:,:,icarbon) = closure_intern(:,:,icarbon) + check_intern(:,:,imbc,icarbon)
1367	ENDDO
1368
1369	!! 4.3 Write outcome of checks
1370	DO ipts=1,npts
1371	DO ivm=1,nvm
1372
1373	! Mass balance closure
1374	IF(ABS(closure_intern(ipts,ivm,icarbon)) .LE. min_stomate)THEN
1375	IF (ld_massbal) WRITE(numout,*) 'Mass balance closure in prescribe_prognostic'
1376	ELSE
1377	WRITE(numout,*) 'Error: mass balance is not closed in prescribe_prognostic'
1378	WRITE(numout,*) ' ipts,ivm; ', ipts,ivm
1379	WRITE(numout,*) ' Difference is, ', closure_intern(ipts,ivm,icarbon)
1380	WRITE(numout,*) ' pool_end,pool_start: ', pool_end(ipts,ivm,icarbon), pool_start(ipts,ivm,icarbon)
1381	WRITE(numout,*) ' check_intern,co2_to_bm,pool_end,veget_max: ', &
1382	check_intern(ipts,ivm,iatm2land,icarbon),co2_to_bm(ipts,ivm), veget_max(ipts,ivm)
1383	IF(ld_stop)THEN
1384	CALL ipslerr_p (3,'prescribe_prognostic', 'Mass balance error.','','')
1385	ENDIF
1386	ENDIF
1387	ENDDO
1388	ENDDO
1389
1390	END SUBROUTINE prescribe_prognostic
1391
1392	END MODULE stomate_prescribe

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

Download in other formats: