;+
;
; @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
;
; @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. Optional, if not specified date2=date1
;
; @keyword FILENAME
; the grads control file name: 'xxxx.ctl'
;
; @keyword TIMESTEP
; to specify that the dates are time steps instead of true calendar.
;
; NOT yet available
;
; @keyword BOX {type=A 4 or 6 elements 1d array, [lon1,lon2,lat1,lat2, depth1, depth2]}
; It specifies the area where data must be read
;
; @keyword EVERYTHING
;
; @keyword NOSTRUCT
;
; @keyword _EXTRA
; Used to pass keywords
;
; @returns
; an array
;
; @uses
; common
;
; @restrictions
; define all the grid parameters (defined in common)
; associated to the data.
;
; 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)
;
; @version
; $Id$
;
;-
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
ras = report( '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
ras = report('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
ras = report('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
ras = report(var+' not found in the variable list 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 file 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 existence 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
ras = report(!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
ras = report(['According to '+filename+' the file size must be '+strtrim(filesize, 1)+' instead of '+strtrim(infof2read.size, 1), $
'jpi: '+strtrim(jpi, 2), $
'jpj: '+strtrim(jpj, 2), $
'jpt: '+strtrim(jpt, 2), $
'format size in byte: '+strtrim(fmtsz, 2), $
'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 extraction
;------------------------
; time arguments
;------------------------
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