
function o2 = load_M1_O2(buoy,dep,ctd)

% buoy = 'm1-slave1';
% dep = '201507';
% ctd = load_M1_CTD('m1',dep);

%read cumulative file
 data = read_M1_O2_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;

%assign to structure
 o2.time = data(:,1);
 o2.o2uMS0 = data(:,2);
 clear data;

%interpolate ctd (T and S) to o2
 for n=1:length(o2.time)
   [m(n),i(n)] = min(abs(o2.time(n) - ctd.time));
 end
 o2temp = ctd.temp(i);
 o2sal = ctd.sal(i);
 o2sigt = ctd.sigt(i);
 f = find(m > 15/1440);
 o2temp(f) = NaN;
 o2sal(f) = NaN;
 o2sigt(f) = NaN;
 clear n m i f;

%reference O2 conc to seawater (S=33.5 psu)
 B0=-6.24097e-3;
 B1=-6.93498e-3;
 B2=-6.90358e-3;
 B3=-4.29155e-3;
 C0=-3.11680e-7;
 S=o2sal;
 Ts = log( (298.15 - o2temp) ./ (273.15 + o2temp));
 o2.o2uM = o2.o2uMS0 .* exp( S.*(B0 + B1*Ts + B2*(Ts.^2) + ...
   B3*(Ts.^3)) + C0*(S.^2) );
 clear B0 B1 B2 B3 C0 S Ts o2uMS0;

%convert uM to umol/kg
 o2.o2umolkg = o2.o2uM*1000./(o2sigt+1000);
 clear ctd;
