;+ ; @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 your 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 exiter 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 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