Changeset 25 for trunk/ToBeReviewed


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

Location:
trunk/ToBeReviewed/CALCULS
Files:
3 added
17 copied

Legend:

Unmodified
Added
Removed

  • trunk/ToBeReviewed/CALCULS/curl.pro

    r23 r25  
    6161   tempsun = systime(1)         ; pour key_performance 
    6262; 
     63   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 
     64     return, report(['This version of curl is based on Arakawa C-grid.' $ 
     65                     , 'U and V grids must therefore be defined']) 
     66; 
    6367   u = litchamp(uu) 
    6468   v = litchamp(vv) 
    … …  
    7175; on trouve les points que u et v ont en communs 
    7276;------------------------------------------------------------ 
    73    indicexu = (lindgen(jpi))[premierxu:premierxu+nxu-1] 
    74    indicexv = (lindgen(jpi))[premierxv:premierxv+nxv-1] 
     77   indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 
     78   indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 
    7579   indicex = inter(indicexu, indicexv) 
    76    indiceyu = (lindgen(jpj))[premieryu:premieryu+nyu-1] 
    77    indiceyv = (lindgen(jpj))[premieryv:premieryv+nyv-1] 
     80   indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 
     81   indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 
    7882   indicey = inter(indiceyu, indiceyv) 
    7983   nx = n_elements(indicex)  
    … …  
    9599             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    96100            case 1 of 
    97                nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
    98                nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    99                nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
    100                nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     101               nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
     102               nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     103               nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
     104               nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    101105               ELSE : 
    102106            endcase 
    … …  
    114118         coefu = (e1u[indice2d])[*]#replicate(1, nzt) 
    115119         coefu = reform(coefu, nx, ny, nzt, /over) 
    116          coefu = coefu*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] 
     120         coefu = coefu*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    117121         terreu = where(coefu EQ 0) 
    118122         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan 
    … …  
    120124         coefv = (e2v[indice2d])[*]#replicate(1, nzt) 
    121125         coefv = reform(coefv, nx, ny, nzt, /over) 
    122          coefv = coefv*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] 
     126         coefv = coefv*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    123127         terrev = where(coefv EQ 0) 
    124128         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan 
    125129 
    126          tabf = (fmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] 
     130         tabf = (fmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    127131         div = (e1f[indice2d]*e2f[indice2d])[*]#replicate(1, nzt) 
    128132         div = reform(div, nx, ny, nzt, /over) 
    … …  
    137141; mise a !values.f_nan de la bordure 
    138142;------------------------------------------------------------ 
    139          if NOT keyword_set(key_periodique)  OR nx NE jpi then begin 
     143         if NOT keyword_set(key_periodic)  OR nx NE jpi then begin 
    140144            psi(0, *, *) = !values.f_nan 
    141145            psi(nx-1, *, *) = !values.f_nan 
    … …  
    150154; pour le trace graphique 
    151155;------------------------------------------------------------ 
    152          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t', 'f'] 
    153          if keyword_set(direc) then psi = moyenne(psi,direc,/nan, boite = boite) 
     156         domdef, (glagmt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f'] 
     157         if keyword_set(direc) then psi = moyenne(psi,direc,/nan) 
    154158 
    155159;---------------------------------------------------------------------------- 
    … …  
    170174             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    171175               if nxu NE nx then $ 
    172                 if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]  
     176                if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]  
    173177               IF nxv NE nx THEN $ 
    174                 if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     178                if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    175179               IF nyu NE ny THEN $ 
    176                 if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]  
     180                if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]  
    177181               IF nyv NE ny THEN $ 
    178                 if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     182                if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    179183            END 
    180184            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    … …  
    191195; calcul du rotationnel 
    192196;---------------------------------------------------------------------------- 
    193          coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*premierzt] 
     197         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 
    194198         terreu = where(coefu EQ 0) 
    195199         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan 
    … …  
    197201         coefu = reform(coefu, nx, ny, jpt, /over) 
    198202; 
    199          coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*premierzt] 
     203         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 
    200204         terrev = where(coefv EQ 0) 
    201205         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan 
    … …  
    203207         coefv = reform(coefv, nx, ny, jpt, /over) 
    204208; 
    205          tabf = (fmask())[indice2d+jpi*jpj*premierzt]/(e1f[indice2d]*e2f[indice2d]) 
     209         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d]) 
    206210         tabf = temporary(tabf[*])#replicate(1, jpt) 
    207211         tabf = reform(tabf, nx, ny, jpt, /over) 
    … …  
    217221; mise a !values.f_nan de la bordure 
    218222;------------------------------------------------------------ 
    219          if NOT keyword_set(key_periodique) OR nx NE jpi then begin 
     223         if NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    220224            psi(0, *, *) = !values.f_nan 
    221225            psi(nx-1, *, *) = !values.f_nan 
    … …  
    227231         if terref[0] NE -1 then psi[temporary(terref)] = valmask 
    228232;---------------------------------------------------------------------------- 
    229          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t', 'f'] 
    230          if keyword_set(direc) then psi = grossemoyenne(psi,direc,/nan, boite = boite) 
     233         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f'] 
     234         if keyword_set(direc) then psi = grossemoyenne(psi,direc,/nan) 
    231235;---------------------------------------------------------------------------- 
    232236      END 
    … …  
    255259             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    256260               if nxu NE nx then $ 
    257                 if indicex[0] EQ premierxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]  
     261                if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]  
    258262               IF nxv NE nx THEN $ 
    259                 if indicex[0] EQ premierxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
     263                if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
    260264               IF nyu NE ny THEN $ 
    261                 if indicey[0] EQ premieryu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]  
     265                if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]  
    262266               IF nyv NE ny THEN $ 
    263                 if indicey[0] EQ premieryv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
     267                if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
    264268            END 
    265269            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    … …  
    273277; calcul du rotationnel 
    274278;------------------------------------------------------------ 
    275          coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*premierzt] 
     279         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 
    276280         terreu = where(coefu EQ 0) 
    277281         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan 
    278          coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*premierzt] 
     282         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 
    279283         terrev = where(coefv EQ 0) 
    280284         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan 
    281          tabf = (fmask())[indice2d+jpi*jpj*premierzt]/(e1f[indice2d]*e2f[indice2d]) 
     285         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d]) 
    282286; 
    283287         zu = u*temporary(coefu) 
    … …  
    290294; mise a !values.f_nan de la bordure 
    291295;------------------------------------------------------------ 
    292          if  NOT keyword_set(key_periodique) OR nx NE jpi then begin 
     296         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    293297            psi(0, *) = !values.f_nan 
    294298            psi(nx-1, *) = !values.f_nan 
    … …  
    303307; pour le trace graphique 
    304308;------------------------------------------------------------ 
    305          domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], prof1, prof2, grille = ['t', 'f'] 
    306          if keyword_set(direc) then psi = moyenne(psi,direc,/nan, boite = boite) 
    307  
     309         domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f'] 
     310         if keyword_set(direc) then psi = moyenne(psi,direc,/nan) 
    308311 
    309312      END 

  • trunk/ToBeReviewed/CALCULS/depth2level.pro

    r23 r25  
    6060;------------------------------------------------------------ 
    6161   in = litchamp(tab) 
    62    grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz 
    63    glam = 1 
    64    gphi = 1 
     62   grille,mask,-1,-1,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz 
    6563;--------------------------------------------------------------- 
    6664; verification de la coherence entre la taille du tableau et le domaine definit par domdef 
    … …  
    7068   if taille[0] NE 2 then return, report('le champ en entree doit contenir un tableau 2d') 
    7169   case 1 of 
    72       taille[1] eq jpi and taille[2] eq jpj:in=in[premierx:dernierx, premiery:derniery] 
     70      taille[1] eq jpi and taille[2] eq jpj:in=in[firstx:lastx, firsty:lasty] 
    7371      taille[1] eq  nx and taille[2] eq  ny: 
    7472      else:return, report('Probleme d''adequation entre les tailles du domaine et celle du champ.') 
    … …  
    8583;------------------------------------------------------------ 
    8684; on passe en tableaux qui ont la taille des tableaux 3d 
    87    prof=replicate(1,nx*ny)#gdep[premierz:dernierz]  
     85   prof=replicate(1,nx*ny)#gdep[firstz:lastz]  
    8886   in = in[*]#replicate(1, nz) 
    8987; 
    … …  
    9492   if keyword_set(upper) then begin 
    9593      levels = levels-1 
    96       notvalid = where(levels EQ 0) 
     94      notvalid = where(levels EQ -1) 
    9795   ENDIF ELSE notvalid = where(levels EQ nz) 
    9896   IF notvalid[0] NE -1 THEN levels[notvalid] = !values.f_nan 

  • trunk/ToBeReviewed/CALCULS/div.pro

    r23 r25  
    5656   tempsun = systime(1)         ; pour key_performance 
    5757@common 
    58    u = uu 
    59    v = vv 
    60  
     58; 
     59   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 
     60     return, report(['This version of div is based on Arakawa C-grid.' $ 
     61                     , 'U and V grids must therefore be defined']) 
     62; 
     63   u = litchamp(uu) 
     64   v = litchamp(vv) 
     65; 
    6166   date1 = time[0] 
    6267   if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1] 
    … …  
    6671; on trouve les points que u et v ont en communs 
    6772;------------------------------------------------------------ 
    68    indicexu = (lindgen(jpi))[premierxu:premierxu+nxu-1] 
    69    indicexv = (lindgen(jpi))[premierxv:premierxv+nxv-1] 
     73   indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 
     74   indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 
    7075   indicex = inter(indicexu, indicexv) 
    71    indiceyu = (lindgen(jpj))[premieryu:premieryu+nyu-1] 
    72    indiceyv = (lindgen(jpj))[premieryv:premieryv+nyv-1] 
     76   indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 
     77   indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 
    7378   indicey = inter(indiceyu, indiceyv) 
    7479   nx = n_elements(indicex)  
    … …  
    8691;------------------------------------------------------------ 
    8792         case 1 of 
    88             (size(u))[0] NE 3 OR (size(v))[0] NE 3: return,  -1 
     93            (size(v))[0] NE 3: return,  -1 
    8994            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    9095             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    9196               case 1 of 
    92                   nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
    93                   nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    94                   nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
    95                   nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     97                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
     98                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     99                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
     100                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    96101                  ELSE : 
    97102               endcase 
    … …  
    114119         zu = reform(zu, nx, ny, nzt, /over) 
    115120         zu = temporary(u)*temporary(zu) $ 
    116           *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] 
     121          *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    117122         terreu = where(zu EQ 0) 
    118123         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan 
    … …  
    121126         zv = reform(zv, nx, ny, nzt, /over) 
    122127         zv = temporary(v)*temporary(zv) $ 
    123           *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] 
     128          *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    124129         terrev = where(zv EQ 0) 
    125130         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan 
    … …  
    129134         zdiv = reform(zdiv, nx, ny, nzt, /over) 
    130135         zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0)) $ 
    131           *tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] 
     136          *tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    132137;------------------------------------------------------------ 
    133138; mise a !values.f_nan de la bordure 
    134139;------------------------------------------------------------ 
    135          if  NOT keyword_set(key_periodique) OR nx NE jpi then begin 
     140         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    136141            zdiv(0, *, *) = !values.f_nan 
    137142            zdiv(nx-1, *, *) = !values.f_nan 
    … …  
    143148; 
    144149         if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    145          terre =  where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] EQ 0) 
     150         terre =  where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0) 
    146151         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask 
    147152;------------------------------------------------------------ 
    … …  
    151156         varname = 'div' 
    152157         varunits = '1e6*s-1' 
    153          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t'] 
    154          if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan, boite = boite) 
     158         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t'] 
     159         if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan) 
    155160      END 
    156161;---------------------------------------------------------------------------- 
    … …  
    168173             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    169174               case 1 of 
    170                   nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
    171                   nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    172                   nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
    173                   nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     175                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
     176                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     177                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
     178                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    174179                  ELSE : 
    175180               endcase 
    … …  
    185190; calcul de la divergence 
    186191;------------------------------------------------------------ 
    187          zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj*premierzt] 
     192         zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 
    188193         terreu = where(zu EQ 0) 
    189194         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan 
    … …  
    192197         zu = temporary(u)*temporary(zu) 
    193198; 
    194          zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj*premierzt] 
     199         zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 
    195200         terrev = where(zv EQ 0) 
    196201         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan 
    … …  
    199204         zv = temporary(v)*temporary(zv) 
    200205; 
    201          zdiv = 1e6*tmask[indice2d+jpi*jpj*premierzt]/(e1t[indice2d]*e2t[indice2d]) 
     206         zdiv = 1e6*tmask[indice2d+jpi*jpj*firstzt]/(e1t[indice2d]*e2t[indice2d]) 
    202207         zdiv = (zdiv)[*]#replicate(1, jpt) 
    203208         zdiv = reform(zdiv, nx, ny, jpt, /over) 
    … …  
    207212; mise a !values.f_nan de la bordure 
    208213;------------------------------------------------------------ 
    209          if  NOT keyword_set(key_periodique) OR nx NE jpi then begin 
     214         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    210215            zdiv(0, *, *) = !values.f_nan 
    211216            zdiv(nx-1, *, *) = !values.f_nan 
    … …  
    222227         varname = 'div' 
    223228         varunits = '1e6*s-1' 
    224          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t'] 
    225          if keyword_set(direc) then  zdiv = grossemoyenne(zdiv,direc,/nan, boite = boite) 
     229         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t'] 
     230         if keyword_set(direc) then  zdiv = grossemoyenne(zdiv,direc,/nan) 
    226231      END 
    227232;---------------------------------------------------------------------------- 
    … …  
    240245      ELSE:BEGIN                ;xy 
    241246         indice3d = lindgen(jpi, jpj, jpk) 
    242          indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, niveau -1] 
     247         indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt] 
    243248;------------------------------------------------------------ 
    244249; extraction de u et v sur le domaine qui convient 
    … …  
    252257             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    253258               case 1 of 
    254                   nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
    255                   nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
    256                   nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
    257                   nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
     259                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
     260                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
     261                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
     262                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
    258263                  ELSE : 
    259264               endcase 
    … …  
    280285; mise a !values.f_nan de la bordure 
    281286;------------------------------------------------------------ 
    282          if  NOT keyword_set(key_periodique) OR nx NE jpi then begin 
     287         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    283288            zdiv(0, *) = !values.f_nan 
    284289            zdiv(nx-1, *) = !values.f_nan 
    … …  
    298303         varname = 'div' 
    299304         varunits = '1e6*s-1' 
    300          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t'] 
    301          if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan, boite = boite) 
     305         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t'] 
     306         if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan) 
    302307 
    303308      END 

  • trunk/ToBeReviewed/CALCULS/grad.pro

    r23 r25  
    3434@common 
    3535;------------------------------------------------------------ 
     36; 
     37   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 
     38     return, report(['This version of grad is based on Arakawa C-grid.' $ 
     39                     , 'U and V grids must therefore be defined']) 
     40; 
    3641   res = litchamp(field) 
    3742   taille=size(res) 
    38    grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz      
     43   grille, mask, glam, gphi, gdep, nx, ny,nz,firstx,firsty,firstz,lastx, lasty, lastz      
    3944   if n_elements(valmask) EQ 0 then valmask = 1e20 
    4045   case strupcase(vargrid) of 
    … …  
    4247         case direc of 
    4348            'x':BEGIN  
    44                divi = e1u[premierx:dernierx, premiery:derniery] 
    45                newmask = (umask())[premierx:dernierx, premiery:derniery, premierz:dernierz] 
     49               divi = e1u[firstx:lastx, firsty:lasty] 
     50               newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] 
    4651               vargrid = 'U' 
    47                domdef, glamt[premierx, 0], glamu[dernierx, 0] $ 
    48                 , gphit[0, premiery], gphiu[0, derniery], grille = ['T','U'] 
     52               domdef, glamt[firstx, 0], glamu[lastx, 0] $ 
     53                , gphit[0, firsty], gphiu[0, lasty], gridtype = ['T','U'] 
    4954            END 
    5055            'y':BEGIN 
    51                divi = e2v[premierx:dernierx, premiery:derniery] 
    52                newmask = (vmask())[premierx:dernierx, premiery:derniery, premierz:dernierz] 
     56               divi = e2v[firstx:lastx, firsty:lasty] 
     57               newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
    5358               vargrid = 'V' 
    54                domdef, glamt[premierx, 0], glamv[dernierx, 0] $ 
    55                 , gphit[0, premiery], gphiv[0, derniery], grille = ['T','V'] 
    56             END 
    57             'z':BEGIN 
    58                divi = e3w[premierz:dernierz] 
     59               domdef, glamt[firstx, 0], glamv[lastx, 0] $ 
     60                , gphit[0, firsty], gphiv[0, lasty], gridtype = ['T','V'] 
     61            END 
     62            'z':BEGIN 
     63               divi = e3w[firstz:lastz] 
    5964               newmask = mask 
    6065               vargrid = 'W' 
    … …  
    6570      'W':BEGIN 
    6671         case direc of 
    67             'x':divi = e1u[premierx:dernierx, premiery:derniery] 
    68             'y':divi = e2v[premierx:dernierx, premiery:derniery] 
    69             'z':BEGIN 
    70                divi = e3t[premierz:dernierz] 
     72            'x':divi = e1u[firstx:lastx, firsty:lasty] 
     73            'y':divi = e2v[firstx:lastx, firsty:lasty] 
     74            'z':BEGIN 
     75               divi = e3t[firstz:lastz] 
    7176               newmask = mask 
    7277               vargrid = 'T' 
    … …  
    7883         case direc of 
    7984            'x':BEGIN 
    80                divi = (shift(e1t, -1, 0))[premierx:dernierx, premiery:derniery] 
    81                newmask = tmask[premierx:dernierx, premiery:derniery, premierz:dernierz] 
     85               divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty] 
     86               newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 
    8287               vargrid = 'T' 
    83                domdef, glamt[premierx, 0], glamu[dernierx] $ 
    84                 , gphit[0, premiery], gphiu[0, derniery], grille = ['T','U'] 
     88               domdef, glamt[firstx, 0], glamu[lastx] $ 
     89                , gphit[0, firsty], gphiu[0, lasty], gridtype = ['T','U'] 
    8590            END 
    8691            'y':BEGIN 
    87                divi = e2f[premierx:dernierx, premiery:derniery] 
    88                newmask = (fmask())[premierx:dernierx, premiery:derniery, premierz:dernierz] 
     92               divi = e2f[firstx:lastx, firsty:lasty] 
     93               newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
    8994               vargrid = 'F' 
    90                domdef, glamu[premierx, 0], glamf[dernierx, 0] $ 
    91                 , gphiu[0, premiery], gphif[0, derniery], grille = ['U','F'] 
    92             END 
    93             'z':BEGIN 
    94                divi = e3w[premierz:dernierz] 
     95               domdef, glamu[firstx, 0], glamf[lastx, 0] $ 
     96                , gphiu[0, firsty], gphif[0, lasty], gridtype = ['U','F'] 
     97            END 
     98            'z':BEGIN 
     99               divi = e3w[firstz:lastz] 
    95100               newmask = mask 
    96101               vargrid = 'W' 
    … …  
    102107         case direc of 
    103108            'x':BEGIN 
    104                divi = e1f[premierx:dernierx, premiery:derniery] 
    105                newmask = (fmask())[premierx:dernierx, premiery:derniery, premierz:dernierz] 
     109               divi = e1f[firstx:lastx, firsty:lasty] 
     110               newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
    106111               vargrid = 'F' 
    107                domdef, glamv[premierx, 0], glamf[dernierx, 0] $ 
    108                 , gphiv[0, premiery], gphif[0, derniery], grille = ['V','F'] 
     112               domdef, glamv[firstx, 0], glamf[lastx, 0] $ 
     113                , gphiv[0, firsty], gphif[0, lasty], gridtype = ['V','F'] 
    109114            END 
    110115            'y':BEGIN 
    111                divi = (shift(e2t, 0, -1))[premierx:dernierx, premiery:derniery] 
    112                newmask = tmask[premierx:dernierx, premiery:derniery, premierz:dernierz] 
     116               divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty] 
     117               newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 
    113118               vargrid = 'T' 
    114                domdef, glamt[premierx, 0], glamv[dernierx, 0] $ 
    115                 , gphit[0, premiery], gphiv[0, derniery], grille = ['T','V'] 
    116             END 
    117             'z':BEGIN 
    118                divi = e3w[premierz:dernierz] 
     119               domdef, glamt[firstx, 0], glamv[lastx, 0] $ 
     120                , gphit[0, firsty], gphiv[0, lasty], gridtype = ['T','V'] 
     121            END 
     122            'z':BEGIN 
     123               divi = e3w[firstz:lastz] 
    119124               newmask = mask 
    120125               vargrid = 'W' 
    … …  
    125130;       'F':BEGIN 
    126131;          case direc of 
    127 ;             'x':divi = (shift(e1v, -1, 0))[premierx:dernierx, premiery:derniery] 
    128 ;             'y':divi = (shift(e2u, 0, -1))[premierx:dernierx, premiery:derniery] 
    129 ;             'z':divi = e3w[premierz:dernierz] 
     132;             'x':divi = (shift(e1v, -1, 0))[firstx:lastx, firsty:lasty] 
     133;             'y':divi = (shift(e2u, 0, -1))[firstx:lastx, firsty:lasty] 
     134;             'z':divi = e3w[firstz:lastz] 
    130135;             ELSE:return, report('Bad definition of direction argument') 
    131136;          endcase 
    … …  
    139144;---------------------------------------------------------------------------- 
    140145      taille[0] EQ 2:BEGIN 
    141          earth = where(mask[*, *, premierz] EQ 0) 
     146         earth = where(mask[*, *, firstz] EQ 0) 
    142147         if earth[0] NE -1 then res[earth] = !values.f_nan 
    143148         case direc of 
    144149            'x':BEGIN  
    145150               res = (shift(res, -1, 0)-res)/divi 
    146                if key_periodique EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan 
     151               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan 
    147152               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0) 
    148153            END 
    … …  
    154159            ELSE:return,  report('Bad definition of direction argument for the type of array') 
    155160         ENDCASE 
    156          earth = where(newmask[*, *, premierz] EQ 0) 
     161         earth = where(newmask[*, *, firstz] EQ 0) 
    157162         if earth[0] NE -1 then res[earth] = valmask 
    158163      END 
    … …  
    161166;---------------------------------------------------------------------------- 
    162167      taille[0] EQ 3 AND jpt NE 1:BEGIN 
    163          earth = where(mask[*, *, premierz] EQ 0) 
     168         earth = where(mask[*, *, firstz] EQ 0) 
    164169         if earth[0] NE -1 then BEGIN 
    165170            earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt)) 
    … …  
    170175            'x':BEGIN  
    171176               res = (shift(res, -1, 0, 0)-res)/divi 
    172                if key_periodique EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 
     177               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 
    173178               if vargrid EQ 'T' OR vargrid EQ 'V'  then res = shift(res, 1, 0, 0) 
    174179            END 
    … …  
    180185            ELSE:return,  report('Bad definition of direction argument for the type of array') 
    181186         ENDCASE 
    182          earth = where(newmask[*, *, premierz] EQ 0) 
    183          if earth[0] NE -1 then BEGIN  
    184             earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt)) 
    185             res[earth] = valmask 
    186          ENDIF 
     187         earth = where(newmask[*, *, firstz] EQ 0) 
     188         if earth[0] NE -1 THEN res[earth] = valmask 
    187189      END 
    188190;---------------------------------------------------------------------------- 
    … …  
    196198               divi = divi[*]#replicate(1, nz) 
    197199               res = (shift(res, -1, 0, 0)-res)/divi 
    198                if key_periodique EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 
     200               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 
    199201               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0) 
    200202            END 
    … …  
    227229         if earth[0] NE -1 then res[earth] = !values.f_nan 
    228230         case direc OF 
     231            'x':BEGIN  
     232               divi = divi[*]#replicate(1, nz*jpt) 
     233               res = (shift(res, -1, 0, 0, 0)-res)/divi 
     234               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *, *] = !values.f_nan 
     235               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0, 0) 
     236            END 
     237            'y':BEGIN  
     238               divi = divi[*]#replicate(1, nz*jpt) 
     239               res = (shift(res, 0, -1, 0, 0)-res)/divi 
     240               res[*, ny-1, *, *] = !values.f_nan 
     241               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0, 0)         
     242            END 
    229243            'z':BEGIN 
    230244               divi = replicate(1, nx*ny)#divi 
    … …  
    238252                  res[*, *, nz-1, *] = !values.f_nan 
    239253               ENDELSE 
    240                if earth[0] NE -1 then res[earth] = valmask 
    241             END 
    242          ENDCASE 
     254            END 
     255         ENDCASE 
     256         if earth[0] NE -1 then res[earth] = valmask 
    243257      END 
    244258;------------------------------------------------------------ 

  • trunk/ToBeReviewed/CALCULS/grossemoyenne.pro

    r23 r25  
    1111; CATEGORY: 
    1212; 
    13 ; CALLING SEQUENCE: result = grossemoyenne(tab,'direc',BOITE=boite) 
     13; CALLING SEQUENCE: result = grossemoyenne(tab,'direc',BOXZOOM=boxzoom) 
    1414; 
    1515; INPUTS:       tab = 3 or 4d field 
    … …  
    1919; KEYWORD PARAMETERS: 
    2020; 
    21 ;               boite = [xmin,xmax,ymin,ymax (,(zmin,)zmax)] pour plus 
     21;               boxzoom = [xmin,xmax,ymin,ymax (,(zmin,)zmax)] pour plus 
    2222;               de detail cf domdef 
    23 ;               boite peut prendre 5 formes:  
    24 ;               [prof2], [prof1, prof2],[lon1, lon2, lat1, lat2],   
    25 ;               [lon1, lon2, lat1, lat2, prof2],[lon1, lon2, lat1, lat2, prof1,prof2] 
     23;               boxzoom peut prendre 5 formes:  
     24;               [vert2], [vert1, vert2],[lon1, lon2, lat1, lat2],   
     25;               [lon1, lon2, lat1, lat2, vert2],[lon1, lon2, lat1, 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 grossemoyenne est appelee via 
    4242;               checkfield) 
    4343; 
    44 ;               INTEGRATION: pour faire une integrale plutot qu''une moyenne 
     44;               INTEGRATION: pour faire une integrale plutot qu''une 
     45;               moyenne 
     46; 
     47;               /SPATIALFIRST, when performing at the same time 
     48;               spatial and temporal mean, grossemoyenne is assuming 
     49;               that the mask is not changing with the time. In 
     50;               consequence, grossemoyenne performs temporal mean 
     51;               first and then call moyenne. Activate /SPATIALFIRST if 
     52;               you want to perform the spatial mean before the 
     53;               temporal mean. Note that if NAN is activated, then 
     54;               SPATIALFIRST is activated automatically. 
     55; 
     56;               /TEMPORALFIRST: to force to perform first temporal 
     57;               mean even if nan is activated (see SPATIALFIRST 
     58;               explanations...) 
     59; 
     60;  
     61;         /WDEPTH: to specify that the field is at W depth instad of T  
     62;         depth (automatically activated if vargrid eq 'W') 
    4563;             
    4664; OUTPUTS:  
    … …  
    90108;------------------------------------------------------------ 
    91109;------------------------------------------------------------ 
    92 function grossemoyenne ,tab,direc,BOITE=boite, INTEGRATION=integration $ 
    93                         , NAN = nan, NODOMDEF = nodomdef, _extra = ex 
    94 @common 
    95    tempsun = systime(1)         ; pour key_performance 
     110function grossemoyenne, tab, direc, BOXZOOM = boxzoom, INTEGRATION = integration $ 
     111                        , NAN = nan, NODOMDEF = nodomdef, WDEPTH = wdepth $ 
     112                        , SPATIALFIRST = spatialfirst, TEMPORALFIRST = temporalfirst $ 
     113                        , _extra = ex 
     114;--------------------------------------------------------- 
     115@cm_4mesh 
     116@cm_4data 
     117@cm_4cal 
     118  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     119@updatenew 
     120@updatekwd 
     121  ENDIF 
     122;--------------------------------------------------------- 
     123  tempsun = systime(1)          ; pour key_performance 
    96124;------------------------------------------------------------ 
    97125;I) preliminaires 
    98126;------------------------------------------------------------ 
    99    dirt = 0 
    100    dirx = 0 
    101    diry = 0 
    102    dirz = 0 
    103    dim  = 'aa' 
     127  dirt = 0 
     128  dirx = 0 
     129  diry = 0 
     130  dirz = 0 
     131  dim  = 'aa' 
    104132;------------------------------------------------------------ 
    105133; I.1) direction(s) suivants lesquelles on integre 
    106134;------------------------------------------------------------ 
    107    if ( strpos(direc,'t') ge 0 ) then dirt = 1 
    108    if ( strpos(direc,'x') ge 0 ) then dirx = 1 
    109    if ( strpos(direc,'y') ge 0 ) then diry = 1 
    110    if ( strpos(direc,'z') ge 0 ) then dirz = 1 
     135  if ( strpos(direc, 't') ge 0 ) then dirt = 1 
     136  if ( strpos(direc, 'x') ge 0 ) then dirx = 1 
     137  if ( strpos(direc, 'y') ge 0 ) then diry = 1 
     138  if ( strpos(direc, 'z') ge 0 ) then dirz = 1 
     139  IF keyword_set(NAN) AND (dirx EQ 1 OR diry EQ 1 OR dirz EQ 1) $ 
     140    THEN spatialfirst = 1 
     141  IF keyword_set(temporalfirst) THEN spatialfirst = 0 
    111142;------------------------------------------------------------ 
    112143; I.2) verification de la taille du tableau d'entree 
    113144;------------------------------------------------------------ 
    114    taille = size(tab) 
    115    case 1 of 
    116       taille[0] eq 1 :return,  report('Le tableau n''a qu''une dimension, cas non traite!') 
    117       taille[0] eq 2 :return,  report('Le tableau n''a qu''deux dimension, cas non traite!') 
    118       taille[0] eq 3 :BEGIN  
    119          dim = '3d' 
    120          if dirx eq 0 and diry eq 0 and dirt eq 0 then return, tab 
    121       END 
    122       taille[0] eq 4 :BEGIN  
    123          dim = '4d' 
    124          if dirx eq 0 and diry eq 0 and dirz eq 0 and dirt eq 0 then return, tab 
    125       END 
    126       else : return, report('Le tableau d entree doit etre a 3 ou 4 dimensions s''il ne contient pas de dim temporelle utiliser moyenne') 
    127    endcase 
     145  taille = size(tab) 
     146  case 1 of 
     147    taille[0] eq 1 :return,  report('Le tableau n''a qu''une dimension, cas non traite!') 
     148    taille[0] eq 2 :return,  report('Le tableau n''a qu''deux dimension, cas non traite!') 
     149    taille[0] eq 3 :BEGIN  
     150      dim = '3d' 
     151      if dirx eq 0 and diry eq 0 and dirt eq 0 then return, tab 
     152    END 
     153    taille[0] eq 4 :BEGIN  
     154      dim = '4d' 
     155      if dirx eq 0 and diry eq 0 and dirz eq 0 and dirt eq 0 then return, tab 
     156    END 
     157    else : return, report('Le tableau d entree doit etre a 3 ou 4 dimensions s''il ne contient pas de dim temporelle utiliser moyenne') 
     158  endcase 
    128159;------------------------------------------------------------ 
    129160;   I.4) obtention des facteurs d''echelles et du masque sur le sous 
    130161;   domaine concerne par la moyenne  
    131 ; redefinition du domaine ajuste a boite (a 6 elements) 
     162; redefinition du domaine ajuste a boxzoom (a 6 elements) 
    132163; ceci va nous permetre de faire les calcules que sur le sous domaine 
    133164; comcerne par la moyenne. domdef, suivit de grille nous donne tous 
    134165; les tableaux de la grille sur le sous domaine  
    135166;------------------------------------------------------------ 
    136    if keyword_set(boite) then BEGIN  
    137       Case 1 Of 
    138          N_Elements(Boite) Eq 1: bte = [lon1, lon2, lat1, lat2, 0.,boite[0]] 
    139          N_Elements(Boite) Eq 2: bte = [lon1, lon2, lat1, lat2, boite[0],boite[1]] 
    140          N_Elements(Boite) Eq 4: bte = [Boite,prof1, prof2] 
    141          N_Elements(Boite) Eq 5: bte = [Boite[0:3], 0, Boite[4]] 
    142          N_Elements(Boite) Eq 6: bte = Boite 
    143          Else: return, report('Mauvaise Definition de Boite') 
    144       endcase 
    145       oldboite = [lon1, lon2, lat1, lat2, prof1, prof2] 
    146       if NOT keyword_set(nodomdef) then domdef, bte,GRILLE=vargrid, _extra = ex 
    147    ENDIF  
     167  if keyword_set(boxzoom) then BEGIN  
     168    Case 1 Of 
     169      N_Elements(Boxzoom) Eq 1: bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]] 
     170      N_Elements(Boxzoom) Eq 2: bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]] 
     171      N_Elements(Boxzoom) Eq 4: bte = [Boxzoom, vert1, vert2] 
     172      N_Elements(Boxzoom) Eq 5: bte = [Boxzoom[0:3], 0, Boxzoom[4]] 
     173      N_Elements(Boxzoom) Eq 6: bte = Boxzoom 
     174      Else: return, report('Wrong Definition of Boxzoom') 
     175    endcase 
     176    if NOT keyword_set(nodomdef) then BEGIN  
     177      savedbox = 1b 
     178      saveboxparam, 'boxparam4grmoyenne.dat' 
     179      domdef, bte, GRIDTYPE = vargrid, _extra = ex 
     180    ENDIF  
     181  ENDIF  
    148182;--------------------------------------------------------------- 
    149183; attribution du mask et des tableaux de longitude et latitude... 
    150184;--------------------------------------------------------------- 
    151    grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz, dernierx, derniery, dernierz,e1,e2,e3 
     185  grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, WDEPTH = wdepth 
    152186;------------------------------------------------------------ 
    153187; I.3) si dirt eq 1 on fait la moyenne temporelle et on envoie ds moyenne 
    154188;------------------------------------------------------------ 
    155    if dirt EQ 1 then begin 
    156       if dim EQ 3d then BEGIN  
     189  if dirt EQ 1 AND NOT keyword_set(spatialfirst) then begin 
     190    if dim EQ 3d then BEGIN  
    157191      case 1 of 
    158          taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $ 
    159           res = tab[premierx:premierx+nx-1 $ 
    160                     ,premiery:premiery+ny-1, *] 
    161          taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpt:res=tab 
    162          else: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)) 
     192        taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $ 
     193          res = tab[firstx:firstx+nx-1 $ 
     194                    , firsty:firsty+ny-1, *] 
     195        taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpt:res = tab 
     196        else:BEGIN  
     197          if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
     198          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)) 
     199        END 
    163200      ENDCASE 
    164          divi = finite(res) 
    165          divi = total(temporary(divi),3, nan = keyword_set(nan)) 
    166          notanum = where(divi EQ 0) 
    167          if keyword_set(integration) then begin 
    168             res = total(res,3, nan = nan) 
    169          ENDIF ELSE BEGIN 
    170             if keyword_set(nan) then BEGIN 
    171                res = total(res,3, nan = keyword_set(nan))/ (1 > divi) 
    172             ENDIF ELSE res = total(res,3)/(1.*taille[3])  
    173          ENDELSE 
     201      if keyword_set(integration) then begin 
     202        res = total(res, 3, nan = nan) 
    174203      ENDIF ELSE BEGIN 
     204        if keyword_set(nan) then BEGIN 
     205          divi = finite(res) 
     206          divi = total(temporary(divi), 3) 
     207          notanum = where(divi EQ 0) 
     208          res = total(res, 3, nan = keyword_set(nan))/ (1 > divi) 
     209          if notanum[0] NE -1 then res[temporary(notanum)] = !values.f_nan 
     210        ENDIF ELSE res = total(res, 3)/(1.*taille[3])  
     211      ENDELSE 
     212    ENDIF ELSE BEGIN 
    175213      case 1 of 
    176          taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $ 
    177           res =tab[premierx:dernierx,premiery:derniery,premierz:dernierz, *] 
    178          taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $ 
    179           res =tab[premierx:dernierx,premiery:derniery,*, *] 
    180          taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq nz and taille[4] eq jpt:res=tab 
    181          taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk and taille[4] eq jpt: $ 
    182           res=tab[*, *, premierz:dernierz, *] 
    183          else:return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $ 
    184                              +strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+'*' $ 
    185                              +strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*' $ 
    186                              +strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)+'*' $ 
    187                              +strtrim(taille[4], 1)) 
    188       endcase 
    189          divi = finite(res) 
    190          divi = total(temporary(divi),4, nan = keyword_set(nan)) 
    191          notanum = where(divi EQ 0) 
    192          if keyword_set(integration) then begin 
    193             res = total(res,4, nan = nan) 
    194          ENDIF ELSE BEGIN 
    195             if keyword_set(nan) then begin 
    196                res = total(res,4, /nan)/(1 > divi) 
    197             ENDIF ELSE res = total(res,4)/(1.*taille[4]) 
    198          ENDELSE 
     214        taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $ 
     215          res = tab[firstx:lastx, firsty:lasty, firstz:lastz, *] 
     216        taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $ 
     217          res = tab[firstx:lastx, firsty:lasty, *, *] 
     218        taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq nz and taille[4] eq jpt:res = tab 
     219        taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk and taille[4] eq jpt: $ 
     220          res = tab[*, *, firstz:lastz, *] 
     221        else:BEGIN  
     222           if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
     223          return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $ 
     224                         +strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+'*' $ 
     225                         +strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*' $ 
     226                         +strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)+'*' $ 
     227                         +strtrim(taille[4], 1)) 
     228        END 
     229        endcase 
     230      if keyword_set(integration) then begin 
     231        res = total(res, 4, nan = nan) 
     232      ENDIF ELSE BEGIN 
     233        if keyword_set(nan) then begin 
     234          divi = finite(res) 
     235          divi = total(temporary(divi), 4) 
     236          notanum = where(divi EQ 0) 
     237          res = total(res, 4, /nan)/(1 > divi) 
     238          if notanum[0] NE -1 then res[temporary(notanum)] = !values.f_nan 
     239        ENDIF ELSE res = total(res, 4)/(1.*taille[4]) 
    199240      ENDELSE 
    200       if notanum[0] NE -1 then res[temporary(notanum)] = !values.f_nan 
    201       return,  moyenne(temporary(res), direc, BOITE = boite, NAN = nan, INTEGRATION = integration, _extra = ex) 
    202    endif 
     241    ENDELSE 
     242    if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'     
     243    return,  moyenne(temporary(res), direc, BOXZOOM = boxzoom, NAN = nan, INTEGRATION = integration, NODOMDEF = nodomdef, WDEPTH = wdepth, _extra = ex) 
     244  ENDIF ELSE res = tab 
     245  IF jpt EQ 1 THEN BEGIN 
     246    if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'     
     247    return, moyenne(reform(res, /over), direc, BOXZOOM = boxzoom, NAN = nan, INTEGRATION = integration, NODOMDEF = nodomdef, WDEPTH = wdepth, _extra = ex) 
     248  END 
    203249;------------------------------------------------------------ 
    204250;------------------------------------------------------------ 
    … …  
    206252;------------------------------------------------------------ 
    207253;------------------------------------------------------------ 
    208    if (dim eq '3d') then begin 
     254  if (dim eq '3d') then begin 
    209255;--------------------------------------------------------------- 
    210256;   II.1) verification de la coherence de la taille du tableau a 
    … …  
    215261; (jpi,jpj,jpt) soit celle du domaine reduit (nx,ny,jpt) 
    216262;--------------------------------------------------------------- 
    217       case 1 of 
    218          taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $ 
    219           res = tab[premierx:premierx+nx-1 $ 
    220                     ,premiery:premiery+ny-1, *] 
    221          taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpt:res=tab 
    222          else: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)) 
    223       endcase 
    224       if keyword_set(nan) NE 0 then BEGIN  
    225          if nan NE 1 then BEGIN ; si nan n''est pas !values.f_nan  
     263    case 1 of 
     264      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $ 
     265        res = tab[firstx:firstx+nx-1 $ 
     266                  , firsty:firsty+ny-1, *] 
     267      taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpt:res = tab 
     268      else:BEGIN  
     269        if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
     270        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)) 
     271      enD 
     272    endcase 
     273    if keyword_set(nan) NE 0 then BEGIN  
     274      if nan NE 1 then BEGIN    ; si nan n''est pas !values.f_nan  
    226275; on le met a !values.f_nan 
    227             if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 
    228             ELSE notanumber = where(abs(res) GT abs(nan)/10.) 
    229             if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 
    230          ENDIF 
     276        if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 
     277        ELSE notanumber = where(abs(res) GT abs(nan)/10.) 
     278        if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 
    231279      ENDIF 
     280    ENDIF 
    232281;--------------------------------------------------------------- 
    233282; rq IL FAUT FAIRE ATTENTION AUX CAS OU LA DIM A MOYENNER = 1, ET 
    … …  
    235284; PEUVENT SEMBLER INUTILE AU DEPART 
    236285;--------------------------------------------------------------- 
    237       if nx EQ 1 OR ny EQ 1 OR jpt EQ 1 then res=reform(res,nx,ny,jpt, /over) 
    238       if nx EQ 1 OR ny EQ 1 then BEGIN  
    239          mask =  reform(mask, nx, ny, nz, /over) 
    240          e1 =  reform(e1, nx, ny, /over) 
    241          e2 =  reform(e2, nx, ny, /over) 
    242       endif 
     286    if nx EQ 1 OR ny EQ 1 then BEGIN  
     287      res = reform(res, nx, ny, jpt, /over) 
     288      e1 =  reform(e1, nx, ny, /over) 
     289      e2 =  reform(e2, nx, ny, /over) 
     290    endif 
     291    if nx EQ 1 OR ny EQ 1 OR nz EQ 1 THEN $ 
     292      mask =  reform(mask, nx, ny, nz, /over) 
    243293;--------------------------------------------------------------- 
    244294; II.3) differents types de moyennes 
    245295;--------------------------------------------------------------- 
    246       if keyword_set(nan) NE 0 then begin 
    247          msknan = replicate(1., nx, ny, jpt) 
    248          notanumber = where(finite(res) EQ 0) 
    249          if notanumber[0] NE -1 then msknan[temporary(notanumber)] = !values.f_nan 
    250       ENDIF ELSE msknan = 1 
    251       if nz eq jpk then mask = mask[*, *,niveau-1] ELSE mask = mask[*, *,nz-1] 
    252       case 1 of 
    253          (dirx eq 1) and (diry eq 0) : begin 
    254             e = temporary(e1)*temporary(mask) 
    255             if nx EQ 1 OR ny EQ 1 then e = reform(e, nx, ny, /over) 
    256             echelle = temporary(e[*])#replicate(1, jpt) 
    257             echelle = reform(echelle, nx, ny, jpt, /over) 
    258             if keyword_set(integration) then divi=1 $ 
    259             ELSE divi = total(echelle*msknan,1, nan = nan) 
    260             res=total(temporary(res)*echelle,1, nan = nan)/(divi > 1.) 
    261             if keyword_set(nan) then begin 
    262                testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 
    263                testnan = total(testnan,1) 
    264             endif 
    265          end 
    266          (dirx eq 0) and (diry eq 1) : begin 
    267             e = temporary(e2)*temporary(mask) 
    268             if nx EQ 1 OR ny EQ 1 then e = reform(e, nx, ny, /over) 
    269             echelle = temporary(e[*])#replicate(1, jpt) 
    270             echelle = reform(echelle, nx, ny, jpt, /over) 
    271             if keyword_set(integration) then divi=1 $ 
    272             ELSE divi = total(echelle*msknan,2, nan = nan) 
    273             res=total(temporary(res)*echelle,2, nan = nan)/(divi > 1.) 
    274             if keyword_set(nan) then begin 
    275                testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 
    276                testnan = total(testnan,2) 
    277             endif 
    278          end 
    279          (dirx eq 1) and (diry eq 1) : begin 
    280             echelle=(temporary(e1)*temporary(e2)*temporary(mask))[*]#replicate(1,jpt) 
    281             echelle=reform(echelle,nx,ny,jpt, /over) 
    282             if keyword_set(integration) then divi=1 $ 
    283             ELSE divi = total(temporary(total(echelle*msknan,1, nan = nan)),1, nan = nan) 
    284             res = total(temporary(total(temporary(res)*echelle,1, nan=nan)),1, nan=nan)/(divi > 1.) 
    285             if keyword_set(nan) then begin 
    286                testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 
    287                testnan = total(total(testnan,1),1) 
    288             endif 
    289          end 
    290       endcase 
    291    endif 
     296    if keyword_set(nan) NE 0 then msknan = finite(res) ELSE msknan = -1 
     297    mask = mask[*, *, 0] 
     298    case 1 of 
     299      (dirx eq 1) and (diry eq 0) : begin 
     300        e = temporary(e1)*temporary(mask) 
     301        echelle = (temporary(e))[*]#replicate(1, jpt) 
     302        echelle = reform(echelle, nx, ny, jpt, /over) 
     303        if keyword_set(integration) then divi = 1 ELSE BEGIN  
     304          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 1) $ 
     305          ELSE divi = total(echelle, 1) 
     306        ENDELSE 
     307        res = total(temporary(res)*echelle, 1, nan = nan)/(divi > 1.) 
     308        if msknan[0] NE -1 then BEGIN 
     309          echelle = temporary(echelle) NE 0 
     310          testnan = temporary(msknan)*echelle 
     311          testnan = total(temporary(testnan), 1) $ 
     312            +(total(temporary(echelle), 1) EQ 0) 
     313        endif 
     314      end 
     315      (dirx eq 0) and (diry eq 1) : begin 
     316        e = temporary(e2)*temporary(mask) 
     317        if nx EQ 1 OR ny EQ 1 then e = reform(e, nx, ny, /over) 
     318        echelle = (temporary(e))[*]#replicate(1, jpt) 
     319        echelle = reform(echelle, nx, ny, jpt, /over) 
     320        if keyword_set(integration) then divi = 1 ELSE BEGIN  
     321          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 2) $ 
     322          ELSE divi = total(echelle, 2) 
     323        ENDELSE  
     324        res = total(temporary(res)*echelle, 2, nan = nan)/(divi > 1.) 
     325        if msknan[0] NE -1 then begin 
     326          echelle = temporary(echelle) NE 0 
     327          testnan = temporary(msknan)*echelle 
     328          testnan = total(temporary(testnan), 2) $ 
     329            +(total(temporary(echelle), 2) EQ 0) 
     330        endif 
     331      end 
     332      (dirx eq 1) and (diry eq 1) : begin 
     333        echelle = (temporary(e1)*temporary(e2)*temporary(mask))[*]#replicate(1, jpt) 
     334        echelle = reform(echelle, nx, ny, jpt, /over) 
     335        if keyword_set(integration) then divi = 1 ELSE BEGIN  
     336          IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 1) $ 
     337          ELSE divi = total(total(echelle, 1), 1) 
     338        ENDELSE 
     339        res = total(temporary(total(temporary(res)*echelle, 1, nan = nan)), 1, nan = nan)/(divi > 1.) 
     340        if msknan[0] NE -1 then begin 
     341          echelle = temporary(echelle) NE 0 
     342          testnan = temporary(msknan)*echelle 
     343          testnan = total(total(temporary(testnan), 1), 1) $ 
     344            +(total(total(temporary(echelle), 1), 1) EQ 0) 
     345        endif 
     346      end 
     347    endcase 
     348  endif 
    292349;------------------------------------------------------------ 
    293350;------------------------------------------------------------ 
    294351; III) Cas serie tableaux 3d (tab4d) 
    295352;------------------------------------------------------------ 
    296    if (dim eq '4d') then begin 
     353  if (dim eq '4d') then begin 
    297354;--------------------------------------------------------------- 
    298355; III.1) verification de la coherence de la taille du tableau a 
    … …  
    303360; (jpi,jpj,jpk,jpt) soit celle du domaine reduit (nx,ny,ny,jpt) 
    304361;--------------------------------------------------------------- 
    305       case 1 of 
    306          taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $ 
    307           res =tab[premierx:dernierx,premiery:derniery,premierz:dernierz, *] 
    308          taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $ 
    309           res =tab[premierx:dernierx,premiery:derniery,*, *] 
    310          taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq nz and taille[4] eq jpt:res=tab 
    311          taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk and taille[4] eq jpt: $ 
    312           res=tab[*, *, premierz:dernierz, *] 
    313          else:return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $ 
    314                              +strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+'*' $ 
    315                              +strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*' $ 
    316                              +strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)+'*' $ 
    317                              +strtrim(taille[4], 1)) 
    318       endcase 
    319       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) 
    320       if keyword_set(nan) NE 0 then BEGIN  
    321          if nan NE 1 then BEGIN ; si nan n''est pas !values.f_nan  
     362    case 1 of 
     363      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $ 
     364        res = tab[firstx:lastx, firsty:lasty, firstz:lastz, *] 
     365      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $ 
     366        res = tab[firstx:lastx, firsty:lasty, *, *] 
     367      taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq nz and taille[4] eq jpt:res = tab 
     368      taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk and taille[4] eq jpt: $ 
     369        res = tab[*, *, firstz:lastz, *] 
     370      else:BEGIN  
     371        if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
     372        return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $ 
     373                       +strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+'*' $ 
     374                       +strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*' $ 
     375                       +strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)+'*' $ 
     376                       +strtrim(taille[4], 1)) 
     377      END 
     378    endcase 
     379    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) 
     380    if keyword_set(nan) NE 0 then BEGIN  
     381      if nan NE 1 then BEGIN    ; si nan n''est pas !values.f_nan  
    322382; on le met a !values.f_nan 
    323             if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 
    324             ELSE notanumber = where(abs(res) GT abs(nan)/10.) 
    325             if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 
    326          ENDIF 
     383        if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 
     384        ELSE notanumber = where(abs(res) GT abs(nan)/10.) 
     385        if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 
    327386      ENDIF 
     387    ENDIF 
    328388;--------------------------------------------------------------- 
    329389; rq IL FAUT FAIRE ATTENTION AUX CAS OU LA DIM A MOYENNER = 1, ET 
    … …  
    331391; PEUVENT SEMBLER INUTILE AU DEPART 
    332392;--------------------------------------------------------------- 
    333       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) 
    334       if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then mask =  reform(mask, nx, ny, nz, /over) 
     393    if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then BEGIN  
     394      res = reform(res, nx, ny, nz, jpt, /over) 
     395      mask =  reform(mask, nx, ny, nz, /over) 
     396    ENDIF 
     397    IF keyword_set(key_partialstep) THEN BEGIN 
     398; the top of the ocean floor is 
     399      IF vargrid EQ 'T' OR vargrid EQ 'W' THEN bottom = total(mask, 3) $ 
     400      ELSE bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)  
     401; we suppress columns with only ocean or land  
     402      good = where(bottom NE 0 AND bottom NE nz) 
     403; the bottom of the ocean in 3D index is: 
     404      bottom = lindgen(nx*ny)+(temporary(bottom)-1L)*nx*ny 
     405      IF good[0] NE -1 THEN bottom = bottom[good] $ 
     406      ELSE bottom = -1 
     407    ENDIF ELSE bottom = -1 
    335408;--------------------------------------------------------------- 
    336409; III.2) differents types de moyennes 
    337410;--------------------------------------------------------------- 
    338       IF keyword_set(nan) NE 0 then begin 
    339          msknan = replicate(1., nx, ny, nz, jpt) 
    340          notanumber = where(finite(res) EQ 0) 
    341          if notanumber[0] NE -1 then msknan[temporary(notanumber)] = !values.f_nan 
    342       ENDIF ELSE msknan = 1 
    343       case 1 of 
    344          (dirx eq 1) and (diry eq 0) and (dirz eq 0) : begin 
    345             e13 = temporary(e1[*])#replicate(1.,nz) 
    346             e13 = reform(e13,nx,ny,nz, /over) 
    347             echelle = (temporary(e13)*temporary(mask))[*]#replicate(1, jpt) 
    348             echelle = reform(echelle, nx, ny, nz, jpt, /over) 
    349             if keyword_set(integration) then divi=1 $ 
    350             ELSE begin  
    351                if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 
    352                divi = total(temporary(divi),1, nan = nan) 
    353             endelse 
    354             res = temporary(res)*echelle 
    355             res = total(temporary(res),1, nan = nan)/(divi > 1) 
    356             if keyword_set(nan) then begin 
    357                testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 
    358                testnan = total(testnan,1) 
    359             endif 
    360          end 
    361          (dirx eq 0) and (diry eq 1) and (dirz eq 0) : begin 
    362             e23 = temporary(e2[*])#replicate(1.,nz) 
    363             e23 = reform(e23,nx,ny,nz, /over) 
    364             echelle = (temporary(e23)*temporary(mask))[*]#replicate(1, jpt) 
    365             echelle = reform(echelle, nx, ny, nz, jpt, /over) 
    366             if keyword_set(integration) then divi=1 $ 
    367             ELSE begin  
    368                if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 
    369                divi = total(temporary(divi), 2, nan = nan) 
    370             endelse 
    371             res = total(temporary(res)*echelle,2, nan = nan)/(divi > 1) 
    372             if keyword_set(nan) then begin 
    373                testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 
    374                testnan = total(testnan,2) 
    375             endif 
    376          end 
    377          (dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin 
    378             e33 = replicate(1,1.*nx*ny)#e3 
    379             e33 = reform(e33,nx,ny,nz, /over) 
    380             echelle =(temporary(e33)*temporary(mask))[*]#replicate(1, jpt) 
    381             echelle = reform(echelle, nx, ny, nz, jpt, /over) 
    382             if keyword_set(integration) then divi=1 $ 
    383             ELSE begin  
    384                if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 
    385                divi = total(temporary(divi), 3, nan = nan) 
    386             endelse 
    387             res = total(temporary(res)*echelle,3, nan = nan)/(divi > 1) 
    388             if keyword_set(nan) then begin 
    389                testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 
    390                testnan = total(testnan,3) 
    391             endif 
    392          end 
    393          (dirx eq 1) and (diry eq 1) and (dirz eq 0) : begin 
    394             e13 = e1[*]#replicate(1.,nz) 
    395             e13 = reform(e13,nx,ny,nz, /over) 
    396             e23 = e2[*]#replicate(1.,nz) 
    397             e23 = reform(e23,nx,ny,nz, /over) 
    398             echelle = (temporary(e13)*temporary(e23)*temporary(mask))[*]#replicate(1, jpt) 
    399             echelle = reform(echelle,nx,ny,nz,jpt, /over) 
    400             if keyword_set(integration) then divi=1 $ 
    401             ELSE begin  
    402                if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 
    403                divi = total(total(temporary(divi),1, nan = nan),1, nan = nan) 
    404             endelse 
    405             res =total(total(temporary(res)*echelle,1, nan = nan),1, nan = nan)/(divi > 1) 
    406             if keyword_set(nan) then begin 
    407                testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 
    408                testnan = total(total(testnan,1),1) 
    409             endif 
    410          end 
    411          (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin 
    412             e133 = e1[*]#e3 
    413             echelle = (temporary(e133[*])*temporary(mask[*]))#replicate(1, jpt) 
    414             echelle = reform(echelle,nx,ny,nz,jpt, /over) 
    415             if keyword_set(integration) then divi=1 $ 
    416             ELSE begin  
    417                if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 
    418                divi = total(total(temporary(divi),1, nan = nan),2, nan = nan) 
    419             endelse 
    420             res = total(total(temporary(res)*echelle,1, nan = nan),2, nan = nan)/(divi > 1) 
    421             if keyword_set(nan) then begin 
    422                testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 
    423                testnan = total(total(testnan,1),2) 
    424             endif 
    425          end 
    426          (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin 
    427             e233 = e2[*]#e3 
    428             echelle = (temporary(e233[*])*temporary(mask[*]))#replicate(1, jpt) 
    429             echelle = reform(echelle,nx,ny,nz,jpt, /over) 
    430             if keyword_set(integration) then divi=1 $ 
    431             ELSE begin  
    432                if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 
    433                divi = total(total(temporary(divi),2, nan = nan),2,nan=nan) 
    434             endelse 
    435             res =total(total(temporary(res)*echelle,2, nan = nan),2, nan = nan)/(divi > 1) 
    436             if keyword_set(nan) then begin 
    437                testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 
    438                testnan = total(total(testnan,2),2) 
    439             endif 
    440          end 
    441          (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin 
    442             e1233 = (e1[*]*e2[*])#e3 
    443             echelle=(temporary(e1233[*])*temporary(mask[*]))#replicate(1,jpt) 
    444             echelle=reform(echelle,nx,ny,nz,jpt, /over) 
    445             if keyword_set(integration) then divi=1 $ 
    446             ELSE begin  
    447                if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 
    448                divi = total(total(total(temporary(divi), 1,nan=nan), 1,nan=nan), 1,nan=nan) 
    449             endelse 
    450             res =total(total(total(temporary(res)*echelle, 1,nan=nan), 1,nan=nan), 1,nan=nan)/(divi > 1) 
    451             if keyword_set(nan) then begin 
    452                testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 
    453                testnan = total(total(total(testnan,1),1),1) 
    454             endif 
    455          end 
    456       endcase 
    457    endif 
     411    IF keyword_set(nan) NE 0 THEN msknan = finite(res) ELSE msknan = -1 
     412    case 1 of 
     413      (dirx eq 1) and (diry eq 0) and (dirz eq 0) : BEGIN 
     414        e13 = (temporary(e1))[*]#replicate(1., nz) 
     415        e13 = reform(e13, nx, ny, nz, /over) 
     416        echelle = (temporary(e13)*temporary(mask))[*]#replicate(1, jpt) 
     417        echelle = reform(echelle, nx, ny, nz, jpt, /over) 
     418        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 
     419          AND nx NE 1 THEN BEGIN 
     420          IF msknan[0] EQ -1 THEN BEGIN  
     421            msknan = replicate(1b, nx, ny, nz, jpt) 
     422            nan = 1 
     423          ENDIF 
     424          bottom = bottom#replicate(1, jpt) $ ; 4D bottom! 
     425            + replicate(1, n_elements(bottom))#(nx*ny*nz*lindgen(jpt)) 
     426          msknan[bottom] = 0 
     427          res[temporary(bottom)] = !values.f_nan 
     428        ENDIF 
     429        if keyword_set(integration) then divi = 1 ELSE begin  
     430          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 1) $ 
     431          ELSE divi = total(echelle, 1) 
     432        endelse 
     433        res = temporary(res)*echelle 
     434        res = total(temporary(res), 1, nan = nan)/(divi > 1) 
     435        if msknan[0] NE -1 then begin 
     436          echelle = temporary(echelle) NE 0 
     437          testnan = temporary(msknan)*echelle 
     438          testnan = total(temporary(testnan), 1) $ 
     439            +(total(temporary(echelle), 1) EQ 0) 
     440        endif 
     441      end 
     442      (dirx eq 0) and (diry eq 1) and (dirz eq 0) : begin 
     443        e23 = temporary(e2[*])#replicate(1., nz) 
     444        e23 = reform(e23, nx, ny, nz, /over) 
     445        echelle = (temporary(e23)*temporary(mask))[*]#replicate(1, jpt) 
     446        echelle = reform(echelle, nx, ny, nz, jpt, /over) 
     447        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 
     448          AND ny NE 1 THEN BEGIN 
     449          IF msknan[0] EQ -1 THEN BEGIN  
     450            msknan = replicate(1b, nx, ny, nz) 
     451            nan = 1 
     452          endif 
     453          bottom = bottom#replicate(1, jpt) $ ; 4D bottom! 
     454            + replicate(1, n_elements(bottom))#(nx*ny*nz*lindgen(jpt)) 
     455          msknan[bottom] = 0 
     456          res[temporary(bottom)] = !values.f_nan 
     457        ENDIF 
     458        if keyword_set(integration) then divi = 1 ELSE begin  
     459          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 2) $ 
     460          ELSE divi = total(echelle, 2) 
     461        endelse 
     462        res = total(temporary(res)*echelle, 2, nan = nan)/(divi > 1) 
     463        if msknan[0] NE -1 then begin 
     464          echelle = temporary(echelle) NE 0 
     465          testnan = temporary(msknan)*echelle 
     466          testnan = total(temporary(testnan), 2) $ 
     467            +(total(temporary(echelle), 2) EQ 0) 
     468        endif 
     469      end 
     470      (dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin 
     471        e33 = replicate(1, 1.*nx*ny)#e3 
     472        e33 = reform(e33, nx, ny, nz, /over) 
     473        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
     474          IF keyword_set(wdepth) THEN $ 
     475            e33[bottom] = (e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
     476          ELSE e33[bottom] = (e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
     477        ENDIF 
     478        echelle = (temporary(e33)*temporary(mask))[*]#replicate(1, jpt) 
     479        echelle = reform(echelle, nx, ny, nz, jpt, /over) 
     480        if keyword_set(integration) then divi = 1 ELSE begin  
     481          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 3) $ 
     482          ELSE divi = total(echelle, 3) 
     483        endelse 
     484        res = total(temporary(res)*echelle, 3, nan = nan)/(divi > 1) 
     485        if msknan[0] NE -1 then begin 
     486          echelle = temporary(echelle) NE 0 
     487          testnan = temporary(msknan)*echelle 
     488          testnan = total(temporary(testnan), 3) $ 
     489            +(total(temporary(echelle), 3) EQ 0) 
     490        endif 
     491      end 
     492      (dirx eq 1) and (diry eq 1) and (dirz eq 0) : begin 
     493        e13 = e1[*]#replicate(1., nz) 
     494        e13 = reform(e13, nx, ny, nz, /over) 
     495        e23 = e2[*]#replicate(1., nz) 
     496        e23 = reform(e23, nx, ny, nz, /over) 
     497        echelle = (temporary(e13)*temporary(e23)*temporary(mask))[*]#replicate(1, jpt) 
     498        echelle = reform(echelle, nx, ny, nz, jpt, /over) 
     499        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 
     500          AND nx*ny NE 1 THEN BEGIN 
     501          IF msknan[0] EQ -1 THEN BEGIN  
     502            msknan = replicate(1b, nx, ny, nz) 
     503            nan = 1 
     504          endif 
     505          bottom = bottom#replicate(1, jpt) $ ; 4D bottom! 
     506            + replicate(1, n_elements(bottom))#(nx*ny*nz*lindgen(jpt)) 
     507          msknan[bottom] = 0 
     508          res[temporary(bottom)] = !values.f_nan 
     509        ENDIF 
     510        if keyword_set(integration) then divi = 1 ELSE begin  
     511          IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 1) $ 
     512          ELSE divi = total(total(echelle, 1), 1) 
     513        endelse 
     514        res = total(total(temporary(res)*echelle, 1, nan = nan), 1, nan = nan)/(divi > 1) 
     515        if msknan[0] NE -1 then begin 
     516          echelle = temporary(echelle) NE 0 
     517          testnan = temporary(msknan)*echelle 
     518          testnan = total(total(temporary(testnan), 1), 1) $ 
     519            +(total(total(temporary(echelle), 1), 1) EQ 0) 
     520        endif 
     521      end 
     522      (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin 
     523        e133 = e1[*]#e3 
     524        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
     525          IF keyword_set(wdepth) THEN $ 
     526            e133[bottom] = (e1*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
     527          ELSE e133[bottom] = (e1*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
     528        ENDIF 
     529        echelle = (temporary(e133[*])*temporary(mask[*]))#replicate(1, jpt) 
     530        echelle = reform(echelle, nx, ny, nz, jpt, /over) 
     531        if keyword_set(integration) then divi = 1 ELSE begin  
     532          IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 2) $ 
     533          ELSE divi = total(total(echelle, 1), 2) 
     534        endelse 
     535        res = total(total(temporary(res)*echelle, 1, nan = nan), 2, nan = nan)/(divi > 1) 
     536        if msknan[0] NE -1 then begin 
     537          echelle = temporary(echelle) NE 0 
     538          testnan = temporary(msknan)*echelle 
     539          testnan = total(total(temporary(testnan), 1), 2) $ 
     540            +(total(total(temporary(echelle), 1), 2) EQ 0) 
     541        endif 
     542      end 
     543      (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin 
     544        e233 = e2[*]#e3 
     545        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
     546          IF keyword_set(wdepth) THEN $ 
     547            e233[bottom] = (e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
     548          ELSE e233[bottom] = (e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
     549        ENDIF 
     550        echelle = (temporary(e233[*])*temporary(mask[*]))#replicate(1, jpt) 
     551        echelle = reform(echelle, nx, ny, nz, jpt, /over) 
     552        if keyword_set(integration) then divi = 1 ELSE begin  
     553          IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 2), 2) $ 
     554          ELSE divi = total(total(echelle, 2), 2) 
     555        endelse 
     556        res = total(total(temporary(res)*echelle, 2, nan = nan), 2, nan = nan)/(divi > 1) 
     557        if msknan[0] NE -1 then begin 
     558          echelle = temporary(echelle) NE 0 
     559          testnan = temporary(msknan)*echelle 
     560          testnan = total(total(temporary(testnan), 2), 2) $ 
     561            +(total(total(temporary(echelle), 2), 2) EQ 0) 
     562        endif 
     563      end 
     564      (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin 
     565        e1233 = (e1[*]*e2[*])#e3 
     566        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
     567          IF keyword_set(wdepth) THEN $ 
     568            e1233[bottom] = (e1*e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
     569          ELSE e1233[bottom] = (e1*e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
     570        ENDIF 
     571        echelle = (temporary(e1233[*])*temporary(mask[*]))#replicate(1, jpt) 
     572        echelle = reform(echelle, nx, ny, nz, jpt, /over) 
     573        if keyword_set(integration) then divi = 1 ELSE begin  
     574          IF msknan[0] NE -1 THEN divi = total(total(total(echelle*msknan, 1), 1), 1) $ 
     575          ELSE divi = total(total(total(echelle, 1), 1), 1) 
     576        endelse 
     577        res = total(total(total(temporary(res)*echelle, 1, nan = nan), 1, nan = nan), 1, nan = nan)/(divi > 1) 
     578        if msknan[0] NE -1 then begin 
     579          echelle = temporary(echelle) NE 0 
     580          testnan = temporary(msknan)*echelle 
     581          testnan = total(total(total(temporary(testnan), 1), 1), 1) $ 
     582            +(total(total(total(temporary(echelle), 1), 1), 1) EQ 0) 
     583        endif 
     584      end 
     585    endcase 
     586  endif 
     587;------------------------------------------------------------ 
     588  if dirt EQ 1 AND keyword_set(spatialfirst) then BEGIN 
     589    IF (reverse(size(res, /dimension)))[0] NE jpt THEN BEGIN 
     590      print, 'the last dimension of res is not equal to jpt: '+strtrim(jpt, 2) 
     591      if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
     592      return, -1 
     593    ENDIF 
     594    tdim = size(res, /n_dimensions) 
     595    if keyword_set(integration) then res = total(res, tdim, nan = nan) $ 
     596    ELSE BEGIN  
     597      if keyword_set(nan) then BEGIN  
     598        testnan = testnan < 1 
     599        testnan = total(temporary(testnan), tdim) 
     600        divi = testnan 
     601      ENDIF ELSE divi = jpt 
     602      res = total(res, tdim, nan = nan)/(1 > divi) 
     603    ENDELSE  
     604  ENDIF 
    458605;------------------------------------------------------------ 
    459606;------------------------------------------------------------ 
    … …  
    463610; IV.1) on masque les terres par une valeur a 1e+20 
    464611;------------------------------------------------------------ 
    465    valmask = 1e+20 
    466    terre = where(divi EQ 0) 
    467    IF terre[0] NE -1 THEN BEGIN  
    468       res(temporary(terre)) = 1e+20 
    469    ENDIF  
     612  valmask = 1e+20 
     613  terre = where(divi EQ 0) 
     614  IF terre[0] NE -1 THEN BEGIN  
     615    res(temporary(terre)) = 1e+20 
     616  ENDIF  
    470617;------------------------------------------------------------ 
    471618; IV.2) on remplace, quand nan ne 1, !values.f_nan par nan 
    472619;------------------------------------------------------------ 
    473    if keyword_set(nan) NE 0 then BEGIN  
    474       puttonan = where(temporary(testnan) EQ 0) 
    475       if puttonan[0] NE -1 then res[temporary(puttonan)] = !values.f_nan 
    476       if nan NE 1 then BEGIN  
    477          notanumber = where(finite(res) eq 0) 
    478          if notanumber[0] NE -1 then res[temporary(notanumber)] = nan 
    479       ENDIF 
    480    ENDIF 
     620  if keyword_set(nan) NE 0 then BEGIN  
     621    puttonan = where(temporary(testnan) EQ 0) 
     622    if puttonan[0] NE -1 then res[temporary(puttonan)] = !values.f_nan 
     623    if nan NE 1 then BEGIN  
     624      notanumber = where(finite(res) eq 0) 
     625      if notanumber[0] NE -1 then res[temporary(notanumber)] = nan 
     626    ENDIF 
     627  ENDIF 
    481628;------------------------------------------------------------ 
    482629; IV.3) on se remplace ds le sous domaine qui etait definit a l''entree de 
    483630; moyenne  
    484631;------------------------------------------------------------ 
    485    if keyword_set(boite) AND NOT keyword_set(nodomdef) then domdef, oldboite,GRILLE=vargrid 
    486 ;------------------------------------------------------------ 
    487    if keyword_set(key_performance) THEN print, 'temps grossemoyenne', systime(1)-tempsun  
    488    return, res 
     632  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 
     633;------------------------------------------------------------ 
     634  if keyword_set(key_performance) THEN print, 'temps grossemoyenne', systime(1)-tempsun  
     635  return, res 
    489636;------------------------------------------------------------ 
    490637;------------------------------------------------------------ 

  • trunk/ToBeReviewed/CALCULS/hdyn.pro

    r23 r25  
    101101         case 1 of 
    102102            tailles[1] eq jpi and tailles[2] eq jpj: BEGIN 
    103                sn = tabsn[premierxt:dernierxt, premieryt:dernieryt, *] 
    104                tn = tabtn[premierxt:dernierxt, premieryt:dernieryt, *] 
     103               sn = tabsn[firstxt:lastxt, firstyt:lastyt, *] 
     104               tn = tabtn[firstxt:lastxt, firstyt:lastyt, *] 
    105105            end 
    106106            tailles[1] eq  nx and tailles[2] eq  ny:BEGIN  
    … …  
    116116         e33d = replicate(1, nx*ny)#e3t 
    117117         e33d = reform(e33d, nx, ny, jpk, /over) 
    118          terre = where(tmask[premierxt:dernierxt, premieryt:dernieryt, *] EQ 0) 
     118         terre = where(tmask[firstxt:lastxt, firstyt:lastyt, *] EQ 0) 
    119119         if terre[0] NE -1 then vol[terre] = !values.f_nan 
    120120         case level of 
    … …  
    128128         case 1 of 
    129129            tailles[1] eq jpi and tailles[2] eq jpj AND tailles[4] EQ jpt: BEGIN 
    130                sn = tabsn[premierxt:dernierxt, premieryt:dernieryt, *, *] 
    131                tn = tabtn[premierxt:dernierxt, premieryt:dernieryt, *, *] 
     130               sn = tabsn[firstxt:lastxt, firstyt:lastyt, *, *] 
     131               tn = tabtn[firstxt:lastxt, firstyt:lastyt, *, *] 
    132132            end 
    133133            tailles[1] eq  nx and tailles[2] eq  ny AND tailles[4] EQ jpt:BEGIN  
    … …  
    144144         e33d = e33d[*]#replicate(1, jpt) 
    145145         e33d = reform(e33d, nx, ny, jpk, jpt, /over) 
    146          mask = tmask[premierxt:dernierxt, premieryt:dernieryt, *] 
     146         mask = tmask[firstxt:lastxt, firstyt:lastyt, *] 
    147147         mask = mask[*]#replicate(1, jpt) 
    148148         terre = where(mask EQ 0) 

  • trunk/ToBeReviewed/CALCULS/level2depth.pro

    r23 r25  
    4545;------------------------------------------------------------ 
    4646   niveaux = litchamp(tab) 
    47    grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz 
     47   grille,mask, -1, -1,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz 
    4848;--------------------------------------------------------------- 
    4949; verification de la coherence entre la taille du tableau et le domaine definit par domdef 
    … …  
    5252   if taille[0] NE 2 then return, report('le champ en entree doit contenir un tableau 2d') 
    5353   case 1 of 
    54       taille[1] eq jpi and taille[2] eq jpj:niveaux=niveaux[premierx:dernierx, premiery:derniery] 
     54      taille[1] eq jpi and taille[2] eq jpj:niveaux=niveaux[firstx:lastx, firsty:lasty] 
    5555      taille[1] eq  nx and taille[2] eq  ny: 
    5656      else:return, report('Probleme d''adequation entre les tailles du domaine et celle du champ.') 

  • trunk/ToBeReviewed/CALCULS/level2mask.pro

    r23 r25  
    66; 
    77; PURPOSE: permet de passer d''un tableau 2d de niveau seuil au 
    8 ; tableau 3d de mask avec des 1 ds les niveaux au dessus (et sur) du 
    9 ; niveau seuil et des 0 en dessous. 
     8; tableau 3d de mask avec des 1 ds les niveaux au dessus du 
     9; niveau seuil et des 0 en dessous (et sur). 
    1010; 
    1111; CATEGORY: SANS BOUCLE 
    … …  
    3232; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr) 
    3333;                       17/6/1999 
     34; Setp 2004: boundary level have 0 values and not 1 (as it was 
     35; explained before in the header). see: 
     36; print, array_equal(niveau, total(level2mask(niveau),3)) 
     37; 
    3438;- 
    3539;------------------------------------------------------------ 
    … …  
    4549;------------------------------------------------------------ 
    4650   niveaux = litchamp(tab) 
    47    grille,maskterre,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz 
     51   grille,maskterre, -1, -1, -1,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz 
    4852;--------------------------------------------------------------- 
    4953; verification de la coherence entre la taille du tableau et le domaine definit par domdef 
    … …  
    5357   if taille[0] NE 2 then return, report('le champ en entree doit contenir un tableau 2d') 
    5458   case 1 of 
    55       taille[1] eq jpi and taille[2] eq jpj:niveaux=niveaux[premierx:dernierx, premiery:derniery] 
     59      taille[1] eq jpi and taille[2] eq jpj:niveaux=niveaux[firstx:lastx, firsty:lasty] 
    5660      taille[1] eq  nx and taille[2] eq  ny: 
    5761      else:return, report('Probleme d''adequation entre les tailles du domaine et celle du champ.') 

  • 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;------------------------------------------------------------ 

  • trunk/ToBeReviewed/CALCULS/norme.pro

    r23 r25  
    2121; KEYWORD PARAMETERS: 
    2222; 
    23 ;       BOITE: boite sur laquelle moyenner (par defaut le domaine 
     23;       BOXZOOM: boxzoom sur laquelle moyenner (par defaut le domaine 
    2424;       selectionner par le dernier domdef effectue) 
    2525; 
    … …  
    6060;     pour calculer la moyenne de la norme des courants sur tout le 
    6161;     dommaine entre 0 et 50: 
    62 ;      IDL> res=norme(un,vn,boite=[0,50],dir='xyz') 
     62;      IDL> res=norme(un,vn,boxzoom=[0,50],dir='xyz') 
    6363; 
    6464; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr) 
    … …  
    6868;------------------------------------------------------------ 
    6969;------------------------------------------------------------ 
    70 FUNCTION norme,composanteu, composantev, BOITE=boite,DIREC=direc 
    71 @common 
     70FUNCTION norme, composanteu, composantev, BOXZOOM = boxzoom, DIREC = direc, _extra = ex 
     71;--------------------------------------------------------- 
     72@cm_4mesh 
     73@cm_4data 
     74@cm_4cal 
     75  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     76@updatenew 
     77@updatekwd 
     78  ENDIF 
     79;--------------------------------------------------------- 
    7280   tempsun = systime(1)         ; pour key_performance 
     81; 
     82   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 
     83     return, report(['This version of norme is based on Arakawa C-grid.' $ 
     84                     , 'U and V grids must therefore be defined']) 
     85; 
     86;------------------------------------------------------------ 
     87  if keyword_set(boxzoom) then BEGIN  
     88    Case 1 Of 
     89      N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]] 
     90      N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]] 
     91      N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2] 
     92      N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]] 
     93      N_Elements(Boxzoom) Eq 6:bte = Boxzoom 
     94      Else: return, report('Mauvaise Definition de Boxzoom') 
     95    ENDCASE 
     96    domdef, boxzoom 
     97  ENDIF  
    7398;------------------------------------------------------------ 
    7499   if NOT keyword_set(direc) then direc = 0 
    … …  
    104129; on trouve les points que u et v ont en communs 
    105130;------------------------------------------------------------ 
    106    indicexu = (lindgen(jpi))[premierxu:premierxu+nxu-1] 
    107    indicexv = (lindgen(jpi))[premierxv:premierxv+nxv-1] 
     131   indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 
     132   indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 
    108133   indicex = inter(indicexu, indicexv) 
    109    indiceyu = (lindgen(jpj))[premieryu:premieryu+nyu-1] 
    110    indiceyv = (lindgen(jpj))[premieryv:premieryv+nyv-1] 
     134   indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 
     135   indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 
    111136   indicey = inter(indiceyu, indiceyv) 
    112137   nx = n_elements(indicex)  
    … …  
    124149         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] 
    125150         indice3d = lindgen(jpi, jpj, jpk) 
    126          indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,premierzt:dernierzt] 
     151         indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt] 
    127152;------------------------------------------------------------ 
    128153; extraction de u et v sur le domaine qui convient 
    … …  
    134159                  nzt:BEGIN  
    135160                     if nxu NE nx then $ 
    136                       if indicex[0] EQ premierxu then u = u[0:nx-1,*,*] ELSE u = u[1: nx, *,*]  
     161                      if indicex[0] EQ firstxu then u = u[0:nx-1,*,*] ELSE u = u[1: nx, *,*]  
    137162                     IF nxv NE nx THEN $ 
    138                       if indicex[0] EQ premierxv then v = v[0:nx-1,*,*] ELSE v = v[1: nx, *,*] 
     163                      if indicex[0] EQ firstxv then v = v[0:nx-1,*,*] ELSE v = v[1: nx, *,*] 
    139164                     IF nyu NE ny THEN $ 
    140                       if indicey[0] EQ premieryu then u = u[*,0:ny-1,*] ELSE u = u[*, 1: ny,*]  
     165                      if indicey[0] EQ firstyu then u = u[*,0:ny-1,*] ELSE u = u[*, 1: ny,*]  
    141166                     IF nyv NE ny THEN $ 
    142                       if indicey[0] EQ premieryv then v = v[*,0:ny-1,*] ELSE v = v[*, 1: ny,*] 
     167                      if indicey[0] EQ firstyv then v = v[*,0:ny-1,*] ELSE v = v[*, 1: ny,*] 
    143168                  end 
    144169                  jpk:BEGIN  
    145170                     if nxu NE nx then $ 
    146                       if indicex[0] EQ premierxu then u = u[0:nx-1, *,premierzt:dernierzt] ELSE u = u[1: nx, *,premierzt:dernierzt]  
     171                      if indicex[0] EQ firstxu then u = u[0:nx-1, *,firstzt:lastzt] ELSE u = u[1: nx, *,firstzt:lastzt]  
    147172                     IF nxv NE nx THEN $ 
    148                       if indicex[0] EQ premierxv then v = v[0:nx-1, *,premierzt:dernierzt] ELSE v = v[1: nx, *,premierzt:dernierzt] 
     173                      if indicex[0] EQ firstxv then v = v[0:nx-1, *,firstzt:lastzt] ELSE v = v[1: nx, *,firstzt:lastzt] 
    149174                     IF nyu NE ny THEN $ 
    150                       if indicey[0] EQ premieryu then u = u[*, 0:ny-1,premierzt:dernierzt] ELSE u = u[*, 1: ny,premierzt:dernierzt]  
     175                      if indicey[0] EQ firstyu then u = u[*, 0:ny-1,firstzt:lastzt] ELSE u = u[*, 1: ny,firstzt:lastzt]  
    151176                     IF nyv NE ny THEN $ 
    152                       if indicey[0] EQ premieryv then v = v[*, 0:ny-1,premierzt:dernierzt] ELSE v = v[*, 1: ny,premierzt:dernierzt] 
     177                      if indicey[0] EQ firstyv then v = v[*, 0:ny-1,firstzt:lastzt] ELSE v = v[*, 1: ny,firstzt:lastzt] 
    153178                  end 
    154179                  ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') 
    … …  
    174199         a=u(0,*,*) 
    175200         u=(u+shift(u,1,0,0))/2. 
    176          if NOT keyword_set(key_periodique) OR nx NE jpi then u(0,*,*)=a 
     201         if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*,*)=a 
    177202         a=v(*,0,*) 
    178203         v=(v+shift(v,0,1,0))/2. 
    179          if NOT keyword_set(key_periodique) OR nx NE jpi then v(*,0,*)=a 
     204         if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0,*)=a 
    180205;---------------------------------------------------------------------------- 
    181206; attribution du mask et des tableau de longitude et latitude 
    … …  
    190215         if landv[0] NE -1 then v[landv] = 0 
    191216         res=sqrt(u^2+v^2) 
    192          if NOT keyword_set(key_periodique) OR nx NE jpi then res(0,*, *)=!values.f_nan 
     217         if NOT keyword_set(key_periodic) OR nx NE jpi then res(0,*, *)=!values.f_nan 
    193218         res(*,0, *)=!values.f_nan 
    194219         mask = where(mask eq 0) 
    195220         IF mask[0] NE -1 THEN res(mask) = valmask 
    196221; moyennes en tous genres 
    197          oldboite = [lon1, lon2, lat1, lat2, prof1, prof2] 
    198          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme 
    199          if keyword_set(direc) then res = moyenne(res,direc,/nan, boite = boite) 
     222         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme 
     223         if keyword_set(direc) then res = moyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef) 
    200224;----------------------------------------------------------- 
    201225      END 
    … …  
    215239             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    216240               if nxu NE nx then $ 
    217                 if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]  
     241                if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]  
    218242               IF nxv NE nx THEN $ 
    219                 if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     243                if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    220244               IF nyu NE ny THEN $ 
    221                 if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]  
     245                if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]  
    222246               IF nyv NE ny THEN $ 
    223                 if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     247                if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    224248            END 
    225249            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    … …  
    235259         a=u(0,*,*) 
    236260         u=(u+shift(u,1,0,0))/2. 
    237          if NOT keyword_set(key_periodique) OR nx NE jpi then u(0,*,*)=a 
     261         if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*,*)=a 
    238262         a=v(*,0,*) 
    239263         v=(v+shift(v,0,1,0))/2. 
    240          if NOT keyword_set(key_periodique) OR nx NE jpi then v(*,0,*)=a 
     264         if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0,*)=a 
    241265;---------------------------------------------------------------------------- 
    242266; attribution du mask et des tableau de longitude et latitude 
    … …  
    245269; compte (faire un petit dessin...) 
    246270;---------------------------------------------------------------------------- 
    247          mask = tmask[indice2d+jpi*jpj*(niveau-1)] 
     271         mask = tmask[indice2d+jpi*jpj*firstzt] 
    248272         if ny EQ 1 then mask = reform(mask, nx, ny, /over) 
    249273;----------------------------------------------------------- 
    … …  
    256280         if landv[0] NE -1 then v[landv] = 0 
    257281         res=sqrt(u^2+v^2) 
    258          if NOT keyword_set(key_periodique) OR nx NE jpi then res(0,*, *)=!values.f_nan 
     282         if NOT keyword_set(key_periodic) OR nx NE jpi then res(0,*, *)=!values.f_nan 
    259283         res(*,0, *)=!values.f_nan 
    260284         mask = where(mask eq 0) 
    … …  
    267291         ENDIF 
    268292; moyennes en tous genres 
    269          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme 
    270          if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boite = boite) 
     293         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme 
     294         if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef) 
    271295      END 
    272296;---------------------------------------------------------------------------- 
    … …  
    279303         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] 
    280304         indice3d = lindgen(jpi, jpj, jpk) 
    281          indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,premierzt:dernierzt] 
     305         indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt] 
    282306;------------------------------------------------------------ 
    283307; extraction de u et v sur le domaine qui convient 
    … …  
    289313                  nzt:BEGIN  
    290314                     if nxu NE nx then $ 
    291                       if indicex[0] EQ premierxu then u = u[0:nx-1,*,*,*] ELSE u = u[1: nx, *,*,*]  
     315                      if indicex[0] EQ firstxu then u = u[0:nx-1,*,*,*] ELSE u = u[1: nx, *,*,*]  
    292316                     IF nxv NE nx THEN $ 
    293                       if indicex[0] EQ premierxv then v = v[0:nx-1,*,*,*] ELSE v = v[1: nx, *,*,*] 
     317                      if indicex[0] EQ firstxv then v = v[0:nx-1,*,*,*] ELSE v = v[1: nx, *,*,*] 
    294318                     IF nyu NE ny THEN $ 
    295                       if indicey[0] EQ premieryu then u = u[*,0:ny-1,*,*] ELSE u = u[*, 1: ny,*,*]  
     319                      if indicey[0] EQ firstyu then u = u[*,0:ny-1,*,*] ELSE u = u[*, 1: ny,*,*]  
    296320                     IF nyv NE ny THEN $ 
    297                       if indicey[0] EQ premieryv then v = v[*,0:ny-1,*,*] ELSE v = v[*, 1: ny,*,*] 
     321                      if indicey[0] EQ firstyv then v = v[*,0:ny-1,*,*] ELSE v = v[*, 1: ny,*,*] 
    298322                  end 
    299323                  jpk:BEGIN  
    300324                     if nxu NE nx then $ 
    301                       if indicex[0] EQ premierxu then u = u[0:nx-1, *,premierzt:dernierzt,*] ELSE u = u[1: nx, *,premierzt:dernierzt,*]  
     325                      if indicex[0] EQ firstxu then u = u[0:nx-1, *,firstzt:lastzt,*] ELSE u = u[1: nx, *,firstzt:lastzt,*]  
    302326                     IF nxv NE nx THEN $ 
    303                       if indicex[0] EQ premierxv then v = v[0:nx-1, *,premierzt:dernierzt,*] ELSE v = v[1: nx, *,premierzt:dernierzt,*] 
     327                      if indicex[0] EQ firstxv then v = v[0:nx-1, *,firstzt:lastzt,*] ELSE v = v[1: nx, *,firstzt:lastzt,*] 
    304328                     IF nyu NE ny THEN $ 
    305                       if indicey[0] EQ premieryu then u = u[*, 0:ny-1,premierzt:dernierzt,*] ELSE u = u[*, 1: ny,premierzt:dernierzt,*]  
     329                      if indicey[0] EQ firstyu then u = u[*, 0:ny-1,firstzt:lastzt,*] ELSE u = u[*, 1: ny,firstzt:lastzt,*]  
    306330                     IF nyv NE ny THEN $ 
    307                       if indicey[0] EQ premieryv then v = v[*, 0:ny-1,premierzt:dernierzt,*] ELSE v = v[*, 1: ny,premierzt:dernierzt,*] 
     331                      if indicey[0] EQ firstyv then v = v[*, 0:ny-1,firstzt:lastzt,*] ELSE v = v[*, 1: ny,firstzt:lastzt,*] 
    308332                  end 
    309333                  ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') 
    … …  
    312336            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $ 
    313337             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN  
    314                u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,premierzt:dernierzt, *] 
    315                v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,premierzt:dernierzt, *] 
     338               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt, *] 
     339               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt, *] 
    316340            END 
    317341            ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') 
    … …  
    322346         a=u(0,*,*,*) 
    323347         u=(u+shift(u,1,0,0,0))/2. 
    324          if NOT keyword_set(key_periodique) OR nx NE jpi then u(0,*,*,*)=a 
     348         if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*,*,*)=a 
    325349         a=v(*,0,*,*) 
    326350         v=(v+shift(v,0,1,0,0))/2. 
    327          if NOT keyword_set(key_periodique) OR nx NE jpi then v(*,0,*,*)=a 
     351         if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0,*,*)=a 
    328352;---------------------------------------------------------------------------- 
    329353; attribution du mask et des tableau de longitude et latitude 
    … …  
    338362         if landv[0] NE -1 then v[landv] = 0 
    339363         res=sqrt(u^2+v^2) 
    340          if NOT keyword_set(key_periodique) OR nx NE jpi then res(0,*, *, *)=!values.f_nan 
     364         if NOT keyword_set(key_periodic) OR nx NE jpi then res(0,*, *, *)=!values.f_nan 
    341365         res(*,0, *, *)=!values.f_nan 
    342366         mask = where(mask eq 0) 
    … …  
    349373         ENDIF 
    350374; moyennes en tous genres 
    351          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme 
    352          if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boite = boite) 
     375         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme 
     376         if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef) 
    353377      END 
    354378;---------------------------------------------------------------------------- 
    … …  
    367391             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    368392               if nxu NE nx then $ 
    369                 if indicex[0] EQ premierxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]  
     393                if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]  
    370394               IF nxv NE nx THEN $ 
    371                 if indicex[0] EQ premierxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
     395                if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
    372396               IF nyu NE ny THEN $ 
    373                 if indicey[0] EQ premieryu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]  
     397                if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]  
    374398               IF nyv NE ny THEN $ 
    375                 if indicey[0] EQ premieryv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
     399                if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
    376400            END 
    377401            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    … …  
    394418         a=u(0,*) 
    395419         u=(u+shift(u,1,0))/2. 
    396          if NOT keyword_set(key_periodique) OR nx NE jpi then u(0,*)=a 
     420         if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*)=a 
    397421         a=v(*,0) 
    398422         v=(v+shift(v,0,1))/2. 
    399          if NOT keyword_set(key_periodique) OR nx NE jpi then v(*,0)=a 
     423         if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0)=a 
    400424;---------------------------------------------------------------------------- 
    401425; attribution du mask et des tableau de longitude et latitude 
    … …  
    404428; compte (faire un petit dessin...) 
    405429;---------------------------------------------------------------------------- 
    406          mask = tmask[indice2d+jpi*jpj*(niveau-1)] 
     430         mask = tmask[indice2d+jpi*jpj*firstzt] 
    407431         if nyt EQ 1 THEN mask = reform(mask, nx, ny, /over) 
    408432;----------------------------------------------------------- 
    … …  
    415439         if landv[0] NE -1 then v[landv] = 0 
    416440         res=sqrt(u^2+v^2) 
    417          if NOT keyword_set(key_periodique) OR nx NE jpi then res(0,*)=!values.f_nan 
     441         if NOT keyword_set(key_periodic) OR nx NE jpi then res(0,*)=!values.f_nan 
    418442         res(*,0)=!values.f_nan 
    419443         mask = where(mask eq 0) 
    420444         IF mask[0] NE -1 THEN res(mask) = valmask 
    421445; moyennes en tous genres 
    422          oldboite = [lon1, lon2, lat1, lat2, prof1, prof2] 
    423          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme 
    424          if keyword_set(direc) then res = moyenne(res,direc,/nan, boite = boite) 
     446         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme 
     447         if keyword_set(direc) then res = moyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef) 
    425448      END 
    426449;---------------------------------------------------------------------------- 

  • trunk/ToBeReviewed/CALCULS/projectondepth.pro

    r23 r25  
    4343;   ->champ nul a 1e-6 pres 
    4444;  
    45 ;  verifcation en projettant la temperature sur la profondeur de la 2o 
     45;  verifcation en projettant la temperature sur la profondeur de la 20 
    4646;  degres par exemple... 
    4747; 
    … …  
    6565; verification de la coherence entre la taille du tableau et le 
    6666; domaine  
    67    grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz 
     67   grille, mask, -1, -1, -1,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz 
    6868   case 1 of 
    69       tailledepth[1] eq jpi and tailledepth[2] eq jpj:depth=depth[premierx:dernierx, premiery:derniery] 
     69      tailledepth[1] eq jpi and tailledepth[2] eq jpj:depth=depth[firstx:lastx, firsty:lasty] 
    7070      tailledepth[1] eq  nx and tailledepth[2] eq  ny: 
    7171      else:return, report('Probleme d''adequation entre les tailles du domaine et celle du tableau de profondeur') 
    … …  
    7373   case 1 OF 
    7474      taillearray[3] NE jpk:return, report('Le tableau 3d doit avoir sa 3eme dimension egale a jpk') 
    75       taillearray[1] eq jpi and taillearray[2] eq jpj:array=array[premierx:dernierx, premiery:derniery, *] 
     75      taillearray[1] eq jpi and taillearray[2] eq jpj:array=array[firstx:lastx, firsty:lasty, *] 
    7676      taillearray[1] eq  nx and taillearray[2] eq  ny: 
    7777      else:return, report('Probleme d''adequation entre les tailles du domaine et celle du tableau de profondeur') 
    … …  
    100100   if out[0] NE -1 then res[out] = !values.f_nan 
    101101; on masque les points terres a valmask 
    102    grille,mask 
    103102   if n_elements(valmask) EQ 0 then valmask = 1e20 
    104103   terre = where((temporary(mask))[*, *, 0] EQ 0) 

  • trunk/ToBeReviewed/CALCULS/remplit.pro

    r23 r25  
    1 FUNCTION remplit, zinput, NAN = nan, NITE = nite, BASIQUE = basique, mask = mask, _extra = ex 
     1FUNCTION remplit, zinput, NAN = nan, NITER = niter, BASIQUE = basique, mask = mask, FILLXDIR = fillxdir, FILLYDIR = fillydir, FILLVAL = fillval, _extra = ex 
    22@common 
    3    tempsun = systime(1)         ; pour key_performance 
     3  tempsun = systime(1)          ; pour key_performance 
    44;+ 
    5    ;; 
    6    ;; Extrapole zinout[jpi,jpj] sur les continents en utilisant les 4 
    7    ;; plus proches voisins masques oceaniquement et construit un nouveau 
     5  ;; 
     6  ;; Extrapole zinout[jpi,jpj] sur les continents en utilisant les 4 
     7  ;; plus proches voisins masques oceaniquement et construit un nouveau 
    88;   masque 
    9    ;; contenant l'ancien masque oceanique PLUSles points extrapoles. 
    10    ;; Reitere le processus nite fois. 
    11    ;; C'est pas clair, essayez ! 
    12    ;; 
    13    ;; 
     9  ;; contenant l'ancien masque oceanique PLUSles points extrapoles. 
     10  ;; Reitere le processus niter fois. 
     11  ;; C'est pas clair, essayez ! 
     12  ;; 
     13  ;; 
    1414; 
    1515;    /Nan: to fill the point which have the value 
    … …  
    1818; 
    1919; 
    20  
    21  
    22  
    2320;- 
    2421; les points non remplis sont masques a valmask 
    25    z = zinput 
    26    if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 
    27    if keyword_set(basique) then begin 
    28       oldkey_gridtype = key_gridtype 
    29       key_gridtype = 'c' 
    30    endif 
    31    IF n_elements(nite) EQ 0 THEN nite = 1 
    32    IF nite EQ 0 THEN return, zinput 
    33    if keyword_set(basique)  then begin 
    34       nx = (size(zinput))[1] 
    35       ny = (size(zinput))[2] 
    36       if NOT keyword_set(mask) then mmmask = basique ELSE mmmask = mask 
    37       if  key_gridtype eq 'e' then begin 
    38          case vargrid of 
    39             'T':glam = glamt[premierxt:dernierxt, premieryt:dernieryt] 
    40             'U':glam = glamu[premierxu:dernierxu, premieryu:dernieryu] 
    41          endcase 
    42       endif 
    43    ENDIF ELSE grille, mmmask, glam, gphi, gdep, nx, ny, nz, _extra = ex 
    44    if keyword_set(mask) then mmmask = mask 
    45 ;--------------------------------------------------------------- 
    46    if (size(mmmask))[0] EQ 3 THEN BEGIN  
    47       if (size(mmmask))[3] EQ jpk THEN mmmask = mmmask[*, *, niveau-1] $ 
    48       ELSE  mmmask = mmmask[*, *, nz-1] 
    49    ENDIF 
    50 ; 
    51    if n_elements(mmmask) EQ 1 then mmmask = replicate(1, nx, ny) 
    52    if keyword_set(nan) then begin 
    53       nanpoint = where(finite(z) EQ 0) 
    54       if nanpoint[0] NE -1 then begin 
    55          mmmask[nanpoint] = 0 
    56          z[nanpoint] = 0 
    57       endif 
    58    endif 
     22  IF n_elements(niter) EQ 0 THEN niter = 1 
     23  IF niter EQ 0 THEN return, zinput 
     24  z = zinput 
     25  if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 
     26  if keyword_set(basique) then begin 
     27    oldkey_gridtype = key_gridtype 
     28    key_gridtype = 'c' 
     29    nx = (size(zinput))[1] 
     30    ny = (size(zinput))[2] 
     31    if NOT keyword_set(mask) then mmmask = basique ELSE mmmask = mask 
     32    if  key_gridtype eq 'e' then begin 
     33      case vargrid of 
     34        'T':glam = glamt[firstxt:lastxt, firstyt:lastyt] 
     35        'U':glam = glamu[firstxu:lastxu, firstyu:lastyu] 
     36      endcase 
     37    endif 
     38  ENDIF ELSE grille, mmmask, glam, gphi, gdep, nx, ny, nz, _extra = ex 
     39  if keyword_set(mask) then mmmask = mask 
     40;--------------------------------------------------------------- 
     41  if (size(mmmask))[0] EQ 3 THEN mmmask = mmmask[*, *, 0] 
     42; 
     43  if n_elements(mmmask) EQ 1 then mmmask = replicate(1b, nx, ny) 
     44  if keyword_set(nan) then begin 
     45    nanpoint = where(finite(z) EQ 0) 
     46    if nanpoint[0] NE -1 then begin 
     47      mmmask[nanpoint] = 0b 
     48      z[nanpoint] = 0 
     49    endif 
     50  ENDIF 
     51  mmmask = byte(mmmask) 
    5952;--------------------------------------------------------------- 
    6053; on ajoute un cadre de zero a z, mask, e1, e2 
    … …  
    6255; soucier des bords du domaine! 
    6356;--------------------------------------------------------------- 
    64    tempdeux = systime(1)        ; pour key_performance =2 
    65    nx2 = nx+2 
    66    case key_gridtype of 
    67       'c':BEGIN 
    68          zero = fltarr(nx+2, ny+2) 
    69          zero[1:nx, 1:ny] = mmmask 
    70          mmmask = zero 
    71          zero = fltarr(nx+2, ny+2) 
    72          zero[1:nx, 1:ny] = z 
    73          if keyword_set(key_periodique) AND nx EQ jpi then begin 
    74             zero[0, 1:ny] = z[jpi-1, *] 
    75             zero[nx+1, 1:ny] = z[0, *] 
    76          endif 
    77          z = zero 
    78       END 
    79       'e':BEGIN 
    80          zero = fltarr(nx+2, ny+4) 
    81          zero[1:nx, 2:ny+1] = mmmask 
    82          mmmask = zero 
    83          zero = fltarr(nx+2, ny+4) 
    84          zero[1:nx, 2:ny+1] = z 
    85          if keyword_set(key_periodique) AND nx EQ jpi then begin 
    86             zero[0, 2:ny+1] = z[jpi-1, *] 
    87             zero[nx+1, 2:ny+1] = z[0, *] 
    88          endif 
    89          z = zero 
    90       END 
    91    endcase 
    92    IF testvar(var = key_performance) EQ 2 THEN $ 
     57  tempdeux = systime(1)         ; pour key_performance =2 
     58  nx2 = nx+2 
     59  case key_gridtype of 
     60    'c':BEGIN 
     61      ztmp = bytarr(nx+2, ny+2) 
     62      ztmp[1:nx, 1:ny] = mmmask 
     63      mmmask = temporary(ztmp) 
     64      ztmp = fltarr(nx+2, ny+2) 
     65      ztmp[1:nx, 1:ny] = z 
     66      if keyword_set(key_periodic) AND nx EQ jpi then begin 
     67        ztmp[0, 1:ny] = z[jpi-1, *] 
     68        ztmp[nx+1, 1:ny] = z[0, *] 
     69      endif 
     70      z = temporary(ztmp) 
     71    END 
     72    'e':BEGIN 
     73      ztmp = bytarr(nx+2, ny+4) 
     74      ztmp[1:nx, 2:ny+1] = mmmask 
     75      mmmask = temporary(ztmp) 
     76      ztmp = fltarr(nx+2, ny+4) 
     77      ztmp[1:nx, 2:ny+1] = z 
     78      if keyword_set(key_periodic) AND nx EQ jpi then begin 
     79        ztmp[0, 2:ny+1] = z[jpi-1, *] 
     80        ztmp[nx+1, 2:ny+1] = z[0, *] 
     81      endif 
     82      z = temporary(ztmp) 
     83    END 
     84  endcase 
     85  IF testvar(var = key_performance) EQ 2 THEN $ 
    9386    print, 'temps remplit: on ajoute un cadre de zero ', systime(1)-tempdeux 
    9487;--------------------------------------------------------------- 
    … …  
    9790;--------------------------------------------------------------- 
    9891;--------------------------------------------------------------- 
    99    FOR n = 1, nite DO BEGIN  
    100 ; on trouve les points cote 
    101       tempdeux = systime(1)     ; pour key_performance =2 
    102 ; les points du bord du cadre ne doivent pas etre dans la cote 
    103       case key_gridtype of 
    104          'c':BEGIN 
    105             mmmask[0, *] = 1 
    106             mmmask[nx+1, *] = 1 
    107             mmmask[*, 0] = 1 
    108             mmmask[*, ny+1] = 1 
    109             if keyword_set(key_periodique) AND nx EQ jpi then begin 
    110                mmmask[0,*] = mmmask[nx-2, *] 
    111                mmmask[nx+1,*] = mmmask[1, *] 
    112             endif 
    113          END 
    114          'e':BEGIN 
    115             mmmask[0, *] = 1 
    116             mmmask[nx+1, *] = 1 
    117             mmmask[*, 0:1] = 1 
    118             mmmask[*, ny+2:ny+3] = 1 
    119             if keyword_set(key_periodique) AND nx EQ jpi then begin 
    120                mmmask[0,*] = mmmask[nx-2, *] 
    121                mmmask[nx+1,*] = mmmask[1, *] 
    122             endif 
    123          END 
    124       endcase 
     92  FOR n = 1, niter DO BEGIN  
     93; on trouve les points coast 
     94    tempdeux = systime(1)       ; pour key_performance =2 
     95; les points du bord du cadre ne doivent pas etre selectionnes comme 
     96; la coast 
     97    case key_gridtype of 
     98      'c':BEGIN 
     99        mmmask[0, *] = 1b 
     100        mmmask[nx+1, *] = 1b 
     101        mmmask[*, 0] = 1b 
     102        mmmask[*, ny+1] = 1b 
     103      END 
     104      'e':BEGIN 
     105        mmmask[0, *] = 1b 
     106        mmmask[nx+1, *] = 1b 
     107        mmmask[*, 0:1] = 1b 
     108        mmmask[*, ny+2:ny+3] = 1b 
     109      END 
     110    endcase 
    125111; liste des points terre restant 
    126       terre = where(mmmask EQ 0) 
    127       if terre[0] EQ -1 then GOTO, fini 
     112    IF keyword_set(fillxdir) THEN BEGIN 
     113; we stop if all the lines, that contains data, have been filled 
     114      test = total(mmmask[1:nx, *], 1) 
     115      IF total((test EQ 0)+(test EQ nx)) EQ ny+2 THEN GOTO, fini 
     116    ENDIF 
     117    IF keyword_set(fillydir) THEN BEGIN 
     118; we stop if all the columns, that contains data, have been filled 
     119      test = total(mmmask[*, 1:ny], 2) 
     120      IF total((test EQ 0)+(test EQ ny)) EQ nx+2 THEN GOTO, fini 
     121    ENDIF 
     122    land = where(mmmask EQ 0) 
     123    if land[0] EQ -1 then GOTO, fini 
    128124; les points du bord du cadre doivent maintenant etre dans la terre 
    129       case key_gridtype of 
    130          'c':BEGIN 
    131             mmmask[0, *] = 0 
    132             mmmask[nx+1, *] = 0 
    133             mmmask[*, 0] = 0 
    134             mmmask[*, ny+1] = 0 
    135          END 
    136          'e':BEGIN 
    137             mmmask[0, *] = 0 
    138             mmmask[nx+1, *] = 0 
    139             mmmask[*, 0:1] = 0 
    140             mmmask[*, ny+2:ny+3] = 0 
    141          END 
    142       endcase 
     125    case key_gridtype of 
     126      'c':BEGIN 
     127        mmmask[0, *] = 0b 
     128        mmmask[nx+1, *] = 0b 
     129        mmmask[*, 0] = 0b 
     130        mmmask[*, ny+1] = 0b 
     131      END 
     132      'e':BEGIN 
     133        mmmask[0, *] = 0b 
     134        mmmask[nx+1, *] = 0b 
     135        mmmask[*, 0:1] = 0b 
     136        mmmask[*, ny+2:ny+3] = 0b 
     137      END 
     138    endcase 
     139    if keyword_set(key_periodic) AND nx EQ jpi then begin 
     140      mmmask[0, *] = mmmask[nx, *] 
     141      mmmask[nx+1, *] = mmmask[1, *] 
     142    endif 
    143143; liste des voisins mer 
    144       case key_gridtype of 
    145          'c':BEGIN 
    146             voisins = shift(mmmask,-1,0)+shift(mmmask,1,0)+shift(mmmask,0,-1)+shift(mmmask,0,1) $ 
    147              +1./sqrt(2)*(shift(mmmask,-1,-1)+shift(mmmask,1,1)+shift(mmmask,1,-1)+shift(mmmask,-1,1)) 
    148             cote = where(voisins[terre] GT 0) 
    149 ; liste des points cote 
    150             cote = terre[cote] 
    151             poids = voisins[cote] 
    152          END 
    153          'e':BEGIN 
    154             shifted = glam[0, 0] LT glam[0, 1] 
    155             oddeven = (terre/nx2+1-shifted) MOD 2 
    156             voisins = mmmask[terre+1]+mmmask[terre-1] $ 
    157              +mmmask[2*nx2+terre]+mmmask[-2*nx2+terre] $ 
    158              +sqrt(2)*(mmmask[terre-nx2+oddeven]+mmmask[terre+nx2+oddeven] $ 
    159                        +mmmask[terre+(-nx2-1)+oddeven]+mmmask[terre+(nx2-1)+oddeven]) 
    160             cote = where(voisins GT 0) 
    161             poids = voisins[cote] 
    162             cote = terre[cote] 
    163          END 
    164       endcase 
    165 ; 
    166       IF testvar(var = key_performance) EQ 2 THEN $ 
    167        print, 'temps remplit: trouver la cote ', systime(1)-tempdeux 
    168 ;--------------------------------------------------------------- 
    169 ; remplissage des points cote 
    170 ;--------------------------------------------------------------- 
    171       tempdeux = systime(1)     ; pour key_performance =2 
     144    case key_gridtype of 
     145      'c':BEGIN 
     146        CASE 1 OF 
     147          keyword_set(fillxdir):weight = mmmask[1+land]+mmmask[-1+land] 
     148          keyword_set(fillydir):weight = mmmask[nx2+land]+mmmask[-nx2+land] 
     149          ELSE:weight = mmmask[1+land]+mmmask[-1+land]+mmmask[nx2+land]+mmmask[-nx2+land] $ 
     150            +1./sqrt(2)*(mmmask[nx2+1+land]+mmmask[nx2-1+land] $ 
     151                         +mmmask[-nx2+1+land]+mmmask[-nx2-1+land]) 
     152        ENDCASE 
     153      END 
     154      'e':BEGIN 
     155        shifted = glam[0, 0] LT glam[0, 1] 
     156        oddeven = (land/nx2+1-shifted) MOD 2 
     157        weight = mmmask[1+land]+mmmask[-1+land] $ 
     158          +mmmask[2*nx2+land]+mmmask[-2*nx2+land] $ 
     159          +sqrt(2)*(mmmask[-nx2+oddeven+land]+mmmask[-nx2-1+oddeven+land] $ 
     160                    +mmmask[nx2+oddeven+land]+mmmask[nx2-1+oddeven+land]) 
     161      END 
     162    endcase 
     163 
     164    ok = where(weight GT 0) 
     165    weight = weight[ok] 
     166    coast = land[temporary(ok)] 
     167; 
     168    IF testvar(var = key_performance) EQ 2 THEN $ 
     169      print, 'temps remplit: trouver la coast ', systime(1)-tempdeux 
     170;--------------------------------------------------------------- 
     171; remplissage des points coast 
     172;--------------------------------------------------------------- 
     173    tempdeux = systime(1)       ; pour key_performance =2 
    172174; on masque z 
    173       z = z*mmmask 
    174 ; 
    175       case key_gridtype of 
    176          'c':BEGIN 
    177             zcote = z[1+cote]+z[-1+cote]+z[nx2+cote]+z[-nx2+cote] $ 
    178              +1./sqrt(2)*(z[nx2+1+cote]+z[nx2-1+cote] $ 
    179                           +z[-nx2+1+cote]+z[-nx2-1+cote]) 
    180          END 
    181          'e':BEGIN 
    182             oddeven = (cote/nx2+1-shifted) MOD 2 
    183             zcote = z[1+cote]+z[-1+cote]+z[2*nx2+cote]+z[-2*nx2+cote] $ 
    184              +sqrt(2)*(z[-nx2+oddeven+cote]+z[-nx2-1+oddeven+cote] $ 
    185                        +z[nx2+oddeven+cote]+z[nx2-1+oddeven+cote]) 
    186          END 
    187       endcase 
    188        
    189       z[cote] = zcote/poids 
     175    z = temporary(z)*mmmask 
     176; 
     177    case key_gridtype of 
     178      'c':BEGIN 
     179        CASE 1 OF 
     180          keyword_set(fillxdir):zcoast = z[1+coast]+z[-1+coast] 
     181          keyword_set(fillydir):zcoast = z[nx2+coast]+z[-nx2+coast] 
     182          ELSE:zcoast = z[1+coast]+z[-1+coast]+z[nx2+coast]+z[-nx2+coast] $ 
     183            +1./sqrt(2)*(z[nx2+1+coast]+z[nx2-1+coast] $ 
     184                         +z[-nx2+1+coast]+z[-nx2-1+coast]) 
     185        ENDCASE 
     186      END 
     187      'e':BEGIN 
     188        oddeven = (coast/nx2+1-shifted) MOD 2 
     189        zcoast = z[1+coast]+z[-1+coast]+z[2*nx2+coast]+z[-2*nx2+coast] $ 
     190          +sqrt(2)*(z[-nx2+oddeven+coast]+z[-nx2-1+oddeven+coast] $ 
     191                    +z[nx2+oddeven+coast]+z[nx2-1+oddeven+coast]) 
     192      END 
     193    endcase 
     194;     
     195    z[coast] =  temporary(zcoast)/ temporary(weight) 
     196; we update the the boundary conditions of z 
     197    if keyword_set(key_periodic) AND nx EQ jpi then begin 
     198      z[0, *] = z[nx, *] 
     199      z[nx+1, *] = z[1, *] 
     200    endif 
    190201;--------------------------------------------------------------- 
    191202; IV) on reduit le masque 
    192203;--------------------------------------------------------------- 
    193       mmmask[cote] = 1 
    194 ; 
    195       IF testvar(var = key_performance) EQ 2 THEN $ 
    196        print, 'temps remplit: une iteration ', systime(1)-tempdeux 
    197    ENDFOR  
     204    mmmask[ temporary(coast)] = 1 
     205; 
     206    IF testvar(var = key_performance) EQ 2 THEN $ 
     207      print, 'temps remplit: une iteration ', systime(1)-tempdeux 
     208  ENDFOR  
    198209fini: 
    199210; 
    200 ; on masque les valeurs sur les terres restantes 
    201 ; 
    202    IF n_elements(valmask) EQ 0 then valmask = 1e20 
    203    z = z*mmmask+valmask*(1-mmmask) 
     211; on masque les valeurs sur les lands restantes 
     212; 
     213  IF n_elements(valmask) EQ 0 then valmask = 1e20 
     214  IF n_elements(fillval) EQ 0 THEN fillval = valmask 
     215  z = temporary(z)*mmmask + fillval*(1b-mmmask) 
    204216;--------------------------------------------------------------- 
    205217; on redecoupe le tableau pour retirer le cadre! 
    206218;--------------------------------------------------------------- 
    207    case key_gridtype of 
    208       'c':BEGIN 
    209          z = z[1:nx, 1:ny] 
    210       END 
    211       'e':BEGIN 
    212          z = z[1:nx, 2:ny+1] 
    213       END 
    214    endcase 
    215 ; 
    216    if keyword_set(basique) then key_gridtype = oldkey_gridtype 
    217 ;--------------------------------------------------------------- 
    218    if keyword_set(key_performance) THEN print, 'temps remplit', systime(1)-tempsun  
    219    return, z 
     219  case key_gridtype of 
     220    'c':BEGIN 
     221      z = z[1:nx, 1:ny] 
     222    END 
     223    'e':BEGIN 
     224      z = z[1:nx, 2:ny+1] 
     225    END 
     226  endcase 
     227; 
     228  if keyword_set(basique) then key_gridtype = oldkey_gridtype 
     229;--------------------------------------------------------------- 
     230  if keyword_set(key_performance) THEN print, 'temps remplit', systime(1)-tempsun  
     231  return, z 
    220232END  
    221233 
Note: See TracChangeset for help on using the changeset viewer.