%
%  This script calculates and graphs pdfs of accumulated tropical cyclone
%  rainfall at specified points of interest.
%       December 2011
%       Copyright 2011 WindRiskTech L.L.C.
%       Revised October 2018...did away with ksdensity.
%       Revised December 2018 using new vertical velocity algorithm that
%       accounts for surface roughness variations from ECMWF roughness
%       data.
%       Revised March, 2020 to include different units choices
%       Revised September, 2020 to use filtyears
%       Revised December 2020 to remove cos(latitude) factor from hy and to
%       allow for higher resolution topography
%-------------------------------------------------------------------------
ratemin=0;     % Minimum rainfal rate (mm/hr) to consider
rmin=20;       % Minimum storm rainfall (mm) to consider
smoothing='y'; % Simple smoothing of frequencies
%------------------------------------------------------------------------
params    %  Load parameters
pfield='rain';
%
clear poi plat plong
clf('reset')
%
if strcmp(randfac,'y')
    %disp('Caution: Random draw for wind radii turned on; results will differ with each call')
    if exist('tmr7','var')
        stop(tmr7)
        delete(tmr7)
        clear tmr7
    end   
    tmr7=timer;
    set(tmr7,'startdelay',6);
    set(tmr7,'timerfcn','try; delete(hms); catch err; end')
    hms=msgbox('Caution: Random draw for wind radii turned on; results will differ with each call');
    start(tmr7)
end 
%
if exist('shearstore','var') == 0
    shape='circ';
    load('temp.mat')
end
fs3='sorted';
load(fs3)
load temp filtyears
%
if strcmpi(topores,'high')
    load bathymetry_high
else
    load bathymetry
end    
[ntopo,~]=size(bathy);
topores=360/ntopo;
toporesi=1/topores;
%
clear rain rtemp rainrate rainratec
sfac=1./(topores.*60.*1852);
pifac=acos(-1)./180;
%
if exist('plat','var') == 0
    if strcmp(shape,'poly') == 1 || city_radius >= 2000
        prompt = {'Enter POI latitude','Enter POI longitude','Enter POI name'};
        dlg_title = 'Enter Point of Interest Data';
        num_lines = 1;
        answer = inputdlg(prompt,dlg_title,num_lines);
        if isempty(answer) || isempty(answer{1})
            return
        end            
        plat=str2double(answer{1});
        plong=str2double(answer{2});
        tmp=answer{3};
        poi=cellstr(tmp);
    else
        prompt = {'Use default POI (y or n; default is y)?'};
        dlg_title = 'Use default POI?';
        num_lines = 1;
        answer = inputdlg(prompt,dlg_title,num_lines);
        pswitch=answer{1};
        %pswitch=input('Use default POI (y or n; default is y)?','s');
        if isempty(pswitch)
            pswitch='y';
        end    
        %pswitch='y';
        if pswitch == 'n'
            %plat=input('Enter POI latitude');
            %plong=input('Enter POI longitude');
            %tmp=input('Enter POI name','s');
            prompt = {'Enter POI latitude','Enter POI longitude','Enter POI name'};
            dlg_title = 'Enter Point of Interest Data';
            num_lines = 1;
            answer = inputdlg(prompt,dlg_title,num_lines);
            plat=str2double(answer{1});
            plong=str2double(answer{2});
            tmp=answer{3};            
            poi=cellstr(tmp);
        else    
            plat=clat;
            plong=clong;
            clear poi
            datea(1)=upper(datea(1));
            poi=cellstr(datea);
        end
    end    
end    
%
knotfac=1852./3600;
%
[n,m]=size(ut);
if plong < 0
    plong=plong+360;
end  
%
ush=zeros(n,m);
vsh=ush;
vdrift=1.5.*3600./1852;
vdrift=vdrift.*latstore(1,1)./(abs(latstore(1,1))+1e-8);
if exist('u850store','var') == 1
    ush=5.*knotfac.*(ut-u850store);
    vsh=5.*knotfac.*(vt-vdrift.*cos(pifac.*latstore)-v850store);
end    
lat=latstore;
long=longstore;
v=vstore;
vse=vsestore;
%  
% Scale and randomize radii of maximum wind
%
[nrm,mrm]=size(rmstore);
rfac=magfac.*(1+zeros(nrm,mrm));
if strcmp(randfac,'y')
    rfac=magfac.*outer(nrm,mrm);
end  
%
rm=rmstore.*rfac;
rmse=rmsestore.*rfac;
%
if max(yearstore) == min(yearstore)
    freqyear(1)=freq;
end    
%
m_to_mm=1000;
rowa_over_rowl=0.00117;
%
bathy=max(bathy,-1);
dhdx=sfac*(circshift(bathy,-1,1)-circshift(bathy,0,1));
dhdy=sfac*(circshift(bathy,-1,2)-circshift(bathy,0,2));
%
ib=1+floor(toporesi.*plong);
ibp=ib+1;
if ibp > ntopo
    ibp=ibp-ntopo;
end    
ibs=1+floor(toporesi.*plong-0.5);
ibsp=ibs+1;
plongs=plong;
if ibs < 0
    ibs=ntopo;
    plongs=plong+360;
end    
if ibsp > ntopo
    ibsp=ibsp-ntopo;
end
jb=1+floor(toporesi.*(plat+90));
jbs=1+floor(toporesi.*(plat+90)-0.5);
b1=bathy(ib,jb);
b2=bathy(ib,jb+1);
b3=bathy(ibp,jb);
b4=bathy(ibp,jb+1);
dely=1+toporesi.*(plat+90)-jb;
delx=1+toporesi.*plong-ib;
d1=(1.-delx).*(1.-dely);
d2=dely.*(1.-delx);
d3=delx.*(1.-dely);
d4=delx.*dely;
h=exp(d1.*log(b1+11)+d2.*log(b2+11)+d3.*log(b3+11)+d4.*log(b4+11))-11;
%
b1=dhdx(ibs,jbs);
b2=dhdx(ibs,jbs+1);
b3=dhdx(ibsp,jbs);
b4=dhdx(ibsp,jbs+1);
dely=0.5+toporesi.*(plat+90)-jbs;
delx=0.5+toporesi.*plongs-ibs;
d1=(1.-delx).*(1.-dely);
d2=dely.*(1.-delx);
d3=delx.*(1.-dely);
d4=delx.*dely;
hx=(b1.*d1+b2.*d2+b3.*d3+b4.*d4)./cos(pifac.*plat);
%
b1=dhdy(ibs,jbs);
b2=dhdy(ibs,jbs+1);
b3=dhdy(ibsp,jbs);
b4=dhdy(ibsp,jbs+1);
hy=(b1.*d1+b2.*d2+b3.*d3+b4.*d4);
%
clear bathy
if exist('T600store','var')
    clear q900 q600
    [ q900, q600 ] = qs900b( T600store, vstore );
    wq = pointwshortnq(lat,long,q900,v,rm,vse,rmse,ut,vt,ush,vsh,plat,plong,h,hx,hy,timeres,wrad);
else
    w = pointwshortn(lat,long,v,rm,vse,rmse,ut,vt,ush,vsh,plat,plong,h,hx,hy,timeres);
    wq=max(w-wrad,0).*q900;
    clear w
end   
%
rainrate=eprecip.*m_to_mm.*3600.*rowa_over_rowl.*wq;
rain=timeres.*sum(rainrate,2); 
year=filtyears;  %  Modified September, 2020
qy=max(size(year));
%
% Create columns sorted by accumulated rain
%
ns=1:n;
dum=[ns',rain];
rainsort=sortrows(dum,2,'descend');
clear dum ns
%
nplot=2;
while nplot > 1
    nplot=menu('Rain Probability Plots','Stop','Peak Rainfall Rate', ...
        'Accumulated Rainfall');
    if nplot == 1    
        nplot=0;
        break
    elseif nplot == 2    
        if strcmp(Docdisp,'y')
            try
                delete(hButton)
            catch    
            end
        end 
        rainratemax=max(rainrate,[],2);
        clear u xi xint
        rmax=10*ceil(0.1*max(rainratemax));
        rmax=max(rmax,1);
        ratrate=sum(histcounts(rainratemax,[ratemin, 10000]))./n;
        xint=0.02*rmax;
        xi=ratemin:xint:rmax;
        xil=xi;
        xil(end+1)=100000; %#ok<SAGROW>
        xil(1)=max(xil(1),0.01);
        qx=max(size(xi));
        u=zeros(1,qx);
        for i=1:qy    
            mask=ismember(yearstore,year(i));
            rainnew=rainratemax.*mask;
            fm=histcounts(rainnew,xil);
            fm=fm*freqyear(i)*ratrate./sum(fm);
            qx=max(size(xi));
            f2=zeros(1,qx);
            for kk=1:qx
                f2(1,kk)=sum(fm(kk:qx));
            end    
            u=u+f2/qy;    
        end 
        %
        if strcmpi(smoothing,'y')
            us=u;
            us(3:end-2)=smooth134(u);
            us(3:end-2)=smooth134(us);
            us(3:end-2)=smooth134(us);
            fu=us;
        end
        %
        rp=1./(1.-exp(-u));
        y=real(rp);
        x=runitsfac*xi;
        h=semilogy(x,y,'.b');
        %set(gca,'ylim',[1 1e5])
        set(h,'markersize',12)
        set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        ylabel('Return Period (Years)','fontweight','bold')
        xlabl=char(strcat('Peak storm rainfall rate  (',runitslabel,'/hr)'));
        rswitch=1;
        xlabel(xlabl,'fontweight','bold')
        if strcmp(rgridx,'y')
             set(gca,'xgrid','on')
        end
        if strcmp(rgridy,'y')
             set(gca,'ygrid','on')
        end  
        title(poi,'fontweight','bold','interpreter','none');
        xadj=x;
        yadj=y;
        %
        % Add button to access documentation of script
        %
        if strcmp(Docdisp,'y')
            tempname=mfilename;
            mname = strcat('graphics_docs/html/',tempname,'_doc.html');
            comstr=['open ',mname,';delete(hButton)'];
            labelStr = '<html><center><a href="">Doc';
            hButton = uicontrol('string',labelStr,'pos',[10,10,50,20],'callback',comstr); 
        end       
        %
    elseif nplot == 3    
        %
        if strcmp(Docdisp,'y')
            try
                delete(hButton)
            catch
            end
        end     
        rswitch=2;
        %rain=rain/25.4;rmin=rmin/25.4;  % Conversion to inches
        xmax=10*ceil(0.1*max(rain));
        xmax=max(xmax,10);
        rat=sum(histcounts(rain,[rmin, 10000]))./n;
        xint=0.02*xmax;
        xi=(rmin:xint:xmax);        % Changed 12/11/13 and again 8/27/2018 and 10/13/2018
        xil=xi;
        xil(end+1)=100000;  %#ok<SAGROW> % This just allows a category in histcounts including all events with rain>rainmax
        xil(1)=max(xil(1),0.01); % This is essential to exclude masked events
        qx=max(size(xi));
        %
        f=zeros(1,qx);
        for i=1:qy
            % 
            % Define rainnew as the array containing the point peak rain totals in each
            %    year i. This was added in October, 2016 and revised in
            %    October 2018.
            %
            mask=ismember(yearstore,year(i));
            rainnew=mask.*rain';
            u=histcounts(rainnew,xil);
            fd=freqyear(i)*rat*u./(0.001+sum(u));
            fdc=zeros(1,qx);
            for k=1:qx
                fdc(k)=sum(fd(k:end));
            end
            f=f+fdc/qy;
        end   
        %
        if strcmpi(smoothing,'y')
            fs=f;
            fs(3:end-2)=smooth134(f);
            fs(3:end-2)=smooth134(fs);
            fs(3:end-2)=smooth134(fs);
            f=fs;
        end
        %
        rp=1./(1-exp(-f));    % Changed 12/10/2013
        x2=runitsfac*xi;
        y2=rp;
        h=semilogy(x2,y2,'.b');
        set(h,'markersize',18)
        hold on
        semilogy(x2,y2,'b','linewidth',1)
        hold off
        set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
        ylabel('Return Period (Years)','fontweight','bold')
        xlabl=char(strcat('Storm total rainfall (',runitslabel,')'));
        xlabel(xlabl,'fontweight','bold')
        if strcmp(rgridx,'y')
             set(gca,'xgrid','on')
        end
        if strcmp(rgridy,'y')
             set(gca,'ygrid','on')
        end  
        title(poi,'fontweight','bold','interpreter','none');
        xadj=x2;
        yadj=y2;
        %
        % Add button to access documentation of script
        %
        if strcmp(Docdisp,'y')
            tempname=mfilename;
            mname = strcat('graphics_docs/html/',tempname,'_doc.html');
            comstr=['open ',mname,';delete(hButton)'];
            labelStr = '<html><center><a href="">Doc';
            hButton = uicontrol('string',labelStr,'pos',[10,10,50,20],'callback',comstr); 
        end               
    end
end       
x=x2;y=y2;
if exist('tmr7','var')
    stop(tmr7)
    delete(tmr7)
    clear tmr7
end
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 
rainsave=rain;
save raint rainsave