Changeset 223 for trunk/SRC/ToBeReviewed/CALCULS/grossemoyenne.pro
- Timestamp:
- 03/14/07 18:13:39 (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SRC/ToBeReviewed/CALCULS/grossemoyenne.pro
r163 r223 4 4 ;+ 5 5 ; 6 ; @file_comments 6 ; @file_comments 7 7 ; averages a 3- or 4-d time series field over a selected 8 8 ; geographical area or along the time axis. For one ore more … … 16 16 ; @param DIREC {in}{required} 17 17 ; 'x' 'y' 'z' 't' 'xy' 'xz' 'yz' 'xyz' 'xt' 'yt' 'zt' 18 ; 'xyt' 'xzt' 'yzt' or 'xyzt' 19 ; 20 ; @keyword BOXZOOM 18 ; 'xyt' 'xzt' 'yzt' or 'xyzt' 19 ; 20 ; @keyword BOXZOOM 21 21 ; [xmin,xmax,ymin,ymax (,(zmin,)zmax)] to more details, see domdef 22 ; boxzoom can have 5 forms: 23 ; [vert2], [vert1, vert2],[lon1, lon2, lat1, lat2], 22 ; boxzoom can have 5 forms: 23 ; [vert2], [vert1, vert2],[lon1, lon2, lat1, lat2], 24 24 ; [lon1, lon2, lat1, lat2, vert2],[lon1, lon2, lat1, lat2, vert1,vert2] 25 ; 26 ; @keyword NAN 27 ; not a number, we activate it if we want to average without considerate some masked values of TAB. 28 ; If masked values of TAB are values consecrated by IDL(!values.f_nan), wr just have to put NAN. 29 ; If masked values of TAB are valued a (a must be different of 1, corresponding to nan = 30 ; !values.f_nan and of 0, which desactivate nan). We have to put NAN=a. 31 ; Comment: In output, rsult points which are NAN will be valued a or !values.f_nan. 32 ; 25 ; 26 ; @keyword NAN 27 ; not a number, we activate it if we want to average without considerate some 28 ; masked values of TAB. 29 ; If masked values of TAB are values consecrated by IDL(!values.f_nan), we 30 ; just have to put NAN. 31 ; If masked values of TAB are valued a (a must be different of 1, 32 ; corresponding to nan = !values.f_nan and of 0, which desactivate nan). 33 ; We have to put NAN=a. 34 ; Comment: In output, rsult points which are NAN will be valued a or 35 ; !values.f_nan. 36 ; 33 37 ; @keyword NODOMDEF 34 ; We activate it if we do not want to pass in domdef even if the keyword boxzoom35 ; is present (like when grossemoyenne is called via checkfield)36 ; 37 ; @keyword INTEGRATION 38 ; We activate it if we do not want to pass in domdef even if the keyword 39 ; boxzoom is present (like when grossemoyenne is called via checkfield) 40 ; 41 ; @keyword INTEGRATION 38 42 ; To make an integral rather than an average 39 43 ; 40 ; @keyword SPATIALFIRST 44 ; @keyword SPATIALFIRST 41 45 ; when performing at the same time 42 46 ; spatial and temporal mean, grossemoyenne is assuming … … 48 52 ; SPATIALFIRST is activated automatically. 49 53 ; 50 ; @keyword TEMPORALFIRST 54 ; @keyword TEMPORALFIRST 51 55 ; to force to perform first temporal 52 ; mean even if nan is activated (see SPATIALFIRST 53 ; explanations...) 54 ; 55 ; 56 ; @keyword WDEPTH 57 ; to specify that the field is at W depth instead of T 58 ; depth (automatically activated if vargrid eq 'W') 56 ; mean even if nan is activated (see SPATIALFIRST explanations...) 57 ; 58 ; @keyword WDEPTH 59 ; to specify that the field is at W depth instead of T 60 ; depth (automatically activated if vargrid eq 'W') 59 61 ; 60 62 ; @uses 61 ; result:un tableau 63 ; result:un tableau 62 64 ; common 63 65 ; domdef 64 66 ; 65 ; @restrictions Put values corresponding to land at 1.e20 66 ; 67 ; @history 67 ; @restrictions 68 ; Put values corresponding to land at 1.e20 69 ; 70 ; @history 68 71 ; Jerome Vialard (jv\@lodyc.jussieu.fr) 69 72 ; 2/7/98 70 73 ; Sebastien Masson (smasson\@lodyc.jussieu.fr) 71 ; adaptation array containing a temporal dimension 74 ; adaptation array containing a temporal dimension 72 75 ; 14/8/98 73 76 ; 15/1/98 … … 118 121 taille = size(tab) 119 122 case 1 of 120 taille[0] eq 1 :return, report(' Le tableau n''a qu''une dimension, cas non traite!')121 taille[0] eq 2 :return, report(' Le tableau n''a qu''deux dimension, cas non traite!')122 taille[0] eq 3 :BEGIN 123 taille[0] eq 1 :return, report('array has only one dimension, not implemented!') 124 taille[0] eq 2 :return, report('array has only two dimensions, not implemented!') 125 taille[0] eq 3 :BEGIN 123 126 dim = '3d' 124 127 if dirx eq 0 and diry eq 0 and dirt eq 0 then return, tab 125 128 END 126 taille[0] eq 4 :BEGIN 129 taille[0] eq 4 :BEGIN 127 130 dim = '4d' 128 131 if dirx eq 0 and diry eq 0 and dirz eq 0 and dirt eq 0 then return, tab 129 132 END 130 else : return, report(' Le tableau d entree doit etre a 3 ou 4 dimensions s''il ne contient pas de dim temporelle utilisermoyenne')133 else : return, report('array must have 3 or 4 dimensions if there is not time dimension use moyenne') 131 134 endcase 132 135 ;------------------------------------------------------------ … … 134 137 ; Redefinition of the domain ajusted at boxzoom (at 6 elements) 135 138 ; This will allowed us to calculate only in the domain concerned by the average. 136 ; Domdef, followed by grid give us all arrays of the grid on the subdomain 137 ;------------------------------------------------------------ 138 if keyword_set(boxzoom) then BEGIN 139 ; Domdef, followed by grid give us all arrays of the grid on the subdomain 140 ;------------------------------------------------------------ 141 if keyword_set(boxzoom) then BEGIN 139 142 Case 1 Of 140 143 N_Elements(Boxzoom) Eq 1: bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]] … … 145 148 Else: return, report('Wrong Definition of Boxzoom') 146 149 endcase 147 if NOT keyword_set(nodomdef) then BEGIN 150 if NOT keyword_set(nodomdef) then BEGIN 148 151 savedbox = 1b 149 152 saveboxparam, 'boxparam4grmoyenne.dat' 150 153 domdef, bte, GRIDTYPE = vargrid, _extra = ex 151 ENDIF 152 ENDIF 154 ENDIF 155 ENDIF 153 156 ;--------------------------------------------------------------- 154 157 ; attribution of the mask and of longitude and latitude arrays... … … 159 162 ;------------------------------------------------------------ 160 163 if dirt EQ 1 AND NOT keyword_set(spatialfirst) then begin 161 if dim EQ 3d then BEGIN 164 if dim EQ 3d then BEGIN 162 165 case 1 of 163 166 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $ … … 165 168 , firsty:firsty+ny-1, *] 166 169 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpt:res = tab 167 else:BEGIN 170 else:BEGIN 168 171 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 169 172 return, report('Probleme d''adequation entre les tailles du domaine nx*ny*jpt '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)) … … 179 182 res = total(res, 3, nan = keyword_set(nan))/ (1 > divi) 180 183 if notanum[0] NE -1 then res[temporary(notanum)] = !values.f_nan 181 ENDIF ELSE res = total(res, 3)/(1.*taille[3]) 184 ENDIF ELSE res = total(res, 3)/(1.*taille[3]) 182 185 ENDELSE 183 186 ENDIF ELSE BEGIN … … 190 193 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpk and taille[4] eq jpt: $ 191 194 res = tab[*, *, firstz:lastz, *] 192 else:BEGIN 195 else:BEGIN 193 196 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 194 197 return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $ … … 211 214 ENDELSE 212 215 ENDELSE 213 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 216 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 214 217 return, moyenne(temporary(res), direc, BOXZOOM = boxzoom, NAN = nan, INTEGRATION = integration, NODOMDEF = nodomdef, WDEPTH = wdepth, _extra = ex) 215 218 ENDIF ELSE res = tab 216 219 IF jpt EQ 1 THEN BEGIN 217 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 220 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 218 221 return, moyenne(reform(res, /over), direc, BOXZOOM = boxzoom, NAN = nan, INTEGRATION = integration, NODOMDEF = nodomdef, WDEPTH = wdepth, _extra = ex) 219 222 END … … 227 230 ; II.1) verification of the coherence of the array's size to average 228 231 ; Verification of the coherence between the array's size and the domain defined by domdef 229 ; The input array must have either the total domain's size (jpi,jpj,jpt) or this230 ; one of the reduced domain (nx,ny,jpt)232 ; The input array must have either the total domain's size (jpi,jpj,jpt) or 233 ; this one of the reduced domain (nx,ny,jpt) 231 234 ;--------------------------------------------------------------- 232 235 case 1 of … … 235 238 , firsty:firsty+ny-1, *] 236 239 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpt:res = tab 237 else:BEGIN 240 else:BEGIN 238 241 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 239 242 return, report('Probleme d''adequation entre les tailles du domaine nx*ny*jpt '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)) 240 243 enD 241 244 endcase 242 if keyword_set(nan) NE 0 then BEGIN 243 if nan NE 1 then BEGIN ; If nan is not !values.f_nan 245 if keyword_set(nan) NE 0 then BEGIN 246 if nan NE 1 then BEGIN ; If nan is not !values.f_nan 244 247 ; we put it at !values.f_nan 245 248 if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ … … 249 252 ENDIF 250 253 ;--------------------------------------------------------------- 251 ; Comment : WE HAVE TO BE CAREFUL ABOUT CASES WHERE THE DIMENSION TO AVERAGE = 1,252 ; A ND MAKE SURE THAT IT EXIST. THAT IS WHY WE USE reform(...,nx,ny,...) WHICH CAN253 ; LOOK USELESS AT THE BEGINNING254 ;--------------------------------------------------------------- 255 if nx EQ 1 OR ny EQ 1 then BEGIN 254 ; Comment : WE HAVE TO BE CAREFUL ABOUT CASES WHERE THE DIMENSION TO 255 ; AVERAGE = 1, AND MAKE SURE THAT IT EXIST. THAT IS WHY WE USE 256 ; reform(...,nx,ny,...) WHICH CAN LOOK USELESS AT THE BEGINNING 257 ;--------------------------------------------------------------- 258 if nx EQ 1 OR ny EQ 1 then BEGIN 256 259 res = reform(res, nx, ny, jpt, /over) 257 260 e1 = reform(e1, nx, ny, /over) … … 270 273 echelle = (temporary(e))[*]#replicate(1, jpt) 271 274 echelle = reform(echelle, nx, ny, jpt, /over) 272 if keyword_set(integration) then divi = 1 ELSE BEGIN 275 if keyword_set(integration) then divi = 1 ELSE BEGIN 273 276 IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 1) $ 274 277 ELSE divi = total(echelle, 1) … … 287 290 echelle = (temporary(e))[*]#replicate(1, jpt) 288 291 echelle = reform(echelle, nx, ny, jpt, /over) 289 if keyword_set(integration) then divi = 1 ELSE BEGIN 292 if keyword_set(integration) then divi = 1 ELSE BEGIN 290 293 IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 2) $ 291 294 ELSE divi = total(echelle, 2) 292 ENDELSE 295 ENDELSE 293 296 res = total(temporary(res)*echelle, 2, nan = nan)/(divi > 1.) 294 297 if msknan[0] NE -1 then begin … … 302 305 echelle = (temporary(e1)*temporary(e2)*temporary(mask))[*]#replicate(1, jpt) 303 306 echelle = reform(echelle, nx, ny, jpt, /over) 304 if keyword_set(integration) then divi = 1 ELSE BEGIN 307 if keyword_set(integration) then divi = 1 ELSE BEGIN 305 308 IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 1) $ 306 309 ELSE divi = total(total(echelle, 1), 1) … … 322 325 if (dim eq '4d') then begin 323 326 ;--------------------------------------------------------------- 324 ; III.1) Verification of the coherence of the array to average size 325 ; Verification of the coherence between the array's size and the domain 326 ; defind by domdef 327 ; The input array must have either the total domain size (jpi,jpj,jpk,jpt) 327 ; III.1) Verification of the coherence of the array to average size 328 ; Verification of the coherence between the array's size and the domain 329 ; defind by domdef 330 ; The input array must have either the total domain size (jpi,jpj,jpk,jpt) 328 331 ; or this one of the reduced domain (nx,ny,ny,jpt) 329 332 ;--------------------------------------------------------------- … … 336 339 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpk and taille[4] eq jpt: $ 337 340 res = tab[*, *, firstz:lastz, *] 338 else:BEGIN 341 else:BEGIN 339 342 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 340 343 return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $ … … 346 349 endcase 347 350 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 OR jpt EQ 1 then res = reform(res, nx, ny, nz, jpt, /over) 348 if keyword_set(nan) NE 0 then BEGIN 349 if nan NE 1 then BEGIN ; if nan is not !values.f_nan 351 if keyword_set(nan) NE 0 then BEGIN 352 if nan NE 1 then BEGIN ; if nan is not !values.f_nan 350 353 ; we put it at !values.f_nan 351 354 if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ … … 355 358 ENDIF 356 359 ;--------------------------------------------------------------- 357 ; Comment : WE HAVE TO BE CAREFUL ABOUT CASES WHERE THE DIMENSION TO AVERAGE = 1,358 ; A ND MAKE SURE THAT IT EXIST. THAT IS WHY WE USE reform(...,nx,ny,...) WHICH CAN359 ; LOOK USELESS AT THE BEGINNING360 ;--------------------------------------------------------------- 361 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then BEGIN 360 ; Comment : WE HAVE TO BE CAREFUL ABOUT CASES WHERE THE DIMENSION TO 361 ; AVERAGE = 1, AND MAKE SURE THAT IT EXIST. THAT IS WHY WE USE 362 ; reform(...,nx,ny,...) WHICH CAN LOOK USELESS AT THE BEGINNING 363 ;--------------------------------------------------------------- 364 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then BEGIN 362 365 res = reform(res, nx, ny, nz, jpt, /over) 363 366 mask = reform(mask, nx, ny, nz, /over) … … 366 369 ; the top of the ocean floor is 367 370 IF vargrid EQ 'T' OR vargrid EQ 'W' THEN bottom = total(mask, 3) $ 368 ELSE bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 369 ; we suppress columns with only ocean or land 371 ELSE bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 372 ; we suppress columns with only ocean or land 370 373 good = where(bottom NE 0 AND bottom NE nz) 371 374 ; the bottom of the ocean in 3D index is: … … 386 389 IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 387 390 AND nx NE 1 THEN BEGIN 388 IF msknan[0] EQ -1 THEN BEGIN 391 IF msknan[0] EQ -1 THEN BEGIN 389 392 msknan = replicate(1b, nx, ny, nz, jpt) 390 393 nan = 1 … … 395 398 res[temporary(bottom)] = !values.f_nan 396 399 ENDIF 397 if keyword_set(integration) then divi = 1 ELSE begin 400 if keyword_set(integration) then divi = 1 ELSE begin 398 401 IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 1) $ 399 402 ELSE divi = total(echelle, 1) … … 415 418 IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 416 419 AND ny NE 1 THEN BEGIN 417 IF msknan[0] EQ -1 THEN BEGIN 420 IF msknan[0] EQ -1 THEN BEGIN 418 421 msknan = replicate(1b, nx, ny, nz) 419 422 nan = 1 … … 424 427 res[temporary(bottom)] = !values.f_nan 425 428 ENDIF 426 if keyword_set(integration) then divi = 1 ELSE begin 429 if keyword_set(integration) then divi = 1 ELSE begin 427 430 IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 2) $ 428 431 ELSE divi = total(echelle, 2) … … 446 449 echelle = (temporary(e33)*temporary(mask))[*]#replicate(1, jpt) 447 450 echelle = reform(echelle, nx, ny, nz, jpt, /over) 448 if keyword_set(integration) then divi = 1 ELSE begin 451 if keyword_set(integration) then divi = 1 ELSE begin 449 452 IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 3) $ 450 453 ELSE divi = total(echelle, 3) … … 467 470 IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 468 471 AND nx*ny NE 1 THEN BEGIN 469 IF msknan[0] EQ -1 THEN BEGIN 472 IF msknan[0] EQ -1 THEN BEGIN 470 473 msknan = replicate(1b, nx, ny, nz) 471 474 nan = 1 … … 476 479 res[temporary(bottom)] = !values.f_nan 477 480 ENDIF 478 if keyword_set(integration) then divi = 1 ELSE begin 481 if keyword_set(integration) then divi = 1 ELSE begin 479 482 IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 1) $ 480 483 ELSE divi = total(total(echelle, 1), 1) … … 497 500 echelle = (temporary(e133[*])*temporary(mask[*]))#replicate(1, jpt) 498 501 echelle = reform(echelle, nx, ny, nz, jpt, /over) 499 if keyword_set(integration) then divi = 1 ELSE begin 502 if keyword_set(integration) then divi = 1 ELSE begin 500 503 IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 2) $ 501 504 ELSE divi = total(total(echelle, 1), 2) … … 518 521 echelle = (temporary(e233[*])*temporary(mask[*]))#replicate(1, jpt) 519 522 echelle = reform(echelle, nx, ny, nz, jpt, /over) 520 if keyword_set(integration) then divi = 1 ELSE begin 523 if keyword_set(integration) then divi = 1 ELSE begin 521 524 IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 2), 2) $ 522 525 ELSE divi = total(total(echelle, 2), 2) … … 539 542 echelle = (temporary(e1233[*])*temporary(mask[*]))#replicate(1, jpt) 540 543 echelle = reform(echelle, nx, ny, nz, jpt, /over) 541 if keyword_set(integration) then divi = 1 ELSE begin 544 if keyword_set(integration) then divi = 1 ELSE begin 542 545 IF msknan[0] NE -1 THEN divi = total(total(total(echelle*msknan, 1), 1), 1) $ 543 546 ELSE divi = total(total(total(echelle, 1), 1), 1) … … 562 565 tdim = size(res, /n_dimensions) 563 566 if keyword_set(integration) then res = total(res, tdim, nan = nan) $ 564 ELSE BEGIN 565 if keyword_set(nan) then BEGIN 567 ELSE BEGIN 568 if keyword_set(nan) then BEGIN 566 569 testnan = testnan < 1 567 570 testnan = total(temporary(testnan), tdim) … … 569 572 ENDIF ELSE divi = jpt 570 573 res = total(res, tdim, nan = nan)/(1 > divi) 571 ENDELSE 574 ENDELSE 572 575 ENDIF 573 576 ;------------------------------------------------------------ … … 580 583 valmask = 1e+20 581 584 terre = where(divi EQ 0) 582 IF terre[0] NE -1 THEN BEGIN 585 IF terre[0] NE -1 THEN BEGIN 583 586 res[temporary(terre)] = 1e+20 584 ENDIF 587 ENDIF 585 588 ;------------------------------------------------------------ 586 589 ; IV.2) We replace, when nan equal 1, !values.f_nan by nan 587 590 ;------------------------------------------------------------ 588 if keyword_set(nan) NE 0 then BEGIN 591 if keyword_set(nan) NE 0 then BEGIN 589 592 puttonan = where(temporary(testnan) EQ 0) 590 593 if puttonan[0] NE -1 then res[temporary(puttonan)] = !values.f_nan 591 if nan NE 1 then BEGIN 594 if nan NE 1 then BEGIN 592 595 notanumber = where(finite(res) eq 0) 593 596 if notanumber[0] NE -1 then res[temporary(notanumber)] = nan … … 599 602 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 600 603 ;------------------------------------------------------------ 601 if keyword_set(key_performance) THEN print, 'temps grossemoyenne', systime(1)-tempsun 604 if keyword_set(key_performance) THEN print, 'temps grossemoyenne', systime(1)-tempsun 602 605 return, res 603 606 ;------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.