source: trunk/SRC/ToBeReviewed/LECTURE/read_ncdf.pro @ 336

;+
;
; @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 <pro>initncdf</pro> 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 ZINVAR {type=named variable}
; Set this keyword to a named variable in which 1 is returned if a
; vertical dimension is found in the variable. Returns 0 otherwise
;
; @keyword _EXTRA
; Used to pass keywords to <pro>isafile</pro>, <pro>initncdf</pro>,
; <pro>ncdf_gettime</pro> and <pro>domdef</pro>
;
; @returns
; Structure readable by <pro>litchamp</pro> 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 $
                    , ZINVAR = zinvar, _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
;------------------------------------------------------------
; we check if the variable has a vertical dimension and/or 
; a record dimension. This is useful to define boxzoom 
; keyword in domdef
;------------------------------------------------------------
  dummy = where(varinq.dim EQ inq.recdim, tinvar) 
; check the presence of a vertical dimension according to the
; number of dimensions and the presence of a record dimension
  zinvar = (varinq.ndims EQ 3 AND tinvar NE 1) OR varinq.ndims EQ 4
;------------------------------------------------------------
; 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 N_Elements(Boxzoom) Of
        1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]]
        2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]]
        4:bte = [Boxzoom, vert1, vert2]
        5:bte = [Boxzoom[0:3], 0, Boxzoom[4]]
        6:bte = Boxzoom
        ELSE: BEGIN
          ncdf_close, cdfid
          return, report('Wrong Definition of Boxzoom')
        END
      ENDCASE
      IF zinvar EQ 1 THEN domdef, bte, GRIDTYPE = ['T', vargrid], _extra = ex $
      ELSE domdef, bte[0:3], 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(abs(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, _extra = ex) $
    ELSE BEGIN
      res = grossemoyenne(temporary(res), direc, _extra = ex)
      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(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