;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; 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:
;
;       BOITE: boite 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,boite=[0,50],dir='xyz')
;
; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
;                       9/6/1999
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION norme,composanteu, composantev, BOITE=boite,DIREC=direc
@common
   tempsun = systime(1)         ; pour key_performance
;------------------------------------------------------------
   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))[premierxu:premierxu+nxu-1]
   indicexv = (lindgen(jpi))[premierxv:premierxv+nxv-1]
   indicex = inter(indicexu, indicexv)
   indiceyu = (lindgen(jpj))[premieryu:premieryu+nyu-1]
   indiceyv = (lindgen(jpj))[premieryv:premieryv+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,premierzt:dernierzt]
;------------------------------------------------------------
; 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 premierxu then u = u[0:nx-1,*,*] ELSE u = u[1: nx, *,*] 
                     IF nxv NE nx THEN $
                      if indicex[0] EQ premierxv then v = v[0:nx-1,*,*] ELSE v = v[1: nx, *,*]
                     IF nyu NE ny THEN $
                      if indicey[0] EQ premieryu then u = u[*,0:ny-1,*] ELSE u = u[*, 1: ny,*] 
                     IF nyv NE ny THEN $
                      if indicey[0] EQ premieryv then v = v[*,0:ny-1,*] ELSE v = v[*, 1: ny,*]
                  end
                  jpk:BEGIN 
                     if nxu NE nx then $
                      if indicex[0] EQ premierxu then u = u[0:nx-1, *,premierzt:dernierzt] ELSE u = u[1: nx, *,premierzt:dernierzt] 
                     IF nxv NE nx THEN $
                      if indicex[0] EQ premierxv then v = v[0:nx-1, *,premierzt:dernierzt] ELSE v = v[1: nx, *,premierzt:dernierzt]
                     IF nyu NE ny THEN $
                      if indicey[0] EQ premieryu then u = u[*, 0:ny-1,premierzt:dernierzt] ELSE u = u[*, 1: ny,premierzt:dernierzt] 
                     IF nyv NE ny THEN $
                      if indicey[0] EQ premieryv then v = v[*, 0:ny-1,premierzt:dernierzt] ELSE v = v[*, 1: ny,premierzt:dernierzt]
                  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_periodique) OR nx NE jpi then u(0,*,*)=a
         a=v(*,0,*)
         v=(v+shift(v,0,1,0))/2.
         if NOT keyword_set(key_periodique) 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_periodique) 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
         oldboite = [lon1, lon2, lat1, lat2, prof1, prof2]
         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme
         if keyword_set(direc) then res = moyenne(res,direc,/nan, boite = boite)
;-----------------------------------------------------------
      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 premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
               IF nxv NE nx THEN $
                if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
               IF nyu NE ny THEN $
                if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
               IF nyv NE ny THEN $
                if indicey[0] EQ premieryv 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_periodique) OR nx NE jpi then u(0,*,*)=a
         a=v(*,0,*)
         v=(v+shift(v,0,1,0))/2.
         if NOT keyword_set(key_periodique) 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*(niveau-1)]
         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_periodique) 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], prof1, prof2, /meme
         if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boite = boite)
      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,premierzt:dernierzt]
;------------------------------------------------------------
; 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 premierxu then u = u[0:nx-1,*,*,*] ELSE u = u[1: nx, *,*,*] 
                     IF nxv NE nx THEN $
                      if indicex[0] EQ premierxv then v = v[0:nx-1,*,*,*] ELSE v = v[1: nx, *,*,*]
                     IF nyu NE ny THEN $
                      if indicey[0] EQ premieryu then u = u[*,0:ny-1,*,*] ELSE u = u[*, 1: ny,*,*] 
                     IF nyv NE ny THEN $
                      if indicey[0] EQ premieryv then v = v[*,0:ny-1,*,*] ELSE v = v[*, 1: ny,*,*]
                  end
                  jpk:BEGIN 
                     if nxu NE nx then $
                      if indicex[0] EQ premierxu then u = u[0:nx-1, *,premierzt:dernierzt,*] ELSE u = u[1: nx, *,premierzt:dernierzt,*] 
                     IF nxv NE nx THEN $
                      if indicex[0] EQ premierxv then v = v[0:nx-1, *,premierzt:dernierzt,*] ELSE v = v[1: nx, *,premierzt:dernierzt,*]
                     IF nyu NE ny THEN $
                      if indicey[0] EQ premieryu then u = u[*, 0:ny-1,premierzt:dernierzt,*] ELSE u = u[*, 1: ny,premierzt:dernierzt,*] 
                     IF nyv NE ny THEN $
                      if indicey[0] EQ premieryv then v = v[*, 0:ny-1,premierzt:dernierzt,*] ELSE v = v[*, 1: ny,premierzt:dernierzt,*]
                  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,premierzt:dernierzt, *]
               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,premierzt:dernierzt, *]
            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_periodique) 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_periodique) 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_periodique) 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], prof1, prof2, /meme
         if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boite = boite)
      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 premierxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
               IF nxv NE nx THEN $
                if indicex[0] EQ premierxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
               IF nyu NE ny THEN $
                if indicey[0] EQ premieryu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
               IF nyv NE ny THEN $
                if indicey[0] EQ premieryv 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_periodique) OR nx NE jpi then u(0,*)=a
         a=v(*,0)
         v=(v+shift(v,0,1))/2.
         if NOT keyword_set(key_periodique) 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*(niveau-1)]
         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_periodique) 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
         oldboite = [lon1, lon2, lat1, lat2, prof1, prof2]
         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme
         if keyword_set(direc) then res = moyenne(res,direc,/nan, boite = boite)
      END
;----------------------------------------------------------------------------
   endcase
;------------------------------------------------------------
   if keyword_set(key_performance) THEN print, 'temps norme', systime(1)-tempsun 
   return, res
end