Index: trunk/ToBeReviewed/CALCULS/curl.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/curl.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/curl.pro	(revision 25)
@@ -0,0 +1,321 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:curl
+;
+; PURPOSE:calcule la composante verticale du rotationnel d''un champ
+; de vecteur horizontaux
+;
+; CATEGORY:calcule sur les matrices
+;
+; CALLING SEQUENCE:res=curl(u,v)
+;
+; INPUTS:
+;       u et v deux matrices representant les coordonnes d''un
+;       champ de vecteur
+;
+; KEYWORD PARAMETERS:
+;
+; OUTPUTS:res: une matrice 2d
+;
+; COMMON BLOCKS:
+;	common.pro
+;
+; SIDE EFFECTS:
+;
+; RESTRICTIONS:
+; les matrices u et v peuvent de 2 a 4 dimensions.
+; attention pour distinger les differents configurations de u et v
+; (xy, xyz, xyt, xyzt), on regarde la variable du common 
+;        -time qui contient le calendrier en jour julien d''IDL auquel
+;        se rapportent u et v ansi que la variable 
+;        -jpt qui est le nombre de pas de temps a considerer ds time.
+; les tableaux u et v sont decoupes sur le meme domaine
+; geographique. A cause du decalage des grilles T, U, V et F il est
+; possiible que ces 2 tableaux n''aient pas la meme taille et se
+; repportent a des indices differents. Si tel est le cas les tableaux
+; sont redecoupes sur les indices qu'ils ont en commun et le dommaine
+; est redefinit pour qu'il colle a ces indices communs.
+; pour eviter ces redecoupes utiliser le mot cles /memeindice ds
+; domdef.pro
+;
+;
+; les points sur le bord du dessin sont mis a !values.f_nan 
+;
+; EXAMPLE:
+;
+; MODIFICATION HISTORY:Guillaume Roullet (grlod@ipsl.jussieu.fr)
+; 
+;                      Sebastien Masson (smasson@lodyc.jussieu.fr)
+;                      adaptation pour marcher avec un domaine reduit
+;
+;                      21/5/1999: valeurs manquantes a !values.f_nan
+;periodicite 
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION curl, uu, vv
+@common
+   tempsun = systime(1)         ; pour key_performance
+;
+   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
+     return, report(['This version of curl is based on Arakawa C-grid.' $
+                     , 'U and V grids must therefore be defined'])
+;
+   u = litchamp(uu)
+   v = litchamp(vv)
+
+   date1 = time[0]
+   if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
+   if (size(u))[0] NE (size(v))[0] then return,  -1
+
+;------------------------------------------------------------
+; on trouve les points que u et v ont en communs
+;------------------------------------------------------------
+   indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
+   indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
+   indicex = inter(indicexu, indicexv)
+   indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
+   indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
+   indicey = inter(indiceyu, indiceyv)
+   nx = n_elements(indicex) 
+   ny = n_elements(indicey)
+   case 1 of
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xyz
+;----------------------------------------------------------------------------
+      (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN 
+         indice2d = lindgen(jpi, jpj)
+         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
+;------------------------------------------------------------
+; extraction de u et v sur le domaine qui convient
+;------------------------------------------------------------
+         case 1 of
+            (size(u))[0] NE 3 OR (size(v))[0] NE 3: return,  -1
+            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
+             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
+            case 1 of
+               nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
+               nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
+               nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
+               nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
+               ELSE :
+            endcase
+            END
+            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
+             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
+               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
+               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
+            END
+            ELSE:return,  -1
+         endcase
+;------------------------------------------------------------
+; calcul du rotationnel
+;------------------------------------------------------------
+         coefu = (e1u[indice2d])[*]#replicate(1, nzt)
+         coefu = reform(coefu, nx, ny, nzt, /over)
+         coefu = coefu*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
+         terreu = where(coefu EQ 0)
+         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
+         
+         coefv = (e2v[indice2d])[*]#replicate(1, nzt)
+         coefv = reform(coefv, nx, ny, nzt, /over)
+         coefv = coefv*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
+         terrev = where(coefv EQ 0)
+         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
+
+         tabf = (fmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
+         div = (e1f[indice2d]*e2f[indice2d])[*]#replicate(1, nzt)
+         div = reform(div, nx, ny, nzt, /over)
+         tabf = tabf/div
+;
+         zu = u*temporary(coefu)
+         zv = v*temporary(coefv)
+
+         psi =  (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
+         psi = tabf*psi
+;------------------------------------------------------------
+; mise a !values.f_nan de la bordure
+;------------------------------------------------------------
+         if NOT keyword_set(key_periodic)  OR nx NE jpi then begin
+            psi(0, *, *) = !values.f_nan
+            psi(nx-1, *, *) = !values.f_nan
+         endif
+         psi(*, 0, *) = !values.f_nan
+         psi(*, ny-1, *) = !values.f_nan
+;
+         if n_elements(valmask) EQ 0 THEN valmask = 1e20
+         terref =  where(tabf EQ 0)
+         if terref[0] NE -1 then psi[temporary(terref)] = valmask
+;------------------------------------------------------------
+; pour le trace graphique
+;------------------------------------------------------------
+         domdef, (glagmt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
+         if keyword_set(direc) then psi = moyenne(psi,direc,/nan)
+
+;----------------------------------------------------------------------------
+      END
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xyt
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+      date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN 
+         indice2d = lindgen(jpi, jpj)
+         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
+;------------------------------------------------------------
+; extraction de u et v sur le domaine qui convient
+;------------------------------------------------------------
+         case 1 of
+            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
+             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
+               if nxu NE nx then $
+                if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
+               IF nxv NE nx THEN $
+                if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
+               IF nyu NE ny THEN $
+                if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
+               IF nyv NE ny THEN $
+                if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
+            END
+            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
+             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
+               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
+               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
+            END
+            ELSE:BEGIN 
+               print, 'problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs'
+               return, -1
+            end
+         endcase
+;----------------------------------------------------------------------------
+; calcul du rotationnel
+;----------------------------------------------------------------------------
+         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
+         terreu = where(coefu EQ 0)
+         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
+         coefu = temporary(coefu[*])#replicate(1, jpt)
+         coefu = reform(coefu, nx, ny, jpt, /over)
+;
+         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
+         terrev = where(coefv EQ 0)
+         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
+         coefv = temporary(coefv[*])#replicate(1, jpt)
+         coefv = reform(coefv, nx, ny, jpt, /over)
+;
+         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
+         tabf = temporary(tabf[*])#replicate(1, jpt)
+         tabf = reform(tabf, nx, ny, jpt, /over)
+;------------------------------------------------------------
+; calcul du rotationnel
+;------------------------------------------------------------
+         zu = u*temporary(coefu)
+         zv = v*temporary(coefv)
+;
+         psi =  (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
+         psi = tabf*psi
+;------------------------------------------------------------
+; mise a !values.f_nan de la bordure
+;------------------------------------------------------------
+         if NOT keyword_set(key_periodic) OR nx NE jpi then begin
+            psi(0, *, *) = !values.f_nan
+            psi(nx-1, *, *) = !values.f_nan
+         endif
+         psi(*, 0, *) = !values.f_nan
+         psi(*, ny-1, *) = !values.f_nan
+         if n_elements(valmask) EQ 0 THEN valmask = 1e20
+         terref =  where(tabf EQ 0)
+         if terref[0] NE -1 then psi[temporary(terref)] = valmask
+;----------------------------------------------------------------------------
+         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
+         if keyword_set(direc) then psi = grossemoyenne(psi,direc,/nan)
+;----------------------------------------------------------------------------
+      END
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xyzt
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+      date1 NE date2 AND (size(u))[0] EQ 4:BEGIN 
+         return, report('non code!')
+      END
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xy
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+      ELSE:BEGIN                ;xy
+         indice2d = lindgen(jpi, jpj)
+         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
+;------------------------------------------------------------
+; extraction de u et v sur le domaine qui convient
+;------------------------------------------------------------
+         case 1 of
+            (size(u))[0] NE 2 OR (size(v))[0] NE 2: return,  -1
+            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
+             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
+               if nxu NE nx then $
+                if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
+               IF nxv NE nx THEN $
+                if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
+               IF nyu NE ny THEN $
+                if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
+               IF nyv NE ny THEN $
+                if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
+            END
+            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
+             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
+               u = u[indice2d]
+               v = v[indice2d]
+            END
+            ELSE:return,  -1
+         endcase
+;------------------------------------------------------------
+; calcul du rotationnel
+;------------------------------------------------------------
+         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
+         terreu = where(coefu EQ 0)
+         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
+         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
+         terrev = where(coefv EQ 0)
+         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
+         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
+;
+         zu = u*temporary(coefu)
+         zv = v*temporary(coefv)
+
+         psi =  (shift(zv, -1, 0)-zv) + (zu-shift(zu, 0, -1))
+         psi = tabf*psi
+
+;------------------------------------------------------------
+; mise a !values.f_nan de la bordure
+;------------------------------------------------------------
+         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
+            psi(0, *) = !values.f_nan
+            psi(nx-1, *) = !values.f_nan
+         endif
+         psi(*, 0) = !values.f_nan
+         psi(*, ny-1) = !values.f_nan
+;
+         if n_elements(valmask) EQ 0 THEN valmask = 1e20
+         terref =  where(tabf EQ 0)
+         if terref[0] NE -1 then psi[temporary(terref)] = valmask
+;------------------------------------------------------------
+; pour le trace graphique
+;------------------------------------------------------------
+         domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
+         if keyword_set(direc) then psi = moyenne(psi,direc,/nan)
+
+      END
+;----------------------------------------------------------------------------
+   endcase
+;------------------------------------------------------------
+   if keyword_set(key_performance) THEN print, 'temps curl', systime(1)-tempsun 
+;------------------------------------------------------------
+   vargrid = 'F'
+   varname = 'vorticity'
+   return, psi
+end
Index: trunk/ToBeReviewed/CALCULS/depth2floatlevel.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/depth2floatlevel.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/depth2floatlevel.pro	(revision 25)
@@ -0,0 +1,74 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:depth2floatlevel
+;
+; PURPOSE: assez comparable a depth2level mais ici le niveau calcule
+; est en float. Par ex le niveau 5.4 correspond a une profondeur egale
+; a gdep[5]+.4*(gdep[6]-gdep[5])
+;
+; CATEGORY: SANS BOUCLE
+;
+; CALLING SEQUENCE:res=depth2floatlevel(depth2d)
+; 
+; INPUTS: depth2d tableau 2d de profondeur (ou une structure repondant
+; aux criteres de litchamp)
+;
+; KEYWORD PARAMETERS: 
+;       /NOMASK: pour ne pas masquer les points terres
+;
+; OUTPUTS: un tableau 2d de float contenant les valeurs des niveaux.
+;
+; COMMON BLOCKS:common.pro
+;
+; SIDE EFFECTS:accepte les vcaleurs a !values.f_nan et masque les
+; points terres a valmask.
+;
+; RESTRICTIONS:
+;
+; EXAMPLE:
+;
+;    IDL> a=(jpk-1)/(1.*jpi*jpj)*findgen(jpi,jpj)
+;    IDL> plt, 1e6*(a-floatlevel2depth(depth2floatlevel(a))),/nocontour
+;
+;  ->champ nul a 1e-6 pres
+;
+; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
+;                      15/06/2000
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION depth2floatlevel, tab, NOMASK = nomask
+   tempsun = systime(1)         ; pour key_performance
+@common
+;------------------------------------------------------------
+   depthin = litchamp(tab)
+;
+   levelup = depth2level(depthin, /UPPER, /nomask)
+   depthup = level2depth(levelup, /nomask)
+;
+   levellow = depth2level(depthin, /lower, /nomask)
+   depthlow = level2depth(levellow, /nomask)
+; calcule de la distance depthlow-depthup et gestion du cas ou cette
+; distance est nulle ou egale a !values.f_nan
+   divi = depthlow-depthup
+   nan = where(finite(divi) EQ 0)
+   if nan[0] NE -1 then divi[nan] = 0
+   nan = where(divi EQ 0)
+   if nan[0] NE -1 then divi[nan] = !values.f_nan
+; calcule du resultat
+   res = levelup+(depthin-depthup)/divi
+; on masque les points terre a valmask
+   if NOT keyword_set(nomask) then begin
+      grille,mask
+      if n_elements(valmask) EQ 0 then valmask = 1e20
+      terre = where((temporary(mask))[*, *, 0] EQ 0)
+      if terre[0] NE -1 then res[terre] = valmask
+   endif
+;------------------------------------------------------------
+   if keyword_set(key_performance) THEN print, 'temps depth2floatlevel', systime(1)-tempsun 
+;
+   return, res
+end
Index: trunk/ToBeReviewed/CALCULS/depth2level.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/depth2level.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/depth2level.pro	(revision 25)
@@ -0,0 +1,120 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME: depth2level
+;
+; PURPOSE: permet de passer d''un tableau 2d de profondeur au tableau
+; 2d correspondant de niveaux.
+;
+; CATEGORY: SANS BOUCLE
+;
+; CALLING SEQUENCE: res=depth2level(depth2d)
+;
+; INPUTS: depth2d tableau 2d de profondeur (ou une structure repondant
+; aux criteres de litchamp)
+;
+; KEYWORD PARAMETERS:
+;
+;      /UPPER: (active par defaut) on selectionne le niveau
+;      directement au dessus de la profondeur
+;
+;      /LOWER: on selectionne le niveau directement au dessous de la
+;      profondeur 
+;
+;      /CLOSER: on selectionne le niveau le plus proche de la
+;      profondeur 
+;
+;      /NOMASK: pour ne pas masquer les points terres
+;
+; OUTPUTS: un tableau 2d contenant les valeurs des niveaux.
+;
+; COMMON BLOCKS:common.pro
+;
+; SIDE EFFECTS:pour les profondeurs hors des valeurs de gdep, la
+; valeur !values.f_nan est retournee.
+; Si la profondeur est superieur a celle du fond, on retourne
+; jpk-1dans le cas upper, et !values.f_nan ds le cas lower.
+;
+; RESTRICTIONS:
+;
+; EXAMPLE:
+;
+; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
+;                       17/6/1999
+;                       15/6/2000 accepte !values.f_nan
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION depth2level, tab, LOWER = lower, UPPER = upper, CLOSER = closer $
+                      , NOMASK = nomask, _extra = ex
+   tempsun = systime(1)         ; pour key_performance
+@common
+;------------------------------------------------------------
+   upper = 1
+   if keyword_set(lower) THEN upper = 0
+;------------------------------------------------------------
+; lecture du champ d''entree et recuperation de la taille du sous
+; domaine utilise
+;------------------------------------------------------------
+   in = litchamp(tab)
+   grille,mask,-1,-1,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz
+;---------------------------------------------------------------
+; verification de la coherence entre la taille du tableau et le domaine definit par domdef
+;---------------------------------------------------------------
+   IF ny EQ 1 THEN in = reform(in, nx, ny, /over)
+   taille = size(in)
+   if taille[0] NE 2 then return, report('le champ en entree doit contenir un tableau 2d')
+   case 1 of
+      taille[1] eq jpi and taille[2] eq jpj:in=in[firstx:lastx, firsty:lasty]
+      taille[1] eq  nx and taille[2] eq  ny:
+      else:return, report('Probleme d''adequation entre les tailles du domaine et celle du champ.')
+   endcase
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+; vire les points a !values.f_nan
+   notanumber = where(finite(in, /nan) EQ 1)
+   if notanumber[0] NE -1 then in[notanumber] = 0
+;------------------------------------------------------------
+; on transforme le tableau 2d de profondeur en tableau 2d de niveaux
+; correspondant aux profondeurs
+;------------------------------------------------------------
+; on passe en tableaux qui ont la taille des tableaux 3d
+   prof=replicate(1,nx*ny)#gdep[firstz:lastz] 
+   in = in[*]#replicate(1, nz)
+;
+   mask01 = (prof[*] LT in[*])
+   mask01 = reform(mask01, nx, ny, nz)
+   levels = total(mask01, 3)
+   notvalid = where(levels EQ nz)
+   if keyword_set(upper) then begin
+      levels = levels-1
+      notvalid = where(levels EQ -1)
+   ENDIF ELSE notvalid = where(levels EQ nz)
+   IF notvalid[0] NE -1 THEN levels[notvalid] = !values.f_nan
+
+; si closer est active
+   if keyword_set(closer) then begin
+      test = [ [[litchamp(tab)-level2depth(levels)]] $
+               , [[level2depth( (levels+1)<(jpk-1) )-litchamp(tab)]] ]
+      test = test[*, *, 0]-test[*, *, 1]
+      changer = where(test GE 0)
+      if changer[0] NE -1 then levels[changer] = (levels[changer]+1) < (jpk-1)
+   endif
+;------------------------------------------------------------
+; on replace les points a !values.f_nan
+   if notanumber[0] NE -1 then levels[notanumber] = !values.f_nan
+; on masque les points terres a valmask
+   if NOT keyword_set(nomask) then begin
+      if n_elements(valmask) EQ 0 then valmask = 1e20
+      terre = where(mask[*, *, 0] EQ 0)
+      if terre[0] NE -1 then levels[terre] = valmask
+   endif
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+   if keyword_set(key_performance) THEN print, 'temps depth2level', systime(1)-tempsun 
+   return, levels
+end
Index: trunk/ToBeReviewed/CALCULS/depth2mask.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/depth2mask.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/depth2mask.pro	(revision 25)
@@ -0,0 +1,64 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME: depth2mask
+;
+; PURPOSE: permet de passer d''un tableau 2d de profondeur seuil au
+; tableau 3d de mask avec des 1 ds les niveaux au dessus de la
+; profondeur seuil et des 0 en dessous.
+;
+; CATEGORY: SANS BOUCLE
+;
+; CALLING SEQUENCE: res=depth2mask(depht2d)
+;
+; INPUTS: depht2d tableau 2d de profondeur seuil (ou une structure repondant
+; aux criteres de litchamp)
+;
+; KEYWORD PARAMETERS:
+;
+;      /UPPER: (active par defaut) on selectionne le niveau
+;      directement au dessus de la profondeur
+;
+;      /LOWER: on selectionne le niveau directement au dessous de la
+;      profondeur 
+;
+;      /CLOSER: on selectionne le niveau le plus proche de la
+;      profondeur 
+;
+; OUTPUTS: un tableau 3d contenant le mask associe au tableau 2d de
+; profondeurs seuil
+;
+; COMMON BLOCKS:
+;	common.pro
+;
+; SIDE EFFECTS:
+;
+; RESTRICTIONS:
+;
+; EXAMPLE:
+;
+; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
+;                       17/6/1999
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION depth2mask, tab, _extra = ex
+;------------------------------------------------------------
+   tempsun = systime(1)         ; pour key_performance
+@common
+;------------------------------------------------------------
+; on transforme le tableau 2d de profondeur en tableau 2d de niveaux
+; correspondant aux profondeurs
+;------------------------------------------------------------
+   niveaux = depth2level(tab, _extra = ex)
+   IF niveaux[0] EQ -1 THEN return,  -1
+;------------------------------------------------------------
+; on transforme le tableau 2d de niveaux en tableau 3d de mask
+;------------------------------------------------------------
+   mask = level2mask(niveaux)
+;------------------------------------------------------------
+   if keyword_set(key_performance) NE 0 THEN print, 'temps depth2mask', systime(1)-tempsun 
+   return, mask
+end
Index: trunk/ToBeReviewed/CALCULS/determ2.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/determ2.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/determ2.pro	(revision 25)
@@ -0,0 +1,42 @@
+;+
+; NAME:determ2
+;
+; PURPOSE: computes the determinant of n 2 by 2 arrays
+;
+; CATEGORY: no DO loops and better accuracy
+;
+; CALLING SEQUENCE: 2 cases:
+;     res = determ2(z2ds)
+;     res = determ2(z1d00,z1d01,z1d10,z1d11)
+;
+; INPUTS:
+;     z2ds: an 2*2*n array
+;     or
+;     z1d00,z1d01,z1d10,z1d11: the four n elements arrays
+;     defined as:
+;              z2ds[0, 0, *] = z1d00
+;              z2ds[0, 1, *] = z1d01
+;              z2ds[1, 0, *] = z1d10
+;              z2ds[1, 1, *] = z1d11
+;
+; OUTPUTS: n elements array, the determinent of each 2*2 arrrays
+;
+; EXAMPLE:
+;
+; a=findgen(2,2,5)
+; print, determ2(a)
+; FOR i=0,4 DO print, determ(a[*,*,i]) ; IDL solution
+;
+; MODIFICATION HISTORY:
+;           S. Masson (smasson@lodyc.jussieu.fr)
+;           July 11th, 2002
+;-
+FUNCTION determ2, a, b, c, d
+  CASE n_params() OF
+    1:res = a[0, 0, *]*a[1, 1, *]-a[0, 1, *]*a[1, 0, *]
+    4:res = a[*]*d[*]-c[*]*b[*]
+    ELSE:stop
+  ENDCASE
+
+  RETURN, res
+END
Index: trunk/ToBeReviewed/CALCULS/determ3.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/determ3.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/determ3.pro	(revision 25)
@@ -0,0 +1,67 @@
+;+
+; NAME:determ3
+;
+; PURPOSE: computes the determinant of n 3 by 3 arrays
+;
+; CATEGORY: no DO loops and better accuracy
+;
+; CALLING SEQUENCE:2 cases:
+;     res = determ2(z2ds)
+;     res = determ2(in00, in01, in02, in10, in11, in12, in20, in21, in22)
+;
+; INPUTS:
+;     z2ds: an 3*3*n array
+;     or
+;     in00, in01, in02, in10, in11, in12, in20, in21, in22: the nine n
+;     elements arrays defined as:
+;           in00 = z2ds[0, 0, *]
+;           in01 = z2ds[0, 1, *]
+;           in02 = z2ds[0, 2, *]
+;           in10 = z2ds[1, 0, *]
+;           in11 = z2ds[1, 1, *]
+;           in12 = z2ds[1, 2, *]
+;           in20 = z2ds[2, 0, *]
+;           in21 = z2ds[2, 1, *]
+;           in22 = z2ds[2, 2, *]
+;
+; OUTPUTS: n elements array, the determinent of each 3*3 arrrays
+;
+; EXAMPLE:
+;
+; a=findgen(3,3,5)
+; print, determ3(a^2)
+; FOR i=0,4 DO print, determ((a[*,*,i])^2) ; IDL solution
+;
+; MODIFICATION HISTORY:
+;           S. Masson (smasson@lodyc.jussieu.fr)
+;           July 11th, 2002
+;-
+FUNCTION determ3, in00, in01, in02, in10, in11, in12, in20, in21, in22
+  IF n_params() EQ 1 THEN BEGIN
+    in00save = temporary(in00)
+    in00 = in00save[0, 0, *]
+    in01 = in00save[0, 1, *]
+    in02 = in00save[0, 2, *]
+    in10 = in00save[1, 0, *]
+    in11 = in00save[1, 1, *]
+    in12 = in00save[1, 2, *]
+    in20 = in00save[2, 0, *]
+    in21 = in00save[2, 1, *]
+    in22 = in00save[2, 2, *]
+  ENDIF
+
+  a01 = -determ2(in10[*], in20[*] $
+                 , in12[*], in22[*])
+
+  a11 = determ2(in00[*], in20[*] $
+                , in02[*], in22[*])
+
+  a21 = -determ2(in00[*], in10[*] $
+                 , in02[*], in12[*])
+
+  res = in01[*]*a01+in11[*]*a11+in21[*]*a21
+
+  IF n_params() EQ 1 THEN in00 = temporary(in00save)
+
+  RETURN, res
+END
Index: trunk/ToBeReviewed/CALCULS/div.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/div.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/div.pro	(revision 25)
@@ -0,0 +1,316 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:div
+;
+; PURPOSE:calcule la divergence d'un champ 2D
+;
+; CATEGORY:calcule sur les matrices
+;
+; CALLING SEQUENCE:res=div(u,v)
+;
+; INPUTS:
+;       u et v deux matrices representant les coordonnes d''un
+;       champ de vecteur
+;
+; KEYWORD PARAMETERS:
+;
+; OUTPUTS:res: une matrice 2d
+;
+; COMMON BLOCKS:
+;	common.pro
+;
+; SIDE EFFECTS:
+;
+; RESTRICTIONS:
+; les matrices u et v peuvent de 2 a 4 dimensions.
+; attention pour distinger les differents configurations de u et v
+; (xy, xyz, xyt, xyzt), on regarde la variable du common 
+;        -time qui contient le calendrier en jour julien d''IDL auquel
+;        se rapportent u et v ansi que la variable 
+;        -jpt qui est le nombre de pas de temps a considerer ds time.
+; les tableaux u et v sont decoupes sur le meme domaine
+; geographique. A cause du decalage des grilles T, U, V et F il est
+; possiible que ces 2 tableaux n''aient pas la meme taille et se
+; repportent a des indices differents. Si tel est le cas les tableaux
+; sont redecoupes sur les indices qu'ils ont en commun et le dommaine
+; est redefinit pour qu'il colle a ces indices communs.
+; pour eviter ces redecoupes utiliser le mot cles /memeindice ds
+; domdef.pro
+;
+; les points sur le bord du dessin sont mis a !values.f_nan
+;
+; EXAMPLE:
+;
+; MODIFICATION HISTORY:Guillaume Roullet (grlod@ipsl.jussieu.fr)
+;                      Creation : printemps 1998
+;                      Sebastien Masson (smasson@lodyc.jussieu.fr)
+;                      adaptation pour marcher avec un domaine reduit
+;                      12/1/2000
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION div, uu, vv
+   tempsun = systime(1)         ; pour key_performance
+@common
+;
+   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
+     return, report(['This version of div is based on Arakawa C-grid.' $
+                     , 'U and V grids must therefore be defined'])
+;
+   u = litchamp(uu)
+   v = litchamp(vv)
+;
+   date1 = time[0]
+   if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
+   if (size(u))[0] NE (size(v))[0] then return,  -1
+
+;------------------------------------------------------------
+; on trouve les points que u et v ont en communs
+;------------------------------------------------------------
+   indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
+   indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
+   indicex = inter(indicexu, indicexv)
+   indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
+   indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
+   indicey = inter(indiceyu, indiceyv)
+   nx = n_elements(indicex) 
+   ny = n_elements(indicey)
+   indice2d = lindgen(jpi, jpj)
+   indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
+   case 1 of
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xyz
+;----------------------------------------------------------------------------
+      (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN 
+;------------------------------------------------------------
+; extraction de u et v sur le domaine qui convient
+;------------------------------------------------------------
+         case 1 of
+            (size(v))[0] NE 3: return,  -1
+            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
+             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
+               case 1 of
+                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
+                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
+                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
+                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
+                  ELSE :
+               endcase
+            END
+            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
+             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
+               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
+               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
+            END
+            ELSE:BEGIN 
+               zdiv = -1
+               GOTO, sortie
+            end
+            
+         endcase
+;------------------------------------------------------------
+; calcul de la divergence
+;------------------------------------------------------------
+         zu = (e2u[indice2d])[*]#replicate(1, nzt)
+         zu = reform(zu, nx, ny, nzt, /over)
+         zu = temporary(u)*temporary(zu) $
+          *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
+         terreu = where(zu EQ 0)
+         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
+;
+         zv = (e1v[indice2d])[*]#replicate(1, nzt)
+         zv = reform(zv, nx, ny, nzt, /over)
+         zv = temporary(v)*temporary(zv) $
+          *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
+         terrev = where(zv EQ 0)
+         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
+;
+         zdiv = 1e6/(e1t[indice2d]*e2t[indice2d])
+         zdiv = (zdiv)[*]#replicate(1, nzt)
+         zdiv = reform(zdiv, nx, ny, nzt, /over)
+         zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0)) $
+          *tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
+;------------------------------------------------------------
+; mise a !values.f_nan de la bordure
+;------------------------------------------------------------
+         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
+            zdiv(0, *, *) = !values.f_nan
+            zdiv(nx-1, *, *) = !values.f_nan
+         endif
+         zdiv(*, 0, *) = !values.f_nan
+         zdiv(*, ny-1, *) = !values.f_nan
+;
+         zdiv = temporary(zdiv)
+;
+         if n_elements(valmask) EQ 0 THEN valmask = 1e20
+         terre =  where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0)
+         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
+;------------------------------------------------------------
+; pour le trace graphique
+;------------------------------------------------------------
+         vargrid = 'T'
+         varname = 'div'
+         varunits = '1e6*s-1'
+         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
+         if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan)
+      END
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xyt
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+      date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN 
+;------------------------------------------------------------
+; extraction de u et v sur le domaine qui convient
+;------------------------------------------------------------
+         case 1 of
+            (size(u))[0] NE 3 OR (size(v))[0] NE 3: return,  -1
+            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
+             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
+               case 1 of
+                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
+                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
+                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
+                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
+                  ELSE :
+               endcase
+            END
+            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
+             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
+               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
+               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
+            END
+            ELSE:return,  -1
+         endcase
+;------------------------------------------------------------
+; calcul de la divergence
+;------------------------------------------------------------
+         zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
+         terreu = where(zu EQ 0)
+         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
+         zu = (zu)[*]#replicate(1, jpt)
+         zu = reform(zu, nx, ny, jpt, /over)
+         zu = temporary(u)*temporary(zu)
+;
+         zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
+         terrev = where(zv EQ 0)
+         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
+         zv = (zv)[*]#replicate(1, jpt)
+         zv = reform(zv, nx, ny, jpt, /over)
+         zv = temporary(v)*temporary(zv)
+;
+         zdiv = 1e6*tmask[indice2d+jpi*jpj*firstzt]/(e1t[indice2d]*e2t[indice2d])
+         zdiv = (zdiv)[*]#replicate(1, jpt)
+         zdiv = reform(zdiv, nx, ny, jpt, /over)
+         terre =  where(zdiv EQ 0)
+         zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0))
+;------------------------------------------------------------
+; mise a !values.f_nan de la bordure
+;------------------------------------------------------------
+         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
+            zdiv(0, *, *) = !values.f_nan
+            zdiv(nx-1, *, *) = !values.f_nan
+         endif
+         zdiv(*, 0, *) = !values.f_nan
+         zdiv(*, ny-1, *) = !values.f_nan
+;
+         if n_elements(valmask) EQ 0 THEN valmask = 1e20
+         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
+;------------------------------------------------------------
+; pour le trace graphique
+;------------------------------------------------------------
+         vargrid = 'T'
+         varname = 'div'
+         varunits = '1e6*s-1'
+         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
+         if keyword_set(direc) then  zdiv = grossemoyenne(zdiv,direc,/nan)
+      END
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xyzt
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+      date1 NE date2 AND (size(u))[0] EQ 4:BEGIN 
+         return, report('non code!')
+      END
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xy
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+      ELSE:BEGIN                ;xy
+         indice3d = lindgen(jpi, jpj, jpk)
+         indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt]
+;------------------------------------------------------------
+; extraction de u et v sur le domaine qui convient
+;------------------------------------------------------------
+         case 1 of
+            (size(u))[0] NE 2 OR (size(v))[0] NE 2: BEGIN 
+               zdiv = -1
+               GOTO, sortie
+            end
+            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
+             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
+               case 1 of
+                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]
+                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
+                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]
+                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
+                  ELSE :
+               endcase
+            END
+            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
+             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
+               u = u[indice2d]
+               v = v[indice2d]
+            END
+            ELSE:return,  -1
+         endcase
+;------------------------------------------------------------
+; calcul de la divergence
+;------------------------------------------------------------
+         zu = temporary(u)*e2u[indice2d]*(umask())[indice3d]
+         terreu = where(zu EQ 0)
+         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
+         zv = temporary(v)*e1v[indice2d]*(vmask())[indice3d]
+         terrev = where(zv EQ 0)
+         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
+         zdiv = zu-shift(zu, 1, 0)+zv-shift(zv, 0, 1)
+         zdiv = temporary(zdiv)*tmask[indice3d]/(e1t[indice2d]*e2t[indice2d])
+;------------------------------------------------------------
+; mise a !values.f_nan de la bordure
+;------------------------------------------------------------
+         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
+            zdiv(0, *) = !values.f_nan
+            zdiv(nx-1, *) = !values.f_nan
+         endif
+         zdiv(*, 0) = !values.f_nan
+         zdiv(*, ny-1) = !values.f_nan
+;
+         zdiv = temporary(zdiv)*1e6
+;
+         if n_elements(valmask) EQ 0 THEN valmask = 1e20
+         terre =  where(tmask[indice3d] EQ 0)
+         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
+;------------------------------------------------------------
+; pour le trace graphique
+;------------------------------------------------------------
+         vargrid = 'T'
+         varname = 'div'
+         varunits = '1e6*s-1'
+         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
+         if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan)
+
+      END
+   endcase
+   
+;------------------------------------------------------------
+sortie:
+   if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun 
+   
+   return, zdiv
+end
Index: trunk/ToBeReviewed/CALCULS/floatlevel2depth.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/floatlevel2depth.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/floatlevel2depth.pro	(revision 25)
@@ -0,0 +1,69 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:floatlevel2depth
+;
+; PURPOSE: assez comparable a level2depth. C'est la fonction inverse
+; de depth2floatlevel
+;
+; CATEGORY:SANS BOUCLE
+;
+; CALLING SEQUENCE:res=floatlevel2depth(niveau)
+; 
+; INPUTS: tableau 2d de niveaux seuil (ou une structure repondant aux
+; criteres de litchamp)
+;
+; KEYWORD PARAMETERS:
+;       /NOMASK: pour ne pas masquer les points terres
+;
+; OUTPUTS: un tableau 2d contenant des profondeurs
+;
+; COMMON BLOCKS:common.pro
+;
+; SIDE EFFECTS:
+;
+; RESTRICTIONS:
+;
+; EXAMPLE:
+;
+;    IDL> a=gdept[jpk-1]/(1.*jpi*jpj)*findgen(jpi,jpj)
+;    IDL> plt, 1e6*(a-floatlevel2depth(depth2floatlevel(a))),/nocontour
+;
+;  ->champ nul a 1e-6 pres
+;
+; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
+;                      15/6/2000
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION floatlevel2depth, tab, NOMASK = nomask
+   tempsun = systime(1)         ; pour key_performance
+@common
+;------------------------------------------------------------
+   flevelin = litchamp(tab)
+; on vire les points a !values.f_nan
+   notanumber = where(finite(flevelin, /nan) EQ 1)
+   if notanumber[0] NE -1 then flevelin[notanumber] = 0
+; on seuil (vire les points terres a valmask par ex)
+   flevelin = 0 > flevelin < (jpk-1)
+; on calcule la profondeur 
+   depthup = level2depth(floor(flevelin), /nomask)
+   depthlow= level2depth(ceil(flevelin), /nomask)
+   weight = flevelin-floor(flevelin)
+   res = depthup+weight*(depthlow-depthup)
+; on replace les points a !values.f_nan
+   if notanumber[0] NE -1 then res[notanumber] = !values.f_nan
+; on masque les points terres a valmask
+   if NOT keyword_set(nomask) then begin
+      grille,mask
+      if n_elements(valmask) EQ 0 then valmask = 1e20
+      terre = where((temporary(mask))[*, *, 0] EQ 0)
+      if terre[0] NE -1 then res[terre] = valmask
+   endif
+;------------------------------------------------------------
+   if keyword_set(key_performance) THEN print, 'temps floatlevel2depth', systime(1)-tempsun 
+;
+   return, res
+end
Index: trunk/ToBeReviewed/CALCULS/fsfzpt.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/fsfzpt.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/fsfzpt.pro	(revision 25)
@@ -0,0 +1,17 @@
+FUNCTION fsfzpt, pfs, pfp
+
+;+
+; Ice freezing point
+; ==================
+;  fsfzpt: freezing point of seawater in degrees celsius
+;       units : salinity        pfs       (ipss-78)
+;               pressure        pfp      decibars
+;               temperature     fszfpt   degrees celsius
+;               freezing pt
+;	reference : unesco tech. papers in the marine science no 28 1978
+;               eigth report jpots
+;               annex 6 freezing point of seawater F.J.Millero pp.29-35
+;	checkvalue: fsfzpt=-2.588567 deg.c,for s=40.0,p=500 decibars
+;-
+  RETURN, ( -0.0575 + 1.710523e-3 * sqrt(pfs) - 2.154996e-4 * pfs  ) * pfs - 7.53e-4 * pfp
+END
Index: trunk/ToBeReviewed/CALCULS/grad.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/grad.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/grad.pro	(revision 25)
@@ -0,0 +1,273 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:
+;
+; PURPOSE:
+;
+; CATEGORY:
+;
+; CALLING SEQUENCE:
+; 
+; INPUTS:
+;
+; KEYWORD PARAMETERS:
+;
+; OUTPUTS:
+;
+; COMMON BLOCKS:common.pro
+;
+; SIDE EFFECTS:
+;
+; RESTRICTIONS:
+;
+; EXAMPLE:
+;
+; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
+;
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION grad, field, direc
+@common
+;------------------------------------------------------------
+;
+   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
+     return, report(['This version of grad is based on Arakawa C-grid.' $
+                     , 'U and V grids must therefore be defined'])
+;
+   res = litchamp(field)
+   taille=size(res)
+   grille, mask, glam, gphi, gdep, nx, ny,nz,firstx,firsty,firstz,lastx, lasty, lastz     
+   if n_elements(valmask) EQ 0 then valmask = 1e20
+   case strupcase(vargrid) of
+      'T':BEGIN
+         case direc of
+            'x':BEGIN 
+               divi = e1u[firstx:lastx, firsty:lasty]
+               newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
+               vargrid = 'U'
+               domdef, glamt[firstx, 0], glamu[lastx, 0] $
+                , gphit[0, firsty], gphiu[0, lasty], gridtype = ['T','U']
+            END
+            'y':BEGIN
+               divi = e2v[firstx:lastx, firsty:lasty]
+               newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
+               vargrid = 'V'
+               domdef, glamt[firstx, 0], glamv[lastx, 0] $
+                , gphit[0, firsty], gphiv[0, lasty], gridtype = ['T','V']
+            END
+            'z':BEGIN
+               divi = e3w[firstz:lastz]
+               newmask = mask
+               vargrid = 'W'
+            END
+            ELSE:return, report('Bad definition of direction argument')
+         ENDCASE
+      END
+      'W':BEGIN
+         case direc of
+            'x':divi = e1u[firstx:lastx, firsty:lasty]
+            'y':divi = e2v[firstx:lastx, firsty:lasty]
+            'z':BEGIN
+               divi = e3t[firstz:lastz]
+               newmask = mask
+               vargrid = 'T'
+            END
+            ELSE:return, report('Bad definition of direction argument')
+         endcase
+      END
+      'U':BEGIN
+         case direc of
+            'x':BEGIN
+               divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty]
+               newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
+               vargrid = 'T'
+               domdef, glamt[firstx, 0], glamu[lastx] $
+                , gphit[0, firsty], gphiu[0, lasty], gridtype = ['T','U']
+            END
+            'y':BEGIN
+               divi = e2f[firstx:lastx, firsty:lasty]
+               newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
+               vargrid = 'F'
+               domdef, glamu[firstx, 0], glamf[lastx, 0] $
+                , gphiu[0, firsty], gphif[0, lasty], gridtype = ['U','F']
+            END
+            'z':BEGIN
+               divi = e3w[firstz:lastz]
+               newmask = mask
+               vargrid = 'W'
+            END
+            ELSE:return, report('Bad definition of direction argument')
+         endcase
+      END
+      'V':BEGIN
+         case direc of
+            'x':BEGIN
+               divi = e1f[firstx:lastx, firsty:lasty]
+               newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
+               vargrid = 'F'
+               domdef, glamv[firstx, 0], glamf[lastx, 0] $
+                , gphiv[0, firsty], gphif[0, lasty], gridtype = ['V','F']
+            END
+            'y':BEGIN
+               divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty]
+               newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
+               vargrid = 'T'
+               domdef, glamt[firstx, 0], glamv[lastx, 0] $
+                , gphit[0, firsty], gphiv[0, lasty], gridtype = ['T','V']
+            END
+            'z':BEGIN
+               divi = e3w[firstz:lastz]
+               newmask = mask
+               vargrid = 'W'
+            END
+            ELSE:return, report('Bad definition of direction argument')
+         endcase
+      END
+;       'F':BEGIN
+;          case direc of
+;             'x':divi = (shift(e1v, -1, 0))[firstx:lastx, firsty:lasty]
+;             'y':divi = (shift(e2u, 0, -1))[firstx:lastx, firsty:lasty]
+;             'z':divi = e3w[firstz:lastz]
+;             ELSE:return, report('Bad definition of direction argument')
+;          endcase
+;       END
+      ELSE:return, report('Bad definition of vargrid')
+   ENDCASE
+   res = fitintobox(res)
+   case 1 of
+;----------------------------------------------------------------------------
+;xy
+;----------------------------------------------------------------------------
+      taille[0] EQ 2:BEGIN
+         earth = where(mask[*, *, firstz] EQ 0)
+         if earth[0] NE -1 then res[earth] = !values.f_nan
+         case direc of
+            'x':BEGIN 
+               res = (shift(res, -1, 0)-res)/divi
+               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan
+               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0)
+            END
+            'y':BEGIN 
+               res = (shift(res, 0, -1)-res)/divi
+               res[*, ny-1] = !values.f_nan
+               if vargrid EQ 'T' OR vargrid EQ 'U' then res =   shift(res, 0, 1)           
+            END
+            ELSE:return,  report('Bad definition of direction argument for the type of array')
+         ENDCASE
+         earth = where(newmask[*, *, firstz] EQ 0)
+         if earth[0] NE -1 then res[earth] = valmask
+      END
+;----------------------------------------------------------------------------
+;xyt
+;----------------------------------------------------------------------------
+      taille[0] EQ 3 AND jpt NE 1:BEGIN
+         earth = where(mask[*, *, firstz] EQ 0)
+         if earth[0] NE -1 then BEGIN
+            earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt))
+            res[earth[*]] = !values.f_nan
+         ENDIF
+         divi = divi[*]#replicate(1, jpt)
+         case direc of
+            'x':BEGIN 
+               res = (shift(res, -1, 0, 0)-res)/divi
+               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
+               if vargrid EQ 'T' OR vargrid EQ 'V'  then res = shift(res, 1, 0, 0)
+            END
+            'y':BEGIN 
+               res = (shift(res, 0, -1, 0)-res)/divi
+               res[*, ny-1, *] = !values.f_nan
+               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0)        
+            END
+            ELSE:return,  report('Bad definition of direction argument for the type of array')
+         ENDCASE
+         earth = where(newmask[*, *, firstz] EQ 0)
+         if earth[0] NE -1 THEN res[earth] = valmask
+      END
+;----------------------------------------------------------------------------
+;xyz
+;----------------------------------------------------------------------------
+      taille[0] EQ 3 AND jpt EQ 1:BEGIN
+         earth = where(mask EQ 0)
+         if earth[0] NE -1 then res[earth] = !values.f_nan
+         case direc OF
+            'x':BEGIN 
+               divi = divi[*]#replicate(1, nz)
+               res = (shift(res, -1, 0, 0)-res)/divi
+               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
+               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0)
+            END
+            'y':BEGIN 
+               divi = divi[*]#replicate(1, nz)
+               res = (shift(res, 0, -1, 0)-res)/divi
+               res[*, ny-1, *] = !values.f_nan
+               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0)        
+            END
+            'z':BEGIN
+               divi = reform(replicate(1, nx*ny)#divi, nx, ny, nz)
+               if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz)
+               if vargrid EQ 'W' THEN BEGIN 
+                  res = (shift(res, 0, 0, 1)-res)/divi
+                  res[*, *, 0] = !values.f_nan
+               ENDIF ELSE BEGIN
+                  res = (res-shift(res, 0, 0, -1))/divi
+                  res[*, *, nz-1] = !values.f_nan
+               ENDELSE
+               if earth[0] NE -1 then res[earth] = valmask
+            END
+         ENDCASE
+      END
+;------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xyzt
+;----------------------------------------------------------------------------
+      taille[0] EQ 4:BEGIN
+         earth = where((mask[*]#replicate(1, jpt)) EQ 0)
+         if earth[0] NE -1 then res[earth] = !values.f_nan
+         case direc OF
+            'x':BEGIN 
+               divi = divi[*]#replicate(1, nz*jpt)
+               res = (shift(res, -1, 0, 0, 0)-res)/divi
+               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *, *] = !values.f_nan
+               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0, 0)
+            END
+            'y':BEGIN 
+               divi = divi[*]#replicate(1, nz*jpt)
+               res = (shift(res, 0, -1, 0, 0)-res)/divi
+               res[*, ny-1, *, *] = !values.f_nan
+               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0, 0)        
+            END
+            'z':BEGIN
+               divi = replicate(1, nx*ny)#divi
+               divi = reform(divi[*]#replicate(1, jpt), nx, ny, nz, jpt, /over)
+               if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt)
+               if vargrid EQ 'W' THEN BEGIN 
+                  res = (shift(res, 0, 0, 1, 0)-res)/divi
+                  res[*, *, 0, *] = !values.f_nan
+               ENDIF ELSE BEGIN
+                  res = (res-shift(res, 0, 0, -1, 0))/divi
+                  res[*, *, nz-1, *] = !values.f_nan
+               ENDELSE
+            END
+         ENDCASE
+         if earth[0] NE -1 then res[earth] = valmask
+      END
+;------------------------------------------------------------
+;------------------------------------------------------------
+   endcase
+   varname = 'grad of '+varname
+   varunit = varunit+'/m'
+
+
+;------------------------------------------------------------
+   return, res
+end
+
+
+
+
+
+
Index: trunk/ToBeReviewed/CALCULS/grossemoyenne.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/grossemoyenne.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/grossemoyenne.pro	(revision 25)
@@ -0,0 +1,638 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; 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',BOXZOOM=boxzoom)
+;
+; 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:
+;
+;		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 grossemoyenne est appelee via
+;               checkfield)
+;
+;               INTEGRATION: pour faire une integrale plutot qu''une
+;               moyenne
+;
+;               /SPATIALFIRST, when performing at the same time
+;               spatial and temporal mean, grossemoyenne is assuming
+;               that the mask is not changing with the time. In
+;               consequence, grossemoyenne performs temporal mean
+;               first and then call moyenne. Activate /SPATIALFIRST if
+;               you want to perform the spatial mean before the
+;               temporal mean. Note that if NAN is activated, then
+;               SPATIALFIRST is activated automatically.
+;
+;               /TEMPORALFIRST: to force to perform first temporal
+;               mean even if nan is activated (see SPATIALFIRST
+;               explanations...)
+;
+; 
+;         /WDEPTH: to specify that the field is at W depth instad of T 
+;         depth (automatically activated if vargrid eq 'W')
+;            
+; 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, BOXZOOM = boxzoom, INTEGRATION = integration $
+                        , NAN = nan, NODOMDEF = nodomdef, WDEPTH = wdepth $
+                        , SPATIALFIRST = spatialfirst, TEMPORALFIRST = temporalfirst $
+                        , _extra = ex
+;---------------------------------------------------------
+@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
+  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
+  IF keyword_set(NAN) AND (dirx EQ 1 OR diry EQ 1 OR dirz EQ 1) $
+    THEN spatialfirst = 1
+  IF keyword_set(temporalfirst) THEN spatialfirst = 0
+;------------------------------------------------------------
+; 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 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('Wrong Definition of Boxzoom')
+    endcase
+    if NOT keyword_set(nodomdef) then BEGIN 
+      savedbox = 1b
+      saveboxparam, 'boxparam4grmoyenne.dat'
+      domdef, bte, GRIDTYPE = vargrid, _extra = ex
+    ENDIF 
+  ENDIF 
+;---------------------------------------------------------------
+; attribution du mask et des tableaux de longitude et latitude...
+;---------------------------------------------------------------
+  grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, WDEPTH = wdepth
+;------------------------------------------------------------
+; I.3) si dirt eq 1 on fait la moyenne temporelle et on envoie ds moyenne
+;------------------------------------------------------------
+  if dirt EQ 1 AND NOT keyword_set(spatialfirst) 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[firstx:firstx+nx-1 $
+                    , firsty:firsty+ny-1, *]
+        taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpt:res = tab
+        else:BEGIN 
+          if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
+          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))
+        END
+      ENDCASE
+      if keyword_set(integration) then begin
+        res = total(res, 3, nan = nan)
+      ENDIF ELSE BEGIN
+        if keyword_set(nan) then BEGIN
+          divi = finite(res)
+          divi = total(temporary(divi), 3)
+          notanum = where(divi EQ 0)
+          res = total(res, 3, nan = keyword_set(nan))/ (1 > divi)
+          if notanum[0] NE -1 then res[temporary(notanum)] = !values.f_nan
+        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[firstx:lastx, firsty:lasty, firstz:lastz, *]
+        taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $
+          res = tab[firstx:lastx, firsty:lasty, *, *]
+        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[*, *, firstz:lastz, *]
+        else:BEGIN 
+           if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
+          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))
+        END
+        endcase
+      if keyword_set(integration) then begin
+        res = total(res, 4, nan = nan)
+      ENDIF ELSE BEGIN
+        if keyword_set(nan) then begin
+          divi = finite(res)
+          divi = total(temporary(divi), 4)
+          notanum = where(divi EQ 0)
+          res = total(res, 4, /nan)/(1 > divi)
+          if notanum[0] NE -1 then res[temporary(notanum)] = !values.f_nan
+        ENDIF ELSE res = total(res, 4)/(1.*taille[4])
+      ENDELSE
+    ENDELSE
+    if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'    
+    return,  moyenne(temporary(res), direc, BOXZOOM = boxzoom, NAN = nan, INTEGRATION = integration, NODOMDEF = nodomdef, WDEPTH = wdepth, _extra = ex)
+  ENDIF ELSE res = tab
+  IF jpt EQ 1 THEN BEGIN
+    if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'    
+    return, moyenne(reform(res, /over), direc, BOXZOOM = boxzoom, NAN = nan, INTEGRATION = integration, NODOMDEF = nodomdef, WDEPTH = wdepth, _extra = ex)
+  END
+;------------------------------------------------------------
+;------------------------------------------------------------
+; 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[firstx:firstx+nx-1 $
+                  , firsty:firsty+ny-1, *]
+      taille[1] EQ  nx and taille[2] eq  ny and taille[3] eq jpt:res = tab
+      else:BEGIN 
+        if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
+        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))
+      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, jpt, /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
+;---------------------------------------------------------------
+    if keyword_set(nan) NE 0 then msknan = finite(res) ELSE msknan = -1
+    mask = mask[*, *, 0]
+    case 1 of
+      (dirx eq 1) and (diry eq 0) : begin
+        e = temporary(e1)*temporary(mask)
+        echelle = (temporary(e))[*]#replicate(1, jpt)
+        echelle = reform(echelle, nx, ny, jpt, /over)
+        if keyword_set(integration) then divi = 1 ELSE BEGIN 
+          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 1) $
+          ELSE divi = total(echelle, 1)
+        ENDELSE
+        res = total(temporary(res)*echelle, 1, nan = nan)/(divi > 1.)
+        if msknan[0] NE -1 then BEGIN
+          echelle = temporary(echelle) NE 0
+          testnan = temporary(msknan)*echelle
+          testnan = total(temporary(testnan), 1) $
+            +(total(temporary(echelle), 1) EQ 0)
+        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 BEGIN 
+          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 2) $
+          ELSE divi = total(echelle, 2)
+        ENDELSE 
+        res = total(temporary(res)*echelle, 2, nan = nan)/(divi > 1.)
+        if msknan[0] NE -1 then begin
+          echelle = temporary(echelle) NE 0
+          testnan = temporary(msknan)*echelle
+          testnan = total(temporary(testnan), 2) $
+            +(total(temporary(echelle), 2) EQ 0)
+        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 BEGIN 
+          IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 1) $
+          ELSE divi = total(total(echelle, 1), 1)
+        ENDELSE
+        res = total(temporary(total(temporary(res)*echelle, 1, nan = nan)), 1, nan = nan)/(divi > 1.)
+        if msknan[0] NE -1 then begin
+          echelle = temporary(echelle) NE 0
+          testnan = temporary(msknan)*echelle
+          testnan = total(total(temporary(testnan), 1), 1) $
+            +(total(total(temporary(echelle), 1), 1) EQ 0)
+        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[firstx:lastx, firsty:lasty, firstz:lastz, *]
+      taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $
+        res = tab[firstx:lastx, firsty:lasty, *, *]
+      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[*, *, firstz:lastz, *]
+      else:BEGIN 
+        if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
+        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))
+      END
+    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 then BEGIN 
+      res = reform(res, nx, ny, nz, jpt, /over)
+      mask =  reform(mask, nx, ny, nz, /over)
+    ENDIF
+    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 = (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(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, jpt)
+            nan = 1
+          ENDIF
+          bottom = bottom#replicate(1, jpt) $ ; 4D bottom!
+            + replicate(1, n_elements(bottom))#(nx*ny*nz*lindgen(jpt))
+          msknan[bottom] = 0
+          res[temporary(bottom)] = !values.f_nan
+        ENDIF
+        if keyword_set(integration) then divi = 1 ELSE begin 
+          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 1) $
+          ELSE divi = total(echelle, 1)
+        endelse
+        res = temporary(res)*echelle
+        res = total(temporary(res), 1, nan = nan)/(divi > 1)
+        if msknan[0] NE -1 then begin
+          echelle = temporary(echelle) NE 0
+          testnan = temporary(msknan)*echelle
+          testnan = total(temporary(testnan), 1) $
+            +(total(temporary(echelle), 1) EQ 0)
+        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(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
+          bottom = bottom#replicate(1, jpt) $ ; 4D bottom!
+            + replicate(1, n_elements(bottom))#(nx*ny*nz*lindgen(jpt))
+          msknan[bottom] = 0
+          res[temporary(bottom)] = !values.f_nan
+        ENDIF
+        if keyword_set(integration) then divi = 1 ELSE begin 
+          IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 2) $
+          ELSE divi = total(echelle, 2)
+        endelse
+        res = total(temporary(res)*echelle, 2, nan = nan)/(divi > 1)
+        if msknan[0] NE -1 then begin
+          echelle = temporary(echelle) NE 0
+          testnan = temporary(msknan)*echelle
+          testnan = total(temporary(testnan), 2) $
+            +(total(temporary(echelle), 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
+        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 msknan[0] NE -1 THEN divi = total(echelle*msknan, 3) $
+          ELSE divi = total(echelle, 3)
+        endelse
+        res = total(temporary(res)*echelle, 3, nan = nan)/(divi > 1)
+        if msknan[0] NE -1 then begin
+          echelle = temporary(echelle) NE 0
+          testnan = temporary(msknan)*echelle
+          testnan = total(temporary(testnan), 3) $
+            +(total(temporary(echelle), 3) EQ 0)
+        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(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
+          bottom = bottom#replicate(1, jpt) $ ; 4D bottom!
+            + replicate(1, n_elements(bottom))#(nx*ny*nz*lindgen(jpt))
+          msknan[bottom] = 0
+          res[temporary(bottom)] = !values.f_nan
+        ENDIF
+        if keyword_set(integration) then divi = 1 ELSE begin 
+          IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 1) $
+          ELSE divi = total(total(echelle, 1), 1)
+        endelse
+        res = total(total(temporary(res)*echelle, 1, nan = nan), 1, nan = nan)/(divi > 1)
+        if msknan[0] NE -1 then begin
+          echelle = temporary(echelle) NE 0
+          testnan = temporary(msknan)*echelle
+          testnan = total(total(temporary(testnan), 1), 1) $
+            +(total(total(temporary(echelle), 1), 1) EQ 0)
+        endif
+      end
+      (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin
+        e133 = e1[*]#e3
+        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
+        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 msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 2) $
+          ELSE divi = total(total(echelle, 1), 2)
+        endelse
+        res = total(total(temporary(res)*echelle, 1, nan = nan), 2, nan = nan)/(divi > 1)
+        if msknan[0] NE -1 then begin
+          echelle = temporary(echelle) NE 0
+          testnan = temporary(msknan)*echelle
+          testnan = total(total(temporary(testnan), 1), 2) $
+            +(total(total(temporary(echelle), 1), 2) EQ 0)
+        endif
+      end
+      (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin
+        e233 = e2[*]#e3
+        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
+        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 msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 2), 2) $
+          ELSE divi = total(total(echelle, 2), 2)
+        endelse
+        res = total(total(temporary(res)*echelle, 2, nan = nan), 2, nan = nan)/(divi > 1)
+        if msknan[0] NE -1 then begin
+          echelle = temporary(echelle) NE 0
+          testnan = temporary(msknan)*echelle
+          testnan = total(total(temporary(testnan), 2), 2) $
+            +(total(total(temporary(echelle), 2), 2) EQ 0)
+        endif
+      end
+      (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin
+        e1233 = (e1[*]*e2[*])#e3
+        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
+        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 msknan[0] NE -1 THEN divi = total(total(total(echelle*msknan, 1), 1), 1) $
+          ELSE divi = total(total(total(echelle, 1), 1), 1)
+        endelse
+        res = total(total(total(temporary(res)*echelle, 1, nan = nan), 1, nan = nan), 1, nan = nan)/(divi > 1)
+        if msknan[0] NE -1 then begin
+          echelle = temporary(echelle) NE 0
+          testnan = temporary(msknan)*echelle
+          testnan = total(total(total(temporary(testnan), 1), 1), 1) $
+            +(total(total(total(temporary(echelle), 1), 1), 1) EQ 0)
+        endif
+      end
+    endcase
+  endif
+;------------------------------------------------------------
+  if dirt EQ 1 AND keyword_set(spatialfirst) then BEGIN
+    IF (reverse(size(res, /dimension)))[0] NE jpt THEN BEGIN
+      print, 'the last dimension of res is not equal to jpt: '+strtrim(jpt, 2)
+      if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
+      return, -1
+    ENDIF
+    tdim = size(res, /n_dimensions)
+    if keyword_set(integration) then res = total(res, tdim, nan = nan) $
+    ELSE BEGIN 
+      if keyword_set(nan) then BEGIN 
+        testnan = testnan < 1
+        testnan = total(temporary(testnan), tdim)
+        divi = testnan
+      ENDIF ELSE divi = jpt
+      res = total(res, tdim, nan = nan)/(1 > divi)
+    ENDELSE 
+  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(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat'
+;------------------------------------------------------------
+  if keyword_set(key_performance) THEN print, 'temps grossemoyenne', systime(1)-tempsun 
+  return, res
+;------------------------------------------------------------
+;------------------------------------------------------------
+end
Index: trunk/ToBeReviewed/CALCULS/hdyn.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/hdyn.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/hdyn.pro	(revision 25)
@@ -0,0 +1,164 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:hdyn
+;
+; PURPOSE:calcule la hauteur dynamique par rapport a un etat de
+; reference pour une profondeur de reference. Cf les mots cles pour
+; les differentes possibilites. 
+; Par defaut l''etat de reference est
+; rho=1020 et la profondeur de reference est gdepw[ka] avec ka le
+; premier niveau W directement au dessus de 1000 m. 
+;
+; CATEGORY: calculs de post-traitement
+;
+; CALLING SEQUENCE:res=hdyn(sn,tn)
+;
+; INPUTS:sn et tn sont des tableaux de meme taille representant la
+; salinite et la temperature.
+;
+; KEYWORD PARAMETERS:
+;
+;        GILL: activer cette cle si on veut faire le calcul de la
+;        hauteur dynamique comme ds le GILL page 215 cad par rapport a
+;        un etat de reference qui varie en profondeur et qui est
+;        determine par une temperature de reference tref a 0 degre et
+;        une salinite de reference sref a 35psu
+;
+;        LEVEL: C''est le niveau de reference a prendre. Ce niveau est
+;        definit tel que gdepw[level] est la profondeur de reference
+;
+;        SREF: donner une valeur a ce mot cle pour changer la salinite
+;        de reference utiliser ds le calcul lorsque GILL est active.
+;   
+;        TREF: donner une valeur a ce mot cle pour changer la temperature
+;        de reference utiliser ds le calcul lorsque GILL est active.
+;
+;        PROFREF: donner a ce mot cle une profondeur qui sera prise comme
+;        la profondeur de reference (ds ce cas LEVEL n'' a aucun
+;        effet). le calcul sera alors effectue jusqu''a cette
+;        profondeur en effectuant une interpolation entre le dernier
+;        niveau W au dessus de PROFREF et PROFREF.
+;
+;        SURFACE_LEVEL: C''est le niveau auquel on veut calculer la
+;        hauteur dynamique. Par defaut c''est le niveau 0.
+;
+; OUTPUTS:un tableau de la meme taille que sn et tn representant la
+; hauteur dynamique calculee a partir d''une profondeur de reference
+; et par rapport a un etat de reference.
+;
+; COMMON BLOCKS:
+;	common.pro
+;
+; SIDE EFFECTS:
+; les points pour lesquels on nje peut calcule la hauteur dynamique
+; (dont la batymetrie est moins profonde que la profondeur de
+; reference) sont mis a la valeur !values.f_nan
+;
+; RESTRICTIONS:approximation: la pression en decibars est egale a la
+; profondeur en metres (la pression augmente de 1bar tous les 10m)
+;
+; EXAMPLE:
+;
+; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
+;
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION hdyn,  tabsn, tabtn, TREF = tref,  SREF = sref, PROFREF = profref, LEVEL = level, GILL = gill, SURFACE_LEVEL = surface_level
+   tempsun = systime(1)         ; pour key_performance
+@common
+
+   if NOT keyword_set(surface_level) then surface_level = 0
+; utile si GILL est active
+   if NOT keyword_set(tref) then tref = 0.
+   if NOT keyword_set(sref) then sref = 35.
+; on determine si besoin est la profondeur de reference et le niveau
+; W situe directement au dessus.
+   if keyword_set(profref) then begin
+      rien = where(gdepw LE profref, level)
+      level = level-1
+      za = gdepw[level]
+   ENDIF ELSE BEGIN 
+      if NOT keyword_set(level) then BEGIN 
+         rien = where(gdepw LE 1000., level)
+         level = level-1
+      ENDIF
+      profref = gdepw[level]
+      za = profref
+   ENDELSE 
+   tailles = size(tabsn)
+   taillet = size(tabtn)
+   if total(tailles[0:tailles[0]] NE taillet[0:taillet[0]]) NE 0 then $
+    return,  report('Les tableaux sn et tn doivent avoir la meme taille')
+   if tailles[3] NE jpk then return, report('La dim verticale des tableaux sn et tn doit etre egalre a jpk')
+   nx = nxt
+   ny = nyt
+   case (size(tabsn))[0] OF
+      3:BEGIN 
+         case 1 of
+            tailles[1] eq jpi and tailles[2] eq jpj: BEGIN
+               sn = tabsn[firstxt:lastxt, firstyt:lastyt, *]
+               tn = tabtn[firstxt:lastxt, firstyt:lastyt, *]
+            end
+            tailles[1] eq  nx and tailles[2] eq  ny:BEGIN 
+               sn = tabsn
+               tn = tabtn
+            end
+            else:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
+         ENDCASE
+         if keyword_set(gill) then $
+          rhonref = rhon(replicate(sref,nx, ny, jpk),replicate(tref,nx, ny, jpk), /insitu) $
+         ELSE rhonref = 1020.
+         vol=(rhonref-rhon(sn,tn, /insitu))/rhonref
+         e33d = replicate(1, nx*ny)#e3t
+         e33d = reform(e33d, nx, ny, jpk, /over)
+         terre = where(tmask[firstxt:lastxt, firstyt:lastyt, *] EQ 0)
+         if terre[0] NE -1 then vol[terre] = !values.f_nan
+         case level of
+            0:hdyn =100.*(profref-gdepw[0])*vol[*, *, 0]
+            1:hdyn =100.*(vol*e33d)[*, *, 0]+(profref-gdepw[1])*vol[*, *, 1]
+            ELSE:hdyn =100.*total( (vol*e33d)[*, *, surface_level: level-1], 3) $
+             +(profref-gdepw[level])*vol[*, *, level]
+         endcase
+      END
+      4:BEGIN 
+         case 1 of
+            tailles[1] eq jpi and tailles[2] eq jpj AND tailles[4] EQ jpt: BEGIN
+               sn = tabsn[firstxt:lastxt, firstyt:lastyt, *, *]
+               tn = tabtn[firstxt:lastxt, firstyt:lastyt, *, *]
+            end
+            tailles[1] eq  nx and tailles[2] eq  ny AND tailles[4] EQ jpt:BEGIN 
+               sn = tabsn
+               tn = tabtn
+            end
+            else:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
+         endcase
+         if keyword_set(gill) then $
+          rhonref = rhon(replicate(sref,nx, ny, jpk, jpt),replicate(tref,nx, ny, jpk, jpt), /insitu) $
+         ELSE rhonref = 1020.
+         vol=(rhonref-rhon(sn,tn, /insitu))/rhonref
+         e33d = replicate(1, nx*ny)#e3t
+         e33d = e33d[*]#replicate(1, jpt)
+         e33d = reform(e33d, nx, ny, jpk, jpt, /over)
+         mask = tmask[firstxt:lastxt, firstyt:lastyt, *]
+         mask = mask[*]#replicate(1, jpt)
+         terre = where(mask EQ 0)
+         if terre[0] NE -1 then vol[terre] = !values.f_nan
+         case level of
+            0:hdyn =100.*(profref-gdepw[0])*vol[*, *, 0, *]
+            1:hdyn =100.*(vol*e33d)[*, *, 0, *]+(profref-gdepw[1])*vol[*, *, 1, *]
+            ELSE:hdyn =100.*total( (vol*e33d)[*, *, surface_level: level-1, *], 3) $
+             +(profref-gdepw[level])*vol[*, *, level, *]
+         endcase
+      END
+      ELSE: return,  report('cas non code')
+   ENDCASE
+   varunit = 'cm'
+   varname = 'Dynamic Height (href='+strtrim(round(profref), 1)+'m)'
+   IF keyword_set(key_performance) THEN print, 'temps hdyn', systime(1)-tempsun 
+
+   return, hdyn
+end
Index: trunk/ToBeReviewed/CALCULS/level2depth.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/level2depth.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/level2depth.pro	(revision 25)
@@ -0,0 +1,77 @@
+;
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:level2depht
+;
+; PURPOSE: permet de passer d''un tableau 2d de niveau au
+; tableau 2d de profondeur correspondant a ces niveaux
+;
+; CATEGORY: SANS BOUCLE
+;
+; CALLING SEQUENCE: res=level2depth(niveau)
+;
+; INPUTS: niveau tableau 2d de niveaux seuil (ou une structure repondant
+; aux criteres de litchamp)
+;
+; KEYWORD PARAMETERS:
+;       /NOMASK: pour ne pas masquer les points terres
+;
+; OUTPUTS: un tableau 2d contenant des profondeurs
+;
+; COMMON BLOCKS:
+;	common.pro
+;
+; SIDE EFFECTS:
+;
+; RESTRICTIONS:
+;
+; EXAMPLE:
+;
+; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
+;                       17/6/1999
+;                       14/6/2000 accepte !values.f_nan
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION level2depth, tab, NOMASK = nomask
+   tempsun = systime(1)         ; pour key_performance
+@common
+;------------------------------------------------------------
+; lecture du champ d''entree et recuperation de la taille du sous
+; domaine utilise
+;------------------------------------------------------------
+   niveaux = litchamp(tab)
+   grille,mask, -1, -1,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz
+;---------------------------------------------------------------
+; verification de la coherence entre la taille du tableau et le domaine definit par domdef
+;---------------------------------------------------------------
+   taille = size(niveaux)
+   if taille[0] NE 2 then return, report('le champ en entree doit contenir un tableau 2d')
+   case 1 of
+      taille[1] eq jpi and taille[2] eq jpj:niveaux=niveaux[firstx:lastx, firsty:lasty]
+      taille[1] eq  nx and taille[2] eq  ny:
+      else:return, report('Probleme d''adequation entre les tailles du domaine et celle du champ.')
+   endcase
+;---------------------------------------------------------------
+
+   wherenan = where(finite(niveaux, /nan) EQ 1)
+   if wherenan[0] NE -1 then niveaux[wherenan] = 0
+   niveaux = 0 > niveaux < (jpk-1)
+;---------------------------------------------------------------
+   gdep = replicate(1, nx*ny)#gdep
+   niveaux = lindgen(nx*ny)+(nx*ny)*niveaux
+   gdep = reform(gdep[niveaux], nx, ny)
+;
+   if wherenan[0] NE -1 then gdep[wherenan] = !values.f_nan
+   if NOT keyword_set(nomask) then begin
+      if n_elements(valmask) EQ 0 then valmask = 1e20
+      terre = where(mask[*, *, 0] EQ 0)
+      if terre[0] NE -1 then gdep[terre] = valmask
+   endif
+;
+   if keyword_set(key_performance) THEN print, 'temps level2depth', systime(1)-tempsun 
+   return, gdep
+end
Index: trunk/ToBeReviewed/CALCULS/level2index.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/level2index.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/level2index.pro	(revision 25)
@@ -0,0 +1,53 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:level2index
+;
+; PURPOSE: on veut ds une matrice 3d extraire un tableau 2d (x,y) dont
+; chacun des elements a ete extrait a un niveau specifie par le
+; tableau 2d level (typiquement, on veut obtenir la salinite le long
+; d''une isopycne que l''on a reperee par son niveau).
+; level2index est une fonction qui donne en fonction de level un
+; tableau 2d d'indice qui permettra d''extraire le tableau 2d du
+; tableau 3d...
+;
+; CATEGORY: SANS BOUCLE
+;
+; CALLING SEQUENCE: index=level2index(level)
+; 
+; INPUTS:level: un tableau 2d de niveaux
+;
+; KEYWORD PARAMETERS:
+;
+; OUTPUTS: untableau 2d d''indices
+;
+; COMMON BLOCKS:
+;
+; SIDE EFFECTS:
+;
+; RESTRICTIONS:
+;
+; EXAMPLE:
+;
+; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
+;                      24/11/1999
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION level2index, level
+; un elements de tableau 3d (dont les 2 premieres dimensions sont nx et
+; ny) dont les coordonnes sont i, j et k a pour indice ds le meme
+; tableau 3d  i + j*nx + k*(nx*ny)
+; level etant donne, pour chaque points de level on connait i, j et k,
+; on peut donc calculer l''indice.
+   taille = size(level)
+   nx = taille[1]
+   ny = taille[2]
+; tableau k*(nx*ny)
+   tabknxny = (nx*ny)*long(level)
+;
+;
+   return, lindgen(nx, ny)+tabknxny
+end
Index: trunk/ToBeReviewed/CALCULS/level2mask.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/level2mask.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/level2mask.pro	(revision 25)
@@ -0,0 +1,75 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:level2mask
+;
+; PURPOSE: permet de passer d''un tableau 2d de niveau seuil au
+; tableau 3d de mask avec des 1 ds les niveaux au dessus du
+; niveau seuil et des 0 en dessous (et sur).
+;
+; CATEGORY: SANS BOUCLE
+;
+; CALLING SEQUENCE: res=level2mask(niveau)
+;
+; INPUTS: niveau tableau 2d de niveaux seuil (ou une structure repondant
+; aux criteres de litchamp)
+;
+; KEYWORD PARAMETERS:
+;
+; OUTPUTS: un tableau 3d contenant le mask associe au tableau 2d de
+; niveaux seuil
+;
+; COMMON BLOCKS:
+;	common.pro
+;
+; SIDE EFFECTS:
+;
+; RESTRICTIONS:
+;
+; EXAMPLE:
+;
+; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
+;                       17/6/1999
+; Setp 2004: boundary level have 0 values and not 1 (as it was
+; explained before in the header). see:
+; print, array_equal(niveau, total(level2mask(niveau),3))
+;
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION level2mask, tab
+;------------------------------------------------------------
+   tempsun = systime(1)         ; pour key_performance
+@common
+;------------------------------------------------------------
+; lecture du champ d''entree et recuperation de la taille du sous
+; domaine utilise
+;------------------------------------------------------------
+   niveaux = litchamp(tab)
+   grille,maskterre, -1, -1, -1,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz
+;---------------------------------------------------------------
+; verification de la coherence entre la taille du tableau et le domaine definit par domdef
+;---------------------------------------------------------------
+   IF ny EQ 1 THEN niveaux = reform(niveaux, nx, ny, /over)
+   taille = size(niveaux)
+   if taille[0] NE 2 then return, report('le champ en entree doit contenir un tableau 2d')
+   case 1 of
+      taille[1] eq jpi and taille[2] eq jpj:niveaux=niveaux[firstx:lastx, firsty:lasty]
+      taille[1] eq  nx and taille[2] eq  ny:
+      else:return, report('Probleme d''adequation entre les tailles du domaine et celle du champ.')
+   endcase
+;------------------------------------------------------------
+;------------------------------------------------------------
+; on transforme le tableau 2d de niveaux en tableau 3d de mask
+;------------------------------------------------------------
+   mask = reform( niveaux[*]#(1./(indgen(nz)+1)), nx, ny, nz )
+   mask = floor(temporary(mask)) < 1
+   mask = temporary(mask)*temporary(maskterre)
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+   if keyword_set(key_performance) THEN print, 'temps level2mask', systime(1)-tempsun 
+   return, mask
+end
Index: trunk/ToBeReviewed/CALCULS/moyenne.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/moyenne.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/moyenne.pro	(revision 25)
@@ -0,0 +1,600 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; 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
+;---------------------------------------------------------
+@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
Index: trunk/ToBeReviewed/CALCULS/norme.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/norme.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/norme.pro	(revision 25)
@@ -0,0 +1,454 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:norme
+;
+; PURPOSE: calcule la norme d''un champ de vecteurs, puis fait une
+; moyenne eventuelle.
+;          Rq1: le champ de vecteur peut etre, 2d:xy, 3d: xyz ou xyt,
+;          4d: xyzt
+;          Rq2: le calcul de la norme est fait avant l''eventuelle
+;          moyenne spatiale ou temporelle car la moyenne de la norme
+;          n''est pas egale a la norme des moyennes.
+;
+; CATEGORY: calcul de post traitement
+;
+; CALLING SEQUENCE:res=norme(champ_de_vecteurs)
+;
+; INPUTS:un tableau 2d, 3d ou 4d
+;
+; KEYWORD PARAMETERS:
+;
+;       BOXZOOM: boxzoom sur laquelle moyenner (par defaut le domaine
+;       selectionner par le dernier domdef effectue)
+;
+;       DIREC:'t' 'x' 'y' 'z' 'xy' 'xz' 'yz' 'xyz' 'xt' 'yt' 'zt' 'xyt'
+;       'xzt' 'yzt' 'xyzt' directions selon lesquelles effectuer les
+;       moyennes 
+;
+; OUTPUTS:tableau a tracer avec plt, pltz ou pltt.
+;
+; COMMON BLOCKS:
+;	common.pro
+;
+; SIDE EFFECTS:
+;
+;      La norme est calculee aux points T. Pour faire ce calcul, on
+;      moyenne les champs U et V aux points T avant de calculer la norme.
+;      Au bord des cotes et du domaine, on ne peut pas calculer les
+;      champs U et V aux points T, ces points sont donc a la valeur
+;      !values.f_nan. 
+; 
+;      lorsqu''on fait le calcul sur un domaine geographique reduit,
+;      les champs U et V ne comprennent pas forcement le meme nombre
+;      de points. Dans ce cas on redecoupe U et V pour ne garder que
+;      les points en commun. Au passage on refait un domdef qui
+;      redefinit un domaine geographique sur lequel les champs U et V
+;      sont extraits sur les meme points.
+;
+; RESTRICTIONS:
+;
+;      pour savoir a quel type de tableau on a a faire, on teste la
+;      taille de celui-ci et les dates donnees par time[0] et
+;      time[jpt-1] pour savoir si il y a une dimension
+;      temporelle. Avant de lancer norme s''assurer que time et jpt
+;      sont bien definis comme il faut! 
+;
+; EXAMPLE:
+;
+;     pour calculer la moyenne de la norme des courants sur tout le
+;     dommaine entre 0 et 50:
+;      IDL> res=norme(un,vn,boxzoom=[0,50],dir='xyz')
+;
+; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
+;                       9/6/1999
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION norme, composanteu, composantev, BOXZOOM = boxzoom, DIREC = direc, _extra = ex
+;---------------------------------------------------------
+@cm_4mesh
+@cm_4data
+@cm_4cal
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+@updatekwd
+  ENDIF
+;---------------------------------------------------------
+   tempsun = systime(1)         ; pour key_performance
+;
+   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
+     return, report(['This version of norme is based on Arakawa C-grid.' $
+                     , 'U and V grids must therefore be defined'])
+;
+;------------------------------------------------------------
+  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
+    domdef, boxzoom
+  ENDIF 
+;------------------------------------------------------------
+   if NOT keyword_set(direc) then direc = 0
+; construction de u et v aux pts T
+   u = litchamp(composanteu)
+   v = litchamp(composantev)
+   date1 = time[0]
+   if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
+
+   if (size(u))[0] NE (size(v))[0] then return,  -1
+
+   vargrid='T'
+   varname = 'norme '
+   valmask = 1e20
+;
+   grilleu = litchamp(composanteu, /grid)
+   if grilleu EQ '' then grilleu = 'U'
+   grillev = litchamp(composantev, /grid)
+   if grillev EQ '' then grillev = 'V'
+   IF grilleu EQ 'V' AND grillev EQ 'U' THEN inverse = 1
+   IF grilleu EQ 'T' AND grillev EQ 'T' THEN BEGIN 
+      interpolle  = 0 
+      return, report('cas non code mais facile a faire!')
+   ENDIF ELSE interpolle = 1
+   if keyword_set(inverse) then begin
+      rien = u
+      u = v
+      v = rien
+   endif
+
+
+;------------------------------------------------------------
+; on trouve les points que u et v ont en communs
+;------------------------------------------------------------
+   indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
+   indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
+   indicex = inter(indicexu, indicexv)
+   indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
+   indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
+   indicey = inter(indiceyu, indiceyv)
+   nx = n_elements(indicex) 
+   ny = n_elements(indicey)
+;----------------------------------------------------------------------------
+   case 1 of
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xyz
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+      (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN 
+;----------------------------------------------------------------------------
+         indice2d = lindgen(jpi, jpj)
+         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
+         indice3d = lindgen(jpi, jpj, jpk)
+         indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt]
+;------------------------------------------------------------
+; extraction de u et v sur le domaine qui convient
+;------------------------------------------------------------
+         case 1 of
+            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
+             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
+               case (size(u))[3] OF
+                  nzt:BEGIN 
+                     if nxu NE nx then $
+                      if indicex[0] EQ firstxu then u = u[0:nx-1,*,*] ELSE u = u[1: nx, *,*] 
+                     IF nxv NE nx THEN $
+                      if indicex[0] EQ firstxv then v = v[0:nx-1,*,*] ELSE v = v[1: nx, *,*]
+                     IF nyu NE ny THEN $
+                      if indicey[0] EQ firstyu then u = u[*,0:ny-1,*] ELSE u = u[*, 1: ny,*] 
+                     IF nyv NE ny THEN $
+                      if indicey[0] EQ firstyv then v = v[*,0:ny-1,*] ELSE v = v[*, 1: ny,*]
+                  end
+                  jpk:BEGIN 
+                     if nxu NE nx then $
+                      if indicex[0] EQ firstxu then u = u[0:nx-1, *,firstzt:lastzt] ELSE u = u[1: nx, *,firstzt:lastzt] 
+                     IF nxv NE nx THEN $
+                      if indicex[0] EQ firstxv then v = v[0:nx-1, *,firstzt:lastzt] ELSE v = v[1: nx, *,firstzt:lastzt]
+                     IF nyu NE ny THEN $
+                      if indicey[0] EQ firstyu then u = u[*, 0:ny-1,firstzt:lastzt] ELSE u = u[*, 1: ny,firstzt:lastzt] 
+                     IF nyv NE ny THEN $
+                      if indicey[0] EQ firstyv then v = v[*, 0:ny-1,firstzt:lastzt] ELSE v = v[*, 1: ny,firstzt:lastzt]
+                  end
+                  ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
+               endcase
+            END
+            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $
+             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN 
+               u = u[indice3d]
+               v = v[indice3d]
+            END
+            ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
+         endcase
+;------------------------------------------------------------------
+; on reform u et v pour s'assurer qu'aucune dimension n'a ete ecrasee
+;------------------------------------------------------------------
+         if nzt EQ 1 then begin
+            u = reform(u, nx, ny, nzt, /over)
+            v = reform(v, nx, ny, nzt, /over)
+         endif
+;------------------------------------------------------------------
+; construction de u et v aux pts T
+;-----------------------------------------------------------
+         a=u(0,*,*)
+         u=(u+shift(u,1,0,0))/2.
+         if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*,*)=a
+         a=v(*,0,*)
+         v=(v+shift(v,0,1,0))/2.
+         if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0,*)=a
+;----------------------------------------------------------------------------
+; attribution du mask et des tableau de longitude et latitude
+;----------------------------------------------------------------------------
+         mask = tmask[indice3d]
+         if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over)
+;-----------------------------------------------------------
+         if n_elements(valmask) EQ 0 THEN valmask = 1e20
+         landu = where(u GE valmask/10)
+         if landu[0] NE -1 then u[landu] = 0
+         landv = where(v GE valmask/10)
+         if landv[0] NE -1 then v[landv] = 0
+         res=sqrt(u^2+v^2)
+         if NOT keyword_set(key_periodic) OR nx NE jpi then res(0,*, *)=!values.f_nan
+         res(*,0, *)=!values.f_nan
+         mask = where(mask eq 0)
+         IF mask[0] NE -1 THEN res(mask) = valmask
+; moyennes en tous genres
+         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
+         if keyword_set(direc) then res = moyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef)
+;-----------------------------------------------------------
+      END
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xyt
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+      date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN 
+         indice2d = lindgen(jpi, jpj)
+         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
+;------------------------------------------------------------
+; extraction de u et v sur le domaine qui convient
+;------------------------------------------------------------
+         case 1 of
+            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
+             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
+               if nxu NE nx then $
+                if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
+               IF nxv NE nx THEN $
+                if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
+               IF nyu NE ny THEN $
+                if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
+               IF nyv NE ny THEN $
+                if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
+            END
+            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
+             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
+               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
+               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
+            END
+            ELSE:return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
+         endcase
+;------------------------------------------------------------------
+; construction de u et v aux pts T
+;-----------------------------------------------------------
+         a=u(0,*,*)
+         u=(u+shift(u,1,0,0))/2.
+         if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*,*)=a
+         a=v(*,0,*)
+         v=(v+shift(v,0,1,0))/2.
+         if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0,*)=a
+;----------------------------------------------------------------------------
+; attribution du mask et des tableau de longitude et latitude
+; on recupere la grille complette pour etablir un grand mask etendu ds les 4
+; directions pour couvrir les points pour lesquels un pt terre a ete pris en
+; compte (faire un petit dessin...)
+;----------------------------------------------------------------------------
+         mask = tmask[indice2d+jpi*jpj*firstzt]
+         if ny EQ 1 then mask = reform(mask, nx, ny, /over)
+;-----------------------------------------------------------
+; construction de terre qui contient tous les point a masquer
+;-----------------------------------------------------------
+         if n_elements(valmask) EQ 0 THEN valmask = 1e20
+         landu = where(u GE valmask/10)
+         if landu[0] NE -1 then u[landu] = 0
+         landv = where(v GE valmask/10)
+         if landv[0] NE -1 then v[landv] = 0
+         res=sqrt(u^2+v^2)
+         if NOT keyword_set(key_periodic) OR nx NE jpi then res(0,*, *)=!values.f_nan
+         res(*,0, *)=!values.f_nan
+         mask = where(mask eq 0)
+         IF mask[0] NE -1 THEN BEGIN 
+            coeftps = lindgen(jpt)*nx*ny
+            coeftps = replicate(1, n_elements(mask))#coeftps
+            mask = (temporary(mask))[*]#replicate(1, jpt)
+            mask =temporary(mask[*]) + temporary(coeftps[*])
+            res(temporary(mask)) = valmask
+         ENDIF
+; moyennes en tous genres
+         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
+         if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef)
+      END
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xyzt
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+      date1 NE date2 AND (size(u))[0] EQ 4:BEGIN 
+         indice2d = lindgen(jpi, jpj)
+         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
+         indice3d = lindgen(jpi, jpj, jpk)
+         indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt]
+;------------------------------------------------------------
+; extraction de u et v sur le domaine qui convient
+;------------------------------------------------------------
+         case 1 of
+            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
+             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
+               case (size(u))[3] OF
+                  nzt:BEGIN 
+                     if nxu NE nx then $
+                      if indicex[0] EQ firstxu then u = u[0:nx-1,*,*,*] ELSE u = u[1: nx, *,*,*] 
+                     IF nxv NE nx THEN $
+                      if indicex[0] EQ firstxv then v = v[0:nx-1,*,*,*] ELSE v = v[1: nx, *,*,*]
+                     IF nyu NE ny THEN $
+                      if indicey[0] EQ firstyu then u = u[*,0:ny-1,*,*] ELSE u = u[*, 1: ny,*,*] 
+                     IF nyv NE ny THEN $
+                      if indicey[0] EQ firstyv then v = v[*,0:ny-1,*,*] ELSE v = v[*, 1: ny,*,*]
+                  end
+                  jpk:BEGIN 
+                     if nxu NE nx then $
+                      if indicex[0] EQ firstxu then u = u[0:nx-1, *,firstzt:lastzt,*] ELSE u = u[1: nx, *,firstzt:lastzt,*] 
+                     IF nxv NE nx THEN $
+                      if indicex[0] EQ firstxv then v = v[0:nx-1, *,firstzt:lastzt,*] ELSE v = v[1: nx, *,firstzt:lastzt,*]
+                     IF nyu NE ny THEN $
+                      if indicey[0] EQ firstyu then u = u[*, 0:ny-1,firstzt:lastzt,*] ELSE u = u[*, 1: ny,firstzt:lastzt,*] 
+                     IF nyv NE ny THEN $
+                      if indicey[0] EQ firstyv then v = v[*, 0:ny-1,firstzt:lastzt,*] ELSE v = v[*, 1: ny,firstzt:lastzt,*]
+                  end
+                  ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
+               endcase
+            END
+            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $
+             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN 
+               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt, *]
+               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt, *]
+            END
+            ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
+         endcase
+;------------------------------------------------------------------
+; construction de u et v aux pts T
+;-----------------------------------------------------------
+         a=u(0,*,*,*)
+         u=(u+shift(u,1,0,0,0))/2.
+         if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*,*,*)=a
+         a=v(*,0,*,*)
+         v=(v+shift(v,0,1,0,0))/2.
+         if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0,*,*)=a
+;----------------------------------------------------------------------------
+; attribution du mask et des tableau de longitude et latitude
+;----------------------------------------------------------------------------
+         mask = tmask[indice3d]
+         if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over)
+;-----------------------------------------------------------
+         if n_elements(valmask) EQ 0 THEN valmask = 1e20
+         landu = where(u GE valmask/10)
+         if landu[0] NE -1 then u[landu] = 0
+         landv = where(v GE valmask/10)
+         if landv[0] NE -1 then v[landv] = 0
+         res=sqrt(u^2+v^2)
+         if NOT keyword_set(key_periodic) OR nx NE jpi then res(0,*, *, *)=!values.f_nan
+         res(*,0, *, *)=!values.f_nan
+         mask = where(mask eq 0)
+         IF mask[0] NE -1 THEN BEGIN 
+            coeftps = lindgen(jpt)*nx*ny*nzt
+            coeftps = replicate(1, n_elements(mask))#coeftps
+            mask = (temporary(mask))[*]#replicate(1, jpt)
+            mask =temporary(mask[*]) + temporary(coeftps[*])
+            res(temporary(mask)) = valmask
+         ENDIF
+; moyennes en tous genres
+         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
+         if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef)
+      END
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+;xy
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+      ELSE:BEGIN                ;xy
+         indice2d = lindgen(jpi, jpj)
+         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
+;------------------------------------------------------------
+; extraction de u et v sur le domaine qui convient
+;------------------------------------------------------------
+         case 1 of
+            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
+             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
+               if nxu NE nx then $
+                if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
+               IF nxv NE nx THEN $
+                if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
+               IF nyu NE ny THEN $
+                if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
+               IF nyv NE ny THEN $
+                if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
+            END
+            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
+             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
+               u = u[indice2d]
+               v = v[indice2d]
+            END
+            ELSE:return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
+         endcase
+;------------------------------------------------------------------
+; on reform u et v pour s'assurer qu'aucune dimension n'a ete ecrasee
+;------------------------------------------------------------------
+         if ny EQ 1 then begin
+            u = reform(u, nx, ny, /over)
+            v = reform(v, nx, ny, /over)
+         endif
+;------------------------------------------------------------------
+; construction de u et v aux pts T
+;-----------------------------------------------------------
+         a=u(0,*)
+         u=(u+shift(u,1,0))/2.
+         if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*)=a
+         a=v(*,0)
+         v=(v+shift(v,0,1))/2.
+         if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0)=a
+;----------------------------------------------------------------------------
+; attribution du mask et des tableau de longitude et latitude
+; on recupere la grille complette pour etablir un grand mask etendu ds les 4
+; directions pour couvrir les points pour lesquels un pt terre a ete pris en
+; compte (faire un petit dessin...)
+;----------------------------------------------------------------------------
+         mask = tmask[indice2d+jpi*jpj*firstzt]
+         if nyt EQ 1 THEN mask = reform(mask, nx, ny, /over)
+;-----------------------------------------------------------
+; construction de terre qui contient tous les point a masquer
+;-----------------------------------------------------------
+         if n_elements(valmask) EQ 0 THEN valmask = 1e20
+         landu = where(u GE valmask/10)
+         if landu[0] NE -1 then u[landu] = 0
+         landv = where(v GE valmask/10)
+         if landv[0] NE -1 then v[landv] = 0
+         res=sqrt(u^2+v^2)
+         if NOT keyword_set(key_periodic) OR nx NE jpi then res(0,*)=!values.f_nan
+         res(*,0)=!values.f_nan
+         mask = where(mask eq 0)
+         IF mask[0] NE -1 THEN res(mask) = valmask
+; moyennes en tous genres
+         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
+         if keyword_set(direc) then res = moyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef)
+      END
+;----------------------------------------------------------------------------
+   endcase
+;------------------------------------------------------------
+   if keyword_set(key_performance) THEN print, 'temps norme', systime(1)-tempsun 
+   return, res
+end
Index: trunk/ToBeReviewed/CALCULS/projectondepth.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/projectondepth.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/projectondepth.pro	(revision 25)
@@ -0,0 +1,108 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:projectondepth
+;
+; PURPOSE: routine permettant de projeter un champ 3d suivant un
+; tableau de profondeurs.
+;
+; CATEGORY: sans boucles
+;
+; CALLING SEQUENCE:res=projectondepth(arrayin, depthin)
+; 
+; INPUTS:
+;       arrayin: un tableau 3d dont la 3eme dimension doit etre egale
+;       a jpk
+;       depthin: un tableau 2d indiquant n chaque point a quel
+;       profondeur projeter
+;       
+; KEYWORD PARAMETERS:none
+;
+; OUTPUTS:res: un tableau 2d projection du tableau 3d suivant les
+; profondeurs indiquees par depthin
+;
+; COMMON BLOCKS:common.pro
+;
+; SIDE EFFECTS: points a !values.f_nan qd calcul impossible. points
+; terres masques a Valmask.
+;
+; RESTRICTIONS:
+;
+; EXAMPLE:
+;
+; on contruit un tableau de profondeurs possibles
+;   IDL> a=gdept[jpk-1]/(1.*jpi*jpj)*findgen(jpi,jpj)
+; on contruit un tableau a projeter sur ces profondeurs. pour le test
+; on construit un tableau 3d dont chaque vecteur suivant z est la
+; profondeur. 
+;   IDL> arraytest=replicate(1,jpi*jpj)#gdept
+;   IDL> arraytest=reform(arraytest,jpi,jpj,jpk, /over)
+; on test la projection du tabeau profondeur sur la profondeur...
+;   IDL> plt, 1e6*(a-projectondepth(arraytest,a)),/nocontour
+;   ->champ nul a 1e-6 pres
+; 
+;  verifcation en projettant la temperature sur la profondeur de la 20
+;  degres par exemple...
+;
+; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
+;                      15/6/2000
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION projectondepth, arrayin, depthin
+   tempsun = systime(1)         ; pour key_performance
+@common
+;------------------------------------------------------------
+   depth = litchamp(depthin)
+   array = litchamp(arrayin)
+; petites verifications
+   tailledepth = size(depth)
+   taillearray = size(array)
+   if tailledepth[0] NE 2 THEN return, report('Depth array must have 2 dimensions')
+   if taillearray[0] NE 3 THEN return, report('Array in must have 3 dimensions')
+; verification de la coherence entre la taille du tableau et le
+; domaine 
+   grille, mask, -1, -1, -1,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz
+   case 1 of
+      tailledepth[1] eq jpi and tailledepth[2] eq jpj:depth=depth[firstx:lastx, firsty:lasty]
+      tailledepth[1] eq  nx and tailledepth[2] eq  ny:
+      else:return, report('Probleme d''adequation entre les tailles du domaine et celle du tableau de profondeur')
+   endcase
+   case 1 OF
+      taillearray[3] NE jpk:return, report('Le tableau 3d doit avoir sa 3eme dimension egale a jpk')
+      taillearray[1] eq jpi and taillearray[2] eq jpj:array=array[firstx:lastx, firsty:lasty, *]
+      taillearray[1] eq  nx and taillearray[2] eq  ny:
+      else:return, report('Probleme d''adequation entre les tailles du domaine et celle du tableau de profondeur')
+   endcase
+;
+; c''est parti
+;
+   flevel = depth2floatlevel(depth)
+; on vire les points a !values.f_nan
+   notanumber = where(finite(flevel, /nan) EQ 1)
+   if notanumber[0] NE -1 then flevel[notanumber] = 0
+; on seuil (vire les points terres a valmask par ex)
+   flevel = 0 > flevel < (jpk-1)
+;
+   indexup = level2index(floor(flevel))
+   indexlow = nx*ny+indexup
+   out = where(indexlow GE nx*ny*jpk-1)
+   if out[0] NE -1 then indexlow[out] = indexlow[out]-nx*ny
+;
+   weight = flevel-floor(flevel)
+   res = array[indexup]
+   res = res+weight*(array[indexlow]-res)
+;
+; on replace les points a !values.f_nan
+   if notanumber[0] NE -1 then res[notanumber] = !values.f_nan
+   if out[0] NE -1 then res[out] = !values.f_nan
+; on masque les points terres a valmask
+   if n_elements(valmask) EQ 0 then valmask = 1e20
+   terre = where((temporary(mask))[*, *, 0] EQ 0)
+   if terre[0] NE -1 then res[terre] = valmask
+;------------------------------------------------------------
+   if keyword_set(key_performance) THEN print, 'temps projectondepth', systime(1)-tempsun 
+   return, res
+end
Index: trunk/ToBeReviewed/CALCULS/remplit.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/remplit.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/remplit.pro	(revision 25)
@@ -0,0 +1,233 @@
+FUNCTION remplit, zinput, NAN = nan, NITER = niter, BASIQUE = basique, mask = mask, FILLXDIR = fillxdir, FILLYDIR = fillydir, FILLVAL = fillval, _extra = ex
+@common
+  tempsun = systime(1)          ; pour key_performance
+;+
+  ;;
+  ;; Extrapole zinout[jpi,jpj] sur les continents en utilisant les 4
+  ;; plus proches voisins masques oceaniquement et construit un nouveau
+;   masque
+  ;; contenant l'ancien masque oceanique PLUSles points extrapoles.
+  ;; Reitere le processus niter fois.
+  ;; C'est pas clair, essayez !
+  ;;
+  ;;
+;
+;    /Nan: to fill the point which have the value
+;    !values.f_nan. Whitout this keyword, these point are not filling
+;    and stays at !values.f_nan.
+;
+;
+;-
+; les points non remplis sont masques a valmask
+  IF n_elements(niter) EQ 0 THEN niter = 1
+  IF niter EQ 0 THEN return, zinput
+  z = zinput
+  if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c'
+  if keyword_set(basique) then begin
+    oldkey_gridtype = key_gridtype
+    key_gridtype = 'c'
+    nx = (size(zinput))[1]
+    ny = (size(zinput))[2]
+    if NOT keyword_set(mask) then mmmask = basique ELSE mmmask = mask
+    if  key_gridtype eq 'e' then begin
+      case vargrid of
+        'T':glam = glamt[firstxt:lastxt, firstyt:lastyt]
+        'U':glam = glamu[firstxu:lastxu, firstyu:lastyu]
+      endcase
+    endif
+  ENDIF ELSE grille, mmmask, glam, gphi, gdep, nx, ny, nz, _extra = ex
+  if keyword_set(mask) then mmmask = mask
+;---------------------------------------------------------------
+  if (size(mmmask))[0] EQ 3 THEN mmmask = mmmask[*, *, 0]
+;
+  if n_elements(mmmask) EQ 1 then mmmask = replicate(1b, nx, ny)
+  if keyword_set(nan) then begin
+    nanpoint = where(finite(z) EQ 0)
+    if nanpoint[0] NE -1 then begin
+      mmmask[nanpoint] = 0b
+      z[nanpoint] = 0
+    endif
+  ENDIF
+  mmmask = byte(mmmask)
+;---------------------------------------------------------------
+; on ajoute un cadre de zero a z, mask, e1, e2
+; comme ca apres on peut faire des shifts ds tous les sens sans se
+; soucier des bords du domaine!
+;---------------------------------------------------------------
+  tempdeux = systime(1)         ; pour key_performance =2
+  nx2 = nx+2
+  case key_gridtype of
+    'c':BEGIN
+      ztmp = bytarr(nx+2, ny+2)
+      ztmp[1:nx, 1:ny] = mmmask
+      mmmask = temporary(ztmp)
+      ztmp = fltarr(nx+2, ny+2)
+      ztmp[1:nx, 1:ny] = z
+      if keyword_set(key_periodic) AND nx EQ jpi then begin
+        ztmp[0, 1:ny] = z[jpi-1, *]
+        ztmp[nx+1, 1:ny] = z[0, *]
+      endif
+      z = temporary(ztmp)
+    END
+    'e':BEGIN
+      ztmp = bytarr(nx+2, ny+4)
+      ztmp[1:nx, 2:ny+1] = mmmask
+      mmmask = temporary(ztmp)
+      ztmp = fltarr(nx+2, ny+4)
+      ztmp[1:nx, 2:ny+1] = z
+      if keyword_set(key_periodic) AND nx EQ jpi then begin
+        ztmp[0, 2:ny+1] = z[jpi-1, *]
+        ztmp[nx+1, 2:ny+1] = z[0, *]
+      endif
+      z = temporary(ztmp)
+    END
+  endcase
+  IF testvar(var = key_performance) EQ 2 THEN $
+    print, 'temps remplit: on ajoute un cadre de zero ', systime(1)-tempdeux
+;---------------------------------------------------------------
+;---------------------------------------------------------------
+; iteration 
+;---------------------------------------------------------------
+;---------------------------------------------------------------
+  FOR n = 1, niter DO BEGIN 
+; on trouve les points coast
+    tempdeux = systime(1)       ; pour key_performance =2
+; les points du bord du cadre ne doivent pas etre selectionnes comme
+; la coast
+    case key_gridtype of
+      'c':BEGIN
+        mmmask[0, *] = 1b
+        mmmask[nx+1, *] = 1b
+        mmmask[*, 0] = 1b
+        mmmask[*, ny+1] = 1b
+      END
+      'e':BEGIN
+        mmmask[0, *] = 1b
+        mmmask[nx+1, *] = 1b
+        mmmask[*, 0:1] = 1b
+        mmmask[*, ny+2:ny+3] = 1b
+      END
+    endcase
+; liste des points terre restant
+    IF keyword_set(fillxdir) THEN BEGIN
+; we stop if all the lines, that contains data, have been filled
+      test = total(mmmask[1:nx, *], 1)
+      IF total((test EQ 0)+(test EQ nx)) EQ ny+2 THEN GOTO, fini
+    ENDIF
+    IF keyword_set(fillydir) THEN BEGIN
+; we stop if all the columns, that contains data, have been filled
+      test = total(mmmask[*, 1:ny], 2)
+      IF total((test EQ 0)+(test EQ ny)) EQ nx+2 THEN GOTO, fini
+    ENDIF
+    land = where(mmmask EQ 0)
+    if land[0] EQ -1 then GOTO, fini
+; les points du bord du cadre doivent maintenant etre dans la terre
+    case key_gridtype of
+      'c':BEGIN
+        mmmask[0, *] = 0b
+        mmmask[nx+1, *] = 0b
+        mmmask[*, 0] = 0b
+        mmmask[*, ny+1] = 0b
+      END
+      'e':BEGIN
+        mmmask[0, *] = 0b
+        mmmask[nx+1, *] = 0b
+        mmmask[*, 0:1] = 0b
+        mmmask[*, ny+2:ny+3] = 0b
+      END
+    endcase
+    if keyword_set(key_periodic) AND nx EQ jpi then begin
+      mmmask[0, *] = mmmask[nx, *]
+      mmmask[nx+1, *] = mmmask[1, *]
+    endif
+; liste des voisins mer
+    case key_gridtype of
+      'c':BEGIN
+        CASE 1 OF
+          keyword_set(fillxdir):weight = mmmask[1+land]+mmmask[-1+land]
+          keyword_set(fillydir):weight = mmmask[nx2+land]+mmmask[-nx2+land]
+          ELSE:weight = mmmask[1+land]+mmmask[-1+land]+mmmask[nx2+land]+mmmask[-nx2+land] $
+            +1./sqrt(2)*(mmmask[nx2+1+land]+mmmask[nx2-1+land] $
+                         +mmmask[-nx2+1+land]+mmmask[-nx2-1+land])
+        ENDCASE
+      END
+      'e':BEGIN
+        shifted = glam[0, 0] LT glam[0, 1]
+        oddeven = (land/nx2+1-shifted) MOD 2
+        weight = mmmask[1+land]+mmmask[-1+land] $
+          +mmmask[2*nx2+land]+mmmask[-2*nx2+land] $
+          +sqrt(2)*(mmmask[-nx2+oddeven+land]+mmmask[-nx2-1+oddeven+land] $
+                    +mmmask[nx2+oddeven+land]+mmmask[nx2-1+oddeven+land])
+      END
+    endcase
+
+    ok = where(weight GT 0)
+    weight = weight[ok]
+    coast = land[temporary(ok)]
+;
+    IF testvar(var = key_performance) EQ 2 THEN $
+      print, 'temps remplit: trouver la coast ', systime(1)-tempdeux
+;---------------------------------------------------------------
+; remplissage des points coast
+;---------------------------------------------------------------
+    tempdeux = systime(1)       ; pour key_performance =2
+; on masque z
+    z = temporary(z)*mmmask
+;
+    case key_gridtype of
+      'c':BEGIN
+        CASE 1 OF
+          keyword_set(fillxdir):zcoast = z[1+coast]+z[-1+coast]
+          keyword_set(fillydir):zcoast = z[nx2+coast]+z[-nx2+coast]
+          ELSE:zcoast = z[1+coast]+z[-1+coast]+z[nx2+coast]+z[-nx2+coast] $
+            +1./sqrt(2)*(z[nx2+1+coast]+z[nx2-1+coast] $
+                         +z[-nx2+1+coast]+z[-nx2-1+coast])
+        ENDCASE
+      END
+      'e':BEGIN
+        oddeven = (coast/nx2+1-shifted) MOD 2
+        zcoast = z[1+coast]+z[-1+coast]+z[2*nx2+coast]+z[-2*nx2+coast] $
+          +sqrt(2)*(z[-nx2+oddeven+coast]+z[-nx2-1+oddeven+coast] $
+                    +z[nx2+oddeven+coast]+z[nx2-1+oddeven+coast])
+      END
+    endcase
+;    
+    z[coast] =  temporary(zcoast)/ temporary(weight)
+; we update the the boundary conditions of z
+    if keyword_set(key_periodic) AND nx EQ jpi then begin
+      z[0, *] = z[nx, *]
+      z[nx+1, *] = z[1, *]
+    endif
+;---------------------------------------------------------------
+; IV) on reduit le masque
+;---------------------------------------------------------------
+    mmmask[ temporary(coast)] = 1
+;
+    IF testvar(var = key_performance) EQ 2 THEN $
+      print, 'temps remplit: une iteration ', systime(1)-tempdeux
+  ENDFOR 
+fini:
+;
+; on masque les valeurs sur les lands restantes
+;
+  IF n_elements(valmask) EQ 0 then valmask = 1e20
+  IF n_elements(fillval) EQ 0 THEN fillval = valmask
+  z = temporary(z)*mmmask + fillval*(1b-mmmask)
+;---------------------------------------------------------------
+; on redecoupe le tableau pour retirer le cadre!
+;---------------------------------------------------------------
+  case key_gridtype of
+    'c':BEGIN
+      z = z[1:nx, 1:ny]
+    END
+    'e':BEGIN
+      z = z[1:nx, 2:ny+1]
+    END
+  endcase
+;
+  if keyword_set(basique) then key_gridtype = oldkey_gridtype
+;---------------------------------------------------------------
+  if keyword_set(key_performance) THEN print, 'temps remplit', systime(1)-tempsun 
+  return, z
+END 
+
Index: trunk/ToBeReviewed/CALCULS/rhon.pro
===================================================================
--- trunk/ToBeReviewed/CALCULS/rhon.pro	(revision 25)
+++ trunk/ToBeReviewed/CALCULS/rhon.pro	(revision 25)
@@ -0,0 +1,147 @@
+FUNCTION rhon, sn, tn, INSITU = insitu, SIGMA_N = sigma_n
+@common
+   tempsun = systime(1)         ; pour key_performance
+;+
+   ;;
+   ;; Calcul de la fonction d'etat (issue de eos.F)
+   ;;
+   ;; Creation : 1997 / G. Roullet
+   ;; adaptation pour les tableaux z,zt,xyz,xyzt
+   ;; par seb.
+   ;;
+;-
+      sn = -1e5 > double(sn) < 1e5
+      tn = -1e5 > double(tn) < 1e5
+
+
+      IF keyword_set(sigma_n) then insitu = 1
+
+   taille = size(sn)
+   case taille[0] of
+      0:BEGIN                   ;z
+         zrhop=0d
+         jkmax = 1
+      END 
+      1:BEGIN                   ;z
+         zrhop=dblarr(taille[1])
+         jkmax = taille[1]
+      END 
+      2:BEGIN                   ;xy (jpt=1) ou zt
+         zrhop=dblarr(taille[1],taille[2])
+         if jpt EQ 1 then jkmax = 1 ELSE jkmax = taille[1]
+      END 
+      3:BEGIN                   ;xyz (jpt=1) ou xyt
+         zrhop=dblarr(taille[1],taille[2],taille[3])
+         if jpt EQ 1 then jkmax = taille[3] ELSE jkmax = 1
+      END 
+      4:BEGIN                   ;xyzt
+         zrhop=dblarr(taille[1],taille[2],taille[3],taille[4] )
+         jkmax = taille[3]
+      END 
+   endcase
+   
+
+
+   FOR jk = 0, jkmax-1 DO BEGIN
+
+      case taille[0] of
+         0:BEGIN                ;z
+            ztt = tn
+            zs = sn
+         END 
+         1:BEGIN                ;z
+            ztt = tn[jk]
+            zs = sn[jk] 
+         END 
+         2:BEGIN                ;xy (jpt=1) ou zt
+            if jpt EQ 1 then begin
+               ztt = tn
+               zs = sn
+            ENDIF ELSE BEGIN
+               ztt = tn[jk, *]
+               zs = sn[jk, *] 
+            ENDELSE
+         END 
+         3:BEGIN                ;xyz (jpt=1) ou xyt
+            if jpt EQ 1 then begin
+               ztt = tn[*, *, jk]
+               zs = sn[*, *,jk] 
+            endif ELSE BEGIN
+               ztt = tn
+               zs = sn
+            ENDELSE
+         END 
+         4:BEGIN                ;xyzt
+            ztt = tn[*, *, jk, *]
+            zs = sn[*, *,jk, *] 
+         END 
+      endcase
+      if n_elements(sigma_n) NE 0 then zh = 1000.*sigma_n ELSE zh = gdept(jk)
+; ...   square root salinity
+      zsr= sqrt(abs(zs))
+; ...   compute density pure water at atm pressure
+      zr1=((((6.536332e-9*ztt-1.120083e-6)*ztt+1.001685e-4)*ztt-9.095290e-3)*ztt+6.793952e-2)*ztt+999.842594
+; ...   seawater density atm pressure
+      zr2= (((5.3875e-9*ztt-8.2467e-7)*ztt+7.6438e-5)*ztt-4.0899e-3)*ztt+0.824493
+      zr3= (-1.6546e-6*ztt+1.0227e-4)*ztt-5.72466e-3
+      zr4= 4.8314e-4
+;
+; ...   potential density (reference to the surface)
+      case taille[0] of
+         0: zrhop    = (zr4*zs + zr3*zsr + zr2)*zs + zr1
+         1: zrhop(jk)= (zr4*zs + zr3*zsr + zr2)*zs + zr1
+         2:BEGIN 
+            if jpt EQ 1 then zrhop = (zr4*zs + zr3*zsr + zr2)*zs + zr1 $
+            ELSE zrhop(jk, *)= (zr4*zs + zr3*zsr + zr2)*zs + zr1
+         END
+         3:BEGIN 
+            if jpt EQ 1 then zrhop(*, *,jk)= (zr4*zs + zr3*zsr + zr2)*zs + zr1 $
+             ELSE zrhop = (zr4*zs + zr3*zsr + zr2)*zs + zr1
+         END
+         4: zrhop(*, *,jk, *)= (zr4*zs + zr3*zsr + zr2)*zs + zr1
+      endcase
+
+      IF n_elements(insitu) EQ 1 THEN BEGIN 
+; ...   add the compression terms
+         ze = (-3.508914e-8*ztt-1.248266e-8)*ztt-2.595994e-6
+         zbw= ( 1.296821e-6*ztt-5.782165e-9)*ztt+1.045941e-4
+         zb = zbw + ze * zs
+;
+         zd = -2.042967e-2
+         zc = (-7.267926e-5*ztt+2.598241e-3)*ztt+0.1571896
+         zaw = ((5.939910e-6*ztt+2.512549e-3)*ztt-0.1028859)*ztt -4.721788
+         za = ( zd*zsr + zc)*zs + zaw
+;
+         zb1= (-0.1909078*ztt+7.390729)*ztt-55.87545
+         za1= ((2.326469e-3*ztt+1.553190)*ztt-65.00517)*ztt+1044.077
+         zkw= (((-1.361629e-4*ztt-1.852732e-2)*ztt-30.41638)*ztt+2098.925)*ztt+190925.6
+         zk0= (zb1*zsr + za1)*zs + zkw
+;
+; ...   masked in situ density
+         case taille[0] of
+            0: zrhop = zrhop / (1.0-zh/(zk0-zh*(za-zh*zb)))
+            1: zrhop(jk) = zrhop(jk) / (1.0-zh/(zk0-zh*(za-zh*zb)))
+            2:BEGIN 
+               if jpt EQ 1 then zrhop = zrhop / (1.0-zh/(zk0-zh*(za-zh*zb))) $
+               ELSE zrhop(jk, *) = zrhop(jk, *) / (1.0-zh/(zk0-zh*(za-zh*zb)))
+            END
+            3:BEGIN 
+               if jpt EQ 1 then zrhop(*, *,jk) = zrhop(*, *,jk) / (1.0-zh/(zk0-zh*(za-zh*zb))) $
+               ELSE zrhop = zrhop / (1.0-zh/(zk0-zh*(za-zh*zb)))
+            END
+            4: zrhop(*, *,jk, *) = zrhop(*, *,jk, *) / (1.0-zh/(zk0-zh*(za-zh*zb)))
+         endcase
+         
+      ENDIF 
+   ENDFOR 
+   
+   terre = where(tn GE 1e6)
+   if terre[0] NE -1 then zrhop[terre] = valmask
+
+   if keyword_set(key_performance) THEN print, 'temps rhon', systime(1)-tempsun 
+
+   
+   return, zrhop
+END 
+
+
