tutorial_global_oce_latlon/diags_matlab/mit_barostream.m

function [psi, psimask] = mit_barostream(u,varargin)
%function [psi, psimask] = mit_barostream(u,umask,dy,dz);
% or
%function [psi, psimask] = mit_barostream(u,gridinformation);
% global barotropic stream function for mitgcm model, time slabs are
% handled as cell objects, integration from north to the south (by
% convention). 

% $Header:  $
% $Name:  $

  if nargin == 2
    g = varargin{1};
    umask = g.umask;
    dy = g.dyg;
    dz = g.dz;
  elseif nargin == 4
    umask = varargin{1};
    dy = varargin{2};
    dz = varargin{3};
  else
    error('need 2 (one of which is the grid structure) or 4 arguments')
  end
  
  [nx ny nz] = size(umask);
  for kz=1:nz
    dydzs(:,:,kz) = (umask(:,:,kz).*dy)*dz(kz);
  end

  % mask for stream function
  pmask = change(squeeze(umask(:,:,1)),'==',NaN,0);
  % add psi-point to the north of all wet points
  pmask(1:nx,2:ny) = pmask(1:nx,2:ny)+pmask(1:nx,1:ny-1);
  pmask = change(pmask,'==',0,NaN);
  pmask = change(pmask,'~=',NaN,1);
  % integrate from the north to the south (by convention), change
  % integration direction by flipping the array ubar (transposed because
  % of MITgcm conventions)
  if iscell(u)
    nt = length(u);
    psi = cell(size(u));
    for k=1:nt
      udxdz = change(u{k}.*dydzs,'==',NaN,0);
      ubar = squeeze(sum(udxdz,3));
      psi{k} = fliplr(cumsum(fliplr(ubar),2)).*pmask;
    end
  else
    nt = size(u,4);
    psi = repmat(NaN,[nx ny nt]);
    udxdz = change(u.*repmat(dydzs,[1 1 1 nt]),'==',NaN,0);
    ubar = squeeze(sum(udxdz,3));
    for kt = 1:nt
      psi(:,:,kt) = fliplr(cumsum(fliplr(squeeze(ubar(:,:,kt))),2)).*pmask;
    end
  end
  if nargout == 2
    psimask = pmask;
  end
    
  return
1	function [psi, psimask] = mit_barostream(u,varargin)
2	%function [psi, psimask] = mit_barostream(u,umask,dy,dz);
3	% or
4	%function [psi, psimask] = mit_barostream(u,gridinformation);
5	% global barotropic stream function for mitgcm model, time slabs are
6	% handled as cell objects, integration from north to the south (by
7	% convention).
8
9	% $Header: $
10	% $Name: $
11
12	if nargin == 2
13	g = varargin{1};
14	umask = g.umask;
15	dy = g.dyg;
16	dz = g.dz;
17	elseif nargin == 4
18	umask = varargin{1};
19	dy = varargin{2};
20	dz = varargin{3};
21	else
22	error('need 2 (one of which is the grid structure) or 4 arguments')
23	end
24
25	[nx ny nz] = size(umask);
26	for kz=1:nz
27	dydzs(:,:,kz) = (umask(:,:,kz).dy)dz(kz);
28	end
29
30	% mask for stream function
31	pmask = change(squeeze(umask(:,:,1)),'==',NaN,0);
32	% add psi-point to the north of all wet points
33	pmask(1:nx,2:ny) = pmask(1:nx,2:ny)+pmask(1:nx,1:ny-1);
34	pmask = change(pmask,'==',0,NaN);
35	pmask = change(pmask,'~=',NaN,1);
36	% integrate from the north to the south (by convention), change
37	% integration direction by flipping the array ubar (transposed because
38	% of MITgcm conventions)
39	if iscell(u)
40	nt = length(u);
41	psi = cell(size(u));
42	for k=1:nt
43	udxdz = change(u{k}.*dydzs,'==',NaN,0);
44	ubar = squeeze(sum(udxdz,3));
45	psi{k} = fliplr(cumsum(fliplr(ubar),2)).*pmask;
46	end
47	else
48	nt = size(u,4);
49	psi = repmat(NaN,[nx ny nt]);
50	udxdz = change(u.*repmat(dydzs,[1 1 1 nt]),'==',NaN,0);
51	ubar = squeeze(sum(udxdz,3));
52	for kt = 1:nt
53	psi(:,:,kt) = fliplr(cumsum(fliplr(squeeze(ubar(:,:,kt))),2)).*pmask;
54	end
55	end
56	if nargout == 2
57	psimask = pmask;
58	end
59
60	return