%  printDynTopomap.m - prints surface dynamic topography relative to
%                      RefPress on a regional map.
%
%   need variables set: ShadeCoast PlotCoast PlotBathy MapNx MapNy
%                       RefPress xrange yrange ctdcastlist xctdcastlist
%                       DataDir CruiseName DepthResolution
%   sets variables: none
%   calls: readctd readxctd sw_svan

global GEOMINFO CRUISEINFO DEFAULTS MASKDONE

%    title for figure
figtitle=['Dynamic Topography relative to ' num2str(RefPress) ' m  for ' CruiseName];
%        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,h]=contour(xdep,ydep,zdep,botcl);
   %clabel(C,h);
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);
%                         calculate specific volume anomaly
   svan=sw_svan(s,t,p);
   i=find(p > RefPress);
%                remove contribution below reference pressure
   if(~isempty(i))
      svan(i)=0;
   end
   [m,n]=size(p);
   avgsvan=.5*(svan(1:m-1)+svan(2:m));%
%              convert pressure from dbar to pascal
   dp=1.e4*(p(2:m)-p(1:m-1));
%       integrate with composite trapezoidal rule
   IntSvan=avgsvan' * dp;

   if(isnan(t))
      val(cast)=NaN;
   else
      val(cast)=IntSvan;
   end
end

if xctdcastlist(1) == 0
   disp('No Xctd data')
else
  for cast=xctdcastlist
    file= ['XCTD' num2str(xcFileno(xccastlist(cast))) '.dat'];
%               read a profile of T, S
     [P,T,S]=readxctd(xcdatadir,file);
%          convert pressure to depth
     Z=-P;
     CastDepth(cast+casts)=max(P);
     X(cast+casts)=xcXloc(xccastlist(cast));
     Y(cast+casts)=xcYloc(xccastlist(cast));
%                         calculate specific volume anomaly
     svan=sw_svan(S,T,P);
     i=find(P > RefPress);
%                remove contribution below reference pressure
     if(~isempty(i))
        svan(i)=0;
     end
     [m,n]=size(P);
     avgsvan=.5*(svan(1:m-1)+svan(2:m));%
%              convert pressure from dbar to pascal
     dp=1.e4*(P(2:m)-P(1:m-1));
%       integrate with composite trapezoidal rule
     IntSvan=avgsvan' * dp;
     if(isnan(T))
        val(cast+casts)=NaN;
     else
        val(cast+casts)=IntSvan;
     end
  end
end

%           remove casts that are shallower than the reference depth
i=find(CastDepth < RefPress);
val(i)=NaN;

%shortval=round(val*100);
%i=find(~isnan(shortval));
%text(X(i)',Y(i)',int2str(shortval(i)));
i=find(~isnan(val));
text(X(i),Y(i),num2str(val(i),3),'Color',[1 0 0], ...
     'FontSize',8,'HorizontalAlignment','center');

%    print numbers
%fid=fopen('dyntopo.asc','w');
%for cast=1:size(i')
%   fprintf(fid,'%9.4f  %9.4f  %9.6f \n',Lat(cast),Lon(cast),val(cast));
%   fprintf(fid,'%6.0f  %6.0f  %9.6f \n',X(cast),Y(cast),val(cast));
%end 
%fclose(fid);


hold off;
drawnow;

%    end of printDynTopomap.m

