function [PARK, HR, INFO]=calc_BGC_anchor(path_data,cast);
% clear all
% close all
addpath(genpath('/Volumes/Chem/ARGO_PROCESSING/MFILES/FLOATS/'));
addpath(genpath('/Users/Mau/Documents/MBARIprojects/CPF/CPFPROCESSING'));
addpath(genpath('/Users/Mau/Documents/MATLABhelp_new/ARGO_PROCESSING-master/')); %

%% PURPOSE: Extract data from .msg, .isus and .dura to 
%           CALCULATE DERIVED VARIABLES for BGC sensors (low res data, LR)
%           Mimicking data output structures in
%           ARGO_PROCESSING/Process_APEX_float.m (and lots of code from
%           that function). Thus, saving HR data in .mat file as well.
%           NOTE: HR structure from this output is analogous to the one produced for Argo
%           floats by Tanya and Josh. HR data is updated in my local folder 
%           with derived CTD variables (e.g. CT, SA, RHO from TEOS and w_ctd), 
%           later when running plot_cpf_HRctd.m

%% INPUT:
% path_data : path to raw data files (.msg, .isus, .dura)
% cast: string with cast number, .i.e, '001'

%% CHANGES:
% -11/17/2021: MC added float_type=CPF to INFO structure (needed for ODV
% coversion)
% -01/06/2021: modify identification of root path to account for a
% different path name (e.g. in a local direcotry);

%%  Path to CPF data, for testing:
% % path_data='/Volumes/CPF/CPF003/Data/TankTests/2021-8-11T20-6-55/'; % use parse_CPFmsg_v3(msg_file);
% % NOTE: AirCalData only for .001.msg file?
% 
% % path_data='/Volumes/CPF/CPF003/Data/TankTests/2021-8-25T19-5-6/'; % max profiles was set to 2, updated firmware
% path_data='/Volumes/CPF/CPF003/Data/TankTests/2021-8-25T21-57-38/'; % latest tank test
% % use parse_CPFmsg_v4(msg_file);

% cast='005'; % select path and cast cast for testing
% path_data='/Volumes/CPF/CPF003/Data/ParagonTests/2021-9-16T16-33-30/';

% cast='001';
% path_data='/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/ParagonTests/MICRORIDER/2021-11-16T19-16-49/';

% cast='006';
% cast='014';
% cast='020';
% path_data='/Volumes/CPF/CPF003/Data/ParagonTests/2022-2-23T16-59-6/'; % Tests wether there is a missalignment for DOXY

% cast='003';
% path_data= '/Volumes/CPF/CPF003/Data/ParagonTests/2022-3-14T16-30-13/'; % profiles 001, 002, 003, 005 and 013 have 6h data sampled at 2'

% dirs.temp      = path to temporary working dir
dirs.temp=path_data;
% msg_list   = ls([dirs.temp,'*.msg']); % get list of file names to process

msg_file = ['CPF003.',cast,'.msg'];
% isus_file = ['CPF003.',cast,'.isus'];

NO3_file = regexprep(msg_file,'msg','isus');
pH_file  = regexprep(msg_file,'msg','dura');

cast_num = regexp(msg_file,'\d+(?=\.msg)','once','match');

lr_ind_chk = 0; % indice toggle - find indices once per float
hr_ind_chk = 0; % toggle

%%  load data from msg files
% d=parse_CPFmsg_v6([dirs.temp,msg_file]); % modifying to extract anchor data
d=parse_CPFmsg_v5([dirs.temp,msg_file]); % should work for Tests done in
anchor=ReadMsgFile_mc(dirs.temp,msg_file); % code from Gene, modified to putput data in right format

data=d; % to match T&J code

% try overwriting lr profile data with anchor data
data.lr_hdr=anchor.park_hdr; %: {16×1 cell}
data.lr_d=anchor.park_d; %: [180×16 double]

d.lr_hdr=anchor.park_hdr; %: {16×1 cell}
d.lr_d=anchor.park_d; %: [180×16 double]

%%  get cal data
get_cpf_cals_mc

%% SET DATA FILL VALUES
fv.bio = 99999;
fv.QC  = 99;

% fv.QC  = 1; % this should fix the QF flag in ODV files

%% 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

% VALID BIO ARGO RANGE CHECKS - QC DATA [MIN MAX]
RC.P     = [0 12000];
RC.S     = [26 38]; % from argo parameter list
RC.T     = [-2.5 40]; % from argo parameter list
RC.O     = [-5 550]; % from argo parameter list
RC.OP    = [10 70]; % optode phase, from argo parameter list
RC.OT    = [-2.5 40]; % optode temperature, from argo parameter list
RC.CHL   = [-0.1 50]; % argoBGC QC manual 09July2016
RC.BB700 = [-0.000025 0.1]; % argoBGC QC manual 09July2016
RC.BB532 = [-0.000005 0.1]; % argoBGC QC manual 09July2016
RC.CDOM  = [-1000 1000]; % Place Holder
RC.NO3   = [-10 55];
RC.PH    = [7.3 8.5];
RC.PHV   = [-1.2 -0.7]; % Range check on pH volts
RC.IB    = [-100 100]; % Range check on pH Ib, nano amps
RC.IK    = [-100 100]; % Range check on pH Ik, nano amps

% INST_ID_str = regexp(MBARI_ID_str,'^\d+', 'once','match');
crazy_val = 99990; %some ridiculous number that is less than the GDAC fill value of 99999

%%
% GATHER INFO FOR BUILDING ODV FILE LATER AND JUST TO MAKE LIFE
% EASIER
INFO.CpActivationP = d.CpActivationP;
INFO.FwRev    = d.FwRev;
INFO.CTDtype  = d.CTDtype;
INFO.CTDsn    = d.CTDsn;
INFO.FlbbMode = d.FlbbMode;
if isfinite(INFO.FlbbMode)  && INFO.FlbbMode >0
    master_FLBB = 1;
end

INFO.sdn      = d.sdn;

%     % CHECK FOR BAD GPS TIME 10/30/19 -jp
%     % WE ARE GETTING TIME FROM PROFILE TERMINATION TIME NOT GPS FIX TIME!!!
%     if ~isempty(first_sdn) && ~isempty(INFO.sdn)
%         dvec = datevec(INFO.sdn);
%         if INFO.sdn > first_sdn + 365*20 && dvec(1) == 2099 % 20 yrs from start?
%             disp(['GPS time for this profile is > 20 years past start ', ...
%                 '- gps week day number bug?!!'])
%             %dvec = datevec(INFO.sdn); % I think this should be moved above 2nd IF? jp 06/20/21
%             dvec(1)  = 1999; % per aoml suggestion don't quite understand jump to 2099
%             INFO.sdn = datenum(dvec) + 1024*7;
%         elseif INFO.sdn < first_sdn % bad gps time fix 10/30/19
%             disp('GPS time for this profile is unreasonable - gps week day number bug!!')
%             disp(['days since first profile = ',num2str((INFO.sdn - ...
%                 first_sdn),'%0.0f')]);
%             INFO.sdn = INFO.sdn + 1024*7;
%         end
%     end
INFO.cast     = d.cast;
INFO.gps      = d.gps;

INFO.float_type='CPF'; % MC, 11/17/2021

%% LOW RES data coolumns

% % 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

% 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

iOCR1 = find(strcmp('OCRCh1_uW/cm2/nm',    d.lr_hdr) == 1); % OCR1
iOCR2 = find(strcmp('OCRCh2',    d.lr_hdr) == 1); % OCR2
iOCR3 = find(strcmp('OCRCh3',    d.lr_hdr) == 1); % OCR3
iOCR4 = find(strcmp('OCRCh4',    d.lr_hdr) == 1); % OCR4

iTime   = find(strcmp('DateTime',      d.lr_hdr) == 1); % ANCHOR data has time associated to it
%% 
% ADD TIME data to the ANCHOR CPF data, MC
LR.TIME      = d.lr_d(:,iTime); % MATLAB Time

%% 
% GET SOME ARGO CORE PARAMTERS
LR.PRES      = d.lr_d(:,iP);

% MAKE AN ARRAY OF ZEROS FOR FILLING ARRAYS LATER
fill0  = ones(size(LR.PRES))* 0;

LR.PRES_QC          = fill0 + fv.QC;
LR.PRES_ADJUSTED    = LR.PRES;
LR.PRES_ADJUSTED_QC = fill0 + fv.QC;

LR.PSAL             = d.lr_d(:,iS);
%LR.PSAL_QC          = fill0 + fv.QC;
LR.PSAL_QC          = fill0 + 1;
LR.PSAL_ADJUSTED    = LR.PSAL;
% LR.PSAL_ADJUSTED_QC = fill0 + fv.QC;
LR.PSAL_ADJUSTED_QC = fill0 + 1;

LR.TEMP             = d.lr_d(:,iT);
% LR.TEMP_QC          = fill0 + fv.QC;
LR.TEMP_QC          = fill0 + 1;
LR.TEMP_ADJUSTED    = LR.TEMP;
% LR.TEMP_ADJUSTED_QC = fill0 + fv.QC;
LR.TEMP_ADJUSTED_QC = fill0 + 1;

% CHECK FOR BAD PRESS VALUES
LRQF_P = LR.PRES < RCR.P(1) | LR.PRES > RCR.P(2);
LR.PRES_QC(LRQF_P)  = 4;  % BAD
LR.PRES_QC(~LRQF_P & LR.PRES ~= fv.bio) = 1; % GOOD

LR.PRES_ADJUSTED_QC(LRQF_P)  = 4;  % BAD
LR.PRES_ADJUSTED_QC(~LRQF_P & LR.PRES_ADJUSTED ~= fv.bio) = 1; % GOOD

% CALCULATE POTENTIAL DENSITY (STILL NEED TO UPDATE TO TEOS 10 !)
potT   = theta(LR.PRES, LR.TEMP, LR.PSAL,0); % potential temp, p=0
% USING POTENTIAL DENSITY ( per: Processing Argo OXYGEN data at the
% DAC level Version 2.0 October 1st 2015, Chap 5)
lr_den = density(LR.PSAL, potT); % potential density, kg /m^3

% % CALCULATE POTENTIAL DENSITY (UPDATED TO USE TEOS 10, by MC)
%     LON=INFO.gps(1); LAT=INFO.gps(2);
%     [SA, in_ocean] = gsw_SA_from_SP(LR.PSAL,LR.PRES,LON,LAT);  % Absolute Salinity from Pactical salinity
%     CT=gsw_CT_from_t(SA,LR.TEMP,LR.PRES);                      % Conservative Temperature
%     lr_den= gsw_rho(SA,CT,0);                                  % In-situ Density, use p_ref=0 to calculate potential density
    
%% ------------------- Calculate CHLA ------------------------------------

% t_nan = isnan(data.lr_d(:,iChl)); % NaN's in data if any
% LR.CHLA= (data.lr_d(~t_nan,iChl) - cal.CHL.ChlDC) .* ...
%                 cal.CHL.ChlScale;

% ****************************************************************
% CALCULATE CHLOROPHYLL CONCENTRATION (µg/L or mg/m^3)
% ****************************************************************

if ~isempty(iChl) % ( ~isempty(iChl) && master_FLBB ~= 0) || ...
    % (~isempty(iChl) && ~isempty(regexp(MBARI_ID_str, bad_chl_filter,'once')))
    
    % Not sure why master_FLBB variable part of the statement, seems having
    % just ~istempty(iChl) would suffice?
    % JP 03/10/20201 - There are a handful of earlier floats where the msg
    % file header indicates that an FLBB exists but there isn't one on the
    % float. For these floats the FLBB_mode = 0 from the start maybe there
    % is a better way...
    
    % Predim variables with fill values then adjust as needed
    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;
    %         INFO.CHLA_SCI_CAL_EQU   = 'not applicable';
    %         INFO.CHLA_SCI_CAL_COEF  = 'not applicable';
    %         INFO.CHLA_SCI_CAL_COM   = 'not applicable';
    %         INFO.CHLA_DATA_MODE  = 'R'; %"not applicable" is not acceptable for this field, should be 'R' (per Coriolis)
    
    t_nan = isnan(d.lr_d(:,iChl)); % NaN's in data if any
    
    LR.FLUORESCENCE_CHLA(~t_nan)    = d.lr_d(~t_nan,iChl);
    %         LR.FLUORESCENCE_CHLA_QC(~t_nan) = fv.QC;
    
    if isfield(cal,'CHL') % Sensor could be bad so maybe no cal info
        LR.CHLA(~t_nan) = (d.lr_d(~t_nan,iChl) - cal.CHL.ChlDC) .* ...
            cal.CHL.ChlScale;
        %             LR.CHLA_QC(~t_nan) =  3; % 3 do not use w/o adjusting
    end
end

%% ------------------- Calculate BBP ------------------------------------

% ****************************************************************
% CALCULATE PARTICLE BACKSCATTER COEFFICIENT FROM VOLUME
% SCATTERING FUNCTION (VSF) (m^-1)
% APEX FLBB
% ****************************************************************
if ~isempty(iBb) % (~isempty(iBb) && master_FLBB ~= 0) || ....
    %     (~isempty(iBb) && ~isempty(regexp(MBARI_ID_str, bad_chl_filter,'once')))
    %     if ~isempty(iBb) && master_FLBB ~= 0
    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;
    %         INFO.BBP700_SCI_CAL_EQU   = 'not applicable';
    %         INFO.BBP700_SCI_CAL_COEF  = 'not applicable';
    %         INFO.BBP700_SCI_CAL_COM   = 'not applicable';
    %         INFO.BBP700_DATA_MODE  = 'R'; %"not applicable" is not acceptable for this field, should be 'R' (per Coriolis)
    
    t_nan = isnan(d.lr_d(:,iBb)); % NaN's in data if any
    
    if isfield(cal,'BB') % Sensor could be bad so maybe no cal info
        % VOLUME SCATTERING FUNCTION (VSF) (m^-1 sr^-1)
        VSF(~t_nan) = (d.lr_d(~t_nan, 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(BETA_SW_ind(ct)), b90sw, bsw] = ...
                    betasw_ZHH2009(LAMBDA,d.lr_d(BETA_SW_ind(ct),iT), ...
                    THETA, d.lr_d(BETA_SW_ind(ct),iS), DELTA);
            end
        end
        %LR.BETA_BACKSCATTERING700 = VSF; % (DOUBLE CHECK SPEC say counts???)
        LR.BETA_BACKSCATTERING700(~t_nan) = d.lr_d(~t_nan,iBb); % counts
        %             LR.BETA_BACKSCATTERING700_QC(~t_nan) = fv.QC; % counts
        LR.BBP700(~t_nan)    = (VSF(~t_nan) - BETA_SW(~t_nan)) * X; %b_bp m^-1
        LR.BBP700_QC(~t_nan) = 1; % 3 do not use w/o adjusting ... 6/10/21 modify qcraw flag from 3 to 2.
        
        %             % CALCULATE ADJUSTED DATA _ NO ADJUSTMENTS AT THIS TIME
        %             % BUT...6/10/21 START POPULATING BBP700_ADJUSTED WITH BBP
        %             % DIRECTLY!!!  KEEP ORIGINAL BLOCK IN CASE AN ACTUAL ADJUSTMENT
        %             % GETS IMPLEMENTED/STORED IN THE QC MATRIX.
        %             LR.BBP700_ADJUSTED(~t_nan) = LR.BBP700(~t_nan);
        %             LR.BBP700_ADJUSTED_QC(~t_nan) = 1;
        %             LR.BBP700_ADJUSTED_ERROR(~t_nan) = fv.bio; % PLACE HOLDER FOR NOW
        %             INFO.BBP700_SCI_CAL_EQU  = 'BBP700_ADJUSTED=BBP700';
        %             INFO.BBP700_SCI_CAL_COEF = [''];
        %             INFO.BBP700_SCI_CAL_COM  =['BBP700_ADJUSTED is being filled with BBP700 directly in real time.  Adjustment method may be enhanced in the future.'];
        %
        %             %             if isfield(QC,'BB')
        %             %                 QCD = [LR.PRES, LR.TEMP, LR.PSAL, LR.BBP700];
        %             %                 LR.BBP700_ADJUSTED(~t_nan) = ...
        %             %                     apply_QC_corr(QCD(~t_nan,:), INFO.sdn, QC.BB);
        %             %                 LR.BBP700_ADJUSTED_QC(~t_nan) =  2;
        %             %                 LR.BBP700_ADJUSTED_ERROR(~t_nan) = fv.bio; % PLACE HOLDER FOR NOW
        %             %                 INFO.BBP700_SCI_CAL_EQU  = 'BBP700_ADJUSTED=BBP700*A-B';
        %             %                 INFO.BBP700_SCI_CAL_COEF = ['A=', ...
        %             %                     num2str(QC.BB.steps(3),'%0.4f'),',B=',...
        %             %                     num2str(QC.BB.steps(4),'%0.4f')];
        %             %                 INFO.BBP700_SCI_CAL_COM  =['A and B determined by comparison', ...
        %             %                     ' to discrete samples from post deployment calibration',...
        %             %                     ' rosette cast'];
        %             %             end
    end
end

%% ------------------- 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 = ?

iPhase=iTPh; % check w/Yui

% ****************************************************************
% CALCULATE OXYGEN CONCENTRATION
% ****************************************************************
if ~isempty(iPhase) % Bphase or Tphase exists so O2 data should exist
    t_nan = isnan(d.lr_d(:,iPhase));
    O2_chk = d.lr_d(:,iPhase) < RCR.OP(1) | ...
        d.lr_d(:,iPhase) > RCR.OP(2) | ...
        d.lr_d(:,iTo) < RCR.OT(1) | ...
        d.lr_d(:,iTo) > RCR.OT(2); % Flag phase, Temp out of range
    if any(O2_chk)
        disp(['Out of range phase or optode T detected for ', ...
            msg_file])
    end
    
    if strcmp(cal.O.type,'3830')
        LR.BPHASE_DOXY            = fill0 + fv.bio; % predim w fill val
        LR.BPHASE_DOXY_QC         = fill0 + fv.QC;
        LR.BPHASE_DOXY(~t_nan)    = d.lr_d(~t_nan,iBPh);
        LR.BPHASE_DOXY_QC(~t_nan) = fv.QC;
    elseif strcmp(cal.O.type,'4330')
        LR.TPHASE_DOXY            = fill0 + fv.bio;
                     LR.TPHASE_DOXY_QC         = fill0 + fv.QC;
        LR.TPHASE_DOXY(~t_nan)    = d.lr_d(~t_nan,iTPh);
                     LR.TPHASE_DOXY_QC(~t_nan) =  fv.QC;
    elseif strcmp(cal.O.type,'SBE83')
        LR.PHASE_DELAY_DOXY            = fill0 + fv.bio;
        LR.PHASE_DELAY_DOXY_QC         = fill0 + fv.QC;
        LR.PHASE_DELAY_DOXY(~t_nan)    = d.lr_d(~t_nan,iPh83);
        LR.PHASE_DELAY_DOXY_QC(~t_nan) =  fv.QC;
    end
    LR.DOXY                = fill0 + fv.bio;
    LR.DOXY_QC             = fill0 + fv.QC;
    LR.TEMP_DOXY           = fill0 + fv.bio;
    LR.TEMP_DOXY_QC        = fill0 + fv.QC;
    %         LR.DOXY_ADJUSTED       = fill0 + fv.bio;
    %         LR.DOXY_ADJUSTED_QC    = fill0 + fv.QC;
    %         LR.DOXY_ADJUSTED_ERROR = fill0 + fv.bio;
    %         INFO.DOXY_SCI_CAL_EQU  = 'not applicable';
    %         INFO.DOXY_SCI_CAL_COEF = 'not applicable';
    %         INFO.DOXY_SCI_CAL_COM  = 'not applicable';
    %         INFO.DOXY_DATA_MODE  = 'R'; %"not applicable" is not acceptable for this field, should be 'R' (per Coriolis)
    
    if strcmp(cal.O.type,'SBE83')
        myOdata = [LR.PRES(~t_nan) LR.TEMP(~t_nan) LR.PSAL(~t_nan) lr_d((~t_nan),iPhase),lr_d((~t_nan),iTo)];
        [ppoxdoxy, O2_uM, O2_T] = Calc_SBE63_O2(myOdata, cal.O);
        LR.DOXY(~t_nan) = O2_uM ./ lr_den(~t_nan) .* 1000; % µmol/kg
    else
        %[S, T, P, Phase, CalPhase, [O2], O2sol, pO2]
        O2 = calc_O2_4ARGO(LR.PSAL(~t_nan), LR.TEMP(~t_nan), ...
            LR.PRES(~t_nan),d.lr_d((~t_nan),iPhase), cal.O); % O2 in µmol/L + more
        LR.DOXY(~t_nan) = O2(:,6) ./ lr_den(~t_nan) .* 1000; % µmol/kg
    end
    tDOXY = abs(LR.DOXY) > crazy_val & ~t_nan; % Unrealistic bad value
    LR.DOXY(tDOXY) = crazy_val; % SET TO crazy bad value
    LR.DOXY_QC(~t_nan)      = 3;
    % JP: QC assignment below not needed - crazy value greater than range limits
    % LR.DOXY_QC(tDOXY) = 4; % set crazy bad value QF to 4
    LR.TEMP_DOXY(~t_nan)    = d.lr_d(~t_nan,iTo);
    LR.TEMP_DOXY_QC(~t_nan) = fv.QC;
    % Save O2sol, pO2?
    clear O2
end

%% AIR OXYGEN

if ~isempty(iPhase) 
    
end


%% ------------------- 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([dirs.temp,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 NITRATE (µmol / kg scale)
% DO DEPTH CORRECTION 1st
% CONVERT TO µmol/kg
% ****************************************************************
if ~isempty(iNO3)
    LR.NITRATE                      = fill0 + fv.bio;
    LR.NITRATE_QC                   = fill0 + fv.QC;
    LR.UV_INTENSITY_DARK_NITRATE    = fill0 + fv.bio;
    LR.UV_INTENSITY_DARK_NITRATE_QC = fill0 + fv.QC;
    LR.UV_INTENSITY_NITRATE         = fill0 + fv.bio;
    LR.UV_INTENSITY_NITRATE_QC      = fill0 + fv.QC;
    
    LR.NITRATE_ADJUSTED             = fill0 + fv.bio;
    LR.NITRATE_ADJUSTED_QC          = fill0 + fv.QC;
    LR.NITRATE_ADJUSTED_ERROR       = fill0 + fv.bio;
    INFO.NITRATE_SCI_CAL_EQU  = 'not applicable';
    INFO.NITRATE_SCI_CAL_COEF = 'not applicable';
    INFO.NITRATE_SCI_CAL_COM  = 'not applicable';
    INFO.NITRATE_DATA_MODE = 'not applicable';
    
    % IF ISUS FILE & CAL INFO EXIST, TRY AND PROCESS IT
    if exist([dirs.temp, NO3_file],'file') && isfield(cal,'N')
        spec = parse_NO3msg([dirs.temp,NO3_file]); % return struct
        UV_INTEN = spec.UV_INTEN;
        if ~isempty(UV_INTEN)
            % [SDN, DarkCur, Pres, Temp, Sal, NO3, BL_int,BL_slope,
            %  RMS_ER, Wl~240, ABS~240] !!! NITRATE STILL µmol/L !!!
            
            NO3  = calc_FLOAT_NO3(spec, cal.N, 1); % ESW P corr
            %NO3  = calc_APEX_NO3_JP(spec, cal.N, 0); % NO P corr
            
            IX = (size(NO3,1):-1:1)'; % FLIP SHALLOW TO DEEP FOR ARGO
            %[B,IX]   = sort(NO3(:,3)); % SORT SHALLOW TO DEEP FOR ARGO
            NO3      = NO3(IX,:); % µmol/L
            UV_INTEN = UV_INTEN(IX,:);
            clear B IX
        end
        
    elseif ~isfield(cal,'N')
        disp('No calibration info to process nitrate with')
    else
        disp([dirs.temp, NO3_file,' not found!!'])
    end
    
    if exist('NO3', 'var') == 1
        % MATCH NITRATE DEPTHS TO CTD DEPTHS IF LENGTHS DIFFERENT
        % MISSING SPECTRA  AND NO3 CONC. SET TO NaN
        % could also  replace missing NO3 from *.isus rwith *.msg
        
        if size(LR.PRES,1) ~= size(NO3,1) % Mismatch in sample lengths
            cS       = size(spec.UV_INTEN,2); % cols of WL's
            rP       = size(LR.PRES,1); % # of CTD samples (*.msg)
            [rN, cN] = size(NO3); % RN = # of isus samples (*.isus)
            disp(['*.isus(',num2str(rN),') & *.msg(',num2str(rP), ...
                ') sample counts not equal for ',msg_file])
            NO3tmp   = ones(rP,cN)* NaN;
            spectmp  = ones(rP,cS)* fv.bio;
            % fill msg file data then replace with NO3 data
            %NO3tmp(:,3:6) = [LR.PRES, LR.TEMP, LR.PSAL lr_d(:,iNO3)];
            NO3tmp(:,3:5) = [LR.PRES, LR.TEMP, LR.PSAL];
            
            for i = 1: rP
                ind  = find(NO3(:,3) == LR.PRES(i),1,'first');
                if ~isempty(ind)
                    NO3tmp(i,:)  = NO3(ind,:);
                    spectmp(i,:) = UV_INTEN(ind,:);
                end
            end
            NO3      = NO3tmp;
            
            UV_INTEN = spectmp; % n x m matrix
            clear rP rN cN NO3tmp i
        end
        
        t_nan  = isnan(NO3(:,6));
        % CHECK FIT ERROR AND BASELINE ABSORBANCE
        tABS08 = NO3(:,11) > 0.8; %ABS240 > 0.8 (QF=3)
        tABS11 = NO3(:,9) > 0.003 | NO3(:,11) > 1.1; % RMS & ABS240 (QF=4)
        
        LR.UV_INTENSITY_DARK_NITRATE    = NO3(:,2);
        LR.UV_INTENSITY_DARK_NITRATE(t_nan) = fv.bio;
        LR.UV_INTENSITY_DARK_NITRATE_QC = fill0 + fv.QC;
        
        LR.UV_INTENSITY_NITRATE    =  UV_INTEN;
        LR.UV_INTENSITY_NITRATE_QC =  UV_INTEN * 0 + fv.QC;
        
        NpotT  = theta(NO3(:,3), NO3(:,4), NO3(:,5),0);
        N_den  = density(NO3(:,5), NpotT);
        %NO3_kg = NO3(:,6) ./ N_den * 1000;
        
        
        % ********************************************************
        % CALCULATE NITRATE IN umol/L FOR PRESSURE = 0
        % (vs in situ p). ISUS calibrated at p = 0 but
        % ISUS MEASURES IN SITU CONCENTRATION SO ISUS NO3 SLIGHTLY
        % GREATER THAN BOTTLE or SURFACE NO3 (p=0) AT DEPTH
        % (more moles NO3 per liter at depth)
        % NO3(p=0) = N03(p=z) *[density(S,T,p=0) /density(S,T,p@z)]
        
        NO3_p0 = NO3(:,6) .* sw_dens(NO3(:,5),NpotT,0) ./ ...
            sw_dens(NO3(:,5),NO3(:,4),NO3(:,3));
        
        NO3_p0_kg         = NO3_p0 ./ N_den * 1000; % µmol/kg
        LR.NITRATE        = NO3_p0_kg;
        LR.NITRATE(t_nan) = fv.bio;
        
        % QC RAW
        LR.NITRATE_QC = fill0 + 3; % questionable to start
        LR.NITRATE_QC(t_nan) = fv.QC;
        LR.NITRATE_QC(~t_nan & tABS11) = 4;
        
%         %             % ********************************************************
%         %             % APPLY QC CORRECTIONS
%         %             % CORRECTIONS DETERMINED ON µmol/L scale so adjust on that
%         %             % scale and then convert
%         %             if isfield(cal,'N') && isfield(QC,'N')
%         %                 QCD = [LR.PRES, LR.TEMP, LR.PSAL,NO3_p0_kg];
%         %                 LR.NITRATE_ADJUSTED = fill0 + fv.bio;
%         %                 LR.NITRATE_ADJUSTED(~t_nan) = ...
%         %                     apply_QC_corr(QCD(~t_nan,:), INFO.sdn, QC.N);
%         %
%         %                 LR.NITRATE_ADJUSTED_QC = fill0 + 1; % set=1 9/27/16
%         %                 LR.NITRATE_ADJUSTED_QC(t_nan) = fv.QC; % nan's to fill value
%         %
%         %                 LR.NITRATE_ADJUSTED_QC(~t_nan & tABS08) = 3; %Bad ABS240
%         %                 if sum(tABS11) > 0
%         %                     disp(['Nitrate ABS@240 > 1.1 or RMS > 0.003 ',...
%         %                         'detected. Setting QF = 4 for flagged values'])
%         %                     LR.NITRATE_ADJUSTED_QC(~t_nan & tABS11) = 4; % Really bad abs or RMS
%         %                 end
%         %
%         %                 LR.NITRATE_ADJUSTED_ERROR = (abs(LR.NITRATE - ...
%         %                     LR.NITRATE_ADJUSTED)) * 0.1 + 0.5;
%         %                 LR.NITRATE_ADJUSTED_ERROR(t_nan) = fv.bio;
%         %
%         %
%         %                 INFO.NITRATE_SCI_CAL_EQU  = ['NITRATE_ADJUSTED=', ...
%         %                     '[NITRATE-SUM(OFFSET(S)+DRIFT(S))]/GAIN'];
%         %                 INFO.NITRATE_SCI_CAL_COEF = ['OFFSET(S) and DRIFT(S) ', ...
%         %                     'from climatology comparisons at 1000m or 1500m. GAIN ',...
%         %                     'from surface/deep comparison where surface values ',...
%         %                     'are known'];
%         %                 INFO.NITRATE_SCI_CAL_COM  =['Contact Tanya Maurer ',...
%         %                     '(tmaurer@mbari.org) or Josh Plant (jplant@mbari.org) ',...
%         %                     'for more information'];
%         %             end
%         %             clear QCD UV_INTEN
%         %         end
%         %
%         %         %         % DO A FINAL RANGE CHECK ON VALUES, IF BAD SET QF = 4
%         %         %         % RAW
%         %         %         [BSLflag, theflag] = isbadsensor(BSL, MBARI_ID_str, INFO.cast, 'N');
%         %         %         t_bio   = LR.NITRATE ~= fv.bio;
%         %         %         %tST     = LR.PSAL_QC == 4 | LR.TEMP_QC == 4; % Bad S or T will affect nitrate
%         %         %         tST     = LR.PSAL_QC == 4 | LR.TEMP_QC == 4 | LR.PRES_QC == 4; % Bad STP will affect nitrate
%         %         %         LR.NITRATE_QC(t_bio) = LR.NITRATE_QC(t_bio) * ~BSLflag + BSLflag*theflag;
%         %         %         LR.UV_INTENSITY_DARK_NITRATE_QC(t_bio) = ...
%         %         %             LR.UV_INTENSITY_DARK_NITRATE_QC(t_bio) * ~BSLflag + BSLflag*theflag;
%         %         %
%         %         %         t_chk1 = t_bio &(LR.NITRATE < RCR.NO3(1)| LR.NITRATE > RCR.NO3(2)); % RANGE
%         %         %         t_chk2 = t_bio & tST; % SALT
%         %         %         LR.NITRATE_QC(t_chk1 | t_chk2) = 4;
%         %         %         LR.UV_INTENSITY_DARK_NITRATE_QC(t_chk1)  = 4;
%         %
%         %         % DO A FINAL RANGE CHECK BASED ON BSL & OTHER DEPENDANT PARAMETERS
%         %         % MODIFIED 04/07/2020 BY JP BECAUSE BSL FLAG OF 3 COULD OVERWRITE A 4
%         %         t_bio = LR.NITRATE ~= fv.bio; % non fill values
%         %         [BSLflag, theflag] = isbadsensor(BSL, MBARI_ID_str, INFO.cast, 'N');
%         %         tBSL  = ones(size(t_bio)) * ~BSLflag + BSLflag*theflag; % flag from BSL
%         %         tSTP  = (LR.PSAL_QC == 4 | LR.TEMP_QC == 4 | LR.PRES_QC == 4)*4; % Bad STP will affect nitrate
%         %         tRC   = (LR.NITRATE < RCR.NO3(1)| LR.NITRATE > RCR.NO3(2))*4; % RAnge check
%         %
%         %         tALL  = max([LR.NITRATE_QC, tBSL, tSTP, tRC],[],2); % get highest flag
%         %         LR.NITRATE_QC(t_bio) = tALL(t_bio);
%         %
%         %         tALL  = max([LR.UV_INTENSITY_DARK_NITRATE_QC, tBSL, tRC],[],2); % get highest flag
%         %         LR.UV_INTENSITY_DARK_NITRATE_QC(t_bio) = tALL(t_bio);
%         %
%         %
%         %
%         %         % *****  ADJUSTED  *****
%         %         %         t_bio   = LR.NITRATE_ADJUSTED ~= fv.bio;
%         %         %         %tST     = LR.PSAL_ADJUSTED_QC == 4 | LR.TEMP_ADJUSTED_QC == 4; % Bad S or T will affect nitrate
%         %         %         tST     = LR.PSAL_ADJUSTED_QC == 4 | LR.TEMP_ADJUSTED_QC == 4 | ...
%         %         %                   LR.PRES_ADJUSTED_QC == 4;
%         %         %         LR.NITRATE_ADJUSTED_QC(t_bio) = LR.NITRATE_ADJUSTED_QC(t_bio) ...
%         %         %             * ~BSLflag + BSLflag*theflag;
%         %         %
%         %         %         t_chk1 = t_bio & (LR.NITRATE_ADJUSTED < RC.NO3(1)| ...
%         %         %                  LR.NITRATE_ADJUSTED > RC.NO3(2));
%         %         %         t_chk2 = t_bio & tST;
%         %         %         LR.NITRATE_ADJUSTED_QC(t_chk1 | t_chk2) = 4;
%         %
%         %         t_bio = LR.NITRATE_ADJUSTED ~= fv.bio;
%         %         tBSL  = ones(size(t_bio)) * ~BSLflag + BSLflag*theflag; % flag from BSL
%         %         tSTP  = (LR.PSAL_ADJUSTED_QC == 4 | LR.TEMP_ADJUSTED_QC == 4 | ...
%         %             LR.PRES_ADJUSTED_QC == 4)*4;
%         %         tRC   = (LR.NITRATE_ADJUSTED < RC.NO3(1)| LR.NITRATE_ADJUSTED > RC.NO3(2))*4;
%         %
%         %         tALL  = max([LR.NITRATE_ADJUSTED_QC, tBSL, tSTP, tRC],[],2); % get highest flag
%         %         LR.NITRATE_ADJUSTED_QC(t_bio) = tALL(t_bio);
%         %
%         %         % double check fill value QC flags - 9634 returns fv's from qc_adjustmet
%         %         % due to bad gps date bug - jp 10/04/19
%         %         LR.NITRATE_ADJUSTED_QC(~t_bio)  = fv.QC;
%         %
%         %         % -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
%         %         % DO A SPIKE TEST ON VALUES, IF SPIKES ARE IDENTIFIED, SET QF TO 4
%         %         % (BAD).  BGC_spiketest screens for nans and fill values.
%         %         %
%         %         % RUN TEST ON RAW NITRATE
%         %         QCscreen_N = LR.NITRATE_QC == 4; % screen for BAD data already assigned.
%         %         [spike_inds, quality_flags] = BGC_spiketest(MBARI_ID_str,str2double(cast_num),[LR.PRES LR.NITRATE],'NO3',dirs.cal,fv.bio,QCscreen_N);
%         %         if ~isempty(spike_inds)
%         %             LR.NITRATE_QC(spike_inds) = quality_flags;
%         %             disp(['LR.NITRATE QUALITY FLAGS ADJUSTED FOR IDENTIFIED SPIKES, PROFILE ',cast_num,' FLOAT ',MBARI_ID_str,'.'])
%         %             %         else
%         %             %             disp(['NO LR.NITRATE SPIKES IDENTIFIED FOR PROFILE ',cast_num,' FLOAT ',MBARI_ID_str,'.'])
%         %         end
%         %         %
%         %         % RUN TEST ON QC NITRATE_ADJUSTED
%         %         QCscreen_Nadj = LR.NITRATE_ADJUSTED_QC == 4; % screen for BAD data already assigned.
%         %         [spike_inds, quality_flags] = BGC_spiketest(MBARI_ID_str,str2double(cast_num),[LR.PRES LR.NITRATE_ADJUSTED],'NO3',dirs.cal,fv.bio,QCscreen_Nadj);
%         %         if ~isempty(spike_inds)
%         %             LR.NITRATE_ADJUSTED_QC(spike_inds) = quality_flags;
%         %             disp(['LR.NITRATE_ADJUSTED QUALITY FLAGS ADJUSTED FOR IDENTIFIED SPIKES, PROFILE ',cast_num,' FLOAT ',MBARI_ID_str,'.'])
%         %             %         else
%         %             %             disp(['NO LR.NITRATE_ADJUSTED SPIKES IDENTIFIED FOR PROFILE ',cast_num,' FLOAT ',MBARI_ID_str,'.'])
    end
    % -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
    clear tchk tABS08 tABS11 NO3 QCscreen_N QCscreen_Nadj spike_inds quality_flags
end

%% ------------------- Calculate pH -------------------------------------
% main functions: parse_pHmsg([data_dir,dura_fns(i).name]);, phcalc

% ****************************************************************
% CALCULATE pH (µmol / kg scale)
% ****************************************************************
if ~isempty(ipH)
    LR.VRS_PH                          = fill0 + fv.bio;
    LR.VRS_PH_QC                       = fill0 + fv.QC;
    LR.PH_IN_SITU_FREE                 = fill0 + fv.bio;
    LR.PH_IN_SITU_FREE_QC              = fill0 + fv.QC;
    LR.PH_IN_SITU_TOTAL                = fill0 + fv.bio;
    LR.PH_IN_SITU_TOTAL_QC             = fill0 + fv.QC;
    LR.IB_PH                           = fill0 + fv.bio;
    LR.IB_PH_QC                        = fill0 + fv.QC;
    IK                                 = fill0 + fv.bio;
    LR.PH_IN_SITU_TOTAL_ADJUSTED       = fill0 + fv.bio;
    LR.PH_IN_SITU_TOTAL_ADJUSTED_QC    = fill0 + fv.QC;
    LR.PH_IN_SITU_TOTAL_ADJUSTED_ERROR = fill0 + fv.bio;
    INFO.PH_SCI_CAL_EQU  = 'not applicable';
    INFO.PH_SCI_CAL_COEF = 'not applicable';
    INFO.PH_SCI_CAL_COM  = 'not applicable';
    INFO.PH_DATA_MODE  = 'R'; %"not applicable" is not acceptable for this field, should be 'R' (per Coriolis)
    
    t_nan = isnan(d.lr_d(:,ipH)); %NaN's?
    
    LR.VRS_PH(~t_nan)      = d.lr_d(~t_nan,ipH); % I param
    LR.VRS_PH_QC(~t_nan)   = fv.QC;
    if isfield(cal,'pH')
        [phfree,phtot] = phcalc(LR.VRS_PH(~t_nan), ...
            LR.PRES(~t_nan), LR.TEMP(~t_nan), LR.PSAL(~t_nan), ...
            cal.pH.k0, cal.pH.k2, cal.pH.pcoefs);
        LR.PH_IN_SITU_FREE(~t_nan)    = phfree; % I param
        LR.PH_IN_SITU_FREE_QC(~t_nan) = fv.QC;
        
        LR.PH_IN_SITU_TOTAL(~t_nan)    = phtot;
        LR.PH_IN_SITU_TOTAL_QC(~t_nan) = 3;
        
        LR_inf = isinf(LR.PH_IN_SITU_FREE); % happens if S = 0
        LR.PH_IN_SITU_FREE(LR_inf)    = 20.1; %UNREAL #
        LR.PH_IN_SITU_FREE_QC(LR_inf) = 4;
        LR.PH_IN_SITU_TOTAL(LR_inf)    = 20.1; %UNREAL #
        LR.PH_IN_SITU_TOTAL_QC(LR_inf) = 4;
    else % data returned  but no cal
        disp(['pH data exists but no cal file for ',INST_ID_str])
    end
    
    % GET  DATA FROM *.dura SO I CAN GRAB Ib & Ik
    dura = parse_pHmsg([dirs.temp,pH_file]);
    myINFO = dir([dirs.temp,pH_file]);
    %         dura = parse_pHmsg([pH_file]); % MC, pH_file has directory in it
    %         myINFO = dir([pH_file]);
    if ~isempty(myINFO)
        ph_filesize = myINFO.bytes;
    else
        ph_filesize = 0;
    end
    if (isempty(dura.data) || ph_filesize< 10000) % usually > 10 Kb for full profile
        disp(['Not enough pH data returned for ',pH_file])
    else
        pH_data = dura.data;
        pH_hdr  = dura.hdr;
        % GET some pH column indices
        ipH_p  = find(strcmp('CTD Pres', pH_hdr) == 1);
        ipH_t  = find(strcmp('CTD Temp', pH_hdr) == 1);
        ipH_s  = find(strcmp('CTD Sal', pH_hdr) == 1);
        ipH_Ib = find(strcmp('Ib', pH_hdr) == 1); % Base current
        ipH_Ik = find(strcmp('Ik', pH_hdr) == 1); % Base current
        
        ph_rows = size(pH_data,1);
        IX      = (ph_rows:-1:1)';
        %[B,IX]  = sort(pH_data(:,ipH_p));
        pH_data = pH_data(IX,:);
        clear B IX ph_rows
        
        % Mismatch in sample lengths
        if size(LR.PRES,1) ~= size(pH_data(:,1),1)
            rP       = size(LR.PRES,1); % # of CTD samples (*.msg)
            [rN, cN] = size(pH_data);       % size pH data
            disp(['*.dura(',num2str(rN),') & *.msg(',num2str(rP), ...
                ') sample counts not equal for ',msg_file])
            
            pHtmp    = ones(rP,cN)* fv.bio; % start with fill values
            % fill msg file data then replace with pH data
            pHtmp(:,[ipH_p,ipH_t,ipH_s]) = [LR.PRES, LR.TEMP, LR.PSAL];
            
            for i = 1: rP
                ind  = find(pH_data(:,ipH_p) == LR.PRES(i),1,'first');
                if ~isempty(ind)
                    pHtmp(i,:)  = pH_data(ind,:);
                end
            end
            pH_data = pHtmp;
            clear rP rN cN NO3tmp i pHtmp
        end
        LR.IB_PH(~t_nan)    = pH_data(~t_nan,ipH_Ib) * 1e9; % nano amps
        LR.IB_PH_QC(~t_nan) = fv.QC;
        IK(~t_nan) = pH_data(~t_nan, ipH_Ik) * 1e9; % nano amps, local not for ARGO
        clear ipH_p ipH_t ipH_t ipH_Ib pH_data pH_hdr
    end
    
    %         if isfield(cal,'pH') && isfield(QC,'pH')
    %             QCD = [LR.PRES(~t_nan), LR.TEMP(~t_nan), LR.PSAL(~t_nan), phtot];
    %             LR.PH_IN_SITU_TOTAL_ADJUSTED(~t_nan) = ...
    %                 apply_QC_corr(QCD, INFO.sdn, QC.pH);
    %             LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(~t_nan) = 1; % set=1 9/27/16 vs
    %             LR.PH_IN_SITU_TOTAL_ADJUSTED_ERROR(~t_nan) = 0.01;
    %
    %             LR.PH_IN_SITU_TOTAL_ADJUSTED(LR_inf)    = 20.1; %UNREAL #
    %             LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(LR_inf) = 4;
    %
    %             INFO.PH_SCI_CAL_EQU  = ['PH_ADJUSTED=PH+[PUMP_OFFSET', ...
    %                 '-SUM(OFFSET(S)+DRIFT(S))]* TCOR'];
    %             INFO.PH_SCI_CAL_COEF = ['PUMP_OFFSET=Pump may add ', ...
    %                 'interfernece in CP mode,OFFSET(S) and DRIFT(S) from ',...
    %                 'climatology comparisons at 1000m or 1500m,','TCOR=', ...
    %                 '(TREF+273.15)./(T+273.15).  TREF = TEMP at 1500m.'];
    %             INFO.PH_SCI_CAL_COM  =['Contact Tanya Maurer ',...
    %                 '(tmaurer@mbari.org) or Josh Plant (jplant@mbari.org) ',...
    %                 'for more information'];
    %
    %             % TEMPORARY ADJUSTED pH FIX 08/02/2016
    %             % FLOATVIZ pH CALCULATED WITH OLDER FUNCTION. QC STEPS
    %             % DETERMINED WITH OLD pH VALUES, BUT A CONSTANT OFFSET
    %             % JP pH - FV pH = 0.0167)
    %             % Commented out 9/28/16 doing Qc on JP files now
    %             %                 disp(['!!!! APPLYING TEMPORARY pH CORRECTION TO ADJUSTED',...
    %             %                     ' VALUES (adj_pH = adj_pH - 0.0167)']);
    %             %                 LR.PH_IN_SITU_TOTAL_ADJUSTED(~t_nan) = ...
    %             %                     LR.PH_IN_SITU_TOTAL_ADJUSTED(~t_nan) - 0.0167;
    %
    %             % FIX 9660 04/02/18 -jp
    %             if strcmpi(MBARI_ID_str, 'ua9660')
    %                 % pH sensor is not in pumped flow stream. Sensor is exposed to
    %                 % ambient light. Use Ib to correct for light sensitivity
    %                 %d(pH)/d(iB), daylight, <50m, Model I R*R > 0.8, = 1.488e6
    %                 gps9660 = mean(INFO.gps,1);
    %                 [~,El] = SolarAzElq(gps9660(1), gps9660(3),gps9660(2), 0);% sdn lat lon al
    %                 if El > 0
    %                     dz1 = LR.PRES <= 50;
    %                     dz2 = LR.PRES <= 120 & LR.PRES >= 70;
    %                     IBo = median(LR.IB_PH(dz2),'omitnan');
    %                     LR.PH_IN_SITU_TOTAL_ADJUSTED(~t_nan & dz1) = ...
    %                         LR.PH_IN_SITU_TOTAL_ADJUSTED(~t_nan & dz1) + ...
    %                         (LR.IB_PH(~t_nan & dz1) - IBo)./1e9 * 1.488e6;
    %                     LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(~t_nan & dz1) = 2;
    %                     clear dz1 dz2 IBo gps9660
    %                 end
    %             end
    %         end
    %
    %         % DO A FINAL RANGE CHECK ON VALUES, IF BAD SET QF = 4
    %         % RAW DATA
    %         %tST     = LR.PSAL_QC == 4 | LR.TEMP_QC == 4; % Bad S or T will affect pH
    %         tST     = LR.PSAL_QC == 4 | LR.TEMP_QC == 4 | LR.PRES_QC == 4;
    %         t_bio   = LR.PH_IN_SITU_TOTAL ~= fv.bio;
    %
    %         t_chk1 = t_bio & (LR.PH_IN_SITU_TOTAL < RCR.PH(1)| ...
    %                  LR.PH_IN_SITU_TOTAL > RCR.PH(2) | tST); % RANGE
    %         t_chk2 = t_bio & LR.IB_PH ~= fv.bio.*1e9 & LR.IB_PH ~= fv.bio ...
    %             & (LR.IB_PH < RCR.IB(1) | LR.IB_PH > RCR.IB(2) | ...
    %             IK < RCR.IK(1) | IK > RCR.IK(2)); % DIAGNOSTIC
    %         %t_chk3 = t_bio & (LRQF_S | LRQF_T);
    %         t_chk4 = LR.VRS_PH ~= fv.bio &  (LR.VRS_PH < RCR.PHV(1) | ...
    %             LR.VRS_PH > RCR.PHV(2));
    %         t_chk5 = t_bio & LR.IB_PH ~= fv.bio.*1e9 & LR.IB_PH ~= fv.bio  ...
    %             & (LR.IB_PH < RCR.IB(1).*25 | LR.IB_PH > RCR.IB(2).*25 | ...
    %             IK < RCR.IK(1).*25 | IK > RCR.IK(2).*25); % DIAGNOSTIC
    %
    %         %LR.PH_IN_SITU_TOTAL_QC(t_chk1 | t_chk2 | t_chk3) = 4;
    %         LR.PH_IN_SITU_TOTAL_QC(t_chk2) = 3; % set pH with out of range Ik/Ib to questionable, unless pH is also out of range (then set to bad)
    %         LR.PH_IN_SITU_TOTAL_QC(t_chk5) = 4;
    % 		LR.PH_IN_SITU_TOTAL_QC(t_chk1) = 4;
    %         LR.VRS_PH_QC(t_chk4) = 4;
    %         LR.IB_PH_QC(t_chk2)  = 3;
    %         LR.IB_PH_QC(t_chk5)  = 4;
    %
    % 		[BSLflag, theflag] = isbadsensor(BSL, MBARI_ID_str, INFO.cast, 'PH');
    %         LR.VRS_PH_QC(t_bio) = LR.VRS_PH_QC(t_bio) * ~BSLflag + BSLflag*theflag;
    %         LR.IB_PH_QC(t_bio)  = LR.IB_PH_QC(t_bio) * ~BSLflag + BSLflag*theflag;
    %         LR.PH_IN_SITU_FREE_QC(t_bio) = LR.PH_IN_SITU_FREE_QC(t_bio) ...
    %             * ~BSLflag + BSLflag*theflag;
    %         LR.PH_IN_SITU_TOTAL_QC(t_bio) = LR.PH_IN_SITU_TOTAL_QC(t_bio) ...
    %             * ~BSLflag + BSLflag*theflag;
    
    
    %         % DO A FINAL QC CHECK
    %         % 04/07/2020 JP - Now that the BSL flag can be 3 or 4 it has the
    %         % potential to overwrite a 4 with a 3 so don't let this happen
    %         t_bio = LR.PH_IN_SITU_TOTAL ~= fv.bio;
    %
    %         %6/11/20 TM need this incase no diagnostic data, 0 if no data.
    %         % added the ".*1e9" to this line.  It was lacking in Josh's update
    %         % from 4/7/20 and causing erroneous flagging for certain cases, ie
    %         % float 9634 surface samples of cycles 67 and 97.
    %         %t_diag = LR.IB_PH ~= fv.bio.*1e9;
    %
    %         % JP 12/17/20 slight modification to TM's fix. If the dura file is
    %         % too small no Ib & Ik tests happen & LR.IB_PH, LR.IK_PH is never
    %         % reassigned so you need to check for fv.bio*1e9 & fv.bio as fill
    %         % values
    %         t_diag = ~(LR.IB_PH == fv.bio.*1e9 | LR.IB_PH == fv.bio); % JP 12/17/20
    %
    %         [BSLflag, theflag] = isbadsensor(BSL, MBARI_ID_str, INFO.cast, 'PH');
    %         tBSL  = ones(size(t_bio)) * ~BSLflag + BSLflag*theflag; % flag from BSL
    %         tSTP  = (LR.PSAL_QC == 4 | LR.TEMP_QC == 4 | LR.PRES_QC == 4)*4; % Bad STP will affect nitrate
    %         tRC   = (LR.PH_IN_SITU_TOTAL < RCR.PH(1)|LR.PH_IN_SITU_TOTAL > RCR.PH(2))*4; % Range check
    %         tVRS  = (LR.VRS_PH < RCR.PHV(1) | LR.VRS_PH > RCR.PHV(2))*4;
    %         tIB3  = (LR.IB_PH < RCR.IB(1) | LR.IB_PH > RCR.IB(2)).*t_diag*3; % DIAGNOSTIC
    %         tIK3  = (IK < RCR.IK(1) | IK > RCR.IK(2)).*t_diag*3; % DIAGNOSTIC
    %         tIB4  = (LR.IB_PH < RCR.IB(1)*25 | LR.IB_PH > RCR.IB(2)*25).*t_diag*4; % DIAGNOSTIC
    %         tIK4  = (IK < RCR.IK(1)*25 | IK > RCR.IK(2)*25).*t_diag*4; % DIAGNOSTIC
    %
    %         % 09/30/20 TM - EXCLUSION BLOCK FOR NEWER EQPAC FLOATS EXHIBITING
    %         % THE ERRONEOUS RAILED DIAG VALUES
    %         if regexp(MBARI_ID_str, skip_ph_diag, 'once')
    %             % keep syntax the same to ensure proper function!
    %             % Set flags in this block to 0 (always good)
    %             disp('Excluding float from pH sensor diagnostic checks (erroneous railed values!)!')
    %             tIB3  = (LR.IB_PH < RCR.IB(1) | LR.IB_PH > RCR.IB(2)).*t_diag*0; % DIAGNOSTIC
    %             tIK3  = (IK < RCR.IK(1) | IK > RCR.IK(2)).*t_diag*0; % DIAGNOSTIC
    %             tIB4  = (LR.IB_PH < RCR.IB(1)*25 | LR.IB_PH > RCR.IB(2)*25).*t_diag*0; % DIAGNOSTIC
    %             tIK4  = (IK < RCR.IK(1)*25 | IK > RCR.IK(2)*25).*t_diag*0; % DIAGNOSTIC
    %         end
    %
    %         tALL  = max([LR.PH_IN_SITU_TOTAL_QC, tBSL, tSTP, tRC, ...
    %             tVRS, tIB3, tIK3, tIB4, tIK4],[],2); % get highest flag
    %         LR.PH_IN_SITU_TOTAL_QC(t_bio) = tALL(t_bio);
    %         LR.PH_IN_SITU_FREE_QC(t_bio)  = tALL(t_bio);
    %
    %         tALL  = max([LR.VRS_PH_QC, tBSL, tVRS, tIB3, tIK3, tIB4, tIK4],[],2); % get highest flag
    %         LR.VRS_PH_QC(t_bio) = tALL(t_bio);
    %
    %         tALL  = max([LR.IB_PH_QC, tBSL, tIB3, tIK3, tIB4, tIK4],[],2); % get highest flag
    %         LR.IB_PH_QC(t_bio)  = tALL(t_bio);
    %
    %
    %         %         % *****   ADJUSTED DATA   *****
    %         %         t_bio = LR.PH_IN_SITU_TOTAL_ADJUSTED ~= fv.bio;
    %         %         %tST     = LR.PSAL_ADJUSTED_QC == 4 | LR.TEMP_ADJUSTED_QC == 4; % Bad S or T will affect pH
    %         %         tST     = LR.PSAL_ADJUSTED_QC == 4 | LR.TEMP_ADJUSTED_QC == 4 |...
    %         %                   LR.PRES_ADJUSTED_QC == 4;
    %         %
    %         %         t_chk1 = t_bio & (LR.PH_IN_SITU_TOTAL_ADJUSTED < RC.PH(1)| ...
    %         %             LR.PH_IN_SITU_TOTAL_ADJUSTED > RC.PH(2) | tST);
    %         %         t_chk2 = t_bio & LR.IB_PH ~= fv.bio.*1e9 & LR.IB_PH ~= fv.bio & ...
    %         %             (LR.IB_PH < RC.IB(1) | LR.IB_PH > RC.IB(2) | ...
    %         %             IK < RC.IK(1) | IK > RC.IK(2));
    %         %         %t_chk3 = t_bio & (LRQF_S | LRQF_T);
    %         %         t_chk5 = t_bio & LR.IB_PH ~= fv.bio.*1e9 & LR.IB_PH ~= fv.bio ...
    %         %             & (LR.IB_PH < RCR.IB(1).*25 | LR.IB_PH > RCR.IB(2).*25 | ...
    %         %             IK < RCR.IK(1).*25 | IK > RCR.IK(2).*25); % DIAGNOSTIC
    %         %
    %         %         %LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(t_chk1 | t_chk2 |t_chk3) = 4;
    %         %         LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(t_chk2) = 3;  % set pH with out of range Ik/Ib to questionable, unless pH is also out of range (then set to bad)
    %         % 		LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(t_chk5) = 4;
    %         %         LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(t_chk1) = 4;
    %         %
    %         % 		% double check fill value QC flags - 9634 returns fv's from
    %         % 		qc_adjustmet
    %         % 		% due to bad gps date bug - jp 10/04/19
    %         % 		LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(~t_bio)  = fv.QC;
    %         %
    %         % 		LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(t_bio) =  ...
    %         %             LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(t_bio) * ~BSLflag + BSLflag*theflag;
    %
    %
    %         % 04/07/2020 JP - Now that the BSL flag can be 3 or 4 it has the
    %         % potential to overwrite a 4 with a 3 so don't let this happen
    %         t_bio = LR.PH_IN_SITU_TOTAL_ADJUSTED ~= fv.bio;
    %         tBSL  = ones(size(t_bio)) * ~BSLflag + BSLflag*theflag; % flag from BSL
    %         tSTP  = (LR.PSAL_ADJUSTED_QC == 4 | LR.TEMP_ADJUSTED_QC == 4 | ....
    %             LR.PRES_ADJUSTED_QC == 4)*4; % Bad STP will affect nitrate
    %         tRC   = (LR.PH_IN_SITU_TOTAL_ADJUSTED < RC.PH(1) | ...
    %             LR.PH_IN_SITU_TOTAL_ADJUSTED > RC.PH(2))*4; % Range check
    %         tVRS  = (LR.VRS_PH < RC.PHV(1) | LR.VRS_PH > RC.PHV(2))*4;
    %
    %         % double check fill value QC flags - 9634 returns fv's from qc_adjustmet
    %         % due to bad gps date bug - jp 10/04/19
    %         LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(~t_bio)  = fv.QC;
    %
    %         tALL  = max([LR.PH_IN_SITU_TOTAL_ADJUSTED_QC, tBSL, tSTP, tRC, ...
    %             tVRS, tIB3, tIK3, tIB4, tIK4],[],2); % get highest flag
    %         LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(t_bio) = tALL(t_bio);
    %
    %         % 09/30/20 TM - EXCLUSION BLOCK FOR NEWER EQPAC FLOATS EXHIBITING
    %         % THE ERRONEOUS RAILED DIAG VALUES
    %         if regexp(MBARI_ID_str, skip_ph_diag, 'once')
    %             if ~isempty(dura.data) && ph_filesize> 10000 && ... % usually > 10000 bytes for full profile...
    %                     any(LR.IB_PH < RC.IB(1) | LR.IB_PH > RC.IB(2))
    %                 disp([pH_file,' has out of range pH Ib values. pH ', ...
    %                     'QF''s for these samples will be set to bad or questionable'])
    %             end
    %             if ~isempty(dura.data) && ph_filesize> 10000 && ... % usually > 10000 bytes for full profile...
    %                     any(IK < RC.IK(1) | IK > RC.IK(2))
    %                 disp([pH_file,' has out of range pH Ik values. pH ', ...
    %                     'QF''s for these samples will be set to bad or questionable'])
    %             end
    %         end
    %
    %         % -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
    %         % FINALLY DO A SPIKE TEST ON VALUES, IF SPIKES ARE IDENTIFIED, SET QF TO 4
    %         % (BAD).  BGC_spiketest screens for nans and fill values.
    %         %
    %         % RUN TEST ON RAW PH_IN_SITU_TOTAL
    %         QCscreen_phT = LR.PH_IN_SITU_TOTAL_QC == 4; % screen for BAD data already assigned.
    %         [spike_inds, quality_flags] = BGC_spiketest(MBARI_ID_str,str2double(cast_num),[LR.PRES LR.PH_IN_SITU_TOTAL],'PH',dirs.cal,fv.bio,QCscreen_phT);
    %         if ~isempty(spike_inds)
    %             % Could add optional code to compare qf already assigned (ie in range
    %             % checks above), and qf assigned during spiketest.  Keep
    %             % whichever is greater.  (At this point, I'm pre-screening the vector using
    %             % current quality flags.  Screening occurs within Function.  Currently assigning 4 (bad) to spikes, although may
    %             % be modified in the future).  5/22/18. TM.
    %             %spikeQF = LR.PH_INSITU_FREE_QC(spike_inds) < quality_flags;
    %             %LR.PH_INSITU_FREE_QC(spike_inds(spikeQF)) = quality_flags(spikeQF);
    %             LR.PH_IN_SITU_TOTAL_QC(spike_inds) = quality_flags;
    %             disp(['LR.PH_IN_SITU_TOTAL QUALITY FLAGS ADJUSTED FOR IDENTIFIED SPIKES, PROFILE ',cast_num,' FLOAT ',MBARI_ID_str,'.'])
    %             %         else
    %             %             disp(['NO LR.PH_IN_SITU_TOTAL SPIKES IDENTIFIED FOR PROFILE ',cast_num,' FLOAT ',MBARI_ID_str,'.'])
    %         end
    %         %
    %         % RUN TEST ON QC PH_IN_SITU_TOTAL_ADJUSTED
    %         QCscreen_phTadj = LR.PH_IN_SITU_TOTAL_ADJUSTED_QC == 4; % screen for BAD data already assigned.
    %         [spike_inds, quality_flags] = BGC_spiketest(MBARI_ID_str,str2double(cast_num),[LR.PRES LR.PH_IN_SITU_TOTAL_ADJUSTED],'PH',dirs.cal,fv.bio,QCscreen_phTadj);
    %         if ~isempty(spike_inds)
    %             LR.PH_IN_SITU_TOTAL_ADJUSTED_QC(spike_inds) = quality_flags;
    %             disp(['LR.PH_IN_SITU_TOTAL_ADJUSTED QUALITY FLAGS ADJUSTED FOR IDENTIFIED SPIKES, PROFILE ',cast_num,' FLOAT ',MBARI_ID_str,'.'])
    %             %         else
    %             %             disp(['NO LR.PH_IN_SITU_TOTAL_ADJUSTED SPIKES IDENTIFIED FOR PROFILE ',cast_num,' FLOAT ',MBARI_ID_str,'.'])
    %         end
    %         % -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
    %
    %         clear phfree phtot QCD dura t_chk1 t_chk2 t_chk3 QCscreen_phF QCscreen_phT QCscreen_phTadj QCscreen_O QCscreenOadj spike_inds quality_flags
    %
    %         %if INFO.cast == 30,pause,end % TESTING
    %
end

%% ------------------- Calculate LIGHT ------------------------------------

% % Just copy output from OCR channels for now
LR.OCR1=d.lr_d(:,iOCR1);
LR.OCR2=d.lr_d(:,iOCR2);
LR.OCR3=d.lr_d(:,iOCR3);
LR.OCR4=d.lr_d(:,iOCR4);

%% ------------------- HR Data structure ----------------------------------

    % ****************************************************************
    % NOW DEAL WITH HIGH RESOLUTION (CP) DATA
    % ****************************************************************
    if isempty(d.hr_d) % CHECK FOR DATA
        disp(['No high res data in message file for ', msg_file]);% MC ...
            % strtrim(msg_list(msg_ct,:))])
        HR.PRES = [];
        HR.PRES_ADJUSTED = [];
        HR.PSAL = [];
        HR.PSAL_ADJUSTED = [];
        HR.TEMP = [];
        HR.TEMP_ADJUSTED = [];
        HR.NBIN_CTD = [];
    else % if data also header variable
        hr_d = d.hr_d; % High resolution data
        % GET VARIABLE INDICES - ONLY NEED TO DO THIS ONCE
        if hr_ind_chk == 0
            hr_ind_chk = 1;
            
            iPP   = find(strcmp('pres_dBar', d.hr_hdr) == 1); % CTD P
            iTT   = find(strcmp('temp_degC', d.hr_hdr) == 1); % CTD T
            iSS   = find(strcmp('sal_ppt', d.hr_hdr) == 1); % CTD S
            iBIN   = find(strcmp('nbin ctd', d.hr_hdr) == 1); % CTD N samples per bin
        end
        HR.PRES = hr_d(:,iPP);
        HR.PRES_ADJUSTED = HR.PRES;
        HR.PSAL = hr_d(:,iSS);
        HR.PSAL_ADJUSTED = HR.PSAL; % NO CORRECTIONS DONE FOR S OR T
        HR.TEMP = hr_d(:,iTT);
        HR.TEMP_ADJUSTED = HR.TEMP;
        HR.NBIN_CTD = hr_d(:,iBIN);
        
        HR.PSAL_QC = HR.PRES * 0 + fv.QC; % Predimmension QF's
        HR.PSAL_ADJUSTED_QC = HR.PSAL_QC;
        HR.TEMP_QC = HR.PSAL_QC;
        HR.TEMP_ADJUSTED_QC = HR.PSAL_QC;
        HR.PRES_QC = HR.PSAL_QC;
        HR.PRES_ADJUSTED_QC = HR.PSAL_QC;
        
        HRQF_P  = HR.PRES < RC.P(1) | HR.PRES > RC.P(2); % BAD HR S
        t_bio  = HR.PRES ~= fv.bio;
        HR.PRES_QC(HRQF_P)  = 4;  % BAD
        HR.PRES_QC(~HRQF_P & HR.PRES ~= fv.bio) = 1; % GOOD
        HR.PRES_ADJUSTED_QC(HRQF_P)  = 4;  % BAD
        HR.PRES_ADJUSTED_QC(~HRQF_P & HR.PRES_ADJUSTED ~= fv.bio) = 1; % GOOD
        
%         [BSLflag, theflag] = isbadsensor(BSL, MBARI_ID_str, INFO.cast, 'S');
%         HRQF_S  = HR.PSAL < RC.S(1) | HR.PSAL > RC.S(2); % BAD HR S
%         t_bio  = HR.PSAL ~= fv.bio;
%         HR.PSAL_QC(HRQF_S)  = 4;  % BAD
%         HR.PSAL_QC(~HRQF_S & HR.PSAL ~= fv.bio) = 1; % GOOD
%         HR.PSAL_QC(t_bio) = HR.PSAL_QC(t_bio) * ~BSLflag + BSLflag*theflag;
%         HR.PSAL_ADJUSTED_QC(HRQF_S)  = 4;  % BAD
%         HR.PSAL_ADJUSTED_QC(~HRQF_S & HR.PSAL_ADJUSTED ~= fv.bio) = 1; % GOOD
%         HR.PSAL_ADJUSTED_QC(t_bio) = HR.PSAL_ADJUSTED_QC(t_bio) * ~BSLflag + BSLflag*theflag;
%         
%         
%         [BSLflag, theflag] = isbadsensor(BSL, MBARI_ID_str, INFO.cast, 'T');
%         HRQF_T  = HR.TEMP < RC.T(1) | HR.TEMP > RC.T(2); % BAD HR T
%         t_bio   = HR.TEMP ~= fv.bio;
%         HR.TEMP_QC(HRQF_T)  = 4;  % BAD
%         HR.TEMP_QC(~HRQF_T & HR.TEMP ~= fv.bio) = 1; % GOOD
%         HR.TEMP_QC(t_bio) = HR.TEMP_QC(t_bio) * ~BSLflag + BSLflag*theflag;
%         HR.TEMP_ADJUSTED_QC(HRQF_T)  = 4;  % BAD
%         HR.TEMP_ADJUSTED_QC(~HRQF_T & HR.TEMP_ADJUSTED ~= fv.bio) = 1; % GOOD
%         HR.TEMP_ADJUSTED_QC(t_bio) = HR.TEMP_ADJUSTED_QC(t_bio) * ~BSLflag + BSLflag*theflag;
%        
%         % CHECK FOR ANY HIGH RESOLUTION RAW QUALITY FLAGS TO ADD
%         if ~isempty(FV_HRdata)
%             
%             iQF = find(strcmp(FV_HRdata.hdr,'QF') == 1); % QF indices
%             tFV         = FV_HRdata.data(:,2) == INFO.cast;
%             FV_cast     = FV_HRdata.data(tFV,:);   % get FV cast
%             FV_QF_sum   = sum(FV_cast(:,iQF),1); % sum of QF columns
%             if sum(FV_QF_sum) > 0 % ANY ODV QF FLAGS GREATER THAN ZERO?
%                 indQF   = iQF(FV_QF_sum > 0); % ODV QF colums w/ non zero flags
%                 for QF_ct = 1 : size(indQF,2)
%                     ind = strcmp(FV_HRdata.hdr{indQF(QF_ct)-1}, ...
%                         QCvars(:,1));
%                     % N col in floatviz but really no data
%                     if ~isfield(HR, [QCvars{ind,2},'_QC'])
%                         continue
%                     end
%                     
%                     % FLOAVIZ VAR MATCHES LIST, GET MATCHING QF's
%                     if sum(ind) > 0 && isfield(HR,QCvars{ind,2})
%                         ODV_QF  = [FV_cast(:,6),FV_cast(:,8), ...
%                             FV_cast(:,indQF(QF_ct))]; % P, T & QC
%                         if strcmp(QCvars{ind,2},'BBP700') == 1
%                             ODV_QF(ODV_QF(:,3) == 4,3) = 2;
%                             ODV_QF(ODV_QF(:,3) == 8,3) = 4;
%                         else
%                             ODV_QF(ODV_QF(:,3) == 4,3) = 3; % CONVERT TO ARGO VALUES
%                             ODV_QF(ODV_QF(:,3) == 8,3) = 4; % CONVERT TO ARGO VALUES
%                         end
%                         
%                         ARGO_QF = [HR.PRES, HR.TEMP, ...
%                             HR.([QCvars{ind,2},'_QC'])];
%                         ct = 0;
%                         for i = 1:size(ARGO_QF,1) % line by line just in case
%                             dP = abs(ODV_QF(:,1) - ARGO_QF(i,1)); % press
%                             dT = abs(ODV_QF(:,2) - ARGO_QF(i,2)); % temp
%                             min_dP = min(dP); % this should be 0 but def < 1m
%                             min_dT = min(dT); % this should be 0 but def < 1m
%                             if min_dP > 1 % poor match
%                                 disp(['Poor HR QF match for ARGO PRES = ', ...
%                                     num2str(ARGO_QF(i,1))])
%                                 continue
%                             end
%                             t1 = dP == min_dP & dT == min_dT; % there should only be one
%                             if ODV_QF(t1,3) > ARGO_QF(i,3) % ODV QF worse than ARGO
%                                 ARGO_QF(i,3) = ODV_QF(t1,3); % replace argo w/ ODV
%                                 ct =ct+1;
%                             end
%                         end
%                         HR.([QCvars{ind,2},'_QC']) = ARGO_QF(:,3); %update structure
%                         if ct > 0
%                             disp([num2str(ct), ' QC flags added from HR ODV file for ', ...
%                                 QCvars{ind,2}, '_QC'])
%                         end
%                     end
%                 end
%                 clear tFV FV_cast FV_QF_sum indQF QF_ct ODV_QF ARGO_QF ct i dP t1
%                 clear min_dP min_dT dT
%             end
%         end
%         
%         % CHECK FOR ANY HIGH RESOLUTION ADJUSTED QUALITY FLAGS TO ADD
%         if ~isempty(FV_HRQCdata)
%             
%             iQF = find(strcmp(FV_HRQCdata.hdr,'QF') == 1); % QF indices
%             tFV         = FV_HRQCdata.data(:,2) == INFO.cast;
%             FV_cast     = FV_HRQCdata.data(tFV,:);   % get FV cast
%             FV_QF_sum   = sum(FV_cast(:,iQF),1); % sum of QF columns
%             if sum(FV_QF_sum) > 0 % ANY ODV QF FLAGS GREATER THAN ZERO?
%                 indQF   = iQF(FV_QF_sum > 0); % ODV QF colums w/ non zero flags
%                 for QF_ct = 1 : size(indQF,2)
%                     ind = strcmp(FV_HRQCdata.hdr{indQF(QF_ct)-1}, ...
%                         QCvars(:,1));
%                     % N col in floatviz but really no data
%                     if ~isfield(HR, [QCvars{ind,2},'_ADJUSTED_QC'])
%                         continue
%                     end
%                     
%                     % FLOAVIZ VAR MATCHES LIST, GET MATCHING QF's
%                     if sum(ind) > 0 && isfield(HR,[QCvars{ind,2},'_ADJUSTED'])
%                         ODV_QF  = [FV_cast(:,6),FV_cast(:,8), ...
%                             FV_cast(:,indQF(QF_ct))]; % P, T & QC
%                         if strcmp(QCvars{ind,2},'BBP700') == 1
%                             ODV_QF(ODV_QF(:,3) == 4,3) = 2;
%                             ODV_QF(ODV_QF(:,3) == 8,3) = 4;
%                         else
%                             ODV_QF(ODV_QF(:,3) == 4,3) = 3; % CONVERT TO ARGO VALUES
%                             ODV_QF(ODV_QF(:,3) == 8,3) = 4; % CONVERT TO ARGO VALUES
%                         end
%                         
%                         ARGO_QF = [HR.PRES, HR.TEMP, ...
%                             HR.([QCvars{ind,2},'_ADJUSTED_QC'])];
%                         ct = 0;
%                         for i = 1:size(ARGO_QF,1) % line by line just in case
%                             dP = abs(ODV_QF(:,1) - ARGO_QF(i,1)); % press
%                             dT = abs(ODV_QF(:,2) - ARGO_QF(i,2)); % temp
%                             min_dP = min(dP); % this should be 0 but def < 1m
%                             min_dT = min(dT); % this should be 0 but def < 1m
%                             if min_dP > 1 % poor match
%                                 disp(['Poor HR QF match for ARGO PRES = ', ...
%                                     num2str(ARGO_QF(i,1))])
%                                 continue
%                             end
%                             t1 = dP == min_dP & dT == min_dT; % there should only be one
%                             if ODV_QF(t1,3) > ARGO_QF(i,3) % ODV QF worse than ARGO
%                                 ARGO_QF(i,3) = ODV_QF(t1,3); % replace argo w/ ODV
%                                 ct =ct+1;
%                             end
%                         end
%                         HR.([QCvars{ind,2},'_ADJUSTED_QC']) = ARGO_QF(:,3); %update structure
%                         if ct > 0
%                             disp([num2str(ct), ' QC flags added from HR ODV file for ', ...
%                                 QCvars{ind,2}, '_ADJUSTED_QC'])
%                         end
%                     end
%                 end
%                 clear tFV FV_cast FV_QF_sum indQF QF_ct ODV_QF ARGO_QF ct i dP t1
%                 clear min_dP min_dT dT
%             end
%         end
    end

%% RENAME data strcuture
PARK=LR;

%% ------------------- Save profile in .mat -------------------------------

% New root directory tracking, allows for local folder path
indsd=regexp(path_data,'Data', 'once'); % find index of D
dirs.root=path_data(1:indsd+4); %'/Volumes/CPF/CPF003/Data/
% dirs.root=path_data(1:25); %'/Volumes/CPF/CPF003/Data/

dirs.mat=[dirs.root,'matfiles/',dirs.temp(length(dirs.root)+1:end)]; % path to save .mat file

% create directory if it doesn't exist already
if ~exist(dirs.mat) & numel(dir(dirs.mat))<=2 % will create directory only if it doesn't exist already!
    mkdir(dirs.mat)
end

if isempty(data.lr_d)~=1 && isempty(data.hr_d)~=1 % MC, added this if to avoid saving a .mat file for last profile (which often has no data)

    display(['Saving .mat file for each profile ... ',dirs.mat])
     
    WMO='CPF003';
    save_str=[dirs.mat, WMO,'.', cast_num,'.mat'];
    
    if exist(save_str) & numel(dir(save_str))<=2 % if a mat file exist in the directory, just append PARK structure
        save(save_str,'PARK','-append');
    else
        save(save_str,'PARK','HR','INFO');
    end

end

%     % *********************************************************************
%     % SAVE THE PROFILE AS WMO_ID#.PROFILE.mat
%     % THE *.mat file will contain 3 structures:
%     %   LR      for low resolution data
%     %   HR      for high resolution data (constant profiling)
%     %   info	float info that may be of use for ARGO data stream, also
%     %           used to build ODV compatible text file from 8.mat files
%     % *********************************************************************
%     if exist('AIR_O2', 'var')
%         INFO.AIR_O2   = AIR_O2; % [p t s phase cor_phase uM uMsat pO2]
%     end
%     if exist('AIRCAL_O2', 'var')
%         % [sdn bladderPress P p t s phase cor_phase uM uMsat pO2]
%         INFO.AIRCAL_O2   = AIRCAL_O2;
%     end
%
%     %         WMO  = INFO.WMO_ID; % string
%     %         if isempty(WMO)
%     %             disp(['NO WMO# FOUND FOR FLOAT! CREATING TEMPORARY ', ...
%     %                   'DATA DIR FOR ', INFO.UW_ID])
%     %             WMO = ['NO_WMO_',INFO.UW_ID]; % CREATE TEMPORARY WMO NAME
%     %         end
%
%     save_str = [dirs.mat, WMO,'\', WMO,'.', cast_num,'.mat'];
%     if regexp(MBARI_ID_str, ocr_testflt, 'once')
%         save(save_str,'OCR','LR','HR','INFO');
%     else
%         save(save_str,'LR','HR','INFO');
%     end
%     %     if msg_ct == 1
%     %         copyfile(fp_cal, [dirs.mat, WMO,'\']); % copy over cal file
%     %     end