clear all;
clf;

% matlab normally only contours when it can interpolate between points
% it will therefore not go outside the areas constrained by obeservation.
% You can force it to close the contours by putting boundaries at the
% edge of the domain.  If 'mod' is specified as on 0 Alexandrium boundaries
% are specidied at the 4 corners of the domain, if it is set to off
% no boundaries are set and the contours are not closed. 

mod = 'on';

% loading in observation data - 93 RMRP, cruise 6a
if strcmp(mod,'on') == 1 
  load ccb_lonlat_mod.mat;
  load ccb_Alexsrf_mod.mat;
else
  load ccb_lonlat.mat;
  load ccb_Alexsrf.mat;
end

% loading coastline data for plotting
load cape_cod_map.dat;
xcoast = cape_cod_map(:,1);
ycoast = cape_cod_map(:,2);

% making grid for interpolation - can change resolution by
% modifying the increment of the x,y vectors.  Right now, it
% is set for a 21x31 grid.

% note that lat/lon are treated as rectangular/cartesian coordinates
% over this domain.  I know this is a bit rough, but figured it would
% be OK for this purpose.
 
X = [-72:0.1:-69];
Y = [41:0.1:43];
[XI, YI] = meshgrid(X,Y);

if strcmp(mod,'on') == 1 
  XX = ccb_lonlat_mod(:,1);
  YY = ccb_lonlat_mod(:,2);
  ZZ = ccb_Alexsrf_mod;
else
  XX = ccb_lonlat(:,1);
  YY = ccb_lonlat(:,2);
  ZZ = ccb_Alexsrf;
end

% interpolate a surface from the data using a cubic interpolation - other
% option include 'linear' and 'nearest' for interpolation method.

ZI = griddata(XX, YY, ZZ, XI, YI, 'cubic');
set(gca,'DataAspectRatio',[1 1 1]);

% plot up the results using contour or pslice

a = contour(XI,YI,ZI,[0:25:250])
hold on
clabel(a,'manual');
plot(xcoast,ycoast,'k-');
cdot(XX,YY,ZZ,jet,15,0,[0 250]);
title('Contoured Cells, 93 RMRP, Cruise 6A')
xlabel('Latitude');
ylabel('Longitude');
pclegend([0 250]', [0.8 0.4 0.06 0.5]);


