function met = read_met_SeaState_2015(dfl1,dfl2,ddd,hhh,adj)
%{
reads met0 and met1 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            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           air pressure, mb
    15  aspir_on        T/RH diagnostic
    16  org_carrier     org diagnostic
    17  Rh337           RH from heated 337 probe
    18  T337            T from 337 probe, C

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

%% initialize

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

nr1 = 19;     % number of variables for met sys#1 in met0
nr2 = 24;     % number of variables for met sys#2 in met1
mean1 = [];   % array for raw met1 data
mean2 = [];   % array for raw met2 data

%% read met0 file
temp = {};
disp(['Reading met0 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 %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';
jd_pc_1 = ddd+(hhh+(mean1(1,:)+(mean1(2,:)+mean1(3,:)/1000)/60)/60)/24;   % timestamp from PSD system
fclose(flist);

%% read met1 file
temp = {};
disp(['Reading met1 file for hour ',int2str(hhh)]);
flist = fopen(dfl2,'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 %*[^\n]','delimiter',', ','headerlines',1,'emptyvalue', NaN);
        mean2 = [mean2 cell2mat(temp)];
    catch
        % if one line is cut short, blank it out totally
        [l8,c8] = size(temp{1,nr2});
        if l8==0 || c8==0; continue; end; % read next batch of lines
        for ll=1:8
            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');
            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(24,:);      % pressure (mbar)
psp1 = mean2(17,:);       % solar radiometers
psp2 = mean2(21,:);
org_carrier = mean1(12,:);
aspir_on = mean1(15,:);
org = mean1(14,:);
Tsea = mean2(23,:) + adj(1);
Tc1 = mean2(7,:) + adj(4);   % case temperature for PIR1 Eppley
Td1 = mean2(9,:) + adj(3);   % dome temperature for PIR1 Eppley
Tc2 = mean2(11,:) + adj(6);  % case temperature for PIR2 Eppley
Td2 = mean2(13,:) + adj(5);  % dome temperature for PIR2 Eppley
therm1 = mean2(15,:);        % PIR thermopile output, W/m2
therm2 = mean2(19,:);
Rhvais = mean1(16,:);
Tvais = mean1(17,:) + adj(2);
Rh337 = mean1(18,:);
T337 = mean1(19,:) + adj(7);

%% Fix Tc1 based on date/time.  Prior to jd 277.6757 TC1 is high  by 1.57 C
% use Tc2 for this period...
ii = find(jd_pc_2 < 277.6757);
if ~isempty(ii)
    Tc1(ii) = Tc2(ii);
end

%% 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>110) = NaN;
Rhvais = replace_NaN_nearest_neighbor(Rhvais);
[Rhvais, N] = despike2(Rhvais);

Rh337(Rh337>110) = NaN;
Rh337 = replace_NaN_nearest_neighbor(Rh337);
[Rh337, N] = despike2(Rh337);

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

T337(T337>30) = NaN;
T337 = replace_NaN_nearest_neighbor(T337);
[T337, N] = despike2(T337);

press = replace_NaN_nearest_neighbor(press);
[press,N] = 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');
met(:,17) = interp1(jd_pc_1',Rh337,jd_ref','nearest','extrap');
met(:,18) = interp1(jd_pc_1',T337,jd_ref','nearest','extrap');

