source: trunk/rseries_ncdf.pro @ 11

Last change on this file since 11 was 11, checked in by pinsard, 17 years ago

introducing rseries_ncdf.pro first step for ORCA025 outputs reading
File size: 10.9 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME: rseries_ncdf
6;
7; PURPOSE: reading timeseries writing on many files.
8;
9; CATEGORY: upper read_ncdf
10;
11; CALLING SEQUENCE: res = rseries_ncdf(var, date1, date2[, exp[, freq]])
12;
13; INPUTS:
14;        var: a scalar string; the name of the variable to be read
15;
16;        date1: a saclar giving the first date of the time
17;        series. Date format is ymmdd, yymmdd, yyymmdd, yyyymmdd,
18;        yyyyymmdd, yyyyyymmdd.
19;
20;        date2: same as date2
21;
22;        exp: a scalar string; the name of the experiment. If exp is
23;        undefined, exp = varexp.
24;
25;
26;        freq: a scalar string; the averaging frequency. if freq is
27;        undefined, freq = '5d'
28;
29; KEYWORD PARAMETERS:
30;
31;        CENTURY: and interger giving the true century of the
32;        calendar. for example if filename is
33;        AA5_5d_920101_921231* and the calendar variable in the file
34;        is 19920101 to 19921231 then CENTURY=19
35;
36;        /NOSTRUCT: Set this keyword to return an array instead of a
37;        structure.
38;
39;       NDAYSPM: developpe par eric, ca veut dire: nombre de jours par mois!
40;                par defaut c''est 30, sinon le specifier en donnant
41;                une valeur a ndayspm
42;                pour passer a un calendrier avec un nombre de jours constant par
43;                mois. Utilise en particulier ds julday et caldat
44;
45; OUTPUTS:
46;
47;        a structure: with the folowing Structure Tags:
48;                     arr:the array output
49;                     grid: the array grid
50;                     units: the array units
51;                     experiment: the name of the experiment
52;                     name: the name of the variable
53;
54; COMMON BLOCKS:common.pro
55;
56; SIDE EFFECTS:
57;
58;        Update the time and jpt common variables. Time is the
59;        calendar in IDL julian days.
60;
61;
62; RESTRICTIONS:
63;
64;        1) The file must contain an unlimited dimension, and the
65;        variable to be read must contain this unlimited dimension.
66;        2) The file must contain a one dimension variable which
67;        dimension is the unlimited dimension. This variable is the
68;        calendar.
69;        3) This variable calendar must have an attribut called "units"
70;        which must have a format similar to
71;                    "seconds since 0001-01-01 00:00:00"
72;                    "hours since 0001-01-01 00:00:00"
73;                    "days since 1979-01-01 00:00:00"
74;                    "months since 1979-01-01 00:00:00"
75;                    "years since 1979-01-01 00:00:00"
76;        4) The calendar must be the Gregorian calendar
77;        5) The name of the file must begining by
78;           exp_freq_datefirst_datelast_*
79;        6) The maximum gap between 2 consecutive files must be one
80;        day.
81;
82;
83; EXAMPLE:res = rseries_ncdf('sozotaux',920501,930410,'AA5','5d')
84;
85; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
86;                      Apr 23 2001
87;-
88;------------------------------------------------------------
89;------------------------------------------------------------
90;------------------------------------------------------------
91FUNCTION rseries_ncdf, var, date1, date2, expin, freqin, CENTURY = century, NOSTRUCT = nostruct, GRIDTYPE = gridtype, _extra = ex
92@common
93; name of the file: exp_freq_datefirst_datelast_*
94  date1 = (long(date1))[0]
95  date2 = (long(date2))[0]
96  if date2 LT date1 then return, report('date2 must be larger than date1')
97  if n_elements(expin) EQ 0 then exp = varexp ELSE BEGIN
98    IF (strpos(expin, '_'))[0] NE -1 AND n_elements(freqin) EQ 0 THEN BEGIN
99      exp = strmid(expin, 0, strpos(expin, '_'))
100      freqin = strmid(expin, strpos(expin, '_')+1)
101    ENDIF ELSE exp = expin
102    varexp = exp
103  ENDELSE
104  if n_elements(freqin) EQ 0 then freq = '5d' ELSE freq = freqin
105;------------------------------------------------
106; determinaition of the filename beginning with exp_freq_*
107;------------------------------------------------
108  IF keyword_set(gridtype) EQ 0 THEN gridtype = ''
109; determination OF datefirst and datelast
110; list of the files beginning by iodir+exp+sep+freq
111  sep = '_'
112  possiblenames = iodir+exp+sep+freq+'*'+gridtype+'.nc'
113  possiblenames = findfile(possiblenames)
114  if possiblenames[0] EQ '' then return, report('nofilename :'+iodir+exp+sep+freq+'*'+gridtype+'.nc'+' found')
115  npos = n_elements(possiblenames)
116; list of the datefirst and datelast of the files beginning by
117; iodir+exp+sep+freq
118  datefirst = lonarr(npos)
119  datelast = lonarr(npos)
120  sepshift = n_elements(str_sep(exp, sep))-1
121  for i = 0, npos-1 do BEGIN
122    separate = str_sep(possiblenames[i], sep)
123    datefirst[i] = separate[2+sepshift]
124    datelast[i] = separate[3+sepshift]
125  ENDFOR
126; selection of the files for which datefirst le date1 and datelast ge
127; date1
128  CASE strmid(freqin, 0, 1, /rever) OF
129    'm':div = 100
130    'y':div = 10000
131    ELSE:div = 1
132  ENDCASE   
133  goodfile = where(datefirst le date1/div and datelast GE date1/div)
134  if goodfile[0] EQ -1 then return, report('filename :'+iodir+exp+sep+freq+'*'+gridtype+'.nc'+' not found with the dates containing '+strtrim(date1, 1))
135  datefirst = datefirst[goodfile]
136  datelast = datelast[goodfile]
137; how many different pairs of dates are present ?
138; find the unique pairs of dates
139  pairofdate = dcomplex(datefirst, datelast)
140  pairofdate = pairofdate[uniq(pairofdate, sort(pairofdate))]
141  fileok = ''
142  datefirstok = 0L
143  datelastok = 0L
144; loop on the files for which date1 is between datefirst and datelast
145  for p = 0, n_elements(pairofdate)-1 do begin
146; name of the file: exp_freq_datefirst_datelast_*
147    datefirst = long(pairofdate[p])
148    datelast = long(imaginary(pairofdate[p]))
149; determination OF the possible filenames
150; test if the year format is i1, i2, i3, i4, i5 or i6
151    namebegin = iodir+exp+sep+freq+sep+strtrim(datefirst, 1)+sep+strtrim(datelast, 1)+'*'+gridtype+'.nc'
152    namebegin = findfile(namebegin)
153    possiblenames = ''
154    if namebegin[0] NE '' then possiblenames = [possiblenames, namebegin]
155    if datefirst LT 1e9 then begin
156      for i = 1, 5-(datefirst GE 1e5)-(datefirst GE 1e6) $
157        -(datefirst GE 1e7)-(datefirst GE 1e8) do begin
158        zero = string(replicate(byte('0'), i))
159        namebegin = iodir+exp+sep+freq+sep $
160          + zero+strtrim(datefirst, 1) +sep+ zero+strtrim(datelast, 1) +'*'+gridtype+'.nc'
161        namebegin = findfile(namebegin)
162        if namebegin[0] NE '' then possiblenames = [possiblenames, namebegin]
163      endfor
164    ENDIF
165    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')
166    possiblenames = possiblenames[1:n_elements(possiblenames)-1]
167; determination OF the filenames
168    ncdf_control, 0, /noverbose
169    i = 0
170    repeat BEGIN
171      cdfid = ncdf_open(possiblenames[i])
172;         print,possiblenames[i],var
173      test = ncdf_varid(Cdfid, var)
174      ncdf_close, cdfid
175      i = i+1
176    endrep until test NE -1 OR i EQ n_elements(possiblenames)
177    ncdf_control, 0, /verbose
178    if test NE -1 then begin
179      fileok = [fileok, possiblenames[i-1]]
180      datefirstok = [datefirstok, datefirst]
181      datelastok = [datelastok, datelast]
182    endif
183  ENDFOR
184  numoffileok = n_elements(fileok)-1
185  if numoffileok EQ 0 then return, report('the variable '+var+' was not fond in the file: '+possiblenames)
186  fileok = fileok[1:numoffileok]
187  datefirstok = datefirstok[1: numoffileok]
188  datelastok = datelastok[1: numoffileok]
189  if numoffileok NE 1 then BEGIN
190; more than one file containts the variable var and has the date1
191; between datefirst and datelast...
192; Is there any file FOR which date2 LE datelast?
193    goodfile = where(date2 LE datelastok)
194    if goodfile[0] NE -1 then begin
195      fileok = fileok[goodfile]
196      datefirstok = datefirstok[goodfile]
197      datelastok = datelastok[goodfile]
198; we choose the file which has the smallest datelast
199      bestfile = sort(datelastok)
200      filename = fileok[bestfile[0]]
201      datefirst = datefirstok[bestfile[0]]
202      datelast = datelastok[bestfile[0]]
203    ENDIF ELSE BEGIN
204; we choose the file which has the bigest datelast
205      bestfile = sort(datelastok)
206      filename = fileok[bestfile[numoffileok-1]]
207      datefirst = datefirstok[bestfile[numoffileok-1]]
208      datelast = datelastok[bestfile[numoffileok-1]]
209    ENDELSE
210  ENDIF ELSE BEGIN
211    filename = fileok[0]
212    datefirst = datefirstok[0]
213    datelast = datelastok[0]
214  ENDELSE
215; now we have find the good file with the good dates and the godd
216; variable inside.
217; we get the number of dimension of the variable
218  cdfid = ncdf_open(filename)
219  varcontient = ncdf_varinq(cdfid, var)
220  ncdf_close, cdfid
221;
222;------------------------------------------------
223; reading of the variables
224;------------------------------------------------
225;
226  if NOT keyword_set(century) then century = 0l
227  IF chkstru(ex, 'filename') THEN ex.filename = filename
228  CASE strmid(freqin, 0, 1, /rever) OF
229    'm':datelast = 100*datelast + daysinmonth(datelast MOD 100, datelast / 100)
230    'y':div = 10000*datelast+360+(5+leapyr(datelast+century))*(key_caltype EQ 'greg')
231    ELSE:div = 1
232  ENDCASE   
233  if date2 GT datelast THEN BEGIN
234; if we need to read more than one file,
235; first we read the first file
236    res1 = read_ncdf(var, date1+century*1000000L, datelast[0]+century*1000000L $
237                     , filename = filename, /nostruct, _extra = ex)
238    time1 = time                ; store the first part of the calendar
239    jpt1 = jpt                  ; store the number time steps already read
240
241; and after we call again rseries_ncdf
242; newdatefirst is defined as datelast+1day
243    newdatefirst = jul2date(date2jul(datelast)+1)
244    res2 = rseries_ncdf(var, newdatefirst, date2, exp, freq, /NOSTRUCT, _extra = ex)
245;
246;
247    CASE 1 OF
248      res1[0] EQ -1 AND res2[0] EQ -1: return, -1
249      res1[0] EQ -1:res = temporary(res2)
250      res2[0] EQ -1:BEGIN
251        time = time1
252        jpt = jpt1
253        res = temporary(res1)
254      END
255      ELSE:BEGIN
256 ; we glue the result of the first read and the result of the new call
257; to rseries_ncdf
258        time = [time1, time]
259        jpt = jpt1+jpt
260        case varcontient.ndims OF
261          1:res = [temporary(res1), temporary(res2)]
262          2:res = [ [temporary(res1)], [temporary(res2)] ]
263          3:res = [ [ [temporary(res1)] ], [ [temporary(res2)] ] ]
264          4:BEGIN
265            res = [(temporary(res1))[*], (temporary(res2))[*]]
266            grille, mask, glam, gphi, gdep, nx, ny, nz
267            res = reform(res, nx, ny, nz, jpt, /over)
268          END
269        ENDCASE
270      ENDELSE
271    ENDCASE
272  endif ELSE BEGIN
273    res = read_ncdf(var, date1+century*1000000L, date2+century*1000000L $
274                    , filename = filename, /nostruct, _extra = ex)
275  ENDELSE
276;------------------------------------------------
277  if keyword_set(nostruct) then return, res $
278  ELSE return, {arr:res, grid:vargrid, units:varunit, experiment:varexp, name:varname}
279
280end
Note: See TracBrowser for help on using the repository browser.