% buoy_comparison_WHOTS_2021.m

% run this after running flux_eval_WHOTS_2021.m
% also, convert hourly WHOTS data from txt to matlab table format
% and save in WHOI data directory

% Sept 2021 BWB

%%  set up parameters, edit as desired and run
clear;
close('all');
fclose('all');

% date range: JD 236 to 244 for complete WHOTS cruise
firstday = 236 + datenum(2020,12,31);  % start datenum
lastday = 244 + datenum(2020,12,31);   % end datenum

plotit = true;  % for doing plots.

graphformat = '.png';  % select graphics format files
graphdevice = '-dpng'; % select graphic device
cruise = 'WHOTS_2021'; % Cruise name
ship = 'RV_Sette';    % Research vessel

% system specific path definitions
sysType = computer;
if strncmp(sysType,'MACI64',7)       % set Mac paths
    data_drive = '/Users/byronb/Raw_Data_Files.noindex';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs','flux');
    path_proc_data = fullfile(data_drive,cruise,ship,'flux','Processed');
    path_whoi_data = fullfile(data_drive,cruise,ship,'flux','WHOI');
    path_images = fullfile(data_drive,cruise,ship,'flux','Processed_Images','buoy_comparison');
    mkdir(path_images);
elseif strncmp(sysType,'PCWIN',7)    % set PC Windows paths
    data_drive = 'C:\Data';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs','flux');
    path_proc_data = fullfile(data_drive,cruise,ship,'flux','Processed');
    path_whoi_data = fullfile(data_drive,cruise,ship,'flux','WHOI');
    path_images = fullfile(data_drive,cruise,ship,'flux','Processed_Images','buoy_comparison');
    mkdir(path_images);
end

% matlab script path
restoredefaultpath
cd(path_prog);
addpath(genpath(path_prog));
savepath;

%% load 10-min and hourly data
% this assumes you have already run flux_eval_WHOTS_2021.m and you have
% converted the hourly mooring data files to matlab table format.

load(fullfile(path_proc_data,'WHOTS_2021_PSL_flux_1hr_all.mat'),'flxhr');
load(fullfile(path_proc_data,'WHOTS_2021_PSL_flux_10min_all.mat'),'flx');

load(fullfile(path_whoi_data,'WHOTS_16_sys1.mat'),'W_16_1');
load(fullfile(path_whoi_data,'WHOTS_16_sys2.mat'),'W_16_2');
load(fullfile(path_whoi_data,'WHOTS_17_sys1.mat'),'W_17_1');
load(fullfile(path_whoi_data,'WHOTS_17_sys2.mat'),'W_17_2');

%%  set up selection criteria

% set up timestamps in datenum format
flxhr.dnum = flxhr.jd + datenum(2020,12,31);
W_16_1.dnum = datenum(W_16_1.Year,W_16_1.mo,W_16_1.day,W_16_1.hr,0,0);
W_16_2.dnum = datenum(W_16_2.Year,W_16_2.mo,W_16_2.day,W_16_2.hr,0,0);
W_17_1.dnum = datenum(W_17_1.Year,W_17_1.mo,W_17_1.day,W_17_1.hr,0,0);
W_17_2.dnum = datenum(W_17_2.Year,W_17_2.mo,W_17_2.day,W_17_2.hr,0,0);

% selector for entire time at HOT site: Aug 27 0600 to Aug 31 1900
tt = find(flxhr.dnum>datenum(2021,8,27,6,0,0) & flxhr.dnum<datenum(2021,8,31,19,0,0));
flxhr = flxhr(tt,:);  % select out ship data for comparison period
N = height(flxhr);
dt = datetime(flxhr.dnum,'convertfrom','datenum'); % datetime formate timestamp for plots...

% select mooring data covering same time period as PSL measurements
% and copy over into the same timebase (some hours will be missing).
% getting rid of the data we don't need simplifies the comparison code to follow.

% first, create empty arrays for the mooring data, same length as PSL table
C = width(W_16_1);
w161 = NaN(N,C);
w161(:,C) = flxhr.dnum;
w162 = w161; w171 = w161; w172 = w161;

% now, find coincident hours for each system and copy data into arrays
for ii = 1:N
    jj = find(abs(W_16_1.dnum - w161(ii,C))<1/1440);
    if ~ isempty(jj)
        w161(ii,1:19) = W_16_1{jj,1:19};
    end
    jj = find(abs(W_16_2.dnum - w162(ii,C))<1/1440);
    if ~ isempty(jj)
        w162(ii,1:19) = W_16_2{jj,1:19};
    end
    jj = find(abs(W_17_1.dnum - w171(ii,C))<1/1440);
    if ~ isempty(jj)
        w171(ii,1:19) = W_17_1{jj,1:19};
    end
    jj = find(abs(W_17_2.dnum - w172(ii,C))<1/1440);
    if ~ isempty(jj)
        w172(ii,1:19) = W_17_2{jj,1:19};
    end
end

% convert arrays to tables
wvars = W_16_1.Properties.VariableNames;
w161 = array2table(w161,'VariableNames',wvars);
w162 = array2table(w162,'VariableNames',wvars);
w171 = array2table(w171,'VariableNames',wvars);
w172 = array2table(w172,'VariableNames',wvars);

% compute hourly distance from the ship to each mooring anchor locations
d16 = NaN(N,1); d17 = d16;
lat16 = 22 + 40.096/60; lon16 = -(157 + 56.788/60); % lat/lon WHOTS 16
lat17 = 22 + 45.045/60; lon17 = -(157 + 53.888/60); % lat/lon WHOTS 17
for ii = 1:N
    [~,d16(ii)] = legs([flxhr.lat(ii),lat16],[flxhr.lon(ii)-360,lon16]);
    [~,d17(ii)] = legs([flxhr.lat(ii),lat17],[flxhr.lon(ii)-360,lon17]);    
end

% selector for period closest to WHOTS17: Aug 30 1200 to Aug 31 1600
tt17 = find(flxhr.dnum>datenum(2021,8,30,12,0,0) & flxhr.dnum<datenum(2021,8,31,16,0,0));

% selector for period closest to WHOTS16: Aug 27 0800 to Aug 28 1600
tt16 = find(flxhr.dnum>datenum(2021,8,27,8,0,0) & flxhr.dnum<datenum(2021,8,28,16,0,0));

% selector for both intercomparison periods, only usable with SCS and WHOTS-17 data
% since WHOTS-16 was recovered on 8/28 1800
all = union(tt16,tt17);

%% met comparison plots

% Distance to Moorings, timeseries
figure('position',[300,300,800,600],'name','Distance_to_moorings');
plot(dt,d16,'.',dt,d17,'.','markersize',10); grid on; grid minor;
ylabel('Distance, N. Miles','fontsize',12);
title('WHOTS 2021, Distance, Ship to Mooring Anchor Location','fontsize',16);
legend('WHOTS-16','WHOTS-17','location','northeast','fontsize',12);
if plotit
    ppath = fullfile(path_images,['01_dist_to_mooring',graphformat]);
    print(graphdevice,ppath);
end

% SL Pressure
padj = 3*0.125; % adjustment of mooring P from 3m to sea level
figure('position',[300,300,800,700],'name','Pressure_all');
plot([1014,1019],[1014,1019],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.Pmb(tt16),w161.BPR(tt16)+padj,'o',flxhr.Pmb(tt16),w162.BPR(tt16)+padj,'o',...
          flxhr.Pmb(all),w171.BPR(all)+padj,'o',flxhr.Pmb(all),w172.BPR(all)+padj,'o',...
          flxhr.Pmb(all),flxhr.Pmb_scs(all),'o','linewidth',1.5); hold off;
xlabel('Pressure, mb (PSL PTB)','fontsize',14);
ylabel('Pressure, mb (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Sea Level Pressure Comparison','fontsize',16);
subtitle('For both intercomparison periods: Aug 27 - Aug 31, 2021','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2','SCS', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['02_pressure_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','Pressure_nearby');
plot([1014,1019],[1014,1019],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.Pmb(tt16),w161.BPR(tt16)+padj,'o',flxhr.Pmb(tt16),w162.BPR(tt16)+padj,'o',...
          flxhr.Pmb(tt17),w171.BPR(tt17)+padj,'o',flxhr.Pmb(tt17),w172.BPR(tt17)+padj,'o',...
          'linewidth',1.5); hold off;
xlabel('Pressure, mb (PSL PTB)','fontsize',14);
ylabel('Pressure, mb (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Sea Level Pressure Comparison','fontsize',16);
subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['03_pressure_nearby',graphformat]);
    print(graphdevice,ppath);
end

% Air Temperature (3m)
figure('position',[300,300,800,700],'name','Tair_all');
plot([24,26],[24,26],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.Trf(tt16),w161.ATMP(tt16),'o',flxhr.Trf(tt16),w162.ATMP(tt16),'o',...
          flxhr.Trf(all),w171.ATMP(all),'o',flxhr.Trf(all),w172.ATMP(all),'o',...
          flxhr.Trf_scs(all),flxhr.ta_scs(all),'o','linewidth',1.5); hold off;
axis([24,26,24,26]);
xlabel('Air Temp, °C (COARE Ref)','fontsize',14);
ylabel('Air Temp, °C (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Air Temperature Comparison','fontsize',16);
subtitle('PSL Ta adjusted to 3m for moorings and 21m for SCS','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2','SCS','location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['04_tair_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','Tair_nearby');
plot([24,26],[24,26],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.Trf(tt16),w161.ATMP(tt16),'o',flxhr.Trf(tt16),w162.ATMP(tt16),'o',...
          flxhr.Trf(tt17),w171.ATMP(tt17),'o',flxhr.Trf(tt17),w172.ATMP(tt17),'o',...
          'linewidth',1.5); hold off;
axis([24,26,24,26]);
xlabel('Air Temp, °C (COARE 3m Ref)','fontsize',14);
ylabel('Air Temp, °C (moorings)','fontsize',14);
title('WHOTS 2021, Hourly 3m Air Temperature Comparison','fontsize',16);
subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['05_tair_nearby',graphformat]);
    print(graphdevice,ppath);
end

% Relative Humidity (3m)
figure('position',[300,300,800,700],'name','RH_all');
plot([66,84],[66,84],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.RHrf(tt16),w161.RH(tt16),'o',flxhr.RHrf(tt16),w162.RH(tt16),'o',...
          flxhr.RHrf(all),w171.RH(all),'o',flxhr.RHrf(all),w172.RH(all),'o',...
          flxhr.RHrf_scs(all),flxhr.rh_scs(all),'o','linewidth',1.5); hold off;
xlabel('RH %, (COARE Ref)','fontsize',14);
ylabel('RH %, (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Relative Humidity Comparison','fontsize',16);
subtitle('PSL RH adjusted to 3m for moorings and 21m for SCS','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2','SCS','location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['06_rh_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','RH_nearby');
plot([68,86],[68,86],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.RHrf(tt16),w161.RH(tt16),'o',flxhr.RHrf(tt16),w162.RH(tt16),'o',...
          flxhr.RHrf(tt17),w171.RH(tt17),'o',flxhr.RHrf(tt17),w172.RH(tt17),'o',...
          'linewidth',1.5); hold off;
xlabel('RH %, (COARE 3m Ref)','fontsize',14);
ylabel('RH %, (moorings)','fontsize',14);
title('WHOTS 2021, Hourly 3m Relative Humidity Comparison','fontsize',16);
subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['07_rh_nearby',graphformat]);
    print(graphdevice,ppath);
end

% Shortwave Flux
figure('position',[300,300,800,700],'name','RS_all');
plot([-200,1200],[-200,1200],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.rs(tt16),w161.SRAD(tt16),'o',flxhr.rs(tt16),w162.SRAD(tt16),'o',...
          flxhr.rs(all),w171.SRAD(all),'o',flxhr.rs(all),w172.SRAD(all),'o',...
          'linewidth',1.5); hold off;
xlabel('Shortwave Radiation, W/m^2 (PSL K&Z)','fontsize',14);
ylabel('Shortwave Radiation, W/m^2 (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Shortwave Radiation Comparison','fontsize',16);
subtitle('For both intercomparison periods: Aug 27 - Aug 31, 2021','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['08_rs_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','RS_nearby');
plot([-200,1200],[-200,1200],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.rs(tt16),w161.SRAD(tt16),'o',flxhr.rs(tt16),w162.SRAD(tt16),'o',...
          flxhr.rs(tt17),w171.SRAD(tt17),'o',flxhr.rs(tt17),w172.SRAD(tt17),'o',...
          'linewidth',1.5); hold off;
xlabel('Shortwave Radiation, W/m^2 (PSL K&Z)','fontsize',14);
ylabel('Shortwave Radiation, W/m^2 (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Shortwave Radiation Comparison','fontsize',16);
subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['09_rs_nearby',graphformat]);
    print(graphdevice,ppath);
end

% Longwave Flux
figure('position',[300,300,800,700],'name','RL_all');
plot([360,440],[360,440],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.rl(tt16),w161.LRAD(tt16),'o',flxhr.rl(tt16),w162.LRAD(tt16),'o',...
          flxhr.rl(all),w171.LRAD(all),'o',flxhr.rl(all),w172.LRAD(all),'o',...
          'linewidth',1.5); hold off;
xlabel('Longwave Radiation, W/m^2 (PSL PIR)','fontsize',14);
ylabel('Longwave Radiation, W/m^2 (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Longwave Radiation Comparison','fontsize',16);
subtitle('For both intercomparison periods: Aug 27 - Aug 31, 2021','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['10_rl_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','RL_nearby');
plot([360,440],[360,440],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.rl(tt16),w161.LRAD(tt16),'o',flxhr.rl(tt16),w162.LRAD(tt16),'o',...
          flxhr.rl(tt17),w171.LRAD(tt17),'o',flxhr.rl(tt17),w172.LRAD(tt17),'o',...
          'linewidth',1.5); hold off;
xlabel('Longwave Radiation, W/m^2 (PSL PIR)','fontsize',14);
ylabel('Longwave Radiation, W/m^2 (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Longwave Radiation Comparison','fontsize',16);
subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['11_rl_nearby',graphformat]);
    print(graphdevice,ppath);
end

% SST
figure('position',[300,300,800,700],'name','SST_all');
plot([26,27],[26,27],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.tsea(tt16),w161.STMP(tt16),'o',flxhr.tsea(tt16),w162.STMP(tt16),'o',...
          flxhr.tsea(all),w171.STMP(all),'o',flxhr.tsea(all),w172.STMP(all),'o',...
          'linewidth',1.5); hold off;
xlabel('Sea Temp, °C, (PSL Seasnake)','fontsize',14);
ylabel('Sea Temp, °C (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Sea Surface Comparison','fontsize',16);
subtitle('For both intercomparison periods: Aug 27 - Aug 31, 2021','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['12_sst_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','SST_nearby');
plot([26,27],[26,27],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.tsea(tt16),w161.STMP(tt16),'o',flxhr.tsea(tt16),w162.STMP(tt16),'o',...
          flxhr.tsea(tt17),w171.STMP(tt17),'o',flxhr.tsea(tt17),w172.STMP(tt17),'o',...
          'linewidth',1.5); hold off;
xlabel('Sea Temp, °C, (PSL Seasnake)','fontsize',14);
ylabel('Sea Temp, °C (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Sea Surface Temperature Comparison','fontsize',16);
subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['13_sst_nearby',graphformat]);
    print(graphdevice,ppath);
end

% Wind Speed (3m)
figure('position',[300,300,800,700],'name', 'WSPD_all');
plot([0,10],[0,10],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.Urf(tt16),w161.WSPD(tt16),'o',flxhr.Urf(tt16),w162.WSPD(tt16),'o',...
          flxhr.Urf(all),w171.WSPD(all),'o',flxhr.Urf(all),w172.WSPD(all),'o',...
          flxhr.Urf_scs(all),flxhr.wspd_scs(all),'o','linewidth',1.5); hold off;
xlabel('Wind Speed, m/s, (COARE Ref)','fontsize',14);
ylabel('Wind Speed, m/s, (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Wind Speed Comparison','fontsize',16);
subtitle('PSL Wspd adjusted to 3m for moorings and 21m for SCS','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2','SCS','location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['14_wspd_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','WSPD_nearby');
plot([0,10],[0,10],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.Urf(tt16),w161.WSPD(tt16),'o',flxhr.Urf(tt16),w162.WSPD(tt16),'o',...
          flxhr.Urf(tt17),w171.WSPD(tt17),'o',flxhr.Urf(tt17),w172.WSPD(tt17),'o',...
          'linewidth',1.5); hold off;
xlabel('Wind Speed, m/s, (COARE 3m Ref)','fontsize',14);
ylabel('Wind Speed, m/s, (moorings)','fontsize',14);
title('WHOTS 2021, Hourly 3m Wind Speed Comparison','fontsize',16);
subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['15_wspd_nearby',graphformat]);
    print(graphdevice,ppath);
end

% Wind Direction
wdir = flxhr.wdir + 180;   % to oceanographic convention 'direction to'
wdir(wdir>360) = wdir(wdir>360) - 360;
wdir(wdir<0) = wdir(wdir<0) + 360;
wdir_scs = flxhr.wdir_scs + 180;   % to oceanographic convention 'direction to'
wdir_scs(wdir_scs>360) = wdir_scs(wdir_scs>360) - 360;
wdir_scs(wdir_scs<0) = wdir_scs(wdir_scs<0) + 360;
figure('position',[300,300,800,700],'name','WDIR_all');
plot([0,360],[0,360],'k-'); grid on; grid minor; hold on;
h1 = plot(wdir(tt16),w161.WDIR(tt16),'o',wdir(tt16),w162.WDIR(tt16),'o',...
          wdir(all),w171.WDIR(all),'o',wdir(all),w172.WDIR(all),'o',...
          wdir(all),wdir_scs(all),'o','linewidth',1.5); hold off;
axis([200,340,200,340]);
xlabel('Wind Direction ''to'', deg (PSL WXT)','fontsize',14);
ylabel('Wind Direction ''to'', deg (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Wind Direction (Ocn Convention)','fontsize',16);
subtitle('For both intercomparison periods: Aug 27 - Aug 31, 2021','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2','SCS','location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['16_wdir_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','WDIR_nearby');
plot([0,360],[0,360],'k-'); grid on; grid minor; hold on;
h1 = plot(wdir(tt16),w161.WDIR(tt16),'o',wdir(tt16),w162.WDIR(tt16),'o',...
          wdir(tt17),w171.WDIR(tt17),'o',wdir(tt17),w172.WDIR(tt17),'o',...
          'linewidth',1.5); hold off;
axis([200,340,200,340]);
xlabel('Wind Direction ''to'', deg (PSL WXT)','fontsize',14);
ylabel('Wind Direction ''to'', deg (moorings)','fontsize',14);
title('WHOTS 2021, Hourly Wind Direction (Ocn Convention)','fontsize',16);
subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'W16-1','W16-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['17_wdir_nearby',graphformat]);
    print(graphdevice,ppath);
end

%% prepare vectors of stats for the ship and buoy comparison:
% comparison is PSL value – ship or mooring value

% for entire period at HOT site:
psl_scs_all_means  = [nanmean1(flxhr.Trf_scs(all) - flxhr.ta_scs(all)),...
                      nanmean1(flxhr.RHrf_scs(all) - flxhr.rh_scs(all)),...
                      NaN,... % no sst
                      nanmean1(flxhr.Urf_scs(all) - flxhr.wspd_scs(all)),...
                      nanmean1(flxhr.wdir(all) - flxhr.wdir_scs(all)),...
                      NaN,... % no rs
                      NaN,... % no rl
                      nanmean1(flxhr.Pmb(all) - flxhr.Pmb_scs(all))];

psl_scs_all_meds  =  [nanmedian1(flxhr.Trf_scs(all) - flxhr.ta_scs(all)),...
                      nanmedian1(flxhr.RHrf_scs(all) - flxhr.rh_scs(all)),...
                      NaN,...
                      nanmedian1(flxhr.Urf_scs(all) - flxhr.wspd_scs(all)),...
                      nanmedian1(flxhr.wdir(all) - flxhr.wdir_scs(all)),...
                      NaN,...
                      NaN,...
                      nanmedian1(flxhr.Pmb(all) - flxhr.Pmb_scs(all))];

psl_scs_all_std  =   [nanstd1(flxhr.Trf_scs(all) - flxhr.ta_scs(all)),...
                      nanstd1(flxhr.RHrf_scs(all) - flxhr.rh_scs(all)),...
                      NaN,...
                      nanstd1(flxhr.Urf_scs(all) - flxhr.wspd_scs(all)),...
                      nanstd1(flxhr.wdir(all) - flxhr.wdir_scs(all)),...
                      NaN,...
                      NaN,...
                      nanstd1(flxhr.Pmb(all) - flxhr.Pmb_scs(all))];

psl_w171_all_means = [nanmean1(flxhr.Trf(all) - w171.ATMP(all)),...
                      nanmean1(flxhr.RHrf(all) - w171.RH(all))...
                      nanmean1(flxhr.tsea(all) - w171.STMP(all))...
                      nanmean1(flxhr.Urf(all) - w171.WSPD(all))...
                      nanmean1(wdir(all) - w171.WDIR(all))...
                      nanmean1(flxhr.rs(all) - w171.SRAD(all))...
                      nanmean1(flxhr.rl(all) - w171.LRAD(all))...
                      nanmean1(flxhr.Pmb(all) - w171.BPR(all))];

psl_w171_all_meds =  [nanmedian1(flxhr.Trf(all) - w171.ATMP(all)),...
                      nanmedian1(flxhr.RHrf(all) - w171.RH(all))...
                      nanmedian1(flxhr.tsea(all) - w171.STMP(all))...
                      nanmedian1(flxhr.Urf(all) - w171.WSPD(all))...
                      nanmedian1(wdir(all) - w171.WDIR(all))...
                      nanmedian1(flxhr.rs(all) - w171.SRAD(all))...
                      nanmedian1(flxhr.rl(all) - w171.LRAD(all))...
                      nanmedian1(flxhr.Pmb(all) - w171.BPR(all))];

psl_w171_all_std =   [nanstd1(flxhr.Trf(all) - w171.ATMP(all)),...
                      nanstd1(flxhr.RHrf(all) - w171.RH(all))...
                      nanstd1(flxhr.tsea(all) - w171.STMP(all))...
                      nanstd1(flxhr.Urf(all) - w171.WSPD(all))...
                      nanstd1(wdir(all) - w171.WDIR(all))...
                      nanstd1(flxhr.rs(all) - w171.SRAD(all))...
                      nanstd1(flxhr.rl(all) - w171.LRAD(all))...
                      nanstd1(flxhr.Pmb(all) - w171.BPR(all))];

psl_w172_all_means = [nanmean1(flxhr.Trf(all) - w172.ATMP(all)),...
                      nanmean1(flxhr.RHrf(all) - w172.RH(all))...
                      nanmean1(flxhr.tsea(all) - w172.STMP(all))...
                      nanmean1(flxhr.Urf(all) - w172.WSPD(all))...
                      nanmean1(wdir(all) - w172.WDIR(all))...
                      nanmean1(flxhr.rs(all) - w172.SRAD(all))...
                      nanmean1(flxhr.rl(all) - w172.LRAD(all))...
                      nanmean1(flxhr.Pmb(all) - w172.BPR(all))];

psl_w172_all_meds =  [nanmedian1(flxhr.Trf(all) - w172.ATMP(all)),...
                      nanmedian1(flxhr.RHrf(all) - w172.RH(all))...
                      nanmedian1(flxhr.tsea(all) - w172.STMP(all))...
                      nanmedian1(flxhr.Urf(all) - w172.WSPD(all))...
                      nanmedian1(wdir(all) - w172.WDIR(all))...
                      nanmedian1(flxhr.rs(all) - w172.SRAD(all))...
                      nanmedian1(flxhr.rl(all) - w172.LRAD(all))...
                      nanmedian1(flxhr.Pmb(all) - w172.BPR(all))];

psl_w172_all_std =   [nanstd1(flxhr.Trf(all) - w172.ATMP(all)),...
                      nanstd1(flxhr.RHrf(all) - w172.RH(all))...
                      nanstd1(flxhr.tsea(all) - w172.STMP(all))...
                      nanstd1(flxhr.Urf(all) - w172.WSPD(all))...
                      nanstd1(wdir(all) - w172.WDIR(all))...
                      nanstd1(flxhr.rs(all) - w172.SRAD(all))...
                      nanstd1(flxhr.rl(all) - w172.LRAD(all))...
                      nanstd1(flxhr.Pmb(all) - w172.BPR(all))];

% Now compute same for comparison intervals near each mooring
psl_w161_nearby_means = [nanmean1(flxhr.Trf(tt16) - w161.ATMP(tt16)),...
                      nanmean1(flxhr.RHrf(tt16) - w161.RH(tt16))...
                      nanmean1(flxhr.tsea(tt16) - w161.STMP(tt16))...
                      nanmean1(flxhr.Urf(tt16) - w161.WSPD(tt16))...
                      nanmean1(wdir(tt16) - w161.WDIR(tt16))...
                      nanmean1(flxhr.rs(tt16) - w161.SRAD(tt16))...
                      nanmean1(flxhr.rl(tt16) - w161.LRAD(tt16))...
                      nanmean1(flxhr.Pmb(tt16) - w161.BPR(tt16))];

psl_w161_nearby_meds =  [nanmedian1(flxhr.Trf(tt16) - w161.ATMP(tt16)),...
                      nanmedian1(flxhr.RHrf(tt16) - w161.RH(tt16))...
                      nanmedian1(flxhr.tsea(tt16) - w161.STMP(tt16))...
                      nanmedian1(flxhr.Urf(tt16) - w161.WSPD(tt16))...
                      nanmedian1(wdir(tt16) - w161.WDIR(tt16))...
                      nanmedian1(flxhr.rs(tt16) - w161.SRAD(tt16))...
                      nanmedian1(flxhr.rl(tt16) - w161.LRAD(tt16))...
                      nanmedian1(flxhr.Pmb(tt16) - w161.BPR(tt16))];

psl_w161_nearby_std =   [nanstd1(flxhr.Trf(tt16) - w161.ATMP(tt16)),...
                      nanstd1(flxhr.RHrf(tt16) - w161.RH(tt16))...
                      nanstd1(flxhr.tsea(tt16) - w161.STMP(tt16))...
                      nanstd1(flxhr.Urf(tt16) - w161.WSPD(tt16))...
                      nanstd1(wdir(tt16) - w161.WDIR(tt16))...
                      nanstd1(flxhr.rs(tt16) - w161.SRAD(tt16))...
                      nanstd1(flxhr.rl(tt16) - w161.LRAD(tt16))...
                      nanstd1(flxhr.Pmb(tt16) - w161.BPR(tt16))];

psl_w162_nearby_means = [nanmean1(flxhr.Trf(tt16) - w162.ATMP(tt16)),...
                      nanmean1(flxhr.RHrf(tt16) - w162.RH(tt16))...
                      nanmean1(flxhr.tsea(tt16) - w162.STMP(tt16))...
                      nanmean1(flxhr.Urf(tt16) - w162.WSPD(tt16))...
                      nanmean1(wdir(tt16) - w162.WDIR(tt16))...
                      nanmean1(flxhr.rs(tt16) - w162.SRAD(tt16))...
                      nanmean1(flxhr.rl(tt16) - w162.LRAD(tt16))...
                      nanmean1(flxhr.Pmb(tt16) - w162.BPR(tt16))];

psl_w162_nearby_meds =  [nanmedian1(flxhr.Trf(tt16) - w162.ATMP(tt16)),...
                      nanmedian1(flxhr.RHrf(tt16) - w162.RH(tt16))...
                      nanmedian1(flxhr.tsea(tt16) - w162.STMP(tt16))...
                      nanmedian1(flxhr.Urf(tt16) - w162.WSPD(tt16))...
                      nanmedian1(wdir(tt16) - w162.WDIR(tt16))...
                      nanmedian1(flxhr.rs(tt16) - w162.SRAD(tt16))...
                      nanmedian1(flxhr.rl(tt16) - w162.LRAD(tt16))...
                      nanmedian1(flxhr.Pmb(tt16) - w162.BPR(tt16))];

psl_w162_nearby_std =   [nanstd1(flxhr.Trf(tt16) - w162.ATMP(tt16)),...
                      nanstd1(flxhr.RHrf(tt16) - w162.RH(tt16))...
                      nanstd1(flxhr.tsea(tt16) - w162.STMP(tt16))...
                      nanstd1(flxhr.Urf(tt16) - w162.WSPD(tt16))...
                      nanstd1(wdir(tt16) - w162.WDIR(tt16))...
                      nanstd1(flxhr.rs(tt16) - w162.SRAD(tt16))...
                      nanstd1(flxhr.rl(tt16) - w162.LRAD(tt16))...
                      nanstd1(flxhr.Pmb(tt16) - w162.BPR(tt16))];

psl_w171_nearby_means = [nanmean1(flxhr.Trf(tt17) - w171.ATMP(tt17)),...
                      nanmean1(flxhr.RHrf(tt17) - w171.RH(tt17))...
                      nanmean1(flxhr.tsea(tt17) - w171.STMP(tt17))...
                      nanmean1(flxhr.Urf(tt17) - w171.WSPD(tt17))...
                      nanmean1(wdir(tt17) - w171.WDIR(tt17))...
                      nanmean1(flxhr.rs(tt17) - w171.SRAD(tt17))...
                      nanmean1(flxhr.rl(tt17) - w171.LRAD(tt17))...
                      nanmean1(flxhr.Pmb(tt17) - w171.BPR(tt17))];

psl_w171_nearby_meds =  [nanmedian1(flxhr.Trf(tt17) - w171.ATMP(tt17)),...
                      nanmedian1(flxhr.RHrf(tt17) - w171.RH(tt17))...
                      nanmedian1(flxhr.tsea(tt17) - w171.STMP(tt17))...
                      nanmedian1(flxhr.Urf(tt17) - w171.WSPD(tt17))...
                      nanmedian1(wdir(tt17) - w171.WDIR(tt17))...
                      nanmedian1(flxhr.rs(tt17) - w171.SRAD(tt17))...
                      nanmedian1(flxhr.rl(tt17) - w171.LRAD(tt17))...
                      nanmedian1(flxhr.Pmb(tt17) - w171.BPR(tt17))];

psl_w171_nearby_std =   [nanstd1(flxhr.Trf(tt17) - w171.ATMP(tt17)),...
                      nanstd1(flxhr.RHrf(tt17) - w171.RH(tt17))...
                      nanstd1(flxhr.tsea(tt17) - w171.STMP(tt17))...
                      nanstd1(flxhr.Urf(tt17) - w171.WSPD(tt17))...
                      nanstd1(wdir(tt17) - w171.WDIR(tt17))...
                      nanstd1(flxhr.rs(tt17) - w171.SRAD(tt17))...
                      nanstd1(flxhr.rl(tt17) - w171.LRAD(tt17))...
                      nanstd1(flxhr.Pmb(tt17) - w171.BPR(tt17))];

psl_w172_nearby_means = [nanmean1(flxhr.Trf(tt17) - w172.ATMP(tt17)),...
                      nanmean1(flxhr.RHrf(tt17) - w172.RH(tt17))...
                      nanmean1(flxhr.tsea(tt17) - w172.STMP(tt17))...
                      nanmean1(flxhr.Urf(tt17) - w172.WSPD(tt17))...
                      nanmean1(wdir(tt17) - w172.WDIR(tt17))...
                      nanmean1(flxhr.rs(tt17) - w172.SRAD(tt17))...
                      nanmean1(flxhr.rl(tt17) - w172.LRAD(tt17))...
                      nanmean1(flxhr.Pmb(tt17) - w172.BPR(tt17))];

psl_w172_nearby_meds =  [nanmedian1(flxhr.Trf(tt17) - w172.ATMP(tt17)),...
                      nanmedian1(flxhr.RHrf(tt17) - w172.RH(tt17))...
                      nanmedian1(flxhr.tsea(tt17) - w172.STMP(tt17))...
                      nanmedian1(flxhr.Urf(tt17) - w172.WSPD(tt17))...
                      nanmedian1(wdir(tt17) - w172.WDIR(tt17))...
                      nanmedian1(flxhr.rs(tt17) - w172.SRAD(tt17))...
                      nanmedian1(flxhr.rl(tt17) - w172.LRAD(tt17))...
                      nanmedian1(flxhr.Pmb(tt17) - w172.BPR(tt17))];

psl_w172_nearby_std =   [nanstd1(flxhr.Trf(tt17) - w172.ATMP(tt17)),...
                      nanstd1(flxhr.RHrf(tt17) - w172.RH(tt17))...
                      nanstd1(flxhr.tsea(tt17) - w172.STMP(tt17))...
                      nanstd1(flxhr.Urf(tt17) - w172.WSPD(tt17))...
                      nanstd1(wdir(tt17) - w172.WDIR(tt17))...
                      nanstd1(flxhr.rs(tt17) - w172.SRAD(tt17))...
                      nanstd1(flxhr.rl(tt17) - w172.LRAD(tt17))...
                      nanstd1(flxhr.Pmb(tt17) - w172.BPR(tt17))];

% save stats vector array as a txt file for import into spreadsheet
out = [psl_scs_all_means;
       psl_scs_all_meds;
       psl_scs_all_std;
       psl_w171_all_means;
       psl_w171_all_meds;
       psl_w171_all_std;
       psl_w172_all_means;
       psl_w172_all_meds;
       psl_w172_all_std;
       psl_w161_nearby_means;
       psl_w161_nearby_meds;
       psl_w161_nearby_std;
       psl_w162_nearby_means;
       psl_w162_nearby_meds;
       psl_w162_nearby_std;
       psl_w171_nearby_means;
       psl_w171_nearby_meds;
       psl_w171_nearby_std;
       psl_w172_nearby_means;
       psl_w172_nearby_meds;
       psl_w172_nearby_std];

sfl = fullfile(data_drive,cruise,ship,'Scientific_analysis','report','comp_table.txt');
save(sfl,'out','-ascii','-double','-tabs');

