function ier = write_roms_wind_stress_forcing(Gname,Fname,wtime,wstrx,wstry)
%
% Create a forcing file of UNIFORM wind stress for ROMS
%
% Input:                                                               %
%                                                                      %
%    Gname       GRID netCDF file name (string).                       %
%    Fname       FORCING metCDF file (string).                         %    
%    wtime       Wind time                                             %
%    wstrx - comp. of wind stress in x-dir, same length as time (N/m2) %
%    wstry - comp. of wind stress in y-dir, same length as time (N/m2) %
%                                                                      %
% On Output:                                                           %
%                                                                      %
%    ier         Error flag structure                                  %
%                                                                      %
% Calls:   MEXCDF (Interface to NetCDF library using Matlab).          %
%          nc_read, wind_stress                                        %
%                                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%TODO - add time units

num_time_steps=length(wtime);
%-----------------------------------------------------------------------
% Read in model grid positions and angle.
%-----------------------------------------------------------------------

nc = netcdf(Gname);
rlon=nc_read(Gname,'lon_rho');
rlat=nc_read(Gname,'lat_rho');
angle=nc_read(Gname,'angle');
 

% 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 eta dir');
nc{'svstr'}.units = ncchar('N/m^2');

close(nc)

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

for n=1:num_time_steps,

% Uniform assignment model grid.

     Ustr =ones(size(rlon))*wstrx(n);
     Vstr =ones(size(rlon))*wstry(n);
    
% 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

    ier.status(1)=nc_write(Fname,'sustr',sustr,n);
    ier.status(2)=nc_write(Fname,'svstr',svstr,n);
    ier.status(3)=nc_write(Fname,'sms_time',wtime(n),n);

end

ncclose