
% -- 
% Jeff Sevadjian
% Monterey Bay Aquarium Research Institute
% 7700 Sandholdt Rd
% Moss Landing, CA 95039, USA
% ph: +1 831-775-1888

% Human-in-the-loop QC procedure to generate final QC'd WG data set.
% Strongly recommend entering commands from this script manually a few
%  lines at a time as opposed to running the entire script all at once.

% All the 'processing' (applying cal's, deriving of engineering units, etc.)
%  is now done without human input in python scripts, either in
%  real-time (makeWGwebQC_wgTiny.py), or post-cruise (loadWGdata_wgTiny.py).
% This Matlab script is used only to identify and remove bad data.

% Outputs final .csv files in the /QC directory
% A final .nc file (merged to a common 10-min time stamp) needs to be
%  created after this script is run using makeWGnetCDF_wgTiny.py.

%depl info
 veh = 'wgTiny';
 dep = '20190114';

%directories for post-cruise and QC'd data
 atlas = '//atlas.shore.mbari.org/';
 mydirPC = [atlas 'WaveGlider/deployment_data/' veh '/' dep '/postCruise/'];
 mydirQC = [atlas 'WaveGlider/deployment_data/' veh '/' dep '/QC/'];
 clear veh dep;

%create /QC dir if doesn't exist
 if ~exist(mydirQC,'dir')
   mkdir(mydirQC);
 end

%paths to WG, M1 code
 addpath([atlas 'jsevadjian/code/WaveGlider/QC code/']);
 addpath([atlas 'jsevadjian/M-buoy/code/']);

%--%

%read-in postCruise .csv files
% NaN zeros, convert time stamp, interp position
 dn = datenum([1970 1 1]);
 %Position
  Position = csvread([mydirPC 'Position.csv'],1,0);
  Position(Position==0) = NaN;
  Position(:,1) = dn + Position(:,1)/86400/1000;
 %CO2 [pH, pCO2water, pCO2air]
  CO2 = csvread([mydirPC 'CO2.csv'],1,0);
  CO2(CO2==0) = NaN;
  CO2(:,1) = dn + CO2(:,1)/86400/1000;
  CO2(:,5) = interp1(Position(:,1),Position(:,2),CO2(:,1));
  CO2(:,6) = interp1(Position(:,1),Position(:,3),CO2(:,1));
 %CTD [T, S, O2]
  CTD = csvread([mydirPC 'CTD.csv'],1,0);
  CTD(CTD==0) = NaN;
  CTD(:,1) = dn + CTD(:,1)/86400/1000;
  CTD(:,5) = interp1(Position(:,1),Position(:,2),CTD(:,1));
  CTD(:,6) = interp1(Position(:,1),Position(:,3),CTD(:,1));
 %Eco-Puck [fluor, bb470, bb650]
  EP = csvread([mydirPC 'EcoPuck.csv'],1,0);
  EP(EP==0) = NaN;
  EP(:,1) = dn + EP(:,1)/86400/1000;
  EP(:,5) = interp1(Position(:,1),Position(:,2),EP(:,1));
  EP(:,6) = interp1(Position(:,1),Position(:,3),EP(:,1));
 %MET [airPress, airTemp]
  MET = csvread([mydirPC 'MET.csv'],1,0);
  MET(MET==0) = NaN;
  MET(:,1) = dn + MET(:,1)/86400/1000;
  MET(:,4) = interp1(Position(:,1),Position(:,2),MET(:,1));
  MET(:,5) = interp1(Position(:,1),Position(:,3),MET(:,1));
 %MET wind [windDir, windSpd]
  METw = csvread([mydirPC 'METwind.csv'],1,0);
  METw(METw==0) = NaN;
  METw(:,1) = dn + METw(:,1)/86400/1000;
  METw(:,4) = interp1(Position(:,1),Position(:,2),METw(:,1));
  METw(:,5) = interp1(Position(:,1),Position(:,3),METw(:,1));
 %TagID [tagID]
%   TagID = csvread([mydirPC 'tagID.csv'],1,0); %errors if no data
 %CO2 Equil
  CO2equil = csvread([mydirPC 'CO2equil.csv'],1,0);
  CO2equil(CO2equil==0) = NaN;
  CO2equil(:,1) = dn + CO2equil(:,1)/86400/1000;
 %CO2 Zero
  CO2zero = csvread([mydirPC 'CO2zero.csv'],1,0);
  CO2zero(CO2zero==0) = NaN;
  CO2zero(:,1) = dn + CO2zero(:,1)/86400/1000;
 %CO2 Air
  CO2air = csvread([mydirPC 'CO2air.csv'],1,0);
  CO2air(CO2air==0) = NaN;
  CO2air(:,1) = dn + CO2air(:,1)/86400/1000;
 %CO2 Stnd
  CO2stnd = csvread([mydirPC 'CO2stnd.csv'],1,0);
  CO2stnd(CO2stnd==0) = NaN;
  CO2stnd(:,1) = dn + CO2stnd(:,1)/86400/1000;
 %CO2 Hum
  CO2hum= csvread([mydirPC 'CO2hum.csv'],1,0);
  CO2hum(CO2hum==0) = NaN;
  CO2hum(:,1) = dn + CO2hum(:,1)/86400/1000;

%Plot all sensor data in one fig
% This is to compare different data sets for simultaneous 'events' -
%  i.e. if there is a spike in O2 and pCO2w at the same time, they 
%  are both likely to be 'real.' This is to aid removal of data which
%  is implemented in the individual sensor break-out figures.
 qcWG_plotAll(Position,CO2,CTD,EP,MET,METw);

%tBeg, tEnd should already be finalized in WG_cal_config.xlsx

%plot map
 coast = dlmread([atlas 'jsevadjian/GIS/montereyBay.txt']);
 figure; hold on;
 plot(coast(:,1),coast(:,2),'k');
 scatter(Position(:,2),Position(:,3),3,Position(:,1),'filled');
 colormap(jet);
 plot(-122.02203,36.75060,'k.','MarkerSize',10); %M1
 plot(-122.378,36.697,'k.','MarkerSize',10); %M2
 plot(-121.847,36.797,'k.','MarkerSize',10); %C1
 set(gca,'XLim',[min(Position(:,2)) max(Position(:,2))] + 0.2*[-1 1]);
 set(gca,'YLim',[min(Position(:,3)) max(Position(:,3))] + 0.2*[-1 1]);
 set(gca,'DataAspectRatio',[1 cos(nanmean(get(gca,'YLim')/180*pi)) 1]);
 box on;
 clear coast;

%--%

%load M1 data
 % See //atlas/jsevadjian/M-Buoy/M1_cal_config.xlsx
  M1dep = '201607'; %'201607' or '201507';
%   M1co2cal = [ ...
%     datenum([2015 7 30 0 0 0]), datenum([2015 10 1 20 0 0]), -12.7455, 0.4090, 0.9983, 0.0012; ...
%     datenum([2015 10 1 23 0 0]), datenum([2015 8 29 0 0 0]), -4.7736, 0.2689, 1.0355, -0.0023 ];
  M1co2cal = [ ...
    datenum([2016 8 29 0 0 0]), now, 3.5927 0.1272 1.0254 -0.003 ];
%   M1pHcal = [ ...
%     datenum([2015 7 30 19 0 0; 2015 10 23 18 50 0])', -0.378917, 284.99; ...
%     datenum([2015 10 23 19 0 0; 2016 8 29 19 0 0])', -0.384977, 279.15 ];
  M1pHcal = [ ...
    datenum([2016 8 29 0 0 0]), datenum([2017 3 14 18 40 0]), -0.404050 294.19; ...
    datenum([2017 3 14 18 40 1]), now, -0.400226, 293.81 ];

 %load data
  M1gps = read_M1_GPS_full('M1',M1dep);
  %CTD
   M1ctd = load_M1_CTD('M1',M1dep);
   M1ctd.lon = interp1(M1gps(:,1), M1gps(:,2), M1ctd.time);
   M1ctd.lat = interp1(M1gps(:,1), M1gps(:,3), M1ctd.time);
  %O2
   M1o2 = load_M1_O2('m1-slave1',M1dep,M1ctd);
   M1o2.lon = interp1(M1gps(:,1), M1gps(:,2), M1o2.time);
   M1o2.lat = interp1(M1gps(:,1), M1gps(:,3), M1o2.time);
  %pH
   M1pH = load_M1_pH('m1-slave1',M1dep,M1ctd,M1pHcal);
   M1pH.lon = interp1(M1gps(:,1), M1gps(:,2), M1pH.time);
   M1pH.lat = interp1(M1gps(:,1), M1gps(:,3), M1pH.time);
  %CO2 %NOW WITH STANDARD?
   M1co2 = load_M1_CO2('m1-slave1',M1dep,M1ctd,M1co2cal);
   M1co2.lon = interp1(M1gps(:,1), M1gps(:,2), M1co2.time);
   M1co2.lat = interp1(M1gps(:,1), M1gps(:,3), M1co2.time);
  clear M1dep M1gps M1*cal;

%index M1 data within 4 km, 10-60 min
 tolSp = 4; %km
 %CTD
  [ia,ib,indBad] = getIndexNear(CTD(:,1), CTD(:,5), CTD(:,6),...
    M1ctd.time, M1ctd.lon, M1ctd.lat, tolSp, 10/1440);
  CTD(:,7:8) = NaN(size(CTD,1),2); %T, S
  CTD(ia,7:8) = [M1ctd.temp(ib(ia)) M1ctd.sal(ib(ia))];
  CTD(indBad,7:8) = NaN;
 %O2
  [ia,ib,indBad] = getIndexNear(CTD(:,1), CTD(:,5), CTD(:,6),...
    M1o2.time, M1o2.lon, M1o2.lat, 2, 10/1440);
  CTD(:,9) = NaN(size(CTD,1),1); %O2
  CTD(ia,9) = M1o2.o2umolkg(ib(ia));
  CTD(indBad,9) = NaN;
 %pH
  [ia,ib,indBad] = getIndexNear(CO2(:,1), CO2(:,5), CO2(:,6),...
    M1pH.time, M1pH.lon, M1pH.lat, tolSp, 15/1440);
  CO2(:,7) = NaN(size(CO2,1),1); %pH
  CO2(ia,7) = M1pH.pH(ib(ia));
  CO2(indBad,7) = NaN;
 %CO2
  [ia,ib,indBad] = getIndexNear(CO2(:,1), CO2(:,5), CO2(:,6),...
    M1co2.time, M1co2.lon, M1co2.lat, tolSp, 60/1440);
  CO2(:,8:9) = NaN(size(CO2,1),2); %pCO2water, pCO2air
  CO2(ia,8:9) = [M1co2.pCO2water(ib(ia)) M1co2.pCO2air(ib(ia))];
  CO2(indBad,8:9) = NaN;
 clear tolSp ia ib indBad;

%--%

% %load underway pCO2 data
%  tBeg = datestr(CTD(1,1),'mm/dd/yyyy');
%  tEnd = datestr(CTD(end,1),'mm/dd/yyyy');
%  sql = [ 'SELECT date_time,dec_lat,dec_long,tmp,sal,pco2,tco2,sio4 '...
%    'FROM underway ' ...
%    'WHERE date_time>''' tBeg ''' AND ' ...
%          'date_time<''' tEnd ''' AND ' ...
%          'pco2>0 ' ...
%    'ORDER BY date_time' ]; 
%  c = database('jdbc:jtds:sqlserver://solstice.shore.mbari.org:1433/BOG',...
%    'everyone','guest');
%  v = jdbcquery(c, sql);
%  close(c);
%  clear tBeg tEnd sql c;
% 
% %index underway data within 4 km, 10-60 min
%  tolSp = 4; %km
%  %CO2
%   [ia,ib,indBad] = getIndexNear(CO2(:,1), CO2(:,5), CO2(:,6),...
%     v.date_time, v.dec_long, v.dec_lat, tolSp, 60/1440);
%   CO2(:,10) = NaN(size(CO2,1),1);
%   CO2(ia,10) = v.pco2(ib(ia));
%   CO2(indBad,10) = NaN;
%  clear tolSp ia ib indBad;

%------------%

%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).
 figCTD = figure('Name','CTD','Units','inches','Position',[0 0 11 8.5]);
 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),'.-');
%    plot(M1ctd.time,M1ctd.temp,'.-','Color',[.8 1 .8]);
%    plot(CTD(:,1),CTD(:,7),'.','Color',[0 .8 0],'MarkerSize',10);
%    le=legend('WG','M1');
%    set(le,'Color','none');
   ylabel(['[' char(186) 'C]']);
   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','CTD(:,2) = qcWG_NaNcb(CTD(:,1),CTD(:,2),ha(1));');
 axes(ha(2));
   plot(CTD(:,1),CTD(:,3),'.-');
%    plot(M1ctd.time,M1ctd.sal,'.-','Color',[.8 1 .8]);
%    plot(CTD(:,1),CTD(:,8),'.','Color',[0 .8 0],'MarkerSize',10);
   ylabel('');
   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','CTD(:,3) = qcWG_NaNcb(CTD(:,1),CTD(:,3),ha(2));');
 axes(ha(3));
   plot(CTD(:,1),CTD(:,4),'.-');
%    plot(M1o2.time,M1o2.o2umolkg,'.-','Color',[.8 1 .8]);
%    plot(CTD(:,1),CTD(:,9),'.','Color',[0 .8 0],'MarkerSize',10);
   ylabel(['[' char(181) 'mol kg^-^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','CTD(:,4) = qcWG_NaNcb(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 = slidingstdev(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.
 f = [];
 x = CTD(:,3);
 y = CTD(:,2);
 figTS = figure('Name','TS','Units','inches','Position',[0 0 11 8.5]);
 plot(x,y,'.');
 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',...
   ['gx=get(gca,''XLim'');' ...
    'gy=get(gca,''YLim'');' ...
    'f=[f; find(x>=gx(1) & x<=gx(2) & y>=gy(1) & y<=gy(2))];' ...
    'figure(figTS); plot(x(f),y(f),''k.'');' ...
    'figure(figCTD); axes(ha(2)); plot(CTD(f,1),CTD(f,3),''k.'');']);
 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.
 f = [];
 x = CTD(:,2);
 y = CTD(:,4);
 figTO = figure('Name','TO','Units','inches','Position',[0 0 11 8.5]);
 plot(x,y,'.');
 hold on;
%  plot(CTD(:,7),CTD(:,9),'g.'); %M1
 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',...
   ['gx=get(gca,''XLim'');' ...
    'gy=get(gca,''YLim'');' ...
    'f=[f; find(x>=gx(1) & x<=gx(2) & y>=gy(1) & y<=gy(2))];' ...
    'figure(figTO); plot(x(f),y(f),''k.'');' ...
    'figure(figCTD); axes(ha(3)); plot(CTD(f,1),CTD(f,4),''k.'');']);
 uiwait(figTO);
 clear f x y g gx gy;

%NaN O2 data where S is NaN
 CTD(isnan(CTD(:,3)),4) = NaN;

%save QC'd .csv file
 time = (CTD(:,1) - dn) * 86400*1000;
 fid = fopen([mydirQC 'CTD.csv'],'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;

close(figCTD);
clear ha le figCTD figTS figTO;

%------------%

%QC pCO2 and pH:

% (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 8.5]);
 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),'.-');
%    plot(CO2(:,1),CO2(:,8),'.','Color',[0 .8 0],'MarkerSize',10);
%    plot(CO2(:,1),CO2(:,10),'.','Color',[255 127 14]/255,'MarkerSize',10);
   ylabel('[ppm]');
   title('pCO2w');
   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','CO2(:,3) = qcWG_NaNcb(CO2(:,1),CO2(:,3),ha(5));');
 axes(ha(6)); %pCO2a
   f = find(~isnan(CO2(:,4)));
   plot(CO2(f,1),CO2(f,4),'.-');
%    plot(CO2(:,1),CO2(:,9),'.','Color',[0 .8 0],'MarkerSize',10);
   ylabel('[ppm]');
   title('pCO2a');
   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','CO2(:,4) = qcWG_NaNcb(CO2(:,1),CO2(:,4),ha(6));');
 axes(ha(7)); %pH
   f = find(~isnan(CO2(:,2)));
   plot(CO2(f,1),CO2(f,2),'.-');
%    plot(M1pH.time,M1pH.pH,'.-','Color',[.8 1 .8]);
%    plot(CO2(:,1),CO2(:,7),'.','Color',[0 .8 0],'MarkerSize',10);
%    le=legend('WG','M1');
%    set(le,'Color','none');
   title('pH');
   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','CO2(:,2) = qcWG_NaNcb(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 = [];
 x = interp1(CTD(:,1),CTD(:,2),CO2(:,1));
 y = CO2(:,2);
 figTpH = figure('Name','TpH','Units','inches','Position',[0 0 11 8.5]);
 plot(x,y,'.');
 hold on;
%  plot(CTD(:,7),CTD(:,9),'g.'); %M1
 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',...
   ['gx=get(gca,''XLim'');' ...
    'gy=get(gca,''YLim'');' ...
    'f=[f; find(x>=gx(1) & x<=gx(2) & y>=gy(1) & y<=gy(2))];' ...
    'figure(figTpH); plot(x(f),y(f),''k.'');' ...
    'figure(figCO2); axes(ha(7)); plot(CO2(f,1),CO2(f,2),''k.'');']);
 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
 sali = [33 34];
 tempi = [8 15];
 alk = 2150.0 + 44.0.*(sali-31.25); %(GF)
 alk = [2100 2200 2300 2400];
%CO2SYS
 sil = 0;
 po4 = 1;
 pH_th = [7.5:0.01:8.5];
 for n=1:length(alk)
%    A = CO2SYS(alk(n),pH_th,1,3,sali(n),tempi(n),NaN,0,NaN,...
%      sil,po4,1,4,1);
   A = CO2SYS(alk(n),pH_th,1,3,33.5,12,NaN,0,NaN,...
     sil,po4,1,4,1);
   xPred(:,n) = A(:,4); %(uatm)
 end
 clear sali tempi alk sil po4 n A;
%plot
 f = [];
 x = deNaN(CO2(:,1),CO2(:,3));
 y = CO2(:,2);
 figPP = figure('Name','PP','Units','inches','Position',[0 0 11 8.5]);
 plot(x,y,'.');
 hold on;
%  plot(CTD(:,8),CTD(:,7),'g.'); %M1
 plot(xPred(:,1),pH_th,'-','Color',[.4 .4 .4],'LineWidth',2);
 plot(xPred(:,2),pH_th,'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',...
   ['gx=get(gca,''XLim'');' ...
    'gy=get(gca,''YLim'');' ...
    'f=[f; find(x>=gx(1) & x<=gx(2) & y>=gy(1) & y<=gy(2))];' ...
    'figure(figPP); plot(x(f),y(f),''k.'');' ...
    'figure(figCO2); axes(ha(5)); plot(CO2(f,1),CO2(f,3),''k.'');' ...
    'axes(ha(7)); plot(CO2(f,1),CO2(f,2),''k.'');']);
 uiwait(figPP);
 clear f x y g gx gy xPred pH_th;

%save QC'd .csv file
 time = (CO2(:,1) - dn) * 86400*1000;
 fid = fopen([mydirQC 'CO2.csv'],'w');
 fprintf(fid, ['Time[GMT/Unix*1000],pH,pCO2water[ppm],' ...
   'pCO2air[ppm]']);
 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;

close(figCO2);
clear ha figCO2 figTpH figPP;

%------------%

%QC Eco-Puck:

% (1) Time-series plot
% This is used to manually NaN data, for example, if the Eco-Puck sensor
%  became fouled or had barnacles attach to sensor body.
 figEP = figure('Name','EP','Units','inches','Position',[0 0 11 8.5]);
 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(EP(:,1),EP(:,2),'.-');
   ylabel(['[' char(181) 'g L^-^1]']);
   title('Chlorophyll');
   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','EP(:,2) = qcWG_NaNcb(EP(:,1),EP(:,2),ha(1));');
 axes(ha(2));
   plot(EP(:,1),EP(:,3),'.-');
   ylabel('[m^-^1]');
   title('bb470');
   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','EP(:,3) = qcWG_NaNcb(EP(:,1),EP(:,3),ha(2));');
 axes(ha(3));
   plot(EP(:,1),EP(:,4),'.-');
   ylabel(['[m^-^1]']);
   title('bb650');
   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','EP(:,4) = qcWG_NaNcb(EP(:,1),EP(:,4),ha(3));');
 set(ha,'XLim',[min(EP(:,1)) max(EP(:,1))] + [-1 1]*0.1);
 linkaxes(ha,'x');

%save QC'd .csv file
 time = (EP(:,1) - dn) * 86400*1000;
 fid = fopen([mydirQC 'EcoPuck.csv'],'w');
 fprintf(fid, ['Time[GMT/Unix*1000],chlorophyll[ug/L],bb470[1/m],' ...
   'bb650[1/m]']);
 for n=1:size(EP,1)
   fprintf(fid, '\n%.0f,%.3f,%.5f,%.5f', [time(n) EP(n,2:4)]);
 end
 fclose(fid);
 clear time fid n;

close(figEP);
clear ha figEP g;

%------------%

%QC MET:

% (1) Time-series plot
% This is used to manually NaN data, for example, if the MET atm pressure
%  sensor started reporting bad data. Typical range for atm pressure is 
%  1000 - 1030 hPa (or mbar).
 figMET = figure('Name','MET','Units','inches','Position',[0 0 11 8.5]);
 ha = tight_subplot(2,1,0.075,0.2,0.1);
 g = get(ha,'Position');
 set(ha,'Box','on','NextPlot','add');
 axes(ha(1));
   plot(MET(:,1),MET(:,2),'.-');
%    plot(M1ctd.time,M1ctd.temp,'.-','Color',[.8 1 .8]);
%    plot(CTD(:,1),CTD(:,7),'.','Color',[0 .8 0],'MarkerSize',10);
%    le=legend('WG','M1');
%    set(le,'Color','none');
   ylabel(['[mbar]']);
   title('Atmospheric Pressure');
   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','MET(:,2) = qcWG_NaNcb(MET(:,1),MET(:,2),ha(1));');
 axes(ha(2));
   plot(MET(:,1),MET(:,3),'.-');
%    plot(M1ctd.time,M1ctd.sal,'.-','Color',[.8 1 .8]);
%    plot(CTD(:,1),CTD(:,8),'.','Color',[0 .8 0],'MarkerSize',10);
   ylabel(['[' char(186) 'C]']);
   title('Air Temperature');
   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','MET(:,3) = qcWG_NaNcb(MET(:,1),MET(:,3),ha(2));');
 set(ha,'XLim',[min(MET(:,1)) max(MET(:,1))] + [-1 1]*0.1);
 linkaxes(ha,'x');

%save QC'd .csv file
 time = (MET(:,1) - dn) * 86400*1000;
 fid = fopen([mydirQC 'MET.csv'],'w');
 fprintf(fid, ['Time[GMT/Unix*1000],airPress[hPa],airTemp[C]']);
 for n=1:size(MET,1)
   fprintf(fid, '\n%.0f,%.1f,%.1f', [time(n) MET(n,2:3)]);
 end
 fclose(fid);
 clear time fid n;

close(figMET);
clear ha figMET g;

%------------%

%QC METw:

% (1) Time-series plot
% This is used to manually NaN data, for example, if the MET wind
%  sensor started reporting bad data.
 figMETw = figure('Name','METw','Units','inches','Position',[0 0 11 8.5]);
 ha = tight_subplot(2,1,0.075,0.2,0.1);
 g = get(ha,'Position');
 set(ha,'Box','on','NextPlot','add');
 axes(ha(1));
   plot(METw(:,1),METw(:,2),'.');
%    plot(M1ctd.time,M1ctd.temp,'.-','Color',[.8 1 .8]);
%    plot(CTD(:,1),CTD(:,7),'.','Color',[0 .8 0],'MarkerSize',10);
%    le=legend('WG','M1');
%    set(le,'Color','none');
   ylabel(['[' char(181) 'True]']);
   title('Wind Direction');
   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','METw(:,2) = qcWG_NaNcb(METw(:,1),METw(:,2),ha(1));');
 axes(ha(2));
   plot(METw(:,1),METw(:,3),'.-');
%    plot(M1ctd.time,M1ctd.sal,'.-','Color',[.8 1 .8]);
%    plot(CTD(:,1),CTD(:,8),'.','Color',[0 .8 0],'MarkerSize',10);
   ylabel(['[m s^-^1]']);
   title('Wind Speed');
   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','METw(:,3) = qcWG_NaNcb(METw(:,1),METw(:,3),ha(2));');
 set(ha,'XLim',[min(METw(:,1)) max(METw(:,1))] + [-1 1]*0.1);
 linkaxes(ha,'x');

%save QC'd .csv file
 time = (METw(:,1) - dn) * 86400*1000;
 fid = fopen([mydirQC 'METwind.csv'],'w');
 fprintf(fid, ['Time[GMT/Unix*1000],windDir[True],windSpeed[m/s]']);
 for n=1:size(METw,1)
   fprintf(fid, '\n%.0f,%.1f,%.1f', [time(n) METw(n,2:3)]);
 end
 fclose(fid);
 clear time fid n;

close(figMETw);
clear ha figMETw g;

%------------%

%QC Position:

% (1) Time-series plot
% This is used to manually NaN data, for example, if the GPS started
%  reporting bad data.
% This has been updated to manually identify BOG Line data
 BogLine = zeros(size(Position(:,1)));
 figPos = figure('Name','Pos','Units','inches','Position',[0 0 11 8.5]);
 ha = tight_subplot(2,1,0.075,0.2,0.1);
 g = get(ha,'Position');
 set(ha,'Box','on','NextPlot','add');
 axes(ha(1));
   plot(Position(:,1),Position(:,2),'.-');
   ylabel(['[' char(186) 'E]']);
   title('Longitude');
   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','Position(:,2) = qcWG_NaNcb(Position(:,1),Position(:,2),ha(1));');
   uicontrol('Units','normalized','Style','pushbutton','String','BOG Line',...
     'Position',[g{1}(1)+g{1}(3)-0.2 g{1}(2)+g{1}(4)+0.003 .08 .03],...
     'Callback','BogLine = qcWG_BogLineCB(Position(:,1),Position(:,2),BogLine,ha(1));');
 axes(ha(2));
   plot(Position(:,1),Position(:,3),'.-');
   ylabel(['[' char(186) 'N]']);
   title('Latitude');
   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','Position(:,3) = qcWG_NaNcb(Position(:,1),Position(:,3),ha(2));');
 set(ha,'XLim',[min(Position(:,1)) max(Position(:,1))] + [-1 1]*0.1);
 linkaxes(ha,'x');

%save QC'd .csv file
 %Position
  time = (Position(:,1) - dn) * 86400*1000;
  fid = fopen([mydirQC 'Position.csv'],'w');
  fprintf(fid, ['Time[GMT/Unix*1000],longitude[E],latitude[N]']);
  for n=1:size(Position,1)
    fprintf(fid, '\n%.0f,%.5f,%.5f', [time(n) Position(n,2:3)]);
  end
  fclose(fid);
  clear time fid n;
 %BOG Line
  time = (Position(:,1) - dn) * 86400*1000;
  fid = fopen([mydirQC 'BogLine.csv'],'w');
  fprintf(fid, ['Time[GMT/Unix*1000],BogLine[N=0,Y=1]']);
  for n=1:size(Position,1)
    fprintf(fid, '\n%.0f,%.0f', [time(n) BogLine(n)]);
  end
  fclose(fid);
  clear time fid n;

%simply copy the Waypoint.csv file without any user-input
 copyfile([mydirPC 'Waypoint.csv'],[mydirQC 'Waypoint.csv']);

close(figPos);
clear ha figPos g;

%------------%

%Final Check:
 qcWG_plotAll(Position,CO2,CTD,EP,MET,METw);
