1 |
function [psi, psimask] = mit_overturning(v,vmask,dx,dz,addlayer); |
2 |
%function [psi, psimask] = mit_overturning(v,vmask,dx,dz); |
3 |
% overturning stream function for mitgcm model, time slabs are handled as |
4 |
% cell objects, integration from bottom to top. |
5 |
|
6 |
if nargin < 5 |
7 |
addlayer = 0; |
8 |
elseif (nargin > 5 & nargin < 4) |
9 |
error('needs 4 or 5 arguments') |
10 |
end |
11 |
[nx ny nz] = size(vmask); |
12 |
for kz=1:nz |
13 |
dxdzs(:,:,kz) = (vmask(:,:,kz).*dx)*dz(kz); |
14 |
end |
15 |
|
16 |
% mask for stream function |
17 |
pmask = change(squeeze(nanmean(vmask)),'~=',NaN,1); |
18 |
% add another psi-point at the bottom of each column |
19 |
for ky=1:ny; |
20 |
iz = min(find(isnan(pmask(ky,:)))); |
21 |
if ~isempty(iz) & iz > 1 |
22 |
pmask(ky,iz) = 1; |
23 |
end |
24 |
end |
25 |
if iscell(v) |
26 |
nt = length(v); |
27 |
psi = cell(size(v)); |
28 |
for k=1:nt |
29 |
vdxdz = v{k}(:,:,:).*dxdzs; %(dxdzs.*vmask); |
30 |
zave = fliplr(change(squeeze(nansum(vdxdz)),'==',NaN,0)); |
31 |
psi{k} = -fliplr(cumsum(zave,2)).*pmask; |
32 |
end |
33 |
else |
34 |
nt = size(v,4); |
35 |
for k=1:nt |
36 |
vdxdz = squeeze(v(:,:,:,k)).*dxdzs; %(dxdzs.*vmask); |
37 |
zave = fliplr(squeeze(sum(change(vdxdz,'==',NaN,0),1))); |
38 |
psi(:,:,k) = -fliplr(cumsum(zave,2)).*pmask; |
39 |
end |
40 |
end |
41 |
% add another layer at the bottom; psi is zero in the layer by definition |
42 |
if addlayer |
43 |
disp('mit_overturning: adding a layer with psi = 0 at the bottom') |
44 |
pmask = [pmask pmask(:,end)]; |
45 |
if iscell(psi) |
46 |
nt = length(psi); |
47 |
for k=1:nt |
48 |
psi{k} = [psi{k} change(psi{k}(:,end),'~=',NaN,0)]; |
49 |
end |
50 |
else |
51 |
psi = cat(2,psi,change(psi(:,end,:),'~=',NaN,0)); |
52 |
end |
53 |
end |
54 |
if nargout == 2 |
55 |
psimask = pmask; |
56 |
end |
57 |
|
58 |
return |
59 |
|
60 |
|
61 |
|
62 |
|
63 |
|