% b2_hourly_met_metcity_8192.m

%{
    Load flux windows (identical to ship run) and compute mean tower met
    for each flux window.  Further compute hourly means.
    
    Table mc_flx will have the following columns:

        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 for 1 hr)
        8     seg_end                   segment end point (36000 for 1 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
%}

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

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

path(pathdef);
path_prog = fullfile(root_data_path,'_matlab');
cd(fullfile(path_prog,'co2_ch4_o3_flux'));
addpath(genpath(path_prog));

save_path = fullfile(root_data_path,'_processed','metcity_flux_v3','flux_merged');
mkdir(save_path);
fig_path = fullfile(save_path,'plots'); % for saving plots
mkdir(fig_path);

met_path = fullfile(root_data_path,'_processed','metcity_tower_level2_v4');

%% load 8192 windows and set up for averaging
load(fullfile(save_path,'wins_8192_all.mat'));  % load flux windows table 'mc_flx'
vars = mc_flx.Properties.VariableNames;
www = table2array(mc_flx);

load(fullfile(met_path,'metcity_1min_met_all.mat'));  % load 1-min met city tower met table 'mc_met'
qc_mat = mc_met{:,94:185}; qc_mat = logical(qc_mat);  % qc variables as mask for bad data
data_mat = mc_met{:,2:93}; data_mat(qc_mat) = NaN;    % apply mask to NaN out bad points
data_mat = [mc_met{:,1},data_mat];                    % add datenum
met_vars = mc_met.Properties.VariableNames; met_vars = met_vars(1:93); % data variables
mc_met = array2table(data_mat,'VariableNames',met_vars); % rebuild table with filtered data

www = [www,NaN(length(www),92)];  % mc_flx array with added columns for met
vars = [vars,met_vars(2:93)];     % variable names for result array

flx_bins = [mc_flx.datenum,mc_flx.datenum + 8192/36000/24];  % overlapping window bins
hr_bins = fix(mc_flx.datenum(1)):1/24:ceil(mc_flx.datenum(end)); hr_bins = hr_bins(:);

%% average all met over flux windows - windows are overlapping, so this loop takes a long time
for ii=1:length(flx_bins)
    display(['window count ',sprintf('%05i',ii)]);
    kk = find(mc_met.datenum>=flx_bins(ii,1) & mc_met.datenum<flx_bins(ii,2));
    if ~isempty(kk)
        temp = mc_met{kk,2:93};  % extract table selection as array
        www(ii,9:100) = nanmean1(temp,1); % save column mean of selection to result array
    end
end

% add 10m true wind speed and direction
[wspd,wdir] = uv_to_sd(www(:,65),-www(:,69));
vars = [vars,{'wspd_10m','wdir_10m'}];
www(:,101:102) = [wspd,wdir];

% save table of flux window data
clear mc_flx
mc_flx = array2table(www,'VariableNames',vars);
sfl = fullfile(save_path,'wins_met_8192_all.mat');
save(sfl,'mc_flx');

%% compute hourly averages

xxx = interval_avg(mc_met{:,1},mc_met{:,2:93},hr_bins);

% add 10m true wind speed and direction
[wspd,wdir] = uv_to_sd(xxx(:,58),-xxx(:,62));
xxx(:,94:95) = [wspd,wdir];
% add datevec columns
temp = datevec(xxx(:,1));
xxx = [xxx(:,1),temp(:,1:4),xxx(:,2:end)];
% save table of hourly met
hr_vars = [{'datenum','yyyy','mm','dd','hr'},met_vars(2:93),{'wspd_10m','wdir_10m'}];
mc_flx_hr = array2table(xxx,'VariableNames',hr_vars);
sfl = fullfile(save_path,'mc_flx_hr.mat');
save(sfl,'mc_flx_hr');

