function [wq,dayk] = pointwshortnqdx(latstore,longstore,dq,vstore,rmstore, ...
    datestore, ...
    vsestore,rmsestore,ut,vt,us,vs,plat,plong,h,hx,hy,timeres,wrad)
%
%  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 account for variable surface roughness
%      Revised August 2023 to correct Omega for use in dM/dr terms
%      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=max(size(plong));  
w=zeros(nn,ntime,sx);
w2=zeros(nn,ntime,sx);
%
for i=1:sx
    if plong(i) < 0
        plong(i)=plong(i) + 360;
    end    
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
%
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(1,sx);
cdx=zeros(1,sx);
cdy=zeros(1,sx);
for i=1:sx
    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(i)+90));
    jbs=1+floor(4.*(plat(i)+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(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;
    cdrag(i)=d1*b1+d2*b2+d3*b3+d4*b4;
    dely=0.5+4.*(plat(i)+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)=(d1*b1x+d2*b2x+d3*b3x+d4*b4x)./cos(pifac.*plat(i));
    cdy(i)=d1*b1y+d2*b2y+d3*b3y+d4*b4y;
end
%
% 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);
dy=dx;
for i=1:sx
    dx(:,:,i)=dfac.*cos(pifac.*plat(i)).*(plong(i)-longstore);
    dy(:,:,i)=dfac.*(plat(i)-latstore);
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);
rmshort=zeros(nn,jtot,sx);
vseshort=zeros(nn,jtot,sx);
rmseshort=zeros(nn,jtot,sx);
rshort=zeros(nn,jtot,sx);
latshort=zeros(nn,jtot,sx);
longshort=zeros(nn,jtot,sx);
dateshort=zeros(nn,jtot,sx); 
utshort=zeros(nn,jtot,sx);
vtshort=zeros(nn,jtot,sx);
usshort=zeros(nn,jtot,sx);
vsshort=zeros(nn,jtot,sx);
dqshort=zeros(nn,jtot,sx);
%
for i=1:sx
    for n=1:nn
        vshort(n,:,i)=vstore(n,jstart(n,1,i):jend(n,1,i));
        rmshort(n,:,i)=rmstore(n,jstart(n,1,i):jend(n,1,i));
        vseshort(n,:,i)=vsestore(n,jstart(n,1,i):jend(n,1,i));
        rmseshort(n,:,i)=rmsestore(n,jstart(n,1,i):jend(n,1,i));
        rshort(n,:,i)=radius(n,jstart(n,1,i):jend(n,1,i),i);
        latshort(n,:,i)=latstore(n,jstart(n,1,i):jend(n,1,i));
        longshort(n,:,i)=longstore(n,jstart(n,1,i):jend(n,1,i));
        dateshort(n,:,i)=datestore(n,jstart(n,1,i):jend(n,1,i));
        utshort(n,:,i)=ut(n,jstart(n,1,i):jend(n,1,i));
        vtshort(n,:,i)=vt(n,jstart(n,1,i):jend(n,1,i));
        usshort(n,:,i)=us(n,jstart(n,1,i):jend(n,1,i));
        vsshort(n,:,i)=vs(n,jstart(n,1,i):jend(n,1,i));
        dqshort(n,:,i)=dq(n,jstart(n,1,i):jend(n,1,i));
    end
end
%
%  Create high time-resolution series
%
vfine=zeros(nn,jfine,sx);
rmfine=zeros(nn,jfine,sx);
vsefine=zeros(nn,jfine,sx);
rmsefine=zeros(nn,jfine,sx);
latfine=zeros(nn,jfine,sx);
longfine=zeros(nn,jfine,sx);
dayk=zeros(nn,jfine,sx); 
utfine=zeros(nn,jfine,sx);
vtfine=zeros(nn,jfine,sx);
usfine=zeros(nn,jfine,sx);
vsfine=zeros(nn,jfine,sx);
hfine=zeros(nn,jfine,sx);
hyfine=zeros(nn,jfine,sx);
hxfine=zeros(nn,jfine,sx);
cdfine=zeros(nn,jfine,sx);
cdyfine=zeros(nn,jfine,sx);
cdxfine=zeros(nn,jfine,sx);
dqfine=zeros(nn,jfine,sx);
%
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,:);
        dayk(:,k,:)=(1-weight).*dateshort(:,j,:)+weight.*dateshort(:,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,:);
        dqfine(:,k,:)=(1-weight).*dqshort(:,j,:)+weight.*dqshort(:,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,:);
dayk(:,k,:)=dateshort(:,jtot,:); 
utfine(:,k,:)=utshort(:,jtot,:);
vtfine(:,k,:)=vtshort(:,jtot,:);
usfine(:,k,:)=usshort(:,jtot,:);
vsfine(:,k,:)=vsshort(:,jtot,:);
dqfine(:,k,:)=dqshort(:,jtot,:);
%
rmsefine=max(rmsefine,0.1);
clear *short dx dy
dx=zeros(nn,jfine,sx);
dy=dx;
for i=1:sx      
    dx(:,:,i)=dfac.*cos(pifac.*plat(i)).*(plong(i)-longfine(:,:,i));
    dy(:,:,i)=dfac.*(plat(i)-latfine(:,:,i));
end    
rfine=sqrt(dx.^2+dy.^2);
%
V=zeros(nn,jfine,sx);
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);
    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);
wq=dqfine.*max(w-wrad,0);
%
end