Context Navigation

source: tags/ORCHIDEE/src_stomate/lpj_establish.f90 @ 6

Last change on this file since 6 was 6, checked in by orchidee, 14 years ago
import first tag equivalent to CVS orchidee_1_9_5 + OOL_1_9_5
File size: 18.1 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_constants
14	USE constantes_veg
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, &
36	leaf_age, leaf_frac, &
37	ind, biomass, age, everywhere, co2_to_bm,veget_max)
38
39	!
40	! 0 declarations
41	!
42
43	! 0.1 input
44
45	! Domain size
46	INTEGER(i_std), INTENT(in) :: npts
47	! Time step of vegetation dynamics (days)
48	REAL(r_std), INTENT(in) :: dt
49	! Is pft there
50	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: PFTpresent
51	! Winter sufficiently cold
52	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: regenerate
53	! 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)
54	INTEGER(i_std), DIMENSION(npts,8), INTENT(in) :: neighbours
55	! resolution at each grid point in m (1=E-W, 2=N-S)
56	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: resolution
57	! in order for this PFT to be introduced, does it have to be present in an
58	! adjacent grid box?
59	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: need_adjacent
60	! time constant of probability of a leaf to be eaten by a herbivore (days)
61	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: herbivores
62	! annual precipitation (mm/year)
63	REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_annual
64	! growing degree days (C)
65	REAL(r_std), DIMENSION(npts), INTENT(in) :: gdd0
66	! last year's maximum leaf mass, for each PFT (gC/(m**2 of ground))
67	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lm_lastyearmax
68	! crown area of individuals (m**2)
69	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: cn_ind
70	! leaf area index OF AN INDIVIDUAL PLANT
71	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lai
72	! space availability for trees
73	REAL(r_std), DIMENSION(npts), INTENT(in) :: avail_tree
74	! space availability for grasses
75	REAL(r_std), DIMENSION(npts), INTENT(in) :: avail_grass
76	! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground
77	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_max
78
79	! 0.2 modified fields
80
81	! leaf age (days)
82	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_age
83	! fraction of leaves in leaf age class
84	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac
85	! Number of individuals / m2
86	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind
87	! biomass (gC/(m**2 of ground))
88	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout) :: biomass
89	! mean age (years)
90	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age
91	! is the PFT everywhere in the grid box or very localized (after its introduction)
92	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: everywhere
93	! biomass uptaken (gC/(m**2 of total ground)/day)
94	!NV passage 2D
95	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: co2_to_bm
96
97	! 0.3 local
98
99	! time during which a sapling can be entirely eaten by herbivores (d)
100	REAL(r_std) :: tau_eatup
101	! new fpc ( foliage protected cover: fractional coverage )
102	REAL(r_std), DIMENSION(npts,nvm) :: fpc_nat
103	! maximum tree establishment rate, based on climate only
104	REAL(r_std), DIMENSION(npts) :: estab_rate_max_climate_tree
105	! maximum grass establishment rate, based on climate only
106	REAL(r_std), DIMENSION(npts) :: estab_rate_max_climate_grass
107	! maximum tree establishment rate, based on climate and fpc
108	REAL(r_std), DIMENSION(npts) :: estab_rate_max_tree
109	! maximum grass establishment rate, based on climate and fpc
110	REAL(r_std), DIMENSION(npts) :: estab_rate_max_grass
111	! total natural fpc
112	REAL(r_std), DIMENSION(npts) :: sumfpc
113	! total woody fpc
114	REAL(r_std), DIMENSION(npts) :: sumfpc_wood
115	! for trees, measures the total concurrence for available space
116	REAL(r_std), DIMENSION(npts) :: spacefight_tree
117	! for grasses, measures the total concurrence for available space
118	REAL(r_std), DIMENSION(npts) :: spacefight_grass
119	! change in number of individuals /m2 per time step (per day in history file)
120	REAL(r_std), DIMENSION(npts,nvm) :: d_ind
121	! biomass increase (gC/(m**2 of ground))
122	REAL(r_std), DIMENSION(npts) :: bm_new
123	! stem diameter (m)
124	REAL(r_std), DIMENSION(npts) :: dia
125	! temporary variable
126	REAL(r_std), DIMENSION(npts) :: b1
127	! new sap mass (gC/(m**2 of ground))
128	REAL(r_std), DIMENSION(npts) :: sm2
129	! woodmass of an individual
130	REAL(r_std), DIMENSION(npts) :: woodmass
131	! ratio of hw(above) to total hw, sm(above) to total sm
132	REAL(r_std), DIMENSION(npts) :: sm_at
133	! reduction factor for establishment if many trees or grasses are present
134	REAL(r_std), DIMENSION(npts) :: factor
135	! from how many sides is the grid box invaded
136	INTEGER(i_std) :: nfrontx
137	INTEGER(i_std) :: nfronty
138	! daily establishment rate is large compared to present number of individuals
139	LOGICAL, DIMENSION(npts) :: many_new
140	! indices
141	INTEGER(i_std) :: i,j,k,m
142
143	! =========================================================================
144
145	IF (bavard.GE.3) WRITE(numout,*) 'Entering establish'
146
147	!
148	! 1 messages and initialization
149	!
150	tau_eatup = one_year/2.
151
152	IF ( firstcall ) THEN
153
154	WRITE(numout,*) 'establish:'
155
156	WRITE(numout,*) ' > time during which a sapling can be entirely eaten by herbivores (d): ', &
157	tau_eatup
158
159	firstcall = .FALSE.
160
161	ENDIF
162
163	!
164	! 2 recalculate fpc
165	!
166
167	!
168	! 2.1 Only natural part of the grid cell
169	!
170
171	DO j = 2,nvm
172
173	IF ( natural(j) ) THEN
174	DO i = 1, npts
175	IF (lai(i,j) == val_exp) THEN
176	fpc_nat(i,j) = cn_ind(i,j) * ind(i,j)
177	ELSE
178	fpc_nat(i,j) = cn_ind(i,j) * ind(i,j) * ( 1. - exp( -lai(i,j) * ext_coeff(j) ) )
179	ENDIF
180	ENDDO
181	ELSE
182
183	fpc_nat(:,j) = 0.0
184
185	ENDIF
186
187	ENDDO
188
189	!
190	! 2.2 total natural fpc on grid
191	!
192
193	sumfpc(:) = SUM( fpc_nat(:,:), DIM=2 )
194
195	!
196	! 2.3 total woody fpc on grid and number of regenerative tree pfts
197	!
198
199	sumfpc_wood(:) = 0.0
200	spacefight_tree(:) = 0.0
201
202	DO j = 2,nvm
203
204	IF ( tree(j) .AND. natural(j) ) THEN
205
206	! total woody fpc
207
208	WHERE ( PFTpresent(:,j) )
209	sumfpc_wood(:) = sumfpc_wood(:) + fpc_nat(:,j)
210	ENDWHERE
211
212	! how many trees are competing? Count a PFT fully only if it is present
213	! on the whole grid box.
214
215	WHERE ( PFTpresent(:,j) .AND. ( regenerate(:,j) .GT. regenerate_crit ) )
216	spacefight_tree(:) = spacefight_tree(:) + everywhere(:,j)
217	ENDWHERE
218
219	ENDIF
220
221	ENDDO
222
223	!
224	! 2.4 number of natural grasses
225	!
226
227	spacefight_grass(:) = 0.0
228
229	DO j = 2,nvm
230
231	IF ( .NOT. tree(j) .AND. natural(j) ) THEN
232
233	! how many grasses are competing? Count a PFT fully only if it is present
234	! on the whole grid box.
235
236	WHERE ( PFTpresent(:,j) )
237	spacefight_grass(:) = spacefight_grass(:) + everywhere(:,j)
238	ENDWHERE
239
240	ENDIF
241
242	ENDDO
243
244	!
245	! 3 establishment rate
246	!
247
248	!
249	! 3.1 maximum establishment rate, based on climate only
250	!
251
252	WHERE ( ( precip_annual(:) .GE. precip_crit ) .AND. ( gdd0(:) .GE. gdd_crit ) )
253
254	estab_rate_max_climate_tree(:) = estab_max_tree
255	estab_rate_max_climate_grass(:) = estab_max_grass
256
257	ELSEWHERE
258
259	estab_rate_max_climate_tree(:) = 0.0
260	estab_rate_max_climate_grass(:) = 0.0
261
262	ENDWHERE
263
264	!
265	! 3.2 reduce maximum tree establishment rate if many trees present.
266	! In the original DGVM, this is done using a step function which yields a
267	! reduction by factor 4 if sumfpc_wood(i) .GT. fpc_crit - 0.05.
268	! This can lead to small oscillations (without consequences however).
269	! Here, a steady linear transition is used between fpc_crit-0.075 and
270	! fpc_crit-0.025.
271	!
272
273	factor(:) = 1. - 15. * ( sumfpc_wood(:) - (fpc_crit-.075) )
274	factor(:) = MAX( 0.25_r_std, MIN( 1._r_std, factor(:) ) )
275
276	estab_rate_max_tree(:) = estab_rate_max_climate_tree(:) * factor(:)
277
278	!
279	! 3.3 Modulate grass establishment rate.
280	! If canopy is not closed (fpc < fpc_crit-0.05), normal establishment.
281	! If canopy is closed, establishment is reduced by a factor 4.
282	! Factor is linear between these two bounds.
283	! This is different from the original DGVM where a step function is
284	! used at fpc_crit-0.05 (This can lead to small oscillations,
285	! without consequences however).
286	!
287
288	factor(:) = 1. - 15. * ( sumfpc(:) - (fpc_crit-.05) )
289	factor(:) = MAX( 0.25_r_std, MIN( 1._r_std, factor(:) ) )
290
291	estab_rate_max_grass(:) = estab_rate_max_climate_grass(:) * factor(:)
292
293	!
294	! 4 do establishment for natural PFTs
295	!
296
297	d_ind(:,:) = 0.0
298
299	DO j = 2,nvm
300
301	! only for natural PFTs
302
303	IF ( natural(j) ) THEN
304
305	!
306	! 4.1 PFT expansion across the grid box. Not to be confused with areal
307	! coverage.
308	!
309
310	IF ( treat_expansion ) THEN
311
312	! only treat plants that are regenerative and present and still can expand
313
314	DO i = 1, npts
315
316	IF ( PFTpresent(i,j) .AND. &
317	( everywhere(i,j) .LT. 1. ) .AND. &
318	( regenerate(i,j) .GT. regenerate_crit ) ) THEN
319
320	! from how many sides is the grid box invaded (separate x and y directions
321	! because resolution may be strongly anisotropic)
322	!
323	! For the moment we only look into 4 direction but that can be extanded (JP)
324	!
325	nfrontx = 0
326	IF ( neighbours(i,3) .GT. 0 ) THEN
327	IF ( everywhere(neighbours(i,3),j) .GT. 1.-min_stomate ) nfrontx = nfrontx+1
328	ENDIF
329	IF ( neighbours(i,7) .GT. 0 ) THEN
330	IF ( everywhere(neighbours(i,7),j) .GT. 1.-min_stomate ) nfrontx = nfrontx+1
331	ENDIF
332
333	nfronty = 0
334	IF ( neighbours(i,1) .GT. 0 ) THEN
335	IF ( everywhere(neighbours(i,1),j) .GT. 1.-min_stomate ) nfronty = nfronty+1
336	ENDIF
337	IF ( neighbours(i,5) .GT. 0 ) THEN
338	IF ( everywhere(neighbours(i,5),j) .GT. 1.-min_stomate ) nfronty = nfronty+1
339	ENDIF
340
341	everywhere(i,j) = &
342	everywhere(i,j) + migrate(j) * dt/one_year * &
343	( nfrontx / resolution(i,1) + nfronty / resolution(i,2) )
344
345	IF ( .NOT. need_adjacent(i,j) ) THEN
346
347	! in that case, we also assume that the PFT expands from places within
348	! the grid box (e.g., oasis).
349
350	everywhere(i,j) = &
351	everywhere(i,j) + migrate(j) * dt/one_year * &
352	2. * SQRT( pieverywhere(i,j)/(resolution(i,1)resolution(i,2)) )
353
354	ENDIF
355
356	everywhere(i,j) = MIN( everywhere(i,j), 1._r_std )
357
358	ENDIF
359
360	ENDDO
361
362	ENDIF ! treat expansion?
363
364	!
365	! 4.2 establishment rate
366	! - Is lower if the PFT is only present in a small part of the grid box
367	! (after its introduction), therefore multiplied by "everywhere".
368	! - Is divided by the number of PFTs that compete ("spacefight").
369	! - Is modulated by space availability (avail_tree, avail_grass).
370	!
371
372	IF ( tree(j) ) THEN
373
374	! 4.2.1 present and regenerative trees
375
376	WHERE ( PFTpresent(:,j) .AND. ( regenerate(:,j) .GT. regenerate_crit ) )
377
378
379	d_ind(:,j) = estab_rate_max_tree(:)everywhere(:,j)/spacefight_tree(:) &
380	avail_tree(:) * dt/one_year
381
382	ENDWHERE
383
384	ELSE
385
386	! 4.2.2 present and regenerative grasses
387
388	WHERE ( PFTpresent(:,j) .AND. ( regenerate(:,j) .GT. regenerate_crit ) )
389
390	d_ind(:,j) = estab_rate_max_grass(:)everywhere(:,j)/spacefight_grass(:) &
391	avail_grass(:) * dt/one_year
392
393	ENDWHERE
394
395	ENDIF ! tree/grass
396
397	!
398	! 4.3 herbivores reduce establishment rate
399	! We suppose that saplings are vulnerable during a given time after establishment.
400	! This is taken into account by preventively reducing the establishment rate.
401	!
402
403	IF ( ok_herbivores ) THEN
404
405	d_ind(:,j) = d_ind(:,j) * EXP( - tau_eatup/herbivores(:,j) )
406
407	ENDIF
408
409	!
410	! 4.4 be sure that ind*cn_ind does not exceed 1
411	!
412
413	WHERE ( ( d_ind(:,j) .GT. 0.0 ) .AND. &
414	( (ind(:,j)+d_ind(:,j))*cn_ind(:,j) .GT. 1. ) )
415
416	d_ind(:,j) = MAX( 1._r_std / cn_ind(:,j) - ind(:,j), 0._r_std )
417
418	ENDWHERE
419
420	!
421	! 4.5 new properties where there is establishment ( d_ind > 0 )
422	!
423
424	! 4.5.1 biomass.
425	! Add biomass only if d_ind, over one year, is of the order of ind.
426	! (If we don't do this, the biomass density can become very low).
427	! In that case, take biomass from the atmosphere.
428
429	! compare establishment rate and present number of inidivuals
430	many_new(:) = ( d_ind(:,j) .GT. 0.1 * ind(:,j) )
431
432	! gives a better vectorization of the VPP
433
434	IF ( ANY( many_new(:) ) ) THEN
435
436	DO k = 1, nparts
437
438	WHERE ( many_new(:) )
439
440	bm_new(:) = d_ind(:,j) * bm_sapl(j,k) / veget_max (:,j)
441
442	biomass(:,j,k) = biomass(:,j,k) + bm_new(:)
443
444	!NV passage 2D
445	co2_to_bm(:,j) = co2_to_bm(:,j) + bm_new(:) / dt
446
447	ENDWHERE
448
449	ENDDO
450
451	! reset leaf ages. Should do a real calculation like in the npp routine,
452	! but this case is rare and not worth messing around.
453
454	WHERE ( many_new(:) )
455	leaf_age(:,j,1) = 0.0
456	leaf_frac(:,j,1) = 1.0
457	ENDWHERE
458
459	DO m = 2, nleafages
460
461	WHERE ( many_new(:) )
462	leaf_age(:,j,m) = 0.0
463	leaf_frac(:,j,m) = 0.0
464	ENDWHERE
465
466	ENDDO
467
468	ENDIF ! establishment rate is large
469
470	WHERE ( d_ind(:,j) .GT. 0.0 )
471
472	! 4.5.2 age decreases
473
474	age(:,j) = age(:,j) * ind(:,j) / ( ind(:,j) + d_ind(:,j) )
475
476	! 4.5.3 new number of individuals
477
478	ind(:,j) = ind(:,j) + d_ind(:,j)
479
480	ENDWHERE
481
482	!
483	! 4.6 eventually convert excess sapwood to heartwood
484	!
485
486	IF ( tree(j) ) THEN
487
488	sm2(:) = 0.0
489
490	WHERE ( d_ind(:,j) .GT. 0.0 )
491
492	! ratio of above / total sap parts
493	sm_at(:) = biomass(:,j,isapabove) / &
494	( biomass(:,j,isapabove) + biomass(:,j,isapbelow) )
495
496	! woodmass of an individual
497
498	woodmass(:) = &
499	( biomass(:,j,isapabove) + biomass(:,j,isapbelow) + &
500	biomass(:,j,iheartabove) + biomass(:,j,iheartbelow) ) / ind(:,j)
501
502	! crown area (m**2) depends on stem diameter (pipe model)
503	dia(:) = ( woodmass(:) / ( pipe_density * pi/4. * pipe_tune2 ) ) &
504	** ( 1. / ( 2. + pipe_tune3 ) )
505
506	b1(:) = pipe_k1 / ( sla(j) * pipe_densitypipe_tune2 dia(:)*pipe_tune3 ) &
507	ind(:,j)
508	sm2(:) = lm_lastyearmax(:,j) / b1(:)
509
510	ENDWHERE
511
512	WHERE ( ( d_ind(:,j) .GT. 0.0 ) .AND. &
513	( biomass(:,j,isapabove) + biomass(:,j,isapbelow) ) .GT. sm2(:) )
514
515	biomass(:,j,iheartabove) = biomass(:,j,iheartabove) + &
516	( biomass(:,j,isapabove) - sm2(:) * sm_at(:) )
517	biomass(:,j,isapabove) = sm2(:) * sm_at(:)
518
519	biomass(:,j,iheartbelow) = biomass(:,j,iheartbelow) + &
520	( biomass(:,j,isapbelow) - sm2(:) * (1. - sm_at) )
521	biomass(:,j,isapbelow) = sm2(:) * (1. - sm_at(:))
522
523	ENDWHERE
524
525	ENDIF ! tree
526
527	ENDIF ! natural
528
529	ENDDO ! loop over pfts
530
531	!
532	! 5 history
533	!
534
535	d_ind = d_ind / dt
536
537	CALL histwrite (hist_id_stomate, 'IND_ESTAB', itime, d_ind, npts*nvm, horipft_index)
538
539	IF (bavard.GE.4) WRITE(numout,*) 'Leaving establish'
540
541	END SUBROUTINE establish
542
543	END MODULE lpj_establish

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

Download in other formats: