Context Navigation

source: branches/publications/ORCHIDEE-GMv3.2/ORCHIDEE/src_stomate/lpj_establish.f90 @ 5816

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

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

Download in other formats: