function met = read_met_HiWinGS_2013_2(dfl1,dfl2,ddd,hhh,adj)
%{
reads met0 and trh0 raw files and returns array of 1-min means with the
following columns:

    1   jd_met          1 min timestamp
    2   Tvais           T vaisala, C
    3   Rhvais          RH vaisala
    4   Tsea            sea snake surface SST
    5   psp1            psp radiation, W/m2
    6   pir1            pir radiation, W/m2
    7   psp2            psp radiation, W/m2
    8   pir2            pir radiation, W/m2
    9   Tc1             pir1 case temperature
    10  Td1             pir1 dome temperature
    11  Tc2             pir1 case temperature
    12  Td2             pir1 dome temperature
    13  org             rain, mm/hr
    14  press           atm pressure, mb
    15  aspir_on        T/RH diagnostic
    16  org_carrier     org diagnostic

input parameters: dfl1 = string path to met0 file
                  dfl2 = string path to me20 file
                  ddd = julian date
                  hhh = hour
                  adj =  array of temperature adjustments
                        [tsea,tair,td1,tc1,td2,tc2]

Input met0 file format:

%}

%% initialize

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

nr1 = 61;     % number of variables in met0
nr2 = 6;      % number of variables in trh0

start = (ddd+hhh/24);
delta = double(1.0/1440);
last = start + 59*delta;
jd_ref = [start:delta:last]; % reference 1-min timestamp for interp

%% read met0 file
mean1 = [];   % array for raw met1 data
temp = {};
disp(['Reading met0 file for hour ',int2str(hhh)]);
flist = fopen(dfl1,'r');
if flist>0
    while ~feof(flist)
        try
            %                       1                                10                            20                            30                            40                            50                            60
            temp = textscan(flist,'%2f%2f%3f %2f/%2f/%4f %2f:%2f:%2f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %*[^\n]','delimiter',', ','headerlines',1,'emptyvalue',NaN);
            mean1 = [mean1 cell2mat(temp)];
        catch
            % if one line is cut short, blank it out totally
            [l8,c8] = size(temp{1,nr1}); % check size of last cell
            if l8==0 || c8==0; continue; end; % read next batch of lines
            for ll=1:9
                if length(temp{1,ll})~=l8
                    temp{1,ll}(l8+1)=[];
                end;
            end;
            for ii=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');
                mean1=[mean1 str2num(char(temp{ii}))];  % 60xnr1 array
            end;
        end;
    end;
    mean1 = mean1';
    % shift data back 1 minute so timestamp is at beginning of averaging interval
    % data from last minute is duplicated to complete time series
    times = mean1(1:3,:);
    data = [mean1(4:end,2:end), mean1(4:end,end)];
    mean1 = [times;data];
    jd_pc_1 = ddd+(hhh+(mean1(1,:)+(mean1(2,:)+mean1(3,:)/1000)/60)/60)/24;   % timestamp from PSD system
    jd_pc_1 = round(jd_pc_1*1440)/1440; % to nearest minute
    fclose(flist);
else
    mean1 = NaN(61,60);
    jd_pc_1 = jd_ref;
end;

%% read trh0 file - this is 1Hz data
mean2 = [];   % array for raw met2 data
temp = {};
disp(['Reading trh0 file for hour ',int2str(hhh)]);
flist = fopen(dfl2,'r');
if flist>0
    while ~feof(flist)
            % read lines until eof or error
            temp = textscan(flist,'%2f%2f%3f %f %f %f %*[^\n]','headerlines',1,'emptyvalue', NaN);
            % if one line is cut short due to error, discard it
            [r1,c1] = size(temp{1,1});
            if r1==0 || c1==0; % no data, read next batch of lines
                continue;
            end;
            test = 1;
            for ii = 2:nr2  % see if variables are all the same length
                test = test * (length(temp{1,ii})==r1);
            end;
            if ~test        % if they aren't , delete last line read from temp
                for jj=1:nr2
                    temp{jj} = temp{jj}(1:r1-1);
                end;
            end;
            [r1,c1] = size(temp{1,1});
            if r1==0 || c1==0; % no more data, read next batch of lines
                continue;
            end;
            if any(any(imag(cell2mat(temp))))
                continue;  % skip this read if imaginary numbers pop up...
            end;
            if r1==1
                mean2 = [mean2;cell2mat(temp)]; % concat - only one line read
            else
                mean2 = [mean2 cell2mat(temp)]; % concat - multiple lines read
            end;
     end;
    fclose(flist);
    jd_pc_2 = ddd+(hhh+(mean2(:,1)+(mean2(:,2)+mean2(:,3)/1000)/60)/60)/24;   % timestamp from PSD system
    jd_bins = [jd_ref, jd_ref(end) + double(1/1440)];
    mean_trh = interval_avg(jd_pc_2,mean2(:,4:6),jd_bins); % avg to jd_ref timestamp
    mean_trh = mean_trh';
else
    mean_trh = NaN(4,60);
    mean_trh(1,:) = jd_ref';
end;

%% format variables
press = mean1(35,:);      % pressure (mbar)
psp1 = -mean1(25,:);      % solar radiometers
psp2 = -mean1(27,:);
org_carrier = mean1(28,:);
aspir_on = ones(size(press)); % not logged on HiWinGS
org = mean1(30,:);
Tsea = mean1(34,:) + adj(1);
Tc1 = mean1(13,:) + adj(4);   % case temperature for PIR1 Eppley
Td1 = mean1(15,:) + adj(3);   % dome temperature for PIR1 Eppley
Tc2 = mean1(17,:) + adj(6);   % case temperature for PIR2 Eppley
Td2 = mean1(19,:) + adj(5);   % dome temperature for PIR2 Eppley
therm1 = mean1(21,:);         % PIR thermopiles
therm2 = mean1(23,:);
Rhvais = mean_trh(2,:);
Tvais = mean_trh(3,:) + adj(2);

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

%% filter out bad values in T, RH and P
Rhvais(Rhvais>200) = NaN;
Rhvais = replace_NaN_nearest_neighbor(Rhvais);
[Rhvais,~] = despike2(Rhvais);

Tvais(Tvais>30) = NaN;
Tvais = replace_NaN_nearest_neighbor(Tvais);
[Tvais,~] = despike2(Tvais);

press = replace_NaN_nearest_neighbor(press);
[press,~] = despike2(press);

%% interp to exact 1-min timestamp
met = zeros(60,16)*NaN;
met(:,1) = jd_ref;
met(:,2) = Tvais;
met(:,3) = Rhvais;
met(:,4) = interp1(jd_pc_1',Tsea,jd_ref','nearest','extrap');
met(:,5) = interp1(jd_pc_1',psp1,jd_ref','nearest','extrap');
met(:,6) = interp1(jd_pc_1',pir1,jd_ref','nearest','extrap');
met(:,7) = interp1(jd_pc_1',psp2,jd_ref','nearest','extrap');
met(:,8) = interp1(jd_pc_1',pir2,jd_ref','nearest','extrap');
met(:,9) = interp1(jd_pc_1',Tc1,jd_ref','nearest','extrap');
met(:,10) = interp1(jd_pc_1',Td1,jd_ref','nearest','extrap');
met(:,11) = interp1(jd_pc_1',Tc2,jd_ref','nearest','extrap');
met(:,12) = interp1(jd_pc_1',Td2,jd_ref','nearest','extrap');
met(:,13) = interp1(jd_pc_1',org,jd_ref','nearest','extrap');
met(:,14) = interp1(jd_pc_1',press,jd_ref','nearest','extrap');
met(:,15) = interp1(jd_pc_1',aspir_on,jd_ref','nearest','extrap')*NaN;
met(:,16) = interp1(jd_pc_1',org_carrier,jd_ref','nearest','extrap');

