%{
 AUTOFLUX STRATUS 2017 :: 2017-03 :: BWB
 Processes daily met files, make daily dianostic plots and save intermediate
 data files. Averages to 10 min. Transmits summary plots via ftp to ESRL.

 Calls:
  auto_read_met_day:     PSD radiation/SST data from logger 2 and WXT520 @ 1 min
  auto_read_gps_day:     PSD GPS data
  auto_read_hed_day:     PSD Heading data
%}

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

%% initialize run

% other options - set as necessary
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_2017';    % Cruise name
ship = 'RHB';               % 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/2017_STRATUS/';
elseif strncmp(sysType,'PCWIN',7)  % set PSD DAC paths
    data_drive = 'D:\Data';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs');
end

% set script paths
restoredefaultpath              % reset matlab path
addpath(genpath(fullfile(path_prog,'auto_flux')));	% add scripts to path
savepath; rehash;
path_python = fullfile(path_prog,'python');

% 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','auto_processed');
mkdir(path_proc_data);
path_ftp_data = fullfile(data_drive,cruise,ship,'flux','ftp');
mkdir(path_ftp_data);  % directory for ftp files
path_raw_images = fullfile(data_drive,cruise,ship,'flux','auto_images');
mkdir(path_raw_images);

% this script auto-runs after 0000 hrs UTC,
% so we define ddd to run the previous day's data
dv = clock; % current date-time vector
ddd = md2yd_vn(dv(1),dv(2),dv(3)) - 1;
[month,day] = yd2md(dv(1),ddd);
Vdate = [dv(1),month,day];      % ddd date vector
yyyy = sprintf('%04i',dv(1));   % year string

% start command diary
logfname = fullfile(path_proc_data,['MlabLog_',sprintf('%03u',ddd),'.log']);
diary(logfname);
disp('Starting AUTOflux Analysis');

% reasonable lat/lon ship track limits
Lonmin = 269;
Lonmax = 281;
Latmin = -19;
Latmax = 33;
PosLims = [Lonmin,Lonmax,Latmin,Latmax];
rdcon = pi/180;

% 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 height
% Riegl is 16.5 above water - waterline to base of tower = 7.16 m
% WXT height (m)
zwxt = 18;

%% Analysis Loop
try
    %% set raw data folder and parse wxt files
    path_working_ddd = fullfile(path_raw_data,[yyyy(3:4),sprintf('%03i',ddd)]);
    if (exist(path_working_ddd,'dir')==0)
        disp('Daily Raw Data Folder Not Found');
        return    % quit if daily folder is not found
    end
    % convert PSD gps and met3 files to gprm and wxt format 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('parseGpsFiles Error');
        return    % quit if python encounters an error
    end
    status = system(['python ',fullfile(path_python,'parseWXTFiles.py')]);
    if status
        disp('parseWXTFiles Error');
        return    % quit if python encounters an error
    end
    
    cd(fullfile(path_prog,'auto_flux'));
    
    %% read data
    % read PSD 10-min avg bulk met data, metm is 144x18
    metm = auto_read_met_day_STRATUS_2017(path_working_ddd,ddd,yyyy,adj);
    % read PSD 10-min avg gps - gprm is 144x7 array
    gprm = auto_read_gps_day_STRATUS_2017(path_working_ddd,ddd,yyyy,PosLims);
    % read PSD 10-min avg heading-roll - hedm is 144x3 array
    hedm = auto_read_hed_day_STRATUS_2017(path_working_ddd,ddd,yyyy);
    
    %{
    metm is 144x18
    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  rwdir_wxt       rel wind direction ,0-360 deg from bow
    12  rwspd_wxt       rel wind speed, m/s
    13  Tair_wxt        air temperature, C
    14  RH_wxt          RH, %
    15  Pmb_wxt         air pressure, mb or hPa
    16  rain_wxt        rain rate, mm/hr
    17  U_wxt           rel U wind component, WXT
    18  V_wxt           rel V wind component, WXT
    
    gprm is 144x7
    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 144x3
    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(:,17),metm(:,18),hedm(:,2),gprm(:,2),gprm(:,3));
    ta = metm(:,13);
    ts = metm(:,2);
    rh = metm(:,14);
    psp1 = metm(:,3);
    psp2 = metm(:,5);
    pir1 = metm(:,4);
    pir2 = metm(:,6);
    rl = (metm(:,4)+metm(:,6))/2;
    rs = (metm(:,3)+metm(:,5))/2;
    Pmb = metm(:,15)+0.125*(zwxt); % @ sea level
    rain = metm(:,16);
    rwspd = metm(:,12);
    rwdir = metm(:,11);
    tc1 = metm(:,7);
    td1 = metm(:,8);
    tc2 = metm(:,9);
    td2 = metm(:,10);
    cog = gprm(:,2);
    sog = gprm(:,3);
    lat = gprm(:,6);
    lon = gprm(:,7);
    hed = hedm(:,2);
    qs = qsea_p([ts Pmb]);
    qa = qair_p([ta rh Pmb]);
    jd_10min = metm(:,1);
    
    %% Some plots...
    
    if plotit
        % PRESSURE
        figure; plot(jd_10min,Pmb,'r-','linewidth',1); 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)');
        axis([ddd ddd+1 median(Pmb(~isnan(Pmb)))-10  median(Pmb(~isnan(Pmb)))+10]);
        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_raw_images,[cruise,'_Pressure_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end
        
        % TEMPERATURES
        figure; plot(jd_10min,tc1,'r',jd_10min,td1,'y',jd_10min,tc2,'m',jd_10min,td2,'k',jd_10min,ta,'b',jd_10min,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_raw_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_10min, rh, 'r','linewidth',1); 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_raw_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_10min, rain, 'r','linewidth',1); 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_raw_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_10min, sog,'linewidth',1);
        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_10min, cog,'b.');
        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_raw_images,[cruise,'_COG_SOG_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end
        
        % HEADING
        figure; plot(jd_10min, hed,'b.');
        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_raw_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_10min, 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_10min, 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_raw_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_10min, 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_10min, 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_raw_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
    % 10-min averaging was done above
    % using 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 Ce Cdn_10 Chn_10 Cen_10 Urf  Trf  RHrf];
    %       1   2   3   4  5   6  7  8   9   10  11   12   13 14   15   16 17 18  19     20     21     22   23   24
    
    % pressure at sealevel
    % reference height set at 3m for buoy intercomparison
    y = coare35vn(wspd,zwxt,ta,zwxt,rh,zwxt,Pmb,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_10min, 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_raw_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_10min, -hsb, jd_10min, -hlb, jd_10min, Rnl, jd_10min, Rns, jd_10min, -rf,'g.','LineWidth',2,'MarkerSize',6);
        hold on; plot(jd_10min, 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','BestOutside'); 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_raw_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_10min, 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_raw_images,[cruise,'_Friction_Velocity_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end
        
        close all;
    end
    
    %% save 10-min bulk results to txt file - historical format
    
    jazz = NaN(28,144);
    nd = NaN(1,144);
    jazz(1,:) = jd_10min;           % Decimal julian day at beginning of line
    jazz(2,:) = wspd;               % psd true wind speed at zu, m/s
    jazz(3,:) = wdir;               % psd true wind direction, deg
    jazz(4,:) = ts;                 % psd seasnake T, C
    jazz(5,:) = nd;                 % ship theromsalinograph T, C  bow
    jazz(6,:) = nd;                 % ship theromsalinograph salinity, psu  bow
    jazz(7,:) = ta;                 % psd air T at zt, C
    jazz(8,:) = qs;                 % psd air specific humidity at sea surface, g/kg
    jazz(9,:) = qa;                 % psd air specific humidity, g/kg
    jazz(10,:) = rs;                % psd solar flux, w/m^2
    jazz(11,:) = rl;                % psd IR flux, w/m^2
    jazz(12,:) = rain;              % WXT precip rate, mm/hr
    jazz(13,:) = cog;               % cog from gps, m/s
    jazz(14,:) = sog;               % sog from gps, m/s
    jazz(15,:) = hed;               % heading from gps, deg
    jazz(16,:) = rwspd;             % psd WXT rel wind speed, m/s
    jazz(17,:) = rwdir;             % psd WXT rel wind direction, deg
    jazz(18,:) = lat;               % decimal latitude, deg
    jazz(19,:) = lon;               % decimal longtude, deg
    jazz(20,:) = nd;	            % Depth of SST sensor used in heat flux calc, m
    jazz(21,:) = usr;               % u-star
    jazz(22,:) = taub;              % wind stress, coare 3.0, N/m^2
    jazz(23,:) = hsb;               % sensible heat flux, coare 3.5, w/m^2
    jazz(24,:) = hlb;               % latent heat flux, coare 3.5, w/m^2
    jazz(25,:) = rf;                % rain heat flux, w/m^2
    jazz(26,:) = Urf;               % PSD wspd at 3m from COARE
    jazz(27,:) = Trf;               % PSD ta at 3m from COARE
    jazz(28,:) = RHrf;              % PSD rh at 3m from COARE
    
    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,['%11.5f ',repmat('%9.3f ',1,27),' \n'],jazz);
    
    zpath = fullfile(path_ftp_data,[cruise,'_PSD_flux_10min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.zip']);
    zip(zpath,fpath);
    
    fclose all;
    clear jazz;
    
    %% 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_10min + double(2.5/1440); % 5-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);
    Rscl = sw;
    
    %% more plots
    
    if plotit
        if ~isempty(rlclr(isfinite(rlclr)))
            ii = find(isfinite(rlclr));
            figure; plot(jd_10min, pir1,'b',jd_10min, pir2,'g',jd_10min(ii),rlclr(ii),'r'); xlim([ddd ddd+1])
            title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD/SCS 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','BestOutside');
            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_raw_images,[cruise,'_IRflux_Comparison_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
                print(graphdevice,ppath);
            end
        end
        
        if ~isempty(Rscl(isfinite(Rscl)))
            ii = find(isfinite(Rscl));
            figure; plot(jd_10min, psp1,'b', jd_10min, psp2,'g',jd_10min(ii), Rscl(ii),'r'); xlim([ddd ddd+1]);
            title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD/SCS 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('PSD', 'SCS','Clearsky','Location','BestOutside');
            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_raw_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
    disp('AUTOflux Run Successful');
catch BadRun
    disp('Abnormal Termination');
    disp(BadRun.identifier);
    disp(getReport(BadRun));
end

close('all');
fclose('all');
disp(['END of AUTOflux_eval_',cruise,' yearday ',sprintf('%03i',ddd),'.']);
diary off;
zpath = fullfile(path_ftp_data,['MlabLog_',sprintf('%03u',ddd),'.zip']);
zip(zpath,logfname);
exit;

