function [MW_B_date_time,MW_E_date_time,MW_date_time, MW_code] = MW(sel_filename,station_id)

if ~exist('sel_filename') || isempty(sel_filename)
    error('No filenames specified.');
end
if ~exist('station_id') || isempty(station_id)
    error('No station code specified.');
end

report = [];
report_time = [];
MW_date_time = [];
MW_code = [];
fake_report = [];
fake_report_time = [];

for i=1:length(sel_filename)
    
    full_filename = sel_filename{i}; % current filename
    
    fid = fopen( full_filename );
    if fid==-1
        err = errordlg({'The following file could not be opened for reading:',full_filename},'Error reading file');
        uiwait(err); return;
    end
    h = waitbar(0,['File name: ' strrep(strrep(full_filename,'\','\\'),'_','\_')],...
        'Name','Extracting manual thunderstorm observations...', 'WindowStyle','modal');
    drawnow;
    file_contents = textscan( fid, '%s', 'whitespace', '', 'delimiter', '\n');
    file_contents = file_contents{1}; 
    
    report_i = zeros(length(file_contents),1);
    report_time_i = zeros(length(file_contents),1);
    MW_date_time_i = zeros(length(file_contents),1);
    MW_code_i = zeros(length(file_contents),1);

    MW_count=0;
    count = 0;
    for i=1:length(file_contents)
        textline = file_contents{i}; 
        if length(textline)<22
            continue
        end
        ac = str2num(textline(5:10));
        if ac ~= station_id
            continue
        end
        
        datetimestr=textline(11:22);
        yr=str2num(datetimestr(1:4));
        mo=str2num(datetimestr(5:6));
        dy=str2num(datetimestr(7:8));
        hr=str2num(datetimestr(9:10));
        mn=str2num(datetimestr(11:12));
        
        count = count+1;
        report_time_i(count) = datenum(yr,mo,dy,hr,mn,0);
        
        n=strfind(textline, 'MW');
                
        for ii=1:length(n)             
            tstm=str2num(textline(n(ii)+3:n(ii)+4));
            qual_check = str2num(textline(n(ii)+5));
            if  ismember(tstm,[17 29 95 96 97 98 99])& ismember(qual_check,[0 1])  % nos. 17, 29, 95-99, designated as thunderstorms by ASOS
                MW_count=MW_count+1;
                MW_date_time_i(MW_count)=datenum(yr,mo,dy,hr,mn,0);
                MW_code_i(MW_count) = tstm;
                report_i(count) = 1;
                if ishandle(h), waitbar(i/length(file_contents),h); drawnow; end

            end
        end
    end
    
    % eliminate extra zeros at the end:
    MW_date_time_i = MW_date_time_i(1:MW_count);
    MW_code_i = MW_code_i(1:MW_count);
    report_i= report_i(1:count);
    report_time_i = report_time_i(1:count);
    
    % add to results from previous files:
    MW_date_time = [MW_date_time; MW_date_time_i];
    MW_code = [MW_code; MW_code_i];
    report = [report;report_i];
    report_time = [report_time;report_time_i];

    if ishandle(h), waitbar(1,h); drawnow; close(h); end;
end

% sort MW codes by increasing date/time and eliminate repetitions (same date/time and same code)
MW_rows = unique([MW_date_time MW_code],'rows');
MW_date_time = MW_rows(:,1);
MW_code = MW_rows(:,2);
% sort reports by increasing date/time and eliminate repetitions (same date/time and same code)
report_rows = unique([report_time report],'rows');
report_time = report_rows(:,1);
report = report_rows(:,2);
% if "report" contains a 0 and a 1 for the same "report_time", then eliminate the 0:
dup_ind = find(diff(report_time)==0); % find repeated times
report(union(dup_ind, dup_ind+1)) = 1; % set both values of report to 1 (one should already by 1)
report_rows = unique([report_time report],'rows'); % eliminate repeated rows again
report_time = report_rows(:,1);
report = report_rows(:,2);

% if a 0 follows a first 1 by less than 15 minutes, then delete the 0 (use the next 0 as the estimated end):
% repeat until no 0 follows a first 1 by less than 15 minutes
for i=1:10
    find((report(1:end-2)==0 & report(2:end-1)==1 & report(3:end)==0) & (report_time(3:end)-report_time(2:end-1))<0.25/24)+2;
    ind_delete = find((diff(report_time)<0.25/24 & report(1:end-1)==1) & report(2:end)==0)+1;
    if isempty(ind_delete), break, end
    ind_keep = setdiff(1:length(report),ind_delete);
    report_time = report_time(ind_keep);
    report = report(ind_keep);
end

% If a reporting gap of > 2h occurs, add a fake non-thunderstorm report 1 h after the last report 
report_gap_ind = find(diff(report_time)>2/24); % find indices of gaps > 2 h between reports
if report(end) % if the last report indicates a thunderstorm, include this as a gap:
    report_gap_ind = union(report_gap_ind,length(report));
end
fake_report = zeros(length(report_gap_ind),1); % add a fake non-thunderstorm report for each gap
for i = 1:length(report_gap_ind)
    % define the time of the fake report as 1 h after the last report:
    fake_report_time(i,1) = report_time(report_gap_ind(i))+1/24;
end

% add fake reports and times to list of real reports and times:
report_time = [report_time;fake_report_time];
report= [report;fake_report];
% sort reports by increasing date/time:
[report_time, report_ind] = sort(report_time);
report = report(report_ind);


ts_report_ind = find(report); % indices of all thunderstorm reports
next_ts_report_ind = find(report)+1; % indices of all reports that follow a thunderstorm report
% Estimated beginnings (date/time of the first in a series of consecutive thunderstorm reports):
MW_B_date_time = report_time(setdiff(ts_report_ind, next_ts_report_ind));
% Estimated ends (date/time of the first non-thunderstorm report after a series of consecutive thunderstorm reports)
MW_E_date_time = report_time(setdiff(next_ts_report_ind, ts_report_ind));