source: trunk/SRC/ToBeReviewed/INIT/initncdf.pro @ 236

;+
;
; @file_comments
; Initfile for Netcdf file. define all the grid parameters through
; an appropriate call to computegid
;
; @categories
; Grid
;
; @param NCFILEIN {in}{required}{type=scalar string}
; A string giving the name of the NetCdf file
;
; @keyword INVMASK {default=0}{type=scalar: 0 or 1}
; Inverse the land/sea mask (that should have 0/1 values for land/sea): mask = 1-mask
;
; @keyword MASKNAME {type=string}
; A string giving the name of the variable in the file
; that contains the land/sea mask
;
; @keyword MISSING_VALUE {type=scalar}
; To define (or redefine if the attribute is
; already existing) the missing values used with USEASMASK
; keyword
;
; @keyword START1 {default=0}{type=scalar: 0 or 1}
; Index the axis from 1 instead of 0 when using
; /xyindex and/or /zindex
;
; @keyword USEASMASK {type=scalar string}
; A string giving the name of the variable in the file
; that will be used to build the land/sea mask. In this case the
; mask is based on the first record (if record dimension
; exists). The mask is build according to :
;    1 the keyword missing_value if existing
;    2 the attribute 'missing_value' if existing
;    3 NaN values if existing
;
; @keyword ZAXISNAME {default='z', 'level', 'lev', 'depth...'}{type=scalar string}
; A string giving the name of the variable in the file
; that contains the [xyz]axis.
;
; @keyword XYINDEX {default=0}{type=scalar: 0 or 1}
; To define the x/y axis with index instead of using
; the values contained in X/YAXISNAME.
; x/yaxis = keyword_set(start1) + findgen(jpi/jpj)
; this forces key_onearth = 0
;
; @keyword ZINDEX {default=0}{type=scalar: 0 or 1}
; To define the z axis with index instead of using
; the values contained in ZAXISNAME.
; zaxis = keyword_set(start1) + findgen(jpk)
;
; @keyword _EXTRA
; Used to pass keywords to <pro>computegrid</pro> and 
; <pro>ncdf_getaxis</pro>
;
; @uses
; common.pro
;
; @restrictions
; Change the grid parameters (see <pro>computegrid</pro>)
;
; @restrictions
; the file must contain an x and an y axis. (1 ou 2 dimensional array)
;
; @examples
; IDL> initncdf,'toto.nc',glam=[-180,180]
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                      8 May 2002
;
; @version
; $Id$
;
;-
;
PRO initncdf, ncfilein $
              , ZAXISNAME = zaxisname, MASKNAME = maskname $
              , INVMASK = invmask, USEASMASK = useasmask $
              , MISSING_VALUE = missing_value, START1 = start1 $
              , XYINDEX = xyindex, ZINDEX = zindex $
              , _EXTRA = ex
;
  compile_opt idl2, strictarrsubs
;
@common
;----------------------------------------------------------
; check the name of the file
  ncfile = isafile(FILENAME = ncfilein, IODIRECTORY = iodir, _extra = ex)
  if size(ncfile, /type) NE 7 then BEGIN
    ras = report( 'initncdf cancelled')
    return
  endif
; if the file is stored on tape
  if !version.os_family EQ 'unix' then spawn, 'file '+ncfile+' > /dev/null'
;----------------------------------------------------------
; open the file
  cdfid = ncdf_open(ncfile)
; what is inside 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
;----------------------------------------------------------
; name of the variables
  namevar = strarr(inside.nvars)
  for varid = 0, inside.nvars-1 do begin
    invar = ncdf_varinq(cdfid, varid)
    namevar[varid] = strlowcase(invar.name)
  ENDFOR
;----------------------------------------------------------
; find the x/yaxis
 ncdf_getaxis, cdfid, dimidx, dimidy, xaxis, yaxis $
               , START1 = start1, XYINDEX = xyindex, ROMSGRID = romsgrid, _extra = ex
;----------------------------------------------------------
; find the zaxis
  IF keyword_set(romsgrid) THEN BEGIN
    FOR i = 0, inside.ndims-1 DO BEGIN
      ncdf_diminq, cdfid, i, name, size
      CASE strlowcase(name) OF
        's_rho':zaxis = reverse(indgen(size))
        's_u':zaxis = reverse(indgen(size))
        's_v':zaxis = reverse(indgen(size))
        's_psi':zaxis = reverse(indgen(size))
        's_w':zaxis = reverse(indgen(size-1))
        ELSE:
      ENDCASE
    ENDFOR
    IF (where(namevar EQ 'h'))[0] NE -1 THEN BEGIN
      ncdf_varget, cdfid, 'h', romsh
    ENDIF ELSE romsh = -1
  ENDIF ELSE BEGIN
    if keyword_set(zaxisname) then zaxisname = strlowcase(zaxisname) ELSE zaxisname = 'z'
    zvarid = (where(namevar EQ 'nav_lev' or namevar EQ zaxisname OR namevar EQ 'level' OR namevar EQ 'lev' OR strmid(namevar, 0, 5) EQ 'depth'))[0]
    if zvarid EQ -1 AND inside.ndims GT 3 then begin
      ras = report( 'initncdf: the zaxis was not found..., check the the use of ZAXISNAME keyword if you whant to find one...')
;     stop
    endif
; read the zaxis
    if zvarid NE -1 THEN ncdf_varget, cdfid, zvarid, zaxis
  ENDELSE
  IF keyword_set(zindex) AND keyword_set(zaxis) THEN $
     zaxis = keyword_set(start1) + findgen(n_elements(zaxis))
;----------------------------------------------------------
; mask
  IF NOT (keyword_set(maskname) OR keyword_set(useasmask)) AND keyword_set(romsgrid) THEN maskname = 'mask_rho'
  CASE 1 OF
    keyword_set(maskname):BEGIN
      mskid = (where(namevar EQ strlowcase(maskname)))[0]
      if mskid NE -1 THEN BEGIN
        mskinq = ncdf_varinq(cdfid, mskid)
; is the mask variable containing the record dimension?
        withrcd = (where(mskinq.dim EQ inside.recdim))[0]
        IF withrcd NE -1 THEN BEGIN
; in order to read only the first record
; we need to get the size of each dimension
          count = replicate(1L, mskinq.ndims)
          FOR d = 0, mskinq.ndims -1 DO BEGIN
            IF d NE withrcd THEN BEGIN
              ncdf_diminq, cdfid, mskinq.dim[d], name, size
              count[d] = size
            ENDIF
          ENDFOR
; read the variable for the first record
          ncdf_varget, cdfid, mskid, tmask, count = count
        ENDIF ELSE ncdf_varget, cdfid, mskid, tmask
; check if we need to applay add_offset and scale factor
        FOR a = 0, mskinq.natts-1 DO BEGIN
          attname = ncdf_attname(cdfid, mskid, a)
          CASE strlowcase(attname) OF
            'add_offset':ncdf_attget, cdfid, mskid, attname, add_offset
            'scale_factor':ncdf_attget, cdfid, mskid, attname, scale_factor
            ELSE:
          ENDCASE
        ENDFOR
        IF n_elements(scale_factor) NE 0 THEN tmask = tmask*scale_factor
        IF n_elements(add_offset) NE 0 THEN tmask = tmask+add_offset
        if keyword_set(invmask) then tmask = 1-tmask
        tmask = byte(round(tmask))
      ENDIF ELSE tmask = -1
    END
;..................
    keyword_set(useasmask):BEGIN
      mskid = (where(namevar EQ strlowcase(useasmask)))[0]
      if mskid NE -1 THEN BEGIN
        mskinq = ncdf_varinq(cdfid, mskid)
; is the mask variable containing the record dimension?
        withrcd = (where(mskinq.dim EQ inside.recdim))[0]
        IF withrcd NE -1 THEN BEGIN
; in order to read only the first record
; we need to get the size of each dimension
          count = replicate(1L, mskinq.ndims)
          FOR d = 0, mskinq.ndims -1 DO BEGIN
            IF d NE withrcd THEN BEGIN
              ncdf_diminq, cdfid, mskinq.dim[d], name, size
              count[d] = size
            ENDIF
          ENDFOR
; read the variable for the first record
          ncdf_varget, cdfid, mskid, tmask, count = count
        ENDIF ELSE ncdf_varget, cdfid, mskid, tmask
; check if we need to applay add_offset and scale factor
        FOR a = 0, mskinq.natts-1 DO BEGIN
          attname = ncdf_attname(cdfid, mskid, a)
          CASE strlowcase(attname) OF
            'add_offset':ncdf_attget, cdfid, mskid, attname, add_offset
            'scale_factor':ncdf_attget, cdfid, mskid, attname, scale_factor
            'missing_value':IF n_elements(missing_value) EQ 0 THEN $
               ncdf_attget, cdfid, mskid, attname, missing_value
            ELSE:
          ENDCASE
        ENDFOR
        IF n_elements(scale_factor) NE 0 THEN tmask = tmask*scale_factor
        IF n_elements(add_offset) NE 0 THEN tmask = tmask+add_offset
        IF n_elements(missing_value) NE 0 THEN BEGIN
; we have to take care of the float accuracy
          CASE 1 OF
            missing_value GE 1.e6:tmask = tmask LT (missing_value-10)
            missing_value LE -1.e6:tmask = tmask GT (missing_value-10)
            abs(missing_value) LE 1.e-6:tmask = abs(tmask) GT 1.e-6
            ELSE:tmask = tmask NE missing_value
          ENDCASE
          if keyword_set(invmask) then tmask = 1-tmask
        ENDIF ELSE BEGIN
          tmask = finite(tmask)
          IF min(tmask) EQ 1 THEN BEGIN
            ras = report( 'missing or nan values not found...')
            tmask = -1
          ENDIF
        ENDELSE
      ENDIF ELSE tmask = -1
    END
;..................
    ELSE:tmask  = -1
  ENDCASE
;
  ncdf_close, cdfid
;
; compute the grid
  if NOT keyword_set(zaxis) then BEGIN
    computegrid, xaxis = xaxis, yaxis = yaxis $
                 , mask = tmask, onearth = 1b - keyword_set(xyindex), ROMSH = romsh, _EXTRA = ex
  ENDIF ELSE BEGIN
    computegrid, xaxis = xaxis, yaxis = yaxis, zaxis = zaxis $
                 , mask = tmask, onearth = 1b - keyword_set(xyindex), ROMSH = romsh, _EXTRA = ex
  ENDELSE
  IF n_elements(time) EQ 0 THEN time = 0
  jpt = n_elements(time)
;----------------------------------------------------------

  return
end