Context Navigation

source: branches/ORCHIDEE_Quest/ORCHIDEE/src_stomate/lpj_fire.f90 @ 7406

Last change on this file since 7406 was 4693, checked in by josefine.ghattas, 7 years ago
Clean in text output, see ticket #394
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 28.8 KB

Line
1	! =================================================================================================================================
2	! MODULE : lpj_fire
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006). This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
7	!
8	!>\BRIEF Calculates the probability of fire and its effects on the carbon cycle
9	!!
10	!!\n DESCRIPTION: None
11	!!
12	!! RECENT CHANGE(S): None
13	!!
14	!! SVN :
15	!! $HeadURL$
16	!! $Date$
17	!! $Revision$
18	!! \n
19	!_ ================================================================================================================================
20
21	MODULE lpj_fire
22
23	! modules used:
24	USE xios_orchidee
25	USE stomate_data
26	USE ioipsl_para
27	USE pft_parameters
28	USE constantes
29
30	IMPLICIT NONE
31
32	! private & public routines
33
34	PRIVATE
35	PUBLIC fire,fire_clear
36
37	LOGICAL, SAVE :: firstcall_fire = .TRUE. !! first call
38	!$OMP THREADPRIVATE(firstcall_fire)
39
40	CONTAINS
41
42
43	!! ================================================================================================================================
44	!! SUBROUTINE : fire_clear
45	!!
46	!>\BRIEF Set the firstcall_fire flag to .TRUE. and activate initialization
47	!!
48	!_ ================================================================================================================================
49
50	SUBROUTINE fire_clear
51	firstcall_fire = .TRUE.
52	END SUBROUTINE fire_clear
53
54
55	!! ================================================================================================================================
56	!! SUBROUTINE : fire
57	!!
58	!>\BRIEF Calculate fire index and fraction of area burned by comparing
59	!! daily litter moisture with prescribed moisture of fire inhibition. If daily
60	!! moisture is below the precribed threshold, allow fire disturbance and
61	!! calculate CO_2 emissions from fire. The main algorithm follows Thonicke et al.
62	!! (2001)
63	!!
64	!! DESCRIPTION : Fire occurs when the three basic prerequisites are met:
65	!! sufficient fuel on the ground, relatively dry fuel (fuel moisture lower than
66	!! the moisture of extinction), and presence of ignition sources. The flag
67	!! ::disable_fire is used to use this fire module or not. While here in the module the
68	!! ignition source is not explicitely represented. It's assumed
69	!! that the ignition source is always available, i.e., if a sufficient amount
70	!! of dead fuel exists with a moisture content below the moisture of fire
71	!! inhibition, then both live and dead fuel will start to burn.\n
72	!!
73	!! The module completes the following tasks:
74	!! 1. Calculates daily fire index, long term fire index and available ground
75	!! litter for burning.\n
76	!! \latexonly
77	!! \input{lpj_fire_1.tex}
78	!! \endlatexonly
79	!! Where, m is the daily litter moisture, m_e is the moisture of extinction, p(m)
80	!! is probability of fire (i.e. the daily fire index).\n
81	!!
82	!! \latexonly
83	!! \input{lpj_fire_2.tex}
84	!! \endlatexonly
85	!! Where, s is the long term fire index and A(s) is the annual fire fraction.\n
86	!! 2. Calculates annual fire fraction, then transform to daily fraction in the
87	!! same time step of stomate.\n
88	!! 3. Ground litter and grass live biomass (in growing season and except root and carbon
89	!! reserve pool) are burned at full fire fraction.\n
90	!! 4. Fire fraction for tree PFT are compromised by prescribed fire resistence.
91	!! Tree live biomass are consumed using this compromised fire fraction and consumption
92	!! fraction for different biomass parts. In the case of activation of dynamic
93	!! vegetation, tree individual density is updated after fire. \n
94	!! 5. For all types of fuel (either ground litter or live biomass) that are
95	!! burned, there is a certain fraction of âresidueâ that are not completely
96	!! burned. The remaining part is transferred to (in case of biomass) or
97	!! remains (in case of ground litter) as litter.\n
98	!! 6. Update the biomass, litter or dead leaves pool after burning.
99	!! 7. If the flag SPINUP_ANALYTIC is set to true, the matrix A is updated following
100	!! Lardy (2011).
101	!!
102	!! RECENT CHANGE(S): April 2015: Black carbon calculations is removed. The black carbon was not conserved.
103	!!
104	!! MAIN OUTPUT VARIABLE(S): ::co2_fire or the carbon emitted into the atmosphere
105	!! by fire, including both living and dead biomass (gC m^{-2} dtslow^{-1})$ @endtex
106	!!
107	!! REFERENCES :
108	!! - Thonicke K., Venevsky S., Sitch S., and Cramer W. (2001) The role of fire
109	!! disturbance for global vegetation dynamics: coupling fire into a Dynamic
110	!! Global Vegetation Model, Global Ecology & Biogeography, 10, 661-677.
111	!! - Kuhlbusch et al. JGR 101, 23651-23665, 1996
112	!! - Kuhlbusch & Crutzen, GBC 9, 491-501, 1995
113	!! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium,
114	!! Environmental Modelling & Software (2011), doi:10.1016\|j.envsoft.2011.05.016
115	!!
116	!! FLOWCHART : None
117	!! \n
118	!_ ================================================================================================================================
119
120	SUBROUTINE fire (npts, dt, litterpart, &
121	litterhum_daily, t2m_daily, lignin_struc,veget_cov_max, &
122	fireindex, firelitter, biomass, ind, &
123	litter, dead_leaves, bm_to_litter, &
124	co2_fire, MatrixA)
125
126	!! 0. Variable and parameter declarations
127
128	!! 0.1 Input variables
129
130	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of pixels
131	!! (unitless)
132	REAL(r_std), INTENT(in) :: dt !! Time step of the simulations for stomate
133	!! (days)
134	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: litterpart !! [DISPENSABLE] Fraction of litter above
135	!! the ground belonging to different PFTs
136	!! ?? this variable is used but might be
137	!! redundant (with value of always 1) ??
138	!! To check but probably litterpart to be
139	!! removed. Probably a residual from
140	!! nat/agri before merge
141	REAL(r_std), DIMENSION(npts), INTENT(in) :: litterhum_daily !! Daily litter moisture (unitless, 0-1)
142	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_daily !! Daily 2 meter temperature (K)
143	REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(in) :: lignin_struc !! Ratio Lignin/Carbon in structural above
144	!! and below ground litter
145	!! @tex $(gC gC^{-1})$ @endtex
146	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_cov_max !! Maximum fraction of vegetation type including
147	!! non-biological fraction (0-1, unitless)
148
149	!! 0.2 Output variables
150
151	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: co2_fire !! Carbon emitted into the atmosphere by
152	!! fire, including both living and dead
153	!! biomass
154	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
155	!! 0.3 Modified variables
156
157	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: fireindex !! Probability of fire; in the code means
158	!! long term fire index (unitless, 0-1)
159	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: firelitter !! Total natural litter
160	!! available for burning above the ground
161	!! @tex $(gC m^{-2})$ @endtex
162	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass !! Biomass @tex $(gC m^{-2})$ @endtex
163	REAl(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind !! Density of individuals
164	!! @tex $(m^{-2})$ @endtex
165	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) :: litter !! Metabolic and structural litter, above
166	!! and below ground
167	!! @tex $(gC m^{-2})$ @endtex
168	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: dead_leaves !! Dead leaves on ground, per PFT,
169	!! metabolic and structural
170	!! @tex $(gC m^{-2})$ @endtex
171	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout):: bm_to_litter !! Biomass entering the litter pool
172	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
173	REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(inout) :: MatrixA !! Matrix containing the fluxes between the carbon pools
174	!! here, it is called only if disable_fire = n
175	!! once a day
176	!! @tex $(gC.m^2.day^{-1})$ @endtex
177
178
179	!! 0.4 Local variables
180
181	REAL(r_std), DIMENSION(npts) :: fire_disturb !! Actual fire disturbance fraction
182	!! after consideration of inherent fire
183	!! resistance of different PFTs
184	!! (unitless, 0-1)
185	REAL(r_std), DIMENSION(npts,nvm) :: firedeath !! In case of activation of dynamic
186	!! vegetation, the daily fraction of
187	!! burned individuals (unitless, 0-1)
188	REAL(r_std), DIMENSION(npts) :: moistlimit !! Moisture for fire inhibition per PFT !!
189	!! (unitless, 0-1) ; temporary variable
190	!! for each PFT in the loop over #PFT ?
191	!! it's not fire inhibition, it's
192	!! "moisture of extinction" by the original
193	!! reference paper, and a more "fire
194	!! professional name" it's better to change
195	!! the name of moistlimit
196	REAL(r_std), DIMENSION(npts) :: litter_above !! Total aboveground litter per PFT !!
197	!! @tex $(gC m^{-2})$ @endtex
198	REAL(r_std), DIMENSION(npts,nvm) :: fireindex_daily !! Daily fire index (unitless, 0-1)
199	REAL(r_std), DIMENSION(npts, nvm) :: firefrac !! Daily fire burning fraction on ground
200	!! (unitless, 0-1)
201	REAL(r_std), DIMENSION(npts) :: struc_residual !! Fraction of structural litter not burned
202	!! (thus residual), depending on strutural
203	!! litter lignin content (unitless, 0-1)
204	REAL(r_std), DIMENSION(npts) :: residue !! Fuel (either biomass or litter) not
205	!! burned @tex $(gC m^{-2})$ @endtex
206	REAL(r_std), DIMENSION(npts) :: x !! Intermediate variable
207	REAL(r_std), DIMENSION(npts) :: aff !! Annual fire fraction (unitless, 0-1)
208	INTEGER(i_std) :: j,k,m !! Indeces
209	!_ ================================================================================================================================
210
211	IF (printlev>=3) WRITE(numout,*) 'Entering fire'
212
213	!! 1. Initialization
214
215	IF ( firstcall_fire ) THEN
216
217	!! 1.1 Fraction to CO_2
218	! What fraction of the plant biomass compartment, if burned, is transformed
219	! into CO_2 released to the atmosphere?
220	!??! the fraction for heartabove seems too big, it's not clear if this value is correct.
221
222	!! 1.2 print control messages
223	IF (printlev >= 2) THEN
224	WRITE(numout,*) 'fire:'
225	WRITE(numout,*) ' > temporal inertia of fire probability (d): ',tau_fire
226
227	WRITE(numout,*) ' > fraction of burned biomass that becomes CO2:'
228	WRITE(numout,*) ' leaves: ', co2frac(ileaf)
229	WRITE(numout,*) ' sap above ground: ', co2frac(isapabove)
230	WRITE(numout,*) ' sap below ground: ', co2frac(isapbelow)
231	WRITE(numout,*) ' heartwood above ground: ', co2frac(iheartabove)
232	WRITE(numout,*) ' heartwood below ground: ', co2frac(iheartbelow)
233	WRITE(numout,*) ' roots: ', co2frac(iroot)
234	WRITE(numout,*) ' fruits: ', co2frac(ifruit)
235	WRITE(numout,*) ' carbohydrate reserve: ', co2frac(icarbres)
236
237	WRITE(numout,) ' > critical litter quantity (gC/m*2): ',litter_crit
238	WRITE(numout,*) ' > We calculate a fire probability on agricultural ground, but'
239	WRITE(numout,*) ' the effective fire fraction is zero.'
240	END IF
241	firstcall_fire = .FALSE.
242
243	ENDIF
244
245	!! 1.4 Initialize output
246	co2_fire(:,:) = zero
247	firedeath(:,:) = zero
248	fireindex_daily(:,:) = zero
249	firefrac(:,:) = zero
250
251	!! 2. Determine fire probability
252
253	! Calculate probability that crops are not burned for the moment
254	! Calculate long-term aboveground litter
255
256	! a long loop over PFT
257	DO j = 2,nvm !loop over #PFT
258
259	!! 2.1 Total above ground litter
260	! Total litter above the ground for the PFT under consideration
261	litter_above(:) = litter(:,imetabolic,j,iabove,icarbon) + &
262	litter(:,istructural,j,iabove,icarbon)
263
264
265	!! 2.2 Calculate the PFT litter amount weighted by moisture
266	! If litter moisture is higher than moisture of extinction, fire is not possible.
267	moistlimit(:) = zero
268
269	WHERE ( (t2m_daily(:) .GT. ZeroCelsius) .AND. (litter_above(:) .GT. min_stomate) )
270	!?? the calculation here is redundant? as the part before flam(j) is 1?
271	!?? see first comment.. litterpart to be removed ?
272	moistlimit(:) = &
273	( litterpart(:,j,imetabolic)*litter(:,imetabolic,j,iabove,icarbon) + &
274	litterpart(:,j,istructural)*litter(:,istructural,j,iabove,icarbon) ) / &
275	litter_above(:) * flam(j)
276
277	ENDWHERE
278
279	!! 2.3 Calculate daily fire index
280	! Calculate daily fire index
281	! \latexonly
282	! \input{lpj_fire_1.tex}
283	! \endlatexonly
284	! Where, m is the daily litter moisture, m_e is the moisture of extinction, p(m)
285	! is probability of fire (i.e. the daily fire index).\n
286	WHERE ( moistlimit(:) .GT. min_stomate )
287	x(:) = litterhum_daily(:)/moistlimit(:)
288	fireindex_daily(:,j) = EXP( - pi * x(:) * x(:) )
289	ELSEWHERE
290	fireindex_daily(:,j) = zero
291	ENDWHERE
292
293	!! 2.4 Calculate long-term fire index
294	! Calculate long-term fire index which is the mean probability of fire
295	fireindex(:,j) = ((tau_fire - dt) * fireindex(:,j) + (dt) * fireindex_daily(:,j)) / tau_fire
296
297	!! 2.5 Calculate long term aboveground litter that are available for burning
298	firelitter(:,j) = &
299	( ( tau_fire-dt ) * firelitter(:,j) + dt * litter_above(:) ) / tau_fire
300
301	!! 3. Calculate fire fraction from litter and fireindex
302
303	!! 3.1 Fire fraction from long term fire index for natural PFTs
304	! Transform the annual fire fraction to daily fraction. This is done by assuming that
305	! each day the fire fraction is the same.
306	aff(:) = firefrac_func (npts, fireindex(:,j))
307
308	! annual_fire_fraction = 1. - ( 1. - daily_fire_fraction )**365
309	! Thus, daily_fire_fraction = 1. - ( 1. - annual_fire_fraction )**(1/365)
310	! If annual firefrac<<1, then firefrac_daily = firefrac * dt/one_year
311	! This approximation avoids numerical problems.
312	! the variable 'un' is 1
313	IF(.NOT.disable_fire.AND.natural(j))THEN
314	WHERE ( aff(:) .GT. 0.1 )
315	firefrac(:,j) = un - ( un - aff(:) ) ** (dt/one_year)
316	ELSEWHERE
317	firefrac(:,j) = aff(:) * dt/one_year
318	ENDWHERE
319	ELSE
320	firefrac(:,j) = zero
321	ENDIF
322
323	! No fire if litter is below critical value
324	WHERE ( firelitter(:,j) .LT. litter_crit )
325	firefrac(:,j) = zero
326	ENDWHERE
327
328	! However, there is a minimum fire extent
329	firefrac(:,j) = MAX( 0.001_r_std * dt/one_year, firefrac(:,j) )
330
331	! if FIRE_DISABLE flag is set no fire
332	IF (disable_fire) firefrac = zero
333
334	!! 3.2 For agricultural ground, no fire is burned for the moment
335	IF ( .NOT. natural(j)) firefrac(:,j) = zero
336
337	!! 4. Determine fire impact
338
339	! Calculate fire disturbance fraction for each PFT, and fire emissions due
340	! to grasses.
341
342	!! 4.1 Tree and grass live biomass
343	! Tree live biomass is not burned in fire.
344	! However, in the dynamic vegetation module, tree individual density
345	! will be updated after fire. The fraction of tree individuals that are
346	! supposed to die in fire is the fire fraction multiplied by the tree PFT fire
347	! resistance which reflect survivorship of different tree PFTs during fire.
348	IF ( is_tree(j) ) THEN
349
350	!! 4.1.1 Disturban,ce factor for trees
351	! Trees are disturbed over the whole year. A resistance factor is
352	! used to reflect survivorship of different tree PFTs during fire.
353	fire_disturb(:) = ( un - resist(j) ) * firefrac(:,j)
354
355	ELSE
356
357	!! 4.1.2 Disturbance factor for grasses
358	! Grass is not disturbed outside the growing season thus grass biomass
359	! is only burned during the growing season.
360	WHERE ( biomass(:,j,ileaf,icarbon) .GT. min_stomate )
361
362	fire_disturb(:) = ( un - resist(j) ) * firefrac(:,j)
363
364	ELSEWHERE
365
366	fire_disturb(:) = zero
367
368	ENDWHERE
369
370	ENDIF
371
372	!! 4.2 Burn live biomass
373	! The burned part goes to either the atmposhere or litter
374	DO k = 1, nparts
375
376	! for tree PFT, all biomass compartments are burned.
377	! for grass biomass compartments, all are burned except root and carbon reserve
378	! IF concerning PFT is tree; OR (the PFT is grass, but not root or carbon reserve biomass); then it's burned.
379	IF ( .NOT. ( ( .NOT. is_tree(j) ) .AND. ( ( k.EQ.iroot ) .OR. ( k.EQ.icarbres ) ) ) ) THEN
380
381	!! 4.2.1 Fraction to the atmosphere.
382	co2_fire(:,j) = co2_fire(:,j)+ biomass(:,j,k,icarbon) * fire_disturb(:) * co2frac(k) / dt
383
384	!! 4.2.2 Determine the residue
385	! Residual is expressed in gC/m^2
386	residue(:) = biomass(:,j,k,icarbon) * fire_disturb(:) * ( un - co2frac(k) )
387
388	!! 4.2.2.3 The rest (largest part) of the residue becomes litter.
389	bm_to_litter(:,j,k,icarbon) = bm_to_litter(:,j,k,icarbon) + residue(:)
390	ENDIF
391
392	ENDDO
393
394
395	!! 4.3 Update biomass pool after burning
396
397	!! 4.3.1 Decrease biomass amount except for grass roots and carbon reserve.
398	DO m = 1, nelements ! Loop over # elements
399
400	DO k = 1, nparts
401
402	! **2
403	IF ( .NOT. ( ( .NOT. is_tree(j) ) .AND. ( ( k.EQ.iroot ) .OR. ( k.EQ.icarbres) ) ) ) THEN
404
405	biomass(:,j,k,m) = ( un - fire_disturb(:) ) * biomass(:,j,k,m)
406
407	ENDIF
408
409	ENDDO
410
411	END DO ! End Loop over # elements
412
413
414	!! 4.3.2 If vegetation is dynamic, then decrease the density of tree individuals.
415	IF ( (ok_dgvm .OR. .NOT.lpj_gap_const_mort) .AND. is_tree(j) ) THEN
416
417	firedeath(:,j) = fire_disturb(:) / dt
418
419	ind(:,j) = ( un - fire_disturb(:) ) * ind(:,j)
420
421	ENDIF
422
423	ENDDO ! loop over #PFT
424
425	!! 5. Burn litter
426
427	! All litter (either forest or grass) is burned throughout the whole year
428	! Metabolic litter is totally burned. Burning fraction of structural litter is related
429	! with its lignin content. The burned part either goes to the atmosphere
430	! or remains in litter as unburned residue.
431	DO j = 2,nvm !loop over #PFT
432
433	!! 5.1 Burn exposed metabolic litter
434	! Exposed metabolic litter burns completely and goes directly into the atmosphere as CO2.
435
436	!! 5.1.1 CO2 flux from litter burning
437	! Flux expressed in gC m^{-2}dtslow^{-1}
438	co2_fire(:,j) = co2_fire(:,j) + litter(:,imetabolic,j,iabove,icarbon) * &
439	firefrac(:,j) / dt
440
441	!! 5.1.2 Decrease metabolic litter
442	DO m = 1,nelements
443	litter(:,imetabolic,j,iabove,m) = litter(:,imetabolic,j,iabove,m) * &
444	( un - firefrac(:,j) )
445	END DO
446
447	!! 5.2 Burning exposed structural litter
448	! The fraction that goes to the atmosphere is related to its lignin content. The remaining unburned
449	! residues remain as litter.
450
451	!! 5.2.1 Incomplete burning of exposed structural litter
452	! Fraction of structural litter that does not burn completly. This fraction depends on lignin
453	! content (lignin_struc).
454	struc_residual(:) = fire_resist_struct * lignin_struc(:,j,iabove)
455
456	!! 5.2.2 CO2 flux from litter burning
457	! Flux expressed in gC m^{-2}dtslow^{-1}
458	co2_fire(:,j) = co2_fire(:,j) + &
459	litter(:,istructural,j,iabove,icarbon) * firefrac(:,j) * &
460	( un - struc_residual(:) )/ dt
461
462	!! 5.2.3 Determine residue
463	! The residue is litter that undergoes fire, but is not transformed into CO2
464	residue(:) = litter(:,istructural,j,iabove,icarbon) * firefrac(:,j) * &
465	struc_residual(:)
466
467	!! 5.2.6 TResidue remaining litter
468	! The rest (largest part) of the residue remains litter. Remaining
469	! litter is the sum of this and of the litter which has not undergone a fire.
470	litter(:,istructural,j,iabove,icarbon) = &
471	litter(:,istructural,j,iabove,icarbon) * ( un - firefrac(:,j) ) + &
472	residue(:)
473
474	!! 5.2.7 Update matrix A for analytical spin-up (only if SPINUP_ANALYTIC is activated)
475
476	IF (spinup_analytic) THEN
477
478	! litter structural above
479	MatrixA(:,j,istructural_above,istructural_above) = MatrixA(:,j,istructural_above,istructural_above) - firefrac(:,j) &
480	+ firefrac(:,j)* struc_residual(:)
481
482	! litter_metabolic above
483	MatrixA(:,j,imetabolic_above,imetabolic_above) = MatrixA(:,j,imetabolic_above,imetabolic_above) - firefrac(:,j)
484
485	!DS! This is the exact formulation. The above is a simplified formulation.
486	! MatrixA(:,j,istructural_above,istructural_above) = MatrixA(:,j,istructural_above,istructural_above)*(1. - firefrac(:,j) + &
487	! firefrac(:,j)* struc_residual(:)
488	! MatrixA(:,j,imetabolic_above,imetabolic_above) = MatrixA(:,j,imetabolic_above,imetabolic_above)*(1 - firefrac(:,j))
489
490	ENDIF !(spinup_analytic)
491
492	ENDDO ! loop over #PFT
493
494	!! 5.3 Update exposed dead leaves (leaf litter) on the ground
495	! Dead leaves are supposed to burn completely in fire even their structural part.
496	DO j = 2,nvm
497
498	DO k = 1, nlitt
499	dead_leaves(:,j,k) = dead_leaves(:,j,k) * ( un - firefrac(:,j) )
500	ENDDO
501
502	ENDDO
503
504	!! 6. write out history variables
505
506	! output in 1/dtslow
507	firefrac(:,:) = firefrac(:,:) / dt
508
509	CALL xios_orchidee_send_field("firefrac",firefrac(:,:))
510	CALL xios_orchidee_send_field("firedeath",firedeath(:,:))
511
512	CALL histwrite_p (hist_id_stomate, 'FIREFRAC', itime, &
513	firefrac(:,:), npts*nvm, horipft_index)
514	CALL histwrite_p (hist_id_stomate, 'FIREDEATH', itime, &
515	firedeath(:,:), npts*nvm, horipft_index)
516
517	IF (printlev>=4) WRITE(numout,*) 'Leaving fire'
518
519	END SUBROUTINE fire
520
521
522	!! ================================================================================================================================
523	!! FUNCTION : firefrac_func
524	!!
525	!>\BRIEF Calculate the annual fire fraction using long term fire index
526	!!
527	!! DESCRIPTION : None
528	!!
529	!! RECENT CHANGE(S): None
530	!!
531	!! RETURN VALUE : annual fire fraction (::firefrac_result)
532	!!
533	!! REFERENCES :
534	!! - Thonicke K., Venevsky S., Sitch S., and Cramer W. (2001) The role of fire
535	!! disturbance for global vegetation dynamics: coupling fire into a Dynamic
536	!! Global Vegetation Model, Global Ecology & Biogeography, 10, 661-677.\n
537	!!
538	!! FLOWCHART : None
539	!! \n
540	!_ ================================================================================================================================
541
542	FUNCTION firefrac_func (npts, x) RESULT (firefrac_result)
543
544	!! 0. Variable and parameter declaration
545
546	!! 0.1 Input variables
547
548	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
549	REAL(r_std), DIMENSION(npts), INTENT(in) :: x !! fire index (unitless, 0-1)
550
551	!! 0.2 Output variables
552
553	REAL(r_std), DIMENSION(npts) :: firefrac_result !! fire fraction (unitless, 0-1)
554
555	!! 0.3 Modified variables
556
557	!! 0.4 Local variables
558	REAL(r_std), DIMENSION(npts) :: xm1 !! intermediate variable
559	!_ ================================================================================================================================
560
561	!! 1. Function
562
563	xm1(:) = x(:) - 1.
564
565	! this equation is from Thonicke et al. (2001) equation (9), it use the fire index as input to calculate annual fire fraction.
566	! but with different parameters with the source literature.
567	!! \latexonly
568	!! \input{lpj_fire_2.tex}
569	!! \endlatexonly
570	!! Where, s is the long term fire index and A(s) is the annual fire fraction.\n
571	!??! here we use a different parameter with K. Thonicke et al. (2001)
572	!??! it's not clear if these parameters are correct.
573	firefrac_result(:) = &
574	& x(:) * EXP( xm1(:) / ( -firefrac_coeff(4)xm1(:)xm1(:)*xm1(:) + &
575	& firefrac_coeff(3)xm1(:)xm1(:) + &
576	& firefrac_coeff(2)*xm1(:) + &
577	& firefrac_coeff(1) ) )
578
579	END FUNCTION firefrac_func
580
581	END MODULE lpj_fire

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

Download in other formats: