function met = read_met_DYNAMO_2011(dfl,ddd,hhh,adj)
%{
reads met0 and me20 raw 1-min mean files and returns array 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   psp             psp radiation, W/m2
6   pir             pir radiation, W/m2
7   Td              pir1 dome temperature
8   org             rain, mm/hr
9   press           Air pressure, mb
10  aspir_on        T/RH diagnostic
11  org_carrier     org diagnostic

input parameters: dfl = string path to met0 file
                  ddd = julian date
                  hhh = hour
                  adj =  array of temperature adjustments
                        [tsea_adj,ta_psd_adj,rh_psd_adj,bp_psd_adj,td1_adj,tc1_adj,td2_adj,tc2_adj]
                             1        2          3          4         5       6       7       8
%}

%% initialize

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

nr1 = 65;     % number of variables for met sys#1 in met0
mean1 = [];   % array for raw me10 data
temp = {};

%% read met0 file
disp(['Reading met0 file for hour ',int2str(hhh)]);
flist = fopen(dfl,'r');
while ~feof(flist)
    try
        temp =textscan(flist,'%2f%2f%3f %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 %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;
% adjust MET for 2min lag issue
mean1 = mean1';
times = mean1(1:7,:);
data = [mean1(8:end,3:end),ones(58,2)*NaN]; % shift data back 2 minutes...
mean1 = [times;data];
jd_pc = ddd + (hhh+(mean1(1,:) + (mean1(2,:)  +mean1(3,:)/1000)/60)/60)/24;   % timestamp from PSD system
fclose(flist);

%% format variables
press = mean1(36,:) + adj(4); % pressure (mbar)
psp1 = mean1(30,:);           % solar radiometers
psp2 = mean1(31,:);
org_carrier = NaN(size(jd_pc));
aspir_on = NaN(size(jd_pc));
org = mean1(32,:)/0.7;        % adjustment based on comparison with UConn rain gauge
Tsea = mean1(25,:) + adj(1);
Tc1 = mean1(21,:) + adj(6);   % case temperature for PIR1 Eppley
Td1 = mean1(22,:) + adj(5);   % dome temperature for PIR1 Eppley
Tc2 = mean1(23,:) + adj(8);   % case temperature for PIR2 Eppley
Td2 = mean1(24,:) + adj(7);   % dome temperature for PIR2 Eppley
therm1 = mean1(28,:);         % PIR thermopiles
therm2 = mean1(29,:);
Rhvais = mean1(27,:) + adj(3);
Tvais = mean1(26,:) + adj(2);

%% compute radiation
pir1 = therm1 + sig_sb*(Tc1+C2K).^4 - 3*sig_sb*((Td1+C2K).^4 - (Tc1+C2K).^4); % adjust B coeff to 3 for leg1
pir2 = therm2 + sig_sb*(Tc2+C2K).^4 - 4*sig_sb*((Td2+C2K).^4 - (Tc2+C2K).^4);

% PIR1 had solar contamination or overheating so use PIR2
pir = pir2;
% average psp
psp = (psp1+psp2)/2;

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

Tvais(Tvais>35) = 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
start = (ddd+hhh/24);
delta = double(1.0/1440);
last = start + 59*delta;
jd_ref = [start:delta:last];
met = zeros(60,11)*NaN;
met(:,1) = jd_ref;
met(:,2) = interp1(jd_pc',Tvais,jd_ref','nearest');
met(:,3) = interp1(jd_pc',Rhvais,jd_ref','nearest');
met(:,4) = interp1(jd_pc',Tsea,jd_ref','nearest');
met(:,5) = interp1(jd_pc',psp,jd_ref','nearest');
met(:,6) = interp1(jd_pc',pir,jd_ref','nearest');
met(:,7) = interp1(jd_pc',Td1,jd_ref','nearest');
met(:,8) = interp1(jd_pc',org,jd_ref','nearest');
met(:,9) = interp1(jd_pc',press,jd_ref','nearest');
met(:,10) = interp1(jd_pc',aspir_on,jd_ref','nearest');
met(:,11) = interp1(jd_pc',org_carrier,jd_ref','nearest');

% NaN last two minutes of data as missing values due to timestamp issue...
met(59:60,2:end) = NaN;

