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