Context Navigation

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

Last change on this file since 7346 was 5691, checked in by sebastiaan.luyssaert, 6 years ago
DEV: tested with 13, 37 and 64 PFTs with LCC on different pixels. Some configuration run for 20 years on a given pixel, other crash on another pixel. There is a mass balance problem in sapiens_lcc (ticket #482). This commit fixes a problem with PFT1 in littercalc. This PFT is now fully integrated in LCC and subsequent litter and soil dynamics. veget_max was changed in veget_cov_max where appropriate, a typo in enerbil was corrected.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 79.6 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.
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, cc_to_lai, &
33	cc_to_biomass, check_surface_area, check_mass_balance, &
34	wood_to_qmdia, Nmax, sort_circ_class_biomass
35
36	IMPLICIT NONE
37
38	! private & public routines
39
40	PRIVATE
41	PUBLIC prescribe, make_saplings, prescribe_clear
42
43	! first call
44	LOGICAL, SAVE :: firstcall_prescribe = .TRUE.
45	!$OMP THREADPRIVATE(firstcall_prescribe)
46
47	CONTAINS
48
49	! =================================================================================================================================
50	!! SUBROUTINE : prescribe_clear
51	!!
52	!>\BRIEF : Set the firstcall_prescribe flag back to .TRUE. to prepare for the next simulation.
53	!_=================================================================================================================================
54
55	SUBROUTINE prescribe_clear
56	firstcall_prescribe=.TRUE.
57	END SUBROUTINE prescribe_clear
58
59
60	!! ================================================================================================================================
61	!! SUBROUTINE : prescribe
62	!!
63	!>\BRIEF Works only with static vegetation and agricultural PFTs.
64	!! Initialize biomass, density, presence in the first call
65	!! and update them in the following.
66	!!
67	!! DESCRIPTION (functional, design, flags): \n
68	!! This module works only with static vegetation and agricultural PFTs.
69	!! In the first call, initialize density of individuals, biomass,
70	!! and leaf age distribution to some reasonable value. In the following calls,
71	!! these variables are updated following a stand replacing disturbance or
72	!! if the conditions are suitable for recruitment.
73	!!
74	!! To fulfill these purposes, pipe model are used:
75	!! \latexonly
76	!! \input{prescribe1.tex}
77	!! \input{prescribe2.tex}
78	!! \endlatexonly
79	!!
80	!! RECENT CHANGE(S): None
81	!!
82	!! MAIN OUTPUT VARIABLES(S): ::ind, ::biomass
83	!!
84	!! REFERENCES :
85	!! - Krinner, G., N. Viovy, et al. (2005). "A dynamic global vegetation model
86	!! for studies of the coupled atmosphere-biosphere system." Global
87	!! Biogeochemical Cycles 19: GB1015, doi:1010.1029/2003GB002199.
88	!! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics,
89	!! plant geography and terrestrial carbon cycling in the LPJ dynamic
90	!! global vegetation model, Global Change Biology, 9, 161-185.
91	!! - McDowell, N., Barnard, H., Bond, B.J., Hinckley, T., Hubbard, R.M., Ishii,
92	!! H., KÃ¶stner, B., Magnani, F. Marshall, J.D., Meinzer, F.C., Phillips, N.,
93	!! Ryan, M.G., Whitehead D. 2002. The relationship between tree height and leaf
94	!! area: sapwood area ratio. Oecologia, 132:12â20.
95	!! - Novick, K., Oren, R., Stoy, P., Juang, F.-Y., Siqueira, M., Katul, G. 2009.
96	!! The relationship between reference canopy conductance and simplified hydraulic
97	!! architecture. Advances in water resources 32, 809-819.
98	!! - Ruger et al. 2009, J. of Ecol. doi: 10.1111/j.1365-2745.2009.01552.x
99	!!
100	!! FLOWCHART : None
101	!! \n
102	!_ ================================================================================================================================
103
104	SUBROUTINE prescribe (npts, veget_cov_max, dt, PFTpresent, &
105	everywhere, when_growthinit, leaf_frac, circ_class_dist, &
106	circ_class_n, circ_class_biomass, atm_to_bm, &
107	forest_managed, KF, plant_status, age, npp_longterm, &
108	lm_lastyearmax, tau_eff_leaf, tau_eff_sap, tau_eff_root, &
109	k_latosa_adapt, light_tran_to_level_season, &
110	EndOfYear, lpft_replant,species_change_map)
111
112	!! 0. Parameters and variables declaration
113
114	!! 0.1 Input variables
115
116	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
117	INTEGER(i_std), DIMENSION (:,:), INTENT(in) :: forest_managed !! forest management flag
118	REAL(r_std), INTENT(in) :: dt !! time step (dt_days)
119	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_cov_max !! "maximal" coverage fraction of a PFT
120	!! (LAI -> infinity) on ground. May sum to
121	!! less than unity if the pixel has
122	!! nobio area. (unitless; 0-1)
123	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_root !! Effective root turnover time that accounts
124	!! waterstress (days)
125	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_sap !! Effective sapwood turnover time that accounts
126	!! waterstress (days)
127	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_leaf !! Effective leaf turnover time that accounts
128	!! waterstress (days)
129	LOGICAL, DIMENSION(:,:), OPTIONAL, INTENT(in) :: lpft_replant !! Set to true if a PFT has been clearcut
130	!! and needs to be replaced by another species
131	INTEGER(i_std), DIMENSION (:,:), INTENT(in) :: species_change_map !! map which gives the PFT number that each PFT
132	REAL(r_std), DIMENSION (:,:,:), INTENT(in) :: light_tran_to_level_season !! Seasonal fraction of light transmitted to canopy levels
133	LOGICAL, INTENT(in) :: EndOfYear !! Flag set on the last day of the year used
134	!! to update "yearly variables". This
135	!! variable must be .TRUE. once a year
136	REAL(r_std), DIMENSION(:), INTENT(in) :: circ_class_dist !! The probability distribution of trees
137	!! in a circ class in case of a
138	!! redistribution (unitless).
139
140	!! 0.2 Output variables
141
142
143	!! 0.3 Modified variables
144
145	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! PFT present (0 or 1)
146	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: plant_status !! Growth and phenological status of the plant
147	!! Different stati are listed in in constantes_var
148	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or
149	!! very localized (after its introduction) (?)
150	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! how many days ago was the beginning of
151	!! the growing season (days)
152	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "long term" net primary productivity
153	!! @tex ($gC m^{-2} year^{-1}$) @endtex
154	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! mean age (years)
155	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
156	!! @tex ($gC m^{-2}$) @endtex
157	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: circ_class_n !! Number of individuals in each circ class
158	!! @tex $(ind m^{-2})$ @endtex
159	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: circ_class_biomass !! Biomass components of the model tree
160	!! within a circumference class
161	!! class @tex $(g C ind^{-1})$ @endtex
162	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: atm_to_bm !! CO2 or N taken from the atmosphere to get
163	!! the seed C and N
164	!! to create the seedlings
165	!! @tex (gC.m^{-2}dt^{-1})$ @endtex
166	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age
167	!! class (unitless;0-1)
168	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: KF !! Scaling factor to convert sapwood mass
169	!! into leaf mass (m) - this variable is
170	!! passed to other routines
171	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: k_latosa_adapt !! Leaf to sapwood area adapted for long
172	!! term water stress (m)
173
174
175	!! 0.4 Local variables
176	REAL(r_std), DIMENSION(nvm) :: c0_alloc !! Root to sapwood tradeoff parameter
177	INTEGER(i_std) :: ipts, ivm !! index (unitless)
178	INTEGER(i_std) :: ipar, ilev !! index (unitless)
179	INTEGER(i_std) :: icir, iele, imbc !! index (unitless)
180	INTEGER(i_std) :: deb,fin, imaxt !! index (unitless)
181	REAL(r_std), DIMENSION(npts,nvm) :: k_latosa !! Height dependent base value to calculate
182	!! KF (-)
183	REAL(r_std), DIMENSION(nvm,ncirc,nparts,nelements):: bm_sapl !! Sapling biomass for the functional
184	!! allocation with a dimension for the
185	!! circumference classes
186	!! @tex $(gC.ind^{-1})$ @endtex
187	REAL(r_std), DIMENSION(npts,nvm) :: LF !! Scaling factor to convert sapwood mass
188	!! into root mass (unitless)
189	REAL(r_std), DIMENSION(npts,nvm) :: lstress_fac !! Light stress factor, based on total
190	!! transmitted light (unitless, 0-1)
191	REAL(r_std), DIMENSION(npts,nvm,ncirc) :: circ_class !! circumference of individual trees
192	REAL(r_std), DIMENSION(nparts) :: bm_init !! Biomass needed to initiate the next
193	!! planting @tex $(gC m^{-2})$ @endtex
194	REAL(r_std), DIMENSION(ncirc) :: nb_trees_i !! Number of trees in each twentith
195	!! circumference quantile of the
196	!! distribution (ind)
197	REAL(r_std) :: excedent !! Number of trees after truncation to be
198	!! reallocated to smaller quantiles of the
199	!! distribution (ind)
200	REAL(r_std) :: max_circ_init !! Maximum initial circumferences of the
201	!! truncated exponential distribution (cm)
202	REAL(r_std), DIMENSION(ncirc) :: ave_tree_height !! The height of the ideal tree in each
203	!! circumference class of the
204	!! distribution (ind)...not saved since
205	!! it should only be used for prescribing
206	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements):: check_intern !! Contains the components of the internal
207	!! mass balance chech for this routine
208	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
209	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
210	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
211	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start and end pool of this routine
212	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
213	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! Start and end pool of this routine
214	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
215	REAL(r_std) :: delta_KF !! Difference between new and old estimate
216	!! of KF while iterating
217	REAL(r_std), DIMENSION(npts,nvm) :: veget_cov_max_begin !! temporary storage of veget_cov_max to check area conservation
218	REAL(r_std) :: sapwood_density !! Sapwood density gC m^{-3}
219	REAL(r_std) :: dia_init !! Diameter of a sapling (m)
220	REAL(r_std) :: qmd_init !! Quadratic mean diameter of the initial stand of seedlings (m)
221	REAL(r_std) :: ind_init !! Maximum number of individuals when establishing a new forest (#/m^{-2})
222	REAL(r_std) :: wood_init !! Wood mass of a sapling (g C tree^{-1})
223	INTEGER(i_std) :: iloop !! Index
224	LOGICAL :: establish !! Logical flag to establish a new young vegetation
225	!! after a stand replacing disturbance
226	LOGICAL :: recruite !! Logical flag to add saplings of a PFT under an
227	!! existing canopy of the same PFT
228	REAL(r_std), DIMENSION(npts,nvm,ncirc) :: old_circ_class_n !! Old number of individuals in each circ class
229	!! @tex $(ind m^{-2})$ @endtex
230	REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements) &
231	:: old_circ_class_biomass !! Old biomass components of the model tree
232	!! within a circumference class
233	!! class @tex $(g C ind^{-1})$ @endtex
234	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements) :: old_biomass !! Old stand level biomass
235	!! tex $(gC.m^{-2})$ @endtex
236	REAL(r_std), DIMENSION(ncirc) :: height_small !! Height of the smallest tree in circ_class_biomass
237	!! tex $(m)$ @endtex
238	REAL(r_std) :: dia_small !! Diameter of the smallest tree in circ_class_biomass
239	!! tex $(m)$ @endtex
240	REAL(r_std) :: height_sapl !! Height of the recruitment sapling
241	!! tex $(m)$ @endtex
242	REAL(r_std) :: new_ind !! density of new individuals
243	REAL(r_std) :: old_ind !! density of old individuals
244	REAL(r_std) :: cn_leaf,cn_root !! CN ratio of leaves, and roots
245	!! (gC/gN)
246	REAL(r_std) :: cn_wood !! CN ratio of wood pool
247	!! (gC/gN)
248	REAL(r_std) :: coeff
249	REAL(r_std), DIMENSION(ncirc) :: dummy
250
251	REAL(r_std) :: lambda !! lambda of the truncated exponential
252	!! distribution of initial diameters (1/cm**2)
253	REAL(r_std) :: tmp_ind !! temporary number of individuals per m-2
254	REAL(r_std), DIMENSION(ncirc) :: circ_class_n_new !! variable used to sort circ_class_n
255	REAL(r_std), DIMENSION(ncirc,nparts,nelements) :: circ_class_biomass_new !! variable used to sort circ_class_biomass
256
257	!_ ================================================================================================================================
258
259	IF (printlev>=2) WRITE(numout,*) 'Entering stomate_prescribe.f90'
260
261	!! 1. Initialize biomass at first call
262
263	!! 1.1 Initialize local printlev
264	!printlev_loc=get_printlev('prescribe')
265
266	!! 1.2 Store for closing the mass balance with recruitment
267	old_biomass(:,:,:,:)= cc_to_biomass(npts,nvm,&
268	circ_class_biomass(:,:,:,:,:), &
269	circ_class_n(:,:,:))
270
271
272	!! 1.3 Initialize check for mass balance closure
273	IF (err_act.GT.1) THEN
274
275	check_intern(:,:,:,:) = zero
276	pool_start(:,:,:) = zero
277	DO iele = 1,nelements
278
279	! atm_to_bm has as intent inout, the variable
280	! accumulates carbon over the course of a day.
281	! Use the difference between start and the end of
282	! this routine
283	check_intern(:,:,iatm2land,iele) = - un * &
284	atm_to_bm(:,:,iele) * veget_cov_max(:,:) * dt
285
286	DO ipar = 1,nparts
287
288	DO icir = 1,ncirc
289
290	! Initial biomass pool. Use circ_class_biomass
291	! because that variable is the prognostic variable
292	! biomass is syncronized to circ_class_biomass. If
293	! many diameter classes are used, each diameter class
294	! separatly passes all criteria but when combined
295	! a mass balance problem occurs.
296	pool_start(:,:,iele) = pool_start(:,:,iele) + &
297	circ_class_biomass(:,:,icir,ipar,iele) * &
298	circ_class_n(:,:,icir) * veget_cov_max(:,:)
299
300	ENDDO
301
302	ENDDO
303
304	ENDDO
305
306	!! 1.3 Initialize check for area conservation
307	veget_cov_max_begin(:,:) = veget_cov_max(:,:)
308
309	ENDIF ! err_act.GT.1
310
311	!! 2. Prescribe carbon and nitrogen biomass
312
313	! We want the vegetation to regrow after a stand replacing disturbance
314	! So loop over every point and PFT. If there is supposed to be
315	! vegetation here (according to veget_cov_max) but there isn't (according
316	! to ind), then establish a new young vegetation in that PFT.
317	DO ivm = 2,nvm ! Loop over PFTs
318
319	DO ipts = 1, npts ! Loop over pixels
320
321	!! 2.1 Wait for the end of the year to replant
322	! If we are regrowing different species after managed forests
323	! are killed/die, we sometimes need to prevent them regrowing.
324	! This is only true if they die in the middle of the year due
325	! to natural causes, in which case they will be replanted at
326	! the end of the year. This loop checks to see if we regrow
327	! this PFT now or later. We can regrow a PFT after lpft_replant
328	! is set back to false so before that we need to process all the
329	! information such that the correct species will be regrown.
330	! lpft_replant is an optional variable. First test whether it is
331	! present.
332	IF(ok_change_species)THEN
333	IF(PRESENT(lpft_replant))THEN
334	IF(lpft_replant(ipts,ivm))THEN
335	! .AND. agec_group(ivm)/=species_change_map(ipts,ivm))THEN
336	CYCLE
337	ENDIF
338	ENDIF
339	ENDIF
340
341	! Debug
342	IF(printlev_loc.GT.4)THEN
343	WRITE(numout,*) 'PFT ',ivm
344	WRITE(numout,*) 'plant_status ', plant_status(ipts,ivm)
345	WRITE(numout,*) 'veget_cov_max ', veget_cov_max(ipts,ivm)
346	WRITE(numout,*) 'recruitment_pft ',recruitment_pft(ivm)
347	WRITE(numout,*) 'ok_dgvm ', ok_dgvm
348	WRITE(numout,*) 'ok_recruitment ',ok_recruitment
349	WRITE(numout,*) 'biomass - C, ', &
350	SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1))
351	WRITE(numout,*) 'biomass - N, ', &
352	SUM(SUM(circ_class_biomass(ipts,ivm,:,:,initrogen),1))
353	WRITE(numout,*) 'EndOfYear ', EndOfYear
354	ENDIF
355	!
356
357	! Decide whether we have to establish new vegetation (::establish)
358	! allow recruitment (::recruite). These two local variables were
359	! introduced to allow maximal flexibility so within a single
360	! run, PFTs with and without recruitment can co-exist.
361	IF( ( (plant_status(ipts,ivm) .EQ. iprescribe) .OR. &
362	(plant_status(ipts,ivm) .EQ. idead) ) .AND. &
363	(SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1)) .LE. min_stomate) .AND. &
364	(veget_cov_max(ipts,ivm) .GT. min_stomate) .AND. &
365	(.NOT. ok_dgvm) ) THEN
366
367	! The vegetation just died from replacing disturbances
368	! (two first conditions). The new vegetation could be
369	! of the same or a different PFT than the previous
370	! vegetation.
371	establish = .TRUE.
372	recruite = .FALSE.
373
374	! Debug
375	IF (printlev_loc>=4) THEN
376	WRITE(numout,*) 'We are calling stomate_prescribe.f90 &
377	& to initialize PFT ',ivm
378	ENDIF
379	!-
380
381	ELSEIF( ( (plant_status(ipts,ivm) .NE. iprescribe) .OR. &
382	(plant_status(ipts,ivm) .NE. idead) ) .AND. &
383	(SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1)) .GT. min_stomate) .AND. &
384	(veget_cov_max(ipts,ivm) .GT. min_stomate) .AND. &
385	(recruitment_pft(ivm)) .AND. &
386	(.NOT. ok_dgvm) .AND. &
387	(ok_recruitment) .AND. &
388	(EndOfYear) ) THEN
389
390	! There is already vegetation (first three conditions)
391	! and this specific PFT has recruitment (condition 5)
392	! Only use recruitment when the recruitment flag is TRUE
393	! (condition 7)
394	establish = .FALSE.
395	recruite = .TRUE.
396
397	! Debug
398	IF (printlev_loc>=4) THEN
399	WRITE(numout,*) 'We are calling stomate_prescribe.f90 &
400	& for sapling recruitment in PFT ',ivm
401	ENDIF
402	!-
403
404	ELSEIF( (veget_cov_max(ipts,ivm).LT. min_stomate) .AND. &
405	(SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1)) .LT. min_stomate) .AND. &
406	(SUM(circ_class_n(ipts,ivm,:)) .LT. min_stomate) )THEN
407
408	! The PFT is not defined. Initialize the variables
409	! to ensure that all goes well when moving and merging
410	! different age classes.
411	circ_class_n(ipts,ivm,:) = zero
412	circ_class_biomass(ipts,ivm,:,:,icarbon) = zero
413
414	! Conditions are not fullfilled for establishing
415	! new vegetation
416	establish = .FALSE.
417	recruite = .FALSE.
418
419	ELSEIF (veget_cov_max(ipts,ivm) .GT. min_stomate) THEN
420
421	! Conditions are not fullfilled for either establishing
422	! a new young vegetation or having recruitment under an
423	! exisisting vegetation of the same PFT.
424	establish = .FALSE.
425	recruite = .FALSE.
426
427	! Debug
428	IF (printlev_loc>=4) THEN
429	WRITE(numout,*) 'We are calling stomate_prescribe.f90 &
430	& but nothing should be done for PFT ',ivm
431	ENDIF
432
433
434	ELSE
435
436	! Unexpected condition
437	CALL ipslerr_p (3,'stomate_prescribe.f90',&
438	'Unexpected condition - time to rethink the', &
439	'IF that choses between establish and recruite','')
440
441	ENDIF
442
443	!! 2.2 Calculate variables that will determine the size of the saplings
444	IF( establish .OR. recruite )THEN
445
446	!! 2.2.1 Calculate c0_alloc
447	! The calculation starts with c0_alloc, it only depends on PFT and
448	! the effective longiveties.
449	c0_alloc(ivm) = calculate_c0_alloc(ipts, ivm, &
450	tau_eff_root(ipts,ivm),tau_eff_sap(ipts,ivm))
451
452	!! 2.2.2 Calculate k_latosa for newly established PFTs
453	! If a new vegetation has to be established we recalculate
454	! k_latosa from the parameter file (see note on physiological
455	! memory of k_latosa below). If we add recruits we use the KF of
456	! the existing vegetation.
457	IF (establish) THEN
458
459	! Initialize tree level variables
460	circ_class_n(ipts,ivm,:) = zero
461	circ_class_biomass(ipts,ivm,:,:,:) = zero
462
463	! Assume that there is no physiological memory for k_latosa
464	! between different generations (this is probably true when
465	! trees are planted but seems less likely for natural
466	! regeneration. The first test that implemented memory has
467	! caused problems with the LAI from the second generation
468	! onwards. Nevertheless, we still use k_latosa_adapt. In case
469	! someone want to try again in the future, the structure of
470	! the code is ready to account for memory
471	k_latosa_adapt(ipts,ivm) = k_latosa_min(ivm)
472
473	! Debugging
474	IF (printlev_loc>=4) THEN
475	WRITE(numout,*) 'Establish a PFT (stomate_prescribe.f90):'
476	WRITE(numout,*) ' Imposing initial biomass for prescribed &
477	& trees, initial reserve mass for prescribed grasses.'
478	WRITE(numout,*) ' Declaring prescribed PFTs present.'
479	WRITE(numout,*) ' Grid point: ',ipts,' PFT Type: ',ivm
480	WRITE(numout,*) ' Establish ',establish,' Recruite ',recruite
481	ENDIF
482	!-
483
484	ENDIF ! establish
485
486	!! 2.2.3 Calculate Lightstress
487	! Note that light and water stress have an opposite effect on KF.
488	! Waterstress will decrease KF because more C should be allocated
489	! to the roots. Light stress will increase KF because more C should
490	! be allocated to the leaves.
491	! Lightstress varies from 0 to 1 and is calculated from the canopy
492	! structure (veget).
493	IF (establish) THEN
494
495	! Given that there is no vegetation at this point, the
496	! lightstress cannot be calculated and is set to 1.
497	! However, as soon as we grow a canopy it will experience
498	! light stress and so KF should be adjusted.
499	lstress_fac(ipts,ivm) = un
500
501	ELSEIF (recruite) THEN
502
503	! We already have vegetation, so we should do recruitment if
504	! we want it for this PFT. First lets compute light stress based
505	! on the fraction of light reaching the ground (level 1),
506	! averaged over all daytime time steps (sinang>0) of the vegetative
507	! period (gpp>0). This is based on Pgap and accounts for canopy
508	! structure.
509	lstress_fac(ipts,ivm)=light_tran_to_level_season(ipts,ivm,1)
510
511	! Debug
512	IF(printlev_loc>=4)THEN
513	WRITE(numout,*) 'Mean seasonal transmitted light used in &
514	& stomate_prescribe.f90 ' // &
515	& '(ivm,ipts,nlevels_tot) ',ivm,ipts,nlevels_tot
516	DO ilev = nlevels_tot,1,-1
517	WRITE(numout,*) ilev,light_tran_to_level_season(ipts,ivm,ilev)
518	ENDDO
519	WRITE(numout,*) 'lstress_fac from CN-CAN is ', &
520	lstress_fac(ipts,ivm)*100, ' (%)'
521	ENDIF
522	!-
523
524	ENDIF ! establish/recruite for light
525
526	ENDIF ! establish .OR. recruite for initialization
527
528	!! 2.3 Initial vegetation has to be prescribed only when the vegetation is
529	! static
530	IF (establish .OR. recruite) THEN
531
532	!! 2.3.1 Initialize woody PFT's
533	! Use veget_cov_max to check whether the PFT is present. If the PFT
534	! is present but it has no biomass, prescribe its biomass (gC m-{2}).
535	IF ( is_tree(ivm) ) THEN
536
537	! prescribe/recruitment branching for actual process
538	IF (establish) THEN
539
540	IF (ncirc .EQ. 1) THEN
541
542	! Calculate the height range based on the prescribed
543	! min/max initial heights. In the general approach
544	! we divide by (ncirc-1) which equals to zero if ncirc
545	! is one. Therefore a special case was defined.
546	ave_tree_height(ncirc) = (height_init_min(ivm) + &
547	height_init_max(ivm))/2
548
549	ELSE
550
551	DO icir = 1,ncirc
552
553	! Calculate the height range based on the prescribed
554	! min/max initial heights. Use the general approach.
555	ave_tree_height(icir) = height_init_min(ivm) + &
556	(icir-1) * (height_init_max(ivm) - &
557	height_init_min(ivm))/(ncirc-1)
558
559	ENDDO
560
561	ENDIF !ncirc.EQ.1
562
563	! Allocation factors
564	! Sapwood to root ratio
565	! Following Magnani et al. 2000 "In order to decreases hydraulic
566	! resistance, the investment of carbon in fine roots or sapwood
567	! yields to the plant very different returns, both because of
568	! different hydraulic conductivities and because of the strong
569	! impact of plant height on shoot resistance. On the other hand,
570	! fine roots and sapwood have markedly different longevities and
571	! the cost of production, discounted for turnover, will differ
572	! accordingly. Optimal growth under hydraulic constraints requires
573	! that the ratio of marginal hydraulic returns to marginal annual
574	! cost for carbon investment in either roots or sapwood be the
575	! same (Bloom, Chapin & Mooney 1985; Case & Fair 1989). This
576	! is formalized in equation (13) and further derived to obtain
577	! equation (17) in Magnani et al 2000. The latter is implemented
578	! here.Pipe_density is given in gC/m-3, convert to kg/m-3. And
579	! apply equation (17) in Magnani et al 2000. Note that c0_alloc
580	! was calculated at the start of this routine. The calculation
581	! itself is done in function_library
582
583	! Calculate leaf area to sapwood area
584	! To be consistent with the hydraulic limitations and pipe theory,
585	! k_latosa is calculated from equation (18) in Magnani et al.
586	! To do so, total hydraulic resistance and tree height need to
587	! known. This poses a problem as the resistance depends on the leaf
588	! area and the leaf area on the resistance. There is no independent
589	! equation and equations 12 and 18 depend on each other and
590	! substitution would be circular. Hence prescribed k_latosa values
591	! were obtained from observational records and are given in
592	! mtc_parameters.f90.
593
594	! The relationship between height and k_latosa as reported in
595	! McDowell et al 2002 and Novick et al 2009 could be implemented
596	! to adjust k_latosa for the height of the stand. This did NOT
597	! result in a realistic model behavior
598	! k_latosa(ipts,ivm) = wstress_fac(ipts,ivm) * &
599	! (k_latosa_max(ivm) - (latosa_height(ivm) * &
600	! (SUM( nb_trees_i(:) * ave_tree_height(:) ) / &
601	! SUM( nb_trees_i(:) ))))
602
603	! Also k_latosa has been reported to be a function of CO2
604	! concentration (Atwell et al. 2003, Tree Physiology, 23:13-21
605	! and Pakati et al. 2000, Global Change Biology, 6:889-897).
606	! This effect is not accounted for in the current code
607
608	! Finally, k_latosa is also reported to be a function of
609	! diameter (i.e. stand thinning, Simonin et al 2006, Tree
610	! Physiology, 26:493-503). Here the relationship with thinning was
611	! interpreted as a realtionship with light stress. Note that light
612	! stress cannot be calculated at this time in the model because
613	! there is no canopy (that's why we are in prescribe!) and
614	! so there is no lightstress. lstress was therefore set to one
615	! (see above). We prefered to keep this redundancy in the code
616	! because it makes it clear that k_latosa is calculated in the
617	! same way in prescribe.f90 and growth_fun_all.f90
618	k_latosa(ipts,ivm) = (k_latosa_adapt(ipts,ivm) + &
619	(lstress_fac(ipts,ivm) * &
620	(k_latosa_max(ivm)-k_latosa_min(ivm))))
621
622	! +++CHECK+++
623	! Ideally waterstress and lightstress should both
624	! affect k_latosa. One of the following approaches could
625	! be used as a starting point.
626	! k_latosa(ipts,ivm) = k_latosa_min(ivm) + &
627	! (wstress_fac(ipts,ivm) * lstress_fac(ipts,ivm) * &
628	! (k_latosa_max(ivm)-k_latosa_min(ivm)))
629	! k_latosa(ipts,ivm) = wstress_fac(ipts,ivm) * (k_latosa_min(ivm) + &
630	! (lstress_fac(ipts,ivm) * &
631	! (k_latosa_max(ivm)-k_latosa_min(ivm))))
632	! ++++++++++++
633
634	! Calculate conversion coefficient for sapwood area to leaf
635	! area
636	! (1) The scaling parameter between leaf and sapwood mass
637	! is derived from LA_ind = k_latosa * SA_ind, where
638	! LA_ind = leaf area of an individual, SA_ind is the sapwood
639	! area of an individual and k_latosa a pipe-model parameter
640	! (2) LA_ind = Cl * sla
641	! (3) Cs = SA_ind * height * wooddensity * tree_ff
642	! Substitute (2) and (3) in (1)
643	! Cl = Cs * k1 / (wooddensity * sla * tree_ff * height)
644	! Cl = Cs*KF/height, where KF is in (m)
645	! KF is passed to the allocation routine and it is saved in the
646	! restart file.
647	KF(ipts,ivm) = k_latosa(ipts,ivm) / &
648	( sla(ivm) * pipe_density(ivm) * tree_ff(ivm))
649
650	! Initialize delta_KF to get the DO WHILE started
651	delta_KF = un
652	iloop=0
653	DO WHILE (delta_KF .GT. max_delta_KF)
654
655	! If there is a WHILE loop, there always needs to be a check
656	! on the number of loops to make sure we don't get stuck in
657	! an infinite loop. This number is completely arbitrary.
658	iloop=iloop+1
659	IF(iloop > 1000)THEN
660
661	WRITE(numout,*) 'Taking too long to converge the '//&
662	'delta_KF loop in prescribe!'
663	WRITE(numout,*) 'iloop,delta_KF,max_delta_KF,ivm,ipts: ',&
664	iloop,delta_KF,max_delta_KF,ivm,ipts
665	CALL ipslerr_p (3,'stomate_prescribe',&
666	'Taking too long to converge the delta_KF','','')
667
668	ENDIF
669
670
671	! Now we create the saplings for each class, based on the
672	! height
673	DO icir = 1,ncirc
674
675	! The assumption we make is that we plant trees of 2 to 3
676	! years old rather than growing trees from seeds. The
677	! allometric relationship between height and diameter is
678	! derived from mature tree and likely unrealistic for
679	! saplings. The height of the saplings is prescribed and
680	! determines the reserves which are especially important
681	! for deciduous species which need to survive on their
682	! reserves for the first year (new phenology scheme requires
683	! annual mean values to get started). These are aboveground
684	! values for dia_init and wood_init.
685	dia_init = ( ave_tree_height(icir) / pipe_tune2(ivm) ) ** &
686	( 1. / pipe_tune3(ivm) )
687	wood_init = ( ave_tree_height(icir) * pi / 4. * &
688	(dia_init) ** 2. ) * pipe_density(ivm) * tree_ff(ivm)
689
690	! Use the variables from the pipe model and the allometric
691	! relationships to calculate the different biomass components
692	! of the saplings in each diameter class.
693	CALL make_saplings(wood_init, bm_sapl, &
694	ave_tree_height(icir), icir, ivm, &
695	ipts, npts, KF, c0_alloc, plant_status, &
696	establish)
697
698	ENDDO ! ncirc
699
700	! Determine the biomass in each circumference class.
701	! Calculations are based on the biomass in each sapling and
702	! the number of trees in each circumference class. Total
703	! leaf biomass is needed for the light stress calculation
704	! below.
705	DO icir=1,ncirc
706	circ_class_biomass(ipts,ivm,icir,:,:) = &
707	bm_sapl(ivm,icir,:,:)
708	ENDDO
709
710	! The light stress should be calculated making use of Pgap so
711	! it accounts for LAI crown dimensions and tree distribution.
712	! However, this would be computationally expensive so we just
713	! use a first order estimate based on light attenuation model
714	! by Lambert-Beer. When LAI is low, a lot of light reaches the
715	! forest floor and so KF should increase to make use of the
716	! available light by growing leaves
717	lstress_fac = &
718	exp(-cc_to_lai(circ_class_biomass(ipts,ivm,:,ileaf,icarbon),&
719	circ_class_n(ipts,ivm,:),ivm))
720
721	! Causing large differences between first and second prescribe
722	delta_KF = ABS (KF(ipts,ivm) - ((k_latosa_adapt(ipts,ivm) + &
723	(lstress_fac(ipts,ivm) * &
724	(k_latosa_max(ivm)-k_latosa_min(ivm))))) / &
725	( sla(ivm) * tree_ff(ivm) * pipe_density(ivm) ))
726	KF(ipts,ivm) = ((k_latosa_adapt(ipts,ivm) + &
727	(lstress_fac(ipts,ivm) * &
728	(k_latosa_max(ivm)-k_latosa_min(ivm)))) / &
729	( sla(ivm) * tree_ff(ivm) * pipe_density(ivm) ))
730
731	! Debug
732	IF(printlev_loc>=4)THEN
733	WRITE(numout,*) 'prescribe delta_KF, ', delta_KF
734	WRITE(numout,*) 'prescribe lstress, ', lstress_fac
735	ENDIF
736	! -
737
738	END DO ! DO WHILE delta_KF
739
740	! The biomass of the saplings is known. Calculate the
741	! quadratic mean diameter and the max number of trees
742	! the stand can carry according to the self-thinning
743	! relationship. Calculate the tree distribution by making
744	! use of the prescribed number of trees
745	circ_class_n(ipts,ivm,:) = circ_class_dist(:) / &
746	SUM(circ_class_dist(:)) * nmaxtrees(ivm) * ha_to_m2
747
748	! Calculate the quadratic mean diameter of the prescribed
749	! stand.
750	qmd_init = &
751	wood_to_qmdia(circ_class_biomass(ipts,ivm,:,:,icarbon),&
752	circ_class_n(ipts,ivm,:),ivm)
753
754	! Calculate the expected density under self-thinning.
755	! unit is trees.ha-1
756	ind_init = Nmax(qmd_init*m_to_cm,ivm)
757
758	! Take the lowest density. The reason for taking the lowest
759	! is that the highest density could be up to 50 10E6 trees
760	! per hectare. Given that our allometric rules are for
761	! mature trees, this then results in an LAI of ~980 which
762	! makes the albedo routine crash (because of quality checks).
763	! Because every individual tree hardly grows, self-thinning
764	! goes way to slow. For these reason we start with a more
765	! reasonable number of trees per hectare and will use a
766	! background mortality to make the number of trees decrease
767	IF (nmaxtrees(ivm) .GT. ind_init) THEN
768
769	! Recalculate stand density and diameter distribution
770	circ_class_n(ipts,ivm,:) = circ_class_dist(:) / &
771	SUM(circ_class_dist(:)) * ind_init * ha_to_m2
772
773	ENDIF
774
775	! Debug
776	IF (printlev_loc>=4) THEN
777	DO icir=1,ncirc
778	WRITE(numout,*)'icir, circ_class_n, ', icir, &
779	circ_class_n(ipts,ivm,icir)
780	WRITE(numout,*)'Initial biomass distribution - carbon'
781	WRITE(numout,*) circ_class_biomass(ipts,ivm,icir,:,icarbon)
782	WRITE(numout,*) circ_class_biomass(ipts,ivm,icir,:,initrogen)
783	ENDDO
784	WRITE(numout,*)' Total number of saplings', &
785	SUM(circ_class_n(ipts,ivm,:))
786	END IF
787	!-
788
789	ELSEIF (recruite) THEN
790
791	!! Calculate the number of recruits
792	! Recruitment is based on the light stress factor [0 1] computed
793	! above. It makes use of the functional response of recruitment
794	! to light described by Ruger et al (2009, J. of Ecol.,
795	! doi: 10.1111/j.1365-2745.2009.01552.x) in the 50-ha plot of Barro
796	! Colorado Island (BCI) in Panama:
797	! Log10(Nrecruits)= A + B * ( Log10(lstress_fac) - MU )
798	!
799	! Where Nrecruits is the number of recruits for an area of 25 m2.
800	! The parameter A is the intercept (mean log10 recruits expected),
801	! B the slope (functional response), and MU=Log10(0.02)=-1.7
802	! Here the recruitment response to light can be negative (B<0),
803	! decelerating (0<B<1), accelerating (B>1) or linear (B=1).
804	! MU is needed here to express the fraction of light between 0 and 1
805	! centred on mean light 2% at BCI. The linearized (log-log) power
806	! model is implemented here for each PFT following the original setup
807	! of Ruger et al 2009 to ensure realistic values. IT COULD EASILY BE
808	! SIMPLIFIED. The predicted number of recruits is for an area of
809	! 25 m2 so we have to compute the antilogarithm and divide by 25 to
810	! get the estimate by 1 m2 which is the unit used in ORCHIDEE.
811	!
812	! Add min_stomate (a very small number) to lstress_fac to avoid numerical
813	! issues with LOG10(zero) being undefined.
814	IF (recruitment_alpha(ivm).LT.min_stomate .AND. &
815	recruitment_beta(ivm).LT.min_stomate) THEN
816
817	! Sanity check. If the users sets both values to zero
818	! new_ind is set to zero. The simulation should be
819	! identical to a simulation without recruitment.
820	new_ind = zero
821
822	ELSE
823
824	! Calculate the number of recruits
825	new_ind = ( 10**( recruitment_alpha(ivm) + &
826	recruitment_beta(ivm) * &
827	( LOG10(lstress_fac(ipts,ivm)+min_stomate) - &
828	LOG10(0.02) ) ) )/25
829
830	ENDIF
831
832	! Update the number of individuals
833	tmp_ind = SUM(circ_class_n(ipts,ivm,:))
834	old_ind = tmp_ind
835	tmp_ind = tmp_ind + new_ind
836
837	! Debug
838	IF(printlev_loc>=4)THEN
839	WRITE(numout,*) 'PFT, ',ivm
840	WRITE(numout,*) 'canopy cover (LAI), ',&
841	cc_to_lai(circ_class_biomass(ipts,ivm,:,ileaf,icarbon),&
842	circ_class_n(ipts,ivm,:),ivm)
843	WRITE(numout,) 'used canopy_gap (%) is, ', lstress_fac(ipts,ivm)100
844	WRITE(numout,) 'Original ind is, ', (tmp_ind - new_ind)m2_to_ha
845	WRITE(numout,) 'new_ind is, ', new_indm2_to_ha
846	WRITE(numout,) 'ind+new_ind is, ', tmp_indm2_to_ha
847	qmd_init = &
848	wood_to_qmdia(circ_class_biomass(ipts,ivm,:,:,icarbon),&
849	circ_class_n(ipts,ivm,:),ivm)
850	WRITE(numout,*) 'Dg (m) ', qmd_init
851	ind_init = Nmax(qmd_init*m_to_cm,ivm)
852	WRITE(numout,*) 'Nmax, ',ind_init
853	ENDIF
854	!-
855
856	!! Calculate the dimensions of the recruits
857	! The height of the saplings is prescribed and determines the reserves
858	! which are especially important for deciduous species, which need to
859	! survive on their reserves for the first couple of months to a whole
860	! year because the phenology scheme requires annual mean values to get
861	! started. Calculate sapling diameter and biomass for prescribed sapling
862	! height
863	dia_init = ( recruitment_height(ivm) / pipe_tune2(ivm) ) ** &
864	( 1. / pipe_tune3(ivm) )
865	wood_init = (recruitment_height(ivm) * pi / 4. * (dia_init) ** 2. ) * &
866	pipe_density(ivm) * tree_ff(ivm)
867
868	! Use the variables from the pipe model and the allometric
869	! relationships to calculate the different biomass components
870	! of the saplings in each diameter class.
871	CALL make_saplings(wood_init, bm_sapl, recruitment_height(ivm), &
872	1, ivm, ipts, npts, KF, c0_alloc, plant_status, establish)
873
874	! Debug
875	IF (printlev_loc>=4) THEN
876	WRITE(numout,*)'Prescribed recruitment height (m): ',recruitment_height(ivm)
877	WRITE(numout,*)'Recruitment diameter (m) for prescribed height: ', dia_init
878	WRITE(numout,*)'Wood_init recruitment: ',wood_init
879	WRITE(numout,*)'ind_init ',ind_init
880	WRITE(numout,*)'bm_sapl after calling make_saplings: ',bm_sapl(ivm,1,:,:)
881	ENDIF
882	!-
883
884	! Update the number of individuals in the first size class. All
885	! recruitment occurs in the circ_class with the smallest diameter
886	old_circ_class_n(ipts,ivm,:) = circ_class_n(ipts,ivm,:)
887	circ_class_n(ipts,ivm,1) = circ_class_n(ipts,ivm,1) + new_ind
888
889	!! Determine the biomass in the first circumference class.
890	! This is a dillution approach were an average tree is made by making use
891	! of the biomass in the available trees and the biomass in the recruits.
892	! This is a simple mean weighted by the number of trees. The model needs
893	! the biomass in an individual tree (gC tree-1)
894	DO iele = 1,nelements
895	old_circ_class_biomass(ipts,ivm,:,:,iele) = &
896	circ_class_biomass(ipts,ivm,:,:,iele)
897	circ_class_biomass(ipts,ivm,1,:,iele) = &
898	(circ_class_biomass(ipts,ivm,1,:,iele) * &
899	old_circ_class_n(ipts,ivm,1) + &
900	bm_sapl(ivm,1,:,iele) * new_ind ) / &
901	circ_class_n(ipts,ivm,1)
902	ENDDO
903
904	! Debug
905	IF (printlev_loc>=4) THEN
906	WRITE(numout,) "Canopy gap (%)", lstress_fac(ipts,ivm)100
907	WRITE(numout,*) "Original density (N/ha): ",&
908	old_circ_class_n(ipts,ivm,:)*m2_to_ha
909	WRITE(numout,) "New density (N/ha): ",circ_class_n(ipts,ivm,:)m2_to_ha
910	WRITE(numout,*) "Recruited saplings per ha in first size class ", &
911	new_ind * m2_to_ha
912	WRITE(numout,*) "Original biomass per individual", &
913	SUM(old_circ_class_biomass(ipts,ivm,:,:,icarbon),2),&
914	SUM(old_circ_class_biomass(ipts,ivm,:,:,initrogen),2)
915	WRITE(numout,*) "New biomass per individual", &
916	SUM( circ_class_biomass(ipts,ivm,:,:,icarbon),2 ),&
917	SUM( circ_class_biomass(ipts,ivm,:,:,initrogen),2 )
918	WRITE(numout,*) "Biomass of new saplings gC/sapling", &
919	bm_sapl(ivm,1,:,icarbon)
920	WRITE(numout,*) "Biomass of new saplings gN/sapling", &
921	bm_sapl(ivm,1,:,initrogen)
922	WRITE(numout,*) "KF ", KF(ipts,ivm)
923	ENDIF
924	!-
925
926	ENDIF ! end prescribe/recruitment block
927
928	IF (establish) THEN
929
930	! Set leaf age classes, all leaves are current year leaves
931	leaf_frac(ipts,ivm,:) = zero
932	leaf_frac(ipts,ivm,1) = un
933
934	! Set time since last beginning of growing season but only
935	! for the first day of the whole simulation. When the model
936	! is initialized when_growthinit is set to undef. In subsequent
937	! time steps it should have a value. For trees without phenology
938	! the growing season starts at the moment the PFT is prescribed
939	IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
940	when_growthinit(ipts,ivm) = 200
941	ENDIF
942
943	ENDIF
944
945	! The carbon and nitrogen to build the saplings is taken from the
946	! atmosphere, keep track of amount to calculate the C-balance
947	! closure
948	IF (establish) THEN
949
950	! The whole biomass of this PFT was taken from the
951	! atmosphere
952	atm_to_bm(ipts,ivm,:) = atm_to_bm(ipts,ivm,:) + &
953	(SUM(cc_to_biomass(npts, ivm, &
954	circ_class_biomass(ipts,ivm,:,:,:), &
955	circ_class_n(ipts,ivm,:)),1) / dt )
956
957	ELSEIF (recruite) THEN
958
959	! Only the increase in biomass of this PFT was taken
960	! from the atmosphere
961	atm_to_bm(ipts,ivm,:) = atm_to_bm(ipts,ivm,:) + &
962	(SUM(cc_to_biomass(npts, ivm, &
963	circ_class_biomass(ipts,ivm,:,:,:), &
964	circ_class_n(ipts,ivm,:)) - &
965	old_biomass(ipts,ivm,:,:),1) / dt )
966
967	ENDIF
968
969	! Check whether the order was preserved. We only want to
970	! preserve the order for the carbon biomass. Nitrogen follows
971	! carbon. We expect that the biomass of an individual tree
972	! increases with an increasing circ class. This is probably
973	! not important for the correct functioning of the model. So,
974	! if this turns out to be computationally expensive it is
975	! worth trying without. We try to preserve the order as
976	! an additional mean to control the quality of the simulations.
977	! A strict order also helps to interpret the output. Hence,
978	! if the order was not preserved we will sort the biomasses.
979	! Note that this subroutine first checks whether the order
980	! was preseverd.
981	circ_class_biomass_new = circ_class_biomass(ipts,ivm,:,:,:)
982	circ_class_n_new = circ_class_n(ipts,ivm,:)
983	CALL sort_circ_class_biomass(circ_class_biomass_new, &
984	circ_class_n_new)
985	circ_class_biomass(ipts,ivm,:,:,:) = circ_class_biomass_new
986	circ_class_n(ipts,ivm,:) = circ_class_n_new
987
988	ENDIF ! tree(ivm)
989
990	!! 2.3.2 Initialize grassy PFTs
991	! Use veget_cov_max to check whether the PFT is present. If the PFT is
992	! present but it has no biomass, prescribe its biomass (gC m-{2}).
993	! It is assumed that at day 1 grasses have all their biomass in the
994	! reserve pool. The criteria exclude crops. Crops are no longer
995	! prescribed but planted the day that begin_leaves is true.
996	IF (establish) THEN
997
998	IF ( ( .NOT. is_tree(ivm) ) .AND. &
999	( natural(ivm) ) .AND. &
1000	( veget_cov_max(ipts,ivm) .GT. min_stomate ) .AND. &
1001	( SUM(circ_class_biomass(ipts,ivm,1,:,icarbon) ) .LE. min_stomate ) ) THEN
1002
1003	! Debug
1004	IF(printlev_loc>=4)THEN
1005	WRITE(numout,*) 'We will prescribe a new grass '//&
1006	'vegetation, the old one died'
1007	ENDIF
1008	! -
1009
1010	! For grasses we assume that an individual grass is 1 m2 of grass.
1011	! This is set in nmaxtrees in pft_parameters.f90. It could be set
1012	! to another value but with the current code this should not have
1013	! any meaning. The grassland does not necessarily covers the whole
1014	! 1 m2 so adjust for the canopy cover
1015	tmp_ind = nmaxtrees(ivm) * ha_to_m2 * canopy_cover(ivm)
1016
1017	!+++CHECK+++
1018	! We do not deal with circumference classes for grasses
1019	! and crops, but we want to keep the arrays the same as for the trees
1020	! so we put the young grass information into the first circumference
1021	! class. Calculate the sapwood to leaf mass in a similar way as has
1022	! been done for trees. For trees this approach had been justified by
1023	! observations. For grasses such justification is not supported by
1024	! observations but we didn't try to find it. Needs more work by
1025	! someone interested in grasses. There might be a more elegant
1026	! solution making use of a well observed parameter.
1027	k_latosa(ipts,ivm) = (k_latosa_adapt(ipts,ivm) + &
1028	(lstress_fac(ipts,ivm) * &
1029	(k_latosa_max(ivm)-k_latosa_min(ivm))))
1030
1031	! The mass of the structural carbon relates to the mass of the leaves
1032	! through a prescribed parameter ::k_latosa
1033	KF(ipts,ivm) = k_latosa(ipts,ivm)
1034	!+++++++++++
1035
1036	! Calculate leaf to root area
1037	LF(ipts,ivm) = c0_alloc(ivm) * KF(ipts,ivm)
1038
1039	! Debug
1040	IF(printlev_loc>=4.AND. test_pft == ivm)THEN
1041	WRITE(numout,*) 'KF, ', k_latosa(ipts,ivm), KF(ipts,ivm)
1042	WRITE(numout,) 'LF, c0_alloc, ', c0_alloc(ivm) KF(ipts,ivm), &
1043	c0_alloc(ivm)
1044	ENDIF
1045	!-
1046
1047	! initialize everything to make sure there are not random values
1048	! floating around for ncirc = 1
1049	bm_sapl(ivm,1,:,:) = val_exp
1050
1051	! Similar as for trees, the initial height of the vegetation was
1052	! defined
1053	bm_sapl(ivm,1,ileaf,icarbon) = height_init_min(ivm) / &
1054	lai_to_height(ivm) / sla(ivm)
1055
1056	! Use allometric relationships to define the root mass based on
1057	! leaf mass. Some 'sapwood mass' is needed to store the reserves.
1058	! An arbitrairy fraction of 5% was used
1059	bm_sapl(ivm,1,iroot,icarbon) = bm_sapl(ivm,1,ileaf,icarbon) / &
1060	LF(ipts,ivm)
1061	bm_sapl(ivm,1,isapabove,icarbon) = bm_sapl(ivm,1,ileaf,icarbon) / &
1062	KF(ipts,ivm)
1063
1064	! Some of the biomass components that exist for trees are undefined
1065	! for grasses
1066	bm_sapl(ivm,1,isapbelow,icarbon) = zero
1067	bm_sapl(ivm,1,iheartabove,icarbon) = zero
1068	bm_sapl(ivm,1,iheartbelow,icarbon) = zero
1069	bm_sapl(ivm,1,ifruit,icarbon) = zero
1070	bm_sapl(ivm,1,icarbres,icarbon) = zero
1071
1072	! Pools that are defined in the same way for trees and grasses
1073	bm_sapl(ivm,1,ilabile,icarbon) = labile_to_total * &
1074	(bm_sapl(ivm,1,ileaf,icarbon) + &
1075	fcn_root(ivm) * bm_sapl(ivm,1,iroot,icarbon) + fcn_wood(ivm) * &
1076	(bm_sapl(ivm,1,isapabove,icarbon) + &
1077	bm_sapl(ivm,1,isapbelow,icarbon) + &
1078	bm_sapl(ivm,1,icarbres,icarbon)))
1079
1080	! Allocate the nitrogen
1081	cn_leaf=cn_leaf_init(ivm)
1082	cn_wood=cn_leaf/fcn_wood(ivm)
1083	cn_root=cn_leaf/fcn_root(ivm)
1084
1085	!+++CHECK+++
1086	! See comments at a similmar block of code for the
1087	! tree saplings
1088	! Use the C/N ratio to calculate the N content for
1089	! all the biomass components.
1090	bm_sapl(ivm,1,:,initrogen) = &
1091	bm_sapl(ivm,1,:,icarbon) / cn_root
1092
1093	! Overwrite the N content for the leaves and wood
1094	! components. Note that only sapabove is defined the
1095	! other wood components are zero
1096	bm_sapl(ivm,1,ileaf,initrogen) = &
1097	bm_sapl(ivm,1,ileaf,icarbon) / cn_leaf
1098	bm_sapl(ivm,1,isapabove,initrogen) = &
1099	bm_sapl(ivm,1,isapabove,icarbon) / cn_wood
1100	!++++++++++++
1101
1102	! Avoid deciduous PFTs to have leaves out at establishment
1103	! Whether the saplings have leaves or don't have leaves the first
1104	! year doesn't really matter. Either the approach is correct in the
1105	! northern hemisphere or in the southern hemisphere. Note
1106	! that the resource limitation approach starts with the sapling having
1107	! leaves. Anyhow, a spin-up is recommended to avoid issues with the
1108	! initial conditions
1109	IF ( pheno_type(ivm) .NE. 1 ) THEN
1110
1111	! Not evergreen. Deciduous PFTs now need to survive half a year
1112	! before bud burst. To ensure survival there are several options:
1113	! (a) either the height of the initial vegetation is increased
1114	! (this results in more reserves) or (b) the reserves could be
1115	! increased. The second option may result in result in numerical
1116	! issues further down the code as the optimal reserve level is
1117	! calculated from the other biomass pools
1118	bm_sapl(ivm,1,icarbres,:) = bm_sapl(ivm,1,icarbres,:) + &
1119	bm_sapl(ivm,1,ileaf,:) + &
1120	bm_sapl(ivm,1,iroot,:) + &
1121	bm_sapl(ivm,1,isapabove,:)
1122	bm_sapl(ivm,1,ileaf,:) = zero
1123	bm_sapl(ivm,1,iroot,:) = zero
1124	bm_sapl(ivm,1,isapabove,:) = zero
1125
1126	! When there are no leaves, the crop/grass is dormant
1127	plant_status(ipts,ivm) = idormant
1128
1129	ELSE
1130
1131	! Initilize carbohydrate reserves for evergreen PFTs
1132	bm_sapl(ivm,1,icarbres,:) = zero
1133
1134	! Evergreen plants always have a canopy
1135	plant_status(ipts,ivm) = icanopy
1136
1137	ENDIF
1138
1139	! Debug
1140	IF (printlev_loc>=4) THEN
1141	DO iele=1,nelements
1142	WRITE(numout,*) ' young grass biomass (gC):',ivm,ipts
1143	WRITE(numout,*) ' bm_sapl(:,1,ileaf,:,) ',&
1144	bm_sapl(ivm,1,ileaf,iele)
1145	WRITE(numout,*) ' bm_sapl(:,1,ispabove,:) ',&
1146	bm_sapl(ivm,1,isapabove,iele)
1147	WRITE(numout,*) ' bm_sapl(:,1,isapbelow,:) ',&
1148	bm_sapl(ivm,1,isapbelow,iele)
1149	WRITE(numout,*) ' bm_sapl(:,1,iheartabove,:) ',&
1150	bm_sapl(ivm,1,iheartabove,iele)
1151	WRITE(numout,*) ' bm_sapl(:,1,iheartbelow,:) ',&
1152	bm_sapl(ivm,1,iheartbelow,iele)
1153	WRITE(numout,*) ' bm_sapl(:,1,iroot,:) ',&
1154	bm_sapl(ivm,1,iroot,iele)
1155	WRITE(numout,*) ' bm_sapl(:,1,ifruit,:)', &
1156	bm_sapl(ivm,1,ifruit,iele)
1157	WRITE(numout,*) ' bm_sapl(:,1,icarbres,:)', &
1158	bm_sapl(ivm,1,icarbres,iele)
1159	WRITE(numout,*) ' bm_sapl(:,1,ilabile,:)', &
1160	bm_sapl(ivm,1,ilabile,iele)
1161	ENDDO
1162	ENDIF
1163	! -
1164
1165	! Synchronize biomass and circ_class_biomass (gC tree-1). This
1166	! allows to more easily check the mass balance closure
1167	circ_class_biomass(ipts,ivm,1,:,:) = bm_sapl(ivm,1,:,:)
1168	circ_class_n(ipts,ivm,1) = tmp_ind
1169
1170	! Set leaf age classes -> all leaves will be current year leaves
1171	leaf_frac(ipts,ivm,:) = zero
1172	leaf_frac(ipts,ivm,1) = un
1173
1174	! Set time since last beginning of growing season but only
1175	! for the first day of the whole simulation. When the model
1176	! is initialized when_growthinit is set to undef. In subsequent
1177	! time steps it should have a value.
1178	IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
1179	when_growthinit(ipts,ivm) = 200
1180	ENDIF
1181
1182	! The carbon and nitrogen biomass to build the saplings is taken
1183	! from the atmosphere, keep track of amount to calculate the mass
1184	! balance closure
1185	DO iele = 1,nelements
1186	atm_to_bm(ipts,ivm,iele) = atm_to_bm(ipts,ivm,iele) + &
1187	((SUM(circ_class_biomass(ipts,ivm,1,:,iele))*&
1188	circ_class_n(ipts,ivm,1)) / dt)
1189	ENDDO
1190
1191	! Debug
1192	IF (printlev_loc>=4) THEN
1193	DO iele=1,nelements
1194	WRITE(numout,*) ' root to sapwood tradeoff (LF) : ', c0_alloc(ivm)
1195	WRITE(numout,*) ' grass sapling biomass: ',bm_sapl(ivm,1,:,iele)
1196	ENDDO
1197	ENDIF
1198	! -
1199
1200	ELSEIF (.NOT. natural(ivm)) THEN
1201
1202	! Initialize croplands - else leaves_begin will never become true
1203	! Set time since last beginning of growing season but only
1204	! for the first day of the whole simulation. When the model
1205	! is initialized when_growthinit is set to undef. In subsequent
1206	! time steps it should have a value.
1207	IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
1208	when_growthinit(ipts,ivm) = 200
1209	ENDIF
1210
1211	! Debug
1212	IF (printlev_loc>=4) WRITE(numout,*) 'Initialize crops',ivm
1213	!-
1214
1215	ENDIF ! .NOT. tree(ivm)
1216
1217	ENDIF ! establish (for .NOT. is_tree(ivm))
1218
1219	!! 2.3.3 Declare PFT present
1220	!! Now that the PFT has biomass it should be declared 'present'
1221	!! everywhere in that grid box. Assign some additional properties
1222	IF (establish) THEN
1223	PFTpresent(ipts,ivm) = .TRUE.
1224	everywhere(ipts,ivm) = un
1225	age(ipts,ivm) = zero
1226	npp_longterm(ipts,ivm) = npp_longterm_init
1227	lm_lastyearmax(ipts,ivm) = zero
1228	ENDIF
1229
1230	ENDIF ! establish or recruite
1231
1232	ENDDO ! loop over pixels
1233
1234	ENDDO ! loop over PFTs
1235
1236	!! 3. Check numerical consistency of this routine
1237
1238	IF (err_act.GT.1) THEN
1239
1240	! 3.1 Check surface area
1241	CALL check_surface_area("stomate_prescribe", npts, veget_cov_max_begin, &
1242	veget_cov_max,'pft')
1243
1244	! 3.2 Mass balance closure
1245	! 3.2.1 Calculate final biomass
1246	pool_end(:,:,:) = zero
1247	DO ipar = 1,nparts
1248	DO iele = 1,nelements
1249	DO icir = 1,ncirc
1250	pool_end(:,:,iele) = pool_end(:,:,iele) + &
1251	(circ_class_biomass(:,:,icir,ipar,iele) * &
1252	circ_class_n(:,:,icir) * veget_cov_max(:,:))
1253	ENDDO
1254	ENDDO
1255	ENDDO
1256
1257	! 3.2.2 Calculate mass balance
1258	! Common processes
1259	DO iele=1,nelements
1260	check_intern(:,:,iatm2land,iele) = check_intern(:,:,iatm2land,iele) + &
1261	atm_to_bm(:,:,iele) * veget_cov_max(:,:) * dt
1262	check_intern(:,:,ipoolchange,iele) = -un * (pool_end(:,:,iele) - &
1263	pool_start(:,:,iele))
1264	ENDDO
1265
1266	closure_intern(:,:,:) = zero
1267	DO imbc = 1,nmbcomp
1268	DO iele=1,nelements
1269	! Debug
1270	IF (printlev_loc>=3) THEN
1271	IF(iele .EQ. 1) WRITE(numout,"(9999(G12.5,:,','))") &
1272	'imbc, test_pft, check_intern', test_grid , imbc, &
1273	test_pft, check_intern(test_grid,test_pft,imbc,iele)
1274	ENDIF
1275	!-
1276	closure_intern(:,:,iele) = closure_intern(:,:,iele) + &
1277	check_intern(:,:,imbc,iele)
1278	ENDDO
1279	ENDDO
1280
1281	! 3.3 Check mass balance closure
1282	CALL check_mass_balance("stomate_prescribe", closure_intern, npts, &
1283	pool_end, pool_start, veget_cov_max, 'pft')
1284
1285	ENDIF ! err_act.GT.1
1286
1287	IF(printlev>=3) WRITE(numout,*) 'Leaving stomate_prescribe.f90'
1288
1289	firstcall_prescribe = .FALSE.
1290
1291	END SUBROUTINE prescribe
1292
1293
1294
1295	!! ================================================================================================================================
1296	!! SUBROUTINE : make_sapling
1297	!!
1298	!>\BRIEF Given the dimensions and the allocation factors, calculate the biomass of a sapling
1299	!!
1300	!! DESCRIPTION (functional, design, flags):
1301	!!
1302	!! RECENT CHANGE(S): None
1303	!!
1304	!! MAIN OUTPUT VARIABLES(S): ::bm_sapl, ::KF
1305	!!
1306	!! REFERENCES :
1307	!!
1308	!! FLOWCHART : None
1309	!! \n
1310	!_ ================================================================================================================================
1311
1312	SUBROUTINE make_saplings(wood_init, bm_sapl, tree_height, icir, &
1313	ivm, ipts, npts, KF, c0_alloc, plant_status, &
1314	establish)
1315
1316	!! 0. Parameters and variables declaration
1317
1318	!! 0.1 Input variables
1319	REAL(r_std), INTENT(in) :: wood_init
1320	REAL(r_std), INTENT(in) :: tree_height !! Tree height (m)
1321	INTEGER(i_std), INTENT(in) :: ivm, icir, ipts !! index (unitless)
1322	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
1323	REAL(r_std), DIMENSION(:), INTENT(in) :: c0_alloc !! Root to sapwood tradeoff parameter
1324	LOGICAL, INTENT(in) :: establish !! yes: establish a new young vegetation.
1325	!! No: add saplings under an existing vegetation
1326
1327	!! 0.2 Output variables
1328
1329	!! 0.3 Modified variables
1330	REAL(r_std), DIMENSION(:,:,:,:),INTENT(inout) :: bm_sapl !! Sapling biomass for the functional
1331	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: KF !! Scaling factor to convert sapwood mass
1332	!! into leaf mass (m) - this variable is
1333	!! passed to other routines
1334	REAL, DIMENSION(:,:), INTENT(inout) :: plant_status !! Growth and phenological status of the plant
1335	!! Different stati defined in constants
1336	!! 0.4 Local variables
1337	INTEGER(i_std) :: iele !! index (unitless)
1338	REAL(r_std) :: cn_leaf,cn_root !! CN ratio of leaves and roots (gC/gN)
1339	REAL(r_std) :: cn_wood !! CN ratio of wood pool (gC/gN)
1340
1341
1342
1343	!_ ================================================================================================================================
1344
1345	! The woody biomass is contained in four components. Thus,
1346	! wood_init = isapabove + isapbelow + iheartabove + iheartbelow.
1347	! Given that isapbelow = isapbelow and iheartabove =
1348	! iheartbelow = bm_sapl_heartabove*isapabove. If
1349	! bm_sapl_heartbelow = bm_sapl_heartabove = 0.2, then
1350	! isapabove = wood_init/2.4
1351	bm_sapl(ivm,icir,isapabove,icarbon) = wood_init / &
1352	(2. + bm_sapl_heartabove + bm_sapl_heartbelow)
1353	bm_sapl(ivm,icir,isapbelow,icarbon) = bm_sapl(ivm,icir,isapabove,icarbon)
1354	bm_sapl(ivm,icir,iheartabove,icarbon) = bm_sapl_heartabove * &
1355	bm_sapl(ivm,icir,isapabove,icarbon)
1356	bm_sapl(ivm,icir,iheartbelow,icarbon) = bm_sapl_heartbelow * &
1357	bm_sapl(ivm,icir,isapbelow,icarbon)
1358
1359	! Use the allometric relationships to calculate initial leaf
1360	! and root mass. Note that wstress should be accounted for
1361	! through c0_alloc. There was an inconsistency in the
1362	! original calculation (OCN) - most pools were in gN but
1363	! leaves were in gC. A correction is proposed. Given that
1364	! icarbes was just set to zero the equation for ilabile is
1365	! a bit overkill and icarbes could be omitted.
1366	bm_sapl(ivm,icir,ileaf,icarbon) = ( bm_sapl(ivm,icir,isapabove,icarbon) + &
1367	bm_sapl(ivm,icir,isapbelow,icarbon) ) * KF(ipts,ivm) / tree_height
1368	bm_sapl(ivm,icir,iroot,icarbon) = bm_sapl(ivm,icir,ileaf,icarbon) / &
1369	( KF(ipts,ivm) * c0_alloc(ivm) )
1370	bm_sapl(ivm,icir,icarbres,icarbon)=zero
1371	bm_sapl(ivm,icir,ifruit,icarbon) = zero
1372	bm_sapl(ivm,icir,ilabile,icarbon) = labile_to_total * &
1373	(bm_sapl(ivm,icir,ileaf,icarbon) + &
1374	fcn_root(ivm) * bm_sapl(ivm,icir,iroot,icarbon) + &
1375	fcn_wood(ivm) * (bm_sapl(ivm,icir,isapabove,icarbon) + &
1376	bm_sapl(ivm,icir,isapbelow,icarbon) + &
1377	bm_sapl(ivm,icir,icarbres,icarbon)))
1378
1379	! Debug
1380	IF (printlev_loc>=4) THEN
1381	WRITE(numout,*) ' +++++ CHECK INSIDE make_sapling ++++ '
1382	WRITE(numout,*) 'Circumference class: ',icir,' PFT type: ',ivm
1383	WRITE(numout,*) ' root to sapwood tradeoff p :', c0_alloc(ivm)
1384	WRITE(numout,*) ' sapling height, sapling diameter, wood_init, ', &
1385	tree_height , ( tree_height/pipe_tune2(ivm) ) ** &
1386	( 1. / pipe_tune3(ivm) ), wood_init
1387	WRITE(numout,*) 'pipe_density, ',pipe_density(ivm)
1388	ENDIF
1389	! -
1390
1391	! Allocate the nitrogen
1392	cn_leaf=cn_leaf_init(ivm)
1393	cn_wood=cn_leaf/fcn_wood(ivm)
1394	cn_root=cn_leaf/fcn_root(ivm)
1395
1396	! Use the C/N ratio to calculate the N content for
1397	! all the biomass components.
1398	bm_sapl(ivm,icir,:,initrogen) = &
1399	bm_sapl(ivm,icir,:,icarbon) / cn_root
1400
1401	! Overwrite the N content for the leaves and wood
1402	! components.
1403	bm_sapl(ivm,icir,ileaf,initrogen) = &
1404	bm_sapl(ivm,icir,ileaf,icarbon) / cn_leaf
1405	bm_sapl(ivm,icir,isapabove,initrogen) = &
1406	bm_sapl(ivm,icir,isapabove,icarbon) / cn_wood
1407	bm_sapl(ivm,icir,isapbelow,initrogen) = &
1408	bm_sapl(ivm,icir,isapbelow,icarbon) / cn_wood
1409	bm_sapl(ivm,icir,iheartabove,initrogen) = &
1410	bm_sapl(ivm,icir,iheartabove,icarbon) / cn_wood
1411	bm_sapl(ivm,icir,iheartbelow,initrogen) = &
1412	bm_sapl(ivm,icir,iheartbelow,icarbon) / cn_wood
1413
1414	! Avoid deciduous PFTs to have leaves out at establishment
1415	! Whether the saplings have leaves or don't have leaves the
1416	! first year doesn't really matter. Either the approach is
1417	! correct in the northern hemisphere or in the southern
1418	! hemisphere. Anyhow, a spin-up is needed to avoid issues
1419	! with the initial conditions
1420	IF ( pheno_type(ivm) .NE. 1 ) THEN
1421
1422	! Not evergreen. Deciduous PFTs now need to survive an extra
1423	! year before bud burst. To ensure survival there are several
1424	! options: (a) either the height of the initial vegetation is
1425	! increased (this results in more reserves) or (b) the reserves
1426	! could be increased. The second option may result in numerical
1427	! issues further down the code as the optimal reserve level is
1428	! calculated from the other biomass pools. Also some initial
1429	! tests showed that higher reserves simply resulted in more
1430	! respiration.
1431	! Tip: use taller trees to start with if they die between
1432	! prescribe and leaf onset
1433	bm_sapl(ivm,icir,icarbres,:) = (bm_sapl(ivm,icir,icarbres,:) + &
1434	bm_sapl(ivm,icir,ileaf,:) + bm_sapl(ivm,icir,iroot,:))
1435	bm_sapl(ivm,icir,ileaf,:) = zero
1436	bm_sapl(ivm,icir,iroot,:) = zero
1437
1438	! When deciduous trees have no leaves they are dormant
1439	! If the code is used for recruitment use the actual
1440	! plant_status of the overstorey trees for the saplings
1441	IF (establish) THEN
1442	plant_status(ipts,ivm) = idormant
1443	ENDIF
1444
1445	ELSE
1446
1447	! Initilize carbohydrate reserves for evergreen PFTs
1448	bm_sapl(ivm,icir,icarbres,:) = zero
1449
1450	! Evergreen trees always have a canopy. If the code
1451	! is used for recruitment use the actual plant_status
1452	! of the overstorey trees for the saplings
1453	IF (establish) THEN
1454	plant_status(ipts,ivm) = icanopy
1455	ENDIF
1456
1457	ENDIF
1458
1459	! Debug
1460	IF (printlev_loc>=4) THEN
1461	DO iele=1,nelements
1462	WRITE(numout,*) ' sapling biomass (gC):',icir,ivm,ipts
1463	WRITE(numout,*) ' bm_sapl(:,:,ileaf,:,) ',&
1464	bm_sapl(ivm,icir,ileaf,iele)
1465	WRITE(numout,*) ' bm_sapl(:,:,ispabove,:) ',&
1466	bm_sapl(ivm,icir,isapabove,iele)
1467	WRITE(numout,*) ' bm_sapl(:,:,isapbelow,:) ',&
1468	bm_sapl(ivm,icir,isapbelow,iele)
1469	WRITE(numout,*) ' bm_sapl(:,:,iheartabove,:) ',&
1470	bm_sapl(ivm,icir,iheartabove,iele)
1471	WRITE(numout,*) ' bm_sapl(:,:,iheartbelow,:) ',&
1472	bm_sapl(ivm,icir,iheartbelow,iele)
1473	WRITE(numout,*) ' bm_sapl(:,:,iroot,:) ',&
1474	bm_sapl(ivm,icir,iroot,iele)
1475	WRITE(numout,*) ' bm_sapl(:,:,ifruit,:)', &
1476	bm_sapl(ivm,icir,ifruit,iele)
1477	WRITE(numout,*) ' bm_sapl(:,:,icarbres,:)', &
1478	bm_sapl(ivm,icir,icarbres,iele)
1479	WRITE(numout,*) ' bm_sapl(:,:,ilabile,:)', &
1480	bm_sapl(ivm,icir,ilabile,iele)
1481	ENDDO
1482	ENDIF
1483
1484	END SUBROUTINE make_saplings
1485
1486	END MODULE stomate_prescribe

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

Download in other formats: