function [rainrate,rain,dayk,rx,f,rp] = raincalc( POI,rmin,rmax,vcirc, ...
    rmw,lat,long,yearst,month,day,hour,T600,u850,v850)
%
%  This script calculates tropical cyclone rain rates and accumulated
%  rainfall at specified points of interest.
%--------------------------------------------------------------------------
%  INPUTS: 
%
%  POI(:,1) is a vector of latitudes, and POI(:,2) is a vector of
%  longitudes of points of interest at which times series of wind speed and
%  direction are calculated. 
%
%  rmin and rmax are the minimum and maximun storm total rainfall (mm) to
%  be considered in constructing the storm total rainfall return periods.
%  (Note that rmin itself is forced to be larger than a critical value set
%  below.)
%
%  All of the following, except year, are n x m arrays where n is the storm
%  number and m is the two-hour point along the track. Year is a a vector
%  dimensioned (1,n).
%
%  vcirc: Maximum azimuthally averaged 10 meter wind speed in m/s
%  rmw:   Radius of maximum winds in km
%  lat:   Latitude in degrees (negative for southern hemisphere)
%  long:  Longitude in degrees east of Greenwich meridian
%  yearst:Year of storm (this is a vector of dimension [1,n])
%  month: Calendar month (1-12)
%  day:   Day of the month 
%  hour:  Hour of the day (Greenwich mean time)
%  T600:  600 hPa environmental temperature (K) interpolated to storm
%           center
%  u850:  850 hPa environmental zonal wind (m/s) interpolated to storm
%           center
%  v850:  850 hPa environmental meridional wind (m/s) interpolated to storm
%           center
%--------------------------------------------------------------------------
%  OUTPUTS:
%
%  rainrate is a matrix containing the rain rates in mm/hr at the POIs. It
%  is dimensioned n x k x s, where n is the number of storms, k is the
%  length of the time series, which is set internally, and s is the number
%  of points of interest. (Note that if s=1, n x k matrices are returned.)
%  The matlab datenumbers corresponding to these series are contained in
%  dayk.
%
%  rain is a vector of dimension n x s containing the storm total rainfall
%  in mm at the POI.
%
%  rx is the range of storm total rainfall (mm) considered for histograms
%  of rainfall at the POI, determined by the paramters rmin and rmax
%  described above.
%
%  f is the annual exceedence frequency, and rp is the return period
%  (inverse of the annual probability) of storm total rainfall at the POI.
%  These are s x m matrices, where s is the number of points of interest and
%  m is the size of the range of rain rates in the histogram.
%
%  Note that this function can be slow if many points of interest are used.
%      
%  Copyright WindRiskTech LLC, 2016
%  Last modified for use in prediction models, February, 2020
%
%-------------------------------------------------------------------------
%  It is very important to put lower bounds on rain for data sets in which
%  the number of storms causing rain at the POI is a fraction of the total
%  number of storms. No harm is done by having such bounds when that
%  fraction is one or close to one. 
%
%-------------------------------------------------------------------------
rminlower=20;     % Minimum storm rainfall (mm) to consider
%------------------------------------------------------------------------
%
rmin=max(rmin,rminlower);
params    %  Load parameters
%
plat=POI(:,1)';
plong=POI(:,2)';
%
sx=max(size(plong));
for i=1:sx
    if plong(i) < 0
        plong(i)=plong(i)+360;
    end
end   
%
% Some dummy arrays for prediction model use
%
vsestore=0.*vcirc;
rmsestore=vsestore;
%
% Convert vcirc, u850 and v850 to knots
%
knotfac=1852./3600;
knotfaci=1/knotfac;
vcirc=vcirc.*knotfaci;
u850=u850.*knotfaci;
v850=v850.*knotfaci;
%
load bathymetry %#ok<LOAD>
clear rain rtemp rainrate rainratec
sfac=1./(0.25.*60.*1852);
pifac=acos(-1)./180;
%

%
[n,m]=size(vcirc);
%
% Calclate translation speeds (this code used in predictive mode)
%
[ut,vt,~,~] = utransfull(lat,long,vcirc,u850,v850);
ush=zeros(n,m);
vsh=ush;
vdrift=1.5.*knotfaci;
vdrift=vdrift.*lat(1,1)./(abs(lat(1,1))+1e-8);
if exist('u850','var') == 1
    ush=5.*knotfac.*(ut-u850);
    vsh=5.*knotfac.*(vt-vdrift.*cos(pifac.*lat)-v850);
end    
%  
% Scale and randomize radii of maximum wind
%
[nrm,mrm]=size(rmw);
rfac=magfac.*(1+zeros(nrm,mrm));
if strcmp(randfac,'y')
    rfac=magfac.*outer(nrm,mrm);
end  
%
rm=rmw.*rfac;
rmse=rmsestore.*rfac;
%
% Define parameters for an array vnew that contains all the peak winds in a
% given year
%
m_to_mm=1000;
rowa_over_rowl=0.00117;
%
bathy=max(bathy,-1); %#ok<NODEF>
dhdx=sfac*(circshift(bathy,-1,1)-circshift(bathy,0,1));
dhdy=sfac*(circshift(bathy,-1,2)-circshift(bathy,0,2));
%
h=zeros(1,sx);
hx=zeros(1,sx);
hy=zeros(1,sx);
%
for i=1:sx
    ib=1+floor(4.*plong(i));
    ibp=ib+1;
    if ibp > 1440
        ibp=ibp-1440;
    end    
    ibs=1+floor(4.*plong(i)-0.5);
    ibsp=ibs+1;
    plongs=plong(i);
    if ibs < 0
        ibs=1440;
        plongs=plong(i)+360;
    end    
    if ibsp > 1440
        ibsp=ibsp-1440;
    end
    jb=1+floor(4.*(plat(i)+90));
    jbs=1+floor(4.*(plat(i)+90)-0.5);
    b1=bathy(ib,jb);
    b2=bathy(ib,jb+1);
    b3=bathy(ibp,jb);
    b4=bathy(ibp,jb+1);
    dely=1+4.*(plat(i)+90)-jb;
    delx=1+4.*plong(i)-ib;
    d1=(1.-delx).*(1.-dely);
    d2=dely.*(1.-delx);
    d3=delx.*(1.-dely);
    d4=delx.*dely;
    h(i)=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+4.*(plat(i)+90)-jbs;
    delx=0.5+4.*plongs-ibs;
    d1=(1.-delx).*(1.-dely);
    d2=dely.*(1.-delx);
    d3=delx.*(1.-dely);
    d4=delx.*dely;
    hx(i)=(b1.*d1+b2.*d2+b3.*d3+b4.*d4)./cos(pifac.*plat(i));
    %
    b1=dhdy(ibs,jbs);
    b2=dhdy(ibs,jbs+1);
    b3=dhdy(ibsp,jbs);
    b4=dhdy(ibsp,jbs+1);
    hy(i)=(b1.*d1+b2.*d2+b3.*d3+b4.*d4)./cos(pifac.*plat(i));
end    
%
clear bathy
clear q900 q600
[ q900, q600 ] = qs900b( T600,vcirc ); %#ok<ASGLU>
[wq,dayk] = pointwnqp(lat,long,q900,vcirc,rm, ...
   vsestore,rmse,yearst,month,day,hour,ut,vt,ush,vsh,plat, ...
   plong,h,hx,hy,timeres,wrad);
%
wpfac=1;
%{
if wprofile == 2
    wpfac=2;
end
%}    
rainrate=wpfac.*eprecip.*m_to_mm.*3600.*rowa_over_rowl.*wq;
rain(:,:)=timeres.*sum(rainrate,2);
%
%xmax=10*ceil(0.1*max(rain));
xmax=rmax;
xmax=max(xmax,10);
%
%maxrain=10*ceil(0.1*max(rain));
%
xint=0.02*xmax;
xi=(rmin:xint:xmax);         % Changed 12/11/13 and 10/14/2018
qx=max(size(xi));
xil=xi;
xil(end+1)=100000;
xil(1)=max(xil(1),0.01); % This is essential to exclude masked events
%
f=zeros(sx,qx);
year=min(yearst):max(yearst);
qy=max(size(year));
for j=1:sx
    rat=sum(histcounts(rain(:,j),[rmin, 10000]))./n;
    for i=1:qy
        % 
        % Define rainnew as the array containing the point accumulated rainfall in each
        %    year i. This was added in October, 2016
        %
        mask=ismember(yearst,year(i));
        rainnew=mask'.*rain(:,j);
        u=histcounts(rainnew,xil);
        fd=rat*u./(0.001+sum(u));
        fdc=zeros(1,qx);
        for k=1:qx
            fdc(k)=sum(fd(k:end));
        end
        f(j,:)=f(j,:)+fdc/qy; 
    end    
end    
rp=1./(1-exp(-f));    % Changed 12/10/2013
rx=xi;
%