[2] | 1 | ;+ |
---|
[231] | 2 | ; |
---|
[150] | 3 | ; @file_comments |
---|
| 4 | ; Reading function for the file net_cdf. |
---|
[495] | 5 | ; This program is less universal than <pro>ncdf_lec</pro> (it appeal to |
---|
| 6 | ; declared variables in <pro>common</pro>) but it is very easier to be used. |
---|
[370] | 7 | ; It considerate |
---|
[150] | 8 | ; the declaration of the different zooms which have been defined |
---|
| 9 | ; (ixminmesh...premierx...), the declaration of the variable key_shift... |
---|
[495] | 10 | ; To put it in a nutshell, the result of read_ncdf can be directly used in |
---|
[370] | 11 | ; <pro>plt</pro> ... |
---|
| 12 | ; |
---|
[226] | 13 | ; This is also this program which is used by default in our reading widgets. |
---|
[2] | 14 | ; |
---|
[150] | 15 | ; @categories |
---|
[157] | 16 | ; Reading |
---|
[226] | 17 | ; |
---|
[163] | 18 | ; @param NAME {in}{required}{type=string} |
---|
| 19 | ; It define the field to be read. |
---|
[2] | 20 | ; |
---|
[150] | 21 | ; @param BEGINNING {in}{required} |
---|
| 22 | ; Relative with the time axis. |
---|
| 23 | ; These can be |
---|
[365] | 24 | ; - 2 dates of the type yyyymmdd and in this case, we select data |
---|
[150] | 25 | ; which are included between these two dates. |
---|
[365] | 26 | ; - 2 indexes which define between which and which time steps we have |
---|
[163] | 27 | ; to extract the temporal dimension. |
---|
[2] | 28 | ; |
---|
[365] | 29 | ; @param ENDING {in}{required} |
---|
[150] | 30 | ; Relative with the time axis. |
---|
| 31 | ; See BEGINNING. |
---|
[226] | 32 | ; |
---|
[174] | 33 | ; @param COMPATIBILITY {in}{optional} |
---|
| 34 | ; Useless, defined for compatibility |
---|
[226] | 35 | ; |
---|
[271] | 36 | ; @keyword ADDSCL_BEFORE {default=0}{type=scalar: 0 or 1} |
---|
[366] | 37 | ; put 1 to apply add_offset and scale factor on data before looking for |
---|
[271] | 38 | ; missing values |
---|
| 39 | ; |
---|
[457] | 40 | ; @keyword ALLRECORDS |
---|
| 41 | ; activate to read all records in a file (no need to specify any |
---|
| 42 | ; parameter other than NAME) |
---|
| 43 | ; |
---|
[226] | 44 | ; @keyword BOXZOOM |
---|
| 45 | ; Contain the boxzoom on which we have to do the reading |
---|
| 46 | ; |
---|
[174] | 47 | ; @keyword CALLITSELF {default=0}{type=scalar: 0 or 1} |
---|
[495] | 48 | ; For ROMS outputs. Use by read_ncdf itself to access auxiliary data |
---|
[366] | 49 | ; (h and zeta). |
---|
[226] | 50 | ; |
---|
[276] | 51 | ; @keyword DIREC |
---|
| 52 | ; a string used to specify the direction along which we want to make |
---|
| 53 | ; spatial and/or temporal mean. It could be: 'x' 'y' 'z' 't' 'xy' 'xz' |
---|
[365] | 54 | ; 'yz' 'xyz' 'xt' 'yt' 'zt' 'xyt' 'xzt' 'yzt' or 'xyzt' |
---|
[276] | 55 | ; |
---|
[174] | 56 | ; @keyword FILENAME {required}{type=string} |
---|
[295] | 57 | ; It contains the file's name. |
---|
[226] | 58 | ; |
---|
[174] | 59 | ; @keyword INIT {default=0}{type=scalar: 0 or 1} |
---|
[365] | 60 | ; To call automatically <pro>initncdf</pro> with filename as input argument |
---|
[366] | 61 | ; and thus redefine all the grid parameters |
---|
[226] | 62 | ; |
---|
[150] | 63 | ; @keyword GRID |
---|
[163] | 64 | ; ='[UTVWF]' to specify the type of grid. Default is (1) |
---|
[150] | 65 | ; based on the name of the file if the file ends by |
---|
| 66 | ; GRID[._][TUVFW].NC (not case sensible) or (2) T if case (1) |
---|
| 67 | ; is not found. |
---|
[226] | 68 | ; |
---|
[174] | 69 | ; @keyword TIMESTEP {default=0}{type=scalar: 0 or 1} |
---|
[495] | 70 | ; Specify that BEGINNING and ENDING refer to indexes of the time axis and not |
---|
[366] | 71 | ; to dates |
---|
[2] | 72 | ; |
---|
[174] | 73 | ; @keyword TOUT {default=0}{type=scalar: 0 or 1} |
---|
[226] | 74 | ; We activate it if we want to read the file on the whole domain without |
---|
| 75 | ; considerate the sub-domain defined by the boxzoom or |
---|
[150] | 76 | ; lon1,lon2,lat1,lat2,vert1,vert2. |
---|
[226] | 77 | ; |
---|
[174] | 78 | ; @keyword NOSTRUCT {default=0}{type=scalar: 0 or 1} |
---|
[226] | 79 | ; We activate it if we do not want that read_ncdf send back a structure |
---|
[163] | 80 | ; but only the array referring to the field. |
---|
[226] | 81 | ; |
---|
[174] | 82 | ; @keyword ZETAFILENAME {default=FILENAME}{type=string} |
---|
[238] | 83 | ; For ROMS outputs. The filename of the file where zeta variable should be read |
---|
[174] | 84 | ; |
---|
| 85 | ; @keyword ZETAZERO {default=0}{type=scalar: 0 or 1} |
---|
| 86 | ; For ROMS outputs. To define zeta to 0. instead of reading it |
---|
| 87 | ; |
---|
[336] | 88 | ; @keyword ZINVAR {type=named variable} |
---|
| 89 | ; Set this keyword to a named variable in which 1 is returned if a |
---|
| 90 | ; vertical dimension is found in the variable. Returns 0 otherwise |
---|
| 91 | ; |
---|
[150] | 92 | ; @keyword _EXTRA |
---|
[271] | 93 | ; Used to pass keywords to <pro>isafile</pro>, <pro>initncdf</pro>, |
---|
[337] | 94 | ; <pro>ncdf_gettime</pro>, <pro>ncdf_getatt</pro> and <pro>domdef</pro> |
---|
[2] | 95 | ; |
---|
[150] | 96 | ; @returns |
---|
[495] | 97 | ; Structure readable by <pro>litchamp</pro> or an array if NOSTRUCT is |
---|
[366] | 98 | ; activated. |
---|
| 99 | ; |
---|
[150] | 100 | ; @uses |
---|
[370] | 101 | ; <pro>common</pro> |
---|
[226] | 102 | ; |
---|
[150] | 103 | ; @restrictions |
---|
| 104 | ; The field must have a temporal dimension. |
---|
[226] | 105 | ; |
---|
[150] | 106 | ; @history |
---|
[157] | 107 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
[226] | 108 | ; 15/10/1999 |
---|
| 109 | ; |
---|
[150] | 110 | ; @version |
---|
[226] | 111 | ; $Id$ |
---|
[325] | 112 | ; |
---|
[2] | 113 | ;- |
---|
[327] | 114 | FUNCTION read_ncdf, name, beginning, ending, compatibility $ |
---|
[457] | 115 | , ALLRECORDS = allrecords, BOXZOOM = boxzoom, FILENAME = filename $ |
---|
[336] | 116 | , PARENTIN = parentin, TIMESTEP = timestep $ |
---|
| 117 | , ADDSCL_BEFORE = addscl_before $ |
---|
[366] | 118 | , TOUT = tout, NOSTRUCT = nostruct $ |
---|
| 119 | , CONT_NOFILL = CONT_NOFILL, INIT = init $ |
---|
[336] | 120 | , GRID = grid, CALLITSELF = callitself, DIREC = direc $ |
---|
| 121 | , ZETAFILENAME = zetafilename, ZETAZERO = zetazero $ |
---|
[501] | 122 | , ZINVAR = zinvar, ISROMS = isroms, _EXTRA = ex |
---|
[114] | 123 | ; |
---|
| 124 | compile_opt idl2, strictarrsubs |
---|
| 125 | ; |
---|
[44] | 126 | @cm_4mesh |
---|
| 127 | @cm_4data |
---|
| 128 | @cm_4cal |
---|
| 129 | IF NOT keyword_set(key_forgetold) THEN BEGIN |
---|
| 130 | @updatenew |
---|
| 131 | @updatekwd |
---|
| 132 | ENDIF |
---|
[2] | 133 | ;------------------------------------------------------------ |
---|
| 134 | ; we find the filename. |
---|
| 135 | ;------------------------------------------------------------ |
---|
[226] | 136 | ; print,filename |
---|
[2] | 137 | ; is parent a valid widget ? |
---|
[271] | 138 | IF keyword_set(parentin) THEN BEGIN |
---|
[44] | 139 | parent = long(parentin) |
---|
| 140 | parent = parent*widget_info(parent, /managed) |
---|
| 141 | ENDIF |
---|
| 142 | filename = isafile(filename = filename, IODIRECTORY = iodir, _EXTRA = ex) |
---|
[2] | 143 | ;------------------------------------------------------------ |
---|
[150] | 144 | ; Opening of the name file |
---|
[2] | 145 | ;------------------------------------------------------------ |
---|
[271] | 146 | IF size(filename, /type) NE 7 THEN $ |
---|
[336] | 147 | return, report('read_ncdf cancelled') |
---|
[44] | 148 | IF !version.OS_FAMILY EQ 'unix' THEN spawn, '\file '+filename+' > /dev/null' |
---|
| 149 | cdfid = ncdf_open(filename) |
---|
[271] | 150 | inq = ncdf_inquire(cdfid) |
---|
[2] | 151 | ;------------------------------------------------------------ |
---|
| 152 | ; we check if the variable name exists in the file. |
---|
| 153 | ;------------------------------------------------------------ |
---|
[271] | 154 | IF ncdf_varid(cdfid, name) EQ -1 THEN BEGIN |
---|
[44] | 155 | ncdf_close, cdfid |
---|
[2] | 156 | return, report('variable '+name+' !C not found in the file '+filename) |
---|
[44] | 157 | ENDIF |
---|
[271] | 158 | varinq = ncdf_varinq(cdfid, name) |
---|
| 159 | IF varinq.ndims LT 2 THEN return, report('read_ncdf cannot read scalar or 1D data') |
---|
[172] | 160 | ; look for the dimension names |
---|
[271] | 161 | dimnames = strarr(varinq.ndims) |
---|
| 162 | FOR i = 0, varinq.ndims-1 DO BEGIN |
---|
| 163 | ncdf_diminq, cdfid, varinq.dim[i], tmp, dimsize |
---|
[172] | 164 | dimnames[i] = tmp |
---|
[226] | 165 | ENDFOR |
---|
[2] | 166 | ;------------------------------------------------------------ |
---|
[150] | 167 | ; check the time axis and the debut and ending dates |
---|
[2] | 168 | ;------------------------------------------------------------ |
---|
[457] | 169 | IF n_elements(beginning) EQ 0 AND NOT keyword_set(allrecords) THEN BEGIN |
---|
[284] | 170 | beginning = 0L |
---|
| 171 | timestep = 1L |
---|
[271] | 172 | ENDIF |
---|
| 173 | ; define time and jpt |
---|
| 174 | CASE 1 OF |
---|
[365] | 175 | keyword_set(timestep):BEGIN |
---|
[284] | 176 | firsttps = long(beginning[0]) |
---|
| 177 | IF n_elements(ending) NE 0 THEN lasttps = long(ending[0]) ELSE lasttps = firsttps |
---|
[495] | 178 | IF NOT keyword_set(callitself) then BEGIN |
---|
[462] | 179 | time = ncdf_gettime(filename, cdfid, caller = 'read_ncdf', err = err, _extra = ex) |
---|
| 180 | jpt = lasttps-firsttps+1 |
---|
| 181 | IF time[0] LT 0 then time = julday(1, 1, 1) + lindgen(jpt) $ |
---|
| 182 | ELSE time = time[firsttps:lasttps] |
---|
| 183 | ENDIF ELSE jpt = lasttps-firsttps+1 |
---|
[271] | 184 | END |
---|
[365] | 185 | keyword_set(parent):BEGIN |
---|
[44] | 186 | widget_control, parent, get_uvalue = top_uvalue |
---|
| 187 | filelist = extractatt(top_uvalue, 'filelist') |
---|
| 188 | IF filelist[0] EQ 'many !' THEN filelist = filename |
---|
| 189 | currentfile = (where(filelist EQ filename))[0] |
---|
| 190 | time = (*(extractatt(top_uvalue, 'fileparameters'))[currentfile]).time_counter |
---|
[298] | 191 | date1 = date2jul(beginning[0]) |
---|
| 192 | if n_elements(ending) NE 0 then date2 = date2jul(ending[0]) ELSE date2 = date1 |
---|
[336] | 193 | firsttps = (where(abs(time - date1) LT 0.9d/86400.d))[0] ; beware of rounding errors... |
---|
[271] | 194 | lasttps = (where(abs(time - date2) LT 0.9d/86400.d))[0] |
---|
| 195 | jpt = lasttps-firsttps+1 |
---|
| 196 | END |
---|
[462] | 197 | keyword_set(allrecords) OR n_elements(beginning) EQ 0:BEGIN |
---|
[457] | 198 | time = ncdf_gettime(filename, cdfid, caller = 'read_ncdf', err = err, _extra = ex) |
---|
[495] | 199 | IF time[0] LT 0 then BEGIN |
---|
[457] | 200 | jpt = -time[0] |
---|
| 201 | time = julday(1, 1, 1) + lindgen(jpt) |
---|
| 202 | ENDIF |
---|
| 203 | firsttps = 0 |
---|
| 204 | lasttps = jpt -1 |
---|
| 205 | END |
---|
[365] | 206 | ELSE:BEGIN |
---|
[271] | 207 | time = ncdf_gettime(filename, cdfid, caller = 'read_ncdf', err = err, _extra = ex) |
---|
| 208 | IF n_elements(err) NE 0 THEN BEGIN |
---|
| 209 | dummy = report(err) |
---|
[44] | 210 | ncdf_close, cdfid |
---|
[271] | 211 | return, -1 |
---|
[44] | 212 | ENDIF |
---|
[271] | 213 | ; date1 |
---|
[298] | 214 | date1 = date2jul(beginning[0]) |
---|
[271] | 215 | ; date2 |
---|
[298] | 216 | if n_elements(ending) NE 0 then date2 = date2jul(ending[0]) ELSE date2 = date1 |
---|
[271] | 217 | ; firsttps |
---|
| 218 | firsttps = where(time GE (date1 - 0.9d/86400.d)) & firsttps = firsttps[0] |
---|
[44] | 219 | if firsttps EQ -1 THEN BEGIN |
---|
| 220 | ncdf_close, cdfid |
---|
| 221 | return, report('date 1: '+strtrim(jul2date(date1), 1)+' is not found in the time axis.') |
---|
[74] | 222 | ENDIF |
---|
[271] | 223 | ; lasttps |
---|
| 224 | lasttps = where(time LE (date2 + 0.9d/86400.d)) & lasttps = lasttps[n_elements(lasttps)-1] |
---|
| 225 | if lasttps EQ -1 THEN BEGIN |
---|
[44] | 226 | ncdf_close, cdfid |
---|
[198] | 227 | return, report('the time axis has no date before date 2: '+strtrim(jul2date(date2), 1)) |
---|
[44] | 228 | endif |
---|
| 229 | if lasttps LT firsttps then BEGIN |
---|
| 230 | ncdf_close, cdfid |
---|
[198] | 231 | return, report('the time axis has no dates between date1 and date 2: '+strtrim(jul2date(date1), 1)+' '+strtrim(jul2date(date2), 1)) |
---|
[44] | 232 | endif |
---|
[271] | 233 | time = time[firsttps:lasttps] |
---|
| 234 | jpt = lasttps-firsttps+1 |
---|
| 235 | END |
---|
| 236 | ENDCASE |
---|
[501] | 237 | IF inq.recdim EQ -1 AND jpt NE 1 THEN inq.recdim = varinq.dim[varinq.ndims-1] |
---|
[2] | 238 | ;------------------------------------------------------------ |
---|
[501] | 239 | ; we check if the variable has a vertical dimension and/or |
---|
| 240 | ; a record dimension. This is useful to define boxzoom |
---|
| 241 | ; keyword in domdef |
---|
| 242 | ;------------------------------------------------------------ |
---|
| 243 | dummy = where(varinq.dim EQ inq.recdim, tinvar) |
---|
| 244 | ; check the presence of a vertical dimension according to the |
---|
| 245 | ; number of dimensions and the presence of a record dimension |
---|
| 246 | zinvar = (varinq.ndims EQ 3 AND tinvar NE 1) OR varinq.ndims EQ 4 |
---|
| 247 | ;------------------------------------------------------------ |
---|
| 248 | ; shall we redefine the grid parameters |
---|
| 249 | ;------------------------------------------------------------ |
---|
| 250 | IF keyword_set(init) THEN initncdf, filename, _extra = ex |
---|
| 251 | ;------------------------------------------------------------ |
---|
[150] | 252 | ; Name of the grid on which the field refer to. |
---|
[2] | 253 | ;------------------------------------------------------------ |
---|
[44] | 254 | IF keyword_set(grid) THEN vargrid = strupcase(grid) ELSE BEGIN |
---|
| 255 | vargrid = 'T' ; default definition |
---|
[69] | 256 | IF finite(glamu[0]) EQ 1 THEN BEGIN |
---|
[172] | 257 | ; are we in one of the case corresponding to ROMS conventions? |
---|
| 258 | CASE 1 OF |
---|
[336] | 259 | dimnames[2 < (varinq.ndims-1)] EQ 's_w':vargrid = 'W' |
---|
[172] | 260 | dimnames[0] EQ 'xi_rho' AND dimnames[1] EQ 'eta_rho':vargrid = 'T' |
---|
| 261 | dimnames[0] EQ 'xi_u' AND dimnames[1] EQ 'eta_u' :vargrid = 'U' |
---|
| 262 | dimnames[0] EQ 'xi_v' AND dimnames[1] EQ 'eta_v' :vargrid = 'V' |
---|
| 263 | dimnames[0] EQ 'xi_psi' AND dimnames[1] EQ 'eta_psi':vargrid = 'F' |
---|
| 264 | dimnames[0] EQ 'xi_rho' AND dimnames[1] EQ 'eta_v' :vargrid = 'V' |
---|
| 265 | dimnames[0] EQ 'xi_u' AND dimnames[1] EQ 'eta_rho':vargrid = 'U' |
---|
| 266 | dimnames[0] EQ 'xi_u' AND dimnames[1] EQ 'eta_v' :vargrid = 'F' |
---|
[226] | 267 | ELSE:BEGIN |
---|
[172] | 268 | ; could we define the grid type from the file name?? |
---|
| 269 | pattern = ['GRID.', 'GRID_', 'GRID', 'UPID_', '30ID_'] |
---|
| 270 | gdtype = ['T', 'U', 'V', 'W', 'F'] |
---|
| 271 | fnametest = strupcase(filename) |
---|
| 272 | FOR i = 0, n_elements(pattern)-1 DO BEGIN |
---|
| 273 | FOR j = 0, n_elements(gdtype)-1 DO BEGIN |
---|
| 274 | substr = pattern[i]+gdtype[j] |
---|
| 275 | pos = strpos(fnametest, substr) |
---|
| 276 | IF pos NE -1 THEN $ |
---|
| 277 | vargrid = strmid(fnametest, pos+strlen(substr)-1, 1) |
---|
| 278 | ENDFOR |
---|
| 279 | ENDFOR |
---|
| 280 | END |
---|
| 281 | ENDCASE |
---|
[336] | 282 | ENDIF |
---|
[44] | 283 | ENDELSE |
---|
[2] | 284 | ;--------------------------------------------------------------- |
---|
[365] | 285 | ; redefinition of the domain |
---|
[2] | 286 | ;--------------------------------------------------------------- |
---|
[271] | 287 | IF keyword_set(tout) THEN BEGIN |
---|
[44] | 288 | nx = jpi |
---|
| 289 | ny = jpj |
---|
| 290 | nz = jpk |
---|
| 291 | firstx = 0 |
---|
| 292 | firsty = 0 |
---|
| 293 | firstz = 0 |
---|
| 294 | lastx = jpi-1 |
---|
| 295 | lasty = jpj-1 |
---|
| 296 | lastz = jpk-1 |
---|
| 297 | case strupcase(vargrid) of |
---|
| 298 | 'T':mask = tmask |
---|
| 299 | 'U':mask = umask() |
---|
| 300 | 'V':mask = vmask() |
---|
| 301 | 'W':mask = tmask |
---|
| 302 | 'F':mask = fmask() |
---|
| 303 | endcase |
---|
| 304 | ENDIF ELSE BEGIN |
---|
[336] | 305 | IF keyword_set(boxzoom) then BEGIN |
---|
| 306 | CASE N_Elements(Boxzoom) Of |
---|
| 307 | 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]] |
---|
| 308 | 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]] |
---|
| 309 | 4:bte = [Boxzoom, vert1, vert2] |
---|
| 310 | 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]] |
---|
| 311 | 6:bte = Boxzoom |
---|
| 312 | ELSE: BEGIN |
---|
[44] | 313 | ncdf_close, cdfid |
---|
| 314 | return, report('Wrong Definition of Boxzoom') |
---|
[336] | 315 | END |
---|
[44] | 316 | ENDCASE |
---|
[336] | 317 | IF zinvar EQ 1 THEN domdef, bte, GRIDTYPE = ['T', vargrid], _extra = ex $ |
---|
| 318 | ELSE domdef, bte[0:3], GRIDTYPE = ['T', vargrid], _extra = ex |
---|
[44] | 319 | ENDIF |
---|
| 320 | grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz |
---|
[226] | 321 | undefine, glam & undefine, gphi & ; We liberate some memory! |
---|
[44] | 322 | ENDELSE |
---|
[2] | 323 | ;--------------------------------------------------------------------- |
---|
[295] | 324 | ; We initialize ixmindta, iymindta if needed |
---|
[2] | 325 | ;--------------------------------------------------------------------- |
---|
[44] | 326 | if n_elements(jpidta) EQ 0 THEN jpidta = jpiglo |
---|
| 327 | if n_elements(jpjdta) EQ 0 THEN jpjdta = jpjglo |
---|
| 328 | if n_elements(jpkdta) EQ 0 THEN jpkdta = jpkglo |
---|
| 329 | if n_elements(ixmindta) EQ 0 THEN ixmindta = 0 |
---|
| 330 | if n_elements(ixmaxdta) EQ 0 then ixmaxdta = jpidta-1 |
---|
| 331 | if ixmindta EQ -1 THEN ixmindta = 0 |
---|
| 332 | IF ixmaxdta EQ -1 then ixmaxdta = jpidta-1 |
---|
| 333 | if n_elements(iymindta) EQ 0 THEN iymindta = 0 |
---|
| 334 | IF n_elements(iymaxdta) EQ 0 then iymaxdta = jpjdta-1 |
---|
| 335 | if iymindta EQ -1 THEN iymindta = 0 |
---|
| 336 | IF iymaxdta EQ -1 then iymaxdta = jpjdta-1 |
---|
| 337 | if n_elements(izmindta) EQ 0 THEN izmindta = 0 |
---|
| 338 | IF n_elements(izmaxdta) EQ 0 then izmaxdta = jpkdta-1 |
---|
| 339 | if izmindta EQ -1 THEN izmindta = 0 |
---|
| 340 | IF izmaxdta EQ -1 then izmaxdta = jpkdta-1 |
---|
[2] | 341 | ;---------------------------------------------------------------- |
---|
[150] | 342 | ; We will read the file |
---|
[2] | 343 | ;--------------------------------------------------------------- |
---|
[44] | 344 | if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1] |
---|
| 345 | key_stride = 1l > long(key_stride) |
---|
| 346 | ;--------------------------------------------------------------------- |
---|
| 347 | ;--------------------------------------------------------------------- |
---|
[2] | 348 | @read_ncdf_varget |
---|
| 349 | ;--------------------------------------------------------------------- |
---|
[44] | 350 | ;--------------------------------------------------------------------- |
---|
[271] | 351 | ; We define common (cm_4data) variables associated with the variable. |
---|
[2] | 352 | ;--------------------------------------------------------------------- |
---|
| 353 | ; varname |
---|
[172] | 354 | IF NOT keyword_set(callitself) THEN varname = name |
---|
[271] | 355 | ; varunit and others... |
---|
[337] | 356 | ncdf_getatt, cdfid, name, add_offset = add_offset, scale_factor = scale_factor, missing_value = missing_value, units = units, _extra = ex |
---|
[271] | 357 | IF NOT keyword_set(callitself) THEN varunit = units |
---|
[2] | 358 | ; vardate |
---|
[398] | 359 | ; We make a legible date |
---|
[399] | 360 | caldat, time[0], month, day, year |
---|
[44] | 361 | vardate = string(format = '(C(CMoA))', 31*(month-1))+' '+strtrim(day, 1)+', '+strtrim(year, 1) |
---|
[74] | 362 | varexp = file_basename(filename) |
---|
[150] | 363 | ; We apply the value valmask on land points. |
---|
[44] | 364 | if NOT keyword_set(cont_nofill) then begin |
---|
[271] | 365 | valmask = 1.e20 |
---|
[44] | 366 | case 1 of |
---|
[275] | 367 | varinq.ndims eq 2: earth = where(mask[*, *, 0] EQ 0) ;xy array |
---|
| 368 | varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] EQ -1:earth = where(mask EQ 0) ;xyz array |
---|
[365] | 369 | varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] NE -1:earth = where(mask[*, *, 0] EQ 0) ;xyt array |
---|
[275] | 370 | varinq.ndims eq 4: earth = where(mask EQ 0) ;xyzt array |
---|
[271] | 371 | ENDCASE |
---|
[44] | 372 | ENDIF ELSE earth = -1 |
---|
| 373 | ; we look for missing_value |
---|
[271] | 374 | ; we apply add_offset, scale_factor and missing_value |
---|
| 375 | IF keyword_set(addscl_before) THEN BEGIN |
---|
[275] | 376 | IF scale_factor NE 1 THEN res = temporary(res)*scale_factor |
---|
| 377 | IF add_offset NE 0 THEN res = temporary(res)+add_offset |
---|
[271] | 378 | ENDIF |
---|
[275] | 379 | IF size(missing_value, /type) NE 7 THEN BEGIN |
---|
[271] | 380 | IF finite(missing_value) EQ 1 THEN BEGIN |
---|
[336] | 381 | CASE 1 OF |
---|
| 382 | missing_value GT 1.e6:missing = where(res GT missing_value/10.) |
---|
| 383 | missing_value LT -1.e6:missing = where(res LT missing_value/10.) |
---|
| 384 | abs(missing_value) LT 1.e-6:missing = where(abs(res) LT 1.e-6) |
---|
| 385 | ELSE:missing = where(res EQ missing_value) |
---|
| 386 | ENDCASE |
---|
[271] | 387 | ENDIF ELSE missing = where(finite(res) EQ 0) |
---|
[44] | 388 | ENDIF ELSE missing = -1 |
---|
[271] | 389 | IF NOT keyword_set(addscl_before) THEN BEGIN |
---|
| 390 | if scale_factor NE 1 then res = temporary(res)*scale_factor |
---|
| 391 | if add_offset NE 0 then res = temporary(res)+add_offset |
---|
| 392 | ENDIF |
---|
[275] | 393 | IF missing[0] NE -1 THEN res[temporary(missing)] = !values.f_nan |
---|
[365] | 394 | IF earth[0] NE -1 THEN BEGIN |
---|
| 395 | IF varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] NE -1 THEN $ ;xyt array |
---|
[275] | 396 | earth = earth#replicate(1, jpt)+replicate(nx*ny, n_elements(earth))#lindgen(jpt) |
---|
[500] | 397 | IF varinq.ndims eq 4 THEN earth = earth#replicate(1L, jpt)+replicate(nx*ny*nz, n_elements(earth))#lindgen(jpt) |
---|
[460] | 398 | CASE size(res, /type) OF |
---|
| 399 | 1:res[temporary(earth)] = 0b |
---|
| 400 | 2:res[temporary(earth)] = 0 |
---|
[466] | 401 | 3:res[temporary(earth)] = 0L |
---|
[460] | 402 | ELSE:res[temporary(earth)] = valmask |
---|
| 403 | ENDCASE |
---|
[365] | 404 | ENDIF |
---|
[2] | 405 | ;--------------------------------------------------------------------- |
---|
[493] | 406 | ; if it is roms outputs, we need to get additional infos... |
---|
[172] | 407 | IF NOT keyword_set(callitself) THEN BEGIN |
---|
[501] | 408 | IF (strmid(dimnames[0], 0, 3) EQ 'xi_' AND strmid(dimnames[1], 0, 4) EQ 'eta_') OR keyword_set(isroms) THEN BEGIN |
---|
[172] | 409 | ncdf_attget, cdfid, 'theta_s', theta_s, /global |
---|
| 410 | ncdf_attget, cdfid, 'theta_b', theta_b, /global |
---|
| 411 | ncdf_attget, cdfid, 'hc', hc, /global |
---|
[174] | 412 | ; look for all variables names |
---|
[271] | 413 | allvarnames = strarr(inq.nvars) |
---|
| 414 | FOR i = 0, inq.nvars-1 DO BEGIN |
---|
[174] | 415 | tmp = ncdf_varinq( cdfid, i) |
---|
| 416 | allvarnames[i] = tmp.name |
---|
| 417 | ENDFOR |
---|
| 418 | CASE 1 OF |
---|
| 419 | keyword_set(zetazero):zeta = fltarr(nx, ny, jpt) |
---|
| 420 | keyword_set(zetafilename): $ |
---|
| 421 | zeta = read_ncdf('zeta', firsttps, lasttps, FILENAME = zetafilename $ |
---|
| 422 | , /TIMESTEP, /NOSTRUCT, CONT_NOFILL = CONT_NOFILL $ |
---|
| 423 | , GRID = vargrid, /CALLITSELF, _EXTRA = ex) |
---|
[181] | 424 | (where(allvarnames EQ 'zeta'))[0] NE -1: $ |
---|
[174] | 425 | zeta = read_ncdf('zeta', firsttps, lasttps, FILENAME = filename $ |
---|
| 426 | , /TIMESTEP, /NOSTRUCT, CONT_NOFILL = CONT_NOFILL $ |
---|
| 427 | , GRID = vargrid, /CALLITSELF, _EXTRA = ex) |
---|
| 428 | ELSE:return, report('The variable zeta was not found in the file, please use the keyword ZETAFILENAME to specify the name of a file containing zeta or use keyword ZETAZERO to define zeta to 0.') |
---|
| 429 | ENDCASE |
---|
[192] | 430 | romszinfos = {h:romszinfos.h, zeta:temporary(zeta), theta_s:theta_s, theta_b:theta_b, hc:hc} |
---|
[172] | 431 | ENDIF ELSE romszinfos = {h:-1, zeta:-1, theta_s:-1, theta_b:-1, hc:-1} |
---|
| 432 | ENDIF |
---|
| 433 | ;--------------------------------------------------------------------- |
---|
[276] | 434 | IF keyword_set(direc) THEN BEGIN |
---|
[310] | 435 | IF jpt EQ 1 THEN res = moyenne(temporary(res), direc, _extra = ex) $ |
---|
[276] | 436 | ELSE BEGIN |
---|
[310] | 437 | res = grossemoyenne(temporary(res), direc, _extra = ex) |
---|
[276] | 438 | IF ( strpos(strlowcase(direc), 't') ge 0 ) THEN BEGIN |
---|
| 439 | vardate = strtrim(jul2date(time[0]), 1)+' - '+strtrim(jul2date(time[jpt-1]), 1) |
---|
| 440 | time = total(time)/float(jpt) |
---|
| 441 | jpt = 1 |
---|
| 442 | ENDIF |
---|
[365] | 443 | ENDELSE |
---|
[276] | 444 | ENDIF |
---|
| 445 | ;--------------------------------------------------------------------- |
---|
[44] | 446 | ncdf_close, cdfid |
---|
[172] | 447 | IF keyword_set(nostruct) THEN return, res |
---|
| 448 | IF keyword_set(key_forgetold) THEN BEGIN |
---|
[226] | 449 | return, {arr:temporary(res), grid:vargrid, unit:varunit, experiment:varexp, name:varname} |
---|
| 450 | ENDIF ELSE BEGIN |
---|
[172] | 451 | return, {tab:temporary(res), grille:vargrid, unite:varunit, experience:varexp, nom:varname} |
---|
| 452 | ENDELSE |
---|
| 453 | |
---|
[44] | 454 | END |
---|