Context Navigation

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

Last change on this file since 7346 was 2860, checked in by nicolas.vuichard, 9 years ago
Merge with the trunk is done. It compiles
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 31.2 KB

Line
1	! =================================================================================================================================
2	! MODULE : lpj_fire
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.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, &
121	litterhum_daily, t2m_daily, lignin_struc,lignin_wood,veget_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), INTENT(in) :: litterhum_daily !! Daily litter moisture (unitless, 0-1)
135	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_daily !! Daily 2 meter temperature (K)
136	REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(in) :: lignin_struc !! Ratio Lignin/Carbon in structural above
137	!! and below ground litter
138	!! @tex $(gC gC^{-1})$ @endtex
139	REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(in) :: lignin_wood !! Ratio Lignine/Carbon in woody litter for above and
140	!! and below ground litter
141	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max !! Maximum fraction of vegetation type including
142	!! non-biological fraction (0-1, unitless)
143
144	!! 0.2 Output variables
145
146	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: co2_fire !! Carbon emitted into the atmosphere by
147	!! fire, including both living and dead
148	!! biomass
149	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
150	!! 0.3 Modified variables
151
152	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: fireindex !! Probability of fire; in the code means
153	!! long term fire index (unitless, 0-1)
154	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: firelitter !! Total natural litter
155	!! available for burning above the ground
156	!! @tex $(gC m^{-2})$ @endtex
157	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout) :: biomass !! Biomass @tex $(gC m^{-2})$ @endtex
158	REAl(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind !! Density of individuals
159	!! @tex $(m^{-2})$ @endtex
160	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) :: litter !! Metabolic and structural litter, above
161	!! and below ground
162	!! @tex $(gC m^{-2})$ @endtex
163	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: dead_leaves !! Dead leaves on ground, per PFT,
164	!! metabolic and structural
165	!! @tex $(gC m^{-2})$ @endtex
166	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(inout):: bm_to_litter !! Biomass entering the litter pool
167	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
168	REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(inout) :: MatrixA !! Matrix containing the fluxes between the carbon pools
169	!! here, it is called only if disable_fire = n
170	!! once a day
171	!! @tex $(gC.m^2.day^{-1})$ @endtex
172
173
174	!! 0.4 Local variables
175
176	REAL(r_std), DIMENSION(npts) :: fire_disturb !! Actual fire disturbance fraction
177	!! after consideration of inherent fire
178	!! resistance of different PFTs
179	!! (unitless, 0-1)
180	REAL(r_std), DIMENSION(npts,nvm) :: firedeath !! In case of activation of dynamic
181	!! vegetation, the daily fraction of
182	!! burned individuals (unitless, 0-1)
183	REAL(r_std), DIMENSION(npts) :: moistlimit !! Moisture for fire inhibition per PFT !!
184	!! (unitless, 0-1) ; temporary variable
185	!! for each PFT in the loop over #PFT ?
186	!! it's not fire inhibition, it's
187	!! "moisture of extinction" by the original
188	!! reference paper, and a more "fire
189	!! professional name" it's better to change
190	!! the name of moistlimit
191	REAL(r_std), DIMENSION(npts) :: litter_above !! Total aboveground litter per PFT !!
192	!! @tex $(gC m^{-2})$ @endtex
193	REAL(r_std), DIMENSION(npts,nvm) :: fireindex_daily !! Daily fire index (unitless, 0-1)
194	REAL(r_std), DIMENSION(npts, nvm) :: firefrac !! Daily fire burning fraction on ground
195	!! (unitless, 0-1)
196	REAL(r_std), DIMENSION(npts) :: struc_residual !! Fraction of structural litter not burned
197	!! (thus residual), depending on strutural
198	!! litter lignin content (unitless, 0-1)
199	REAL(r_std), DIMENSION(npts) :: wood_residual !! Fraction of woody litter not burned
200	!! (thus residual), depending on strutural
201	!! litter lignin content (unitless, 0-1)
202	REAL(r_std), DIMENSION(npts) :: residue_frac !! Fraction of Fuel (either biomass or litter) not
203	!! burned (unitless, 0-1)
204	REAL(r_std), DIMENSION(npts) :: x !! Intermediate variable
205	REAL(r_std), DIMENSION(npts) :: aff !! Annual fire fraction (unitless, 0-1)
206	INTEGER(i_std) :: j,k,m !! Indeces
207	REAL(r_std), DIMENSION(npts,nvm) :: nox_emission !! NOx emitted into the atmosphere by fire (living and dead biomass)
208	!! (in gN/m**2/day)
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	WRITE(numout,*) 'fire:'
224
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	WRITE(numout,*) ' labile reserve: ', co2frac(ilabile)
237
238	WRITE(numout,) ' > critical litter quantity (gC/m*2): ',litter_crit
239	WRITE(numout,*) ' > We calculate a fire probability on agricultural ground, but'
240	WRITE(numout,*) ' the effective fire fraction is zero.'
241
242	firstcall_fire = .FALSE.
243
244	ENDIF
245
246	!! 1.4 Initialize output
247	co2_fire(:,:) = zero
248	firedeath(:,:) = zero
249	fireindex_daily(:,:) = zero
250	firefrac(:,:) = zero
251	nox_emission(:,:)=zero
252	!! 2. Determine fire probability
253
254	! Calculate probability that crops are not burned for the moment
255	! Calculate long-term aboveground litter
256
257	! a long loop over PFT
258	DO j = 2,nvm !loop over #PFT
259
260	!! 2.1 Total above ground litter
261	! Total litter above the ground for the PFT under consideration
262	litter_above(:) = litter(:,imetabolic,j,iabove,icarbon) + &
263	litter(:,istructural,j,iabove,icarbon)
264
265
266	!! 2.2 Calculate the PFT litter amount weighted by moisture
267	! If litter moisture is higher than moisture of extinction, fire is not possible.
268	moistlimit(:) = zero
269
270	WHERE ( (t2m_daily(:) .GT. ZeroCelsius) .AND. (litter_above(:) .GT. min_stomate) )
271	!?? the calculation here is redundant? as the part before flam(j) is 1?
272	!?? see first comment.. litterpart to be removed ?
273	moistlimit(:) = flam(j)
274
275	ENDWHERE
276
277	!! 2.3 Calculate daily fire index
278	! Calculate daily fire index
279	! \latexonly
280	! \input{lpj_fire_1.tex}
281	! \endlatexonly
282	! Where, m is the daily litter moisture, m_e is the moisture of extinction, p(m)
283	! is probability of fire (i.e. the daily fire index).\n
284	WHERE ( moistlimit(:) .GT. min_stomate )
285	x(:) = litterhum_daily(:)/moistlimit(:)
286	fireindex_daily(:,j) = EXP( - pi * x(:) * x(:) )
287	ELSEWHERE
288	fireindex_daily(:,j) = zero
289	ENDWHERE
290
291	!! 2.4 Calculate long-term fire index
292	! Calculate long-term fire index which is the mean probability of fire
293	fireindex(:,j) = ((tau_fire - dt) * fireindex(:,j) + (dt) * fireindex_daily(:,j)) / tau_fire
294
295	!! 2.5 Calculate long term aboveground litter that are available for burning
296	firelitter(:,j) = &
297	( ( tau_fire-dt ) * firelitter(:,j) + dt * litter_above(:) ) / tau_fire
298
299	!! 3. Calculate fire fraction from litter and fireindex
300
301	!! 3.1 Fire fraction from long term fire index for natural PFTs
302	! Transform the annual fire fraction to daily fraction. This is done by assuming that
303	! each day the fire fraction is the same.
304	aff(:) = firefrac_func (npts, fireindex(:,j))
305
306	! annual_fire_fraction = 1. - ( 1. - daily_fire_fraction )**365
307	! Thus, daily_fire_fraction = 1. - ( 1. - annual_fire_fraction )**(1/365)
308	! If annual firefrac<<1, then firefrac_daily = firefrac * dt/one_year
309	! This approximation avoids numerical problems.
310	! the variable 'un' is 1
311	IF(.NOT.disable_fire.AND.natural(j))THEN
312	WHERE ( aff(:) .GT. 0.1 )
313	firefrac(:,j) = un - ( un - aff(:) ) ** (dt/one_year)
314	ELSEWHERE
315	firefrac(:,j) = aff(:) * dt/one_year
316	ENDWHERE
317	ELSE
318	firefrac(:,j) = zero
319	ENDIF
320
321	! No fire if litter is below critical value
322	WHERE ( firelitter(:,j) .LT. litter_crit )
323	firefrac(:,j) = zero
324	ENDWHERE
325
326	! However, there is a minimum fire extent
327	firefrac(:,j) = MAX( 0.001_r_std * dt/one_year, firefrac(:,j) )
328
329	! if FIRE_DISABLE flag is set no fire
330	IF (disable_fire) firefrac = zero
331
332	!! 3.2 For agricultural ground, no fire is burned for the moment
333	IF ( .NOT. natural(j)) firefrac(:,j) = zero
334
335	!! 4. Determine fire impact
336
337	! Calculate fire disturbance fraction for each PFT, and fire emissions due
338	! to grasses.
339
340	!! 4.1 Tree and grass live biomass
341	! Tree live biomass is not burned in fire.
342	! However, in the dynamic vegetation module, tree individual density
343	! will be updated after fire. The fraction of tree individuals that are
344	! supposed to die in fire is the fire fraction multiplied by the tree PFT fire
345	! resistance which reflect survivorship of different tree PFTs during fire.
346	IF ( is_tree(j) ) THEN
347
348	!! 4.1.1 Disturban,ce factor for trees
349	! Trees are disturbed over the whole year. A resistance factor is
350	! used to reflect survivorship of different tree PFTs during fire.
351	fire_disturb(:) = ( un - resist(j) ) * firefrac(:,j)
352
353	ELSE
354
355	!! 4.1.2 Disturbance factor for grasses
356	! Grass is not disturbed outside the growing season thus grass biomass
357	! is only burned during the growing season.
358	WHERE ( biomass(:,j,ileaf,icarbon) .GT. min_stomate )
359
360	fire_disturb(:) = ( un - resist(j) ) * firefrac(:,j)
361
362	ELSEWHERE
363
364	fire_disturb(:) = zero
365
366	ENDWHERE
367
368	ENDIF
369
370	!! 4.2 Burn live biomass
371	! The burned part goes to either the atmposhere or litter
372	DO k = 1, nparts
373
374	! for tree PFT, all biomass compartments are burned.
375	! for grass biomass compartments, all are burned except root and carbon reserve
376	! IF concerning PFT is tree; OR (the PFT is grass, but not root or carbon reserve biomass); then it's burned.
377	IF ( .NOT. ( ( .NOT. is_tree(j) ) .AND. ( ( k.EQ.iroot ) .OR. ( k.EQ.icarbres ) ) ) ) THEN
378
379	!! 4.2.1 Fraction to the atmosphere.
380	co2_fire(:,j) = co2_fire(:,j)+ biomass(:,j,k,icarbon) * fire_disturb(:) * co2frac(k) / dt
381
382	nox_emission(:,j) = nox_emission(:,j) + biomass(:,j,k,initrogen) &
383	* fire_disturb(:) * co2frac(k) / dt
384
385	!! 4.2.2 Determine the residue
386	! Residual is expressed as a fraction
387	residue_frac(:) = fire_disturb(:) * ( un - co2frac(k) )
388
389	!! 4.2.2.3 The rest (largest part) of the residue becomes litter.
390	bm_to_litter(:,j,k,icarbon) = bm_to_litter(:,j,k,icarbon) + residue_frac(:) &
391	* biomass(:,j,k,icarbon)
392
393	bm_to_litter(:,j,k,initrogen) = bm_to_litter(:,j,k,initrogen) + residue_frac(:) &
394	* biomass(:,j,k,initrogen)
395	ENDIF
396
397	ENDDO
398
399
400	!! 4.3 Update biomass pool after burning
401
402	!! 4.3.1 Decrease biomass amount except for grass roots and carbon reserve.
403	DO m = 1, nelements ! Loop over # elements
404
405	DO k = 1, nparts
406
407	! **2
408	IF ( .NOT. ( ( .NOT. is_tree(j) ) .AND. ( ( k.EQ.iroot ) .OR. ( k.EQ.icarbres) ) ) ) THEN
409
410	biomass(:,j,k,m) = ( un - fire_disturb(:) ) * biomass(:,j,k,m)
411
412	ENDIF
413
414	ENDDO
415
416	END DO ! End Loop over # elements
417
418
419	!! 4.3.2 If vegetation is dynamic, then decrease the density of tree individuals.
420	IF ( (ok_dgvm .OR. .NOT.lpj_gap_const_mort) .AND. is_tree(j) ) THEN
421
422	firedeath(:,j) = fire_disturb(:) / dt
423
424	ind(:,j) = ( un - fire_disturb(:) ) * ind(:,j)
425
426	ENDIF
427
428	ENDDO ! loop over #PFT
429
430	!! 5. Burn litter
431
432	! All litter (either forest or grass) is burned throughout the whole year
433	! Metabolic litter is totally burned. Burning fraction of structural litter is related
434	! with its lignin content. The burned part either goes to the atmosphere
435	! or remains in litter as unburned residue.
436	DO j = 2,nvm !loop over #PFT
437
438	!! 5.1 Burn exposed metabolic litter
439	! Exposed metabolic litter burns completely and goes directly into the atmosphere as CO2.
440
441	!! 5.1.1 CO2 flux from litter burning
442	! Flux expressed in gC m^{-2}dtslow^{-1}
443	co2_fire(:,j) = co2_fire(:,j) + litter(:,imetabolic,j,iabove,icarbon) * &
444	firefrac(:,j) / dt
445
446	nox_emission(:,j) = nox_emission(:,j) + litter(:,imetabolic,j,iabove,initrogen) * &
447	firefrac(:,j) / dt
448
449	!! 5.1.2 Decrease metabolic litter
450	DO m = 1,nelements
451	litter(:,imetabolic,j,iabove,m) = litter(:,imetabolic,j,iabove,m) * &
452	( un - firefrac(:,j) )
453	END DO
454
455	!! 5.2 Burning exposed structural litter
456	! The fraction that goes to the atmosphere is related to its lignin content. The remaining unburned
457	! residues remain as litter.
458
459	!! 5.2.1 Incomplete burning of exposed structural litter
460	! Fraction of structural litter that does not burn completly. This fraction depends on lignin
461	! content (lignin_struc).
462	struc_residual(:) = fire_resist_lignin * lignin_struc(:,j,iabove)
463
464	!! 5.2.2 CO2 flux from litter burning
465	! Flux expressed in gC m^{-2}dtslow^{-1}
466	co2_fire(:,j) = co2_fire(:,j) + &
467	litter(:,istructural,j,iabove,icarbon) * firefrac(:,j) * &
468	( un - struc_residual(:) )/ dt
469
470	nox_emission(:,j) = nox_emission(:,j) + &
471	litter(:,istructural,j,iabove,initrogen) * firefrac(:,j) * &
472	( 1. - struc_residual(:) )/ dt
473
474	!! 5.2.3 Determine residue
475	! The residue is litter that undergoes fire, but is not transformed into CO2
476	residue_frac(:) = firefrac(:,j) * struc_residual(:)
477
478	!! 5.2.6 TResidue remaining litter
479	! The rest (largest part) of the residue remains litter. Remaining
480	! litter is the sum of this and of the litter which has not undergone a fire.
481	litter(:,istructural,j,iabove,icarbon) = &
482	litter(:,istructural,j,iabove,icarbon) * ( un - firefrac(:,j) ) + &
483	residue_frac(:) * litter(:,istructural,j,iabove,icarbon)
484
485	litter(:,istructural,j,iabove,initrogen) = &
486	litter(:,istructural,j,iabove,initrogen) * ( un - firefrac(:,j) ) + &
487	residue_frac(:) * litter(:,istructural,j,iabove,initrogen)
488
489	!! 5.3 Burning exposed woody litter is not totally transformed into CO2.
490
491	!! 5.3.1 Fraction of exposed woody litter that does not
492	! burn totally should depend on lignin content (lignin_wood).
493
494	wood_residual(:) = fire_resist_lignin * lignin_wood(:,j,iabove)
495
496	!! 5.3.2 CO2 flux, in gC/(m**2 of total ground)/day.
497
498	co2_fire(:,j) = co2_fire(:,j) + &
499	litter(:,iwoody,j,iabove,icarbon) * firefrac(:,j) * &
500	( 1. - wood_residual(:) )/ dt
501
502	nox_emission(:,j) = nox_emission(:,j) + &
503	litter(:,iwoody,j,iabove,initrogen) * firefrac(:,j) * &
504	( 1. - wood_residual(:) )/ dt
505
506	!! 5.3.3 determine residue (litter that undergoes fire, but is not transformed
507	! into CO2)
508
509	residue_frac(:) = firefrac(:,j) * wood_residual(:)
510
511	!! 5.3.5 The rest (largest part) of the residue remains litter. Remaining litter
512	! is the sum of this and of the litter which has not undergone a fire.
513
514	litter(:,iwoody,j,iabove,icarbon) = &
515	litter(:,iwoody,j,iabove,icarbon) * ( 1. - firefrac(:,j) ) + &
516	residue_frac(:) * litter(:,iwoody,j,iabove,icarbon)
517
518	litter(:,iwoody,j,iabove,initrogen) = &
519	litter(:,iwoody,j,iabove,initrogen) * ( 1. - firefrac(:,j) ) + &
520	residue_frac(:) * litter(:,iwoody,j,iabove,initrogen)
521
522	!! 5.2.7 Update matrix A for analytical spin-up (only if SPINUP_ANALYTIC is activated)
523
524	IF (spinup_analytic) THEN
525
526	! litter structural above
527	MatrixA(:,j,istructural_above,istructural_above) = MatrixA(:,j,istructural_above,istructural_above) - firefrac(:,j) &
528	+ firefrac(:,j)* struc_residual(:)
529
530	! litter_metabolic above
531	MatrixA(:,j,imetabolic_above,imetabolic_above) = MatrixA(:,j,imetabolic_above,imetabolic_above) - firefrac(:,j)
532
533	! litter woody above
534	MatrixA(:,j,iwoody_above,iwoody_above) = MatrixA(:,j,iwoody_above,iwoody_above) - firefrac(:,j) &
535	+ firefrac(:,j)* wood_residual(:)
536
537	!DS! This is the exact formulation. The above is a simplified formulation.
538	! MatrixA(:,j,istructural_above,istructural_above) = MatrixA(:,j,istructural_above,istructural_above)*(1. - firefrac(:,j) + &
539	! firefrac(:,j)* struc_residual(:)
540	! MatrixA(:,j,imetabolic_above,imetabolic_above) = MatrixA(:,j,imetabolic_above,imetabolic_above)*(1 - firefrac(:,j))
541
542	ENDIF !(spinup_analytic)
543
544	ENDDO ! loop over #PFT
545
546	!! 5.3 Update exposed dead leaves (leaf litter) on the ground
547	! Dead leaves are supposed to burn completely in fire even their structural part.
548	DO j = 2,nvm
549
550	DO k = 1, nlitt
551	dead_leaves(:,j,k) = dead_leaves(:,j,k) * ( un - firefrac(:,j) )
552	ENDDO
553
554	ENDDO
555
556	!! 6. write out history variables
557
558	! output in 1/dtslow
559	firefrac(:,:) = firefrac(:,:) / dt
560
561	CALL xios_orchidee_send_field("FIREFRAC",firefrac(:,:))
562	CALL xios_orchidee_send_field("FIREDEATH",firedeath(:,:))
563	CALL xios_orchidee_send_field("NOX_FIRE_EMISSION",nox_emission(:,:))
564
565	CALL histwrite_p (hist_id_stomate, 'FIREFRAC', itime, &
566	firefrac(:,:), npts*nvm, horipft_index)
567	CALL histwrite_p (hist_id_stomate, 'FIREDEATH', itime, &
568	firedeath(:,:), npts*nvm, horipft_index)
569	CALL histwrite_p (hist_id_stomate, 'NOX_FIRE_EMISSION', itime, &
570	nox_emission(:,:), npts*nvm, horipft_index)
571
572	IF (printlev>=4) WRITE(numout,*) 'Leaving fire'
573
574	END SUBROUTINE fire
575
576
577	!! ================================================================================================================================
578	!! FUNCTION : firefrac_func
579	!!
580	!>\BRIEF Calculate the annual fire fraction using long term fire index
581	!!
582	!! DESCRIPTION : None
583	!!
584	!! RECENT CHANGE(S): None
585	!!
586	!! RETURN VALUE : annual fire fraction (::firefrac_result)
587	!!
588	!! REFERENCES :
589	!! - Thonicke K., Venevsky S., Sitch S., and Cramer W. (2001) The role of fire
590	!! disturbance for global vegetation dynamics: coupling fire into a Dynamic
591	!! Global Vegetation Model, Global Ecology & Biogeography, 10, 661-677.\n
592	!!
593	!! FLOWCHART : None
594	!! \n
595	!_ ================================================================================================================================
596
597	FUNCTION firefrac_func (npts, x) RESULT (firefrac_result)
598
599	!! 0. Variable and parameter declaration
600
601	!! 0.1 Input variables
602
603	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
604	REAL(r_std), DIMENSION(npts), INTENT(in) :: x !! fire index (unitless, 0-1)
605
606	!! 0.2 Output variables
607
608	REAL(r_std), DIMENSION(npts) :: firefrac_result !! fire fraction (unitless, 0-1)
609
610	!! 0.3 Modified variables
611
612	!! 0.4 Local variables
613	REAL(r_std), DIMENSION(npts) :: xm1 !! intermediate variable
614	!_ ================================================================================================================================
615
616	!! 1. Function
617
618	xm1(:) = x(:) - 1.
619
620	! this equation is from Thonicke et al. (2001) equation (9), it use the fire index as input to calculate annual fire fraction.
621	! but with different parameters with the source literature.
622	!! \latexonly
623	!! \input{lpj_fire_2.tex}
624	!! \endlatexonly
625	!! Where, s is the long term fire index and A(s) is the annual fire fraction.\n
626	!??! here we use a different parameter with K. Thonicke et al. (2001)
627	!??! it's not clear if these parameters are correct.
628	firefrac_result(:) = &
629	& x(:) * EXP( xm1(:) / ( -firefrac_coeff(4)xm1(:)xm1(:)*xm1(:) + &
630	& firefrac_coeff(3)xm1(:)xm1(:) + &
631	& firefrac_coeff(2)*xm1(:) + &
632	& firefrac_coeff(1) ) )
633
634	END FUNCTION firefrac_func
635
636	END MODULE lpj_fire

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

Download in other formats: