%meteorStatsNY.m
%
%Program to calculate various statistics (# of points, mean, st. deviation, 
%maximum, and minimum) for available meteorological data from stations near
%the New York Bight study area
%
%Soupy Alexander, 2/6/2

station = input('Select Ambrose Light (A) or Buoy 44025 (B): ');

if station == 'A';
    file_name1 = 'alsn6-a.nc';
elseif station == 'B';
    file_name1 = '44025A-a.nc';
    file_name2 = '44025B-a.nc';
else error('Improper station ID, must be capital letter in single quotes.');
end


%Pull out the data.
file = netcdf(file_name1,'nowrite');
t = file{'time'}(:);
t2 = file{'time2'}(:);
time = singleJD(t,t2);
air_temp = file{'AT_21'}(:);
ss_temp = file{'T_25'}(:);
windu = file{'WU_422'}(:);
windv = file{'WV_423'}(:);

%Concatenate the 44025 data.

if station == 'B';
    file2 = netcdf(file_name2,'nowrite');
    t_2 = file2{'time'}(:);
    t2_2 = file2{'time2'}(:);
    time2 = singleJD(t_2,t2_2);
    air_temp2 = file2{'AT_21'}(:);
    ss_temp2 = file2{'T_25'}(:);
    windu2 = file2{'WU_422'}(:);
    windv2 = file2{'WV_423'}(:);
    
    time = ([time' time2'])';
    air_temp = ([air_temp' air_temp2'])';
    ss_temp = ([ss_temp' ss_temp2'])';
    windu = ([windu' windu2'])';
    windv = ([windv' windv2'])';
end

air_temp = ridfill_nan(air_temp);
ss_temp = ridfill_nan(ss_temp);
windu = ridfill_nan(windu);
windv = ridfill_nan(windv);
windsp = sqrt(windu.^2 + windv.^2);

[taux,tauy,u10,v10] = wstress(windu,windv);

taux = ridfill_nan(taux);
tauy = ridfill_nan(tauy);

tauxy = sqrt(taux.^2 + tauy.^2);


%Now do month by month
greg_time = gregorian(time);

index_times = [12 1 2 3 4];

for index = 1:5;
    in_month = find(greg_time(:,2)==index_times(index));
    if length(in_month) > 0;
    n_windu(1,index) = length(find(isnan(windu(in_month))==0));
    n_windv(1,index) = length(find(isnan(windv(in_month))==0));
    n_airtemp(1,index) = length(find(isnan(air_temp(in_month))==0));
    n_sstemp(1,index) = length(find(isnan(ss_temp(in_month))==0));
    
    mean_windu(1,index) = nanmean(windu(in_month));
    mean_windv(1,index) = nanmean(windv(in_month));
    mean_airtemp(1,index) = nanmean(air_temp(in_month));
    mean_sstemp(1,index) = nanmean(ss_temp(in_month));
    mean_taux(1,index) = nanmean(taux(in_month));
    mean_tauy(1,index) = nanmean(tauy(in_month));
    mean_tauxy(1,index) = nanmean(tauxy(in_month));
    mean_windsp(1,index) = nanmean(windsp(in_month));
    
    
    std_windu(1,index) = nanstd(windu(in_month));
    std_windv(1,index) = nanstd(windv(in_month));
    std_airtemp(1,index) = nanstd(air_temp(in_month));
    std_sstemp(1,index) = nanstd(ss_temp(in_month));
    std_taux(1,index) = nanstd(taux(in_month));
    std_tauy(1,index) = nanstd(tauy(in_month));
    std_tauxy(1,index) = nanstd(tauxy(in_month));
    std_windsp(1,index) = nanstd(windsp(in_month));    

    max_windu(1,index) = nanmax(windu(in_month));
    max_windv(1,index) = nanmax(windv(in_month));
    max_airtemp(1,index) = nanmax(air_temp(in_month));
    max_sstemp(1,index) = nanmax(ss_temp(in_month));
    max_taux(1,index) = nanmax(taux(in_month));
    max_tauy(1,index) = nanmax(tauy(in_month));
    max_tauxy(1,index) = nanmax(tauxy(in_month));
    max_windsp(1,index) = nanmax(windsp(in_month));    

    min_windu(1,index) = nanmin(windu(in_month));
    min_windv(1,index) = nanmin(windv(in_month));
    min_airtemp(1,index) = nanmin(air_temp(in_month));
    min_sstemp(1,index) = nanmin(ss_temp(in_month));
    min_taux(1,index) = nanmin(taux(in_month));
    min_tauy(1,index) = nanmin(tauy(in_month));
    min_tauxy(1,index) = nanmin(tauxy(in_month));
    min_windsp(1,index) = nanmin(windsp(in_month));    
end
end
    
n_windu(1,6) = length(find(isnan(windu)==0))
n_windv(1,6) = length(find(isnan(windv)==0))
n_airtemp(1,6) = length(find(isnan(air_temp)==0))
n_sstemp(1,6) = length(find(isnan(ss_temp)==0))

mean_windu(1,6) = nanmean(windu)
mean_windv(1,6) = nanmean(windv)
mean_airtemp(1,6) = nanmean(air_temp)
mean_sstemp(1,6) = nanmean(ss_temp)
mean_taux(1,6) = nanmean(taux)
mean_tauy(1,6) = nanmean(tauy)
mean_tauxy(1,6) = nanmean(tauxy)
mean_windsp(1,6) = nanmean(windsp)

std_windu(1,6) = nanstd(windu)
std_windv(1,6) = nanstd(windv)
std_airtemp(1,6) = nanstd(air_temp)
std_sstemp(1,6) = nanstd(ss_temp)
std_taux(1,6) = nanstd(taux)
std_tauy(1,6) = nanstd(tauy)
std_tauxy(1,6) = nanstd(tauxy)
std_windsp(1,6) = nanstd(windsp)

max_windu(1,6) = nanmax(windu)
max_windv(1,6) = nanmax(windv)
max_airtemp(1,6) = nanmax(air_temp)
max_sstemp(1,6) = nanmax(ss_temp)
max_taux(1,6) = nanmax(taux)
max_tauy(1,6) = nanmax(tauy)
max_tauxy(1,6) = nanmax(tauxy)
max_windsp(1,6) = nanmax(windsp)

min_windu(1,6) = nanmin(windu)
min_windv(1,6) = nanmin(windv)
min_airtemp(1,6) = nanmin(air_temp)
min_sstemp(1,6) = nanmin(ss_temp)    
min_taux(1,6) = nanmin(taux)
min_tauy(1,6) = nanmin(tauy)
min_tauxy(1,6) = nanmin(tauxy)
min_windsp(1,6) = nanmin(windsp)