function monthly_stats2(dataDirectory, plotDirectory, siteA, siteB)
%
%   monthly_stats2(dataDirectory, siteA, siteB)
%
%   Code to create the plots/tables of monthly statistics for the
%   meteorological and oceanographic observations.  A continuation of
%   monthly_stats.m, which was starting to get unwieldy; this does the
%   calculations for wind stress and ADCP current.

%Soupy Alexander, 11/5/03

%Search the data archives to establish the time series of wind stress and
%current

currTime = [];
currTime_lp = [];
c5 = [];
c10 = [];
c1 = [];
c5_lp = [];
c10_lp = [];
c1_lp = [];

%Load in the site A data
for indexFile = 1:length(siteA)
    fileName = siteA(indexFile).fileList;
    if isempty(strfind(fileName, 'adc'))
        continue
    end
    
    currFileA = netcdf(fullfile(dataDirectory, siteA(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;
    
    hit5 = value2Index(fileDepth, 5);
    cU5 = fileU(:, hit5);
    try
        [fcU5_lp, fileTime_lp] = plfilt(cU5, fileTime, 6); 
    catch
        continue
    end
    cV5 = fileV(:, hit5);
    [fcV5_lp, fileTime_lp] = plfilt(cV5, fileTime, 6);
    fc5 = cU5 + i.*cV5;
    fc5_lp = fcU5_lp + i.*fcV5_lp;
    
    hit10 = value2Index(fileDepth, 15);
    cU10 = fileU(:, hit10);
    [fcU10_lp, fileTime_lp] = plfilt(cU10, fileTime, 6); 
    cV10 = fileV(:, hit10);
    [fcV10_lp, fileTime_lp] = plfilt(cV10, fileTime, 6);
    fc10 = cU10 + i.*cV10;
    fc10_lp = fcU10_lp + i.*fcV10_lp;
    
    hit1 = find(fileDepth == max(fileDepth));
    cU1 = fileU(:, hit1);
    [fcU1_lp, fileTime_lp] = plfilt(cU1, fileTime, 6); 
    cV1 = fileV(:, hit1);
    [fcV1_lp, fileTime_lp] = plfilt(cV1, fileTime, 6);
    fc1 = cU1 + i.*cV1;
    fc1_lp = fcU1_lp + i.*fcV1_lp;
    
    [currTime] = [currTime(:); fileTime(:)];
    [currTime_lp] = [currTime_lp(:); fileTime_lp(:)];
    [c5] = [c5(:); fc5(:)];
    [c10] = [c10(:); fc10(:)];
    [c1] = [c1(:); fc1(:)];
    [c5_lp] = [c5_lp(:); fc5_lp(:)];
    [c10_lp] = [c10_lp(:); fc10_lp(:)];
    [c1_lp] = [c1_lp(:); fc1_lp(:)];
end

BcurrTime = [];
BcurrTime_lp = [];
Bc5 = [];
Bc10 = [];
Bc1 = [];
Bc5_lp = [];
Bc10_lp = [];
Bc1_lp = [];

%Load in the site B data
for indexFile = 1:length(siteB)
    fileName = siteB(indexFile).fileList;
    if isempty(strfind(fileName, 'adc'))
        continue
    end
    
    currFileA = netcdf(fullfile(dataDirectory, siteB(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;
    
    hit5 = value2Index(fileDepth, 5);
    cU5 = fileU(:, hit5);
    try
        [fcU5_lp, fileTime_lp] = plfilt(cU5, fileTime, 6); 
    catch
        continue
    end
    cV5 = fileV(:, hit5);
    [fcV5_lp, fileTime_lp] = plfilt(cV5, fileTime, 6);
    fc5 = cU5 + i.*cV5;
    fc5_lp = fcU5_lp + i.*fcV5_lp;
    
    hit10 = value2Index(fileDepth, 10);
    cU10 = fileU(:, hit10);
    [fcU10_lp, fileTime_lp] = plfilt(cU10, fileTime, 6); 
    cV10 = fileV(:, hit10);
    [fcV10_lp, fileTime_lp] = plfilt(cV10, fileTime, 6);
    fc10 = cU10 + i.*cV10;
    fc10_lp = fcU10_lp + i.*fcV10_lp;
    
    hit1 = find(fileDepth == max(fileDepth));
    cU1 = fileU(:, hit1);
    [fcU1_lp, fileTime_lp] = plfilt(cU1, fileTime, 6); 
    cV1 = fileV(:, hit1);
    [fcV1_lp, fileTime_lp] = plfilt(cV1, fileTime, 6);
    fc1 = cU1 + i.*cV1;
    fc1_lp = fcU1_lp + i.*fcV1_lp;
    
    [BcurrTime] = [currTime(:); fileTime(:)];
    [BcurrTime_lp] = [currTime_lp(:); fileTime_lp(:)];
    [Bc5] = [Bc5(:); fc5(:)];
    [Bc10] = [Bc10(:); fc10(:)];
    [Bc1] = [Bc1(:); fc1(:)];
    [Bc5_lp] = [Bc5_lp(:); fc5_lp(:)];
    [Bc10_lp] = [Bc10_lp(:); fc10_lp(:)];
    [Bc1_lp] = [Bc1_lp(:); fc1_lp(:)];
end

%Wind stress data
windTime = [];
windTime_lp = [];
w = [];
w_lp = [];

%Load in the wind data
windFile = netcdf(fullfile(dataDirectory, 'lt_wind.nc'), 'nowrite');
windTime = singlejd(windFile{'time'}(:), windFile{'time2'}(:));
windU = windFile{'WU_422'}(:);
windV = windFile{'WV_423'}(:);

%Estimate the 10 m wind speed (buoy 44013 is at 5 m)
[tauX, tauY, windU_10, windV_10] = wstress(windU, windV, 5);
[tauU, tauV] = wstress(windU_10, windV_10);
w = tauU + i.*tauV;

%Set up the new time base & interpolate the data
newTime = min(windTime):1/24:max(windTime);
[windU_i] = smart_interp(windTime, windU, newTime, 5);
[windV_i] = smart_interp(windTime, windV, newTime, 5);
[windU_lp, windTime_lp] = plfilt(windU_i, newTime, 6);
[windV_lp] = plfilt(windV_i, newTime, 6);

%Estimate the 10 m wind speed (buoy 44013 is at 5 m)
[tauX, tauY, windU_10, windV_10] = wstress(windU_lp, windV_lp, 5);
[tauU, tauV] = wstress(windU_10, windV_10);
w_lp = tauU + i.*tauV;

clear windV_lp windV_i windV_10 windV windU_lp windU_i windU_10 windU
clear windFile tauY tauX tauV tauU newTime mbltArchive indexFile hit5
clear hit10 hit1 fileU fileV fileTime_lp fileTime fileName fileDepth
clear fcV5_lp fcV1_lp fcV10_lp fcU5_lp fcU1_lp fcU10_lp fc5_lp fc1_lp fc10_lp
clear fc1 fc5 fc10 currFileA cV5 cV10 cV1 cU5 cU10 cU1 bads

timeVars = {'currTime'; 'currTime'; 'currTime'; 'BcurrTime'; 'BcurrTime'; ...
        'BcurrTime'; 'windTime'};
theVars = {'c5'; 'c10'; 'c1'; 'Bc5'; 'Bc10'; 'Bc1'; 'w'};
fileNames = {'a_c5a'; 'a_c15a'; 'a_c1a'; 'b_c5a'; 'b_c10a'; 'b_c1a'; 'stress'};
varNames = {'Site A, ADCP Current Speed, 5 mbs'; 'Site A, ADCP Current Speed, 15 mbs'; ...
        'Site A, ADCP Current Speed, Deepest Bin'; 'Site B, ADCP Current Speed, 5 mbs'; ...
        'Site B, ADCP Current Speed, 10 mbs'; 'Site B, ADCP Current Speed, Deepest Bin'; 'Wind Stress'};
theUnits = {'cm/s'; 'cm/s'; 'cm/s'; 'cm/s'; 'cm/s'; 'cm/s'; 'dynes/cm^2'};
theLimits = {[0 20]; [0 20]; [0 15]; [0 20]; [0 15]; [0 15]; [0 3]};
theLimits2 = {[0 10]; [0 10]; [0 6]; [0 10]; [0 10]; [0 6]; [0 3]}; 
monthList = {'Jan'; 'Feb'; 'Mar'; 'Apr'; 'May'; 'Jun'; 'Jul'; 'Aug'; 'Sep'; 'Oct'; ...
        'Nov'; 'Dec'};

for indexFile = 1:length(timeVars)
    eval(['theTime = ' timeVars{indexFile} ';']);
    eval(['xData = real(' theVars{indexFile} ');']);
    eval(['yData = imag(' theVars{indexFile} ');']);
        
    gregTime = gregorian(theTime);
    post1990 = find(gregTime(:,1) >=1990);
    theTime = theTime(post1990);
    gregTime = gregTime(post1990, :);
    xData = xData(post1990);
    yData = yData(post1990);
    good = find(~isnan(xData + yData));
    xData = xData(good);
    yData = yData(good);
    theTime = theTime(good);
    gregTime = gregTime(good,:);
    
    eval(['theTimeLP = ' timeVars{indexFile} '_lp;']);
    eval(['xDataLP = real(' theVars{indexFile} '_lp);']);
    eval(['yDataLP = imag(' theVars{indexFile} '_lp);']);
    
    theData = sqrt(xData.^2 + yData.^2);
    theDataLP = sqrt(xDataLP.^2 + yDataLP.^2);
    gregTimeLP = gregorian(theTimeLP);
    
    yearList = min(gregTime(:,1)):max(gregTime(:,1));
    
    %Figure out the monthly means
    for indexMonth = 1:12
        inMonth = find(gregTime(:,2) == indexMonth);
        inMonthLP = find(gregTimeLP(:,2) == indexMonth);
       
        if isempty(inMonth)
            monthPointsAll(indexMonth) = 0;
            continue
        end
        
        eastMeanAll(indexMonth) = nanmean(xData(inMonth));
        eastSTDAll(indexMonth) = nanstd(xData(inMonth));
        eastSTD_lpAll(indexMonth) = nanstd(xDataLP(inMonthLP));
        
        northMeanAll(indexMonth) = nanmean(yData(inMonth));
        northSTDAll(indexMonth) = nanstd(yData(inMonth));
        northSTD_lpAll(indexMonth) = nanstd(yDataLP(inMonthLP));
        
        [vecDirI, vecSpdI] = cart2pol(eastMeanAll(indexMonth), northMeanAll(indexMonth));
        vecSpdAll(indexMonth) = vecSpdI;
        vecDirAll(indexMonth) = anglextocompass(vecDirI*180/pi);
        
        monthMeanAll(indexMonth) = nanmean(theData(inMonth));
        monthAnomAll(indexMonth) = 0;
        monthSTDAll(indexMonth) = nanstd(theData(inMonth));
        monthSTD_lpAll(indexMonth) = nanstd(theDataLP(inMonthLP));
        monthMinAll(indexMonth) = nanmin(theData(inMonth));
        monthMaxAll(indexMonth) = nanmax(theData(inMonth));
        monthPointsAll(indexMonth) = length(find(~isnan(theData(inMonth))));
    end
    
    yearMeanAll = nanmean(theData);
    
    eastMeanAll(13) = nanmean(xData);
    eastSTDAll(13) = nanstd(xData);
    eastSTD_lpAll(13) = nanstd(xDataLP);
    
    northMeanAll(13) = nanmean(yData);
    northSTDAll(13) = nanstd(yData);
    northSTD_lpAll(13) = nanstd(yDataLP);
    
    [vecDirI, vecSpdI] = cart2pol(eastMeanAll(13), northMeanAll(13));
    vecSpdAll(13) = vecSpdI;
    vecDirAll(13) = anglextocompass(vecDirI*180/pi);
    
    monthMeanAll(13) = nanmean(theData);
    monthAnomAll(13) = 0;
    monthSTDAll(13) = nanstd(theData);
    monthSTD_lpAll(13) = nanstd(theDataLP);
    monthMinAll(13) = nanmin(theData);
    monthMaxAll(13) = nanmax(theData);
    monthPointsAll(13) = length(find(~isnan(theData)));
    
    eastMean = repmat(NaN, length(yearList), 13);
    eastSTD = repmat(NaN, length(yearList), 13);
    eastSTD_lp = repmat(NaN, length(yearList), 13);
    northMean = repmat(NaN, length(yearList), 13);
    northSTD = repmat(NaN, length(yearList), 13);
    northSTD_lp = repmat(NaN, length(yearList), 13);    
    vecSpd = repmat(NaN, length(yearList), 13);
    vecDir = repmat(NaN, length(yearList), 13);
    monthMean = repmat(NaN, length(yearList), 13);
    monthAnom = repmat(NaN, length(yearList), 13);
    monthSTD = repmat(NaN, length(yearList), 13);
    monthSTD_lp = repmat(NaN, length(yearList), 13);
    monthMin = repmat(NaN, length(yearList), 13);
    monthMax = repmat(NaN, length(yearList), 13);
    monthPoints = repmat(NaN, length(yearList), 13);
    
    
    %Do the monthly averages
    for indexYear = 1:length(yearList)
        targetYear = yearList(indexYear);
        
        inYear = find(gregTime(:,1) == targetYear);
        
        if isempty(inYear)
            monthPoints(indexYear, 13) = 0;
            continue
        end
        
        inYearLP = find(gregTimeLP(:,1) == targetYear);
        
        timeInYear = theTime(inYear);
        gregTimeInYear = gregorian(timeInYear);
        dataInYear = theData(inYear);
        timeInYearLP = theTimeLP(inYearLP);
        gregTimeInYearLP = gregorian(timeInYearLP);
        dataInYearLP = theDataLP(inYearLP);
        
        xDataInYear = xData(inYear);
        xDataInYearLP = xDataLP(inYearLP);
        yDataInYear = yData(inYear);
        yDataInYearLP = yDataLP(inYearLP);
        
        eastMean(indexYear, 13) = nanmean(xDataInYear);
        eastSTD(indexYear, 13) = nanstd(xDataInYear);
        eastSTD_lp(indexYear, 13) = nanstd(xDataInYearLP);
        
        northMean(indexYear, 13) = nanmean(yDataInYear);
        northSTD(indexYear, 13) = nanstd(yDataInYear);
        northSTD_lp(indexYear, 13) = nanstd(yDataInYearLP);
        
        [vecDirI, vecSpdI] = cart2pol(eastMean(indexYear, 13), northMean(indexYear, 13));
        vecSpd(indexYear, 13) = vecSpdI;
        vecDir(indexYear, 13) = anglextocompass(vecDirI*180/pi);
        
        monthMean(indexYear, 13) = nanmean(dataInYear);
        monthAnom(indexYear, 13) = nanmean(dataInYear) - yearMeanAll;
        monthSTD(indexYear, 13) = nanstd(dataInYear);
        monthSTD_lp(indexYear, 13) = nanstd(dataInYearLP);
        monthMin(indexYear, 13) = nanmin(dataInYear);
        monthMax(indexYear, 13) = nanmax(dataInYear);
        monthPoints(indexYear, 13) = length(find(~isnan(dataInYear)));
        
        for indexMonth = 1:12
            inMonth = find(gregTimeInYear(:,2) == indexMonth);
            if isempty(inMonth) | (length(inMonth) < 6)
                monthPoints(indexYear, indexMonth) = 0;
                continue
            end
            inMonthLP = find(gregTimeInYearLP(:,2) == indexMonth);
            
            eastMean(indexYear, indexMonth) = nanmean(xDataInYear(inMonth));
            eastSTD(indexYear, indexMonth) = nanstd(xDataInYear(inMonth));
            eastSTD_lp(indexYear, indexMonth) = nanstd(xDataInYearLP(inMonthLP));
            
            northMean(indexYear, indexMonth) = nanmean(yDataInYear(inMonth));
            northSTD(indexYear, indexMonth) = nanstd(yDataInYear(inMonth));
            northSTD_lp(indexYear, indexMonth) = nanstd(yDataInYearLP(inMonthLP));
            
            [vecDirI, vecSpdI] = cart2pol(eastMean(indexYear, indexMonth), ...
                northMean(indexYear, indexMonth));
            vecSpd(indexYear, indexMonth) = vecSpdI;
            vecDir(indexYear, indexMonth) = anglextocompass(vecDirI*180/pi);
            
            monthMean(indexYear,indexMonth) = nanmean(dataInYear(inMonth));
            monthAnom(indexYear, indexMonth) = nanmean(dataInYear(inMonth)) - ...
                monthMeanAll(indexMonth);
            monthSTD(indexYear, indexMonth) = nanstd(dataInYear(inMonth));
            monthSTD_lp(indexYear, indexMonth) = nanstd(dataInYearLP(inMonthLP));
            monthMin(indexYear, indexMonth) = nanmin(dataInYear(inMonth));
            monthMax(indexYear, indexMonth) = nanmax(dataInYear(inMonth));
            monthPoints(indexYear, indexMonth) = length(find(~isnan(dataInYear(inMonth))));
        end
    end
    
    %Do the actual table creation/plotting
    fID = fopen(fullfile(plotDirectory, [fileNames{indexFile} '.txt']), 'w');
    fprintf(fID, [varNames{indexFile} '\n\n']);
    spaces = '\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t';
    titles = '# Points\tU Mean\tU STD (H)\tU STD(L)\tV Mean\tV STD (H)\tV STD (L)\t';
    titles = [titles 'Vec. Speed\tVec. Dir.\tSc. Mean\tSc. Anom.\tSc. STD (H)\tSc.STD(L)\t'];
    titles = [titles 'Sc. Min\tSc. Max\t'];
    fprintf(fID, ['\t' repmat('Jan\t', 1, 15) repmat('Feb\t', 1, 15) ...
            repmat('Mar\t', 1, 15) repmat('Apr\t', 1, 15) repmat('May\t', 1, 15) ...
            repmat('Jun\t', 1, 15) repmat('Jul\t', 1, 15) repmat('Aug\t', 1, 15) ...
            repmat('Sep\t', 1, 15) repmat('Oct\t', 1, 15) repmat('Nov\t', 1, 15) ...
            repmat('Dec\t', 1, 15) repmat('Year\t', 1, 15) '\n']);
    fprintf(fID, ['\t' repmat(titles, 1, 13) '\n']);
    
    f = figure; orient landscape
    f2 = figure; orient landscape
    colList = {'b'; 'g'; 'y'; 'm'; 'k'; 'c'};
    colList = repmat(colList(:), 10, 1);
    symList = {'.'; 'o'; 'x'; '+'; 's'; 'd'; 'v'; 'p'};
    symList = repmat(symList(:), 10, 1);
    
    figure(f)
    A = subplot(2,1,1); hold on
    B = subplot(2,1,2); hold on
    
    figure(f2)
    A2 = subplot(2,1,1); hold on
    B2 = subplot(2,1,2); hold on
    
    for indexYear = 1:length(yearList)
        
        if monthPoints(indexYear, 13) == 0
            continue
        end
        
        fprintf(fID, [num2str(yearList(indexYear), '%.0f') '\t']);
        
        for indexMonth = 1:13
            fprintf(fID, [num2str(monthPoints(indexYear, indexMonth), '%.0f') '\t']);
            if monthPoints(indexYear, indexMonth) ~= 0
                fprintf(fID, [num2str(eastMean(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(eastSTD(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(eastSTD_lp(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(northMean(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(northSTD(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(northSTD_lp(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(vecSpd(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(vecDir(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(monthMean(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(monthAnom(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(monthSTD(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(monthSTD_lp(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(monthMin(indexYear, indexMonth), '%.2f') '\t' ...
                        num2str(monthMax(indexYear, indexMonth), '%.2f') '\t']);
            else
                fprintf(fID, '---\t---\t---\t---\t---\t---\t---\t---\t---\t---\t---\t---\t---\t---\t');
            end
            
            if indexMonth ~= 13   
                figure(f)
                subplot(A)
                plot(indexMonth, monthMean(indexYear, indexMonth), ...
                    [colList{indexYear} symList{indexYear}]);
                subplot(B)
                plot(indexMonth, monthAnom(indexYear, indexMonth), ...
                    [colList{indexYear} symList{indexYear}]);
                
                figure(f2)
                subplot(A2)
                plot(indexMonth, vecSpd(indexYear, indexMonth), ...
                    [colList{indexYear} symList{indexYear}]);
                subplot(B2)
                plot(indexMonth, vecDir(indexYear, indexMonth), ...
                    [colList{indexYear} symList{indexYear}]);
            end
        end
        fprintf(fID, '\n');
    end
    
    
    fprintf(fID, 'All Years\t');
    
    for indexMonth = 1:13
        fprintf(fID, [num2str(monthPointsAll(indexMonth), '%.0f') '\t']);
        if monthPointsAll(indexMonth) ~= 0
            fprintf(fID, [num2str(eastMeanAll(indexMonth), '%.2f') '\t' ...
                    num2str(eastSTDAll(indexMonth), '%.2f') '\t' ...
                    num2str(eastSTD_lpAll(indexMonth), '%.2f') '\t' ...
                    num2str(northMeanAll(indexMonth), '%.2f') '\t' ...
                    num2str(northSTDAll(indexMonth), '%.2f') '\t' ...
                    num2str(northSTD_lpAll(indexMonth), '%.2f') '\t' ...
                    num2str(vecSpdAll(indexMonth), '%.2f') '\t' ...
                    num2str(vecDirAll(indexMonth), '%.0f') '\t' ...
                    num2str(monthMeanAll(indexMonth), '%.2f') '\t' ...
                    num2str(monthAnomAll(indexMonth), '%.2f') '\t' ...
                    num2str(monthSTDAll(indexMonth), '%.2f') '\t' ...
                    num2str(monthSTD_lpAll(indexMonth), '%.2f') '\t' ...
                    num2str(monthMinAll(indexMonth), '%.2f') '\t' ...
                    num2str(monthMaxAll(indexMonth), '%.2f') '\t']);
        else
            fprintf(fID, '---\t---\t---\t---\t---\t---\t---\t---\t---\t---\t---\t---\t---\t---\t');
        end
        
        if indexMonth ~= 13
            figure(f)
            subplot(A)
            plot(indexMonth, monthMeanAll(indexMonth), 'r*');
            figure(f2)
            subplot(A2)
            plot(indexMonth, vecSpdAll(indexMonth), 'r*');
            subplot(B2)
            plot(indexMonth, vecDirAll(indexMonth), 'r*');
        end
    end
    
    fclose(fID);
    
    figure(f)
    subplot(A)
    xlim([0 13])
    ylim(theLimits{indexFile})
    set(gca, 'xtick', [1:12], 'xticklabel', monthList)
    ylabel(theUnits{indexFile})
    theSize = get(gca, 'position');
    theSize(3) = theSize(3)*.89;
    set(gca, 'position', theSize);
    title('Monthly Mean')
    box('on')
    
    subplot(B)
    xlim([0 13])
    ylim([-1*max(abs(ylim)) max(abs(ylim))])
    set(gca, 'xtick', [1:12], 'xticklabel', monthList)
    ylabel(theUnits{indexFile})
    theSize1 = get(gca, 'position');
    theSize1(3) = theSize1(3)*.89;
    set(gca, 'position', theSize1);
    title('Monthly Anomaly')
    box('on')
    
    legSize(1) = theSize1(1) + theSize1(3)*1.01;
    legSize(2) = theSize1(2);
    legSize(3) = theSize1(3)*0.09;
    legSize(4) = (theSize(2)+theSize(4)) - theSize1(2);
    theLeg = subplot('position', legSize);
    hold on
    
    for index = 1:length(yearList)
        plot(1, -1* yearList(index), [colList{index} symList{index}]);
        t = text(3, -1* yearList(index), num2str(yearList(index), '%.0f'));
        set(t, 'fontsize', 8);
    end
    plot(1, -1*(max(yearList)+1), 'r*')
    t = text(3, -1*(max(yearList)+1), 'All Years');
    set(t, 'fontsize', 8);
    axis('off')
    box('off')
    xlim([0 5])
    ylim([-1*(max(yearList)+2) -1*(min(yearList)-1)])
    
    suptitle(varNames{indexFile})
    
    print('-dpdf', fullfile(plotDirectory, fileNames{indexFile}))
    close
    
    
    figure(f2)
    subplot(A2)
    xlim([0 13])
    ylim(theLimits2{indexFile})
    set(gca, 'xtick', [1:12], 'xticklabel', monthList)
    ylabel(theUnits{indexFile})
    theSize = get(gca, 'position');
    theSize(3) = theSize(3)*.89;
    set(gca, 'position', theSize);
    title('Vector Mean Speed')
    box('on')
    
    subplot(B2)
    xlim([0 13])
    ylim([0 360])
    set(gca, 'xtick', [1:12], 'xticklabel', monthList)
    ylabel('Degrees CW from North')
    theSize1 = get(gca, 'position');
    theSize1(3) = theSize1(3)*.89;
    set(gca, 'position', theSize1);
    title('Vector Mean Direction')
    box('on')
    
    legSize(1) = theSize1(1) + theSize1(3)*1.01;
    legSize(2) = theSize1(2);
    legSize(3) = theSize1(3)*0.09;
    legSize(4) = (theSize(2)+theSize(4)) - theSize1(2);
    theLeg = subplot('position', legSize);
    hold on
    
    for index = 1:length(yearList)
        plot(1, -1* yearList(index), [colList{index} symList{index}]);
        t = text(3, -1* yearList(index), num2str(yearList(index), '%.0f'));
        set(t, 'fontsize', 8);
    end
    plot(1, -1*(max(yearList)+1), 'r*')
    t = text(3, -1*(max(yearList)+1), 'All Years');
    set(t, 'fontsize', 8);
    axis('off')
    box('off')
    xlim([0 5])
    ylim([-1*(max(yearList)+2) -1*(min(yearList)-1)])
    
    suptitle(varNames{indexFile})
    
    print('-dpdf', fullfile(plotDirectory, ['p2_' fileNames{indexFile}]))
    close
end