MITgcm_contrib/mitgcm_tools/overturning.m

function [Psi] = overturning(GRID,v,varargin);
%psi=overturning(GRID,v);
%psi=overturning(GRID,v,mask);
%
%Calculates the volume overturning stream function (m^3/s).
%e.g.
% >> GRID=loadgrid;
% >> STATE=loadstate;
% >> psi=overturning(GRID,STATE.V,GRID.matl);
% >> contourf(GRID.yg(1,:),GRID.rf,psi'/1e6);colorbar
%
%Written by adcroft@mit.edu, 2002
%$Header:

if nargin==3
 mapmsk=varargin{:};
elseif nargin==2
 mapmsk=max( GRID.mskc(:,:,1), GRID.mskc(:,:,end) );
else
 error('Wrong number of arguments')
end
mapmsk=mapmsk.*mapmsk(:,[1 1:end-1]);

N=size(GRID.hfacs);
nxy=prod(size(GRID.rac));
nr=prod(size(GRID.drf));

DRF=spdiags(GRID.drf',0,nr,nr);
dz=reshape(GRID.hfacs,[nxy nr])*DRF;
area=dz.*( (GRID.dxg(:).*mapmsk(:))*ones(1,nr) );
area=reshape(area, N);

V=squeeze( sum(v.*area,1) );

Psi=zeros(N(2),nr+1);
for k=nr:-1:1;
 Psi(:,k)=Psi(:,k+1)-V(:,k);
end

vmsk=squeeze( sum(GRID.hfacs.*area,1) );
vmsk( find(vmsk~=0) )=1;
pmsk=vmsk(:,[1 1:end]).*vmsk(:,[1:end end]);

pmsk( find(pmsk==0) )=NaN;
Psi=Psi.*pmsk(:,[1 1:end-1]);
1	function [Psi] = overturning(GRID,v,varargin);
2	%psi=overturning(GRID,v);
3	%psi=overturning(GRID,v,mask);
4	%
5	%Calculates the volume overturning stream function (m^3/s).
6	%e.g.
7	% >> GRID=loadgrid;
8	% >> STATE=loadstate;
9	% >> psi=overturning(GRID,STATE.V,GRID.matl);
10	% >> contourf(GRID.yg(1,:),GRID.rf,psi'/1e6);colorbar
11	%
12	%Written by adcroft@mit.edu, 2002
13	%$Header:
14
15	if nargin==3
16	mapmsk=varargin{:};
17	elseif nargin==2
18	mapmsk=max( GRID.mskc(:,:,1), GRID.mskc(:,:,end) );
19	else
20	error('Wrong number of arguments')
21	end
22	mapmsk=mapmsk.*mapmsk(:,[1 1:end-1]);
23
24	N=size(GRID.hfacs);
25	nxy=prod(size(GRID.rac));
26	nr=prod(size(GRID.drf));
27
28	DRF=spdiags(GRID.drf',0,nr,nr);
29	dz=reshape(GRID.hfacs,[nxy nr])*DRF;
30	area=dz.( (GRID.dxg(:).mapmsk(:))*ones(1,nr) );
31	area=reshape(area, N);
32
33	V=squeeze( sum(v.*area,1) );
34
35	Psi=zeros(N(2),nr+1);
36	for k=nr:-1:1;
37	Psi(:,k)=Psi(:,k+1)-V(:,k);
38	end
39
40	vmsk=squeeze( sum(GRID.hfacs.*area,1) );
41	vmsk( find(vmsk~=0) )=1;
42	pmsk=vmsk(:,[1 1:end]).*vmsk(:,[1:end end]);
43
44	pmsk( find(pmsk==0) )=NaN;
45	Psi=Psi.*pmsk(:,[1 1:end-1]);