% a2_motcorr_flux_FATIMA_2022.m

%{
    2020/07  BWB
    
    Driver script for the FATIMA_2022 flux calculations on MV Atlantic Condor
    
    Fluxes are processed in 10-min segments.
    
    Set run parameters in the first code block below as desired and run.
    
    V2: - With updated and vectorized code, May 2015.
        - 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, 2018
    V5 : Using Intertial Labs MRU for motion, reformatted output da table, July 2022
%}

% Output da table format:  not all variables are computed on this cruise
%{
                        TIMESTAMP VARIABLES
    1     jd              decimal_day   decimal day of year timestamp
    2     dnum            decimal_day   matlab datenum timestamp
    3     year            none          year (cols 3:8 are date;time vector)
    4     month           none          month
    5     day             none          day of month
    6     hour            none          hour of day; 0;23
    7     minute          none          minute of hour
    8     second          none          second (zeros for this file)
                        COARE 3.5 SHIP (SCS) INPUT FIELDS, cols 9:21
    9     wspd_ship       m/s           true wind speed; ship; m/s
    10    zu_ship         m             wspd measurement height; ship; m
    11    ta_ship         degC          air temperature; ship; degC
    12    zt_ship         m             ta measurement height; ship; m
    13    rh_ship         percent       relative humidity; ship; %
    14    zq_ship         m             rh measurement height; ship; m
    15    Pmb_ship        mb            surface air pressure; ship; mb
    16    sst_ship        degC          SST ship TSG; ¬degC
    17    rs_ship         W/m2          solar radiative flux; ship; W/m2
    18    rl_ship         W/m2          IR radiative flux; ship; W/m2
    19    lat_ship        degN          latitude; ship; degN
    20    zi_ship         m             MBL depth; NaN for default 600m
    21    rain_ship       mm/hr         rain rate; ship; mm/hr
                        COARE 3.5 NOAA PSL INPUT FIELDS, cols 22:34
    22    wspd_psl        m/s           true wind speed; NOAA PSL sonic; m/s
    23    zu_psl          m             wspd measurement height; m
    24    ta_psl          degC          air temperature; NOAA PSL Vaisala HMT; degC
    25    zt_psl          m             ta measurement height; m
    26    rh_psl          percent       relative humidity; NOAA PSL Vaisala HMT; %
    27    zq_psl          m             rh measurement height; m
    28    Pmb_psl         mb            surface air pressure; NOAA PSL Vaisala PTB; mb
    29    sst_psl         degC          SST; seasnake; ¬degC
    30    rs_psl          W/m2          solar radiative flux; NOAA PSL psp; W/m2
    31    rl_psl          W/m2          IR radiative flux; NOAA PSL pir; W/m2
    32    lat_psl         degN          latitude; NOAA PSL gps; degN
    33    zi_psl          m             MBL depth; NaN for default 600m
    34    rain_psl        mm/hr         rain rate; org; NOAA PSL; mm/hr
                        OTHER WIND VARIABLES
    35    wdir_ship       deg           ship mean true wind direction; deg
    36    wdir_ship_std   deg           std dev ship true wind direction; deg
    37    rwspd_ship      m/s           ship rel wind speed
    38    rwspd_ship_std  m/s           ship rwspd std dev
    39    rwdir_ship      deg           ship rel wind dir; +/; 180 deg
    40    rwdir_ship_std  deg           ship rwdir std dev
    41    wdir_psl        deg           true wind direction from N;E components; NOAA PSL sonic; deg
    42    rwspd_psl       m/s           streamwise relative wind speed from NOAA PSL sonic; m/s
    43    rwspd_psl_std   m/s           std dev rwspd from NOAA PSL sonic; m/s
    44    rwdir_psl       deg           relative wind direction from NOAA PSL sonic; deg
    45    rwdir_psl_std   deg           std dev rwdir from NOAA PSL sonic; deg
    46    Un              m/s           N component; corrected true windspeed; NOAA PSL sonic; m/s
    47    Ue              m/s           E component; corrected true windspeed; NOAA PSL sonic; m/s
    48    Ubar            m/s           streamwise true wspd; NOAA PSL sonic; no flow dist. correction; m/s
                        NOAA PSL WXT VARIABLES
    49    rwdir_wxt       deg           rel wind direction ;0;360 deg from bow; NOAA PSL WXT
    50    rwspd_wxt       m/s           wind speed; m/s; NOAA PSL WXT
    51    ta_wxt          degC          air temperature; degC; NOAA PSL WXT
    52    rh_wxt          percent       RH; %; NOAA PSL WXT
    53    Pmb_wxt         mb            sea level air pressure; mb or hPa; NOAA PSL WXT
    54    rain_wxt        mm/hr         rain rate; mm/hr; NOAA PSL WXT
    55    wspd_wxt        m/s           true wind speed; m/s; NOAA PSL WXT
    56    wdir_wxt        deg           true wind dir; m/s; NOAA PSL WXT
    57    Urel_wxt        m/s           U relative; m/s; NOAA PSL WXT
    58    Vrel_wxt        m/s           V relative; m/s; NOAA PSL WXT
                        OTHER MET/OCEAN VARIABLES
    59    qa_psl          g/kg          air specific humidity; NOAA PSL sensors; g/kg
    60    qa_ship         g/kg          air specific humidity; ship sensors; g/kg
    61    qa_wxt          g/kg          WXT specific humidity; g/kg
    62    qs_psl          g/kg          sea surface saturated specific humidity; seasnake SST; g/kg
    63    qs_ship         g/kg          sea surface saturated specific humidity; ship SST; g/kg
    64    qs_rosr         g/kg          sea surface saturated specific humidity; ROSR SST; g/kg
    65    rhoa_psl        kg/m3         moist air density from NOAA PSL sensors; kg/m3
    66    sal_ship        ppth          surface seawater salinity; ship TSG; ppth
    67    sst_rosr        degC          ROSR radiometric SST; degC
    68    aspir_on        none          Vaisala HMT ventilation on/off parameter
    69    org_carrier     volts         ORG rainrate QA parameter; volts
                       NAVIGATION VARIABLES
    70    sog_ship        m/s           mean sog; ship gps; m/s
    71    sog_ship_std    m/s           std sog; ship gps; m/s
    72    sog_ship_min    m/s           minimum sog; ship gps; m/s
    73    sog_ship_max    m/s           maximum sog; ship gps; m/s
    74    cog_ship        deg           mean cog; ship gps; deg
    75    cog_ship_std    deg           std cog; ship gps; deg
    76    sog_psl         m/s           mean sog; m/s; NOAA PSL gps; deg
    77    sog_psl_std     m/s           std sog; NOAA PSL gps; m/s
    78    sog_psl_min     m/s           minimum sog; NOAA PSL gps; deg
    79    sog_psl_max     m/s           maximum sog; NOAA PSL gps; deg
    80    cog_psl         deg           mean cog; deg; NOAA PSL gps; deg
    81    cog_psl_std     deg           std cog; NOAA PSL gps; deg
    82    lrg             deg           heading; ship gyro; deg
    83    lrg_std         deg           std dev heading; ship gyro; deg
    84    hed_psl         deg           heading; NOAA PSL gps; deg
    85    hed_psl_std     deg           std dev heading; NOAA PSL gps; deg
    86    lon_ship        degE          longitude; ship gps
    87    lon_psl         degE          latitiude; NOAA PSL gps
    88    odec            m/s           odec (ship speed through water); m/s
    89    odec_std        m/s           stdev odec; m/s
    90    odec_min        m/s           minimum odec; m/s
    91    odec_max        m/s           maximum odec; m/s
                        COARE 3.5 ship data output fields
    92    usb_ship        m/s           friction velocity; ship; m/s
    93    taub_ship       N/m2          surface wind stress; ship; N/m2
    94    hsb_ship        W/m2          sensible heat flux; ship; W/m2
    95    hlb_ship        W/m2          latent heat flux; ship; W/m2
    96    hlwebb_ship     W/m2          hlb 'webb' correction; W/m2
    97    tsb_ship        degK          t;star; ship; ¬degC
    98    qsb_ship        g/kg          q;star; ship; g/kg
    99    zo_ship         m             surface roughness length; ship; m
    100   zot_ship        m             surface temperature roughness length; ship; m
    101   zoq_ship        m             surface moisture roughness length; ship; m
    102   Cd_ship         none          drag coefficient; ship
    103   Ch_ship         none          sensible heat transfer coefficient; ship
    104   Ce_ship         none          water vapor transfer coerfficient; ship
    105   L_ship          m             Obukhov length; ship; m
    106   zet_ship        none          z/L; ship @ zu_ship
    107   dter_ship       degC          cool skin temperature depression; ship; degC
    108   dqer_ship       degC          cool skin saturation humidity depression; g/kg
    109   tkt_ship        m             cool skin thickness; ship; m
    110   RF_ship         W/m2          rain heat flux; ship; W/m2
    111   Cdn_10_ship     none          10m neutral drag coefficient; ship
    112   Chn_10_ship     none          10m neutral heat transfer coefficient; ship
    113   Cen_10_ship     none          10m neutral moisture transfer coefficient; ship
                        COARE 3.5 NOAA PSL data output fields
    114   usb_psl         m/s           friction velocity; NOAA PSL; m/s
    115   taub_psl        N/m2          surface wind stress; NOAA PSL; N/m2
    116   hsb_psl         W/m2          sensible heat flux; NOAA PSL; W/m2
    117   hlb_psl         W/m2          latent heat flux; NOAA PSL; W/m2
    118   hlwebb_psl      W/m2          hlb 'webb' correction; NOAA PSL; W/m2
    119   tsb_psl         degK          t;star; NOAA PSL; degK
    120   qsb_psl         g/kg          q;star; NOAA PSL; g/kg
    121   zo_psl          m             surface roughness length; NOAA PSL; m
    122   zot_psl         m             surface temperature roughness length; NOAA PSL; m
    123   zoq_psl         m             surface moisture roughness length; NOAA PSL; m
    124   Cd_psl          none          drag coefficient; NOAA PSL
    125   Ch_psl          none          sensible heat transfer coefficient; NOAA PSL
    126   Ce_psl          none          water vapor transfer coerfficient; NOAA PSL
    127   L_psl           m             Obukhov length; NOAA PSL; m
    128   zet_psl         none          z/L; NOAA PSL @ zu_psl
    129   dter_psl        degC          cool skin temperature depression; NOAA PSL; degC
    130   dqer_psl        degC          cool skin humidity depression; NOAA PSL; degC
    131   tkt_psl         m             cool skin thickness; NOAA PSL; m
    132   RF_psl          W/m2          rain heat flux; NOAA PSL; W/m2
    133   Cdn_10_psl      none          10m neutral drag coefficient; NOAA PSL
    134   Chn_10_psl      none          10m neutral heat transfer coefficient; NOAA PSL
    135   Cen_10_psl      none          10m neutral moisture transfer coefficient; NOAA PSL
                        FLUX VARIABLES
    136   wu_cov          m2/s2         wu covariance; Edson correction; m2/s2
    137   wu_csp          m2/s2         wu cospectrum integral; Edson correction; m2/s2
    138   wv_cov          m2/s2         wv covariance; Edson correction; m2/s2
    139   wv_csp          m2/s2         wv cospectrum integral; Edson correction; m2/s2
    140   wTs_cov         degC_m/s      wTsonic covariance; Edson correction; ¬degC m/s
    141   wTs_csp         degC_m/s      wTsonic cospectrum integral; Edson correction; m2/s2
    142   wT_cov          degC_m/s      wT covariance; ¬degC m/s
    143   wT_csp          degC_m/s      wT cospectrum integral; Edson correction; m2/s2
    144   wq_cov          g/kg_m/s      wq covariance; Licor; g/kg m/s
    145   wq_csp          g/kg_m/s      wq cospectrum integral; Licor; g/kg m/s
    146   hs_cov          W/m2          sensible heat flux from covariance; W/m2
    147   hs_csp          W/m2          sensible heat flux from cospectrum; W/m2
    148   hl_cov          W/m2          latent heat flux from covariance; uncorrected; W/m2
    149   hl_csp          W/m2          latent heat flux from cospectrum; uncorrected; W/m2
    150   hlc_cov         W/m2          webb corrected hl_cov; W/m2
    151   hlc_csp         W/m2          webb corrected hl_csp; W/m2
    152   Cu2a            none          Cu structure func parameter; method a
    153   Cu2b            none          Cu structure func parameter; method b
    154   Cw2a            none          Cw structure func parameter; method a
    155   Cw2b            none          Cw structure func parameter; method b
    156   Ct2a            none          Ct structure func parameter; method a
    157   Ct2b            none          Ct structure func parameter; method b
    158   Cq2a            none          Cq structure func parameter; method a
    159   Cq2b            none          Cq structure func parameter; method b
    160   usid_a          m/s           ustar from inertial diss.; method a
    161   usid_b          m/s           ustar from inertial diss.; method b
    162   tsid_a          degK          tstar from inertial diss.; method a
    163   tsid_b          degK          tstar from inertial diss.; method b
    164   qsid_a          g/kg          qstar from inertial diss.; method a
    165   qsid_b          g/kg          qstar from inertial diss.; method b
    166   Lid_a           m             L from inertial diss.; method a
    167   Lid_b           m             L from inertial diss.; method b
                        TURBULENCE STATISTICS AND QA/QC PARAMETERS
    168   uvar            m2/s2         u variance; Edson correction; m2/s2
    169   vvar            m2/s2         v variance; Edson correction; m2/s2
    170   wvar            m2/s2         w variance; Edson correction; m2/s2
    171   tvar            degC^2        Tsonic variance; degC^2
    172   qvar            g2/kg2        Licor q variance; g/kg
    173   azm             deg           azimuth angle from first rotation; radians
    174   tilt            deg           tilt angle from first rotation; radians
    175   fsonic          Hz            sonic anemometer sampling frequency; Hz
    176   missing_son     none          missing sonic wind data lines per 10 min interval @ fsonic
    177   badPnts_son     none          #NaNs in sonic data per 10 min interval @ fsonic
    178   flicor          Hz            Licor 7500 sampling frequency; Hz
    179   missing_lic     none          missing licor data lines per 10 min interval @ flicor
    180   fmotion         Hz            motion system sampling frequency; Hz
    181   missing_mot     none          missing motion data lines per 10 min interval @ fmotion
    182   fheading        Hz            Hemisphere heading system sampling frequency; Hz
    183   missing_hed     none          missing heading data lines per 10 min interval @ fheading
    184   uplat_std       m/s           U platform velocity std dev; Edson correction; m2/s2
    185   vplat_std       m/s           V platform velocity std dev; Edson correction; m2/s2
    186   wplat_std       m/s           W platform velocity std dev; Edson correction; m2/s2
    187   Ts_noise        degC^2        Sonic temperature noise variance; degC^2
                    OTHER LICOR 7500 QA/QC VARIABLES - 'good' is ADC<65
    188   h2o_lic_all     mmol/m3       mean H2O; all data; mmol/m3
    189   h2o_lic_good    mmol/m3       mean H2O; good data; mmol/m3
    190   co2_lic_all     ppm           mean CO2; all data; ppm mole fraction
    191   co2_lic_good    ppm           mean CO2; good data; ppm mole fractio
    192   Pmb_lic         mb            pressure from Licor box; mb
    193   T_lic           degC          temperature from Licor box; degC
    194   h2o_std_all     mmol/m3       std dev H2O; all data; mmol/m3
    195   h2o_std_good    mmol/m3       std dev H2O; good data; mmol/m3
    196   co2_std_all     ppm           std dev CO2; all data; ppm mole fractio
    197   co2_std_good    ppm           std dev CO2; good data; ppm mole fractio
    198   AGC_lic         none          Licor diagnostic; 50-100; 50 = best
    199   q_lic           g/kg          Licor specific humidity
%}

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

% set start and end dates for the flux calcs
dnum_start = datenum(2022,7,4);
dnum_end = datenum(2022,8,1);

read_raw = false;  % set true to force read of raw 10Hz data files

% cruise ID and paths
cruise = 'FATIMA_2022';  % cruise name
ship = 'Atl_Condor';     % research vessel
plotit = true;           % create plots and save to .png format

% system specific path defs
data_drive = '/Users/byronb/Raw_Data_Files.noindex';

% matlab script path
path_prog = fullfile(data_drive,cruise,'_matlab');
cd(path_prog);
addpath(genpath(path_prog));

% path to excel spreadsheet of da variable names and attributes
path_da_params = fullfile(path_prog,'da_vars.xlsx');

% lat/lon limits and temperature adjustments
Lonmin = 292; Lonmax = 313; Latmin = 35; Latmax = 52;
PosLims = [Lonmin, Lonmax, Latmin, Latmax];

% no temp adjustments this run
tsea_adj = 0; td1_adj = 0; tc1_adj = 0;
td2_adj = 0; tc2_adj = 0; ta_wxt_adj = 0;
adj = [tsea_adj,td1_adj,tc1_adj,td2_adj,tc2_adj,ta_wxt_adj];

% longest gap in met data to fill with est values using fillgaps()
met_gap_lim = 5; % met0 data gaps of 5 min or less will be filled

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

fsons = 20; % sonic and motion data rate, Hz

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

% build paths relative to data_drive
% these 3 paths should be preexisting
path_raw_data_flux = fullfile(data_drive,cruise,'raw');
path_proc_data_flux = fullfile(data_drive,cruise,'processed');
path_proc_images_flux = fullfile(data_drive,cruise,'processed_images');
rosr_path = fullfile(path_raw_data_flux,'rosr','sst_rosr_all.mat');

% 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,'true_wind'));
% subdirectory for saving corrected wind timeseries plots
path_windPlots = fullfile(path_proc_images_flux,['wind_plots_',dcname]);
mkdir(path_windPlots);
mkdir(fullfile(path_windPlots,'motcorr_segment'));
mkdir(fullfile(path_windPlots,'uvw'));
mkdir(fullfile(path_windPlots,'Tson'));
% subdirectory for saving corrected wind spectra plots
path_spectraPlots = fullfile(path_proc_images_flux,['spectra_',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 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

% flux data sample rate
fsam = 10;            % frequency (Hz) for all fast data
Npts = fsam*600;      % data points in 10 min window @ frequency fsam
% raw data acquisition rates
fsonic = 20; flicor = 10; fmotion = 20; fheading = 2;

% as measured in July 2022, Halifax installation
zu = 18.6;      % sonic wind speed measurement height (m)
z_wxt = 18;     % wxt T/RH measurement height (m)
z_PTB = 16.25;  % PTB pressure sensor height (m)

% sonic-MotionPak displacement vector [Rx Ry Rz]
% sonic is 0.914 m above and 0.15m forward of MRU
sens_disp = [0.15,0,0.914]; %[X,Y,Z]

% constants
d2r = pi/180;  % degrees to radians
r2d = 180/pi;  % radians to degrees
C2K = 273.15;  % °C to °K
Rgas = 287.1;              % Pa m-3 kg-1 K-1 for dry air
Rgas_universal = 8.314472; % Pa m3 mol-1 K-1
Mw = 18.01528; % molar mass of H20 g/mol
Md = 28.964;   % molar mass dry air g/mol
epsilon = Mw/Md; % mass of water to mass of dry air = 0.622

% Load variable names, units and descriptions for output da table.
% Variable names and attributes are maintained in an excel spreadsheet
% Edit the spreadsheet to update variable names, etc.
da_params = readcell(path_da_params);
da_vars = da_params(:,2);
da_units = da_params(:,3);
da_descript = da_params(:,4);

%% ----------Main loop----------
for dxx = dnum_start:dnum_end  % iterate over days
    % Date variables for this day
    Vdate = datevec(dxx); Vdate = Vdate(1:3);
    yyyy = Vdate(1); YYYY = sprintf('%4i',yyyy); YY = YYYY(3:4);
    mm = Vdate(2); MM = sprintf('%02i',mm);
    dd = Vdate(3); DD = sprintf('%02i',dd);
    ddd = md2yd_vn(yyyy, mm, dd);
    DDD = sprintf('%03i',ddd);

    % 10-min timestamp for this day
    dnum_10min = dxx:10/1440:dxx+1-10/1440;  dnum_10min = dnum_10min(:);

    % empty da table for 10-min results
    da = array2table(NaN(144,length(da_vars)),'VariableNames',da_vars);
    da.Properties.VariableUnits = da_units;
    da.Properties.VariableDescriptions = da_descript;
    da.dnum = dnum_10min;
    da.jd = da.dnum - datenum(yyyy-1,12,31);
    da.year(:) = yyyy; da.month(:) = mm; da.day(:) = dd;
    da.second(:) = 0;

    % set constants in da table for this day
    % using wxt for psl T/RH - no HMT sensor this cruise
    da.fsonic(:) = fsonic; da.flicor(:) = flicor;
    da.fmotion(:) = fmotion; da.fheading(:) = fheading;
    da.zu_psl(:) = zu; da.zt_psl(:) = z_wxt;
    da.zq_psl(:) = z_wxt;
    
    % array for spectra
    sp = NaN(144,481); sp(:,1) = dnum_10min;
    
    if decorr
        % table for W and Tsonic decorrelation coefficients
        mu_vars = {'dnum','muu_accx0','muu_uplat0','muu_accx','muu_uplat',...
            'muv_accy0','muv_vplat0','muv_accy','muv_vplat',...
            'muw_accz0','muw_wplat0','muw_accz','muw_wplat',...
            'muTs_accx0','muTs_accy0','muTs_accz0',...
            'muTs_uplat0','muTs_vplat0','muTs_wplat0',...
            'muTs_accx','muTs_accy','muTs_accz',...
            'muTs_uplat','muTs_vplat','muTs_wplat'};
        mu_decorr = array2table(NaN(144,25),'VariableNames',mu_vars);
        mu_decorr.dnum = dnum_10min;

        % table for q Licor decorrelation coefficients
        muq_vars = {'dnum','muq_accx0','muq_accy0','muq_accz0','muq_uplat0','muq_vplat0',...
            'muq_wplat0','muq_accx','muq_accy','muq_accz','muq_uplat','muq_vplat',...
            'muq_wplat'};
        mu_q_decorr = array2table(NaN(144,13),'VariableNames',muq_vars);
        mu_q_decorr.dnum = dnum_10min;
    end
    
    % table for hourly lag output
    lags = array2table(NaN(24,2),'VariableNames',{'dnum','lag_points'});
    lags.dnum = (dxx:1/24:dxx+1-1/24)';

    %%%%%%%%%% loop over 24 hours of data %%%%%%%%%%

    for hh = 0:23   % iterate over hours
        %% Set timestamp variables
        HH = sprintf('%02i',hh);
        disp(['%%%% DAY ',DDD,' HOUR ',HH,' %%%%']);
        % start of hour
        start = datenum(yyyy,mm,dd,hh,0,0);
        % ref 10 Hz timestamp
        dnum_10Hz = start:1/864000:(start + 1/24 - 1/864000);
        dnum_10Hz = dnum_10Hz(:);
        % 10-min bin edges for averaging and computing flux segments
        dnum_10min_bin = (start:10/1440:start+1/24);
        dnum_10min_bin = dnum_10min_bin(:);
        % select indices in daily tables for this hour's results
        iii = find(abs(da.dnum-start)<1/86400);
        iii = iii:1:iii+5; % six indices for this hour

        da.hour(iii) = hour(datetime(start,'ConvertFrom','datenum'));
        da.minute(iii) = [0:10:50]';

        %% read raw data files...
        
        % read 1-min PSL met table with the following variables:
        %{
            1   dnum            matlab datenum timestamp
            2   Tsea            sea snake surface SST, °C
            3   psp1            psp1 radiation, W/m2
            4   pir1            pir1 radiation, W/m2
            5   psp2            psp2 radiation, W/m2
            6   pir2            pir2 radiation, W/m2
            7   Tc1             pir1 case temperature, °C
            8   Td1             pir1 dome temperature, °C
            9   Tc2             pir2 case temperature, °C
            10  Td2             pir2 dome temperature, °C
            11  Pmb             sea level air pressure, mb(hPa), Vaisala PTB330
            12  rwdir_wxt       rel wind direction ,0-360 deg from bow
            13  rwspd_wxt       wind speed, m/s
            14  Tair_wxt        air temperature, °C
            15  RH_wxt          RH, %
            16  Pmb_wxt         sea level air pressure, mb(hPa)
            17  rain_wxt        rain rate, mm/hr
            18  U_wxt           rel U wind component, WXT, m/s
            19  V_wxt           rel V wind component, WXT, m/s
            20  sst_rosr        radiometric sst, °C
        %}
        dfl1 = fullfile(path_raw_data_flux,'MeanMet',[YY,DDD],...
                        ['met0',YY,DDD,HH,'_raw.txt']);
        dfl2 = fullfile(path_raw_data_flux,'WXT520',[YY,DDD],...
                        ['wxt_parsed_',YY,DDD,HH,'_raw.txt']);
        met = read_met_FATIMA_2022(dfl1,dfl2,rosr_path,z_wxt,z_PTB,...
                                   adj,met_gap_lim);
        % sea snake SST bad after 17 Jul 06:00
        ii = met.dnum >= datenum(2022,7,17,6,0,0);
        met.Tsea(ii) = NaN;

        % read 10 Hz mot table with following variables:
        %{
            1   dnum     matlab datenum timestamp
            2   yyyy     year
            3   mm       month
            4   dd       day of month
            5   hr       hour of day
            6   min      minute
            7   sec      decimal seconds
            8   hed      heading, deg
            9   pitch    pitch angle, deg
            10  roll     roll angle, deg
            11  ratex    x axis rotational rate, rad/s
            12  ratey    y axis rotational rate, rad/s
            13  ratez    z axis rotational rate, rad/s
            14  accx     x axis acceleration, m/s2
            15  accy     y axis acceleration, m/s2
            16  accz     z axis acceleration, m/s2
            17  magx     x axis magnetic field strength, nT
            18  magy     y axis magnetic field strength, nT
            19  magz     z axis magnetic field strength, nT
            20  lat      latitude N, deg
            21  lon      longitude E, deg
            22  heave    vertical heave, m
            23  uplat    forward(N) velocity, m/s, includes ship speed...
            24  vplat    to-port(W) velocity, m/s, includes ship speed...
            25  wplat    vertical velocity, m/s

            PSL standard right-handed axis/sign conventions:
             +x velocity is forward
             +y velocity is to port
             +z velocity is up
             +rotation about the x-axis is port up (phi = roll)
             +rotation about the y-axis is bow down (theta = pitch)
             +rotation about the z-axis is bow rotation to port (psi = yaw)
             vplat is positive W in earth coords, doesn't include ship speed...
             uplat is positive N in earth coords, doesn't include ship speed...
        %}
        dfl3 = fullfile(path_raw_data_flux,'MRUdata',[YY,DDD],...
                        [YYYY,DDD,HH,'00MRU_P.txt']);
        [mot,dT_hist,dT_edges,missing_mot] = read_motion_mru_FATIMA_2022(dfl3,read_raw);
        da.missing_mot(iii) = missing_mot;

        % read sonic table with the following variables:
        %{
            1   dnum        matlab datenum 10 Hz timestamp
            2   u           bow-stern axis wind velocity, m/s
            3   v           port-starboard axis wind velocity, m/s
            4   w           z axis wind velocity, m/s
            5   Tsonic      sonic temperature, C
        %}
        if start< datenum(2022,7,18,0,0,0)
            rotationsonic = true;  % leg 1 sonic
        else
            rotationsonic = false; % leg 2 sonic - no rotation
        end
        dfl4 = fullfile(path_raw_data_flux,'SonicGillR3',[YY,DDD],...
                        ['son0',YY,DDD,HH,'_raw.txt']);

        [son,~,~,badPnts,missing_son] = read_sonic_FATIMA_2022(dfl4,rotationsonic,read_raw,fsons);
        xxx = missing_son>100 | missing_son>100;  % flag for bad 10-min sonic segments
        da.missing_son(iii) = missing_son;
        da.badPnts_son(iii) = badPnts;

        % read gps table - PYTHON-processed $GPRMC lines from Hemisphere gps
        %{
            1   tref        10 Hz timestamp
            2   cog         cog, deg
            3   sog         sog, m/s
            4   Ngps        filtered N component of ship speed, m/s
            5   Egps        filtered E component of ship speed, m/s
            6   lat         latitude N, deg
            7   lon         longitude E, deg
        %}
        dfl5 = fullfile(path_raw_data_flux,'GPS_Hemisphere',[YY,DDD],...
                        ['gprm' YY,DDD,HH,'_raw.txt']);
        dfl5b = fullfile(path_raw_data_flux,'GPS_MRUP',[YY,DDD],...
                        ['gprm' YY,DDD,HH,'_raw.txt']);
        gps = read_gprm_FATIMA_2022(dfl5,PosLims);
        if sum(isnan(gps.lat)) > 0.05*3600
            gps = read_gprm_FATIMA_2022(dfl5b,PosLims);
        end

        % read Hemisphere heading and roll angle table
        %{
            1   dnum        matlab datenum timestamp
            2   hed         heading, deg
            3   roll        roll angle, deg
        %}
        dfl6 = fullfile(path_raw_data_flux,'HED_Hemisphere',[YY,DDD],...
                        ['hed0' YY,DDD,HH,'_raw.txt']);
        [hed,~,~,missing_hed] = read_hed_roll_FATIMA_2022(dfl6);
        da.missing_hed(iii) = missing_hed;

        % read licor table
        %{
            1   dnum        10 Hz matlab datenum timestamp
            2   h20         H2O, mmol/m3
            3   co2         CO2, approx mmol/m3
            5   Ta_lic      Licor box temperature, °C
            6   Pmb_lic     Licor box pressure, mb
            7   AGC         Licor diagnostic, 0-100
        %}
        dfl7 = fullfile(path_raw_data_flux,'GasLicor',[YY,DDD],...
                        ['lic0' YY,DDD,HH,'_raw.txt']);
        [lic,~,~,missing_lic] = read_licor_FATIMA_2022(dfl7,read_raw);
        da.missing_lic(iii) = missing_lic;

        %% interpolate 10 Hz heading
        hed0 = unwrap(hed.hed*d2r); % in radians
        hed1 = fillgaps(hed0,1800); % fill in data gaps
        heading = interp1(hed.dnum,hed1,dnum_10Hz,'linear','extrap'); % 10 Hz heading
        shead = sin(heading);       % 10Hz vector components
        chead = cos(heading);
        
        %% run motcorr if NaNs are less than 5% of total number of wind or
        % motion data points
        if sum(isnan(mot.accx))<0.1*36000 && sum(isnan(son.u)) < 0.1*36000 && ...
           sum(isnan(hed.hed))<0.05*36000
            % E98 motion corrections, with decorr if desired
            % use motcorr4, decorr in 10 min segments
            % prepare input arrays for motcorr4()
            motion(:,1) = mot{:,1};
            motion(:,2:4) = fillgaps(mot{:,14:16},1800);    % fill gaps in accs, m/s2
            motion(:,5:7) = fillgaps(mot{:,11:13},1800);    % fill gaps in rates, rad/s
            motion(:,8:10) = fillgaps(mot{:,23:25},1800);   % fill gaps in mru uvwplat, m/s
            sonic(:,1) = son{:,1};
            sonic(:,2:5) = fillgaps(son{:,2:5},1800);       % fix gaps in uvwt

            % remove filtered gps speed components from mru uplat & vplat
            motion(:,8) = motion(:,8) + gps.Ngps;
            motion(:,9) = motion(:,9) + gps.Egps;

            % run motion corrections
            [uvw,T_10Hz,accplat,uvwplat,xyzplat,euler,plat_rate,lagPnts,mu,uvw_raw,...
                uvw_motcorr,f1,f2,f3] = motcorr4(sonic,motion,heading,sens_disp,fsam,decorr);

            % save motcorr wind plots every 6 hours
            if ~all(isnan(uvw(:,1))) && plotit && (hh==0 || hh==6 || hh==12 || hh==18)
                ppath = fullfile(path_windPlots,'motcorr_segment',...
                    ['uvw_motcorr_',DDD,'_',HH,'_',dcname,graphformat]);
                print(f1,graphdevice,ppath);
                ppath = fullfile(path_windPlots,'uvw',['uvw_',DDD,'_',HH,'_',dcname,graphformat]);
                print(f2,graphdevice,ppath);
                ppath = fullfile(path_windPlots,'Tson',['Tson_',DDD,'_',HH,'_',dcname,graphformat]);
                print(f3,graphdevice,ppath);
                close all;
            else
                close all;
            end

            % reformat output of motcorr4 to column arrays
            uvw = uvw';                    % uvw rel winds - Edson motcorr and decorr
            T_10Hz = T_10Hz';              % sonic temperature
            accplat = accplat';            % platform acc,vel and displacement in earth frame
            uvwplat = uvwplat';
            xyzplat = xyzplat';
            euler = euler';                % euler angles
            plat_rate = plat_rate';        % platform rotational rates
            uvw_raw = uvw_raw';            % sonic uvw rotated to earth frame
            uvw_motcorr = uvw_motcorr';    % uvw with Edson correction only, no decorr

            lll = find(abs(lags.dnum-start)<1/86400);
            lags.lag_points (lll) = lagPnts; % save lag adjustment

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

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

            % save azm and tilt
            da.azm(iii) = azm*r2d; da.tilt(iii) = tilt*r2d;

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

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

            % add timestamp and sonicT to uvw_str table
            % these are motion corrected winds for flux calcs.
            uvw_str = array2table([dnum_10Hz, uvw_str, T_10Hz],'VariableNames',...
                      {'dnum','ustr','vstr','wstr','Tson'});
            clear uvw_rot azm tilt azm2
        else
            % too many NaNs, so return empty motcorr output arrays
            uvw_str = array2table([dnum_10Hz,NaN(36000,4)],'VariableNames',...
                      {'dnum','ustr','vstr','wstr','Tson'});
            accplat = NaN(36000,3);
            uvwplat = NaN(36000,3);
            xyzplat = NaN(36000,3);
            euler = NaN(36000,3);
        end

        %% compute 10-min statistics

        % 10-min psd gps statistics
        % 1:dnum, 2:cog 3:sog 4:Ngps 5: Egps 6:lat 7:lon
        gps_10min_avg = interval_avg(gps.dnum, gps{:,2:end},dnum_10min_bin);
        da.lat_psl(iii) = gps_10min_avg(:,6);
        da.lon_psl(iii) = gps_10min_avg(:,7);
        Ngps = gps_10min_avg(:,4); Egps = gps_10min_avg(:,5);
        % cog from speed components
        gps_10min_avg(:,2) = atan2(gps_10min_avg(:,5),gps_10min_avg(:,4))*r2d;
        gps_10min_avg(:,2) = mod(gps_10min_avg(:,2)+360, 360);
        da.cog_psl(iii) = gps_10min_avg(:,2);
        % sog from speed components
        gps_10min_avg(:,3) = sqrt(gps_10min_avg(:,4).^2 + gps_10min_avg(:,5).^2);
        da.sog_psl(iii) = gps_10min_avg(:,3);
        clear gps_10min_avg

        gps_10min_std = interval_std(gps.dnum, gps{:,2:end},dnum_10min_bin);
        % cog std from unwrapped raw cog
        temp = interval_std(gps.dnum,unwrap(gps.cog*d2r),dnum_10min_bin)*r2d;
        da.cog_psl_std(iii) = temp(:,2);
        da.sog_psl_std(iii) = gps_10min_std(:,3);
        clear gps_10min_std temp

        gps_10min_min = interval_min(gps.dnum, gps{:,2:end},dnum_10min_bin);
        da.sog_psl_min(iii) = gps_10min_min(:,3);
        clear gps_10min_min

        gps_10min_max = interval_max(gps.dnum, gps{:,2:end},dnum_10min_bin);
        da.sog_psl_max(iii) = gps_10min_max(:,3);
        clear gps_10min_max

        % 10-min psd heading
        chead_10 = interval_avg(dnum_10Hz, chead, dnum_10min_bin);
        shead_10 = interval_avg(dnum_10Hz, shead, dnum_10min_bin);
        [hed_10,~] = cart2pol(chead_10(:,2), shead_10(:,2));
        hed_10 = hed_10*r2d;
        da.hed_psl(iii) = mod(hed_10+360, 360);
        temp = interval_std(dnum_10Hz,unwrap(heading),dnum_10min_bin);
        da.hed_psl_std(iii) = temp(:,2)*r2d;
        clear chead chead_10 shead shead_10 hed_10 temp

        % 10-min psd met statistics is 6x19 array:
        %{
            1   dnum            matlab datenum timestamp
            2   Tsea            sea snake surface SST, °C
            3   psp1            psp1 radiation, W/m2
            4   pir1            pir1 radiation, W/m2
            5   psp2            psp2 radiation, W/m2
            6   pir2            pir2 radiation, W/m2
            7   Tc1             pir1 case temperature, °C
            8   Td1             pir1 dome temperature, °C
            9   Tc2             pir2 case temperature, °C
            10  Td2             pir2 dome temperature, °C
            11  Pmb             sea level air pressure, mb(hPa), Vaisala PTB330
            12  rwdir_wxt       rel wind direction ,0-360 deg from bow
            13  rwspd_wxt       wind speed, m/s
            14  Tair_wxt        air temperature, °C
            15  RH_wxt          RH, %
            16  Pmb_wxt         sea level air pressure, mb(hPa)
            17  rain_wxt        rain rate, mm/hr
            18  U_wxt           rel U wind component, WXT, m/s
            19  V_wxt           rel V wind component, WXT, m/s
            20  sst_rosr        radiometric sst, °C
        %}
        met_10min_avg = interval_avg(met.dnum, met{:,2:end},dnum_10min_bin);
        da.sst_psl(iii) = met_10min_avg(:,2);
        da.rs_psl(iii) = (met_10min_avg(:,3) + met_10min_avg(:,3))/2;
        da.rl_psl(iii) = (met_10min_avg(:,4) + met_10min_avg(:,6))/2;
        da.Pmb_psl(iii) = met_10min_avg(:,11);
        da.ta_wxt(iii) = met_10min_avg(:,14);
        da.rh_wxt(iii) = met_10min_avg(:,15);
        da.Pmb_wxt(iii) = met_10min_avg(:,16);
        da.rain_wxt(iii) = met_10min_avg(:,17);
        da.Urel_wxt(iii) = met_10min_avg(:,18);
        da.Vrel_wxt(iii) = met_10min_avg(:,19);
        da.sst_rosr(iii) = met_10min_avg(:,20);
        [rwspd_wxt,rwdir_wxt] = uv_to_sd(da.Urel_wxt(iii),da.Vrel_wxt(iii));
        da.rwspd_wxt(iii) = rwspd_wxt;
        da.rwdir_wxt(iii) = rwdir_wxt;
        [spd_wxt,dir_wxt,~,~] = uv_rel_to_sd_true(da.Urel_wxt(iii),da.Vrel_wxt(iii),...
                                da.hed_psl(iii),da.cog_psl(iii),da.sog_psl(iii));
        da.wspd_wxt(iii) = spd_wxt;
        da.wdir_wxt(iii) = dir_wxt;
        clear met_10min_avg rwspd_wxt rwdir_wxt spd_wxt dir_wxt

        % 10-min motion-corrected wind statistics: 6x5 array as follows:
        % 1:jd, 2:ustr, 3:vstr, 4:wstr, 5:Tsonic
        uvwt_10min_avg = interval_avg(dnum_10Hz, uvw_str{:,2:end},dnum_10min_bin);
        uvwt_10min_avg(xxx,2:5) = NaN; % remove bad 10-min periods of son data
        da.Ubar(iii) = uvwt_10min_avg(:,2); % mean streamline U, no flow dist. correction

        uvwt_10min_std = interval_std(dnum_10Hz, uvw_str{:,2:end},dnum_10min_bin);
        uvwt_10min_std(xxx,2:5) = NaN;
        da.uvar(iii) = uvwt_10min_std(:,2).^2;
        da.vvar(iii) = uvwt_10min_std(:,3).^2;
        da.wvar(iii) = uvwt_10min_std(:,4).^2;
        da.tvar(iii) = uvwt_10min_std(:,5).^2;
        clear uvwt_10min_avg uvwt_10min_std

        % 10-min platform velocity std deviations: 6x4 array as follows:
        % 1:jd, 2:uplat, 3:vplat, 4:wplat
        plat_10min_std = interval_std(dnum_10Hz, uvwplat(:,:),dnum_10min_bin);
        da.uplat_std(iii) = plat_10min_std(:,2);
        da.vplat_std(iii) = plat_10min_std(:,3);
        da.wplat_std(iii) = plat_10min_std(:,4);
        clear plat_10min_std

        % compute psd true wind spd & dir from N/E sonic components
        % first, apply flow distortion correction to the raw wind speeds
        U = son.u/0.95;  % +5% for head-on wind
        V = son.v/1.15;  % -15% for side-on wind
        Un_10min = zeros(6,1)*NaN; % empty arrays
        Ue_10min = zeros(6,1)*NaN;
        ii = 1;
        for jj=1:6
            Un_10min(jj) = (nanmean1(U(ii:ii+Npts-1).*cos(heading(ii:ii+Npts-1)) + ...
                    V(ii:ii+Npts-1).*sin(heading(ii:ii+Npts-1)) + Ngps(jj)));
            Ue_10min(jj) = (nanmean1(U(ii:ii+Npts-1).*sin(heading(ii:ii+Npts-1)) - ...
                    V(ii:ii+Npts-1).*cos(heading(ii:ii+Npts-1)) + Egps(jj)));
            ii = ii + Npts;
        end
        Un_10min(xxx) = NaN; Ue_10min(xxx) = NaN;
        da.Un(iii) = Un_10min; da.Ue(iii) = Ue_10min;
        da.wdir_psl(iii) = atan2(-da.Ue(iii), -da.Un(iii))*r2d;
        ww = find(da.wdir_psl(iii)<0);
        da.wdir_psl(iii(ww)) = da.wdir_psl(iii(ww)) + 360;
        da.wspd_psl(iii) = sqrt(da.Un(iii).^2 + da.Ue(iii).^2);
        clear U V Un_10min Ue_10min ww

        % 10-min relative wind statistics
        rwdir = unwrap(atan2(son.v,-son.u)); % unwrapped rwdir in radians
        % average and convert to +/- 180 deg format
        rwdir_psd = interval_avg(dnum_10Hz,rwdir,dnum_10min_bin)*r2d;
        rwdir_psd = rwdir_psd(:,2); rwdir_psd(xxx) = NaN;
        rwdir_psd = mod(rwdir_psd + 360, 360);
        rwdir_psd(rwdir_psd>180) = rwdir_psd(rwdir_psd>180) - 360;
        da.rwdir_psl(iii) = rwdir_psd;
        rwdir_psd_std = interval_std(dnum_10Hz,rwdir,dnum_10min_bin)*r2d;
        rwdir_psd_std = rwdir_psd_std(:,2); rwdir_psd_std(xxx) = NaN;
        da.rwdir_psl_std(iii) = rwdir_psd_std;
        clear rwdir_psd rwdir_psd_std

        % 10-min raw sonic wind averages: 6x5 array as follows:
        % 1:jd, 2:u, 3:v, 4:w, 5:Tsonic
        son_10min_avg = interval_avg(dnum_10Hz,son{:,2:end},dnum_10min_bin);
        son_10min_avg(xxx,2:5) = NaN;
        da.rwspd_psl(iii) = sqrt(son_10min_avg(:,2).^2 + son_10min_avg(:,3).^2 + ...
                            son_10min_avg(:,4).^2);
        rwspd = sqrt(son.u.^2 + son.v.^2 + son.w.^2);
        son_10min_std = interval_std(dnum_10Hz,rwspd,dnum_10min_bin);
        son_10min_std(xxx,2) = NaN;
        da.rwspd_psl_std(iii) = son_10min_std(:,2);
        clear son_10min_avg rwspd son_10min_std

        % 10-min licor H2O and CO2 statistics: 6x5 array output
        % stats are for 'all' data, not screened for bad AGC
        % 1 dnum, 2 Licor_H2O(mmol/m3), 3 Licor_CO2(ppm), 4 Licor_T(°C), 5 Licor_P(mb), 6 AGC
        lic_10min_avg_all = interval_avg(dnum_10Hz,lic{:,2:end},dnum_10min_bin);
        da.AGC_lic(iii) = lic_10min_avg_all(:,6);
        ww = find(da.AGC_lic(iii)>65);
        da.h2o_lic_all(iii) = lic_10min_avg_all(:,2);
        da.h2o_lic_good(iii) = lic_10min_avg_all(:,2); da.h2o_lic_good(iii(ww)) = NaN;
        da.co2_lic_all(iii) = lic_10min_avg_all(:,3);
        da.co2_lic_good(iii) = lic_10min_avg_all(:,3); da.co2_lic_good(iii(ww)) = NaN;
        da.T_lic(iii) = lic_10min_avg_all(:,4);
        da.Pmb_lic(iii) = lic_10min_avg_all(:,5);
        
        
        lic_10min_std_all = interval_std(dnum_10Hz,lic{:,2:end},dnum_10min_bin);
        da.h2o_std_all(iii) = lic_10min_std_all(:,2);
        da.h2o_std_good(iii) = lic_10min_std_all(:,2); da.h2o_std_good(iii(ww)) = NaN;
        da.co2_std_all(iii) = lic_10min_std_all(:,3);
        da.co2_std_good(iii) = lic_10min_std_all(:,3); da.co2_std_good(iii(ww)) = NaN;
        clear lic_10min_avg_all lic_10min_std_all ww

        % 10-min specific humidities (q = g/kg)
        da.qa_psl(iii) = qair_p([da.ta_psl(iii),da.rh_psl(iii),da.Pmb_psl(iii)]);
        da.qa_ship(iii) = qair_p([da.ta_ship(iii),da.rh_ship(iii),da.Pmb_ship(iii)]);
        da.qa_wxt(iii) = qair_p([da.ta_wxt(iii),da.rh_wxt(iii),da.Pmb_wxt(iii)]);
        da.qs_psl(iii) = qsea_p([da.sst_psl(iii),da.Pmb_psl(iii)]);
        da.qs_ship(iii) = qsea_p([da.sst_ship(iii),da.Pmb_ship(iii)]);
        da.qs_rosr(iii) = qsea_p([da.sst_rosr(iii),da.Pmb_psl(iii)]);

        % 10-min mean air density (kg/m3) and H2O number density (mmol/m3) from bulk met
        da.rhoa_psl(iii) = rhoa(da.ta_wxt(iii),da.Pmb_wxt(iii),da.rh_wxt(iii));

        % compute 10-min abnd fast q from licor mmol/m3, using 10-min moist air density
        da.q_lic(iii) = (da.h2o_lic_all(iii)*1e-3).*Mw./da.rhoa_psl(iii); 
        q_10Hz_raw = NaN(Npts*6,1);
        lic_h2o = fillgaps(lic.h2o,1800);  
        ii=1;
        for jj=1:6
            q_10Hz_raw(ii:ii+Npts-1) = lic_h2o(ii:ii+Npts-1)*18/1000/da.rhoa_psl(iii(jj));
            ii = ii + Npts;
        end

        % decorrelate licor water vapor w/respect to 3-axis motion
        if decorr
            X = [ones(length(q_10Hz_raw),1),accplat,uvwplat];
            [q,mu_q] = interval_decorr(dnum_10Hz,q_10Hz_raw,dnum_10min_bin,X);
            mu_q_decorr(iii,2:13) = array2table(mu_q,'VariableNames',muq_vars(2:13));
        else
            q = q_10Hz_raw;
        end
        % compute q variance
        q_10min_std = interval_std(dnum_10Hz, q, dnum_10min_bin);
        da.qvar(iii) = q_10min_std(:,2).^2;
        % add q to uvw_str table
        uvwtq_str = addvars(uvw_str,q);
        clear uvw_str q_10min_std mu_q q X lic_h2o

        % finally, compute 10-min wu, wv, wt and wq covariances
        wu_cov_10min = interval_cov(dnum_10Hz,uvwtq_str.wstr,uvwtq_str.ustr,dnum_10min_bin);
        wv_cov_10min = interval_cov(dnum_10Hz,uvwtq_str.wstr,uvwtq_str.vstr,dnum_10min_bin);
        wTs_cov_10min = interval_cov(dnum_10Hz,uvwtq_str.wstr,uvwtq_str.Tson,dnum_10min_bin);
        wq_cov_10min = interval_cov(dnum_10Hz,uvwtq_str.wstr,uvwtq_str.q,dnum_10min_bin);
        wu_cov_10min(xxx,2) = NaN; da.wu_cov(iii) = wu_cov_10min(:,2);
        wv_cov_10min(xxx,2) = NaN; da.wv_cov(iii) = wv_cov_10min(:,2);
        wTs_cov_10min(xxx,2) = NaN; da.wTs_cov(iii) = wTs_cov_10min(:,2);
        wq_cov_10min(xxx,2) = NaN; da.wq_cov(iii) = wq_cov_10min(:,2);
        clear wu_cov_10min wv_cov_10min wTs_cov_10min wq_cov_10min

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

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

        output: 6 x 25 array:
         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      zo      roughness length (m)
         9      zot     thermal roughness length (m)
         10     zoq     moisture roughness length (m)
         11     Cd      wind stress transfer (drag) coefficient at height zu
         12     Ch      sensible heat transfer coefficient (Stanton number) at ht zu
         13     Ce      latent heat transfer coefficient (Dalton number) at ht zq
         14     L       Obukhov length scale (m)
         15     zet     Monin-Obukhov stability parameter zu/L
         16     dter    cool-skin temperature depression (degC)
         17     dqer    cool-skin humidity depression (degC)
         18     tkt     cool-skin thickness (m)
         19     RF      rain heat flux (W/m2)
         20     Cdn_10  neutral value of drag coefficient at 10m
         21     Chn_10  neutral value of Stanton number at 10m
         22     Cen_10  neutral value of Dalton number at 10m
         23     U10N    10m neutral wind speed, m/s
         24     T10N    10m neutral air temperature, degC
         25     Q10N    10m neutral air humidity, g/kg
        %}
        
        % ADJUST COARE INPUT VARIABLES AS REQUIRED
        % NO WAVE INPUTS THIS RUN
        coare_psd = coare35vn(da.wspd_psl(iii), da.zu_psl(iii),...
                    da.ta_wxt(iii), da.zt_psl(iii), da.rh_wxt(iii),...
                    da.zq_psl(iii), da.Pmb_wxt(iii), da.sst_rosr(iii),...
                    da.rs_psl(iii), da.rl_psl(iii), da.lat_psl(iii),...
                    600, da.rain_wxt(iii), NaN, NaN, 10);
        bb = find(strcmp(da_vars,'usb_psl'));
        da(iii,bb:bb+21) = coare_psd(:,1:22); % just using vars 1:22
        clear coare_psd

        %% compute covariance heat fluxes
        % correct sonic wT for humidity flux -  use C35 for <wq>=-usb*qsb
        % <wT> = <wTs> - 0.51*T*<wq>, <wq> = -u*q*
        da.wT_cov(iii) = da.wTs_cov(iii) - 0.51*(da.ta_wxt(iii)+C2K).* ...
                         -da.usb_psl(iii).*da.qsb_psl(iii)/1000;
        % 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
        da.hs_cov(iii) = da.wT_cov(iii).* da.rhoa_psl(iii) * 1004.67;

        % latent heat vap for water and seawater from SST
        % using default salinity...
        [Le_w_10min, Le_sw_10min] = Le_water(da.sst_rosr(iii), ones(6,1)*33);
        % hl = wq * rhoair * Le_w = g/kg m/s * kg/m3 * J/g = W/m2
        da.hl_cov(iii) = da.wq_cov(iii) .* da.rhoa_psl(iii) .* Le_sw_10min;
        % Webb corrected latent heat flux, using C35 webb correction
        da.hlc_cov(iii) = da.hl_cov(iii) + da.hlwebb_psl(iii);

        %% --------------------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;
        Ts_noise = NaN(6,1);

        df = fsam/(Npts-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 process 'good' 10-min son segments
            ii = 1 + Npts*(jj-1);
            [Puu(jj,:),Fxx(1,:)] = psd2(detrend(uvwtq_str.ustr(ii:ii+Npts-1))',Npts,fsam,hamming(Npts));
            Puu = real(Puu);
            [Pu(jj,:),Fx(1,:)] = specsmoo(Puu(jj,:),fsam);
            [Pvv(jj,:),~] = psd2(detrend(uvwtq_str.vstr(ii:ii+Npts-1))',Npts,fsam,hamming(Npts));
            Pvv = real(Pvv);
            [Pv(jj,:),~] = specsmoo(Pvv(jj,:),fsam);
            [Pww(jj,:),~] = psd2(detrend(uvwtq_str.wstr(ii:ii+Npts-1))',Npts,fsam,hamming(Npts));
            Pww = real(Pww);
            [Pw(jj,:),~] = specsmoo(Pww(jj,:),fsam);
            [Ptt(jj,:),~] = psd2(detrend(uvwtq_str.Tson(ii:ii+Npts-1))',Npts,fsam,hamming(Npts));
            Ptt = real(Ptt);
            [Pt(jj,:),~] = specsmoo(Ptt(jj,:),fsam); % This is Tsonic...
            [Pqq(jj,:),~] = psd2(detrend(uvwtq_str.q(ii:ii+Npts-1))',Npts,fsam,hamming(Npts));
            Pqq = real(Pqq);
            [Pq(jj,:),~] = specsmoo(Pqq(jj,:),fsam);

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

            % cospectra
            % wu
            [temp,~] = csd2(detrend(uvwtq_str.wstr(ii:ii+Npts-1)),detrend(uvwtq_str.ustr(ii:ii+Npts-1)),...
                                    Npts,fsam,hamming(Npts)); temp = real(temp);
            da.wu_csp(iii(jj)) = sum(temp*df); % integrate for wu_csp
            [Cuw(jj,:),~] = specsmoo(temp,fsam);
            % wv
            [temp,~] = csd2(detrend(uvwtq_str.wstr(ii:ii+Npts-1)),detrend(uvwtq_str.vstr(ii:ii+Npts-1)),...
                                    Npts,fsam,hamming(Npts)); temp = real(temp);
            da.wv_csp(iii(jj)) = sum(temp*df); % integrate for wv_csp
            [Cvw(jj,:),~] = specsmoo(temp,fsam);
            % wTs
            [temp,~] = csd2(detrend(uvwtq_str.wstr(ii:ii+Npts-1)),detrend(uvwtq_str.Tson(ii:ii+Npts-1)),...
                                    Npts,fsam,hamming(Npts)); temp = real(temp);
            da.wTs_csp(iii(jj)) = sum(temp*df); % integrate for wTs_csp
            [Ctw(jj,:),~] = specsmoo(temp,fsam);
            % wTs_raw
            [temp,~] = csd2(detrend(uvwtq_str.wstr(ii:ii+Npts-1)),detrend(son.Tsonic(ii:ii+Npts-1)),...
                                    Npts,fsam,hamming(Npts)); temp = real(temp);
            [Ctw_raw(jj,:),~] = specsmoo(temp,fsam);
            % wq
            [temp,~] = csd2(detrend(uvwtq_str.wstr(ii:ii+Npts-1)),detrend(uvwtq_str.q(ii:ii+Npts-1)),...
                                    Npts,fsam,hamming(Npts)); temp = real(temp);
            da.wq_csp(iii(jj)) = sum(temp*df); % integrate for wq_csp
            [Cqw(jj,:),~] = specsmoo(temp,fsam);
            % wq_raw
            [temp,~] = csd2(detrend(uvwtq_str.wstr(ii:ii+Npts-1)),detrend(q_10Hz_raw(ii:ii+Npts-1)),...
                                    Npts,fsam,hamming(Npts)); temp = real(temp);
            [Cqw_raw(jj,:),~] = specsmoo(temp,fsam);
            clear temp
            ii = ii + Npts;
        end

        da.Ts_noise(iii) = Ts_noise;

        % Decide if noise subtraction from Ts spectrum is justified.
        % Compare Pt*Fx magnitudes at 0.8-1.2 Hz and 3.2-5 Hz for the hourly mean
        % smoothed spectrum. If magnitude at the higher frequencies is greater,
        % then keep the adjusted spectrum.  Otherwise, replace Pt2 with the original
        % Pt spectrum.
        hi = Fx>3.2; lo = Fx>0.8 & Fx<1.2; testPt = nanmean1(Pt,1).*Fx;
        hi_mean = nanmean1(testPt(hi)); HIM = sprintf('%6.4f',hi_mean);
        lo_mean = nanmean1(testPt(lo)); LOM = sprintf('%6.4f',lo_mean);
        disp(['Pt-hi-f = ',HIM,'  Pt-lo-f = ',LOM]);
        if hi_mean < lo_mean % no apparent noise tail this hour...
            Pt2 = Pt;        % so don't keep the adjusted spectra
        end

        % compute wT, hs, hl, hlc for cospectral integral fluxes, as above
        da.wT_csp(iii) = da.wTs_csp(iii) + 0.51*(da.ta_wxt(iii)+C2K).* ...
                         da.usb_psl(iii).*da.qsb_psl(iii)*1e-3;
        da.hs_csp(iii) = da.wT_csp(iii).* da.rhoa_psl(iii) * 1004.67;
        da.hl_csp(iii) = da.wq_csp(iii).* da.rhoa_psl(iii).* Le_sw_10min;
        da.hlc_csp(iii) = da.hl_csp(iii) + da.hlwebb_psl(iii);

        %% save spectra to sp array
        %{
            10-min spec output arrays for spectra and cospectra: 144x481
    
            column designations as follows, all spectra are 48 bin
    
            COLUMNS     VARIABLE        DESCRIPTION
            1           dnum            10-min timestamp
            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(iii,2:end) = [Fx2, Pu, Pv, Pw, Pt, Pq, Cuw, Cvw, Ctw, Cqw];


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

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

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

%         fminT = find(Fxx>0.28,1,'first');
%         fmaxT = find(Fxx<0.4,1,'last');
%         da.Ct2a(iii) = idiss_struc_func(Ptt2, Fxx, fminT, fmaxT, da.rwspd_psl(iii), 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
        da.Cq2a(iii) = idiss_struc_func(Pq, Fx, fminq, fmaxq, da.rwspd_psl(iii), 0);

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

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

        da.Cu2a(iii) = da.Cu2a(iii)./f1;
        da.Cw2a(iii) = da.Cw2a(iii)./f2;
        da.Ct2a(iii) = da.Ct2a(iii)./f3;
        da.Cq2a(iii) = da.Cq2a(iii)./f3;

        % compute usid, tsid, qsid, Lid
        % (Cu2, Ct2, Cq2, qbar, Tbar, SST, zu, zq, L_bulk, usb, tsb, qsb, dter)
        [da.usid_a(iii),da.tsid_a(iii),da.qsid_a(iii),da.Lid_a(iii)] = ...
            idiss(da.Cu2a(iii),da.Ct2a(iii),da.Cq2a(iii),da.qa_wxt(iii),...
            da.ta_wxt(iii),da.sst_rosr(iii),zu,z_wxt,da.L_psl(iii),...
            da.usb_psl(iii),da.tsb_psl(iii),da.qsb_psl(iii),da.dter_psl(iii));
        % qsid converted from kg/kg to g/kg
        da.qsid_a(iii) = da.qsid_a(iii)*1000;

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

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

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

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

        % compute usid, tsid, qsid, Lid
        % (Cu2, Ct2, Cq2, qbar, Tbar, SST, zu, zq, L_bulk, usb, tsb, qsb, dter)
        [da.usid_b(iii),da.tsid_b(iii),da.qsid_b(iii),da.Lid_b(iii)] = ...
            idiss(da.Cu2b(iii),da.Ct2b(iii),da.Cq2b(iii),da.qa_wxt(iii),...
            da.ta_wxt(iii),da.sst_rosr(iii),zu,z_wxt,da.L_psl(iii),...
            da.usb_psl(iii),da.tsb_psl(iii),da.qsb_psl(iii),da.dter_psl(iii));
        % qsid converted from kg/kg to g/kg
        da.qsid_b(iii) = da.qsid_b(iii)*1000;

        clear fminU fmaxU fminq fmaxq fminT fmaxT f1 f2 f3

        %% save 10_Hz wind and motion data as tables

        rwdir = mod(rwdir*r2d + 360, 360);
        rwdir(rwdir>180) = rwdir(rwdir>180) - 360;
        flx = addvars(uvwtq_str,rwdir);
        fpath = fullfile(path_uvwStreamData,['uvwtq_',DDD,'_',HH,'_',dcname,'.mat']);
        save (fpath,'flx');

        mot=[dnum_10Hz,uvwplat,xyzplat,euler];
        mot_vars = {'dnum','uplat','vplat','wplat','xplat','yplat','zplat','phi','theta','psi'};
        mot = array2table(mot,'VariableNames',mot_vars);
        fpath = fullfile(path_motionData,['motion_',DDD,'_',HH,'_',dcname,'.mat']);
        save (fpath,'mot');

        %% selected wind correction and spectra plots every 6 hrs
        if plotit
            if hh==0 || hh==6 || hh==12 || hh==18
                % cospectra plots
                if any(isfinite(Cuw))
                    wu = nanmean1(Cuw,1);
                    figure(4);
                    semilogx(Fx,Fx.*wu,'b-',[1e-3 10 ],[0 0],'LineWidth',1);
                    title([DDD,' ',HH,' wu Cospectrum']); grid;
                    ylabel('C_{uw}(f) x f'); xlabel('f (Hz)');
                    ppath = fullfile(path_spectraPlots,'wu_cospectrum',['wuCospec_',DDD,...
                        '_',HH,'_',dcname,graphformat]);
                    print(graphdevice,ppath);
                    close all;
                end

                if any(isfinite(Cqw)) | any(isfinite(Cqw_raw))
                    wq = nanmean1(Cqw,1); wq_raw = nanmean1(Cqw_raw,1);
                    figure(5);
                    semilogx(Fx,Fx.*wq,'b-',Fx,Fx.*wq_raw,'b--',[1e-3 10 ],[0 0],'LineWidth',1);
                    title([DDD,' ',HH,' wq Cospectrum: raw(dashed), decorr(solid)']); grid;
                    ylabel('C_{wq}(f) x f'); xlabel('f (Hz)');
                    ppath = fullfile(path_spectraPlots,'wq_cospectrum',['wqCospec_',DDD,'_',...
                        HH,'_',dcname,graphformat]);
                    print(graphdevice,ppath);
                    close all;
                end

                if any(isfinite(Ctw)) | any(isfinite(Ctw_raw))
                    wt = nanmean1(Ctw,1); wt_raw = nanmean1(Ctw_raw,1);
                    figure(6);
                    semilogx(Fx,Fx.*wt,'b-',Fx,Fx.*wt_raw,'r--',[1e-3 10 ],[0 0],'LineWidth',1);
                    title([DDD,' ',HH,' wT Cospectrum: raw(dashed), decorr(solid)']); grid;
                    ylabel('C_{wt}(f) x f'); xlabel('f (Hz)');
                    ppath = fullfile(path_spectraPlots,'wt_cospectrum',['wtCospec_',DDD,...
                        '_',HH,'_',dcname,graphformat]);
                    print(graphdevice,ppath);
                    close all;
                end

                % spectra plots
                if any(isfinite(Pw))
                    Sw = nanmean1(Pw,1);
                    [Sw_raw,~] = psd2(detrend(sonic(:,4)),Npts*6,10,hamming(Npts*6));
                    [Sw_raw,Fraw] = specsmoo(Sw_raw,10);
                    figure(7);
                    loglog(Fx,Fx.*Sw,'b-',Fraw,Fraw.*Sw_raw,'b--','LineWidth',1); hold on;
                    loglog(Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*Sw(35).*Fx(35),'r',...
                        'LineWidth',1);hold off;
                    ylabel('S_{ww} x f'); xlabel('f (Hz)'); legend('corrected','raw','location','northeast');
                    title([DDD,' ',HH,' Spectrum of vertical wind velocity, W']); grid;
                    ppath = fullfile(path_spectraPlots,'w_spectrum',['wSpectrum_',DDD,...
                        '_',HH,'_',dcname,graphformat]);
                    print(graphdevice,ppath);
                    close all;
                end

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

                if any(isfinite(Pq))
                    Sq = nanmean1(Pq,1);
                    figure(9);
                    loglog(Fx,Fx.*Sq,'b-',Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*Sq(35).*Fx(35),...
                        'r','LineWidth',1);
                    ylabel('S_{qq} x f'); xlabel('f (Hz)');
                    title([DDD,' ',HH,' Spectrum of specific humidity, q']); grid;
                    ppath = fullfile(path_spectraPlots,'q_spectrum',['qSpectrum_',DDD,...
                        '_',HH,'_',dcname,graphformat]);
                    print(graphdevice,ppath);
                    close all;
                end

                if any(isfinite(Pu)) | any(isfinite(Pv))
                    Su = nanmean1(Pu,1);
                    Sv = nanmean1(Pv,1);
                    [Su_raw,~] = psd2(detrend(sonic(:,2)),Npts*6,10,hamming(Npts*6));
                    [Su_raw,~] = specsmoo(real(Su_raw),10);
                    [Sv_raw,~] = psd2(detrend(sonic(:,3)),Npts*6,10,hamming(Npts*6));
                    [Sv_raw,~] = specsmoo(real(Sv_raw),10);
                    figure(10);
                    loglog(Fx,Fx.*Su,'b-',Fraw,Fraw.*Su_raw,'b--',Fx,Fx.*Sv,'g-',...
                        Fraw,Fraw.*Sv_raw,'g--','LineWidth',1); hold on;
                    loglog(Fx(17:48),((Fx(17:48)).^(-2/3))./(Fx(35).^(-2/3)).*Su(35).*Fx(35),...
                        'r','LineWidth',1); hold off
                    ylabel('S_{uu} x f'); xlabel('f (Hz)');
                    legend('U','U raw','V','V raw','-2/3','location','northeast');
                    title([DDD,' ',HH,' U V Wind Spectra']); grid;
                    ppath = fullfile(path_spectraPlots,'uv_spectrum',['uvSpectrum_',DDD,...
                        '_',HH,'_',dcname,graphformat]);
                    print(graphdevice,ppath);
                    close all;
                end
            end
        end
    end % end of 24 hour loop

    %% save da, sp, lags and mu files
    da_fname = ['da_',DDD,'_',dcname,'.mat'];
    dfl = fullfile(path_da,da_fname);
    save (dfl,'da');

    sp = real(sp);
    sp_fname = ['sp_',DDD,'_',dcname,'.mat'];
    dfl = fullfile(path_sp,sp_fname);
    save (dfl,'sp');

    lag_fname = ['lag_',DDD,'_',dcname,'.mat'];
    dfl = fullfile(path_lag,lag_fname);
    save (dfl,'lags');

    mu_fname = ['mu_',DDD,'_',dcname,'.mat'];
    dfl = fullfile(path_mu,'w_Ts',mu_fname);
    save (dfl,'mu_decorr');

    mu_q_fname = ['mu_q_',DDD,'_',dcname,'.mat'];
    dfl = fullfile(path_mu,'q',mu_q_fname);
    save (dfl,'mu_q_decorr');

    %% daily plots
    
    % a few parameters...
    wtv = -da.usb_psl.*(da.tsb_psl+0.61.*(da.ta_wxt+C2K).*da.qsb_psl);
    ws = (9.83./(da.ta_wxt+C2K).*abs(wtv)*600).^0.333;  % convective velocity
    ug = sqrt(da.usb_psl.^2 + ws.^2);
    ugw = sqrt(1.5*da.usb_psl.^2 + ws.^2);
    ugu = sqrt(8*da.usb_psl.^2 + 1.*ws.^2);
    gg = sqrt(1 + 1.6*ws.^2./da.wspd_psl.^2);

    % filters
    % kk is for rwdir, 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.
    % hh is for kk & licor AGC < 65
    % jj is for kk & isfinite(usid) & (abs(-wu-usid)^2<.02+.006*U)
    % mm; both hh and jj for hlid

    kk = find(abs(da.rwdir_psl)<60 & da.hed_psl_std<5 &...
               da.sog_psl_std<0.6 & da.vplat_std<0.8 &...
               (sqrt(da.wvar)./ugw < .7+.0015*da.wspd_psl.^2) &...
               (sqrt(da.vvar)./ugu<1.4) & da.missing_son<100 & ...
               da.badPnts_son<100 & da.rain_wxt<5);
    hh = intersect(kk, find(da.AGC_lic < 65));
    jj = intersect(kk, find(isfinite(da.usid_a)));
    mm = intersect(hh,jj);

    dt_10min = datetime(da.dnum,'ConvertFrom','datenum'); % datetime for plots

    % Air/Sea Temperatures
    figure;
    subplot(2,1,1);
        plot(dt_10min,da.ta_ship,'r-',dt_10min,da.ta_psl,'b-',dt_10min,...
             da.ta_wxt,'g-','LineWidth',2);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Air/Sea Temperatures',cruise,Vdate(1),Vdate(2),...
            Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        ylabel('Air Temp (C)'); grid on; grid minor;
        legend('Ship','PSL','WXT','Location','Best');
    subplot(2,1,2); 
        plot(dt_10min,da.sst_ship,'r-',dt_10min,da.sst_psl,'b-',dt_10min,...
             da.sst_rosr,'g-','LineWidth',2);
        ylabel('SST (C)'); grid on; grid minor;
        legend('Ship','PSL','ROSR','Location','Best');
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',...
            {'NOAA/ESRL/PSL/BoundaryLayerObs'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,'Air_Sea_Temps',['Air_Sea_Temps_',cruise,'_',sprintf('%02i_%02i_%03i',...
            Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end
    close all;

    % RH/q
    Licor_H2O_a = (da.h2o_lic_all*1e-3).*Rgas_universal.*(da.ta_wxt+C2K)./(da.Pmb_wxt*1e2); % mmol/m3 to mol/mol
    Licor_H2O_m = (da.h2o_lic_all*1e-3).*Mw./da.rhoa_psl; % mmol/m3 to g/kg
    licor_rh = relhum5([da.ta_wxt,Licor_H2O_m,da.Pmb_wxt]);
    dd = find(da.AGC_lic < 65);
    figure;
    subplot(2,1,1); plot(dt_10min,da.rh_ship,'r-',dt_10min,da.rh_psl,'b-',dt_10min,da.rh_wxt,'g-',...
        dt_10min(dd),licor_rh(dd),'ok','LineWidth',2);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  RH & Specific Humidity',cruise,Vdate(1),...
            Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        ylabel('RH');  grid on; grid minor; ylim([40,100]);
        legend('Ship','PSL','WXT','Licor','Location','Best');
    subplot(2,1,2); plot(dt_10min,da.qa_ship,'r-',dt_10min,da.qa_wxt,'b-',dt_10min(dd),Licor_H2O_m(dd),...
        'ok','LineWidth',2);
        ylabel('q (g/kg)'); grid on; grid minor;
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',...
            {'NOAA/ESRL/PSL/BoundaryLayerObs'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,'RH',['RH_q_',cruise,'_',sprintf('%02i_%02i_%03i',...
            Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end
    close all;

    % True Wspd/Wdir
    figure;
    subplot(2,1,1); plot(dt_10min,da.wspd_ship,'r-',dt_10min,da.wspd_psl,'b-',dt_10min,da.wspd_wxt,'g-',...
        'LineWidth',2);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  True Wind',cruise,Vdate(1),Vdate(2),...
            Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        ylabel('True Wind Speed (m/s)'); grid on; grid minor;
        legend('Ship','PSL','WXT','Location','Best');
    subplot(2,1,2); plot(dt_10min,da.wdir_ship,'r.',dt_10min,da.wdir_psl,'b.',dt_10min,da.wdir_wxt,'g.',...
        'LineWidth',2);
        ylabel('True Wind Dir'); grid on; grid minor;
        ylim([0 360]); yticks([0,45,90,135,180,225,270,315,360]);
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',...
            {'NOAA/ESRL/PSL/BoundaryLayerObs'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,'True_Wind',['True_Wind_',cruise,'_',sprintf('%02i_%02i_%03i',...
            Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end
    close all;

    % COG / SOG
    figure;
    subplot(2,1,1); plot(dt_10min,da.sog_ship,'r.',dt_10min,da.sog_psl,'b.','MarkerSize',8);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  SOG COG',cruise,Vdate(1),Vdate(2),...
            Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        ylabel('SOG (m/s)'); grid on; grid minor;
        legend('Ship','PSL','Location','Best');
    subplot(2,1,2); plot(dt_10min,da.cog_ship,'r.',dt_10min,da.cog_psl,'b.','MarkerSize',8);
        xlabel('Hour (UTC)');ylabel('COG'); grid on; grid minor; ylim([0 360]);
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',...
            {'NOAA/ESRL/PSL/BoundaryLayerObs'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,'SOG_COG',['SOG_COG_',cruise,'_',sprintf('%02i_%02i_%03i',...
            Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end
    close all;

    % Relative Wspd/Wdir
    figure;
    Uplim = ones(144,1)*60;
    LoLim = ones(144,1)*-60;
    subplot(2,1,1); plot(dt_10min,da.rwspd_psl,'b-',dt_10min,da.rwspd_ship,'r-','LineWidth',2);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSL Relative Wind',cruise,Vdate(1),...
            Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        ylabel('Rel Wind Speed (m/s)'); grid on; grid minor;
        legend('PSL','Ship','Location','Best');
    subplot(2,1,2); plot(dt_10min,da.rwdir_ship,'r.',dt_10min,da.rwdir_psl,'b.','MarkerSize',8); hold on;
        plot(dt_10min,Uplim,'k--',dt_10min,LoLim,'k--');
        ylabel('Rel Wind Dir'); grid on; grid minor;
        ylim([-180 180]); yticks([-180,-90,-60,0,60,90,180]);
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',...
            {'NOAA/ESRL/PSL/BoundaryLayerObs'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,'Rel_Wind',['Rel_Wind_',cruise,'_',sprintf('%02i_%02i_%03i',...
            Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end
    close all;

    % Rain
    temp = da.rain_wxt;
    temp(isnan(temp)) = 0;
    accm_wxt = cumsum(temp)/6; % mm/hr to mm/10 min
    temp = da.rain_psl;
    temp(isnan(temp)) = 0;
    accm_psd = cumsum(temp)/6;
    figure;
    subplot(2,1,1); plot(dt_10min,da.rain_psl,'b',dt_10min,da.rain_wxt,'g','LineWidth',3);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Rain',cruise,Vdate(1),Vdate(2),...
            Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        ylabel('Rain Rate (mm/hr)'); grid on; grid minor;
        legend('ORG','WXT','Location','Best');
    subplot(2,1,2); plot(dt_10min,accm_psd,'b',dt_10min,accm_wxt,'g','LineWidth',3);
        ylabel('Accumulated Rain (mm)'); grid;
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',...
            {'NOAA/ESRL/PSL/BoundaryLayerObs'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,'Rain',['Rain_',cruise,'_',sprintf('%02i_%02i_%03i',...
            Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end
    close all;

    % Std Dev Platform Motion & Tilt Angle
    figure;
    subplot(2,1,1); plot(dt_10min,da.uplat_std,'b-',dt_10min,da.vplat_std,'r-',...
        dt_10min,da.wplat_std,'g-','LineWidth',2);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Motion Std Dev & Tilt',cruise,Vdate(1),...
            Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        ylabel('Motion Std Dev (m/s)'); grid on; grid minor;
        legend('U plat','V plat','W plat','Location','Best');
    subplot(2,1,2); plot(dt_10min,da.tilt,'b-','LineWidth',2);
        ylabel('Wind Tilt Angle (deg)'); grid on; grid minor;
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',...
            {'NOAA/ESRL/PSL/BoundaryLayerObs'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,'Motion_Tilt',['Motion_Tilt_',cruise,'_',sprintf('%02i_%02i_%03i',...
            Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end
    close all;

    % Heat Fluxes
    % compute ID hs and hl
    hsid = -da.rhoa_psl*1004.67.*da.tsid_a.*da.usid_a;               % idiss method a
    [Le_w, ~] = Le_water(da.sst_rosr, ones(height(da))*33);           % using default salinity
    hlid = -da.rhoa_psl.*Le_w.*da.qsid_a.*da.usid_a - da.hlwebb_psl; % Webb corrected, g/kg
    maxhl = max(max(da.hlb_psl),max(da.hlc_cov)); minhl = min(min(da.hlb_psl),min(da.hlc_cov));
    maxhs = max(max(da.hsb_psl),max(da.hs_cov)); minhs = min(min(da.hsb_psl),min(da.hs_cov));
    figure;
    subplot(2,1,1); plot(dt_10min,zeros(size(dt_10min)),'k-','LineWidth',1); hold on
        h1 = plot(dt_10min,da.hlb_ship,'r-',dt_10min,da.hlb_psl,'b-',dt_10min(hh),...
        da.hlc_cov(hh),'go',dt_10min(mm),hlid(mm),'mo','LineWidth',2); hold off;
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Heat Fluxes',cruise,Vdate(1),...
            Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        ylabel('Latent Heat (W/m^2)'); grid on; grid minor;
        ylim([minhl-0.1*abs(minhl), maxhl+0.1*abs(maxhl)]);
        legend(h1,'Bulk-ship','Bulk-psd','Cov-Webb','ID-a','Location','Best');
    subplot(2,1,2); plot(dt_10min,zeros(size(dt_10min)),'k-','LineWidth',1); hold on
        plot(dt_10min,da.hsb_ship,'r-',dt_10min,da.hsb_psl,'b-',dt_10min(kk),...
        da.hs_cov(kk),'go',dt_10min(kk),hsid(kk),'mo','LineWidth',2); hold off;
        ylabel('Sensible Heat (W/m^2)'); grid on; grid minor;
        ylim([minhs-0.1*abs(minhs), maxhs+0.1*abs(maxhs)]);
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',...
            {'NOAA/ESRL/PSL/BoundaryLayerObs'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,'Heat_Fluxes',['Heat_Fluxes_',cruise,'_',sprintf('%02i_%02i_%03i',...
            Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end
    close all;

    % Solar/IR radiation
    figure;
    subplot(2,1,1); plot(dt_10min,da.rs_ship,'r-',dt_10min,da.rs_psl,'b-','linewidth',2);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Solar/IR Radiation',cruise,Vdate(1),...
            Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
        ylabel('Solar Downwelling (W/m^2)'); grid on; grid minor;
        legend('Ship','PSL','Location','Best');
    subplot(2,1,2); plot(dt_10min,da.rl_ship,'r-',dt_10min,da.rl_psl,'b-','linewidth',2);
        ylabel('IR Downwelling (W/m^2)'); grid;
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',...
            {'NOAA/ESRL/PSL/BoundaryLayerObs'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
    orient tall
    if plotit
        ppath = fullfile(path_dailyPlots,'Radiation',['Radiation_',cruise,'_',sprintf('%02i_%02i_%03i',...
            Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end
    close all;

    % Stress
    figure; plot(dt_10min,da.taub_ship,'r-',dt_10min,da.taub_psl,'b-',dt_10min(kk),...
        -(da.rhoa_psl(kk).*da.wu_cov(kk)),'ko',dt_10min(jj),(da.rhoa_psl(jj).*da.usid_a(jj).^2),...
        'go',dt_10min(jj),(da.rhoa_psl(jj).*da.usid_b(jj).^2),'mo','LineWidth',2);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Streamwise Stress',cruise,Vdate(1),...
        Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    ylabel('Stress (N/m^2)'); grid on; grid minor; ylim([-0.02 max(0,max(da.taub_psl))*1.25]);
    legend('Bulk-ship','Bulk-psd','Cov','ID-A','ID-B','Location','Best');
    annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',...
        {'NOAA/ESRL/PSL/BoundaryLayerObs'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
    if plotit
        ppath = fullfile(path_dailyPlots,'Stress',['Stress_',cruise,'_',sprintf('%02i_%02i_%03i',...
            Vdate(2),Vdate(3),ddd),graphformat]);
        print(graphdevice,ppath);
    end
    close all;

    disp(['FINISHED motcorr ',cruise,' yearday ',DDD,'.']);

end % end daily loop

%% concatenate daily files into one whole-cruise file
cd(path_da);
flist =  dir('da*.mat');
for ii = 1:length(flist)
    load(flist(ii).name,'da');
    if ii==1
        temp = da;
    else
        temp = [temp;da];
    end
end
da = temp;
fname = ['FATIMA_2022_da_',dcname,'.mat'];
dfl = fullfile(path_da,fname);
save(dfl,'da');

cd(path_lag);
flist =  dir('lag*.mat');
for ii = 1:length(flist)
    load(flist(ii).name,'lags');
    if ii==1
        temp = lags;
    else
        temp = [temp;lags];
    end
end
lags = temp;
fname = ['FATIMA_2022_lags_',dcname,'.mat'];
dfl = fullfile(path_lag,fname);
save(dfl,'lags');

cd(fullfile(path_mu,'w_ts'));
flist =  dir('mu*.mat');
for ii = 1:length(flist)
    load(flist(ii).name,'mu_decorr');
    if ii==1
        temp = mu_decorr;
    else
        temp = [temp;mu_decorr];
    end
end
mu_decorr = temp;
fname = ['FATIMA_2022_mu_w_Ts_',dcname,'.mat'];
dfl = fullfile(path_mu,'w_ts',fname);
save(dfl,'mu_decorr');

cd(fullfile(path_mu,'q'));
flist =  dir('mu*.mat');
for ii = 1:length(flist)
    load(flist(ii).name,'mu_q_decorr');
    if ii==1
        temp = mu_q_decorr;
    else
        temp = [temp;mu_q_decorr];
    end
end
mu_q_decorr = temp; clear temp;
fname = ['FATIMA_2022_mu_q_',dcname,'.mat'];
dfl = fullfile(path_mu,'q',fname);
save(dfl,'mu_q_decorr');

cd(path_sp);
flist =  dir('sp*.mat');
for ii = 1:length(flist)
    load(flist(ii).name,'sp');
    if ii==1
        temp = sp;
    else
        temp = [temp;sp];
    end
end
sp = temp; clear temp;
fname = ['FATIMA_2022_sp_',dcname,'.mat'];
dfl = fullfile(path_sp,fname);
save(dfl,'sp');

cd(path_prog);

disp('END OF motcorr PROCESSING')
