-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathcalc_just_bi.ncl
178 lines (136 loc) · 9.89 KB
/
calc_just_bi.ncl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
/;
2020, Copyright University Corporation for Atmospheric Research
;/
/;
mc1, mc10, mc100, are the percent moisture content of 1-, 10-, 100-hour timelag
climat is NFDRS climate class
erc is energy release component
bi is burning index
;/
load "retrieve_constants.ncl"
function calc_just_bi(mc1, mc10, mc100, erc, mcherb, mcwood, ws, fuelmod, climat)
local w1,w10,w100,w1000,wherb,wwood,depth,sg1,sg10,sg100,sg1000,sgherb,sgwood,mxd,hd,wndfc,hl,rhol,rhod,std,stl,sd,sl,fctcur,wherbc,w1p,sa1,sa10,sa100,sawood,saherb,sadead,salive,f1,f10,f100,wherbp,fdead,flive,fherb,fwood,sgbrd,sgbrl,wtotd,wtotl,wtot,sa1,sa10,sa100,sawood,saherb,w1n,w10n,w100n,wherbn,wwoodn,hn1,hn10,hn100,hnherb,hnwood,wrat,mclfe,sgbrt,rhobed,rhobar,wtmcd,wtmcl,mxl,gmamx,betbar,betop,ad,dedrt,fherb,livrt,gmaop,wdeadn,etasd,etamd,wliven,etasl,etaml,slope_angle,slpfct,c,e,b,ufact,ir,zeta,phislp,phiwnd,htsink,sc,bi
begin
constants = retrieve_constants(fuelmod)
w1 = constants(0) ; fuel model specified 1-hour class fuel loading
w10 = constants(1) ; fuel model specified 10-hour class fuel loading
w100 = constants(2) ; fuel model specified 100-hour class fuel loading
w1000 = constants(3) ; fuel model specified 1000-hour class fuel loading
wwood = constants(4) ; fuel model specified woody shrub fuel loading
wherb = constants(5) ; fuel model specified herbaceous fuel loading
sg1 = constants(6) ; specified surface area-to-volume ratio for 1-hour class of fuel model
sg10 = constants(7) ; specified surface area-to-volume ratio for 10-hour class of fuel model
sg100 = constants(8) ; specified surface area-to-volume ratio for 100-hour class of fuel model
sgwood = constants(10) ; specified surface area-to-volume ratio for woody class of fuel model
sgherb = constants(11) ; specified surface area-to-volume ratio for herbaceous class of fuel model
hd = constants(12) ; specified dead-fuel heat of combustion of fuel model
mxd = constants(14) ; specified dead-fuel moisture of extinction for fuel model
depth = constants(15) ; effective fuel-bed depth measured, in feet
wndfc = constants(16) ; fuel model specified wind reduction factor. wndfc is calculated under certain conditions.
c1 = 0.046 ; conversion factor for T/Ac to lbs/ft^2
w1 = w1 * c1
w10 = w10 * c1
w100 = w100 * c1
w1000 = w1000 * c1
wwood = wwood * c1
wherb = wherb * c1 ; conversions from T/Ac to lbs/ft^2
hl = hd ; specified live-fuel heat of combustion of fuel model
rhol = 32.0 ; particle density of live fuel
rhod = 32.0 ; particle density of dead fuel
std = 0.0555 ; proportion of inert mineral content of dead fuels
stl = 0.0555 ; proportion of inert mineral content of live fuels
sd = 0.01 ; proportion of silica-free mineral content of dead fuels
sl = 0.01 ; proportion of silica-free mineral content of live fuels
fctcur = 1.33 - 0.0111 * mcherb ; fraction of fuel model herbaceous fuel loading transferred to 1-hour fuel class
fctcur = fctcur < 1.0
fctcur = fctcur > 0.0
wherbc = fctcur * wherb ; amount of herbaceous fuel loading transferred to 1-hour class
w1p = w1 + wherbc ; 1-hour fuel loading and transferred herbaceous loading
sa1 = (w1p / rhod) * sg1 ; surface area of 1-hour class fuels
sa10 = (w10 / rhod) * sg10 ; surface area of 10-hour class fuels
sa100 = (w100 / rhod) * sg100 ; surface area of 100-hour class fuels
sawood = (wwood / rhol) * sgwood ; surface area of woody fuel class
saherb = (wherb / rhol) * sgherb ; surface area of herbaceous fuel class
sadead = sa1 + sa10 + sa100 ; surface area of dead-fuel classes
salive = saherb + sawood ; surface area of live-fuel classes
f1 = sa1 / sadead ; proportion of dead-fuel surface area in 1-hour class
f10 = sa10 / sadead ; proportion of dead-fuel surface area in 10-hour class
f100 = sa100 / sadead ; proportion of dead-fuel surface area in 100-hour class
; f1, f10, and f100 are used as weighting factors in the ros/sc calculation
wherbp = wherb - wherbc ; amount of herbaceous floading left after transfer to 1-hour fuel loading
fdead = sadead / (sadead + salive) ; proportion of total surface area in dead-fuel classes
flive = salive / (sadead + salive) ; proportion of total surface area in live-fuel classes
fherb = saherb / salive ; proportion of live surface area in herbaceous class
fwood = sawood / salive ; proportion of live surface area in woody class
sgbrd = (f1 * sg1) + (f10 * sg10) + (f100 * sg100) ; characteristic surface area-to-volume ratio of dead fuel, surface area weighted
sgbrl = (fherb * sgherb) + (fwood * sgwood) ; characteristic surface area-to-volume ratio of live fuel, surface area weighted
wtotd = w1p + w10 + w100 + w1000 ; total dead-fuel loading
wtotl = wherbp + wwood ; total live-fuel loading
wtot = wtotd + wtotl ; total fuel loading
w1n = w1p * (1.0 - std) ; net fuel loading of 1-hour class
w10n = w10 * (1.0 - std) ; net fuel loading of 10-hour class
w100n = w100 * (1.0 - std) ; net fuel loading of 100-hour class
wherbn = wherbp * (1.0 - stl) ; net fuel loading of herbaceous class
wwoodn = wwood * (1.0 - stl) ; net fuel loading of woody class
wtotln = wtotl * (1 - stl) ; total live-fuel loading
hn1 = w1n * exp( -138.0 / sg1) ; heating number of 1-hour class
hn10 = w10n * exp( -138.0 / sg10) ; heating number of 10-hour class
hn100 = w100n * exp( -138.0 / sg100) ; heating number of 100-hour class
if((-500. / sgherb) .lt. -180.218) then
hnherb = conform_dims(dimsizes(wherbn), 0., -1) ; heating number of herbaceous class
else
hnherb = wherbn * exp(-500. / sgherb)
end if
hnwood = where((-500. / sgwood) .lt. -180.218, 0., wwoodn * exp(-500. / sgwood)) ; heating number of woody class
wrat = where((hnherb + hnwood) .eq. 0., 0., (hn1 + hn10 + hn100) / (hnherb + hnwood)) ; ratio of dead-to-live heating numbers for calculation of live moisture of extinction (mxl). note that contributions from the 1000-hour class are omitted because their contributions are negligible
mclfe = ((mc1 * hn1) + (mc10 * hn10) + (mc100 * hn100)) / (hn1 + hn10 + hn100) ; dead-fuel moisture content, weighted by heating number, for calculation of live moisture of extinction (mxl)
sgbrt = (fdead * sgbrd) + (flive * sgbrl) ; characteristic surface area-to-volume ratio of fuel bed, surface area weighted
rhobed = (wtot - w1000) / depth ; bulk density of fuel bed
rhobar = ((wtotl * rhol) + (wtotd * rhod)) / wtot ; particle density of weighted fuel
wtmcd = (f1 * mc1) + (f10 * mc10) + (f100 * mc100) ; surface area weighted dead-fuel moisture content
wtmcl = (fherb * mcherb) + (fwood * mcwood) ; surface area weighted live-fuel moisture content
mxl = where(wtotln .gt. 0, (2.9 * wrat * (1.0 - mclfe / mxd) - 0.226) * 100.0 > mxd, 0) ; calculated live-fuel moisture of extinction
gmamx = (sgbrt^1.5) / (495.0 + 0.0594 * sgbrt^1.5) ; weighted maximum reaction velocity of surface area
betbar = rhobed/rhobar ; packing ratio
betop = 3.348 * sgbrt^(-0.8189) ; optimum packing ratio, surface area weighted
ad = 133.0 * sgbrt^(-0.7913) ; exponent in gmaop equation
dedrt = wtmcd / mxd ; ratio in calculation of etamd
fherb = saherb / salive ; proportion of live surface area in herbaceous class
livrt = wtmcl / mxl ; ratio in calculation of etaml
gmaop = gmamx * (betbar / betop)^ad * exp(ad * (1.0 + betbar / betop)) ; weighted optimum reaction velocity of surface area
wdeadn = (f1 * w1n) + (f10 * w10n) + (f100 * w100n) ; net loading of dead fuels, surface-area weighted
etasd = 0.174 * sd^(-0.19) ; mineral damping coefficient of dead fuel
etamd = 1.0 - (2.59 * dedrt) + (5.11 * (dedrt^2.0)) - (3.52 * (dedrt^3.0)) ; surface area weighted dead-fuel moisture damping coefficient
etamd = etamd > 0 ;NEW LINE 3/9/2020
etamd = etamd < 1 ;NEW LINE 3/9/2020
wliven = (fwood * wwoodn) + (fherb * wherbn) ; surface area weighted net loading of live fuels
etasl = 0.174 * sl^(-0.19) ; live fuel mineral damping coefficient
etaml = 1.0 - (2.59 * livrt) + (5.11 * (livrt^2.0)) - (3.52 * (livrt^3.0)) ; surface area weighted live-fuel moisture damping coefficient
etaml = etaml > 0 ;NEW LINE 3/9/2020
etaml = etaml < 1 ;NEW LINE 3/9/2020 these four new lines are what make BI not break :)
;slope angle depends on nfdrs slope class:
slope_angles = (/12.67, 17.63, 24.23, 32.46, 41.99/) ; corresponding to NFDRS slope class (/1,2,3,4,5/) ****Not the same as NFDRS climate class****
slope_class = 1 ;default taken from MATLAB code. If better information can be found, replace this line.
slope_angle = slope_angles(slope_class - 1)
slpfct = 5.275 * (tan(slope_angle * get_pi("float") / 180.0)) ^ 2 ; slope effect multiplier coefficient
;wind effect multiplier coefficients/exponents
c = 7.47 * exp( -0.133 * sgbrt^0.55)
e = 0.715 * exp(-3.59 * 10^(-4) * sgbrt)
b = 0.02526 * sgbrt^0.54
ufact = c * (betbar / betop)^(-e)
ir = gmaop * ((wdeadn * hd * etasd * etamd) + (wliven * hl * etasl * etaml)) ; surface area weighted reaction intensity
zeta = exp((0.792 + 0.681 * sgbrt^0.5) * (betbar + 0.1)) / (192.0 + 0.2595 * sgbrt) ; no wind propogating flux ratio
phislp = slpfct * betbar^(-0.3) ; multiplier for slope effect
; if((ws * 88.0 *wndfc) .gt. 0.9 * ir) then
; phiwnd = ufact * ((0.9 * ir)^b) ; wind effect multiplier
; else
; phiwnd = ufact * ((ws * 88.0 * wndfc)^b)
; end if
phiwnd = where((ws * 88.0 * wndfc) .gt. (0.9 * ir), ufact * ((0.9 * ir)^b), ufact * ((ws * 88.0 * wndfc)^b))
; phiwnd = where((ws * 88.0 * wndfc) .gt. (0.9 * ir), ufact * ((0.9 * ir)), 1)
htsink = rhobed * (fdead * ( f1 * exp(-138.0 / sg1) * (250.0 * 11.16 * mc1) + f10 * exp(-138.0 / sg10) * (250.0 * 11.16 * mc10) + f100 * exp(-138.0 / sg100) * (250.0 * 11.16 * mc100))) + (flive * (fherb * exp(-138.0 / sgherb) * ( 250.0 + 11.16 * mcherb) + fwood * exp(-138.0 / sgwood) * ( 250.0 + 11.16 * mcwood))) ; heat sink
;the following calculation is technically a variable called ros, or Rate of Spread. However, sc is just a rounding of ros.
sc = ir * zeta * (1 + phislp + phiwnd) / htsink ; spread component
bi = 3.01 * ((sc * erc)^0.46) ; burning index
return(bi)
end