function ier = write_roms_wave_stress_forcing(Gname,Fname,wtime,Hwave,Dwave,Pwave,Lwave)
%
% Create a forcing file of waves for ROMS
%
% Input:                                                               %
%                                                                      %
%    Gname       GRID netCDF file name (string).                       %
%    Fname       FORCING metCDF file (string).                         %    
%    Hwave       Wave height                                           %
%    Dwave       Wave direction                                        %
%    Pwave       Wave period                                           %
%                                                                      %
% On Output:                                                           %
%                                                                      %
%    ier         Error flag structure                                  %
%                                                                      %
% Calls:   MEXCDF (Interface to NetCDF library using Matlab).          %
%          nc_read                                                     %
%                                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-----------------------------------------------------------------------
% Read in model grid positions and angle.
%-----------------------------------------------------------------------

nc = netcdf(Gname);
rlon=nc_read(Gname,'lon_rho');
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
num_time_steps= length(wtime);
% nc('ocean_time')= length(wtime);
nc('ocean_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{'ocean_time'} = ncdouble('ocean_time');
nc{'ocean_time'}.long_name = ncchar('surface stress time');
nc{'ocean_time'}.units = ncchar('s');
nc{'ocean_time'}.field = ncchar('scalar');

nc{'Hwave'} = ncdouble('ocean_time', 'eta_rho', 'xi_rho'); 
nc{'Hwave'}.long_name = ncchar('wind-induced significant wave height');
nc{'Hwave'}.units = ncchar('meter');
nc{'Hwave'}.time = ncchar('ocean_time');
nc{'Hwave'}.coordinates = ncchar('x_rho y_rho');
nc{'Hwave'}.field = ncchar('Hwave, scalar, series');
 
nc{'Dwave'} = ncdouble('ocean_time', 'eta_rho', 'xi_rho'); 
nc{'Dwave'}.long_name = ncchar('wind-induced wave direction');
nc{'Dwave'}.units = ncchar('degrees');
nc{'Dwave'}.time = ncchar('ocean_time');
nc{'Dwave'}.coordinates = ncchar('x_rho y_rho');
nc{'Dwave'}.field = ncchar('Dwave, scalar, series');

nc{'Pwave_bot'} = ncdouble('ocean_time', 'eta_rho', 'xi_rho'); 
nc{'Pwave_bot'}.long_name = ncchar('wind-induced wave period');
nc{'Pwave_bot'}.units = ncchar('degrees');
nc{'Pwave_bot'}.time = ncchar('ocean_time');
nc{'Pwave_bot'}.coordinates = ncchar('x_rho y_rho');
nc{'Pwave_bot'}.field = ncchar('Pwave, scalar, series');

nc{'Lwave'} = ncdouble('ocean_time', 'eta_rho', 'xi_rho'); 
nc{'Lwave'}.long_name = ncchar('wind-induced wave length');
nc{'Lwave'}.units = ncchar('m');
nc{'Lwave'}.time = ncchar('ocean_time');
nc{'Lwave'}.coordinates = ncchar('x_rho y_rho');
nc{'Lwave'}.field = ncchar('Lwave, scalar, series');

close(nc)

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

for n=1:num_time_steps,
    Hw = Hwave(n).*ones(size(rlon));
    Dw = Dwave(n).*ones(size(rlon));
    Pw = Pwave(n).*ones(size(rlon));
    Lw = Lwave(n).*ones(size(rlon));
    ier.status(1)=nc_write(Fname,'Hwave',Hw,n);
    ier.status(2)=nc_write(Fname,'Dwave',Dw,n);
    ier.status(3)=nc_write(Fname,'Pwave_bot',Pw,n);
    ier.status(4)=nc_write(Fname,'Lwave',Lw,n);
end

% Dwave = Dwave.* ones(size(rlon));
% Pwave = Pwave.* ones(size(rlon));
% Lwave = Lwave.* ones(size(rlon));
% ier.status(2)=nc_write(Fname,'Dwave',Dwave);
% ier.status(3)=nc_write(Fname,'Pwave_bot',Pwave);
% ier.status(4)=nc_write(Fname,'Lwave',Lwave);
ier.status(5)=nc_write(Fname,'ocean_time',wtime);

ncclose