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	cnh	1.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