;+
;
; @file_comments
; compute the 3D e3t from e3t and e3t_ps
;
; @categories
; diagnostics
;
; @keyword e1 {default=0}{type=scalar: 0 or 1}
; activate to compute e1t * e3t instead of e3t
; note that of both keywords e1 and e2 are activated we compute e1t * e2t * e3t instead of e3t
;
; @keyword e2 {default=0}{type=scalar: 0 or 1}
; activate to compute e2t * e3t instead of e3t
; note that of both keywords e1 and e2 are activated we compute e1t * e2t * e3t instead of e3t
;
; @returns
; e3t
;
; @uses
; cm_4mesh.pro
;
; @history
; Sebastien Masson, Jan 2011
;
; @version
; $Id$
;
;-
function e3t_3d, e1 = e1, e2 = e2, fstx = fstx, lstx = lstx, fsty = fsty, lsty = lsty
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
;
  IF n_elements(fstx) EQ 0 THEN fstx = firstxt
  IF n_elements(lstx) EQ 0 THEN lstx = lastxt
  IF n_elements(fsty) EQ 0 THEN fsty = firstyt
  IF n_elements(lsty) EQ 0 THEN lsty = lastyt
  nx = lstx - fstx + 1
  ny = lsty - fsty + 1
;
  CASE 1 OF
    keyword_set(e1) AND keyword_set(e2):arr2d = e1t[fstx:lstx, fsty:lsty] * e2t[fstx:lstx, fsty:lsty]
    keyword_set(e1)                    :arr2d = e1t[fstx:lstx, fsty:lsty]
                        keyword_set(e2):arr2d =                             e2t[fstx:lstx, fsty:lsty]
    ELSE:arr2d = replicate(1.d, nx*ny)
  ENDCASE
  e3t3d = arr2d[*] # e3t[firstzt:lastzt]
  e3t3d = reform(e3t3d, nx, ny, nzt, /overwrite)
;
  IF keyword_set(key_partialstep) THEN BEGIN
; Change bottom values with e3t_ps
; level of the bottom of the ocean
    bottom = total(tmask[fstx:lstx, fsty:lsty, *], 3)
    sea = where(bottom NE 0)
    bottom2 = long(temporary(bottom)) - firstzt
    bottom  = bottom2 - 1L
    ok  = inter(sea, where(bottom  GE 0 AND bottom  LT nzt))
    ok2 = inter(sea, where(bottom2 GE 0 AND bottom2 LT nzt))
; apply e3t_ps to e3t_3D at the bottom of the ocean
    IF ok[0]  NE -1 THEN BEGIN
; the bottom of the ocean in 3D index is:
      bottom  = (lindgen(nx*ny))[ok ] + nx*ny*(temporary(bottom ))[ok ]
      e3t3d[temporary(bottom )] = arr2d[ok ] * (e3t_ps[fstx:lstx, fsty:lsty])[ok ]
    ENDIF
    IF ok2[0] NE -1 THEN BEGIN
      bottom2 = (lindgen(nx*ny))[ok2] + nx*ny*(temporary(bottom2))[ok2]
      e3t3d[temporary(bottom2)] = arr2d[ok2] * (e3t_ps[fstx:lstx, fsty:lsty])[ok2]
    ENDIF
;
  ENDIF
;
  return, e3t3d
END