disp('save_flux_hr_FATIMA_2022');

%{
   Saves filtered 10-min output variable table flx10 as follows:
              Variable    Units     Description
    1         dnum        none      matlab datenum timestamp (decimal day)
    2         wspd        m/s       true wind speed
    3         wdir        deg       true wind direction
    4         rwspd       m/s       relative wind speed
    5         rwdir       deg       relative wind direction
    6         ta          degC      air temperature
    7         ta_wxt      degC      air temperature, WXT sensor
    8         rh          percent   relative humidity
    9         rh_wxt      percent   relative humidity, WXT sensor
    10        sst         degC      sea surface temperature
    11        sst_rosr    degC      radiometruc sea surface temperature
    12        Pmb         mb        sea level air pressure
    13        qs          g/kg      surface saturation specific humidity
    14        qa          g/kg      specific humidity
    15        rs          W/m2      downwelling solar radiative flux
    16        rl          W/m2      downwelling IR flux
    17        rain        mm/hr     rain rate
    18        hsc         W/m2      turbulent sensible heat flux
    19        hsib        W/m2      inertial dissipation sensible heat flux
    20        hsb         W/m2      bulk sensible heat flux
    21        hlc         W/m2      turbulent latent heat flux
    22        hlib        W/m2      inertial dissipation latent heat flux
    23        hlb         W/m2      bulk latent heat flux
    24        tauc        N/m2      turbulent streamwise wind stress
    25        taucx       N/m2      turbulent cross-stream wind stress
    26        tauib       N/m2      inertial dissipation wind stress
    27        taub        N/m2      bulk wind stress
    28        ct          none      Ct structure func parameter
    29        cq          none      Cq structure func parameter
    30        cu          none      Cu structure func parameter
    31        cw          none      Cw structure func parameter
    32        jplume      none      exhaust plume filter flag
    33        jmanuv      none      ship maneuver filter flag
    34        tilt        deg       streamline tilt in wind at location of sonic anemometer
    35        RF          W/m2      rain heat flux
    36        hlwebb      W/m2      latent heat 'Webb' correction
    37        rlclr       W/m2      clear sky IR downwelling radiative flux model
    38        rsclr       W/m2      clear sky solar downwelling radiative flux model
    39        u10n        m/s       10m neutral wind speed
    40        t10n        degC      10m neutral air temperature
    41        q10n        g/kg      10m neutral specific humidity
    42        usid_flag   none      quality flag for ID ustar
    43        qsid_flag   none      quality flag for ID qstar
    44        tsid_flag   none      quality flag for ID tstar
    45        zu          m         wind measurement height
    46        zt          m         air temperature measurement height
    47        zq          m         humidity measurement height
    48        hed         deg       ship heading
    49        sog         m/s       ship speed over ground from gps
    50        cog         deg       ship course over ground from gps
    51        lat         degN      latitude
    52        lon         degE      longitude
    53        wspd_ship   m/s       true wind speed from ship sensor
    54        wdir_ship   deg       true wind direction from ship sensor
    55        ta_ship     degC      air temperature from ship sensor
    56        rh_ship     percent   relative humidity from ship sensor
    57        sst_ship    degC      sea surface temperature from ship sensor
    58        Pmb_ship    mb        sea level air pressure from ship sensor
    59        qs_ship     g/kg      surface saturation specific humidity from ship sensors
    60        qa_ship     g/kg      specific humidity from ship sensor
    61        rs_ship     W/m2      downwelling solar radiative flux from ship sensor
    62        rl_ship     W/m2      downwelling IR flux from ship sensor
    63        q_lic       g/kg      specific humidity from Licor 7500
    64        q_lic_std   g/kg      std dev specific humidity from Licor 7500
    65        co2_lic     ppm       co2 ppm from Licor 7500
    66        co2_lic_std ppm       std dev co2 ppm from Licor 7500
    67        rnl         W/m2      net longwave radiation (into ocean)
    68        hnetb       W/m2      net heat flux into ocean, bulk hs & hl (into ocean)
    69        hnetc       W/m2      net heat flux into ocean, turb hs & hl (into ocean)
    70        hnetib      W/m2      net heat flux into ocean, ID hs & hl (into ocean)
%}
fclose('all');
clear flux1hr;

%% initialize
path_f60_vars = fullfile(path_prog,'f60_vars.xlsx');
f60_params = readcell(path_f60_vars); % load variable names, units and descriptions
f60_vars = f60_params(:,2);
f60_units = f60_params(:,3);
f60_descript = f60_params(:,4);

% prepare empty table for output
% names, units and descriptions are read in from an excel spreadsheet
% to add variables or adjust metadata, edit the excel sheet
npz = height(da_hr);    % number of rows of data
nnn = length(f60_vars); % number of columns of output
flxhr = array2table(NaN(npz,nnn),'VariableNames',f60_vars);
flxhr.Properties.VariableUnits = f60_units;
flxhr.Properties.VariableDescriptions = f60_descript;

%% prepare a few of the output variables
reldir2 = rwdir_best;
reldir2(reldir2<0) = reldir2(reldir2<0)+360;
% define good ID data flags
yfilt_u = isfinite(yusib);                  % wind
yfilt_q = isfinite(yqsib_lic);              % Licor
yfilt_t = isfinite(ytsib);                  % sonic T
ytaucx = -(da_hr.rhoa_psl.*da_hr.wv_csp);  % x-stream stress
ytauid = (da_hr.rhoa_psl.*yusib.^2);       % ID stress

% copy variables into table
flxhr.dnum = da_hr.dnum;
flxhr.wspd = wspd_hr;
flxhr.wdir = wdir_hr;
flxhr.rwspd = rwspd_hr;
flxhr.rwdir = rwdir_hr;
flxhr.ta = best.ta;
flxhr.ta_wxt = da_hr.ta_wxt;
flxhr.rh = best.rh;
flxhr.rh_wxt = da_hr.rh_wxt;
flxhr.sst = best.ts;
flxhr.sst_rosr = da_hr.sst_rosr;
flxhr.Pmb = best.Pmb;
flxhr.qs = best.qs;
flxhr.qa = best.qa;
flxhr.rs = best.rs;
flxhr.rl = best.rl;
flxhr.rain = best.rain;
flxhr.hsc = yhsc;
flxhr.hsib = yhsib;
flxhr.hsb = da_hr.hsb_psl;
flxhr.hlc = yhlc;
flxhr.hlib = yhlib_lic;
flxhr.hlb = da_hr.hlb_psl;
flxhr.tauc = ytauc;
flxhr.taucx = ytaucx;
flxhr.tauib = ytauib;
flxhr.taub = da_hr.taub_psl;
flxhr.ct = yct;
flxhr.cq = ycql;
flxhr.cu = ycu;
flxhr.cw = ycw;
flxhr.jplume = qc_hr.jplume;
flxhr.jmanuv = qc_hr.jmanuv;
flxhr.tilt = da_hr.tilt;
flxhr.RF = da_hr.RF_psl;
flxhr.hlwebb = da_hr.hlwebb_psl;
flxhr.rlclr = yrlclr;
flxhr.rsclr = yrsclr;
flxhr.u10n = c35_hr.U10N;
flxhr.t10n = c35_hr.T10N;
flxhr.q10n = c35_hr.Q10N;
flxhr.usid_flag = yfilt_u;
flxhr.qsid_flag = yfilt_q;
flxhr.tsid_flag = yfilt_t;
flxhr.zu = da_hr.zu_psl;
flxhr.zt = da_hr.zt_psl;
flxhr.zq = da_hr.zq_psl;
flxhr.hed = hed_hr;
flxhr.sog = sog_hr;
flxhr.cog = cog_hr;
flxhr.lat = da_hr.lat_psl;
flxhr.lon = da_hr.lon_psl;
flxhr.wspd_ship = wspd_ship_hr;
flxhr.wdir_ship = wdir_ship_hr;
flxhr.ta_ship = da_hr.ta_ship;
flxhr.rh_ship = da_hr.rh_ship;
flxhr.sst_ship = da_hr.sst_ship;
flxhr.Pmb_ship = da_hr.Pmb_ship;
flxhr.qs_ship = da_hr.qs_ship;
flxhr.qa_ship = da_hr.qa_ship;
flxhr.rs_ship = da_hr.rs_ship;
flxhr.rl_ship = da_hr.rl_ship;
flxhr.q_lic = da_hr.q_lic;
flxhr.q_lic_std = sqrt(da_hr.qvar);
flxhr.co2_lic = da_hr.co2_lic_good;
flxhr.co2_lic_std = da_hr.co2_std_good;
flxhr.rnl = yrnl;
flxhr.hnetb = yhnet_C35;
flxhr.hnetc = yhnet_cov;
flxhr.hnetib = yhnet_ID;

% save in .mat format
dfl = fullfile(da_path,['FATIMA_2022_flux_hr_',dcname,'.mat']);
save(dfl,'flxhr');

%% save output in ascii format with variable name and unit header lines
% format header strings
hdr1 = 'FATIMA hourly flux Rev 1, July, 2022, BWB\n';
hdr2 = ''; hdr3 = '';
for ii = 1:length(f60_vars)
    n = f60_vars{ii}; ln = length(n); padn = repmat(' ',1,21-ln);
    hdr2 = [hdr2,n,padn];
    u = f60_units{ii}; lu = length(u); padu = repmat(' ',1,21-lu);
    hdr3 = [hdr3,u,padu];
end
hdr2 = [hdr2,'\n']; hdr3 = [hdr3,'\n'];

temp = table2array(flxhr);
% save results to ascii format with header lines
fout = fopen(fullfile(da_path,['FATIMA_2022_flux_hr_',dcname,'.txt']),'w');
fprintf(fout,hdr1);
fprintf(fout,hdr2);
fprintf(fout,hdr3);
fprintf(fout,[repmat('%20.10e ',1,nnn),' \n'],temp');
fclose all;


