function adcp_contour(theYears, dataDirectory, plotDirectory, siteArchive, theSite)
%
%   adcp_contour(theYears, dataDirectory, plotDirectory, siteArchive, theSite)
%
%   Function to create a contour plot of the low passed wind stress and current
%   observations from the profiling current meters for theYears.
%   Files are output as SITE_contour_YYYY.

%Soupy Alexander, 10/29/03

%Load in the wind stress
[windTime, tauU, tauV] = mblt_windstress(dataDirectory);

%Load in the site A data
depthVec = [4:2:30];
timeVec = [julian(1996, 1, 1, 00):6/24:julian(max(theYears), 12, 31, 24)];
allU = repmat(NaN, length(depthVec), length(timeVec));
allV = repmat(NaN, length(depthVec), length(timeVec));

warning off MATLAB:divideByZero

for indexFile = 1:length(siteArchive)
    if rem(indexFile,10) == 0
        fprintf(['On file number ' num2str(indexFile, '%.0f') ' of ' ...
                num2str(length(siteArchive), '%.0f') '.\n']);
    end
    fileName = siteArchive(indexFile).fileList;
    if isempty(strfind(fileName, 'adc'))
        continue
    end
    currFileA = netcdf(fullfile(dataDirectory, siteArchive(indexFile).fileList), 'nowrite');
    fileDepth = currFileA{'depth'}(:);
    fileTime = singlejd(currFileA{'time'}(:), currFileA{'time2'}(:));
    fileU = currFileA{'u_1205'}(:);
    fileV = currFileA{'v_1206'}(:);
    bads = find(fileU > 1e3 | fileV > 1e3);
    fileU(bads) = NaN;
    fileV(bads) = NaN;
    
    for indexDepth = 1:length(depthVec)
        hitDepth = value2index(fileDepth, depthVec(indexDepth), 2.1);
        
        if isnan(hitDepth)
            continue
        end
        
        fileUHit = fileU(:,hitDepth);
        fileVHit = fileV(:,hitDepth);
        bads = find(isnan(fileUHit + fileVHit));
        
        fileTimeI = fileTime;
        fileTimeI(bads) = [];
        fileUHit(bads) = [];
        fileVHit(bads) = [];
        
        if isempty(fileTimeI)
            continue
        end
        
        try
            [fileIU, fileTime_lp] = plfilt(fileUHit, fileTimeI, 6);
            [fileIV, fileTime_lp] = plfilt(fileVHit, fileTimeI, 6);
        catch
            continue
        end
        
        timeStart = max(find(timeVec <= min(fileTime_lp)));
        timeEnd = min(find(timeVec >= max(fileTime_lp)));
        
        fileIU2 = smart_interp(fileTime_lp, fileIU, timeVec(timeStart:timeEnd), 12);
        fileIV2 = smart_interp(fileTime_lp, fileIV, timeVec(timeStart:timeEnd), 12);
        
        allU(indexDepth, timeStart:timeEnd) = fileIU2;
        allV(indexDepth, timeStart:timeEnd) = fileIV2;
    end    
end  

for indexYear = 1:length(theYears)
    
    targetYear = theYears(indexYear);
    if (targetYear < 1996 & strcmp(lower(theSite), 'a')) | ...
            (targetYear < 1997 & strcmp(lower(theSite), 'b'))
        continue
    end
    
    figure, orient landscape
    
    xTick = [];
    for indexMonth = 1:12
        xTick = [xTick(:); datenum(targetYear, indexMonth, 1)];
    end
    
    %Do the plots
    
    hitTime = windTime;
    hitU = tauU;
    hitV = tauV;
    
    inYear = find((hitTime >= julian(targetYear, 1, 1, 00)) & ...
        hitTime <= julian(targetYear, 12, 31, 24));
    hitTime = hitTime(inYear); 
    
    [theAngleRad, theLength] = cart2pol(hitU(inYear), hitV(inYear));
    theAngle = theAngleRad * 180/pi;
    
    goods = find(~isnan(theAngle));
    theAngle = theAngle(goods);
    theLength = theLength(goods);
    
    plotX = julian2datenum(hitTime(goods));
    plotY = zeros(size(plotX));
    
    A = subplot(3,1,1);
    sticksafe(plotX, plotY, theLength, theAngle)
    ylim([-4 4])
    set(gca, 'xtick', xTick, 'ytick', [-4:2:4])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    ylabel('dynes/cm^2')
    title('Low-Passed Wind Stress')
    l = line([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)], [0 0]);
    set(l, 'color', 'k')
    box('on')
    sticksafe
    
    B = subplot(3,1,2);
    gregTime = gregorian(timeVec);
    inTime = find(gregTime(:,1) == targetYear);
    contourf(julian2datenum(timeVec(inTime)), depthVec, allU(:,inTime))
    caxis([-40 40])
    if strcmp(lower(theSite), 'a')
        ylim([0 28])
    else
        ylim([0 18])
    end
    axis ij
    shading flat
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    set(gca, 'xtick', xTick)
    ylabel('Depth (m)')
    title('U Component of Current')
    box('on')
    datetick('x', 3, 'keeplimits', 'keepticks')
    
    C = subplot(3,1,3);
    gregTime = gregorian(timeVec);
    inTime = find(gregTime(:,1) == targetYear);
    contourf(julian2datenum(timeVec(inTime)), depthVec, allV(:,inTime))
    caxis([-40 40])
    if strcmp(lower(theSite), 'a')
        ylim([0 28])
    else
        ylim([0 18])
    end
    axis ij
    shading flat
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    set(gca, 'xtick', xTick)
    ylabel('Depth (m)')
    title('V Component of Current')
    box('on')
    datetick('x', 3, 'keeplimits', 'keepticks')
    D = colorbar('horiz');
    
    Dparam = get(D, 'position');
    dWidth = Dparam(4);
    
    Bparam = get(B, 'position');
    plotHeight = Bparam(4);
    Bparam(2) = Bparam(2) + dWidth/2;
    set(B, 'position', Bparam)
    
    Cparam = get(C, 'position');
    Cparam(4) = plotHeight;
    Cparam(2) = Cparam(2) - 1.25*dWidth;
    set(C, 'position', Cparam);
    
    Dparam(2) = 0.05;
    set(D, 'position', Dparam)
    
    suptitle({'Low-Pass Filtered Wind Stress and Profiling Current Meter Observations'; ...
            ['Site ' upper(theSite) ', ' num2str(targetYear, '%.0f')]})
    
    print(fullfile(plotDirectory, ...
        [lower(theSite) '_contour_' num2str(targetYear, '%.0f')]), '-dpdf')
    
    close
    
end





