global_ocean.128x64x15/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.

  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	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