flux2dIsen.m

function flux2dIsen = calIsenFlux(vcomp3d,PT3d,pfull,isenCoord,deltaP,minP,maxP,ps2d)
% calculate mass flux in isentropic coordinate
%----INPUT----
% vcomp3d:  meridonal velocity at one time,   3d matrix (lon,lat,press)
% PT3d:     potential temperature,            3d matrix
% pfull:    temperature coordinate, unit mb,  1d array, for PT3d and
%   vcomp3d
% isenCoord:isentropic coordinate that converts to, increases with index, 1d array
%  usually should have regular spacing, or unexpect problem can happen
% deltaP:   usually 10 mb
% minP, maxP and deltaP describe the pressure coordinate to interpolate to
% ps2d:     pressure at surface, unit mb,     2d matrix (lon,lat)
%  if ps2d is not present, does not place a bound on surface pressure

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

vcomp3d = double(vcomp3d);
PT3d    = double(PT3d);
pfull   = double(pfull);

G = 9.8; % gravitational constant

if ~ exist('minP','var')
minP = 100; %unit mb
end
if ~ exist('maxP','var')
maxP = 1050;
end

pInterp = minP:deltaP:maxP; % interpolated coordinate for pressure
numLon = size(vcomp3d,1);
numLat = size(vcomp3d,2);

flux2dIsen = zeros(numLat,length(isenCoord));
deltaPT    = (max(isenCoord)-min(isenCoord))/length(isenCoord);
isenCoordEdge = isenCoord-deltaPT/2; % edges to collect vcomp into isentropic coord
isenCoordEdge(length(isenCoord)+1) = max(isenCoord)+deltaPT/2;

% -------- interpolate PT and vcomp to higher resolution pressure grid -----------------
PT3dInterp    = nakeinterp3(numLon,numLat,pfull(:),PT3d,   pInterp(:));
vcomp3dInterp = nakeinterp3(numLon,numLat,pfull(:),vcomp3d,pInterp(:));
% -------- end interpolation -----------------------------------------------------------

% -------- eliminate vcomp outside surface pressure range ------------------------------
if exist('ps2d','var')
     ps3d      = repmat(ps2d,[1 1 length(pInterp)]);
     pInterp3d        = zeros(size(ps3d));
     pInterp3d(1,1,:) = pInterp(:);
     pInterp3d = repmat(pInterp3d(1,1,:),[numLon numLat 1]);
     outsideIndex = pInterp3d(:) > ps3d(:);
     vcomp3dInterp(outsideIndex) = 0;
end
% -------- end eliminating ------------------------------------------------------------


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

flux2dIsen = p2isen(vcomp3dInterp,bin,...
            [numLon numLat size(vcomp3dInterp,3) length(isenCoord)]);

flux2dIsen = flux2dIsen*deltaP*100/G/deltaPT/numLon;