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