source: trunk/SRC/ToBeReviewed/LECTURE/GRIB/read_grib.pro @ 495

;+
; @file_comments
;
; @categories
;
; @param VARCODE
;
; @param DATE1 {in}{optional}
; Date of the beginning (yyyymmdd if TIMESTEP is not activate)
;
; @param DATE2 {in}{optional}
; Last date. Optional, if not specified date2=date1
;
; @keyword FILE{type=array or string}
; A single filename or an array of filenames to be retrieved.
;
; @returns
;
; @restrictions
;
; @examples
;
; @history
;
; @version
; $Id$
;
;-
FUNCTION read_grib, varcode, date1, date2, FILE=file
;
  compile_opt idl2, strictarrsubs
;
@common
; http://www.wmo.ch/web/www/WDM/Guides/Guide-binary-2.html
;
;  gribfile =
;  '/d1fes2-raid6/SINTEX-common/ES10.d.00/atm/5d/ES10.d.00_5d_00911201_00911230.grib'
  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'
;
  openr, num, gribfile, /GET_LUN, ERROR = err, /SWAP_IF_LITTLE_ENDIAN
  if err ne 0 then begin
    ras = report(!err_string)
    return, -1
  ENDIF
;
  recstart = scan_grib_recstart(num)
;
;    messize = scan_grib_messize(num, recstart)
;    addoff = lonarr(n_elements(recstart))
;    FOR i = 1L, n_elements(recstart)-1 DO $
;      addoff[i] = recstart[i]-recstart[i-1]-messize[i-1]
;
;    nbits = scan_grib_nbits(num, recstart)
;    print,nbits[uniq(nbits,sort(nbits))]
;
  codes = scan_grib_code(num, recstart)
  nbcodes =  uniq(codes, sort(codes))
;
  dates = scan_grib_date(num, recstart)
  nbdates = uniq(dates, sort(dates))
;
  goodvar = where(codes EQ varcode)
  IF goodvar[0] EQ -1 THEN BEGIN
    ras = report( 'no var code '+strtrim(varcode, 2)+' in the file')
    return, -1
  ENDIF
;
  recstart = recstart[goodvar]
  dates = dates[goodvar]
;
  gooddate = where(dates GE date1 AND dates LE date2)
  IF gooddate[0] EQ -1 THEN BEGIN
    ras = report( 'no dates between '+strtrim(date1, 2)+' and '+strtrim(date2, 2)+' in the file')
    return, -1
  ENDIF
  recstart = recstart[gooddate]
  dates = dates[gooddate]
  key_caltype = '360d'
  time = date2jul(dates)
  jpt = n_elements(time)
  IF jpt EQ 1 THEN vardate = strtrim(dates[0], 2) ELSE vardate = strtrim(dates[0], 2)+' - '+strtrim(dates[jpt-1], 2)
;  varname
  vargrid = 'T'
;  varexp
;  varunit
;
  grib_pds = read_grib_pds(num, recstart[0])
; grid parameters
  IF grib_pds.gdsnotomitted THEN BEGIN
    grib_gds = read_grib_gds(num, recstart[0])
; min, max of the latitude with a precision of 10^-2
    lat1 = fix(100*grib_gds.la1)/100.
    lat2 = fix(100*grib_gds.la2)/100.
;
    CASE grib_gds.gridtype OF
; Latitude/Longitude Grid
      0:BEGIN
        computegrid, grib_gds.lo1, grib_gds.la1, grib_gds.di, -grib_gds.dj $
          , grib_gds.ni, grib_gds.nj
       END
; Gaussian Latitude/Longitude Grid
      4:BEGIN
; find the latitude axis
        CASE 1 OF
; n48
          grib_gds.n EQ 48 AND lat1 EQ 88.57 AND lat2 EQ -88.57:$
            gphit = n48gaussian()
; n80
           grib_gds.n EQ 80 AND lat1 EQ 89.14 AND lat2 EQ -89.14:$
             gphit = n80gaussian()
; n128
           grib_gds.n EQ 128 AND lat1 EQ 89.46 AND lat2 EQ -89.46:$
             gphit = n128gaussian()
; n160
           grib_gds.n EQ 160 AND lat1 EQ 89.57 AND lat2 EQ -89.57:$
             gphit = n160gaussian()
; n256
           grib_gds.n EQ 256 AND lat1 EQ 89.73 AND lat2 EQ -89.73:$
             gphit = n256gaussian()
; part of one of the gaussian grids defined above
           ELSE:BEGIN
             cnt = 0
             REPEAT BEGIN
               CASE cnt OF
                 0:gphit = n48gaussian()
                 1:gphit = n80gaussian()
                 2:gphit = n128gaussian()
                 3:gphit = n160gaussian()
                 4:gphit = n256gaussian()
                 5:BEGIN
                   gphit = n80gaussian()
                   lat1 = 29.71
                   lat2 = -19.62
                 END
                 ELSE:stop
               ENDCASE
               nfix = fix(gphit*100)/100.
               nlat1 = (where(nfix EQ lat1))[0]
               nlat2 = (where(nfix EQ lat2))[0]
               IF nlat1 NE -1 AND  nlat2 NE -1 $
                 AND nlat2-nlat1+1 EQ grib_gds.nj $
                 THEN gphit = gphit[nlat1:nlat2] ELSE gphit = -1
               cnt = cnt+1
             ENDREP UNTIL gphit[0] NE -1
           END
         ENDCASE
         computegrid, grib_gds.lo1, -1, grib_gds.di, -1, grib_gds.ni, -1, YAXIS = gphit
       END
; Mercator Projection Grid
       gridtype EQ 1:
; Gnomonic Projection Grid
       gridtype EQ 2:
; Lambert Conformal, secant or tangent, conical or bipolar (normal or
; oblique) Projection Grid
       gridtype EQ 3:
; Polar Stereographic Projection Grid
       gridtype EQ 5:
; Oblique Lambert conformal, secant or tangent, conical or bipolar,
; projection
       gridtype EQ 13:
; Spherical Harmonic Coefficients
       gridtype EQ 50:
; Space view perspective or orthographic grid
       gridtype EQ 90:
; reserved - see Manual on Codes
       ELSE:
     ENDCASE
   ENDIF ELSE stop
;
   res = fltarr(grib_gds.ni, grib_gds.nj, n_elements(recstart))
   FOR i = 0, n_elements(recstart)-1 DO BEGIN
     res[*, *, i] = read_grib_bds(num, recstart[i], grib_gds.ni, grib_gds.nj)
   ENDFOR
;
   free_lun, num
;
   IF keyword_set(key_yreverse) THEN res = reverse(res, 2)

   RETURN, res
 END