%{
 FLUX EVAL STRATUS 2019 :: 2018-03 :: BWB
 Process daily met files, make daily dianostic plots and save intermediate
 data files. Saves data at 1 min resolution.

 Calls:
  auto_read_met_day:     mean2 radiation/SST/pressure and WXT520 @ 1 min
  auto_read_gps_day:     GPS @ 1 min
  auto_read_hed_day:     Heading @ 1 min
%}

%%  initialize parameters
clear;
close('all');
fclose('all');

%% other options - set as necessary

% date range: 93 to 107 for complete cruise leg
stjd = 114;  % start day
endjd = 114; % end day
yyyy = '2019';

plotit = true; % for doing plots.
prtit = true;  % for saving plots--plotit must also be true
graphformat = '.png';  % select graphics format files
graphdevice = '-dpng'; % select graphic device
cruise = 'STRATUS_2019'; % Cruise name
ship = 'CaboDeHornos';    % Research vessel

% system specific path definitions
sysType = computer;
%if strncmp(sysType,'MACI64',7)       % set Mac paths
 %   data_drive = '/Users/byronb/Raw_Data_Files.noindex';
 %   path_prog = '/Users/byronb/Documents/MATLAB/2019_STRATUS/auto_flux';
 %   path_python = '/Users/byronb/Documents/MATLAB/2019_STRATUS/python';
%elseif strncmp(sysType,'PCWIN',7)  % set PSD DAC paths
    data_drive = 'D:\Data';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs','auto_flux');
    path_python = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs','python');
%end

% matlab script path
restoredefaultpath
userpath reset
cd(path_prog);
addpath(genpath(path_prog));
savepath;
rehash toolboxcache;

% paths relative to data_drive - first one should be preexisting
path_raw_data = fullfile(data_drive,cruise,ship,'flux','Raw');
path_proc_data = fullfile(data_drive,cruise,ship,'flux','Processed');
mkdir(path_proc_data);
path_images = fullfile(data_drive,cruise,ship,'flux','Processed_Images');
mkdir(path_images);

% Track map limits
Lonmin = 250;
Lonmax = 295;
Latmin = -28;
Latmax = -10;
PosLims = [Lonmin,Lonmax,Latmin,Latmax];

% temperature adjustments, all deg C
tsea_adj = 0;
ta_adj = 0;
td1_adj = 0;
tc1_adj = 0;
td2_adj = 0;
tc2_adj = 0;
adj = [tsea_adj,ta_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 Cabo de Hornos
% WXT height (wind,P,T,RH,rain), top of bow tower
zwxt = 15.4;
% mean2 pressure sensor height (m), top of wheelhouse
zp = 16.4;
% 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 met3 files to gprm and wxt0 files with python script
    % python language must be installed on the computer.
    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
    cd (path_prog);
    
    %% read data
    % read PSD 1-min avg bulk met data, metm is 1440x19
    metm = auto_read_met_day_STRATUS_2019(path_working_ddd,ddd,yyyy,adj);
    % read PSD 1-min avg gps - gprm is 1440x7 array
    gprm = auto_read_gps_day_STRATUS_2019(path_working_ddd,ddd,yyyy,PosLims);
    % read PSD 1-min avg heading-roll - hedm is 1440x3 array
    hedm = auto_read_hed_day_STRATUS_2019(path_working_ddd,ddd,yyyy);
    
    %{
    metm is 1440x19
    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
    12  rwdir_wxt       rel wind direction ,0-360 deg from bow
    13  rwspd_wxt       rel wind speed, m/s
    14  Tair_wxt        air temperature, C
    15  RH_wxt          RH, %
    16  Pmb_wxt         air pressure, mb or hPa
    17  rain_wxt        rain rate, mm/hr
    18  U_wxt           U wind component, WXT
    19  V_wxt           V wind component, WXT
    
    gprm is 1440x7
    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
    
    hedm is 1440x3
    1   jd_ref      10 min timestamp
    2   hed         heading
    3   roll        roll angle
    %}
    
    %% recompute true winds from relative U,V components
    %  and define variables
    [wspd,wdir,~,~] = uv_rel_to_sd_true(metm(:,18),metm(:,19),hedm(:,2),gprm(:,2),gprm(:,3));
    ta = metm(:,14);
    ts = metm(:,2);
    rh = metm(:,15);
    psp1 = metm(:,3);
    psp2 = metm(:,5);
    pir1 = metm(:,4)*-1;
    pir2 = metm(:,6)*-1;
    rl = (metm(:,4)+metm(:,6))/2;      % using mean PIR
    rs = (metm(:,3)+metm(:,5))/2;      % using mean PSP
    Pmb = metm(:,11)+0.125*(zp);       % @ sea level
    Pmb_wxt = metm(:,16)+0.125*(zwxt); % @ sea level
    rain = metm(:,17);
    rwspd = metm(:,13);
    rwdir = metm(:,12);
    tc1 = metm(:,7);
    td1 = metm(:,8);
    tc2 = metm(:,9);
    td2 = metm(:,10);
    cog = gprm(:,2);
    sog = gprm(:,3);
    Ngps = gprm(:,4);
    Egps = gprm(:,5);
    lat = gprm(:,6);
    lon = gprm(:,7);
    hed = hedm(:,2);
    qs = qsea_p([ts Pmb_wxt]);
    qa = qair_p([ta rh Pmb_wxt]);
    jd_min = metm(:,1);
    
    %% Some plots...
    [month,day] = yd2md(2019,ddd);
    Vdate = [2019,month,day]; 
    if plotit
        % PRESSURE
        figure; plot(jd_min,Pmb,'r-',jd_min,Pmb_wxt,'b-','linewidth',2); grid;
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Pressure',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        xlabel('Hour (UTC)'); ylabel('Sea Level Pressure (mb)'); legend('Mean2','WXT');
        axis([ddd ddd+1 min(Pmb(~isnan(Pmb)))-1  max(Pmb(~isnan(Pmb)))+1]);
        set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'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(jd_min,tc1,'r',jd_min,td1,'y',jd_min,tc2,'m',jd_min,td2,'k',jd_min,ta,'b',jd_min,ts,'k:','linewidth',1);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD temperatures',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        legend('Tc1','Td1','Tc2','Td2','Ta','SST','Location','Best');
        xlabel('Hour (UTC)'); ylabel('Temperature (degC)');
        xlim([ddd ddd+1]); grid;
        set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'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(jd_min, rh, 'r','linewidth',2); grid;
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Relative Humidity',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlabel('Hour (UTC)'); ylabel('Relative Humidity (%)');xlim([ddd ddd+1]);
        set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'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(jd_min, rain, 'r','linewidth',2); grid;
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Rain Rate ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlabel('Hour (UTC)'); ylabel('Rain Rate (mm/hr)'); xlim([ddd ddd+1]);
        set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'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(jd_min, sog,'linewidth',2);
        xlabel('Hour (UTC)');ylabel('SOG (m/s)'); grid;
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  GPS SOG/COG',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlim([ddd ddd+1]); ylim([0 8]); set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        subplot(2,1,2); plot(jd_min, cog,'b.','linewidth',2);
        xlabel('Hour (UTC)'); ylabel('COG (deg)'); axis([ddd ddd+1 0 360]); grid;
        set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        orient tall;
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'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_STRATUS_2019(lon,lat,cruise,PosLims,ddd);
        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(jd_min, hed,'b.','linewidth',2);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Heading',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlabel('Hour (UTC)');ylabel('Heading (deg)');
        axis([ ddd ddd+1 0 360]); grid;
        set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'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;
        subplot(2,1,1);
        plot(jd_min, rwdir,'b.',...	% data points
            [ddd ddd+1],[60 60],'k:', ...
            [ddd ddd+1],[-60 -60],'k:', ...
            [ddd ddd+1],[90 90],'k-', ...
            [ddd ddd+1],[-90 -90],'k-','linewidth',1);
        xlabel('Hour (UTC)'); ylabel('Relative Wind Direction (deg)'); axis([ddd ddd+1 -180 180]); grid;
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Relative Wind Direction',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        subplot(2,1,2); plot(jd_min, rwspd,'b','linewidth',1); grid;
        xlabel('Hour (UTC)'); ylabel('Relative Wind Speed (m/s)'); xlim([ddd ddd+1]);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Relative Wind Speed',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        orient tall;
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'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;
        subplot(2,1,1);
        plot(jd_min, wdir,'b.','linewidth',1); grid;
        xlabel('Hour (UTC)'); ylabel('Wind Direction (deg)'); axis([ddd ddd+1 -5 365]);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  True Wind Direction',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
        set(gca(gcf),'YTick',0:45:360);
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        subplot(2,1,2); plot(jd_min, wspd,'b','linewidth',1); grid;
        xlabel('Hour (UTC)'); ylabel('Wind Speed (m/s)'); xlim([ddd ddd+1]);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  True Wind Speed',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        orient tall;
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'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;
    
    %% 1-min 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 in coare35vn.m for buoy intercomparison
    y = coare35vn(wspd,zwxt,ta,zwxt,rh,zwxt,Pmb_wxt,ts,rs,rl,lat,600,rain,NaN,NaN);
    
    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);
    
    Rns = rs*0.955;
    Rnl = 0.97*(rl - 5.67e-8*(ts+273.15).^4);
    hnet = Rns + Rnl - hsb - hlb - rf;
    
    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 && ~isempty(metm(isfinite(metm(:,2:end))))
        % NET HEAT FLUX
        figure; plot(jd_min, hnet,'linewidth',1); xlabel('Hour (UTC)'); ylabel('Heat flux (W/m^2)'); xlim([ddd ddd+1])
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Net Heat flux = %04.2f',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hm),'FontWeight','Bold','Interpreter','none');
        set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks'); grid;
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_images,[cruise,'_Net_Heat_Flux_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end
        
        % HEAT FLUX COMPONENTS - SEA
        figure; plot(jd_min, -hsb, jd_min, -hlb, jd_min, Rnl, jd_min, Rns, jd_min, -rf,'g.','LineWidth',2,'MarkerSize',6);
        hold on; plot(jd_min, hnet,'k','linewidth',2); xlabel('Hour (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'); grid;
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Heat Flux Components (sea)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlim([ddd ddd+1]); set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'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(jd_min, usr,'linewidth',1); xlabel('Hour (UTC)'); ylabel('u_*(m/s)'); xlim([ddd ddd+1]);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Friction Velocity from COARE algorithm 3.5',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks'); grid;
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'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
    rlclr = (0.52+0.13/60*abs(lat)+(0.082-0.03/60.*abs(lat)).*sqrt(qa)).*(5.67e-8*(ta+273.15).^4);
    
    % Solar model
    latd = nanmean1(lat);
    lond = nanmean1(lon);
    qrat = 4.0;
    watvap = iv;
    jdx = jd_min + 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
    
    flux10 = NaN(32,144);
    nd = NaN(1,144);
    zts = ones(1,144)*zsst;
    flux10(1,:) = jd_min;             % Decimal julian day at beginning of interval
    flux10(2,:) = wspd;               % psd true wind speed at zu, m/s
    flux10(3,:) = wdir;               % psd true wind direction, deg
    flux10(4,:) = ts;                 % psd seasnake T, C
    flux10(5,:) = nd;                 % ship theromsalinograph T, C  bow - NOT USED
    flux10(6,:) = nd;                 % ship theromsalinograph salinity, psu  bow - NOT USED
    flux10(7,:) = ta;                 % psd air T at zt, C
    flux10(8,:) = qs;                 % psd air specific humidity at sea surface, g/kg
    flux10(9,:) = qa;                 % psd air specific humidity, g/kg
    flux10(10,:) = rs;                % psd solar flux, w/m^2
    flux10(11,:) = rl;                % psd IR flux, w/m^2
    flux10(12,:) = rain;              % WXT precip rate, mm/hr
    flux10(13,:) = cog;               % cog from gps, m/s
    flux10(14,:) = sog;               % sog from gps, m/s
    flux10(15,:) = hed;               % heading from gps, deg
    flux10(16,:) = rwspd;             % psd WXT rel wind speed, m/s
    flux10(17,:) = rwdir;             % psd WXT rel wind direction, deg
    flux10(18,:) = lat;               % decimal latitude, deg
    flux10(19,:) = lon;               % decimal longtude, deg
    flux10(20,:) = zts;	             % Depth of SST sensor used in heat flux calc, m
    flux10(21,:) = usr;               % u-star, m/s
    flux10(22,:) = taub;              % wind stress, coare 3.5, N/m^2
    flux10(23,:) = hsb;               % sensible heat flux, coare 3.5, w/m^2
    flux10(24,:) = hlb;               % latent heat flux, coare 3.5, w/m^2
    flux10(25,:) = rf;                % rain heat flux, w/m^2
    flux10(26,:) = Urf;               % neutral wspd at 3m from COARE
    flux10(27,:) = Trf;               % neutral ta at 3m from COARE
    flux10(28,:) = RHrf;              % rh at 3m from COARE
    flux10(29,:) = Pmb_wxt;           % sea level pressure
    flux10(30,:) = rh;                % PSD RH at height zp
    flux10(31,:) = rsclr;             % shortwave clearsky model
    flux10(32,:) = rlclr;             % longwave clearsky model

    
    fpath = fullfile(path_proc_data,[cruise,'_PSD_flux_10min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
    fout = fopen(fpath,'w');
    fprintf(fout,[repmat('%20.10e ',1,32),' \n'],flux10);
    fclose all;
    flux10 = flux10';
        
    %% more plots
    
    if plotit
        if ~isempty(rlclr(isfinite(rlclr)))
            ii = find(isfinite(rlclr));
            figure; plot(jd_min, pir1,'b',jd_min, pir2,'g',jd_min(ii),rlclr(ii),'r'); xlim([ddd ddd+1])
            title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Infrared Flux',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
            xlabel('Hour (UTC)'); ylabel('IR Flux (W/m^2)'); grid; legend('PIR1', 'PIR2', 'Clearsky','Location','Best');
            set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
            xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
            annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'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
        end
        
        if ~isempty(rsclr(isfinite(rsclr)))
            ii = find(isfinite(rsclr));
            figure; plot(jd_min, psp1,'b', jd_min, psp2,'g',jd_min(ii), rsclr(ii),'r'); xlim([ddd ddd+1]);
            title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD solar flux',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
            xlabel('Hour (UTC)');ylabel('Solar Flux (W/m^2)'); grid; legend('PSP1', 'PSP2','Clearsky','Location','Best');
            set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
            xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2));
            annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'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
        end
        close all;
    end
    
    %% average to 1-hr timebase and run bulk model
    
    % straight average of everything
    jd_hr_bin = ddd:1/24:ddd+1; jd_hr_bin = jd_hr_bin';
    flux1hr = interval_avg(flux10(:,1),flux10(:,2:end),jd_hr_bin);
    
    % 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(flux1hr(:,2),zwxt,flux1hr(:,7),zwxt,flux1hr(:,30),zwxt,flux1hr(:,29),...
        flux1hr(:,4),flux1hr(:,10),flux1hr(:,11),flux1hr(:,18),600,flux1hr(:,12),NaN,NaN);
    
    % copy new bulk results into flux1hr array
    flux1hr(:,21) = y(:,8);     % u-star, m/s
    flux1hr(:,22) = y(:,3);     % wind stress, N/m^2
    flux1hr(:,23) = y(:,1);     % sensible heat flux, w/m^2
    flux1hr(:,24) = y(:,2);     % latent heat flux, w/m^2
    flux1hr(:,25) = y(:,14);    % rain heat flux, w/m^2
    flux1hr(:,26) = y(:,22);	% neutral wspd at 3m
    flux1hr(:,27) = y(:,23);	% neutral ta at 3m
    flux1hr(:,28) = y(:,24);	% rh at 3m
    
    % recompute angles from 10-min data
    [U,V] = sd_to_uv(wspd,wdir);
    ii = rwdir<0; rwdir(ii) = rwdir(ii) + 360;
    [Ur,Vr] = sd_to_uv(rwspd,rwdir);
    Nhed = cos(hed./180.*pi); % 1x3600 arrays
    Ehed = sin(hed./180.*pi);
    xx = [Ngps,Egps,U,V,Ur,Vr,Nhed,Ehed];
    xxx = interval_avg(jd_min,xx,jd_hr_bin); xxx = xxx(:,2:end);
    % wdir
    [~,flux1hr(:,3)] = uv_to_sd(xxx(:,3),xxx(:,4));
    % rwdir
    [~,temp] = uv_to_sd(xxx(:,5),xxx(:,6));
    temp = mod(temp+360,360);
    ii = temp >180; temp(ii) = temp(ii) -360;
    flux1hr(:,17) = temp;
    % cog
    flux1hr(:,13) = atan2(xxx(:,2), xxx(:,1))*180/pi;
    flux1hr(:,13) = mod(flux1hr(:,13) + 360*3, 360);
    % hed
    flux1hr(:,15) = atan2(xxx(:,8), xxx(:,7))*180/pi;
    flux1hr(:,15) = mod(flux1hr(:,15) + 360*3, 360);
    
    % add std deviation rwdir to 1-hr array
    xx = flux10(:,17); xx(xx<0) = xx(xx<0)+360;
    xx = unwrap(xx*pi/180);
    temp = interval_std(flux10(:,1),xx,jd_hr_bin); temp = temp(:,2);
    flux1hr(:,33) = temp*180/pi;
    
    %% save 1-hr data file
    
    fpath = fullfile(path_proc_data,[cruise,'_PSD_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,33),' \n'],flux1hr');
    fclose all;
    
    filecat(fullfile(path_proc_data,'STRATUS_2019_PSD_flux_10min_2019_*.txt'),...
    fullfile(path_proc_data,'STRATUS_2019_PSD_flux_10min_all.txt'),'Overwrite');
    filecat(fullfile(path_proc_data,'STRATUS_2019_PSD_flux_1hr_2019_*.txt'),...
    fullfile(path_proc_data,'STRATUS_2019_PSD_flux_1hr_all.txt'),'Overwrite');
    
end

