;+
;
; @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 addtion 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 definig 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)
;
; @returns {type=2D yz array}
; meridional stream function in Sv
;
; @uses
; cm_4mesh
; cm_4data
; grille
; litchamp
; fitintobox
;
; @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
;
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 bondaries must be land point')
IF keyword_set(key_partialstep) AND total(msk[*, ny-1, *]) NE 0 THEN flagdata = 1
;
;
IF keyword_set(key_partialstep) THEN BEGIN
; Rebuild the T-point 3D partial steps array from 2D bottom e3t_ps
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; 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
e3t_3D = replicate(1, nx*ny) # e3t[firstz:lastz]
; apply e3t_ps to e3t_3D at the bottom of the ocean
e3t_3D[temporary(bottom)] = e3t_ps[firstx:lastx, firsty:lasty]
; 3D e3t_3D array
e3t_3D = reform(e3t_3D, nx, ny, nz, /overwrite)
; Rebuild the V-point 3D partial steps array from T-point 3D e3t_3D array
tmp = shift(e3t_3D, 0, -1, 0)
e3v_3D = [ [ (temporary(e3t_3D))[*] ], [ (temporary(tmp))[*] ] ]
e3v_3D = min(temporary(e3v_3D), dimension = 2)
; 3D e1v array
e1_3D = (e1v[firstx:lastx, firsty:lasty])[*] # replicate(1., nz)
; e1v*e3v
e13 = temporary(e1_3D) * temporary(e3v_3D)
ENDIF ELSE BEGIN
e13 = (e1v[firstx:lastx, firsty:lasty])[*] # e3t[firstz:lastz]
ENDELSE
;
; mask the array
vn = 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(e13) * temporary(vn), 1), 2, /cumulative)
; but we force the computation from bottom to up.
msf = total(temporary(vn) * temporary(e13), 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