function ier = write_roms_wave_ub_forcing(...
   Gname,Fname,wtime,Hwave,Dwave,Pwave,Lwave,Ub)
% write_roms_wave_ub_forcing - Create wave forcing file for ROMS
% ier = write_roms_wave_ub_forcing(...
%          Gname,Fname,wtime,Hwave,Dwave,Pwave,Lwave,Ub)
%
% Input:                                                               %
%                                                                      %
%    Gname        GRID netCDF file name (string).                      %
%    Fname        FORCING metCDF file (string).                        %
%    Hwave        Wave height (m)                                      %
%    Dwave        Wave direction                                       %
%    Pwave        Wave period (s)                                      %
%    Ub(optional) Wave-orbital velocity (m/s)                          %
%                                                                      %
% 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('wave_time')= length(wtime);
nc('wave_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{'wave_time'} = ncdouble('wave_time');
nc{'wave_time'}.long_name = ncchar('wave time');
nc{'wave_time'}.units = ncchar('days');
nc{'wave_time'}.field = ncchar('scalar');

nc{'Hwave'} = ncdouble('wave_time', 'eta_rho', 'xi_rho'); 
nc{'Hwave'}.long_name = ncchar('wave height');
nc{'Hwave'}.units = ncchar('meters');
nc{'Hwave'}.time = ncchar('wave_time');
nc{'Hwave'}.coordinates = ncchar('x_rho y_rho');
nc{'Hwave'}.field = ncchar('Hwave, scalar, series');

nc{'Dwave'} = ncdouble('wave_time', 'eta_rho', 'xi_rho'); 
nc{'Dwave'}.long_name = ncchar('wave direction');
nc{'Dwave'}.units = ncchar('degrees');
nc{'Dwave'}.time = ncchar('wave_time');
nc{'Dwave'}.coordinates = ncchar('x_rho y_rho');
nc{'Dwave'}.field = ncchar('Dwave, scalar, series');

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

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

if(exist('Ub','var'))
nc{'Ub_swan'} = ncdouble('wave_time', 'eta_rho', 'xi_rho'); 
nc{'Ub_swan'}.long_name = ncchar('Wave-orbital velocity');
nc{'Ub_swan'}.units = ncchar('m/s');
nc{'Ub_swan'}.time = ncchar('wave_time');
nc{'Ub_swan'}.coordinates = ncchar('x_rho y_rho');
nc{'Ub_swan'}.field = ncchar('Ub, scalar, series');
end
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);
    if(exist('Ub','var'))
       Ub_swan = Ub(n).*ones(size(rlon));
       ier.status(5)=nc_write(Fname,'Ub_swan',Ub_swan,n);
    end
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(6)=nc_write(Fname,'wave_time',wtime);

ncclose