1 | function read_grib, varcode, date1, date2, file = file |
---|
2 | ; |
---|
3 | compile_opt idl2, strictarrsubs |
---|
4 | ; |
---|
5 | @common |
---|
6 | ; http://www.wmo.ch/web/www/WDM/Guides/Guide-binary-2.html |
---|
7 | ; |
---|
8 | ; gribfile = |
---|
9 | ; '/d1fes2-raid6/SINTEX-common/ES10.d.00/atm/5d/ES10.d.00_5d_00911201_00911230.grib' |
---|
10 | IF keyword_set(file) THEN gribfile = isafile(file = file, iodir = iodir) ELSE gribfile = '/d1fes2-raid6/SINTEX-common/ES10/atm/5d/ZOOM_IND/ES10_5d_00210101_00301230.grib' |
---|
11 | ; |
---|
12 | openr, num, gribfile, /GET_LUN, ERROR = err, /SWAP_IF_LITTLE_ENDIAN |
---|
13 | if err ne 0 then begin |
---|
14 | print, !err_string |
---|
15 | return, -1 |
---|
16 | ENDIF |
---|
17 | ; |
---|
18 | recstart = scan_grib_recstart(num) |
---|
19 | ; |
---|
20 | ; messize = scan_grib_messize(num, recstart) |
---|
21 | ; addoff = lonarr(n_elements(recstart)) |
---|
22 | ; FOR i = 1L, n_elements(recstart)-1 DO $ |
---|
23 | ; addoff[i] = recstart[i]-recstart[i-1]-messize[i-1] |
---|
24 | ; |
---|
25 | ; |
---|
26 | ; nbits = scan_grib_nbits(num, recstart) |
---|
27 | ; print,nbits[uniq(nbits,sort(nbits))] |
---|
28 | ; |
---|
29 | codes = scan_grib_code(num, recstart) |
---|
30 | nbcodes = uniq(codes, sort(codes)) |
---|
31 | ; |
---|
32 | dates = scan_grib_date(num, recstart) |
---|
33 | nbdates = uniq(dates, sort(dates)) |
---|
34 | ; |
---|
35 | goodvar = where(codes EQ varcode) |
---|
36 | IF goodvar[0] EQ -1 THEN BEGIN |
---|
37 | print, 'no var code '+strtrim(varcode, 2)+' in the file' |
---|
38 | return, -1 |
---|
39 | ENDIF |
---|
40 | ; |
---|
41 | recstart = recstart[goodvar] |
---|
42 | dates = dates[goodvar] |
---|
43 | ; |
---|
44 | gooddate = where(dates GE date1 AND dates LE date2) |
---|
45 | IF gooddate[0] EQ -1 THEN BEGIN |
---|
46 | print, 'no dates between '+strtrim(date1, 2)+' and '+strtrim(date2, 2)+' in the file' |
---|
47 | return, -1 |
---|
48 | ENDIF |
---|
49 | recstart = recstart[gooddate] |
---|
50 | dates = dates[gooddate] |
---|
51 | key_caltype = '360d' |
---|
52 | time = date2jul(dates) |
---|
53 | jpt = n_elements(time) |
---|
54 | IF jpt EQ 1 THEN vardate = strtrim(dates[0], 2) ELSE vardate = strtrim(dates[0], 2)+' - '+strtrim(dates[jpt-1], 2) |
---|
55 | ; varname |
---|
56 | vargrid = 'T' |
---|
57 | ; varexp |
---|
58 | ; varunit |
---|
59 | ; |
---|
60 | grib_pds = read_grib_pds(num, recstart[0]) |
---|
61 | ; grid parameters |
---|
62 | IF grib_pds.gdsnotomitted THEN BEGIN |
---|
63 | grib_gds = read_grib_gds(num, recstart[0]) |
---|
64 | ; min, max of the latitude with a precision of 10^-2 |
---|
65 | lat1 = fix(100*grib_gds.la1)/100. |
---|
66 | lat2 = fix(100*grib_gds.la2)/100. |
---|
67 | ; |
---|
68 | CASE grib_gds.gridtype OF |
---|
69 | ; Latitude/Longitude Grid |
---|
70 | 0:BEGIN |
---|
71 | computegrid, grib_gds.lo1, grib_gds.la1, grib_gds.di, -grib_gds.dj $ |
---|
72 | , grib_gds.ni, grib_gds.nj |
---|
73 | END |
---|
74 | ; Gaussian Latitude/Longitude Grid |
---|
75 | 4:BEGIN |
---|
76 | ; find the latitude axis |
---|
77 | CASE 1 OF |
---|
78 | ; n48 |
---|
79 | grib_gds.n EQ 48 AND lat1 EQ 88.57 AND lat2 EQ -88.57:$ |
---|
80 | gphit = n48gaussian() |
---|
81 | ; n80 |
---|
82 | grib_gds.n EQ 80 AND lat1 EQ 89.14 AND lat2 EQ -89.14:$ |
---|
83 | gphit = n80gaussian() |
---|
84 | ; n128 |
---|
85 | grib_gds.n EQ 128 AND lat1 EQ 89.46 AND lat2 EQ -89.46:$ |
---|
86 | gphit = n128gaussian() |
---|
87 | ; n160 |
---|
88 | grib_gds.n EQ 160 AND lat1 EQ 89.57 AND lat2 EQ -89.57:$ |
---|
89 | gphit = n160gaussian() |
---|
90 | ; n256 |
---|
91 | grib_gds.n EQ 256 AND lat1 EQ 89.73 AND lat2 EQ -89.73:$ |
---|
92 | gphit = n256gaussian() |
---|
93 | ; part of one of the gaussian grids defined above |
---|
94 | ELSE:BEGIN |
---|
95 | cnt = 0 |
---|
96 | REPEAT BEGIN |
---|
97 | CASE cnt OF |
---|
98 | 0:gphit = n48gaussian() |
---|
99 | 1:gphit = n80gaussian() |
---|
100 | 2:gphit = n128gaussian() |
---|
101 | 3:gphit = n160gaussian() |
---|
102 | 4:gphit = n256gaussian() |
---|
103 | 5:BEGIN |
---|
104 | gphit = n80gaussian() |
---|
105 | lat1 = 29.71 |
---|
106 | lat2 = -19.62 |
---|
107 | END |
---|
108 | ELSE:stop |
---|
109 | ENDCASE |
---|
110 | nfix = fix(gphit*100)/100. |
---|
111 | nlat1 = (where(nfix EQ lat1))[0] |
---|
112 | nlat2 = (where(nfix EQ lat2))[0] |
---|
113 | IF nlat1 NE -1 AND nlat2 NE -1 $ |
---|
114 | AND nlat2-nlat1+1 EQ grib_gds.nj $ |
---|
115 | THEN gphit = gphit[nlat1:nlat2] ELSE gphit = -1 |
---|
116 | cnt = cnt+1 |
---|
117 | ENDREP UNTIL gphit[0] NE -1 |
---|
118 | END |
---|
119 | ENDCASE |
---|
120 | computegrid, grib_gds.lo1, -1, grib_gds.di, -1, grib_gds.ni, -1, YAXIS = gphit |
---|
121 | END |
---|
122 | ; Mercator Projection Grid |
---|
123 | gridtype EQ 1: |
---|
124 | ; Gnomonic Projection Grid |
---|
125 | gridtype EQ 2: |
---|
126 | ; Lambert Conformal, secant or tangent, conical or bipolar (normal or |
---|
127 | ; oblique) Projection Grid |
---|
128 | gridtype EQ 3: |
---|
129 | ; Polar Stereographic Projection Grid |
---|
130 | gridtype EQ 5: |
---|
131 | ; Oblique Lambert conformal, secant or tangent, conical or bipolar, |
---|
132 | ; projection |
---|
133 | gridtype EQ 13: |
---|
134 | ; Spherical Harmonic Coefficients |
---|
135 | gridtype EQ 50: |
---|
136 | ; Space view perspective or orthographic grid |
---|
137 | gridtype EQ 90: |
---|
138 | ; reserved - see Manual on Codes |
---|
139 | ELSE: |
---|
140 | ENDCASE |
---|
141 | ENDIF ELSE stop |
---|
142 | ; |
---|
143 | res = fltarr(grib_gds.ni, grib_gds.nj, n_elements(recstart)) |
---|
144 | FOR i = 0, n_elements(recstart)-1 DO BEGIN |
---|
145 | res[*, *, i] = read_grib_bds(num, recstart[i], grib_gds.ni, grib_gds.nj) |
---|
146 | ENDFOR |
---|
147 | ; |
---|
148 | free_lun, num |
---|
149 | ; |
---|
150 | IF keyword_set(key_yreverse) THEN res = reverse(res, 2) |
---|
151 | |
---|
152 | RETURN, res |
---|
153 | END |
---|