Context Navigation

source: branches/publications/ORCHIDEE-Clateral/src_stomate/lpj_establish.f90 @ 7346

Last change on this file since 7346 was 2944, checked in by josefine.ghattas, 9 years ago

Ticket #182

Cleaning after commit update of DGVM.

Calculated veget_max_tree and nbtree to make independent of the number of pfts.
Added new parameters used instead of numbers in the code.
Added more comments.

Property svn:keywords set to HeadURL Date Author Revision

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

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

Download in other formats: