Index: trunk/SRC/Computation/norm.pro
===================================================================
--- trunk/SRC/ToBeReviewed/CALCULS/norme.pro	(revision 296)
+++ trunk/SRC/Computation/norm.pro	(revision 314)
@@ -2,24 +2,20 @@
 ;
 ; @file_comments
-; calculate the norm of a field of vectors, then make a possible average.
-;   Comment 1: The field of vector can be, 2d:xy, 3d: xyz or xyt,
-; 4d: xyzt
-;   Comment 2:
-; The calculation of the norm is made before the possible spatial or
-; temporal average because the average of the norm is not equal to the
-; norm of averages
+; calculate the norm of vectors field located on Arakawa C-grid
 ;
 ; @categories
 ; Calculation
 ;
-; @param COMPOSANTEU {in}{required}
-; an 2d, 3d or 4d array
-;
-; @param COMPOSANTEV {in}{required}
-; an 2d, 3d or 4d array
-;
-; @keyword BOXZOOM
-; boxzoom on which do the average (by default the domain selected
-; by the last <pro>domdef</pro> done)
+; @param UU {in}{required}
+; Matrix representing the zonal coordinates (at U/V point) of a field of vectors
+; A 2D (xy), 3D (xyz or yt), 4D (xyzt) or a structure readable by
+; <pro>litchamp</pro> and containing a 2D (xy), 3D (xyz or yt), 4D (xyzt) array.
+; Note that the dimension of the array must suit the domain dimension.
+;
+; @param VV {in}{required}
+; Matrix representing the meridional coordinates (at V/U point) of a field of vectors
+; A 2D (xy), 3D (xyz or yt), 4D (xyzt) or a structure readable by
+; <pro>litchamp</pro> and containing a 2D (xy), 3D (xyz or yt), 4D (xyzt) array.
+; Note that the dimension of the array must suit the domain dimension.
 ;
 ; @keyword DIREC
@@ -27,13 +23,18 @@
 ;       'xzt' 'yzt' 'xyzt' Direction on which do averages
 ;
+; @keyword _EXTRA
+; Used to declare that this routine accepts inherited keywords
+;
 ; @returns
-; Array to trace with <pro>plt</pro>, <pro>pltz</pro> or <pro>pltt</pro>.
+; A 2D (xy), 3D (xyz or yt), 4D (xyzt) Array
 ;
 ; @uses
-; common.pro
+; @cm_4mesh
+; @cm_4data
+; @cm_4cal
 ;
 ; @restrictions
 ; The norm is calculated on points T. To do this calculation, we average
-; field U and V on points T before calculate the norme. At the edge of
+; field U and V on points T before calculate the norm. At the edge of
 ; coast and of domain, we can not calculate fields U and V at points T,
 ; that is why these points are at value !values.f_nan.
@@ -46,16 +47,16 @@
 ; To know what type of array we work with, we  test its size and dates
 ; gave by time[0] and time[jpt-1] to know if thee is a temporal dimension.
-; Before to start norme, make sure that time and jpt are defined how
+; Before to start norm, make sure that time and jpt are defined how
 ; they have to!
 ;
 ; @examples
-; To calculate the average of the norme of streams on all the domain
+; To calculate the average of the norm of streams on all the domain
 ; between 0 and 50:
-;      IDL> res=norme(un,vn,boxzoom=[0,50],dir='xyz')
+;      IDL> domdef, 0, 50
+;      IDL> res = norm(un, vn, dir = 'xyz')
 ;
 ; @history
 ; Sebastien Masson (smasson\@lodyc.jussieu.fr)
 ;                       9/6/1999
-;
 ; @version
 ; $Id$
@@ -63,5 +64,5 @@
 ;-
 ;
-FUNCTION norme, composanteu, composantev, BOXZOOM = boxzoom, DIREC = direc, _EXTRA = ex
+FUNCTION norm, uu, vv, DIREC = direc, _EXTRA = ex
 ;
   compile_opt idl2, strictarrsubs
@@ -75,65 +76,52 @@
   ENDIF
 ;---------------------------------------------------------
-   tempsun = systime(1)         ; To 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.' $
+  time1 = systime(1)          ; To key_performance
+;
+  IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
+     return, report(['This version of norm 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
+;------------------------------------------------------------
+  if NOT keyword_set(direc) then direc = 0
 ; construction of u and v at points T
-   u = litchamp(composanteu)
-   v = litchamp(composantev)
-   date1 = time[0]
-   if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
+  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
-
-   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
-
-
+  if (size(u))[0] NE (size(v))[0] then return,  -1
+;
+  grilleu = litchamp(uu, /grid)
+  if grilleu EQ '' then grilleu = 'U'
+  grillev = litchamp(vv, /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('Case not coded, but easy to do...!')
+  ENDIF ELSE interpolle = 1
+  if keyword_set(inverse) then begin
+    tmp = u
+    u = temporary(v)
+    v = temporary(tmp)
+  endif
 ;------------------------------------------------------------
 ; We find common points between u and v
 ;------------------------------------------------------------
-   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
+  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)
+;----------------------------------------------------------------------------
+  vargrid = 'T'
+  varname = 'norm'
+  if n_elements(valmask) EQ 0 THEN valmask = 1e20
+  firstxt = indicex[0] & lastxt = indicex[0]+nx-1 & nxt = nx
+  firstyt = indicey[0] & lastyt = indicey[0]+ny-1 & nyt = ny
+;----------------------------------------------------------------------------
+  case 1 of
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
@@ -141,84 +129,88 @@
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
-      (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]
+    (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 of u and v on the appropriated domain
 ;------------------------------------------------------------
-         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
+      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(['the third dimension of u (' + strtrim((size(u))[3], 1) $
+                                  +') must be equal to nzt (' + strtrim(nzt, 1) + ') or jpk ('+strtrim(jpk, 1)+')'])
+          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(['We cannot find re-cut on common indexes (see restriction un the header.).' $
+                     , 'To avoid this problem, read the full domain' $
+                     , 'or use the keyword /memeindice in domdef when defining the zoom area'])
+      endcase
 ;------------------------------------------------------------------
 ; We reshape u and v to make sure that no dimension has been erased
 ;------------------------------------------------------------------
-         if nzt EQ 1 then begin
-            u = reform(u, nx, ny, nzt, /over)
-            v = reform(v, nx, ny, nzt, /over)
-         endif
+      if nzt EQ 1 then begin
+        u = reform(u, nx, ny, nzt, /over)
+        v = reform(v, nx, ny, nzt, /over)
+      endif
 ;------------------------------------------------------------------
 ; construction of u and v at points 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
+      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 of the mask and of longitude and latitude arrays
 ;----------------------------------------------------------------------------
-         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
+      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
 ; All kind of average
-         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
+      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
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
@@ -226,38 +218,41 @@
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
-      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]
+    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 of u and v on the appropriated domain
 ;------------------------------------------------------------
-         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
+      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(['We cannot find re-cut on common indexes (see restriction un the header.).' $
+                     , 'To avoid this problem, read the full domain' $
+                     , 'or use the keyword /memeindice in domdef when defining the zoom area'])
+      endcase
 ;------------------------------------------------------------------
 ; construction of u and v at points 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
+      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 of the mask and of longitude and latitude arrays.
@@ -266,29 +261,29 @@
 ; considerated (make a small drawing)
 ;----------------------------------------------------------------------------
-         mask = tmask[indice2d+jpi*jpj*firstzt]
-         if ny EQ 1 then mask = reform(mask, nx, ny, /over)
+      mask = tmask[indice2d+jpi*jpj*firstzt]
+      if ny EQ 1 then mask = reform(mask, nx, ny, /over)
 ;-----------------------------------------------------------
 ; construction of land containing all points to mask
 ;-----------------------------------------------------------
-         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
+      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
+      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
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
@@ -296,81 +291,85 @@
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
-      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]
+    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 of u and v on the appropriated domain
 ;------------------------------------------------------------
-         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
+      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: report, (['the third dimension of u (' + strtrim((size(u))[3], 1) $
+                            +') must be equal to nzt (' + strtrim(nzt, 1) + ') or jpk ('+strtrim(jpk, 1)+')'])
+          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(['We cannot find re-cut on common indexes (see restriction un the header.).' $
+                     , 'To avoid this problem, read the full domain' $
+                     , 'or use the keyword /memeindice in domdef when defining the zoom area'])
+      endcase
 ;------------------------------------------------------------------
 ; construction of u and v at points 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
+      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 of the mask and of longitude and latitude arrays
 ;----------------------------------------------------------------------------
-         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
+      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
 ; All kind of average
-         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
+      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
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
@@ -378,45 +377,48 @@
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
-      ELSE:BEGIN                ;xy
-         indice2d = lindgen(jpi, jpj)
-         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
+    ELSE:BEGIN                  ;xy
+      indice2d = lindgen(jpi, jpj)
+      indice2d = indice2d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1]
 ;------------------------------------------------------------
 ; extraction of u and v on the appropriated domain
 ;------------------------------------------------------------
-         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
+      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(['We cannot find re-cut on common indexes (see restriction un the header.).' $
+                     , 'To avoid this problem, read the full domain' $
+                     , 'or use the keyword /memeindice in domdef when defining the zoom area'])
+      endcase
 ;------------------------------------------------------------------
 ; We reshape u and v to make sure that no dimension has been erased
 ;------------------------------------------------------------------
-         if ny EQ 1 then begin
-            u = reform(u, nx, ny, /over)
-            v = reform(v, nx, ny, /over)
-         endif
+      if ny EQ 1 then begin
+        u = reform(u, nx, ny, /over)
+        v = reform(v, nx, ny, /over)
+      endif
 ;------------------------------------------------------------------
 ; construction of u and v at points 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
+      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 of the mask and of longitude and latitude arrays.
@@ -425,27 +427,27 @@
 ; considerated (make a small drawing)
 ;----------------------------------------------------------------------------
-         mask = tmask[indice2d+jpi*jpj*firstzt]
-         if nyt EQ 1 THEN mask = reform(mask, nx, ny, /over)
+      mask = tmask[indice2d+jpi*jpj*firstzt]
+      if nyt EQ 1 THEN mask = reform(mask, nx, ny, /over)
 ;-----------------------------------------------------------
 ; construction of land containing all points to mask
 ;-----------------------------------------------------------
-         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
+      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
 ; All kind of average
-         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
+      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, 'time norm', systime(1)-time1
+  return, res
 end
