function [S,T,TSvol]=TSvolumetric(theta,salt,gamma,lon,lat,depth,gammarange,thetarange,saltrange);

% Use [S,T,TSvol]=TSvolumetric(theta,salt,gamma,lon,lat,depth,gammarange,thetarange,saltrange)
%
% This routine calculates the volumetric theta-S diagram for waters with neutral densities between gammarange(1) and gammarange(2). The routine will calculate the volume span by theta-S values that fall between the elements in the thetarange and saltrange values
%
% theta is the 3D potential temperature field 
% salt is the 3D salinity field 
% gamma is the 3D neutral density field
% lon is the 1D longitude field
% lat is the 1D latitude field
% depth is the 1D pressure field
% gammarange is a two element vector that gives the maximum and minimum neutral density to be considered 
% thetarange is a 1D vector with the potential temperature bins
% saltrange is a 1D vector with the temperature bins
%
% Example
% [S,T,TSvol]=TSvolumetric(t,s,g,lon,lat,press,[28.1,28.5],[-1:.1:4],[34.4:.01:35]);

disp('Setting up parameters');
% Compute the lateral area span by each point in the data set
[X,Y]=meshgrid(lon,lat);
R=sw_dist([-90 90]/pi,[0 0],'km')*1e3;
ll=size(X);
dlat=mean(diff(lat(:)))*pi/180;
dlon=mean(diff(lon(:)))*pi/180;
xx=X*pi/180;yy=Y*pi/180;
al(2:ll(1),:)=2*pi*(sin((yy(1:ll(1)-1,:)+yy(2:ll(1),:))/2)-sin(lat(1)*pi/180));
al(1,:)=2*pi*(sin((yy(1,:)-dlat+yy(1,:))/2)-sin(lat(1)*pi/180));
al(ll(1)+1,:)=2*pi*(sin((yy(end,:)+dlat+yy(end,:))/2)-sin(lat(1)*pi/180));
al=diff(al)/ll(2)*R^2;

% Compute the height span by  each point in the data set
press=[0;diff(depth(:))];

% Compute the volume span by  each point in the data set
for i=1:length(press)     
     volume(:,:,i)=al*press(i);
end

% Compute the range of neutral densities
g1=gammarange(1);
g2=gammarange(2);

kk=find((gamma>g1) & (gamma<g2));
ss=salt(kk);
tt=theta(kk);
vol=volume(kk);

disp('Calculating T-S volumetric diagram');
x=ss(:);
y=tt(:);
xbins=saltrange(:);
ybins=thetarange(:);
nxbins=length(saltrange);
nybins=length(thetarange);
TSvol=zeros(nxbins,nybins);
ii=~isnan(x) & ~isnan(y);
x=x(ii);
y=y(ii);
temp=[xbins(1);(xbins(1:end-1)+xbins(2:end))/2;xbins(end)];
intx=diff(temp);
intx(1)=2*intx(1);intx(end)=2*intx(end);
posx=[xbins(1);(temp(2:end-2)+temp(3:end-1))/2;xbins(end)];
temp=[ybins(1);(ybins(1:end-1)+ybins(2:end))/2;ybins(end)];
inty=diff(temp);
inty(1)=2*inty(1);inty(end)=2*inty(end);
posy=[ybins(1);(temp(2:end-2)+temp(3:end-1))/2;ybins(end)];
norm=intx*inty';

num=length(x);
for n=1:num
   [dum,i]=min(abs(x(n)-saltrange));
   [dum,j]=min(abs(y(n)-thetarange));
   TSvol(i,j)=TSvol(i,j)+vol(n);
end
TSvol(1,:)=zeros(size(TSvol(1,:)));
TSvol(:,end)=zeros(size(TSvol(:,end)));
S=posx;
T=posy;