function [x,y,windspeed] = windfieldx(nt,monthplot,dayplot,hourplot)
%
%  This script calculates the distribution of surface wind speed, in knots,
%  for a given storm at a specified time.
%
%  nt is the track number of the storm, monthplot is the month (0-12), day
%  plot is the day, and hourplot is the GMT time.
%
%  x and y are vectors containing the longitudes and latitudes of the grid
%
%  windspeed is the wind speed (knots) at each grid point.
%
%  The grid resolution and dimensions can be set automatically or
%  specified in the windfield map section of params.m.
%
%       Copyright 2011 WindRiskTech L.L.C.
%       Last modified May, 2019
%
%------------------------------------------------------------------------
%
params    %  Load parameters
%
load temp vstore vsestore latstore longstore rmstore rmsestore uinc vinc
load temp monthstore hourstore daystore
%
[nrm,mrm]=size(rmstore);
rfac=magfac.*(1+zeros(nrm,mrm));
%  
% Randomize radii of maximum wind
%
if strcmp(randfac,'y')
    rfac=magfac.*outer(nrm,mrm);
end
%
[~,jmaxd]=min(vstore,[],2);
jmaxd=jmaxd-1;
%
dum=zeros(1,jmaxd(nt));
V=[dum+2000; monthstore(nt,1:jmaxd(nt)); daystore(nt,1:jmaxd(nt));...
    hourstore(nt,1:jmaxd(nt)); dum; dum];
V=transpose(V);
timev=datenum(V);
timeplot=datenum([2000 monthplot dayplot hourplot 0 0]);
[~,jplot]=min(max(timeplot-timev,0));
jplot=max(jplot,2);
weight=(timeplot-timev(jplot-1))./(timev(jplot)-timev(jplot-1));
%
latstorm=weight.*latstore(nt,jplot)+(1-weight).*latstore(nt,jplot-1);
longstorm=weight.*longstore(nt,jplot)+(1-weight).*longstore(nt,jplot-1);
vstorm=weight.*vstore(nt,jplot)+(1-weight).*vstore(nt,jplot-1);
rmstorm=(weight.*rfac(nt,jplot).*rmstore(nt,jplot)+(1-weight).* ...
    rfac(nt,jplot-1).*rmstore(nt,jplot-1));
if exist('vsestore','var')
    vsestorm=weight.*vsestore(nt,jplot)+(1-weight).*vsestore(nt,jplot-1);
    rmsestorm=(weight.*rfac(nt,jplot).*rmsestore(nt,jplot)+(1-weight).* ...
        rfac(nt,jplot-1).*rmsestore(nt,jplot-1));
end
utstorm=weight.*uinc(nt,jplot)+(1-weight).*uinc(nt,jplot-1);
vtstorm=weight.*vinc(nt,jplot)+(1-weight).*vinc(nt,jplot-1);
%
pifac=acos(-1)/180;
dfac=60.*1.852;
%
if strcmp(bound,'auto')
    bxmin=floor(longstorm-deltax);
    bxmax=ceil(longstorm+deltax);
    bymin=floor(latstorm-deltay);
    bymax=ceil(latstorm+deltay);
end
%
x=bxmin:dellatlong:bxmax;
y=bymin:dellatlong:bymax;
sx=max(size(x));
sy=max(size(y));
%
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(sx,sy);
    for i=1:sx
        for j=1:sy
            ib=1+floor(4.*x(i));
            ibp=ib+1;
            if ibp > 1440
                ibp=1;
            end    
            jb=1+floor(4.*(y(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.*(y(j)+90)-jb;
            delx=1+4.*x(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 
%
% Calculate distance of each POI from track
%
dx=zeros(sx,sy);
dy=dx;
for i=1:sx
    for j=1:sy
        dx(i,j)=dfac.*cos(pifac.*latstorm).*(x(i)-longstorm);
        dy(i,j)=dfac.*(y(j)-latstorm);
    end
end    
r=sqrt(dx.^2+dy.^2);
if exist('vsestore','var')
    vcirc=windprofilem(vstorm,rmstorm,vsestorm,rmsestorm,r,wprofile);
else
    vcirc=windprofiles(vstorm,rmstorm,r,wprofile);
end 
vcirc=vcirc.*latstorm./(abs(latstorm)+1e-8);
un=-vcirc.*dy./max(r,1e-8)+utstorm;
vn=vcirc.*dx./max(r,1e-8)+vtstorm;
%
rfac=r./sqrt(r.^2+4.*radcity.^2);
windspeed=sqrt(un.^2+vn.^2)-rfac.*sqrt(utstorm.^2+vtstorm.^2);
%
if strcmp(wreduc,'y')
    cdfac=max(1+sqrt(cdrag).*logfac,0);
    windspeed=windspeed.*rat.*cdfac;
end
windspeed=max(windspeed,0);
%        
windspeed=transpose(windspeed);
% 
end

