function river_ts(theYears, dataDirectory, plotDirectory)
%
%   river_ts(theYears, dataDirectory, plotDirectory
%
%   Plots of the streamflow from the Charles, Parker, Ipswich, and
%   Merrimack rivers and the temperature and salinity from various depths
%   at sites A and B.  Plots are named SITE_river_YYYY.

%Soupy Alexander, 10/29/03

%Load in the river data
charlesData = xlsread(fullfile(dataDirectory, 'charles_modified.xls'));
charlesTime = charlesData(:,1) + datenum('30-Dec-1899');    %Convert from Excel
charlesFlow = charlesData(:,2);
[charlesTime, charlesFlow] = nanfill_2(charlesTime, charlesFlow);

parkerData = xlsread(fullfile(dataDirectory, 'parker_modified.xls'));
parkerTime = parkerData(:,1) + datenum('30-Dec-1899');    %Convert from Excel
parkerFlow = parkerData(:,2);
[parkerTime, parkerFlow] = nanfill_2(parkerTime, parkerFlow);

ipswichData = xlsread(fullfile(dataDirectory, 'ipswich_modified.xls'));
ipswichTime = ipswichData(:,1) + datenum('30-Dec-1899');    %Convert from Excel
ipswichFlow = ipswichData(:,2);
[ipswichTime, ipswichFlow] = nanfill_2(ipswichTime, ipswichFlow);

merrimackData = xlsread(fullfile(dataDirectory, 'merrimack_modified.xls'));
merrimackTime = merrimackData(:,1) + datenum('30-Dec-1899');    %Convert from Excel
merrimackFlow = merrimackData(:,2);
[merrimackTime, merrimackFlow] = nanfill_2(merrimackTime, merrimackFlow);

%Load in the temperature data
tempFile1 = netcdf(fullfile(dataDirectory, 'lt_t_0m.nc'), 'nowrite');
tempTime1 = singlejd(tempFile1{'time'}(:), tempFile1{'time2'}(:));
tempData1 = tempFile1{'T_25'}(:);
[tempTime1, tempData1] = nanfill_2(tempTime1, tempData1);

tempFile2 = netcdf(fullfile(dataDirectory, 'a_temp_5mbs.nc'), 'nowrite');
tempTime2 = singlejd(tempFile2{'time'}(:), tempFile2{'time2'}(:));
tempData2 = tempFile2{'T_20'}(:);
[tempTime2, tempData2] = nanfill_2(tempTime2, tempData2);

tempFile3 = netcdf(fullfile(dataDirectory, 'a_temp_10mab.nc'), 'nowrite');
tempTime3 = singlejd(tempFile3{'time'}(:), tempFile3{'time2'}(:));
tempData3 = tempFile3{'T_20'}(:);
[tempTime3, tempData3] = nanfill_2(tempTime3, tempData3);

tempFile4 = netcdf(fullfile(dataDirectory, 'a_temp_1mab.nc'), 'nowrite');
tempTime4 = singlejd(tempFile4{'time'}(:), tempFile4{'time2'}(:));
tempData4 = tempFile4{'T_20'}(:);
[tempTime4, tempData4] = nanfill_2(tempTime4, tempData4);

tempFile3_b = netcdf(fullfile(dataDirectory, 'b_temp_10mab.nc'), 'nowrite');
tempTime3_b = singlejd(tempFile3_b{'time'}(:), tempFile3_b{'time2'}(:));
tempData3_b = tempFile3_b{'T_20'}(:);
[tempTime3_b, tempData3_b] = nanfill_2(tempTime3_b, tempData3_b);

tempFile4_b = netcdf(fullfile(dataDirectory, 'b_temp_1mab.nc'), 'nowrite');
tempTime4_b = singlejd(tempFile4_b{'time'}(:), tempFile4_b{'time2'}(:));
tempData4_b = tempFile4_b{'T_20'}(:);
[tempTime4_b, tempData4_b] = nanfill_2(tempTime4_b, tempData4_b);

%Load in the salinity data
salFile2 = netcdf(fullfile(dataDirectory, 'a_sal_5mbs.nc'), 'nowrite');
salTime2 = singlejd(salFile2{'time'}(:), salFile2{'time2'}(:));
salData2 = salFile2{'S_40'}(:);
[salTime2, salData2] = nanfill_2(salTime2, salData2);

salFile3 = netcdf(fullfile(dataDirectory, 'a_sal_10mab.nc'), 'nowrite');
salTime3 = singlejd(salFile3{'time'}(:), salFile3{'time2'}(:));
salData3 = salFile3{'S_40'}(:);
[salTime3, salData3] = nanfill_2(salTime3, salData3);

salFile4 = netcdf(fullfile(dataDirectory, 'a_sal_1mab.nc'), 'nowrite');
salTime4 = singlejd(salFile4{'time'}(:), salFile4{'time2'}(:));
salData4 = salFile4{'S_40'}(:);
[salTime4, salData4] = nanfill_2(salTime4, salData4);

salFile3_b = netcdf(fullfile(dataDirectory, 'b_sal_10mab.nc'), 'nowrite');
salTime3_b = singlejd(salFile3_b{'time'}(:), salFile3_b{'time2'}(:));
salData3_b = salFile3_b{'S_40'}(:);
[salTime3_b, salData3_b] = nanfill_2(salTime3_b, salData3_b);

salFile4_b = netcdf(fullfile(dataDirectory, 'b_sal_1mab.nc'), 'nowrite');
salTime4_b = singlejd(salFile4_b{'time'}(:), salFile4_b{'time2'}(:));
salData4_b = salFile4_b{'S_40'}(:);
[salTime4_b, salData4_b] = nanfill_2(salTime4_b, salData4_b);

for indexYear = 1:length(theYears)
    targetYear = theYears(indexYear);
    figure, orient landscape
    
    xTick = [];
    for indexMonth = 1:12
        xTick = [xTick(:); datenum(targetYear, indexMonth, 1)];
    end
       
    subplot(4,1,1)
    plot(charlesTime, charlesFlow)
    hold on
    plot(ipswichTime, ipswichFlow, 'r')
    plot(parkerTime, parkerFlow, 'g')
    ylabel('cfs')
    ylim([0 3000])
    title('Streamflow')
    set(gca, 'xtick', xTick, 'ytick', [0:500:3000])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])  
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.9*diff(ylim) min(ylim)+.9*diff(ylim)]);
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.9*diff(ylim), 'Charles');
    set(t, 'fontsize', 8);
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.8*diff(ylim) min(ylim)+.8*diff(ylim)]);
    set(l, 'color', 'r')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.8*diff(ylim), 'Ipswich');
    set(t, 'fontsize', 8)
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.7*diff(ylim) min(ylim)+.7*diff(ylim)]);
    set(l, 'color', 'g')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.7*diff(ylim), 'Parker');
    set(t, 'fontsize', 8)
    datetick('x', 3, 'keeplimits', 'keepticks')
    box('on')
    
    subplot(4,1,2)
    plot(merrimackTime, merrimackFlow)
    ylabel('cfs')
    ylim([0 5e4])
    title('Streamflow:  Merrimack River')
    set(gca, 'xtick', xTick, 'ytick', [0:1e4:5e4])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    box('on')
    
    subplot(4,1,3)
    plot(julian2datenum(salTime2), salData2)
    hold on
    plot(julian2datenum(salTime3), salData3, 'r')
    plot(julian2datenum(salTime4), salData4, 'g')
    ylabel('psu')
    ylim([28 34])
    title('Salinity')
    set(gca, 'xtick', xTick, 'ytick', [28:2:34])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.3*diff(ylim) min(ylim)+.3*diff(ylim)]);
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.3*diff(ylim), '5 mbs');
    set(t, 'fontsize', 8);
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.2*diff(ylim) min(ylim)+.2*diff(ylim)]);
    set(l, 'color', 'r')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.2*diff(ylim), '10 mab');
    set(t, 'fontsize', 8)
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.1*diff(ylim) min(ylim)+.1*diff(ylim)]);
    set(l, 'color', 'g')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.1*diff(ylim), '1 mab');
    set(t, 'fontsize', 8)
    datetick('x', 3, 'keeplimits', 'keepticks')
    box('on')
    
    subplot(4,1,4)
    plot(julian2datenum(tempTime1), tempData1, 'm')
    hold on
    plot(julian2datenum(tempTime2), tempData2)
    plot(julian2datenum(tempTime3), tempData3, 'r')
    plot(julian2datenum(tempTime4), tempData4, 'g')
    ylabel('^oC')
    ylim([0 25])
    title('Temperature')
    set(gca, 'xtick', xTick, 'ytick', [0:5:25])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.9*diff(ylim) min(ylim)+.9*diff(ylim)]);
    set(l, 'color', 'm')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.9*diff(ylim), '0 m');
    set(t, 'fontsize', 8);
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.8*diff(ylim) min(ylim)+.8*diff(ylim)]);
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.8*diff(ylim), '5 mbs');
    set(t, 'fontsize', 8);
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.7*diff(ylim) min(ylim)+.7*diff(ylim)]);
    set(l, 'color', 'r')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.7*diff(ylim), '10 mab');
    set(t, 'fontsize', 8)
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.6*diff(ylim) min(ylim)+.6*diff(ylim)]);
    set(l, 'color', 'g')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.6*diff(ylim), '1 mab');
    set(t, 'fontsize', 8)
    datetick('x', 3, 'keeplimits', 'keepticks')
    box('on')
    
    suptitle({'River Streamflow, Salinity, and Temperature'; ...
            ['Site A, ' num2str(targetYear, '%.0f')]}) 
    
    print(fullfile(plotDirectory, ...
        ['a_river_' num2str(targetYear, '%.0f')]), '-dpdf');
    close
    
    %%%%%%%%%%%%%%%%%%%
    
    if targetYear < 1998
        continue
    end
    
    figure, orient landscape
       
    subplot(4,1,1)
    plot(charlesTime, charlesFlow)
    hold on
    plot(ipswichTime, ipswichFlow, 'r')
    plot(parkerTime, parkerFlow, 'g')
    ylabel('cfs')
    ylim([0 3000])
    title('Streamflow')
    set(gca, 'xtick', xTick, 'ytick', [0:500:3000])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])  
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.9*diff(ylim) min(ylim)+.9*diff(ylim)]);
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.9*diff(ylim), 'Charles');
    set(t, 'fontsize', 8);
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.8*diff(ylim) min(ylim)+.8*diff(ylim)]);
    set(l, 'color', 'r')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.8*diff(ylim), 'Ipswich');
    set(t, 'fontsize', 8)
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.7*diff(ylim) min(ylim)+.7*diff(ylim)]);
    set(l, 'color', 'g')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.7*diff(ylim), 'Parker');
    set(t, 'fontsize', 8)
    datetick('x', 3, 'keeplimits', 'keepticks')
    box('on')
    
    subplot(4,1,2)
    plot(merrimackTime, merrimackFlow)
    ylabel('cfs')
    ylim([0 5e4])
    title('Streamflow:  Merrimack River')
    set(gca, 'xtick', xTick, 'ytick', [0:1e4:5e4])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    box('on')
    
    subplot(4,1,3)
    hold on
    plot(julian2datenum(salTime3_b), salData3_b, 'r')
    plot(julian2datenum(salTime4_b), salData4_b, 'g')
    ylabel('psu')
    ylim([28 34])
    title('Salinity')
    set(gca, 'xtick', xTick, 'ytick', [28:2:34])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.2*diff(ylim) min(ylim)+.2*diff(ylim)]);
    set(l, 'color', 'r')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.2*diff(ylim), '10-16 mbs');
    set(t, 'fontsize', 8);
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.1*diff(ylim) min(ylim)+.1*diff(ylim)]);
    set(l, 'color', 'g')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.1*diff(ylim), '20-25 mbs');
    set(t, 'fontsize', 8)
    datetick('x', 3, 'keeplimits', 'keepticks')
    box('on')
    
    subplot(4,1,4)
    plot(julian2datenum(tempTime1), tempData1, 'm')
    hold on
    plot(julian2datenum(tempTime3_b), tempData3_b, 'r')
    plot(julian2datenum(tempTime4_b), tempData4_b, 'g')
    ylabel('^oC')
    ylim([0 25])
    title('Temperature')
    set(gca, 'xtick', xTick, 'ytick', [0:5:25])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
        l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.9*diff(ylim) min(ylim)+.9*diff(ylim)]);
    set(l, 'color', 'm')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.9*diff(ylim), '0 m');
    set(t, 'fontsize', 8);
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.8*diff(ylim) min(ylim)+.8*diff(ylim)]);
    set(l, 'color', 'r')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.8*diff(ylim), '10-16 mbs');
    set(t, 'fontsize', 8);
    l = line([datenum(targetYear, 1, 15, 00, 00, 00) datenum(targetYear, 1, 25, 00, 00, 00)], ...
        [min(ylim)+.7*diff(ylim) min(ylim)+.7*diff(ylim)]);
    set(l, 'color', 'g')
    t = text(datenum(targetYear, 1, 30, 00, 00, 00), min(ylim)+.7*diff(ylim), '20-25 mbs');
    set(t, 'fontsize', 8)
    datetick('x', 3, 'keeplimits', 'keepticks')
    box('on')
    
    suptitle({'River Streamflow, Salinity, and Temperature'; ...
            ['Site B, ' num2str(targetYear, '%.0f')]}) 
    
    print(fullfile(plotDirectory, ...
        ['b_river_' num2str(targetYear, '%.0f')]), '-dpdf');
    close
    
end