% b6_compute_fluxes_metcity_8192.m
%{
    Computes 13.65 min (2^13 / 8192 point) stress, heat, CO2 and CH4 EC fluxes for the metcity tower.

    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                      methane wet mole fraction, ppm
        106   ch4d                      methane dry mole fraction, ppm
        107   co2w                      carbon dioxide wet mole fraction, ppm
        108   co2d                      carbon dioxide wet mole fraction, ppm
        109   h20                       water vapor mole fraction, %
        110   dco2_dt                   co2 rate of change, ppm/hr
        111   dch4_dt                   ch4 rate of change, ppb/hr
        112   sig_co2                   std dev co2 dry mole fraction, ppm
        113   sig_ch4                   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                   std dev co2 (despiked), ppm
        118   ch4_lag                   'optimal' lag for ch4, points @ 10Hz
        119   ch4_std                   std dev ch4 (despiked), ppb
        120   max_Corr_co2              maximum correlation coef, w:co2
        121   max_Corr_ch4              maximum correlation coef, w:ch4
        122   q_pic                     specific humidity, Picarro, g/kg
        123   q_pic_tau                 inlet attn time const for Picarro q
        124   q_pic_lag                 'optimal' lag for Picarro q, points @ 10Hz
        125   q_pic_std                 std dev Picarro q (despiked), g/kg
        126   puff_stpt                 array index for start of puff segment @ 36000 pnts/hr
        127   puff_endpt                array index for end of puff segment @ 36000 pnts/hr

        This script appends the following 54 variables for each flux time segment in the mc_flx table:
        128   wu_cov                    w'u' covariance
        129   wu_csp                    w'u' cospectral integral
        130   wv_cov                    w'v' covariance
        131   wv_csp                    w'v' cospectral integral
        132   wTs_cov                   w'Ts' sonic temperature covariance
        133   wTs_csp                   w'Ts' sonic temperature cospectral integral
        134   wT_cov                    w'T' temperature covariance
        135   wT_csp                    w'T' temperature cospectral integral
        136   std_wT                    std dev wT covariance (Finlkelstein et al.)
        137   Ts_noise                  Sonic temperature noise variance
        138   hsens_cov                 sensible heat flux from wT_cov, W/m2
        139   hsens_csp                 sensible heat flux from wT_csp, W/m2
        140   std_hs                    std dev sensible heat flux (Finlkelstein et al.)
        141   wq_cov_g_kg               w'q' covariance, g/kg m/s
        142   wq_csp_g_kg               w'q' cospectral integral, g/kg m/s
        143   wq_cov_corr_g_kg          w'q' cov attenuation corrected, g/kg m/s
        144   wq_csp_corr_g_kg          w'q' csp attenuation corrected, g/kg m/s
        145   std_wq_g_kg               std dev wq covariance (Finlkelstein et al.)
        146   hlat_cov_W_m2             cov latent heat flux, W/m2
        147   hlat_csp_W_m2             csp latent heat flux, W/m2
        148   std_hl_W_m2               std dev latent heat flux (Finlkelstein et al.)
        149   DLwq_g_kg                 w'q' detection limit, g/kg m/s
        150   DLhlat_W_m2               latent heat detection limit, W/m2
        151   Aq                        q attenuation correction factor
        152   wco2_cov_ppm              w'co2' covariance, ppm m/s
        153   wco2_csp_ppm              w'co2' cospectral integral, ppm m/s
        154   wco2_cov_corr_ppm         w'co2' cov attenuation corrected ppm m/s
        155   wco2_csp_corr_ppm         w'co2' csp attenuation corrected ppm m/s
        156   wco2_cov_corr_mM          w'co2' cov attn corrected mMoles/m2/day
        157   wco2_csp_corr_mM          w'co2' csp attn corrected mMoles/m2/day
        158   std_wco2_mM               std dev wco2 covariance (Finlkelstein et al.)
        159   DLwco2_ppm                w'co2' detection limit, ppm m/s
        160   DLwco2_mM                 w'co2' detection limit, mM/m2/d
        161   uco2_mM                   u'co2' streamwise covariance attn corrected, mM/m2/d
        162   vco2_mM                   v'co2' crosswind covariance attn corrected, mM/m2/d
        163   Aco2                      attenuation factor, CO2
        164   wch4_cov_ppb              w'ch4' covariance ppm m/s
        165   wch4_csp_ppb              w'ch4' cospectral integral ppm m/s
        166   wch4_cov_corr_ppb         w'ch4' cov attenuation corrected ppm m/s
        167   wch4_csp_corr_ppb         w'ch4' csp attenuation corrected ppm m/s
        168   wch4_cov_corr_uM          w'ch4' cov attn corrected, uMoles/m2/d
        169   wch4_csp_corr_uM          w'ch4' csp attn corrected, uMoles/m2/d
        170   std_wch4_uM               std dev wch4 covariance (Finlkelstein et al.)
        171   DLwch4_ppb                w'ch4' detection limit, ppm m/s
        172   DLwch4_uM                 w'ch4' detection limit, uMoles/m2/d
        173   uch4_uM                   u'ch4', streamwise covariance attn corrected, uM/m2/d
        174   vch4_uM                   v'ch4' crosswind covariance attn corrected, uM/m2/d
        175   Ach4                      attenuation factor, CH4
        176   std_w                     std dev w vertical wind velocity
        177   std_u                     std dev u streamwise wind velocity
        178   std_v                     std dev v crosswind velocity
        179   std_wdir                  std dev wdir
        180   L_mean                    hourly Obukhov length from mean fluxes @ 2/6/10m levels for intervals passing mc wind filtering, m
        181   zet10                     hourly z/L from mean flux @10m tower level for intervals passing mc wind filtering, unitless
 
%}

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

%% adjust these parameters as desired, then run

plotit = true;     % set true to save plots - script runs faster if plotting is turned off
len = 8192;        % number of points in each flux window
wq_cutoff = datenum(2020,04,16);    % last datenum for picarro wq calcs - dryer installed
low_pass = true;  % use the low pass filtered co2/ch4 concentrations for computing flux?

% other constants
z = 10;        % measurement height, m, met city tower
fson = 10;     % sampling frequency

%{
    Define start/stop indices for 7 overlapping 8192 point windows in 1 hr of 10Hz data
    (36000 total points/hr) Windows are 13.653 minutes (2^13 data points each).
             600: 8791
            4696:12887
            8792:16983
           12888:21079
           16984:25175
           21080:29271
           25176:33367
    Flux will be computed for these windows each hour.
    The windows specifically exclude the last two minutes and first minute of each hour
    due to the 'puff' which happens some time in this time range.
%}
Nbins = 7;
bins = zeros(Nbins,2);
bins(:,1) = 600:4096:34800-8192;
bins(:,2) = bins(:,1) + 8191;

%% paths
root_data_path = '/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019';
path_prog = fullfile(root_data_path,'_matlab');

picarro_path = fullfile(root_data_path,'crds_ghg_ice_station_tower');
tower_path = fullfile(root_data_path,'_processed','metcity_tower_level2_v3');
wind_path = fullfile(tower_path,'10Hz_turb_L2V3');

if low_pass
    save_path = fullfile(root_data_path,'_processed','metcity_flux_v2','flux_8192_lp');
else
    save_path = fullfile(root_data_path,'_processed','metcity_flux_v2','flux_8192');
end
mkdir(save_path);

plot_path = fullfile(save_path,'plots'); mkdir(plot_path);
path_spectraPlots = fullfile(plot_path,'spectra_cospectra'); mkdir(path_spectraPlots);
path_wucsp = fullfile(path_spectraPlots,'wucsp'); mkdir(path_wucsp);
path_wvcsp = fullfile(path_spectraPlots,'wvcsp'); mkdir(path_wvcsp);
path_wTcsp = fullfile(path_spectraPlots,'wTcsp'); mkdir(path_wTcsp);
path_wco2csp = fullfile(path_spectraPlots,'wco2csp'); mkdir(path_wco2csp);
path_wch4csp = fullfile(path_spectraPlots,'wch4csp'); mkdir(path_wch4csp);
path_wqcsp = fullfile(path_spectraPlots,'wqcsp'); mkdir(path_wqcsp);
path_uvspc = fullfile(path_spectraPlots,'uvspc'); mkdir(path_uvspc);
path_Tspc = fullfile(path_spectraPlots,'Tspc'); mkdir(path_Tspc);
path_wspc = fullfile(path_spectraPlots,'wspc'); mkdir(path_wspc);
path_co2spc = fullfile(path_spectraPlots,'co2spc'); mkdir(path_co2spc);
path_ch4spc = fullfile(path_spectraPlots,'ch4spc'); mkdir(path_ch4spc);
path_qspc = fullfile(path_spectraPlots,'qspc'); mkdir(path_qspc);
path_attn = fullfile(path_spectraPlots,'attn_csp'); mkdir(path_attn);

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

%% load mc_flx table and append new columns

% load flux windows table 'mc_flx' with 127 variables
load(fullfile(save_path,'wins_met_gases_lags_8192_all.mat'));
N = height(mc_flx);

% load cell array of variable names, units and long names,
% including new variables to be added by this script.  Variable
% names and attributes are maintained in an external excel file.
vars = readcell(fullfile(save_path,'flux_vars_V2.xlsx'));

% update existing variable names in mc_flx from excel file
mc_flx.Properties.VariableNames = vars(1:127,1);

% names of 54 new variables added to mc_flux in columns 128-181
% add new variables to mc_flux table
temp = array2table(NaN(N,54),'VariableNames',vars(128:181,1));
mc_flx = [mc_flx,temp]; clear temp;

% define units and long names for all variables
% long names will be copied into the standard 'VariableDescriptions'
% property.
mc_flx.Properties.VariableUnits = vars(:,2);
mc_flx.Properties.VariableDescriptions = vars(:,3);

% load hourly tower flux file for MetCity, will use this for z/L estimate
load(fullfile(tower_path,'metcity_1hr_corA10_all.mat'));
mc_tower_hr = mc_flux_hr; clear mc_flux_hr;

% freq array for computing model inlet attenuation cospectra:
% constant delta log(f), Npts=50.
dlogf = abs(log10(5)-log10(1e-4))/50;
log_f = log10(1e-4):dlogf:log10(5);
fx = 10.^log_f;

% load co2/ch4 low-pass filter transfer function data if needed
% this will be used to remove the Picarro oscillation in the observed
% co2 and ch4 concentrations
if low_pass
    load(fullfile(root_data_path,'_processed','ship_flux_V3','TF_lowpass.mat'),'TFlo');
    Hlp = interp1(TFlo.f,TFlo.h,fx); % this is the lp filter transfer function applied to co2/ch4
    Hlp = replace_NaN_nearest_neighbor(Hlp);  % in case there's a NaN in the interp
end

%% set up structures to store smoothed spectra/cospectra

spec_sm = struct('datenum',NaN(N,1));
spec_sm.datenum = mc_flx.datenum;
[spec_sm(1).Fx] = NaN(1,50);
[spec_sm(1).Pu] = NaN(N,50); [spec_sm(1).Pv] = NaN(N,50);
[spec_sm(1).Pw] = NaN(N,50); [spec_sm(1).Pt] = NaN(N,50);
[spec_sm(1).Pt2] = NaN(N,50); [spec_sm(1).Pq] = NaN(N,50);
[spec_sm(1).Pc] = NaN(N,50); [spec_sm(1).Pm] = NaN(N,50);
[spec_sm(1).Cuw] = NaN(N,50); [spec_sm(1).Cvw] = NaN(N,50);
[spec_sm(1).Ctw] = NaN(N,50); [spec_sm(1).Cqw] = NaN(N,50);
[spec_sm(1).Cco2w] = NaN(N,50); [spec_sm(1).Cch4w] = NaN(N,50);

%%  flux loop - process fluxes in 8192 point data windows, 7 windows per hour

current_day = 0; % initialize variable identifying the current day
for kk = 35155:7:N  % indices of whole hours in flx table to process

    %% setup this hour

    yyyy = mc_flx.yyyy(kk); mm = mc_flx.mm(kk); dd = mc_flx.dd(kk); hh = mc_flx.hh(kk);
    YYYY = sprintf('%4i',yyyy); YY = YYYY(3:4);
    MM = sprintf('%02i',mm);
    DD = sprintf('%02i',dd);
    HH = sprintf('%02i',hh);
    ddd = md2yd_vn(yyyy,mm,dd);
    DDD = sprintf('%03i',ddd);
    dnum_this_hour = datenum(mc_flx.yyyy(kk),mc_flx.mm(kk),mc_flx.dd(kk),mc_flx.hh(kk),0,0);
    this_day = floor(dnum_this_hour);
    LEG = sprintf('%01i',mc_flx.leg(kk));

    % mc_flux indices for this hour
    zz = kk:1:kk+6;

    % get observed for 10m z/L for this hour, compute L
    % hourly mc tower met is filtered for basic wind direction
    aa = find(abs(mc_tower_hr.datenum  - dnum_this_hour)<5/86400); % find hour, 5 sec tolerance
    if isempty(aa)
        zet10 = NaN;
        L_mean = NaN;
    else
        zet10 = mc_tower_hr.zeta_level_n_10m(aa); % leg 1 z_sonic = 10.34m
        L_mean = mc_tower_hr.L_all_levels(aa); 
    end
    mc_flx.L_mean(zz) = L_mean;
    mc_flx.zet10(zz) = zet10;

    if this_day~=current_day  % if new day load 10Hz sonic wind - one file for entire day
        %{
        son is 864000 x 42 table with the following variables
        1    dt                     datetime timestamp
        2    dnum                   datenum timestamp
        3    metek_2m_heatstatus    2m sonic heater status flag
        4    metek_2m_u             2m sonic U (+W) m/s
        5    metek_2m_v             2m sonic V (+N) m/s
        6    metek_2m_w             2m sonic W (+up) m/s
        7    metek_2m_T             2m sonic temperature, °C
        8    metek_2m_hspd          2m horizontal wind speed, m/s
        9    metek_2m_ts            2m wind direction, deg
        10   metek_2m_incx          2m pitch angle, + is antiCW about sonic E-W axis viewing E to W'
        11   metek_2m_incy          2m roll angle, + is antiCW about sonic S-N axis viewing S to N'
        12   metek_6m_heatstatus    same for 60 sonic
        13   metek_6m_u             ...
        14   metek_6m_v             ...
        15   metek_6m_w             ...
        16   metek_6m_T             ...
        17   metek_6m_hspd          ...
        18   metek_6m_ts            ...
        19   metek_6m_incx          ...
        20   metek_6m_incy          ...
        21   metek_10m_heatstatus   same for 10m sonic
        22   metek_10m_u            ...
        23   metek_10m_v            ...
        24   metek_10m_w            ...
        25   metek_10m_T            ...
        26   metek_10m_hspd         ...
        27   metek_10m_ts           ...
        28   metek_10m_incx         ...
        29   metek_10m_incy         ...
        30   metek_mast_u           32/23m mast U (+W) m/s
        31   metek_mast_v           32/23m mast V (+N) m/s
        32   metek_mast_w           32/23m mast W (+up) m/s
        33   metek_mast_T           32/23m mast sonic temperature, °C
        34   licor_diag             Licor 7500 diagnostic flag, bits 0-3 = signal strength; bit 5 = PLL; bit 6 = detector temp; bit 7 = chopper temp
        35   licor_co2              Licor CO2, mmol/m3
        36   licor_h2o              Licor h2o, mmol/m3
        37   licor_T                Licor temperature at strut, °C
        38   licor_pr               Licor pressure at heatd, hPa
        39   licor_co2_str          Licor CO2 signal strength, raw co2 reference signal relative to expected value
        40   licor_pll              Licor phase lock loop, 1 = optical filter wheel rotating at correct rate (boolean)
        41   licor_dt               Licor detector temperature, 1 = temperature near set point (boolean)
        42   licor_ct               Licor chopper temperature, 1 = temperature near set point (boolean)
        %}
        dfl_wind = fullfile(wind_path,['metcity_wind_',YYYY,MM,DD,'.mat']);
        if exist(dfl_wind,'file')
            load(dfl_wind);
            current_day = this_day;
        else
            continue % no wind data, so continue to next hour
        end
    end

    % load various parameters for each 8192 point segment this hour
    T_K = mc_flx.temp_10m(zz) + 273.15;
    RH = mc_flx.rh_10m(zz);
    PPa = mc_flx.atmos_pressure_2m(zz)*100;
    rwspd = mc_flx.wspd_10m(zz);
    puff_stpt = mc_flx.puff_stpt(zz);
    puff_endpt = mc_flx.puff_endpt(zz); if puff_endpt==1; puff_endpt=2; end
    tau_q = mc_flx.q_pic_tau(zz);
    q_lag = mc_flx.q_pic_lag(zz);
    tau_co2_ch4 = mc_flx.co2_ch4_tau(zz);
    co2_lag = mc_flx.co2_lag(zz);
    ch4_lag = mc_flx.ch4_lag(zz);

    %% process wind, compute stress and sensible heat flux
        
    % index for wind this hour
    st_idx = find(abs(son.dnum - dnum_this_hour)<0.5/864000);
    ww = st_idx:st_idx + 35999;
    
    % if no wind data this hour, continue to next hour
    if all(isnan(son.metek_10m_v(ww)))
        continue
    end
    
    % switch POP uv wind convention to BLO standard
    u = son.metek_10m_v(ww);      % U wind positive to N
    v = -son.metek_10m_u(ww);     % V wind positive to W
    w = son.metek_10m_w(ww);      % vertical wind
    Ts = son.metek_10m_T(ww);     % sonic temperature, °C
    badw = length(find(isnan(w)));
    % if too many bad points, skip this hour and set gap flag
    if badw>100
        continue  % skip this hour if more than 10 sec of NaN in wind data
    end
    
    % Split wind arrays into overlapping, despiked, detrended flux window bins.
    % Rotate uv into azimuth so u is streamwise wind speed and mean v = 0.
    % We do not do a double rotation to get mean w = 0 because there is a
    % question about persistent W bias in the metek sonic data even though
    % the tilt angle is very small.
    % For now we are ignoring the W bias and will detrend it away before computing fluxes.
    % Also, compute std dev of the 10Hz wdir.
    uvwt_win = NaN(len,4,Nbins);  % len x uvwt x #windows
    for ii = 1:Nbins
        [u_str,v_str] = uv_rot(u(bins(ii,1):bins(ii,2)),v(bins(ii,1):bins(ii,2)));
        uvwt_win(:,1,ii) = despike2(detrend(u_str,'linear'));
        uvwt_win(:,2,ii) = despike2(detrend(v_str,'linear'));
        uvwt_win(:,3,ii) = despike2(detrend(w(bins(ii,1):bins(ii,2)),'linear'));
        uvwt_win(:,4,ii) = despike2(detrend(Ts(bins(ii,1):bins(ii,2)),'linear'));
        [~,wd] = uv_to_sd(u(bins(ii,1):bins(ii,2)),v(bins(ii,1):bins(ii,2)));
        wd = unwrap(wd*pi/180)*180/pi;   % unwrapped 10Hz wind angle
        mc_flx.std_wdir(zz(ii)) = nanstd1(wd);   % save std wdir to output array
    end

    % compute sonic covariances and save to output array
    for ii = 1:Nbins
        % wu_cov
        temp = nancov1(uvwt_win(:,1,ii),uvwt_win(:,3,ii));
        mc_flx.wu_cov(zz(ii)) = temp(1,2);
        % wv_cov
        temp = nancov1(uvwt_win(:,2,ii),uvwt_win(:,3,ii));
        mc_flx.wv_cov(zz(ii)) = temp(1,2);
        % wTs_cov
        temp = nancov1(uvwt_win(:,4,ii),uvwt_win(:,3,ii));
        mc_flx.wTs_cov(zz(ii)) = temp(1,2);
        % sigma wTs (after Finkelstein et al.)
        mc_flx.std_wT(zz(ii)) = sqrt(var_of_cov(uvwt_win(:,4,ii),uvwt_win(:,3,ii),200));
        % sigma w, u, v
        mc_flx.std_w(zz(ii)) = nanstd1(uvwt_win(:,3,ii));
        mc_flx.std_u(zz(ii)) = nanstd1(uvwt_win(:,1,ii));
        mc_flx.std_v(zz(ii)) = nanstd1(uvwt_win(:,2,ii));
    end

    %% compute sonic spectra & cospectra
    % initialize arrays for this hour's spectral data
    Fx = NaN(1,50);        % empty arrays for 6x50-bin smoothed spectra
    Pu = NaN(Nbins,50);
    Pv = Pu; Pw = Pu; Pt = Pu; Pt2 = Pu;
    Cuw = Pu; Cvw = Pu; Ctw = Pu;
    Ts_noise = NaN(Nbins,1);

    Fxx = NaN(1,len/2+1);  % empty arrays for raw spectra
    Puu = NaN(Nbins,len/2+1);
    Pvv = Puu; Pww = Puu; Ptt = Puu; Ptt2 = Puu;
    Cuuw = Puu; Cvvw = Puu; Cttw = Puu;

    df = fson/(len-1);

    for ii = 1:Nbins
        % Ustr spectrum
        [Puu(ii,:),Fxx(1,:)] = psd2(uvwt_win(:,1,ii),len,fson,hamming(len));
        Puu(ii,:) = real(Puu(ii,:));
        [Pu(ii,:),Fx(1,:)] = specsmoo(Puu(ii,:),fson);
        % V_str spectrum
        [Pvv(ii,:),~] = psd2(uvwt_win(:,2,ii),len,fson,hamming(len));
        Pvv(ii,:) = real(Pvv(ii,:));
        [Pv(ii,:),~] = specsmoo(Pvv(ii,:),fson);
        % W_str spectrum
        [Pww(ii,:),~] = psd2(uvwt_win(:,3,ii),len,fson,hamming(len));
        Pww(ii,:) = real(Pww(ii,:));
        [Pw(ii,:),~] = specsmoo(Pww(ii,:),fson);
        % Ts spectrum
        [Ptt(ii,:),~] = psd2(uvwt_win(:,4,ii),len,fson,hamming(len));
        Ptt(ii,:) = real(Ptt(ii,:));
        [Pt(ii,:),~] = specsmoo(Ptt(ii,:),fson);

        % Subtract noise from Tsonic spectra, assuming sonic spectrum at ~5 Hz is all noise
        % (assume last 3 points in smoothed spectrum are just noise)
        Ts_noise(ii) = median(Pt(ii,end-2:end));
        Ptt2(ii,:) = Ptt(ii,:) - Ts_noise(ii);
        [Pt2(ii,:),~] = specsmoo(Ptt2(ii,:),fson);
        jj = find(Pt2(ii,:)<0); Pt2(ii,jj) = 0; % be sure there are no negative variances

        % cospectra
        [Cuuw(ii,:),~] = csd2(uvwt_win(:,3,ii),uvwt_win(:,1,ii),len,fson,hamming(len));
        Cuuw(ii,:) = real(Cuuw(ii,:));
        % integrate cospectrum for wu_csp
        mc_flx.wu_csp(zz(ii)) = sum(Cuuw(ii,:)*df);
        [Cuw(ii,:),~] = specsmoo(Cuuw(ii,:),fson);
        [Cvvw(ii,:),~] = csd2(uvwt_win(:,3,ii),uvwt_win(:,2,ii),len,fson,hamming(len));
        Cvvw(ii,:) = real(Cvvw(ii,:));
        % integrate cospectrum for wv_csp, col 169
        mc_flx.wv_csp(zz(ii)) = sum(Cvvw(ii,:)*df);
        [Cvw(ii,:),~] = specsmoo(Cvvw(ii,:),fson);
        [Cttw(ii,:),~] = csd2(uvwt_win(:,3,ii),uvwt_win(:,4,ii),len,fson,hamming(len));
        Cttw(ii,:) = real(Cttw(ii,:));
        % integrate cospectrum for wTs_csp, col 170
        mc_flx.wTs_csp(zz(ii)) = sum(Cttw(ii,:)*df);
        [Ctw(ii,:),~] = specsmoo(Cttw(ii,:),fson);
        % not making q correction, so wT = wTs for now
        mc_flx.wT_csp(zz(ii)) = mc_flx.wTs_csp(zz(ii));
    end

    % save Ts_noise for segments this hour
    mc_flx.Ts_noise(zz) = Ts_noise;

    %% save smoothed spectra/cospectra to spec structure
    if all(isnan(spec_sm(1).Fx)); spec_sm(1).Fx = Fx; end

    spec_sm(1).Pu(zz,:) = Pu;
    spec_sm(1).Pv(zz,:) = Pv;
    spec_sm(1).Pw(zz,:) = Pw;
    spec_sm(1).Pt(zz,:) = Pt;
    spec_sm(1).Pt2(zz,:) = Pt2;

    spec_sm(1).Cuw(zz,:) = Cuw;
    spec_sm(1).Cvw(zz,:) = Cvw;
    spec_sm(1).Ctw(zz,:) = Ctw;

    %% plot wind spectra & cospectra

    if plotit
        % cospectra plots
        if exist('h1','var')
    		clf(h1)
            set(0,'CurrentFigure', h1);
        else
            h1 = figure('position',[500,100,700,500],'Name','wu_csp');
        end
        wu = nanmean1(Cuw,1);
        semilogx(Fx,Fx.*wu,'b-',[1e-3 10 ],[0 0],'k-','linewidth',2);
        title([MM,'/',DD,'/',YYYY,' ',HH,'00',' wu Cospectrum']); grid;
        ylabel('C_{uw}(f) x f'); xlabel('f (Hz)');
        ppath = fullfile(path_wucsp,[YYYY,'_',MM,'_',DD,'_',HH,'00_wuCospec.png']);
        print(ppath,'-dpng');

        if exist('h2','var')
        	clf(h2)
            set(0,'CurrentFigure', h2);
        else
            h2 = figure('position',[500,100,700,500],'Name','wv_csp');
        end
        wv = nanmean1(Cvw,1);
        semilogx(Fx,Fx.*wv,'b-',[1e-3 10 ],[0 0],'k-','linewidth',2);
        title([MM,'/',DD,'/',YYYY,' ',HH,'00',' wv Cospectrum']); grid;
        ylabel('C_{vw}(f) x f'); xlabel('f (Hz)');
        ppath = fullfile(path_wvcsp,[YYYY,'_',MM,'_',DD,'_',HH,'00_wvCospec.png']);
        print(ppath,'-dpng');

        if exist('h3','var')
        	clf(h3)
            set(0,'CurrentFigure', h3);
        else
            h3 = figure('position',[500,100,700,500],'Name','wt_csp');
        end
        wt = nanmean1(Ctw,1);
        semilogx(Fx,Fx.*wt,'b-',[1e-3 10 ],[0 0],'k-','linewidth',2);
        title([MM,'/',DD,'/',YYYY,' ',HH,'00',' wTs Cospectrum']); grid;
        ylabel('C_{wt}(f) x f'); xlabel('f (Hz)');
        ppath = fullfile(path_wTcsp,[YYYY,'_',MM,'_',DD,'_',HH,'00_wtCospec.png']);
        print(ppath,'-dpng');

        % spectra plots
        if exist('h4','var')
        	clf(h4)
            set(0,'CurrentFigure', h4);
        else
            h4 = figure('position',[500,100,700,500],'Name','w_spc');
        end
        Sw = nanmean1(Pw,1);
        loglog(Fx,Fx.*Sw,'b-','linewidth',2); hold on;
        loglog(Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*Sw(35).*Fx(35),'r');
        hold off; ylabel('S_{ww} x f'); xlabel('f (Hz)');
        title([MM,'/',DD,'/',YYYY,' ',HH,'00',' Spectrum of vertical wind velocity, W']); grid;
        ppath = fullfile(path_wspc,[YYYY,'_',MM,'_',DD,'_',HH,'00_wSpec.png']);
        print(ppath,'-dpng');

        if exist('h5','var')
        	clf(h5)
            set(0,'CurrentFigure', h5);
        else
            h5 = figure('position',[500,100,700,500],'Name','t_spc');
        end
        St = nanmean1(Pt,1); St2 = nanmean1(Pt2,1);
        WN = zeros(size(St)) + nanmean1(Ts_noise); WN(1:24) = NaN;
        loglog(Fx,Fx.*St,'b--',Fx(1:end-5),Fx(1:end-5).*St2(1:end-5),'b-','linewidth',2); hold on;
        loglog(Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*St(35).*Fx(35),'r-', Fx,Fx.*WN,'k--');
        hold off; ylabel('S_{tt} x f'); xlabel('f (Hz)'); title([MM,'/',DD,'/',YYYY,' ',HH,'00',' Spectrum of sonic temperature, Ts']); grid;
        legend('raw','minus noise','-2/3','white noise','location','northwest');
        ppath = fullfile(path_Tspc,[YYYY,'_',MM,'_',DD,'_',HH,'00_tSpec.png']);
        print(ppath,'-dpng');

        if exist('h6','var')
        	clf(h6)
            set(0,'CurrentFigure', h6);
        else
            h6 = figure('position',[500,100,700,500],'Name','uv_spc');
        end
        Su = nanmean1(Pu,1);
        Sv = nanmean1(Pv,1);
        loglog(Fx,Fx.*Su,'b-',Fx,Fx.*Sv,'g-','linewidth',2); hold on;
        loglog(Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*Su(35).*Fx(35),'r');
        hold off; ylabel('S_{uu} x f'); xlabel('f (Hz)');
        legend('U','V','-2/3','location','northeast');
        title([MM,'/',DD,'/',YYYY,' ',HH,'00',' U V Wind Spectra']); grid;
        ppath = fullfile(path_uvspc,[YYYY,'_',MM,'_',DD,'_',HH,'00_uvSpec.png']);
        print(ppath,'-dpng');
    end

    %% process Picarro CO2/CH4/H2O file
    
    co2_dfl = fullfile(picarro_path,[YY,DDD],['pic1',YY,DDD,HH,'_raw.mat']);
    if ~exist(co2_dfl,'file')
        continue % no Picarro data this hour, so continue to next hour
    end
    
    %{
        Reads serial feed file from DAS100 specified by dfl and
        return 36000 x 12 array with the following columns:
        1     jd   decimal day-of-year timestamp
        2    sec   Picarro timestamp, seconds since 1/1/1970 00:00
        3     Pc   cavity pressure, torr
        4     Tc   cavity temperature, °C
        5     V1   solenoid valve flag
        6   ch4w   methane conc, ppm, wet
        7   ch4d   methane conc, ppm, dry
        8   co2w   carbon dioxide conc, ppm, wet
        9   co2d   carbon dioxide conc, ppm, dry
        10   h20   water vapor conc, %
        11  ch4lp  low-pass filtered ch4d, 0.8 Hz passband
        12  co2lp  low-pass filtered co2d, 0.8 Hz passband
    %}
    pic = read_picarro_mosaic_2019(co2_dfl,yyyy,ddd,hh);
    q = pic(:,10)*0.01*0.622*1000; % percent mole fraction to g/kg
    % select raw or low-pass filtered co2/ch4 data
    if low_pass
        co2 = pic(:,12);        % low pass filtered co2 dry
        ch4 = pic(:,11);        % loq pass filtered ch4 dry
    else
        co2 = pic(:,9);         % 'dry' gas mole fraction, ppm
        ch4 = pic(:,7);         % 'dry' gas mole fraction, ppm
    end
    
    
    % process Picarro data if fewer than 50 bad points this hour
    if sum(isnan(q))<50 && sum(isnan(co2))<50 && sum(isnan(ch4))<50 && isfinite(puff_stpt(1)) && co2_lag(1)>1
        
        %% pad puff region with constant value
        if (puff_stpt(1)>0 && puff_endpt(1)>0)
            q(1:puff_endpt(1)) = q(puff_endpt(1)+1);
            co2(1:puff_endpt(1)) = co2(puff_endpt(1)+1);
            ch4(1:puff_endpt(1)) = ch4(puff_endpt(1)+1);
        elseif (puff_stpt(1)<0 && puff_endpt(1)<0)
            q(puff_stpt(1)+36000:end) = q(puff_stpt(1)+36000-1);
            co2(puff_stpt(1)+36000:end) = co2(puff_stpt(1)+36000-1);
            ch4(puff_stpt(1)+36000:end) = ch4(puff_stpt(1)+36000-1);
        elseif (puff_stpt(1)<0 && puff_endpt(1)>0)
            q(1:puff_endpt(1)) = q(puff_endpt(1)+1);
            co2(1:puff_endpt(1)) = co2(puff_endpt(1)+1);
            ch4(1:puff_endpt(1)) = ch4(puff_endpt(1)+1);
            q(puff_stpt(1)+36000:end) = q(puff_stpt(1)+36000-1);
            co2(puff_stpt(1)+36000:end) = co2(puff_stpt(1)+36000-1);
            ch4(puff_stpt(1)+36000:end) = ch4(puff_stpt(1)+36000-1);
        end
        
        % synchronize picarro and wind arrays
        % shift trace gas arrays by indicated lag (in points)
        if isfinite(q_lag(1))
            q_seg = q;
            q_seg(1:end-q_lag(1)+1) = q(q_lag(1):end);
            q_seg = despike2(q_seg);
        else
            q_seg = NaN(36000,1);
        end
        
        if isfinite(co2_lag(1))
            co2_seg = co2;
            co2_seg(1:end-co2_lag(1)+1) = co2(co2_lag(1):end);
            co2_seg = despike2(co2_seg);
        else
            co2_seg = NaN(36000,1);
        end
        
        if isfinite(ch4_lag(1))
            ch4_seg = ch4;
            ch4_seg(1:end-ch4_lag(1)+1) = ch4(ch4_lag(1):end);
            ch4_seg = detrend(despike2(ch4_seg));
        else
            ch4_seg = NaN(36000,1);
        end
        
        %% compute covariances and detection limits
        for ii=1:Nbins
            % unit conversions for gas fluxes:
            % ppm m/s to mM/m2/d conversion = 86400*Pa/(8.314*T_K)/1000
            % this same constant also converts concentration units ppb m/s to uM/m2/d
            mM_m2_d_conv = 86400*PPa(ii)/8.31451./T_K(ii)/1000;
            
            % wq_cov and detection limit estimate
            temp = nancov1(detrend(q_seg(bins(ii,1):bins(ii,2)),'linear'),...
                detrend(uvwt_win(:,3,ii),'linear'));
            mc_flx.wq_cov_g_kg(zz(ii)) = temp(1,2); % save wq_cov, g/kg m/s
            mc_flx.DLwq_g_kg(zz(ii)) = flux_DL(detrend(uvwt_win(:,3,ii),'linear'),...
                detrend(q_seg(bins(ii,1):bins(ii,2)),'linear'),20);
            vcv_q = var_of_cov(q_seg(bins(ii,1):bins(ii,2)),uvwt_win(:,3,ii),200);
            mc_flx.std_wq_g_kg(zz(ii)) = sqrt(vcv_q); % wq_cov std dev, Finkelstein 2001
            
            % wco2_cov/uco2_cov/vco2_cov and detection limit
            temp = nancov1(detrend(co2_seg(bins(ii,1):bins(ii,2)),'linear'),...
                detrend(uvwt_win(:,3,ii),'linear'));
            mc_flx.DLwco2_ppm(zz(ii)) = flux_DL(detrend(uvwt_win(:,3,ii),'linear'),...
                detrend(co2_seg(bins(ii,1):bins(ii,2)),'linear'),20);
            mc_flx.wco2_cov_ppm(zz(ii)) = temp(1,2); % save wco2_cov, ppm/m/s
            mc_flx.DLwco2_mM(zz(ii)) = mc_flx.DLwco2_ppm(zz(ii))*mM_m2_d_conv; % save DL in mM units

            temp = nancov1(detrend(co2_seg(bins(ii,1):bins(ii,2)),'linear'),...
                detrend(uvwt_win(:,1,ii),'linear'));
            mc_flx.uco2_mM(zz(ii)) = temp(1,2)*mM_m2_d_conv; % save uco2_cov, mM/m2/d
            temp = nancov1(detrend(co2_seg(bins(ii,1):bins(ii,2)),'linear'),...
                detrend(uvwt_win(:,2,ii),'linear'));
            mc_flx.vco2_mM(zz(ii)) = temp(1,2)*mM_m2_d_conv; % save vco2_cov, mM/m2/d
            vcv_co2 = var_of_cov(detrend(co2_seg(bins(ii,1):bins(ii,2)),'linear'),detrend(uvwt_win(:,3,ii),'linear'),200);
            mc_flx.std_wco2_mM(zz(ii)) = sqrt(vcv_co2)*mM_m2_d_conv; % wco2_cov std dev, Finkelstein 2001
            
            % wch4_cov/uch4_cov/vch4_cov and detection limit
            temp = nancov1(detrend(ch4_seg(bins(ii,1):bins(ii,2)),'linear'),...
                detrend(uvwt_win(:,3,ii),'linear'));
            mc_flx.DLwch4_ppb(zz(ii)) = flux_DL(detrend(uvwt_win(:,3,ii),'linear'),...
                detrend(ch4_seg(bins(ii,1):bins(ii,2)),'linear'),20);
            mc_flx.wch4_cov_ppb(zz(ii)) = temp(1,2); % save wch4_cov  ppm m/s
            mc_flx.DLwch4_uM(zz(ii)) = mc_flx.DLwch4_ppb(zz(ii))*mM_m2_d_conv*1000; % save DL in uM units
            
            temp = nancov1(detrend(ch4_seg(bins(ii,1):bins(ii,2)),'linear'),...
                detrend(uvwt_win(:,1,ii),'linear'));
            mc_flx.uch4_uM(zz(ii)) = temp(1,2)*mM_m2_d_conv*1000; % save uch4_cov, uM/m2/d
            temp = nancov1(detrend(ch4_seg(bins(ii,1):bins(ii,2)),'linear'),...
                detrend(uvwt_win(:,2,ii),'linear'));
            mc_flx.vch4_uM(zz(ii)) = temp(1,2)*mM_m2_d_conv*1000; % save vch4_cov, uM/m2/d
            vcv_ch4 = var_of_cov(detrend(ch4_seg(bins(ii,1):bins(ii,2)),'linear'),detrend(uvwt_win(:,3,ii),'linear'),200);
            mc_flx.std_wch4_uM(zz(ii)) = sqrt(vcv_ch4)*mM_m2_d_conv*1000; % wch4_cov std dev, Finkelstein 2001
            
            % compute Kaimal cospectrum
            if isnan(L_mean)
                Km = kaimal_neutral_scalar_Co(fx,z,rwspd(ii)); % neutral stability model
            else
                Km = moore_scalar_Co(fx,z,L_mean,rwspd(ii)); % stability-dependent cospectral model
            end

            % 1st order linear transfer functions
            Hfq = (1 + (2*pi*fx*tau_q(ii)).^2).^(-1);
            if low_pass
                Hfcc_p = (1 + (2*pi*fx*mc_flx.co2_ch4_tau(zz(ii))).^2).^(-1);  % from puff time const
                Hfcc = Hfcc_p.*Hlp; % combined with lp filter used on co2/ch4
%                 figure; semilogx(fx,Hfcc_p,'r-',fx,Hlp,'b-',fx,Hfcc,'k--','linewidth',1.5); grid on; grid minor;
            else
                Hfcc = (1 + (2*pi*fx*tau_co2_ch4(ii)).^2).^(-1);
            end
            
            % filtered cospectra
            % use sqrt(Hf) since only gas conc is low-pass filtered by inlet, w is not.
            Kq = Km.*sqrt(Hfq);
            Kcc = Km.*sqrt(Hfcc);
            
            % compute attenuation factors as ratio of cospectral areas
            Aq = trapz(fx,Kq./fx)/trapz(fx,Km./fx);  %AQ = sprintf('%5.4f',Aq);
            Ac = trapz(fx,Kcc./fx)/trapz(fx,Km./fx); %AC = sprintf('%5.4f',Ac);
            mc_flx.Aq(zz(ii)) = Aq;
            mc_flx.Aco2(zz(ii)) = Ac; % same value for CO2 and CH4
            mc_flx.Ach4(zz(ii)) = Ac;
            
            % computed corrected fluxes
            mc_flx.wq_cov_corr_g_kg(zz(ii)) = mc_flx.wq_cov_g_kg(zz(ii))/Aq;    % wq_cov_corr, g/kg m/s
            mc_flx.wco2_cov_corr_ppm(zz(ii)) = mc_flx.wco2_cov_ppm(zz(ii))/Ac;  % wco2_cov_corr, ppm/m/s
            mc_flx.wco2_cov_corr_mM(zz(ii)) = mc_flx.wco2_cov_corr_ppm(zz(ii))*mM_m2_d_conv; % mM units
            mc_flx.uco2_mM(zz(ii)) = mc_flx.uco2_mM(zz(ii))/Ac;                 % uco2_mM, corrected
            mc_flx.vco2_mM(zz(ii)) = mc_flx.vco2_mM(zz(ii))/Ac;                 % vco2_mM, corrected
            mc_flx.wch4_cov_corr_ppb(zz(ii)) = mc_flx.wch4_cov_ppb(zz(ii))/Ac;  % wch4_cov_corr, ppm/m/s
            mc_flx.wch4_cov_corr_uM(zz(ii)) = mc_flx.wch4_cov_corr_ppb(zz(ii))*mM_m2_d_conv*1000; %uM units
            mc_flx.uch4_uM(zz(ii)) = mc_flx.uch4_uM(zz(ii))/Ac;                 % uch4_uM, corrected
            mc_flx.vch4_uM(zz(ii)) = mc_flx.vch4_uM(zz(ii))/Ac;                 % vch4_uM, corrected
            
            % latent heat vap for fresh water
            [Le_w, ~] = Le_water(-1.7, 33); % use approx SST and salinity
            % air density kg/m3
            rhoair = rhoa(T_K(ii)-273.15,PPa(ii)/100,RH(ii));
            % compute latent heat flux: hlat = wq * rhoair * Le_w = g/kg m/s * kg/m3 * J/g = W/m2
            mc_flx.hlat_cov_W_m2(zz(ii)) = mc_flx.wq_cov_corr_g_kg(zz(ii)).*rhoair*Le_w;
            % latent heat detection limit
            mc_flx.DLhlat_W_m2(zz(ii)) = mc_flx.DLwq_g_kg(zz(ii)).*rhoair*Le_w;
            % latent heat uncertainty
            mc_flx.std_hl_W_m2(zz(ii)) = mc_flx.std_wq_g_kg(zz(ii)).*rhoair*Le_w;
            
            % correct wTs for humidity flux & compute sensible heat flux
            % this is a very small correction when q is small
            % <wT> = <wTs> - 0.51*T*<wq>, where wq is kg/kg units
%             if isfinite(mc_flx.wq_cov_g_kg(zz(ii)))
%                 mc_flx.wT_cov(zz(ii)) = mc_flx.wq_cov_g_kg(zz(ii)) + 0.51*T_K(ii).*mc_flx.wq_cov_g_kg(zz(ii))/1000;
%                 % convert std_wTs to std_wT
%                 mc_flx.std_wT(zz(ii)) = mc_flx.std_wT(zz(ii)) + 0.51*T_K(ii).*mc_flx.wq_cov_g_kg(zz(ii))/1000;
%             else
                mc_flx.wT_cov(zz(ii)) = mc_flx.wTs_cov(zz(ii)); % but skip q correction if no q covariance...
%             end
            % cov sensible heat flux and uncertainty
            mc_flx.hsens_cov(zz(ii)) = mc_flx.wT_cov(zz(ii)).*rhoair*1004.67;   % hs: Cp_air = 1004.67 J/kg K
            mc_flx.std_hs(zz(ii)) = mc_flx.std_wT(zz(ii)).*rhoair*1004.67;   % std hs
            % csp sensible heat flux
            mc_flx.hsens_csp(zz(ii)) = mc_flx.wT_csp(zz(ii)).*rhoair*1004.67;
        end
        
        %% compute spectra and cospectra
        % initialize arrays for this hour's spectral data
        Fx = NaN(1,50);      % empty arrays for 6x48-bin smoothed spectra
        Pq = NaN(Nbins,50); Pc = Pq; Pm = Pq;
        Cqw = Pq; Cco2w = Pq; Cch4w = Pq;

        df = fson/(len-1);
        
        Fxx = NaN(1,len/2+1);  % empty arrays for raw spectra
        Pqq = NaN(Nbins,len/2+1); Pcc = Pqq; Pmm = Pqq;
        Cqqw = Pqq; Cccw = Pqq; Cmmw = Pqq;
        
        for ii = 1:Nbins
            % q spectrum
            [Pqq(ii,:),Fxx] = psd2(detrend(q_seg(bins(ii,1):bins(ii,2)),'linear'),len,fson,hamming(len));
            Pqq(ii,:) = real(Pqq(ii,:));
            [Pq(ii,:),Fx] = specsmoo(Pqq(ii,:),fson);

            % co2 spectrum
            [Pcc(ii,:),~] = psd2(detrend(co2_seg(bins(ii,1):bins(ii,2)),'linear'),len,fson,hamming(len));
            Pcc(ii,:) = real(Pcc(ii,:));
            [Pc(ii,:),~] = specsmoo(Pcc(ii,:),fson);

            % ch4 spectrum
            [Pmm(ii,:),~] = psd2(detrend(ch4_seg(bins(ii,1):bins(ii,2)),'linear'),len,fson,hamming(len));
            Pmm(ii,:) = real(Pmm(ii,:));
            [Pm(ii,:),~] = specsmoo(Pmm(ii,:),fson);

            % w:q cospectrum
            [Cqqw(ii,:),~] = csd2(uvwt_win(:,3,ii),detrend(q_seg(bins(ii,1):bins(ii,2)),'linear'),...
                                  len,fson,hamming(len));
            Cqqw(ii,:) = real(Cqqw(ii,:));
            % integrate cospectrum for wq_csp, wq_csp_corr & hlat_csp
            mc_flx.wq_csp_g_kg(zz(ii)) = sum(Cqqw(ii,:)*df);
            mc_flx.wq_csp_corr_g_kg(zz(ii)) = mc_flx.wq_csp_g_kg(zz(ii))/Aq;
            mc_flx.hlat_csp_W_m2(zz(ii)) = mc_flx.wq_csp_corr_g_kg(zz(ii)).*rhoair*Le_w;
            [Cqw(ii,:),~] = specsmoo(Cqqw(ii,:),fson);

            % w:co2 cospectrum
            [Cccw(ii,:),~] = csd2(uvwt_win(:,3,ii),detrend(co2_seg(bins(ii,1):bins(ii,2)),'linear'),...
                                  len,fson,hamming(len));
            Cccw(ii,:) = real(Cccw(ii,:));
            % integrate cospectrum for wco2_csp, wco2_csp_corr & wco2_csp_mM
            mc_flx.wco2_csp_ppm(zz(ii)) = sum(Cccw(ii,:)*df);
            mc_flx.wco2_csp_corr_ppm(zz(ii)) = mc_flx.wco2_csp_ppm(zz(ii))/Ac;
            mc_flx.wco2_csp_corr_mM(zz(ii)) = mc_flx.wco2_csp_corr_ppm(zz(ii))*mM_m2_d_conv; % mM/m2/d units
            [Cco2w(ii,:),~] = specsmoo(Cccw(ii,:),fson);

            % w:ch4 cospectrum
            [Cmmw(ii,:),~] = csd2(uvwt_win(:,3,ii),detrend(ch4_seg(bins(ii,1):bins(ii,2)),'linear'),...
                                  len,fson,hamming(len));
            Cmmw(ii,:) = real(Cmmw(ii,:));
            % integrate cospectrum for wch4_csp, wch4_csp_corr & wch4_csp_uM
            mc_flx.wch4_csp_ppb(zz(ii)) = sum(Cmmw(ii,:)*df);
            mc_flx.wch4_csp_corr_ppb(zz(ii)) = mc_flx.wch4_csp_ppb(zz(ii))/Ac;
            mc_flx.wch4_csp_corr_uM(zz(ii)) = mc_flx.wch4_csp_corr_ppb(zz(ii))*mM_m2_d_conv*1000; % uM/m2/d units
            [Cch4w(ii,:),~] = specsmoo(Cmmw(ii,:),fson);
        end
        
        %% save smoothed spectra/cospectra to spec structure
        if all(isnan(spec_sm(1).Fx)); spec_sm(1).Fx = Fx; end
        
        spec_sm(1).Pq(zz,:) = Pq;
        spec_sm(1).Pc(zz,:) = Pc;
        spec_sm(1).Pm(zz,:) = Pm;
        
        spec_sm(1).Cqw(zz,:) = Cqw;
        spec_sm(1).Cco2w(zz,:) = Cco2w;
        spec_sm(1).Cch4w(zz,:) = Cch4w;
        
        %% plot Picarro spectra & cospectra

        if plotit
            % wq cospectrum
            if exist('h7','var')
                clf(h7)
                set(0,'CurrentFigure', h7);
            else
                h7 = figure('position',[500,100,700,500],'Name','wq_csp');
            end
            wq = nanmean1(Cqw,1);
            semilogx(Fx,Fx.*wq,'b-',[1e-3 10 ],[0 0],'k-','linewidth',2);
            title([MM,'/',DD,'/',YYYY,' ',HH,'00 wq Cospectrum']); grid;
            ylabel('C_{wq}(f) x f'); xlabel('f (Hz)');
            ppath = fullfile(path_wqcsp,[YYYY,'_',MM,'_',DD,'_',HH,'00_wqCospec.png']);
            print(ppath,'-dpng');
            
            % wco2 cospectrum
            if exist('h8','var')
                clf(h8)
                set(0,'CurrentFigure', h8);
            else
                h8 = figure('position',[500,100,700,500],'Name','wco2_csp');
            end
            wco2 = nanmean1(Cco2w,1);
            semilogx(Fx,Fx.*wco2,'b-',[1e-3 10 ],[0 0],'k-','linewidth',2);
            title([MM,'/',DD,'/',YYYY,' ',HH,'00 wco2 Cospectrum']); grid;
            ylabel('C_{wco2}(f) x f'); xlabel('f (Hz)');
            ppath = fullfile(path_wco2csp,[YYYY,'_',MM,'_',DD,'_',HH,'00_wco2Cospec.png']);
            print(ppath,'-dpng');
            
            % wch4 cospectrum
            if exist('h9','var')
                clf(h9)
                set(0,'CurrentFigure', h9);
            else
                h9 = figure('position',[500,100,700,500],'Name','wch4_csp');
            end
            wch4 = nanmean1(Cch4w,1);
            semilogx(Fx,Fx.*wch4,'b-',[1e-3 10 ],[0 0],'k-','linewidth',2);
            title([MM,'/',DD,'/',YYYY,' ',HH,'00 wch4 Cospectrum']); grid;
            ylabel('C_{wch4}(f) x f'); xlabel('f (Hz)');
            ppath = fullfile(path_wch4csp,[YYYY,'_',MM,'_',DD,'_',HH,'00_wch4Cospec.png']);
            print(ppath,'-dpng');
            
            % q spectrum
            if exist('h10','var')
                clf(h10)
                set(0,'CurrentFigure', h10);
            else
                h10 = figure('position',[500,100,700,500],'Name','q_spc');
            end
            Sq = nanmean1(Pq,1);
            loglog(Fx(3:end),Fx(3:end).*Sq(3:end),'b-','linewidth',2); hold on;
            loglog(Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*Sq(35).*Fx(35),'r');
            hold off; ylabel('S_{qq} x f'); xlabel('f (Hz)');
            title([MM,'/',DD,'/',YYYY,' ',HH,'00 Spectrum of specific humidity, q']); grid;
            ppath = fullfile(path_qspc,[YYYY,'_',MM,'_',DD,'_',HH,'00_qSpec.png']);
            print(ppath,'-dpng');
            
            % co2 spectrum
            if exist('h11','var')
                clf(h11)
                set(0,'CurrentFigure', h11);
            else
                h11 = figure('position',[500,100,700,500],'Name','co2_spc');
            end
            Sco2 = nanmean1(Pc,1);
            loglog(Fx(3:end),Fx(3:end).*Sco2(3:end),'b-','linewidth',2); hold on;
            loglog(Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*Sco2(35).*Fx(35),'r');
            hold off; ylabel('S_{co2} x f'); xlabel('f (Hz)');
            title([MM,'/',DD,'/',YYYY,' ',HH,'00 Spectrum of carbon dioxide, CO_2']); grid;
            ppath = fullfile(path_co2spc,[YYYY,'_',MM,'_',DD,'_',HH,'00_co2Spec.png']);
            print(ppath,'-dpng');
            
            % ch4 spectrum
            if exist('h12','var')
                clf(h12)
                set(0,'CurrentFigure', h12);
            else
                h12 = figure('position',[500,100,700,500],'Name','ch4_spc');
            end
            Sch4 = nanmean1(Pm,1);
            loglog(Fx(3:end),Fx(3:end).*Sch4(3:end),'b-','linewidth',2); hold on;
            loglog(Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*Sch4(35).*Fx(35),'r');
            hold off; ylabel('S_{ch4} x f'); xlabel('f (Hz)');
            title([MM,'/',DD,'/',YYYY,' ',HH,'00 Spectrum of methane, CH_4']); grid;
            ppath = fullfile(path_ch4spc,[YYYY,'_',MM,'_',DD,'_',HH,'00_ch4Spec.png']);
            print(ppath,'-dpng');
        end
    end
    
    
    %% plot theoretical cospectrum and transfer functions for this hour
    if plotit && exist('Km','var')
        if exist('h15','var')
            clf(h15)
            set(0,'CurrentFigure', h15);
        else
            h15 = figure('position',[500,100,700,500],'Name','attn_csp');
        end
        semilogx(fx,Km,'k-','linewidth',2); hold on; grid;
        if exist('Kq','var'); semilogx(fx,Kq,'r-','linewidth',1); end
        if exist('Kcc','var'); semilogx(fx,Kcc,'b-','linewidth',1); end
        hold off; ylabel('f C_{wx}(f) / u_*X_*'); xlabel('Frequency, Hz');
        title([MM,'/',DD,'/',YYYY,' ',HH,'00 Estimated Cospectral Attenuation, Kaimal (Neutral, Scalar)']);
        subtitle('black = Kaimal; red = q; blue = CO_2/CH_4');
        ppath = fullfile(path_attn,[YYYY,'_',MM,'_',DD,'_',HH,'00_attnCospec.png']);
        print(ppath,'-dpng');
    end
    
    % plot inlet attenuation transfer functions
    if plotit && (exist('Hfq','var') || exist('Hfcc','var'))
        if exist('h16','var')
            clf(h16)
            set(0,'CurrentFigure', h16);
        else
            h16 = figure('position',[500,100,700,500],'Name','attn_transf');
        end
        if exist('Hfq','var'); semilogx(fx,sqrt(Hfq),'r-','linewidth',2); hold on; end
        if exist('Hfcc','var'); semilogx(fx,sqrt(Hfcc),'b-','linewidth',2); hold on; end
        hold off; grid; xlabel('Frequency, Hz'); ylabel('\surd H_f');
        title([MM,'/',DD,'/',YYYY,' ',HH,'00 Attenuation Transfer Functions']);
        subtitle('red = q; blue = CO_2/CH_4');
        ppath = fullfile(path_attn,[YYYY,'_',MM,'_',DD,'_',HH,'00_attnTransfunc.png']);
        print(ppath,'-dpng');
    end
    clear Km Kq Kcc Aq Ac Ao Hfq Hfcc
    
end
clear h1 h2 h3 h4 h5 h6 h7 h8 h9 h10 h11 h12 h13 h14 h15 h16
close all

%% save the mc_flx data table with additional columns

sfl = fullfile(save_path,'wins_met_gases_lags_flux_8192_v2.mat');
save(sfl,'mc_flx');

sfl = fullfile(save_path,'spectra_mc_8192_v2.mat');
save(sfl,'spec_sm','-v7.3');

%%............. END OF 8192 POINT FLUX CALCULATIONS...................%%

display('DONE')

