%  MakeSigGrid.m
%     make a 3d grid of values from a cruise
%     interpolated onto density surfaces
%
%    define to load depth
%  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];
siglevels=[26.8 26.9 27 27.1 27.2 27.3 27.4 27.5 27.6 27.65 ...
27.7 27.725 27.75 27.775 27.8 27.825 ];

makeplot=1;
skipplot=0;
RefPress=0;
PlotRange=[-50 550 -50 350];
CruiseName='NBP0103';
%CruiseName='NBP0104';

load([CruiseName '.depgrid.mat'],'depgrid');
im=find( depgrid >  0);
load([CruiseName '.gridstns.mat']);

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

if loaddep
%                    load etopo2 for west Ant Pen region
%                    loads files londep, latdep, zdep, xdep, ydep
  load wap.etopo2.mat;
  loaddep=0;
end

[m nsig]=size(siglevels);
xi=20:20:500;
yi=-20:20:300;
[n,nx]=size(xi);
[n,ny]=size(yi);
[xxi,yyi]=meshgrid(xi,yi);

PTonSigGrid=ones(ny,nx,nsig);
SonSigGrid=ones(ny,nx,nsig);
DepSigGrid=ones(ny,nx,nsig);
%        station data on sigma surfaces
Depsig=zeros(size(pt));
PTsig=zeros(size(pt));
Tsig=zeros(size(pt));
Ssig=zeros(size(pt));
O2sig=zeros(size(pt));
[m,nstn]=size(pt);

p=stddep*ones(1,nstn);

%           make potential density
PT=potentialtemp(s,t,p,RefPress);
[SVAN,Pden]=swstate(s,PT,RefPress);
for i=1:nstn
   Pdenmin=min(Pden(:,i));
   Pdenmax=max(Pden(:,i));
   j=find( siglevels >= Pdenmin & siglevels <=Pdenmax);
   zi=NaN*ones(33,1);
   val=Pden(:,i);
   k=find(~isnan(val));
   zi(j)=interp1(val(k),stddep(k),siglevels(j)');
   Depsig(:,i)=zi;
%         interpolate pt and s on sigma depths
   val=PT(:,i);
   k=find(~isnan(val));
   vali=NaN*ones(33,1);
   vali(j)=interp1(stddep(k),val(k),zi(j));
   PTsig(:,i)=vali;
%
   val=s(:,i);
   k=find(~isnan(val));
   vali=NaN*ones(33,1);
   vali(j)=interp1(stddep(k),val(k),zi(j));
   Ssig(:,i)=vali;
%
   val=t(:,i);
   k=find(~isnan(val));
   vali=NaN*ones(33,1);
   vali(j)=interp1(stddep(k),val(k),zi(j));
   Tsig(:,i)=vali;
%
   val=o2(:,i);
   k=find(~isnan(val));
   vali=NaN*ones(33,1);
   vali(j)=interp1(stddep(k),val(k),zi(j));
   O2sig(:,i)=vali;
end

%      grid a sigma surface
for isig=1:nsig
%for isig=1:2
   disp(['Gridding sigma level: ' num2str(isig) ' Sigma: ' num2str(siglevels(isig))])
   xv=Depsig(isig,:);
   i=find(~isnan(xv));
   if(isempty(i))
     DepSigGrid(:,:,isig)=NaN;
   else
     vi=griddata(xgrid(i),ygrid(i),xv(i),xxi,yyi);
     vi(im)=NaN;
     DepSigGrid(:,:,isig)=vi;
     i=find(DepSigGrid(:,:,isig) > -depgrid);
     DepSigGrid(i)=NaN;
   end
%
   xv=PTsig(isig,:);
   i=find(~isnan(xv));
   if(isempty(i))
     PTonSigGrid(:,:,isig)=NaN;
   else
     vi=griddata(xgrid(i),ygrid(i),xv(i),xxi,yyi);
     vi(im)=NaN;
     PTonSigGrid(:,:,isig)=vi;
   end
%
   xv=Tsig(isig,:);
   i=find(~isnan(xv));
   if(isempty(i))
     TonSigGrid(:,:,isig)=NaN;
   else
     vi=griddata(xgrid(i),ygrid(i),xv(i),xxi,yyi);
     vi(im)=NaN;
     TonSigGrid(:,:,isig)=vi;
   end
%
   xv=Ssig(isig,:);
   i=find(~isnan(xv));
   if(isempty(i))
     SonSigGrid(:,:,isig)=NaN;
   else
     vi=griddata(xgrid(i),ygrid(i),xv(i),xxi,yyi);
     vi(im)=NaN;
     SonSigGrid(:,:,isig)=vi;
   end
%
   xv=O2sig(isig,:);
   i=find(~isnan(xv));
   if(isempty(i))
     O2onSigGrid(:,:,isig)=NaN;
   else
     vi=griddata(xgrid(i),ygrid(i),xv(i),xxi,yyi);
     vi(im)=NaN;
     O2onSigGrid(:,:,isig)=vi;
   end
end

save([CruiseName '.siggrid.mat'],'DepSigGrid','siglevels','xi','yi', ...
   'PTsig','Tsig','Ssig','O2sig','RefPress', ...
   'PTonSigGrid','TonSigGrid','SonSigGrid','O2onSigGrid', ...
   'depgrid');

%if skipplot
if makeplot
   for isig=1:nsig
     figure
     i=find(~isnan(DepSigGrid(:,:,isig)));
     if(~isempty(i))
       [Con,Han]=contour(xxi,yyi,DepSigGrid(:,:,isig));
       clabel(Con,Han);
     else
       text(250,200,'Density surface does not occur')
     end
     title(['Depth of sigma ' num2str(siglevels(isig))]);
     axis(PlotRange);
     hold on;
     i=find(~isnan(PTsig(isig,:)));
     plot(xgrid(i),ygrid(i),'r*');
     i=find(isnan(DepSigGrid(:,:,isig)));
     plot(xxi(i),yyi(i),'go');
     text(50,0,CruiseName)
     hold off
  end
end
if skipplot
%if makeplot
   for isig=1:nsig
     figure
     [Con,Han]=contour(xxi,yyi,TonSigGrid(:,:,isig));
     clabel(Con,Han);
     title(['Gridded Temp on sigma ' num2str(siglevels(isig))]);
     axis(PlotRange);
     hold on;
     i=find(~isnan(PTsig(isig,:)));
     plot(xgrid(i),ygrid(i),'r*');
     text(50,0,CruiseName)
     hold off
  end
end
if skipplot
%if makeplot
   for isig=1:nsig
     figure
     [Con,Han]=contour(xxi,yyi,PTonSigGrid(:,:,isig));
     clabel(Con,Han);
     title(['Gridded Pot Temp on sigma ' num2str(siglevels(isig))]);
     axis(PlotRange);
     hold on;
     i=find(~isnan(PTsig(isig,:)));
     plot(xgrid(i),ygrid(i),'r*');
     text(50,0,CruiseName)
     hold off
  end
end
if skipplot
%if makeplot
   for isig=1:nsig
     figure
     [Con,Han]=contour(xxi,yyi,SonSigGrid(:,:,isig));
     clabel(Con,Han);
     title(['Gridded Salinity on sigma ' num2str(siglevels(isig))]);
     axis(PlotRange);
     hold on;
     i=find(~isnan(Ssig(isig,:)));
     plot(xgrid(i),ygrid(i),'r*');
     text(50,0,CruiseName)
     hold off
  end
end
if skipplot
%if makeplot
   for isig=1:nsig
     figure
     [Con,Han]=contour(xxi,yyi,O2onSigGrid(:,:,isig));
     clabel(Con,Han);
     title(['Gridded Oxygen on sigma ' num2str(siglevels(isig))]);
     axis(PlotRange);
     hold on;
     i=find(~isnan(Ssig(isig,:)));
     plot(xgrid(i),ygrid(i),'r*');
     text(50,0,CruiseName)
     hold off
  end
end

