% fixes applied after evalflux.m and cat_mats.m, before run_motcorr.m and da_red.m
% EJT 2021

clear all;
close all;
setup_cruise; % set default conversion and calibration factors, and cruise specific names

% system specific path defs
sysType = computer;
username=char(java.lang.System.getProperty('user.name'));
if strncmp(sysType,'MACI64',7)     % set Mac paths
    data_drive = '/Users/eliz/DATA/';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs');
end

% matlab script path
restoredefaultpath
cd(fullfile(path_prog,'flux'));
addpath(genpath(fullfile(path_prog,'flux')));
% addpath('/Users/eliz/Documents/MATLAB/');
rehash toolboxcache;

% make directory for plots of fixes
warning ('off','MATLAB:MKDIR:DirectoryExists');
path_fix_plots = fullfile(data_drive,cruise,ship,'flux','Raw_Images','fixes');
mkdir(path_fix_plots);

%% load concatenated cruise data
in_version = 'v1'; % v1 for Jan 2022 tests
out_version = 'v2'; % v2 for Jan 2022 tests

% v0 met sea flux nav data
indir_1 = [data_drive cruise '/' ship '/flux/Processed/v0_1min/'];
indir_10 = [data_drive cruise '/' ship '/flux/Processed/v0_10min/'];
outdir = [data_drive cruise '/' ship '/flux/Processed/' out_version '/'];
mkdir(outdir);

load([indir_1 cruise '_1min_data_' in_version '.mat']); %% structure b1
load([indir_10 cruise '_10min_data_' in_version '.mat']); %% structure b10


%%% v1 files to be saved in this program
m_outfile_1 = [outdir cruise '_1min_nav_met_sea_flux_' out_version '.mat'];
m_outfile_10 = [outdir cruise '_10min_nav_met_sea_flux_' out_version '.mat'];

%%% adcp data from read_adcp.m in flux code directory 
adir = '/Users/eliz/DATA/PISTON_2019/Sally_Ride/ship/adcp/';
matname = 'piston-WH300_RV-Sally-Ride_20190904_R0.mat';
load([adir matname]);

%%% load otter data for QC
o1 = load([data_drive cruise '/Otter/piston-SurfOtter_RV-Sally-Ride_20190907_R0.mat']);
o2 = load([data_drive cruise '/Otter/piston-SurfOtter_RV-Sally-Ride_20190913_R0.mat']);
o3 = load([data_drive cruise '/Otter/piston-SurfOtter_RV-Sally-Ride_20190921_R0.mat']);

%%% ROSR data
load('/Users/eliz/DATA/PISTON_2019/Sally_Ride/SST/ROSR/all_data/skinT_PISTON_2019.mat');
% There isn't a great way to do this interp. The ROSR outputs a value at
% irregularly spaced time steps because its averaging period is 5 min, and
% on top of that it does not report when rain or excessive sea spray occurs
% automatically (i.e. when it detects enough liquid for voltage to drop on
% its primitive optical sensor). In practice, the rain sensor is cleaned
% off every day briefly to make sure salt or spray are artificially
% obscuring the sensor. 
b1.Tskin_ir = interp1(sea.t, sea.SSST, b1.t);
b10.Tskin_ir = interp1(sea.t, sea.SSST, b10.t);

min10 = datenum(2009,12,25,0,10,0) - datenum(2009,12,25,0,0,0);
for i = 1:length(b1.t)
    dtime = abs(b1.t(i) - sea.t(closetime(sea.t, b1.t(i))));
    if dtime > min10
        b1.Tskin_ir(i) = nan;
    end
end
for i = 1:length(b10.t)
    dtime = abs(b10.t(i) - sea.t(closetime(sea.t, b10.t(i))));
    if dtime > min10
        b10.Tskin_ir(i) = nan;
    end
end

plot_rosr_interp = 0;
if plot_rosr_interp == 1
    figure;
    subplot(2,1,1); hold on;
    plot(b1.t, b1.Tskin_ir,'o');
    plot(sea.t, sea.SSST,'x');
    legend('interp to 1 min','orig');
    grid on;
    xlim([min(b1.t) max(b1.t)]);
    datetick('x','mm/dd','keeplimits');

    subplot(2,1,2); hold on;
    plot(b10.t, b10.Tskin_ir,'o');
    plot(sea.t, sea.SSST,'x');
    legend('interp to 10 min','orig');
    grid on;
    xlim([min(b1.t) max(b1.t)]);
    datetick('x','mm/dd','keeplimits');
    
    print(graphdevice,[path_fix_plots '/rosr_check_' cruise graphformat]);
end


%% add or mask data from prior / new structures if needed

if nanmean(b10.hs > 0) && nanmean(b10.hs_s) < 0
    disp('***need to flip ship or PSL turb fluxes');
%     flips = {'hs_s';'hl_s';'hrain_s'};
%     for k = 1:length(flips)
%         eval(['b10.' flips{k} ' = -b10.' flips{k} ';']);
%         eval(['b1.' flips{k} ' = -b1.' flips{k} ';']);
end
    
    
% if isfield(b10,'q10') ==1
% %   %%% name changes and var changes unless the new evalflux program
% %   hasn't been run yet
%     b1.qa10 = b1.q10;
%     b1.qa10N = b1.q10N;
%     b1.qa10_s = b1.q10_s;
%     b1.qa10N_s = b1.q10N_s;
%     b1.qa2 = b1.q2;
%     b1.qa2_s = b1.q2_s;
%     b10.qa10 = b10.q10;
%     b10.qa10N = b10.q10N;
%     b10.qa10_s = b10.q10_s;
%     b10.qa10N_s = b10.q10N_s;
%     b10.qa2 = b10.q2;
%     b10.qa2_s = b10.q2_s;
% 
%     fields_to_remove = {'q10';'q10N';'q10_s';'q10N_s';'q2';'q2_s'};
%     b1 = rmfield(b1,fields_to_remove);
%     b10 = rmfield(b10,fields_to_remove);
% %    
% end

%% if you need to test output from evalflux.m... 
test_plot = 0;
if test_plot == 1
    figure;
    subplot(1,3,1);
    plot(b10.t, b10.hrain, b10.t, b10.hrain_s);
    title('hrain');
    datetick('x','mm/dd'); axis tight; grid on;
    subplot(1,3,2);
    plot(b10.t, b10.hl, b10.t, b10.hl_s);
    title('hl');
    datetick('x','mm/dd'); axis tight; grid on;
    legend('PSL','ship','location','best');
    subplot(1,3,3);
    plot(b10.t, b10.hs, b10.t, b10.hs_s);
    title('hs');
    datetick('x','mm/dd'); axis tight; grid on;
end

%% adjustments to temps... this is how it used to be done but has been phased out. 
% Even if used, it was only called in evalflux.
% Currently, corrections get done instead in fix_met_sea, not here.
% but potentially, this could be used if you knew ahead of time what the
% corrections should be. We usually don't and have to do tests to figure it
% out. So we assume here that no corrections are needed initially.

% td1_adj = 0;        % PIR dome 1
% tc1_adj = 0;        % PIR case 1
% td2_adj = 0;        % PIR dome 2
% tc2_adj = 0;        % PIR case 2
% td_ship_adj = 0;    % PIR dome ship
% tc_ship_adj = 0;    % PIR case ship
% tsea_adj = 0;       % sea snake
% ta_adj = 0;         % air temp
% 
% adj = [tsea_adj,ta_adj,td1_adj,tc1_adj,td2_adj,tc2_adj];
% ship_adj = [td_ship_adj,tc_ship_adj];

%% bad sea snake data... was left on deck and heated while ship within the EEZ
% wh_no_snk_1 =find(b1.t < 249.21);
% b1.Tsnk(wh_no_snk_1) = nan;
% wh_no_snk_10 =find(b10.t < 249.21);
% b10.Tsnk(wh_no_snk_10) = nan;

%% EEZ: NaN the data while in EEZ
EEZ = find(b1.t < datenum(2019,9,6,5,0,0) | b1.t >= datenum(2019,9,25,12,0,0));
notEEZ = find(b1.t >= datenum(2019,9,6,5,0,0) & b1.t < datenum(2019,9,25,12,0,0));
EEZ_10 = find(b10.t < datenum(2019,9,6,5,0,0) | b10.t >= datenum(2019,9,25,12,0,0));
notEEZ_10 = find(b10.t >= datenum(2019,9,6,5,0,0) & b10.t < datenum(2019,9,25,12,0,0));
EEZ_ID = zeros(length(b1.t),1);
EEZ_ID_10 = zeros(length(b10.t),1);
EEZ_ID(EEZ) = 1; % 1 = in EEZ... bad data
EEZ_ID_10(EEZ_10) = 1; % 1 = in EEZ... bad data


fmet = fields(b1);
field1 = zeros(1,length(fmet));

%%% if you want to just make missing data nan
% for j = 1:length(fmet)
%     x1 = b1.(fmet{j});
%     x10 = b10.(fmet{j});
%     if length(x1) > 100 
%         if strcmp(fmet(j),'t') ~= 1 && strcmp(fmet(j),'jd') ~= 1 && strcmp(fmet(j),'year') ~= 1 ...
%                 && strcmp(fmet(j),'month') ~= 1 && strcmp(fmet(j),'day') ~= 1 && strcmp(fmet(j),'hour') ~= 1 ...
%                 && strcmp(fmet(j),'minute') ~= 1 && strcmp(fmet(j),'second') ~= 1 ...
%                 && strcmp(fmet(j),'lat') ~= 1 && strcmp(fmet(j),'lon') ~= 1
% % %             disp(['applying EEZ for ' fmet{j}]);
%             x1(EEZ) = nan;
%             x10(EEZ_10) = nan;
%             b1.(fmet{j}) = x1;
%             b10.(fmet{j}) = x10; 
%          end
%     else
%         % track this as a 1-member array
%         field1(j) = 1;
%     end
% end

%%% if you want to remove bad data altogether
for j = 1:length(fmet)
    x1 = b1.(fmet{j});
    x10 = b10.(fmet{j});
%     if length(x1) > 100  %%% not needed because all fields are greater than 100; no readme fields
%         if strcmp(fmet(j),'t') ~= 1 && strcmp(fmet(j),'jd') ~= 1 && strcmp(fmet(j),'year') ~= 1 ...
%                 && strcmp(fmet(j),'month') ~= 1 && strcmp(fmet(j),'day') ~= 1 && strcmp(fmet(j),'hour') ~= 1 ...
%                 && strcmp(fmet(j),'minute') ~= 1 && strcmp(fmet(j),'second') ~= 1 ...
%                 && strcmp(fmet(j),'lat') ~= 1 && strcmp(fmet(j),'lon') ~= 1
% %             disp(['applying EEZ for ' fmet{j}]);
            x1(EEZ) = [];
            x10(EEZ_10) = [];
            b1.(fmet{j}) = x1;
            b10.(fmet{j}) = x10; 
%          end
%     else
%         % track this as a 1-member array
%         field1(j) = 1;
%     end
end

%% jd_10bin for averaging
jd_10bin = [b10.jd; b10.jd(end) + datenum(0,0,0,0,10,0)];
jd_1bin = [b1.jd; b1.jd(end) + datenum(0,0,0,0,1,0)];

%% speed through water... check to make sure it's in m/s not knots
% convert kt to m/s
% b1.spd_log_u_s = b1.spd_log_u_s/1.944;
% b1.spd_log_v_s = b1.spd_log_v_s/1.944;
% b1.spd_log_s = sqrt(b1.spd_log_u_s.^2 + b1.spd_log_v_s.^2);
% 
% b10.spd_log_u_s = b10.spd_log_u_s/1.944;
% b10.spd_log_v_s = b10.spd_log_v_s/1.944;
% b10.spd_log_s = sqrt(b10.spd_log_u_s.^2 + b10.spd_log_v_s.^2);

plot_spd_log = 0;
if plot_spd_log == 1
    figure('position',[1,1, 1200, 400]);
    subplot(1,3,1);hold on;
    plot(b1.t, b1.spd_log_u_s, b1.t, b1.spd_log_v_s); 
    ylim([-1 6]);
    grid on;
    ylabel('m s^{-1}');
    legend('spd log u','spd log v');
    datetick('x','mm/dd');
    title(ptitle)

    subplot(1,3,2);hold on;
    plot(b1.t, b1.spd_log_s, b1.t, b1.sog_s); 
    ylim([-1 6]);
    grid on;
    ylabel('m s^{-1}');
    legend('spd log','sog');
    datetick('x','mm/dd');

    subplot(1,3,3);hold on;
    plot(b1.spd_log_s, b1.sog_s,'.'); 
    ylim([-1 6]); xlim([-1 6]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    grid on;
    xlabel('spd log, m s^{-1}');
    ylabel('sog, m s^{-1}');
    axis square;

    print(graphdevice,[path_fix_plots '/spd_log_' cruise graphformat]);
end

%% determine best nav data... fill gaps to have continuous time series
% for PISTON 2019, PSL and ship data are nearly equivalent. Just need to
% fill a few gaps. Either with fill_ga

disp(['missing 1 min sog = ' sprintf('%i',length(find(isnan(b1.sog) == 1)))]);
disp(['missing 1 min sog_s = ' sprintf('%i',length(find(isnan(b1.sog_s) == 1)))]);
disp(['missing 1 min hed = ' sprintf('%i',length(find(isnan(b1.hed) == 1)))]);
disp(['missing 1 min hed_s = ' sprintf('%i',length(find(isnan(b1.hed_s) == 1)))]);
disp(['missing 1 min cog = ' sprintf('%i',length(find(isnan(b1.cog) == 1)))]);
disp(['missing 1 min cog_s = ' sprintf('%i',length(find(isnan(b1.cog_s) == 1)))]);
disp(['missing 1 min lat = ' sprintf('%i',length(find(isnan(b1.lat) == 1)))]);
disp(['missing 1 min lat_s = ' sprintf('%i',length(find(isnan(b1.lat_s) == 1)))]);
disp(['missing 1 min lon = ' sprintf('%i',length(find(isnan(b1.lon) == 1)))]);
disp(['missing 1 min lon_s = ' sprintf('%i',length(find(isnan(b1.lon_s) == 1)))]);

plot_nav = 0;
if plot_nav == 1
figure;
subplot(2,3,1); hold on;
plot(b10.sog, b10.sog_s,'x'); axis square; xlim([0 10]);
    plot(xlim, xlim,'--y'); title('SOG'); grid on;
    axis square; xlabel('NOAA'); ylabel('ship');
subplot(2,3,2); hold on;
plot(b10.cog, b10.cog_s,'x'); axis square; xlim([0 360]);
    plot(xlim, xlim,'--y'); title('COG'); grid on;
    axis square; xlabel('NOAA'); ylabel('ship');
subplot(2,3,3); hold on;
plot(b10.hed, b10.hed_s,'x'); axis square; xlim([0 360]);
    plot(xlim, xlim,'--y'); title('HEADING'); grid on;
    axis square; xlabel('NOAA'); ylabel('ship');
subplot(2,3,4); hold on;
plot(b10.lat, b10.lat_s,'x'); axis square; xlim([min(xlim)-1 max(xlim)+1]);...
    ylim(xlim); plot(xlim, xlim,'--y'); title('LAT'); grid on;
    axis square; xlabel('NOAA'); ylabel('ship');
subplot(2,3,5); hold on;
plot(b10.lon, b10.lon_s,'x'); axis square; xlim([min(xlim)-1 max(xlim)+1]); ...
    ylim(xlim); plot(xlim, xlim,'--y'); title('LON'); grid on;
    axis square; xlabel('NOAA'); ylabel('ship');

print(graphdevice,[path_fix_plots '/nav_' cruise graphformat]);

end

%% snake and TSG corrections; loading otter data
%%% decide whether to correct snake to the ROSR or the otter
%%% and whether to correct the TSG to the snake or Otter or ROSR
cor_type_snk = 'o';  % correct to otter
cor_type_tsg = 's';  % correct to snake... meaning correct to otter :)

% cor_type = 'o';

%   *** otter 1:
%                   T: [21×7321 double]
%                   z: [21×1 double]
%                  sn: {21×1 cell}
%                   S: [5×7321 double]
%                  zS: [5×1 double]
%                   p: [5×7321 double]
%                  zp: [5×1 double]
%                   t: [1×7321 double]
%      S_uncalibrated: [5×7321 double]
%     S_calibrate_idx: [1×7321 logical]
%     T_calibrate_idx: [1×7321 logical]
%      T_uncalibrated: [21×7321 double]
%            T_spline: [21×7321 double]
%                dTdz: [21×7321 double]
%            S_spline: [21×7321 double]
%                dSdz: [21×7321 double]
%         pden_spline: [21×7321 double]
%                  N2: [21×7321 double]
% readme: 'Processed by Ken Hughes's script create_otter_ctd_matrix_piston19 on 20-Sep-2019 23:19:09.
%           ... see below for more readme repeated

%   *** otter 2:
%                   T: [21×5831 double]
%                   z: [21×1 double]
%                  sn: {21×1 cell}
%                   S: [5×5831 double]
%                  zS: [5×1 double]
%                   p: [5×5831 double]
%                  zp: [5×1 double]
%                   t: [1×5831 double]
%      S_uncalibrated: [5×5831 double]
%     S_calibrate_idx: [1×5831 logical]
%     T_calibrate_idx: [1×5831 logical]
%      T_uncalibrated: [21×5831 double]
%            T_spline: [21×5831 double]
%                dTdz: [21×5831 double]
%            S_spline: [21×5831 double]
%                dSdz: [21×5831 double]
%         pden_spline: [21×5831 double]
%                  N2: [21×5831 double]
% readme: 'Processed by Ken Hughes's script create_otter_ctd_matrix_piston19 on 21-Sep-2019 00:35:42.
% Temperature, salinity, and pressure data from all of the sensors on the surf otter, excluding fast thermistors (GusTs, TPods).
% Temperature and salinity are 'calibrated' by minimizing the offset to their respective means at night, 
% when the surface layer is assumed isothermal/isohaline unless there is an obvious, true gradient in T or S.
% Uncalibrated T and S are included for reference. As are the time indices at which data were included in the calibration.
% Sensors included some or all of RBR Solo T (sn=10???? or 076???), Idronaut CTDs (sn=021663?), 
% RBR Concertos (sn=06????), and Seabird SBE37s (sn=371660?). T and S gradients calculated from spline fits.
% Splines fit with Matlab's csaps function with a smoothing value (p) dependent on the size of the signal'

% find otter data closest to TSG # 2 (port aft) level... between 2.3 - 3.4 m. Nominally, we're calling it 3 m.
wh_otter_Tsea_depth_1 = find(o1.otter.z > 2.3 & o1.otter.z < 3.4);
wh_otter_Tsea_depth_2 = find(o2.otter.z > 2.3 & o2.otter.z < 3.4);
wh_otter_Tsea_depth_3 = find(o3.otter.z > 2.3 & o3.otter.z < 3.4);

wh_otter_Ssea_depth_1 = find(o1.otter.zS > 2.3 & o1.otter.zS < 3.4);
wh_otter_Ssea_depth_2 = find(o2.otter.zS > 2.3 & o2.otter.zS < 3.4);
wh_otter_Ssea_depth_3 = find(o3.otter.zS > 2.3 & o3.otter.zS < 3.4);

% do a mean along depth indicies
o1_mean_at_Tsea = nanmean(o1.otter.T(wh_otter_Tsea_depth_1,:), 1);
o1_mean_at_Ssea = o1.otter.S(wh_otter_Ssea_depth_1,:);
o2_mean_at_Tsea = nanmean(o2.otter.T(wh_otter_Tsea_depth_2,:), 1);
o2_mean_at_Ssea = o2.otter.S(wh_otter_Ssea_depth_2,:);
o3_mean_at_Tsea = nanmean(o3.otter.T(wh_otter_Tsea_depth_3,:), 1);
o3_mean_at_Ssea = o3.otter.S(wh_otter_Ssea_depth_3,:);

o1_top_at_Tsea = o1.otter.T(wh_otter_Tsea_depth_1(1),:);
o2_top_at_Tsea = o2.otter.T(wh_otter_Tsea_depth_2(1),:);
o3_top_at_Tsea = o3.otter.T(wh_otter_Tsea_depth_3(1),:);

o1_base_at_Tsea = o1.otter.T(wh_otter_Tsea_depth_1(end),:);
o2_base_at_Tsea = o2.otter.T(wh_otter_Tsea_depth_2(end),:);
o3_base_at_Tsea = o3.otter.T(wh_otter_Tsea_depth_3(end),:);

o_mean_at_Tsea = [o1_mean_at_Tsea o2_mean_at_Tsea o3_mean_at_Tsea];
o_mean_at_Ssea = [o1_mean_at_Ssea o2_mean_at_Ssea o3_mean_at_Ssea];

o_top_at_Tsea = [o1_top_at_Tsea o2_top_at_Tsea o3_top_at_Tsea];

o_base_at_Tsea = [o1_base_at_Tsea o2_base_at_Tsea o3_base_at_Tsea];

%%% grab data
oT_at_snake = [o1.otter.T(1,:) o2.otter.T(1,:) o3.otter.T(1,:)];  % from 0.10 m below vehicle == 0.6 m total
oS_at_snake = [o1.otter.S(1,:) o2.otter.S(1,:) o3.otter.S(1,:)];  % from 0.35 m below vehicle == 0.85 m total

ot = [o1.otter.t o2.otter.t o3.otter.t];
no_otter_1 = find(b1.t < o1.otter.t(1));
no_otter_2 = find(b1.t > o1.otter.t(end) & b1.t < o2.otter.t(1));
no_otter_3 = find(b1.t > o2.otter.t(end) & b1.t < o3.otter.t(1));
no_otter_4 = find(b1.t > o3.otter.t(end));
no_otter = [no_otter_1; no_otter_2; no_otter_3; no_otter_4];
ID_no_otter = zeros(length(b1.t), 1);
ID_no_otter(no_otter) = 1;
ID_otter = double(~ID_no_otter);
wh_otter = find(ID_otter == 1);

oSsea_mean = interp1(ot, o_mean_at_Ssea, b1.t);
oTsea_mean = interp1(ot, o_mean_at_Tsea, b1.t);
oTsea_top = interp1(ot, o_top_at_Tsea, b1.t);
oTsea_base = interp1(ot, o_base_at_Tsea, b1.t);
oTsnake = interp1(ot, oT_at_snake, b1.t);
oSsnake = interp1(ot, oS_at_snake, b1.t);

oSsea_mean(no_otter) = nan;
oTsea_mean(no_otter) = nan;
oTsea_top(no_otter) = nan;
oTsea_base(no_otter) = nan;
oTsnake(no_otter) = nan;
oSsnake(no_otter) = nan;

%%% add a few otter variables to structures
b1.ID_otter = ID_otter;
b1.ID_no_otter = ID_no_otter;
b1.S_0p85 = oSsnake;
b10.S_0p85 = interval_avg_var(b1.jd, b1.S_0p85, jd_10bin);

wh_otter_10 = find(isfinite(b10.S_0p85) == 1);
b10.ID_otter = zeros(length(b10.t), 1);
b10.ID_otter(wh_otter_10) = 1;
b10.ID_no_otter = double(~b10.ID_otter);
wh_no_otter_10 = find(b10.ID_no_otter == 1);

%% infrared radiation
sig_sb = 5.67e-8; % Stefan Boltzmann constant
C2K = 273.15; % K
c3_1 = 3.63; % coefficient k1 for PIR1 from June 2017 calibration... not the correct one but it's all we've got so far
c3_2 = 3.37; % coefficient k2 for PIR2 from June 2017 calibration... not the correct one but it's all we've got so far
k2_corr = 0; % start off as zero (not sure what this is)

% IDEAS: try different c3 values
% 18 W = how much it matters if case and dome are 1 deg different. 
% an extreme error is on the order of 10 W. 

% note: PIR2 domeT was all bad because KZ unand PIR1 caseT aren't trustworthy. 
%%% this is perhaps a solution. Thermopile outputs look good. Then
%%% PIR2 had best case T, but PIR2 had best case T. For these equations,
%%% measured therm is positive and should be negative. So we add a negative
%%% sign to the first term of this equation.
pir1 = -b1.lw_therm_1 + sig_sb*(b1.lw_T_case_2+C2K).^4 - c3_1*sig_sb*((b1.lw_T_dome_1+C2K).^4-(b1.lw_T_case_2+C2K).^4);
% pir2 = -b1.lw_therm_2 + sig_sb*(b1.lw2_T_case+C2K).^4 - k1*sig_sb*((b1.lw1_T_dome+C2K).^4-(b1.lw2_T_case+C2K).^4);
pir2 = -b1.lw_therm_2 + sig_sb*(b1.lw_T_case_2+C2K).^4;
% pir2_corr = -b1.lw_therm_2 + sig_sb*(b1.lw_T_case_2+C2K).^4 - k2_corr*b1.sw1; 

%%% ... correcting as a function of solar forcing, which is causing the dome to heat up. 
% k2_corr is a new calibration coefficient for adjusting the IR for solar
% heating when the unit doesn't measure dome T
plot_sw_ir = 0;
if plot_sw_ir == 1
    figure; plot(b1.sw_dn_1, pir2-pir1,'.');
    ylabel('corrected pir2 - pir1');
    xlabel('sw dn 1');
    print(graphdevice,[path_fix_plots '/IR_Solar_' cruise graphformat]);
end

plot_ir = 0;
if plot_ir == 1
    figure;
    plot(b1.t, b1.lw_dn_1, b1.t, pir1, b1.t, pir2);
    legend('orig PIR1','new PIR1','new PIR2');%, 'ship');
    datetick;
    ylabel('W m^{-2}');
    grid on;
    title(['original and corrected IR_{down} ' cruise]);
    print(graphdevice,[path_fix_plots '/IR_' cruise graphformat]);
end

% save new 
b1.lw_dn_1 = pir1;
b1.lw_dn_2 = pir2;
b10.lw_dn_1 = interval_avg_var(b1.jd, b1.lw_dn_1, jd_10bin);
b10.lw_dn_2 = interval_avg_var(b1.jd, b1.lw_dn_2, jd_10bin);

% check PIR temp offsets from night air temps
wh_night = find(b1.hour < 21 | b1.hour > 9); %% for PISTON
tc1_adj2 = nanmean(b1.Ta(wh_night)-b1.lw_T_case_1(wh_night)); disp(['tc1_adj = ',sprintf('%4.2f',tc1_adj2)]);
td1_adj2 = nanmean(b1.Ta(wh_night)-b1.lw_T_dome_1(wh_night)); disp(['td1_adj = ',sprintf('%4.2f',td1_adj2)]);
tc2_adj2 = nanmean(b1.Ta(wh_night)-b1.lw_T_case_2(wh_night)); disp(['tc2_adj = ',sprintf('%4.2f',tc2_adj2)]);
td2_adj2 = nanmean(b1.Ta(wh_night)-b1.lw_T_dome_2(wh_night)); disp(['td2_adj = ',sprintf('%4.2f',td2_adj2)]);
tc_sh_adj = nanmean(b1.Ta(wh_night)-b1.lw_T_dome_s(wh_night)); disp(['tc_sh_adj = ',sprintf('%4.2f',tc_sh_adj)]);
td_sh_adj = nanmean(b1.Ta(wh_night)-b1.lw_T_dome_s(wh_night)); disp(['td_sh_adj = ',sprintf('%4.2f',td_sh_adj)]);

%% bad and moderate wind flags
wh_mod_wind = find(abs(b1.rdir) > 90);
wh_bad_wind = find(abs(b1.rdir) > 120);

b1.ID_mod_wind = b1.t*0;
b1.ID_bad_wind = b1.t*0;
b1.ID_mod_wind(wh_mod_wind) = 1;
b1.ID_bad_wind(wh_bad_wind) = 1;

wh_mod_wind_10 = find(abs(b10.rdir) > 90);
wh_bad_wind_10 = find(abs(b10.rdir) > 120);

b10.ID_mod_wind = b10.t*0;
b10.ID_bad_wind = b10.t*0;
b10.ID_mod_wind(wh_mod_wind_10) = 1;
b10.ID_bad_wind(wh_bad_wind_10) = 1;


%% wind speed
% readjust wind if needed

%% snake T corrections
       
%%% Tsnake T correction... originally estimated at 0.12 C in the field
%%% Check otter and rosr at night again to be sure
dT_snk_est = 0.12;

wh_pm_r = find(b1.hour >= 16 & b1.hour <= 21 & b1.Tskin_ir > 28.72);
wh_pm_o = find(b1.hour >= 16 & b1.hour <= 21 & ID_otter == 1);
dT_snk_otter = b1.Tsnk(wh_pm_o) - oTsnake(wh_pm_o);
dT_sst_rosr = b1.Tskin(wh_pm_r) - b1.Tskin_ir(wh_pm_r);

offset_snk_otter = prctilex(dT_snk_otter, 50);
offset_sst_rosr = prctilex(dT_sst_rosr, 50);

Tsnk_cor_o = b1.Tsnk-offset_snk_otter;
Tsnk_cor_r = b1.Tsnk-offset_sst_rosr;

% look at resulting SST changes... but don't apply them because we'll
% recalculate SST with COARE later anyway.
plot_sst = 0;
if plot_sst == 1
    figure('position',[1,1, 1200, 800]);
    subplot(2,3,1); hold on;
    plot(b1.Tsnk(wh_pm_o), oTsnake(wh_pm_o),'.');
    ylabel('60 cm T_{otter}');
    xlabel('50 cm T_{snake}');
    ylim([28.5 29.3]);xlim([28.5 29.3]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('snake offset');
    axis square;
    grid on;
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.05*diff(ylim), ['med = ' sprintf('%3.3f',offset_snk_otter)],...
        'fontsize',16);
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.1*diff(ylim), ...
        ['mean = ' sprintf('%3.3f',nanmean(dT_snk_otter))],'fontsize',16);


    subplot(2,3,2); hold on;
    plot(Tsnk_cor_o(wh_pm_o), oTsnake(wh_pm_o),'.');
    ylabel('60 cm T_{otter}');
    xlabel('50 cm T_{snake}');
    ylim([28.5 29.3]);xlim([28.5 29.3]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('Corrected snake: otter');
    axis square;
    grid on;
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.05*diff(ylim), ...
        ['rmse = ' sprintf('%3.3f',rmse(Tsnk_cor_o(wh_pm_o), oTsnake(wh_pm_o)))],'fontsize',16);

    subplot(2,3,3); hold on;
    plot(Tsnk_cor_r(wh_pm_o), oTsnake(wh_pm_o),'.');
    ylabel('60 cm T_{otter}');
    xlabel('50 cm T_{snake}');
    ylim([28.5 29.3]);xlim([28.5 29.3]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('Corrected snake: rosr');
    axis square;
    grid on;
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.05*diff(ylim), ...
        ['rmse = ' sprintf('%3.3f',rmse(Tsnk_cor_r(wh_pm_o), oTsnake(wh_pm_o)))],'fontsize',16);


    subplot(2,3,4); hold on;
    plot(b1.Tskin(wh_pm_r), b1.Tskin_ir(wh_pm_r),'.');
    ylabel('0 cm SST');
    xlabel('0 cm T_{rosr}');
    ylim([28.5 29.3]);xlim([28.5 29.3]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('SST offset');
    axis square;
    grid on;
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.05*diff(ylim), ['med = ' sprintf('%3.3f',offset_sst_rosr)],...
        'fontsize',16);
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.1*diff(ylim), ...
        ['mean = ' sprintf('%3.3f',nanmean(dT_sst_rosr))],'fontsize',16);

    subplot(2,3,5); hold on;
    plot(b1.Tskin(wh_pm_r)-offset_sst_rosr, b1.Tskin_ir(wh_pm_r),'.');
    ylabel('0 cm SST');
    xlabel('0 cm T_{rosr}');
    ylim([28.5 29.3]);xlim([28.5 29.3]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('Corrected SST: rosr');
    axis square;
    grid on;
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.05*diff(ylim), ...
        ['rmse = ' sprintf('%3.3f',rmse(b1.Tskin(wh_pm_r)-offset_sst_rosr, b1.Tskin_ir(wh_pm_r)))],'fontsize',16);

    subplot(2,3,6); hold on;
    plot(b1.Tskin(wh_pm_r)-offset_snk_otter, b1.Tskin_ir(wh_pm_r),'.');
    ylabel('0 cm SST');
    xlabel('0 cm T_{rosr}');
    ylim([28.5 29.3]);xlim([28.5 29.3]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('Corrected SST: otter');
    axis square;
    grid on;
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.05*diff(ylim), ...
        ['rmse = ' sprintf('%3.3f',rmse(b1.Tskin(wh_pm_r)-offset_snk_otter, b1.Tskin_ir(wh_pm_r)))],'fontsize',16);
    print(graphdevice,[path_fix_plots '/snake_T_' cruise graphformat]);
end

plot_snake_before_after = 0;
if plot_snake_before_after == 1
    figure;
    plot(b1.t, b1.Tsnk,'o', b1.t, Tsnk_cor_o,'-');
    legend('before','after');
    grid on;
    axis tight;
    datetick('x','mm/dd','keeplimits');
    title([ptitle ' snake T']);
    print(graphdevice,[path_fix_plots '/snake_T_orig_corr_' cruise graphformat]);
end

%% choose corrected snake value and save
if cor_type_snk == 'r'
    % correct Tsnk to ROSR
    b1.Tsnk = Tsnk_cor_r; 
%     b1.Tskin = b1.Tskin - offset_sst_rosr; %%% this is going to get redone
%     anyway
elseif cor_type_snk == 'o'
    % correct Tsnk to otter
    b1.Tsnk = Tsnk_cor_o; 
%     b1.Tskin = b1.Tskin - offset_snk_otter; %%% this is going to get
%     redone anyway
end

% redo qsnk from Tsnk... uses same code as COARE
b1.qsnk = qsea_sp(b1.Tsnk, b1.psealevel, b1.Ssea2_s);
b10.Tsnk = interp1(b1.t, b1.Tsnk, b10.t);
b10.qsnk = qsea_sp(b10.Tsnk, b10.psealevel, b10.Ssea2_s);

%% QC TSG data
% switch in TSGs lasted from 9/11 23:39 to 9/12 4 UTC
switchstart = datenum(2019,9,11,23,39,0);
switchend = datenum(2019,9,12,4,0,0);
bad_TSG_switch = find(b1.t > switchstart & b1.t < switchend); % primary flow was switched from bow to stern intake
more_bad_TSG_a = find(b1.Fsea2_s < 2.2 & b1.t < datenum(2019,9,12,0,0,0)); % too low before switch
more_bad_TSG_b = find(b1.Fsea2_s < 1.8 & b1.t > datenum(2019,9,15,0,0,0)); % too low after switch
more_bad_TSG_c = find(b1.t >= datenum(2019,9,11,4,0,0) & b1.t < datenum(2019,9,11,5,15,0)); % flow rate too high
more_bad_TSG_d = find(b1.Fsea2_s > 2.85 & b1.t > datenum(2019,9,15,0,0,0)); % too high after switch
more_bad_TSG_e = find(b1.t > datenum(2019,9,7,4,0,0) & b1.t < datenum(2019,9,7,6,45,0)); % too high before switch
bad_TSG = [bad_TSG_switch; more_bad_TSG_a; more_bad_TSG_b;  more_bad_TSG_c; more_bad_TSG_d; more_bad_TSG_e; ...
    more_bad_TSG_a-2; more_bad_TSG_a-1;more_bad_TSG_a+2;more_bad_TSG_a+1;...
    more_bad_TSG_b-2; more_bad_TSG_b-1;more_bad_TSG_b+2;more_bad_TSG_b+1;...
    more_bad_TSG_d-2; more_bad_TSG_d-1;more_bad_TSG_d+2;more_bad_TSG_d+1;...
    (more_bad_TSG_d(end)+3:more_bad_TSG_d(end)+22)'];

b1.ID_bad_TSG = zeros(length(b1.t),1);
b1.ID_bad_TSG(bad_TSG) = 1;
b10.ID_bad_TSG = ceil(interval_avg_var(b1.jd, b1.ID_bad_TSG, jd_10bin));

plot_t_s_orig = 0;
if plot_t_s_orig == 1
    
    figure('position',[1,1,1400,800]);
    subplot(3,1,1);hold on;
    plot(b1.t, b1.Ssea1_s,'-',b1.t, b1.Ssea2_s, '-'); 
    plot(b1.t(bad_TSG), b1.Ssea1_s(bad_TSG),'xr',b1.t(bad_TSG), b1.Ssea2_s(bad_TSG), 'xr'); 
    grid on;
    xlim([min(b1.t) max(b1.t)]);
    datetick('x','mm/dd','keeplimits');
    ylabel('psu');
    title([ptitle ' TSG S']);
    legend('TSG 1','TSG 2','location','south');

    subplot(3,1,2);hold on;
    plot(b1.t, b1.Tsea1_s,'-',b1.t, b1.Tsea2_s, '-'); 
    plot(b1.t(bad_TSG), b1.Tsea1_s(bad_TSG),'xr',b1.t(bad_TSG), b1.Tsea2_s(bad_TSG), 'xr'); 
    grid on;
    xlim([min(b1.t) max(b1.t)]);
    datetick('x','mm/dd','keeplimits');
    ylabel('^oC');
    title([ptitle ' TSG T']);
    legend('TSG 1','TSG 2','location','south');

    subplot(3,1,3);
    yyaxis left;
    hold on;
    plot(b1.t, b1.ID_tsg_s,'.k'); 
    plot(b1.t(bad_TSG), b1.ID_tsg_s(bad_TSG), 'xr');
    grid on;
    ylim([0 2.5]);
    xlim([min(b1.t) max(b1.t)]);
    text(min(b1.t)+1, 0.2,'bad','fontweight','bold','fontsize',16);
    text(min(b1.t)+1, 1.2,'TSG 1 main','fontweight','bold','fontsize',16);
    text(min(b1.t)+1, 2.2,'TSG 2 main','fontweight','bold','fontsize',16);
    datetick('x','mm/dd','keeplimits');
    title([ptitle ' TSG ID']);
    yyaxis right;
    hold on;
    plot(b1.t, b1.Fsea1_s,'-','color',rgb('purple'));
    plot(b1.t, b1.Fsea2_s,'-','color',rgb('coral')); 
    plot(b1.t(bad_TSG), b1.Fsea1_s(bad_TSG),'xr');
    plot(b1.t(bad_TSG), b1.Fsea2_s(bad_TSG),'xr'); 
    ylabel('flow rate');
    legend('ID','bad','TSG 1','TSG 2','bad','bad');

    print(graphdevice,[path_fix_plots '/TSG_TS_orig_' cruise graphformat]);

end

b1.Tsea1_s(bad_TSG) = nan;
b1.Tsea2_s(bad_TSG) = nan;
b1.Ssea1_s(bad_TSG) = nan;
b1.Ssea2_s(bad_TSG) = nan;
b1.ID_tsg_s(bad_TSG) = 0;

%% TSG T correction
Tsea_best = b1.Tsea2_s;  % temp
Ssea_best = b1.Ssea2_s;  % salinity

%%% determine lat night - early morning offset between TSG and Tsnk
wh_pre_pm = find(b1.t <= switchstart & b1.hour >= 16 & b1.hour <= 21);
wh_pre_pm_o = find(b1.t <= switchstart & b1.hour >= 16 & b1.hour <= 21 & ID_otter == 1);
wh_post_pm_A = find(b1.t >= switchend & b1.t <= datenum(2019,9,18,0,0,0) & b1.hour >= 16 & b1.hour <= 21);
wh_post_pm_B = find(b1.t > datenum(2019,9,18,0,0,0) & b1.hour >= 16 & b1.hour <= 21);
wh_post_pm_A_o = find(b1.t >= switchend & b1.t <= datenum(2019,9,18,0,0,0) & b1.hour >= 16 & b1.hour <= 21 & ID_otter == 1);
wh_post_pm_B_o = find(b1.t > datenum(2019,9,18,0,0,0) & b1.hour >= 16 & b1.hour <= 21 & ID_otter == 1);

wh_fix_pre = find(b1.t <= switchstart);
wh_fix_post_A = find(b1.t >= switchend & b1.t <= datenum(2019,9,18,0,0,0));
wh_fix_post_B = find(b1.t > datenum(2019,9,18,0,0,0));

wh_test = [wh_pre_pm; wh_post_pm_A; wh_post_pm_B];
wh_test_o = [wh_pre_pm_o; wh_post_pm_A_o; wh_post_pm_B_o];


%%% try this with corrected b1.Tsnk from above and with o_mean_at_Tsea
dT_pre_snk      = Tsea_best(wh_pre_pm) - b1.Tsnk(wh_pre_pm);
dT_post_A_snk   = Tsea_best(wh_post_pm_A) - b1.Tsnk(wh_post_pm_A);
dT_post_B_snk   = Tsea_best(wh_post_pm_B) - b1.Tsnk(wh_post_pm_B);

dT_pre_o        = Tsea_best(wh_pre_pm_o) - oTsea_mean(wh_pre_pm_o);
dT_post_A_o     = Tsea_best(wh_post_pm_A_o) - oTsea_mean(wh_post_pm_A_o);
dT_post_B_o     = Tsea_best(wh_post_pm_B_o) - oTsea_mean(wh_post_pm_B_o);

t_pre = b1.t(wh_pre_pm);
t_post = b1.t([wh_post_pm_A; wh_post_pm_B]);

offset_pre_snk = prctilex(dT_pre_snk, 50);
offset_post_A_snk = prctilex(dT_post_A_snk, 50);
offset_post_B_snk = prctilex(dT_post_B_snk, 50);

offset_pre_o_T = prctilex(dT_pre_o, 50);
offset_post_A_o_T = prctilex(dT_post_A_o, 50);
offset_post_B_o_T = prctilex(dT_post_B_o, 50);

% correct Tsea to Tsnk, which was corrected to otter
Tsea_cor_o = Tsea_best;
Tsea_cor_s = Tsea_best;

% pick the best offset to apply
% if cor_type_tsg == 'r'
    Tsea_cor_s(wh_fix_pre) = Tsea_best(wh_fix_pre) - offset_pre_snk;
    Tsea_cor_s(wh_fix_post_A) = Tsea_best(wh_fix_post_A) - offset_post_A_snk;
    Tsea_cor_s(wh_fix_post_B) = Tsea_best(wh_fix_post_B) - offset_post_B_snk;
% elseif cor_type_tsg == 'o'
    Tsea_cor_o(wh_fix_pre) = Tsea_best(wh_fix_pre) - offset_pre_o_T;
    Tsea_cor_o(wh_fix_post_A) = Tsea_best(wh_fix_post_A) - offset_post_A_o_T;
    Tsea_cor_o(wh_fix_post_B) = Tsea_best(wh_fix_post_B) - offset_post_B_o_T;
% end

plot_t = 0;
if plot_t == 1
    
    figure('position',[1,1, 1200, 800]);
    subplot(2,3,1); hold on;
    plot(Tsea_best(wh_pre_pm), b1.Tsnk(wh_pre_pm),'.');
    plot(Tsea_best(wh_post_pm_A), b1.Tsnk(wh_post_pm_A),'.');
    plot(Tsea_best(wh_post_pm_B), b1.Tsnk(wh_post_pm_B),'.');
    legend('pre','post A','post B','location','southeast');
    ylabel('50 cm T_{snake}');
    xlabel('3 m T_{TSG}');
    ylim([28.6 29.4]);xlim([28.6 29.4]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('TSG T offset');
    axis square;
    grid on;

    subplot(2,3,2); hold on;
    plot(Tsea_cor_s(wh_pre_pm), b1.Tsnk(wh_pre_pm),'.');
    plot(Tsea_cor_s(wh_post_pm_A), b1.Tsnk(wh_post_pm_A),'.');
    plot(Tsea_cor_s(wh_post_pm_B), b1.Tsnk(wh_post_pm_B),'.');
    legend('pre','post A','post B','location','southeast');
    ylabel('50 cm T_{snake}');
    xlabel('3 m T_{TSG}');
    ylim([28.6 29.4]);xlim([28.6 29.4]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('corrected TSG T: snake');
    axis square;
    grid on;
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.05*diff(ylim), ...
        ['rmse = ' sprintf('%3.3f',rmse(Tsea_cor_s(wh_test), b1.Tsnk(wh_test)))],'fontsize',16);

    subplot(2,3,3); hold on;
    plot(Tsea_cor_o(wh_pre_pm), b1.Tsnk(wh_pre_pm),'.');
    plot(Tsea_cor_o(wh_post_pm_A), b1.Tsnk(wh_post_pm_A),'.');
    plot(Tsea_cor_o(wh_post_pm_B), b1.Tsnk(wh_post_pm_B),'.');
    legend('pre','post A','post B','location','southeast');
    ylabel('50 cm T_{snake}');
    xlabel('3 m T_{TSG}');
    ylim([28.6 29.4]);xlim([28.6 29.4]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('corrected TSG T: otter');
    axis square;
    grid on;
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.05*diff(ylim), ...
        ['rmse = ' sprintf('%3.3f',rmse(Tsea_cor_o(wh_test), b1.Tsnk(wh_test)))],'fontsize',16);


    subplot(2,3,4); hold on;
    plot(Tsea_best(wh_pre_pm_o), oTsea_mean(wh_pre_pm_o),'.');
    plot(Tsea_best(wh_post_pm_A_o), oTsea_mean(wh_post_pm_A_o),'.');
    plot(Tsea_best(wh_post_pm_B_o), oTsea_mean(wh_post_pm_B_o),'.');
    legend('pre','post A','post B','location','southeast');
    ylabel('3 m T_{otter}');
    xlabel('3 m T_{TSG}');
    ylim([28.6 29.4]);xlim([28.6 29.4]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('TSG T offset');
    axis square;
    grid on;

    subplot(2,3,5); hold on;
    plot(Tsea_cor_o(wh_pre_pm_o), oTsea_mean(wh_pre_pm_o),'.');
    plot(Tsea_cor_o(wh_post_pm_A_o), oTsea_mean(wh_post_pm_A_o),'.');
    plot(Tsea_cor_o(wh_post_pm_B_o), oTsea_mean(wh_post_pm_B_o),'.');
    legend('pre','post A','post B','location','southeast');
    ylabel('3 m T_{otter}');
    xlabel('3 m T_{TSG} corrected');
    ylim([28.6 29.4]);xlim([28.6 29.4]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('corrected TSG T: otter');
    axis square;
    grid on;
    print(graphdevice,[path_fix_plots '/TSG_T_' cruise graphformat]);
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.05*diff(ylim), ...
        ['rmse = ' sprintf('%3.3f',rmse(Tsea_cor_o(wh_test_o), oTsea_mean(wh_test_o)))],'fontsize',16);

    subplot(2,3,6); hold on;
    plot(Tsea_cor_s(wh_pre_pm_o), oTsea_mean(wh_pre_pm_o),'.');
    plot(Tsea_cor_s(wh_post_pm_A_o), oTsea_mean(wh_post_pm_A_o),'.');
    plot(Tsea_cor_s(wh_post_pm_B_o), oTsea_mean(wh_post_pm_B_o),'.');
    legend('pre','post A','post B','location','southeast');
    ylabel('3 m T_{otter}');
    xlabel('3 m T_{TSG} corrected');
    ylim([28.6 29.4]);xlim([28.6 29.4]);
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('corrected TSG T: snake');
    axis square;
    grid on;
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.05*diff(ylim), ...
        ['rmse = ' sprintf('%3.3f',rmse(Tsea_cor_s(wh_test_o), oTsea_mean(wh_test_o)))],'fontsize',16);
    print(graphdevice,[path_fix_plots '/TSG_T_' cruise graphformat]);

end

%% fill missing holes in Tsea and then save
Tsea_cor_fill = fillmissing(Tsea_cor_s, 'linear');

%%% The snake is corrected to the sparse otter data, and the snake data is
%%% then used to correct the TSG since its available for longer record

b1.Tsea_s = Tsea_cor_fill;
b10.Tsea_s = interval_avg_var(b1.jd, b1.Tsea_s, jd_10bin);

%% seawater S
%%% determine late night - early morning offset between TSG and Tsnk
dS_pre_o        = Ssea_best(wh_pre_pm_o) - oSsea_mean(wh_pre_pm_o);
dS_post_A_o     = Ssea_best(wh_post_pm_A_o) - oSsea_mean(wh_post_pm_A_o);
dS_post_B_o     = Ssea_best(wh_post_pm_B_o) - oSsea_mean(wh_post_pm_B_o);

offset_pre_o_S = prctilex(dS_pre_o, 50);
offset_post_A_o_S = prctilex(dS_post_A_o, 50);
offset_post_B_o_S = prctilex(dS_post_B_o, 50);

Ssea_cor = Ssea_best;
Ssea_cor(wh_fix_pre) = Ssea_best(wh_fix_pre) - offset_pre_o_S;
Ssea_cor(wh_fix_post_A) = Ssea_best(wh_fix_post_A) - offset_post_A_o_S;
Ssea_cor(wh_fix_post_B) = Ssea_best(wh_fix_post_B) - offset_post_B_o_S;

Ssea_cor_fill = fillmissing(Ssea_cor, 'linear');
b1.Ssea_s = Ssea_cor_fill;
b10.Ssea_s = interval_avg_var(b1.jd, b1.Ssea_s, jd_10bin);

plot_S = 0;
if plot_S == 1
    figure('position',[1,1, 800, 400]);
    subplot(1,2,1); hold on;
    plot(Ssea_best(wh_test_o), oSsea_mean(wh_test_o),'.');
    ylabel('3 m S_{otter}');
    xlabel('3 m S_{TSG}');
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    grid on;
    title('TSG S offset');
    axis square;

    subplot(1,2,2); hold on;
    plot(Ssea_cor(wh_test_o), oSsea_mean(wh_test_o),'.');
    ylabel('3 m S_{otter}');
    xlabel('3 m S_{TSG} corrected');
    plot(ylim, ylim, '--','color',rgb('skyblue'));
    title('Corrected TSG S');
    grid on;
    axis square;
    text(min(xlim)+0.1*diff(xlim), max(ylim)-0.05*diff(ylim), ...
        ['rmse = ' sprintf('%3.3f',rmse(Ssea_cor(wh_test_o), oSsea_mean(wh_test_o)))],'fontsize',16);
    print(graphdevice,[path_fix_plots '/TSG_S_' cruise graphformat]);
end


plot_t_s_after = 0;
if plot_t_s_after == 1
    figure('position',[1,1,1400,800]);
    subplot(3,1,1);
    plot(b1.t, b1.Ssea1_s,'x',b1.t, b1.Ssea2_s, 's',b1.t, Ssea_cor,'.'); 
    grid on;
    xlim([min(b1.t) max(b1.t)]);
    datetick('x','mm/dd','keeplimits');
    ylabel('psu');
    title([ptitle ' TSG S']);
    legend('TSG 1','TSG 2','TSG 2 QCed','location','south');

    subplot(3,1,2);
    plot(b1.t, b1.Tsea1_s,'x',b1.t, b1.Tsea2_s, 's',b1.t, Tsea_cor_s,'.'); 
    grid on;
    xlim([min(b1.t) max(b1.t)]);
    datetick('x','mm/dd','keeplimits');
    ylabel('^oC');
    title([ptitle ' TSG T']);
    legend('TSG 1','TSG 2','TSG 2 QCed','location','south');

    subplot(3,1,3);
    yyaxis left;
    plot(b1.t, b1.ID_tsg_s,'ok'); 
    grid on;
    ylim([0 2.5]);
    xlim([min(b1.t) max(b1.t)]);
    text(min(b1.t)+1, 0.2,'bad','fontweight','bold','fontsize',16);
    text(min(b1.t)+1, 1.2,'TSG 1 main','fontweight','bold','fontsize',16);
    text(min(b1.t)+1, 2.2,'TSG 2 main','fontweight','bold','fontsize',16);
    datetick('x','mm/dd','keeplimits');
    title([ptitle ' TSG ID']);
    yyaxis right;
    hold on;
    plot(b1.t, b1.Fsea1_s,'-','color',rgb('purple'));
    plot(b1.t, b1.Fsea2_s,'-','color',rgb('coral')); 
    ylabel('flow rate');
    legend('ID','TSG 1','TSG 2');

    print(graphdevice,[path_fix_plots '/TSG_TS_orig_corr_' cruise graphformat]);

    figure;
    subplot(2,1,1);
    plot(b1.t, b1.Tsea_s,'o',b1.t, Tsea_cor_s,'.');
    legend('filled','unfilled','location','best');
    grid on;
    axis tight;
    datetick('x','mm/dd','keeplimits');

    subplot(2,1,2);
    plot(b1.t, b1.Ssea_s,'o',b1.t, Ssea_cor,'.');
    legend('filled','unfilled','location','best');
    grid on;
    axis tight;
    datetick('x','mm/dd','keeplimits');

    print(graphdevice,[path_fix_plots '/TSG_TS_corr_fill_' cruise graphformat]);

end

%% seawater sigma
b1.psea_s = sw_pres(ones(length(b1.t),1)*3,b1.lat); % dbar
b1.rhosea_s = sw_dens(b1.Ssea_s,b1.Tsea_s,b1.psea_s); % kg / m3
b10.psea_s = sw_pres(ones(length(b10.t),1)*3,b10.lat); % dbar
b10.rhosea_s = sw_dens(b10.Ssea_s,b10.Tsea_s,b10.psea_s); % kg / m3

%% rain rate 
% a prior threshold of 0.03 mm/hr was used for ORG... why? Don't think that's necessary. 
% prate(prate<0.03) = 0;     % threshold of 0.03 mm/hr

%%% reset nan to zero... should already be taken care of in evalflux.m
bad_prate = find(isnan(b1.prate) == 1);
disp(['resetting NaN prate 1-min of length = ' sprintf('%i',length(bad_prate)) ' to zero']);
b1.prate(bad_prate) = 0; % reset nan to 0


%% adjust licor variables if updates to qa or Ta are performed, or if bad data are found

%% q licor plot

% For PISTON 2019, the licor stopped washing itself after this point and so
% its data should just not be used. There is a drifting offset in qa past this point.
bad_licor_day = datenum(2019,9,7,22,0,0);

%%% from da_red... do not change coefficients
AGC_lim     = 60;   % max licor AGC
CO2_LoLim   = 350;  % min CO2 ppm
CO2_HiLim   = 450;  % max CO2 ppm
sigCO2_lim  = 20;   % max sigma CO2
q_lic_LoLim = 0.1;  % min q licor
q_lic_HiLim = 26;   % max q licor
q_std_LoLim  = 0;    % min sigma q licor
q_std_HiLim  = 3;    % max sigma q licor

clear ii; 
ii = find(b10.licor_AGC>AGC_lim | b10.licor_qa<q_lic_LoLim | b10.licor_qa>q_lic_HiLim); 
%     usually we also look for.... find(b10.licor_qa_std>q_std_HiLim | b10.licor_qa_std<q_std_LoLim);
if strcmp(cruise,'PISTON_2019') == 1
   iii = find(b10.t > bad_licor_day);
   ii = unique([ii; iii]);
end

% %%% special --> save these updates too
b10.licor_qa(ii) = nan;
% save(b10_infile, 'b10');

t0 = min(b10.t);
tN = max(b10.t);
figure;
subplot(2,1,1); 
    plot(b10.t,b10.qa,'or',b10.t,b10.licor_qa,'xk'); grid; xlim([t0 tN]);
    ylabel('q, g kg^{-1}'); datetick('x','dd','keeplimits');
    title([ptitle,' specific humidity']);
    legend('Vaisala','Licor','location','best');
subplot(2,1,2);
    plot(b10.t,b10.licor_AGC,'.r',[t0,tN],[AGC_lim,AGC_lim],'g--',[t0,tN],[50,50],'g-');
    grid; ylabel('AGC'); datetick('x','dd','keeplimits');
    ylim([48,100]); xlim([t0 tN]);
    if strcmp(cruise,'PISTON_2019') == 1
       hold on;
       plot([bad_licor_day bad_licor_day], ylim,'--k');
       text(bad_licor_day-0.25, AGC_lim, 'bad after this point',...
           'rotation',90,'fontweight','bold','fontsize',16);
    end
    legend('Licor AGC','AGC Limit','AGC = 50','location','northoutside','orientation','horizontal');
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA PSL plot 27'},'FontSize',10,'FitBoxToText','off','LineStyle','none');
print(graphdevice,[path_fix_plots '/licor_qa_' cruise graphformat]);


%% fill missing holes in snake with adjusted ROSR or TSG if needed
% % would need to use something like the below code, except with
% corrections for cool skin and warm layer as applicable
% b1.Tsnk_best = b1.Tsnk;
% b10.Tsnk_best = b10.Tsnk;

% % choose what to fill the holes in the snake with ROSR - cool skin
% fill_snk = 'rosr';
% if strcmp(fill_snk,'rosr') == 1
%     wh_bad_snk_good_rosr = find(isnan(b1.Tsnk) == 1 & isfinite(SSST1) == 1);
%     Tsnk_best_rosr = SSST1 + b1.dT_skin;
%     b1.Tsnk_best(wh_bad_snk_good_rosr) = Tsnk_best_rosr(wh_bad_snk_good_rosr);
%     
%     wh_bad_snk_good_rosr = find(isnan(b10.Tsnk) == 1 & isfinite(b10.Tskin) == 1);
%     Tsnk_best_rosr = b10.Tskin + b10.dT_skin;
%     b10.Tsnk_best(wh_bad_snk_good_rosr) = Tsnk_best_rosr(wh_bad_snk_good_rosr);
%     
% elseif strcmp(fill_snk,'tsgs') == 1
%     wh_bad_snk_good_tsg = find(isnan(b1.Tsnk) == 1 & isfinite(b1.Ttsg_s) == 1);
%     b1.Tsnk_best(wh_bad_snk_good_tsg) = b1.Ttsg_s(wh_bad_snk_good_tsg);
%     
%     wh_bad_snk_good_tsg = find(isnan(b10.Tsnk) == 1 & isfinite(b10.Ttsg_s) == 1);
%     b10.Tsnk_best(wh_bad_snk_good_tsg) = b10.Ttsg_s(wh_bad_snk_good_tsg);
% end

%% best radiative fluxes and met data... for both b10 and b1
b1.sw_dn = (b1.sw_dn_1 + b1.sw_dn_2)/2; % if averaging solar sensors together
b1.lw_dn = (b1.lw_dn_1 + b1.lw_dn_2)/2; % if averaging IR sensors together

b10.sw_dn = (b10.sw_dn_1 + b10.sw_dn_2)/2; % if averaging solar sensors together
b10.lw_dn = (b10.lw_dn_1 + b10.lw_dn_2)/2; % if averaging IR sensors together 

%% redo clear sky values if changes are made to Ta or qa, or to tune it locally
%   function coefficients can be tuned for location and aerosol loading
% b10.lw_dn_clr = (.52+.13/60*abs(f10.lat)+(.082-.03/60*abs(f10.lat)).*...
%     sqrt(f10.qa)).*(5.67e-8*(f10.Ta+C2K).^4);

b10.sw_dn_clr = rs_clear(b10.jd,nanmean(b10.psealevel),b10.qa,b10.lat,b10.lon,k1,k2,oz);
b1.sw_dn_clr = rs_clear(b1.jd,nanmean(b1.psealevel),b1.qa,b1.lat,b1.lon,k1,k2,oz);

b10.lw_dn_clr = (0.52+0.13/60*abs(b10.lat)+(0.082-0.03/60.*abs(b10.lat)).*sqrt(b10.qa)).*(5.67e-8*(b10.Ta+273.15).^4);
b1.lw_dn_clr = (0.52+0.13/60*abs(b1.lat)+(0.082-0.03/60.*abs(b1.lat)).*sqrt(b1.qa)).*(5.67e-8*(b1.Ta+273.15).^4);

%% redo paccum for full cruise of data
temp = b1.prate;
temp(isnan(temp)) = 0;
b1.paccum = cumsum(temp)/60; % 60 min per hour

temp = b1.prate_s;
temp(isnan(temp)) = 0;
b1.paccum_s = cumsum(temp)/60; % 60 min per hour

temp = b10.prate;
temp(isnan(temp)) = 0;
b10.paccum = cumsum(temp)/6; % 6 x 10 min per hour

temp = b10.prate_s;
temp(isnan(temp)) = 0;
b10.paccum_s = cumsum(temp)/6; % 6 x 10 min per hour

plot_check_rain = 0;
if plot_check_rain == 1
    
figure;
subplot(2,2,1);
plot(b1.t, b1.paccum,'o', b10.t, b10.paccum,'x'); 
legend('1-min','10-min','location','northwest');
title('NOAA');
ymax_noaa = max(b1.paccum);
datetick('x','mm/dd'); grid on; ylabel('\Sigma Precip (mm)');
subplot(2,2,2);
plot(b1.t, b1.paccum_s,'o', b10.t, b10.paccum_s,'x'); 
ylim([0 ymax_noaa]);
legend('1-min','10-min','location','northwest');
title('ship');
datetick('x','mm/dd'); grid on; ylabel('\Sigma Precip (mm)');

subplot(2,2,3);
plot(b1.t, b1.prate,'o', b10.t, b10.prate,'x'); 
legend('1-min','10-min','location','northwest');
title('NOAA');
ymax_noaa = max(b1.prate);
datetick('x','mm/dd'); grid on; ylabel('P rate (mm/hr)');
subplot(2,2,4);
plot(b1.t, b1.prate_s,'o', b10.t, b10.prate_s,'x'); 
ylim([0 ymax_noaa]);
legend('1-min','10-min','location','northwest');
title('ship');
datetick('x','mm/dd'); grid on; ylabel('P rate (mm/hr)');

print(graphdevice,[path_fix_plots '/rain_' cruise graphformat]);

end

%% bad L value from ship data COARE output. Reset to nan.
bad_L1 = find(b1.L_s > 5e5);
bad_L10 = find(b10.L_s > 5e5);
b1.L_s(bad_L1) = nan;
b10.L_s(bad_L10) = nan;

%% missing PSL data.... no missing PSL data 
wh_no_Tsnk = find(isnan(b1.Tsnk) == 1 & isnan(b1.Tsea_s) == 0);
wh_no_Ta = find(isnan(b1.Ta) == 1 & isnan(b1.Ta_s) == 0);
wh_no_rh = find(isnan(b1.rh) == 1 & isnan(b1.rh_s) == 0);
wh_no_psealevel = find(isnan(b1.psealevel) == 1 & isnan(b1.psealevel_s) == 0);
wh_no_wspd = find(isnan(b1.wspd) == 1 & isnan(b1.wspd_s) == 0);
wh_no_wdir = find(isnan(b1.wdir) == 1 & isnan(b1.wdir_s) == 0);
wh_no_rspd = find(isnan(b1.rspd) == 1 & isnan(b1.rspd_s) == 0);
wh_no_rdir = find(isnan(b1.rdir) == 1 & isnan(b1.rdir_s) == 0);
wh_no_lat = find(isnan(b1.lat) == 1 & isnan(b1.lat_s) == 0);
wh_no_lon = find(isnan(b1.lon) == 1 & isnan(b1.lon_s) == 0);
wh_no_hed = find(isnan(b1.hed) == 1 & isnan(b1.hed_s) == 0);
wh_no_cog = find(isnan(b1.cog) == 1 & isnan(b1.cog_s) == 0);
wh_no_sog = find(isnan(b1.sog) == 1 & isnan(b1.sog_s) == 0);
 
% b1.Ta_best = b1.Ta;
% b1.Ta_best(wh_no_Ta) = b1.Ta_s;
% 
% b1.rh_best = b1.rh;
% b1.rh_best(wh_no_rh) = b1.rh_s;
% 
% b1.slp_best = b1.slp;
% b1.slp_best(wh_no_slp) = b1.slp_s;
% 
% b1.wspd_best = b1.wspd;
% b1.wspd_best(wh_no_wspd) = b1.wspd_s;
% 
% b1.rspd_best = b1.rspd;
% b1.rspd_best(wh_no_rspd) = b1.rspd_s;
% 
% b1.wdir_best = b1.wdir;
% b1.wdir_best(wh_no_wdir) = b1.wdir_s;
% 
% b1.rdir_best = b1.rdir;
% b1.rdir_best(wh_no_rdir) = b1.rdir_s;

%% currents
%%% true = relative to fixed earth coords (from gps)
%%% true_sfc = relative to ocean (from spd_log and/or currents)
%%% odec was a true ship log for speed, not used anymore. New ships have
%%% speed log though.

% wind direction toward, for reference
wdir_toward = b10.wdir+180;
wdir_toward(wdir_toward>360) = wdir_toward(wdir_toward> 360) - 360;

% relative wind direction in 0-360 format
rdir2 = b10.rdir;
rdir2(rdir2<0)=rdir2(rdir2<0)+360;

cu = interp1(a.t, a.usfc_fill, b1.t);
cv = interp1(a.t, a.vsfc_fill, b1.t);
cu_10 = interp1(a.t, a.usfc_fill, b10.t);
cv_10 = interp1(a.t, a.vsfc_fill, b10.t);

% compute surface current from odec, head, sog_ship, cog_ship
% current direction is 'toward' = oceanographic convention, or like the
% vector would appear

%%% check to make sure the signs and equations are correct? 

[cspd_adcp, cdir_adcp] = uv_to_sd_ocean(cu, cv);
[cspd_adcp_10, cdir_adcp_10] = uv_to_sd_ocean(cu_10, cv_10);

[cspd_spd_log,cdir_spd_log]         = current(b1.spd_log_s,b1.hed,b1.sog,b1.cog);
[cspd_spd_log_10,cdir_spd_log_10]   = current(b10.spd_log_s,b10.hed,b10.sog,b10.cog);

wh_bad_cspd_spd_log = find(cspd_spd_log > 5);
wh_bad_cspd_spd_log_10 = find(cspd_spd_log_10 > 5);
cspd_spd_log(wh_bad_cspd_spd_log) = nan;
cdir_spd_log(wh_bad_cspd_spd_log) = nan;
cspd_spd_log_10(wh_bad_cspd_spd_log_10) = nan;
cdir_spd_log_10(wh_bad_cspd_spd_log_10) = nan;

% use spd_log as main version, but fill gaps from adcp. If further gaps are
% there, could just use regular wind speed and direction... ? 
cspd = cspd_spd_log; 
wh_bad_spd_log = find(isnan(cspd_spd_log) == 1);
cspd(wh_bad_spd_log) = cspd_adcp(wh_bad_spd_log);

cspd_10 = cspd_spd_log_10; 
wh_bad_spd_log_10 = find(isnan(cspd_spd_log_10) == 1);
cspd_10(wh_bad_spd_log_10) = cspd_adcp_10(wh_bad_spd_log_10);

cdir = cdir_spd_log; 
wh_bad_spd_log_dir = find(isnan(cdir_spd_log) == 1);
cdir(wh_bad_spd_log_dir) = cdir_adcp(wh_bad_spd_log_dir);

cdir_10 = cdir_spd_log_10; 
wh_bad_spd_log_dir_10 = find(isnan(cdir_spd_log_10) == 1);
cdir_10(wh_bad_spd_log_dir_10) = cdir_adcp_10(wh_bad_spd_log_dir_10);

disp(['missing 1 min cspd = ' sprintf('%i',length(find(isnan(cspd) == 1)))]);
disp(['missing 1 min cdir = ' sprintf('%i',length(find(isnan(cdir) == 1)))]);
disp(['missing 10 min cspd = ' sprintf('%i',length(find(isnan(cspd_10) == 1)))]);
disp(['missing 10 min cdir = ' sprintf('%i',length(find(isnan(cdir_10) == 1)))]);

plot_current = 0;
if plot_current == 1
    figure;
    subplot(1,2,1);
    plot(b1.t, cspd_adcp,'o',b1.t, cspd_spd_log,'x');
    datetick; grid on;
    legend('ADCP','SPEED LOG');
    ylabel('current speed, m/s');
    title(ptitle)

    subplot(1,2,2);
    plot(b1.t, cdir_adcp,'o',b1.t, cdir_spd_log,'x');
    datetick; grid on;
    legend('ADCP','SPEED LOG');
    ylabel('current dir, ^o');
    print(graphdevice,[path_fix_plots '/currents_' cruise graphformat]);

    figure;
    subplot(1,2,1);
    plot(cspd_adcp,cspd_spd_log,'.k');
    hold on; plot([0 1.4],[0 1.4],'--y');
    grid on; axis square;
    xlabel('ADCP');
    ylabel('SPEED LOG');
    title('current speed, m/s');

    subplot(1,2,2);
    plot(cdir_adcp,cdir_spd_log,'.k');
    hold on; plot([0 360],[0 360],'--y');
    grid on; axis square;
    xlabel('ADCP');
    ylabel('SPEED LOG');
    title('current dir, ^o');
    print(graphdevice,[path_fix_plots '/currents_scatter' cruise graphformat]);
    
    figure;
    subplot(1,2,1);
    plot(b10.t, cdir_10, 'o',b10.t, rdir2, 'x', b10.t, wdir_toward,'.');
    legend('cdir toward','rdir','wdir toward');
    grid on;
    datetick('x','mm/dd');
    
    subplot(1,2,2);
    plot(cdir_10,rdir2, 'x', cdir_10, wdir_toward,'o');
    legend('rdir','wdir toward');
    xlabel('cdir ^o');
    grid on;
    
    print(graphdevice,[path_fix_plots '/currents_dir_' cruise graphformat]);
    
end


% sub function current wasnt working because of dimensions so I changed the function by: * -> .*
% this is how we get cspd and cdir from topmost level of ADCP
% compute wind w/respect to moving surface
[wspd_sfc,wdir_sfc] = wind_current_adjust(b1.wspd,b1.wdir,cspd,cdir);
[wspd_sfc_spd_log,wdir_sfc_spd_log] = wind_current_adjust(b1.wspd,b1.wdir,cspd_spd_log,cdir_spd_log);
[wspd_sfc_adcp,wdir_sfc_adcp] = wind_current_adjust(b1.wspd,b1.wdir,cspd_adcp,cdir_adcp);

[wspd_sfc_10,wdir_sfc_10] = wind_current_adjust(b10.wspd,b10.wdir,cspd_10,cdir_10);
[wspd_sfc_spd_log_10,wdir_sfc_spd_log_10] = wind_current_adjust(b10.wspd,b10.wdir,cspd_spd_log_10,cdir_spd_log_10);
[wspd_sfc_adcp_10,wdir_sfc_adcp_10] = wind_current_adjust(b10.wspd,b10.wdir,cspd_adcp_10,cdir_adcp_10);

plot_wind_sfc = 0;
if plot_wind_sfc == 1
    
    figure;
    subplot(1,2,1);
    plot(b1.t, wspd_sfc_adcp,'o',b1.t, wspd_sfc_spd_log,'x',...
        b1.t, wspd_sfc,'o',b1.t, b1.wspd,'.');
    grid on; datetick;
    legend('ADCP adjusted','SPEED LOG adjusted','hybrid adjusted','original','location','southoutside');
    ylabel('wind speed, m/s');
    title(ptitle);

    subplot(1,2,2);
    plot(b1.t, wdir_sfc_adcp,'o',b1.t, wdir_sfc_spd_log,'x',...
        b1.t, wdir_sfc,'o',b1.t, b1.wdir,'.');
    grid on; datetick;
    legend('ADCP adjusted','SPEED LOG adjusted','hybrid adjusted','original','location','southoutside');
    ylabel('wind dir, ^o');
    print(graphdevice,[path_fix_plots '/wind_sfc_' cruise graphformat]);

    figure;
    subplot(1,2,1);
    plot(b1.wspd, wspd_sfc_adcp, 'o', b1.wspd, wspd_sfc_spd_log,'x', ... 
        b1.wspd, wspd_sfc, '.');
    hold on;
    plot([0 18],[0 18],'--k');
    grid on;
    xlabel('original wind speed m/s');
    ylabel('adjusted');
    legend('ADCP adjusted','SPEED LOG adjusted','hybrid adjusted','location','northwest');
    title(ptitle);
    subplot(1,2,2);
    plot(b1.wdir, wdir_sfc_adcp, 'o', b1.wdir, wdir_sfc_spd_log,'x', ...
                b1.wdir, wdir_sfc, '.');
    hold on;
    plot([0 360],[0 360],'--k');
    grid on;
    xlabel('original wind dir, ^o');
    ylabel('adjusted');
    legend('ADCP adjusted','SPEED LOG adjusted','hybrid adjusted','location','northoutside','orientation','horizontal');
    print(graphdevice,[path_fix_plots '/wind_sfc_scatter_' cruise graphformat]);


    figure;
    plot(b1.t, cspd_adcp,'o', b1.t, cspd_spd_log,'x',...
        b1.t, cspd,'s', b1.t, b1.wspd/10,'.');
    datetick;
    grid on;
    title(ptitle);
    legend('ADCP current speed','SPEED LOG current speed','hybrid current speed','wind speed / 10','location','southoutside');
    print(graphdevice,[path_fix_plots '/wind_current_' cruise graphformat]);

end


% add adcp data to structure
b1.cspd_adcp = cspd_adcp;
b10.cspd_adcp = cspd_adcp_10;
b1.cdir_adcp = cdir_adcp;
b10.cdir_adcp = cdir_adcp_10;
b1.cspd_u_adcp = cu;
b1.cspd_v_adcp = cv;
b10.cspd_u_adcp = cu_10;
b10.cspd_v_adcp = cv_10;

% and/or add spd log current to structure
b1.cspd_spdlog = cspd_spd_log;
b10.cspd_spdlog = cspd_spd_log_10;
b1.cdir_spdlog = cdir_spd_log;
b10.cdir_spdlog = cdir_spd_log_10;

b1.cspd = cspd;
b10.cspd = cspd_10;
b1.cdir = cdir;
b10.cdir = cdir_10;

%%% We think cspd/cdir from spd_log is superior for final dataset because
%%% it is more representative of the surface waters than the water sensed
%%% by the ADCP at 10-14 m (10.1, 12.1, 14.1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%% wspd sfc: use spd log version, with gaps
%%%%%%%%%%%%%%%%%%%%%%%%%%%% filled by ADCP.
b1.wspd_sfc = wspd_sfc; 
b1.wdir_sfc = wdir_sfc;
b10.wspd_sfc = wspd_sfc_10; 
b10.wdir_sfc = wdir_sfc_10;

%% redo   and SST and save

disp(['missing 1 min wspd_sfc = ' sprintf('%i',length(find(isnan(b1.wspd_sfc) == 1)))]);
disp(['missing 1 min wdir_sfc = ' sprintf('%i',length(find(isnan(b1.wdir_sfc) == 1)))]);
disp(['missing 1 min Ta = ' sprintf('%i',length(find(isnan(b1.Ta) == 1)))]);
disp(['missing 1 min rh = ' sprintf('%i',length(find(isnan(b1.rh) == 1)))]);
disp(['missing 1 min psealevel = ' sprintf('%i',length(find(isnan(b1.psealevel) == 1)))]);
disp(['missing 1 min Tsnk = ' sprintf('%i',length(find(isnan(b1.Tsnk) == 1)))]);
disp(['missing 1 min sw_dn = ' sprintf('%i',length(find(isnan(b1.sw_dn) == 1)))]);
disp(['missing 1 min lw_dn = ' sprintf('%i',length(find(isnan(b1.lw_dn) == 1)))]);
disp(['missing 1 min prate = ' sprintf('%i',length(find(isnan(b1.prate) == 1)))]);
disp(['missing 1 min Ssea_s = ' sprintf('%i',length(find(isnan(b1.Ssea_s) == 1)))]);
disp(['missing 1 min wspd_sfc = ' sprintf('%i',length(find(isnan(b1.wspd_sfc) == 1)))]);

disp(['missing 10 min wdir_sfc = ' sprintf('%i',length(find(isnan(b10.wdir_sfc) == 1)))]);
disp(['missing 10 min Ta = ' sprintf('%i',length(find(isnan(b10.Ta) == 1)))]);
disp(['missing 10 min rh = ' sprintf('%i',length(find(isnan(b10.rh) == 1)))]);
disp(['missing 10 min psealevel = ' sprintf('%i',length(find(isnan(b10.psealevel) == 1)))]);
disp(['missing 10 min Tsnk = ' sprintf('%i',length(find(isnan(b10.Tsnk) == 1)))]);
disp(['missing 10 min sw_dn = ' sprintf('%i',length(find(isnan(b10.sw_dn) == 1)))]);
disp(['missing 10 min lw_dn = ' sprintf('%i',length(find(isnan(b10.lw_dn) == 1)))]);
disp(['missing 10 min prate = ' sprintf('%i',length(find(isnan(b10.prate) == 1)))]);
disp(['missing 10 min Ssea_s = ' sprintf('%i',length(find(isnan(b10.Ssea_s) == 1)))]);

%% albedo
% check albedo function for proper year day, units/sign of lon, and lat

% vector version of function
[alb,T_sw,solarmax_sw,psi_sw] = albedo_vector_test(b10.sw_dn,b10.jd,b10.lon,b10.lat,'E');

% single value version of function
alb_k = nan(length(b10.t),1);
T_sw_k = nan(length(b10.t),1);
solarmax_sw_k = nan(length(b10.t),1);
psi_sw_k = nan(length(b10.t),1);

% single value version of function
for k = 1:length(b10.t)
    [alb_k(k),T_sw_k(k),solarmax_sw_k(k),psi_sw_k(k)] = albedo_test(b10.sw_dn(k),b10.jd(k),b10.lon(k),b10.lat(k),'E');
end

plot_albedo = 0;
if plot_albedo == 1
    figure; 
    subplot(2,2,1); hold on;
    plot(b10.t, alb); title('albedo');
    plot(b10.t, alb_k,'--');
    plot(b10.t, ones(length(b10.t), 1)*0.055);
    legend('vector function','single value function','constant value 0.0550','position',[0.5,0.5,0.01,0.01]);
    datetick('x','mm/dd'); axis tight;
    subplot(2,2,2); hold on;
    plot(b10.t, T_sw); title('T');
    plot(b10.t, T_sw_k,'--');
    datetick('x','mm/dd'); axis tight;
    subplot(2,2,3); hold on;
    plot(b10.t, solarmax_sw); title('solarmax');
    plot(b10.t, solarmax_sw_k,'--');
    datetick('x','mm/dd'); axis tight;
    subplot(2,2,4); hold on;
    plot(b10.t, psi_sw); title('psi');
    plot(b10.t, psi_sw_k,'--');
    datetick('x','mm/dd'); axis tight;
    print(graphdevice,[path_fix_plots '/albedo_' cruise graphformat]);
end

%% run coare

% testing air-sea temp difference and theta expressions
grav = grv(b10.lat);
lapse=grav/cpa;

P_tq=(b10.psealevel - (0.125*zt));   % P at tq measurement height (mb)

% Ta converted to K from sensor height, p in mb, theta in K
theta = (b10.Ta+C2K).*(1000./P_tq).^(Rgas/cpa); % potential temp at instrument height in K
TadjK = b10.Ta+C2K + lapse*zt; % adjusted temp at instrument height in K
Tadj = b10.Ta + lapse*zt; % adjusted temp at instrument height in C
theta_sfc = (b10.Tskin+C2K).*(1000./b10.psealevel).^(Rgas/cpa); % potential temp surface in K
TadjK_sfc = b10.Tskin+C2K; % adjusted temp surface in K
Tadj_sfc = b10.Tskin; % adjusted temp surface in C

% different versions of differences, of the same type
dtheta = theta_sfc - theta;
dTadjK = TadjK_sfc - TadjK;
dTadj = Tadj_sfc - Tadj;

% just for fun, mix and match the different types
dtheta_TadjK = theta_sfc - TadjK;
dTadjK_theta = TadjK_sfc - theta;

%%% test plot
plot_theta = 0;
if plot_theta == 1
    figure;
    subplot(2,2,1);
    plot(b10.t, theta, b10.t, TadjK,'--');
    legend('\theta zt','TadjK zt');
    ylabel('K');
    datetick('x','mm/dd');
    grid on;
    
    subplot(2,2,2);
    plot(b10.t, theta_sfc, b10.t, TadjK_sfc,'--');
    legend('\theta SFC','TadjK SFC');
    ylabel('K');
    datetick('x','mm/dd');
    grid on;   
    
    subplot(2,2,3);
    plot(b10.t, dtheta, b10.t, dTadjK,'--',b10.t, dTadj,'.');
    legend('theta sfc - theta zt','TadjK sfc - TadjK zt', 'Tadj sfc - Tadj zt');
    ylabel('K');
    datetick('x','mm/dd');
    grid on;     
    
    subplot(2,2,4);
    plot(b10.t, dtheta_TadjK, b10.t, dTadjK_theta,'--');
    legend('theta sfc - TadjK zt','TadjK sfc - theta zt');
    ylabel('K');
    datetick('x','mm/dd');
    grid on;         
    
    print(graphdevice,[path_fix_plots '/theta_' cruise graphformat]);
end

ref = 2;
%%% v 3.6 for PSL data
B=coare36vnWarm_et(b10.jd, b10.wspd_sfc, zu, b10.Ta, zt, b10.rh, zq, b10.psealevel, b10.Tsnk, b10.sw_dn, b10.lw_dn, b10.lat, b10.lon, 600, b10.prate, zsnk, b10.Ssea_s, nan, nan, ref, ref, ref);
%%% v 3.6 for ship data
Bs=coare36vnWarm_et(b10.jd, b10.wspd_s, zu_ship, b10.Ta_s, zt_ship, b10.rh_s, zq_ship, b10.psealevel_s, b10.Tsea_s, b10.sw_dn_s, b10.lw_dn_s, b10.lat, b10.lon, 600, b10.prate, zsea2_ship, b10.Ssea_s, nan, nan, ref, ref, ref);
%%% v 3.5 for PSL data
C=coare35vnWarm_et(b10.jd, b10.wspd_sfc, zu, b10.Ta, zt, b10.rh, zq, b10.psealevel, b10.Tsnk, b10.sw_dn, b10.lw_dn, b10.lat, b10.lon, 600, b10.prate, zsnk, ref);

cfields = {'usr';'tau';'hs';'hl';'hb';'hb_son';'hlwebb';'tsr';'qsr';'zo';'zot';'zoq';...
    'Cd';'Ch';'Ce';'L';'zeta';'dT_skin';'dq_skin';'dz_skin';'U2';'Ta2';'qa2';...
    'rh2';'U2N';'Ta2N';'qa2N';'lw_net';'sw_net';'Le';'rhoa';'UN';'U10';'U10N';'CdN10';'ChN10';'CeN10';...
    'hrain';'qs';'erate';'Ta10';'Ta10N';'qa10';'qa10N';'rh10';'pa10';'rhoa10';'gust';'wc_frac';'Edis';...
    'dT_warm';'dz_warm';'dT_warm_to_skin';'du_warm'};

% % testing buoyancy flux v3.6
% cfields = {'usr';'tau';'hs';'hl';'hb';'hb1';'hb_son';'hb_son1';'hlwebb';'tsr';'qsr';'zo';'zot';'zoq';...
%     'Cd';'Ch';'Ce';'L';'zeta';'dT_skin';'dq_skin';'dz_skin';'U2';'T2';'qa2';...
%     'rh2';'U2N';'T2N';'qa2N';'lw_net';'sw_net';'Le';'rhoa';'UN';'U10';'U10N';'CdN10';'ChN10';'CeN10';...
%     'hrain';'qs';'erate';'T10';'T10N';'qa10';'qa10N';'rh10';'pa10';'rhoa10';'gust';'wc_frac';'Edis';...
%     'dT_warm';'dz_warm';'dT_warm_to_skin';'du_warm'};

cfields_35 = {'usr';'tau';'hs';'hl';'hb';'hb_son';'hlwebb';'tsr';'qsr';'zo';'zot';'zoq';...
    'Cd';'Ch';'Ce';'L';'zeta';'dT_skin';'dq_skin';'dz_skin';'U2';'Ta2';'qa2';...
    'rh2';'U2N';'lw_net';'sw_net';'Le';'rhoa';'UN';'U10';'U10N';'CdN10';'ChN10';'CeN10';...
    'hrain';'qs';'erate';'Ta10';'Ta10N';'qa10';'qa10N';'rh10';'pa10';'rhoa10';'gust';...
    'dT_warm';'dz_warm';'dT_warm_to_skin'};

missing_35 = {'du_warm';'wc_frac';'Edis';'Ta2N';'qa2N'};

for i = 1:length(cfields)
    % just name PSL fields for now, add to array later
    eval([cfields{i} '_10 = B(:,i);']);
    eval([cfields{i} ' = interp1(b10.t,' cfields{i} '_10, b1.t);' ]);
    
    % add ship fields to array
    eval(['b10.' cfields{i} '_s = Bs(:,i);']);
    eval(['b1.' cfields{i} '_s = interp1(b10.t, b10.' cfields{i} '_s, b1.t);' ]);
end

% do any coare checks of 3.6 that are needed
plot_new_buoy_flux = 0;
if plot_new_buoy_flux == 1
    figure;
    
    subplot(1,2,1); hold on;
    plot(b10.t, hb1_10, b10.t, hb_10,'--');
    legend('original','new');
    title('buoyancy flux');
    ylabel('W m^{-2}');
    datetick('x','mm/dd');
    grid on;
    
    subplot(1,2,2); hold on;
    plot(b10.t, hb_son1_10, b10.t, hb_son_10,'--');
    legend('original','new');
    title('sonic buoyancy flux');
    ylabel('W m^{-2}');
    datetick('x','mm/dd');
    grid on;
    
    print(graphdevice,[path_fix_plots '/newflux_' cruise graphformat]);
end


% interpolate from 10 min to 1 min
for i = 1:length(cfields_35)
    eval([cfields_35{i} '_10_35 = C(:,i);']);
    eval([cfields_35{i} '_35 = interp1(b10.t,' cfields_35{i} '_10_35, b1.t);' ]);
end

nanfill = nan(length(tau),1);
nanfill_10 = nan(length(tau_10),1);

for i = 1:length(missing_35)
    eval([missing_35{i} '_10_35 = nanfill_10;']);
    eval([missing_35{i} '_35 = nanfill;']);
end

%%% option: sign changes of COARE fields... the bulk sensible, latent, and rain
%%% heat fluxes are defined positive by COARE. 
% flips = {'hs';'hl';'hrain'};
% for k = 1:length(flips)
%     eval([flips{k} ' = - ' flips{k} ';']);
%     eval([flips{k} '_35 = - ' flips{k} '_35;']);
%     eval([flips{k} '_10 = - ' flips{k} '_10;']);
%     eval([flips{k} '_10_35 = - ' flips{k} '_10_35;']);
% end

%% recalculate some fields using COARE output

% interface
Tskin       = b1.Tsnk + dT_warm_to_skin - dT_skin;
Tskin_10    = b10.Tsnk + dT_warm_to_skin_10 - dT_skin_10;

% recalculate upwelling radiative fluxes
% lw_up       = lw_net - b1.lw_dn; %%% this is REALLY NOISY
lw_up_10    = lw_net_10 - b10.lw_dn;
% sw_up       = sw_net - b1.sw_dn; %%% this is REALLY NOISY
sw_up_10    = sw_net_10 - b10.sw_dn;

sw_up = interp1(b10.jd, sw_up_10, b1.jd);
lw_up = interp1(b10.jd, lw_up_10, b1.jd);

% net heat flux: heating into the ocean
hnet        = sw_net + lw_net - hs - hl - hrain;
hnet_10     = sw_net_10 + lw_net_10 - hs_10 - hl_10 - hrain_10;

% this could be dTadj or dTadjK, but might as well do it as theta to be
% absolutely correct. The theta version is only 0.0022 greater. 
p_tq = b1.psealevel - 0.125*zt;
p_tq_10 = b10.psealevel - 0.125*zt;

theta10     = (b1.Ta+C2K).*(1000./p_tq).^(Rgas/cpa);  
theta0      = (Tskin+C2K).*(1000./b1.psealevel).^(Rgas/cpa);  
dtheta      = theta0 - theta10;

theta10_10  = (b10.Ta+C2K).*(1000./p_tq_10).^(Rgas/cpa);  
theta0_10   = (Tskin_10+C2K).*(1000./b10.psealevel).^(Rgas/cpa);  
dtheta_10   = theta0_10 - theta10_10;

tsr_son = tsr + 0.51*(b1.Ta+C2K).*qsr;
tsr_son_10 = tsr_10 + 0.51*(b10.Ta+C2K).*qsr_10;

% fields to replace now that they are checked, these already appeared in
% the original v0 b1 and b10 datasets. 
more_fields = {'hnet';'lw_up';'sw_up';'Tskin';'dtheta';'theta0';'theta10';'tsr_son'};
af = sort([cfields; more_fields]);
 
%% Example: name changing if needed...
%%% if the old version had a different name but you want to compare it to the current version now
% if isfield(b1,'SSQ') == 1
%     b1.qs = b1.SSQ;
%     b10.qs = b10.SSQ;
%     b1 = rmfield(b1,'SSQ');
%     b10 = rmfield(b10,'SSQ');
% end

%% plots
plot_checks = 0;
if plot_checks == 1
%%% now check to see before after of fluxes. 
figure;
counter = 1;
type_plot = 'coare';

cfields_s = cell(length(cfields),1);
for i = 1:length(cfields)
    cfields_s(i) = {[char(cfields(i)) '_s']};
end

if  strcmp(type_plot,'coare') == 1
    thefields = sort(cfields); %%% if you want to plot only coare output
elseif strcmp(type_plot,'all') == 1
    thefields = af; %%% if you want to plot all vars for entire experiment
end
nc = length(thefields);
for i = 1:4:nc
    clf;
    for j = 1:4
       if (i+j-1) <= nc
           subplot(2,2,j); hold on;
           eval(['newvar = ' thefields{i+j-1} '_10;']);
           plot(b10.t, newvar,'o');
           grid on;
           var_name = {strrep(thefields{i+j-1},'_',' ')};
           title(var_name);
           xlim([min(b10.t) max(b10.t)]);
           datetick('x','mm/dd','keeplimits');
           grid on;
           % plot coare v 3.5 if doing the coare plots
           if strcmp(type_plot,'coare') == 1
               eval(['oldvar = ' thefields{i+j-1} '_10_35;']);
               plot(b10.t, oldvar,'x');
               legend('3.6','3.5');
           % plot prior b10 fields before any corrections were made herein
           elseif strcmp(type_plot,'all') == 1
                if isfield(b10,thefields{i+j-1}) == 1
                plot(b10.t, b10.(thefields{i+j-1}),'.k','markersize',5);
                legend('3.6','3.6 prior','location','best');
                end
           % if neither coare 3.5 fields or prior b10 fields exist, just 
           %      continue plotting only new coare 3.6 fields     
           else
                legend('3.6');
           end
       end
    end
   print(graphdevice,[path_fix_plots '/' type_plot '_' sprintf('%i',counter) 'b_' cruise graphformat]);
   counter = counter + 1;
end  % for all vars

end %% if make plots

%% save new COARE and post-COARE derived values in stuctures now that they've been checked
for j = 1:length(af)
    eval(['b1.(af{j}) = ' af{j} ';']);
    eval(['b10.(af{j}) = ' af{j} '_10;']);
end

% save new ship data fields for completeness
b10.lw_up_s = b10.lw_net_s - b10.lw_dn_s;
b10.sw_up_s = b10.sw_net_s - b10.sw_dn_s;
b10.hnet_s = b10.sw_net_s + b10.lw_net_s - b10.hs_s - b10.hl_s - b10.hrain_s;
b1.lw_up_s = b1.lw_net_s - b1.lw_dn_s;
b1.sw_up_s = b1.sw_net_s - b1.sw_dn_s;
b1.hnet_s = b1.sw_net_s + b1.lw_net_s - b1.hs_s - b1.hl_s - b1.hrain_s;
b1.Tskin_s = b1.Tsea_s + b1.dT_warm_to_skin_s - b1.dT_skin_s;
b10.Tskin_s = b10.Tsea_s + b10.dT_warm_to_skin_s - b10.dT_skin_s;

%% fix a few PISTON cruise-specific things
b1.Ta_s(b1.Ta_s < 20) = nan;
b1.psealevel_s(b10.psealevel_s < 999) = nan;

%% elapsed time since Jan 1 of beginning of year... computed instead in save_nc.m
% t_day1 = datenum(yr,1,1, 0,0,0);
% t_day0 = datenum(yr,0,0, 0,0,0);
% 
% b10.t_since     = (b10.t - t_day1)*24*3600;     % subtracts Jan 1 of this year from t and converts days to seconds.
% b1.t_since      = (b1.t - t_day1)*24*3600;      % subtracts Jan 1 of this year from t and converts days to seconds.
% 
% % jd is vector PSL decimal year day, where Jan 1  = day 1.
% t_check = jd + t_day0;  % t = matlab date time, in units of decimal days since January 1, 0000
% te1 = (t - t_day1)*24*3600; % te is seconds elapsed since Jan 1 of that year, in this case 2019   
% te2 = (jd + t_day0 - t1)*24*3600; % same
% te3 = (jd - 1)*24*3600; % same
% jd_check = (te3/(24*3600)) +1 ; % PSL year date in our normal programs
% 
% disp(datestr(t(1)),0);
% disp(datestr(jd(1)-1),0);

%% save corrected data
% organize structures for saving
% put all variable names in alphabetical order;
b1 = orderfields(b1);
b10 = orderfields(b10);
save(m_outfile_1, 'b1');
save(m_outfile_10, 'b10');

%%% make some plots:
plot_more = 0;
if plot_more == 1

%%% skin T

figure;
subplot(2,1,1); hold on;
plot(b1.t, b1.dT_skin,'-','color',rgb('gray'));
legend('dT cool skin','location','southwest')
grid on;
ylabel('^o C');
title('PISTON 2019');
datetick('x','mm/dd');
xlim([min(b1.t) max(b1.t)]);

subplot(2,1,2);
hold on;
plot(sea.t, sea.SSST,'xb');
plot(b1.t, b1.Tskin,'--r',b1.t, b1.Tsnk,'-g');
ylabel('^o C');
grid on;
datetick('x','mm/dd');
xlim([min(b1.t) max(b1.t)]);
legend('ROSR T_{skin}','COARE T_{skin}','snake T_{0.05m}','location','northwest')

print(graphdevice,[path_fix_plots '/warm_skin_1_' cruise graphformat]);


case1end    = datenum(2019,9,9,20,0,0);
case1start  = datenum(2019,9,8,15,0,0);
case2end    = datenum(2019,9,23,20,0,0);
case2start  = datenum(2019,9,22,20,0,0);
case1 = find(b1.t <= case1end & b1.t >= case1start);
case2 = find(b1.t <= case2end & b1.t >= case2start);
case1_rosr = find(sea.t <= case1end & sea.t >= case1start);
case2_rosr = find(sea.t <= case2end & sea.t >= case2start);

case1end_f      = datenum(2019,9,9,6,0,0);
case1start_f    = datenum(2019,9,9,4,0,0);
case2end_f      = datenum(2019,9,23,6,0,0);
case2start_f    = datenum(2019,9,23,4,0,0);
case1_f = find(b1.t <= case1end_f & b1.t >= case1start_f);
case2_f = find(b1.t <= case2end_f & b1.t >= case2start_f);
case1_rosr_f = find(sea.t <= case1end_f & sea.t >= case1start_f);
case2_rosr_f = find(sea.t <= case2end_f & sea.t >= case2start_f);


figure;
subplot(2,2,1);
yyaxis right; hold on;
plot(sea.t(case1_rosr), sea.SSST(case1_rosr),'xb');
plot(b1.t(case1), b1.Tskin(case1),'--r',b1.t(case1), b1.Tsnk(case1),'-g');
set(gca,'ycolor','k')
ylabel('^o C');
yyaxis left;
hold on;
plot(b1.t(case1), b1.dT_skin(case1),'-','color',rgb('gray'));
set(gca,'ycolor',rgb('gray'))
legend('dT cool skin','ROSR T_{skin}','COARE T_{skin}','snake T_{0.05m}','location','southeast')
grid on;
ylabel('^o C');
title('PISTON 2019 Sept 9-10, DOY 252-253');
datetick('x','hh');
axis tight;
xlabel('UTC')

subplot(2,2,3);
yyaxis left;
plot(b1.t(case1), b1.prate(case1), 'o','color',rgb('gray'))
set(gca,'ycolor',rgb('gray'))
ylabel('mm hr^{-1}');
yyaxis right; hold on;
plot(sea.t(case1_rosr), sea.SSST(case1_rosr),'xb');
plot(b1.t(case1), b1.Tskin(case1),'--r',b1.t(case1), b1.Tsnk(case1),'-g');
ylabel('^o C');
set(gca,'ycolor','k')
grid on;
datetick('x','hh');
axis tight;
legend('rain rate','ROSR T_{skin}','COARE T_{skin}','snake T_{0.05m}','position',[0.3,0.15, 0.1,0.1])
xlabel('UTC');

subplot(2,2,2);
yyaxis right; hold on;
plot(sea.t(case2_rosr), sea.SSST(case2_rosr),'xb');
plot(b1.t(case2), b1.Tskin(case2),'--r',b1.t(case2), b1.Tsnk(case2),'-g');
set(gca,'ycolor','k')
ylabel('^o C');
yyaxis left;
hold on;
plot(b1.t(case2), b1.dT_skin(case2),'-','color',rgb('gray'));
set(gca,'ycolor',rgb('gray'))
% legend('dT cool skin','ROSR T_{skin}','COARE T_{skin}','snake T_{0.05m}','location','southeast')
grid on;
ylabel('^o C');
title('PISTON 2019 Sept 22-23, DOY 265-266');
datetick('x','hh');
axis tight;
xlabel('UTC')


subplot(2,2,4);
yyaxis left;
plot(b1.t(case2), b1.prate(case2), 'o','color',rgb('gray'))
set(gca,'ycolor',rgb('gray'))
ylabel('mm hr^{-1}');
yyaxis right; hold on;
plot(sea.t(case2_rosr), sea.SSST(case2_rosr),'xb');
plot(b1.t(case2), b1.Tskin(case2),'--r',b1.t(case2), b1.Tsnk(case2),'-g');
ylabel('^o C');
set(gca,'ycolor','k')
grid on;
datetick('x','hh');
axis tight;
% legend('rain rate','ROSR T_{skin}','COARE T_{skin}','snake T_{0.05m}','position',[0.3,0.15, 0.1,0.1])
xlabel('UTC');

print(graphdevice,[path_fix_plots '/warm_skin_2_' cruise graphformat]);


%%%% focused in on time when dT_skin < 0 for warm skin

figure;
subplot(1,2,1);
yyaxis right; hold on;
plot(sea.t(case1_rosr_f), sea.SSST(case1_rosr_f),'xb');
plot(b1.t(case1_f), b1.Tskin(case1_f),'--r',b1.t(case1_f), b1.Tsnk(case1_f),'-g');
set(gca,'ycolor','k')
ylabel('^o C');
yyaxis left;
hold on;
plot(b1.t(case1_f), b1.dT_skin(case1_f),'-','color',rgb('gray'));
set(gca,'ycolor',rgb('gray'))
legend('dT cool skin','ROSR T_{skin}','COARE T_{skin}','snake T_{0.05m}','location','southeast')
grid on;
ylabel('^o C');
title('PISTON 2019 Sept 10, DOY 253 zoom');
datetick('x','hh:MM');
axis tight;
xlabel('UTC')

subplot(1,2,2);
yyaxis right; hold on;
plot(sea.t(case2_rosr_f), sea.SSST(case2_rosr_f),'xb');
plot(b1.t(case2_f), b1.Tskin(case2_f),'--r',b1.t(case2_f), b1.Tsnk(case2_f),'-g');
set(gca,'ycolor','k')
ylabel('^o C');
yyaxis left;
hold on;
plot(b1.t(case2_f), b1.dT_skin(case2_f),'-','color',rgb('gray'));
set(gca,'ycolor',rgb('gray'))
% legend('dT cool skin','ROSR T_{skin}','COARE T_{skin}','snake T_{0.05m}','location','southeast')
grid on;
ylabel('^o C');
title('PISTON 2019 Sept 23, DOY 266 zoom');
datetick('x','hh:MM');
axis tight;
xlabel('UTC')

print(graphdevice,[path_fix_plots '/warm_skin_3_' cruise graphformat]);



end % end more plots

plot_warm_layer = 0;
if plot_warm_layer == 1
    
figure;
subplot(2,1,1);
hold on;
plot(sea.t, sea.SSST,'xb');
plot(b1.t, b1.Tskin,'--r',b1.t, b1.Tsnk,'-g');
ylabel('^o C');
grid on;
datetick('x','mm/dd');
xlim([min(b1.t) max(b1.t)]);
legend('ROSR T_{skin}','COARE T_{skin}','snake T_{0.05m}','location','northwest')

subplot(2,1,2);
hold on;
yyaxis left;
plot(sea.t, sea.SSST,'xb');
set(gca,'ycolor',rgb('blue'))
ylabel('^o C');
yyaxis right;
plot(b1.t, b1.U10,'-','color',rgb('gray'));
set(gca,'ycolor',rgb('gray'))
ylabel('U_{10} m s^{-1}');
grid on;
datetick('x','mm/dd');
xlim([min(b1.t) max(b1.t)]);
legend('ROSR T_{skin}','U_{10}','location','northwest')

print(graphdevice,[path_fix_plots '/warm_layer_1_' cruise graphformat]);


casetime = [datenum(2019,9,22,0,0,0) datenum(2019,9,24,23,59,0)];

figure;
subplot(2,1,1);
hold on;
plot(sea.t, sea.SSST,'xb');
plot(b1.t, b1.Tskin,'--r',b1.t, b1.Tsnk,'-g');
xlim([casetime]);
ylabel('^o C');
grid on;
datetick('x','mm/dd','keeplimits');
legend('ROSR T_{skin}','COARE T_{skin}','snake T_{0.05m}','location','northwest')

subplot(2,1,2);
hold on;
yyaxis left;
plot(sea.t, sea.SSST,'xb');
set(gca,'ycolor',rgb('blue'))
xlim([casetime]);
ylabel('^o C');
yyaxis right;
plot(b1.t, b1.U10,'-','color',rgb('gray'));
set(gca,'ycolor',rgb('gray'))
xlim([casetime]);
ylabel('U_{10} m s^{-1}');
grid on;
datetick('x','mm/dd','keeplimits');
legend('ROSR T_{skin}','U_{10}','location','northwest')

print(graphdevice,[path_fix_plots '/warm_layer_2_' cruise graphformat]);

end

function g=grv(lat)
% computes g [m/sec^2] given lat in deg
gamma=9.7803267715;
c1=0.0052790414;
c2=0.0000232718;
c3=0.0000001262;
c4=0.0000000007;
phi=lat*pi/180;
x=sin(phi);
g=gamma*(1+c1*x.^2+c2*x.^4+c3*x.^6+c4*x.^8);
end