function [w] = pointwfield(latstore,longstore,vstore,rmstore, ...
    vsestore,rmsestore,ut,vt,us,vs,plat,plong,h,hx,hy)
%
%  This function takes the nn X m 'store' matrices  and, 
%  using the latitude (plat) and longitude (plong) vectors of points of 
%  interest, calculates the spatial distribution of vertical velocity 
%  centered at plong, plat. The routine utrans.m should be run before this 
%  one. The arrays us and vs are the vertical shears (m/s) used to estimate 
%  the baroclinic components of the vertical motion. The matrices h, hx, 
%  and hy contain topographic heights and their derivatives in x and y. 
%      Copyright 2009 WindRiskTech L.L.C.
%      Revised August, 2017 to include shear dotted with storm entropy
%      in hxmod and hymod definitions
%      Revised December, 2018 to use ECMWF neutral drag coefficients
%      Revised August 2023 to correct Omega for use in dM/dr terms
%      Revised October 23 2024 to correct line 328 and 423 with factor of
%      0.5. To compensate for radar/gauge comparisons, defualt precip
%      efficiency increased to 1.0. 
%-------------------------------------------------------------------------
%
params    %  Load parameters
%
se=max(max(vsestore));  % Test for secondary eyewalls
deltari=1./deltar;
timereswi=1./(3600.*timeresw);
Hi=1./Htrop;
%
[nn,m]=size(ut);
[~,sx]=size(plong);
[~,sy]=size(plat);
ngrid=0;
if sx == sy 
    ngrid=1;
end    
if ngrid == 1
    sy=1;
else    
    sy=max(size(plat));
end
w=zeros(nn,3,sx,sy);
w2=zeros(nn,3,sx,sy);
%
if plong(1) < 0
    plong=plong + 360;
end    
%
pifac=acos(-1)/180;
dfac=60.*1.852;
sfac=1./(0.25.*60.*1852);
knotfac=1852./3600;
ut=ut.*knotfac;
vt=vt.*knotfac;
Omega=acos(-1)./(12.*3600);    % Earth angular velocity parameter
%
% Load neutral drag coefficients
%
load C_Drag500 cd
cd=0.9*cd./(1+50.*cd); % This corrects the drag coefficient to be better applied to gradient wind;
                       % see Esau et al. (2004)
cd=max(cd,1e-3); % Align over-water values with Fortran program (including some wave drag effect)
%
% Interpolate drag coefficient and its gradients to POI
%
sy=max(size(plat));
sx=max(size(plong));
dcddx=sfac*(circshift(cd,-1,1)-circshift(cd,0,1));
dcddy=sfac*(circshift(cd,-1,2)-circshift(cd,0,2));
%
cdrag=zeros(sx,sy);
cdx=zeros(sx,sy);
cdy=zeros(sx,sy);
for i=1:sx
    for j=1:sy
        ib=1+floor(4.*plong(i));
        ibp=ib+1;
        if ibp > 1440
            ibp=1;
        end   
        ibs=1+floor(4.*plong(i)-0.5);
        plongs=plong;
        if ibs < 1
            ibs=1440;
            plongs=plong+360;
        end    
        ibsp=ibs+1;
        jb=1+floor(4.*(plat(j)+90));
        jbs=1+floor(4.*(plat(j)+90)-0.5);
        b1=cd(ib,jb);
        b2=cd(ib,jb+1);
        b3=cd(ibp,jb);
        b4=cd(ibp,jb+1);
        b1x=dcddx(ibs,jbs);
        b2x=dcddx(ibs,jbs+1);
        b3x=dcddx(ibsp,jbs);
        b4x=dcddx(ibsp,jbs+1);
        b1y=dcddy(ibs,jbs);
        b2y=dcddy(ibs,jbs+1);
        b3y=dcddy(ibsp,jbs);
        b4y=dcddy(ibsp,jbs+1);
        dely=1+4.*(plat(j)+90)-jb;
        delx=1+4.*plong(i)-ib;
        d1=(1.-delx).*(1.-dely);
        d2=dely.*(1.-delx);
        d3=delx.*(1.-dely);
        d4=delx.*dely;
        cdrag(i,j)=d1*b1+d2*b2+d3*b3+d4*b4;
        dely=0.5+4.*(plat(j)+90)-jbs;
        delx=0.5+4.*plongs(i)-ibs;
        d1=(1.-delx).*(1.-dely);
        d2=dely.*(1.-delx);
        d3=delx.*(1.-dely);
        d4=delx.*dely;
        cdx(i,j)=(d1*b1x+d2*b2x+d3*b3x+d4*b4x)./cos(pifac.*plat(j));
        cdy(i,j)=d1*b1y+d2*b2y+d3*b3y+d4*b4y;
    end
end
cdrag=max(cdrag,0.0015); % Added December, 2018 to put lower bound on drag coefficeint over water
%
% Three statements below added Dec. 2018 to reduce drag coefficient of
% near-coastal locations to reflect response time of landfalling TCs. 
%
%disp([num2str(cdrag(1,1)),' ',num2str(cdx(1,1)/sfac),' ',num2str(cdy(1,1)/sfac)]);
%cdrag=min(cdrag,1.5e-3.*(1+max(h,0)./100));
%cdx=cdx.*0.01.*min(max(h,0),100);
%cdy=cdy.*0.01.*min(max(h,0),100);
%
% Calculate distance of each POI from track
%
dx=zeros(nn,m,sx,sy);
dy=dx;
for i=1:sx
    for jj=1:sy
        if ngrid == 1
            j=i;
        else
            j=jj;
        end    
        dx(:,:,i,j)=dfac.*cos(pifac.*plat(j)).*(plong(i)-longstore);
        dy(:,:,i,j)=dfac.*(plat(j)-latstore);
    end
end    
rfull=sqrt(dx.*dx+dy.*dy);
rfull=max(rfull,0.5);
denom=1./rfull;
costhetafull=dx.*denom;
sinthetafull=dy.*denom;
%
vfull=zeros(nn,3,sx,sy);
rmfull=zeros(nn,3,sx,sy);
vsefull=zeros(nn,3,sx,sy);
rmsefull=zeros(nn,3,sx,sy);
latfull=zeros(nn,3,sx,sy);
longfull=zeros(nn,3,sx,sy);
utfull=zeros(nn,3,sx,sy);
vtfull=zeros(nn,3,sx,sy);
usfull=zeros(nn,3,sx,sy);
vsfull=zeros(nn,3,sx,sy);
hfull=zeros(nn,3,sx,sy);
hyfull=zeros(nn,3,sx,sy);
hxfull=zeros(nn,3,sx,sy);
cdfull=zeros(nn,3,sx,sy);
cdyfull=zeros(nn,3,sx,sy);
cdxfull=zeros(nn,3,sx,sy);
%
for i=1:sx
    for jj=1:sy
        if ngrid == 1
            j=i;
        else
            j=jj;
        end          
        for n=1:nn
            vfull(n,:,i,j)=vstore(n,:);
            rmfull(n,:,i,j)=rmstore(n,:);
            vsefull(n,:,i,j)=vsestore(n,:);
            rmsefull(n,:,i,j)=rmsestore(n,:);            
            latfull(n,:,i,j)=latstore(n,:);
            longfull(n,:,i,j)=longstore(n,:);
            utfull(n,:,i,j)=ut(n,:);
            vtfull(n,:,i,j)=vt(n,:);
            usfull(n,:,i,j)=us(n,:);
            vsfull(n,:,i,j)=vs(n,:);
        end
    end
end
%
V=zeros(nn,3,sx,sy);
Vd=V;
Vpp=V;
Vmp=V;
Vpm=V;
Vmm=V;
Vrp=V;
Vrm=V;
up=V;
um=V;
%
[ V(:,:,:,:) ] = windprofiles( vfull(:,:,:,:), rmfull(:,:,:,:), ...
    rfull(:,:,:,:), wprofile);
[ Vd(:,:,:,:) ] = windprofilem( vfull(:,:,:,:), rmfull(:,:,:,:), ...
    vsefull(:,:,:,:), rmsefull(:,:,:,:), rfull(:,:,:,:), wprofile);
[ Vpp(:,2,:,:) ] = windprofiles( vfull(:,3,:,:), rmfull(:,3,:,:), ...
    rfull(:,2,:,:)+deltar, wprofile);
[ Vpm(:,2,:,:) ] = windprofiles( vfull(:,3,:,:), rmfull(:,3,:,:), ...
    max(rfull(:,2,:,:)-deltar,0), wprofile);
[ Vmp(:,2,:,:) ] = windprofiles( vfull(:,1,:,:), rmfull(:,1,:,:), ...
    rfull(:,2,:,:)+deltar, wprofile);
[ Vmm(:,2,:,:) ] = windprofiles( vfull(:,1,:,:), rmfull(:,1,:,:), ...
    max(rfull(:,2,:,:)-deltar,0), wprofile);
[ Vrp(:,:,:,:) ] = windprofiles( vfull(:,:,:,:), rmfull(:,:,:,:), ...
    rfull(:,:,:,:)+deltar, wprofile);
[ Vrm(:,:,:,:) ] = windprofiles( vfull(:,:,:,:), rmfull(:,:,:,:), ...
    max(rfull(:,:,:,:)-deltar,0), wprofile);
%
% Convert to meters per second. Done now because windprofile expects knots.
%
V=knotfac*V;
Vd=knotfac*Vd;
Vpp=knotfac*Vpp;
Vpm=knotfac*Vpm;
Vmp=knotfac*Vmp;
Vmm=knotfac*Vmm;
Vrp=knotfac*Vrp;
Vrm=knotfac*Vrm;
%
latfac=latstore(1,1)./(abs(latstore(1,1))+1e-8);
vph=0.5.*(V+Vrp);
vmh=0.5.*(V+Vrm);
u1temp=vtfull.*costhetafull-utfull.*sinthetafull;
u2temp=vtfull.^2+utfull.^2;
vnetp=sqrt(vph.^2+2.*vph.*latfac.*u1temp+u2temp);
vnetm=sqrt(vmh.^2+2.*vmh.*latfac.*u1temp+u2temp);
%
for n=1:nn
    for j=1:3
         hfull(n,j,:,:)=h(:,:);
         hyfull(n,j,:,:)=hy(:,:);
         hxfull(n,j,:,:)=hx(:,:);
         cdfull(n,j,:,:)=cdrag(:,:);
         cdyfull(n,j,:,:)=cdy(:,:);
         cdxfull(n,j,:,:)=cdx(:,:);
    end 
end
%
cdfac=1e3.*0.5.*deltar.*(cdxfull.*costhetafull+cdyfull.*sinthetafull);
Cdp=cdfull+cdfac;
Cdm=cdfull-cdfac;
Cdp=max(Cdp,0);
Cdp=min(Cdp,0.005);
Cdm=max(Cdm,0);
Cdm=min(Cdm,0.005);
%
%  These lines of code added/modified  December 2018 to account for
%  roughness changes over water
%
facp=1+0.0193*vnetp;
facm=1+0.0193*vnetm;
facp(hfull>0)=1;  %  Do not let roughness change wind wind over land
facm(hfull>0)=1;
uekp=-Hi.*Cdp.*facp.*vph.*vnetp;
uekm=-Hi.*Cdm.*facm.*vmh.*vnetm;
%
cp=1000.*2.*Omega.*sin(pifac.*abs(latfull));
dMdrp=cp.*(rfull+0.5.*deltar)+(rfull+0.5.*deltar).*deltari.*(Vrp-V)+ ...
    0.5.*(Vrp+V);
dMdrp=max(dMdrp,2);
dMdrm=cp.*(rfull-0.5.*deltar)+(rfull-0.5.*deltar).*deltari.*(V-Vrm)+ ...
    0.5.*(Vrm+V);
dMdrm=max(dMdrm,2);
efacp=min((-1+2.*((rfull(:,2,:,:)+deltar)./rmfull(:,2,:,:)).^2),1);
efacm=min((-1+2.*((rfull(:,2,:,:)-deltar)./rmfull(:,2,:,:)).^2),1);
up(:,2,:,:)=(rfull(:,2,:,:)+deltar).*(-0.5.*timereswi.*efacp.* ...
    (Vpp(:,2,:,:)-Vmp(:,2,:,:))+uekp(:,2,:,:))./dMdrp(:,2,:,:);
um(:,2,:,:)=(rfull(:,2,:,:)-deltar).*(-0.5.*timereswi.*efacm.* ...
    (Vpm(:,2,:,:)-Vmm(:,2,:,:))+uekm(:,2,:,:))./dMdrm(:,2,:,:);
w(:,2,:,:)=-Htrop.*0.5*deltari.*((rfull(:,2,:,:)+deltar).*up(:,2,:,:)- ...
    (rfull(:,2,:,:)-deltar).*um(:,2,:,:))./max(rfull(:,2,:,:),1);
%
if se > 0   %  If secondary eyewalls present, add in their contribution to w
    %
    V=zeros(nn,3,sx,sy);
    Vpp=V;
    Vmp=V;
    Vpm=V;
    Vmm=V;
    Vrp=V;
    Vrm=V;
    up=V;
    um=V;
    [ V(:,:,:,:) ] = windprofiles( vsefull(:,:,:,:), rmsefull(:,:,:,:), ...
        rfull(:,:,:,:), wprofile);
    [ Vpp(:,2,:,:) ] = windprofiles( vsefull(:,3,:,:), rmsefull(:,3,:,:), ...
        rfull(:,2,:,:)+deltar, wprofile);
    [ Vpm(:,2,:,:) ] = windprofiles( vsefull(:,3,:,:), rmsefull(:,3,:,:), ...
        max(rfull(:,2,:,:)-deltar,0), wprofile);
    [ Vmp(:,2,:,:) ] = windprofiles( vsefull(:,1,:,:), rmsefull(:,1,:,:), ...
        rfull(:,2,:,:)+deltar, wprofile);
    [ Vmm(:,2,:,:) ] = windprofiles( vsefull(:,1,:,:), rmsefull(:,1,:,:), ...
        max(rfull(:,2,:,:)-deltar,0), wprofile);
    [ Vrp(:,:,:,:) ] = windprofiles( vsefull(:,:,:,:), rmsefull(:,:,:,:), ...
        rfull(:,:,:,:)+deltar, wprofile);
    [ Vrm(:,:,:,:) ] = windprofiles( vsefull(:,:,:,:), rmsefull(:,:,:,:), ...
        max(rfull(:,:,:,:)-deltar,0), wprofile); 
    %
    % Convert to meters per second. Done now because windprofile expects knots.
    %
    V=knotfac*V;
    Vpp=knotfac*Vpp;
    Vpm=knotfac*Vpm;
    Vmp=knotfac*Vmp;
    Vmm=knotfac*Vmm;
    Vrp=knotfac*Vrp;
    Vrm=knotfac*Vrm;    
    %
    vph=0.5.*(V+Vrp);
    vmh=0.5.*(V+Vrm);
    vnetp=sqrt(vph.^2+2.*vph.*latfac.*u1temp+u2temp);
    vnetm=sqrt(vmh.^2+2.*vmh.*latfac.*u1temp+u2temp);
    clear u1temp u2temp
    uekp=-Hi.*Cdp.*vph.*vnetp;
    uekm=-Hi.*Cdm.*vmh.*vnetm;
    %
    cp=1000.*2.*Omega.*sin(pifac.*abs(latfull));
    dMdrp=cp.*(rfull+0.5.*deltar)+(rfull+0.5.*deltar).*deltari.*(Vrp-V)+ ...
    0.5.*(Vrp+V);
    dMdrp=max(dMdrp,2);
    dMdrm=cp.*(rfull-0.5.*deltar)+(rfull-0.5.*deltar).*deltari.*(V-Vrm)+ ...
        0.5.*(Vrm+V);
    dMdrm=max(dMdrm,2);
    efacp=min((-1+2.*((rfull(:,2,:,:)+deltar)./rmsefull(:,2,:,:)).^2),1);
    efacm=min((-1+2.*((rfull(:,2,:,:)-deltar)./rmsefull(:,2,:,:)).^2),1);
    %
    %  Do not consider time rate of change when either velocity is zero
    %
    efacp=efacp.*min(Vpp(:,2,:,:),1).*min(Vmp(:,2,:,:),1);
    efacm=efacm.*min(Vpm(:,2,:,:),1).*min(Vmm(:,2,:,:),1);
    %
    up(:,2,:,:)=(rfull(:,2,:,:)+deltar).*(-0.5.*timereswi.*efacp.* ...
    (Vpp(:,2,:,:)-Vmp(:,2,:,:))+uekp(:,2,:,:))./dMdrp(:,2,:,:);
    um(:,2,:,:)=(rfull(:,2,:,:)-deltar).*(-0.5.*timereswi.*efacm.* ...
        (Vpm(:,2,:,:)-Vmm(:,2,:,:))+uekm(:,2,:,:))./dMdrm(:,2,:,:);
    %
    %  Do not calculate vertical velocities if either radial velocity is
    %  zero. This is necessary because secondary wind maximum can vanish
    %  from one time step to the next, so that time rate of change blows up
    %
    ufac=min(abs(30*up(:,2,:,:)),1).*min(abs(30*um(:,2,:,:)),1);
    %
    w2(:,2,:,:)=-Htrop*ufac*0.5*deltari.*((rfull(:,2,:,:)+deltar).*up(:,2,:,:)- ...
        (rfull(:,2,:,:)-deltar).*um(:,2,:,:))./max(rfull(:,2,:,:),1);
    %
    w=max(w,w2);
end    
%
%  Now add in topographic and shear components
%
clear um up efacm efacp dMdrp dMdrm
%hxmod=-0.0005.*cp.*vsfull;
%hymod=0.0005.*cp.*usfull;
%hxmod=-0.0005.*(cp+2.*V./(0.1+rfull)).*vsfull; % Modified August 8 2017 to include shear dotted with storm entropy
%hymod=0.0005.*(cp+2.*V./(0.1+rfull)).*usfull;
hxmod=-0.0005.*(2.*cp+2.*V./(0.1+rfull)+deltari*(Vrp-Vrm)).*vsfull; % Modified August 8 2017 to include shear dotted with storm entropy
hymod=0.0005.*(2.*cp+2.*V./(0.1+rfull)+deltari*(Vrp-Vrm)).*usfull;
ufunc=((radcity-rfull)./50);
ufunc=min(ufunc,1);
ufunc=max(ufunc,0);
utfull=utfull.*ufunc;
vtfull=vtfull.*ufunc;
w(:,2,:,:)=w(:,2,:,:)+(Vd(:,2,:,:).*latfac.*costhetafull(:,2,:,:)+vtfull(:,2,:,:)).* ...
    hyfull(:,2,:,:)+(utfull(:,2,:,:)-Vd(:,2,:,:).*latfac.*sinthetafull(:,2,:,:)).* ...
    hxfull(:,2,:,:)+Vd(:,2,:,:).*costhetafull(:,2,:,:).*hymod(:,2,:,:)-Vd(:,2,:,:).* ...
    sinthetafull(:,2,:,:).*hxmod(:,2,:,:);
w=min(w,7);
end