%  MakeGrid.m
%     make a 3d grid of values from a cruise
%     1000 m is depth interval 19
%
CruiseName='NBP0103';
%CruiseName='NBP0104';

%      define these once at beginning
%makemask=1;
%loaddep=1;

ptcl=[-2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.45 -0.2 0 .2 .4 .6 .8 1 1.2 1.4 1.6 1.8 2];
scl=[32.2 32.4 32.6 32.8 33 33.2 33.4 33.6 33.8 34 34.2 34.4 34.6 34.7 34.8];
makeplot=1;
skipplot=0;

load([CruiseName '.gridstns.mat']);

load([CruiseName '.depgrid.mat'],'depgrid');
makemask=0;

i=find(xgrid < 0);
pt(:,i) = NaN;
t(:,i)  = NaN;
s(:,i)  = NaN;
o2(:,i) = NaN;

ndep=19;
xi=20:20:500;
yi=-20:20:300;
[n,nx]=size(xi);
[n,ny]=size(yi);
[xxi,yyi]=meshgrid(xi,yi);
stddepgrid=stddep(1:ndep);

nzgrid=ones(ny,nx);
ptgrid=ones(ny,nx,ndep);
tgrid=ones(ny,nx,ndep);
sgrid=ones(ny,nx,ndep);
o2grid=ones(ny,nx,ndep);

if loaddep
%                    load etopo2 for west Ant Pen region
%                    loads files londep, latdep, zdep, xdep, ydep
  load wap.etopo2.mat;
  loaddep=0;
end
if makemask
  disp(' Making land mask. ')
  depgrid=ones(ny,nx);
  depgrid=griddata(xdep,ydep,zdep,xxi,yyi);
  makemask=0;
  im=find( depgrid >  0);
  if skipplot
%  if makeplot
    figure
    [Con,Han]=contour(xxi,yyi,depgrid);
    clabel(Con,Han);
    title('depth on 20 km grid ');
    axis([0 500 -50 300]);
    hold on;
    plot(xgrid,ygrid,'r*');
    hold off
  end
else
  im=find( depgrid >  0);
end
%           extract a depth interval
for idep=1:ndep
%for idep=1:2
   disp(['Extracting depth level: ' num2str(idep) ' Depth: ' num2str(stddep(idep))])
%
   xpt=pt(idep,:);
   xt = t(idep,:);
   xs = s(idep,:);
   xo2=o2(idep,:);
%
   i=find(~isnan(xt));
   ti=griddata(xgrid(i),ygrid(i),xt(i),xxi,yyi);
   ti(im)=NaN;
   tgrid(:,:,idep)=ti;
%
   i=find(~isnan(xpt));
   pti=griddata(xgrid(i),ygrid(i),xpt(i),xxi,yyi);
   pti(im)=NaN;
   ptgrid(:,:,idep)=pti;
%
   i=find(~isnan(xs));
   si=griddata(xgrid(i),ygrid(i),xs(i),xxi,yyi);
   si(im)=NaN;
   sgrid(:,:,idep)=si;
%
   i=find(~isnan(xo2));
   o2i=griddata(xgrid(i),ygrid(i),xo2(i),xxi,yyi);
   o2i(im)=NaN;
   o2grid(:,:,idep)=o2i;
%
end
for i=1:nx
  for j=1:ny
    for idep=1:ndep
       if stddep(idep) < -depgrid(j,i) 
         nzgrid(j,i)=idep;
       end
    end
  end
end
%         set land to zero
nzgrid(im)=0;
if makeplot
  figure
%  cl=[1 5 10 11 12 13 14 15 16 17 18 19];
%  [Con,Han]=contour(xxi,yyi,nzgrid,cl);
  [Con,Han]=contour(xxi,yyi,nzgrid);
  clabel(Con,Han);
  title('deepest index on 20 km grid ');
  axis([0 500 -50 300]);
  hold on;
  plot(xgrid,ygrid,'r*');
     text(50,0,CruiseName)
  hold off
end

save([CruiseName '.depgrid.mat'],'stddepgrid','xi','yi', ...
     'ptgrid','tgrid','sgrid','o2grid', ...
     'depgrid','nzgrid');

if skipplot
%if makeplot
   for idep=1:ndep
     figure
     [Con,Han]=contour(xxi,yyi,ptgrid(:,:,idep));
     clabel(Con,Han);
     title(['Pot Temp on 20 km grid at depth ' num2str(stddep(idep))]);
     axis([0 500 -50 300]);
     hold on;
     i=find(~isnan(pt(idep,:)));
     plot(xgrid(i),ygrid(i),'r*');
     text(50,0,CruiseName)
     hold off
  end
end
if skipplot
%if makeplot
   for idep=1:ndep
     figure
     [Con,Han]=contour(xxi,yyi,sgrid(:,:,idep));
     clabel(Con,Han);
     title(['Salinity on 20 km grid at depth ' num2str(stddep(idep))]);
     axis([0 500 -50 300]);
     hold on;
     i=find(~isnan(s(idep,:)));
     plot(xgrid(i),ygrid(i),'r*');
     text(50,0,CruiseName)
     hold off
  end
end

