% b13_houly_fluxes_metcity_merged_dbl_rot.m

%     Loads 13.65 min (2^13 / 8192 point) met and flux vars for the
%     metcity tower.
%{
    Output table mc_flx will have the following columns @ 13 min time base:
    
    1     leg                       cruise leg number
    2     datenum                   matlab datenum, start of interval
    3     yyyy                      year
    4     mm                        month
    5     dd                        day of month
    6     hh                        hour of day, 0-23
    7     seg_begin                 segment begin point (@36000 per hr)
    8     seg_end                   segment end point (@36000 per hr)
    9     lat_tower                 columns 9-100 are met city tower data averaged over flux windows
    10    lon_tower                 ...
    11    heading_tower
    12    heading_mast
    13    lat_mast
    14    lon_mast
    15    zenith_true
    16    zenith_apparent
    17    azimuth
    18    ship_bearing
    19    ship_distance
    20    sr50_dist
    21    temp_2m
    22    temp_6m
    23    temp_10m
    24    temp_mast
    25    dew_point_2m
    26    dew_point_6m
    27    dew_point_10m
    28    dew_point_mast
    29    rh_2m
    30    rh_6m
    31    rh_10m
    32    rh_mast
    33    atmos_pressure_2m
    34    atmos_pressure_mast
    35    brightness_temp_surface
    36    subsurface_heat_flux_A
    37    subsurface_heat_flux_B
    38    snow_depth
    39    mixing_ratio_2m
    40    mixing_ratio_6m
    41    mixing_ratio_10m
    42    mixing_ratio_mast
    43    rhi_2m
    44    rhi_6m
    45    rhi_10m
    46    rhi_mast
    47    vapor_pressure_2m
    48    vapor_pressure_6m
    49    vapor_pressure_10m
    50    vapor_pressure_mast
    51    wspd_vec_mean_2m
    52    wspd_vec_mean_6m
    53    wspd_vec_mean_10m
    54    wspd_vec_mean_mast
    55    wdir_vec_mean_2m
    56    wdir_vec_mean_6m
    57    wdir_vec_mean_10m
    58    wdir_vec_mean_mast
    59    h2o_licor
    60    co2_licor
    61    temp_licor
    62    pressure_licor
    63    wspd_u_mean_2m
    64    wspd_u_mean_6m
    65    wspd_u_mean_10m
    66    wspd_u_mean_mast
    67    wspd_v_mean_2m
    68    wspd_v_mean_6m
    69    wspd_v_mean_10m
    70    wspd_v_mean_mast
    71    wspd_w_mean_2m
    72    wspd_w_mean_6m
    73    wspd_w_mean_10m
    74    wspd_w_mean_mast
    75    temp_acoustic_mean_2m
    76    temp_acoustic_mean_6m
    77    temp_acoustic_mean_10m
    78    temp_acoustic_mean_mast
    79    wspd_u_std_2m
    80    wspd_v_std_2m
    81    wspd_w_std_2m
    82    temp_acoustic_std_2m
    83    wspd_u_std_6m
    84    wspd_v_std_6m
    85    wspd_w_std_6m
    86    temp_acoustic_std_6m
    87    wspd_u_std_10m
    88    wspd_v_std_10m
    89    wspd_w_std_10m
    90    temp_acoustic_std_10m
    91    wspd_u_std_mast
    92    wspd_v_std_mast
    93    wspd_w_std_mast
    94    temp_acoustic_std_mast
    95    skin_temp_surface
    96    down_long_hemisp
    97    down_short_hemisp
    98    up_long_hemisp
    99    up_short_hemisp
    100   net_radiation
    101   wspd_10m                  10m true windspeed, m/s
    102   wdir_10m                  10m true wind direction, deg
    103   Pcav                      Picarro cavity pressure, torr
    104   Tcav                      Picarro cavity temperature, °C
    105   ch4w_ppb                  methane wet mole fraction, ppb
    106   ch4d_ppb                  methane dry mole fraction, ppb
    107   co2w_ppm                  carbon dioxide wet mole fraction, ppm
    108   co2d_ppm                  carbon dioxide wet mole fraction, ppm
    109   h20_pct                   water vapor mole fraction, %
    110   dco2_dt_ppm               co2 rate of change, ppm/hr
    111   dch4_dt_ppb               ch4 rate of change, ppb/hr
    112   sig_co2_ppm               std dev co2 dry mole fraction, ppm
    113   sig_ch4_ppb               std dev ch4 dry mole fraction, ppb
    114   puff_lag_co2              lag between valve and start of Picarro gas puff, points @ 10Hz
    115   co2_ch4_tau               inlet attn time const for co2 & ch4
    116   co2_lag                   'optimal' lag for co2, points @ 10Hz
    117   co2_std_ppm               std dev co2 (despiked), ppm
    118   ch4_lag                   'optimal' lag for ch4, points @ 10Hz
    119   ch4_std_ppb               std dev ch4 (despiked), ppb
    120   max_Corr_co2              maximum correlation coef, w:co2
    121   max_Corr_ch4              maximum correlation coef, w:ch4
    122   puff_stpt                 array index for start of puff segment @ 36000 pnts/hr
    123   puff_endpt                array index for end of puff segment @ 36000 pnts/hr
    124   wu_cov                    w'u' covariance
    125   wu_csp                    w'u' cospectral integral
    126   wv_cov                    w'v' covariance
    127   wv_csp                    w'v' cospectral integral
    128   wTs_cov                   w'Ts' sonic temperature covariance
    129   wTs_csp                   w'Ts' sonic temperature cospectral integral
    130   wT_cov                    w'T' temperature covariance
    131   wT_csp                    w'T' temperature cospectral integral
    132   std_wT                    std dev wT covariance (Finlkelstein et al.)
    133   Ts_noise                  Sonic temperature noise variance
    134   hsens_cov                 sensible heat flux from wT_cov, W/m2
    135   hsens_csp                 sensible heat flux from wT_csp, W/m2
    136   std_hs                    std dev sensible heat flux (Finlkelstein et al.)
    137   wco2_cov_ppm              w'co2' covariance, ppm m/s
    138   wco2_csp_ppm              w'co2' cospectral integral, ppm m/s
    139   wco2_cov_corr_ppm         w'co2' cov attenuation corrected ppm m/s
    140   wco2_csp_corr_ppm         w'co2' csp attenuation corrected ppm m/s
    141   wco2_cov_mM               w'co2' cov attn corrected mMoles/m2/day
    142   wco2_csp_mM               w'co2' csp attn corrected mMoles/m2/day
    143   std_wco2_mM               std dev wco2 covariance (Finlkelstein et al.)
    144   DLwco2_ppm                w'co2' detection limit, ppm m/s
    145   DLwco2_mM                 w'co2' detection limit, mM/m2/d
    146   uco2_mM                   u'co2' streamwise covariance attn corrected, mM/m2/d
    147   vco2_mM                   v'co2' crosswind covariance attn corrected, mM/m2/d
    148   Aco2                      attenuation factor, CO2
    149   wch4_cov_ppb              w'ch4' covariance ppm m/s
    150   wch4_csp_ppb              w'ch4' cospectral integral ppm m/s
    151   wch4_cov_corr_ppb         w'ch4' cov attenuation corrected ppm m/s
    152   wch4_csp_corr_ppb         w'ch4' csp attenuation corrected ppm m/s
    153   wch4_cov_uM               w'ch4' cov attn corrected, uMoles/m2/d
    154   wch4_csp_uM               w'ch4' csp attn corrected, uMoles/m2/d
    155   std_wch4_uM               std dev wch4 covariance (Finlkelstein et al.)
    156   DLwch4_ppb                w'ch4' detection limit, ppm m/s
    157   DLwch4_uM                 w'ch4' detection limit, uMoles/m2/d
    158   uch4_uM                   u'ch4', streamwise covariance attn corrected, uM/m2/d
    159   vch4_uM                   v'ch4' crosswind covariance attn corrected, uM/m2/d
    160   Ach4                      attenuation factor, CH4
    161   std_w                     std dev w vertical wind velocity
    162   std_u                     std dev u streamwise wind velocity
    163   std_v                     std dev v crosswind velocity
    164   std_wdir                  std dev wdir
    165   L_mean                    hourly Obukhov length from mean fluxes @ 2/6/10m levels for intervals passing mc wind filtering, m
    166   zet10                     hourly z/L from mean flux @10m tower level for intervals passing mc wind filtering, unitless
%}

% Computes filtered hourly averages, plots results, and saves hourly mean
% flux table with the following variables:
%{
    1    datenum               matlab datenum
    2    yyyy                  year
    3    mm                    month
    4    dd                    day of month
    5    hh                    hour of day
    6    wspd_hr               vector mean true wind speed, 10m, m/s
    7    wdir_hr               vector mean true wind direction, 10m
    8    Ta_hr                 air temperature, 10m, °C
    9    Pmb_hr                air pressure, 2m, mbar
    10   rh_hr                 relative humidity, 10m, %
    11   sflag_sum_hr          number of segments/hr passing basic wdir filter (0-7)
    12   hs_csp_hr             weighted hourly mean sensible heat flux, W/m2/s
    13   SEhs_csp_hr           std error of weighted mean HS
    14   Fwco2_csp_mM_hr       weighted hourly mean corrected w'co2', mM/m2/d
    15   SEwco2_csp_mM_hr      std error of weighted mean w'co2', mM/m2/d
    16   Fwco2_cov_mM_hr       weighted hourly mean corrected w'co2', mM/m2/d
    17   SEwco2_cov_mM_hr      std error of weighted mean w'co2', mM/m2/d
    18   DLwco2_mM_hr          w'co2' det. lim. (mean DL/sqrt(N)), Finkelstein 2001, mM/m2/d
    19   Fuco2_mM_hr           weighted hourly mean streamwise corrected u'co2', mM/m2/d
    20   sflag_co2_sum_hr      number of segments/hr passing wind and co2 filters (0-7)
    21   co2d_ppm_hr           hourly mean dry co2 mole fraction, ppm
    22   dco2_dt_ppm_hr        mean dco2/dt for intervals passing co2 filter, ppm/hr
    23   Fwch4_csp_uM_hr       weighted hourly mean corrected w'ch4', uM/m2/d
    24   SEwch4_csp_uM_hr      std error of weighted mean w'ch4', uM/m2/d
    25   Fwch4_cov_uM_hr       weighted hourly mean corrected w'ch4', uM/m2/d
    26   SEwch4_cov_uM_hr      std error of weighted mean w'ch4', uM/m2/d
    27   DLwch4_uM_hr          w'ch4' det. lim. (mean DL/sqrt(N)), Finkelstein 2001, uM/m2/d
    28   Fuch4_uM_hr           weighted hourly mean streamwise corrected u'ch4', uM/m2/d
    29   sflag_ch4_sum_hr      number of segments/hr passing wind and ch4 filters (0-7)
    30   ch4d_ppb_hr           hourly mean dry ch4 mole fraction, ppb
    31   dch4_dt_ppb_hr        mean dch4/dt for intervals passing ch4 filter, ppb/hr
    32   rhoair_hr             dry air density, kg/m3
    33   wu_csp_hr             mean w'u' csp integral for intervals passing wind filter, m2/s2
    34   wv_csp_hr             mean w'v' csp integral for intervals passing wind filter, m2/s2
    35   tau_hr                mean streamwise stress for intervals passing wind filter, kg/m/s2
    36   std_w_hr              std dev vertical wind velocity for intervals passing wind filter, m/s
    37   usr_csp_hr            mean ustar for intervals passing wind filter, m/s
    38   wTs_csp_hr            weighted mean w'Tsonic' for intervals passing wind filter, °C m/s
    39   wT_csp_hr             weighted mean w'T' for intervals passing wind filter, °C m/s
    40   tsr_csp_hr            weighted mean tstar for intervals passing wind filter, °C
    41   std_u_hr              std dev u (sqrt(mean(std_u^2)) for intervals passing wind filter, m/s
    42   std_v_hr              std dev v (sqrt(mean(std_v^2)) for intervals passing wind filter, m/s
    43   L_hr                  hourly Obukhov length recomputed from mean fluxes @ 2/6/10m levels for intervals passing mc wind filtering, m
    44   zet10_hr              z/L from hourly mean flux @10m level for intervals passing mc wind filtering, unitless
    45   Rib_hr                'bulk' Richardson number for 10m level w/respect to surface temperature and wind speed
    46   Ri_L1                 'gradient' Richardson number, 2m–6m levels
    47   z_L1                  estimated height for Ri_L1, m
    48   Ri_L2                 as above but for 6m–10m levels
    49   z_L2                  estimated height for Ri_L2, m
    50   Ri_L3                 as above but for 10m to 20m bow tower level
    51   z_L3                  estimated height for Ri_L3, m
    52   Ri_L4                 as above but for 10m to 23-32m mast level
    53   z_L4                  estimated height for Ri_L4, m
%}

% June 2020 BWB

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

plotit = true;     % set true to save plots - script runs faster if plotting is turned off

%% paths
root_data_path = '/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019';
path_prog = fullfile(root_data_path,'_matlab');
save_path = fullfile(root_data_path,'_processed','metcity_flux_v3','flux_merged');
plot_path = fullfile(save_path,'plots','hourly_analysis');
mkdir(plot_path);

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

%% load data

% 13-min flux table 'mc_flx' from b5_compute_fluxes_metcity_8192.m
load(fullfile(save_path,'wins_met_gases_lags_flux_merged_v3.mat'));

% load variable names, units and long names for output hourly mean table
vars = readcell(fullfile(save_path,'flux_vars_hr_V3.xlsx'));

% hourly Ri and flux profile arrays from h_analyze_profiles.m (bulk model
% and winds script)
load(fullfile(root_data_path,'_processed','metcity_tower_level2_v3','Ri_profile_4levels_hourly.mat'));

%% compute hourly means

% bins for 1-hour averaging
dm_bins = floor(mc_flx.datenum(1)):1/24:ceil(mc_flx.datenum(end)); dm_bins = dm_bins(:);

% basic wind filter for turbulence parameters: wdir outside of +/- 10° of
% bearing to ship and when U10 > 4 m/s
test = abs(mc_flx.ship_bearing - mc_flx.wdir_vec_mean_10m);
sflag = test>10 & test<350 & mc_flx.wspd_vec_mean_10m>=4;

% set trace gas selection filters for non-stationary conditions
% limits are the 99 percentile of the Weibull distribution for each
% parameter.  See plots from stationarity_distributioins.m
dco2_dt_lim = 1.38;
uco2_lim = 40.50;
sig_co2_lim = 0.27; % ppm
dch4_dt_lim = 7.15;
uch4_lim = 197;
sig_ch4_lim = 0.74; % ppb
sco2flag = abs(mc_flx.dco2_dt_ppm)<=dco2_dt_lim & abs(mc_flx.uco2_mM)<=uco2_lim & mc_flx.sig_co2_ppm<=sig_co2_lim;
sch4flag = abs(mc_flx.dch4_dt_ppb)<=dch4_dt_lim & abs(mc_flx.uch4_uM)<=uch4_lim & mc_flx.sig_ch4_ppb<=sig_ch4_lim;

% selection filters for each gas
cflag = sflag & sco2flag;   % filter co2, in-sector, stationary
mflag = sflag & sch4flag ;  % filter ch4, in-sector, stationary

% filtered hourly mean stress and wind stats
wu_cov_hr = interval_avg(mc_flx.datenum(sflag),mc_flx.wu_cov(sflag),dm_bins); wu_cov_hr = wu_cov_hr(:,2);
wv_cov_hr = interval_avg(mc_flx.datenum(sflag),mc_flx.wv_cov(sflag),dm_bins); wv_cov_hr = wv_cov_hr(:,2);
wu_csp_hr = interval_avg(mc_flx.datenum(sflag),mc_flx.wu_csp(sflag),dm_bins); wu_csp_hr = wu_csp_hr(:,2);
wv_csp_hr = interval_avg(mc_flx.datenum(sflag),mc_flx.wv_csp(sflag),dm_bins); wv_csp_hr = wv_csp_hr(:,2);
std_w_hr = interval_avg(mc_flx.datenum(sflag),mc_flx.std_w(sflag).^2,dm_bins); std_w_hr = std_w_hr(:,2);
std_w_hr = sqrt(std_w_hr);
std_u_hr = interval_avg(mc_flx.datenum(sflag),mc_flx.std_w(sflag).^2,dm_bins); std_u_hr = std_u_hr(:,2);
std_u_hr = sqrt(std_u_hr);
std_v_hr = interval_avg(mc_flx.datenum(sflag),mc_flx.std_w(sflag).^2,dm_bins); std_v_hr = std_v_hr(:,2);
std_v_hr = sqrt(std_v_hr);

% number of good segments per hour
sflag_sum_hr = interval_sum(mc_flx.datenum,sflag,dm_bins);
dt_1hr = datetime(sflag_sum_hr(:,1)+0.5/24,'convertfrom','datenum'); % midpoint timestamp
sflag_sum_hr = sflag_sum_hr(:,2); % # good segments per hour

sflag_co2_sum_hr = interval_sum(mc_flx.datenum,cflag,dm_bins);
sflag_co2_sum_hr = sflag_co2_sum_hr(:,2); % # good segments per hour

sflag_ch4_sum_hr = interval_sum(mc_flx.datenum,mflag,dm_bins);
sflag_ch4_sum_hr = sflag_ch4_sum_hr(:,2); % # good segments per hour

% hourly mean met variables
Pmb_hr = interval_avg(mc_flx.datenum,mc_flx.atmos_pressure_2m,dm_bins); Pmb_hr = Pmb_hr(:,2);
Ta_hr = interval_avg(mc_flx.datenum,mc_flx.temp_10m,dm_bins); Ta_hr = Ta_hr(:,2);
[~,wdir_hr,wspd_hr] = interval_avg_vectors(mc_flx.datenum,mc_flx.wdir_vec_mean_10m,mc_flx.wspd_vec_mean_10m,dm_bins);
rh_hr = interval_avg(mc_flx.datenum,mc_flx.rh_10m,dm_bins); rh_hr = rh_hr(:,2);
Tsfc_hr = interval_avg(mc_flx.datenum,mc_flx.skin_temp_surface,dm_bins); Tsfc_hr = Tsfc_hr(:,2);
% L and zeta are constant for each hour so mean will just return the hourly values...
L_hr = interval_avg(mc_flx.datenum,mc_flx.L_mean,dm_bins); L_hr = L_hr(:,2);
zet10_hr = interval_avg(mc_flx.datenum,mc_flx.zet10,dm_bins); zet10_hr = zet10_hr(:,2);

% compute hourly bulk Richardson number between 10m and surface
dU = wspd_hr;                        % delta U is ~= U10
Tv = Ta_hr + 273.15 + 0.0098*10;     % 10m virtual T ~= 10m Theta
dz = 10;                             % delta z is 10m to surface
dT = Tv - (Tsfc_hr + 273.15);        % delta Theta, 10m to sfc
Rib_hr = (9.82*dz.*dT)./(Tv.*dU.^2); % Ri = (g*dz*dT)/(Tv*dU)

% copy over Ri, zRI, L and zeta from profile analysis
Ri_L1 = R_prof_3_5_hr(:,2);
Ri_L2 = R_prof_7_5_hr(:,2);
Ri_L3 = R_prof_14_2_hr(:,2);
Ri_L4 = R_prof_16_hr(:,2);
z_L1 = R_prof_3_5_hr(:,3);
z_L2 = R_prof_7_5_hr(:,3);
z_L3 = R_prof_14_2_hr(:,3);
z_L4 = R_prof_16_hr(:,3);

% compute hourly weighted means and std error for fluxes
% from 13 min fluxes that pass selection filters.  Weighting uses
% 13min flux variances estimated from Fink4elstein et al. method
% to compute standard error of the hourly mean.
[hs_cov_hr,SEhs_cov_hr,~] = interval_weighted_avg(mc_flx.datenum(sflag),mc_flx.hsens_cov(sflag),...
    mc_flx.std_hs(sflag), dm_bins);
[hs_csp_hr,SEhs_csp_hr,~] = interval_weighted_avg(mc_flx.datenum(sflag),mc_flx.hsens_csp(sflag),...
    mc_flx.std_hs(sflag), dm_bins);
[wTs_cov_hr,SEwTs_cov_hr,~] = interval_weighted_avg(mc_flx.datenum(sflag),mc_flx.wTs_cov(sflag),...
    mc_flx.std_wT(sflag), dm_bins);
[wTs_csp_hr,SEwTs_csp_hr,~] = interval_weighted_avg(mc_flx.datenum(sflag),mc_flx.wTs_csp(sflag),...
    mc_flx.std_wT(sflag), dm_bins);
[wT_cov_hr,SEwT_cov_hr,~] = interval_weighted_avg(mc_flx.datenum(sflag),mc_flx.wT_cov(sflag),...
    mc_flx.std_wT(sflag), dm_bins);
[wT_csp_hr,SEwT_csp_hr,~] = interval_weighted_avg(mc_flx.datenum(sflag),mc_flx.wT_csp(sflag),...
    mc_flx.std_wT(sflag), dm_bins);
[Fwco2_csp_mM_hr,SEwco2_csp_mM_hr,~] = interval_weighted_avg(mc_flx.datenum(cflag),mc_flx.wco2_csp_corr_mM(cflag),...
    mc_flx.std_wco2_mM(cflag), dm_bins);
[Fwco2_cov_mM_hr,SEwco2_cov_mM_hr,~] = interval_weighted_avg(mc_flx.datenum(cflag),mc_flx.wco2_cov_corr_mM(cflag),...
    mc_flx.std_wco2_mM(cflag), dm_bins);
[Fwch4_csp_uM_hr,SEwch4_csp_uM_hr,~] = interval_weighted_avg(mc_flx.datenum(mflag),mc_flx.wch4_csp_corr_uM(mflag),...
    mc_flx.std_wch4_uM(mflag), dm_bins);
[Fwch4_cov_uM_hr,SEwch4_cov_uM_hr,~] = interval_weighted_avg(mc_flx.datenum(mflag),mc_flx.wch4_cov_corr_uM(mflag),...
    mc_flx.std_wch4_uM(mflag), dm_bins);

% hourly detection limits: simple average 13-min DL divided by root N,
% where N is determined by the number of valid segments per hour
DLwco2_mM_hr = interval_avg(mc_flx.datenum(cflag),mc_flx.DLwco2_mM(cflag),dm_bins); DLwco2_mM_hr = DLwco2_mM_hr(:,2);
fff = isfinite(sflag_co2_sum_hr) & sflag_co2_sum_hr>0;
DLwco2_mM_hr(fff) = DLwco2_mM_hr(fff)./sqrt(sflag_co2_sum_hr(fff));

DLwch4_uM_hr = interval_avg(mc_flx.datenum(mflag),mc_flx.DLwch4_uM(mflag),dm_bins); DLwch4_uM_hr = DLwch4_uM_hr(:,2);
fff = isfinite(sflag_ch4_sum_hr) & sflag_ch4_sum_hr>0;
DLwch4_uM_hr(fff) = DLwch4_uM_hr(fff)./sqrt(sflag_ch4_sum_hr(fff));

% mean concentrations and dC/dt for each gas
co2d_ppm_hr = interval_avg(mc_flx.datenum(cflag),mc_flx.co2d_ppm(cflag),dm_bins); co2d_ppm_hr = co2d_ppm_hr(:,2);
dco2_dt_ppm_hr = interval_avg(mc_flx.datenum(cflag),mc_flx.dco2_dt_ppm(cflag),dm_bins); dco2_dt_ppm_hr = dco2_dt_ppm_hr(:,2);
ch4d_ppm_hr = interval_avg(mc_flx.datenum(mflag),mc_flx.ch4d_ppb(mflag),dm_bins); ch4d_ppm_hr = ch4d_ppm_hr(:,2);
dch4_dt_ppb_hr = interval_avg(mc_flx.datenum(mflag),mc_flx.dch4_dt_ppb(mflag),dm_bins); dch4_dt_ppb_hr = dch4_dt_ppb_hr(:,2);

% compute hourly turbulence scale parameters
% these have sflag applied, so +/-20° sector toward ship is excluded
usr_cov_hr = NaN(size(wu_cov_hr)); usr_csp_hr = usr_cov_hr;
ii = isfinite(wu_cov_hr);
% ustar from filtered streamwise wu
usr_cov_hr(ii) = -sign(wu_cov_hr(ii)).*sqrt(abs(wu_cov_hr(ii)));
usr_csp_hr(ii) = -sign(wu_csp_hr(ii)).*sqrt(abs(wu_csp_hr(ii)));
% usr_cov = (wu_hr.^2 + wv_hr.^2).^(1/4);
% eh2o = 6.112.*exp(17.67.*Ta_hr./(Ta_hr+243.5)).*rh_hr/100;
% wT2_hr = wTs_hr./(1 + 0.32*eh2o./Pmb_hr);
tsr_cov_hr = -wT_cov_hr./usr_cov_hr;
tsr_csp_hr = -wT_csp_hr./usr_csp_hr;

%% load hourly mean flux file from POP group analysis of tower data
% same timestamp as mc_flx_hr table
dfl = '/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019/_processed/metcity_tower_level2_v3/metcity_1hr_corA10_all.mat';
load (dfl); mc_tower_hr = mc_flux_hr; clear mc_flux_hr; % renaming tower met/flux table...
dt_pop = datetime(mc_tower_hr.datenum,'convertfrom','datenum');
% basic wind filter for stress
test = abs(mc_tower_hr.ship_bearing - mc_tower_hr.wdir_vec_mean_10m);
sflag_hr_pop = test>10 & test<350;  % excluding when wdir is within +/- 10° of ship bearing.

%% plots
plotit = 1; % setting plotit to 1 in case it was set to zero at start of script

pname = 'usr_cov_vs_csp';
plot_one2one(pname,usr_cov_hr,usr_csp_hr);
title('Hourly ustar, cov vs. csp','fontsize',18);
xlabel('ustar by covariance, m/s','fontsize',14);
ylabel('ustar by cospec integral, m/s','fontsize',14);
ppath = fullfile(plot_path,'ustar_cov_vs_csp_hourly_metcity.png');
print(ppath,'-dpng');

pname = 'tsr_cov_vs_csp';
plot_one2one(pname,tsr_cov_hr,tsr_csp_hr);
ylim([-0.2,0.12]);
title('Hourly tstar, cov vs. csp','fontsize',18);
xlabel('tstar by covariance, °C','fontsize',14);
ylabel('tstar by cospec integral, °C','fontsize',14);
ppath = fullfile(plot_path,'tstar_cov_vs_csp_hourly_metcity.png');
print(ppath,'-dpng');

pname = 'Fco2_cov_vs_csp';
plot_one2one(pname,Fwco2_cov_mM_hr,Fwco2_csp_mM_hr);
title('Hourly w''co2'', cov vs. csp','fontsize',18);
xlabel('w''co2'' by covariance, mM/m2/d','fontsize',14);
ylabel('w''co2'' by cospec integral, mM/m2/d','fontsize',14);
ppath = fullfile(plot_path,'wco2_cov_vs_csp_hourly_metcity.png');
print(ppath,'-dpng');

pname = 'Fch4_cov_vs_csp';
plot_one2one(pname,Fwch4_cov_uM_hr,Fwch4_csp_uM_hr);
title('Hourly w''ch4'', cov vs. csp','fontsize',18);
xlabel('w''ch4'' by covariance, uM/m2/d','fontsize',14);
ylabel('w''ch4'' by cospec integral, uM/m2/d','fontsize',14);
ppath = fullfile(plot_path,'wch4_cov_vs_csp_hourly_metcity.png');
print(ppath,'-dpng');

% plot a few inlet attenuation factors for trace gases
dt = datetime(mc_flx.datenum,'convertfrom','datenum');
figure('position',[500,100,1000,700],'Name','attn_factor_timeseries');
plot(dt,mc_flx.Aco2,'.','markersize',10); grid on; grid minor;
xlabel('Date / Time'); ylabel('CO_2/CH_4 Inlet Attn Factors');
title('MOSAiC MetCity Fluxes, Fractional Inlet Attenuation Factors','fontsize',16);
ppath = fullfile(plot_path,'inlet_attn_factors_metcity.png');
print(ppath,'-dpng');
ppath = fullfile(plot_path,'inlet_attn_factors_metcity.fig');
savefig(ppath);

% sensible heat flux
figure('position',[400,400,650,400],'Name','sensible_heat_flux');
plot([dt_1hr(1),dt_1hr(end)],[0,0],'k-','linewidth',1); hold on; grid on; grid minor;
plot(dt_1hr,hs_csp_hr,'.','markersize',15); hold off;
xlim([dt_1hr(1),dt_1hr(end)]); %ylim([-40,40]);
xlabel('Time / Date'); ylabel('Sensible Heat Flux, W/m^2');
title('MOSAiC, MetCity Hourly Mean Sensible Heat Flux','fontsize',16);
ppath = fullfile(plot_path,'hsFlux_metcity.png');
print(ppath,'-dpng');
ppath = fullfile(plot_path,'hsFlux_metcity.fig');
savefig(ppath);

% sensible heat flux distribution
figure('position',[400,400,650,400],'Name','sensible_heat_flux_dist');
histogram(hs_csp_hr,[-30:1:30]); grid on; grid minor;
ylabel('count','fontsize',12); xlabel('Sensible Heat Flux, W/m^2','fontsize',12);
title('MOSAiC, MetCity Hourly Sensible Heat Flux Distribution','fontsize',18);
ppath = fullfile(plot_path,'hsFlux_dist_metcity.png');
print(ppath,'-dpng');

% CO2 flux
Fuco2_mM_hr = interval_avg(mc_flx.datenum(cflag),mc_flx.uco2_mM(cflag),dm_bins); Fuco2_mM_hr = Fuco2_mM_hr(:,2);
cc = sflag_co2_sum_hr>2 & abs(Fwco2_csp_mM_hr)>DLwco2_mM_hr;
ccc = sflag_co2_sum_hr>2 & abs(Fwco2_csp_mM_hr)<=DLwco2_mM_hr;
figure('position',[100,750,600,800],'Name','co2_flux_timeseries');
subplot(2,1,1)
    %     yyaxis left;
    plot(dt_1hr(ccc),Fwco2_csp_mM_hr(ccc),'b.','markersize',15); hold on; grid on; grid minor;
    plot(dt_1hr(cc),Fwco2_csp_mM_hr(cc),'r.','markersize',15);
    plot([dt_1hr(1),dt_1hr(end)],[0,0],'k-','linewidth',1); hold off;
    ylabel('CO_2 flux, w''c'', mM/m^2/d');
    xlim([dt_1hr(1),dt_1hr(end)]); ylim([-15,15]);
    %     yyaxis right; plot(dt_co2,UFco2_mean(:,2),'x','markersize',6);
    %     ylabel('advective CO_2 flux, u''c'',mM/m^2/d');
    legend('<DL','>DL','location','north'); xlabel('Time / Date');
    title('MOSAiC, MetCity Hourly Mean Carbon Dioxide Flux','fontsize',16);
subplot(2,1,2)
    plot(dt_1hr,DLwco2_mM_hr,'.','markersize',15); grid minor;
    xlabel('Time / Date'); ylabel('Hourly Flux Det Lim, mM/m^2/d');
    xlim([dt_1hr(1),dt_1hr(end)]); ylim([0,6]);
ppath = fullfile(plot_path,'co2Flux_metcity.png');
print(ppath,'-dpng');
ppath = fullfile(plot_path,'co2Flux_metcity.fig');
savefig(ppath)

% CH4 flux
Fuch4_uM_hr = interval_avg(mc_flx.datenum(sco2flag),mc_flx.uch4_uM(sco2flag),dm_bins); Fuch4_uM_hr = Fuch4_uM_hr(:,2);
mm = sflag_ch4_sum_hr>2 & abs(Fwch4_csp_uM_hr)>DLwch4_uM_hr;
mmm = sflag_ch4_sum_hr>2 & abs(Fwch4_csp_uM_hr)<=DLwch4_uM_hr;
figure('position',[100,750,600,800],'Name','ch4_flux_timeseries');
subplot(2,1,1)
    %     yyaxis left;
    plot(dt_1hr(mmm),Fwch4_csp_uM_hr(mmm),'b.','markersize',15); grid on; grid minor; hold on;
    plot(dt_1hr(mm),Fwch4_csp_uM_hr(mm),'r.','markersize',15); grid minor;
    plot([dt_1hr(1),dt_1hr(end)],[0,0],'k-','linewidth',1);
    ylabel('CH_4 flux, w''c'', uM/m^2/d');
    xlim([dt_1hr(1),dt_1hr(end)]); ylim([-100,100]);
    %     yyaxis right; plot(dt_ch4,UFch4_mean(:,2),'x','markersize',6);
    %     ylabel('advective CH_4 flux, u''c'', uM/m^2/d');
    legend('<DL','>DL','location','north'); xlabel('Time / Date');
    title('MOSAiC, MetCity Hourly Mean Methane Flux','fontsize',16);
subplot(2,1,2)
    plot(dt_1hr,DLwch4_uM_hr,'.','markersize',15); grid minor;
    xlabel('Time / Date'); ylabel('Hourly Flux Det Lim, mM/m^2/d');
    xlim([dt_1hr(1),dt_1hr(end)]);
ppath = fullfile(plot_path,'ch4Flux_metcity.png');
print(ppath,'-dpng');
ppath = fullfile(plot_path,'ch4Flux_metcity.fig');
savefig(ppath);

% u'c'/v'c' vs. dC/dt
figure('position',[100,750,600,800],'Name','stationarity_params');
subplot(2,1,1)
    plot([-6,6],[0,0],'k-',[0,0],[-200,200],'k-','linewidth',1); hold on;
    h1 = plot(mc_flx.dco2_dt_ppm(sflag),mc_flx.uco2_mM(sflag),'r.',mc_flx.dco2_dt_ppm(sflag),mc_flx.vco2_mM(sflag),'b.','markersize',5);
    xlabel('dCO_2/dt, ppm/hr'); ylabel('CO_2 adv. turb. flux, mM/m^2/d'); grid;
    xlim([-6,6]); ylim([-100,100]); yticks(-100:25:100);
    rectangle('Position',[-dco2_dt_lim,-uco2_lim,dco2_dt_lim*2,uco2_lim*2],'linewidth',2,'linestyle','-');
    legend(h1,'u''co_2''','v''co_2''','location','northwest');
    title('MOSAiC, MetCity CO_2/CH_4 Stationarity Test','fontsize',16);
subplot(2,1,2)
    plot([-30,30],[0,0],'k-',[0,0],[-1000,1000],'k-','linewidth',1); hold on;
    h2 = plot(mc_flx.dch4_dt_ppb(sflag),mc_flx.uch4_uM(sflag),'r.',mc_flx.dch4_dt_ppb(sflag),mc_flx.vch4_uM(sflag),'b.','markersize',5);
    xlabel('dCH_4/dt, ppb/hr'); ylabel('CH_4 adv. turb. flux, uM/m^2/d'); grid;
    xlim([-30,30]); ylim([-500,500]); yticks(-500:100:500);
    rectangle('Position',[-dch4_dt_lim,-uch4_lim,dch4_dt_lim*2,uch4_lim*2],'linewidth',2,'linestyle','-');
    legend(h2,'u''ch_4''','v''ch_4''','location','northwest');
ppath = fullfile(plot_path,'horiz_flux_v_dcdt_metcity.png');
print(ppath,'-dpng');
ppath = fullfile(plot_path,'horiz_flux_v_dcdt_metcity.fig');
savefig(ppath);

% CO2/CH4 flux detection limit distribution
fff = sflag_co2_sum_hr>2;
ggg = sflag_ch4_sum_hr>2;
figure('position',[500,100,800,900],'Name','co2_ch4_fluxDL_dist');
subplot(2,1,1);
    histogram(DLwco2_mM_hr(fff),[0:0.1:4]); grid on; grid minor;
    xlabel('CO_2 Flux Det.Lim., mM/m^2/d','fontsize',12);
    ylabel('N_{tot}','fontsize',12);
    title('MOSAiC, MetCity CO_2/CH_4 1\sigma Flux Detection Limit Distributions','fontsize',18);
    subtitle('N > 2 good segments / hr','fontsize',12);
subplot(2,1,2);
    histogram(DLwch4_uM_hr(ggg),[0:1:30]); grid on; grid minor;
    xlabel('CH_4 Flux Det.Lim., \muM/m^2/d','fontsize',12);
    ylabel('N_{tot}','fontsize',12);
ppath = fullfile(plot_path,'co2_ch4_fluxDL_dist_metcity.png');
print(ppath,'-dpng');

% CO2/CH4 flux distribution
figure('position',[500,100,800,900],'Name','co2_ch4_flux_dist');
subplot(2,1,1);
    histogram(Fwco2_csp_mM_hr(fff),[-15:0.2:15]); grid on; grid minor;
    xlabel('CO_2 Flux, mM/m^2/d','fontsize',12);
    ylabel('N_{tot}','fontsize',12); ylim([0,250]);
    title('MOSAiC, MetCity CO_2/CH_4 Flux Distributions','fontsize',18);
    subtitle('N > 2 good segments / hr','fontsize',12);
subplot(2,1,2);
    histogram(Fwch4_csp_uM_hr(ggg),[-80:2.5:80]); grid on; grid minor;
    xlabel('CH_4 Flux, \muM/m^2/d','fontsize',12);
    ylabel('N_{tot}','fontsize',12); ylim([0,500]);
ppath = fullfile(plot_path,'o2_ch4_flux_dist_metcity.png');
print(ppath,'-dpng');

% CO2/CH4 flux distribution - above DL measurements
fff = sflag_co2_sum_hr>2 & abs(Fwco2_csp_mM_hr)>DLwco2_mM_hr;
ggg = sflag_ch4_sum_hr>2 & abs(Fwch4_csp_uM_hr)>DLwch4_uM_hr;
figure('position',[500,100,800,900],'Name','co2_ch4_aboveDL_dist');
subplot(2,1,1);
    histogram(Fwco2_csp_mM_hr(fff),[-15:0.2:15]); grid on; grid minor;
    xlabel('CO_2 Flux, mM/m^2/d','fontsize',12);
    ylabel('N_{tot}','fontsize',12); ylim([0,100]);
    title('MOSAiC, MetCity CO_2/CH_4 Flux Distributions','fontsize',18);
    subtitle('N > 2 good segments / hr & above 1\sigma DL','fontsize',12);
subplot(2,1,2);
    histogram(Fwch4_csp_uM_hr(ggg),[-100:2.5:100]); grid on; grid minor;
    xlabel('CH_4 Flux, \muM/m^2/d','fontsize',12);
    ylabel('N_{tot}','fontsize',12); ylim([0,200]);
ppath = fullfile(plot_path,'co2_ch4_flux_dist_aboveDL_metcity.png');
print(ppath,'-dpng');

% # good segments per hr, frequency distribution
figure('position',[500,100,800,900],'Name','Ngood_dist');
subplot(2,1,1);
    histogram(sflag_co2_sum_hr); grid on; grid minor;
    xlabel('N_{segs} F_{co2} per hr','fontsize',12);
    ylabel('N_{tot}','fontsize',12);
    annotation('textbox',[0.3,0.8,0.1,0.1],'String',['CO_2 total (N>2) = ',sprintf('%4i',sum(sflag_co2_sum_hr>2))],...
        'fontsize',12);
    title('MOSAiC, MetCity CO_2/CH_4/O_3 Flux Segs per hr, Distribution','fontsize',18);
subplot(2,1,2);
    histogram(sflag_ch4_sum_hr); grid on; grid minor;
    xlabel('N_{segs} F_{ch4} per hr','fontsize',12);
    ylabel('N_{tot}','fontsize',12);
annotation('textbox',[0.3,0.3,0.1,0.1],'String',['CH_4 total (N>2) = ',sprintf('%4i',sum(sflag_ch4_sum_hr>2))],...
    'fontsize',12);
ppath = fullfile(plot_path,'co2_ch4_Ngood_per_hr_dist_metcity.png');
print(ppath,'-dpng');

% stream-wise stress vs. wspd
rhoair_hr = rhoa(Ta_hr,Pmb_hr,rh_hr);
tau_hr = -rhoair_hr.*wu_csp_hr;
[rw_bin, tau_bin, med_tau_bin, se_bin, n_bin] = binave3(wspd_hr,tau_hr,1,0,20);
figure('position',[400,400,650,400],'Name','stress_v_u10_hr');
plot([0,20],[0,0],'k-','linewidth',1.5); grid on; grid minor; hold on;
plot(wspd_hr,tau_hr,'r.','markersize',4);
errorbar(rw_bin,tau_bin,se_bin,'b-o','linewidth',2); ylim([-0.1,0.4]);
ylabel('Stress, N/m^2','fontsize',12); xlabel('U_{10} m/s','fontsize',12);
title('MOSAiC MetCity, Wind Stress vs. Wind Speed','fontsize',18);
ppath = fullfile(plot_path,'ustress_v_wspd_hourly_metcity.png');
print(ppath,'-dpng');

rhoair  = rhoa(mc_flx.temp_10m,mc_flx.atmos_pressure_2m,mc_flx.rh_10m);
tau = NaN(size(mc_flx.wu_csp));
% tflag = mc_flx.rwdir>-90 & mc_flx.rwdir<90 & mc_flx.hed_std<20 & mc_flx.rwspd>1 &...
%         mc_flx.son_tilt>-5 & mc_flx.son_tilt<10;
tau(sflag) = -mc_flx.wu_csp(sflag).*rhoair(sflag);
figure('position',[400,400,650,400],'Name','stress_v_u10_13min');
plot([0,20],[0,0],'k-','linewidth',1.5); grid on; grid minor; hold on;
plot(mc_flx.wspd_vec_mean_10m,tau,'r.','markersize',4); ylim([-0.5,0.6]);
ylabel('Stress, N/m^2','fontsize',12); xlabel('U_{10} m/s','fontsize',12);
title('MOSAiC MetCity, Wind Stress vs. Wind Speed','fontsize',18);
ppath = fullfile(plot_path,'ustress_v_wspd_raw_metcity.png');
print(ppath,'-dpng');

% cross-stream stress vs. wspd
taux_hr = -rhoair_hr.*wv_csp_hr;
[rwx_bin, taux_bin, med_taux_bin, setx_bin, nx_bin] = binave3(wspd_hr,taux_hr,1,0,20);
figure('position',[400,400,650,400],'Name','xtream_stress_v_u10_hr');
plot([0,20],[0,0],'k-','linewidth',1.5); grid on; grid minor; hold on;
plot(wspd_hr,taux_hr,'r.','markersize',4);  ylim([-0.2,1.4]);
errorbar(rwx_bin,taux_bin,setx_bin,'b-o','linewidth',2); ylim([-0.25,0.25]);
ylabel('X-stream Stress, N/m^2','fontsize',12); xlabel('U_{10} m/s','fontsize',12);
title('MOSAiC MetCity, Cross-stream Wind Stress vs. Wind Speed','fontsize',18);
ppath = fullfile(plot_path,'vstress_v_wspd_hourly_metcity.png');
print(ppath,'-dpng');

% sigma-w vs wdir
sigwflag = mc_flx.wspd_vec_mean_10m>1;
figure('position',[400,400,650,400],'Name','sigw_v_wdir');
plot(mc_flx.wdir_vec_mean_10m(sigwflag),mc_flx.std_w(sigwflag),'.','markersize',4); grid on; grid minor;
xlim([0,360]); ylim([0,3]);
xlabel('Wind Dir','fontsize',12); ylabel('\sigma_w  m/s','fontsize',12);
title('MOSAiC, MetCity Sigma-W vs Wind Direction','fontsize',18);
ppath = fullfile(plot_path,'sigw_v_wdir_metcity.png');
print(ppath,'-dpng');

% hourly sigma-w/usr(z/L) vs wdir
usr = -sign(mc_flx.wu_csp).*sqrt(abs(mc_flx.wu_csp));
Nsigw = NaN(size(mc_flx.std_w));
ii = mc_flx.zet10<0; % for the unstable cases
Nsigw(ii) = mc_flx.std_w(ii)./(usr(ii).*(1 - 3*mc_flx.zet10(ii)).^(1/3));
ii = mc_flx.zet10>0; % for the stable cases
Nsigw(ii) = mc_flx.std_w(ii)./(usr(ii).*(1 + 0.2*mc_flx.zet10(ii)));

figure('position',[400,400,650,400],'Name','Norm_sigw_v_wdir');
plot(mc_flx.wdir_vec_mean_10m(sigwflag),Nsigw(sigwflag),'.','markersize',4); hold on;
plot([0,360],[1.25,1.25],'k-','linewidth',1); hold off;
xlim([0,360]); ylim([0,4]); grid on; grid minor;
xlabel('Wind Dir','fontsize',12); ylabel('\sigma_w / (u_* \psi(z/l)/1.25)','fontsize',12);
title('MOSAiC, MetCity Normalized Sigma-W vs Wind Dir','fontsize',18);
ppath = fullfile(plot_path,'Norm_sigw_v_wdir_metcity.png');
print(ppath,'-dpng');

%% save 1-hr mean values

dm_hr = dm_bins(1:end-1);
dm_vec = datevec(dm_hr);

% add variables to table using var name as defined above
mc_flx_hr = array2table([dm_hr,dm_vec(:,1:4)],'VariableNames',{'datenum','yyyy','mm','dd','hh'});
mc_flx_hr = addvars(mc_flx_hr,wspd_hr,wdir_hr,Ta_hr,Pmb_hr,rh_hr,sflag_sum_hr,hs_csp_hr,SEhs_csp_hr,...
    Fwco2_csp_mM_hr,SEwco2_csp_mM_hr,Fwco2_cov_mM_hr,SEwco2_cov_mM_hr,DLwco2_mM_hr,Fuco2_mM_hr,...
    sflag_co2_sum_hr,co2d_ppm_hr,dco2_dt_ppm_hr,Fwch4_cov_uM_hr,SEwch4_cov_uM_hr,Fwch4_csp_uM_hr,...
    SEwch4_csp_uM_hr,DLwch4_uM_hr,Fuch4_uM_hr,sflag_ch4_sum_hr,ch4d_ppm_hr,dch4_dt_ppb_hr,rhoair_hr,...
    wu_csp_hr,wv_csp_hr,std_w_hr,usr_csp_hr,wTs_csp_hr,wT_csp_hr,tsr_csp_hr,std_u_hr,std_v_hr,Tsfc_hr,...
    L_hr,zet10_hr,Rib_hr,Ri_L1,z_L1,Ri_L2,z_L2,Ri_L3,z_L3,Ri_L4,z_L4);

% rename vars using names from Excel spreadsheet, just in case some names
% have been updated.  Also add units and long names.
mc_flx_hr.Properties.VariableNames = vars(:,1);
mc_flx_hr.Properties.VariableUnits = vars(:,2);
mc_flx_hr.Properties.VariableDescriptions = vars(:,3);

sfl = fullfile(save_path,'MOSAiC_metcity_flux_co2_ch4_hr.mat');
save(sfl,'mc_flx_hr');

%% other plots

% compare usr with POP group's hourly value
[~,ia,ib] = intersect(mc_flx_hr.datenum,mc_tower_hr.datenum);
plot_one2one('ustar_blo_v_pop',mc_flx_hr.usr_csp_hr(ia),mc_tower_hr.ustar_10m(ib));
ylim([-0.2,0.8]);
xlabel('ustar (BLO)','fontsize',14); ylabel('ustar (POP)','fontsize',14);
title('MOSAiC MetCity, Hourly ustar comparison','fontsize',18);
ppath = fullfile(plot_path,'ustar_blo_v_pop_metcity.png');
print(ppath,'-dpng');

% compare usr with corA10
[~,ia,ib] = intersect(mc_flx_hr.datenum,corA10hr_10m.dnum);
plot_one2one('ustar_blo_v_bulk',mc_flx_hr.usr_csp_hr(ia),corA10hr_10m.usr(ib));
xlabel('ustar (BLO)','fontsize',14); ylabel('ustar (bulk, 10m)','fontsize',14);
title('MOSAiC MetCity, Hourly ustar comparison','fontsize',18);
ppath = fullfile(plot_path,'ustar_blo_v_bulk_metcity.png');
print(ppath,'-dpng');

% timeseries of POP and BLO ustar
Nblo = sum(isfinite(mc_flx_hr.usr_csp_hr));
Npop = sum(isfinite(mc_tower_hr.ustar_10m(sflag_hr_pop)));
figure('position',[500,100,1000,700],'Name','blo_pop_ustar_timeseries');
plot(dt_1hr,mc_flx_hr.usr_csp_hr,'.',dt_pop(sflag_hr_pop),mc_tower_hr.ustar_10m(sflag_hr_pop),'.','markersize',10);
grid on; grid minor; ylabel('ustar, m/s');
legend(['BLO, N = ',sprintf('%5i',Nblo)],['POP, N = ',sprintf('%5i',Npop)],'location','northwest')';
title('MOSAiC MetCity, ustar timeseries','fontsize',18);
ppath = fullfile(plot_path,'ustar_blo_pop_timeseries_metcity.png');
print(ppath,'-dpng');

% CO2 / CH4 hr flux timeseries - Legs 4/5
figure('position',[400,400,800,800],'Name','L4_L5_co2_ch4_flux_timeseries');
subplot(2,1,1)
    plot(dt_1hr(sflag_co2_sum_hr>2),Fwco2_csp_mM_hr(sflag_co2_sum_hr>2),'b.','markersize',15); hold on; grid on; grid minor;
    plot(dt_1hr(cc),Fwco2_csp_mM_hr(cc),'r.','markersize',15);
    plot([dt_1hr(1),dt_1hr(end)],[0,0],'k-','linewidth',1); hold off;
    ylabel('CO_2 flux, w''c'', mM/m^2/d','fontsize',14); xlabel('Time / Date','fontsize',14);
    % ylim([-15,15]);
    % xlim([dt_3hr(1),dt_3hr(end)]);
    xlim([datetime(2020,6,15),datetime(2020,9,30)]);
    title('MetCity CO_2 / CH_4 flux, hr mean, Legs 4-5','fontsize',20);
subplot(2,1,2)
    plot(dt_1hr(sflag_ch4_sum_hr>2),Fwch4_csp_uM_hr(sflag_ch4_sum_hr>2),'b.','markersize',15); grid on; grid minor; hold on;
    plot(dt_1hr(mm),Fwch4_csp_uM_hr(mm),'r.','markersize',15);
    plot([dt_1hr(1),dt_1hr(end)],[0,0],'k-','linewidth',1); hold off;
    ylabel('CH_4 flux, w''c'', uM/m^2/d','fontsize',14);  xlabel('Time / Date','fontsize',14);
    % ylim([-100,100]);
    % xlim([dt_3hr(1),dt_3hr(end)]);
    xlim([datetime(2020,6,15),datetime(2020,9,30)]);
ppath = fullfile(plot_path,'co2_ch4_flux_legs4-5_metcity.png');
print(ppath,'-dpng');

% CO2/CH4 flux distribution
fff = sflag_co2_sum_hr>2;
ggg = sflag_ch4_sum_hr>2;
figure('position',[500,100,800,900],'Name','co2_ch4_flux_dist');
subplot(2,1,1);
histogram(Fwco2_csp_mM_hr(fff),[-20:0.2:20]); grid on; grid minor;
xlabel('CO_2 Flux, mM/m^2/d','fontsize',12);
ylabel('N_{tot}','fontsize',12); %ylim([0,60]);
title('MOSAiC, MetCity CO_2/CH_4 Hourly Flux Distributions','fontsize',18);
subtitle('N > 2 good segments / hr','fontsize',12);
subplot(2,1,2);
histogram(Fwch4_csp_uM_hr(ggg),[-80:2.5:80]); grid on; grid minor;
xlabel('CH_4 Flux, \muM/m^2/d','fontsize',12);
ylabel('N_{tot}','fontsize',12); %ylim([0,120]);
ppath = fullfile(plot_path,'co2_ch4_flux_dist_metcity.png');
print(ppath,'-dpng');

figure('position',[500,100,800,900],'Name','Rib_dist');
subplot(2,1,1);
histogram(mc_flx_hr.Rib_hr,[-0.3:0.02:0.3]); grid on; grid minor;
xlabel('Ri bulk @ 10m','fontsize',12);
ylabel('N_{tot}','fontsize',12); ylim([0,2500]);
title('Ri Distributions','fontsize',18);
subplot(2,1,2);
histogram(mc_flx_hr.Ri_L2,[-0.3:0.02:0.3]); grid on; grid minor;
xlabel('Ri (6m-10m levels)','fontsize',12);
ylabel('N_{tot}'); ylim([0,2500]);
ppath = fullfile(plot_path,'Ri_distributions_metcity.png');
print(ppath,'-dpng');

% CO2 flux and Ri_bulk - upper limit for Ri = 0.1
cc = mc_flx_hr.sflag_co2_sum_hr>3 & mc_flx_hr.Rib_hr <= 0.1;
ccc = mc_flx_hr.sflag_co2_sum_hr>3 & mc_flx_hr.Rib_hr > 0.1;
figure('position',[100,750,800,500],'Name','co2_flux_Rib_timeseries');
plot(dt_1hr(cc),mc_flx_hr.Fwco2_csp_mM_hr(cc),'b.','markersize',15); hold on; grid on; grid minor;
plot(dt_1hr(ccc),mc_flx_hr.Fwco2_csp_mM_hr(ccc),'r.','markersize',15);
plot([dt_1hr(1),dt_1hr(end)],[0,0],'k-','linewidth',1); hold off;
ylabel('CO_2 flux, w''c'', mM/m^2/d','fontsize',14);
xlim([dt_1hr(1),dt_1hr(end)]); ylim([-7,10]);
legend('Ri $\le$ 0.1','Ri $>$ 0.1','location','northeast','interpreter','latex','fontsize',14);
title('MOSAiC, MetCity Hourly Mean Carbon Dioxide Flux','fontsize',16);
ppath = fullfile(plot_path,'co2Flux_Rib_limit_metcity.png');
print(ppath,'-dpng');
ppath = fullfile(plot_path,'co2Flux_Rib_limit_metcity.fig');
savefig(ppath)

% CH4 flux and Ri_bulk - upper limit for Ri = 0.1
mm = mc_flx_hr.sflag_ch4_sum_hr>3 & mc_flx_hr.Rib_hr <= 0.1;
mmm = mc_flx_hr.sflag_ch4_sum_hr>3 & mc_flx_hr.Rib_hr > 0.1;
figure('position',[100,750,800,500],'Name','ch4_flux_Rib_timeseries');
plot(dt_1hr(mm),mc_flx_hr.Fwch4_csp_uM_hr(mm),'b.','markersize',15); hold on; grid on; grid minor;
plot(dt_1hr(mmm),mc_flx_hr.Fwch4_csp_uM_hr(mmm),'r.','markersize',15);
plot([dt_1hr(1),dt_1hr(end)],[0,0],'k-','linewidth',1); hold off;
ylabel('CH_4 flux, w''c'', uM/m^2/d','fontsize',14);
xlim([dt_1hr(1),dt_1hr(end)]); ylim([-40,50]);
legend('Ri $\le$ 0.1','Ri $>$ 0.1','location','northeast','interpreter','latex','fontsize',14);
title('MOSAiC, MetCity Hourly Mean Methane Flux','fontsize',16);
ppath = fullfile(plot_path,'ch4Flux_Rib_limit_metcity.png');
print(ppath,'-dpng');
ppath = fullfile(plot_path,'ch4Flux_Rib_limit_metcity.fig');
savefig(ppath)

% CO2/CH4 flux distribution, Ri < 0.1 limit
fff = mc_flx_hr.sflag_co2_sum_hr>2;
ggg = mc_flx_hr.sflag_ch4_sum_hr>2;
kkk = mc_flx_hr.Rib_hr<0.1;
figure('position',[500,100,800,900]);
subplot(2,1,1);
    histogram(mc_flx_hr.Fwco2_csp_mM_hr(fff),[-20:0.5:20]); grid on; grid minor; hold on;
    histogram(mc_flx_hr.Fwco2_csp_mM_hr(fff&kkk),[-20:0.5:20]); hold off;
    xlabel('CO_2 Flux, mM/m^2/d','fontsize',14);
    ylabel('N_{tot}','fontsize',14); %ylim([0,60]);
    legend('all','Ri < 0.1','Location','northwest','fontsize',12);
    title('MOSAiC, MetCity CO_2/CH_4 Hourly Flux Distributions','fontsize',18);
    subtitle('N > 2 good segments / hr','fontsize',14);
subplot(2,1,2);
    histogram(mc_flx_hr.Fwch4_csp_uM_hr(ggg),[-80:2.5:80]); grid on; grid minor; hold on;
    histogram(mc_flx_hr.Fwch4_csp_uM_hr(ggg&kkk),[-80:2.5:80]); hold off;
    xlabel('CH_4 Flux, \muM/m^2/d','fontsize',14);
    ylabel('N_{tot}','fontsize',14); %ylim([0,120]);
    legend('all','Ri < 0.1','Location','northwest','fontsize',12);
ppath = fullfile(plot_path,'co2_ch4_flux_dist_metcity.png');
print(ppath,'-dpng');

% heat flux distribution, Ri < 0.1 limit
figure('position',[500,100,800,500]);
% subplot(2,1,1);
    histogram(mc_tower_hr.Hs_10m(fff),[-50:1:50]); grid on; grid minor; hold on;
    histogram(mc_tower_hr.Hs_10m(fff&kkk),[-50:1:50]); hold off;
    xlabel('10m Sensible Heat, W/m^2','fontsize',14);
    ylabel('N_{tot}','fontsize',14); %ylim([0,60]);
    legend('all','Ri < 0.1','Location','northwest','fontsize',12);
    title('MOSAiC, MetCity Heat Flux Hourly Flux Distributions','fontsize',18);
    subtitle('N > 2 good segments / hr','fontsize',14);
% subplot(2,1,2);
%     histogram(mc_tower_hr.Hl(ggg),[-4:0.25:10]); grid on; grid minor; hold on;
%     histogram(mc_tower_hr.Hl(ggg&kkk),[-4:0.25:10]); hold off;
%     xlabel('2m Latent Heat, W/m^2, Licor','fontsize',14);
%     ylabel('N_{tot}','fontsize',14); %ylim([0,120]);
%     legend('all','Ri < 0.1','Location','northwest','fontsize',12);
ppath = fullfile(plot_path,'heat_flux_dist_metcity.png');
print(ppath,'-dpng');

