Context Navigation

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

Last change on this file since 7852 was 7518, checked in by sebastiaan.luyssaert, 2 years ago
Contributes to tickets #633 and #651
File size: 50.2 KB

Line
1	! ============================================================================================================================
2	! MODULE : stomate_pest
3	!
4	! CONTACT : orchidee-help _at_ ipsl.listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006).
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF : Calculates tree damage due to bark-beetle outbreaks
10	!!
11	!!\n DESCRIPTION: See subroutine bark_beetle_damage for a more elaborate description of the
12	!! principles.
13	!!
14	!! RECENT CHANGE(S): None
15	!!
16	!! REFERENCES: C. TEMPERLI, H. BUGMANN, AND C. ELKIN. Cross-scale interactions
17	!! among bark beetles, climate change, and wind disturbances: a
18	!! landscape modeling approach. Ecological Monographs, 83(3), 2013, pp. 383â402.
19	!!
20	!! SVN :
21	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_stomate/stomate_pest.f90
22	!! $Date: 2019-01-22 13:41:13 +0100 (mar., 06 aout. 2019)
23	!! $Revision: $
24	!! \n
25	!_ ===========================================================================================================================
26
27	MODULE stomate_pest
28
29	! modules used:
30	USE ioipsl_para
31	USE xios_orchidee
32	USE stomate_data
33	USE constantes
34	USE constantes_soil
35	USE pft_parameters
36	USE pft_parameters_var
37	USE function_library, ONLY: wood_to_ba, wood_to_qmdia, Nmax, wood_to_dia
38	IMPLICIT NONE
39
40	! private & public routines
41
42	PRIVATE
43	PUBLIC pest_clear,bark_beetle_damage
44
45	LOGICAL, SAVE :: firstcall_bark_beetle_damage = .TRUE. !! First call
46	!$OMP THREADPRIVATE(firstcall_bark_beetle_damage)
47	INTEGER(i_std), SAVE :: printlev_loc !! Local level of text output
48	!$OMP THREADPRIVATE(printlev_loc)
49
50	CONTAINS
51
52	!!
53	!==================================================================================================================================
54	!! SUBROUTINE : pest_clear
55	!!
56	!>\BRIEF Set the firstcall flag to .TRUE. and activate initialization
57	!!
58	!_
59	!================================================================================================================================
60
61	SUBROUTINE pest_clear
62	firstcall_bark_beetle_damage = .TRUE.
63	END SUBROUTINE pest_clear
64
65
66	!================================================================================================================================
67	!! SUBROUTINE : bark_beetle_damage
68	!!
69	!>\BRIEF Calculates tree damage due to bark-beetle outbreaks
70	!!
71	!!\n DESCRIPTION : This subroutine is used to calculate the annual amount of trees kill by
72	!! bark beetles. This module is strongly inspired by the work of Temperli
73	!! et al. 2013 in Ecological Monographs. The bark beetle module developped
74	!! for LANDCLIM landscape model was the basis of this work. The module was
75	!! adapted to fit the constraints of ORCHIDEE in regard to its spatial and
76	!! temporal resolution. In short, this module calculate a suceptibility of
77	!! of a forest to be attack by bark beetles at the pixel/PFT level. At
78	!! present the model has been parameterized for Picea abies age. The risk
79	!! probability of a pixel to be attacked depends on pixel suceptibility and
80	!! bark beetle presure where pressure is calculated as an index that
81	!! indicates the activity of the beetles in this pixel). Then, a mortality
82	!! rate for each Picea abies PFT is calculated based on PFT suceptibility
83	!! and beetle pressure. Finally, the wood biomass kill by bark beetle of a
84	!! given PFT is the product of the pixel risk * mortality rate * the PFT
85	!! wood biomass. For each PFT we first kill the small trees and if more
86	!! biomass needs to be killed, ORCHIDEE will continue with killing the
87	!! larger diameter classes.
88	!! The most important difference beetween this module and the original one
89	!! in LANDCLIM is that ORCHIDEE lost all information about the spatialization
90	!! of the attack. In other words the attack is no longer spatialy explicit.
91	!!
92	!! RECENT CHANGE(S): Added in December 2019
93	!!
94	!! MAIN OUTPUT VARIABLE(S) : circ_class_kill
95	!!
96	!! REFERENCES: CHRISTIAN TEMPERLI, HARALD BUGMANN, AND CHE ELKIN. Cross-scale
97	!! interactions among bark beetles, climate change, and wind
98	!! disturbances: a landscape modeling approach. Ecological Monographs,
99	!! 83(3), 2013, pp. 383â402.
100	!!
101	!! FLOWCHART :
102	!! \n
103	!_================================================================================================================================
104
105	SUBROUTINE bark_beetle_damage( &
106	npts, circ_class_biomass, beetle_generation_index,&
107	veget_max, circ_class_n, circ_class_dist, &
108	age_stand, season_drought_legacy, resolution, &
109	circ_class_kill, forest_managed, wood_leftover_legacy, &
110	beetle_pop_legacy, risk_index_legacy, &
111	beetle_damage_legacy, risk_index, epidemic_monitor, &
112	epidemic, beetle_flyaway, windthrow_suscept_monitor, &
113	beetle_pressure_monitor, suscept_index_monitor )
114
115	IMPLICIT NONE
116
117	!! 0. Variable declarations
118	!! 0.1 Input variables
119	INTEGER(i_std), INTENT(in) :: npts !! number of pixels
120	INTEGER(i_std), DIMENSION (:, :), INTENT(in) :: forest_managed !! forest management flag
121	!! (0 : unmanaged, >0 : managed)
122	REAL(r_std), DIMENSION(:, :, :, :, :), INTENT(in) :: circ_class_biomass !! Biomass of an individual tree (gC.m-2)
123	!! in each circ class
124	REAL(r_std), DIMENSION(:, :, :), INTENT(in) :: circ_class_n !! Number of trees(tree.m-2)
125	!! within each circumference
126	REAL(r_std), DIMENSION(:), INTENT(in) :: circ_class_dist !! The probability distribution of trees
127	REAL(r_std), DIMENSION(:, :), INTENT(in) :: veget_max !! "maximal" coverage fraction of
128	!! a PFT on the ground
129	REAL(r_std), DIMENSION(:, :), INTENT(in) :: resolution !! The length in m of each grid-box in
130	!! X and Y axis
131	REAL(r_std), DIMENSION(:, :, :), INTENT(in) :: beetle_generation_index !! number of generation that BB can
132	!! achieved in one year
133	INTEGER(i_std), DIMENSION(:, :), INTENT(in) :: age_stand !! Age of stand (years)
134	REAL(r_std), DIMENSION(:, :, :), INTENT(in) :: season_drought_legacy !! mean growing season moisture availability
135	REAL(r_std), DIMENSION(:, :, :), INTENT(in) :: wood_leftover_legacy !! Coarse woody bedris, snags and logs
136	!! (gC.m-2) left on the PFT
137	!! during previous timestep
138	REAL(r_std), DIMENSION(:, :, :), INTENT(in) :: risk_index_legacy !! risk index stored every year, n time
139	!! 0.3 Modified variables
140	REAL(r_std), DIMENSION(:, :, :, :, :), INTENT(inout) :: circ_class_kill !! Number of trees within a circ that needs
141	REAL(r_std), DIMENSION(:, :, :),INTENT(inout) :: beetle_pop_legacy !! Index (0-1) Biomass of tree from the
142	!! same PFT that was
143	REAL(r_std), DIMENSION(:, :), INTENT(inout) :: epidemic !! logical flag indicationg if beetle
144	!! outbreak is in a epidemic stage
145	REAL(r_std), DIMENSION(:, :),INTENT(inout) :: risk_index !! risk that a given pixel experience an
146	!! outbreak (unitless)
147	REAL(r_std), DIMENSION(:, :), INTENT(inout) :: beetle_flyaway
148	REAL(r_std), DIMENSION(:, :), INTENT(inout) :: epidemic_monitor !! variable used to monitor epidemic flag
149	REAL(r_std), DIMENSION(:, :, :), INTENT(inout) :: beetle_damage_legacy !! temporary variable(gC.m-2)
150	!! infected during the previous timestep
151	REAL(r_std), DIMENSION(:, :), INTENT(inout) :: windthrow_suscept_monitor!! variable used to monitor windthrow susceptibility (unitless)
152	REAL(r_std), DIMENSION(:, :), INTENT(inout) :: beetle_pressure_monitor !! variable used to monitor beetle pressure index (unitless)
153	REAL(r_std), DIMENSION(:, :), INTENT(inout) :: suscept_index_monitor !! variable used to monitor overall susceptibility index (unitless)
154
155	!! 0.4 Local variables
156	INTEGER(i_std) :: ipts, ileg, ivm, icir, &
157	ivmap, iagec, iyear, xvm, &
158	ncirc_loc
159	INTEGER(i_std), DIMENSION(npts, nvmap, nagec) :: nc !! number of PFT counter
160	!! within a species group that are
161	!! infected by bark beetles (unitless)
162	REAL(r_std), DIMENSION(npts, nvm) :: season_drought_mean !! average drought index on this PFT
163	!! during previous time steps.
164	REAL(r_std), DIMENSION(npts, nvm) :: wood_leftover_mean !! average biomass (gC m-2)
165	!! of dead wood left on this PFT
166	!! during previous time steps.
167	REAL(r_std), DIMENSION(npts, nvmap) :: windthrow_suscept_gp !! Species group-level susceptibility
168	!! to bark beetle attacks due to the
169	!! leftover of windthrow debris (unitless)
170	REAL(r_std), DIMENSION(npts, nvmap) :: drought_suscept_gp !! Species group-level susceptibility
171	!! to bark beetle attacks due to weakening
172	!! from droughts (unitless)
173	REAL(r_std), DIMENSION(npts, nvmap) :: standage_suscept_gp !! Species group-level susceptibility
174	!! to bark beetle attacks
175	REAL(r_std), DIMENSION(npts, nvmap) :: rdi_suscept_gp !! Species group-level susceptibility
176	!! to relative density of stem
177	!! due to the stand age (unitless)
178	REAL(r_std), DIMENSION(npts, nvmap) :: share_suscept_gp !! Species group-level susceptibility
179	!! to bark beetle attacks due to
180	!! the landscape heterogeneity (unitless)
181	REAL(r_std), DIMENSION(npts, nvmap) :: suscept_index_gp !! Species group-level susceptibility
182	!! to bark beetle attacks due to
183	!! the combined effect of windthrow,
184	!! stand age and
185	REAL(r_std), DIMENSION(npts, nvmap) :: risk_index_gp !! risk index averaged for all
186	!! age classes of the same PFT
187	!! droughts (unitless)
188	REAL(r_std), DIMENSION(npts, nvmap) :: beetle_pressure_index !! Beetle activity of a given pixel
189	REAL(r_std), DIMENSION(npts, nvmap) :: biomass_infected !! Total biomass of a species group
190	!! that is infected by bark
191	!! beetles (gC or gN m-2).
192	REAL(r_std), DIMENSION(npts, nvmap) :: vegetmax_gp !! share of all PFTs of the same species
193	!! (called species group) within
194	REAL(r_std), DIMENSION(npts) :: vegetmax_deciduous !! share of deciduous PFTs
195	!! a pixel (0-1, unitless)
196	REAL(r_std), DIMENSION(npts, nvm) :: standage_suscept !! Pft-level susceptibility
197	!! due to stand age (0-1, unitless)
198	REAL(r_std), DIMENSION(npts, nvm) :: drought_suscept !! Pft-level susceptibility
199	!! due to weakening from droughts
200	!! (0-1, unitless)
201	REAL(r_std), DIMENSION(npts, nvm) :: suscept_index !! Pft-level susceptibility to bark beetle
202	!! attacks due to the combined effect of
203	!! windthrow, stand age and
204	!! droughts (0-1, unitless)
205	REAL(r_std), DIMENSION(npts, nvm) :: mortality_rate !! probability that an infected tree die
206	REAL(r_std), DIMENSION(npts, nvm) :: beetle_damage !! bark beetle damage in gC m-2 year-1
207	REAL(r_std), DIMENSION(npts, nvmap) :: G !! Average bark beetle pressure during
208	!! the three recent years
209	REAL(r_std) :: wght_sirdi !! weigth rdi_susceptibility in the calculation
210	!! of susceptibility_index
211	REAL(r_std) :: wght_siw !! weigth windthrow_susceptibility in the
212	!! calculation of susceptibility_index
213	REAL(r_std) :: qm_dia !! quadratic diameter
214	REAL(r_std), DIMENSION(npts, nvmap) :: woody_frac !! woody/total biomass (unitless)
215	REAL(r_std), DIMENSION(ncirc) :: qdia !! quadratic diameter per circumference class
216	REAL(r_std), DIMENSION(npts, nvm) :: rdi !! relative density index
217	REAL(r_std), DIMENSION(npts, nvmap) :: rdi_gp !! rdi for all age-classe of a same PFT
218	REAL(r_std), DIMENSION(npts, nvmap) :: epidemic_risk !! with the highest bark beetle pressure
219	!! (0-1, unitless)
220	REAL(r_std), DIMENSION(legacy_years) :: bgi_current !! temporary variable
221	REAL(r_std), DIMENSION(npts, nvmap) :: biomass_current !! Actual biomass of the PFT (gC or gN m-2)
222	REAL(r_std), DIMENSION(npts, nvmap) :: woody_current !! Actual aboveground woody biomass of the PFT
223	!! (gC or gN m-2)
224	REAL(r_std) :: share !! Deciduous/coniferous share in a pixel
225	REAL(r_std), DIMENSION(npts, nvm, legacy_years) :: beetle_temp1 !! temporary variable
226	REAL(r_std), DIMENSION(npts, nvm, beetle_legacy) :: beetle_temp2
227	REAL(r_std), DIMENSION(nvm) :: rdi_limit !! Rdi a which epidemic come back to endimic (unitless)
228	LOGICAL, DIMENSION(npts, nvmap) :: beetle_gp !! Flag that indicates whether or not
229	!! bark beetle damage should
230	!! be calculated for that species group.
231	!================================================================================================================================
232
233	IF (firstcall_bark_beetle_damage) THEN
234	!! Initialize local printlev
235	printlev_loc = get_printlev('barkbeetle_damage')
236	firstcall_bark_beetle_damage = .FALSE.
237	END IF
238
239	IF (printlev_loc .GE. 2) WRITE(numout, *) 'Entering stomate pest damage'
240	! Initialize variables
241	vegetmax_gp(:, :) = zero
242	standage_suscept_gp(:, :) = zero
243	windthrow_suscept_gp(:, :) = zero
244	drought_suscept_gp(:, :) = zero
245	share_suscept_gp(:, :) = zero
246	rdi_suscept_gp(:, :) = zero
247	risk_index(:, :) = zero
248	rdi(:, :) = zero
249	rdi_gp(:, :) = zero
250	biomass_current(:, :) = zero
251	woody_current(:, :) = zero
252	woody_frac(:, :) = un
253	season_drought_mean(:, :) = zero
254	wood_leftover_mean(:, :) = zero
255	beetle_damage(:, :) = zero
256	mortality_rate(:, :) = zero
257	G(:, :) = zero
258	nc(:, :, :) = zero
259	! THe end of an epidemic is trigger by the forest density (RDI)
260	! rdi_limit must be above rdi_susceptibility_b in order to avoid
261	! long tail in epidemic phase the rdi_target_suscept must be ranged
262	! between 0.5 < rdi_target_suscept < 0.7.
263	rdi_limit(:) = log(1 / rdi_target_suscept(:) - 1) / &
264	rdi_susceptibility_a(:) + &
265	rdi_susceptibility_b(:)
266
267	! Assume that we are not interested in calculating bark beetle
268	! damage for this species group.
269	beetle_gp(:, :) = .FALSE.
270
271	! This is the main loop for wich the beetle damage will be calculated
272	DO ipts = 1, npts
273
274	DO ivmap = 1, nvmap
275
276	! a fist loop over PFT to setup crusial variable: vegetmax_group,
277	! vegetmax_deciduous, rdi_group, beetle_group and nc
278	DO iagec = 1, nagec_pft(ivmap)
279
280	! Determine the PFT index for the age classes within the species
281	! group under consideration
282	ivm = start_index(ivmap) + (iagec - 1)
283
284	WRITE(numout,*) 'PFT, ivmap, ipts', ivm, ivmap, ipts
285	IF(printlev_loc >= 4)THEN
286	WRITE(numout,*) 'PFT', ivm
287	WRITE(numout,*) 'beetle_pft', beetle_pft(ivm)
288	WRITE(numout,*) 'beetle_group', beetle_gp(ipts, ivmap)
289	WRITE(numout,*) 'veget_max', veget_max(ipts, ivm)
290	END IF
291	!-
292
293	IF (.NOT. beetle_pft(ivm) .OR. (veget_max(ipts, ivm) .LE. min_stomate)) THEN
294	! We don't want to calculate bark beetle infections for this PFT
295	! or there is no vegetation present in the PFT. nc is a counter.
296	! Because at least one PFT of the age class group has no beetles in it,
297	! the counter is decreased. It will be used later in other calculations.
298	nc(ipts, ivmap, iagec) = 0
299	ELSE
300	nc(ipts, ivmap, iagec) = 1
301	!Relative density index (RDI) is needed here in order to
302	!calculate susceptibility to RDI in the next circ_class loop
303	qm_dia = wood_to_qmdia( &
304	circ_class_biomass(ipts, ivm, :, :, icarbon), &
305	circ_class_n(ipts, ivm, :), ivm)
306
307	rdi(ipts, ivm) = (SUM(circ_class_n(ipts, ivm, :)) * m2_to_ha) / &
308	Nmax(qm_dia * m_to_cm, ivm)
309	! Calculate the total area of this species group across all age
310	! classes. Indicate that there bark beetles damage needs to be
311	! calculated in this species group.
312	vegetmax_gp(ipts, ivmap) = vegetmax_gp(ipts, ivmap) + veget_max(ipts, ivm)
313	beetle_gp(ipts, ivmap) = .TRUE.
314	ENDIF
315
316	! Deciduous tree prevent picea trees to be infested by bark beetle
317	! by acting "smoke screen". We include that effect in the code
318	! later through share susceptibility.
319	IF (is_deciduous(ivm) .AND. is_tree(ivm)) THEN
320	vegetmax_deciduous(ipts) = vegetmax_deciduous(ipts) + veget_max(ipts, ivm)
321	ENDIF
322
323	ENDDO ! iagec
324
325	IF (beetle_gp(ipts, ivmap) .AND. vegetmax_gp(ipts, ivmap) .GT. min_stomate) THEN
326	DO iagec = 1, nagec_pft(ivmap)
327	IF (nc(ipts, ivmap, iagec) .EQ. 0) THEN
328	CYCLE
329	ENDIF
330	ivm = start_index(ivmap) + (iagec - 1)
331	! All picea pft class are conrcerned by suceptibility to RDI
332	rdi_gp(ipts, ivmap) = rdi_gp(ipts, ivmap) + rdi(ipts, ivm) * un / SUM(nc(ipts, ivmap, :))
333
334	! Actual biomass of the PFT (gC or gN m-2)
335	biomass_current(ipts, ivmap) = biomass_current(ipts, ivmap) + &
336	SUM(SUM(circ_class_biomass(ipts, ivm, :, :, icarbon), 2) * &
337	circ_class_n(ipts, ivm, :))
338	! Woody biomass of the PFT (gC or gN m-2)
339	woody_current(ipts, ivmap) = woody_current(ipts, ivmap) + &
340	SUM((circ_class_biomass(ipts, ivm, :, isapabove,icarbon) + &
341	circ_class_biomass(ipts, ivm, :, iheartabove, icarbon)) * &
342	circ_class_n(ipts, ivm, :))
343	woody_frac(ipts, ivmap) = SUM((circ_class_biomass(ipts, ivm, :, isapabove,icarbon) + &
344	circ_class_biomass(ipts, ivm, :, iheartabove, icarbon))* &
345	circ_class_n(ipts, ivm, :)) / &
346	SUM(SUM(circ_class_biomass(ipts, ivm, :, :, icarbon), 2) * &
347	circ_class_n(ipts, ivm, :))
348	! Update the average biomass (gC or gN m-2)
349	! of dead wood left on this PFT during previous time steps.
350	! Only biomass
351	! added during the legacy period is accounted for.
352	wood_leftover_mean(ipts, ivmap) = MAX(wood_leftover_mean(ipts, ivmap) + &
353	(SUM(wood_leftover_legacy(ipts, ivm, :)) / legacy_years - &
354	SUM(beetle_damage_legacy(ipts, ivm,:) * &
355	woody_frac(ipts, ivmap))), zero)
356	season_drought_mean(ipts, ivm) = &
357	SUM(season_drought_legacy(ipts, ivm, :)) / legacy_years
358
359	ENDDO
360
361	ENDIF
362
363	ENDDO
364
365	! The main PFT loop :
366	! Bark beetle damage is calculated across all age classes of the
367	! same species group. Hence the loop over nvmap instead of nvm.
368	DO ivmap = 1, nvmap
369
370	! The fist age class PFT use when PFT parameter is needed in
371	! equation outside iagec loop
372	xvm = start_index(ivmap)
373
374	! In the original model "share" represent the purity of the stand
375	! (monospecific vs mixte). This index reflect the fact that beetle are
376	! sensitive to deciduous (AGE_GROUP 6 and 8 for 15 PFT) species
377	! which act as natural "smoke screen".
378	IF(vegetmax_deciduous(ipts) .GT. zero) THEN
379	DO iagec = 1, nagec_pft(ivmap)
380	IF (nc(ipts, ivmap, iagec) .EQ. 0) THEN
381	CYCLE
382	ENDIF
383	ivm = start_index(ivmap) + (iagec - 1)
384	share = vegetmax_deciduous(ipts) / &
385	(vegetmax_gp(ipts, ivmap) + vegetmax_deciduous(ipts))
386	share_suscept_gp(ipts, ivmap) = &
387	un / (un + exp(share_susceptibility_a(ivm) * &
388	(share - share_susceptibility_b(ivm))))
389	ENDDO
390	ELSE
391	share_suscept_gp(ipts, ivmap) = un
392	END IF
393
394	!Debug
395	IF(printlev_loc >= 4)THEN
396	WRITE(numout,*) 'vegetmax_group', vegetmax_gp(ipts, ivmap)
397	WRITE(numout,*) 'rdi_group', rdi_gp(ipts, ivmap)
398	WRITE(numout,*) 'vegetmax_deciduous', vegetmax_deciduous(ipts)
399	WRITE(numout,*) 'beetle_group', beetle_gp(ipts, ivmap)
400	WRITE(numout,*) 'nc', nc(ipts, ivmap, :)
401	END IF
402	!-
403
404	IF (beetle_gp(ipts, ivmap) .AND. vegetmax_gp(ipts, ivmap) .GT. min_stomate) THEN
405	! Bark beetle damage needs to be calculated for this species group and
406	! there is vegetation present.
407	! Epidemic risk is an important concept of this model. It give the
408	! probability that an epidemic phase will happen. An Epidemic state
409	! is define by diferences in the importance of the specific
410	! susceptibility in the calculation of the infested rate.
411	DO iagec = 1, nagec_pft(ivmap)
412	IF (nc(ipts, ivmap, iagec) .EQ. 0) THEN
413	CYCLE
414	ENDIF
415	ivm = start_index(ivmap) + (iagec - 1)
416	! risk_index_legacy are averaged is across legacy_year. It is a minimum
417	! and can be increase up to the whole isk_index_legacy vector
418	epidemic_risk(ipts, ivmap) = epidemic_risk(ipts, ivmap) + &
419	SUM(risk_index_legacy(ipts, ivm, :)) / legacy_years * &
420	veget_max(ipts, ivm) / vegetmax_gp(ipts, ivmap)
421	epidemic_monitor(ipts, ivm) = SUM(risk_index_legacy(ipts, ivm, :)) / legacy_years
422	ENDDO ! iagec
423
424	! The start of an epidemic is only triggered by risk_index over the
425	! last legacy_years years. wght_sirdi_b(xvm) relate is a risk index
426	! value at wich wght_sirdi (the importance the rdi
427	! susceptibility )become dominant over woodleftover susceptibility
428	! this trick ensure a smooth transition from endemic to epidemic
429	! state.
430	IF(epidemic_risk(ipts, ivmap) > wght_sirdi_b(xvm)) THEN
431	! There is a minimum requierement to go out of an epidemic phase
432	! fist : the stand need to reach a RDI thresold meaning that the
433	! there is not enough stems available for the actual beetle
434	! population which will fly away to another location.
435	IF(rdi_gp(ipts,ivmap) .LT. rdi_limit(xvm)) THEN
436	DO iagec = 1, nagec_pft(ivmap)
437	IF (nc(ipts, ivmap, iagec) .EQ. 0) THEN
438	CYCLE
439	ENDIF
440	ivm = start_index(ivmap) + (iagec - 1)
441	! The fist time we go out of the epidemic, a lot of beetle
442	! will flayaway to find new source of food. This reduce the
443	! beetle pressure in the pixel and help the stand to come
444	! back to a recover phase.
445	IF(epidemic(ipts, ivm) .EQ. 1) THEN
446	beetle_flyaway(ipts, ivm) = remaining_beetles(ivm)
447	ENDIF
448	epidemic(ipts, ivm) = 0
449	ENDDO
450	ELSE
451	! There is a bark beetle outbreak going on !
452	DO iagec = 1, nagec_pft(ivmap)
453	IF (nc(ipts, ivmap, iagec) .EQ. 0) THEN
454	CYCLE
455	ENDIF
456	ivm = start_index(ivmap) + (iagec - 1)
457	epidemic(ipts, ivm) = 1
458	ENDDO
459	ENDIF
460	! second : when the risk is low (<wght_sirdi_b), epidemic is
461	! impossible
462	ELSE
463	DO iagec = 1, nagec_pft(ivmap)
464	IF (nc(ipts, ivmap, iagec) .EQ. 0) THEN
465	CYCLE
466	ENDIF
467	ivm = start_index(ivmap) + (iagec - 1)
468	epidemic(ipts, ivm) = 0
469	ENDDO
470	ENDIF
471
472	DO iagec = 1, nagec_pft(ivmap)
473	IF (nc(ipts, ivmap, iagec) .EQ. 0) THEN
474	CYCLE
475	ENDIF
476	ivm = start_index(ivmap) + (iagec - 1)
477
478	! Beetle_flyaway alway come back to 1 (neutral) a rate define by
479	! legacy_years.
480	IF( beetle_flyaway(ipts,ivm) .LT. un ) THEN
481	IF(rdi_gp(ipts,ivmap) .LT. rdi_limit(xvm)) THEN
482	beetle_flyaway(ipts,ivm) = remaining_beetles(ivm)
483	ELSE
484	beetle_flyaway(ipts,ivm) = MIN(beetle_flyaway(ipts,ivm) + &
485	remaining_beetles(ivm)/(legacy_years*6),un)
486	ENDIF
487	ENDIF
488
489	! This is where the importance of each specific susceptibilities
490	! are updated according to epidemic state and risk level. A S shape
491	! relationship with risk level is used to ensure a smooth
492	! transition between epidemic and endemic states.
493	IF( epidemic(ipts, ivm) .EQ. 1 ) THEN
494	!When an epidemic phase is triggered beetle feed on all the
495	!trees available. Beetles don't need woodleftover anymore.
496	wght_siw = zero
497	wght_sirdi = un - (wght_sis(xvm) + wght_sid(xvm) + wght_siw)
498	ELSE
499	wght_sirdi = un / ( un + exp(wght_sirdi_a(xvm) * (epidemic_risk(ipts, ivmap) - &
500	wght_sirdi_b(xvm)))) * (un - (wght_sis(xvm) + wght_sid(xvm)))
501	wght_siw = un - (wght_sirdi + wght_sis(xvm) + wght_sid(xvm))
502	ENDIF
503
504	! Calculate the average bark beetle pressure. The average of three
505	! highest
506	! values from within the total number of legacy_years is used.
507	! bgi (number of beetle generation in one year)
508	bgi_current(:) = beetle_generation_index(ipts, ivm, :)
509
510	! Sort the values to retrieve the three highest values.
511	CALL quicksort(bgi_current, 1, legacy_years)
512
513	! Average the three recent years with the highest bark beetle
514	! pressure. The
515	! same bark beetle pressure is experienced by the whole pixel.
516	! G (number of beetle generation in one year)
517	G(ipts, ivmap) = G(ipts, ivmap) + SUM(bgi_current(legacy_years - &
518	(nb_years_bgi - 1):legacy_years)) / &
519	nb_years_bgi * veget_max(ipts, ivm) / vegetmax_gp(ipts, ivmap)
520
521	IF (.NOT. beetle_pft(ivm) .OR. veget_max(ipts, ivm) .LT. min_stomate) THEN
522	! No reason to calculate bark beetle damage
523	CYCLE
524	ELSE
525	IF(printlev_loc >= 4) THEN
526	WRITE(numout,*) 'bgi_current', bgi_current(:)
527	WRITE(numout,*) 'wood_leftover', wood_leftover_legacy(ipts, ivm, :)
528	WRITE(numout,*) 'beetle_damage_legacy', beetle_damage_legacy(ipts, ivm, :)
529	WRITE(numout,*) 'woody_frac', woody_frac(ipts,ivm)
530	WRITE(numout,*) 'wood_leftover_mean', &
531	SUM(wood_leftover_legacy(ipts, ivm, :)) / legacy_years - &
532	SUM(beetle_damage_legacy(ipts, ivm,:) * &
533	woody_frac(ipts,ivmap))
534	WRITE(numout,*) 'season_drought_mean',season_drought_mean(ipts, ivm)
535	WRITE(numout,*) 'epidemic_risk', epidemic_risk(ipts, ivmap)
536	WRITE(numout,*) 'epidemic flag', epidemic(ipts, ivm)
537	WRITE(numout,*) 'biomass_current', biomass_current(ipts,ivmap)
538	WRITE(numout,*) 'woody_current', woody_current(ipts, ivmap)
539	WRITE(numout,*) 'RDI group', rdi_gp(ipts, ivmap)
540	WRITE(numout,*) 'RDI', rdi(ipts, ivm)
541	WRITE(numout,*) 'wght_sirdi', wght_sirdi
542	WRITE(numout,*) 'wght_siw', wght_siw
543	WRITE(numout,*) 'beetle_flyaway', beetle_flyaway(ipts,ivm)
544	WRITE(numout,*) 'rdi_limit', rdi_limit(ivm)
545	ENDIF
546
547
548	! The self-thinning and yield relationships are for diameters
549	! expressed in cm
550	rdi_suscept_gp(ipts, ivmap) = rdi_suscept_gp(ipts, ivmap) + &
551	(un / (un + exp(rdi_susceptibility_a(ivm) * &
552	(rdi(ipts, ivm) - rdi_susceptibility_b(ivm)))))* &
553	veget_max(ipts, ivm) / vegetmax_gp(ipts, ivmap)
554
555	!!!!!!! NOT USED ANYMORE !!!!!!!!
556	! Stand age effect on bark beetle vulnerability. First calculaed at the
557	! pft level and then accumulated for all age classes (thus accumulated
558	! for the species group)
559	!standage_suscept(ipts, ivm) = age_susceptibility_a(ivm) + &
560	! ( un - age_susceptibility_a(ivm)) / &
561	! ( un + exp(-(age_susceptibility_b(ivm) * &
562	! (age_stand(ipts, ivm) - age_susceptibility_c(ivm)))))
563	!standage_suscept_gp(ipts, ivmap) = standage_suscept_gp(ipts, ivmap) + &
564	! standage_suscept(ipts ,ivm) * &
565	! veget_max(ipts, ivm) / vegetmax_gp(ipts, ivmap)
566
567	! LONG term drought susceptibility is an indirect measurement
568	! tree defence to beetle colonisation. Multiple drought event
569	! will decrease the amount of repelant produce by the tree.
570	drought_suscept(ipts,ivm) = un / (un + exp(drought_susceptibility_a(ivm) * &
571	((un - season_drought_mean(ipts, ivm)) - drought_susceptibility_b(ivm))))
572	drought_suscept_gp(ipts, ivmap) = drought_suscept_gp(ipts, ivmap) + &
573	drought_suscept(ipts, ivm) * veget_max(ipts, ivm) / vegetmax_gp(ipts, ivmap)
574
575	! Susceptibility_index is used to calculate the mortality
576	! rate of infested trees. It is conceptually different from
577	! susceptibility_index_group use to calculate trees that will
578	! be infested by bark beetle.
579	! Mortality rate only consider susceptibilty
580	! which affect physiological state of a tree like long term
581	! drought or intraspecific comptition (RDI). windtrhow and
582	! share are not present in the equation because rely on
583	! mechanism involve only in the infesteation stage.
584	suscept_index(ipts, ivm) = (un / (un + exp(rdi_susceptibility_a(ivm) * &
585	(rdi(ipts, ivm) - rdi_susceptibility_b(ivm))))) * wght_sirdi + &
586	drought_suscept(ipts, ivm) * (un - wght_sirdi)
587	ENDIF
588	ENDDO ! end loop nagec
589
590	! The calculation of the windthrow suceptibility differs from the
591	! original paper of Temperli et al. 2013. First, the wood biomass
592	! lost by windthrow is replace by the wood leftover. This change
593	! implied that windthrow suceptibility is only representing beetle x
594	! windthrow interaction anymore. Of course, this index is still
595	! increasing drasticaly after a windthrow event. Second, the
596	! equation is using the actual woody biomass instead of the maximal
597	! potential woody biomass, a fix value of 300 t/ha in the paper, in
598	! ORCHIDEE we never reach this value and it will
599	! underestimate the contribution of windthrows.
600	! The amount of dead wood left on-site by previous wind storms is
601	! smaller than the living aboveground biomass. This will be accounted
602	! for to decrease the pft-level susceptibility to bark beetle attacks
603	! due to the left-over of windthrow debris. When 30% of
604	! the current wooody biomass is on the floor
605	! susceptibility reach its maximum(1)(unitless)
606	windthrow_suscept_gp(ipts, ivmap) = MIN(un, (wood_leftover_mean(ipts, ivmap) / &
607	woody_current(ipts, ivmap)) / windthrow_susceptibility_tune(start_index(ivmap)))
608
609	! Calculate pixel sucseptibility at the level of the whole species group. It is
610	! assumed that half of the vulnerability comes from the relative density index.
611	! But here RDI is not use has a proxy for competition but rather
612	! than a proxy for bettle spead rate wich is stimulatate when tree
613	! canopy close from each other. Long term drought (proxy of plant health)
614	! and pixel deciduous share also account in the calculation.
615	! wood leftover is also use to trigger the epidemic phase by increasing
616	! the risk index right after a storm but don't have any impact
617	! during the epidemic phase.
618	! susceptibility_index_group is multiplied with vegetmax to account for area in the
619	! pixel where susceptibility_index is zero because there is no suceptible forest
620	! i.e. picea abies. Finally susceptibility_index_group will be used
621	! to calculate the amount of picea tree that will be infested by
622	! bark beetle.
623
624	suscept_index_gp(ipts, ivmap) = &
625	(drought_suscept_gp(ipts, ivmap) * wght_sid(ivm) + &
626	rdi_suscept_gp(ipts, ivmap) * wght_sirdi + &
627	share_suscept_gp(ipts, ivmap) * wght_sis(ivm) + &
628	windthrow_suscept_gp(ipts, ivmap) * wght_siw) * &
629	vegetmax_gp(ipts, ivmap)
630
631	DO iagec = 1, nagec_pft(ivmap)
632	IF (nc(ipts, ivmap, iagec) .EQ. 0) THEN
633	CYCLE
634	ENDIF
635	ivm = start_index(ivmap) + (iagec - 1)
636	! Estimate the risk index at the Age class level.
637	suscept_index_monitor(ipts, ivm) = suscept_index_gp(ipts, ivmap)
638	windthrow_suscept_monitor(ipts, ivm) = windthrow_suscept_gp(ipts,ivmap)
639
640	END DO
641
642
643	!Debug
644	IF(printlev_loc >= 4) THEN
645	WRITE(numout,*) 'drought_susceptibility', drought_suscept_gp(ipts, ivmap)
646	WRITE(numout,*) 'rdi_supceptibilit', rdi_suscept_gp(ipts, ivmap)
647	WRITE(numout,*) 'share_susceptibility', share_suscept_gp(ipts, ivmap)
648	WRITE(numout,*) 'windthrow_susceptibility', windthrow_suscept_gp(ipts, ivmap)
649	END IF
650	!-
651
652	! Calculate as the real pressure of the beetles. Avegare of the average bark
653	! beetle of the past three years and the beetle legacy (which is a ratio between
654	! affected dead wood (breeding grounds) and fresh biomass (food)). This average
655	! is an index for the available beetles when multiplied with the susceptibility
656	! of the species group we get the beetle pressure (which could be seen as a proxy
657	! for beetle biomass). Make sure the index stays between zero and one. Here the
658	! lower bound was set to 0,05 to reflect the fact that there is always an endemic
659	! population of bark beetles. If we would allow this index to go to zero beetles
660	! could go extint. Given that both G and beetle_legacy_index can reach zero it is
661	! slightly easier to control the value of the beetle_pressure_index here.
662	beetle_pressure_index(ipts, ivmap) = MIN(MAX(suscept_index_gp(ipts, ivmap) * &
663	(G(ipts, ivmap) + beetle_pop_legacy(ipts, xvm, 1)) * 0.5* &
664	beetle_flyaway(ipts, ivmap), 0.05), un)
665
666	! Estimate the risk index at the pixel level.
667	risk_index_gp(ipts, ivmap) = suscept_index_gp(ipts, ivmap) * &
668	beetle_pressure_index(ipts, ivmap)
669
670	! This loop is use for monitoring intermediate result.
671	! ivmap ---> ivm, because there is no ivmap dimension in the
672	! output files
673	DO iagec = 1, nagec_pft(ivmap)
674	IF (nc(ipts, ivmap, iagec) .EQ. 0) THEN
675	CYCLE
676	ENDIF
677	ivm = start_index(ivmap) + (iagec - 1)
678	risk_index(ipts, ivm) = risk_index_gp(ipts, ivmap)
679	suscept_index_monitor(ipts, ivm) = suscept_index_gp(ipts, ivmap)
680	beetle_pressure_monitor(ipts, ivm) = beetle_pressure_index(ipts, ivmap)
681	windthrow_suscept_monitor(ipts, ivm) = windthrow_suscept_gp(ipts, ivmap)
682	END DO
683
684	!Debug
685	IF(printlev_loc >= 4) THEN
686	WRITE(numout,*) 'G', G(ipts, ivmap)
687	WRITE(numout,*) 'beetle_pressure_index', beetle_pressure_index(ipts, ivmap)
688	WRITE(numout,*) 'beetle_pop_legacy', beetle_pop_legacy(ipts, xvm, 1)
689	WRITE(numout,*) 'risk_index_group', risk_index_gp(ipts, :)
690	WRITE(numout,*) 'susceptibility_index_group', suscept_index_gp(ipts, ivmap)
691	END IF
692	!-
693
694	! Estimate the total biomass of spruce infected by bark beetle. The
695	! calculation diverge from the original equation which used a fix
696	! value corresponding the half of the maximun potential biomass (also fix).
697	! In our case maximun potential biomass is not a fix value and depend
698	! climate forcing, soil properties, interspecies water competition. Seems
699	! easier to use the actual biomass and recalibrate suceptibility parameters
700	! in order to decrease the overall risk index values.
701	biomass_infected(ipts, ivmap) = biomass_current(ipts, ivmap) * risk_index_gp(ipts, ivmap)
702
703	! Calculate Beetle pressure index. The actual biomass(biomass_current) is
704	! used instead of the maximal potential biomass, i.e., 300 t/ha, as was done
705	! in temperli et al. 2013. This approach differs from Temperli. The ratio
706	! between dead wood from previous beetle outbreaks (breeding grounds) and
707	! uninfected biomass (fresh food) is seen as a driver of new infections.
708	! Make sure the legacy stays inbetween zero and one.
709	beetle_temp1(ipts, ivm, :) = zero
710
711	DO iyear = 1, legacy_years - 1
712	beetle_temp1(ipts, ivm, iyear + 1) = beetle_pop_legacy(ipts, ivm, iyear)
713	ENDDO
714
715	beetle_temp1(ipts, ivm, 1) = MIN(biomass_infected(ipts, ivmap) / &
716	(biomass_current(ipts, ivmap) * pressure_feedback(ivm)), un)
717	beetle_pop_legacy(ipts, ivm, :) = beetle_temp1(ipts, ivm, :)
718
719	! Distribute the total infected biomass over the different age classes
720	! of the same species
721	DO iagec = 1, nagec_pft(ivmap)
722	IF (nc(ipts, ivmap, iagec) .EQ. 0) THEN
723	CYCLE
724	ENDIF
725	! Determine the PFT index for the age classes within the species
726	! group under consideration
727	ivm = start_index(ivmap) + (iagec - 1)
728
729	! Estimate mortality rate at the PFT level (unitless).
730	mortality_rate(ipts, ivm) = (suscept_index(ipts, ivm) + &
731	beetle_pressure_index(ipts, ivmap)) * 0.5
732
733	! Estimate living biomass killed by bark beetles (gC or gN m-2)
734	beetle_damage(ipts, ivm) = mortality_rate(ipts, ivm) * &
735	biomass_infected(ipts, ivmap) * &
736	veget_max(ipts, ivm) / &
737	vegetmax_gp(ipts, ivmap)
738
739	beetle_temp2(ipts, ivm, :) = zero
740	DO iyear = 1, beetle_legacy - 1
741	beetle_temp2(ipts, ivm, iyear + 1) = &
742	beetle_damage_legacy(ipts, ivm, iyear)
743	ENDDO
744	beetle_temp2(ipts, ivm, 1) = beetle_damage(ipts, ivm)
745	beetle_damage_legacy(ipts, ivm, :) = beetle_temp2(ipts, ivm, :)
746
747	!Debug
748	IF(printlev_loc >= 4) THEN
749	WRITE(numout,*) 'biomass_current', biomass_current(ipts,ivmap)
750	WRITE(numout,*) 'biomass_infected', biomass_infected(ipts, ivmap)
751	WRITE(numout,*) 'beetle_damage', beetle_damage(ipts, ivm)
752	WRITE(numout,*) 'effective mortality in %', beetle_damage(ipts, ivm) / &
753	biomass_current(ipts,ivmap) * 100
754	WRITE(numout,*) 'mortality_rate', mortality_rate(ipts, ivm)
755	WRITE(numout,*) 'vegetmax_group', vegetmax_gp(ipts, ivmap)
756	WRITE(numout,*) 'veget_max', veget_max(ipts, ivm)
757	END IF
758	!-
759	! Select which circ_class would be damaged based on the thier
760	! diameters
761	ncirc_loc = 0
762	qdia(:) = wood_to_dia(circ_class_biomass(ipts, ivm, :, :, icarbon), ivm)
763	DO icir = 1, ncirc
764	IF (qdia(icir) > 0.075) THEN
765	ncirc_loc = ncirc_loc + 1
766	WRITE(numout,) 'diameter in cm:', qdia(icir) 100
767	ENDIF
768	ENDDO
769	WRITE(numout,*) 'ncirc_loc: ', ncirc_loc
770	! By using an decreasing index, I made the assumption that the
771	! first
772	! tree which die from bark beetle infection are the bigest one.
773	DO icir = ncirc_loc, 1, -1
774	IF(SUM(circ_class_biomass(ipts, ivm, icir, :, icarbon)) * &
775	circ_class_n(ipts, ivm, icir) .LT. beetle_damage(ipts, ivm)) THEN
776	! All individuals in the circumference class need to be killed to satisfy
777	! the total amount of biomass that needs to be killed.
778	circ_class_kill(ipts, ivm, icir, forest_managed(ipts, ivm), icut_beetle) = &
779	circ_class_n(ipts, ivm, icir)
780	beetle_damage(ipts, ivm) = beetle_damage(ipts, ivm) - &
781	SUM(circ_class_biomass(ipts,ivm, icir, :, icarbon)) * circ_class_n(ipts, ivm, icir)
782	ELSE
783	IF(beetle_damage(ipts, ivm) .GT. zero) THEN
784	! Only part of the trees in this circumference class should be killed
785	circ_class_kill(ipts, ivm, icir, forest_managed(ipts, ivm), icut_beetle) = &
786	circ_class_n(ipts, ivm, icir) * beetle_damage(ipts, ivm) / &
787	(SUM(circ_class_biomass(ipts, ivm, icir, :, icarbon)) * circ_class_n(ipts, ivm, icir))
788	beetle_damage(ipts, ivm) = zero
789	ELSE
790	! The total demand of mortality has been satisfied. No more killing is
791	! required.
792	circ_class_kill(ipts, ivm, icir, forest_managed(ipts, ivm), icut_beetle) = zero
793	ENDIF
794	ENDIF
795	ENDDO ! End loop ncirc
796
797	ENDDO ! End loop nagec
798
799	ELSE
800	! No mortality from bark beetles. Set the mortality to zero to avoid problems
801	! with uninitialized values
802	DO iagec = 1, nagec_pft(ivmap)
803	IF (nc(ipts, ivmap, iagec) .EQ. 0) THEN
804	CYCLE
805	ENDIF
806	! Determine the PFT index for the age classes within the species
807	! group under consideration
808	ivm = start_index(ivmap) + (iagec - 1)
809	! No trees are to be killed
810	circ_class_kill(ipts, ivm, :, forest_managed(ipts, ivm), icut_beetle) = zero
811	ENDDO ! nagec
812	ENDIF ! beetle_group .EQV. .TRUE. .AND. vegetmax_group(ipts,ivmap) .GT. min_stomate
813
814	ENDDO ! End loop nvmap
815	ENDDO ! End loop npts
816
817	END SUBROUTINE bark_beetle_damage
818
819
820	! +++MOVE TO FUNCTION_LIBRARY.F90
821	! +++COMPLETE HEADER
822	!=============================================================================================================================
823	!! SUBROUTINE : quicksort
824	!!
825	!>\BRIEF
826	!!
827	!!\n DESCRIPTION :
828	!!
829	!! RECENT CHANGE(S):
830	!!
831	!! MAIN OUTPUT VARIABLE(S) :
832	!!
833	!! REFERENCES:
834	!!
835	!! FLOWCHART :
836	!! \n
837	!_
838	!=============================================================================================================================
839
840	! quicksort.f --f90--
841	! Author: t-nissie
842	! License: GPLv3
843	! Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
844	RECURSIVE SUBROUTINE quicksort(a, first, last)
845	IMPLICIT NONE
846	REAL8 a(), x, t
847	INTEGER first, last
848	INTEGER i, j
849	!_
850	!=============================================================================================================================
851	x = a( (first + last) * 0.5 )
852	i = first
853	j = last
854	DO
855	DO WHILE (a(i) < x)
856	i = i + 1
857	ENDDO
858	DO WHILE (x < a(j))
859	j = j - 1
860	ENDDO
861	IF (i >= j) EXIT
862	t = a(i); a(i) = a(j); a(j) = t
863	i = i + 1
864	j = j - 1
865	ENDDO
866	IF (first < i - 1) call quicksort(a, first, i - 1)
867	IF (j + 1 < last) call quicksort(a, j + 1, last)
868	END SUBROUTINE quicksort
869
870	END MODULE stomate_pest
871

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

Download in other formats: