%% Run mode
% UPDATE_MODE = true:  only fetch/process profiles newer than the latest
%   entry already saved in for_distribution/CORE_ARGO_ROMS.csv, then
%   append the new rows to that file. Skips reprocessing years of old
%   data (minutes instead of hours).
% UPDATE_MODE = false: full historical reprocess from 2018-01-01,
%   overwriting for_distribution/CORE_ARGO_ROMS.csv (the original,
%   slow behavior of this script).
UPDATE_MODE = true;
distfile = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
    'for_distribution','CORE_ARGO_ROMS.csv');

if UPDATE_MODE
    % Pull just the date/time columns from the existing distribution file
    % (avoid loading all ~5M rows x 12 vars) to find the last profile
    % already saved, and use that as the cutoff for this run.
    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','dd-MMM-yyyy HH:mm');
    cutoff_dt = max(dt_existing);
    fprintf('UPDATE_MODE: latest profile already saved is %s. Fetching newer data only.\n',...
        datestr(cutoff_dt));
    datestart = str2double(datestr(cutoff_dt,'yyyymmddHHMMSS')); % index-file date format
    datestart_dnum = datenum(cutoff_dt); % datenum format, used later for T.SDN
else
    datestart = 20180101000000; % format matching index file
    datestart_dnum = datenum(2018,1,1);
end
%% Import Argo profile index
Sprof_tbl = import_prof_index();
%% Filter index file for space and time bounds
basepath = 'https://data-argo.ifremer.fr/dac/';
% Set lat/lon bounds for study area
latbounds = [30 48];
lonbounds = [-134 -115];
% bounding box
ibounds = Sprof_tbl.latitude > latbounds(1) &...
    Sprof_tbl.latitude < latbounds(2) &...
    Sprof_tbl.longitude > lonbounds(1) &...
    Sprof_tbl.longitude < lonbounds(2) &...
    Sprof_tbl.date > datestart;
% All CORE Argo profs in the bounding time/space window
filtered_profiles = Sprof_tbl(ibounds, :); 
%% Get all float WMOs from index file
parts = split(filtered_profiles.file, '/');
WMO = str2double(parts(:,2));
filtered_profiles.WMO = WMO;
floats = unique(WMO);
%% Map of all matching profiles
n_floats = length(unique(WMO));
n_profiles = height(filtered_profiles);
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,filtered_profiles.longitude,...
    'magenta','filled')
hold on;
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_floats, n_profiles);
legend(lg1)
if UPDATE_MODE
    title(sprintf('Float profiles since %s',datestr(datestart_dnum,'dd-mmm-yyyy')));
    fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
        'CORE_ARGO_OUTPUT','MAP_CORE_ARGO_UPDATE.jpg');
else
    title('Float Profiles since 2018');
    fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
        'CORE_ARGO_OUTPUT','MAP_CORE_ARGO_2018_2025.jpg');
end
saveas(fig,fsave)
clear fsave % reset the variable name
%% Download NC files from the GDAC server and convert data to csv format
% DOI for 12/08/2025 snapshot: 10.17882/42182
% Define variable names
linvars =...
    {'JULD','LATITUDE','LONGITUDE','CYCLE_NUMBER','JULD_QC','POSITION_QC'};
matvars =...
    {'PRES_ADJUSTED','TEMP_ADJUSTED','PSAL_ADJUSTED',...
    'PRES_ADJUSTED_QC','TEMP_ADJUSTED_QC','PSAL_ADJUSTED_QC'};

% Create output variable names with SDN instead of JULD
linvars_output = {'SDN','LATITUDE','LONGITUDE',...
    'CYCLE_NUMBER','SDN_QC','POSITION_QC'};
allVars_output = horzcat({'WMO_ID'},linvars_output, matvars);

% Initialize empty table with proper variable names and types
T = table('Size', [0, length(allVars_output)], ...
'VariableTypes', repmat({'double'}, 1, length(allVars_output)), ...
'VariableNames', allVars_output);

% Reference date for JULD (1950-01-01 for Argo floats)
reference_date = datenum(1950, 1, 1, 0, 0, 0);

% Keep track of floats without data
floatswithoutdata = cell(length(floats), 1);

% Loop through floats
for i = 1:length(floats)
    float = num2str(floats(i));
    fname = [float,'_prof.nc'];
    fsave = fullfile('C:','Users','bwerb','Documents',...
        'CUGNROMS','CORE_NC_QC',fname);
    furl = fullfile(basepath,parts{i,1},num2str(floats(i)),fname);
    
    % Check if file already exists, if not download it. In UPDATE_MODE we
    % always re-download: each float's *_prof.nc is a cumulative file, so
    % a locally cached copy from a prior run can't contain profiles added
    % to the float since then.
    if ~isfile(fsave) || UPDATE_MODE
        try
            websave(fsave,furl);
            fprintf('Downloaded: %s\n', fname);
        catch
            floatswithoutdata{i,:} = fsave;
            continue
        end
    else
        fprintf('File already exists, skipping download: %s\n', fname);
    end
    
    % Read data into struct
    data = struct;
    for j = 1:length(matvars)
        data.(matvars{j}) = ncread(fsave,matvars{j});
        
        % Convert QC variables from char to double if needed
        % QC variables are stored as char arrays in NetCDF
        if contains(matvars{j}, '_QC')
            if ischar(data.(matvars{j})) || iscell(data.(matvars{j}))
                % For 2D char arrays, convert each element
                data.(matvars{j}) = double(data.(matvars{j}))...
                    - double('0');
                % Convert -16 to NaN
                data.(matvars{j})(data.(matvars{j}) == -16) = NaN;
                % Also convert anything outside 0-9 to NaN
                data.(matvars{j})(data.(matvars{j}) < 0 |...
                    data.(matvars{j}) > 9) = NaN;
            end
        end
    end
    
    % Add lin (array) variables to struct
    % Get N_LEVELS first (lenth dimension of matrix)
    info = ncinfo(fsave);
    N_LEVELS = info.Dimensions(11).Length;
    for k = 1:length(linvars)
        data.(linvars{k}) = ncread(fsave,linvars{k});
        
        % Convert JULD to SDN (MATLAB Serial Date Number)
        if strcmp(linvars{k}, 'JULD')
            data.SDN = reference_date + data.JULD;
        elseif strcmp(linvars{k}, 'JULD_QC')
            % Convert JULD_QC from char to double and assign to SDN_QC
            if ischar(data.JULD_QC)
                data.SDN_QC = double(data.JULD_QC) - double('0');
                % Convert -16 to NaN
                data.SDN_QC(data.SDN_QC == -16) = NaN;
                data.SDN_QC(data.SDN_QC < 0 | data.SDN_QC > 9) = NaN;
            else
                data.SDN_QC = data.JULD_QC;
            end
        end
        
        % Convert other QC variables from char to double if needed
        if contains(linvars{k}, '_QC') && ~strcmp(linvars{k}, 'JULD_QC')
            if ischar(data.(linvars{k}))
                data.(linvars{k}) = double(data.(linvars{k}))...
                    - double('0');
                % Convert -16 to NaN
                data.(linvars{k})(data.(linvars{k}) == -16) = NaN;
                data.(linvars{k})(data.(linvars{k}) < 0 |...
                    data.(linvars{k}) > 9) = NaN;
            end
        end
        
        % Replicate linear variables to match N_LEVELS
        data.(linvars{k}) = repmat(data.(linvars{k})', N_LEVELS, 1);
    end
    
    % Replicate SDN and SDN_QC to match N_LEVELS
    data.SDN = repmat(data.SDN', N_LEVELS, 1);
    data.SDN_QC = repmat(data.SDN_QC', N_LEVELS, 1);
    
    % Create WMO_ID column with the current float ID for all rows
    data.WMO_ID = repmat(floats(i), numel(data.SDN), 1);

    % Create table for this float using SDN instead of JULD
    T1 = table(data.WMO_ID(:),data.SDN(:), data.LATITUDE(:),...
        data.LONGITUDE(:), data.CYCLE_NUMBER(:), data.SDN_QC(:),...
        data.POSITION_QC(:), data.PRES_ADJUSTED(:),...
        data.TEMP_ADJUSTED(:), data.PSAL_ADJUSTED(:),...
        data.PRES_ADJUSTED_QC(:), data.TEMP_ADJUSTED_QC(:),...
        data.PSAL_ADJUSTED_QC(:), 'VariableNames', allVars_output);
    
    % Append to main table
    T = [T; T1];
end

% Clean up floatswithoutdata (remove empty cells)
floatswithoutdata =...
    floatswithoutdata(~cellfun('isempty', floatswithoutdata));
if UPDATE_MODE
    fsavefloatswithoutdata = fullfile('C:','Users','bwerb','Documents',...
        'CUGNROMS','CORE_ARGO_OUTPUT','floats_without_data_UPDATE.mat');
else
    fsavefloatswithoutdata = fullfile('C:','Users','bwerb','Documents',...
        'CUGNROMS','CORE_ARGO_OUTPUT','floats_without_data.mat');
end
save(fsavefloatswithoutdata,"floatswithoutdata")
% Display results
fprintf('Total rows in combined table: %d\n', height(T));
fprintf('Number of floats successfully processed: %d\n',...
    i - length(floatswithoutdata));
fprintf('Number of floats with missing/failed data: %d\n',...
    length(floatswithoutdata));

if ~isempty(floatswithoutdata)
    fprintf('\nFloats without data:\n');
    disp(floatswithoutdata);
end

% Calculate depth from pressure
T.DEPTH_ADJUSTED = sw_dpth(T.PRES_ADJUSTED, T.LATITUDE);

clear fsave;
%% Create QC table: QC flags = 1 'Good data'.
good_Data = 1;
idx_QC = T.SDN_QC == good_Data &...
    T.POSITION_QC == good_Data &...
    T.PRES_ADJUSTED_QC == good_Data &...
    T.TEMP_ADJUSTED_QC == good_Data &...
    T.PSAL_ADJUSTED_QC == good_Data;

% Set lat/lon bounds for study area. _prof downloads all profiles. Just
% want ones inside area (and, in UPDATE_MODE, only ones newer than
% datestart_dnum -- these floats' _prof.nc files hold their full history,
% not just the new profiles).
latbounds = [30 48];
lonbounds = [-134 -115];
% bounding box
ibounds = T.LATITUDE > latbounds(1) &...
    T.LATITUDE < latbounds(2) &...
    T.LONGITUDE > lonbounds(1) &...
    T.LONGITUDE < lonbounds(2) &...
    T.SDN > datestart_dnum;
T_QC = T(idx_QC & ibounds,:);

if isempty(T_QC)
    fprintf('No new QC''d profiles found since %s. Nothing to append.\n',...
        datestr(datestart_dnum));
    return
end
%% Save to CSV
if UPDATE_MODE
    fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
        'CORE_ARGO_OUTPUT','CORE_ARGO_UPDATE_QC.csv');
else
    fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
        'CORE_ARGO_OUTPUT','CORE_ARGO_2018_2025_QC.csv');
end
writetable(T_QC,fsave);
clear fsave;
%% T,S -> TA
DesiredVariables = 1;  % 1 = Total Alkalinity (TA)
% Define measurement locations [longitude, latitude, depth]
OutputCoordinates = [T_QC.LONGITUDE, T_QC.LATITUDE, T_QC.DEPTH_ADJUSTED];

% Define predictor measurements [Salinity, Temperature]
PredictorMeasurements = [T_QC.PSAL_ADJUSTED, T_QC.TEMP_ADJUSTED];

% Specify which parameters are in which columns
PredictorTypes = [1, 2];  % 1=Salinity, 2=Temperature

% Specify to use equation 8 (S, T only)
Equations = 8;

% Call the function
[Estimates, Uncertainties] = ...
    ESPER_Mixed(DesiredVariables, OutputCoordinates, ...
              PredictorMeasurements, PredictorTypes, ...
              'Equations', Equations);
T_QC.TALK_ESPER_MIXED__mol_kg_T_S_LAT_LON = Estimates.TA;
T_QC.TALK_ESPER_MIXED__mol_kg_T_S_LAT_LON_Uncertainty = Uncertainties.TA;

%% Normalize table
% Final table, no QC vars since all = 1 'good data'
T_CORE_ARGO = table();
T_CORE_ARGO.WMO_ID = T_QC.WMO_ID;
T_CORE_ARGO.Prof_num = T_QC.CYCLE_NUMBER;
T_CORE_ARGO.mon_day_yr = datetime(T_QC.SDN,...
    'ConvertFrom', 'datenum','Format','dd-MMM-yyyy'); % matches existing CSV
T_CORE_ARGO.hh_mm = datetime(T_QC.SDN,...
    'ConvertFrom', 'datenum','Format','HH:mm');
T_CORE_ARGO.Lon_E = T_QC.LONGITUDE;
T_CORE_ARGO.Lat_N = T_QC.LATITUDE;
T_CORE_ARGO.Pressure_dbar = T_QC.PRES_ADJUSTED;
T_CORE_ARGO.Temperature_C = T_QC.TEMP_ADJUSTED;
T_CORE_ARGO.Salinity_pss = T_QC.PSAL_ADJUSTED;
T_CORE_ARGO.Depth_m = T_QC.DEPTH_ADJUSTED;
T_CORE_ARGO.TALK_ESPER_umol_kg = T_QC.TALK_ESPER_MIXED__mol_kg_T_S_LAT_LON;
% Note: the existing distribution CSV does not carry a per-row
% uncertainty column (see metadata header: fixed +/-15 umol/kg), so it's
% left out here to keep the column layout matching for append.

% pot den anom (sigma-theta)
potT   = theta(T_CORE_ARGO.Pressure_dbar, T_CORE_ARGO.Temperature_C,...
    T_CORE_ARGO.Salinity_pss, 0);
T_CORE_ARGO.Sigma_theta_kg_m3    = density(T_CORE_ARGO.Salinity_pss, potT) - 1000; 
%% Map location of all profiles
% Draw a box for reference
n_floats = length(unique(T_CORE_ARGO.WMO_ID));
combinedID =...
    string(T_CORE_ARGO.WMO_ID) + "_" + string(T_CORE_ARGO.Prof_num);
[C,IA,IC] = unique(combinedID);
n_profiles = length(C);
lat_box =...
    [latbounds(1), latbounds(1), latbounds(2), latbounds(2), latbounds(1)];
lon_box =...
    [lonbounds(1), lonbounds(2), lonbounds(2), lonbounds(1), lonbounds(1)];
T_plot = T_CORE_ARGO(IA,:);
% Create figure
fig = figure;
geobasemap('satellite')
geoscatter(T_plot.Lat_N,T_plot.Lon_E,'magenta','filled')
hold on;
geoplot(lat_box, lon_box, 'b-', 'LineWidth', 2);
geolimits(...
    [latbounds(1)-2, latbounds(2)+2], [lonbounds(1)-2, lonbounds(2)+2]);
lg1 = sprintf('CORE floats, profiles: %d, %d',n_floats, n_profiles);
legend(lg1)
if UPDATE_MODE
    title(sprintf('Float profiles since %s',datestr(datestart_dnum,'dd-mmm-yyyy')));
    fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
        'CORE_ARGO_OUTPUT','CORE_ARGO_MAP_UPDATE.jpg');
else
    title('Float Profiles since 2018');
    fsave = fullfile('C:','Users','bwerb','Documents','CUGNROMS',...
        'CORE_ARGO_OUTPUT','CORE_ARGO_MAP_2018_2025.jpg');
end
saveas(fig,fsave)
clear fsave;
%% Save
% distfile already points at for_distribution/CORE_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_CORE_ARGO,distfile,"WriteMode","append",...
        "WriteVariableNames",false);
    fprintf('Appended %d rows (%d profiles, %d floats) to %s\n',...
        height(T_CORE_ARGO), n_profiles, n_floats, distfile);
else
    % Define metadata lines
    metadata = {
        '// CORE_ARGO_ROMS_v6'
        ['// Date: ', datestr(now,'yyyy-mm-dd')]
        '// Estimated parameters are average ESPER-LIR, and ESPER-NN'
        '// Uncertainty: TALK_ALGORITHM_umol_kg = 15 umol/kg'
    };

    % 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_CORE_ARGO,distfile,"WriteMode","append","WriteVariableNames",true);
end
%%
close all;