function [xx,yy,z3,lon,lat,var] = dospgrid ( variable, ...
                                             gridfile, ...
                                             cyear, ...
					     crnum, ...
					     legnum, ...
					     depth, ...
					     vmin, ...
					     vmax, ...
					     nrange);
% DOSPGRID  plots color shaded horizontal sections in the 
%          Western Gulf of Maine from a .mat data file for a cruise:
%
%    XXgomYz.mat (XX=yr=cyear, Y=cruise number=crnum, z=leg (a or b)=legnum)
%
%................................................................................
%  INPUT VARIABLE DESCRIPTIONS:
% 
%   varstr= variable to be plotted
%		(e.g., 'salt' for salinity, 'temp' for temperature, 'alex' for
%                      alexandrium tameranse, 'flor' for flourometry data,
%   gridfile:  Previously existing grid.
%		
%   cyear = cruise year (93 or 94)
%   crnum = cruise number (2-9 for 94 or 2-6 for 93)
%   legnum = 1 or 2 for leg number
%   depth    = depth of horizontal section to be displayed
%   vmin = minimum value for plot, (e.g. 28 for salinity)
%   vmax = maximum value for plot, (e.g. 32 for salinity)
%   nrange = range of data influence in grid cells.  Grid
%          cells farther than nrange cells away from data
%          are masked.  Recommended value=5-10.
%
%  dospgrid plots a color shaded horizontal section 
%  of a specified 
%  variable at a specified water depth
%
%  Rich Signell and Derek Fong  April 1994
%  
%
% modified 11/4/94 DAF
% modified 12/15/95 JGE
% 

% load the grid and associated parameters, ie, inner r1, r2, theta1, theta2,
% etc.
eval ( sprintf ('load %s', gridfile) );

%set limits on color axis
cax=[vmin vmax];


% Construct the filename string for the proper matlab file.
% Then extract the specified variable at the specified depth.
filename = sprintf ( '%igom%i%c', cyear, crnum, legnum+96 );
[var, lon, lat] = readmat ( filename, variable, depth );
if ( length(var) == 0 )
  fprintf(1, 'dospgrid:  could not read %s\n', filename );
  return;
end


% Convert data into polar coordinates. 
centered_lon = lon - c_lon;
centered_lat = lat - c_lat;
data_r = sqrt ( centered_lon.^2 + centered_lat.^2 );
data_theta = atan2( centered_lat, centered_lon );

% atan2 returns values between -pi and pi.  We need the values between 0 and
% 2pi.  Otherwise we get confused when trying to put the data into the grid. 
ifind = find(data_theta<0);
data_theta(ifind) = data_theta(ifind) + 2*pi;


% Put the data into the grid.
r_ind = 1 + round((data_r - r1)/dr);
theta_ind = 1 + round((data_theta - theta1)/dtheta);
n = length(data_r);
for i = 1:n
   polar_grid(theta_ind(i),r_ind(i)) = var(i);
end



% parameters for gridding
T=0.9;
tol=0.1;
omega=1.4;
alpha=3.0;
nrange=10;
% new, supposedly working algorithm
% Oh, it works all right.
z2=ngrid(polar_grid,T,alpha,tol,omega,nrange);  


% put up a little plot showin' what we done
fig = figure;
surf ( x_grid, y_grid, z2 );
set(gca, 'view', [0 90] );
hold on
plot3 ( lon, lat, 1e6*ones(size(lon)), '+w' );
set ( gca, 'aspectratio', [cos(43*pi/180) NaN] );
shading interp





return;











% Derek's stuff, I'm not going to mess around with it.
%....................................................................
% Plotting stuff

z3=flipud(rot90(z2));
ind=find(finite(z3));

% mark of -inf values with nans so that some versions of matlab
% won't goof and plot these points as a rust red.
%
ind2=find(z3==-inf);
z3(ind2)=nan*ones(1,length(ind2));

%
% peg data out of range by max and min colors
%
z3(ind)=max(z3(ind),cax(1));
z3(ind)=min(z3(ind),cax(2));

%figure(1);
colormap(jet);

%.................................................................
% need to decide if we want a shaded contour map or a line contour
% plot...conflag=1 (shaded map), conflag=0 (line contour)
%
conflag=1;
if conflag
%   clegend4(cax');
%   hmesh=pcolor2(xx(1,:),yy(:,1),z3);caxis(cax);
   zsize=size(get(hmesh,'zdata'));
   set(hmesh,'zdata',ones(zsize)*(-3));
else
   cs=contour(xx(1,:),yy(:,1),z3);
   vc=26:.5:33;
   clabel(cs,'manual');
end

set(gca,'aspectratio',[NaN, cos(42*pi/180)]);shading('flat');
%..................................................................
% now plot topography based on gomcoast.mat data set
%
hold on;
doland;
hold off;


% plots a dot for each CTD station
line(lon,lat,'linestyle','.','color','black',...
  'erasemode','none');
set(gcf,'inverthardcopy','off')

% makes a color map such that blue is high values and red are low values
% (except of alexandrium tamarense or temperature)
if (varstr~='alex' & varstr~='temp')
  cc=colormap;
  colormap(flipud(cc));
end
if conflag
    hh=gca;
    diff=max(cax)-min(cax);
    mid=floor(diff/4);
    inc=floor(diff/4)/10;
%    hb=cbar(hh,'v',inc,mid,'fontsize',9);
    axes(hh);
    set(gca,'fontsize',9);
end





