function [ ] = prepf(filename)
%
%  This script reads in the raw tracks files, runs "utransfull" to determine
%  translation speed, adds increments to the circular wind speed based,
%  orders the tracks by maximum lifetime wind and also
%  by maximum wind at a point of interest, and stores the results. 
%  Note that this function clears main event set matrices from base
%  workspace.
%
%  Copyright 2012 WindRiskTech L.L.C.
%  Last modified October, 2016
%  Modified January, 2019 to store index of location of maximum winds
%  Modified May, 2019 (code cleaned up)
%  Modified January, 2020 to include outer radius calculation
%  Modified September 2020 to write out filtyears
%  Modified December, 2020 to accept netcdf files as input
%--------------------------------------------------------------------------
%
params    %  Load parameters
%
filename=char(filename);
if isempty(filename) 
    return
end    
datea=filename;
storm=[]; 
%
% Copy event set files to temporary files
%
if exist('poly,in','file')
    delete poly.in
end 
%
fs=[edir,datea,'.zip'];
%
try
    if exist('7z.exe','file')
        comm=['7z e -aoa -bso0 ',fs];
        system(comm);
    else
        unzip(fs)
    end    
catch
    disp('Unzip error');
    return
end   
%
%  Determine date file was created for use in correcting surface pressure
%
rawfile=dir(fs);
filedate=rawfile.datenum;
pfiledate=736116;     %   Datenumber corresponding to 1 June 2015, when pressure algorithm was changed in fortran codes
%
if strcmpi(datea(end-2:end),'.nc')
    file=[datea storm];
    convertnc(file)
else    
    file=[datea storm '.mat'];
    movefile(file,'temp.mat','f') 
end   
%
% Load event sets
%
shape='circ';  % (default value used for old event sets)
load('temp.mat') %#ok<LOAD>
%
[nn,m]=size(vstore);
%
% This section of code added on September 27 2020 to deal with tracks over
% Himalayas. This section should not be needed for tracks prior to
% implementation of CD or for after Sept. 27 2020 when FORTRAN code
% corrected
%
[~,jmaxv]=min(vstore,[],2);          % Find when circular wind speeds drop to zero
[~,jmaxlat]=min(abs(latstore),[],2); %#ok<NODEF>
% 
if sum(jmaxv) < sum(jmaxlat)
    for i=1:nn
        latstore(i,jmaxv(i):end)=0;  
        daystore(i,jmaxv(i):end)=0;  
        hourstore(i,jmaxv(i):end)=0;  
        longstore(i,jmaxv(i):end)=0;  
        monthstore(i,jmaxv(i):end)=0;  
        pstore(i,jmaxv(i):end)=0;  
        rhstore(i,jmaxv(i):end)=0;  
        rmsestore(i,jmaxv(i):end)=0;  
        rmstore(i,jmaxv(i):end)=0;  
        shearstore(i,jmaxv(i):end)=0;  
        T600store(i,jmaxv(i):end)=0;  %#ok<*AGROW>
        u850store(i,jmaxv(i):end)=0; 
        v850store(i,jmaxv(i):end)=0; 
        vsestore(i,jmaxv(i):end)=0; 
        vpstore(i,jmaxv(i):end)=0; 
    end    
    save temp latstore daystore hourstore longstore monthstore pstore rhstore -append
    save temp rmsestore rmstore shearstore T600store u850store v850store vsestore vpstore -append
end
%
if exist('freqyear','var') == 0
    freqyear=freq;   %#ok<USENS> %  For single period event sets, set freqyear equal to freq
    save temp freqyear -append
end    
%
if exist('yearstore','var') == 0
    %
    yearstore=zeros(1,nn)+1;   %  For single period event sets, set yearstore = 1
    save temp yearstore -append
end    
%
% Define filtyears 
%
filtyears=min(yearstore):max(yearstore); % Use prepfilterf if filtering of years desired   
%
if exist('T600store','var')
    [router] = roquick(T600store,rmstore,vstore);
    save temp router -append
end    
%
%  Correct central surface pressures for pre June 1 2015 event sets
%
if filedate < pfiledate
    vms=vstore.*1852/3600;  %  Convert wind speed to m/s
    pstore=7+pstore.*exp(0.8.*vms.^2./(287*300))-0.6.*vms;   
    clear vms
    save temp pstore -append
end
%
if exist('vsestore','var') == 0 || strcmp(seyewall,'n')
    vsestore=zeros(nn,m);
    rmsestore=zeros(nn,m);
    save temp vsestore rmsestore -append
end    
%
%  Super-gradient wind pre-processor
%
%[~,jmax]=min(vstore,[],2);
%jmax=jmax-1;
%vtemp=vstore;
%vupper=105;
%vlower=65;
%mag=0.25;
%for i=1:nn,
%    for j=1:jmax(i);
%        if vstore(i,j) >= vlower && vstore(i,j) < vupper
%            coef=4.*(vstore(i,j)-vlower).*(vupper-vstore(i,j))./(vupper-vlower)^2;
%            vtemp(i,j)=(1+mag*coef).*vstore(i,j);
%        end
%    end
%end
%vstore=vtemp;
%clear vtemp
%
%  Calculate translation velocity and wind speed increments
%
if exist('u850store','var')
    [ut,vt,uinc,vinc ] = utransfull(latstore,longstore,vstore,u850store,v850store);
    [utse,vtse,uincse,vincse ] = utransfull(latstore,longstore,vsestore,u850store,v850store);
else
    [ut,vt,uinc,vinc ] = utransfull(latstore,longstore);
    utse=ut;
    vtse=vt;
    uincse=uinc;
    vincse=vinc;
end
%
vnet=vstore+sqrt(uinc.^2+vinc.^2); % Maximum ground-relative wind speed
[vnetmax,ivmax]=max(vnet,[],2);    % Storm lifetime maximum wind
vnetmax=transpose(vnetmax);
%
% Create mask for secondary eyewall events and find maximum ground relative
% wind for them
%
semask=min(vsestore,1);
vnetse=semask.*(vsestore+sqrt(uincse.^2+vincse.^2));
[vnetmaxse,~]=max(vnetse,[],2);    % Storm lifetime maximum wind
vnetmaxse=transpose(vnetmaxse);
%
passi=zeros(nn,m)+1;
if strcmp(shape,'poly')
    [polytype, passi, passe, xint, yint, jint, kint, kfrac] = polyfilter(latstore,longstore);
    if strcmp(polytype,'closed')
        poly=load('poly.in'); 
        x=poly(:,1); 
        y=poly(:,2);
        A=polyarea(x',y');
        if A < poly_area
            polytype='open';
        end
    end        
    %
    % Find maximum ground-relative wind at time track crosses line segements
    %
    [~,xsize]=size(xint);
    vseg=zeros(nn,xsize); % Added 11/28/17
    for nx=1:xsize % Find intersection with maximum wind
        for i=1:nn
            vseg(i,nx)=passe(i).*vnet(i,max(kint(i,nx),1)); % Changed 11/28/17            
        end
    end    
    [vmax,iimax]=max(vseg,[],2);
    kmax2=zeros(1,nn);
    for i=1:nn
        kmax2(i)=kint(i,iimax(i));
    end      
    vmax=vmax';
    if strcmp(polytype,'closed')
        vmaxcl=max(vnet.*passi,[],2);  % Maximum ground-relative wind speed within closed polygon
        vmaxcl=transpose(vmaxcl);
        vmax=max(vmaxcl,vmax);         % For narrow closed polygons, max wind may be on boundary
        clear vmaxcl
    end      
else  % Case of circular filter, find maximum wind at point of interest
    if city_radius < 2000 % Exclude cases where city radius encompasses whole basin
        [~, vmax, kmax2] = best_filter(vnet,latstore,longstore,city_radius,clat,clong);
        yearbig=repmat(yearstore',[1,m]);
        datestore=datenum(double(yearbig),monthstore,daystore,hourstore,0,0); %#ok<DATNM>
        [vs, dir1, dayk] = pointseriesn(latstore,longstore,vstore,rmstore.*magfac, ...
                 vsestore,rmsestore.*magfac,datestore,uinc,vinc, ...
                 clat,clong,timeres);     
        vsmax=max(vs,[],2);
        vsmax=transpose(vsmax);
    else
        vmax=vnetmax;
        kmax2=ivmax;
    end    
end
[mf,nf]=size(vmax);
if mf > nf
    vmax=transpose(vmax);
end    
%
%  Sort by local and by liftime maximum wind speeds
%
[vsort,nsort]=sortrows(transpose(vnetmax),-1);
vsort=transpose(vsort); 
nsort=transpose(nsort); 
%
[vsortp,nsortp]=sortrows(transpose(vmax),-1);
vsortp=transpose(vsortp); 
nsortp=transpose(nsortp); 
%
% Save to files
%
file='sorted';
%
save(file,'vsort','nsort','vsortp','nsortp','freq',...
    'datea','storm')
if strcmp(shape,'circ') && city_radius < 2000
    [vsortc,nsortc]=sortrows(transpose(vsmax),-1);
    vsortc=transpose(vsortc);
    nsortc=transpose(nsortc);
    save sorted vsortc nsortc -append
end    
eset=datea;
if strcmpi(datea(end-2:end),'.nc')
    eset=datea(1:end-3);
end
%
save temp uinc vinc ut vt vnet vmax vnetmax passi kmax2 eset filtyears -append
if strcmp(shape,'poly')
    save temp polytype passe xint yint jint kint kfrac -append
elseif city_radius < 2000
    save temp vs dir1 dayk vsmax -append
end    
%
save temp uincse vincse utse vtse vnetse vnetmaxse -append
%
% Create list of protected variables that will not be cleared by individual
% scripts
%
protected_variables=who('-file', 'temp.mat'); % All variables in 'temp'
%
% Additional protected variables:
%
additional_variables={'protected_variables','x','x2','x3','x4','y','y2','y3', ...
    'y4','z','z2','gplot','tmr','maptoolpresent','xmin*','xmax','ymin','ymax', ...
    'imin','imax','kmin','kmax2', 'histsize', 'sumtotal', 'datea', ...
    'fr','nx','ny','nt','rplot','ax','ay','aspect','nplot','storm'};
%
protected_variables=[protected_variables;additional_variables'];
save temp protected_variables -append
%
evalin( 'base', 'clear *store' )  %  This command clears main arrays from base workspace
disp('prep finished')
end
