Context Navigation

lpj_fire.f90 @ 7193

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

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

Context Navigation

source: branches/publications/ORCHIDEE-MICT-OP-r6850/src_stomate/lpj_fire.f90 @ 7193

Download in other formats: