;+
;
; @file_comments
;
; @categories
;
; @param NAMEFILE
;
; @keyword GRID {default='T'}{type=scalar string}
; Used to specify on which grid type are located the data
;
; @keyword _EXTRA
; Used to pass keywords to isafile
; and ncdf_getaxis
;
; @returns
;
; @uses
;
; @restrictions
;
; @examples
;
; @history
;
; @version
; $Id$
;
; @todo
; seb : I don't know what to do with that...
;
;-
;
; liste des presupposes:
; 1) le fichier a lire est un fichier netcdf
; 2) le nom de ce fichier finit
; par U.nc, V.nc, W.nc, T.nc ou F.nc, la lettre avant le
; nc designant la grille a laquelle se rapporte la champ.
; Si tel n''est pas la cas, le fichier est attribue a la grille
; T.
; 3) ce fichier contient une dimension infinie qui doit etre
; celle qui se rapporte au temps et au mois 2 autres dimensions
; dont les noms sont 'x','lon...','xi_...' et 'y','lat...' ou
; 'eta_...' ou bien en majuscule.
; 4) il doit exister ds ce fichier une unique variable n''ayant
; qu''une dimension et etant la dimension temporelle. cette
; variable sera prise comme axe des temps. Rq: si plusieurs
; variables verifient ces criteres on considere la premiere
; variable
; 5) Cette variable axe des temps doit contenir l''attribut
; 'units' qui doit etre ecrit suivant la syntaxe:
; "seconds since 0001-01-01 00:00:00"
; "hours since 0001-01-01 00:00:00"
; "days since 1979-01-01 00:59:59"
; "months since 1979-01-01 00:59:59"
; "years since 1979-01-01 00:59:59"
;
; je crois que c''est tout!
;
; GRID='[UTVWF]' to specify the type of grid. Defaut is (1)
; based on the name of the file if the file ends by
; GRID[._][TUVFW].NC (not case sensible) or (2) T if case (1)
; is not found.
;
;
FUNCTION scanfile, namefile, GRID = GRID, _EXTRA = ex
;
compile_opt idl2, strictarrsubs
;
@common
;------------------------------------------------------------
; filename
;------------------------------------------------------------
fullname = isafile(filename = namefile, IODIRECTORY = iodir, _extra = ex)
IF size(fullname, /type) NE 7 THEN return, -1
;------------------------------------------------------------
; open file
;------------------------------------------------------------
cdfid = ncdf_open(fullname)
;------------------------------------------------------------
; What contains the file?
;------------------------------------------------------------
inside = ncdf_inquire(cdfid) ;
;------------------------------------------------------------
; name of all dimensions
;------------------------------------------------------------
namedim = strarr(inside.ndims)
for dimiq = 0, inside.ndims-1 do begin
ncdf_diminq, cdfid, dimiq, tmpname, value
namedim[dimiq] = strlowcase(tmpname)
ENDFOR
;------------------------------------------------------------
; x/y dimensions id
;------------------------------------------------------------
ncdf_getaxis, cdfid, dimidx, dimidy, _extra = ex
;------------------------------------------------------------
; name of all variables
;------------------------------------------------------------
; we keep only the variables containing at least x, y and time dimension (if existing...)
namevar = strarr(inside.nvars)
for varid = 0, inside.nvars-1 do begin
invar = ncdf_varinq(cdfid, varid) ; what contains the variable?
if (inter(invar.dim, dimidx))[0] NE -1 AND $
(inter(invar.dim, dimidy))[0] NE -1 AND $
((where(invar.dim EQ inside.recdim))[0] NE -1 OR inside.recdim EQ -1) $
THEN namevar[varid] = invar.name
ENDFOR
namevar = namevar[where(namevar NE '')]
;------------------------------------------------------------
; find vargrid for each selected variable...
;------------------------------------------------------------
listgrid = strarr(n_elements(namevar))
; default definitions
IF keyword_set(grid) THEN vargrid = strupcase(grid) ELSE vargrid = 'T'
; look for values of vargrid for each variable
IF finite(glamu[0]) EQ 1 AND NOT keyword_set(grid) THEN BEGIN
; for each variable, look if we in one of the case corresponding to ROMS conventions?
FOR i = 0, n_elements(namevar)-1 do begin
invar = ncdf_varinq(cdfid, namevar[i])
tmpnm = namedim[invar.dim]
; are we in one of the case corresponding to ROMS conventions?
CASE 1 OF
tmpnm[2 <(invar.ndims-1)] EQ 's_w':vargrid = 'W'
tmpnm[0] EQ 'xi_rho' AND tmpnm[1] EQ 'eta_rho':listgrid[i] = 'T'
tmpnm[0] EQ 'xi_u' AND tmpnm[1] EQ 'eta_u' :listgrid[i] = 'U'
tmpnm[0] EQ 'xi_v' AND tmpnm[1] EQ 'eta_v' :listgrid[i] = 'V'
tmpnm[0] EQ 'xi_psi' AND tmpnm[1] EQ 'eta_psi':listgrid[i] = 'F'
tmpnm[0] EQ 'xi_rho' AND tmpnm[1] EQ 'eta_v' :listgrid[i] = 'V'
tmpnm[0] EQ 'xi_u' AND tmpnm[1] EQ 'eta_rho':listgrid[i] = 'U'
tmpnm[0] EQ 'xi_u' AND tmpnm[1] EQ 'eta_v' :listgrid[i] = 'F'
ELSE:
ENDCASE
ENDFOR
empty = where(listgrid EQ '')
IF empty[0] NE -1 THEN BEGIN
; could we define the grid type from the file name??
pattern = ['GRID.', 'GRID_', 'GRID', 'UPID_', '30ID_']
gdtype = ['T', 'U', 'V', 'W', 'F']
fnametest = strupcase(fullname)
FOR i = 0, n_elements(pattern)-1 DO BEGIN
FOR j = 0, n_elements(gdtype)-1 DO BEGIN
substr = pattern[i]+gdtype[j]
pos = strpos(fnametest, substr)
IF pos NE -1 THEN $
vargrid = strmid(fnametest, pos+strlen(substr)-1, 1)
ENDFOR
ENDFOR
listgrid[empty] = vargrid
ENDIF
ENDIF ELSE listgrid[*] = vargrid
;------------------------------------------------------------
; time axis
;------------------------------------------------------------
date0fk = date2jul(19000101)
IF inside.recdim EQ -1 THEN BEGIN
jpt = 1
time = date0fk
fakecal = 1
ENDIF ELSE BEGIN
ncdf_diminq, cdfid, inside.recdim, timedimname, jpt
; we look for the variable containing the time axis
; we look for the first variable having for only dimension inside.recdim
varid = 0
repeat BEGIN
IF varid LT inside.nvars THEN BEGIN
invar = ncdf_varinq(cdfid, varid)
varid = varid+1
ENDIF ELSE varid = 0
endrep until (n_elements(invar.dim) EQ 1 AND invar.dim[0] EQ inside.recdim) OR (varid EQ 0)
varid = varid-1
;
CASE 1 OF
varid EQ -1:BEGIN
dummy = report('the file '+fullname+' has no time axis.!C we create a fake calendar ...')
fakecal = 1
time = date0fk + lindgen(1>jpt)
END
invar.natts EQ 0:BEGIN
dummy = report('the variable '+invar.name+' has no attribut.!C we create a fake calendar ...')
fakecal = 1
time = date0fk + lindgen(1>jpt)
END
ELSE:BEGIN
;
; we want to know which attributes are attached to the time variable...
;
attnames = strarr(invar.natts)
for attiq = 0, invar.natts-1 do attnames[attiq] = ncdf_attname(cdfid, varid, attiq)
if (where(attnames EQ 'units'))[0] EQ -1 then BEGIN
dummy = report('Attribut ''units'' not found for the variable '+invar.name+'!C we create a fake calendar ...')
fakecal = 1
time = date0fk + lindgen(1>jpt)
ENDIF ELSE BEGIN
; we read the time axis
ncdf_varget, cdfid, varid, time
time = double(time)
ncdf_attget, cdfid, varid, 'units', value
; time_counter:units = "seconds since 0001-01-01 00:00:00" ;
; time_counter:units = "hours since 0001-01-01 00:00:00" ;
; time_counter:units = "days since 1979-01-01 00:00:00" ;
; time_counter:units = "months since 1979-01-01 00:00:00" ;
; time_counter:units = "years since 1979-01-01 00:00:00" ;
value = strtrim(strcompress(string(value)), 2)
mots = str_sep(value, ' ')
unite = mots[0]
unite = strlowcase(unite)
IF strpos(unite, 's', strlen(unite)-1) NE -1 THEN unite = strmid(unite, 0, strlen(unite)-1)
IF strpos(unite, 'julian_') NE -1 THEN unite = strmid(unite, 7)
err = 0
IF unite NE 'second' AND unite NE 'hour' AND unite NE 'day' $
AND unite NE 'month' AND unite NE 'year' THEN BEGIN
dummy = report('time units does not start with seconds/hours/days/months/years')
err = 1
ENDIF
IF stregex(value, '[^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*', /boolean) EQ 0 THEN BEGIN
dummy = report('attribut units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*')
err = 1
ENDIF
IF err GT 0 THEN BEGIN
fakecal = 1
time = date0fk + lindgen(1>jpt)
ENDIF ELSE BEGIN
debut = str_sep(mots[2], '-')
;
; now we try to find the attribut called calendar...
; the attribute "calendar" exists?
; If no, we suppose that the calendar is gregorian calendar
;
if (where(attnames EQ 'calendar'))[0] NE -1 then BEGIN
ncdf_attget, cdfid, varid, 'calendar', value
value = string(value)
CASE value OF
'noleap':key_caltype = 'noleap'
'360d':key_caltype = '360d'
'greg':IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg'
ELSE:BEGIN
; notused = report('Unknown calendar: '+value+', we use greg calendar.')
key_caltype = 'greg'
END
ENDCASE
ENDIF ELSE BEGIN
; notused = report('Unknown calendar, we use '+key_caltype+' calendar.')
IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg'
ENDELSE
;
; BEWARE we have to get back the calendar attribute and ajust time by consequence...
;
;
; We pass time in IDL julian days
;
case unite of
'second':time = julday(debut[1], debut[2], debut[0])+time/86400.d
'hour':time = julday(debut[1], debut[2], debut[0])+time/24.d
'day':time = julday(debut[1], debut[2], debut[0])+time
'month':BEGIN
if total(fix(time) NE time) NE 0 then $ ; we switch to days with 30d/m
time = julday(debut[1], debut[2], debut[0])+round(time*30) $
ELSE for t = 0, n_elements(time)-1 DO $
time[t] = julday(debut[1]+time[t], debut[2], debut[0])
END
'year':BEGIN
if total(fix(time) NE time) NE 0 then $ ; we switch to days with 365d/y
time = julday(debut[1], debut[2], debut[0])+round(time*365) $
ELSE for t = 0, n_elements(time)-1 do $
time[t] = julday(debut[1], debut[2], debut[0]+time[t])
END
ENDCASE
;
; high frequency calendar: more than one element per day
IF max(histogram([long(time-time[0])])) GT 1 THEN fakecal = 1 ELSE fakecal = 0
date0fk = date2jul(19000101)
IF keyword_set(fakecal) THEN time = date0fk+lindgen(1>jpt) $
ELSE time = long(time)
;
ENDELSE
ENDELSE
END
ENDCASE
ENDELSE
;------------------------------------------------------------
ncdf_close, cdfid
;------------------------------------------------------------
return, {filename:fullname, time_counter:time, listvar:namevar $
, listgrid:strupcase(listgrid), caltype:key_caltype $
, fakecal:date0fk*fakecal}
end