-
Notifications
You must be signed in to change notification settings - Fork 0
/
get_atcf_output.ncl
executable file
·123 lines (90 loc) · 3.22 KB
/
get_atcf_output.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
;
; Alicia Bentley
; Last Updated: 25 September 2018
;
; Basic Description: Obtains ATCF output needed for a script
;
;
;###############################################################################
undef("get_atcf")
function get_atcf(name,model,time_str,opt)
begin
atcf_url = "./"
; INPUTS ARE:
; name = "aal062018"
; model = "AP01"
; time_str = "2018090100"
; opt = True
;;;; parse data file to obtain GFS track
trk_read = asciiread(atcf_url+name+".dat",-1,"string")
delim = ", "
nfilds = str_fields_count(trk_read(0), delim)
;;;; okay lets get the one tcg pathway letter
time = str_get_field(trk_read, 3, delim)
slat = str_get_field(trk_read, 7, delim)
slon = str_get_field(trk_read, 8, delim)
model_type = str_get_field(trk_read, 5, delim)
fhr = str_get_field(trk_read,6, delim)
wind = str_get_field(trk_read,9,delim)
mslp = str_get_field(trk_read,10,delim)
radii = str_get_field(trk_read,12,delim)
;;;; now that we have done that, lets find a way to discriminate
;;;; pick a model, pick a time and pick a radii
cv_yr = toint(str_get_cols(time_str,0,3))
cv_mo = toint(str_get_cols(time_str,4,5))
cv_day = toint(str_get_cols(time_str,6,7))
cv_hr = toint(str_get_cols(time_str,8,9))
base_time = cd_inv_calendar(cv_yr,cv_mo,cv_day,cv_hr,0,0,"hours since 1800-01-01 00:00:00",0)
;print(cd_string(base_time,""))
atcf_ind = ind(model_type .eq. model .and. time .eq. time_str .and. radii .eq. "34")
atcf_lat = tofloat(slat(atcf_ind))
atcf_lon = tofloat(slon(atcf_ind))
atcf_wind = tofloat(wind(atcf_ind))
atcf_mslp = tofloat(mslp(atcf_ind))
atcf_time = toint(time(atcf_ind))
atcf_fhr = toint(fhr(atcf_ind))
atcf_radii = tofloat(radii(atcf_ind))
;;;; fix lat and lon to coordinates we can use!!!
;print(atcf_lon)
hemi = str_get_cols(slon(atcf_ind),4,4)
;print(hemi)
atcf_lat = atcf_lat / 10
do k = 0,dimsizes(atcf_lon)-1
if hemi(k) .eq. "W" then
atcf_lon(k) = 360 + ((atcf_lon(k) / 10) * -1)
else if hemi(k) .eq. "E" then
atcf_lon(k) = (atcf_lon(k) / 10)
else
atcf_lon(k) = 360 + ((atcf_lon(k) / 10) * -1)
end if
end if
end do
;print(atcf_lon)
;;;; okay last part is put lat lon and time all in one file
atcf_actual_ftime = atcf_fhr + base_time
atcf_actual_ftime@units = "hours since 1800-01-01 00:00:00"
atcf_lat_lon = new((/dimsizes(atcf_lat),2/),"float")
atcf_lat_lon(:,0) = atcf_lat
atcf_lat_lon(:,1) = atcf_lon
atcf_lat_lon!0 = "time"
;atcf_lat_lon&time = cd_convert(atcf_actual_ftime,time_tot@units)
atcf_lat_lon&time = atcf_actual_ftime
;print(atcf_lat_lon(:,0))
;print(atcf_lat_lon(:,1))
;print(atcf_lat_lon&time)
atcf_hr_times = cd_convert(atcf_lat_lon&time,"hours since 1800-01-01 00:00:00")
;print(atcf_hr_times)
;-------------------------------------------------------
atcf_fnl_output := new((/5,dimsizes(atcf_lat_lon&time)/),"float")
atcf_fnl_output(0,:) = (/tofloat(atcf_lat_lon&time)/)
atcf_fnl_output(1,:) = (/atcf_mslp/)
atcf_fnl_output(2,:) = (/atcf_lat_lon(:,0)/)
atcf_fnl_output(3,:) = (/atcf_lat_lon(:,1)/)
atcf_fnl_output(4,:) = (/atcf_wind/)
;printVarSummary(atcf_fnl_output)
atcf_fnl_output!1 = "time"
printVarSummary(atcf_fnl_output)
atcf_fnl_output&time = atcf_hr_times
atcf_fnl_output&time@units = opt@timeUnits
return(atcf_fnl_output)
end