;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME: moyenne
;
; PURPOSE:  averages a 2- or 3-d field over a selected
;           geographical area and along one ore more
;	    selected axes (x, y or z)
;
; CATEGORY:
;
; CALLING SEQUENCE: result = moyenne(tab,'direc',BOXZOOM=boxzoom)
;
; INPUTS:	tab = 2 or 3d field
;		direc = 'x' 'y' 'z' 'xy' 'xz' 'yz' or 'xyz'
;
; KEYWORD PARAMETERS:
;
;		BOXZOOM = [xmin,xmax,ymin,ymax (,(zmin,)zmax)] pour plus
;		de detail cf domdef.
;               boxzoom peut prendre 5 formes: 
;               [vert2], [vert1, vert2],[lon1, lon2, lat1, lat2],  
;               [lon1, lon2, lat1, lat2, vert2],[lon1, lon2, lat1,
;               lat2, vert1,vert2]
;
;               NAN: not a number, a activer si l''on peut faire veut
;               faire une moyenne sans tenir compte de certaines
;               valeurs masques de tab.
;               si les valeurs masques de tab sont la valeur consacree
;               par IDL (!values.f_nan), il suffit de mettre /NAN.
;               si les valeurs masques de tab on pour valeur a (a doit
;               etre differente de 1 (correspond a nan =
;               !values.f_nan) et de 0 (qui desactive nan)
;               il faut mettre NAN=a. 
;               Rq: en sorties les points de result qui sont NAN
;               auront pour valeur a ou !values.f_nan.
;
;               NODOMDEF: activer si l''on ne veut pas passer ds
;               domdef bien que le mot cle boxzoom soit present (comme
;               c''est le cas qd moyenne est appelee via checkfield) 
;            
;               INTEGRATION: pour faire une integrale plutot qu''une moyenne
; 
;         /WDEPTH: to specify that the field is at W depth instad of T 
;         depth (automatically activated if vargrid eq 'W')
;            
; OUTPUTS:  result:un tableau 
;
; COMMON BLOCKS:
;	common   domdef
;
; SIDE EFFECTS:met les valeurs correspondants a la terre a 1e20
;
; RESTRICTIONS:
;
; EXAMPLE:
;
; MODIFICATION HISTORY: Jerome Vialard (jv@lodyc.jussieu.fr)
;			2/7/98
;                       Sebastien Masson (smasson@lodyc.jussieu.fr)
;                       mise au propre de certains truc (les terres...)
;                       14/8/98
;                       15/1/98
;                       11/3/99 adaptation pour NAN
;                       28/7/99 moyennes tableaux 1d
;-
;------------------------------------------------------------
;------------------------------------------------------------
; PLAN DU PROGRAMME:
;------------------------------------------------------------
;------------------------------------------------------------
;I) preliminaires
;   I.1) determination des directions de moyennes d''apres direc
;   I.2) verification de la taille du tableau d''entree
;   I.3) obtention des facteurs d''echelles et du masque sur le sous
;   domaine concerne par la moyenne 
;  ) moyennes pour les tableaux 1d (x,y)
;II) moyennes pour les tableaux 2d (x,y)
;   II.1) verification de la coherence de la taille du tableau a
;   moyenner 
;   II.2) differents types de moyennes possibles
;III) moyennes pour les tableaux 3d (x,y,z)
;   III.1) verification de la coherence de la taille du tableau a
;   moyenner 
;   III.2) differents types de moyennes possibles
;IV ) finitions
;   IV.1) on masque les terres par une valeur a 1e+20
;   IV.2) on remplace, quand nan ne 1, !values.f_nan par nan
;   IV.3) on revient au sous domaine initial.
;------------------------------------------------------------
;------------------------------------------------------------
function moyenne, tab, direc, BOXZOOM = boxzoom, INTEGRATION = integration $
                  , NAN = nan, NODOMDEF = nodomdef, WDEPTH = wdepth $
                  , _extra = ex
;---------------------------------------------------------
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
@cm_4cal
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
@updatekwd
  ENDIF
;---------------------------------------------------------
  tempsun = systime(1)          ; pour key_performance
;------------------------------------------------------------
;I) preliminaires
;------------------------------------------------------------
  dirt = 0
  dirx = 0
  diry = 0
  dirz = 0
;------------------------------------------------------------
; I.1) direction(s) suivants lesquelles on integre
;------------------------------------------------------------
  if ( strpos(direc, 't') ge 0 ) then dirt = 1
  if ( strpos(direc, 'x') ge 0 ) then dirx = 1
  if ( strpos(direc, 'y') ge 0 ) then diry = 1
  if ( strpos(direc, 'z') ge 0 ) then dirz = 1
  if (dirx eq 0 and diry eq 0 and dirz eq 0) then return, tab
;------------------------------------------------------------
; I.2) verification de la taille du tableau d'entree
;------------------------------------------------------------
  taille = size(tab)
  case 1 of
    taille[0] eq 1 :dim = '1d'
    taille[0] eq 2 :BEGIN
      dim = '2d'
      if dirx eq 0 and diry eq 0 then return, tab
    END
    taille[0] eq 3 :BEGIN
      dim = '3d'
      if dirx eq 0 and diry eq 0 and dirz eq 0 then return, tab
    END
    else : return, report('Le tableau d''entree doit etre a 2 ou 3 dimensions s''il contient une dim temporelle utiliser grossemoyenne')
  endcase
;------------------------------------------------------------
;   I.3) obtention des facteurs d''echelles et du masque sur le sous
;   domaine concerne par la moyenne 
; redefinition du domaine ajuste a boxzoom (a 6 elements)
; ceci va nous permetre de faire les calcules que sur le sous domaine
; comcerne par la moyenne. domdef, suivit de grille nous donne tous
; les tableaux de la grille sur le sous domaine 
;------------------------------------------------------------
  if keyword_set(boxzoom) then BEGIN 
    Case 1 Of
      N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]]
      N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]]
      N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2]
      N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]]
      N_Elements(Boxzoom) Eq 6:bte = Boxzoom
      Else: return, report('Mauvaise Definition de Boxzoom')
    endcase
    if NOT keyword_set(nodomdef) then BEGIN 
      savedbox = 1b
      saveboxparam, 'boxparam4moyenne.dat'
      domdef, bte, GRIDTYPE = vargrid, _extra = ex
    ENDIF 
  ENDIF 
;---------------------------------------------------------------
; attribution du mask et des tableaux de longitude et latitude...
;---------------------------------------------------------------
  IF vargrid EQ 'W' THEN wdepth = 1
  grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, WDEPTH = wdepth
;------------------------------------------------------------
;------------------------------------------------------------
; II) Cas du tableau 1d
;------------------------------------------------------------
;------------------------------------------------------------
  if dim EQ '1d' then BEGIN
    if n_elements(tab) NE nx*ny AND n_elements(tab) NE nx*ny*nz then BEGIN
      if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
      return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.')
    ENDIF
    case 1 of
      nx EQ 1 AND ny EQ 1:BEGIN ;vecteur suivant z
        case n_elements(tab) of
          jpk:res = tab[firstz:lastz]
          nz:res = tab
          ELSE:BEGIN 
            if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
            return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.')
          END
        ENDCASE
        if dirz EQ 1 then BEGIN 
          dim = '3d' 
          taille = size(reform(res, nx, ny, nz))
        ENDIF ELSE BEGIN
          if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
          return, res
        ENDELSE 
      END
      ny EQ 1:BEGIN             ;vecteur suivant x
        case n_elements(tab) of
          jpi:res = tab[firstx:lastx]
          nx:res = tab
          ELSE:BEGIN 
            if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
            return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.')
          END
        ENDCASE
        if dirx EQ 1 then BEGIN 
          dim = '2d' 
          taille = size(reform(res, nx, ny))
        ENDIF ELSE BEGIN 
          if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
          return, res
        ENDELSE
      END
      nx EQ 1:BEGIN             ;vecteur suivant y
        case n_elements(tab) of
          jpj:res = tab[firsty:lasty]
          ny:res = tab
          ELSE:BEGIN 
            if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
            return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.')
          END
        ENDCASE
        if diry EQ 1 then BEGIN
          dim = '2d'
          taille = size(reform(res, nx, ny))
        ENDIF ELSE BEGIN 
          if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
          return, res
        ENDELSE 
      END
    endcase
  endif
;------------------------------------------------------------
;------------------------------------------------------------
; II) Cas du tableau 2d
;------------------------------------------------------------
;------------------------------------------------------------
  if (dim eq '2d') then begin
;---------------------------------------------------------------
;   II.1) verification de la coherence de la taille du tableau a
;   moyenner 
; verification de la coherence entre la taille du tableau et le
; domaine definit par domdef
; le tableau en entree doit avoir soit la taille du domaine total
; (jpi,jpj) soit celle du domaine reduit (nx,ny)
;---------------------------------------------------------------
    case 1 of
      taille[1] eq jpi and taille[2] eq jpj: $
        res = tab[firstx:lastx, firsty:lasty]
      taille[1] eq  nx and taille[2] eq  ny:res = tab
      else:BEGIN 
        if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
        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))
      END
    ENDCASE
    if keyword_set(nan) NE 0 then BEGIN 
      if nan NE 1 then BEGIN    ; si nan n''est pas !values.f_nan 
; on le met a !values.f_nan
        if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $
        ELSE notanumber = where(abs(res) GT abs(nan)/10.)
        if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan
      ENDIF
    ENDIF
;---------------------------------------------------------------
; rq IL FAUT FAIRE ATTENTION AUX CAS OU LA DIM A MOYENNER = 1, ET
; S'ASSURER QU'ELLE EXISTE BIEN. D'OU LES reform(...,nx,ny,...) QUI
; PEUVENT SEMBLER INUTILE AU DEPART
;---------------------------------------------------------------
    if nx EQ 1 OR ny EQ 1 then BEGIN 
      res = reform(res, nx, ny, /over)
      e1 =  reform(e1, nx, ny, /over)
      e2 =  reform(e2, nx, ny, /over)
    endif
    if nx EQ 1 OR ny EQ 1 OR nz EQ 1 THEN $
      mask =  reform(mask, nx, ny, nz, /over)
;---------------------------------------------------------------
; II.3) differents types de moyennes
;---------------------------------------------------------------
    mask = mask[*, *, 0]
    if keyword_set(nan) NE 0 then msknan = finite(res) ELSE msknan = -1
    case 1 of
      (dirx eq 1) and (diry eq 0) : begin
        e = e1*mask
        if keyword_set(integration) then divi = 1 $
        else begin 
          divi = e
          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan
          if ny EQ 1 then divi = reform(divi, nx, ny, /over)
          divi = total(divi, 1)
        endelse
        res = res*e
        if ny EQ 1 then res = reform(res, nx, ny, /over)
        res = total(res, 1, nan = nan)/(divi > 1.)
        if msknan[0] NE -1 then begin
          testnan  = msknan*mask
          if ny EQ 1 then testnan = reform(testnan, nx, ny, /over)
          testnan = total(testnan, 1)+(total(mask, 1) EQ 0)
        endif
      end
      (dirx eq 0) and (diry eq 1) : begin
        e = e2*mask
        if keyword_set(integration) then divi = 1 $
        else begin 
          divi = e
          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan
          if ny EQ 1 then divi = reform(divi, nx, ny, /over)
          divi = total(divi, 2)
        endelse
        res = res*e
        if ny EQ 1 then res = reform(res, nx, ny, /over)
        res = total(res, 2, nan = nan)/(divi > 1.)
        if msknan[0] NE -1 then begin
          testnan  = msknan*mask
          if ny EQ 1 then testnan  = reform(testnan, nx, ny, /over)
          testnan = total(testnan, 2)+(total(mask, 2) EQ 0)
        endif
      end
      (dirx eq 1) and (diry eq 1) : begin
        if keyword_set(integration) then divi = 1 else BEGIN 
          IF msknan[0] NE -1 THEN divi = total(e1*e2*mask*msknan) $
          ELSE divi = total(e1*e2*mask)
        ENDELSE 
        res = total(res*e1*e2*mask, nan = nan)/(divi > 1.)
        if msknan[0] NE -1 then begin
          testnan  = msknan*mask
          testnan = total(testnan)+(total(mask) EQ 0)
        endif
      end
    endcase
  endif
;------------------------------------------------------------
;------------------------------------------------------------
; III) Cas du tableau 3d
;------------------------------------------------------------
;------------------------------------------------------------
  if (dim eq '3d') then begin
;---------------------------------------------------------------
; III.1) verification de la coherence de la taille du tableau a
; moyenner 
; verification de la coherence entre la taille du tableau et le
; domaine definit par domdef 
; le tableau en entree doit avoir soit la taille du domaine total
; (jpi,jpj,jpk) soit celle du domaine reduit (nx,ny,ny)
;---------------------------------------------------------------
    case 1 of
      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk: $
        res = tab[firstx:lastx, firsty:lasty, firstz:lastz]
      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz: $
        res = tab[firstx:lastx, firsty:lasty, *]
      taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq nz :res = tab
      taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk : $
        res = tab[*, *, firstz:lastz]
      else:BEGIN 
        if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
        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))
      END
    endcase
    if keyword_set(nan) NE 0 then BEGIN 
      if nan NE 1 then BEGIN    ; si nan n''est pas !values.f_nan 
; on le met a !values.f_nan
        if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $
        ELSE notanumber = where(abs(res) GT abs(nan)/10.)
        if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan
      ENDIF
    ENDIF
;---------------------------------------------------------------
; rq IL FAUT FAIRE ATTENTION AUX CAS OU LA DIM A MOYENNER = 1, ET
; S'ASSURER QU'ELLE EXISTE BIEN. D'OU LES reform(...,nx,ny,...) QUI
; PEUVENT SEMBLER INUTILE AU DEPART
;---------------------------------------------------------------
    if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then BEGIN 
      res = reform(res, nx, ny, nz, /over)
      e1 =  reform(e1, nx, ny, /over)
      e2 =  reform(e2, nx, ny, /over)
    endif
    if nx EQ 1 OR ny EQ 1 OR nz EQ 1 THEN $
      mask =  reform(mask, nx, ny, nz, /over)
    IF keyword_set(key_partialstep) THEN BEGIN
; the top of the ocean floor is
      IF vargrid EQ 'T' OR vargrid EQ 'W' THEN bottom = total(mask, 3) $
      ELSE bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 
; we suppress columns with only ocean or land 
      good = where(bottom NE 0 AND bottom NE nz)
; the bottom of the ocean in 3D index is:
      bottom = lindgen(nx*ny)+(temporary(bottom)-1L)*nx*ny
      IF good[0] NE -1 THEN bottom = bottom[good] $
      ELSE bottom = -1
    ENDIF ELSE bottom = -1
;---------------------------------------------------------------
; III.2) differents types de moyennes
;---------------------------------------------------------------
    if keyword_set(nan) NE 0 then msknan = finite(res) ELSE msknan = -1
    case 1 of
      (dirx eq 1) and (diry eq 0) and (dirz eq 0) : begin
        e13 = e1[*]#replicate(1., nz)
        e13 = reform(e13, nx, ny, nz, /over)
        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $
          AND nx NE 1 THEN BEGIN
          IF msknan[0] EQ -1 THEN BEGIN 
            msknan = replicate(1b, nx, ny, nz)
            nan = 1
          endif
          msknan[bottom] = 0
          res[bottom] = !values.f_nan
        ENDIF
        if keyword_set(integration) then divi = 1 else begin
          divi = e13*mask
          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan
          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over)
          divi = total(divi, 1)
        ENDELSE
        res = res*e13*mask
        if nz EQ 1 then res = reform(res, nx, ny, nz, /over)
        res = total(res, 1, nan = nan)/(divi > 1.)
        e13 = 1
        if msknan[0] NE -1 then begin
          testnan  = msknan*mask
          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over)
          testnan = total(testnan, 1)+(total(mask, 1) EQ 0)
        endif
      end
      (dirx eq 0) and (diry eq 1) and (dirz eq 0) : begin
        e23 = e2[*]#replicate(1., nz)
        e23 = reform(e23, nx, ny, nz, /over)
        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $
          AND ny NE 1 THEN BEGIN
          IF msknan[0] EQ -1 THEN BEGIN 
            msknan = replicate(1b, nx, ny, nz)
            nan = 1
          endif
          msknan[bottom] = 0
          res[bottom] = !values.f_nan
        ENDIF
        if keyword_set(integration) then divi = 1 else begin
          divi = e23*mask
          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan
          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over)
          divi = total(divi, 2)
        ENDELSE
        res = res*e23*mask
        if nz EQ 1 then res = reform(res, nx, ny, nz, /over)
        res = total(res, 2, nan = nan)/(divi > 1.)
        e23 = 1
        if msknan[0] NE -1 then begin
          testnan  = msknan*mask
          if nz EQ 1 then testnan = reform(testnan, nx, ny, nz, /over)
          testnan = total(testnan, 2)+(total(mask, 2) EQ 0)
        endif
      end
      (dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin
        e33 = replicate(1, 1.*nx*ny)#e3
        e33 = reform(e33, nx, ny, nz, /over)
        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN
          IF keyword_set(wdepth) THEN $
            e33[bottom] = (e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $
          ELSE e33[bottom] = (e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)]
        ENDIF
        if keyword_set(integration) then divi = 1 else begin
          divi = e33*mask
          if msknan[0] NE -1 then divi = temporary(divi)*msknan
          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over)
          divi = total(divi, 3)
        ENDELSE
        res = res*e33*mask
        if nz EQ 1 then res = reform(res, nx, ny, nz, /over)
        res = total(res, 3, nan = nan)/(divi > 1.)
        e33 = 1
        if msknan[0] NE -1 then begin
          testnan  = msknan*mask
          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over)
          testnan = total(testnan, 3)+(total(mask, 3) EQ 0)
        endif
      end
      (dirx eq 1) and (diry eq 1) and (dirz eq 0) : begin
        e123 = (e1*e2)[*]#replicate(1., nz)
        e123 = reform(e123, nx, ny, nz, /over)
        IF keyword_set(key_partialstep) AND bottom[0] NE -1 $
          AND nx*ny NE 1 THEN BEGIN
          IF msknan[0] EQ -1 THEN BEGIN 
            msknan = replicate(1b, nx, ny, nz)
            nan = 1
          endif
          msknan[bottom] = 0
          res[bottom] = !values.f_nan
        ENDIF
        if keyword_set(integration) then divi = 1 else BEGIN 
          divi = e123*mask
          IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan
          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over)
          divi = total(total(divi, 1), 1)
        ENDELSE
        res = res*e123*mask
        if nz EQ 1 then res = reform(res, nx, ny, nz, /over)
        res = total(total(res, 1, nan = nan), 1, nan = nan) / (divi > 1.)
        e123 = 1
        if msknan[0] NE -1 then begin
          testnan  = msknan*mask
          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over)
          testnan = total(total(testnan, 1), 1)+(total(total(mask, 1), 1) EQ 0)
        endif
      end
      (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin
        e133 = e1[*]#e3
        e133 = reform(e133, nx, ny, nz, /over)
        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN
          IF keyword_set(wdepth) THEN $
            e133[bottom] = (e1*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $
          ELSE e133[bottom] = (e1*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)]
        ENDIF
        if keyword_set(integration) then divi = 1 else BEGIN 
          divi = e133*mask
          if msknan[0] NE -1 then divi = temporary(divi)*msknan
          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over)
          divi = total(total(divi, 1), 2)
        ENDELSE
        res = res*e133*mask
        if nz EQ 1 then res = reform(res, nx, ny, nz, /over)
        res = total(total(res, 1, nan = nan), 2, nan = nan) / (divi > 1.)
        e133 = 1
        if msknan[0] NE -1 then begin
          testnan  = msknan*mask
          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over)
          testnan = total(total(testnan, 1), 2)+(total(total(mask, 1), 2) EQ 0)
        endif
      end
      (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin
        e233 = e2[*]#e3
        e233 = reform(e233, nx, ny, nz, /over)
        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN
          IF keyword_set(wdepth) THEN $
            e233[bottom] = (e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $
          ELSE e233[bottom] = (e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)]
        ENDIF
        if keyword_set(integration) then divi = 1 else BEGIN 
          divi = e233*mask
          if msknan[0] NE -1 then divi = temporary(divi)*msknan
          if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over)
          divi = total(total(divi, 2), 2)
        ENDELSE
        res = res*e233*mask
        if nz EQ 1 then res = reform(res, nx, ny, nz, /over)
        res = total(total(res, 2, nan = nan), 2, nan = nan) / (divi > 1.)
        e233 = 1
        if msknan[0] NE -1 then begin
          testnan  = msknan*mask
          if nz EQ 1 then testnan  = reform(testnan, nx, ny, nz, /over)
          testnan = total(total(testnan, 2), 2)+(total(total(mask, 2), 2) EQ 0)
        endif
      end
      (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin
        e1233 = (e1*e2)[*]#e3
        e1233 = reform(e1233, nx, ny, nz, /over)
        IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN
          IF keyword_set(wdepth) THEN $
            e1233[bottom] = (e1*e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $
          ELSE e1233[bottom] = (e1*e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)]
        ENDIF
        if keyword_set(integration) then divi = 1 else BEGIN 
          if msknan[0] NE -1 then divi = total(e1233*mask*msknan) $
          ELSE divi = total(e1233*mask)
        ENDELSE
        res = total(res*e1233*mask, nan = nan) / (divi > 1.)
        e1233 = 1
        if msknan[0] NE -1 then begin
          testnan  = msknan*mask
          testnan = total(testnan)+(total(mask) EQ 0)
        endif
      end
    endcase
  endif
;------------------------------------------------------------
;------------------------------------------------------------
;IV ) finitions
;------------------------------------------------------------
;------------------------------------------------------------
; IV.1) on masque les terres par une valeur a 1e+20
;------------------------------------------------------------
  valmask = 1e+20
  terre = where(divi EQ 0)
  IF terre[0] NE -1 THEN BEGIN
    res[terre] = 1e+20
  ENDIF  
;------------------------------------------------------------
; IV.2) on remplace, quand nan ne 1, !values.f_nan par nan
;------------------------------------------------------------
  if keyword_set(nan) NE 0 then BEGIN 
    puttonan = where(testnan EQ 0)
    if puttonan[0] NE -1 then res[puttonan] = !values.f_nan
    if nan NE 1 then BEGIN 
      notanumber = where(finite(res) eq 0)
      if notanumber[0] NE -1 then res[notanumber] = nan
    ENDIF
  ENDIF
;------------------------------------------------------------
; IV.3) on se remplace ds le sous domaine qui etait definit a l''entree de
; moyenne 
;------------------------------------------------------------
  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat'
;------------------------------------------------------------
  if keyword_set(key_performance) THEN print, 'temps moyenne', systime(1)-tempsun 
  return, res
;------------------------------------------------------------
;------------------------------------------------------------
end