;+
;
; @file_comments
; calcule la divergence d'un champ 2D
;
; @categories
; matrix
;
; @examples
; IDL> res=divfred(u,v)
;
; @param uu {in}{required}{type=2D array}
; Matrix representing the zonal coordinates (U point)
;
; @param vv {in}{required}{type=2D array}
; Matrix representing the zonal coordinates (V point)
;
; @returns
; the divergence of the input data (with the same size)
;
; @uses
; common
;
; @restrictions
; - Requires SAXO tools
; - les matrices u et v peuvent de 2 a 4 dimensions.
; attention pour distinguer les différentes configurations de u et v
; (xy, xyz, xyt, xyzt), on regarde la variable du common
; -time qui contient le calendrier en jour julien d''IDL auquel
; se rapportent u et v ansi que la variable
; -jpt qui est le nombre de pas de temps à considérer dans time.
; les tableaux u et v sont découpés sur le même domaine
; géographique. A cause du décalage des grilles T, U, V et F il est
; possible que ces 2 tableaux n''aient pas la même taille et se
; reportent à des indices différents. Si tel est le cas les tableaux
; sont redécoupés sur les indices qu'ils ont en commun et le dommaine
; est redéfini pour qu'il colle à ces indices communs.
; pour éviter ces redécoupes utiliser le mot clé /memeindice dans
; domdef
;
; les points sur le bord du dessin sont mis à !values.f_nan
;
; @todo
; ++ pas fini de comprendre, tester (compatibilité saxo), adapter, commenter
; nettoyer
; ++ a comparer et merger avec SAXO_DIR/ToBeReviewed/CALCULS/div.pro
; Written by: fplod 2006-03-27T13:45:14Z aedon.lodyc.jussieu.fr (Darwin)
;
; @history
; fplod 2007-07-30T10:11:33Z aedon.locean-ipsl.upmc.fr (Darwin)
; add header for idldoc
; created from
; /usr/work/sur/fvi/OPA/geomag/divfred.pro
; Guillaume Roullet (grlod\@ipsl.jussieu.fr)
; Creation : printemps 1998
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
; adaptation pour marcher avec un domaine reduit
; 12/1/2000
; @version
; $Id$
;
;-
;
FUNCTION divfred, uu, vv
tempsun = systime(1) ; pour key_performance
@common
u = uu
v = 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
;------------------------------------------------------------
; on trouve les points que u et v ont en communs
;------------------------------------------------------------
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]
case 1 of
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyz
;----------------------------------------------------------------------------
(size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN
;------------------------------------------------------------
; extraction de u et v sur le domaine qui convient
;------------------------------------------------------------
case 1 of
(size(u))[0] NE 3 OR (size(v))[0] NE 3: return, -1
(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 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
(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 :BEGIN
zdiv = -1
GOTO, sortie
end
endcase
;------------------------------------------------------------
; calcul de la divergence
;------------------------------------------------------------
zu = (e2u[indice2d])[*]#replicate(1, nzt)
zu = reform(zu, nx, ny, nzt, /over)
zu = temporary(u)*temporary(zu) $
*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
terreu = where(zu EQ 0)
;;; if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
;
zv = (e1v[indice2d])[*]#replicate(1, nzt)
zv = reform(zv, nx, ny, nzt, /over)
zv = temporary(v)*temporary(zv) $
*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
terrev = where(zv EQ 0)
;;; if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
;
zdiv = 1e6/(e1t[indice2d]*e2t[indice2d])
zdiv = (zdiv)[*]#replicate(1, nzt)
zdiv = reform(zdiv, nx, ny, nzt, /over)
zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0)) $
*tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
;------------------------------------------------------------
; mise a !values.f_nan de la bordure
;------------------------------------------------------------
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
;
zdiv = temporary(zdiv)
;
if n_elements(valmask) EQ 0 THEN valmask = 1e20
terre = where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0)
if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
;------------------------------------------------------------
; pour le trace graphique
;------------------------------------------------------------
vargrid = 'T'
varname = 'div'
varunits = '1e6*s-1'
domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, grille = ['t']
if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan, boite = boite)
END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyt
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN
;------------------------------------------------------------
; extraction de u et v sur le domaine qui convient
;------------------------------------------------------------
case 1 of
(size(u))[0] NE 3 OR (size(v))[0] NE 3: return, -1
(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 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
(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, -1
endcase
;------------------------------------------------------------
; calcul de la divergence
;------------------------------------------------------------
zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
terreu = where(zu EQ 0)
;;; if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
zu = (zu)[*]#replicate(1, jpt)
zu = reform(zu, nx, ny, jpt, /over)
zu = temporary(u)*temporary(zu)
;
zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
terrev = where(zv EQ 0)
;;; if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
zv = (zv)[*]#replicate(1, jpt)
zv = reform(zv, nx, ny, jpt, /over)
zv = temporary(v)*temporary(zv)
;
zdiv = 1e6*tmask[indice2d+jpi*jpj*firstzt]/(e1t[indice2d]*e2t[indice2d])
zdiv = (zdiv)[*]#replicate(1, jpt)
zdiv = reform(zdiv, nx, ny, jpt, /over)
terre = where(zdiv EQ 0)
zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0))
;------------------------------------------------------------
; mise a !values.f_nan de la bordure
;------------------------------------------------------------
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
;
if n_elements(valmask) EQ 0 THEN valmask = 1e20
if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
;------------------------------------------------------------
; pour le trace graphique
;------------------------------------------------------------
vargrid = 'T'
varname = 'div'
varunits = '1e6*s-1'
domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, grille = ['t']
if keyword_set(direc) then zdiv = grossemoyenne(zdiv,direc,/nan, boite = boite)
END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyzt
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
date1 NE date2 AND (size(u))[0] EQ 4:BEGIN
return, report('non code!')
END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xy
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
ELSE : BEGIN ;xy
indice3d = lindgen(jpi, jpj, jpk)
indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt]
;------------------------------------------------------------
; extraction de u et v sur le domaine qui convient
;------------------------------------------------------------
case 1 of
(size(u))[0] NE 2 OR (size(v))[0] NE 2: BEGIN
zdiv = -1
GOTO, sortie
end
(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 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
(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, -1
endcase
;------------------------------------------------------------
; calcul de la divergence
;------------------------------------------------------------
zu = temporary(u)*e2u[indice2d]*(umask())[indice3d]
terreu = where(zu EQ 0)
;;; if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
zv = temporary(v)*e1v[indice2d]*(vmask())[indice3d]
terrev = where(zv EQ 0)
;;; if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
zdiv = zu-shift(zu, 1, 0)+zv-shift(zv, 0, 1)
zdiv = temporary(zdiv)*tmask[indice3d]/(e1t[indice2d]*e2t[indice2d])
;------------------------------------------------------------
; mise a !values.f_nan de la bordure
;------------------------------------------------------------
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
;
zdiv = temporary(zdiv)*1e6
;
if n_elements(valmask) EQ 0 THEN valmask = 1e20
terre = where(tmask[indice3d] EQ 0)
if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
;------------------------------------------------------------
; pour le trace graphique
;------------------------------------------------------------
vargrid = 'T'
varname = 'div'
varunits = '1e6*s-1'
domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, grille = ['t']
if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan, boite = boite)
END
endcase
;------------------------------------------------------------
sortie:
if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun
return, zdiv
end