function [lic,dT_hist,dT_edges,missing] = read_licor_FATIMA_2022(dfl,read_raw)
%{
    reads licor file specified by dfl and returns a table, lic, with the
    following columns:
    
    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

    Also returns histogram counts and edges for delta T between
    data points and the number of missing points per hour.
    
    input parameters: dfl = string path to lic file

    raw format: MMSSsss AGC co2 h2o Ta P_hPa
    07/03/2022 08:00 Licor 7500 (SN: 0772)
    0000071 249	26.9842	1099.435	21.73	100.8	
    0000041 249	26.9816	1099.876	21.73	100.9	
    0000130 249	26.9819	1099.682	21.72	100.9
%}

fclose all;

[fpath,fname,~] = fileparts(dfl);
YY = fname(5:6); yyyy = 2000 + str2double(YY);
DDD = fname(7:9); ddd = str2double(DDD);
HH = fname(10:11);  hh = str2double(HH);

dflm = [fname,'.mat'];
% load lic mat table if it exists, unless force read_raw is indicated
if ~read_raw && exist(fullfile(fpath,dflm),'file')
    disp(['Reading Licor file ',fname]);
    load(fullfile(fpath,dflm),'lic','dT_hist','dT_edges','missing');
    return
end

%% 10 Hz reference timestamp and 1 Hz averaging bins
start = datenum(yyyy-1,12,31) + ddd + hh/24;
delta = double(1.0/864000);
last = start + 35999*delta;
dnum = start:delta:last;  dnum = dnum(:); % exact hourly 10Hz timestamp
dnum_bin = start:double(1.0/86400):last+double(1.0/86400); % 1 sec bins
dnum_bin = dnum_bin(:);
dnum_bin_10min = start:10/1440:start+1/24; dnum_bin_10min = dnum_bin_10min(:); % 10 min

missing = ones(6,1)*36000/6; %initialize to all missing
dT_hist = NaN(1,64);
dT_edges = NaN(1,65);

%% read data if file exists, if not return NaNs
if exist(dfl,'file')==2
    %% read data
    disp(['Reading Licor file ',fname]);
    flist = fopen(dfl,'r');
    temp = {};
    z = [];

    while feof(flist)==0
        try    % read file into cell array. temp is 1x8 cell array
            temp = textscan(flist,'%2f%2f%3f %f %f %f %f %f %*[^\n]','delimiter', ', ','headerlines', 1,'emptyvalue',NaN,'treatAsEmpty','Overflow');
            z = [z cell2mat(temp)'];
        catch
            for ii=1:10  % length of last cell reflects missing values from any/all fields
                if length(temp{1,ii})~=length(temp{1,10})
                    temp{1,ii}(length(temp{1,10})+1)=[]; %truncate length of all cells to length(temp{1,10})
                end
            end
            if ~isempty(temp{1,1})
                z = [z cell2mat(temp)'];
            end
        end
    end
    fclose(flist);
    z(:,any(isnan(z),1))=[]; % kill lines beginning with NaN
    z = z';
    
    % raw datemun timestamp
    tlic =  start + z(:,1)/1440 + (z(:,2) + z(:,3)/1000)/86400;
    % histogram of delta T for the raw data, should all be 0.1 sec
    [dT_hist,dT_edges] = histcounts(diff(tlic)*86400,-0.2:0.05:3);
    % count missing points in each 10-min interval
    % should be 6000 points per 10min
    for ii=1:6
        missing(ii) = 6000 - numel(tlic(tlic>=dnum_bin_10min(ii) & tlic<dnum_bin_10min(ii+1)));
    end

    Licor_diag = z(:,4);
    co2 = z(:,5);         % mmol/m^3
    h2o = z(:,6);         % mmol/m^3
    Ta_lic = z(:,7);      % °C
    Pmb_lic = z(:,8)*10;  % kPa to mb

    %% convert CO2 units from mol/m3 to best estimate of umol/mol or ppm
    Rgas_universal = 8.314472; % Pa m3 /K /mol
    co2 = ((co2*1e-3).*Rgas_universal.*(Ta_lic+273.15)./(Pmb_lic*1e2))*1e6; % mmol/m3 to umol/mol (ppm)

    %% interpret status diagnostic
    xx = Licor_diag;
    b7 = floor(xx/128);
    y = b7*128;
    b6 = floor((xx-y)/64);
    y = y+b6*64;
    b5 = floor((xx-y)/32);
    y = y+b5*32;
    b4 = floor((xx-y)/16);
    y = y+b4*16;
    AGC = (xx-y)*6.25;

    %% Format output to exactly 10 Hz
    lic = NaN(36000,6);
    lic(:,1) = dnum;

    % interpolate to 10 Hz
%     [~,ii,~] = unique(tlic);
%     lic(:,2) = interp1(tlic(ii),h2o(ii),tref,'linear','extrap');
%     lic(:,3) = interp1(tlic(ii),co2(ii),tref,'linear','extrap');
%     lic(:,4) = interp1(tlic(ii),Tlic(ii),tref,'linear','extrap');
%     lic(:,5) = interp1(tlic(ii),Plic(ii),tref,'linear','extrap');
%     lic(:,6) = interp1(tlic(ii),AGC(ii),tref,'linear','extrap');

    % or, copy over original points to the 10 Hz grid - this is slower
    temp = dnum;
    for ii = 1:length(tlic)                          % loop over all data points
        if any(isfinite(temp))                       % if there are unused grid points
            [~,jj] = nanmin1(abs(temp - tlic(ii)));  % find closest grid point for this datum
            dT = (tlic(ii)-temp(jj))*86400*1000;     % compute deltaT in mSec
            if dT>=-90 && dT<90                      % only if dT is within range of +/- 90 ms...
                lic(jj,2) = h2o(ii);                 % copy data over to gridded array point
                lic(jj,3) = co2(ii);
                lic(jj,4) = Ta_lic(ii);
                lic(jj,5) = Pmb_lic(ii);
                lic(jj,6) = AGC(ii);
                temp(jj) = NaN;            % then remove this grid point from further consideration
            end
        end
    end

    % fill in NaNs with nearest neighbor
%     lic(:,2) = replace_NaN_nearest_neighbor(lic(:,2));
%     lic(:,3) = replace_NaN_nearest_neighbor(lic(:,3));
%     lic(:,4) = replace_NaN_nearest_neighbor(lic(:,4));
%     lic(:,5) = replace_NaN_nearest_neighbor(lic(:,5));
%     lic(:,6) = replace_NaN_nearest_neighbor(lic(:,6));

else
    lic = NaN(36000,6);
    lic(:,1) = dnum;
end

vars = {'dnum','h2o','co2','Ta_lic','Pmb_lic','agc'};
units = {'none','mmol/m3','ppm','degC','mb','none'};
lic = array2table(lic,'VariableNames',vars);
lic.Properties.VariableUnits = units;

% save as .mat file for quicker loading next time
save(fullfile(fpath,dflm),'lic','dT_hist','dT_edges','missing');

end
