source: trunk/SRC/Grid/romsdepth.pro @ 504

;+
;
; @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, new_s_coord = new_s_coord
;
  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[firstx:lastx, firsty:lasty, *])/nx/ny
;
  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
;
  IF keyword_set(new_s_coord) THEN BEGIN

    IF theta_s GT 0. THEN BEGIN
      csrf = ( 1. - cosh(theta_s*sc) ) / ( cosh(theta_s) - 1. )
    ENDIF ELSE BEGIN
      csrf = - sc*sc
    ENDELSE 
    if theta_b GT 0. THEN BEGIN 
      cs = ( exp(theta_b*csrf) - 1. ) / ( 1. - exp(-theta_b) )  
    ENDIF ELSE BEGIN
      cs = csrf
    ENDELSE
    
    hh = (hroms[*]#replicate(1., jpk))
    z0 = hh * ( replicate(1., nx*ny)#(hc*sc)[*] + hh * replicate(1., nx*ny)#cs[*] ) / ( hc + hh )
    
    
  ENDIF ELSE BEGIN
    
    cff1 = 1./sinh(theta_s)
    cff2 = 0.5/tanh(0.5*theta_s)

    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

    z0 = replicate(1., nx*ny)#cff[*] + (replicate(1., nx*ny)#cff1[*]) * (hroms[*]#replicate(1., jpk))

  ENDELSE

;
  hinv = 1./hroms
  hinv = hinv[*]#replicate(1., jpk)
; put a z dimension to zeta
  zeta = transpose((temporary(zeta))[firstx:lastx, firsty:lasty, *])
  zeta = reform((temporary(zeta))[*]#replicate(1., jpk), nt, ny, nx, jpk, /overwrite)
  zeta = transpose(temporary(zeta), [2, 1, 3, 0])

  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