% b20_mc_co2_flux_comparison
% compare Licor and Picarro CO2 fluxes at MetCity tower
% compute Webb correction factors and compare to POP Webb term

clear
close all

root_data_path = '/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019';
path_plots = fullfile(root_data_path,'_processed','metcity_flux_v3','co2_flux_comparison_plots');
path_prog = fullfile(root_data_path,'_matlab');
path(pathdef);
cd(path_prog);
addpath(genpath(path_prog));

% load mc_flux table with 10-min POP flux data
load('/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019/_processed/metcity_tower_level2_v4/metcity_10min_flux_all.mat');
dt = datetime(mc_flux.datenum,'convertfrom','datenum');

% constants and bulk met parameters
Rgas = 287.1;   % J/kg/K  dry air gas constant
cpa  = 1004.67; % J/kg/K  dry air heat capacity
T = mc_flux.temp_2m;
T_K = T + 273.15;
Le   = (2.501 - 0.00237*0)*1e6; % latent heat evap for pure water at 0°C, J/kg
RH = mc_flux.rh_2m;
Pmb = mc_flux.atmos_pressure_2m; % Pair in mb
mu = 1.6077; % molar mass ratio air/water = 28.9647/18.01528, unitless
[Q,Pv]  = qsat26air(T,Pmb,RH);  % specific humidity Q in g/kg, Pv in mb
Q = Q/1000; % kg/kg units
rhodry = (Pmb-Pv)*100./(Rgas*T_K);        % rho air dry, kg/m3
rhoa = Pmb*100./(Rgas*(T_K).*(1+0.61*Q)); % rho moist air, kg/m3

% CO2 flux conversion factor : mg/m2/s to mmol/m2/day
fconvert = 86400/44; % 86400 s/day x 1 mmol/44mg;
co2_mg_m3 = mc_flux.co2_licor*44; % co2 concentration : convert mmol/m3 to mg/m3
% co2_ppm = mmol/m3 * R * T_K / P(kPa)
co2_ppm = (mc_flux.co2_licor*8.314463.*T_K./(Pmb/10)); % co2 molar mixing ratio in ppm
fco2 = mc_flux.CO2_flux*fconvert;        % POP-computed CO2 flux
fco2w = mc_flux.CO2_flux_Webb*fconvert;  % POP-computed CO2 flux, Webb adjusted

% Webb method 1) compute wbar with C36 bulk model
% BE SURE TO USE ICE OUTPUT VALUEs, e.g. wbari, not wbar
% wbari = 1.61*hlbi./Le./(1+1.61*Q)./rhoa+hsbi./rhoa./cpa./ta;
%       = J/s/m2 * kg/J * kg/kg * m3/kg + J/s/m2 * m3/kg * kg*K/J * 1/K
%       = m/s + m/s
% run coare36vn_ice_zrf_f(u,zu,t,zt,rh,zq,P,ts,Rs,Rl,lat,zi,rain,Ss,cp,sigH,zrf_u,zrf_t,zrf_q,zoi_0)
[~,C36] = coare36vn_ice_zrf_f(mc_flux.wspd_vec_mean_2m,ones(height(mc_flux),1)*2,... % u,zu
                             mc_flux.temp_2m,ones(height(mc_flux),1)*2,...           % t,zt
                             mc_flux.rh_2m,ones(height(mc_flux),1)*2,...             % rh,zq
                             mc_flux.atmos_pressure_2m,...                           % P
                             mc_flux.skin_temp_surface,...                           % ts 
                             mc_flux.down_short_hemisp ,...                          % Rs
                             mc_flux.down_long_hemisp,...                            % Rl
                             mc_flux.lat_tower,100,zeros(height(mc_flux),1),...      % lat, zi, rain
                             ones(height(mc_flux),1)*32,...                          % Ss
                             NaN,NaN,10,10,10,...                                    % cp, sigH, zrf_u, zrf_t, zrf_q
                             ones(height(mc_flux),1)*1.5e-3);                        % zoi_0
% now compute Webb term from wbar and co2 concenttation in mg units
Webb_bulk = co2_mg_m3.*C36.wbari;  % mg/m3 * m/s = mg/m2/s

% convert POP w'Ts' to w'T'
% method 1) use usb/qsb to convert WTs to WT: wT = wTs - (0.51 * ta_K * -(usb * qsb*1e-3));  qsb to kg/kg units
% WTb = mc_flux.WT_csp_2m + (0.51.*T_K.*C36.usri.*C36.qsri.*1e-3);
% method 2) use observed w'q': WT = WTs - (0.51 * ta_K * wq)
% NEED TO USE UPDATED wq WITH CORRECTED q...
wq = mc_flux.Wq_csp;   % kg/kg m/s
% WT = mc_flux.WT_csp_2m - (0.51.*T_K.*wq);

% OR, just use WTs and ignore moisture effect
WT = mc_flux.WT_csp_2m;

% Webb method 2) compute moisture and thermal web terms from observed w'q' and w'T'
% Webb water vapor dilution term from observed wq
% (mg/m3 / mg/m3) * mg/m2/s = mg/m2/s
Webb_v = mu*co2_mg_m3./(rhodry*1e6).*wq*1e6; % 1e6 = kg to mg
% Webb thermal expansion term from observed WT
% (1 + mu*(mg/m3 / mg/m3)) * (K*m/s)/K * (mg/m3) = mg/m2/s
Webb_t = (1 + mu*(co2_mg_m3./(rhodry*1e6))).*(WT./T_K).*co2_mg_m3;

% now compute Webb-adjusted fco2 in mmol/m2/d units
fco2w_bulk_mg = mc_flux.CO2_flux + Webb_bulk;         % using Webb method 1
fco2w_bulk = fco2w_bulk_mg*fconvert;                  % mmol/m2/d units
fco2w_recalc_mg = mc_flux.CO2_flux + Webb_v + Webb_t; % using Webb method 2
fco2w_recalc = fco2w_recalc_mg*fconvert;              % mmol/m2/d units


% now some plots...
% selection filter for good 2m fluxes
ff = mc_flux.turbulence_2m_qc == 0; % & dt>datetime(2020,6,1);

%1
figure; plot(dt(ff),fco2(ff),'.',dt(ff),fco2w(ff),'.','markersize',6);
grid on; grid minor; legend('raw','webb','fontsize',12,'fontsize',12,'location','south');
ylabel('CO_2 flux, mmol/m^2/d'); ylim([-500,500]); ax = gca; ax.FontSize = 12;
title('CO_2 Flux, POP Results','FontSize',16);
ppath = fullfile(path_plots,'co2_flux_timeseries_POP.png');
print(ppath,'-dpng');

%2
figure('Position',[300,300,1500,800]);
subtightplot(1,1,1,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
plot(dt(ff),Webb_t(ff)*fconvert,'.',dt(ff),Webb_v(ff)*fconvert,'.',...
     'MarkerSize',12); grid on; grid minor;
xlim([datetime(2019,10,26),datetime(2020,9,20)]); ax = gca; ax.FontSize = 14;
legend('webb-T_{obs}','webb-v_{obs}','fontsize',12,'location','southeast');
legendmarkeradjust(20); ylabel('Webb Terms, mmol/m^2/d'); ylim([-200,200]); 
title('CO_2 Flux Webb Adjustment Terms','FontSize',20);
ppath = fullfile(path_plots,'co2_flux_webb_terms.png');
print(ppath,'-dpng');

%3
figure('Position',[300,300,600,800]);
subtightplot(2,1,1,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
plot([-400,400],[-400,400],'k-','LineWidth',1.5); hold on;
plot(fco2w(ff),fco2w_recalc(ff),'.','markersize',6); hold off;
axis([-400,400,-400,400]); grid on; grid minor; ax = gca; ax.FontSize = 12;
ylabel('fco2-webb mmol/m2^/d Recalc'); xlabel('fco2-webb, mmol/m^2/d POP');
title('Webb Corrected LI-7500 CO_2 Flux Comparison','FontSize',16);
subtightplot(2,1,2,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
plot([-400,400],[-400,400],'k-','LineWidth',1.5); hold on;
plot(fco2w(ff),fco2w_bulk(ff),'.','markersize',6); hold off;
axis([-400,400,-400,400]); grid on; grid minor; ax = gca; ax.FontSize = 12;
ylabel('fco2-webb mmol/m^2/d Bulk'); xlabel('fco2-webb, mmol/m^2/d POP');
ppath = fullfile(path_plots,'co2_flux_webb_recalc_vs_POP.png');
print(ppath,'-dpng');


% compute filtered hourly average fluxes
% exclude periods when 2m turbulence is bad
% load Picarro flux table mc_flx_hr
load('/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019/_processed/metcity_flux_v3/flux_merged/MOSAiC_metcity_flux_co2_ch4_hr.mat');
bins = mc_flx_hr.datenum;
bins = [bins;bins(end)+1/24];
dt_hr = datetime(mc_flx_hr.datenum,'ConvertFrom','datenum');

fco2w_hr = interval_avg(mc_flux.datenum(ff),fco2w(ff),bins); fco2w_hr = fco2w_hr(:,2);
fco2w_recalc_hr = interval_avg(mc_flux.datenum(ff),fco2w_recalc(ff),bins); fco2w_recalc_hr = fco2w_recalc_hr(:,2);
fco2w_bulk_hr = interval_avg(mc_flux.datenum(ff),fco2w_bulk(ff),bins); fco2w_bulk_hr = fco2w_bulk_hr(:,2);

%4
figure; plot(dt_hr,fco2w_hr,'.',dt_hr,fco2w_recalc_hr,'.',dt_hr,fco2w_bulk_hr,'.',dt_hr,mc_flx_hr.Fwco2_csp_mM_hr,'.','markersize',8);
grid on; grid minor; ax = gca; ax.FontSize = 12;
legend('pop-Webb','recalc-Webb','bulk-Webb','crds','fontsize',12,'location','south');
ylabel('CO_2 flux, mmol/m^2/d'); ylim([-100,100]);
title('CO_2 Flux','FontSize',16);
ppath = fullfile(path_plots,'co2_flux_timeseries_licor_crds.png');
print(ppath,'-dpng');

%5
figure('Position',[300,300,1500,800]);
subtightplot(2,1,1,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
plot(dt_hr,fco2w_recalc_hr,'.',dt_hr,mc_flx_hr.Fwco2_csp_mM_hr,'.','markersize',12);
legend('LI-7500','CRDS','fontsize',12,'location','south');  legendmarkeradjust(20);
ylabel('CO_2 flux, mmol/m^2/d','FontSize',15); ylim([-100,100]);grid on; grid minor;
xlim([datetime(2019,10,26),datetime(2020,9,20)]); ax = gca; ax.FontSize = 14;
title('Met CIty CO_2 Flux, LI-7500 Open Path and CRDS Closed Path','FontSize',20);
subtightplot(2,1,2,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
plot(dt_hr,fco2w_recalc_hr,'.',dt_hr,mc_flx_hr.Fwco2_csp_mM_hr,'.','markersize',16);
legend('LI-7500','CRDS','fontsize',12,'location','south'); legendmarkeradjust(20);
ylabel('CO_2 flux, mmol/m^2/d','FontSize',15); ylim([-10,10]);grid on; grid minor;
xlim([datetime(2019,10,26),datetime(2020,9,20)]); ax = gca; ax.FontSize = 14;
ppath = fullfile(path_plots,'co2_flux_hourly_timeseries_licor_crds.png');
print(ppath,'-dpng');

%6
hl_hr = interval_avg(mc_flux.datenum(ff),mc_flux.Hl_Webb(ff),bins); hl_hr = hl_hr(:,2);
ii = hl_hr > -5 & hl_hr < 10 & fco2w_recalc_hr > -100 & fco2w_recalc_hr < 200;
coef_hl = polyfit(hl_hr(ii),fco2w_recalc_hr(ii),1);
fco2_Hl_line = polyval(coef_hl,-5:1:10);
figure('Position',[300,300,800,800]);
plot(hl_hr,fco2w_recalc_hr,'o','MarkerSize',6); hold on
plot(-5:1:10,fco2_Hl_line,'r-',LineWidth=2); hold off;
axis([-5,10,-100,200]); grid on; grid minor;
ax = gca; ax.FontSize = 12;
xlabel('Latent Heat, W/m2'); ylabel('Fco2 Licor, mM/m^2/d');
title('Licor CO_2 Flux (Webb Corr) vs. Latent Heat Flux','FontSize',18);
ppath = fullfile(path_plots,'co2_flux_licor_vs_hl.png');
print(ppath,'-dpng');

%7
hs_hr = interval_avg(mc_flux.datenum(ff),mc_flux.Hs_2m(ff),bins); hs_hr = hs_hr(:,2);
ii = hs_hr > -10 & hs_hr < 50 & fco2w_recalc_hr > -100 & fco2w_recalc_hr < 200;
coef_hs = polyfit(hs_hr(ii),fco2w_recalc_hr(ii),1);
fco2_Hs_line = polyval(coef_hs,-10:1:50);
figure('Position',[300,300,800,800]);
plot(hs_hr,fco2w_recalc_hr,'o','MarkerSize',6); hold on
plot(-10:1:50,fco2_Hs_line,'r-',LineWidth=2); hold off;
axis([-10,50,-100,200]); grid on; grid minor;
ax = gca; ax.FontSize = 12;
xlabel('Sensible Heat, W/m2'); ylabel('Fco2 Licor, mM/m^2/d');
title('Licor CO_2 Flux (Webb Corr) vs. Sensible Heat Flux','FontSize',18);
ppath = fullfile(path_plots,'co2_flux_licor_vs_hs.png');
print(ppath,'-dpng');

