function [met1, met10] = read_met_WHOTS_2022(dfl1,dfl2,dfl3,ddd,hh,adj)
%{
    reads met2 and wxp1/wxp2 files and returns 1 min and 10 min avg tables for the
    hour with the following variables:

        1   jd_met          1 min timestamp, decimal DOY @ start of interval
        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             vaisala pressure from mean2 box (no height adjustment)
        12  rwdir_wxt1      rel wind direction, WXT1, 0-360 deg from bow
        13  rwspd_wxt1      rel wind speed, WXT1, m/s
        14  Tair_wxt1       air temperature, WXT1, C
        15  RH_wxt1         RH, WXT1, %
        16  Pmb_wxt1        air pressure, WXT1, mb or hPa
        17  rain_wxt1       rain rate, WXT1, mm/hr
        18  U_wxt1          U (bow-stern) wind component/0.95, WXT1
        19  V_wxt1          V (port-stbd) wind component/1.15, WXT1
        20  rwdir_wxt2      rel wind direction, WXT2, 0-360 deg from bow
        21  rwspd_wxt2      rel wind speed, WXT2, m/s
        22  Tair_wxt2       air temperature, WXT2, C
        23  RH_wxt2         RH, WXT2, %
        24  Pmb_wxt2        air pressure, WXT2, mb or hPa
        25  rain_wxt2       rain rate, WXT2, mm/hr
        26  U_wxt2          U (bow-stern) wind component/0.95, WXT2
        27  V_wxt2          V (port-stbd) wind component/1.15, WXT2


    input parameters: dfl1 = string path to met2 file
                      dfl2 = string path to wxt1 file
                      dfl3 = string path to wxt2 file
                      ddd = julian date
                      hh = hour
                      adj =  1x7 array of temperature adjustments
                            [tsea,tair1,tair2,td1,tc1,td2,tc2]

        Mean 2 raw file format:
         1:3    PSD timecode: MMSSsss
         4      Batt_Volt
         5      PTemp
         6      Case1_R_Avg
         7      Case1_T_Avg
         8      Dome1_R_Avg
         9      Dome1_T_Avg
         10     Case2_R_Avg
         11     Case2_T_Avg
         12     Dome2_R_Avg
         13     Dome2_T_Avg
         14     PIR1_V_Avg
         15     PIR1_wm2_Avg
         16     PSP1_V_Avg
         17     PSP1_wm2_Avg
         18     PIR2_V_Avg
         19     PIR2_wm2_Avg
         20     PSP2_V_Avg
         21     PSP2_wm2_Avg
         22     Snake_R_Med
         23     Snake_C_Med
         24     hpat_Avg

        wxt 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
        15  heater temperature, C
        16  heater voltage

        NOTE: wxpN*.txt FILE IS GENERATED FROM RAW wxtN*.txt FILES WITH PYTHON SCRIPT
%}

%% initialize
sig_sb = 5.67e-8;    % Stefan Boltzmann constant
C2K = 273.15;        % temp conversion constant

nr2 = 24;     % number of raw fields in met2 file, as above
mean2 = [];   % empty array for raw met2 data
temp = {};

start = (ddd+hh/24); stop = start + 1/24;
delta1 = double(1.0/1440);
delta10 = double(1.0/144);
bins1 = start:delta1:stop; bins1 = bins1(:);
tref1 = bins1(1:60);
bins10 = start:delta10:stop; bins10 = bins10(:);

met1m = NaN(60,27);
met1m(:,1) = tref1;
met10m = NaN(6,27);
met10m(:,1) = bins10(1:6);

%% read met2 file
disp(['Reading met 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);
        mean2 = [mean2; cell2mat(temp)];
    catch
        % if one line is cut short, blank it out totally
        [r,c] = size(temp{1,nr2});
        if r==0 || c==0; continue; 
        end % read next batch of lines
        for ll=1:nr2
            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');
% %             mean2 = [mean2; str2num(char(temp{ii}))];  % 60x nr2 array
%         end
        mean2 = [mean2; cell2mat(temp)];
    end
end
fclose(flist);

% compute timestamp for campbell logger output
jd_met2 = ddd+(hh+(mean2(:,1)+(mean2(:,2)+mean2(:,3)/1000)/60)/60)/24;   % timestamp from PSD system
% subtract one minute from campbell timestamp, since it is at the end of averaging interval
% this means the last min of data this hour is actually the first point from next hour
% so, the output from this routine will be missing the last point each hour in the 1m array
jd_met2 = jd_met2 - 1/1440;
% and round time to nearest minute
jd_met2 = round(jd_met2*1440)/1440;

%% temperature adjustments
ii = mean2(:,23)<0;  mean2(ii,23) = NaN;  % remove bad SST (snake unplugged)
mean2(:,23) = mean2(:,23) + adj(1);  % sea snake
mean2(:,7) = mean2(:,7) + adj(5);    % case temperature for PIR1
mean2(:,9) = mean2(:,9) + adj(4);    % dome temperature for PIR1
mean2(:,11) = mean2(:,11) + adj(7);  % case temperature for PIR2
mean2(:,13) = mean2(:,13) + adj(6);  % dome temperature for PIR2

%% compute IR radiative flux (downwelling, so sign is negative)
pir1 = -1*(mean2(:,15) + sig_sb*(mean2(:,7)+C2K).^4 - 4*sig_sb*((mean2(:,9)+C2K).^4 - (mean2(:,7)+C2K).^4));
pir2 = -1*(mean2(:,19) + sig_sb*(mean2(:,11)+C2K).^4 - 4*sig_sb*((mean2(:,13)+C2K).^4 - (mean2(:,11)+C2K).^4));

%% format psd met array [jd,ts,psp1,pir1,psp2,pir2,tc1,td1,tc2,td2,Pmb]
% negate PIR values here so heat flux is w/respect to ocean
met2 = [jd_met2,mean2(:,23),mean2(:,17),-pir1,mean2(:,21),-pir2,mean2(:,7),mean2(:,9),...
         mean2(:,11),mean2(:,13),mean2(:,24)];

%% read wxt files
disp(['Reading wxp1 file for hour ',int2str(hh)]); 
flist = fopen(dfl2,'r');
wxt2=[];
while ~feof(flist)
    try
        temp = textscan(flist,['%2f%2f%3f ',repmat('%f ',1,13),'%*[^\n]'],'delimiter',',');
        wxt1 = cell2mat(temp);
    catch
        % if one line is cut short, blank it out totally
            [r,c] = size(temp{1,16});
            if r==0 || c==0; continue; 
            end % read next batch of lines
            for ll=1:16
                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');
                wxt1 = [wxt1 cell2mat(temp)];  % 60x nr2 array
                dump=fgetl(flist);
%             end
    end
end
fclose(flist);


wxt1 = cell2mat(temp);
wxt1(:,6) = wxt1(:,6) + adj(2);  % air temp adjustment
[wxtU,wxtV] = sd_to_uv(wxt1(:,5),wxt1(:,4));  % compute wind components
% add UV cols 17,18 to array, adjust for flow distortion effects
wxt1 = [wxt1,wxtU/0.95,wxtV/1.15];
jd_wxt1 = ddd+(hh+(wxt1(:,1)+(wxt1(:,2)+wxt1(:,3)/1000)/60)/60)/24;   % timestamp from PSD system

% format wxt1 array as [jd,ta,rh,Pmb,rain_rate,U,V]
wxt1 = [jd_wxt1,wxt1(:,6:8),wxt1(:,11),wxt1(:,17:18)];

disp(['Reading wxp2 file for hour ',int2str(hh)]); 
flist = fopen(dfl3,'r');
wxt2=[];
while ~feof(flist)
    try
        temp = textscan(flist,['%2f%2f%3f ',repmat('%f ',1,13),'%*[^\n]'],'delimiter',',');
        wxt2 = cell2mat(temp);
    catch
        % if one line is cut short, blank it out totally
            [r,c] = size(temp{1,16});
            if r==0 || c==0; continue; 
            end % read next batch of lines
            for ll=1:16
                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');
                wxt2 = [wxt2 cell2mat(temp)];  % 60x nr2 array
                dump=fgetl(flist);
%             end
    end
end
fclose(flist);


wxt2(:,6) = wxt2(:,6) + adj(2);  % air temp adjustment
[wxtU,wxtV] = sd_to_uv(wxt2(:,5),wxt2(:,4));  % compute wind components
%% adjusts wxt2 by 5deg offset compared to WXT1 (see install pictures).
wxtU=wxtU.*cos(deg2rad(-5))-wxtV.*sin(deg2rad(-5));
wxtV=wxtV.*cos(deg2rad(-5))+wxtU.*sin(deg2rad(-5));
% add UV cols 17,18 to array, adjust for flow distortion effects
wxt2 = [wxt2,wxtU/0.95,wxtV/1.15];
jd_wxt2 = ddd+(hh+(wxt2(:,1)+(wxt2(:,2)+wxt2(:,3)/1000)/60)/60)/24;   % timestamp from PSD system

% format wxt1 array as [jd,ta,rh,Pmb,rain_rate,U,V]
wxt2 = [jd_wxt2,wxt2(:,6:8),wxt2(:,11),wxt2(:,17:18)];

%% prepare 1m average output

% copy met2 data into output array
for ii = 1:60
    idx = find(abs(met2(:,1)-met1m(ii,1))<1/86400); % to nearest second
    if isempty(idx); continue; end
    if length(idx)>1
        met1m(ii,2:11) = met2(idx(1),2:11);
    else
        met1m(ii,2:11) = met2(idx,2:11);
    end
end

% average wxt arrays and compute rwspd/rwdir
% wind is adjusted for flow distortion!!
temp1 = interval_avg(wxt1(:,1),wxt1(:,2:end),bins1);
[spd1,dir1] = uv_to_sd(temp1(:,6),temp1(:,7));
temp2 = interval_avg(wxt2(:,1),wxt2(:,2:end),bins1);
[spd2,dir2] = uv_to_sd(temp2(:,6),temp2(:,7));

% merge to output array
met1m(:,12:19) = [spd1,dir1,temp1(:,2:7)];
met1m(:,20:27) = [spd2,dir2,temp2(:,2:7)];

% output table
vars = {'jd','Tsea','psp1','pir1','psp2','pir2','Tc1','Td1','Tc2','Td2','Pmb',...
        'rwspd_wxt1','rwdir_wxt1','Tair_wxt1','RH_wxt1','Pmb_wxt1','rain_wxt1','U_wxt1',...
        'V_wxt1','rwspd_wxt2','rwdir_wxt2','Tair_wxt2','RH_wxt2','Pmb_wxt2','rain_wxt2',...
        'U_wxt2','V_wxt2'};
met1 = array2table(met1m,'VariableNames',vars);

%% prepare 10m average output as above

% average met and wxt arrays and compute rwspd/rwdir
% wind is adjusted for flow distortion!!
temp1 = interval_avg(met2(:,1),met2(:,2:end),bins10);
temp2 = interval_avg(wxt1(:,1),wxt1(:,2:end),bins10);
[spd1_10,dir1_10] = uv_to_sd(temp2(:,6),temp2(:,7));
temp3 = interval_avg(wxt2(:,1),wxt2(:,2:end),bins10);
[spd2_10,dir2_10] = uv_to_sd(temp3(:,6),temp3(:,7));

% add std deviation of 1-min rwdir for both wxt units
rw1 = unwrap(dir1*pi/180)*180/pi;
rw2 = unwrap(dir2*pi/180)*180/pi;
temp4 = interval_std(met1m(:,1),rw1,bins10);
temp5 = interval_std(met1m(:,1),rw2,bins10);

met10m(:,2:29) = [temp1(:,2:11),spd1_10,dir1_10,temp2(:,2:7),spd2_10,dir2_10,temp3(:,2:7),...
                  temp4(:,2),temp5(:,2)];

% output table
vars10 = [vars,'std_rwdir1','std_rwdir2'];
met10 = array2table(met10m,'VariableNames',vars10);

end
