% Derive_empirical_parameters.m
%
% Yui Takeshita
% MBARI
% 4/11/2023


clear
close all

% .mat files produced by Extract_CUGN_to_mat.m
ddir = 'C:\Users\bwerb\Documents\CUGNROMS\CUGN\data';

% list of filenames for CUGN 2019 data. CUGN_line_93 is deliberately not
% included here -- it has no oxygen data, so ESPER/CANYONB (both of which
% require oxygen) can't run on it.
fnames = {'CUGN_along.mat';'CUGN_line_56.mat';'CUGN_line_66.mat';'CUGN_line_80.mat';'CUGN_line_90.mat'};

% loop through each one and derive various parameters

for i = 1:length(fnames)

    % load 2019 data
    disp(['Starting ',fnames{i},'...']);
    load(fullfile(ddir,fnames{i}),'s');
    
    % create inputs for algorithms
    latin = ones(size(s.tc,1),1) * s.lat;
    lonin = ones(size(s.tc,1),1) * s.lon;
    sdnin = ones(size(s.tc,1),1) * s.sdn;
    
    % *** commented 1/7/2026; LIRs are replaced by ESPERLIR****
    % % ========== LIRS =====================
    % 
    % % LIARv2
    % disp('Starting LIARv2...');
    % ta_liar = LIAR_TSO(latin(:), lonin(:),  s.depth(:), s.psal(:), s.tc(:), s.doxy(:));
    % s.ta_liar = reshape(ta_liar, size(s.tc));
    % 
    % % LIPHR no oA
    % disp('Starting LIPHR no OA...');
    % pH_liphr = LIPHR_TSO(latin(:), lonin(:),  s.depth(:), s.psal(:), s.tc(:), s.doxy(:), false, sdnin(:));
    % s.pH_liphr = reshape(pH_liphr, size(s.tc));
    % 
    % % LINR
    % disp('Starting LINR...');
    % no3_linr = LINR_TSO(latin(:), lonin(:),  s.depth(:), s.psal(:), s.tc(:), s.doxy(:));
    % s.no3_linr = reshape(no3_linr, size(s.tc));
    % ************************************************* 

    % ============ CANYONB ================
    
    % TA, Nitrate
    % Estimate TA, DIC, pH, NO3 from CANYONB
    disp('Starting CANYONB...');
    trex = CANYONB(sdnin(:), latin(:), lonin(:), s.depth(:), s.tc(:), s.psal(:), s.doxy(:), {'NO3', 'AT', 'CT', 'pH'});
    % adjust pH so that it matches spec pH. carter_corr isn't on this
    % machine's path -- this is its documented equivalent (same formula
    % used inline for the same purpose in BGC_ARGO.m).
    pHcanb = trex.pH + (trex.pH.*0.0404 - 0.3168);
    disp('CANYONB finished');

    % code to save these into s structure
    s.ta_canb = reshape(trex.AT, size(s.tc));
    s.dic_canb = reshape(trex.CT, size(s.tc));
    s.pH_canb = reshape(pHcanb, size(s.tc));
    s.no3_canb = reshape(trex.NO3, size(s.tc));

    % ============ ESPER ====================
    
    DesireVar = [1, 2, 3, 5]; %TA, DIC, pH, Nitrate
    OutCoords = horzcat(lonin(:), latin(:), s.depth(:)); %lon, lat, depth
    PredictorTypes = [1 2 6]; % PSAL, TEMP, DOXY_ADJ, 
    MeasEsper = [s.psal(:), s.tc(:), s.doxy(:)];
    refyear = year(sdnin(:)); % 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  = reshape(EspLir.TA, size(s.tc));
    s.dic_esplir = reshape(EspLir.DIC, size(s.tc));
    s.pH_esplir  = reshape(EspLir.pH, size(s.tc));
    s.no3_esplir = reshape(EspLir.nitrate, size(s.tc));
    
    % 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  = reshape(EspNN.TA, size(s.tc));
    s.dic_espnn = reshape(EspNN.DIC, size(s.tc));
    s.pH_espnn  = reshape(EspNN.pH, size(s.tc));
    s.no3_espnn = reshape(EspNN.nitrate, size(s.tc));
    
    % ESPER Mix, which is just an average 
    s.ta_espmix = reshape(mean(horzcat(EspLir.TA, EspNN.TA),2), size(s.tc));
    s.dic_espmix = reshape(mean(horzcat(EspLir.DIC, EspNN.DIC),2), size(s.tc));
    s.pH_espmix = reshape(mean(horzcat(EspLir.pH, EspNN.pH),2), size(s.tc));
    s.no3_espmix = reshape(mean(horzcat(EspLir.nitrate, EspNN.nitrate),2), size(s.tc));
    
    % average the non-LIR parameters
    s.ta_avg   = mean(cat(3,s.ta_esplir,  s.ta_espnn,  s.ta_canb), 3,'omitnan');
    s.dic_avg  = mean(cat(3,s.dic_esplir, s.dic_espnn, s.dic_canb),3,'omitnan');
    s.no3_avg  = mean(cat(3,s.no3_esplir, s.no3_espnn, s.no3_canb),3,'omitnan');
    s.pH_avg   = mean(cat(3,s.pH_esplir,  s.pH_espnn,  s.pH_canb), 3,'omitnan');
    
    % write QF flags
    
    s.position_QC = ones(size(s.lat));
    % pres, temp, sal, density, depth, doxy, derived parameters
    vars = {'pres', 'tc', 'psal', 'pdens', 'depth', 'doxy', 'pH_avg', 'ta_avg', 'dic_avg', 'no3_avg'};
    % basically good, or missing data
    for v=1:length(vars)
        s.([vars{v},'_QC']) = ones(size(s.(vars{v})));
        s.([vars{v},'_QC'])(isnan(s.([vars{v}]))) = 9;
    end
    
    % save data
    save(fullfile(ddir,[fnames{i}(1:end-4),'_withEmp.mat']), 's');
    clearvars -except ddir fnames

end































