source: trunk/procs/nc_read.pro @ 92

Last change on this file since 92 was 92, checked in by ericg, 16 years ago

Misc EG updates following merge with new nc_read.pro (r86)
File size: 5.2 KB
Line 
1;--------------------------- 
2; Reading any netcdf file
3;--------------------------- 
4FUNCTION nc_read, file_name, var_name, ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, NO_MEAN = no_mean
5
6; arguments = file_name, varname
7; ncdf_db=<location>:<path> or just <path>
8; mot clef = iodir time_1 time_2
9
10@common
11@com_eg
12 
13; inits
14   IF debug_w THEN print, '   '
15   IF debug_w THEN print, '   ENTER nc_read...'
16
17   IF NOT keyword_set(TIME_1) THEN time_1 =  1
18   IF NOT keyword_set(TIME_2) THEN time_2 =  time_1
19   IF debug_w THEN print, '    key_yreverse : ', key_yreverse 
20;
21
22; decide which subdomain of data to read
23   IF debug_w THEN print, '     keyword_set(ALL_DATA) : ', keyword_set(ALL_DATA)
24   IF keyword_set(ALL_DATA) THEN BEGIN
25      tout = 1
26   ENDIF ELSE tout = 0
27   CASE vert_type OF
28      'z' : zindex = 0
29      'level' : zindex = 1
30      ELSE: zindex = 0
31   ENDCASE
32
33; local directory
34   IF strpos(ncdf_db, ':') GE 1 THEN directory = (strsplit(ncdf_db, ':', /EXTRACT))[1] $
35    ELSE directory = ncdf_db
36
37; Find the variable's attribute
38   ncdf_getatt, directory+file_name, var_name, add_offset = add_offset, scale_factor = scale_factor, $
39    missing_value = missing_value, units = units, calendar = calendar, long_name = long_name
40
41; By default units for the Z axis are meters
42   IF n_elements(gdept) GT 1 THEN BEGIN
43      IF name_level NE '-' THEN $
44       ncdf_getatt, directory+file_name, name_level, units = units_depth ELSE units_depth = 'm'
45   ENDIF ELSE units_depth = ''
46
47; Check consistency between time_2 computed from param in cmdline and the max number
48; of time steps found in the file.
49   jpt_max = find_jptmax(file_name, IODIRECTORY = directory)
50   IF time_2-time_1+1 GT jpt_max AND jpt_max NE -1 THEN time_2 = time_1+jpt_max-1
51
52; Read the data with a call to read_ncdf
53   IF debug_w THEN print, '     '
54   IF debug_w THEN print, '     Reading raw data from ',  file_name
55
56   lec_data = read_ncdf(var_name, time_1-1, time_2-1, BOXZOOM = boxzoom, FILENAME = directory+file_name, $
57                        /TIMESTEP, /ADDSCL_BEFORE, TOUT = tout, /NOSTRUCT, /CONT_NOFILL, GRID = vargrid, $
58                        ZINVAR = zinvar, ZINDEX = zindex)
59   field_dim = (size(lec_data))[0]
60   IF field_dim LE 0 THEN stop,  '  Something wrong happened in read_ncdf '
61
62   IF debug_w THEN print, '     '
63   IF debug_w THEN print,  '    boxzoom=', boxzoom
64   IF debug_w THEN print,  '    firstxt,firstyt,firstzt=', firstxt, firstyt, firstzt
65   IF debug_w THEN print,  '    lastxt,lastyt,lastzt=', firstxt, firstyt, firstzt   
66   IF debug_w THEN print,  '    zinvar=', zinvar
67   IF debug_w THEN print, '     ixmindta,iymindta,izmindta  =', ixmindta,iymindta,izmindta
68   IF debug_w THEN print, '     ixminmesh,iyminmesh=', ixminmesh,iyminmesh
69   IF debug_w THEN print, '     ixmaxmesh,iymaxmesh=', ixmaxmesh,iymaxmesh
70   IF debug_w THEN print, '     izminmesh,izmaxmesh=', izminmesh,izmaxmesh
71   IF debug_w THEN print, '     jpt=', jpt
72   IF debug_w THEN print, '     key_shift=', key_shift
73   IF debug_w THEN print, '     key_yreverse=',key_yreverse
74   IF debug_w THEN print, '     '
75   
76; Average data along vertical if needed and update some features
77; needed for plt (data_direc, name_suff)
78   name_suff = ''
79   data_direc = ''
80   update_data, TAB = lec_data, VNAME = var_name, BOXZOOM = boxzoom, ZUNITS = units_depth, $
81    ZINVAR = zinvar, SUFFIX = name_suff, D_DIREC = data_direc, $
82    TIME_1 = time_1, TIME_2 =  time_2, NO_MEAN = no_mean
83   field_dim = (size(lec_data))[0]
84
85; Field attributes
86   field = {name: '', data: lec_data, legend: '', units: '', origin: '', dim: 0, direc: data_direc}
87   field.name = var_name
88   field.origin = directory+file_name
89   field.legend = long_name+name_suff
90   field.units = units
91   field.dim = field_dim
92
93; replace NaN with valmask
94
95   idx_t = where (~finite(field.data))
96   IF idx_t(0) NE -1 THEN field.data(idx_t) = valmask
97
98; get valmask (might need valmask = float(string(valmask))
99
100; valmask = 1.e20
101   IF size(missing_value,  /TYPE) EQ 4 OR size(missing_value,  /TYPE) EQ 5 THEN BEGIN
102      valmask = missing_value
103; ensure valmask is positive
104      IF valmask LT 0 THEN BEGIN
105         print, '      *** Warning valmask is negative - changing sign: ', valmask
106         idx_t = where (field.data EQ valmask)
107         IF idx_t(0) NE -1 THEN field.data(idx_t) = -valmask
108         valmask = -valmask
109         idx_t=0                ; free memory
110      ENDIF
111   ENDIF
112
113; if pseudo 3d mask, read mask file and set masked points to valmask
114
115; still to do
116
117; min/max
118
119   chardim = ' - dim = '
120   FOR i = 1, (size(field.data))[0] DO BEGIN
121      chardim = chardim+string((size(field.data))[i], format = '(I5)')
122   ENDFOR 
123
124   index_test = (where (field.data EQ valmask))
125   IF index_test(0) NE -1 THEN BEGIN
126      minf = min(field.data(where (field.data NE valmask)))
127      maxf = max(field.data(where (field.data NE valmask)))
128   ENDIF ELSE BEGIN
129      minf =  min(field.data)
130      maxf =  max(field.data)
131   ENDELSE
132
133   print, '      = ',field.legend, '    [min/max = ',minf , maxf,'  ', field.units,' - masked values = ',valmask, chardim, ']'
134
135   IF keyword_set(key_yreverse) THEN print, '    key_yreverse active : ',  key_yreverse
136
137   IF debug_w THEN print,  '  ...EXIT nc_read'
138   IF debug_w THEN print,  '  '
139   
140   return, field
141
142END
Note: See TracBrowser for help on using the repository browser.