;+
;
; @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
;
; @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 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 norm 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
;
; 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 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
; 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 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
; 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