Context Navigation

lpj_establish.f90 @ 6737

Last change on this file since 6737 was 4719, checked in by albert.jornet, 7 years ago

Merge: from revisions [4491:4695/trunk/ORCHIDEE]

Merge done in [4671:4718/perso/albert.jornet/MICT_MERGE]

Property svn:keywords set to HeadURL Date Author Revision

File size: 42.0 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_stomate/lpj_establish.f90 @ 6737

Download in other formats: