function adcp_contour(theSite)
%
%   adcp_contour(theSite)
%
%   Function to create a contour plot of the low passed wind stress and current
%   observations from the profiling current meters.

%Soupy Alexander, 6/23/04

%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.


% %Full Data Record
% startTime = [1998 08 07 00 00 00];
% endTime = [1998 09 01 00 00 00];
% timeVec = julian(startTime):1/24:julian(endTime);
% theTickSpace = 2;
% tickCommand = 'datetick(''x'', 6, ''keepticks'', ''keeplimits'')';
% theGo = 1;
% titleLabel = 'Full Record';
% theExt = 'all';

% August 19 - August 24
% startTime = [1998 08 19 00 00 00];
% endTime = [1998 08 25 00 00 00];
% timeVec = julian(startTime):5/(60*24):julian(endTime);
% theTickSpace = 1;
% tickCommand = 'datetick(''x'', 6, ''keepticks'', ''keeplimits'')';
% theGo = 1;
% titleLabel = '8/19 - 8/24'; 
% theExt = '0819_0824';

% 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'')';
theGo = 0;
titleLabel = '8/22';
theExt = '0822';

tickTime = julian2datenum(julian(startTime)):theTickSpace:...
    julian2datenum(julian(endTime));

plotDirectory = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\PLOTS\TIME_SERIES';

%Load in the data
if strcmp(lower(theSite), 'b')
    theFile = '5372adc-a_wave.nc';
elseif strcmp(lower(theSite), 'c')
    theFile = '5352adc-a_wave.nc';
end

ncID = netcdf(fullfile('C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\DATAFILES', ...
    theFile), 'nowrite');
theTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
theDepth = ncID{'depth'}(:);
theU = ncID{'u_1205'}(:);
theV = ncID{'v_1206'}(:);
theU(:,end-1:end) = NaN;
theV(:,end-1:end) = NaN;

[theTheta, theLength] = cart2pol(theU, theV);
newTheta = theTheta - 20*pi/180;
[newU, newV] = pol2cart(newTheta, theLength);


if theGo == 1
    for indexTime = 1:length(timeVec)-1;
        inTime = find((theTime >= timeVec(indexTime)) & ...
            (theTime < timeVec(indexTime+1)));
        thisTime(indexTime) = nanmean(theTime(inTime));
        newU2(:,indexTime) = nanmean(newU(inTime,:));
        newV2(:,indexTime) = nanmean(newV(inTime,:));
    end
else
    inTime = tind(theTime, startTime, endTime);
    thisTime = theTime(inTime);
    newU2 = (newU(inTime,:))';
    newV2 = (newV(inTime,:))';
end

figure, orient landscape

B = subplot(2,1,1);
gregTime = gregorian(theTime);
contourf(julian2datenum(thisTime), theDepth(1:end-2), newU2(1:end-2,:), 32)
caxis([-40 40])
axis ij
shading flat
xlim([julian2datenum(julian(startTime)) julian2datenum(julian(endTime))])
ylim([0 max(theDepth)])
set(gca, 'xtick', tickTime)
ylabel('Depth (m)')
title('70^o/250^o Component of Current')
box('on')
eval(tickCommand);

C = subplot(2,1,2);
gregTime = gregorian(theTime);
contourf(julian2datenum(thisTime), theDepth(1:end-2), newV2(1:end-2,:), 32)
caxis([-40 40])
axis ij
shading flat
xlim([julian2datenum(julian(startTime)) julian2datenum(julian(endTime))])
ylim([0 max(theDepth)])
set(gca, 'xtick', tickTime)
ylabel('Depth (m)')
title('-20^o/160^o Component of Current')
box('on')
eval(tickCommand);

D = colorbar('southoutside');

set(gcf, 'paperpositionmode', 'auto')
set(gcf, 'units', 'inches')
set(gcf, 'position', [1 1 10 7.5])

suptitle({'Profiling Current Meter Observations'; ...
    ['Site ' upper(theSite) ', ' titleLabel]})

Cposit = get(C, 'position');
Bposit = get(B, 'position');
Bposit(4) = Cposit(4);
set(B, 'position', Bposit);

Cposit(2) = Cposit(2) + 0.05;
set(C, 'position', Cposit)

print(fullfile(plotDirectory, ...
    [lower(theSite) '_contour_' theExt]), '-dpdf')

close