% fixit_da1_FATIMA_2022.m
% This routine eliminates known bad segments of raw 10-min data

% Don't run this script directly.  It is called from a3_final_filtered_FATIMA_2022.m

% You can edit the lines below to make crude corrections to raw data.

disp('fixit_da1_FATIMA_2022')

% define crude wdir filter for bad periods (e.g. wind from stern sector))
ee = abs(da.rwdir_psl)>120; % +/- 120° wind sector

% filter for bad periods of ship met data
% there's no ship data for FATIMA, so set this entire flag to 'true'
pp = true(height(da),1);

%% wind speed
% remove known bad wind data...
da.wspd_ship(pp) = NaN;
da.wspd_psl(ee) = NaN; da.wspd_ship(ee) = NaN; da.wspd_wxt(ee) = NaN;

figure('position',[1,1,800,475],'name','wspd');
plot(dt,da.wspd_psl,'b.',dt,da.wspd_ship,'r.',dt,da.wspd_wxt,'g.','MarkerSize',8);
legend('GillR3','Ship','WXT','location','best'); grid on; grid minor;
ylabel('True Wspd m/s'); title([ptitle,' True Wind Speed (m/s)']);
if savefigs
    ppath = fullfile(png_path,['001_',cruise,'_True_Wind_Spd.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['001_',cruise,'_True_Wind_Spd.fig']));
end

[xbins, ybins1, ~, SD1, ~] = binave2(da.wspd_psl,da.wspd_ship,2,0,18);
[~, ybins0, ~, SD0, ~] = binave2(da.wspd_psl,da.wspd_wxt,2,0,18);
figure('name','wspd_comp'); plot([0,20],[0,20],'k-'); hold on;
p1 = errorbar(xbins,ybins1,SD1,'rd-','LineWidth',1); grid on; grid minor;
p2 = errorbar(xbins,ybins0,SD0,'bo-','LineWidth',1); hold  off;
xlim([0 20]); ylim([0 20]);
legend([p1 p2],{'Ship','WXT'},'location','southeast');
ylabel('Wspd, various m/s'); xlabel('Wspd, NOAA sonic, m/s');
title([ptitle,' True Wind Speed Comparison vs Gill R3 Sonic']);
if savefigs
    ppath = fullfile(png_path,['002_',cruise,'_True_Wind_Spd_Comparison.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['002_',cruise,'_True_Wind_Spd_Comparison.fig']));
end

figure('position',[1,1,800,475],'name','sonic_qc');
subplot(2,1,1);
    plot(dt,da.missing_son,'b.','markersize',8); grid on; grid minor;
    ylabel('missing pnts'); ylim([-100,12500]);
    title([ptitle,' Missing Sonic Data Points per 10min @ ',...
                  sprintf('%02i',da.fsonic(1)),' Hz']);
subplot(2,1,2);
    plot(dt,da.badPnts_son,'b.','markersize',8); grid on; grid minor;
    ylabel('bad pnts');  ylim([-100,12500]);
    title([ptitle,' Bad Sonic Data Points per 10min @ ',...
                  sprintf('%02i',da.fsonic(1)),' Hz']);
if savefigs
    ppath = fullfile(png_path,['003_',cruise,'_Sonic_QA.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['003_',cruise,'_Sonic_QA.fig']));
end

%% wind direction
% remove bad wdir, as above
da.wdir_ship(pp) = NaN;
da.wdir_psl(ee) = NaN; da.wdir_ship(ee) = NaN;

figure('position',[1,1,800,475],'name','wdir');
plot(dt,da.wdir_psl,'b.',dt,da.wdir_ship,'r.',dt,da.wdir_wxt,'g.','markersize',8);
ylim([0,360]); yticks([0,90,180,270,360]); grid on; grid minor;
ylabel('True Wdir, deg');
legend('GillR3','Ship','WXT','location','best');
title([ptitle,' True Wind Direction (relative to earth)']);
if savefigs
    ppath = fullfile(png_path,['004_',cruise,'_True_Wind_Dir.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['004_',cruise,'_True_Wind_Dir.fig']));
end

figure('name','wdir_hist'); histogram(da.wdir_psl,0:15:360); grid;
xlim([0,360]); xticks(0:60:360); xlabel('True Wind Direction');
ylabel('Frequency Histogram, 10-min averages');
title([ptitle,' Wind Direction Frequency Historgam, NOAA sonic']);
if savefigs
    ppath = fullfile(png_path,['005_',cruise,'_True_Wind_Dir_Dist.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['005_',cruise,'_True_Wind_Dir_Dist.fig']));
end

%% relative wind direction
da.rwdir_ship(pp) = NaN;

figure('position',[1,1,800,475],'name','rwdir');
plot([stdt,endt],[90,90],'k--',[stdt,endt],[-90,-90],'k--','LineWidth',1); hold on;
plot([stdt,endt],[60,60],'r--',[stdt,endt],[-60,-60],'r--','LineWidth',1);
plot([stdt,endt],[0,0],'k-');
p1 = plot(dt,da.rwdir_psl,'b.',dt,da.rwdir_ship,'r.','markersize',8);hold off;
ylim([-180,180]); yticks([-180,-90,0,90,180]); grid on; grid minor;
ylabel('Rel Wdir, deg');
% legend(p1,{'GillR3','Ship'},'location','best');
title([ptitle,' Relative Wind Direction']);
if savefigs
    ppath = fullfile(png_path,['006_',cruise,'_Rel_Wind_Dir.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['006_',cruise,'_Rel_Wind_Dir.fig']));
end

figure('name','rwdir_hist'); histogram(da.rwdir_psl,-180:15:180); grid on; grid minor;
xlim([-180,180]); xticks(-180:45:180); xlabel('Relative Wind Direction');
ylabel('Frequency Histogram, 10-min averages');
title([ptitle,' Relative Wind Direction Frequency Historgam, NOAA sonic']);
if savefigs
    ppath = fullfile(png_path,['007_',cruise,'_Rel_Wind_Dir_Dist.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['007_',cruise,'_Rel_Wind_Dir_Dist.fig']));
end

%% sea level pressure
da.Pmb_ship(pp) = NaN;

figure('position',[1,1,800,475],'name','pressure');
plot(dt,da.Pmb_psl,'b.',dt,da.Pmb_ship,'r.',dt,da.Pmb_wxt,'g.','markersize',8);
ylabel('Sea Level pressure, mb'); grid on; grid minor;
legend('PTB','Ship','WXT','location','best');
title([ptitle,' Sea Level Pressure']);
if savefigs
    ppath = fullfile(png_path,['008_',cruise,'_Pressure.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['008_',cruise,'_Pressure.fig']));
end

[bpbinNOAA, bpbinSCS, ~, SDbpSCS, ~] = binave2(da.Pmb_psl,da.Pmb_ship,2,998,1012);
[~, bpbinWXT, ~, SDbpWXT, ~] = binave2(da.Pmb_psl,da.Pmb_wxt,2,998,1012);
figure('name','press_comp'); plot([998,1018],[998,1018],'k-'); hold on;
h1 = errorbar(bpbinNOAA,bpbinSCS,SDbpSCS,'rd-','LineWidth',1);
h2 = errorbar(bpbinNOAA,bpbinWXT,SDbpWXT,'gd-','LineWidth',1); grid on; grid minor;
legend([h1,h2],'Ship','WXT','location','best'); % axis([1000,1022,1006,1022]);
ylabel('Sea Level BP, mb'); xlabel('Sea Level BP, NOAA PTB, mb');
title([ptitle,' Pressure Comparison vs PTB220']);
if savefigs
    ppath = fullfile(png_path,['009_',cruise,'_Pressure_Comparison.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['009_',cruise,'_Pressure_Comparison.fig']));
end

disp(['Mean Diff BP, PTB - WXT = ',sprintf('%6.2f',nanmean1(da.Pmb_psl-da.Pmb_wxt))]);

%% tilt angle - flow distortion angle derived from the motion correction script
% this is tilt angle prior to ship speed correction

figure('position',[1,1,800,475],'name','tilt_rwdir');
plot(da.rwdir_psl,da.tilt,'b.',[-60 -60],[-5 25],'k--',[60 60],[-5 25],'k--',...
      [-180 180],[0 0],'k-','MarkerSize',8);
axis([-180 180 -5 25]); grid on; grid minor;
ylabel('Streamline Tilt Angle, deg'); xlabel('Relative Wind Direction');
title([ptitle,' Wind Streamline Tilt Angle, NOAA sonic']);
if savefigs
    ppath = fullfile(png_path,['010_',cruise,'_Tilt.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['010_',cruise,'_Tilt.fig']));
end

[xbins, ybins, ~, SD, ~] = binave2(da.rwdir_psl,da.tilt,20,-180,180);
figure('position',[1,1,800,475],'name','tilt_rwdir_binned');
plot([-60 -60],[-5 25],'k--',[60 60],[-5 25],'k--',[-180 180],[0 0],'k-'); hold on;
errorbar(xbins,ybins,SD,'rd-','LineWidth',1); hold off; axis([-180 180 -5 25]);
ylabel('Streamline Tilt Angle, deg'); xlabel('Relative Wind Direction');
title([ptitle,' Binned Streamline Tilt Angle, NOAA Sonic']);  grid on; grid minor;
if savefigs
    ppath = fullfile(png_path,['011_',cruise,'_Tilt_Binned.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['011_',cruise,'_Tilt_Binned.fig']));
end

[xbins, ybinswv, ~, ~, ~] = binave2(da.rwdir_psl,da.wvar,20,-180,180);
[~, ybinsusb, ~, ~, ~] = binave2(da.rwdir_psl,(da.usb_psl .* psi_wvar(da.zet_psl)),20,-180,180);
sigW_ustar = sqrt(ybinswv)./ybinsusb;
figure('position',[1,1,800,475],'name','sigW_usr_binned');
plot([-75 -75],[-5 25],'k--',[75 75],[-5 25],'k--',[-180 180],[0 0],'k-',...
     [-180 180],[1.25 1.25],'b--'); hold on;
plot(xbins,sigW_ustar,'rd-','MarkerSize',10,'LineWidth',1.5); hold off;
axis([-180,180,-1,5.0]); ylabel('\sigma_w / (u_* x psi(z/L))'); 
xlabel('Relative Wind Direction'); grid on; grid minor;
title([ptitle,' Sigma W/ustar vs Rel Wind Direction']);
if savefigs
    ppath = fullfile(png_path,['012_',cruise,'_SigWusr_Binned.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['012_',cruise,'_SigW_Binned.fig']));
end

%% sea surface temperature
% filter out bad rosr data
aa = false(height(da),1);  % ROSR is all good for now
da.sst_rosr(aa) = NaN;
bb = false(height(da),1);  % bad snake data already removed in met reader
da.sst_psl(bb) = NaN;
cc = true(height(da),1);   % no ship data, so set to all bad
da.sst_ship(cc) = NaN;

maxx = max(da.sst_rosr)+0.5; minx = min(da.sst_rosr)-0.5;
figure('position',[1,1,800,800],'name','sst');
subtightplot(2,1,1,[0.08,0.08],[0.05,0.05],[0.08,0.05]);
    plot(dt,da.sst_psl,'r-',dt,da.sst_ship,'b-',dt,da.sst_rosr,...
         'g-','LineWidth',1);
    ylim([minx,maxx]); grid on; grid minor;
    ylabel('SST, °C'); title([ptitle,' Sea Surface Temperature']);
    legend('Sea Snake','Ship','ROSR','location','best');
subtightplot(2,1,2,[0.08,0.08],[0.05,0.05],[0.08,0.05]);
    histogram(da.sst_rosr-da.sst_psl,[-1.5:0.1:1.5]); grid on; grid minor;
    xlabel('SST-rosr - SST-psl'); ylabel('count, 10-min periods');
if savefigs
    ppath = fullfile(png_path,['013_',cruise,'_Tsea.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['013_',cruise,'_Tsea.fig']));
end

ii = isfinite(da.sst_psl);
[xbins0, ybins0, ~, SD0, ~] = binave2(da.hour(ii),da.sst_psl(ii),2,0,24);
[xbins1, ybins1, ~, SD1, ~] = binave2(da.hour(ii),da.sst_ship(ii),2,0,24);
[xbins2, ybins2, ~, SD2, ~] = binave2(da.hour(ii),da.sst_rosr(ii),2,0,24);
figure('position',[1,1,800,475],'name','sst_diel');
p1 = errorbar(xbins0,ybins0,SD0,'rd-','LineWidth',1); hold on;
p2 = errorbar(xbins1,ybins1,SD1,'bo-','LineWidth',1);
p3 = errorbar(xbins2,ybins2,SD2,'gs-','LineWidth',1); hold off
xticks(0:1:24); xlabel('Hour of Day, UTC'); grid on; grid minor;
legend([p1 p2 p3],{'Sea Snake','Ship','ROSR'},'location','best');
ylabel('SST, °C'); title([ptitle,' SST Diel Cycle']);
if savefigs
    ppath = fullfile(png_path,['014_',cruise,'_Tsea_Diel.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['014_',cruise,'_Tsea_Diel.fig']));
end

%% relative humidity
da.rh_ship(pp) = NaN;

figure('position',[1,1,800,475],'name','rh');
plot(dt,da.rh_psl,'b.',dt,da.rh_ship,'r.',dt,da.rh_wxt,'g.','MarkerSize',8);
ylabel('RH, %'); grid on; grid minor;
% legend('HMT','Ship','WXT','location','best');
title([ptitle,'  Relative Humidity']);
if savefigs
    ppath = fullfile(png_path,['015_',cruise,'_RH.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['015_',cruise,'_RH.fig']));
end

%% air temperature
da.ta_ship(pp) = NaN;

figure('position',[1,1,800,475],'name','ta');
plot(dt,da.ta_psl,'b.',dt,da.ta_ship,'r.',dt,da.ta_wxt,'g.','MarkerSize',8);
ylabel('T air, °C'); grid on; grid minor;
% legend('HMT','Ship','WXT','location','best');
title([ptitle,'  Air temperature']);
if savefigs
    ppath = fullfile(png_path,['016_',cruise,'_Tair.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['016_',cruise,'_Tair.fig']));
end

%% specific humidity
% recompute q in case T/RH has changed
da.qa_psl = qair_p([da.ta_psl,da.rh_psl,da.Pmb_psl]);
da.qa_ship = qair_p([da.ta_ship,da.rh_ship,da.Pmb_ship]);
da.qa_wxt = qair_p([da.ta_wxt,da.rh_wxt,da.Pmb_wxt]);

figure('position',[1,1,800,475],'name','qa');
plot(dt,da.qa_psl,'b.',dt,da.qa_ship,'r.',dt,da.qa_wxt,'g.','MarkerSize',8); 
ylabel('Q air, g/kg'); grid on; grid minor;
% legend('HMT','Ship','WXT','location','best');
title([ptitle,' Air specific humidity']);
if savefigs
    ppath = fullfile(png_path,['017_',cruise,'_Qa.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['017_',cruise,'_Qa.fig']));
end

%% surface saturation specific humidity
% recompute qs in case SST has changed
da.qs_psl = qsea_p([da.sst_psl,da.Pmb_psl]);
da.qs_ship = qsea_p([da.sst_ship,da.Pmb_ship]);
da.qs_rosr = qsea_p([da.sst_rosr,da.Pmb_psl]);

figure('position',[1,1,800,475],'name','qs');
plot(dt,da.qs_psl,'b.',dt,da.qs_ship,'r.',dt,da.qs_rosr,'g.','MarkerSize',8);
ylabel('Qsea, g/kg'); grid on; grid minor;
legend('Sea Snake','Ship','ROSR','location','best');
title([ptitle,' Saturation Specific humidity at Ocean Surface (g/kg)']);
if savefigs
    ppath = fullfile(png_path,['018_',cruise,'_Qsea.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['018_',cruise,'_Qsea.fig']));
end

%% solar radiative flux
gg = true(height(da),1);  % no ship solar flux
da.rs_ship(gg) = NaN;

figure('position',[1,1,800,475],'name','rs');
plot(dt,da.rs_psl,'b.',dt,da.rs_ship,'r.','MarkerSize',6);
ylabel('Solar Flux, W/m^2'); grid on; grid minor;
legend('PSD','Ship','location','best');
title([ptitle,' Downwelling Solar Radiation']);
if savefigs
    ppath = fullfile(png_path,['019_',cruise,'_SolarRad.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['019_',cruise,'_SolarRad.fig']));
end

%% longwave radiative flux
gg = true(height(da),1);  % no ship IR flux
da.rl_ship(gg) = NaN;

figure('position',[1,1,800,475],'name','rl');
plot(dt,da.rl_psl,'b.',dt,da.rl_ship,'r.','MarkerSize',6);
ylabel('IR Flux, W/m^2'); grid on; grid minor;
legend('PSD','Ship','location','best');
title([ptitle,' Downwelling IR Radiation']);
if savefigs
    ppath = fullfile(png_path,['020_',cruise,'_IR_Rad.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['020_',cruise,'_IR_Rad.fig']));
end

%% rain rate
figure('position',[1,1,800,475],'name','rain');
plot(dt,da.rain_psl,'r.',dt,da.rain_wxt,'b.','MarkerSize',8);
ylabel('Rain, mm/hr'); grid on; grid minor;
legend('ORG','WXT','location','best');
title([ptitle,' Rainrate']);
if savefigs
    ppath = fullfile(png_path,['021_',cruise,'_Rain.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['021_',cruise,'_Rain.fig']));
end

%% SOG/COG/Heading/Map plots
figure('position',[1,1,800,475],'name','sog');
plot(dt,da.sog_psl,'b.',dt,da.sog_ship,'r.','MarkerSize',6);
ylabel('SOG, m/s'); grid on; grid minor;
legend('NOAA','Ship','location','best');
title([ptitle,' Ship Speed Over Ground']);
if savefigs
    ppath = fullfile(png_path,['022_',cruise,'_SOG.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['022_',cruise,'_SOG.fig']));
end

figure('position',[1,1,800,475],'name','cog');
plot(dt,da.cog_psl,'b.',dt,da.cog_ship,'r.');
yticks([0,45,90,135,180,225,270,315,360]); ylim([0,360]);
ylabel('COG, deg'); grid on; grid minor;
legend('NOAA','Ship','location','best');
title([ptitle,' Ship Course Over Ground']);
if savefigs
    ppath = fullfile(png_path,['023_',cruise,'_COG.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['023_',cruise,'_COG.fig']));
end

figure('position',[1,1,800,475],'name','hed');
plot(dt,da.hed_psl,'b.',dt,da.lrg,'r.','MarkerSize',6);
yticks([0,45,90,135,180,225,270,315,360]); ylim([0,360]);
ylabel('Heading,deg'); grid on; grid minor;
legend('NOAA','Ship','location','best');
title([ptitle,' Heading']);
if savefigs
    ppath = fullfile(png_path,['024_',cruise,'_Heading.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['024_',cruise,'_Heading.fig']));
end

figure('position',[1,1,800,475],'name','std_hed');
plot(dt,da.hed_psl_std,'b.','MarkerSize',6);
ylabel('Std Dev Heading, deg'); grid on; grid minor;
title([ptitle,' Standard Deviation Heading']);
if savefigs
    ppath = fullfile(png_path,['025a_',cruise,'_StdDev_Heading.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['025a_',cruise,'_StdDev_Heading.fig']));
end

figure('position',[1,1,800,475],'name','std_sog');
plot(dt,da.sog_psl_std,'b.','MarkerSize',6);
ylabel('Std Dev SOG, m/s'); grid on; grid minor;
title([ptitle,' Standard Deviation SOG']);
if savefigs
    ppath = fullfile(png_path,['025b_',cruise,'_StdDev_SOG.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['025b_',cruise,'_StdDev_SOG.fig']));
end

map_FATIMA_2022(da.lon_psl,da.lat_psl,LatLonLim)
if savefigs
    ppath = fullfile(png_path,['026_',cruise,'_TrackMap.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['026_',cruise,'_TrackMap.fig']));
end

%% Licor fast humidity & co2
% lots of clean up required here...

% wind from stern of ship or AGC is bad or CO2 is bad or points missing
% could add more here based on std devs, etc.
xx = ee | da.AGC_lic>AGC_lim | da.co2_lic_all<CO2_LoLim | ...
     da.co2_lic_all>CO2_HiLim | da.missing_lic > da.flicor*600*0.05;
da.h2o_lic_all(xx) = NaN; da.h2o_lic_good(xx) = NaN;
da.co2_lic_good(xx) = NaN;
da.wq_cov(xx) = NaN; da.wq_csp(xx) = NaN;
da.hl_cov(xx) = NaN; da.hl_csp(xx) = NaN;
da.hlc_cov(xx) = NaN; da.hlc_csp(xx) = NaN;
da.Cq2a(xx) = NaN; da.Cq2b(xx) = NaN;

% remove remaining outlier wq covariances
da.wq_cov = da.wq_cov*1.045; % wq Licor with correction for sensor sep: 1.045
da.wq_csp = da.wq_csp*1.045;
ii = da.wq_cov<wq_LoLim | da.wq_cov>wq_HiLim;
da.wq_cov(ii) = NaN;
ii = da.wq_csp<wq_LoLim | da.wq_csp>wq_HiLim;
da.wq_csp(ii) = NaN;

% recompute licor q from filtered data
da.q_lic = (da.h2o_lic_all*1e-3).*Mw./da.rhoa_psl;

 % q licor plot
figure('position',[1,1,800,475],'name','q_lic_agc');
subplot(2,1,1); % q licor plot
plot(dt,da.qa_wxt,'b-',dt,da.q_lic,'r-','LineWidth',1.5);
ylabel('Qair, g/kg'); grid on; grid minor;
title([ptitle,' Licor Specific Humidity and AGC (red dashed = limit)']);
legend('WXT','Licor','location','northwest');
subplot(2,1,2); % Licor AGC
plot(dt,da.AGC_lic,[stdt,endt],[AGC_lim,AGC_lim],'r--',[stdt,endt],...
     [50,50],'k-','LineWidth',1.5);
ylabel('AGC (arbitrary units)'); ylim([48,100]); grid on; grid minor;
legend('AGC','Limit','location','northwest');
if savefigs
    ppath = fullfile(png_path,['027_',cruise,'_Licor_Qa_AGC.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['027_',cruise,'_Licor_Qa_AGC.fig']));
end

% Licor CO2 (ship plume indicator)
figure('position',[1,1,800,475],'name','co2_lic');
plot(dt,da.co2_lic_all,[stdt,endt],[CO2_LoLim,CO2_LoLim],'r--',[stdt,endt],...
     [CO2_HiLim,CO2_HiLim],'r--','LineWidth',1.5);
grid on; grid minor; ylabel('CO_2 (ppm)'); ylim([0,1000]);
title([ptitle,' LICOR CO2, Contamination Indicator, red dashed = limits']);
if savefigs
    ppath = fullfile(png_path,['028_',cruise,'_Licor_CO2.png']);
    print('-dpng',ppath);
%     saveas(gcf,fullfile(fig_path,['028_',cruise,'_Licor_CO2.fig']));
end

close all

%% filter covariances and ID variables for outlier values

% first, apply crude wind sector and missing son/nav data criteria.
% discard periods with more than 5% missing son, mot or hed
kk = ee & (da.missing_mot>da.fmotion*600*0.05 |...
           da.missing_hed>da.fheading*600*0.05 | ...
           da.missing_son>da.fsonic*600*0.05 |...
           da.badPnts_son>badSon_lim);

da.wu_cov(kk) = NaN; da.wu_csp(kk) = NaN; da.wv_cov(kk) = NaN; da.wv_csp(kk) = NaN;
da.wTs_cov(kk) = NaN; da.wTs_csp(kk) = NaN; da.wT_cov(kk) = NaN; da.wT_csp(kk) = NaN;
da.wq_cov(kk) = NaN; da.wq_csp(kk) = NaN; da.hs_cov(kk) = NaN; da.hs_csp(kk) = NaN;
da.hl_cov(kk) = NaN; da.hl_csp(kk) = NaN; da.hlc_cov(kk) = NaN; da.hlc_csp(kk) = NaN;
da.Cu2a(kk) = NaN; da.Cu2b(kk) = NaN; da.Cw2a(kk) = NaN; da.Cw2b(kk) = NaN;
da.Ct2a(kk) = NaN; da.Ct2b(kk) = NaN; da.Cq2a(kk) = NaN; da.Cq2b(kk) = NaN;
da.usid_a(kk) = NaN; da.usid_b(kk) = NaN; da.tsid_a(kk) = NaN; da.tsid_b(kk) = NaN;
da.qsid_a(kk) = NaN; da.qsid_b(kk) = NaN; da.Lid_a(kk) = NaN; da.Lid_b(kk) = NaN;

% apply general covariance outlier limits
ii = da.wu_cov<wu_LoLim | da.wu_cov>wu_HiLim; da.wu_cov(ii) = NaN;
ii = da.wu_csp<wu_LoLim | da.wu_csp>wu_HiLim; da.wu_csp(ii) = NaN;
ii = da.wv_cov<wv_LoLim | da.wv_cov>wv_HiLim; da.wv_cov(ii) = NaN;
ii = da.wv_csp<wv_LoLim | da.wv_csp>wv_HiLim; da.wv_csp(ii) = NaN;
ii = da.wTs_cov<wT_LoLim | da.wTs_cov>wT_HiLim; da.wTs_cov(ii) = NaN;
ii = da.wTs_csp<wT_LoLim | da.wTs_csp>wT_HiLim; da.wTs_csp(ii) = NaN;
ii = da.wT_cov<wT_LoLim | da.wT_cov>wT_HiLim; da.wT_cov(ii) = NaN;
ii = da.wT_csp<wT_LoLim | da.wT_csp>wT_HiLim; da.wT_csp(ii) = NaN;

% apply general cu2,cw2,ct2 and cq2 outlier limits
ii = find(da.Cu2a<cu_LoLim |da.Cu2a>cu_HiLim); da.Cu2a(ii) = NaN;
ii = find(da.Cu2b<cu_LoLim |da.Cu2b>cu_HiLim); da.Cu2b(ii) = NaN;
ii = find(da.Cw2a<cw_LoLim |da.Cw2a>cw_HiLim); da.Cw2a(ii) = NaN;
ii = find(da.Cw2b<cw_LoLim |da.Cw2b>cw_HiLim); da.Cw2b(ii) = NaN;
ii = find(da.Ct2a<ct_LoLim |da.Ct2a>ct_HiLim); da.Ct2a(ii) = NaN;
ii = find(da.Ct2b<ct_LoLim |da.Ct2b>ct_HiLim); da.Ct2b(ii) = NaN;
ii = find(da.Cq2a<cq_LoLim |da.Cq2a>cq_HiLim); da.Cq2a(ii) = NaN;
ii = find(da.Cq2b<cq_LoLim |da.Cq2b>cq_HiLim); da.Cq2b(ii) = NaN;

% apply filter for Ts noise and variance limits
ii = (da.Ts_noise>Tnoise_lim) | (da.tvar>Tson_var_lim);
da.Ct2a(ii) = NaN; da.Ct2b(ii) = NaN;
da.wTs_cov(ii) = NaN; da.wTs_csp(ii) = NaN;
da.wT_cov(ii) = NaN; da.wT_csp(ii) = NaN;

disp('DONE DA_FIXIT')

