
% This will be run manually as a script. User-input and response
%  are an essential part of the processing!

% This will be good enough for the Wave Glider comparisons, but
%  one lingering issue/thing to check is, why are there some time
%  stamps not spaced exactly 1 hour apart,
%  e.g., plot(diff(ctd.time)*1440,'.');

% At some point in the future, a fully-documented NetCDF file
%  will be needed. For now, I am outputting .mat and .txt files only.

dep = '201407';

% Will define co2cal below, as we will first check whether the Licor
%  was ever switched (it wasn't).

addpath('//atlas.shore.mbari.org/jsevadjian/M-Buoy/code/');

%read & convert ctd data
 ctd = load_M1_CTD('m1',dep);
 ctd.sal(ctd.sal<15) = NaN; %first-order!

%read co2 data
 data = read_M1_CO2_full('m1-slave1',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)
% For initial QC plot only - will re-define later
 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);

%manually QC data

%axes param's
 ax.NextPlot = 'add';
 ax.Box = 'on';
 ax.XLim = [min(data(:,1)) max(data(:,1))] + [-3 3];

%plot
 myCol = jet(6);
 figure('Units','Inches','Position',[0 0 11 8.5]);
 ha = tight_subplot(3,1,0.03,0.1,[0.1 0.15]);
 set(ha,ax);
 %temp
  axes(ha(1));
  for n=1:6
    plot(data(f1+n-1,1),data(f1+n-1,2),'.-','Color',myCol(n,:));
  end
  ylabel({'Temperature'; '[C]'});
 %press
  axes(ha(2));
  for n=1:6
    plot(data(f1+n-1,1),data(f1+n-1,3),'.-','Color',myCol(n,:));
  end
  ylabel({'Pressure'; '[kPa]'});
 %xCO2
  axes(ha(3));
  for n=1:6
    plot(data(f1+n-1,1),data(f1+n-1,4),'.-','Color',myCol(n,:));
  end
  ylabel({'xCO2'; '[ppm]'});
 %labels, etc.
  linkaxes(ha,'x');
  for n=1:length(ha)
    le = legend(ha(n),{'Equil ON','Equil OFF','Zero ON','Zero OFF',...
      'Air ON','Air OFF'});
    ga=get(ha(n),'Position');
    gl = get(le,'Position');
    set(le,'Position',[ga(1)+ga(3)+0.01 gl(2:4)]);
  end
  clear myCol n le ga gl;

% Hit the water 7/16/2014
% First good data point is 19:13 GMT
% Can tell by the equil and air xCO2's diverging; equil press dropping

% % Was the Licor ever switched? Plot zero T vs CO2.
%  zeroInd = f1+3;
%  figure; hold on; plot(data(zeroInd,2),data(zeroInd,4),'.');
%  f=find(data(zeroInd,1)>datenum([2015 1 4 12 0 0]));
%  plot(data(zeroInd(f),2),data(zeroInd(f),4),'g.');

% Zero T vs CO2 does change ~mid deployment, but is gradual, not a sharp
%  change. No OASIS log notes about a change in the CO2 system.

%demarcate separate systems (if any) used for 201407 (ginput)
 oa4dates = datenum([2014 7 16 19 13 0; 2015 7 30 15 46 0])';

%Use these dates also to eliminate data before/after deployment
 f = find(data(:,1) < oa4dates(1) | data(:,1) > oa4dates(2));
 data(f,:) = [];
 clear f;

%Now have to re-index cycles
 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);

%enter calibrations from M1_cal_config.xlsx
 co2cal = [oa4dates -7.89500	0.32795	1.03580	-0.00229	1];

%define time according to EquilOn time
% round to nearest minute (cycle takes ~5 min to complete)
 dt = 1/1440;
 co2.time = round(data(f1,1)/dt) * dt;
 clear dt;

%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
  xCO2airFinal=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));
    xCO2airFinal(f)=( co2cal(n,3) + ...
      co2cal(n,4)*AirOffTemp(f) + ...
      co2cal(n,5)*xCO2airAdj(f) + ...
      co2cal(n,6)*AirOffTemp(f).*xCO2airAdj(f) ) * ...
      co2cal(n,7);
    xCO2waterFinal(f)=( co2cal(n,3) + ...
      co2cal(n,4)*EquilOffTemp(f) + ...
      co2cal(n,5)*xCO2waterAdj(f) + ...
      co2cal(n,6)*EquilOffTemp(f).*xCO2waterAdj(f) ) * ...
      co2cal(n,7);
  end
  clear Equil* Air* Zero* n f;

%convert xCO2 to pCO2
 %interp T, S to CO2 time indices
  tempi=interp1(ctd.time,ctd.temp,co2.time);
  sali=interp1(ctd.time,ctd.sal,co2.time);
 %check that you don't have too many NaNs from CTD
  goodInd = sum([~isnan(tempi)'; ~isnan(sali)'])==2;
  disp(length(goodInd) - sum(goodInd)); %2 out of 9087, great!
 %saturate the atmosphere with water vapor,
 % assuming 1 atm of pressure, at the in-situ T and S
 % **IS THIS THE SAME FOR AIR AND WATER? WHY A SALINITY TERM IN AIR?** %
  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.time(~goodInd)=[];
  co2.pCO2air = xCO2airFinal(goodInd).*VPFac(goodInd);
  co2.pCO2water = xCO2waterFinal(goodInd).*VPFac(goodInd);
  clear tempi sali goodInd TempK VP*;
 %check
  axes(ha(3));
  plot(co2.time,co2.pCO2air,'.-','Color',[.8 .8 .8]);
  plot(co2.time,co2.pCO2water,'.-','Color',[0 0 0]);
 clear xCO2*;

%save final data as .mat file
 mydir = '/Volumes/jsevadjian/M-Buoy/final CO2 data/';
 save([mydir dep '.mat'],'co2');

%write final data to .txt file
 dv = datevec(co2.time);
 fid = fopen([mydir dep '.txt'],'w+');
 fprintf(fid,'%s\n',['# GMT Year, Month, Day, Hr, Min, ' ...
   'Air pCO2 [ppm], Water pCO2 [ppm]']);
 fmt = [repmat('%02d%c',1,5) '%.2f%c%.2f\n'];
 for n=1:length(co2.time)
   fprintf(fid,fmt,...
     dv(n,1),',', dv(n,2),',', dv(n,3),',', dv(n,4),',', dv(n,5), ',', ...
     co2.pCO2air(n), ',', ...
     co2.pCO2water(n));
 end
 fclose(fid);
 clear dv fid fmt n;

%--%

% JS added 2017.03.20 to check against Gernot's daily-avg dpCO2 xls file

 u = unique(floor(co2.time));
 for n=1:length(u)
   f = find( floor(co2.time) == u(n) );
   pCO2wAvg(n) = nanmean(co2.pCO2water(f));
   pCO2aAvg(n) = nanmean(co2.pCO2air(f));
 end
 clear n f;

 dpCO2avg = [pCO2wAvg - pCO2aAvg]';