
function co2 = load_M1_CO2(buoy,dep,ctd,co2cal)

% buoy = 'm1-slave1';
% dep = '201507';
% ctd = load_M1_CTD('m1',dep);
% co2cal = [ datenum([2015 7 30 18 1 0]), datenum([2015 10 1 19 4 0]), ...
%   -12.74550	0.40896	0.99828	0.00122; ...
%   datenum([2015 10 1 23 34 0]), datenum(now+1), ...
%   -4.77360	0.26887	1.03550	-0.00227 ];

%read cumulative file
 data = read_M1_CO2_full(buoy,dep);

%get rid of nonsense time values
 f = find( data(:,1)<datenum([1900 1 0]) | ...
   data(:,1)>now +30 );
 data(f,:)=[];
 clear f;

%take unique time values only
 [~,i]=unique(data(:,1));
 data=data(i,:);
 clear i;

%sort by time (may not be nec)
 [~,i]=sort(data(:,1));
 data=data(i,:);
 clear i;

%index cycles (equil, zero, air; on, off)
 f1=find(data(1:end-5,end)==1 & data(2:end-4,end)==2 & ...
   data(3:end-3,end)==3 & data(4:end-2,end)==4 & ...
   data(5:end-1,end)==5 & data(6:end,end)==6);

%assign to structure
 co2.time = data(f1,1);

%apply calibrations
 %define data needed for cal
  EquilOffCO2=data(f1+1,4);
  EquilOffTemp=data(f1+1,2);
  ZeroOffCO2=data(f1+3,4);
  AirOffCO2=data(f1+5,4);
  AirOffTemp=data(f1+5,2);
  clear f1 dt;

 %subtract zeros
  xCO2airAdj = AirOffCO2-ZeroOffCO2;
  xCO2waterAdj = EquilOffCO2-ZeroOffCO2;

 %apply T-curve fits
  co2.pCO2air=NaN(size(co2.time));
  xCO2waterFinal=NaN(size(co2.time));
  for n=1:size(co2cal,1)
    f=find(co2.time>=co2cal(n,1) & co2.time<=co2cal(n,2));
    co2.pCO2air(f) = co2cal(n,3) + ...
      co2cal(n,4)*AirOffTemp(f) + ...
      co2cal(n,5)*xCO2airAdj(f) + ...
      co2cal(n,6)*AirOffTemp(f).*xCO2airAdj(f);
    xCO2waterFinal(f) = co2cal(n,3) + ...
      co2cal(n,4)*EquilOffTemp(f) + ...
      co2cal(n,5)*xCO2waterAdj(f) + ...
      co2cal(n,6)*EquilOffTemp(f).*xCO2waterAdj(f);
  end
  clear Equil* Air* Zero* n f;

%convert xCO2 to pCO2 (pCO2water only)
 %interpolate ctd (T and S) to o2
  for n=1:length(co2.time)
    [m(n),i(n)] = min(abs(co2.time(n) - ctd.time));
  end
  tempi = ctd.temp(i);
  sali = ctd.sal(i);
  f = find(m > 15/1440);
  tempi(f) = NaN;
  sali(f) = NaN;
  clear n m i f;
 %saturate the atmosphere with water vapor,
 % assuming 1 atm of pressure, at the in-situ T and S
  TempK = tempi + 273.15;
  VPWP = exp(24.4543 - 67.4509 * (100./TempK) - 4.8489 * ...
    log(TempK./100));
  VPCorrWP = exp(-0.000544 * sali);
  VPSWWP = VPWP .* VPCorrWP;
  VPFac = 1 - VPSWWP;
 %define final data
  co2.pCO2water = xCO2waterFinal.*VPFac;
  clear tempi sali TempK VP*;
