[2] | 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 | ; |
---|
[10] | 23 | |
---|
[2] | 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 | |
---|
[27] | 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 |
---|
[2] | 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) |
---|
[27] | 111 | inq=ncdf_inquire(cdfid) |
---|
[2] | 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 |
---|
[27] | 125 | varinq=ncdf_varinq(cdfid, var_name) |
---|
[2] | 126 | |
---|
| 127 | ; test sur la dimension |
---|
| 128 | err_mess = '' |
---|
[27] | 129 | field_dim = n_elements(varinq.dim) |
---|
[2] | 130 | |
---|
| 131 | ; get unlimited record variable |
---|
[27] | 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 |
---|
[2] | 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 |
---|
[27] | 153 | inq.recdim = dimidt |
---|
[2] | 154 | ENDIF ELSE BEGIN |
---|
| 155 | print, 'No record dimension considered in the file' |
---|
| 156 | nb_time = 0 |
---|
| 157 | ENDELSE |
---|
| 158 | |
---|
[27] | 159 | IF varinq.dim[varinq.ndims-1] NE dimidt THEN STOP, $ |
---|
[2] | 160 | 'Post_it cannot handle variables whose record dimension is not the last one' |
---|
| 161 | |
---|
[27] | 162 | ; ncdf_diminq,cdfid,(n_elements(varinq.dim)-1), name_time, nb_time |
---|
[2] | 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 | |
---|
[27] | 204 | CASE n_elements(varinq.dim) OF |
---|
[2] | 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 |
---|
[27] | 217 | ;; si varinq.dim contient la dim infinie (no 3) -> temps |
---|
[2] | 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 |
---|
[4] | 456 | ; idx = where(tmask EQ 0) |
---|
| 457 | ; lec_data(idx) = 1.e20 |
---|
[2] | 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 ? |
---|
[27] | 497 | FOR i = 0, varinq.natts-1 DO BEGIN |
---|
[2] | 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 = '???' |
---|
[27] | 513 | FOR i = 0, varinq.natts-1 DO BEGIN |
---|
[2] | 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 |
---|
[27] | 535 | FOR i = 0, varinq.natts-1 DO BEGIN |
---|
[2] | 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 | |
---|
[10] | 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 | |
---|
[2] | 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 |
---|