function met = read_met_FATIMA_2022(dfl1,dfl2,dfl3,z_wxt,z_PTB,adj,gap_lim)
%{
reads met0 and wxt_parsed files and returns 1 min avg table with the
following variables:

    1   dnum            matlab datenum timestamp
    2   Tsea            sea snake surface SST
    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
    8   Td1             pir1 dome temperature
    9   Tc2             pir2 case temperature
    10  Td2             pir2 dome temperature
    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
    19  V_wxt           rel V wind component, WXT
    20  sst_rosr        radiometric surface skin temperature, °C

input parameters: dfl1 = string path to met0 file
                  dfl2 = string path to wxt0 file
                  dfl3 = path to rosr sst mat file
                  z_wxt = wxt height above SL, m
                  z_PTB = PTB330 height above SL, m
                  adj =  1x6 array of temperature adjustments
                        [tsea,td1,tc1,td2,tc2,wxt_t]
                  path_python: path to python folder with WXT parse script

    Mean 1 raw file format:
    1:3 PSD timecode: MMSSsss
    4   Batt_Volt
    5   PTemp
    6   PIR1 case R
    7   PIR1 case T, °C
    8   PIR1 dome R
    9   PIR1 dome T, °C
    10  PIR2 case R
    11  PIR2 case T, °C
    12  PIR2 dome R
    13  PIR2 dome T, °C
    14  PIR1 volts
    15  PIR1 W/m2, uncorrected
    16  PSP1 volts
    17  PSP1 W/m2
    18  PIR2 volts
    19  PIR2 W/m2, uncorrected
    20  PSP2 volts
    21  PSP2 W/m2
    22  Seasnake R
    23  Seasnake T, °C
    24  Pmb, Vaisala PTB330

    wxt_parsed file format:
    1:3 PSD timecode: MMSSsss
    4   wind direction, degrees
    5   wind speed, m/s
    6   air temperature, C
    7   RH, %
    8   pressure, mb or hPa
    9   cumulative rain, mm
    10  rain event duration, sec
    11  rain intensity, mm/hr
    12  cumulative hail, hits/cm^2
    13  hail event duration, s
    14  hail intensity, hits/cm^2/hr

    NOTE: wxt_parsed FILE IS GENERATED FROM RAW wxt0 FILES WITH PYTHON SCRIPT
%}

%% initialize

g = 0.0098;          % dry adiabatic lapse rate, C/m
sig_sb = 5.67e-8;    % Stefan Boltzmann constant
C2K = 273.15;        % temp conversion constant

nr1 = 24;     % number of variables in met0
nr2 = 14;     % number of variables in wxt0
mean0 = [];   % array for raw met0 data
wxt = [];
temp = {};

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

% reference timestamp and 1-min averaging bins for this hour
start = datenum(yyyy-1,12,31) + (ddd+hh/24);
delta = double(1.0/1440);
last = start + 59*delta;
dnum = start:delta:last; dnum = dnum(:);
dnum_bin = start:delta:last+delta; dnum_bin = dnum_bin(:);

%% read met0 file
if exist(dfl1,'file')==2
    disp(['Reading met0 file for hour ',int2str(hh)]);
    flist = fopen(dfl1,'r');
    while ~feof(flist)
        try
            temp = textscan(flist,['%2f%2f%3f ',repmat('%f ',1,21),'%*[^\n]'],'delimiter',', ','headerlines',1);
            mean0 = [mean0 cell2mat(temp)];
        catch
            [r,~] = size(temp{1,nr1});
            error(['read error ',fname,' line ',sprintf('%2i',r+2)]);
            % if one line is cut short, blank it out totally
%             [r,c] = size(temp{1,nr1}); % check size of last cell
%             if r==0 || c==0; continue; end % read next batch of lines
%             for ll=1:nr1
%                 if length(temp{1,ll})~=r
%                     temp{1,ll}(r+1)=[];
%                 end
%             end
%             for ii=1:length(temp)  % Look for NaNs and blanks in data fields
%                 temp{ii}(ismember(temp{ii},'NAN'))=cellstr('NaN');
%                 temp{ii}(ismember(temp{ii},''))=cellstr('NaN');
%                 mean0 = [mean0 str2num(char(temp{ii}))];
%             end
        end
    end
    dnum_1 = datenum(yyyy-1,12,31) + ddd + hh/24 + mean0(:,1)/1440 + ...
        (mean0(:,2) + mean0(:,3)/1000)/86400;   % timestamp from PSD system
    fclose(flist);
else
    mean0 = NaN(60,nr1);
    dnum_1 = dnum;
end

%% read wxt0 file
if exist(dfl2,'file')==2
    disp(['Reading wxt_parsed file for hour ',int2str(hh)]);
    flist = fopen(dfl2,'r');
    try
        temp = textscan(flist,['%2f%2f%3f ',repmat('%f ',1,11),'%*[^\n]'],'delimiter',',');
        wxt = [wxt cell2mat(temp)];
    catch
        [~,fname,~] = fileparts(dfl2);
        [r,~] = size(temp{1,nr2});
        error(['read error ',fname,' line ',sprintf('%2i',r+1)]);
    end

    [wxtU,wxtV] = sd_to_uv(wxt(:,5),wxt(:,4));  % compute wind components
    wxt = [wxt,wxtU,wxtV];                      % add to array

    dnum_2 = datenum(yyyy-1,12,31) + ddd + hh/24 + wxt(:,1)/1440 + ...
        (wxt(:,2) + wxt(:,3)/1000)/86400;  % timestamp from PSD system
    fclose(flist);
    wxt(:,6) = wxt(:,6) + adj(6);               % Ta WXT air temp adjustment
else
    wxt = NaN(60,nr2+2);
    dnum_2 = dnum';
end

%% format variables
% copy over original points to the nearest time on grid
press = NaN(length(dnum),1); psp1 = press; psp2 = press; Tsea = press;
Tc1 = press; Tc2 = press; Td1 = press; Td2 = press; therm1 = press; therm2 = press;
temp = dnum;
for ii = 1:length(dnum_1)                          % loop over all data points
    if any(isfinite(temp))                         % if there are unused grid points
        [~,jj] = nanmin1(abs(temp - dnum_1(ii)));  % find closest grid point for this value
        dT = (dnum_1(ii)-temp(jj))*86400;          % compute deltaT in Sec
        if abs(dT)<=30                             % only if dT is within range of +/- 30 sec...
            press(jj) = mean0(ii,24) + 0.125*z_PTB;% pressure (mbar)
            psp1(jj) = mean0(ii,17);               % solar radiometers
            psp2(jj) = mean0(ii,21);
            Tsea(jj) = mean0(ii,23) + adj(1);      % sea snake
            Tc1(jj) = mean0(ii,7) + adj(3);        % case temperature for PIR1 Eppley
            Td1(jj) = mean0(ii,9) + adj(2);        % dome temperature for PIR1 Eppley
            Tc2(jj) = mean0(ii,11) + adj(5);       % case temperature for PIR2 Eppley
            Td2(jj) = mean0(ii,13) + adj(4);       % dome temperature for PIR2 Eppley
            therm1(jj) = mean0(ii,15);             % PIR1 thermopile output, W/m2
            therm2(jj) = mean0(ii,19);             % PIR2 thermopile output, W/m2
            temp(jj) = NaN;                        % then remove this grid point from further consideration
        end
    end
end

%% despike mean met variables
press = despike(press);
psp1 = despike(psp1);
psp2 = despike(psp2);
Tsea = despike(Tsea);
Tc1 = despike(Tc1);
Td1 = despike(Td1);
Tc2 = despike(Tc2);
Td2 = despike(Td2);
therm1 = despike(therm1);
therm2 = despike(therm2);

%% compute IR radiation
pir1 = -therm1 + sig_sb*(Tc1+C2K).^4 - 3.63*sig_sb*((Td1+C2K).^4-(Tc1+C2K).^4);
pir2 = -therm2 + sig_sb*(Tc2+C2K).^4 - 3.37*sig_sb*((Td2+C2K).^4-(Tc2+C2K).^4);

%% average wxt to 1 min timestamp
if exist(dfl2,'file')==2
    wxt1 = interval_avg(dnum_2,wxt(:,4:16),dnum_bin);
    [spd,dir] = uv_to_sd(wxt1(:,13),wxt1(:,14));    % recompute wspd/wdir
    wxt1(:,2) = dir;
    wxt1(:,3) = spd;
    wxt1(:,6) = wxt1(:,6) + 0.125*z_wxt;
else
    wxt1 = NaN(60,16);
    wxt1(:,1) = dnum;
end

%% read rosr and interp to 1-min timestamp
rosr_sst = read_rosr_FATIMA_2022(dfl3,dnum);

%% format output table
mvars = {'dnum','Tsea','psp1','pir1','psp2','pir2','Tc1',...
         'Td1','Tc2','Td2','press','rwdir_wxt','rwspd_wxt','Tair_wxt',...
         'RH_wxt','Pmb_wxt','rain_wxt','U_wxt','V_wxt','rosr_sst'};

met = [dnum,Tsea,psp1,pir1,psp2,pir2,Tc1,Td1,Tc2,Td2,press,wxt1(:,2:6),...
       wxt1(:,9),wxt1(:,13:14),rosr_sst];

% fill short gaps in data of 5 min or less
% don't gap fill last column of rosr_sst for now
[~,c] = size(met);
for ii = 2:c-1
    jj = NaN_runs(met(:,ii),gap_lim);   % indices for NaN runs >= gap_lim
    met(:,ii) = fillgaps(met(:,ii),60); % fill NaN gaps, use 60 point window
    met(jj,ii) = NaN;                   % restore NaN to longer gaps
end

met = array2table(met,'VariableNames',mvars);

end
