%{
run_motcorr.m

Includes hl / hs computed with Simon's method (wpl_schotanus.m).

Driver script for the flux calculations

Fluxes are processed in 10-min segments.

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

The additional flux utility functions usually don't require
project-specific modifications.

Set run parameters in the first block below as desired.

V2: May 2015 with updated and vectorized code, BRB.
- Motion corrections are Edson98 only - Hare code removed.
- CO2 flux FROM LI 7500 is not computed.
V4 : TSONIC & Q_licor DECORRELATED W/ RESPECT TO 3-AXIS MOTION.

%}

%% Initialize run parameters
close all;
fclose all;
clear;
warning ('off','MATLAB:MKDIR:DirectoryExists');

% cruise ID and paths
cruise = 'PISTON_2019';  % cruise name
cruise_str = 'PISTON';  % cruise name without _
ship = 'Sally';             % research vessel
plotit = true;                   % create plots and save to .png format
yyyy = '2019';

% Dates
% for atomic, use 9 to 44
% for piston 2019 use 248 to 268
stjd = 248; endjd = 248;

%%% testing increase in wind speed based on suspected flow distortion.
%%% gets multiplied on wspd, which goes into bulk fluxes
%%% gest multiplied on rspd, which goes into ID fluxes

perc_inc = 98;
if perc_inc == 0        %% none
    folder_suffix = '';
    u_mult = 1.0;
elseif perc_inc == 10   %% 10% increase
    folder_suffix = '_uinc_10';
    u_mult = 1.10;
elseif perc_inc == 15   %% 15% increase
    folder_suffix = '_uinc_15';
    u_mult = 1.15;
elseif perc_inc == 30   %%% 30% increase
    folder_suffix = '_uinc_30';
    u_mult = 1.30;
elseif perc_inc == 98
    folder_suffix = '_uinc_var2_check';
elseif perc_inc == 99
    folder_suffix = '_uinc_var';
end

% system specific path defs
sysType = computer;
username=char(java.lang.System.getProperty('user.name'));
if strncmp(sysType,'MACI64',7)     % set Mac paths
    data_drive = '/Users/eliz/DATA/';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs');
elseif strncmp(sysType,'PCWIN64',7)  % set PSL DAC paths
    data_drive = 'D:\DATA\';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs');
end

% matlab script path
restoredefaultpath
cd(fullfile(path_prog,'flux'));
addpath(genpath(fullfile(path_prog,'flux')));
rehash toolboxcache;

% select lat/lon limits and adjustments based on leg number
Lonmin = 119; Lonmax = 139; Latmin = 4; Latmax = 23;
td1_adj = 0;        % ?
tc1_adj = 0;        % ?
td2_adj = 0;        % ?
tc2_adj = 0;        % ?
td_s_adj = 0;     % ?
tc_s_adj = 0;     % ?
ts_adj = 0;       % sea snake
ta_adj = 0;         % air temp
ta_s_adj = 0;  % ship air T
ta_vais_adj = 0;  % psl vaisala air T
wxt_adj = 0;

PosLims = [Lonmin, Lonmax, Latmin, Latmax];
adj = [ts_adj,ta_vais_adj,td1_adj,tc1_adj,td2_adj,tc2_adj,wxt_adj];
SCS_adj = [td_s_adj,tc_s_adj];

%% types of processing: motion corrected or not, using WXT or not
use_wxt = 0;
sens_sep_factor = 1.0; % this was 1.045 in a prior set up... not sure what is best here. 
decorr = true;      % set true decorrelate winds & motion...
if decorr
    dcname = 'decorr';
else
    dcname = 'std';
end

%% 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');
path_proc_data_flux = fullfile(data_drive,cruise,ship,'flux',['Processed' folder_suffix]);
path_proc_images_flux = fullfile(data_drive,cruise,ship,'flux',['Processed_Images' folder_suffix]);

path_output = fullfile(data_drive,cruise,ship,'flux',['Field_Processed' folder_suffix],'flux_met_sea_10min_v0');
path_output_final = fullfile(data_drive,cruise,ship,'flux',['Field_Processed' folder_suffix],'flux_met_sea_10min');
mkdir(path_output);
mkdir(path_output_final);

% subdirectory for saving daily met/flux plots
path_dailyPlots = fullfile(path_proc_images_flux,['Daily_',dcname]);
mkdir(path_dailyPlots);
mkdir(fullfile(path_dailyPlots,'Air_Sea_Temps'));
mkdir(fullfile(path_dailyPlots,'SOG_COG'));
mkdir(fullfile(path_dailyPlots,'Motion_Tilt'));
mkdir(fullfile(path_dailyPlots,'Radiation'));
mkdir(fullfile(path_dailyPlots,'Rain'));
mkdir(fullfile(path_dailyPlots,'Rel_Wind'));
mkdir(fullfile(path_dailyPlots,'RH'));
mkdir(fullfile(path_dailyPlots,'Heat_Fluxes'));
mkdir(fullfile(path_dailyPlots,'Stress'));
mkdir(fullfile(path_dailyPlots,'Track_Plot'));
mkdir(fullfile(path_dailyPlots,'True_Wind'));
% subdirectory for saving corrected wind timeseries plots
path_windPlots = fullfile(path_proc_images_flux,['WindPlots_',dcname]);
mkdir(path_windPlots);
mkdir(fullfile(path_windPlots,'motcorr_segment'));
mkdir(fullfile(path_windPlots,'uvw'));
% subdirectory for saving corrected wind spectra plots
path_spectraPlots = fullfile(path_proc_images_flux,['SpectraPlots_',dcname]);
mkdir(path_spectraPlots);
mkdir(fullfile(path_spectraPlots,'q_spectrum'));
mkdir(fullfile(path_spectraPlots,'t_spectrum'));
mkdir(fullfile(path_spectraPlots,'uv_spectrum'));
mkdir(fullfile(path_spectraPlots,'w_spectrum'));
mkdir(fullfile(path_spectraPlots,'wq_cospectrum'));
mkdir(fullfile(path_spectraPlots,'wt_cospectrum'));
mkdir(fullfile(path_spectraPlots,'wu_cospectrum'));
% subdirectory for saving corrected wind data
path_uvwStreamData = fullfile(path_proc_data_flux,['uvwStream_',dcname]);
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]);
mkdir(path_da);
% subdirectory for saving hourly sp files
path_sp = fullfile(path_proc_data_flux,['sp_',dcname]);
mkdir(path_sp);
% subdirectory for saving hourly cr files
path_cr = fullfile(path_proc_data_flux,['cr_',dcname]);
mkdir(path_cr);
% subdirectory for saving hourly lag times
path_lag = fullfile(path_proc_data_flux,['lag_',dcname]);
mkdir(path_lag);
% subdirectory for saving hourly decorr coefs
path_mu = fullfile(path_proc_data_flux,['mu_',dcname]);
mkdir(path_mu);
mkdir(fullfile(path_mu,'w_Ts'));
mkdir(fullfile(path_mu,'q'));
% initialize plot format
graphformat = '.png';  % select graphics format
graphdevice = '-dpng'; % select graphic device

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

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

%% fetch variables and strings related to this particular cruise
setup_cruise;


%% constants
d2r = pi/180;
r2d = 180/pi;

C2K = 273.15;               % T conversion constant
d2r = pi/180;               % angle conversion constants
r2d = 180/pi;
Rgas = 287.1;               % Pa m-3 kg-1 K-1 for dry air
Rgas_universal = 8.314472;	% Pa m3 K-1 mol-1
Mw = 18.01528;              % molar mass of H2O g/mol
Md = 28.964;                % molar mass dry air
epsilon = Mw/Md;            % mass of water to mass of dry air = 0.622... epsilon = 0.9715; % mean emissivity of the ocean
cpa  = 1004.67;             % heat capacity air
sigma = 5.67E-8;            % stephan boltzmann constant W m-2 K-4

%% other ship data
% load ROSR data as structure "sea"
load('/Users/eliz/DATA/PISTON_2019/Sally/SST/ROSR/all_data/seawater_PISTON_2019.mat');
% load rest of ship data as structure "MET"
load('/Users/eliz/DATA/PISTON_2019/Sally/ship/MET_SallyRide_PISTON_2019.mat');

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

    % 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,184);
    sp = NaN(144,481);
    cr = NaN(144,17);
    lag = NaN(24,2); lag(:,1) = ddd:1/24:ddd+1-1/24;
    mu_decorr = NaN(144,17); mu_decorr(:,1) = ddd:10/1440:ddd+1-10/1440;
    mu_q_decorr = NaN(144,13); mu_q_decorr(:,1) = ddd:10/1440:ddd+1-10/1440;

    if strcmp(cruise,'PISTON_2019') == 1 && ddd == 268
        hr_end = 11;  %%% there are only 11 hours in day 268 for PISTON
    else
        hr_end = 23;
    end
    
    for hhh = 00 : hr_end   % iterate over hours
        %% read raw data files...
        
        hr = sprintf('%02i',hhh);

        % read psl 1 min met files: met1, met2, wxt
        dfl1 = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['met1',yyyy(3:4),jd,hr,'_raw.txt']);
        dfl2 = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['met2',yyyy(3:4),jd,hr,'_raw.txt']);
        if use_wxt == 0
            met = read_met(dfl1,dfl2,ddd,hhh,adj,zp);
            % add in nans for WXT data that doesn't exist and wasn't read in
            met = [met nan(60,8)];       
        elseif use_wxt == 1
            dfl3 = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['wxt0',yyyy(3:4),jd,hr,'_raw.txt']);
            wxt = read_wxt(dfl3, ddd, hhh, zwxt);
            % Add wxt to met, but don't include the same 1-min time array
            met = [met; wxt(:,(2:8))]; 
        end
        
                % 10-min psl met statistics is 6x25 array:
        %{
            1   jd_met          1 min timestamp
            2   Tvais           T vaisala, C
            3   Rhvais          RH vaisala
            4   Tsnk            sea snake surface SST
            5   psp1            psp1 radiation, W/m2
            6   pir1            pir1 radiation, W/m2
            7   psp2            psp2 radiation, W/m2
            8   pir2            pir2 radiation, W/m2
            9   Tc1             pir1 case temperature
            10  Td1             pir1 dome temperature
            11  Tc2             pir2 case temperature
            12  Td2             pir2 dome temperature
            13  org             rain, mm/hr - none for this cruise
            14  slp             air slp, mb
            15  aspir_on        T/RH diagnostic
            16  org_carrier     org diagnostic
            17  org_V           ORG raw measurement, Volts
            18  therm1          thermopile reading PIR1 W/m2
            19  therm2          thermopile reading PIR2 W/m2
            20  P               measured pressure, mb            
            21  rdir_wxt       rel wind direction ,0-360 deg from bow
            22  rspd_wxt       rel wind speed, m/s
            23  Ta_wxt        air temperature, C
            24  rh_wxt          RH, %
            25  slp_wxt         air pressure, mb or hPa
            26  prate_wxt        rain rate, mm/hr
            27  U_wxt           U wind component, WXT
            28  V_wxt           V wind component, WXT
        %}
            
        % read motion file
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['mot0' yyyy(3:4),jd,hr,'_raw.txt']);
        mot = read_motion_real(dfl,ddd,hhh);

        % read sonic file
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['son0',yyyy(3:4),jd,hr,'_raw.txt']);
        [son,badPnts,missing] = read_sonic(dfl,ddd,hhh,sonicmodel,rotationsonic);
        badPnts = badPnts(:); % # NaNs in raw file (water or ice on transducers, etc.)
        missing = missing(:); % # missing data points in raw file - timestamp gaps
        xxx = missing>100 | badPnts>100;  % bad sonic data

        % read gprm file - PYTHON processed output from raw Hemosphere gps file
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['gprm',yyyy(3:4),jd,hr,'_raw.txt']);
        gps = read_gprm(dfl,ddd,hhh,PosLims);

        % read Hemisphere heading and roll angle
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['hed0',yyyy(3:4),jd,hr,'_raw.txt']);
        hed = read_hed_pitch(dfl,ddd,hhh);

        % read licor
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['lic0',yyyy(3:4),jd,hr,'_raw.txt']);
        lic = read_licor(dfl,ddd,hhh);

        % read 1 Hz ship data
        dfl = fullfile(path_raw_data_flux,[yyyy(3:4),jd],['scs0',yyyy(3:4),jd,hr,'_raw.txt']);
        
        %%% FIX - get the data that are needed just to run this code out of
        %%% structure. then worry about the rest later. Run it as is for
        %%% now. 
		ship_1min = read_scs_v2(dfl,ddd,hhh,PosLims,SCS_adj,zp_ship);
        
        % FIX: write an array of expected values of type scs from scs_hr
        % structure
        
%         ship_nan = nan(3600,1);
%         
%         ship = [ship_hr.jd, ship_hr.lat, ship_hr.lon, ship_hr.cog, ship_hr.sog, ...  % 1-5
%             ship_nan, ship_nan, ship_hr.rh, ship_hr.Ta, ship_hr.sw_dn, ....
%             ship_nan, ship_hr.hed, ship_hr.wspd, ship_hr.wdir, ship_hr.lw_dn, ...
%             ship_hr.slp, ship_hr.qa, ship_hr.paccum, ship_nan, ship_nan, ...
%             ship_nan, ship_nan, ship_hr.rspd, ship_hr.rdir, ship_hr.sogN, ...
%             ship_hr.sogE, ship_nan];

%%% expected order of an ship array size 86400 x 27 ... :
%     1    jd = Decimal DOY
%     2    lat
%     3    lon
%     4    cog      deg
%     5    sog      m/s
%     6    Ttsg     TSG temperature
%     7    Stsg     TSG salinity
%     8    rh       %
%     9    Ta       air temp, C
%     10   sw_dn    downwelling solar radiation, W/m2
%     11   prate    rain rate, mm/hr
%     12   hed      heading, deg
%     13   wspd     true wind speed, m/s
%     14   wdir     true wind direction
%     15   lw_dn    downwelling IR, W/m2
%     16   slp      sealevel pressure, mb
%     17   qa       air specific humidity, g/kg
%     18   NAN
%     19   NAN
%     20   NAN
%     21   NAN
%     22   NAN
%     23   rspd     relative wind speed, m/s
%     24   rdir     relative wind dir, +/- 180 deg from bow
%     25   sogN     N ship speed component, m/s
%     26   sogE     E ship speed component, m/s
%     27   NAN 
%     28   un       *** added here: comp of true wind + south to north, like met convention v
%     29   uw       *** added here: comp of true wind + east to west, like met convention -u
%     30   hedN     *** added in this program, N comp of hed
%     31   hedE     *** added in this program, E comp of hed
%     32   run      *** added here: comp of rel wind + south to north, like met convention v
%     33   ruw      *** added here: comp of rel wind + east to west, like met convention -u
        
        %% copy out selected variables
        jd_1min = met(:,1);   % ref 1 min timestamp
        jd_1Hz = ship_1min.jd;    % ref 1 Hz timestamp
        jd_10Hz = son(:,1);   % ref 10 Hz timestamp

        % 1 sec ship met variables used for motion correction
%         lat = ship(:,2);      % not used
%         lon = ship(:,3);      % not used
%         cog_s = ship(:,4);
% %         sogm = ship(:,5);     % not used
% %         sogNm = ship(:,25);   % not used
% %         sogEm = ship(:,26);   % not used
%         hed_s = ship(:,12);
        
        % 1-min psl met variables



        t_1min = datenum(str2double(yyyy),0,0,0,0,0) + jd_1min;

        %% correct ORG voltage offset
        org_V = met(:,17);                              % output voltage
        org_bkgd = nanmedian1(despike2(org_V));     % median zero offset
        org_adj = 0.06484 - org_bkgd;                   % adjustment to desired offset of 0.06484 V
        org_rate_raw = (org_V+org_adj).^1.87*25 - 0.15; % recompute rain rate
        org_rate = org_rate_raw;                        % use that as first guess
        org_rate(org_rate<0.05)= 0;                     % remove data that is noisy ~ rates < 0.05 mm/hr threshold

        %% compute 10 Hz heading
        heading_interp = replace_NaN_nearest_neighbor(hed(:,2)*d2r);  % to radians
        heading = heading_interp;
        wh_heading_missing = find(isnan(heading_interp) ==1);
        wh_heading_valid = find(isfinite(heading_interp) ==1);
%         heading(wh_missing_hed) = [];
        shed = sin(heading);    % 10Hz components
        ched = cos(heading);

        %% E98 motion corrections, with decorr if desired
        % motcorr can handle missing heading because there is a check in it
        % before it runs angles.m, which adjusts the length of all the
        % arrays to still run properly.
        % use motcorr3, decorr in 10 min segments
%         [uvw,T_10Hz,accplat,uvwplat,xyzplat,euler,plat_rate,lagPnts,mu,uvw_raw,uvw_motcorr,f1,ff2] = ...
%                 motcorr3(son(wh_missing_hed,:),mot(wh_missing_hed,:),heading,sens_disp,fsonic,decorr);
        [uvw,T_10Hz,accplat,uvwplat,xyzplat,euler,plat_rate,lagPnts,mu,uvw_raw,uvw_motcorr,f1,ff2] = ...
                motcorr3(son,mot,heading,sens_disp,fsonic,decorr);

        if ~all(isnan(uvw(:,1))) && plotit && (hhh==0 || hhh==6 || hhh==12 || hhh==18)
            ppath = fullfile(path_windPlots,'motcorr_segment',['uvw_motcorr','_',jd,'_',hr,'_',dcname,graphformat]);
            print(f1,graphdevice,ppath);
            ppath = fullfile(path_windPlots,'uvw',['uvw_',jd,'_',hr,'_',dcname,graphformat]);
            print(ff2,graphdevice,ppath);
            close all;
        else
            close all;
        end

        uvw = uvw';                    % uvw - Edson motcorr and decorr
        T_10Hz = T_10Hz';              % this is sonic temperature
        accplat = accplat';            % platform acc,vel and disp in earth frame
        uvwplat = uvwplat';
        xyzplat = xyzplat';
        euler = euler';
        plat_rate = plat_rate';
        uvw_raw = uvw_raw';            % sonic uvw rotated to earth frame
        uvw_motcorr = uvw_motcorr';    % uvw with Edson correction only
        lagmat = [son(1,1), lagPnts];

        %% 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);

        % use filtered psl gps velocity @ 10 Hz for speed correction
        sogN = gps(:,4);
        sogE = gps(:,5);

        % Other option: interp ship velocity to 10 Hz for speed correction
%         sogN = interp1(jd_1min_s,sogNm,jd_10Hz,'linear','extrap');
%         sogE = interp1(jd_1min_s,sogEm,jd_10Hz,'linear','extrap');

        % 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) + sogN;
        uvw_rot(:,2) = uvw_rot(:,2) - sogE;

        % 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
        jd_10bin = (jd_start:10/1440:jd_end)';

        % Compute 10-min psl gps statistics: 6x5 array as follows:
        % 1:jd, 2:cog, 3:sog, m/s, 4:sogN, 5: sogE
%         gps_10min_avg = interval_avg(gps(:,1), gps(:,2:end),jd_10bin);
%         gps_10min_std = interval_std(gps(:,1), gps(:,2:end),jd_10bin);
%         gps_10min_min = interval_min(gps(:,1), gps(:,2:end),jd_10bin);
%         gps_10min_max = interval_max(gps(:,1), gps(:,2:end),jd_10bin);

        xavg = interval_avg(gps(:,1), gps(:,2:end),jd_10bin);
        cog_avg = xavg(:,2); sog_avg = xavg(:,3); sogN_avg = xavg(:,4); sogE_avg = xavg(:,5);
        xstd = interval_std(gps(:,1), gps(:,2:end),jd_10bin);
        cog_std = xstd(:,2); sog_std = xstd(:,3); sogN_std = xstd(:,4); sogE_std = xstd(:,5);
        xmin = interval_min(gps(:,1), gps(:,2:end),jd_10bin);
        cog_min = xmin(:,2); sog_min = xmin(:,3); sogN_min = xmin(:,4); sogE_min = xmin(:,5);
        xmax = interval_max(gps(:,1), gps(:,2:end),jd_10bin);
        cog_max = xmax(:,2); sog_max = xmax(:,3); sogN_max = xmax(:,4); sogE_max = xmax(:,5);
        
        % redo sog from speed components
        sog_avg = sqrt(sogN_avg.^2 + sogE_avg.^2);
        % redo cog from speed components, then unwrap
        cog_avg = atan2(sogE_avg,sogN_avg)*r2d;
        cog_avg = mod(cog_avg+360, 360);
        % redo std of cog from 10 Hz data
        cog_10Hz = gps(:,2);
        cog_std = interval_std_var(gps(:,1),unwrap(cog_10Hz*d2r),jd_10bin)*r2d;

        % 10-min heading
        ched_avg = interval_avg(hed(:,1), ched, jd_10bin);
        shed_avg = interval_avg(hed(:,1), shed, jd_10bin);
        [hed_avg,~] = cart2pol(ched_avg(:,2), shed_avg(:,2));
        hed_avg = hed_avg*r2d;
        hed_avg = mod(hed_avg+360, 360);
        
        % 10 min stats of met data
%         met_10min_avg = interval_avg(jd_1min, met(:,2:end),jd_10bin);
%         met_10min_std = interval_std(jd_1min, met(:,2:end),jd_10bin);
%         met_10min_min = interval_min(jd_1min, met(:,2:end),jd_10bin);
%         met_10min_max = interval_max(jd_1min, met(:,2:end),jd_10bin);
      
        % put met into a structure called mett
        mfields = {'jd';'Ta';'rh';'Tsnk';'psp1';'pir1';'psp2';'pir2';'Tc1';'Td1';'Tc2';'Td2';...  % 1-12
            'org';'slp';'aspir_on';'org_carrier';'org_V';'therm1';'therm2';'P';...                   % 13-20
            'rdir_wxt';'rspd_wxt';'Ta_wxt';'rh_wxt';'slp_wxt';'prate_wxt';'run_wxt';'ruw_wxt'}  ;        % 21-28
        
        for i = 1:length(mfields)
            mett.(mfields{i}) = met(:,i);
        end

        clear met_avg met_std met_min met_max
        for j = 1:length(mfields) 
            x = mett.(mfields{j});
            met_avg.(mfields{j}) = interval_avg_var(jd_1min, x, jd_10bin);
            met_std.(mfields{j}) = interval_std_var(jd_1min, x, jd_10bin);
            met_min.(mfields{j}) = interval_min_var(jd_1min, x, jd_10bin);
            met_max.(mfields{j}) = interval_max_var(jd_1min, x, jd_10bin);
        end

        %% Define a few 10-min PSL met variables for easier use below
        jd_10min = jd_10bin(1:end-1);
        t_10min = datenum(str2double(yyyy),0,0,0,0,0) + jd_10min;
         
        Ta = met_avg.Ta;
        rh = met_avg.rh;
        Tsnk = met_avg.Tsnk;
        P = met_avg.P;                  % measured pressure, mb
        slp = met_avg.slp;      % sea level pressure, mb

        %%% solar - choose how to do it based on sensor performance; 
        %%% from 1 sensor or a combination of 2
        rs = -(met_avg.psp1 + met_avg.psp2)/2; % if averaging solar sensors together
        %%% rs should be positive according to all the read statements that come before it.
        rs(rs<0) = 0.0;

        %%% IR - choose how to do it based on sensor performance;
        %%% from 1 sensor or a combination of 2
        rl = (met_avg.pir1 + met_avg.pir2)/2; % if averaging IR sensors together
        %%% rl should be positive according to all the read statements that come before it.
        rl(rl>500) = NaN;
        
        %%% rain at 10 min frequemncy, calculated from 1 min data
        rain = interval_avg_var(jd_1min,org_rate,jd_10bin);
        rain(rain<0.03) = 0;     % threshold of 0.03 mm/hr
                
        %% ROSR skin T
        Tskin = interp1(sea.t, sea.SSST, t_10min,'linear');
        
        % NOTE: load ROSR data
        
        %% remaining ship met data
        % [C,IA,IC] = unique(A) also returns index vectors IA and IC such that C = A(IA)
        [u_dnum, ID_u] = unique(MET.dnum);

        MET_out_vars = {'Ttsg1';'Ttsg2';'Stsg1';'Stsg2';...
            'Ftsg1';'Ftsg2';...
            'ID_Ftsg1';'ID_Ftsg2';'spd_log_u';'spd_log_v'};

        MET_in_vars = {'TSG_temperature_1';'TSG_temperature_2';'salinity_1';'salinity_2';
            'USW_flow_meter_1';'USW_flow_meter_2';...
            'TSG_goodflow1';'TSG_goodflow2';'ship_speed_transverse';'ship_speed_longitudinal'};

        for i = 1:length(MET_in_vars)
            thisvar = MET.(MET_in_vars{i});
            outvar = interp1(MET.dnum(ID_u), thisvar(ID_u), t_10min,'linear');
            eval([MET_out_vars{i} ' = outvar;']); 
        end

        ID_tsg  = nan(length(t_10min),1);      % 1 = use bow, 2 = use sea chest, ship, 0 if shouldn't be used at all.

        %%% on ??? day / hour / min, we switched from bow water to seachest water.
        %%% You can see the change in the flow rates and the data output. Use
        %%% seachest after this date. Use bow before.
        tsgswitchtime1 = datenum(2019,9,11, 23, 40, 0);
        tsgswitchtime2 = datenum(2019,9,12, 01, 20, 0);

        % remove bad data from bubbles
        Stsg1(Stsg1 < 30) = nan;
        Stsg2(Stsg2 < 30) = nan;

        ID_tsg(t_10min < tsgswitchtime1) = 1;  % use bow before switch
        ID_tsg(t_10min > tsgswitchtime2) = 2;  % use seachest after switch
        ID_tsg(t_10min >= tsgswitchtime1 & t_10min <= tsgswitchtime2) = 0;

        Ttsg1(ID_tsg == 0) = nan;
        Ttsg2(ID_tsg == 0) = nan;
        Stsg1(ID_tsg == 0) = nan;
        Stsg2(ID_tsg == 0) = nan;
        Ftsg1(ID_tsg == 0) = nan;
        Ftsg2(ID_tsg == 0) = nan;

        %%% ship precip ... my best but still not great attempt at rain rate
        paccum_s = interp1(MET.dnum1min, MET.paccum1min, t_10min);
        prate_s = cum2rate(paccum_s', t_10min')';

        %%% all the MET variables we want to save
        MET_vars = {'Ttsg1'; 'Ttsg2'; 'Stsg1'; 'Stsg2'; 'Ftsg1'; 'Ftsg2'; 'ID_Ftsg1'; 'ID_Ftsg2'; 'ID_tsg';...
        'spd_log_u'; 'spd_log_v';...
        'prate_s'; 'paccum_s'};


        %% work on met data from wxt: 
        
        % rain rate
        met_avg.prate_wxt(met_avg.prate_wxt < 0) = 0;
        
        % calculate true wind, and recalculate relative wind, from the averaged relative wind components
        % ... NOTE: these are NOT u and v components in meteorological convention
        [met_avg.wspd_wxt,met_avg.wdir_wxt,~,~] = uv_rel_to_sd_true(met_avg.run_wxt,met_avg.ruw_wxt,hed_avg,cog_avg,sog_avg);
        [met_avg.rspd_wxt,met_avg.rdir_wxt] = uv_to_sd(met_avg.run_wxt,met_avg.ruw_wxt);
        met_avg.rdir_wxt(met_avg.rdir_wxt>180) = met_avg.rdir_wxt(met_avg.rdir_wxt>180)-360;
        
        
        %% 10-min stats for ship data
        %%%% start with 1 min data
        jd_1min_s = ship_1min.jd;

        % NOTE: we already have relative wind direction and speed from ship
        % data, rdir and rspd. Here we adjust the winds for suspected flow
        % distortion by U + 5%, V - 15%
        % NOTE: are rspd2_s and ship_1min.rspd the same? If so, remove. 
        % NOTE: the u and v components here are NOT in meteorological
        % convention. They are positive to north (un) and to west (uw)
        [ship_1min.rspd2,ship_1min.rdir2,run_s,ruw_s] = sd_true_to_uv_rel(ship_1min.wspd,ship_1min.wdir,ship_1min.hed,ship_1min.cog,ship_1min.sog);
        % apply flow distortion correction to mean SCS rel wind speed components
        ship_1min.run = run_s/0.95; % U + 5%
        ship_1min.ruw = ruw_s/1.15; % V - 15%
        % recompute SCS true wind and true wind components from adjusted
        % relative u and v components. % NOTE: the u and v components here 
        % are NOT in meteorological convention. They are positive to north (un) and to west (uw)
        % recompute SCS true wind and true wind components
        [ship_1min.wspd2,ship_1min.wdir2,ship_1min.un,ship_1min.uw] = ...
            uv_rel_to_sd_true(ship_1min.run,ship_1min.ruw,ship_1min.hed,ship_1min.cog,ship_1min.sog);
        % use these corrected values by reassaigning wspd = wspd2 and wdir = wdir2
        ship_1min.wspd = ship_1min.wspd2;
        ship_1min.wdir = ship_1min.wdir2;
                
%         % original
%         [imrwspd,imrwdir,imUrel,imVrel] = sd_true_to_uv_rel(scs(:,13),scs(:,14),scs(:,12),scs(:,4),scs(:,5));
%         % apply flow distortion correction to mean SCS rel wind speeds
%         imUrel = imUrel/0.95; % U + 5%
%         imVrel = imVrel/1.15; % V - 15%
%         [imum2,imdm2,imU,imV] = uv_rel_to_sd_true(imUrel,imVrel,scs(:,12),scs(:,4),scs(:,5));
%         scs(:,13) = imum2;
%         scs(:,14) = imdm2;
        
        % compute N/E heading components
        [ship_1min.hedN, ship_1min.hedE] = pol2cart(ship_1min.hed*d2r,1);
        
        % these variables were added to former ship array as columns 28-33
%         ship_1min = [ship_1min,u,v,hedN,hedE,ru,rv];

        %% stats for everything at 10 min 
        
        clear ship_avg ship_std ship_min ship_max

        %%% for averaging 1-min averages to 10 min... probably not advisable.
        sfields = fields(ship_1min);

        for j = 1:length(sfields) 
            x = ship_1min.(sfields{j});
            ship_avg.(sfields{j}) = interval_avg_var(jd_1min_s, x, jd_10bin);
            ship_std.(sfields{j}) = interval_std_var(jd_1min_s, x, jd_10bin);
            ship_min.(sfields{j}) = interval_min_var(jd_1min_s, x, jd_10bin);
            ship_max.(sfields{j}) = interval_max_var(jd_1min_s, x, jd_10bin);
        end
    
        
%         ship_10min_avg = interval_avg(jd_1min_s, ship_1min(:,2:end),jd_10bin);
%         ship_10min_std = interval_std(jd_1min_s, ship_1min(:,2:end),jd_10bin);
%         ship_10min_min = interval_min(jd_1min_s, ship_1min(:,2:end),jd_10bin);
%         ship_10min_max = interval_max(jd_1min_s, ship_1min(:,2:end),jd_10bin);
%         ship_10min_avg(xxx,28:29) = NaN; % if psl sonic is bad, assume same for ship
%         ship_10min_std(xxx,28:29) = NaN;

        %%% FIX
        % recompute sog/cog from components
        ship_avg.cog = atan2(ship_avg.sogE, ship_avg.sogN)*r2d; % cog from components
        ship_avg.cog = mod(ship_avg.cog+360, 360);
        ship_avg.sog = sqrt(ship_avg.sogN.^2 + ship_avg.sogE.^2); % sog from components
        % recompute std cog from 1 min data
        ship_std.cog = interval_std_var(jd_1min_s, unwrap(ship_1min.cog*d2r),jd_10bin)*r2d; 
        % redo heading from components and std heading from 1 min data; same reason
        ship_avg.hed = atan2(ship_avg.hedE, ship_avg.hedN)*r2d; % hed from components
        ship_avg.hed = mod(ship_avg.hed+360, 360);
        ship_std.hed = r2d*(interval_std_var(jd_1min_s,unwrap(ship_1min.hed*d2r),jd_10bin));
        % redo wdir from components and std from 1 Hz data
        [ship_avg.wspd,ship_avg.wdir] = uv_to_sd(ship_avg.un,ship_avg.uw);
        % redo wdir std with 1 min adjusted data
        ship_std.wdir = r2d*(interval_std_var(jd_1min_s,unwrap(ship_1min.wdir*d2r),jd_10bin));
        % redo rdir from components, then unwrap
        [ship_avg.rspd,ship_avg.rdir] = uv_to_sd(ship_avg.run,ship_avg.ruw);
        clear xx; xx = ship_avg.rdir>180; 
        ship_avg.rdir(xx) = ship_avg.rdir(xx)-360;
        % redo std rdir and std rspd from 1 min data
        ship_std.rspd = interval_std_var(jd_1min_s,ship_1min.rspd,jd_10bin); % NOTE: use rspd or rspd2?... use adjusted or not? rspd2 isn't adjusted
        ship_std.rdir = r2d*(interval_std_var(jd_1min_s,unwrap(ship_1min.rdir*d2r),jd_10bin)); % NOTE: use rdir or rdir2?... use adjusted or not? rdir2 isn't adjusted

        %%% FIX
        % Define a few 10-min variables for use below
        Ttsg_s = Ttsg2;  % there are 2 TSGs so I had to enter this a different way. otherwise, Ttsg_s = ship_avg.Ttsg;  
        Stsg_s = Stsg2; % there are 2 TSGs so I had to enter this a different way. otherwise, Stsg_s = ship_avg.Stsg;  
        ship_avg.S = nan(6,1); % NOTE: I didn't do this yet
        ship_max.S = nan(6,1); % NOTE: I didn't do this yet
        ship_min.S = nan(6,1); % NOTE: I didn't do this yet
        ship_std.S = nan(6,1); % NOTE: I didn't do this yet
        
        lat_s = ship_avg.lat;
        lon_s = ship_avg.lon;
        cog_s = ship_avg.cog;
        sog_s = ship_avg.sog;
        rh_s = ship_avg.rh; % rh_s(rh_s<70 | rh_s>100) = NaN;
        Ta_s = ship_avg.Ta; % Ta_s(Ta_s<20 | Ta_s>35) = NaN;
        rs_s = ship_avg.sw_dn;
        rain_s = nan(6,1); % No ship rain data... otherwise rain_s = ship_avg.prate;
        hed_s = ship_avg.hed;
        wspd_s = ship_avg.wspd;
        wdir_s = ship_avg.wdir;
        rspd_s = ship_avg.rspd;
        rdir_s = ship_avg.rdir;
        rl_s = ship_avg.lw_dn; % rl_s(rl_s<380 | rl_s>500) = NaN;
        slp_s = ship_avg.slp;
        rosr_s = nan(6,1); % ship didn't have rosr... we brought our own. otherwise rain_s = ship_avg.rosr;
        rspd_s_std = ship_std.rspd;
        rdir_s_std = ship_std.rdir;   
        
        %%% FIX
        % std dev for sin/cos of hed heading, in radians
        shed_s_std = interval_std_var(jd_1min_s, sin(ship_1min.hed*d2r),jd_10bin);
        ched_s_std = interval_std_var(jd_1min_s, cos(ship_1min.hed*d2r),jd_10bin);

        %% 10-min motion-corrected wind statistics: 6x5 array as follows:
        % 1:jd, 2:ustr, 3:vstr, 4:wstr, 5:Tson
        uvwt_10min_avg = interval_avg(uvw_str(:,1), uvw_str(:,2:end),jd_10bin);
        uvwt_10min_std = interval_std(uvw_str(:,1), uvw_str(:,2:end),jd_10bin);
        uvwt_10min_avg(xxx,2:5) = NaN;
        uvwt_10min_std(xxx,2:5) = NaN;
        % define Ubar, the streamwise true windspeed from Edson correction
        Ubar = uvwt_10min_avg(:,2);
        
        % 10-min raw sonic wind averages: 6x5 array as follows:
        % 1:jd, 2:u, 3:v, 4:w, 5:Tson
        son_10min_avg = interval_avg(jd_10Hz, son(:,2:end),jd_10bin);
        son_10min_avg(xxx,2:5) = NaN;
        uson_avg = son_10min_avg(:,2);
        vson_avg = son_10min_avg(:,3);
        wson_avg = son_10min_avg(:,4);
        Tson_avg = son_10min_avg(:,5);

        % 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_10bin);
        uplat_std = plat_10min_std(:,2);
        vplat_std = plat_10min_std(:,3);
        wplat_std = plat_10min_std(:,4);

        % compute psl 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  %%% WTF is this???
        V = son(:,3)/1.15;  % -15% for side-on wind  %%% WTF is this???
        Un = zeros(6,1)*NaN; % empty arrays
        Ue = zeros(6,1)*NaN;
        ii = 1;
        for jj=1:6
            % puts data into earth coordinates N/S, E/W in average sense
            Un(jj) = (nanmean1(U(ii:ii+Npts-1).*cos(heading(ii:ii+Npts-1)) + ...
                    V(ii:ii+Npts-1).*sin(heading(ii:ii+Npts-1)) + sogN_avg(jj)));
            Ue(jj) = (nanmean1(U(ii:ii+Npts-1).*sin(heading(ii:ii+Npts-1)) - ...
                    V(ii:ii+Npts-1).*cos(heading(ii:ii+Npts-1)) + sogE_avg(jj)));
            wdir = atan2(-Ue, -Un)*r2d;
            wdir(wdir<0) = wdir(wdir<0)+360;
            wspd = sqrt(Un.^2 + Ue.^2);
            ii = ii + Npts;
        end
        
        % adjust wind for PISTON 2019 due to suspected flow distortion
        if perc_inc > 90
            wh_u0 = find(wspd < 0.5);
            wh_u1 = find(wspd >= 0.5 & wspd < 1.5);
            wh_u2 = find(wspd >= 1.5 & wspd < 2.5);
            wh_u3 = find(wspd >= 2.5 & wspd < 3.5);
            wh_u4 = find(wspd >= 3.5 & wspd < 4.5);
            wh_u5 = find(wspd >= 4.5 & wspd < 5.5);
            wh_u6 = find(wspd >= 5.5 & wspd < 6.5);
            wh_u7 = find(wspd >= 6.5 & wspd < 7.5);
            wh_u8 = find(wspd >= 7.5 & wspd < 8.5);
            wh_u9 = find(wspd >= 8.5 & wspd < 9.5);
            wh_u10 = find(wspd >= 9.5 & wspd < 10.5);
            wh_u11 = find(wspd >= 10.5 & wspd < 11.5);
            wh_u12 = find(wspd >= 11.5 & wspd < 12.5);
            wh_u13 = find(wspd >= 12.5 & wspd < 13.5);
            wh_u14 = find(wspd >= 13.5 & wspd < 14.5);
            wh_u15 = find(wspd >= 14.5 & wspd < 15.5);
            wh_u16 = find(wspd >= 15.5 & wspd < 16.5);
            wh_u17 = find(wspd >= 16.5);
            u_mult_var = ones(length(wspd),1);
            if perc_inc == 98
                u_mult_var(wh_u0)   = 1.17;
                u_mult_var(wh_u1)   = 1.163;
                u_mult_var(wh_u2)   = 1.156;
                u_mult_var(wh_u3)   = 1.149;
                u_mult_var(wh_u4)   = 1.142;
                u_mult_var(wh_u5)   = 1.135;
                u_mult_var(wh_u6)   = 1.128;
                u_mult_var(wh_u7)   = 1.121;
                u_mult_var(wh_u8)   = 1.114;
                u_mult_var(wh_u9)   = 1.107;
                u_mult_var(wh_u10)  = 1.100;
                u_mult_var(wh_u11)  = 1.093;
                u_mult_var(wh_u12)  = 1.086;
                u_mult_var(wh_u13)  = 1.079;
                u_mult_var(wh_u14)  = 1.072;
                u_mult_var(wh_u15)  = 1.065;
                u_mult_var(wh_u16)  = 1.065;
                u_mult_var(wh_u17)  = 1.065;
            elseif perc_inc == 99
                u_mult_var(wh_u0)   = 1.165;
                u_mult_var(wh_u1)   = 1.160;
                u_mult_var(wh_u2)   = 1.155;
                u_mult_var(wh_u3)   = 1.150;
                u_mult_var(wh_u4)   = 1.145;
                u_mult_var(wh_u5)   = 1.140;
                u_mult_var(wh_u6)   = 1.135;
                u_mult_var(wh_u7)   = 1.130;
                u_mult_var(wh_u8)   = 1.125;
                u_mult_var(wh_u9)   = 1.120;
                u_mult_var(wh_u10)  = 1.115;
                u_mult_var(wh_u11)  = 1.110;
                u_mult_var(wh_u12)  = 1.105;
                u_mult_var(wh_u13)  = 1.100;
                u_mult_var(wh_u14)  = 1.095;
                u_mult_var(wh_u15)  = 1.090;
                u_mult_var(wh_u16)  = 1.085;
                u_mult_var(wh_u17)  = 1.080;
            end
            wspd = u_mult_var.*wspd;
        else
            wspd = u_mult*wspd;
        end
        
        Un(xxx) = NaN; Ue(xxx) = NaN;
        wdir(xxx) = NaN; wspd(xxx) = NaN;

        % 10-min relative wind statistics
        % compute unwrapped rdir in radians
        rdir_10Hz = unwrap(atan2(son(:,3),-son(:,2))); 
        % average and convert to +/- 180 deg format
        rdir = interval_avg_var(jd_10Hz,rdir_10Hz,jd_10bin)*r2d;
        rdir = mod(rdir + 360, 360);
        rdir(rdir>180) = rdir(rdir>180) - 360;
        rdir_std = interval_std_var(jd_10Hz,rdir_10Hz,jd_10bin)*r2d;
        rdir(xxx) = NaN;
        rdir_std(xxx) = NaN;
        % note: this is 10 min streamwise wind speed UNCORRECTED for ship velocity
        rspd = sqrt(uson_avg.^2 + vson_avg.^2 + wson_avg.^2);
        % adjust relative wind for PISTON 2019 due to suspected flow distortion
        if perc_inc > 90
            wh_u0 = find(rspd < 0.5);
            wh_u1 = find(rspd >= 0.5 & rspd < 1.5);
            wh_u2 = find(rspd >= 1.5 & rspd < 2.5);
            wh_u3 = find(rspd >= 2.5 & rspd < 3.5);
            wh_u4 = find(rspd >= 3.5 & rspd < 4.5);
            wh_u5 = find(rspd >= 4.5 & rspd < 5.5);
            wh_u6 = find(rspd >= 5.5 & rspd < 6.5);
            wh_u7 = find(rspd >= 6.5 & rspd < 7.5);
            wh_u8 = find(rspd >= 7.5 & rspd < 8.5);
            wh_u9 = find(rspd >= 8.5 & rspd < 9.5);
            wh_u10 = find(rspd >= 9.5 & rspd < 10.5);
            wh_u11 = find(rspd >= 10.5 & rspd < 11.5);
            wh_u12 = find(rspd >= 11.5 & rspd < 12.5);
            wh_u13 = find(rspd >= 12.5 & rspd < 13.5);
            wh_u14 = find(rspd >= 13.5 & rspd < 14.5);
            wh_u15 = find(rspd >= 14.5 & rspd < 15.5);
            wh_u16 = find(rspd >= 15.5 & rspd < 16.5);
            wh_u17 = find(rspd >= 16.5);
            u_mult_var = ones(length(wspd),1);
            if perc_inc == 98
                u_mult_var(wh_u0)   = 1.17;
                u_mult_var(wh_u1)   = 1.163;
                u_mult_var(wh_u2)   = 1.156;
                u_mult_var(wh_u3)   = 1.149;
                u_mult_var(wh_u4)   = 1.142;
                u_mult_var(wh_u5)   = 1.135;
                u_mult_var(wh_u6)   = 1.128;
                u_mult_var(wh_u7)   = 1.121;
                u_mult_var(wh_u8)   = 1.114;
                u_mult_var(wh_u9)   = 1.107;
                u_mult_var(wh_u10)  = 1.100;
                u_mult_var(wh_u11)  = 1.093;
                u_mult_var(wh_u12)  = 1.086;
                u_mult_var(wh_u13)  = 1.079;
                u_mult_var(wh_u14)  = 1.072;
                u_mult_var(wh_u15)  = 1.065;
                u_mult_var(wh_u16)  = 1.065;
                u_mult_var(wh_u17)  = 1.065;
            elseif perc_inc == 99
                u_mult_var(wh_u0)   = 1.165;
                u_mult_var(wh_u1)   = 1.160;
                u_mult_var(wh_u2)   = 1.155;
                u_mult_var(wh_u3)   = 1.150;
                u_mult_var(wh_u4)   = 1.145;
                u_mult_var(wh_u5)   = 1.140;
                u_mult_var(wh_u6)   = 1.135;
                u_mult_var(wh_u7)   = 1.130;
                u_mult_var(wh_u8)   = 1.125;
                u_mult_var(wh_u9)   = 1.120;
                u_mult_var(wh_u10)  = 1.115;
                u_mult_var(wh_u11)  = 1.110;
                u_mult_var(wh_u12)  = 1.105;
                u_mult_var(wh_u13)  = 1.100;
                u_mult_var(wh_u14)  = 1.095;
                u_mult_var(wh_u15)  = 1.090;
                u_mult_var(wh_u16)  = 1.085;
                u_mult_var(wh_u17)  = 1.080;
            end
            rspd = u_mult_var.*rspd;
        else
            rspd = u_mult*rspd;
        end
        
       %% Wind fix: sonic anemometer is dropping out at times. So fluxes don't get
       % calculated. Fill it with SCS wind for now. 
       wspd_best = wspd;
       wdir_best = wdir;
       rspd_best = rspd;
       rdir_best = rdir;
       wspd_best(xxx == 1) = wspd_s(xxx == 1);
       wdir_best(xxx == 1) = wdir_s(xxx == 1);
       rspd_best(xxx == 1) = rspd_s(xxx == 1);
       rdir_best(xxx == 1) = rdir_s(xxx == 1);

       %% seawater T fix
       % note: use snake T if possible, use ship to fill hole in snake time series if present?
        Ts_best = Tsnk;
        
        %% 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_T(C), 4 Licor_P(kpa), 5 AGC
        lic_jd = lic(:,1);
        lic_10min_avg_all = interval_avg(lic_jd,lic(:,2:end),jd_10bin);
        lic_10min_std_all = interval_std(lic_jd,lic(:,2:end),jd_10bin);
        lic_10min_max_all = interval_max(lic_jd,lic(:,2:end),jd_10bin);
        
        licor_H2O = lic_10min_avg_all(:,2); % Licor_H20, H2O concentration, mmol/m3 
        licor_Tbox = lic_10min_avg_all(:,3); % Licor box T, C
        licor_Pbox = lic_10min_avg_all(:,4); % Licor box P, kPa
        licor_AGC = lic_10min_avg_all(:,5); % Licor AGC, mean diagnostic value
        
        licor_sig_H2O = lic_10min_std_all(:,2); % standard deviation of Licor H2O, mmol/m3
        licor_sig_Pbox = lic_10min_std_all(:,4); % standard deviation of Licor box P, kPa

        %%% make computations and name variables that we'll save
        licor_H2O_MR = (licor_H2O*1e-3).*Rgas_universal.*(licor_Tbox+273.15)./(licor_Pbox*1e3); % mixing ratio, converted mmol/m3 to mol/mol
        licor_rhoa = (licor_Pbox*1e3).*(1-licor_H2O_MR*(1-epsilon))./(Rgas.*(licor_Tbox+273.15)); % air density in licor [kg/m3]
        licor_qa = (licor_H2O*1e-3).*Mw./licor_rhoa; % mmol/m3 to g/kg
        licor_rh = relhum([Ta,licor_qa,slp]);
        
%         % OG
%         lic_h20_10 = lic_10min_avg_all(:,2); % Licor_H20, mmol/m3
%         lic_boxT_10 = lic_10min_avg_all(:,3); % Licor box T, C
%         lic_boxP_10 = lic_10min_avg_all(:,4); % Licor box P, kPa
%         lic_AGC_10 = lic_10min_avg_all(:,5); % Licor AGC, mean diagnostic value
%         
%         %%% make computations and name variables that we'll save
%         licor_H2O_a = (lic_h20_10*1e-3).*Rgas_universal.*(lic_boxT_10+273.15)./(lic_boxP_10*1e3); % mmol/m3 to mol/mol
%         rhoa_licor = (lic_boxP_10*1e3).*(1-licor_H2O_a*(1-epsilon))./(Rgas.*(lic_boxT_10+273.15)); % air density in licor [kg/m3]
%         qa_licor = (lic_h20_10*1e-3).*Mw./rhoa_licor; % mmol/m3 to g/kg
%         rh_licor = relhum([ta_psd,qa_licor,slp_psd]);
 
        
        %%% quality control licor data based on AGC. Good value where AGC < 60
        % note: the high frequency licor data are used in computing eddy
        % covariance and inertial dissipation latent heat flux, not these
        % 10-min values computed and saved here. 
        wh_good_licor = find(licor_AGC < 60);
        wh_bad_licor = find(licor_AGC >= 60);
        licor_rhoa(wh_bad_licor) = nan; % kg / m3
        licor_qa(wh_bad_licor) = nan; % g / kg
        licor_rh(wh_bad_licor) = nan; % rh
        licor_H2O(wh_bad_licor) = nan;  % mmol/m3
        licor_H2O_MR(wh_bad_licor) = nan;  % mol/mol
        licor_Tbox(wh_bad_licor) = nan; % C
        licor_Pbox(wh_bad_licor) = nan; % kPa
        
        % 10-min specific humidities (q = g/kg)
        qa      = qair_p(Ta,rh,slp - (0.125*zq));
        qa_s    = qair_p(Ta_s,rh_s,slp_s - (0.125*zq_ship));
        qsnk    = qsea_p(Tsnk,slp);
        q_s     = qsea_p(Ttsg_s,slp_s);

        % 10-min mean air density (kg/m3) and H2O number density (mmol/m3) from bulk met
        rhoa = air_density(Ta,slp,rh);
        % number density from vaisala RH, mmol/m3
        H2O_10min_avg = qa.*rhoa/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(jj);
            ii = ii + Npts;
        end

        % decorrelate licor water vapor vars w/respect to 3-axis motion
        if decorr
%             q_10Hz = decorr_q(q_10Hz_raw,uvwplat,mot(:,2:4));
            jd_bins = lic(1,1):10/1440:son(1,1)+60/1440; jd_bins = jd_bins';
            X = [ones(length(q_10Hz_raw),1),accplat,uvwplat];
            [q_10Hz,mu_q] = interval_decorr(lic_jd,q_10Hz_raw,jd_bins,X);
            [H2O_10Hz,~] = interval_decorr(lic_jd,lic(:,2),jd_bins,X);
            mu_q = [jd_bins(1:end-1),mu_q];
        else
            q_10Hz = q_10Hz_raw;
            H2O_10Hz = lic(:,2);
        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_10bin);
        wv_cov_10min = interval_cov(uvw_str(:,1),uvw_str(:,4),uvw_str(:,3),jd_10bin);
        wT_cov_10min = interval_cov(uvw_str(:,1),uvw_str(:,4),T_10Hz,jd_10bin); % This is for sonic T
        wH2O_cov_10min = interval_cov(uvw_str(:,1),uvw_str(:,4),H2O_10Hz,jd_10bin);
        wq_cov_10min = interval_cov(uvw_str(:,1),uvw_str(:,4),q_10Hz,jd_10bin);
        wu_cov_10min(xxx,2) = NaN;
        wv_cov_10min(xxx,2) = NaN;
        wT_cov_10min(xxx,2) = NaN;
        wH2O_cov_10min(xxx,2) = NaN;
        wq_cov_10min(xxx,2) = NaN;

        wu_cov = wu_cov_10min(:,2);
        wv_cov = wv_cov_10min(:,2);
        wq_cov = wq_cov_10min(:,2)*sens_sep_factor;
        wH2O_cov = wH2O_cov_10min(:,2);
        wT_cov = wT_cov_10min(:,2);
        
        % Simon's alternate computation of wT and wq, Webb corrected
        [wT_sds, wq_sds] = wpl_schotanus(Ta+273.15,qa/1000,rhoa,wT_cov,wH2O_cov);
        wq_sds = wq_sds*sens_sep_factor*1000; % kg/kg -> g/kg, and with correction for sensor sep
        
        %% ---------INITIAL COARE BULK MODEL RUNS WITH SHIP AND PSL DATA--------
        %{
        input format: A = coare35vn(u,zu,t,zt,rh,zq,P,Tsnk,Rs,Rl,lat,zi,rain,cp,sigH)
        Rs and Rl are downwelling and positive. 
        Rs_down = positive, abs value - dont direction dependent
        Rs_up = positive, abs value - dont direction dependent
        Rl_down = positive, abs value - dont direction dependent
        Rl_up = positive, abs value - dont direction dependent
        anytime a radiometer is pointed at sky, it gives a positive number.
        If you pointed it at ground (measuring reflection), it would give a positive number. 
        Rnet = down  - up, both up and down are positive, so that positive net values go into the surface and
        warm it. Net IR tends to be negative (cooling the ocean, warming
        atm). Net solar tends to be positive 
        

        cp, sigH set to NaN
        zi=600 and jcool=1 (default)
        FOR MISOBOB only: using ship sst for all, psl snake not working
        no rain data

        A = [usr tau hsb hlb hlwebb tsr qsr zot zoq Cd Ch Ce L zet dter dqer tkt RF Cdn_10 Chn_10 Cen_10];
        %     1   2   3   4     5    6   7   8   9  10 11 12 13 14  15   16   17 18   19     20     21

        usr (1), tau (2), hsb (3), hlb (4), L (13), zet (14), dter (15), tkt (17),
        RH (18), Cdn_10 (19), Chn_10 (20), Cen_10 (21)
        
        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
        %}

        % COARE 3.5 - using psl rain rate, Tair, RH, PIR, SLP, Tsnk for both runs
  
       %%% this version of code has output array of 26 columns.       
       coare_psl = coare35vn_motcorr(wspd_best,zu,Ta,zt,rh,zq,slp,Ts_best,...
                     rs,rl,lat_s,600,rain,NaN,NaN,10);
       coare_s = coare35vn_motcorr(wspd_s,zu_ship,Ta_s,zt,rh_s,zq,slp_s,Ttsg_s,...
                     rs_s,rl_s,lat_s,600,rain,NaN,NaN,10);

       %%% how I want to name them or explain them to others
%        cfields_named =      {'ustar'; 'tau'; 'hs'; 'hl';'Tstar'; 'qstar';'L';'zet';'dT_skin';'dz_skin';'hr';'Cd';       'Ch';  'Ce';    'U10';   'evap';'T10';'q10';  'rh10'};
       cfields = {'usr';'tau';'hsb';'hlb';'hlwebb';'tsr';'qsr';'zot';'zoq';...  %1-9
           'Cd';'Ch';'Ce';'L';'zet';'dter';'dqer';'tkt';'RF';'Cdn_10';...       %10-19
           'Chn_10';'Cen_10';'U10';'Evap';'T10';'Q10';'RH10'};                  %20-26         
       
        for i = 1:length(cfields)
            eval([cfields{i} ' = coare_psl(:,i);']);
%             cfields(i) =
            eval([cfields{i} '_s = coare_s(:,i);']);
        end

        
        % IF USING OLDER COARE 2.5 CODE...
        % input = [U odec Tsnk 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_psl(ii,:) = coul_lkb([wspd(ii),0,Tsnk(ii),Ta(ii),qsnk(ii),qa(ii),...
%                                   rs(ii),rl(ii),rain_all(ii),600,slp(ii),zu,zt,zq]);
%             coare_s(ii,:) = coul_lkb([wspd_s(ii),0,Ttsg_s(ii),Ta_s(ii),qs_s(ii),qa_s(ii),...
%                                   rs_s(ii),rl_s(ii),rain_all(ii),600,slp_s(ii),zu_ship,zt_ship,zq_ship]);
%         end;


        %% compute covariance heat fluxes
        % correct sonic wT for humidity -  use C35 for <wq>=-usb*qsb
        % <wT> = <wTs> - 0.51*T*<wq>
        wT_cov = wT_cov + 0.51*(Ta+273.15).*usr.*qsr*1e-3;
        %%% NOTE - correct Simon's too? 
        
        % 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 = wT_cov.* rhoa * 1004.67;
        hs_cov_sds = wT_sds.* rhoa * 1004.67;
        
        % latent heat vap [J/g] for water and seawater from ship SST and ship salinity
        [Le_w_10min, Le_sw_10min] = Le_water(Ttsg_s, Stsg_s);
        % hl = wq * rhoair * Le_w = g/kg m/s * kg/m3 * J/g = W/m2
        hl_cov = wq_cov .* rhoa .* Le_sw_10min;
        % Webb corrected hl
        hl_covW = hl_cov + hlwebb; % using C35 webb correction
        hl_cov_sds = wq_sds .* rhoa .* Le_sw_10min; % Simon's method, Le in J/kg
         
        %% --------------------COMPUTE SPECTRA AND COSPECTRA--------------------
        % initialize arrays for hourly spectral data
        Fx = NaN(1,48);      % empty arrays for 6x48-bin smoothed spectra
        Pu = NaN(6,48);
        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 = NaN(1,1+Npts/2);  % empty arrays for raw spectra
        Puu = NaN(6,1+Npts/2);  % used for idiss calculations...
        Pvv = Puu; Pww = Puu; Ptt = Puu; Ptt2 = Puu; Pqq = Puu;
        Tson_noise = NaN(6,1);

        % compute variance spectra and cospectra in 10-min windows
        % but only for periods when sonic data is good
        segs = 1:6;
        for jj = segs(~xxx)  % just good 10-min intervals
            ii = 1 + Npts*(jj-1);
            [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 Tson...
            [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 Tson spectra (for idiss)
            % assumes sonic spectrum at ~5 Hz is noise (from last 3 points in smoothed spectrum)
            Tson_noise(jj) = median(Pt(jj,end-2:end));
            Ptt2(jj,:) = Ptt(jj,:) - Tson_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 rspd
        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, rspd, 0);
        Cw2a = idiss_struc_func(Pw, Fx, fminU, fmaxU, rspd, 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, rspd, 0);
%         Cw2a = idiss_struc_func(Pww, Fxx, fminU, fmaxU, rspd, 0);

        % Ct2 - this is Tson...
        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, rspd, 0);

%         fminT = find(Fxx>0.28,1,'first');
%         fmaxT = find(Fxx<0.4,1,'last');
%         Ct2a = idiss_struc_func(Ptt2, Fxx, fminT, fmaxT, rspd, 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, rspd, 0);

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


        % apply Wyngaard & Clifford correction
        % f1, f2, f3 are same length as uvar, etc.
        ustr2 = rspd.^2;
        uvar = uvwt_10min_std(:,2).^2;
        vvar = uvwt_10min_std(:,3).^2;
        wvar = uvwt_10min_std(:,4).^2;
        Tvar = uvwt_10min_std(:,5).^2;
        f1 = 1 - (uvar/9.0 - 2*vvar/3.0 - 2*wvar/3.0)./ustr2;
        ff2 = 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./ff2;
        Ct2a = Ct2a./f3;
        Cq2a = Cq2a./f3;

        % compute usid, tsid, qsid, Lid
        % qsid output is kg/kg
        [usid_a,tsid_a,qsid_a,Lid_a] = idiss_sst(Cu2a,Ct2a,Cq2a,qa,Ta,Tsnk,zu,zq,L,usr,tsr,qsr,dter);

        %% 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, rspd);
        Cw2b = idiss_spec_fit(Pw, Fx, fminU, fmaxU, rspd);

        % 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, rspd);

        % 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, rspd);

        % apply Wyngaard & Clifford correction
        Cu2b = Cu2b./f1;
        Cw2b = Cw2b./ff2;
        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_sst(Cu2b,Ct2b,Cq2b,qa,Ta,Tsnk,zu,zq,L,usr,tsr,qsr,dter);

        %% save 10_Hz wind and motion data
        % wrap rdir_10Hz back to original format
        rdir_10Hz_wrap = mod(rdir_10Hz*r2d + 360, 360); 
        rdir_10Hz_wrap(rdir_10Hz>180) = rdir_10Hz_wrap(rdir_10Hz>180) - 360;
        uvwmat = [uvw_str,rdir_10Hz_wrap,q_10Hz]; %[jd_10Hz,u,v,w,T,rdir,q_10Hz]
        fpath = fullfile(path_uvwStreamData,['UVWstream_',jd,'_',hr,'_',dcname,'.txt']);
        save (fpath,'uvwmat','-ascii');

        motmatall=[jd_10Hz,uvwplat,xyzplat,euler];
        fpath = fullfile(path_motionData,['motion_',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

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

                wq = nanmean1(Cqw,1); wq_raw = nanmean1(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,'wq_cospectrum',['wqCospec_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                wt = nanmean1(Ctw,1); wT_raw = nanmean1(Ctw_raw,1);
                figure(6); clf; semilogx(Fx,Fx.*wt,'b-',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,'wt_cospectrum',['wtCospec_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                % spectra plots
                Sw = nanmean1(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','northeast');
                title([jd,' ',hr,' Spectrum of vertical wind velocity, W']); grid;
                ppath = fullfile(path_spectraPlots,'w_spectrum',['wSpectrum_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

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

                Sq = nanmean1(Pq,1);
                figure(9);clf;loglog(Fx,Fx.*Sq,'b-',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,'q_spectrum',['qSpectrum_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                Su = nanmean1(Pu,1);
                Sv = nanmean1(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','northeast');
                title([jd,' ',hr,' U V Wind Spectra']); grid;
                ppath = fullfile(path_spectraPlots,'uv_spectrum',['uvSpectrum_',jd,'_',hr,'_',dcname,graphformat]);
                print(graphdevice,ppath);

                close all;
            end
        end
        
        %% save each hour of 10-min data in daily mat and text files
        
        % coare_psl variables:         
        save_coare_fields = [1          2     3     4      6       7     13    14     15    17:26];
        %%% how they are labeled in COARE algorithm
%         my_cfields =            {'usr'; 'tau'; 'hsb'; 'hlb'; 'tstr'; 'qstr'; 'L'; 'zet'; 'dter'; 'tkt'; 'RF'; 'Cdn_10'; 'Chn_10'; 'Cen_10'; 'U10N'; 'Evap'; 'T10'; 'Q10'; 'RH10'};
        %%% how I want to name them or explain them to others
%         cfields_named =      {'ustar'; 'tau'; 'hs'; 'hl';'Tstar'; 'qstar';'L';'zet';'dT_skin';'dz_skin';'hr';'Cd';       'Ch';  'Ce';    'U10';   'evap';'T10';'q10';  'rh10'};

       cfields_named = {'ustar';'tau';'hs';'hl';'hlwebb';'tstar';'qstar';'zot';'zoq';...    %1-9
       'Cd';'Ch';'Ce';'L';'zet';'dter';'dqer';'dT_skin';'hr';'Cdn_10';...                   %10-19
       'Chn_10';'Cen_10';'U10';'Evap';'T10';'q10';'rh10'};                                  %20-26

        
        % coare_s variables:         
        save_coare_s_fields = [             2     3     4   ];
        %%% how they are labeled in COARE algorithm
%         cfields_s =            {      'tau'; 'hsb'; 'hlb'};
        %%% how I want to name them or explain them to others
%         cfields_s_named =      {     'tau_s'; 'hs_s'; 'hl_s'};
        
        fields_to_save = {'rain'; 'Ta'; 'rh'; 'Ts_best'; 'rs_dn'; 'rl_dn'; 'slp'; 'P'; 'qa'; 'rhoa';...
        'cog_avg'; 'sog_avg'; 'hed_avg'; 'lat_s'; 'lon_s'; 'Un'; 'Ue'; 'wspd_best'; 'wdir_best'; 'rspd_best'; 'rdir_best';...
        'Ttsg_s'; 'Stsg_s';'Tskin';'uplat_std';'vplat_std';'wplat_std';'the_tilt';...
        'wu_cov';'wv_cov';'wq_cov';'wT_cov';'wq_sds';'wT_sds';...
        'licor_qa';'licor_rh';'licor_rhoa';'licor_AGC';'licor_H2O_MR';'licor_H2O';...
            'licor_Tbox';'licor_Pbox';'licor_sig_Pbox';'licor_sig_H2O';'missing';'badPnts';...
        'usid_a';'usid_b';'tsid_a';'qsid_a';'hlwebb';'hs_cov';'hs_cov_sds';'hl_cov';'hl_covW';'hl_cov_sds';...
        'rdir_s';'rspd_s';'wdir_s';'wspd_s';'rl_s_dn';'rs_s_dn';'Ta_s';'qa_s';'rh_s';'sog_s';'cog_s';...
        'Ttsg1'; 'Ttsg2'; 'Stsg1'; 'Stsg2'; 'Ftsg1'; 'Ftsg2'; 'ID_Ftsg1'; 'ID_Ftsg2'; 'ID_tsg';...
        'spd_log_u'; 'spd_log_v';...
        'prate_s'; 'paccum_s';...
        'uvar';'vvar';'wvar';'Tvar'};
    
    % other stuff to save needed in dared
        cu   =   (Cu2a+Cu2b)/2;
        cw   =   (Cw2a+Cw2b)/2;
        ct   =   (Ct2a+Ct2b)/2;
        cq  =    (Cq2a+Cq2b)/2;
        usi  =   (usid_a+usid_b)/2;
        tsi  =   (tsid_a+tsid_b)/2;
        qsi  =   (qsid_a+qsid_b)/2;
    Ts_noise
    sig_h   =   da(:,113);   % std dev of heading, >5 indicates ship maneuver
sig_sog =   da(:,98);    % std dev of ship speed, m/s (>.2 indicates maneuver)
if any(isfinite(odec))
    % compute surface current from odec, head, sog_ship, cog_ship
    % current direction is 'to' - oceanographic convention
    [cspd,cdir] = current(odec,hed_best,sog_best,cog_best);
% first square all standard deviations
da(:,33) = da(:,33).^2;           % std wspd, ship
da(:,37) = da(:,37).^2;           % std wdir, ship
da(:,95:96) = da(:,95:96).^2;     % std Rl & Rs
da(:,98) = da(:,98).^2;           % std SOG, ship
da(:,102) = da(:,102).^2;         % std COG, ship
da(:,106) = da(:,106).^2;         % std COG, PSD
da(:,113) = da(:,113).^2;         % std heading, ship
da(:,129:133) = da(:,129:133).^2; % std, sp1, sp2, sp3, sin/cos heading
da(:,154) = da(:,154).^2;         % std odec
da(:,165) = da(:,165).^2;         % std rwdir, PSD
cr(:,12:15) = cr(:,12:15).^2;     % std licor vars
% for turbulence variables, insert filtered means from y2 into y1
y1(:,116:121) = y2(:,116:121);    % covariances, PSD
y1(:,134:149) = y2(:,134:149);    % ID variables, PSD
y1(:,160:162) = y2(:,160:162);    % hs/hl, PSD
y1(:,181:184) = y2(:,181:184);    % hs/hl, Simon's function
y1(:,44) = y2(:,44);              % tilt
y3(:,8:9) = ycr(:,8:9);           % wq and hlc
[~,ywspdx,ywdirx] = timave_vec1(jdy,dta,dax(:,16),dax(:,41));
[~,ywspd_ship,ywdir_ship] = timave_vec1(jdy,dta,da(:,30),da(:,34));

        % variables in this program:
        % Simplify labels and make consistent with met/sea file
        % note: lat, lon, hed are vars that exist prior to this
        fields_to_name = {'prate'; 'Ta'; 'rh'; 'Tsnk'; 'sw_dn'; 'lw_dn'; 'slp'; 'P'; 'qa'; 'rhoa';...
        'cog'; 'sog'; 'hed'; 'lat'; 'lon'; 'Un'; 'Ue'; 'wspd'; 'wdir'; 'rspd'; 'rdir';...
        'Ttsg'; 'Stsg';'Tskin';'uplat_std';'vplat_std';'wplat_std';'tilt';...
        'wu_cov';'wv_cov';'wq_cov';'wT_cov';'wq_sds';'wT_sds';...
        'licor_qa';'licor_rh';'licor_rhoa';'licor_AGC';'licor_H2O_MR';'licor_H2O';...
            'licor_Tbox';'licor_Pbox';'licor_sig_Pbox';'licor_sig_H2O';'missing';'badPnts';...
        'usid_a';'usid_b';'tsid_a';'qsid_a';'hlwebb';'hs_cov';'hs_cov_sds';'hl_cov';'hl_covW';'hl_cov_sds';...
        'rdir_s';'rspd_s';'wdir_s';'wspd_s';'lw_dn_s';'sw_dn_s';'Ta_s';'qa_s';'rh_s';'sog_s';'cog_s';...
        'Ttsg1_s'; 'Ttsg2_s'; 'Stsg1_s'; 'Stsg2_s'; 'Ftsg1_s'; 'Ftsg2_s'; 'ID_Ftsg1_s'; 'ID_Ftsg2_s'; 'ID_tsg_s';...
        'spd_log_u_s'; 'spd_log_v_s';...
        'prate_s'; 'paccum_s';...
        'uvar';'vvar';'wvar';'Tvar'};
 
        % address sign or orientation issues with these arrays for saving
        the_tilt = tilt';

        %%% radiation is measured as positive value when downwelling.
        rl_s_dn = rl_s;
        rs_s_dn = rs_s;
        rs_dn = rs;
        rl_dn = rl;
        
        %%% start a structure
        if hhh == 0
            clear ff;
            
            ff.jd = jd_10min;
            ff.t = datenum(str2double(yyyy),0,0,0,0,0) + jd_10min;
            ff.sst = Ts_best - dter;

            for i = save_coare_fields
%                 ff.(cfields_named{i}) = coare_psl(:,save_coare_fields(i));
                eval(['ff.' cfields_named{i} ' = ' cfields{i}]);
            end

            for i = save_coare_s_fields
%                 ff.(fields_s_named{i}) = coare_s(:,save_coare_s_fields(i));
                eval(['ff.' cfields_named{i} '_s = ' cfields{i} '_s;']);
            end
            
            for i = 1:length(fields_to_save)
                % add to structure as newly named var
                eval(['ff.(fields_to_name{i}) = ' fields_to_save{i} ';']);
%                 ff.(fields_to_name{i}) = fields_to_save{i};
            end
            
        %%% add to the structure
        else
            ff.jd = [ff.jd; jd_10min];
            ff.t = [ff.t; datenum(str2double(yyyy),0,0,0,0,0) + jd_10min];
            ff.sst = [ff.sst; Ts_best - dter];
            
            for i = save_coare_fields
%                 ff.(cfields_named{i}) = [ff.(cfields_named{i}); cfields_named{i}];
                eval(['ff.' cfields_named{i} ' = [ff.' cfields_named{i} ' ; ' cfields{i} '];']);
            end
            
            for i = save_coare_s_fields
%                 ff.(cfields_s_named{i}) = [ff.(cfields_s_named{i}); ];
                eval(['ff.' cfields_named{i} '_s = [ff.' cfields_named{i} '_s; ' cfields{i} '_s];']);
            end
            
            for i = 1:length(fields_to_save)
                % add to structure as newly named var
                eval(['ff.(fields_to_name{i}) = [ff.(fields_to_name{i}); ' fields_to_save{i} '];']);
%                 ff.(fields_to_name{i}) = [ff.(fields_to_name{i}); fields_to_save{i}];
            end
            
        end
        
        
    %% format da, sp, cr output and concatenate to daily array
%{
    da 10-min variable designations as follows:  144 x 176 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    $ship           COARE 3.5 input, ship:[Utrue,Tair,RH,P_mb,Tsnk,
                            Rs,Rl,lat,rain]
    11:15   no_data         To maintain historical column numbering...
    16:24   $psl            COARE 3.5 input, PSL, 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      rspd            streamwise relative wind speed (m/s)from PSL sonic
    39      Un              N component, corrected true windspeed, PSL sonic
    40      Ue              E component, corrected true windspeed, PSL sonic
    41      wdir            true wind direction from N,E components, PSL sonic
    42      rdir            relative wind direction from PSL sonic
    43      azm             azimuth angle from first rotation, radians
    44      tilt            tilt angle from first rotation, radians
    45      Tsnk            seasnake SST, degC
    46      Tair            air temperature, PSL vaisala, degC
    47      RH              RH, PSL vaisala
    48      RH_ship         RH, ship sensor
    49      Q               specific humidity, PSL sensors, g/kg
    50      Q_ship          specific humidity, ship sensors, g/kg
    51      Qsea            sea surface specific humidity from seasnake SST
    52      rhoa            moist air density from PSL sensors (kg/m3)
    53:73   $ship           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   $psl            COARE 3.5 output, PSL data, as above
    95      Rl_std          std PIR flux, PSL
    96      Rs_std          std PSP flux, PSL
    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             mean cog, deg, PSL gps
    106     cog_std         std cog, PSL gps
    107     cog_min         minimum cog, PSL gps
    108     cog_max         maximum cog, PSL gps
    109     sog             mean sog, m/s, PSL gps
    110     sog_min         minimum sog, PSL gps
    111     sog_max         maximum sog, PSL gps
    112     hed             mean heading, deg, ship gyro
    113     hed_std         std heading, ship gyro
    114     hed_min         minimum heading, ship gyro
    115     hed_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 hed heading from ship (NOTE: correct?)
    133     hd_cos_std      stdev cosine hed heading (NOTE: correct?)
    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, ship gps
    151     lat             latitiude, ship 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     slp_mb          sea level atmospheric pressure, PSL sensor, mb
    158     missing         missing sonic wind data lines per 10 min interval
    159     badPnts         # NaNs in sonic data per 10 min inteerval
    160     hs_cov          sensible heat flux from covariance
    161     hl_cov          latent heat flux from covariance, uncorrected
    162     hl_covW         webb corrected hl_cov
    163     Stsg_s          salinity, ship system, psu
    164     Tson_noise      Sonic temperature noise variance
    165     rdir_std        std dev relative wind direction
    166     rdir_wxt        rel wind direction ,0-360 deg from bow, WXT
    167     rspd_wxt        wind speed, m/s, WXT
    168     Ta_wxt          air temperature, C, WXT
    169     rh_wxt          RH, %, WXT
    170     slp_wxt         sea level air pressure, mb or hPa, WXT
    171     prate_wxt       rain rate, mm/hr, WXT
    172     wspd_wxt        true wind speed, m/s, WXT
    173     wdir_wxt        true wind dir, m/s, WXT
    174     run_wxt          Un rel, m/s, WXT + south-to-north like the met convention "v"
    175     ruw_wxt          Uw rel, m/s, WXT + east-to-west like the met convention "-u"
    176     rosr_s          ship ROSR radiometric SST
    177     rspd_s          ship rel wind speed
    178     rdir_s          ship rel wind dir, +/- 180 deg
    179     rspd_s_std      ship rspd std dev
    180     rdir_s_std      ship rdir std dev
    181     wT_sds          heat flux covariance from Simon's function
    182     wq_sds          latent flux covariance from Simon's function, with CWF instrument separation correction
    183     hs_cov_sds      simon's covariance hs
    184     hl_cov_sds      simon's covariance hl
%}
        no_data = zeros(6,1)*NaN;
        StRow = hhh*6+1;
        EndRow = StRow + 5;
        da(StRow:EndRow,:) = [jd_10min,wspd_s,Ta_s,rh_s,slp_s,...        % 1-5
            Ttsg_s,rs_s,rl_s,lat_s,rain_s,no_data,wspd_best,...          % 6-12
            wdir_best,rspd_best,rdir_best,wspd,Ta,rh,slp,...   % 13-19
            Ts_best,rs,rl,lat_s,rain,no_data,no_data,...           % 20-26
            no_data,no_data,no_data,wspd_s,ship_min.S,...             % 27-31
            ship_max.S,ship_std.S,wdir_s,...                  % 32-34
            ship_min.rh,ship_max.rh,ship_std.rh,...        % 35-37
            rspd,Un,Ue,wdir,rdir,azm',tilt',...       % 38-44
            Ts_best,Ta,rh,rh_s,qa,qa_s,qsnk,...             % 45-51
            rhoa,...                                                        % 52
            usr_s, tau_s, hsb_s, hlb_s, hlwebb_s, tsr_s, qsr_s, zot_s, zoq_s, Cd_s, ...  % 53-62
            Ch_s, Ce_s, L_s, zet_s, dter_s, dqer_s, tkt_s, RF_s, Cdn_10_s, Chn_10_s, Cen_10_s,...  % 53-73
            usr, tau, hsb, hlb, hlwebb, tsr, qsr, zot, zoq, Cd,... % 74-83
            Ch, Ce, L, zet, dter, dqer, tkt, RF, Cdn_10, Chn_10, Cen_10,... % 74-94
            met_std.pir1,...                                              % 95
            met_std.psp1,sog_s,ship_std.sog,...                   % 96-98
            ship_min.sog,ship_max.sog,cog_s,...                   % 99-101
            ship_std.cog,ship_min.cog,ship_max.cog,...        % 102-104
            cog_avg,cog_std,cog_min,...        % 105-107
            cog_max,sog_avg,sog_min,...        % 108-110
            sog_max,hed_s,ship_std.hed,...                  % 111-113
            ship_min.hed,ship_max.hed,no_data,...                 % 114-116
            wu_cov,no_data,wv_cov,no_data,...                                   % 117-120
            wT_cov,no_data,no_data,no_data,uvar,...                         % 121-125
            vvar,wvar,...                                                    % 126-127
            Tvar,uplat_std,vplat_std,wplat_std,...                                      % 128-131
            shed_s_std,ched_s_std,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_s,lat_s,Ubar,...          % 147-152
            no_data,no_data,no_data,no_data,slp,missing,badPnts,...          % 153-159
            hs_cov,hl_cov,hl_covW,Stsg_s,Tson_noise,...                       % 160-164
            rdir_std, met_avg.rdir_wxt, met_avg.rspd_wxt, met_avg.Ta_wxt, met_avg.rh_wxt, met_avg.slp_wxt,...         % 165-170
            met_avg.prate_wxt, met_avg.wspd_wxt, met_avg.wdir_wxt, met_avg.run_wxt, met_avg.ruw_wxt, rosr_s,...             % 171-176
            rspd_s,rdir_s,rspd_s_std,rdir_s_std,...                     % 177-180
            wT_sds, wq_sds,hs_cov_sds,hl_cov_sds];       % 181-184

            %%% this chunk was broken out by variable for tracking individual fields
            % rhoa,coare_s(:,1:21),coare_psl(:,1:21),met_10min_std(:,6),... % 52-95 

% output of array A from coare35vn_motcorr:
% coare_psl = [usr tau hsb hlb hlwebb tsr qsr zot zoq  Cd Ch Ce L zet dter dqer tkt RF Cdn_10 Chn_10 Cen_10 U10 Evap T10 Q10 RH10];
%    coare_psl:  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]
%   da for psl: 74  75  76  77    78   79  80  81  82  83 84 85 86 87  88   89   90 91   92     93     94  ... the rest is not saved in da]
%  da for ship: 53  54  55  56    57   58  59  60  61  62 63 64 65 66  67   68   69 70   71     72     73  ... the rest is not saved in da]

%     %%% calculate inertial dissipation fluxes for plotting
%     da(:,52) = rhoa; % air density PSL
%     da(:,144) = tsid_a; % sensible heat flux from inertial dissipation method A
%     da(:,142) = usid_a; % stress from inertial dissipation method A
%     da(:,20) = Ts_best; % best surface temp
%     da(:,163) = Stsg_s; % salinity from ship
%     da(:,146) = qsid_a; % latent heat flux from inertial dissipation method A
%     
%     da(:,78) = coare_psl(:,5) = hlwebb; % webb latent flux from coare_psl
%     da(:,56) = coare_s(:,4) = hlb_s; % bulk latent flux from coare_s
%     da(:,77) = coare_psl(:,4) = hlb; % bulk latent flux from coare_psl
%     da(:,55) = coare_s(:,3) = hsb_s; % bulk sensible flux from coare_s
%     da(:,76) = coare_psl(:,3) = hsb; % bulk sensbile flux from coare_psl
%         
%     hlwebb = coare_psl(:,5);
%     hlb_s = coare_s(:,4);
%     hlb = coare_psl(:,4);
%     hsb_s = coare_s(:,3);
%     hsb = coare_psl(:,3);
%     
    
%{
        10-min output array for CO2/H2O statistics: 144x17
        follows historical format, but CO2 covariance is not computed.

        cr column designations as follows:  'good' variables are NaN
            for this project.  Filtering is done in the da routine.

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

        cr(StRow:EndRow,:) = [jd_10min,...
            licor_H2O,... 
                    no_data,...
                    no_data,...
                    no_data,...
                licor_Pbox,...
                licor_Tbox,...
                wq_cov,...
                hl_cov,...
                    no_data,...
                    no_data,...
                licor_sig_H2O,...
                    no_data,...
                    no_data,...
                    no_data,...
                licor_sig_Pbox,...
                licor_AGC];

%{
        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];

        % concatenate lag times and 10-min decorr coefficients
        lag(hhh+1,:) = lagmat;
        mu_decorr(StRow:EndRow,:) = mu;
        mu_q_decorr(StRow:EndRow,:) = mu_q;

    end % end of 24 hour loop

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    %% save daily 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_',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_',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_',jd,'_',dcname,'.txt'];
    fl = fullfile(path_sp,sp_fname);
    save (fl,'sp','-ascii','-tabs');

    lag_fname = ['lag_',jd,'_',dcname,'.txt'];
    fl = fullfile(path_lag,lag_fname);
    save (fl,'lag','-ascii','-tabs');

    mu_fname = ['mu_',jd,'_',dcname,'.txt'];
    fl = fullfile(path_mu,'w_Ts',mu_fname);
    save (fl,'mu_decorr','-ascii','-tabs');

    mu_q_fname = ['mu_q_',jd,'_',dcname,'.txt'];
    fl = fullfile(path_mu,'q',mu_q_fname);
    save (fl,'mu_q_decorr','-ascii','-tabs');

    %% Prep for daily plots

    %%% Compute good data flags for calculating inertial dissipation and eddy
    %%% covariance fluxes of sensible heating, evaporation, and stress

    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(:,42)<90 & da(:,42)>-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) & da(:,158)<100 & ...
               da(:,159)<100 & da(:,24)<5);
    
    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);
            
    % wh_good_motion = kk 
    %      general filter for rdir, std heading, std sog, vplat variance,
    %     (sqrt(std_w)/ugw < .7*.0015*U^2) & (sqrt(std_v)/ugu < 1.4),
    %     BadPnts < 100, Missing < 100 and rain<5.
    good_motion = zeros(6,1);
    good_motion(kk) = 1;
    
    % wh_good_motion_licor = hh 
    %     filter for kk & licor AGC < 60 
    good_motion_licor = zeros(6,1);
    good_motion_licor(hh) = 1;
    
    % wh_good_motion_id = jj
    %     additional filter for kk & isfinite(usid) & (abs(-wu-usid)^2<.02+.006*U)
    good_motion_id = zeros(6,1);
    good_motion_id(jj) = 1;
    
    % wh_good_motion_id_licor = mm;     
    %     filter for both hh (motion and licor) and jj (motion and id) for hl_id
    good_motion_id_licor = zeros(6,1);
    good_motion_id_licor(mm) = 1;
   

    %%% calculate inertial dissipation (id) and eddy covariance (cov) fluxes for saving and plotting
    
%     to get here, must save following through hourly loops: 
%     {'usid_a';'usid_b';'tsid_a';'qsid_a';'hlwebb';'hs_cov';'hs_cov_sds';'hl_cov';'hl_cov_sds'};
    
    % idiss method a for latent and sensible heat fluxes
    hsid = -ff.rhoa*1004.67.*ff.tsid_a.*ff.usid_a; 
    [Le_w, ~] = Le_water(ff.Tsnk, ff.Stsg);
    hlid = -ff.rhoa.*Le_w.*ff.qsid_a*1000.*ff.usid_a - ff.hlwebb; % Webb corrected, g/kg
%     maxhl = max(max(hlb_s),max(hlb))+20; minhl = min(min(hlb_s),min(hlb))-20;
%     maxhs = max(max(hsb_s),max(hsb))+5; minhs = min(min(hsb_s),min(hsb))-5;
    ff.hs_id = nan(length(ff.t),1);
    ff.hl_id = nan(length(ff.t),1);
    ff.hs_id(kk) = hsid(kk);
    ff.hl_id(mm) = hlid(mm);
    
    % id a and b methods for stress
    ff.tau_ida = nan(length(ff.t),1);
    ff.tau_idb = nan(length(ff.t),1);
    ff.tau_ida(jj) = (ff.rhoa(jj).*ff.usid_a(jj).^2);
    ff.tau_idb(jj) = (ff.rhoa(jj).*ff.usid_b(jj).^2);
    
    % eddy covariance method for stress
    ff.tau_cov = nan(length(ff.t),1);
    ff.tau_cov(kk) = -(ff.rhoa(kk).*ff.wu_cov(kk));

    % eddy covariance method for sensible and latent heat
    ff.hs_covW = nan(length(ff.t),1);
    ff.hs_covS = nan(length(ff.t),1);
    ff.hs_covW(kk) = ff.hs_cov(kk);
    ff.hs_covS(kk) = ff.hs_cov_sds(kk);
    
    ff.hl_covW = nan(length(ff.t),1);
    ff.hl_covS = nan(length(ff.t),1);
    ff.hl_covW(hh) = ff.hl_cov(hh);
    ff.hl_covS(hh) = ff.hl_cov_sds(hh);
   
                 
    %%% fields now saved in daily array: tau_cov, tau_ida, tau_idb, hs_id, hl_id, hs_covW, hl_covW, hs_covS, hl_covS
    %%% So, remove the fields that we don't really need to carry along, now
    %%% that the quality controlled and correct-unit id and cov fluxes are
    %%% included
%     fields_to_remove = {'wu_cov';'usid_a';'usid_b';'tsid_a';'qsid_a';'hlwebb';'hs_cov';'hs_cov_sds';'hl_cov';'hl_cov_sds'};
%     ff = rmfield(ff,fields_to_remove);     
    
    %%% Rain
    
        % Daily wxt rain accum
            temp = da(:,171);
        %     temp = ff.prate_wxt;
        %     temp = nan(1440,1);
        %     temp(isnan(temp)) = 0;
            accm_wxt = cumsum(temp)/6; % mm/hr to mm/10 min
        %     ff.paccum_wxt = accum_wxt;
    
        % Daily psl org_rate rain accum
        %     temp = da(:,24);
            temp = ff.prate;
            temp(isnan(temp)) = 0;
            accm = cumsum(temp)/6;
            ff.paccum = accm;
            
        % Daily ship rain accum
        %     temp = da(:,11);
        %     accm_s = cumsum(temp)/6;
    
        
    %%% Bulk surface fluxes
    ff.hr = real(ff.hr); % for some reason these had imaginary components??
    %%% prior code should have already taken care of sign conventions for
    %%% this version, in which negative = down
%     ff.lw_dn = -ff.lw_dn; % change to sign convention where negative = down
%     ff.lw_s = ff.lw_s; % this sign was already negative... not sure where that happened but whatever. 
        
    % OG yrnl = 0.97*(5.67e-8*(yts + tdk).^4 - yrl);

    ff.lw_net = 0.97*(5.67e-8*(ff.sst+C2K).^4 - ff.lw_dn); % net longwave...    lw_net = lw_up - lw_dn ... net is + upward cooling ocean , - downward heating ocean 
    ff.lw_up = ff.lw_dn + ff.lw_net;                       % upward longwave... lw_up = lw_dn + lw_net

%     ff.sw_dn = -ff.sw_dn; % change to sign convention where negative = down
%     ff.sw_s = -ff.sw_s; % change to sign convention where negative = down
    ff.sw_net = 0.955*ff.sw_dn;                             % net shortwave     sw_net = sw_up - sw_dn .... net is + downward heating ocean, - upward cooling ocean
    ff.sw_up = ff.sw_dn - ff.sw_net;                        % upward shortwave  sw_up = sw_dn - sw_net
    
    % OG hnet = 0.955*rs_best - rnl - hsb - hlb - RF;           % original net heat calc from coare output
    ff.hnet = ff.sw_net - ff.lw_net - ff.hs - ff.hl - ff.hr;    % net

    % compute sea surface specific humidity with SST
    ff.qs = qsea_p(ff.sst,ff.slp);
    
    % compute air-sea T and q gradients
    ff.dT = ff.sst - ff.T10;
    ff.dq = ff.qs - ff.q10;
    
    %%% meta data for ff plotting later
    ff.cruise_str = cruise_str;
    ff.PosLims = PosLims;
    ff.ztsg1_ship = ztsg1_ship;
    ff.ztsg2_ship = ztsg2_ship;
        
    disp(['FINISHED motcorr ',cruise,' yearday ',sprintf('%03i',ddd),'.']);

            
    %% DAILY flux file prep to save = all hours
    
    % put all variable names in alphabetical order;
    ff = orderfields(ff);
    
    ff.readme = {'PISTON 2019 field project 10-min data of flux, meteorology, seawater';...
                'from sensors by NOAA ESRL PSL, APL-UW, R/V Sally Ride';
                'Produced by Elizabeth Thompson with run_motcorr.m';....
                datetime;
                'Contact: elizabeth.thompson@noaa.gov';
                };


    % base file name
%     path_output = fullfile(data_drive,cruise,ship,'flux','Field_Processed','flux_met_sea_v0');
    
    fname = [cruise '_10min_flux_met_sea_data_' sprintf('%02i_%02i_jd%03i',Vdate(2),Vdate(3),ddd)];
    dfl_mat = [path_output '/' fname '.mat'];
    save(dfl_mat,'ff');    
    allfields = fields(ff);
    nff = length(allfields);
    
    %%% save text file of all these fields, minus the 1D fields: readme and sensor heights, and matlab t.
    %%% save text file of all these fields, minus the 1D fields: readme and sensor heights, and matlab t.
    f_fields = fields(ff);
    nf1 = zeros(1,length(f_fields));
    for i = 1:length(f_fields)
        nf1(i) = length(ff.(f_fields{i}));
    end
    wh_1D = find(nf1 < 50);
    ff2 = rmfield(ff,vertcat('t',f_fields(wh_1D)));
    

    %%% add in a human readable time array
    [ff2.year, ff2.month, ff2.day, ff2.hour, ff2.minute] = datevec(ff.t);
    ff2 = orderfields(ff2);

    %%% save text file
    k = mat2txt(fname, 'n', [path_output '/'] , ff2);
    %%% save header for text file
    
    %%% put entire structure into same ordered array with columns and rows
    fields_ff2 = fields(ff2);
    nff2 = size(fields_ff2);
    k = repmat('                    ',length(ff2.jd),1);
    for j = 1:nff2
    k(j) = fprintf('%20s',char(fields_ff2(j)));
    end
    disp('');
    disp(['num vars in txt file: ' sprintf('%i',nff2)]);
    disp(['num vars in mat file: ' sprintf('%i',nff)]);

    %%% save that in a text file manually and call it the header. 
    
%     %%% open .txt file with same name as .mat, then print vectors
%     dfl = [[path_output '/'] 'header_' fname '.txt'];
% %     fprintf('%20s',char(cfields_named(1)))
%     
%     flist = fopen(dfl,'w');
%     fmtStr = '';
% %     fmtStr = strcat('\r\n', fmtStr);
%     for i = 1:nff2
%        fmtStr = strcat('%20s ', fmtStr);
%     end
%     fprintf(flist,fmtStr,fields_ff2);
%     fclose('all');
%     disp(['Text file written as ' dfl]);
    
    %% Daily plots 
        
    %%% test plotting program made from this program
%     plot_motcorr(ff,path_dailyPlots);

    
end % end daily loop

% concatenate daily files into one whole-cruise file
filecat(fullfile(path_da,['da_*_',dcname,'.txt']), fullfile(path_proc_data_flux,...
    [cruise,'_da_',dcname,'.txt']),'Overwrite');
filecat(fullfile(path_cr,['cr_*_',dcname,'.txt']), fullfile(path_proc_data_flux,...
    [cruise,'_cr_',dcname,'.txt']),'Overwrite');
filecat(fullfile(path_sp,['sp_*_',dcname,'.txt']), fullfile(path_proc_data_flux,...
    [cruise,'_sp_',dcname,'.txt']),'Overwrite');
filecat(fullfile(path_lag,['lag_*_',dcname,'.txt']), fullfile(path_proc_data_flux,...
    [cruise,'_lag_',dcname,'.txt']),'Overwrite');
filecat(fullfile(path_mu,'w_Ts','mu_*.txt'), fullfile(path_proc_data_flux,...
    [cruise,'_mu_',dcname,'.txt']),'Overwrite');
filecat(fullfile(path_mu,'q','mu_*.txt'), fullfile(path_proc_data_flux,...
    [cruise,'_mu_q_',dcname,'.txt']),'Overwrite');

disp('END OF motcorr PROCESSING')
