tutorial_global_oce_latlon/diags_matlab/mit_overturning.m

function [psi, psimask] = mit_overturning(v,vmask,dx,dz,addlayer);
%function [psi, psimask] = mit_overturning(v,vmask,dx,dz);
% overturning stream function for mitgcm model, time slabs are handled as 
% cell objects, integration from bottom to top.

% $Header:  $
% $Name:  $

  if nargin < 5
    addlayer = 0;
  elseif (nargin > 5 & nargin < 4)
    error('needs 4 or 5 arguments')
  end
  [nx ny nz] = size(vmask);
  for kz=1:nz
    dxdzs(:,:,kz) = (vmask(:,:,kz).*dx)*dz(kz);
  end

  % mask for stream function
  pmask = change(squeeze(nanmean(vmask)),'~=',NaN,1);
  % add another psi-point at the bottom of each column
  for ky=1:ny;
    iz = min(find(isnan(pmask(ky,:))));
    if ~isempty(iz) & iz > 1
      pmask(ky,iz) = 1;
    end
  end
  if iscell(v)
    nt = length(v);
    psi = cell(size(v));
    for k=1:nt
      vdxdz  = v{k}(:,:,:).*dxdzs; %(dxdzs.*vmask);
      zave   = fliplr(change(squeeze(nansum(vdxdz)),'==',NaN,0));
      psi{k} = -fliplr(cumsum(zave,2)).*pmask;
    end
  else
    nt = size(v,4);
    for k=1:nt
      vdxdz = squeeze(v(:,:,:,k)).*dxdzs; %(dxdzs.*vmask);
      zave  = fliplr(squeeze(sum(change(vdxdz,'==',NaN,0),1)));
      psi(:,:,k)   = -fliplr(cumsum(zave,2)).*pmask;
    end
  end
  % add another layer at the bottom; psi is zero in the layer by definition
  if addlayer
    disp('mit_overturning: adding a layer with psi = 0 at the bottom')
    pmask = [pmask pmask(:,end)];
    if iscell(psi)
      nt = length(psi);
      for k=1:nt
        psi{k} = [psi{k} change(psi{k}(:,end),'~=',NaN,0)];
      end
    else
      psi = cat(2,psi,change(psi(:,end,:),'~=',NaN,0));
    end
  end
  if nargout == 2
    psimask = pmask;
  end
    
  return


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