Context Navigation

source: tags/ORCHIDEE_1_9_6/ORCHIDEE/src_stomate/lpj_establish.f90 @ 880

Last change on this file since 880 was 837, checked in by didier.solyga, 12 years ago
Add a test to prevent error for the choice of the Metaclasses. Correct fVegLitter ipcc output. Replace the last 1 by un.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 31.2 KB

Line
1	! establishment routine
2	! Suppose seed pool >> establishment rate.
3	!
4	!< $HeadURL$
5	!< $Date$
6	!< $Author$
7	!< $Revision$
8	! IPSL (2006)
9	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
10	!
11	MODULE lpj_establish
12
13	! modules used:
14
15	USE ioipsl
16	USE stomate_data
17	USE constantes
18
19	IMPLICIT NONE
20
21	! private & public routines
22
23	PRIVATE
24	PUBLIC establish,establish_clear
25
26	! first call
27	LOGICAL, SAVE :: firstcall = .TRUE.
28	CONTAINS
29
30
31	SUBROUTINE establish_clear
32	firstcall = .TRUE.
33	END SUBROUTINE establish_clear
34
35	SUBROUTINE establish (npts, dt, PFTpresent, regenerate, &
36	neighbours, resolution, need_adjacent, herbivores, &
37	precip_annual, gdd0, lm_lastyearmax, &
38	cn_ind, lai, avail_tree, avail_grass, npp_longterm, &
39	leaf_age, leaf_frac, &
40	ind, biomass, age, everywhere, co2_to_bm,veget_max, woodmass_ind)
41	!
42	! 0 declarations
43	!
44
45	! 0.1 input
46
47	! Domain size
48	INTEGER(i_std), INTENT(in) :: npts
49	! Time step of vegetation dynamics (days)
50	REAL(r_std), INTENT(in) :: dt
51	! Is pft there
52	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: PFTpresent
53	! Winter sufficiently cold
54	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: regenerate
55	! indices of the 8 neighbours of each grid point (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
56	INTEGER(i_std), DIMENSION(npts,8), INTENT(in) :: neighbours
57	! resolution at each grid point in m (1=E-W, 2=N-S)
58	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: resolution
59	! in order for this PFT to be introduced, does it have to be present in an
60	! adjacent grid box?
61	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: need_adjacent
62	! time constant of probability of a leaf to be eaten by a herbivore (days)
63	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: herbivores
64	! annual precipitation (mm/year)
65	REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_annual
66	! growing degree days (C)
67	REAL(r_std), DIMENSION(npts), INTENT(in) :: gdd0
68	! last year's maximum leaf mass, for each PFT (gC/(m**2 of ground))
69	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lm_lastyearmax
70	! crown area of individuals (m**2)
71	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: cn_ind
72	! leaf area index OF AN INDIVIDUAL PLANT
73	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lai
74	! space availability for trees
75	REAL(r_std), DIMENSION(npts), INTENT(in) :: avail_tree
76	! space availability for grasses
77	REAL(r_std), DIMENSION(npts), INTENT(in) :: avail_grass
78	! longterm NPP, for each PFT (gC/(m**2 of ground))
79	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: npp_longterm
80	! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground
81	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max
82
83	! 0.2 modified fields
84
85	! leaf age (days)
86	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_age
87	! fraction of leaves in leaf age class
88	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac
89	! Number of individuals / m2
90	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind
91	! biomass (gC/(m**2 of ground))
92	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout) :: biomass
93	! mean age (years)
94	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age
95	! is the PFT everywhere in the grid box or very localized (after its introduction)
96	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: everywhere
97	! biomass uptaken (gC/(m**2 of total ground)/day)
98	!NV passage 2D
99	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: co2_to_bm
100	! woodmass of the individual, needed to calculate crownarea in lpj_crownarea
101	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: woodmass_ind
102
103	! 0.3 local
104
105	! time during which a sapling can be entirely eaten by herbivores (d)
106	REAL(r_std) :: tau_eatup
107	! new fpc ( foliage protected cover: fractional coverage )
108	REAL(r_std), DIMENSION(npts,nvm) :: fpc_nat
109	! maximum tree establishment rate, based on climate only
110	REAL(r_std), DIMENSION(npts) :: estab_rate_max_climate_tree
111	! maximum grass establishment rate, based on climate only
112	REAL(r_std), DIMENSION(npts) :: estab_rate_max_climate_grass
113	! maximum tree establishment rate, based on climate and fpc
114	REAL(r_std), DIMENSION(npts) :: estab_rate_max_tree
115	! maximum grass establishment rate, based on climate and fpc
116	REAL(r_std), DIMENSION(npts) :: estab_rate_max_grass
117	! total natural fpc
118	REAL(r_std), DIMENSION(npts) :: sumfpc
119	! total fraction occupied by natural vegetation
120	REAL(r_std), DIMENSION(npts) :: fracnat
121	! total woody fpc
122	REAL(r_std), DIMENSION(npts) :: sumfpc_wood
123	! for trees, measures the total concurrence for available space
124	REAL(r_std), DIMENSION(npts) :: spacefight_tree
125	! for grasses, measures the total concurrence for available space
126	REAL(r_std), DIMENSION(npts) :: spacefight_grass
127	! change in number of individuals /m2 per time step (per day in history file)
128	REAL(r_std), DIMENSION(npts,nvm) :: d_ind
129	! biomass increase (gC/(m**2 of ground))
130	REAL(r_std), DIMENSION(npts) :: bm_new
131	! stem diameter (m)
132	REAL(r_std), DIMENSION(npts) :: dia
133	! temporary variable
134	REAL(r_std), DIMENSION(npts) :: b1
135	! new sap mass (gC/(m**2 of ground))
136	REAL(r_std), DIMENSION(npts) :: sm2
137	! woodmass of an individual
138	REAL(r_std), DIMENSION(npts) :: woodmass
139	! carbon mass in youngest leaf age class (gC/m**2 PFT)
140	REAL(r_std), DIMENSION(npts) :: leaf_mass_young
141	! ratio of hw(above) to total hw, sm(above) to total sm
142	REAL(r_std), DIMENSION(npts) :: sm_at
143	! reduction factor for establishment if many trees or grasses are present
144	REAL(r_std), DIMENSION(npts) :: factor
145	! Total carbon mass for all pools
146	REAL(r_std), DIMENSION(npts) :: total_bm_c
147	! Total sappling biomass for all pools
148	REAL(r_std), DIMENSION(npts) :: total_bm_sapl
149	! from how many sides is the grid box invaded
150	INTEGER(i_std) :: nfrontx
151	INTEGER(i_std) :: nfronty
152	! daily establishment rate is large compared to present number of individuals
153	!LOGICAL, DIMENSION(npts) :: many_new
154	! flow due to new individuals
155	! veget_max after establishment, to get a proper estimate of carbon and nitrogen
156	REAL(r_std), DIMENSION(npts) :: vn
157	! lai on each PFT surface
158	REAL(r_std), DIMENSION(npts) :: lai_ind
159
160	! indices
161	INTEGER(i_std) :: i,j,k,m
162
163	! =========================================================================
164
165	IF (bavard.GE.3) WRITE(numout,*) 'Entering establish'
166
167	!
168	! 1 messages and initialization
169	!
170	tau_eatup = one_year/2.
171
172	IF ( firstcall ) THEN
173
174	WRITE(numout,*) 'establish:'
175
176	WRITE(numout,*) ' > time during which a sapling can be entirely eaten by herbivores (d): ', &
177	tau_eatup
178
179	firstcall = .FALSE.
180
181	ENDIF
182
183
184	IF (control%ok_dgvm) THEN
185	!
186	! 2 recalculate fpc
187	!
188
189	!
190	! 2.1 Only natural part of the grid cell
191	!
192
193	fracnat(:) = un
194	do j = 2,nvm
195	IF ( .NOT. natural(j) ) THEN
196	fracnat(:) = fracnat(:) - veget_max(:,j)
197	ENDIF
198	ENDDO
199
200	!
201	! 2.2 total natural fpc on grid
202	!
203	sumfpc(:) = zero
204	DO j = 2,nvm
205
206	IF ( natural(j) ) THEN
207	WHERE(fracnat(:).GT.min_stomate)
208	WHERE (lai(:,j) == val_exp)
209	fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) / fracnat(:)
210	ELSEWHERE
211	fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) / fracnat(:) &
212	* ( un - exp( - lm_lastyearmax(:,j) * sla(j) * ext_coeff(j) ) )
213	ENDWHERE
214	ENDWHERE
215
216	WHERE ( PFTpresent(:,j) )
217	sumfpc(:) = sumfpc(:) + fpc_nat(:,j)
218	ENDWHERE
219	ELSE
220
221	fpc_nat(:,j) = zero
222
223	ENDIF
224
225	ENDDO
226
227	!
228	! 2.3 total woody fpc on grid and number of regenerative tree pfts
229	!
230
231	sumfpc_wood(:) = zero
232	spacefight_tree(:) = zero
233
234	DO j = 2,nvm
235
236	IF ( tree(j) .AND. natural(j) ) THEN
237
238	! total woody fpc
239
240	WHERE ( PFTpresent(:,j) )
241	sumfpc_wood(:) = sumfpc_wood(:) + fpc_nat(:,j)
242	ENDWHERE
243
244	! how many trees are competing? Count a PFT fully only if it is present
245	! on the whole grid box.
246
247	WHERE ( PFTpresent(:,j) .AND. ( regenerate(:,j) .GT. regenerate_crit ) )
248	spacefight_tree(:) = spacefight_tree(:) + everywhere(:,j)
249	ENDWHERE
250
251	ENDIF
252
253	ENDDO
254
255	!
256	! 2.4 number of natural grasses
257	!
258
259	spacefight_grass(:) = zero
260
261	DO j = 2,nvm
262
263	IF ( .NOT. tree(j) .AND. natural(j) ) THEN
264
265	! how many grasses are competing? Count a PFT fully only if it is present
266	! on the whole grid box.
267
268	WHERE ( PFTpresent(:,j) )
269	spacefight_grass(:) = spacefight_grass(:) + everywhere(:,j)
270	ENDWHERE
271
272	ENDIF
273
274	ENDDO
275
276	!
277	! 3 establishment rate
278	!
279
280	!
281	! 3.1 maximum establishment rate, based on climate only
282	!
283
284	WHERE ( ( precip_annual(:) .GE. precip_crit ) .AND. ( gdd0(:) .GE. gdd_crit_estab ) )
285
286	estab_rate_max_climate_tree(:) = estab_max_tree
287	estab_rate_max_climate_grass(:) = estab_max_grass
288
289	ELSEWHERE
290
291	estab_rate_max_climate_tree(:) = zero
292	estab_rate_max_climate_grass(:) = zero
293
294	ENDWHERE
295
296	!
297	! 3.2 reduce maximum tree establishment rate if many trees present.
298	! In the original DGVM, this is done using a step function which yields a
299	! reduction by factor 4 if sumfpc_wood(i) .GT. fpc_crit - 0.05.
300	! This can lead to small oscillations (without consequences however).
301	! Here, a steady linear transition is used between fpc_crit-0.075 and
302	! fpc_crit-0.025.
303	!
304
305	! factor(:) = 1. - 15. * ( sumfpc_wood(:) - (fpc_crit - 0.075) )
306	! factor(:) = MAX( 0.25_r_std, MIN( un, factor(:) ) )
307
308	!SZ modified according to Smith et al. 2001, 080806
309	factor(:)=(un - exp(- establish_scal_fact * (un - sumfpc_wood(:))))*(un - sumfpc_wood(:))
310
311	estab_rate_max_tree(:) = estab_rate_max_climate_tree(:) * factor(:)
312
313	!
314	! 3.3 Modulate grass establishment rate.
315	! If canopy is not closed (fpc < fpc_crit-0.05), normal establishment.
316	! If canopy is closed, establishment is reduced by a factor 4.
317	! Factor is linear between these two bounds.
318	! This is different from the original DGVM where a step function is
319	! used at fpc_crit-0.05 (This can lead to small oscillations,
320	! without consequences however).
321	!
322
323	! factor(:) = 1. - establish_scal_fact * ( sumfpc(:) - (fpc_crit - 0.05) )
324	! factor(:) = MAX( 0.25_r_std, MIN( un, factor(:) ) )
325	! estab_rate_max_grass(:) = estab_rate_max_climate_grass(:) * factor(:)
326
327	!SZ modified to true LPJ formulation, grasses are only allowed in the
328	!fpc fraction not occupied by trees..., 080806
329	!NVmodif estab_rate_max_grass(:)=MAX(0.98-sumfpc(:),zero)
330	estab_rate_max_grass(:) = MAX(MIN(estab_rate_max_climate_grass(:), max_tree_coverage - sumfpc(:)),zero)
331
332	! SZ: longterm grass NPP for competition between C4 and C3 grasses
333	! to avoid equal veget_max, the idea is that more reestablishment
334	! is possible for the more productive PFT
335	factor(:) = min_stomate
336	DO j = 2,nvm
337	IF ( natural(j) .AND. .NOT.tree(j)) &
338	factor(:) = factor(:) + npp_longterm(:,j) * &
339	lm_lastyearmax(:,j) * sla(j)
340	ENDDO
341	!
342	!
343	!
344	! 4 do establishment for natural PFTs
345	!
346
347	d_ind(:,:) = zero
348
349	DO j = 2,nvm
350
351	! only for natural PFTs
352
353	IF ( natural(j) ) THEN
354
355	!
356	! 4.1 PFT expansion across the grid box. Not to be confused with areal
357	! coverage.
358	!
359
360	IF ( treat_expansion ) THEN
361
362	! only treat plants that are regenerative and present and still can expand
363
364	DO i = 1, npts
365
366	IF ( PFTpresent(i,j) .AND. &
367	( everywhere(i,j) .LT. un ) .AND. &
368	( regenerate(i,j) .GT. regenerate_crit ) ) THEN
369
370	! from how many sides is the grid box invaded (separate x and y directions
371	! because resolution may be strongly anisotropic)
372	!
373	! For the moment we only look into 4 direction but that can be extanded (JP)
374	!
375	nfrontx = 0
376	IF ( neighbours(i,3) .GT. 0 ) THEN
377	IF ( everywhere(neighbours(i,3),j) .GT. 1.-min_stomate ) nfrontx = nfrontx+1
378	ENDIF
379	IF ( neighbours(i,7) .GT. 0 ) THEN
380	IF ( everywhere(neighbours(i,7),j) .GT. 1.-min_stomate ) nfrontx = nfrontx+1
381	ENDIF
382
383	nfronty = 0
384	IF ( neighbours(i,1) .GT. 0 ) THEN
385	IF ( everywhere(neighbours(i,1),j) .GT. 1.-min_stomate ) nfronty = nfronty+1
386	ENDIF
387	IF ( neighbours(i,5) .GT. 0 ) THEN
388	IF ( everywhere(neighbours(i,5),j) .GT. 1.-min_stomate ) nfronty = nfronty+1
389	ENDIF
390
391	everywhere(i,j) = &
392	everywhere(i,j) + migrate(j) * dt/one_year * &
393	( nfrontx / resolution(i,1) + nfronty / resolution(i,2) )
394
395	IF ( .NOT. need_adjacent(i,j) ) THEN
396
397	! in that case, we also assume that the PFT expands from places within
398	! the grid box (e.g., oasis).
399
400	everywhere(i,j) = &
401	everywhere(i,j) + migrate(j) * dt/one_year * &
402	2. * SQRT( pieverywhere(i,j)/(resolution(i,1)resolution(i,2)) )
403
404	ENDIF
405
406	everywhere(i,j) = MIN( everywhere(i,j), un )
407
408	ENDIF
409
410	ENDDO
411
412	ENDIF ! treat expansion?
413
414	!
415	! 4.2 establishment rate
416	! - Is lower if the PFT is only present in a small part of the grid box
417	! (after its introduction), therefore multiplied by "everywhere".
418	! - Is divided by the number of PFTs that compete ("spacefight").
419	! - Is modulated by space availability (avail_tree, avail_grass).
420	!
421
422	IF ( tree(j) ) THEN
423
424	! 4.2.1 present and regenerative trees
425
426	WHERE ( PFTpresent(:,j) .AND. ( regenerate(:,j) .GT. regenerate_crit ) )
427
428
429	d_ind(:,j) = estab_rate_max_tree(:)everywhere(:,j)/spacefight_tree(:) &
430	avail_tree(:) * dt/one_year
431
432	ENDWHERE
433
434	ELSE
435
436	! 4.2.2 present and regenerative grasses
437
438	WHERE ( PFTpresent(:,j) .AND. ( regenerate(:,j) .GT. regenerate_crit ) &
439	.AND.factor(:).GT.min_stomate .AND. spacefight_grass(:).GT. min_stomate)
440
441	d_ind(:,j) = estab_rate_max_grass(:)everywhere(:,j)/spacefight_grass(:) &
442	MAX(min_stomate,npp_longterm(:,j)lm_lastyearmax(:,j)sla(j)/factor(:)) * fracnat(:) * dt/one_year
443
444	ENDWHERE
445
446	ENDIF ! tree/grass
447
448	ENDIF ! if natural
449	ENDDO ! PFTs
450
451	ELSE ! lpj establishment in static case, SZ 080806, account for real LPJ dynamics in
452	! prescribed vegetation, i.e. population dynamics within a given area of the
453	! grid cell
454
455	d_ind(:,:) = zero
456
457	DO j = 2,nvm
458
459	! only for natural PFTs
460
461	WHERE(ind(:,j)*cn_ind(:,j).GT.min_stomate)
462	lai_ind(:) = sla(j) * lm_lastyearmax(:,j)/(ind(:,j)*cn_ind(:,j))
463	ELSEWHERE
464	lai_ind(:) = zero
465	ENDWHERE
466
467	IF ( natural(j) .AND. tree(j)) THEN
468
469	fpc_nat(:,j) = MIN(un, cn_ind(:,j) * ind(:,j) * &
470	MAX( ( un - exp( - ext_coeff(j) * lai_ind(:) ) ), min_cover ) )
471	!fpc_nat(:,j) = max(fpc_nat(:,j),1.-exp(-0.5sla(j) lm_lastyearmax(:,j)))
472
473
474	WHERE (veget_max(:,j).GT.min_stomate.AND.ind(:,j).LE.2.)
475
476	! only establish into growing stands, ind can become very
477	! large in the static mode because LAI is very low in poor
478	! growing conditions, favouring continuous establishment. To avoid this
479	! a maximum IND is set. BLARPP: This should be replaced by a
480	! better stand density criteria
481	!
482	factor(:)=(un - exp(-establish_scal_fact * (un - fpc_nat(:,j))))*(un - fpc_nat(:,j))
483
484	estab_rate_max_tree(:) = estab_max_tree * factor(:)
485	!
486	! 4 do establishment for natural PFTs
487	!
488	d_ind(:,j) = MAX( zero, estab_rate_max_tree(:) * dt/one_year)
489
490	ENDWHERE
491
492	!SZ: quickfix: to simulate even aged stand, uncomment the following lines...
493	!where (ind(:,j) .LE. min_stomate)
494	!d_ind(:,j) = 0.1 !MAX( 0.0, estab_rate_max_tree(:) * dt/one_year)
495
496	WHERE (veget_max(:,j).GT.min_stomate .AND. ind(:,j).EQ.zero)
497	d_ind(:,j) = ind_0_estab
498	! elsewhere
499	!d_ind(:,j) =0.0
500	ENDWHERE
501
502	ELSEIF ( natural(j) .AND. .NOT.tree(j)) THEN
503
504	WHERE (veget_max(:,j).GT.min_stomate)
505
506	fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) * &
507	MAX( ( un - exp( - ext_coeff(j) * lai_ind(:) ) ), min_cover )
508
509	d_ind(:,j) = MAX(zero , (un - fpc_nat(:,j)) * dt/one_year )
510
511	ENDWHERE
512
513	WHERE (veget_max(:,j).GT.min_stomate .AND. ind(:,j).EQ. zero)
514	d_ind(:,j) = ind_0_estab
515	ENDWHERE
516
517	ENDIF
518
519	ENDDO
520
521	ENDIF ! DGVM OR NOT
522
523	DO j = 2,nvm
524
525	! only for natural PFTs
526
527	IF ( natural(j) ) THEN
528
529	!
530	! 4.3 herbivores reduce establishment rate
531	! We suppose that saplings are vulnerable during a given time after establishment.
532	! This is taken into account by preventively reducing the establishment rate.
533	!
534
535	IF ( ok_herbivores ) THEN
536
537	d_ind(:,j) = d_ind(:,j) * EXP( - tau_eatup/herbivores(:,j) )
538
539	ENDIF
540
541	!
542	! 4.4 be sure that ind*cn_ind does not exceed 1
543	!SZ This control is now moved to lpj_cover.f90
544	!SZ
545
546	!The aim is to control for sum(veget)=1., irrespective of ind*cnd (crowns can overlap as long as
547	! there is enough light
548	!
549	!SZ: This could be part of the dynamic vegetation problem of Orchidee
550	!in conjunction with the wrong formulation of establishment response
551	!to tree fpc above...
552	! WHERE ( ( d_ind(:,j) .GT. zero ) .AND. &
553	! ( (ind(:,j)+d_ind(:,j))*cn_ind(:,j) .GT. un ) )
554	!
555	! d_ind(:,j) = MAX( 1._stnd / cn_ind(:,j) - ind(:,j), zero )
556	!
557	! ENDWHERE
558
559	!
560	! 4.5 new properties where there is establishment ( d_ind > 0 )
561	!
562
563	! 4.5.1 biomass.
564	! Add biomass only if d_ind, over one year, is of the order of ind.
565	! (If we don't do this, the biomass density can become very low).
566	! In that case, take biomass from the atmosphere.
567
568	! compare establishment rate and present number of inidivuals
569	!many_new(:) = ( d_ind(:,j) .GT. 0.1 * ind(:,j) )
570
571	! gives a better vectorization of the VPP
572
573	!IF ( ANY( many_new(:) ) ) THEN
574
575	! save old leaf mass to calculate leaf age
576	leaf_mass_young(:) = leaf_frac(:,j,1) * biomass(:,j,ileaf)
577	! total biomass of existing PFT to limit biomass added from establishment
578	total_bm_c(:) = zero
579
580	DO k = 1, nparts
581	total_bm_c(:)=total_bm_c(:)+biomass(:,j,k)
582	ENDDO
583	IF(control%ok_dgvm) THEN
584	vn(:) = veget_max(:,j)
585	ELSE
586	vn(:) = un
587	ENDIF
588	total_bm_sapl(:)=zero
589	DO k = 1, nparts
590	WHERE(d_ind(:,j).GT.min_stomate.AND.vn(:).GT.min_stomate)
591
592	total_bm_sapl(:) = total_bm_sapl(:) + &
593	bm_sapl(j,k) * d_ind(:,j) / vn(:)
594	ENDWHERE
595	ENDDO
596
597	IF(control%ok_dgvm) THEN
598	! SZ calculate new woodmass_ind and veget_max after establishment (needed for correct scaling!)
599	! essential correction for MERGE!
600	IF(tree(j))THEN
601	DO i=1,npts
602	IF((d_ind(i,j)+ind(i,j)).GT.min_stomate) THEN
603
604	IF((total_bm_c(i).LE.min_stomate) .OR. (veget_max(i,j) .LE. min_stomate)) THEN
605
606	! new wood mass of PFT
607	woodmass_ind(i,j) = &
608	& (((biomass(i,j,isapabove) + biomass(i,j,isapbelow) &
609	& + biomass(i,j,iheartabove) + biomass(i,j,iheartbelow))*veget_max(i,j)) &
610	& + (bm_sapl(j,isapabove) + bm_sapl(j,isapbelow) &
611	& + bm_sapl(j,iheartabove) + bm_sapl(j,iheartbelow))*d_ind(i,j))/(ind(i,j) + d_ind(i,j))
612
613	ELSE
614	! new biomass is added to the labile pool, hence there is no change in CA associated with establishment
615	woodmass_ind(i,j) = &
616	& (biomass(i,j,isapabove) + biomass(i,j,isapbelow) &
617	& +biomass(i,j,iheartabove) + biomass(i,j,iheartbelow))*veget_max(i,j) &
618	& /(ind(i,j) + d_ind(i,j))
619
620	ENDIF
621
622	! new diameter of PFT
623	dia(i) = (woodmass_ind(i,j)/(pipe_densitypi/4.pipe_tune2)) &
624	& **(1./(2.+pipe_tune3))
625	vn(i) = (ind(i,j) + d_ind(i,j))pipe_tune1MIN(dia(i),maxdia(j))**pipe_tune_exp_coeff
626
627	ENDIF
628	ENDDO
629	ELSE ! for grasses, cnd=1, so the above calculation cancels
630	vn(:) = ind(:,j) + d_ind(:,j)
631	ENDIF
632	ELSE ! static
633	DO i=1,npts
634	IF(tree(j).AND.(d_ind(i,j)+ind(i,j)).GT.min_stomate) THEN
635	IF(total_bm_c(i).LE.min_stomate) THEN
636
637	! new wood mass of PFT
638	woodmass_ind(i,j) = &
639	& (((biomass(i,j,isapabove) + biomass(i,j,isapbelow) &
640	& + biomass(i,j,iheartabove) + biomass(i,j,iheartbelow))) &
641	& + (bm_sapl(j,isapabove) + bm_sapl(j,isapbelow) &
642	& + bm_sapl(j,iheartabove) + bm_sapl(j,iheartbelow))*d_ind(i,j))/(ind(i,j)+d_ind(i,j))
643
644	ELSE ! new biomass is added to the labile pool, hence there is no change in CA associated with establishment
645
646	woodmass_ind(i,j) = &
647	& (biomass(i,j,isapabove) + biomass(i,j,isapbelow) &
648	& + biomass(i,j,iheartabove) + biomass(i,j,iheartbelow)) &
649	& /(ind(i,j) + d_ind(i,j))
650
651	ENDIF
652	ENDIF
653	ENDDO
654
655	vn(:) = un ! cannot change in static!, and veget_max implicit in d_ind
656
657	ENDIF
658	! total biomass of PFT added by establishment defined over veget_max ...
659	total_bm_sapl(:) = zero
660	DO k = 1, nparts
661	WHERE(d_ind(:,j).GT.min_stomate.AND.total_bm_c(:).GT.min_stomate.AND.vn(:).GT.min_stomate)
662
663	total_bm_sapl(:) = total_bm_sapl(:) + &
664	bm_sapl(j,k) * d_ind(:,j) / vn(:)
665	ENDWHERE
666	ENDDO
667
668	DO k = 1, nparts
669
670	bm_new(:) = zero
671
672	! first ever establishment, C flows
673	WHERE( d_ind(:,j).GT.min_stomate .AND. &
674	total_bm_c(:).LE.min_stomate .AND. &
675	vn(:).GT.min_stomate)
676	! WHERE ( many_new(:) )
677
678	!bm_new(:) = d_ind(:,j) * bm_sapl(j,k) / veget_max (:,j)
679	bm_new(:) = d_ind(:,j) * bm_sapl(j,k) / vn(:)
680
681	biomass(:,j,k) = biomass(:,j,k) + bm_new(:)
682
683	co2_to_bm(:,j) = co2_to_bm(:,j) + bm_new(:) / dt
684
685	ENDWHERE
686
687	! establishment into existing population, C flows
688	WHERE(d_ind(:,j).GT.min_stomate.AND.total_bm_c(:).GT.min_stomate)
689
690	bm_new(:) = total_bm_sapl(:) * biomass(:,j,k) / total_bm_c(:)
691
692	biomass(:,j,k) = biomass(:,j,k) + bm_new(:)
693	co2_to_bm(:,j) = co2_to_bm(:,j) + bm_new(:) / dt
694
695	ENDWHERE
696	ENDDO
697
698	! reset leaf ages. Should do a real calculation like in the npp routine,
699	! but this case is rare and not worth messing around.
700	! SZ 080806, added real calculation now, because otherwise leaf_age/leaf_frac
701	! are not initialised for the calculation of vmax, and hence no growth at all.
702	! logic follows that of stomate_npp.f90, just that it's been adjusted for the code here
703	!
704	! 4.5.2 Decrease leaf age in youngest class if new leaf biomass is higher than old one.
705	!
706
707	!!$ WHERE ( many_new(:) )
708	!!$ leaf_age(:,j,1) = zero
709	!!$ leaf_frac(:,j,1) = un
710	!!$ ENDWHERE
711	!!$
712	!!$ DO m = 2, nleafages
713	!!$
714	!!$ WHERE ( many_new(:) )
715	!!$ leaf_age(:,j,m) = zero
716	!!$ leaf_frac(:,j,m) = zero
717	!!$ ENDWHERE
718	!!$
719	!!$ ENDDO
720
721	WHERE ( d_ind(:,j) * bm_sapl(j,ileaf) .GT. min_stomate )
722
723	leaf_age(:,j,1) = leaf_age(:,j,1) * leaf_mass_young(:) / &
724	( leaf_mass_young(:) + d_ind(:,j) * bm_sapl(j,ileaf) )
725
726	ENDWHERE
727
728	leaf_mass_young(:) = leaf_mass_young(:) + d_ind(:,j) * bm_sapl(j,ileaf)
729
730	!
731	! new age class fractions (fraction in youngest class increases)
732	!
733
734	! youngest class: new mass in youngest class divided by total new mass
735
736	WHERE ( biomass(:,j,ileaf) .GT. min_stomate )
737
738	leaf_frac(:,j,1) = leaf_mass_young(:) / biomass(:,j,ileaf)
739
740	ENDWHERE
741
742	! other classes: old mass in leaf age class divided by new mass
743
744	DO m = 2, nleafages
745
746	WHERE ( biomass(:,j,ileaf) .GT. min_stomate )
747
748	leaf_frac(:,j,m) = leaf_frac(:,j,m) * &
749	( biomass(:,j,ileaf) + d_ind(:,j) * bm_sapl(j,ileaf) ) / biomass(:,j,ileaf)
750
751	ENDWHERE
752
753	ENDDO
754
755	!ENDIF ! establishment rate is large
756
757	WHERE ( d_ind(:,j) .GT. min_stomate )
758
759	! 4.5.3 age decreases
760
761	age(:,j) = age(:,j) * ind(:,j) / ( ind(:,j) + d_ind(:,j) )
762
763	! 4.5.4 new number of individuals
764
765	ind(:,j) = ind(:,j) + d_ind(:,j)
766
767	ENDWHERE
768
769	!
770	! 4.6 eventually convert excess sapwood to heartwood
771	!
772
773	!SZ to clarify with Gerhard Krinner: This is theoretically inconsistent because
774	! the allocation to sapwood and leaves do not follow the LPJ logic in stomate_alloc.f90
775	! hence imposing this here not only solves for the uneveness of age (mixing new and average individual)
776	! but also corrects for the discrepancy between SLAVE and LPJ logic of allocation, thus leads to excess heartwood
777	! and thus carbon accumulation!
778	! should be removed.
779
780	IF ( tree(j) ) THEN
781
782	!!$ sm2(:) = 0.0
783	!!$ WHERE ( d_ind(:,j) .GT. 0.0 )
784	!!$
785	!!$ ! ratio of above / total sap parts
786	!!$ sm_at(:) = biomass(:,j,isapabove) / &
787	!!$ ( biomass(:,j,isapabove) + biomass(:,j,isapbelow) )
788	!!$
789	!!$ ! woodmass of an individual
790	!!$
791	!!$ woodmass(:) = &
792	!!$ ( biomass(:,j,isapabove) + biomass(:,j,isapbelow) + &
793	!!$ biomass(:,j,iheartabove) + biomass(:,j,iheartbelow) ) / ind(:,j)
794	!!$
795	!!$ ! crown area (m**2) depends on stem diameter (pipe model)
796	!!$ dia(:) = ( woodmass(:) / ( pipe_density * pi/4. * pipe_tune2 ) ) &
797	!!$ ** ( 1. / ( 2. + pipe_tune3 ) )
798	!!$
799	!!$ b1(:) = pipe_k1 / ( sla(j) * pipe_densitypipe_tune2 dia(:)*pipe_tune3 ) &
800	!!$ ind(:,j)
801	!!$ sm2(:) = lm_lastyearmax(:,j) / b1(:)
802	!!$
803	!!$ ENDWHERE
804
805	sm2(:) = biomass(:,j,isapabove) + biomass(:,j,isapbelow)
806
807	WHERE ( ( d_ind(:,j) .GT. min_stomate ) .AND. &
808	( biomass(:,j,isapabove) + biomass(:,j,isapbelow) ) .GT. sm2(:) )
809
810	sm_at(:) = biomass(:,j,isapabove) / &
811	( biomass(:,j,isapabove) + biomass(:,j,isapbelow) )
812
813	biomass(:,j,iheartabove) = biomass(:,j,iheartabove) + &
814	( biomass(:,j,isapabove) - sm2(:) * sm_at(:) )
815	biomass(:,j,isapabove) = sm2(:) * sm_at(:)
816
817	biomass(:,j,iheartbelow) = biomass(:,j,iheartbelow) + &
818	( biomass(:,j,isapbelow) - sm2(:) * (un - sm_at) )
819	biomass(:,j,isapbelow) = sm2(:) * (un - sm_at(:))
820
821	ENDWHERE
822
823	ENDIF ! tree
824
825	ENDIF ! natural
826
827	ENDDO ! loop over pfts
828
829	!
830	! 5 history
831	!
832
833	d_ind = d_ind / dt
834
835	CALL histwrite (hist_id_stomate, 'IND_ESTAB', itime, d_ind, npts*nvm, horipft_index)
836	CALL histwrite (hist_id_stomate, 'ESTABTREE', itime, estab_rate_max_tree, npts, hori_index)
837	CALL histwrite (hist_id_stomate, 'ESTABGRASS', itime, estab_rate_max_grass, npts, hori_index)
838
839	IF (bavard.GE.4) WRITE(numout,*) 'Leaving establish'
840
841	END SUBROUTINE establish
842
843	END MODULE lpj_establish

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

Download in other formats: