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