[2] | 1 | ;+ |
---|
[150] | 2 | ; @file_comments |
---|
| 3 | ; Reading function for the file net_cdf. |
---|
[226] | 4 | ; This program is less universal than ncdf_lec (it appeal to declared |
---|
[150] | 5 | ; variables in common.pro) but it is very easier to be used. It considerate |
---|
| 6 | ; the declaration of the different zooms which have been defined |
---|
| 7 | ; (ixminmesh...premierx...), the declaration of the variable key_shift... |
---|
| 8 | ; To put it in a nutshell, the result of read_ncdf can be directly used in plt... |
---|
[226] | 9 | ; This is also this program which is used by default in our reading widgets. |
---|
[2] | 10 | ; |
---|
[150] | 11 | ; @categories |
---|
[157] | 12 | ; Reading |
---|
[226] | 13 | ; |
---|
[163] | 14 | ; @param NAME {in}{required}{type=string} |
---|
| 15 | ; It define the field to be read. |
---|
[2] | 16 | ; |
---|
[150] | 17 | ; @param BEGINNING {in}{required} |
---|
| 18 | ; Relative with the time axis. |
---|
| 19 | ; These can be |
---|
[226] | 20 | ; - 2 date of the type yyyymmdd and in this case, we select dates |
---|
[150] | 21 | ; which are included between these two dates. |
---|
[226] | 22 | ; - 2 indexes which define between which and which time step we have |
---|
[163] | 23 | ; to extract the temporal dimension. |
---|
[2] | 24 | ; |
---|
[150] | 25 | ; @param ENDING {in}{required} |
---|
| 26 | ; Relative with the time axis. |
---|
| 27 | ; See BEGINNING. |
---|
[226] | 28 | ; |
---|
[174] | 29 | ; @param COMPATIBILITY {in}{optional} |
---|
| 30 | ; Useless, defined for compatibility |
---|
[226] | 31 | ; |
---|
| 32 | ; @keyword BOXZOOM |
---|
| 33 | ; Contain the boxzoom on which we have to do the reading |
---|
| 34 | ; |
---|
[174] | 35 | ; @keyword CALLITSELF {default=0}{type=scalar: 0 or 1} |
---|
| 36 | ; For ROMS outputs. Use by read_ncdf itself to access auxilliary data (h and zeta). |
---|
[226] | 37 | ; |
---|
[174] | 38 | ; @keyword FILENAME {required}{type=string} |
---|
[163] | 39 | ; It contains he file's name. |
---|
[226] | 40 | ; |
---|
[174] | 41 | ; @keyword INIT {default=0}{type=scalar: 0 or 1} |
---|
[150] | 42 | ; To call automatically initncdf, filename and thus |
---|
| 43 | ; redefine all the grid parameters |
---|
[226] | 44 | ; |
---|
[150] | 45 | ; @keyword GRID |
---|
[163] | 46 | ; ='[UTVWF]' to specify the type of grid. Default is (1) |
---|
[150] | 47 | ; based on the name of the file if the file ends by |
---|
| 48 | ; GRID[._][TUVFW].NC (not case sensible) or (2) T if case (1) |
---|
| 49 | ; is not found. |
---|
[226] | 50 | ; |
---|
[174] | 51 | ; @keyword TIMESTEP {default=0}{type=scalar: 0 or 1} |
---|
| 52 | ; Specify that BEGINNING and ENDING refer to indexes of the time axis and not to dates |
---|
[2] | 53 | ; |
---|
[174] | 54 | ; @keyword TOUT {default=0}{type=scalar: 0 or 1} |
---|
[226] | 55 | ; We activate it if we want to read the file on the whole domain without |
---|
| 56 | ; considerate the sub-domain defined by the boxzoom or |
---|
[150] | 57 | ; lon1,lon2,lat1,lat2,vert1,vert2. |
---|
[226] | 58 | ; |
---|
[174] | 59 | ; @keyword NOSTRUCT {default=0}{type=scalar: 0 or 1} |
---|
[226] | 60 | ; We activate it if we do not want that read_ncdf send back a structure |
---|
[163] | 61 | ; but only the array referring to the field. |
---|
[226] | 62 | ; |
---|
[163] | 63 | ; @keyword TIMEVAR {type=string} |
---|
| 64 | ; It define the name of the variable that |
---|
| 65 | ; contains the time axis. This keyword can be useful if there |
---|
| 66 | ; is no unlimited dimension or if the time axis selected by default |
---|
[150] | 67 | ; (the first 1D array with unlimited dimension) is not the good one. |
---|
[2] | 68 | ; |
---|
[174] | 69 | ; @keyword ZETAFILENAME {default=FILENAME}{type=string} |
---|
| 70 | ; For ROMS outputs. The filename of the file where zeta vriable should be read |
---|
| 71 | ; |
---|
| 72 | ; @keyword ZETAZERO {default=0}{type=scalar: 0 or 1} |
---|
| 73 | ; For ROMS outputs. To define zeta to 0. instead of reading it |
---|
| 74 | ; |
---|
[150] | 75 | ; @keyword _EXTRA |
---|
[226] | 76 | ; Used to pass your keywords |
---|
[2] | 77 | ; |
---|
[150] | 78 | ; @returns |
---|
[163] | 79 | ; Structure readable by litchamp.pro or an array if NOSTRUCT is activated. |
---|
[226] | 80 | ; |
---|
[150] | 81 | ; @uses |
---|
| 82 | ; common.pro |
---|
[226] | 83 | ; |
---|
[150] | 84 | ; @restrictions |
---|
| 85 | ; The field must have a temporal dimension. |
---|
[226] | 86 | ; |
---|
[150] | 87 | ; @history |
---|
[157] | 88 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
[226] | 89 | ; 15/10/1999 |
---|
| 90 | ; |
---|
[150] | 91 | ; @version |
---|
[226] | 92 | ; $Id$ |
---|
[2] | 93 | ;- |
---|
[150] | 94 | ;--------------------------------------------------------- |
---|
| 95 | ;--------------------------------------------------------- |
---|
| 96 | ;--------------------------------------------------------- |
---|
| 97 | |
---|
| 98 | FUNCTION read_ncdf, name, beginning, ending, compatibility, BOXZOOM = boxzoom, FILENAME = filename $ |
---|
[44] | 99 | , PARENTIN = parentin, TIMESTEP = timestep, TIMEVAR = timevar $ |
---|
| 100 | , TOUT = tout, NOSTRUCT = nostruct, CONT_NOFILL = CONT_NOFILL, INIT = init $ |
---|
[209] | 101 | , GRID = grid, CALLITSELF = callitself $ |
---|
[192] | 102 | , ZETAFILENAME = zetafilename, ZETAZERO = zetazero $ |
---|
[174] | 103 | , _EXTRA = ex |
---|
[44] | 104 | ;--------------------------------------------------------- |
---|
[114] | 105 | ; |
---|
| 106 | compile_opt idl2, strictarrsubs |
---|
| 107 | ; |
---|
[44] | 108 | @cm_4mesh |
---|
| 109 | @cm_4data |
---|
| 110 | @cm_4cal |
---|
| 111 | IF NOT keyword_set(key_forgetold) THEN BEGIN |
---|
| 112 | @updatenew |
---|
| 113 | @updatekwd |
---|
| 114 | ENDIF |
---|
[2] | 115 | ;------------------------------------------------------------ |
---|
| 116 | ; we find the filename. |
---|
| 117 | ;------------------------------------------------------------ |
---|
[226] | 118 | ; print,filename |
---|
[2] | 119 | ; is parent a valid widget ? |
---|
[44] | 120 | if keyword_set(parentin) then BEGIN |
---|
| 121 | parent = long(parentin) |
---|
| 122 | parent = parent*widget_info(parent, /managed) |
---|
| 123 | ENDIF |
---|
| 124 | filename = isafile(filename = filename, IODIRECTORY = iodir, _EXTRA = ex) |
---|
[2] | 125 | ;------------------------------------------------------------ |
---|
[150] | 126 | ; Opening of the name file |
---|
[2] | 127 | ;------------------------------------------------------------ |
---|
[44] | 128 | if size(filename, /type) NE 7 then $ |
---|
[2] | 129 | return, report('read_ncdf cancelled') |
---|
[44] | 130 | IF !version.OS_FAMILY EQ 'unix' THEN spawn, '\file '+filename+' > /dev/null' |
---|
| 131 | cdfid = ncdf_open(filename) |
---|
| 132 | contient = ncdf_inquire(cdfid) |
---|
[2] | 133 | ;------------------------------------------------------------ |
---|
| 134 | ; we check if the variable name exists in the file. |
---|
| 135 | ;------------------------------------------------------------ |
---|
[44] | 136 | if ncdf_varid(cdfid, name) EQ -1 then BEGIN |
---|
| 137 | ncdf_close, cdfid |
---|
[2] | 138 | return, report('variable '+name+' !C not found in the file '+filename) |
---|
[44] | 139 | ENDIF |
---|
| 140 | varcontient = ncdf_varinq(cdfid, name) |
---|
[172] | 141 | IF varcontient.ndims LT 2 THEN return, report('read_ncdf cannot read scalar or 1D data') |
---|
| 142 | ; look for the dimension names |
---|
| 143 | dimnames = strarr(varcontient.ndims) |
---|
| 144 | FOR i = 0, varcontient.ndims-1 DO BEGIN |
---|
| 145 | ncdf_diminq, cdfid, varcontient.dim[i], tmp, dimsize |
---|
| 146 | dimnames[i] = tmp |
---|
[226] | 147 | ENDFOR |
---|
[2] | 148 | ;------------------------------------------------------------ |
---|
[44] | 149 | ; shall we redefine the grid parameters |
---|
| 150 | ;------------------------------------------------------------ |
---|
| 151 | if keyword_set(init) THEN initncdf, filename, _extra = ex |
---|
| 152 | ;------------------------------------------------------------ |
---|
[150] | 153 | ; check the time axis and the debut and ending dates |
---|
[2] | 154 | ;------------------------------------------------------------ |
---|
[150] | 155 | if n_elements(beginning) EQ 0 then begin |
---|
| 156 | beginning = 0 |
---|
[44] | 157 | timestep = 1 |
---|
| 158 | endif |
---|
| 159 | if keyword_set(timestep) then begin |
---|
[150] | 160 | firsttps = beginning[0] |
---|
| 161 | if n_elements(ending) NE 0 then lasttps = ending[0] ELSE lasttps = firsttps |
---|
[44] | 162 | jpt = lasttps-firsttps+1 |
---|
[172] | 163 | IF NOT keyword_set(callitself) then time = julday(1, 1, 1) + lindgen(jpt) |
---|
[44] | 164 | ENDIF ELSE BEGIN |
---|
| 165 | if keyword_set(parent) then BEGIN |
---|
| 166 | widget_control, parent, get_uvalue = top_uvalue |
---|
| 167 | filelist = extractatt(top_uvalue, 'filelist') |
---|
| 168 | IF filelist[0] EQ 'many !' THEN filelist = filename |
---|
| 169 | currentfile = (where(filelist EQ filename))[0] |
---|
| 170 | time = (*(extractatt(top_uvalue, 'fileparameters'))[currentfile]).time_counter |
---|
[150] | 171 | date1 = date2jul(beginning[0]) |
---|
| 172 | if n_elements(ending) NE 0 then date2 = date2jul(ending[0]) ELSE date2 = date1 |
---|
[44] | 173 | firsttps = where(time EQ date1) & firsttps = firsttps[0] |
---|
| 174 | lasttps = where(time EQ date2) & lasttps = lasttps[0] |
---|
| 175 | ENDIF ELSE BEGIN |
---|
[226] | 176 | IF keyword_set(timevar) THEN BEGIN |
---|
[44] | 177 | timeid = ncdf_varid(cdfid, timevar) |
---|
| 178 | IF timeid EQ -1 THEN BEGIN |
---|
| 179 | ncdf_close, cdfid |
---|
| 180 | return, report('the file '+filename+' as no variable ' + timevar $ |
---|
| 181 | + '. !C Use the TIMESTEP keyword') |
---|
| 182 | endif |
---|
| 183 | timecontient = ncdf_varinq(cdfid, timeid) |
---|
| 184 | contient.recdim = timecontient.dim[0] |
---|
[226] | 185 | ENDIF ELSE BEGIN |
---|
[2] | 186 | ; we find the infinite dimension |
---|
[44] | 187 | timedim = contient.recdim |
---|
| 188 | if timedim EQ -1 then BEGIN |
---|
| 189 | ncdf_close, cdfid |
---|
| 190 | return, report('the file '+filename+' as no infinite dimension. !C Use TIMESTEP or TIMEVAR keyword') |
---|
| 191 | endif |
---|
[226] | 192 | ; we find the FIRST time axis |
---|
[44] | 193 | timeid = 0 |
---|
[150] | 194 | repeat BEGIN ; As long as we have not find a variable having only one dimension: the infinite one |
---|
| 195 | timecontient = ncdf_varinq(cdfid, timeid) ; that the variable contain. |
---|
[44] | 196 | timeid = timeid+1 |
---|
| 197 | endrep until (n_elements(timecontient.dim) EQ 1 $ |
---|
| 198 | AND timecontient.dim[0] EQ contient.recdim) $ |
---|
[2] | 199 | OR timeid EQ contient.nvars+1 |
---|
| 200 | ; |
---|
[44] | 201 | if timeid EQ contient.nvars+1 then BEGIN |
---|
| 202 | ncdf_close, cdfid |
---|
[2] | 203 | return, report('the file '+filename+' as no time axis variable. !C Use the TIMESTEP keyword') |
---|
[44] | 204 | endif |
---|
| 205 | timeid = timeid-1 |
---|
[226] | 206 | ENDELSE |
---|
[2] | 207 | ; we must found the time origin of the julian calendar used in the |
---|
[226] | 208 | ; time axis. |
---|
[44] | 209 | ; does the attribut units an dcalendar exist for the variable time axis? |
---|
| 210 | if timecontient.natts EQ 0 then BEGIN |
---|
| 211 | ncdf_close, cdfid |
---|
| 212 | return, report('the variable '+timecontient.name+' has no attribut.!C Use the TIMESTEP keyword or add the attribut ''units'' to the variable') |
---|
| 213 | endif |
---|
| 214 | attnames = strarr(timecontient.natts) |
---|
| 215 | for attiq = 0, timecontient.natts-1 do attnames[attiq] = ncdf_attname(cdfid, timeid, attiq) |
---|
| 216 | if (where(attnames EQ 'units'))[0] EQ -1 then BEGIN |
---|
| 217 | ncdf_close, cdfid |
---|
| 218 | return, report('Attribut ''units'' not found for the variable '+timecontient.name+' !C Use the TIMESTEP keyword') |
---|
| 219 | ENDIF |
---|
[2] | 220 | ; |
---|
[44] | 221 | ; now we try to find the attribut called calendar... |
---|
[150] | 222 | ; the attribute "calendar" exists? |
---|
[44] | 223 | ; If no, we suppose that the calendar is gregorian calendar |
---|
| 224 | ; |
---|
| 225 | if (where(attnames EQ 'calendar'))[0] NE -1 then BEGIN |
---|
| 226 | ncdf_attget, cdfid, timeid, 'calendar', value |
---|
| 227 | value = string(value) |
---|
| 228 | CASE value OF |
---|
[69] | 229 | 'noleap':key_caltype = 'noleap' |
---|
[44] | 230 | '360d':key_caltype = '360d' |
---|
| 231 | 'greg':IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg' |
---|
| 232 | ELSE:BEGIN |
---|
[226] | 233 | ; notused = report('Unknown calendar: '+value+', we use greg calendar.') |
---|
[44] | 234 | key_caltype = 'greg' |
---|
| 235 | END |
---|
| 236 | ENDCASE |
---|
| 237 | ENDIF ELSE BEGIN |
---|
[226] | 238 | ; notused = report('Unknown calendar, we use '+key_caltype+' calendar.') |
---|
[44] | 239 | IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg' |
---|
| 240 | ENDELSE |
---|
| 241 | ; |
---|
| 242 | ; now we take acre of units attribut |
---|
| 243 | ncdf_attget, cdfid, timeid, 'units', value |
---|
| 244 | ; |
---|
[2] | 245 | ; time_counter:units = "seconds since 0001-01-01 00:00:00" ; |
---|
| 246 | ; time_counter:units = "hours since 0001-01-01 00:00:00" ; |
---|
| 247 | ; time_counter:units = "days since 1979-01-01 00:00:00" ; |
---|
| 248 | ; time_counter:units = "months since 1979-01-01 00:00:00" ; |
---|
| 249 | ; time_counter:units = "years since 1979-01-01 00:00:00" ; |
---|
| 250 | ; |
---|
| 251 | ; we decript the "units" attribut to find the time origin |
---|
[44] | 252 | value = strtrim(strcompress(string(value)), 2) |
---|
| 253 | mots = str_sep(value, ' ') |
---|
| 254 | unite = mots[0] |
---|
[198] | 255 | unite = strlowcase(unite) |
---|
| 256 | IF strpos(unite, 's', strlen(unite)-1) NE -1 THEN unite = strmid(unite, 0, strlen(unite)-1) |
---|
| 257 | IF strpos(unite, 'julian_') NE -1 THEN unite = strmid(unite, 7) |
---|
| 258 | IF unite NE 'second' AND unite NE 'hour' AND unite NE 'day' $ |
---|
| 259 | AND unite NE 'month' AND unite NE 'year' THEN BEGIN |
---|
[172] | 260 | ncdf_close, cdfid |
---|
| 261 | return, report('time units does not start with seconds/hours/days/months/years') |
---|
| 262 | ENDIF |
---|
[199] | 263 | IF stregex(value, '[^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*', /boolean) EQ 0 THEN BEGIN |
---|
[172] | 264 | ncdf_close, cdfid |
---|
[199] | 265 | return, report('attribut units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*') |
---|
[172] | 266 | ENDIF |
---|
[44] | 267 | depart = str_sep(mots[2], '-') |
---|
| 268 | ncdf_varget, cdfid, timeid, time |
---|
| 269 | time = double(time) |
---|
| 270 | case unite of |
---|
[69] | 271 | 'second':time = julday(depart[1], depart[2], depart[0])+time/86400.d |
---|
| 272 | 'hour':time = julday(depart[1], depart[2], depart[0])+time/24.d |
---|
[44] | 273 | 'day':time = julday(depart[1], depart[2], depart[0])+time |
---|
[226] | 274 | 'month':BEGIN |
---|
[44] | 275 | if total(fix(time) NE time) NE 0 then $ ; we switch to days with 30d/m |
---|
| 276 | time = julday(depart[1], depart[2], depart[0])+round(time*30) $ |
---|
| 277 | ELSE for t = 0, n_elements(time)-1 DO $ |
---|
| 278 | time[t] = julday(depart[1]+time[t], depart[2], depart[0]) |
---|
| 279 | END |
---|
| 280 | 'year':BEGIN |
---|
| 281 | if total(fix(time) NE time) NE 0 then $ ; we switch to days with 365d/y |
---|
| 282 | time = julday(depart[1], depart[2], depart[0])+round(time*365) $ |
---|
| 283 | ELSE for t = 0, n_elements(time)-1 do $ |
---|
| 284 | time[t] = julday(depart[1], depart[2], depart[0]+time[t]) |
---|
| 285 | END |
---|
| 286 | ELSE:BEGIN |
---|
| 287 | ncdf_close, cdfid |
---|
[172] | 288 | return, report('The "units" attribute of the time axis must be something like: !C "seconds since 0001-01-01 ..." !C "days since 1979-01-01 ..." !C "months since 1979-01-01 ..." !C "years since 1979-01-01 ..."') |
---|
[44] | 289 | end |
---|
| 290 | ENDCASE |
---|
[150] | 291 | date1 = date2jul(beginning[0]) |
---|
| 292 | if n_elements(ending) NE 0 then date2 = date2jul(ending[0]) ELSE date2 = date1 |
---|
[69] | 293 | time = double(time) |
---|
[44] | 294 | firsttps = where(time GE date1) & firsttps = firsttps[0] |
---|
| 295 | if firsttps EQ -1 THEN BEGIN |
---|
| 296 | ncdf_close, cdfid |
---|
| 297 | return, report('date 1: '+strtrim(jul2date(date1), 1)+' is not found in the time axis.') |
---|
[74] | 298 | ENDIF |
---|
[44] | 299 | lasttps = where(time LE date2) |
---|
| 300 | if lasttps[0] EQ -1 THEN BEGIN |
---|
| 301 | ncdf_close, cdfid |
---|
[198] | 302 | return, report('the time axis has no date before date 2: '+strtrim(jul2date(date2), 1)) |
---|
[44] | 303 | endif |
---|
| 304 | lasttps = lasttps[n_elements(lasttps)-1] |
---|
| 305 | if lasttps LT firsttps then BEGIN |
---|
| 306 | ncdf_close, cdfid |
---|
[198] | 307 | return, report('the time axis has no dates between date1 and date 2: '+strtrim(jul2date(date1), 1)+' '+strtrim(jul2date(date2), 1)) |
---|
[44] | 308 | endif |
---|
| 309 | ENDELSE |
---|
| 310 | time = time[firsttps:lasttps] |
---|
| 311 | jpt = lasttps-firsttps+1 |
---|
| 312 | ENDELSE |
---|
[2] | 313 | ;------------------------------------------------------------ |
---|
[150] | 314 | ; Name of the grid on which the field refer to. |
---|
[2] | 315 | ;------------------------------------------------------------ |
---|
[44] | 316 | IF keyword_set(grid) THEN vargrid = strupcase(grid) ELSE BEGIN |
---|
| 317 | vargrid = 'T' ; default definition |
---|
[69] | 318 | IF finite(glamu[0]) EQ 1 THEN BEGIN |
---|
[172] | 319 | ; are we in one of the case corresponding to ROMS conventions? |
---|
| 320 | CASE 1 OF |
---|
| 321 | dimnames[2 <(varcontient.ndims-1)] EQ 's_w':vargrid = 'W' |
---|
| 322 | dimnames[0] EQ 'xi_rho' AND dimnames[1] EQ 'eta_rho':vargrid = 'T' |
---|
| 323 | dimnames[0] EQ 'xi_u' AND dimnames[1] EQ 'eta_u' :vargrid = 'U' |
---|
| 324 | dimnames[0] EQ 'xi_v' AND dimnames[1] EQ 'eta_v' :vargrid = 'V' |
---|
| 325 | dimnames[0] EQ 'xi_psi' AND dimnames[1] EQ 'eta_psi':vargrid = 'F' |
---|
| 326 | dimnames[0] EQ 'xi_rho' AND dimnames[1] EQ 'eta_v' :vargrid = 'V' |
---|
| 327 | dimnames[0] EQ 'xi_u' AND dimnames[1] EQ 'eta_rho':vargrid = 'U' |
---|
| 328 | dimnames[0] EQ 'xi_u' AND dimnames[1] EQ 'eta_v' :vargrid = 'F' |
---|
[226] | 329 | ELSE:BEGIN |
---|
[172] | 330 | ; could we define the grid type from the file name?? |
---|
| 331 | pattern = ['GRID.', 'GRID_', 'GRID', 'UPID_', '30ID_'] |
---|
| 332 | gdtype = ['T', 'U', 'V', 'W', 'F'] |
---|
| 333 | fnametest = strupcase(filename) |
---|
| 334 | FOR i = 0, n_elements(pattern)-1 DO BEGIN |
---|
| 335 | FOR j = 0, n_elements(gdtype)-1 DO BEGIN |
---|
| 336 | substr = pattern[i]+gdtype[j] |
---|
| 337 | pos = strpos(fnametest, substr) |
---|
| 338 | IF pos NE -1 THEN $ |
---|
| 339 | vargrid = strmid(fnametest, pos+strlen(substr)-1, 1) |
---|
| 340 | ENDFOR |
---|
| 341 | ENDFOR |
---|
| 342 | END |
---|
| 343 | ENDCASE |
---|
| 344 | ENDIF |
---|
[44] | 345 | ENDELSE |
---|
[2] | 346 | ;--------------------------------------------------------------- |
---|
[44] | 347 | ; call the init function ? |
---|
| 348 | ;--------------------------------------------------------------- |
---|
| 349 | |
---|
| 350 | ;--------------------------------------------------------------- |
---|
[150] | 351 | ; redefinition of the domain |
---|
[2] | 352 | ;--------------------------------------------------------------- |
---|
[44] | 353 | if keyword_set(tout) then begin |
---|
| 354 | nx = jpi |
---|
| 355 | ny = jpj |
---|
| 356 | nz = jpk |
---|
| 357 | firstx = 0 |
---|
| 358 | firsty = 0 |
---|
| 359 | firstz = 0 |
---|
| 360 | lastx = jpi-1 |
---|
| 361 | lasty = jpj-1 |
---|
| 362 | lastz = jpk-1 |
---|
| 363 | case strupcase(vargrid) of |
---|
| 364 | 'T':mask = tmask |
---|
| 365 | 'U':mask = umask() |
---|
| 366 | 'V':mask = vmask() |
---|
| 367 | 'W':mask = tmask |
---|
| 368 | 'F':mask = fmask() |
---|
| 369 | endcase |
---|
| 370 | ENDIF ELSE BEGIN |
---|
[226] | 371 | if keyword_set(boxzoom) then BEGIN |
---|
[44] | 372 | Case 1 Of |
---|
| 373 | N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]] |
---|
| 374 | N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]] |
---|
| 375 | N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2] |
---|
| 376 | N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]] |
---|
| 377 | N_Elements(Boxzoom) Eq 6:bte = Boxzoom |
---|
| 378 | Else: BEGIN |
---|
| 379 | ncdf_close, cdfid |
---|
| 380 | return, report('Wrong Definition of Boxzoom') |
---|
| 381 | end |
---|
| 382 | ENDCASE |
---|
| 383 | savedbox = 1b |
---|
| 384 | saveboxparam, 'boxparam4rdncdf.dat' |
---|
| 385 | domdef, bte, GRIDTYPE = ['T', vargrid], _extra = ex |
---|
| 386 | ENDIF |
---|
| 387 | grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz |
---|
[226] | 388 | undefine, glam & undefine, gphi & ; We liberate some memory! |
---|
[44] | 389 | ENDELSE |
---|
[2] | 390 | ;--------------------------------------------------------------------- |
---|
[150] | 391 | ; We initializate ixmindta, iymindta if needed |
---|
[2] | 392 | ;--------------------------------------------------------------------- |
---|
[44] | 393 | if n_elements(jpidta) EQ 0 THEN jpidta = jpiglo |
---|
| 394 | if n_elements(jpjdta) EQ 0 THEN jpjdta = jpjglo |
---|
| 395 | if n_elements(jpkdta) EQ 0 THEN jpkdta = jpkglo |
---|
| 396 | if n_elements(ixmindta) EQ 0 THEN ixmindta = 0 |
---|
| 397 | if n_elements(ixmaxdta) EQ 0 then ixmaxdta = jpidta-1 |
---|
| 398 | if ixmindta EQ -1 THEN ixmindta = 0 |
---|
| 399 | IF ixmaxdta EQ -1 then ixmaxdta = jpidta-1 |
---|
| 400 | if n_elements(iymindta) EQ 0 THEN iymindta = 0 |
---|
| 401 | IF n_elements(iymaxdta) EQ 0 then iymaxdta = jpjdta-1 |
---|
| 402 | if iymindta EQ -1 THEN iymindta = 0 |
---|
| 403 | IF iymaxdta EQ -1 then iymaxdta = jpjdta-1 |
---|
| 404 | if n_elements(izmindta) EQ 0 THEN izmindta = 0 |
---|
| 405 | IF n_elements(izmaxdta) EQ 0 then izmaxdta = jpkdta-1 |
---|
| 406 | if izmindta EQ -1 THEN izmindta = 0 |
---|
| 407 | IF izmaxdta EQ -1 then izmaxdta = jpkdta-1 |
---|
[2] | 408 | ;---------------------------------------------------------------- |
---|
[150] | 409 | ; We will read the file |
---|
[2] | 410 | ;--------------------------------------------------------------- |
---|
[44] | 411 | if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1] |
---|
| 412 | key_stride = 1l > long(key_stride) |
---|
| 413 | ;--------------------------------------------------------------------- |
---|
| 414 | ;--------------------------------------------------------------------- |
---|
[2] | 415 | @read_ncdf_varget |
---|
| 416 | ;--------------------------------------------------------------------- |
---|
[44] | 417 | ;--------------------------------------------------------------------- |
---|
[150] | 418 | ; We define global variable joined with the variable. |
---|
[2] | 419 | ;--------------------------------------------------------------------- |
---|
| 420 | ; varname |
---|
[172] | 421 | IF NOT keyword_set(callitself) THEN varname = name |
---|
[2] | 422 | ; varunit |
---|
[44] | 423 | if varcontient.natts NE 0 then begin |
---|
| 424 | attnames = strarr(varcontient.natts) |
---|
| 425 | for attiq = 0, varcontient.natts-1 do attnames[attiq] = ncdf_attname(cdfid, name, attiq) |
---|
| 426 | lowattnames = strlowcase(attnames) |
---|
[2] | 427 | ; |
---|
[44] | 428 | found = (where(lowattnames EQ 'units'))[0] |
---|
| 429 | IF found NE -1 then ncdf_attget, cdfid, name, attnames[found], value ELSE value = '' |
---|
[172] | 430 | IF NOT keyword_set(callitself) THEN varunit = strtrim(string(value), 2) |
---|
[2] | 431 | ; |
---|
[44] | 432 | found = (where(lowattnames EQ 'add_offset'))[0] |
---|
| 433 | if found NE -1 then ncdf_attget, cdfid, name, attnames[found], add_offset ELSE add_offset = 0. |
---|
[2] | 434 | ; |
---|
[44] | 435 | found = (where(lowattnames EQ 'scale_factor'))[0] |
---|
| 436 | if found NE -1 then ncdf_attget, cdfid, name, attnames[found], scale_factor ELSE scale_factor = 1. |
---|
[2] | 437 | ; |
---|
[44] | 438 | missing_value = 'no' |
---|
| 439 | found = (where(lowattnames EQ '_fillvalue'))[0] |
---|
| 440 | if found NE -1 then ncdf_attget, cdfid, name, attnames[found], missing_value |
---|
| 441 | found = (where(lowattnames EQ 'missing_value'))[0] |
---|
| 442 | if found NE -1 then ncdf_attget, cdfid, name, attnames[found], missing_value |
---|
[2] | 443 | ; |
---|
[226] | 444 | ENDIF ELSE BEGIN |
---|
[172] | 445 | IF NOT keyword_set(callitself) THEN varunit = '' |
---|
[44] | 446 | add_offset = 0. |
---|
| 447 | scale_factor = 1. |
---|
| 448 | missing_value = 'no' |
---|
| 449 | ENDELSE |
---|
[2] | 450 | ; vardate |
---|
[150] | 451 | ; We make a legible date in function of the specified language. |
---|
| 452 | year = long(beginning[0])/10000 |
---|
| 453 | month = (long(beginning[0])/100) MOD 100 |
---|
| 454 | day = (long(beginning[0]) MOD 100) |
---|
[44] | 455 | vardate = string(format = '(C(CMoA))', 31*(month-1))+' '+strtrim(day, 1)+', '+strtrim(year, 1) |
---|
[74] | 456 | varexp = file_basename(filename) |
---|
[44] | 457 | |
---|
[150] | 458 | ; We apply the value valmask on land points. |
---|
[44] | 459 | if NOT keyword_set(cont_nofill) then begin |
---|
| 460 | valmask = 1e20 |
---|
| 461 | case 1 of |
---|
| 462 | varcontient.ndims eq 2:BEGIN ;xy array |
---|
| 463 | mask = mask[*, *, 0] |
---|
| 464 | earth = where(mask EQ 0) |
---|
| 465 | END |
---|
| 466 | varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] EQ -1:BEGIN ;xyz array |
---|
| 467 | earth = where(mask EQ 0) |
---|
| 468 | END |
---|
| 469 | varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] NE -1:BEGIN ;xyt array |
---|
| 470 | mask = mask[*, *, 0] |
---|
| 471 | earth = where(mask EQ 0) |
---|
| 472 | if earth[0] NE -1 then BEGIN |
---|
| 473 | earth = earth#replicate(1, jpt)+replicate(nx*ny, n_elements(earth))#lindgen(jpt) |
---|
| 474 | END |
---|
| 475 | END |
---|
| 476 | varcontient.ndims eq 4:BEGIN ;xyzt array |
---|
| 477 | earth = where(mask EQ 0) |
---|
| 478 | if earth[0] NE -1 then BEGIN |
---|
| 479 | earth = earth#replicate(1, jpt)+replicate(nx*ny*nz, n_elements(earth))#lindgen(jpt) |
---|
| 480 | END |
---|
| 481 | END |
---|
| 482 | endcase |
---|
| 483 | ENDIF ELSE earth = -1 |
---|
| 484 | ; we look for missing_value |
---|
| 485 | IF size(missing_value, /type) NE 7 then BEGIN |
---|
[226] | 486 | IF size(missing_value, /type) EQ 1 THEN BEGIN |
---|
[206] | 487 | missing_value = strlowcase(string(missing_value)) |
---|
[226] | 488 | IF strmid(missing_value, 0, 1, /reverse_offset) EQ 'f' THEN $ |
---|
[206] | 489 | missing_value = strmid(missing_value, 0, strlen(missing_value)-1) |
---|
[226] | 490 | IF isnumber(string(missing_value), tmp) EQ 1 THEN missing_value = tmp ELSE BEGIN |
---|
[206] | 491 | print, 'Warning: missing value is not a number: ', missing_value |
---|
| 492 | missing_value = - 1 |
---|
| 493 | ENDELSE |
---|
[226] | 494 | ENDIF |
---|
[44] | 495 | ; if missing_value NE valmask then begin |
---|
| 496 | if abs(missing_value) LT 1e6 then missing = where(res EQ missing_value) $ |
---|
| 497 | ELSE missing = where(abs(res) gt abs(missing_value)/10.) |
---|
| 498 | ; ENDIF ELSE missing = -1 |
---|
| 499 | ENDIF ELSE missing = -1 |
---|
[226] | 500 | ; we apply add_offset, scale_factor and missing_value |
---|
[44] | 501 | if scale_factor NE 1 then res = temporary(res)*scale_factor |
---|
| 502 | if add_offset NE 0 then res = temporary(res)+add_offset |
---|
| 503 | if missing[0] NE -1 then res[temporary(missing)] = !values.f_nan |
---|
[172] | 504 | if earth[0] NE -1 then res[temporary(earth)] = 1.e20 |
---|
[2] | 505 | ;--------------------------------------------------------------------- |
---|
[172] | 506 | ; if it is roms outputs, we need to get additionals infos... |
---|
| 507 | IF NOT keyword_set(callitself) THEN BEGIN |
---|
| 508 | IF strmid(dimnames[0], 0, 3) EQ 'xi_' AND strmid(dimnames[1], 0, 4) EQ 'eta_' THEN BEGIN |
---|
| 509 | ncdf_attget, cdfid, 'theta_s', theta_s, /global |
---|
| 510 | ncdf_attget, cdfid, 'theta_b', theta_b, /global |
---|
| 511 | ncdf_attget, cdfid, 'hc', hc, /global |
---|
[174] | 512 | ; look for all variables names |
---|
[181] | 513 | allvarnames = strarr(contient.nvars) |
---|
| 514 | FOR i = 0, contient.nvars-1 DO BEGIN |
---|
[174] | 515 | tmp = ncdf_varinq( cdfid, i) |
---|
| 516 | allvarnames[i] = tmp.name |
---|
| 517 | ENDFOR |
---|
| 518 | CASE 1 OF |
---|
| 519 | keyword_set(zetazero):zeta = fltarr(nx, ny, jpt) |
---|
| 520 | keyword_set(zetafilename): $ |
---|
| 521 | zeta = read_ncdf('zeta', firsttps, lasttps, FILENAME = zetafilename $ |
---|
| 522 | , /TIMESTEP, /NOSTRUCT, CONT_NOFILL = CONT_NOFILL $ |
---|
| 523 | , GRID = vargrid, /CALLITSELF, _EXTRA = ex) |
---|
[181] | 524 | (where(allvarnames EQ 'zeta'))[0] NE -1: $ |
---|
[174] | 525 | zeta = read_ncdf('zeta', firsttps, lasttps, FILENAME = filename $ |
---|
| 526 | , /TIMESTEP, /NOSTRUCT, CONT_NOFILL = CONT_NOFILL $ |
---|
| 527 | , GRID = vargrid, /CALLITSELF, _EXTRA = ex) |
---|
| 528 | 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.') |
---|
| 529 | ENDCASE |
---|
[192] | 530 | romszinfos = {h:romszinfos.h, zeta:temporary(zeta), theta_s:theta_s, theta_b:theta_b, hc:hc} |
---|
[172] | 531 | ENDIF ELSE romszinfos = {h:-1, zeta:-1, theta_s:-1, theta_b:-1, hc:-1} |
---|
| 532 | ENDIF |
---|
| 533 | ;--------------------------------------------------------------------- |
---|
[44] | 534 | ncdf_close, cdfid |
---|
[172] | 535 | |
---|
| 536 | IF keyword_set(savedbox) THEN restoreboxparam, 'boxparam4rdncdf.dat' |
---|
| 537 | |
---|
| 538 | IF keyword_set(nostruct) THEN return, res |
---|
| 539 | IF keyword_set(key_forgetold) THEN BEGIN |
---|
[226] | 540 | return, {arr:temporary(res), grid:vargrid, unit:varunit, experiment:varexp, name:varname} |
---|
| 541 | ENDIF ELSE BEGIN |
---|
[172] | 542 | return, {tab:temporary(res), grille:vargrid, unite:varunit, experience:varexp, nom:varname} |
---|
| 543 | ENDELSE |
---|
| 544 | |
---|
[44] | 545 | END |
---|
[2] | 546 | |
---|
| 547 | |
---|
| 548 | |
---|
| 549 | |
---|
| 550 | |
---|
| 551 | |
---|