function mbiwe_weather
%
%   mbiwe_weather
%
%   Function to plot the wind speed, vector wind stress, and significant
%   wave height from the period of the Mass. Bay Internal Wave Experiment
%

%Soupy Alexander, 1/13/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.

theLimits = [datenum([1998 07 29 00 00 00]) datenum([1998 09 02 00 00 00])];
theTicks = min(theLimits):4:max(theLimits);

dataDirectory = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\DATAFILES';
plotDirectory = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\PLOTS';
windFile = 'lt_wind.nc';
waveFile = 'lt_wave.nc';

%Wave data
ncWave = netcdf(fullfile(dataDirectory, waveFile), 'nowrite');
waveTime = singlejd(ncWave{'time'}(:), ncWave{'time2'}(:));
waveData = ncWave{'wh_4061'}(:);
waveData(find(waveData >1e3 | waveData < -1e3)) = NaN;

%Wind data
ncWind = netcdf(fullfile(dataDirectory, windFile), 'nowrite');
windTime = singlejd(ncWind{'time'}(:), ncWind{'time2'}(:));
windU = ncWind{'WU_422'}(:);
windV = ncWind{'WV_423'}(:);
theBads = find(windU > 1e3 | windU < -1e3 | windV > 1e3 | windV < -1e3);
windU(theBads) = NaN;
windV(theBads) = NaN;
windSpeed = sqrt(windU.^2 + windV.^2);

%Set up the new time base & interpolate the data
windTime_new = min(windTime):1/24:max(windTime);
[windU_i] = smart_interp(windTime, windU, windTime_new, 5);
[windV_i] = smart_interp(windTime, windV, windTime_new, 5);
[windU_lp, windTime_lp] = plfilt(windU_i, windTime_new, 6);
[windV_lp] = plfilt(windV_i, windTime_new, 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);

figure, orient landscape
subplot(3, 1, 1)
plot(julian2datenum(waveTime), waveData)
xlim(theLimits)
ylim([0 4])
set(gca, 'xtick', theTicks)
datetick('x', 6, 'keepticks', 'keeplimits')
ylabel('m')
title('Significant Wave Height')

subplot(3,1,2)
plot(julian2datenum(windTime), windSpeed)
xlim(theLimits)
ylim([0 10])
set(gca, 'xtick', theTicks)
datetick('x', 6, 'keepticks', 'keeplimits')
ylabel('m/s')
title('Wind Speed')

subplot(3,1,3)
[theAngles, theLengths] = cart2pol(tauU, tauV);
theAngles = theAngles*180/pi;
sticksafe(julian2datenum(windTime_lp), ...
    zeros(size(theLengths)), theLengths, theAngles)
xlim(theLimits)
ylim([-2 2])
set(gca, 'xtick', theTicks)
datetick('x', 6, 'keepticks', 'keeplimits')
ylabel('dynes/cm^2')
title('Low-Passed Wind Stress')
box('on')
sticksafe

suptitle('Wind/Wave Data')
print('-dpdf', fullfile(plotDirectory, 'weather.pdf'))
close