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	adcroft	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]);