function [vari, vari_qc]=cleanvar(vari, vari_qc, varpres,btlvar,btlpres,ttlstr, p_auto_stat,p_othercast)
% visualize ctd variable (o2 or sali) that was calibrated; finds spikes and cut/ interpolates
% From PMEL G. Johnson's routine; adapted by A. Ganachaud (2014)
% INPUTS
%   p_auto_stat is for automated interpolation (0 or 1 for each station ; use if many small gaps)
%   btlvar /btlpres : bottles to overlay
%   p_othercast=1 if we do not plot the bottles from current cast but from
%   same station / different cast

%WITHIN THIS ROUTINE, ALL VARIABLES ARE FOR ONE SINGLE STATION

disp('*******************************')
disp('* running cleanvar')
disp('*******************************')
global pl ipl %this is just to manage plots from previous profiles


% KEEPS ONLY THE FOUR PREVIOUS PROFILES ON THE GRAPHIC
if ipl>=4
    delete(pl(1))
    pl(1)=[];ipl=ipl-1;
end

figure(4); %contains the profiles
pl0=plot(vari,-varpres,'r'); hold on;grid on;
title(ttlstr)

%overlay bottles
if p_othercast
    plo2=plot(btlvar,-btlpres,'bx',btlvar,-btlpres,'bo');
else
    plo2=plot(btlvar,-btlpres,'rx',btlvar,-btlpres,'bo');
end

if ipl>1
    if ~isempty(plo2)
        legend([pl0,pl(ipl-1), plo2(2)],'current station','previous stations', 'samples','Location','southeast')
    else
        legend([pl0,pl(ipl-1)],'current station','previous stations', 'Location','southeast')
    end
else
    if ~isempty(plo2)
        legend([pl0, plo2(2)],'current station', 'samples','Location','southeast')
    else
        legend([pl0],'current station', 'Location','southeast')
    end
end
%save axis
axc=axis; %save current axis

%automated interpolation (use only if many small gaps)
if p_auto_stat
    ipl=ipl+1;
    ginan=find(isnan(vari));
    if ~isempty(ginan)
        gilagend=[find(diff(ginan)~=1); length(ginan)];
        gilagdeb=[1;gilagend+1];
        if varpres(max(ginan))==max(varpres) %do not deal with ends
            gilagend(length(gilagend))=[];
            gilagdeb(length(gilagdeb))=[];
        end
        for ilag=1:length(gilagend)
            x1=ginan(gilagdeb(ilag))-1;
            x2=ginan(gilagend(ilag))+1;
            kk=x1+1:x2-1;
            ll=[x1,x2];
            vari(kk)=interp1(varpres(ll),vari(ll),varpres(kk));
            vari_qc(kk)=6;  %flag 6 as interpolated
        end %for ilag
        pl(ipl)=plot(vari, -varpres,'b');
        title([st ' interpolated profile (blue) OK ?'])
        ppause
        %axis([0 max(varpres) min(ox) max(ox)]);
        delete(pl0);
    else %if isempty(ginan)
        set(pl0,'color','b')
        title([st ' profile (blue) OK ?'])
    end
else %there was an interactive interpolation code, but now it is in the while loop below
end %p_auto

%second pass to remove points if needed

loop1=menu('Removal /interpolation ?','yes','no');

while loop1==1
    title('Zoom in on current area of interest, type ''return'' when done');
    set(pl0,'color','b')
    zoom on
    disp('enter instructions or type ''''return'''' on keyboard to continue')
    keyboard
    zoom off
    title('Click on end-points for removal/interpolation');
    legend off
    [x,y]=ginput(2);
    y=sort(-round(y));
    x1=find(varpres==y(1))
    x2=find(varpres==y(2))
    if isempty(x1)
        x1=1;
    elseif isempty(x2)
        x2=length(varpres);
    else %both are empty (double click)
        kk=[];
    end
    kk=x1+1:x2-1;
    if ~isempty(kk)
        ll=[x1,x2];
        yn=input('Interpolation ? y or n=removal or f=flag ','s');
        switch yn
            case 'y'
                vari(kk)=interp1(varpres(ll),vari(ll),varpres(kk));
                vari_qc(kk)=6;  %flag as interpolated
            case 'n'
                vari(kk)=NaN;
                vari_qc(kk)=4;  %flag as bad
            case 'f'
                fn=input('flag value (1 to 9)? ');
                vari_qc(kk)=fn;  %flag
        end
    end %if ~isempty(kk)
    pl1=pl0;
    set(pl1,'color','b')
    pl0=plot(vari,-varpres,'r');
    axis(axc)
    title('More Interpolating? (red = new profile)')
    disp(sprintf('treated between %d and %d db',varpres(min(kk)),varpres(max(kk))))
    loop1=menu('More Interpolating/removal?','yes','no');
    delete(pl1)
end
set(pl0,'color','g')
ipl=ipl+1;
pl(ipl)=pl0;
%delete bottles
delete(plo2)
end %end function

% interactive interpolation code (removed from the first if/else)
function dummy
if 0    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %interactive interpolation
    % remove artifacts interactively, interpolate, change WOCE flag to 6
    loop1=menu('Interpolate?','yes','no');
    while loop1==1
        title('Zoom in on current area of interest, type ''return'' when done');
        zoom on
        keyboard
        zoom off
        title('Click on end-points for interpolation');
        [x,y]=ginput(2);
        y=sort(-round(y));
        x1=find(varpres==y(1));
        x2=find(varpres==y(2));
        kk=x1+1:x2-1;
        if   (x(1)>min(varpres)) & (x(2)<max(varpres))
            ll=[x1,x2];
            vari(kk)=interp1(varpres(ll),vari(ll),varpres(kk));
            vari_qc(kk)=6;  %flag as interpolated
        else
            if (x(2)>max(varpres)) %last depths
                kk=x1+1:length(varpres);
            else %first depths
                kk=1:x2-1;
            end
            vari(kk)=NaN;
            vari_qc(kk)=4;  %flag as "bad"
        end
        pl1=pl0;
        set(pl1,'color','b')
        pl0=plot(vari,-varpres,'r');
        axis(axc)
        title('More Interpolating? (blue=before interpolation)')
        loop1=menu('More Interpolating?','yes','no');
        delete(pl1)
    end
end
end