% flux_eval_WHOTS_2022.m
%{
    Process daily met files, make daily dianostic plots, and save
    data files at 10 min and 1 hr resolution.

    Calls python programs parseGpsFiles.py and parseWXTFiles.py to convert
    gps and wxt data into matlab-readable text format.  Python 3.x must be
    installed on the computer and executable from the command line.

    Calls these to load daily tables for each system:
    read_met_day:     mean2 and wxt data @ 10min
    read_gps_day:     GPS @ 10 min
    read_hed_day:     Heading @ 10 min
    read_scs_day;     SCS @ 10 min

    Adjusted headers and switch to COARE 3.6
    July 2022 LB
%}

%%  set up parameters, edit as desired, and run
clear;
close('all');
fclose('all');

% date range: 203:208 for all days, XXX:XXX for WHOTS only
stjd = 205 ;  % start day
endjd = 208; % end day
yyyy = 2022; YYYY = sprintf('%04i',yyyy);

plotit = true;      % for doing plots.
prtit = true;       % for saving plots--plotit must also be true
run_python = true; % set true if raw gps/wxt files are not already converted

graphformat = '.png';  % select graphics format files
graphdevice = '-dpng'; % select graphic device
cruise = 'WHOTS_2022'; % Cruise name
ship = 'RV_Sette';     % Research vessel

% system specific path definitions
sysType = computer;
if strncmp(sysType,'MACI64',7)       % set Mac paths
    data_drive = '/Users/ludo/Documents/Work/';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs','flux');
    path_python = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs','python');
    path_raw_data = fullfile(data_drive,cruise,ship,'flux','Raw');
    path_proc_data = fullfile(data_drive,cruise,ship,'flux','Processed');
    mkdir(path_proc_data); mkdir(fullfile(path_proc_data,'daily'));
    path_images = fullfile(data_drive,cruise,ship,'flux','Processed_Images','daily');
    mkdir(path_images);
elseif strncmp(sysType,'PCWIN',7)    % set PC Windows paths
    data_drive = 'C:\Data';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs','flux');
    path_python = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs','python');
    path_raw_data = fullfile(data_drive,cruise,ship,'flux','Raw');
    path_proc_data = fullfile(data_drive,cruise,ship,'flux','Processed');
    mkdir(path_proc_data); mkdir(fullfile(path_proc_data,'daily'));
    path_images = fullfile(data_drive,cruise,ship,'flux','Processed_Images','daily');
    mkdir(path_images);
end

% matlab script path
restoredefaultpath
cd(path_prog);
addpath(genpath(path_prog));
savepath;

% Track map limits
Lonmin = 201;
Lonmax = 203;
Latmin = 21;
Latmax = 23;
PosLims = [Lonmin,Lonmax,Latmin,Latmax];

% temperature adjustments, all deg C
tsea_adj = 0.06;  %based on comparison with TSG at night
ta1_adj = 0.2;  %in lab comparison
ta2_adj = 0.2;  %in lab comparison
td1_adj = -.2;  
tc1_adj = -.2;
td2_adj = -.2;
tc2_adj = -.2;
adj = [tsea_adj,ta1_adj,ta2_adj,td1_adj,tc1_adj,td2_adj,tc2_adj];

% Clear sky model configuration
k1 = 0.05;	% aerosol optical depth, band 1
k2 = 0.05;	% aerosol optical depth, band 2
oz = 0.2;	% column ozone
iv = 3.5;	% column water vapor (cm), if not calculated from obs.

% sensor heights for RV_Sette
% z_pdeck = 22+6/12 + 8+6/12 + 8+10/12; % measured height of P deck, ft. --> deck seems to be pilot house
% % WXT height (wind,P,T,RH,rain), top of bow tower
% zwxt = (z_pdeck + 5)*0.3048; % wxt height in meters
% % Vaisala PTB sensor height (m), top of wheelhouse
% zp = (z_pdeck + 2.5)*0.3048; % ptb sensor height in meters
%07/2022 - Ludo and Alex new measurements
z_pdeck=41.375;  %pilot house height (ft) = 12.6m
zwxt = (21 + 24)*0.3048; % wxt height in meters = 13.7m
% % Vaisala PTB sensor height (m), top of wheelhouse
zp = (z_pdeck + 3)*0.3048; % ptb sensor height in meters ~ 13.1m

% ship pressure sensor height
zimp = 0; % SCS pressure reported as SL pressure
% ship T/RH and wind height(m) on main mast (70 ft)
zscs = 70*0.3048;
% depth of sst, m
zsst = 0.05;  % using sea snake

%% Analysis Loop

for ddd = stjd:endjd
    %% set path to daily directory of raw files
    path_working_ddd = fullfile(path_raw_data,[YYYY(3:4),sprintf('%03i',ddd)]);
    if (exist(path_working_ddd,'dir')==0)
        continue
    end

    %% convert raw gps and wxt files to gprm and wxp files with python scripts
    % python 3.x must be installed on the computer.
    if run_python
        cd (path_working_ddd);
        status = system(['python ',fullfile(path_python,'parseGpsFiles.py')]);
        if status
            disp(['GPS FILE PROCESSING ERROR FOR ',path_working_ddd]);
            continue    % go to next day if python encounters an error
        end
        status = system(['python ',fullfile(path_python,'parseWXTFiles.py')]);
        if status
            disp(['WXT FILE PROCESSING ERROR FOR ',path_working_ddd]);
            continue    % go to next day if python encounters an error
        end
        status = system(['python ',fullfile(path_python,'parseSCSFiles4SST.py')]);
        if status
            disp(['SCS FILE PROCESSING ERROR FOR ',path_working_ddd]);
            continue    % go to next day if python encounters an error
        end
        cd (path_prog);
    end

    %% read 10-min data tables for entire day
    %{
    met is 144x27 table
        1   jd_met          10 min timestamp, decimal DOY @ start of interval
        2   Tsea            sea snake surface SST
        3   psp1            psp1 radiation, W/m2
        4   pir1            pir1 radiation, W/m2
        5   psp2            psp2 radiation, W/m2
        6   pir2            pir2 radiation, W/m2
        7   Tc1             pir1 case temperature
        8   Td1             pir1 dome temperature
        9   Tc2             pir2 case temperature
        10  Td2             pir2 dome temperature
        11  Pmb             vaisala pressure from mean2 box (no height adjustment)
        12  rwdir_wxt1      rel wind direction, WXT1, 0-360 deg from bow
        13  rwspd_wxt1      rel wind speed, WXT1, m/s
        14  Tair_wxt1       air temperature, WXT1, C
        15  RH_wxt1         RH, WXT1, %
        16  Pmb_wxt1        air pressure, WXT1, mb or hPa
        17  rain_wxt1       rain rate, WXT1, mm/hr
        18  U_wxt1          U wind component, WXT1, flow distortion corrected
        19  V_wxt1          V wind component, WXT1, flow distortion corrected
        20  rwdir_wxt2      rel wind direction, WXT2, 0-360 deg from bow
        21  rwspd_wxt2      rel wind speed, WXT2, m/s
        22  Tair_wxt2       air temperature, WXT2, C
        23  RH_wxt2         RH, WXT2, %
        24  Pmb_wxt2        air pressure, WXT2, mb or hPa
        25  rain_wxt2       rain rate, WXT2, mm/hr
        26  U_wxt2          U wind component, WXT2, flow distortion corrected
        27  V_wxt2          V wind component, WXT2, flow distortion corrected

    gps is 144x7 table
        1   jd_ref      10 min timestamp
        2   cog         cog, deg
        3   sog         sog, m/s
        4   Ngps        filtered N component of ship speed, m/s
        5   Egps        filtered E component of ship speed, m/s
        6   lat         latitude, deg
        7   lon         longitude E, 0-360

    hed is 144x3 table
        1   jd_ref      10 min timestamp
        2   hed         heading
        3   pitch       pitch angle
        4   roll        roll angle

    scs is 144 x 11 table
        1   jd_ref      1 Hz timestamp
        2   lat         latitude deg, N positive
        3   lon         longitude deg, E positive
        4   sogm        gps speed over ground, m/s
        5   imrh        RH humidity, g/kg
        6   imta        Air temp, C
        7   lrg         heading, deg
        8   imum        true wind speed, m/s
        9   imdm        true wind direction, deg
        10  imPress     air pressure, mb
    %}
    % read PSD 10-min avg bulk met data, metm is 144x29 table
    met = read_met_day_WHOTS_2022(path_working_ddd,ddd,yyyy,adj);
    % read PSD 10-min avg gps - gprm is 144x7 table
    gps = read_gps_day_WHOTS_2022(path_working_ddd,ddd,yyyy,PosLims);
    % read PSD 10-min avg heading-pitch - hed is 144x5 table
    hed = read_hed_day_WHOTS_2022(path_working_ddd,ddd,yyyy);
    % read SCS 10-min avg SCS scsm is 144x14 table
    scs = read_scs_day_WHOTS_2022(path_working_ddd,ddd,yyyy,PosLims);

    %adjust offsets based on lab comparison
%     met.Tair_wxt1=met.Tair_wxt1+.2;
%     met.Tair_wxt2=met.Tair_wxt2+.2;
    met.RH_wxt1=met.RH_wxt1-2;
    met.RH_wxt2=met.RH_wxt2-2;
  
    
    
    
    %% compute a few adjusted 10-min variables
    dt = datetime((met.jd + datenum(yyyy-1,12,31)),'convertfrom','datenum'); % timestamp for plots
    
    % recompute true winds from WXT relative U,V components - these are corrected for flow dist.
    [wspd1,wdir1,~,~] = uv_rel_to_sd_true(met.U_wxt1,met.V_wxt1,hed.hed,gps.cog,gps.sog);
    [wspd2,wdir2,~,~] = uv_rel_to_sd_true(met.U_wxt2,met.V_wxt2,hed.hed,gps.cog,gps.sog);
    rl = nanmean1([met.pir1,met.pir2],2);   % use mean PIR for rl
    rs = nanmean1([met.psp1,met.psp2],2);   % use mean PSP for rs
    Pmb = met.Pmb + 0.125*(zp);             % ptb @ sea level
    Pmb_wxt1 = met.Pmb_wxt1 + 0.125*zwxt;   % wxt1 @ sea level
    Pmb_wxt2 = met.Pmb_wxt2 + 0.125*zwxt;   % wxt2 @ sea level
    Pmb_scs = scs.Pmb + 0.125*zimp;         % scs @ sea level
    rwdir1 = met.rwdir_wxt1;
    rwdir1(rwdir1>180) = rwdir1(rwdir1>180) - 360; % wxt1 rwdir in +/- 180° format
    rwdir2 = met.rwdir_wxt2;
    rwdir2(rwdir2>180) = rwdir2(rwdir2>180) - 360; % wxt2 rwdir in +/- 180° format
    qs = qsea_p([met.Tsea, met.Pmb]);
    qa1 = qair_p([met.Tair_wxt1, met.RH_wxt1, met.Pmb_wxt1]);  % qa wxt1
    qa2 = qair_p([met.Tair_wxt2, met.RH_wxt2, met.Pmb_wxt2]);  % qa wxt2

    %% Some plots...

    Vdate = datevec(dt(1)); % [yyyy,mm,dd,hh,mm,ss] for file names
    dstr = datestr(Vdate);  % for plot titles

    if plotit
        % PRESSURE
        figure; plot(dt,Pmb,dt,Pmb_wxt1,dt,Pmb_wxt2,dt,Pmb_scs,'linewidth',2); grid on; grid minor;
        title(['WHOTS ',dstr,'  Sea Level Pressure'],'FontWeight','Bold','FontSize',14)
        xlabel('Date / Time UTC'); ylabel('Sea Level Pressure (mb)');
        legend('ptb','wxt1','wxt2','scs','Location','Best');
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/PSL/Boundary Layer Obs.'},...
                   'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_images,[cruise,'_Pressure_',sprintf('%04i_%02i_%02i_%03i',...
                             Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        % TEMPERATURES
        figure; plot(dt,met.Tc1,dt,met.Td1,dt,met.Tc2,dt,met.Td2,dt,met.Tair_wxt1,dt,met.Tair_wxt2,...
                     dt,met.Tsea,'.',dt,scs.ta,dt,scs.sst2,':','linewidth',1); grid on; grid minor;
                         ylim([23 33])
        title(['WHOTS ',dstr,'  Temperatures'],'FontWeight','Bold','FontSize',14);
        xlabel('Date / Time UTC'); ylabel('Temperature (degC)');
        legend('tc1','td1','tc2','td2','wxt1','wxt2','sea-snake','scs tair','scs tsg','Location','Best');
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/PSL/Boundary Layer Obs.'},...
                   'FontSize',6,'FitBoxToText','off','LineStyle','none');
       
        if prtit
            ppath = fullfile(path_images,[cruise,'_Temperatures_',sprintf('%04i_%02i_%02i_%03i',...
                             Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        % RELATIVE HUMIDITY
        figure; plot(dt,met.RH_wxt1,dt,met.RH_wxt2,dt,scs.rh,'linewidth',2); grid on; grid minor;
        title(['WHOTS ',dstr,'  Relative Humidity'],'FontWeight','Bold','FontSize',14);
        legend('wxt1','wxt2','scs','Location','Best');
        xlabel('Date / Time UTC'); ylabel('Relative Humidity (%)');
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/PSL/Boundary Layer Obs.'},...
                   'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_images,[cruise,'_RH_',sprintf('%04i_%02i_%02i_%03i',...
                             Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        % RAIN RATE
        figure; plot(dt,met.rain_wxt1,dt,met.rain_wxt2,'linewidth',2); grid on; grid minor;
        title(['WHOTS ',dstr,'  Rain Rate'],'FontWeight','Bold','FontSize',14);
        xlabel('Date / Time UTC'); ylabel('Rain Rate (mm/hr)');
        legend('wxt1','wxt2','Location','Best');
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/PSL/Boundary Layer Obs.'},...
                   'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_images,[cruise,'_Rainrate_',sprintf('%04i_%02i_%02i_%03i',...
                             Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        % GPS SOG/COG
        figure;
        subplot(2,1,1);
            plot(dt,gps.sog,dt,scs.sog,'linewidth',2); grid on; grid minor;
            xlabel('Date / Time UTC'); ylabel('SOG (m/s)');
            legend('psl','scs','Location','Best');
            title(['WHOTS ',dstr,'  SOG / COG'],'FontWeight','Bold','FontSize',14);
        subplot(2,1,2);
            plot(dt,gps.cog,dt,scs.cog,'linewidth',2); grid on; grid minor;
            xlabel('Date / Time UTC'); ylabel('COG (deg)');
        orient tall;
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/PSL/Boundary Layer Obs.'},...
                   'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_images,[cruise,'_COG_SOG_',sprintf('%04i_%02i_%02i_%03i',...
                             Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        % Track Map
        map_WHOTS_2022(gps.lon,gps.lat,PosLims,dstr);
        if prtit
            ppath = fullfile(path_images,[cruise,'_Track_Map_',sprintf('%04i_%02i_%02i_%03i',...
                             Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        % HEADING
        figure; plot(dt,hed.hed,'.',dt,scs.hed,'.-','markersize',10); grid on; grid minor;
        title(['WHOTS ',dstr,'  Heading'],'FontWeight','Bold','FontSize',14);
        xlabel('Date / Time UTC'); ylabel('Heading (deg)');
        legend('psl','scs','Location','Best');
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/PSL/Boundary Layer Obs.'},...
                   'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_images,[cruise,'_Heading_',sprintf('%04i_%02i_%02i_%03i',...
                             Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        % RELATIVE WIND SPEED & DIRECTION
        figure('position',[500,500,600,800])
        subplot(3,1,1)
            plot([dt(1),dt(end)],[60 60],'k:',[dt(1),dt(end)],[-60 -60],'k:',[dt(1),dt(end)],[90 90],...
                 'k-',[dt(1),dt(end)],[-90 -90],'k-','linewidth',1); hold on;
            h1 = plot(dt,rwdir1,dt,rwdir2,dt,scs.rwdir,'linewidth',2); grid on; grid minor; hold off;
            ylim([-180,180]); legend(h1,'wxt1','wxt2','scs','location','best');
            title(['WHOTS ',dstr,'  Relative Wind'],'FontWeight','Bold','FontSize',14);
            xlabel('Date / Time UTC'); ylabel('Rel Wind Dir (deg)');
        subplot(3,1,2)
            plot(dt,met.std_rwdir1,dt,met.std_rwdir2,'linewidth',2); grid on; grid minor;
            xlabel('Date / Time UTC'); ylabel('Std Dev 1-min Wdir (deg)');
            legend('wxt1','wxt2','location','best');
        subplot(3,1,3)
            plot(dt,met.rwspd_wxt1,dt,met.rwspd_wxt2,dt,scs.rwspd,'linewidth',2); grid on; grid minor;
            xlabel('Date / Time UTC'); ylabel('Rel Wind Spd (m/s)');
            legend('wxt1','wxt2','scs','location','best');
        orient tall;
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/PSL/Boundary Layer Obs.'},...
                   'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_images,[cruise,'_Relative_Winds_',sprintf('%04i_%02i_%02i_%03i',...
                             Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        % TRUE WIND SPEED & DIRECTION
        figure('position',[500,500,600,600])
        subplot(2,1,1);
            plot(dt,wdir1,dt,wdir2,dt,scs.wdir','linewidth',2); grid on; grid minor;
            xlabel('Date / Time UTC'); ylabel('Wind Dir (deg)'); ylim([0,360]);
            legend('wxt1','wxt2','scs','Location','Best');
            title(['WHOTS ',dstr,'  True Wind'],'FontWeight','Bold','FontSize',14);
        subplot(2,1,2);
            plot(dt,wspd1,dt,wspd2,dt,scs.wspd','linewidth',2); grid on; grid minor;
            xlabel('Hour (UTC)'); ylabel('Wind Speed (m/s)');
            legend('wxt1','wxt2','scs','Location','Best');
        orient tall;
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/PSL/Boundary Layer Obs.'},...
                   'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_images,[cruise,'_True_Winds_',sprintf('%04i_%02i_%02i_%03i',...
                             Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end
    end

    close all;

    %% 10-min bulk flux run

    % use 'best' spd/dir & rwdir from  wxt1 (starboard) and wxt2 (port) sensors;
    % if rwd > 60, use wxt1; rwd <-60: wxt2; rwd +/- 60: use mean of wxt1 & wxt2
    % first, compute mean spd/dir and rwspd/rwdir for bow sector case
    [um,dm,~,~] = uv_rel_to_sd_true(nanmean1([met.U_wxt1,met.U_wxt2],2),...
                  nanmean1([met.V_wxt1,met.V_wxt2],2),hed.hed,gps.cog,gps.sog);
    [rwm,rdm] = uv_to_sd(nanmean1([met.U_wxt1,met.U_wxt2],2),nanmean1([met.V_wxt1,met.V_wxt2],2));
    % starboard case - use wxt1
    rwdir = rwdir1;
    rwspd = met.rwspd_wxt1;
    u = wspd1;
    dir = wdir1;
    % bow sector case - use mean computed above
    rwdir(abs(rdm)<60) = rdm(abs(rdm)<60);
    rwspd(abs(rdm)<60) = rwm(abs(rdm)<60);
    u(abs(rdm)<=60) = um(abs(rdm)<=60);
    dir(abs(rdm)<=60) = dm(abs(rdm)<=60);
    % port case - use wxt2
    rwdir(rdm<-60) = rwdir2(rdm<-60);
    rwspd(rdm<-60) = met.rwspd_wxt2(rdm<-60);
    u(rdm<-60) = wspd2(rdm<-60);
    dir(rdm<-60) = wdir2(rdm<-60);

    % for COARE use mean u,ta,rh,rain values from the WXT sensors
    % use SL pressure from ptb sensor, SST from sea snake
    % Reference height is set at 3m in coare35vn.m for mooring intercomparison
    ta = nanmean1([met.Tair_wxt1,met.Tair_wxt2],2);
    rh = nanmean1([met.RH_wxt1,met.RH_wxt2],2);
    qa = qair_p([ta,rh,Pmb]);
    rain = nanmean1([met.rain_wxt1,met.rain_wxt2],2);
    
    % Reference height of 3m in coare35vn.m for mooring intercomparison
%     y = coare35vn(u,zwxt,ta,zwxt,rh,zwxt,Pmb,met.Tsea,rs,rl,gps.lat,600,rain,NaN,NaN,3);
    y=coare36vn_zrf_et(u,zwxt,ta,zwxt,rh,zwxt,Pmb,met.Tsea,rs,rl,gps.lat,gps.lon,met.jd,600,rain,rain.*35,NaN,NaN,3,3,3);
    y = real(y); % in case any complex values are returned...
%     A=[usr tau hsb hlb hbb hsbb hlwebb tsr qsr zo  zot zoq Cd Ch Ce  L  zeta dT_skinx dq_skinx dz_skin Urf Trf Qrf RHrf UrfN TrfN QrfN  lw_net sw_net Le rhoa UN U10 U10N Cdn_10 Chn_10 Cen_10 hrain Qs Evap T10 T10N Q10 Q10N  RH10 P10 rhoa10 gust wc_frac Edis];

    % Reference height of 21m in coare35vn.m for mooring intercomparison
%     yscs = coare35vn(u,zwxt,ta,zwxt,rh,zwxt,Pmb,met.Tsea,rs,rl,gps.lat,600,rain,NaN,NaN,zscs);
    yscs = coare36vn_zrf_et(u,zwxt,ta,zwxt,rh,zwxt,Pmb,met.Tsea,rs,rl,gps.lat,gps.lon,met.jd,600,rain,rain.*35,NaN,NaN,zscs,zscs,zscs);
    yscs = real(yscs);
    
    % select desired bulk model parameters...
%     hsb = y(:,1); hlb = y(:,2); taub = y(:,3); usr = y(:,8); tcol = y(:,11);
%     qcol = y(:,12); rf = y(:,14); zo = y(:,4);
%     Urf = y(:,22); Trf = y(:,23); RHrf = y(:,24);
%     Urf_scs = yscs(:,22); Trf_scs = yscs(:,23); RHrf_scs = yscs(:,24);
    hsb = y(:,3); hlb = y(:,4); taub = y(:,2); usr = y(:,1); tcol = y(:,18);
    qcol = y(:,19)./1000; rf = y(:,38); zo = y(:,10);
    Urf = y(:,21); Trf = y(:,22); RHrf = y(:,24);
    Urf_scs = yscs(:,21); Trf_scs = yscs(:,22); RHrf_scs = yscs(:,24);

    Rns = rs*0.955;  % downwelling solar flux
    Rnl = 0.97*(rl - 5.67e-8*(met.Tsea+273.15).^4);  % net IR flux
    hnet = Rns + Rnl - hsb - hlb - rf; % hnet with respect to ocean

    hm = nanmean1(hnet); sm = nanmean1(hsb); lm = nanmean1(hlb);
    Sm = nanmean1(Rns);  Lm = nanmean1(Rnl); rfm = nanmean1(rf);

    %% 10-min bulk model plots

    if plotit && any(isfinite(usr))
        % HEAT FLUX COMPONENTS - SEA
        figure; plot(dt,-hsb,dt,-hlb,dt,Rnl,dt,Rns,dt,-rf,'LineWidth',2); grid on; grid minor; hold on;
        plot(dt, hnet,'k:','linewidth',2); hold off; xlabel('Date / Time UTC'); ylabel('Heat Flux (W/m^2)');
        legend(sprintf('H_{sens}= %4.1f',-sm),sprintf('H_{lat}= %4.1f',-lm),sprintf('R_{long}= %4.1f', Lm), sprintf('R_{short}= %4.1f',Sm),...
            sprintf('Rain= %4.1f',-rfm),sprintf('Net= %4.1f',nanmean1(hnet)),'Location','Best');
        title(['WHOTS ',dstr,'  Heat Fluxes (w/respect to Sea)'],'FontWeight','Bold','FontSize',14);
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/PSL/Boundary Layer Obs.'},...
                   'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_images,[cruise,'_Heat_Flux_Components_sea_',sprintf('%04i_%02i_%02i_%03i',...
                             Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        % BULK FRICTION VELOCITY
        figure; plot(dt,usr,'linewidth',2); grid on; grid minor;
        xlabel('Date / Time UTC'); ylabel('u_*(m/s)');
        title(['WHOTS ',dstr,'  Friction Velocity, COARE 3.6'],'FontWeight','Bold','FontSize',14);
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/PSL/Boundary Layer Obs.'},...
                   'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_images,[cruise,'_Friction_Velocity_',sprintf('%04i_%02i_%02i_%03i',...
                             Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        close all;
    end

    %% radiation models @ 10 min resolution
    % Downward IR
    rlb = 5;  % CF location specific model bias for WHOTS
    rlclr = (0.52+0.13/60*abs(gps.lat)+(0.082-0.03/60.*abs(gps.lat)).*sqrt(qa)).*(5.67e-8*(ta+273.15).^4) + rlb;

    % Solar model
    latd = nanmean1(gps.lat);
    lond = nanmean1(gps.lon);
    qrat = 4.0;
    watvap = iv;
    jdx = met.jd + double(5/1440); % 10-min jd bin centers
    tutc = (jdx - floor(jdx)).*24;  % decimal hour of the day, 0-24
    clear sw sz saz dirs sky;
    [sw,sz,saz,dirs,sky] = SolarRadiancex(latd,lond,jdx,tutc,watvap,nanmean1(Pmb),k1,k2,oz);
    rsclr = sw;

    %% save 10-min bulk results to txt file - historical format with extra columns

    flx10 = NaN(144,70);
    nd = NaN(144,1);
    zts = ones(144,1)*zsst;
    zu = ones(144,1)*zwxt;
    zship = ones(144,1)*zscs; % height of scs U,Ta,RH

    flx10(:,1) = met.jd;             % decimal julian day at beginning of interval
    flx10(:,2) = u;                  % mean true wind speed at zwxt, m/s
    flx10(:,3) = dir;                % mean true wind direction, deg
    flx10(:,4) = met.Tsea;           % psd seasnake T, °C
    flx10(:,5) = nd;                 % ship theromsalinograph T, °C  bow - NOT USED
    flx10(:,6) = nd;                 % ship theromsalinograph salinity, psu  bow - NOT USED
    flx10(:,7) = ta;                 % mean air T at zu, °C
    flx10(:,8) = qs;                 % sea surface air specific humidity, g/kg
    flx10(:,9) = qa;                 % mean air specific humidity, g/kg
    flx10(:,10) = rs;                % mean solar flux, w/m^2
    flx10(:,11) = rl;                % mean IR flux, w/m^2
    flx10(:,12) = rain;              % mean WXT precip rate, mm/hr
    flx10(:,13) = gps.cog;           % cog from gps, m/s
    flx10(:,14) = gps.sog;           % sog from gps, m/s
    flx10(:,15) = hed.hed;           % heading from gps, deg
    flx10(:,16) = rwspd;             % blended psd WXT rel wind speed, m/s
    flx10(:,17) = rwdir;             % blended psd WXT rel wind direction, deg
    flx10(:,18) = gps.lat;           % decimal latitude, deg
    flx10(:,19) = gps.lon;           % decimal longtude, deg
    flx10(:,20) = zts;               % Depth of SST sensor used in heat flux calc, m
    flx10(:,21) = usr;               % u-star, m/s
    flx10(:,22) = taub;              % wind stress, COARE 3.6, N/m^2
    flx10(:,23) = hsb;               % sensible heat flux, COARE 3.6, w/m^2
    flx10(:,24) = hlb;               % latent heat flux, COARE 3.6, w/m^2
    flx10(:,25) = rf;                % rain heat flux, w/m^2
    flx10(:,26) = Urf;               % wspd at 3m from COARE
    flx10(:,27) = Trf;               % ta at 3m from COARE
    flx10(:,28) = RHrf;              % rh at 3m from COARE
    flx10(:,29) = Pmb;               % sea level pressure, ptb sensor, mb
    flx10(:,30) = rh;                % mean rh at height zu
    flx10(:,31) = rsclr;             % shortwave clearsky model
    flx10(:,32) = rlclr;             % longwave clearsky model
    flx10(:,33) = scs.ta;            % scs Tair, °C
    flx10(:,34) = scs.rh;            % scs RH, %
    flx10(:,35) = Pmb_scs;           % scs SL pressure, mb
    flx10(:,36) = scs.qa;            % scs specific humidity, g/kg
    flx10(:,37) = scs.wdir;          % scs true wind dir, degrees
    flx10(:,38) = scs.wspd;          % scs true wind spd, m/s
    flx10(:,39) = zu;                % PSL wxt sensor height, m
    flx10(:,40) = met.std_rwdir1;    % std dev 1-min rwdir1, deg
    flx10(:,41) = Urf_scs;           % wspd at height zscs, m/s
    flx10(:,42) = Trf_scs;           % ta at height zscs, °C
    flx10(:,43) = RHrf_scs;          % rh at height zscs, %
    flx10(:,44) = zship;             % height of ship wind,T and RH sensors, m
    flx10(:,45) = met.std_rwdir2;    % std dev 1-min rwdir2, deg
    flx10(:,46) = met.rwdir_wxt1;    % next 8 vars are raw wxt1 ouptut
    flx10(:,47) = met.rwspd_wxt1;
    flx10(:,48) = met.Tair_wxt1;
    flx10(:,49) = met.RH_wxt1;
    flx10(:,50) = met.Pmb_wxt1;
    flx10(:,51) = met.rain_wxt1;
    flx10(:,52) = met.U_wxt1;        % U (bow-stern) wind component/0.95
    flx10(:,53) = met.V_wxt1;        % V (port-stbd) wind component/1.15
    flx10(:,54) = met.rwdir_wxt2;    % next 8 vars are raw wxt2 output
    flx10(:,55) = met.rwspd_wxt2;
    flx10(:,56) = met.Tair_wxt2;
    flx10(:,57) = met.RH_wxt2;
    flx10(:,58) = met.Pmb_wxt2;
    flx10(:,59) = met.rain_wxt2;
    flx10(:,60) = met.U_wxt2;        % U (bow-stern) wind component/0.95
    flx10(:,61) = met.V_wxt2;        % V (port-stbd) wind component/1.15
    flx10(:,62) = met.psp1;          % next 8 vars are raw radiometer output
    flx10(:,63) = met.pir1;
    flx10(:,64) = met.psp2;
    flx10(:,65) = met.pir2;
    flx10(:,66) = met.Tc1;
    flx10(:,67) = met.Td1;
    flx10(:,68) = met.Tc2;
    flx10(:,69) = met.Td2;
    flx10(:,70) = hed.hed_std;       % std dev heading, deg
    flx10(:,71) = scs.sst2;       % tsg sst, C
    flx10(:,72) = scs.sal;       % salinity, %
    
    % we will save the hourly flux file as headerless .txt for now
    % and concat all the hourly files when daily processing is complete
    fpath = fullfile(path_proc_data,'daily',[cruise,'_PSL_flux_10min_',sprintf('%04i_%02i_%02i_%03i',...
                     Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
    % write out to file
    fout = fopen(fpath,'w');
    fprintf(fout,[repmat('%20.10e ',1,72),' \n'],flx10');
    fclose all;

    %% IR / Solar flux plots

    if plotit
        ii = find(isfinite(rlclr));
        figure; plot(dt,met.pir1,dt,met.pir2,dt(ii),rlclr(ii),'linewidth',2); grid on; grid minor;
        title(['WHOTS ',dstr,'  IR Flux'],'FontWeight','Bold','FontSize',14);
        xlabel('Date / Time UTC'); ylabel('IR Flux (W/m^2)'); legend('pir1', 'pir2', 'clearsky','Location','Best');
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/PSL/Boundary Layer Obs.'},...
            'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_images,[cruise,'_IRflux_Comparison_',sprintf('%04i_%02i_%02i_%03i',...
                Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        ii = find(isfinite(rsclr));
        figure; plot(dt,met.psp1,dt,met.psp2,dt(ii),rsclr(ii),'linewidth',2); grid on; grid minor;
        title(['WHOTS ',dstr,'  Solar Flux'],'FontWeight','Bold','FontSize',14);
        xlabel('Hour (UTC)');ylabel('Solar Flux (W/m^2)'); grid; legend('PSP1', 'PSP2','Clearsky','Location','Best');
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/PSL/Boundary Layer Obs.'},...
            'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_images,[cruise,'_Solarflux_Comparison_',sprintf('%04i_%02i_%02i_%03i',...
                Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        close all;
    end

    %% average to 1-hr timebase and rerun bulk model

    % straight average of cols 2-44 from 10-min output
    jd_hr_bin = ddd:1/24:ddd+1; jd_hr_bin = jd_hr_bin';
    flx1hr = interval_avg(flx10(:,1),[flx10(:,2:44) flx10(:,71:72)],jd_hr_bin);
%     flx1hr(:,45) = interval_avg(flx10(:,1),flx10(:,71),jd_hr_bin); %sst2
%     flx1hr(:,46) = interval_avg(flx10(:,1),flx10(:,72),jd_hr_bin); %salinity
    
    % 1-hr bulk flux run
    % Vectorized COARE 3.5 with output adjusted to match historical cor30b
    % y = [hsb hlb tau zo zot zoq L usr tsr qsr dter dqer tkt RF hlwebb Cd Ch
    %       1   2   3   4  5   6  7  8   9   10  11   12   13 14   15   16 17
    %      Ce Cdn_10 Chn_10 Cen_10 Urf  Trf  RHrf];
    %      18  19     20     21     22   23   24
    % pressure at sealevel
    % reference height set at 3m for buoy intercomparison
%     y = coare35vn(flx1hr(:,2),zwxt,flx1hr(:,7),zwxt,flx1hr(:,30),zwxt,flx1hr(:,29),...
%          flx1hr(:,4),flx1hr(:,10),flx1hr(:,11),flx1hr(:,18),600,flx1hr(:,12),NaN,NaN,3);
    y=coare36vn_zrf_et(flx1hr(:,2),zwxt,flx1hr(:,7),zwxt,flx1hr(:,30),zwxt,flx1hr(:,29),flx1hr(:,4),flx1hr(:,10),flx1hr(:,11),flx1hr(:,18),flx1hr(:,19),flx1hr(:,1),600,flx1hr(:,12),flx1hr(:,12).*35,NaN,NaN,3,3,3);
    y = real(y);
    
    % Reference height of zscs in coare35vn.m for ship intercomparison
%     yscs = coare35vn(flx1hr(:,2),zwxt,flx1hr(:,7),zwxt,flx1hr(:,30),zwxt,flx1hr(:,29),...
%          flx1hr(:,4),flx1hr(:,10),flx1hr(:,11),flx1hr(:,18),600,flx1hr(:,12),NaN,NaN,zscs);
    yscs=coare36vn_zrf_et(flx1hr(:,2),zwxt,flx1hr(:,7),zwxt,flx1hr(:,30),zwxt,flx1hr(:,29),flx1hr(:,4),flx1hr(:,10),flx1hr(:,11),flx1hr(:,18),flx1hr(:,19),flx1hr(:,1),600,flx1hr(:,12),flx1hr(:,12).*35,NaN,NaN,zscs,zscs,zscs);
    yscs = real(yscs);

    % copy new bulk results into flux1hr array
    flx1hr(:,21) = y(:,1);     % u-star, m/s
    flx1hr(:,22) = y(:,2);     % wind stress, N/m^2
    flx1hr(:,23) = y(:,3);     % sensible heat flux, w/m^2
    flx1hr(:,24) = y(:,4);     % latent heat flux, w/m^2
    flx1hr(:,25) = y(:,38);    % rain heat flux, w/m^2
    flx1hr(:,26) = y(:,21);	   % wspd at 3m
    flx1hr(:,27) = y(:,22);	   % ta at 3m
    flx1hr(:,28) = y(:,24);	   % rh at 3m
    flx1hr(:,41) = yscs(:,21); % wspd at zscs
    flx1hr(:,42) = yscs(:,22); % ta at zscs
    flx1hr(:,43) = yscs(:,24); % rh at zscs

    % recompute angles from the 10-min data
    [U,V] = sd_to_uv(u,dir);
    ii = rwdir<0; rwdir(ii) = rwdir(ii) + 360; % back to 0-360 format
    [Ur,Vr] = sd_to_uv(rwspd,rwdir);
    Nhed = cos(hed.hed./180.*pi);
    Ehed = sin(hed.hed./180.*pi);
    xx = [gps.Ngps,gps.Egps,U,V,Ur,Vr,Nhed,Ehed];
    xxx = interval_avg(flx10(:,1),xx,jd_hr_bin); xxx = xxx(:,2:end);
    % wdir
    [~,flx1hr(:,3)] = uv_to_sd(xxx(:,3),xxx(:,4));
    % rwdir
    [~,temp] = uv_to_sd(xxx(:,5),xxx(:,6));
    temp = mod(temp+10*360,360);
    ii = temp >180; temp(ii) = temp(ii) -360;
    flx1hr(:,17) = temp;
    % cog
    flx1hr(:,13) = atan2(xxx(:,2), xxx(:,1))*180/pi;
    flx1hr(:,13) = mod(flx1hr(:,13) + 360*3, 360);
    % hed
    flx1hr(:,15) = atan2(xxx(:,8), xxx(:,7))*180/pi;
    flx1hr(:,15) = mod(flx1hr(:,15) + 360*3, 360);

    % update col 40 to hourly std dev of 10-min rwdir
    xx = flx10(:,17); xx(xx<0) = xx(xx<0)+360;
    xx = unwrap(xx*pi/180);
    temp = interval_std(flx10(:,1),xx,jd_hr_bin); temp = temp(:,2);
    flx1hr(:,40) = temp*180/pi;

    %% save 1-hr data file as headerless txt

    fpath = fullfile(path_proc_data,'daily',[cruise,'_PSL_flux_1hr_',sprintf('%04i_%02i_%02i_%03i',...
                     Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
    fout = fopen(fpath,'w');
    fprintf(fout,[repmat('%20.10e ',1,46),' \n'],flx1hr');
    fclose all;

end

%% concat all the daily 10-min and hourly result files and reload
filecat(fullfile(path_proc_data,'daily','WHOTS_2022_PSL_flux_10min_2022_*.txt'),...
        fullfile(path_proc_data,'WHOTS_2022_PSL_flux_10min_all.txt'),'Overwrite');
filecat(fullfile(path_proc_data,'daily','WHOTS_2022_PSL_flux_1hr_2022_*.txt'),...
        fullfile(path_proc_data,'WHOTS_2022_PSL_flux_1hr_all.txt'),'Overwrite');

clear flx10 flx1hr
flx10 = load(fullfile(path_proc_data,'WHOTS_2022_PSL_flux_10min_all.txt'),'-ascii');
flx1hr = load(fullfile(path_proc_data,'WHOTS_2022_PSL_flux_1hr_all.txt'),'-ascii');

%% save concatenated daily results in matlab table format with var name & unit attributes

vars = {'jd','wspd','wdir','tsea','nd1','nd2','ta','qs','qa','rs','rl','rain',...
        'cog','sog','hed','rwspd','rwdir','lat','lon','zts_psl','usr','taub',...
        'hsb','hlb','rf','Urf','Trf','RHrf','Pmb','rh','rsclr','rlclr','ta_scs',...
        'rh_scs','Pmb_scs','qa_scs','wdir_scs','wspd_scs','zu_psl','std_rwdir1',...
        'Urf_scs','Trf_scs','RHrf_scs','zu_scs','std_rwdir2','rwdir_wxt1','rwspd_wxt1',...
        'Tair_wxt1','RH_wxt1','Pmb_wxt1','rain_wxt1','U_wxt1','V_wxt1','rwdir_wxt2',...
        'rwspd_wxt2','Tair_wxt2','RH_wxt2','Pmb_wxt2','rain_wxt2','U_wxt2','V_wxt2',...
        'psp1','pir1','psp2','pir2','Tc1','Td1','Tc2','Td2','std_hed','tsg2','sal'};
units = {'dec day','m/s','deg','degC','none','none','degC','g/kg','g/kg','W/m2',...
        'W/m2','mm/hr','deg','m/s','deg','m/s','deg','deg','deg','m','m/s','N/m2',...
        'W/m2','W/m2','W/m2','m/s','degC','pct','mbar','pct','W/m2','W/m2','degC',...
        'pct','mb','g/kg','deg','m/s','m','deg','m/s','degC','pct','m','deg','deg',...
        'm/s','degC','pct','mb','mm/hr','m/s','m/s','deg','m/s','degC','pct','mb',...
        'mm/hr','m/s','m/s','W/m2','W/m2','W/m2','W/m2','degC','degC','degC','degC','deg','degC','none'};

spath10 = fullfile(path_proc_data,'WHOTS_2022_PSL_flux_10min_all.mat');
spath1hr = fullfile(path_proc_data,'WHOTS_2022_PSL_flux_1hr_all.mat');

flx = array2table(flx10,'VariableNames',vars);
flx.Properties.VariableUnits = units;
save(spath10,'flx');

flxhr = array2table(flx1hr,'VariableNames',[vars(1:44) vars(71:end)]);
flxhr.Properties.VariableUnits = [units(1:44) units(71:end)];
save(spath1hr,'flxhr');

%% re-save concatenated txt data with 3 header lines

% format header strings
hdr1 = 'WHOTS 2022, NOAA/PSL flux/met 10-min result file, Rev 1, July, 2022, LWB\n';
hdr2 = ''; hdr3 = '';
for ii = 1:length(vars)
    n = vars{ii}; ln = length(n); padn = repmat(' ',1,21-ln);
    hdr2 = [hdr2,n,padn];
    u = units{ii}; lu = length(u); padu = repmat(' ',1,21-lu);
    hdr3 = [hdr3,u,padu];
end
hdr2 = [hdr2,'\n']; hdr3 = [hdr3,'\n'];

% save 10min results to ascii format
fout = fopen(fullfile(path_proc_data,'WHOTS_2022_PSL_flux_10min_all.txt'),'w');
fprintf(fout,hdr1);
fprintf(fout,hdr2);
fprintf(fout,hdr3);
fprintf(fout,[repmat('%20.10e ',1,72),' \n'],flx10');
fclose all;

% save 1hr results to same ascii format - only first 44 variables
hdr1 = 'WHOTS 2022, NOAA/PSL flux/met 1-hr result file, Rev 1, July, 2022, LB\n';
hdr2 = ''; hdr3 = '';
for ii = [1:44 71:72]
    n = vars{ii}; ln = length(n); padn = repmat(' ',1,21-ln);
    hdr2 = [hdr2,n,padn];
    u = units{ii}; lu = length(u); padu = repmat(' ',1,21-lu);
    hdr3 = [hdr3,u,padu];
end
hdr2 = [hdr2,'\n']; hdr3 = [hdr3,'\n'];
fout = fopen(fullfile(path_proc_data,'WHOTS_2022_PSL_flux_1hr_all.txt'),'w');
fprintf(fout,hdr1);
fprintf(fout,hdr2);
fprintf(fout,hdr3);
fprintf(fout,[repmat('%20.10e ',1,46),' \n'],flx1hr');
fclose all;

