%% Read in table properly
fpath = 'C:\Users\bwerb\Documents\CUGNROMS\IOOS glider data 2021';
files = dir(fullfile(fpath, '*.csv'));   % filter to CSVs, or use '*.*' for all
% i = 17;
for i = 1
    fprintf('file %s', string(i));
    fname = files(i).name;
    
    if contains(fname, {'ce','ocg'})
        opts = detectImportOptions(fname);
        opts = setvaropts(opts, opts.VariableNames{1}, 'Type', 'char');   % time (UTC)
        opts = setvaropts(opts, opts.VariableNames{2}, 'Type', 'double'); % depth (m)
        opts = setvaropts(opts, opts.VariableNames{3}, 'Type', 'double'); % latitude
        opts = setvaropts(opts, opts.VariableNames{4}, 'Type', 'double'); % longitude
        opts = setvaropts(opts, opts.VariableNames{5}, 'Type', 'double'); % temperature
        opts = setvaropts(opts, opts.VariableNames{6}, 'Type', 'double'); % salinity
        opts = setvaropts(opts, opts.VariableNames{7}, 'Type', 'double'); % dissolved_oxygen
        opts = setvaropts(opts, opts.VariableNames{8}, 'Type', 'double'); % density
        % opts.SelectedVariableNames = opts.VariableNames(2:end);
        opts.datalines = [3 Inf]; % skip first 2 lines of metadata
        opts.VariableNames = {'time_UTC_','depth','latn','lone','tempc','psal','do','sigma'};
        T_RAW = readtable(fname, opts);
        T_RAW.time_UTC_ = datetime(T_RAW.time_UTC_, 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss''Z''', 'TimeZone', 'UTC');    
        T_RAW.mon_day_yr = T_RAW.time_UTC_;
        T_RAW.mon_day_yr.Format = 'MM/dd/uuuu';
        T_RAW.hh_mm = T_RAW.time_UTC_;
        T_RAW.hh_mm.Format = 'HH:mm';
    
        T = table();
        T.WMO_ID = repmat({fname(4:6)}, height(T_RAW), 1);
        T.Prof_num = NaN(height(T_RAW), 1);
        T.mon_day_yr = T_RAW.mon_day_yr;
        T.hh_mm      = T_RAW.hh_mm;
        T.Lon_E = T_RAW.lone;
        T.Lat_N = T_RAW.latn;
        T.Position_QF = NaN(height(T_RAW), 1);
        T.Pressure_dbar = gsw_p_from_z(-T_RAW.depth, T_RAW.latn);
        T.Pressure_dbar_QF = NaN(height(T_RAW), 1);
        T.Temperature_C = T_RAW.tempc;
        T.Temperature_C_QF = NaN(height(T_RAW), 1);
        T.Salinity_pss = T_RAW.psal;
        T.Salinity_pss_QF = NaN(height(T_RAW), 1);
        T.Sigma_theta_kg_m3 = T_RAW.sigma;
        T.Sigma_theta_kg_m3_QF = NaN(height(T_RAW), 1);
        T.Depth_m = T_RAW.depth;
        T.Depth_m_QC = NaN(height(T_RAW), 1);
        T.Oxygen_umol_kg = T_RAW.do;
        T.Oxygen_umol_kg_QF = NaN(height(T_RAW), 1);
        T.OxygenSat = NaN(height(T_RAW), 1);
        T.OxygenSat_QF = NaN(height(T_RAW), 1);
        T.Nitrate_umol_kg = NaN(height(T_RAW), 1);
        T.Nitrate_umol_kg_QF = NaN(height(T_RAW), 1);
        T.pHinsitu_Total = NaN(height(T_RAW), 1);
        T.pHinsitu_Total_QF = NaN(height(T_RAW), 1);
        T.TALK_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.DIC_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.pH_ALGORITHM_insitu_total = NaN(height(T_RAW), 1);
        T.NO3_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.PO4_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.SILICATE_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.DIC_pHTA_umol_kg = NaN(height(T_RAW), 1);
        T.PCO2_pHTA_uatm = NaN(height(T_RAW), 1);
    
    elseif contains(fname, 'gp_')
        % gp-specific readtable setup
        opts = detectImportOptions(fname);
        opts = setvaropts(opts, opts.VariableNames{1}, 'Type', 'char');   % time (UTC)
        opts = setvaropts(opts, opts.VariableNames{2}, 'Type', 'double'); % depth (m)
        opts = setvaropts(opts, opts.VariableNames{3}, 'Type', 'double'); % latitude
        opts = setvaropts(opts, opts.VariableNames{4}, 'Type', 'double'); % longitude
        opts = setvaropts(opts, opts.VariableNames{5}, 'Type', 'double'); % temperature
        opts = setvaropts(opts, opts.VariableNames{6}, 'Type', 'double'); % salinity
        opts = setvaropts(opts, opts.VariableNames{7}, 'Type', 'double'); % dissolved_oxygen
        opts = setvaropts(opts, opts.VariableNames{8}, 'Type', 'double'); % density
        opts.datalines = [3 Inf]; % skip first 2 lines of metadata
        % opts.SelectedVariableNames = opts.VariableNames(2:end);
        opts.VariableNames = {'time_UTC_','depth','latn','lone','tempc','psal','do','sigma'};
        T_RAW = readtable(fname, opts);
        T_RAW.time_UTC_ = datetime(T_RAW.time_UTC_, 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss''Z''', 'TimeZone', 'UTC');    
        T_RAW.mon_day_yr = T_RAW.time_UTC_;
        T_RAW.mon_day_yr.Format = 'MM/dd/uuuu';
        T_RAW.hh_mm = T_RAW.time_UTC_;
        T_RAW.hh_mm.Format = 'HH:mm';
    
        T = table();
        T.WMO_ID = repmat({fname(4:6)}, height(T_RAW), 1);
        T.Prof_num = NaN(height(T_RAW), 1);
        T.mon_day_yr = T_RAW.mon_day_yr;
        T.hh_mm      = T_RAW.hh_mm;
        T.Lon_E = T_RAW.lone;
        T.Lat_N = T_RAW.latn;
        T.Position_QF = NaN(height(T_RAW), 1);
        T.Pressure_dbar = gsw_p_from_z(-T_RAW.depth, T_RAW.latn);
        T.Pressure_dbar_QF = NaN(height(T_RAW), 1);
        T.Temperature_C = T_RAW.tempc;
        T.Temperature_C_QF = NaN(height(T_RAW), 1);
        T.Salinity_pss = T_RAW.psal;
        T.Salinity_pss_QF = NaN(height(T_RAW), 1);
        T.Sigma_theta_kg_m3 = T_RAW.sigma;
        T.Sigma_theta_kg_m3_QF = NaN(height(T_RAW), 1);
        T.Depth_m = T_RAW.depth;
        T.Depth_m_QC = NaN(height(T_RAW), 1);
        T.Oxygen_umol_kg = T_RAW.do;
        T.Oxygen_umol_kg_QF = NaN(height(T_RAW), 1);
        T.OxygenSat = NaN(height(T_RAW), 1);
        T.OxygenSat_QF = NaN(height(T_RAW), 1);
        T.Nitrate_umol_kg = NaN(height(T_RAW), 1);
        T.Nitrate_umol_kg_QF = NaN(height(T_RAW), 1);
        T.pHinsitu_Total = NaN(height(T_RAW), 1);
        T.pHinsitu_Total_QF = NaN(height(T_RAW), 1);
        T.TALK_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.DIC_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.pH_ALGORITHM_insitu_total = NaN(height(T_RAW), 1);
        T.NO3_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.PO4_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.SILICATE_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.DIC_pHTA_umol_kg = NaN(height(T_RAW), 1);
        T.PCO2_pHTA_uatm = NaN(height(T_RAW), 1);
    
    elseif contains(fname, 'osu')
        % osu-specific readtable setup
        opts = detectImportOptions(fname);
        opts = setvaropts(opts, opts.VariableNames{1}, 'Type', 'char');   % time (UTC)
        opts = setvaropts(opts, opts.VariableNames{2}, 'Type', 'double'); % depth (m)
        opts = setvaropts(opts, opts.VariableNames{3}, 'Type', 'double'); % latitude
        opts = setvaropts(opts, opts.VariableNames{4}, 'Type', 'double'); % longitude
        opts = setvaropts(opts, opts.VariableNames{5}, 'Type', 'double'); % temperature
        opts = setvaropts(opts, opts.VariableNames{6}, 'Type', 'double'); % salinity
        opts = setvaropts(opts, opts.VariableNames{7}, 'Type', 'double'); % dissolved_oxygen
        opts = setvaropts(opts, opts.VariableNames{8}, 'Type', 'double'); % density
        % opts.SelectedVariableNames = opts.VariableNames(2:end);
        opts.datalines = [3 Inf]; % skip first 2 lines of metadata
        opts.VariableNames = {'time_UTC_','depth','latn','lone','tempc','psal','do','sigma'};
        T_RAW = readtable(fname, opts);
        T_RAW.time_UTC_ = datetime(T_RAW.time_UTC_, 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss''Z''', 'TimeZone', 'UTC');    
        T_RAW.mon_day_yr = T_RAW.time_UTC_;
        T_RAW.mon_day_yr.Format = 'MM/dd/uuuu';
        T_RAW.hh_mm = T_RAW.time_UTC_;
        T_RAW.hh_mm.Format = 'HH:mm';
    
        T = table();
        T.WMO_ID = repmat({fname(4:6)}, height(T_RAW), 1);
        T.Prof_num = NaN(height(T_RAW), 1);
        T.mon_day_yr = T_RAW.mon_day_yr;
        T.hh_mm      = T_RAW.hh_mm;
        T.Lon_E = T_RAW.lone;
        T.Lat_N = T_RAW.latn;
        T.Position_QF = NaN(height(T_RAW), 1);
        T.Pressure_dbar = gsw_p_from_z(-T_RAW.depth, T_RAW.latn);
        T.Pressure_dbar_QF = NaN(height(T_RAW), 1);
        T.Temperature_C = T_RAW.tempc;
        T.Temperature_C_QF = NaN(height(T_RAW), 1);
        T.Salinity_pss = T_RAW.psal;
        T.Salinity_pss_QF = NaN(height(T_RAW), 1);
        T.Sigma_theta_kg_m3 = T_RAW.sigma;
        T.Sigma_theta_kg_m3_QF = NaN(height(T_RAW), 1);
        T.Depth_m = T_RAW.depth;
        T.Depth_m_QC = NaN(height(T_RAW), 1);
        T.Oxygen_umol_kg = T_RAW.do;
        T.Oxygen_umol_kg_QF = NaN(height(T_RAW), 1);
        T.OxygenSat = NaN(height(T_RAW), 1);
        T.OxygenSat_QF = NaN(height(T_RAW), 1);
        T.Nitrate_umol_kg = NaN(height(T_RAW), 1);
        T.Nitrate_umol_kg_QF = NaN(height(T_RAW), 1);
        T.pHinsitu_Total = NaN(height(T_RAW), 1);
        T.pHinsitu_Total_QF = NaN(height(T_RAW), 1);
        T.TALK_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.DIC_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.pH_ALGORITHM_insitu_total = NaN(height(T_RAW), 1);
        T.NO3_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.PO4_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.SILICATE_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.DIC_pHTA_umol_kg = NaN(height(T_RAW), 1);
        T.PCO2_pHTA_uatm = NaN(height(T_RAW), 1);
    
    elseif contains(fname, 'UW')
        % A. Quality control values
        % https://iop.apl.washington.edu/Seaglider_Quality_Control_Manual.html
        % These are the available quality control names and their numeric equivalents. They are taken from Argo, 2010 with the addition of QC_UNSAMPLED. Not all values are currently used.
        % 
        % QC_NO_CHANGE	0 - No QC was performed
        % QC_GOOD	1 - Value is ok
        % QC_PROBABLY_GOOD	2 - Value is likely good
        % QC_PROBABLY_BAD	3 - Potentially correctable
        % QC_BAD	4 - Untrustworthy and uncorrectable
        % QC_CHANGED	5 - Explicit manual change
        % QC_UNSAMPLED	6 - Explicitly not sampled (vs. expected but QC_MISSING)
        % QC_INTERPOLATED	8 - Interpolated value
        % QC_MISSING	9 - Value missing; instrument timed out
        opts = detectImportOptions(fname);
        opts = setvaropts(opts, opts.VariableNames{1}, 'Type', 'char');   % time (UTC)
        opts = setvaropts(opts, opts.VariableNames{2}, 'Type', 'double'); % depth (m)
        opts = setvaropts(opts, opts.VariableNames{3}, 'Type', 'double'); % latitude
        opts = setvaropts(opts, opts.VariableNames{4}, 'Type', 'double'); % longitude
        opts = setvaropts(opts, opts.VariableNames{5}, 'Type', 'double'); % temperature
        opts = setvaropts(opts, opts.VariableNames{6}, 'Type', 'double'); % salinity
        opts = setvaropts(opts, opts.VariableNames{7}, 'Type', 'double'); % oxygen
        opts = setvaropts(opts, opts.VariableNames{8}, 'Type', 'double'); % density
        % opts.SelectedVariableNames = opts.VariableNames(2:end);
        opts.datalines = [3 Inf]; % skip first 2 lines of metadata
        opts.VariableNames = {'time_UTC_','depth','latn','lone','tempc','psal','do','sigma'};
        T_RAW = readtable(fname, opts);
        T_RAW.time_UTC_ = datetime(T_RAW.time_UTC_, 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss''Z''', 'TimeZone', 'UTC');    
        T_RAW.mon_day_yr = T_RAW.time_UTC_;
        T_RAW.mon_day_yr.Format = 'MM/dd/uuuu';
        T_RAW.hh_mm = T_RAW.time_UTC_;
        T_RAW.hh_mm.Format = 'HH:mm';
        QC_idx = find(T_RAW.temp_QC == 1 & T_RAW.psal_QC == 1 & T_RAW.do_QC == 1);
        T_RAW = T_RAW(QC_idx,:);
        
        T = table();
        T.WMO_ID = repmat({fname(4:6)}, height(T_RAW), 1);
        T.Prof_num = NaN(height(T_RAW), 1);
        T.mon_day_yr = T_RAW.mon_day_yr;
        T.hh_mm      = T_RAW.hh_mm;
        T.Lon_E = T_RAW.lone;
        T.Lat_N = T_RAW.latn;
        T.Position_QF = NaN(height(T_RAW), 1);
        T.Pressure_dbar = gsw_p_from_z(-T_RAW.depth, T_RAW.latn);
        T.Pressure_dbar_QF = NaN(height(T_RAW), 1);
        T.Temperature_C = T_RAW.tempc;
        T.Temperature_C_QF = T_RAW.temp_QC;
        T.Salinity_pss = T_RAW.psal;
        T.Salinity_pss_QF = T_RAW.psal_QC;
        T.Sigma_theta_kg_m3 = T_RAW.sigma;
        T.Sigma_theta_kg_m3_QF = NaN(height(T_RAW), 1);
        T.Depth_m = T_RAW.depth;
        T.Depth_m_QC = NaN(height(T_RAW), 1);
        T.Oxygen_umol_kg = T_RAW.do;
        T.Oxygen_umol_kg_QF = T_RAW.do_QC;
        T.OxygenSat = NaN(height(T_RAW), 1);
        T.OxygenSat_QF = NaN(height(T_RAW), 1);
        T.Nitrate_umol_kg = NaN(height(T_RAW), 1);
        T.Nitrate_umol_kg_QF = NaN(height(T_RAW), 1);
        T.pHinsitu_Total = NaN(height(T_RAW), 1);
        T.pHinsitu_Total_QF = NaN(height(T_RAW), 1);
        T.TALK_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.DIC_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.pH_ALGORITHM_insitu_total = NaN(height(T_RAW), 1);
        T.NO3_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.PO4_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.SILICATE_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.DIC_pHTA_umol_kg = NaN(height(T_RAW), 1);
        T.PCO2_pHTA_uatm = NaN(height(T_RAW), 1);
    
    elseif contains(fname, 'dfo')
        opts = detectImportOptions(fname);
        opts = setvaropts(opts, opts.VariableNames{1}, 'Type', 'char');   % time (UTC)
        opts = setvaropts(opts, opts.VariableNames{2}, 'Type', 'double'); % depth (m)
        opts = setvaropts(opts, opts.VariableNames{3}, 'Type', 'double'); % latitude
        opts = setvaropts(opts, opts.VariableNames{4}, 'Type', 'double'); % longitude
        opts = setvaropts(opts, opts.VariableNames{5}, 'Type', 'double'); % temperature
        opts = setvaropts(opts, opts.VariableNames{6}, 'Type', 'double'); % salinity
        opts = setvaropts(opts, opts.VariableNames{7}, 'Type', 'double'); % oxygen_concentration ml/l need to convert to do umol/kg!!!!!
        opts = setvaropts(opts, opts.VariableNames{8}, 'Type', 'double'); % density
        % opts.SelectedVariableNames = opts.VariableNames(2:end);
        opts.datalines = [3 Inf]; % skip first 2 lines of metadata
        opts.VariableNames = {'time_UTC_','depth','latn','lone','tempc','psal','do','sigma'};
        T_RAW = readtable(fname, opts);
        T_RAW.time_UTC_ = datetime(T_RAW.time_UTC_, 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss''Z''', 'TimeZone', 'UTC');    
        T_RAW.mon_day_yr = T_RAW.time_UTC_;
        T_RAW.mon_day_yr.Format = 'MM/dd/uuuu';
        T_RAW.hh_mm = T_RAW.time_UTC_;
        T_RAW.hh_mm.Format = 'HH:mm';
        QC_idx = find(T_RAW.temp_QC == 1 & T_RAW.psal_QC == 1 & T_RAW.do_QC == 1);
        T_RAW = T_RAW(QC_idx,:);
        
        T = table();
        T.WMO_ID = repmat({fname(10:12)}, height(T_RAW), 1);
        T.Prof_num = NaN(height(T_RAW), 1);
        T.mon_day_yr = T_RAW.mon_day_yr;
        T.hh_mm      = T_RAW.hh_mm;
        T.Lon_E = T_RAW.lone;
        T.Lat_N = T_RAW.latn;
        T.Position_QF = NaN(height(T_RAW), 1);
        T.Pressure_dbar = gsw_p_from_z(-T_RAW.depth, T_RAW.latn);
        T.Pressure_dbar_QF = NaN(height(T_RAW), 1);
        T.Temperature_C = T_RAW.tempc;
        T.Temperature_C_QF = T_RAW.temp_QC;
        T.Salinity_pss = T_RAW.psal;
        T.Salinity_pss_QF = T_RAW.psal_QC;
        T.Sigma_theta_kg_m3 = T_RAW.sigma;
        T.Sigma_theta_kg_m3_QF = NaN(height(T_RAW), 1);
        T.Depth_m = T_RAW.depth;
        T.Depth_m_QC = NaN(height(T_RAW), 1);
        T.Oxygen_umol_kg = T_RAW.do;
        T.Oxygen_umol_kg_QF = T_RAW.do_QC;
        T.OxygenSat = NaN(height(T_RAW), 1);
        T.OxygenSat_QF = NaN(height(T_RAW), 1);
        T.Nitrate_umol_kg = NaN(height(T_RAW), 1);
        T.Nitrate_umol_kg_QF = NaN(height(T_RAW), 1);
        T.pHinsitu_Total = NaN(height(T_RAW), 1);
        T.pHinsitu_Total_QF = NaN(height(T_RAW), 1);
        T.TALK_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.DIC_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.pH_ALGORITHM_insitu_total = NaN(height(T_RAW), 1);
        T.NO3_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.PO4_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.SILICATE_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.DIC_pHTA_umol_kg = NaN(height(T_RAW), 1);
        T.PCO2_pHTA_uatm = NaN(height(T_RAW), 1);

    elseif contains(fname, 'sg')
        % gp-specific readtable setup
        opts = detectImportOptions(fname);
        opts = setvaropts(opts, opts.VariableNames{1}, 'Type', 'char');   % time (UTC)
        opts = setvaropts(opts, opts.VariableNames{2}, 'Type', 'double'); % depth (m)
        opts = setvaropts(opts, opts.VariableNames{3}, 'Type', 'double'); % latitude
        opts = setvaropts(opts, opts.VariableNames{4}, 'Type', 'double'); % longitude
        opts = setvaropts(opts, opts.VariableNames{5}, 'Type', 'double'); % temperature
        opts = setvaropts(opts, opts.VariableNames{6}, 'Type', 'double'); % salinity
        opts = setvaropts(opts, opts.VariableNames{7}, 'Type', 'double'); % density
        opts.datalines = [3 Inf]; % skip first 2 lines of metadata
        % opts.SelectedVariableNames = opts.VariableNames(2:end);
        opts.VariableNames = {'time_UTC_','depth','latn','lone','tempc','psal','do','sigma'};
        T_RAW = readtable(fname, opts);
        T_RAW.time_UTC_ = datetime(T_RAW.time_UTC_, 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss''Z''', 'TimeZone', 'UTC');    
        T_RAW.mon_day_yr = T_RAW.time_UTC_;
        T_RAW.mon_day_yr.Format = 'MM/dd/uuuu';
        T_RAW.hh_mm = T_RAW.time_UTC_;
        T_RAW.hh_mm.Format = 'HH:mm';
    
        T = table();
        T.WMO_ID = repmat({fname(4:6)}, height(T_RAW), 1);
        T.Prof_num = NaN(height(T_RAW), 1);
        T.mon_day_yr = T_RAW.mon_day_yr;
        T.hh_mm      = T_RAW.hh_mm;
        T.Lon_E = T_RAW.lone;
        T.Lat_N = T_RAW.latn;
        T.Position_QF = NaN(height(T_RAW), 1);
        T.Pressure_dbar = gsw_p_from_z(-T_RAW.depth, T_RAW.latn);
        T.Pressure_dbar_QF = NaN(height(T_RAW), 1);
        T.Temperature_C = T_RAW.tempc;
        T.Temperature_C_QF = NaN(height(T_RAW), 1);
        T.Salinity_pss = T_RAW.psal;
        T.Salinity_pss_QF = NaN(height(T_RAW), 1);
        T.Sigma_theta_kg_m3 = T_RAW.sigma;
        T.Sigma_theta_kg_m3_QF = NaN(height(T_RAW), 1);
        T.Depth_m = T_RAW.depth;
        T.Depth_m_QC = NaN(height(T_RAW), 1);
        T.Oxygen_umol_kg = T_RAW.do;
        T.Oxygen_umol_kg_QF = NaN(height(T_RAW), 1);
        T.OxygenSat = NaN(height(T_RAW), 1);
        T.OxygenSat_QF = NaN(height(T_RAW), 1);
        T.Nitrate_umol_kg = NaN(height(T_RAW), 1);
        T.Nitrate_umol_kg_QF = NaN(height(T_RAW), 1);
        T.pHinsitu_Total = NaN(height(T_RAW), 1);
        T.pHinsitu_Total_QF = NaN(height(T_RAW), 1);
        T.TALK_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.DIC_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.pH_ALGORITHM_insitu_total = NaN(height(T_RAW), 1);
        T.NO3_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.PO4_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.SILICATE_ALGORITHM_umol_kg = NaN(height(T_RAW), 1);
        T.DIC_pHTA_umol_kg = NaN(height(T_RAW), 1);
        T.PCO2_pHTA_uatm = NaN(height(T_RAW), 1);
    end
    
    clear T_RAW;
    %% Compute ESPER
    sdn = datenum(T.mon_day_yr + timeofday(T.hh_mm));
    DesireVar = [1, 2, 3, 4, 5, 6]; %TA, DIC, pH, phosphate, nitrate, silicate
    OutCoords = horzcat(T.Lon_E, T.Lat_N, T.Depth_m); %lon, lat, depth
    PredictorTypes = [1 2 6]; % PSAL, TEMP, DOXY_ADJ, 
    MeasEsper = [T.Salinity_pss, T.Temperature_C, T.Oxygen_umol_kg];
    refyear = year(sdn)+month(sdn)/12; % reference year for OA adjustment
    Equations = 7; % S, T, O2
    
    % ESPER-LIR for TA, DIC, pH, and Nitrate
    disp('Starting ESPER_LIR...');
    [EspLir,~,~] = ESPER_LIR(DesireVar, OutCoords, MeasEsper,...
        PredictorTypes, 'Equations', Equations, 'EstDates', refyear); 
    disp('Finished ESPER_LIR');
    % save variables to 's' structure
    s.ta_esplir  = EspLir.TA;
    s.dic_esplir = EspLir.DIC;
    s.pH_esplir  = EspLir.pH;
    s.no3_esplir = EspLir.nitrate;
    s.po4_esplir = EspLir.phosphate;
    s.sioh4_esplir = EspLir.silicate;
    % Esper NN
    disp('Starting ESPER_NN...');
    [EspNN,~] = ESPER_NN(DesireVar, OutCoords, MeasEsper,...
        PredictorTypes, 'Equations', Equations, 'EstDates', refyear); 
    disp('Finished ESPER_NN');
    % save variables to 's' structure
    s.ta_espnn  = EspNN.TA;
    s.dic_espnn = EspNN.DIC;
    s.pH_espnn  = EspNN.pH;
    s.no3_espnn = EspNN.nitrate;
    s.po4_espnn = EspNN.phosphate;
    s.sioh4_espnn = EspNN.silicate;
    
    clear DesireVar OutCoords PredictorTypes MeasEsper refyear Equations
    %% CANB
    % canb = CANYONB(sdn, T.Lat_N, T.Lon_E, T.Pressure_dbar, T.Temperature_C, T.Salinity_pss, T.Oxygen_umol_kg, {'NO3', 'AT', 'CT', 'pH', 'SiOH4', 'PO4'}); % make sure same weights
    % 
    % s.pH_canb = canb.pH + (canb.pH*0.0404 - 0.3168); % additional correction to make it consistent to 'spec pH'
    % s.ta_canb  = canb.AT;
    % s.dic_canb  = canb.CT;
    % s.no3_canb  = canb.NO3;
    % s.po4_canb = canb.PO4;
    % s.sioh4_canb = canb.SiOH4;
    %% CANB fast
    % --- find rows with no NaNs across all CANYONB inputs ---
    measCanb   = [T.Pressure_dbar, T.Temperature_C, T.Salinity_pss, T.Oxygen_umol_kg];
    validIdx   = find(all(~isnan(measCanb), 2));
    nTotal     = height(T);
    nValid     = length(validIdx);
    fprintf('Running CANYONB on %d of %d rows\n', nValid, nTotal);
    
    % --- run CANYONB on valid rows only ---
    canb = CANYONB(...
        sdn(validIdx), ...
        T.Lat_N(validIdx), ...
        T.Lon_E(validIdx), ...
        T.Pressure_dbar(validIdx), ...
        T.Temperature_C(validIdx), ...
        T.Salinity_pss(validIdx), ...
        T.Oxygen_umol_kg(validIdx), ...
        {'NO3', 'AT', 'CT', 'pH', 'SiOH4', 'PO4'});
    %%
    % --- reindex outputs to full table size ---
    pH_canb_raw  = expandToFull(canb.pH,    validIdx, nTotal);
    s.pH_canb    = pH_canb_raw + (pH_canb_raw * 0.0404 - 0.3168);
    s.ta_canb    = expandToFull(canb.AT,    validIdx, nTotal);
    s.dic_canb   = expandToFull(canb.CT,    validIdx, nTotal);
    s.no3_canb   = expandToFull(canb.NO3,   validIdx, nTotal);
    s.po4_canb   = expandToFull(canb.PO4,   validIdx, nTotal);
    s.sioh4_canb = expandToFull(canb.SiOH4, validIdx, nTotal);
    %% Add to table: compute the average between ESPERLIR,ESPERNN,CANB
    T.TALK_ALGORITHM_umol_kg   = mean(cat(3,s.ta_esplir,  s.ta_espnn,  s.ta_canb), 3,'omitnan');
    T.DIC_ALGORITHM_umol_kg  = mean(cat(3,s.dic_esplir, s.dic_espnn, s.dic_canb),3,'omitnan');
    T.pH_ALGORITHM_insitu_total  = mean(cat(3,s.pH_esplir,  s.pH_espnn,  s.pH_canb), 3,'omitnan');
    T.NO3_ALGORITHM_umol_kg   = mean(cat(3,s.no3_esplir, s.no3_espnn, s.no3_canb),3,'omitnan');
    T.PO4_ALGORITHM_umol_kg = mean(cat(3,s.po4_esplir, s.po4_espnn, s.po4_canb),3,'omitnan');
    T.SILICATE_ALGORITHM_umol_kg = mean(cat(3,s.sioh4_esplir, s.sioh4_espnn, s.sioh4_canb),3,'omitnan');
    
    %% Now use CO2SYS with TA_ALGORITHM, pH, nutrients to estimate DIC and PCO2.
    % need everything in vectors
    ta = T.TALK_ALGORITHM_umol_kg;
    pHin = T.pHinsitu_Total;
    tc = T.Temperature_C;
    psal = T.Salinity_pss;
    pres = T.Pressure_dbar;
    sil = T.SILICATE_ALGORITHM_umol_kg;
    po4 = T.PO4_ALGORITHM_umol_kg; 
    
    trex = CO2SYSv3(ta, pHin, 1, 3, psal, tc, 25, pres, 0, sil, po4, 0, 0, 1, 10, 1, 2, 2);
    
    % extract variables
    T.DIC_pHTA_umol_kg = trex(:,2);
    T.PCO2_pHTA_uatm = trex(:,4);
    
    clear ta pHin tc psal pres sil po4 trex
    %% -10000000000, -999 -> nan
    T = standardizeMissing(T, [-10000000000, -999]);%% save files
    %% SAVE
    fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS','IOOS glider data 2020',...
    'IOOS_DERIVED_PARAMS',[fname(1:end-4), '_ROMS.csv']);
    % Define metadata lines
    metadata = {
        '// IOOS GLIDER WITH DERIVED PARAMETERS'
        '// Date: 2026-03-27'
        '// Estimated parameters are average of CANYONB, ESPER-LIR, and ESPER-NN'
        '// Uncertainty: TALK_ALGORITHM_umol_kg = 15 umol/kg, DIC_ALGORITHM_umol_kg = 20 umol/kg, pH_ALGORITHM_insitu_total = 0.05, NO3_ALGORITHM_umol_kg = 3 umol/kg, pHinsitu_Total = 0.01, Nitrate_umol_kg = 1 umol/kg, DIC_pHTA_umol_kg = 10 umol/kg (DIC calculated from measured pH + algorithm TA), Oxygen_umol_kg = 1%'
    };
    
    % Open file and write metadata
    fid = fopen(fsave, 'w');
    for i = 1:length(metadata)
        fprintf(fid, '%s\n', metadata{i});
    end
    fclose(fid);
    writetable(T,fsave,"WriteMode","append",WriteVariableNames=true);
    clear fsave;
end
%%
% --- helper (if not already defined) ---
function out = expandToFull(vals, idx, n)
    out = nan(n, 1);
    out(idx) = vals;
end