Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE/src_stomate/lpj_establish.f90 @ 257

Last change on this file since 257 was 257, checked in by didier.solyga, 13 years ago
Externalized version merged with the trunk
File size: 31.9 KB

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

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

Download in other formats: