function [vs, dir1, dayfine] = pointseriesn(latstore,longstore,vstore,rmstore, ...
    vsestore,rmsestore,yearstore,monthstore,daystore,hourstore,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);
[~,vend]=min(vstore,[],2);
vend=vend-1;
%
[nn,m]=size(ut);
[~,sx]=size(plong);
for i=1:sx
    if plong(i) < 0
        plong(i)=plong(i)+360;
    end
end    
%
vendr=repmat(vend,[1,sx]);
pifac=acos(-1)/180;
dfac=60*1.852;
%
if strcmp(wreduc,'y')
    logfac=log(wheight/500)/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./(1+sqrt(mincd)*logfac); % Factor to ensure correct wind speeds over water    
    %
    % Interpolate drag coefficient to POI
    %
    cdrag=zeros(1,sx);
    for i=1:sx
            ib=1+floor(4.*plong(i));
            ibp=ib+1;
            if ibp > 1440
                ibp=1;
            end    
            jb=1+floor(4.*(plat(i)+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(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)=1./(d1./b1+d2./b2+d3./b3+d4./b4); % Conservative estimate; weights small roughness
            %
    end
end 
%
% 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=min(jmin,vendr-delj);
jmin=max(jmin,1+delj);
jstart=jmin-delj;
jend=jmin+delj;
jtot=2.*delj+1;
clear vend vendr
%
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);
latshort=zeros(nn,jtot,sx);
longshort=zeros(nn,jtot,sx);
utshort=zeros(nn,jtot,sx);
vtshort=zeros(nn,jtot,sx);
monthshort=zeros(nn,jtot,sx);
dayshort=zeros(nn,jtot,sx);
hourshort=zeros(nn,jtot,sx);
days=zeros(nn,jtot,sx);
%
for i=1:sx
    for n=1:nn
        vshort(n,:,i)=vstore(n,jstart(n,i):jend(n,i));
        vseshort(n,:,i)=vsestore(n,jstart(n,i):jend(n,i));
        rmshort(n,:,i)=rmstore(n,jstart(n,i):jend(n,i));
        rmseshort(n,:,i)=rmsestore(n,jstart(n,i):jend(n,i));
        latshort(n,:,i)=latstore(n,jstart(n,i):jend(n,i));
        longshort(n,:,i)=longstore(n,jstart(n,i):jend(n,i));
        utshort(n,:,i)=ut(n,jstart(n,i):jend(n,i));
        vtshort(n,:,i)=vt(n,jstart(n,i):jend(n,i));
        monthshort(n,:,i)=monthstore(n,jstart(n,i):jend(n,i));
        dayshort(n,:,i)=daystore(n,jstart(n,i):jend(n,i));
        hourshort(n,:,i)=hourstore(n,jstart(n,i):jend(n,i));
    end
end
for i=1:sx
    for n=1:nn
        for k=1:jtot
            days(n,k,i)=datenum(double(yearstore(n)),monthshort(n,k,i),dayshort(n,k,i),hourshort(n,k,i),0,0);
        end
        [daymax,djmax]=max(days(n,:,i));
        if days(n,end,i) < days(n,1,i)
            days(n,djmax:end,i)=daymax;
        end    
    end
end    
%
%  Create high time=resolution series
%
vfine=zeros(nn,jfine,sx);
vsefine=zeros(nn,jfine,sx);
rmfine=zeros(nn,jfine,sx);
rmsefine=zeros(nn,jfine,sx);
latfine=zeros(nn,jfine,sx);
longfine=zeros(nn,jfine,sx);
utfine=zeros(nn,jfine,sx);
vtfine=zeros(nn,jfine,sx);
dayfine=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,:);
        vsefine(:,k,:)=(1-weight).*vseshort(:,j,:)+weight.*vseshort(:,j+1,:);
        rmfine(:,k,:)=(1-weight).*rmshort(:,j,:)+weight.*rmshort(:,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,:);
        dayfine(:,k,:)=(1-weight).*days(:,j,:)+weight.*days(:,j+1,:);
    end
end
k=k+1;
vfine(:,k,:)=vshort(:,jtot,:);
vsefine(:,k,:)=vseshort(:,jtot,:);
rmfine(:,k,:)=rmshort(:,jtot,:);
rmsefine(:,k,:)=rmseshort(:,jtot,:);
latfine(:,k,:)=latshort(:,jtot,:);
longfine(:,k,:)=longshort(:,jtot,:);
utfine(:,k,:)=utshort(:,jtot,:);
vtfine(:,k,:)=vtshort(:,jtot,:);
dayfine(:,k,:)=days(:,jtot,:);
clear *short
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 ] = windprofilem( vfine, rmfine, vsefine, rmsefine, rfine, wprofile); 
%
V=V.*latfine./(abs(latfine)+1e-8);
%
if strcmp(wreduc,'y')
    cdfine=zeros(nn,jfine,sx);
    for i=1:nn
        for j=1:jfine
            cdfine(i,j,:)=cdrag(1,:);
        end
    end    
    cdfac=max(1+sqrt(cdfine).*logfac,0);
    V=V.*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