%  Tmaxmap.m   - contour maximum potential temperature below a given depth
%                at horizontal locations on a regional map.
%
%   need variables set: ShadeCoast PlotCoast PlotBathy xrange yrange MinDepth
%                       UsePotTemp RefPress DataDir CruiseName
%                       DepthResolution ctdcastlist xctdcastlist
%   sets variables: MaskDone
%   calls: readctd readxbt sw_ptmp

global GEOMINFO CRUISEINFO DEFAULTS MASKDONE

figtitle=['Max Potential Temperature below ' num2str(MinDepth) ' 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) 
   [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);
   Lat(cast)=ctdlat(cast);
   Lon(cast)=ctdlon(cast);
   if UsePotTemp 
%   calculate the potential temperature
      pt=ones(size(t))*NaN;
      pt = sw_ptmp(s,t,p,RefPress);
      t=pt;
   end
   if CastDepth(cast) > MinDepth
     i=find(p > MinDepth);
     val(cast)=max(t(i));   
   else
     val(cast)=NaN;
   end
end

if xctdcastlist(1) == 0
   disp('No Xctd data')
else
%                 need to fix this part
   for cast=xctdcastlist
%          count number of casts
     Casts=Casts+1;
     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(Casts)=max(P);
      X(Casts)=xctdxgrid(xccastlist(cast));
      Y(Casts)=xctdygrid(xccastlist(cast));
      Lat(Casts)=xctdlat(xccastlist(cast));
      Lon(Casts)=xctdlon(xccastlist(cast));
      if UsePotTemp 
%   calculate the potential temperature
         PT=ones(size(T))*NaN;
         PT = sw_ptmp(S,T,P,RefPress);
         T=PT;
      end
   end
end

%           remove casts that are shallower than the plot depth
i=find(CastDepth < MinDepth);
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);

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

if MaskLand == 1
%            mask grid values over land
  if MaskDone == 1
    [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;

%               print values to a file
fid=fopen('MaxTemp.asc','w')
for cast=1:Casts
   fprintf(fid,'%9.4f  %9.4f  %9.6f \n',Lat(cast),Lon(cast),val(cast));
end 
fclose(fid)

drawnow
%    end of Tmaxmap.m

