function CO2out=QC_pCO2(CO2,CO2equil,CO2zero,CO2air,CO2stnd,CO2hum,CTD,savefile)
dn = datenum([1970 1 1]);

%QC pCO2 and pH:

load('CO2SYS_2100.mat');
load('CO2SYS_2400.mat');

% (1) Time-series plot
% This is used to manually NaN data, for example, when the diagnostics
%  data indicate a leak or blockage in the gas flow path, or the presence
%  of water vapor in the OASIS can.
 figCO2 = figure('Name','CO2','Units','inches','Position',[0 0 11 7.5],...
     'CloseRequestFcn', @(h,e) fprintf(1, 'Not allowed, use "close CO2 force"\n.'));
 ha=tight_subplot(7,1,0.05,0.1,0.1);
 g = get(ha,'Position');
 set(ha,'Box','on','NextPlot','add');
 axes(ha(1)); %press
   plot(CO2equil(:,1),CO2equil(:,5),'.-','Color',[.8 0 0]);
   plot(CO2zero(:,1),CO2zero(:,5),'.-','Color',[0 0 0]);
   plot(CO2air(:,1),CO2air(:,5),'.-','Color',[0 0 .8]);
   f = find(~isnan(CO2stnd(:,7)));
   plot(CO2stnd(f,1),CO2stnd(f,5),'.-','Color',[0 .8 0]);
   ylabel('[kPa]');
   title('Pressure');
 axes(ha(2)); %temp
   plot(CO2equil(:,1),CO2equil(:,3),'.-','Color',[.8 0 0]);
   plot(CO2zero(:,1),CO2zero(:,3),'.-','Color',[0 0 0]);
   plot(CO2air(:,1),CO2air(:,3),'.-','Color',[0 0 .8]);
   plot(CO2stnd(f,1),CO2stnd(f,3),'.-','Color',[0 .8 0]);
   ylabel('[C]');
   title('Temperature');
 axes(ha(3)); %humidity
   plot(CO2hum(:,1),CO2hum(:,2),'.-','Color',[.8 0 0]);
   plot(CO2hum(:,1),CO2hum(:,3),'.-','Color',[0 0 .8]);
   plot(CO2hum(f,1),CO2hum(f,4),'.-','Color',[0 .8 0]);
   ylabel('[V]');
   title('Humidity');
 axes(ha(4)); %xCO2
   plot(CO2equil(:,1),CO2equil(:,7),'.-','Color',[.8 0 0]);
   plot(CO2zero(:,1),CO2zero(:,7),'.-','Color',[0 0 0]);
   plot(CO2air(:,1),CO2air(:,7),'.-','Color',[0 0 .8]);
   plot(CO2stnd(f,1),CO2stnd(f,7),'.-','Color',[0 .8 0]);
   ylabel('[ppm]');
   title('xCO2');
   
 axes(ha(5)); %pCO2w
   f = find(~isnan(CO2(:,3)));
   plot(CO2(f,1),CO2(f,3),'.-');
   ylabel('[ppm]');
   title('pCO2w');
   fUnNaN=@qcWG_UnNaNcb;
   fNaN=@qcWG_NaNcb;
   ha(5).UserData=zeros(length(CO2(:,2)),1); 
   
   uicontrol('Units','normalized','Style','pushbutton','String','NaN',...
     'Position',[g{5}(1)+g{5}(3)-0.08 g{5}(2)+g{5}(4)+0.003 .08 .03],...
     'Callback',{fNaN,CO2(:,1),CO2(:,3),ha(5)});
 
   uicontrol('Units','normalized','Style','pushbutton','String','Un-NaN',...
     'Position',[g{5}(1)+g{5}(3)-0.16 g{5}(2)+g{5}(4)+0.003 .08 .03],...
     'Callback',{fUnNaN,CO2(:,1),CO2(:,3),ha(5)});

 axes(ha(6)); %pCO2a
   f = find(~isnan(CO2(:,4)));
   plot(CO2(f,1),CO2(f,4),'.-');
   ylabel('[ppm]');
   title('pCO2a');
   ha(6).UserData=zeros(length(CO2(:,2)),1);
   
   uicontrol('Units','normalized','Style','pushbutton','String','NaN',...
     'Position',[g{6}(1)+g{6}(3)-0.08 g{6}(2)+g{6}(4)+0.003 .08 .03],...
     'Callback',{fNaN,CO2(:,1),CO2(:,4),ha(6)});
 
   uicontrol('Units','normalized','Style','pushbutton','String','Un-NaN',...
     'Position',[g{6}(1)+g{6}(3)-0.16 g{6}(2)+g{6}(4)+0.003 .08 .03],...
     'Callback',{fUnNaN,CO2(:,1),CO2(:,4),ha(6)});
 
 axes(ha(7)); %pH
   f = find(~isnan(CO2(:,2)));
   plot(CO2(f,1),CO2(f,2),'.-');
   title('pH');
   ha(7).UserData=zeros(length(CO2(:,2)),1);
   
   uicontrol('Units','normalized','Style','pushbutton','String','NaN',...
     'Position',[g{7}(1)+g{7}(3)-0.08 g{7}(2)+g{7}(4)+0.003 .08 .03],...
     'Callback',{fNaN,CO2(:,1),CO2(:,2),ha(7)});
 
  uicontrol('Units','normalized','Style','pushbutton','String','Un-NaN',...
     'Position',[g{7}(1)+g{7}(3)-0.16 g{7}(2)+g{7}(4)+0.003 .08 .03],...
     'Callback',{fUnNaN,CO2(:,1),CO2(:,2),ha(7)});
 
 set(ha,'XLim',[min(CO2(:,1)) max(CO2(:,1))] + [-1 1]*0.1);
 set(ha(1:end-1),'XTickLabel',[]);
 linkaxes(ha,'x');

% (2) Flow-path pressure < 90 or > 104 kPa
% These thresholds tend to do a good job identifying leaks or plugs in 
%  the CO2 system flow path. I have been implementing this as a rule for
%  every WG mission. It is also worth manually checking Humidity for values 
%  greater than 1.5 V.
 fe = find( CO2equil(:,5) < 90 | CO2equil(:,5) > 104 );
 fz = find( CO2zero(:,5) < 90 | CO2zero(:,5) > 104 );
 fa = find( CO2air(:,5) < 90 | CO2air(:,5) > 104 );
 %index to 'CO2' matrix
  [~,fei] = intersect(CO2(:,1), CO2equil(fe,1));
  [~,fzi] = intersect(CO2(:,1), CO2zero(fz,1));
  [~,fai] = intersect(CO2(:,1), CO2air(fa,1));
  clear fe fz fa;
%check
 axes(ha(5));
 plot(CO2([fei fzi],1),CO2([fei fzi],3),'m.');
 axes(ha(6));
 plot(CO2([fai fzi],1),CO2([fai fzi],4),'m.');
%NaN bad data
 CO2([fei fzi],3) = NaN;
 CO2([fai fzi],4) = NaN;
 clear fei fzi fai;

% (4) T vs pH
% This is used to identify bad pH data by looking for outliers in
%  water-mass properties. 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.
 f = [];
 [c,ia,ic]=unique(CTD(:,1));
 x = interp1(CTD(ia,1),CTD(ia,2),CO2(:,1));
 y = CO2(:,2);
 figTpH = figure('Name','TpH','Units','inches','Position',[0 0 11 7.5]);
 plot(x,y,'.');
 hold on;
 xlabel(['T [' char(186) 'C]']);
 ylabel(['pH']);
 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',...
   {@PlotTpH,figTpH,figCO2,x,y,CO2,ha});
 uiwait(figTpH);
 clear f x y g gx gy;

% (5) pCO2w vs pH
% This is used to identify bad pCO2w and/or pH data by comparing to the
%  theoretical relation obtained from CO2SYS using an estimated alkalinity
%  derived from S. 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.
%calc theoretical relation

%plot
 figPP = figure('Name','PP','Units','inches','Position',[0 0 11 7.5]);
 scatter(CO2(:,3),CO2(:,2),15,'filled');
 hold on;
 plot(CO2SYS2100(:,1),CO2SYS2100(:,2),'-','Color',[.4 .4 .4],'LineWidth',2);
 plot(CO2SYS2400(:,1),CO2SYS2400(:,2),'k-','LineWidth',2);
 xlabel('pCO2 [ppm]');
 ylabel('pH');
 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',...
   {@plotPP,figPP,figCO2,CO2,ha});
 uiwait(figPP);
 clear f x y g gx gy xPred pH_th;

 sprintf('Paused, Hit any key to save data and continue')
 pause;
 for n1=5:7
     bool=zeros(size(CO2,1),1);
     for n=1:size(CO2,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
     
     if n1==5
        filtered=CO2(:,3);
        filtered(bool==1)=NaN;
        CO2(:,n1-2)=filtered; %POS 3
     elseif n1==6
         filtered=CO2(:,4);
         filtered(bool==1)=NaN;
         CO2(:,n1-2)=filtered; %POS 4
     elseif n1==7
         filtered=CO2(:,2);
         filtered(bool==1)=NaN;
        CO2(:,n1-5)=filtered; %POS 2
     end
 end

 close CO2 force

 CO2out=CO2;
 
% save QC'd .csv file
 time = (CO2(:,1) - dn) * 86400*1000;
 fid = fopen(savefile,'w');
 fprintf(fid, ['Time[GMT/Unix*1000],pCO2water[ppm],' ...
   'pCO2air[ppm],pH']);
 for n=1:size(CO2,1)
   fprintf(fid, '\n%.0f,%.3f,%.3f,%.3f', [time(n) CO2(n,2:4)]); 
 end
 fclose(fid);
 clear time fid n;

clear ha figCO2 figTpH figPP;

end