;+
;
; @file_comments
; Reading function for the file net_cdf.
; This program is less universal than ncdf_lec (it appeal to declared
; variables in common.pro) but it is very easier to be used. It considerate
; the declaration of the different zooms which have been defined
; (ixminmesh...premierx...), the declaration of the variable key_shift...
; To put it in a nutshell, the result of read_ncdf can be directly used in plt...
; This is also this program which is used by default in our reading widgets.
;
; @categories
; Reading
;
; @param NAME {in}{required}{type=string}
; It define the field to be read.
;
; @param BEGINNING {in}{required}
; Relative with the time axis.
; These can be
; - 2 date of the type yyyymmdd and in this case, we select dates
; which are included between these two dates.
; - 2 indexes which define between which and which time step we have
; to extract the temporal dimension.
;
; @param ENDING {in}{required}
; Relative with the time axis.
; See BEGINNING.
;
; @param COMPATIBILITY {in}{optional}
; Useless, defined for compatibility
;
; @keyword ADDSCL_BEFORE {default=0}{type=scalar: 0 or 1}
; put 1 to apply add_offset ad scale factor on data before looking for
; missing values
;
; @keyword BOXZOOM
; Contain the boxzoom on which we have to do the reading
;
; @keyword CALLITSELF {default=0}{type=scalar: 0 or 1}
; For ROMS outputs. Use by read_ncdf itself to access auxilliary data (h and zeta).
;
; @keyword DIREC
; a string used to specify the direction along which we want to make
; spatial and/or temporal mean. It could be: 'x' 'y' 'z' 't' 'xy' 'xz'
; 'yz' 'xyz' 'xt' 'yt' 'zt' 'xyt' 'xzt' 'yzt' or 'xyzt'
;
; @keyword FILENAME {required}{type=string}
; It contains the file's name.
;
; @keyword INIT {default=0}{type=scalar: 0 or 1}
; To call automatically initncdf with filename as input argument
; and thus ; redefine all the grid parameters
;
; @keyword GRID
; ='[UTVWF]' to specify the type of grid. Default 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.
;
; @keyword TIMESTEP {default=0}{type=scalar: 0 or 1}
; Specify that BEGINNING and ENDING refer to indexes of the time axis and not to dates
;
; @keyword TOUT {default=0}{type=scalar: 0 or 1}
; We activate it if we want to read the file on the whole domain without
; considerate the sub-domain defined by the boxzoom or
; lon1,lon2,lat1,lat2,vert1,vert2.
;
; @keyword NOSTRUCT {default=0}{type=scalar: 0 or 1}
; We activate it if we do not want that read_ncdf send back a structure
; but only the array referring to the field.
;
; @keyword ZETAFILENAME {default=FILENAME}{type=string}
; For ROMS outputs. The filename of the file where zeta variable should be read
;
; @keyword ZETAZERO {default=0}{type=scalar: 0 or 1}
; For ROMS outputs. To define zeta to 0. instead of reading it
;
; @keyword _EXTRA
; Used to pass keywords to isafile, initncdf,
; ncdf_gettime and domdef
;
; @returns
; Structure readable by litchamp or an array if NOSTRUCT is activated.
; @uses
; common.pro
;
; @restrictions
; The field must have a temporal dimension.
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
; 15/10/1999
;
; @version
; $Id$
;-
;
FUNCTION read_ncdf, name, beginning, ending, compatibility, BOXZOOM = boxzoom, FILENAME = filename $
, PARENTIN = parentin, TIMESTEP = timestep, ADDSCL_BEFORE = addscl_before $
, TOUT = tout, NOSTRUCT = nostruct, CONT_NOFILL = CONT_NOFILL, INIT = init $
, GRID = grid, CALLITSELF = callitself, DIREC = direc $
, ZETAFILENAME = zetafilename, ZETAZERO = zetazero $
, _EXTRA = ex
;
compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
@cm_4cal
IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
@updatekwd
ENDIF
;------------------------------------------------------------
; we find the filename.
;------------------------------------------------------------
; print,filename
; is parent a valid widget ?
IF keyword_set(parentin) THEN BEGIN
parent = long(parentin)
parent = parent*widget_info(parent, /managed)
ENDIF
filename = isafile(filename = filename, IODIRECTORY = iodir, _EXTRA = ex)
;------------------------------------------------------------
; Opening of the name file
;------------------------------------------------------------
IF size(filename, /type) NE 7 THEN $
return, report('read_ncdf cancelled')
IF !version.OS_FAMILY EQ 'unix' THEN spawn, '\file '+filename+' > /dev/null'
cdfid = ncdf_open(filename)
inq = ncdf_inquire(cdfid)
;------------------------------------------------------------
; we check if the variable name exists in the file.
;------------------------------------------------------------
IF ncdf_varid(cdfid, name) EQ -1 THEN BEGIN
ncdf_close, cdfid
return, report('variable '+name+' !C not found in the file '+filename)
ENDIF
varinq = ncdf_varinq(cdfid, name)
IF varinq.ndims LT 2 THEN return, report('read_ncdf cannot read scalar or 1D data')
; look for the dimension names
dimnames = strarr(varinq.ndims)
FOR i = 0, varinq.ndims-1 DO BEGIN
ncdf_diminq, cdfid, varinq.dim[i], tmp, dimsize
dimnames[i] = tmp
ENDFOR
;------------------------------------------------------------
; shall we redefine the grid parameters
;------------------------------------------------------------
IF keyword_set(init) THEN initncdf, filename, _extra = ex
;------------------------------------------------------------
; check the time axis and the debut and ending dates
;------------------------------------------------------------
IF n_elements(beginning) EQ 0 THEN BEGIN
beginning = 0L
timestep = 1L
ENDIF
; define time and jpt
CASE 1 OF
keyword_set(timestep):BEGIN
firsttps = long(beginning[0])
IF n_elements(ending) NE 0 THEN lasttps = long(ending[0]) ELSE lasttps = firsttps
jpt = lasttps-firsttps+1
IF NOT keyword_set(callitself) then time = julday(1, 1, 1) + lindgen(jpt)
END
keyword_set(parent):BEGIN
widget_control, parent, get_uvalue = top_uvalue
filelist = extractatt(top_uvalue, 'filelist')
IF filelist[0] EQ 'many !' THEN filelist = filename
currentfile = (where(filelist EQ filename))[0]
time = (*(extractatt(top_uvalue, 'fileparameters'))[currentfile]).time_counter
date1 = date2jul(beginning[0])
if n_elements(ending) NE 0 then date2 = date2jul(ending[0]) ELSE date2 = date1
firsttps = (where(abs(time - date1) LT 0.9d/86400.d))[0] ; beware of rounding errors...
lasttps = (where(abs(time - date2) LT 0.9d/86400.d))[0]
jpt = lasttps-firsttps+1
END
ELSE:BEGIN
time = ncdf_gettime(filename, cdfid, caller = 'read_ncdf', err = err, _extra = ex)
IF n_elements(err) NE 0 THEN BEGIN
dummy = report(err)
ncdf_close, cdfid
return, -1
ENDIF
; date1
date1 = date2jul(beginning[0])
; date2
if n_elements(ending) NE 0 then date2 = date2jul(ending[0]) ELSE date2 = date1
; firsttps
firsttps = where(time GE (date1 - 0.9d/86400.d)) & firsttps = firsttps[0]
if firsttps EQ -1 THEN BEGIN
ncdf_close, cdfid
return, report('date 1: '+strtrim(jul2date(date1), 1)+' is not found in the time axis.')
ENDIF
; lasttps
lasttps = where(time LE (date2 + 0.9d/86400.d)) & lasttps = lasttps[n_elements(lasttps)-1]
if lasttps EQ -1 THEN BEGIN
ncdf_close, cdfid
return, report('the time axis has no date before date 2: '+strtrim(jul2date(date2), 1))
endif
if lasttps LT firsttps then BEGIN
ncdf_close, cdfid
return, report('the time axis has no dates between date1 and date 2: '+strtrim(jul2date(date1), 1)+' '+strtrim(jul2date(date2), 1))
endif
time = time[firsttps:lasttps]
jpt = lasttps-firsttps+1
END
ENDCASE
;------------------------------------------------------------
; Name of the grid on which the field refer to.
;------------------------------------------------------------
IF keyword_set(grid) THEN vargrid = strupcase(grid) ELSE BEGIN
vargrid = 'T' ; default definition
IF finite(glamu[0]) EQ 1 THEN BEGIN
; are we in one of the case corresponding to ROMS conventions?
CASE 1 OF
dimnames[2 <(varinq.ndims-1)] EQ 's_w':vargrid = 'W'
dimnames[0] EQ 'xi_rho' AND dimnames[1] EQ 'eta_rho':vargrid = 'T'
dimnames[0] EQ 'xi_u' AND dimnames[1] EQ 'eta_u' :vargrid = 'U'
dimnames[0] EQ 'xi_v' AND dimnames[1] EQ 'eta_v' :vargrid = 'V'
dimnames[0] EQ 'xi_psi' AND dimnames[1] EQ 'eta_psi':vargrid = 'F'
dimnames[0] EQ 'xi_rho' AND dimnames[1] EQ 'eta_v' :vargrid = 'V'
dimnames[0] EQ 'xi_u' AND dimnames[1] EQ 'eta_rho':vargrid = 'U'
dimnames[0] EQ 'xi_u' AND dimnames[1] EQ 'eta_v' :vargrid = 'F'
ELSE: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(filename)
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
END
ENDCASE
ENDIF
ENDELSE
;---------------------------------------------------------------
; redefinition of the domain
;---------------------------------------------------------------
IF keyword_set(tout) THEN BEGIN
nx = jpi
ny = jpj
nz = jpk
firstx = 0
firsty = 0
firstz = 0
lastx = jpi-1
lasty = jpj-1
lastz = jpk-1
case strupcase(vargrid) of
'T':mask = tmask
'U':mask = umask()
'V':mask = vmask()
'W':mask = tmask
'F':mask = fmask()
endcase
ENDIF ELSE BEGIN
if keyword_set(boxzoom) then BEGIN
Case 1 Of
N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]]
N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]]
N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2]
N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]]
N_Elements(Boxzoom) Eq 6:bte = Boxzoom
Else: BEGIN
ncdf_close, cdfid
return, report('Wrong Definition of Boxzoom')
end
ENDCASE
savedbox = 1b
saveboxparam, 'boxparam4rdncdf.dat'
domdef, bte, GRIDTYPE = ['T', vargrid], _extra = ex
ENDIF
grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz
undefine, glam & undefine, gphi & ; We liberate some memory!
ENDELSE
;---------------------------------------------------------------------
; We initialize ixmindta, iymindta if needed
;---------------------------------------------------------------------
if n_elements(jpidta) EQ 0 THEN jpidta = jpiglo
if n_elements(jpjdta) EQ 0 THEN jpjdta = jpjglo
if n_elements(jpkdta) EQ 0 THEN jpkdta = jpkglo
if n_elements(ixmindta) EQ 0 THEN ixmindta = 0
if n_elements(ixmaxdta) EQ 0 then ixmaxdta = jpidta-1
if ixmindta EQ -1 THEN ixmindta = 0
IF ixmaxdta EQ -1 then ixmaxdta = jpidta-1
if n_elements(iymindta) EQ 0 THEN iymindta = 0
IF n_elements(iymaxdta) EQ 0 then iymaxdta = jpjdta-1
if iymindta EQ -1 THEN iymindta = 0
IF iymaxdta EQ -1 then iymaxdta = jpjdta-1
if n_elements(izmindta) EQ 0 THEN izmindta = 0
IF n_elements(izmaxdta) EQ 0 then izmaxdta = jpkdta-1
if izmindta EQ -1 THEN izmindta = 0
IF izmaxdta EQ -1 then izmaxdta = jpkdta-1
;----------------------------------------------------------------
; We will read the file
;---------------------------------------------------------------
if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1]
key_stride = 1l > long(key_stride)
;---------------------------------------------------------------------
;---------------------------------------------------------------------
@read_ncdf_varget
;---------------------------------------------------------------------
;---------------------------------------------------------------------
; We define common (cm_4data) variables associated with the variable.
;---------------------------------------------------------------------
; varname
IF NOT keyword_set(callitself) THEN varname = name
; varunit and others...
ncdf_getatt, cdfid, name, add_offset = add_offset, scale_factor = scale_factor, missing_value = missing_value, units = units
IF NOT keyword_set(callitself) THEN varunit = units
; vardate
; We make a legible date in function of the specified language.
year = long(beginning[0])/10000
month = (long(beginning[0])/100) MOD 100
day = (long(beginning[0]) MOD 100)
vardate = string(format = '(C(CMoA))', 31*(month-1))+' '+strtrim(day, 1)+', '+strtrim(year, 1)
varexp = file_basename(filename)
; We apply the value valmask on land points.
if NOT keyword_set(cont_nofill) then begin
valmask = 1.e20
case 1 of
varinq.ndims eq 2: earth = where(mask[*, *, 0] EQ 0) ;xy array
varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] EQ -1:earth = where(mask EQ 0) ;xyz array
varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] NE -1:earth = where(mask[*, *, 0] EQ 0) ;xyt array
varinq.ndims eq 4: earth = where(mask EQ 0) ;xyzt array
ENDCASE
ENDIF ELSE earth = -1
; we look for missing_value
; we apply add_offset, scale_factor and missing_value
IF keyword_set(addscl_before) THEN BEGIN
IF scale_factor NE 1 THEN res = temporary(res)*scale_factor
IF add_offset NE 0 THEN res = temporary(res)+add_offset
ENDIF
IF size(missing_value, /type) NE 7 THEN BEGIN
IF finite(missing_value) EQ 1 THEN BEGIN
CASE 1 OF
missing_value GT 1.e6:missing = where(res GT missing_value/10.)
missing_value LT -1.e6:missing = where(res LT missing_value/10.)
abs(missing_value) LT 1.e-6:missing = where(res LT 1.e-6)
ELSE:missing = where(res EQ missing_value)
ENDCASE
ENDIF ELSE missing = where(finite(res) EQ 0)
ENDIF ELSE missing = -1
IF NOT keyword_set(addscl_before) THEN BEGIN
if scale_factor NE 1 then res = temporary(res)*scale_factor
if add_offset NE 0 then res = temporary(res)+add_offset
ENDIF
IF missing[0] NE -1 THEN res[temporary(missing)] = !values.f_nan
IF earth[0] NE -1 THEN BEGIN
IF varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] NE -1 THEN $ ;xyt array
earth = earth#replicate(1, jpt)+replicate(nx*ny, n_elements(earth))#lindgen(jpt)
IF varinq.ndims eq 4 THEN earth = earth#replicate(1, jpt)+replicate(nx*ny*nz, n_elements(earth))#lindgen(jpt)
res[temporary(earth)] = 1.e20
ENDIF
;---------------------------------------------------------------------
; if it is roms outputs, we need to get additionals infos...
IF NOT keyword_set(callitself) THEN BEGIN
IF strmid(dimnames[0], 0, 3) EQ 'xi_' AND strmid(dimnames[1], 0, 4) EQ 'eta_' THEN BEGIN
ncdf_attget, cdfid, 'theta_s', theta_s, /global
ncdf_attget, cdfid, 'theta_b', theta_b, /global
ncdf_attget, cdfid, 'hc', hc, /global
; look for all variables names
allvarnames = strarr(inq.nvars)
FOR i = 0, inq.nvars-1 DO BEGIN
tmp = ncdf_varinq( cdfid, i)
allvarnames[i] = tmp.name
ENDFOR
CASE 1 OF
keyword_set(zetazero):zeta = fltarr(nx, ny, jpt)
keyword_set(zetafilename): $
zeta = read_ncdf('zeta', firsttps, lasttps, FILENAME = zetafilename $
, /TIMESTEP, /NOSTRUCT, CONT_NOFILL = CONT_NOFILL $
, GRID = vargrid, /CALLITSELF, _EXTRA = ex)
(where(allvarnames EQ 'zeta'))[0] NE -1: $
zeta = read_ncdf('zeta', firsttps, lasttps, FILENAME = filename $
, /TIMESTEP, /NOSTRUCT, CONT_NOFILL = CONT_NOFILL $
, GRID = vargrid, /CALLITSELF, _EXTRA = ex)
ELSE:return, report('The variable zeta was not found in the file, please use the keyword ZETAFILENAME to specify the name of a file containing zeta or use keyword ZETAZERO to define zeta to 0.')
ENDCASE
romszinfos = {h:romszinfos.h, zeta:temporary(zeta), theta_s:theta_s, theta_b:theta_b, hc:hc}
ENDIF ELSE romszinfos = {h:-1, zeta:-1, theta_s:-1, theta_b:-1, hc:-1}
ENDIF
;---------------------------------------------------------------------
IF keyword_set(direc) THEN BEGIN
IF jpt EQ 1 THEN res = moyenne(temporary(res), direc) $
ELSE BEGIN
res = grossemoyenne(temporary(res), direc)
IF ( strpos(strlowcase(direc), 't') ge 0 ) THEN BEGIN
vardate = strtrim(jul2date(time[0]), 1)+' - '+strtrim(jul2date(time[jpt-1]), 1)
time = total(time)/float(jpt)
jpt = 1
ENDIF
ENDELSE
ENDIF
;---------------------------------------------------------------------
ncdf_close, cdfid
IF keyword_set(savedbox) THEN restoreboxparam, 'boxparam4rdncdf.dat'
IF keyword_set(nostruct) THEN return, res
IF keyword_set(key_forgetold) THEN BEGIN
return, {arr:temporary(res), grid:vargrid, unit:varunit, experiment:varexp, name:varname}
ENDIF ELSE BEGIN
return, {tab:temporary(res), grille:vargrid, unite:varunit, experience:varexp, nom:varname}
ENDELSE
END