;+
;
; @file_comments
; compute the horizontal divergence of a vectors field
;
; @categories
; Calculation
;
; @param UU
; Matrix representing the zonal coordinates (U point) of a field of vectors
; A 2D (xy), 3D (xyz or yt) or a structure readable by litchamp and containing
; a 2D (xy), 3D (xyz or yt) array (4D case is not coded yet).
; note that the dimension of the array must suit the domain dimension.
;
; @param VV
; Matrix representing the meridional coordinates (V point) of a field of vectors
; A 2D (xy), 3D (xyz or yt) or a structure readable by litchamp and containing
; a 2D (xy), 3D (xyz or yt) array (4D case is not coded yet).
; note that the dimension of the array must suit the domain dimension.
;
; @keyword DIREC {type=scalar string}
; Use if you want to call moyenne or
; grossemoyenne after the div computation
; (stupid ?) with a mean done in the DIREC direction
;
; @returns
; the divergence of the input data (with the same size)
;
; @uses
; cm_4cal
; cm_4data
; cm_4mmesh
;
; @restrictions
;
; - Works only for Arakawa C-grid.
; - UU must be on U grid, VV must be on V grid
; - 4D case is not coded yet
; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases.
; - U and V arrays are cut in the same geographic domain. Because of the shift between
; T, U, V and F grids, it is possible that these two arrays do not have the same
; size and refer to different indexes. In this case, arrays are re-cut on
; common indexes. To avoid these re-cuts, use the keyword /memeindice in
; domdef
; - When computing the divergence, we update, vargrid, varname, varunits and the
; grid position parameters (firstxt, lastxt, nxt, firstyt, lastyt, nyt).
; - points that cannot be computed (domain boundaries, coastline) are set to NaN
;
; @examples
; IDL> \@tst_initorca2
; IDL> plt, div(dist(jpi,jpj), dist(jpi,jpj))
;
; @history
; Guillaume Roullet (grlod\@ipsl.jussieu.fr): creation; spring 1998
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
; adaptation to work with a reduce domain; 12/1/2000
;
; @version
; $Id$
;
; @todo
; code the 4D case
;
;-
;
FUNCTION div, uu, vv, DIREC = direc
;
compile_opt idl2, strictarrsubs
;
@cm_4cal ; for jpt
@cm_4data ; for varname, vargrid, vardate, varunit, valmask
@cm_4mesh
;
tempsun = systime(1) ; For key_performance
;
IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
return, report(['This version of div is based on Arakawa C-grid.' $
, 'U and V grids must therefore be defined'])
;
u = litchamp(uu)
v = litchamp(vv)
;
szu = size(u)
szv = size(v)
if szu[0] NE szv[0] then return, report('U and V input data must have the same number of dimensions')
;------------------------------------------------------------
; 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)
indice2d = lindgen(jpi, jpj)
indice2d = indice2d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1]
;----------------------------------------------------------------------------
vargrid = 'T'
varname = 'div'
varunits = '1.e6*s-1'
if n_elements(valmask) EQ 0 THEN valmask = 1.e20
firstxt = indicex[0] & lastxt = indicex[0]+nx-1 & nxt = nx
firstyt = indicey[0] & lastyt = indicey[0]+ny-1 & nyt = ny
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
case 1 of
;----------------------------------------------------------------------------
;xyz
;----------------------------------------------------------------------------
szu[0] EQ 3 AND jpt EQ 1:BEGIN
;------------------------------------------------------------
; extraction of U and V on the appropriated domain
;------------------------------------------------------------
case 1 of
szu[1] EQ nxu AND szu[2] EQ nyu AND $
szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN
case 1 of
nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
ELSE :
endcase
END
szu[1] EQ jpi AND szu[2] EQ jpj AND $
szv[1] EQ jpi AND szv[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, -1
endcase
;------------------------------------------------------------
; divergence computation
;------------------------------------------------------------
zu = (e2u[indice2d])[*]#replicate(1., nzt)
landu = where((umask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0)
if landu[0] NE -1 then zu[temporary(landu)] = !values.f_nan
zu = temporary(u) * temporary(zu)
;
zv = (e1v[indice2d])[*]#replicate(1., nzt)
landv = where((vmask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0)
if landv[0] NE -1 then zv[temporary(landv)] = !values.f_nan
zv = temporary(v) * temporary(zv)
;
zdiv = (e1t[indice2d]*e2t[indice2d])[*]#replicate(1.e6, nzt)
zdiv = ( zu - shift(zu, 1, 0, 0) + zv - shift(zv, 0, 1, 0) ) * temporary(zdiv)
;------------------------------------------------------------
; Edging put at !values.f_nan
;------------------------------------------------------------
if NOT keyword_set(key_periodic) OR nx NE jpi then begin
zdiv[0, *, *] = !values.f_nan
zdiv[nx-1, *, *] = !values.f_nan
endif
zdiv[*, 0, *] = !values.f_nan
zdiv[*, ny-1, *] = !values.f_nan
;
land = where(tmask[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0)
if land[0] NE -1 then zdiv[temporary(land)] = valmask
if keyword_set(direc) then zdiv = moyenne(zdiv, direc, /nan)
END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyt
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
szu[0] EQ 3 AND jpt GT 1:BEGIN
;------------------------------------------------------------
; extraction of U and V on the appropriated domain
;------------------------------------------------------------
case 1 of
szu[1] EQ nxu AND szu[2] EQ nyu AND $
szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN
case 1 of
nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
ELSE :
endcase
END
szu[1] EQ jpi AND szu[2] EQ jpj AND $
szv[1] EQ jpi AND szv[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, -1
endcase
;------------------------------------------------------------
; divergence computation
;------------------------------------------------------------
zu = e2u[indice2d]
landu = where((umask())[indice2d+jpi*jpj*firstzt] EQ 0)
if landu[0] NE -1 then zu[temporary(landu)] = !values.f_nan
zu = (temporary(zu))[*]#replicate(1., jpt)
zu = temporary(u) * temporary(zu)
;
zv = e1v[indice2d]
landv = where((vmask())[indice2d+jpi*jpj*firstzt] EQ 0)
if landv[0] NE -1 then zv[temporary(landv)] = !values.f_nan
zv = (temporary(zv))[*]#replicate(1., jpt)
zv = temporary(v) * temporary(zv)
;
zdiv = (e1t[indice2d]*e2t[indice2d])[*]#replicate(1.e6, jpt)
zdiv = ( zu - shift(zu, 1, 0, 0) + zv - shift(zv, 0, 1, 0) ) * temporary(zdiv)
;------------------------------------------------------------
; Edging put at !values.f_nan
;------------------------------------------------------------
if NOT keyword_set(key_periodic) OR nx NE jpi then begin
zdiv[0, *, *] = !values.f_nan
zdiv[nx-1, *, *] = !values.f_nan
endif
zdiv[*, 0, *] = !values.f_nan
zdiv[*, ny-1, *] = !values.f_nan
;
land = where(tmask[indice2d+jpi*jpj*firstzt] EQ 0, cnt)
if land[0] NE -1 then BEGIN
land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt))
zdiv[temporary(land)] = valmask
ENDIF
if keyword_set(direc) then zdiv = grossemoyenne(zdiv, direc, /nan)
END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyzt
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
szu[0] EQ 4:BEGIN
return, report('Case not coded contact saxo team or make a do loop!')
END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xy
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
szu[0] EQ 2:BEGIN
;------------------------------------------------------------
; extraction of U and V on the appropriated domain
;------------------------------------------------------------
case 1 of
szu[1] EQ nxu AND szu[2] EQ nyu AND $
szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN
case 1 of
nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]
nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]
nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
ELSE :
endcase
END
szu[1] EQ jpi AND szu[2] EQ jpj AND $
szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN
u = u[indice2d]
v = v[indice2d]
END
ELSE:return, -1
endcase
;------------------------------------------------------------
; divergence computation
;------------------------------------------------------------
zu = e2u[indice2d]
landu = where((umask())[indice2d+jpi*jpj*firstzt] EQ 0)
if landu[0] NE -1 then zu[temporary(landu)] = !values.f_nan
zu = temporary(u) * temporary(zu)
zv = e1v[indice2d]
landv = where((vmask())[indice2d+jpi*jpj*firstzt] EQ 0)
if landv[0] NE -1 then zv[temporary(landv)] = !values.f_nan
zv = temporary(v) * temporary(zv)
zdiv = 1.e6 / (e1t[indice2d]*e2t[indice2d])
zdiv = ( zu - shift(zu, 1, 0) + zv - shift(zv, 0, 1) ) * temporary(zdiv)
;------------------------------------------------------------
; Edging put at !values.f_nan
;------------------------------------------------------------
if NOT keyword_set(key_periodic) OR nx NE jpi then begin
zdiv[0, *] = !values.f_nan
zdiv[nx-1, *] = !values.f_nan
endif
zdiv[*, 0] = !values.f_nan
zdiv[*, ny-1] = !values.f_nan
;
land = where(tmask[indice2d+jpi*jpj*firstzt] EQ 0)
if land[0] NE -1 then zdiv[temporary(land)] = valmask
if keyword_set(direc) then zdiv = moyenne(zdiv, direc, /nan)
END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
ELSE:return, report('U and V input arrays must have 2, 3 or 4 dimensions')
ENDCASE
;------------------------------------------------------------
if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun
return, zdiv
end