% a10_build_OOT_datafiles_ship.m

%{
    Merge synchronized sonic and gas analyzer data into hourly files
    for the Ogive Optimization Toolbox routines.  Output files are matlab
    tables.
    
    table columns (variables) are are:
    DATA = [datenum,u,v,w,T,Pmb,co2_ppm,ch4_ppm,o3_ppm,dms_ppm] 
    where: 
      - time is the matlab-native datenum format
      - uvwT are three-dimensional rotated wind, and temperature output
        from the Metek USA-1 sonic
      - Pmb is pressure from the Licor 7500
      - co2_ppm, ch4_ppm, o3_ppb, dms_ppm are dry air umol/mol gas
        concentrations, lag adjusted using output from picarro and APIMS
        'N2 puff' lag analysis.Gas concentration data during the hourly
        N2 puff is set to the mean value, as are all NaNs.
      - Gas concentrations are all ppm because the smallest scaling factor
        in OOT is 10^-6 (micro).
%}

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

%% adjust these parameters as if necesssary, then run
% anemometer, APIMS flow and dms isotopomer abundance constants
fson = 10;              % sampling frequency, sonic and all gases
orifice_flow = 0.3;     % LPM flow into vacuum chamber
a63s = 0.0081;          % labeled DMS std abundance mass 63
a66s = 0.86124;         % labeled DMS std abundance mass 66
a63 = 0.9298;           % natural DMS abundance mass 63
a66 = 0.0004;           % natural DMS abundance mass 66 (can ignore this)
% output table variable names
vars = {'datenum','u','v','w','Ts','Pmb','co2_ppm','ch4_ppm',...
        'o3_ppm','dms_ppm'};

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

wind_path = fullfile(root_data_path,'sonic_anemometer_3D_ship_bow');
picarro_path = fullfile(root_data_path,'crds_ghg_ship_bow');
ozone_path = fullfile(root_data_path,'ozone_fast_ship_bow');
dms_path = fullfile(root_data_path,'dms_apims_ship_bow','SIM');
licor_path = fullfile(root_data_path,'licor_7500ds_ship_bow');
inlet_flow_path = fullfile(root_data_path,'CO2_Alicat');
o3_flow_path = fullfile(ozone_path,'FastO3#2_control_files');

save_path = fullfile(root_data_path,'_processed','ship_flux_V3',...
                     'merged_10Hz_turb');
mkdir(save_path);

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

%% load flx table with lag times for CO2, CH4, O3 and DMS
% for CO2,CH4 and O3
load(fullfile(root_data_path,'_processed','ship_flux_V3','flux_8192',...
              'wins_gases_lags_8192_all.mat'),'flx');
% for DMS
load(fullfile(root_data_path,'_processed','dms_flux_v2','flux_8192',...
              'dms_wins_lags_8192_all.mat'),'dms_flx');

% load fo3 calibration coefs for all legs
load(fullfile(merge_path,'fo3_cal_coefs_all_legs.mat'));

% load TEI ozone for entire cruise - table slo3_ship: col1 is datenum, col 2 is O3_ppb
load(fullfile(root_data_path,'_processed','ship_flux_V3','o3slow_ship_1min_all_legs.mat'));

%% set up hourly timestamp table for entire cruise
times_hr = flx(1:7:end,1:6);
times_hr.ddd = md2yd_vn(times_hr.yyyy,times_hr.mm,times_hr.dd);
N = height(times_hr);

for ii = 1:N
    YYYY = sprintf('%04i',times_hr.yyyy(ii));
    YY = YYYY(3:4);
    MM = sprintf('%02i',times_hr.mm(ii));
    DD = sprintf('%02i',times_hr.dd(ii));
    DDD= sprintf('%03i',times_hr.ddd(ii));
    HH = sprintf('%02i',times_hr.hh(ii));
    wind_fldr = fullfile(wind_path,[YY,DDD]);
    co2_fldr = fullfile(picarro_path,[YY,DDD]);
    o3_fldr = fullfile(ozone_path,[YY,DDD]);
    lic_fldr = fullfile(licor_path,[YY,DDD]);
    dms_fldr = fullfile(dms_path,[YY,DDD]);
    co2flow_fldr = fullfile(inlet_flow_path,[YY,DDD]);
    
    wind_dfl = fullfile(wind_fldr,['son0',YY,DDD,HH,'_raw.txt']);
    co2_dfl = fullfile(co2_fldr,['pic0',YY,DDD,HH,'_raw.txt']);
    o3_dfl = fullfile(o3_fldr,['fo30',YY,DDD,HH,'_raw.txt']);
    dms_dfl = fullfile(dms_fldr,['DMS_',YYYY,'_',DDD,'_',HH,'00.txt']);
    lic_dfl = fullfile(lic_fldr,['lic0',YY,DDD,HH,'_raw.txt']);
    co2flow_dfl = fullfile(co2flow_fldr,['ali0',YY,DDD,HH,'_raw.txt']);

    % load inlet flow rate data for co2/dms and o3 systems
    co2flow_dir = dir(fullfile(o3_flow_path,['control_',YY,'-',MM,'-',DD,'*.lvm']));
    % build cell array of file paths to O3 lvm files for this day
    o3flow_flist = {};
    for ff = 1:length(co2flow_dir)
        o3flow_flist{ff} = fullfile(co2flow_dir(ff).folder,co2flow_dir(ff).name);
    end
    
    % select DMS standard concentration based on date
    if times_hr.yyyy(ii) == 2020 && times_hr.ddd(ii) + times_hr.hh(ii)/24 >= (241+4/24)
        conc_d3_std_ppb = 532.2; % EB0125416 in use after 8/28/20 0400
    elseif times_hr.yyyy(ii) == 2020 && ...
           times_hr.ddd(ii) + times_hr.hh(ii)/24 >= (218+15/24) && ...
           times_hr.ddd(ii) + times_hr.hh(ii)/24 < (241+4/24)
        conc_d3_std_ppb = (542.9+547.6+548.7)/3; % CA05700 in use after 8/5/20 1500 UTC
    else
        conc_d3_std_ppb = (563.8+576.3)/2; % CA05698 used prior to 8/5/20 1500
    end

    if exist(wind_dfl,'file') % execute only if there's wind data for this hour...
        %% load 10Hz wind and gas data for this hour

        % these are dummy variables required for son reader, but not used
        % in the output result we use here
        hed0 = ones(1440,2); hed0(:,1) = times_hr.ddd(ii); hed0(:,2) = 0:1:1439;
        Nbins = 7; bins = zeros(Nbins,2);
        bins(:,1) = 600:4096:34800-8192; bins(:,2) = bins(:,1) + 8191;
        len = 8192;
        
        [son,~,~,~,~,~,~] = read_sonic3D_ship_mosaic_2019_bins(wind_dfl,...
                            hed0,times_hr.ddd(ii),times_hr.hh(ii),bins,len,0);
        pic = read_picarro_mosaic_2019(co2_dfl,times_hr.yyyy(ii),times_hr.ddd(ii),times_hr.hh(ii));
        fo3_raw = read_fo3_cnts_mosaic_2019(o3_dfl,times_hr.ddd(ii),times_hr.hh(ii));
        lic = read_licor_mosaic_2019(lic_dfl,times_hr.ddd(ii),times_hr.hh(ii));
        dms = read_dms_mosaic_2019(dms_dfl,times_hr.ddd(ii),times_hr.hh(ii));
        
        % read 1-min mean inlet flow rates for co2/dms and o3 systems
        co2flow = read_co2_alicat_mosaic_2019(co2flow_dfl,times_hr.ddd(ii),times_hr.hh(ii));
        co2flow = [co2flow(:,1),co2flow(:,5)];  % select inlet flow slpm
        o3flow = read_fo3_control_mosaic_2019(o3flow_flist,times_hr.ddd(ii)); % reads 1 day of data
        zz = o3flow(:,1)>=co2flow(1,1) & o3flow(:,1)<=co2flow(end,1); % select minutes for this hour
        o3flow = [o3flow(zz,1),o3flow(zz,11)]; % select inlet flow slpm this hour

        %% process licor for pressure data
        Pmb = lic(:,6)*10; % kPa to hPa (mb)

        %% process CO2/CH4
        % load hourly co2/ch4 parameters from flx table
        zz = find(flx{:,3}==times_hr.yyyy(ii) & ...
                  flx{:,4}==times_hr.mm(ii) & ...
                  flx{:,5}==times_hr.dd(ii) & ...
                  flx{:,6}==times_hr.hh(ii),1,'first');
        if ~isempty(zz)
            lag_co2 = flx.co2_lag(zz);
            puff_endpt = flx.puff_endpt(zz); if puff_endpt==1; puff_endpt=2; end

            % Find puff - might occurr at start or end of hour.
            % Puff indices in flx indicate the puff nearest top of hour,
            % so indices may be negative if nearest puff was just prior
            % to the start of the hour (i.e. in the prior hour data file).
            %
            % If puff index is positive, then it's at start of this hour, so no problem.
            % If this hour puff index is negative, then the true puff index for 
            % this hour's data will be saved in the next hour in the mc_flx table.
            if flx.puff_stpt(zz)>0                % puff starts after top of hour
                puff_stpt = flx.puff_stpt(zz);
            elseif ii<N && flx.puff_stpt(zz+1)<0  % puff starts before end of hour
                puff_stpt = flx.puff_stpt(zz+1) + 36000;
            else
                puff_stpt = NaN;
            end
            if flx.puff_endpt(zz)>0               % puff ends after top of hour
                puff_endpt = flx.puff_endpt(zz);
            elseif ii<N && flx.puff_endpt(zz+1)<0 % puff ends before end of hour
                puff_endpt = flx.puff_endpt(zz+1) + 36000;
            else
                puff_endpt = NaN;
            end
        else
            lag_co2 = NaN;
            puff_stpt = NaN;
            puff_endpt = NaN;
        end

        % process and save cos/ch4 if lag is good and inlet flow is >25 slpm for
        % the entire hour
        if isfinite(lag_co2) && lag_co2>0 && all(co2flow(:,2)>25)
            co2 = pic(:,9);
            ch4 = pic(:,7);

            % Remove puff periods from co2/ch4 conc data
            % when puff_stpt and puff_endpt are not NaN.
            if (puff_stpt < puff_endpt) % entire puff is within this hour
                co2(puff_stpt:puff_endpt) = nanmean1(co2);
                ch4(puff_stpt:puff_endpt) = nanmean1(ch4);
            elseif (puff_stpt >puff_endpt) % puff straddles top of hour
                co2(1:puff_endpt) = nanmean1(co2);
                ch4(1:puff_endpt) = nanmean1(ch4);
                co2(puff_stpt:end) = nanmean1(co2);
                ch4(puff_stpt:end) = nanmean1(ch4);
            end

            % lag-shift co2/ch4 data
            co2_shift = co2;
            co2_shift(1:end-lag_co2+1) = co2(lag_co2:end);
            co2_shift = despike(co2_shift);
            ch4_shift = ch4;
            ch4_shift(1:end-lag_co2+1) = ch4(lag_co2:end);
            ch4_shift = despike(ch4_shift);
        else
            % don't use co2/ch4 this hour if  shift or flow is bad
            co2_shift = NaN(36000,1);
            ch4_shift = NaN(36000,1);
        end

        %% process O3
        if ~isempty(zz)
            lag_o3 = flx.o3_lag(zz);
        else
            lag_o3 = NaN;
        end
        % process and save o3 if lag is good and inlet flow is >25 slpm for
        % the entire hour
        if isfinite(lag_o3) && lag_o3>0 && all(o3flow(isfinite(o3flow(:,2)),2)>25)
            % recompute O3 from counts using cal coefficient for cruise leg
            temp = fo3_raw(:,2)*fo3_coef(flx.leg(zz),1) + fo3_coef(flx.leg(zz),2);
            % adjust mean concentration to the TEI value
            oo = slo3_ship.datenum>=flx.datenum(zz) & slo3_ship.datenum<flx.datenum(zz)+1/24;
            slo3 = nanmean1(slo3_ship.O3_ppb(oo)); % mean slow ozone conc for this hour
            temp = temp*(slo3/nanmean1(fo3_raw(1000:end-1000,3)));
            o3 = temp;  % adjusted raw fast ozone cocnentration

            % Remove puff periods from o3 conc data
            % when puff_stpt and puff_endpt are not NaN.
            if (puff_stpt < puff_endpt) % entire puff is within this hour
                o3(puff_stpt:puff_endpt) = nanmean1(o3);
            elseif (puff_stpt >puff_endpt) % puff straddles top of hour
                o3(1:puff_endpt) = nanmean1(o3);
                o3(puff_stpt:end) = nanmean1(o3);
            end

            % lag-shift o3 data
            o3_shift = o3;
            o3_shift(1:end-lag_o3+1) = o3(lag_o3:end);
            o3_shift = despike(o3_shift);
        else
            % don't use o3 this hour if shift or flow is bad
            o3_shift = NaN(36000,1);
        end

        %% process DMS
        % load hourly dms parameters from dms flx table
        zz = find(dms_flx{:,3}==times_hr.yyyy(ii) & ...
                  dms_flx{:,4}==times_hr.mm(ii) & ...
                  dms_flx{:,5}==times_hr.dd(ii) & ...
                  dms_flx{:,6}==times_hr.hh(ii),1,'first');
        if ~isempty(zz)
            lag_dms = dms_flx.lag(zz);
            tau_dms = dms_flx.tau(zz);
            m66b = dms_flx.m66b(zz);
            m63b = dms_flx.m63b(zz);
        else
            lag_dms = NaN;
            m66b = NaN;
            m63b = NaN;
        end
        % process and save dms if lag is good and inlet flow is >25 slpm for
        % the entire hour
        if isfinite(lag_dms) && lag_dms>0 && all(co2flow(:,2)>25)
            % compute 'dry' 10Hz dms concentration in pptv
            cal_flow_lpm = calFlow(dms(:,6));
            api_flow_lpm = apiFlow(dms(:,7));
            man_flow_lpm = manFlow(dms(:,8));
            total_flow = man_flow_lpm + api_flow_lpm + orifice_flow;
            std_conc_ppt = conc_d3_std_ppb*1000.*cal_flow_lpm./total_flow;
            m63 = dms(:,3) - m63b*0.80; % background cnts ~ 20% lower when d3-DMS is added
            m66 = dms(:,2) - m66b*0.80;
            R = m63./m66;
            dms_ppt = std_conc_ppt.*((R*a66s - a63s)./(a63));
            % remove blank and puff periods from dms conc data
            dms(dms(:,5)>4.5,5) = NaN; dms(:,5) = replace_NaN_nearest_neighbor(dms(:,5));
            bad1 = find(dms(:,5)>4.05,1,'first');
            bad2 = find(dms(:,5)>4.05,1,'last'); bad2 = bad2+100;
            dms_ppt_clean = dms_ppt; dms_ppt_clean(bad1:bad2) = NaN;
            % lag-shift dms data
            dms_shift = dms_ppt_clean;
            dms_shift(1:end-lag_dms+1) = dms_ppt_clean(lag_dms:end);
            dms_shift = despike(dms_shift);
        else
            % don't use dms this hour if lag shift or flow bad
            dms_shift = NaN(36000,1);
        end

        %% format output table and save as .mat
        dnum = datenum(times_hr.yyyy(ii)-1,12,31) + son(:,1);
        % using raw u,v,w in ship frame (relative) coordinates, no rotation
        % converting O3 and DMS to ppm
        outmat = [dnum,son(:,2:5),Pmb,co2_shift,ch4_shift,o3_shift/1000,dms_shift/1e6];
        flx_10Hz = array2table(outmat,'VariableNames',vars);
        fname = [YYYY,'_',MM,'_',DD,'_',HH,'_MOSAiC_10Hz_winds_gases_ship.mat'];
        save(fullfile(save_path,fname),'flx_10Hz');
        disp(['Done with file ',fname]);
    end
end


