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	jmc	1.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	jmc	1.3	% $Header: $
10	jmc	1.1	% $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