Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

lpj_establish.f90

Last change on this file was 1902, checked in by matthew.mcgrath, 11 years ago
DEV: Trunk merges up to and including r1892
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 36.0 KB

Line
1	! =================================================================================================================================
2	! MODULE : lpj_establish
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 Establish pft's
10	!!
11	!!\n DESCRIPTION: None
12	!!
13	!! RECENT CHANGE(S): None
14	!!
15	!! REFERENCE(S) :
16	!! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics,
17	!! plant geography and terrestrial carbon cycling in the LPJ dynamic
18	!! global vegetation model, Global Change Biology, 9, 161-185.\n
19	!! - Haxeltine, A. and I. C. Prentice (1996), BIOME3: An equilibrium
20	!! terrestrial biosphere model based on ecophysiological constraints,
21	!! resource availability, and competition among plant functional types,
22	!! Global Biogeochemical Cycles, 10(4), 693-709.\n
23	!! - Smith, B., I. C. Prentice, et al. (2001), Representation of vegetation
24	!! dynamics in the modelling of terrestrial ecosystems: comparing two
25	!! contrasting approaches within European climate space,
26	!! Global Ecology and Biogeography, 10, 621-637.\n
27	!!
28	!! SVN :
29	!! $HeadURL$
30	!! $Date$
31	!! $Revision$
32	!! \n
33	!_ ================================================================================================================================
34
35	MODULE lpj_establish
36
37	! modules used:
38	USE xios_orchidee
39	USE ioipsl_para
40	USE stomate_data
41	USE constantes
42	USE function_library
43
44	IMPLICIT NONE
45
46	! private & public routines
47	PRIVATE
48	PUBLIC establish_clear, establish_diagnostic
49
50	LOGICAL, SAVE :: firstcall = .TRUE. !! first call
51	!$OMP THREADPRIVATE(firstcall)
52	CONTAINS
53
54
55	!! ================================================================================================================================
56	!! SUBROUTINE : fire_clear
57	!!
58	!>\BRIEF Set the firstcall flag to .TRUE. and activate initialization
59	!_ ================================================================================================================================
60
61	SUBROUTINE establish_clear
62	firstcall = .TRUE.
63	END SUBROUTINE establish_clear
64
65
66	! =================================================================================================================================
67	! SUBROUTINE : establish_diagnostic
68	!
69	!>\BRIEF Calculate sstablishment of new woody PFT and herbaceous PFTs
70	!!
71	!! DESCRIPTION : Establishments of new woody and herbaceous PFT are simulated.
72	!! Maximum establishment rate (0.12) declines due to competition for light (space).
73	!! There are two establishment estimates: one for the for DGVM and one for the
74	!! static cases.\n
75	!! In the case of DGVM, competitive process of establishment for the area of
76	!! available space is represented using more detailed description compared with static
77	!! one. Biomass and distribution of plant age are updated on the basis of changes
78	!! in number of individuals. Finally excess sapwood of is converted to heartwood.
79	!!
80	!! \latexonly
81	!! \input{equation_lpj_establish.tex}
82	!! \endlatexonly
83	!! \n
84	!!
85	!! RECENT CHANGE(S): None
86	!!
87	!! REFERENCE(S) :
88	!! Smith, B., I. C. Prentice, et al. (2001), Representation of vegetation
89	!! dynamics in the modelling of terrestrial ecosystems: comparing two
90	!! contrasting approaches within European climate space,
91	!! Global Ecology and Biogeography, 10, 621-637.
92	!!
93	!! FLOWCHART :
94	!! \latexonly
95	!! \includegraphics[scale = 0.7]{establish.png}
96	!! \endlatexonly
97	!! \n
98	!_ ================================================================================================================================
99
100	SUBROUTINE establish_diagnostic (npts, dt, PFTpresent, regenerate, &
101	neighbours, resolution, need_adjacent, herbivores, &
102	precip_annual, gdd0, lm_lastyearmax, &
103	cn_ind, lai, avail_tree, avail_grass, npp_longterm, &
104	leaf_age, leaf_frac, &
105	ind, biomass, age, everywhere, co2_to_bm,veget_max, woodmass_ind)
106
107	!! 0. Variable and parameter declaration
108
109	!! 0.1 Input variables
110
111	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of pixels (dimensionless)
112	REAL(r_std), INTENT(in) :: dt !! Time step of vegetation dynamics for stomate
113	!! (days)
114	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: PFTpresent !! Is pft there (unitless)
115	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: regenerate !! Winter sufficiently cold (unitless)
116	INTEGER(i_std), DIMENSION(npts,8), INTENT(in) :: neighbours !! indices of the 8 neighbours of each grid point
117	!! (unitless);
118	!! [1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW]
119	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: resolution !! resolution at each grid point (m); 1=E-W, 2=N-S
120	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: need_adjacent !! in order for this PFT to be introduced, does it
121	!! have to be present in an adjacent grid box?
122	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: herbivores !! time constant of probability of a leaf to
123	!! be eaten by a herbivore (days)
124	REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_annual !! annual precipitation (mm year^{-1})
125	REAL(r_std), DIMENSION(npts), INTENT(in) :: gdd0 !! growing degree days (degree C)
126	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
127	!! (gC m^{-2 })
128	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: cn_ind !! crown area of individuals (m^2)
129	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lai !! leaf area index OF an individual plant
130	!! (m^2 m^{-2})
131	REAL(r_std), DIMENSION(npts), INTENT(in) :: avail_tree !! space availability for trees (unitless)
132	REAL(r_std), DIMENSION(npts), INTENT(in) :: avail_grass !! space availability for grasses (unitless)
133	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: npp_longterm !! longterm NPP, for each PFT (gC m^{-2})
134	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT
135	!! (LAI -> infinity) on ground (unitless)
136	!! 0.2 Output variables
137
138	!! 0.3 Modified variables
139
140	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_age !! leaf age (days)
141	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class (unitless)
142	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind !! Number of individuals (individuals m^{-2})
143	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout):: biomass !! biomass (gC m^{-2 })
144	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age !! mean age (years)
145	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: everywhere !! is the PFT everywhere in the grid box or very
146	!! localized (unitless)
147	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: co2_to_bm !! biomass up take for establishment i.e.
148	!! pseudo-photosynthesis (gC m^{-2} day^{-1})
149	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: woodmass_ind !! woodmass of the individual, needed to calculate
150	!! crownarea in lpj_crownarea (gC m^{-2 })
151
152	!! 0.4 Local variables
153
154	REAL(r_std) :: tau_eatup !! time during which a sapling can be entirely
155	!! eaten by herbivores (days)
156	REAL(r_std), DIMENSION(npts,nvm) :: fpc_nat !! new fpc, foliage projective cover: fractional
157	!! coverage (unitless)
158	REAL(r_std), DIMENSION(npts) :: estab_rate_max_climate_tree !! maximum tree establishment rate,
159	!! based on climate only (unitless)
160	REAL(r_std), DIMENSION(npts) :: estab_rate_max_climate_grass !! maximum grass establishment rate,
161	!! based on climate only (unitless)
162	REAL(r_std), DIMENSION(npts) :: estab_rate_max_tree !! maximum tree establishment rate,
163	!! based on climate and fpc
164	!! (unitless)
165	REAL(r_std), DIMENSION(npts) :: estab_rate_max_grass !! maximum grass establishment rate,
166	!! based on climate and fpc
167	!! (unitless)
168	REAL(r_std), DIMENSION(npts) :: sumfpc !! total natural fpc (unitless)
169	REAL(r_std), DIMENSION(npts) :: fracnat !! total fraction occupied by natural
170	!! vegetation (unitless)
171	REAL(r_std), DIMENSION(npts) :: sumfpc_wood !! total woody fpc (unitless)
172	REAL(r_std), DIMENSION(npts) :: spacefight_tree !! for trees, measures the total concurrence for
173	!! available space (unitless)
174	REAL(r_std), DIMENSION(npts) :: spacefight_grass!! for grasses, measures the total concurrence
175	!! for available space (unitless)
176	REAL(r_std), DIMENSION(npts,nvm) :: d_ind !! change in number of individuals per time step
177	!! (individuals m^{-2} day{-1})
178	REAL(r_std), DIMENSION(npts) :: bm_new !! biomass increase (gC m^{-2 })
179	REAL(r_std), DIMENSION(npts) :: dia !! stem diameter (m)
180	REAL(r_std), DIMENSION(npts) :: b1 !! temporary variable
181	REAL(r_std), DIMENSION(npts) :: woodmass !! woodmass of an individual (gC m^{-2})
182	REAL(r_std), DIMENSION(npts) :: leaf_mass_young !! carbon mass in youngest leaf age class
183	!! (gC m^{-2})
184	REAL(r_std), DIMENSION(npts) :: factor !! reduction factor for establishment if many
185	!! trees or grasses are present (unitless)
186	REAL(r_std), DIMENSION(npts) :: total_bm_c !! Total carbon mass for all pools (gC m^{-2})
187	REAL(r_std), DIMENSION(npts,nelements) :: total_bm_sapl !! Total sappling biomass for all pools
188	!! (gC m^{-2})
189	INTEGER(i_std) :: nfrontx !! from how many sides is the grid box invaded
190	!! (unitless?)
191	INTEGER(i_std) :: nfronty !! from how many sides is the grid box invaded
192	!! (unitless?)
193	!LOGICAL, DIMENSION(npts) :: many_new !! daily establishment rate is large compared to
194	!! present number of individuals
195	REAL(r_std), DIMENSION(npts) :: vn !! flow due to new individuals veget_max after
196	!! establishment, to get a proper estimate of
197	!! carbon and nitrogen
198	REAL(r_std), DIMENSION(npts) :: lai_ind !! lai on each PFT surface (m^2 m^{-2})
199	INTEGER(i_std) :: i,j,k,m !! indices (unitless)
200	!_ ================================================================================================================================
201
202	IF (bavard.GE.3) WRITE(numout,*) 'Entering establish'
203
204	!! 1. messages and initialization
205
206	! Assumption: time durasion that sapling is completely eaten by hervioures is a half year?
207	! No reference
208	tau_eatup = one_year/2.
209
210	!! 1.1 First call only
211	IF ( firstcall ) THEN
212
213	WRITE(numout,*) 'establish:'
214
215	WRITE(numout,*) ' > time during which a sapling can be entirely eaten by herbivores (d): ', &
216	tau_eatup
217
218	firstcall = .FALSE.
219
220	ENDIF
221
222	!! 2. recalculate fpc
223
224	IF (control%ok_dgvm) THEN
225	fracnat(:) = un
226
227	!! 2.1 Only natural part of the grid cell
228	do j = 2,nvm ! Loop over # PFTs
229
230	IF ( .NOT. natural(j) ) THEN
231	fracnat(:) = fracnat(:) - veget_max(:,j)
232	ENDIF
233	ENDDO ! Loop over # PFTs
234
235	sumfpc(:) = zero
236
237	!! 2.2 Total natural fpc on grid
238	! The overall fractional coverage of a PFT in a grid is calculated here.
239	! FPC is related to mean individual leaf area index by the Lambert-Beer law.
240	! See Eq. (1) in tex file.\n
241	DO j = 2,nvm ! Loop over # PFTs
242	IF ( natural(j) ) THEN
243	WHERE(fracnat(:).GT.min_stomate)
244	WHERE (lai(:,j) == val_exp)
245	fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) / fracnat(:)
246	ELSEWHERE
247	fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) / fracnat(:) &
248	* ( un - exp( - lm_lastyearmax(:,j) * sla(j) * ext_coeff(j) ) )
249	ENDWHERE
250	ENDWHERE
251
252	WHERE ( PFTpresent(:,j) )
253	sumfpc(:) = sumfpc(:) + fpc_nat(:,j)
254	ENDWHERE
255	ELSE
256
257	fpc_nat(:,j) = zero
258
259	ENDIF
260
261	ENDDO ! Loop over # PFTs
262
263	!! 2.3 Total woody fpc on grid and number of regenerative tree pfts
264	! Total woody FPC increases by adding new FPC.
265	! Under the condition that temperature in last winter is higher than a threshold,
266	! woody plants is exposed in higher competitive environment.
267	sumfpc_wood(:) = zero
268	spacefight_tree(:) = zero
269
270	DO j = 2,nvm ! Loop over # PFTs
271
272	IF ( is_tree(j) .AND. natural(j) ) THEN
273
274	! total woody fpc
275	WHERE ( PFTpresent(:,j) )
276	sumfpc_wood(:) = sumfpc_wood(:) + fpc_nat(:,j)
277	ENDWHERE
278
279	! how many trees are competing? Count a PFT fully only if it is present
280	! on the whole grid box.
281	WHERE ( PFTpresent(:,j) .AND. ( regenerate(:,j) .GT. regenerate_crit ) )
282	spacefight_tree(:) = spacefight_tree(:) + everywhere(:,j)
283	ENDWHERE
284
285	ENDIF
286
287	ENDDO ! Loop over # PFTs
288
289	!! 2.4 Total number of natural grasses on grid\n
290	! Grass increment equals 'everywhere'\n
291	spacefight_grass(:) = zero
292
293	DO j = 2,nvm ! Loop over # PFTs
294
295	IF ( .NOT. is_tree(j) .AND. natural(j) ) THEN
296
297	! Count a PFT fully only if it is present on a grid.
298	WHERE ( PFTpresent(:,j) )
299	spacefight_grass(:) = spacefight_grass(:) + everywhere(:,j)
300	ENDWHERE
301
302	ENDIF
303
304	ENDDO ! Loop over # PFTs
305
306	!! 2.5 Maximum establishment rate, based on climate only\n
307	WHERE ( ( precip_annual(:) .GE. precip_crit ) .AND. ( gdd0(:) .GE. gdd_crit_estab ) )
308
309	estab_rate_max_climate_tree(:) = estab_max_tree ! 'estab_max_*'; see 'stomate_constants.f90'
310	estab_rate_max_climate_grass(:) = estab_max_grass
311
312	ELSEWHERE
313
314	estab_rate_max_climate_tree(:) = zero
315	estab_rate_max_climate_grass(:) = zero
316
317	ENDWHERE
318
319	!! 2.6 Reduce maximum tree establishment rate if many trees are present.
320	! In the original DGVM, this is done using a step function which yields a
321	! reduction by factor 4 if sumfpc_wood(i) .GT. fpc_crit - 0.05.
322	! This can lead to small oscillations (without consequences however).
323	! Here, a steady linear transition is used between fpc_crit-0.075 and
324	! fpc_crit-0.025.
325	! factor(:) = 1. - 15. * ( sumfpc_wood(:) - (fpc_crit-.075))
326	! factor(:) = MAX( 0.25_r_std, MIN( 1._r_std, factor(:)))
327	! SZ modified according to Smith et al. 2001
328	! See Eq. (2) in header
329	factor(:)=(un - exp(- establish_scal_fact * (un - sumfpc_wood(:))))*(un - sumfpc_wood(:))
330	estab_rate_max_tree(:) = estab_rate_max_climate_tree(:) * factor(:)
331
332	!! 2.7 Modulate grass establishment rate.
333	! If canopy is not closed (fpc < fpc_crit-0.05), normal establishment.
334	! If canopy is closed, establishment is reduced by a factor 4.
335	! Factor is linear between these two bounds.
336	! This is different from the original DGVM where a step function is
337	! used at fpc_crit-0.05 (This can lead to small oscillations,
338	! without consequences however).
339	! factor(:) = 1. - 15. * ( sumfpc(:) - (fpc_crit-.05))
340	! factor(:) = MAX( 0.25_r_std, MIN( 1._r_std, factor(:)))
341	! estab_rate_max_grass(:) = estab_rate_max_climate_grass(:) * factor(:)
342	! SZ modified to true LPJ formulation, grasses are only allowed in the
343	! fpc fraction not occupied by trees..., 080806
344	! estab_rate_max_grass(:)=MAX(0.98-sumfpc(:),zero)
345	! See Eq. (3) in header
346	estab_rate_max_grass(:) = MAX(MIN(estab_rate_max_climate_grass(:), max_tree_coverage - sumfpc(:)),zero)
347
348	!! 2.8 Longterm grass NPP for competition between C4 and C3 grasses
349	! to avoid equal veget_max, the idea is that more reestablishment
350	! is possible for the more productive PFT
351	factor(:) = min_stomate
352	DO j = 2,nvm ! Loop over # PFTs
353	IF ( natural(j) .AND. .NOT.is_tree(j)) &
354	factor(:) = factor(:) + npp_longterm(:,j) * &
355	lm_lastyearmax(:,j) * sla(j)
356	ENDDO ! Loop over # PFTs
357
358	!! 2.9 Establish natural PFTs
359	d_ind(:,:) = zero
360
361	DO j = 2,nvm ! Loop over # PFTs
362
363	IF ( natural(j) ) THEN ! only for natural PFTs
364
365	!! 2.9.1 PFT expansion across the grid box. Not to be confused with areal coverage.
366	IF ( treat_expansion ) THEN
367
368	! only treat plants that are regenerative and present and still can expand
369	DO i = 1, npts ! Loop over # pixels - domain size
370
371	IF ( PFTpresent(i,j) .AND. &
372	( everywhere(i,j) .LT. un ) .AND. &
373	( regenerate(i,j) .GT. regenerate_crit ) ) THEN
374
375	! from how many sides is the grid box invaded (separate x and y directions
376	! because resolution may be strongly anisotropic)
377	! For the moment we only look into 4 direction but that can be expanded (JP)
378	nfrontx = 0
379	IF ( neighbours(i,3) .GT. 0 ) THEN
380	IF ( everywhere(neighbours(i,3),j) .GT. 1.-min_stomate ) nfrontx = nfrontx+1
381	ENDIF
382	IF ( neighbours(i,7) .GT. 0 ) THEN
383	IF ( everywhere(neighbours(i,7),j) .GT. 1.-min_stomate ) nfrontx = nfrontx+1
384	ENDIF
385
386	nfronty = 0
387	IF ( neighbours(i,1) .GT. 0 ) THEN
388	IF ( everywhere(neighbours(i,1),j) .GT. 1.-min_stomate ) nfronty = nfronty+1
389	ENDIF
390	IF ( neighbours(i,5) .GT. 0 ) THEN
391	IF ( everywhere(neighbours(i,5),j) .GT. 1.-min_stomate ) nfronty = nfronty+1
392	ENDIF
393
394	everywhere(i,j) = &
395	everywhere(i,j) + migrate(j) * dt/one_year * &
396	( nfrontx / resolution(i,1) + nfronty / resolution(i,2) )
397
398	IF ( .NOT. need_adjacent(i,j) ) THEN
399
400	! in that case, we also assume that the PFT expands from places within
401	! the grid box (e.g., oasis).
402	! What is this equation? No reference.
403	everywhere(i,j) = &
404	everywhere(i,j) + migrate(j) * dt/one_year * &
405	2. * SQRT( pieverywhere(i,j)/(resolution(i,1)resolution(i,2)) )
406
407	ENDIF
408
409	everywhere(i,j) = MIN( everywhere(i,j), un )
410
411	ENDIF
412
413	ENDDO ! Loop over # pixels - domain size
414
415	ENDIF ! treat expansion?
416
417	!! 2.9.2 Establishment rate
418	! - Is lower if the PFT is only present in a small part of the grid box
419	! (after its introduction), therefore multiplied by "everywhere".
420	! - Is divided by the number of PFTs that compete ("spacefight").
421	! - Is modulated by space availability (avail_tree, avail_grass).
422
423	!! 2.9.2.1 present and regenerative trees
424	IF ( is_tree(j) ) THEN
425
426	WHERE ( PFTpresent(:,j) .AND. ( regenerate(:,j) .GT. regenerate_crit ) )
427
428
429	d_ind(:,j) = estab_rate_max_tree(:)everywhere(:,j)/spacefight_tree(:) &
430	avail_tree(:) * dt/one_year
431
432	ENDWHERE
433
434	!! 2.9.2.2 present and regenerative grasses
435	ELSE
436
437	WHERE ( PFTpresent(:,j) .AND. ( regenerate(:,j) .GT. regenerate_crit ) &
438	.AND.factor(:).GT.min_stomate .AND. spacefight_grass(:).GT. min_stomate)
439
440	d_ind(:,j) = estab_rate_max_grass(:)everywhere(:,j)/spacefight_grass(:) &
441	MAX(min_stomate,npp_longterm(:,j)lm_lastyearmax(:,j)sla(j)/factor(:)) * fracnat(:) * dt/one_year
442
443	ENDWHERE
444
445	ENDIF ! tree/grass
446
447	ENDIF ! if natural
448	ENDDO ! Loop over # PFTs
449
450	!! 3. Lpj establishment in static case
451
452	! Lpj establishment in static case, SZ 080806, account for real LPJ dynamics in
453	! prescribed vegetation, i.e. population dynamics within a given area of the grid cell.
454	ELSE
455
456	d_ind(:,:) = zero
457
458	DO j = 2,nvm ! Loop over # PFTs
459
460	WHERE(ind(:,j)*cn_ind(:,j).GT.min_stomate)
461	lai_ind(:) = sla(j) * lm_lastyearmax(:,j)/(ind(:,j)*cn_ind(:,j))
462	ELSEWHERE
463	lai_ind(:) = zero
464	ENDWHERE
465
466	!! 3.1 For natural woody PFTs
467	IF ( natural(j) .AND. is_tree(j)) THEN
468
469	! See Eq. (4) in tex file.
470	fpc_nat(:,j) = MIN(un, cn_ind(:,j) * ind(:,j) * &
471	MAX( ( un - exp( - ext_coeff(j) * lai_ind(:) ) ), min_cover ) )
472
473
474	WHERE (veget_max(:,j).GT.min_stomate.AND.ind(:,j).LE.2.)
475
476	!! 3.1.1 Only establish into growing stands
477	! Only establish into growing stands, ind can become very
478	! large in the static mode because LAI is very low in poor
479	! growing conditions, favouring continuous establishment.
480	! To avoid this a maximum IND is set. BLARPP: This should be
481	! replaced by a better stand density criteria.
482	factor(:)=(un - exp(-establish_scal_fact * (un - fpc_nat(:,j))))*(un - fpc_nat(:,j))
483
484	estab_rate_max_tree(:) = estab_max_tree * factor(:)
485
486	!! 3.1.2 do establishment for natural PFTs\n
487	d_ind(:,j) = MAX( zero, estab_rate_max_tree(:) * dt/one_year)
488
489	ENDWHERE
490
491	!SZ: quickfix: to simulate even aged stand, uncomment the following lines...
492	!where (ind(:,j) .LE. min_stomate)
493	!d_ind(:,j) = 0.1 !MAX( 0.0, estab_rate_max_tree(:) * dt/one_year)
494	WHERE (veget_max(:,j).GT.min_stomate .AND. ind(:,j).EQ.zero)
495	d_ind(:,j) = ind_0_estab
496	ENDWHERE
497
498	!! 3.2 For natural grass PFTs
499	ELSEIF ( natural(j) .AND. .NOT.is_tree(j)) THEN
500
501	WHERE (veget_max(:,j).GT.min_stomate)
502
503	fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) * &
504	MAX( ( un - exp( - ext_coeff(j) * lai_ind(:) ) ), min_cover )
505
506	d_ind(:,j) = MAX(zero , (un - fpc_nat(:,j)) * dt/one_year )
507
508	ENDWHERE
509
510	WHERE (veget_max(:,j).GT.min_stomate .AND. ind(:,j).EQ. zero)
511	d_ind(:,j) = ind_0_estab
512	ENDWHERE
513
514	ENDIF
515
516	ENDDO ! Loop over # PFTs
517
518	ENDIF ! DGVM OR NOT
519
520	!! 4. Biomass calculation
521
522	DO j = 2,nvm ! Loop over # PFTs
523
524	IF ( natural(j) ) THEN ! only for natural PFTs
525
526	!! 4.1 Herbivores reduce establishment rate
527	! We suppose that saplings are vulnerable during a given time after establishment.
528	! This is taken into account by preventively reducing the establishment rate.
529	IF ( control%ok_herbivory ) THEN
530
531	d_ind(:,j) = d_ind(:,j) * EXP( - tau_eatup/herbivores(:,j) )
532
533	ENDIF
534
535	!! 4.2 Total biomass.
536	! Add biomass only if d_ind, over one year, is of the order of ind.
537	! save old leaf mass to calculate leaf age
538	leaf_mass_young(:) = leaf_frac(:,j,1) * biomass(:,j,ileaf,icarbon)
539
540	! total biomass of existing PFT to limit biomass added from establishment
541	total_bm_c(:) = zero
542
543	DO k = 1, nparts
544	total_bm_c(:) = total_bm_c(:) + biomass(:,j,k,icarbon)
545	ENDDO
546	IF(control%ok_dgvm) THEN
547	vn(:) = veget_max(:,j)
548	ELSE
549	vn(:) = un
550	ENDIF
551
552	total_bm_sapl(:,:)=zero
553	DO k = 1, nparts
554	WHERE(d_ind(:,j).GT.min_stomate.AND.vn(:).GT.min_stomate)
555
556	total_bm_sapl(:,icarbon) = total_bm_sapl(:,icarbon) + &
557	bm_sapl_old(j,k,icarbon) * d_ind(:,j) / vn(:)
558
559	ENDWHERE
560	ENDDO
561
562	!! 4.3 Woodmass calculation
563
564	!! 4.3.1 with DGVM
565	IF(control%ok_dgvm) THEN
566
567	! SZ calculate new woodmass_ind and veget_max after establishment (needed for correct scaling!)
568	! essential correction for MERGE!
569	IF (is_tree(j)) THEN
570	DO i = 1,npts ! Loop over # pixels - domain size
571	IF((d_ind(i,j)+ind(i,j)).GT.min_stomate) THEN
572
573	IF((total_bm_c(i).LE.min_stomate) .OR. (veget_max(i,j) .LE. min_stomate)) THEN
574
575	! new wood mass of PFT
576	woodmass_ind(i,j) = &
577	(((biomass(i,j,isapabove,icarbon) + biomass(i,j,isapbelow,icarbon) &
578	+ biomass(i,j,iheartabove,icarbon) + biomass(i,j,iheartbelow,icarbon))*veget_max(i,j)) &
579	+ ((bm_sapl_old(j,isapabove,icarbon)) + (bm_sapl_old(j,isapbelow,icarbon)) &
580	+ (bm_sapl_old(j,iheartabove,icarbon)) &
581	+ (bm_sapl_old(j,iheartbelow,icarbon)))*d_ind(i,j))/(ind(i,j) + d_ind(i,j))
582
583	ELSE
584
585	! new biomass is added to the labile pool, hence there is no change
586	! in CA associated with establishment
587	woodmass_ind(i,j) = &
588	& (biomass(i,j,isapabove,icarbon) + biomass(i,j,isapbelow,icarbon) &
589	& +biomass(i,j,iheartabove,icarbon) + biomass(i,j,iheartbelow,icarbon))*veget_max(i,j) &
590	& /(ind(i,j) + d_ind(i,j))
591
592	ENDIF
593
594	! new diameter of PFT
595	dia(i) = (woodmass_ind(i,j)/(pipe_density(j)pi/4.pipe_tune2(j))) &
596	& **(1./(2.+pipe_tune3(j)))
597	vn(i) = (ind(i,j) + d_ind(i,j))pipe_tune1(j)MIN(dia(i),maxdia(j))**pipe_tune_exp_coeff(j)
598
599	ENDIF
600	ENDDO ! Loop over # pixels - domain size
601	ELSE ! for grasses, cnd=1, so the above calculation cancels
602	vn(:) = ind(:,j) + d_ind(:,j)
603	ENDIF
604
605	!! 4.3.2 without DGVM (static)\n
606	ELSE
607	DO i=1,npts ! Loop over # pixels - domain size
608	IF(is_tree(j).AND.(d_ind(i,j)+ind(i,j)).GT.min_stomate) THEN
609	IF(total_bm_c(i).LE.min_stomate) THEN
610
611	! new wood mass of PFT
612	woodmass_ind(i,j) = &
613	& (((biomass(i,j,isapabove,icarbon) + biomass(i,j,isapbelow,icarbon) &
614	& + biomass(i,j,iheartabove,icarbon) + biomass(i,j,iheartbelow,icarbon))) &
615	& + ((bm_sapl_old(j,isapabove,icarbon)) + &
616	(bm_sapl_old(j,isapbelow,icarbon)) &
617	& + (bm_sapl_old(j,iheartabove,icarbon)) + &
618	(bm_sapl_old(j,iheartbelow,icarbon)))*d_ind(i,j))/(ind(i,j)+d_ind(i,j))
619
620	ELSE
621
622	! new biomass is added to the labile pool, hence there is no change
623	! in CA associated with establishment
624	woodmass_ind(i,j) = &
625	& (biomass(i,j,isapabove,icarbon) + biomass(i,j,isapbelow,icarbon) &
626	& + biomass(i,j,iheartabove,icarbon) + biomass(i,j,iheartbelow,icarbon)) &
627	& /(ind(i,j) + d_ind(i,j))
628
629	ENDIF
630	ENDIF
631	ENDDO ! Loop over # pixels - domain size
632
633	vn(:) = un ! cannot change in static!, and veget_max implicit in d_ind
634
635	ENDIF
636
637	!! 4.4 total biomass of PFT added by establishment defined over veget_max ...
638	total_bm_sapl(:,:) = zero
639	DO k = 1, nparts ! Loop over # litter tissues (nparts=8); see 'stomate_constants.f90'
640	WHERE(d_ind(:,j).GT.min_stomate.AND.total_bm_c(:).GT.min_stomate.AND.vn(:).GT.min_stomate)
641
642	total_bm_sapl(:,icarbon) = total_bm_sapl(:,icarbon) + &
643	(bm_sapl_old(j,k,icarbon)) * d_ind(:,j) / vn(:)
644	ENDWHERE
645	ENDDO ! Loop over # litter tissues
646
647	!! 4.5 Update biomass at each component
648	DO k = 1, nparts ! Loop over # litter tissues
649
650	bm_new(:) = zero
651
652	! first ever establishment, C flows
653	WHERE( d_ind(:,j).GT.min_stomate .AND. &
654	total_bm_c(:).LE.min_stomate .AND. &
655	vn(:).GT.min_stomate)
656	! WHERE ( many_new(:) )
657
658	!bm_new(:) = d_ind(:,j) * bm_sapl(j,k) / veget_max (:,j)
659	bm_new(:) = d_ind(:,j) * (bm_sapl_old(j,k,icarbon)) / vn(:)
660
661	biomass(:,j,k,icarbon) = biomass(:,j,k,icarbon) + bm_new(:)
662
663	co2_to_bm(:,j) = co2_to_bm(:,j) + bm_new(:) / dt
664
665	ENDWHERE
666
667	! establishment into existing population, C flows
668	WHERE(d_ind(:,j).GT.min_stomate.AND.total_bm_c(:).GT.min_stomate)
669
670	bm_new(:) = total_bm_sapl(:,icarbon) * biomass(:,j,k,icarbon) / total_bm_c(:)
671
672	biomass(:,j,k,icarbon) = biomass(:,j,k,icarbon) + bm_new(:)
673	co2_to_bm(:,j) = co2_to_bm(:,j) + bm_new(:) / dt
674
675	ENDWHERE
676	ENDDO ! Loop over # litter tissues
677
678
679	!! 4.6 Decrease leaf age in youngest class if new leaf biomass is higher than old one.
680	WHERE ( d_ind(:,j) * bm_sapl_old(j,ileaf,icarbon) .GT. min_stomate )
681
682	! reset leaf ages. Should do a real calculation like in the npp routine,
683	! but this case is rare and not worth messing around.
684	! SZ 080806, added real calculation now, because otherwise leaf_age/leaf_frac
685	! are not initialised for the calculation of vmax, and hence no growth at all.
686	! logic follows that of stomate_npp.f90, just that it's been adjusted for the code here
687	leaf_age(:,j,1) = leaf_age(:,j,1) * leaf_mass_young(:) / &
688	( leaf_mass_young(:) + d_ind(:,j) * bm_sapl_old(j,ileaf,icarbon) )
689
690	ENDWHERE
691
692	leaf_mass_young(:) = leaf_mass_young(:) + d_ind(:,j) * bm_sapl_old(j,ileaf,icarbon)
693
694	!! 4.7 Youngest class: new mass in youngest class divided by total new mass
695	WHERE ( biomass(:,j,ileaf,icarbon) .GT. min_stomate )
696	! new age class fractions (fraction in youngest class increases)
697	leaf_frac(:,j,1) = leaf_mass_young(:) / biomass(:,j,ileaf,icarbon)
698
699	ENDWHERE
700
701	!! 4.8 Other classes: old mass in leaf age class divided by new mass
702	DO m = 2, nleafages
703
704	WHERE ( biomass(:,j,ileaf,icarbon) .GT. min_stomate )
705
706	leaf_frac(:,j,m) = leaf_frac(:,j,m) * &
707	( biomass(:,j,ileaf,icarbon) + d_ind(:,j) * &
708	bm_sapl_old(j,ileaf,icarbon) ) / biomass(:,j,ileaf,icarbon)
709
710	ENDWHERE
711
712	ENDDO
713
714	!! 4.9 Update age and number of individuals
715	WHERE ( d_ind(:,j) .GT. min_stomate )
716
717	age(:,j) = age(:,j) * ind(:,j) / ( ind(:,j) + d_ind(:,j) )
718
719	ind(:,j) = ind(:,j) + d_ind(:,j)
720
721	ENDWHERE
722
723	!! 4.10 Convert excess sapwood to heartwood
724	!! No longer done : supressed by S. Zaehle given that the LPJ logic of carbon allocation was
725	!! contradictory to SLAVE allocation. See CVS tag 1_5 for initial formulation.
726
727
728	ENDIF ! natural
729
730	ENDDO ! Loop over # PFTs
731
732	!! 5. history
733
734	d_ind = d_ind / dt
735
736	CALL xios_orchidee_send_field("IND_ESTAB",d_ind)
737	CALL xios_orchidee_send_field("ESTABTREE",estab_rate_max_tree)
738	CALL xios_orchidee_send_field("ESTABGRASS",estab_rate_max_grass)
739
740	CALL histwrite_p (hist_id_stomate, 'IND_ESTAB', itime, d_ind, npts*nvm, horipft_index)
741	CALL histwrite_p (hist_id_stomate, 'ESTABTREE', itime, estab_rate_max_tree, npts, hori_index)
742	CALL histwrite_p (hist_id_stomate, 'ESTABGRASS', itime, estab_rate_max_grass, npts, hori_index)
743
744	IF (bavard.GE.4) WRITE(numout,*) 'Leaving establish'
745
746	END SUBROUTINE establish_diagnostic
747
748	END MODULE lpj_establish

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

Context Navigation

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

Download in other formats: