 %   CALIBRATION OF CTD-O2 DATA : script MATCH_SGN.M
%   - READ THE _CAL BOTTLE FILE
%   - READ DOWNCAST PROFILES
%   - FIND BOTTLE DENSITIES ON DOWNCAST PROFILE
%   - ASSIGN O2 SAMPLE DATA TO CTD ON SAME DENSITY FROM DOWNCAST
%   - GRAPHICS WITH EACH PROFILE THEN ENSEMBLE RESIDUALS
%   - SAVE A COMPLETE .CLO FILE AND MATLAB FOR STEP 3
%----------------------------------

% read uncalibrated bottle data file
%----------------------------------
clear;
p_match=1;
o2_params %get parameters

sig.lev=[0:1000:6000]; %define sigma reference pressure to use
nlev=length(sig.lev);

% Figure root names
namefig1='match_sig_pres';
namefig2='match_sig_oxyg';

fmergemat = strcat(dirmerge, cname, '_mergedbtl_cal.mat');
eval(['load ',fmergemat]);
%load variable btl.xxxx %with all samples and ctd_rosette info
% btl.upres: bottle depth ('U'pcast pressure)
% btl.dpres: new bottle depth on 'D'owncast
% btl.uoxv: ctd .btl (upcast)
% btl.doxv: ctd interpolated (downcast)

%
% compute bottle sigmas
%
for ilev=1:nlev
    btl.csg(:,ilev)=sw_pden(btl.usali,btl.utemp68,btl.upres,sig.lev(ilev))-1000;
end

% initialize output vectors (d= down ctd)

btl.dpres=NaN*zeros(size(btl.upres));
btl.dtemp68=NaN*zeros(size(btl.upres));
btl.dsali=NaN*zeros(size(btl.upres));
btl.doxv=NaN*zeros(size(btl.upres));
btl.ddvdt=NaN*zeros(size(btl.upres));
btl.dsdvdt=NaN*zeros(size(btl.upres));

% same but with "P" suffix is for pressure-matched only (for plot purpose only to see the
%difference
btl.dpresP=NaN*zeros(size(btl.upres));
btl.dtemp68P=NaN*zeros(size(btl.upres));
btl.dsaliP=NaN*zeros(size(btl.upres));
btl.doxvP=NaN*zeros(size(btl.upres));
btl.ddvdtP=NaN*zeros(size(btl.upres));
btl.dsdvdtP=NaN*zeros(size(btl.upres));

% read uncalibrated data files and match local sigmas
% ctd files
read_moorspice_nc

%
% start loop on profiles
%
for ii=1:length(gprof)
    iprof=find(actd.profileid==gprof(ii)); %current profile index in ctd vectors
    %ctd. variable in this loop is for each profile while actd contails all
    %profiles (see read_pandora_nc)
    clear ctd
    disp(        '*********************************')
    disp(sprintf('**** PROFILE %03d            ****',gprof(ii)))
    disp(        '*********************************')
    %get profile and remove useless data at the end
    gigood=1:max(find(~isnan(actd.dsali(:,iprof))));
    gibad=(max(gigood)+1):size(actd.dsali,1);
    if any(~isnan(actd.dsali(gibad,iprof))),error('data remaining below NaN'),end
    for ivar=1:length(varstr)
        eval(sprintf('%s=a%s(gigood,iprof);',varstr{ivar},varstr{ivar}))
        eval(sprintf('%s(%s>=1.000e+36)=NaN;',varstr{ivar},varstr{ivar}))
    end
    if any(diff(ctd.dpres)~=1)
        warning('ctd.dpres is not regularily monotonic: check!');ppause
    end
    
    % compute dv/dt
    if ~isfield(ctd,'ddvdt') % il faut time !
        ctd.ddvdt=conv(diff(ctd.doxv)./(diff(ctd.dtime)),[.5 .5]);
        ctd.ddvdt(1)=2*ctd.ddvdt(1);
        lmax=length(ctd.ddvdt);
        ctd.ddvdt(lmax)=2*ctd.ddvdt(lmax);
    end
    
    if p_filter
        % smooth ctd.ddvdt a bit with a filter
        % this should be a bit longer than the expected lag,
        % so a half-width of 5 dbar (9 points) for this 1-dbar
        % binned 1-m/s lowering rate cast.     
        ctd.dsdvdt=filtend(ctd.ddvdt,p_filter);     % shortened from 8 dbar (15 points), omit altogether?
    else
        ctd.dsdvdt=ctd.ddvdt;
    end
    if p_plot
        figure(1);clf;plot(ctd.ddvdt,-ctd.dpres,'r');hold on;
        plot(ctd.dsdvdt,-ctd.dpres,'b');set(gca,'xlim',[-0.08 0.02]);title('dv/dt');ppause
    end
    
    % find the deepest pressure so that we are just dealing with the downcast
    
    ipmax=find(ctd.dpres==max(ctd.dpres));
    ctd=keep_data(ctd,1:ipmax);
    
    % get rid of NaNs from delooping
    
    %ii=find(finite(ctd.dpres.*ctd.dtemp68.*ctd.dsali.*ctd.doxv.*ctd.dsdvdt)==0);
    gigood=find(~isnan(ctd.dpres.*ctd.dtemp68.*ctd.dsali.*ctd.doxv.*ctd.dsdvdt));
    ctd=keep_data(ctd,gigood);
        
    % compute CTD sigmas
    for ilev=1:nlev
        ctd.sg(:,ilev)=sw_pden(ctd.dsali,ctd.dtemp68,ctd.dpres,sig.lev(ilev))-1000;
        
        % add some random noise to make sure that no two sigmas are
        % the same and the sorted series is monotonic
        ctd.sg(:,ilev)=ctd.sg(:,ilev)+1e-8*randn(size(ctd.sg(:,ilev)));
    
        % now sort the sigmas
        [gisort,gjsort]=sort(ctd.sg(:,ilev)); 
        ctd.sgsorted(:,ilev)=ctd.sg(gjsort,ilev);
    
        % some plots to explore the size of the sorting and
        % its effect on oxygen volts...
        if p_plot
            figure(1);clf
            plot(ctd.dpres(gjsort)-ctd.dpres,-ctd.dpres);title('dpres sorted')
            pause
            plot(ctd.doxv,ctd.sg(:,ilev),ctd.doxv(gjsort),ctd.sg(gjsort,ilev),'r');
            title('effect dpres sorted on oxygen voltage')
            pause
        end
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % now find all bottle data for the correct station and cast number within
        % 500 dbar of the reference pressure for a given sigma
        gid500=find(btl.stat==Stat(iprof)&btl.cast==Cast(iprof)&...
            btl.upres>(sig.lev(ilev)-500) & btl.upres<=(sig.lev(ilev)+500));
        %disp(btl.upres(gid500));
        
        % now interpolate by sorted sigmas at
        % WILL WE MISS SOME TOP AND BOTTOM POINTS?  WE MIGHT, SO WE FORCE THE
        % ISSUE WITH fsgP and lsgP at either end....
        fsg=min([min(ctd.sgsorted(:,ilev)),min(btl.csg(gid500,ilev))])-1e-4;
        lsg=max([max(ctd.sgsorted(:,ilev)),max(btl.csg(gid500,ilev))])+1e-4;
        xsg=[fsg;ctd.sgsorted(:,ilev);lsg];
        btl.dpres(gid500)  =interp1(xsg,[ctd.dpres(1);  ctd.dpres;  ctd.dpres(ipmax)],  btl.csg(gid500,ilev));
        btl.dtemp68(gid500)=interp1(xsg,[ctd.dtemp68(1);ctd.dtemp68;ctd.dtemp68(ipmax)],btl.csg(gid500,ilev));
        btl.dsali(gid500)  =interp1(xsg,[ctd.dsali(1);  ctd.dsali;  ctd.dsali(ipmax)],  btl.csg(gid500,ilev));
        btl.doxv(gid500)   =interp1(xsg,[ctd.doxv(1);   ctd.doxv;   ctd.doxv(ipmax)],   btl.csg(gid500,ilev));
        btl.ddvdt(gid500)  =interp1(xsg,[ctd.ddvdt(1);  ctd.ddvdt;  ctd.ddvdt(ipmax)],  btl.csg(gid500,ilev));
        btl.dsdvdt(gid500) =interp1(xsg,[ctd.dsdvdt(1); ctd.dsdvdt; ctd.dsdvdt(ipmax)], btl.csg(gid500,ilev));
    end %ilev
    clear ilev %(just in case)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CANCEL SIGMA MATCH BELOW A CERTAIN DEPTH IF SPECIFIED
    % (p_gidmax_match)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    gid2p=find(btl.stat==Stat(iprof)&btl.cast==Cast(iprof)&...
        btl.upres>=p_gidmax_match(ii));
    if ~isempty(gid2p)
        disp('No match for these depths')
        disp(btl.upres(gid2p));ppause
        disp('No match for these depths')
        btl.dpres(gid2p)=btl.upres(gid2p);
        xpr=[ctd.dpres;ctd.dpres(ipmax)+1e-4];
        btl.dtemp68(gid2p)=interp1(xpr,[ctd.dtemp68;ctd.dtemp68(ipmax)],btl.upres(gid2p));
        btl.dsali(gid2p)  =interp1(xpr,[ctd.dsali;  ctd.dsali(ipmax)],  btl.upres(gid2p));
        btl.doxv(gid2p)   =interp1(xpr,[ctd.doxv;   ctd.doxv(ipmax)],   btl.upres(gid2p));
        btl.ddvdt(gid2p)  =interp1(xpr,[ctd.ddvdt;  ctd.ddvdt(ipmax)],  btl.upres(gid2p));
        btl.dsdvdt(gid2p) =interp1(xpr,[ctd.dsdvdt; ctd.dsdvdt(ipmax)], btl.upres(gid2p));
    end
    
    % now CHECK and if difference in pressure between upcast and sigma-matched
    % downcast values exceeds a limit, interpolate those values in pressure
    % WILL WE MISS SOME TOP AND BOTTOM POINTS?  WE MIGHT, SO WE FORCE THE
    % ISSUE WITH fpri and lpri at either end....
    gidmx50=find((btl.stat==Stat(iprof))&btl.cast==Cast(iprof)&...
        (abs(btl.upres-btl.dpres)>50)); %50 m limit
    if ~isempty(gidmx50)
        disp(sprintf(...
            'Station %d WARNING: NO DENSITY MATCH BECAUSE > 50M DEVIATION (MAY BE DUE TO HOMOGENEOUS REGION)',gprof(ii)));
        disp(sprintf(' AT %g DBARS',btl.upres(gidmx50)))
        for idd=1:length(gidmx50)
            gidep=find(abs(ctd.dpres-btl.upres(gidmx50(idd)))<200);
            figure(6);clf;subplot(1,2,1);plot(ctd.dsali(gidep),-ctd.dpres(gidep));
            title(sprintf('PROFILE %d , bottle depth %1.0f, matched %1.0f',gprof(ii),btl.upres(gidmx50(idd)),btl.dpres(gidmx50(idd))))
            subplot(1,2,2);plot(ctd.dtemp68(gidep),-ctd.dpres(gidep));ppause
        end
    end
    fpri=min([min(ctd.dpres),min(btl.upres(gidmx50))])-1;
    lpri=max([max(ctd.dpres),max(btl.upres(gidmx50))])+1;
    btl.dpres(gidmx50)=btl.upres(gidmx50);
    xpr50=[fpri;ctd.dpres;lpri];
    btl.dtemp68(gidmx50)=interp1(xpr50,[ctd.dtemp68(1);ctd.dtemp68;ctd.dtemp68(ipmax)],btl.upres(gidmx50));
    btl.dsali(gidmx50)  =interp1(xpr50,[ctd.dsali(1);  ctd.dsali;  ctd.dsali(ipmax)],  btl.upres(gidmx50));
    btl.doxv(gidmx50)   =interp1(xpr50,[ctd.doxv(1);   ctd.doxv;   ctd.doxv(ipmax)],   btl.upres(gidmx50));
    btl.ddvdt(gidmx50)  =interp1(xpr50,[ctd.ddvdt(1);  ctd.ddvdt;  ctd.ddvdt(ipmax)],  btl.upres(gidmx50));
    btl.dsdvdt(gidmx50) =interp1(xpr50,[ctd.dsdvdt(1); ctd.dsdvdt; ctd.dsdvdt(ipmax)], btl.upres(gidmx50));    
     
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % DO INTERPOLATION ON DOWNCAST IN PRESSURE TOO FOR PLOTS
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    gid2p=find(btl.stat==Stat(iprof)&btl.cast==Cast(iprof));
    btl.dpresP(gid2p)=btl.upres(gid2p);
    xpr=[ctd.dpres;ctd.dpres(ipmax)+1e-4];
    btl.dtemp68P(gid2p)=interp1(xpr,[ctd.dtemp68;ctd.dtemp68(ipmax)],btl.upres(gid2p));
    btl.dsaliP(gid2p)  =interp1(xpr,[ctd.dsali;  ctd.dsali(ipmax)],  btl.upres(gid2p));
    btl.doxvP(gid2p)   =interp1(xpr,[ctd.doxv;   ctd.doxv(ipmax)],   btl.upres(gid2p));
    btl.ddvdtP(gid2p)  =interp1(xpr,[ctd.ddvdt;  ctd.ddvdt(ipmax)],  btl.upres(gid2p));
    btl.dsdvdtP(gid2p) =interp1(xpr,[ctd.dsdvdt; ctd.dsdvdt(ipmax)], btl.upres(gid2p));
    
    %View downcast profile with matched data
    if p_plot 
        jj1=find(btl.stat==Stat(iprof)&btl.cast==Cast(iprof));
        figure(3);clf
        plot(ctd.doxv,-ctd.dpres,'r-',btl.doxv(jj1),-btl.dpres(jj1),'ro',btl.uoxv(jj1),-btl.upres(jj1),'go',btl.doxvP(jj1),-btl.upres(jj1),'bo')
        % plot(ctd.doxv,-ctd.dpres,'r-',btl.doxv(gid500),-btl.dpres(gid500),'b+',btl.uoxv(gid500),-btl.upres(gid500),'mo')
        title(sprintf('O2 voltage: correction of P, Profile %d', gprof(ii)));grid on
        xlabel('oxygen volts'); ylabel('Pressure [dbar]');
        legend('CTD downcast','sigma-matched','ctd up (.btl)','pressure only','location','southeast')
        
        figure(4);clf
        %ctd.doxum is the SEABIRD oxygen; btl.wnklo2umum is bottle data; btl.dpres is sigma
        %adjusted ctd.dpres on downcast; btl.upres is upcast bottle pressure
        plot(ctd.doxum,-ctd.dpres,'r-',btl.wnklo2um(jj1),-btl.dpres(jj1),'ro',btl.wnklo2um(jj1),-btl.upres(jj1),'go')
        % plot(ctd.doxv,-ctd.dpres,'r-',btl.doxv(gid500),-btl.dpres(gid500),'b+',btl.uoxv(gid500),-btl.upres(gid500),'mo')
        title(sprintf('Sample O2: correction of P, Profile %d', gprof(ii)));grid on
        xlabel('oxygen umol/kg'); ylabel('Pressure [dbar]');
        legend('CTD downcast','sigma matched','upcast pressure','location','southeast')
        pause
    end
%
% END LOOP ON PROFILES **********************************
%
end %ii (loop on profiles/ iprof)

% compute oxsat in umol/kg according to Benson & Krausse (1984)
% compute oxsol in umol/kg according to Garcia & Gordon (1992)
% check value is 274.610 at 10C and salinity 35 psu
%
% 1-on density-mached o2 (in principle that is what we want)
%
th = sw_ptmp(btl.dsali,btl.dtemp68,btl.dpres,0);
btl.dosat=o2sol_gg92(btl.dsali,btl.dtemp68,btl.dpres);

% compute CTD oxygen in umol/kg 
btl.doxum=oxyg_eq(p_equation, 0, sbecoef0,sbecoefd, btl);

if p_plot
    figure(5);clf
    plot(btl.uoxum-btl.doxum,-btl.dpres,'o');title('CTD oxygen UP - Oxygen recalculated DOWN');
    ppause
end
btl.ddoxum=btl.doxum-btl.wnklo2um; %oxygen difference with matched sigma
%
% 2-SAME CALCULATION BUT FOR PRESSURE-MATCHED OXYGEN (suffix P)
%     For some diagnostics
btl.dosatP=o2sol_gg92(btl.dsaliP,btl.dtemp68P,btl.dpresP);
params.stat=btl.stat;
params.dpres=btl.dpresP;
params.dtemp68=btl.dtemp68P;
params.dsali=btl.dsaliP;
params.doxv=btl.doxvP;
params.ddvdt=btl.ddvdtP;
params.dsdvdt=btl.dsdvdtP;
params.dosat=btl.dosatP;
btl.doxumP=oxyg_eq(p_equation, 0, sbecoef0,sbecoefd,params);
clear params

figure(1);clf
subplot(2,1,1)
plot(btl.dpres-btl.upres,-btl.upres,'ro');
title('Downcast - upcast sigma-matched pressures')
xlabel('Pressure difference [dbar]')
ylabel('Upcast pressures [dbar]')
grid on;
subplot(2,1,2)
plot(btl.stat,btl.dpres-btl.upres,'ro');
xlabel('Station number')
ylabel('Pressure difference [dbar]')
grid on;

figure(2); clf
subplot(2,1,1)
ii=find(btl.wnklo2um>0);
if p_plot_upcastbottles
    pl1=plot(btl.uoxum(ii)-btl.wnklo2um(ii),-btl.upres(ii),'go');
end
%    title('Upcast CTD - upcast bottle with pressure (green) or sigma-matched (red) oxygens')
title('Downcast CTD - upcast bottle with pressure or sigma-matched oxygens')
hold on
pl2=plot(btl.doxumP(ii)-btl.wnklo2um(ii),-btl.upres(ii),'bo'); %PRESSURE-MATCHED
pl3=plot(btl.doxum(ii)-btl.wnklo2um(ii),-btl.dpres(ii),'ro');
ylabel('Upcast pressures [dbar]')
xlabel('Oxygen difference [\mu mol  kg^{-1}]')
grid on;
if p_plot_upcastbottles
    legend([pl1,pl2,pl3],{'Upcast CTD','Pressure-match','Sigma-match'},'location','southeast')
else
    legend([pl2,pl3],{'Pressure-match','Sigma-match'},'location','southeast')
end 
%
subplot(2,1,2)
if p_plot_upcastbottles
    pl1=plot(btl.stat(ii),btl.uoxum(ii)-btl.wnklo2um(ii),'go');
end
hold on
plot(btl.stat(ii),btl.doxumP(ii)-btl.wnklo2um(ii),'bo'); 
plot(btl.stat(ii),btl.doxum(ii)-btl.wnklo2um(ii),'ro');
xlabel('Station number')
ylabel('Oxygen difference [\mu mol  kg^{-1}]')
grid on;

%Verify volume calibration for flasks (to be re-done after calibration too)
figure(5)
hold on
plot(btl.flask(ii),btl.doxumP(ii)-btl.wnklo2um(ii),'bo'); 
plot(btl.flask(ii),btl.doxum(ii)-btl.wnklo2um(ii),'ro');
xlabel('Flask id')
ylabel('Oxygen difference [\mu mol  kg^{-1}]')
grid on;
%
% SAVE OUTPUT
%
eval(['save ' dirmerge 'oxydata_' cname 'p btl'])

if p_final
    set(1,'paperunits','cent','paperorient','port')
    set(1,'position',[778   175   540   740])
    nfopjpg=sprintf('%s%s%d.png',dirjpg,namefig1,gisensorO)
    exportfig(1,nfopjpg,'format','png','LineMode','scaled','color','rgb',...
        'FontMode','fixed','FontSize',10,'height',25);
    
    nfopeps=sprintf('%s%s%d.eps',direps,namefig1,gisensorO)
    exportfig(1,nfopeps,'preview','tiff','LineMode','scaled','color','rgb',...
        'FontMode','fixed','FontSize',10,'height',25);
    
    set(2,'paperunits','cent','paperorient','port')
    set(2,'position',[778   175   540   740])
    nfopjpg=sprintf('%s%s%d.png',dirjpg,namefig2,gisensorO)
    exportfig(2,nfopjpg,'format','png','LineMode','scaled','color','rgb',...
        'FontMode','fixed','FontSize',10,'height',25);
    
    nfopeps=sprintf('%s%s%d.eps',direps,namefig2,gisensorO)
    exportfig(2,nfopeps,'preview','tiff','LineMode','scaled','color','rgb',...
        'FontMode','fixed','FontSize',10,'height',25);
    
end


