function [jd_ceil,cloudcode,cloud1,cloud2,cloud3,ceil_range,ceil_b] = load_CL31(ddd,dfl)

%{
    reads and decodes raw CL31 data files
    
    INPUTS:
       ddd = day of year
       dfl is the full-path string to the CL31 file
    OUTPUTS:
       jd_ceil is the decimal day of year
       cloudcode is detetection parameter, defined in code below
       cloud1, cloud2 and cloud3 are cloud base heights in m
       ceil_range is the range of each gate in meters
       ceil_b is the backscatter coefficient in the units (srad km)^(-1)
%}

nn = 7000;          % oversize default array index
nbits = 20;         % number of bits for reflectivity
start_height = 0;   % in m
gates_number = 770; % predefine CL31 gates and resolution
resolution = 10;

% Declare Arrays
tempcloud1 = NaN(nn,1);
tempcloud2 = NaN(nn,1);
tempcloud3 = NaN(nn,1);
tempcloudcode = NaN(nn,1);
tempceil_sec = NaN(nn,1);
tempceil_min = NaN(nn,1);
tempceil_hour =NaN(nn,1);
tempceil_b = NaN(gates_number,nn);

% open file
fid  =fopen(dfl,'r');
if fid == -1
    error(['Cannot read ' dfl]);
end

% dump header
dmp = fgetl(fid);   %dump line -Ceilometer Logfile
dmp1 = fgetl(fid);  %dump line -File created: 17.11.2008 0:00:00

%start read loop
k = 0; % count time samples
while ~feof(fid)
    k = k+1;
    time = fscanf(fid,'%*11c %2d:%2d:%2d',[3 1]); %dump date, read hh:mm:ss
    if ~isempty(time)
        tempceil_hour(k) = time(1);
        tempceil_min(k) = time(2);
        tempceil_sec(k) = time(3);

        fgetl(fid);       % go to end of line
        fgetl(fid);       % get rid of CL31

        tempht = fgetl(fid); % cloud base heights  20 00940 01235 ///// 000000008080
        if tempht(1)=='/' || tempht(2)=='A' % Alarms present!
            tempcloudcode(k) = -1; % flag error
        else % no alarms present
            %{
                tempcloudcode code:
                    0  No significant backscatter
                    1  1 cloudbase detected
                    2  2 cloudbase detected
                    3  3 cloudbase detected
                    4  full obscuration, but no cloud base detected
                    5  some obscuration, but determined to be transparent
            %}
            tempcloudcode(k) = sscanf(tempht(1:2),'%1d%*1c');
            if tempcloudcode(k) < 4 % codes 4:fully obscured, no base detected; 5:some obscuration, but determined transparent
                if tempcloudcode(k) > 0
                    tempcloud1(k) = sscanf(tempht(4:8),'%5d',1);
                    if tempcloudcode(k) > 1
                        tempcloud2(k) = sscanf(tempht(10:14),'%5d',1);
                        if tempcloudcode(k) == 3
                            tempcloud3(k) = sscanf(tempht(16:20),'%5d',1);
                        end
                    end
                end
            end
        end % no alarms present

        % position 21 is a space before the detection status bits in 12
        if k<2 % only check detection status once per file
            Detection_status = tempht(22:22+12-1); % hexadecimal digits
            units_bit = sscanf(Detection_status(11),'%1x',1)>=8;
            % units_bit true if meters
            % see table of hexadecimal digits and detection bits in ceilo manual
            if ~units_bit
                mm = 0.3048; % conversion coefficient from feet to meters
            end
        end

        fgetl(fid);         % skip next line:  -1 ///  0 ///  0 ///  0 ///  0 ///
        line3 = fgetl(fid); % 100 N 100 +21  74  219 +13    0 LF7HN1  38
        if k<2 % update scale and resolution once per file
            scaleres = sscanf(line3,'%5d');
            scale = scaleres(1)/100;
%             resolution = scaleres(2);
%             gates_number = scaleres(3);
        end

        line = fgets(fid);
        tempceil_b(:,k) = sscanf(line,'%5x');

    end

    fgetl(fid); % dump 
    fgetl(fid);  % get rid of blank line
end
fclose(fid);
k = k-1; % length of timeseries
ceil_b = tempceil_b(:,1:k);

% wrap 2s complement for nbits = 20
ceil_b2 = ceil_b;
index = find(ceil_b>=2^(nbits-1));
if ~isempty(index); ceil_b2(index) = ceil_b(index)-2^nbits; end

if scale==1
    ceil_b = ceil_b2*1e-5; % scaling factor of 1e5/(km srad), p50 of CL31 manual
else
    ceil_b = ceil_b2*scale*1e-5;
end
% ceil_b is now VOLUME BACKSCATTER COEFFICIENT [units: 1/(km*steradians)].

% convert units to meters if necessary
if units_bit
    cloud1 = tempcloud1(1:k);      % meters
    cloud2 = tempcloud2(1:k);
    cloud3 = tempcloud3(1:k);
else
    cloud1 = tempcloud1(1:k)*mm;   % ft to m
    cloud2 = tempcloud2(1:k)*mm;
    cloud3 = tempcloud3(1:k)*mm;
end

cloudcode = tempcloudcode(1:k);
ceil_hr = tempceil_hour(1:k);
ceil_mn = tempceil_min(1:k);
ceil_sc = tempceil_sec(1:k);
jd_ceil(1:k,1) = ddd + ceil_hr/24 + ceil_mn/1440 + ceil_sc/86400;

ceil_range = (start_height+(0:gates_number-1)*resolution); % profile range is 0 to 7690 m

% do some tests to see if last points are not from the next hour...
tt = find(diff(ceil_hr)<0);
cloudcode(tt+1) = [];
cloud1(tt+1) = [];
cloud2(tt+1) = [];
cloud3(tt+1) = [];
ceil_b(:,tt+1) = [];
jd_ceil(tt+1) = [];

% does the same for extra NaNs...
tt = find(~isfinite(ceil_hr));
cloudcode(tt) = [];
cloud1(tt) = [];
cloud2(tt) = [];
cloud3(tt) = [];
ceil_b(:,tt) = [];
jd_ceil(tt) = [];

end
