%{
 Process daily met files, make diagnostic plots, and save intermediate
 1-min and 10-min data files.

 Updated May 2021 EJT

 Calls:
  read_met_day:     PSL met data from loggers 1 & 2 @ 1 min     from 1 min data
  read_gps_day:     PSL gps course & speed @ 1 Hz               from 10 Hz data
  read_hed_day:     PSL gps heading & pitch @ 1 Hz              from 10 Hz data
  read_motion_day:  PSL MotionPak @ 1 Hz                        from 10 Hz data
  read_sonic_day:   PSL sonic wind @ 1 Hz                       from 10 Hz data
  read_licor_day:   PSL licor fast water vapor / CO2 @ 1 Hz     from 10 Hz data
  read_ship_day:    ship data @ 1 Hz                            from 1 Hz or sometimes faster data
  read_wxt_day      PSL WXT met data @ 1 Hz                     from 1 min data
%}

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

%% set day-of-year range for this run
% for atomic 2020, use 9 to 44
% for piston 2019, use 248 to 268
% jdStart = 248; jdStop = 268;
jdStart = 268; jdStop = 268;

min10 = datenum(2018,8,27,0,10,0) - datenum(2018,8,27,0,0,0);

%% set paths
plotit = true; % for doing plots.
prtit = true;  % for saving plots (plotit must also be true)
graphformat = '.png';  % select graphics format files
graphdevice = '-dpng'; % select graphic device

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

% paths relative to data_drive - these should be preexisting
path_readme = fullfile(data_drive,cruise,ship,'flux','readme');
mkdir(path_readme);
path_python = fullfile(path_prog,'python');
path_raw_data = fullfile(data_drive,cruise,ship,'flux','Raw');

% define folders for saving data and plots
path_proc_data = fullfile(data_drive,cruise,ship,'flux','Processed');
mkdir(path_proc_data);
mkdir(fullfile(path_proc_data,'legacy','flux'));
mkdir(fullfile(path_proc_data,'legacy','gps'));
mkdir(fullfile(path_proc_data,'legacy','met'));
mkdir(fullfile(path_proc_data,'legacy','rad'));
mkdir(fullfile(path_proc_data,'v0_1min'));
mkdir(fullfile(path_proc_data,'v0_10min'));

path_raw_images = fullfile(data_drive,cruise,ship,'flux','Raw_Images');
mkdir(path_raw_images);
mkdir(fullfile(path_raw_images,'qSpectra'));
mkdir(fullfile(path_raw_images,'Rwspd'));
mkdir(fullfile(path_raw_images,'sonicQC'));
mkdir(fullfile(path_raw_images,'tSpectra'));
mkdir(fullfile(path_raw_images,'uSpectra'));
mkdir(fullfile(path_raw_images,'wSpectra'));
mkdir(fullfile(path_raw_images,'SL_pressure'));
mkdir(fullfile(path_raw_images,'Temps'));
mkdir(fullfile(path_raw_images,'RH'));
mkdir(fullfile(path_raw_images,'Wind'));
mkdir(fullfile(path_raw_images,'Rainrate'));
mkdir(fullfile(path_raw_images,'T_RH_fan'));
mkdir(fullfile(path_raw_images,'Track_plot'));
mkdir(fullfile(path_raw_images,'COG_SOG'));
mkdir(fullfile(path_raw_images,'Heading'));
mkdir(fullfile(path_raw_images,'pitch'));
mkdir(fullfile(path_raw_images,'Relative_Wind'));
mkdir(fullfile(path_raw_images,'True_Wind'));
mkdir(fullfile(path_raw_images,'Motion'));
mkdir(fullfile(path_raw_images,'Licor_AGC'));
mkdir(fullfile(path_raw_images,'Heat_Fluxes'));
mkdir(fullfile(path_raw_images,'u-star'));
mkdir(fullfile(path_raw_images,'IR_flux'));
mkdir(fullfile(path_raw_images,'Solar_flux'));
mkdir(fullfile(path_raw_images,'check'));


%% other initialization settings

%%%% Note on ship seawater measurements for PISTON 2019:
% '  Tsnake = NOAA sea snake at 0.05 m depth, 5 cm. Deployed at 5 Z on 6 Sept.';... 
% '  ship bow intake = starboard hole at 3.6 m depth, 0.3-0.6 m forward and below bow thruster vent,';...
% '         feeds into bow thruster room.';...
% '  seachest intake = port vent between 2.3-3.4 m depth, 0.6-0.9 m aft and just above bilge keel,';...
% '         feeds into engine room.';...
% '  TSG 1 = TSG located in bow thruster room (SBE45).';...
% '         At beginning of cruise, this TSG sampled bow water through a direct line. ';...
% '         After switch, this TSG re-sampled sea chest water after TSG #2.';...
% '  TSG 2 = TSG located in main lab.';...
% '         At beginning of cruise, this TSG re-sampled bow water after TSG #1.';...
% '         After switch, this TSG sampled seachest water through a direct line.';...
% '  The switch in plumbing happened between 00:00 and 00:15 on 12 Sept because the bow was pitching';...
% '         so much in rough seas that the bow pump was losing its prime and taking in air.';...
% '         For remainder of cruise, we decided to use sea chest water (TSG#2) instead of the bow water (TSG#1).';...

% grab remaining ship data "MET" that weren't saved in data aquisition
% files (they didn't contain all the needed data automatically)
load('/Users/eliz/DATA/PISTON_2019/Sally_Ride/ship/MET_SallyRide_PISTON_2019.mat');

%%% compute MET.jd time step, if it hasn't already been calculated
% MET.jd = MET.dum - datenum(yr, 0,0,0,0,0);
%%% or this way ... hard way
% if ~isfield(MET, 'jd')
%     [yy0, mm0, dd0, hh0, min0, ss0] = datevec(MET.dnum);
%     jd0 = md2yd_vn(yy0,mm0,dd0);
%     dayfraction = days(duration(hh0, mm0, ss0));
%     % dayfraction_2 = (ihr + (imin + (isec./60))./60)./24;  %% equivalent
%     % sprintf('%10.5f',dayfraction(2))
%     MET.jd = jd0 + dayfraction;
% end

%% loop through days
for ddd = jdStart:jdStop
    %% run python scripts on gps and WXT raw files
    [m,d] = yd2md(yr, ddd);
    Vdate = [yr, m, d];
    path_working_ddd = fullfile(path_raw_data,[yr_st(3:4),sprintf('%03i',ddd)]);
    % convert PSL gps and met3 files to gprm and wxt format with python script
    % python language must be installed on the computer.
    cd(path_working_ddd);

    files = dir('gps*.txt');
    if ~isempty(files)
        system(['python ',fullfile(path_python,'parseGpsFiles.py')]);
    end

    % why is this commented out?
%     files = dir('hed0*.txt');
%     if ~isempty(files)
%         system(['python ',fullfile(path_python,'parseHedFiles.py')]);
%     end
    
    files = dir('met3*.txt');
    if ~isempty(files)
        system(['python ',fullfile(path_python,'parseWXTFiles.py')]);
    end
    
    %%%=====================================================================
    %% read data - this takes time...
    
    % read ship (used to be called "SCS") data - from ship at 1 Hz 86400 by numvars: 1 sec
    % # fields changes each cruise depending on how much ship data is provided and what type
    ship_day = read_ship_day(path_working_ddd,ddd,yr_st,PosLims,ship_adj,zp_ship,zq_ship);
    % test with read_ship_practice or with this command:
    %  ship_day = read_ship_day('/Users/eliz/DATA/ATOMIC/Brown/flux/Raw/20010',10,'2020',[45 63 5 15],[0 0],15.6337)
        
    % read PSL 1-min avg bulk met data - metm is 1440x17 array: 1 min
    metm = read_met_day(path_working_ddd,ddd,yr_st,zp);

    % read PSL gps - gprm is 1Hz 86400x7 array: 1 sec
    gprm = read_gps_day(path_working_ddd,ddd,yr_st,PosLims);

    % read PSL heading-pitch - hedm is 1Hz 86400x3 array: 1 sec
    % array for MISOBOB - heading not working
    hedm = read_hed_day(path_working_ddd,ddd,yr_st);

    % read PSL motionPak - motm is 1Hz 86400x13: 1 sec
    motm = read_motion_day(path_working_ddd,ddd,yr_st);
   
    % MOTION PLOT - 1 Hz
        figure;
        subplot(3,2,1);plot(motm(:,1),motm(:,2),'.','markersize',2);ylabel('accx (m.s^-^2)');
        ylim([-5 5]);set(gca(gcf),'XTick',ddd:2/24:ddd+1); datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        grid;
        subplot(3,2,2);plot(motm(:,1),motm(:,3),'.','markersize',2);ylabel('accy (m.s^-^2)');
        ylim([-5 5]);set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        grid;
        subplot(3,2,3);plot(motm(:,1),motm(:,4),'.','markersize',2);ylabel('accz (m.s^-^2)');
        ylim([5 15]);set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        grid;
        subplot(3,2,4);plot(motm(:,1),motm(:,5),'.','markersize',2);ylabel('ratex (rad/s)');
        ylim([-5 5]*.02);set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        grid;
        subplot(3,2,5);plot(motm(:,1),motm(:,6),'.','markersize',2);ylabel('ratey (rad/s)');xlabel('Hour (UTC)');
        ylim([-5 5]*.02);set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        grid;
        subplot(3,2,6);plot(motm(:,1),motm(:,7),'.','markersize',2);ylabel('ratez (rad/s)');xlabel('Hour (UTC)');
        ylim([-5 5]*.02);set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        grid;
        annotation(gcf,'textbox',[0.35 0.96 0.8 0.02462],'String',sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSL motion pack',cruise_str,Vdate(1),Vdate(2),Vdate(3),ddd),'Fontsize',18,'FontWeight','Bold','FitBoxToText','off','LineStyle','none','Interpreter','none');
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA ESRL PSL'},'FontSize',14,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_raw_images,'Motion',['MotionPak_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end
    
    % read PSL sonic - sonm is 1Hz 86400x5: 1 sec data
    sonm = read_sonic_day(path_working_ddd,ddd,yr_st,sonicmodel,...
         rotationsonic,fsonic,prtit,path_raw_images,graphformat,graphdevice); 

    % read PSL Licor - licm is 1Hz 86400x8: 1 sec data 
    licm = read_licor_day(path_working_ddd,ddd,yr_st,...
                flicor,prtit,path_raw_images,graphformat,graphdevice);  
     
    % read wxt data - wxtm is 1440 x 9 array: 1 min data
    wxtm = read_wxt_day(path_working_ddd,ddd,yr_st,zwxt);

    close all
           
    %% 1-min averages for some variables
    disp('doing work on 1 min vars');
    
    %%% time
    jd_1min = metm(:,1); % ref 1-min timestamp
    delta = double(1.0/1440);
    last = ddd + 1440*delta;
    jd_1bin = (ddd:delta:last);	
    jd_1bin = jd_1bin'; % ref 1-min bin edges
    t_1min = datenum(yr,0,0,0,0,0) + jd_1min;
    
    delta_10 = double(10.0/1440);
    last_10 = ddd + 144*delta_10;
    jd_10bin = (ddd:delta_10:last_10);	
    jd_10bin = jd_10bin'; % ref 10-min bin edges 
    jd_10min = jd_10bin(1:end-1); % reference 10-min timestamp
    t_10min = datenum(yr,0,0,0,0,0) + jd_10min;
    
    
    %% 1-min PSL data
    
%     %%% motion data...saved for high resolution run_motcorr program later
%     motm1 = interval_avg(motm(:,1), motm(:,2:7), jd_1bin); % average everything
%     accX = motm1(:,2);
%     accY = motm1(:,3);
%     accZ = motm1(:,4);
%     rateX = motm1(:,5);
%     rateY = motm1(:,6);
%     rateZ = motm1(:,7);

    %%% PSL gps data, sog, cog
    gprm1 = interval_avg(gprm(:,1), gprm(:,2:7), jd_1bin); % average everything
    lat = gprm1(:,6); %
    lon = gprm1(:,7); %
    sogN = gprm1(:,4); % mean ship speed toward N
    sogE = gprm1(:,5); % mean ship speed toward W
    sog = sqrt(sogN.^2 + sogE.^2);  % recompute sog from speed components
    cog = atan2(sogE, sogN)*r2d; % recompute cog from speed components
    cog = mod(cog + 360, 360);
    gprm1(:,2) = cog; % replace with updated 1-min cog/sog
    gprm1(:,3) = sog;
    
    %%% PSL heading: from 1 Hz to 1 min sampling. Keep in hedm structure
    %%% since that's what is averaged to 10 min later on.
    [hedN_10Hz, hedE_10Hz] = pol2cart(hedm(:,2)*d2r,1);  % N/E heading components at 1 Hz in radians
    hedm = [hedm, hedN_10Hz, hedE_10Hz];                   % add to array of 1 Hz data
    hedm1 = interval_avg(hedm(:,1), hedm(:,2:5), jd_1bin); % average everything to 1 min
    hedN1 = hedm1(:,4);
    hedE1 = hedm1(:,5);
    hed = atan2(hedE1, hedN1)*r2d; % recompute hed from avg N/E components in radians, then convert back to degrees
    hed = mod(hed + 360, 360);
    hedm1(:,2) = hed; % replace with updated 1-min hed
    [hedN, hedE] = pol2cart(hed*d2r,1);  % recomputed N/E heading components with 1 min avgs to be consistent
    
    %%% PSL pitch from heading
    pitch = hedm1(:,3); % this is at 1 min
    ii = find(isnan(pitch)); 
    pitch(ii) = nanmean1(pitch);
    
    %%% sonic
    sonm1 = interval_avg(sonm(:,1), sonm(:,2:5), jd_1bin); % average all sonm
    Tsonic = sonm1(:,5);
    
    % correct ORG offset
    orgV_bkgd = nanmedian1(despike2(metm(:,17))); % single value
    orgV_adj = 0.06484 - orgV_bkgd; % single value
    org_recalc = (metm(:,17)+orgV_adj).^1.87*25 - 0.15; % recompute rain rate after adjustment
    org_recalc(org_recalc<0.05) = 0; % set lower limit on org rain rate sensitivity
    prate = org_recalc;
    prate(isnan(org_recalc) == 1) = 0; % reset nan to zero
    orgV = metm(:,17); % voltage 
    orgV_desp = despike2(metm(:,17)); % voltage despiked
    prate_orig = metm(:,13); % orgiginal rain rate estimated from ORG
    temp = prate;
    temp(isnan(temp)) = 0;
    paccum = cumsum(temp)/60; % 60 min per hour
    
    %%% relative and true wind from sonic and PSL gps, heading, sog, cog
    % Apply flow distortion corrections for mean spd using 0.95 & 1.15 
    % coefficients. Note: these are based on global class vessels like
    % R/V Revelle, R/V Thompson, and NOAA Ship Brown. They could need to be
    % adjusted for smaller or different profile ships. Also use ship
    % heading, sog, and cog, for adjusting relative wind. We don't know 
    % whether the w component or w^2 in vector average/sum is needed 
    % to compute wspd... we assume not though because we are instead 
    % correcting for flow distortion with the coefficients. Coefficients
    % were reduced to account for increase in flow distortion on R/V ride.
    % The change was arrived at by inspecting da_red plots to come up with
    % a percentage increase/decrease ~ 15% decrease in coefficients 
    
%     [rspd,rdir] = uv_to_sd(sonm1(:,2)/0.95,sonm1(:,3)/1.15); % for global class ships
%    %% chris had writfdsa  Eadzy that wasn't 15% c kj
    [rspd,rdir] = uv_to_sd(sonm1(:,2)/0.86,sonm1(:,3)/1.12); % for Sally Ride Ocean Class vessel
    [rUn, rUw] = sd_to_uv(rspd, rdir);
    rdir(rdir>180) = rdir(rdir>180)-360;    % +/- 180 deg format
    [wspd,wdir,Un, Uw] = uv_rel_to_sd_true(sonm1(:,2)/0.86,sonm1(:,3)/1.12,hed,cog,sog);
   

    %% other met data

    Ta = metm(:,2);
    rh = metm(:,3);
    Tsnk = metm(:,4);
    sw_dn_1 = -metm(:,5); % should be positive value... PISTON IR and solar fluxes were flipped, so apply neg sign
    lw_dn_1 = metm(:,6); % should be positive value... 1st guess is from met reader and wrong for now. It is re-fixed in fix_met_sea program after this
    sw_dn_2 = -metm(:,7); % should be positive value... PISTON was flipped, so apply neg sign
    lw_dn_2 = metm(:,8); % should be positive value... 1st guess is from met reader and wrong for now. It is re-fixed in fix_met_sea program after this
    lw_T_case_1 = metm(:,9);
    lw_T_dome_1 = metm(:,10);
    lw_T_case_2 = metm(:,11);
    lw_T_dome_2 = metm(:,12);
    psealevel = metm(:,14);  % sea level press
    pa = metm(:,20);  % measured press at height zp
    pa_at_zq = psealevel - 0.125*zq;
    aspir_TRH = metm(:,15);
    lw_therm_1 = metm(:,18); % measured value should be negative voltage. If not, flip sign when using this value to compute IR downwelling radiation in fix_met_sea.
    lw_therm_2 = metm(:,19); % measured value should be negative voltage. If not, flip sign when using this value to compute IR downwelling radiation in fix_met_sea.
    qa = qair_p(Ta,rh,pa_at_zq);
    qsnk = qsea_p(Tsnk,psealevel);
    % air density (kg/m3) and H2O number density (mmol/m3)
    rhoa = air_density(Ta,pa_at_zq,rh); %%% this gets recalculated in COARE below
    H2O = qa.*rhoa/0.018016;
            
    % average licm
    licm1 =  interval_avg(licm(:,1), licm(:,2:end), jd_1bin); 
%     licor_AGC = licm1(:,5);
    
    % 1-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 = licm1(:,1);
    
    licor_H2O = licm1(:,2); % Licor_H20, H2O density, mmol/m3 
    licor_Tbox = licm1(:,3); % Licor box T, C
    licor_Pbox = licm1(:,4); % Licor box pa, kPa
    licor_AGC = licm1(:,5); % Licor AGC, mean diagnostic value
    licor_CO2 = licm1(:,6); % Licor CO2 density, mmol/m3
    
    % there are SO many bad values in Tbox and Pbox. They correspond to bad
    % values in CO2 but not H2O, q, or rhoa... we use Ta and Pa from PSL
    % instead to compute licor met quantities below. Don't bother QC-ing
    % Pbox and Tbox. 
%     licor_Tbox(licor_Tbox > 40) = nan;
%     licor_Pbox(licor_Pbox*10 > 1030) = nan;
    
    %%% these are redone with 10 Hz data in motcorr, so are commented out
%     licor_H2O_std = interval_std_var(licm(:,1), licm(:,2), jd_1bin); % standard deviation of Licor H2O, mmol/m3
%     licor_Pbox_std = interval_std_var(licm(:,1), licm(:,4), jd_1bin); % standard deviation of Licor box pa, kPa

    licor_H2O_MR = (licor_H2O*1e-3).*Rgas_universal.*(Ta+C2K)./(pa_at_zq*1e-1); % WV MR converted from mmol/m3 to mol water / mol dry air
    licor_rhoa_dry = (pa_at_zq*1e3).*(1-licor_H2O_MR*(1-epsilon))./(Rgas.*(Ta+C2K)); % dry air density in licor [kg/m3]... conserved
    licor_qa_dry = (licor_H2O*1e-3).*Mw./licor_rhoa_dry; % mmol/m3 to g/kg dry air only
    licor_qa = (licor_H2O*1e-3).*Mw./rhoa; % mmol/m3 to g/kg moist air mix
    licor_CO2_MR = (licor_CO2).*Rgas_universal.*(Ta+C2K)./(pa_at_zq*1e-1); % C02 MR converted from mmol/m3 to ppm... micromol C02 / mol dry air
    licor_rh = relhum([Ta,licor_qa,pa_at_zq]);
    
    %%% 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 that are computed and saved here. 
    wh_good_licor = find(licor_AGC < 60);
    wh_bad_licor = find(licor_AGC >= 60);
    licor_rhoa_dry(wh_bad_licor) = nan; % kg / m3
    licor_qa_dry(wh_bad_licor) = nan; % g / kg
    licor_qa(wh_bad_licor) = nan; % g / kg
    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
    licor_CO2_MR(wh_bad_licor) = nan; % ppm
    licor_rh(wh_bad_licor) = nan; % %

    %% 1-min WXT data
    rUn_wxt = wxtm(:,8);
    rUw_wxt = wxtm(:,9);
    [rspd_wxt,rdir_wxt] = uv_to_sd(rUn_wxt,rUw_wxt); % recompute wspd/wdir with cog, sog, hed (which WXT doesn't do in its reader, since it doesn't measure it)
    rr = rdir_wxt>180; 
    rdir_wxt(rr) = rdir_wxt(rr)-360;
    [wspd_wxt,wdir_wxt,Un_wxt,Uw_wxt] = uv_rel_to_sd_true(wxtm(:,8),wxtm(:,9),hed,cog,sog);
    Ta_wxt = wxtm(:,4);
    rh_wxt = wxtm(:,5);
    pa_wxt = wxtm(:,6);
    qa_wxt = qair_p(Ta_wxt, rh_wxt, pa_wxt); % spec humidity
    % air density (kg/m3) and H2O number density (mmol/m3)
    rhoa_wxt = air_density(Ta_wxt,pa_wxt,rh_wxt);
    H2O_wxt = qa_wxt.*rhoa_wxt/0.018016;
    psealevel_wxt = wxtm(:,6)+0.125*zwxt;
    prate_wxt = wxtm(:,7);  % mm/hr
    prate_wxt(isnan(prate_wxt) == 1) = 0; % reset nan to zero
    temp = prate_wxt;
    temp(isnan(temp)) = 0;
    paccum_wxt = cumsum(temp)/60; % 60 X 1 min segments per hour
    
    %% 1-min ship data 
    
%     The reader checks to make sure sog is in m/s. This could also be done
%     here instead or later in fixit.m

    %%% for PISTON only, change units of T from K to C, eliminate noisy
    %%% (negative) values of lw_dn. NOTE: do most of the other fixes later
    %%% in fixit.m afterwards. 
    ship_day.lw_T_case = ship_day.lw_T_case-273.15;
    ship_day.lw_T_dome = ship_day.lw_T_dome-273.15;
    ship_day.lw_dn(ship_day.lw_dn < 0 | ship_day.lw_dn > 500) = nan;        

    %%% ship heading, cog, sog, wind: from 1 Hz to 1 min sampling, into ship_avg structure.
    % recomputed N/E heading components from 1 min avg magnitude and direction
    [ship_day.hedN, ship_day.hedE] = pol2cart(ship_day.hed*d2r,1);  
    [ship_day.Un, ship_day.Uw] = sd_to_uv(ship_day.wspd, ship_day.wdir);
    [ship_day.rUn, ship_day.rUw] = sd_to_uv(ship_day.rspd, ship_day.rdir);
    
    clear ship_avg;
    mfields = fields(ship_day);
    for j = 1:length(mfields) 
        x = ship_day.(mfields{j});
        ship_avg.(mfields{j}) = interval_avg_var(ship_day.jd, x, jd_1bin);
    end

    lat_s = ship_avg.lat;
    lon_s = ship_avg.lon;
    sog_s = sqrt(ship_avg.sogN.^2 + ship_avg.sogE.^2);  % recompute sog from speed components
    cog_s = atan2(ship_avg.sogE, ship_avg.sogN)*r2d;  % recompute cog from speed components
    cog_s = mod(cog_s + 360, 360);
    ship_avg.cog = cog_s; % replace with updated 1-min cog/sog
    ship_avg.sog = sog_s;
    hed_s = atan2(ship_avg.hedE, ship_avg.hedN)*r2d; % recompute avg hed from avg N/E components
    hed_s = mod(hed_s + 360, 360);
    ship_avg.hed = hed_s; % replace with updated 1-min hed

    [rspd_s,rdir_s] = uv_to_sd(ship_avg.rUn, ship_avg.rUw); 
    rdir_s(rdir_s>180) = rdir_s(rdir_s>180)-360;    % +/- 180 deg format
    %%% use same rUn and rUw components, and the recently recalculated 10-min avgs of hed, sog, cog, to do 10-min true wind
    [wspd_s,wdir_s,Un_s,Uw_s] = uv_rel_to_sd_true(ship_avg.rUn,ship_avg.rUw,ship_avg.hed,ship_avg.cog,ship_avg.sog);
    ship_avg.wspd = wspd_s;
    ship_avg.rspd = rspd_s;
    ship_avg.wdir = wdir_s;
    ship_avg.rdir = rdir_s;
    ship_avg.Un = Un_s;
    ship_avg.Uw = Uw_s;  
    
    %% 1 min remaining ship data
    % MET file with missing ship variables that are not present DAQ file.
    % [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 = {'Tsea1';'Tsea2';'Ssea1';'Ssea2';...
        'Fsea1';'Fsea2';...
        'ID_Fsea1';'ID_Fsea2';'spd_log_u';'spd_log_v';...
        'lw_therm'};

    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';...
        'LWR_theropile'};
    
    disp('start MET vars');
    for i = 1:length(MET_in_vars)
        thisvar = MET.(MET_in_vars{i});
        outvar = interval_avg_var(MET.jd(ID_u), thisvar(ID_u), jd_1bin);
        outvar_10 = interval_avg_var(MET.jd(ID_u), thisvar(ID_u), jd_10bin);
        eval([MET_out_vars{i} ' = outvar;']);  
        eval([MET_out_vars{i} '_10 = outvar_10;']);  
    end
    disp('end MET vars');

    ID_Fsea1_10 = round(ID_Fsea1_10);
    ID_Fsea2_10 = round(ID_Fsea2_10);
    ID_Fsea1 = round(ID_Fsea1);
    ID_Fsea2 = round(ID_Fsea2);
    
    %%% get std of spd_log at high res
    MET_spd_log_u = MET.ship_speed_transverse/1.944;
    MET_spd_log_v = MET.ship_speed_longitudinal/1.944;
    MET_spd_log = sqrt(MET_spd_log_u.^2 + MET_spd_log_v.^2);
    spd_log_std_s = interval_std_var(MET.jd(ID_u), MET_spd_log(ID_u), jd_1bin);
    spd_log_std_s_10 = interval_std_var(MET.jd(ID_u), MET_spd_log(ID_u), jd_10bin);

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

    %%% on this 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
    % there are no real matlab warnings here. these are already defined, matlab just doesn't know it
    Ssea1(Ssea1 < 30) = nan;
    Ssea2(Ssea2 < 30) = nan;
    
    ID_tsg(t_1min < tsgswitchtime1) = 1;  % use bow before switch
    ID_tsg(t_1min > tsgswitchtime2) = 2;  % use seachest after switch
    ID_tsg(t_1min >= tsgswitchtime1 & t_1min <= tsgswitchtime2) = 0;

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

    % there are no real matlab warnings here. these are already defined, matlab just doesn't know it
    Tsea1(ID_tsg == 0) = nan;
    Tsea2(ID_tsg == 0) = nan;
    Ssea1(ID_tsg == 0) = nan;
    Ssea2(ID_tsg == 0) = nan;
    Fsea1(ID_tsg == 0) = nan;
    Fsea2(ID_tsg == 0) = nan;
    
    % PISTON: speed log. Should be in m/s to match sog. It was originally in kt.
    spd_log_u = spd_log_u./1.944;
    spd_log_v = spd_log_v./1.944;
    spd_log_u_10 = spd_log_u_10./1.944;
    spd_log_v_10 = spd_log_v_10./1.944;
    
    %%% add remaining MET vars not accounted for in MET_in_vars
    %%% ship precip ... this is my best but still good attempt at rain rate
    %%% computed as a function of their accumulation gauge.
    paccum_s = interp1(MET.dnum1min, MET.paccum1min, t_1min);
    prate_s = cum2rate(paccum_s', t_1min')';
    prate_s(isnan(prate_s) == 1) = 0; % reset nan to zero
    spd_log_s = sqrt(spd_log_u.^2 + spd_log_v.^2);
    ID_tsg_s = ID_tsg;

    paccum_s_10 = interval_avg_var(jd_1min, paccum_s, jd_10bin);
    prate_s_10 = cum2rate(paccum_s_10', t_10min')';
    prate_s_10(isnan(prate_s_10) == 1) = 0; % reset nan to zero
    spd_log_s_10 = sqrt(spd_log_u_10.^2 + spd_log_v_10.^2);
    ID_tsg_s_10 = ID_tsg_10;

    
    %%% all the MET variables we want to save
    
    %%% these last ones aren't part of MET_in_vars, and already have _s names
    MET_more_vars = {'paccum_s';'prate_s';'spd_log_s';'spd_log_std_s';'ID_tsg_s'}; 
    
    %%% These are done but need _s names
    MET_vars_s = strcat(MET_out_vars, '_s');

    %%% all with _s now, combine:
    MET_all = [MET_vars_s; MET_more_vars];

    % save 1-min averages of these qc-ed 1-min data
    % Not included in MET_out_vars: spd_log, prate, paccum
    for i = 1:length(MET_vars_s)
        eval([MET_vars_s{i} ' = ' MET_out_vars{i} ';']); 
        eval([MET_vars_s{i} '_10 = ' MET_out_vars{i} '_10;']); 
    end
    
    %% 10-min ship data
    %%% all the MET variables we want to save ... redo some of the 10-min
    %%% vars for averaging / time / component reasons

    clear ship_avg_10;
    mfields = fields(ship_day);
    for j = 1:length(mfields) 
        x = ship_day.(mfields{j});
        ship_avg_10.(mfields{j}) = interval_avg_var(ship_day.jd, x, jd_10bin);
    end
        
    ship_avg_10.sog_std = interval_std_var(ship_day.jd, ship_day.sog, jd_10bin);    % sog std dev in ship speed
    ship_avg_10.hed_std = interval_std_var(ship_day.jd, unwrap(ship_day.hed*d2r), jd_10bin)*r2d;    % sog std dev in ship speed
    ship_avg_10.cog_std = interval_std_var(ship_day.jd, unwrap(ship_day.cog*d2r), jd_10bin)*r2d;    % sog std dev in ship speed
    ship_avg_10.ched_std = interval_std_var(ship_day.jd, sin(ship_day.hed*d2r), jd_10bin);    % sin of sog std dev in ship speed, in radians
    ship_avg_10.shed_std = interval_std_var(ship_day.jd, cos(ship_day.hed*d2r), jd_10bin);    % cos of sog std dev in ship speed, in radians
    
    ship_avg_10.sog_min = interval_min_var(ship_day.jd, ship_day.sog, jd_10bin);    % sog min dev in ship speed
    ship_avg_10.hed_min = interval_min_var(ship_day.jd, unwrap(ship_day.hed*d2r), jd_10bin)*r2d;    % sog min dev in ship speed
    ship_avg_10.cog_min = interval_min_var(ship_day.jd, unwrap(ship_day.cog*d2r), jd_10bin)*r2d;    % sog min dev in ship speed
    
    ship_avg_10.sog_max = interval_max_var(ship_day.jd, ship_day.sog, jd_10bin);    % sog max dev in ship speed
    ship_avg_10.hed_max = interval_max_var(ship_day.jd, unwrap(ship_day.hed*d2r), jd_10bin)*r2d;    % sog max dev in ship speed
    ship_avg_10.cog_max = interval_max_var(ship_day.jd, unwrap(ship_day.cog*d2r), jd_10bin)*r2d;    % sog max dev in ship speed
    
    ship_avg_10.wspd_max = interval_max_var(ship_day.jd, ship_day.wspd, jd_10bin);
    ship_avg_10.wspd_min = interval_min_var(ship_day.jd, ship_day.wspd, jd_10bin);
    ship_avg_10.wspd_std = interval_std_var(ship_day.jd, ship_day.wspd, jd_10bin);
    ship_avg_10.wdir_max = interval_max_var(ship_day.jd, ship_day.wdir, jd_10bin);
    ship_avg_10.wdir_min = interval_min_var(ship_day.jd, ship_day.wdir, jd_10bin);
    ship_avg_10.wdir_std = interval_std_var(ship_day.jd, ship_day.wdir, jd_10bin);

    ship_avg_10.rdir_std = interval_std_var(ship_day.jd, unwrap(ship_day.rdir*d2r), jd_10bin)*r2d;
    ship_avg_10.rspd_std = interval_std_var(ship_day.jd, ship_day.rspd, jd_10bin);
    
    ship_avg.sog_std = interval_std_var(ship_day.jd, ship_day.sog, jd_1bin);    
    ship_avg.hed_std = interval_std_var(ship_day.jd, unwrap(ship_day.hed*d2r), jd_1bin)*r2d;   
    ship_avg.cog_std = interval_std_var(ship_day.jd, unwrap(ship_day.cog*d2r), jd_1bin)*r2d;   
    ship_avg.ched_std = interval_std_var(ship_day.jd, sin(ship_day.hed*d2r), jd_1bin);    % cosine of heading, in radians
    ship_avg.shed_std = interval_std_var(ship_day.jd, cos(ship_day.hed*d2r), jd_1bin);    % sine of heading, in radians

    ship_avg.sog_min = interval_min_var(ship_day.jd, ship_day.sog, jd_1bin);    
    ship_avg.hed_min = interval_min_var(ship_day.jd, unwrap(ship_day.hed*d2r), jd_1bin)*r2d;    
    ship_avg.cog_min = interval_min_var(ship_day.jd, unwrap(ship_day.cog*d2r), jd_1bin)*r2d;   
    
    ship_avg.sog_max = interval_max_var(ship_day.jd, ship_day.sog, jd_1bin);   
    ship_avg.hed_max = interval_max_var(ship_day.jd, unwrap(ship_day.hed*d2r), jd_1bin)*r2d;   
    ship_avg.cog_max = interval_max_var(ship_day.jd, unwrap(ship_day.cog*d2r), jd_1bin)*r2d;   
    
    ship_avg.wspd_max = interval_max_var(ship_day.jd, ship_day.wspd, jd_1bin);
    ship_avg.wspd_min = interval_min_var(ship_day.jd, ship_day.wspd, jd_1bin);
    ship_avg.wspd_std = interval_std_var(ship_day.jd, ship_day.wspd, jd_1bin);
    ship_avg.wdir_max = interval_max_var(ship_day.jd, ship_day.wdir, jd_1bin);
    ship_avg.wdir_min = interval_min_var(ship_day.jd, ship_day.wdir, jd_1bin);
    ship_avg.wdir_std = interval_std_var(ship_day.jd, ship_day.wdir, jd_1bin);   
    
    ship_avg.rdir_std = interval_std_var(ship_day.jd, unwrap(ship_day.rdir*d2r), jd_1bin)*r2d;
    ship_avg.rspd_std = interval_std_var(ship_day.jd, ship_day.rspd, jd_1bin);
   
    % rest of ship met and nav data
    lat_s_10 = ship_avg_10.lat;
    lon_s_10 = ship_avg_10.lon;
    sog_s_10 = sqrt(ship_avg_10.sogN.^2 + ship_avg_10.sogE.^2);  % recompute sog from speed components
    cog_s_10 = atan2(ship_avg_10.sogE, ship_avg_10.sogN)*r2d;  % recompute cog from speed components
    cog_s_10 = mod(cog_s_10 + 360, 360);
    ship_avg_10.cog = cog_s_10; % replace with updated 1-min cog/sog
    ship_avg_10.sog = sog_s_10;
    hedE_s_10 = ship_avg_10.hedE;
    hedN_s_10 = ship_avg_10.hedN;
    hed_s_10 = atan2(hedE_s_10, hedN_s_10)*r2d; % recompute avg hed from avg N/E components
    hed_s_10 = mod(hed_s_10 + 360, 360);
    ship_avg_10.hed = hed_s_10; % replace with updated 1-min hed
    
    rh_s_10 = ship_avg_10.rh; % rh
    Ta_s_10 = ship_avg_10.Ta; % air t
    psealevel_s_10 = ship_avg_10.psealevel; % psealevel
    pa_at_zq_s_10 = psealevel_s_10 - (0.125*zq_ship); % 
    ship_avg_10.pa = psealevel_s_10 - (0.125*zp_ship); % p measured
    sw_dn_s_10 = ship_avg_10.sw_dn; % solar radiation
    lw_dn_s_10 = ship_avg_10.lw_dn; % ir radiation
    qa_s_10 = qair_p(Ta_s_10, rh_s_10, pa_at_zq_s_10); % spec humidity
    ship_avg_10.qa = qa_s_10;
    
    % air density (kg/m3) and H2O number density (mmol/m3)
    rhoa_s_10 = air_density(Ta_s_10,pa_at_zq_s_10,rh_s_10);
    H2O_s_10 = qa_s_10.*rhoa_s_10/0.018016;
    ship_avg_10.rhoa = rhoa_s_10;
    ship_avg_10.H2O = H2O_s_10;
    
    [rspd_s_10,rdir_s_10] = uv_to_sd(ship_avg_10.rUn, ship_avg_10.rUw); 
    rdir_s_10(rdir_s_10>180) = rdir_s_10(rdir_s_10>180)-360;    % +/- 180 deg format
    %%% use same rUn and rUw components, and the recently recalculated 10-min avgs of hed, sog, cog, to do 10-min true wind
    [wspd_s_10,wdir_s_10,Un_s_10,Uw_s_10] = uv_rel_to_sd_true(ship_avg_10.rUn,ship_avg_10.rUw,ship_avg_10.hed,ship_avg_10.cog,ship_avg_10.sog);
    ship_avg_10.wspd = wspd_s_10;
    ship_avg_10.rspd = rspd_s_10;
    ship_avg_10.wdir = wdir_s_10;
    ship_avg_10.rdir = rdir_s_10;
    ship_avg_10.Un = Un_s_10;
    ship_avg_10.Uw = Uw_s_10;   
    
    
    %% 10 min PSL data 
    % 10-min PSL gps, sog, cog from 1-Hz data
    gprm10 = interval_avg(gprm(:,1), gprm(:,2:7), jd_10bin); % average everything
    sogN_10 = gprm10(:,4); % ship speed toward N
    sogE_10 = gprm10(:,5); % ship speed toward W
    lat_10 = gprm10(:,6);
    lon_10 = gprm10(:,7);
    sog_10 = sqrt(sogN_10.^2 + sogE_10.^2);  % sog from speed components from GPS
    cog_10 = atan2(sogE_10, sogN_10)*r2d; % cog from speed components from GPS
    cog_10 = mod(cog_10 + 360, 360); % cog fom GPS
      
    %%% these are redone with 10 Hz data in motcorr, so are commented out
%     cog_std_10 = interval_std_var(gprm(:,1),unwrap(gprm(:,2)*d2r),jd_10bin)*r2d;
%     sog_std_10 = interval_std_var(gprm(:,1),gprm(:,3),jd_10bin);
%     cog_std = interval_std_var(gprm(:,1),unwrap(gprm(:,2)*d2r),jd_1bin)*r2d;
%     sog_std = interval_std_var(gprm(:,1),gprm(:,3),jd_1bin);
    
    % 10-min PSL heading from 1Hz raw data
    hedm10 = interval_avg(hedm(:,1), hedm(:,2:5), jd_10bin); % average everything
    hedN_10 = hedm10(:,4);
    hedE_10 = hedm10(:,5);
    hed_10 = atan2(hedE_10, hedN_10)*r2d; % recompute hed from avg N/E components
    hed_10 = mod(hed_10 + 360, 360);
    hedm10(:,2) = hed_10; % replace with updated 10-min hed
    
    %%% these are redone with 10 Hz data in motcorr, so are commented out
%     hed_std_10 = interval_std_var(hedm(:,1), unwrap(hedm(:,2)*d2r), jd_10bin)*r2d;   
%     hed_std = interval_std_var(hedm(:,1), unwrap(hedm(:,2)*d2r), jd_1bin)*r2d;   
    
    pitch_10 = hedm10(:,3); % this is at 1 min
    ii_10 = find(isnan(pitch_10)); 
    pitch_10(ii_10) = nanmean1(pitch_10);

%     % Another way to do the 10-min heading from 1 Hz components, used in motcorr
%     heading_interp = replace_NaN_nearest_neighbor(hed(:,2)*d2r);  % to radians
%     heading = heading_interp;
%     shed = sin(heading);    
%     ched = cos(heading);
%     ched = interval_avg(hed(:,1), ched, jd_10bin);
%     shed_avg = interval_avg(hed(:,1), shed, jd_10bin);
%     [hed_avg,~] = cart2pol(ched(:,2), shed_avg(:,2));
%     hed_avg = hed_avg*r2d;
%     hed_avg = mod(hed_avg+360, 360);
 
    % 10-min PSL met data from 1 Hz data
    metm10 = interval_avg(metm(:,1), metm(:,2:end), jd_10bin);
    Ta_10 = metm10(:,2);
    rh_10 = metm10(:,3);
    Tsnk_10 = metm10(:,4); 
    sw_dn_1_10 = -metm10(:,5); % should be positive. PISTON installation was flipped, so apply neg sign. 
    lw_dn_1_10 = metm10(:,6); % should be positive.
    sw_dn_2_10 = -metm10(:,7); % should be positive. PISTON installation was flipped, so apply neg sign. 
    lw_dn_2_10 = metm10(:,8); % should be positive.
    lw_T_case_1_10 = metm10(:,9);
    lw_T_dome_1_10 = metm10(:,10);
    lw_T_case_2_10 = metm10(:,11);
    lw_T_dome_2_10 = metm10(:,12);
    org_10 = metm10(:,13);
    psealevel_10 = metm10(:,14);  % sea level press
    pa_at_zq_10 = psealevel_10-(0.125*zq);
    pa_10 = psealevel_10-(0.125*zp);  % measured press
    aspir_TRH_10 = metm10(:,15);
    lw_therm_1_10 = metm10(:,18); % should be measured as negative value. If not, flip sign.
    lw_therm_2_10 = metm10(:,19); % should be measured as negative value. If not, flip sign.
    qa_10 = qair_p(Ta_10,rh_10,pa_at_zq_10);
    qsnk_10 = qsea_p(Tsnk_10,psealevel_10);
    % air density (kg/m3) and H2O number density (mmol/m3)
    rhoa_10 = air_density(Ta_10,pa_at_zq_10,rh_10); 
    H2O_10 = qa_10.*rhoa_10/0.018016;

    % recalculate ORG 10 min from 1 min org data QC'd above
    if all(isnan(org_10)); org_10 = 0; end
    org_carrier_10 = metm10(:,16);
    org_10_recalc = interval_avg_var(metm(:,1), org_recalc, jd_10bin);
    org_10_recalc(org_10_recalc<0.05) = 0;  % set lower limit on org rain rate sensitivity
    prate_10 = org_10_recalc;
    prate_10(isnan(org_10_recalc) == 1) = 0; % reset nan to 0
    temp = prate_10;
    temp(isnan(temp)) = 0;
    paccum_10 = cumsum(temp)/6; % 6 X 10 min segments per hour

    orgV_bkgd_10    = nanmedian1(despike2(interval_avg_var(metm(:,1), metm(:,17), jd_10bin))); % single value
    orgV_adj_10     = 0.06484 - orgV_bkgd_10; % single value
    orgV_10         = interval_avg_var(metm(:,1), metm(:,17), jd_10bin); 
    orgV_desp_10    = despike2(interval_avg_var(metm(:,1), metm(:,17), jd_10bin));
    prate_orig_10   = interval_avg_var(metm(:,1), metm(:,13), jd_10bin);

    % compute 10-min PSL winds from sonic 10Hz
    sonm10 = interval_avg(sonm(:,1), sonm(:,2:end), jd_10bin); % average everything
    Tsonic_10 = sonm10(:,5);
    
    % Apply flow distortion corrections for mean spd, using ship heading for corrections
    % increases wind speed for bow-on winds and decreases wind speed for
    % cross-ship winds. Agrees with results we have and from other ships. 
    [rspd_10,rdir_10] = uv_to_sd(sonm10(:,2)/0.86,sonm10(:,3)/1.12);
    rdir_10(rdir_10>180) = rdir_10(rdir_10>180)-360; % +/- 180 deg format
    [wspd_10,wdir_10,Un_10,Uw_10] = uv_rel_to_sd_true(sonm10(:,2)/0.86,sonm10(:,3)/1.12,hed_10,cog_10,sog_10);

    %% other met data at 10-min

    % 10-min PSL licor
    licm10 = interval_avg(licm(:,1), licm(:,2:end), jd_10bin);

    % 10-min wxt met 
    wxt10 = interval_avg(wxtm(:,1), wxtm(:,2:9), jd_10bin);      % average everything
    rUn_wxt_10 = wxt10(:,8);
    rUw_wxt_10 = wxt10(:,9);
    [rspd_wxt_10,rdir_wxt_10] = uv_to_sd(rUn_wxt_10,rUw_wxt_10); % recompute wspd/wdir with cog, sog, hed (which WXT doesn't do in its reader, since it doesn't measure it)
    rr = rdir_wxt_10>180; 
    rdir_wxt_10(rr) = rdir_wxt_10(rr)-360; % +/- 180 deg format
    [wspd_wxt_10,wdir_wxt_10,Un_wxt_10,Uw_wxt_10] = uv_rel_to_sd_true(wxt10(:,8),wxt10(:,9),hed_10,cog_10,sog_10);
    Ta_wxt_10 = wxt10(:,4);
    rh_wxt_10 = wxt10(:,5);
    pa_wxt_10 = wxt10(:,6);
    qa_wxt_10 = qair_p(Ta_wxt_10, rh_wxt_10, pa_wxt_10); % spec humidity
    % air density (kg/m3) and H2O number density (mmol/m3)
    rhoa_wxt_10 = air_density(Ta_wxt_10,pa_wxt_10,rh_wxt_10);
    H2O_wxt_10 = qa_wxt_10.*rhoa_wxt_10/0.018016;
    psealevel_wxt_10 = pa_wxt_10+0.125*zwxt;
    prate_wxt_10 = wxt10(:,7);  % mm/hr
    prate_wxt_10(isnan(prate_wxt_10) == 1) = 0; % reset nan to 0
    temp = prate_wxt_10;
    temp(isnan(temp)) = 0;
    paccum_wxt_10 = cumsum(temp)/6; % 6 X 10 min segments per hour

    %% 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_10_jd = licm10(:,1);

    licor_H2O_10 = licm10(:,2); % Licor_H20, WV density, mmol/m3 
    licor_Tbox_10 = licm10(:,3); % Licor box T, C
    licor_Pbox_10 = licm10(:,4); % Licor box pa, kPa
    licor_AGC_10 = licm10(:,5); % Licor AGC, mean diagnostic value
    licor_CO2_10 = licm10(:,6); % Licor CO2, density, mmol/m3

    % there are SO many bad values in Tbox and Pbox. They correspond to bad
    % values in CO2 but not H2O, q, or rhoa...
%     licor_Tbox_10(licor_Tbox_10 > 40) = nan;
%     licor_Pbox_10(licor_Pbox_10*10 > 1030) = nan;
    
    %%% these are redone with 10 Hz data in motcorr, so are commented out
%     licor_H2O_std_10 = interval_std_var(licm(:,1), licm(:,2), jd_10bin); % standard deviation of Licor H2O, mmol/m3
%     licor_Pbox_std_10 = interval_std_var(licm(:,1), licm(:,4), jd_10bin); % standard deviation of Licor box pa, kPa

    licor_H2O_MR_10 = (licor_H2O_10*1e-3).*Rgas_universal.*(Ta_10+C2K)./(pa_at_zq_10*1e3); % WV MR converted from mmol/m3 to mol water / mol dry air
    licor_rhoa_dry_10 = (pa_at_zq_10*1e3).*(1-licor_H2O_MR_10*(1-epsilon))./(Rgas.*(Ta_10+C2K)); % dry air density in licor [kg/m3]... conserved
    licor_qa_dry_10 = (licor_H2O_10*1e-3).*Mw./licor_rhoa_dry_10; % mmol/m3 to g/kg dry air only
    licor_qa_10 = (licor_H2O_10*1e-3).*Mw./rhoa_10; % mmol/m3 to g/kg moist air mix
    licor_CO2_MR_10 = (licor_CO2_10*1e3).*Rgas_universal.*(Ta_10+C2K)./(pa_at_zq_10*1e3); % C02 MR converted from mmol/m3 to ppm... micromol C02 / mol dry air
    licor_rh_10 = relhum([Ta_10,licor_qa_10,pa_at_zq_10]);

    %%% 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 in motcorr
    % program. These means are not used later. They are just here for
    % reference.
    wh_good_licor_10 = find(licor_AGC_10 < 60);
    wh_bad_licor_10 = find(licor_AGC_10 >= 60);
    licor_rhoa_dry_10(wh_bad_licor_10) = nan; % kg / m3
    licor_qa_dry_10(wh_bad_licor_10) = nan; % g / kg
    licor_qa_10(wh_bad_licor_10) = nan; % g / kg
    licor_H2O_10(wh_bad_licor_10) = nan;  % mmol/m3
    licor_H2O_MR_10(wh_bad_licor_10) = nan;  % mol/mol
    licor_CO2_10(wh_bad_licor_10) = nan;  % mmol/m3
    licor_CO2_MR_10(wh_bad_licor_10) = nan;  % micromol/m3 or ppm
    licor_Tbox_10(wh_bad_licor_10) = nan; % C
    licor_Pbox_10(wh_bad_licor_10) = nan; % kPa
    
    
    %% save 1-min results to text files, historical format
    % with added columns in ship, rad and met files for ship and WXT variables

    % 1-min PSL met data
    % jd pir1 psp1 Tc1 Td1 Tsnk Ta rh rain org_diag aspir_TRH Pmb rdir_wxt rspd_wxt ta_wxt rh_wxt Pmb_wxt rain_wxt U_wxt V_wxt
    fpath = fullfile(path_proc_data,'legacy','met',[cruise,'_met_1min_',]%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
    fout = fopen(fpath,'w');
    jazz = zeros(12,1440)*NaN;
    jazz(1,:) = jd_1min;
    jazz(2,:) = metm(:,6);      % pir1 W/m2
    jazz(3,:) = metm(:,5);      % psp1 W/m2
    jazz(4,:) = metm(:,9);      % Tc1 C
    jazz(5,:) = metm(:,10);     % Td1 C
    jazz(6,:) = metm(:,4);      % Tsnake C
    jazz(7,:) = metm(:,2);      % Ta vaisala C
    jazz(8,:) = metm(:,3);      % RH vaisala
    jazz(9,:) = prate;          % org mm/hr
    jazz(10,:) = metm(:,16);    % org carrier
    jazz(11,:) = metm(:,15);    % aspir_TRH on
    jazz(12,:) = metm(:,14);	% psealevel mb

    fprintf(fout,['%11.5f ',repmat('%9.2f ',1,12),' \n'],jazz);
    fclose all;

    % 1-min radiation data
    % jd pir1 pir2 psp1 psp2 Tc1 Td1 Tc2 Td2 Ta rh Tsnk imir imsol impir imtc imtd imir2
    fpath = fullfile(path_proc_data,'legacy','rad',[cruise,'_rad_1min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
    fout = fopen(fpath,'w');
    jazz = zeros(17,1440)*NaN;
    jazz(1,:) = jd_1min;
    jazz(2,:) = metm(:,6);           % pir1 W/m2
    jazz(3,:) = metm(:,8);           % pir2 W/m2
    jazz(4,:) = metm(:,5);           % psp1 W/m2
    jazz(5,:) = metm(:,7);           % psp2 W/m2
    jazz(6,:) = metm(:,9);           % Tc1 C
    jazz(7,:) = metm(:,10);          % Td1 C
    jazz(8,:) = metm(:,11);          % Tc2 C
    jazz(9,:) = metm(:,12);          % Td2 C
    jazz(10,:) = metm(:,2);          % Ta vaisala C
    jazz(11,:) = metm(:,3);          % RH vaisala
    jazz(12,:) = metm(:,4);          % Tsnake C
    jazz(13,:) = ship_avg.lw_dn;     % ship pir
    jazz(14,:) = ship_avg.sw_dn;     % ship psp
    jazz(15,:) = lw_therm;           % ship thermopile mV
    jazz(16,:) = ship_avg.lw_T_case; % ship T case
    jazz(17,:) = ship_avg.lw_T_dome; % ship T dome
    fprintf(fout,['%11.5f ',repmat('%9.3f ',1,17),' \n'],jazz);
    fclose('all');

    % 1-min gps data
    % jd lat lon sog cog hed pitch lat lon head cog sog
    fpath = fullfile(path_proc_data,'legacy','gps',[cruise,'_gpsnav_1min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
    fout = fopen(fpath,'w');
    jazz = zeros(12,1440)*NaN;
    jazz(1,:) = jd_1min;
    jazz(2,:) = gprm1(:,6);     % PSL decimal latitude, deg
    jazz(3,:) = gprm1(:,7);     % PSL decimal longitude, deg
    jazz(4,:) = gprm1(:,3);     % PSL GPS SOG, m/s
    jazz(5,:) = gprm1(:,2);     % PSL GPS COG, deg
    jazz(6,:) = hedm1(:,2);     % PSL GPS heading, deg
    jazz(7,:) = hedm1(:,3);     % PSL GPS angle (pitch), deg
    jazz(8,:) = ship_avg.lat;   % ship Lat
    jazz(9,:) = ship_avg.lon;   % ship Lon
    jazz(10,:) = ship_avg.hed;  % ship Heading
    jazz(11,:) = ship_avg.cog;  % ship COG
    jazz(12,:) = ship_avg.sog;  % ship SOG
    fprintf(fout,['%11.5f ',repmat('%9.3f ',1,11),' \n'],jazz);
    fclose('all');

    %% Bulk met data and surface fluxes for now

    %%% for PISTON, choose TSG2 for now
    Tsea_10 = Tsea2_s_10;
    Ssea_10 = Ssea2_s_10;
    zsea_ship = zsea2_ship;
    
    % tentatively choose the 1st PIR and PSP for first round of fluxes
    lw_dn_10 = lw_dn_1_10;
    sw_dn_10 = sw_dn_1_10;
    lw_dn = lw_dn_1;
    sw_dn = sw_dn_1;
    
    lw_dn_std_10 = interval_std_var(jd_1min, lw_dn, jd_10bin);
    sw_dn_std_10 = interval_std_var(jd_1min, sw_dn, jd_10bin);
    %%% lw and sw are not really offered at time intervals fast enough to
    %%% produce an accurate standard deviation. So, just interpolate the
    %%% 10-min down to 1-min. This has to be redone anyway later when the
    %%% corrections to the radiometers are completed in fixit.m
    lw_dn_std = interp1(jd_10min, lw_dn_std_10, jd_1min);
    sw_dn_std = interp1(jd_10min, sw_dn_std_10, jd_1min);

    
    %%% coare

    coare_psl = coare36vnWarm_et(jd_10min, wspd_10,zu,Ta_10,zt,rh_10,zq,psealevel_10,Tsnk_10,...
                     sw_dn_10,lw_dn_10,lat_10,lon_10, 600,prate_10,zsnk,...
                     Ssea_10, nan, nan, 2, 2, 2);
    coare_s = coare36vnWarm_et(jd_10min, wspd_s_10,zu_ship,Ta_s_10,zt,rh_s_10,zq_ship,psealevel_s_10,Tsea_10,...
                     sw_dn_s_10,lw_dn_s_10,lat_s_10,lon_s_10,600,prate_s_10,zsea_ship,...
                     Ssea_10, nan, nan, 2, 2, 2);
                 
    cfields = {'usr';'tau';'hs';'hl';'hb';'hb_son';'hlwebb';'tsr';'qsr';'zo';'zot';'zoq';...
        'Cd';'Ch';'Ce';'L';'zeta';'dT_skin';'dq_skin';'dz_skin';'U2';'Ta2';'qa2';...
        'rh2';'U2N';'Ta2N';'qa2N';'lw_net';'sw_net';'Le';'rhoa';'UN';'U10';'U10N';'CdN10';'ChN10';'CeN10';...
        'hrain';'qs';'erate';'Ta10';'Ta10N';'qa10';'qa10N';'rh10';'pa10';'rhoa10';'gust';'wc_frac';'Edis';...
        'dT_warm';'dz_warm';'dT_warm_to_skin';'du_warm'};

    cfields_s = strcat(cfields,'_s');

    for i = 1:length(cfields)
        eval([cfields{i} '_10 = coare_psl(:,i);']);
        eval([cfields{i} ' = interp1(jd_10min,' cfields{i} '_10, jd_1min);' ]);
        eval([cfields{i} '_s_10 = coare_s(:,i);']);
        eval([cfields{i} '_s = interp1(jd_10min,' cfields{i} '_s_10, jd_1min);' ]);
    end

    
    %%% in case you want signs to be + heating ocean now, vs later... the code
    %%% right now doesn't flip the signs until save_nc.m
    %%% The bulk sensible, latent, and rain heat fluxes are defined + warming 
    %%% the atmosphere instead by COARE. 
    % flips = {'hs';'hl';'hrain'};
    % for k = 1:length(flips)
    %     eval([flips{k} ' = - ' flips{k} ';']);
    %     eval([flips{k} '_s = - ' flips{k} '_s;']);
    %     eval([flips{k} '_s_10 = - ' flips{k} '_s_10;']);
    %     eval([flips{k} '_10 = - ' flips{k} '_10;']);
    % end

    %% recalculate some fields using COARE output

    % interface
    Tskin = Tsnk + dT_warm_to_skin - dT_skin;
    Tskin_10 = Tsnk_10 + dT_warm_to_skin_10 - dT_skin_10;

    % this is redone in fixit.m once TSGs are checked
    % Tskin_s = Tsea + dT_warm_to_skin_s - dT_skin_s;
    % Tskin_s_10 = Tsea_10 + dT_warm_to_skin_s_10 - dT_skin_s_10;

    % recalculate upwelling radiative fluxes
    lw_up       = lw_net - lw_dn;
    lw_up_10    = lw_net_10 - lw_dn_10;
    sw_up       = sw_net - sw_dn;
    sw_up_10    = sw_net_10 - sw_dn_10;

    % net heat flux now that signs are all consistent: heating into the ocean
    hnet        = sw_net + lw_net - hs - hl - hrain;
    hnet_10     = sw_net_10 + lw_net_10 - hs_10 - hl_10 - hrain_10;
    
    grav = grv(lat);
    lapse=grav/cpa;
    dtheta = Tskin - (Ta + lapse*zt);
    
    grav_10 = grv(lat_10);
    lapse_10=grav_10/cpa;
    dtheta_10 = Tskin_10 - (Ta_10 + lapse_10*zt);
    
%     % potential temperature test of prior approximation in Ta
%     p_tq = psealevel - 0.125*zt;
%     p_tq_10 = psealevel_10 - 0.125*zt;
%     theta0   = (Tskin+C2K).*(1000./psealevel).^(Rgas/cpa);  
%     theta10  = (Ta+C2K).*(1000./p_tq).^(Rgas/cpa);  
%     dtheta   = theta0 - theta10;
%     theta10_10 = (Ta_10+C2K).*(1000./p_tq_10).^(Rgas/cpa);  
%     theta0_10  = (Tskin_10+C2K).*(1000./psealevel_10).^(Rgas/cpa);  
%     dtheta_10  = theta0_10 - theta10_10;

    % tstar from sonic without humidity correction... the bulk value corrects 
    % for it automatically so we have to add it back in (tstr < 0 for SST > Ta)
    tsr_son = tsr + 0.51*(Ta+C2K).*qsr*1e-3;
    tsr_son_10 = tsr_10 + 0.51*(Ta_10+C2K).*qsr_10*1e-3;

    %% BULK FRICTION VELOCITY
        figure; plot(jd_10min, usr_10,'b','linewidth',1); xlabel('Hour (UTC)'); ylabel('u_*(m/s)'); xlim([ddd ddd+1]);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSL Friction Velocity from COARE algorithm',cruise_str,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks'); grid;
        xax = get(gca(gcf),'XTickLabel'); xax(end,:)='24:00'; set(gca(gcf),'XTickLabel',xax(:,1:2));
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA PSL'},'FontSize',14,'FitBoxToText','off','LineStyle','none');
        if prtit
            ppath = fullfile(path_raw_images,'ustar',['Friction_Velocity_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),graphformat]);
            print(graphdevice,ppath);
        end

        close all;
 
    %% save 10-min met and bulk model results to txt file - historical format with WXT data at end

    fpath = fullfile(path_proc_data,'legacy','flux',[cruise,'_flux_10min_',sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt']);
    fout = fopen(fpath,'w');
    jazz = NaN(46,144);
    jazz(1,:) = jd_10min;            % Decimal julian day at beginning of line
    jazz(2,:) = wspd_10;             % PSL true wind speed at zu_PSL, m/s
    jazz(3,:) = wdir_10;             % PSL true wind direction, deg
    jazz(4,:) = Tsnk_10;             % PSL seasnake T, C
    jazz(5,:) = Tsea_10;            % ship theromsalinograph T, C  
    jazz(6,:) = Ssea_10;            % ship theromsalinograph salinity, psu 
    jazz(7,:) = Ta_10;               % PSL air T at zt_PSL, C
    jazz(8,:) = qsnk_10;             % PSL air specific humidity at sea surface, g/kg
    jazz(9,:) = qa_10;               % PSL air specific humidity, g/kg
    jazz(10,:) = sw_dn_10;      % PSL solar flux, w/m^2
    jazz(11,:) = lw_dn_10;      % PSL IR flux, w/m^2
    jazz(12,:) = prate_10;           % PSL org precip rate, mm/hr
    jazz(13,:) = sog_s_10;           % ship sog from gps, m/s
    jazz(14,:) = hed_s_10;           % ship heading from gyrocompass, deg
    jazz(15,:) = rspd_10;           % PSL rel wind speed, m/s
    jazz(16,:) = rdir_10;           % PSL rel wind direction, deg
    jazz(17,:) = lat_10;             % decimal latitude, deg
    jazz(18,:) = lon_10;             % decimal longtude, deg
    jazz(19,:) = ones(size(jd_10min))*zsea_ship;	% Depth of seawater T sensor used in heat flux calc, m
    jazz(20,:) = ship_avg_10.sog_std;        % standard deviation of ship speed, m/s
    jazz(21,:) = tau_10;               % wind stress, coare 3.0, N/m^2
    jazz(22,:) = hs_10;                % sensible heat flux, coare, w/m^2
    jazz(23,:) = hl_10;                % latent heat flux, coare, w/m^2
    jazz(24,:) = hrain_10;                 % rain heat flux, w/m^2
    jazz(25,:) = Ta_s_10;            % ship air temp, C
    jazz(26,:) = qa_s_10;            % ship air specific humidity, g/kg
    jazz(27,:) = wspd_s_10;          % ship true wind speed, m/s
    jazz(28,:) = wdir_s_10;          % ship true wind direction, deg
    jazz(29,:) = sw_dn_s_10;            % ship solar flux, w/m^2
    jazz(30,:) = lw_dn_s_10;            % ship IR flux, w/m^2
    jazz(31,:) = psealevel_s_10;           % ship SLP, mb
    jazz(32,:) = rh_s_10;            % ship RH (%)
    jazz(33,:) = U10_10;                % PSL wspd at 10m from COARE
    jazz(34,:) = Ta10_10;                % PSL Ta at 10m from COARE
    jazz(35,:) = rh10_10;               % PSL rh at 10m from COARE
    jazz(36,:) = U10_s_10;              % ship wspd at 10m from COARE
    jazz(37,:) = Ta10_s_10;              % ship Ta at 10m from COARE
    jazz(38,:) = rh10_s_10;             % ship rh at 10m from COARE
    jazz(39,:) = rspd_wxt_10;       % WXT relative wind speed, m/s
    jazz(40,:) = rdir_wxt_10;       % WXT relative wind direction, deg
    jazz(41,:) = wspd_wxt_10;        % WXT true wind speed, m/s
    jazz(42,:) = wdir_wxt_10;        % WXT true wind direction, deg
    jazz(43,:) = Ta_wxt_10;            % WXT air temp, C
    jazz(44,:) = rh_wxt_10;           % WXT RH, %
    jazz(45,:) = pa_wxt_10;            % WXT sea level pressure, mb
    jazz(46,:) = prate_wxt_10;         % WXT rain rate, mm/hr

    fprintf(fout,['%11.5f ',repmat('%9.3f ',1,47),' \n'],jazz);
    fclose all;
    clear jazz;

    %% radiation models @ 1 and 10 min resolution
    
    % Downward IR 10 min
    if all(isnan(qa_10))
        lw_dn_clr_10 = (0.52+0.13/60*abs(lat_s_10)+(0.082-0.03/60.*abs(lat_s_10)).*sqrt(qa_s_10)).*(5.67e-8*(Ta_s_10+273.15).^4);
    else
        lw_dn_clr_10 = (0.52+0.13/60*abs(lat_10)+(0.082-0.03/60.*abs(lat_10)).*sqrt(qa_10)).*(5.67e-8*(Ta_10+273.15).^4);
    end

    % Downward IR 1 min
    if all(isnan(qa))
        lw_dn_clr = (0.52+0.13/60*abs(lat_s)+(0.082-0.03/60.*abs(lat_s)).*sqrt(qa)).*(5.67e-8*(Ta+273.15).^4);
    else
        lw_dn_clr = (0.52+0.13/60*abs(lat)+(0.082-0.03/60.*abs(lat)).*sqrt(qa)).*(5.67e-8*(Ta+273.15).^4);
    end
    
    % Solar model 10 min
    sw_dn_clr_10 = rs_clear(jd_10min,nanmean(psealevel_10),qa_10,lat_10,lon_10,k1,k2,oz);

    % Solar model 1 min
    sw_dn_clr = rs_clear(jd_1min,nanmean(psealevel),qa,lat,lon,k1,k2,oz);

    disp(['FINISHED evalflux.m ',cruise,' yearday ',sprintf('%03i',ddd),'.']);

    %% save daily 1- and 10-min data to matalb structures f and g

    %%% Do 1-min averages of 1-sec data, then put into daily structure. 
    %%% Do this for fields 3-end, so ignore t and jd since they'll be
    %%% handled separately here in this loop.
    clear a1 a10 
   
    % add these directly. The interval_avg times are centered on 30 sec marks.
    a1.jd    = jd_1min;
    a1.t     = t_1min;
    a10.jd   = jd_10min;
    a10.t    = t_10min;
         
    %%% if there is an issue with t, check rounding.
%     [yyyy, mmm, dd, HH, MM, SS] = datevec(to);
%     t = datenum(yyyy, mmm, dd, HH, roundn(MM,1), 0);

    % ship_vars are same for 1-min and 10-min averages (ship_avg and ship_avg_10)
    ship_avg_no_wind_cmps = rmfield(ship_avg,{'rUn';'rUw';'Un';'Uw';'jd';'t'});
    ship_vars = fields(ship_avg_no_wind_cmps);
    ship_vars_s = strcat(ship_vars,'_s');
    % save 1- and 10-min averages of ship data
    for i = 1:length(ship_vars_s)
        eval([ship_vars_s{i} ' = ship_avg.' ship_vars{i} ';']); 
        eval([ship_vars_s{i} '_10 = ship_avg_10.' ship_vars{i} ';']); 
    end
    
    wxt_vars = {'wspd_wxt';'wdir_wxt';'rspd_wxt';'rdir_wxt';...
            'Ta_wxt';'qa_wxt';'rh_wxt';'rhoa_wxt';'H2O_wxt';...
            'prate_wxt';'paccum_wxt';'psealevel_wxt';'pa_wxt';...
            'rUn_wxt';'rUw_wxt'};
            % 'Un_wxt';'Uw_wxt';
    
    % 1-min vars to save from PSL
    f_psl = {'lat';'lon';'Ta';'Tsonic';'Tsnk';'Tskin';'rh';'qa';'qsnk';...
            'psealevel';'pa';'rhoa';'H2O';'prate';'paccum';...
            'prate_orig';'orgV';'orgV_desp';...
            'sw_dn_1';'sw_dn_2';'lw_dn_1';'lw_dn_2';'sw_dn_clr';'lw_dn_clr';...
            'sw_dn_std';'lw_dn_std';...
            'lw_T_dome_1';'lw_T_case_1';'lw_T_dome_2';'lw_T_case_2';...
            'lw_dn';'lw_up';'sw_dn';'sw_up';...
            'lw_therm_1';'lw_therm_2';'wspd';'wdir';'rspd';'rdir';...
%             'wspd_raw';'wdir_raw';'rspd_raw';'rdir_raw';
            'sog';'cog';'sogE';'sogN';'hed';'hedE';'hedN';...
            'Tskin';'aspir_TRH';'pitch';...
            'licor_AGC';'licor_rhoa_dry';'licor_qa';'licor_qa_dry';'licor_H2O';'licor_H2O_MR';...
            'licor_Tbox';'licor_Pbox';'licor_CO2';'licor_CO2_MR';'licor_rh';...
            'hnet';'lw_net';'sw_net';'dtheta';'tsr_son'};
%              'rUn';'rUw';'Un';'Uw';'rUn_raw';'rUw_raw';'Un_raw';'Uw_raw';...

                
    % add MET, ship, and coare fields (osl + ship) 
    f_fields = vertcat(f_psl, ship_vars_s, MET_all, cfields, cfields_s, wxt_vars);
    
    % save in structure
    for j = 1:length(f_fields)
        eval(['a1.(f_fields{j}) = ' f_fields{j} ';']);
        thisvar = a1.(f_fields{j});
        if length(thisvar) > 50
            eval(['a10.(f_fields{j}) = ' f_fields{j} '_10;']);
        else
           eval(['a10.(f_fields{j}) = ' f_fields{j} ';']);
        end
    end

    % organize structures for saving
    % put all variable names in alphabetical order;
    a1 = orderfields(a1);
    a10 = orderfields(a10);
    
    [a10.year, a10.month, a10.day, a10.hour, a10.minute] = datevec(a10.t);
    [a1.year, a1.month, a1.day, a1.hour, a1.minute] = datevec(a1.t);

    %% 1 min and 10 min File I/O  
    % base file name
    a1_name = [cruise '_1min_data_' sprintf('%02i_%02i_jd%03i',Vdate(2),Vdate(3),ddd)];
    a10_name = [cruise '_10min_data_' sprintf('%02i_%02i_jd%03i',Vdate(2),Vdate(3),ddd)];
    a1_mat = [path_proc_data '/v0_1min/' a1_name '.mat'];
    a10_mat = [path_proc_data '/v0_10min/' a10_name '.mat'];
    save(a1_mat,'a1');    
    save(a10_mat,'a10'); 
    fields_1 = fields(a1);
    fields_10 = fields(a10);
    n1 = length(fields_1);
    n10 = length(fields_10);
    
    % save text file of all these fields; no single value fields are in here
    mat2txt(a1_name,   'n', [path_proc_data '/v0_1min/']  , a1);
    mat2txt(a10_name,  'n', [path_proc_data '/v0_10min/'] , a10);

    make_readme = 0;
    if make_readme == 1
        % assign long variable name or text to each field name
        [l_1,  u_1,  ln_1]   = assign_long_info_met_PISTON(fields_1, n1);
        [l_10, u_10, ln_10]  = assign_long_info_met_PISTON(fields_10, n10);
        for j = 1:n1
            if isempty(l_1{j}) == 1
                disp(['missing info for field ' sprintf('%i',j) ' = ' string(fields_1(j))]);
            end
            if isempty(u_1{j}) == 1
                disp(['missing unit for field ' sprintf('%i',j) ' = ' string(fields_1(j))]);
            end
            if isempty(ln_1{j}) == 1
                disp(['missing long name for field ' sprintf('%i',j) ' = ' string(fields_1(j))]);
            end
        end

        % save readme and header automatically
        readme_full_1     = [path_readme '/' '_readme_1min_' cruise '_met_sea_flux_data_v0.txt'];
        readme_full_10    = [path_readme '/' '_readme_10min_' cruise '_met_sea_flux_data_v0.txt'];
        readme_data_1     = [path_readme '/' '_fields_1min_' cruise '_met_sea_flux_data_v0.txt'];
        readme_data_10    = [path_readme '/' '_fields_10min_' cruise '_met_sea_flux_data_v0.txt'];
        header_data_1     = [path_readme '/' '_header_1min_' cruise '_met_sea_flux_data_v0.txt'];
        header_data_10    = [path_readme '/' '_header_10min_' cruise '_met_sea_flux_data_v0.txt'];

        fileID_data1 = fopen(readme_data_1,'w');
        fileID_data10 = fopen(readme_data_10,'w');
        fileID_head1 = fopen(header_data_1,'w');
        fileID_head10 = fopen(header_data_10,'w');

        for j = 1:n1
            formatSpec = 'i(%i) = %15s %-60s\n';
            fprintf(fileID_data1, formatSpec, j, string(fields_1(j)), string(l_1(j)) );
            fprintf(fileID_head1,'%20s',string(fields_1(j)) );
        end
        fclose(fileID_data1);
        fclose(fileID_head1);

        for j = 1:n10
            formatSpec = 'i(%i) = %15s %-60s\n';
            fprintf(fileID_data10, formatSpec, j, string(fields_10(j)), string(lg(j)) );
            fprintf(fileID_head10,'%20s',string(fields_g_t(j)) );
        end
        fclose(fileID_data10);
        fclose(fileID_head10);

        % concatenate these data fields to the end of readme narrative 
        readme_narrative = [path_readme '/' '_narrative_PISTON_2019_met_sea_flux_v0.txt'];
        command1= ['cat ' readme_narrative ' ' readme_data_1 ' > ' readme_full_1];
        system(command1);

        command10= ['cat ' readme_narrative ' ' readme_data_10 ' > ' readme_full_10];
        system(command10);

    end
    
    % print to screen
    disp('');
    disp(['num vars in  1 min mat file: ' sprintf('%i',n1)]);
    disp(['num vars in 10 min mat file: ' sprintf('%i',n10)]);
    %%% save that in a text file manually and call it the header. 
    
    %% make plots of 1 and 10 min data

    plot_eval(a1, a10, path_prog, path_raw_images);
    %     if ddd == 250
    %         plot_1min_10min(f, a10, path_raw_images);
    %     end
    %     close all; 

end %%% daily for loop


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OUTPUTS:

% a10 and a1 have the same fields:

%                       Cd: [1440×1 double]
%                    CdN10: [1440×1 double]
%                  CdN10_s: [1440×1 double]
%                     Cd_s: [1440×1 double]
%                       Ce: [1440×1 double]
%                    CeN10: [1440×1 double]
%                  CeN10_s: [1440×1 double]
%                     Ce_s: [1440×1 double]
%                       Ch: [1440×1 double]
%                    ChN10: [1440×1 double]
%                  ChN10_s: [1440×1 double]
%                     Ch_s: [1440×1 double]
%                     Edis: [1440×1 double]
%                   Edis_s: [1440×1 double]
%                  Fsea1_s: [1440×1 double]
%                  Fsea2_s: [1440×1 double]
%                      H2O: [1440×1 double]
%                  H2O_wxt: [1440×1 double]
%               ID_Fsea1_s: [1440×1 double]
%               ID_Fsea2_s: [1440×1 double]
%                 ID_tsg_s: [1440×1 double]
%                        L: [1440×1 double]
%                      L_s: [1440×1 double]
%                       Le: [1440×1 double]
%                     Le_s: [1440×1 double]
%                  Ssea1_s: [1440×1 double]
%                  Ssea2_s: [1440×1 double]
%                  Stsg1_s: [1440×1 double]
%                       Ta: [1440×1 double]
%                     Ta10: [1440×1 double]
%                    Ta10N: [1440×1 double]
%                  Ta10N_s: [1440×1 double]
%                   Ta10_s: [1440×1 double]
%                      Ta2: [1440×1 double]
%                     Ta2N: [1440×1 double]
%                   Ta2N_s: [1440×1 double]
%                    Ta2_s: [1440×1 double]
%                     Ta_s: [1440×1 double]
%                   Ta_wxt: [1440×1 double]
%                  Tsea1_s: [1440×1 double]
%                  Tsea2_s: [1440×1 double]
%                    Tskin: [1440×1 double]
%                     Tsnk: [1440×1 double]
%                   Tsonic: [1440×1 double]
%                  Ttsg1_s: [1440×1 double]
%                      U10: [1440×1 double]
%                     U10N: [1440×1 double]
%                   U10N_s: [1440×1 double]
%                    U10_s: [1440×1 double]
%                       U2: [1440×1 double]
%                      U2N: [1440×1 double]
%                    U2N_s: [1440×1 double]
%                     U2_s: [1440×1 double]
%                       UN: [1440×1 double]
%                     UN_s: [1440×1 double]
%                aspir_TRH: [1440×1 double]
%               ched_std_s: [1440×1 double]
%                      cog: [1440×1 double]
%                cog_max_s: [1440×1 double]
%                cog_min_s: [1440×1 double]
%                    cog_s: [1440×1 double]
%                cog_std_s: [1440×1 double]
%                  dT_skin: [1440×1 double]
%                dT_skin_s: [1440×1 double]
%                  dT_warm: [1440×1 double]
%                dT_warm_s: [1440×1 double]
%          dT_warm_to_skin: [1440×1 double]
%        dT_warm_to_skin_s: [1440×1 double]
%                  dq_skin: [1440×1 double]
%                dq_skin_s: [1440×1 double]
%                   dtheta: [1440×1 double]
%                  du_warm: [1440×1 double]
%                du_warm_s: [1440×1 double]
%                  dz_skin: [1440×1 double]
%                dz_skin_s: [1440×1 double]
%                  dz_warm: [1440×1 double]
%                dz_warm_s: [1440×1 double]
%                    erate: [1440×1 double]
%                  erate_s: [1440×1 double]
%                     gust: [1440×1 double]
%                   gust_s: [1440×1 double]
%                       hb: [1440×1 double]
%                     hb_s: [1440×1 double]
%                   hb_son: [1440×1 double]
%                 hb_son_s: [1440×1 double]
%                      hed: [1440×1 double]
%                     hedE: [1440×1 double]
%                   hedE_s: [1440×1 double]
%                     hedN: [1440×1 double]
%                   hedN_s: [1440×1 double]
%                hed_max_s: [1440×1 double]
%                hed_min_s: [1440×1 double]
%                    hed_s: [1440×1 double]
%                hed_std_s: [1440×1 double]
%                       hl: [1440×1 double]
%                     hl_s: [1440×1 double]
%                   hlwebb: [1440×1 double]
%                 hlwebb_s: [1440×1 double]
%                     hnet: [1440×1 double]
%                    hrain: [1440×1 double]
%                  hrain_s: [1440×1 double]
%                       hs: [1440×1 double]
%                     hs_s: [1440×1 double]
%                       jd: [1440×1 double]
%                      lat: [1440×1 double]
%                    lat_s: [1440×1 double]
%                licor_AGC: [1440×1 double]
%                licor_CO2: [1440×1 double]
%             licor_CO2_MR: [1440×1 double]
%                licor_H2O: [1440×1 double]
%             licor_H2O_MR: [1440×1 double]
%               licor_Pbox: [1440×1 double]
%               licor_Tbox: [1440×1 double]
%                 licor_qa: [1440×1 double]
%             licor_qa_dry: [1440×1 double]
%                 licor_rh: [1440×1 double]
%           licor_rhoa_dry: [1440×1 double]
%                      lon: [1440×1 double]
%                    lon_s: [1440×1 double]
%              lw_T_case_1: [1440×1 double]
%              lw_T_case_2: [1440×1 double]
%              lw_T_case_s: [1440×1 double]
%              lw_T_dome_1: [1440×1 double]
%              lw_T_dome_2: [1440×1 double]
%              lw_T_dome_s: [1440×1 double]
%                    lw_dn: [1440×1 double]
%                  lw_dn_1: [1440×1 double]
%                  lw_dn_2: [1440×1 double]
%                lw_dn_clr: [1440×1 double]
%                  lw_dn_s: [1440×1 double]
%                lw_dn_std: [1440×1 double]
%                   lw_net: [1440×1 double]
%                 lw_net_s: [1440×1 double]
%               lw_therm_1: [1440×1 double]
%               lw_therm_2: [1440×1 double]
%               lw_therm_s: [1440×1 double]
%                    lw_up: [1440×1 double]
%                     orgV: [1440×1 double]
%                orgV_desp: [1440×1 double]
%                       pa: [1440×1 double]
%                     pa10: [1440×1 double]
%                   pa10_s: [1440×1 double]
%                     pa_s: [1440×1 double]
%                   pa_wxt: [1440×1 double]
%                   paccum: [1440×1 double]
%                 paccum_s: [1440×1 double]
%               paccum_wxt: [1440×1 double]
%                    pitch: [1440×1 double]
%                    prate: [1440×1 double]
%               prate_orig: [1440×1 double]
%                  prate_s: [1440×1 double]
%                prate_wxt: [1440×1 double]
%                psealevel: [1440×1 double]
%     psealevel_reported_s: [1440×1 double]
%              psealevel_s: [1440×1 double]
%            psealevel_wxt: [1440×1 double]
%                       qa: [1440×1 double]
%                     qa10: [1440×1 double]
%                    qa10N: [1440×1 double]
%                  qa10N_s: [1440×1 double]
%                   qa10_s: [1440×1 double]
%                      qa2: [1440×1 double]
%                     qa2N: [1440×1 double]
%                   qa2N_s: [1440×1 double]
%                    qa2_s: [1440×1 double]
%                     qa_s: [1440×1 double]
%                   qa_wxt: [1440×1 double]
%                       qs: [1440×1 double]
%                     qs_s: [1440×1 double]
%                     qsnk: [1440×1 double]
%                      qsr: [1440×1 double]
%                    qsr_s: [1440×1 double]
%                  rUn_wxt: [1440×1 double]
%                  rUw_wxt: [1440×1 double]
%                     rdir: [1440×1 double]
%                   rdir_s: [1440×1 double]
%               rdir_std_s: [1440×1 double]
%                 rdir_wxt: [1440×1 double]
%                       rh: [1440×1 double]
%                     rh10: [1440×1 double]
%                   rh10_s: [1440×1 double]
%                      rh2: [1440×1 double]
%                    rh2_s: [1440×1 double]
%                     rh_s: [1440×1 double]
%                   rh_wxt: [1440×1 double]
%                     rhoa: [1440×1 double]
%                   rhoa10: [1440×1 double]
%                 rhoa10_s: [1440×1 double]
%                   rhoa_s: [1440×1 double]
%                 rhoa_wxt: [1440×1 double]
%                     rspd: [1440×1 double]
%                   rspd_s: [1440×1 double]
%               rspd_std_s: [1440×1 double]
%                 rspd_wxt: [1440×1 double]
%               shed_std_s: [1440×1 double]
%                      sog: [1440×1 double]
%                     sogE: [1440×1 double]
%                   sogE_s: [1440×1 double]
%                     sogN: [1440×1 double]
%                   sogN_s: [1440×1 double]
%                sog_max_s: [1440×1 double]
%                sog_min_s: [1440×1 double]
%                    sog_s: [1440×1 double]
%                sog_std_s: [1440×1 double]
%                spd_log_s: [1440×1 double]
%            spd_log_std_s: [1440×1 double]
%              spd_log_u_s: [1440×1 double]
%              spd_log_v_s: [1440×1 double]
%                    sw_dn: [1440×1 double]
%                  sw_dn_1: [1440×1 double]
%                  sw_dn_2: [1440×1 double]
%                sw_dn_clr: [1440×1 double]
%                  sw_dn_s: [1440×1 double]
%                sw_dn_std: [1440×1 double]
%                   sw_net: [1440×1 double]
%                 sw_net_s: [1440×1 double]
%                    sw_up: [1440×1 double]
%                        t: [1440×1 double]
%                      tau: [1440×1 double]
%                    tau_s: [1440×1 double]
%                      tsr: [1440×1 double]
%                    tsr_s: [1440×1 double]
%                  tsr_son: [1440×1 double]
%                      usr: [1440×1 double]
%                    usr_s: [1440×1 double]
%                  wc_frac: [1440×1 double]
%                wc_frac_s: [1440×1 double]
%                     wdir: [1440×1 double]
%               wdir_max_s: [1440×1 double]
%               wdir_min_s: [1440×1 double]
%                   wdir_s: [1440×1 double]
%               wdir_std_s: [1440×1 double]
%                 wdir_wxt: [1440×1 double]
%                     wspd: [1440×1 double]
%               wspd_max_s: [1440×1 double]
%               wspd_min_s: [1440×1 double]
%                   wspd_s: [1440×1 double]
%               wspd_std_s: [1440×1 double]
%                 wspd_wxt: [1440×1 double]
%                     zeta: [1440×1 double]
%                   zeta_s: [1440×1 double]
%                       zo: [1440×1 double]
%                     zo_s: [1440×1 double]
%                      zoq: [1440×1 double]
%                    zoq_s: [1440×1 double]
%                      zot: [1440×1 double]
%                    zot_s: [1440×1 double]
%                     year: [1440×1 double]
%                    month: [1440×1 double]
%                      day: [1440×1 double]
%                     hour: [1440×1 double]
%                   minute: [1440×1 double]
