;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
;
; @file_comments 
; calculate the norme 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 norme is made before the possible spatial or 
; temporal average because the average of the norme is not equal to the 
; norme of averages

;
; @categories
; Post traitment 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 domdef done)
;
; @keyword DIREC
; 't' 'x' 'y' 'z' 'xys' 'xz' 'yz' 'xyz' 'xt' 'yt' 'zt' 'xyt'
;       'xzt' 'yzt' 'xyzt' Direction on which do averages 
;
; @returns
; Array to trace with plt, pltz or pltt.
;
; @uses
; common.pro
;
; @restrictions 
; The norme is calculated on points TTo do this calculation, we average 
; field U and Von points T before calculate the norme. 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. 
; 
; When we calculate on a reduce geographic domain, field U and V have not 
; necessarily the same number of point. In this case, we recut U and V to 
; keep only common points. We profit of this to redo a domdef which redefine 
; a geographic domain on which fields U and V are extracted on same points
;
; @restrictions
; 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 
; they have to! 
;
; @examples
; To calculate the average of the norme of streams on all the domain 
; between 0 et 50:
;      IDL> res=norme(un,vn,boxzoom=[0,50],dir='xyz')
;
; @history
; Sebastien Masson (smasson@lodyc.jussieu.fr)
;                       9/6/1999
;
; @version
; $Id$
;
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION norme, composanteu, composantev, BOXZOOM = boxzoom, DIREC = direc, _extra = ex
;---------------------------------------------------------
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
@cm_4cal
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
@updatekwd
  ENDIF
;---------------------------------------------------------
   tempsun = systime(1)         ; 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.' $
                     , '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 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]

   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


;------------------------------------------------------------
; 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
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;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 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
;------------------------------------------------------------------
; 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
;------------------------------------------------------------------
; 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
;----------------------------------------------------------------------------
; attribution of the mask and of logitude 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
; 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
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;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 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
;------------------------------------------------------------------
; 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
;----------------------------------------------------------------------------
; attribution of the mask and of longitude and latitude arrays.
; We recover the complete grid to establish a big mask extent in the four 
; direction to cover pointsfor which a land point has been 
; considerated (make a small drawing)
;----------------------------------------------------------------------------
         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
; 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 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
;------------------------------------------------------------------
; 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
;----------------------------------------------------------------------------
; attribution of the mask and of logitude 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
; 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
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xy
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
      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
;------------------------------------------------------------------
; 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
;------------------------------------------------------------------
; 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
;----------------------------------------------------------------------------
; attribution of the mask and of longitude and latitude arrays.
; We recover the complete grid to establish a big mask extent in the four 
; direction to cover pointsfor which a land point has been 
; considerated (make a small drawing)
;----------------------------------------------------------------------------
         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
; 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
end