;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME: grossemoyenne
;
; PURPOSE:  averages a 3- or 4-d time serie field over a selected
;           geographical area or along the time axis. For one ore more
;	    selected axes (x, y, z, t)
;
; CATEGORY:
;
; CALLING SEQUENCE: result = grossemoyenne(tab,'direc',BOITE=boite)
;
; INPUTS:	tab = 3 or 4d field
;		direc ='x' 'y' 'z' 't' 'xy' 'xz' 'yz' 'xyz' 'xt' 'yt' 'zt'
;		'xyt' 'xzt' 'yzt' or 'xyzt' 
;
; KEYWORD PARAMETERS:
;
;		boite = [xmin,xmax,ymin,ymax (,(zmin,)zmax)] pour plus
;		de detail cf domdef
;               boite peut prendre 5 formes: 
;               [prof2], [prof1, prof2],[lon1, lon2, lat1, lat2],  
;               [lon1, lon2, lat1, lat2, prof2],[lon1, lon2, lat1, lat2, prof1,prof2]
;             
;               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 boite soit present (comme
;               c''est le cas qd grossemoyenne est appelee via
;               checkfield)
;
;               INTEGRATION: pour faire une integrale plutot qu''une moyenne
;            
; OUTPUTS: 
;
; COMMON BLOCKS:   result:un tableau 
;	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)
; adaptation pour les tableaux comportants une dimension temporelle
;                       14/8/98
;                       15/1/98
;                       12/3/99 adaptation pour NAN et utilisation de TEMPORARY
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
; 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 
;II) moyennes pour les tableaux 3d (x,y,t)
;   II.1) verification de la coherence de la taille du tableau a
;   moyenner 
;   II.2) renvoie sur moyenne qd une moyenne sur t est demandee
;   II.3) differents types de moyennes possibles
;III) moyennes pour les tableaux 4d (x,y,z,t)
;   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 grossemoyenne ,tab,direc,BOITE=boite, INTEGRATION=integration $
                        , NAN = nan, NODOMDEF = nodomdef, _extra = ex
@common
   tempsun = systime(1)         ; pour key_performance
;------------------------------------------------------------
;I) preliminaires
;------------------------------------------------------------
   dirt = 0
   dirx = 0
   diry = 0
   dirz = 0
   dim  = 'aa'
;------------------------------------------------------------
; 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
;------------------------------------------------------------
; I.2) verification de la taille du tableau d'entree
;------------------------------------------------------------
   taille = size(tab)
   case 1 of
      taille[0] eq 1 :return,  report('Le tableau n''a qu''une dimension, cas non traite!')
      taille[0] eq 2 :return,  report('Le tableau n''a qu''deux dimension, cas non traite!')
      taille[0] eq 3 :BEGIN 
         dim = '3d'
         if dirx eq 0 and diry eq 0 and dirt eq 0 then return, tab
      END
      taille[0] eq 4 :BEGIN 
         dim = '4d'
         if dirx eq 0 and diry eq 0 and dirz eq 0 and dirt eq 0 then return, tab
      END
      else : return, report('Le tableau d entree doit etre a 3 ou 4 dimensions s''il ne contient pas de dim temporelle utiliser moyenne')
   endcase
;------------------------------------------------------------
;   I.4) obtention des facteurs d''echelles et du masque sur le sous
;   domaine concerne par la moyenne 
; redefinition du domaine ajuste a boite (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(boite) then BEGIN 
      Case 1 Of
         N_Elements(Boite) Eq 1: bte = [lon1, lon2, lat1, lat2, 0.,boite[0]]
         N_Elements(Boite) Eq 2: bte = [lon1, lon2, lat1, lat2, boite[0],boite[1]]
         N_Elements(Boite) Eq 4: bte = [Boite,prof1, prof2]
         N_Elements(Boite) Eq 5: bte = [Boite[0:3], 0, Boite[4]]
         N_Elements(Boite) Eq 6: bte = Boite
         Else: return, report('Mauvaise Definition de Boite')
      endcase
      oldboite = [lon1, lon2, lat1, lat2, prof1, prof2]
      if NOT keyword_set(nodomdef) then domdef, bte,GRILLE=vargrid, _extra = ex
   ENDIF 
;---------------------------------------------------------------
; attribution du mask et des tableaux de longitude et latitude...
;---------------------------------------------------------------
   grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz, dernierx, derniery, dernierz,e1,e2,e3
;------------------------------------------------------------
; I.3) si dirt eq 1 on fait la moyenne temporelle et on envoie ds moyenne
;------------------------------------------------------------
   if dirt EQ 1 then begin
      if dim EQ 3d then BEGIN 
      case 1 of
         taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $
          res = tab[premierx:premierx+nx-1 $
                    ,premiery:premiery+ny-1, *]
         taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpt:res=tab
         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))
      ENDCASE
         divi = finite(res)
         divi = total(temporary(divi),3, nan = keyword_set(nan))
         notanum = where(divi EQ 0)
         if keyword_set(integration) then begin
            res = total(res,3, nan = nan)
         ENDIF ELSE BEGIN
            if keyword_set(nan) then BEGIN
               res = total(res,3, nan = keyword_set(nan))/ (1 > divi)
            ENDIF ELSE res = total(res,3)/(1.*taille[3]) 
         ENDELSE
      ENDIF ELSE BEGIN
      case 1 of
         taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $
          res =tab[premierx:dernierx,premiery:derniery,premierz:dernierz, *]
         taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $
          res =tab[premierx:dernierx,premiery:derniery,*, *]
         taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq nz and taille[4] eq jpt:res=tab
         taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk and taille[4] eq jpt: $
          res=tab[*, *, premierz:dernierz, *]
         else:return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $
                             +strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+'*' $
                             +strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*' $
                             +strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)+'*' $
                             +strtrim(taille[4], 1))
      endcase
         divi = finite(res)
         divi = total(temporary(divi),4, nan = keyword_set(nan))
         notanum = where(divi EQ 0)
         if keyword_set(integration) then begin
            res = total(res,4, nan = nan)
         ENDIF ELSE BEGIN
            if keyword_set(nan) then begin
               res = total(res,4, /nan)/(1 > divi)
            ENDIF ELSE res = total(res,4)/(1.*taille[4])
         ENDELSE
      ENDELSE
      if notanum[0] NE -1 then res[temporary(notanum)] = !values.f_nan
      return,  moyenne(temporary(res), direc, BOITE = boite, NAN = nan, INTEGRATION = integration, _extra = ex)
   endif
;------------------------------------------------------------
;------------------------------------------------------------
; II) Cas serie tableaux 2d (tab3d)
;------------------------------------------------------------
;------------------------------------------------------------
   if (dim eq '3d') 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,jpt) soit celle du domaine reduit (nx,ny,jpt)
;---------------------------------------------------------------
      case 1 of
         taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $
          res = tab[premierx:premierx+nx-1 $
                    ,premiery:premiery+ny-1, *]
         taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpt:res=tab
         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))
      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 jpt EQ 1 then res=reform(res,nx,ny,jpt, /over)
      if nx EQ 1 OR ny EQ 1 then BEGIN 
         mask =  reform(mask, nx, ny, nz, /over)
         e1 =  reform(e1, nx, ny, /over)
         e2 =  reform(e2, nx, ny, /over)
      endif
;---------------------------------------------------------------
; II.3) differents types de moyennes
;---------------------------------------------------------------
      if keyword_set(nan) NE 0 then begin
         msknan = replicate(1., nx, ny, jpt)
         notanumber = where(finite(res) EQ 0)
         if notanumber[0] NE -1 then msknan[temporary(notanumber)] = !values.f_nan
      ENDIF ELSE msknan = 1
      if nz eq jpk then mask = mask[*, *,niveau-1] ELSE mask = mask[*, *,nz-1]
      case 1 of
         (dirx eq 1) and (diry eq 0) : begin
            e = temporary(e1)*temporary(mask)
            if nx EQ 1 OR ny EQ 1 then e = reform(e, nx, ny, /over)
            echelle = temporary(e[*])#replicate(1, jpt)
            echelle = reform(echelle, nx, ny, jpt, /over)
            if keyword_set(integration) then divi=1 $
            ELSE divi = total(echelle*msknan,1, nan = nan)
            res=total(temporary(res)*echelle,1, nan = nan)/(divi > 1.)
            if keyword_set(nan) then begin
               testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0)
               testnan = total(testnan,1)
            endif
         end
         (dirx eq 0) and (diry eq 1) : begin
            e = temporary(e2)*temporary(mask)
            if nx EQ 1 OR ny EQ 1 then e = reform(e, nx, ny, /over)
            echelle = temporary(e[*])#replicate(1, jpt)
            echelle = reform(echelle, nx, ny, jpt, /over)
            if keyword_set(integration) then divi=1 $
            ELSE divi = total(echelle*msknan,2, nan = nan)
            res=total(temporary(res)*echelle,2, nan = nan)/(divi > 1.)
            if keyword_set(nan) then begin
               testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0)
               testnan = total(testnan,2)
            endif
         end
         (dirx eq 1) and (diry eq 1) : begin
            echelle=(temporary(e1)*temporary(e2)*temporary(mask))[*]#replicate(1,jpt)
            echelle=reform(echelle,nx,ny,jpt, /over)
            if keyword_set(integration) then divi=1 $
            ELSE divi = total(temporary(total(echelle*msknan,1, nan = nan)),1, nan = nan)
            res = total(temporary(total(temporary(res)*echelle,1, nan=nan)),1, nan=nan)/(divi > 1.)
            if keyword_set(nan) then begin
               testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0)
               testnan = total(total(testnan,1),1)
            endif
         end
      endcase
   endif
;------------------------------------------------------------
;------------------------------------------------------------
; III) Cas serie tableaux 3d (tab4d)
;------------------------------------------------------------
   if (dim eq '4d') 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,jpt) soit celle du domaine reduit (nx,ny,ny,jpt)
;---------------------------------------------------------------
      case 1 of
         taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $
          res =tab[premierx:dernierx,premiery:derniery,premierz:dernierz, *]
         taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $
          res =tab[premierx:dernierx,premiery:derniery,*, *]
         taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq nz and taille[4] eq jpt:res=tab
         taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpk and taille[4] eq jpt: $
          res=tab[*, *, premierz:dernierz, *]
         else:return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $
                             +strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+'*' $
                             +strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*' $
                             +strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)+'*' $
                             +strtrim(taille[4], 1))
      endcase
      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)
      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 OR jpt EQ 1 then res=reform(res,nx,ny,nz,jpt, /over)
      if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then mask =  reform(mask, nx, ny, nz, /over)
;---------------------------------------------------------------
; III.2) differents types de moyennes
;---------------------------------------------------------------
      IF keyword_set(nan) NE 0 then begin
         msknan = replicate(1., nx, ny, nz, jpt)
         notanumber = where(finite(res) EQ 0)
         if notanumber[0] NE -1 then msknan[temporary(notanumber)] = !values.f_nan
      ENDIF ELSE msknan = 1
      case 1 of
         (dirx eq 1) and (diry eq 0) and (dirz eq 0) : begin
            e13 = temporary(e1[*])#replicate(1.,nz)
            e13 = reform(e13,nx,ny,nz, /over)
            echelle = (temporary(e13)*temporary(mask))[*]#replicate(1, jpt)
            echelle = reform(echelle, nx, ny, nz, jpt, /over)
            if keyword_set(integration) then divi=1 $
            ELSE begin 
               if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle
               divi = total(temporary(divi),1, nan = nan)
            endelse
            res = temporary(res)*echelle
            res = total(temporary(res),1, nan = nan)/(divi > 1)
            if keyword_set(nan) then begin
               testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0)
               testnan = total(testnan,1)
            endif
         end
         (dirx eq 0) and (diry eq 1) and (dirz eq 0) : begin
            e23 = temporary(e2[*])#replicate(1.,nz)
            e23 = reform(e23,nx,ny,nz, /over)
            echelle = (temporary(e23)*temporary(mask))[*]#replicate(1, jpt)
            echelle = reform(echelle, nx, ny, nz, jpt, /over)
            if keyword_set(integration) then divi=1 $
            ELSE begin 
               if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle
               divi = total(temporary(divi), 2, nan = nan)
            endelse
            res = total(temporary(res)*echelle,2, nan = nan)/(divi > 1)
            if keyword_set(nan) then begin
               testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0)
               testnan = total(testnan,2)
            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)
            echelle =(temporary(e33)*temporary(mask))[*]#replicate(1, jpt)
            echelle = reform(echelle, nx, ny, nz, jpt, /over)
            if keyword_set(integration) then divi=1 $
            ELSE begin 
               if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle
               divi = total(temporary(divi), 3, nan = nan)
            endelse
            res = total(temporary(res)*echelle,3, nan = nan)/(divi > 1)
            if keyword_set(nan) then begin
               testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0)
               testnan = total(testnan,3)
            endif
         end
         (dirx eq 1) and (diry eq 1) and (dirz eq 0) : begin
            e13 = e1[*]#replicate(1.,nz)
            e13 = reform(e13,nx,ny,nz, /over)
            e23 = e2[*]#replicate(1.,nz)
            e23 = reform(e23,nx,ny,nz, /over)
            echelle = (temporary(e13)*temporary(e23)*temporary(mask))[*]#replicate(1, jpt)
            echelle = reform(echelle,nx,ny,nz,jpt, /over)
            if keyword_set(integration) then divi=1 $
            ELSE begin 
               if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle
               divi = total(total(temporary(divi),1, nan = nan),1, nan = nan)
            endelse
            res =total(total(temporary(res)*echelle,1, nan = nan),1, nan = nan)/(divi > 1)
            if keyword_set(nan) then begin
               testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0)
               testnan = total(total(testnan,1),1)
            endif
         end
         (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin
            e133 = e1[*]#e3
            echelle = (temporary(e133[*])*temporary(mask[*]))#replicate(1, jpt)
            echelle = reform(echelle,nx,ny,nz,jpt, /over)
            if keyword_set(integration) then divi=1 $
            ELSE begin 
               if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle
               divi = total(total(temporary(divi),1, nan = nan),2, nan = nan)
            endelse
            res = total(total(temporary(res)*echelle,1, nan = nan),2, nan = nan)/(divi > 1)
            if keyword_set(nan) then begin
               testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0)
               testnan = total(total(testnan,1),2)
            endif
         end
         (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin
            e233 = e2[*]#e3
            echelle = (temporary(e233[*])*temporary(mask[*]))#replicate(1, jpt)
            echelle = reform(echelle,nx,ny,nz,jpt, /over)
            if keyword_set(integration) then divi=1 $
            ELSE begin 
               if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle
               divi = total(total(temporary(divi),2, nan = nan),2,nan=nan)
            endelse
            res =total(total(temporary(res)*echelle,2, nan = nan),2, nan = nan)/(divi > 1)
            if keyword_set(nan) then begin
               testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0)
               testnan = total(total(testnan,2),2)
            endif
         end
         (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin
            e1233 = (e1[*]*e2[*])#e3
            echelle=(temporary(e1233[*])*temporary(mask[*]))#replicate(1,jpt)
            echelle=reform(echelle,nx,ny,nz,jpt, /over)
            if keyword_set(integration) then divi=1 $
            ELSE begin 
               if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle
               divi = total(total(total(temporary(divi), 1,nan=nan), 1,nan=nan), 1,nan=nan)
            endelse
            res =total(total(total(temporary(res)*echelle, 1,nan=nan), 1,nan=nan), 1,nan=nan)/(divi > 1)
            if keyword_set(nan) then begin
               testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0)
               testnan = total(total(total(testnan,1),1),1)
            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(temporary(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(temporary(testnan) EQ 0)
      if puttonan[0] NE -1 then res[temporary(puttonan)] = !values.f_nan
      if nan NE 1 then BEGIN 
         notanumber = where(finite(res) eq 0)
         if notanumber[0] NE -1 then res[temporary(notanumber)] = nan
      ENDIF
   ENDIF
;------------------------------------------------------------
; IV.3) on se remplace ds le sous domaine qui etait definit a l''entree de
; moyenne 
;------------------------------------------------------------
   if keyword_set(boite) AND NOT keyword_set(nodomdef) then domdef, oldboite,GRILLE=vargrid
;------------------------------------------------------------
   if keyword_set(key_performance) THEN print, 'temps grossemoyenne', systime(1)-tempsun 
   return, res
;------------------------------------------------------------
;------------------------------------------------------------
end