Ignore:
Timestamp:
03/14/07 18:13:39 (17 years ago)
Author:
pinsard
Message:

improvements of some *.pro

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/SRC/ToBeReviewed/CALCULS/grossemoyenne.pro

    r163 r223  
    44;+ 
    55; 
    6 ; @file_comments   
     6; @file_comments 
    77; averages a 3- or 4-d time series field over a selected 
    88; geographical area or along the time axis. For one ore more 
     
    1616; @param DIREC {in}{required} 
    1717; '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 
    2121; [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], 
    2424; [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; 
    3337; @keyword NODOMDEF 
    34 ; We activate it if we do not want to pass in domdef even if the keyword boxzoom  
    35 ; 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 
    3842; To make an integral rather than an average 
    3943; 
    40 ; @keyword SPATIALFIRST  
     44; @keyword SPATIALFIRST 
    4145; when performing at the same time 
    4246; spatial and temporal mean, grossemoyenne is assuming 
     
    4852; SPATIALFIRST is activated automatically. 
    4953; 
    50 ; @keyword TEMPORALFIRST  
     54; @keyword TEMPORALFIRST 
    5155; 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') 
    5961; 
    6062; @uses 
    61 ; result:un tableau  
     63; result:un tableau 
    6264; common 
    6365; domdef 
    6466; 
    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 
    6871; Jerome Vialard (jv\@lodyc.jussieu.fr) 
    6972;                       2/7/98 
    7073;                       Sebastien Masson (smasson\@lodyc.jussieu.fr) 
    71 ; adaptation array containing a temporal dimension  
     74; adaptation array containing a temporal dimension 
    7275;                       14/8/98 
    7376;                       15/1/98 
     
    118121  taille = size(tab) 
    119122  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 
    123126      dim = '3d' 
    124127      if dirx eq 0 and diry eq 0 and dirt eq 0 then return, tab 
    125128    END 
    126     taille[0] eq 4 :BEGIN  
     129    taille[0] eq 4 :BEGIN 
    127130      dim = '4d' 
    128131      if dirx eq 0 and diry eq 0 and dirz eq 0 and dirt eq 0 then return, tab 
    129132    END 
    130     else : return, report('Le tableau d entree doit etre a 3 ou 4 dimensions s''il ne contient pas de dim temporelle utiliser moyenne') 
     133    else : return, report('array must have 3 or 4 dimensions if there is not time dimension use moyenne') 
    131134  endcase 
    132135;------------------------------------------------------------ 
     
    134137; Redefinition of the domain ajusted at boxzoom (at 6 elements) 
    135138; 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 
    139142    Case 1 Of 
    140143      N_Elements(Boxzoom) Eq 1: bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]] 
     
    145148      Else: return, report('Wrong Definition of Boxzoom') 
    146149    endcase 
    147     if NOT keyword_set(nodomdef) then BEGIN  
     150    if NOT keyword_set(nodomdef) then BEGIN 
    148151      savedbox = 1b 
    149152      saveboxparam, 'boxparam4grmoyenne.dat' 
    150153      domdef, bte, GRIDTYPE = vargrid, _extra = ex 
    151     ENDIF  
    152   ENDIF  
     154    ENDIF 
     155  ENDIF 
    153156;--------------------------------------------------------------- 
    154157; attribution of the mask and of longitude and latitude arrays... 
     
    159162;------------------------------------------------------------ 
    160163  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 
    162165      case 1 of 
    163166        taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $ 
     
    165168                    , firsty:firsty+ny-1, *] 
    166169        taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpt:res = tab 
    167         else:BEGIN  
     170        else:BEGIN 
    168171          if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
    169172          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)) 
     
    179182          res = total(res, 3, nan = keyword_set(nan))/ (1 > divi) 
    180183          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]) 
    182185      ENDELSE 
    183186    ENDIF ELSE BEGIN 
     
    190193        taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk and taille[4] eq jpt: $ 
    191194          res = tab[*, *, firstz:lastz, *] 
    192         else:BEGIN  
     195        else:BEGIN 
    193196           if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
    194197          return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $ 
     
    211214      ENDELSE 
    212215    ENDELSE 
    213     if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'     
     216    if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
    214217    return,  moyenne(temporary(res), direc, BOXZOOM = boxzoom, NAN = nan, INTEGRATION = integration, NODOMDEF = nodomdef, WDEPTH = wdepth, _extra = ex) 
    215218  ENDIF ELSE res = tab 
    216219  IF jpt EQ 1 THEN BEGIN 
    217     if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'     
     220    if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
    218221    return, moyenne(reform(res, /over), direc, BOXZOOM = boxzoom, NAN = nan, INTEGRATION = integration, NODOMDEF = nodomdef, WDEPTH = wdepth, _extra = ex) 
    219222  END 
     
    227230;   II.1) verification of the coherence of the array's size to average 
    228231; 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 this  
    230 ; 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) 
    231234;--------------------------------------------------------------- 
    232235    case 1 of 
     
    235238                  , firsty:firsty+ny-1, *] 
    236239      taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpt:res = tab 
    237       else:BEGIN  
     240      else:BEGIN 
    238241        if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
    239242        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)) 
    240243      enD 
    241244    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 
    244247; we put it at !values.f_nan 
    245248        if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 
     
    249252    ENDIF 
    250253;--------------------------------------------------------------- 
    251 ; Comment : WE HAVE TO BE CAREFUL ABOUT CASES WHERE THE DIMENSION TO AVERAGE = 1,  
    252 ; AND MAKE SURE THAT IT EXIST. THAT IS WHY WE USE reform(...,nx,ny,...) WHICH CAN  
    253 ; LOOK USELESS AT THE BEGINNING 
    254 ;--------------------------------------------------------------- 
    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 
    256259      res = reform(res, nx, ny, jpt, /over) 
    257260      e1 =  reform(e1, nx, ny, /over) 
     
    270273        echelle = (temporary(e))[*]#replicate(1, jpt) 
    271274        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 
    273276          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 1) $ 
    274277          ELSE divi = total(echelle, 1) 
     
    287290        echelle = (temporary(e))[*]#replicate(1, jpt) 
    288291        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 
    290293          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 2) $ 
    291294          ELSE divi = total(echelle, 2) 
    292         ENDELSE  
     295        ENDELSE 
    293296        res = total(temporary(res)*echelle, 2, nan = nan)/(divi > 1.) 
    294297        if msknan[0] NE -1 then begin 
     
    302305        echelle = (temporary(e1)*temporary(e2)*temporary(mask))[*]#replicate(1, jpt) 
    303306        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 
    305308          IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 1) $ 
    306309          ELSE divi = total(total(echelle, 1), 1) 
     
    322325  if (dim eq '4d') then begin 
    323326;--------------------------------------------------------------- 
    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) 
    328331; or this one of the reduced domain (nx,ny,ny,jpt) 
    329332;--------------------------------------------------------------- 
     
    336339      taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk and taille[4] eq jpt: $ 
    337340        res = tab[*, *, firstz:lastz, *] 
    338       else:BEGIN  
     341      else:BEGIN 
    339342        if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
    340343        return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $ 
     
    346349    endcase 
    347350    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 
    350353; we put it at !values.f_nan 
    351354        if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 
     
    355358    ENDIF 
    356359;--------------------------------------------------------------- 
    357 ; Comment : WE HAVE TO BE CAREFUL ABOUT CASES WHERE THE DIMENSION TO AVERAGE = 1,  
    358 ; AND MAKE SURE THAT IT EXIST. THAT IS WHY WE USE reform(...,nx,ny,...) WHICH CAN  
    359 ; LOOK USELESS AT THE BEGINNING 
    360 ;--------------------------------------------------------------- 
    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 
    362365      res = reform(res, nx, ny, nz, jpt, /over) 
    363366      mask =  reform(mask, nx, ny, nz, /over) 
     
    366369; the top of the ocean floor is 
    367370      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 
    370373      good = where(bottom NE 0 AND bottom NE nz) 
    371374; the bottom of the ocean in 3D index is: 
     
    386389        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 
    387390          AND nx NE 1 THEN BEGIN 
    388           IF msknan[0] EQ -1 THEN BEGIN  
     391          IF msknan[0] EQ -1 THEN BEGIN 
    389392            msknan = replicate(1b, nx, ny, nz, jpt) 
    390393            nan = 1 
     
    395398          res[temporary(bottom)] = !values.f_nan 
    396399        ENDIF 
    397         if keyword_set(integration) then divi = 1 ELSE begin  
     400        if keyword_set(integration) then divi = 1 ELSE begin 
    398401          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 1) $ 
    399402          ELSE divi = total(echelle, 1) 
     
    415418        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 
    416419          AND ny NE 1 THEN BEGIN 
    417           IF msknan[0] EQ -1 THEN BEGIN  
     420          IF msknan[0] EQ -1 THEN BEGIN 
    418421            msknan = replicate(1b, nx, ny, nz) 
    419422            nan = 1 
     
    424427          res[temporary(bottom)] = !values.f_nan 
    425428        ENDIF 
    426         if keyword_set(integration) then divi = 1 ELSE begin  
     429        if keyword_set(integration) then divi = 1 ELSE begin 
    427430          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 2) $ 
    428431          ELSE divi = total(echelle, 2) 
     
    446449        echelle = (temporary(e33)*temporary(mask))[*]#replicate(1, jpt) 
    447450        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 
    449452          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 3) $ 
    450453          ELSE divi = total(echelle, 3) 
     
    467470        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 
    468471          AND nx*ny NE 1 THEN BEGIN 
    469           IF msknan[0] EQ -1 THEN BEGIN  
     472          IF msknan[0] EQ -1 THEN BEGIN 
    470473            msknan = replicate(1b, nx, ny, nz) 
    471474            nan = 1 
     
    476479          res[temporary(bottom)] = !values.f_nan 
    477480        ENDIF 
    478         if keyword_set(integration) then divi = 1 ELSE begin  
     481        if keyword_set(integration) then divi = 1 ELSE begin 
    479482          IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 1) $ 
    480483          ELSE divi = total(total(echelle, 1), 1) 
     
    497500        echelle = (temporary(e133[*])*temporary(mask[*]))#replicate(1, jpt) 
    498501        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 
    500503          IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 2) $ 
    501504          ELSE divi = total(total(echelle, 1), 2) 
     
    518521        echelle = (temporary(e233[*])*temporary(mask[*]))#replicate(1, jpt) 
    519522        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 
    521524          IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 2), 2) $ 
    522525          ELSE divi = total(total(echelle, 2), 2) 
     
    539542        echelle = (temporary(e1233[*])*temporary(mask[*]))#replicate(1, jpt) 
    540543        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 
    542545          IF msknan[0] NE -1 THEN divi = total(total(total(echelle*msknan, 1), 1), 1) $ 
    543546          ELSE divi = total(total(total(echelle, 1), 1), 1) 
     
    562565    tdim = size(res, /n_dimensions) 
    563566    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 
    566569        testnan = testnan < 1 
    567570        testnan = total(temporary(testnan), tdim) 
     
    569572      ENDIF ELSE divi = jpt 
    570573      res = total(res, tdim, nan = nan)/(1 > divi) 
    571     ENDELSE  
     574    ENDELSE 
    572575  ENDIF 
    573576;------------------------------------------------------------ 
     
    580583  valmask = 1e+20 
    581584  terre = where(divi EQ 0) 
    582   IF terre[0] NE -1 THEN BEGIN  
     585  IF terre[0] NE -1 THEN BEGIN 
    583586    res[temporary(terre)] = 1e+20 
    584   ENDIF  
     587  ENDIF 
    585588;------------------------------------------------------------ 
    586589; IV.2) We replace, when nan equal 1, !values.f_nan by nan 
    587590;------------------------------------------------------------ 
    588   if keyword_set(nan) NE 0 then BEGIN  
     591  if keyword_set(nan) NE 0 then BEGIN 
    589592    puttonan = where(temporary(testnan) EQ 0) 
    590593    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 
    592595      notanumber = where(finite(res) eq 0) 
    593596      if notanumber[0] NE -1 then res[temporary(notanumber)] = nan 
     
    599602  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
    600603;------------------------------------------------------------ 
    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 
    602605  return, res 
    603606;------------------------------------------------------------ 
Note: See TracChangeset for help on using the changeset viewer.