;+
;
; @file_comments
; compute the meridional (along the model grid) stream function in Sv using
; the following equation: total(total(e1v*e3v*vn),1), 2, /cumulative)/1.e6
;
; @categories
; diagnostics
;
; @param z3d {in} {type=3D xyz array}
; meridional velocity
;
; @param mask2d {in} {type=2D xy array} {optional} {default=1 everywhere}
; land/sea (0/1) mask to be applied (in addition to the 3D mask
; of the model) to mask some parts of the domain (and keep for example
; the northend atlantic).
;
; @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 INDEXBOXZOOM
; Set this keyword to a named variable in which msf will return the x
; and y indexes defining the zoom-box that should be used by pltz to do
; the plot (see example)
;
; @keyword MASKOUT
; Set this keyword to a named variable in which msf will return the 2d
; array (yz) should be used by pltz to do plot the land/sea mask(see example)
;
; @keyword TRANSPORT {type=scalar: 0 or 1} {default=0.}
; activate to specify that z3d is not a zonal current but a zonal
; transport (e1v*e3v*vn)
;
; @returns {type=2D yz array}
; meridional 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 V 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 = 'MSF'
;       varunit = 'Sv'
;       vargrid = 'W' ; -> should introduce WF point...
;
; @examples
;
;   IDL> initorca05
;   IDL> file = '/Volumes/pim/F31/oce/1y/F31_1y_0040_0040_grid_V.nc'
;   IDL> vn = read_ncdf('vomecrty', 0, 0, /timestep, file = file)
;   IDL> mm = msf(vn, mskatl, indexboxzoom = ind, maskout = ma)
;   IDL> pltz, mm, -20, 20, int = 1, boxzoom = [ind, 0, 5200], /xindex  $
;            , /yindex, zoom = 1000, maskdta = ma, /no_partial, style = 'so0so'
;
; @history
; Gurvan Madec, Christian Ethe, Claude Talandier and Sebastien Masson, Dec 2008
;
; @version
; $Id$
;
;-
;
FUNCTION msf, z3d, mask2d, INDEXBOXZOOM = indexboxzoom, MASKOUT = maskout, NOSTRUCTURE = nostructure, TRANSPORT = transport
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
;
  CASE 1 OF
    firstxt NE firstxv OR lastxt NE lastxv: return, report('T and V box must be defined over the same x indexes, use /memeindices when calling domdef')
    firstyt NE firstyv OR lastyt NE lastyv: return, report('T and V 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'
;
  vn = litchamp(z3d)
  IF size(vn, /n_dimensions) NE 3 THEN return, report('input data must be a 3D array')
  vn = fitintobox(temporary(vn))
;
  IF n_elements(mask2d) EQ 0 THEN mask2d = replicate(1b, nx, ny)
  IF size(mask2d, /n_dimensions) NE 2 THEN return, report('input mask2d must be a 2D array')
  msk = fitintobox(mask2d)
  msk = (temporary(msk))[*]#replicate(1b, nz)
  msk = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] * temporary(msk)
;
  IF total(msk[*, *, nz-1]) NE 0 THEN return, report('The bottom of the box (defined by lastzt) must be set to land value (=0)')
  IF (total(msk[0, *, *]) NE 0 OR total(msk[nx-1, *, *]) NE 0) AND NOT (keyword_set(key_periodic) AND nx EQ jpi) THEN $
     return, report('for non-periodic domains, eastern and western boundaries must be land point')
  IF keyword_set(key_partialstep) AND total(msk[*, ny-1, *]) NE 0 THEN flagdata = 1
;
; current -> transport
  IF NOT keyword_set(transport) THEN vn = temporary(vn) * e3v_3d(/e1) 
  vtr = temporary(vn) * msk
;
  IF arg_present(maskout) THEN BEGIN
    msk = total(temporary(msk), 1)
    maskout = temporary(msk) NE 0
    maskout = shift(maskout, 0, 1) ; bottom value of W should be taken into account
    maskout[*, 0] = maskout[*, 1]
    maskout = reverse(temporary(maskout), 2)
  ENDIF
;  
; defaut computation:
;  msf = 1.e-6 * total(total(temporary(vtr), 1), 2, /cumulative)
; but we force the computation from bottom to up.
  msf = total(temporary(vtr), 1) ; -> yz array
  msf = reverse(temporary(msf), 2)
  msf = 1.e-6 * total(temporary(msf), 2, /cumulative)
  msf = -1.*reverse(temporary(msf), 2)
  IF keyword_set(flagdata) THEN msf[ny-1, *] = !values.f_nan

; set msf to 0 in the largest continent... no done...

; update data informations
  varname = 'MSF'
  varunit = 'Sv'
  vargrid = 'W' ; -> should introduce WF point...
;
; we must define a domain related to the section
;
; are we using an ORCA grid???

  IF (jpiglo EQ 182 AND jpjglo EQ 149) OR (jpiglo EQ 722 AND jpjglo EQ 511) $
     OR (jpiglo EQ 1442 AND jpjglo EQ 1021) OR (jpiglo EQ 4322 AND jpjglo EQ 3059) THEN BEGIN 
; find the i index with the northest latitude
    ind = (where(gphit[firstx:lastx, firsty:lasty] EQ max(gphit[firstx:lastx, firsty:lasty])))[0] mod nx
  ENDIF ELSE ind = 0
;
  indexboxzoom = [ind, ind, firsty, lasty]
;
  IF keyword_set(nostructure) THEN return, msf

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

END