function licm = read_licor_day_STRATUS_2017(path_working_ddd,ddd,yyyy,...
                flicor,prtit,path_raw_images,graphformat,graphdevice)
%{
Reads hourly licor files for one day
Loops through 24 hours, calling read_licor_STRATUS_2017 for each hour
Raw 10Hz data is averaged to 1Hz output

Inputs: path_working_ddd: path to daily folder containing hourly files
        ddd: day-of-year variable
        yyyy: year string

Output: licm: 86400 x 10 array of 1-Hz avg gps data
        Returns NaN for data gaps

        1   jd_ref      10 Hz timestamp
        2   Licor_H2O   H2O, mmol/m3
        3   Licor_CO2   CO2, approx mmol/m3
        4   Licor_CO2a  CO2, approx mol/mol
        5   Licor_T     Licor box temperature
        6   Licor_P     Licor box pressure kPa
        7   AGC         Licor diagnostic
        8   Licor_mr    H2O, mol/mol
        9   Rho_air     Air density, kg/m3
        10  Licor_q     H2O, g/kg
%}

fclose all;

licm = zeros(86400,10)*NaN;
delta = double(1.0/86400);
last = ddd + 86400*delta;
jd_ref = ddd:delta:last;  jd_ref = jd_ref(:); % ref 1 Hz timestamp
licm(:,1) = jd_ref(1:end-1)';
jd = sprintf('%03i',ddd);
Rgas = 287.1;
Rgas_universal = 8.314472; % Pa m3 /K /mol
Mw = 18.01528;  % molar mass of H20 g/mol
Md = 28.964;    % molar mass dry air
epsilon = Mw/Md; % mass of water to mass of dry air = 0.622

for hhh = 0:23              % cycle thru 24 hourly gprm files
    hr = sprintf('%02i',hhh);
    dfl = fullfile(path_working_ddd,['lic0' yyyy(3:4),jd,hr,'_raw.txt']);
    if exist(dfl,'file')==2
        lic = read_licor_STRATUS_2017(dfl,ddd,hhh);
    else
        continue % no data this hour
    end
    jdlic = lic(:,1);

    % convert to mol/mol and specific humidity
    h2o_mr = (lic(:,2)*1e-3).*Rgas_universal.*(lic(:,5)+273.15)./(lic(:,6)*1e3); % mmol/m3 to mol/mol
    rhoa_lic = (lic(:,6)*1e3).*(1-h2o_mr*(1-epsilon))./(Rgas.*(lic(:,5)+273.15));% air density [kg/m3]
    h2o_q = (lic(:,2)*1e-3).*Mw./rhoa_lic;  % mmol/m3 to g/kg
    lic = [lic,h2o_mr,rhoa_lic,h2o_q];      % add columns to lic array

    % average this hour's data into 1Hz licm array
    start = ddd + hhh/24.0; % jd start time this hour
    diff = jd_ref - start;  % look for closest time stamp to start
    [~,ii] = min(abs(diff));
    temp = interval_avg(jdlic, lic(:,2:10), jd_ref(ii:ii+3600));
    licm(ii:ii+3599,2:10) = temp(:,2:10);

    if prtit && mod(hhh,6)==0 % spectra plot
        Npts = length(h2o_q);
        [Pqq,~] = psd2(detrend(h2o_q),Npts,flicor,hamming(Npts));
        [Pqq,Fxx] = specsmoo(Pqq,flicor);

        figure;loglog(Fxx,Fxx.*Pqq,'b-');
        hold on;
        loglog(Fxx(17:48),((Fxx(17:48)).^(-2/3))./(Fxx(35).^(-2/3)).*Pqq(35).*Fxx(35),'r');
        title([jd,' ',hr,' Spectrum of Specific Humidity']); grid;
        ppath = fullfile(path_raw_images,['qSpectrum_',jd,'_',hr,'_',graphformat]);
        print(graphdevice,ppath);

        close all;
    end
end
