Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE/src_stomate/lpj_light.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: 23.0 KB

Line
1	! Light competition
2	!
3	! If canopy is almost closed (fpc > fpc_crit), then trees outcompete grasses.
4	! fpc_crit is normally fpc_crit.
5	! Here, fpc ("foilage protected cover") also takes into account the minimum fraction
6	! of space covered by trees through branches etc. This is done to prevent strong relative
7	! changes of fpc from one day to another for deciduous trees at the beginning of their
8	! growing season, which would yield to strong cutbacks (see 3.2.1.1.2)
9	! No competition between woody pfts (height of individuals is not considered) !
10	! Exception: when one woody pft is overwhelming (i.e. fpc > fpc_crit). In that
11	! case, eliminate all other woody pfts and reduce dominant pft to fpc_crit.
12	! Age of individuals is not considered. In reality, light competition would more
13	! easily kill young individuals, thus increasing the mean age of the stand.
14	! Exclude agricultural pfts from competition
15	!
16	! SZ: added light competition for the static case if the mortality is not
17	! assumed to be constant.
18	! other modifs:
19	! -1 FPC is now always calculated from lm_lastyearmax*sla, since the aim of this DGVM is
20	! to represent community ecology effects; seasonal variations in establishment related to phenology
21	! may be relevant, but beyond the scope of a 1st generation DGVM
22	! -2 problem, if agriculture is present, fpc can never reach 1.0 since natural veget_max < 1.0. To
23	! correct for this, ind must be recalculated to correspond to the natural density...
24	! since ind is 1/m2 grid cell, this can be achived by dividing ind by the agricultural fraction
25
26	!
27	! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/lpj_light.f90,v 1.8 2009/01/06 15:01:25 ssipsl Exp $
28	! IPSL (2006)
29	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
30	!
31	MODULE lpj_light
32
33	! modules used:
34
35	USE ioipsl
36	USE constantes
37	USE stomate_data
38
39	IMPLICIT NONE
40
41	! private & public routines
42
43	PRIVATE
44	PUBLIC light, light_clear
45
46	! first call
47	LOGICAL, SAVE :: firstcall = .TRUE.
48
49	CONTAINS
50
51	SUBROUTINE light_clear
52	firstcall=.TRUE.
53	END SUBROUTINE light_clear
54
55	SUBROUTINE light (npts, dt, &
56	veget_max, fpc_max, PFTpresent, cn_ind, lai, maxfpc_lastyear, &
57	lm_lastyearmax, ind, biomass, veget_lastlight, bm_to_litter, mortality)
58
59	!
60	! 0 declarations
61	!
62
63	! 0.1 input
64
65	! Domain size
66	INTEGER(i_std), INTENT(in) :: npts
67	! Time step (days)
68	REAL(r_std), INTENT(in) :: dt
69	! Is pft there
70	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: PFTpresent
71	! crown area of individuals (m**2)
72	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: cn_ind
73	! leaf area index OF AN INDIVIDUAL PLANT
74	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lai
75	! last year's maximum fpc for each natural PFT, on ground
76	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: maxfpc_lastyear
77	! last year's maximum leafmass for each natural PFT, on ground
78	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lm_lastyearmax
79	! last year's maximum fpc for each natural PFT, on ground
80	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max
81	! last year's maximum fpc for each natural PFT, on ground
82	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: fpc_max
83
84	! 0.2 modified fields
85
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	! Vegetation cover after last light competition
91	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: veget_lastlight
92	! biomass taken away (gC/m**2)
93	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout) :: bm_to_litter
94	! fraction of individuals that died this time step
95	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: mortality
96
97	! 0.3 local
98
99	! index
100	INTEGER(i_std) :: i,j,k,m
101	! total natural fpc
102	REAL(r_std), DIMENSION(npts) :: sumfpc
103	! fraction of natural vegetation at grid cell level
104	REAL(r_std), DIMENSION(npts) :: fracnat
105	! total natural woody fpc
106	REAL(r_std) :: sumfpc_wood
107	! change in total woody fpc
108	REAL(r_std) :: sumdelta_fpc_wood
109	! maximum wood fpc
110	REAL(r_std) :: maxfpc_wood
111	! which woody pft is maximum
112	INTEGER(i_std) :: optpft_wood
113	! total natural grass fpc
114	REAL(r_std) :: sumfpc_grass
115	! this year's foliage protected cover on natural part of the grid cell
116	REAL(r_std), DIMENSION(npts,nvm) :: fpc_nat
117	! fpc change within last year
118	REAL(r_std), DIMENSION(nvm) :: deltafpc
119	! Relative change of number of individuals for trees
120	REAL(r_std) :: reduct
121	! Fraction of plants that survive
122	REAL(r_std), DIMENSION(nvm) :: survive
123	! FPC for static mode
124	REAL(r_std), DIMENSION(npts) :: fpc_real
125	! FPC mortality for static mode
126	REAL(r_std), DIMENSION(npts) :: lai_ind
127	! number of grass PFTs present in the grid box
128	! INTEGER(i_std) :: num_grass
129	! New total grass fpc
130	REAL(r_std) :: sumfpc_grass2
131	! fraction of plants that dies each day (1/day)
132	REAL(r_std), DIMENSION(npts,nvm) :: light_death
133	! Relative change of number of individuals for trees
134	REAL(r_std) :: fpc_dec
135
136	! =========================================================================
137
138	IF (bavard.GE.3) WRITE(numout,*) 'Entering light'
139
140	!
141	! 1 first call
142	!
143
144	IF ( firstcall ) THEN
145
146	WRITE(numout,*) 'light:'
147
148	WRITE(numout,*) ' > Maximum total number of grass individuals in'
149	WRITE(numout,*) ' a closed canopy: ', grass_mercy
150
151	WRITE(numout,*) ' > Minimum fraction of trees that survive even in'
152	WRITE(numout,*) ' a closed canopy: ', tree_mercy
153
154	WRITE(numout,*) ' > For trees, minimum fraction of crown area covered'
155	WRITE(numout,*) ' (due to its branches etc.)', min_cover
156
157	WRITE(numout,*) ' > for diagnosis of fpc increase, compare today''s fpc'
158	IF ( annual_increase ) THEN
159	WRITE(numout,*) ' to last year''s maximum.'
160	ELSE
161	WRITE(numout,*) ' to fpc of the last time step.'
162	ENDIF
163
164	firstcall = .FALSE.
165
166	ENDIF
167
168	IF (control%ok_dgvm) THEN
169	!
170	! 2 fpc characteristics
171	!
172
173	! 2.0 Only natural part of the grid cell:
174	! calculate fraction of natural and agricultural (1-fracnat) surface
175
176	fracnat(:) = 1.
177	DO j = 2,nvm
178	IF ( .NOT. natural(j) ) THEN
179	fracnat(:) = fracnat(:) - veget_max(:,j)
180	ENDIF
181	ENDDO
182	!
183	! 2.1 calculate fpc on natural part of grid cell.
184	!
185	fpc_nat(:,:)=zero
186	fpc_nat(:,ibare_sechiba)=un
187
188	DO j = 2, nvm
189
190	IF ( natural(j) ) THEN
191
192	! 2.1.1 natural PFTs
193
194	IF ( tree(j) ) THEN
195
196	! 2.1.1.1 trees: minimum cover due to stems, branches etc.
197
198	! DO i = 1, npts
199	! IF (lai(i,j) == val_exp) THEN
200	! fpc_nat(i,j) = cn_ind(i,j) * ind(i,j)
201	! ELSE
202	! fpc_nat(i,j) = cn_ind(i,j) * ind(i,j) * &
203	! MAX( ( 1._r_std - exp( -lai(i,j) * ext_coeff(j) ) ), min_cover )
204	! ENDIF
205	! ENDDO
206
207	!NV : modif from SZ version : fpc is based on veget_max, not veget.
208	WHERE(fracnat(:).GE.min_stomate)
209	! WHERE(LAI(:,j) == val_exp)
210	! fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) / fracnat(:)
211	! ELSEWHERE
212	! fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) / fracnat(:) * &
213	! MAX( ( 1._r_std - exp( - lm_lastyearmax(:,j) * sla(j) * ext_coeff(j) ) ), min_cover )
214	! ENDWHERE
215	fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) / fracnat(:)
216	ENDWHERE
217
218	ELSE
219
220	!NV : modif from SZ version : fpc is based on veget_max, not veget.
221	WHERE(fracnat(:).GE.min_stomate)
222	! WHERE(LAI(:,j) == val_exp)
223	! fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) / fracnat(:)
224	! ELSEWHERE
225	! fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) / fracnat(:) * &
226	! ( 1._r_std - exp( - lm_lastyearmax(:,j) * sla(j) * ext_coeff(j) ) )
227	! ENDWHERE
228	fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) / fracnat(:)
229	ENDWHERE
230
231	!!$ ! 2.1.1.2 bare ground
232	!!$ IF (j == ibare_sechiba) THEN
233	!!$ fpc_nat(:,j) = cn_ind(:,j) * ind(:,j)
234	!!$
235	!!$ ! 2.1.1.3 grasses
236	!!$ ELSE
237	!!$ DO i = 1, npts
238	!!$ IF (lai(i,j) == val_exp) THEN
239	!!$ fpc_nat(i,j) = cn_ind(i,j) * ind(i,j)
240	!!$ ELSE
241	!!$ fpc_nat(i,j) = cn_ind(i,j) * ind(i,j) * &
242	!!$ ( 1._r_std - exp( -lai(i,j) * ext_coeff(j) ) )
243	!!$ ENDIF
244	!!$ ENDDO
245	!!$ ENDIF
246
247	ENDIF ! tree/grass
248
249	ELSE
250
251	! 2.1.2 agricultural PFTs: not present on natural part
252
253	fpc_nat(:,j) = zero
254
255	ENDIF ! natural/agricultural
256
257	ENDDO
258
259	!
260	! 2.2 sum natural fpc for every grid point
261	!
262
263	sumfpc(:) = zero
264	DO j = 2,nvm
265	!SZ bug correction MERGE: need to subtract agricultural area!
266	sumfpc(:) = sumfpc(:) + fpc_nat(:,j)
267	ENDDO
268
269	!
270	! 3 Light competition
271	!
272
273	light_death(:,:) = zero
274
275	DO i = 1, npts ! SZ why this loop and not a vector statement ?
276
277	! Only if vegetation cover is dense
278
279	IF ( sumfpc(i) .GT. fpc_crit ) THEN
280
281	! fpc change for each pft
282	! There are two possibilities: either we compare today's fpc with the fpc after the last
283	! time step, or we compare it to last year's maximum fpc of that PFT. In the first case,
284	! the fpc increase will be strong for seasonal PFTs at the beginning of the growing season.
285	! As for trees, the cutback is proportional to this increase, this means that seasonal trees
286	! will be disadvantaged compared to evergreen trees. In the original LPJ model, with its
287	! annual time step, the second method was used (this corresponds to annual_increase=.TRUE.)
288
289	IF ( annual_increase ) THEN
290	deltafpc(:) = MAX( (fpc_nat(i,:)-maxfpc_lastyear(i,:)), zero )
291	ELSE
292	deltafpc(:) = MAX( (fpc_nat(i,:)-veget_lastlight(i,:)), zero )
293	ENDIF
294
295	! default: survive
296
297	survive(:) = 1.0
298
299	!
300	! 3.1 determine some characteristics of the fpc distribution
301	!
302
303	sumfpc_wood = zero
304	sumdelta_fpc_wood = zero
305	maxfpc_wood = zero
306	optpft_wood = 0
307	sumfpc_grass = zero
308	! num_grass = 0
309
310	DO j = 2,nvm
311
312	! only natural pfts
313
314	IF ( natural(j) ) THEN
315
316	IF ( tree(j) ) THEN
317
318	! trees
319
320	! total woody fpc
321
322	sumfpc_wood = sumfpc_wood + fpc_nat(i,j)
323
324	! how much did the woody fpc increase
325
326	sumdelta_fpc_wood = sumdelta_fpc_wood + deltafpc(j)
327
328	! which woody pft is preponderant
329
330	IF ( fpc_nat(i,j) .GT. maxfpc_wood ) THEN
331
332	optpft_wood = j
333
334	maxfpc_wood = fpc_nat(i,j)
335
336	ENDIF
337
338	ELSE
339
340	! grasses
341
342	! total (natural) grass fpc
343
344	sumfpc_grass = sumfpc_grass + fpc_nat(i,j)
345
346	! number of grass PFTs present in the grid box
347
348	! IF ( PFTpresent(i,j) ) THEN
349	! num_grass = num_grass + 1
350	! ENDIF
351
352	ENDIF ! tree or grass
353
354	ENDIF ! natural
355
356	ENDDO ! loop over pfts
357
358	!
359	! 3.2 light competition: assume wood outcompetes grass
360	!
361	!SZ
362	!!$ IF (sumfpc_wood .GE. fpc_crit ) THEN
363
364	!
365	! 3.2.1 all allowed natural space is covered by wood:
366	! cut back trees to fpc_crit.
367	! Original DGVM: kill grasses. Modified: we let a very
368	! small fraction of grasses survive.
369	!
370
371	DO j = 2,nvm
372
373	! only present and natural pfts compete
374
375	IF ( PFTpresent(i,j) .AND. natural(j) ) THEN
376
377	IF ( tree(j) ) THEN
378
379	!
380	! 3.2.1.1 tree
381	!
382
383	! no single woody pft is overwhelming
384	! (original DGVM: tree_mercy = 0.0 )
385	! The reduction rate is proportional to the ratio deltafpc/fpc.
386
387	IF (sumfpc_wood .GE. fpc_crit .AND. fpc_nat(i,j) .GT. min_stomate .AND. &
388	sumdelta_fpc_wood .GT. min_stomate) THEN
389
390	! reduct = MIN( ( ( deltafpc(j)/sumdelta_fpc_wood * &
391	! (sumfpc_wood-fpc_crit) ) / fpc_nat(i,j) ), &
392	! ( 1._r_std - tree_mercy ) )
393	reduct = un - MIN((fpc_nat(i,j)-(sumfpc_wood-fpc_crit) &
394	* deltafpc(j)/sumdelta_fpc_wood)/fpc_nat(i,j), un )
395
396	ELSE
397
398	! tree fpc didn't icrease or it started from nothing
399
400	reduct = zero
401
402	ENDIF
403
404	survive(j) = un - reduct
405
406	ELSE
407
408	!
409	! 3.2.1.2 grass: let a very small fraction survive (the sum of all
410	! grass individuals may make up a maximum cover of
411	! grass_mercy [for lai -> infinity]).
412	! In the original DGVM, grasses were killed in that case,
413	! corresponding to grass_mercy = 0.
414	!
415
416	! survive(j) = ( grass_mercy / REAL( num_grass,r_std ) ) / ind(i,j)
417
418	! survive(j) = MIN( 1._r_std, survive(j)
419
420	IF(sumfpc_grass .GE. 1.0-MIN(fpc_crit,sumfpc_wood).AND. &
421	sumfpc_grass.GE.min_stomate) THEN
422
423	fpc_dec=(sumfpc_grass-1.+MIN(fpc_crit,sumfpc_wood))*fpc_nat(i,j)/sumfpc_grass
424
425	reduct=fpc_dec
426	ELSE
427	reduct = zero
428	ENDIF
429	survive(j) = ( un - reduct )
430
431	ENDIF ! tree or grass
432
433	ENDIF ! pft there and natural
434
435	ENDDO ! loop over pfts
436
437	!SZ
438	!!$ ELSE
439	!!$
440	!!$ !
441	!!$ ! 3.2.2 not too much wood so that grasses can subsist
442	!!$ !
443	!!$
444	!!$ ! new total grass fpc
445	!!$ sumfpc_grass2 = fpc_crit - sumfpc_wood
446	!!$
447	!!$ DO j = 2,nvm
448	!!$
449	!!$ ! only present and natural PFTs compete
450	!!$
451	!!$ IF ( PFTpresent(i,j) .AND. natural(j) ) THEN
452	!!$
453	!!$ IF ( tree(j) ) THEN
454	!!$
455	!!$ ! no change for trees
456	!!$
457	!!$ survive(j) = 1.0
458	!!$
459	!!$ ELSE
460	!!$
461	!!$ ! grass: fractional loss is the same for all grasses
462	!!$
463	!!$ IF ( sumfpc_grass .GT. min_stomate ) THEN
464	!!$ survive(j) = sumfpc_grass2 / sumfpc_grass
465	!!$ ELSE
466	!!$ survive(j)= zero
467	!!$ ENDIF
468	!!$
469	!!$ ENDIF
470	!!$
471	!!$ ENDIF ! pft there and natural
472	!!$
473	!!$ ENDDO ! loop over pfts
474	!!$
475	!!$ ENDIF ! sumfpc_wood > fpc_crit
476
477	!
478	! 3.3 update output variables
479	!
480
481	DO j = 2,nvm
482
483	IF ( PFTpresent(i,j) .AND. natural(j) ) THEN
484
485	bm_to_litter(i,j,:) = bm_to_litter(i,j,:) + &
486	biomass(i,j,:) * ( un - survive(j) )
487
488	biomass(i,j,:) = biomass(i,j,:) * survive(j)
489
490	IF ( control%ok_dgvm ) THEN
491	ind(i,j) = ind(i,j) * survive(j)
492	ENDIF
493
494	! fraction of plants that dies each day.
495	! exact formulation: light_death(i,j) = 1. - survive(j) ** (1/dt)
496	light_death(i,j) = ( un - survive(j) ) / dt
497
498	ENDIF ! pft there and natural
499
500	ENDDO ! loop over pfts
501
502	ENDIF ! sumfpc > fpc_crit
503
504	ENDDO ! loop over grid points
505
506	!
507	! 4 recalculate fpc on natural part of grid cell (for next light competition)
508	!
509
510	DO j = 2,nvm
511
512	IF ( natural(j) ) THEN
513
514	!
515	! 4.1 natural PFTs
516	!
517
518	IF ( tree(j) ) THEN
519
520	! 4.1.1 trees: minimum cover due to stems, branches etc.
521
522	DO i = 1, npts
523	!NVMODIF
524	! IF (lai(i,j) == val_exp) THEN
525	! veget_lastlight(i,j) = cn_ind(i,j) * ind(i,j)
526	! ELSE
527	! veget_lastlight(i,j) = &
528	! cn_ind(i,j) * ind(i,j) * &
529	! MAX( ( un - exp( -lai(i,j) * ext_coeff(j) ) ), min_cover )
530	! ENDIF
531	!! veget_lastlight(i,j) = cn_ind(i,j) * ind(i,j)
532	IF (lai(i,j) == val_exp) THEN
533	veget_lastlight(i,j) = cn_ind(i,j) * ind(i,j)
534	ELSE
535	veget_lastlight(i,j) = &
536	cn_ind(i,j) * ind(i,j) * &
537	MAX( ( un - EXP( - lm_lastyearmax(i,j) * sla(j) * ext_coeff(j) ) ), min_cover )
538	ENDIF
539	ENDDO
540
541	ELSE
542
543	! 4.1.2 grasses
544	DO i = 1, npts
545	!NVMODIF
546	! IF (lai(i,j) == val_exp) THEN
547	! veget_lastlight(i,j) = cn_ind(i,j) * ind(i,j)
548	! ELSE
549	! veget_lastlight(i,j) = cn_ind(i,j) * ind(i,j) * &
550	! ( un - exp( -lai(i,j) * ext_coeff(j) ) )
551	! ENDIF
552	!!veget_lastlight(i,j) = cn_ind(i,j) * ind(i,j)
553	IF (lai(i,j) == val_exp) THEN
554	veget_lastlight(i,j) = cn_ind(i,j) * ind(i,j)
555	ELSE
556	veget_lastlight(i,j) = cn_ind(i,j) * ind(i,j) * &
557	( un - exp( - lm_lastyearmax(i,j) * sla(j) * ext_coeff(j) ) )
558	ENDIF
559	ENDDO
560	ENDIF ! tree/grass
561
562	ELSE
563
564	!
565	! 4.2 agricultural PFTs: not present on natural part
566	!
567
568	veget_lastlight(:,j) = zero
569
570	ENDIF ! natural/agricultural
571
572	ENDDO
573
574	ELSE ! static
575
576	light_death(:,:) = zero
577
578	DO j = 2, nvm
579
580	IF ( natural(j) ) THEN
581
582	! 2.1.1 natural PFTs, in the one PFT only case there needs to be no special case for grasses,
583	! neither a redistribution of mortality (delta fpc)
584
585	WHERE( ind(:,j)*cn_ind(:,j) .GT. min_stomate )
586	lai_ind(:)=sla(j) * lm_lastyearmax(:,j) / ( ind(:,j) * cn_ind(:,j) )
587	ELSEWHERE
588	lai_ind(:)=zero
589	ENDWHERE
590
591	fpc_nat(:,j) = cn_ind(:,j) * ind(:,j) * &
592	MAX( ( 1._r_std - exp( - ext_coeff(j) * lai_ind(:) ) ), min_cover )
593
594	WHERE(fpc_nat(:,j).GT.fpc_max(:,j))
595
596	light_death(:,j)=MIN(un,un-fpc_max(:,j)/fpc_nat(:,j))
597
598	ENDWHERE
599
600	DO k=1,nparts
601
602	bm_to_litter(:,j,k)=bm_to_litter(:,j,k)+light_death(:,j)*biomass(:,j,k)
603	biomass(:,j,k)=biomass(:,j,k)-light_death(:,j)*biomass(:,j,k)
604
605	ENDDO
606	ind(:,j)=ind(:,j)-light_death(:,j)*ind(:,j)
607	! if (j==10) print ,'ind10bis=',ind(:,j),light_death(:,j)ind(:,j)
608	ENDIF
609	ENDDO
610
611	light_death(:,:)=light_death(:,:)/dt
612
613	ENDIF
614
615	!
616	! 5 history
617	!
618
619	CALL histwrite (hist_id_stomate, 'LIGHT_DEATH', itime, &
620	light_death, npts*nvm, horipft_index)
621
622	IF (bavard.GE.4) WRITE(numout,*) 'Leaving light'
623
624	END SUBROUTINE light
625
626	END MODULE lpj_light

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

Download in other formats: