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 |