- Timestamp:
- 04/28/06 14:18:03 (18 years ago)
- Location:
- trunk
- Files:
-
- 10 added
- 1 deleted
- 16 copied
Legend:
- Unmodified
- Added
- Removed
-
trunk/Grid/computegrid.pro
r12 r13 5 5 ; NAME:computegrid 6 6 ; 7 ; PURPOSE:calcule les parametres (necessaires aux traces) d''une 8 ; grille reguliere avec 9 ; - les points T, U et V sont confondus 10 ; - un seul niveau vertical. 11 ; 12 ; CATEGORY:grille 13 ; 14 ; CALLING SEQUENCE:computegrid, departx, departy, pasx, pasy, nx, ny 7 ; PURPOSE:compute the grid parameters from cm_4mesh common: 8 ; 9 ; horizontal parameters: 10 ; glam[tf], gphi[tf], e1t and e2t 11 ; and if FULLCGRID keyword is defined: 12 ; glam[uv], gphi[uv], e1[uvf] and e2[uvf] 13 ; 14 ; verticals parameters: 15 ; gdep[tw], e3[tw] 16 ; 17 ; masks: tmask 18 ; and if FULLCGRID keyword is defined:[uv]maskred fmaskred[xy] 19 ; 20 ; triangulation: triangles_list 21 ; 22 ; key_ parameters: 23 ; key_shift, key_periodic, key_zreverse, key_yreverse, 24 ; key_stride, key_onearth, key_partialstep 25 ; 26 ; CATEGORY:grid 27 ; 28 ; CALLING SEQUENCE: 29 ; 30 ; computegrid, startx, starty, stepx, stepy, nx, ny 31 ; computegrid, startx, starty, stepx, stepy 32 ; computegrid, xaxis = xaxis, yaxis = yaxis 33 ; or a suitable mix... 15 34 ; 16 ; INPUTS:des saclaires. 17 ; departx:la longitude minimum 18 ; departy:la latitude minimum 19 ; pasx:la pas suivant la longitude 20 ; pasy:la pas suivant la latitude 21 ; nx:le nombre de points en longitude 22 ; ny:le nombre de points en latitude 35 ; INPUTS: 36 ; startx:scalar, x starting point 37 ; starty:scalar, y starting point 38 ; stepx:scalar or vector: x direction step, must be > 0 39 ; if vector nx is not used 40 ; stepy:scalar or vector: y direction step, 41 ; could be > 0 (south to north) or < 0 (north to south) 42 ; if vector ny is not used 43 ; nx:scalar, number of points in x direction 44 ; ny:scalar, number of points in y direction 23 45 ; 24 46 ; KEYWORD PARAMETERS: 25 47 ; 48 ; /FULLCGRID: activate to specify that you want to compute 49 ; all the paremeters of a C grid. Computation of glam[uv], 50 ; gphi[uv], e1[uvf], e2[uvf], [uv]maskred and fmaskred[xy] 51 ; will be add to the default computations 52 ; 53 ; GLAMBOUNDARY: a 2 elements vector, [lon1,lon2], the longitute 54 ; boundaries that should be used to visualize the data. 55 ; we must have lon2 > lon1 and lon2 - lon1 le 360 56 ; key_shift will be defined automaticaly computed according to 57 ; glamboundary by using the FIRST LINE of glamt but 58 ; key_shift will /= 0 only if key_periodic = 1 59 ; 60 ; MASK: to specify the mask with a 2 or 3 dimension array 61 ; 62 ; ONEARTH = 0 or 1: to force the manual definition of 63 ; key_onearth (to specify if the data are on earth -> use longitude 64 ; /latitude etc...). By default, key_onearth = 1. 65 ; note that ONEARTH = 0 forces PERIODIC = 0, SHIFT = 0, 66 ; and is cancelling GLAMBOUNDARY 67 ; 68 ; PERIODIC = 0 or 1: to force the manual definition of 69 ; key_periodic. By default, key_periodic is automaticaly 70 ; computed by using the first line of glamt. 71 ; 72 ; SHIFT = scalar to force the manual definition of key_shift. By 73 ; debault, key_shift is automaticaly computed according to 74 ; glamboundary (when defined) by using the FIRST LINE of glamt. if 75 ; key_periodic=0 then in any case key_shift = 0. 76 ; 77 ; STRIDE = : a 3 elements vector to specify the stride in x, y, z 78 ; direction. Default definition is [1, 1, 1]. The resulting value 79 ; will be stored in the common (cm_4mesh) variable key_stride 80 ; 81 ; XAXIS: to specify longitude1 with a 1 or 2 dimension array, in 82 ; this case startx, stepx and nx are not used but could be 83 ; necessary if the y axis is not defined with yaxis. It must be 84 ; possible to sort the first line of xaxis in the increasing 85 ; order by shifting its elements. 86 ; 87 ; YAXIS: to specify latitudes with a 1 or 2 dimension array, in 88 ; this case starty, stepy and ny are not used but starty and 89 ; stepy could be necessary if the x axis is not defined with xaxis. 90 ; It must be sorted in the increasing or deceasing order 91 ; (along each column if 2d array). 92 ; 93 ; /XYINDEX: activate to specify that the horizontal grid should 94 ; be simply defined by using the index of the points 95 ; (xaxis = findgen(nx) and yaxis = findgen(ny)) 96 ; using this keyword forces key_onearth=0 97 ; 98 ; [XYZ]MINMESH: to define the common variables i[xyz]minmesh 99 ; used to define the grid only in a zoomed part of the original 100 ; grid. Defaut values are 0L, max value is [XYZ]MAXMESH 101 ; 102 ; [XYZ]MAXMESH: to define the common variables i[xyz]maxmesh 103 ; used to define the grid only in a zoomed part of the original 104 ; grid. Defaut values are jp[ijk]glo-1, max value is 105 ; jp[ijk]glo-1. if [XYZ]MAXMESH is negative, then we define 106 ; i[xyz]maxmesh as jp[ijk]glo - 1 + [XYZ]MAXMESH instead of 107 ; [XYZ]MAXMESH 108 ; 109 ; ZAXIS: to specify the vertical axis with a 1 dimension 110 ; array. Must be sorted in the increasing or deceasing order 111 ; 26 112 ; OUTPUTS: 27 113 ; 28 ; COMMON BLOCKS:common.pro 29 ; 30 ; SIDE EFFECTS:definit aussi les parametres pour les grilles U et V 31 ; meme si on n''en a pas besoin. 32 ; 33 ; RESTRICTIONS:programme vite fait, avec les latitudes qui doivent 34 ; aller du sud vers le nord. 114 ; COMMON BLOCKS: cm_4mesh cm_4data cm_4cal 115 ; 116 ; SIDE EFFECTS: 117 ; 118 ; RESTRICTIONS:FUV points definition... 35 119 ; 36 120 ; EXAMPLE: … … 38 122 ; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr) 39 123 ; 2000-04-20 124 ; Sept 2004, several bug fixs to suit C grid type... 125 ; Aug 2005, rewritte almost everything... 40 126 ;- 41 127 ;------------------------------------------------------------ 42 128 ;------------------------------------------------------------ 43 129 ;------------------------------------------------------------ 44 PRO computegrid, departx, departy, pasx, pasy, nx, ny 45 @common 130 PRO computegrid, startx, starty, stepxin, stepyin, nxin, nyin $ 131 , XAXIS = xaxis, YAXIS = yaxis, ZAXIS = zaxis $ 132 , MASK = mask, GLAMBOUNDARY = glamboundary $ 133 , XMINMESH = xminmesh, XMAXMESH = xmaxmesh $ 134 , YMINMESH = yminmesh, YMAXMESH = ymaxmesh $ 135 , ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $ 136 , ONEARTH = onearth, PERIODIC = periodic $ 137 , SHIFT = shift, STRIDE = stride $ 138 , FULLCGRID = fullcgrid, XYINDEX = xyindex $ 139 , FBASE2TBASE = fbase2tbase, _extra = ex 140 ;--------------------------------------------------------- 141 @cm_4mesh 142 @cm_4data 143 @cm_4cal 144 IF NOT keyword_set(key_forgetold) THEN BEGIN 145 @updatenew 146 @updatekwd 147 ENDIF 148 ;--------------------------------------------------------- 46 149 ;------------------------------------------------------------ 47 jpiglo = long(nx) 48 jpjglo = long(ny) 49 jpkglo = 1l 50 ; 51 jpi = jpiglo 52 jpj = jpjglo 53 jpk = jpkglo 54 ; 55 ixminmesh =0l 56 ixmaxmesh =long(jpi-1) 57 iyminmesh =0l 58 iymaxmesh =long(jpj-1) 59 izminmesh =0l 60 izmaxmesh =long(jpk-1) 61 ; 62 glamt = departx+findgen(jpi)*pasx 63 glamt = glamt#replicate(1, jpj) 64 glamf = glamt+pasx/2. 65 gphit = departy+findgen(jpj)*pasy 66 gphit = replicate(1, jpi)#gphit 67 gphif = gphit+pasy/2. 68 glamu = glamf 69 glamv = glamt 70 gphiu = gphit 71 gphiv = gphif 72 gdept = 5. 73 gdepw = 5. 74 ; 75 pasx = abs(pasx) 76 pasy = abs(pasy) 77 ; facteurs d'echelle: 78 r = 6371000. 79 rprojete = r*cos(gphit[0,*]*!pi/180.) 80 e1t = rprojete*pasx*!pi/180. 81 e1t = reform(replicate(1, jpi)#e1t[*], jpi, jpj) 82 e2t = replicate(r*pasy*!pi/180., jpi, jpj) 83 ; 84 e1u = e1t 85 e2u = e2t 86 ; 87 rprojete = r*cos(gphiv[0,*]*!pi/180.) 88 e1v = rprojete*pasx*!pi/180. 89 e1v = reform(replicate(1, jpi)#e1v[*], jpi, jpj) 90 e2v = e2t 91 ; 92 e1f = e1v 93 e2f = e2t 94 ; 95 e3t = 10. 96 e3w = 10. 97 ; mask par defauyt a 1 98 tmask = replicate(1, jpi, jpj) 99 umaskred = replicate(1, jpj) 100 vmaskred = replicate(1, jpi) 101 fmaskredy = replicate(1, jpj) 102 fmaskredx = replicate(1, jpi) 150 time1 = systime(1) ; for key_performance 151 ;------------------------------------------------------------ 152 ; 153 ;==================================================== 154 ; Check input parameters 155 ;==================================================== 156 ; 157 ; xaxis related parameters 158 ; 159 if n_elements(xaxis) NE 0 then BEGIN 160 CASE (size(xaxis))[0] OF 161 0:nx = 1L 162 1:nx = (size(xaxis))[1] 163 2:nx = (size(xaxis))[1] 164 ENDCASE 165 ENDIF ELSE BEGIN 166 IF n_elements(startx) EQ 0 THEN BEGIN 167 dummy = report('If xaxis is not given, startx must be defined') 168 return 169 ENDIF 170 CASE n_elements(stepxin) OF 171 0:BEGIN 172 dummy = report('If xaxis is not given, stepxin must be defined') 173 return 174 END 175 1:BEGIN 176 IF n_elements(nxin) EQ 0 THEN BEGIN 177 dummy = report('If xaxis is not given and stepxin has only one element, nx must be defined') 178 return 179 ENDIF ELSE nx = nxin 180 END 181 ELSE:nx = n_elements(stepxin) 182 ENDCASE 183 ENDELSE 184 ; 185 ; yaxis related parameters 186 ; 187 if n_elements(yaxis) NE 0 then BEGIN 188 CASE (size(yaxis))[0] OF 189 0:ny = 1L 190 1:ny = (size(yaxis))[1] 191 2:ny = (size(yaxis))[2] 192 ENDCASE 193 ENDIF ELSE BEGIN 194 IF n_elements(starty) EQ 0 THEN BEGIN 195 dummy = report('If yaxis is not given, starty must be defined') 196 return 197 ENDIF 198 CASE n_elements(stepyin) OF 199 0:BEGIN 200 dummy = report('If yaxis is not given, stepyin must be defined') 201 return 202 END 203 1:BEGIN 204 IF n_elements(nyin) EQ 0 THEN BEGIN 205 dummy = report('If yaxis is not given and stepyin has only one element, ny must be defined') 206 return 207 ENDIF ELSE ny = nyin 208 END 209 ELSE:ny = n_elements(stepyin) 210 ENDCASE 211 ENDELSE 212 ; 213 ; zaxis related parameters 214 ; 215 if n_elements(zaxis) NE 0 then BEGIN 216 CASE (size(zaxis))[0] OF 217 0:nz = 1L 218 1:nz = (size(zaxis))[1] 219 ELSE:BEGIN 220 print, 'not coded' 221 stop 222 END 223 ENDCASE 224 ENDIF ELSE nz = 1L 225 ; 226 ;==================================================== 227 ; Others automatic definitions... 228 ;==================================================== 229 ; 230 jpiglo = long(nx) 231 jpjglo = long(ny) 232 jpkglo = long(nz) 233 ; 234 IF n_elements(xminmesh) NE 0 THEN ixminmesh = long(xminmesh[0]) ELSE ixminmesh = 0l 235 IF n_elements(xmaxmesh) NE 0 THEN ixmaxmesh = long(xmaxmesh[0]) ELSE ixmaxmesh = jpiglo-1 236 IF n_elements(yminmesh) NE 0 THEN iyminmesh = long(yminmesh[0]) ELSE iyminmesh = 0l 237 IF n_elements(ymaxmesh) NE 0 THEN iymaxmesh = long(ymaxmesh[0]) ELSE iymaxmesh = jpjglo-1 238 IF n_elements(zminmesh) NE 0 THEN izminmesh = long(zminmesh[0]) ELSE izminmesh = 0l 239 IF n_elements(zmaxmesh) NE 0 THEN izmaxmesh = long(zmaxmesh[0]) ELSE izmaxmesh = jpkglo-1 240 ; 241 iymaxmesh = iymaxmesh-keyword_set(fbase2tbase) 242 ; 243 IF ixmaxmesh LT 0 THEN ixmaxmesh = jpiglo -1 + ixmaxmesh 244 IF iymaxmesh LT 0 THEN iymaxmesh = jpjglo -1 + iymaxmesh 245 IF izmaxmesh LT 0 THEN izmaxmesh = jpkglo -1 + izmaxmesh 246 ; avoid basics errors... 247 ixmaxmesh = 0 > ixmaxmesh < (jpiglo-1) 248 ixminmesh = 0 > ixminmesh < ixmaxmesh 249 iymaxmesh = 0 > iymaxmesh < (jpjglo-1) 250 iyminmesh = 0 > iyminmesh < iymaxmesh 251 izmaxmesh = 0 > izmaxmesh < (jpkglo-1) 252 izminmesh = 0 > izminmesh < izmaxmesh 253 ; 254 jpi = ixmaxmesh-ixminmesh+1 255 jpj = iymaxmesh-iyminmesh+1 256 jpk = izmaxmesh-izminmesh+1 257 ; 258 jpidta = jpiglo 259 jpjdta = jpjglo 260 jpkdta = jpkglo 261 ixmindta = 0 262 ixmaxdta = jpidta-1 263 iymindta = 0 264 iymaxdta = jpjdta-1 265 izmindta = 0 266 izmaxdta = jpkdta-1 267 ; 268 key_partialstep = 0 269 if n_elements(stride) eq 3 then key_stride = stride $ 270 ELSE key_stride = [1, 1, 1] 271 key_gridtype = 'c' 272 ; 273 ; check xyindex and its consequences 274 ; 275 if keyword_set(xyindex) then onearth = 0 276 ; 277 ; check onearth and its consequences 278 ; 279 IF n_elements(onearth) EQ 0 THEN key_onearth = 1b $ 280 ELSE key_onearth = keyword_set(onearth) 281 IF NOT key_onearth THEN BEGIN 282 periodic = 0 283 shift = 0 284 ENDIF; 103 285 286 r = 6371000. 287 ; 288 ;==================================================== 289 ; X direction : glamt 290 ;==================================================== 291 ; 292 ; def of glamt 293 ; 294 if n_elements(xaxis) NE 0 then BEGIN 295 if keyword_set(xyindex) THEN glamt = findgen(jpiglo) ELSE glamt = xaxis 296 ENDIF ELSE BEGIN 297 if keyword_set(xyindex) THEN stepx = 1. ELSE stepx = stepxin 298 CASE 1 OF 299 n_elements(stepx):glamt = startx + findgen(jpiglo)*stepx 300 size(stepx, /n_dimensions):glamt = startx + total(stepx, /cumulative) 301 ELSE:BEGIN 302 dummy = report('Wrong definition of stepx...') 303 return 304 END 305 ENDCASE 306 ENDELSE 307 ; 308 ; apply glamboundary 309 ; 310 IF keyword_set(glamboundary) AND key_onearth THEN BEGIN 311 IF glamboundary[0] GE glamboundary[1] THEN stop 312 IF glamboundary[1]-glamboundary[0] GT 360 THEN stop 313 glamt = glamt MOD 360 314 smaller = where(glamt LT glamboundary[0]) 315 if smaller[0] NE -1 then glamt[smaller] = glamt[smaller]+360. 316 bigger = where(glamt GE glamboundary[1]) 317 if bigger[0] NE -1 then glamt[bigger] = glamt[bigger]-360. 318 ENDIF 319 ; 320 ; force glamt to have 2 dimensions 321 ; 322 CASE size(reform(glamt), /n_dimensions) OF 323 0:glamt = replicate(glamt, jpi, jpj) 324 1:glamt = glamt[ixminmesh:ixmaxmesh]#replicate(1, jpj) 325 2:glamt = glamt[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh] 326 ENDCASE 327 ; keep 2d array even with degenereted dimension 328 IF jpj EQ 1 THEN glamt = reform(glamt, jpi, jpj, /over) 329 ; 330 ;==================================================== 331 ; Y direction : gphit 332 ;==================================================== 333 ; 334 ; def of gphit 335 ; 336 if n_elements(yaxis) NE 0 THEN BEGIN 337 if keyword_set(xyindex) THEN gphit = findgen(jpjglo) ELSE gphit = yaxis 338 ENDIF ELSE BEGIN 339 if keyword_set(xyindex) THEN stepy = 1. ELSE stepy = stepyin 340 CASE 1 OF 341 n_elements(stepy):gphit = starty + findgen(jpjglo)*stepy 342 size(stepy, /n_dimensions):gphit = starty + total(stepy, /cumulative) 343 ELSE:BEGIN 344 dummy = report('Wrong definition of stepy...') 345 return 346 END 347 ENDCASE 348 ENDELSE 349 ; 350 ; force gphit to have 2 dimensions 351 ; 352 CASE size(reform(gphit), /n_dimensions) OF 353 0:gphit = replicate(gphit, jpi, jpj) 354 1:gphit = replicate(1, jpi)#gphit[iyminmesh:iymaxmesh] 355 2:gphit = gphit[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh] 356 ENDCASE 357 ; keep 2d array even with degenereted dimension 358 IF jpj EQ 1 THEN gphit = reform(gphit, jpi, jpj, /over) 359 ; 360 ;==================================================== 361 ; def of key_shift 362 ;==================================================== 363 ; 364 ; definition of key_shift by shiftting the array to have the min 365 ; values of glamt[*, 0] in glamt[0, 0] 366 ; 367 IF n_elements(shift) EQ 0 THEN BEGIN 368 IF jpi GT 1 then BEGIN 369 xtest = glamt[*, 0] 370 key_shift = (where(xtest EQ min(xtest)))[0] 371 IF key_shift NE 0 THEN key_shift = jpi - key_shift 372 ENDIF ELSE key_shift = 0 373 ENDIF ELSE key_shift = shift 374 ; 375 ;==================================================== 376 ; def of key_periodic 377 ;==================================================== 378 ; 379 IF n_elements(periodic) EQ 0 THEN BEGIN 380 IF jpi GT 1 THEN BEGIN 381 xtest = shift(glamt[*, 0], key_shift) 382 ; check that xtest is now sorted in the increasing order 383 IF array_equal(sort(xtest), lindgen(jpi)) EQ 0 THEN BEGIN 384 print, 'WARNING: we cannot sort the xaxis with a simple shift...' 385 print, 'we force key_periodic = 0 and key_shift = 0' 386 print, 'only horizontal plot may be ok...' 387 key_periodic = 0 388 xnotsorted = 1 389 ENDIF ELSE BEGIN 390 key_periodic = (xtest[jpi-1]+2*(xtest[jpi-1]-xtest[jpi-2])) $ 391 GE (xtest[0]+360) 392 ENDELSE 393 ENDIF ELSE key_periodic = 0 394 ENDIF ELSE key_periodic = keyword_set(periodic) 395 ; 396 ; update key_shift 397 ; 398 key_shift = key_shift * (key_periodic EQ 1) 399 ; 400 IF NOT keyword_set(key_periodic) AND keyword_set(fbase2tbase) THEN BEGIN 401 ixmaxmesh = ixmaxmesh-1 402 jpi = jpi-1 403 ENDIF 404 ; 405 ;==================================================== 406 ; apply key_shift 407 ;==================================================== 408 ; 409 if keyword_set(key_shift) then BEGIN 410 glamt = shift(glamt, key_shift, 0) 411 gphit = shift(gphit, key_shift, 0) 412 IF jpj EQ 1 THEN BEGIN 413 glamt = reform(glamt, jpi, jpj, /over) 414 gphit = reform(gphit, jpi, jpj, /over) 415 ENDIF 416 ENDIF 417 ; 418 ;==================================================== 419 ; def key_yreverse 420 ;==================================================== 421 ; 422 IF jpj GT 1 THEN BEGIN 423 if gphit[0, 1] LT gphit[0, 0] then begin 424 key_yreverse = 1 425 gphit = reverse(gphit, 2) 426 glamt = reverse(glamt, 2) 427 ENDIF ELSE key_yreverse = 0 428 ENDIF ELSE key_yreverse = 0 429 ; 430 ;==================================================== 431 ; Are we using a "regular" grid (that can be described 432 ; with x vector and y vector)? 433 ;==================================================== 434 ; 435 ; to get faster, we first test the most basic cases before 436 ; testing the full array. 437 ; 438 CASE 1 OF 439 keyword_set(xyindex):key_irregular = 0b 440 jpi EQ 1 OR jpj EQ 1:key_irregular = 0b 441 n_elements(xaxis) EQ 0 AND n_elements(yaxis) EQ 0:key_irregular = 0b 442 size(reform(xaxis), /n_dimensions) EQ 1 AND size(reform(xaxis), /n_dimensions) EQ 1:key_irregular = 0b 443 n_elements(xaxis) EQ 0 AND size(reform(yaxis), /n_dimensions) EQ 1:key_irregular = 0b 444 n_elements(yaxis) EQ 0 AND size(reform(xaxis), /n_dimensions) EQ 1:key_irregular = 0b 445 array_equal(glamt[*, 0], glamt[*, jpj-1]) EQ 0:key_irregular = 1b 446 array_equal(gphit[0, *], gphit[jpi-1, *]) EQ 0:key_irregular = 1b 447 array_equal(glamt, glamt[*, 0]#replicate(1, jpj)) EQ 0:key_irregular = 1b 448 array_equal(gphit, replicate(1, jpi)#(gphit[0, *])[*]) EQ 0:key_irregular = 1b 449 ELSE:key_irregular = 0b 450 ENDCASE 451 ; 452 ;==================================================== 453 ; def of glamf: defined as the middle of T(i,j) T(i+1,j+1) 454 ;==================================================== 455 ; 456 IF jpi GT 1 THEN BEGIN 457 ; we must compute stepxf: x distance between T(i,j) T(i+1,j+1) 458 CASE 1 OF 459 n_elements(stepx):stepxf = stepx 460 size(stepx, /n_dimensions):stepxf = stepx#replicate(1, jpj) 461 ELSE:BEGIN 462 if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $ 463 OR (keyword_set(key_periodic) AND key_irregular) then BEGIN 464 stepxf = (glamt + 720) MOD 360 465 IF jpj EQ 1 THEN stepxf = reform(stepxf, jpi, jpj, /over) 466 stepxf = shift(stepxf, -1, -1) - stepxf 467 stepxf = [ [[stepxf]], [[stepxf + 360]], [[stepxf - 360]] ] 468 stepxf = min(abs(stepxf), dimension = 3) 469 IF NOT keyword_set(key_periodic) THEN $ 470 stepxf[jpi-1, *] = stepxf[jpi-2, *] 471 ENDIF ELSE BEGIN 472 stepxf = shift(glamt, -1, -1) - glamt 473 IF keyword_set(key_periodic) THEN $ 474 stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $ 475 ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *] 476 ENDELSE 477 IF jpj GT 1 THEN BEGIN 478 stepxf[*, jpj-1] = stepxf[*, jpj-2] 479 stepxf[jpi-1, jpj-1] = stepxf[jpi-2, jpj-2] 480 ENDIF 481 END 482 ENDCASE 483 glamf = glamt + 0.5 * stepxf 484 IF jpj EQ 1 THEN glamf = reform(glamf, jpi, jpj, /over) 485 ENDIF ELSE glamf = glamt + 0.5 486 IF keyword_set(glamboundary) AND key_onearth THEN $ 487 glamf = glamboundary[0] > temporary(glamf) < glamboundary[1] 488 ; 489 ;==================================================== 490 ; def of gphif: defined as the middle of T(i,j) T(i+1,j+1) 491 ;==================================================== 492 ; 493 IF jpj GT 1 THEN BEGIN 494 ; we must compute stepyf: y distance between T(i,j) T(i+1,j+1) 495 CASE 1 OF 496 n_elements(stepy):stepyf = stepy 497 size(stepy, /n_dimensions):stepyf = replicate(1, jpi)#stepy 498 ELSE:BEGIN 499 stepyf = shift(gphit, -1, -1) - gphit 500 stepyf[*, jpj-1] = stepyf[*, jpj-2] 501 IF jpi GT 1 THEN BEGIN 502 if NOT keyword_set(key_periodic) THEN $ 503 stepyf[jpi-1, *] = stepyf[jpi-2, *] 504 stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2] 505 ENDIF 506 END 507 ENDCASE 508 gphif = gphit + 0.5 * stepyf 509 ENDIF ELSE gphif = gphit + 0.5 510 IF key_onearth THEN gphif = -90. > gphif < 90. 511 ; 512 ;==================================================== 513 ; e1t: x distance between U(i-1,j) and U(i,j) 514 ;==================================================== 515 ; 516 ; *-|-*---|---*---| 517 ; 518 IF jpi GT 1 THEN BEGIN 519 IF n_elements(stepx) NE 1 THEN BEGIN 520 IF keyword_set(irregular) THEN BEGIN 521 ; we must compute stepxu: x distance between T(i,j) T(i+1,j) 522 IF keyword_set(key_periodic) THEN BEGIN 523 stepxu = (glamt + 720) MOD 360 524 stepxu = shift(stepxu, -1, 0) - stepxu 525 stepxu = [ [[stepxu]], [[stepxu + 360]], [[stepxu - 360]] ] 526 stepxu = min(abs(stepxu), dimension = 3) 527 ENDIF ELSE BEGIN 528 stepxu = shift(glamt, -1, 0) - glamt 529 stepxu[jpi-1, *] = stepxf[jpi-2, *] 530 ENDELSE 531 ENDIF ELSE stepxu = stepxf 532 IF jpj EQ 1 THEN stepxu = reform(stepxu, jpi, jpj, /over) 533 e1t = 0.5*(stepxu+shift(stepxu, 1, 0)) 534 IF NOT keyword_set(key_periodic) THEN $ 535 e1t[0, *] = e1t[1, *] 536 IF jpj EQ 1 THEN e1t = reform(e1t, jpi, jpj, /over) 537 ENDIF ELSE e1t = replicate(stepx, jpi, jpj) 538 ENDIF ELSE e1t = replicate(1b, jpi, jpj) 539 ; 540 ;==================================================== 541 ; e2t: y distance between V(i,j-1) and V(i,j) 542 ;==================================================== 543 ; 544 IF jpj GT 1 THEN BEGIN 545 ; we must compute stepyv: y distance between T(i,j) T(i,j+1) 546 IF n_elements(stepy) NE 1 THEN BEGIN 547 IF keyword_set(key_irregular) THEN BEGIN 548 stepyv = shift(gphit, 0, -1) - gphit 549 stepyv[*, jpj-1] = stepyv[*, jpj-2] 550 ENDIF ELSE stepyv = stepyf 551 e2t = 0.5*(stepyv+shift(stepyv, 0, 1)) 552 e2t[*, 0] = e2t[*, 1] 553 ENDIF ELSE e2t = replicate(stepy, jpi, jpj) 554 ENDIF ELSE e2t = replicate(1b, jpi, jpj) 555 ; 556 IF key_onearth THEN e2t = r * !pi/180. * temporary(e2t) 557 ; 558 ;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 559 IF keyword_set(fullcgrid) THEN BEGIN 560 ; 561 ;==================================================== 562 ; def of glamu: defined as the middle of T(i,j) T(i+1,j) 563 ;==================================================== 564 ; 565 IF keyword_set(irregular) THEN BEGIN 566 glamu = glamt + 0.5 * stepxu 567 IF keyword_set(glamboundary) AND key_onearth THEN $ 568 glamu = glamboundary[0] > temporary(glamu) < glamboundary[1] 569 ENDIF ELSE glamu = glamf 570 ; 571 ;==================================================== 572 ; def of gphiu: defined as the middle of T(i,j) T(i+1,j) 573 ;==================================================== 574 ; 575 IF jpi GT 1 THEN BEGIN 576 ; we must compute stepyu: y distance between T(i+1,j) T(i,j) 577 IF keyword_set(key_irregular) THEN BEGIN 578 stepyu = shift(gphit, -1, 0) - gphit 579 IF NOT keyword_set(key_periodic) THEN $ 580 stepyu[jpi-1, *] = stepyu[jpi-2, *] 581 gphiu = gphit + 0.5 * stepyu 582 IF jpj EQ 1 THEN gphiu = reform(gphiu, jpi, jpj, /over) 583 ENDIF ELSE gphiu = gphit 584 ENDIF ELSE gphiu = gphit 585 IF key_onearth THEN gphiu = -90. > gphiu < 90. 586 ; 587 ;==================================================== 588 ; def of glamv: defined as the middle of T(i,j) T(i,j+1) 589 ;==================================================== 590 ; 591 IF jpj GT 1 THEN BEGIN 592 ; we must compute stepxv: x distance between T(i,j) T(i,j+1) 593 IF keyword_set(irregular) THEN BEGIN 594 IF keyword_set(key_periodic) THEN BEGIN 595 stepxv = (glamt + 720) MOD 360 596 stepxv = shift(stepxv, 0, -1) - stepxv 597 stepxv = [ [[stepxv]], [[stepxv + 360]], [[stepxv - 360]] ] 598 stepxv = min(abs(stepxv), dimension = 3) 599 ENDIF ELSE stepxv = shift(glamt, 0, -1) - glamt 600 stepxv[*, jpj-1] = stepxv[*, jpj-2] 601 glamv = glamt + 0.5 * stepxv 602 IF keyword_set(glamboundary) AND key_onearth THEN $ 603 glamv = glamboundary[0] > temporary(glamv) < glamboundary[1] 604 ENDIF ELSE glamv = glamt 605 ENDIF ELSE glamv = glamt 606 ; 607 ;==================================================== 608 ; def of gphiv: defined as the middle of T(i,j) T(i,j+1) 609 ;==================================================== 610 ; 611 IF keyword_set(key_irregular) THEN $ 612 gphiv = gphit + 0.5 * stepyv $ 613 ELSE gphiv = gphif 614 IF key_onearth THEN gphiv = -90. > gphiv < 90. 615 ; 616 ;==================================================== 617 ; e1u: x distance between T(i,j) and T(i+1,j) 618 ;==================================================== 619 ; 620 IF jpi GT 1 AND n_elements(stepx) NE 1 THEN $ 621 e1u = stepxu ELSE e1u = e1t 622 ; 623 ;==================================================== 624 ; e2u: y distance between F(i,j-1) and F(i,j) 625 ;==================================================== 626 ; 627 IF keyword_set(key_irregular) THEN BEGIN 628 e2u = gphif - shift(gphif, 0, 1) 629 e2u[*, 0] = e2u[*, 1] 630 IF key_onearth THEN e2u = r * !pi/180. * temporary(e2u) 631 ENDIF ELSE e2u = e2t 632 ; 633 ;==================================================== 634 ; e1v: x distance between F(i-1,j) and F(i,j) 635 ;==================================================== 636 ; 637 IF keyword_set(irregular) THEN BEGIN 638 IF keyword_set(key_periodic) THEN BEGIN 639 e1v = (glamf + 720) MOD 360 640 e1v = e1v - shift(e1v, 1, 0) 641 e1v = [ [[e1v]], [[e1v + 360]], [[e1v - 360]] ] 642 e1v = min(abs(e1v), dimension = 3) 643 ENDIF ELSE BEGIN 644 e1v = glamf - shift(glamf, 1, 0) 645 e1v[0, *] = stepxf[1, *] 646 ENDELSE 647 ENDIF ELSE e1v = e1t 648 ; 649 ;==================================================== 650 ; e2v: y distance between T(i,j) and T(i+1,j) 651 ;==================================================== 652 ; 653 IF jpj GT 1 and n_elements(stepy) NE 1 THEN BEGIN 654 e2v = stepyv 655 IF key_onearth THEN e2v = r * !pi/180. * temporary(e2v) 656 ENDIF ELSE e2v = e2t 657 ; 658 ;==================================================== 659 ; e1f: x distance between V(i,j) and V(i+1,j) 660 ;==================================================== 661 ; 662 IF keyword_set(irregular) THEN BEGIN 663 IF keyword_set(key_periodic) THEN BEGIN 664 e1f = (glamv + 720) MOD 360 665 e1f = shift(e1f, -1, 0) - e1f 666 e1f = [ [[e1f]], [[e1f + 360]], [[e1f - 360]] ] 667 e1f = min(abs(e1f), dimension = 3) 668 ENDIF ELSE BEGIN 669 e1f = shift(glamv, -1, 0) - glamt 670 e1f[jpi-1, *] = stepxf[jpi-2, *] 671 ENDELSE 672 ENDIF ELSE e1f = e1u 673 ; 674 ;==================================================== 675 ; e2f: y distance between U(i,j) and U(i,j+1) 676 ;==================================================== 677 ; 678 IF keyword_set(key_irregular) THEN BEGIN 679 e2f = shift(gphiu, 0, -1) - gphiu 680 e2f[*, jpj-1] = e2f[*, jpj-2] 681 IF key_onearth THEN e2f = r * !pi/180. * temporary(e2f) 682 ENDIF ELSE e2f = e2v 683 ; 684 ENDIF 685 ;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 686 ; 687 ; 688 ;==================================================== 689 ; e1[tuvf] from degree to meters 690 ;==================================================== 691 ; 692 IF keyword_set(key_onearth) THEN BEGIN 693 e1t = r * !pi/180. * temporary(e1t) * cos(!pi/180.*gphit) 694 IF keyword_set(fullcgrid) THEN BEGIN 695 e1u = r * !pi/180. * temporary(e1u) * cos(!pi/180.*gphiu) 696 e1v = r * !pi/180. * temporary(e1v) * cos(!pi/180.*gphiv) 697 e1f = r * !pi/180. * temporary(e1f) * cos(!pi/180.*gphif) 698 ENDIF 699 ENDIF 700 ; 701 ;==================================================== 702 ; if not fullcgrid: make sure we don't use glam[uv], gphi[uv], e[12][uvf] 703 ;==================================================== 704 ; 705 IF NOT keyword_set(fullcgrid) THEN BEGIN 706 glamu = !values.f_nan & glamv = !values.f_nan 707 gphiu = !values.f_nan & gphiv = !values.f_nan 708 e1u = !values.f_nan & e1v = !values.f_nan & e1f = !values.f_nan 709 e2u = !values.f_nan & e2v = !values.f_nan & e2f = !values.f_nan 710 firstxu = !values.f_nan & lastxu = !values.f_nan & nxu = !values.f_nan 711 firstyu = !values.f_nan & lastyu = !values.f_nan & nyu = !values.f_nan 712 firstxv = !values.f_nan & lastxv = !values.f_nan & nxv = !values.f_nan 713 firstyv = !values.f_nan & lastyv = !values.f_nan & nyv = !values.f_nan 714 ENDIF 715 ; 716 ;==================================================== 717 ; Z direction 718 ;==================================================== 719 ; 720 ; z axis 721 ; 722 CASE n_elements(zaxis) OF 723 0:BEGIN 724 gdept = 0. 725 key_zreverse = 0 726 END 727 1:BEGIN 728 gdept = zaxis 729 key_zreverse = 0 730 END 731 ELSE:BEGIN 732 gdept = zaxis[izminmesh:izmaxmesh] 733 IF jpk GT 1 THEN BEGIN 734 if gdept[0] GT gdept[1] then begin 735 gdept = reverse(gdept) 736 key_zreverse = 1 737 ENDIF ELSE key_zreverse = 0 738 ENDIF ELSE key_zreverse = 0 739 END 740 ENDCASE 741 ; 742 if n_elements(gdept) GT 1 then BEGIN 743 stepz = shift(gdept, -1)-gdept 744 stepz[jpk-1] = stepz[jpk-2] 745 gdepw = 0. > (gdept-stepz/2.) 746 ENDIF ELSE BEGIN 747 stepz = 1. 748 gdepw = gdept 749 ENDELSE 750 ; 751 ;==================================================== 752 ; e3[tw]: 753 ;==================================================== 754 ; 755 e3t = stepz 756 IF n_elements(stepz) GT 1 THEN BEGIN 757 e3w = 0.5*(stepz+shift(stepz, 1)) 758 e3w[0] = 0.5*e3t[0] 759 ENDIF ELSE e3w = e3t 760 ; 761 ;==================================================== 762 ; Mask 763 ;==================================================== 764 ; 765 ; defaut mask eq 1 766 if NOT keyword_set(mask) then mask = -1 767 ; 768 if mask[0] NE -1 then BEGIN 769 tmask = byte(mask[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh, izminmesh:izmaxmesh]) 770 tmask = reform(tmask, jpi, jpj, jpk, /over) 771 if key_shift NE 0 then tmask = shift(tmask, key_shift, 0, 0) 772 ; because tmask = reverse(tmask, 2) is not working if the 3rd 773 ; dimension of tmask = 1, we call reform. 774 IF jpk EQ 1 THEN tmask = reform(tmask, /over) 775 IF key_yreverse EQ 1 THEN tmask = reverse(tmask, 2) 776 IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over) 777 IF key_zreverse EQ 1 THEN tmask = reverse(tmask, 3) 778 IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over) 779 IF keyword_set(fullcgrid) THEN BEGIN 780 IF keyword_set(key_periodic) THEN BEGIN 781 msk = tmask*shift(tmask, -1, 0, 0) 782 umaskred = msk[jpi-1, *, *] 783 ENDIF ELSE umaskred = tmask[jpi-1, *, *] 784 vmaskred = tmask[*, jpj-1, *] 785 fmaskredy = tmask[jpi-1, *, *] 786 fmaskredx = tmask[*, jpj-1, *] 787 ENDIF 788 ENDIF ELSE BEGIN 789 tmask = replicate(1b, jpi, jpj, jpk) 790 IF keyword_set(fullcgrid) THEN BEGIN 791 umaskred = replicate(1b, jpj, jpk) 792 vmaskred = replicate(1b, jpi, jpk) 793 fmaskredy = replicate(1b, jpj, jpk) 794 fmaskredx = replicate(1b, jpi, jpk) 795 ENDIF 796 ENDELSE 797 ; 798 IF NOT keyword_set(fullcgrid) THEN BEGIN 799 umaskred = !values.f_nan 800 vmaskred = !values.f_nan 801 fmaskredy = !values.f_nan 802 fmaskredx = !values.f_nan 803 ENDIF 804 ; 805 ;==================================================== 806 ; stride... 807 ;==================================================== 808 ; 809 IF total(key_stride) GT 3 THEN BEGIN 810 IF key_shift NE 0 THEN BEGIN 811 ; for explanation, see header of read_ncdf_varget.pro 812 jpiright = key_shift 813 jpileft = jpi - key_shift - ( (key_stride[0]-1)-((key_shift-1) MOD key_stride[0]) ) 814 jpi = ((jpiright-1)/key_stride[0]+1) + ((jpileft-1)/key_stride[0]+1) 815 ENDIF ELSE jpi = (jpi-1)/key_stride[0]+1 816 jpj = (jpj-1)/key_stride[1]+1 817 jpk = (jpk-1)/key_stride[2]+1 818 ; 819 glamt = (temporary(glamt))[0:*:stride[0], 0:*:stride[1]] 820 gphit = (temporary(gphit))[0:*:stride[0], 0:*:stride[1]] 821 e1t = (temporary(e1t))[0:*:stride[0], 0:*:stride[1]] 822 e2t = (temporary(e2t))[0:*:stride[0], 0:*:stride[1]] 823 tmask = (temporary(tmask))[0:*:stride[0], 0:*:stride[1], 0:*:stride[2]] 824 gdept = gdept[0:*:stride[2]] 825 gdepw = gdepw[0:*:stride[2]] 826 e3t = e3t[0:*:stride[2]] 827 e3w = e3w[0:*:stride[2]] 828 ; we must recompute glamf and gphif... 829 IF jpi GT 1 THEN BEGIN 830 if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $ 831 OR (keyword_set(key_periodic) AND key_irregular) then BEGIN 832 stepxf = (glamt + 720) MOD 360 833 stepxf = shift(stepxf, -1, -1) - stepxf 834 stepxf = [ [[stepxf]], [[stepxf + 360]], [[stepxf - 360]] ] 835 stepxf = min(abs(stepxf), dimension = 3) 836 IF NOT keyword_set(key_periodic) THEN $ 837 stepxf[jpi-1, *] = stepxf[jpi-2, *] 838 ENDIF ELSE BEGIN 839 stepxf = shift(glamt, -1, -1) - glamt 840 IF keyword_set(key_periodic) THEN $ 841 stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $ 842 ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *] 843 ENDELSE 844 IF jpj GT 1 THEN BEGIN 845 stepxf[*, jpj-1] = stepxf[*, jpj-2] 846 stepxf[jpi-1, jpj-1] = stepxf[jpi-2, jpj-2] 847 ENDIF 848 glamf = glamt + 0.5 * stepxf 849 IF jpj EQ 1 THEN glamf = reform(glamf, jpi, jpj, /over) 850 ENDIF ELSE glamf = glamt + 0.5 851 IF jpj GT 1 THEN BEGIN 852 ; we must compute stepyf: y distance between T(i,j) T(i+1,j+1) 853 stepyf = shift(gphit, -1, -1) - gphit 854 stepyf[*, jpj-1] = stepyf[*, jpj-2] 855 IF jpi GT 1 THEN BEGIN 856 if NOT keyword_set(key_periodic) THEN $ 857 stepyf[jpi-1, *] = stepyf[jpi-2, *] 858 stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2] 859 ENDIF 860 gphif = gphit + 0.5 * stepyf 861 ENDIF ELSE gphif = gphit + 0.5 862 ; 863 IF keyword_set(fullcgrid) THEN BEGIN 864 glamu = (temporary(glamu))[0:*:stride[0], 0:*:stride[1]] 865 gphiu = (temporary(gphiu))[0:*:stride[0], 0:*:stride[1]] 866 e1u = (temporary(e1u))[0:*:stride[0], 0:*:stride[1]] 867 e2u = (temporary(e2u))[0:*:stride[0], 0:*:stride[1]] 868 glamv = (temporary(glamv))[0:*:stride[0], 0:*:stride[1]] 869 gphiv = (temporary(gphiv))[0:*:stride[0], 0:*:stride[1]] 870 e1v = (temporary(e1v))[0:*:stride[0], 0:*:stride[1]] 871 e2v = (temporary(e2v))[0:*:stride[0], 0:*:stride[1]] 872 e1f = (temporary(e1f))[0:*:stride[0], 0:*:stride[1]] 873 e2f = (temporary(e2f))[0:*:stride[0], 0:*:stride[1]] 874 umaskred = (temporary(umaskred))[0, 0:*:stride[1], 0:*:stride[2]] 875 vmaskred = (temporary(vmaskred))[0:*:stride[0], 0, 0:*:stride[2]] 876 fmaskredy = (temporary(fmaskredy))[0, 0:*:stride[1], 0:*:stride[2]] 877 fmaskredx = (temporary(fmaskredx))[0:*:stride[0], 0, 0:*:stride[2]] 878 ENDIF 879 ENDIF 880 ; 881 ;==================================================== 882 ; apply all the grid parameters 883 ;==================================================== 884 ; 885 @updateold 886 domdef 887 ; 888 ;==================================================== 889 ; Triangulation 890 ;==================================================== 891 ; 892 IF total(mask) EQ jpi*jpj*jpk $ 893 AND NOT keyword_set(key_irregular) THEN triangles_list = -1 $ 894 ELSE triangles_list = triangule(/keep_cont) 895 ; 896 ;==================================================== 897 ; time axis (default definition) 898 ;==================================================== 899 ; 900 jpt = 1 901 time = 0 902 ; 903 IF NOT keyword_set(key_forgetold) THEN BEGIN 904 @updateold 905 ENDIF 104 906 ;------------------------------------------------------------ 105 return 907 IF keyword_set(key_performance) EQ 1 THEN $ 908 print, 'time computegrid', systime(1)-time1 909 ;------------------------------------------------------------ 910 return 106 911 end 107 912 108 -
trunk/Obsolete/meshlec.pro
r12 r13 330 330 endif 331 331 ;------------------------------------------------------- 332 key_yreverse = 0 333 key_zreverse = 0 334 key_partialstep = 0 335 key_stride = [1, 1, 1] 336 key_gridtype = 'c' 337 ; 332 338 if not keyword_set(pasblabla) then print,'lecture '+nomfich+' finie' 333 339 widget_control, noticebase, bad_id = toto, /destroy -
trunk/Obsolete/ncdf_meshlec.pro
r12 r13 5 5 ; NAME:ncdf_meshlec 6 6 ; 7 ; PURPOSE:lit le fichier meshmask au format netcdf cree par OPA ds sa 8 ; version tout NETCDF! 7 ; PURPOSE:obsolete, use ncdf_meshread instead... 9 8 ; 10 ; CATEGORY:lecture de grille9 ; MODIFICATION HISTORY: 11 10 ; 12 ; CALLING SEQUENCE:ncdf_meshlec [,' nomfich'] 13 ; 14 ; INPUTS: 15 ; 16 ; nomfich: un string definissant le nom du fichier a lire. Si 17 ; ce string ne contient pas l''adresse entiere du fichier nomfich est 18 ; complete par iodir. Par defaut nomfich=iodir+'meshmask.nc' 19 ; 20 ; KEYWORD PARAMETERS: 21 ; 22 ; GLAMBOUNDARY: un vecteur de 2 elements specifaint le min et le 23 ; max qui doivent etre imposes en longitude (obligatoire si le 24 ; tableau depasse 360 degres) 25 ; 26 ; /CHECKDAT: cherche si il existe un ficher du meme nom que le 27 ; fichier meshmask mais terminant par .dat au lieu de .dat. 28 ; Si ce fichier n''existe pas, ncdf_meshlec lit le fichier 29 ; NetCdf et sauve tous les tableaux lus ds le fichier 30 ; ...dat. Rq: Ce fichier est bcp plus petit que le fichier 31 ; NetCdf d''origine car par ex on ne sauve pas les tableaux 32 ; umask, vmask, fmask mais simplement certains de leurs bords 33 ; (cf. umask.pro vmask.pro et fmask.pro) et ces masks sont 34 ; sauves en bytes et toues les autres tableaux en 35 ; float. Cependant ce fichier .dat est cree en fonction des 36 ; parametres rentres ds le init... (ixminmesh.....) et Je ne 37 ; sais rien sur la portabilite du fichier .dat. 38 ; Si ce fichier .dat existe on le lit avec restore. 39 ; 40 ; OUTPUTS:mise a jours des variables de common glam, common 41 ; gphi,common e12,common mask,common tab (cf. common.pro) 42 ; 43 ; COMMON BLOCKS:common.pro 44 ; 45 ; SIDE EFFECTS:prend en compte (si il sont definits) les 46 ; ixminmesh,...et definit jpiglo,jpi,.... 47 ; 48 ; RESTRICTIONS: 49 ; 50 ; La definition de ixminmesh,ixmaxmesh,iyminmesh,iymaxmesh 51 ; ,izminmesh,izmaxmesh doit etre faite avant l''entree dans cette 52 ; routine. pour attribuer automatiquement ces valeurs au maximum 53 ; possible les mettre toutes a -1 et ncdf_meshlec les calculera. 54 ; 55 ; EXAMPLE: 56 ; 57 ; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr) 58 ; 12/1999 11 ; Aug. 2005, Sebastien Masson: switch to ncdf_meshread 59 12 ;- 60 13 ;------------------------------------------------------------ 61 14 ;------------------------------------------------------------ 62 15 ;------------------------------------------------------------ 63 PRO ncdf_meshlec, nomfich, GLAMBOUNDARY = glamboundary, GETDIMENSIONS = getdimensions, CHECKDAT = checkdat 64 @common 65 tempsun = systime(1) ; pour key_performance 66 ;---------------------------------------------------------- 67 ; constitution de l''adresse s_fichier et 68 ; ouverture du fichier a l''adresse s_fichier 69 ;------------------------------------------------------- 70 ; def de nomfich par defaut 71 IF n_params() EQ 0 then nomfich = 'meshmask.nc' 72 ;------------------ 73 if keyword_set(CHECKDAT) then begin 74 nomfichdat = strmid(nomfich, 0, strlen(nomfich)-2)+'dat' 75 thisOS = strupcase(strmid(!version.os_family, 0, 3)) 76 CASE thisOS OF 77 'MAC':sep = ':' 78 'WIN':sep = '\' 79 ELSE:sep = '/' 80 ENDCASE 81 if strpos(nomfichdat, sep) EQ -1 then begin 82 ; si iodir ne finit pas par sep on le rajoute 83 if rstrpos(iodir, sep) NE strlen(iodir)-1 then iodir = iodir+sep 84 nomfichdat = iodir+nomfichdat 85 endif 86 test = findfile(nomfichdat) ; le nom cherche correspond bien a un fichier? 87 if test[0] NE '' then begin 88 noticebase=xnotice('Lecture du fichier !C '+nomfichdat+'!C ...') 89 restore, nomfichdat 90 widget_control, noticebase, bad_id = nothing, /destroy 91 if keyword_set(key_performance) THEN print, 'temps ncdf_meshlec', systime(1)-tempsun 92 return 93 ENDIF 94 endif 95 ;------------------ 96 s_fichier = isafile(file = nomfich, iodir = iodir) 16 PRO ncdf_meshlec, filename, _EXTRA = ex 97 17 ; 98 noticebase=xnotice('Lecture du fichier !C '+s_fichier+'!C ...') 18 CASE n_params() OF 19 0:ncdf_meshread, _EXTRA = ex 20 1:ncdf_meshread, filename, _EXTRA = ex 21 ENDCASE 99 22 ; 100 cdfid=ncdf_open(s_fichier) 101 contient=ncdf_inquire(cdfid) 102 103 ;------------------------------------------------------------ 104 ; differentes dimensions 105 ;------------------------------------------------------------ 106 ncdf_diminq,cdfid,'x',name,jpiglo 107 ncdf_diminq,cdfid,'y',name,jpjglo 108 ncdf_diminq,cdfid,'z',name,jpkglo 109 ; 110 if keyword_set(getdimensions) then begin 111 ncdf_close, cdfid 112 return 113 endif 114 ;------------------------------------------------------- 115 ;------------------------------------------------------- 116 if n_elements(ixminmesh) EQ 0 THEN ixminmesh = 0 117 if n_elements(ixmaxmesh) EQ 0 then ixmaxmesh = jpiglo-1 118 if ixminmesh EQ -1 THEN ixminmesh = 0 119 IF ixmaxmesh EQ -1 then ixmaxmesh = jpiglo-1 120 if n_elements(iyminmesh) EQ 0 THEN iyminmesh = 0 121 IF n_elements(iymaxmesh) EQ 0 then iymaxmesh = jpjglo-1 122 if iyminmesh EQ -1 THEN iyminmesh = 0 123 IF iymaxmesh EQ -1 then iymaxmesh = jpjglo-1 124 if n_elements(izminmesh) EQ 0 THEN izminmesh = 0 125 IF n_elements(izmaxmesh) EQ 0 then izmaxmesh = jpkglo-1 126 if izminmesh EQ -1 THEN izminmesh = 0 127 IF izmaxmesh EQ -1 then izmaxmesh = jpkglo-1 128 ; 129 ixmindtasauve = testvar(var = ixmindta) 130 iymindtasauve = testvar(var = iymindta) 131 izmindtasauve = testvar(var = izmindta) 132 ; 133 ixmindta = 0l 134 iymindta = 0l 135 izmindta = 0l 136 ; 137 jpi = long(ixmaxmesh-ixminmesh+1) 138 jpj = long(iymaxmesh-iyminmesh+1) 139 jpk = long(izmaxmesh-izminmesh+1) 140 ; 141 if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1] 142 key_stride = 1l > long(key_stride) 143 ; 144 jpitotal = jpi 145 key_shift = long(testvar(var = key_shift)) 146 key = long(key_shift MOD jpitotal) 147 if key LT 0 then key = key+jpitotal 148 ixmin = ixminmesh-ixmindta 149 ; 150 jpi = jpi/key_stride[0] 151 jpj = jpj/key_stride[1] 152 jpk = jpk/key_stride[2] 153 ; 154 jpt = 1 155 premiertps = 0 156 ; 157 premierx = 0 158 dernierx = jpi-1 159 premiery = 0 160 derniery = jpj-1 161 premierz = 0 162 dernierz = jpk-1 163 nx = jpi 164 ny = jpj 165 nz = 1 166 izminmeshsauve = izminmesh 167 izminmesh = 0 168 ; 169 ;------------------------------------------------------- 170 ; tableaux 2d: 171 ;------------------------------------------------------- 172 if total(key_stride) EQ 3 then $ 173 namevar = ['glamt','glamu','glamv','glamf','gphit','gphiu','gphiv','gphif','e1t','e1u','e1v','e1f','e2t','e2u','e2v','e2f'] $ 174 ELSE namevar = ['glamt','glamu','glamv','gphit','gphiu','gphiv'] 175 ; 176 for i = 0, n_elements(namevar)-1 do begin 177 varcontient=ncdf_varinq(cdfid, namevar[i]) 178 name = varcontient.name 179 @read_ncdf_varget 180 commande = namevar[i]+'=float(res)' 181 rien = execute(commande) 182 ENDFOR 183 ; 184 if total(key_stride) NE 3 then BEGIN 185 glamf = glamt 186 glamf = [glamf, reform(2.*reform(glamt[jpi-1, *])-reform(glamt[jpi-2, *]), 1, jpj)] 187 glamf = [ [glamf], [reform(2*glamf[*, jpj-1]-glamf[*, jpj-2], jpi+1, 1)] ] 188 glamf = glamf+shift(glamf, -1, -1) 189 glamf = .5*glamf[0:jpi-1, 0:jpj-1] 190 ; 191 gphif = gphit 192 gphif = [gphif, reform(2.*reform(gphit[jpi-1, *])-reform(gphit[jpi-2, *]), 1, jpj)] 193 gphif = [ [gphif], [reform(2*gphif[*, jpj-1]-gphif[*, jpj-2], jpi+1, 1)] ] 194 gphif = gphif+shift(gphif, -1, -1) 195 gphif = .5*gphif[0:jpi-1, 0:jpj-1] 196 ; 197 e1t = replicate(1, jpi, jpj) 198 e2t = e1t 199 e1u = e1t 200 e2u = e1t 201 e1v = e1t 202 e2v = e1t 203 e1f = e1t 204 e2f = e1t 205 ENDIF 206 ;------------------------------------------------------- 207 ; tableaux 3d: 208 ;------------------------------------------------------- 209 nz = jpk 210 izminmesh = izminmeshsauve 211 ; 212 varcontient=ncdf_varinq(cdfid,'tmask') 213 name = varcontient.name 214 @read_ncdf_varget 215 tmask = byte(res) 216 ; 217 varcontient=ncdf_varinq(cdfid,'umask') 218 name = varcontient.name 219 nx = 1 220 premierx = jpi-1 221 @read_ncdf_varget 222 umaskred = byte(res) 223 umaskred = reform(umaskred, /over) 224 ; 225 varcontient=ncdf_varinq(cdfid,'vmask') 226 name = varcontient.name 227 nx = jpi 228 premierx = 0 229 ny = 1 230 premiery = jpj-1 231 @read_ncdf_varget 232 vmaskred = byte(res) 233 vmaskred = reform(vmaskred, /over) 234 ; 235 varcontient=ncdf_varinq(cdfid,'fmask') 236 name = varcontient.name 237 nx = 1 238 premierx = jpi-1 239 ny = jpj 240 premiery = 0 241 @read_ncdf_varget 242 fmaskredy = byte(res) 243 coast = where(fmaskredy NE 0 and fmaskredy NE 1) 244 IF coast[0] NE -1 THEN fmaskredy[coast] = 0b 245 fmaskredy= reform(fmaskredy, /over) 246 nx = jpi 247 premierx = 0 248 ny = 1 249 premiery = jpj-1 250 @read_ncdf_varget 251 fmaskredx = byte(res) 252 coast = where(fmaskredx NE 0 and fmaskredx NE 1) 253 IF coast[0] NE -1 THEN fmaskredx[coast] = 0b 254 fmaskredx = reform(fmaskredx, /over) 255 ;------------------------------------------------------- 256 ; tableaux 1d: 257 ;------------------------------------------------------- 258 namevar = ['e3t','e3w','gdept','gdepw'] 259 for i = 0, n_elements(namevar)-1 do begin 260 varcontient=ncdf_varinq(cdfid, namevar[i]) 261 if n_elements(varcontient.dim) EQ 4 THEN BEGIN 262 commande = 'ncdf_varget,cdfid,namevar[i],'+namevar[i] $ 263 +',offset = [0,0,izminmesh,0], count = [1,1,jpk,1]' 264 if key_stride[2] NE 1 then commande = commande+', stride=[1,1,key_stride[2],1]' 265 ENDIF ELSE BEGIN 266 commande = 'ncdf_varget,cdfid,namevar[i],'+namevar[i] $ 267 +',offset = [izminmesh], count = [jpk]' 268 if key_stride[2] NE 1 then commande = commande+', stride=key_stride[2]' 269 ENDELSE 270 rien = execute(commande) 271 commande = namevar[i]+'=float('+namevar[i]+')' 272 rien = execute(commande) 273 commande = namevar[i]+' = reform('+namevar[i]+', /over)' 274 rien = execute(commande) 275 ENDFOR 276 ;------------------------------------------------------- 277 ncdf_close, cdfid 278 ;------------------------------------------------------- 279 ; bornes de glam qui ne doivent pas depasser 360 degres... 280 ;------------------------------------------------------- 281 ; minglam = min(glamt, max = maxglam) 282 ; if maxglam-minglam GE 360 AND NOT keyword_set(glamboundary) then $ 283 ; nothing = execute('glamboundary = '+xquestion('What are the longitudes boundary?', '[-180,180]',/chkwidget)) 284 ;------------------------------------------------------- 285 if keyword_set(glamboundary) then BEGIN 286 if glamboundary[0] NE glamboundary[1] then BEGIN 287 glamt = glamt MOD 360 288 smaller = where(glamt LT glamboundary[0]) 289 if smaller[0] NE -1 then glamt[smaller] = glamt[smaller]+360 290 bigger = where(glamt GE glamboundary[1]) 291 if bigger[0] NE -1 then glamt[bigger] = glamt[bigger]-360 292 glamu = glamu MOD 360 293 smaller = where(glamu LT glamboundary[0]) 294 if smaller[0] NE -1 then glamu[smaller] = glamu[smaller]+360 295 bigger = where(glamu GE glamboundary[1]) 296 if bigger[0] NE -1 then glamu[bigger] = glamu[bigger]-360 297 glamv = glamv MOD 360 298 smaller = where(glamv LT glamboundary[0]) 299 if smaller[0] NE -1 then glamv[smaller] = glamv[smaller]+360 300 bigger = where(glamv GE glamboundary[1]) 301 if bigger[0] NE -1 then glamv[bigger] = glamv[bigger]-360 302 glamf = glamf MOD 360 303 smaller = where(glamf LT glamboundary[0]) 304 if smaller[0] NE -1 then glamf[smaller] = glamf[smaller]+360 305 bigger = where(glamf GE glamboundary[1]) 306 if bigger[0] NE -1 then glamf[bigger] = glamf[bigger]-360 307 endif 308 endif 309 ;------------------------------------------------------- 310 ixmindta = ixmindtasauve 311 iymindta = iymindtasauve 312 izmindta = izmindtasauve 313 ;------------------------------------------------------- 314 widget_control, noticebase, bad_id = nothing, /destroy 315 ; 316 if keyword_set(checkdat) then BEGIN 317 noticebase=xnotice('Ecriture du fichier !C '+nomfichdat+'!C ...') 318 ; 319 save, jpi, jpj, jpk, ixminmesh, ixmaxmesh, iyminmesh, iymaxmesh, izminmesh, izmaxmesh, glamt,glamu,glamv,glamf,gphit,gphiu,gphiv,gphif,e1t,e1u,e1v,e1f,e2t,e2u,e2v,e2f, tmask, umaskred, vmaskred, fmaskredx, fmaskredy, gdept, gdepw, e3t, e3w, filename = nomfichdat 320 ; 321 widget_control, noticebase, bad_id = nothing, /destroy 322 endif 323 ; 324 if keyword_set(key_performance) THEN print, 'temps ncdf_meshlec', systime(1)-tempsun 325 326 ;------------------------------------------------------- 327 return 328 end 23 return 24 END -
trunk/ToBeReviewed/GRILLE/changegrid.pro
r12 r13 28 28 if newgrid.filetype EQ 'batch file' THEN BEGIN 29 29 createpro, '@'+strmid(newgrid.filename[0], 0, strlen(newgrid.filename)-4) $ 30 , filename = isadirectory(io = homedir, title = 'Bad definition of Homedir') $ 31 +'for_createpro.pro' 30 , filename = myuniquetmpdir +'for_createpro.pro' 32 31 return, 1 33 32 ENDIF ELSE BEGIN … … 36 35 ; 37 36 ; 38 key_periodique = newgrid.key_periodique 37 key_periodic = newgrid.key_periodic 38 ; 39 @updateold 39 40 domdef 40 if newgrid.triangulation EQ 1 then triangles = triangule() ELSE triangles = -1 41 if newgrid.triangulation EQ 1 then triangles_list = triangule() $ 42 ELSE triangles_list = -1 43 ; 44 IF NOT keyword_set(key_forgetold) THEN BEGIN 45 @updateold 46 ENDIF 41 47 ; 42 48 return, 1 -
trunk/ToBeReviewed/GRILLE/cmpgrid.pro
r12 r13 14 14 ; we compare the structure which caracterise the grid whith 15 15 ; ccmeshparameters 16 ; 16 17 ; 17 18 case 1 of -
trunk/ToBeReviewed/GRILLE/domdef.pro
r12 r13 10 10 ; CATEGORY: 11 11 ; 12 ; CALLING SEQUENCE:domdef [,lon1, lon2, lat1, lat2[, prof1,prof2]] ou12 ; CALLING SEQUENCE:domdef [,lon1, lon2, lat1, lat2[,vert1,vert2]] ou 13 13 ; bien domdef,vecteur 14 14 ; 15 15 ; INPUTS:(facultatif), [vecteur a] 2, 4 ou 6 elements:; 16 16 ; sans l''utilisation des mots cles index,xindex,yindex,zindex: 17 ; * prof1, prof2: pour un domaine 3D dont la partie horizontale couvre tout17 ; *vert1, vert2: pour un domaine 3D dont la partie horizontale couvre tout 18 18 ; glam et gphi 19 19 ; *lon1, lon2, lat1, lat2: 20 20 ; definissant les longitudes min. max et les latitudes min, max du domaine a 21 21 ; etudier (tous les niveaux sont selectiones) 22 ; *lon1,lon2,lat1,lat2, prof1,prof2 pour specifier les profondeurs.22 ; *lon1,lon2,lat1,lat2,vert1,vert2 pour specifier les profondeurs. 23 23 ; 24 24 ; KEYWORD PARAMETERS: 25 ; 26 ; ENDPOINTS: a four elements vector [x1,y1,x2,y2] used to specify 27 ; that domdef must define the box used to make a plot (pltz, pltt, 28 ; plt1d) done strictly along the line (that can have any direction) 29 ; starting at (x1, y1) ending at (x2, y2). When defining endpoints, 30 ; you must also define TYPE which define the type of plots 31 ; ('pltz', 'xt', 'yt', 'zt', 'x', 'y', 'z', 't') will used 32 ; ENDPOINTS keywords 25 33 ; 26 34 ; FINDALWAYS:oblige a redefinir une boite meme qd auqun point … … 28 36 ; grille. 29 37 ; 30 ; GRI LLE:un string ou un vecteur de string contennant le nom des38 ; GRIDTYPE:un string ou un vecteur de string contennant le nom des 31 39 ; grilles (determinees uniquement par : 'T','U','V','W','F') pour 32 40 ; lesquelles le calcul doit etre fait. par … … 61 69 ; -nzt,w:entier qui contient le nombre de pts en profondeur de 62 70 ; la grille reduite au domaine 3D 63 ; - premierxt,u,v,f: le premierindice qui delimite le sous domaine71 ; -firstxt,u,v,f: le first indice qui delimite le sous domaine 64 72 ; suivant x 65 ; - premieryt,u,v,f: le premierindice qui delimite le sous domaine73 ; -firstyt,u,v,f: le first indice qui delimite le sous domaine 66 74 ; suivant y 67 ; - premierzt,w: le premierindice qui delimite le sous domaine75 ; -firstzt,w: le first indice qui delimite le sous domaine 68 76 ; suivant z 69 ; - dernierxt,u,v,f: le dernierindice qui delimite le sous domaine77 ; -lastxt,u,v,f: le last indice qui delimite le sous domaine 70 78 ; suivant x 71 ; - dernieryt,u,v,f: le dernierindice qui delimite le sous domaine79 ; -lastyt,u,v,f: le last indice qui delimite le sous domaine 72 80 ; suivant y 73 ; - dernierzt,w: le dernierindice qui delimite le sous domaine81 ; -lastzt,w: le last indice qui delimite le sous domaine 74 82 ; suivant z 75 83 ; … … 85 93 ; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr) 86 94 ; 8/2/98 95 ; rewrite everything, debug and spee-up Sebastien Masson April 2005 87 96 ;- 88 97 ;------------------------------------------------------------ 89 98 ;------------------------------------------------------------ 90 99 ;------------------------------------------------------------ 91 pro domdef,x1,x2,y1,y2,z1,z2, FINDALWAYS = findalways, GRILLE=grille, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 92 @common 93 tempsun = systime(1) ; pour key_performance 94 ;------------------------------------------------------------------- 95 ; determination en fonction du nombre d''arguments de 96 ; lon1,lon2,lat1,lat2,prof1,prof2 97 ;------------------------------------------------------------------- 98 if not keyword_set(grille) then grille=['T', 'U', 'V', 'W', 'F'] ELSE $ 99 for i = 0, n_elements(grille)-1 do grille[i] = strupcase(grille[i]) 100 if keyword_set(memeindices) then grille=['T', grille] 101 CASE N_PARAMS() OF 102 0: BEGIN 100 pro domdef, x1, x2, y1, y2, z1, z2, FINDALWAYS = findalways $ 101 , GRIDTYPE = gridtype, MEMEINDICES = memeindices $ 102 , XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex $ 103 , ENDPOINTS = endpoints, TYPE = type $ 104 , INDEX = index, _extra = ex 105 ;------------------------------------------------------------ 106 ; include commons 107 @cm_4mesh 108 IF NOT keyword_set(key_forgetold) THEN BEGIN 109 @updatenew 110 @updatekwd 111 ENDIF 112 ;--------------------- 113 tempsun = systime(1) ; pour key_performance 114 ; 115 CASE N_PARAMS() OF 116 0: 117 1: 118 2: 119 4: 120 6: 121 ELSE:BEGIN 122 ras = report('Bad number of parameter in the call of domdef.') 123 RETURN 124 END 125 ENDCASE 126 ; 127 IF keyword_set(endpoints) THEN BEGIN 128 IF NOT keyword_set(type) THEN BEGIN 129 dummy = report(['If domdef is used do find the box associated' $ 130 , 'to endpoints, you must also specify type keyword']) 131 return 132 ENDIF 133 CASE N_PARAMS() OF 134 0: 135 1:boxzoom = [x1] 136 2:boxzoom = [x1, x2] 137 4:boxzoom = [x1, x2, y1, y2] 138 6:boxzoom = [x1, x2, y1, y2, z1, z2] 139 ENDCASE 140 section, BOXZOOM = boxzoom, ENDPOINTS = endpoints, TYPE = type, /ONLYBOX 141 return 142 ENDIF 143 144 ;------------------------------------------------------------------- 145 ; recall domdef when there is only one input parameter 146 ;------------------------------------------------------------------- 147 ; 148 IF N_PARAMS() EQ 1 THEN BEGIN 149 CASE n_elements(x1) OF 150 2:domdef, x1[0], x1[1], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 151 4:domdef, x1[0], x1[1], x1[2], x1[3], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 152 6:domdef, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 153 ELSE:BEGIN 154 ras = report('Bad number of elements in x1') 155 RETURN 156 END 157 ENDCASE 158 RETURN 159 ENDIF 160 ;------------------------------------------------------------------- 161 ; default definitions and checks 162 ;------------------------------------------------------------------- 163 IF NOT keyword_set(gridtype) THEN gridtype = ['T', 'U', 'V', 'W', 'F'] $ 164 ELSE gridtype = strupcase(gridtype) 165 IF keyword_set(memeindices) THEN gridtype = ['T', gridtype] 166 ; default definitions 167 lon1t = 99999. & lon2t = -99999. & lat1t = 99999. & lat2t = -99999. 168 lon1u = 99999. & lon2u = -99999. & lat1u = 99999. & lat2u = -99999. 169 lon1v = 99999. & lon2v = -99999. & lat1v = 99999. & lat2v = -99999. 170 lon1f = 99999. & lon2f = -99999. & lat1f = 99999. & lat2f = -99999. 171 vert1t = 99999. & vert2t = -99999. & vert1w = 99999. & vert2w = -99999. 172 ; 173 IF N_PARAMS() EQ 2 THEN GOTO, vertical 174 ; 175 ;------------------------------------------------------------------- 176 ;------------------------------------------------------------------- 177 ; define all horizontal parameters ... 103 178 ; lon1 et lon2 104 test = 1105 IF (where(grille eq 'T'))[0] NE -1 THEN test = [test, glamt[*]]106 IF (where(grille eq 'W'))[0] NE -1 AND (where(grille eq 'T'))[0] EQ -1 THEN $107 test = [test, glamt[*]]108 IF (where(grille eq 'U'))[0] NE -1 AND n_elements(glamu) NE 0 THEN test = [test, glamu[*]]109 IF (where(grille eq 'V'))[0] NE -1 AND n_elements(glamv) NE 0 THEN test = [test, glamv[*]]110 IF (where(grille eq 'F'))[0] NE -1 AND n_elements(glamf) NE 0 THEN test = [test, glamf[*]]111 test = test[1:n_elements(test)-1]112 lon1=min([test], max = max)113 lon2=max114 179 ; lat1 et lat2 115 test = 1 116 IF (where(grille eq 'T'))[0] NE -1 THEN test = [test, gphit[*]] 117 IF (where(grille eq 'W'))[0] NE -1 AND (where(grille eq 'T'))[0] EQ -1 THEN $ 118 test = [test, gphit[*]] 119 IF (where(grille eq 'U'))[0] NE -1 AND n_elements(gphiu) NE 0 THEN test = [test, gphiu[*]] 120 IF (where(grille eq 'V'))[0] NE -1 AND n_elements(gphiv) NE 0 THEN test = [test, gphiv[*]] 121 IF (where(grille eq 'F'))[0] NE -1 AND n_elements(gphif) NE 0 THEN test = [test, gphif[*]] 122 test = test[1:n_elements(test)-1] 123 lat1=min([test], max = max) 124 lat2=max 125 ; prof1 et prof2 126 prof1=0 127 test = 1 128 IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN test = [test, gdept[jpk-1]] 129 IF (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN test = [test, gdepw[jpk-1]] 130 test = test[1:n_elements(test)-1] 131 prof2=max([test]) 132 end 133 1: BEGIN 134 xx1 = reform(x1) 135 taille = size(xx1) 136 CASE taille[1] OF 137 2:domdef, xx1[0], xx1[1] , FINDALWAYS = findalways,GRILLE=grille, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 138 4:domdef, xx1[0], xx1[1], xx1[2], xx1[3] , FINDALWAYS = findalways,GRILLE=grille, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 139 6:domdef, xx1[0], xx1[1], xx1[2], xx1[3], xx1[4], xx1[5], FINDALWAYS = findalways,GRILLE=grille, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 140 else: begin 141 ras = report('Mauvais nombre de parametre dans l''appel de DOMDEF') 142 return 143 end 144 ENDCASE 145 return 146 end 147 2: begin 148 ; lon1 et lon2 149 IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN 150 ouca = where(gdept ge prof1 and gdept le prof2,nzt) 151 if ouca[0] NE -1 then zt=gdept[ouca] ELSE zt = -1 152 ENDIF 153 if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then BEGIN 154 ouca = where(gdepw ge prof1 and gdepw le prof2,nzw) 155 if ouca[0] NE -1 then zw=gdepw[ouca] ELSE zw = -1 156 ENDIF 157 test = 1 158 IF (where(grille eq 'T'))[0] NE -1 THEN test = [test, glamt[*]] 159 IF (where(grille eq 'W'))[0] NE -1 AND (where(grille eq 'T'))[0] EQ -1 THEN $ 160 test = [test, glamt[*]] 161 IF (where(grille eq 'U'))[0] NE -1 AND n_elements(glamu) NE 0 THEN test = [test, glamu[*]] 162 IF (where(grille eq 'V'))[0] NE -1 AND n_elements(glamv) NE 0 THEN test = [test, glamv[*]] 163 IF (where(grille eq 'F'))[0] NE -1 AND n_elements(glamf) NE 0 THEN test = [test, glamf[*]] 164 test = test[1:n_elements(test)-1] 165 lon1=min([test], max = max) 166 lon2=max 167 ; lat1 et lat2 168 test = 1 169 IF (where(grille eq 'T'))[0] NE -1 THEN test = [test, gphit[*]] 170 IF (where(grille eq 'W'))[0] NE -1 AND (where(grille eq 'T'))[0] EQ -1 THEN $ 171 test = [test, gphit[*]] 172 IF (where(grille eq 'U'))[0] NE -1 AND n_elements(gphiu) NE 0 THEN test = [test, gphiu[*]] 173 IF (where(grille eq 'V'))[0] NE -1 AND n_elements(gphiv) NE 0 THEN test = [test, gphiv[*]] 174 IF (where(grille eq 'F'))[0] NE -1 AND n_elements(gphif) NE 0 THEN test = [test, gphif[*]] 175 test = test[1:n_elements(test)-1] 176 lat1=min([test], max = max) 177 lat2=max 178 ; prof1 et prof2 179 if keyword_set(zindex) OR keyword_set(index) then BEGIN 180 z1 = x1 181 z2 = x2 182 IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN prof1 = gdept[z1] ELSE prof1=gdepw[z1] 183 if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then prof2=gdepw[z2] ELSE prof2=gdept[z2] 184 ENDIF ELSE BEGIN 185 prof1=x1 186 prof2=x2 187 ENDELSE 188 end 189 4: begin 190 ; lon1 et lon2 191 if keyword_set(xindex) OR keyword_set(index) then BEGIN 192 if keyword_set(yindex) OR keyword_set(index) then BEGIN 193 if n_elements(glamf) NE 0 then BEGIN 194 lon1 = min([glamt[x1:x2, y1:y2], glamf[x1:x2, y1:y2]], max = lon2) 195 ENDIF ELSE BEGIN 196 lon1 = min(glamt[x1:x2, y1:y2], max = lon2) 197 ENDELSE 198 ENDIF ELSE BEGIN 199 lon1 = glamt[x1, 0] 200 if n_elements(glamf) NE 0 then lon2 = glamf[x2, 0] $ 201 ELSE lon2 = glamt[x2, 0] 180 ; firstx[tuvf], lastx[tuvf], nx[tuvf] 181 ;------------------------------------------------------------------- 182 ; check if the grid is defined for U and V points. If not, take care 183 ; of the cases gridtype eq 'U' or 'V' 184 ;------------------------------------------------------------------- 185 errstatus = 0 186 IF (finite(glamu[0]*gphiu[0]) EQ 0 OR n_elements(glamu) EQ 0 OR n_elements(gphiu) EQ 0) AND (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN 187 firstxu = !values.f_nan 188 lastxu = !values.f_nan 189 nxu = !values.f_nan 190 okgrid = where(gridtype NE 'U', count) 191 IF count NE 0 THEN gridtype = gridtype[okgrid] $ 192 ELSE errstatus = report('U grid is undefined. Impossible to call domdef with vargid = ''U''') 193 ENDIF 194 ; 195 IF (finite(glamv[0]*gphiv[0]) EQ 0 OR n_elements(glamv) EQ 0 OR n_elements(gphiv) EQ 0) AND (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN 196 firstxv = !values.f_nan 197 lastxv = !values.f_nan 198 nxv = !values.f_nan 199 okgrid = where(gridtype NE 'V', count) 200 IF count NE 0 THEN gridtype = gridtype[okgrid] $ 201 ELSE errstatus = report('V grid is undefined. Impossible to call domdef with vargid = ''V''') 202 ENDIF 203 IF errstatus EQ -1 THEN return 204 ;------------------------------------------------------------------- 205 ;------------------------------------------------------------------- 206 ; horizontal domain defined with lon1, lon2, lat1 and lat2 207 ;------------------------------------------------------------------- 208 ;------------------------------------------------------------------- 209 IF N_PARAMS() EQ 0 $ 210 OR ( (N_PARAMS() EQ 4 OR N_PARAMS() EQ 6) $ 211 AND NOT keyword_set(xindex) AND NOT keyword_set(yindex) AND NOT keyword_set(index) ) THEN BEGIN 212 IF N_PARAMS() EQ 0 THEN BEGIN 213 ; find lon1 and lon2 the longitudinal boudaries of the full domain 214 IF (where(gridtype eq 'T'))[0] NE -1 THEN lon1t = min(glamt, max = lon2t) 215 IF (where(gridtype eq 'W'))[0] NE -1 AND (where(gridtype eq 'T'))[0] EQ -1 THEN lon1t = min(glamt, max = lon2t) 216 IF (where(gridtype eq 'U'))[0] NE -1 THEN lon1u = min(glamu, max = lon2u) 217 IF (where(gridtype eq 'V'))[0] NE -1 THEN lon1v = min(glamv, max = lon2v) 218 IF (where(gridtype eq 'F'))[0] NE -1 THEN lon1f = min(glamf, max = lon2f) 219 lon1 = min([lon1t, lon1u, lon1v, lon1f]) 220 lon2 = max([lon2t, lon2u, lon2v, lon2f]) 221 ; find lat1 and lat2 the latitudinal boudaries of the full domain 222 IF (where(gridtype eq 'T'))[0] NE -1 THEN lat1t = min(gphit, max = lat2t) 223 IF (where(gridtype eq 'W'))[0] NE -1 AND (where(gridtype eq 'T'))[0] EQ -1 THEN lat1t = min(gphit, max = lat2t) 224 IF (where(gridtype eq 'U'))[0] NE -1 THEN lat1u = min(gphiu, max = lat2u) 225 IF (where(gridtype eq 'V'))[0] NE -1 THEN lat1v = min(gphiv, max = lat2v) 226 IF (where(gridtype eq 'F'))[0] NE -1 THEN lat1f = min(gphif, max = lat2f) 227 lat1 = min([lat1t, lat1u, lat1v, lat1f]) 228 lat2 = max([lat2t, lat2u, lat2v, lat2f]) 229 ENDIF ELSE BEGIN 230 lon1 = min([x1, x2], max = lon2) 231 lat1 = min([y1, y2], max = lat2) 232 ENDELSE 233 ; find firstxt, firstxt, nxt and nyt according to lon1, lon2, lat1 and lat2 234 IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN 235 dom = where( (glamt GE lon1) AND (glamt LE lon2) $ 236 AND (gphit GE lat1) AND (gphit LE lat2) ) 237 IF (dom[0] EQ -1) THEN BEGIN 238 IF keyword_set(findalways) THEN BEGIN 239 domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 240 IF N_PARAMS() EQ 6 THEN $ 241 domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 242 RETURN 243 ENDIF ELSE BEGIN 244 ras = report('WARNING, The box does not contain any T points...') 245 firstxt = -1 & lastxt = -1 & nxt = 0 246 firstyt = -1 & lastyt = -1 & nyt = 0 247 ENDELSE 248 ENDIF ELSE BEGIN 249 jyt = dom / jpi 250 ixt = temporary(dom) MOD jpi 251 firstxt = min(temporary(ixt), max = lastxt) 252 firstyt = min(temporary(jyt), max = lastyt) 253 nxt = lastxt - firstxt + 1 254 nyt = lastyt - firstyt + 1 255 ENDELSE 256 ENDIF 257 ; find firstxu, firstxu, firstyu, firstyu, nxu and nyu 258 ; according to lon1, lon2, lat1 and lat2 259 IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN 260 IF keyword_set(memeindices) THEN BEGIN 261 firstxu = firstxt & lastxu = lastxt & nxu = nxt 262 firstyu = firstyt & lastyu = lastyt & nyu = nyt 263 ENDIF ELSE BEGIN 264 dom = where( (glamu GE lon1) AND (glamu LE lon2) $ 265 AND (gphiu GE lat1) AND (gphiu LE lat2) ) 266 IF (dom[0] EQ -1) THEN BEGIN 267 IF keyword_set(findalways) THEN BEGIN 268 domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 269 IF N_PARAMS() EQ 6 THEN $ 270 domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 271 RETURN 272 ENDIF ELSE BEGIN 273 ras = report('WARNING, The box does not contain any U points...') 274 firstxu = -1 & lastxu = -1 & nxu = 0 275 firstyu = -1 & lastyu = -1 & nyu = 0 276 ENDELSE 277 ENDIF ELSE BEGIN 278 jyu = dom / jpi 279 ixu = temporary(dom) MOD jpi 280 firstxu = min(temporary(ixu), max = lastxu) 281 firstyu = min(temporary(jyu), max = lastyu) 282 nxu = lastxu - firstxu + 1 283 nyu = lastyu - firstyu + 1 284 ENDELSE 285 ENDELSE 286 ENDIF 287 ; find firstxv, firstxv, firstyv, firstyv, nxv and nyv 288 ; according to lon1, lon2, lat1 and lat2 289 IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN 290 IF keyword_set(memeindices) THEN BEGIN 291 firstxv = firstxt & lastxv = lastxt & nxv = nxt 292 firstyv = firstyt & lastyv = lastyt & nyv = nyt 293 ENDIF ELSE BEGIN 294 dom = where( (glamv GE lon1) AND (glamv LE lon2) $ 295 AND (gphiv GE lat1) AND (gphiv LE lat2) ) 296 IF (dom[0] EQ -1) THEN BEGIN 297 IF keyword_set(findalways) THEN BEGIN 298 domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 299 IF N_PARAMS() EQ 6 THEN $ 300 domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 301 RETURN 302 ENDIF ELSE BEGIN 303 ras = report('WARNING, The box does not contain any V points...') 304 firstxv = -1 & lastxv = -1 & nxv = 0 305 firstyv = -1 & lastyv = -1 & nyv = 0 306 ENDELSE 307 ENDIF ELSE BEGIN 308 jyv = dom / jpi 309 ixv = temporary(dom) MOD jpi 310 firstxv = min(temporary(ixv), max = lastxv) 311 firstyv = min(temporary(jyv), max = lastyv) 312 nxv = lastxv - firstxv + 1 313 nyv = lastyv - firstyv + 1 314 ENDELSE 315 ENDELSE 316 ENDIF 317 ; find firstxf, firstxf, firstyf, firstyf, nxf and nyf 318 ; according to lon1, lon2, lat1 and lat2 319 IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN 320 IF keyword_set(memeindices) THEN BEGIN 321 firstxf = firstxt & lastxf = lastxt & nxf = nxt 322 firstyf = firstyt & lastyf = lastyt & nyf = nyt 323 ENDIF ELSE BEGIN 324 dom = where( (glamf GE lon1) AND (glamf LE lon2) $ 325 AND (gphif GE lat1) AND (gphif LE lat2) ) 326 IF (dom[0] EQ -1) THEN BEGIN 327 IF keyword_set(findalways) THEN BEGIN 328 domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 329 IF N_PARAMS() EQ 6 THEN $ 330 domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 331 RETURN 332 ENDIF ELSE BEGIN 333 ras = report('WARNING, The box does not contain any F points...') 334 firstxf = -1 & lastxf = -1 & nxf = 0 335 firstyf = -1 & lastyf = -1 & nyf = 0 336 ENDELSE 337 ENDIF ELSE BEGIN 338 jyf = dom / jpi 339 ixf = temporary(dom) MOD jpi 340 firstxf = min(temporary(ixf), max = lastxf) 341 firstyf = min(temporary(jyf), max = lastyf) 342 nxf = lastxf - firstxf + 1 343 nyf = lastyf - firstyf + 1 344 ENDELSE 345 ENDELSE 346 ENDIF 347 ; 348 ENDIF ELSE BEGIN 349 CASE 1 OF 350 ;------------------------------------------------------------------- 351 ;------------------------------------------------------------------- 352 ; horizontal domain defined with the X and Y indexes 353 ;------------------------------------------------------------------- 354 ;------------------------------------------------------------------- 355 (keyword_set(xindex) AND keyword_set(yindex)) OR keyword_set(index):BEGIN 356 fstx = min([x1, x2], max = lstx) 357 fsty = min([y1, y2], max = lsty) 358 IF fstx LT 0 OR lstx GE jpi THEN BEGIN 359 ras = report('Bad definition of X1 or X2') 360 return 361 ENDIF 362 IF fsty LT 0 OR lsty GE jpj THEN BEGIN 363 ras = report('Bad definition of Y1 or Y2') 364 return 365 ENDIF 366 nx = lstx - fstx + 1 367 ny = lsty - fsty + 1 368 ; find lon1t, lon2t, lat1t, lat2t, firstxt, firstxt, nxt and nyt 369 ; according to x1, x2, y1, y2 370 IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype eq 'W'))[0] NE -1 THEN BEGIN 371 lon1t = min(glamt[fstx:lstx, fsty:lsty], max = lon2t) 372 lat1t = min(gphit[fstx:lstx, fsty:lsty], max = lat2t) 373 firstxt = fstx & lastxt = lstx 374 firstyt = fsty & lastyt = lsty 375 nxt = nx & nyt = ny 376 ENDIF 377 ; find lon1u, lon2u, lat1u, lat2u, firstxu, firstxu, nxu and nyu 378 ; according to x1, x2, y1, y2 379 IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN 380 lon1u = min(glamu[fstx:lstx, fsty:lsty], max = lon2u) 381 lat1u = min(gphiu[fstx:lstx, fsty:lsty], max = lat2u) 382 firstxu = fstx & lastxu = lstx 383 firstyu = fsty & lastyu = lsty 384 nxu = nx & nyu = ny 385 ENDIF 386 ; find lon1v, lon2v, lat1v, lat2v, firstxv, firstxv, nxv and nyv 387 ; according to x1, x2, y1, y2 388 IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN 389 lon1v = min(glamv[fstx:lstx, fsty:lsty], max = lon2v) 390 lat1v = min(gphiv[fstx:lstx, fsty:lsty], max = lat2v) 391 firstxv = fstx & lastxv = lstx 392 firstyv = fsty & lastyv = lsty 393 nxv = nx & nyv = ny 394 ENDIF 395 ; find lon1f, lon2f, lat1f, lat2f, firstxf, firstxf, nxf and nyf 396 ; according to x1, x2, y1, y2 397 IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN 398 lon1f = min(glamf[fstx:lstx, fsty:lsty], max = lon2f) 399 lat1f = min(gphif[fstx:lstx, fsty:lsty], max = lat2f) 400 firstxf = fstx & lastxf = lstx 401 firstyf = fsty & lastyf = lsty 402 nxf = nx & nyf = ny 403 ENDIF 404 lon1 = min([lon1t, lon1u, lon1v, lon1f]) 405 lon2 = max([lon2t, lon2u, lon2v, lon2f]) 406 lat1 = min([lat1t, lat1u, lat1v, lat1f]) 407 lat2 = max([lat2t, lat2u, lat2v, lat2f]) 408 END 409 ;------------------------------------------------------------------- 410 ;------------------------------------------------------------------- 411 ; horizontal domain defined with the X index and lat1/lat2 412 ;------------------------------------------------------------------- 413 ;------------------------------------------------------------------- 414 keyword_set(xindex):BEGIN 415 fstx = min([x1, x2], max = lstx) 416 IF fstx LT 0 OR lstx GE jpi THEN BEGIN 417 ras = report('Bad definition of X1 or X2') 418 return 419 ENDIF 420 nx = lstx - fstx + 1 421 lat1 = min([y1, y2], max = lat2) 422 ; find lon1t, lon2t, firstxt, firstxt, firstyt, firstyt, nxt 423 ; and nyt according to x1, x2, lat1 and lat2 424 IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN 425 firstxt = fstx & lastxt = lstx & nxt = nx 426 dom = where( (gphit[fstx:lstx, *] GE lat1) AND (gphit[fstx:lstx, *] LE lat2) ) 427 IF (dom[0] EQ -1) THEN BEGIN 428 IF keyword_set(findalways) THEN BEGIN 429 domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 430 IF N_PARAMS() EQ 6 THEN $ 431 domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 432 RETURN 433 ENDIF ELSE BEGIN 434 ras = report('WARNING, The box does not contain any T points...') 435 firstyt = -1 & lastyt = -1 & nyt = 0 202 436 ENDELSE 203 ENDIF ELSE BEGIN 204 lon1=x1 205 lon2=x2 206 ENDELSE 207 ; lat1 et lat2 208 if keyword_set(yindex) OR keyword_set(index) then BEGIN 209 if keyword_set(xindex) OR keyword_set(index) then begin 210 lat1 = min([gphit[x1, y1], gphit[x2, y1]]) 211 if n_elements(gphif) NE 0 then lat2 = max([gphif[x1, y2], gphif[x2, y2]]) $ 212 ELSE lat2 = max([gphit[x1, y2], gphit[x2, y2]]) 213 ENDIF ELSE BEGIN 214 lat1 = gphit[0, y1] 215 if n_elements(gphif) NE 0 then lat2 = gphif[0, y2] $ 216 ELSE lat2 = gphit[0, y2] 437 ENDIF ELSE BEGIN 438 jyt = temporary(dom) / nx 439 firstyt = min(temporary(jyt), max = lastyt) 440 nyt = lastyt - firstyt + 1 441 ENDELSE 442 IF nyt NE 0 THEN lon1t = min(glamt[firstxt:lastxt, firstyt:lastyt], max = lon2t) 443 ENDIF 444 ; find lon1u, lon2u, firstxu, firstxu, firstyu, firstyu, nxu 445 ; and nyu according to x1, x2, lat1 and lat2 446 IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN 447 firstxu = fstx & lastxu = lstx & nxu = nx 448 IF keyword_set(memeindices) THEN BEGIN 449 firstyu = firstyt & lastyu = lastyt & nyu = nyt 450 ENDIF ELSE BEGIN 451 dom = where( (gphiu[fstx:lstx, *] GE lat1) AND (gphiu[fstx:lstx, *] LE lat2) ) 452 IF (dom[0] EQ -1) THEN BEGIN 453 IF keyword_set(findalways) THEN BEGIN 454 domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 455 IF N_PARAMS() EQ 6 THEN $ 456 domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 457 RETURN 458 ENDIF ELSE BEGIN 459 ras = report('WARNING, The box does not contain any U points...') 460 firstyu = -1 & lastyu = -1 & nyu = 0 461 ENDELSE 462 ENDIF ELSE BEGIN 463 jyu = temporary(dom) / nx 464 firstyu = min(temporary(jyu), max = lastyu) 465 nyu = lastyu - firstyu + 1 217 466 ENDELSE 218 ENDIF ELSE BEGIN 219 lat1=y1 220 lat2=y2 221 ENDELSE 222 ; prof1 et prof2 223 prof1=0 224 test = 1 225 IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN test = [test, gdept[jpk-1]] 226 IF (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN test = [test, gdepw[jpk-1]] 227 test = test[1:n_elements(test)-1] 228 prof2=max([test]) 229 end 230 6: begin 231 ; lon1 et lon2 232 if keyword_set(xindex) OR keyword_set(index) then BEGIN 233 if keyword_set(yindex) OR keyword_set(index) then begin 234 lon1 = min([glamt[x1, y1], glamt[x1, y2]]) 235 if n_elements(glamf) NE 0 then lon2 = max([glamf[x2, y1], glamf[x2, y2]]) $ 236 ELSE lon2 = max([glamt[x2, y1], glamt[x2, y2]]) 237 ENDIF ELSE BEGIN 238 lon1 = glamt[x1, 0] 239 if n_elements(glamf) NE 0 then lon2 = glamf[x2, 0] $ 240 ELSE lon2 = glamt[x2, 0] 467 ENDELSE 468 IF nyu NE 0 THEN lon1u = min(glamu[firstxu:lastxu, firstyu:lastyu], max = lon2u) 469 ENDIF 470 ; find lon1v, lon2v, firstxv, firstxv, firstyv, firstyv, nxv 471 ; and nyv according to x1, x2, lat1 and lat2 472 IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN 473 firstxv = fstx & lastxv = lstx & nxv = nx 474 IF keyword_set(memeindices) THEN BEGIN 475 firstyv = firstyt & lastyv = lastyt & nyv = nyt 476 ENDIF ELSE BEGIN 477 dom = where( (gphiv[fstx:lstx, *] GE lat1) AND (gphiv[fstx:lstx, *] LE lat2) ) 478 IF (dom[0] EQ -1) THEN BEGIN 479 IF keyword_set(findalways) THEN BEGIN 480 domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 481 IF N_PARAMS() EQ 6 THEN $ 482 domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 483 RETURN 484 ENDIF ELSE BEGIN 485 ras = report('WARNING, The box does not contain any V points...') 486 firstyv = -1 & lastyv = -1 & nyv = 0 487 ENDELSE 488 ENDIF ELSE BEGIN 489 jyv = temporary(dom) / nx 490 firstyv = min(temporary(jyv), max = lastyv) 491 nyv = lastyv - firstyv + 1 241 492 ENDELSE 242 ENDIF ELSE BEGIN 243 lon1=x1 244 lon2=x2 245 ENDELSE 246 ; lat1 et lat2 247 if keyword_set(yindex) OR keyword_set(index) then BEGIN 248 if keyword_set(xindex) OR keyword_set(index) then begin 249 lat1 = min([gphit[x1, y1], gphit[x2, y1]]) 250 if n_elements(gphif) NE 0 then lat2 = max([gphif[x1, y2], gphif[x2, y2]]) $ 251 ELSE lat2 = max([gphit[x1, y2], gphit[x2, y2]]) 252 ENDIF ELSE BEGIN 253 lat1 = gphit[0, y1] 254 if n_elements(gphif) NE 0 then lat2 = gphif[0, y2] $ 255 ELSE lat2 = gphit[0, y2] 493 ENDELSE 494 IF nyv NE 0 THEN lon1v = min(glamv[firstxv:lastxv, firstyv:lastyv], max = lon2v) 495 ENDIF 496 ; find lon1f, lon2f, firstxf, firstxf, firstyf, firstyf, nxf 497 ; and nyf according to x1, x2, lat1 and lat2 498 IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN 499 firstxf = fstx & lastxf = lstx & nxf = nx 500 IF keyword_set(memeindices) THEN BEGIN 501 firstyf = firstyt & lastyf = lastyt & nyf = nyt 502 ENDIF ELSE BEGIN 503 dom = where( (gphif[fstx:lstx, *] GE lat1) AND (gphif[fstx:lstx, *] LE lat2) ) 504 IF (dom[0] EQ -1) THEN BEGIN 505 IF keyword_set(findalways) THEN BEGIN 506 domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 507 IF N_PARAMS() EQ 6 THEN $ 508 domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 509 RETURN 510 ENDIF ELSE BEGIN 511 ras = report('WARNING, The box does not contain any F points...') 512 firstyf = -1 & lastyf = -1 & nyf = 0 513 ENDELSE 514 ENDIF ELSE BEGIN 515 jyf = temporary(dom) / nx 516 firstyf = min(temporary(jyf), max = lastyf) 517 nyf = lastyf - firstyf + 1 256 518 ENDELSE 257 ENDIF ELSE BEGIN 258 lat1=y1 259 lat2=y2 260 ENDELSE 261 ; prof1 et prof2 262 if keyword_set(zindex) OR keyword_set(index) then begin 263 IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN prof1=gdept[z1] ELSE prof1=gdepw[z1] 264 if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then prof2=gdepw[z2] ELSE prof2=gdept[z2] 265 ENDIF ELSE BEGIN 266 prof1=z1 267 prof2=z2 268 ENDELSE 269 end 270 else: begin 271 ras = report('Mauvais nombre de parametre dans l''appel de DOMDEF') 272 return 273 end 274 endcase 275 ;---------------------------------------------------------------------- 276 ; determination de nzt et nwt 277 ;---------------------------------------------------------------------- 278 IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN 279 if keyword_set(zindex) OR keyword_set(index) then nzt = z2-z1+1 $ 280 ELSE BEGIN 281 ouca = where(gdept ge prof1 and gdept le prof2,nzt) 282 if ouca[0] NE -1 then zt=gdept[ouca] ELSE zt = -1 283 ENDELSE 284 ENDIF 285 if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then BEGIN 286 if keyword_set(zindex) OR keyword_set(index) then nzw = z2-z1+1 $ 287 ELSE BEGIN 288 ouca = where(gdepw ge prof1 and gdepw le prof2,nzw) 289 if ouca[0] NE -1 then zw=gdepw[ouca] ELSE zw = -1 290 ENDELSE 291 ENDIF 292 ;---------------------------------------------------------------------- 293 ; determination de premierzt, dernierzt, premierzw, dernierzw 294 ;---------------------------------------------------------------------- 295 IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN 296 if keyword_set(zindex) OR keyword_set(index) then BEGIN 297 premierzt = z1 298 dernierzt = z2 519 ENDELSE 520 IF nyf NE 0 THEN lon1f = min(glamf[firstxf:lastxf, firstyf:lastyf], max = lon2f) 521 ENDIF 522 lon1 = min([lon1t, lon1u, lon1v, lon1f]) 523 lon2 = max([lon2t, lon2u, lon2v, lon2f]) 524 END 525 ;------------------------------------------------------------------- 526 ;------------------------------------------------------------------- 527 ; horizontal domain defined with the Y index and lon1/lon2 528 ;------------------------------------------------------------------- 529 ;------------------------------------------------------------------- 530 keyword_set(yindex):BEGIN 531 fsty = min([y1, y2], max = lsty) 532 IF fsty LT 0 OR lsty GE jpj THEN BEGIN 533 ras = report('Bad definition of Y1 or Y2') 534 return 535 ENDIF 536 ny = lsty - fsty + 1 537 lon1 = min([x1, x2], max = lon2) 538 ; find lat1t, lat2t, firstxt, firstxt, firstyt, firstyt, nxt 539 ; and nyt according to x1, x2, lon1 and lon2 540 IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN 541 firstyt = fsty & lastyt = lsty & nyt = ny 542 dom = where( (glamt[*, fsty:lsty] GE lon1) AND (glamt[*, fsty:lsty] LE lon2) ) 543 IF (dom[0] EQ -1) THEN BEGIN 544 IF keyword_set(findalways) THEN BEGIN 545 domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 546 IF N_PARAMS() EQ 6 THEN $ 547 domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 548 RETURN 549 ENDIF ELSE BEGIN 550 ras = report('WARNING, The box does not contain any T points...') 551 firstxt = -1 & lastxt = -1 & nxt = 0 552 ENDELSE 553 ENDIF ELSE BEGIN 554 jxt = temporary(dom) MOD jpi 555 firstxt = min(temporary(jxt), max = lastxt) 556 nxt = lastxt - firstxt + 1 557 ENDELSE 558 IF nxt NE 0 THEN lat1t = min(gphit[firstxt:lastxt, firstyt:lastyt], max = lat2t) 559 ENDIF 560 ; find lat1u, lat2u, firstxu, firstxu, firstyu, firstyu, nxu 561 ; and nyu according to x1, x2, lon1 and lon2 562 IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN 563 firstyu = fsty & lastyu = lsty & nyu = ny 564 IF keyword_set(memeindices) THEN BEGIN 565 firstxu = firstyt & lastxu = lastyt & nxu = nxt 566 ENDIF ELSE BEGIN 567 dom = where( (glamu[*, fsty:lsty] GE lon1) AND (glamu[*, fsty:lsty] LE lon2) ) 568 IF (dom[0] EQ -1) THEN BEGIN 569 IF keyword_set(findalways) THEN BEGIN 570 domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 571 IF N_PARAMS() EQ 6 THEN $ 572 domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 573 RETURN 574 ENDIF ELSE BEGIN 575 ras = report('WARNING, The box does not contain any U points...') 576 firstxu = -1 & lastxu = -1 & nxu = 0 577 ENDELSE 578 ENDIF ELSE BEGIN 579 jxu = temporary(dom) MOD jpi 580 firstxu = min(temporary(jxu), max = lastxu) 581 nxu = lastxu - firstxu + 1 582 ENDELSE 583 ENDELSE 584 IF nxu NE 0 THEN lat1u = min(gphiu[firstxu:lastxu, firstyu:lastyu], max = lat2u) 585 ENDIF 586 ; find lat1v, lat2v, firstxv, firstxv, firstyv, firstyv, nxv 587 ; and nyv according to x1, x2, lon1 and lon2 588 IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN 589 firstyv = fsty & lastyv = lsty & nyv = ny 590 IF keyword_set(memeindices) THEN BEGIN 591 firstxv = firstyt & lastxv = lastyt & nxv = nxt 592 ENDIF ELSE BEGIN 593 dom = where( (glamv[*, fsty:lsty] GE lon1) AND (glamv[*, fsty:lsty] LE lon2) ) 594 IF (dom[0] EQ -1) THEN BEGIN 595 IF keyword_set(findalways) THEN BEGIN 596 domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 597 IF N_PARAMS() EQ 6 THEN $ 598 domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 599 RETURN 600 ENDIF ELSE BEGIN 601 ras = report('WARNING, The box does not contain any V points...') 602 firstxv = -1 & lastxv = -1 & nxv = 0 603 ENDELSE 604 ENDIF ELSE BEGIN 605 jxv = temporary(dom) MOD jpi 606 firstxv = min(temporary(jxv), max = lastxv) 607 nxv = lastxv - firstxv + 1 608 ENDELSE 609 ENDELSE 610 IF nxv NE 0 THEN lat1v = min(gphiv[firstxv:lastxv, firstyv:lastyv], max = lat2v) 611 ENDIF 612 ; find lat1f, lat2f, firstxf, firstxf, firstyf, firstyf, nxf 613 ; and nyf according to x1, x2, lon1 and lon2 614 IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN 615 firstyf = fsty & lastyf = lsty & nyf = ny 616 IF keyword_set(memeindices) THEN BEGIN 617 firstxf = firstyt & lastxf = lastyt & nxf = nxt 618 ENDIF ELSE BEGIN 619 dom = where( (glamf[*, fsty:lsty] GE lon1) AND (glamf[*, fsty:lsty] LE lon2) ) 620 IF (dom[0] EQ -1) THEN BEGIN 621 IF keyword_set(findalways) THEN BEGIN 622 domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 623 IF N_PARAMS() EQ 6 THEN $ 624 domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 625 RETURN 626 ENDIF ELSE BEGIN 627 ras = report('WARNING, The box does not contain any F points...') 628 firstxf = -1 & lastyf = -1 & nxf = 0 629 ENDELSE 630 ENDIF ELSE BEGIN 631 jxf = temporary(dom) MOD jpi 632 firstxf = min(temporary(jxf), max = lastxf) 633 nxf = lastxf - firstxf + 1 634 ENDELSE 635 ENDELSE 636 IF nxf NE 0 THEN lat1f = min(gphif[firstxf:lastxf, firstyf:lastyf], max = lat2f) 637 ENDIF 638 lat1 = min([lat1t, lat1u, lat1v, lat1f]) 639 lat2 = max([lat2t, lat2u, lat2v, lat2f]) 640 END 641 ENDCASE 642 ENDELSE 643 ;------------------------------------------------------------------- 644 ; The extracted domain is it regular or not? 645 ;------------------------------------------------------------------- 646 CASE 1 OF 647 ((where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype eq 'W'))[0] NE -1) AND nxt NE 0 AND nyt NE 0:BEGIN 648 ; to get faster, we first test the most basic cases befor testing the 649 ; full array. 650 CASE 0 OF 651 array_equal(glamt[firstxt:lastxt, firstyt], glamt[firstxt:lastxt, lastyt]):key_irregular = 1b 652 array_equal(gphit[firstxt, firstyt:lastyt], gphit[lastxt, firstyt:lastyt]):key_irregular = 1b 653 array_equal(glamt[firstxt:lastxt, firstyt:lastyt] $ 654 , glamt[firstxt:lastxt, firstyt]#replicate(1, nyt)):key_irregular = 1b 655 array_equal(gphit[firstxt:lastxt, firstyt:lastyt] $ 656 , replicate(1, nxt)#(gphit[firstxt, firstyt:lastyt])[*]):key_irregular = 1b 657 ELSE:key_irregular = 0b 658 ENDCASE 659 END 660 (where(gridtype eq 'U'))[0] NE -1 AND nxu NE 0 AND nyu NE 0:BEGIN 661 CASE 0 OF 662 array_equal(glamu[firstxu:lastxu, firstyu], glamu[firstxu:lastxu, lastyu]):key_irregular = 1b 663 array_equal(gphiu[firstxu, firstyu:lastyu], gphiu[lastxu, firstyu:lastyu]):key_irregular = 1b 664 array_equal(glamu[firstxu:lastxu, firstyu:lastyu] $ 665 , glamu[firstxu:lastxu, firstyu]#replicate(1, nyu)):key_irregular = 1b 666 array_equal(gphiu[firstxu:lastxu, firstyu:lastyu] $ 667 , replicate(1, nxu)#(gphiu[firstxu, firstyu:lastyu])[*]):key_irregular = 1b 668 ELSE:key_irregular = 0b 669 ENDCASE 670 END 671 (where(gridtype eq 'V'))[0] NE -1 AND nxv NE 0 AND nyv NE 0:BEGIN 672 CASE 0 OF 673 array_equal(glamv[firstxv:lastxv, firstyv], glamv[firstxv:lastxv, lastyv]):key_irregular = 1b 674 array_equal(gphiv[firstxv, firstyv:lastyv], gphiv[lastxv, firstyv:lastyv]):key_irregular = 1b 675 array_equal(glamv[firstxv:lastxv, firstyv:lastyv] $ 676 , glamv[firstxv:lastxv, firstyv]#replicate(1, nyv)):key_irregular = 1b 677 array_equal(gphiv[firstxv:lastxv, firstyv:lastyv] $ 678 , replicate(1, nxv)#(gphiv[firstxv, firstyv:lastyv])[*]):key_irregular = 1b 679 ELSE:key_irregular = 0b 680 ENDCASE 681 END 682 (where(gridtype eq 'F'))[0] NE -1 AND nxf NE 0 AND nyf NE 0:BEGIN 683 CASE 0 OF 684 array_equal(glamf[firstxf:lastxf, firstyf], glamf[firstxf:lastxf, lastyf]):key_irregular = 1b 685 array_equal(gphif[firstxf, firstyf:lastyf], gphif[lastxf, firstyf:lastyf]):key_irregular = 1b 686 array_equal(glamf[firstxf:lastxf, firstyf:lastyf] $ 687 , glamf[firstxf:lastxf, firstyf]#replicate(1, nyf)):key_irregular = 1b 688 array_equal(gphif[firstxf:lastxf, firstyf:lastyf] $ 689 , replicate(1, nxf)#(gphif[firstxf, firstyf:lastyf])[*]):key_irregular = 1b 690 ELSE:key_irregular = 0b 691 ENDCASE 692 END 693 ELSE: 694 ENDCASE 695 ; 696 ;------------------------------------------------------------------- 697 ;------------------------------------------------------------------- 698 ; define all vertical parameters ... 699 ; vert1, vert2 700 ; firstz[tw], lastz[tw], nz[tw] 701 ;------------------------------------------------------------------- 702 ;------------------------------------------------------------------- 703 ; 704 vertical: 705 ; 706 ;------------------------------------------------------------------- 707 ; vertical domain defined with vert1, vert2 708 ;------------------------------------------------------------------- 709 IF NOT (keyword_set(zindex) OR keyword_set(index)) THEN BEGIN 710 ; define vert1 et vert2 711 CASE N_PARAMS() OF 712 2:vert1 = min([x1, x2], max = vert2) 713 6:vert1 = min([z1, z2], max = vert2) 714 ELSE:BEGIN 715 IF (inter(byte(gridtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN $ 716 vert1t = min(gdept, max = vert2t) 717 IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN $ 718 vert1w = min(gdepw, max = vert2w) 719 vert1 = min([vert1t, vert1w]) 720 vert2 = max([vert2t, vert2w]) 721 END 722 ENDCASE 723 ; define firstzt, firstzt, nzt 724 IF (inter(byte(gridtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN 725 domz = where(gdept ge vert1 and gdept le vert2, nzt) 726 IF nzt NE 0 THEN BEGIN 727 firstzt = domz[0] 728 lastzt = domz[nzt-1] 299 729 ENDIF ELSE BEGIN 300 if zt[0] NE -1 then begin 301 premierzt = (where(gdept eq zt[0]))[0] 302 dernierzt = (where(gdept eq zt[nzt-1]))[0] 303 ENDIF ELSE BEGIN 304 premierzt = -1 305 dernierzt = -1 306 ENDELSE 307 ENDELSE 308 endif 309 if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then begin 310 if keyword_set(zindex) OR keyword_set(index) then BEGIN 311 premierzw = z1 312 dernierzw = z2 730 ras = report('WARNING, The box does not contain any T level...') 731 firstzt = -1 732 lastzt = -1 733 ENDELSE 734 ENDIF 735 ; define firstzw, firstzw, nzw 736 IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN 737 IF keyword_set(memeindices) THEN BEGIN 738 firstzw = firstzt 739 lastzw = lastzt 740 nzw = nzt 313 741 ENDIF ELSE BEGIN 314 if zw[0] NE -1 then begin 315 premierzw = (where(gdepw eq zw[0]))[0] 316 dernierzw = (where(gdepw eq zw[nzw-1]))[0] 317 ENDIF ELSE BEGIN 318 premierzw = -1 319 dernierzw = -1 320 ENDELSE 321 ENDELSE 322 endif 323 ;------------------------------------------------------------------- 324 ;------------------------------------------------------------------- 325 ; calcul pour les pts t et w 326 ;------------------------------------------------------------------- 327 ;------------------------------------------------------------------- 328 tempdeux = systime(1) ; pour key_performance =2 329 if (where(grille eq 'T'))[0] NE -1 or (where(grille EQ 'W'))[0] NE -1 then begin 330 if (keyword_set(xindex) AND keyword_set(yindex)) OR keyword_set(index) then begin 331 premierxt = x1 332 premieryt = y1 333 dernierxt = x2 334 dernieryt = y2 335 ENDIF ELSE BEGIN 336 dom = where( (glamt ge lon1) $ 337 and (glamt le lon2) $ 338 and (gphit ge lat1) $ 339 and (gphit le lat2) ) 340 if (dom[0] eq -1) then BEGIN 341 if keyword_set(findalways) then BEGIN 342 domdef, prof1,prof2, GRILLE=grille, MEMEINDICES = memeindices, _extra = ex 343 return 344 ENDIF ELSE BEGIN 345 ras = report('ATTENTION, la boite ne contient pas de points t') 346 goto, ptu 347 ENDELSE 348 endif 349 IF testvar(var = key_performance) EQ 2 THEN $ 350 print, 'temps domdef: trouver le sous domaine', systime(1)-tempdeux 351 ;------------------------------------------------------------------- 352 ; on recupere les abscisses et les ordonees des points selectionnes 353 ;------------------------------------------------------------------- 354 tempdeux = systime(1) ; pour key_performance =2 355 ordonnee = dom/jpi 356 abscisse = dom-ordonnee*jpi 357 ;------------------------------------------------------------------- 358 ; on les ordonnes et on elimine les elements en double 359 ;------------------------------------------------------------------- 360 ordonnee = ordonnee 361 abscisse = abscisse 362 ordonnee = ordonnee(uniq(ordonnee, sort(ordonnee))) 363 abscisse = abscisse(uniq(abscisse, sort(abscisse))) 364 ;------------------------------------------------------------------- 365 ; indices des bornes du domaine 366 ;------------------------------------------------------------------- 367 premierxt = abscisse[0] 368 premieryt = ordonnee[0] 369 dernierxt = abscisse[n_elements(abscisse) -1] 370 dernieryt = ordonnee[n_elements(ordonnee) -1] 742 domz = where(gdepw ge vert1 and gdepw le vert2, nzw) 743 IF nzw NE 0 THEN BEGIN 744 firstzw = domz[0] 745 lastzw = domz[nzw-1] 746 ENDIF ELSE BEGIN 747 ras = report('WARNING, The box does not contain any W level...') 748 firstzw = -1 749 lastzw = -1 750 ENDELSE 371 751 ENDELSE 372 ;------------------------------------------------------------------- 373 ; nombre de points selectionnes 374 ;------------------------------------------------------------------- 375 nxt = dernierxt-premierxt+1 376 nyt = dernieryt-premieryt+1 377 IF testvar(var = key_performance) EQ 2 THEN $ 378 print, 'temps domdef: calculs de premier, dernier,...', systime(1)-tempdeux 379 ;------------------------------------------------------------------- 380 ; la grille est reguliere ou non ?? 381 ;------------------------------------------------------------------- 382 tempdeux = systime(1) ; pour key_performance =2 383 if (total(glamt[premierxt:dernierxt,premieryt:dernieryt] $ 384 NE glamt[premierxt:dernierxt,0]#replicate(1, nyt)) eq 0 $ 385 AND total(gphit[premierxt:dernierxt,premieryt:dernieryt] $ 386 NE replicate(1, nxt)#(gphit[0,premieryt:dernieryt])[*]) eq 0 ) $ 387 then key_irregular = 0 ELSE key_irregular = 1 388 IF testvar(var = key_performance) EQ 2 THEN $ 389 print, 'temps domdef: la grille est reguliere ou non ', systime(1)-tempdeux 390 ;------------------------------------------------------------------- 391 if keyword_set(memeindices) then begin 392 premierxu = premierxt & premierxv = premierxt & premierxf = premierxt 393 dernierxu = dernierxt & dernierxv = dernierxt & dernierxf = dernierxt 394 premieryu = premieryt & premieryv = premieryt & premieryf = premieryt 395 dernieryu = dernieryt & dernieryv = dernieryt & dernieryf = dernieryt 396 premierzw = premierzt 397 dernierzw = dernierzt 398 nxu = nxt & nxv = nxt & nxf = nxt 399 nyu = nyt & nyv = nyt & nyf = nyt 400 nzw = nzt 401 lon1 = glamt[premierxt, 0] 402 lon2 = glamf[dernierxf, 0] 403 lat1 = gphit[0, premieryt] 404 lat2 = gphif[0, dernieryf] 405 return 406 endif 407 ptu: 408 endif 409 ;------------------------------------------------------------------- 410 ; calcul pour les pts u 411 ;------------------------------------------------------------------- 412 if (where(grille eq 'U'))[0] NE -1 then begin 413 if (keyword_set(xindex) AND keyword_set(xindex)) OR keyword_set(index) then begin 414 premierxu = x1 415 premieryu = y1 416 dernierxu = x2 417 dernieryu = y2 418 ENDIF ELSE BEGIN 419 dom = where( (glamu ge lon1) $ 420 and (glamu le lon2) $ 421 and (gphiu ge lat1) $ 422 and (gphiu le lat2) ) 423 if (dom[0] eq -1) then begin 424 if keyword_set(findalways) then begin 425 domdef, prof1,prof2, GRILLE=grille, MEMEINDICES = memeindices, _extra = ex 426 return 427 ENDIF ELSE BEGIN 428 ras = report('ATTENTION, la boite ne contient pas de points u') 429 goto, ptv 430 ENDELSE 431 ENDIF 432 ordonnee = dom/jpi 433 abscisse = dom-ordonnee*jpi 434 ordonnee = ordonnee(uniq(ordonnee, sort(ordonnee))) 435 abscisse = abscisse(uniq(abscisse, sort(abscisse))) 436 premierxu = abscisse[0] 437 premieryu = ordonnee[0] 438 dernierxu = abscisse[n_elements(abscisse) -1] 439 dernieryu = ordonnee[n_elements(ordonnee) -1] 440 ENDELSE 441 nxu = dernierxu-premierxu+1 442 nyu = dernieryu-premieryu+1 443 if (total(glamu[premierxu:dernierxu,premieryu:dernieryu] $ 444 NE glamu[premierxu:dernierxu,0]#replicate(1, nyu)) eq 0 $ 445 AND total(gphiu[premierxu:dernierxu,premieryu:dernieryu] $ 446 NE replicate(1, nxu)#(gphiu[0,premieryu:dernieryu])[*]) eq 0 ) $ 447 then key_irregular = 0 ELSE key_irregular = 1 448 ptv: 449 endif 450 ;------------------------------------------------------------------- 451 ; calcul pour les pts v 452 ;------------------------------------------------------------------- 453 if (where(grille eq 'V'))[0] NE -1 then begin 454 if (keyword_set(xindex) AND keyword_set(xindex)) OR keyword_set(index) then begin 455 premierxv = x1 456 premieryv = y1 457 dernierxv = x2 458 dernieryv = y2 459 ENDIF ELSE BEGIN 460 dom = where( (glamv ge lon1) $ 461 and (glamv le lon2) $ 462 and (gphiv ge lat1) $ 463 and (gphiv le lat2) ) 464 if (dom[0] eq -1) then begin 465 if keyword_set(findalways) then begin 466 domdef, prof1,prof2, GRILLE=grille, MEMEINDICES = memeindices, _extra = ex 467 return 468 ENDIF ELSE BEGIN 469 ras = report('ATTENTION, la boite ne contient pas de points v') 470 goto, ptf 471 ENDELSE 472 ENDIF 473 ordonnee = dom/jpi 474 abscisse = dom-ordonnee*jpi 475 ordonnee = ordonnee(uniq(ordonnee, sort(ordonnee))) 476 abscisse = abscisse(uniq(abscisse, sort(abscisse))) 477 premierxv = abscisse[0] 478 premieryv = ordonnee[0] 479 dernierxv = abscisse[n_elements(abscisse) -1] 480 dernieryv = ordonnee[n_elements(ordonnee) -1] 481 ENDELSE 482 nxv = dernierxv-premierxv+1 483 nyv = dernieryv-premieryv+1 484 if (total(glamv[premierxv:dernierxv,premieryv:dernieryv] $ 485 NE glamv[premierxv:dernierxv,0]#replicate(1, nyv)) eq 0 $ 486 AND total(gphiv[premierxv:dernierxv,premieryv:dernieryv] $ 487 NE replicate(1, nxv)#(gphiv[0,premieryv:dernieryv])[*]) eq 0 ) $ 488 then key_irregular = 0 ELSE key_irregular = 1 489 ptf: 490 endif 491 ;------------------------------------------------------------------- 492 ; calcul pour les pts f 493 ;------------------------------------------------------------------- 494 if (where(grille eq 'F'))[0] NE -1 then begin 495 if (keyword_set(xindex) AND keyword_set(xindex)) OR keyword_set(index) then begin 496 premierxf = x1 497 premieryf = y1 498 dernierxf = x2 499 dernieryf = y2 500 ENDIF ELSE BEGIN 501 dom = where( (glamf ge lon1) $ 502 and (glamf le lon2) $ 503 and (gphif ge lat1) $ 504 and (gphif le lat2) ) 505 if (dom[0] eq -1) then begin 506 if keyword_set(findalways) then begin 507 domdef, prof1,prof2, GRILLE=grille, MEMEINDICES = memeindices, _extra = ex 508 return 509 ENDIF ELSE BEGIN 510 ras = report('ATTENTION, la boite ne contient pas de points f') 511 return 512 ENDELSE 513 ENDIF 514 ordonnee = dom/jpi 515 abscisse = dom-ordonnee*jpi 516 ordonnee = ordonnee(uniq(ordonnee, sort(ordonnee))) 517 abscisse = abscisse(uniq(abscisse, sort(abscisse))) 518 premierxf = abscisse[0] 519 premieryf = ordonnee[0] 520 dernierxf = abscisse[n_elements(abscisse) -1] 521 dernieryf = ordonnee[n_elements(ordonnee) -1] 522 ENDELSE 523 nxf = dernierxf-premierxf+1 524 nyf = dernieryf-premieryf+1 525 if (total(glamf[premierxf:dernierxf,premieryf:dernieryf] $ 526 NE glamf[premierxf:dernierxf,0]#replicate(1, nyf)) eq 0 $ 527 AND total(gphif[premierxf:dernierxf,premieryf:dernieryf] $ 528 NE replicate(1, nxf)#(gphif[0,premieryf:dernieryf])[*]) eq 0 ) $ 529 then key_irregular = 0 ELSE key_irregular = 1 530 endif 531 ;------------------------------------------------------------------- 532 if keyword_set(key_performance) THEN print, 'temps domdef', systime(1)-tempsun 752 ENDIF 753 ;------------------------------------------------------------------- 754 ; vertical domain defined with the Z index 755 ;------------------------------------------------------------------- 756 ENDIF ELSE BEGIN 757 CASE N_PARAMS() OF 758 2:fstz = min([x1, x2], max = lstz) 759 4:return 760 6:fstz = min([z1, z2], max = lstz) 761 ENDCASE 762 IF fstz LT 0 OR lstz GE jpk THEN BEGIN 763 ras = report('Bad definition of X1, X2, Z1 or Z2') 764 return 765 ENDIF 766 nz = lstz - fstz + 1 767 ; find vert1t, vert2t, firstzt, firstzt, nzt 768 ; according to (x1, x2) or (z1, z2) 769 IF (where(gridtype eq 'T'))[0] NE -1 THEN BEGIN 770 vert1t = min(gdept[fstz:lstz], max = vert2t) 771 firstzt = fstz & lastzt = lstz & nzt = nz 772 ENDIF 773 ; find vert1w, vert2w, firstzw, firstzw, nzw 774 ; according to (x1, x2) or (z1, z2) 775 IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN 776 vert1w = min(gdepw[fstz:lstz], max = vert2w) 777 firstzw = fstz & lastzw = lstz & nzw = nz 778 ENDIF 779 vert1 = min([vert1t, vert1w]) 780 vert2 = max([vert2t, vert2w]) 781 ENDELSE 782 ; 783 ;------------------------------------------------------------------- 784 IF NOT keyword_set(key_forgetold) THEN BEGIN 785 @updateold 786 ENDIF 787 ;------------------------------------------------------------------- 788 if keyword_set(key_performance) THEN print, 'temps domdef', systime(1)-tempsun 533 789 ;------------------------------------------------------------------- 534 790 535 791 ;------------------------------------------------------------------- 536 792 return 537 793 end -
trunk/ToBeReviewed/GRILLE/f2v.pro
r12 r13 38 38 ;------------------------------------------------------------ 39 39 FUNCTION f2v, temp 40 @common 40 ;--------------------------------------------------------- 41 @cm_4mesh 42 @cm_4data 43 @cm_4cal 44 IF NOT keyword_set(key_forgetold) THEN BEGIN 45 @updatenew 46 ENDIF 47 ;--------------------------------------------------------- 41 48 res = temp 42 49 ;on force nxt=nxf, etc ... 43 premierxv = premierxf44 dernierxv = dernierxf45 premieryv = premieryf46 dernieryv = dernieryf50 firstxv = firstxf 51 lastxv = lastxf 52 firstyv = firstyf 53 lastyv = lastyf 47 54 nxv = nxf 48 55 nyv = nyf 49 56 vargrid = 'V' 50 57 if NOT keyword_set(valmask) then valmask = 1e20 51 lon1 = glamv[ premierxv, 0]52 lon2 = glamf[ dernierxf, 0]58 lon1 = glamv[firstxv, 0] 59 lon2 = glamf[lastxf, 0] 53 60 54 61 ; cas sur la taille du tableau et application … … 60 67 taille[1] eq nxf and taille[2] eq nyf: 61 68 taille[1] eq jpi and taille[2] eq jpj: $ 62 res=res[ premierxf:dernierxf, premieryf:dernieryf]69 res=res[firstxf:lastxf, firstyf:lastyf] 63 70 else: $ 64 71 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 65 72 endcase 66 mask = (fmask())[ premierxf:dernierxf, premieryf:dernieryf, 0]73 mask = (fmask())[firstxf:lastxf, firstyf:lastyf, 0] 67 74 terre = where(mask EQ 0) 68 75 IF terre[0] NE -1 THEN res[terre] = !values.f_nan 69 76 res = 0.5*(res + shift(res, 1, 0)) 70 if NOT (keyword_set(key_periodi que) AND nxf EQ jpi) then res[0, *] = !values.f_nan71 mask = (vmask())[ premierxf:dernierxf, premieryf:dernieryf, 0]77 if NOT (keyword_set(key_periodic) AND nxf EQ jpi) then res[0, *] = !values.f_nan 78 mask = (vmask())[firstxf:lastxf, firstyf:lastyf, 0] 72 79 terre = where(mask EQ 0) 73 80 IF terre[0] NE -1 THEN res[terre] = valmask … … 77 84 taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ nzt: 78 85 taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ jpk: $ 79 res=res[*, *, premierzt:dernierzt]86 res=res[*, *, firstzt:lastzt] 80 87 taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ jpt: 81 88 taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk: $ 82 res=res[ premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt]89 res=res[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt] 83 90 taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpt: $ 84 res=res[ premierxf:dernierxf, premieryf:dernieryf, *]91 res=res[firstxf:lastxf, firstyf:lastyf, *] 85 92 else: $ 86 93 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 87 94 ENDCASE 88 95 if taille[3] EQ jpt then begin 89 mask = (fmask())[ premierxf:dernierxf, premieryf:dernieryf, dernierzt*(nzt NE jpk)]96 mask = (fmask())[firstxf:lastxf, firstyf:lastyf, lastzt*(nzt NE jpk)] 90 97 mask = temporary(mask[*])#replicate(1, jpt) 91 98 mask = reform(mask, nxf, nyf, jpt, /over) 92 ENDIF ELSE mask = (fmask())[ premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt]99 ENDIF ELSE mask = (fmask())[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt] 93 100 terre = where(temporary(mask) EQ 0) 94 101 IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan 95 102 res = 0.5*(res + shift(res, 1, 0, 0)) 96 if NOT (keyword_set(key_periodi que) AND nxf EQ jpi) then res[0, *, *] = !values.f_nan103 if NOT (keyword_set(key_periodic) AND nxf EQ jpi) then res[0, *, *] = !values.f_nan 97 104 if taille[3] EQ jpt then BEGIN 98 mask = tmask[ premierxf:dernierxf, premieryf:dernieryf, dernierzt*(nzt NE jpk)]105 mask = tmask[firstxf:lastxf, firstyf:lastyf, lastzt*(nzt NE jpk)] 99 106 mask = temporary(mask[*])#replicate(1, jpt) 100 107 mask = reform(mask, nxf, nyf, jpt, /over) 101 ENDIF ELSE mask = (vmask())[ premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt]108 ENDIF ELSE mask = (vmask())[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt] 102 109 terre = where(temporary(mask) EQ 0) 103 110 IF terre[0] NE -1 THEN res[temporary(terre)] = valmask … … 107 114 taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ nzt AND taille[4] EQ jpt: 108 115 taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 109 res=res[*, *, premierzt:dernierzt, *]116 res=res[*, *, firstzt:lastzt, *] 110 117 taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 111 res=res[ premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt, *]118 res=res[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt, *] 112 119 else: $ 113 120 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 114 121 ENDCASE 115 mask = (fmask())[ premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt]122 mask = (fmask())[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt] 116 123 mask = temporary(mask[*])#replicate(1, jpt) 117 124 mask = reform(mask, nxf, nyf, nzt, jpt, /over) … … 119 126 IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan 120 127 res = 0.5*(res + shift(res, 1, 0, 0, 0)) 121 if NOT (keyword_set(key_periodi que) AND nxf EQ jpi) then res[0, *, *, *] = !values.f_nan122 mask = (vmask())[ premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt]128 if NOT (keyword_set(key_periodic) AND nxf EQ jpi) then res[0, *, *, *] = !values.f_nan 129 mask = (vmask())[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt] 123 130 mask = temporary(mask[*])#replicate(1, jpt) 124 131 mask = reform(mask, nxf, nyf, nzt, jpt, /over) … … 128 135 endcase 129 136 137 IF NOT keyword_set(key_forgetold) THEN BEGIN 138 @updateold 139 ENDIF 140 130 141 return, res 131 142 END -
trunk/ToBeReviewed/GRILLE/fmask.pro
r12 r13 25 25 ;------------------------------------------------------------ 26 26 FUNCTION fmask 27 @common 28 tempsun = systime(1) ; pour key_performance 29 if jpk EQ 1 then begin 30 res = tmask*shift(tmask, -1, 0)*shift(tmask, 0, -1)*shift(tmask, -1, -1) 31 if NOT keyword_set(key_periodique) then res[jpi-1, *] = fmaskredy 32 res[*, jpj-1] = fmaskredx 33 ENDIF ELSE BEGIN 34 res = tmask*shift(tmask, -1, 0, 0)*shift(tmask, 0, -1, 0)*shift(tmask, -1, -1, 0) 35 if NOT keyword_set(key_periodique) then res[jpi-1, *, *] = fmaskredy 36 res[*, jpj-1, *] = fmaskredx 37 ENDELSE 38 if keyword_set(key_performance) THEN print, 'temps fmask', systime(1)-tempsun 39 40 return, res 27 ;--------------------------------------------------------- 28 @cm_4mesh 29 IF NOT keyword_set(key_forgetold) THEN BEGIN 30 @updatenew 31 ENDIF 32 ;--------------------------------------------------------- 33 tempsun = systime(1) ; pour key_performance 34 ; 35 CASE size(tmask, /n_dimensions) OF 36 2:res = tmask*shift(tmask, -1, 0)*shift(tmask, 0, -1)*shift(tmask, -1, -1) 37 3:res = tmask*shift(tmask, -1, 0, 0)*shift(tmask, 0, -1, 0)*shift(tmask, -1, -1, 0) 38 ENDCASE 39 ; 40 if NOT keyword_set(key_periodic) then res[jpi-1, *, *] = fmaskredy 41 res[*, jpj-1, *] = fmaskredx 42 ; 43 if keyword_set(key_performance) THEN print, 'temps fmask', systime(1)-tempsun 44 45 return, res 41 46 end -
trunk/ToBeReviewed/GRILLE/grille.pro
r12 r13 13 13 ; 14 14 ; CALLING SEQUENCE: 15 ; grille,mask,glam,gphi,gdep,nx,ny,nz, premierx,premiery,premierz,dernierx,derniery,dernierz,e1,e2,e315 ; grille,mask,glam,gphi,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz,e1,e2,e3 16 16 ; 17 17 ; INPUTS:rien. ATTENTION les choix de la grille se fait a partir de la … … 32 32 ; /NOTRI: utile seulement qd TRI est active. dans ce cas 33 33 ; grille retourne -1 ds la variable tri meme si la variable du 34 ; common triangles est definie et differente de -1 35 ; 36 ; OUTPUTS:mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz, 37 ; dernierx,derniery,dernierz,e1,e2,e3 34 ; common triangles_list est definie et differente de -1 35 ; 36 ; /WDEPTH: to specify that the field is at W depth instad of T 37 ; depth (automatically activated if vargrid eq 'W') 38 ; 39 ; OUTPUTS:mask,glam,gphi,gdep,nx,ny,nz,firstx,firsty,firstz, 40 ; lastx,lasty,lastz,e1,e2,e3 38 41 ; 39 42 ; pour leur definition cf domdef et la gestion des sous … … 43 46 ; mask, glam et gphi il suffit de taper grille, mask, glam, gphi 44 47 ; 45 ; COMMON BLOCKS: 46 ; common.pro congridseb.pro 48 ; COMMON BLOCKS: cm_4mesh and cm_4data 47 49 ; 48 50 ; SIDE EFFECTS: utilise la variable globale vargird … … 59 61 ;------------------------------------------------------------ 60 62 ;------------------------------------------------------------ 61 pro grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz, e1,e2,e3, TRI = tri, NOTRI = notri, TOUT = tout, FORPLT = forplt, _EXTRA = ex 62 @common 63 tempsun = systime(1) ; pour key_performance 64 ;------------------------------------------------------------ 65 if keyword_set(tout) then begin 66 oldboite = [lon1, lon2, lat1, lat2, prof1, prof2] 67 domdef, grille = vargrid, _EXTRA = ex 68 endif 69 tempdeux = systime(1) ; pour key_performance =2 70 CASE 1 OF 71 ;------------------------------------------------------------ 72 ; grille T 73 ;------------------------------------------------------------ 74 vargrid eq 'OPAPTDHT' or vargrid eq 'OPAPT3DT' $ 75 or vargrid eq 'T': begin 63 pro grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, TRI = tri, NOTRI = notri, TOUT = tout, FORPLT = forplt, IFPLTZ = ifpltz, WDEPTH = wdepth, _EXTRA = ex 64 ;------------------------------------------------------------ 65 ; include commons 66 @cm_4mesh 67 @cm_4data 68 IF NOT keyword_set(key_forgetold) THEN BEGIN 69 @updatenew 70 ENDIF 71 ;--------------------- 72 tempsun = systime(1) ; pour key_performance 73 ;------------------------------------------------------------ 74 vargrid = strupcase(strmid(vargrid,0,/reverse_offset)) 75 ; 76 if vargrid eq 'W' then wdepth = 1 77 if keyword_set(tout) then begin 78 savedbox = 1b 79 saveboxparam, 'boxparam4grille.dat' 80 domdef, gridtype = vargrid, _EXTRA = ex 81 endif 82 tempdeux = systime(1) ; pour key_performance =2 83 ;------------------------------------------------------------ 84 ;------------------------------------------------------------ 85 IF keyword_set(wdepth) THEN BEGIN 86 firstz = firstzw 87 lastz = lastzw 88 nz = nzw 89 ENDIF ELSE BEGIN 90 firstz = firstzt 91 lastz = lastzt 92 nz = nzt 93 ENDELSE 94 ;------------------------------------------------------------ 95 ;------------------------------------------------------------ 96 CASE 1 OF 97 ;------------------------------------------------------------ 98 ; grille T and W 99 ;------------------------------------------------------------ 100 vargrid eq 'T' OR vargrid eq 'W' : begin 76 101 ;scalaires 77 nx=nxt 78 ny=nyt 79 nz=nzt 80 premierx = premierxt 81 premiery = premieryt 82 premierz = premierzt 83 dernierx = dernierxt 84 derniery = dernieryt 85 dernierz = dernierzt 102 nx = nxt 103 ny = nyt 104 firstx = firstxt 105 firsty = firstyt 106 lastx = lastxt 107 lasty = lastyt 86 108 ;vecteurs 2d 87 glam=glamt[premierx:dernierx, premiery:derniery]88 gphi=gphit[premierx:dernierx, premiery:derniery]89 e1 =e1t[premierx:dernierx, premiery:derniery]90 e2 =e2t[premierx:dernierx, premiery:derniery]109 IF n_elements(glam) NE 1 THEN glam = glamt[firstx:lastx, firsty:lasty] 110 IF n_elements(gphi) NE 1 THEN gphi = gphit[firstx:lastx, firsty:lasty] 111 IF n_elements(e1) NE 1 THEN e1 = e1t[firstx:lastx, firsty:lasty] 112 IF n_elements(e2) NE 1 THEN e2 = e2t[firstx:lastx, firsty:lasty] 91 113 ;vecteurs 3d 92 mask=tmask[premierx:dernierx, premiery:derniery, premierz:dernierz] 93 end 94 ;------------------------------------------------------------ 95 ; grille W 96 ;------------------------------------------------------------ 97 vargrid eq 'OPAPT3DW' or vargrid eq 'W': begin 114 IF keyword_set(forplt) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz] $ 115 ELSE IF n_elements(mask) NE 1 THEN mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 116 end 117 ;------------------------------------------------------------ 118 ; grille U 119 ;------------------------------------------------------------ 120 vargrid eq 'U': begin 98 121 ;scalaires 99 nx=nxt 100 ny=nyt 101 nz=nzw 102 premierx = premierxt 103 premiery = premieryt 104 premierz = premierzw 105 dernierx = dernierxt 106 derniery = dernieryt 107 dernierz = dernierzw 122 nx = nxu 123 ny = nyu 124 firstx = firstxu 125 firsty = firstyu 126 lastx = lastxu 127 lasty = lastyu 108 128 ;vecteurs 2d 109 terre = where(tmask[*, *, 0] EQ 0) 110 glam=glamt[premierx:dernierx, premiery:derniery] 111 gphi=gphit[premierx:dernierx, premiery:derniery] 112 e1 =e1t[premierx:dernierx, premiery:derniery] 113 e2 =e2t[premierx:dernierx, premiery:derniery] 129 IF n_elements(glam) NE 1 THEN glam = glamu[firstx:lastx, firsty:lasty] 130 IF n_elements(gphi) NE 1 THEN gphi = gphiu[firstx:lastx, firsty:lasty] 131 if keyword_set(forplt) then BEGIN 132 mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz] 133 eastboarder = mask-shift(mask, 1, 0)*mask 134 westboarder = mask-shift(mask, -1, 0)*mask 135 if key_periodic NE 1 OR nx NE jpi then westboarder[nx-1, *] = 0b 136 tmp1 = shift(eastboarder, 0, 1) 137 tmp1[*, 0] = 0b 138 tmp2 = shift(eastboarder, 0, -1) 139 tmp2[*, ny-1] = 0b 140 add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder)) 141 eastboarder = temporary(eastboarder)+temporary(add) 142 tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b 143 tmp1[*, ny-1] = 1b 144 tmp1[*, 0] = 1b 145 tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b 146 if key_periodic NE 1 OR nx NE jpi then begin 147 tmp2[nx-1, *] = 1b 148 tmp2[0, *] = 0b 149 endif 150 no1 = temporary(tmp1)*temporary(tmp2) 151 tmp = temporary(eastboarder)*temporary(no1)*mask 152 mask[0:nx-2, *] = 0b 153 tmp = temporary(tmp)+temporary(mask) 154 tmp = where(tmp GE 1) 155 if tmp[0] NE -1 then begin 156 glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp] 157 gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp] 158 endif 159 ENDIF 160 IF n_elements(e1) NE 1 THEN e1 = e1u[firstx:lastx, firsty:lasty] 161 IF n_elements(e2) NE 1 THEN e2 = e2u[firstx:lastx, firsty:lasty] 114 162 ;vecteurs 3d 115 mask=tmask[premierx:dernierx, premiery:derniery, premierz:dernierz] 116 end 117 ;------------------------------------------------------------ 118 ; grille U 119 ;------------------------------------------------------------ 120 vargrid eq 'OPAPTDHU' or vargrid eq 'OPAPT3DU' $ 121 or vargrid eq 'U': begin 163 IF keyword_set(forplt) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz] $ 164 ELSE IF n_elements(mask) NE 1 THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] 165 end 166 ;------------------------------------------------------------ 167 ; grille V 168 ;------------------------------------------------------------ 169 vargrid eq 'OPAPTDHV' or vargrid eq 'OPAPT3DV' $ 170 or vargrid eq 'V': begin 122 171 ;scalaires 123 nx=nxu 124 ny=nyu 125 nz=nzt 126 premierx = premierxu 127 premiery = premieryu 128 premierz = premierzt 129 dernierx = dernierxu 130 derniery = dernieryu 131 dernierz = dernierzt 172 nx = nxv 173 ny = nyv 174 firstx = firstxv 175 firsty = firstyv 176 lastx = lastxv 177 lasty = lastyv 132 178 ;vecteurs 2d 133 glam=glamu[premierx:dernierx, premiery:derniery] 134 gphi=gphiu[premierx:dernierx, premiery:derniery] 135 if keyword_set(forplt) then BEGIN 136 mask = tmask[premierx:dernierx, premiery:derniery, niveau-1] 137 terre = where(mask EQ 0) 138 if terre[0] NE -1 then begin 139 glam[terre] = (glamt[premierx:dernierx, premiery:derniery])[terre] 140 gphi[terre] = (gphit[premierx:dernierx, premiery:derniery])[terre] 141 endif 142 ENDIF 143 e1 =e1u[premierx:dernierx, premiery:derniery] 144 e2 =e2u[premierx:dernierx, premiery:derniery] 179 IF n_elements(glam) NE 1 THEN glam = glamv[firstx:lastx, firsty:lasty] 180 IF n_elements(gphi) NE 1 THEN gphi = gphiv[firstx:lastx, firsty:lasty] 181 if keyword_set(forplt) then BEGIN 182 mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz] 183 northboarder = mask-shift(mask, 0, 1)*mask 184 southboarder = mask-shift(mask, 0, -1)*mask 185 southboarder[*, ny-1] = 0b 186 tmp1 = shift(northboarder, -1, 0) 187 if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b 188 tmp2 = shift(northboarder, 1, 0) 189 if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b 190 add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder) 191 northboarder = temporary(northboarder)+temporary(add) 192 tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b 193 tmp1[*, ny-1] = 1b 194 tmp1[*, 0] = 0b 195 tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b 196 if key_periodic NE 1 OR nx NE jpi then begin 197 tmp2[nx-1, *] = 1b 198 tmp2[0, *] = 1b 199 endif 200 no1 = temporary(tmp1)*temporary(tmp2) 201 tmp = temporary(northboarder)*mask*temporary(no1) 202 mask[*, 0:ny-2] = 0b 203 tmp = temporary(tmp)+temporary(mask) 204 tmp = where(tmp GE 1) 205 if tmp[0] NE -1 then begin 206 glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp] 207 gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp] 208 endif 209 ENDIF 210 IF n_elements(e1) NE 1 THEN e1 = e1v[firstx:lastx, firsty:lasty] 211 IF n_elements(e2) NE 1 THEN e2 = e2v[firstx:lastx, firsty:lasty] 145 212 ;vecteurs 3d 146 mask=(umask())[premierx:dernierx, premiery:derniery, premierz:dernierz] 147 end 148 ;------------------------------------------------------------ 149 ; grille V 150 ;------------------------------------------------------------ 151 vargrid eq 'OPAPTDHV' or vargrid eq 'OPAPT3DV' $ 152 or vargrid eq 'V': begin 213 IF keyword_set(forplt) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz] $ 214 ELSE IF n_elements(mask) NE 1 THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] 215 end 216 ;------------------------------------------------------------ 217 ; grille F 218 ;------------------------------------------------------------ 219 vargrid eq 'OPAPTDHF' or vargrid eq 'OPAPT3DF' $ 220 or vargrid eq 'F': begin 153 221 ;scalaires 154 nx=nxv 155 ny=nyv 156 nz=nzt 157 premierx = premierxv 158 premiery = premieryv 159 premierz = premierzt 160 dernierx = dernierxv 161 derniery = dernieryv 162 dernierz = dernierzt 222 nx = nxf 223 ny = nyf 224 firstx = firstxf 225 firsty = firstyf 226 lastx = lastxf 227 lasty = lastyf 163 228 ;vecteurs 2d 164 glam=glamv[premierx:dernierx, premiery:derniery] 165 gphi=gphiv[premierx:dernierx, premiery:derniery] 166 if keyword_set(forplt) then BEGIN 167 mask = tmask[premierx:dernierx, premiery:derniery, niveau-1] 168 terre = where(mask EQ 0) 169 if terre[0] NE -1 then begin 170 glam[terre] = (glamt[premierx:dernierx, premiery:derniery])[terre] 171 gphi[terre] = (gphit[premierx:dernierx, premiery:derniery])[terre] 172 endif 173 ENDIF 174 e1 =e1v[premierx:dernierx, premiery:derniery] 175 e2 =e2v[premierx:dernierx, premiery:derniery] 229 IF n_elements(glam) NE 1 THEN glam = glamf[firstx:lastx, firsty:lasty] 230 IF n_elements(gphi) NE 1 THEN gphi = gphif[firstx:lastx, firsty:lasty] 231 if keyword_set(forplt) then BEGIN 232 mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz] 233 eastboarder = mask-shift(mask, 1, 0)*mask 234 westboarder = mask-shift(mask, -1, 0)*mask 235 westboarder[nx-1, *] = 0b 236 northboarder = mask-shift(mask, 0, 1)*mask 237 southboarder = mask-shift(mask, 0, -1)*mask 238 southboarder[*, ny-1] = 0b 239 tmp1 = shift(northboarder, -1, 0) 240 if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b 241 tmp2 = shift(northboarder, 1, 0) 242 if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b 243 add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder) 244 northboarder = temporary(northboarder)+temporary(add) 245 tmp1 = shift(eastboarder, 0, 1) 246 tmp1[*, 0] = 0b 247 tmp2 = shift(eastboarder, 0, -1) 248 tmp2[*, ny-1] = 0b 249 add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder)) 250 eastboarder = temporary(eastboarder)+temporary(add) 251 tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b 252 tmp1[*, ny-1] = 1b 253 tmp1[*, 0] = 1b 254 tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b 255 if key_periodic NE 1 OR nx NE jpi then begin 256 tmp2[nx-1, *] = 1b 257 tmp2[0, *] = 1b 258 endif 259 no1 = temporary(tmp1)*temporary(tmp2) 260 tmp = (temporary(northboarder)+temporary(eastboarder))*mask*temporary(no1) 261 mask[0:nx-2, *] = 0b 262 mask[*, 0:ny-2] = 0b 263 tmp = temporary(tmp)+temporary(mask) 264 tmp = where(tmp GE 1) 265 if tmp[0] NE -1 then begin 266 glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp] 267 gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp] 268 endif 269 ENDIF 270 IF n_elements(e1) NE 1 THEN e1 = e1f[firstx:lastx, firsty:lasty] 271 IF n_elements(e2) NE 1 THEN e2 = e2f[firstx:lastx, firsty:lasty] 176 272 ;vecteurs 3d 177 mask=(vmask())[premierx:dernierx, premiery:derniery, premierz:dernierz] 178 end 179 ;------------------------------------------------------------ 180 ; grille F 181 ;------------------------------------------------------------ 182 vargrid eq 'OPAPTDHF' or vargrid eq 'OPAPT3DF' $ 183 or vargrid eq 'F': begin 184 ;scalaires 185 nx=nxf 186 ny=nyf 187 nz=nzt 188 premierx = premierxf 189 premiery = premieryf 190 premierz = premierzt 191 dernierx = dernierxf 192 derniery = dernieryf 193 dernierz = dernierzt 194 ;vecteurs 2d 195 glam=glamf[premierx:dernierx, premiery:derniery] 196 gphi=gphif[premierx:dernierx, premiery:derniery] 197 if keyword_set(forplt) then BEGIN 198 mask = tmask[premierx:dernierx, premiery:derniery, niveau-1] 199 terre = where(mask EQ 0) 200 if terre[0] NE -1 then begin 201 glam[terre] = (glamt[premierx:dernierx, premiery:derniery])[terre] 202 gphi[terre] = (gphit[premierx:dernierx, premiery:derniery])[terre] 203 endif 204 ENDIF 205 e1 =e1f[premierx:dernierx, premiery:derniery] 206 e2 =e2f[premierx:dernierx, premiery:derniery] 207 ;vecteurs 3d 208 mask=(fmask())[premierx:dernierx, premiery:derniery, premierz:dernierz] 273 IF keyword_set(forplt) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz] $ 274 ELSE IF n_elements(mask) NE 1 THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 275 END 276 ;------------------------------------------------------------ 277 ELSE:BEGIN 278 ras = report('Wrong definition of vargrid = '+vargrid+'. Only T, U, V, W or F are acceptable') 279 stop 280 END 281 ENDCASE 282 IF testvar(var = key_performance) EQ 2 THEN $ 283 print, 'temps grille: attribution des scalaires, vecteurs et tableaux ', systime(1)-tempdeux 284 ; 285 ;------------------------------------------------------------ 286 ;------------------------------------------------------------ 287 ;------------------------------------------------------------ 288 ; Variables se rapportant a la dimension verticale 289 ;------------------------------------------------------------ 290 ;------------------------------------------------------------ 291 ;------------------------------------------------------------ 292 ; 293 ; 294 tempdeux = systime(1) ; pour key_performance =2 295 if keyword_set(wdepth) then begin 296 gdep = gdepw[firstz:lastz] 297 e3 = e3w[firstz:lastz] 298 endif else begin 299 gdep = gdept[firstz:lastz] 300 e3 = e3t[firstz:lastz] 301 ENDELSE 302 ; for the vertical sections with partial steps 303 IF keyword_set(ifpltz) AND keyword_set(key_partialstep) THEN BEGIN 304 CASE 1 OF 305 ifpltz EQ 'xz' AND ny EQ 1:BEGIN 306 bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 307 good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth)) 308 bottom = lindgen(nx)+(bottom-1l+keyword_set(wdepth))*nx 309 IF good[0] NE -1 THEN BEGIN 310 bottom = bottom[good] 311 IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw) 312 gdep = replicate(1, nx)#gdep 313 if keyword_set(wdepth) THEN $ 314 truegdep = hdepw[firstx:lastx, firsty:lasty] $ 315 ELSE truegdep = hdept[firstx:lastx, firsty:lasty] 316 gdep[bottom] = truegdep[good] 317 ENDIF 209 318 END 210 ;------------------------------------------------------------ 211 ELSE:BEGIN 212 ras = report('Vauvaise definition de la variable vargrid: '+vargrid+'ceete variable doit etre T, U, V, W ou F') 213 stop 319 ifpltz EQ 'yz' AND nx EQ 1:BEGIN 320 bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 321 good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth)) 322 bottom = lindgen(ny)+(bottom-1l+keyword_set(wdepth))*ny 323 IF good[0] NE -1 THEN BEGIN 324 bottom = bottom[good] 325 IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw) 326 gdep = replicate(1, ny)#gdep 327 if keyword_set(wdepth) THEN $ 328 truegdep = hdepw[firstx:lastx, firsty:lasty] $ 329 ELSE truegdep = hdept[firstx:lastx, firsty:lasty] 330 gdep[bottom] = truegdep[good] 331 ENDIF 214 332 END 215 ENDCASE 216 IF testvar(var = key_performance) EQ 2 THEN $ 217 print, 'temps grille: attribution des scalaires, vecteurs et tableaux ', systime(1)-tempdeux 218 ;------------------------------------------------------------ 219 ; Variables se rapportant a la dimension verticale 220 ;------------------------------------------------------------ 221 tempdeux = systime(1) ; pour key_performance =2 222 if vargrid eq 'OPAPT3DW' or vargrid eq 'W' then begin 223 gdep = gdepw[premierz:dernierz] 224 e3=e3w[premierz:dernierz] 225 endif else begin 226 gdep = gdept[premierz:dernierz] 227 e3=e3t[premierz:dernierz] 228 endelse 229 IF testvar(var = key_performance) EQ 2 THEN $ 333 ELSE: 334 ENDCASE 335 ENDIF 336 IF testvar(var = key_performance) EQ 2 THEN $ 230 337 print, 'temps grille: Variables se rapportant a la dimension verticale ', systime(1)-tempdeux 231 338 ;------------------------------------------------------------ 232 339 ; vecteur triangulation Qd TRI est active 233 340 ;------------------------------------------------------------ 234 235 if triangles [0] EQ -1 OR keyword_set(notri) then tri = -1 ELSE BEGIN236 tempdeux = systime(1); pour key_performance =2237 238 msk[premierx:dernierx,premiery:derniery] = 1239 ind = where( msk[triangles[0, *]]*msk[triangles[1, *]]*msk[triangles[2, *]] EQ 1 )240 tri =triangles[*, ind]-(premierx+premiery*jpi)241 242 243 244 245 246 341 if arg_present(TRI) then $ 342 if triangles_list[0] EQ -1 OR keyword_set(notri) then tri = -1 ELSE BEGIN 343 tempdeux = systime(1) ; pour key_performance =2 344 msk = bytarr(jpi, jpj) 345 msk[firstx:lastx, firsty:lasty] = 1 346 ind = where( msk[triangles_list[0, *]]*msk[triangles_list[1, *]]*msk[triangles_list[2, *]] EQ 1 ) 347 tri = triangles_list[*, ind]-(firstx+firsty*jpi) 348 y = tri/jpi 349 x = tri-y*jpi 350 tri = x+y*nx 351 IF testvar(var = key_performance) EQ 2 THEN $ 352 print, 'temps grille: decoupage de la triangulation ', systime(1)-tempdeux 353 ENDELSE 247 354 ;------------------------------------------------------------------ 248 355 ; pour s'assurer qu'il n'y a pas de dimension degenerees (=1) … … 256 363 ; e3=reform(e3, /over) 257 364 258 if keyword_set(tout) then domdef, oldboite, grille = vargrid 259 if keyword_set(key_performance) THEN print, 'temps grille', systime(1)-tempsun 260 261 return 365 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grille.dat' 366 if keyword_set(key_performance) THEN print, 'temps grille', systime(1)-tempsun 367 368 ;------------------------------------------------------------ 369 IF NOT keyword_set(key_forgetold) THEN BEGIN 370 @updateold 371 ENDIF 372 ;--------------------- 373 return 262 374 263 375 end -
trunk/ToBeReviewed/GRILLE/t2v.pro
r12 r13 38 38 ;------------------------------------------------------------ 39 39 FUNCTION t2v, temp 40 @common 40 ;--------------------------------------------------------- 41 @cm_4mesh 42 @cm_4data 43 @cm_4cal 44 IF NOT keyword_set(key_forgetold) THEN BEGIN 45 @updatenew 46 ENDIF 47 ;--------------------------------------------------------- 41 48 res = temp 42 49 43 50 ;on force nxt=nxv, etc ... 44 premierxv = premierxt45 dernierxv = dernierxt46 premieryv = premieryt47 dernieryv = dernieryt51 firstxv = firstxt 52 lastxv = lastxt 53 firstyv = firstyt 54 lastyv = lastyt 48 55 nxv = nxt 49 56 nyv = nyt 50 57 vargrid = 'V' 51 58 if NOT keyword_set(valmask) then valmask = 1e20 52 lat1 = gphit[0, premieryt]53 lat2 = gphiv[0, dernieryv]59 lat1 = gphit[0, firstyt] 60 lat2 = gphiv[0, lastyv] 54 61 55 62 ; cas sur la taille du tableau et application … … 61 68 taille[1] eq nxt and taille[2] eq nyt: 62 69 taille[1] eq jpi and taille[2] eq jpj: $ 63 res=res[ premierxt:dernierxt, premieryt:dernieryt]70 res=res[firstxt:lastxt, firstyt:lastyt] 64 71 else: $ 65 72 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 66 73 endcase 67 mask = tmask[ premierxt:dernierxt, premieryt:dernieryt, 0]74 mask = tmask[firstxt:lastxt, firstyt:lastyt, 0] 68 75 terre = where(mask EQ 0) 69 76 IF terre[0] NE -1 THEN res[terre] = !values.f_nan 70 77 res = 0.5*(res + shift(res, 0, -1)) 71 78 res[*, nyt-1] = !values.f_nan 72 mask = (vmask())[ premierxt:dernierxt, premieryt:dernieryt, 0]79 mask = (vmask())[firstxt:lastxt, firstyt:lastyt, 0] 73 80 terre = where(mask EQ 0) 74 81 IF terre[0] NE -1 THEN res[terre] = valmask … … 78 85 taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ nzt: 79 86 taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ jpk: $ 80 res=res[*, *, premierzt:dernierzt]87 res=res[*, *, firstzt:lastzt] 81 88 taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ jpt: 82 89 taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk: $ 83 res=res[ premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt]90 res=res[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt] 84 91 taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpt: $ 85 res=res[ premierxt:dernierxt, premieryt:dernieryt, *]92 res=res[firstxt:lastxt, firstyt:lastyt, *] 86 93 else: $ 87 94 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 88 95 ENDCASE 89 96 if taille[3] EQ jpt then begin 90 mask = tmask[ premierxt:dernierxt, premieryt:dernieryt, dernierzt*(nzt NE jpk)]97 mask = tmask[firstxt:lastxt, firstyt:lastyt, lastzt*(nzt NE jpk)] 91 98 mask = temporary(mask[*])#replicate(1, jpt) 92 99 mask = reform(mask, nxt, nyt, jpt, /over) 93 ENDIF ELSE mask = tmask[ premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt]100 ENDIF ELSE mask = tmask[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt] 94 101 terre = where(temporary(mask) EQ 0) 95 102 IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan … … 97 104 res[*, nyt-1, *] = !values.f_nan 98 105 if taille[3] EQ jpt then BEGIN 99 mask = (vmask())[ premierxt:dernierxt, premieryt:dernieryt, dernierzt*(nzt NE jpk)]106 mask = (vmask())[firstxt:lastxt, firstyt:lastyt, lastzt*(nzt NE jpk)] 100 107 mask = temporary(mask[*])#replicate(1, jpt) 101 108 mask = reform(mask, nxt, nyt, jpt, /over) 102 ENDIF ELSE mask = (vmask())[ premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt]109 ENDIF ELSE mask = (vmask())[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt] 103 110 terre = where(temporary(mask) EQ 0) 104 111 IF terre[0] NE -1 THEN res[temporary(terre)] = valmask … … 108 115 taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ nzt AND taille[4] EQ jpt: 109 116 taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 110 res=res[*, *, premierzt:dernierzt, *]117 res=res[*, *, firstzt:lastzt, *] 111 118 taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 112 res=res[ premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt, *]119 res=res[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt, *] 113 120 else: $ 114 121 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 115 122 ENDCASE 116 mask = tmask[ premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt]123 mask = tmask[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt] 117 124 mask = temporary(mask[*])#replicate(1, jpt) 118 125 mask = reform(mask, nxt, nyt, nzt, jpt, /over) … … 121 128 res = 0.5*(res + shift(res, 0, -1, 0, 0)) 122 129 res[*, nyt-1, *, *] = !values.f_nan 123 mask = (vmask())[ premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt]130 mask = (vmask())[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt] 124 131 mask = temporary(mask[*])#replicate(1, jpt) 125 132 mask = reform(mask, nxt, nyt, nzt, jpt, /over) … … 127 134 IF terre[0] NE -1 THEN res[temporary(terre)] = valmask 128 135 END 129 endcase 136 ENDCASE 137 138 IF NOT keyword_set(key_forgetold) THEN BEGIN 139 @updateold 140 ENDIF 130 141 131 142 return, res -
trunk/ToBeReviewed/GRILLE/tracegrille.pro
r12 r13 11 11 ; CALLING SEQUENCE:tracegrille 12 12 ; 13 ; INPUTS:none 13 ; INPUTS:glam et gphi, les tableaux 1d ou 2d des position en 14 ; longitude/latitude des points de la grille. Si glam et gphi ne sont 15 ; pas specifies, trace la grille specifiee par vargrid, sur le domaine 16 ; definit par le dernier domdef. 14 17 ; 15 18 ; KEYWORD PARAMETERS: … … 21 24 ; qu''une ligne de j constant tout les ystride points 22 25 ; 23 ; OCEAN: pour ne tracer la grille que sur les points oceans26 ; /OCEAN: pour ne tracer la grille que sur les points oceans 24 27 ; 25 ; EARTH: pour ne tracer la grille que sur les points terre 28 ; /EARTH: pour ne tracer la grille que sur les points terre 29 ; 30 ; /RMOUT:select to remove all cell having one corner out of the 31 ; plot boundaries (!x.range, !y.range) 26 32 ; 27 33 ; + tous les mots clefs de la procedure PLOTS … … 31 37 ; COMMON BLOCKS:common.pro 32 38 ; 33 ; SIDE EFFECTS:trace la grille specifiee par vargrid, sur le domaine 34 ; definit par le dernier domdef. 39 ; SIDE EFFECTS: 35 40 ; 36 41 ; RESTRICTIONS: … … 38 43 ; EXAMPLE: 39 44 ; 40 ; IDL> plt,indgen(jpi,jpj),/nocontour,/nocouleur 45 ; IDL> plt,indgen(jpi,jpj),/nocontour,/nofill 46 ; IDL> vargrid='T' 41 47 ; IDL> tracegrille,/ocean,color=20 42 48 ; IDL> tracegrille,/earth,color=80 … … 49 55 ;------------------------------------------------------------ 50 56 ;------------------------------------------------------------ 51 PRO tracegrille, OCEAN = ocean, EARTH = earth, XSTRIDE = xstride, YSTRIDE = ystride, _extra = extra 52 @common 53 tempsun = systime(1) ; pour key_performance 54 if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 57 PRO tracegrille, glamin, gphiin, OCEAN = ocean, EARTH = earth $ 58 , XSTRIDE = xstride, YSTRIDE = ystride, RMOUT = rmout $ 59 , _extra = extra 60 ;--------------------------------------------------------- 61 @cm_4mesh 62 @cm_4data 63 IF NOT keyword_set(key_forgetold) THEN BEGIN 64 @updatenew 65 ENDIF 66 ;--------------------------------------------------------- 67 tempsun = systime(1) ; pour key_performance 68 ; to avoid warning message 69 oldexcept = !except 70 !except = 0 71 if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 72 ; 73 if n_elements(glamin) * n_elements(gphiin) EQ 0 then BEGIN 74 grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz 75 IF keyword_set(ocean) AND key_gridtype EQ 'c' THEN BEGIN 76 ; we reduce the mask to take into account the point located ON the coastline. 77 CASE vargrid OF 78 'U':BEGIN 79 mask = tmask[firstx:lastx, firsty:lasty] 80 IF NOT keyword_set(key_periodic) OR nx NE jpi $ 81 THEN tmpx = mask[nx-1, *] 82 mask = (mask+shift(mask, -1, 0)) < 1 83 IF NOT keyword_set(key_periodic) OR nx NE jpi $ 84 THEN mask[nx-1, *] = temporary(tmpx) 85 END 86 'V':BEGIN 87 mask = tmask[firstx:lastx, firsty:lasty] 88 tmpy = mask[*, ny-1] 89 mask = (mask+shift(mask, 0, -1)) < 1 90 mask[*, ny-1] = temporary(tmpy) 91 END 92 'F':BEGIN 93 mask = tmask[firstx:lastx, firsty:lasty] 94 IF NOT keyword_set(key_periodic) OR nx NE jpi $ 95 THEN tmpx = mask[nx-1, *] 96 tmpy = mask[*, ny-1] 97 mask = (mask+shift(mask, -1, 0)+shift(mask, 0, -1)+shift(mask, -1, -1)) < 1 98 mask[*, ny-1] = temporary(tmpy) 99 IF NOT keyword_set(key_periodic) OR nx NE jpi $ 100 THEN mask[nx-1, *] = temporary(tmpx) 101 END 102 ELSE: 103 ENDCASE 104 ENDIF 105 ENDIF ELSE BEGIN 106 glam = glamin 107 gphi = gphiin 108 IF (size(glam))[0] EQ 1 AND (size(gphi))[0] EQ 1 THEN BEGIN 109 nx = n_elements(glam) 110 ny = n_elements(gphi) 111 glam = glam#replicate(1, ny) 112 gphi = replicate(1, nx)#gphi 113 ENDIF ELSE BEGIN 114 nx = (size(glam))[1] 115 ny = (size(glam))[2] 116 ENDELSE 117 ENDELSE 118 if n_elements(mask) EQ 0 then mask = replicate(1b, nx, ny) 119 if (size(mask))[0] EQ 3 then mask = mask[*, *, 0] 120 ; 121 IF keyword_set(RMOUT) THEN BEGIN 122 out = where(glam GT max(!x.range) OR glam LT min(!x.range) $ 123 OR gphi GT max(!y.range) OR gphi LT min(!y.range)) 124 IF out[0] NE -1 THEN BEGIN 125 glam[out] = !values.f_nan 126 gphi[out] = !values.f_nan 127 ENDIF 128 ENDIF 129 ; 130 IF keyword_set(ocean) then BEGIN 131 earth = where(mask EQ 0) 132 if earth[0] NE -1 then begin 133 glam[earth] = !values.f_nan 134 gphi[earth] = !values.f_nan 135 ENDIF 136 earth = 0 137 ENDIF 138 ; 139 IF keyword_set(earth) THEN BEGIN 140 ocean = where(mask EQ 1) 141 if ocean[0] NE -1 then begin 142 glam[ocean] = !values.f_nan 143 gphi[ocean] = !values.f_nan 144 ENDIF 145 ocean = 0 146 ENDIF 147 ; 148 if NOT keyword_set(xstride) then xstride = 1 149 if NOT keyword_set(ystride) then ystride = 1 150 case key_gridtype of 151 'c':BEGIN 152 for i = 0, ny-1, ystride do begin 153 plots, glam[*, i], gphi[*, i], _extra = extra 154 endfor 155 for i = 0, nx-1, xstride do begin 156 plots, glam[i, *], gphi[i, *], _extra = extra 157 endfor 158 END 159 'e':BEGIN 160 shifted = glam[0, 0] LT glam[0, 1] 161 glam2 = glam+(glam[1]-glam[0])/2. 162 if shifted then begin 163 for i = 0, ny-2 do BEGIN 164 xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*] 165 yy = (transpose([[gphi[*, i]], [gphi[*, i+1]]]))[*] 166 plots, xx[0:2*nx-2], yy[0:2*nx-2], _extra = extra 167 ENDFOR 168 ENDIF ELSE BEGIN 169 for i = 1, ny-1 do BEGIN 170 xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*] 171 yy = (transpose([[gphi[*, i]], [gphi[*, i-1]]]))[*] 172 plots, xx[0:2*nx-2], yy[0:2*nx-2], _extra = extra 173 ENDFOR 174 ENDELSE 175 for i = 1, (ny-1)/2 do $ 176 plots, [glam[0, 2*i-1], glam[0, 2*i]] $ 177 , [gphi[0, 2*i-1], gphi[0, 2*i]], _extra = extra 178 for i = 0, (ny-2)/2 do $ 179 plots, [glam[nx-1, 2*i], glam[nx-1, 2*i+1]] $ 180 , [gphi[nx-1, 2*i], gphi[nx-1, 2*i+1]], _extra = extra 181 END 182 endcase 55 183 56 grille, mask, glam, gphi, gdep, nx, ny 57 if (size(mask))[0] EQ 3 then mask = mask[*, *, 0] 58 case 1 of 59 keyword_set(ocean):BEGIN 60 earth = where(mask EQ 0) 61 if earth[0] NE -1 then begin 62 glam[earth] = !values.f_nan 63 gphi[earth] = !values.f_nan 64 endif 65 END 66 keyword_set(earth):BEGIN 67 ocean = where(mask EQ 1) 68 if ocean[0] NE -1 then begin 69 glam[ocean] = !values.f_nan 70 gphi[ocean] = !values.f_nan 71 endif 72 END 73 ELSE: 74 endcase 75 if NOT keyword_set(xstride) then xstride = 1 76 if NOT keyword_set(ystride) then ystride = 1 77 case key_gridtype of 78 'c':BEGIN 79 for i = 0, ny-1, ystride do begin 80 plots, glam[*, i], gphi[*, i], color = c_cote, _extra = extra 81 endfor 82 for i = 0, nx-1, xstride do begin 83 plots, glam[i, *], gphi[i, *], color = c_cote, _extra = extra 84 endfor 85 END 86 'e':BEGIN 87 shifted = glam[0, 0] LT glam[0, 1] 88 glam2 = glam+(glam[1]-glam[0])/2. 89 if shifted then begin 90 for i = 0, ny-2 do BEGIN 91 xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*] 92 yy = (transpose([[gphi[*, i]], [gphi[*, i+1]]]))[*] 93 plots, xx[0:2*nx-2], yy[0:2*nx-2], color = c_cote, _extra = extra 94 ENDFOR 95 ENDIF ELSE BEGIN 96 for i = 1, ny-1 do BEGIN 97 xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*] 98 yy = (transpose([[gphi[*, i]], [gphi[*, i-1]]]))[*] 99 plots, xx[0:2*nx-2], yy[0:2*nx-2], color = c_cote, _extra = extra 100 ENDFOR 101 ENDELSE 102 for i = 1, (ny-1)/2 do $ 103 plots, [glam[0, 2*i-1], glam[0, 2*i]] $ 104 , [gphi[0, 2*i-1], gphi[0, 2*i]], color = c_cote, _extra = extra 105 for i = 0, (ny-2)/2 do $ 106 plots, [glam[nx-1, 2*i], glam[nx-1, 2*i+1]] $ 107 , [gphi[nx-1, 2*i], gphi[nx-1, 2*i+1]], color = c_cote, _extra = extra 108 END 109 endcase 184 if keyword_set(key_performance) THEN print, 'temps trace grille', systime(1)-tempsun 185 !except = oldexcept 110 186 111 if keyword_set(key_performance) THEN print, 'temps trace grille', systime(1)-tempsun 112 return 187 return 113 188 end -
trunk/ToBeReviewed/GRILLE/u2t.pro
r12 r13 38 38 ;------------------------------------------------------------ 39 39 FUNCTION u2t, temp 40 @common 40 ;--------------------------------------------------------- 41 @cm_4mesh 42 @cm_4data 43 @cm_4cal 44 IF NOT keyword_set(key_forgetold) THEN BEGIN 45 @updatenew 46 ENDIF 47 ;--------------------------------------------------------- 41 48 res = temp 42 49 ;on force nxt=nxu, etc ... 43 premierxt = premierxu44 dernierxt = dernierxu45 premieryt = premieryu46 dernieryt = dernieryu50 firstxt = firstxu 51 lastxt = lastxu 52 firstyt = firstyu 53 lastyt = lastyu 47 54 nxt = nxu 48 55 nyt = nyu 49 56 vargrid = 'T' 50 57 if NOT keyword_set(valmask) then valmask = 1e20 51 lon1 = glamt[ premierxt, 0]52 lon2 = glamu[ dernierxu, 0]58 lon1 = glamt[firstxt, 0] 59 lon2 = glamu[lastxu, 0] 53 60 ; 54 61 ; cas sur la taille du tableau et application … … 60 67 taille[1] eq nxu and taille[2] eq nyu: 61 68 taille[1] eq jpi and taille[2] eq jpj: $ 62 res=res[ premierxu:dernierxu, premieryu:dernieryu]69 res=res[firstxu:lastxu, firstyu:lastyu] 63 70 else: $ 64 71 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 65 72 endcase 66 mask = (umask())[ premierxu:dernierxu, premieryu:dernieryu, 0]73 mask = (umask())[firstxu:lastxu, firstyu:lastyu, 0] 67 74 terre = where(mask EQ 0) 68 75 IF terre[0] NE -1 THEN res[terre] = !values.f_nan 69 76 res = 0.5*(res + shift(res, 1, 0)) 70 if NOT (keyword_set(key_periodi que) AND nxu EQ jpi) then res[0, *] = !values.f_nan71 mask = tmask[ premierxu:dernierxu, premieryu:dernieryu, 0]77 if NOT (keyword_set(key_periodic) AND nxu EQ jpi) then res[0, *] = !values.f_nan 78 mask = tmask[firstxu:lastxu, firstyu:lastyu, 0] 72 79 terre = where(mask EQ 0) 73 80 IF terre[0] NE -1 THEN res[terre] = valmask … … 77 84 taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ nzt: 78 85 taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ jpk: $ 79 res=res[*, *, premierzt:dernierzt]86 res=res[*, *, firstzt:lastzt] 80 87 taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ jpt: 81 88 taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk: $ 82 res=res[ premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt]89 res=res[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt] 83 90 taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpt: $ 84 res=res[ premierxu:dernierxu, premieryu:dernieryu, *]91 res=res[firstxu:lastxu, firstyu:lastyu, *] 85 92 else: $ 86 93 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 87 94 ENDCASE 88 95 if taille[3] EQ jpt then begin 89 mask = (umask())[ premierxu:dernierxu, premieryu:dernieryu, dernierzt*(nzt NE jpk)]96 mask = (umask())[firstxu:lastxu, firstyu:lastyu, lastzt*(nzt NE jpk)] 90 97 mask = temporary(mask[*])#replicate(1, jpt) 91 98 mask = reform(mask, nxu, nyu, jpt, /over) 92 ENDIF ELSE mask = (umask())[ premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt]99 ENDIF ELSE mask = (umask())[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt] 93 100 terre = where(temporary(mask) EQ 0) 94 101 IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan 95 102 res = 0.5*(res + shift(res, 1, 0, 0)) 96 if NOT (keyword_set(key_periodi que) AND nxu EQ jpi) then res[0, *, *] = !values.f_nan103 if NOT (keyword_set(key_periodic) AND nxu EQ jpi) then res[0, *, *] = !values.f_nan 97 104 if taille[3] EQ jpt then BEGIN 98 mask = tmask[ premierxu:dernierxu, premieryu:dernieryu, dernierzt*(nzt NE jpk)]105 mask = tmask[firstxu:lastxu, firstyu:lastyu, lastzt*(nzt NE jpk)] 99 106 mask = temporary(mask[*])#replicate(1, jpt) 100 107 mask = reform(mask, nxu, nyu, jpt, /over) 101 ENDIF ELSE mask = tmask[ premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt]108 ENDIF ELSE mask = tmask[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt] 102 109 terre = where(temporary(mask) EQ 0) 103 110 IF terre[0] NE -1 THEN res[temporary(terre)] = valmask … … 107 114 taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ nzt AND taille[4] EQ jpt: 108 115 taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 109 res=res[*, *, premierzt:dernierzt, *]116 res=res[*, *, firstzt:lastzt, *] 110 117 taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 111 res=res[ premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt, *]118 res=res[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt, *] 112 119 else: $ 113 120 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 114 121 ENDCASE 115 mask = (umask())[ premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt]122 mask = (umask())[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt] 116 123 mask = temporary(mask[*])#replicate(1, jpt) 117 124 mask = reform(mask, nxu, nyu, nzt, jpt, /over) … … 119 126 IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan 120 127 res = 0.5*(res + shift(res, 1, 0, 0, 0)) 121 if NOT (keyword_set(key_periodi que) AND nxu EQ jpi) then res[0, *, *, *] = !values.f_nan122 mask = tmask[ premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt]128 if NOT (keyword_set(key_periodic) AND nxu EQ jpi) then res[0, *, *, *] = !values.f_nan 129 mask = tmask[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt] 123 130 mask = temporary(mask[*])#replicate(1, jpt) 124 131 mask = reform(mask, nxu, nyu, nzt, jpt, /over) … … 128 135 endcase 129 136 137 IF NOT keyword_set(key_forgetold) THEN BEGIN 138 @updateold 139 ENDIF 140 130 141 return, res 131 142 END -
trunk/ToBeReviewed/GRILLE/umask.pro
r12 r13 39 39 ;------------------------------------------------------------ 40 40 FUNCTION umask 41 tempsun = systime(1) ; pour key_performance 42 @common 43 if jpk EQ 1 then begin 44 res = tmask*shift(tmask, -1, 0) 45 if NOT keyword_set(key_periodique) then res[jpi-1, *] = umaskred 46 ENDIF ELSE BEGIN 47 res = tmask*shift(tmask, -1, 0, 0) 48 if NOT keyword_set(key_periodique) then res[jpi-1, *, *] = umaskred 49 ENDELSE 41 ;--------------------------------------------------------- 42 @cm_4mesh 43 IF NOT keyword_set(key_forgetold) THEN BEGIN 44 @updatenew 45 ENDIF 46 ;--------------------------------------------------------- 47 tempsun = systime(1) ; pour key_performance 50 48 ; 51 if keyword_set(key_performance) THEN print, 'temps umask', systime(1)-tempsun 49 CASE size(tmask, /n_dimensions) OF 50 2:res = tmask*shift(tmask, -1, 0) 51 3:res = tmask*shift(tmask, -1, 0, 0) 52 ENDCASE 53 ; 54 if NOT keyword_set(key_periodic) then res[jpi-1, *, *] = umaskred 55 if keyword_set(key_performance) THEN print, 'temps umask', systime(1)-tempsun 52 56 ; 53 57 return, res 54 58 end -
trunk/ToBeReviewed/GRILLE/v2t.pro
r12 r13 38 38 ;------------------------------------------------------------ 39 39 FUNCTION v2t, temp 40 @common 40 ;--------------------------------------------------------- 41 @cm_4mesh 42 @cm_4data 43 @cm_4cal 44 IF NOT keyword_set(key_forgetold) THEN BEGIN 45 @updatenew 46 ENDIF 47 ;--------------------------------------------------------- 41 48 res = temp 42 49 ;on force nxt=nxv, etc ... 43 premierxt = premierxv44 dernierxt = dernierxv45 premieryt = premieryv46 dernieryt = dernieryv50 firstxt = firstxv 51 lastxt = lastxv 52 firstyt = firstyv 53 lastyt = lastyv 47 54 nxt = nxv 48 55 nyt = nyv 49 56 vargrid = 'T' 50 57 if NOT keyword_set(valmask) then valmask = 1e20 51 lat1 = gphit[0, premieryt]52 lat2 = gphiv[0, dernieryv]58 lat1 = gphit[0, firstyt] 59 lat2 = gphiv[0, lastyv] 53 60 54 61 ; cas sur la taille du tableau et application … … 60 67 taille[1] eq nxv and taille[2] eq nyv: 61 68 taille[1] eq jpi and taille[2] eq jpj: $ 62 res=res[ premierxv:dernierxv, premieryv:dernieryv]69 res=res[firstxv:lastxv, firstyv:lastyv] 63 70 else: $ 64 71 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 65 72 endcase 66 mask = (vmask())[ premierxv:dernierxv, premieryv:dernieryv, 0]73 mask = (vmask())[firstxv:lastxv, firstyv:lastyv, 0] 67 74 terre = where(mask EQ 0) 68 75 IF terre[0] NE -1 THEN res[terre] = !values.f_nan 69 76 res = 0.5*(res + shift(res, 0, +1)) 70 77 res[*, 0] = !values.f_nan 71 mask = tmask[ premierxv:dernierxv, premieryv:dernieryv, 0]78 mask = tmask[firstxv:lastxv, firstyv:lastyv, 0] 72 79 terre = where(mask EQ 0) 73 80 IF terre[0] NE -1 THEN res[terre] = valmask … … 77 84 taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ nzt: 78 85 taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ jpk: $ 79 res=res[*, *, premierzt:dernierzt]86 res=res[*, *, firstzt:lastzt] 80 87 taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ jpt: 81 88 taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk: $ 82 res=res[ premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt]89 res=res[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt] 83 90 taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpt: $ 84 res=res[ premierxv:dernierxv, premieryv:dernieryv, *]91 res=res[firstxv:lastxv, firstyv:lastyv, *] 85 92 else: $ 86 93 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 87 94 ENDCASE 88 95 if taille[3] EQ jpt then begin 89 mask = (vmask())[ premierxv:dernierxv, premieryv:dernieryv, dernierzt*(nzt NE jpk)]96 mask = (vmask())[firstxv:lastxv, firstyv:lastyv, lastzt*(nzt NE jpk)] 90 97 mask = temporary(mask[*])#replicate(1, jpt) 91 98 mask = reform(mask, nxv, nyv, jpt, /over) 92 ENDIF ELSE mask = (vmask())[ premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt]99 ENDIF ELSE mask = (vmask())[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt] 93 100 terre = where(temporary(mask) EQ 0) 94 101 IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan … … 96 103 res[*, 0, *] = !values.f_nan 97 104 if taille[3] EQ jpt then BEGIN 98 mask = tmask[ premierxv:dernierxv, premieryv:dernieryv, dernierzt*(nzt NE jpk)]105 mask = tmask[firstxv:lastxv, firstyv:lastyv, lastzt*(nzt NE jpk)] 99 106 mask = temporary(mask[*])#replicate(1, jpt) 100 107 mask = reform(mask, nxv, nyv, jpt, /over) 101 ENDIF ELSE mask = tmask[ premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt]108 ENDIF ELSE mask = tmask[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt] 102 109 terre = where(temporary(mask) EQ 0) 103 110 IF terre[0] NE -1 THEN res[temporary(terre)] = valmask … … 107 114 taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ nzt AND taille[4] EQ jpt: 108 115 taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 109 res=res[*, *, premierzt:dernierzt, *]116 res=res[*, *, firstzt:lastzt, *] 110 117 taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 111 res=res[ premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt, *]118 res=res[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt, *] 112 119 else: $ 113 120 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 114 121 ENDCASE 115 mask = (vmask())[ premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt]122 mask = (vmask())[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt] 116 123 mask = temporary(mask[*])#replicate(1, jpt) 117 124 mask = reform(mask, nxv, nyv, nzt, jpt, /over) … … 120 127 res = 0.5*(res + shift(res, 0, +1, 0, 0)) 121 128 res[*, 0, *, *] = !values.f_nan 122 mask = tmask[ premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt]129 mask = tmask[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt] 123 130 mask = temporary(mask[*])#replicate(1, jpt) 124 131 mask = reform(mask, nxv, nyv, nzt, jpt, /over) … … 128 135 endcase 129 136 130 return, res 137 IF NOT keyword_set(key_forgetold) THEN BEGIN 138 @updateold 139 ENDIF 140 141 return, res 131 142 END 132 143 -
trunk/ToBeReviewed/GRILLE/vmask.pro
r12 r13 27 27 FUNCTION vmask 28 28 @common 29 tempsun = systime(1) ; pour key_performance 30 if jpk EQ 1 then begin 31 res = tmask*shift(tmask, 0, -1) 32 res[*, jpj-1] = vmaskred 33 ENDIF ELSE BEGIN 34 res = tmask*shift(tmask, 0, -1, 0) 35 res[*, jpj-1, *] = vmaskred 36 ENDELSE 37 if keyword_set(key_performance) THEN print, 'temps vmask', systime(1)-tempsun 29 tempsun = systime(1) ; pour key_performance 38 30 ; 39 return, res 31 CASE size(tmask, /n_dimensions) OF 32 2:res = tmask*shift(tmask, 0, -1) 33 3:res = tmask*shift(tmask, 0, -1, 0) 34 ENDCASE 35 ; 36 res[*, jpj-1, *] = vmaskred 37 if keyword_set(key_performance) THEN print, 'temps vmask', systime(1)-tempsun 38 ; 39 return, res 40 40 end -
trunk/ToBeReviewed/GRILLE/whichgrid.pro
r12 r13 4 4 case filetype of 5 5 'OPA C-Grid.Binary IEEE: Meshmask':BEGIN 6 mesh lec, filename, /pasblabla, _extra = ex6 meshread, filename, /pasblabla, _extra = ex 7 7 END 8 8 'OPA C-Grid.Net Cdf: Meshmask':BEGIN 9 ncdf_mesh lec, filename, _extra = ex9 ncdf_meshread, filename, _extra = ex 10 10 END 11 11 'Regular 2D.Binary IEEE: mask': … … 28 28 ;********************************************************************* 29 29 FUNCTION getgridparameter, top 30 @common 31 30 ;--------------------------------------------------------- 31 @cm_4mesh 32 IF NOT keyword_set(key_forgetold) THEN BEGIN 33 @updatenew 34 ENDIF 35 ;--------------------------------------------------------- 32 36 widget_control, widget_info(top, find_by_uname = 'xmesh'), get_value = answer 33 37 jpiglo = long(answer[0]) … … 48 52 , get_value = answer 49 53 key_shift = long(answer[0]) 50 widget_control, widget_info(top, find_by_uname = 'key_periodi que') $54 widget_control, widget_info(top, find_by_uname = 'key_periodic') $ 51 55 , get_value = answer 52 key_periodi que= long(answer[0])56 key_periodic = long(answer[0]) 53 57 widget_control, widget_info(top, find_by_uname = 'triangulation') $ 54 58 , get_value = answer … … 61 65 , iyminmesh:iyminmesh, iymaxmesh:iymaxmesh $ 62 66 , izminmesh:izminmesh, izmaxmesh:izmaxmesh $ 63 , key_shift:key_shift, key_periodi que:key_periodique$67 , key_shift:key_shift, key_periodic:key_periodic $ 64 68 , triangulation:triangulation, boundary:boundary} 69 70 @updateold 65 71 66 72 return, res … … 68 74 ;********************************************************************* 69 75 pro showgridparameter, basetop, EDITABLE = editable, _EXTRA = ex 70 ; 71 @common 72 ;------------------------------------------------------------ 73 ; 76 ;--------------------------------------------------------- 77 @cm_4mesh 78 IF NOT keyword_set(key_forgetold) THEN BEGIN 79 @updatenew 80 ENDIF 81 ;--------------------------------------------------------- 74 82 widget_control, basetop, update = 0 75 83 base=widget_base(basetop, /COLUMN, /align_center, _EXTRA = ex) … … 82 90 if NOT keyword_set(key_shift) then key_shift = 0 83 91 nothing = widget_text(basea, value = strtrim(key_shift*(1-keyword_set(editable)),1), uname = 'key_shift', xsize = 4, EDITABLE = editable) 84 nothing = widget_label(basea, value = 'key_periodi que')85 if NOT keyword_set(key_periodi que) then key_periodique= 086 nothing = widget_text(basea, value = strtrim(key_periodi que*(1-keyword_set(editable)),1), uname = 'key_periodique', xsize = 4, EDITABLE = editable)92 nothing = widget_label(basea, value = 'key_periodic') 93 if NOT keyword_set(key_periodic) then key_periodic = 0 94 nothing = widget_text(basea, value = strtrim(key_periodic*(1-keyword_set(editable)),1), uname = 'key_periodic', xsize = 4, EDITABLE = editable) 87 95 baseb=widget_base(base, /row, /align_center) 88 96 nothing = widget_label(baseb, value = 'use a triangulation (y/n) ?') … … 159 167 ;********************************************************************* 160 168 PRO whichgrid_event, event 161 @common 169 ;--------------------------------------------------------- 170 @cm_4mesh 171 IF NOT keyword_set(key_forgetold) THEN BEGIN 172 @updatenew 173 ENDIF 174 ;--------------------------------------------------------- 162 175 widget_control, event.id, get_uvalue = eventuvalue 163 176 IF chkstru(eventuvalue,'name') EQ 0 THEN return … … 181 194 if rstrpos(filename, '.pro') EQ strlen(filename)-4 then begin 182 195 createpro, '@'+strmid(filename, 0, strlen(filename)-4) $ 183 , filename = isadirectory(io = homedir, title = 'Bad definition of Homedir') $ 184 +'for_createpro.pro' 196 , filename = myuniquetmpdir +'for_createpro.pro' 185 197 showgridparameter, event.handler, group_leader = event.handler,/frame, uname = 'showgridparameter' 186 198 ENDIF ELSE BEGIN … … 204 216 if rstrpos(filename, '.pro') EQ strlen(filename)-4 then begin 205 217 createpro, '@'+strmid(filename, 0, strlen(filename)-4) $ 206 , filename = isadirectory(io = homedir, title = 'Bad definition of Homedir') $ 207 +'for_createpro.pro' 218 , filename = myuniquetmpdir +'for_createpro.pro' 208 219 showgridparameter, event.handler, group_leader = event.handler,/frame, uname = 'showgridparameter' 209 220 ENDIF ELSE BEGIN … … 223 234 case filetype of 224 235 'OPA C-Grid.Binary IEEE: Meshmask':BEGIN 225 mesh lec, filename, /pasblabla, /getdimensions236 meshread, filename, /pasblabla, /getdimensions 226 237 showgridparameter, event.handler, /editable, group_leader = event.handler,/frame, uname = 'showgridparameter' 227 238 END 228 239 'OPA C-Grid.Net Cdf: Meshmask':BEGIN 229 ncdf_mesh lec, filename, /getdimensions240 ncdf_meshread, filename, /getdimensions 230 241 showgridparameter, event.handler, /editable, group_leader = event.handler,/frame, uname = 'showgridparameter' 231 242 END … … 269 280 END 270 281 ;********************************************************************* 271 FUNCTION whichgrid, name, IODIRECTORY = iodirectory, PARENT = parent, _EXTRA = ex 272 @common 273 ; 282 FUNCTION whichgrid, name, PARENT = parent, _EXTRA = ex 283 ;--------------------------------------------------------- 284 @cm_4mesh 285 IF NOT keyword_set(key_forgetold) THEN BEGIN 286 @updatenew 287 ENDIF 288 ;--------------------------------------------------------- 274 289 ; 275 290 if n_elements(name) NE 0 then begin 276 filename = isafile(filename = name )291 filename = isafile(filename = name, IODIRECTORY = iodir, _extra = ex) 277 292 if size(filename, /type) NE 7 then return, -1 278 293 ENDIF ELSE filaname = 'no file' … … 301 316 if rstrpos(filename, '.pro') EQ strlen(filename)-4 then begin 302 317 createpro, '@'+strmid(filename, 0, strlen(filename)-4) $ 303 , filename = isadirectory(io = homedir, title = 'Bad definition of Homedir') $ 304 +'for_createpro.pro' 318 , filename = myuniquetmpdir +'for_createpro.pro' 305 319 showgridparameter, base, group_leader = base,/frame, uname = 'showgridparameter' 306 320 ENDIF ELSE BEGIN
Note: See TracChangeset
for help on using the changeset viewer.