function flag = AUTOflux_eval_CAPRICORN_2016(ddd)
%{
 AUTOflux_eval_CAPRICORN_2016.m :: 2016-03 :: BWB
 Processes daily flux/met files, saves daily dianostic plots and
 intermediate data files, averaged to 1, 5 and 30 minutes.

 Calls:
  read_scs_day:     Ship data
  read_met_day:     PSD 1-min met data loggers 1 and 2
  read_gps_day:     PSD gps course & speed @ 1 Hz
  read_hed_day:     PSD gps heading & roll @ 10 Hz
  read_son_day:     PSD sonic uvw and Ts @ 10 Hz
  read_motion_day:  PSD MotionPak @ 10 Hz
  read_licor_day:   PSD licor water vapor @ 10 Hz
%}

cruise = 'CAPRICORN_2016';  % Cruise name, e.g. 'STRATUS_2009', etc...
yyyy = '2016';              % year string
ship = 'Investigator';      % Research vessel
data_drive = 'D:Data';      % define base path - DAS system
% data_drive = '/Users/byronb/Raw_Data_Files.noindex';   % base path - mac
% Lacie_drive = 'G:';
% Lacie_drive = '/Users/byronb/Raw_Data_Files.noindex/Lacie'; % test mac directory

% paths relative to data_drive - these 4 paths should be preexisting
path_raw_data = fullfile(data_drive,cruise,ship,'flux','Raw');
% Investigator 5 sec underway data is single netCDF file for entire cruise
path_scs_data = fullfile(data_drive,cruise,ship,'flux','SCS','in2016_v02uwy.nc');
path_proc_data = fullfile(data_drive,cruise,ship,'flux','Processed_Raw');
path_raw_images = fullfile(data_drive,cruise,ship,'flux','Raw_Images');
% path_backup_raw = fullfile(Lacie_drive,cruise,ship,'flux','Raw');
% mkdir(path_backup_raw); % make backup directory if it doesn't exist
% path_backup_proc = fullfile(Lacie_drive,cruise,ship,'flux','Processed');
% mkdir(path_backup_proc);
% path_backup_rawIm = fullfile(Lacie_drive,cruise,ship,'flux','Raw_Images');
% mkdir(path_backup_rawIm);
% path_backup_prog = fullfile(Lacie_drive,cruise,ship,'Scientific_analysis','programs');
% mkdir(path_backup_prog);

% construct data path variables
[m,d] = yd2md(str2double(yyyy), ddd);
Vdate = [str2double(yyyy),m,d];
path_working_ddd = fullfile(path_raw_data,[yyyy(3:4),sprintf('%03i',ddd)]);

% Adjustments & postcalibrations for PSD temperatures, all deg C
tsea_adj = -0.5;  % sea snake high by 0.5 deg - 3/13/16 test
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];

% Sonic configuration
% sonicmodel = 'WindMasterPro';	% Sonic model
fsonic = 10;                    % frequency (Hz)
% rotationsonic = false;	    % rotate 30 deg = true
% outputsonic = 1;              % Sonic output. 0=speed of sound, 1=Temperature (degC)

% Licor configuration
flicor = 10;            % frequency (Hz)

% 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.
iv = 1.25;  % approx column water vapor (cm), from local radiometer obs

% Track map limits
Lonmin = 130;
Lonmax = 160;
Latmin = -56;
Latmax = -40;

% psd sensor heights
zu = 23;	% wind speed measurement height (m)
zt = 21.5;	% air T measurement height (m)
zq = 21.5;	% air q measurement height (m)
dtsg = 7;   % depth of ship tsgm, tssm
zp = 16;    % psd pressure
zp_ship = 20; % ship pressure

% other options
plotit = true; % for making plots.
prtit = true;  % for printing plots--plotit must also be true
graphformat = '.png';  % select graphics format files
graphdevice = '-dpng'; % select graphic device

%% load data & average 10 Hz to 1 Hz
% read ship data
scsm = read_scs_day_CAPRICORN_2016(path_scs_data,ddd);
scsm = scsm(:,1:37);  % ignore last column of cum rain

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

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

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

% read PSD sonic - sonm is 1Hz 86400x5 array
sonm = read_sonic_day_CAPRICORN_2016(path_working_ddd,ddd,yyyy,...
         fsonic,prtit,path_raw_images,graphformat,graphdevice);

% read PSD motionPak - motm is 1Hz 86400x7
motm = read_motion_day_CAPRICORN_2016(path_working_ddd,ddd,yyyy);

% read PSD licor for fast water vapor - licm is 1Hz 86400x7 array
licm = read_licor_day_CAPRICORN_2016(path_working_ddd,ddd,yyyy,...
        flicor,prtit,path_raw_images,graphformat,graphdevice);

%% compute 1-min averages
rdcon = pi/180;

% compte 1-min averages
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 psd gps data
gprm1 = interval_avg(gprm(:,1), gprm(:,2:end), 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 psd 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:end), 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

% 1-min psd sonic winds
sonm1 = interval_avg(sonm(:,1), sonm(:,2:end), jd_1bin); % average everything
% compute wind spd/dir - apply flow distortion corrections
[rwspd_1,rwdir_1] = uv_to_sd(sonm1(:,2)/0.95,sonm1(:,3)/1.15);
rwdir_1(rwdir_1>180) = rwdir_1(rwdir_1>180)-360;    % +/- 180 deg format
[wspd_1,wdir_1,~,~] = uv_rel_to_sd_true(sonm1(:,2)/0.95,sonm1(:,3)/1.15,hed_1,cog_1,sog_1);

%1-min psd licor
licm1 = interval_avg(licm(:,1), licm(:,2:end), jd_1bin);

% recompute scs winds and cog/sog from components
[Uport,Vport] = sd_to_uv(scsm(:,7),scsm(:,8));
[Ustbd,Vstbd] = sd_to_uv(scsm(:,11),scsm(:,12));
[Uson,Vson] = sd_to_uv(scsm(:,15),scsm(:,16));
[Urelport,Vrelport] = sd_to_uv(scsm(:,9),scsm(:,10));
[Urelstbd,Vrelstbd] = sd_to_uv(scsm(:,13),scsm(:,14));
[Urelson,Vrelson] = sd_to_uv(scsm(:,17),scsm(:,18));
Nspath = scsm(:,4).*cos(scsm(:,5)./180.*pi);
Espath = scsm(:,4).*sin(scsm(:,5)./180.*pi);
[Nhed, Ehed] = pol2cart(scsm(:,6)*rdcon,1);    % N/E heading components
scsm(:,38:41) = [Nspath, Espath, Nhed, Ehed];  % add N/E components to 5 sec array
    
% compute 1-min avg scs
scsm1 = interval_avg(scsm(:,1), scsm(:,2:end), jd_1bin); % average everything

% recompute cumulative rain
scsm1(:,24) = cumsum(scsm1(:,25))*1/60;

% recompute average velocity & direction components
% first, define array of uv wind & heading componenets

% scs_uv = [scsm(:,1),Uport,Vport,Ustbd,Vstbd,Uson,Vson,Urelport,...
%         Vrelport,Urelstbd,Vrelstbd,Urelson,Vrelson,Nspath,Espath,Nhed,Ehed];

% OR, define array with flow distortion correction to ship winds
scs_uv = [scsm(:,1),Uport,Vport,Ustbd,Vstbd,Uson,Vson,Urelport/0.95,...
        Vrelport/1.15,Urelstbd/0.95,Vrelstbd/1.15,Urelson/0.95,...
        Vrelson/1.15,Nspath,Espath,Nhed,Ehed];
scs_uv1 = interval_avg(scs_uv(:,1), scs_uv(:,2:end), jd_1bin);

% compute true winds
% [scsm1(:,7),scsm1(:,8)] = uv_to_sd(scs_uv1(:,2),scs_uv1(:,3));     % port true winds
% scsm1(:,8) = mod(scsm1(:,8) + 360, 360);
% [scsm1(:,11),scsm1(:,12)] = uv_to_sd(scs_uv1(:,4),scs_uv1(:,5));   % stbd true winds
% scsm1(:,12) = mod(scsm1(:,12) + 360, 360);
% [scsm1(:,15),scsm1(:,16)] = uv_to_sd(scs_uv1(:,6),scs_uv1(:,7));   % son true winds
% scsm1(:,16) = mod(scsm1(:,16) + 360, 360);

% compute relative winds and heading,sog,cog
[scsm1(:,9),scsm1(:,10)] = uv_to_sd(scs_uv1(:,8),scs_uv1(:,9));    % port rel winds
zz = scsm1(:,10)>180; scsm1(zz,10) = scsm1(zz,10)-360;             % +/-180 deg format
[scsm1(:,13),scsm1(:,14)] = uv_to_sd(scs_uv1(:,10),scs_uv1(:,11)); % stbd rel winds
zz = scsm1(:,14)>180; scsm1(zz,14) = scsm1(zz,14)-360;             % +/-180 deg format
[scsm1(:,17),scsm1(:,18)] = uv_to_sd(scs_uv1(:,12),scs_uv1(:,13)); % son rel winds
zz = scsm1(:,18)>180; scsm1(zz,18) = scsm1(zz,18)-360;             % +/-180 deg format
scsm1(:,4) = sqrt(scs_uv1(:,14).^2 + scs_uv1(:,15).^2); % recompute sog from speed components
scsm1(:,5) = atan2(scs_uv1(:,15), scs_uv1(:,14))/rdcon; % recompute cog from speed components
scsm1(:,5) = mod(scsm1(:,5) + 360, 360);
scsm1(:,6) = atan2(scs_uv1(:,17), scs_uv1(:,16))/rdcon; % recompute heading from components
scsm1(:,6) = mod(scsm1(:,6) + 360, 360);

% OR scs true winds from flow distortion corrected relative winds
[scsm1(:,7),scsm1(:,8)] = uv_rel_to_sd_true(scs_uv1(:,8),scs_uv1(:,9),scsm1(:,6),scsm1(:,5),scsm1(:,4)); % port
[scsm1(:,11),scsm1(:,12)] = uv_rel_to_sd_true(scs_uv1(:,10),scs_uv1(:,11),scsm1(:,6),scsm1(:,5),scsm1(:,4)); % stbd
[scsm1(:,15),scsm1(:,16)] = uv_rel_to_sd_true(scs_uv1(:,12),scs_uv1(:,13),scsm1(:,6),scsm1(:,5),scsm1(:,4)); % sonic

%% compute 5-min averages
delta = double(5.0/1440);
last = ddd + 288*delta;
jd_5bin = ddd:delta:last; jd_5bin = jd_5bin(:); % ref 5-min bin edges

% 5-min SCS
scsm5 = interval_avg(scsm(:,1), scsm(:,2:end), jd_5bin);     % average everything
% recompute cumulative rain
scsm5(:,24) = cumsum(scsm5(:,25))*5/60;
scsm5_std = interval_std(scsm(:,1), scsm(:,2:end), jd_5bin); % std dev of everything
scsm_uv5 = interval_avg(scs_uv(:,1), scs_uv(:,2:end), jd_5bin); 
jd_5min = scsm5(:,1); % reference 5-min timestamp
latm_5 = scsm5(:,2);
lonm_5 = scsm5(:,3);
lrg_5n = scsm5(:,40);
lrg_5e = scsm5(:,41);
lrgm_5 = atan2(lrg_5e, lrg_5n)/rdcon; % from avg components
lrgm_5 = mod(lrgm_5 + 360, 360);
SOG_5n = scsm5(:,38);
SOG_5e = scsm5(:,39);
cogm_5 = atan2(SOG_5e, SOG_5n)/rdcon; % from avg components
cogm_5 = mod(cogm_5 + 360, 360);
sogm_5 = sqrt(SOG_5e.^2 + SOG_5n.^2); % from avg components
sogm_5_std = scsm5_std(:,4);    % std dev in ship speed

imta_5p = scsm5(:,19);
imta_5s = scsm5(:,20);
imrh_5p = scsm5(:,21);
imrh_5s = scsm5(:,22);
impress_5 = scsm5(:,23);
imrain_5 = scsm5(:,25);
% imrain_opt_5 = scsm5(:,27);
imsol_5p = scsm5(:,28);
imsol_5s = scsm5(:,29);
imir_5p = scsm5(:,30);
imir_5s = scsm5(:,31);
tsir_5 = scsm5(:,32);
tsgm_5 = scsm5(:,33);
tssm_5 = scsm5(:,34);
imqa_5p = qair3([impress_5, imta_5p, imrh_5p]);
imqa_5s = qair3([impress_5, imta_5s, imrh_5s]);

% recompute winds from components
[imum_5p,imdm_5p] = uv_to_sd(scsm_uv5(:,2),scsm_uv5(:,3));
imdm_5p = mod(imdm_5p + 360, 360);
[imum_5s,imdm_5s] = uv_to_sd(scsm_uv5(:,4),scsm_uv5(:,5));
imdm_5s = mod(imdm_5s + 360, 360);
[imum_5u,imdm_5u] = uv_to_sd(scsm_uv5(:,6),scsm_uv5(:,7));
imdm_5u = mod(imdm_5u + 360, 360);

[imrspd_5p,imrdir_5p] = uv_to_sd(scsm_uv5(:,8),scsm_uv5(:,9));
zz = imrdir_5p>180; imrdir_5p(zz) = imrdir_5p(zz) - 360;
[imrspd_5s,imrdir_5s] = uv_to_sd(scsm_uv5(:,10),scsm_uv5(:,11));
zz = imrdir_5s>180; imrdir_5s(zz) = imrdir_5s(zz) - 360;
[imrspd_5u,imrdir_5u] = uv_to_sd(scsm_uv5(:,12),scsm_uv5(:,13));
zz = imrdir_5u>180; imrdir_5u(zz) = imrdir_5u(zz) - 360;

% 5-min PSD met & radiation
metm5 = interval_avg(metm(:,1), metm(:,2:end), jd_5bin);
Tvais_5 = metm5(:,2);
Rhvais_5 = metm5(:,3);
Tsea_5 = metm5(:,4);
psp1_5 = metm5(:,5);
pir1_5 = metm5(:,6);
psp2_5 = metm5(:,7);
pir2_5 = metm5(:,8);
Tc1_5 = metm5(:,9);
Td1_5 = metm5(:,10);
Tc2_5 = metm5(:,11);
Td2_5 = metm5(:,12);
org_5 = metm5(:,13);
press_5 = metm5(:,14);
aspir_on_5 = metm5(:,15);
org_carrier_5 = metm5(:,16);
qair_5 = qair3([press_5,Tvais_5,Rhvais_5]);
qsea_5 = qsea(Tsea_5);

% 5-min PSD gps
gprm5 = interval_avg(gprm(:,1), gprm(:,2:end), jd_5bin); % average everything
Ngps_5 = gprm5(:,4); % ship speed toward N
Egps_5 = gprm5(:,5); % ship speed toward W
lat_5 = gprm5(:,6);
lon_5 = gprm5(:,7);
sog_5 = sqrt(Ngps_5.^2 + Egps_5.^2);  % sog from speed components
cog_5 = atan2(Egps_5, Ngps_5)/rdcon;  % cog from speed components
cog_5 = mod(cog_5 + 360, 360);

% 5-min PSD heading
hedm5 = interval_avg(hedm(:,1), hedm(:,2:end), jd_5bin); % average everything
hed_5n = hedm5(:,4);
hed_5e = hedm5(:,5);
hed_5 = atan2(hed_5e, hed_5n)/rdcon; % recompute hed from avg N/E components
hed_5 = mod(hed_5 + 360, 360);
hedm5(:,2) = hed_5;

% 5-min Licor
licm5 = interval_avg(licm(:,1), licm(:,2:end), jd_5bin);
qLic_5 = licm5(:,10);
Lic_T_5 = licm5(:,5);
Lic_P_5 = licm5(:,6)*10;  % kPa to mb
Lic_AGC_5 = licm5(:,7);

% 5-min Sonic Winds
sonm5 = interval_avg(sonm(:,1), sonm(:,2:end), jd_5bin);
U_5 = sonm5(:,2)/0.95; % flow distortion corrected rel wind components
V_5 = sonm5(:,3)/1.15;
W_5 = sonm5(:,4);      % mean W
Ts_5 = sonm5(:,5);     % mean Tsonic

[~,rwdir_5] = uv_to_sd(U_5,V_5);
zz = rwdir_5>180; rwdir_5(zz) = rwdir_5(zz)  -360; % +/- 180 deg format
rwspd_5 = sqrt(U_5.^2 + V_5.^2 + W_5.^2);

% compute true wind
[wspd_5,wdir_5,Un_5,Ue_5] = uv_rel_to_sd_true(U_5,V_5,hed_5,cog_5,sog_5);

%% compute 30-min averages
delta = double(30.0/1440);
last = ddd + 48*delta;
jd_30bin = (ddd:delta:last); jd_30bin = jd_30bin(:); % ref 30-min bin edges

% 30-min SCS
scsm30 = interval_avg(scsm(:,1), scsm(:,2:end), jd_30bin);     % average everything
% recompute cumulative rain
scsm30(:,24) = cumsum(scsm30(:,25))*30/60;
scsm30_std = interval_std(scsm(:,1), scsm(:,2:end), jd_30bin); % std dev of everything
scsm_uv30 = interval_avg(scs_uv(:,1), scs_uv(:,2:end), jd_30bin); 
jd_30min = scsm30(:,1);
latm_30 = scsm30(:,2);
lonm_30 = scsm30(:,3);
lrg_30n = scsm30(:,40);
lrg_30e = scsm30(:,41);
lrgm_30 = atan2(lrg_30e, lrg_30n)/rdcon;  % from avg components
lrgm_30 = mod(lrgm_30 + 360, 360);
SOG_30n = scsm30(:,38);
SOG_30e = scsm30(:,39);
cogm_30 = atan2(SOG_30e, SOG_30n)/rdcon;  % from avg components
cogm_30 = mod(cogm_30 + 360, 360);
sogm_30 = sqrt(SOG_30e.^2 + SOG_30n.^2);  % from avg components
sogm_30_std = scsm30_std(:,6);     % std dev in ship speed

imta_30p = scsm30(:,19);
imta_30s = scsm30(:,20);
imrh_30p = scsm30(:,21);
imrh_30s = scsm30(:,22);
impress_30 = scsm30(:,23);
imrain_30 = scsm30(:,25);
% imrain_opt_30 = scsm30(:,27);
imsol_30p = scsm30(:,28);
imsol_30s = scsm30(:,29);
imir_30p = scsm30(:,30);
imir_30s = scsm30(:,31);
tsir_30 = scsm30(:,32);
tsgm_30 = scsm30(:,33);
tssm_30 = scsm30(:,34);
imqa_30p = qair3([impress_30, imta_30p, imrh_30p]);
imqa_30s = qair3([impress_30, imta_30s, imrh_30s]);

% recompute winds from components
[imum_30p,imdm_30p] = uv_to_sd(scsm_uv30(:,2),scsm_uv30(:,3));
imdm_30p = mod(imdm_30p/rdcon + 360, 360);
[imum_30s,imdm_30s] = uv_to_sd(scsm_uv30(:,4),scsm_uv30(:,5));
imdm_30s = mod(imdm_30s/rdcon + 360, 360);
[imum_30u,imdm_30u] = uv_to_sd(scsm_uv30(:,6),scsm_uv30(:,7));
imdm_30u = mod(imdm_30u/rdcon + 360, 360);

[imrspd_30p,imrdir_30p] = uv_to_sd(scsm_uv30(:,8),scsm_uv30(:,9));
zz = imrdir_30p>180; imrdir_30p(zz) = imrdir_30p(zz) - 360;
[imrspd_30s,imrdir_30s] = uv_to_sd(scsm_uv30(:,10),scsm_uv30(:,11));
zz = imrdir_30s>180; imrdir_30s(zz) = imrdir_30s(zz) - 360;
[imrspd_30u,imrdir_30u] = uv_to_sd(scsm_uv30(:,12),scsm_uv30(:,13));
zz = imrdir_30u>180; imrdir_30u(zz) = imrdir_30u(zz) - 360;

% 30-min PSD gps
gprm30 = interval_avg(gprm(:,1), gprm(:,2:end), jd_30bin);
Ngps_30 = gprm30(:,4); % ship speed toward N
Egps_30 = gprm30(:,5); % ship speed toward W
lat_30 = gprm30(:,6);
lon_30 = gprm30(:,7);
sog_30 = sqrt(Ngps_30.^2 + Egps_30.^2);  % sog from speed components
cog_30 = atan2(Egps_30, Ngps_30)/rdcon;  % cog from speed components
cog_30 = mod(cog_30 + 360, 360);

% 30-min PSD heading
hedm30 = interval_avg(hedm(:,1), hedm(:,2:end), jd_30bin); % average everything
hed_30n = hedm30(:,4);
hed_30e = hedm30(:,5);
hed_30 = atan2(hed_30e, hed_30n)/rdcon; % recompute hed from avg N/E components
hed_30 = mod(hed_30 + 360, 360);
hedm30(:,2) = hed_30;

% 30-min PSD met & radiation
metm30 = interval_avg(metm(:,1), metm(:,2:end), jd_30bin);
Tvais_30 = metm30(:,2);
Rhvais_30 = metm30(:,3);
Tsea_30 = metm30(:,4);
psp1_30 = metm30(:,5);
pir1_30 = metm30(:,6);
psp2_30 = metm30(:,7);
pir2_30 = metm30(:,8);
org_30 = metm30(:,13);
press_30 = metm30(:,14);
qair_30 = qair3([press_30,Tvais_30,Rhvais_30]);
qsea_30 = qsea(Tsea_30);

% 30-min Sonic Winds
sonm30 = interval_avg(sonm(:,1), sonm(:,2:end), jd_30bin);
U_30 = sonm30(:,2)/0.95; % flow distortion corrected rel wind components
V_30 = sonm30(:,3)/1.15;
W_30 = sonm30(:,4);

[~,rwdir_30] = uv_to_sd(U_30,V_30);
zz = rwdir_30>180; rwdir_30(zz) = rwdir_30(zz)  -360; % +/- 180 deg format
rwspd_30 = sqrt(U_30.^2 + V_30.^2 + W_30.^2);

% compute true wind
[wspd_30,wdir_30,Un_30,Ue_30] = uv_rel_to_sd_true(U_30,V_30,hed_30,cog_30,sog_30);

%% radiation models @ 5 min resolution
% Downwelling IR
rlclr_5 = (0.52+0.13/60*abs(latm_5)+(0.082-0.03/60.*abs(latm_5)).*sqrt(qair_5)).*(5.67e-8*(Tvais_5+273.15).^4);

% Downwelling Solar
latd = nanmean1(latm_5);
lond = nanmean1(lonm_5);
% qa = qair_5;	% specific humidity
% qrat = 4.0;
% watvap = qa/qrat;
watvap = iv;
jdx = jd_5min + 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;
[Rscl_5,sz,saz,dirs,sky] = SolarRadiance(latd,lond,jdx,tutc,watvap,nanmean1(press_5),k1,k2,oz);

%% Some plots...

if plotit
    % PRESSURE
    figure; plot(jd_5min,press_5+zp*0.125,'r',jd_5min,impress_5+zp_ship*0.125,'b'); hold on;
    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('Sea Level Pressure (mb)'); grid;
    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));
    legend('NOAA','Ship','location','best');
    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;

    % AIR TEMPERATURES
    figure; plot(jd_5min,Tc1_5,'r',jd_5min,Td1_5,'b',jd_5min,Tc2_5,'c',jd_5min,Td2_5,'k',jd_5min,Tvais_5,'g.',jd_5min,imta_5p,'m.',jd_5min,imta_5s,'k.');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD temperatures',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    legend('PIR1 Tc','PIR1 Td','PIR2 Tc','PIR2 Td','Ta noaa','Ta ship port','Ta ship stbd','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,'_Air_Temperatures_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % SOLAR RADIATON
    figure; plot(jd_5min,Rscl_5,'k--',jd_5min, psp1_5,'r', jd_5min, psp2_5,'b',jd_5min,imsol_5p,'g',jd_5min,imsol_5s,'m');
    xlabel('Hour (UTC)'); ylabel('Solar Flux (W/m^2)'); legend('NOAA 1','NOAA 2','Ship port','Ship stbd','Location','best');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Spectral Pyranometers',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    axis([ddd ddd+1 -5 1200]); grid; legend('Clearsky','NOAA 1','NOAA 2','Ship port','Ship stbd','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_raw_images,[cruise,'_Eppley_Solarflux_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % IR RADIATION
    figure; plot(jd_5min,rlclr_5,'k--',jd_5min, pir1_5,'r', jd_5min, pir2_5,'b',jd_5min,imir_5p,'g',jd_5min,imir_5s,'m');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Infrared Radiometers',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    xlabel('Hour (UTC)');ylabel('IR Flux (W/m^2)'); xlim([ddd ddd+1])
    set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks'); grid;
    legend('Clearsky','PIR1','PIR2','Ship port','Ship stbd','location','best');
    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,'_Eppley_IRflux_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % RELATIVE HUMIDITY
    figure; plot(jd_5min, Rhvais_5, 'r',jd_5min, imrh_5p,'b',jd_5min, imrh_5s,'g');
    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 (%)');
    axis([ddd ddd+1 0 105]); grid; legend('NOAA','Ship port','Ship stbd','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_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_5min, org_5, 'r',jd_5min, imrain_5, 'b');
    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'); grid;
    legend('NOAA ORG','Ship RM Young','location','best');
    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;
    
    % CUMULATIVE RAIN
    org_cum = cumsum(org_5)*5/60;
    figure; plot(jd_5min, org_cum, 'r-',jd_5min, scsm5(:,24), 'b-',scsm1(:,1), scsm1(:,24), 'b.');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Cumulative Rain ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    xlabel('Hour (UTC)'); ylabel('Rain (mm)'); xlim([ddd ddd+1]);
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks'); grid;
    legend('NOAA ORG','Ship RM Young','location','best');
    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,'_RainCum_',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_5min, aspir_on_5);
    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,'_PSD_TRH_fan_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % PSD ORG CARRIER
    figure; plot(jd_5min, org_carrier_5);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD ORG signal carrier',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    xlabel('Hour (UTC)');ylabel('Rain Gauge Voltage (volt)'); axis([ddd ddd+1 4.0 5.5]);
    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,'_PSD_ORG_carrier_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % SST
    figure; plot(jd_5min, Tsea_5, 'r',jd_5min, tsgm_5, 'b',jd_5min, tsir_5, 'g');
    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('SeaSnake','Ship SBE','Ship IR','location','best');
    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,'_SST_',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_5min, sog_5,'b.-', jd_5min, sogm_5,'r');
    xlabel('Hour (UTC)');ylabel('SOG (m/s)'); legend('NOAA','Ship','location','best'); 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_5min, cog_5,'b.-', jd_5min, cogm_5,'r'); legend('NOAA','Ship','location','best');
    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;

    % MAP PLOT - m_map folder should be in matlab path
    figure;
    m_proj('mercator','longitudes',[Lonmin Lonmax],'latitudes',[Latmin Latmax]);
    m_grid;
    m_coast('linewidth',2,'color',0.5+[0 0 0]);
    m_coast('patch',[.7 .7 .7],'edgecolor','none');
    m_line(lon_5,lat_5,'color','r','linestyle','none','marker','.','markersize',6);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  R/V %s track',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,ship),'FontWeight','Bold','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,'_shiptrack_map_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % HEADING
    figure; plot(jd_5min, hed_5,'b.-', jd_5min, lrgm_5, '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('NOAA','Ship','location','best');
    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); % 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,'_PSD_Roll_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % SONIC TEMPERATURE
    figure; plot(jd_5min, Ts_5); xlabel('Hour (UTC)'); ylabel('Temperature (degC)'); xlim([ddd ddd+1]);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Sonic Air Temperature',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,'_PSD_Sonic_Temperature_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % PSD SONIC WIND SPEEDS
    figure; plot(jd_5min, U_5,'r', jd_5min, V_5,'b', jd_5min, W_5,'g', jd_5min, rwspd_5,'k');
    xlabel('Hour (UTC)'); ylabel('Wind speed (m/s)'); xlim([ddd ddd+1]);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Sonic Relative Wind Components',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    legend('U','V','W','speed','Location','best'); 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,'_PSD_Sonic_Winds_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % RELATIVE WIND SPEED & DIRECTION
    figure; plot(jd_5min,rwdir_5,'r',jd_5min,imrdir_5p,'b',jd_5min,imrdir_5s,'g',jd_5min,imrdir_5u,'m',...
        [ddd ddd+1],[60 60],'k', ...  % 60 degree sector
        [ddd ddd+1],[90 90],'k--', ...% 90 degree sector
        [ddd ddd+1],[-60 -60],'k', ...
        [ddd ddd+1],[-90 -90],'k--');
    xlabel('Hour (UTC)'); ylabel('Relative Wind Direction (deg)'); axis([ddd ddd+1 -180 180]); grid;
    legend('PSD','Ship port','Ship stbd','Ship sonic','+/- 60 deg','+/- 90 deg','location','best');
    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));
    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_Wdir_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;
    
    figure; plot(jd_5min,rwspd_5,'r',jd_5min,imrspd_5p,'b',jd_5min,imrspd_5s,'g',jd_5min,imrspd_5u,'m'); grid;
    xlabel('Hour (UTC)'); ylabel('Relative Wind Speed (m/s)'); xlim([ddd ddd+1]);
    legend('PSD','Ship port','Ship stbd','Ship sonic','location','best');
    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));
    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_Wspd_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;
    
    % GILL SONIC BIAS W/ RELWDIR
    rwdirp = scsm(:,10); rwdirp(rwdirp>180) = rwdirp(rwdirp>180)-360;
    rwdiru = scsm(:,18); rwdiru(rwdiru>180) = rwdiru(rwdiru>180)-360;
    ylo = min(scsm(:,17)-scsm(:,9))-1;
    yhi = max(scsm(:,17)-scsm(:,9))+1;
    figure; plot([0,0],[-15,15],'k-',[-180,180],[0,0],'k-',rwdirp,scsm(:,17)-scsm(:,9),'r.'); grid;
    xlabel('Rel Wdir, port'); ylabel('Delta Wspd, sonic-port (m/s)'); axis([-180,180,ylo,yhi]);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Gill 2D Bias ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','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,'_Gill_Rwspd_Bias_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;
    
    figure; plot([-180,180],[-180,180],'k-',rwdirp,rwdiru,'r.'); grid;
    xlabel('Rel Wdir, port'); ylabel('Rel Wdir, usonic'); axis([-180,180,-180,180]);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Gill 2D Bias ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','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,'_Gill_Rwdir_Bias_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    
    % TRUE WIND SPEED & DIRECTION
    figure; plot(jd_5min, wdir_5,'r',jd_5min,imdm_5p,'b',jd_5min,imdm_5s,'g',jd_5min,imdm_5u,'m'); grid;
    xlabel('Hour (UTC)'); ylabel('Wind Direction (deg)'); axis([ddd ddd+1 0 360]);
    legend('PSD','Ship port','Ship stbd','Ship sonic','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');
    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,'_True_Wdir_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;
    
    figure; plot(jd_5min,wspd_5,'r',jd_5min,imum_5p,'b',jd_5min, imum_5s,'g',jd_5min,imum_5u,'m'); grid;
    legend('PSD','Ship port','Ship stbd','Ship sonic','location','best');
    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));
    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_Wspd_',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,'_PSD_MotionPak_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % SPECIFIC HUMIDITY
    figure; plot(jd_5min, qLic_5, 'b', jd_5min, qair_5,'r');
    xlabel('Hour (UTC)'); ylabel('Specific humidity (g/kg)');
    legend('Licor 7500','Vaisala aspirated','Location','best'); grid;
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Specific Humidity',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    axis([ddd ddd+1 nanmedian1(qair_5(~isnan(qair_5)))-10 nanmedian1(qair_5(~isnan(qair_5)))+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,'_Specific_Humidity_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % LICOR CO2, 1Hz
    co2 = licm(:,4)*1e6;
    figure; plot(licm(:,1), co2,'b');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Licor CO2',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    xlabel('Hour (UTC)');ylabel('CO2 (ppm)'); grid;
    axis([ddd ddd+1 nanmedian1(co2(~isnan(co2)))-60 nanmedian1(co2(~isnan(co2)))+60]);
    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,'_PSD_Licor_CO2_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % LICOR BOX TEMPERATURE
    figure;plot(jd_5min, Lic_T_5);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Licor Temperature ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    xlabel('Hour (UTC)');ylabel('Temperature (degC)'); grid;
    axis([ddd ddd+1 nanmedian1(Lic_T_5(~isnan(Lic_T_5)))-10 nanmedian1(Lic_T_5(~isnan(Lic_T_5)))+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,'_PSD_Licor_Temp_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % LICOR PRESSURE
    figure;plot(jd_5min, Lic_P_5);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Licor Pressure ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    xlabel('Hour (UTC)');ylabel('Pressure (hPa)'); grid;
    axis([ddd ddd+1 nanmedian1(Lic_P_5(~isnan(Lic_P_5)))-10 nanmedian1(Lic_P_5(~isnan(Lic_P_5)))+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,'_PSD_Licor_Pressure_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % LICOR AGC
    figure;plot(jd_5min, Lic_AGC_5,'b.', [ddd ddd+1],[55 55],'r--', [ddd ddd+1],[60 60], 'r');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Licor AGC value ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    axis([ddd ddd+1 40 100]);
    xlabel('Hour (UTC)');ylabel('AGC (%)'); 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,'_PSD_Licor_AGC_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;
    
    % Sea Snake / Air Temp comp - only useful when snake is on deck at night
    kk = find(psp1_5<1); % && jd_5min<x && jd_5min>y);
    Tavg = (Tc1_5+Tc2_5+Tvais_5)/3; % average air temp from PIRs and Vaisala T/RH
    lo = min([Tavg(kk);Tsea_5(kk)])-1;
    hi = max([Tavg(kk);Tsea_5(kk)])+1;
    figure; plot([0,30],[0,30],'k-',Tavg(kk), Tsea_5(kk),'bo'); axis([lo,hi,lo,hi]); grid;
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Snake-Tair Comparison ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    xlabel('Tair(avg), C');ylabel('T snake, C');
    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,'_Snake_Tair_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
        close all;
    end;
end;

%% save 1-min results to text files, historical format

% 1-min SCS met/nav data
fpath = fullfile(path_proc_data,[cruise,'_proc_scs_1min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
fout = fopen(fpath,'w');
jazz = scsm1';
fprintf(fout,['%11.6f ',repmat('%10.4f ',1,40),'\r\n'],jazz);
fclose all;

% 1-min PSD met data
fpath = fullfile(path_proc_data,[cruise,'_proc_met_1min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
fout = fopen(fpath,'w');
jazz = zeros(12,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
fprintf(fout,['%11.6f ',repmat('%10.4f ',1,11),'\r\n'],jazz);
fclose all;

% 1-min radiation data
fpath = fullfile(path_proc_data,[cruise,'_proc_rad_1min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
fout = fopen(fpath,'w');
jazz = zeros(12,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,:) = scsm1(:,30);    % ship pir port W/m2
jazz(7,:) = scsm1(:,31);    % ship pir stbd W/m2
jazz(8,:) = scsm1(:,28);    % ship psp port W/m2
jazz(9,:) = scsm1(:,29);    % ship psp stbd W/m2
jazz(10,:) = metm(:,2);     % Tvais C
jazz(11,:) = metm(:,3);     % RHvais
jazz(12,:) = metm(:,4);     % Tsnake C
fprintf(fout,['%11.6f ',repmat('%10.4f ',1,11),'\r\n'],jazz);
fclose('all');

% 1-min gps data
fpath = fullfile(path_proc_data,[cruise,'_proc_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.6f ',repmat('%10.4f ',1,6),'\r\n'],jazz);
fclose('all');

% 1-min sonic wind results - flow distortion corrected true & relative
fpath = fullfile(path_proc_data,[cruise,'_proc_son_1min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
fout = fopen(fpath,'w');
jazz = zeros(9,1440)*NaN;
jazz(1,:) = jd_min;
jazz(2,:) = sonm1(:,2);      % U
jazz(3,:) = sonm1(:,3);      % V
jazz(4,:) = sonm1(:,4);      % W
jazz(5,:) = sonm1(:,5);      % Tsonic
jazz(6,:) = rwdir_1;         % rel wind dir
jazz(7,:) = rwspd_1;         % rel wind spd
jazz(8,:) = wdir_1;          % true wind dir
jazz(9,:) = wspd_1;          % true wind spd
fprintf(fout,['%11.6f ',repmat('%10.4f ',1,8),'\r\n'],jazz);
fclose('all');

% 1-min licor results
% 1   jd          1-min timestamp
% 2   Licor_H2O   H2O, mmol/m3
% 3   Licor_CO2   CO2, approx mmol/m3
% 4   Licor_CO2a  CO2, approx mol/mol
% 5   Licor_T     Licor box temperature
% 6   Licor_P     Licor box pressure kPa
% 7   AGC         Licor diagnostic
% 8   Licor_mr    H2O, mol/mol
% 9   Rho_air     Air density, kg/m3
% 10  Licor_q     H2O, g/kg
fpath = fullfile(path_proc_data,[cruise,'_proc_lic_1min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
fout = fopen(fpath,'w');
jazz = licm1';
fprintf(fout,['%11.6f ',repmat('%10.4f ',1,9),'\r\n'],jazz);
fclose('all');
clear jazz;

%% timatch_5
% 5-min averaging was done above
% run bulk model with 5 min data
% 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 wbar Cd Ch Ce Cdn_10 Chn_10 Cen_10 ug ];
%     1   2   3   4  5   6  7  8   9  10   11   12  13  14  15  16 17 18    19    20     21   22
% 3/16 added U10N, Cdn_10 to output for CAPRICORN

% choose best radiation & SST for bulk model and heat budget
pir_5 = (pir1_5 + pir2_5)/2;     % use psd pir
% psp_5 = (psp1_5 + psp2_5)/2;    % use psd psp
psp_5 = max(imsol_5p,imsol_5s);  % use max of ship PSPs
% SST_5 = Tsea_5;                 % use snake
SST_5 = tsgm_5;                  % use ship SST

[y,U10N] = coare35vn(wspd_5,zu,Tvais_5,zt,Rhvais_5,zq,press_5,SST_5,psp_5,pir_5,latm_5,600,org_5,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);
Cdn_10 = y(:,19);

Rns = 0.955*psp_5;
Rnl = 0.97*(pir_5 - 5.67e-8*(SST_5+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);

%% 5-min bulk model plots

if plotit && ~isempty(sonm5(isfinite(sonm5(:,2:end)))) && ~isempty(metm5(isfinite(metm5(:,2:end))))
    % NET HEAT FLUX
    figure; plot(jd_5min, hnet); 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,'_PSD_Net_Heat_Flux_5min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % HEAT FLUX COMPONENTS - AIR
    figure; plot(jd_5min, hsb, jd_5min, hlb, jd_5min, Rnl, jd_5min, Rns, jd_5min, rf,'g.','LineWidth',2,'MarkerSize',12);
    hold on; plot(jd_5min, 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).  Heat Flux Components (air)',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,'_PSD_Heat_Flux_Components_5min_air_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % HEAT FLUX COMPONENTS - SEA
    figure; plot(jd_5min, -hsb, jd_5min, -hlb, jd_5min, Rnl, jd_5min, Rns, jd_5min, -rf,'g.','LineWidth',2,'MarkerSize',12);
    hold on; plot(jd_5min, 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).  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,'_PSD_Heat_Flux_Components_5min_sea_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % BULK FRICTION VELOCITY
    figure; plot(jd_5min, usr); xlabel('Hour (UTC)'); ylabel('u_*(m/s)'); xlim([ddd ddd+1]);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Friction Velocity, COARE 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);
        close all;
    end;
end;

%% save 5-min bulk results to txt file - historical format + U10N

fpath = fullfile(path_proc_data,[cruise,'_PSD_flux_5min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
fout = fopen(fpath,'w');
jazz = NaN(34,288);
jazz(1,:) = jd_5min;            % Decimal julian day at beginning of line
jazz(2,:) = wspd_5;             % psd true wind speed at 23 m ht, m/s
jazz(3,:) = wdir_5;             % psd true wind direction, deg
jazz(4,:) = Tsea_5;             % psd seasnake T, C
jazz(5,:) = tsgm_5;             % ship theromsalinograph T, C
jazz(6,:) = tssm_5;             % ship theromsalinograph salinity, psu
jazz(7,:) = Tvais_5;            % psd air T at 18 m, C
jazz(8,:) = qsea_5;             % psd air specific humidity at sea surface, g/kg
jazz(9,:) = qair_5;             % psd air specific humidity at 18 m, g/kg
jazz(10,:) = psp_5;             % psd solar flux, w/m^2
jazz(11,:) = pir_5;             % psd IR flux, w/m^2
jazz(12,:) = org_5;             % psd optical raingage precip rate, mm/hr
jazz(13,:) = sogm_5;            % ship sog from gps, m/s
jazz(14,:) = lrgm_5;            % ship heading from gyrocompass, deg
jazz(15,:) = rwspd_5;           % rel wind speed, m/s
jazz(16,:) = rwdir_5;           % rel wind direction, deg
jazz(17,:) = latm_5;            % decimal latitude, deg
jazz(18,:) = lonm_5;            % decimal longtude, deg
jazz(19,:) = scsm5(:,37);	    % depth of SST sensor used in heat flux calc, m
jazz(20,:) = sogm_5_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.0, w/m^2
jazz(23,:) = hlb;               % latent heat flux, coare 3.0, w/m^2
jazz(24,:) = rf;                % rain heat flux, w/m^2
jazz(25,:) = imta_5p;           % IMET air temp, C (31.4 m)
jazz(26,:) = imqa_5p;           % IMET air specific humidity, g/kg (31.4 m)
jazz(27,:) = imum_5p;           % IMET true wind speed, m/s (25 m)
jazz(28,:) = imdm_5p;           % IMET true wind direction, deg (25 m)
jazz(29,:) = imsol_5p;          % IMET solar flux, w/m^2
jazz(30,:) = imir_5p;           % IMET IR flux, w/m^2
jazz(31,:) = impress_5;         % BP, mb @ 20 m
jazz(32,:) = imrh_5p;           % RH (%) @ 31.4 m
jazz(33,:) = U10N;              % 10m neutral windspeed, COARE 3.5
jazz(34,:) = Cdn_10;            % 10m neutral drag coefficient, COARE 3.5

fprintf(fout,['%11.6f ',repmat('%10.4f ',1,33),'\r\n'],jazz);
fclose all;
clear jazz;

%% timatch_30 - 30-min bulk model run
% 30 min averages were computed above
% run bulk model with 30 min data
% 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 wbar Cd Ch Ce Cdn_10 Chn_10 Cen_10 ug ];
%     1   2   3   4  5   6  7  8   9  10   11   12  13  14  15  16 17 18    19    20     21   22

% choose best radiation & SST for bulk model and heat budget
pir_30 = (pir1_30 + pir2_30)/2;     % use psd pir
% psp_30 = (psp1_30 + psp2_30)/2;    % use psd psp
psp_30 = max(imsol_30p,imsol_30s);  % use max of ship PSPs
% SST_30 = Tsea_30;                  % use snake
SST_30 = tsgm_30;                   % use ship SST

[y,~] = coare35vn(wspd_30,zu,Tvais_30,zt,Rhvais_30,zq,press_30,SST_30,psp_30,pir_30,latm_30,600,org_30,NaN,NaN);

hsb = y(:,1);
hlb = y(:,2);
taub = y(:,3);
tcol = y(:,11);
qcol = y(:,12);
rf = y(:,14);
Lmo = y(:,7);
usr = y(:,8);
tsr = y(:,9);
qsr = y(:,10);
zo = y(:,4);

Rns = 0.955*psp_30;
Rnl = 0.97*(pir_30 - 5.67e-8*(SST_30+273.15).^4);
hnet = Rns + Rnl - hsb - hlb - rf;

%% save 30 min bulk results to txt file

fpath = fullfile(path_proc_data,[cruise,'_PSD_flux_30min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
fout = fopen(fpath,'w');
jazz = zeros(32,48)*NaN;

jazz(1,:) = jd_30min;         % Decimal julian day at beginning of line
jazz(2,:) = wspd_30;          % psd true wind speed at 23 m ht, m/s
jazz(3,:) = wdir_30;          % psd true wind direction, deg
jazz(4,:) = Tsea_30;          % psd seasnake T, C
jazz(5,:) = tsgm_30;          % ship theromsalinograph T, C
jazz(6,:) = tssm_30;          % ship theromsalinograph salinity, psu
jazz(7,:) = Tvais_30;         % psd air T at 18 m, C
jazz(8,:) = qsea_30;          % psd air specific humidity at sea surface, g/kg
jazz(9,:) = qair_30;          % psd air specific humidity at 18 m, g/kg
jazz(10,:) = psp_30;          % psd solar flux, w/m^2
jazz(11,:) = pir_30;          % psd IR flux, w/m^2
jazz(12,:) = org_30;          % psd optical raingage precip rate, mm/hr
jazz(13,:) = sogm_30;         % ship sog from gps, m/s
jazz(14,:) = lrgm_30;         % ship heading from gyrocompass, deg
jazz(15,:) = rwspd_30;        % rel wind speed, m/s
jazz(16,:) = rwdir_30;        % rel wind direction, deg
jazz(17,:) = latm_30;         % decimal latitude, deg
jazz(18,:) = lonm_30;         % decimal longtude, deg
jazz(19,:) = scsm30(:,37);    % depth of SST sensor used in heat flux calc, m
jazz(20,:) = sogm_30_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.0, w/m^2
jazz(23,:) = hlb;             % latent heat flux, coare 3.0, w/m^2
jazz(24,:) = rf;              % rain heat flux, w/m^2
jazz(25,:) = imta_30p;        % IMET air temp, C (31.4 m)
jazz(26,:) = imqa_30p;        % IMET air specific humidity, g/kg (31.4 m)
jazz(27,:) = imum_30p;        % IMET true wind speed, m/s (25 m)
jazz(28,:) = imdm_30p;        % IMET true wind direction, deg (25 m)
jazz(29,:) = imsol_30p;       % IMET solar flux, w/m^2
jazz(30,:) = imir_30p;        % IMET IR flux, w/m^2
jazz(31,:) = impress_30;      % BP, mb @ 20 m
jazz(32,:) = Rhvais_30;       % RH (%) @ 31.4 m

fprintf(fout,['%11.6f ',repmat('%10.4f ',1,31),'\r\n'],jazz);
fclose all;
clear jazz;

disp('END of MANUALflux_eval')
flag = 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


