%{
 STRATUS 2017 :: 2017-03 :: BWB
 Processes daily met files, make daily diagnostic plots and save intermediate
 data files. Averages to 10 min

 Calls:
  read_met_day:     PSD met data from loggers 1, 2, and WXT520 @ 1 min
  read_gps_day:     PSD gps course & speed @ 1 Hz
  read_hed_day:     PSD gps heading & roll @ 1 Hz
  read_motion_day:  PSD MotionPak @ 1 Hz
  read_sonic_day:   PSD sonic wind @ 1 Hz
  read_licor_day:   PSD licor fast water vapor / CO2 @ 1 Hz
  read_scs_day:     RHB ship data @ 1 Hz
%}

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

% set date range for run, whole STRATUS 2017 cruise is 116:139
jdStart = 122;
jdStop = 122;

% other options
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
fsonic = 10;
flicor = 10;
rotationsonic = false;
sonicmodel = 'WindMasterPro';
cruise = 'STRATUS_2017';    % Cruise name
yyyy = '2017';              % year string
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 matlab script path
restoredefaultpath              % reset matlab path
addpath(genpath(fullfile(path_prog,'flux')));	% add scripts to path
savepath; rehash;

% paths relative to data_drive - first 2 should be preexisting
path_python = fullfile(path_prog,'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);
path_raw_images = fullfile(data_drive,cruise,ship,'flux','Raw_Images');
mkdir(path_raw_images);

% Adjustments & postcalibrations for PSD temperatures
tsea_adj = 0;	% all deg C;
ta_adj = 0;
td1_adj = 0;
tc1_adj = 0;
td2_adj = 0;
tc2_adj = 0;
wxt_t_adj = 0;
adj = [tsea_adj,ta_adj,td1_adj,tc1_adj,td2_adj,tc2_adj,wxt_t_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.

% Track map limits
Lonmin = 269;
Lonmax = 284;
Latmin = -19;
Latmax = 33;
PosLims = [Lonmin,Lonmax,Latmin,Latmax];

% sensor heights - provisional until installation is complete
% as measured on 01/13/2015 for CalWater: waterline to base of tower = 7.16 m.
% PSD wind height (m), top of bow tower
zu = 18.6;
% vaisala air T measurement height (m), bow tower
zt = 17;
% vaisala air q measurement height (m), bow tower
zq = 17;
% PSD pressure sensor height (m), top of wheelhouse
zp = 16.5;
% SCS wind height (m)
zuim = 15;
% SCS air T height (m)
ztim = 15;
% SCS air q height (m)
zqim = 15;
% SCS P sensor height (m)
zpim = 15;
% depth of ship tsgm, tssm
dtsg = 5;

%% loop through days
for ddd = jdStart:jdStop

    [m,d] = yd2md(str2num(yyyy), ddd);
    Vdate = [str2num(yyyy),m,d];
    path_working_ddd = fullfile(path_raw_data,[yyyy(3:4),sprintf('%03i',ddd)]);
    if (exist(path_working_ddd,'dir')==0)
        disp(['PATH ',path_working_ddd,' NOT FOUND']);
        continue    % go to next day if 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(['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 (fullfile(path_prog,'flux'));

    %% read data
    % read PSD 1-min avg bulk met data, metm is 1440x24 array
    metm = read_met_day_STRATUS_2017(path_working_ddd,ddd,yyyy,adj);

    % read PSD gps - gprm is 1Hz 86400x7 array
    gprm = read_gps_day_STRATUS_2017(path_working_ddd,ddd,yyyy,PosLims);

    % read PSD heading-roll - hedm is 1Hz 86400x3 array
    hedm = read_hed_day_STRATUS_2017(path_working_ddd,ddd,yyyy);

    % read PSD motionPak - motm is 1Hz 86400x7
    motm = read_motion_day_STRATUS_2017(path_working_ddd,ddd,yyyy);
    
    % read PSD sonic - sonm is 1Hz 86400x5
    sonm = read_sonic_day_STRATUS_2017(path_working_ddd,ddd,yyyy,sonicmodel,...
         rotationsonic,fsonic,prtit,path_raw_images,graphformat,graphdevice);
    
    % read PSD Licor - licm is 1Hz 86400x10
    licm = read_licor_day_STRATUS_2017(path_working_ddd,ddd,yyyy,...
                flicor,prtit,path_raw_images,graphformat,graphdevice);
    
    % read ship SCS data - scsm is 1Hz 86400x17
    scsm = read_scs_day_STRATUS_2017(path_working_ddd,ddd,yyyy,PosLims);

    %% compute 1-min averages for PSD & SCS variables
    rdcon = pi/180;

    % compte 1-min averages of some 10 Hz variables
    jd_min = metm(:,1); % ref 1-min timestamp
    delta = double(1.0/1440);
    last = ddd + 1440*delta;
    jd_1bin = (ddd:delta:last);	jd_1bin = jd_1bin'; % ref 1-min bin edges
    % 1-min gps data
    gprm1 = interval_avg(gprm(:,1), gprm(:,2:7), jd_1bin); % average everything
    Ngps_1 = gprm1(:,4); % mean ship speed toward N
    Egps_1 = gprm1(:,5); % mean ship speed toward W
    sog_1 = sqrt(Ngps_1.^2 + Egps_1.^2);  % recompute sog from speed components
    cog_1 = atan2(Egps_1, Ngps_1)/rdcon; % recompute cog from speed components
    cog_1 = mod(cog_1 + 360, 360);
    gprm1(:,2) = cog_1; % replace with updated 1-min cog/sog
    gprm1(:,3) = sog_1;
    % 1-min heading
    [hed_n, hed_e] = pol2cart(hedm(:,2)*rdcon,1);  % N/E heading components
    hedm = [hedm, hed_n, hed_e];                   % add to array
    hedm1 = interval_avg(hedm(:,1), hedm(:,2:5), jd_1bin); % average everything
    hed_1 = atan2(hedm1(:,5), hedm1(:,4))/rdcon; % recompute hed from avg N/E components
    hed_1 = mod(hed_1 + 360, 360);
    hedm1(:,2) = hed_1; % replace with updated 1-min hed
    
    scsm1 = interval_avg(scsm(:,1), scsm(:,2:17), jd_1bin); % average everything
    
    licm1 =  interval_avg(licm(:,1), licm(:,2:10), jd_1bin); % average everything

    %% compute 10-min averages
    delta = double(10.0/1440);
    last = ddd + 144*delta;
    jd_10bin = (ddd:delta:last);	jd_10bin = jd_10bin'; % ref 10-min bin edges

    % 10-min SCS
    %{
    1 = jday    scsm column IDs
    2 = latm
    3 = lonm
    4 = lrgm
    5 = cogm
    6 = sogm
    7 = imrh
    8 = imta
    9 = imum
    10 = imdm
    11 = impress
    12 = imrspd
    13 = imrdir
    14 = tsgm
    15 = tssm
    %}
    [lrg_1n, lrg_1e] = pol2cart(scsm(:,12)*rdcon,1);         % N/E heading components
    [SOG_1n, SOG_1e] = pol2cart(scsm(:,4)*rdcon, scsm(:,5)); % N/E ship velocity components
    scsm = [scsm, lrg_1n, lrg_1e, SOG_1n, SOG_1e];           % add to scsm array, now 86400x21

    scsm10 = interval_avg(scsm(:,1), scsm(:,2:21), jd_10bin); % average everything
    scsm10_std = interval_std(scsm(:,1), scsm(:,2:21), jd_10bin); % std dev of everything
    jd_10min = scsm10(:,1); % reference10-min timestamp
    latm_10 = scsm10(:,2);
    lonm_10 = scsm10(:,3);
    lrg_10n = scsm10(:,18);
    lrg_10e = scsm10(:,19);
    lrgm_10 = atan2(lrg_10e, lrg_10n)/rdcon; % from avg components
    lrgm_10 = mod(lrgm_10 + 360, 360);
    SOG_10n = scsm10(:,20);
    SOG_10e = scsm10(:,21);
    cogm_10 = atan2(SOG_10e, SOG_10n)/rdcon; % from avg components
    cogm_10 = mod(cogm_10 + 360, 360);
    sogm_10 = sqrt(SOG_10e.^2 + SOG_10n.^2);  % from avg components
    sogm_10_std = scsm10_std(:,5);    % std dev in ship speed
    imrh_10 = scsm10(:,8);
    imta_10 = scsm10(:,9);
    imum_10 = scsm10(:,13);
    % recompute wdir
    imdm_10 = interval_avg(scsm(:,1), unwrap(scsm(:,14)*rdcon), jd_10bin);
    imdm_10 = imdm_10(:,2);
    imdm_10 = mod(imdm_10/rdcon + 360*6, 360);
    impress_10 = scsm10(:,16)+0.125*(zpim); % @ sea level
    tsgm_10 = scsm10(:,6);
    tssm_10 = scsm10(:,7);
%     ts38_10 = scsm10(:,6);
    imsol_10 = scsm10(:,10);
    imir_10 = scsm10(:,15);
    imqa_10 = qair_p([impress_10, imta_10, imrh_10]);
    % compute relative winds
    [imrspd_10,imrdir_10,~,~] = sd_true_to_uv_rel(imum_10,imdm_10,lrgm_10,cogm_10,sogm_10);
    imrdir_10 = mod(imrdir_10 + 720, 360);
    imrdir_10(imrdir_10>180) = imrdir_10(imrdir_10>180) - 360; % +/- 180 deg format

    % 10-min PSD gps
    gprm10 = interval_avg(gprm(:,1), gprm(:,2:7), jd_10bin); % average everything
    Ngps_10 = gprm10(:,4); % ship speed toward N
    Egps_10 = gprm10(:,5); % ship speed toward W
    lat_10 = gprm10(:,6);
    lon_10 = gprm10(:,7);
    sog_10 = sqrt(Ngps_10.^2 + Egps_10.^2);  % sog from speed components
    cog_10 = atan2(Egps_10, Ngps_10)/rdcon; % cog from speed components
    cog_10 = mod(cog_10 + 360, 360);

    % 10-min PSD heading
    hedm10 = interval_avg(hedm(:,1), hedm(:,2:5), jd_10bin); % average everything
    hed_10n = hedm10(:,4);
    hed_10e = hedm10(:,5);
    hed_10 = atan2(hed_10e, hed_10n)/rdcon; % recompute hed from avg N/E components
    hed_10 = mod(hed_10 + 360, 360);
    hedm10(:,2) = hed_10; % replace with updated 10-min hed

    % 10-min PSD met & radiation
    %{
    1   jd_met          1 min timestamp
    2   Tvais           T vaisala, C
    3   Rhvais          RH vaisala
    4   Tsea            sea snake surface SST
    5   psp1            psp1 radiation, W/m2
    6   pir1            pir1 radiation, W/m2
    7   psp2            psp2 radiation, W/m2
    8   pir2            pir2 radiation, W/m2
    9   Tc1             pir1 case temperature
    10  Td1             pir1 dome temperature
    11  Tc2             pir2 case temperature
    12  Td2             pir2 dome temperature
    13  org             rain, mm/hr - none for this cruise
    14  press           air pressure, mb
    15  aspir_on        T/RH diagnostic
    16  org_carrier     org diagnostic - none for this cruise
    17  rwdir_wxt       rel wind direction ,0-360 deg from bow
    18  rwspd_wxt       rel wind speed, m/s
    19  Tair_wxt        air temperature, C
    20  RH_wxt          RH, %
    21  Pmb_wxt         air pressure, mb or hPa
    22  rain_wxt        rain rate, mm/hr
    23  U_wxt           U wind component, WXT
    24  V_wxt           V wind component, WXT
    %}

    metm10 = interval_avg(metm(:,1), metm(:,2:24), jd_10bin);
    Tvais_10 = metm10(:,2);
    Rhvais_10 = metm10(:,3);
    Tsea_10 = metm10(:,4);
    psp1_10 = metm10(:,5);
    pir1_10 = metm10(:,6);
    psp2_10 = metm10(:,7);
    pir2_10 = metm10(:,8);
    Tc1_10 = metm10(:,9);
    Td1_10 = metm10(:,10);
    Tc2_10 = metm10(:,11);
    Td2_10 = metm10(:,12);
    org_10 = metm10(:,13);
    press_10 = metm10(:,14)+0.125*(zp); % @ sea level
    aspir_on_10 = metm10(:,15);
    org_carrier_10 = metm10(:,16);
    qair_10 = qair_p([Tvais_10,Rhvais_10,press_10]);
    qsea_10 = qsea_p([Tsea_10,press_10]);
    
    % compute 10-min PSD winds from sonic
    sonm10 = interval_avg(sonm(:,1), sonm(:,2:end), jd_10bin); % average everything
    % compute wind spd/dir - apply flow distortion corrections
    [rwspd_10,rwdir_10] = uv_to_sd(sonm10(:,2)/0.95,sonm10(:,3)/1.15);
    rwdir_10(rwdir_10>180) = rwdir_10(rwdir_10>180)-360;    % +/- 180 deg format
    [wspd_10,wdir_10,~,~] = uv_rel_to_sd_true(sonm10(:,2)/0.95,sonm10(:,3)/1.15,hed_10,cog_10,sog_10);
    
    %10-min psd licor
    licm10 = interval_avg(licm(:,1), licm(:,2:end), jd_10bin);
    
    % recompute WXT winds from U,V components
    [rwspd_wxt_10,rwdir_wxt_10] = uv_to_sd(metm10(:,23),metm10(:,24));
    rwdir_wxt_10(rwdir_wxt_10>180) = rwdir_wxt_10(rwdir_wxt_10>180) - 360; % +/- 180 format
    [wspd_wxt_10,wdir_wxt_10,~,~] = uv_rel_to_sd_true(metm10(:,23),metm10(:,24),hed_10,cog_10,sog_10);

    Tair_wxt_10 = metm10(:,19);
    RH_wxt_10 = metm10(:,20);
    Pmb_wxt_10 = metm10(:,21)+0.125*(zu); % @ sea level
    rain_wxt_10 = metm10(:,22);

 
    %% Some plots...

    if plotit
        % PRESSURE
        figure; plot(jd_10min,press_10,'r-',jd_10min,impress_10,'g-',jd_10min,Pmb_wxt_10,'b-','linewidth',1); grid;
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Pressures',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        xlabel('Hour (UTC)'); ylabel('3m Pressure (mb)');
        axis([ddd ddd+1 median(press_10(~isnan(press_10)))-10  median(press_10(~isnan(press_10)))+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));
        legend('PSD','SCS','WXT','location','BestOutside');
        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_10,'r',jd_10min,Td1_10,'y',jd_10min,Tc2_10,'m',jd_10min,Td2_10,'k',jd_10min,Tair_wxt_10,'b',jd_10min,Tvais_10,'g.',jd_10min,imta_10,'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','WXT','Tvais','SCS','Location','BestOutside');
        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, Rhvais_10, 'r',jd_10min, RH_wxt_10,'b',jd_10min, imrh_10,'g','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]); legend('PSD','WXT','SCS','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,'_RH_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        % RAIN RATE
        figure; plot(jd_10min, org_10, 'r',jd_10min, rain_wxt_10, 'b','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');
        legend('ORG','WXT','location','BestOutside');
        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

        % PSD T/RH ASPIRATOR FAN
        figure; plot(jd_10min, aspir_on_10,'b','linewidth',1);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD T/RH aspirator fan',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlabel('Hour (UTC)');ylabel('Aspirator Fan Speed (rpm)'); xlim([ddd,ddd+1]);
        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,'_backflow_indicator_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        % SST
        figure; plot(jd_10min, Tsea_10,'r',jd_10min, tsgm_10,'b','linewidth',1); grid;
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  SST',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlabel('Hour (UTC)');ylabel('Sea Surface Temperature (degC)'); xlim([ddd ddd+1]);
        legend('tsnk','tsgm','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,'_SST_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        % LAT/LON TRACK PLOT
        figure; plot(lon_10, lat_10, 'b.', lonm_10, latm_10, 'r.','linewidth',1); grid;
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  GPS Lat/Lon',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlabel('Longitude (deg E)'); ylabel('Latitude (deg N)');
        legend('PSD','SCS','location','BestOutside'); axis([Lonmin Lonmax Latmin Latmax]);
        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,'_GPS_track_LAT_LON_',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_10, jd_10min, sogm_10,'linewidth',1);
        xlabel('Hour (UTC)');ylabel('SOG (m/s)'); legend('PSD','SCS','location','BestOutside'); 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_10,'b.', jd_10min, cogm_10,'r.'); legend('PSD','SCS','location','BestOutside');
        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_10,'b.', jd_10min, lrgm_10, 'r.');
        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)'); legend('PSD','SCS','location','BestOutside');
        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

        % Roll ANGLE
        roll = hedm(:,3); % this is 1Hz
        ii = find(isnan(roll)); roll(ii) = nanmean1(roll);
        figure; plot(hedm(:,1),roll,'b-'); grid;
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Roll Angle ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlabel('Hour (UTC)'); ylabel('Roll (deg)'); 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,'_Roll_',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_10,'r.',jd_10min, rwdir_wxt_10,'b.', jd_10min, imrdir_10, 'g.',...	% 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;
        legend('WXT','SCS','location','BestOutside');
        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_10, 'r',jd_10min, rwspd_wxt_10, 'b', jd_10min, imrspd_10, 'g','linewidth',1);
        legend('PSD','WXT','SCS','location','BestOutside'); 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_10,'r.',jd_10min, wdir_wxt_10,'b.', jd_10min, imdm_10, 'g.','linewidth',1); grid;
        xlabel('Hour (UTC)'); ylabel('Wind Direction (deg)'); axis([ddd ddd+1 -5 365]);
        legend('PSD','WXT','SCS','location','Best');
        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_10, 'r',jd_10min, wspd_wxt_10, 'b', jd_10min, imum_10, 'g','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


        % MOTION PLOT - 1 Hz
        figure;
        subplot(3,2,1);plot(motm(:,1),motm(:,2),'.','markersize',2);ylabel('accx (m.s^-^2)','FontSize',8);
        ylim([-5 5]);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));
        grid;
        subplot(3,2,2);plot(motm(:,1),motm(:,3),'.','markersize',2);ylabel('accy (m.s^-^2)','FontSize',8);
        ylim([-5 5]);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));
        grid;
        subplot(3,2,3);plot(motm(:,1),motm(:,4),'.','markersize',2);ylabel('accz (m.s^-^2)','FontSize',8);
        ylim([5 15]);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));
        grid;
        subplot(3,2,4);plot(motm(:,1),motm(:,5),'.','markersize',2);ylabel('ratex (rad/s)','FontSize',8);
        ylim([-5 5]*.02);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));
        grid;
        subplot(3,2,5);plot(motm(:,1),motm(:,6),'.','markersize',2);ylabel('ratey (rad/s)','FontSize',8);xlabel('Hour (UTC)','FontSize',8);
        ylim([-5 5]*.02);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));
        grid;
        subplot(3,2,6);plot(motm(:,1),motm(:,7),'.','markersize',2);ylabel('ratez (rad/s)','FontSize',8);xlabel('Hour (UTC)','FontSize',8);
        ylim([-5 5]*.02);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));
        grid;
        annotation(gcf,'textbox',[0.25 0.96 0.8 0.02462],'String',sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD motion pack',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','FitBoxToText','off','LineStyle','none','Interpreter','none');
        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,'_MotionPak_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end
    
        % LICOR AGC
        figure; plot(jd_min, licm1(:,7), 'r.','linewidth',1); grid;
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Licor AGC ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlabel('Hour (UTC)'); ylabel('AGC'); xlim([ddd ddd+1]); ylim([0,100]);
        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,'_AGC_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end
    end

    close all;

    %% save 1-min results to text files, historical format
    % with added columns in scs and met files for nav and WXT variables

    % 1-min PSD met data
    % jd pir1 psp1 Tc1 Td1 tsnk ta rh rain org_diag aspir Pmb rwdir_wxt rwspd_wxt ta_wxt rh_wxt Pmb_wxt rain_wxt U_wxt V_wxt
    fpath = fullfile(path_proc_data,[cruise,'_met_1min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
    fout = fopen(fpath,'w');
    jazz = zeros(20,1440)*NaN;
    jazz(1,:) = jd_min;
    jazz(2,:) = metm(:,6);      % pir1 W/m2
    jazz(3,:) = metm(:,5);      % psp1 W/m2
    jazz(4,:) = metm(:,9);      % Tc1 C
    jazz(5,:) = metm(:,10);     % Td1 C
    jazz(6,:) = metm(:,4);      % Tsnake C
    jazz(7,:) = metm(:,2);      % Tvais C
    jazz(8,:) = metm(:,3);      % RHvais
    jazz(9,:) = metm(:,13);     % org mm/hr
    jazz(10,:) = metm(:,16);    % org carrier
    jazz(11,:) = metm(:,15);    % aspir on
    jazz(12,:) = metm(:,14);	% press mb at zp
    jazz(13,:) = metm(:,17);    % WXT rel wind direction ,0-360 deg from bow
    jazz(14,:) = metm(:,18);    % WXT rel wind speed, m/s
    jazz(15,:) = metm(:,19);    % WXT air temperature, C
    jazz(16,:) = metm(:,20);    % WXT RH, %
    jazz(17,:) = metm(:,21);    % WXT air pressure, mb or hPa at zu
    jazz(18,:) = metm(:,22);    % WXT rain rate, mm/hr
    jazz(19,:) = metm(:,23);    % U wind component, WXT
    jazz(20,:) = metm(:,24);    % V wind component, WXT
    fprintf(fout,['%11.5f ',repmat('%9.2f ',1,19),' \n'],jazz);
    fclose all;

    % 1-min radiation data
    % jd pir1 pir2 psp1 psp2 Tc1 Td1 Tc2 Td2 ta rh tsnk imir imsol
    fpath = fullfile(path_proc_data,[cruise,'_rad_1min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
    fout = fopen(fpath,'w');
    jazz = zeros(14,1440)*NaN;
    jazz(1,:) = jd_min;
    jazz(2,:) = metm(:,6);      % pir1 W/m2
    jazz(3,:) = metm(:,8);      % pir2 W/m2
    jazz(4,:) = metm(:,5);      % psp1 W/m2
    jazz(5,:) = metm(:,7);      % psp2 W/m2
    jazz(6,:) = metm(:,9);      % Tc1 C
    jazz(7,:) = metm(:,10);     % Td1 C
    jazz(8,:) = metm(:,11);     % Tc2 C
    jazz(9,:) = metm(:,12);     % Td2 C
    jazz(10,:) = metm(:,2);     % Tvais C
    jazz(11,:) = metm(:,3);     % RHvais
    jazz(12,:) = metm(:,4);     % Tsnake C
    jazz(13,:) = scsm1(:,15);    % imet pir
    jazz(14,:) = scsm1(:,10);    % imet psp
    fprintf(fout,['%11.5f ',repmat('%9.3f ',1,13),' \n'],jazz);
    fclose('all');

    % 1-min gps data
    % jd lat lon sog cog hed roll
    fpath = fullfile(path_proc_data,[cruise,'_gpsnav_1min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
    fout = fopen(fpath,'w');
    jazz = zeros(7,1440)*NaN;
    jazz(1,:) = jd_min;
    jazz(2,:) = gprm1(:,6);     % decimal latitude, deg
    jazz(3,:) = gprm1(:,7);     % decimal longitude, deg
    jazz(4,:) = gprm1(:,3);     % GPS SOG, m/s
    jazz(5,:) = gprm1(:,2);     % GPS COG, deg
    jazz(6,:) = hedm1(:,2);     % Crescent GPS heading, deg
    jazz(7,:) = hedm1(:,3);     % Crescent GPS angle (roll), deg
    fprintf(fout,['%11.5f ',repmat('%9.3f ',1,6),' \n'],jazz);
    fclose('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

    % using PSD data, pressure at sealevel
    y = coare35vn(wspd_10,zu,Tvais_10,zt,Rhvais_10,zq,press_10,Tsea_10,psp1_10,pir1_10,lat_10,600,org_10,NaN,NaN);
    
    % using SCS met data, pressure at sealevel, using PSD rain
    yim = coare35vn(imum_10,zuim,imta_10,ztim,imrh_10,zqim,impress_10,tsgm_10,imsol_10,imir_10,latm_10,600,org_10,NaN,NaN);

    hsb = y(:,1);
    hlb = y(:,2);
    taub = y(:,3);
    tcol = y(:,11);
    qcol = y(:,12);
    rf = y(:,14);
    usr = y(:,8);
    zo = y(:,4);
    Urf = y(:,22);
    Trf = y(:,23);
    RHrf = y(:,24);
    imUrf = yim(:,22);
    imTrf = yim(:,23);
    imRHrf = yim(:,24);

    Rns = psp1_10*0.955;
    Rnl = 0.97*(pir1_10 - 5.67e-8*(Tsea_10+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(metm10(isfinite(metm10(:,2:end))))
        % NET HEAT FLUX
        figure; plot(jd_10min, hnet,'b-','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',1,'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,'b','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

    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');
    jazz = NaN(38,144);
    jazz(1,:) = jd_10min;            % Decimal julian day at beginning of line
    jazz(2,:) = wspd_10;             % psd true wind speed at zu, m/s
    jazz(3,:) = wdir_10;             % psd true wind direction, deg
    jazz(4,:) = Tsea_10;             % psd seasnake T, C
    jazz(5,:) = tsgm_10;             % ship theromsalinograph T, C  bow
    jazz(6,:) = tssm_10;             % ship theromsalinograph salinity, psu  bow
    jazz(7,:) = Tvais_10;            % psd air T at zt, C
    jazz(8,:) = qsea_10;             % psd air specific humidity at sea surface, g/kg
    jazz(9,:) = qair_10;             % psd air specific humidity, g/kg
    jazz(10,:) = psp1_10;            % psd solar flux, w/m^2
    jazz(11,:) = pir1_10;            % psd IR flux, w/m^2
    jazz(12,:) = org_10;             % psd org precip rate, mm/hr
    jazz(13,:) = sogm_10;            % ship sog from gps, m/s
    jazz(14,:) = lrgm_10;            % ship heading from gyrocompass, deg
    jazz(15,:) = rwspd_10;           % psd WXT rel wind speed, m/s
    jazz(16,:) = rwdir_10;           % psd WXT rel wind direction, deg
    jazz(17,:) = latm_10;            % decimal latitude, deg
    jazz(18,:) = lonm_10;            % decimal longtude, deg
    jazz(19,:) = ones(size(jd_10min))*dtsg;	% Depth of SST sensor used in heat flux calc, m
    jazz(20,:) = sogm_10_std;        % standard deviation of ship speed, m/s
    jazz(21,:) = taub;               % wind stress, coare 3.0, N/m^2
    jazz(22,:) = hsb;                % sensible heat flux, coare 3.5, w/m^2
    jazz(23,:) = hlb;                % latent heat flux, coare 3.5, w/m^2
    jazz(24,:) = rf;                 % rain heat flux, w/m^2
    jazz(25,:) = imta_10;            % IMET air temp, C
    jazz(26,:) = imqa_10;            % IMET air specific humidity, g/kg
    jazz(27,:) = imum_10;            % IMET true wind speed, m/s
    jazz(28,:) = imdm_10;            % IMET true wind direction, deg
    jazz(29,:) = imsol_10;           % IMET solar flux, w/m^2
    jazz(30,:) = imir_10;            % IMET IR flux, w/m^2
    jazz(31,:) = impress_10;         % IMET BP, mb @ SL
    jazz(32,:) = imrh_10;            % IMET RH (%)
    jazz(33,:) = Urf;               % PSD wspd at 3m from COARE
    jazz(34,:) = Trf;               % PSD ta at 3m from COARE
    jazz(35,:) = RHrf;              % PSD rh at 3m from COARE
    jazz(36,:) = imUrf;             % SCS wspd at 3m from COARE
    jazz(37,:) = imTrf;             % SCS ta at 3m from COARE
    jazz(38,:) = imRHrf;            % SCS rh at 3m from COARE
    fprintf(fout,['%11.5f ',repmat('%9.3f ',1,37),' \n'],jazz);
    fclose all;
    clear jazz;

    %% radiation models @ 10 min resolution
    % Downward IR
    rlclr_10 = (0.52+0.13/60*abs(latm_10)+(0.082-0.03/60.*abs(latm_10)).*sqrt(qair_10)).*(5.67e-8*(Tvais_10+273.15).^4);

    % Solar model
    latd = nanmean1(latm_10);
    lond = nanmean1(lonm_10);
    qa = qair_10; % humidity needed for solar model
    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(press_10),k1,k2,oz);
    Rscl_10 = sw;

    %% more plots

    if plotit
        if ~isempty(rlclr_10(isfinite(rlclr_10)))
            ii = find(isfinite(rlclr_10));
            figure; plot(jd_10min, pir1_10,'b',jd_10min, imir_10,'g',jd_10min(ii),rlclr_10(ii),'r','linewidth',1); 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('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,'_IRflux_Comparison_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
                print(graphdevice,ppath);
            end
        end

        if ~isempty(Rscl_10(isfinite(Rscl_10)))
            ii = find(isfinite(Rscl_10));
            figure; plot(jd_10min, psp1_10,'b', jd_10min, imsol_10,'g',jd_10min(ii), Rscl_10(ii),'r','linewidth',1); 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

    close('all');
    fclose('all');
    disp(['END of MANUALflux_eval_',cruise,' yearday ',sprintf('%03i',ddd),'.']);

end
