Context Navigation

lpj_establish.f90 @ 7746

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

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

Context Navigation

source: branches/publications/ORCHIDEE-MICT-OP-r6850/src_stomate/lpj_establish.f90 @ 7746

Download in other formats: