Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-261/src_stomate/stomate_prescribe.f90 @ 8692

Last change on this file since 8692 was 4333, checked in by josefine.ghattas, 8 years ago
Merge with trunk for revision r3623 - r3977 Some corrections for gfortran
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 55.9 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, crown area.
10	!!
11	!!\n DESCRIPTION: None
12	!!
13	!! RECENT CHANGE(S): None
14	!!
15	!! REFERENCE(S) :
16	!!
17	!! SVN :
18	!! $HeadURL$
19	!! $Date$
20	!! $Revision$
21	!! \n
22	!_ ================================================================================================================================
23
24	MODULE stomate_prescribe
25
26	! modules used:
27
28	USE ioipsl_para
29	USE stomate_data
30	USE pft_parameters
31	USE constantes
32	USE function_library, ONLY: calculate_c0_alloc, biomass_to_lai
33
34	IMPLICIT NONE
35
36	! private & public routines
37
38	PRIVATE
39	PUBLIC prescribe,prescribe_clear
40
41	! first call
42	LOGICAL, SAVE :: firstcall_prescribe = .TRUE.
43	!$OMP THREADPRIVATE(firstcall_prescribe)
44
45	CONTAINS
46
47	! =================================================================================================================================
48	!! SUBROUTINE : prescribe_clear
49	!!
50	!>\BRIEF : Set the firstcall_prescribe flag back to .TRUE. to prepare for the next simulation.
51	!_=================================================================================================================================
52
53	SUBROUTINE prescribe_clear
54	firstcall_prescribe=.TRUE.
55	END SUBROUTINE prescribe_clear
56
57
58	!! ================================================================================================================================
59	!! SUBROUTINE : prescribe
60	!!
61	!>\BRIEF Works only with static vegetation and agricultural PFT. Initialize biomass,
62	!! density, presence in the first call and update them in the following.
63	!!
64	!! DESCRIPTION (functional, design, flags): \n
65	!! This module works only with static vegetation and agricultural PFT.
66	!! In the first call, initialize density of individuals, biomass,
67	!! and leaf age distribution to some reasonable value. In the following calls,
68	!! these variables are updated.
69	!!
70	!! To fulfill these purposes, pipe model are used:
71	!! \latexonly
72	!! \input{prescribe1.tex}
73	!! \input{prescribe2.tex}
74	!! \endlatexonly
75	!!
76	!! RECENT CHANGE(S): None
77	!!
78	!! MAIN OUTPUT VARIABLES(S): ::ind, ::cn_ind, ::leaf_frac
79	!!
80	!! REFERENCES :
81	!! - Krinner, G., N. Viovy, et al. (2005). "A dynamic global vegetation model
82	!! for studies of the coupled atmosphere-biosphere system." Global
83	!! Biogeochemical Cycles 19: GB1015, doi:1010.1029/2003GB002199.
84	!! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics,
85	!! plant geography and terrestrial carbon cycling in the LPJ dynamic
86	!! global vegetation model, Global Change Biology, 9, 161-185.
87	!! - McDowell, N., Barnard, H., Bond, B.J., Hinckley, T., Hubbard, R.M., Ishii,
88	!! H., KÃ¶stner, B., Magnani, F. Marshall, J.D., Meinzer, F.C., Phillips, N.,
89	!! Ryan, M.G., Whitehead D. 2002. The relationship between tree height and leaf
90	!! area: sapwood area ratio. Oecologia, 132:12â20.
91	!! - Novick, K., Oren, R., Stoy, P., Juang, F.-Y., Siqueira, M., Katul, G. 2009.
92	!! The relationship between reference canopy conductance and simplified hydraulic
93	!! architecture. Advances in water resources 32, 809-819.
94	!!
95	!! FLOWCHART : None
96	!! \n
97	!_ ================================================================================================================================
98
99	SUBROUTINE prescribe (npts, veget_max, veget, dt, PFTpresent, &
100	everywhere, when_growthinit, biomass, leaf_frac, &
101	ind, co2_to_bm, &
102	KF, &
103	senescence, age, npp_longterm, lm_lastyearmax, k_latosa_adapt)
104
105	!! 0. Parameters and variables declaration
106
107	!! 0.1 Input variables
108
109	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
110	REAL(r_std), INTENT(in) :: dt !! time step (dt_days)
111	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT
112	!! (LAI -> infinity) on ground. May sum to
113	!! less than unity if the pixel has
114	!! nobio area. (unitless; 0-1)
115	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget !! Fraction of vegetation type including
116	!! non-biological fraction (unitless)
117
118	!! 0.2 Output variables
119
120
121	!! 0.3 Modified variables
122
123	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! PFT present (0 or 1)
124	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Flag for setting senescence stage (only
125	!! for deciduous trees)
126	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
127	!! very localized (after its introduction) (?)
128	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! how many days ago was the beginning of
129	!! the growing season (days)
130	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Density of individuals at the stand level
131	!! @tex $(m^{-2})$ @endtex
132	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "long term" net primary productivity
133	!! @tex ($gC m^{-2} year^{-1}$) @endtex
134	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
135	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
136	!! @tex ($gC m^{-2}$) @endtex
137	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Stand level biomass
138	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C
139	!! to create the seedlings
140	!! @tex (gC.m^{-2}dt^{-1})$ @endtex
141	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age
142	!! class (unitless;0-1)
143	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: KF !! Scaling factor to convert sapwood mass
144	!! into leaf mass (m) - this variable is
145	!! passed to other routines
146	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: k_latosa_adapt !! Leaf to sapwood area adapted for long
147	!! term water stress (m)
148
149
150	!! 0.4 Local variables
151
152	REAL(r_std), DIMENSION(nvm) :: c0_alloc !! Root to sapwood tradeoff parameter
153	INTEGER(i_std) :: ipts, ivm, ipar, k !! index (unitless)
154	INTEGER(i_std) :: icir, iele, imbc !! index (unitless)
155	INTEGER(i_std) :: deb,fin, imaxt !! index (unitless)
156	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements) :: sync_biomass !! Temporary stand level biomass
157	!! @tex $(gC.m^{-2})$ @endtex
158	REAL(r_std), DIMENSION(npts,nvm) :: sync_ind !! Temporary density of individuals at the
159	!! stand level @tex $(m^{-2})$ @endtex
160	REAL(r_std), DIMENSION(npts,nvm) :: k_latosa !! Height dependent base value to calculate
161	!! KF (-)
162	! REAL(r_std), DIMENSION(nvm,nparts,nelements) :: bm_sapl !! Sapling biomass for the functional
163	!! allocation with a dimension for the
164	!! circumference classes
165	!! @tex $(gC.ind^{-1})$ @endtex
166	REAL(r_std), DIMENSION(npts,nvm) :: LF !! Scaling factor to convert sapwood mass
167	!! into root mass (unitless)
168	REAL(r_std), DIMENSION(npts,nvm) :: lstress_fac !! Light stress factor, based on total
169	!! transmitted light (unitless, 0-1)
170	REAL(r_std), DIMENSION(nparts) :: bm_init !! Biomass needed to initiate the next
171	!! planting @tex $(gC m^{-2})$ @endtex
172	REAL(r_std) :: nb_trees_i !! Number of trees in each twentith
173	!! circumference quantile of the
174	!! distribution (ind)
175	REAL(r_std) :: excedent !! Number of trees after truncation to be
176	!! reallocated to smaller quantiles of the
177	!! distribution (ind)
178	REAL(r_std) :: ave_tree_height !! The height of the ideal tree in each
179	!! circumference class of the
180	!! distribution (ind)...not saved since
181	!! it should only be used for prescribing
182	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
183	!! mass balance chech for this routine
184	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
185	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
186	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
187	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start and end pool of this routine
188	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
189	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! Start and end pool of this routine
190	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
191	REAL(r_std) :: delta_KF !! Difference between new and old estimate
192	!! of KF while iterating
193	REAL(r_std) :: min_height_init
194	REAL(r_std) :: max_height_init
195	REAL(r_std) :: sapwood_density
196	REAL(r_std) :: dia_init
197	REAL(r_std) :: wood_init
198	INTEGER(i_std) :: iloop
199
200	!++++ temp ++++
201	!!$ REAL(r_std), DIMENSION(npts,nvm,ncirc) :: height,dia,cn !! Tree height calculated from allometric
202	!!$ !! relationships (m)
203	REAL(r_std) :: cn_leaf,cn_root,cn_wood
204	!! CN ratio of leaves, root and wood pool
205	!! (gC/gN)
206	!_ ================================================================================================================================
207
208	!! 1. Initialize biomass at first call
209
210	!! 1.1.1 Initialize check for mass balance closure
211	! The mass balance is calculated at the end of this routine
212	! in section 4
213	! Initial biomass pool
214	pool_start(:,:,:) = zero
215	DO ipar = 1,nparts
216	DO iele = 1,nelements
217	pool_start(:,:,iele) = pool_start(:,:,iele) + &
218	(biomass(:,:,ipar,iele) * veget_max(:,:))
219	ENDDO
220	ENDDO
221
222	! co2_to_bm is has as intent inout, the variable accumulates
223	! carbon over the course of a day. Use the difference between
224	! start and the end of this routine
225	check_intern(:,:,iatm2land,icarbon) = - un * co2_to_bm(:,:) * veget_max(:,:) * dt
226
227	! Prescribe was taken out of the first call loop. Because of the firstcall here,
228	! nothing was done when the vegetation was killed in lpj_gap and the vegetation
229	! thus stayed dead forever. We want it to regrow. So instead, let's loop over
230	! every point and PFT. If there is supposed to be vegetation here (according to
231	! veget_max) but there isn't (according to ind), then we grow it.
232	DO ivm = 2,nvm ! Loop over PFTs
233
234	DO ipts = 1, npts ! Loop over pixels
235
236
237	! We are supposed to have vegetation, but we don't have any: prescribe some.
238	IF((veget_max(ipts,ivm) .GT. min_stomate) .AND. (ind(ipts,ivm) .LE. min_stomate))THEN
239
240	! Initilaize tree level variables
241	ind(ipts,ivm) = zero
242	biomass(ipts,ivm,:,:) = zero
243
244	! Assume that there is no memory for k_latosa between
245	! different generations (this is probably true when trees
246	! are planted but seems less likely for natural regeneration
247	! Allowing for memory has caused problems with the LAI from
248	! the second generation onwards
249	k_latosa_adapt(ipts,ivm) = k_latosa_min(ivm)
250
251	! Write output
252	WRITE(numout,*) 'Prescribe (stomate_prescribe.f90):'
253	WRITE(numout,*) ' > Imposing initial biomass for prescribed trees, '// &
254	'initial reserve mass for prescribed grasses.'
255	WRITE(numout,*) ' > Declaring prescribed PFTs present.'
256	WRITE(numout,*) ' > Grid point: ',ipts,' PFT Type: ',ivm
257
258	!! 2. Calculate the vegetation characteristics of a newly established vegetation
259
260	!! 2.1 Stress factors
261
262	! Note that light and water stress have an opposite effect on KF. Waterstress
263	! will decrease KF because more C should be allocated to the roots. Light
264	! stress will increase KF because more C should be allocated to the leaves.
265
266	! We will need the c0_alloc, it only depends on PFT and the effective
267	! longiveties
268	c0_alloc(ivm) = calculate_c0_alloc(ipts, ivm, tau_root(ivm), &
269	tau_sap(ivm))
270
271	! Lightstress varies from 0 to 1 and is calculated from the canopy structure (veget)
272	! Given that there is no vegetation at this point, the lightstress cannot be calculated
273	! and is set to 1. However as soon as we grow a canopy it will experience light stress
274	! and so KF should be adjusted. If this is not done, we can't prescribe tall vegetation
275	! which is a useful feature to speed up optimisation and testing. We will have iterate
276	! over light stress and KF to prescribe vegetation that is in balance with its
277	! light environment.
278	lstress_fac(ipts,ivm) = un
279
280	! Initial vegetation has to be prescribed only when the vegetation is static
281	IF ( ( .NOT. ok_dgvm ) .AND. &
282	( veget_max(ipts,ivm) .GT. min_stomate ) ) THEN
283
284	!! 2.2 Initialize woody PFT's
285	! Use veget_max to check whether the PFT is present. If the PFT is present but it has no
286	! biomass, prescribe its biomass (gC m-{2}).
287	IF ( is_tree(ivm) .AND. &
288	( veget_max(ipts,ivm) .GT. min_stomate ) .AND. &
289	( SUM( biomass(ipts,ivm,:,icarbon) ) .LE. min_stomate ) ) THEN
290
291	! PFT is present so it needs to be initialized
292	IF (veget_max(ipts,ivm) .GT. min_stomate) THEN
293
294	! The forest stand starts with the number of individuals as prescribed in constantes_mtc.f90
295	! This number is defined per hectare whereas these calculations are per m2
296	ind(ipts,ivm) = nmaxtrees(ivm) * ha_to_m2
297
298	!! 2.3 Initialize height distribution
299	! Initializing height distribution...I've taken this from the circumference initilization of Valentin
300	! we want ncirc bins between height_init_min and height_init_max
301	min_height_init = height_init_min(ivm)
302	max_height_init = height_init_max(ivm)
303
304	!! 2.4 Calculate height distribution
305
306	!! 2.4.1 Height distribution of a high stand
307	! Deterministic initial distribution following a truncated exponential law
308	! if they are coppices, we handle this differently
309
310
311	! this array is going to be used to create the sapling for this class so that we can create
312	! the total amount of biomass
313	! in this class...after this we will always calculate the height/circ/diam/whatever
314	! from the biomass in the class
315	ave_tree_height=(0.5_r_std)*&
316	(max_height_init+min_height_init)
317
318	! I am changing this so that our lowest class is min_height_init and our biggest class
319	! is max_height_init. With this change, we can use the same code the redistribute
320	! the trees if the biomass in one of our height classes is equal to zero later on due to mortality or
321	! thinning.
322	!+++ This caused some strange behavior, so change it back for now.
323	! ave_tree_height(icir)=min_height_init+REAL(icir-1,r_std)/REAL(ncirc-1,r_std)*&
324	! (max_height_init-min_height_init)
325
326
327	IF (printlev>=4) THEN
328	WRITE(numout,*)"min initial height ",min_height_init,'max initial height ',max_height_init
329	ENDIF
330
331	! ENDIF ! check FM type
332
333	!! 2.4.2 Circumference distribution of a coppice
334	! I'm not see why this should be any different than a standard prescribe.
335	! Ideally, the plantings should be stems 20-25 cm long stuck into the soil,
336	! but the allocation scheme will have the same issues with this that it does
337	! with planting a small sapling, i.e. it won't like it. The stems do not
338	! follow standard allocation rules. Then again, neither do coppices.
339
340	!!$ IF (forest_managed(ipts,ivm) == ifm_src) THEN
341	!!$ CALL ipslerr(3,'stomate_prescribe.f90','Problem with SRC','','')
342	!!$ ENDIF
343
344	!! 2.5 Allocation factors
345	! Sapwood to root ratio
346	! Following Magnani et al. 2000 "In order to decreases hydraulic resistance,
347	! the investment of carbon in fine roots or sapwood yields to the plant very
348	! different returns, both because of different hydraulic conductivities and
349	! because of the strong impact of plant height on shoot resistance. On the
350	! other hand, fine roots and sapwood have markedly different longevities and
351	! the cost of production, discounted for turnover, will differ accordingly.
352	! Optimal growth under hydraulic constraints requires that the ratio of marginal
353	! hydraulic returns to marginal annual cost for carbon investment in either roots
354	! or sapwood be the same (Bloom, Chapin & Mooney 1985; Case & Fair 1989). This
355	! is formalized in equation (13) and further derived to obtain equation (17)
356	! in Magnani et al 2000. The latter is implemented here.
357	! Pipe_density is given in gC/m-3, convert to kg/m-3. And apply equation (17)
358	! in Magnani et al 2000. Note that c0_alloc was calculated at the start of this
359	! routine. The calculation itself is done in function_library
360
361	! Calculate leaf area to sapwood area
362	! To be consistent with the hydraulic limitations and pipe theory,
363	! k_latosa is calculated from equation (18) in Magnani et al.
364	! To do so, total hydraulic resistance and tree height need to known. This
365	! poses a problem as the resistance depends on the leaf area and the leaf
366	! area on the resistance. There is no independent equation and equations 12
367	! and 18 depend on each other and substitution would be circular. Hence
368	! prescribed k_latosa values were obtained from observational records
369	! and are given in mtc_parameters.f90.
370
371	! The relationship between height and k_latosa as reported in McDowell
372	! et al 2002 and Novick et al 2009 is implemented to adjust k_latosa for
373	! the height of the stand. This did NOT result in a realistic model behavior
374	!!$ k_latosa(ipts,ivm) = wstress_fac(ipts,ivm) * &
375	!!$ (k_latosa_max(ivm) - (latosa_height(ivm) * &
376	!!$ (SUM( nb_trees_i(:) * ave_tree_height(:) ) / SUM( nb_trees_i(:) ))))
377
378	! Alternatively, k_latosa is also reported to be a function of diameter
379	! (i.e. stand thinning, Simonin et al 2006, Tree Physiology, 26:493-503).
380	! Here the relationship with thinning was interpreted as a realtionship with
381	! light stress. Note that light stress cannot be calculated at this time in
382	! the model because there is no canopy (that's why we are in prescribe!) and
383	! so there is no lightstress. lstress was therefore set to one (see above).
384	! We prefered this redundancy in the code because it makes it clear that
385	! k_latosa is calculated in the same way in prescribe.f90 and growth_fun_all.f90
386	! +++CHECK+++
387	! How do we want to deal with waterstress
388	!!$ k_latosa(ipts,ivm) = k_latosa_min(ivm) + &
389	!!$ (wstress_fac(ipts,ivm) * lstress_fac(ipts,ivm) * &
390	!!$ (k_latosa_max(ivm)-k_latosa_min(ivm)))
391	!!$ k_latosa(ipts,ivm) = wstress_fac(ipts,ivm) * (k_latosa_min(ivm) + &
392	!!$ (lstress_fac(ipts,ivm) * &
393	!!$ (k_latosa_max(ivm)-k_latosa_min(ivm))))
394	k_latosa(ipts,ivm) = (k_latosa_adapt(ipts,ivm) + &
395	(lstress_fac(ipts,ivm) * &
396	(k_latosa_max(ivm)-k_latosa_min(ivm))))
397	! ++++++++++++
398
399
400	! Also k_latosa has been reported to be a function of CO2 concentration
401	! (Atwell et al. 2003, Tree Physiology, 23:13-21 and Pakati et al. 2000,
402	! Global Change Biology, 6:889-897). This effect is not accounted for in
403	! the current code
404
405	! Calculate conversion coefficient for sapwood area to leaf area
406	! (1) The scaling parameter between leaf and sapwood mass is derived from
407	! LA_ind = k_latosa * SA_ind, where LA_ind = leaf area of an individual, SA_ind is the
408	! sapwood area of an individual and k_latosa a pipe-model parameter
409	! (2) LA_ind = Cl * sla
410	! (3) Cs = SA_ind * height * wooddensity * tree_ff
411	! Substitute (2) and (3) in (1)
412	! Cl = Cs * k1 / (wooddensity * sla * tree_ff * height)
413	! Cl = Cs*KF/height, where KF is in (m)
414	! KF is passed to the allocation routine and it is saved in the restart file.
415	KF(ipts,ivm) = k_latosa(ipts,ivm) / ( sla(ivm) * pipe_density(ivm) * tree_ff(ivm))
416
417	! Initialize delta_KF to get the DO WHILE started
418	delta_KF = un
419
420	iloop=0
421	DO WHILE (delta_KF .GT. max_delta_KF)
422
423	! If there is a WHILE loop, there always needs to be a check on the number
424	! of loops to make sure we don't get stuck in an infinite loop. This
425	! number is completely arbitrary.
426	iloop=iloop+1
427	IF(iloop > 1000)THEN
428	WRITE(numout,*) 'Taking too long to converge the delta_KF loop in prescribe!'
429	WRITE(numout,*) 'iloop,delta_KF,max_delta_KF,ivm,ipts: ',&
430	iloop,delta_KF,max_delta_KF,ivm,ipts
431	CALL ipslerr_p (3,'stomate_prescribe',&
432	'Taking too long to converge the delta_KF','','')
433
434	ENDIF
435
436	!! 2.6 Create saplings for each height class
437	! Now we create the saplings for each class, based on the height
438
439	! The assumption we make is that we plant trees of 2 to 3 years old rather than
440	! growing trees from seeds. The allometric relationship between height and
441	! diameter is derived from mature tree and likely unrealistic for saplings.
442	! The height of the saplings is prescribed and determines the reserves which are
443	! especially important for deciduous species which need to survive on their
444	! reserves for the first year (new phenology scheme requires annual mean values
445	! to get started)
446	dia_init = ( ave_tree_height / pipe_tune2(ivm) ) ** ( 1. / pipe_tune3(ivm) )
447	wood_init = ( ave_tree_height * pi / 4. * (dia_init) ** 2. ) * &
448	pipe_density(ivm) * tree_ff(ivm)
449
450	! The woody biomass is contained in four components. Thus, wood_init = isapabove +
451	! isapbelow + iheartabove + iheartbelow. Given that isapbelow = isapbelow and
452	! iheartabove = iheartbelow = bm_sapl_heartabove*isapabove. If bm_sapl_heartbelow =
453	! bm_sapl_heartabove = 0.2, then isapabove = wood_init/2.4
454	bm_sapl(ivm,isapabove,icarbon) = wood_init / (2. + bm_sapl_heartabove + bm_sapl_heartbelow)
455	bm_sapl(ivm,isapbelow,icarbon) = bm_sapl(ivm,isapabove,icarbon)
456	bm_sapl(ivm,iheartabove,icarbon) = bm_sapl_heartabove * bm_sapl(ivm,isapabove,icarbon)
457	bm_sapl(ivm,iheartbelow,icarbon) = bm_sapl_heartbelow * bm_sapl(ivm,isapbelow,icarbon)
458
459	! Use the allometric relationships to calculate initial leaf and root mass
460	bm_sapl(ivm,ileaf,icarbon) = ( bm_sapl(ivm,isapabove,icarbon) + &
461	bm_sapl(ivm,isapbelow,icarbon) ) * KF(ipts,ivm) / ave_tree_height
462	!+++CHECK+++
463	!How do we want to deal with water stress? wstress is accounted for through c0
464	!!$ bm_sapl(ivm,iroot,icarbon) = bm_sapl(ivm,ileaf,icarbon) / ( KF(ipts,ivm) * c0_alloc(ivm) )
465	bm_sapl(ivm,iroot,icarbon) = bm_sapl(ivm,ileaf,icarbon) / &
466	( KF(ipts,ivm) * c0_alloc(ivm) )
467	!+++++++++++
468
469	! Write initial values
470	IF (printlev>=4) THEN
471	WRITE(numout,*) ' Circumference class: ',icir,' PFT type: ',ivm
472	WRITE(numout,*) ' root to sapwood tradeoff p :', c0_alloc(ivm)
473	WRITE(numout,*) 'height_init, dia_init, wood_init, ', &
474	ave_tree_height , dia_init, wood_init
475	WRITE(numout,*) 'pipe_density, ',pipe_density(ivm)
476	ENDIF
477
478	!++++++ CHECK ++++++
479	! The carbohydrate reserves do not seem to be set before this line. This
480	! is a problem since it then uses an unitililized value. Therefore, I
481	! will initilize it.
482	! Should this really be zero? If so the code below is nonesensical and icarbres
483	! could be omitted. nonesensical code = 2 * (bm_sapl(ivm,icir,icarbres,icarbon) + &
484	! bm_sapl(ivm,icir,ileaf,icarbon) + bm_sapl(ivm,icir,iroot,icarbon))
485	bm_sapl(ivm,icarbres,icarbon)=zero
486	!++++++++++++++++++
487
488	! Pools that are defined in the same way for trees and grasses
489	bm_sapl(ivm,ifruit,icarbon) = zero
490
491	!+++CHECK+++
492	! There is an inconsistency in the calculation - most pools are in gN
493	! but leaves is in gC. The correction is proposed, that implies that
494	! the parameter labile_reserve will need to be tuned
495	!!$ bm_sapl(ivm,ilabile,icarbon) = labile_to_total * &
496	!!$ (bm_sapl(ivm,ileaf,icarbon) / cn_leaf_prescribed(ivm) + &
497	!!$ fcn_root(ivm) * bm_sapl(ivm,iroot,icarbon) + fcn_wood(ivm) * &
498	!!$ (bm_sapl(ivm,isapabove,icarbon) + bm_sapl(ivm,isapbelow,icarbon) + &
499	!!$ bm_sapl(ivm,icarbres,icarbon)))
500
501	bm_sapl(ivm,ilabile,icarbon) = labile_to_total * (bm_sapl(ivm,ileaf,icarbon) + &
502	fcn_root(ivm) * bm_sapl(ivm,iroot,icarbon) + fcn_wood(ivm) * &
503	(bm_sapl(ivm,isapabove,icarbon) + bm_sapl(ivm,isapbelow,icarbon) + &
504	bm_sapl(ivm,icarbres,icarbon)))
505	!+++++++++++
506
507	! Avoid deciduous PFTs to have leaves out at establishment
508	! Whether the saplings have leaves or don't have leaves the first year doesn't really matter
509	! Either the approach is correct in the northern hemisphere or in the southern hemisphere. Note
510	! that the resource limitation approach starts with the sapling having leaves.
511	! Anyhow, a spin-up is needed to avoid issues with the initial conditions
512	IF ( pheno_type(ivm) .NE. 1 ) THEN
513
514	! Not evergreen. Deciduous PFTs now need to survive an extra year before bud burst. To
515	! ensure survival there are several options: (a) either the height of the initial
516	! vegetation is increased (this results in more reserves) or (b) the reserves could be
517	! increased. The second option may result in numerical issues further down
518	! the code as the optimal reserve level is calculated from the other biomass pools.
519	! Also some initial tests showed that higher results simply resulted in more respiration.
520	! Use taller trees to start with.
521	bm_sapl(ivm,icarbres,icarbon) = (bm_sapl(ivm,icarbres,icarbon) + &
522	bm_sapl(ivm,ileaf,icarbon) + bm_sapl(ivm,iroot,icarbon))
523	bm_sapl(ivm,ileaf,icarbon) = zero
524	bm_sapl(ivm,iroot,icarbon) = zero
525
526	! When deciduous trees have no leaves they are senescent
527	senescence(ipts,ivm) = .TRUE.
528
529	ELSE
530
531	! Initilize carbohydrate reserves for evergreen PFTs
532	bm_sapl(ivm,icarbres,icarbon) = zero
533
534	! Evergreen trees never go into senescence
535	senescence(ipts,ivm) = .FALSE.
536
537	ENDIF
538
539	IF (printlev>=4) THEN
540	WRITE(numout,*) ' sapling biomass (gC):',icir,ivm,ipts
541	WRITE(numout,*) ' leaves: (::bm_sapl(ivm,ileaf,icarbon))',&
542	bm_sapl(ivm,ileaf,icarbon)
543	WRITE(numout,*) ' sap above ground: (::bm_sapl(ivm,ispabove,icarbon)):',&
544	bm_sapl(ivm,isapabove,icarbon)
545	WRITE(numout,*) ' sap below ground: (::bm_sapl(ivm,isapbelow,icarbon))',&
546	bm_sapl(ivm,isapbelow,icarbon)
547	WRITE(numout,*) ' heartwood above ground: (::bm_sapl(ivm,iheartabove,icarbon))',&
548	bm_sapl(ivm,iheartabove,icarbon)
549	WRITE(numout,*) ' heartwood below ground: (::bm_sapl(ivm,iheartbelow,icarbon))',&
550	bm_sapl(ivm,iheartbelow,icarbon)
551	WRITE(numout,*) ' roots: (::bm_sapl(ivm,iroot,icarbon))',&
552	bm_sapl(ivm,iroot,icarbon)
553	WRITE(numout,*) ' fruits: (::bm_sapl(ivm,ifruit,icarbon))',&
554	bm_sapl(ivm,ifruit,icarbon)
555	WRITE(numout,*) ' carbohydrate reserve: (::bm_sapl(ivm,icarbres,icarbon))',&
556	bm_sapl(ivm,icarbres,icarbon)
557	WRITE(numout,*) ' labile reserve: (::bm_sapl(ivm,ilabile,icarbon))',&
558	bm_sapl(ivm,ilabile,icarbon)
559	ENDIF
560
561
562	cn_leaf=cn_leaf_init(ivm)
563	cn_wood=cn_leaf/fcn_wood(ivm)
564	cn_root=cn_leaf/fcn_root(ivm)
565
566	DO k=1,nparts
567	IF (k.EQ.ileaf) THEN
568	bm_sapl(ivm,k,initrogen) = bm_sapl(ivm,k,icarbon) / cn_leaf
569	ELSE IF (k.LT.iroot) THEN
570	bm_sapl(ivm,k,initrogen) = bm_sapl(ivm,k,icarbon) / cn_wood
571	ELSE
572	bm_sapl(ivm,k,initrogen) = bm_sapl(ivm,k,icarbon) / cn_root
573	ENDIF
574	ENDDO
575
576
577	IF (printlev>=4) THEN
578	WRITE(numout,*)'Initial distribution, method 2',ivm
579	WRITE(numout,*)'Average trees height (m): ',ave_tree_height
580	ENDIF
581
582	!! 2.7 Determine the biomass
583	! I do this based on the biomass in each sapling and the number of trees in each
584	! circumference class...we need the biomass in an average tree
585	biomass(ipts,ivm,:,:)= &
586	bm_sapl(ivm,:,:)*ind(ipts,ivm)
587
588	! The light stress should be calculated making use of Pgap so it accounts for LAI
589	! crown dimensions and tree distribution. However, this would be computationally
590	! expensive so we just use a first order estimate based on light attenuation model
591	! by Lambert-Beer. When LAI is low, a lot of light reaches the forest floor and so
592	! KF should increase to make use of the available light by growing leaves
593	lstress_fac = exp(-biomass_to_lai(biomass(ipts,ivm,ileaf,icarbon),ivm) * 0.5)
594	! Causing large differences between first and second prescribe
595	delta_KF = ABS (KF(ipts,ivm) - ((k_latosa_adapt(ipts,ivm) + &
596	(lstress_fac(ipts,ivm) * (k_latosa_max(ivm)-k_latosa_min(ivm))))) / &
597	( sla(ivm) * tree_ff(ivm) * pipe_density(ivm) ))
598	KF(ipts,ivm) = ((k_latosa_adapt(ipts,ivm) + &
599	(lstress_fac(ipts,ivm) * &
600	(k_latosa_max(ivm)-k_latosa_min(ivm)))) / &
601	( sla(ivm) * tree_ff(ivm) * pipe_density(ivm) ))
602
603	IF(printlev>=4)THEN
604	WRITE(numout,*) 'prescribe delta_KF, ', delta_KF
605	WRITE(numout,*) 'prescribe lstress, ', lstress_fac
606	ENDIF
607	END DO
608
609	IF (printlev>=4) THEN
610	WRITE(numout,*)'Initial biomass distribution'
611	WRITE(numout,*)'End initial biomass distribution',ind(ipts,ivm)
612	WRITE(numout,*)"biomass(ipts,ivm,:,icarbon)",biomass(ipts,ivm,:,icarbon)
613	WRITE(numout,*)"biomass(ipts,ivm,:,initrogen)",biomass(ipts,ivm,:,initrogen)
614	END IF
615
616	IF (ivm .EQ. test_pft .AND. printlev>=4) THEN
617	WRITE(numout,*) 'Check prescribe'
618	WRITE(numout,*) 'stomate_prescribe::init ind ',&
619	ipts,ivm,ind(ipts,ivm)
620	ENDIF
621
622	! PFT is not present
623	ELSE
624
625	! At the stand level
626	biomass(ipts,ivm,:,:) = zero
627	ind(ipts,ivm) = zero
628
629	ENDIF
630
631	! Set leaf age classes, all leaves are current year leaves
632	leaf_frac(ipts,ivm,:) = zero
633	leaf_frac(ipts,ivm,1) = un
634
635	!+++CHECK+++
636	! Set time since last beginning of growing season but only
637	! for the first day of the whole simulation. When the model
638	! is initialized when_growthinit is set to undef. In subsequent
639	! time steps it should have a value. For trees without phenology
640	! the growing season starts at the moment the PFT is prescribed
641	IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
642	when_growthinit(ipts,ivm) = 200
643	ENDIF
644	!+++++++++++
645
646	! Seasonal trees have no leaves at beginning
647	! Saplings of evergreen trees have a leaf mass on day 1 and mass in the other components
648	! saplings of deciduous trees have no leaves on day 1 but mass in the other components.
649	IF ( pheno_model(ivm) .NE. 'none' ) THEN
650
651	! Add the carbon from the leaves to the reserve pool
652	biomass(ipts,ivm,icarbres,:) = biomass(ipts,ivm,icarbres,:) + biomass(ipts,ivm,ileaf,:)
653	biomass(ipts,ivm,ileaf,:) = zero
654	leaf_frac(ipts,ivm,1) = zero
655
656	!+++CHECK+++
657	! Set time since last beginning of growing season but only
658	! for the first day of the whole simulation. When the model
659	! is initialized when_growthinit is set to undef. In subsequent
660	! time steps it should have a value. The phenology module
661	! prevents leaf onset soon after senescence, by setting
662	! ::when_growthinit to a value, leaf offset
663	! will occur at the first opportunity
664	!!$ when_growthinit(ipts,ivm) = large_value
665	IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
666	when_growthinit(ipts,ivm) = 200
667	ENDIF
668	!++++++++++++
669
670	! Redundant, flag has already been set. When there are no leaves, the tree is in senescence
671	senescence(ipts,ivm) = .TRUE.
672
673	ENDIF ! pheno_model(ivm)
674
675	! The biomass to build the saplings is taken from the atmosphere, keep track of
676	! amount to calculate the C-balance closure
677	co2_to_bm(ipts,ivm) = co2_to_bm(ipts,ivm) + ( SUM(biomass(ipts,ivm,:,icarbon)) / dt )
678
679	ENDIF ! tree(ivm)
680
681	!! 2.8 Initialize grassy PFTs
682	!! Use veget_max to check whether the PFT is present. If the PFT is present but it
683	!! has no biomass, prescribe its biomass (gC m-{2}). It is assumed that at day 1
684	!! grasses have all their biomass in the reserve pool. The criteria exclude crops.
685	!! Crops are no longer prescribed but planted the day that begin_leaves is true.
686
687	!+++TEMP+++
688	IF(printlev>=4 .AND. test_pft == ivm)THEN
689	WRITE(numout,*) 'prescribe - total biomass, ',SUM(biomass(ipts,ivm,:,icarbon))
690	ENDIF
691	!++++++++++
692
693	IF ( ( .NOT. is_tree(ivm) ) .AND. &
694	! ( natural(ivm) ) .AND. &
695	( veget_max(ipts,ivm) .GT. min_stomate ) .AND. &
696	( SUM( biomass(ipts,ivm,:,icarbon) ) .LE. min_stomate ) ) THEN
697
698	!+++TEMP+++
699	IF(printlev>=4 .AND. test_pft == ivm)THEN
700	WRITE(numout,*) 'We will prescribe a new vegetation, the old one died'
701	ENDIF
702	!++++++++++
703
704	! For grasses we assume that an individual grass is 1 m2 of grass. This is set
705	! in nmaxtrees in pft_parameters.f90. It could be set to another value but
706	! with the current code this should not have any meaning. The grassland
707	! does not necessarily covers the whole 1 m2 so adjust for the canopy cover
708	ind(ipts,ivm) = nmaxtrees(ivm) * ha_to_m2 * canopy_cover(ivm)
709
710	! now we generate the size of a single grass sapling...this was all taken from
711	! stomate_data.f90...we do not deal with circumference classes for grasses
712	! and crops, but we want to keep the arrays the same as for the trees so
713	! we put the sapling information into the first circumference class
714
715	!+++CHECK+++
716	! Calculate the sapwood to leaf mass in a similar way as has been done for trees.
717	! For trees this approach had been justified by observations. For grasses such
718	! justification is not supported by observations but we didn't try to find it.
719	! Needs more work by someone interested in grasses. There might be a more elegant
720	! solution making use of a well observed parameter.
721	!!$ k_latosa(ipts,ivm) = k_latosa_min(ivm) + &
722	!!$ (wstress_fac(ipts,ivm) * lstress_fac(ipts,ivm) * &
723	!!$ (k_latosa_max(ivm)-k_latosa_min(ivm)))
724	!!$ k_latosa(ipts,ivm) = wstress_fac(ipts,ivm) * (k_latosa_min(ivm) + &
725	!!$ (lstress_fac(ipts,ivm) * &
726	!!$ (k_latosa_max(ivm)-k_latosa_min(ivm))))
727	k_latosa(ipts,ivm) = (k_latosa_adapt(ipts,ivm) + &
728	(lstress_fac(ipts,ivm) * &
729	(k_latosa_max(ivm)-k_latosa_min(ivm))))
730
731	! The mass of the structural carbon relates to the mass of the leaves through
732	! a prescribed parameter ::k_latosa
733	KF(ipts,ivm) = k_latosa(ipts,ivm)
734	!+++++++++++
735
736	! Calculate leaf to root area
737	LF(ipts,ivm) = c0_alloc(ivm) * KF(ipts,ivm)
738
739	!---TEMP---
740	IF(printlev>=4.AND. test_pft == ivm)THEN
741	WRITE(numout,*) 'KF, ', k_latosa(ipts,ivm), KF(ipts,ivm)
742	WRITE(numout,) 'LF, c0_alloc, ', c0_alloc(ivm) KF(ipts,ivm), &
743	c0_alloc(ivm)
744	ENDIF
745	!----------
746
747	! initialize everything to make sure there are not random values floating around
748	! for ncirc != 1
749	bm_sapl(ivm,:,:) = val_exp
750
751	! Similar as for trees, the initial height of the vegetation was defined
752	bm_sapl(ivm,ileaf,icarbon) = height_init_min(ivm) / lai_to_height(ivm) / sla(ivm)
753
754	! Use allometric relationships to define the root mass based on leaf mass. Some
755	! sapwood mass is needed to store the reserves. An arbitrairy fraction of 5% was
756	! used
757	bm_sapl(ivm,iroot,icarbon) = bm_sapl(ivm,ileaf,icarbon) / LF(ipts,ivm)
758	bm_sapl(ivm,isapabove,icarbon) = bm_sapl(ivm,ileaf,icarbon) / KF(ipts,ivm)
759
760	! Some of the biomass components that exist for trees are undefined for grasses
761	bm_sapl(ivm,isapbelow,icarbon) = zero
762	bm_sapl(ivm,iheartabove,icarbon) = zero
763	bm_sapl(ivm,iheartbelow,icarbon) = zero
764	bm_sapl(ivm,ifruit,icarbon) = zero
765	bm_sapl(ivm,icarbres,icarbon) = zero
766
767	! Pools that are defined in the same way for trees and grasses
768	bm_sapl(ivm,ilabile,icarbon) = labile_to_total * (bm_sapl(ivm,ileaf,icarbon) + &
769	fcn_root(ivm) * bm_sapl(ivm,iroot,icarbon) + fcn_wood(ivm) * &
770	(bm_sapl(ivm,isapabove,icarbon) + bm_sapl(ivm,isapbelow,icarbon) + &
771	bm_sapl(ivm,icarbres,icarbon)))
772
773	! Avoid deciduous PFTs to have leaves out at establishment
774	! Whether the saplings have leaves or don't have leaves the first year doesn't really matter
775	! Either the approach is correct in the northern hemisphere or in the southern hemisphere. Note
776	! that the resource limitation approach starts with the sapling having leaves.
777	! Anyhow, a spin-up is needed to avoid issues with the initial conditions
778	IF ( pheno_type(ivm) .NE. 1 ) THEN
779
780	! Not evergreen. Deciduous PFTs now need to survive an extra year before bud burst. To
781	! ensure survival there are several options: (a) either the height of the initial
782	! vegetation is increased (this results in more reserves) or (b) the reserves could be
783	! increased. The second option may result in result in numerical issues further down
784	! the code as the optimal reserve level is calculated from the other biomass pools
785	bm_sapl(ivm,icarbres,icarbon) = bm_sapl(ivm,icarbres,icarbon) + bm_sapl(ivm,ileaf,icarbon) + &
786	bm_sapl(ivm,iroot,icarbon) + bm_sapl(ivm,isapabove,icarbon)
787	bm_sapl(ivm,ileaf,icarbon) = zero
788	bm_sapl(ivm,iroot,icarbon) = zero
789	bm_sapl(ivm,isapabove,icarbon) = zero
790	! When there are no leaves, the crop/grass is in senescence
791	senescence(ipts,ivm) = .TRUE.
792
793	ELSE
794
795	! Initilize carbohydrate reserves for evergreen PFTs
796	bm_sapl(ivm,icarbres,icarbon) = zero
797
798	! Evergreen plants never go into senescence
799	senescence(ipts,ivm) = .FALSE.
800
801	ENDIF
802
803	IF (printlev>=4) THEN
804	WRITE(numout,*) ' sapling biomass (gC):',1,ivm,ipts
805	WRITE(numout,*) ' leaves: (::bm_sapl(ivm,ileaf,icarbon))',&
806	bm_sapl(ivm,ileaf,icarbon)
807	WRITE(numout,*) ' sap above ground: (::bm_sapl(ivm,ispabove,icarbon)):',&
808	bm_sapl(ivm,isapabove,icarbon)
809	WRITE(numout,*) ' sap below ground: (::bm_sapl(ivm,isapbelow,icarbon))',&
810	bm_sapl(ivm,isapbelow,icarbon)
811	WRITE(numout,*) ' heartwood above ground: (::bm_sapl(ivm,iheartabove,icarbon))',&
812	bm_sapl(ivm,iheartabove,icarbon)
813	WRITE(numout,*) ' heartwood below ground: (::bm_sapl(ivm,iheartbelow,icarbon))',&
814	bm_sapl(ivm,iheartbelow,icarbon)
815	WRITE(numout,*) ' roots: (::bm_sapl(ivm,iroot,icarbon))',&
816	bm_sapl(ivm,iroot,icarbon)
817	WRITE(numout,*) ' fruits: (::bm_sapl(ivm,ifruit,icarbon))',&
818	bm_sapl(ivm,ifruit,icarbon)
819	WRITE(numout,*) ' carbohydrate reserve: (::bm_sapl(ivm,icarbres,icarbon))',&
820	bm_sapl(ivm,icarbres,icarbon)
821	WRITE(numout,*) ' labile reserve: (::bm_sapl(ivm,ilabile,icarbon))',&
822	bm_sapl(ivm,ilabile,icarbon)
823	ENDIF
824
825
826
827	cn_leaf=cn_leaf_init(ivm)
828	cn_wood=cn_leaf/fcn_wood(ivm)
829	cn_root=cn_leaf/fcn_root(ivm)
830
831	DO k=1,nparts
832	IF (k.EQ.ileaf) THEN
833	bm_sapl(ivm,k,initrogen) = bm_sapl(ivm,k,icarbon) / cn_leaf
834	ELSE IF (k.LT.iroot) THEN
835	bm_sapl(ivm,k,initrogen) = bm_sapl(ivm,k,icarbon) / cn_wood
836	ELSE
837	bm_sapl(ivm,k,initrogen) = bm_sapl(ivm,k,icarbon) / cn_root
838	ENDIF
839	ENDDO
840
841	! Write initial values
842	IF (printlev>=4) THEN
843	WRITE(numout,*) ' root to sapwood tradeoff (LF) : ', c0_alloc(ivm)
844	WRITE(numout,*) ' grass sapling biomass: ',bm_sapl(ivm,:,icarbon)
845	ENDIF
846
847	! Initial biomass (g C m-2)
848	biomass(ipts,ivm,:,:) = bm_sapl(ivm,:,:) * ind(ipts,ivm)
849	IF (printlev>=4) THEN
850	WRITE(numout,*) 'bm_grass, prescribe, ', biomass(ipts,ivm,:,icarbon)
851	ENDIF
852
853
854	! Set leaf age classes -> all leaves will be current year leaves
855	leaf_frac(ipts,ivm,:) = zero
856	leaf_frac(ipts,ivm,1) = un
857
858	! Set time since last beginning of growing season but only
859	! for the first day of the whole simulation. When the model
860	! is initialized when_growthinit is set to undef. In subsequent
861	! time steps it should have a value.
862	!!$ when_growthinit(ipts,ivm) = large_value
863	IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
864	when_growthinit(ipts,ivm) = 200
865	ENDIF
866
867	! The biomass to build the saplings is taken from the atmosphere, keep track of
868	! amount to calculate the C-balance closure
869	co2_to_bm(ipts,ivm) = co2_to_bm(ipts,ivm) + ( SUM(biomass(ipts,ivm,:,icarbon)) / dt )
870
871	! ELSEIF (.NOT. natural(ivm)) THEN
872	!
873	! ! Initialize croplands - else leaves_begin will never become true
874	! ! Set time since last beginning of growing season but only
875	! ! for the first day of the whole simulation. When the model
876	! ! is initialized when_growthinit is set to undef. In subsequent
877	! ! time steps it should have a value.
878	! IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
879	! when_growthinit(ipts,ivm) = 200
880	! ENDIF
881
882	ENDIF ! .NOT. tree(ivm)
883
884	!! 2.3 Declare PFT present
885	!! Now that the PFT has biomass it should be declared 'present'
886	!! everywhere in that grid box. Assign some additional properties
887	PFTpresent(ipts,ivm) = .TRUE.
888	everywhere(ipts,ivm) = un
889	age(ipts,ivm) = zero
890	npp_longterm(ipts,ivm) = npp_longterm_init
891	lm_lastyearmax(ipts,ivm) = zero
892
893	ENDIF ! not ok_dgvm or agricultural
894
895	ENDIF ! IF (veget_max .GT. zero .AND. ind .EQ. zero)
896
897	ENDDO ! loop over pixels
898
899	ENDDO ! loop over PFTs
900
901	!! 4. Calculate components of the mass balance
902
903	!! 4.1 Calculate final biomass
904	pool_end(:,:,:) = zero
905	DO ipar = 1,nparts
906	DO iele = 1,nelements
907	pool_end(:,:,iele) = pool_end(:,:,iele) + &
908	(biomass(:,:,ipar,iele) * veget_max(:,:))
909	ENDDO
910	ENDDO
911
912	!! 4.2 Calculate mass balance
913	check_intern(:,:,iatm2land,icarbon) = check_intern(:,:,iatm2land,icarbon) + &
914	co2_to_bm(:,:) * veget_max(:,:) * dt
915	check_intern(:,:,iland2atm,icarbon) = -un * zero
916	check_intern(:,:,ilat2out,icarbon) = zero
917	check_intern(:,:,ilat2in,icarbon) = -un * zero
918	check_intern(:,:,ipoolchange,icarbon) = -un * (pool_end(:,:,icarbon) - pool_start(:,:,icarbon))
919	closure_intern = zero
920	DO imbc = 1,nmbcomp
921	closure_intern(:,:,icarbon) = closure_intern(:,:,icarbon) + check_intern(:,:,imbc,icarbon)
922	ENDDO
923
924	!! 4.3 Write outcome
925	DO ipts=1,npts
926	DO ivm=1,nvm
927	IF(ABS(closure_intern(ipts,ivm,icarbon)) .LE. min_stomate)THEN
928	IF (ld_massbal) WRITE(numout,*) 'Mass balance closure in prescribe_prognostic'
929	ELSE
930	WRITE(numout,*) 'Error: mass balance is not closed in prescribe_prognostic'
931	WRITE(numout,*) ' ipts,ivm; ', ipts,ivm
932	WRITE(numout,*) ' Difference is, ', closure_intern(ipts,ivm,icarbon)
933	WRITE(numout,*) ' pool_end,pool_start: ', pool_end(ipts,ivm,icarbon), pool_start(ipts,ivm,icarbon)
934	WRITE(numout,*) ' check_intern,co2_to_bm,pool_end,veget_max: ', &
935	check_intern(ipts,ivm,iatm2land,icarbon),co2_to_bm(ipts,ivm), veget_max(ipts,ivm)
936	IF(ld_stop)THEN
937	CALL ipslerr_p (3,'prescribe_prognostic', 'Mass balance error.','','')
938	ENDIF
939	ENDIF
940	ENDDO
941	ENDDO
942
943	firstcall_prescribe = .FALSE.
944
945	END SUBROUTINE prescribe
946
947
948	END MODULE stomate_prescribe

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

Download in other formats: