Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_stomate/stomate_mark_kill.f90 @ 7852

Last change on this file since 7852 was 7555, checked in by sebastiaan.luyssaert, 2 years ago
Contributes to ticket #837. Simplified forest management for young stands a bit
File size: 52.7 KB

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

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

Download in other formats: