function [vs, dir1] = pointshort(latstore,longstore,vstore,rmstore, ...
   vsestore,rmsestore,ut,vt,plat,plong,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 wind speed and direction at the 
%  points of interest. The routine utrans.m should be run before this one. 
%      Copyright 2009 WindRiskTech L.L.C.
%      Last revised October, 2014
%-------------------------------------------------------------------------
%
params    %  Load parameters
%
nsteps=round(2./timeres);
nstepsi=1./nsteps;
delj=floor(timelength./4);
%
[nn,m]=size(ut);
[~,sx]=size(plong);
[~,sy]=size(plat);
ngrid=0;
if sx == sy 
    ngrid=1;
    sy=1;
end    
%
if plong(1) < 0
    plong=plong + 360;
end    
%
pifac=acos(-1)/180;
dfac=60.*1.852;
%
% 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);
latshort=zeros(nn,jtot,sx,sy);
longshort=zeros(nn,jtot,sx,sy);
utshort=zeros(nn,jtot,sx,sy);
vtshort=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));
            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));
        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);
%
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,:,:);
    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,:,:);
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 ] = windprofiledouble( vfine, rmfine, vsefine, rmsefine, rfine, wprofile); 
%
V=V.*latfine./(abs(latfine)+1e-8);
vn=vtfine+V.*dx./max(rfine,0.5);
un=utfine-V.*dy./max(rfine,0.5);
rfac=rfine./sqrt(rfine.^2+radcity.^2);
vs=sqrt(un.^2+vn.^2)-rfac.*sqrt(utfine.^2+vtfine.^2);
vs=max(vs,0);
tempd=360+atan2(-un,-vn)./pifac;
dir1=mod(tempd,360);
%
end