Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lpj_establish.f90 @ 7541

Last change on this file since 7541 was 7541, checked in by fabienne.maignan, 2 years ago
Zhang publication on coupling factor
File size: 40.0 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE_2.2_r7266/ORCHIDEE/src_stomate/lpj_establish.f90 @ 7541

Download in other formats: