%% Run mode
% UPDATE_MODE = true:  only recompute the expensive ESPER-LIR/ESPER-NN/
%   CANYONB/CO2SYS estimates for profiles newer than the latest entry
%   already saved in for_distribution/BGC_ARGO_ROMS.csv, then append
%   just those new rows (minutes instead of hours).
% UPDATE_MODE = false: full historical reprocess, overwriting
%   for_distribution/BGC_ARGO_ROMS.csv (the original, slow behavior).
%
% NOTE: unlike CORE_ARGO.m, this script never downloads anything itself
% -- it reads whatever float files are already sitting in
% soccom_snapshot_dir (below) and MBARI_floatviz_QC/. Which floats/files
% get read is deliberately left unfiltered by the cutoff, since that's
% driven by the (separately stale) local argo_synthetic-profile_index.txt
% and must keep covering every float that has a local ODV file, not just
% recently-active ones. UPDATE_MODE only skips ESPER/CANYONB/CO2SYS for
% rows already in the CSV -- it will not find a profile that isn't
% already sitting in one of those two folders. Refresh those folders
% (new SOCCOM/GO-BGC ODV snapshot, updated FloatViz QC pull) first if you
% need genuinely new float data.
UPDATE_MODE = true;

% SOCCOM/GO-BGC ODV snapshot directory -- THIS IS A MANUAL DOWNLOAD STEP.
% Unlike everything else in this script, this data is not fetched
% automatically. It's a dated, citable snapshot release (has its own
% DOI), not a live feed, so getting a newer one means going to the UCSD
% Library archive, downloading the latest dated release, and pointing
% this variable at the new folder name:
%   https://library.ucsd.edu/dc/object/bb1310816p
soccom_snapshot_dir = "SOCCOM_GO-BGC_LoResQC_LIAR_01Jul2026_odvtxt";

distfile = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
    'for_distribution','BGC_ARGO_ROMS.csv');

if UPDATE_MODE
    optsDate = detectImportOptions(distfile,'NumHeaderLines',5);
    optsDate.SelectedVariableNames = {'mon_day_yr','hh_mm'};
    optsDate = setvartype(optsDate,{'mon_day_yr','hh_mm'},'string');
    T_existing_dates = readtable(distfile,optsDate);
    dt_existing = datetime(...
        T_existing_dates.mon_day_yr + " " + T_existing_dates.hh_mm,...
        'InputFormat','MM/dd/uuuu HH:mm');
    cutoff_dt = max(dt_existing);
    fprintf('UPDATE_MODE: latest profile already saved is %s. Only newer profiles will be (re)computed.\n',...
        datestr(cutoff_dt));
end
%% Read in the synthetic profile index file that includes the basic meta data for all BGC-Argo profiles
bgc_profiles = readtable("argo_synthetic-profile_index.txt");
basepath = 'https://usgodae.org/ftp/outgoing/argo/dac/';

% Set lat/lon bounds for study area
latbounds = [30 48];
lonbounds = [-134 -115];
datestart = 20180101000000; % format matching index file
% bounding box
ibounds = bgc_profiles.latitude > latbounds(1) &...
    bgc_profiles.latitude < latbounds(2) &...
    bgc_profiles.longitude > lonbounds(1) &...
    bgc_profiles.longitude < lonbounds(2) &...
    bgc_profiles.date > datestart;

filtered_profiles = bgc_profiles(ibounds, :); % All BGC-Argo profs in the bounding time/space window
%% Get all float WMOs from index file
parts = split(filtered_profiles.file, '/');
WMO = str2double(parts(:,2));
filtered_profiles.WMO = WMO;
%% Create map with MBARI managed and non-MBARI managed floats
% Get MBARI float WMOs
MBARI_float_list = readtable(fullfile(soccom_snapshot_dir,"MBARI_float_list.txt"));
i = ismember(MBARI_float_list.WMO,WMO);
MBARI_float_list_filtered = MBARI_float_list(i,:);
imbari = ismember(WMO,MBARI_float_list.WMO); % Floats from the index file that match the MBARI float list
n_mbari_floats = length(MBARI_float_list_filtered.WMO);
n_non_mbari_floats = length(unique(filtered_profiles.WMO(~imbari)));
n_mbari = length(filtered_profiles.longitude(imbari));
n_non_mbari = length(filtered_profiles.longitude(~imbari));

% Draw a box for reference
lat_box = [latbounds(1), latbounds(1), latbounds(2), latbounds(2), latbounds(1)];
lon_box = [lonbounds(1), lonbounds(2), lonbounds(2), lonbounds(1), lonbounds(1)];
% Create figure
fig = figure;
geobasemap('satellite')
geoscatter(filtered_profiles.latitude(imbari),filtered_profiles.longitude(imbari),'magenta')
hold on;
geoscatter(filtered_profiles.latitude(~imbari),filtered_profiles.longitude(~imbari),'cyan')
geoplot(lat_box, lon_box, 'b-', 'LineWidth', 2);
geolimits([latbounds(1)-2, latbounds(2)+2], [lonbounds(1)-2, lonbounds(2)+2]);
lg1 = sprintf('MBARI floats, profiles: %d, %d',n_mbari_floats, n_mbari);
lg2 = sprintf('Non-MBARI floats, profiles: %d, %d',n_non_mbari_floats, n_non_mbari);
legend(lg1,lg2)
title('Float Profiles since 2018');
fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
    'BGC_ARGO_OUTPUT','MAP_BGC_ARGO_2018_2025.jpg');
saveas(fig,fsave)
clear fsave;
%% Create map with MBARI managed and non-MBARI managed floats by year
MBARI_float_list = readtable(fullfile(soccom_snapshot_dir,"MBARI_float_list.txt"));
parts = split(filtered_profiles.file, '/');
WMO = str2double(parts(:,2));
imbari = ismember(WMO,MBARI_float_list.WMO); % Floats from the index file that match the MBARI float list
lat_box = [latbounds(1), latbounds(1), latbounds(2), latbounds(2), latbounds(1)];
lon_box = [lonbounds(1), lonbounds(2), lonbounds(2), lonbounds(1), lonbounds(1)];
yrs = 2018:1:2025;
for i = 1:length(yrs)
    datestart = str2double(sprintf('%d%010d', yrs(i), 101000000));  % Jan 1 00:00:00
    dateend = str2double(sprintf('%d%010d', yrs(i), 1231235959));   % Dec 31 23:59:59
    
    ibounds = bgc_profiles.latitude > latbounds(1) &...
        bgc_profiles.latitude < latbounds(2) &...
        bgc_profiles.longitude > lonbounds(1) &...
        bgc_profiles.longitude < lonbounds(2) &...
        bgc_profiles.date >= datestart &...
        bgc_profiles.date <= dateend;
    
    filtered_profiles_yrs = bgc_profiles(ibounds, :);
    
    % Recalculate imbari for the year-filtered data
    parts_yrs = split(filtered_profiles_yrs.file, '/');
    WMO_yrs = str2double(parts_yrs(:,2));
    imbari_yrs = ismember(WMO_yrs, MBARI_float_list.WMO);
    
    fig = figure;
    geobasemap('satellite')
    geoscatter(filtered_profiles_yrs.latitude(imbari_yrs), filtered_profiles_yrs.longitude(imbari_yrs), 'magenta')
    hold on;
    geoscatter(filtered_profiles_yrs.latitude(~imbari_yrs), filtered_profiles_yrs.longitude(~imbari_yrs), 'cyan')
    geoplot(lat_box, lon_box, 'b-', 'LineWidth', 2);
    geolimits([latbounds(1)-2, latbounds(2)+2], [lonbounds(1)-2, lonbounds(2)+2]);
    
    n_mbari = sum(imbari_yrs);
    n_non_mbari = sum(~imbari_yrs);
    
    legend(['MBARI profiles: ',num2str(n_mbari)], ['Non-MBARI profiles: ',num2str(n_non_mbari)])
    title(['Float Profiles in ', num2str(yrs(i))]);
    fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
    'BGC_ARGO_OUTPUT',['map_', num2str(yrs(i)), '.jpg']);
    saveas(fig, fsave)
    close(fig)
    clear fsave;
end
%% Get profile data from ODV text files for MBARI sprof floats
% Create list of filenames from sprof directory: https://library.ucsd.edu/dc/object/bb1310816p
dirPath_snapshot = soccom_snapshot_dir;
files = dir(dirPath_snapshot);
% Remove '.' and '..' entries
files = files(~ismember({files.name}, {'.', '..'}));
% Get just file names
fileNames = {files.name}';
MBARI_WMO = unique(WMO(imbari)); % Unique WMOs
% Filenames -> strings
numericFromCells = cellfun(@(x) str2double(regexp(x, '^\d+', 'match', 'once')), fileNames);
% Find  MBARI WMOs that match any of the WMOs from the filenames
[tf, loc] = ismember(numericFromCells, MBARI_WMO);
indices = find(tf);
% Files with MBARI float data
fnames_mbari_floats = fileNames(indices); % 15 matching files. The remaining 4 are: 3 PMEL floats not processed by mbari. 1 UW original OCR test float w/o pH/NO3 or floatviz files.
MBARI_WMO_hasfile = numericFromCells(indices);
non_MBARI_managed_floats = setdiff(MBARI_float_list_filtered.WMO,MBARI_WMO_hasfile); % see comment above
%% Get profile data from ODV text files for nonMBARI floats (managed by mbari and downloaded via chem\ARGO_PROCESSING\DATA\FLOATVIZ\QC)
% Out of 24 nonMBARI floats, 12 are managed by MBARI with complete ODV
% files
% The remaining 12 are: 8 Argo Canada (1 with OXYGEN_ADJUSTED), 3 TS
% floats, 1 triple O2 float. Further described below
% Create list of filenames from sprof directory: https://library.ucsd.edu/dc/object/bb1310816p
% Absolute path -- websave unreliably resolves relative output paths.
dirPath_floatviz = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
    'MBARI_floatviz_QC');
NONMBARI_WMO = unique(WMO(~imbari)); % Unique WMOs

% Refresh just these WMOs' files from MBARI's live FloatViz server before
% reading the directory. dirPath_floatviz is a shared cache used by other
% things too (700+ files, most outside this study), so we only touch the
% files we're actually about to read -- not the whole folder. Each file
% is a per-float cumulative QC pull (like CORE_ARGO.m's *_prof.nc), so we
% always overwrite rather than skip existing ones, to pick up new
% profiles. Floats with no FloatViz presence (Argo Canada, TS-only, etc)
% 404/403 and are silently skipped, same as before.
% websave creates/truncates its output file before the HTTP request is
% known to succeed, so a failed request (403/404 for floats with no
% FloatViz presence) would otherwise leave a 0-byte file behind -- or
% worse, clobber a good file from a previous run if a live float's
% request happens to fail transiently. Download to a temp file first and
% only move it into place once it's confirmed good.
floatviz_url_base = 'https://www3.mbari.org/lobo/Data/FloatVizData/QC/';
floatviz_download_failed = [];
for w = 1:length(NONMBARI_WMO)
    wmo_str = num2str(NONMBARI_WMO(w));
    fsave = fullfile(dirPath_floatviz, [wmo_str, 'QC.TXT']);
    fsave_tmp = [fsave, '.part'];
    furl = [floatviz_url_base, wmo_str, 'QC.txt'];
    try
        websave(fsave_tmp, furl);
        movefile(fsave_tmp, fsave);
    catch
        floatviz_download_failed = [floatviz_download_failed; NONMBARI_WMO(w)]; %#ok<AGROW>
        if isfile(fsave_tmp)
            delete(fsave_tmp);
        end
    end
end
if ~isempty(floatviz_download_failed)
    fprintf('No FloatViz QC file on the server for %d WMO(s) (expected for Argo Canada/TS-only/etc): %s\n',...
        length(floatviz_download_failed), mat2str(floatviz_download_failed'));
end

files = dir(dirPath_floatviz);
% Remove '.' and '..' entries
files = files(~ismember({files.name}, {'.', '..'}));
% Get just file names
fileNames = {files.name}';
% Filenames -> strings
numericFromCells = cellfun(@(x) str2double(regexp(x, '^\d+', 'match', 'once')), fileNames);
% Find  MBARI WMOs that match any of the WMOs from the filenames
[tf, loc] = ismember(numericFromCells, NONMBARI_WMO);
indices = find(tf);
% Files with MBARI float data
fnames_nonmbari_floats = fileNames(indices);
Non_MBARI_WMO_hasfile = numericFromCells(indices);
% Found 12/24 floats. Remaining floats are:
% Might need these:
% Argo Canada:
% 4902663
% 4902643
% 4902642
% 4902624
% 4902554
% 4902551
% 4902549
% 4901784
% params = readtable("nonmbari_params.txt");
% i = ismember(params.WMO,[4902663, 4902643, 4902642, 4902624, 4902554, 4902551, 4902549, 4901784]); % Cross referenced with https://biogeochemical-argo.org/new-floats-data-table.php?search=4900874 to find these floats
% writetable(params(i,:),'ARGOCanada_params.txt');
% Do not need these
% No Match, P,T,S Floats:
% 4902541 - TS
% 4902474 - TS
% 4902472 - TS
% Triple O2:
% 5906445
%% Concatenate all MBARI floatviz files into a singular table, adding nan columns when there is a missing column
% Store all tables in a cell array
tables = cell(length([fnames_mbari_floats; fnames_nonmbari_floats]), 1);

for i = 1:length(fnames_mbari_floats)
    fname = fullfile(dirPath_snapshot, fnames_mbari_floats(i));
    tables{i} = readtable(fname);
end
for i = 1:length(fnames_nonmbari_floats)
    fname = fullfile(dirPath_floatviz, fnames_nonmbari_floats(i));
    tables{i+length(fnames_mbari_floats)} = readtable(fname);
end

% Get all unique column names
allVars = {};
for i = 1:length(tables)
    allVars = [allVars, tables{i}.Properties.VariableNames];
end
allVars = unique(allVars, 'stable');

% Add missing columns to each table
for i = 1:length(tables)
    missingVars = setdiff(allVars, tables{i}.Properties.VariableNames);
    for j = 1:length(missingVars)
        % Add column filled with NaN
        tables{i}.(missingVars{j}) = nan(height(tables{i}), 1);
    end
    % Reorder columns to match
    tables{i} = tables{i}(:, allVars);
end

% Concatenate all tables
T_combined = vertcat(tables{:});
T_combined.Lon__E_ = T_combined.Lon__E_ - 360;
% ODV files have lon in degrees East, +360 to correct the bounds set above
ibounds = T_combined.Lat__N_ > latbounds(1) &...
        T_combined.Lat__N_ < latbounds(2) &...
        T_combined.Lon__E_ > lonbounds(1) &...
        T_combined.Lon__E_ < lonbounds(2);

% GPS,P,T,S,O2 are 'good data'
T_combined_QC = T_combined(T_combined.QF == 0 & T_combined.QF_1 == 0 &...
    T_combined.QF_2 == 0 & T_combined.QF_3 == 0 &...
    T_combined.QF_6 == 0 & ibounds,:); 

% Replace questionable and bad pH with NAN
ibadpH = T_combined_QC.QF_12 > 0; 
T_combined_QC.pHinsitu_Total_(ibadpH) = NaN;
T_combined_QC.pH25C_Total_(ibadpH) = NaN;

%% Create map from the real data set not meta data
lat_box = [latbounds(1), latbounds(1), latbounds(2), latbounds(2), latbounds(1)];
lon_box = [lonbounds(1), lonbounds(2), lonbounds(2), lonbounds(1), lonbounds(1)];
fig = figure;
geobasemap('satellite')
geoscatter(T_combined_QC.Lat__N_,T_combined_QC.Lon__E_,'magenta')
hold on;
geoplot(lat_box, lon_box, 'b-', 'LineWidth', 2);
geolimits([latbounds(1)-2, latbounds(2)+2], [lonbounds(1)-2, lonbounds(2)+2]);
n_mbari_floats = length(unique(T_combined_QC.Cruise));
n_mbari = length(unique([T_combined_QC.Cruise T_combined_QC.Station],'rows')); % cruise+station gives unique profile identifier
% n_non_mbari = length(filtered_profiles.latitude(~imbari));
lg = sprintf('MBARI floats, profiles: %d, %d',n_mbari_floats, n_mbari);
legend(lg)
title('MBARI Managed Float Profiles since 2018');
fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
'BGC_ARGO_OUTPUT','MBARI_MANAGED_FLOATS_2018_2025.jpg');
saveas(fig, fsave)
clear fsave;
%% Standardize column names
% Define the transformations (original name -> new name, or 'REMOVE')
transforms = {
    'Cruise', 'WMO_ID';
    'Station', 'Prof_num';
    'Type', 'REMOVE';
    'mon_day_yr', '';
    'hh_mm', '';
    'Lon__E_', 'Lon_E';
    'Lat__N_', 'Lat_N';
    'QF', 'Position_QF';
    'IceEvRec', 'REMOVE';
    'Pressure_dbar_', 'Pressure_dbar';
    'QF_1', 'Pressure_dbar_QF';
    'Temperature__C_', 'Temperature_C';
    'QF_2', 'Temperature_C_QF';
    'Salinity_pss_', 'Salinity_pss';
    'QF_3', 'Salinity_pss_QF';
    'Sigma_theta_kg_m_3_', 'Sigma_theta_kg_m3';
    'QF_4', 'Sigma_theta_kg_m3_QF';
    'Depth_m_', 'Depth_m';
    'QF_5', 'Depth_m_QC';
    'Oxygen__mol_kg_', 'Oxygen_umol_kg';
    'QF_6', 'Oxygen_umol_kg_QF';
    'OxygenSat___', 'OxygenSat';
    'QF_7', 'OxygenSat_QF';
    'Nitrate__mol_kg_', 'Nitrate_umol_kg';
    'QF_8', 'Nitrate_umol_kg_QF';
    'Chl_a_mg_m_3_', 'REMOVE';
    'QF_9', 'REMOVE';
    'b_bp700_1_m_', 'REMOVE';
    'QF_10', 'REMOVE';
    'POC_mmol_m_3_', 'REMOVE';
    'QF_11', 'REMOVE';
    'pHinsitu_Total_', 'pHinsitu_Total';
    'QF_12', 'pHinsitu_Total_QF';
    'pH25C_Total_', 'REMOVE';
    'QF_13', 'REMOVE';
    'TALK_LIAR__mol_kg_', 'REMOVE';
    'QF_14', 'REMOVE';
    'DIC_LIAR__mol_kg_', 'REMOVE';
    'QF_15', 'REMOVE';
    'pCO2_LIAR__atm_', 'REMOVE';
    'QF_16', 'REMOVE';
    'CDOM_ppb_', 'REMOVE';
    'QF_17', 'REMOVE';
    'DOWN_IRRAD380_W_m_2_nm_', 'REMOVE';
    'QF_18', 'REMOVE';
    'DOWN_IRRAD443_W_m_2_nm_', 'REMOVE';
    'QF_19', 'REMOVE';
    'DOWN_IRRAD490_W_m_2_nm_', 'REMOVE';
    'QF_20', 'REMOVE';
    'DOWN_IRRAD555_W_m_2_nm_', 'REMOVE';
    'QF_21', 'REMOVE';
    'DOWN_IRRAD412_W_m_2_nm_', 'REMOVE';
    'DOWNWELL_PAR__molQuanta_m_2_sec_', 'REMOVE';
    'Chl_a435_mg_m_3_', 'REMOVE';
};

% Apply transformations
for i = 1:size(transforms, 1)
    oldName = transforms{i, 1};
    newName = transforms{i, 2};
    
    % Check if the variable exists in the table
    if ismember(oldName, T_combined_QC.Properties.VariableNames)
        if strcmp(newName, 'REMOVE')
            % Remove the variable
            T_combined_QC.(oldName) = [];
        elseif ~isempty(newName)
            % Rename the variable
            T_combined_QC.Properties.VariableNames{strcmp(T_combined_QC.Properties.VariableNames, oldName)} = newName;
        end
        % If newName is empty (''), do nothing (keep original name)
    end
end
%% Convert mon_day_yr and hh:mm
T_combined_QC.mon_day_yr.Format = 'MM/dd/uuuu';
T_combined_QC.hh_mm = datetime(T_combined_QC.hh_mm, 'InputFormat', 'HH:mm', 'Format', 'HH:mm');
%% Keep only profiles newer than what's already saved (UPDATE_MODE)
% This is what actually skips the expensive ESPER/CANYONB/CO2SYS work
% below for old profiles -- everything upstream still reads the full
% local file set.
if UPDATE_MODE
    dt_combined = T_combined_QC.mon_day_yr + timeofday(T_combined_QC.hh_mm);
    T_combined_QC = T_combined_QC(dt_combined > cutoff_dt, :);
    if isempty(T_combined_QC)
        fprintf('No new profiles found since %s. Nothing to append.\n',...
            datestr(cutoff_dt));
        return
    end
end
%% Save file
if UPDATE_MODE
    fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
    'BGC_ARGO_OUTPUT','BGC_ARGO_N0_DERIVE_UPDATE.csv');
else
    fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
    'BGC_ARGO_OUTPUT','BGC_ARGO_N0_DERIVE.csv');
end
writetable(T_combined_QC,fsave);
clear fsave;
%% Add ESPER params to MBARI floats
T = T_combined_QC;
%% Calculate all params with ESPER and CANYONB
% sdn = datenum(T.mon_day_yr + timeofday(T.hh_mm));
% [EspMixed, uncertainty] = ESPER_TSO(T.Lat_N, T.Lon_E, T.Depth_m, T.Salinity_pss, T.Temperature_C, T.Oxygen_umol_kg, sdn);
%% 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;
%% 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;
%% 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');

%% Add to table
% T.TALK_ESPER_umol_kg = EspMixed.TA;
% T.DIC_ESPER_umol_kg = EspMixed.DIC;
% T.pH_ESPER_insitu_total = EspMixed.pH;
% T.PO4_ESPER_umol_kg = EspMixed.phosphate;
% T.NO3_ESPER_umol_kg = EspMixed.nitrate;
% T.SILICATE_ESPER_umol_kg = EspMixed.silicate;

% T.TALK_ESPER_umol_kg_Uncertainty = uncertainty.TA;
% T.DIC_ESPER_umol_kg_Uncertainty = uncertainty.DIC;
% T.pH_ESPER_insitu_total_Uncertainty = uncertainty.pH;
% T.PO4_ESPER_umol_kg_Uncertainty = uncertainty.phosphate;
% T.NO3_ESPER_umol_kg_Uncertainty = uncertainty.nitrate;
% T.SILICATE_ESPER_umol_kg_Uncertainty = uncertainty.silicate;
%% 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);
%% -10000000000, -999 -> nan
T = standardizeMissing(T, [-10000000000, -999]);
%% SAVE
% distfile already points at for_distribution/BGC_ARGO_ROMS.csv (set at
% the top of the script, where it's also used to find the update cutoff)
if UPDATE_MODE
    % Append only -- keep the existing metadata/header lines untouched.
    writetable(T,distfile,"WriteMode","append","WriteVariableNames",false);
    fprintf('Appended %d rows to %s\n', height(T), distfile);
else
    % Define metadata lines
    metadata = {
        '// BGC_ARGO_ROMS_v7'
        ['// Date: ', datestr(now,'yyyy-mm-dd')]
        '// 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(distfile, 'w');
    for i = 1:length(metadata)
        fprintf(fid, '%s\n', metadata{i});
    end
    fclose(fid);
    writetable(T,distfile,"WriteMode","append","WriteVariableNames",true);
end
%%
close all;