clear all
close all

%% Calculate DERIVED VARIABLES for BGC sensors

% Path to CPF data:
path_data='/Volumes/CPF/CPF003/Data/TankTests/2021-8-11T20-6-55/';

cast='002'; %
msg_file = [path_data,'CPF003.',cast,'.msg'];
isus_file = [path_data,'CPF003.',cast,'.isus'];
% NOTE: AirCalData only for .001.msg file?

%%  load data from msg files
d=parse_CPFmsg_v3(msg_file);
data=d;

%%  get cal data
get_cpf_cals_mc
% cal.type='4330'; %otherwise, calc_O2_4ARGO.m gets stuck

%% LOW RES data coolumns

% iP=1;
% iT=2;
% iS=3;
% iNO3 = 8;
% iChl=10;
% iBb=11;

% Identify columns in header of low res data
% (adapted to include units and lr header names in CPF .msg files)
iP   = find(strcmp('pres(dBar)',      d.lr_hdr) == 1); % CTD P
iT   = find(strcmp('temp(degC)',      d.lr_hdr) == 1); % CTD T
iS   = find(strcmp('sal(ppt)',      d.lr_hdr) == 1); % CTD S
% iPh83  = find(strcmp('Phase',   d.lr_hdr) == 1); % SBE83 O2 phase delay
% iT83 = find(strcmp('T83', d.lr_hdr) == 1); % SBE83 Temp (degC)
iTo  = find(strcmp('Topt(degC)',   d.lr_hdr) == 1); % optode temp
iTPh = find(strcmp('TPhase(deg)', d.lr_hdr) == 1); % T Phase 4330
iRPh = find(strcmp('RPhase(deg)', d.lr_hdr) == 1); % R Phase 4330
% iBPh = find(strcmp('BPhase', d.lr_hdr) == 1); % B Phase 3830
ipH  = find(strcmp('pH(V)',  d.lr_hdr) == 1); % pH
iChl = find(strcmp('FSig(counts)',   d.lr_hdr) == 1); % CHL fluor
iBb  = find(strcmp('BbSig(counts)',  d.lr_hdr) == 1); % Backscatter
iNO3 = find(strcmp('NO3(V)',    d.lr_hdr) == 1); % Nitrate

%% ------------------- Calculate CHLA ------------------------------------
LR.CHLA= (data.lr_d(:,iChl) - cal.CHL.ChlDC) .* ...
                cal.CHL.ChlScale;

%         LR.FLUORESCENCE_CHLA    = fill0 + fv.bio;
%         LR.FLUORESCENCE_CHLA_QC = fill0 + fv.QC;
%         LR.CHLA                 = fill0 + fv.bio;
%         LR.CHLA_QC              = fill0 + fv.QC;
%         LR.CHLA_ADJUSTED        = fill0 + fv.bio;
%         LR.CHLA_ADJUSTED_QC     = fill0 + fv.QC;
%         LR.CHLA_ADJUSTED_ERROR  = fill0 + fv.bio;

%% ------------------- Calculate BBP ------------------------------------
            % VOLUME SCATTERING FUNCTION (VSF) (m^-1 sr^-1)
            VSF(:) = (data.lr_d(:, iBb) - cal.BB.BetabDC) ...
                .* cal.BB.BetabScale;

            % NOW CALCULATE PARTICLE BACKSCATTER COEFFICIENT
            %X       = 1.13*2*pi; % (Barnes and Antoine, 2014)
            %X       = 1.17*2*pi; % (email from E. Boss 24 May 2016)
            
            % (Proc. Bio-Argo particle backscattering at the DAC level
            % Version 1.2, July 21th 2016
            X      = 1.097*2*pi; % FLBB,BBP processing doc July 2016
            LAMBDA = 700;
            DELTA  = 0.039;      % depolarization ratio
            THETA  = 142;        % FLBBAP2 , DAC manual says 120
%            BETA_SW_ind = find(t_nan == 0);
%            if ~isempty(BETA_SW_ind)
%                for ct = 1:size(BETA_SW_ind,1)
                    [BETA_SW(:), b90sw, bsw] = ...
                        betasw_ZHH2009(LAMBDA,data.lr_d(:,iT), ...
                        THETA, data.lr_d(:,iS), DELTA);
%                end
%            end
            
            %LR.BETA_BACKSCATTERING700 = VSF; % (DOUBLE CHECK SPEC say counts???)
            LR.BETA_BACKSCATTERING700(:) = data.lr_d(:,iBb); % counts
%             LR.BETA_BACKSCATTERING700_QC(:) = fv.QC; % counts
            LR.BBP700(:)    = (VSF(:) - BETA_SW(:)) * X; %b_bp m^-1
%            LR.BBP700_QC(:) = 2; % 3 do not use w/o adjusting ... 6/10/21 modify qcraw flag from 3 to 2.
            
%         VSF                          = fill0 + fv.bio;
%         BETA_SW                      = fill0 + fv.bio;
%         LR.BETA_BACKSCATTERING700    = fill0 + fv.bio;
%         LR.BETA_BACKSCATTERING700_QC = fill0 + fv.QC;
%         LR.BBP700                    = fill0 + fv.bio;
%         LR.BBP700_QC                 = fill0 + fv.QC;
%         LR.BBP700_ADJUSTED           = fill0 + fv.bio;
%         LR.BBP700_ADJUSTED_QC        = fill0 + fv.QC;
%         LR.BBP700_ADJUSTED_ERROR     = fill0 + fv.bio;

%% ------------------- Calculate O2 --------------------------------------

P=data.lr_d(:,iP);
T=data.lr_d(:,iT);
S=data.lr_d(:,iS);
Phase=data.lr_d(:,5);

O2data = calc_O2_4ARGO(S, T, P, Phase, cal.O); % from ARGOPROCESSING/
% OUTPUT O2data    =  [S, T, P, Phase, CalPhase, [O2], O2sol, pO2] 
%                     units: [O2 ] = [mmol/m^3 =uM/L]
%                             O2 sol= % Sat ; [mmol/m^3 =uM/L],
%                             pO2 = ?

%% ------------------- Calculate NITRATE ---------------------------------
% CHECK FOR WL_OFFSET IN ISUS NITRATE CAL
% IF IT EXISTS AND IS NOT THE DEFAULT (210), COMMENT
if isfield(cal, 'N')
    if isfield(cal.N, 'WL_offset') && cal.N.WL_offset ~= 210
        disp(['NON STANDARD WAVE LENGTH OFFSET IN NITRATE CAL = ',...
            num2str(cal.N.WL_offset)]);
    end
end
spec = parse_NO3msg(isus_file); 

Pcor_flag=0; % Pcor_flag = 1 to pressure correct ESW, otherwise don't
ncal=cal.N;
NO3 = calc_FLOAT_NO3(spec, ncal,Pcor_flag); 

% OUTPUT: NO3 =   a matrix of nitrate concentrations (umol/L) and other info
%        [SDN, DarkCur,Pres,Temp,Sal,NO3,BL_B,BL_M,RMS ERROR,WL~240,ABS~240]
%% ------------------- Calculate LIGHT ----------------------------------

%% VALID BIO ARGO RANGE CHECKS - RAW DATA [MIN MAX]
% values and code above from Process_APEX_float.m
% last modification 07/13/21 JP

RCR.P     = [0 12000];
RCR.S     = [26 38]; % from argo parameter list
RCR.T     = [-2.5 40]; % from argo parameter list
RCR.O     = [-5 550]; % from argo parameter list
RCR.PO2   = [-5 5000]; % from argo parameter list; added to this code 5/25/21, TM: is this range right?  Should be max 500??
RCR.OP    = [10 70]; % optode phase, from argo parameter list
RCR.OT    = [-2.5 40]; % optode temperature, from argo parameter list
RCR.CHL   = [-0.1 150]; % argoBGC QC manual 09July2016
RCR.BB700 = [-0.000025 0.1]; % argoBGC QC manual 09July2016
RCR.BB532 = [-0.000005 0.1]; % argoBGC QC manual 09July2016
RCR.CDOM = [-1000 1000]; % Place Holder
RCR.NO3  = [-15 65];
RCR.PH     = [7.0 8.8];
RCR.PHV  = [-1.2 -0.7]; % Range check on pH volts
RCR.IB   = [-100 100]; % Range check on pH Ib, nano amps
RCR.IK   = [-100 100]; % Range check on pH Ik, nano amps
