function [vs, dir1] = pointshortn(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 January, 2019
%-------------------------------------------------------------------------
%
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;
%
if strcmp(wreduc,'y')
    logfac=log(500/wheight)/0.35; % Factor to Reduce gradient wind to 10 m wind
    %
    % Load neutral drag coefficients
    %
    load C_Drag500 cd
    load bathymetry bathy
    mincd=min(min(cd));
    cd(bathy<0)=mincd; % This statement added to set ocean point drag coefficients to 
                       % their proper value; this treats a few bogus points
                       % in the drag coefficient data set resulting from
                       % ECMWF's interpolation from a 1 x 1 grid
    clear bathy
    rat=1+sqrt(mincd)*logfac; % Factor to ensure correct wind speeds over water    
    %
    % Interpolate drag coefficient to POI
    %
    cdrag=zeros(sx,sy);
    for i=1:sx
        for j=1:sy
            ib=1+floor(4.*plong(i));
            ibp=ib+1;
            if ib > 1440
                ib=ib-1440;
            end
            if ibp > 1440
                ibp=ibp-1440;
            end    
            jb=1+floor(4.*(plat(j)+90));
            b1=cd(ib,jb);
            b2=cd(ib,jb+1);
            b3=cd(ibp,jb);
            b4=cd(ibp,jb+1);
            %
            %  This block interpolates drag coefficient to plong, plat
            %
            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)=1./(d1./b1+d2./b2+d3./b3+d4./b4); % Conservative estimate; weights small roughness
            %
        end
    end
end 
%disp(num2str(cdrag(1,1)))
%
% 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;
jstart=max(jstart,1);
jend=jmin+delj;
jtot=jend(1)-jstart(1)+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 ] = windprofilem( vfine, rmfine, vsefine, rmsefine, rfine, wprofile); 
%
V=V.*sign(latfine);
%
if strcmp(wreduc,'y')
    cdfine=zeros(nn,jfine,sx,sy);
    for i=1:nn
        for j=1:jfine
            cdfine(i,j,:,:)=cdrag(:,:);
        end
    end    
    cdfac=1./(1+max(sqrt(cdfine).*logfac,0));
    V=V.*rat.*cdfac;
    utfine=utfine.*rat.*cdfac;
    vtfine=vtfine.*rat.*cdfac;
end    
%
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