%{
 SeaState 2015 :: 2015-10 :: BWB
 Processes daily flux files, makes daily dianostic plots and intermediate
 data files. Averages to 5 and 30 minutes and saves plots and output files.

 Calls:
  read_SKQ1min_SeaState_2015:  1-min ship data from the SKQ raw files
  read_met_day_SeaState_2015:  PSD met data loggers 1 and 2
  read_gps_day_SeaState_2015:  PSD gps course & speed @ 10 Hz
  read_hed_day_SeaState_2015:  PSD gps heading & pitch @ 10 Hz
  read_sonic_day_SeaState_2015:  PSD sonic uvw and Ts
  read_motion_day_SeaState_2015:  PSD MotionPak
  read_licor_day_SeaState_2015:  PSD licor water vapor
  read_paroP_day_SeaState_2015:  Paro fast pressure sensors
%}

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

% set 'false' if rerun for ddd - no need to repeat copy & clean
data_copy = false;
% set 'true' to activate ftp to Boulder server - if bandwidth allows
data_ftp = false;
% set 'true' to save data, results and plots to external drive
laCie_backup = false;

cruise = 'SeaState_2015';      % Cruise name, e.g. 'STRATUS_2009', etc...
yyyy = '2015';              % year string
ship = 'Sikuliaq';         % Research vessel
% data_drive = 'D:';   % 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 - first 3 should be preexisting
path_raw_data = fullfile(data_drive,cruise,ship,'flux','Raw');
path_proc_data = fullfile(data_drive,cruise,ship,'flux','Processed');
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);

% matlab script path
% path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs','flux');
% cd(path_prog);                  % cd to SeaState_2015 matlab directory
% restoredefaultpath              % reset matlab path
% addpath(genpath(path_prog));	  % add scripts to path
% savepath; rehash;

% Adjustments & postcalibrations for PSD temperatures
tsea_adj = 0;	% all deg C;
ta_adj = 0;
td1_adj = 0;
tc1_adj = 0;  % set tc_adj based on time in read met
td2_adj = 0;
tc2_adj = 0;
ta337_adj = 0;
adj = [tsea_adj,ta_adj,td1_adj,tc1_adj,td2_adj,tc2_adj,ta337_adj];

% Sampling frequencies
fsonic = 20;
flicor = 10;
fpressure = 10;

% 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.
% Ola's values for Arctic
k1 = 0.001;	% aerosol optical depth, band 1
k2 = 0.001;	% aerosol optical depth, band 2
oz = 0.001;	% column ozone
iv = 0.2;	% column water vapor (cm), if not calculated from obs.


% Track map limits
% SeaState_2015
Lonmin = 190;
Lonmax = 220;
Latmin = 63;
Latmax = 76;

% sensor heights
% measured on 09/28/2015
zu_hi = 16.49;	% tower wind speed measurement height (m)
zu_lo = 10.63;  % bow wind speed measurement height (m)
zt = 15.16;	    % air T measurement height (m)
zq = 15.16;	    % air RH measurement height (m)
zt337 = 15.41;  % Vais 337 T sensor height (m)
zq337 = 15.41;  % Vais 337 RH sensor height (m)
zp_psd = 17;    % pressure sensor heights
zp_ship = 15.41;
zp_bow = 9.33;
zp_tower = 15.02;
dtsg = 6;       % depth of ship tsgm, tssm

% 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
saveit = true;  % for saving data files going to a temp directory
clearit = true; % for clearing data sometimes to save memory

%% copy data from raw directory to working directory & clean data files

% enter date and construct data path variables
if ~exist('ddd','var'); ddd=input('Input yearday to evaluate as ddd : '); end;
[m,d] = yd2md(str2num(yyyy), ddd);
Vdate = [str2num(yyyy),m,d];
path_working_ddd = fullfile(path_raw_data,[yyyy(3:4),sprintf('%03i',ddd)]);
path_raw_ddd = fullfile(data_drive,'Data',[yyyy(3:4),sprintf('%03i',ddd)]);

% copy data to working directories and clean up files
if data_copy
    % transfer from local raw data directory to working directory
    copyfile(path_raw_ddd, path_working_ddd);

    % now fix NaNs in the raw data files...
    % these statements will search every text file in the path and replace
    % NAN with NaN.  Perl language must be installed on the computer
    file_ext = 'txt';
    srch_str = 'NAN';
    repl_str = 'NaN';
    disp('Start cleaning NAN');
    f = replace_str_in_files(path_working_ddd, file_ext, srch_str, repl_str);
    disp('Finished cleaning NAN');

    % convert PSD gps files to gprm format with python script
    % gprm files are easier/faster for matlab to parse
    % python language must be installed on the computer
    cd (path_working_ddd); % cd to data directory
    pyscript = fullfile(path_prog,'parseGpsFiles.py');
    system(['python ',pyscript]); % run python script
    cd(path_prog);  % switch back to flux script directory
    disp('Finished generating gprm files');
end;

%% read data
% read SKQ data for entire day
% scs is 1-min avg data, 1440x46
% column designations:
% jd lat lon sog cog hed roll pitch hed heave lat lon depth Pmb Tair rh
%  1  2   3   4   5   6   7     8    9   10    11  12   13   14  15  16  
% sst45 cond sal sos sst38 rwdir rwspd dir spd tpile rl Tcase Tdome rs par
%  17    18  19  20   21    22    23   24  25   26   27  28    29   30  31
% ct15 flowSpd flowVol flowT flowCond lrg rot odecL odecX sogL sogX dirPort
%  32     33     34      35     36     37  38  39    40    41   42    43
% spdPort dirStbd spdStbd
%   44      45      46
scsm = read_SKQ1min_SeaState_2015(ddd,path_working_ddd);

% read PSD 1-min avg bulk met data, metm is 1440x11 array
% jd tvais rhvais sst psp1 pir1 psp2 pir2 tc1 td1 tc2 td2 org press
%  1   2      3    4   5    6    7    8    9   10  11  12  13  14
% aspir_on org_carrier rh337 t337
%    15        16       17    18
metm = read_met_day_SeaState_2015(path_working_ddd,ddd,yyyy,adj);

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

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

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

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

% read PSD licor for fast water vapor - licm is 1Hz 86400x7 array
licm = read_licor_day_SeaState_2015(path_working_ddd,ddd,yyyy,...
           flicor,prtit,path_raw_images,graphformat,graphdevice);
       
% read Paro fast pressure sensor data, fpsm is 86400x3 array
fpsm = read_paroP_day_SeaState_2015(path_working_ddd,ddd,yyyy,...
         fpressure,prtit,path_raw_images,graphformat,graphdevice);

%% compute 1-min averages and define 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);	% 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
% 1-min sonic winds
sonm1 = interval_avg(sonm(:,1), sonm(:,2:9), jd_1bin); % average everything
% compute relative wind spd/dir
rwdir_hi_1 = unwrap(atan2(sonm1(:,3)/1.15,-sonm1(:,2)/0.95))/rdcon;  % rel wind dir
rwdir_hi_1 = mod(rwdir_hi_1 + 360*10, 360);               % +/- 360 deg format
rwdir_hi_1(rwdir_hi_1>180) = rwdir_hi_1(rwdir_hi_1>180)-360;    % +/- 180 deg format
rwspd_hi_1 = sqrt((sonm1(:,2)/0.95).^2 + (sonm1(:,3)/1.15).^2 + sonm1(:,4).^2);
rwdir_lo_1 = unwrap(atan2(sonm1(:,7)/1.15,-sonm1(:,6)/0.95))/rdcon;  % rel wind dir
rwdir_lo_1 = mod(rwdir_lo_1 + 360*10, 360);               % +/- 360 deg format
rwdir_lo_1(rwdir_lo_1>180) = rwdir_lo_1(rwdir_lo_1>180)-360;    % +/- 180 deg format
rwspd_lo_1 = sqrt((sonm1(:,6)/0.95).^2 + (sonm1(:,7)/1.15).^2 + sonm1(:,8).^2);

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

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

% 5-min SCS
[lrg_1n, lrg_1e] = pol2cart(scsm(:,37)*rdcon,1);            % N/E heading components
[SOG_1n, SOG_1e] = pol2cart(scsm(:,5)*rdcon, scsm(:,4));    % N/E ship velocity components
scsm = [scsm, lrg_1n, lrg_1e, SOG_1n, SOG_1e];              % add to array
%               47      48      49      50
scsm5 = interval_avg(scsm(:,1), scsm(:,2:50), jd_5bin);     % average everything
scsm5_std = interval_std(scsm(:,1), scsm(:,2:50), jd_5bin); % std dev of everything
jd_5min = scsm5(:,1); % reference 5-min timestamp
latm_5 = scsm5(:,2);
lonm_5 = scsm5(:,3);
lonm_5(lonm_5<0) = lonm_5(lonm_5<0)+360;
lrg_5n = scsm5(:,47);
lrg_5e = scsm5(:,48);
lrgm_5 = atan2(lrg_5e, lrg_5n)/rdcon; % from avg components
lrgm_5 = mod(lrgm_5 + 360, 360);
SOG_5n = scsm5(:,49);
SOG_5e = scsm5(:,50);
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
imrh_5 = scsm5(:,16);
imta_5 = scsm5(:,15);
ct15_5 = scsm5(:,32) + 0.9;     % corrected to sea snake
impress_5 = scsm5(:,14);
tsgm_5 = scsm5(:,17);
tssm_5 = scsm5(:,19);
sst38_5 = scsm5(:,21);
imqa_5 = qair3([impress_5, imta_5, imrh_5]);
imsol_5 = scsm5(:,30);
% imir_5 = scsm5(:,27);  % recompute below
imtcase_5 = scsm5(:,28);
imtdome_5 = scsm5(:,29);

% ship winds
imrspd_5 = scsm5(:,23);
imum_5 = scsm5(:,25);
imum_port_5 = scsm5(:,44);
imum_stbd_5 = scsm5(:,46);
% recompute wdir
imdm_5 = interval_avg(scsm(:,1), unwrap(scsm(:,24)*rdcon), jd_5bin);
imdm_5 = imdm_5(:,2);
imdm_5 = mod(imdm_5/rdcon + 360*6, 360);
imdm_port_5 = interval_avg(scsm(:,1), unwrap(scsm(:,43)*rdcon), jd_5bin);
imdm_port_5 = imdm_port_5(:,2);
imdm_port_5 = mod(imdm_port_5/rdcon + 360*6, 360);
imdm_stbd_5 = interval_avg(scsm(:,1), unwrap(scsm(:,45)*rdcon), jd_5bin);
imdm_stbd_5 = imdm_stbd_5(:,2);
imdm_stbd_5 = mod(imdm_stbd_5/rdcon + 360*6, 360);
% recompute rwdir
imrdir_5 = interval_avg(scsm(:,1), unwrap(scsm(:,22)*rdcon), jd_5bin);
imrdir_5 = imrdir_5(:,2);
imrdir_5 = mod(imrdir_5/rdcon + 720, 360);
imrdir_5(imrdir_5>180) = imrdir_5(imrdir_5>180) - 360; % +/- 180 deg format

% 5-min PSD met & radiation
metm5 = interval_avg(metm(:,1), metm(:,2:18), jd_5bin);
Tvais_5 = metm5(:,2);
Rhvais_5 = metm5(:,3);
Tsea_5 = metm5(:,4);
Tsea_5(jd_5min>301.229166666) = NaN; % snake dead after this date/time
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)*0;  % rain gauge not working
press_5 = metm5(:,14);
aspir_on_5 = metm5(:,15);
org_carrier_5 = metm5(:,16);
rh337_5 = metm5(:,17);
t337_5 = metm5(:,18);
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:7), 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:5), 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; % replace with updated 1-min hed

% 5-min Licor
licm5 = interval_avg(licm(:,1), licm(:,2:10), 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 fast pressure
fps5 = interval_avg(fpsm(:,1), fpsm(:,2:3), jd_5bin);
fps0_5 = fps5(:,2);
fps1_5 = fps5(:,3);

% 5-min Sonic Winds, tower (high) sonic
sonm5 = interval_avg(sonm(:,1), sonm(:,2:9), jd_5bin);
U_hi_5 = sonm5(:,2)/0.95; % flow distortion corrected rel wind components
V_hi_5 = sonm5(:,3)/1.15;
W_hi_5 = sonm5(:,4);      % mean W
Ts_hi_5 = sonm5(:,5);     % mean Tsonic
rwdir_hi_5 = unwrap(atan2(V_hi_5,-U_hi_5))/rdcon;     % rel wind dir
rwdir_hi_5 = mod(rwdir_hi_5 + 360*6, 360);            % +/- 360 deg format
rwdir_hi_5(rwdir_hi_5>180) = rwdir_hi_5(rwdir_hi_5>180)-360;% +/- 180 deg format
rwspd_hi_5 = sqrt(U_hi_5.^2 + V_hi_5.^2 + W_hi_5.^2);
Un_hi_5 = U_hi_5.*lrg_5n + V_hi_5.*lrg_5e + SOG_5n;    % true wind components...
Ue_hi_5 = U_hi_5.*lrg_5e - V_hi_5.*lrg_5n + SOG_5e;    % corrected for ship speed
wdir_hi_5 = mod(atan2(-Ue_hi_5,-Un_hi_5)/rdcon + 360, 360);  % true wind dir
wspd_hi_5 = sqrt(Un_hi_5.^2 + Ue_hi_5.^2);                   % true wind speed

% 5-min Sonic Winds, bow (low) sonic
U_lo_5 = sonm5(:,6)/0.95; % flow distortion corrected rel wind components
V_lo_5 = sonm5(:,7)/1.15; % not sure this correction is valid for lower sonic
W_lo_5 = sonm5(:,8);      % mean W
Ts_lo_5 = sonm5(:,9);     % mean Tsonic
rwdir_lo_5 = unwrap(atan2(V_lo_5,-U_lo_5))/rdcon;     % rel wind dir
rwdir_lo_5 = mod(rwdir_lo_5 + 360*6, 360);            % +/- 360 deg format
rwdir_lo_5(rwdir_lo_5>180) = rwdir_lo_5(rwdir_lo_5>180)-360;% +/- 180 deg format
rwspd_lo_5 = sqrt(U_lo_5.^2 + V_lo_5.^2 + W_lo_5.^2);
Un_lo_5 = U_lo_5.*lrg_5n + V_lo_5.*lrg_5e + SOG_5n;    % true wind components...
Ue_lo_5 = U_lo_5.*lrg_5e - V_lo_5.*lrg_5n + SOG_5e;    % corrected for ship speed
wdir_lo_5 = mod(atan2(-Ue_lo_5,-Un_lo_5)/rdcon + 360, 360);  % true wind dir
wspd_lo_5 = sqrt(Un_lo_5.^2 + Ue_lo_5.^2);                   % true wind speed

% recompute imir from raw data, using Tc/Tc relationship derived from NOAA PIRs
imir_5 = (scsm5(:,26)+4.9)/1000/838.6/2.83e-6;  % using 2.83e-6 V per W/m2 calibration
psp_ave = (psp1_5+psp2_5)/2;        % avg NOAA psp value
imTc = imta_5 - 0.10 + 0.00914*psp_ave - 9.13e-6*psp_ave.^2;  % empirical fit
imTd = imta_5 - 0.11 + 0.00966*psp_ave - 1.40e-5*psp_ave.^2;  % empirical fit
imir_5 = imir_5 + 5.67e-8*(imTc+273.15).^4 - 4*5.67e-8*((imTd+273.15).^4-(imTc+273.15).^4);

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

% 30-min SCS
scsm30 = interval_avg(scsm(:,1), scsm(:,2:50), jd_30bin);     % average everything
scsm30_std = interval_std(scsm(:,1), scsm(:,2:50), jd_30bin); % std dev of everything
jd_30min = scsm30(:,1);
latm_30 = scsm30(:,2);
lonm_30 = scsm30(:,3);
lrg_30n = scsm30(:,47);
lrg_30e = scsm30(:,48);
lrgm_30 = atan2(lrg_30e, lrg_30n)/rdcon; % from avg components
lrgm_30 = mod(lrgm_30 + 360, 360);
SOG_30n = scsm30(:,49);
SOG_30e = scsm30(:,50);
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(:,4);     % std dev in ship speed
imrh_30 = scsm30(:,16);
imta_30 = scsm30(:,15);
tsgm_30 = scsm30(:,17);
tssm_30 = scsm30(:,19);
ct15_30 = scsm30(:,32) + 0.9; % corrected to sea snake
impress_30 = scsm30(:,14);
imqa_30 = qair3([impress_30, imta_30, imrh_30]);
imsol_30 = scsm30(:,30);
% imir_30 = scsm30(:,27);  % recompute below

% ship winds
imum_30 = scsm30(:,25);

% recompute wdir
imdm_30 = interval_avg(scsm(:,1), unwrap(scsm(:,23)*rdcon), jd_30bin);
imdm_30 = imdm_30(:,2);
imdm_30 = mod(imdm_30/rdcon + 360*6, 360);

% 30-min PSD gps
gprm30 = interval_avg(gprm(:,1), gprm(:,2:7), jd_30bin);

% 30-min PSD met
metm30 = interval_avg(metm(:,1), metm(:,2:18), jd_30bin);
Tvais_30 = metm30(:,2);
Rhvais_30 = metm30(:,3);
press_30 = metm30(:,14);
Tsea_30 = metm30(:,4);
psp1_30 = metm30(:,5);
pir1_30 = metm30(:,6);
psp2_30 = metm30(:,7);
pir2_30 = metm30(:,8);
org_30 = metm30(:,13)*0;  % rain gauge not working
rh337_30 = metm30(:,17);
t337_30 = metm30(:,18);
qair_30 = qair3([press_30,Tvais_30,Rhvais_30]);
qsea_30 = qsea(Tsea_30);

% 30-min Sonic Winds for tower (high) sonic
sonm30 = interval_avg(sonm(:,1), sonm(:,2:9), jd_30bin);
U_hi_30 = sonm30(:,2)/0.95; % flow distortion corrected rel wind components
V_hi_30 = sonm30(:,3)/1.15;
W_hi_30 = sonm30(:,4);
rwdir_hi_30 = unwrap(atan2(V_hi_30,-U_hi_30))/rdcon;         % rel wind dir
rwdir_hi_30 = mod(rwdir_hi_30 + 360*6, 360);                 % +/- 360 deg format
rwdir_hi_30(rwdir_hi_30>180) = rwdir_hi_30(rwdir_hi_30>180)-360;   % +/- 180 deg format
rwspd_hi_30 = sqrt(U_hi_30.^2 + V_hi_30.^2 + W_hi_30.^2);
Un_hi_30 = U_hi_30.*lrg_30n + V_hi_30.*lrg_30e + SOG_30n;    % true wind components...
Ue_hi_30 = U_hi_30.*lrg_30e - V_hi_30.*lrg_30n + SOG_30e;    % corrected for ship speed
wdir_30 = mod(atan2(-Ue_hi_30,-Un_hi_30)/rdcon + 360, 360);  % true wind dir
wspd_30 = sqrt(Un_hi_30.^2 + Ue_hi_30.^2);                   % true wind speed

% 30-min Sonic Winds for bow (lower) sonic
U_lo_30 = sonm30(:,6)/0.95; % flow distortion corrected rel wind components
V_lo_30 = sonm30(:,7)/1.15; % not sure this correction is valid for lower sonic
W_lo_30 = sonm30(:,8);
rwdir_lo_30 = unwrap(atan2(V_lo_30,-U_lo_30))/rdcon;         % rel wind dir
rwdir_lo_30 = mod(rwdir_lo_30 + 360*6, 360);                 % +/- 360 deg format
rwdir_lo_30(rwdir_lo_30>180) = rwdir_lo_30(rwdir_lo_30>180)-360;   % +/- 180 deg format
rwspd_lo_30 = sqrt(U_lo_30.^2 + V_lo_30.^2 + W_lo_30.^2);
Un_lo_30 = U_lo_30.*lrg_30n + V_lo_30.*lrg_30e + SOG_30n;    % true wind components...
Ue_lo_30 = U_lo_30.*lrg_30e - V_lo_30.*lrg_30n + SOG_30e;    % corrected for ship speed
wdir_30 = mod(atan2(-Ue_lo_30,-Un_lo_30)/rdcon + 360, 360);  % true wind dir
wspd_30 = sqrt(Un_lo_30.^2 + Ue_lo_30.^2);                   % true wind speed

% recompute imir from raw data, using Tc/Tc relationship derived from NOAA PIRs
imir_30 = (scsm30(:,26)+4.9)/1000/838.6/2.83e-6;  % using 2.83e-6 V per W/m2 calibration
psp_ave = (psp1_30+psp2_30)/2;        % avg NOAA psp value
imTc = imta_30 - 0.10 + 0.00914*psp_ave - 9.13e-6*psp_ave.^2;  % empirical fit
imTd = imta_30 - 0.11 + 0.00966*psp_ave - 1.40e-5*psp_ave.^2;  % empirical fit
imir_30 = imir_30 + 5.67e-8*(imTc+273.15).^4 - 4*5.67e-8*((imTd+273.15).^4-(imTc+273.15).^4);

%% Some plots...

if plotit
    % PRESSURE
    figure; plot(jd_5min,press_5+zp_psd*0.125,'r',jd_5min,fps0_5+zp_tower*0.125,'b',jd_5min,fps1_5+zp_bow*0.125,'g',jd_5min,impress_5+zp_ship*0.125,'k');
    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('Atmospheric Pressure (mb), adjusted to sea level'); grid;
    axis([ddd ddd+1 min(press_5(~isnan(press_5)))+zp_psd*0.125-2  max(press_5(~isnan(press_5)))+zp_psd*0.125+2]);
    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','Paro-tower','Paro-bow','Ship','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;

    % AIR TEMPERATURES
    figure; plot(jd_5min,Tc1_5,'r',jd_5min,Td1_5,'y',jd_5min,Tc2_5,'m',jd_5min,Td2_5,'k',jd_5min,Tsea_5,'b',jd_5min,Tvais_5,'g.',jd_5min,imta_5,'k.',jd_5min,t337_5,'r.');
    % hold on; plot(jd_5min,imtcase_5,'b.',jd_5min,imtdome_5,'m.');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Temperatures',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    legend('Tc1','Td1','Tc2','Td2','Seasnake','Tvais','Ship','T337','Tc_ship','Td_ship','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,'_PSD_Temperatures_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % SOLAR RADIATON
    figure; plot(jd_5min, psp1_5, 'r', jd_5min, psp2_5,'b',jd_5min,imsol_5,'g');
    xlabel('Hour (UTC)'); ylabel('Solar Flux (W/m^2)'); legend('PSP1','PSP2','SCS','Location','BestOutside');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Precision Spectral Pyranometers',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    axis([ddd ddd+1 -50 150]); grid; legend('PSP1','PSP2','Ship','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,'_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, pir1_5, 'r', jd_5min, pir2_5,'b',jd_5min,imir_5,'g');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Precision 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('PIR1','PIR2','Ship','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,'_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, rh337_5,'b',jd_5min, imrh_5,'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 115]); grid; legend('RHvias','RH337','Ship','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_5min, org_5, 'r');
    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;
    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_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_backflow_indicator_',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_rain_gauge_function_',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, sst38_5,'g',jd_5min, ct15_5, 'k');
    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)'); axis([ddd,ddd+1,min(ct15_5)-1,max(tsgm_5)+1]);
    legend('SeaSnake','SBE45','SBE38','CT15','location','BestOutside');
    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;

    % LAT/LON TRACK PLOT
    figure; plot(lon_5, lat_5, 'b', lonm_5, latm_5, 'r');
    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)'); grid;
    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_5min, sog_5, jd_5min, sogm_5);
    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_5min, cog_5,'b.', jd_5min, cogm_5,'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;

    % MAP PLOT - m_map folder should be in matlab path
    figure;
    m_proj('Albers Equal-Area Conic','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('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;

    % PITCH ANGLE
    pitch = hedm(:,3); % this is 1Hz
    ii = find(isnan(pitch)); pitch(ii) = nanmean1(pitch);
    figure; plot(hedm(:,1),pitch,'b-'); grid;
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Pitch Angle ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    xlabel('Hour (UTC)'); ylabel(['Pitch (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_Pitch_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % SONIC TEMPERATURE
    figure; plot(jd_5min, Ts_hi_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 RELATIVE WIND SPEED COMPONENTS
    figure; plot(jd_5min, U_hi_5, jd_5min, V_hi_5, jd_5min, W_hi_5, jd_5min, rwspd_hi_5);
    xlabel('Hour (UTC)'); ylabel('Wind speed (m/s)'); xlim([ddd ddd+1]);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Tower Sonic Rel Wind Components',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    legend('U','V','W','speed','Location','BestOutside'); 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,'_Tower_Sonic_Winds_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;
    
    figure; plot(jd_5min, U_lo_5, jd_5min, V_lo_5, jd_5min, W_lo_5, jd_5min, rwspd_lo_5);
    xlabel('Hour (UTC)'); ylabel('Wind speed (m/s)'); xlim([ddd ddd+1]);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Bow Sonic Rel Wind Components',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    legend('U','V','W','speed','Location','BestOutside'); 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,'_Bow_Sonic_Winds_',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_5min, rwdir_hi_5,'r.',jd_5min,rwdir_lo_5,'b.', jd_5min, imrdir_5, 'k.',...	% data points
        [ddd ddd+1],[45 45],'g', ...                        % 45 degree sector
        [ddd ddd+1],[-45 -45],'g', ...
        [ddd ddd+1],[60 60],'b', ...                        % 60 degree sector
        [ddd ddd+1],[-60 -60],'b');
    xlabel('Hour (UTC)'); ylabel('Relative Wind Direction (deg)'); axis([ddd ddd+1 -180 180]); grid;
    legend('Tower','Bow','Ship','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_5min, rwspd_hi_5, 'r',jd_5min,rwspd_lo_5, 'b', jd_5min, imrspd_5, 'k');
    legend('Tower','Bow','Ship','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_5min, wdir_hi_5,'r.',jd_5min, wdir_lo_5,'b.', jd_5min, imdm_5, 'k.');
    hold on; plot(jd_5min, imdm_port_5, 'g.',jd_5min, imdm_stbd_5, 'c.');
    xlabel('Hour (UTC)'); ylabel('Wind Direction (deg)'); grid;
    legend('Tower','Bow','Ship','Ship-port','Ship-stbd','location','BestOutside'); axis([ddd ddd+1 0 360]);
    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));
    subplot(2,1,2);
    plot(jd_5min, wspd_hi_5, 'r',jd_5min, wspd_lo_5,'b',jd_5min, imum_5, 'k');
    hold on; plot(jd_5min, imum_port_5, 'g.',jd_5min, imum_stbd_5, 'c.');
    legend('Tower','Bow','Ship','Ship-port','Ship-stbd','location','BestOutside'); 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,'_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('7500','Vaisala aspirated','Location','BestOutside'); 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 min(co2(~isnan(co2)))-10 max(co2(~isnan(co2)))+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_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;
end;

close all;

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

% 1-min SCS met/nav data, traditional format
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 = [scsm(:,1:3)';scsm(:,33)';scsm(:,5)';scsm(:,4)';scsm(:,16)';scsm(:,15)';scsm(:,24)';scsm(:,23)';scsm(:,14)';scsm(:,22)';scsm(:,21)';scsm(:,17)';scsm(:,19)'];
fprintf(fout,'%11.5f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f \r\n',jazz);
%               jd    lat   lon   lrg   cog   sog    rh    ta   wspd  wdir press rwspd  rwdir  sst  sal
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(14,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
jazz(13,:) = metm(:,18);	% t337
jazz(14,:) = metm(:,17);	% rh337
fprintf(fout,'%11.5f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f  %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f \r\n',jazz);
%                jd    pir1  psp1   tc1   td1   ts    ta     rh   org  car   aspr  Pmb  t337  rh337
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,:) = 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
fprintf(fout,'%11.5f %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f \r\n',jazz);
%               jd    pir1  pir2  psp1  psp2  tc1   td1   tc2   td2   ta    rh    tsea
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 (pitch), deg
fprintf(fout,'%11.5f %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f \r\n',jazz);
%                jd    lat   lon   sog   cog   hed   pit
fclose('all');

% 1-min sonic wind results
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(6,1440)*NaN;
jazz(1,:) = jd_min;
jazz(2,:) = sonm1(:,2);      % U high
jazz(3,:) = sonm1(:,3);      % V high
jazz(4,:) = sonm1(:,4);      % W high
jazz(5,:) = sonm1(:,5);      % Tsonic high
jazz(6,:) = rwdir_hi_1;      % rel wind dir high
jazz(7,:) = sonm1(:,6);      % U low
jazz(8,:) = sonm1(:,7);      % V low
jazz(9,:) = sonm1(:,8);      % W low
jazz(10,:) = sonm1(:,9);     % Tsonic low
jazz(11,:) = rwdir_lo_1;     % rel wind dir low
fprintf(fout,'%11.5f %9.2f %9.2f %9.2f %9.2f %9.2f   %9.2f %9.2f %9.2f %9.2f %9.2f \r\n',jazz);
%                jd   U_hi  V_hi  W_hi  T_hi rdir_hi  U_lo  V_lo  W_lo T_lo rdir_lo
fclose('all');

% 1-min licor results
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.5f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f \r\n',jazz);
%                jd   H2O   CO2   CO2a   T     P    AGC    mr   rhoa   q
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

% using PSD data
y = coare35vn(wspd_hi_5,zu_hi,Tvais_5,zt,Rhvais_5,zq,press_5,ct15_5,psp1_5,pir2_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);

Rns = psp1_5*0.7;  % 0.7 est for ice, 0.92 for water
Rnl = 0.985*(pir2_5 - 5.67e-8*(ct15_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','BestOutside'); grid;
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD 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','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,'_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).  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,'_PSD_Friction_Velocity_5min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    close all;
end;

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

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 = zeros(32,288)*NaN;
jazz(1,:) = jd_5min;            % Decimal julian day at beginning of line
jazz(2,:) = wspd_hi_5;          % psd true wind speed at 16.5 m ht, m/s
jazz(3,:) = wdir_hi_5;          % psd true wind direction, deg
jazz(4,:) = Tsea_5;             % psd seasnake T, C
jazz(5,:) = tsgm_5;             % ship theromsalinograph T, C  bow
jazz(6,:) = tssm_5;             % ship theromsalinograph salinity, psu  bow
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,:) = psp1_5;            % psd solar flux, w/m^2
jazz(11,:) = pir2_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_hi_5;           % rel wind speed, m/s
jazz(16,:) = rwdir_hi_5;           % rel wind direction, deg
jazz(17,:) = latm_5;            % decimal latitude, deg
jazz(18,:) = lonm_5;            % decimal longtude, deg
jazz(19,:) = ones(size(jd_5min))*dtsg;	% 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_5;            % IMET air temp, C (21 m)
jazz(26,:) = imqa_5;            % IMET air specific humidity, g/kg (21 m)
jazz(27,:) = imum_5;            % IMET true wind speed, m/s (21 m)
jazz(28,:) = imdm_5;            % IMET true wind direction, deg (21 m)
jazz(29,:) = imsol_5;           % IMET solar flux, w/m^2
jazz(30,:) = imir_5;            % IMET IR flux, w/m^2
jazz(31,:) = impress_5;         % BP, mb
jazz(32,:) = imrh_5;            % RH (%)
%               1        2       3       4       5       6       7       8       9      10      11      12      13      14      15      16      17       18       19      20      21      22      23      24       25      26      27      28      29      30      31       32
fprintf(fout,'%11.5f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %12.4f\t %12.4f\t %9.2f\t %9.2f\t %8.3f\t %9.2f\t %9.2f\t %9.2f\t  %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t  %9.2f\t \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

% using PSD data

y = coare35vn(wspd_30,zu_hi,Tvais_30,zt,Rhvais_30,zq,press_30,ct15_30,psp1_30,pir2_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 = psp1_30*0.7;  % 0.7 est for ice, 0.92 for water
Rnl = 0.985*(pir2_30 - 5.67e-8*(ct15_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;       %etl true wind speed at 16.5 m ht, m/s
jazz(3,:) = wdir_30;       %etl true wind direction, deg
jazz(4,:) = Tsea_30;       %etl seasnake T, C
jazz(5,:) = tsgm_30;       %ship theromsalinograph T, C  bow
jazz(6,:) = tssm_30;       %ship theromsalinograph salinity, psu  bow
jazz(7,:) = Tvais_30;      %etl air T at 18 m, C
jazz(8,:) = qsea_30;       %etl air specific humidity at sea surface, g/kg
jazz(9,:) = qair_30;       %etl air specific humidity at 18 m, g/kg
jazz(10,:) = psp1_30;      %etl solar flux, w/m^2
jazz(11,:) = pir2_30;      %etl IR flux, w/m^2
jazz(12,:) = org_30;       %etl 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_hi_30;  %rel wind speed, m/s
jazz(16,:) = rwdir_hi_30;  %rel wind direction, deg
jazz(17,:) = latm_30;      %decimal latitude, deg
jazz(18,:) = lonm_30;      %decimal longtude, deg
jazz(19,:) = ones(size(jd_30min))*zt;       %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_30;      %IMET air temp, C (21 m)
jazz(26,:) = imqa_30;      %IMET air specific humidity, g/kg (21 m)
jazz(27,:) = imum_30;      %IMET true wind speed, m/s (21 m)
jazz(28,:) = imdm_30;      %IMET true wind direction, deg (21 m)
jazz(29,:) = imsol_30;     %IMET solar flux, w/m^2
jazz(30,:) = imir_30;      %IMET IR flux, w/m^2
jazz(31,:) = impress_30;   %BP, mb
jazz(32,:) = Rhvais_30;    %RH (%)

%vectorized print
%                 1        2       3      4       5       6       7       8       9      10      11       12     13       14      15     16        17     18       19      20      21       22      23      24       25      26      27      28      29      30      31      32
fprintf(fout,'%11.5f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %12.4f\t %12.4f\t %9.2f\t %9.2f\t %8.3f\t %9.2f\t %9.2f\t %9.2f\t  %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t %9.2f\t \r\n',jazz);
fclose all;
clear jazz;

%% radiation models @ 5 min resolution
% Downward 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);

% Solar model
latd = nanmean1(latm_5);
lond = nanmean1(lonm_5);
qa = qair_5; % humidity needed for solar model
qrat = 4.0;
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;
[sw,sz,saz,dirs,sky] = SolarRadiancex(latd,lond,jdx,tutc,watvap,nanmean1(press_5),k1,k2,oz);
Rscl = sw;

%% more plots

if plotit
    if ~isempty(rlclr_5(isfinite(rlclr_5)))
        ii = find(isfinite(rlclr_5));
        figure; plot(jd_5min, pir1_5,'b',jd_5min, pir2_5,'r',jd_5min, imir_5,'g',jd_5min(ii),rlclr_5(ii),'k--'); 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('NOAA-1','NOAA-2','Ship','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,'_PSD-SCS_IRflux_Comparison_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end;
    end;

    if ~isempty(Rscl(isfinite(Rscl)))
        ii = find(isfinite(Rscl));
        figure; plot(jd_5min, psp1_5,'b', jd_5min, psp2_5,'r',jd_5min, imsol_5,'g',jd_5min(ii), Rscl(ii),'k--'); 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('NOAA-1','NOAA-2','Ship','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,'_PSD-SCS_Solarflux_Comparison_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end;
    end;
    close all;
end;

%% backup raw and processed data to external Lacie hard drive
if laCie_backup
    %Check if some folders have not been backed-up
    dirraw = dir(path_raw_data);
    dirLacie = dir(path_backup_raw);
    if ~isempty(dirLacie) && ~isempty(dirraw)
        % find files which are not backed up
        ii = find(~ismember(char(dirraw(:).name),char(dirLacie(:).name),'rows'));
        for ll=1:length(ii)
            copyfile(fullfile(path_raw_data, dirraw(ii(ll)).name), fullfile(path_backup_raw, dirraw(ii(ll)).name));
        end;
    end;
    % copy over processed data and plots
    copyfile(fullfile(path_proc_data,['*',num2str(ddd),'*.txt']), path_backup_proc);
    copyfile(fullfile(path_raw_images,['*',num2str(ddd),'*.png']),path_backup_rawIm);
end;

%% ftp selected data to our server
if data_ftp
    ftp_to_Boulder_v2015;
end

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


