Context Navigation

source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_stomate/lpj_establish.f90 @ 8398

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

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

Download in other formats: