-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathcalc_cfwi.ncl
145 lines (115 loc) · 4.03 KB
/
calc_cfwi.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
/;
2018, Copyright University Corporation for Atmospheric Research
;/
/;
tmax: daily maximum temperature
prec: 24 hour precipitation
hurs: daily relative humidity
spd: daily wind speed
time: time series of all variables
lat: latitude coordinates for all variables
lon: longitude coordinates for all variables
iounits: integer array corresponding to units of tmax, prec, hurs, spd
iounit(0)=-1 input tmax (metadata)
iounit(0)=0 input tmax (degC)
iounit(0)=1 input tmax (degK)
iounit(0)=2 input tmax (degF)
iounit(1)=-1 input prec (metadata)
iounit(1)=0 input prec (mm/day)
iounit(1)=1 input prec (inches/day)
iounit(1)=2 input prec (kg m-2 s-1)
iounit(2)=-1 input rhmax (metadata)
iounit(2)=0 input rhmax (%)
iounit(2)=1 input rhmax (1)
iounit(3)=-1 input spd (metadata)
iounit(3)=0 input spd (m s-1)
iounit(3)=1 input spd (mph)
;/
load "calc_ffmc.ncl"
load "calc_dmc.ncl"
load "calc_dc.ncl"
load "calc_isi.ncl"
load "calc_bui.ncl"
load "convert_temp.ncl"
load "convert_prec.ncl"
load "convert_humid.ncl"
load "convert_wind.ncl"
load "unit_handling.ncl"
load "check_time_length.ncl"
function calc_cfwi(tmax, prec, hurs, spd, time, lat, lon, iounits, opt)
local fd, B, cfwi, ffmc, dmc, dc, Lf, latlen, lonlen, ndays, cal, cal1, day_year, Le, Le1, Le2, tmaxdmc, tmaxdc, m, ffmc, dmc, dc, isi, bui, fd, B, varatts, size
begin
size = dimsizes(iounits)
if(size .eq. 1) then
if(iounits .eq. -1) then
convert_humid("%", hurs)
convert_prec("mm/day", prec)
convert_temp("degC", tmax)
convert_wind("m/s", spd)
else
print("The iounits value "+ iounits +" is not recognized. To indicate that all input units should be taken from metadata, input the value -1. Otherwise, indicate units according to documentation in a vector of (/tmax, prec, hurs, spd/).")
end if
end if
if(size .gt. 1) then
units = unit_handling((/"tmax", "prec", "hurs", "spd"/), iounits)
if(units(0) .ne. "metadata") then
tmax@units = units(0)
end if
if(units(1) .ne. "metadata") then
prec@units=units(1)
end if
if(units(2) .ne. "metadata") then
hurs@units = units(2)
end if
if(units(3) .ne. "metadata") then
spd@units = units(3)
end if
convert_humid("%", hurs)
convert_prec("mm/day", prec)
convert_temp("degC", tmax)
convert_wind("m/s", spd)
end if
ffmc = prec(0, :, :)
dmc = prec(0, :, :)
dc = prec(0, :, :)
m = prec(0, :, :)
ffmc = 85.0
dmc = 6.0
dc = 15.0
m = 0.0
Lf = (/-1.6, -1.6, -1.6, .9, 3.8, 5.8, 5.8, 6.4, 5.0, 2.4, .4, -1.6, -1.6/) ;"Day length factors" for deeper in the soil... according to month, Jan-Dec
cfwi = prec
cfwi = 0.0
latlen = dimsizes(lat)
lonlen = dimsizes(lon)
ndays = dimsizes(time)
check_time_length(ndays, 365, opt)
do i = 0, ndays-1
cal = cd_calendar(time(i), 0)
cal1 = tointeger(cal(0,:))
cal1@calendar = time@calendar
day_year = day_of_year(cal1(0), cal1(1), cal1(2))
Le = daylight_fao56(day_year, lat)
Le1 = Le(0, :)
Le2 = conform_dims((/latlen, lonlen/), Le1, 0)
m = (/ calc_ffmc(ffmc, hurs(i, :, :), tmax(i, :, :), spd(i, :, :), prec(i, :, :) ) /)
ffmc = 59.5 * (250. - m) / (147.2 + m)
dmc = (/ calc_dmc(dmc, tmax(i, :, :), prec(i, :, :), hurs(i, :, :), Le2) /)
dc = (/ calc_dc(prec(i, :, :), tmax(i, :, :), dc, Lf(cal1(1)-1)) /)
isi = (/ calc_isi(spd(i, :, :), m) /)
bui = (/ calc_bui(dmc, dc) /)
fd = where(bui .gt. 80., 1000. / (25. + 108.64 * exp(-0.023 * bui)), 0.626 * bui^0.809 + 2)
B = 0.1 * isi * fd
cfwi(i, :, :) = where(B .gt. 1., exp(2.72 * (.434 * log(B))^.647), B)
cfwi(i, :, :) = cfwi(i, :, :) > 0.0
end do
delete_VarAtts(cfwi, -1)
cfwi@long_name = "Canadian Fire Weather Index"
cfwi@standard_name = "canadian_fire_weather_index"
varatts = (/"units", "missing_value", "_FillValue"/)
cfwi@$varatts(0)$ = "1" ; cfwi is unitless
do i = 1, dimsizes(varatts)-1
cfwi@$varatts(i)$ = prec@$varatts(i)$
end do
return(cfwi)
end