;+
;
; @file_comments
; compute the barotropic (along the model grid) stream function in Sv using
; the following equation: total(total(e2u*e3u*un),3), 2, /cumulative)/1.e6
;
; @categories
; diagnostics
;
; @param z3d {in} {type=3D xyz array}
; zonal velocity
;
; @keyword NOSTRUCTURE {default=0}{type=scalar: 0 or 1}
; activate if you do not want that msf returns a structure
; but only the array referring to the field.
;
; @keyword REFPOINT {type=2 element vector}
; a 2 element vector [x,y] giving the (approximative) position of the
; point that should be taken as a reference. This position should be
; given with the same coordinates system as the one used in
; glam/gphi. see example
;
; @keyword REFVALUE {type=scalar} {default=0.}
; the bsf value that we want to specify at the position defined by refpoint
;
; @keyword TRANSPORT {type=scalar: 0 or 1} {default=0.}
; activate to specify that z3d is not a zonal current but a zonal
; transport (e2u*e3u*un)
;
; @returns {type=2D xy array}
; barotropic stream function in Sv
;
; @uses
; <pro>cm_4mesh</pro>
; <pro>cm_4data</pro>
; <pro>grille</pro>
; <pro>litchamp</pro>
; <pro>fitintobox</pro>
;
; @restrictions
; - T and U boxzoom parameters must be the same
; - to be valid, mask must be equal to 0 at the bottom and on each
;   side of the domain along x direction
; - define the common variables (of cm_4data)
;       varname = 'BSF'
;       varunit = 'Sv'
;       vargrid = 'F'
;
; @examples
;
;   IDL> initorca05
;   IDL> file = '/Volumes/pim/F31/oce/1y/F31_1y_0040_0040_grid_U.nc'
;   IDL> un = read_ncdf('vozocrtx', 0, 0, /timestep, file = file)
;   IDL> bar = bsf(un, refvalue = 0., refpoint = [0, 20])
;   IDL> plt, bar, -150, 150, int = 10, style = 'so0so', format = '(i4)', lct = 64
;
; @history
; Gurvan Madec, Christian Ethe, Claude Talandier and Sebastien Masson, Dec 2008
;
; @version
; $Id$
;
;-
;
FUNCTION bsf, z3d, NOSTRUCTURE = nostructure, REFPOINT = refpoint, REFVALUE = refvalue, TRANSPORT = transport
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
;
  CASE 1 OF
    firstxt NE firstxu OR lastxt NE lastxu: return, report('T and U box must be defined over the same x indexes, use /memeindices when calling domdef')
    firstyt NE firstyu OR lastyt NE lastyu: return, report('T and U box must be defined over the same y indexes, use /memeindices when calling domdef')
    ELSE:
  ENDCASE
;
  grille, -1, -1, -1, -1, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, grid = 'T'
  umsk = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
  IF total(umsk[*, *, nz-1]) NE 0 THEN return, report('The bottom of the box (defined by lastzt) must be set to land value (=0)')
  IF total(umsk[*, 0, *]) NE 0 THEN return, report('southern boundary must be land point if you want to compute the bsf with this function')
;
  un = litchamp(z3d)
  IF size(un, /n_dimensions) NE 3 THEN return, report('input data must be a 3D array')
  un = fitintobox(temporary(un))
;
; current -> transport
  IF NOT keyword_set(transport) THEN un = temporary(un) * e3u_3d(/e2)
  utr = temporary(umsk) * temporary(un)
; compute the bsf
  bsf = 1.e-6 * total(total(temporary(utr), 3), 2, /cumulative)
; set bsf to 0 in the largest continent... no done...
  IF keyword_set(refpoint) THEN BEGIN
    refind = neighbor(refpoint[0], refpoint[1], glamf[firstx:lastx, firsty:lasty], gphif[firstx:lastx, firsty:lasty], SPHERE = key_onearth)
    IF n_elements(refvalue) EQ 0 THEN refval = - bsf[refind[0]] ELSE refval = refvalue - bsf[refind[0]]
    bsf = temporary(bsf) + refval
  ENDIF
  IF keyword_set(flagdata) THEN bsf[nx-1, *] = !values.f_nan
;
; update data informations
  varname = 'BSF'
  varunit = 'Sv'
  vargrid = 'F'
;
  IF keyword_set(nostructure) THEN return, bsf

  IF keyword_set(key_forgetold) THEN BEGIN
    return, {arr:temporary(bsf), grid:vargrid, unit:varunit, experiment:varexp, name:varname}
  ENDIF ELSE BEGIN
    return, {tab:temporary(bsf), grille:vargrid, unite:varunit, experience:varexp, nom:varname}
  ENDELSE

END