function [status]=create_wind_stress_forcing(Gname,Fname,wtime, Lon, Lat, windx, windy)

%function [status]=create_wind_stress_forcing(Gname,Fname,wtime, Lon, Lat, windx, windy)
%
% this function creates a forcing file of wind stress for ROMS
% from a user supplied wind field
% On Input:                                                            %
%                                                                      %
%    Gname       GRID NetCDF file name (string).                       %
%    Fname       FORCING NetCDF file (string).                         %    
%    wtime
%    lon  - not used 
%    lat  - not used 
%    windx - magnitude of wind in x-dir, same length as time
%    windy - magnitude of wind in y-dir, same length as time
%                                                                      %
% On Output:                                                           %
%                                                                      %
%    status      Error flag (integer).                                 %
%                                                                      %
% Calls:   MEXCDF (Interface to NetCDF library using Matlab).          %
%          nc_read, wind_stress                                        %
%                                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

method = 'linear';
num_time_steps=length(wtime);

%-----------------------------------------------------------------------
% Read in model grid positions and angle.
%-----------------------------------------------------------------------

%rlon=nc_read(Gname,'lon_rho');
%rlat=nc_read(Gname,'lat_rho');
%angle=nc_read(Gname,'angle');
nc = netcdf(Gname);
rlon=nc{'lon_rho'}(:);
rlat=nc{'lat_rho'}(:);
angle=nc{'angle'}(:);
close(nc)

[Lp,Mp]=size(rlon);
L=Lp-1;
M=Mp-1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create FORCING NetCDF file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nc = netcdf(Fname, 'clobber');
if isempty(nc), return, end
nc.type = ncchar('Gridpak file');
nc.gridid = 'combined grid';
nc.history = ncchar(['Created by "' mfilename '" on ' datestr(now)]);

nc.CPP_options = ncchar('DCOMPLEX, DBLEPREC, NCARG_32, PLOTS,');
name(nc.CPP_options, 'CPP-options')

%define dimensions
%nc('num_time_steps')= length(wtime);
nc('sms_time')= length(wtime);
nc('xi_psi') = L;
nc('xi_rho') = Lp;
nc('xi_u') = L;
nc('xi_v') = Lp;

nc('eta_psi') = M;
nc('eta_rho') = Mp;
nc('eta_u') = Mp;
nc('eta_v') = M;


%define variables
%nc{'sms_time'} = ncdouble('num_time_steps');
nc{'sms_time'} = ncdouble('sms_time');
nc{'sms_time'}.long_name = ncchar('surface stress time');
nc{'sms_time'}.units = ncchar('s');
nc{'sms_time'}.field = ncchar('scalar');

%nc{'sustr'} = ncdouble('num_time_steps','eta_u','xi_u'); %% 2 elements.
nc{'sustr'} = ncdouble('sms_time','eta_u','xi_u'); %% 2 elements.
nc{'sustr'}.long_name = ncchar('surface stress in xi dir');
nc{'sustr'}.units = ncchar('N/m^2');

%nc{'svstr'} = ncdouble('num_time_steps','eta_v','xi_v'); %% 2 elements.
nc{'svstr'} = ncdouble('sms_time','eta_v','xi_v'); %% 2 elements.
nc{'svstr'}.long_name = ncchar('surface stress in etq dir');
nc{'svstr'}.units = ncchar('N/m^2');

close(nc)

%-----------------------------------------------------------------------
% Process wind.
%-----------------------------------------------------------------------
nc = netcdf(Fname, 'write');

for n=1:num_time_steps,

% Compute wind stress.
%    [TauX,TauY]=wind_stress(windx(n),windy(n));
    [TauX,TauY]=wind_stress(squeeze(windx(n,:,:))',squeeze(windy(n,:,:))');

% Interpolate to model grid.

%    [Ustr]=interp2(Lon,Lat,TauX,rlon,rlat,method);
%    [Vstr]=interp2(Lon,Lat,TauY,rlon,rlat,method);
    [Ustr]=zeros(size(rlon))+TauX;
    [Vstr]=zeros(size(rlon))+TauY;
    
% Rotate to model grid.

    taux=Ustr.*cos(angle)+Vstr.*sin(angle);
    tauy=Vstr.*cos(angle)-Ustr.*sin(angle);
    
% Average to U- and V-points.

    sustr=0.5.*(taux(1:L,1:Mp)+taux(2:Lp,1:Mp));
    svstr=0.5.*(tauy(1:Lp,1:M)+tauy(1:Lp,2:Mp));

% Fill the arrays with data

    [status]=nc_write(Fname,'sustr',sustr,n);
    [status]=nc_write(Fname,'svstr',svstr,n);
    [status]=nc_write(Fname,'sms_time',wtime(n),n);

end

ncclose