%  printO2onMap.m  - print oxygen on a given depth or density
%                    level at horizontal locations on a regional map.
%
%   need variables set: ShadeCoast PlotCoast PlotBathy UseDepth Depth2Plot
%                       Dens2plot RefPress xrange yrange o2units
%                       DataDir CruiseName DepthResolution
%   sets variables: none
%   calls: readctd swstate

global GEOMINFO CRUISEINFO DEFAULTS

%    title for figure
if(UseDepth == 1)
   zi=[-Depth2Plot -Depth2Plot];
   figtitle=['Oxygen at ' num2str(Depth2Plot) ' m'];
else
   RefPress=0;
   zi=[Dens2Plot Dens2Plot];
   figtitle=['Oxygen on ' num2str(Dens2Plot) ' sigma'];
end

%        open an new graphics window.
figure;
if(ShadeCoast)
   [Cc,hc]=contourf(xdep,ydep,zdep,[0 0]);
end
if(PlotCoast) 
   plot(xcoast,ycoast);
end
hold on;
if(PlotBathy) 
%      contour intervals for bottom topography
   [C,hc]=contour(xdep,ydep,zdep,botcl);
   %clabel(C,hc);
end
axis([xrange yrange]);
title(figtitle);
xlabel('distance (km)');
ylabel('distance (km)');

val=NaN*ones(size(ctdxgrid));

Casts=0;
for cast=ctdcastlist
%           count number of casts
   Casts=Casts+1;
  file= [DataDir '/' CruiseName '/Ctd/' DepthResolution '/ctd' castnum(cast) '.dat'];

  [nlines,ndepths,desc,p,t,s,pt,r0,r1,bv,t0,t1,s0,s1,o2,par,tr,fl] = readctd(file);

%          convert pressure to depth
   Z=-p;
   CastDepth(cast)=max(p);
   X(cast)=ctdxgrid(cast);
   Y(cast)=ctdygrid(cast);
%                  reduce to 10 m resolution
   [m,n]=size(p);
   if(UseDepth == 1)
      if(isnan(t))
         val(cast)=NaN;
      else
         ti=interp1(Z(1:10:m),o2(1:10:m),zi,'linear');
         val(cast)=ti(1);
      end
   else
      [SVAN,Pden]=swstate(s(1:10:m),pt(1:10:m),RefPress);
%        remove unstable densities by pushing values
%        deeper if next deeper value is less dense
%        add a small amount to have a vertical difference
      i=find(diff(Pden) <= 0);
      while  (~isempty(i))
         Pden(i+1)=Pden(i)+1.e-5;
         i=find(diff(Pden) <= 0);
      end
      if(isnan(t))
         val(cast)=NaN;
      else
         ti=interp1(Pden,o2(1:10:m),zi,'linear');
         val(cast)=ti(1);
      end
   end
end

i=find(CastDepth < Depth2Plot);
val(i)=NaN;

i=find(~isnan(val));
%    convert O2 from micro mol per kg to ml per l
val(i) = val(i)*o2units;

text(X(i),Y(i),num2str(val(i),3),'Color',[1 0 0], ...
     'FontSize',8,'HorizontalAlignment','center');

hold off;

%    end of printO2onMap.m

