% a19_compare_oot_ship.m
% compare ship flux results computed with OOT and PSL EC methods

%% initialize
clear
close all
warning('off','MATLAB:MKDIR:DirectoryExists');

% paths
root_data_path = '/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019';
path_prog = fullfile(root_data_path,'_matlab');
path_oot = fullfile(root_data_path,'_processed','ship_flux_v3','flux_merged_OOT_v1');
path_blo = fullfile(root_data_path,'_processed','ship_flux_v3','flux_merged');
path_plots = fullfile(path_oot,'comprison plots'); mkdir(path_plots);

% set matlab script path
cd(path_prog);
addpath(genpath(path_prog));
cd(fullfile(path_prog,'co2_ch4_o3_flux'));

%% load data and average to 1 hr
load(fullfile(path_oot,'flx_ship_oot.mat'));               % loads oot 30-min table
load(fullfile(path_blo,'flux_co2_ch4_o3_dms_hr_all.mat')); % loads BLO flx_hr

% average OOT result to hourly timescale
bins_hr = [flx_hr.datenum;flx_hr.datenum(end)+1/24];

% average OOT - first unwrap wdir, then rewrap after averaging
temp = table2array(oot);
temp(:,6) = unwrap(temp(:,6)*pi/180);
ship_oot = interval_avg(temp(:,1),temp(:,2:end),bins_hr);
ship_oot(:,6) = mod(ship_oot(:,6)*180/pi,360);
ii = ship_oot(:,6)<0; ship_oot(ii,6) = ship_oot(ii,6) + 360;
oot_vars = oot.Properties.VariableNames;
ship_oot = array2table(ship_oot,'VariableNames',oot_vars);

% load hourly metcity file for surface temperature data...
load(fullfile(root_data_path,'_processed','metcity_tower_level2_v4','metcity_1hr_flux_all.mat'));

%% plots

% 1 compare mean wind
figure('Position',[300,300,800,800],'Name','wspd');
plot([0,20],[0,20],'k-','LineWidth',1); hold on; grid on; grid minor;
plot(flx_hr.wspd_hr,ship_oot.u,'.','MarkerSize',10); hold off;
xlabel('wspd, BLO, m/s','FontSize',14);
ylabel('wspd OOT, m/s','FontSize',14);
pname = fullfile(path_plots,'Fig01_wspd_comparison.png');
print(pname,'-dpng');

% 2 wind binned
[wspd_blo_oot, wspd_oot_mean, wspd_oot_med, SD_wspd_oot, N_wspd_oot] = binave3(flx_hr.wspd_hr,ship_oot.u,2,0,20);
figure('Position',[300,300,800,800],'Name','wspd_bin');
plot([0,20],[0,20],'k-','LineWidth',1); hold on; grid on; grid minor;
errorbar(wspd_blo_oot,wspd_oot_mean,SD_wspd_oot,'square','MarkerSize',12,'LineWidth',2); hold off;
xlabel('wspd, BLO, m/s','FontSize',14);
ylabel('wspd, OOT, m/s','FontSize',14);
pname = fullfile(path_plots,'Fig02_wspd_comparison_bin.png');
print(pname,'-dpng');

% 3 compare temperature
figure('Position',[300,300,800,800],'Name','tson');
plot([-40,5],[-40,5],'k-','LineWidth',1); hold on; grid on; grid minor;
h12 = plot(flx_hr.Ta_Licor_hr,ship_oot.Tair,'.','MarkerSize',10); hold off;
xlabel('Tlicor, BLO, °C','FontSize',14);
ylabel('Tsonic, OOT, °C','FontSize',14);
pname = fullfile(path_plots,'Fig03_Tsonic_comparison.png');
print(pname,'-dpng');

% 4 compare momentum flux
% selector for data where both blo and oot return non-NaN Tau result
jj = isfinite(ship_oot.tau_OO) & isfinite(flx_hr.wu_csp_hr);
tau_blo = -flx_hr.rhoair_hr.*flx_hr.wu_csp_hr;
figure('Position',[300,300,800,800],'Name','tau');
plot([-0.1,0.5],[-0.1,0.5],'k-','LineWidth',1); hold on;
h1 = plot(tau_blo(jj),-ship_oot.tau(jj),'.',tau_blo(jj),-ship_oot.tau_OO(jj),'.','MarkerSize',10);
grid on; grid minor; xlabel('Tau, BLO, N/m2','FontSize',14); hold off;
ylabel('Tau, OOT, N/m2','FontSize',14);
legend(h1,'OOT-cov','OOT','Location','northwest','fontsize',14);
pname = fullfile(path_plots,'Fig04_tau_comparison.png');
print(pname,'-dpng');

% 5 tau binned
[tau_blo_oot, tau_oot_mean, tau_oot_med, SD_tau_oot, N_tau_oot] = binave3(tau_blo(jj),-ship_oot.tau(jj),0.04,-0.02,0.48);
[tau_blo_ootOO, tau_oot_meanOO, tau_oot_medOO, SD_tau_ootOO, N_tau_ootOO] = binave3(tau_blo(jj),-ship_oot.tau_OO(jj),0.04,-0.02,0.48);
figure('Position',[300,300,800,800],'Name','tau_bin');
plot([-0.1,0.5],[-0.1,0.5],'k-','LineWidth',1); hold on;
h3 = errorbar(tau_blo_oot,tau_oot_mean,SD_tau_oot,'square','MarkerSize',12,'LineWidth',2);
h4 = errorbar(tau_blo_ootOO,tau_oot_meanOO,SD_tau_ootOO,'diamond','MarkerSize',12,'LineWidth',2);
grid on; grid minor; xlabel('Tau, BLO, N/m2','FontSize',14); hold off;
ylabel('Tau, OOT, N/m2','FontSize',14);
legend([h3,h4],'OOT-cov','OOT','Location','northwest','fontsize',14);
pname = fullfile(path_plots,'Fig05_tau_comparison_bin.png');
print(pname,'-dpng');

% 6 compare sensible heat flux
% selector for data where both blo and oot return non-NaN Hs result
ii = isfinite(ship_oot.Hs_OO) & isfinite(flx_hr.hs_csp_hr);
figure('Position',[300,300,800,800],'Name','hs');
plot([-60,60],[-60,60],'k-','LineWidth',1); hold on;
h5 = plot(flx_hr.hs_csp_hr(ii),ship_oot.Hs(ii),'.',...
          flx_hr.hs_csp_hr(ii),ship_oot.Hs_OO(ii),'.','MarkerSize',10);
grid on; grid minor; xlabel('Hs, BLO, W/m2','FontSize',14); hold off;
ylabel('Hs, W/m2','FontSize',14); axis([-60,60,-60,60]);
legend(h5,'OOT-cov','OOT','Location','northwest','fontsize',14);
pname = fullfile(path_plots,'Fig06_hs_comparison.png');
print(pname,'-dpng');

% 7 hs binned
[hs_blo_oot, hs_oot_mean, hs_oot_med, SD_hs_oot, N_hs_oot] = binave3(flx_hr.hs_csp_hr(ii),ship_oot.Hs(ii),5,-30,30);
[hs_blo_ootOO, hs_oot_meanOO, hs_oot_medOO, SD_hs_ootOO, N_hs_ootOO] = binave3(flx_hr.hs_csp_hr(ii),ship_oot.Hs_OO(ii),5,-30,30);
figure('Position',[300,300,800,800],'Name','hs_bin');
plot([-40,40],[-40,40],'k-','LineWidth',1); hold on;
h7 = errorbar(hs_blo_oot,hs_oot_mean,SD_hs_oot,'square','MarkerSize',12,'LineWidth',2);
h8 = errorbar(hs_blo_ootOO,hs_oot_meanOO,SD_hs_ootOO,'diamond','MarkerSize',12,'LineWidth',2);
grid on; grid minor; xlabel('Hs, BLO, W/m2','FontSize',14); hold off;
ylabel('Hs, OOT, W/m2','FontSize',14); axis([-40,40,-40,40]);
legend([h7,h8],'OOT-cov','OOT','Location','northwest','fontsize',14);
pname = fullfile(path_plots,'Fig07_hs_comparison_bin.png');
print(pname,'-dpng');

% 8 OOT Hs Diff
figure('Name','OOT_Hs_Diff');
plot(ship_oot.Hs(ii),ship_oot.Hs(ii) - ship_oot.Hs_OO(ii),'o'); grid on; grid minor;
xlabel('Hs, OOTcov, W/m2','FontSize',14); ylabel('Hs Difference (OOTcov - OOT)','FontSize',14);
axis([-30,30,-10,10]);
pname = fullfile(path_plots,'Fig08_OOT_Hs_Diff.png');
print(pname,'-dpng');

% 9 Hs diff vs. Hs
figure('Name','OOT_BLO_Diff');
plot(flx_hr.hs_csp_hr(ii),flx_hr.hs_csp_hr(ii) - ship_oot.Hs_OO(ii),'o'); grid on; grid minor;
xlabel('Hs, BLO, W/m2','FontSize',14); ylabel('Hs Difference (BLO-OOT)','FontSize',14);
axis([-30,30,-10,10]);
pname = fullfile(path_plots,'Fig09_OOT_BLO_Diff.png');
print(pname,'-dpng');

% 10 normalized Hs distributions - using surface skin temp from metcity
delT = mc_pop.skin_temp_surface - flx_hr.Ta_Licor_hr;
Ch_ootcv = ship_oot.Hs(ii)./flx_hr.wspd_hr(ii)./delT(ii);
Ch_oot = ship_oot.Hs_OO(ii)./flx_hr.wspd_hr(ii)./delT(ii);
pctile_oot = prctilex(Ch_oot,[2.5 25 50 75 97.5]);
Noot = sum(isfinite(Ch_oot));
oot_str = ['percentiles:  N = ',sprintf('%4i',Noot),char(10),'     2.5     25    50     75     97.5  ',char(10),sprintf('%3.2f  ',pctile_oot)];
Ch_blo = flx_hr.hs_csp_hr(ii)./flx_hr.wspd_hr(ii)./delT(ii);
pctile_blo = prctilex(Ch_blo,[2.5 25 50 75 97.5]);
Nblo = sum(isfinite(Ch_blo));
blo_str = ['percentiles:  N = ',sprintf('%4i',Nblo),char(10),'     2.5     25    50     75     97.5  ',char(10),sprintf('%3.2f  ',pctile_blo)];

figure('Position',[300,300,600,600],'Name','norm_Hs_distributions');
subplot(2,1,1); histogram(Ch_blo,-20:0.25:20); grid on;
ylabel('BLO','FontSize',14);
annotation('textbox',[.2 .5 .3 .3],'String',blo_str,'FitBoxToText','on');
title('Normalized Hs Distributions','FontSize',16);
subtitle('Hs_{norm}  = Hs_{10} / (U_{10} \DeltaT)','FontSize',14);
subplot(2,1,2); histogram(Ch_oot,-20:0.25:20); grid on;
ylabel('OOT','FontSize',14); xlabel('Hs_{norm}','FontSize',14);
annotation('textbox',[.2 0 .3 .3],'String',oot_str,'FitBoxToText','on');
pname = fullfile(path_plots,'Fig10_norm_Hs_distributions.png');
print(pname,'-dpng');

% 11 Hs variance vs Rib binned
[Rib_bins_blo, hs_rib_blo_mean, hs_rib_blo_med, SD_hs_rib_blo, N_hs_rib_blo] = binave3(flx_hr.Rib_hr(ii),flx_hr.hs_csp_hr(ii),0.025,-0.1,0.1);
[Rib_bins_oot, hs_rib_oot_mean, hs_rib_oot_med, SD_hs_rib_oot, N_hs_rib_oot] = binave3(flx_hr.Rib_hr(ii),ship_oot.Hs_OO(ii),0.025,-0.1,0.1);
[Rib_bins_ootcv, hs_rib_ootcv_mean, hs_rib_ootcv_med, SD_hs_rib_ootcv, N_hs_rib_ootcv] = binave3(flx_hr.Rib_hr(ii),ship_oot.Hs(ii),0.025,-0.1,0.1);
figure('Position',[300,300,800,800],'Name','hs_rib_bin');
subtightplot(3,1,1,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
errorbar(Rib_bins_blo,hs_rib_blo_mean,SD_hs_rib_blo,'o','MarkerSize',12,'LineWidth',2);
ylabel('Hs (BLO), W/m2','FontSize',14); grid on; grid minor;
title('Hs vs Ri_b','FontSize',18);
subtightplot(3,1,2,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
errorbar(Rib_bins_ootcv,hs_rib_ootcv_mean,SD_hs_rib_ootcv,'square','MarkerSize',12,'LineWidth',2);
ylabel('Hs (OOT-cov), W/m2','FontSize',14); grid on; grid minor; 
subtightplot(3,1,3,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
errorbar(Rib_bins_oot,hs_rib_oot_mean,SD_hs_rib_oot,'diamond','MarkerSize',12,'LineWidth',2);
ylabel('Hs (OOT), W/m2','FontSize',14); grid on; grid minor; 
xlabel('Ri_b','FontSize',14);
pname = fullfile(path_plots,'Fig11_hs_vs_Rib_bin.png');
print(pname,'-dpng');

% 12 Hs_norm variance vs Rib
[Rib_bins_bloN, hs_rib_blo_meanN, hs_rib_blo_medN, SD_hs_rib_bloN, N_hs_rib_bloN] = binave3(flx_hr.Rib_hr(ii),Ch_blo,0.025,-0.1,0.1);
[Rib_bins_ootN, hs_rib_oot_meanN, hs_rib_oot_medN, SD_hs_rib_ootN, N_hs_rib_ootN] = binave3(flx_hr.Rib_hr(ii),Ch_oot,0.025,-0.1,0.1);
[Rib_bins_ootcvN, hs_rib_ootcv_meanN, hs_rib_ootcv_medN, SD_hs_rib_ootcvN, N_hs_rib_ootcv] = binave3(flx_hr.Rib_hr(ii),Ch_ootcv,0.025,-0.1,0.1);
figure('Position',[300,300,800,800],'Name','hsNorm_rib_bin');
subtightplot(3,1,1,[0.05,0.05],[0.05,0.07],[0.08,0.05]);
errorbar(Rib_bins_bloN,hs_rib_blo_meanN,SD_hs_rib_bloN,'o','MarkerSize',12,'LineWidth',2);
ylabel('Hs_{norm} (BLO)','FontSize',14); ylim([-5,7]); grid on; grid minor;
title('Hs_{norm} vs Ri_b','FontSize',18);
subtitle('Hs_{norm}  = Hs_{10} / (U_{10} \DeltaT)','FontSize',14);
subtightplot(3,1,2,[0.05,0.05],[0.05,0.07],[0.08,0.05]);
errorbar(Rib_bins_ootcvN,hs_rib_ootcv_meanN,SD_hs_rib_ootcvN,'square','MarkerSize',12,'LineWidth',2);
ylabel('Hs_{norm} (OOT-cov)','FontSize',14); ylim([-5,7]);  grid on; grid minor; 
subtightplot(3,1,3,[0.05,0.05],[0.05,0.07],[0.08,0.05]);
errorbar(Rib_bins_ootN,hs_rib_oot_meanN,SD_hs_rib_ootN,'diamond','MarkerSize',12,'LineWidth',2);
ylabel('Hs_{norm} (OOT)','FontSize',14); ylim([-5,7]);  grid on; grid minor; 
xlabel('Ri_b','FontSize',14);
pname = fullfile(path_plots,'Fig12_hsNorm_vs_Rib_bin.png');
print(pname,'-dpng');

% 13 compare OOT and BLO CO2 flux
ii = isfinite(flx_hr.Fwco2_csp_mM_hr) & isfinite(ship_oot.Fco2_OO);
co2_coefs_13min = polyfit(flx_hr.Fwco2_csp_mM_hr(ii),ship_oot.Fco2_OO(ii),1);
xref = -20:1:20;
fco2_13min_line = polyval(co2_coefs_13min,xref);
figure; plot(xref,xref,'k-'); hold on;
plot(flx_hr.Fwco2_csp_mM_hr(ii),ship_oot.Fco2_OO(ii),'bo','MarkerSize',8);
plot(xref,fco2_13min_line,'r-','LineWidth',1.5); hold off; grid on; grid minor;
axis([-20,20,-20,20]);
ylabel('F_{co2} OOT, mM/m2/d','fontsize',14);
xlabel('F_{co2} PSL, mM/m2/d','fontsize',14);
pname = fullfile(path_plots,'Fig13_fco2OOT_vs_fco2PSL.png');
print(pname,'-dpng');

% 14 compare OOT versions of CO2 flux
ii = isfinite(ship_oot.Fco2) & isfinite(ship_oot.Fco2_OO) & ...
              ship_oot.Fco2 < 20 & ship_oot.Fco2 > -20;
co2_coefs_30min = polyfit(ship_oot.Fco2(ii),ship_oot.Fco2_OO(ii),1);
fco2_30min_line = polyval(co2_coefs_30min,xref);
figure; plot(xref,xref,'k-'); hold on;
plot(ship_oot.Fco2(ii),ship_oot.Fco2_OO(ii),'bo','MarkerSize',8);
plot(xref,fco2_30min_line,'r-','LineWidth',1.5); hold off; grid on; grid minor;
axis([-20,20,-20,20]);
ylabel('F_{co2} OOT, mM/m2/d','fontsize',14);
xlabel('F_{co2} 30min, mM/m2/d','fontsize',14);
pname = fullfile(path_plots,'Fig14_fco2OOT_vs_fco2_30min.png');
print(pname,'-dpng');
