% b22_daily_fluxes_metcity.m
% compute weighted daily means for met city gas fluxes
% compare co2 flux with project means from other cruises

% saves the following weekly mean variables in the mc_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_dd       mean CO2 flux from cospectrum, mM/m2/d
7   SEwco2_csp_mM_dd      std error of weekly mean Fwco2_csp_mM_dd
8   Fwch4_csp_uM_dd       mean CH4 flux from cospectrum, uM/m2/d
9   SEwch4_csp_uM_dd      std error of weekly mean Fwch4_csp_uM_dd
10  hs_csp_dd             mean sensible heat flux from cospectrum, W/m2
11  SEhs_csp_dd           std error of weekly mean hs_csp_dd
12  Ngood_co2             number of 'good' hours of Fco2 data per day
13  Ngood_ch4             number of 'good' hours of Fch4 data per day
%}

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

plotit = 1;

root_data_path = '/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019';

path(pathdef);
path_prog = fullfile(root_data_path,'_matlab');
plot_path = fullfile(root_data_path,'_processed','metcity_flux_V3','flux_merged','plots','daily_means');
mkdir(plot_path);

cd(path_prog);
addpath(genpath(path_prog));

% load mc_flx_hr table and comparison cruise data tables
load(fullfile(root_data_path,'_processed','metcity_flux_v3','flux_merged','MOSAiC_metcity_flux_co2_ch4_hr.mat'));
load(fullfile(root_data_path,'_processed','ship_flux_V2','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 weekly avg bins
dnum_bins = datenum(2019,10,27):1:datenum(2020,10,4); dnum_bins = dnum_bins(:); % daily bin edges

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

%% co2 flux
% filter results for cases with 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 = mc_flx_hr.sflag_co2_sum_hr>=3 & ~(isfinite(mc_flx_hr.Rib_hr) & mc_flx_hr.Rib_hr>0.1);
Ngood_co2 = interval_sum(mc_flx_hr.datenum(ii),ii(ii),dnum_bins);
Ngood_co2 = Ngood_co2(:,2);

[fco2,se_fco2,t_dd] = interval_weighted_avg(mc_flx_hr.datenum(ii),mc_flx_hr.Fwco2_csp_mM_hr(ii),...
    mc_flx_hr.SEwco2_csp_mM_hr(ii),dnum_bins);
dnum_dd = t_dd + 12/24; % mid-day timestamp
dt_dd = datetime(dnum_dd,'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_daily');
plot(x_lims,[0,0],'k-','linewidth',1.5); grid on; grid minor; hold on;
errorbar(dt_dd(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); ylim([-1.5,2.5]);
ylabel('F_{co2}  mM/m^2/d \pm SEM','fontsize',14);
title('MOSAiC Daily Mean CO_2 Turbulent Flux, Met City','fontsize',18);
if plotit
    ppath = fullfile(plot_path,'fco2_daily_metcity.png');
    print(ppath,'-dpng');
    fpath = fullfile(plot_path,'fco2_daily_metcity.fig');
    savefig(fpath);
end

%% co2 flux - comp cruises

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

figure('position',[100,750,1000,500],'Name','fco2_comp_cruises');
plot(x_lims,[0,0],'k-','linewidth',1.5); grid on; grid minor; hold on;
errorbar(dt_dd(xx),fco2(xx),fco2_sigma(xx),'o','markersize',8,'linewidth',1.5);
p1=errorbar(dt_amt,amtFc02_mean,amtFc02_std,'rd','markersize',12,'linewidth',2);
p2=errorbar(dt_torN,torN_mean,torN_std,'bd','markersize',12,'linewidth',2);
p3=errorbar(dt_torS,torS_mean,torS_std,'gd','markersize',12,'linewidth',2);
p4=errorbar(dt_hw,hwco2_mean,hwco2_std,'kd','markersize',12,'linewidth',2);
p5=errorbar(dt_JRbar,JRbar_mean,JRbar_std,'md','markersize',12,'linewidth',2);
p6=errorbar(dt_JRfra,JRfra_mean,JRfra_std,'cd','markersize',12,'linewidth',2);
hold off; xlim(x_lims); %datetick('x', 'mmm','keeplimits');
xlabel('Month, 2019-2020','fontsize',12); ylabel('F_{co2}  mM/m^2/d \pm 1\sigma','fontsize',14);
title('CO_2 Turbulent Flux +/- 1\sigma, MetCity, Comparison Cruises','fontsize',18);
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_daily_comp_cruises_1sig_metcity.png');
    print(ppath,'-dpng');
    fpath = fullfile(plot_path,'fco2_daily_comp_cruises_1sig_metcity.fig');
    savefig(fpath);
end

%% ch4 flux
% filter results for cases with 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 = mc_flx_hr.sflag_ch4_sum_hr>=3 & ~(isfinite(mc_flx_hr.Rib_hr) & mc_flx_hr.Rib_hr>0.1);
Ngood_ch4 = interval_sum(mc_flx_hr.datenum(jj),jj(jj),dnum_bins);
Ngood_ch4 = Ngood_ch4(:,2);

[fch4,se_fch4,~] = interval_weighted_avg(mc_flx_hr.datenum(jj),mc_flx_hr.Fwch4_csp_uM_hr(jj),...
    mc_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_dd(xx),fch4(xx),se_fch4(xx),'o','markersize',8,'linewidth',1.5); hold off;
ylim([-20,25]); 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 Daily Mean CH_4 Turbulent Flux, MetCity','fontsize',18);
if plotit
    ppath = fullfile(plot_path,'fch4_daily_metcity.png');
    print(ppath,'-dpng');
    fpath = fullfile(plot_path,'fch4_daily_metcity.fig');
    savefig(fpath);
end

%% sensible heat
hh = mc_flx_hr.sflag_sum_hr>=3;
[hs,se_hs,~] = interval_weighted_avg(mc_flx_hr.datenum(hh),mc_flx_hr.hs_csp_hr(hh),...
    mc_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_dd,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 Daily Mean Sensible Heat Flux, MetCity','fontsize',18);
if plotit
    ppath = fullfile(plot_path,'hs_daily_metcity.png');
    print(ppath,'-dpng');
    fpath = fullfile(plot_path,'hs_daily_metcity.fig');
    savefig(fpath);
end

%% save table of daily mean results
dv = datevec(t_dd);
vars = {'datenum','yyyy','mm','dd','hh','Fwco2_csp_mM_dd','SEwco2_csp_mM_dd','Fwch4_csp_uM_dd',...
    'SEwch4_csp_uM_dd','hs_csp_dd','SEhs_csp_dd','Ngood_co2','Ngood_ch4'};
mc_flx_dd = array2table([t_dd,dv(:,1:4),fco2,se_fco2,fch4,se_fch4,hs,se_hs,Ngood_co2,Ngood_ch4],'VariableNames',vars);

sfl = fullfile(root_data_path,'_processed','metcity_flux_v3','flux_merged','MOSAiC_metcity_flux_co2_ch4_daily.mat');
save(sfl,'mc_flx_dd');

