[67] | 1 | function read_grib, varcode, date1, date2, file = file |
---|
[114] | 2 | ; |
---|
| 3 | compile_opt idl2, strictarrsubs |
---|
| 4 | ; |
---|
[67] | 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' |
---|
[69] | 52 | time = date2jul(dates) |
---|
[67] | 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 |
---|