function met = read_met_CAPRICORN_2016(dfl1,dfl2,ddd,hhh,adj)
%{
reads met1 and met2 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   psp1            psp1 radiation, W/m2
    6   pir1            pir1 radiation, W/m2
    7   psp2            psp2 radiation, W/m2
    8   pir2            pir2 radiation, W/m2
    9   Tc1             pir1 case temperature
    10  Td1             pir1 dome temperature
    11  Tc2             pir2 case temperature
    12  Td2             pir2 dome temperature
    13  org             rain, mm/hr
    14  press           air pressure, mb
    15  aspir_on        T/RH diagnostic
    16  org_carrier     org diagnostic

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

    Mean 1 raw file format:
    1:3 PSD timecode: MMSSsss
    4:6 Campbell mm/dd/yyyy
    7:9 Campbell hh:mm:ss
    10 Average(1,Batt_Volt,FP2,0)
    11 Average (1,PTemp,FP2,0)
    12 Average(1,ORG_Car_v,FP2,0)
    13 Average(1,ORG_Sig_v,FP2,0)
    14 Average(1,ORG_mmhr,FP2,0)
    15 Average(1, airflow, FP2,0)
    16 Sample (1,TRH(2),FP2): RH %
    17 Sample (1,TRH(1),FP2): Temp C

    Mean 2 raw file format:
    1:3 PSD timecode: MMSSsss
    4 record number
	5 Sample(1,Batt_Volt,FP2)
	6 Sample (1,PTemp,FP2)
	7 Average (1,Case1_R,IEEE4,0)
	8 Average (1,Case1_T,FP2,0)
	9 Average (1,Dome1_R,IEEE4,0)
	10 Average (1,Dome1_T,FP2,0)
	11 Average (1,Case2_R,IEEE4,0)
	12 Average (1,Case2_T,FP2,0)
	13 Average (1,Dome2_R,IEEE4,0)
	14 Average (1,Dome2_T,FP2,0)
	15 Average (1,PIR1_V,FP2,0)
	16 Average (1,PIR1_wm2,FP2,0)
	17 Average (1,PSP1_V,FP2,0)
	18 Average (1,PSP1_wm2,FP2,0)
	19 Average (1,PIR2_V,FP2,0)
	20 Average (1,PIR2_wm2,FP2,0)
	21 Average (1,PSP2_V,FP2,0)
	22 Average (1,PSP2_wm2,FP2,0)
	23 Average(1,Snake_r4w,IEEE4,0)
    24 Average(1,SnakeC_4W,IEEE4,0)
	25 Pressure, mb
%}

%% initialize

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

nr1 = 17;     % number of variables for met sys#1 in met1
nr2 = 25;     % number of variables for met sys#2 in met2
mean1 = [];   % array for raw met1 data
mean2 = [];   % array for raw met2 data
temp = {};

%% read met0 file
disp(['Reading met1 file for hour ',int2str(hhh)]);
flist = fopen(dfl1,'r');
while ~feof(flist)
    try
        temp = textscan(flist,['%2f%2f%3f %f/%f/%f %f:%f:%f ',repmat('%f ',1,8),'%*[^\n]'],'delimiter',', ','headerlines',1);
        mean1 = [mean1 cell2mat(temp)];
    catch
        % 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=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}))];
        end;
    end;
end;
mean1 = mean1';
jd_pc_1 = ddd+(hhh+(mean1(1,:)+(mean1(2,:)+mean1(3,:)/1000)/60)/60)/24;   % timestamp from PSD system
fclose(flist);

%% read me20 file
disp(['Reading met2 file for hour ',int2str(hhh)]);
flist = fopen(dfl2,'r');
while ~feof(flist)
    try
        temp = textscan(flist,['%2f%2f%3f ',repmat('%f ',1,22),'%*[^\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=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}))];  % 60xnr2 array
        end;
    end; %end try/catch
end;
mean2 = mean2';
jd_pc_2 = ddd+(hhh+(mean2(1,:)+(mean2(2,:)+mean2(3,:)/1000)/60)/60)/24;   % timestamp from PSD system
fclose(flist);

%% format variables
press = mean2(25,:);      % pressure (mbar)
psp1 = mean2(18,:);       % solar radiometers
psp2 = mean2(22,:);
org_carrier = mean1(12,:);
aspir_on = mean1(15,:);
org = mean1(14,:);
org(org<0) = 0;
Tsea = mean2(24,:) + adj(1);
Tc1 = mean2(8,:) + adj(4);   % case temperature for PIR1 Eppley
Td1 = mean2(10,:) + adj(3);  % dome temperature for PIR1 Eppley
Tc2 = mean2(12,:) + adj(6);  % case temperature for PIR2 Eppley
Td2 = mean2(14,:) + adj(5);  % dome temperature for PIR2 Eppley
therm1 = mean2(16,:);        % PIR thermopile output, W/m2
therm2 = mean2(20,:);
Rhvais = mean1(16,:);
Tvais = mean1(17,:) + adj(2);

%% compute IR 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);

% Select data to be output.
% pir = pir2;
% pir = (pir1+pir2)/2;    % use averages
% psp = (psp1+psp2)/2;

%% 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
start = (ddd+hhh/24);
delta = double(1.0/1440);
last = start + 59*delta;
jd_ref = start:delta:last;
met = zeros(60,16)*NaN;
met(:,1) = jd_ref;
met(:,2) = interp1(jd_pc_1',Tvais,jd_ref','nearest','extrap');
met(:,3) = interp1(jd_pc_1',Rhvais,jd_ref','nearest','extrap');
met(:,4) = interp1(jd_pc_2',Tsea,jd_ref','nearest','extrap');
met(:,5) = interp1(jd_pc_2',psp1,jd_ref','nearest','extrap');
met(:,6) = interp1(jd_pc_2',pir1,jd_ref','nearest','extrap');
met(:,7) = interp1(jd_pc_2',psp2,jd_ref','nearest','extrap');
met(:,8) = interp1(jd_pc_2',pir2,jd_ref','nearest','extrap');
met(:,9) = interp1(jd_pc_2',Tc1,jd_ref','nearest','extrap');
met(:,10) = interp1(jd_pc_2',Td1,jd_ref','nearest','extrap');
met(:,11) = interp1(jd_pc_2',Tc2,jd_ref','nearest','extrap');
met(:,12) = interp1(jd_pc_2',Td2,jd_ref','nearest','extrap');
met(:,13) = interp1(jd_pc_1',org,jd_ref','nearest','extrap');
met(:,14) = interp1(jd_pc_2',press,jd_ref','nearest','extrap');
met(:,15) = interp1(jd_pc_1',aspir_on,jd_ref','nearest','extrap');
met(:,16) = interp1(jd_pc_1',org_carrier,jd_ref','nearest','extrap');

