Ignore:
Timestamp:
05/02/06 14:59:12 (18 years ago)
Author:
pinsard
Message:

upgrade of CALCULS according to cerbere.lodyc.jussieu.fr: /usr/home/smasson/SAXO_RD/ : files

File:
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/ToBeReviewed/CALCULS/moyenne.pro

    r23 r25  
    1111; CATEGORY: 
    1212; 
    13 ; CALLING SEQUENCE: result = moyenne(tab,'direc',BOITE=boite) 
     13; CALLING SEQUENCE: result = moyenne(tab,'direc',BOXZOOM=boxzoom) 
    1414; 
    1515; INPUTS:       tab = 2 or 3d field 
    … …  
    1818; KEYWORD PARAMETERS: 
    1919; 
    20 ;               BOITE = [xmin,xmax,ymin,ymax (,(zmin,)zmax)] pour plus 
     20;               BOXZOOM = [xmin,xmax,ymin,ymax (,(zmin,)zmax)] pour plus 
    2121;               de detail cf domdef. 
    22 ;               boite peut prendre 5 formes:  
    23 ;               [prof2], [prof1, prof2],[lon1, lon2, lat1, lat2],   
    24 ;               [lon1, lon2, lat1, lat2, prof2],[lon1, lon2, lat1, 
    25 ;               lat2, prof1,prof2] 
     22;               boxzoom peut prendre 5 formes:  
     23;               [vert2], [vert1, vert2],[lon1, lon2, lat1, lat2],   
     24;               [lon1, lon2, lat1, lat2, vert2],[lon1, lon2, lat1, 
     25;               lat2, vert1,vert2] 
    2626; 
    2727;               NAN: not a number, a activer si l''on peut faire veut 
    … …  
    3838; 
    3939;               NODOMDEF: activer si l''on ne veut pas passer ds 
    40 ;               domdef bien que le mot cle boite soit present (comme 
     40;               domdef bien que le mot cle boxzoom soit present (comme 
    4141;               c''est le cas qd moyenne est appelee via checkfield)  
    4242;             
    4343;               INTEGRATION: pour faire une integrale plutot qu''une moyenne 
     44;  
     45;         /WDEPTH: to specify that the field is at W depth instad of T  
     46;         depth (automatically activated if vargrid eq 'W') 
    4447;             
    4548; OUTPUTS:  result:un tableau  
    … …  
    8891;------------------------------------------------------------ 
    8992;------------------------------------------------------------ 
    90 function moyenne ,tab,direc,BOITE=boite, INTEGRATION=integration $ 
    91                   , NAN = nan, NODOMDEF = nodomdef, _extra = ex 
    92 @common 
    93    tempsun = systime(1)         ; pour key_performance 
     93function moyenne, tab, direc, BOXZOOM = boxzoom, INTEGRATION = integration $ 
     94                  , NAN = nan, NODOMDEF = nodomdef, WDEPTH = wdepth $ 
     95                  , _extra = ex 
     96;--------------------------------------------------------- 
     97@cm_4mesh 
     98@cm_4data 
     99@cm_4cal 
     100  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     101@updatenew 
     102@updatekwd 
     103  ENDIF 
     104;--------------------------------------------------------- 
     105  tempsun = systime(1)          ; pour key_performance 
    94106;------------------------------------------------------------ 
    95107;I) preliminaires 
    96108;------------------------------------------------------------ 
    97    dirt = 0 
    98    dirx = 0 
    99    diry = 0 
    100    dirz = 0 
     109  dirt = 0 
     110  dirx = 0 
     111  diry = 0 
     112  dirz = 0 
    101113;------------------------------------------------------------ 
    102114; I.1) direction(s) suivants lesquelles on integre 
    103115;------------------------------------------------------------ 
    104    if ( strpos(direc,'t') ge 0 ) then dirt = 1 
    105    if ( strpos(direc,'x') ge 0 ) then dirx = 1 
    106    if ( strpos(direc,'y') ge 0 ) then diry = 1 
    107    if ( strpos(direc,'z') ge 0 ) then dirz = 1 
    108    if (dirx eq 0 and diry eq 0 and dirz eq 0) then BEGIN  
    109       IF keyword_set(bavard) then $ 
    110        IF dirt NE 1 THEN ras = report('MOYENNE: aucune valeur de direc ne convient, le champ reste inchange') 
    111       return, tab 
    112    ENDIF  
     116  if ( strpos(direc, 't') ge 0 ) then dirt = 1 
     117  if ( strpos(direc, 'x') ge 0 ) then dirx = 1 
     118  if ( strpos(direc, 'y') ge 0 ) then diry = 1 
     119  if ( strpos(direc, 'z') ge 0 ) then dirz = 1 
     120  if (dirx eq 0 and diry eq 0 and dirz eq 0) then return, tab 
    113121;------------------------------------------------------------ 
    114122; I.2) verification de la taille du tableau d'entree 
    115123;------------------------------------------------------------ 
    116    taille = size(tab) 
    117    case 1 of 
    118       taille[0] eq 1 :dim = '1d' 
    119       taille[0] eq 2 :BEGIN 
    120          dim = '2d' 
    121          if dirx eq 0 and diry eq 0 then return, tab 
    122       END 
    123       taille[0] eq 3 :BEGIN 
    124          dim = '3d' 
    125          if dirx eq 0 and diry eq 0 and dirz eq 0 then return, tab 
    126       END 
    127       else : return, report('Le tableau d''entree doit etre a 2 ou 3 dimensions s''il contient une dim temporelle utiliser grossemoyenne') 
    128    endcase 
     124  taille = size(tab) 
     125  case 1 of 
     126    taille[0] eq 1 :dim = '1d' 
     127    taille[0] eq 2 :BEGIN 
     128      dim = '2d' 
     129      if dirx eq 0 and diry eq 0 then return, tab 
     130    END 
     131    taille[0] eq 3 :BEGIN 
     132      dim = '3d' 
     133      if dirx eq 0 and diry eq 0 and dirz eq 0 then return, tab 
     134    END 
     135    else : return, report('Le tableau d''entree doit etre a 2 ou 3 dimensions s''il contient une dim temporelle utiliser grossemoyenne') 
     136  endcase 
    129137;------------------------------------------------------------ 
    130138;   I.3) obtention des facteurs d''echelles et du masque sur le sous 
    131139;   domaine concerne par la moyenne  
    132 ; redefinition du domaine ajuste a boite (a 6 elements) 
     140; redefinition du domaine ajuste a boxzoom (a 6 elements) 
    133141; ceci va nous permetre de faire les calcules que sur le sous domaine 
    134142; comcerne par la moyenne. domdef, suivit de grille nous donne tous 
    135143; les tableaux de la grille sur le sous domaine  
    136144;------------------------------------------------------------ 
    137    if keyword_set(boite) then BEGIN  
    138       Case 1 Of 
    139          N_Elements(Boite) Eq 1:bte=[lon1, lon2, lat1, lat2, 0.,boite[0]] 
    140          N_Elements(Boite) Eq 2:bte=[lon1, lon2, lat1, lat2, boite[0],boite[1]] 
    141          N_Elements(Boite) Eq 4:bte=[Boite, prof1, prof2] 
    142          N_Elements(Boite) Eq 5:bte=[Boite[0:3], 0, Boite[4]] 
    143          N_Elements(Boite) Eq 6:bte=Boite 
    144          Else: return, report('Mauvaise Definition de Boite') 
    145       endcase 
    146       oldboite = [lon1, lon2, lat1, lat2, prof1, prof2] 
    147       if NOT keyword_set(nodomdef) then domdef, bte,GRILLE=vargrid, _extra = ex 
    148    ENDIF  
     145  if keyword_set(boxzoom) then BEGIN  
     146    Case 1 Of 
     147      N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]] 
     148      N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]] 
     149      N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2] 
     150      N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]] 
     151      N_Elements(Boxzoom) Eq 6:bte = Boxzoom 
     152      Else: return, report('Mauvaise Definition de Boxzoom') 
     153    endcase 
     154    if NOT keyword_set(nodomdef) then BEGIN  
     155      savedbox = 1b 
     156      saveboxparam, 'boxparam4moyenne.dat' 
     157      domdef, bte, GRIDTYPE = vargrid, _extra = ex 
     158    ENDIF  
     159  ENDIF  
    149160;--------------------------------------------------------------- 
    150161; attribution du mask et des tableaux de longitude et latitude... 
    151162;--------------------------------------------------------------- 
    152    grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz,e1,e2,e3 
     163  IF vargrid EQ 'W' THEN wdepth = 1 
     164  grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, WDEPTH = wdepth 
    153165;------------------------------------------------------------ 
    154166;------------------------------------------------------------ 
    … …  
    156168;------------------------------------------------------------ 
    157169;------------------------------------------------------------ 
    158    if dim EQ '1d' then BEGIN 
    159       if n_elements(tab) NE nx*ny AND n_elements(tab) NE nx*ny*nz then $ 
    160        return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    161       case 1 of 
    162          nx EQ 1 AND ny EQ 1:BEGIN ;vecteur suivant z 
    163             case n_elements(tab) of 
    164                jpk:res = tab[premierz:dernierz] 
    165                nz:res = tab 
    166                ELSE:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    167             ENDCASE 
    168             if dirz EQ 1 then BEGIN  
    169                dim = '3d'  
    170                taille = size(reform(res, nx, ny, nz)) 
    171             ENDIF ELSE return, res 
    172          END 
    173          ny EQ 1:BEGIN          ;vecteur suivant x 
    174             case n_elements(tab) of 
    175                jpi:res = tab[premierx:dernierx] 
    176                nx:res = tab 
    177                ELSE:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    178             ENDCASE 
    179             if dirx EQ 1 then BEGIN  
    180                dim = '2d'  
    181                taille = size(reform(res, nx, ny)) 
    182             ENDIF ELSE return, res 
    183          END 
    184          nx EQ 1:BEGIN          ;vecteur suivant y 
    185             case n_elements(tab) of 
    186                jpj:res = tab[premiery:derniery] 
    187                ny:res = tab 
    188                ELSE:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    189             ENDCASE 
    190             if diry EQ 1 then BEGIN 
    191                dim = '2d' 
    192                taille = size(reform(res, nx, ny)) 
    193             ENDIF ELSE return, res 
    194          END 
    195       endcase 
    196    endif 
     170  if dim EQ '1d' then BEGIN 
     171    if n_elements(tab) NE nx*ny AND n_elements(tab) NE nx*ny*nz then BEGIN 
     172      if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 
     173      return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.') 
     174    ENDIF 
     175    case 1 of 
     176      nx EQ 1 AND ny EQ 1:BEGIN ;vecteur suivant z 
     177        case n_elements(tab) of 
     178          jpk:res = tab[firstz:lastz] 
     179          nz:res = tab 
     180          ELSE:BEGIN  
     181            if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 
     182            return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.') 
     183          END 
     184        ENDCASE 
     185        if dirz EQ 1 then BEGIN  
     186          dim = '3d'  
     187          taille = size(reform(res, nx, ny, nz)) 
     188        ENDIF ELSE BEGIN 
     189          if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 
     190          return, res 
     191        ENDELSE  
     192      END 
     193      ny EQ 1:BEGIN             ;vecteur suivant x 
     194        case n_elements(tab) of 
     195          jpi:res = tab[firstx:lastx] 
     196          nx:res = tab 
     197          ELSE:BEGIN  
     198            if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 
     199            return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.') 
     200          END 
     201        ENDCASE 
     202        if dirx EQ 1 then BEGIN  
     203          dim = '2d'  
     204          taille = size(reform(res, nx, ny)) 
     205        ENDIF ELSE BEGIN  
     206          if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 
     207          return, res 
     208        ENDELSE 
     209      END 
     210      nx EQ 1:BEGIN             ;vecteur suivant y 
     211        case n_elements(tab) of 
     212          jpj:res = tab[firsty:lasty] 
     213          ny:res = tab 
     214          ELSE:BEGIN  
     215            if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 
     216            return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.') 
     217          END 
     218        ENDCASE 
     219        if diry EQ 1 then BEGIN 
     220          dim = '2d' 
     221          taille = size(reform(res, nx, ny)) 
     222        ENDIF ELSE BEGIN  
     223          if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 
     224          return, res 
     225        ENDELSE  
     226      END 
     227    endcase 
     228  endif 
    197229;------------------------------------------------------------ 
    198230;------------------------------------------------------------ 
    … …  
    200232;------------------------------------------------------------ 
    201233;------------------------------------------------------------ 
    202    if (dim eq '2d') then begin 
     234  if (dim eq '2d') then begin 
    203235;--------------------------------------------------------------- 
    204236;   II.1) verification de la coherence de la taille du tableau a 
    … …  
    209241; (jpi,jpj) soit celle du domaine reduit (nx,ny) 
    210242;--------------------------------------------------------------- 
    211       case 1 of 
    212          taille[1] eq jpi and taille[2] eq jpj: $ 
    213           res = tab[premierx:dernierx, premiery:derniery] 
    214          taille[1] eq  nx and taille[2] eq  ny:res = tab 
    215          else:return, report('Probleme d''adequation entre les tailles du domaine nx*ny '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)) 
    216       ENDCASE 
    217       if keyword_set(nan) NE 0 then BEGIN  
    218          if nan NE 1 then BEGIN ; si nan n''est pas !values.f_nan  
     243    case 1 of 
     244      taille[1] eq jpi and taille[2] eq jpj: $ 
     245        res = tab[firstx:lastx, firsty:lasty] 
     246      taille[1] eq  nx and taille[2] eq  ny:res = tab 
     247      else:BEGIN  
     248        if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 
     249        return, report('Probleme d''adequation entre les tailles du domaine nx*ny '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)) 
     250      END 
     251    ENDCASE 
     252    if keyword_set(nan) NE 0 then BEGIN  
     253      if nan NE 1 then BEGIN    ; si nan n''est pas !values.f_nan  
    219254; on le met a !values.f_nan 
    220             if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 
    221             ELSE notanumber = where(abs(res) GT abs(nan)/10.) 
    222             if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 
    223          ENDIF 
     255        if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 
     256        ELSE notanumber = where(abs(res) GT abs(nan)/10.) 
     257        if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 
    224258      ENDIF 
     259    ENDIF 
    225260;--------------------------------------------------------------- 
    226261; rq IL FAUT FAIRE ATTENTION AUX CAS OU LA DIM A MOYENNER = 1, ET 
    … …  
    228263; PEUVENT SEMBLER INUTILE AU DEPART 
    229264;--------------------------------------------------------------- 
    230       if nx EQ 1 OR ny EQ 1 then BEGIN  
    231          res = reform(res, nx, ny, /over) 
    232          mask =  reform(mask, nx, ny, nz, /over) 
    233          e1 =  reform(e1, nx, ny, /over) 
    234          e2 =  reform(e2, nx, ny, /over) 
    235       endif 
     265    if nx EQ 1 OR ny EQ 1 then BEGIN  
     266      res = reform(res, nx, ny, /over) 
     267      e1 =  reform(e1, nx, ny, /over) 
     268      e2 =  reform(e2, nx, ny, /over) 
     269    endif 
     270    if nx EQ 1 OR ny EQ 1 OR nz EQ 1 THEN $ 
     271      mask =  reform(mask, nx, ny, nz, /over) 
    236272;--------------------------------------------------------------- 
    237273; II.3) differents types de moyennes 
    238274;--------------------------------------------------------------- 
    239       if nz eq jpk then mask = mask[*, *,niveau-1] ELSE mask = mask[*, *,nz-1] 
    240       if keyword_set(nan) NE 0 then begin 
    241          msknan = replicate(1., nx, ny) 
    242          notanumber = where(finite(res) EQ 0) 
    243          if notanumber[0] NE -1 then msknan[notanumber] = !values.f_nan 
    244       ENDIF ELSE msknan = 1 
    245       case 1 of 
    246          (dirx eq 1) and (diry eq 0) : begin 
    247             e = e1*mask 
    248             if keyword_set(integration) then divi=1 $ 
    249             else begin  
    250                divi = e*msknan 
    251                if ny EQ 1 then divi = reform(divi,nx,ny, /over) 
    252                divi = total(divi,1, nan = nan) 
    253             endelse 
    254             res = res*e 
    255             if ny EQ 1 then res = reform(res,nx,ny, /over) 
    256             res = total(res,1, nan = nan)/(divi > 1.) 
    257             if keyword_set(nan) then begin 
    258                testnan  = finite(msknan)+(1-mask) 
    259                if ny EQ 1 then testnan  = reform(testnan,nx,ny, /over) 
    260                testnan = total(testnan,1) 
    261             endif 
    262          end 
    263          (dirx eq 0) and (diry eq 1) : begin 
    264             e = e2*mask 
    265             if keyword_set(integration) then divi=1 $ 
    266             else begin  
    267                divi = e*msknan 
    268                if ny EQ 1 then divi = reform(divi,nx,ny, /over) 
    269                divi = total(divi,2, nan = nan) 
    270             endelse 
    271             res = res*e 
    272             if ny EQ 1 then res = reform(res,nx,ny, /over) 
    273             res = total(res,2, nan = nan)/(divi > 1.) 
    274             if keyword_set(nan) then begin 
    275                testnan  = finite(msknan)+(1-mask) 
    276                if ny EQ 1 then testnan  = reform(testnan,nx,ny, /over) 
    277                testnan = total(testnan,2) 
    278             endif 
    279          end 
    280          (dirx eq 1) and (diry eq 1) : begin 
    281             if keyword_set(integration) then divi=1 else divi = total(e1*e2*mask*msknan, nan = nan) 
    282             res = total(res*e1*e2*mask, nan = nan)/(divi > 1.) 
    283             if keyword_set(nan) then begin 
    284                testnan  = finite(msknan)+(1-mask) 
    285                testnan = total(testnan) 
    286             endif 
    287          end 
    288       endcase 
    289    endif 
     275    mask = mask[*, *, 0] 
     276    if keyword_set(nan) NE 0 then msknan = finite(res) ELSE msknan = -1 
     277    case 1 of 
     278      (dirx eq 1) and (diry eq 0) : begin 
     279        e = e1*mask 
     280        if keyword_set(integration) then divi = 1 $ 
     281        else begin  
     282          divi = e 
     283          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan 
     284          if ny EQ 1 then divi = reform(divi, nx, ny, /over) 
     285          divi = total(divi, 1) 
     286        endelse 
     287        res = res*e 
     288        if ny EQ 1 then res = reform(res, nx, ny, /over) 
     289        res = total(res, 1, nan = nan)/(divi > 1.) 
     290        if msknan[0] NE -1 then begin 
     291          testnan  = msknan*mask 
     292          if ny EQ 1 then testnan = reform(testnan, nx, ny, /over) 
     293          testnan = total(testnan, 1)+(total(mask, 1) EQ 0) 
     294        endif 
     295      end 
     296      (dirx eq 0) and (diry eq 1) : begin 
     297        e = e2*mask 
     298        if keyword_set(integration) then divi = 1 $ 
     299        else begin  
     300          divi = e 
     301          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan 
     302          if ny EQ 1 then divi = reform(divi, nx, ny, /over) 
     303          divi = total(divi, 2) 
     304        endelse 
     305        res = res*e 
     306        if ny EQ 1 then res = reform(res, nx, ny, /over) 
     307        res = total(res, 2, nan = nan)/(divi > 1.) 
     308        if msknan[0] NE -1 then begin 
     309          testnan  = msknan*mask 
     310          if ny EQ 1 then testnan  = reform(testnan, nx, ny, /over) 
     311          testnan = total(testnan, 2)+(total(mask, 2) EQ 0) 
     312        endif 
     313      end 
     314      (dirx eq 1) and (diry eq 1) : begin 
     315        if keyword_set(integration) then divi = 1 else BEGIN  
     316          IF msknan[0] NE -1 THEN divi = total(e1*e2*mask*msknan) $ 
     317          ELSE divi = total(e1*e2*mask) 
     318        ENDELSE  
     319        res = total(res*e1*e2*mask, nan = nan)/(divi > 1.) 
     320        if msknan[0] NE -1 then begin 
     321          testnan  = msknan*mask 
     322          testnan = total(testnan)+(total(mask) EQ 0) 
     323        endif 
     324      end 
     325    endcase 
     326  endif 
    290327;------------------------------------------------------------ 
    291328;------------------------------------------------------------ 
    … …  
    293330;------------------------------------------------------------ 
    294331;------------------------------------------------------------ 
    295    if (dim eq '3d') then begin 
     332  if (dim eq '3d') then begin 
    296333;--------------------------------------------------------------- 
    297334; III.1) verification de la coherence de la taille du tableau a 
    … …  
    302339; (jpi,jpj,jpk) soit celle du domaine reduit (nx,ny,ny) 
    303340;--------------------------------------------------------------- 
    304       case 1 of 
    305          taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk: $ 
    306           res=tab[premierx:dernierx, premiery:derniery, premierz:dernierz] 
    307          taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz: $ 
    308           res=tab[premierx:dernierx, premiery:derniery, *] 
    309          taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq nz :res=tab 
    310          taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk : $ 
    311           res=tab[*, *, premierz:dernierz] 
    312          else:return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)) 
    313       endcase 
    314       if keyword_set(nan) NE 0 then BEGIN  
    315          if nan NE 1 then BEGIN ; si nan n''est pas !values.f_nan  
     341    case 1 of 
     342      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk: $ 
     343        res = tab[firstx:lastx, firsty:lasty, firstz:lastz] 
     344      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz: $ 
     345        res = tab[firstx:lastx, firsty:lasty, *] 
     346      taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq nz :res = tab 
     347      taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk : $ 
     348        res = tab[*, *, firstz:lastz] 
     349      else:BEGIN  
     350        if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 
     351        return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)) 
     352      END 
     353    endcase 
     354    if keyword_set(nan) NE 0 then BEGIN  
     355      if nan NE 1 then BEGIN    ; si nan n''est pas !values.f_nan  
    316356; on le met a !values.f_nan 
    317             if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 
    318             ELSE notanumber = where(abs(res) GT abs(nan)/10.) 
    319             if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 
    320          ENDIF 
     357        if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 
     358        ELSE notanumber = where(abs(res) GT abs(nan)/10.) 
     359        if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 
    321360      ENDIF 
     361    ENDIF 
    322362;--------------------------------------------------------------- 
    323363; rq IL FAUT FAIRE ATTENTION AUX CAS OU LA DIM A MOYENNER = 1, ET 
    … …  
    325365; PEUVENT SEMBLER INUTILE AU DEPART 
    326366;--------------------------------------------------------------- 
    327       if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then BEGIN  
    328          res = reform(res, nx, ny, nz, /over) 
    329          mask =  reform(mask, nx, ny, nz, /over) 
    330          e1 =  reform(e1, nx, ny, /over) 
    331          e2 =  reform(e2, nx, ny, /over) 
    332       endif 
     367    if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then BEGIN  
     368      res = reform(res, nx, ny, nz, /over) 
     369      e1 =  reform(e1, nx, ny, /over) 
     370      e2 =  reform(e2, nx, ny, /over) 
     371    endif 
     372    if nx EQ 1 OR ny EQ 1 OR nz EQ 1 THEN $ 
     373      mask =  reform(mask, nx, ny, nz, /over) 
     374    IF keyword_set(key_partialstep) THEN BEGIN 
     375; the top of the ocean floor is 
     376      IF vargrid EQ 'T' OR vargrid EQ 'W' THEN bottom = total(mask, 3) $ 
     377      ELSE bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)  
     378; we suppress columns with only ocean or land  
     379      good = where(bottom NE 0 AND bottom NE nz) 
     380; the bottom of the ocean in 3D index is: 
     381      bottom = lindgen(nx*ny)+(temporary(bottom)-1L)*nx*ny 
     382      IF good[0] NE -1 THEN bottom = bottom[good] $ 
     383      ELSE bottom = -1 
     384    ENDIF ELSE bottom = -1 
    333385;--------------------------------------------------------------- 
    334386; III.2) differents types de moyennes 
    335387;--------------------------------------------------------------- 
    336       if keyword_set(nan) NE 0 then begin 
    337          msknan = replicate(1., nx, ny, nz) 
    338          notanumber = where(finite(res) EQ 0) 
    339          if notanumber[0] NE -1 then msknan[notanumber] = !values.f_nan 
    340       ENDIF ELSE msknan = 1 
    341       case 1 of 
    342          (dirx eq 1) and (diry eq 0) and (dirz eq 0) : begin 
    343             e13 = e1[*]#replicate(1.,nz) 
    344             e13 = reform(e13,nx,ny,nz, /over) 
    345             if keyword_set(integration) then divi = 1 else begin 
    346                divi = e13*mask*msknan 
    347                if nz EQ 1 then divi = reform(divi,nx,ny,nz, /over) 
    348                divi = total(divi,1, nan = nan) 
    349             endelse 
    350             res = res*e13*mask 
    351             if nz EQ 1 then res = reform(res,nx,ny,nz, /over) 
    352             res = total(res,1, nan = nan)/(divi > 1.) 
    353             e13 = 1 
    354             if keyword_set(nan) then begin 
    355                testnan  = finite(msknan)+(1-mask) 
    356                if nz EQ 1 then testnan  = reform(testnan,nx,ny,nz, /over) 
    357                testnan = total(testnan,1) 
    358             endif 
    359          end 
    360          (dirx eq 0) and (diry eq 1) and (dirz eq 0) : begin 
    361             e23 = e2[*]#replicate(1.,nz) 
    362             e23 = reform(e23,nx,ny,nz, /over) 
    363             if keyword_set(integration) then divi =1 else begin 
    364                divi = e23*mask*msknan 
    365                if nz EQ 1 then divi = reform(divi,nx,ny,nz, /over) 
    366                divi = total(divi,2, nan = nan) 
    367             endelse 
    368             res = res*e23*mask 
    369             if nz EQ 1 then res = reform(res,nx,ny,nz, /over) 
    370             res = total(res,2, nan = nan)/(divi > 1.) 
    371             e23 = 1 
    372             if keyword_set(nan) then begin 
    373                testnan  = finite(msknan)+(1-mask) 
    374                if nz EQ 1 then testnan  = reform(testnan,nx,ny,nz, /over) 
    375                testnan = total(testnan,2) 
    376             endif 
    377          end 
    378          (dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin 
    379             e33 = replicate(1,1.*nx*ny)#e3 
    380             e33 = reform(e33,nx,ny,nz, /over) 
    381             if keyword_set(integration) then divi =1 else begin 
    382                divi = e33*mask*msknan 
    383                if nz EQ 1 then divi = reform(divi,nx,ny,nz, /over) 
    384                divi = total(divi,3, nan = nan) 
    385             endelse 
    386             res = res*e33*mask 
    387             if nz EQ 1 then res = reform(res,nx,ny,nz, /over) 
    388             res = total(res,3, nan = nan)/(divi > 1.) 
    389             e33 = 1 
    390             if keyword_set(nan) then begin 
    391                testnan  = finite(msknan)+(1-mask) 
    392                if nz EQ 1 then testnan  = reform(testnan,nx,ny,nz, /over) 
    393                testnan = total(testnan,3) 
    394             endif 
    395          end 
    396          (dirx eq 1) and (diry eq 1) and (dirz eq 0) : begin 
    397             e123 = (e1*e2)[*]#replicate(1.,nz) 
    398             e123 = reform(e123,nx,ny,nz, /over) 
    399             if keyword_set(integration) then divi =1 else $ 
    400              divi = e123*mask*msknan 
    401             if nz EQ 1 then divi = reform(divi,nx,ny,nz, /over) 
    402             divi = total(total(divi,1, nan = nan),1, nan = nan) 
    403             res = res*e123*mask 
    404             if nz EQ 1 then res = reform(res,nx,ny,nz, /over) 
    405             res = total(total(res,1, nan = nan),1, nan = nan) / (divi > 1.) 
    406             e123 = 1 
    407             if keyword_set(nan) then begin 
    408                testnan  = finite(msknan)+(1-mask) 
    409                if nz EQ 1 then testnan  = reform(testnan,nx,ny,nz, /over) 
    410                testnan = total(total(testnan,1),1) 
    411             endif 
    412          end 
    413          (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin 
    414             e133 = e1[*]#e3 
    415             e133 = reform(e133,nx,ny,nz, /over) 
    416             divi = e133*mask*msknan 
    417             if nz EQ 1 then divi = reform(divi,nx,ny,nz, /over) 
    418             if keyword_set(integration) then divi =1 else $ 
    419              divi = total(total(divi,1, nan = nan),2, nan = nan) 
    420             res = res*e133*mask 
    421             if nz EQ 1 then res = reform(res,nx,ny,nz, /over) 
    422             res = total(total(res,1, nan = nan),2, nan = nan) / (divi > 1.) 
    423             e133 = 1 
    424             if keyword_set(nan) then begin 
    425                testnan  = finite(msknan)+(1-mask) 
    426                if nz EQ 1 then testnan  = reform(testnan,nx,ny,nz, /over) 
    427                testnan = total(total(testnan,1),2) 
    428             endif 
    429          end 
    430          (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin 
    431             e233 = e2[*]#e3 
    432             e233 = reform(e233,nx,ny,nz, /over) 
    433             divi = e233*mask*msknan 
    434             if nz EQ 1 then divi = reform(divi,nx,ny,nz, /over) 
    435             if keyword_set(integration) then divi =1 else $ 
    436              divi = total(total(divi,2, nan = nan),2, nan = nan) 
    437             res = res*e233*mask 
    438             if nz EQ 1 then res = reform(res,nx,ny,nz, /over) 
    439             res = total(total(res,2, nan = nan),2, nan = nan) / (divi > 1.) 
    440             e233 = 1 
    441             if keyword_set(nan) then begin 
    442                testnan  = finite(msknan)+(1-mask) 
    443                if nz EQ 1 then testnan  = reform(testnan,nx,ny,nz, /over) 
    444                testnan = total(total(testnan,2),2) 
    445             endif 
    446          end 
    447          (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin 
    448             e1233 = (e1*e2)[*]#e3 
    449             e1233 = reform(e1233,nx,ny,nz, /over) 
    450             if keyword_set(integration) then divi =1 else divi = total(e1233*mask*msknan, nan = nan) 
    451             res = total(res*e1233*mask, nan = nan) / (divi > 1.) 
    452             e1233 = 1 
    453             if keyword_set(nan) then begin 
    454                testnan  = finite(msknan)+(1-mask) 
    455                testnan = total(testnan) 
    456             endif 
    457          end 
    458       endcase 
    459    endif 
     388    if keyword_set(nan) NE 0 then msknan = finite(res) ELSE msknan = -1 
     389    case 1 of 
     390      (dirx eq 1) and (diry eq 0) and (dirz eq 0) : begin 
     391        e13 = e1[*]#replicate(1., nz) 
     392        e13 = reform(e13, nx, ny, nz, /over) 
     393        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 
     394          AND nx NE 1 THEN BEGIN 
     395          IF msknan[0] EQ -1 THEN BEGIN  
     396            msknan = replicate(1b, nx, ny, nz) 
     397            nan = 1 
     398          endif 
     399          msknan[bottom] = 0 
     400          res[bottom] = !values.f_nan 
     401        ENDIF 
     402        if keyword_set(integration) then divi = 1 else begin 
     403          divi = e13*mask 
     404          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan 
     405          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over) 
     406          divi = total(divi, 1) 
     407        ENDELSE 
     408        res = res*e13*mask 
     409        if nz EQ 1 then res = reform(res, nx, ny, nz, /over) 
     410        res = total(res, 1, nan = nan)/(divi > 1.) 
     411        e13 = 1 
     412        if msknan[0] NE -1 then begin 
     413          testnan  = msknan*mask 
     414          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over) 
     415          testnan = total(testnan, 1)+(total(mask, 1) EQ 0) 
     416        endif 
     417      end 
     418      (dirx eq 0) and (diry eq 1) and (dirz eq 0) : begin 
     419        e23 = e2[*]#replicate(1., nz) 
     420        e23 = reform(e23, nx, ny, nz, /over) 
     421        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 
     422          AND ny NE 1 THEN BEGIN 
     423          IF msknan[0] EQ -1 THEN BEGIN  
     424            msknan = replicate(1b, nx, ny, nz) 
     425            nan = 1 
     426          endif 
     427          msknan[bottom] = 0 
     428          res[bottom] = !values.f_nan 
     429        ENDIF 
     430        if keyword_set(integration) then divi = 1 else begin 
     431          divi = e23*mask 
     432          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan 
     433          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over) 
     434          divi = total(divi, 2) 
     435        ENDELSE 
     436        res = res*e23*mask 
     437        if nz EQ 1 then res = reform(res, nx, ny, nz, /over) 
     438        res = total(res, 2, nan = nan)/(divi > 1.) 
     439        e23 = 1 
     440        if msknan[0] NE -1 then begin 
     441          testnan  = msknan*mask 
     442          if nz EQ 1 then testnan = reform(testnan, nx, ny, nz, /over) 
     443          testnan = total(testnan, 2)+(total(mask, 2) EQ 0) 
     444        endif 
     445      end 
     446      (dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin 
     447        e33 = replicate(1, 1.*nx*ny)#e3 
     448        e33 = reform(e33, nx, ny, nz, /over) 
     449        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
     450          IF keyword_set(wdepth) THEN $ 
     451            e33[bottom] = (e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
     452          ELSE e33[bottom] = (e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
     453        ENDIF 
     454        if keyword_set(integration) then divi = 1 else begin 
     455          divi = e33*mask 
     456          if msknan[0] NE -1 then divi = temporary(divi)*msknan 
     457          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over) 
     458          divi = total(divi, 3) 
     459        ENDELSE 
     460        res = res*e33*mask 
     461        if nz EQ 1 then res = reform(res, nx, ny, nz, /over) 
     462        res = total(res, 3, nan = nan)/(divi > 1.) 
     463        e33 = 1 
     464        if msknan[0] NE -1 then begin 
     465          testnan  = msknan*mask 
     466          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over) 
     467          testnan = total(testnan, 3)+(total(mask, 3) EQ 0) 
     468        endif 
     469      end 
     470      (dirx eq 1) and (diry eq 1) and (dirz eq 0) : begin 
     471        e123 = (e1*e2)[*]#replicate(1., nz) 
     472        e123 = reform(e123, nx, ny, nz, /over) 
     473        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 
     474          AND nx*ny NE 1 THEN BEGIN 
     475          IF msknan[0] EQ -1 THEN BEGIN  
     476            msknan = replicate(1b, nx, ny, nz) 
     477            nan = 1 
     478          endif 
     479          msknan[bottom] = 0 
     480          res[bottom] = !values.f_nan 
     481        ENDIF 
     482        if keyword_set(integration) then divi = 1 else BEGIN  
     483          divi = e123*mask 
     484          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan 
     485          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over) 
     486          divi = total(total(divi, 1), 1) 
     487        ENDELSE 
     488        res = res*e123*mask 
     489        if nz EQ 1 then res = reform(res, nx, ny, nz, /over) 
     490        res = total(total(res, 1, nan = nan), 1, nan = nan) / (divi > 1.) 
     491        e123 = 1 
     492        if msknan[0] NE -1 then begin 
     493          testnan  = msknan*mask 
     494          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over) 
     495          testnan = total(total(testnan, 1), 1)+(total(total(mask, 1), 1) EQ 0) 
     496        endif 
     497      end 
     498      (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin 
     499        e133 = e1[*]#e3 
     500        e133 = reform(e133, nx, ny, nz, /over) 
     501        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
     502          IF keyword_set(wdepth) THEN $ 
     503            e133[bottom] = (e1*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
     504          ELSE e133[bottom] = (e1*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
     505        ENDIF 
     506        if keyword_set(integration) then divi = 1 else BEGIN  
     507          divi = e133*mask 
     508          if msknan[0] NE -1 then divi = temporary(divi)*msknan 
     509          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over) 
     510          divi = total(total(divi, 1), 2) 
     511        ENDELSE 
     512        res = res*e133*mask 
     513        if nz EQ 1 then res = reform(res, nx, ny, nz, /over) 
     514        res = total(total(res, 1, nan = nan), 2, nan = nan) / (divi > 1.) 
     515        e133 = 1 
     516        if msknan[0] NE -1 then begin 
     517          testnan  = msknan*mask 
     518          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over) 
     519          testnan = total(total(testnan, 1), 2)+(total(total(mask, 1), 2) EQ 0) 
     520        endif 
     521      end 
     522      (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin 
     523        e233 = e2[*]#e3 
     524        e233 = reform(e233, nx, ny, nz, /over) 
     525        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
     526          IF keyword_set(wdepth) THEN $ 
     527            e233[bottom] = (e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
     528          ELSE e233[bottom] = (e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
     529        ENDIF 
     530        if keyword_set(integration) then divi = 1 else BEGIN  
     531          divi = e233*mask 
     532          if msknan[0] NE -1 then divi = temporary(divi)*msknan 
     533          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over) 
     534          divi = total(total(divi, 2), 2) 
     535        ENDELSE 
     536        res = res*e233*mask 
     537        if nz EQ 1 then res = reform(res, nx, ny, nz, /over) 
     538        res = total(total(res, 2, nan = nan), 2, nan = nan) / (divi > 1.) 
     539        e233 = 1 
     540        if msknan[0] NE -1 then begin 
     541          testnan  = msknan*mask 
     542          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over) 
     543          testnan = total(total(testnan, 2), 2)+(total(total(mask, 2), 2) EQ 0) 
     544        endif 
     545      end 
     546      (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin 
     547        e1233 = (e1*e2)[*]#e3 
     548        e1233 = reform(e1233, nx, ny, nz, /over) 
     549        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
     550          IF keyword_set(wdepth) THEN $ 
     551            e1233[bottom] = (e1*e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
     552          ELSE e1233[bottom] = (e1*e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
     553        ENDIF 
     554        if keyword_set(integration) then divi = 1 else BEGIN  
     555          if msknan[0] NE -1 then divi = total(e1233*mask*msknan) $ 
     556          ELSE divi = total(e1233*mask) 
     557        ENDELSE 
     558        res = total(res*e1233*mask, nan = nan) / (divi > 1.) 
     559        e1233 = 1 
     560        if msknan[0] NE -1 then begin 
     561          testnan  = msknan*mask 
     562          testnan = total(testnan)+(total(mask) EQ 0) 
     563        endif 
     564      end 
     565    endcase 
     566  endif 
    460567;------------------------------------------------------------ 
    461568;------------------------------------------------------------ 
    … …  
    465572; IV.1) on masque les terres par une valeur a 1e+20 
    466573;------------------------------------------------------------ 
    467    valmask = 1e+20 
    468    terre = where(divi EQ 0) 
    469    IF terre[0] NE -1 THEN BEGIN 
    470       res(terre) = 1e+20 
    471    ENDIF   
     574  valmask = 1e+20 
     575  terre = where(divi EQ 0) 
     576  IF terre[0] NE -1 THEN BEGIN 
     577    res(terre) = 1e+20 
     578  ENDIF   
    472579;------------------------------------------------------------ 
    473580; IV.2) on remplace, quand nan ne 1, !values.f_nan par nan 
    474581;------------------------------------------------------------ 
    475    if keyword_set(nan) NE 0 then BEGIN  
    476       puttonan = where(testnan EQ 0) 
    477       if puttonan[0] NE -1 then res[puttonan] = !values.f_nan 
    478       if nan NE 1 then BEGIN  
    479          notanumber = where(finite(res) eq 0) 
    480          if notanumber[0] NE -1 then res[notanumber] = nan 
    481       ENDIF 
    482    ENDIF 
     582  if keyword_set(nan) NE 0 then BEGIN  
     583    puttonan = where(testnan EQ 0) 
     584    if puttonan[0] NE -1 then res[puttonan] = !values.f_nan 
     585    if nan NE 1 then BEGIN  
     586      notanumber = where(finite(res) eq 0) 
     587      if notanumber[0] NE -1 then res[notanumber] = nan 
     588    ENDIF 
     589  ENDIF 
    483590;------------------------------------------------------------ 
    484591; IV.3) on se remplace ds le sous domaine qui etait definit a l''entree de 
    485592; moyenne  
    486593;------------------------------------------------------------ 
    487    if keyword_set(boite) AND NOT keyword_set(nodomdef) then domdef, oldboite,GRILLE=vargrid 
    488 ;------------------------------------------------------------ 
    489    if keyword_set(key_performance) THEN print, 'temps moyenne', systime(1)-tempsun  
    490    return, res 
     594  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 
     595;------------------------------------------------------------ 
     596  if keyword_set(key_performance) THEN print, 'temps moyenne', systime(1)-tempsun  
     597  return, res 
    491598;------------------------------------------------------------ 
    492599;------------------------------------------------------------ 
Note: See TracChangeset for help on using the changeset viewer.