function [w] = pointwshortn(latstore,longstore,vstore,rmstore, ...
    vsestore,rmsestore,ut,vt,us,vs,plat,plong,h,hx,hy,timeres)
%
%  This function take the nn X m 'store' matrices  and,
%  using the latitude (plat) and longitude (plong) vectors of points of 
%  interest, calculates time series of vertical velocity at the points of
%  interest. 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. 
%      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 include variabel surface roughness over
%      land
%      Revised August 2023 to correct Omega for use in dM/dr terms
%      Revised October 23 2024 to correct w and w2 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
ntime=timelength./timeres+1;
nsteps=round(2./timeres);
nstepsi=1./nsteps;
delj=floor(timelength./4);
deltari=1./deltar;
timeresi=1./(3600.*timeres);
Hi=1./Htrop;
%
[nn,m]=size(ut);
[~,sx]=size(plong);
[~,sy]=size(plat);
ngrid=0;
if sx == sy 
    ngrid=1;
    sy=1;
end    
w=zeros(nn,ntime,sx,sy);
w2=zeros(nn,ntime,sx,sy);
%
if plong(1) < 0
    plong=plong + 360;
end    
%
pifac=acos(-1)/180;
dfac=60.*1.852;
sfac=1./(0.25.*60.*1852);
%
%  Convert knots to meters per second. (Do not convert wind speeds until
%  later.)
%
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
%
% Three statements below added Dec. 2018 to reduce drag coefficient of
% near-coastal locations to reflect response time of landfalling TCs. 
%
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,jj)=dfac.*cos(pifac.*plat(j)).*(plong(i)-longstore);
        dy(:,:,i,jj)=dfac.*(plat(j)-latstore);
    end
end    
radius=sqrt(dx.*dx+dy.*dy);
radius=max(radius,0.5);
clear dx dy
[~,jmin]=min(radius,[],2);
jmin=max(jmin,1+delj);
jmin=min(jmin,m-delj);
jstart=jmin-delj;
jend=jmin+delj;
jtot=2.*delj+1;
%
jfine=1+nsteps.*(jtot-1);
%
%  Create reduced length time series of each quantity
%
vshort=zeros(nn,jtot,sx,sy);
rmshort=zeros(nn,jtot,sx,sy);
vseshort=zeros(nn,jtot,sx,sy);
rmseshort=zeros(nn,jtot,sx,sy);
rshort=zeros(nn,jtot,sx,sy);
latshort=zeros(nn,jtot,sx,sy);
longshort=zeros(nn,jtot,sx,sy);
utshort=zeros(nn,jtot,sx,sy);
vtshort=zeros(nn,jtot,sx,sy);
usshort=zeros(nn,jtot,sx,sy);
vsshort=zeros(nn,jtot,sx,sy);
%
for i=1:sx
    for j=1:sy
        for n=1:nn
            vshort(n,:,i,j)=vstore(n,jstart(n,1,i,j):jend(n,1,i,j));
            rmshort(n,:,i,j)=rmstore(n,jstart(n,1,i,j):jend(n,1,i,j));
            vseshort(n,:,i,j)=vsestore(n,jstart(n,1,i,j):jend(n,1,i,j));
            rmseshort(n,:,i,j)=rmsestore(n,jstart(n,1,i,j):jend(n,1,i,j));
            rshort(n,:,i,j)=radius(n,jstart(n,1,i,j):jend(n,1,i,j),i,j);
            latshort(n,:,i,j)=latstore(n,jstart(n,1,i,j):jend(n,1,i,j));
            longshort(n,:,i,j)=longstore(n,jstart(n,1,i,j):jend(n,1,i,j));
            utshort(n,:,i,j)=ut(n,jstart(n,1,i,j):jend(n,1,i,j));
            vtshort(n,:,i,j)=vt(n,jstart(n,1,i,j):jend(n,1,i,j));
            usshort(n,:,i,j)=us(n,jstart(n,1,i,j):jend(n,1,i,j));
            vsshort(n,:,i,j)=vs(n,jstart(n,1,i,j):jend(n,1,i,j));
        end
    end
end
%
%  Create high time-resolution series
%
vfine=zeros(nn,jfine,sx,sy);
rmfine=zeros(nn,jfine,sx,sy);
vsefine=zeros(nn,jfine,sx,sy);
rmsefine=zeros(nn,jfine,sx,sy);
latfine=zeros(nn,jfine,sx,sy);
longfine=zeros(nn,jfine,sx,sy);
utfine=zeros(nn,jfine,sx,sy);
vtfine=zeros(nn,jfine,sx,sy);
usfine=zeros(nn,jfine,sx,sy);
vsfine=zeros(nn,jfine,sx,sy);
hfine=zeros(nn,jfine,sx,sy);
hyfine=zeros(nn,jfine,sx,sy);
hxfine=zeros(nn,jfine,sx,sy);
cdfine=zeros(nn,jfine,sx,sy);
cdyfine=zeros(nn,jfine,sx,sy);
cdxfine=zeros(nn,jfine,sx,sy);
%
k=0;
for j=1:jtot-1
    for n=1:nsteps
        k=k+1;
        weight=(n-1).*nstepsi;
        vfine(:,k,:,:)=(1-weight).*vshort(:,j,:,:)+weight.*vshort(:,j+1,:,:);
        rmfine(:,k,:,:)=(1-weight).*rmshort(:,j,:,:)+weight.*rmshort(:,j+1,:,:);
        vsefine(:,k,:,:)=(1-weight).*vseshort(:,j,:,:)+weight.*vseshort(:,j+1,:,:);
        rmsefine(:,k,:,:)=(1-weight).*rmseshort(:,j,:,:)+weight.*rmseshort(:,j+1,:,:);
        latfine(:,k,:,:)=(1-weight).*latshort(:,j,:,:)+weight.*latshort(:,j+1,:,:);
        longfine(:,k,:,:)=(1-weight).*longshort(:,j,:,:)+weight.*longshort(:,j+1,:,:);
        utfine(:,k,:,:)=(1-weight).*utshort(:,j,:,:)+weight.*utshort(:,j+1,:,:);
        vtfine(:,k,:,:)=(1-weight).*vtshort(:,j,:,:)+weight.*vtshort(:,j+1,:,:);
        usfine(:,k,:,:)=(1-weight).*usshort(:,j,:,:)+weight.*usshort(:,j+1,:,:);
        vsfine(:,k,:,:)=(1-weight).*vsshort(:,j,:,:)+weight.*vsshort(:,j+1,:,:);
    end
end
k=k+1;
vfine(:,k,:,:)=vshort(:,jtot,:,:);
rmfine(:,k,:,:)=rmshort(:,jtot,:,:);
vsefine(:,k,:,:)=vseshort(:,jtot,:,:);
rmsefine(:,k,:,:)=rmseshort(:,jtot,:,:);
latfine(:,k,:,:)=latshort(:,jtot,:,:);
longfine(:,k,:,:)=longshort(:,jtot,:,:);
utfine(:,k,:,:)=utshort(:,jtot,:,:);
vtfine(:,k,:,:)=vtshort(:,jtot,:,:);
usfine(:,k,:,:)=usshort(:,jtot,:,:);
vsfine(:,k,:,:)=vsshort(:,jtot,:,:);
%
rmsefine=max(rmsefine,0.1);
clear *short dx dy
dx=zeros(nn,jfine,sx,sy);
dy=dx;
for i=1:sx
    for jj=1:sy
        if ngrid == 1
            j=i;
        else
            j=jj;
        end    
        dx(:,:,i,jj)=dfac.*cos(pifac.*plat(j)).*(plong(i)-longfine(:,:,i,jj));
        dy(:,:,i,jj)=dfac.*(plat(j)-latfine(:,:,i,jj));
    end
end    
rfine=sqrt(dx.^2+dy.^2);
%
V=zeros(nn,jfine,sx,sy);
Vd=V;
Vpp=V;
Vmp=V;
Vpm=V;
Vmm=V;
Vrp=V;
Vrm=V;
up=V;
um=V;
[ V(:,:,:,:) ] = windprofiles( vfine(:,:,:,:), rmfine(:,:,:,:), ...
    rfine(:,:,:,:), wprofile);
[ Vd(:,:,:,:) ] = windprofilem( vfine(:,:,:,:), rmfine(:,:,:,:), ...
        vsefine(:,:,:,:), rmsefine(:,:,:,:), rfine(:,:,:,:), wprofile);
[ Vpp(:,2:jfine-1,:,:) ] = windprofiles( vfine(:,3:jfine,:,:), rmfine(:,3:jfine,:,:), ...
    rfine(:,2:jfine-1,:,:)+deltar, wprofile);
[ Vpm(:,2:jfine-1,:,:) ] = windprofiles( vfine(:,3:jfine,:,:), rmfine(:,3:jfine,:,:), ...
    max(rfine(:,2:jfine-1,:,:)-deltar,0), wprofile);
[ Vmp(:,2:jfine-1,:,:) ] = windprofiles( vfine(:,1:jfine-2,:,:), rmfine(:,1:jfine-2,:,:), ...
    rfine(:,2:jfine-1,:,:)+deltar, wprofile);
[ Vmm(:,2:jfine-1,:,:) ] = windprofiles( vfine(:,1:jfine-2,:,:), rmfine(:,1:jfine-2,:,:), ...
    max(rfine(:,2:jfine-1,:,:)-deltar,0), wprofile);
[ Vrp(:,:,:,:) ] = windprofiles( vfine(:,:,:,:), rmfine(:,:,:,:), ...
    rfine(:,:,:,:)+deltar, wprofile);
[ Vrm(:,:,:,:) ] = windprofiles( vfine(:,:,:,:), rmfine(:,:,:,:), ...
    max(rfine(:,:,:,:)-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);
rfinei=1./max(rfine,1);
u1temp=vtfine.*dx.*rfinei-utfine.*dy.*rfinei;
u2temp=vtfine.^2+utfine.^2;
vnetp=sqrt(vph.^2+2.*vph.*latfac.*u1temp+u2temp);
vnetm=sqrt(vmh.^2+2.*vmh.*latfac.*u1temp+u2temp);
for i=1:nn
    for j=1:jfine
        hfine(i,j,:,:)=h(:,:);
        hxfine(i,j,:,:)=hx(:,:);
        hyfine(i,j,:,:)=hy(:,:);
        cdfine(i,j,:,:)=cdrag(:,:);
        cdxfine(i,j,:,:)=cdx(:,:);
        cdyfine(i,j,:,:)=cdy(:,:);

    end
end    
%
cdfac=1e3.*0.5.*deltar.*(cdxfine.*dx.*rfinei+cdyfine.*dy.*rfinei);
Cdp=cdfine+cdfac;
Cdm=cdfine-cdfac;
Cdp=max(Cdp,0);
Cdp=min(Cdp,0.005);
Cdm=max(Cdm,0);
Cdm=min(Cdm,0.005);
%{
facp=1+0.0193*vnetp;
facm=1+0.0193*vnetm;
facp(hfine>0)=1;  %  Do not let roughness change wind wind over land
facm(hfine>0)=1;
uekp=-Hi.*Cdp.*facp.*vph.*vnetp;
uekm=-Hi.*Cdm.*facm.*vmh.*vnetm;
%}
uekp=-Hi.*Cdp.*vph.*vnetp;
uekm=-Hi.*Cdm.*vmh.*vnetm;
%
cp=1000.*2.*Omega.*sin(pifac.*abs(latfine));
dMdrp=cp.*(rfine+0.5.*deltar)+(rfine+0.5.*deltar).*deltari.*(Vrp-V)+ ...
    0.5.*(Vrp+V);
dMdrp=max(dMdrp,10);
dMdrm=cp.*(rfine-0.5.*deltar)+(rfine-0.5.*deltar).*deltari.*(V-Vrm)+ ...
    0.5.*(Vrm+V);
dMdrm=max(dMdrm,10);
efacp=min((-1+2.*((rfine(:,2:jfine-1,:,:)+deltar)./rmfine(:,2:jfine-1,:,:)).^2),1);
efacm=min((-1+2.*((rfine(:,2:jfine-1,:,:)-deltar)./rmfine(:,2:jfine-1,:,:)).^2),1);
up(:,2:jfine-1,:,:)=(rfine(:,2:jfine-1,:,:)+deltar).*(-0.5.*timeresi.* ...
    efacp.*(Vpp(:,2:jfine-1,:,:)-Vmp(:,2:jfine-1,:,:))+uekp(:,2:jfine-1,:,:))./ ...
    dMdrp(:,2:jfine-1,:,:);
um(:,2:jfine-1,:,:)=(rfine(:,2:jfine-1,:,:)-deltar).*(-0.5.*timeresi.* ...
    efacm.*(Vpm(:,2:jfine-1,:,:)-Vmm(:,2:jfine-1,:,:))+uekm(:,2:jfine-1,:,:))./ ...
    dMdrm(:,2:jfine-1,:,:);
w(:,2:jfine-1,:,:)=-Htrop.*0.5*deltari.*((rfine(:,2:jfine-1,:,:)+deltar).*up(:,2:jfine-1,:,:)- ...
    (rfine(:,2:jfine-1,:,:)-deltar).*um(:,2:jfine-1,:,:)).*rfinei(:,2:jfine-1,:,:);
%
if se > 0   %  If secondary eyewalls present, add in their contribution to w
    %
    V=zeros(nn,jfine,sx,sy);
    Vpp=V;
    Vmp=V;
    Vpm=V;
    Vmm=V;
    Vrp=V;
    Vrm=V;
    up=V;
    um=V;
    [ V(:,:,:,:) ] = windprofiles( vsefine(:,:,:,:), rmsefine(:,:,:,:), ...
        rfine(:,:,:,:), wprofile);
    [ Vpp(:,2:jfine-1,:,:) ] = windprofiles( vsefine(:,3:jfine,:,:), rmsefine(:,3:jfine,:,:), ...
        rfine(:,2:jfine-1,:,:)+deltar, wprofile);
    [ Vpm(:,2:jfine-1,:,:) ] = windprofiles( vsefine(:,3:jfine,:,:), rmsefine(:,3:jfine,:,:), ...
        max(rfine(:,2:jfine-1,:,:)-deltar,0), wprofile);
    [ Vmp(:,2:jfine-1,:,:) ] = windprofiles( vsefine(:,1:jfine-2,:,:), rmsefine(:,1:jfine-2,:,:), ...
        rfine(:,2:jfine-1,:,:)+deltar, wprofile);
    [ Vmm(:,2:jfine-1,:,:) ] = windprofiles( vsefine(:,1:jfine-2,:,:), rmsefine(:,1:jfine-2,:,:), ...
        max(rfine(:,2:jfine-1,:,:)-deltar,0), wprofile);
    [ Vrp(:,:,:,:) ] = windprofiles( vsefine(:,:,:,:), rmsefine(:,:,:,:), ...
        rfine(:,:,:,:)+deltar, wprofile);
    [ Vrm(:,:,:,:) ] = windprofiles( vsefine(:,:,:,:), rmsefine(:,:,:,:), ...
        max(rfine(:,:,:,:)-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(latfine));
    dMdrp=cp.*(rfine+0.5.*deltar)+(rfine+0.5.*deltar).*deltari.*(Vrp-V)+ ...
        0.5.*(Vrp+V);
    dMdrp=max(dMdrp,10);
    dMdrm=cp.*(rfine-0.5.*deltar)+(rfine-0.5.*deltar).*deltari.*(V-Vrm)+ ...
        0.5.*(Vrm+V);
    dMdrm=max(dMdrm,10);
    efacp=min((-1+2.*((rfine(:,2:jfine-1,:,:)+deltar)./rmsefine(:,2:jfine-1,:,:)).^2),1);
    efacm=min((-1+2.*((rfine(:,2:jfine-1,:,:)-deltar)./rmsefine(:,2:jfine-1,:,:)).^2),1);
    %
    %  Do not consider time rate of change when either velocity is zero
    %
    efacp=efacp.*min(Vpp(:,2:jfine-1,:,:),1).*min(Vmp(:,2:jfine-1,:,:),1);
    efacm=efacm.*min(Vpm(:,2:jfine-1,:,:),1).*min(Vmm(:,2:jfine-1,:,:),1);
    %
    up(:,2:jfine-1,:,:)=(rfine(:,2:jfine-1,:,:)+deltar).*(-0.5.*timeresi.* ...
        efacp.*(Vpp(:,2:jfine-1,:,:)-Vmp(:,2:jfine-1,:,:))+uekp(:,2:jfine-1,:,:))./ ...
        dMdrp(:,2:jfine-1,:,:);
    um(:,2:jfine-1,:,:)=(rfine(:,2:jfine-1,:,:)-deltar).*(-0.5.*timeresi.* ...
        efacm.*(Vpm(:,2:jfine-1,:,:)-Vmm(:,2:jfine-1,:,:))+uekm(:,2:jfine-1,:,:))./ ...
        dMdrm(:,2:jfine-1,:,:);
    %
    %  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:jfine-1,:,:)),1).*min(abs(30*um(:,2:jfine-1,:,:)),1);
    %
    w2(:,2:jfine-1,:,:)=-ufac.*Htrop.*0.5*deltari.*((rfine(:,2:jfine-1,:,:)+deltar).*up(:,2:jfine-1,:,:)- ...
        (rfine(:,2:jfine-1,:,:)-deltar).*um(:,2:jfine-1,:,:)).*rfinei(:,2:jfine-1,:,:);
    %
    w=max(w,w2);
end    
%
%  Now add in topographic and shear components
%
clear um up efacm efacp dMdrp dMdrm
%hxmod=-0.0005.*cp.*vsfine;
%hymod=0.0005.*cp.*usfine;
hxmod=-0.0005*(2.*cp+2.*V./(0.1+rfine)+deltari*(Vrp-Vrm)).*vsfine; % Modified August 8 2017 to include shear dotted with storm entropy
hymod=0.0005.*(2.*cp+2.*V./(0.1+rfine)+deltari*(Vrp-Vrm)).*usfine;

%
% Reduce effect of translation speed outside of radcity
%
ufunc=((radcity-rfine)./50);
ufunc=min(ufunc,1);
ufunc=max(ufunc,0);
utfine=utfine.*ufunc;
vtfine=vtfine.*ufunc;
%
% Reduce effect of orography outside of storm core
%
ufunc=(150-rfine)./30;
ufunc=min(ufunc,0.6);
ufunc=max(ufunc,0.2);
hxfine=hxfine.*ufunc;
hyfine=hyfine.*ufunc;
%
w=w+(Vd.*latfac.*dx.*rfinei+vtfine).*hyfine+(utfine-Vd.*latfac.*dy.*rfinei).* ...
    hxfine+Vd.*dx.*rfinei.*hymod-Vd.*dy.*rfinei.*hxmod;
w=min(w,7);
%
end