[time,M]=read_RWS(filein,station_id);

function [time,M]=read_RWS(filein,station_id);

% % This function gives the datenum 'time' and the matrix 'M' of the 
% % radiosoundings read in the input file, obtained as ascii text by
% % 'cut and paste' from 'http://weather.uwyo.edu/upperair/sounding.html'.
% % The file can contain an arbitrary number of radiosoundings.
% % The 'station_id' is the number of the station, required to locate the
% % beginning of each data block. It must be in string format.
% % "M" has three dimensions: 
% % 1. Vertical measurements of every radiosounding (50 rows)
% % 2. Columns of the different quantities (11 columns):
% %       PRES   HGHT   TEMP   DWPT   RELH   MIXR   DRCT   SKNT   THTA   THTE   THTV
% %        hPa     m      C      C      %    g/kg    deg   knot     K      K      K 
% 
% % 3. Number of radiosoundings in the input file (variable).
% % "time" is the vector of datenum of each radiosounding, as read from the file. 
% % This variable has the dimension of the third dimension of M.
% 
% % Created by Federico Angelini, ISAC-CNR 
% % 30 Oct 2009


%first step to find the number of empty rows at the beginning of data
%block (considered header lines):
fid=fopen(filein);
nheader=0;
endlabel=[];
while 1
    tline = fgetl(fid);
    if length(tline)>4
        if tline(1:5)==station_id
                while 1
                    tline = fgetl(fid);
                    endlabel=strfind(tline,'.');
                    if length(endlabel)>3,break, end
                    nheader=nheader+1;
                end
        end
    end
    if length(endlabel)>3,break, end
end
fclose(fid);
disp([' Found ' num2str(nheader) ' header lines...'])

fid=fopen(filein);
M=nan*zeros(50,11,1);
j=1;
while 1
    tline = fgetl(fid);
    if length(tline)>4
        if tline(1:5)==station_id
            disp('Found a block of data ...')
            %Read the date:
            time(j)=datenum([tline(end-10:end-9) '-' tline(end-7:end-5) '-' tline(end-3:end) ], 'dd-mmm-yyyy')+str2num(tline(end-14:end-13))/24; 
            disp(datestr(time(j)))
            %Read 5 header lines:
            tline = fgetl(fid);
            tline = fgetl(fid);
            tline = fgetl(fid);
            tline = fgetl(fid);
            tline = fgetl(fid);
            %tline = fgetl(fid);
            %Read the data block:
            d=1;
            while 1
                tline = fgetl(fid);
                endlabel=strfind(tline,'        ');%Stop at the first lack of data
                if (~isempty(endlabel) | d==50 | isempty(tline)) & d>1 %Stop at the first lack of data, or anyway after 50 measurements...
                        disp('End block of data')
                        break
                end
                if isempty(endlabel)%put the data only if the line is complete (may happen only at the first line).
                    M(d,:,j)=str2num(tline);%M(50 levels, 11 variables, number of radiosounding)
                end
                d=d+1;
            end%while single data block
            M(d:50,:,j)=nan*zeros(50-d+1,11,1);
            j=j+1;
        end%if found data block
    end%if tline non empty
    if ~ischar(tline),   break,   end
end
fclose(fid);

%put nans if the first lines of data blocks are incomplete...
indnan=find(isnan(M(1:nheader-5,1,1)));
M(indnan,:,:)=nan;