Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-261/src_stomate/lpj_establish.f90 @ 7474

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

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

Download in other formats: