function for calculating mass flux on isentropic surfaces

flux2dIsen.m

function flux2dIsen = calIsenFlux(vcomp3d,PT3d,pfull,isenCoord,deltaP)
% calculate mass flux in isentropic coordinate
%----INPUT----
% vcomp3d:  meridonal velocity at one time
% PT3d:     potential temperature
% pfull:    temperature coordinate, unit mb
% isenCoord:isentropic coordinate that converts to, increases with index,
% should have regular spacing
% deltaP:   usually 10 mb

%----OUTPUT----
% flux3dIsen: unit kg/(mKs). It is the V(theta) defined in Held and
% Schneider 1999 

G = 9.8; % gravitational constant

minP = 100; %unit mb
maxP = 1050; 
pInterp = minP:deltaP:maxP; % interpolated coordinate for pressure

%flux3dIsen = zeros(size(vcomp3d,1),size(vcomp3d,2),length(isenCoord));
flux2dIsen = zeros(size(vcomp3d,2),length(isenCoord));
deltaPT    = (max(isenCoord)-min(isenCoord))/length(isenCoord);
isenCoordEdge = isenCoord-deltaPT/2;
isenCoordEdge(length(isenCoord)+1) = max(isenCoord)+deltaPT/2;

[x y z]    = ndgrid(1:1:size(vcomp3d,1),1:1:size(vcomp3d,2),pfull); % meshgrid generate X(j,i)
[xi yi zi] = ndgrid(1:1:size(vcomp3d,1),1:1:size(vcomp3d,2),pInterp);

PT3dInterp    = interp3(x,y,z,PT3d,   xi,yi,zi,'cubic');
vcomp3dInterp = interp3(x,y,z,vcomp3d,xi,yi,zi,'cubic');

[~, bin] = histc(PT3dInterp(:),isenCoordEdge);
bin3d    = reshape(bin,size(PT3dInterp));

for i = 1:size(vcomp3dInterp,1)
    for j = 1:size(vcomp3dInterp,2)
        for k = 1:size(vcomp3dInterp,3)
            if bin3d(i,j,k) ~= 0
                flux2dIsen(j,bin3d(i,j,k)) = flux2dIsen(j,bin3d(i,j,k)) ...
                                           + vcomp3dInterp(i,j,k);
            end
        end
    end
end

%{
% using flux1d, very slow; 1000 day takes more than 1 hour
for i = 1:size(vcomp3d,1)
    for j = 1:size(vcomp3d,2)
        vcomp1d = reshape(vcomp3d(i,j,:),size(vcomp3d,3),1);
        PT1d    = reshape(PT3d(i,j,:),   size(PT3d,3),1);
        flux1d  = calIsenFlux1d(vcomp1d,PT1d,pfull,isenCoordEdge,pInterp);
        flux3dIsen(i,j,:) = flux1d(:);
    end
end
%}

%flux3dIsen = flux3dIsen*deltaP*1000/G/deltaPT;
flux2dIsen = flux2dIsen*deltaP*1000/G/deltaPT;