% a15_weekly_fluxes_merged_ship.m
% compute weighted weekly means for gas fluxes
% compare co2 flux with project means from other cruises

% saves the following weekly mean variables in the flx_wk table
%{
1   datenum               matlab datenum
2   yyyy                  year
3   mm                    month
4   dd                    day of month
5   hh                    hour of day
6   Fwco2_csp_mM_wk       mean CO2 flux from cospectrum, mM/m2/d
7   SEwco2_csp_mM_wk      std error of weekly mean Fwco2_csp_mM_wk
8   Fwch4_csp_uM_wk       mean CH4 flux from cospectrum, uM/m2/d
9   SEwch4_csp_uM_wk      std error of weekly mean Fwch4_csp_uM_wk
10  Fwo3_csp_uM_wk        mean O3 flux from cospectrum, uM/m2/d
11  SEwo3_csp_uM_wk       std error of weekly mean Fwo3_csp_uM_wk
12  Fwdms_csp_uM_wk       mean DMS flux from cospectrum, uM/m2/d
13  SEwdms_csp_uM_wk      std error of weekly mean Fwdms_csp_uM_wk
14  Vd_o3_wk              mean O3 deopsition velocity, cm/s
15  hs_csp_wk             mean sensible heat flux from cospectrum, W/m2
16  SEhs_csp_wk           std error of weekly mean hs_csp_wk
17  Ngood_co2             number of 'good' hours of Fco2 data per week
18  Ngood_ch4             number of 'good' hours of Fch4 data per week
19  Ngood_o3              number of 'good' hours of Fo3 data per week
20  Ngood_dms             number of 'good' hours of Fdms data per week
%}

%% initialize
clear
close all
plotit = 1;

root_data_path = '/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019';
path_prog = fullfile(root_data_path,'_matlab');
plot_path = fullfile(root_data_path,'_processed','ship_flux_V3','flux_merged','plots','weekly_means');
mkdir(plot_path);

path(pathdef);
cd(path_prog);
addpath(genpath(path_prog));

% load flx_hr table and comparison cruise data tables
load(fullfile(root_data_path,'_processed','ship_flux_V3','flux_merged','MOSAiC_ship_flux_co2_ch4_o3_dms_hr.mat'));
load(fullfile(root_data_path,'_processed','ship_flux_V3','comparison_cruises','comp_cruises.mat'));

%% compute mean fluxes for comparison cruises
% AMT28 / AMT29: Eq.S.Atl., 0° to -10° lat
% TORERO: Eq.E.Pac. Cold Tongue, <0° lat
% TORERO: Eq.E.Pac. S. of Hawaii, 13° to 19° lat
% HiWinGS: Lab.Sea, Oct. 2013, 51-55° lat
% JR18006: Barents Sea Jun-Jul 2019, > 75° lat
% JR18007: Fram Strait Aug-Sept 2019, > 75° lat

amt_all = [AMT28.co2flux(AMT28.lat<0 & AMT28.lat>-10);AMT29.co2flux(AMT29.lat<0 & AMT29.lat>-10)];
amtFc02_mean = nanmean1(amt_all);
amtFc02_std = nanstd1(amt_all);
dt_amt = datetime(2019,11,1);

torN_mean = nanmean1(TORERO.CO2flux_uM(TORERO.lat>13));
torN_std = nanstd1(TORERO.CO2flux_uM(TORERO.lat>13));
dt_torN = datetime(2020,1,15);

torS_mean = nanmean1(TORERO.CO2flux_uM(TORERO.lat<0));
torS_std = nanstd1(TORERO.CO2flux_uM(TORERO.lat<0));
dt_torS = datetime(2020,2,15);

hwco2_uM = HiWinGSCO2.CO2flux_ppm_m_s.*86400.*HiWinGSCO2.Press_Pa./8.31451./(HiWinGSCO2.Ta_C + 273.15)./1000;
hwco2_mean = nanmean1(hwco2_uM(HiWinGSMetNavC35.lat>51 & HiWinGSMetNavC35.lat<55));
hwco2_std = nanstd1(hwco2_uM(HiWinGSMetNavC35.lat>51 & HiWinGSMetNavC35.lat<55));
dt_hw = datetime(2020,10,1);

JRbar_mean = nanmean1(JR18006.co2flux(JR18006.lat>75));
JRbar_std = nanstd1(JR18006.co2flux(JR18006.lat>75));
dt_JRbar = datetime(2020,7,1);

JRfra_mean = nanmean1(JR18007.co2flux(JR18007.lat>75));
JRfra_std = nanstd1(JR18007.co2flux(JR18007.lat>75));
dt_JRfra = datetime(2020,8,1);

%% set up for weekly averages
dnum_bins = datenum(2019,10,27):7:datenum(2020,10,4); dnum_bins = dnum_bins(:); % weekly bin edges

figure; histogram(flx_hr.Rib_hr,-0.3:0.02:0.5); grid on; grid minor;
xlabel('Ri bulk, 20m','fontsize',14); title('MOSAiC Hourly Ri distribution (20m), Ship','fontsize',18);
if plotit
    ppath = fullfile(plot_path,'Ri_dist_ship.png');
    print(ppath,'-dpng');
end

%% co2 flux
% filter for at least 3 good segments per hour and Rib < 0.1;
% hours where Rib is NaN (no sfc temp data) are passed as 'good'
ii = flx_hr.sflag_co2_sum_hr>=3 & ~(isfinite(flx_hr.Rib_hr) & flx_hr.Rib_hr>0.1);
Ngood_co2 = interval_sum(flx_hr.datenum(ii),ii(ii),dnum_bins);
Ngood_co2 = Ngood_co2(:,2);

[fco2,se_fco2,t_wk] = interval_weighted_avg(flx_hr.datenum(ii),flx_hr.Fwco2_csp_mM_hr(ii),...
    flx_hr.SEwco2_csp_mM_hr(ii),dnum_bins);
dnum_wk = t_wk + 7/2; % mid-week timestamp
dt_wk = datetime(dnum_wk,'convertfrom','datenum');
x_lims = [datetime(2019,10,1),datetime(2020,10,15)];

xx = Ngood_co2>6;  % at least 6 hours of good data for a weekly mean...
figure('position',[100,750,1000,500],'Name','fco2_weekly');
plot(x_lims,[0,0],'k-','linewidth',1.5); grid on; grid minor; hold on;
errorbar(dt_wk(xx),fco2(xx),se_fco2(xx),'o','markersize',8,'linewidth',1.5);
hold off; xlim(x_lims); datetick('x', 'mmm','keeplimits');
xlabel('Month, 2019-2020','fontsize',14);
ylabel('F_{co2}  mM/m^2/d \pm SEM','fontsize',14);
title('MOSAiC Weekly Mean CO_2 Turbulent Flux, Ship','fontsize',18);
if plotit
    ppath = fullfile(plot_path,'fco2_weekly_ship.png');
    print(ppath,'-dpng');
end

%% co2 flux - compared to other cruises

fco2_sigma = interval_std(flx_hr.datenum(ii),flx_hr.Fwco2_csp_mM_hr(ii),dnum_bins);
fco2_sigma = fco2_sigma(:,2);

% compute mean flux for winter/spring and summer seasons
jj = dt_wk < datetime(2020,06,01); % winter/spring dates
kk = xx & jj;  % good co2 fluxes in winter/spring
mm = xx & ~jj; % good co2 fluxes in summer
fco2_winter = nanmean1(fco2(kk));
fco2_summer = nanmean1(fco2(mm));
dt_winter = [datetime(2019,11,1),datetime(2020,5,1)];
dt_summer = [datetime(2020,6,1),datetime(2020,10,1)];

dim = [0.2 0.2 0.3 0.3]; % text box dimensions
str = {['fco2 winter = ',sprintf('%3.2f',fco2_winter),' mM/m^2/d'],...
       ['fco2 summer = ',sprintf('%3.2f',fco2_summer),' mM/m^2/d']};

figure('position',[100,750,1000,500]);
plot(x_lims,[0,0],'k-','linewidth',1.5); grid on; grid minor; hold on;
plot(dt_winter,[fco2_winter,fco2_winter],'k--','LineWidth',1.5);
plot(dt_summer,[fco2_summer,fco2_summer],'k--','LineWidth',1.5);
errorbar(dt_wk(xx),fco2(xx),fco2_sigma(xx),'o','markersize',8,'linewidth',1.5,...
    'Color','#D95319');
p1=errorbar(dt_amt,amtFc02_mean,amtFc02_std,'rd','markersize',12,'linewidth',2.5);
p2=errorbar(dt_torN,torN_mean,torN_std,'bd','markersize',12,'linewidth',2.5);
p3=errorbar(dt_torS,torS_mean,torS_std,'gd','markersize',12,'linewidth',2.5);
p4=errorbar(dt_hw,hwco2_mean,hwco2_std,'kd','markersize',12,'linewidth',2.5);
p5=errorbar(dt_JRbar,JRbar_mean,JRbar_std,'md','markersize',12,'linewidth',2.5);
p6=errorbar(dt_JRfra,JRfra_mean,JRfra_std,'cd','markersize',12,'linewidth',2.5);
hold off; xlim(x_lims); datetick('x', 'mmm','keeplimits');
xlabel('Month, 2019-2020','fontsize',14); ylabel('F_{co2}  mM/m^2/d \pm 1\sigma','fontsize',14);
title('CO_2 Turbulent Flux +/- 1\sigma, Ship, Comparison Cruises','fontsize',18);
annotation('textbox',dim,'String',str,'FitBoxToText','on','FontSize',14);
legend([p1,p2,p3,p4,p5,p6],'AMT Oct 2018/Nov 2019 Eq.S.Atl, 0 to -10° lat, Dong et al. 2021',...
    'TORERO Jan/Feb 2012, C.Pac., 13-19°N lat, Blomquist et al. 2014',...
    'TORERO Jan/Feb 2012, E.Pac. Cold Tongue, 0-9°S lat, Blomquist et al. 2014',...
    'HiWinGS Oct 2013, Labrador Sea, 51-55° lat, Blomquist et al. 2017',...
    'JCR Jun/Jul 2019, Barents Sea, >75°N, Dong et al. 2021',...
    'JCR Aug/Sept 2019, Fram Strait, >75°N, Dong et al. 2021','location','southwest','fontsize',14);
if plotit
    ppath = fullfile(plot_path,'fco2_weekly_comp_cruises_1sig_ship.png');
    print(ppath,'-dpng');
end

%% ch4 flux
% filter for at least 3 good segments per hour and Rib < 0.1;
% hours where Rib is NaN (no sfc temp data) are passed as 'good'
jj = flx_hr.sflag_ch4_sum_hr>=3 & ~(isfinite(flx_hr.Rib_hr) & flx_hr.Rib_hr>0.1);
Ngood_ch4 = interval_sum(flx_hr.datenum(jj),jj(jj),dnum_bins); Ngood_ch4 = Ngood_ch4(:,2);

[fch4,se_fch4,~] = interval_weighted_avg(flx_hr.datenum(jj),flx_hr.Fwch4_csp_uM_hr(jj),...
    flx_hr.SEwch4_csp_uM_hr(jj),dnum_bins);

xx = Ngood_ch4>6;
figure('position',[100,750,1000,500]);
plot(x_lims,[0,0],'k-','linewidth',1.5); grid on; grid minor; hold on;
errorbar(dt_wk(xx),fch4(xx),se_fch4(xx),'o','markersize',8,'linewidth',1.5); hold off;
xlim(x_lims); datetick('x', 'mmm','keeplimits');
xlabel('Month, 2019-2020','fontsize',14); ylabel('F_{ch4}  \muM/m^2/d \pm SEM','fontsize',14);
title('MOSAiC Weekly Mean CH_4 Turbulent Flux, Ship System','fontsize',18);
if plotit
    ppath = fullfile(plot_path,'fch4_weekly_ship.png');
    print(ppath,'-dpng');
end

%% dms flux
% filter for at least 3 good segments per hour and Rib < 0.1;
% hours where Rib is NaN (no sfc temp data) are passed as 'good'
jj = flx_hr.sflag_dms_sum_hr>=3 & ~(isfinite(flx_hr.Rib_hr) & flx_hr.Rib_hr>0.1);
Ngood_dms = interval_sum(flx_hr.datenum(jj),jj(jj),dnum_bins); Ngood_dms = Ngood_dms(:,2);

[fdms,se_fdms,~] = interval_weighted_avg(flx_hr.datenum(jj),flx_hr.Fwdms_csp_uM_hr(jj),...
    flx_hr.SEwdms_csp_uM_hr(jj),dnum_bins);

x_lims_dms = [datetime(2020,7,1),datetime(2020,10,15)];

xx = Ngood_dms>6;
figure('position',[100,750,1000,500]);
plot(x_lims_dms,[0,0],'k-','linewidth',1.5); grid on; grid minor; hold on;
errorbar(dt_wk(xx),fdms(xx),se_fdms(xx),'o','markersize',8,'linewidth',1.5); hold off;
xlim(x_lims_dms); datetick('x', 'mmm','keeplimits');
xlabel('Month, 2019-2020','fontsize',14); ylabel('F_{dms}  \muM/m^2/d \pm SEM','fontsize',14);
title('MOSAiC Weekly Mean DMS Turbulent Flux, Ship System','fontsize',18);
if plotit
    ppath = fullfile(plot_path,'fdms_weekly_ship.png');
    print(ppath,'-dpng');
end

%% o3 flux
% filter for at least 3 good segments per hour and Rib < 0.1;
% hours where Rib is NaN (no sfc temp data) are passed as 'good'
kk = flx_hr.sflag_o3_sum_hr>=3 & ~(isfinite(flx_hr.Rib_hr) & flx_hr.Rib_hr>0.1);
Ngood_o3 = interval_sum(flx_hr.datenum(kk),kk(kk),dnum_bins); Ngood_o3 = Ngood_o3(:,2);

[fo3,se_fo3,~] = interval_weighted_avg(flx_hr.datenum(kk),flx_hr.Fwo3_csp_uM_hr(kk),...
    flx_hr.SEwo3_csp_uM_hr(kk),dnum_bins);
Vd = interval_avg(flx_hr.datenum(kk),flx_hr.Vd_o3_hr(kk),dnum_bins);
Vd = Vd(:,2);

xx = Ngood_o3>6;
figure('position',[100,750,1000,500]);
plot(x_lims,[0,0],'k-','linewidth',1.5); grid on; grid minor; hold on;
errorbar(dt_wk(xx),fo3(xx),se_fo3(xx),'o','markersize',8,'linewidth',1.5); hold off;
xlim(x_lims); datetick('x', 'mmm','keeplimits');
xlabel('Month, 2019-2020','fontsize',14); ylabel('F_{o3}  \muM/m^2/d \pm SEM','fontsize',14);
title('MOSAiC Weekly Mean O_3 Turbulent Flux, Ship System','fontsize',18);
ppath = fullfile(plot_path,'fo3_weekly_ship.png');
print(ppath,'-dpng');

figure('position',[100,750,1000,500]);
plot(x_lims,[0,0],'k-','linewidth',1.5); grid on; grid minor; hold on;
plot(dt_wk(xx),Vd(xx),'o','markersize',8,'linewidth',1.5); hold off;
xlim(x_lims); datetick('x', 'mmm','keeplimits');
xlabel('Month, 2019-2020','fontsize',14); ylabel('V_d (O_3) cm/s','fontsize',14);
title('MOSAiC Weekly Mean O_3 Deposition Velocity, Ship System','fontsize',18);
if plotit
    ppath = fullfile(plot_path,'Vdo3_weekly_ship.png');
    print(ppath,'-dpng');
end

%% sensible heat
% filter for at least 3 good segments per hour and Rib < 0.1;
% hours where Rib is NaN (no sfc temp data) are passed as 'good'
hh = flx_hr.sflag_sum_hr>=3 & ~(isfinite(flx_hr.Rib_hr) & flx_hr.Rib_hr>0.1);

[hs,se_hs,~] = interval_weighted_avg(flx_hr.datenum(hh),flx_hr.hs_csp_hr(hh),...
    flx_hr.SEhs_csp_hr(hh),dnum_bins);

figure('position',[100,750,1000,500]);
plot(x_lims,[0,0],'k-','linewidth',1.5); grid on; grid minor; hold on;
errorbar(dt_wk,hs,se_hs,'o','markersize',8,'linewidth',1.5); hold off;
xlim(x_lims); datetick('x', 'mmm','keeplimits');
xlabel('Month, 2019-2020','fontsize',14); ylabel('Sensible Heat Flux, W/m^2','fontsize',14);
title('MOSAiC Weekly Mean Sensible Heat Flux, Ship System','fontsize',18);
if plotit
    ppath = fullfile(plot_path,'hs_weekly_ship.png');
    print(ppath,'-dpng');
end

%% save table of weekly mean results
dv = datevec(t_wk);
vars = {'datenum','yyyy','mm','dd','hh','Fwco2_csp_mM_wk','SEwco2_csp_mM_wk','Fwch4_csp_uM_wk',...
    'SEwch4_csp_uM_wk','Fwo3_csp_uM_wk','SEwo3_csp_uM_wk','Vd_o3_wk','hs_csp_wk','SEhs_csp_wk',...
    'Fwdms_csp_uM_wk','SEwdms_csp_uM_wk','Ngood_co2','Ngood_ch4','Ngood_o3','Ngood_dms'};
flx_wk = array2table([t_wk,dv(:,1:4),fco2,se_fco2,fch4,se_fch4,fo3,se_fo3,Vd,hs,se_hs,...
    fdms,se_fdms,Ngood_co2,Ngood_ch4,Ngood_o3,Ngood_dms],'VariableNames',vars);

sfl = fullfile(root_data_path,'_processed','ship_flux_V3','flux_merged','MOSAiC_ship_flux_co2_ch4_o3_dms_wk.mat');
save(sfl,'flx_wk');
