Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lpj_establish.f90 @ 5491

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

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

Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_stomate/lpj_establish.f90 @ 5491

Download in other formats: