function ier = write_roms_init_forcing(Gname,Fname,ubar,vbar,u,v)
%
% Create initial condition file for ROMS
%
% Input:                                                               %
%                                                                      %
%    Gname       GRID netCDF file name (string).                       %
%    Fname       FORCING metCDF file (string).                         %    
%    u,v         3D u-v momentum components (m/s)                      %
%    ubar,vbar   vertically integrated u-v momentum components (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.
%-----------------------------------------------------------------------
ncG = netcdf(Gname);
rlon=ncG{'lon_rho'}(:);
close(ncG)

[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('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{'ubar'} = ncdouble('ocean_time', 'eta_u', 'xi_u');
nc{'ubar'}.long_name = ncchar('vertically integrated u-momentum component');
nc{'ubar'}.units = ncchar('meter second-1');
nc{'ubar'}.time = ncchar('ocean_time');
nc{'ubar'}.coordinates = ncchar('x_u y_u');
nc{'ubar'}.field = ncchar(''ubar-velocity, scalar, series'');
 
nc{'vbar'} = ncdouble('ocean_time', 'eta_v', 'xi_v');
nc{'vbar'}.long_name = ncchar('vertically integrated v-momentum component');
nc{'vbar'}.units = ncchar('meter second-1');
nc{'vbar'}.time = ncchar('ocean_time');
nc{'vbar'}.coordinates = ncchar('x_v y_v');
nc{'vbar'}.field = ncchar('vbar-velocity, scalar, series');
 
nc{'u'} = ncdouble('ocean_time', 's_rho', 'eta_u', 'xi_u');
nc{'u'}.long_name = ncchar('u-momentum component');
nc{'u'}.units = ncchar('meter second-1');
nc{'u'}.time = ncchar('ocean_time');
nc{'u'}.coordinates = ncchar('x_u y_u');
nc{'u'}.field = ncchar(''u-velocity, scalar, series'');
 
nc{'v'} = ncdouble('ocean_time', 's_rho', 'eta_v', 'xi_v'); 
nc{'v'}.long_name = ncchar('v-momentum component');
nc{'v'}.units = ncchar('meter second-1');
nc{'v'}.time = ncchar('ocean_time');
nc{'v'}.coordinates = ncchar('x_v y_v');
nc{'v'}.field = ncchar('v-velocity, scalar, series');

close(nc)

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

Hwave = Hwave.* ones(size(rlon));
Dwave = Dwave.* ones(size(rlon));
Pwave = Pwave.* ones(size(rlon));
Lwave = Lwave.* ones(size(rlon));
ier.status(1)=nc_write(Fname,'Hwave',Hwave);
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