%  DynTopomap.m - contour dynamic topography on a regional map.
%                 plots surface dynamic topography relative to RefPress.
%
%   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;

xi=linspace(min(xrange),max(xrange),MapNx);
yi=linspace(min(yrange),max(yrange),MapNy);
[Xi,Yi]=meshgrid(xi,yi);
i=find(~isnan(val));
Vali=griddata(X(i),Y(i),val(i),Xi,Yi);

%    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);


MaxVal=ceil(10*max(val(i)));
MinVal=floor(10*min(val(i)));
DiffVal=MaxVal-MinVal;
if DiffVal < 5 
  cl=MinVal:.5:MaxVal;
elseif DiffVal < 15 
  cl=MinVal:MaxVal;
elseif DiffVal < 30
  cl=MinVal:2:MaxVal;
else
  cl=MinVal:3:MaxVal;
end
cl=.1*cl;

plot(X(i),Y(i),'r.')

if MaskLand == 1
%            mask grid values over land
  if MaskDone == 1
%            check existing mask for size
    [m n]= size(MapMask);
    if n ~= MapNx | m ~= MapNy
%         map not on same grid as MapNx, MapNy
%         need to recalculate the map
       MaskDone=0;
    end
  end
%      
  if MaskDone == 0
     disp(' Making land mask. ');
     MapMask=griddata(xdep,ydep,zdep,Xi,Yi);
     MaskDone = 1;
  end
  i=find( MapMask >  0);
  Vali(i)=NaN;
end

[C,hc]=contour(Xi,Yi,Vali,cl);
clabel(C,hc);

hold off;
drawnow;

%    end of DynTopomap.m


