source: trunk/rseries_ncdf.pro

;+
;
; @file_comments
; reading timeseries writing on many files
;
; @categories
; Reading
;
; @param var {in} {type=string}
; the name of the variable to be read
;
; @param date1 {in} {type=scalar}
; give the first date of the time series.
; Date format is ymmdd, yymmdd, yyymmdd, yyyymmdd, yyyyymmdd, yyyyyymmdd.
;
; @param date2  {in} {type=scalar}
; give the last date of the time series.
; Date format is ymmdd, yymmdd, yyymmdd, yyyymmdd, yyyyymmdd, yyyyyymmdd.
;
; @param expin {in} {type=string} {default=varexp}
; the name of the experiment.
;
; @param freqin {in} {type=string} {default='5d'}
; the averaging frequency.
;
; @keyword CENTURY {type=integer}
; give the true century of the calendar.
; for example if filename is AA5_5d_920101_921231* and the calendar variable
; in the file is 19920101 to 19921231 then CENTURY=19
;
; @keyword NOSTRUCT {default=1}
; Set this keyword to return an array instead of a structure.
;
; @keyword GRIDTYPE
; name of the grid to be read 'grid_T' ++'gridT'
;
; @keyword _EXTRA
; Used to pass read_ncdf keywords
;
; @returns
; a structure with the following Structure Tags:
;  arr:the array output
;  grid: the array grid
;  units: the array units
;  experiment: the name of the experiment
;  name: the name of the variable
;
; @uses
; <pro>common</pro>
;
; @restrictions
; Update the time and jpt common variables.
; Time is the calendar in IDL julian days.
; 1) The file must contain an unlimited dimension, and the
; variable to be read must contain this unlimited dimension.
; 2) The file must contain a one dimension variable which
; dimension is the unlimited dimension. This variable is the
; calendar.
; 3) This variable calendar must have an attribut called "units"
; which must have a format similar to
; "seconds since 0001-01-01 00:00:00"
; "hours since 0001-01-01 00:00:00"
; "days since 1979-01-01 00:00:00"
; "months since 1979-01-01 00:00:00"
; "years since 1979-01-01 00:00:00"
; 4) The calendar must be the Gregorian calendar
; 5) The name of the file must beginning by
; exp_freq_datefirst_datelast_*
; 6) The maximum gap between 2 consecutive files must be one day.
;
; @examples
; IDL> res = rseries_ncdf('sozotaux',920501,930410,'AA5','5d')
;
; for ORCA2 outputs
; IDL> iodir='/usr/work/sur/fvi/OPA/ORCA2/'
; IDL> res= rseries_ncdf('votemper',010101,061231,'ESS','5d')
; IDL> help,res,/structure
;
;** Structure <310ebc0>, 5 tags, length=1800, data length=1800, refs=1:
;   ARR             FLOAT     Array[1, 1, 1, 438]
;   GRID            STRING    'T'
;   UNITS           STRING    'C'
;   EXPERIMENT      STRING    'ESS_5d_060101_061231_grid_T.nc'
;   NAME            STRING    'votemper'
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr) Apr 23 2001
;
; @todo
; adaptation to Drakkar terminologie
; ORCA2 terminology of output files
; <i>experience</i>_5d_<i>yymmdd</i>_<i>yymmdd</i>_<i>grid</i>.nc
; with experience in 'ESS'
; with grid in 'grid_T'
;
; ORCA025 terminology of output files
; /u/rech/cli/rcli002/ORCA025/ORCA025-<i>experience_drakkar<i>-S/<i>year_drakkar</i>/ORCA025-<i>experience_drakkar</i>-y<i>year_drakkar</i>n<i>mois</i>d<i>jour</i>_<i>grid</i>.nc
; with experience_drakkar in 'G70',
; with grid in 'gridT'
;
; @version
; $Id$
;
;-
FUNCTION rseries_ncdf, var, date1, date2, expin, freqin, CENTURY = century, NOSTRUCT = nostruct, GRIDTYPE = gridtype, _EXTRA = ex
@common
; name of the file: exp_freq_datefirst_datelast_*
  date1 = (long(date1))[0]
  date2 = (long(date2))[0]
  if date2 LT date1 then return, report('date2 must be larger than date1')
  if n_elements(expin) EQ 0 then exp = varexp ELSE BEGIN
    IF (strpos(expin, '_'))[0] NE -1 AND n_elements(freqin) EQ 0 THEN BEGIN
      exp = strmid(expin, 0, strpos(expin, '_'))
      freqin = strmid(expin, strpos(expin, '_')+1)
    ENDIF ELSE exp = expin
    varexp = exp
  ENDELSE
  if n_elements(freqin) EQ 0 then freq = '5d' ELSE freq = freqin
;------------------------------------------------
; determination of the filename beginning with exp_freq_*
;------------------------------------------------
  IF keyword_set(gridtype) EQ 0 THEN gridtype = ''
; determination OF datefirst and datelast
; list of the files beginning by iodir+exp+sep+freq
  sep = '_'
  possiblenames = iodir+exp+sep+freq+'*'+gridtype+'.nc'
  possiblenames = findfile(possiblenames)
  if possiblenames[0] EQ '' then return, report('nofilename :'+iodir+exp+sep+freq+'*'+gridtype+'.nc'+' found')
  npos = n_elements(possiblenames)
; list of the datefirst and datelast of the files beginning by
; iodir+exp+sep+freq
  datefirst = lonarr(npos)
  datelast = lonarr(npos)
  sepshift = n_elements(str_sep(exp, sep))-1
  for i = 0, npos-1 do BEGIN
    separate = str_sep(possiblenames[i], sep)
    datefirst[i] = separate[2+sepshift]
    datelast[i] = separate[3+sepshift]
  ENDFOR
; selection of the files for which datefirst le date1 and datelast ge
; date1
  CASE strmid(freqin, 0, 1, /rever) OF
    'm' : div = 100
    'y' : div = 10000
    ELSE : div = 1
  ENDCASE
  goodfile = where(datefirst le date1/div and datelast GE date1/div)
  if goodfile[0] EQ -1 then return, report('filename :'+iodir+exp+sep+freq+'*'+gridtype+'.nc'+' not found with the dates containing '+strtrim(date1, 1))
  datefirst = datefirst[goodfile]
  datelast = datelast[goodfile]
; how many different pairs of dates are present ?
; find the unique pairs of dates
  pairofdate = dcomplex(datefirst, datelast)
  pairofdate = pairofdate[uniq(pairofdate, sort(pairofdate))]
  fileok = ''
  datefirstok = 0L
  datelastok = 0L
; loop on the files for which date1 is between datefirst and datelast
  for p = 0, n_elements(pairofdate)-1 do begin
; name of the file: exp_freq_datefirst_datelast_*
    datefirst = long(pairofdate[p])
    datelast = long(imaginary(pairofdate[p]))
; determination OF the possible filenames
; test if the year format is i1, i2, i3, i4, i5 or i6
    namebegin = iodir+exp+sep+freq+sep+strtrim(datefirst, 1)+sep+strtrim(datelast, 1)+'*'+gridtype+'.nc'
    namebegin = findfile(namebegin)
    possiblenames = ''
    if namebegin[0] NE '' then possiblenames = [possiblenames, namebegin]
    if datefirst LT 1e9 then begin
      for i = 1, 5-(datefirst GE 1e5)-(datefirst GE 1e6) $
        -(datefirst GE 1e7)-(datefirst GE 1e8) do begin
        zero = string(replicate(byte('0'), i))
        namebegin = iodir+exp+sep+freq+sep $
          + zero+strtrim(datefirst, 1) +sep+ zero+strtrim(datelast, 1) +'*'+gridtype+'.nc'
        namebegin = findfile(namebegin)
        if namebegin[0] NE '' then possiblenames = [possiblenames, namebegin]
      endfor
    ENDIF
    if n_elements(possiblenames) eq 1 then return, report('filename :'+iodir+exp+sep+freq+sep+'*'+strtrim(datefirst, 1)+sep+'*'+strtrim(datelast, 1)+'*'+gridtype+'.nc'+' not found')
    possiblenames = possiblenames[1:n_elements(possiblenames)-1]
; determination OF the filenames
    ncdf_control, 0, /noverbose
    i = 0
    repeat BEGIN
      cdfid = ncdf_open(possiblenames[i])
;         print,possiblenames[i],var
      test = ncdf_varid(Cdfid, var)
      ncdf_close, cdfid
      i = i+1
    endrep until test NE -1 OR i EQ n_elements(possiblenames)
    ncdf_control, 0, /verbose
    if test NE -1 then begin
      fileok = [fileok, possiblenames[i-1]]
      datefirstok = [datefirstok, datefirst]
      datelastok = [datelastok, datelast]
    endif
  ENDFOR
  numoffileok = n_elements(fileok)-1
  if numoffileok EQ 0 then return, report('the variable '+var+' was not fond in the file: '+possiblenames)
  fileok = fileok[1:numoffileok]
  datefirstok = datefirstok[1: numoffileok]
  datelastok = datelastok[1: numoffileok]
  if numoffileok NE 1 then BEGIN
; more than one file containts the variable var and has the date1
; between datefirst and datelast...
; Is there any file FOR which date2 LE datelast?
    goodfile = where(date2 LE datelastok)
    if goodfile[0] NE -1 then begin
      fileok = fileok[goodfile]
      datefirstok = datefirstok[goodfile]
      datelastok = datelastok[goodfile]
; we choose the file which has the smallest datelast
      bestfile = sort(datelastok)
      filename = fileok[bestfile[0]]
      datefirst = datefirstok[bestfile[0]]
      datelast = datelastok[bestfile[0]]
    ENDIF ELSE BEGIN
; we choose the file which has the bigest datelast
      bestfile = sort(datelastok)
      filename = fileok[bestfile[numoffileok-1]]
      datefirst = datefirstok[bestfile[numoffileok-1]]
      datelast = datelastok[bestfile[numoffileok-1]]
    ENDELSE
  ENDIF ELSE BEGIN
    filename = fileok[0]
    datefirst = datefirstok[0]
    datelast = datelastok[0]
  ENDELSE
; now we have find the good file with the good dates and the godd
; variable inside.
; we get the number of dimension of the variable
  cdfid = ncdf_open(filename)
  varcontient = ncdf_varinq(cdfid, var)
  ncdf_close, cdfid
;
;------------------------------------------------
; reading of the variables
;------------------------------------------------
;
  if NOT keyword_set(century) then century = 0l
  IF chkstru(ex, 'filename') THEN ex.filename = filename
  CASE strmid(freqin, 0, 1, /rever) OF
    'm':datelast = 100*datelast + daysinmonth(datelast MOD 100, datelast / 100)
    'y':div = 10000*datelast+360+(5+leapyr(datelast+century))*(key_caltype EQ 'greg')
    ELSE : div = 1
  ENDCASE
  if date2 GT datelast THEN BEGIN
; if we need to read more than one file,
; first we read the first file
;++    print,var,date1+century*1000000L, datelast[0]+century*1000000L
;++READ, B, PROMPT='Enter Name: '
    res1 = read_ncdf(var, date1+century*1000000L, datelast[0]+century*1000000L $
                     , filename = filename, /nostruct, _extra = ex)
    time1 = time                ; store the first part of the calendar
    jpt1 = jpt                  ; store the number time steps already read

; and after we call again rseries_ncdf
; newdatefirst is defined as datelast+1day
    newdatefirst = jul2date(date2jul(datelast)+1)
    res2 = rseries_ncdf(var, newdatefirst, date2, exp, freq, /NOSTRUCT, _extra = ex)
;
;
    CASE 1 OF
      res1[0] EQ -1 AND res2[0] EQ -1: return, -1
      res1[0] EQ -1:res = temporary(res2)
      res2[0] EQ -1:BEGIN
        time = time1
        jpt = jpt1
        res = temporary(res1)
      END
      ELSE : BEGIN
 ; we glue the result of the first read and the result of the new call
; to rseries_ncdf
        time = [time1, time]
        jpt = jpt1+jpt
        case varcontient.ndims OF
          1:res = [temporary(res1), temporary(res2)]
          2:res = [ [temporary(res1)], [temporary(res2)] ]
          3:res = [ [ [temporary(res1)] ], [ [temporary(res2)] ] ]
          4:BEGIN
            res = [(temporary(res1))[*], (temporary(res2))[*]]
            grille, mask, glam, gphi, gdep, nx, ny, nz
            res = reform(res, nx, ny, nz, jpt, /over)
          END
        ENDCASE
      ENDELSE
    ENDCASE
  endif ELSE BEGIN
    res = read_ncdf(var, date1+century*1000000L, date2+century*1000000L $
                    , filename = filename, /nostruct, _extra = ex)
  ENDELSE
;------------------------------------------------
  if keyword_set(nostruct) then return, res $
  ELSE return, {arr:res, grid:vargrid, units:varunit, experiment:varexp, name:varname}

end