source: trunk/SRC/ReadWrite/read_grads.pro @ 121

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; @file_comments reading grads file (except "data type station" or "grib")
; from the grads control file even if there is multiple data files.
;
; @categories reading function
;
;       @param var {in}{required} the variable name
;       @param date1 {in}{required} date of the beginning (yyyymmdd if TIMESTEP is not activate)
;       @param date2 {in}{optional} last date. Optionnal, if not scpecified date2=date1
;      
; @keyword FILENAME the grads control file name: 'xxxx.ctl'
;
; @file_comments keyword GLAMBOUNDARY (via computegrid.pro) a 2 elements vector,
;       {lon1,lon2], giving the longitude boundaries that should be
;       used to visualize the data.
;         lon2 > lon1
;         lon2 - lon1 eq 360
;       key_shift will be automatically defined according to
;       GLAMBOUNDARY.
;
;       @keyword /TIMESTEP to specify that the dates are time steps instead of 
;       true calendar.
;
;        @file_comments keyword IODIRECTORY a string giving the name of iodirectory (see
;       isafile.pro for all possibilities). default value is common 
;       variable iodir
;
; @file_comments
;---------------
; NOT yet available
;---------------
;
;       @hidden BOX a 4 or 6 elements 1d array, [lon1,lon2,lat1,lat2, depth1, 
;       depth2], that specifies the area where data must be read
;
;       @hidden EVERYTHING
; 
;       @hidden NOSTRUCTURE
;       
; @returns an array
;
; @uses common
;
; @restriction define all the grid parameters (defined in common.pro)
; associated to the data.
;
; @restrictions this function call the procedure scanfile that use the 
; unix commands grep and sed 
;
; @examples
;    IDL> a=read_grads('sst',19900101,19900131,filename='outputs.ctl')
;    IDL> plt, a
;
; @history Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------

FUNCTION read_grads, var, date1, date2, FILENAME = filename, BOX=box, TIMESTEP = timestep, EVERYTHING = everything, NOSTRUCT = nostruct, _EXTRA = ex
;---------------------------------------------------------
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
@cm_4cal
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
  ENDIF
;------------------------------------------------------------
; we find the filename.
;------------------------------------------------------------
   filename = isafile(FILENAME = filename, IODIRECTORY = iodir, _EXTRA = ex)
   if size(filename, /type) NE 7 then $
    return, report('read_ncdf cancelled')
;------------------------------------------------------------
; we scan the control file called filename
;------------------------------------------------------------
   scanctl, filename, filesname, jpt1file, varsname, varslev, swapbytes, bigendian, littleendian, f77sequential, fileheader, theader, xyheader, VARFMT = varfmt, _EXTRA = ex
   if n_elements(varfmt) EQ 0 then varfmt = 'float'
;------------------------
;------------------------
; check date1 and date2 and found the starting index "t1" and the
; ending index "t2" that corresponds to the time series specified by
; date1 and date2 for the time axis defined in the .ctl file.
;------------------------
;------------------------
   if n_elements(date1) EQ 0 then begin
      t0 = 0
      t1 = 0
   ENDIF
   if n_elements(date2) EQ 0 then date2 = date1
   if keyword_set(timestep) then BEGIN
      if date1 GT date2 then begin
         print, 'date2 must be larger than date1'
         return, -1
      endif
      t1 = 0 > long(date1) < (jpt-1)
      t2 = t1 > long(date2) < (jpt-1)
   ENDIF ELSE BEGIN
      jdate1 = date2jul(date1, /grads) < time[jpt-1]
      jdate2 = time[0] > date2jul(date2, /grads)
      if jdate1 GT jdate2 then begin
         print, 'date2 must be larger than date1'
         return, -1
      endif
      t1 = (where(time GE jdate1))[0]
      tmp = where(time LE jdate2, t2)
      t2 = t2-1
   ENDELSE
   if t2 LT t1 then begin
      print, 'There is no date between date1 and date2'
      return, -1
   endif
   jpt2read = t2-t1+1
;------------------------
;------------------------
; index of the variable
;------------------------
;------------------------
   varid = where(strlowcase(varsname) EQ strlowcase(var))
   varid = varid[0]
   if varid EQ -1 then begin
      print, var+' not found in the variable liste of '+filename
      return,  -1
   ENDIF
   varname = var
   if varslev[varid] EQ 1 then res = fltarr(jpi, jpj, jpt2read, /nozero) $
   ELSE res = fltarr(jpi, jpj, varslev[varid], jpt2read, /nozero)
;------------------------
;------------------------
; find the first file to be read according to the lile list,, the
; number of time step in each file and t1 and t2
;------------------------
;------------------------
   indf2read = t1/jpt1file
   startread = t1-indf2read*jpt1file
   alreadyread = 0
readagain:
   jpt2read1file = min([jpt1file-startread, jpt2read])
   f2read = filesname[indf2read]
;------------------------
;------------------------
; opening
;------------------------
;------------------------;
; check the existance of the file
   f2read = isafile(filename = f2read, iodirectory = iodir, _EXTRA = ex)
; if the file is stored on tape
   if !version.os_family EQ 'unix' then spawn, 'file '+f2read+' > /dev/null'
; open the file
   openr, unit, f2read, /get_lun, error=err $
    , swap_if_little_endian = bigendian $
    , swap_if_big_endian = littleendian $
    , swap_endian = swapbytes
   if err ne 0 then begin
      print,!err_string
      return, -1
   endif
;------------------------
   case varfmt of
      'byte':fmtsz = 1l
      'uint':fmtsz = 2l
      'int':fmtsz = 2l
      'long':fmtsz = 4l
      'float':fmtsz = 4l
   endcase
;------------------------
; check its size
   addf77sec = long(4*2*f77sequential)
   xyblocsize = xyheader + addf77sec*(xyheader NE 0) + jpi*jpj*fmtsz +addf77sec
   nxybloc = long(total(varslev))
   filesize = (fileheader + addf77sec*(fileheader NE 0)) $
    +(theader+addf77sec*(theader NE 0) + nxybloc*xyblocsize)*jpt1file
   infof2read=fstat(unit)
   if infof2read.size NE filesize then begin
      print, 'According to '+filename+' the file size must be '+strtrim(filesize, 1)+' instead of '+strtrim(infof2read.size, 1)
      print, 'jpi: '+strtrim(jpi, 2)
      print, 'jpj: '+strtrim(jpj, 2)
      print, 'jpt: '+strtrim(jpt, 2)
      print, 'format size in byte: '+strtrim(fmtsz, 2)
      print, 'number of xy arrays: '+strtrim(nxybloc, 2)
      return, -1
   endif
;------------------------
;------------------------
; reading
;------------------------
;------------------------;
;------------------------;
; loop on the time steps to be read in one file
   for i = 0, jpt2read1file-1 do begin
; computing the offset
      offset = (fileheader + addf77sec*(fileheader NE 0)) $
       +(theader+addf77sec*(theader NE 0) + nxybloc*xyblocsize)*(startread+i) $
       +theader+addf77sec*(theader NE 0)
      if varid NE 0 THEN $
       offset = offset + long(total(varslev[0:varid-1]))*xyblocsize
; if there is only one level
      IF varslev[varid] EQ 1 then begin
         case varfmt of
            'byte':a=assoc(unit, bytarr(jpi,jpj,/nozero), offset+4*f77sequential)
            'uint':a=assoc(unit,uintarr(jpi,jpj,/nozero), offset+4*f77sequential)
            'int':a=assoc(unit,  intarr(jpi,jpj,/nozero), offset+4*f77sequential)
            'long':a=assoc(unit, lonarr(jpi,jpj,/nozero), offset+4*f77sequential)
            'float':a=assoc(unit,fltarr(jpi,jpj,/nozero), offset+4*f77sequential)
         endcase
         res[*, *, i+alreadyread]=a[0] 
      ENDIF ELSE BEGIN ; more than 1 level to be read
         if f77sequential then BEGIN ; sequential access
            case varfmt of
               'byte':a=assoc(unit, bytarr(jpi*jpj+8, varslev[varid],/nozero), offset)
               'uint':a=assoc(unit,uintarr(jpi*jpj+4, varslev[varid],/nozero), offset)
               'int':a=assoc(unit,  intarr(jpi*jpj+4, varslev[varid],/nozero), offset)
               'long':a=assoc(unit, lonarr(jpi*jpj+2, varslev[varid],/nozero), offset)
               'float':a=assoc(unit,fltarr(jpi*jpj+2, varslev[varid],/nozero), offset)
            endcase
            tmp = a[0]
            case varfmt OF ; we cut the headers and tailers of f77 write
               'byte': tmp = tmp[4:jpi*jpj+3, *]
               'uint': tmp = tmp[2:jpi*jpj+1, *]
               'int':  tmp = tmp[2:jpi*jpj+1, *]
               'long': tmp = tmp[1:jpi*jpj+0, *]
               'float':tmp = tmp[1:jpi*jpj+0, *]
            endcase
            if keyword_set(key_zreverse) then res[*, *, *, i+alreadyread]=reverse(reform(tmp,  jpi, jpj, varslev[varid], /over), 3) ELSE res[*, *, *, i+alreadyread]=reform(tmp,  jpi, jpj, varslev[varid], /over)
         ENDIF ELSE BEGIN  ; direct acces
            case varfmt of
               'byte':a=assoc(unit, bytarr(jpi,jpj,varslev[varid],/nozero),offset)
               'uint':a=assoc(unit,uintarr(jpi,jpj,varslev[varid],/nozero),offset)
               'int':a=assoc(unit,  intarr(jpi,jpj,varslev[varid],/nozero),offset)
               'long':a=assoc(unit, lonarr(jpi,jpj,varslev[varid],/nozero),offset)
               'float':a=assoc(unit,fltarr(jpi,jpj,varslev[varid],/nozero),offset)
            endcase
            if keyword_set(key_zreverse) then res[*, *, *, i+alreadyread]=reverse(a[0], 3) ELSE res[*, *, *, i+alreadyread]=a[0] 
         ENDELSE
      ENDELSE
   endfor
;------------------------
; close the file
   free_lun,unit
   close,unit
;------------------------
;------------------------
; do we need to read a new file to complete the time series ?
;------------------------
;------------------------
   if jpt2read1file NE jpt2read then BEGIN
      indf2read = indf2read+1
      startread = 0
      alreadyread = alreadyread+jpt2read1file
      jpt2read = jpt2read-jpt2read1file
      GOTO, readagain
   ENDIF
;------------------------
;------------------------
; post processing
;------------------------
;------------------------
   if keyword_set(key_yreverse) then res = reverse(res, 2)
   if keyword_set(key_shift) then begin
      case (size(res))[0] of
         2:res = shift(res, key_shift, 0)
         3:res = shift(res, key_shift, 0, 0)
         4:res = shift(res, key_shift, 0, 0, 0)
      endcase
   endif
;------------------------
;   mask
;   IF varslev[varid] EQ 1 then begin
;      if abs(valmask) LE 1e5 then notgood = where(res[*, *, 0] EQ valmask) $
;      ELSE notgood = where(abs(res[*, *, 0]) GE abs(valmask)/10)
;      if notgood[0] NE -1 then tmask[notgood] = 0b
;   ENDIF ELSE BEGIN
;      if abs(valmask) LE 1e5 then notgood = where(res[*, *, *, 0] EQ valmask) $
;      ELSE notgood = where(abs(res[*, *, *, 0]) GE abs(valmask)/10)
;      if notgood[0] NE -1 then tmask[notgood] = 0b
;   ENDELSE
   if abs(valmask) LE 1e5 then notgood = where(res EQ valmask) $
   ELSE notgood = where(abs(res) GE abs(valmask)/10)
   if notgood[0] NE -1 THEN res[notgood] = !values.f_nan
;   valmask = 1e20
;   if abs(valmask) LE 1e5 then notgood = where(res EQ valmask) $
;   ELSE notgood = where(abs(res) GE abs(valmask)/10)
;   if notgood[0] NE -1 THEN res[notgood] = 1e20
;   valmask = 1e20
   triangles_list = triangule()
;------------------------
; subdomain extration

;------------------------
; time aguments
;------------------------
   time = time[t1:t2]
   jpt = t2-t1+1
   if keyword_set(timestep) then vardate = strtrim(time[0], 2) $
   ELSE vardate = date2string(vairdate(time[0]))
;------------------------
   @updateold
;------------------------

   return, res
end