source: trunk/SRC/ReadWrite/ncdf_gettime.pro @ 325

;+
;
; @file_comments
; get the time axis from a netcdf_file and transforms it in
; Julian days of IDL.
;
; @categories
; Read NetCDF file
;
; @param filename {in}{required}{type=scalar string}
; the name of the ncdf_file
;
; @param cdfid {in}{required}{type=scalar}
; the ID of the ncdf_file, which is already open
;
; @keyword TIMEVAR {type=string}
; It define the name of the variable that
; contains the time axis. This keyword can be useful if there
; is no unlimited dimension or if the time axis selected by default
; (the first 1D array with unlimited dimension) is not the good one.
;
; @keyword CALLER {required}{type=scalar}
; Used to specify the error messages. Give the name of the calling
; procedure. It can be only 'read_ncdf' or 'scanfile'.
;
; @keyword ERR
; Set this keyword to a named variable in which the value of the error
; message will be returned
;
; @keyword _EXTRA
; _EXTRA to be able to call ncdf_getmask with _extra keyword
;
; @returns
; a double 1D array of IDL Julian days
; In case of error return -1 if the time dimension was not found
;               or return -jpt if it as been found that the time dimension size is jpt
;
; @restrictions
; the calendar variable must have the units attribute
; following the syntax bellow:
;
; 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" ;
;
; @history
; August 2007: Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;-
FUNCTION ncdf_gettime, filename, cdfid, TIMEVAR = timevar, CALLER = caller, ERR = err, _EXTRA = ex
;
  compile_opt idl2, strictarrsubs
;
  @cm_4cal                      ;   needed for key_caltype
;
  inq = ncdf_inquire(cdfid)
;----------------------------------------------------
;`find the variable containing the time axis
;----------------------------------------------------
; we get its name through the keyword timevar
  IF keyword_set(timevar) THEN BEGIN
    timeid = ncdf_varid(cdfid, timevar)
    IF timeid EQ -1 THEN BEGIN   ;   the variable is not found
      CASE caller OF
        'read_ncdf':err = 'No variable ' + timevar + 'found in '+filename+'. Use the TIMESTEP keyword'
        'scanfile':err = 'No variable ' + timevar + 'found in '+filename+'. We create a fake calendar'
      ENDCASE
      return, -1
    ENDIF
    timeinq = ncdf_varinq(cdfid, timeid)
    inq.recdim = timeinq.dim[0]
    ncdf_diminq, cdfid, inq.recdim, timedimname, jpt
  ENDIF ELSE BEGIN
; we try to find the time axis automatically
; we look for the infinite dimension
    IF inq.recdim EQ -1 THEN BEGIN
      CASE caller OF
        'read_ncdf':err = 'the file '+filename+' as no infinite dimension. !C Use TIMESTEP or TIMEVAR keyword'
        'scanfile':err = 'the file '+filename+' as no infinite dimension. We create a fake calendar'
      ENDCASE
      return, -1
    ENDIF
    ncdf_diminq, cdfid, inq.recdim, timedimname, jpt
; we look for the variable containing the time axis
; we look for the first variable having for only dimension inq.recdim
    timeid = 0
    REPEAT BEGIN                ; As long as we have not find a variable having only one dimension: the infinite one
      timeinq = ncdf_varinq(cdfid, timeid) ; that the variable contain.
      timeid = timeid+1
    ENDREP UNTIL (n_elements(timeinq.dim) EQ 1 AND timeinq.dim[0] EQ inq.recdim) $
       OR timeid EQ inq.nvars+1
    IF timeid EQ inq.nvars+1 THEN BEGIN
      CASE caller OF
        'read_ncdf':err = 'the file '+filename+' as no time axis variable. !C Use the TIMESTEP keyword'
        'scanfile':err = 'the file '+fullname+' has no time axis.!C we create a fake calendar ...'
      ENDCASE
      return, -jpt
    ENDIF
    timeid = timeid-1
  ENDELSE
;----------------------------------------------------
; look for attribute units and calendar to know how to compte the calendar
;----------------------------------------------------
; no attribute for time variable...
  IF timeinq.natts EQ 0 then begin
    CASE caller OF
      'read_ncdf':err = 'the variable '+timeinq.name+' has no attribute.!C Use the TIMESTEP keyword'
      'scanfile':err = 'the variable '+timeinq.name+' has no attribute.!C we create a fake calendar ...'
    ENDCASE
    return, -jpt
  ENDIF
; get attributes names
  attnames = strarr(timeinq.natts)
  for attiq = 0, timeinq.natts-1 do attnames[attiq] = strlowcase(ncdf_attname(cdfid, timeid, attiq))
; do we find units attribute?
  IF (where(attnames EQ 'units'))[0] EQ -1 then BEGIN
    CASE caller OF
      'read_ncdf':err = 'Attribute ''units'' not found for the variable '+timeinq.name+' !C Use the TIMESTEP keyword'
      'scanfile':err = 'Attribute ''units'' not found for the variable '+timeinq.name+'!C we create a fake calendar ...'
    ENDCASE
    return, -jpt
  ENDIF
; Is attribute "calendar" existing?
; 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 = strlowcase(strtrim(value, 2))
; unpade cm_4cal variable key_caltype (used by julday and caldat)
    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
;----------------------------------------------------
; decode units attribute
;----------------------------------------------------
  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" attribute to find the time origin
  value = strtrim(strcompress(string(value)), 2)
  words = str_sep(value, ' ')
  units = words[0]
  units = strlowcase(units)
  IF strpos(units, 's', strlen(units)-1) NE -1 THEN units = strmid(units, 0, strlen(units)-1)
  IF strpos(units, 'julian_') NE -1 THEN units = strmid(units, 7)
  IF units NE 'second' AND units NE 'hour' AND units NE 'day' $
     AND units NE 'month' AND units NE 'year' THEN BEGIN
    CASE caller OF
      'read_ncdf':err = 'time units does not start with seconds/hours/days/months/years !C Use the TIMESTEP keyword'
      'scanfile':err = 'time units does not start with seconds/hours/days/months/years !C we create a fake calendar ...'
    ENDCASE
    return, -jpt
  ENDIF
  IF stregex(value, '[^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*', /boolean) EQ 0 THEN BEGIN
    CASE caller OF
      'read_ncdf':err = 'attribute units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*!C Use the TIMESTEP keyword'
      'scanfile':err = 'attribute units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*!C we create a fake calendar ...'
    ENDCASE
    return, -jpt
  ENDIF
  start = str_sep(words[2], '-')
;----------------------------------------------------
; compute time axis
;----------------------------------------------------
  ncdf_varget, cdfid, timeid, time
  time = double(time)
  case units of
    'second':time = julday(start[1], start[2], start[0], 0, 0, 0)+time/86400.d
    'hour':time = julday(start[1], start[2], start[0], 0, 0, 0)+time/24.d
    'day':time = julday(start[1], start[2], start[0], 0, 0, 0)+time
    'month':BEGIN
      if total(fix(time) NE time) NE 0 then $ ; we switch to days with 30d/m
         time = julday(start[1], start[2], start[0])+round(time*30) $
      ELSE for t = 0, n_elements(time)-1 DO $
         time[t] = julday(start[1]+time[t], start[2], start[0])
    END
    'year':BEGIN
      if total(fix(time) NE time) NE 0 then $ ; we switch to days with 365d/y
         time = julday(start[1], start[2], start[0])+round(time*365) $
      ELSE for t = 0, n_elements(time)-1 do $
         time[t] = julday(start[1], start[2], start[0]+time[t])
    END
  ENDCASE
  time = double(time)
;
  return, time
END