Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_stomate/lpj_spitfire.f90 @ 6737

Last change on this file since 6737 was 4080, checked in by albert.jornet, 7 years ago
New: add spitfire history outputs. Done by Chao
File size: 115.0 KB

Line
1	! ==============================================================================================================================
2	! MODULE : lpj_spitfire
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF This module simulates fires in
10	!! natural vegetations in a prognostic way.
11	!!
12	!!
13	!!\n DESCRIPTION : Use daily meterological
14	!! information combined with information for fuel availability to simulate
15	!! burned area and gas emissions in fire. Only fires in natural PFTs have been
16	!! simulated.
17	!!
18	!! REFERENCE(S)
19	!! - K.Thonicke, A.Spessa, I.C. Prentice, S.P.Harrison, L.Dong, and
20	!! C.Carmona-Moreno, 2010: The influence of vegetation, fire spread and fire
21	!! behaviour on biomass burning and trace gas emissions: results from a
22	!! process-based model. Biogeosciences, 7, 1991-2011.
23	!!
24	!! SVN :
25	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/DOC/ORCHIDEE/src_sechiba/enerbil.f90 $
26	!! $Date: 2012-00-00 12:22:59 +0100 $
27	!! $Revision: ??? $
28	!! \n
29	!_ ================================================================================================================================
30
31	!! Below are some remarks that further document model behavior or existing minor
32	!! issues in the model, they're minor in nature and do not affect large-scale
33	!! model performance.
34	!! ----------------------- Remarks ------------------------------------------------
35	!! Date 2015/04/13
36	! 1. There are few places where clarifications are needed, mainmy minor, they're
37	! marked as ??, or [UNCLEAR]
38	! 2. some varialbes are not actually used in some SUBROUTINEs or the whole module.
39	! and they're designated as "NOT USED"
40	! 3. some sources of parameters need to be check. they're marked as "SOURCE"
41	! 4. it's been tested that observed BA will be read as they are (the monthly
42	! ba will be repeated by day number of concerning month) if the external observed
43	! burned area is forced into the model. The only strange behavior is that the
44	! last day of the year seems to repeat the first month value (it's also the case
45	! for lightning), (cf./home/orchidee03/ychao/OUTPUT_A5F_STD_WITHFIRE_OLDmap_NA5070
46	! /FIRETEST/TEST_BA_LIGHTN/BATEST/observed_ba.txt)
47	! When foring BA, if we assume that BA is distributed equally across the month,
48	! we should divide the monthly BA in GFED3 by number of days in that month
49	! before feed the data into model.
50	! 4. it's been test (2012/06/03) for 0.5dX0.5d lightning data, it behaves the same
51	! way exactly as burned area. this means lightning data must be provided as daily
52	! value.
53	! 5. as popdensity is given one value for one year. this value is read directly into
54	! model (model use constant popdens value throughout the whole simulation year)
55	!! --------------------------------------------------------------------------------
56
57	!!------------ Explanations for forcing external Burned area ----------------------
58	!! Date 2015/04/13
59	! 1. External burned area could be forced into model, currently only monthly external
60	! BA input is configuered, this monthly BA will be treated as to occur at the
61	! beginning day of each month. The flag to activate the external BA forcing is
62	! the read_observed_ba and ratio_flag variable, for regions where ratio_flag is bigger than 0,
63	! BA will be externally forced.
64	! 2. When external BA is not forced, baadjust_ratio could be used to adjust the
65	! simulated BA.
66	! 3. The consumption fraction in fires for fuels could also be imposed from
67	! external input, read_cf_coarse control the input of coarse fuel (corresponding
68	! to 1000hr fuel in the model); read_cf_fine control the input for the fine
69	! fuel (corresponding to live grass, 1hr/10hr/100hr fuel in the model).
70	!
71	! This behaviour is very flexible, one has to read the code to fully understand it.
72	!!--------------------------------------------------------------------------------
73
74	MODULE lpj_spitfire
75
76	! modules used:
77	USE xios_orchidee
78	USE stomate_data
79	USE pft_parameters
80	USE constantes
81	USE ioipsl_para
82
83	IMPLICIT NONE
84
85	! private & public routines
86	PRIVATE
87	PUBLIC spitfire,spitfire_clear
88
89	! first call
90	LOGICAL, SAVE :: firstcall = .TRUE.
91
92	CONTAINS
93
94
95	SUBROUTINE spitfire_clear
96	firstcall = .TRUE.
97	END SUBROUTINE spitfire_clear
98
99	!! ================================================================================================================================
100	!! SUBROUTINE : SPITFIRE
101	!!
102	!>\BRIEF A prognostic fire module. It uses daily temperature and
103	!! precipitation to calculate fire danger index, combined with ignition sources to
104	!! derive number of fires. Then Burned area, fuel consumption and tree crown
105	!! scorching, tree mortality were calculated.
106	!!
107	!! DESCRIPTION :
108	!!
109	!! The module completes the following tasks:
110	!! 1. Calculate fire danger index from daily temperature and precipitation using
111	!! Nesterov Index (NI). The fuel mositure and moisture of extinction are compared
112	!! to derive a daily Fire Danger Index (FDI). Fuel moiture is calculated using
113	!! Nesterov Index by weighting among differnet types of fuels. The fuel type with
114	!! a bigger surface-area-to-volume ratio receives less weight in determining the
115	!! fuel moisture state.
116	!! 2. Use daily FDI and ignition source data to derive daily fire numbers.
117	!! 3. Apply Rothermel's equation to calculate fire Rate of Spread (ROS) for both
118	!! forward and backward. Rate of spread is determined in a way related with fire
119	!! reaction intensity(+), propagating flux ratio(+), wind speed(+), fuel bulk
120	!! density(-),effective heating number(-), and heat of pre-ignition(-).
121	!! 4. From daily FDI derive fire duration time, combined with rate of spread to
122	!! derive burned area.
123	!! 5. Fuel mositure calculated in the step of daily FDI can be again used to
124	!! determin consumption fraction of different types of fuels in fire.
125	!! 6. Calculate fire surface intensity using ROSf, fuel consumption, and burned
126	!! area, thus further to determin flame height and fire damage to tree crown.
127	!! 7. Fire caused tree mortality are determined based on 2 considerations, the
128	!! damage of flame to tree crown, and damage of fire residence time to tree barks.
129	!! crowning scorching biomass release as emissions, while unburned dead tree
130	!! biomass part transfer to litter pool.
131	!!
132	!! In summary, fire weather (temperature, preciptation, wind speed) plays a central
133	!! role in this module, it determines fuel moisture, which futher determines
134	!! daily fire danger index (FDI) --> nubmer of fires --> fire duration --> fire
135	!! size --> burned area --> fire intensity, residence time --> tree demage
136	!!
137	!! Important Note:
138	!! 1. over one day and within one pixel, fire has unified size and BA=Size X
139	!! Number.
140	!! 2. Fuel consumption and burned area are calculated independently.
141	!!
142	!! RECENT CHANGE(S): None
143	!!
144	!! MAIN OUTPUT VARIABLE(S): fire_frac, co2_fire, and trace gas emissions.
145	!!
146	!! REFERENCES :
147	!! - K.Thonicke, A.Spessa, I.C. Prentice, S.P.Harrison, L.Dong, and
148	!! C.Carmona-Moreno, 2010: The influence of vegetation, fire spread and fire
149	!! behaviour on biomass burning and trace gas emissions: results from a
150	!! process-based model. Biogeosciences, 7, 1991-2011.
151	!! - Rothermel, R. C.: A Mathematical Model for Predicting Fire
152	!! Spread in Wildland Fuels, Intermountain Forest and Range Experiment
153	!! Station, Forest Service, US Dept. of Agriculture, Ogden, UtahUSDA
154	!! Forest Service General Technical Report INT-115, 1972.
155	!! - Pyne, S. J., Andrews, P. L., and Laven, R. D.: Introduction to wildland
156	!! fire, 2 edition, Wiley, New York, 769 pp., 1996
157	!!
158	!! FLOWCHART : None
159	!! \n
160	!_ ================================================================================================================================
161
162
163	SUBROUTINE spitfire(npts, dt_days, veget_max_org,resolution,contfrac, &
164	PFTpresent,t2m_min_daily,t2m_max_daily, &
165	precip_daily,wspeed_daily,soilhum_daily, &
166	lightn,litter,ni_acc,fire_numday, &
167	fuel_1hr,fuel_10hr,fuel_100hr, &
168	fuel_1000hr,ind,biomass,popd,a_nd,height, &
169	read_observed_ba, observed_ba,read_cf_fine,cf_fine, &
170	read_cf_coarse,cf_coarse,read_ratio_flag, &
171	ratio_flag,read_ratio,baadjust_ratio,date, &
172	bm_to_litter,co2_fire, &
173	lcc,bafrac_deforest_accu,emideforest_litter_accu,emideforest_biomass_accu,&
174	deforest_litter_remain,deforest_biomass_remain,&
175	def_fuel_1hr_remain,def_fuel_10hr_remain,&
176	def_fuel_100hr_remain,def_fuel_1000hr_remain)
177
178
179
180	! PARAMETERS
181	INTEGER, parameter :: ntrace=6 !!number of different trace gas types included
182	REAL(r_std), parameter :: MINER_TOT=0.055 !!total mineral content, source: Thonicke et al., 2010. p1010 Table A1 Row 5
183	REAL(r_std), parameter :: H=18000.0 !!calorific heat content (kJ kg^{-1}), source: Thonicke et al. (2010) Appendix A
184	REAL(r_std), parameter :: sigma_1hr=66.0 !!surface-area-to-volume ratio (cm^{2} cm^{-3}), source: (USDA 1978, Harvey et al. 1997);
185	!! These two sources are given in a working doc by K. Thonicke, the exact sources unknown.
186	REAL(r_std), parameter :: sigma_10hr=3.58 !!same as above
187	REAL(r_std), parameter :: sigma_100hr=0.98 !!same as above
188	REAL(r_std), parameter :: sigma_livegrass=80.0 !! NOT USED
189	REAL(r_std), parameter :: me_livegrass=0.2 !!moisture of extinction for live grass fuel (unitless)
190
191	! SOURCE this two parameters are said to follow Brown et al. 1981 but I don't
192	!have the paper.
193	REAL(r_std), parameter :: fbd_a = 1.2 !!scalar used to transform FBD of 10hr fuels to 1hr equivs, FBD: fuel bulk density (kg m^{-3})
194	REAL(r_std), parameter :: fbd_b = 1.4 !!scalar used to transform FBD of 100hr fuels to 1hr equivs
195	!SOURCE
196	REAL(r_std), parameter :: fbd_C3_livegrass=4.0 !!kg/m3
197	REAL(r_std), parameter :: fbd_C4_livegrass=4.0 !!kg/m3
198
199	!INPUT VARIABLES
200	INTEGER, INTENT(in) :: npts !! Domain size
201	REAL(r_std), INTENT(in) :: dt_days !! time step of Stomate in days
202	INTEGER(i_std),INTENT(in) :: date !! Date (days)
203	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max_org !! veget_max
204	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: resolution !! resolution at each grid point in m (1=E-W, 2=N-S)
205	REAL(r_std),DIMENSION (npts), INTENT (in) :: contfrac !! fraction of continent in the grid
206	LOGICAL, DIMENSION(npts,nvm), INTENT(inout) :: PFTpresent !! logical variables to check whether a PFT present, defined in stomate.f90
207	REAL(r_std),DIMENSION(npts), INTENT(in) :: t2m_min_daily !! Daily minimum 2 meter temperatures (K)
208	REAL(r_std),DIMENSION(npts), INTENT(in) :: t2m_max_daily !! Daily maximum 2 meter temperatures (K)
209	REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_daily !! daily precip (mm/day)
210	REAL(r_std), DIMENSION(npts), INTENT(in) :: wspeed_daily !! wind speed [m/s]
211	REAL(r_std), DIMENSION(npts), INTENT(in) :: lightn !! number of daily lightning [flashes per km2 per day]]
212	!! lightn declared and allocated in slowproc.f90
213	LOGICAL, INTENT (in) :: read_observed_ba !! Flag for read in observed burned area, i.e., model is forced by external BA
214	REAL(r_std),DIMENSION (npts), INTENT (in) :: observed_ba !! Observed burned area, currently must be in monthly timestep
215
216
217	LOGICAL, INTENT (in) :: read_cf_coarse !! Flag to read external forced consumption fraction for coarse fuels, will be used
218	!! on the 1000hr fuel in the model.
219	REAL(r_std),DIMENSION (npts), INTENT (in) :: cf_coarse !! Externally forced fire consumption fraction for coarse fuel, applied on 1000hr fuel.
220	LOGICAL, INTENT (in) :: read_cf_fine !! Flag to read external forced consumption fraction for fine fuel.
221	REAL(r_std),DIMENSION (npts), INTENT (in) :: cf_fine !! Externally forced fire consumption fraction for fine fuel, applied on live grass,
222	!! and 1hr/10hr/100hr dead fuel.
223	LOGICAL, INTENT (in) :: read_ratio !! Flag for read in ratio to adjust simulated burned area.
224	REAL(r_std),DIMENSION (npts), INTENT (in) :: baadjust_ratio !! gridded input data used to adjust the simulated burned area, unitless.
225	LOGICAL, INTENT (in) :: read_ratio_flag !! Flag for read in observed burned area
226	REAL(r_std),DIMENSION (npts), INTENT (in) :: ratio_flag !! gridded flag value indicate where the burned area, conumption fractions
227	!! will be forced in the model. When bigger than 0, the burned area and
228	!! consumption fractions will be forced in the model.
229
230
231	REAL(r_std), DIMENSION(npts), INTENT(in) :: popd !! human population density [person per sqkm]
232	!! popd declared and allocated in slowproc.f90
233	REAL(r_std), DIMENSION(npts), INTENT(in) :: a_nd !! parameter for potential human-caused ignitions;
234	!! a_nd is parameterized in stomate_lpj before call of SPITFIRE
235	REAL(r_std), DIMENSION(npts), INTENT(in) :: soilhum_daily !! daily top layer soil humidity, here we need only the top soil
236	!! layer humdity, and this is addressed by
237	!! passing soilhum_daily(:,1) as soilhum_daily when we
238	!! call SPITFIRE in stomate_lpj.f90
239	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: height !! tree height in m
240	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lcc !! gross forest loss as fraction of griecell area
241
242
243	! MODIFIED VARIABLES
244	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs), INTENT(inout) :: litter !! metabolic and structural litter, above and belowground
245	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout) :: biomass !! biomass
246	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: fuel_1hr !! [gC^{m-2}]
247	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: fuel_10hr !! [gC^{m-2}]
248	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: fuel_100hr !! [gC^{m-2}]
249	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: fuel_1000hr !! [gC^{m-2}]
250	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: def_fuel_1hr_remain !! [gC^{m-2}]
251	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: def_fuel_10hr_remain !! [gC^{m-2}]
252	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: def_fuel_100hr_remain !! [gC^{m-2}]
253	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: def_fuel_1000hr_remain !! [gC^{m-2}]
254	REAL(r_std), DIMENSION(npts),INTENT(INOUT) :: ni_acc !! ni_acc declared and allocated in stomate.f90
255	REAL(r_std), DIMENSION(npts),INTENT(INOUT) :: fire_numday !! declared and allocated in stomate.f90
256	REAL(r_std), DIMENSION(npts,nvm),INTENT(inout) :: co2_fire !! total C emission in fire, including live crown scorching and
257	!! dead ground litter(tree/grass) + live grass leaf&fruit
258	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind !! Number of individuals / m2; defined in stomate.f90
259	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout) :: bm_to_litter !! Biomass transfer to litter
260	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: bafrac_deforest_accu !!accumulated deforestation fire fraction
261	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout) :: emideforest_litter_accu !!cumulative deforestation fire emissions from litter
262	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout) :: emideforest_biomass_accu !!cumulative deforestation fire emissions from live biomass burning
263	REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs), INTENT(inout) :: deforest_litter_remain
264	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout) :: deforest_biomass_remain
265
266
267
268	! LOCAL SPITFIRE VARIABLES
269	REAL(r_std), DIMENSION(npts) :: dfm_lg !! Daily fuel moisture for live grass, unitless, (0,1)
270	REAL(r_std), DIMENSION(npts) :: dfm !! Daily fuel moisture of dead fuels on ground excluding 1000hr fuel,(0,1)
271	REAL(r_std), DIMENSION(npts) :: ignition_efficiency !! Ignition efficiency as limited by fuel load, (0,1), unitless.
272	REAL(r_std), DIMENSION(npts) :: dfm_1hr !! daily fuel moisture for 1hr-fuel
273	REAL(r_std), DIMENSION(npts) :: fuel_1hr_total !! fuel_1hr_total is different from fuel_1hr(unit in C) and
274	!! it combines the fuel of all types and PFTs. unit: dry mass
275	REAL(r_std), DIMENSION(npts) :: fuel_10hr_total !! fuel_1hr/10hr/100hr/_total unit in dry mass for all natural PFTs
276	REAL(r_std), DIMENSION(npts) :: fuel_100hr_total
277	REAL(r_std), DIMENSION(npts) :: fuel_1000hr_total
278	REAL(r_std), DIMENSION(npts) :: natural_litter_ag_total !! total above-ground litter of natural PFTs
279	REAL(r_std), DIMENSION(npts) :: dead_fuel !! sum of 1hr/10hr/100hr dead fuel (g dry mass/m2) for all natural PFTs
280	REAL(r_std), DIMENSION(npts) :: dead_fuel_all !! sum of 1hr/10hr/100hr/1000hr dead fuel (g dry mass/m2) for all natural PFTs
281	REAL(r_std), DIMENSION(npts) :: livegrass !! sum the leaf bimoass of natural grassland PFTs. unit is in dry mass but not carbon!
282	REAL(r_std), DIMENSION(npts) :: dead_fuel_all_livegrass !! sum of 1hr/10hr/100hr/1000hr dead fuel (g dry mass/m2)
283	!! for all natural PFTs plus livegrass biomass
284	REAL(r_std), DIMENSION(npts) :: ratio_dead_fuel
285	REAL(r_std), DIMENSION(npts) :: ratio_live_fuel
286	REAL(r_std), DIMENSION(npts) :: me_litterfuel !! the moisture of extinction weighted by the
287	!! dead groud litter fuel of different natural PFTs, unitless
288	REAL(r_std), DIMENSION(npts) :: moist_extinction !! the moisture of extinction weighted by ground dead fuel and livegrass,
289	!! m_{e} in Eq.(8)
290	REAL(r_std), DIMENSION(npts) :: d_fdi !! daily fire danger index,climatic fire risk, unitless between 0. and 1.
291	REAL(r_std), DIMENSION(npts) :: net_fuel !! mineral content adjusted dead_fuel (in kg biomass!)
292	REAL(r_std), DIMENSION(npts) :: char_net_fuel !! total net fuel(mineral content corrected) by summing the dead fuel and live grass,
293	!! , in unit of kg dry mass per m2.
294	REAL(r_std), DIMENSION(npts) :: wind_speed !! wind_speed is Uforward in Thonicke et al. 2010; unit: m/min
295	REAL(r_std), DIMENSION(npts) :: wind_forward !! the name here is somewhat misleading; it's only for used in Eq(A5)
296	REAL(r_std), DIMENSION(npts) :: lb !! length-to-breadth ratio of the fire ellipse, unitless, range[1,8]
297	REAL(r_std), DIMENSION(npts) :: df
298	REAL(r_std), DIMENSION(npts) :: db
299	REAL(r_std), DIMENSION(npts) :: fire_durat !! fire duration (min)
300	REAL(r_std), DIMENSION(npts) :: ros_f !! forward fire spred rate(m min^{-1})
301	REAL(r_std), DIMENSION(npts) :: ros_b !! backward fire spread rate; unit: m min^{-1}
302	REAL(r_std), DIMENSION(npts) :: ratio_fbd !! proportion of dead fuel in class i PFT j as a proportion of total dead fuel
303	REAL(r_std), DIMENSION(npts) :: ratio_C3_livegrass !! ratio of C3 and C4 livegrass in total livegrass
304	REAL(r_std), DIMENSION(npts) :: ratio_C4_livegrass
305	REAL(r_std), DIMENSION(npts) :: dens_fuel_ave !! weighted fuel bulk density by different fuel types except 1000h fuel,
306	!! including only natural PFTs, in unit of kg dry mass per m3.
307	REAL(r_std), DIMENSION(npts) :: dens_livegrass_ave !! C3/C4 weighted livegrass fuel bulk density, kg m^{-3}
308	REAL(r_std), DIMENSION(npts) :: char_dens_fuel_ave !! livegrass and dead fule weighted fuel bulk density
309	!! (excluding 1000hr dead fuel), in unit of kg dry mass per m3.
310	REAL(r_std), DIMENSION(npts) :: sigma !! mass weighted mean fuel surface-area-to-volume ratio by different fuel types,
311	!! in unit of cm^{-1}, given that area in unit of cm^{2} and volume in cm^{3}
312	REAL(r_std), DIMENSION(npts) :: gamma !! tau_maxita_mita_s; used in calculation of tau_l for postfire mortality.
313	REAL(r_std), DIMENSION(npts) :: fpc_tree_total !! total tree fraction
314	REAL(r_std), DIMENSION(npts) :: fpc_grass_total !! total grass fraction (including natural and agricultural grassland)
315	REAL(r_std), DIMENSION(npts) :: area
316	REAL(r_std), DIMENSION(npts) :: wetness !! wetness of fuel, calculated as dfm/moisture_of_extinction
317	REAL(r_std), DIMENSION(npts) :: fuel_consum !! total fuel consumption including 1hr/10hr/100hr (not 1000hr!) fuel. [g mass m^{-2}]
318	REAL(r_std), DIMENSION(npts) :: d_i_surface !! surface fire intensity [kW m^{-1}]
319	REAL(r_std), DIMENSION(npts,nvm) :: disturb_crown !! a temporary variable, used to store the C emissions
320	!! from crown scorch of live tree crown.
321	REAL(r_std), DIMENSION(npts) :: fc_crown !! PFT weighted carbon emissions caused by crown scorching (gCm^{-2}).
322	!! only for history file writing.
323	!! this variable is PFT vcmax weighted sum of above one
324	REAL(r_std), DIMENSION(npts) :: d_area_burnt !! daily burned area (ha); local variable, only for history file writing.
325	REAL(r_std), DIMENSION(npts) :: ba_escape !! Escaped deforestation fire burned area (ha), it's included in d_area_burnt,
326	!! However is separated for being put into history file.
327	REAL(r_std), DIMENSION(npts) :: d_numfire !! local variable, only for history file writing
328	REAL(r_std), DIMENSION(npts) :: human_numfire !! local variable, only for history file writing.
329	REAL(r_std), DIMENSION(npts) :: lightn_numfire !! local variable, only for history file writing.
330	REAL(r_std), DIMENSION(npts) :: fire_frac !! burned fraction, unitless
331	REAL(r_std), DIMENSION(npts) :: area_land !! Land area within gridcel excluding water body,[hectar=10000 m2]
332	REAL(r_std), DIMENSION(npts) :: area_natural_veg !! Land area covered with vegetation within gridcell [hectar]
333	REAL(r_std), DIMENSION(npts) :: human_ign !! human ignitions, [1 day^{-1} km^{-2}]
334	REAL(r_std), DIMENSION(npts) :: lightn_efficiency !! ligtning fractions that reach ground with sufficient energy to ignite,unitless
335	REAL(r_std), DIMENSION(npts) :: human_suppression !! huamn suppression effects on fires, unitless
336	REAL(r_std), DIMENSION(npts) :: lightn_ign !! lightning ignitions, [1 day^{-1} km^{-2}]
337	REAL(r_std), DIMENSION(npts,ntrace) :: dcflux_trace
338	REAL(r_std), DIMENSION(npts,nvm,nlitt) :: fc_1hr !! fule consumption for 1hr fuel, [g biomass^m{-2}]
339	REAL(r_std), DIMENSION(npts,nvm,nlitt) :: fc_10hr !! fule consumption for 10hr fuel, [g biomass^m{-2}]
340	REAL(r_std), DIMENSION(npts,nvm,nlitt) :: fc_100hr !! fule consumption for 100hr fuel, in the unit of g dry mass/m**2 (C/0.45)
341	REAL(r_std), DIMENSION(npts,nvm,nlitt) :: fc_1000hr !! fule consumption for 1000hr fuel, [gC^{m-2}]
342	REAL(r_std), DIMENSION(npts,nvm) :: fc_lg !! leaf biomass consumption for GRASS, [gC m^{-2}]
343	REAL(r_std), DIMENSION(npts,nvm) :: fc_lf !! fruit biomass consumption for GRASS, [gC m^{-2}]
344	REAL(r_std), DIMENSION(npts) :: cf_ave !! the average fuel fire consumption fraction of 1hr/10hr/100hr fuel, unitless
345	REAL(r_std), DIMENSION(npts) :: cf_1hr !! fire fuel consumption fraction for 1hr fuel, unitless
346	REAL(r_std), DIMENSION(npts) :: cf_10hr !! fire fuel consumption fraction for 10hr fuel, unitless
347	REAL(r_std), DIMENSION(npts) :: cf_100hr !! fire fuel consumption fraction for 100hr fuel, unitless
348	REAL(r_std), DIMENSION(npts) :: cf_1000hr !! fire fuel consumption fraction for 1000hr fuel, unitless
349	REAL(r_std), DIMENSION(npts) :: cf_lg !! fire fule consumption for livegrass leaf and fruit biomass, unitless
350	REAL(r_std), DIMENSION(npts,nvm) :: tau_l !! the residence time of fire; used in Eq(19)
351	REAL(r_std), DIMENSION(npts,nvm) :: ck !!the proportion of crown affected by fire ( Eq17: CK=(SH-H+CL)/CL )
352	REAL(r_std), DIMENSION(npts,nvm) :: postf_mort !!
353	REAL(r_std), DIMENSION(npts,nvm) :: pm_tau
354	REAL(r_std), DIMENSION(npts,nvm) :: pm_ck
355	REAL(r_std), DIMENSION(npts,nvm) :: sh
356	REAL(r_std), DIMENSION(npts,nvm) :: nind_fa !! # trees affected by fire
357	REAL(r_std), DIMENSION(npts,nvm) :: nind_kill !! # trees killed by fire
358	REAL(r_std), DIMENSION(npts,nvm) :: prop_fa_nk !! propn of indivs fire affected but not killed
359	REAL(r_std), DIMENSION(npts,nvm) :: dcflux_fire_pft !! fire carbon flux to atmosphere; including both crown emssions
360	!! and litter consumption emissions
361	REAL(r_std), DIMENSION(npts,nvm,ntrace) :: dcflux_trace_pft
362	REAL(r_std), DIMENSION(npts,nvm) :: litter_consump_pft !! C emissions from ground litter burning (tree&grass) and grass leaf&fruit.
363	REAL(r_std), DIMENSION(npts) :: litter_consump !! Consumed litter during fire for each PFT,including ground dead litter for
364	!! trees and grasses, and live grass leaf and fruit burning. (gCm^{-2})
365	REAL(r_std), DIMENSION(npts,nvm) :: tau_c !! critical time to cambial kill
366	REAL(r_std), DIMENSION(npts) :: dia !! stem diameter (cm)
367	REAL(r_std), DIMENSION(npts,nvm) :: bt !! bark thickness (cm)
368	REAL(r_std), DIMENSION(npts) :: tree_kill_frac !! tree fraction that are killed in fire
369	REAL(r_std), DIMENSION(npts) :: mean_fire_size_original !! mean fire size before the correction of surface fire intensity
370	REAL(r_std), DIMENSION(npts) :: mean_fire_size !! mean fire size
371	REAL(r_std), DIMENSION(npts) :: natural_vegfrac !! the fraction of the grid cell that are occupied by natural vegetation.
372	REAL(r_std), DIMENSION(nvm,ntrace) :: ef_trace !! emission factors for trace gases
373	REAL(r_std), DIMENSION(npts,nvm) :: fc_1hr_carbon !! [gCm^{-2}]
374	REAL(r_std), DIMENSION(npts,nvm) :: fc_10hr_carbon !! [gCm^{-2}]
375	REAL(r_std), DIMENSION(npts,nvm) :: fc_100hr_carbon !! [gCm^{-2}]
376	REAL(r_std), DIMENSION(npts,nvm) :: fc_1000hr_carbon !! [gCm^{-2}]
377	REAL(r_std), DIMENSION(npts,nvm,nlitt) :: emideforest_fuel_1hr !! fire emissions from deforestation fires for 1hr fuel [gC^{m-2}]
378	REAL(r_std), DIMENSION(npts,nvm,nlitt) :: emideforest_fuel_10hr !! [gC^{m-2}]
379	REAL(r_std), DIMENSION(npts,nvm,nlitt) :: emideforest_fuel_100hr !! [gC^{m-2}]
380	REAL(r_std), DIMENSION(npts,nvm,nlitt) :: emideforest_fuel_1000hr !! [gC^{m-2}]
381	REAL(r_std), DIMENSION(npts,nvm) :: emi_nondef !! Fire carbon emissions not including deforestation fires, based on PFT area
382	REAL(r_std), DIMENSION(npts,nvm) :: emideforest_all !! Fire carbon emissions from deforestation fires, based on PFT area
383	REAL(r_std), DIMENSION(npts,nvm) :: bafrac_deforest !! Deforestation fire burned fraction as gridcell area
384	REAL(r_std), DIMENSION(npts,nvm,nlitt) :: emideforest_litter !! Deforestation fire emissions from litter burning [gCm^{-2}], based on ground area
385	REAL(r_std), DIMENSION(npts,nvm,nparts) :: emideforest_biomass !! live biomass emissions from deforestation fires [gCm^{-2}], based on ground area
386	REAL(r_std), DIMENSION(npts,nvm) :: veget_max !!
387	INTEGER :: j,x,m,k,n,i,ivm,ipts
388
389	REAL(r_std), DIMENSION(npts) :: proc_ba !! a proxy variable used to put the monthly observed BA at the first day
390	!! of each month.
391	REAL(r_std) :: surface_threshold=50.0
392	REAL(r_std) :: fuel_low_bound=444. !! below this ignition efficiency is 0 [dry matter/m2],
393	!! assuming a 0.45 carbon fraction.
394	REAL(r_std) :: fuel_high_bound=2222. !! above this ignition efficiency is 1 [dry matter/m2],
395	!! assuming a 0.45 carbon fraction.
396	!! the bounds sources are Arora & Boer (2005)
397	INTEGER(r_std),DIMENSION(12) :: monthday
398	monthday = (/1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335/)
399
400
401
402	! Initialise variables
403	disturb_crown(:,:)=0.0
404	sigma(:)=0.0
405	lb(:)=0.0
406
407	cf_lg(:)=0.0
408	cf_1hr(:)=0.0
409	cf_10hr(:)=0.0
410	cf_100hr(:)=0.0
411	cf_1000hr(:)=0.0
412	fpc_tree_total(:)=0.0
413	fpc_grass_total(:)=0.0
414	net_fuel(:)=0.0
415	dead_fuel(:)=0.0
416	ratio_dead_fuel(:)=0.0
417	ratio_live_fuel(:)=0.0
418	livegrass(:)=0.0
419	dead_fuel_all_livegrass(:)=0.0
420	fuel_1hr_total(:)=0.0
421	fuel_10hr_total(:)=0.0
422	fuel_100hr_total(:)=0.0
423	fuel_1000hr_total(:)=0
424	dens_fuel_ave(:)=0.0
425	d_area_burnt(:)=0.0
426	ba_escape(:)=0.0
427	fire_frac(:)=0.0
428	d_i_surface(:)=0.0
429	ck(:,:)=0.0
430	tau_l(:,:)=0.0
431	pm_tau(:,:)=0.0
432	postf_mort(:,:)=0.0
433	prop_fa_nk(:,:)=0.0
434	nind_fa(:,:)=0.0
435	nind_kill(:,:)=0.0
436	dcflux_fire_pft(:,:)=0.0
437	dcflux_trace(:,:)=0.0
438	dcflux_trace_pft(:,:,:)=0.0
439	fc_crown(:)=0.0
440	litter_consump_pft(:,:)=0.0
441	litter_consump(:)=0.0
442	emi_nondef(:,:)=0.0
443	co2_fire(:,:)=0.0
444	human_ign(:)=0.0
445	tree_kill_frac(:)=0.0
446	fc_1hr_carbon(:,:)=0.0
447	fc_10hr_carbon(:,:)=0.0
448	fc_100hr_carbon(:,:)=0.0
449	fc_1000hr_carbon(:,:)=0.0
450	human_suppression(:)=0.0
451	natural_vegfrac(:) = zero
452	veget_max(:,:) = veget_max_org(:,:) - lcc(:,:)
453	WHERE(veget_max(:,:) <= zero) veget_max(:,:)=zero
454
455	emideforest_fuel_1hr(:,:,:) = zero
456	emideforest_fuel_10hr(:,:,:) = zero
457	emideforest_fuel_100hr(:,:,:) = zero
458	emideforest_fuel_1000hr(:,:,:) = zero
459	emideforest_litter(:,:,:) = zero
460	emideforest_biomass(:,:,:) = zero
461	emideforest_all(:,:) = zero
462	bafrac_deforest(:,:) = zero
463
464	ratio_C3_livegrass(:) = 1.0
465	ratio_C4_livegrass(:) = zero
466	char_dens_fuel_ave(:) = zero
467
468	!we assume 0.03 of lightnning flashes input has the power to ignite a fire
469	lightn_efficiency(:)=0.03
470
471
472	! Assign emission factors for trace gas emissions resulting from acflux_fire_pft
473	! ef_CO2 defined in stomate_constants.f90; the unit of ef_trace(:) is in g/(g mass burned)
474	DO j=1,nvm
475	ef_trace(j,1)=ef_CO2(j)/1000.0 !CO2
476	ef_trace(j,2)=ef_CO(j)/1000.0 !CO
477	ef_trace(j,3)=ef_CH4(j)/1000.0 !CH4
478	ef_trace(j,4)=ef_VOC(j)/1000.0 !VOC
479	ef_trace(j,5)=ef_TPM(j)/1000.0 !TPM
480	ef_trace(j,6)=ef_NOx(j)/1000.0 !NOx
481	ENDDO
482
483	! this is done by assuming that in forced-BA simulations, the input BA data
484	! is at monthly time step.
485	IF (ANY(monthday.eq.date)) THEN
486	proc_ba(:)=observed_ba(:)
487	ELSE
488	proc_ba(:)=0.
489	ENDIF
490
491	!----------------------------------------------------------------
492	! Prepare multiple variables needed for fire simulation
493	!----------------------------------------------------------------
494
495	!! 0. Prepare some parameters for the simulation of fires.
496
497	!0.1 calculate the fraction of area occupied by natural PFTs.
498	DO j=2,nvm
499	IF (natural(j)) THEN
500	natural_vegfrac(:)=natural_vegfrac(:)+veget_max(:,j)
501	ENDIF
502	ENDDO
503	!calculate pixelarea(lat,area)
504	area_land(:) = (resolution(:,1)resolution(:,2)contfrac(:))/10000.0 !unit:hectar (10000m2)
505	area_natural_veg(:) = area_land(:)*natural_vegfrac(:) !unit:hectar (10000m2)
506
507	! Calculate total above-ground litter
508	natural_litter_ag_total(:)=0.0 !natural_litter_ag_total include litter from all natural PFTs.
509	DO j=1,nvm
510	IF (natural(j)) THEN
511	natural_litter_ag_total(:)=natural_litter_ag_total(:)+(litter(:,imetabolic,j,iabove)+ &
512	litter(:,istructural,j,iabove))*veget_max(:,j)
513	ENDIF
514	ENDDO
515
516	!0.2 Calculate the moisture of extinction weighted by the fuel presence
517	! of different natural PFTs (only dead fuel expressed by sum of litter)
518	! me: moisture of extinction; me is defined in constantes_mtc.f90
519	! as 'flammability threshold' and is PFT-specific
520	me_litterfuel(:)=0.0 !PFT weighted moisture of extinction.
521	DO j=1,nvm
522	IF (natural(j)) THEN
523	WHERE (natural_litter_ag_total(:).gt.min_stomate)
524	me_litterfuel(:)=me_litterfuel(:)+((litter(:,imetabolic,j,iabove)+litter(:,istructural,j,iabove)) * &
525	veget_max(:,j)/natural_litter_ag_total(:))*me(j)
526	ELSEWHERE
527	me_litterfuel(:)=0.0
528	ENDWHERE
529	ENDIF
530	ENDDO
531
532	! calculate tree and grass coverage respectively, here the grass should
533	! not include pasture, so if the pasture PFT is included, the `natural`
534	! variable for the pasture PFTs should be FALSE
535	fpc_tree_total(:)=0.0; fpc_grass_total(:)=0.0
536	DO j=1,nvm
537	IF (natural(j)) THEN
538	IF (is_tree(j)) THEN
539	fpc_tree_total(:)=fpc_tree_total(:)+veget_max(:,j)
540	ELSE
541	fpc_grass_total(:)=fpc_grass_total(:)+veget_max(:,j)
542	ENDIF
543	ENDIF
544	ENDDO
545
546	!0.3 Introduce reduction factor for wind speed with respect to forest or grass ground coverage.
547	! wind_speed is @tex $U_{forward}$ @endtex in Thonicke 2010.
548	! Note that the unit of input wind speed (wspeed_daily) should be @tex $m s^{-1}$ @endtex,
549	! so that ROS is in unit of @tex $m min^{-1}$ @endtex.
550	wind_speed(:)=(fpc_tree_total(:)wspeed_daily(:)60.0*0.4)+ &
551	(fpc_grass_total(:)wspeed_daily(:)60.0*0.6)
552	wind_forward(:)=3.281wind_speed(:) !wind_forward is used in Eq(A5) as 3.281Uforward
553
554	!0.4 Calculate the length-to-breadth ratio of the elliptical fire scar.
555	! the minimum value of lb is 1 and maximum value 8
556	! [UNCLEAR] Notice that for Eq.(12) and Eq.(13) in Thonicke 2010,
557	! 0.06*wind_speed is used here as U_{forward}. Compared with original Eq.(12),
558	! an extra 0.06 was multiplied. Not yet know why this 0.06 is used, it's
559	! in the original code given to me by Patricia.
560	WHERE (wind_speed(:).lt.16.67)
561	lb(:) = 1
562	ELSEWHERE (wind_speed(:).ge.16.67)
563	lb(:)=min(8.,fpc_tree_total(:)* &
564	(1.0+(8.729((1.0-(exp(-0.030.06wind_speed(:))))*2.155)))+ &
565	(fpc_grass_total(:)(1.1+((0.06wind_speed(:))**0.0464))))
566	ENDWHERE
567
568	!0.5 Calculate fuel characteristics, fuel amount, bulk density.
569	DO k=1,nlitt
570	DO j=1,nvm
571	IF (natural(j)) THEN
572	WHERE (fuel_1hr(:,j,k).gt.0.0)
573	fuel_1hr_total(:)=fuel_1hr_total(:) + (fuel_1hr(:,j,k)/0.45)*veget_max(:,j) !fuel_1hr devided by 0.45 to change into the unit of dry mass(ranther than C)
574	ENDWHERE
575
576	WHERE (fuel_10hr(:,j,k).gt.0.0)
577	fuel_10hr_total(:)=fuel_10hr_total(:) + (fuel_10hr(:,j,k)/0.45)*veget_max(:,j)
578	ENDWHERE
579
580	WHERE (fuel_100hr(:,j,k).gt.0.0)
581	fuel_100hr_total(:)=fuel_100hr_total(:) + (fuel_100hr(:,j,k)/0.45)*veget_max(:,j)
582	ENDWHERE
583
584	WHERE (fuel_1000hr(:,j,k).gt.0.0)
585	fuel_1000hr_total(:)=fuel_1000hr_total(:) + (fuel_1000hr(:,j,k)/0.45)*veget_max(:,j)
586	ENDWHERE
587
588	! k.eq.1 so that the grass leaf biomass would not be counted twice in the livegrass
589	! the PFT should be natural but grass(i.e., not tree or pasture),
590	! note livegrass is in unit of dry mass but not carbon!
591	IF (k.eq.1 .and. natural(j) .and. (.not. is_tree(j))) THEN
592	livegrass(:)=livegrass(:) + (biomass(:,j,ileaf)/0.45)*veget_max(:,j)
593	ENDIF !k
594	ENDIF!natural
595	ENDDO
596	ENDDO
597
598	! calculate dead fuel load on the groud, note that dead fuel load include
599	! 1_hr/10_hr/100_hr fuel class of natural PFTs, and 1000_hr fuel is not
600	! included in deaf_fuel and net_fuel
601	dead_fuel(:) = fuel_1hr_total(:) + fuel_10hr_total(:) + &
602	fuel_100hr_total(:) !total dead fuel g dry mass per m2
603	dead_fuel_all(:) = fuel_1hr_total(:) + fuel_10hr_total(:) + &
604	fuel_100hr_total(:)+fuel_1000hr_total(:)
605	dead_fuel_all_livegrass(:)=dead_fuel_all(:)+livegrass(:)
606	net_fuel(:)=(1.0-MINER_TOT)*(dead_fuel(:)/1000.0) ! in unit of kg dry matter per m2
607	! total net fuel accounting for mineral content, in kg dry mass per m2.
608	char_net_fuel(:) = net_fuel(:) +(1.0-MINER_TOT)*livegrass(:)/1000.0 ! in kg dry mass
609
610	!fuel bulk density, weighted per fuel class and fuel load
611	!dens_fuel is PFT-specific defined and initiated in constantes_mtc.f90, unit: kg/m3
612	dens_fuel_ave(:)=0.0
613	DO j=1,nvm
614	IF (natural(j)) THEN
615	WHERE (dead_fuel(:).gt.min_stomate)
616	ratio_fbd(:) = ( ( (fuel_1hr(:,j,istructural) + fuel_1hr(:,j,imetabolic)) &
617	+ fbd_a * (fuel_10hr(:,j,istructural)+ fuel_10hr(:,j,imetabolic)) &
618	+ fbd_b * (fuel_100hr(:,j,istructural) + fuel_100hr(:,j,imetabolic)) ) *veget_max(:,j) / 0.45) &
619	/ dead_fuel(:)
620	dens_fuel_ave(:) = dens_fuel_ave(:) + dens_fuel(j) * ratio_fbd(:) !kg/m3
621	ELSEWHERE
622	dens_fuel_ave(:)=0.0
623	ENDWHERE
624	ENDIF
625	ENDDO
626
627	!calculate livegrass fuel moisture
628	WHERE (livegrass(:).GT.min_stomate)
629	!Eq(B2); soilhum_daily is \omiga_{s} in Eq.(B2); dfm_lg is live grass moisture
630	dfm_lg(:) = max(0.0,((10.0/9.0)*soilhum_daily(:)-(1.0/9.0)))
631	ratio_C3_livegrass(:) = biomass(:,8,ileaf)/0.45*veget_max(:,8) / livegrass(:)
632	ratio_C4_livegrass(:) = biomass(:,9,ileaf)/0.45*veget_max(:,9) / livegrass(:)
633	ELSEWHERE
634	dfm_lg(:)=0.0
635	ENDWHERE
636
637	! influence of livegrass on FBD(fuel bulk density)
638	dens_livegrass_ave(:) = &
639	fbd_C3_livegrass * ratio_C3_livegrass(:) + &
640	fbd_C4_livegrass * ratio_C4_livegrass(:)
641
642	!calculate average fuel density weighted by different PFT and deal fuel & live grass
643	WHERE (dead_fuel(:).gt.min_stomate.or. livegrass(:).gt.min_stomate)
644	ratio_dead_fuel(:) = dead_fuel(:) / (dead_fuel(:) + livegrass(:))
645	ratio_live_fuel(:) = livegrass(:) / (dead_fuel(:) + livegrass(:))
646	char_dens_fuel_ave(:) = dens_fuel_ave(:)* ratio_dead_fuel(:) + &
647	dens_livegrass_ave(:) * ratio_live_fuel(:)
648	ENDWHERE
649
650	! calculate deal fuel and livegrass weighted average moisture of extinction.
651	! moist_extinction is me in Eq(8); moisture of extinction.
652	moist_extinction(:) = me_litterfuel(:) * ratio_dead_fuel(:) + &
653	me_livegrass * ratio_live_fuel(:)
654
655	!calculate mass weighted surface-area-to-volume ratio by fuel types.
656	WHERE (dead_fuel(:).gt.min_stomate)
657	sigma(:)=(fuel_1hr_total(:) * sigma_1hr + &
658	fuel_10hr_total(:) * sigma_10hr + &
659	fuel_100hr_total(:) * sigma_100hr) / dead_fuel(:)
660	ELSEWHERE
661	sigma(:)=0.00001
662	ENDWHERE
663
664	!-----------------------------------------------------------
665	! end preparing daily variables - start of simulation
666	!-----------------------------------------------------------
667
668
669	!! 1 Calculate daily fire burned area
670	!! 1.1 Calculate fire danger index
671	! output variables: dfm(\omiga_{0} in Eq.6); dfm_1hr(moisture content for 1hr fuel);
672	! d_fdi(fire danger index);
673	! modified variables: ni_acc(accumulative Nesterov Index)
674	CALL fire_danger_index(npts,d_fdi,dfm,dfm_1hr,t2m_min_daily, &
675	t2m_max_daily, precip_daily,me_litterfuel, ni_acc, &
676	fuel_1hr_total, fuel_10hr_total,fuel_100hr_total, &
677	dead_fuel,moist_extinction, ratio_dead_fuel,ratio_live_fuel)
678
679	!calculate human ignitions
680	CALL human_ign_fct(npts,popd,a_nd,human_ign)
681	!CALL human_suppression_func(npts,popd,human_suppression)
682
683	!apply the ignition efficiency as adjusted by fuel load
684	WHERE (dead_fuel_all_livegrass(:) .LE. fuel_low_bound)
685	ignition_efficiency(:)=0.
686	ELSEWHERE (dead_fuel_all_livegrass(:) .GT. fuel_low_bound .AND. dead_fuel_all_livegrass(:) .LT. fuel_high_bound)
687	ignition_efficiency(:)=(dead_fuel_all_livegrass(:) - fuel_low_bound)/(fuel_high_bound - fuel_low_bound)
688	ELSEWHERE
689	ignition_efficiency(:)=1.
690	ENDWHERE
691
692	!! 1.2 Calculate number of fires from lightning and human ignitions
693
694	!dervie potential fire number by lightnings
695	lightn_ign(:)=lightn(:)lightn_efficiency(:)(1-human_suppression(:))
696
697	!Here we apply Eq.(2), d_fdi: unitless; lightning ignition: 1/day/km**2 (as in the original nc file);
698	!human_ign: unit as 1/day/km**2
699	!0.01area_land: grid cell total land area in km*2;
700	!final unit: 1/day in the concered grid cell
701	WHERE (net_fuel(:).gt.0.001)
702	lightn_numfire(:)=d_fdi(:)ignition_efficiency(:) lightn_ign(:) * (0.01 * area_natural_veg(:))
703	human_numfire(:)=d_fdi(:)ignition_efficiency(:) human_ign(:) * (0.01 * area_natural_veg(:))
704	d_numfire(:)=lightn_numfire(:) + human_numfire(:)
705	ELSEWHERE !not enough fuel
706	d_numfire(:)=0.0
707	ENDWHERE
708
709	!! 1.3 Calculate daily area burnt
710
711	! calculate fire spread rate
712	!output variables: ros_f(ROS_{f,surface} in Eq.9); gamma(used in calculation of tau_l);
713	! wetness (dfm/me).
714	! note when call rate_of_spread, ros_f has replaced U_front in the definition of rate_of_spread SUBROUTINE
715	CALL rate_of_spread(npts,ros_f,wind_forward,sigma, &
716	dfm,me_litterfuel,H,char_dens_fuel_ave, &
717	char_net_fuel,moist_extinction,gamma,wetness)
718
719	! Calculate ROS_{b,surface}; Eq(10); wind_speed is U_{forward} in Eq.(10)
720	ros_b(:) = ros_f(:) * exp(-0.012 * wind_speed(:)) !Can FBP System; Eq.(10) in Thonicke et al. (2010)
721	WHERE (ros_b(:).lt. 0.05)
722	ros_b(:) = 0.0
723	ENDWHERE
724
725	! fire duration as a function of d_fdi; Eq(14)
726	WHERE (d_fdi(:).ge.min_stomate)
727	fire_durat(:)=241.0/(1.0+(240.0exp(-11.06d_fdi(:)))) !unit:min, Eq.(14) in Thonicke et al. (2010)
728	ELSEWHERE
729	fire_durat(:)=0.0
730	ENDWHERE
731
732	db(:)=ros_b(:)fire_durat(:) !unit: m/min min = m
733	df(:)=ros_f(:)*fire_durat(:)
734
735	! Here we call deforestation fire process
736	IF (allow_deforest_fire) THEN
737	! Calculate deforestation fire burned area fraction and emissions
738	CALL deforestation_fire_proxy(npts,d_fdi,lcc, &
739	deforest_biomass_remain, &
740	def_fuel_1hr_remain,def_fuel_10hr_remain, &
741	def_fuel_100hr_remain,def_fuel_1000hr_remain, &
742	bafrac_deforest,emideforest_fuel_1hr,emideforest_fuel_10hr, &
743	emideforest_fuel_100hr, emideforest_fuel_1000hr,emideforest_biomass)
744
745	! If the accumulated deforestation fire burned area is higher than two times
746	! of deforested area, we treat the surpassing part as natural fires.
747	DO ivm = 1,nvm
748	IF (is_tree(ivm)) THEN
749	DO ipts = 1,npts
750	IF (bafrac_deforest(ipts,ivm) > min_stomate) THEN
751	IF (bafrac_deforest_accu(ipts,ivm)+bafrac_deforest(ipts,ivm) > 3*lcc(ipts,ivm)) THEN
752	d_area_burnt(ipts) = d_area_burnt(ipts) + bafrac_deforest(ipts,ivm) * area_land(ipts)
753	ba_escape(ipts) = ba_escape(ipts) + bafrac_deforest(ipts,ivm) * area_land(ipts)
754	emideforest_fuel_1hr(ipts,ivm,:) = zero
755	emideforest_fuel_10hr(ipts,ivm,:) = zero
756	emideforest_fuel_100hr(ipts,ivm,:) = zero
757	emideforest_fuel_1000hr(ipts,ivm,:) = zero
758	emideforest_biomass(ipts,ivm,:) = zero
759	bafrac_deforest(ipts,ivm) = zero
760	ELSE
761	bafrac_deforest_accu(ipts,ivm) = bafrac_deforest_accu(ipts,ivm) + bafrac_deforest(ipts,ivm)
762	emideforest_litter(ipts,ivm,:) = emideforest_fuel_1hr(ipts,ivm,:) + emideforest_fuel_10hr(ipts,ivm,:) &
763	+ emideforest_fuel_100hr(ipts,ivm,:) + emideforest_fuel_1000hr(ipts,ivm,:)
764
765	emideforest_litter_accu(ipts,ivm,:) = emideforest_litter_accu(ipts,ivm,:) + emideforest_litter(ipts,ivm,:)
766	emideforest_biomass_accu(ipts,ivm,:) = emideforest_biomass_accu(ipts,ivm,:) + emideforest_biomass(ipts,ivm,:)
767	!update the remaining fuel
768	def_fuel_1hr_remain(ipts,ivm,:) = def_fuel_1hr_remain(ipts,ivm,:)-emideforest_fuel_1hr(ipts,ivm,:)
769	def_fuel_10hr_remain(ipts,ivm,:) = def_fuel_10hr_remain(ipts,ivm,:)-emideforest_fuel_10hr(ipts,ivm,:)
770	def_fuel_100hr_remain(ipts,ivm,:) = def_fuel_100hr_remain(ipts,ivm,:)-emideforest_fuel_100hr(ipts,ivm,:)
771	def_fuel_1000hr_remain(ipts,ivm,:) = def_fuel_1000hr_remain(ipts,ivm,:)-emideforest_fuel_1000hr(ipts,ivm,:)
772	deforest_biomass_remain(ipts,ivm,:) = deforest_biomass_remain(ipts,ivm,:)-emideforest_biomass(ipts,ivm,:)
773	deforest_litter_remain(ipts,:,ivm,iabove) = deforest_litter_remain(ipts,:,ivm,iabove)-emideforest_litter(ipts,ivm,:)
774	END IF
775	END IF
776	END DO
777	END IF
778	END DO ! loop of nvm
779
780	WHERE(veget_max(:,:)>min_stomate)
781	emideforest_all(:,:) = (SUM(emideforest_litter,DIM=3) + SUM(emideforest_biomass,DIM=3))/veget_max(:,:)/dt_days
782	ENDWHERE
783	END IF
784
785	! area burnt from caused by ignitions on "this (simulation)" day [sptifire on daily timestep]
786	WHERE (lb(:).gt.min_stomate)
787	mean_fire_size_original(:)=(pi/(4.0lb(:)))((df(:)+db(:))**2.0)/10000.0 !mean_fire_size with unit as ha
788	d_area_burnt(:)=d_numfire(:)* mean_fire_size_original(:) ! d_area_burnt unit in ha
789	mean_fire_size(:)=mean_fire_size_original(:)
790
791	!set the flag to 1 when forcing BA, otherwise set the flag as -1 or 0
792	WHERE (ratio_flag(:).gt.0.0)
793	d_area_burnt(:)=proc_ba(:)
794	ELSEWHERE
795	d_area_burnt(:)=d_area_burnt(:)*baadjust_ratio(:) !set the ratio alwasy to 1 when not forcing BA.
796	ENDWHERE
797
798	! Prevent burned area within the limit of vegetated area
799	WHERE (d_area_burnt(:).gt.area_natural_veg(:))
800	d_area_burnt(:)=area_natural_veg(:)
801	ENDWHERE
802
803	fire_frac(:) = d_area_burnt(:)/area_land(:)
804
805	ENDWHERE ! lb.gt.0
806
807	!! 1.4 Calculate daily fire emissions
808
809	fuel_consum(:)=0.0
810	CALL fuel_consumption(npts,fuel_consum,fire_frac, &
811	fuel_1hr,fuel_10hr,fuel_100hr,fuel_1000hr,livegrass,biomass, &
812	fuel_1hr_total,dfm,dfm_1hr,wetness,MINER_TOT,fc_1hr,fc_lg, &
813	fc_lf,fc_10hr,fc_100hr,fc_1000hr,cf_ave, &
814	moist_extinction,dfm_lg,veget_max, &
815	ratio_flag,cf_coarse,cf_fine,cf_lg,cf_1hr,cf_10hr,cf_100hr,cf_1000hr)
816
817	! calculate surface fire intensity
818	! Notice the fuel_consum doesn't include the 1000hr fuel
819	WHERE ((fire_frac(:).gt.min_stomate).and.(fuel_consum(:).gt.min_stomate))
820	d_i_surface(:)=H(fuel_consum(:)/1000.0/fire_frac(:))(ros_f(:)/60.0) !Eq(15); unit of d_i_surface: kW/m
821	ENDWHERE
822
823	! we suppress the simulated fires if their intensity is too small
824	WHERE (d_i_surface(:) .lt. surface_threshold)
825	lightn_numfire(:)=0.0
826	human_numfire(:)=0.0
827	d_numfire(:)=0.0
828	d_area_burnt(:)=0.0
829	fire_frac(:)=0.0
830	d_i_surface(:)=0.0
831	mean_fire_size(:)=0.0
832	cf_lg(:)=0.0
833	cf_1hr(:)=0.0
834	cf_ave(:)=0.0
835	cf_10hr(:)=0.0
836	cf_100hr(:)=0.0
837	cf_1000hr(:)=0.0
838	ELSEWHERE
839	fire_numday(:)=fire_numday(:)+1.
840	d_numfire(:)=d_numfire(:)/area_land(:) !change the unit of d_numfire to 1ha^{-1}day^{-1}
841	lightn_numfire(:)=lightn_numfire(:)/area_land(:) !change the unit of d_numfire to 1ha^{-1}day^{-1}
842	human_numfire(:)=human_numfire(:)/area_land(:) !change the unit of d_numfire to 1ha^{-1}day^{-1}
843	ENDWHERE
844
845	! Implement fire consequences, remove biomass being burned, mortality etc.
846	DO j=1,nvm
847	IF (natural(j)) THEN
848	DO k = 1, nlitt
849	! surface fire: update fuel load per dead fuel class and ground litter
850	WHERE (d_i_surface(:) .ge. surface_threshold)
851	fuel_1hr(:,j,k)=fuel_1hr(:,j,k) - fc_1hr(:,j,k)*0.45 ! convert back to gC/m2; fuel_1h(npts,nvm,nlitt) unit shoule be gC/m2
852	fuel_10hr(:,j,k)=fuel_10hr(:,j,k) - fc_10hr(:,j,k)*0.45
853	fuel_100hr(:,j,k)=fuel_100hr(:,j,k) - fc_100hr(:,j,k)*0.45
854	fuel_1000hr(:,j,k)=fuel_1000hr(:,j,k) - fc_1000hr(:,j,k) !always in gC/m2
855	litter(:,k,j,iabove)=litter(:,k,j,iabove)-(fc_1hr(:,j,k)+fc_10hr(:,j,k)+fc_100hr(:,j,k))*0.45-fc_1000hr(:,j,k)
856	ELSEWHERE
857	fc_1hr(:,j,k)=0.
858	fc_10hr(:,j,k)=0.
859	fc_100hr(:,j,k)=0.
860	fc_1000hr(:,j,k)=0.
861	ENDWHERE
862	ENDDO !nlitt
863
864	!include ground litter emission to litter_consump_pft
865	litter_consump_pft(:,j)=litter_consump_pft(:,j) &
866	+ (fc_1hr(:,j,imetabolic)+fc_10hr(:,j,imetabolic)+fc_100hr(:,j,imetabolic))*0.45 &
867	+ fc_1000hr(:,j,imetabolic) &
868	+ (fc_1hr(:,j,istructural)+fc_10hr(:,j,istructural)+fc_100hr(:,j,istructural))*0.45 &
869	+ fc_1000hr(:,j,istructural)
870
871	!! Apply fire effects for trees, when surface fire intensity is greater than 50 W/m
872
873	!! start of crown scorching
874	! crown_length defined in constantes_mtc.f90 indicating the fraction of crown length to tree height.
875	IF (is_tree(j)) THEN
876
877	!when I_{surface} < 50kWm-1, ignitions are extinguished
878	WHERE (d_i_surface(:).ge. surface_threshold .and. PFTpresent(:,j))
879	sh(:,j)=f_sh(j)(d_i_surface(:)*0.667) !scorch height per PFT Eq.(16);
880	!f_sh(j) defined in constantes_mtc.f90
881	ELSEWHERE
882	sh(:,j)=0.0
883	ENDWHERE
884
885	WHERE (gamma(:).gt.min_stomate)
886	tau_l(:,j)=2.0*(cf_ave(:)/gamma(:)) ! Equation(46) in SPITFIRE technical document.
887	ELSEWHERE
888	tau_l(:,j)=0.0
889	ENDWHERE
890
891	ck(:,j)=0.0
892	pm_ck(:,j)=0.0
893
894	WHERE (height(:,j).gt.min_stomate .and. crown_length(j).gt.min_stomate &
895	.and. sh(:,j).ge.(height(:,j) - (height(:,j) * crown_length(j))))
896	ck(:,j) = (sh(:,j) - height(:,j) + (height(:,j) * crown_length(j))) &
897	/(height(:,j) * crown_length(j))
898	ENDWHERE
899
900	WHERE (sh(:,j).GE.height(:,j))
901	ck(:,j) = 1.0
902	ENDWHERE
903
904	! post-fire mortality from crown scorching
905	WHERE (ck(:,j).gt.min_stomate)
906	pm_ck(:,j)=r_ck(j)(ck(:,j)*p_ck(j)) !Eq.(22); r_ck & p_ck defined in constantes_mtc.f90
907	ENDWHERE
908
909	! post-fire mortality from cambial damage
910	! calculate Bark Thickness (bt) and tau_c
911	dia(:) = (height(:,j)/pipe_tune2)*(1./pipe_tune3)100 !height is defined as height_presc in src_parameters/constantes_veg.f90; unit:m
912	! because in lpj_crown.f90: height(:,j) = pipe_tune2(dia(:)*pipe_tune3)
913	bt(:,j)=BTpar1(j) * dia(:)+BTpar2(j) !Eq(21)
914	tau_c(:,j)=2.9 * (bt(:,j)**2.0)
915
916	! tau_l/tau_c is the ratio of the residence time of the fire to the critical time for cambial damage.
917	WHERE (tau_c(:,j).ge.min_stomate)
918	WHERE ( (tau_l(:,j)/tau_c(:,j)).ge.2.0)
919	pm_tau(:,j)=1.0
920	ELSEWHERE ((tau_l(:,j)/tau_c(:,j)).gt.0.22)
921	pm_tau(:,j)= (0.563*(tau_l(:,j)/tau_c(:,j)))-0.125
922	ELSEWHERE ((tau_l(:,j)/tau_c(:,j)).le.0.22)
923	pm_tau(:,j)=0.0
924	ENDWHERE
925	ELSEWHERE
926	pm_tau(:,j)=0.0
927	ENDWHERE
928
929	! Calculate total post-fire mortality from crown scorching AND cambial kill; Eq(18)
930	postf_mort(:,j)=pm_tau(:,j)+pm_ck(:,j)-(pm_tau(:,j)*pm_ck(:,j))
931
932	! number of individuals affected by fire in grid cell
933	nind_fa(:,j)=fire_frac(:)*ind(:,j)
934	! the number of individuals killed by crown scorching is
935	! derived by multiplying mortality rate with number of individuals affected by fire.
936	nind_kill(:,j) = postf_mort(:,j) * nind_fa(:,j)
937
938
939	! Here explains how consequence of fire-killed trees are handled.
940	! the fraction of trees that are kill in fire is `postf_mort*fire_frac`,
941	! Of the killed trees, `ck` of them have been crown-scroached by fire
942	! and the 1-,10- and 100-h live fuel and 5% of the 1000-h fuel are combusted.
943	! the biomass of trees that are crown-scoarched but not consumed in fire
944	! would go to bm_to_litter. Then, the killed but not crown-scroached tree
945	! biomass (1-ck) go directly to bm_to_litter.
946
947	! In summary, all the biomass of killed tress should be removed from
948	! live biomass pool as they are either combusted or transferred to bm_to_litter.
949	! Of the fraction (postf_mort * fire_frac) that are killed in fire, of which
950	! --- ck (fraction) have been crown scorached, of which
951	! -- 1-,10- and 100-h live fule and 5% of the 1000-h fuel are combusted;
952	! -- 95% 1000-h fuel go to bm_to_litter;
953	! --- (1-ck) (fraction) go to bm_to_litter.
954
955	tree_kill_frac(:)=postf_mort(:,j)*fire_frac(:)
956	!disturb_crown is carbon emssions from crown scorching.
957	disturb_crown(:,j) = tree_kill_frac(:) * ck(:,j) * &
958	(biomass(:,j,ileaf) + biomass(:,j,ifruit) + &
959	0.0450 * biomass(:,j,isapabove) + 0.0450 * biomass(:,j,iheartabove) + &
960	0.0450 * biomass(:,j,icarbres) + &
961	0.0750 * biomass(:,j,isapabove) + 0.0750 * biomass(:,j,iheartabove) + &
962	0.0750 * biomass(:,j,icarbres) + &
963	0.2100 * biomass(:,j,isapabove) + 0.2100 * biomass(:,j,icarbres) + &
964	0.2100 * biomass(:,j,iheartabove) + &
965	0.6700 * 0.050 * biomass(:,j,icarbres) + 0.6700 * 0.050 * biomass(:,j,iheartabove) + &
966	0.6700 * 0.050 * biomass(:,j,isapabove) &
967	)
968
969	!fire flux to atmosphere, now biomass emissions from trees are included.
970	dcflux_fire_pft(:,j)=dcflux_fire_pft(:,j)+disturb_crown(:,j)
971	fc_crown(:)=fc_crown(:)+disturb_crown(:,j)*veget_max(:,j)
972
973	! biomass to litter transfer; (the killed but not crown scorching) AND (crown scorching but not combusted)
974	! 1-ck(:,j)---killed but not crown scorching;
975	! ck(:,j)(1-0.045-0.075-0.21-0.670.05)---crown scorching but not combusted.
976	bm_to_litter(:,j,ileaf) = bm_to_litter(:,j,ileaf) + &
977	(tree_kill_frac(:) * (1-ck(:,j)) * biomass(:,j,ileaf))
978	bm_to_litter(:,j,ifruit) = bm_to_litter(:,j,ifruit) + &
979	(tree_kill_frac(:) * (1-ck(:,j)) * biomass(:,j,ifruit))
980	bm_to_litter(:,j,isapabove) = bm_to_litter(:,j,isapabove) + &
981	(tree_kill_frac(:) * (1-ck(:,j)+ck(:,j)* &
982	(1-0.045-0.075-0.21-0.670.05)) biomass(:,j,isapabove))
983	bm_to_litter(:,j,icarbres) = bm_to_litter(:,j,icarbres) + &
984	(tree_kill_frac(:) * (1-ck(:,j)+ck(:,j)* &
985	(1-0.045-0.075-0.21-0.670.05)) biomass(:,j,icarbres))
986	bm_to_litter(:,j,iheartabove) = bm_to_litter(:,j,iheartabove) + &
987	(tree_kill_frac(:) * (1-ck(:,j)+ck(:,j)* &
988	(1-0.045-0.075-0.21-0.670.05)) biomass(:,j,iheartabove))
989	! move root/sapbelow/heartbelow biomass to belowground litter pool.
990	bm_to_litter(:,j,iroot) = bm_to_litter(:,j,iroot) + &
991	(tree_kill_frac(:) * biomass(:,j,iroot))
992	bm_to_litter(:,j,isapbelow) = bm_to_litter(:,j,isapbelow) + &
993	(tree_kill_frac(:) * biomass(:,j,isapbelow))
994	bm_to_litter(:,j,iheartbelow) = bm_to_litter(:,j,iheartbelow) + &
995	(tree_kill_frac(:) * biomass(:,j,iheartbelow))
996
997	! biomass update after combustion and litter transfer.
998	biomass(:,j,1) = biomass(:,j,1) - (tree_kill_frac(:) * biomass(:,j,1))
999	biomass(:,j,2) = biomass(:,j,2) - (tree_kill_frac(:) * biomass(:,j,2))
1000	biomass(:,j,3) = biomass(:,j,3) - (tree_kill_frac(:) * biomass(:,j,3))
1001	biomass(:,j,4) = biomass(:,j,4) - (tree_kill_frac(:) * biomass(:,j,4))
1002	biomass(:,j,5) = biomass(:,j,5) - (tree_kill_frac(:) * biomass(:,j,5))
1003	biomass(:,j,6) = biomass(:,j,6) - (tree_kill_frac(:) * biomass(:,j,6))
1004	biomass(:,j,7) = biomass(:,j,7) - (tree_kill_frac(:) * biomass(:,j,7))
1005	biomass(:,j,8) = biomass(:,j,8) - (tree_kill_frac(:) * biomass(:,j,8))
1006
1007
1008	! Update number of individuals surviving fire, ready for the next day.
1009	ind(:,j)= ind(:,j) - nind_kill(:,j)
1010
1011	ELSE !grass
1012
1013	WHERE (d_i_surface(:).ge.surface_threshold .and. PFTpresent(:,j))
1014	!grass leaf & fruit consumption are also included in litter consumption
1015	litter_consump_pft(:,j)=litter_consump_pft(:,j)+fc_lg(:,j)+fc_lf(:,j) !leaf & fruit consumption in fire for GRASS also included in litter consumption
1016	biomass(:,j,ileaf)=biomass(:,j,ileaf)-fc_lg(:,j)
1017	biomass(:,j,ifruit)=biomass(:,j,ifruit)-fc_lf(:,j)
1018	ENDWHERE
1019	ENDIF !IF (tree(j))
1020
1021	! add ground fuel combustion to the output variable of fire carbon flux per PFT
1022	! dcflux_fire_pft:total fire emissions, including crown fuels and ground litter consumption
1023	dcflux_fire_pft(:,j)=dcflux_fire_pft(:,j) + litter_consump_pft(:,j)
1024	litter_consump(:)=litter_consump(:)+(litter_consump_pft(:,j) *veget_max(:,j)/dt_days)
1025	ENDIF !IF natural
1026	ENDDO !fire effects per PFTj
1027
1028
1029	! total carbon emission this day; co2_fire should be equal to litter_consump+fc_crown
1030
1031	emi_nondef(:,:)= dcflux_fire_pft(:,:)/dt_days
1032	co2_fire(:,:)= emi_nondef(:,:) + emideforest_all(:,:)
1033	fc_crown(:)=fc_crown(:)/dt_days
1034
1035	!!1.6 Calculate trace gas emissions
1036	DO x=1,6 !trace_species
1037	DO j=1,nvm
1038	dcflux_trace_pft(:,j,x)=dcflux_fire_pft(:,j) * ef_trace(j,x)/0.45 !
1039	dcflux_trace(:,x)=dcflux_trace(:,x)+ (dcflux_trace_pft(:,j,x) / dt_days)veget_max(:,j) ! dcflux_trace has the same unit with dcflux_fire_pft, in g/m*2/day
1040	ENDDO
1041	ENDDO
1042
1043	DO k=1,nlitt
1044	fc_1hr_carbon(:,:)=fc_1hr_carbon(:,:)+fc_1hr(:,:,k)*0.45
1045	fc_10hr_carbon(:,:)=fc_10hr_carbon(:,:)+fc_10hr(:,:,k)*0.45
1046	fc_100hr_carbon(:,:)=fc_100hr_carbon(:,:)+fc_100hr(:,:,k)*0.45
1047	fc_1000hr_carbon(:,:)=fc_1000hr_carbon(:,:)+fc_1000hr(:,:,k)
1048	ENDDO
1049
1050	CALL histwrite_p (hist_id_stomate, 'CO2_FIRE_NonDef', itime, &
1051	emi_nondef(:,:), npts*nvm, horipft_index)
1052	CALL histwrite_p (hist_id_stomate, 'CO2_FIRE_Def', itime, &
1053	emideforest_all(:,:), npts*nvm, horipft_index)
1054	CALL histwrite_p (hist_id_stomate, 'D_FDI', itime, &
1055	d_fdi(:), npts, hori_index)
1056	CALL histwrite_p (hist_id_stomate, 'ROS_F', itime, &
1057	ros_f(:), npts, hori_index)
1058	CALL histwrite_p (hist_id_stomate, 'LIGHTN_NUMFIRE', itime, &
1059	lightn_numfire(:), npts, hori_index)
1060	CALL histwrite_p (hist_id_stomate, 'HUMAN_NUMFIRE', itime, &
1061	human_numfire(:), npts, hori_index)
1062	CALL histwrite_p (hist_id_stomate, 'D_NUMFIRE', itime, &
1063	d_numfire(:), npts, hori_index)
1064	CALL histwrite_p (hist_id_stomate, 'OBSERVED_BA', itime, &
1065	observed_ba(:), npts, hori_index)
1066	CALL histwrite_p (hist_id_stomate, 'D_AREA_BURNT', itime, &
1067	d_area_burnt(:), npts, hori_index)
1068	CALL histwrite_p (hist_id_stomate, 'BA_ESCAPE', itime, &
1069	ba_escape(:), npts, hori_index)
1070	CALL histwrite_p (hist_id_stomate, 'FIREFRAC_SPITFIRE', itime, &
1071	fire_frac(:), npts, hori_index)
1072	CALL histwrite_p (hist_id_stomate, 'CROWN_CONSUMP', itime, &
1073	fc_crown(:), npts, hori_index)
1074	CALL histwrite_p (hist_id_stomate, 'LITTER_CONSUMP', itime, &
1075	litter_consump(:), npts, hori_index)
1076	CALL histwrite_p (hist_id_stomate, 'LIGHTN_IGN_TOTAL', itime, &
1077	lightn(:), npts, hori_index)
1078	CALL histwrite_p (hist_id_stomate, 'TRACE_GAS_CO2', itime, &
1079	dcflux_trace(:,1), npts, hori_index)
1080	CALL histwrite_p (hist_id_stomate, 'TRACE_GAS_CO', itime, &
1081	dcflux_trace(:,2), npts, hori_index)
1082	CALL histwrite_p (hist_id_stomate, 'TRACE_GAS_CH4', itime, &
1083	dcflux_trace(:,3), npts, hori_index)
1084	CALL histwrite_p (hist_id_stomate, 'TRACE_GAS_VOC', itime, &
1085	dcflux_trace(:,4), npts, hori_index)
1086	CALL histwrite_p (hist_id_stomate, 'TRACE_GAS_TPM', itime, &
1087	dcflux_trace(:,5), npts, hori_index)
1088	CALL histwrite_p (hist_id_stomate, 'TRACE_GAS_NOx', itime, &
1089	dcflux_trace(:,6), npts, hori_index)
1090	CALL histwrite_p (hist_id_stomate, 'FIRE_NUMDAY', itime, &
1091	fire_numday(:), npts, hori_index)
1092
1093	CALL histwrite_p (hist_id_stomate, 'bafrac_deforest', itime, &
1094	bafrac_deforest(:,:), npts*nvm, horipft_index)
1095	CALL histwrite_p (hist_id_stomate, 'bafrac_deforest_accu', itime, &
1096	bafrac_deforest_accu(:,:), npts*nvm, horipft_index)
1097
1098	CALL histwrite_p (hist_id_stomate, 'EDlitMET', itime, &
1099	emideforest_litter(:,:,imetabolic), npts*nvm, horipft_index)
1100	CALL histwrite_p (hist_id_stomate, 'EDlitSTR', itime, &
1101	emideforest_litter(:,:,istructural), npts*nvm, horipft_index)
1102	CALL histwrite_p (hist_id_stomate, 'AccEDlitMET', itime, &
1103	emideforest_litter_accu(:,:,imetabolic), npts*nvm, horipft_index)
1104	CALL histwrite_p (hist_id_stomate, 'AccEDlitSTR', itime, &
1105	emideforest_litter_accu(:,:,istructural), npts*nvm, horipft_index)
1106
1107
1108
1109	CALL histwrite_p (hist_id_stomate, 'EDbioLEAF', itime, &
1110	emideforest_biomass(:,:,ileaf), npts*nvm, horipft_index)
1111	CALL histwrite_p (hist_id_stomate, 'EDbioRESERVE', itime, &
1112	emideforest_biomass(:,:,icarbres), npts*nvm, horipft_index)
1113	CALL histwrite_p (hist_id_stomate, 'EDbioFRUIT', itime, &
1114	emideforest_biomass(:,:,ifruit), npts*nvm, horipft_index)
1115	CALL histwrite_p (hist_id_stomate, 'EDbioSapABOVE', itime, &
1116	emideforest_biomass(:,:,isapabove), npts*nvm, horipft_index)
1117	CALL histwrite_p (hist_id_stomate, 'EDbioHeartABOVE', itime, &
1118	emideforest_biomass(:,:,iheartabove), npts*nvm, horipft_index)
1119	CALL histwrite_p (hist_id_stomate, 'EDbioSapBELOW', itime, &
1120	emideforest_biomass(:,:,isapbelow), npts*nvm, horipft_index)
1121	CALL histwrite_p (hist_id_stomate, 'EDbioHeartBELOW', itime, &
1122	emideforest_biomass(:,:,iheartbelow), npts*nvm, horipft_index)
1123	CALL histwrite_p (hist_id_stomate, 'EDbioROOT', itime, &
1124	emideforest_biomass(:,:,iroot), npts*nvm, horipft_index)
1125
1126	CALL histwrite_p (hist_id_stomate, 'AccEDbioLEAF', itime, &
1127	emideforest_biomass_accu(:,:,ileaf), npts*nvm, horipft_index)
1128	CALL histwrite_p (hist_id_stomate, 'AccEDbioRESERVE', itime, &
1129	emideforest_biomass_accu(:,:,icarbres), npts*nvm, horipft_index)
1130	CALL histwrite_p (hist_id_stomate, 'AccEDbioFRUIT', itime, &
1131	emideforest_biomass_accu(:,:,ifruit), npts*nvm, horipft_index)
1132	CALL histwrite_p (hist_id_stomate, 'AccEDbioSapABOVE', itime, &
1133	emideforest_biomass_accu(:,:,isapabove), npts*nvm, horipft_index)
1134	CALL histwrite_p (hist_id_stomate, 'AccEDbioHeartABOVE', itime, &
1135	emideforest_biomass_accu(:,:,iheartabove), npts*nvm, horipft_index)
1136	CALL histwrite_p (hist_id_stomate, 'AccEDbioSapBELOW', itime, &
1137	emideforest_biomass_accu(:,:,isapbelow), npts*nvm, horipft_index)
1138	CALL histwrite_p (hist_id_stomate, 'AccEDbioHeartBELOW', itime, &
1139	emideforest_biomass_accu(:,:,iheartbelow), npts*nvm, horipft_index)
1140	CALL histwrite_p (hist_id_stomate, 'AccEDbioROOT', itime, &
1141	emideforest_biomass_accu(:,:,iroot), npts*nvm, horipft_index)
1142
1143
1144	!debuging only
1145	!CALL histwrite_p (hist_id_stomate, 'cf_lg', itime, &
1146	! cf_lg(:), npts, hori_index)
1147	!CALL histwrite_p (hist_id_stomate, 'cf_1hr', itime, &
1148	! cf_1hr(:), npts, hori_index)
1149	!CALL histwrite_p (hist_id_stomate, 'cf_10hr', itime, &
1150	! cf_10hr(:), npts, hori_index)
1151	!CALL histwrite_p (hist_id_stomate, 'cf_100hr', itime, &
1152	! cf_100hr(:), npts, hori_index)
1153	!CALL histwrite_p (hist_id_stomate, 'cf_ave', itime, &
1154	! cf_ave(:), npts, horipft_index)
1155	!CALL histwrite_p (hist_id_stomate, 'cf_1000hr', itime, &
1156	! cf_1000hr(:), npts, hori_index)
1157	!CALL histwrite_p (hist_id_stomate, 'fc_1hr_carbon', itime, &
1158	! fc_1hr_carbon(:,:), npts*nvm, horipft_index)
1159	!CALL histwrite_p (hist_id_stomate, 'fc_10hr_carbon', itime, &
1160	! fc_10hr_carbon(:,:), npts*nvm, horipft_index)
1161	!CALL histwrite_p (hist_id_stomate, 'fc_100hr_carbon', itime, &
1162	! fc_100hr_carbon(:,:), npts*nvm, horipft_index)
1163	!CALL histwrite_p (hist_id_stomate, 'fc_1000hr_carbon', itime, &
1164	! fc_1000hr_carbon(:,:), npts*nvm, horipft_index)
1165	CALL histwrite_p (hist_id_stomate, 'fuel_1hr_met', itime, &
1166	fuel_1hr(:,:,imetabolic), npts*nvm, horipft_index)
1167	CALL histwrite_p (hist_id_stomate, 'fuel_1hr_str', itime, &
1168	fuel_1hr(:,:,istructural), npts*nvm, horipft_index)
1169	CALL histwrite_p (hist_id_stomate, 'fuel_10hr_met', itime, &
1170	fuel_10hr(:,:,imetabolic), npts*nvm, horipft_index)
1171	CALL histwrite_p (hist_id_stomate, 'fuel_10hr_str', itime, &
1172	fuel_10hr(:,:,istructural), npts*nvm, horipft_index)
1173	CALL histwrite_p (hist_id_stomate, 'fuel_100hr_met', itime, &
1174	fuel_100hr(:,:,imetabolic), npts*nvm, horipft_index)
1175	CALL histwrite_p (hist_id_stomate, 'fuel_100hr_str', itime, &
1176	fuel_100hr(:,:,istructural), npts*nvm, horipft_index)
1177	CALL histwrite_p (hist_id_stomate, 'fuel_1000hr_met', itime, &
1178	fuel_1000hr(:,:,imetabolic), npts*nvm, horipft_index)
1179	CALL histwrite_p (hist_id_stomate, 'fuel_1000hr_str', itime, &
1180	fuel_1000hr(:,:,istructural), npts*nvm, horipft_index)
1181	!CALL histwrite_p (hist_id_stomate, 'dead_fuel', itime, &
1182	! dead_fuel(:), npts, hori_index)
1183	!CALL histwrite_p (hist_id_stomate, 'dead_fuel_all', itime, &
1184	! dead_fuel_all(:), npts, hori_index)
1185	!CALL histwrite_p (hist_id_stomate, 'ck', itime, &
1186	! ck(:,:), npts*nvm, horipft_index)
1187	!CALL histwrite_p (hist_id_stomate, 'sh', itime, &
1188	! sh(:,:), npts*nvm, horipft_index)
1189	!CALL histwrite_p (hist_id_stomate, 'postf_mort', itime, &
1190	! postf_mort(:,:), npts*nvm, horipft_index)
1191	!CALL histwrite_p (hist_id_stomate, 'pm_tau', itime, &
1192	! pm_tau(:,:), npts*nvm, horipft_index)
1193	!CALL histwrite_p (hist_id_stomate, 'pm_ck', itime, &
1194	! pm_ck(:,:), npts*nvm, horipft_index)
1195
1196	!CALL histwrite_p (hist_id_stomate, 'topsoilhum_daily', itime, &
1197	! soilhum_daily(:), npts, horipft_index)
1198
1199	!CALL histwrite_p (hist_id_stomate, 'fire_durat', itime, &
1200	! fire_durat(:), npts, hori_index)
1201	!CALL histwrite_p (hist_id_stomate, 'ros_b', itime, &
1202	! ros_b(:), npts, hori_index)
1203	!CALL histwrite_p (hist_id_stomate, 'ros_f', itime, &
1204	! ros_f(:), npts, hori_index)
1205	!CALL histwrite_p (hist_id_stomate, 'wind_speed', itime, &
1206	! wind_speed(:), npts, hori_index)
1207	!CALL histwrite_p (hist_id_stomate, 'mean_fire_size_or', itime, &
1208	! mean_fire_size_original(:), npts, hori_index)
1209	!CALL histwrite_p (hist_id_stomate, 'mean_fire_size', itime, &
1210	! mean_fire_size(:), npts, hori_index)
1211	!CALL histwrite_p (hist_id_stomate, 'd_i_surface', itime, &
1212	! d_i_surface(:), npts, hori_index)
1213	!CALL histwrite_p (hist_id_stomate, 'sigma', itime, &
1214	! sigma(:), npts, hori_index)
1215	!CALL histwrite_p (hist_id_stomate, 'char_dens_fuel_ave', itime, &
1216	! char_dens_fuel_ave(:), npts, hori_index)
1217	CALL xios_orchidee_send_field("D_AREA_BURNT",d_area_burnt)
1218	CALL xios_orchidee_send_field("FIREFRAC_SPITFIRE",fire_frac)
1219
1220	END SUBROUTINE SPITFIRE
1221
1222
1223	!! ================================================================================================================================
1224	!! SUBROUTINE : fire_danger_index
1225	!!
1226	!>\BRIEF Calculation of the fire danger index using the Nesterov Index
1227	!!
1228	!! DESCRIPTION : None
1229	!!
1230	!! RECENT CHANGE(S): None
1231	!!
1232	!! RETURN VALUE : dfm, dfm_1hr, d_fdi
1233	!!
1234	!! MODIFIED VALUE: ni_acc
1235	!!
1236	!! REFERENCE(S) :
1237	!! - K.Thonicke, A.Spessa, I.C. Prentice, S.P.Harrison, L.Dong, and
1238	!! C.Carmona-Moreno, 2010: The influence of vegetation, fire spread and fire
1239	!! behaviour on biomass burning and trace gas emissions: results from a
1240	!! process-based model. Biogeosciences, 7, 1991-2011.
1241	!!
1242	!! FLOWCHART : None
1243	!! \n
1244	!_ ================================================================================================================================
1245
1246	SUBROUTINE fire_danger_index(npts,d_fdi,dfm,dfm_1hr,t2m_min_daily, &
1247	t2m_max_daily, precip_daily,me_litterfuel,ni_acc, &
1248	fuel_1hr_total,fuel_10hr_total,fuel_100hr_total, &
1249	dead_fuel,moist_extinction,ratio_dead_fuel, ratio_live_fuel)
1250
1251
1252	implicit none
1253
1254	! Domain size
1255	INTEGER, INTENT(in) :: npts
1256
1257	! INPUT VARIABLES
1258	REAL(r_std),DIMENSION(npts), INTENT(in) :: t2m_min_daily !! Daily minimum 2 meter temperatures (K)
1259	REAL(r_std),DIMENSION(npts), INTENT(in) :: t2m_max_daily !! Daily maximum 2 meter temperatures (K)
1260	REAL(r_std), DIMENSION(npts), INTENT(in) :: precip_daily !! daily precip (mm/day)
1261	REAL(r_std), DIMENSION(npts), INTENT(in) :: me_litterfuel !! the moisture of extinction weighted by the
1262	!! dead groud litter fuel of different natural PFTs, unitless
1263	REAL(r_std), DIMENSION(npts), INTENT(in) :: fuel_1hr_total !! in the unit of dry mass rather than C
1264	REAL(r_std), DIMENSION(npts), INTENT(in) :: fuel_10hr_total
1265	REAL(r_std), DIMENSION(npts), INTENT(in) :: fuel_100hr_total
1266	REAL(r_std), DIMENSION(npts), INTENT(in) :: dead_fuel
1267	REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_extinction
1268	REAL(r_std), DIMENSION(npts), INTENT(in) :: ratio_dead_fuel
1269	REAL(r_std), DIMENSION(npts), INTENT(in) :: ratio_live_fuel
1270
1271	! OUTPUT VARIABLES
1272	REAL(r_std), DIMENSION(npts), INTENT(out) :: dfm !! Daily fuel moisture (0,1), unitless.
1273	REAL(r_std), DIMENSION(npts), INTENT(out) :: dfm_1hr !! Fuel moisture for the 1hr fuel, used only in calculating the dfm_lg_d1hr
1274	REAL(r_std), DIMENSION(npts), INTENT(out) :: d_fdi
1275
1276	! MODIFIED VARIABLES
1277	REAL(r_std), DIMENSION(npts),INTENT(INOUT) :: ni_acc !!
1278
1279	! LOCAL VARIABLES
1280	! alpha values are in proportion to their surface-to-volume ratios and reflects their impact on the overall dryness of fuel.
1281	! lower suface-area-to-volume ratio leads to low alpha value, meaning that fuels needs a long time to be dry has less importance
1282	! in determinning overall fuel mositure.
1283	REAL(r_std), parameter :: alpha=0.0015 !! this variables is NOT USED.
1284	REAL(r_std), parameter :: alpha_1hr=0.001 !! from alpha_1hr --> alpha_1000hr: source SupInfo Table 4 and Thonicke 2010 P1996
1285	REAL(r_std), parameter :: alpha_10hr=0.00005424
1286	REAL(r_std), parameter :: alpha_100hr=0.00001485
1287	REAL(r_std), parameter :: alpha_1000hr = 0.000001 !! NOT USED
1288	REAL(r_std), parameter :: alpha_livegrass = 0.0005 !! SOURCE
1289	REAL(r_std), DIMENSION(npts) :: d_NI !Nesterov Index (.C^2)
1290	REAL(r_std), DIMENSION(npts) :: alpha_fuel
1291	REAL(r_std), DIMENSION(npts) :: char_alpha_fuel
1292
1293
1294
1295	! Initialise
1296	d_NI(:)=0.
1297	alpha_fuel(:)=0.
1298
1299	! calculate Nesterov Index Eq.(5); ZeroCelsius=273.15 in stomate_constant.f90
1300	WHERE ( (precip_daily(:) .LE. 3.) .AND. ((t2m_min_daily(:) -4.) .GE. ZeroCelsius) )
1301	d_NI(:) = (t2m_max_daily(:)-273.16)*( t2m_max_daily(:) - (t2m_min_daily(:)-4.) )
1302	ENDWHERE
1303
1304	WHERE (d_NI(:) .GT. 0.)
1305	ni_acc(:) = ni_acc(:) + d_NI(:)
1306	ELSEWHERE
1307	ni_acc(:)=0.
1308	ENDWHERE
1309
1310	! litter moisture index, weighted per dead fuel class and
1311
1312	WHERE (dead_fuel(:).gt.min_stomate)
1313	alpha_fuel(:)=(alpha_1hr*fuel_1hr_total(:)+ &
1314	alpha_10hr*fuel_10hr_total(:)+ &
1315	alpha_100hr*fuel_100hr_total(:))/dead_fuel(:)
1316	ELSEWHERE
1317	alpha_fuel(:)=0.00001
1318	ENDWHERE
1319
1320	dfm(:)=exp(-alpha_fuel(:)*ni_acc(:)) ! daily fuel moisture; \omega_{o} in Eq(6);
1321	! used for fire spread rate calculation. 1000hr fuel ignored.
1322	dfm_1hr(:) = exp(-alpha_1hr * ni_acc(:)) !moisture for the fuel_1hr.
1323
1324	! calculate Fire Danger Index d_fdi, Eq.(8)
1325	! moist_extinction acts as the 'me' in the FDI equation (Equation 8 in Thonicke 2010 Biogeoscience)
1326	WHERE (moist_extinction(:).le.min_stomate)
1327	d_fdi(:)=0.0
1328	ELSEWHERE
1329	d_fdi(:)=max(0.0,(1.0 - 1.0/moist_extinction(:)*dfm(:) ))
1330	ENDWHERE
1331
1332	!Fire danger index becomes zero when Nesterov index is lower than 0.
1333	WHERE (d_NI(:).le.0.)
1334	d_fdi(:)=0.0
1335	ENDWHERE
1336
1337
1338	! output for testing only
1339	!CALL histwrite_p (hist_id_stomate, 'moist_extinction', itime, &
1340	! moist_extinction(:), npts, horipft_index)
1341	!CALL histwrite_p (hist_id_stomate, 'alpha_fuel', itime, &
1342	! alpha_fuel(:), npts, hori_index)
1343	!CALL histwrite_p (hist_id_stomate, 'ni_acc', itime, &
1344	! ni_acc(:), npts, hori_index)
1345	!CALL histwrite_p (hist_id_stomate, 't2m_min_daily', itime, &
1346	! t2m_min_daily(:), npts, hori_index)
1347	!CALL histwrite_p (hist_id_stomate, 't2m_max_daily', itime, &
1348	! t2m_max_daily(:), npts, hori_index)
1349	!CALL histwrite_p (hist_id_stomate, 'precip_daily', itime, &
1350	! precip_daily(:), npts, hori_index)
1351	!CALL histwrite_p (hist_id_stomate, 'dfm_1hr', itime, &
1352	! dfm_1hr(:), npts, hori_index)
1353
1354	END SUBROUTINE fire_danger_index
1355
1356	!! ================================================================================================================================
1357	!! SUBROUTINE : human_ign_fct
1358	!!
1359	!>\BRIEF Calculation human ignitions
1360	!!
1361	!! DESCRIPTION : None
1362	!!
1363	!! RECENT CHANGE(S): None
1364	!!
1365	!! RETURN VALUE : Nonw
1366	!!
1367	!! REFERENCE(S) :
1368	!! - K.Thonicke, A.Spessa, I.C. Prentice, S.P.Harrison, L.Dong, and
1369	!! C.Carmona-Moreno, 2010: The influence of vegetation, fire spread and fire
1370	!! behaviour on biomass burning and trace gas emissions: results from a
1371	!! process-based model. Biogeosciences, 7, 1991-2011.
1372	!!
1373	!! FLOWCHART : None
1374	!! \n
1375	!_ ================================================================================================================================
1376
1377	! Algorithm in this function is very carefully checked (2012/05/31)
1378	SUBROUTINE human_ign_fct(npts,popd,a_nd,human_ign)
1379
1380	implicit none
1381
1382	! Domain size
1383	INTEGER, INTENT(in) :: npts
1384	! INPUT Variables
1385	REAL(r_std), DIMENSION(npts), INTENT(in) :: popd !individuals km^-2
1386	REAL(r_std), DIMENSION(npts), INTENT(in) :: a_nd !ignitions individual^-1 d^-1
1387	! OUTPUT
1388	REAL(r_std), DIMENSION(npts), INTENT(out) :: human_ign
1389
1390	WHERE(popd(:).ge.min_stomate)
1391	!compared with Eq(3) in Thonicke et al. (2010), we have to again divide
1392	!by 100 to change the unit from ignitions d^{-1}ha^{-1} to ignitions d^{-1}km^{-1}
1393	human_ign(:)=30.0(exp(-0.5(popd(:)*0.5)))a_nd(:)*popd(:)/10000.0 !unit: ignitions d^{-1}km^{-1}
1394	ELSEWHERE
1395	human_ign(:)=0.0
1396	ENDWHERE
1397
1398	END SUBROUTINE human_ign_fct
1399
1400
1401	!! ================================================================================================================================
1402	!! SUBROUTINE : human_ign_fct
1403	!!
1404	!>\BRIEF Calculation human ignitions
1405	!!
1406	!! DESCRIPTION : None
1407	!!
1408	!! RECENT CHANGE(S): None
1409	!!
1410	!! RETURN VALUE : Nonw
1411	!!
1412	!! REFERENCE(S) :
1413	!! Li et al. (2011) A process-based fire parameterization of intermediate
1414	!! complexity in a Dynamic Global Vegetation Model.
1415	!!
1416	!! FLOWCHART : None
1417	!! \n
1418	!_ ================================================================================================================================
1419
1420	SUBROUTINE human_suppression_func(npts,popd,human_suppression)
1421
1422	implicit none
1423
1424	!Domain size
1425	INTEGER, INTENT(in) :: npts
1426
1427	!INPUT Variables
1428	REAL(r_std), DIMENSION(npts), INTENT(in) :: popd !individuals km^-2
1429	!OUTPUT
1430	REAL(r_std), DIMENSION(npts), INTENT(out) :: human_suppression
1431
1432	WHERE(popd(:).ge.min_stomate)
1433	! source: Li et al. (2011) A process-based fire parameterization of intermediate complexity in a Dynamic Global
1434	! Vegetation Model. Page 2766
1435	human_suppression(:) = 0.99-0.98exp(-0.025popd(:))
1436	ELSEWHERE
1437	human_suppression(:)=0.0
1438	ENDWHERE
1439
1440	END SUBROUTINE human_suppression_func
1441
1442
1443	!! ================================================================================================================================
1444	!! SUBROUTINE : rate_of_spread
1445	!!
1446	!>\BRIEF Calculation fire spread rate, This is the application of Eq.(9) in
1447	!! Thonicke et al. (2010)
1448	!!
1449	!! DESCRIPTION : None
1450	!!
1451	!! RECENT CHANGE(S): None
1452	!!
1453	!! RETURN VALUE : Nonw
1454	!!
1455	!! REFERENCE(S) :
1456	!! - K.Thonicke, A.Spessa, I.C. Prentice, S.P.Harrison, L.Dong, and
1457	!! C.Carmona-Moreno, 2010: The influence of vegetation, fire spread and fire
1458	!! behaviour on biomass burning and trace gas emissions: results from a
1459	!! process-based model. Biogeosciences, 7, 1991-2011.
1460	!!
1461	!! FLOWCHART : None
1462	!! \n
1463	!_ ================================================================================================================================
1464
1465	SUBROUTINE rate_of_spread(npts,U_front,wind_forward,sigma, &
1466	dfm,me_litterfuel,H, char_dens_fuel_ave, &
1467	char_net_fuel, moist_extinction,gamma,wetness)
1468
1469	implicit none
1470
1471	! Domain size
1472	INTEGER, INTENT(in) :: npts
1473
1474	! INPUT Variables
1475	REAL(r_std), DIMENSION(npts), INTENT(in) :: wind_forward
1476	REAL(r_std), DIMENSION(npts), INTENT(in) :: char_dens_fuel_ave
1477	REAL(r_std), DIMENSION(npts), INTENT(in) :: sigma
1478	REAL(r_std), DIMENSION(npts), INTENT(in) :: char_net_fuel
1479	REAL(r_std), DIMENSION(npts), INTENT(in) :: dfm
1480	REAL(r_std), DIMENSION(npts), INTENT(in) :: me_litterfuel
1481	REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_extinction
1482	REAL(r_std), INTENT(in) :: H !! heat content of the fuel (18000kJ/kg)
1483
1484	! OUTPUT Variables
1485	REAL(r_std), DIMENSION(npts), INTENT(out) :: U_front !! Forward fire spread rate (m min^{-1}), ROS_{f,surface} in Thonicke et al. 2010
1486	REAL(r_std), DIMENSION(npts), INTENT(out) :: gamma
1487	REAL(r_std), DIMENSION(npts), INTENT(out) :: wetness !! wetness, unitless, the ratio of dead fuel moisture to the moisture of extinction
1488
1489	! LOCAL Variables
1490	REAL(r_std), DIMENSION(npts) :: dummy
1491	REAL(r_std), DIMENSION(npts) :: dummy2
1492	REAL(r_std), DIMENSION(npts) :: beta !! Packing ratio, unitless, the ratio of fuel_bulk_density(kg m^{-3}) against
1493	!! oven-dry particle_density (513kg m^{-3}).
1494	REAL(r_std), DIMENSION(npts) :: beta_op
1495	REAL(r_std), DIMENSION(npts) :: q_ig
1496	REAL(r_std), DIMENSION(npts) :: eps
1497	REAL(r_std), DIMENSION(npts) :: b
1498	REAL(r_std), DIMENSION(npts) :: c
1499	REAL(r_std), DIMENSION(npts) :: e
1500	REAL(r_std), DIMENSION(npts) :: phi_wind
1501	REAL(r_std), DIMENSION(npts) :: xi !!propagating flux ratio used in the calculation of ROS_{f,surface}
1502	REAL(r_std), DIMENSION(npts) :: a
1503	REAL(r_std), DIMENSION(npts) :: gamma_max
1504	REAL(r_std), DIMENSION(npts) :: gamma_aptr !!gamma'
1505	REAL(r_std), DIMENSION(npts) :: ir
1506	REAL(r_std), DIMENSION(npts) :: bet
1507	REAL(r_std), DIMENSION(npts) :: moist_damp
1508
1509
1510	! mineral dampening coefficient
1511	REAL(r_std), parameter :: MINER_DAMP=0.41739 !Thonicke 2010 Table A1 Row 7
1512	REAL(r_std), parameter :: part_dens=513.0 !Thonicke 2010 Table A1 Row 5
1513	INTEGER :: m
1514
1515
1516
1517	! Initialise
1518	xi(:)=0.0
1519	gamma_aptr=0.0
1520	dummy2(:)=0.0
1521
1522	beta(:) = char_dens_fuel_ave(:) / part_dens !packing ratio
1523	beta_op(:) = 0.200395(sigma(:)*(-0.8189)) !Eq.(A6)
1524	bet(:) = beta(:) / beta_op(:) !prepare for Eq.(A5)
1525
1526	! heat of pre-ignition
1527	q_ig(:)=581.0+2594.0*dfm(:)
1528	! effective heating number
1529	eps(:)=exp(-4.528/sigma(:))
1530	! influence of wind speed
1531	b(:)=0.15988(sigma(:)*0.54)
1532	c(:)=7.47(exp(-0.8711(sigma(:)**0.55)))
1533	e(:)=0.7515(exp(-0.01094sigma(:)))
1534
1535	WHERE(bet(:).GT.0.00001.AND.wind_forward(:).GT.min_stomate)
1536	phi_wind(:)=c(:)(wind_forward(:)b(:))(bet(:)**(-e(:))) !Eq.(A5)
1537	ENDWHERE
1538
1539	WHERE(bet(:).LE.0.00001)
1540	phi_wind(:)=c(:)(wind_forward(:)b(:))0.00001 !??UNIMPORTANT not sure this will pose a problem
1541	!but when bet --> 0+, bet**(-e) can be very big, which
1542	!is on the contrary with the simplied form.
1543	ENDWHERE
1544
1545	WHERE(wind_forward(:).LE.min_stomate)
1546	phi_wind(:)=zero
1547	ENDWHERE
1548
1549	! propagating flux
1550	WHERE (sigma(:).le.0.00001)
1551	xi(:)=0.0
1552	ELSEWHERE
1553	xi(:) = (exp((0.792 + 3.7597 * (sigma(:)*0.5)) & !Eq.(A2)
1554	(beta(:) + 0.1))) / (192 + 7.9095 * sigma(:))
1555	ENDWHERE
1556
1557	! reaction intensity
1558	a(:)=8.9033(sigma(:)*(-0.7913))
1559	WHERE (a(:).le.0.00001 .OR. bet(:).le.0.00001)
1560	dummy(:)=0.0
1561	ELSEWHERE
1562	dummy(:)=exp(a(:)*(1.0-bet(:)))
1563	ENDWHERE
1564	gamma_max(:)=1.0/(0.0591+2.926(sigma(:)*(-1.5))) !Thonicke 2010 Table A1 Row 2
1565
1566	WHERE(bet(:).GT.0.0001)
1567	gamma_aptr(:)=gamma_max(:)(bet(:)a(:))dummy(:)
1568	ELSEWHERE
1569	gamma_aptr(:)=gamma_max(:)(0.0001)dummy(:)
1570	ENDWHERE
1571
1572	! wn/me in Eq. row 6 in Table A1 in Thonicke 2010.
1573	WHERE (moist_extinction(:).gt.min_stomate)
1574	wetness(:)=dfm(:)/moist_extinction(:)
1575	ELSEWHERE
1576	wetness(:)=0.0
1577	ENDWHERE
1578
1579	WHERE (wetness(:).gt.min_stomate)
1580	moist_damp(:)=max(0.0,(1.0-(2.59*wetness(:))+ &
1581	(5.11(wetness(:)2.0))-(3.52(wetness(:)**3.0))))
1582	ELSEWHERE
1583	moist_damp(:)=0.0
1584	ENDWHERE
1585
1586	ir(:)=gamma_aptr(:) * char_net_fuel(:) * H * moist_damp(:) * MINER_DAMP !Eq.(A1)
1587
1588	!for use in calculating tau_l for postfire mortality.
1589	gamma(:)=gamma_aptr(:)moist_damp(:)MINER_DAMP
1590
1591	! reaction intensity end
1592
1593	WHERE ((char_dens_fuel_ave(:).le.min_stomate).or. &
1594	(eps(:).le.0.0).or.(q_ig(:).le.min_stomate))
1595	U_front(:)=0.0
1596	ELSEWHERE
1597	U_front(:)=(ir(:) * xi(:) * (1.0 + phi_wind(:))) / &
1598	(char_dens_fuel_ave(:) * eps(:) * q_ig(:))
1599	ENDWHERE
1600
1601	END SUBROUTINE rate_of_spread
1602
1603	!! ================================================================================================================================
1604	!! SUBROUTINE : fuel_consumption
1605	!!
1606	!>\BRIEF Calculation fire fuel consumption
1607	!!
1608	!! DESCRIPTION : None
1609	!!
1610	!! RECENT CHANGE(S): None
1611	!!
1612	!! RETURN VALUE : Nonw
1613	!!
1614	!! REFERENCE(S) :
1615	!! - K.Thonicke, A.Spessa, I.C. Prentice, S.P.Harrison, L.Dong, and
1616	!! C.Carmona-Moreno, 2010: The influence of vegetation, fire spread and fire
1617	!! behaviour on biomass burning and trace gas emissions: results from a
1618	!! process-based model. Biogeosciences, 7, 1991-2011.
1619	!!
1620	!! FLOWCHART : None
1621	!! \n
1622	!_ ================================================================================================================================
1623
1624	!fuel_consum--total fuel consumed excluding the 1000hr fuel; fc_1hr/10hr/100hr--fuel consumed; in the unit of g biomass/m2 (C/0.45); fc_1000hr (unit gC/m2)
1625	!fuel_consum: unit g mass/m**2 weighted by veget_max among PFTs, includes only 1hr/10hr/100hr dead fuel, used for calculation of d_i_surface
1626	SUBROUTINE fuel_consumption(npts,fuel_consum,fire_frac, &
1627	fuel_1hr,fuel_10hr,fuel_100hr,fuel_1000hr,livegrass,biomass, &
1628	fuel_1hr_total,dfm,dfm_1hr,wetness,MINER_TOT,fc_1hr,fc_lg, &
1629	fc_lf,fc_10hr,fc_100hr,fc_1000hr,cf_ave, &
1630	moist_extinction,dfm_lg,veget_max, &
1631	ratio_flag,cf_coarse,cf_fine,cf_lg,cf_1hr,cf_10hr,cf_100hr,cf_1000hr)
1632
1633
1634	implicit none
1635
1636	integer :: d,j,k
1637
1638	! Domain size
1639	INTEGER, INTENT(in) :: npts
1640
1641	REAL(r_std), DIMENSION(npts), INTENT(in) :: livegrass
1642	REAL(r_std), DIMENSION(npts), INTENT(in) :: fire_frac
1643	REAL(r_std), DIMENSION(npts), INTENT(in) :: dfm
1644	REAL(r_std), DIMENSION(npts), INTENT(in) :: dfm_1hr !! Daily fuel moisture for 1hr fuel, unitless, (0,1)
1645	REAL(r_std), DIMENSION(npts), INTENT(in) :: dfm_lg !! Daily live grass fuel moisture, unitless, (0,1)
1646	REAL(r_std), DIMENSION(npts), INTENT(in) :: wetness
1647	REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_extinction
1648	REAL(r_std), DIMENSION(npts), INTENT(in) :: fuel_1hr_total
1649	REAL(r_std), DIMENSION(npts), INTENT(in) :: cf_coarse !! Externally forced consumption fraction for coarse fuel, used
1650	!! on 1000hr fuel.
1651	REAL(r_std), DIMENSION(npts), INTENT(in) :: cf_fine !! Externally forced consumption fraction for fine fuel, applied
1652	!! on live grass fuel, 1hr/10hr/100hr fuel.
1653	REAL(r_std), DIMENSION(npts), INTENT(in) :: ratio_flag !! gridded value used to indicate where the external consumption
1654	!! fractions will be applied in the model, when ratio_flag is larger
1655	!! than 0, the external consumption fractions will be used.
1656
1657	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: fuel_1hr
1658	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: fuel_10hr
1659	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: fuel_100hr
1660	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: fuel_1000hr
1661	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: biomass
1662	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max
1663	REAL(r_std), INTENT(in) :: MINER_TOT
1664
1665	REAL(r_std), DIMENSION(npts), INTENT(inout) :: fuel_consum !! total dead ground fuel consumed excluding the 1000hr fuel, as
1666	!! weighted by presence of PFTs (g dry mass per m2)
1667	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(out) :: fc_1hr !! fc_*hr:fuel consumed, in the unit of biomass (C/0.45).
1668	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(out) :: fc_10hr
1669	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(out) :: fc_100hr
1670	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(out) :: fc_1000hr
1671	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: fc_lg !! livegrass leaf fuel cosumed (g C m^{-2})
1672	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: fc_lf !! livegrass fruit fuel consumed (gC m^{-2})
1673	REAL(r_std), DIMENSION(npts), INTENT(out) :: cf_ave !! the average consumption fraction of 1hr/10hr/100hr fuel, unitless, (0,1)
1674
1675	!LOCAL VARIABLES
1676	REAL(r_std), DIMENSION(npts) :: dfm_lg_d1hr !! daily fuel moisture for 1hr dead fuel and livegrass moisture combined
1677	REAL(r_std), DIMENSION(npts), INTENT(out) :: cf_1hr !! cf_*hr:consumption fraction, unitless, (0,1)
1678	REAL(r_std), DIMENSION(npts), INTENT(out) :: cf_lg
1679	REAL(r_std), DIMENSION(npts), INTENT(out) :: cf_10hr
1680	REAL(r_std), DIMENSION(npts), INTENT(out) :: cf_100hr
1681	REAL(r_std), DIMENSION(npts), INTENT(out) :: cf_1000hr
1682	REAL(r_std), DIMENSION(npts) :: wetness_lg
1683	REAL(r_std), DIMENSION(npts) :: wetness_1hr !! the wetness of 1hr-fuel in comparison with moisture of extinction
1684	!! is influenced by live grass fuel moisture
1685	REAL(r_std), parameter :: me_livegrass=0.2 !! moisture of extinction for live grass (unitless)
1686
1687
1688	! Initialise
1689	fuel_consum(:)=0.0
1690	dfm_lg_d1hr(:)=0.0
1691	wetness_lg(:)=0.0
1692	wetness_1hr(:)=0.0
1693
1694	!dfm_lg is the wlg in Eq.(B2) and Eq.(B3)
1695	WHERE (fuel_1hr_total(:).gt.min_stomate)
1696	dfm_lg_d1hr(:)=dfm_1hr(:)+(dfm_lg(:)*(livegrass(:)/fuel_1hr_total(:))) !Eq.(B3)
1697	ENDWHERE
1698
1699	WHERE (moist_extinction(:).gt.min_stomate)
1700	wetness_lg(:)=dfm_lg(:)/me_livegrass
1701	wetness_1hr(:)=dfm_lg_d1hr(:)/moist_extinction(:)
1702	ELSEWHERE
1703	wetness_lg(:)=1.0
1704	wetness_1hr(:)=1.0
1705	ENDWHERE
1706
1707	!Note that the fuel consumption fractions depend on only the moisture state and
1708	!thus are PFT independent.
1709	DO j=1,nvm
1710	IF (natural(j)) THEN
1711	! livegrass consumption
1712	IF (.not.is_tree(j)) THEN
1713	fc_lg(:,j)=0.0
1714	fc_lf(:,j)=0.0
1715	WHERE (wetness_lg(:).le.0.18)
1716	cf_lg(:)=1.0
1717	ELSEWHERE (wetness_lg(:).gt.0.18.and.wetness_lg(:).le.0.73)
1718	cf_lg(:)=1.11-0.62*wetness_lg(:)
1719	ELSEWHERE (wetness_lg(:).gt.0.73.and.wetness_lg(:).le.1.0)
1720	cf_lg(:)=2.45-2.45*wetness_lg(:)
1721	ELSEWHERE (wetness_lg(:).gt.1.0)
1722	cf_lg(:)=0.0
1723	ENDWHERE
1724
1725	WHERE (ratio_flag(:).gt.0.0)
1726	cf_lg(:)=cf_fine(:)
1727	ENDWHERE
1728
1729	fc_lg(:,j)=cf_lg(:) (1.0-MINER_TOT) biomass(:,j,ileaf)*fire_frac(:)
1730	fc_lf(:,j)=cf_lg(:) (1.0-MINER_TOT) biomass(:,j,ifruit)*fire_frac(:)
1731	ENDIF !not a tree
1732
1733	!dead fuel consumption
1734	DO k=1,nlitt
1735	!1hr fuel consumption
1736	fc_1hr(:,j,k)=0.0
1737	WHERE (fuel_1hr(:,j,k).gt.0.0)
1738
1739	WHERE (wetness_1hr(:).le.0.18)
1740	cf_1hr(:)=1.0
1741	ELSEWHERE (wetness_1hr(:).gt.0.18.and.wetness_1hr(:).le.0.73)
1742	cf_1hr(:)=1.11-0.62*wetness_1hr(:)
1743	ELSEWHERE (wetness_1hr(:).gt.0.73.and.wetness_1hr(:).le.1.0)
1744	cf_1hr(:)=2.45-2.45*wetness_1hr(:)
1745	ELSEWHERE (wetness_1hr(:).gt.1.0)
1746	cf_1hr(:)=0.0
1747	ENDWHERE
1748
1749	WHERE (cf_1hr(:) .gt. 0.9)
1750	cf_1hr(:) = 0.9
1751	ENDWHERE
1752
1753	WHERE (ratio_flag(:).gt.0.0)
1754	cf_1hr(:)=cf_fine(:)
1755	ENDWHERE
1756
1757	fc_1hr(:,j,k)=cf_1hr(:)(1.0-MINER_TOT)(fuel_1hr(:,j,k)/0.45)*fire_frac(:)
1758	ENDWHERE
1759
1760
1761	!10hr fuel consumption
1762	fc_10hr(:,j,k)=0.0
1763	WHERE (fuel_10hr(:,j,k).gt.0.0)
1764	WHERE (wetness(:).le.0.12)
1765	cf_10hr(:)=1.0
1766	ELSEWHERE (wetness(:).gt.0.12.and.wetness(:).le.0.51)
1767	cf_10hr(:)=1.09-0.72*wetness(:)
1768	ELSEWHERE (wetness(:).gt.0.51.and.wetness(:).le.1.0)
1769	cf_10hr(:)=1.47-1.47*wetness(:)
1770	ELSEWHERE (wetness(:).gt.1.0)
1771	cf_10hr(:)=0.0
1772	ENDWHERE
1773
1774	WHERE (cf_10hr(:) .gt. 0.9)
1775	cf_10hr(:) = 0.9
1776	ENDWHERE
1777
1778	WHERE (ratio_flag(:).gt.0.0)
1779	cf_10hr(:)=cf_fine(:)
1780	ENDWHERE
1781
1782	fc_10hr(:,j,k)=cf_10hr(:)* (1.0-MINER_TOT)(fuel_10hr(:,j,k)/0.45)fire_frac(:)
1783	ENDWHERE
1784
1785
1786	!100hr fuel consumption
1787	fc_100hr(:,j,k)=0.0
1788	WHERE (fuel_100hr(:,j,k).gt.0.0)
1789	WHERE (wetness(:).le.0.38)
1790	cf_100hr(:)=0.98-0.85*wetness(:)
1791	ELSEWHERE (wetness(:).gt.0.38 .and. wetness(:).le.1.0)
1792	cf_100hr(:)=1.06-1.06*wetness(:)
1793	ELSEWHERE (wetness(:).gt.1.0)
1794	cf_100hr(:)=0.0
1795	ENDWHERE
1796
1797	WHERE (cf_100hr(:) .gt. fire_max_cf_100hr(j))
1798	cf_100hr(:) = fire_max_cf_100hr(j)
1799	ENDWHERE
1800
1801	WHERE (ratio_flag(:).gt.0.0)
1802	cf_100hr(:)=cf_fine(:)
1803	ENDWHERE
1804
1805	fc_100hr(:,j,k)=cf_100hr(:)(1.0-MINER_TOT) (fuel_100hr(:,j,k)/0.45)*fire_frac(:)
1806	ENDWHERE
1807
1808	!calculating average fuel consumption fraction
1809	!the 1000hr fuel consumption fraction is not included in the average consumption fraction.
1810	cf_ave(:)=(cf_1hr(:)+cf_10hr(:)+cf_100hr(:))/3.0
1811
1812	!1000hr fuel consumption, not influencing rate of spread or surface fire intensity (Rothermel 1972)
1813	fc_1000hr(:,j,k)=0.0
1814	WHERE (fuel_1000hr(:,j,k).gt.0.0 .and. wetness(:).le.1.0)
1815	cf_1000hr(:) = -0.8*wetness(:)+0.8
1816
1817	WHERE (cf_1000hr(:) .gt. fire_max_cf_1000hr(j))
1818	cf_1000hr(:) = fire_max_cf_1000hr(j)
1819	ENDWHERE
1820
1821	WHERE (ratio_flag(:).gt.0.0)
1822	cf_1000hr(:)=cf_coarse(:)
1823	ENDWHERE
1824	fc_1000hr(:,j,k)=cf_1000hr(:) * (1.0-MINER_TOT)* fuel_1000hr(:,j,k)*fire_frac(:)
1825	ENDWHERE
1826
1827	!total fuel consumption (without 1000hr fuel) in g biomass per m**2
1828	!Used to calculate fire frontline intensity
1829	fuel_consum(:)=fuel_consum(:)+(fc_1hr(:,j,k)+fc_10hr(:,j,k)+fc_100hr(:,j,k))*veget_max(:,j)
1830	ENDDO !nlitt
1831	ENDIF !natural
1832	ENDDO !pft
1833
1834	! Only for debug purposes
1835	!CALL histwrite_p (hist_id_stomate, 'cf_fine', itime, &
1836	! cf_fine(:), npts, hori_index)
1837	!CALL histwrite_p (hist_id_stomate, 'cf_coarse', itime, &
1838	! cf_coarse(:), npts, hori_index)
1839	!CALL histwrite_p (hist_id_stomate, 'dfm_lg_d1hr', itime, &
1840	! dfm_lg_d1hr(:), npts, hori_index)
1841	!CALL histwrite_p (hist_id_stomate, 'dfm_lg', itime, &
1842	! dfm_lg(:), npts, hori_index)
1843	!CALL histwrite_p (hist_id_stomate, 'dfm', itime, &
1844	! dfm(:), npts, hori_index)
1845	!CALL histwrite_p (hist_id_stomate, 'wetness', itime, &
1846	! wetness(:), npts, hori_index)
1847	!CALL histwrite_p (hist_id_stomate, 'wetness_lg', itime, &
1848	! wetness_lg(:), npts, hori_index)
1849	!CALL histwrite_p (hist_id_stomate, 'wetness_1hr', itime, &
1850	! wetness_1hr(:), npts, hori_index)
1851
1852	END SUBROUTINE fuel_consumption
1853
1854
1855	!! ================================================================================================================================
1856	!! SUBROUTINE : deforestation_fire
1857	!!
1858	!>\BRIEF Simulate deforestation processes from land cover change.
1859	!!
1860	!! DESCRIPTION : None
1861	!!
1862	!! RECENT CHANGE(S): None
1863	!!
1864	!! RETURN VALUE : Nonw
1865	!!
1866	!! REFERENCE(S) :
1867	!!
1868	!! FLOWCHART : None
1869	!! \n
1870	!_ ================================================================================================================================
1871
1872	SUBROUTINE deforestation_fire(npts,d_fdi,lcc,biomass,fuel_1hr,&
1873	fuel_10hr,fuel_100hr,fuel_1000hr,bafrac_deforest,&
1874	emideforest_fuel_1hr,emideforest_fuel_10hr,emideforest_fuel_100hr,&
1875	emideforest_fuel_1000hr,emideforest_biomass)
1876
1877	IMPLICIT NONE
1878
1879	!! Input variables
1880	! Domain size
1881	INTEGER, INTENT(in) :: npts
1882	REAL(r_std), DIMENSION(npts) :: d_fdi !! daily fire danger index,climatic fire risk, unitless (0,1)
1883	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lcc !! gross forest cover loss
1884	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: biomass
1885	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: fuel_1hr !! [gC^{m-2}]
1886	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: fuel_10hr !! [gC^{m-2}]
1887	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: fuel_100hr !! [gC^{m-2}]
1888	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: fuel_1000hr !! [gC^{m-2}]
1889
1890	!! Output variables
1891	REAL(r_std), DIMENSION(npts,nvm),INTENT(inout) :: bafrac_deforest !! Deforestation fire burned fraction,unitless
1892	REAL(r_std), DIMENSION(npts,nvm,nlitt),INTENT(inout) :: emideforest_fuel_1hr !! fuel consumption for 1hr fuel (gCm^{-2})
1893	REAL(r_std), DIMENSION(npts,nvm,nlitt),INTENT(inout) :: emideforest_fuel_10hr !!
1894	REAL(r_std), DIMENSION(npts,nvm,nlitt),INTENT(inout) :: emideforest_fuel_100hr !!
1895	REAL(r_std), DIMENSION(npts,nvm,nlitt),INTENT(inout) :: emideforest_fuel_1000hr !!
1896	REAL(r_std), DIMENSION(npts,nvm,nparts),INTENT(inout) :: emideforest_biomass !! live biomass emissions from deforestation fires [gCm^{-2}]
1897
1898	!! Local variables
1899	REAL(r_std), DIMENSION(npts) :: deforest_factor !!Factor to scale the deforestation fire fraction, unitless
1900	REAL(r_std), DIMENSION(npts) :: vartemp
1901	REAL(r_std), DIMENSION(npts) :: cc_litter !! Fire consumption fraction for litter, leaf,fruit and reserve biomass
1902	REAL(r_std), DIMENSION(npts) :: cc_wood !! Fire consumption fraction for heart and sapwood
1903	INTEGER(i_std) :: ilit,j !! Indices (unitless)
1904
1905
1906	deforest_factor(:) = 0.017
1907	cc_litter(:) = 0.7 + 0.3*d_fdi(:)
1908	cc_wood(:) = 0.3 + d_fdi(:)*0.3
1909
1910	DO j=2,nvm
1911	IF (is_tree(j)) THEN
1912	WHERE (lcc(:,j) > min_stomate)
1913	!The equation of y=0.58*x-0.006 is derived by regressing GFED3 reported
1914	!deforestation fire burned fraction to the gross forest loss by Hansen 2010
1915	!for the closed canopy Amazonian forest
1916	vartemp(:) = MAX(0.0001,lcc(:,j)*0.58-0.006)
1917	bafrac_deforest(:,j) = deforest_factor(:) * vartemp(:) * d_fdi(:)
1918
1919	!We calculate fire emissions from litter and biomass
1920	emideforest_fuel_1hr(:,j,imetabolic) = cc_litter(:) * bafrac_deforest(:,j) * fuel_1hr(:,j,imetabolic)
1921	emideforest_fuel_1hr(:,j,istructural) = cc_litter(:) * bafrac_deforest(:,j) * fuel_1hr(:,j,istructural)
1922	emideforest_fuel_10hr(:,j,imetabolic) = cc_litter(:) * bafrac_deforest(:,j) * fuel_10hr(:,j,imetabolic)
1923	emideforest_fuel_10hr(:,j,istructural) = cc_litter(:) * bafrac_deforest(:,j) * fuel_10hr(:,j,istructural)
1924	emideforest_fuel_100hr(:,j,imetabolic) = cc_litter(:) * bafrac_deforest(:,j) * fuel_100hr(:,j,imetabolic)
1925	emideforest_fuel_100hr(:,j,istructural) = cc_litter(:) * bafrac_deforest(:,j) * fuel_100hr(:,j,istructural)
1926	emideforest_fuel_1000hr(:,j,imetabolic) = cc_litter(:) * bafrac_deforest(:,j) * fuel_1000hr(:,j,imetabolic)
1927	emideforest_fuel_1000hr(:,j,istructural) = cc_litter(:) * bafrac_deforest(:,j) * fuel_1000hr(:,j,istructural)
1928
1929	emideforest_biomass(:,j,icarbres) = cc_litter(:) * bafrac_deforest(:,j) * biomass(:,j,icarbres)
1930	emideforest_biomass(:,j,ifruit) = cc_litter(:) * bafrac_deforest(:,j) * biomass(:,j,ifruit)
1931	emideforest_biomass(:,j,ileaf) = cc_litter(:) * bafrac_deforest(:,j) * biomass(:,j,ileaf)
1932	emideforest_biomass(:,j,iheartabove) = cc_wood(:) * bafrac_deforest(:,j) * biomass(:,j,iheartabove)
1933	emideforest_biomass(:,j,isapabove) = cc_wood(:) * bafrac_deforest(:,j) * biomass(:,j,isapabove)
1934	ENDWHERE
1935	END IF
1936	END DO
1937
1938	END SUBROUTINE deforestation_fire
1939
1940	!! ================================================================================================================================
1941	!! SUBROUTINE : deforestation_fire_proxy
1942	!!
1943	!>\BRIEF Simulate deforestation processes from land cover change.
1944	!!
1945	!! DESCRIPTION : None
1946	!!
1947	!! RECENT CHANGE(S): None
1948	!!
1949	!! RETURN VALUE : Nonw
1950	!!
1951	!! REFERENCE(S) :
1952	!!
1953	!! FLOWCHART : None
1954	!! \n
1955	!_ ================================================================================================================================
1956
1957	SUBROUTINE deforestation_fire_proxy(npts,d_fdi,lcc, &
1958	deforest_biomass_remain, &
1959	def_fuel_1hr_remain,def_fuel_10hr_remain, &
1960	def_fuel_100hr_remain,def_fuel_1000hr_remain, &
1961	bafrac_deforest, emideforest_fuel_1hr,emideforest_fuel_10hr, &
1962	emideforest_fuel_100hr, emideforest_fuel_1000hr,emideforest_biomass)
1963
1964
1965	IMPLICIT NONE
1966
1967	!! Input variables
1968	! Domain size
1969	INTEGER, INTENT(in) :: npts
1970	REAL(r_std), DIMENSION(npts) :: d_fdi !! daily fire danger index,climatic fire risk, unitless (0,1)
1971	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lcc !! gross forest cover loss
1972	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: def_fuel_1hr_remain !! [gC^{m-2}]
1973	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: def_fuel_10hr_remain !! [gC^{m-2}]
1974	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: def_fuel_100hr_remain !! [gC^{m-2}]
1975	REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in) :: def_fuel_1000hr_remain !! [gC^{m-2}]
1976	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in) :: deforest_biomass_remain
1977
1978	!! Output variables
1979	REAL(r_std), DIMENSION(npts,nvm),INTENT(inout) :: bafrac_deforest !! Deforestation fire burned fraction,unitless
1980	REAL(r_std), DIMENSION(npts,nvm,nlitt),INTENT(inout) :: emideforest_fuel_1hr !! fuel consumption for 1hr fuel (gCm^{-2})
1981	REAL(r_std), DIMENSION(npts,nvm,nlitt),INTENT(inout) :: emideforest_fuel_10hr !!
1982	REAL(r_std), DIMENSION(npts,nvm,nlitt),INTENT(inout) :: emideforest_fuel_100hr !!
1983	REAL(r_std), DIMENSION(npts,nvm,nlitt),INTENT(inout) :: emideforest_fuel_1000hr !!
1984	REAL(r_std), DIMENSION(npts,nvm,nparts),INTENT(inout) :: emideforest_biomass !! live biomass emissions from deforestation fires [gCm^{-2}]
1985
1986	!! Local variables
1987	REAL(r_std), DIMENSION(npts,nvm) :: bafrac_net !! Deforestation fire burned fraction,unitless
1988	REAL(r_std), DIMENSION(npts) :: deforest_factor !!Factor to scale the deforestation fire fraction, unitless
1989	REAL(r_std), DIMENSION(npts) :: vartemp
1990	REAL(r_std), DIMENSION(npts) :: cc_litter !! Fire consumption fraction for litter, leaf,fruit and reserve biomass
1991	REAL(r_std), DIMENSION(npts) :: cc_wood !! Fire consumption fraction for heart and sapwood
1992	INTEGER(i_std) :: ilit,j !! Indices (unitless)
1993
1994
1995	! the 7.2 is tested. The test result could be found at
1996	! /home/cyue/Documents/ORCHIDEE/TEST_MICTv6bis/Deforestation_test/Deforestation_fire_annual_2001_2010_55e0e3e_T2.nc
1997	! with the notebook be found at: /home/cyue/Documents/ORCHIDEE/TEST_MICTv6bis/Test_tag_MICTv6bisT2
1998	! the readme file is at: /home/orchidee01/ychao/MICTv6/bin/readme.txt
1999	deforest_factor(:) = 0.017*7.2
2000	cc_litter(:) = 0.7 + 0.3*d_fdi(:)
2001	cc_wood(:) = 0.3 + d_fdi(:)*0.3
2002
2003	DO j=2,nvm
2004	IF (is_tree(j)) THEN
2005	WHERE (lcc(:,j) > min_stomate)
2006	!The equation of y=0.58*x-0.006 is derived by regressing GFED3 reported
2007	!deforestation fire burned fraction to the gross forest loss by Hansen 2010
2008	!for the closed canopy Amazonian forest
2009	vartemp(:) = MAX(0.0001,lcc(:,j)*0.58-0.006)
2010	bafrac_deforest(:,j) = deforest_factor(:) * vartemp(:) * d_fdi(:)
2011	bafrac_net(:,j)=bafrac_deforest(:,j)/lcc(:,j)
2012
2013	!We calculate fire emissions from litter and biomass
2014	!note here because cc_litter*bafrac_net < 1, the fuel will never be
2015	!exhausted.
2016	emideforest_fuel_1hr(:,j,imetabolic) = cc_litter(:) * bafrac_net(:,j) * def_fuel_1hr_remain(:,j,imetabolic)
2017	emideforest_fuel_1hr(:,j,istructural) = cc_litter(:) * bafrac_net(:,j) * def_fuel_1hr_remain(:,j,istructural)
2018	emideforest_fuel_10hr(:,j,imetabolic) = cc_litter(:) * bafrac_net(:,j) * def_fuel_10hr_remain(:,j,imetabolic)
2019	emideforest_fuel_10hr(:,j,istructural) = cc_litter(:) * bafrac_net(:,j) * def_fuel_10hr_remain(:,j,istructural)
2020	emideforest_fuel_100hr(:,j,imetabolic) = cc_litter(:) * bafrac_net(:,j) * def_fuel_100hr_remain(:,j,imetabolic)
2021	emideforest_fuel_100hr(:,j,istructural) = cc_litter(:) * bafrac_net(:,j) * def_fuel_100hr_remain(:,j,istructural)
2022	emideforest_fuel_1000hr(:,j,imetabolic) = cc_litter(:) * bafrac_net(:,j) * def_fuel_1000hr_remain(:,j,imetabolic)
2023	emideforest_fuel_1000hr(:,j,istructural) = cc_litter(:) * bafrac_net(:,j) * def_fuel_1000hr_remain(:,j,istructural)
2024
2025	emideforest_biomass(:,j,icarbres) = cc_litter(:) * bafrac_net(:,j) * deforest_biomass_remain(:,j,icarbres)
2026	emideforest_biomass(:,j,ifruit) = cc_litter(:) * bafrac_net(:,j) * deforest_biomass_remain(:,j,ifruit)
2027	emideforest_biomass(:,j,ileaf) = cc_litter(:) * bafrac_net(:,j) * deforest_biomass_remain(:,j,ileaf)
2028	emideforest_biomass(:,j,iheartabove) = cc_wood(:) * bafrac_net(:,j) * deforest_biomass_remain(:,j,iheartabove)
2029	emideforest_biomass(:,j,isapabove) = cc_wood(:) * bafrac_net(:,j) * deforest_biomass_remain(:,j,isapabove)
2030	ENDWHERE
2031	END IF
2032	END DO
2033
2034	END SUBROUTINE deforestation_fire_proxy
2035	END MODULE lpj_spitfire

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

Download in other formats: