;+
;
; @file_comments
; compute depth of ROMS outputs according to ROMS parameters
; stored in the common variable (cm_4mesh) romszinfos
;
; @categories
; Grid
;
; @returns
; the depth of the points (or -1 if error)
;
; @uses
; cm_4mesh
; cm_4data
;
; @restrictions
; common variable (cm_4mesh) romszinfos must be correctly defined
;
; @history
; Sept 2006 Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
FUNCTION romsdepth
;
compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
;
theta_s = romszinfos.theta_s
IF theta_s EQ -1 THEN return, -1
theta_b = romszinfos.theta_b
IF theta_b EQ -1 THEN return, -1
hc = romszinfos.hc
IF hc EQ -1 THEN return, -1
hroms = romszinfos.h
IF hroms[0] EQ -1 THEN return, -1
zeta = romszinfos.zeta
IF zeta[0] EQ -1 THEN return, -1
type = vargrid
;
grille, -1, -1, -1, -1, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz
hroms = hroms[firstx:lastx, firsty:lasty]
nt = n_elements(zeta)/nx/ny
;
cff1 = 1./sinh(theta_s)
cff2 = 0.5/tanh(0.5*theta_s)
;
IF type EQ 'W' THEN BEGIN
sc = (findgen(jpk)-jpk)/jpk
; sc = (dindgen(jpk+1)-jpk)/jpk
; jpk = jpk+1
ENDIF ELSE BEGIN
sc = (findgen(jpk)-jpk-0.5+1)/jpk
ENDELSE
;
cs = (1.-theta_b)*cff1*sinh(theta_s*sc)+theta_b*(cff2*tanh(theta_s*(sc+0.5))-0.5)
cff = hc*(sc-cs)
cff1 = cs
;
hinv = 1./hroms
hinv = hinv[*]#replicate(1., jpk)
; put a z dimension to zeta
zeta = transpose(temporary(zeta))
zeta = reform((temporary(zeta))[*]#replicate(1., jpk), nt, ny, nx, jpk, /overwrite)
zeta = transpose(temporary(zeta), [2, 1, 3, 0])
z0 = replicate(1., nx*ny)#cff[*] + (replicate(1., nx*ny)#cff1[*]) * (hroms[*]#replicate(1., jpk))
z = (z0[*]#replicate(1., nt)) + temporary(zeta) * ((1.+z0*hinv)[*]#replicate(1., nt))
z = reform(z, nx, ny, jpk, nt, /overwrite)
z = -1*reverse(temporary(z), 3)
return, z
end