;+
;
; @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
; <pro>cm_4mesh</pro>
; <pro>cm_4data</pro>
;
; @restrictions
; common variable (<pro>cm_4mesh</pro>) 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