% calctd.m applies calibrations to  oxygen
clear all
o2_params
%for quality flags see
%http://woce.nodc.noaa.gov/woce_v3/wocedata_1/whp/exchange/exchange_format_desc.htm

global pl ipl %this is just for cleanox graphics
ipl=0;

% read uncalibrated data files and match local sigmas
%ctd files
% Modif glazeO - 16/01/2015
read_moorspice_nc

if p_onlyonecast %takes only the casts for which o2 samples were made
    if ~exist('giprof2cal')
        iprofmin=find(floor(gprof2get/100)==min(gistat2cal));
        iprofmax=find(floor(gprof2get/100)==max(gistat2cal));
        gprof2cal=gprof2get(iprofmin:iprofmax);
    end
    giprof=[];
    for iproff=1:length(gprof2cal)
        giprof=[giprof;find(actd.profileid==gprof2cal(iproff))];
    end
    giprof=giprof(:)
    profbeg=min(giprof);
    profend=max(giprof);
else %otherwise takes all casts
    gipbeg=find(actd.profileid==100*statbeg+castbeg);
    gipend=find(actd.profileid==100*statend+castend);
    giprof=gipbeg:gipend;
    if isempty(giprof),error('check castbeg and castend exist'),end
end

% load oxygen coefficients (and btl variable)
strf=sprintf('load Califiles/ox_cals_%s_is%d_ws%d_%s.mat',cname,iswitch,wswitch,csuf);
disp(strf)
eval(strf) %loads calib matrix 'final'

%allocate final/interpolated variables
if p_cleanox
    actd.oxcali=NaN*ones(size(actd.doxv));
    actd.oxcali_qc=9*ones(size(actd.doxv));
end
if p_cleansalt
    if p_cleanox | p_cleantemp , error('first use p_cleanox then p_cleansalt, not simultaneously'),end
    actd.dsalii=NaN*ones(size(actd.dsali));
    actd.dsalii_qc=9*ones(size(actd.dsali));
end
if p_cleantemp
    if p_cleanox | p_cleansalt , error('first use p_cleansalt then p_cleantemp, not simultaneously'),end
    actd.dtemp68i=NaN*ones(size(actd.dsali));
    actd.dtemp68i_qc=9*ones(size(actd.dsali));
end

%
% giprof is the profile index in the CTD file.
figure(1);figure(2);figure(3);figure(4);
clf(1);clf(2);clf(3);clf(4) %need to clear figs

for ipp=1:length(giprof) %17
    iprof=giprof(ipp);
    %ctd. variable contains current ctd profile while actd contails all
    %profiles (see read_pandora_nc)
    disp(        '*********************************')
    disp(sprintf('**** STATION %03d, CAST %02d ****',Stat(iprof),Cast(iprof)))
    disp(        '*********************************')
    
    % selects bottle data for graphics, later on
    %select bottles for that profile (for graphic purpose)
    gistatb=find((ismember(btl0.stat,Stat(iprof))==1)&(ismember(btl0.cast,Cast(iprof))==1));
    if isempty(gistatb)
        gistatb=find((ismember(btl0.stat,Stat(iprof))==1));
        %get bottles from other cast of that station
        p_othercast=1;
    else
        p_othercast=0;
    end
    %variable btl contains only bottles for current ctd profile
    btl=keep_data(btl0, gistatb);
    %find flagged bottles for graphics
    %good bottles
    gigoodbtl=find(btl.wnklo2_flag==2);
    %bad bottles
    gibadbtl=find(btl.wnklo2_flag==3);
    %measured salinit bottles:
    gigoodbtls=find(btl.bsali~=35); %35 exact means no data
    
    if ~(p_cleansalt | p_cleantemp)
        %
        % first selects the current ctd profiles
        %
        clear ctd
        ctd.stat=Stat(iprof); %current station
        ctd.cast=Cast(iprof); %current cast
        ctd.dayd=actd.dayd(iprof);
        ctd.latx=actd.latx(iprof);
        ctd.lonx=actd.lonx(iprof);
        ctd.iswitch=iswitch;
        ctd.wswitch=wswitch;
        ctd.csuf=csuf;
        for ivar=1:length(varstr)
            eval(sprintf('%s=a%s(:,iprof);',varstr{ivar},varstr{ivar}))
            eval(sprintf('%s(%s>=1.000e+36)=NaN;',varstr{ivar},varstr{ivar}))
        end
        %
        % osat and sdvdt
        %
        ctd.dosat=o2sol_gg92(ctd.dsali,ctd.dtemp68,ctd.dpres); %using potential temperature inside as for calibration
        
        % compute dv/dt
        if ~isfield(ctd,'ddvdt') % need variable "time" in this case
            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?
            disp(sprintf('********** SDVDT FILTER AT %d M FOR O2 DATA **********',p_filter))
        else
            ctd.dsdvdt=ctd.ddvdt;
        end
        if p_plot
            figure(12);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
        
        if p_debug
            figure(10);clf;plot(ctd.ddvdt,-ctd.dpres,'r',ctd.dsdvdt,-ctd.dpres,'b');legend('dvdt','filtered');
            figure(11);clf;plotxx(ctd.doxv,-ctd.dpres,ctd.dsdvdt,-ctd.dpres);title('ctd.doxv and ctd.sdvdt');ppause
        end
        
        icalmat=find(final.stat==ctd.stat); 
        if isempty(icalmat) %This profile had no sample. Where do we get the coefficients from?
            icalmat=find(final.stat==gistat_nosample{ctd.stat+1}); %+1 for station 0
        end
        if length(icalmat)>1, error('need only one station in calmatrix, per coefficient'), end
        %
        %can use equation with p_slope=0 because sbecoef (1) already contains linear slope dependency
        %
        p_slope=0;
        ctd.oxcal=oxyg_eq(p_equation, p_slope, final.sbecoefs(icalmat,1:5),sbecoefd, ctd);
        ctd.oxcal_qc=2*ones(size(ctd.oxcal)); %Acceptable measurement
        ctd.oxcal_qc(find(isnan(ctd.oxcal)))=4; %unless NaN: "Bad measurement
        
        % keep only downcast data (normally useless with netcdf output)
        imaxp=find(ctd.dpres==max(ctd.dpres));
        ctd=keep_data(ctd,1:imaxp);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Interpolation with cleanvar
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %
        % Oxygen
        %
        if p_cleanox
            %start with calibrated ox
            ctd.oxcali=ctd.oxcal;
            ctd.oxcali_qc=ctd.oxcal_qc;
            % plot ctd temp and sali to possibly detect problems
            figure(3);clf;plotxx(ctd.dtemp68,-ctd.dpres,ctd.dsali,-ctd.dpres);grid on;
            title(sprintf('temp/salt station %d cast %d',ctd.stat, ctd.cast))
            %ctd=cleanox(ctd,p_auto(ctd.stat),keep_data(btl, gigoodbtl),p_othercast);
            ttlstr=sprintf('oxyg station %d cast %d',ctd.stat, ctd.cast);
            [ctd.oxcali, ctd.oxcali_qc]=cleanvar(ctd.oxcali, ctd.oxcali_qc, ctd.dpres, ...
                btl.wnklo2um(gigoodbtl),btl.dpres(gigoodbtl),ttlstr,p_auto(iprof),p_othercast);
            %        else
            %            ctd.oxcali_qc=9*ones(size(ctd.oxcal));
        end
        
    else %~(p_cleansalt | p_cleantemp)
        %in this case p_cleanox was already run so we load the individual
        %ctd profile
        strfe=sprintf('load Caliprofiles/ox_cal_%s_STN%03d_CAST%02d.mat ',...
            cname,Stat(iprof),Cast(iprof));
        disp(strfe);
        eval(strfe)
        imaxp=length(ctd.dpres);
        title(sprintf('temp/salt station %d cast %d',ctd.stat, ctd.cast))
        %
        % Salinity (interpolation only)
        %
        if p_cleansalt
            figure(3);clf;plotxx(ctd.dtemp68,-ctd.dpres,ctd.dsali,-ctd.dpres);grid on;
            if exist('p_getfrompreviousdir') & p_getfrompreviousdir
                % temporary function to pass on data between two version: get
                % unspiked oxygen from that directory
                ctd=getfromprevious(ctd,p_getfrompreviousdir,cname);
                %if ctd.stat==13 & ctd.cast==2
                %    disp('quick oxygen fix for 1302')
                %    ctd.oxcali(1)=NaN;
                %    ctd.oxcali_qc(1)=4;
                %    ppause
                %end
            end
            ctd.dsali_qc=2*ones(size(ctd.dsali));
            ctd.dsali_qc(find(isnan(ctd.dsali)))=9;
            %temporary variable to modify and send to cleanvar routine
            ctd_dsali=ctd.dsali;
            ctd_dsali_qc=ctd.dsali_qc;
            % find where CTD O2 has been interpolated or removed
            ginan=find(isnan(ctd.oxcali));
            ginan(ginan>(length(ctd.dpres)-floor(p_filter/2)))=[]; %dismiss last NaNs due to filter
            giinterp=find(ctd.oxcali_qc==6);
            giflag=find(~isnan(ctd.oxcali)&(ctd.oxcali_qc==3|ctd.oxcali_qc==4|ctd.oxcali_qc==9));
            if ~isempty(ginan)
                gisegments=[1,(1+find(diff(ginan)~=1)');find(diff(ginan)~=1)',length(ginan)];
                disp('O2 Data removed between these depths')
                disp(ctd.dpres(ginan(gisegments)))
                yn=input('Remove also salinity ? (y/n)','s');
                switch yn
                    case 'y'
                        ctd_dsali(ginan)=NaN;
                        ctd_dsali_qc(ginan)=4;
                        disp('removed')
                        ppause
                end
            end
            if ~isempty(giinterp)
                gisegments=[1,(1+find(diff(giinterp)~=1)');find(diff(giinterp)~=1)',length(giinterp)];
                disp('O2 Data interpolated between these depths')
                disp(ctd.dpres(giinterp(gisegments)))
                yn=input('Interpolate also salinity ? (y/n)','s');
                switch yn
                    case 'y'
                        for isegment=1:size(gisegments,2)
                            kk=giinterp(gisegments(1,isegment):gisegments(2,isegment));
                            ll=[kk(1)-1,(max(kk)+1)];
                            ctd_dsali(kk)=interp1(ctd.dpres(ll),ctd.dsali(ll),ctd.dpres(kk));
                            ctd_dsali_qc(kk)=6;
                        end
                        disp('interpolated')
                        ppause
                end
            end
            if ~isempty(giflag)
                gisegments=[1,(1+find(diff(giflag)~=1)');find(diff(giflag)~=1)',length(giflag)];
                disp('O2 Data flagged between these depths (second matrix is flag)')
                disp(ctd.dpres(giflag(gisegments)))
                disp(ctd.oxcali_qc(giflag(gisegments)))
                yn=input('Flag also salinity ? (y/n)','s');
                switch yn
                    case 'y'
                        ctd_dsali_qc(giflag)=ctd.oxcali_qc(giflag);
                        disp('flagged')
                        ppause
                end
            end
            % plot ctd temp and sali to possibly detect problems
            ttlstr=sprintf('Salinity station %d cast %d',ctd.stat, ctd.cast);
            [ctd.dsalii, ctd.dsalii_qc]=cleanvar(ctd_dsali, ctd_dsali_qc, ctd.dpres, ...
                btl.bsali(gigoodbtls),btl.dpres(gigoodbtls),ttlstr,0,p_othercast);
        end %if p_cleansalt
        %
        % Temperature (interpolation only)
        %
        if p_cleantemp
            figure(3);clf;plotxx(ctd.dtemp68,-ctd.dpres,ctd.dsalii,-ctd.dpres);grid on;
            ctd.dtemp68_qc=2*ones(size(ctd.dtemp68));
            ctd.dtemp68_qc(find(isnan(ctd.dtemp68)))=9;
            ctd_dtemp68=ctd.dtemp68;
            ctd_dtemp68_qc=ctd.dtemp68_qc;
            % find where CTD salinity has been interpolated or removed
            ginan=find(isnan(ctd.dsalii));
            giinterp=find(ctd.dsalii_qc==6);
            giflag=find(~isnan(ctd.dsalii)&(ctd.dsalii_qc==3|ctd.dsalii_qc==4|ctd.dsalii_qc==9));
            if ~isempty(ginan)
                gisegments=[1,(1+find(diff(ginan)~=1)');find(diff(ginan)~=1)',length(ginan)];
                disp('Salinity Data removed between these depths')
                disp(ctd.dpres(ginan(gisegments)))
                yn=input('Remove also temperature ? (y/n)','s');
                switch yn
                    case 'y'
                        ctd_dtemp68(ginan)=NaN;
                        ctd_dtemp68_qc(ginan)=4;
                        disp('removed')
                        ppause
                end
            end
            if ~isempty(giinterp)
                gisegments=[1,(1+find(diff(giinterp)~=1)');find(diff(giinterp)~=1)',length(giinterp)];
                disp('Salinity Data interpolated between these depths')
                disp(ctd.dpres(giinterp(gisegments)))
                yn=input('Interpolate also temperature ? (y/n)','s');
                switch yn
                    case 'y'
                        for isegment=1:size(gisegments,2)
                            kk=giinterp(gisegments(1,isegment):gisegments(2,isegment));
                            ll=[kk(1)-1,(max(kk)+1)];
                            ctd_dtemp68(kk)=interp1(ctd.dpres(ll),ctd.dtemp68(ll),ctd.dpres(kk));
                            ctd_dtemp68_qc(kk)=6;
                        end
                        disp('interpolated')
                        ppause
                end
            end
            if ~isempty(giflag)
                gisegments=[1,(1+find(diff(giflag)~=1)');find(diff(giflag)~=1)',length(giflag)];
                disp('Salinity Data flagged between these depths (second matrix is flag)')
                disp(ctd.dpres(giflag(gisegments)))
                disp(ctd.dsalii_qc(giflag(gisegments)))
                yn=input('Flag also temperature ? (y/n)','s');
                switch yn
                    case 'y'
                        ctd_dtemp68_qc(giflag)=ctd.dsalii_qc(giflag);
                        disp('flagged')
                        ppause
                end
            end
            % plot ctd temp and sali to possibly detect problems
            figure(3);clf;plotxx(ctd_dtemp68,-ctd.dpres,ctd.dsalii,-ctd.dpres);grid on;
            title(sprintf('temp/salt station %d cast %d',ctd.stat, ctd.cast))
            ttlstr=sprintf('Temperature station %d cast %d',ctd.stat, ctd.cast);
            [ctd.dtemp68i, ctd.dtemp68i_qc]=cleanvar(ctd_dtemp68, ctd_dtemp68_qc, ctd.dpres, ...
                [],[],ttlstr,0,p_othercast);
            ctd.dtemp90i=ctd.dtemp68i/1.00024; %save only as ITS90
        end %if p_cleantemp
        
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Save as individual files per profile
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ctd.dtemp90=ctd.dtemp68/1.00024; %save only as ITS90
    
    if p_saveoutput
        eval(sprintf(...
            ['save Caliprofiles/ox_cal_%s_STN%03d_CAST%02d.mat ' ...
            'ctd iprof iswitch wswitch csuf'],...
            cname,ctd.stat,ctd.cast))
    end
    %Graphic to save
    %Use interpolated if available
    if isfield(ctd,'oxcali');
        ctdox=ctd.oxcali;
    else
        ctdox=ctd.oxcal;
    end
    if isfield(ctd,'dsalii')
        ctdsal=ctd.dsalii;
    else
        ctdsal=ctd.dsali;
    end
    if isfield(ctd,'dtemp68i')
        ctdtemp=ctd.dtemp68i;
    else
        ctdtemp=ctd.dtemp68;
    end
    
    figure(1);clf;
    %RED crosses for bottles correspond to current cast; blue to another
    %cast
    subplot(1,2,1);
    [ax1,ax2,hl1,hl2]=plotxx(ctdtemp,-ctd.dpres,ctdsal,-ctd.dpres);grid on;
    hold on;
    if p_othercast
        plo2s=plot(btl.bsali(gigoodbtls),-btl.dpres(gigoodbtls),'bx',btl.bsali(gigoodbtls),-btl.dpres(gigoodbtls),'bo');
    else
        plo2s=plot(btl.bsali(gigoodbtls),-btl.dpres(gigoodbtls),'rx',btl.bsali(gigoodbtls),-btl.dpres(gigoodbtls),'bo');
    end
    gibadbtls=intersect(gigoodbtls,gibadbtl); %this will report flaged oxygen but existing salinities
    plo2sb=plot(btl.bsali(gibadbtls),-btl.dpres(gibadbtls),'kd');
    for itx=1:length(gigoodbtls)
        if btl.dpres(gigoodbtls(itx)) <= max(ctd.dpres)
            text(btl.bsali(gigoodbtls(itx))+0.02,-btl.dpres(gigoodbtls(itx)),sprintf('%d',btl.btlid(gigoodbtls(itx))))
        end
    end
    set(ax2,'YTickLabel',[])
    set(ax2,'position',[0.1300    0.1100    0.3347    0.800],'ylim',[-max(ctd.dpres) 0]); %to avoid having PSU label outside
    set(ax1,'position',[0.1300    0.1100    0.3347    0.800],'ylim',[-max(ctd.dpres) 0])
    xlabel('PSU')
    axes(ax1);xlabel('deg C');axes(ax2)
    
    subplot(1,2,2);
    pl0=plot(ctdox,-ctd.dpres,'r'); hold on;grid on;
    title(sprintf('Oxygen station %d cast %d',ctd.stat, ctd.cast))
    set(gca,'position',[0.5703    0.1100    0.3347    0.800],'ylim',[-max(ctd.dpres) 0])
  % glazeo - for station 40 bifurcation :  set(gca,'xlim',[140 220],'XTick',[140:20:220])
    plo3=plot(btl.wnklo2um(gibadbtl),-btl.dpres(gibadbtl),'kd');
    if p_othercast
        plo2=plot(btl.wnklo2um(gigoodbtl),-btl.dpres(gigoodbtl),'bx',btl.wnklo2um(gigoodbtl),-btl.dpres(gigoodbtl),'bo');
    else
        plo2=plot(btl.wnklo2um(gigoodbtl),-btl.dpres(gigoodbtl),'rx',btl.wnklo2um(gigoodbtl),-btl.dpres(gigoodbtl),'bo');
        for itx=1:length(btl.wnklo2um)
            if btl.dpres(itx)<= max(ctd.dpres)
                text(btl.wnklo2um(itx)+2,-btl.dpres(itx),sprintf('%d',btl.btlid(itx)))
            end
        end
    end
    if ~isempty(plo2)
        if ~isempty(plo3)
            legend([pl0, plo2(2), plo3],'CTD', 'sig-matched samples', 'o2-flagged','Location','southeast')
        else
            legend([pl0, plo2(2)],'CTD', 'sig-matched samples', 'Location','southeast')
        end
    else
        legend([pl0],'CTD', 'Location','southeast')
    end
    xlabel('umol/kg')
    
    if p_final
                
        namefig=sprintf('ox_cal_%s_STN%03d_CAST%02d',cname,ctd.stat,ctd.cast);
        set(1,'paperunits','inch','paperorient','port')
        set(1,'paperposition',[0.3 0.5 5 6])
        
        nfopeps=sprintf('%s%s.eps',direps,namefig)
        exportfig(1,nfopeps,'preview','tiff','LineMode','scaled','color','rgb',...
            'FontMode','fixed','FontSize',13,'height',12);
        
        nfopjpg=sprintf('%s%s.png',dirjpg,namefig)
        exportfig(1,nfopjpg,'format','png','LineMode','scaled','color','rgb',...
            'FontMode','fixed','FontSize',13,'height',12);
        
        pause(2)
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Adds current profiles to variable
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    actd.oxcal(1:imaxp,iprof)=ctd.oxcal;
    if p_cleanox
        actd.oxcali(1:imaxp,iprof)=ctd.oxcali;
        actd.oxcali_qc(1:imaxp,iprof)=ctd.oxcali_qc;
    end
    if p_cleansalt
        actd.dsalii(1:imaxp,iprof)=ctd.dsalii;
        actd.dsalii_qc(1:imaxp,iprof)=ctd.dsalii_qc;
    end
    
end %ist0
if p_saveoutput
    if p_onlyonecast
        eval(sprintf(...
        ['save Califiles/all-o2-calibrated_%s_STN%03d_STN%03d.mat ' ...
        'actd iswitch wswitch csuf gisensorT gisensorS gisensorO'],...
        cname,profbeg,profend))
    else
        eval(sprintf(...
        ['save Califiles/all-o2-calibrated_%s_STN%03d_STN%03d.mat ' ...
        'actd iswitch wswitch csuf gisensorT gisensorS gisensorO'],...
        cname,gipbeg,gipend))
    end
end
