function [F]=ts_diag2(gname,fname,tindex,Smin,Smax,Tmin,Tmax,Rval,...
                      Xmin,Xmax,Ymin,Ymax,color);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2007 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% function [F]=ts_diag2(gname,fname,tindex,Smin,Smax,Tmin,Tmax,Rval,...     %
%                       Xmin,Xmax,Ymin,Ymax,color);                         %
%                                                                           %
% This function plots a TS-diagram from the specified time-record of        %
% temperature and salinity in NetCDF file over the selected region          %
% bounded by coordinates (Xmin,Ymin) to (Xmax,Ymax) for grid depth          %
% greater than Zmin.                                                        %
%                                                                           %
% On Input:                                                                 %
%                                                                           %
%    fname       NetCDF file name (character string).                       %
%    tindex      Time index (integer).                                      %
%    Smin        Minimum salinity to consider.                              %
%    Smax        Maximum salinity to consider.                              %
%    Tmin        Minimum temperature to consider.                           %
%    Tmax        Maximum temperature to consider.                           %
%    Rval        Density contours to draw.                                  %
%    Xmin        Analysis box X-left   coordinate.                          %
%    Xmax        Analysis box X-righy  coordinate.                          %
%    Ymin        Analysis box Y-bottom coordinate.                          %
%    Ymax        Analysis box Y-top    coordinate.                          %
%    Zmin        Minimum local grid bottom to consider (negative, meter).   %                 
%                                                                           %
% On Input:                                                                 %
%                                                                           %
%    F           TS data (structure array):                                 %
%                  F.temp  =>  temperature.                                 %
%                  F.salt  =>  salinity.                                    %
%                  F.den   =>  density anomaly.                             %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Deactivate printing information when reading data from NetCDF file.

global IPRINT

IPRINT=0;

%----------------------------------------------------------------------------
%  Set equation of state expansion coefficients.
%----------------------------------------------------------------------------

A00=+19092.56D0;  A01=+209.8925D0;   A02=-3.041638D0;   A03=-1.852732D-3;
A04=-1.361629D-5; B00=+104.4077D0;   B01=-6.500517D0;   B02=+0.1553190D0;
B03=+2.326469D-4; D00=-5.587545D0;   D01=+0.7390729D0;  D02=-1.909078D-2;
E00=+4.721788D-1; E01=+1.028859D-2;  E02=-2.512549D-4;  E03=-5.939910D-7;
F00=-1.571896D-2; F01=-2.598241D-4;  F02=+7.267926D-6;  G00=+2.042967D-3;
G01=+1.045941D-5; G02=-5.782165D-10; G03=+1.296821D-7;  H00=-2.595994D-7;
H01=-1.248266D-9; H02=-3.508914D-9;  Q00=+999.842594D0; Q01=+6.793952D-2;
Q02=-9.095290D-3; Q03=+1.001685D-4;  Q04=-1.120083D-6;  Q05=+6.536332D-9;
U00=+0.824493D0;  U01=-4.08990D-3;   U02=+7.64380D-5;   U03=-8.24670D-7;
U04=+5.38750D-9;  V00=-5.72466D-3;   V01=+1.02270D-4;   V02=-1.65460D-6;
W00=+4.8314D-4; 

%----------------------------------------------------------------------------
% Compute sigma-t.
%----------------------------------------------------------------------------

Sdelta=(Smax-Smin)/100;
Tdelta=(Tmax-Tmin)/100;

[s,t]=meshgrid(Smin:Sdelta:Smax,Tmin:Tdelta:Tmax);

sqrtS=sqrt(s);

%  Compute density (kg/m3) at standard one atmosphere pressure.

den1 = Q00 + Q01.*t + Q02.*t.^2 + Q03.*t.^3 + Q04.*t.^4 + Q05.*t.^5 + ...
       U00.*s + U01.*s.*t + U02.*s.*t.^2 + U03.*s.*t.^3 + U04.*s.*t.^4 + ...
       V00.*s.*sqrtS + V01.*s.*sqrtS.*t + V02.*s.*sqrtS.*t.^2 + ...
       W00.*s.^2;

den1=den1-1000;

%----------------------------------------------------------------------------
% Read in temperature and salinity.
%----------------------------------------------------------------------------

T=nc_read(fname,'temp',tindex);
S=nc_read(fname,'salt',tindex);

[Lp,Mp,N]=size(T);

%----------------------------------------------------------------------------
% Read in grid coordinates and  Land/Sea mask.
%----------------------------------------------------------------------------

X=nc_read(gname,'lon_rho');
Y=nc_read(gname,'lat_rho');

rmask=nc_read(gname,'mask_rho');

X=repmat(X,[1 1 N]);
Y=repmat(Y,[1 1 N]);
mask=repmat(rmask,[1 1 N]);

%----------------------------------------------------------------------------
% Plot TS-diagram.
%----------------------------------------------------------------------------

ind=find(X < Xmin | X > Xmax | Y < Ymin | Y > Ymax);
T(ind)=NaN;
S(ind)=NaN;

ind=find(mask < 0.5);
T(ind)=NaN;
S(ind)=NaN;

T=reshape(T,1,Lp*Mp*N);
S=reshape(S,1,Lp*Mp*N);

ind=isnan(T);
T(ind)=[];
S(ind)=[];

F.temp=T;
F.salt=S;
F.den=den1;

[c,h]=contour(s,t,den1,Rval); clabel(c,h);
set(h,'edgecolor','k');
set(gcf,'Units','Normalized',...
        'Position',[0.2 0.1 0.6 0.8],...
        'PaperUnits','Normalized',...
        'PaperPosition',[0.2 0.1 0.6 0.8]);
box on;
grid on;
hold on;

plot(S,T,color);
xlabel('Salinity');
ylabel('Temperature');
hold off

return
