;+
;
; @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 speficy at the position defined by refpoint
;
; @returns {type=2D xy array}
; barotropic stream function in Sv
;
; @uses
; cm_4mesh
; cm_4data
; grille
; litchamp
; fitintobox
;
; @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
;
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 boudary 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))
;
IF keyword_set(key_partialstep) THEN BEGIN
; 2D e3u at the bottom of the ocean
; see zgr_zps: e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
e3u_ps = [[[e3t_ps[firstx:lastx, firsty:lasty]]], [[shift(e3t_ps[firstx:lastx, firsty:lasty], -1, 0)]]]
e3u_ps = min(temporary(e3u_ps), dimension = 3)
flagdata = NOT (keyword_set(key_periodic) AND nx EQ jpi) AND total(umsk[nx-1, *, *]) NE 0
; level of the bottom of the ocean
bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
bottom = 0L > ( long(temporary(bottom)) - 1L )
; the bottom of the ocean in 3D index is:
bottom = lindgen(nx*ny) + temporary(bottom)*nx*ny
; 3D e3t array
e33 = replicate(1., nx*ny) # e3t[firstz:lastz]
; apply e3u_ps to e33 at the bottom of the ocean
e33[temporary(bottom)] = e3u_ps
; 3D e2u array
e23d = (e2u[firstx:lastx, firsty:lasty])[*] # replicate(1., nzt)
; e2u*e3u
e23 = temporary(e23d) * temporary(e33)
ENDIF ELSE BEGIN
e23 = (e2u[firstx:lastx, firsty:lasty])[*] # e3t[firstz:lastz]
ENDELSE
;
; mask the array
un = temporary(umsk) * temporary(un)
; compute the bsf
bsf = 1.e-6 * total(total(temporary(un) * temporary(e23), 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