CLASSIC
Canadian Land Surface Scheme including Biogeochemical Cycles
disturbance_scheme Module Reference

Calculates disturbance as both natural and human-influenced fires. More...

Functions/Subroutines

subroutine, public disturb (thliq, THLW, THFC, uwind, useTracer, vwind, lightng, fcancmx, isand, rmatctem, ilg, il1, il2, sort, grclarea, thice, popdin, lucemcom, dofire, currlat, iday, fsnow, stemmass, rootmass, gleafmas, bleafmas, litrmass, tracerStemMass, tracerRootMass, tracerGLeafMass, tracerBLeafMass, tracerLitrMass, stemltdt, rootltdt, glfltrdt, blfltrdt, glcaemls, rtcaemls, stcaemls, blcaemls, ltrcemls, burnfrac, pstemmass, pgleafmass, emit_co2, emit_ch4, emit_co, emit_nmhc, emit_h2, emit_nox, emit_n2o, emit_pm25, emit_tpm, emit_tc, emit_oc, emit_bc, burnvegf, bterm_veg, mterm_veg, lterm, smfunc_veg)
 Calculates whether fire occurs, burned area, amount of C emitted and litter generated. More...
 
subroutine, public burntobare (il1, il2, nilg, sort, pvgbioms, pgavltms, pgavscms, burnvegf, pstemmass, pgleafmass, useTracer, fcancmx, stemmass, rootmass, gleafmas, bleafmas, litrmass, soilcmas, nppveg, tracerLitrMass, tracerSoilCMass, tracerGLeafMass, tracerBLeafMass, tracerStemMass, tracerRootMass)
 Update fractional coverages of pfts to take into account the area burnt by fire. Adjust all pools with new densities in their new areas and increase bare fraction. And while we are doing this also run a small check to make sure grid averaged quantities do not get messed up. More...
 

Detailed Description

Calculates disturbance as both natural and human-influenced fires.

CTEM represents disturbance as both natural and human-influenced fires. The original fire parametrization corresponding to CTEM v. 1.0 is described in Arora and Boer (2005) [5]. The parametrization has since been adapted and used in several other DGVMs (Kloster et al. 2010 [49], Kloster et al. 2012 [50], Migliavacca et al. 2013 [69], and Li et al. 2012 [61]). CTEM v. 2.0 incorporates changes suggested in these studies as well as several new improvements.

Fire in CTEM is simulated using a process-based scheme of intermediate complexity that accounts for all elements of the fire triangle: fuel load, combustibility of fuel, and an ignition source. CTEM represents the probability of a fire occurrence ( \(P_\mathrm{f}\)), for a representative area of \(500\, km^2\) ( \(a_{rep}\)), as

\[ \label{fieya} P_\mathrm{f} = P_\mathrm{b}P_\mathrm{i}P_\mathrm{m}, \qquad (Eqn 1)\]

where the right hand side terms represent the fire probabilities that are conditioned on (i) the availability of biomass as a fuel source ( \(P_\mathrm{b}\)), (ii) the combustibility of the fuel based on its moisture content ( \(P_\mathrm{m}\)), and (iii) the presence of an ignition source ( \(P_\mathrm{i}\)). The probability of fire and the subsequent calculations are performed for each PFT present in a grid cell (but the PFT index \(\alpha\) is omitted for clarity in Eqn 1). Since the CTEM parametrization is based on one fire per day per representative area, the representative area has to be sufficiently small that the requirement of only one fire per day is reasonable, yet sufficiently large such that it is not possible to burn the entire representative area in 1 day. Based on MODIS observed fire counts in Fig. 1 of Li et al. (2012) [61], \(500\, km^2\) is an appropriate size to not have more than one fire per day and still be a large enough area to be assumed representative of the grid cell as a whole.

The \(P_\mathrm{b}\) term depends on the aboveground biomass ( \(B_{ag}\)) available for sustaining a fire (which includes the green and brown leaf mass, stem mass and litter mass, \(B_{ag} = C_\mathrm{L} + C_\mathrm{S} + C_\mathrm{D}\)). Below a lower threshold of aboveground biomass ( \(B_{low}\); \(0.2\, kg\, C\, m^{-2}\) similar to Moorcroft et al. (2001) [71], and Kucharik et al. (2000) [53]), fire is not sustained and thus has a probability of 0. Above a biomass of \(1.0\, kg\, C\, m^{-2}\) ( \(B_{high}\)), \(P_\mathrm{b}\) is set to 1 as the amount of fuel available is assumed sufficient for fire. \(P_\mathrm{b}\) is then calculated using the aboveground biomass, \(B_{ag}\) ( \(kg\, C\, m^{-2}\)) with a linear variation between the upper and lower thresholds as

\[ \label{eqn:Pb} P_\mathrm{b}=\max\left[0, \min\left(1, \frac{B_{ag}-B_{low}} {B_{high} - B_{low}}\right)\right]. \qquad (Eqn 2)\]

The linear decrease of \(P_\mathrm{b}\) from \(B_{high}\) to \(B_{low}\) reflects the fragmentation of fuel that occurs as biomass decreases. Fuel fragmentation impacts upon area burned as it impedes the fire spread rate (Guyette et al. 2002) [43].

The probability of fire based on the presence of ignition sources ( \(P_\mathrm{i}\)) is influenced by both natural (lightning) and anthropogenic agents (either intentional or accidental). An initial lightning scalar, \(\vartheta_F\), that varies between 0 and 1 is found as

\[ \vartheta_F = \max\left[0, \min \left(1, \frac{F_c2g - F_{low}}{F_{high} - F_{low}} \right)\right], \qquad (Eqn 3)\]

where \(F_{low}\) and \(F_{high}\) represent lower and upper thresholds of cloud-to-ground lightning strikes ( \(F_c2g\), \(flashes\, km^{-2}\, month^{-1}\)) , respectively. Similar to Eqn 2, below the lower threshold ( \(F_{low}\); \(0.25\, flashes\, km^{-2}\, month^{-1}\)), \(\vartheta_F\) is 0 implying lightning strikes are not sufficient to cause fire ignition, above the upper threshold ( \(F_{high}\); \(10.0\, flashes\, km^{-2}\, month^{-1}\)) \(\vartheta_F\) is 1, as ignition sources now do not pose a constraint on fire. The amount of cloud-to-ground lightning, \(F_c2g\), is a fraction of the total lightning based on the relationship derived by Price and Rind (1993) [77] (approximation of their Eqs. 1 and 2) as

\[ F_c2g = 0.22 \exp (0.0059 \times \vert {\Phi}\vert) F_{tot}, \qquad (Eqn 4)\]

where \(\Phi\) is the grid cell latitude in degrees and \(F_{tot}\) is the total number of lightning \(flashes\, km^{-2}\, month^{-1}\) (both cloud-to-cloud and cloud-to-ground). The probability of fire due to natural ignition, \(P_i, n\), depends on the lightning scalar, \(\vartheta_F\), as

\[ P_i, n = y(\vartheta_F) - y(0)(1 - \vartheta_F) + \vartheta_F[1-y(1)] \nonumber\\ y(\vartheta_F) = \frac{1}{1 + \exp\left(\frac{0.8 - \vartheta_F}{0.1}\right)}. \qquad (Eqn 5)\]

Fire probability due to ignition caused by humans, \(P_i, h\), is parametrized following Kloster et al. (2010) [49] with a dependence on population density, \(p_\mathrm{d}\) ( \(number of people\, km^{-2}\))

\[ \label{eqn:Ph} P_i, h = \min\left[1, \left(\frac{p_\mathrm{d}}{p_{thres}}\right)^{0.43}\right], \qquad (Eqn 6)\]

where \(p_{thres}\) is a population threshold ( \(300\, people\, km^{-2}\)) above which \(P_{i, h}\) is 1. The probability of fire conditioned on ignition, \(P_\mathrm{i}\), is then the total contribution from both natural and human ignition sources

\[ \label{eqn:Pi} P_\mathrm{i} = \max[0, \min\{1, P_{i, n} + (1 - P_{i, n})P_{i, h}\}].\qquad (Eqn 7) \]

The population data used to calculate probability of fire ignition caused by humans and anthropogenic fire suppression (discussed further down in this section) is read in as a model geophysical field.

The probability of fire due to the combustibility of the fuel, \(P_\mathrm{m}\), is dependent on the soil moisture in vegetation's root zone and in the litter layer. The root-zone soil wetness ( \(\phi_{root}\), allocate.f90 Eqn. 1) is used as a surrogate for the vegetation moisture content and the soil wetness of the top soil layer as a surrogate for the litter moisture content. If a grid cell is covered by snow, \(P_\mathrm{m}\) is set to zero. The probability of fire conditioned on soil wetness in vegetation's rooting zone, \(P_{m, V}\), is then

\[ P_{m, V} = 1-\tanh \left[\left( \frac{1.75\ \phi_{root}} {E_\mathrm{V}}\right)^2\right], \qquad (Eqn 8)\]

where \(E_\mathrm{V}\) is the extinction soil wetness above which \(P_{f, V}\) is reduced to near zero and is set to 0.30.

The probability of fire based on the moisture content in the \(\textit{duff}\) layer, \(P_{m, D}\), which includes the brown leaf mass (grasses only) and litter mass ( \(B_{duff} = C_{L, b} + C_\mathrm{D}\); \(kg\, C\, m^{-2}\)), is calculated in a similar way but uses the soil wetness of the first soil layer, ( \(\phi_1\), photosynCanopyConduct.f90 Eqn. 7), as a surrogate for the moisture in the duff layer itself as

\[ P_{m, D} = 1 -\tanh\left[\left(\frac{1.75 \phi_1}{E_{\mathrm{D}}}\right)^2\right], \qquad (Eqn 9)\]

where the extinction soil wetness for the litter layer, \(E_{\mathrm{D}}\), is set to 0.50, which yields a higher probability of fire for the litter layer than for the vegetation for the same soil wetness. \(P_\mathrm{m}\) is then the weighted average of \(P_{m, V}\) and \(P_{m, D}\) given by

\[ \label{eqn:Pf} P_\mathrm{m} = P_{m, V} (1-f_{duff}) + P_{m, D} f_{duff} \nonumber \\ f_{duff}=\frac{B_{duff}}{B_{ag}}\qquad (Eqn 10)\]

where \(f_{duff}\) is the duff fraction of aboveground combustible biomass.

The area burned ( \(a\)) is assumed to be elliptical in shape for fires based upon the wind speed and properties of an ellipse

\[ a(t)=\pi \frac{l}{2}\frac{w}{2}= \frac{\pi}{2} (v_\mathrm{d}+v_\mathrm{u})v_\mathrm{p}t^2, \qquad (Eqn 11)\]

where \(l\) ( \(m\)) and \(w\) ( \(m\)) are the lengths of major and minor axes of the elliptical area burnt; \(v_\mathrm{d}\) ( \(km\, h^{-1}\)) and \(v_\mathrm{u}\) ( \(km\, h^{-1}\)) are the fire spread rates in the downwind and upwind directions, respectively; \(v_\mathrm{p}\) ( \(km\, h^{-1}\)) is the fire spread rate perpendicular to the wind direction and \(t\) is the time ( \(h\)).

The fire spread rate in the downwind direction ( \(v_\mathrm{d}\)) is represented as

\[ \label{firespreadrate} v_\mathrm{d} = v_{d, max}\, g(u)\, h(\phi_{r, d})\qquad (Eqn 12)\]

where \(v_{d, max}\) ( \(km\, h^{-1}\)) is the PFT-specific maximum fire spread rate from Li et al. (2012)[61], which is set to zero for crop PFTs (see also classicParams.f90). The functions \(g(u)\) accounts for the effect of wind speed and \( h(\phi_{r, d})\) accounts for the effect of rooting zone and duff soil wetness on the fire spread rate, as discussed below.

The wind speed ( \(u\); \(km\, h^{-1}\)) is used to determine the length ( \(l\)) to breadth ( \(w\)) ratio, \(L_\mathrm{b}\), of the elliptical area burned by fire

\[ \label{lb} L_\mathrm{b}= \frac{l}{w} = \frac{v_\mathrm{d} + v_\mathrm{u}}{2v_\mathrm{p}} = 1 + 10 [1 -\exp(-0.06 u)] \qquad (Eqn 13)\]

and its head to back ratio, \(H_\mathrm{b}\), following Li et al. (2012) [61], as

\[ \label{hb} H_\mathrm{b} = \frac{v_\mathrm{d}}{v_\mathrm{u}} = \frac{L_\mathrm{b} + (L_\mathrm{b}^2 - 1)^{0.5}}{L_\mathrm{b} - (L_\mathrm{b}^2 - 1)^{0.5}}, \qquad (Eqn 14)\]

which help determine the fire spread rate in the direction perpendicular to the wind speed and in the downward direction. Equations 13 and 14 are combined to estimate the wind scalar \(g(u)\) as

\[ g(u)= g(0) \frac{2.0 L_\mathrm{b}}{(1 + 1/H_\mathrm{b})} \nonumber\\ \frac{g(u)}{g(0)}=\frac{v_\mathrm{d}}{v_\mathrm{p}} = \frac{2.0 L_\mathrm{b}} {(1 + 1/H_\mathrm{b})}, \qquad (Eqn 15) \]

which varies between 0.05 and 1. The lower limit is imposed by the \(g(0)\) term, which has a value of 0.05 and represents the fire spread rate in the absence of wind ( \(u = 0\)); the upper limit is assigned a maximum value of 1. The fire spread rate in the absence of wind is essentially the spread rate in the direction perpendicular to the wind speed ( \(v_\mathrm{p}\)). The value of the \(g(0)\) term is derived by considering the case where the wind speed becomes very large. As \(u\) \(\rightarrow \infty\) then \(L_\mathrm{b} \rightarrow 11\) and \(H_\mathrm{b} \rightarrow 482\), while \(g(\infty)=1\) due to its definition, which yields \(g(0) = 0.0455 \approx 0.05\).

The dependence of fire spread rate on the rooting zone and duff soil wetness, \(h(\phi_{r, d})\) is represented as

\[ h(\phi_{r, d})= h(\phi_{root})(1-f_{duff}) + h(\phi_{1})f_{duff}\nonumber \\ h(\phi_{root})= \left(1-min \left(1, \frac{\phi_{root}}{E_\mathrm{V}} \right) \right)^2\nonumber \\ h(\phi_{1})= \left(1-min \left(1, \frac{\phi_{1}}{E_\mathrm{D}} \right) \right)^2. \qquad (Eqn 16)\]

Both \(h(\phi_{root})\) and \(h(\phi_{1})\) gradually decrease from 1 (when soil wetness is 0 and soil moisture does not constrain fire spread rate) to 0 when soil wetness exceeds the respective extinction wetness thresholds, \(E_\mathrm{V}\) and \(E_\mathrm{D}\).

With fire spread rate determined, and the geometry of the burned area defined, the area burned in 1 day, \(a_{1{\mathrm{d}}}\) ( \(km^2\, day^{-1}\)), following Li et al. (2012) [61], is calculated as

\[ a_{1{\mathrm{d}}} = \frac{\pi v_\mathrm{d}^2 t^2}{4L_\mathrm{b}}\left(1 + \frac{1}{H_\mathrm{b}}\right)^2 \nonumber\\ = \frac{\pi v_\mathrm{d}^2 (24^2)} {4L_\mathrm{b}}\left(1 + \frac{1}{H_\mathrm{b}}\right)^2\label{aburned} \qquad (Eqn 17)\]

by setting \(t\) equal to \(24\, h\).

The fire extinguishing probability, \(q\), is used to calculate the duration ( \(\tau\), \(days\)) of the fire, which in turn is used to calculated the area burned over the duration of the fire, \({a_{\tau d}}\). \(q\) is represented following Kloster et al. (2010) [49] as

\[ q = 0.5 + \frac{\max\left[0, 0.9 - \exp(-0.025\, p_\mathrm{d})\right]}{2}, \qquad (Eqn 18) \]

which yields a value of \(q\) that varies from 0.5 to 0.95 as population density, \(p_\mathrm{d}\) ( \(\frac{\text{number of people}}{km^{2}}\)), increases from zero to infinity. Higher population density thus implies a higher probability of fire being extinguished. \(q\) represents the probability that a fire will be extinguished on the same day it initiated and the probability that it will continue to the next day is ( \(1-q\)). Assuming individual days are independent, the probability that the fire will still be burning on day \(\tau\) is \((1-q)^\tau\). The probability that a fire will last exactly \(\tau\) days, \(P(\tau)\), is the product of the probability that the fire still exists at day \(\tau\) and the probability it will be extinguished on that day hence \(P(\tau) = q(1-q)^\tau\). This yields an exponential distribution of fire duration whose expected value is

\[ \overline{\tau} = E(\tau) = \sum_{\tau=0}^\infty\, \tau\, q(1-q)^{\tau} = \frac{1-q}{q}.\qquad (Eqn 19) \]

Based on this fire duration and the area burned in 1 day (Eqn. 17), the area burned over the duration of the fire ( \(a_{\tau \mathrm{d}}\)) (but still implemented in 1 day since the model does not track individual fires over their duration, \(km^2\, day^{-1}\)) is calculated as

\[ a_{\tau \mathrm{d}} =E(a_{1{\mathrm{d}}} \tau^2)=\sum_{\tau=0}^ \infty\, a_{1{\mathrm{d}}}\, \tau^2 q(1-q)^{\tau} \\ = a_{1{\mathrm{d}}}\, \frac{(1-q) (2-q)}{q^2}.\nonumber\qquad (Eqn 20) \]

Finally, and reintroducing the PFT index \(\alpha\), the area burned is extrapolated for a PFT \(\alpha\) ( \(A_{\mathrm{b}, \alpha}\), \(km^2\, day^{-1}\)) to the whole grid cell as

\[A_{\mathrm{b}, \alpha}=P_{f, \alpha}\, a_{\tau \mathrm{d}, \alpha} \frac{A_\mathrm{g}f_\alpha}{a_{rep}}, \qquad (Eqn 21) \]

where \(A_\mathrm{g}\) is area of a grid cell ( \(km^2\)), \(f_\alpha\) the fractional coverage of PFT \(\alpha\) and \(a_{rep}\) the representative area of \(500\, km^2\), as mentioned earlier. Area burned over the whole grid cell ( \(A_\mathrm{b}\), \(km^2\, day^{-1}\)) is then calculated as the sum of area burned for individual PFTs,

\[ A_\mathrm{b}=\sum_{\alpha=1}^{N}A_{\mathrm{b}, \alpha}.\qquad (Eqn 22)\]

Fire emits \(CO_2\), other trace gases, and aerosols as biomass is burned while plant mortality and damage due to fire contribute to the litter pool. The emissions of a trace gas/aerosol species \(j\) from PFT \(\alpha\), \(E_{\alpha, j}\) ( \(g species (m^{-2} grid cell area) day^{-1}\)) are obtained from a vector of carbon densities \(\vec{C}_{\alpha} = (C_\mathrm{L}, C_\mathrm{S}, C_\mathrm{R}, C_\mathrm{D})_\alpha\) ( \(kg\, C\, m^{-2}\)) for its leaf, stem, root and litter components, multiplied by a vector of combustion factors \(mho_{\alpha} = (mho_\mathrm{L}, mho_\mathrm{S}, mho_\mathrm{R}, mho_\mathrm{D})_\alpha\), which determines what fraction of leaf, stem, root and litter components gets burned, multiplied by a vector of emissions factors \(\Upsilon_{j} = (\Upsilon_\mathrm{L}, \Upsilon_\mathrm{S}, \Upsilon_\mathrm{R}, \Upsilon_\mathrm{D})_j\) ( \(g species (kg\, C\, dry organic matter)^{-1}\)), and by the area burned \(A_{\mathrm{b}, \alpha}\) for that PFT.

The dot product of \(\vec{C}_{\alpha}\), \(\Upsilon_{j}\) and \(mho_{\alpha} \) thus yields emissions per unit grid cell area of species \(j\) from PFT \(\alpha\),

\[ \label{emiss_combust_factor} {E_{\alpha, j}}= ((\vec{C}_\alpha\cdot mho_{\alpha} )\cdot \Upsilon_{j}) \frac{A_{\mathrm{b}, \alpha}}{A_\mathrm{g}}\frac{1000}{450}, \qquad (Eqn 23) \]

where the constant 1000 converts \(\vec{C}_\alpha\) from \(kg\, C\, m^{-2}\) to \(g\, C\, m^{-2}\) and the constant 450 ( \(g\, C\, (kg dry organic matter)^{-1} \)) converts biomass from carbon units to dry organic matter (Li et al. 2012)[61]. The corresponding loss of carbon ( \(kg\, C\, m^{-2}\, day^{-1}\)) from the three live vegetation components (L, S, R) and the litter pool (D) of PFT \(\alpha\) is given by

\[ \label{emiss_combust_loss} H_{\alpha, i}= C_{\alpha, i} mho_i\left(\frac{A_{\mathrm{b}, \alpha}}{A_\mathrm{g}}\right)\quad i={L, S, R, D}.\qquad (Eqn 24) \]

The PFT-specific combustion factors for leaf ( \(mho_\mathrm{L}\)), stem ( \(mho_{\mathrm{S}}\)), root ( \(mho_{\mathrm{R}}\)) and litter ( \(mho_{\mathrm{D}}\)) components are summarized in classicParams.f90. Emission factors for all species of trace gases and aerosols ( \(CO_2\), \(CO\), \(CH_4\), \(H_2\), \(NHMC\), \(NO_x\), \(N_2O\), total particulate matter, particulate matter less than \(2.5\, \mu m\) in diameter, and black and organic carbon) are read in from the CLASSIC parameters namelist file.

Litter generated by fire is based on similar mortality factors, which reflect a PFT's susceptibility to damage due to fire \(\vec{\Theta}_{\alpha} = (\Theta_\mathrm{L}, \Theta_\mathrm{S}, \Theta_\mathrm{R})_\alpha\) (fraction). The contribution to litter pool of each PFT due to plant mortality associated with fire ( \(kg\, C\, m^{-2}\, day^{-1}\)) is calculated as

\[ \label{eqn_using_mort_factors} {M_{\alpha}}= (\vec{C}_\alpha \cdot \Theta_{\alpha} ) \frac{A_{\mathrm{b}, \alpha}}{A_\mathrm{g}}, \qquad (Eqn 25) \]

which is the sum of contribution from individual live vegetation pools

\[ \label{eqn_using_mort_factors_individual} M_{\alpha, i}= C_{\alpha, i} \Theta_{\alpha, i} \left(\frac{A_{\mathrm{b}, \alpha}}{A_\mathrm{g}} \right) \quad i={L, S, R}.\qquad (Eqn 26) \]

The carbon loss terms associated with combustion of vegetation components and litter ( \(H_{\alpha, i}, i={L, S, R, D}\)) and mortality of vegetation components ( \(M_{\alpha, i}, i={L, S, R}\)) due to fire are used in Rate change equations for carbon pools Eqns 1 and 3, which describe the rate of change of carbon in model's five pools (however, listed there without the PFT subscript \(\alpha\)). The PFT-specific mortality factors for leaf ( \(\Theta_\mathrm{L}\)), stem ( \(\Theta_{\mathrm{S}}\)) and root ( \(\Theta_\mathrm{R}\)) components are listed in classicParams.f90.

When CTEM is run with prescribed PFT fractional cover, the area of PFTs does not change and the fire-related emissions of \(CO_2\), other trace gases and aerosols , and generation of litter act to thin the remaining biomass. When competition between PFTs for space is allowed, fire both thins the remaining biomass and through plant mortality creates bare ground, which is subsequently available for colonization. The creation of bare ground depends on the susceptibility of each PFT to stand replacing fire ( \(\zeta_\mathrm{r}\), fraction) (see also classicParams.f90) and the PFT area burned. The fire-related mortality rate, \(m_{dist}\) ( \(day^{-1}\)), used in mortality.f90 Eqn. 1, is then

\[ \label{m_dist} m_{dist, \alpha} = \zeta_{\mathrm{r}, \alpha} \frac{A_{\mathrm{b}, \alpha}}{f_\alpha A_\mathrm{g}}.\qquad (Eqn 27) \]

After bare ground generation associated with fire, the thinned biomass is spread uniformly over the remaining fraction of a PFT. However, it is ensured that the carbon density of the remaining biomass does not increase to a value above what it was before the fire occurred.

Function/Subroutine Documentation

◆ burntobare()

subroutine, public disturbance_scheme::burntobare ( integer, intent(in)  il1,
integer, intent(in)  il2,
integer, intent(in)  nilg,
integer, dimension(icc), intent(in)  sort,
real, dimension(nilg), intent(in)  pvgbioms,
real, dimension(nilg), intent(in)  pgavltms,
real, dimension(nilg), intent(in)  pgavscms,
real, dimension(nilg,icc), intent(in)  burnvegf,
real, dimension(nilg,icc), intent(in)  pstemmass,
real, dimension(nilg,icc), intent(in)  pgleafmass,
integer, intent(in)  useTracer,
real, dimension(nilg,icc), intent(inout)  fcancmx,
real, dimension(nilg,icc), intent(inout)  stemmass,
real, dimension(nilg,icc), intent(inout)  rootmass,
real, dimension(nilg,icc), intent(inout)  gleafmas,
real, dimension(nilg,icc), intent(inout)  bleafmas,
real, dimension(nilg,iccp2), intent(inout)  litrmass,
real, dimension(nilg,iccp2), intent(inout)  soilcmas,
real, dimension(nilg,icc), intent(inout)  nppveg,
real, dimension(:,:,:), intent(inout)  tracerLitrMass,
real, dimension(:,:,:), intent(inout)  tracerSoilCMass,
real, dimension(:,:), intent(inout)  tracerGLeafMass,
real, dimension(:,:), intent(inout)  tracerBLeafMass,
real, dimension(:,:), intent(inout)  tracerStemMass,
real, dimension(:,:), intent(inout)  tracerRootMass 
)

Update fractional coverages of pfts to take into account the area burnt by fire. Adjust all pools with new densities in their new areas and increase bare fraction. And while we are doing this also run a small check to make sure grid averaged quantities do not get messed up.

Author
Joe Melton
Parameters
[in]nilgno. of grid cells in latitude circle (this is passed in as either ilg or nlat depending on mos/comp)
[in]usetracerSwitch for use of a model tracer. If useTracer is 0 then the tracer code is not used. useTracer = 1 turns on a simple tracer that tracks pools and fluxes. The simple tracer then requires that the tracer values in the init_file and the tracerCO2file are set to meaningful values for the experiment being run. useTracer = 2 means the tracer is 14C and will then call a 14C decay scheme. useTracer = 3 means the tracer is 13C and will then call a 13C fractionation scheme.
[in]sortindex for correspondence between 9 ctem pfts and size 12 of parameter vectors
[in]pvgbiomsinitial veg biomass
[in]pgavltmsinitial litter mass
[in]pgavscmsinitial soil c mass
[in]burnvegfper PFT fraction burned of that PFTs area
[in]pstemmassgrid averaged stemmass prior to disturbance, \(kg c/m^2\)
[in]pgleafmassgrid averaged rootmass prior to disturbance, \(kg c/m^2\)
[in,out]fcancmxinitial fractions of the ctem pfts
[in,out]gleafmasgreen leaf carbon mass for each of the 9 ctem pfts, \(kg c/m^2\)
[in,out]bleafmasbrown leaf carbon mass for each of the 9 ctem pfts, \(kg c/m^2\)
[in,out]stemmassstem carbon mass for each of the 9 ctem pfts, \(kg c/m^2\)
[in,out]rootmassroots carbon mass for each of the 9 ctem pfts, \(kg c/m^2\)
[in,out]nppvegnpp for individual pfts, \(u-mol co_2/m^2.sec\)
[in,out]soilcmassoil carbon mass for each of the 9 ctem pfts + bare, \(kg c/m^2\)
[in,out]litrmasslitter carbon mass for each of the 9 ctem pfts + bare, \(kg c/m^2\)
[in,out]tracergleafmassTracer mass in the green leaf pool for each of the CTEM pfts, \(tracer C units/m^2\)
[in,out]tracerbleafmassTracer mass in the brown leaf pool for each of the CTEM pfts, \(tracer C units/m^2\)
[in,out]tracerstemmassTracer mass in the stem for each of the CTEM pfts, \(tracer C units/m^2\)
[in,out]tracerrootmassTracer mass in the roots for each of the CTEM pfts, \(tracer C units/m^2\)
[in,out]tracerlitrmassTracer mass in the litter pool for each of the CTEM pfts + bareground and LUC products, \(tracer C units/m^2\)
[in,out]tracersoilcmassTracer mass in the soil carbon pool for each of the CTEM pfts + bareground and LUC products, \(kg c/m^2\)

Do some initializations

Account for disturbance creation of bare ground. This occurs with relatively low frequency and is PFT dependent. We must adjust the amount of bare ground created to ensure that we do not increase the density of the remaining vegetation.

Test the pftfraca to ensure it does not cause densification of the exisiting biomass Trees compare the stemmass while grass compares the root mass.

the pstemmass is from before the fire occurred,i.e. no thinning !

adjust the bare frac to accomodate for the changes

the pgleafmass is from before the fire occurred,i.e. no thinning !

adjust the bare frac to accomodate for the changes

Soil and litter carbon are treated such that we actually transfer the carbon to the bare fraction since it would remain in place as a location was devegetated In doing so we do not adjust the litter or soilc density on the remaining vegetated fraction. But we do adjust it on the bare fraction to ensure our carbon balance works out.

only do checks if we actually shifted fractions here.

check if total grid average biomass density is same before and after adjusting fractions

only do checks if we actually shifted fractions here.

◆ disturb()

subroutine, public disturbance_scheme::disturb ( real, dimension(ilg,ignd), intent(in)  thliq,
real, dimension(ilg,ignd), intent(in)  THLW,
real, dimension(ilg,ignd), intent(in)  THFC,
real, dimension(ilg), intent(in)  uwind,
integer, intent(in)  useTracer,
real, dimension(ilg), intent(in)  vwind,
real, dimension(ilg), intent(in)  lightng,
real, dimension(ilg,icc), intent(in)  fcancmx,
integer, dimension(ilg,ignd), intent(in)  isand,
real, dimension(ilg,icc,ignd), intent(in)  rmatctem,
integer, intent(in)  ilg,
integer, intent(in)  il1,
integer, intent(in)  il2,
integer, dimension(icc), intent(in)  sort,
real, dimension(ilg)  grclarea,
real, dimension(ilg,ignd), intent(in)  thice,
real, dimension(ilg), intent(in)  popdin,
real, dimension(ilg), intent(in)  lucemcom,
logical, intent(in)  dofire,
real, dimension(ilg), intent(in)  currlat,
integer, intent(in)  iday,
real, dimension(ilg), intent(in)  fsnow,
real, dimension(ilg,icc), intent(inout)  stemmass,
real, dimension(ilg,icc), intent(inout)  rootmass,
real, dimension(ilg,icc), intent(inout)  gleafmas,
real, dimension(ilg,icc), intent(inout)  bleafmas,
real, dimension(ilg,iccp2), intent(inout)  litrmass,
real, dimension(:,:), intent(inout)  tracerStemMass,
real, dimension(:,:), intent(inout)  tracerRootMass,
real, dimension(:,:), intent(inout)  tracerGLeafMass,
real, dimension(:,:), intent(inout)  tracerBLeafMass,
real, dimension(:,:,:), intent(inout)  tracerLitrMass,
real, dimension(ilg,icc), intent(out)  stemltdt,
real, dimension(ilg,icc), intent(out)  rootltdt,
real, dimension(ilg,icc), intent(out)  glfltrdt,
real, dimension(ilg,icc), intent(out)  blfltrdt,
real, dimension(ilg,icc), intent(out)  glcaemls,
real, dimension(ilg,icc), intent(out)  rtcaemls,
real, dimension(ilg,icc), intent(out)  stcaemls,
real, dimension(ilg,icc), intent(out)  blcaemls,
real, dimension(ilg,icc), intent(out)  ltrcemls,
real, dimension(ilg), intent(out)  burnfrac,
real, dimension(ilg,icc), intent(out)  pstemmass,
real, dimension(ilg,icc), intent(out)  pgleafmass,
real, dimension(ilg,icc), intent(out)  emit_co2,
real, dimension(ilg,icc), intent(out)  emit_ch4,
real, dimension(ilg,icc), intent(out)  emit_co,
real, dimension(ilg,icc), intent(out)  emit_nmhc,
real, dimension(ilg,icc), intent(out)  emit_h2,
real, dimension(ilg,icc), intent(out)  emit_nox,
real, dimension(ilg,icc), intent(out)  emit_n2o,
real, dimension(ilg,icc), intent(out)  emit_pm25,
real, dimension(ilg,icc), intent(out)  emit_tpm,
real, dimension(ilg,icc), intent(out)  emit_tc,
real, dimension(ilg,icc), intent(out)  emit_oc,
real, dimension(ilg,icc), intent(out)  emit_bc,
real, dimension(ilg,icc), intent(out)  burnvegf,
real, dimension(ilg,icc), intent(out)  bterm_veg,
real, dimension(ilg,icc), intent(out)  mterm_veg,
real, dimension(ilg), intent(out)  lterm,
real, dimension(ilg,icc), intent(out)  smfunc_veg 
)

Calculates whether fire occurs, burned area, amount of C emitted and litter generated.

Author
Vivek Arora and Joe Melton
Parameters
[in]il1il1=1
[in]il2il2=ilg
[in]sortindex for correspondence between 9 pfts and size 12 of parameters vectors
[in]dofireboolean, if true allow fire, if false no fire
[in]usetracerSwitch for use of a model tracer. If useTracer is 0 then the tracer code is not used. useTracer = 1 turns on a simple tracer that tracks pools and fluxes. The simple tracer then requires that the tracer values in the init_file and the tracerCO2file are set to meaningful values for the experiment being run. useTracer = 2 means the tracer is 14C and will then call a 14C decay scheme. useTracer = 3 means the tracer is 13C and will then call a 13C fractionation scheme.
[in]thliqliquid soil moisture content
[in]thlwwilting point soil moisture content
[in]thfcfield capacity soil moisture content
[in]uwindwind speed, \(m/s\)
[in]vwindwind speed, \(m/s\)
[in]fcancmxfractional coverages of ctem's 9 pfts
[in]lightngtotal \(lightning, flashes/(km^2 . year)\) it is assumed that cloud to ground lightning is some fixed fraction of total lightning.
[in]rmatctemfraction of roots in each soil layer for each pft
[in]thicefrozen soil moisture content over canopy fraction
[in]popdinpopulation density \((people / km^2)\)
[in]lucemcomland use change (luc) related combustion emission losses, \(u-mol co2/m2.sec\)
[in]fsnowfraction of snow simulated by class
[in,out]stemmassstem mass for each of the 9 ctem pfts, \(kg c/m^2\)
[in,out]rootmassroot mass for each of the 9 ctem pfts, \(kg c/m^2\)
[in,out]gleafmasgreen leaf mass for each of the 9 ctem pfts, \(kg c/m^2\)
[in,out]bleafmasbrown leaf mass
[in,out]litrmasslitter mass for each of the CTEM pfts
[in,out]tracergleafmassTracer mass in the green leaf pool for each of the CTEM pfts, \(tracer C units/m^2\)
[in,out]tracerbleafmassTracer mass in the brown leaf pool for each of the CTEM pfts, \(tracer C units/m^2\)
[in,out]tracerstemmassTracer mass in the stem for each of the CTEM pfts, \(tracer C units/m^2\)
[in,out]tracerrootmassTracer mass in the roots for each of the CTEM pfts, \(tracer C units/m^2\)
[in,out]tracerlitrmassTracer mass in the litter pool for each of the CTEM pfts + bareground and LUC products, \(tracer C units/m^2\)
[out]stemltdtstem litter generated due to disturbance \((kg c/m^2)\)
[out]rootltdtroot litter generated due to disturbance \((kg c/m^2)\)
[out]glfltrdtgreen leaf litter generated due to disturbance \((kg c/m^2)\)
[out]burnvegfper PFT fraction burned of that PFTs area
[out]bterm_vegbiomass fire probability term, Vivek
[out]mterm_vegmoisture fire probability term, Vivek
[out]ltermlightning fire probability term
[out]smfunc_vegsoil moisture dependence on fire spread rate, Vivek
[out]glcaemlsgreen leaf carbon emission losses, \(kg c/m^2\)
[out]rtcaemlsroot carbon emission losses, \(kg c/m^2\)
[out]stcaemlsstem carbon emission losses, \(kg c/m^2\)
[out]ltrcemlslitter carbon emission losses, \(kg c/m^2\)
[out]blfltrdtbrown leaf litter generated due to disturbance \((kg c/m^2)\)
[out]blcaemlsbrown leaf carbon emission losses, \(kg c/m^2\)
[out]burnfractotal areal :: fraction burned, (%)
[out]emit_co2carbon dioxide
[out]emit_cocarbon monoxide
[out]emit_ch4methane
[out]emit_nmhcnon-methane hydrocarbons
[out]emit_h2hydrogen gas
[out]emit_noxnitrogen oxides
[out]emit_n2onitrous oxide
[out]emit_pm25particulate matter less than 2.5 um in diameter
[out]emit_tpmtotal particulate matter
[out]emit_tctotal carbon
[out]emit_ocorganic carbon
[out]emit_bcblack carbon
grclareagcm grid cell area, \(km^2\)

initialize required arrays to zero, or assign value

initialization ends

Find pft areas before

area in \(km^2\)

Find the probability of fire as a product of three functions with dependence on total biomass, soil moisture, and lightning

1 Dependence on total biomass

don't allow it to bring in crops since they are not allowed to burn.

Root biomass is not used to initiate fire. For example if the last fire burned all grass leaves, and some of the roots were left, its unlikely these roots could catch fire. Here we ignore the LUC litrmass on iccp2 and the litter on the bare ground (iccp1). We only consider the litrmass on layer 1 as the rest are buried.

Find average biomass over the vegetated fraction

Sum up the vegetated area

calculate bterm for individual PFTs as well

2 Dependence on soil moisture

This is calculated in a way such that the more dry the root zone of a pft type is, and more fractional area is covered with that pft, the more likely it is that fire will get started. that is the dryness factor is weighted by fraction of roots in soil layers, as well as according to the fractional coverage of different pfts. the assumption here is that if there is less moisture in root zone, then it is more likely the vegetation will be dry and thus the likeliness of fire is more.

First find the dryness factor for each soil layer.

If there is snow on the ground, similarly if not soil, do not allow fire so set betadrgt to 0 for all soil layers otherwise calculate as per normal.

Now find weighted value of this dryness factor averaged over the rooting depth, for each pft

Next find this dryness factor averaged over the vegetated fraction \(avgdryns(i) = avgdryns(i) + drgtstrs(i,j)*fcancmx(i,j)\)

The litter and brown leaves are not affected by the soil water potential therefore they will react only to the moisture conditions (our proxy here is the upper soil moisture). If it is dry they increase the probability of fire corresponding to the proportion of total C they contribute. Only allow if there is no snow.

Use average root zone vegetation dryness to find likelihood of fire due to moisture.

calculate mterm for each PFT

duff fraction for each PFT, Vivek We only consider the first layer litter

\(drgtstrs(i,j)\) is \(\phi_{root}\) in Melton and Arora GMDD (2015) paper

no fire likelihood due to moisture if no vegetation

3 dependence on lightning

Dependence on lightning is modelled in a simple way which implies that a large no. of lightning flashes are more likely to cause fire than few lightning flashes.

New approximation of Price and Rind equation. It was developed from a more complete dataset than Prentice and Mackerras. Lightning comes in in units of \(flashes/km^2/yr\) so divide by 12 to make per month.

Determine the probability of fire due to human causes this is based upon the population density from the .popd read-in file

–——————— Number of fire calculations -———————\

This is not used in CTEM in general but we keep the code here for testing purposes

calculate natural and anthorpogenic ignitions/km2.day the constant 0.25 assumes not all c2g lightning hits causes ignitions, only 0.25 of them do the constant (1/30.4) converts c2g lightning from flashes/km2.month to flashes/km2.day MAKE SURE LIGHTNING IS IN UNITS OF FLASHES/KM2.MONTH

Eqs. (4) and (5) of Li et al. 2012 doi:10.5194/bg-9-2761-2012 + also see corrigendum

calculate fire suppression also as a function of population density. Li et al. (2012) formulation is quite similar to what we already use based on Kloster et al. (2010, I think) but we use Kloster's formulation together with our fire extingishing probability. Here, the fire suppression formulation is just by itself

–——————— Number of fire calculations -———————//

calculate fire probability for each PFT. Recall that lightning term is still grid averaged

–——————— Number of fire calculations -———————\

This is not used in CTEM in general but we keep the code here for testing purposes

calculate total number of ignitions based natural and anthorpogenic ignitions for the whole grid cell

finally calculate number of fire,noting that not all ignitions turn into fire because moisture and biomass may now allow that to happen,and some of those will be suppressed due to fire fighting efforts

–——————— Number of fire calculations -———————//

Calculate area burned for each PFT, make sure it's not greater than the area available, then find area and fraction burned for the whole gridcell

soil moisture dependence on fire spread rate

wind speed, which is gridcell specific

Length to breadth ratio from Li et al. (2012)

head to back ratio from Li et al. (2012).

dependence of spread rate on wind

fire spread rate per PFT

area burned in 1 day for that PFT

fire extinguishing probability as a function of grid-cell averaged population density

account for low suppression in Savanna regions, see above for increase in ignition due to cultural practices

area multipler to calculate area burned over the duration of the fire

per PFT area burned, \(km^2\)

Calculate gridcell area burned and fraction

Finally estimate amount of litter generated from each pft, and each vegetation component (leaves, stem, and root) based on their resistance to combustion. Update the veg pools due to combustion.

Set aside these pre-disturbance stem and root masses for use in burntobare subroutine.

Update the pools:

Output the burned area per PFT (the units here are burned fraction of each PFTs area. So if a PFT has 50% gridcell cover and 50% of that burns it will have a burnvegf of 0.5 (which then translates into a gridcell fraction of 0.25). This units is for consistency outside of this subroutine.

Even if no fire we still enter here and perform the calculations for the emissions since we might have some from luc.

We also estimate \(CO_2\) emissions from each of these components. Note that the litter which is generated due to disturbance is uniformly distributed over the entire area of a given pft, and this essentially thins the vegetation biomass. If PFTCompetition is not on, this does not change the vegetation fractions, if competition is on a fraction will become bare. That is handled in burntobare subroutine called from competition subroutine.

Calculate the emissions of trace gases and aerosols based upon how much plant matter was burnt

Sum all pools that will be converted to emissions/aerosols \((g c/m^2 /day)\)

Add in the emissions due to luc fires (deforestation) the luc emissions are converted from \(umol co_2 m-2 s-1 to g c m-2\) (day-1) before adding to tot_emit. NOTE: This tot_emit is not what is used to calculate NBP so the LUC emissions are accounted for in ctem.f90, not here. This here is for model outputting.

Convert burnt plant matter from carbon to dry organic matter using a conversion factor, assume all parts of the plant has the same ratio of carbon to dry organic matter. units: \(kg dom / m^2 / day\)

Convert the dom to emissions/aerosols using emissions factors units: \(g species / m^2 /d\)

Also convert units to \(kg species / m^2 / s^{-1}\)