1 | ;--------------------------- |
---|
2 | ; Reading ORCA netcdf files |
---|
3 | ; EG 24/2/99 |
---|
4 | ;--------------------------- |
---|
5 | FUNCTION nc_read, file_name, var_name, ncdf_db, TIME_1 = time_1, TIME_2 = time_2, ALL_DATA = all_data, NO_MEAN = no_mean |
---|
6 | |
---|
7 | ; arguments = file_name, varname |
---|
8 | ; ncdf_db=<location>:<path> or just <path> |
---|
9 | ; mot clef = iodir time_1 time_2 |
---|
10 | |
---|
11 | @common |
---|
12 | @com_eg |
---|
13 | |
---|
14 | ;; print, 'key_yreverse : ', key_yreverse |
---|
15 | |
---|
16 | ; inits |
---|
17 | IF debug_w THEN print, 'DEBUG: Entering nc_read...' |
---|
18 | IF NOT keyword_set(TIME_1) THEN time_1 = 1 |
---|
19 | IF NOT keyword_set(TIME_2) THEN time_2 = time_1 |
---|
20 | ; |
---|
21 | ; decide which subdomain of data to read |
---|
22 | ; |
---|
23 | |
---|
24 | IF debug_w THEN print, 'keyword_set(ALL_DATA) : ', keyword_set(ALL_DATA) |
---|
25 | |
---|
26 | IF keyword_set(ALL_DATA) THEN BEGIN |
---|
27 | print, ' Reading whole domain' |
---|
28 | firstx = 0 |
---|
29 | firsty = 0 |
---|
30 | firstz = 0 |
---|
31 | nx = jpi |
---|
32 | ny = jpj |
---|
33 | nz = jpk |
---|
34 | lastx = jpi-1 |
---|
35 | lasty = jpj-1 |
---|
36 | lastz = jpk-1 |
---|
37 | ;Rajout MK |
---|
38 | ;Trouver une meilleure place |
---|
39 | ixminmesh = 0 |
---|
40 | iyminmesh = 0 |
---|
41 | ixmaxmesh = jpi-1 |
---|
42 | iymaxmesh = jpj-1 |
---|
43 | ;Fin rajout |
---|
44 | ENDIF ELSE BEGIN |
---|
45 | grille,mask,glam,gphi,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz |
---|
46 | mask = 1 |
---|
47 | ENDELSE |
---|
48 | |
---|
49 | IF debug_w THEN print, 'nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz' |
---|
50 | IF debug_w THEN print, nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz |
---|
51 | |
---|
52 | IF debug_w THEN print, 'key_offset =', key_offset |
---|
53 | IF debug_w THEN print, 'jpi, jpj, jpk =', jpi, jpj, jpk |
---|
54 | |
---|
55 | |
---|
56 | ; define reading boundaries |
---|
57 | |
---|
58 | x_count = nx |
---|
59 | y_count = ny |
---|
60 | z_count = nz |
---|
61 | |
---|
62 | IF debug_w THEN print, 'nx,ny,nz =', nx,ny,nz |
---|
63 | |
---|
64 | ; force izmaxmesh to lastz |
---|
65 | jpk = lastz+1 |
---|
66 | |
---|
67 | ; dealing with longitude periodicity |
---|
68 | IF x_count NE jpi THEN BEGIN |
---|
69 | key_shift_read = 0 |
---|
70 | ENDIF ELSE key_shift_read = key_shift |
---|
71 | IF debug_w THEN print, 'key_shift_read=', key_shift_read |
---|
72 | |
---|
73 | x_offset = key_offset[0]+ firstx+(key_shift_read-key_shift) |
---|
74 | y_offset = key_offset[1]+ firsty |
---|
75 | z_offset = key_offset[2]+ firstz |
---|
76 | IF debug_w THEN print, 'offset =', x_offset, y_offset, z_offset |
---|
77 | |
---|
78 | ; pour trouver un fichier, tu as |
---|
79 | ; findfile (tres pratique) et dialog_pickfile (interactif) |
---|
80 | |
---|
81 | ; pour verifier si il y a une variable je fais |
---|
82 | |
---|
83 | ; inq=ncdf_inquire(cdfid) |
---|
84 | ; for varid=0,inq.nvars-1 do BEGIN |
---|
85 | ; varinq=ncdf_varinq(cdfid,varid) |
---|
86 | ; if varinq.name eq nom then goto, variabletrouvee |
---|
87 | ; endfor |
---|
88 | ; return, -1 |
---|
89 | ; variabletrouvee: |
---|
90 | |
---|
91 | ; local directory |
---|
92 | |
---|
93 | IF strpos(ncdf_db, ':') GE 1 THEN directory = (str_sep(ncdf_db, ':'))[1] $ |
---|
94 | ELSE directory = ncdf_db |
---|
95 | |
---|
96 | ; existence of file |
---|
97 | |
---|
98 | list = findfile(directory+file_name, count = nb_file) |
---|
99 | |
---|
100 | IF nb_file EQ 0 THEN BEGIN |
---|
101 | print, '' |
---|
102 | print, ' ***** file ', file_name, ' not found in ', $ |
---|
103 | directory |
---|
104 | print, '' |
---|
105 | field = {data: -1.0} |
---|
106 | return, field |
---|
107 | ENDIF |
---|
108 | |
---|
109 | ; ouverture fichier netCDF + contenu |
---|
110 | cdfid=ncdf_open(directory+file_name) |
---|
111 | inq=ncdf_inquire(cdfid) |
---|
112 | ; que contient la variable |
---|
113 | varid = ncdf_varid(cdfid, var_name) |
---|
114 | |
---|
115 | name_suff = '' |
---|
116 | |
---|
117 | IF varid EQ -1 THEN BEGIN |
---|
118 | print, '' |
---|
119 | print, ' ***** field ', var_name, ' not found in file ', $ |
---|
120 | file_name |
---|
121 | print, '' |
---|
122 | field = {data: -1.0} |
---|
123 | return, field |
---|
124 | ENDIF |
---|
125 | varinq=ncdf_varinq(cdfid, var_name) |
---|
126 | |
---|
127 | ; test sur la dimension |
---|
128 | err_mess = '' |
---|
129 | field_dim = n_elements(varinq.dim) |
---|
130 | |
---|
131 | ; get unlimited record variable |
---|
132 | IF inq.recdim NE -1 THEN BEGIN |
---|
133 | ncdf_diminq, cdfid, inq.recdim, name_time, nb_time |
---|
134 | ;;ncdf_varget, cdfid, inq.recdim, time_array |
---|
135 | ;;nb_time = (size(time_array))(1) |
---|
136 | ENDIF ELSE BEGIN |
---|
137 | ; Look for a potential record dimension |
---|
138 | IF debug_w THEN print, ' Warning : no unlimited record in netCDF file' |
---|
139 | ; Look for the names of all dimensions and all variables |
---|
140 | list_dims = ncdf_listdims(cdfid) |
---|
141 | list_vars = ncdf_listvars(cdfid) |
---|
142 | ; Propose all the names and make the user choose the right one |
---|
143 | ; for the record dimension |
---|
144 | print, 'In the file ', file_name, ', the dimensions are named :' |
---|
145 | print, list_dims |
---|
146 | print, 'In the file ', file_name, ', the variables are named :' |
---|
147 | print, list_vars |
---|
148 | name_time = xquestion('Among the preceding names, which one could refer to the record dimension ?', /chkwid) |
---|
149 | IF name_time NE '-' THEN BEGIN |
---|
150 | dimidt = ncdf_dimid(cdfid, name_time) |
---|
151 | ncdf_diminq, cdfid, dimidt, name_time, nb_time |
---|
152 | print, 'You chose ', name_time, ' as a record dimension and its size is ', nb_time |
---|
153 | inq.recdim = dimidt |
---|
154 | ENDIF ELSE BEGIN |
---|
155 | print, 'No record dimension considered in the file' |
---|
156 | nb_time = 0 |
---|
157 | ENDELSE |
---|
158 | |
---|
159 | IF varinq.dim[varinq.ndims-1] NE dimidt THEN STOP, $ |
---|
160 | 'Post_it cannot handle variables whose record dimension is not the last one' |
---|
161 | |
---|
162 | ; ncdf_diminq,cdfid,(n_elements(varinq.dim)-1), name_time, nb_time |
---|
163 | ; dimidl = ncdf_dimid(cdfid, name_time) |
---|
164 | ; ncdf_diminq,cdfid,dimidl, name_time, nb_time |
---|
165 | ; varidl = ncdf_varid(cdfid, name_time) |
---|
166 | ; ncdf_varget, cdfid, varidl, time_array |
---|
167 | |
---|
168 | ENDELSE |
---|
169 | |
---|
170 | IF jpt GT nb_time THEN BEGIN |
---|
171 | jpt = nb_time |
---|
172 | print, 'Warning : the specified jpt does not correspond to the real time dimension in the file !' |
---|
173 | print, 'New jpt =', jpt |
---|
174 | ENDIF |
---|
175 | ; defs for @read_ncdf_varget |
---|
176 | |
---|
177 | ; key_stride = time_stride |
---|
178 | key_stride = 1 |
---|
179 | if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1] |
---|
180 | key_stride = 1l > long(key_stride) |
---|
181 | jpitotal = long(ixmaxmesh-ixminmesh+1) |
---|
182 | key_shift = long(testvar(var = key_shift)) |
---|
183 | |
---|
184 | key = long(key_shift MOD jpitotal) |
---|
185 | if key LT 0 then key = key+jpitotal |
---|
186 | ixmin = ixminmesh-ixmindta |
---|
187 | iymin = iyminmesh-iymindta |
---|
188 | firsttps = time_1-1 |
---|
189 | lasttps = time_2-1 |
---|
190 | |
---|
191 | name = varid |
---|
192 | |
---|
193 | IF debug_w THEN print, 'key=', key |
---|
194 | IF debug_w THEN print, 'jpitotal=', jpitotal |
---|
195 | IF debug_w THEN print, 'ixmindta,iymindta,izmindta =', ixmindta,iymindta,izmindta |
---|
196 | IF debug_w THEN print, 'ixminmesh,iyminmesh=', ixminmesh,iyminmesh |
---|
197 | IF debug_w THEN print, 'ixmaxmesh,iymaxmesh=', ixmaxmesh,iymaxmesh |
---|
198 | IF debug_w THEN print, 'izminmesh,izmaxmesh=', izminmesh,izmaxmesh |
---|
199 | IF debug_w THEN print, 'ixmin,iymin=', ixmin,iymin |
---|
200 | IF debug_w THEN print, 'firsttps,lasttps=', firsttps, lasttps |
---|
201 | IF debug_w THEN print, 'key_shift in nc_read=', key_shift |
---|
202 | IF debug_w THEN print, 'key_yreverse, firsty, lasty',key_yreverse, firsty, lasty |
---|
203 | |
---|
204 | CASE n_elements(varinq.dim) OF |
---|
205 | |
---|
206 | ;; fichier 2d : surface |
---|
207 | 2: BEGIN |
---|
208 | print, ' Reading ', var_name, ' (2D data) from ', file_name |
---|
209 | @read_ncdf_varget |
---|
210 | lec_data = res |
---|
211 | data_direc = 'xy' |
---|
212 | niveau = 1 |
---|
213 | END |
---|
214 | |
---|
215 | ;; fichier 3d : 2 cas = 1/ 2d espace + temps 2/ 3d espace |
---|
216 | 3: BEGIN |
---|
217 | ;; si varinq.dim contient la dim infinie (no 3) -> temps |
---|
218 | dim_3 = '3d' |
---|
219 | IF nb_time GE 1 THEN dim_3 = '2d' |
---|
220 | IF dim_3 EQ '2d' THEN BEGIN |
---|
221 | |
---|
222 | IF time_2 EQ time_1 THEN BEGIN |
---|
223 | print, ' Reading ', var_name, ' (2D data - index ', strcompress(string(time_1)),') from ', file_name |
---|
224 | @read_ncdf_varget |
---|
225 | lec_data = res |
---|
226 | data_direc = 'xy' |
---|
227 | ENDIF ELSE BEGIN |
---|
228 | print, ' Reading ', var_name, ' (2D data time serie)', strcompress(string(time_1)),'-', strtrim(string(time_2), 2), ' [every ',strtrim(string(time_stride), 2), '] from ', file_name |
---|
229 | IF debug_w THEN print, 'x_offset et y_offset et time_1 :', x_offset, y_offset, time_1 |
---|
230 | @read_ncdf_varget |
---|
231 | lec_data = res |
---|
232 | data_direc = 'xyt' |
---|
233 | ENDELSE |
---|
234 | |
---|
235 | ENDIF ELSE BEGIN |
---|
236 | print, ' Reading ', var_name, ' (3D data)', ' from ', file_name |
---|
237 | @read_ncdf_varget |
---|
238 | lec_data = res |
---|
239 | data_direc = 'xyz' |
---|
240 | ENDELSE |
---|
241 | END |
---|
242 | ;; fichier 4d : volume + temps |
---|
243 | 4: BEGIN |
---|
244 | IF time_2 EQ time_1 THEN BEGIN |
---|
245 | print, ' Reading ', var_name, ' (3D data - index ', strcompress(string(time_1)),') from ', file_name |
---|
246 | ; read vertical levels |
---|
247 | IF mesh_type EQ 'atm' THEN BEGIN |
---|
248 | ncdf_diminq,cdfid,2, name_level, nb_level |
---|
249 | |
---|
250 | ; get name of vertical level from metadata |
---|
251 | file_grid_config = hom_def+'grid_config.def' |
---|
252 | spawn, 'grep -i "\ '+meshlec_type+' " '+file_grid_config, line |
---|
253 | line = strcompress(strtrim(line[0], 2)) |
---|
254 | length = strlen(line) |
---|
255 | |
---|
256 | IF length EQ 0 THEN BEGIN |
---|
257 | print, ' *** nc_read : define name_level from grid ', meshlec_type, $ |
---|
258 | ' in file ', file_grid_config |
---|
259 | name_level = '' |
---|
260 | return, -1 |
---|
261 | ENDIF ELSE BEGIN |
---|
262 | argvar = str_sep(line, ' ') |
---|
263 | name_level = argvar[4] |
---|
264 | ENDELSE |
---|
265 | dimidl = ncdf_dimid(cdfid, name_level) |
---|
266 | ncdf_diminq,cdfid,dimidl, name_level, nb_level |
---|
267 | |
---|
268 | varidl = ncdf_varid(cdfid, name_level) |
---|
269 | ; make sure name_level is in hPa |
---|
270 | ncdf_attget, cdfid, varidl, 'units', val_unit |
---|
271 | IF string(val_unit) NE 'hPa' AND string(val_unit) NE 'millibar' AND string(val_unit) NE 'mbar' THEN BEGIN |
---|
272 | print, ' vertical levels coordinate not obvious = ', name_level |
---|
273 | print, ' ... using <levels> ...' |
---|
274 | varidl = ncdf_varid(cdfid, 'levels') |
---|
275 | ENDIF |
---|
276 | |
---|
277 | ncdf_varget, cdfid, varidl, gdept |
---|
278 | gdept = gdept |
---|
279 | gdepw = gdept |
---|
280 | e3t = shift(gdept, 1)-gdept |
---|
281 | e3t[0] = e3t[1] |
---|
282 | e3w = e3t |
---|
283 | jpk = nb_level |
---|
284 | tmask = reform((reform(tmask(*, *), jpi*jpj)#replicate(1, jpk)), jpi, jpj, jpk) |
---|
285 | ENDIF |
---|
286 | IF debug_w THEN print, ' going into @read_ncdf_varget ' |
---|
287 | @read_ncdf_varget |
---|
288 | lec_data = res |
---|
289 | ; lec_data = reform(lec_data,x_count, y_count, z_count) |
---|
290 | data_direc = 'xyz' |
---|
291 | IF debug_w THEN help, lec_data |
---|
292 | |
---|
293 | ; vertical average ? |
---|
294 | IF vert_switch ge 1 AND NOT keyword_set(no_mean) THEN BEGIN |
---|
295 | old_boite = [lon1, lon2, lat1, lat2, prof1, prof2] |
---|
296 | print, ' Average in vertical domain ', vert_type, vert_mean |
---|
297 | IF mesh_type EQ 'atm' THEN BEGIN |
---|
298 | CASE atmos_msk OF |
---|
299 | 0: print, ' [atmos grid : take all points] ' ; take all points |
---|
300 | 1: BEGIN |
---|
301 | ; take ocean points only |
---|
302 | idx = where(tmask EQ 0) |
---|
303 | lec_data(idx) = 1.e20 |
---|
304 | print, ' [atmos grid : take ocean points only] ' |
---|
305 | END |
---|
306 | 2: BEGIN |
---|
307 | ; take land points only |
---|
308 | idx = where(tmask GT 0) |
---|
309 | lec_data(idx) = 1.e20 |
---|
310 | print, ' [atmos grid : take land points only] ' |
---|
311 | END |
---|
312 | ENDCASE |
---|
313 | CASE vert_type OF |
---|
314 | 'z': BEGIN ; depth range |
---|
315 | IF time_1 EQ time_2 THEN BEGIN |
---|
316 | zmean = moyenne(lec_data, 'z', boite = vert_mean, NAN =1.e20) |
---|
317 | END ELSE zmean = grossemoyenne(lec_data, 'z', boite = vert_mean, NAN =1.e20) |
---|
318 | IF vert_mean[0] EQ vert_mean[1] THEN BEGIN |
---|
319 | name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0)))), 2)+' hPa' |
---|
320 | ENDIF ELSE BEGIN |
---|
321 | name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' |
---|
322 | ENDELSE |
---|
323 | END |
---|
324 | ELSE: BEGIN ; levels range |
---|
325 | IF vert_mean[0] EQ vert_mean[1] THEN BEGIN |
---|
326 | zmean = lec_data(*, *, vert_mean[0]) |
---|
327 | name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0)))), 2)+' hPa' |
---|
328 | ENDIF ELSE BEGIN |
---|
329 | zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex, NAN =1.e20) |
---|
330 | name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' |
---|
331 | ENDELSE |
---|
332 | END |
---|
333 | ENDCASE |
---|
334 | ; tmask = reform(tmask(*, *, 0), jpi, jpj) |
---|
335 | ENDIF ELSE BEGIN |
---|
336 | ; ocean = always mask |
---|
337 | ; idx = where(tmask EQ 0) |
---|
338 | ; lec_data(idx) = 1.e20 |
---|
339 | CASE vert_type OF |
---|
340 | 'z': BEGIN |
---|
341 | IF vert_mean[0] EQ vert_mean[1] THEN BEGIN |
---|
342 | zmean = lec_data ; right depth already selected |
---|
343 | name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' |
---|
344 | ENDIF ELSE BEGIN |
---|
345 | zmean = moyenne(lec_data, 'z', boite = vert_mean) |
---|
346 | name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' |
---|
347 | ENDELSE |
---|
348 | END |
---|
349 | ELSE: BEGIN ; case level (zindex) |
---|
350 | IF vert_mean[0] EQ vert_mean[1] THEN BEGIN |
---|
351 | zmean = lec_data ; right depth already selected |
---|
352 | name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' |
---|
353 | ENDIF ELSE BEGIN |
---|
354 | zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex) |
---|
355 | name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' |
---|
356 | ENDELSE |
---|
357 | END |
---|
358 | ENDCASE |
---|
359 | ENDELSE |
---|
360 | lec_data = zmean |
---|
361 | domdef, old_boite |
---|
362 | field_dim = field_dim - 1 |
---|
363 | data_direc = 'xy' |
---|
364 | nzt = 1 |
---|
365 | firstz = 0 |
---|
366 | lastz = 0 |
---|
367 | ENDIF |
---|
368 | |
---|
369 | ENDIF ELSE BEGIN |
---|
370 | print, ' Reading ', var_name, ' (3D data time serie)', strcompress(string(time_1)),'-', strtrim(string(time_2), 2), ' from ', file_name |
---|
371 | ; read vertical levels |
---|
372 | IF mesh_type EQ 'atm' THEN BEGIN |
---|
373 | ncdf_diminq,cdfid,2, name_level, nb_level |
---|
374 | |
---|
375 | ; get name of vertical level from metadata |
---|
376 | file_grid_config = hom_def+'grid_config.def' |
---|
377 | spawn, 'grep -i "\ '+meshlec_type+' " '+file_grid_config, line |
---|
378 | line = strcompress(strtrim(line[0], 2)) |
---|
379 | length = strlen(line) |
---|
380 | |
---|
381 | IF length EQ 0 THEN BEGIN |
---|
382 | print, ' *** nc_read : define name_level from grid ', meshlec_type, $ |
---|
383 | ' in file ', file_grid_config |
---|
384 | name_level = '' |
---|
385 | return, -1 |
---|
386 | ENDIF ELSE BEGIN |
---|
387 | argvar = str_sep(line, ' ') |
---|
388 | name_level = argvar[1] |
---|
389 | ENDELSE |
---|
390 | dimidl = ncdf_dimid(cdfid, name_level) |
---|
391 | ncdf_diminq,cdfid,dimidl, name_level, nb_level |
---|
392 | |
---|
393 | varidl = ncdf_varid(cdfid, name_level) |
---|
394 | ; make sure name_level is in hPa |
---|
395 | ncdf_attget, cdfid, varidl, 'units', val_unit |
---|
396 | IF string(val_unit) NE 'hPa' AND string(val_unit) NE 'millibar' AND string(val_unit) NE 'mbar' THEN BEGIN |
---|
397 | print, ' vertical levels coordinate not obvious' |
---|
398 | print, ' ... using <levels> ...' |
---|
399 | varidl = ncdf_varid(cdfid, 'levels') |
---|
400 | ENDIF |
---|
401 | ncdf_varget, cdfid, varidl, gdept |
---|
402 | gdepw = gdept |
---|
403 | e3t = shift(gdept, 1)-gdept |
---|
404 | e3t[0] = e3t[1] |
---|
405 | e3w = e3t |
---|
406 | jpk = nb_level |
---|
407 | tmask = reform((reform(tmask(*, *), jpi*jpj)#replicate(1, jpk)), jpi, jpj, jpk) |
---|
408 | ENDIF |
---|
409 | @read_ncdf_varget |
---|
410 | lec_data = res |
---|
411 | data_direc = 'xyzt' |
---|
412 | |
---|
413 | ; vertical average ? |
---|
414 | IF vert_switch ge 1 AND NOT keyword_set(no_mean) THEN BEGIN |
---|
415 | old_boite = [lon1, lon2, lat1, lat2, prof1, prof2] |
---|
416 | print, ' Average in vertical domain ', vert_type, vert_mean |
---|
417 | IF mesh_type EQ 'atm' THEN BEGIN |
---|
418 | CASE atmos_msk OF |
---|
419 | 0: print, ' [take all points] ' ; take all points |
---|
420 | 1: BEGIN |
---|
421 | ; take ocean points only |
---|
422 | idx = where(tmask EQ 0) |
---|
423 | lec_data(idx) = 1.e20 |
---|
424 | print, ' [take ocean points only] ' |
---|
425 | END |
---|
426 | 2: BEGIN |
---|
427 | ; take land points only |
---|
428 | idx = where(tmask GT 0) |
---|
429 | lec_data(idx) = 1.e20 |
---|
430 | print, ' [take land points only] ' |
---|
431 | END |
---|
432 | ENDCASE |
---|
433 | CASE vert_type OF |
---|
434 | 'z': BEGIN ; levels |
---|
435 | |
---|
436 | IF vert_mean[0] EQ vert_mean[1] THEN BEGIN |
---|
437 | name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0)))), 2)+' hPa' |
---|
438 | ENDIF ELSE BEGIN |
---|
439 | zmean = grossemoyenne(lec_data, 'z', boite = vert_mean) |
---|
440 | name_suff = ' averaged in '+vert_type+'['+strtrim(string(vert_mean(0)), 2)+','+strtrim(string(vert_mean(1)), 2)+']' |
---|
441 | ENDELSE |
---|
442 | END |
---|
443 | ELSE: BEGIN ; levels |
---|
444 | IF vert_mean[0] EQ vert_mean[1] THEN BEGIN |
---|
445 | zmean = lec_data(*, *, vert_mean[0]) |
---|
446 | name_suff = ' at '+vert_type+' '+strtrim(string(long(vert_mean(0))), 2)+' ['+strtrim(string(gdept(vert_mean(0))), 2)+' hPa]' |
---|
447 | ENDIF ELSE BEGIN |
---|
448 | zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex) |
---|
449 | name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' |
---|
450 | ENDELSE |
---|
451 | END |
---|
452 | ENDCASE |
---|
453 | ; tmask = reform(tmask(*, *, 0), jpi, jpj) |
---|
454 | ENDIF ELSE BEGIN |
---|
455 | ; ocean = always mask |
---|
456 | ; idx = where(tmask EQ 0) |
---|
457 | ; lec_data(idx) = 1.e20 |
---|
458 | CASE vert_type OF |
---|
459 | 'z': BEGIN |
---|
460 | IF vert_mean[0] EQ vert_mean[1] THEN BEGIN |
---|
461 | zmean = lec_data |
---|
462 | name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' |
---|
463 | ENDIF ELSE BEGIN |
---|
464 | zmean = grossemoyenne(lec_data, 'z', boite = vert_mean) |
---|
465 | name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' |
---|
466 | ENDELSE |
---|
467 | END |
---|
468 | ELSE: BEGIN |
---|
469 | IF vert_mean[0] EQ vert_mean[1] THEN BEGIN |
---|
470 | zmean = lec_data |
---|
471 | name_suff = ' at '+strtrim(string(long(gdept(vert_mean(0))+.1)), 2)+' m' |
---|
472 | ENDIF ELSE BEGIN |
---|
473 | zmean = moyenne(lec_data, 'z', boite = vert_mean, /zindex, NAN =1.e20) |
---|
474 | name_suff = ' averaged in '+vert_type+'['+strtrim(string(long(vert_mean(0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']' |
---|
475 | ENDELSE |
---|
476 | END |
---|
477 | ENDCASE |
---|
478 | ENDELSE |
---|
479 | lec_data = zmean |
---|
480 | domdef, old_boite |
---|
481 | field_dim = field_dim - 1 |
---|
482 | data_direc = 'xyt' |
---|
483 | nzt = 1 |
---|
484 | firstz = 0 |
---|
485 | lastz = 0 |
---|
486 | ENDIF |
---|
487 | ENDELSE |
---|
488 | END |
---|
489 | ;; erreur si dim > 4 |
---|
490 | ELSE: BEGIN |
---|
491 | err_mess = ' *** nc_read : ERROR dimension > 4' |
---|
492 | lec_data = -1.0 |
---|
493 | END |
---|
494 | ENDCASE |
---|
495 | |
---|
496 | ; scaling ? |
---|
497 | FOR i = 0, varinq.natts-1 DO BEGIN |
---|
498 | att_txt = ncdf_attname(cdfid, varid, i) |
---|
499 | IF att_txt EQ 'scale_factor' THEN BEGIN |
---|
500 | ncdf_attget, cdfid, varid, att_txt, valscale |
---|
501 | lec_data = lec_data*valscale |
---|
502 | ENDIF |
---|
503 | ENDFOR |
---|
504 | |
---|
505 | ; Attribut du champ |
---|
506 | |
---|
507 | field = {name: '', data: lec_data, legend: '', units: '', origin: '', dim: 0, direc: data_direc} |
---|
508 | field.name = var_name |
---|
509 | field.origin = directory+file_name |
---|
510 | ; |
---|
511 | ; get long name |
---|
512 | ; result = '???' |
---|
513 | FOR i = 0, varinq.natts-1 DO BEGIN |
---|
514 | att_txt = ncdf_attname(cdfid, varid, i) |
---|
515 | IF att_txt EQ 'long_name' OR att_txt EQ 'title' THEN ncdf_attget, cdfid, varid, att_txt, result |
---|
516 | ENDFOR |
---|
517 | field.legend = strtrim(string(result),1)+name_suff |
---|
518 | |
---|
519 | ; get unit |
---|
520 | ; if it exists |
---|
521 | isunits = ncdf_attinq(cdfid, varid, 'units') |
---|
522 | IF isunits.datatype NE 'UNKNOWN' THEN BEGIN |
---|
523 | ncdf_attget, cdfid, varid, 'units', result |
---|
524 | field.units = strtrim(string(result),1) |
---|
525 | ENDIF ELSE BEGIN |
---|
526 | print, 'No units for the variable ', field.name |
---|
527 | print, ' Or units were forgotten in the file !' |
---|
528 | ENDELSE |
---|
529 | field.dim = field_dim |
---|
530 | |
---|
531 | |
---|
532 | ; get valmask (might need valmask = float(string(valmask)) |
---|
533 | |
---|
534 | valmask = 1.e20 |
---|
535 | FOR i = 0, varinq.natts-1 DO BEGIN |
---|
536 | att_txt = ncdf_attname(cdfid, varid, i) |
---|
537 | IF att_txt EQ 'missing_value' OR att_txt EQ 'mask value' OR att_txt EQ '_FillValue' THEN ncdf_attget, cdfid, varid, att_txt, valmask |
---|
538 | ENDFOR |
---|
539 | |
---|
540 | ; set valmask to 1.e20 |
---|
541 | |
---|
542 | ; ensure valmask is positive |
---|
543 | |
---|
544 | IF valmask LT 0 THEN BEGIN |
---|
545 | print, ' *** Warning valmask is negative - changing sign: ', valmask |
---|
546 | idx_t = where (field.data EQ valmask) |
---|
547 | IF idx_t(0) NE -1 THEN field.data(idx_t) = -valmask |
---|
548 | valmask = -valmask |
---|
549 | idx_t=0 ; free memory |
---|
550 | ENDIF |
---|
551 | |
---|
552 | ; min/max |
---|
553 | |
---|
554 | chardim = ' - dim = ' |
---|
555 | FOR i = 1, (size(field.data))[0] DO BEGIN |
---|
556 | chardim = chardim+string((size(field.data))[i], format = '(I5)') |
---|
557 | ENDFOR |
---|
558 | |
---|
559 | index_test = (where (field.data NE valmask)) |
---|
560 | IF index_test(0) NE -1 THEN BEGIN |
---|
561 | minf = min(field.data(where (field.data NE valmask))) |
---|
562 | maxf = max(field.data(where (field.data NE valmask))) |
---|
563 | ENDIF ELSE BEGIN |
---|
564 | minf = min(field.data) |
---|
565 | maxf = max(field.data) |
---|
566 | ENDELSE |
---|
567 | |
---|
568 | print, ' = ',field.legend, ' [min/max = ',minf , maxf,' ', field.units,' - masked values = ',valmask, chardim, ']' |
---|
569 | |
---|
570 | |
---|
571 | ncdf_close, cdfid |
---|
572 | |
---|
573 | ;key_offset = key_offset_orig |
---|
574 | ;jpi = jpi_orig |
---|
575 | ;jpj = jpj_orig |
---|
576 | ;jpk = jpk_orig |
---|
577 | |
---|
578 | IF keyword_set(key_yreverse) THEN print, ' key_yreverse active : ', key_yreverse |
---|
579 | |
---|
580 | IF debug_w THEN print, 'DEBUG: Exiting nc_read...' |
---|
581 | |
---|
582 | return, field |
---|
583 | |
---|
584 | |
---|
585 | END |
---|