% a3_final_filtered_FATIMA_2022.m

% Load motcorr da table, compute fluxes, output filtered results.
% Save filtered output in 10 min and hourly flux tables.
% Create analysis and diagnostic plots.

% This script creates placeholder variables in the output tables
% for ship measurements (e.g. SCS) but none of the code to read
% ship data or compute the output fields is implemented in this
% version.

% RUN THIS SCRIPT AFTER RUNNING a2_motcorr_flux_FATIMA_2022.m

% BWB 7/2022

%% initialize & set matlab paths for your system
close all;
fclose all;
clear;
warning ('off','MATLAB:MKDIR:DirectoryExists');

saveflx = true; % set true to save flux result tables in mat and ascii formats
savefigs = true;        % set true to save plots
ptitle = 'FATIMA 2022'; % title string for plots (no underscore chr)
cruise = 'FATIMA_2022'; % string for file names

% system specific path defs
data_drive = '/Users/byronb/Raw_Data_Files.noindex';
% matlab script path
path_prog = fullfile(data_drive,cruise,'_matlab');
cd(path_prog);
addpath(genpath(path_prog));

% COARE options
jcool = 0;  % set to 1 to compute cool skin
% warm layer is only applied if tsg is used for tsea in COARE loop
jwarm = 0;  % set to 1 for warm layer estimate if using ship 5m SST
zref = 10;  % user defined reference height

% set flag for decorrelation status of raw da file
decorr = true;
if decorr
    dcname ='decorr';
else
    dcname = 'std';
end

% define path to 'processed' folder with da subfolder(s)
da_path = fullfile(data_drive,cruise,'processed');

% define paths to folders for saving da plots
png_path = fullfile(data_drive,cruise,'processed_images','final_plots_png');
mkdir(png_path); % create if it doesn't exist

% set limits for various filtering criteria
rdir_lo     = -90;  % rel wind direction limits, deg from bow
rdir_hi     = 90;
rain_lim    = 1;    % max rain rate, mm/hr
vplat_lim   = 0.8;  % max vplat std dev  m/s
sigu_lim    = 1.5;  % max ship speed std dev  m/s
sigh_lim    = 5;    % max ship heading std dev, deg
badSon_lim  = 100;  % max# bad sonic points per 10min interval
missing_lim = 100;  % max# missing sonic points per 10min interval
ushp_lim    = 10;   % max ship speed, m/s
tilt_hi     = 10;   % max flow tilt angle, degrees
tilt_lo     = 0;    % min flow tilt angle, degrees
sig_rwd_lim = 15;   % max std dev rwdir
AGC_lim     = 60;   % max licor AGC
CO2_LoLim   = 200;  % min CO2 ppm
CO2_HiLim   = 700;  % max CO2 ppm
sigCO2_lim  = 20;   % max sigma CO2
q_lic_LoLim = 0.1;  % min q licor
q_lic_HiLim = 26;   % max q licor
sigq_LoLim  = 0;    % min sigma q licor
sigq_HiLim  = 3;    % max sigma q licor
wq_LoLim    = -0.2; % min wq licor
wq_HiLim    = 0.2;  % max wq licor
wu_LoLim    = -0.8; % min wu
wu_HiLim    = 0.1;  % max wu
wv_LoLim    = -0.2; % min wv
wv_HiLim    = 0.2;  % max wv
wT_LoLim    = -0.25;% min wT
wT_HiLim    = 0.25; % max wT
Tnoise_lim  = 0.001;% max Ts noise
Tson_var_lim = 0.15;% max variance in Tsonic
cu_LoLim    = 0;    % min cu2
cu_HiLim    = 1;    % max cu2
cw_LoLim    = 0;    % min cw2
cw_HiLim    = 1;    % max cw2
ct_LoLim    = 1e-6; % min ct2
ct_HiLim    = 0.6;  % max ct2
cq_LoLim   = 1e-4; % min cql2
cq_HiLim   = 2;    % max cql2

% choose ID method, 'a', 'b', or 'ab'
% a: compute Cx2 from median over inertial subrange in smoothed spectra
% b: compute Cx2 from polynomial fit over inertial subrange in smoothed spectra
% ab: average of a and b methods
IDflag = 'b';

% measurement heights as measured in July 2022 installation
zu = 18.60;        % NOAA R3 sonic height (m)
zt = 18;           % air T measurement height (m), (WXT for this cruise)
zq = 18;           % air q measurement height (m), (WXT for this cruise)
zP = 16.25;        % PSL PTB pressure sensor height on bridge roof (m)
zWXT = 18;         % PSL WXT on bow tower (m)
da.zu_ship = NaN;  % ship sensor heights (m) (none this cruise)
da.zt_ship = NaN;
da.zq_ship = NaN;
ztsg = NaN;        % no ship sst/sal sensor this cruise
zts = 0.05;        % approximate depth of sea snake, m
zi = 600;          % default boundary layer depth, m

% constants
C2K = 273.15;               % Temperature conversion constant
d2r = pi/180;               % degrees to radians
r2d = 180/pi;               % radians to degrees
Rgas = 287.1;               % Gas Const: Pa m-3 kg-1 K-1 for dry air
Rgas_universal = 8.314472;	% Gas Const: Pa m3 K-1 mol-1
Mw = 18.01528;              % molar mass of H20, g/mole
Md = 28.964;                % molar mass dry air, g/mole
epsilon = Mw/Md;            % Mw/Md ratio = 0.622
cpa  = 1004.67;             % heat capacity dry air

% lat/lon limits based on for filtering gps outliers
Lonmin = 292; Lonmax = 313; Latmin = 38; Latmax = 50;
LatLonLim = [Lonmin, Lonmax, Latmin, Latmax];

%% load 10-min da data table from motcorr processing
%{
                        TIMESTAMP VARIABLES
    1     jd              decimal_day   decimal day of year timestamp
    2     dnum            decimal_day   matlab datenum timestamp
    3     year            none          year (cols 3:8 are date;time vector)
    4     month           none          month
    5     day             none          day of month
    6     hour            none          hour of day; 0;23
    7     minute          none          minute of hour
    8     second          none          second (zeros for this file)
                        COARE 3.5 SHIP (SCS) INPUT FIELDS, cols 9:21
    9     wspd_ship       m/s           true wind speed; ship; m/s
    10    zu_ship         m             wspd measurement height; ship; m
    11    ta_ship         degC          air temperature; ship; ¬degC
    12    zt_ship         m             ta measurement height; ship; m
    13    rh_ship         percent       relative humidity; ship; %
    14    zq_ship         m             rh measurement height; ship; m
    15    Pmb_ship        mb            surface air pressure; ship; mb
    16    sst_ship        degC          SST ship TSG; ¬degC
    17    rs_ship         W/m2          solar radiative flux; ship; W/m2
    18    rl_ship         W/m2          IR radiative flux; ship; W/m2
    19    lat_ship        degN          latitude; ship; degN
    20    zi_ship         m             MBL depth; NaN for default 600m
    21    rain_ship       mm/hr         rain rate; ship; mm/hr
                        COARE 3.5 NOAA PSL INPUT FIELDS, cols 22:34
    22    wspd_psl        m/s           true wind speed; NOAA PSL sonic; m/s
    23    zu_psl          m             wspd measurement height; m
    24    ta_psl          degC          air temperature; NOAA PSL Vaisala HMT; ¬degC
    25    zt_psl          m             ta measurement height; m
    26    rh_psl          percent       relative humidity; NOAA PSL Vaisala HMT; %
    27    zq_psl          m             rh measurement height; m
    28    Pmb_psl         mb            surface air pressure; NOAA PSL Vaisala PTB; mb
    29    sst_psl         degC          SST; seasnake; ¬degC
    30    rs_psl          W/m2          solar radiative flux; NOAA PSL psp; W/m2
    31    rl_psl          W/m2          IR radiative flux; NOAA PSL pir; W/m2
    32    lat_psl         degN          latitude; NOAA PSL gps; degN
    33    zi_psl          m             MBL depth; NaN for default 600m
    34    rain_psl        mm/hr         rain rate; org; NOAA PSL; mm/hr
                        OTHER WIND VARIABLES
    35    wdir_ship       deg           ship mean true wind direction; deg
    36    wdir_ship_std   deg           std dev ship true wind direction; deg
    37    rwspd_ship      m/s           ship rel wind speed
    38    rwspd_ship_std  m/s           ship rwspd std dev
    39    rwdir_ship      deg           ship rel wind dir; +/; 180 deg
    40    rwdir_ship_std  deg           ship rwdir std dev
    41    wdir_psl        deg           true wind direction from N;E components; NOAA PSL sonic; deg
    42    rwspd_psl       m/s           streamwise relative wind speed from NOAA PSL sonic; m/s
    43    rwspd_psl_std   m/s           std dev rwspd from NOAA PSL sonic; m/s
    44    rwdir_psl       deg           relative wind direction from NOAA PSL sonic; deg
    45    rwdir_psl_std   deg           std dev rwdir from NOAA PSL sonic; deg
    46    Un              m/s           N component; corrected true windspeed; NOAA PSL sonic; m/s
    47    Ue              m/s           E component; corrected true windspeed; NOAA PSL sonic; m/s
    48    Ubar            m/s           streamwise true wspd; NOAA PSL sonic; no flow dist. correction; m/s
                        NOAA PSL WXT VARIABLES
    49    rwdir_wxt       deg           rel wind direction ;0;360 deg from bow; NOAA PSL WXT
    50    rwspd_wxt       m/s           wind speed; m/s; NOAA PSL WXT
    51    ta_wxt          degC          air temperature; ¬degC; NOAA PSL WXT
    52    rh_wxt          percent       RH; %; NOAA PSL WXT
    53    Pmb_wxt         mb            sea level air pressure; mb or hPa; NOAA PSL WXT
    54    rain_wxt        mm/hr         rain rate; mm/hr; NOAA PSL WXT
    55    wspd_wxt        m/s           true wind speed; m/s; NOAA PSL WXT
    56    wdir_wxt        deg           true wind dir; m/s; NOAA PSL WXT
    57    Urel_wxt        m/s           U relative; m/s; NOAA PSL WXT
    58    Vrel_wxt        m/s           V relative; m/s; NOAA PSL WXT
                        OTHER MET/OCEAN VARIABLES
    59    qa_psl          g/kg          air specific humidity; NOAA PSL sensors; g/kg
    60    qa_ship         g/kg          air specific humidity; ship sensors; g/kg
    61    qa_wxt          g/kg          WXT specific humidity; g/kg
    62    qs_psl          g/kg          sea surface saturated specific humidity; seasnake SST; g/kg
    63    qs_ship         g/kg          sea surface saturated specific humidity; ship SST; g/kg
    64    qs_rosr         g/kg          sea surface saturated specific humidity; ROSR SST; g/kg
    65    rhoa_psl        kg/m3         moist air density from NOAA PSL sensors; kg/m3
    66    sal_ship        ppth          surface seawater salinity; ship TSG; ppth
    67    sst_rosr        degC          ROSR radiometric SST; ¬degC
    68    aspir_on        none          Vaisala HMT ventilation on/off parameter
    69    org_carrier     volts         ORG rainrate QA parameter; volts
                       NAVIGATION VARIABLES
    70    sog_ship        m/s           mean sog; ship gps; m/s
    71    sog_ship_std    m/s           std sog; ship gps; m/s
    72    sog_ship_min    m/s           minimum sog; ship gps; m/s
    73    sog_ship_max    m/s           maximum sog; ship gps; m/s
    74    cog_ship        deg           mean cog; ship gps; deg
    75    cog_ship_std    deg           std cog; ship gps; deg
    76    sog_psl         m/s           mean sog; m/s; NOAA PSL gps; deg
    77    sog_psl_std     m/s           std sog; NOAA PSL gps; m/s
    78    sog_psl_min     m/s           minimum sog; NOAA PSL gps; deg
    79    sog_psl_max     m/s           maximum sog; NOAA PSL gps; deg
    80    cog_psl         deg           mean cog; deg; NOAA PSL gps; deg
    81    cog_psl_std     deg           std cog; NOAA PSL gps; deg
    82    lrg             deg           heading; ship gyro; deg
    83    lrg_std         deg           std dev heading; ship gyro; deg
    84    hed_psl         deg           heading; NOAA PSL gps; deg
    85    hed_psl_std     deg           std dev heading; NOAA PSL gps; deg
    86    lon_ship        degE          longitude; ship gps
    87    lon_psl         degE          latitiude; NOAA PSL gps
    88    odec            m/s           odec (ship speed through water); m/s
    89    odec_std        m/s           stdev odec; m/s
    90    odec_min        m/s           minimum odec; m/s
    91    odec_max        m/s           maximum odec; m/s
                        COARE 3.5 ship data output fields
    92    usb_ship        m/s           friction velocity; ship; m/s
    93    taub_ship       N/m2          surface wind stress; ship; N/m2
    94    hsb_ship        W/m2          sensible heat flux; ship; W/m2
    95    hlb_ship        W/m2          latent heat flux; ship; W/m2
    96    hlwebb_ship     W/m2          hlb 'webb' correction; W/m2
    97    tsb_ship        degK          t;star; ship; ¬degC
    98    qsb_ship        g/kg          q;star; ship; g/kg
    99    zo_ship         m             surface roughness length; ship; m
    100   zot_ship        m             surface temperature roughness length; ship; m
    101   zoq_ship        m             surface moisture roughness length; ship; m
    102   Cd_ship         none          drag coefficient; ship
    103   Ch_ship         none          sensible heat transfer coefficient; ship
    104   Ce_ship         none          water vapor transfer coerfficient; ship
    105   L_ship          m             Obukhov length; ship; m
    106   zet_ship        none          z/L; ship @ zu_ship
    107   dter_ship       degC          cool skin temperature depression; ship; ¬degC
    108   dqer_ship       degC          cool skin saturation humidity depression; g/kg
    109   tkt_ship        m             cool skin thickness; ship; m
    110   RF_ship         W/m2          rain heat flux; ship; W/m2
    111   Cdn_10_ship     none          10m neutral drag coefficient; ship
    112   Chn_10_ship     none          10m neutral heat transfer coefficient; ship
    113   Cen_10_ship     none          10m neutral moisture transfer coefficient; ship
                        COARE 3.5 NOAA PSL data output fields
    114   usb_psl         m/s           friction velocity; NOAA PSL; m/s
    115   taub_psl        N/m2          surface wind stress; NOAA PSL; N/m2
    116   hsb_psl         W/m2          sensible heat flux; NOAA PSL; W/m2
    117   hlb_psl         W/m2          latent heat flux; NOAA PSL; W/m2
    118   hlwebb_psl      W/m2          hlb 'webb' correction; NOAA PSL; W/m2
    119   tsb_psl         degK          t;star; NOAA PSL; ¬degK
    120   qsb_psl         g/kg          q;star; NOAA PSL; g/kg
    121   zo_psl          m             surface roughness length; NOAA PSL; m
    122   zot_psl         m             surface temperature roughness length; NOAA PSL; m
    123   zoq_psl         m             surface moisture roughness length; NOAA PSL; m
    124   Cd_psl          none          drag coefficient; NOAA PSL
    125   Ch_psl          none          sensible heat transfer coefficient; NOAA PSL
    126   Ce_psl          none          water vapor transfer coerfficient; NOAA PSL
    127   L_psl           m             Obukhov length; NOAA PSL; m
    128   zet_psl         none          z/L; NOAA PSL @ zu_psl
    129   dter_psl        degC          cool skin temperature depression; NOAA PSL; ¬degC
    130   dqer_psl        degC          cool skin humidity depression; NOAA PSL; ¬degC
    131   tkt_psl         m             cool skin thickness; NOAA PSL; m
    132   RF_psl          W/m2          rain heat flux; NOAA PSL; W/m2
    133   Cdn_10_psl      none          10m neutral drag coefficient; NOAA PSL
    134   Chn_10_psl      none          10m neutral heat transfer coefficient; NOAA PSL
    135   Cen_10_psl      none          10m neutral moisture transfer coefficient; NOAA PSL
                        FLUX VARIABLES
    136   wu_cov          m2/s2         wu covariance; Edson correction; m2/s2
    137   wu_csp          m2/s2         wu cospectrum integral; Edson correction; m2/s2
    138   wv_cov          m2/s2         wv covariance; Edson correction; m2/s2
    139   wv_csp          m2/s2         wv cospectrum integral; Edson correction; m2/s2
    140   wTs_cov         degC_m/s      wTsonic covariance; Edson correction; ¬degC m/s
    141   wTs_csp         degC_m/s      wTsonic cospectrum integral; Edson correction; m2/s2
    142   wT_cov          degC_m/s      wT covariance; ¬degC m/s
    143   wT_csp          degC_m/s      wT cospectrum integral; Edson correction; m2/s2
    144   wq_cov          g/kg_m/s      wq covariance; Licor; g/kg m/s
    145   wq_csp          g/kg_m/s      wq cospectrum integral; Licor; g/kg m/s
    146   hs_cov          W/m2          sensible heat flux from covariance; W/m2
    147   hs_csp          W/m2          sensible heat flux from cospectrum; W/m2
    148   hl_cov          W/m2          latent heat flux from covariance; uncorrected; W/m2
    149   hl_csp          W/m2          latent heat flux from cospectrum; uncorrected; W/m2
    150   hlc_cov         W/m2          webb corrected hl_cov; W/m2
    151   hlc_csp         W/m2          webb corrected hl_csp; W/m2
    152   Cu2a            none          Cu structure func parameter; method a
    153   Cu2b            none          Cu structure func parameter; method b
    154   Cw2a            none          Cw structure func parameter; method a
    155   Cw2b            none          Cw structure func parameter; method b
    156   Ct2a            none          Ct structure func parameter; method a
    157   Ct2b            none          Ct structure func parameter; method b
    158   Cq2a            none          Cq structure func parameter; method a
    159   Cq2b            none          Cq structure func parameter; method b
    160   usid_a          m/s           ustar from inertial diss.; method a
    161   usid_b          m/s           ustar from inertial diss.; method b
    162   tsid_a          degK          tstar from inertial diss.; method a
    163   tsid_b          degK          tstar from inertial diss.; method b
    164   qsid_a          g/kg          qstar from inertial diss.; method a
    165   qsid_b          g/kg          qstar from inertial diss.; method b
    166   Lid_a           m             L from inertial diss.; method a
    167   Lid_b           m             L from inertial diss.; method b
                        TURBULENCE STATISTICS AND QA/QC PARAMETERS
    168   uvar            m2/s2         u variance; Edson correction; m2/s2
    169   vvar            m2/s2         v variance; Edson correction; m2/s2
    170   wvar            m2/s2         w variance; Edson correction; m2/s2
    171   tvar            degC^2        Tsonic variance; ¬degC^2
    172   qvar            g2/kg2        Licor q variance; g/kg
    173   azm             deg           azimuth angle from first rotation; radians
    174   tilt            deg           tilt angle from first rotation; radians
    175   fsonic          Hz            sonic anemometer sampling frequency; Hz
    176   missing_son     none          missing sonic wind data lines per 10 min interval @ fsonic
    177   badPnts_son     none          #NaNs in sonic data per 10 min interval @ fsonic
    178   flicor          Hz            Licor 7500 sampling frequency; Hz
    179   missing_lic     none          missing licor data lines per 10 min interval @ flicor
    180   fmotion         Hz            motion system sampling frequency; Hz
    181   missing_mot     none          missing motion data lines per 10 min interval @ fmotion
    182   fheading        Hz            Hemisphere heading system sampling frequency; Hz
    183   missing_hed     none          missing heading data lines per 10 min interval @ fheading
    184   uplat_std       m/s           U platform velocity std dev; Edson correction; m2/s2
    185   vplat_std       m/s           V platform velocity std dev; Edson correction; m2/s2
    186   wplat_std       m/s           W platform velocity std dev; Edson correction; m2/s2
    187   Ts_noise        degC^2        Sonic temperature noise variance; ¬degC^2
                    OTHER LICOR 7500 QA/QC VARIABLES - 'good' is ADC<65
    188   h2o_lic_all     mmol/m3       mean H2O; all data; mmol/m3
    189   h2o_lic_good    mmol/m3       mean H2O; good data; mmol/m3
    190   co2_lic_all     ppm           mean CO2; all data; ppm
    191   co2_lic_good    ppm           mean CO2; good data; ppm
    192   Pmb_lic         mb            pressure from Licor box; mb
    193   T_lic           degC          temperature from Licor box; ¬degC
    194   h2o_std_all     mmol/m3       std dev H2O; all data; mmol/m3
    195   h2o_std_good    mmol/m3       std dev H2O; good data; mmol/m3
    196   co2_std_all     ppm           std dev CO2; all data; ppm
    197   co2_std_good    ppm           std dev CO2; good data; ppm
    198   AGC_lic         none          Licor diagnostic; 50-100; 50=good
    199   q_lic           g/kg          Licor specific humidity
%}
load(fullfile(da_path,['da_',dcname],['FATIMA_2022_da_',dcname,'.mat']),'da');
da_vars = da.Properties.VariableNames;
da_units = da.Properties.VariableUnits;
da_descript = da.Properties.VariableDescriptions;

%% define variables from da array
%  date limits for plots etc.
dt = datetime(da.dnum,'ConvertFrom','datenum');
stdt = dt(1); endt = dt(end) + minutes(10);

% fix up a few 10-min variables
da.lon_psl(da.lon_psl<0) = da.lon_psl(da.lon_psl<0) + 360;  % reformat negative W longitudes
da.cog_psl(da.cog_psl<0) = da.cog_psl(da.cog_psl<0)+360;
da.cog_ship(da.cog_ship<0) = da.cog_ship(da.cog_ship<0)+360;
da.missing_son(isnan(da.missing_son)) = 6000;

%% make adjustments, filter outliers, etc.
% Apply fixes to data and select best quality measurements
% for use in bulk flux calculations.  Edit the fixit subroutine
% to make corrections and remove outliers.
fixit_da1_FATIMA_2022

%% select inertial dissipation results from da file
switch IDflag
    case 'a'
        cu   =   da.Cu2a;   % cu^2 U method a
        cw   =   da.Cw2a;   % cw^2 W method a
        ct   =   da.Ct2a;   % ct^2 temperature method a
        cql  =   da.Cq2a;   % cq^2 H2O Licor method a
        usi  =   da.usid_a; % inertial dissipation scaling paramters, all method a
        tsi  =   da.tsid_a;
        qsi  =   da.qsid_a;
    case 'b'
        cu   =   da.Cu2b;   % cu^2 U method b
        cw   =   da.Cw2b;   % cw^2 W method b
        ct   =   da.Ct2b;   % ct^2 temperature method b
        cql  =   da.Cq2b;   % cq^2 H2O Licor method b
        usi  =   da.usid_b; % inertial dissipation scaling paramters, all method b
        tsi  =   da.tsid_b;
        qsi  =   da.qsid_b;
    case 'ab'                % average of a and b methods
        cu   =   (da.Cu2a + da.Cu2b)/2;
        cw   =   (da.Cw2a + da.Cw2b)/2;
        ct   =   (da.Ct2a + da.Ct2b)/2;
        cql  =   (da.Cq2a + da.Cq2b)/2;
        usi  =   (da.usid_a + da.usid_b)/2;
        tsi  =   (da.tsid_a + da.tsid_b)/2;
        qsi  =   (da.qsid_a + da.qsid_b)/2;
    otherwise               % default to case 'a' for any other input flag
        cu   =   da.Cu2a;   % cu^2 U method a
        cw   =   da.Cw2a;   % cw^2 W method a
        ct   =   da.Ct2a;   % ct^2 temperature method a
        cql  =   da.Cq2a;   % cq^2 H2O Licor method a
        usi  =   da.usid_a; % inertial dissipation scaling paramters, all method a
        tsi  =   da.tsid_a;
        qsi  =   da.qsid_a;
end

% average cu and cw for smoother structure function
cs = (cu + 0.75*cw)/2;
ii = find(isnan(cu));
cs(ii) = .75*cw(ii);
ii = find(isnan(cw));
cs(ii) = cu(ii);

%% define best bulk met values from filtered raw data
% edit this as needed, based on quality of bulk parameters
% reviewed in the fixit routine
wspd_best = da.wspd_psl; wspd_best(isnan(wspd_best)) = da.wspd_wxt(isnan(wspd_best));
wdir_best = da.wdir_psl; wdir_best(isnan(wdir_best)) = da.wdir_wxt(isnan(wdir_best));
rwdir_best = da.rwdir_psl; rwdir_best(isnan(rwdir_best)) = da.rwdir_wxt(isnan(rwdir_best));
ts_best = da.sst_rosr;
ta_best = da.ta_wxt;
rh_best = da.rh_wxt;
rs_best = da.rs_psl;
rl_best = da.rl_psl;
ii = isnan(da.Pmb_psl); % replace any missing PTB pressures with WXT
da.Pmb_psl(ii) = da.Pmb_wxt(ii);
P_best = da.Pmb_psl;
rain_best = da.rain_wxt;
qa_best = qair_p([ta_best,rh_best,P_best]);
qs_best = qsea_p([ts_best P_best]);
da.rhoa_psl = P_best.*100./(Rgas*(ta_best + C2K).*(1 + 0.61*qa_best/1000));
Le = (2.501-.00237*nanmedian1(ts_best))*1e6;
% set best nav variables
sog_best = da.sog_psl; std_sog_best = da.sog_psl_std;
cog_best = da.cog_psl;
hed_best = da.hed_psl; std_hed_best = da.hed_psl_std;
Ngps_best = sog_best.*cos(cog_best*d2r); % N component ship speed
Egps_best = sog_best.*sin(cog_best*d2r); % E component ship speed

%% correct U for ocean surface currents
if any(isfinite(da.odec))
    % compute surface current from da.odec, head, da.sog_ship, da.cog_ship
    % current direction is 'to' - oceanographic convention
    [cspd,cdir] = current(da.odec,hed_best,sog_best,cog_best);
    % compute wind w/respect to moving surface
    [wspd_sfc,wdir_sfc] = wind_current_adjust(wspd_best,wdir_best,cspd,cdir);
    % use current-relative winds for 'best'
    wspd_best = wspd_sfc; wdir_best = wdir_sfc;
    % append sfc-relative winds to da array for later averaging
    da = addvars(da,wspd_sfc,wdir_sfc);
else    % if no da.odec, set sfc-relative winds to NaN in da array
    wspd_sfc = NaN(size(wspd_best));
    wdir_sfc = NaN(size(wspd_best));
    da = addvars(da,wspd_sfc,wdir_sfc);
end

da_vars = [da_vars,{'wspd_sfc','wdir_sfc'}];
da_units = [da_units,{'m/s','deg'}];
da_descript = [da_descript,{'surface wspd corrected for ocean current',...
                            'surface wdir corrected for ocean current'}];
da.Properties.VariableUnits = da_units;
da.Properties.VariableDescriptions = da_descript;

%% compute estimate for clear sky radiative fluxes & plot
rlclr = (0.52 + 0.13/60*abs(da.lat_psl) + (0.082 - 0.03/60*abs(da.lat_psl)).*...
    sqrt(qa_best)).*(5.67e-8*(ta_best + C2K).^4);

rsclr = rs_clear(da.jd,nanmean1(P_best),qa_best,da.lat_psl,da.lon_psl);

figure('position',[1,1,800,475],'name','rsclr');
plot(dt,rsclr,'k--',dt,rs_best,'r-','LineWidth',1);
ylabel('Solar Flux (W/m^2)'); grid on; grid minor;
legend('clear sky','psp','location','best');
title('Downwelling Solar Radiative Flux');
if savefigs
    ppath = fullfile(png_path,['029_',cruise,'_Solar_Best_Clearsky.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['029_',cruise,'_Solar_Best_Clearsky.fig']));
end

figure('position',[1,1,800,475],'name','rlclr');
plot(dt,rlclr,'k--',dt,rl_best,'r-','LineWidth',1);
ylabel('IR Flux (W/m^2)'); grid on; grid minor;
legend('clear sky','pir','location','best');
title('Downwelling Longwave Radiative Flux');
if savefigs
    ppath = fullfile(png_path,['030_',cruise,'_IR_Best_Clearsky.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['030_',cruise,'_IR_Best_Clearsky.fig']));
end

%% run COARE 3.5 for 10-min bulk model estimates
%{
        input format:
        A = coare35vn(u,zu,t,zt,rh,zq,P,ts,Rs,Rl,lat,zi,rain,cp,sigH,jcool)

        cp, sigH set to NaN
        zi=600

        output: 6 x 21 array: conforms to legacy output format
         1      usr     friction velocity that includes gustiness (m/s)
         2      tau     wind stress (N/m^2)
         3      hsb     sensible heat flux into ocean (W/m^2)
         4      hlb     latent heat flux into ocean (W/m^2)
         5      hlwebb  Webb correction for latent heat flux, add this to directly
                        measured eddy covariance latent heat flux from Licor
         6      tsr     temperature scaling parameter (K)
         7      qsr     specific humidity scaling parameter (g/Kg)
         8      zo      roughness length (m)
         9      zot     thermal roughness length (m)
         10     zoq     moisture roughness length (m)
         11     Cd      wind stress transfer (drag) coefficient at height zu
         12     Ch      sensible heat transfer coefficient (Stanton number) at ht zu
         13     Ce      latent heat transfer coefficient (Dalton number) at ht zq
         14     L       Obukhov length scale (m)
         15     zet     Monin-Obukhov stability parameter zu/L
         16     dter    cool-skin temperature depression (degC)
         17     dqer    cool-skin humidity depression (degC)
         18     tkt     cool-skin thickness (m)
         19     RF      rain heat flux (W/m2)
         20     Cdn_10  neutral value of drag coefficient at 10m
         21     Chn_10  neutral value of Stanton number at 10m
         22     Cen_10  neutral value of Dalton number at 10m
         23     U10N    10m neutral wind speed, m/s
         24     T10N    10m neutral air temperature, degC
         25     Q10N    10m neutral air humidity, g/kg
         26     wbar    mean vertical velocity, m/s
         27     dt_wrm  warming across entire warm layer degC
         28     tk_pwp  warm layer thickness m
         29     dsea    heating at depth of tsg sst
%}

if jwarm && ~all(isnan(da.sst_ship)) % run warm layer version of C35 with Ship SST first
    % need lon in +/- 180 format for warm layer...
    Lon2 = da.lon_psl; Lon2(Lon2>180) = Lon2(Lon2>180) -360;
    temp = coare35vnWarm(da.jd,wspd_best,zu,ta_best,zt,rh_best,zq,P_best,...
        da.sst_ship,rs_best,rl_best,da.lat_psl,Lon2,zi,rain_best,ztsg,zref,jcool);
    wrm = temp(:,end-2:end);
    % now run c35 again with 'best' sst
    c35 = coare35vn(wspd_best,zu,ta_best,zt,rh_best,zq,P_best,ts_best,rs_best,...
        rl_best,da.lat_psl,zi,rain_best,NaN,NaN,zref,jcool);
    % add warm layer variables to 'best' output
    c35 = [c35, wrm];
else
    % just run with 'best' sst
    c35 = coare35vn(wspd_best,zu,ta_best,zt,rh_best,zq,P_best,ts_best,rs_best,...
        rl_best,da.lat_psl,zi,rain_best,NaN,NaN,zref,jcool);
    temp = array2table(NaN(height(da),3),'VariableNames',{'dt_wrm','tk_pwp','dsea'});
    % add NaN warm layer variables
    c35 = [c35, temp];
end

% update c35 variables in da table
da.hsb_psl = c35.hsb;        % sensible heat flux W/m2
da.hlb_psl = c35.hlb;        % latent heat flux W/m2
da.hlwebb_psl = c35.hlwebb;  % hl webb correction
da.taub_psl = c35.tau;       % stress N/m2
da.zo_psl = c35.zo;          % vel roughness length
da.zot_psl = c35.zot;        % temp "
da.zoq_psl = c35.zoq;        % hum  "
da.L_psl = c35.L;            % Obukhov Length, m
da.usb_psl = c35.usr;        % ustar
da.tsb_psl = c35.tsr;        % tstar
da.qsb_psl = c35.qsr;        % qstar (g/kg)
da.zet_psl = c35.zet;        % z/L
da.dter_psl = c35.dter;      % cool skin delta T
da.dqer_psl = c35.dqer;      %  "   "      "   q
da.tkt_psl = c35.tkt;        % thickness of cool skin
da.RF_psl = c35.RF;          % rain heat flux
da.Cd_psl = c35.Cd;          % drag @ zu
da.Ch_psl = c35.Ch;          % Stanton
da.Ce_psl = c35.Ce;          % Dalton
da.Cdn_10_psl = c35.Cdn_10;  % neutral drag @ 10m [includes gustiness]
da.Chn_10_psl = c35.Chn_10;  %
da.Cen_10_psl = c35.Cen_10;  %

% could add additional variables to C35 output...
% ug      =    c35.;      % gustiness
% Urf     =    c35.;      % U,Ta,Qa & RH at user ref height
% Trf     =    c35.;
% Qrf     =    c35.;
% RHrf    =    c35.;
% dtw     =    c35.;      % warm layer T diff
% zw      =    c35.;      % warm layer thickness

% latent heat vap for seawater
% using default salinity of 33...
[~,Le] = Le_water(ts_best, ones(size(ts_best))*33); % Le in J/g
Le = Le*1000; % J/kg

%% ********* define filter criteria for 10-min turbulence results ***********

tsbs = da.tsb_psl + 0.51*(ta_best+C2K).*da.qsb_psl/1000; 
wtv = -da.usb_psl.*(da.tsb_psl + 0.61.*(ta_best+C2K).*da.qsb_psl/1000);  % w'Tsonic'
ws = (9.83./(ta_best+C2K).*abs(wtv)*600).^0.333;  % convective velocity
ug = sqrt(da.usb_psl.^2 + ws.^2);
ugw = sqrt(1.5*da.usb_psl.^2 + ws.^2);
ugu = sqrt(8*da.usb_psl.^2 + 1.*ws.^2);
gg = sqrt(1+1.6*ws.^2./wspd_best.^2);

% filter with basic turbulence criteria
ttt = (rwdir_best<rdir_hi & rwdir_best>rdir_lo) &...
       std_hed_best<sigh_lim & std_sog_best<sigu_lim & rain_best<rain_lim & ...
       da.vplat_std<vplat_lim & (sqrt(da.wvar)./ugw < 2+0.0015*wspd_best.^2) & ...
       (sqrt(da.vvar)./ugu < 1.4);

%% 10min covariance stress & heat fluxes
% streamwise stress = rhoair * w'u'
tauc = -(da.rhoa_psl.*da.wu_csp); % NOAA sonic, csp

% compute hsc again - u* and q* have been updated by new COARE run
% first, correct w'Tsonic' for humidity contribution
da.wT_csp = da.wTs_csp - 0.51*(ta_best+C2K).*-da.usb_psl.*da.qsb_psl/1000;

% then convert to heat flux units
% hsc = w't' * rhoair * Cp_air = K*m/s * kg/m3 * J/kg*K = W/m2
hsc = da.wT_csp*cpa.*da.rhoa_psl; % NOAA sonic

% compute hlc from licor w'q', no Webb correction here
hlc_lic = da.wq_csp/1000.*Le.*da.rhoa_psl;  % kg/kg m/s * J/kg * kg/m3 = J/s/m2 = W/m2

%% prepare heat flux QC plots
% plot of w't'/U10N vs deltaT should be linear, slope is ~Ch10n, and
% x-intercept is measure of bias in deltaT measurement
grav = grv(nanmean1(da.lat_psl));
lapse = grav/cpa;                   % lapse rate for potential temp calc
Theta10 = c35.T10N + lapse*10;      % Pot temp @ 10m
qs_adj = qs_best; %- da.dqer_psl*1000;   % qs - dqer, no adjust for cool skin
ts_adj = ts_best; %- da.dter_psl;        % ts - dter, no adjust for cool skin

wt_U = da.wT_csp./c35.U10N;
dTs = ts_adj - Theta10;  % Delta T, use potential temp at 10m
% wT plot criteria
idx = (isfinite(wt_U) & isfinite(dTs) & dTs<3 & dTs>-2 & wt_U>-0.01 & wt_U<0.01);
minx = min(dTs(idx)); maxx = max(dTs(idx));
coef = polyfit(dTs(idx),wt_U(idx),1);
liney = coef(1)*[minx,maxx] + coef(2); xIntcpt = coef(2)/-coef(1);

figure('position',[1,1,800,475],'name','wT_v_dT');
plot(dTs(idx),wt_U(idx),'r.',[minx,maxx],liney,'b','MarkerSize',8); ylim([-0.01,0.01]);
hold on; plot([0,0],[-0.01,0.01],'k',[minx,maxx],[0,0],'k'); xlim([minx,maxx]);
ylabel('wT/U_{10n}'); xlabel('(Ts - \Theta_{10})'); grid on; grid minor;
title(['Ch10N x 1000 = ',num2str(coef(1)*1000,'%6.5f'),'   y-Intcpt = ',...
    num2str(coef(2),'%6.5f'),'   x-Intcpt = ',num2str(xIntcpt,'%6.5f')],'FontSize',16);
if savefigs
    ppath = fullfile(png_path,['031_',cruise,'_wt_U_vs_dTs.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['031_',cruise,'_wt_U_vs_dTs.fig']));
end

wq_U = (da.wq_csp - (c35.wbar.*qa_best))./c35.U10N; % (wq - wbar*qbar)/u10n - Webb correction
dQs = qs_adj - c35.Q10N;
idx = (isfinite(wq_U) & isfinite(dQs) & wq_U>-0.005 & wq_U<0.02); % wq plot criteria
minx = min(dQs(idx))-2; maxx = max(dQs(idx))+2;
coef = polyfit(dQs(idx),wq_U(idx),1);
liney = coef(1)*[minx,maxx] + coef(2); xIntcpt = coef(2)/-coef(1);

figure('position',[1,1,800,475],'name','wq_v_dQ');
plot(dQs(idx),wq_U(idx),'r.',[minx,maxx],liney,'b','MarkerSize',8); ylim([-0.005,0.005]);
hold on; plot([0,0],[-0.02,0.02],'k',[minx,maxx],[0,0],'k'); xlim([-4,8]);
ylabel('wq/U_{10n}'); xlabel('(Qs - Q_{10N})'); grid on; grid minor;
title(['Ce10N x 1000 = ',num2str(coef(1)*1000,'%6.5f'),'   y-Intcpt = ',...
    num2str(coef(2),'%6.5f'),'   x-Intcpt = ',num2str(xIntcpt,'%6.5f')],'FontSize',16);
if savefigs
    ppath = fullfile(png_path,['032_',cruise,'_wq_U_vs_dQs.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['032_',cruise,'_wq_U_vs_dQs.fig']));
end

%% inertial dissipation processing
zm = -100:.1:20;
cuxx = psi_fu(zm);
ctxx = psi_ft(zm);
%figure; semilogy(zm,cuxx,zm,ctxx);

% ship plume contam index based on wind variances, 1=bad, 0=good
jplume = ones(length(da.dnum),1);
ii = (sqrt(da.wvar)./ugw < 2 + 0.0015*wspd_best.^2) & (sqrt(da.vvar)./ugu < 2);
jplume(ii) = 0;
jplume(isnan(jplume)) = 1;

% ship maneuver index, <3 good
jmanuv = 3*ones(length(da.dnum),1);
ii = (std_hed_best<sigh_lim) & (std_sog_best<sigu_lim) & (da.vplat_std<vplat_lim);
jmanuv(ii) = std_sog_best(ii);
jmanuv(isnan(jmanuv)) = 3;

% examine wind speed structure function parameter z/L dependence
zzl = [-60 -30 -20 -8 -5 -4 -3 -2 -1 -.7 -.5 -.3 -.2 -.1 -.05 0 .05 .1 1 2 3 4 5 8 20];
nn = length(zzl)-1;
zetxm = NaN(nn,1); csm = zetxm; cnt = zetxm;
for i=1:nn
   ii = find(ttt & da.zet_psl>zzl(i) & da.zet_psl<zzl(i+1) & isfinite(cs));
   if ~isempty(ii)
      zetxm(i) = nanmedian1(da.zet_psl(ii));
      csm(i) = nanmedian1(cs(ii)./da.usb_psl(ii).^2*zu^0.667);
      cnt(i) = length(ii);
   else
      zetxm(i) = NaN;
      csm(i) = NaN;
      cnt(i) = NaN;
   end
end

figure('position',[1,1,800,475],'name','cs_psi_fu');
semilogy(zm,cuxx,'k-'); hold on;
semilogy(zetxm,csm,'rO','LineWidth',1);
grid on; grid minor; axis([-20 20 1 100]);
xlabel('z/L'); ylabel('Psi-fu');
legend('Psi-fu (z/L)','NOAA','location','best');
title([ptitle,'  Cs^2/usr * zu^{2/3}, Psi-fu (z/L)']);
if savefigs
    ppath = fullfile(png_path,['033_',cruise,'_Psi_fu.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['033_',cruise,'_Psi_fu.fig']));
end

figure('position',[1,1,800,475],'name','cs_psi_fu_ratio');
plot(zm,ones(length(zm)),'k-','LineWidth',1); hold on;
plot(zetxm,csm./psi_fu(zetxm),'bd','LineWidth',1); hold off; grid on; grid minor;
xlim([min(zetxm)-5,max(zetxm)+5]); xlabel('z/L'); ylabel('Ratio Csm/Psi-fu');
title([ptitle,'  Csm/Psi-fu Ratio']);
if savefigs
    ppath = fullfile(png_path,['034_',cruise,'_Csm_Psi_fu_ratio.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['034_',cruise,'_Csm_Psi_fu_ratio.fig']));
end

% compute ID usr,tsr,qsr
usib = sqrt(cs.*zu^0.667./psi_fu(da.zet_psl)); usib = real(usib);
% tsibs = -sqrt(ct.*zu^0.667./psi_ft(zetx)).*sign(wtv); tsibs = real(tsibs);
DT = ts_best - ta_best + da.dter_psl;  % using delta T rather than wTsonic for sign
tsibs = -sqrt(ct.*zu^0.667./psi_ft(da.zet_psl)).*sign(DT); tsibs = real(tsibs);

qsib_lic = sqrt(cql.*zq^0.667./psi_ft(da.zet_psl)).*sign(da.qsb_psl)/1000; % convert to kg/kg
qsib_lic = real(qsib_lic);

% update filter for ID stress
sss = ttt & (abs(-da.wu_csp - usib.^2) < 0.02 + 0.006*wspd_best);

% Plot ID estimate for 10-min Cd10N vs U
% filter out cdn outliers
mmm = ttt & (usib./c35.U10N).^2 < 0.004;
figure('position',[1,1,800,475],'name','id_cd10n_v_u');
plot(c35.U10N(mmm),c35.Cdn_10(mmm)/1000,'k.',c35.U10N(mmm),...
     (usib(mmm)./c35.U10N(mmm)).^2,'r.','LineWidth',1,'MarkerSize',8);
axis([0 14 4e-4 2e-3]); xlabel('U10n m/s'); ylabel('ID Cd_{10n}');
legend('COARE 3.5','NOAA','location','best'); grid on; grid minor;
title([ptitle,' ID Cd10n, 10-min']);
if savefigs
    ppath = fullfile(png_path,['035_',cruise,'_ID_Cd10N_v_U.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['035_',cruise,'_ID_Cd10N_v_U.fig']));
end

% binned Cd10N vs U
[Ubins, Cdbin, Cdbinm, SD_Cdbin, nn] = binave2(c35.U10N(mmm),(usib(mmm)./c35.U10N(mmm)).^2,2,1,15);
cc = isfinite(Ubins);
figure('position',[1,1,800,475],'name','id_cd10n_v_u_bin');
p0 = plot(c35.U10N(ttt),c35.Cdn_10(ttt)/1000,'k.','MarkerSize',8); hold on;
ii = nn>0;
p1 = errorbar(Ubins(ii),Cdbin(ii),SD_Cdbin(ii)./sqrt(nn(ii)),'rd-',...
              'LineWidth',1,'MarkerSize',8); grid on; grid minor;
axis([0 14 4e-4 2e-3]); xlabel('U_{10n} m/s'); ylabel('Median ID Cd_{10n}');
legend([p0,p1],{'COARE 3.5','NOAA'},'location','best');
title([ptitle,' Binned ID Cd10n, Median +/- SEM']);
if savefigs
    ppath = fullfile(png_path,['036_',cruise,'_bin_ID_Cd10N_v_U.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['036_',cruise,'_bin_ID_Cd10N_v_U.fig']));
end

% ID sensible heat fluxes
tsib = tsibs - 0.51*(ta_best + C2K).*da.qsb_psl/1000; % moisture corrrected tsr
hsib = -da.rhoa_psl*cpa.*(tsib).*usib;

minx = max(min(da.hsb_psl(ttt)),-40)-10; maxx = max(da.hsb_psl(ttt))+10;
figure('position',[1,1,800,475],'name','hsib_v_hsb');
plot([minx maxx],[minx maxx],'k--'); hold on;
plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k--');
plot(da.hsb_psl(ttt),hsib(ttt),'r.','LineWidth',1,'MarkerSize',8); grid on; grid minor;
axis([minx maxx minx maxx]); xlabel('Hs COARE 3.5 (W/m^2))');
ylabel('Hs ID (W/m^2))'); title([ptitle,' ID Sensible Heat Flux, 10-min']);
if savefigs
    ppath = fullfile(png_path,['037_',cruise,'_ID_Hs_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['037_',cruise,'_ID_Hs_v_bulk.fig']));
end

% binned ID Hs
[hsbbin, hsib_bin, hsib_binm, SD_hsib, ~] = binave2(da.hsb_psl(ttt),hsib(ttt),5,-30,30);
minx = min(hsbbin)-10; maxx = max(hsbbin)+10;
figure('position',[1,1,800,475],'name','hsib_v_hsb_bin');
p0 = plot([minx maxx],[minx maxx],'k--'); hold on; grid on; grid minor;
p1 = plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k--');
p2 = errorbar(hsbbin,hsib_bin,SD_hsib,'rd-','LineWidth',1,'MarkerSize',8);
axis([minx maxx minx maxx]); xlabel('Hs COARE 3.5 (W/m^2))');
ylabel(' Mean Hs ID (W/m^2))'); title([ptitle,' binned ID Sensible Heat Flux, Mean +/- SD']);
if savefigs
    ppath = fullfile(png_path,['038_',cruise,'_bin_ID_Hs_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['038_',cruise,'_bin_ID_Hs_v_bulk.fig']));
end

% compare hsib and hsc
minx = max(min(da.hsb_psl(ttt)),-20)-10; maxx = max(da.hsb_psl(ttt))+10;
figure('position',[1,1,800,475],'name','hsib_hsb_v_hsb_');
p1 = plot([minx maxx],[minx maxx],'k--'); hold on;
p2 = plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k--');
p3 = plot(da.hsb_psl(ttt),hsib(ttt),'b.',da.hsb_psl(ttt),hsc(ttt),'r.','LineWidth',1,'MarkerSize',8);
axis([minx maxx minx maxx]); grid on; grid minor; xlabel('Hs COARE 3.5 (W/m^2))');
ylabel('Hs (W/m^2))'); title([ptitle,' NOAA hsib & hsc, 10-min']);
legend(p3,{'hsib','hsc'},'location','best');
if savefigs
    ppath = fullfile(png_path,['039_',cruise,'_hsib_hsc_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['039_',cruise,'_hsib_hsc_v_bulk.fig']));
end

% ID latent heat flux
hlib_lic = -da.rhoa_psl.*Le.*(qsib_lic).*usib; % kg/m3 * J/kg * kg/kg * m/s = J/s/m2 = W/m2
hlib_lic(hlib_lic>1000) = NaN;  hlib_lic(hlib_lic<-1000) = NaN;% kill bad id hl values if present
%ii = find(sgq_lic<5e-3 | sgq_lic>.7 | hlib_lic<-200); hlib_lic(ii)=NaN;
minx = min(hlib_lic(ttt))-10; maxx = max(hlib_lic(ttt))+10;
figure('position',[1,1,800,475],'name','hlib_v_hlb');
plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k--'); hold on;
plot(da.hlb_psl(ttt),hlib_lic(ttt) + da.hlwebb_psl(ttt),'.','MarkerSize',8);
axis([-60 60 -60 60]);  grid on; grid minor;
xlabel('Bulk Hl (W/m^2)'); ylabel('Hlib + webb (W/m^2)');
title([ptitle,' ID Latent Heat Flux (Webb Corrected), 10-min']);
if savefigs
    ppath = fullfile(png_path,['040_',cruise,'_ID_Hl.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['040_',cruise,'_ID_Hl.fig']));
end

% plot covariance latent heat
minx = min(hlc_lic(ttt))-10; maxx = max(hlc_lic(ttt))+10;
figure('position',[1,1,800,475],'name','hlc_v_hlb');
plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k--'); hold on;
plot(da.hlb_psl(ttt),hlc_lic(ttt) + da.hlwebb_psl(ttt),'.'); axis([minx maxx minx maxx]);
xlabel('Bulk Hl (W/m^2)'); ylabel('Hlc + webb (W/m^2)'); grid on; grid minor;
title([ptitle,' Covariance Latent Heat Flux (Webb Corrected), 10-min']);
if savefigs
    ppath = fullfile(png_path,['041_',cruise,'_Cov_Hl.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['041_',cruise,'_Cov_Hl.fig']));
end

% using ID Hl as best latent heat measurement...
hlm = NaN(length(hlib_lic),1);
hlm(ttt) = hlib_lic(ttt);
hlm(hlm>1000) = NaN;

% wind speed dependence of Hl
uu = 0:20; % wind speed bins
hlmm = NaN(length(uu),1);  hlbm = hlmm; dqm = hlmm; nn = hlmm;
for i=1:length(uu)
    ii = find(ttt & floor(c35.U10N)==uu(i));
    if isempty(ii)
        hlmm(i) = NaN;
        hlbm(i) = NaN;
        nn(i) = 0;
        dqm(i) = NaN;
    else
        hlmm(i) = nanmedian1(hlm(ii) + da.hlwebb_psl(ii));
        hlbm(i) = nanmedian1(da.hlb_psl(ii));
        nn(i) = length(ii);
        dqm(i) = nanmedian1(qs_best(ii)-qa_best(ii))-1e3*nanmean1(da.dqer_psl(ii));
    end
end

miny = min(hlmm)-10; maxy = max(hlmm)+10;
figure('position',[1,1,800,475],'name','hl_ubin');
plot(uu+0.5,hlbm,'go-',uu+0.5,hlmm,'bo-','LineWidth',1,'MarkerSize',8);
ylabel('Hlm + webb (W/m^2)');  xlabel('U_{10N} (m/s)'); grid on; grid minor;
if isfinite(miny); ylim([miny maxy]); end
legend('COARE 3.5 median','Hlm + webb, median','location','northwest');
title([ptitle,' 10-min Wind Binned Latent Heat Flux']);
if savefigs
    ppath = fullfile(png_path,['042_',cruise,'_Hlm_WindBin.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['042_',cruise,'_Hlm_WindBin.fig']));
end

figure('position',[1,1,800,475],'name','hlc_v_hlb_u_bin');
plot(hlbm,hlmm,'o',[miny maxy],[miny maxy],'LineWidth',1,'MarkerSize',8);
ylabel('Hlm + webb (W/m^2)');  xlabel('Bulk Hl (W/m^2)'); grid on; grid minor;
if isfinite(maxy); axis([miny maxy miny maxy]); end
title([ptitle,' 10-min Wind Binned Hlib vs COARE 3.5']);
if savefigs
    ppath = fullfile(png_path,['043_',cruise,'_Hlm_Bin_v_C35.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['043_',cruise,'_Hlm_Bin_v_C35.fig']));
end

% filtered ID stress - NOAA sonic
minx = min(-da.wu_csp(ttt))-0.02; maxx = max(-da.wu_csp(ttt))+0.02;
figure('position',[1,1,800,475],'name','wu_filt');
p1 = plot([stdt endt],[0 0]); hold on;
p2 = plot(dt(ttt),-da.wu_csp(ttt),'.',dt(ttt),usib(ttt).^2,'.','MarkerSize',8);
xlim([stdt endt]); ylim([minx maxx]); ylabel('-wu & usib^2 (m^2/s^2)');
legend(p2,{'-wu','usib^2'},'location','northwest'); grid on; grid minor;
title([ptitle,'  WU Std Filtering: ID and Cov, NOAA sonic']);
if savefigs
    ppath = fullfile(png_path,['044_',cruise,'_WU_StdFiltered.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['044_',cruise,'_WU_StdFiltered.fig']));
end

minx = min(-da.wu_csp(sss))-0.02; maxx = max(-da.wu_csp(sss))+0.02;
figure('position',[1,1,800,475],'name','wu_filt2');
p1 = plot([stdt endt],[0 0]); hold on;
p2 = plot(dt(sss),-da.wu_csp(sss),'.',dt(sss),usib(sss).^2,'.','MarkerSize',8);
xlim([stdt endt]); ylim([minx maxx]); xlabel('Day-Of-year'); ylabel('-wu & usib^2 (m^2/s^2)');
legend(p2,{'-wu','usib^2'},'location','northwest'); grid on; grid minor;
title([ptitle,'  WU Additional Filtering: ID and Cov, NOAA sonic']);
if savefigs
    ppath = fullfile(png_path,['045_',cruise,'_WU_AddFiltered.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['045_',cruise,'_WU_AddFiltered.fig']));
end

% Other ID QC plots
figure('position',[1,1,800,475],'name','cq2');
semilogy(da.zet_psl(ttt),cql(ttt)./(da.qsb_psl(ttt)).^2.*zu^0.67,'.',...
    da.zet_psl(ttt),psi_ft(da.zet_psl(ttt)),'.','MarkerSize',8);
axis([-2 2 0.1 5e2])
xlabel('z/L'); ylabel('Cq^2/q_*^2/z^{2/3} Licor'); grid on; grid minor;
title([ptitle,'  Cq^2 and Psi-ft(z/L)']);
if savefigs
    ppath = fullfile(png_path,['046_',cruise,'_cq2_lic.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['046_',cruise,'_cq2_lic.fig']));
end

figure('position',[1,1,800,475],'name','ct2');
semilogy(da.zet_psl(ttt),psi_ft(da.zet_psl(ttt)),'.',...
    da.zet_psl(ttt),ct(ttt)./(tsbs(ttt)).^2.*zu^0.67,'.','MarkerSize',8);
axis([-2,2,0.1,5e2]); xlabel('z/L'); ylabel('Ct^2/t_*^2/z^{2/3}');
legend('Psi-ft(z/L)','NOAA','location','northwest');
title([ptitle,'  Ct^2 and Psi-ft(z/L)']); grid on; grid minor;
if savefigs
    ppath = fullfile(png_path,['047_',cruise,'_ct2.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['047_',cruise,'_ct2.fig']));
end

figure('position',[1,1,800,475],'name','sigQ');
semilogy(da.zet_psl(ttt),phi_tt(da.zet_psl(ttt)),'.',da.zet_psl(ttt),...
         sqrt(da.qvar(ttt))./abs(da.qsb_psl(ttt)),'.','MarkerSize',8);
axis([-2,2,0.1,5e1]); xlabel('z/L'); ylabel('Dimensionless \sigma q Licor');
title([ptitle,'  Dimensionless \sigma q Licor and Psi-tt(z/L)']); grid on; grid minor;
if savefigs
    ppath = fullfile(png_path,['048_',cruise,'_sigmaQ.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['048_',cruise,'_sigmaQ.fig']));
end

sigTs = sqrt(da.tvar);
figure('position',[1,1,800,475],'name','sigTs');
semilogy(da.zet_psl(ttt),phi_tt(da.zet_psl(ttt)),'.',...
    da.zet_psl(ttt),sigTs(ttt)./abs(tsbs(ttt)),'.','MarkerSize',8);
axis([-2,2,0.1,1e3]); xlabel('z/L'); ylabel('Dimensionless \sigma Tosnic');
legend('Phi-tt(z/L)','NOAA','location','best');
title([ptitle,'  Dimensionless \sigma Ts and Psi-tt(z/L)']); grid on; grid minor;
if savefigs
    ppath = fullfile(png_path,['049_',cruise,'_sigmaTs.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['049_',cruise,'_sigmaTs.fig']));
end

figure('position',[1,1,800,475],'name','sigQ_v_sigTs');
semilogy([-10 2],[.85 .85],da.zet_psl(ttt),(sqrt(da.qvar(ttt))./abs(da.qsb_psl(ttt))./...
         (sigTs(ttt)./abs(tsbs(ttt)))),'.','MarkerSize',8); grid on; grid minor;
axis([-2 2 0.01 1e2]); xlabel('z/L'); ylabel('Dimensionless \sigma qLic/ \sigma t');
title([ptitle,'  \sigma q Licor / \sigma Tosnic, NOAA Sonic']);
if savefigs
    ppath = fullfile(png_path,['050_',cruise,'_sigmaQ_Ts_Ratio.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['050_',cruise,'_sigmaQ_Ts_Ratio.fig']));
end

%% compute 10min u*, t* and q* from covariance measurements
% compare cov and ID results with bulk model

usr_obs = sqrt(abs(da.wu_csp)).* -sign(da.wu_csp);
qsr_obs = -da.wq_csp./usr_obs;
eh2o = 6.112.*exp(17.67.*ta_best./(ta_best+243.5)).*rh_best/100;  % VP in mb
wt_obs = da.wTs_csp./(1+0.32*eh2o./P_best);  % humidity correction to Tsonic, Kaimal 1991
tsr_obs = -wt_obs./usr_obs;

% cov ustar
figure('position',[1,1,800,475],'name','usr_v_usb');
plot([-0.2,0.75],[0,0],'k-',[0,0],[-0.2,0.75],'k-'); hold on;
plot([-1,5],[-1,5],'k-',da.usb_psl(sss),usr_obs(sss),'.','MarkerSize',8); hold off;
grid on; grid minor; axis([-0.2,0.75,-0.2,0.75]);
xlabel('u_*, COARE 3.5, m/s'); ylabel('u_* cov, m/s');
title([ptitle,'  u_*, Cov vs. COARE']);
if savefigs
    ppath = fullfile(png_path,['051_',cruise,'_ustar_cov_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['051_',cruise,'_ustar_cov_v_bulk.fig']));
end

% binned cov ustar
[usrbin, usrc_bin, usrc_binm, SD_usrc, ~] = binave2(da.usb_psl(sss),...
                                            usr_obs(sss),0.05,-0.1,0.75);
minx = min(usrbin)-0.15; maxx = max(usrbin)+0.15;
figure('position',[1,1,800,475],'name','usr_v_usb_bin');
plot([minx maxx],[minx maxx],'k--'); hold on; grid on; grid minor;
plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k--');
errorbar(usrbin,usrc_bin,SD_usrc,'rd-','MarkerSize',8,'LineWidth',1);
axis([minx maxx minx maxx]); xlabel('u_* COARE 3.5 (m/s)');
ylabel(' Mean Cov u_* (m/s)'); title([ptitle,' Binned Cov u_*, Mean +/- SD']);
if savefigs
    ppath = fullfile(png_path,['052_',cruise,'_bin_ustar_cov_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['052_',cruise,'_bin_ustar_cov_v_bulk.fig']));
end

% ID ustar
figure('position',[1,1,800,475],'name','usib_v_usb');
plot([-1,5],[-1,5],'k-',da.usb_psl(sss),usib(sss),'.','MarkerSize',8);
xlabel('u_*, COARE 3.5, m/s'); ylabel('u_* ID, m/s'); grid on; grid minor;
title([ptitle,'  u_* ID vs. COARE']); axis([0,0.75,0,0.75]);
if savefigs
    ppath = fullfile(png_path,['053_',cruise,'_ustar_ID_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['053_',cruise,'_ustar_ID_v_bulk.fig']));
end

% binned ID ustar
[usrbin, usib_bin, usib_binm, SD_usib, ~] = binave2(da.usb_psl(sss),usib(sss),0.05,-0.1,0.5);
minx = min(usrbin)-0.1; maxx = max(usrbin)+0.1;
figure('position',[1,1,800,475],'name','usib_v_usb_bin');
plot([minx maxx],[minx maxx],'k--'); hold on; grid on; grid minor;
plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k--');
errorbar(usrbin,usib_bin,SD_usib,'rd-','LineWidth',1,'MarkerSize',8);
axis([minx maxx minx maxx]); xlabel('u_* COARE 3.5 (m/s))');
ylabel(' Mean ID u_* (m/s))'); title([ptitle,' Binned ID u_*, Mean +/- SD']);
if savefigs
    ppath = fullfile(png_path,['054_',cruise,'_bin_ustar_ID_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['054_',cruise,'_bin_ustar_ID_v_bulk.fig']));
end

% cov tstar
figure('position',[1,1,800,475],'name','tsr_v_tsb');
plot([-1,5],[-1,5],'k-',[-1,1],[0,0],'k-',[0,0],[-1,1],'k-'); hold on;
plot(da.tsb_psl(ttt),tsr_obs(ttt),'.','MarkerSize',8); hold off; axis([-0.3,0.3,-0.3,0.3]);
xlabel('t_*, COARE 3.5, °C'); ylabel('t_*, turb, °C');
title([ptitle,'  t_*, Turbulent vs. COARE']); grid on; grid minor;
if savefigs
    ppath = fullfile(png_path,['055_',cruise,'_tstar_cov_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['055_',cruise,'_tstar_cov_v_bulk.fig']));
end

% binned cov tstar
[tsrbin, tsrc_bin, tsrc_binm, SD_tsrc, ~] = binave2(da.tsb_psl(ttt),tsr_obs(ttt),0.025,-0.15,0.15);
minx = min(tsrbin)-0.05; maxx = max(tsrbin)+0.05;
figure('position',[1,1,800,475],'name','tsr_v_tsb_bin');
plot([minx maxx],[minx maxx],'k--'); hold on;
plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k--');
errorbar(tsrbin,tsrc_bin,SD_tsrc,'rd-','MarkerSize',8,'LineWidth',1);
axis([minx maxx minx maxx]); xlabel('t_* COARE 3.5 (°C))'); grid on; grid minor;
ylabel(' Mean Turb t_* (°C)'); title([ptitle,' Binned Turbulent t_*, Mean +/- SD']);
if savefigs
    ppath = fullfile(png_path,['056_',cruise,'_bin_tstar_cov_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['056_',cruise,'_bin_tstar_cov_v_bulk.fig']));
end

% ID tstar
figure('position',[1,1,800,475],'name','tsib_v_tsb');
plot([-1,5],[-1,5],'k-',[-1,1],[0,0],'k-',[0,0],[-1,1],'k-'); hold on;
plot(da.tsb_psl(ttt),tsib(ttt),'.','MarkerSize',8); axis([-0.3,0.3,-0.3,0.3]);
xlabel('t_*, COARE 3.5, °C'); ylabel('t_*, ID, °C');
title([ptitle,'  t_*, ID vs. COARE']); grid on; grid minor;
if savefigs
    ppath = fullfile(png_path,['057_',cruise,'_tstar_ID_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['057_',cruise,'_tstar_ID_v_bulk.fig']));
end

% binned ID tstar
[tsrbin, tsib_bin, tsib_binm, SD_tsib, ~] = binave2(da.tsb_psl(ttt),tsib(ttt),0.025,-0.15,0.15);
minx = min(tsrbin)-0.15; maxx = max(tsrbin)+0.05;
figure('position',[1,1,800,475],'name','tsib_v_tsb_bin');
plot([minx maxx],[minx maxx],'k--'); hold on; grid on; grid minor;
plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k--');
errorbar(tsrbin,tsib_bin,SD_tsib,'rd-','LineWidth',1,'MarkerSize',8);
axis([minx maxx minx maxx]); xlabel('t_* COARE 3.5 (°C))');
ylabel(' Mean ID t_* (°C))'); title([ptitle,' Binned ID t_*, Mean +/- SD']);
if savefigs
    ppath = fullfile(png_path,['058_',cruise,'_bin_tstar_ID_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['058_',cruise,'_bin_tstar_ID_v_bulk.fig']));
end

figure('position',[1,1,800,475],'name','qsr_v_qsb');
plot([-1,5],[-1,5],'k-',[-1,1],[0,0],'k-',[0,0],[-1,1],'k-'); hold on;
plot(da.qsb_psl(ttt),qsr_obs(ttt),'.','MarkerSize',8); axis([-0.1,0.05,-0.1,0.05]);
xlabel('q_*, COARE 3.5, g/kg'); ylabel('q_*, turb, g/kg');
title([ptitle,'  q_*, Turbulent vs. COARE, NOAA sonic']); grid on; grid minor;
if savefigs
    ppath = fullfile(png_path,['059_',cruise,'_qstar_cov_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['059_',cruise,'_qstar_cov_v_bulk.fig']));
end

figure('position',[1,1,800,475],'name','qsib_v_qsb');
plot([-1,5],[-1,5],'k-',[-1,1],[0,0],'k-',[0,0],[-1,1],'k-'); hold on;
plot(da.qsb_psl(ttt),qsib_lic(ttt)*1000,'.','MarkerSize',8); axis([-0.1,0.05,-0.1,0.05]);
xlabel('q_*, COARE 3.5, g/kg'); ylabel('q_*, ID, g/kg'); grid on; grid minor;
title([ptitle,'  q_*, ID vs. COARE, NOAA sonic']);
if savefigs
    ppath = fullfile(png_path,['060_',cruise,'_qstar_ID_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['060_',cruise,'_qstar_ID_v_bulk.fig']));
end

close all

%% ocean warm layer profile
% this only works if the warm layer bulk model is run using tsg sst...
% otherwise zw = 0 and tmxd = tsg
if jwarm
    np = length(da.dnum);
    dd = round(ztsg);  % meters depth for tsg inlet
    dtz = zeros(np,dd+1);
    zz=0:dd;
    for i=1:dd+1
        ii = find(zw>zz(i));
        dtz(ii,i) = dtw(ii).*(1-zz(i)./zw(ii));  % warming from mixed layer to z
    end
    dt_dd = dtz(:,dd+1);
    tmxd = da.sst_ship - dtz(:,dd+1);  % ML temp: tsg with with 5-m warm layer removed
end

%% 10-minute net heat flux from bulk and turbulence measurements
% net longwave heat flux - upward positive
rnl = 0.97*(5.67e-8*(ts_best+C2K).^4) - rl_best; % rl_up - rl_down

% for accurate net heat, we need to fill gaps in the rs, rnl, hs, hl time series
% first, interp rs...
rs_best_interp = fillgaps(rs_best,36);
% NaN-out time specific segments with too much missing data...
bads = datenum(2022,7,16,1,30,00):10/1440:datenum(2022,7,16,8,50,00);
bads = [bads, datenum(2022,7,17,03,10,00):10/1440:datenum(2022,7,17,10,00,00)];
bads = [bads, datenum(2022,7,17,17,00,00):10/1440:datenum(2022,7,18,14,00,00)];
bads = [bads, datenum(2022,7,19,01,40,00):10/1440:datenum(2022,7,20,08,00,00)];
bads = [bads, datenum(2022,7,23,00,00,00):10/1440:datenum(2022,7,24,02,00,00)];
bads = [bads, datenum(2022,8,1,12,20,00):10/1440:datenum(2022,8,2,00,00,00)];
badi = [];
for ii = 1:length(bads)
    idx = find(abs(da.dnum - bads(ii))<1/1440);
    badi = [badi,idx];
end
rs_best_interp(badi) = NaN;
% set a few other nighttime spikes to zero
zeros = datenum(2022,7,17,00,30,00):10/1440:datenum(2022,7,17,01,00,00);
zeros = [zeros,datenum(2022,7,20,06,00,00):10/1440:datenum(2022,7,20,08,00,00)];
zeros = [zeros,datenum(2022,7,24,03,30,00):10/1440:datenum(2022,7,24,04,30,00)];
zeroi = [];
for ii = 1:length(zeros)
    idx = find(abs(da.dnum - zeros(ii))<1/1440);
    zeroi = [zeroi,idx];
end
rs_best_interp(zeroi) = 0;
rsclr_interp = fillgaps(rsclr); rsclr_interp(badi) = NaN;
% figure; plot(dt,rs_best,'r.',dt,rs_best_interp,'b-',dt,rsclr_interp,'k-',...
%              'markersize',8); grid on; grid minor;


% now interp rl, hsb, hlb, rlclr, sst & RF
% and NaN-out large gaps
rl_best_interp = fillgaps(rl_best,36); rl_best_interp(badi) = NaN;
hsb_interp = fillgaps(da.hsb_psl); hsb_interp(badi) = NaN;
hlb_interp = fillgaps(da.hlb_psl); hlb_interp(badi) = NaN;
rlclr_interp = fillgaps(rlclr); rlclr_interp(badi) = NaN;
sst_interp = fillgaps(ts_best); sst_interp(badi) = NaN;
% and set missing RF = zero
RF_interp = da.RF_psl;
ii = isnan(da.RF_psl); RF_interp(ii) = 0;
% figure; plot(dt,rl_best,'r.',dt,rl_interp,'b-',dt,rlclr_interp,'k-',...
%              'markersize',8); grid on; grid minor;
% figure; plot(dt,da.hsb_psl,'r.',dt,hsb_interp,'b-','markersize',8);
% grid on; grid minor;
% figure; plot(dt,da.hlb_psl,'r.',dt,hlb_interp,'b-','markersize',8);
% grid on; grid minor;
% figure; plot(dt,da.RF_psl,'r.',dt,RF_interp,'b-','markersize',8);
% grid on; grid minor;
% figure; plot(dt,ts_best,'r.',dt,sst_interp,'b-','markersize',8);
% grid on; grid minor;

% compute net IR from interpolated SST and rl
rnl_interp = 0.97*(5.67e-8*(sst_interp+C2K).^4) - rl_best_interp; % rl_up - rl_down

% best net heat flux estimate, from interpolated rs/rnl, hsb, hlb
hnet = 0.955*rs_best_interp - rnl_interp - hsb_interp - hlb_interp - RF_interp;
% net heat from cov fluxes hsc & hlc
hnet_c = 0.955*rs_best - rnl - hsc - (hlc_lic + da.hlwebb_psl) - da.RF_psl;
% net heat from id fluxes hsib & hlib_lic
hnet_id = 0.955*rs_best - rnl - hsib - (hlib_lic + da.hlwebb_psl) - da.RF_psl;

%% save ascii file of 10-min flux results if desired
if saveflx
    save_flux_10_FATIMA_2022;      % save 10-min results
end

%% daily average of heat budget variables and cloud forcing analysis
% using interpolated radiation and heat fluxes, since missing data
% is a big problem for computing daily means of the heat budget
dnum_bins = floor(da.dnum(1)):ceil(da.dnum(end)); dnum_bins = dnum_bins';
rsd = interval_avg(da.dnum,rs_best_interp,dnum_bins); rsd = rsd(:,2);
rld = interval_avg(da.dnum,rl_best_interp,dnum_bins); rld = rld(:,2);
hnd = interval_avg(da.dnum,hnet,dnum_bins); hnd = hnd(:,2);  % bulk net heat
latd = interval_avg(da.dnum,da.lat_psl,dnum_bins); latd = latd(:,2);
rsclrd = interval_avg(da.dnum,rsclr_interp,dnum_bins); rsclrd = rsclrd(:,2);
rlclrd = interval_avg(da.dnum,rlclr_interp,dnum_bins);
dt_dd = datetime(rlclrd(:,1),'ConvertFrom','datenum');
rlclrd = rlclrd(:,2);

% count number of good 10-min rs & hs values per day
Nday_rs = interval_sum(da.dnum,isfinite(rs_best_interp),dnum_bins);
Nday_hs = interval_sum(da.dnum,isfinite(hsb_interp),dnum_bins);
% identify days with 100% data coverage (144 10-min values)
ii = (Nday_rs(:,2)==144 & Nday_hs(:,2)==144);

figure('name','cld_trans_coef');
plot(latd(ii),rsd(ii)./rsclrd(ii),'o','LineWidth',1,'MarkerSize',8);
axis([Latmin,Latmax,0,1.1]); set(gca,'fontsize',14); grid on; grid minor;
xlabel('Latitude (deg)'); ylabel('Cloud Trans. Coef.');
title([ptitle,' Cloud Transmission Coef: rs/rsclr']);
if savefigs
    ppath = fullfile(png_path,['061_',cruise,'_cloud_trans_coef_v_lat.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['061_',cruise,'_cloud_trans_coef_v_lat.fig']));
end

jj = ii & (latd>Latmin);
figure('name','cld_ir_flux');
plot(1-rsd(jj)./rsclrd(jj),rld(jj)-rlclrd(jj),'o','LineWidth',1,'MarkerSize',8);
axis([0,1,0,200]); set(gca,'fontsize',14); grid on; grid minor;
xlabel('1 - Tran. Coeff.'); ylabel('Measured - Clear IR Flux (W/m^2)');
title([ptitle,' IR Flux from Cloud']);
if savefigs
    ppath = fullfile(png_path,['062_',cruise,'_cloud_IR_v_lat.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['062_',cruise,'_cloud_IR_v_lat.fig']));
end

figure('name','cld_ir_v_solar');
plot([0 300],[0 200],'k-',rsclrd(jj)-rsd(jj),rld(jj)-rlclrd(jj),'o',...
     'LineWidth',1,'MarkerSize',8);
axis([0 300 0 200]); set(gca,'fontsize',14); grid on; grid minor;
xlabel('SFC Solar Flux (W/m^2)'); ylabel('SFC IR Flux (W/m^2)');
title([ptitle,' Surface IR vs Solar Flux']);
if savefigs
    ppath = fullfile(png_path,['063_',cruise,'_IR_v_solar_flux.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['063_',cruise,'_IR_v_solar_flux.fig']));
end

figure('name','cld_forcing_v_lat'); 
plot(latd(ii),rsd(ii)-rsclrd(ii),'d',latd(ii),rld(ii)-rlclrd(ii),'x',[Latmin Latmax],[0 0 ],...
     'LineWidth',1,'MarkerSize',8);
xlim([Latmin Latmax]); set(gca,'fontsize',14); grid on; grid minor;
xlabel('Latitude (deg)'); ylabel('Cloud Forcing (W/m^2)');
legend('shortwave', 'longwave','location','best');
title([ptitle,' Cloud Forcing vs Latitude']);
if savefigs
    ppath = fullfile(png_path,['064_',cruise,'_cloud_forcing_v_lat.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['064_',cruise,'_cloud_forcing_v_lat.fig']));
end

figure('name','cld_forcing');
plot(dt_dd(ii)+hours(12),rsd(ii)-rsclrd(ii),'o-',dt_dd(ii)+hours(12),rld(ii)-rlclrd(ii),...
     'x-',[dt_dd(1),dt_dd(end) + days(1)],[0 0 ],'k-','LineWidth',1,'MarkerSize',8);
xlim([stdt endt]); set(gca,'fontsize',14); grid on; grid minor;
ylabel('Cloud Forcing (W/m^2)');
legend('shortwave', 'longwave','location','best');
title([ptitle,' Cloud Forcing']);
if savefigs
    ppath = fullfile(png_path,['065_',cruise,'_cloud_forcing.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['065_',cruise,'_cloud_forcing.fig']));
end

figure('name','daily_hnet');
plot(dt_dd(ii)+hours(12), hnd(ii),'bo-',[stdt,endt],[0,0],'k-','LineWidth',1,'MarkerSize',8);
xlim([stdt endt]); ylim([0,350]);
set(gca,'fontsize',14); grid on; grid minor;
ylabel('Daily Avg Net Heat Flux (W/m^2)');
title([ptitle,' Daily Net Heat Flux']);
if savefigs
    ppath = fullfile(png_path,['066_',cruise,'_daily_net_heat.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['066_',cruise,'_daily_net_heat.fig']));
end

close all

%% now compute 1-hr bulk averages & rerun COARE

% hourly timestamp and bins
dnum_hr_bin = floor(da.dnum(1)):1/24:ceil(da.dnum(end));
dnum_hr_bin = dnum_hr_bin(:);
dnum_hr_bin = floor(da.dnum(1)):1/24:ceil(da.dnum(end)); dnum_hr_bin = dnum_hr_bin';
dnum_hr = dnum_hr_bin(1:end-1);
dt_hr = datetime(dnum_hr,'ConvertFrom','datenum');

% first, square all raw standard deviations
da.wdir_ship_std = da.wdir_ship_std.^2;    % std wdir, ship
da.rwspd_ship_std = da.rwspd_ship_std.^2;  % std rwspd, psl
da.rwdir_ship_std = da.rwdir_ship_std.^2;  % std rwdir, psl
da.rwspd_psl_std = da.rwspd_psl_std .^2;   % std rwspd, psl
da.rwdir_psl_std = da.rwdir_psl_std.^2;    % std rwdir, psl
da.sog_ship_std = da.sog_ship_std .^2;     % std SOG, ship
da.cog_ship_std = da.cog_ship_std.^2;      % std COG, ship
da.sog_psl_std = da.sog_psl_std .^2;       % std SOG, psl
da.cog_psl_std = da.cog_psl_std.^2;        % std COG, ship
da.hed_psl_std = da.hed_psl_std.^2;        % std heading, psl
da.lrg_std = da.lrg_std.^2;                % std heading, ship
da.uplat_std = da.uplat_std.^2;            % std, uplat_std
da.vplat_std = da.vplat_std.^2;            % std, vplat_std        
da.wplat_std = da.wplat_std.^2;            % std, wplat_std        
da.odec_std = da.odec_std.^2;              % std da.odec
da.h2o_std_all = da.h2o_std_all.^2;        % std licor vars
da.h2o_std_good = da.h2o_std_good.^2; 
da.co2_std_all = da.co2_std_all.^2; 
da.co2_std_good = da.co2_std_good.^2; 

% da_hr = straight hourly average of all data in da
da_hr = interval_avg(da.dnum,da{:,:},dnum_hr_bin);

% insert date/time columns and convert to a table
da_hr(:,3) = da_hr(:,1); da_hr(:,2) = da_hr(:,1) - datenum(da_hr(1,4)-1,12,31);
da_hr = da_hr(:,2:end); da_hr(:,3:8) = datevec(da_hr(:,2));
da_hr = array2table(da_hr,'VariableNames',da_vars);
da_hr.Properties.VariableUnits = da_units;
da_hr.Properties.VariableDescriptions = da_descript;

% restore hourly standard deviations
da_hr.wdir_ship_std = sqrt(da_hr.wdir_ship_std);    % std wdir, ship
da_hr.rwspd_ship_std = sqrt(da_hr.rwspd_ship_std);  % std rwspd, psl
da_hr.rwdir_ship_std = sqrt(da_hr.rwdir_ship_std);  % std rwdir, psl
da_hr.rwspd_psl_std = sqrt(da_hr.rwspd_psl_std);    % std rwspd, psl
da_hr.rwdir_psl_std = sqrt(da_hr.rwdir_psl_std);    % std rwdir, psl
da_hr.sog_ship_std = sqrt(da_hr.sog_ship_std);      % std SOG, ship
da_hr.cog_ship_std = sqrt(da_hr.cog_ship_std);      % std COG, ship
da_hr.sog_psl_std = sqrt(da_hr.sog_psl_std );       % std SOG, psl
da_hr.cog_psl_std = sqrt(da_hr.cog_psl_std);        % std COG, ship
da_hr.hed_psl_std = sqrt(da_hr.hed_psl_std);        % std heading, psl
da_hr.lrg_std = sqrt(da_hr.lrg_std);                % std heading, ship
da_hr.uplat_std = sqrt(da_hr.uplat_std);            % std, uplat_std
da_hr.vplat_std = sqrt(da_hr.vplat_std);            % std, vplat_std        
da_hr.wplat_std = sqrt(da_hr.wplat_std);            % std, wplat_std        
da_hr.odec_std = sqrt(da_hr.odec_std);              % std da_hr.odec
da_hr.h2o_std_all = sqrt(da_hr.h2o_std_all);        % std licor vars
da_hr.h2o_std_good = sqrt(da_hr.h2o_std_good); 
da_hr.co2_std_all = sqrt(da_hr.co2_std_all); 
da_hr.co2_std_good = sqrt(da_hr.co2_std_good); 

% build 'best' table = straight hourly average of 'best' bulk met data
temp = [wspd_best,ts_best,qs_best,ta_best,rh_best,qa_best,rs_best_interp,...
        rl_best_interp,P_best,rain_best];
best = interval_avg(da.dnum,temp,dnum_hr_bin); best = best(:,2:end);
best = array2table(best,'VariableNames',{'wspd','ts','qs','ta','rh',...
                                         'qa','rs','rl','Pmb','rain'});
% hourly mean of QC variables
temp2 = [jplume,jmanuv];
temp = interval_avg(da.dnum,temp2,dnum_hr_bin);
qc_hr = array2table(temp,'VariableNames',{'dnum','jplume','jmanuv'});
qc_hr.jplume(isnan(qc_hr.jplume)) = 1;
qc_hr.jmanuv(isnan(qc_hr.jmanuv)) = 3;

% hourly mean of interpolated rs/rl clear sky models
temp = interval_avg(da.dnum,[rsclr_interp,rlclr_interp],dnum_hr_bin);
yrsclr = temp(:,2);
yrlclr = temp(:,3);

% vector average winds and compass angles:
rdir2 = rwdir_best; rdir2(rdir2<0) = rdir2(rdir2<0) + 360;
[~,rwdir_hr,rwspd_hr] = interval_avg_vectors(da.dnum, rdir2, da.rwspd_psl, dnum_hr_bin);
rwdir_hr(rwdir_hr>180) = rwdir_hr(rwdir_hr>180) + 360;
[~,wdir_hr,wspd_hr] = interval_avg_vectors(da.dnum, wdir_best, wspd_best, dnum_hr_bin);
[~,wdir_ship_hr,wspd_ship_hr] = interval_avg_vectors(da.dnum, da.wdir_ship, da.wspd_ship, dnum_hr_bin);
[~,wdir_wxt_hr,wspd_wxt_hr] = interval_avg_vectors(da.dnum, da.wdir_wxt, da.wspd_wxt, dnum_hr_bin);
if all(isnan(da.odec))
    [~,hed_hr,~] = interval_avg_vectors(da.dnum, hed_best, ones(size(hed_best)), dnum_hr_bin);
else
    [~,hed_hr,odec_hr] = interval_avg_vectors(da.dnum, hed_best, da.odec, dnum_hr_bin);
end
[~,cog_hr,sog_hr] = interval_avg_vectors(da.dnum, cog_best, sog_best, dnum_hr_bin);

%% recomputing COARE 3.5 with 'best' hourly data

% input format: c35_hr = coare35vn(u,zu,ta,zt,rh,zq,P,ts,Rs,Rl,lat,zi,rain,cp,sigH,jcool)

c35_hr = coare35vn(best.wspd,zu,best.ta,zt,best.rh,zq,best.Pmb,best.ts,best.rs,best.rl,...
    da_hr.lat_psl,zi,best.rain,NaN,NaN,10,jcool);

ywtvb = -c35_hr.usr.*(c35_hr.tsr + 0.61.*(best.ta + C2K).*c35_hr.qsr);
ywt_sonb = -c35_hr.usr.*(c35_hr.tsr + 0.51.*(best.ta + C2K).*c35_hr.qsr);
yws = (9.83./(best.ta + C2K).*abs(ywtvb)*600).^0.333; % convective velocity
ygg = sqrt(1 + 1.44*yws.^2./best.wspd.^2);  % hrly average gust factor
yhsbs = da_hr.rhoa_psl*cpa.*ywt_sonb;

[~,Le_hr] = Le_water(best.ts, ones(size(best.ts))*33); % Le in J/g, default salinity
Le_hr = Le_hr*1000; % in J/kg units

%% recompute covariance and ID fluxes

% da_hr_ttt = hourly average of turbulence QC filtered data
da_hr_ttt = interval_avg(da.dnum(ttt),da{ttt,:},dnum_hr_bin);
% insert date/time columns and convert to table
da_hr_ttt(:,3) = da_hr_ttt(:,1); da_hr_ttt(:,2) = da_hr_ttt(:,1) - datenum(da_hr_ttt(1,4)-1,12,31);
da_hr_ttt = da_hr_ttt(:,2:end); da_hr_ttt(:,3:8) = datevec(da_hr_ttt(:,2));
da_hr_ttt = array2table(da_hr_ttt,'VariableNames',da_vars);
da_hr_ttt.Properties.VariableUnits = da_units;
da_hr_ttt.Properties.VariableDescriptions = da_descript;

% da_hr_sss = hourly average of ID stress QC filtered data
da_hr_sss = interval_avg(da.dnum(sss),da{sss,:},dnum_hr_bin);
% insert date/time columns and convert to table
da_hr_sss(:,3) = da_hr_sss(:,1); da_hr_sss(:,2) = da_hr_sss(:,1) - datenum(da_hr_sss(1,4)-1,12,31);
da_hr_sss = da_hr_sss(:,2:end); da_hr_sss(:,3:8) = datevec(da_hr_sss(:,2));
da_hr_sss = array2table(da_hr_sss,'VariableNames',da_vars);
da_hr_sss.Properties.VariableUnits = da_units;
da_hr_sss.Properties.VariableDescriptions = da_descript;

zety = zu./c35_hr.L;  % hourl bulk z/L
% for C2 vars, take mean of methods a and b
ycu = nanmean1([da_hr_sss.Cu2a,da_hr_sss.Cu2b],2);
ycw = nanmean1([da_hr_sss.Cw2a,da_hr_sss.Cw2b],2);
yct = nanmean1([da_hr_ttt.Ct2a,da_hr_ttt.Ct2b],2);
ycql = nanmean1([da_hr_ttt.Cq2a,da_hr_ttt.Cq2b],2);

% avg cu and cw for smoother structure function
ycs = (ycu + 0.75*ycw)/2;
ycs(isnan(ycu)) = 0.75*ycw(isnan(ycu));
ycs(isnan(ycw)) = ycu(isnan(ycw));

% fix any complex numbers from psi_fu(zety)
ypsifu  = psi_fu(zety);
ypsifu(isnan(ypsifu)) = nanmedian1(ypsifu);

% now compute usib, tsib, qsib, etc.
yusib = sqrt(ycs.*zu^0.667./ypsifu);
ytauib = da_hr.rhoa_psl.*yusib.*yusib./ygg;
% ytsibs = -sqrt(yct.*zu^0.667./psi_ft(zety)).*sign(ywt_sonb);
yDT = best.ts - best.ta; % + c35_hr.dter;  % using delta T rather than wTsonic, ts is skin temp
ytsibs = -sqrt(yct.*zu^0.667./psi_ft(zety)).*sign(yDT);
ytsib = ytsibs - 0.51*(best.ta + C2K).*c35_hr.qsr/1000;
yhsibs = -da_hr.rhoa_psl*cpa.*yusib.*ytsibs;
yhsib = -da_hr.rhoa_psl*cpa.*yusib.*ytsib;
yqsib_lic = sign(c35_hr.qsr).*sqrt(ycql.*zu^0.667./psi_ft(zety))/1000;
ii = yqsib_lic<-3e-4 | yqsib_lic>3e-4; yqsib_lic(ii) = NaN;
yhlib_lic = -da_hr.rhoa_psl.*Le_hr.*yusib.*yqsib_lic;

% covariance fluxes - turbulence variables are filtered...
ywt = da_hr_ttt.wTs_csp + 0.51*(best.ta + C2K).*c35_hr.usr.*c35_hr.qsr/1000;
yhsc = da_hr.rhoa_psl*cpa.*ywt;   % hsc NOAA sonic
yhscs = da_hr.rhoa_psl*cpa.*da_hr_ttt.wTs_csp;

% Licor csp Hl, with separation correction fator, 1.045, no Webb correction
yhlc = da_hr.rhoa_psl.*Le_hr.*da_hr_ttt.wq_csp/1000;

ytauc = -da_hr.rhoa_psl.*da_hr_sss.wu_csp;    % streamwise stress
ytaucx = -da_hr.rhoa_psl.*da_hr_sss.wv_csp;   % cross stream stress

%% Hourly net heat flux
% net longwave relative to sea
yrnl = best.rl - 0.97*(5.67e-8*(best.ts + C2K).^4); % IR heat flux into ocean
% hnet from COARE 3.5 using ship bulk met
yhnet_C35 = 0.955*best.rs + yrnl - c35_hr.hsb - c35_hr.hlb - best.rain;
% net heat from covariance measurements
yhnet_cov = 0.955*best.rs + yrnl - yhsc - (yhlc + c35_hr.hlwebb) - best.rain;
% net heat from ID measurements
yhnet_ID = 0.955*best.rs + yrnl - yhsib - (yhlib_lic + c35_hr.hlwebb) - best.rain;

%% webb correction from hourly bulk hs and hl
yrhobar = best.Pmb.*100./(Rgas*(best.ta + C2K).*(1 + 0.61*best.qa/1000));  % hr mean air density
ywbar_t = (1 + 1.61*best.qa/1000).*c35_hr.hsb./da_hr.rhoa_psl./cpa./(best.ta + C2K);  % webb correction
ywbar_w = 1.61*c35_hr.hlb./Le_hr;

%% hourly heat flux and stress plots
% timeseries for various verisons of hourly Hl
miny = -100; maxy = 200;
figure('position',[1,1,800,475],'name','hr_hl');
% yhlc is hr ave heat flux computed from hourly filtered 10-min wq_lic
subtightplot(2,1,1,[0.08,0.08],[0.05,0.05],[0.08,0.05]);
    plot([dt_hr(1),dt_hr(end)],[0,0],'k-','LineWidth',1); hold on;
    h1 = plot(dt_hr,yhlc+c35_hr.hlwebb,'bo',dt_hr,c35_hr.hlb,'r-',...
              'MarkerSize',8,'LineWidth',1);
    xlim([stdt endt]); ylim([ miny maxy]); grid on; grid minor;
    ylabel('HL , W/m2'); legend(h1,'cov','C35','location','north');
    title('Hourly filtered covariance latent heat flux');
% yhlib_lic is hr ave heat flux from filtered hourly cq2
subtightplot(2,1,2,[0.08,0.08],[0.05,0.05],[0.08,0.05]);
    plot([dt_hr(1),dt_hr(end)],[0,0],'k-','LineWidth',1); hold on;
    h2 = plot(dt_hr,yhlib_lic+c35_hr.hlwebb,'bo',dt_hr,c35_hr.hlb,'r-',...
              'MarkerSize',8,'LineWidth',1);
    xlim([stdt endt]); ylim([ miny maxy]); grid on; grid minor;
    ylabel('HL, W/m2');  legend(h2,'ID','C35','location','north');
    title('Hourly filtered ID latent heat flux');
if savefigs
    ppath = fullfile(png_path,['067_',cruise,'_Hl_timeseries.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['067_',cruise,'_Hl_timeseries.fig']));
end

% Hl cov vs bulk
minx = -40; maxx = 60;
figure('position',[1,1,800,475],'name','hr_hlc_v_hlb');
plot([0,0],[minx,maxx],'k-',[minx,maxx],[0,0],'k-','LineWidth',1); hold on;
plot(c35_hr.hlb,yhlc+c35_hr.hlwebb,'bo',[minx maxx],[minx maxx],...
     'k--','MarkerSize',8,'LineWidth',1);
grid on; grid minor; axis([minx maxx minx maxx]);
xlabel('COARE HL (W/m^2)'); ylabel('HL Cov + Webb (W/m^2)');
title([ptitle,' Hourly Covariance Latent Heat (Webb Corr) vs Bulk']);
if savefigs
    ppath = fullfile(png_path,['068_',cruise,'_Hl_cov_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['068_',cruise,'_Hl_cov_v_bulk.fig']));
end

% Hl ID vs bulk
minx = -40; maxx = 60;
figure('position',[1,1,800,475],'name','hr_hlib_v_hlb');
plot([0,0],[minx,maxx],'k-',[minx,maxx],[0,0],'k-','LineWidth',1); hold on;
plot(c35_hr.hlb,yhlib_lic+c35_hr.hlwebb,'bo',[minx maxx],[minx maxx],...
     'k--','MarkerSize',8,'LineWidth',1);
grid on; grid minor; axis([minx maxx minx maxx]);
xlabel('COARE HL (W/m^2)'); ylabel('HL ID + Webb (W/m^2)');
title([ptitle,' Hourly ID Latent Heat (Webb Corr) vs Bulk']);
if savefigs
    ppath = fullfile(png_path,['069_',cruise,'_Hl_ID_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['069_',cruise,'_Hl_ID_v_bulk.fig']));
end

% timeseries of Hs cov
minx = min(yhsc)-20; maxx = max(yhsc)+20;
figure('position',[1,1,800,475],'name','hr_hsc');
plot(dt_hr,yhsc,'bo',dt_hr,c35_hr.hsb,'r-',[stdt endt],[0 0],'k-',...
     'MarkerSize',8,'LineWidth',1.5);
grid on; grid minor; xlim([stdt endt]); ylim([ minx maxx]);
ylabel('Cov Sensible Heat Flux (W/m^2)');
legend('NOAA','C35','location','best');
title([ptitle,' Hourly Filtered Covariance Hs']);
if savefigs
    ppath = fullfile(png_path,['070_',cruise,'_Hs_cov_timeseries.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['070_',cruise,'_Hs_cov_timeseries.fig']));
end

% timeseries of Hs ID
minx = -50; maxx = 50;
figure('position',[1,1,800,475],'name','hr_hsib');
plot(dt_hr,yhsib,'bo',dt_hr,c35_hr.hsb,'r-',[stdt endt],[0 0],'k-',...
     'MarkerSize',8,'LineWidth',1.5);
grid on; grid minor; xlim([stdt endt]); ylim([ minx maxx]);
ylabel('ID Sensible Heat Flux (W/m^2)');
legend('NOAA','C35','location','best');
title([ptitle,' Hourly Filtered ID Hs']);
if savefigs
    ppath = fullfile(png_path,['071_',cruise,'_Hs_ID_timeseries.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['071_',cruise,'_Hs_ID_timeseries.fig']));
end

% Hs cov vs bulk
% minx = min(yhsc)-20; maxx = max(yhsc)+20;
figure('position',[1,1,800,475],'name','hr_hsc_v_hsb');
plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k-'); hold on;
plot(c35_hr.hsb,yhsc,'bo','MarkerSize',8,'LineWidth',1.5); grid on; grid minor;
axis([-50 50 -50 50]); xlabel('COARE Hs (W/m^2)'); ylabel('Cov Hs (W/m^2)');
title([ptitle,' Hourly Cov Sensible Heat vs Bulk']);
if savefigs
    ppath = fullfile(png_path,['072_',cruise,'_Hs_cov_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['072_',cruise,'_Hs_cov_v_bulk.fig']));
end

% Hs ID vs bulk
minx = -50; maxx = 50;
figure('position',[1,1,800,475],'name','hr_hsib_v_hsb');
plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k-'); hold on;
plot(c35_hr.hsb,yhsib,'bo','MarkerSize',8,'LineWidth',1.5); grid on; grid minor;
axis([minx maxx minx maxx]); xlabel('COARE Hs (W/m^2)'); ylabel('ID Hs (W/m^2)');
title([ptitle,' Hourly ID Sensible Heat vs Bulk']);
if savefigs
    ppath = fullfile(png_path,['073_',cruise,'_Hs_ID_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['073_',cruise,'_Hs_ID_v_bulk.fig']));
end

% Stress timeseries, cov
minx = -0.05; maxx = 0.35;
figure('position',[1,1,800,475],'name','hr_tauc');
plot(dt_hr,ytauc,'bo',dt_hr,c35_hr.tau,'r-',[stdt endt],[0 0],'k-',...
     'MarkerSize',8,'LineWidth',1.5);
xlim([stdt endt]); ylim([ minx maxx]); ylabel('Covariance Stress (N/m^2)');
legend('NOAA','C35','location','best'); grid on; grid minor;
title([ptitle ' Hourly Filtered Covariance Stress']);
if savefigs
    ppath = fullfile(png_path,['074_',cruise,'_Tau_cov_timeseries.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['074_',cruise,'_Tau_cov_timeseries.fig']));
end

% Stress timeseries, ID
minx = -0.05; maxx = 0.35;
figure('position',[1,1,800,475],'name','hr_tauib');
plot(dt_hr,ytauib,'bo',dt_hr,c35_hr.tau,'r-',[stdt endt],[0 0],'k-',...
     'MarkerSize',8,'LineWidth',1.5);
xlim([stdt endt]); ylim([ minx maxx]); ylabel('ID Stress (N/m^2)');
legend('NOAA','C35','location','best'); grid on; grid minor;
title([ptitle ' Hourly Filtered ID Stress']);
if savefigs
    ppath = fullfile(png_path,['075_',cruise,'_Tau_ID_timeseries.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['075_',cruise,'_Tau_ID_timeseries.fig']));
end

% Stress, cov vs bulk
minx = -0.05; maxx = 0.35;
figure('position',[1,1,800,475],'name','hr_tauc_v_taub');
plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k-'); hold on;
plot(c35_hr.tau,ytauc,'bo','MarkerSize',8,'LineWidth',1.5); grid on; grid minor;
axis([minx maxx minx maxx]);
xlabel('COARE Stress (N/m^2)'); ylabel('Covariance Stress (N/m^2)');
title([ptitle,' Hourly Covariance Stress vs Bulk']);
if savefigs
    ppath = fullfile(png_path,['076_',cruise,'_Tau_cov_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['076_',cruise,'_Tau_cov_v_bulk.fig']));
end

% Stress, ID vs bulk
minx = -0.05; maxx = 0.35;
figure('position',[1,1,800,475],'name','hr_tauib_v_taub');
plot([0 0],[minx maxx],'k-',[minx maxx],[0 0],'k-',[minx maxx],[minx maxx],'k-'); hold on;
plot(c35_hr.tau,ytauib,'bo','MarkerSize',8,'LineWidth',1.5); grid on; grid minor;
axis([minx maxx minx maxx]);
xlabel('COARE Stress (N/m^2)'); ylabel('ID Stress (N/m^2)');
title([ptitle,' Hourly ID Stress vs Bulk']);
if savefigs
    ppath = fullfile(png_path,['077_',cruise,'_Tau_ID_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['077_',cruise,'_Tau_ID_v_bulk.fig']));
end

% close all

%% median hourly hl & stress vs bulk & vs U
% Wind bin latent heat
uu=0:15;
hlibm = NaN(length(uu),1); hlbm = hlibm; nn = hlibm; dqm = hlibm;
hlcm = hlibm; hlcmx = hlibm; hlibmmx = hlibm; hlbmx = hlibm; dum = hlibm;
for i=1:length(uu)
    ii = find(floor(best.wspd)==uu(i));
    jj = find(isfinite(yhlib_lic(ii)) & isfinite(c35_hr.hlb(ii)) );
    if isempty(jj)
        hlibm(i) = NaN;
        hlbm(i) = NaN;
        hlibmmx(i) = NaN;
        hlbmx(i) = NaN;
        hlcm(i) = NaN;
        hlcmx(i) = NaN;
        nn(i) = 0;
        dqm(i) = NaN;
        dum(i) = NaN;
    else
        hlibm(i) = nanmedian1(yhlib_lic(ii(jj))+c35_hr.hlwebb(ii(jj)));
        hlcm(i) = nanmedian1(yhlc(ii(jj))+c35_hr.hlwebb(ii(jj)));
        hlbm(i) = nanmedian1(c35_hr.hlb(ii(jj)));
        nn(i) = length(jj);
        dqm(i) = nanmedian1(best.qs(ii(jj))-best.qa(ii(jj)));
        hlibmmx(i) = nanmean1(yhlib_lic(ii(jj))+c35_hr.hlwebb(ii(jj)));
        hlcmx(i) = nanmean1(yhlc(ii(jj))+c35_hr.hlwebb(ii(jj)));
        hlbmx(i) = nanmean1(c35_hr.hlb(ii(jj)));
        dum(i) = nanmean1(best.wspd(ii(jj)));
    end
end

maxx = nanmax1(hlbm)+40; minx = nanmin1(hlbm)-20;
figure('position',[1,1,800,475],'name','med_hl_v_u');
plot(dum,hlbm,'rx-',dum,hlibm,'ro-',dum,hlcm,'r*-',dum,hlbmx,'bx-',...
     dum,hlibmmx,'bo-',dum,hlcmx,'b*-',[0,max(dum)+1],[0,0],...
     'k-','MarkerSize',8,'LineWidth',1);
grid on; grid minor; axis([0 max(dum)+1 minx maxx]);
xlabel('Wspd (m/s)'); ylabel('HL (W/m^2)');
legend('median HL C35','median HL ID','median HL cov','mean HL C35',...
    'mean HL ID','mean HL cov','location','north');
title([ptitle,' Wind-bin HL + Webb vs Wspd']);
if savefigs
    ppath = fullfile(png_path,['078_',cruise,'_med_Hl_v_U.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['078_',cruise,'_med_Hl_v_U.fig']));
end

figure('position',[1,1,800,475],'name','med_hl_v_hlb');
plot([0,0],[-100,300],'k-',[-100,300],[0,0],'k-',[-100,300],[-100,300],...
     'k--','LineWidth',1); hold on;
h1 = plot(hlbm,hlibm,'ro',hlbm,hlcm,'r*',hlbmx,hlibmmx,'bo',...
     hlbmx,hlcmx,'b*','MarkerSize',8,'LineWidth',1);
grid on; grid minor; axis([minx maxx minx maxx]);
xlabel('COARE HL (W/m^2)'); ylabel('Turb HL (W/m^2)');
legend(h1,'median HL ID','median HL cov','mean HL ID','mean HL cov','location','northwest');
title([ptitle,' Wind-bin HL + Webb vs bulk']);
if savefigs
    ppath = fullfile(png_path,['079_',cruise,'_med_Hl_v_bulk.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['079_',cruise,'_med_Hl_v_bulk.fig']));
end

%% Wind bin stress for NOAA sonic
taubm = NaN(length(uu),1); taubmm = taubm; taucm = taubm; taucxm = taubm;
tauibm = taubm; taucmm = taubm; taucxmm = taubm; tauibmm = taubm;
wum = taubm; wvm = taubm; nm = taubm; dum = taubm;

for i=1:length(uu)
    ii = find(floor(best.wspd)==uu(i));
    www = find(isfinite(ytauc(ii)) & isfinite(ytaucx(ii)) & isfinite(c35_hr.tau(ii)) & ...
              (da_hr.rwdir_psl(ii)>rdir_lo & da_hr.rwdir_psl(ii)<rdir_hi) & ...
               da_hr.tilt(ii)<tilt_hi & da_hr.tilt(ii)>tilt_lo);
    if isempty(www)
      taubm(i) = NaN;
      taucm(i) = NaN;
      taucxm(i) = NaN;
      tauibm(i) = NaN;
      wum(i) = NaN;
      wvm(i) = NaN;
      nm(i) = 0;
      dum(i) = NaN;
      taubmm(i) = NaN;
      taucmm(i) = NaN;
      taucxmm(i) = NaN;
      tauibmm(i) = NaN;
    else
      taubm(i) = nanmedian1(c35_hr.tau(ii(www)));
      taucm(i) = nanmedian1(ytauc(ii(www)));
      taucxm(i) = nanmedian1(ytaucx(ii(www)));
      tauibm(i) = nanmedian1(ytauib(ii(www)));
      wum(i) = nanmedian1(da_hr.wu_csp(ii(www)));
      wvm(i) = nanmedian1(da_hr.wv_csp(ii(www)));
      nm(i) = length(www);
      dum(i) = nanmedian1(best.wspd(ii(www)));
      taubmm(i) = nanmean1(c35_hr.tau(ii(www)));
      taucmm(i) = nanmean1(ytauc(ii(www)));
      taucxmm(i) = nanmean1(ytaucx(ii(www)));
      tauibmm(i) = nanmean1(ytauib(ii(www)));
    end
end

ii = nm>2;
figure('position',[1,1,800,475],'name','med_tau_v_u');
plot(dum,taubm,'k-',dum,taubmm,'k--',dum(ii),taucm(ii),'bo',dum(ii),taucmm(ii),'ro',...
    dum(ii),tauibm(ii),'go',[0 15],[0 0],'k-','MarkerSize',8,'LineWidth',1);
grid on; grid minor; xlabel('Wspd (m/s)'); ylabel('Tau (N/m^2)'); % axis([0 16 -.025 0.5]);
legend('bulk median','bulk mean','cov median','cov mean','ID median','location','northwest')
title([ptitle,' Wind-bin Averaged Hourly Streamwise Stress vs. Wind Speed']);
if savefigs
    ppath = fullfile(png_path,['080_',cruise,'_med_Tau_v_U.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['080_',cruise,'_med_Tau_v_U.fig']));
end

figure('position',[1,1,800,475],'name','med_taux_v_u');
plot(dum(ii),-0.004-tauibm(ii)/10,'k--',dum(ii),taucxm(ii),'bo',dum(ii),taucxmm(ii),'ro',...
    [0 16],[0 0],'k-','MarkerSize',8,'LineWidth',1);
grid on; grid minor; xlabel('Wspd (m/s)'); ylabel('Taux (N/m^2)'); % axis([0 16 -.05 0.05]);
legend('-0.004 - tauib/10','Taux mean','Taux median','location','southwest')
title([ptitle,' Wind-bin Averaged Hourly Cross-stream Stress vs. Wind Speed']);
if savefigs
    ppath = fullfile(png_path,['081_',cruise,'_med_Taux_v_U.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['081_',cruise,'_med_Taux_v_U.fig']));
end

figure('position',[1,1,800,475],'name','med_wu_wv_v_u');
plot(dum(ii),-wum(ii),'-bd',dum(ii),-wvm(ii),'-rd',[0 16],[0 0],'k-',...
     'MarkerSize',8,'LineWidth',1); grid on; grid minor;
xlabel('Wspd (m/s)'); ylabel('<w''u''> & <w''v''>, m^2/s^2'); %axis([0 13 -.1 0.25]);
legend('-<w''u''> median','-<w''v''> median','location','northwest')
title([ptitle,' Wind-bin Stress Components vs. Wind Speed']);
if savefigs
    ppath = fullfile(png_path,['082_',cruise,'_wu_wv_v_U.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['082_',cruise,'_wu_wv_v_U.fig']));
end

%% net heat timeseries
figure('position',[1,1,800,475],'name','hnet');
plot(dt_hr',0.955*best.rs,'r-',dt_hr,yrnl,'m-',dt_hr,-c35_hr.hsb,'b-',dt_hr,...
     -c35_hr.hlb,'g-',dt_hr,yhnet_C35,'k-','LineWidth',1.5); grid on; grid minor;
xlim([stdt,endt]); ylabel('Heat Budget Terms, (W/m^2)');
title([ptitle,' Heat Budget Terms (into Ocean): Hnet = 0.955*rs + net rl - Hs - Hl']);
legend('Solar','net IR','Hs','Hl','Net','location','best');
if savefigs
    ppath = fullfile(png_path,['083_',cruise,'_heat_budget.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['083_',cruise,'_heat_budget.fig']));
end

%% normalized sigma-W vs. z/L from filtered hourly results
psi_w_cov = sqrt(da_hr_ttt.wvar)./sqrt(-da_hr_ttt.wu_csp);  % nondim sigma-W/u*(cov)
psi_w_covb = sqrt(da_hr_ttt.wvar)./c35_hr.usr;              % nondim sigma-W/usb
aa = zety>1e-4; bb = zety<-1e-4;     % selectors for stable & convective
xzet_a = 10.^(-3:0.25:2);            % stable z/L for K&F function
xzet_b = -1*10.^(-3:0.25:2);         % convective z/L for K&F function
phi_w_a = 1.25*(1+0.2*xzet_a);            % Kaimal & Finnigan, 1994, Eq.1.33)
phi_w_b = 1.25*(1+3*abs(xzet_b)).^(1/3);  % Kaimal & Finnigan, 1994, Eq.1.33)

figure('position',[1,1,800,800],'name','sigWc_v_zet');
subtightplot(2,1,1,[0.08,0.08],[0.05,0.05],[0.08,0.05]);
    loglog(zety(bb),real(psi_w_cov(bb)),'rd','LineWidth',1.5,'MarkerSize',8); hold on;
    loglog(xzet_b,phi_w_b,'k--','linewidth',2); grid on;
    ylim([5e-1,50]);
    ylabel('\Phi_w = \sigma_w / sqrt(-w''u'')'); xlabel('z/L');
    title([ptitle,'  Non-dimensional \sigma_w vs z/L']);
subtightplot(2,1,2,[0.08,0.08],[0.05,0.05],[0.08,0.05]);
    loglog(zety(aa),real(psi_w_cov(aa)),'rd','LineWidth',1.5,'MarkerSize',8); hold on;
    loglog(xzet_a,phi_w_a,'k--','linewidth',2); grid on;
    ylim([5e-1,50]);
    ylabel('\Phi_w = \sigma_w / sqrt(-w''u'')'); xlabel('z/L');
if savefigs
    ppath = fullfile(png_path,['084_',cruise,'_nondim_sigma_W_v_zet.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['084_',cruise,'_nondim_sigma_W_v_zet.fig']));
end

figure('position',[1,1,800,800],'name','sigWb_v_zet');
subtightplot(2,1,1,[0.08,0.08],[0.05,0.05],[0.08,0.05]);
    loglog(zety(bb),real(psi_w_covb(bb)),'rd','LineWidth',1.5,'MarkerSize',8); hold on;
    loglog(xzet_b,phi_w_b,'k--','linewidth',2); grid on;
    ylim([5e-1,50]);
    ylabel('\Phi_w = \sigma_w / u_* bulk'); xlabel('z/L');
    title([ptitle,'  Non-dimensional \sigma_w vs z/L']);
subtightplot(2,1,2,[0.08,0.08],[0.05,0.05],[0.08,0.05]);
    loglog(zety(aa),real(psi_w_covb(aa)),'rd','LineWidth',1.5,'MarkerSize',8); hold on;
    loglog(xzet_a,phi_w_a,'k--','linewidth',2); grid;
    ylim([5e-1,50]);
    ylabel('\Phi_w = \sigma_w / u_* bulk'); xlabel('z/L');
if savefigs
    ppath = fullfile(png_path,['085_',cruise,'_nondim_sigma_Wb_v_zet.png']);
    print('-dpng',ppath);
%    saveas(gcf,fullfile(fig_path,['085_',cruise,'_nondim_sigma_Wb_v_zet.fig']));
end

% if any(isfinite(da.sst_rosr))
%     ii = da.sog_psl>1;
%     [xbins, ybins, mbins, SD, N] = binave2(da.dter_psl(ii),da.sst_rosr(ii)-da.sst_psl(ii),...
%                                0.05,-0.5,0.1);
%     figure;
%     plot([-1,0.1],[-1,0.1],'k--',[0,0],[-1,0.1],'k-',[-1,0.1],[0,0],'k-',...
%         da.dter_psl(ii),(da.sst_rosr(ii)-da.sst_psl(ii)),'r.','MarkerSize',8); hold on;
%     errorbar(xbins,ybins,SD./sqrt(N),'bd-','linewidth',2); grid on; grid minor;
%     xlabel('dter, °C COARE 3.5'); ylabel('rosr - tsnk, °C'); axis([-1,0.1,-1,0.1]);
%     title('Sea Surface Cool Skin Temperature Depression Comparison');
%     if savefigs
%         ppath = fullfile(png_path,['086_',cruise,'_cool_skin_delta_T.png']);
%         print('-dpng',ppath);
%         %    saveas(gcf,fullfile(fig_path,['086_',cruise,'_cool_skin_delta_T.fig']));
%     end
% end

%% save hourly results
if saveflx
    save_flux_hr_FATIMA_2022;     % save hourly flux to file if desired
end

%% END
close all
disp(['FINISHED da_red PROCESSING ',cruise]);



