;+
;
; @file_comments
; calculate the norm of vectors field located on Arakawa C-grid
;
; @categories
; Calculation
;
; @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
; litchamp 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
; litchamp 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
; 't' 'x' 'y' 'z' 'xys' 'xz' 'yz' 'xyz' 'xt' 'yt' 'zt' 'xyt'
; 'xzt' 'yzt' 'xyzt' Direction on which do averages
;
; @keyword _EXTRA
; Used to declare that this routine accepts inherited keywords
;
; @returns
; A 2D (xy), 3D (xyz or yt), 4D (xyzt) Array
;
; @uses
; 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 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.
;
; 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 norm, make sure that time and jpt are defined how
; they have to!
;
; @examples
; To calculate the average of the norm of streams on all the domain
; between 0 and 50:
; IDL> domdef, 0, 50
; IDL> res = norm(un, vn, dir = 'xyz')
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
; 9/6/1999
; @version
; $Id$
;
;-
FUNCTION norm, uu, vv, 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
;---------------------------------------------------------
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 NOT keyword_set(direc) then direc = 0
; construction of u and v at points T
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
;
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)
;----------------------------------------------------------------------------
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
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;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(['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
;------------------------------------------------------------------
; 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(['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
;----------------------------------------------------------------------------
; 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: 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
;----------------------------------------------------------------------------
; 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(['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
;------------------------------------------------------------------
; 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, 'time norm', systime(1)-time1
return, res
end