function CTDout=QC_CTD(CTD,savefile)
    %QC CTD:

    % (1) Time-series plot
    % This is used to manually NaN data, for example, if there are intermittent
    %  issues with the pump flow rate of the pH/O2 stream, or salinity
    %  spikes not caught by the filter (below).
    
    dn = datenum([1970 1 1]);
    
     figCTD = figure('Name','CTD','Units','inches','Position',[0 0 11 8.5],'CloseRequestFcn', @(h,e) fprintf(1, 'Not allowed, use "close CTD force"\n.'));
     ha = tight_subplot(3,1,0.075,0.1,0.1);
     g = get(ha,'Position');
     set(ha,'Box','on','NextPlot','add');
     
     axes(ha(1));
     plot(CTD(:,1),CTD(:,2),'.-');
     ylabel(['[' char(186) 'C]']);
     fUnNaN=@qcWG_UnNaNcb;
     fNaN=@qcWG_NaNcb;
     ha(1).UserData=zeros(length(CTD(:,2)),1); 
     title('Temperature');
     
     uicontrol('Units','normalized','Style','pushbutton','String','NaN',...
         'Position',[g{1}(1)+g{1}(3)-0.08 g{1}(2)+g{1}(4)+0.003 .08 .03],...
         'Callback',{fNaN,CTD(:,1),CTD(:,2),ha(1)});
     uicontrol('Units','normalized','Style','pushbutton','String','Un-NaN',...
         'Position',[g{1}(1)+g{1}(3)-0.16 g{1}(2)+g{1}(4)+0.003 .08 .03],...
         'Callback',{fUnNaN,CTD(:,1),CTD(:,2),ha(1)});
     
     axes(ha(2));
       plot(CTD(:,1),CTD(:,3),'.-');
       ylabel('');
       ha(2).UserData=zeros(length(CTD(:,2)),1);
       title('Salinity');
       
       uicontrol('Units','normalized','Style','pushbutton','String','NaN',...
         'Position',[g{2}(1)+g{2}(3)-0.08 g{2}(2)+g{2}(4)+0.003 .08 .03],...
         'Callback',{fNaN,CTD(:,1),CTD(:,3),ha(2)});
     uicontrol('Units','normalized','Style','pushbutton','String','Un-NaN',...
         'Position',[g{2}(1)+g{2}(3)-0.16 g{2}(2)+g{2}(4)+0.003 .08 .03],...
         'Callback',{fUnNaN,CTD(:,1),CTD(:,3),ha(2)});
     
     axes(ha(3));
       plot(CTD(:,1),CTD(:,4),'.-');
       ylabel(['[' char(181) 'mol kg^-^1]']);
       ha(3).UserData=zeros(length(CTD(:,2)),1);
       title('Oxygen');
       
       uicontrol('Units','normalized','Style','pushbutton','String','NaN',...
         'Position',[g{3}(1)+g{3}(3)-0.08 g{3}(2)+g{3}(4)+0.003 .08 .03],...
         'Callback',{fNaN,CTD(:,1),CTD(:,4),ha(3)});
     uicontrol('Units','normalized','Style','pushbutton','String','Un-NaN',...
         'Position',[g{3}(1)+g{3}(3)-0.16 g{3}(2)+g{3}(4)+0.003 .08 .03],...
         'Callback',{fUnNaN,CTD(:,1),CTD(:,4),ha(3)});
     
     set(ha,'XLim',[min(CTD(:,1)) max(CTD(:,1))] + [-1 1]*0.1);
     linkaxes(ha,'x');

    % (2) Sliding variance filter on salinity
    % These parameters tend to do a good job identifying spikes in salinity,
    %  which occur in pretty much every deployment. I have been implementing
    %  this as a rule for every WG mission. Additional spikes can be manually
    %   removed in the CTD time-series plot.
     win = 100; %window                                                                  
     nst = 5; %number of stdev's above median
     sm = slidingmed(CTD(:,3), win);
     ss = slidingstddev(CTD(:,3), win);
     fs = find(abs(CTD(:,3) - sm) > nst*ss);
     clear win nst sm ss;
    %check
     axes(ha(2));
     plot(CTD(fs,1),CTD(fs,3),'m.');
    %NaN bad data
     CTD(fs,3) = NaN;
     clear fs;                                                                   

    % (3) T vs S
    % This is used to further identify bad salinity data by looking for
    %  outliers in water-mass properties. Single-point outliers are often
    %  easily identified here. This only 'flags' the outliers; you still
    %  need to zoom-in and manually click the 'NaN' button in the time-series
    %  plot to NaN the data.
    
     figTS = figure('Name','TS','Units','inches','Position',[0 0 11 8.5]);
     plot(CTD(:,3),CTD(:,2),'.');
     hold on;
    %  plot(CTD(:,8),CTD(:,7),'g.'); %M1
     xlabel('S');
     ylabel(['T [' char(186) 'C]']);
     g = get(gca,'Position');
     uicontrol(figTS,'Units','normalized','Style','pushbutton','String','FLAG',...
       'Position',[g(1)+g(3)-0.08 g(2)+g(4)+0.003 .08 .03],...
       'Callback',...
       {@TSPlot,figTS,figCTD,CTD,ha});
     uiwait(figTS);
     clear f x y g gx gy;

    % (4) T vs O2
    % This is used to identify bad oxygen data by looking for outliers in
    %  water-mass properties. Note, according to Booth et al., O2 in MB is
    %  basically never < 62 umol/kg (2.0 mg/l), and is rarely below 143 umol/kg.
    %  This only 'flags' the outliers; you still need to zoom-in and manually
    %  click the 'NaN' button in the time-series plot to NaN the data.

     figTO = figure('Name','TO','Units','inches','Position',[0 0 11 8.5]);
     plot(CTD(:,2),CTD(:,4),'.');
     hold on;
     xlabel(['T [' char(186) 'C]']);
     ylabel(['O2 [' char(181) 'mol kg^-^1]']);
     g = get(gca,'Position');
     
     uicontrol('Units','normalized','Style','pushbutton','String','FLAG',...
       'Position',[g(1)+g(3)-0.08 g(2)+g(4)+0.003 .08 .03],...
       'Callback',...
       {@TOPlot,figTO,figCTD,CTD,ha});
     uiwait(figTO);
     clear f x y g gx gy;
     
     sprintf('Paused, Hit any key to save data and continue')
     pause;
     for n1=1:3
         bool=zeros(size(CTD,1),1);
         for n=1:size(CTD,1)
             samp=ha(n1).UserData(n,:);
             f=samp>0;
             if sum(f)==0
                 bool(n)=0;
                 continue
             end
             samp=samp(f);
             if samp(end)==2
                 bool(n)=0;
             else
                 bool(n)=1;
             end
         end
         filtered=CTD(:,n1+1);
         filtered(bool==1)=NaN;
         CTD(:,n1+1)=filtered;
     end
     
     close CTD force
    %NaN O2 data where S is NaN
     CTD(isnan(CTD(:,3)),4) = NaN;
    
    CTDout=CTD;
    
    %save QC'd .csv file
     time = (CTD(:,1) - dn) * 86400*1000;
     fid = fopen(savefile,'w');
     fprintf(fid, ['Time[GMT/Unix*1000],temperature[C],salinity,' ...
       'oxygen[umol/kg]']);
     for n=1:size(CTD,1)
       fprintf(fid, '\n%.0f,%.3f,%.3f,%.3f', [time(n) CTD(n,2:4)]);
     end
     fclose(fid);
     clear time fid n;

    clear ha le figCTD figTS figTO;
end