Bad contour lines when plotting gridded data

Question

I have heavily edited data from an underwater glider that measures temperature, conductivity, and pressure. From that we get salinity, density and sigma (sigma = density - 1000). The instrument did not sample very well, hence the heavy editing. I need to make "pretty" contour plots for a poster. This is not for a publication so I am not worried about the amount of filtering or smoothing.

I am having trouble with the contour lines of sigma (density), please see below.

The black contour lines should trace the filled contours of sigma but they look very bad. The data was binned by 1 m before any gridding is done. here is the code used to generate the plots.

Here is the code used to generate this image

 % load data
 load Matlab_data
 maxy = 50;

t = 1;
Tstep = t./24./60; % Grid data on t minute time intervals
X1 = time(1)-Tstep:Tstep:time(end)+Tstep; 
Y1 = 0:1:maxy; % Grid data on 1 meter depth intervals

ygrid = [' Depth grid: ' num2str(diff(Y1(1:2)))];
xgrid = [' Time grid: ' num2str(diff(X1(1:2)))];

[X,Y]= meshgrid(X1,Y1);
bad_vals = isnan(sal) .* ~isfinite(sal) .* isnan(press); % don't include NaNs or funky imaginary salinities in the contours
vals = find(bad_vals == 0);

Zd = griddata(time(vals),depth(vals),density(vals),X,Y);
Zt = griddata(time(vals),depth(vals),temp(vals),X,Y);
Zs = griddata(time(vals),depth(vals),sal(vals),X,Y);
Zst = griddata(time(vals),depth(vals),sigmat(vals),X,Y); 

% Interpolate over gaps
vst = interp1gap(sigmat);
vs = interp1gap(sal);
% Grid interpolated salinity and sigma data
Zst_interp = griddata(time(vals),depth(vals),vst(vals),X,Y);
Zs_interp = griddata(time(vals),depth(vals),vs(vals),X,Y);

 %% Contour Plot
  % Set up the figure
  figure3=figure('Position', [2000, 50, 1500, 500]);
  clf
  colormap(jet);

% Temperature
ax1 = subplot(2, 1,1);
  [c,h] = contourf(X,Y,Zst,50,'linestyle','none'); %,[4:.5:9],'linestyle','none');
  cRange = caxis;
   hold on
  [c2,h2] = contour(X,Y,Zst,[24 24.5 25],'color','k'); %,[22:.5:26.5],'linewidth',1.5,'color','k');
  clabel(c2,h2,'fontsize',10,'labelspacing',150);
  set(h2,'linewidth',1)
  hc = colorbar;
  colormap(jet);

  datetick('x','mm/dd','keeplimits','keepticks');

  grid on;
  box on
  pos = get(gca,'position');
  set(gca,'YDir','reverse')%,'position',[pos(1) pos(2) pos(3)-.06 pos(4)]);
  set(gca,'xlim',[time(1)+.5./24 time(end)-.5./24],...
      'ylim',[0 maxy],'fontsize',8,'xminortick','on','yminortick','on');
  set(get(hc,'ylabel'),'string','Sigma-Theta (kg m^-^3)','rotation',270,'verticalalignment','bottom');
  ylabel('Ocean Depth (m)');
  xlabel('Date');
  %title(['Sigma Theta (kg m^-^3) -' strrep(tag,'_','\_')], 'fontweight', 'bold','FontSize',12)
  title(['Sigma Theta (kg m^-^3) :' ygrid xgrid ], 'fontweight', 'bold','FontSize',12)    

ax2 = subplot(2, 1,2);
  %h=pcolor(X,Y,Zst);
  [c,h] = contourf(X,Y,Zst_interp,50,'linestyle','none'); %,[4:.5:9],'linestyle','none');
  shading interp
  hc = colorbar;
  cRange = caxis;
  hold on
  [c2,h2] = contour(X,Y,Zst_interp,[24 24.5 25],'color','k'); %,[22:.5:26.5],'linewidth',1.5,'color','k');
  clabel(c2,h2,'fontsize',10,'labelspacing',150);
  set(h2,'linewidth',1)
  hc = colorbar;
  colormap(jet);
  caxis(cRange);

  datetick('x','mm/dd','keeplimits','keepticks');

  grid on;
  box on
  pos = get(gca,'position');
  set(gca,'YDir','reverse')%,'position',[pos(1) pos(2) pos(3)-.06 pos(4)]);
  set(gca,'xlim',[time(1)+.5./24 time(end)-.5./24],...
      'ylim',[0 maxy],'fontsize',8,'xminortick','on','yminortick','on');
  set(get(hc,'ylabel'),'string','Sigma-Theta (kg m^-^3)','rotation',270,'verticalalignment','bottom');
  ylabel('Ocean Depth (m)');
  xlabel('Date');
 % title(['Sigma Theta (kg m^-^3) -' strrep(tag,'_','\_')], 'fontweight', 'bold','FontSize',12)

  title(['Sigma Theta interp1gap (kg m^-^3) :' ygrid xgrid ], 'fontweight', 'bold','FontSize',12)

Please help! I have TWO issues... the first and most obvious is the ugly black contour lines, these should "flow" around the filled contours nicely. If anyone has experience or suggestions of how to smooth these please pass it along, with code if possible.

The second issue lies in the bottom plot, I need to fill the gaps in data (due to editing original bad data). I used the function interp1gap, available on the File exchange but it also interpolates the data to deeper depths, which I do not want. I just want the gaps to be filled in, such as by choosing their horizontal neighbors and filling in.

Please let me know if you have any suggestions for fixing this! I have attached the data (Matlab_data.mat) and it includes time, depth, sal (salinity), sigmat, press, temp, and density.

Is this an issue of gridding? Please be as specific as possible and any code and figures would be greatly appreciated if you have time to look into this.

The data is available drop box https://www.dropbox.com/s/mjd8f9bzdvwddk5/Matlab_data.mat?dl=0

Thank you very much in advance!