function contour_temp(theSite)
%
%   contour_temp(theSite)
%
%   Function to create contour plots of temperature for theSite.

%Soupy Alexander, 4/2004

%This Matlab m-file was used to create portions of U.S. Geological Survey 
%Data Series 85. Although this program has been used by the USGS, no 
%warranty, expressed or implied, is made by the USGS or the United States 
%Government as to the accuracy and functioning of the program and related 
%program material nor shall the fact of distribution constitute any such 
%warranty, and no responsibility is assumed by the USGS in connection 
%therewith.

%Citation:  Butman, Bradford, Alexander, P. Soupy, Anderson, S.P., 
%Lightsom, F.L., Scotti, Alberto, and Beardsley, R.C., 2004, The 
%Massachusetts Bay Internal Wave Experiment, August 1998:  Data Report:  
%U.S. Geological Survey Data Series 85, 1 DVD-ROM.

filePath = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\DATAFILES\';

%Full Data Record
startTime = [1998 08 07 00 00 00];
endTime = [1998 09 02 00 00 00];
theStep = 1/24;
theTickSpace = 2;
tickCommand = 'datetick(''x'', 6, ''keepticks'', ''keeplimits'')';

% %August 19 - August 24
% startTime = [1998 08 19 00 00 00];
% endTime = [1998 08 25 00 00 00];
% theStep = 5/(24*60);
% theTickSpace = 1;
% tickCommand = 'datetick(''x'', 6, ''keepticks'', ''keeplimits'')';

% %August 22
% startTime = [1998 08 22 00 00 00];
% endTime = [1998 08 23 00 00 00];
% theStep = 1/(24*60);
% theTickSpace = 2/24;
% tickCommand = 'datetick(''x'', 15, ''keepticks'', ''keeplimits'')';

if strcmp(upper(theSite), 'B')
    theFiles = {'tpod4494-a.nc'; 'tpod4491-a.nc'; 'tpod4483-a.nc'; 'tpod3662-a.nc'; ...
            'mcat0684-a.nc'; '5383sc-a.cdf'; '5371tcp-a.nc'};
    theDepths = [0 10 20 30 45 70 85];
elseif strcmp(upper(theSite), 'C')
    theFiles = {'5331mc-a.cdf'; '5341mc-a.cdf'; '5342mc-a.cdf'; '5343sc-a.cdf'; '5344-a.nc'};
    theDepths = [10 21 31 41 45];
elseif strcmp(upper(theSite), 'A')
    theFiles = {'5391ts-a.nc'; '5392ts-a.nc'; '5393ts-a.nc'; '5394ts-a.nc'; ...
            '5395ts-a.nc'};
    theDepths = [10 20 30 40 50];
end

depthRange = [0:5:(max(theDepths)+5)];
timeRange = julian(startTime):theStep:julian(endTime);
tickTime = julian2datenum(floor(min(timeRange))):theTickSpace:...
    julian2datenum(ceil(max(timeRange)));

theTitle = ['Temperature, Site ' upper(theSite)];

tempMat = repmat(NaN, length(depthRange), length(timeRange));

for index = 1:length(theFiles)
    theFile = theFiles{index};
    
    ncID = netcdf(fullfile(filePath, theFile), 'nowrite');
    theTimes = singleJD(ncID{'time'}(:), ncID{'time2'});
    theTemp = ncID{'T_20'}(:);
    
    if isempty(theTemp)
        theTemp = ncID{'temp'}(:);
    end
    
    if isempty(theTemp)
        theTemp = ncID{'temperature'}(:);
    end
    if isempty(theTemp)
        warning(['No temperature data found in file ' theFile '.'])
    end
    
    theTemp = ridfill_nan(theTemp);

    hitDepth = value2Index(depthRange, theDepths(index), 5);
    tempMat(hitDepth, :) = interp1(theTimes, theTemp, timeRange);
end

    
goods = find(~isnan(nanmean(tempMat')));
iTemp = interp2(timeRange(:)', depthRange(goods), ...
    tempMat(goods, :), timeRange(:)', depthRange(:),'linear');

figure, orient landscape

%Do the temperature plotting
[c, h] = contourf(julian2datenum(timeRange(:))', depthRange, iTemp, [4:1:16]);
axis ij
set(gca, 'xtick', tickTime)
xlim([julian2datenum(min(timeRange)-(max(timeRange)-min(timeRange))*0.05) ...
        julian2datenum(max(timeRange)+(max(timeRange)-min(timeRange))*0.05)])
eval(tickCommand);
orient landscape
l = line(repmat(([min(julian2datenum(timeRange)) max(julian2datenum(timeRange))])', ...
    1, length(theDepths)), repmat(theDepths(:)', 2, 1));
shading flat
xlabel('Time')
ylabel('Depth (m)')
title(theTitle)
colorbar
ylim([max([min(theDepths)-5 0]) max(depthRange)])
%clabel(c, h, [5:5:25], 'labelspacing', 200)
