Context Navigation

stomate_mark_kill.f90 @ 7346

Last change on this file since 7346 was 5691, checked in by sebastiaan.luyssaert, 6 years ago

DEV: tested with 13, 37 and 64 PFTs with LCC on different pixels. Some configuration run for 20 years on a given pixel, other crash on another pixel. There is a mass balance problem in sapiens_lcc (ticket #482). This commit fixes a problem with PFT1 in littercalc. This PFT is now fully integrated in LCC and subsequent litter and soil dynamics. veget_max was changed in veget_cov_max where appropriate, a typo in enerbil was corrected.

File size: 44.6 KB

Line
1
2	! =================================================================================================================================
3	! MODULE : stomate_mark_kill
4	!
5	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
6	!
7	! LICENCE : IPSL (2006)
8	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
9	!
10	!>\BRIEF Marks vegetation for killing by natural causes. Replaces
11	! lpj_kill. Does not work with the DGVM.
12	!!
13	!!\n DESCRIPTION: None
14	!!
15	!! RECENT CHANGE(S): None
16	!!
17	!! REFERENCE(S) :
18	!!
19	!! SVN :
20	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_stomate/stomate_mark_kill.f90 $
21	!! $Date: 2013-09-13 17:14:52 +0200 (Fri, 13 Sep 2013) $
22	!! $Revision: 1470 $
23	!! \n
24	!_ ================================================================================================================================
25
26
27	MODULE stomate_mark_kill
28
29	! modules used:
30
31	USE ioipsl_para
32	USE stomate_data
33	USE pft_parameters
34	USE constantes
35	USE function_library, ONLY : nmax, wood_to_qmdia, check_surface_area, &
36	check_mass_balance
37
38	IMPLICIT NONE
39
40	! private & public routines
41
42	PRIVATE
43	PUBLIC mark_to_kill, distribute_mortality_biomass
44
45	CONTAINS
46
47	!! ================================================================================================================================
48	!! SUBROUTINE : mark_to_kill
49	!!
50	!>\BRIEF Kills natural pfts that have low biomass or low number of individuals.
51	!! Does not change the biomass or litter pools.
52	!!
53	!! DESCRIPTION : Kills natural PFTS. Forests can die through self-thinning or
54	!! environmental conditions. Crops can die through environmental
55	!! conditions, and grasses don't die at all. The total number
56	!! of individuals to be killed in each circ class are outputted
57	!! in circ_class_kill, and the biomass is actually killed in
58	!! gap_prognostic.
59	!!
60	!! RECENT CHANGE(S): None
61	!!
62	!! MAIN OUTPUT VARIABLE(S): ::circ_class_kill
63	!!
64	!! REFERENCE(S) : None
65	!!
66	!! FLOWCHART : None
67	!! \n
68	!_ ================================================================================================================================
69
70	SUBROUTINE mark_to_kill (npts, whichroutine, lm_lastyearmax, &
71	PFTpresent, npp_longterm, circ_class_biomass, &
72	circ_class_n, circ_class_kill, rue_longterm, turnover_longterm, &
73	dt, veget_cov_max, st_dist, forest_managed)
74
75	!! 0. Variable and parameter description
76
77	!! 0.1 Input variables
78
79	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
80	CHARACTER(LEN=10), INTENT(in) :: whichroutine !! Message (unitless)
81	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_cov_max !! "maximal" coverage fraction of a PFT on
82	!! the ground (unitless, 0-1)
83	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lm_lastyearmax !! Last year's maximum leaf mass, for each
84	!! PFT @tex $(gC.m^{-2})$ @endtex
85	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: rue_longterm !! Longterm radiation use efficiency
86	!! (??units??)
87	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
88	INTENT(in) :: turnover_longterm !! "Long term" (default 3-year) turnover
89	!! rate @tex $(gC m^{-2} year^{-1})$ @endtex
90	REAL(r_std), INTENT(in) :: dt !! Time step (days)
91	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: npp_longterm !! "Long term" (default = 3-year) net primary
92	!! productivity
93	!! @tex $(gC.m^{-2} year^{-1})$ @endtex
94	REAL(r_std), DIMENSION(:), INTENT(in) :: st_dist !! The probability distribution of trees
95	!! removed from a circ class in case of
96	!! self thinning (unitless).
97
98	!! 0.2 Output variables
99
100
101	!! 0.3 Modified variables
102
103	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: PFTpresent !! Is pft there (true/false)
104	REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), &
105	INTENT(inout) :: circ_class_biomass !! Biomass of the componets of the model
106	!! tree within a circumference
107	!! class @tex $(gC ind^{-1})$ @endtex
108	REAL(r_std), DIMENSION(npts,nvm,ncirc), INTENT(inout) :: circ_class_n !! Number of individuals in each circ class
109	!! @tex $(m^{-2})$ @endtex
110	REAL(r_std), DIMENSION(npts,nvm,ncirc,nfm_types,ncut_times),&
111	INTENT(inout) :: circ_class_kill !! Number of trees within a circ that needs
112	!! to be killed @tex $(ind m^{-2})$ @endtex
113	INTEGER(i_std), DIMENSION (:,:), INTENT(inout) :: forest_managed !! forest management flag (is the forest
114	!! being managed?)
115
116	!! 0.4 Local variables
117
118	INTEGER(i_std) :: ipts, ivm,icir,ipar !! Indices (unitless)
119	INTEGER(i_std) :: iele, imbc, istep !! Indices (unitless)
120	INTEGER(i_std) :: ifm !! Indices (unitless)
121	LOGICAL, DIMENSION(npts) :: was_killed !! Bookkeeping (true/false)
122	REAL(r_std), DIMENSION(ncirc) :: diameters_temp !! Trunk diameter calculated from allometry
123	REAL(r_std) :: Dg !! Quadratic mean diameter of the stand (cm)
124	REAL(r_std) :: bm_difference !! the difference between self-thinning biomass
125	!! and that due to environmental mortality
126	REAL(r_std) :: difference !! the difference between the self-thinning
127	!! density
128	!! and the true density (individuals per m**2)
129	REAL(r_std), DIMENSION(npts,nvm) :: st_dens !! the self-thinning density of the stand
130	!! (individuals per m**2)
131	REAL(r_std), DIMENSION(npts,nvm) :: st_mortality !! biomass killed due to self-thinning density
132	REAL(r_std), DIMENSION(npts,nvm) :: d_mortality !! biomass killed due to the environment
133	REAL(r_std) :: delta_biomass !! Net biomass increase for the previous
134	!! year @tex $(gC m^{-2} year^{-1})$ @endtex
135	REAL(r_std) :: vigour !! Growth efficiency, an indicator of tree
136	!! vitality, used to calculate mortality
137	REAL(r_std) :: mortality_greff !! Mortality rate derived by growth
138	!! efficiency @tex $(year^{-1})$ @endtex
139	REAL(r_std) :: mortality !! Mortality (fraction of trees that is
140	!! dying per time step)
141	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
142	!! mass balance chech for this routine
143	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
144	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
145	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
146	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start, pool_end!! Start and end pool of this routine
147	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
148	REAL(r_std) :: temp
149	LOGICAL :: lconverged !! Checking for loop convergence
150	REAL(r_std), DIMENSION(ncirc) :: target_st_kill !! How many trees from each circumference
151	!! class we would like to kill.
152	!! @tex $(trees m^{-2})$ @endtex
153	REAL(r_std), DIMENSION(ncirc) :: true_st_kill !! How many trees from each circumference
154	!! class that we actually killed.
155	!! @tex $(trees m^{-2})$ @endtex
156	REAL(r_std) :: excess_ind !! THe difference between the amount of
157	!! trees we want to kill and how many we
158	!! can actually kill based on our
159	!! circ_class distribution.
160	!! @tex $(trees m^{-2})$ @endtex
161	REAL(r_std), DIMENSION(npts,nvm) :: veget_cov_max_begin !! temporary storage of veget_cov_max to check area conservation
162
163	!_ ================================================================================================================================
164
165
166	IF (printlev.GE.3) WRITE(numout,*) 'Entering mark_to_kill'
167
168	!! 1. Initialize
169
170	!! 1.1 Initialize local printlev
171	!printlev_loc=get_printlev('mark_to_kill')
172
173	!! 1.2 Initialize check for mass balance closure
174	IF (err_act.GT.1) THEN
175
176	check_intern = zero
177	pool_start = zero
178	DO ipar = 1,nparts
179	DO iele = 1,nelements
180	DO icir = 1,ncirc
181
182	! Initial biomass pool
183	pool_start(:,:,iele) = pool_start(:,:,iele) + &
184	circ_class_biomass(:,:,icir,ipar,iele) * &
185	circ_class_n(:,:,icir) * veget_cov_max(:,:)
186
187	ENDDO
188	ENDDO
189	ENDDO
190
191	!! 1.3 Initialize check for area conservation
192	veget_cov_max_begin(:,:) = veget_cov_max(:,:)
193
194	ENDIF ! err_act.GT.1
195
196	!! 1.4 Initialize variables
197	d_mortality(:,:)=zero
198	st_mortality(:,:)=zero
199
200
201	!! 2. Kill PFTs
202
203	! Kill plants if number of individuals or last year's leaf mass is close to zero.
204	! the "was_killed" business is necessary for a more efficient code on the VPP
205	! processor
206	DO ivm = 2,nvm ! loop plant functional types
207
208	was_killed(:) = .FALSE.
209	IF ( natural(ivm) ) THEN
210
211	!! 2.2 Kill natural PFTs when running STOMATE without DGVM
212
213	!! 2.2.1 Kill PFTs
214	! Kill PFTs but not in the first year because one year
215	! is needed to fill the long-term and seasonal variables.
216	! Therefore no leaves are grown in the first year. As
217	! long as the plant has never grown leaves, ::rue_longterm
218	! keeps its initial value of 1. Note that the logic of
219	! this statement differs from the original verison that
220	! was used for the resource limited allocation where the
221	! carbohydrate reserves pool was treated differently than
222	! in the functional allocation.
223	WHERE ( PFTpresent(:,ivm) .AND. &
224	(rue_longterm(:,ivm) .NE. un) .AND. &
225	(SUM(circ_class_biomass(:,ivm,:,icarbres,icarbon)*circ_class_n(:,ivm,:),2) .LE. min_stomate .AND. &
226	SUM(circ_class_biomass(:,ivm,:,ilabile,icarbon)*circ_class_n(:,ivm,:),2) .LE. min_stomate .AND. &
227	SUM(circ_class_biomass(:,ivm,:,ileaf,icarbon)*circ_class_n(:,ivm,:),2) .LE. min_stomate ) .AND. &
228	SUM(circ_class_n(:,ivm,:),2) .GT. zero)
229
230	was_killed(:) = .TRUE.
231
232	ENDWHERE
233
234	! If the density of a forest goes below a certain level,
235	! the forest dies off. The density limit is referred to
236	! in Bellasen et al 2010.
237	! NOTE: this is currently the only mechanism that can kill
238	! unmanaged forest.
239	IF(is_tree(ivm))THEN
240
241	! The set-up does not allow for species/management
242	! changes. Simulate a stand replacing disturbance.
243	WHERE ( PFTpresent(:,ivm) .AND. &
244	(SUM(circ_class_n(:,ivm,:),2) .LE. dens_target(ivm)*ha_to_m2))
245
246	was_killed(:) = .TRUE.
247
248	ENDWHERE
249
250
251	ENDIF ! checking for trees
252
253	! Debug
254	IF(printlev_loc>4)THEN
255	IF(is_tree(ivm))THEN
256	DO ipts=1,npts
257	IF( PFTpresent(ipts,ivm) .AND. &
258	(SUM(circ_class_n(ipts,ivm,:)) .LE. dens_target(ivm)*ha_to_m2))THEN
259
260	WRITE(numout,*) 'Killing PFT, ipts: ',ivm,ipts
261	WRITE(numout,*) 'Forest density too low'
262	WRITE(numout,*) 'circ_class_n(ipts,ivm,:)',circ_class_n(ipts,ivm,:)
263	WRITE(numout,) 'dens_target(ivm)ha_to_m2',dens_target(ivm)*ha_to_m2
264	ENDIF
265
266
267	ENDDO
268	ENDIF
269
270	DO ipts=1,npts
271	IF ( PFTpresent(ipts,ivm) .AND. &
272	(rue_longterm(ipts,ivm) .NE. un) .AND. &
273	(SUM(circ_class_biomass(ipts,ivm,:,icarbres,icarbon)) .LE. min_stomate .AND. &
274	SUM(circ_class_biomass(ipts,ivm,:,ilabile,icarbon)) .LE. min_stomate .AND. &
275	SUM(circ_class_biomass(ipts,ivm,:,ileaf,icarbon)) .LE. min_stomate ) .AND. &
276	SUM(circ_class_n(ipts,ivm,:)) .GT. zero)THEN
277
278	WRITE(numout,*) 'Killing PFT, ipts: ',ivm,ipts
279	WRITE(numout,*) 'No roots, leaves, labile pool.'
280	WRITE(numout,*) 'circ_class_n(ipts,ivm,:)',SUM(circ_class_n(ipts,ivm,:))
281	WRITE(numout,*) 'Reserve pool',SUM(circ_class_biomass(ipts,ivm,:,icarbres,icarbon))
282	WRITE(numout,*) 'Labile pool',SUM(circ_class_biomass(ipts,ivm,:,ilabile,icarbon))
283	WRITE(numout,*) 'Leaves',SUM(circ_class_biomass(ipts,ivm,:,ileaf,icarbon))
284
285	ENDIF
286	ENDDO
287	ENDIF
288	!-
289
290	!! 2.3 Bookkeeping for PFTs that were killed
291	! For PFTs that were killed, return biomass pools to litter
292	! and update biomass pools to zero
293	DO ipts = 1,npts
294
295	IF ( was_killed(ipts) ) THEN
296
297	! If a PFT is killed, this just means that we mark all the
298	! individuals for a natural death (all debris stays on site).
299	! The killing actually takes place in lpj_gap.f90.
300	! We have to do this differently for trees and non-trees, since
301	! circ_class_n is not defined for non-trees.
302	IF(is_tree(ivm))THEN
303
304	! In reality, if the forest is managed we will not lose this
305	! wood. The forest would be harvested. So, let's put these into
306	! different pools depending on the management type.
307	ifm=forest_managed(ipts,ivm)
308
309	DO icir=1,ncirc
310	circ_class_kill(ipts,ivm,icir,ifm,icut_clear)=&
311	circ_class_n(ipts,ivm,icir)
312	ENDDO
313
314	ELSE
315
316	! Since lpj_gap only knows about circ_class_kill and circ_class_biomass,
317	! we need to give values to icir=1 for these variables for grasses
318	! and crops so that the correct amount of biomass is killed.
319	circ_class_kill(ipts,ivm,1,ifm_none,icut_clear) = &
320	SUM(circ_class_n(ipts,ivm,:))
321
322	ENDIF ! is_tree(ivm)
323
324	!! 2.7 Print sub routine messages
325	IF (printlev_loc>=4) THEN
326
327	WRITE(numout,*) 'kill: eliminated ',PFT_name(ivm)
328	WRITE(numout,*) ' pft number: ',ivm
329	WRITE(numout,*) ' at grid: ',ipts
330
331	ENDIF
332
333	ENDIF ! was_killed
334
335	ENDDO ! loop over points
336
337	ENDIF ! PFT is natural
338
339	!! We need to enforce two different types of killing at the moment,
340	!! self-thinning and environmental mortality. Self-thinning is due
341	!! to resouce competition in dense stands, and environmental
342	!! mortality is due to environmental stress. Self-thinning is
343	!! only applicable to forests.
344	DO ipts=1,npts
345
346	! Don't waste time if the whole PFT was already scheduled to die,
347	! or if the PFT is not here.
348	IF(was_killed(ipts))CYCLE
349	IF(.NOT.PFTpresent(ipts,ivm))CYCLE
350
351	!! let's calculate the self-thinning limit first
352	! Number of indiviudals for trees and grasses are calculated
353	! differently
354	IF ( is_tree(ivm) ) THEN
355
356	! Calculate the self-thinning density. This is the target density
357	! for which we will kill small trees until we reach.
358	! there are unit conversions here because the forestry routines
359	! prefer to work with trees per hectare. I convert everything to
360	! individuals per m**2 when I can. Self-thinning relationship
361	! is fitted for the aboveground biomass/diameter so wood_to_qmdia
362	! should be used here (wood_to_dia_eff accounts for both the above
363	! and belowground biomass).
364	Dg = wood_to_qmdia(circ_class_biomass(ipts,ivm,:,:,icarbon), &
365	circ_class_n(ipts,ivm,:),ivm)
366
367	!!$ Dg=zero
368	!!$ diameters_temp(:) = &
369	!!$ wood_to_dia(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
370	!!$ diameters_temp(:) = &
371	!!$ diameters_temp(:)*m_to_cm
372	!!$
373	!!$ DO icir=1,ncirc
374	!!$ Dg=Dg+circ_class_n(ipts,ivm,icir)m2_to_hadiameters_temp(icir)**2
375	!!$ ENDDO
376
377	! After a coppice at the end of the year, it's possible that all
378	! the aboveground biomass is killed, but we still have individuals
379	! because the roots are alive. This means we don't have
380	! a diameter, though, and therefore we will divide by zero
381	! in the next block of code. So we skip it, which is okay because
382	! there will be no self-thinning in this case anyways.
383	IF(Dg .LT. min_stomate) CYCLE
384
385	!!$ Dg=SQRT(Dg/(ind(ipts,ivm)*m2_to_ha))
386	st_dens(ipts,ivm)=Nmax(Dgm_to_cm,ivm)ha_to_m2
387
388	IF(printlev_loc>=4 .AND. ivm == test_pft .AND. &
389	ipts == test_grid)THEN
390	WRITE(numout,*) 'Does self thinning happen?'
391	WRITE(numout,*) 'ipts,ivm: ',ipts,ivm
392	WRITE(numout,*) 'circ_class_n(ipts,ivm,:),st_dens(ipts,ivm): ',&
393	SUM(circ_class_n(ipts,ivm,:)),st_dens(ipts,ivm)
394	ENDIF
395
396	IF(SUM(circ_class_n(ipts,ivm,:)) .GT. st_dens(ipts,ivm))THEN
397
398	IF(printlev_loc>=4 .AND. ivm == test_pft)THEN
399	WRITE(numout,*) 'Self thinning is taking place.'
400	WRITE(numout,*) 'ipts,ivm: ',ipts,ivm
401	WRITE(numout,*) 'circ_class_n(ipts,ivm,:),st_dens(ipts,ivm): ',&
402	SUM(circ_class_n(ipts,ivm,:)),st_dens(ipts,ivm)
403	WRITE(numout,*) 'difference: ',&
404	SUM(circ_class_n(ipts,ivm,:))-st_dens(ipts,ivm)
405	WRITE(numout,*) 'st_dist: ',&
406	st_dist(:)
407	ENDIF
408
409	! Here is where we actually do the self thinning. The number of
410	! trees that we want to kill is the differnce between the
411	! self thinning density and the true density. Originally, we
412	! killed only small trees until the densities were equal, the
413	! reason being that self-thinning represents competition
414	! driven mortality, and small trees are not able to compete for
415	! resources (e.g. light) as well as larger trees. However, this
416	! seemed to cause strange results for both temperate and tropical
417	! PFTs, so now we kill trees based on a user-defined distribution.
418	difference=SUM(circ_class_n(ipts,ivm,:)) - st_dens(ipts,ivm)
419
420	! This is what we want to kill.
421	target_st_kill(:)=st_dist(:)*difference
422
423	! However, we may not have a sufficient number of trees in each class,
424	! so we might have to adjust things a bit.
425	true_st_kill(:)=target_st_kill(:)
426
427	istep=1
428	DO
429
430	lconverged=.FALSE.
431
432	! Find out how many indivuals we cannot kill with the current
433	! distribution.
434	excess_ind=zero
435	DO icir=1,ncirc
436	IF(circ_class_n(ipts,ivm,icir) .LT. true_st_kill(icir))THEN
437	excess_ind=excess_ind+true_st_kill(icir)-circ_class_n(ipts,ivm,icir)
438	true_st_kill(icir)=circ_class_n(ipts,ivm,icir)
439	ENDIF
440	ENDDO
441
442	IF(excess_ind .GT. zero)THEN
443	! We have something to redistribute.
444	! Put all of the excess in the first class which is not already completely
445	! killed. If we add too much, we'll catch it in the next loop.
446	DO icir=1,ncirc
447	IF(circ_class_n(ipts,ivm,icir) .GT. true_st_kill(icir))THEN
448	true_st_kill(icir)=true_st_kill(icir)+excess_ind
449	excess_ind=zero
450	ENDIF
451	ENDDO
452	ELSE
453	! Done.
454	lconverged=.TRUE.
455	ENDIF
456
457	IF(lconverged)EXIT
458
459	istep=istep+1
460	IF(istep .GT. 100)THEN
461
462	! This is an arbitrary exit here, just to make sure
463	! we don't get stuck in an infinite loop. We should
464	! exit before we ever get here. Otherwise, you can increase
465	! the number in the if statement and try again.
466	WRITE (numout,*) 'ipts, ivm : ',ipts,ivm
467	CALL ipslerr_p (3,'stomate_mark_kill', &
468	'Was not able to distribute self-thinning individuals!', &
469	'See the comments at this part of the code for a solution','')
470
471	END IF
472	ENDDO
473
474	circ_class_kill(ipts,ivm,:,ifm_none,icut_thin)=true_st_kill(:)
475
476	! Debug
477	IF(printlev_loc>=4 .AND. ivm == test_pft)THEN
478	WRITE(numout,*) 'After self thinning.'
479	WRITE(numout,*) 'circ_class_kill: ',&
480	circ_class_kill(ipts,ivm,:,ifm_none,icut_thin)
481	ENDIF
482	!-
483
484	ENDIF ! SUM(circ_class_n(ipts,ivm,:)) .GT. st_dens(ipts,ivm)
485
486	! now let's look at environmental mortality. The calculation of this
487	! rate depends on if it is fixed for the whole simulation or if it's
488	! allowed to respond to things like NPP.
489
490	!! 3.1 Use growth efficiency or constant mortality?
491	IF ( .NOT. ok_constant_mortality ) THEN
492
493	! Was the PFT alive last year?
494	IF ( ( lm_lastyearmax(ipts,ivm) .GT. min_stomate ) ) THEN
495
496	!! 3.1.1 Estimate net biomass increment
497	! Estimate net biomass increment by subtracting turnover from
498	! last year's NPP.
499	! Note that npp_longterm is the longterm growth efficiency
500	! (NPP/LAI) to be fair to deciduous trees
501	delta_biomass = MAX( npp_longterm(ipts,ivm) - &
502	( turnover_longterm(ipts,ivm,ileaf,icarbon) + &
503	turnover_longterm(ipts,ivm,iroot,icarbon) + &
504	turnover_longterm(ipts,ivm,ifruit,icarbon) + &
505	turnover_longterm(ipts,ivm,isapabove,icarbon) + &
506	turnover_longterm(ipts,ivm,isapbelow,icarbon) ), zero)
507
508	!! 3.1.2 Calculate growth efficiency
509	! Calculate growth efficiency by dividing net biomass increment
510	! by last year's maximum LAI. (corresponding actually to the
511	! maximum LAI of the previous year)
512	vigour = delta_biomass / (lm_lastyearmax(ipts,ivm) * sla(ivm))
513
514
515	ELSE
516
517	! PFT does not exist or is dead
518	vigour = zero
519
520	END IF ! ( lm_lastyearmax(ipts,ivm) .GT. min_stomate)
521
522	!! 3.1.3 Calculate growth efficiency mortality rate
523	mortality_greff = mortality_max(ivm) / &
524	( un + ref_mortality(ivm) * vigour )
525
526	! Scale mortality rate (fraction, 0-1) by timesteps per year
527	mortality = MAX(mortality_min(ivm), mortality_greff) * dt/&
528	one_year
529
530	ELSE ! .NOT. ok_constant_mortality
531
532	!! 3.2 Use constant mortality accounting for the residence time
533	!! of each tree PFT
534
535	! Scale mortality rate (fraction, 0-1) by timesteps per year
536	! NOTE: the meaning of residence_time is very different between
537	! the DOFOCO branch and the trunk. In the trunk biomass has no age
538	! thus the residence time accounts for all forest dynamics including
539	! self-thinning, pests, diseases and windthrow. In the DOFOCO branch
540	! biomass does have an age and self-thinning is explicitly accounted
541	! for, hence, the residence time should be much higher as it only
542	! accounts for pest, diseases and windthrow. Even the latter is not
543	! exact because as long as those disturbances are small scale they
544	! are probably accounted for in the parametrization of self-thinning.
545	mortality = dt/(residence_time(ivm)*one_year)
546
547	ENDIF ! .NOT. ok_constant_mortality
548
549	! Debug
550	IF (printlev_loc>=4 .AND. ivm .EQ. test_pft) THEN
551	WRITE(numout,*) 'Mortality rate (0-1), ', mortality
552	ENDIF
553	!-
554
555	! Now the mortality is just a percentage of the total biomass
556	! on the site. Since the main variable we output here is
557	! circ_class_kill, we need to convert that biomass into the
558	! number of individuals that are killed, and from which
559	! circumference classes. Mortality in excess of the self-thinning
560	! mortality will preferenctially kill large trees, since they are
561	! assumed to be less resistant to climatic change.
562
563	! Aggregate over the different biomass components
564	DO ipar = 1,nparts
565
566	! Aggregate over the different circumference classes
567	DO icir = 1,ncirc
568
569	! Mortality from self-thinning
570	! Note that the unit of st_mortaility (g C m-2 dt-1) differs
571	! from that of mortality (rate dt-1)
572	st_mortality(ipts,ivm) = st_mortality(ipts,ivm) + &
573	circ_class_kill(ipts,ivm,icir,ifm_none,icut_thin) * &
574	circ_class_biomass(ipts,ivm,icir,ipar,icarbon)
575
576	ENDDO
577
578	! Ageing and environmental mortality
579	! Note that the unit of d_mortaility (g C m-2 dt-1) differs from
580	! that of mortality (rate dt-1)
581	d_mortality(ipts,ivm) = d_mortality(ipts,ivm) + &
582	mortality * SUM(circ_class_biomass(ipts,ivm,:,ipar,icarbon)*&
583	circ_class_n(ipts,ivm,:))
584
585	ENDDO ! ipar = 1,nparts
586
587	IF(printlev_loc>=4 .AND. test_pft .EQ. ivm .AND. &
588	test_grid == ipts)THEN
589	WRITE(numout,*) 'Does mortality happen?'
590	WRITE(numout,*) 'ipts,ivm: ',ipts,ivm
591	WRITE(numout,*) 'd_mortality(ipts,ivm): ',d_mortality(ipts,ivm)
592	WRITE(numout,*) 'st_mortality(ipts,ivm): ',st_mortality(ipts,ivm)
593	ENDIF
594
595
596	!+++CHECK+++
597	! There are different ways to combine these two sources of
598	! mortality. First one could argue that self-thinning
599	! includes the background mortality. If background mortality
600	! exceeds self-thinning this is due to the simplified way
601	! background mortality is calculated.
602	IF (st_mortality(ipts,ivm) .GT. min_stomate) THEN
603
604	! There is self-thinning mortality, no need to account
605	! for environmental mortality
606	bm_difference = zero
607
608	ELSE
609
610	! There is no self_thinning mortality, use the
611	! environmental mortality
612	bm_difference = d_mortality(ipts,ivm)
613
614	ENDIF
615
616	!!$ ! The second way of looking at this is that the calculation of
617	!!$ ! background mortality is reliable and should therefore always
618	!!$ ! be taken into account. In this reasoning self-thinning is
619	!!$ ! considered part of the background mortality. Account for both
620	!!$ ! but avoid double counting.
621	!!$ bm_difference=d_mortality(ipts,ivm) - st_mortality(ipts,ivm)
622	!+++++++++++
623
624
625	!+++CHECK +++
626	! The question still stands of which trees we kill with the
627	! environmental mortality. If we only kill the biggest, we
628	! lose a circ class within the first year sometimes.
629	IF(bm_difference .GT. zero )THEN
630
631	! We know the total biomass that we want to kill. Let's partition
632	! this biomass among all our circumference classes based on an
633	! exponential distribution which takes more biomass from
634	! the largest class than the smallest. There is a chance that we
635	! will not have enough biomass to kill in some classes, in which
636	! case we need to rearrange it a bit.
637
638	! Debug
639	IF(printlev_loc>=4 .AND. ivm == test_pft .AND. &
640	ipts == test_grid)THEN
641	WRITE(numout,*) 'Mortality is taking place.'
642	WRITE(numout,*) 'ipts,ivm: ',ipts,ivm
643	WRITE(numout,*) 'bm_difference: ',bm_difference
644	WRITE(numout,*) 'circ_class_kill(ipts,ivm,:,ifm_none,icut_thin): ', &
645	circ_class_kill(ipts,ivm,:,ifm_none,icut_thin)
646	WRITE(numout,*) 'circ_class_n(ipts,ivm,:): ', &
647	circ_class_n(ipts,ivm,:)
648	ENDIF
649	!-
650
651	CALL distribute_mortality_biomass( bm_difference, &
652	death_distribution_factor(ivm), circ_class_n(ipts,ivm,:), &
653	circ_class_biomass(ipts,ivm,:,:,icarbon), &
654	circ_class_kill(ipts,ivm,:,ifm_none,icut_thin) )
655
656	IF(printlev_loc>=4 .AND. ivm == test_pft .AND. &
657	ipts == test_grid)THEN
658	WRITE(numout,*) 'after circ_class_kill(ipts,ivm,:,ifm_none,icut_thin): ', &
659	circ_class_kill(ipts,ivm,:,ifm_none,icut_thin)
660	WRITE(numout,*) 'circ_class_n(ipts,ivm,:): ', &
661	circ_class_n(ipts,ivm,:)
662	ENDIF
663
664	ENDIF ! bm_difference .GT. zero
665
666	ELSE ! IF(is_tree(ivm))
667
668	! we don't seem to use self thinning for grasses, since grasses have
669	! a constant density of 1/m**2. This is the original code for the
670	! mortality. Notice that it doesn't actually do anything, due to
671	! Sebastiaan's changes and CHECK statement, so I will leave it
672	! commented out for now.
673
674	! one thing I will set is
675	! circ_class_kill(ipts,ivm,1,ifm_none,icut_thin) and
676	! circ_class_biomass(:,:,1,:,:), since if we decide to
677	! use grass/crop mortality later on we will use these variables
678	! to tell lpj_gap how much biomass to actually kill.
679	circ_class_kill(ipts,ivm,1,ifm_none,icut_thin)=zero
680
681	!!$ IF ( ok_dgvm .OR. .NOT. ok_constant_mortality) THEN
682	!!$
683	!!$ ! +++CHECK+++
684	!!$ ! For grasses, if last year's NPP is very small (less than 10 gCm^{-2}year{-1})
685	!!$ ! the grasses completely die. Grasses and crops are deciduous so in the first
686	!!$ ! year they have no leaves. Consequently npp_longterm is less then
687	!!$ ! ::npp_longterm_init. At the beginning of the second year when rue_longterm is
688	!!$ ! no longer 1, the PFT is killed immediatly. The IF statement needs to be refined
689	!!$ ! to properly allow for this type of PFT killing.
690	!!$ ! The more fundamental and conceptual question is whether we want to kill grasses?
691	!!$ ! Unless we simulate desertification, it is not clear why we would kill the grass
692	!!$ ! most likely to replace it by ... grass. Because grass has hardly any biomass it
693	!!$ ! seem unlikely that killing the grass will affect the C-fluxes.
694	!!$ ! For both reasons i.e. the technical and conceptual, it was decided NOT to kill
695	!!$ ! the grass PFTs
696	!!$ IF ( PFTpresent(ipts,ivm) .AND. ( npp_longterm(ipts,ivm) .LE. npp_longterm_init ) ) THEN
697	!!$
698	!!$ ! ORIGINAL code - see above to justify the change
699	!!$ !!$mortality = un
700	!!$ mortality = zero
701	!!$
702	!!$ END IF
703	!!$ ! +++++++++++
704	!!$
705	!!$ ! Update biomass and litter pools
706	!!$ DO ielem = 1,nelements
707	!!$
708	!!$ DO ipart = 1, nparts
709	!!$
710	!!$ IF ( PFTpresent(ipts,ivm) ) THEN
711	!!$
712	!!$ d_mortality(ipts,ivm) = mortality * biomass(ipts,ivm,ipart,ielem)
713	!!$ bm_to_litter(ipts,ivm,ipart,ielem) = bm_to_litter(ipts,ivm,ipart,ielem) + &
714	!!$ d_mortality(ipts,ivm)
715	!!$ biomass(ipts,ivm,ipart,ielem) = biomass(ipts,ivm,ipart,ielem) - &
716	!!$ d_mortality(ipts,ivm)
717	!!$
718	!!$ END IF ! PFTpresent(ipts,ivm)
719	!!$
720	!!$ ENDDO ! nparts
721	!!$
722	!!$ END DO ! nelements
723	!!$
724	!!$ ENDIF ! .NOT.ok_dgvm .AND. .NOT.ok_constant_mortality
725
726	ENDIF ! end check for trees
727
728	ENDDO ! loop over grid points
729
730	ENDDO ! loop over PFTs
731
732
733	!! 4. Check numerical consistency of this routine
734
735	IF (err_act.GT.1) THEN
736
737	! 4.2 Check surface area
738	CALL check_surface_area("stomate_mark_to_kill", npts, veget_cov_max_begin, &
739	veget_cov_max,'pft')
740
741	! 4.3 Mass balance closure
742	! 4.3.1 Calculate final biomass
743	pool_end = zero
744	DO ipar = 1,nparts
745	DO iele = 1,nelements
746	DO icir = 1,ncirc
747	pool_end(:,:,iele) = pool_end(:,:,iele) + &
748	circ_class_biomass(:,:,icir,ipar,iele) * &
749	circ_class_n(:,:,icir) * veget_cov_max(:,:)
750	ENDDO
751	ENDDO
752	ENDDO
753
754	!! 4.2 Calculate components of the mass balance
755	! 4.3.2 Calculate mass balance
756	! Common processes
757	DO iele=1,nelements
758	check_intern(:,:,ipoolchange,iele) = -un * (pool_end(:,:,iele) - &
759	pool_start(:,:,iele))
760	ENDDO
761
762	closure_intern = zero
763	DO imbc = 1,nmbcomp
764	DO iele=1,nelements
765	! Debug
766	IF (printlev_loc>=4) WRITE(numout,*) &
767	'check_intern, ivm, imbc, iele, ', imbc, &
768	iele, check_intern(:,test_pft,imbc,iele)
769	!-
770	closure_intern(:,:,iele) = closure_intern(:,:,iele) + &
771	check_intern(:,:,imbc,iele)
772	ENDDO
773	ENDDO
774
775	! 4.3.3 Check mass balance closure
776	CALL check_mass_balance("stomate_mark_to_kill", closure_intern, npts, &
777	pool_end, pool_start, veget_cov_max, 'pft')
778
779	ENDIF ! err_act.GT.1
780
781	IF (printlev>=4) WRITE(numout,*) 'Leaving mark_to_kill'
782
783	END SUBROUTINE mark_to_kill
784
785	!! ================================================================================================================================
786	!! SUBROUTINE : distribute_mortality_biomass
787	!!
788	!>\BRIEF Distributes biomass that is going to be killed by natural
789	!! causes (not self thinning) over circ classes.
790	!!
791	!! DESCRIPTION : Mortality is going to kill a certain amount of biomass
792	!! in forests every day. Since we have circumference classes
793	!! now for our forests, we need to determine which classes
794	!! of trees will suffer from this environmental mortality.
795	!! Right now we are taking an exponential distribution.
796	!! Notice that this is NOT the same as redistributing biomass
797	!! after one of the circ classes becomes empty.
798	!!
799	!! RECENT CHANGE(S): None
800	!!
801	!! MAIN OUTPUT VARIABLE(S): ::circ_class_kill
802	!!
803	!! REFERENCE(S) : None
804	!!
805	!! FLOWCHART : None
806	!! \n
807	!_ ================================================================================================================================
808	SUBROUTINE distribute_mortality_biomass ( bm_difference, ddf_temp, circ_class_n_temp, &
809	circ_class_biomass_temp, circ_class_kill_temp )
810
811	!! 0. Variable and parameter description
812
813	!! 0.1 Input variables
814	REAL(r_std),INTENT(in) :: bm_difference !! the biomass to distribute
815	REAL(r_std),INTENT(in) :: ddf_temp !! the death_distribution_factor for this pft
816	REAL(r_std),DIMENSION(:),INTENT(in) :: circ_class_n_temp !! circ_class_n for this point/PFT
817	REAL(r_std),DIMENSION(:,:),INTENT(in) :: circ_class_biomass_temp !! circ_class_biomass for this point/PFT
818
819	!! 0.2 Output variables
820
821	!! 0.3 Modified variables
822	REAL(r_std), DIMENSION(ncirc), INTENT(OUT) :: circ_class_kill_temp !! Number of trees within a circ that needs
823	!! to be killed @tex $(ind m^{-2})$ @endtex
824
825	!! 0.4 Local variables
826	REAL(r_std), DIMENSION(ncirc) :: biomass_desired !! The biomass that dies naturally
827	REAL(r_std), DIMENSION(ncirc) :: death_distribution !! The fraction of biomass taken from
828	!! each circ class for mortality.
829	LOGICAL :: ldone !! Flag to exit a loop
830	REAL(r_std) :: scale_factor !!
831	REAL(r_std) :: sum_total !!
832	REAL(r_std) :: leftover_bm !! excess biomass that we need to kill
833	REAL(r_std) :: living_biomass !! summed biomass
834	REAL(r_std) :: living_trees !! summed biomass
835	INTEGER :: icir
836
837	!_ ================================================================================================================================
838
839	IF(ncirc == 1)THEN
840
841	! This is the easy case. All of our biomass will be taken from the only
842	! circumference class that we have.
843	circ_class_kill_temp(1)=circ_class_kill_temp(1)+&
844	bm_difference/SUM(circ_class_biomass_temp(1,:))
845
846	RETURN
847	ENDIF
848
849	! Here we assume an exponential distribution, arranged so that
850	! ddf_temp times more biomass is taken from the largest circ class
851	! compared to the smallest.
852	biomass_desired(:)=zero
853	death_distribution(:)=un
854	scale_factor=ddf_temp**(un/REAL(ncirc-1))
855	DO icir=2,ncirc
856	death_distribution(icir)=death_distribution(icir-1)*scale_factor
857	ENDDO
858	! Normalize it
859	sum_total=SUM(death_distribution(:))
860	death_distribution(:)=death_distribution(:)/sum_total
861
862	! Ideally, how much is killed from each class? Be careful to include
863	! what was killed in self-thinning here! If we don't, we may try to kill
864	! more biomass than is available in the loop below.
865	DO icir=1,ncirc
866	biomass_desired(icir)=death_distribution(icir)*bm_difference+&
867	circ_class_kill_temp(icir)*SUM(circ_class_biomass_temp(icir,:))
868	ENDDO
869
870	! Right now, we know how much biomass we want to kill in each
871	! class (biomass_desired). What we will do now is to loop through
872	! all the circ_classes and see if we have this much biomass
873	! in each class still alive. The total amount of vegetation still
874	! alive is circ_class_n_temp(icir)-circ_class_kill_temp(icir). If
875	! this number of individuals cannot give us the total biomass that
876	! we need, we keep track of a residual quantity, leftover_bm, which
877	! we try to take from other circ_classes. It's possible that we
878	! will have to loop several times, which makes things more complicated.
879
880	ldone=.FALSE.
881	leftover_bm=zero
882	DO
883	DO icir=ncirc,1,-1
884
885	living_trees=circ_class_n_temp(icir)-circ_class_kill_temp(icir)
886	living_biomass=&
887	SUM(circ_class_biomass_temp(icir,:))*living_trees
888	biomass_desired(icir)=biomass_desired(icir)+leftover_bm
889
890	IF(living_biomass .LE. biomass_desired(icir))THEN
891
892	! We can't get everything from this class, so we kill whatever is left.
893	leftover_bm=biomass_desired(icir)-living_biomass
894	biomass_desired(icir)=zero
895	circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+&
896	living_trees
897
898	ELSE
899
900	! We have enough in this class.
901	circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+&
902	biomass_desired(icir)/SUM(circ_class_biomass_temp(icir,:))
903	biomass_desired(icir)=zero
904	leftover_bm=zero
905
906	ENDIF ! living_biomass .LE. biomass_desired(icir)
907
908	ENDDO ! loop over circ classes
909
910	IF(leftover_bm .LE. min_stomate) EXIT
911
912	! it's possible that we don't have enough biomass left to kill what needs to be
913	! killed, so everything just dies. I cannot think of a case where this
914	! would happen, though, since mortality should always be a percentage of the
915	! total biomass.
916	ldone=.TRUE.
917	DO icir=1,ncirc
918
919	IF( circ_class_kill_temp(icir) .LT. circ_class_n_temp(icir) ) &
920	ldone=.FALSE.
921
922	ENDDO
923
924	IF(ldone)EXIT ! All our biomass is dead, and we still want to kill more!
925
926	ENDDO
927
928	END SUBROUTINE distribute_mortality_biomass
929
930	END MODULE stomate_mark_kill

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

Context Navigation

source: branches/publications/ORCHIDEE_CN_CAN_r5698/src_stomate/stomate_mark_kill.f90 @ 7346

Download in other formats: