%{
run_motcorr_flux_DYNAMO_2011_v3.m

Driver script for the DYNAMO 2011 RHB flux calculations
Fluxes are processed in 10-min windows

***************************************************************************
08/18/2015 BWB

TSONIC & Q_licor DECORRELATED W/ RESPECT TO 3-AXIS MOTION.
***************************************************************************

This file and the associated read data functions generally require
project-specific editing and alterations.

The assiciated sub-folder of flux utility functions should be in
the Matlab path.  These generally don't require project-specific mods.

Set run parameters in the first block below as desired and execute.

V2 NOTES:
- Version: May 2015 with updated and vectorized code, BWB.
- Motion corrections are Edson98 only - Hare code removed.
- CO2 flux FROM LI 7500 is not computed.
%}

%% Initialize run parameters
close all;
fclose all;
clear;

%% Set cruise parameters
% cruise ID and paths
cruise = 'DYNAMO_2011';      % Cruise name
ship = 'Revelle';            % Research vessel
plotit = true;               % create plots and save to .png format
yyyy = '2011';

% datasetv = 'v1';% legacy method, not yet supported in this version.=
% datasetv = 'v2';  % v2 is tilt correction prior to ship speed correction
datasetv = 'v3';    % v2 + Tsonic & q_licor decorrelated w/ respect to
                    %  3-axis motion

decorr = true;      % set true decorrelate winds & motion...
if decorr
    dcname = 'decorr';
else
    dcname = 'std';
end;

% cruise data root directory - edit as required for your file system and OS
% note: windows uses the '\' separator and mac-unix uses '/'
data_drive = '/Users/byronb/Raw_Data_Files.noindex/';

% build paths relative to data_drive
% these and subsequent path constructs are OS-independent due to use of fullfile()
% these first 3 paths should be preexisting
path_raw_data_flux = fullfile(data_drive,cruise,ship,'flux','Raw','leg2');
path_proc_data_flux = fullfile(data_drive,cruise,ship,'flux','Processed_new','leg2');
path_proc_images_flux = fullfile(data_drive,cruise,ship,'flux','Processed_Images_new','leg2');
% subdirectory for saving daily met/flux plots
path_dailyPlots = fullfile(path_proc_images_flux,['Daily_',dcname,'_',datasetv,'_c25']);
mkdir(path_dailyPlots);
% subdirectory for saving corrected wind timeseries plots
path_windPlots = fullfile(path_proc_images_flux,['WindPlots_',dcname,'_',datasetv]);
mkdir(path_windPlots);
% subdirectory for saving corrected wind spectra plots
path_spectraPlots = fullfile(path_proc_images_flux,['SpectraPlots_',dcname,'_',datasetv]);
mkdir(path_spectraPlots);
% subdirectory for saving corrected wind data
path_uvwStreamData = fullfile(path_proc_data_flux,['uvwStream_',dcname,'_',datasetv]);
mkdir(path_uvwStreamData);
% subdirectory for saving motion data
path_motionData = fullfile(path_proc_data_flux,['motion_',dcname]);
mkdir(path_motionData);
% subdirectory for saving hourly da files
path_da = fullfile(path_proc_data_flux,['da_',dcname,'_',datasetv,'_c25']);
mkdir(path_da);
% subdirectory for saving hourly sp files
path_sp = fullfile(path_proc_data_flux,['sp_',dcname,'_',datasetv]);
mkdir(path_sp);
% subdirectory for saving hourly cr files
path_cr = fullfile(path_proc_data_flux,['cr_',dcname,'_',datasetv,'_c25']);
mkdir(path_cr);

% initialize plot format
graphformat = '.png';  % select graphics format
graphdevice = '-dpng'; % select graphic device

% adjustments to temperature measurements
tsea_adj = 0;        % sea snake adjustment
ta_scs_adj = 0;      % scs air T
ta_psd_adj = 0.05;   % psd vaisala air T
rh_psd_adj = 3;      % psd rh adjustment
bp_psd_adj = -0.58;  % psd pressure adjustment
tc1_adj = -0.15;     % psd pir1 and pir 2 temps
td1_adj = 0.05;
tc2_adj = 0;
td2_adj = 0.1;
adj = [tsea_adj,ta_psd_adj,rh_psd_adj,bp_psd_adj,td1_adj,tc1_adj,td2_adj,tc2_adj];

% Sonic configuration
% valid models: 'WindMasterPro', 'R3', 'R3A', 'R2', 'R2A'
sonicmodel = 'R3A';
fsonic = 10;                  % frequency (Hz) for all fast data
rotationsonic = true;        % true = rotate 30 deg for some Gill sonics
outputsonic = 1;    % defines Ts output: 0 = speed of sound, 1 = degC
Npts = fsonic*600;  % data points in 10 min window @ frequency fsonic

% clear sky model configuration
k1 = 0.1;     % aerosol optical depth, band 1
k2 = 0.1;     % aerosol optical depth, band 2
oz = 0.2;     % column ozone
iv = 3.5;     % column water vapor (cm), if not calculated from obs.

% cruise track map limits
Lonmin = 68;
Lonmax = 110;
Latmin = -15;
Latmax = 10;
PosLims = [Lonmin, Lonmax, Latmin, Latmax];

% sensor heights
% as measured on 01/13/2015: waterline to base of tower = 7.16 m.
zu = 18.16;     % wind speed measurement height (m): 7.16m + 11m
zt = 16.5;      % air T measurement height (m)
zq = 16.5;      % air q measurement height (m)
zP = 10;        % psd pressure sensor height
zu_ship = 15;   % ship SCS sensor heights
zt_ship = 15;
zq_ship = 15;
zP_ship = 15;

% sonic-MotionPak displacement vector [Rx Ry Rz]
sens_disp = [0,0,0.7]; % sonic is 0.7 m above MotionPak

% start and end julian dates for the main flux loop
stjd = 275;
endjd = 278;

% arrays to keep track of bad sonic data
jday_badPnts = [];
badSonPnts = [];
missingPnts = [];

% empty output arrays for daily results: rows 144 = 24 x 6 per hr
% except for sp, file formats compatible with older flux code
% no Hare correction spectra in this version
da = NaN(144,164);
sp = NaN(144,481);
cr = NaN(144,17);

%**************************************************************************
% load mean met file for leg3
% MMF = load('/Users/byronb/Raw_Data_Files.noindex/DYNAMO_2011/Revelle/flux/Processed_old/Revelle10minutesLeg3_r2.mat');
% ocean current components
% CN = -MMF.cspd.*cos((MMF.cdir)./180.*pi);  %'toward' North
% CE = -MMF.cspd.*sin((MMF.cdir)./180.*pi);  %'toward' East
%**************************************************************************

%% ----------Main loop----------
for ddd=stjd:endjd  % iterate over days
    [m,d] = yd2md(str2num(yyyy), ddd);
    Vdate = [str2num(yyyy) m d];
    jd = sprintf('%03i',ddd);

    for hhh=00:23;   % iterate over hours
        %% read raw data files...
        hr = sprintf('%02i',hhh);

        % read scs and interp to 1 Hz
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['scs0',yyyy(3:4),jd,hr,'_raw.txt']);
        if exist(dfl,'file')==0
            continue  % if none continue to next hour
        end;
        scs = read_scs_DYNAMO_2011(dfl,ddd,hhh,ta_scs_adj,PosLims);

        % read psd slow met files
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['met0',yyyy(3:4),jd,hr,'_raw.txt']);
        if exist(dfl,'file')==0
            continue  % if none continue to next hour
        end;
        met = read_met_DYNAMO_2011(dfl,ddd,hhh,adj);

        % read mot and interp to 10 Hz
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['mot0',yyyy(3:4),jd,hr,'_raw.txt']);
        if exist(dfl,'file')==0
            continue  % if none continue to next hour
        end;
        mot = read_motion_DYNAMO_2011(dfl,ddd,hhh);

        % read son and interp to 10 Hz
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['son0',yyyy(3:4),jd,hr,'_raw.txt']);
        if exist(dfl,'file')==0
            continue  % if none continue to next hour
        end;
        [son,badPnts,missing] = read_sonic_DYNAMO_2011(dfl,ddd,hhh,sonicmodel,rotationsonic);
        badPnts = badPnts(:); % convert to col vectors
        missing = missing(:);

        % read gps and interp to 10 Hz
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['gprm',yyyy(3:4),jd,hr,'_raw.txt']);
        if exist(dfl,'file')==0;
            continue; % if none continue to next hour
        end;
        gps = read_gprm_DYNAMO_2011(dfl,ddd,hhh,PosLims);

        % read hed and interp to 10 Hz
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['hed0',yyyy(3:4),jd,hr,'_raw.txt']);
        if exist(dfl,'file')==0;
            continue; % if none continue to next hour
        end;
        hed = read_hed_pitch_DYNAMO_2011(dfl,ddd,hhh);

        % read licor
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['lic0',yyyy(3:4),jd,hr,'_raw.txt']);
        if exist(dfl,'file')==0;
            lic = NaN(Npts*6,7);   % return empty array if no file
            lic(:,1) = son(:,1);
        else
            lic = read_licor_DYNAMO_2011(dfl,ddd,hhh);
        end;

        %% copy out selected variables
        jd_1min = met(:,1);   % ref 1 min timestamp
        jd_1Hz = scs(:,1);    % ref 1 Hz timestamp
        jd_10Hz = son(:,1);   % ref 10 Hz timestamp

        % 1-Hz scs met variables
        lat = scs(:,2);
        lon = scs(:,3);
        cogm = scs(:,4);
        sogm = scs(:,5);
        tsgm = scs(:,6);
        tssm = scs(:,7);
        imrh = scs(:,8);
        imta = scs(:,9);
        imsol = scs(:,10);
        orgm = scs(:,11);
        lrg = scs(:,12);
        imum = scs(:,13);
        imdm = scs(:,14);
        imir = scs(:,15);
        imPress = scs(:,16);
        imqa = scs(:,17);

        % 1-min psd met variables
        Tvais = met(:,2);
        Rhvais = met(:,3);
        Tsea = met(:,4);
        psp = met(:,5);
        pir = met(:,6);
        Td = met(:,7);
        org = met(:,8);
        press = met(:,9);
        aspir_on = met(:,10);
        org_carrier = met(:,11);

        %% compute 10 Hz heading from 15 sec ship lrg heading
%         rdcon = pi/180;
%         lrgc = cos(lrg.*rdcon);    % convert to N
%         lrgs = sin(lrg.*rdcon);    % convert to E
%         % interpolate heading to 10 Hz for motion correction
%         lrgcm = interp1(jd_1Hz,lrgc,jd_10Hz,'linear','extrap');
%         lrgsm = interp1(jd_1Hz,lrgs,jd_10Hz,'linear','extrap');
%         % be sure NaNs are gone
%         lrgcm = replace_NaN_nearest_neighbor(lrgcm);
%         lrgsm = replace_NaN_nearest_neighbor(lrgsm);
%         % compute heading in radians - add 2pi so all are positive
%         heading = mod(atan2(lrgsm,lrgcm)+2*pi, 2*pi);
%         shead = sin(heading);  % 10Hz components
%         chead = cos(heading);

        % 10Hz heading from psd gps
        heading = replace_NaN_nearest_neighbor(hed(:,2));
        heading = heading*pi/180;   % to radians
        shead = sin(heading);       % 10Hz components
        chead = cos(heading);

        %% E98 motion corrections, with decorr if desired
        % using motcorr2.m - decorrelates Tsonic w/respect to 3-axis motion
        [uvw,T_10Hz,uvwplat,xyzplat,euler,lagPnts] = motcorr2(son,mot,heading,sens_disp,fsonic,decorr);
        lagmat = [ddd hhh lagPnts];
        uvw = uvw'; % transpose to column vectors for this code
        uvwplat = uvwplat';
        xyzplat = xyzplat';
        euler = euler';
        T_10Hz = T_10Hz';   % this is sonic temperature

        %% correct uv for ship velocity
        uvw_str = uvw.*0;   % preallocate output arrays
        uvw_rot = uvw.*0;
        tilt = zeros(1,6);	% for saving tilt and azm this hr
        azm = zeros(1,6);

        % first remove tilt
        % double rotation into streamlines over 10-min segs of 10Hz data
        ii=1;
        for jj=1:6
            [uvw_str(ii:ii+Npts-1,:),azm(jj),tilt(jj)] = dbl_rot(uvw(ii:ii+Npts-1,:));
            ii = ii + Npts;
        end;
        % single-rotate azimuth ccw, back to earth coordinates
        ii = 1;
        for jj=1:6
            uvw_rot(ii:ii+Npts-1,:) = azm_rot(uvw_str(ii:ii+Npts-1,:),azm(jj),1);
            ii = ii + Npts;
        end;

        % then sum with N/E ship speed components to remove platform velocity
        uvw_rot(:,1) = uvw_rot(:,1) + gps(:,4); % Ngps
        uvw_rot(:,2) = uvw_rot(:,2) - gps(:,5); % Egps


        % finally compute new azimuth and rotate uv clockwise into streamlines
        % mean v & w should now be zero
        ii = 1;
        for jj=1:6
            azm2 = atan2(mean(uvw_rot(ii:ii+Npts-1,2)), mean(uvw_rot(ii:ii+Npts-1,1)));
            uvw_str(ii:ii+Npts-1,:) = azm_rot(uvw_rot(ii:ii+Npts-1,:),azm2,0);
            ii = ii + Npts;
        end;
        disp(['mean W = ',num2str(mean(uvw_str(:,3)),'%3.2f'),'  mean V = ',num2str(mean(uvw_str(:,2)),'%3.2f')]);

        % add timestamp and sonicT to uvw_str array
        uvw_str = [jd_10Hz, uvw_str, T_10Hz];
        clear uvw_rot

        %% compute 10-min statistics

        % Set up to process data in 10-min blocks
        jd_start = floor(son(1,1)*24)/24;   % start of hour
        jd_end = jd_start + 60/1440;        % start of next hour
        % 10-min bin edges
        %  of bin edges = 7 for 6 10-min bins
        jd_10min = (jd_start:10/1440:jd_end);

        % 10-min psd met statistics - 6x11 output array:
        % 1: jd, 2:Tvais, 3:Rhvais, 4:Tsea, 5:psp, 6:pir, 7:Td
        % 8:org, 9:press, 10:aspir_on, 11:org_carrier
        met_10min_avg = interval_avg(met(:,1), met(:,2:end),jd_10min);
        met_10min_std = interval_std(met(:,1), met(:,2:end),jd_10min);
        met_10min_min = interval_min(met(:,1), met(:,2:end),jd_10min);
        met_10min_max = interval_max(met(:,1), met(:,2:end),jd_10min);

        % Define a few 10-min PSD met variables for use below
        ta_psd = met_10min_avg(:,2);
        rh_psd = met_10min_avg(:,3);
        ts_psd = met_10min_avg(:,4);
        rs_psd = met_10min_avg(:,5);
        rl_psd = met_10min_avg(:,6);
        rain_psd = met_10min_avg(:,8);
        bp_psd = met_10min_avg(:,9);

        % 10-min stats for scs data; #cols depends on ship set up.
        % for DYNAMO2011 output is 6x17 array as follows:
        % 1:jd, 2:lat, 3:lon, 4:cogm, 5:sogm-m/s, 6:tsgm, 7:tssm, 8:imrh
        % 9:imta, 10:imsol, 11:orgm, 12:lrg, 13:imum, 14:imdm, 15:imir
        % 16:impress, 17:imqa
        scs_10min_avg = interval_avg(scs(:,1), scs(:,2:end),jd_10min);
        scs_10min_std = interval_std(scs(:,1), scs(:,2:end),jd_10min);
        scs_10min_min = interval_min(scs(:,1), scs(:,2:end),jd_10min);
        scs_10min_max = interval_max(scs(:,1), scs(:,2:end),jd_10min);
        % redo cog - unwrap before averaging, then compute mod
        temp = interval_avg(scs(:,1),unwrap(scs(:,4)*pi/180.0),jd_10min)*180.0/pi;
        scs_10min_avg(:,4) = mod(temp(:,2)+360, 360);
        temp = interval_std(scs(:,1),unwrap(scs(:,4)*pi/180.0),jd_10min)*180.0/pi;
        scs_10min_std(:,4) = temp(:,2);
        % redo lrg, same reason
        temp = interval_avg(scs(:,1),unwrap(scs(:,12)*pi/180.0),jd_10min)*180.0/pi;
        scs_10min_avg(:,12) = mod(temp(:,2)+360, 360);
        temp = interval_std(scs(:,1),unwrap(scs(:,12)*pi/180.0),jd_10min)*180.0/pi;
        scs_10min_std(:,12) = temp(:,2);
        % redo wdir, same reason
        temp = interval_avg(scs(:,1),unwrap(scs(:,14)*pi/180.0),jd_10min)*180.0/pi;
        scs_10min_avg(:,14) = mod(temp(:,2)+360, 360);
        temp = interval_std(scs(:,1),unwrap(scs(:,14)*pi/180.0),jd_10min)*180.0/pi;
        scs_10min_std(:,14) = temp(:,2);

        % Define a few 10-min variables for use below
        lat_scs = scs_10min_avg(:,2);
        lon_scs = scs_10min_avg(:,3);
        cog_scs = scs_10min_avg(:,4);
        sog_scs = scs_10min_avg(:,5);
        ts_scs = scs_10min_avg(:,6);
        sal_scs = scs_10min_avg(:,7);
        rh_scs = scs_10min_avg(:,8);
        ta_scs = scs_10min_avg(:,9);
        rs_scs = scs_10min_avg(:,10);
        rain_scs = scs_10min_avg(:,11);
        lrg_scs = scs_10min_avg(:,12);
        wspd_scs = scs_10min_avg(:,13);
        wdir_scs = scs_10min_avg(:,14);
        rl_scs = scs_10min_avg(:,15);
        bp_scs = scs_10min_avg(:,16);

        % std dev for sin/cos of lrg heading
        shead_10min_std = interval_std(scs(:,1), sin(scs(:,12)),jd_10min);
        chead_10min_std = interval_std(scs(:,1), cos(scs(:,12)),jd_10min);

        % 10-min psd gps statistics: 6x5 array as follows:
        % 1:jd, 2:cog, 3:sog, m/s, 4:Ngps, 5: Egps
        gps_10min_avg = interval_avg(gps(:,1), gps(:,2:end),jd_10min);
        gps_10min_std = interval_std(gps(:,1), gps(:,2:end),jd_10min);
        gps_10min_min = interval_min(gps(:,1), gps(:,2:end),jd_10min);
        gps_10min_max = interval_max(gps(:,1), gps(:,2:end),jd_10min);
        % redo cog as above
        temp = interval_avg(gps(:,1),unwrap(gps(:,2)*pi/180.0),jd_10min)*180.0/pi;
        gps_10min_avg(:,2) = mod(temp(:,2)+360, 360);
        temp = interval_std(gps(:,1),unwrap(gps(:,2)*pi/180.0),jd_10min)*180.0/pi;
        gps_10min_std(:,2) = temp(:,2);

        % 10-min motion-corrected wind statistics: 6x5 array as follows:
        % 1:jd, 2:ustr, 3:vstr, 4:wstr, 5:Tsonic
        uvwt_10min_avg = interval_avg(uvw_str(:,1), uvw_str(:,2:end),jd_10min);
        uvwt_10min_std = interval_std(uvw_str(:,1), uvw_str(:,2:end),jd_10min);

        % 10-min raw sonic wind averages: 6x5 array as follows:
        % 1:jd, 2:u, 3:v, 4:w, 5:Tsonic
        son_10min_avg = interval_avg(son(:,1), son(:,2:end),jd_10min);

        % 10-min platform velocity std deviations: 6x4 array as follows:
        % 1:jd, 2:uplat, 3:vplat, 4:wplat
        plat_10min_std = interval_std(mot(:,1), uvwplat(:,:),jd_10min);

        % compute true wind spd & dir from N/E sonic components
        % first, apply flow distortion correction to the raw wind speeds
        U = son(:,2)/0.95;  % -5% for head-on wind
        V = son(:,3)/1.15;  % +15% for side-on wind
        Un_10min = zeros(6,1)*NaN; % empty arrays
        Ue_10min = Un_10min;
        ii = 1;
        for jj=1:6
            Un_10min(jj) = (nanmean1(U(ii:ii+Npts-1).*cos(heading(ii:ii+Npts-1)) + ...
                    V(ii:ii+Npts-1).*sin(heading(ii:ii+Npts-1)) + gps_10min_avg(jj,4)));
            Ue_10min(jj) = (nanmean1(U(ii:ii+Npts-1).*sin(heading(ii:ii+Npts-1)) - ...
                    V(ii:ii+Npts-1).*cos(heading(ii:ii+Npts-1)) + gps_10min_avg(jj,5)));
            wdir_psd = atan2(-Ue_10min, -Un_10min)*180/pi;
            wdir_psd(wdir_psd<0) = wdir_psd(wdir_psd<0)+360;
            wspd_psd = sqrt(Un_10min.^2 + Ue_10min.^2);
            ii = ii + Npts;
        end;

        % 10-min relative wind statistics
        rwdir = unwrap(atan2(son(:,3),-son(:,2))); % unwrapped rwdir in radians
        % average and convert to +/- 180 deg format
        rwdir_psd = interval_avg(son(:,1),rwdir,jd_10min)*180/pi;
        rwdir_psd = rwdir_psd(:,2);
        rwdir_psd = mod(rwdir_psd + 360, 360);
        rwdir_psd(rwdir_psd>180) = rwdir_psd(rwdir_psd>180) - 360;
        rwdir_psd_std = interval_std(son(:,1),rwdir,jd_10min)*180/pi;
        rwdir_psd_std = rwdir_psd_std(:,2);
        % note: this is streamwise wind speed UNCORRECTED for ship velocity
        rwspd_psd = sqrt(son_10min_avg(:,2).^2 + son_10min_avg(:,3).^2 + son_10min_avg(:,4).^2);

        % 10-min licor H2O and CO2 statistics: 6x5 array output
        % stats are for 'all' data, not screened for bad AGC
        % 1 jd_ref, 2 Licor_H2O(mmol/m3), 3 Licor_CO2(mmol/m3), 4 Licor_CO2a(mol/mol), 5 Licor_T(C), 6 Licor_P(kpa), 7 AGC
        lic_10min_avg_all = interval_avg(lic(:,1),lic(:,2:end),jd_10min);
        lic_10min_std_all = interval_std(lic(:,1),lic(:,2:end),jd_10min);
        lic_10min_max_all = interval_max(lic(:,1),lic(:,2:end),jd_10min);

        % 10-min specific humidities (q = g/kg)
        qa_psd = qair_p([ta_psd,rh_psd,bp_psd]);
        qa_scs = qair_p([ta_scs,rh_scs,bp_scs]);
        qsea_psd = qsea_p([ts_psd,bp_psd]);
        qsea_scs = qsea_p([ts_scs,bp_scs]);

        % 10-min mean air density (kg/m3) and H2O number density (mmol/m3) from bulk met
        rhoa_psd = rhoa(ta_psd,bp_psd,rh_psd);
        % number density from vaisala RH, mmol/m3
        H2O_10min_avg = qa_psd.*rhoa_psd/0.018016;

        % compute fast q from licor mmol/m3, using 10-min moist air density
        q_10Hz_raw = NaN(Npts*6,1);
        ii=1;
        for jj=1:6
            q_10Hz_raw(ii:ii+Npts-1) = lic(ii:ii+Npts-1,2)*18/1000/rhoa_psd(jj);
            ii = ii + Npts;
        end;
        % decorrelate licor q w/respect to 3-axis motion
%         if decorr
%             q_10Hz = decorr_q(q_10Hz_raw,uvwplat,mot(:,2:4));
%         else
            q_10Hz = q_10Hz_raw;
%         end;

        % finally, compute 10-min wu, wv, wt and wq covariances
        wu_cov_10min = interval_cov(uvw_str(:,1),uvw_str(:,4),uvw_str(:,2),jd_10min);
        wv_cov_10min = interval_cov(uvw_str(:,1),uvw_str(:,4),uvw_str(:,3),jd_10min);
        wt_cov_10min = interval_cov(uvw_str(:,1),uvw_str(:,4),T_10Hz,jd_10min); % This is sonic T
        wq_cov_10min = interval_cov(uvw_str(:,1),uvw_str(:,4),q_10Hz,jd_10min);

        % redefine jd_10min from bin edges to timestamp
        jd_10min = met_10min_avg(:,1);

        %% ---------INITIAL COARE BULK MODEL RUNS WITH SHIP AND PSD DATA--------
        %{
        input format: A = coare35vn(u,zu,t,zt,rh,zq,P,ts,Rs,Rl,lat,zi,rain,cp,sigH)

        cp, sigH set to NaN
        zi=600 and jcool=1 (default)

        output: 6 x 21 array: conforms to legacy output format
         1      usr     friction velocity that includes gustiness (m/s)
         2      tau     wind stress (N/m^2)
         3      hsb     sensible heat flux into ocean (W/m^2)
         4      hlb     latent heat flux into ocean (W/m^2)
         5      hlwebb  Webb correction for latent heat flux, add this to directly
                        measured eddy covariance latent heat flux from Licor
         6      tsr     temperature scaling parameter (K)
         7      qsr     specific humidity scaling parameter (g/Kg)
         8      zot     thermal roughness length (m)
         9      zoq     moisture roughness length (m)
         10     Cd      wind stress transfer (drag) coefficient at height zu
         11     Ch      sensible heat transfer coefficient (Stanton number) at ht zu
         12     Ce      latent heat transfer coefficient (Dalton number) at ht zq
         13     L       Obukhov length scale (m)
         14     zet     Monin-Obukhov stability parameter zu/L
         15     dter    cool-skin temperature depression (degC)
         16     dqer    cool-skin humidity depression (degC)
         17     tkt     cool-skin thickness (m)
         18     RF      rain heat flux (W/m2)
         19     Cdn_10  neutral value of drag coefficient at 10m
         20     Chn_10  neutral value of Stanton number at 10m
         21     Cen_10  neutral value of Dalton number at 10m
        %}

        % using rain_psd for both runs and warm layer model
%         coare_psd = coare35vnWarm(jd_10min,wspd_psd,zu,ta_psd,zt,rh_psd,zq,bp_psd,...
%                       ts_psd,rs_psd,rl_psd,lat_scs,lon_scs,600,rain_psd,0.05,10);
%         coare_scs = coare35vnWarm(jd_10min,wspd_scs,zu_ship,ta_scs,zt_ship,rh_scs,zq_ship,...
%                       bp_scs,ts_scs,rs_scs,rl_scs,lat_scs,lon_scs,600,rain_psd,5,10);

%         coare_psd = coare35vn(wspd_psd,zu,ta_psd,zt,rh_psd,zq,bp_psd,ts_psd,...
%                       rs_psd,rl_psd,lat_scs,600,rain_psd,NaN,NaN);
%         coare_scs = coare35vn(wspd_scs,zu_ship,ta_scs,zt_ship,rh_scs,zq_ship,...
%                       bp_scs,ts_scs,rs_scs,rl_scs,lat_scs,600,rain_psd,NaN,NaN);

        %  USE OLDER COARE 2.5 CODE...
        % input = [U odec ts ta qs qa rs rl rain zi P zu zt zq];
        % set odec to zero since U is ship speed corrected above...
        % coul_lkb output format adjusted to match c35 above

        for ii=1:6
            coare_psd(ii,:) = coul_lkb([wspd_psd(ii),0,ts_psd(ii),ta_psd(ii),qsea_psd(ii),qa_psd(ii),...
                                  rs_psd(ii),rl_psd(ii),rain_psd(ii),600,bp_psd(ii),zu,zt,zq]);
            coare_scs(ii,:) = coul_lkb([wspd_scs(ii),0,ts_scs(ii),ta_scs(ii),qsea_scs(ii),qa_scs(ii),...
                                  rs_scs(ii),rl_scs(ii),rain_psd(ii),600,bp_scs(ii),zu_ship,zt_ship,zq_ship]);
        end;

        %% compute covariance heat fluxes
        % correct sonic wT for humidity -  use usb, qsb
        hs_cov_10min = wt_cov_10min(:,2) + 0.51*(ta_psd+273.15).*coare_psd(:,1).*coare_psd(:,7)*1e-3;
        % then convert to heat flux units
        % hs = wt * rhoair * Cp_air = K*m/s * kg/m3 * J/kg*K = W/m2
        % Cp_air = 1004.67 J/kg K
        hs_cov_10min = hs_cov_10min.* rhoa_psd * 1004.67;

        % latent heat vap for water and seawater from ship SST and ship salinity
        [Le_w_10min, Le_sw_10min] = Le_water(ts_scs, sal_scs);
        % hl = wq * rhoair * Le_w = g/kg m/s * kg/m3 * J/g = W/m2
        hl_cov_10min = wq_cov_10min(:,2) .* rhoa_psd .* Le_w_10min;
        % Webb corrected hl
        hlc_cov_10min = hl_cov_10min + coare_psd(:,5);

        %% --------------------COMPUTE SPECTRA AND COSPECTRA--------------------
        % initialize arrays for hourly spectral data
        Fx = zeros(1,48);      % empty arrays for 6x48-bin smoothed spectra
        Pu = zeros(6,48)*NaN;
        Pv = Pu; Pw = Pu; Pt = Pu; Pt2 = Pu; Pq = Pu;
        Cuw = Pu; Cvw = Pu; Ctw = Pu; Cqw = Pu;
        Ctw_raw = Pu; Cqw_raw = Pu;

        Fxx = zeros(1,1+Npts/2);      % empty arrays for raw spectra
        Puu = zeros(6,1+Npts/2)*NaN;  % used for idiss calculations...
        Pvv = Puu; Pww = Puu; Ptt = Puu; Ptt2 = Puu; Pqq = Puu;
        Ts_noise = zeros(6,1);

        % compute variance spectra and cospectra in 10-min windows
        ii = 1;
        for jj=1:6
            [Puu(jj,:),Fxx(1,:)] = psd2(detrend(uvw_str(ii:ii+Npts-1,2)),Npts,fsonic,hamming(Npts));
            Puu = real(Puu);
            [Pu(jj,:),Fx(1,:)] = specsmoo(Puu(jj,:),fsonic);
            [Pvv(jj,:),~] = psd2(detrend(uvw_str(ii:ii+Npts-1,3))',Npts,fsonic,hamming(Npts));
            Pvv = real(Pvv);
            [Pv(jj,:),~] = specsmoo(Pvv(jj,:),fsonic);
            [Pww(jj,:),~] = psd2(detrend(uvw_str(ii:ii+Npts-1,4))',Npts,fsonic,hamming(Npts));
            Pww = real(Pww);
            [Pw(jj,:),~] = specsmoo(Pww(jj,:),fsonic);
            [Ptt(jj,:),~] = psd2(detrend(T_10Hz(ii:ii+Npts-1))',Npts,fsonic,hamming(Npts));
            Ptt = real(Ptt);
            [Pt(jj,:),~] = specsmoo(Ptt(jj,:),fsonic); % This is Tsonic...
            [Pqq(jj,:),~] = psd2(detrend(q_10Hz(ii:ii+Npts-1))',Npts,fsonic,hamming(Npts));
            Pqq = real(Pqq);
            [Pq(jj,:),~] = specsmoo(Pqq(jj,:),fsonic);

            % subtract noise from Tsonic spectra (for idiss)
            % assumes sonic spectrum above 2.5 Hz is noise only (last 5 points in smoothed spectrum)
            Ts_noise(jj) = median(Pt(jj,end-5:end));
            Ptt2(jj,:) = Ptt(jj,:) - Ts_noise(jj);
            [Pt2(jj,:),~] = specsmoo(Ptt2(jj,:),fsonic);
            zz = find(Pt2(jj,:)<0); Pt2(jj,zz) = 0; % be sure there are no negative variances

            [temp,~] = csd2(detrend(uvw_str(ii:ii+Npts-1,4)),detrend(uvw_str(ii:ii+Npts-1,2)),Npts,fsonic,hamming(Npts));
            [Cuw(jj,:),~] = specsmoo(real(temp),fsonic);
            [temp,~] = csd2(detrend(uvw_str(ii:ii+Npts-1,4)),detrend(uvw_str(ii:ii+Npts-1,3)),Npts,fsonic,hamming(Npts));
            [Cvw(jj,:),~] = specsmoo(real(temp),fsonic);
            [temp,~] = csd2(detrend(uvw_str(ii:ii+Npts-1,4)),detrend(T_10Hz(ii:ii+Npts-1)),Npts,fsonic,hamming(Npts));
            [Ctw(jj,:),~] = specsmoo(real(temp),fsonic);
            [temp,~] = csd2(detrend(uvw_str(ii:ii+Npts-1,4)),detrend(son(ii:ii+Npts-1,5)),Npts,fsonic,hamming(Npts));
            [Ctw_raw(jj,:),~] = specsmoo(real(temp),fsonic);
            [temp,~] = csd2(detrend(uvw_str(ii:ii+Npts-1,4)),detrend(q_10Hz(ii:ii+Npts-1)),Npts,fsonic,hamming(Npts));
            [Cqw(jj,:),~] = specsmoo(real(temp),fsonic);
            [temp,~] = csd2(detrend(uvw_str(ii:ii+Npts-1,4)),detrend(q_10Hz_raw(ii:ii+Npts-1)),Npts,fsonic,hamming(Npts));
            [Cqw_raw(jj,:),~] = specsmoo(real(temp),fsonic);
            clear temp
            ii = ii + Npts;
        end;

        %% inertial dissipation calculations
        % METHOD A: compute Cx2 from median over inertial subrange in
        % smoothed spectra
        % Cu2 / Cw2: array length same as rwspd_psd
        fminU = find(Fx>0.8,1,'first');
        fminU = fminU - 2;
        if fminU < 1
            fminU = 1;
        end;
        fmaxU = fminU + 2;
        if fmaxU > length(Fx)
            fmaxU = length(Fx);
        end;
        Cu2a = idiss_struc_func(Pu, Fx, fminU, fmaxU, rwspd_psd, 0);
        Cw2a = idiss_struc_func(Pw, Fx, fminU, fmaxU, rwspd_psd, 0);

         % compute as median from raw spectra
%         fminU = find(Fxx>0.7,1,'first');
%         fmaxU = find(Fxx<1.0,1,'last');
%         Cu2a = idiss_struc_func(Puu, Fxx, fminU, fmaxU, rwspd_psd, 0);
%         Cw2a = idiss_struc_func(Pww, Fxx, fminU, fmaxU, rwspd_psd, 0);

        % Ct2 - this is Tsonic...
        fminT = find(Fx>0.35,1,'first');
        fminT = fminT - 2;
        if fminT < 1
            fminT = 1;
        end;
        fmaxT = fminT + 2;
        if fmaxU > length(Fx)
            fmaxT = length(Fx);
        end;
        Ct2a = idiss_struc_func(Pt2, Fx, fminT, fmaxT, rwspd_psd, 0);

%         fminT = find(Fxx>0.28,1,'first');
%         fmaxT = find(Fxx<0.4,1,'last');
%         Ct2a = idiss_struc_func(Ptt2, Fxx, fminT, fmaxT, rwspd_psd, 0);

        % Cq2 - q licor
        fminq = find(Fx>0.7,1,'first');
        fminq = fminq - 2;
        if fminq < 1
            fminq = 1;
        end;
        fmaxq = fminq + 2;
        if fmaxq > length(Fx)
            fmaxq = length(Fx);
        end;
        Cq2a = idiss_struc_func(Pq, Fx, fminq, fmaxq, rwspd_psd, 0);

%         fminq = find(Fxx>0.6,1,'first');
%         fmaxq = find(Fxx<0.8,1,'last');
%         Cq2a = idiss_struc_func(Pqq, Fxx, fminq, fmaxq, rwspd_psd, 0);

        % apply Wyngaard & Clifford correction
        % f1, f2, f3 are same length as uvar, etc.
        ustr2 = rwspd_psd.^2;
        uvar = uvwt_10min_std(:,2).^2;
        vvar = uvwt_10min_std(:,3).^2;
        wvar = uvwt_10min_std(:,4).^2;
        f1 = 1 - (uvar/9.0 - 2*vvar/3.0 - 2*wvar/3.0)./ustr2;
        f2 = 1 - (uvar/9.0 - vvar/12.0 - wvar/3.0)./ustr2;
        f3 = 1 - (uvar/9.0 - vvar/3.0 - wvar/3.0)./ustr2;

        Cu2a = Cu2a./f1;
        Cw2a = Cw2a./f2;
        Ct2a = Ct2a./f3;
        Cq2a = Cq2a./f3;

        % compute usid, tsid, qsid, Lid
        L_bulk = coare_psd(:,13);
        usr_bulk = coare_psd(:,1);
        tsr_bulk = coare_psd(:,6);
        qsr_bulk = coare_psd(:,7); % this is g/kg
        % qsid output is kg/kg
        [usid_a,tsid_a,qsid_a,Lid_a] = idiss(Cu2a,Ct2a,Cq2a,qa_psd,ta_psd,zu,zq,L_bulk,usr_bulk,tsr_bulk,qsr_bulk);

        %% inertial dissipation calculations
        % METHOD B: compute Cx2 from spectral fit over inertial subrange in
        % smoothed spectra
        % Cu2 / Cw2
        fminU = find(Fx>0.6*fsonic/20,1,'first');
        if fminU<1
            fminU = 1;
        end;
        fmaxU = find(Fx<3.5*fsonic/20,1,'last');
        if fmaxU>length(Fx)
            fmaxU = length(Fx);
        end;
        Cu2b = idiss_spec_fit(Pu, Fx, fminU, fmaxU, rwspd_psd);
        Cw2b = idiss_spec_fit(Pw, Fx, fminU, fmaxU, rwspd_psd);

        % Ct2
        fminT = find(Fx>0.35*fsonic/20,1,'first');
        if fminT<1
            fminT = 1;
        end;
        fmaxT = find(Fx<3*fsonic/20,1,'last');
        if fmaxT>length(Fx)
            fmaxT = length(Fx);
        end;
        Ct2b = idiss_spec_fit(Pt2, Fx, fminT, fmaxT, rwspd_psd);

        % Cq2
        fminq = find(Fx>0.7*fsonic/20,1,'first');
        if fminq<1
            fminq = 1;
        end;
        fmaxq = find(Fx<3.6*fsonic/20,1,'last');
        if fmaxq>length(Fx)
            fmaxq = length(Fx);
        end;
        Cq2b = idiss_spec_fit(Pq, Fx, fminq, fmaxq, rwspd_psd);

        % apply Wyngaard & Clifford correction
        Cu2b = Cu2b./f1;
        Cw2b = Cw2b./f2;
        Ct2b = Ct2b./f3;
        Cq2b = Cq2b./f3;

        % compute usid, tsid, qsid, Lid
        % qsid output is kg/kg
        [usid_b,tsid_b,qsid_b,Lid_b] = idiss(Cu2b,Ct2b,Cq2b,qa_psd,ta_psd,zu,zq,L_bulk,usr_bulk,tsr_bulk,qsr_bulk);


        %% save 10_Hz wind and motion data

        rwdir = mod(rwdir*180/pi + 360, 360);
        rwdir(rwdir>180) = rwdir(rwdir>180) - 360;
        uvwmat = [uvw_str,rwdir,q_10Hz]; %[jd_10Hz,u,v,w,T,rwdir,q_10Hz]
        fpath = fullfile(path_uvwStreamData,['UVWstream_',datasetv,'_',jd,'_',hr,'_',dcname,'.txt']);
        save (fpath,'uvwmat','-ascii');

        motmatall=[jd_10Hz,uvwplat,xyzplat,euler];
        fpath = fullfile(path_motionData,['motion_',datasetv,'_',jd,'_',hr,'_',dcname,'.txt']);
        save (fpath,'motmatall','-ascii','-double');

        %% selected wind correction and spectra plots
        if plotit
            if hhh==0 || hhh==6 || hhh==12 || hhh==18
                % wind timeseries plots
                hr=sprintf('%02i',hhh);
                jd=sprintf('%03i',ddd);
                figure(1);clf;plot(jd_10Hz(1:1000),uvw_str(1:1000,2),'g-',jd_10Hz(1:1000),uvwplat(1:1000,1),'r-',jd_10Hz(1:1000),son(1:1000,2),'b-');
                title([jd,' ',hr,' blue = u sonic, green = u corrected, red = u platform']); grid;
                ppath = fullfile(path_windPlots,['u_',datasetv,'_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                figure(2);clf;plot(jd_10Hz(1:1000),uvw_str(1:1000,3),'g-',jd_10Hz(1:1000),uvwplat(1:1000,2),'r-',jd_10Hz(1:1000),son(1:1000,3),'b-');
                title([jd,' ',hr,' blue = v sonic, green = v corrected, red = v platform']); grid;
                ppath = fullfile(path_windPlots,['v_',datasetv,'_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                figure(3);clf;plot(jd_10Hz(1:1000),uvw_str(1:1000,4),'g-',jd_10Hz(1:1000),uvwplat(1:1000,3),'r-',jd_10Hz(1:1000),son(1:1000,4),'b-');
                hold on; plot([jd_10Hz(1) jd_10Hz(1000)],[0 0],'k--')
                title([jd,' ',hr,' blue = w sonic, green = w corrected, red = w platform']); grid;
                ppath = fullfile(path_windPlots,['w_',datasetv,'_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                % cospectra plots
                wu = mean(Cuw,1);
                figure(4); clf; semilogx(Fx,Fx.*wu,'g-',[1e-3 10 ],[0 0]);
                title([jd,' ',hr,' wu Cospectrum']); grid;
                ylabel('C_{uw}(f) x f'); xlabel('f (Hz)');
                ppath = fullfile(path_spectraPlots,['wuCospec_',datasetv,'_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                wq = mean(Cqw,1); wq_raw = mean(Cqw_raw,1);
                figure(5); clf; semilogx(Fx,Fx.*wq,'b',Fx,Fx.*wq_raw,'b--',[1e-3 10 ],[0 0]);
                title([jd,' ',hr,' wq Cospectrum: raw(dashed), decorr(solid)']); grid;
                ylabel('C_{wq}(f) x f'); xlabel('f (Hz)');
                ppath = fullfile(path_spectraPlots,['wqCospec_',datasetv,'_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                wt = mean(Ctw,1); wt_raw = mean(Ctw_raw,1);
                figure(6); clf; semilogx(Fx,Fx.*wt,'r-',Fx,Fx.*wt_raw,'r--',[1e-3 10 ],[0 0]);
                title([jd,' ',hr,' wT Cospectrum: raw(dashed), decorr(solid)']); grid;
                ylabel('C_{wt}(f) x f'); xlabel('f (Hz)');
                ppath = fullfile(path_spectraPlots,['wtCospec_',datasetv,'_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                % spectra plots
                Sw = mean(Pw,1);
                [Sw_raw,~] = psd2(detrend(son(:,4)),Npts*6,10,hamming(Npts*6));
                [Sw_raw,Fraw] = specsmoo(Sw_raw,10);
                figure(7);clf;loglog(Fx,Fx.*Sw,'b-',Fraw,Fraw.*Sw_raw,'b--');
                hold on; loglog(Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*Sw(35).*Fx(35),'r');
                ylabel('S_{ww} x f'); xlabel('f (Hz)'); legend('corrected','raw','location','northwest');
                title([jd,' ',hr,' Spectrum of vertical wind velocity, W']); grid;
                ppath = fullfile(path_spectraPlots,['wSpectrum_',datasetv,'_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                St = mean(Pt,1); St2 = mean(Pt2,1);
                WN = zeros(size(St)) + nanmean1(Ts_noise); WN(1:24) = NaN;
                figure(8);clf;loglog(Fx,Fx.*St,'g--',Fx,Fx.*St2,'g-',Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*St(35).*Fx(35),'r-');
                hold on; loglog(Fx,Fx.*WN,'k--'); ylabel('S_{tt} x f'); xlabel('f (Hz)');
                title([jd,' ',hr,' Spectrum of sonic temperature, Ts']); grid;
                legend('raw','minus noise','-2/3','white noise','location','northwest');
                ppath = fullfile(path_spectraPlots,['tSpectrum_',datasetv,'_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                Sq = mean(Pq,1);
                figure(9);clf;loglog(Fx,Fx.*Sq,'g-',Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*Sq(35).*Fx(35),'r');
                ylabel('S_{qq} x f'); xlabel('f (Hz)');
                title([jd,' ',hr,' Spectrum of specific humidity, q']); grid;
                ppath = fullfile(path_spectraPlots,['qSpectrum_',datasetv,'_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                Su = mean(Pu,1);
                Sv = mean(Pv,1);
                [Su_raw,~] = psd2(detrend(son(:,2)),Npts*6,10,hamming(Npts*6));
                [Su_raw,~] = specsmoo(real(Su_raw),10);
                [Sv_raw,~] = psd2(detrend(son(:,3)),Npts*6,10,hamming(Npts*6));
                [Sv_raw,~] = specsmoo(real(Sv_raw),10);
                figure(10);clf;loglog(Fx,Fx.*Su,'b-',Fraw,Fraw.*Su_raw,'b--',Fx,Fx.*Sv,'g-',Fraw,Fraw.*Sv_raw,'g--');
                hold on; loglog(Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*Su(35).*Fx(35),'r');
                ylabel('S_{uu} x f'); xlabel('f (Hz)');
                legend('U','U raw','V','V raw','-2/3','location','northwest');
                title([jd,' ',hr,' U V Wind Spectra']); grid;
                ppath = fullfile(path_spectraPlots,['uvSpectrum_',datasetv,'_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                close all;
            end;
        end;

        %% format da, sp, cr output and concatenate to daily array
%{
    da 10-min variable designations as follows:  144 x 164 array

    Note:  historical da format is retained, but Hare corrections not computed
    in this version of the flux code & therefore set to NaN in the da output.

    Also note: COARE input format is for vectorized version 3.5.



    COL     VARIABLE        DESCRIPTION
    1       doy             decimal day of year
    2:10    xship           COARE 3.5 input, scs:[Utrue,Tair,RH,P_mb,Tsea,
                            Rs,Rl,lat,rain]
    11:15   no_data         To maintain historical column numbering...
    16:24   xpsd            COARE 3.5 input, PSD, as above
    25:29   no_data         To maintain historical column numbering...
    30      u_ship          ship mean true wind speed (m/s)
    31      u_ship_min      ship minimum true wind speed (m/s)
    32      u_ship_max      ship maximum true wind speed (m/s)
    33      u_ship_std      ship true wind speed stdev (m/s)
    34      wdir_ship       ship mean true wind direction
    35      wdir_ship_min   ship minimum true wind direction
    36      wdir_ship_max   ship maximum true wind direction
    37      wdir_ship_std   ship true wind direction stdev
    38      rwspd_psd       streamwise relative wind speed (m/s)from PSD sonic
    39      Un              N component, corrected true windspeed, PSD sonic
    40      Ue              E component, corrected true windspeed, PSD sonic
    41      wdir_psd        true wind direction from N,E components, PSD sonic
    42      rwdir_psd       relative wind direction from PSD sonic
    43      azm             azimuth angle from first rotation, radians
    44      tilt            tilt angle from first rotation, radians
    45      Tsea_psd        seasnake SST, degC
    46      Tair_psd        air temperature, PSD vaisala, degC
    47      RH_psd          RH, PSD vaisala
    48      RH_ship         RH, ship sensor
    49      Q_psd           specific humidity, PSD sensors, g/kg
    50      Q_ship          specific humidity, ship sensors, g/kg
    51      Qsea_psd        sea surface specific humidity from seasnake SST
    52      rhoa_psd        moist air density from PSD sensors (kg/m3)
    53:73   yship           COARE 3.5 output, ship data:[usr,tau,hsb,hlb,hlwebb,
                         tsr,qsr,zot,zoq,Cd,Ch,Ce,L,zet,dter,dqer,tkt,RF,
                         Cdn_10,Chn_10,Cen_10]
    74:94   ypsd            COARE 3.5 output, PSD data, as above
    95      Rl_std_psd      std PIR flux, PSD
    96      Rs_std_psd      std PSP flux, PSD
    97      sog_ship        mean sog, ship gps, m/s
    98      sog_ship_std    std sog, ship gps
    99      sog_ship_min    minimum sog, ship gps
    100     sog_ship_max    maximum sog, ship gps
    101     cog_ship        mean cog, deg, ship gps
    102     cog_ship_std    std cog, ship gps
    103     cog_ship_min    minimum cog, ship gps
    104     cog_ship_max    maximum cog, ship gps
    105     cog_psd         mean cog, deg, PSD gps
    106     cog_psd_std     std cog, PSD gps
    107     cog_psd_min     minimum cog, PSD gps
    108     cog_psd_max     maximum cog, PSD gps
    109     sog_psd         mean sog, m/s, PSD gps
    110     sog_psd_min     minimum sog, PSD gps
    111     sog_psd_max     maximum sog, PSD gps
    112     lrg             mean heading, deg, ship gyro
    113     lrg_std         std heading, ship gyro
    114     lrg_min         minimum heading, ship gyro
    115     lrg_max         maximum heading, ship gyro
    116     no_data         wu, Hare correction, not computed
    117     wu              wu covariance, Edson correction
    118     no_data         wv, Hare correction, not computed
    119     wv              wv covariance, Edson correction
    120     no_data         wT, Hare correction, not computed
    121     wt              wT covariance, Edson correction
    122     no_data         uvar, Hare correction, not computed
    123     no_data         vvar, Hare correction, not computed
    124     no_data         wvar, Hare correction, not computed
    125     uvar            u variance, Edson correction
    126     vvar            v variance, Edson correction
    127     wvar            w variance, Edson correction
    128     tvar            T variance
    129     Sp1             U platform velocity std, Edson correction
    130     Sp2             V platform velocity std, Edson correction
    131     Sp3             W platform velocity std, Edson correction
    132     hd_sin_std      stdev sine lrg heading
    133     hd_cos_std      stdev cosine lrg heading
    134     Cu2a            Cu structure func parameter, method a
    135     Cu2b            Cu structure func parameter, method b
    136     Cw2a            Cw structure func parameter, method a
    137     Cw2b            Cw structure func parameter, method b
    138     Ct2a            Ct structure func parameter, method a
    139     Ct2b            Ct structure func parameter, method b
    140     Cq2a            Ch2o structure func parameter, method a
    141     Cq2b            Ch2o structure func parameter, method b
    142     usid_a          ustar from idiss, method a
    143     usid_b          ustar from idiss, method b
    144     tsid_a          tstar from idiss, method a
    145     tsid_b          tstar from idiss, method b
    146     qsid_a          qstar from idiss, method a, kg/kg
    147     qsid_b          qstar from idiss, method b, kg/kg
    148     Lid_a           L from idiss, method a
    149     Lid_b           L from idiss, method b
    150     lon             longitude, scs gps
    151     lat             latitiude, scs gps
    152     Ubar            streamwise true windspeed from Edson correction
    153     odec            odec (ship speed through water), m/s
    154     odec_std        stdev odec
    155     odec_min        minimum odec
    156     odec_max        maximum odec
    157     P_mb            atmospheric pressure, PSD sensor, mb
    158     missing         missing data lines in son file per 10 min
    159     badPnts         # bad data points (spikes, etc.) in sonic data per 10 min
    160     hs_cov          sensible heat flux from covariance
    161     hl_cov          latent heat flux from covariance, uncorrected
    162     hlc_cov         webb corrected hl_cov
    163     sal_scs         salinity, ship system, ppt
    164     Ts_noise        Sonic temperature noise variance
%}
        no_data = zeros(6,1)*NaN;
        StRow = hhh*6+1;
        EndRow = StRow + 5;
        da(StRow:EndRow,:) = [jd_10min,wspd_scs,ta_scs,rh_scs,bp_scs,...        % 1-5
            ts_scs,rs_scs,rl_scs,lat_scs,rain_scs,no_data,no_data,...           % 6-12
            no_data,no_data,no_data,wspd_psd,ta_psd,rh_psd,bp_psd,...           % 13-19
            ts_psd,rs_psd,rl_psd,lat_scs,rain_psd,no_data,no_data,...           % 20-26
            no_data,no_data,no_data,wspd_scs,scs_10min_min(:,13),...            % 27-31
            scs_10min_max(:,13),scs_10min_std(:,13),wdir_scs,...                % 32-34
            scs_10min_min(:,14),scs_10min_max(:,14),scs_10min_std(:,14),...     % 35-37
            rwspd_psd,Un_10min,Ue_10min,wdir_psd,rwdir_psd,azm',tilt',...       % 38-44
            ts_psd,ta_psd,rh_psd,rh_scs,qa_psd,qa_scs,qsea_psd,...              % 45-51
            rhoa_psd,coare_scs,coare_psd,met_10min_std(:,6),...             % 52-95
            met_10min_std(:,5),sog_scs,scs_10min_std(:,5),...                   % 96-98
            scs_10min_min(:,5),scs_10min_max(:,5),cog_scs,...                   % 99-101
            scs_10min_std(:,4),scs_10min_min(:,4),scs_10min_max(:,4),...        % 102-104
            gps_10min_avg(:,2),gps_10min_std(:,2),gps_10min_min(:,2),...        % 105-107
            gps_10min_max(:,2),gps_10min_avg(:,3),gps_10min_min(:,3),...        % 108-110
            gps_10min_max(:,3),lrg_scs,scs_10min_std(:,12),...                  % 111-113
            scs_10min_min(:,12),scs_10min_max(:,12),no_data,...                 % 114-116
            wu_cov_10min(:,2),no_data,wv_cov_10min(:,2),no_data,...             % 117-120
            wt_cov_10min(:,2),no_data,no_data,no_data,uvwt_10min_std(:,2).^2,...% 121-125
            uvwt_10min_std(:,3).^2,uvwt_10min_std(:,4).^2,...                   % 126-127
            uvwt_10min_std(:,5).^2,plat_10min_std(:,2:4),...                    % 128-131
            shead_10min_std(:,2),chead_10min_std(:,2),Cu2a,Cu2b,Cw2a,...        % 132-136
            Cw2b,Ct2a,Ct2b,Cq2a,Cq2b,usid_a,usid_b,tsid_a,tsid_b,qsid_a,...     % 137-146
            qsid_b,Lid_a,Lid_b,lon_scs,lat_scs,uvwt_10min_avg(:,2),...          % 147-152
            no_data,no_data,no_data,no_data,bp_psd,missing,badPnts,...          % 153-159
            hs_cov_10min,hl_cov_10min,hlc_cov_10min,sal_scs,Ts_noise];          % 160-164
%{
        10-min output array for CO2/H2O statistics: 144x17
        follows historical format, but CO2 covariance is not computed.

        cr column designations as follows:

        COLUMN  VARIABLE        DESCRIPTION
        1       doy             decimal day of year
        2       H2O_Lic_all     mean H2O, all data, mmol/m3
        3       H2O_Lic_good    mean H2O, good data, mmol/m3
        4       CO2_Lic_all     mean CO2, all data, mmol/m3
        5       CO2_Lic_good    mean CO2, good data, mmol/m3
        6       P_Licor         atm pressure from Licor box, kPa
        7       T_Licor         air temperature from Licor box, C
        8       wq              water vapor flux
        9       hlc             latent heat flux = wq*Rho*Cp
        10      not_used1       CO2 covariance, not computed
        11      not_used2       CO2 covariance, not computed
        12      std_H2O_all     std H2O, all data
        13      std_H2O_good    std H2O, good data
        14      std_CO2_all     std CO2, all data
        15      std_CO2_good    std CO2, good data
        16      P_std           std Licor box pressure
        17      Licor_AGC       mean diagnostic value
%}

        cr(StRow:EndRow,:) = [jd_10min,lic_10min_avg_all(:,2),...
                no_data,lic_10min_avg_all(:,3),no_data,lic_10min_avg_all(:,6),...
                lic_10min_avg_all(:,5),wq_cov_10min(:,2),hl_cov_10min,no_data,...
                no_data,lic_10min_std_all(:,2),no_data,lic_10min_std_all(:,3),...
                no_data,lic_10min_std_all(:,6),lic_10min_avg_all(:,7)];

%{
        10-min spec output arrays for spectra and cospectra: 144x481

        column designations as follows, all spectra are 48 bin

        COLUMNS     VARIABLE        DESCRIPTION
        1           doy             decimal day of year
        2:49        Fxx             frequency bins
        50:97       Puu             U spectrum
        98:145      Pvv             V spectrum
        146:193     Pww             W spectrum
        194:241     Ptt             T spectrum
        242:289     Pqq             q spectrum, LI7500
        290:337     Cuw             wu cospectrum
        338:385     Cvw             wv cospectrum
        386:433     Ctw             wt cospectrum
        434:481     Cqw             wq cospectrum
%}

        Fx2 = zeros(6,48);
        for ii=1:6; Fx2(ii,:) = Fx; end; % reshape Fx to 6x48 to match other arrays
        sp(StRow:EndRow,:) = [jd_10min, Fx2, Pu, Pv, Pw, Pt, Pq, Cuw, Cvw, Ctw, Cqw];

    end; % end of 24 hour loop

    %% save da, sp, cr files
    % reset first col to be sure missing hour timestamps are not NaN
    da = real(da);
    da(:,1) = (ddd:10/1440:ddd+1-10/1440);
    da_fname = ['da_',datasetv,'_',jd,'_',dcname,'.txt'];
    fl = fullfile(path_da,da_fname);
    save (fl,'da','-ascii','-tabs');

    cr(:,1) = (ddd:10/1440:ddd+1-10/1440);
    cr_fname = ['cr_',datasetv,'_',jd,'_',dcname,'.txt'];
    fl = fullfile(path_cr,cr_fname);
    save (fl,'cr','-ascii','-tabs');

    sp = real(sp);
    sp(:,1) = (ddd:10/1440:ddd+1-10/1440);
    sp_fname = ['sp_',datasetv,'_',jd,'_',dcname,'.txt'];
    fl = fullfile(path_sp,sp_fname);
    save (fl,'sp','-ascii','-tabs');

    %% daily plots

    % FILTERS
    % kk is general filter for rwdir, std heading, std sog, vplat variance,
    %     (sqrt(std_w)/ugw < .7*.0015*U^2) & (sqrt(std_v)/ugu < 1.4)
    % hh is filter for kk & licor AGC < 60
    % jj is additional filter for kk & isfinite(usid) & (abs(-wu-usid)^2<.02+.006*U)
    % mm; both hh and jj for hlid

    wtv = -da(:,74).*(da(:,79)+.61.*(da(:,17)+273.16).*da(:,80));
    ws = (9.83./(da(:,17)+273.16).*abs(wtv)*600).^.333;  % convective velocity
    ug = sqrt(da(:,74).^2+ws.^2);
    ugw = sqrt(1.5*da(:,74).^2+ws.^2);
    ugu = sqrt(8*da(:,74).^2+1.*ws.^2);
    gg = sqrt(1+1.6*ws.^2./da(:,16).^2);
    kk = find((da(:,56)<90 & da(:,56)>-90) & da(:,113)<5 &...
               da(:,98)<0.6 & da(:,130)<0.8 &...
               (sqrt(da(:,127))./ugw < .7+.0015*da(:,16).^2) &...
               (sqrt(da(:,126))./ugu<1.4));
    hh = intersect(kk, find(cr(:,17) < 60));
%     jj = intersect(kk, find(isfinite(da(:,142)) &...
%                   (abs(-da(:,117)-da(:,142).^2) < .02+.006*da(:,16))));
    jj = intersect(kk, find(isfinite(da(:,142))));
    mm = intersect(hh,jj);

    % Air/Sea Temperatures
    figure(1);
    subplot(2,1,1); plot(da(:,1),da(:,3),'bo-',da(:,1),da(:,17),'ro-');
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Air/Sea Temperatures',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlabel('Hour (UTC)');ylabel('Air Temp (C)'); 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));
        legend('Ship','PSD','Location','Best');
    subplot(2,1,2); plot(da(:,1),da(:,6),'bo-',da(:,1),da(:,20),'ro-');
        xlabel('Hour (UTC)');ylabel('SST (C)'); 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');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,['Air_Sea_Temps_',cruise,'_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % RH/q
    Rgas = 287.1;
    Rgas_universal = 8.314472; % Pa m3 /K /mol
    Licor_H2O_a=(cr(:,2)*1e-3).*Rgas_universal.*(cr(:,7)+273.15)./(cr(:,6)*1e3); %mmol/m3 to mol/mol
    Mw = 18.01528; % molar mass of H20 g/mol
    Md = 28.964;
    epsilon = Mw/Md; % mass of water to mass of dry air 0.622
    rhoa_lic = (cr(:,6)*1e3).*(1-Licor_H2O_a*(1-epsilon))./(Rgas.*(cr(:,7)+273.15));%air density in licor [kg/m3]
    Licor_H2O_m = (cr(:,2)*1e-3).*Mw./rhoa_lic; %mmol/m3 to g/kg
    licor_rh = relhum5([cr(:,7),Licor_H2O_m,cr(:,6)*10]);
    dd = find(cr(:,17) < 60);
    figure(2);
    subplot(2,1,1); plot(da(:,1),da(:,4),'bo-',da(:,1),da(:,18),'ro-',da(dd,1),licor_rh(dd),'go');
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  RH & Specific Humidity',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlabel('Hour (UTC)');ylabel('RH'); grid; ylim([40,100]);
        set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
        xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2));
        legend('Ship','PSD-Vais','Licor','Location','Best');
    subplot(2,1,2); plot(da(:,1),da(:,50),'bo-',da(:,1),da(:,49),'ro-',da(dd,1),Licor_H2O_m(dd),'go');
        xlabel('Hour (UTC)');ylabel('q (g/kg)'); 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');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,['RH_q_',cruise,'_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % True Wspd/Wdir
    figure(3);
    subplot(2,1,1); plot(da(:,1),da(:,30),'bo-',da(:,1),da(:,16),'ro-');
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  True Wind',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        xlabel('Hour (UTC)');ylabel('True Wind Speed (m/s)'); 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));
        legend('Ship','PSD','Location','Best');
    subplot(2,1,2); plot(da(:,1),da(:,34),'bo-',da(:,1),da(:,41),'ro-');
        xlabel('Hour (UTC)');ylabel('True Wind Dir'); grid; ylim([0 360]);
        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');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,['True_Wind_',cruise,'_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % Relative Wspd/Wdir
    figure(4);
    Uplim = ones(144,1)*60;
    LoLim = ones(144,1)*-60;
    subplot(2,1,1); plot(da(:,1),da(:,38),'ro-');
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Relative Wind',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        xlabel('Hour (UTC)');ylabel('Rel Wind Speed (m/s)'); 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));
    subplot(2,1,2); plot(da(:,1),da(:,42),'ro-',da(:,1),Uplim,'b-',da(:,1),LoLim,'b-');
        xlabel('Hour (UTC)');ylabel('Rel Wind Dir'); grid; ylim([-180 180]);
        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');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,['Rel_Wind_',cruise,'_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % Rain
    figure(5);
%     temp = da(:,10);
%     temp(isnan(temp)) = 0;
%     accm_scs = cumsum(temp)/6; % mm/hr to mm/10 min
    temp = da(:,24);
    temp(isnan(temp)) = 0;
    accm_psd = cumsum(temp)/6;
%    subplot(2,1,1); plot(da(:,1),da(:,10),'b-',da(:,1),da(:,24),'r-');
    subplot(2,1,1); plot(da(:,1),da(:,24),'r-');
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Rain',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        xlabel('Hour (UTC)');ylabel('Rain Rate (mm/hr)'); 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));
%        legend('Ship','PSD','Location','Best');
%    subplot(2,1,2); plot(da(:,1),accm_scs,'b-',da(:,1),accm_psd,'r-');
    subplot(2,1,2); plot(da(:,1),accm_psd,'r-');
        xlabel('Hour (UTC)');ylabel('Accumulated Rain (mm)'); 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');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,['Rain_',cruise,'_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % Std Dev Platform Motion & Tilt Angle
    figure(6);
    subplot(2,1,1); plot(da(:,1),da(:,129),'b-',da(:,1),da(:,130),'r-',da(:,1),da(:,131),'g-');
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Motion Std Dev & Tilt',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        xlabel('Hour (UTC)');ylabel('Motion Std Dev (m/s)'); grid;
        set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
        legend('U plat','V plat','W plat','Location','Best');
        xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2));
    subplot(2,1,2); plot(da(:,1),da(:,44)*180/pi,'b-');
        xlabel('Hour (UTC)');ylabel('Wind Tilt Angle (deg)'); 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');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,['Motion_Tilt_',cruise,'_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % Heat Fluxes
    hsid = -da(:,52)*1004.67.*da(:,144).*da(:,142); % idiss method a
    [Le_w, ~] = Le_water(da(:,5), da(:,163));
    hlid = -da(:,52).*Le_w.*da(:,146)*1000.*da(:,142) - da(:,78); % Webb corrected, g/kg
    figure(7);
    subplot(2,1,1); plot(da(:,1),da(:,56),'bo-',da(:,1),da(:,77),'ro-',da(hh,1),da(hh,162),'go',da(mm,1),hlid(mm),'mo');
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Heat Fluxes',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        xlabel('Hour (UTC)');ylabel('Latent Heat (W/m^2)'); grid;
        ylim([min(-5,min(da(:,77)-0.25*abs(min(da(:,77))))), max(50,max(da(:,77))*1.25)]);
        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('Bulk-scs','Bulk-psd','Cov-Webb','ID-a','Location','Best');
        subplot(2,1,2); plot(da(:,1),da(:,55),'bo-',da(:,1),da(:,76),'ro-',da(kk,1),da(kk,160),'go',da(kk,1),hsid(kk),'mo');
        xlabel('Hour (UTC)');ylabel('Sensible Heat (W/m^2)'); grid;
        ylim([min(-5,min(da(:,76))-0.25*abs(min(da(:,76)))), max(5,max(da(:,76))+0.5*abs(max(da(:,76))))]);
        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');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,['Heat_Fluxes_',cruise,'_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % Stress
    figure(8); plot(da(:,1),da(:,54),'bo-',da(:,1),da(:,75),'ro-',da(kk,1),-(da(kk,52).*da(kk,117)),'ko',da(jj,1),(da(jj,52).*da(jj,142).^2),'go',da(jj,1),(da(jj,52).*da(jj,143).^2),'mo');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Streamwise Stress',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    xlabel('Hour (UTC)');ylabel('Stress (N/m^2)'); grid; ylim([-0.02 max(0,max(da(:,75)))*1.25]);
    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('Bulk-scs','Bulk-psd','Cov','ID-A','ID-B','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 plotit
        ppath = fullfile(path_dailyPlots,['Stress_',cruise,'_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    % GPS Track Plot
    [lon2,~] = despike2(da(:,150));
    [lat2,~] = despike2(da(:,151));
    figure(9); plot(lon2,lat2,'k-.');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Track Plot',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    xlabel('Longitude');ylabel('Latitude'); grid;
    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 plotit
        ppath = fullfile(path_dailyPlots,['TrackPlot_',cruise,'_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end;

    close all;

    %% reinitialize output arrays and loop back to next day

    da = zeros(144,164)*NaN;
    sp = zeros(144,481)*NaN;
    cr = zeros(144,17)*NaN;

end; % end daily loop

% concatenate daily files into one whole-cruise file
filecat(fullfile(path_da,['da_',datasetv,'_*_',dcname,'.txt']), fullfile(path_proc_data_flux,[cruise,'_da_',datasetv,'_',dcname,'_c25.txt']),'Overwrite');
filecat(fullfile(path_cr,['cr_',datasetv,'_*_',dcname,'.txt']), fullfile(path_proc_data_flux,[cruise,'_cr_',datasetv,'_',dcname,'_c25.txt']),'Overwrite');
filecat(fullfile(path_sp,['sp_',datasetv,'_*_',dcname,'.txt']), fullfile(path_proc_data_flux,[cruise,'_sp_',datasetv,'_',dcname,'.txt']),'Overwrite');


