%merge_ctdbtl_chimie
%FOR GIVEN STATION GROUP:
%  READ ROSETTE FILE (.BTL)
%  READ OXYGEN EXCEL SHEET
%  GROUP THEM IN A .cal (formerly .CLO) FILE
%0-READ ROSETTE DATA
%1-READ BOTTLE DATA
%2-FIND BOTTLE CLOSE TO ROSETTE1 BOTTLE and ASSOCIATE O2 (WARNING IF 2 CLOSE BOTTLESOR NO MATCH
%3- WRITE _cal.mat and .CAL bottle file with all stations for G. Johnsons calib routine
%(match_sgn_xx)
%oxygen is converted from ml/l to umol/kg here using TS from ctd rosette
%data; but also keep ml/l
%A. Ganachaud, Pandora Cruise, 2012
%revision Dec 2014: use profiles instead of stations;

clear all
p_merge=1;
o2_params;  %get parameters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Output .cal text file for compatibility with PMEL software
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fidmerge    = fopen(strcat(dirmerge, cname, '_mergedbtl.cal'),'w');
fmergemat = strcat(dirmerge, cname, '_mergedbtl_cal.mat');

% READ ROSETTE DATA
%-------------------
% Modif filename, glazeO - 15/01/2014
% fros=sprintf('%s%s_btl.nc',dirctd,cname);
fros=sprintf('%s%s%s-ALL_btl.nc',dirctd,'OS_',upper(cname));
if ~exist(fros,'file')
    error(sprintf('file does not exist %s',fros))
end
ncdisp(fros)
profileid=ncread(fros,'PROFILE');
Stat=floor(profileid/100);
Cast=profileid-100*Stat;

%check profile order; re-order if necessary so that giprof contains indexes
%corresponding to profiles given by gprof
for ii=1:length(gprof)
    giprof(ii)=find(profileid==gprof(ii));
end
giprof=giprof(:);
% giprof contains profile indexes in the CTD file.
%
disp(['  Converting temperature from ITS-90 to IPTS-68'])
allbtls=ncread(fros,'BOTL');
cbtlpresall=ncread(fros,'PRES');
% LOAD OXYGEN FOR A SINGLE CHANNEL FOR CALIBRATION
eval(sprintf(['cbtloxVall=ncread(fros,' '''' 'DOV%d' '''' ');'],gisensorO));
eval(sprintf(['cbtloxVdtall=ncread(fros,' '''' 'DVT%d' '''' ');'],gisensorO));
eval(sprintf(['cbtluoxmlall=ncread(fros,' '''' 'DO%d1' '''' ');'],gisensorO));%in ml/l
eval(sprintf(['cbtluoxumall=ncread(fros,' '''' 'DO%d2' '''' ');'],gisensorO));%in umol/kg
% LOAD T AND S FOR EACH SPECIFIED CHANNEL
%eval(sprintf(['btltemp68all=1.00024*ncread(fros,' '''' 'TE0%d' '''' ');'],gisensorT));
%eval(sprintf(['btlusaliall=ncread(fros,' '''' 'PSA%d' '''' ');'],gisensorS));
eval(sprintf(['btltemp68all1=1.00024*ncread(fros,' '''' 'TE0%d' '''' ');'],1));
eval(sprintf(['btlusaliall1=ncread(fros,' '''' 'PSA%d' '''' ');'],1));
eval(sprintf(['btltemp68all2=1.00024*ncread(fros,' '''' 'TE0%d' '''' ');'],2));
eval(sprintf(['btlusaliall2=ncread(fros,' '''' 'PSA%d' '''' ');'],2));
btltemp68all=nan(size(btltemp68all1));
btlusaliall=nan(size(btlusaliall1));


% FIRST READS XLSX OXYGEN FILE
%-----------------------------
if p_read_xls_all %read xls for all station once
    nfsamplexlsx=sprintf('%s%s','Oxygen-',lower(cruiseName), excelExtension);
    disp(sprintf('Reading chemistry stat in file %s',nfsamplexlsx))
    [Sbtl]=read_xlsx_o2_all([diro2 nfsamplexlsx],nbtl,Stat,Cast);
else %o2 chemistry has multiple excel sheets
    clear sbtl_xls
    for ii=1:length(giprof)
        iprof=giprof(ii);
        
        nfsamplexlsx=sprintf('%s_O2_Stations_%03d_a_%03d.%s',cruiseName,10*floor((Stat(iprof)-1)/10)+1,...
            10*(1+floor((Stat(iprof)-1)/10)), excelExtension);
        disp('')
        disp('***********************************************')
        
        disp(sprintf('Reading chemistry stat %d cast %d in file %s',Stat(iprof),Cast(iprof),nfsamplexlsx))
        if ~exist('sbtl_xls')
            [sbtl_xls]=read_xlsx_o2_all([diro2 nfsamplexlsx],nbtl);
        end
        iprof_xls_found=[];
        for iprof_xls=1:length(sbtl_xls)
            if sbtl_xls(iprof_xls).profileid==profileid(iprof)
                iprof_xls_found=iprof_xls;
                break
            end
        end
        if isempty(iprof_xls_found) %it was not there: need to load xlsx
            disp(['PROFILE ', num2str(profileid(iprof)), ' NOT FOUND: LOADING ', nfsamplexlsx])
            [sbtl_xls]=read_xlsx_o2_all([diro2 nfsamplexlsx],nbtl);
        end
        iprof_xls_found=[];
        for iprof_xls=1:length(sbtl_xls)
            if sbtl_xls(iprof_xls).profileid==profileid(iprof)
                iprof_xls_found=iprof_xls;
                break
            end
        end
        if isempty(iprof_xls_found)
            error('Profile not in xls sheet!')
        end
        
        if any(diff(sbtl_xls(iprof_xls).btlid)<1)
            error('check bottle order')
        end
        Sbtl(iprof).stat=sbtl_xls(iprof_xls_found).stat;
        Sbtl(iprof).cast=sbtl_xls(iprof_xls_found).cast;
        Sbtl(iprof).btlid=sbtl_xls(iprof_xls_found).btlid;
        Sbtl(iprof).pres=sbtl_xls(iprof_xls_found).pres;
        Sbtl(iprof).sali=sbtl_xls(iprof_xls_found).sali;
        Sbtl(iprof).temp=sbtl_xls(iprof_xls_found).temp;
        Sbtl(iprof).flask=sbtl_xls(iprof_xls_found).flask;
        Sbtl(iprof).o2ml=sbtl_xls(iprof_xls_found).o2ml;
    end
end
% Initialize output variable btl (dimension=as many total samples)
%-------------------------------
btl.stat=[];%station
btl.cast=[];%cast
btl.btlid=[];%bottle id
btl.upres=[];%u = up ctd pressure
btl.utemp68=[];% up ctd temperature
btl.bsali=[];%salinity from samples
btl.usali=[];% up ctd salinity
btl.flask=[];%flask (flacon) id
btl.uoxv=[];%up ctd o2 voltage
btl.uoxum=[];%up ctd o2 umol/kg
btl.uoxml=[];%up ctd o2 ml/l
btl.wnklo2um=[];%oxygen from samples (umol/kg)
btl.wnklo2ml=[];%oxygen from samples (ml/l)

for ii=1:length(giprof)
    iprof=giprof(ii);
    clear gibtl dsapr btltemp68 roxyg roxV rdoxVdt ro2id ro2 sbtl
    disp(sprintf('************************************************'))
    disp(sprintf('PROFILE %d',profileid(iprof)))
    disp(sprintf('************************************************'))
    %SWITCH TO THE CHOOSEN CHANNEL
    isensor=find(gprof==profileid(iprof)); % because gisensorTu/Su are on gprof index
    switch gisensorTu(isensor)
        case 1
            btltemp68all(:,iprof)=btltemp68all1(:,iprof);
        case 2
            btltemp68all(:,iprof)=btltemp68all2(:,iprof);
            disp(sprintf('Use Secondary temperature; type keyboard'));pause
    end
    switch gisensorSu(isensor)
        case 1
            btlusaliall(:,iprof)=btlusaliall1(:,iprof);
        case 2
            btlusaliall(:,iprof)=btlusaliall2(:,iprof);
            disp(sprintf('Use Secondary Salinity; type keyboard'));pause
    end
    gibtli=1:max(find(~isnan(allbtls(:,iprof))));
    cgibtl=allbtls(gibtli,iprof); %index of all valid bottles in -ctd
    cbtlpres=cbtlpresall(gibtli,iprof);
    btltemp68=btltemp68all(gibtli,iprof);
    btlusali=btlusaliall(gibtli,iprof);
    if any(abs(btlusali-35))>2
        error('check salinity first!')
    end
    btluoxum=cbtluoxumall(gibtli,iprof); %umol/kg
    btluoxml=cbtluoxmlall(gibtli,iprof); %ml/l
    btluoxv=cbtloxVall(gibtli,iprof);
    btldoxVdt=cbtloxVdtall(gibtli,iprof);
    btltpot=sw_ptmp(btlusali,btltemp68,cbtlpres,0);
    % compute oxsol in umol/kg according to Garcia & Gordon (1992)
    % check value is 274.610 at 10C and salinity 35 psu
    btlosat=o2sol_gg92(btlusali, btltemp68,cbtlpres);
    disp(sprintf('************************************************'))
    disp(sprintf('FROM CTD .btl FILE'))
    disp(sprintf('************************************************'))
    disp('btlid     pres             te68   th68    salt    oxst    oxV     doVdt       ox   ')
    disp(sprintf('%2d\t %7.1f \t%4.1f  \t%4.1f  \t%4.1f  \t%4.1f  \t%4.1f  \t%4.5f  \t%4.1f\n',...
        [cgibtl';cbtlpres';btltemp68';btltpot';btlusali';btlosat';btluoxv';btldoxVdt';btluoxum']))
    cbtldens=sw_dens(btlusali,btltemp68,cbtlpres)-1000;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % READ BOTTLE (CHEMISTRY) DATA
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % variables with prefix s (= sample)
    % recovers the chemistry data by profile
    sbtl.stat=Sbtl(iprof).stat(:);
    sbtl.cast=Sbtl(iprof).cast(:);
    sbtl.btlid=Sbtl(iprof).btlid(:);
    sbtl.pres=Sbtl(iprof).pres(:);
    sbtl.sali=Sbtl(iprof).sali(:);
    sbtl.temp=Sbtl(iprof).temp(:);
    sbtl.flask=Sbtl(iprof).flask(:);
    sbtl.o2ml=Sbtl(iprof).o2ml(:);
    if any(diff(sbtl.btlid)<1)
        error('check bottle order')
    end
    gigood=find(~isnan(sbtl.o2ml));
    if ~p_keep_all_btl %keep only bottles with Winkler data
         sbtl.stat=sbtl.stat(gigood);
        sbtl.cast=sbtl.cast(gigood);
        sbtl.pres=sbtl.pres(gigood);
        sbtl.o2ml=sbtl.o2ml(gigood);
        sbtl.btlid=sbtl.btlid(gigood);
        sbtl.sali=sbtl.sali(gigood);
        sbtl.temp=sbtl.temp(gigood);
        sbtl.flask=sbtl.flask(gigood);
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % FIND BOTTLE ID CORRESPONDING TO ROSETTE DEPTH
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %looks for the bottle ID and verifies that depth agrees
    
    pdevmax=0.08; %deviation in percent between rosette depth and excel bottle
    nbtlstat=length(cbtlpres);
    btlwklo2um=NaN*ones(nbtlstat,1);
    btlwklo2ml=NaN*ones(nbtlstat,1);
    btlid=NaN*ones(nbtlstat,1);
    btlflask=NaN*ones(nbtlstat,1);
    btlbsali=NaN*ones(nbtlstat,1);
    ibtlros=1; %ibtlros is index in rosette (.btl) file
%    for ibs=1:length(sbtl.pres) %ibs is index is sample variables
    for ibs=1:nbtlstat %ibs is index is sample variables
        while cgibtl(ibtlros)<sbtl.btlid(ibs) %FIND CORRESPONDING BOTTLE ID this way to manage 2 casts
            ibtlros=ibtlros+1;
            if ibtlros>nbtlstat
                error('chemistry bottle not found in .ros file')
            end
        end
        if cgibtl(ibtlros)~=sbtl.btlid(ibs)
            error(sprintf('bottle %d not found in rosette file',sbtl.btlid(ibs)))
        end
        if abs(sbtl.pres(ibs)-cbtlpres(ibtlros)) >= max(5,(min(cbtlpres(ibtlros),500)*pdevmax))
            error(sprintf('large pressure deviation between chemistry %4.0f db and rosette %4.0f db files',...
                sbtl.pres(ibs),cbtlpres(ibtlros)))
        end
        if 1 %comparison in umol/kg
            disp(sprintf('Bottle %2d found (ros=%6.1f db btl=%6.1f db ; o2ros-o2btl = %4.1f',...
                cgibtl(ibtlros),cbtlpres(ibtlros),sbtl.pres(ibs),btluoxum(ibtlros)-ox_units(sbtl.o2ml(ibs),btlusali(ibtlros),20*ones(length(ibtlros),1))))
        else
            disp(sprintf('Bottle %2d found (ros=%6.1f db btl=%6.1f db ; o2ros-o2btl = %4.1f',...
                cgibtl(ibtlros),cbtlpres(ibtlros),sbtl.pres(ibs),btluoxml(ibtlros)-sbtl.o2ml(ibs)))
        end
        
        %
        % Now add the sample data to rosette variables
        %
        %btlwklo2um(ibtlros)=ox_units(sbtl.o2ml(ibs),btlusali(ibtlros),btltpot(ibtlros)); %with conversion to umol/kg
        %Change AG (F. Baurand) Aug 2015: conversion is for lab conditions
        %at 20?C. Difference reaches 0.4umol/kg if Tpot at 1.8?C
        btlwklo2um(ibtlros)=ox_units(sbtl.o2ml(ibs),btlusali(ibtlros),20*ones(length(ibtlros),1)); %with conversion to umol/kg
        btlwklo2ml(ibtlros)=sbtl.o2ml(ibs);
        btlid(ibtlros)=sbtl.btlid(ibs);
        btlflask(ibtlros)=sbtl.flask(ibs);
        btlbsali(ibtlros)=sbtl.sali(ibs);
    end
    
    % REMOVE NaN BOTTLES
    if ~p_keep_all_btl %keep only bottles with Winkler data
        gibtlnan=find(isnan(btlwklo2um));
        cbtlpres(gibtlnan)=[];
        btlwklo2um(gibtlnan)=[];
        btlwklo2ml(gibtlnan)=[];
        btlid(gibtlnan)=[];
        btlflask(gibtlnan)=[];
        btlbsali(gibtlnan)=[];
        btltemp68(gibtlnan)=[];
        btltpot(gibtlnan)=[];
        btlusali(gibtlnan)=[];
        cbtldens(gibtlnan)=[];
        btlosat(gibtlnan)=[];
        btluoxv(gibtlnan)=[];
        btldoxVdt(gibtlnan)=[];
        btluoxum(gibtlnan)=[];
        btluoxml(gibtlnan)=[];
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 4-Write .cal using structure variable btl.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    nbtlstat=length(btlid);
    btl.stat=[btl.stat;Stat(iprof)*ones(nbtlstat,1)];
    btl.cast=[btl.cast;sbtl.cast(1:nbtlstat)];
    btl.btlid=[btl.btlid;btlid(:)];
    btl.upres=[btl.upres;cbtlpres(:)];
    btl.utemp68=[btl.utemp68;btltemp68(:)];
    btl.bsali=[btl.bsali;sbtl.sali(1:nbtlstat)];
    btl.usali=[btl.usali;btlusali(:)];
    btl.flask=[btl.flask;btlflask(:)];
    btl.uoxv=[btl.uoxv;btluoxv(:)];
    btl.uoxum=[btl.uoxum;btluoxum(:)];
    btl.uoxml=[btl.uoxml;btluoxml(:)];
    btl.wnklo2um=[btl.wnklo2um;btlwklo2um(:)];
    btl.wnklo2ml=[btl.wnklo2ml;btlwklo2ml(:)];
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 4bis-Write .cal file for compatibility with PMEL software
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    uu=[Stat(iprof)*ones(nbtlstat,1),sbtl.cast(1:nbtlstat), btlid(:), cbtlpres(:), btltemp68(:), ...
        sbtl.sali(1:nbtlstat), btlusali(:), btlflask(:),...
        btluoxv(:),btluoxum(:),btlwklo2um(:)];
    fprintf(fidmerge,'%3d %3d %3d %10.3f %8.4f %11.6f %9.4f %5d %9.5f %9.3f %9.3f\n',uu');
end %for iprof
fclose(fidmerge);
%add profile ID to btl (usefull thereafter)
btl.profileid=100*btl.stat+btl.cast;
eval(sprintf('save %s btl',fmergemat))
%next step is match_sgn_xxx.m

