;+
;
; @file_comments
; compute the gradient of a variable located on Arakawa C-grid.
;
; @categories
; Calculation
;
; @param FIELD
; The field for which we want to compute the gradient. A 2D (xy),
; 3D (xyz or yt) or 4D (xyzt) array or a structure readable by
; litchamp and containing a 2D (xy), 3D (xyz or yt) or 4D (xyzt) array.
; Note that the dimension of the array must suit the domain dimension.
;
; @param DIREC {type=scalar string}
; the gradient direction: 'x', 'y', 'z'
;
; @keyword MILLION {default=0}{type=scalar: 0 or 1}
; Activate to multiply returned array by 1.e6
;
; @keyword _EXTRA
; Used to declare that this routine accepts inherited keywords
;
; @returns
; the gradient of the input data with the same size 2D, 3D or 4D
; array, located on the appropriate grid (see restrictions)
;
; @uses
; cm_4cal
; cm_4data
; cm_4mesh
;
; @restrictions
; - Works only for Arakawa C-grid.
; - When computing the gradient, the result is not on the same grid point
; than the input data. In consequence, we update, vargrid and the grid position
; parameters (firstx[tuvf], lastx[tuvf], nx[tuvf], firsty[tuvf], lasty[tuvf],
; ny[tuvf], firstz[tw], lastz[tw], nz[tw]).
; - points that cannot be computed (domain boundaries, coastline) are set to NaN
; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases.
;
; @examples
;
; IDL> \@tst_initorca2
; IDL> plt, grad({arr:gphit,g:'T'}, 'x')
; IDL> plt, grad({arr:gphit,g:'T'}, 'y')
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
FUNCTION grad, field, direc, MILLION=million, _EXTRA=ex
;
compile_opt idl2, strictarrsubs
;
@cm_4cal ; for jpt
@cm_4data ; for varname, vargrid, vardate, varunit, valmask
@cm_4mesh
;
time1 = 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 grad is based on Arakawa C-grid.' $
, 'U and V grids must therefore be defined'])
;
res = litchamp(field)
szres = size(res)
grille, mask, glam, gphi, gdep, nx, ny, nz $
, firstx, firsty, firstz, lastx, lasty, lastz
;
if n_elements(valmask) EQ 0 then valmask = 1.e20
varname = 'grad of '+varname
IF keyword_set(million) THEN varunit = '1.e6*'+varunit+'/m' ELSE varunit = varunit+'/m'
case strupcase(vargrid) of
'T':BEGIN
case direc of
'x':BEGIN
divi = e1u[firstx:lastx, firsty:lasty]
newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
vargrid = 'U'
firstxu = firstxt & lastxu = lastxt & nxu = nxt
firstyu = firstyt & lastyu = lastyt & nyu = nyt
END
'y':BEGIN
divi = e2v[firstx:lastx, firsty:lasty]
newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
vargrid = 'V'
firstxv = firstxt & lastxv = lastxt & nxv = nxt
firstyv = firstyt & lastyv = lastyt & nyv = nyt
END
'z':BEGIN
divi = e3w[firstz:lastz]
newmask = mask
vargrid = 'W'
firstzw = firstzt & lastzw = lastzt & nzw = nzt
END
ELSE:return, report('Bad definition of direction argument')
ENDCASE
END
'W':BEGIN
case direc of
'x':BEGIN
divi = e1u[firstx:lastx, firsty:lasty]
newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
vargrid = 'U'
firstxu = firstxt & lastxu = lastxt & nxu = nxt
firstyu = firstyt & lastyu = lastyt & nyu = nyt
END
'y':BEGIN
divi = e2v[firstx:lastx, firsty:lasty]
newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
vargrid = 'V'
firstxv = firstxt & lastxv = lastxt & nxv = nxt
firstyv = firstyt & lastyv = lastyt & nyv = nyt
END
'z':BEGIN
divi = e3t[firstz:lastz]
newmask = mask
vargrid = 'T'
firstzt = firstzw & lastzt = lastzw & nzt = nzw
END
ELSE:return, report('Bad definition of direction argument')
endcase
END
'U':BEGIN
case direc of
'x':BEGIN
divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty]
newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
vargrid = 'T'
firstxt = firstxu & lastxt = lastxu & nxt = nxu
firstyt = firstyu & lastyt = lastyu & nyt = nyu
END
'y':BEGIN
divi = e2f[firstx:lastx, firsty:lasty]
newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
vargrid = 'F'
firstxf = firstxu & lastxf = lastxu & nxf = nxu
firstyf = firstyu & lastyf = lastyu & nyf = nyu
END
'z':BEGIN
divi = e3w[firstz:lastz]
newmask = mask
vargrid = 'W'
firstzw = firstzt & lastzw = lastzt & nzw = nzt
END
ELSE:return, report('Bad definition of direction argument')
endcase
END
'V':BEGIN
case direc of
'x':BEGIN
divi = e1f[firstx:lastx, firsty:lasty]
newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
vargrid = 'F'
firstxf = firstxv & lastxf = lastxv & nxf = nxv
firstyf = firstyv & lastyf = lastyv & nyf = nyv
END
'y':BEGIN
divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty]
newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
vargrid = 'T'
firstxt = firstxv & lastxt = lastxv & nxt = nxv
firstyt = firstyv & lastyt = lastyv & nyt = nyv
END
'z':BEGIN
divi = e3w[firstz:lastz]
newmask = mask
vargrid = 'W'
firstzw = firstzt & lastzw = lastzt & nzw = nzt
END
ELSE:return, report('Bad definition of direction argument')
endcase
END
'F':BEGIN
; case direc of
; 'x':divi = (shift(e1v, -1, 0))[firstx:lastx, firsty:lasty]
; 'y':divi = (shift(e2u, 0, -1))[firstx:lastx, firsty:lasty]
; 'z':divi = e3w[firstz:lastz]
; ELSE:return, report('Bad definition of direction argument')
; endcase
return, report('F grid: case not coded, please contact SAXO team...')
END
ELSE:return, report('Bad definition of vargrid')
ENDCASE
IF keyword_set(million) THEN divi = temporary(divi)*1.e-6
res = fitintobox(temporary(res))
IF n_elements(res) EQ 1 AND res[0] EQ -1 THEN return, res
case 1 of
;----------------------------------------------------------------------------
;xy
;----------------------------------------------------------------------------
szres[0] EQ 2:BEGIN
land = where((temporary(mask))[*, *, firstz] EQ 0)
if land[0] NE -1 then res[temporary(land)] = !values.f_nan
case direc of
'x':BEGIN
res = (shift(res, -1, 0)-res)/temporary(divi)
if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan
if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0)
END
'y':BEGIN
res = (shift(res, 0, -1)-res)/temporary(divi)
res[*, ny-1] = !values.f_nan
if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1)
END
ELSE:return, report('Bad definition of direction argument for the type of array')
ENDCASE
land = where((temporary(newmask))[*, *, firstz] EQ 0)
if land[0] NE -1 then res[temporary(land)] = valmask
END
;----------------------------------------------------------------------------
;xyt
;----------------------------------------------------------------------------
szres[0] EQ 3 AND jpt NE 1:BEGIN
land = where((temporary(mask))[*, *, firstz] EQ 0, cnt)
if land[0] NE -1 then BEGIN
land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt))
res[temporary(land)] = !values.f_nan
ENDIF
divi = (temporary(divi))[*]#replicate(1., jpt)
case direc of
'x':BEGIN
res = (shift(res, -1, 0, 0)-res)/temporary(divi)
if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0)
END
'y':BEGIN
res = (shift(res, 0, -1, 0)-res)/temporary(divi)
res[*, ny-1, *] = !values.f_nan
if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0)
END
ELSE:return, report('Bad definition of direction argument for the type of array')
ENDCASE
land = where((temporary(newmask))[*, *, firstz] EQ 0, cnt)
if land[0] NE -1 then BEGIN
land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt))
res[temporary(land)] = valmask
ENDIF
END
;----------------------------------------------------------------------------
;xyz
;----------------------------------------------------------------------------
szres[0] EQ 3 AND jpt EQ 1:BEGIN
land = where(mask EQ 0)
if land[0] NE -1 then res[temporary(land)] = !values.f_nan
case direc OF
'x':BEGIN
divi = (temporary(divi))[*]#replicate(1., nz)
res = (shift(res, -1, 0, 0)-res)/temporary(divi)
if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0)
END
'y':BEGIN
divi = (temporary(divi))[*]#replicate(1., nz)
res = (shift(res, 0, -1, 0)-res)/temporary(divi)
res[*, ny-1, *] = !values.f_nan
if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0)
END
'z':BEGIN
divi = replicate(1., nx*ny)#(temporary(divi))[*]
if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, /overwrite)
if vargrid EQ 'W' THEN BEGIN
res = (shift(res, 0, 0, 1)-res)/temporary(divi)
res[*, *, 0] = !values.f_nan
ENDIF ELSE BEGIN
res = (res-shift(res, 0, 0, -1))/temporary(divi)
res[*, *, nz-1] = !values.f_nan
ENDELSE
END
ENDCASE
land = where(temporary(newmask) EQ 0)
if land[0] NE -1 then res[temporary(land)] = valmask
END
;----------------------------------------------------------------------------
;xyzt
;----------------------------------------------------------------------------
szres[0] EQ 4:BEGIN
land = where(temporary(mask) EQ 0, cnt)
if land[0] NE -1 then BEGIN
land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*nz*lindgen(jpt))
res[temporary(land)] = !values.f_nan
ENDIF
case direc OF
'x':BEGIN
divi = (temporary(divi))[*]#replicate(1., nz*jpt)
res = (shift(res, -1, 0, 0, 0)-res)/temporary(divi)
if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *, *] = !values.f_nan
if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0, 0)
END
'y':BEGIN
divi = (temporary(divi))[*]#replicate(1., nz*jpt)
res = (shift(res, 0, -1, 0, 0)-res)/temporary(divi)
res[*, ny-1, *, *] = !values.f_nan
if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0, 0)
END
'z':BEGIN
divi = replicate(1., nx*ny)#(temporary(divi))[*]
divi = (temporary(divi))[*]#replicate(1L, jpt)
if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt, /overwrite)
if vargrid EQ 'W' THEN BEGIN
res = (shift(res, 0, 0, 1, 0)-res)/temporary(divi)
res[*, *, 0, *] = !values.f_nan
ENDIF ELSE BEGIN
res = (res-shift(res, 0, 0, -1, 0))/temporary(divi)
res[*, *, nz-1, *] = !values.f_nan
ENDELSE
END
ENDCASE
land = where(newmask EQ 0, cnt)
if land[0] NE -1 then BEGIN
land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*nz*lindgen(jpt))
res[temporary(land)] = valmask
ENDIF
END
ELSE:return, report('input array must have 2, 3 or 4 dimensions')
ENDCASE
if keyword_set(key_performance) THEN print, 'time curl', systime(1)-time1
return, res
END