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

;+
; @file_comments
; Reading function for the file net_cdf.
; This program is less universal than ncdf_lec (it appeal to declarated 
; 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 inour reading widgets.  
;
; @categories
; reading
;
; 
; @param NAME {in}{required}
; It is a string defining 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 telporal dimension.
;
; @param ENDING  {in}{required}
; Relative with the time axis.
; See BEGINNING.
; 
; @param COMPATIBILITY {in}{required}
; Useless
;  
; @keyword BOXZOOM 
; Contain the boxzoom on which we have to do the reading 
; 
; @keyword FILENAME
; It is a string containing the file's name.
; 
; @keyword INIT
; To call automatically initncdf, filename and thus
; redefine all the grid parameters
; 
; @keyword 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.
; 
; @keyword  IODIRECTORY
; It isa string giving the name of iodirectory (see isafile.pro 
; for all possibilities). default value is common variable iodir.
; 
; @keyword TIMESTEP
; We activate to specify that BEGINNING and ENDING refer to indexes
; of the time axis and not to dates
;
; @keyword TOUT
; 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
; We activate it if we do not want that read_ncdf send back a structure 
; but only the array refering to the field.
; 
; @keyword TIMEVAR
; It is a string which define the name of the variable that
; contains the time axis. This keyword can be usefull if there
; is no unlimited dimension or if the time axis selected by defaut
; (the first 1D array with unlimited dimension) is not the good one.
;
; @keyword _EXTRA
; Used to pass your keywords 
;
; @returns
; Structure readable by litchamp.pro or an array if NOSRUCT 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, TIMEVAR = timevar $
                    , TOUT = tout, NOSTRUCT = nostruct, CONT_NOFILL = CONT_NOFILL, INIT = init $
                    , GRID = grid, FBASE2TBASE = fbase2tbase, _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)
  contient = 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
  varcontient = ncdf_varinq(cdfid, name)
;------------------------------------------------------------
; 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 = 0
    timestep = 1
  endif
  if keyword_set(timestep) then begin
    firsttps = beginning[0]
    if n_elements(ending) NE 0 then lasttps = ending[0] ELSE lasttps = firsttps
    jpt = lasttps-firsttps+1
    time = julday(1, 1, 1) + lindgen(jpt)
  ENDIF ELSE BEGIN
    if keyword_set(parent) then 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(time EQ date1) & firsttps = firsttps[0]
      lasttps = where(time EQ date2) & lasttps = lasttps[0]
    ENDIF ELSE BEGIN
      IF keyword_set(timevar) THEN BEGIN 
        timeid = ncdf_varid(cdfid, timevar)
        IF timeid EQ -1 THEN BEGIN
          ncdf_close, cdfid
          return, report('the file '+filename+' as no variable ' + timevar $
                         + '. !C Use the TIMESTEP keyword')
        endif
        timecontient = ncdf_varinq(cdfid, timeid)
        contient.recdim = timecontient.dim[0]
      ENDIF ELSE BEGIN 
; we find the infinite dimension
        timedim = contient.recdim
        if timedim EQ -1 then BEGIN
          ncdf_close, cdfid
          return, report('the file '+filename+' as no infinite dimension. !C Use TIMESTEP or TIMEVAR keyword')
        endif
; we find the FIRST time axis      
        timeid = 0
        repeat BEGIN       ; As long as we have not find a variable having only one dimension: the infinite one
          timecontient = ncdf_varinq(cdfid, timeid) ; that the variable contain.
          timeid = timeid+1
        endrep until (n_elements(timecontient.dim) EQ 1 $
                      AND timecontient.dim[0] EQ contient.recdim) $
          OR timeid EQ contient.nvars+1
;
        if timeid EQ contient.nvars+1 then BEGIN
          ncdf_close, cdfid
          return, report('the file '+filename+' as no time axis variable. !C Use the TIMESTEP keyword')
        endif
        timeid = timeid-1
      ENDELSE 
; we must found the time origin of the julian calendar used in the
; time axis. 
; does the attribut units an dcalendar exist for the variable time axis?
      if timecontient.natts EQ 0 then BEGIN
        ncdf_close, cdfid
        return, report('the variable '+timecontient.name+' has no attribut.!C Use the TIMESTEP keyword or add the attribut ''units'' to the variable')
      endif
      attnames = strarr(timecontient.natts)
      for attiq = 0, timecontient.natts-1 do attnames[attiq] = ncdf_attname(cdfid, timeid, attiq)
      if (where(attnames EQ 'units'))[0] EQ -1 then BEGIN
        ncdf_close, cdfid
        return, report('Attribut ''units'' not found for the variable '+timecontient.name+' !C Use the TIMESTEP keyword')
      ENDIF
;
; 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, timeid, '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
;
; now we take acre of units attribut
      ncdf_attget, cdfid, timeid, '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" ;
;
; we decript the "units" attribut to find the time origin
      value = strtrim(strcompress(string(value)), 2)
      mots = str_sep(value, ' ')
      unite = mots[0]
      depart = str_sep(mots[2], '-')
      ncdf_varget, cdfid, timeid, time
      time = double(time)
      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)
      case unite of
        'second':time = julday(depart[1], depart[2], depart[0])+time/86400.d
        'hour':time = julday(depart[1], depart[2], depart[0])+time/24.d
        'day':time = julday(depart[1], depart[2], depart[0])+time
        'month':BEGIN 
          if total(fix(time) NE time) NE 0 then $ ; we switch to days with 30d/m
            time = julday(depart[1], depart[2], depart[0])+round(time*30) $
            ELSE for t = 0, n_elements(time)-1 DO $
            time[t] = julday(depart[1]+time[t], depart[2], depart[0])
        END
        'year':BEGIN
          if total(fix(time) NE time) NE 0 then $ ; we switch to days with 365d/y
            time = julday(depart[1], depart[2], depart[0])+round(time*365) $
            ELSE for t = 0, n_elements(time)-1 do $
            time[t] = julday(depart[1], depart[2], depart[0]+time[t])
        END
        ELSE:BEGIN
          ncdf_close, cdfid
          return, report('The "units" attribu of the time axis must be something like: !C "seconds since 0001-01-01 ..." !C "days since 1979-01-01 ..." !C "months since 1979-01-01 ..." !C "years since 1979-01-01 ..."')
        end
      ENDCASE
      date1 = date2jul(beginning[0])
      if n_elements(ending) NE 0 then date2 = date2jul(ending[0]) ELSE date2 = date1
      time = double(time)
      firsttps = where(time GE date1) & 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 = where(time LE date2)
      if lasttps[0] EQ -1 THEN BEGIN
        ncdf_close, cdfid
        return, report('the time axis as no date before date 2: '+strtrim(jul2date(date2), 1))
      endif
      lasttps = lasttps[n_elements(lasttps)-1]
      if lasttps LT firsttps then BEGIN
        ncdf_close, cdfid
        return, report('the time axis as no dates between date1 and  date 2: '+strtrim(jul2date(date1), 1)+' '+strtrim(jul2date(date2), 1))
      endif
    ENDELSE
    time = time[firsttps:lasttps]
    jpt = lasttps-firsttps+1
  ENDELSE
;------------------------------------------------------------
; 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
      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
    ENDIF
  ENDELSE
;---------------------------------------------------------------
; call the init function ?
;---------------------------------------------------------------

;---------------------------------------------------------------
; 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 memoty!
  ENDELSE
;---------------------------------------------------------------------
; We initializate 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)
  key_shift = long(testvar(var = key_shift))
;
  IF n_elements(key_yreverse) EQ 0 THEN key_yreverse = 0
  IF keyword_set(key_yreverse) THEN BEGIN
    tmp = jpj-1-firsty
    firsty = jpj-1-lasty
    lasty = tmp
  ENDIF
;
  IF keyword_set(fbase2tbase) THEN BEGIN
    case strupcase(vargrid) of
      'U':BEGIN
        IF NOT keyword_set(key_periodic) THEN BEGIN
          firstx = firstx+1
          lastx = lastx+1
        ENDIF
      END
      'V':BEGIN
        firsty = firsty+1
        lasty = lasty+1
      END
      'F':BEGIN
        firsty = firsty+1
        lasty = lasty+1
        IF NOT keyword_set(key_periodic) THEN BEGIN
          firstx = firstx+1
          lastx = lastx+1
        ENDIF
      END
      ELSE:
    endcase
  ENDIF
;
  IF keyword_set(fbase2tbase) AND keyword_set(key_periodic) $
    AND (strupcase(vargrid) EQ 'U' OR strupcase(vargrid) EQ 'F') THEN key_shift =  key_shift-1
;
;---------------------------------------------------------------------
;---------------------------------------------------------------------
@read_ncdf_varget
;---------------------------------------------------------------------
;---------------------------------------------------------------------
;
  IF keyword_set(fbase2tbase) AND keyword_set(key_periodic) $
    AND (strupcase(vargrid) EQ 'U' OR strupcase(vargrid) EQ 'F') THEN key_shift =  key_shift+1
;---------------------------------------------------------------------
; We define global variable joined with the variable.
;---------------------------------------------------------------------
; varname
  varname = name
; varunit
  if varcontient.natts NE 0 then begin
    attnames = strarr(varcontient.natts)
    for attiq = 0, varcontient.natts-1 do attnames[attiq] = ncdf_attname(cdfid, name, attiq)
    lowattnames = strlowcase(attnames)
;
    found = (where(lowattnames EQ 'units'))[0]
    IF found NE -1 then ncdf_attget, cdfid, name, attnames[found], value ELSE value = ''
    varunit = strtrim(string(value), 2)
;
    found = (where(lowattnames EQ 'add_offset'))[0]
    if found NE -1 then ncdf_attget, cdfid, name, attnames[found], add_offset ELSE add_offset = 0.
;
    found = (where(lowattnames EQ 'scale_factor'))[0]
    if found NE -1 then ncdf_attget, cdfid, name, attnames[found], scale_factor ELSE scale_factor = 1.
;
    missing_value = 'no'
    found = (where(lowattnames EQ '_fillvalue'))[0]
    if found NE -1 then ncdf_attget, cdfid, name, attnames[found], missing_value
    found = (where(lowattnames EQ 'missing_value'))[0]
    if found NE -1 then ncdf_attget, cdfid, name, attnames[found], missing_value
;
  ENDIF ELSE BEGIN 
    varunit = ''
    add_offset = 0.
    scale_factor = 1.
    missing_value = 'no'
  ENDELSE
; 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 reverse
  if keyword_set(key_yreverse) then res = reverse(temporary(res),  2)
  if keyword_set(key_zreverse) AND (size(res))[0] EQ 3 AND jpt EQ 1 then res = reverse(temporary(res), 3)
  if keyword_set(key_zreverse) AND (size(res))[0] EQ 4 THEN res = reverse(temporary(res), 3)

; We apply the value valmask on land points.
  if NOT keyword_set(cont_nofill) then begin
    valmask = 1e20
    case 1 of
      varcontient.ndims eq 2:BEGIN ;xy array
        mask = mask[*, *, 0]
        earth = where(mask EQ 0)
      END
      varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] EQ -1:BEGIN ;xyz array
        earth = where(mask EQ 0)
      END
      varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] NE -1:BEGIN ;xyt array
        mask = mask[*, *, 0]
        earth = where(mask EQ 0)
        if earth[0] NE -1 then BEGIN
          earth = earth#replicate(1, jpt)+replicate(nx*ny, n_elements(earth))#lindgen(jpt)
        END
      END
      varcontient.ndims eq 4:BEGIN ;xyzt array
        earth = where(mask EQ 0)
        if earth[0] NE -1 then BEGIN
          earth = earth#replicate(1, jpt)+replicate(nx*ny*nz, n_elements(earth))#lindgen(jpt)
        END
      END
    endcase
  ENDIF ELSE earth = -1
; we look for  missing_value
  IF size(missing_value, /type) NE 7 then BEGIN
    IF size(missing_value, /type) EQ 1 THEN BEGIN 
      IF isnumber(string(missing_value), tmp) EQ 1 THEN missing_value = tmp
    ENDIF 
;    if missing_value NE valmask then begin
    if abs(missing_value) LT 1e6 then missing = where(res EQ missing_value) $
    ELSE missing = where(abs(res) gt abs(missing_value)/10.)
;    ENDIF ELSE missing = -1
  ENDIF ELSE missing = -1
; we apply add_offset, scale_factor and missing_value 
  if scale_factor NE 1 then res = temporary(res)*scale_factor
  if add_offset NE 0 then res = temporary(res)+add_offset
  if missing[0] NE -1 then res[temporary(missing)] = !values.f_nan
  if earth[0] NE -1 then res[temporary(earth)] = 1e20
;---------------------------------------------------------------------
  ncdf_close, cdfid
;---------------------------------------------------------------------
  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4rdncdf.dat'
  if keyword_set(nostruct) then return, res $
  ELSE BEGIN 
    IF keyword_set(key_forgetold) THEN BEGIN
      return, {arr:res, grid:vargrid, unit:varunit, experiment:varexp, name:varname} 
    ENDIF ELSE BEGIN 
      return, {tab:res, grille:vargrid, unite:varunit, experience:varexp, nom:varname}
    ENDELSE 
  ENDELSE 
END