%function wavebuoy2forceroms(grdfile,f291file,timestruct,plot_key);
grdfile = 'D:\crs\proj\pv\PV2004_data\roms_stuff\pv_grd2b.nc';

%TODO - add ability to vary water level
%TODO - change f291file to netCDF

% timestruct shoud contain
%   ts.start
%   ts.end
%   ts.dt 
ts.start = julian(2002,11,1,0);
ts.end = julian(2002,12,1,0);
ts.dt = 1/24;

wave_jd = (ts.start:ts.dt:ts.end)'; % this still needs to be converted to ROMS time
nt = length(wave_jd)
gregorian( wave_jd(1) )
gregorian( wave_jd(2) )
gregorian( wave_jd(end) )

% these are required
write_hsig = 1;
write_xp = 1;
write_yp = 1;
write_depth = 1;

% these are standard
% existence of these files is not tested
write_rtp = 1;
write_wdir = 1;
write_ubot = 1;
write_tmbot = 1;

% these are optional and only written if the corresponding SWAN output file
% exists
write_wlen = 0;
write_wind = 0;
write_dissip = 0;
write_force = 0;
write_setup = 0;
write_vel = 0;

%% Open the .grd file for dimensions
% ncload(grdfile);
%%
ncg = netcdf(grdfile,'r');
rlat = ncg{'lat_rho'}(:,:);
rlon = ncg{'lon_rho'}(:,:);
rmask = ncg{'mask_rho'}(:,:);
h = ncg{'h'}(:,:);
[ysize,xsize]=size(h);

%% Open the f291 file
% TODO  - this should read a .nc file instead
load nov_2002 % contains structures S and m
f = S(100).f;
df = S(100).df;

jd = [S.jd]';
hsig = [S.hsig]';
tdom = [S.tdom]';
if(any(isnan(jd))),
   fprintf(1,'Warning % NaNs found in wave time.',sum(isnan(jd)))
end
hsig = fillnan(jd,hsig);
tdom = fillnan(jd,tdom);
ok = find( ~isnan( jd + hsig + tdom ) );
jd = jd(ok);
hsig = hsig(ok);
tdom = tdom(ok);
hsig = interp1(jd,hsig,wave_jd,'nearest');
tdom = interp1(jd,tdom,wave_jd,'nearest');

%%
depth = h;
dum = ones(nt,ysize,xsize);
%hsig = reshape(hsig,[1 1 nt]);
Hsig = repmat( hsig,[1 ysize xsize] );
%tdom = reshape(hsig,[1 1 nt]);
Tdom = repmat( tdom, [1 ysize xsize] );
Tmbot = Tdom;
Ubot = 0.*dum;
Dwave = 270.*dum;
Lwave = 0.1*dum;

%%
for j=1:xsize,
   for i=1:ysize,
      if( (h(i,j)<200) &  ( rmask(i,j)==1 ) ),
         for k=1:nt
            if( sum( S(k).sd > 0 ) ),
               %            Ubot(k,i,j) = ubqf( Tdom(k,i,j), Hsig(k,i,j), h(i,j) );
               kk = tind( [S.jd], gregorian( wave_jd(k) ) );
               ubs(k) = ubspec( h(i,j), S(kk).sd, f, 'default', df );
               Ubot(k,i,j) = ubs.ubr;
               Tmbot(k,i,j)= ubs.Tr;
            end
         end
      end
   end
end

%% Create netcdf file.
fname = 'fake_wave.nc';
nc=netcdf(fname,'clobber');
if isempty(nc), return, end

%% Global attributes:

disp(' ## Defining Global Attributes...')
nc.history = ncchar(['Created by wavebuoy2forceroms on ' datestr(now)]);
nc.type = ncchar('forcing file from buoy data via wavebuoy2forceroms.m');

%% Dimensions:

disp(' ## Defining Dimensions...')

LP=xsize;
MP=ysize;
L=LP-1;
M=MP-1;
N=0;

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;

%nc('s_rho') = N;
%nc('s_w') = N+1;
nc('wave_time')=nt;
nc('one') = 1;

%% Variables and attributes:
disp(' ## Defining Dimensions, Variables, and Attributes...')

nc{'wave_time'} = ncdouble('wave_time');
nc{'wave_time'}.long_name = ncchar('wave field time');
nc{'wave_time'}.units = ncchar('days since 1968-05-23 00:00 UTC');
nc{'wave_time'}.field = ncchar('wave_time, scalar, series');
if(write_depth)
   nc{'depth'} = ncfloat('eta_rho', 'xi_rho');
   nc{'depth'}.long_name = ncchar('swan_depth');
   nc{'depth'}.units = ncchar('meter');
   nc{'depth'}.field = ncchar('depth, scalar, series');
end
if(write_dissip)
   nc{'Wave_dissip'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'Wave_dissip'}.long_name = ncchar('Wave_dissip');
   nc{'Wave_dissip'}.units = ncchar('Watts meter-2');
   nc{'Wave_dissip'}.field = ncchar('Wave_dissip, scalar, series');
end
if(write_force)
   nc{'Wave_forcex'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'Wave_forcex'}.long_name = ncchar('Wave_forcex');
   nc{'Wave_forcex'}.units = ncchar('Newtons meter-2');
   nc{'Wave_forcex'}.field = ncchar('Wave_forcex, scalar, series');

   nc{'Wave_forcey'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'Wave_forcey'}.long_name = ncchar('Wave_forcey');
   nc{'Wave_forcey'}.units = ncchar('Newtons meter-2');
   nc{'Wave_forcey'}.field = ncchar('Wave_forcey, scalar, series');
end
if(write_hsig)
   nc{'Hwave'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'Hwave'}.long_name = ncchar('Hwave');
   nc{'Hwave'}.units = ncchar('meter');
   nc{'Hwave'}.field = ncchar('Hwave, scalar, series');
end
if(write_setup),
   nc{'Wave_setup'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'Wave_setup'}.long_name = ncchar('Wave_setup');
   nc{'Wave_setup'}.units = ncchar('meter');
   nc{'Wave_setup'}.field = ncchar('Wave_setup, scalar, series');
end
if(write_rtp)
   nc{'Pwave_top'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'Pwave_top'}.long_name = ncchar('Pwave_top');
   nc{'Pwave_top'}.units = ncchar('second');
   nc{'Pwave_top'}.field = ncchar('Pwave_top, scalar, series');
end
if(write_tmbot)
   nc{'Pwave_bot'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'Pwave_bot'}.long_name = ncchar('Pwave_bot');
   nc{'Pwave_bot'}.units = ncchar('second');
   nc{'Pwave_bot'}.field = ncchar('Pwave_bot, scalar, series');
end
if(write_ubot)
   nc{'Ub_swan'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'Ub_swan'}.long_name = ncchar('Ub_swan');
   nc{'Ub_swan'}.units = ncchar('meter second-1');
   nc{'Ub_swan'}.field = ncchar('Ub_swan, scalar, series');
end
if(write_vel)
   nc{'velx'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'velx'}.long_name = ncchar('velx');
   nc{'velx'}.units = ncchar('meter second-1');
   nc{'velx'}.field = ncchar('velx, scalar, series');

   nc{'vely'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'vely'}.long_name = ncchar('vely');
   nc{'vely'}.units = ncchar('meter second-1');
   nc{'vely'}.field = ncchar('vely, scalar, series');
end
if(write_wdir)
   nc{'Dwave'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'Dwave'}.long_name = ncchar('Dwave');
   nc{'Dwave'}.units = ncchar('degrees');
   nc{'Dwave'}.field = ncchar('Dwave, scalar, series');
end
if(write_wlen)
   nc{'Lwave'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'Lwave'}.long_name = ncchar('Lwave');
   nc{'Lwave'}.units = ncchar('meter');
   nc{'Lwave'}.field = ncchar('Lwave, scalar, series');
end
if(write_xp)
   nc{'lon'} = ncfloat('eta_rho', 'xi_rho');
   nc{'lon'}.long_name = ncchar('lon');
   nc{'lon'}.units = ncchar('degree_east');
   nc{'lon'}.FillValue_ = ncfloat(100000.);
   nc{'lon'}.missing_value = ncfloat(100000.);
end
if(write_yp)
   nc{'lat'} = ncfloat('eta_rho', 'xi_rho');
   nc{'lat'}.long_name = ncchar('yp');
   nc{'lat'}.units = ncchar('degree_north');
   nc{'lat'}.missing_value = ncfloat(100000.);
   nc{'lat'}.FillValue_ = ncfloat(100000.);
end
if(write_wind)
   nc{'Windx'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'Windx'}.long_name = ncchar('Windx');
   nc{'Windx'}.units = ncchar('meter second-1');
   nc{'Windx'}.field = ncchar('Windx, scalar, series');

   nc{'Windy'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
   nc{'Windy'}.long_name = ncchar('Windy');
   nc{'Windy'}.units = ncchar('meter second-1');
   nc{'Windy'}.field = ncchar('Windy, scalar, series');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%now write the data from the arrays to the netcdf file

if(write_hsig)
% already read in above
disp('hsig -> Hwave')
nc{'Hwave'}(:)  = Hsig;
end
nc{'depth'}(:) = depth;
if(write_dissip)
   nc{'Wave_dissip'}(:) = dissip;
end

if(write_force)
   disp('force => Wave_forcex, Wave_forcey')
   nc{'Wave_forcex'}(:) = force(:,:,:,1);
   nc{'Wave_forcey'}(:) = force(:,:,:,2);
end
if(write_setup)
   disp('setup -> Wave_setup')
   nc{'Wave_setup'}(:) = setup;
end
if(write_rtp)
    disp(' - > Pwave_top')
   nc{'Pwave_top'}(:) = Tdom;
end
if(write_tmbot)
   disp('tmbot - > Pwave_bot')
   nc{'Pwave_bot'}(:) = Tmbot;
   clear tmbot
end
if(write_ubot)
   disp('ubot -> Ub_swan')
   nc{'Ub_swan'}(:) = Ubot;
end
if(write_vel)
   disp('-> velx, vely')
   nc{'velx'}(:) = vel(:,:,:,1);
   nc{'vely'}(:) = vel(:,:,:,2);
end
if(write_wdir)
   disp('-> Dwave')
   nc{'Dwave'}(:) = Dwave;
end
if(write_wlen)
   disp('-> Lwave')
   nc{'Lwave'}(:) = Lwave;
end
if(write_xp)
   disp('-> rlon')
   nc{'lon'}(:) = rlon;
end
if(write_yp)
   disp('-> rlat')
   nc{'lat'}(:) = rlat;
end
if(write_wind)
   disp('wind -> Windx, Windy')
   nc{'Windx'}(:) = wind(:,:,:,1);
   nc{'Windy'}(:) = wind(:,:,:,2);
end

% write time
nc{'wave_time'}(:) = wave_jd-2440000;   % Modified Julian Days = (Julian Day - 2440000.) = days since 1968-05-23 00:00 UTC
close(nc)

%% Plot the values, if this was requested - in the plots, defaults are replaced with NaN's.
if exist('plot_key')
if plot_key >= 1;

   if nt == 1
      ncload(fname)
   else
      tidx =plot_key;
      nc=netcdf(fname);
      depth=squeeze(nc{'depth'}(:,:));
      Wave_dissip=squeeze(nc{'Wave_dissip'}(tidx,:,:));
      Wave_forcex=squeeze(nc{'Wave_forcex'}(tidx,:,:));
      Wave_forcey=squeeze(nc{'Wave_forcey'}(tidx,:,:));
      Hwave=squeeze(nc{'Hwave'}(tidx,:,:));
      Wave_setup=squeeze(nc{'Wave_setup'}(tidx,:,:));
      Pwave_top=squeeze(nc{'Pwave_top'}(tidx,:,:));
      Pwave_bot=squeeze(nc{'Pwave_bot'}(tidx,:,:));
      Ub_swan=squeeze(nc{'Ub_swan'}(tidx,:,:));
      velx=squeeze(nc{'velx'}(tidx,:,:));
      vely=squeeze(nc{'vely'}(tidx,:,:));
      Dwave=squeeze(nc{'Dwave'}(tidx,:,:));
      Lwave=squeeze(nc{'Lwave'}(tidx,:,:));
      windx=squeeze(nc{'Windx'}(tidx,:,:));
      windy=squeeze(nc{'Windy'}(tidx,:,:));
      xp=nc{'lon'}(:);
      yp=nc{'lat'}(:);
   end
   figure;
   swan_depth(find(depth == 0)) = NaN;
   pcolor(xp,yp,squeeze(depth(:,:))); shading interp; colorbar;
   title('Depth Output of Swan');

   figure;
   %   hsig(find(hsig == 0)) = NaN;
   pcolor(xp,yp,squeeze(Hwave(:,:))); shading interp; colorbar;
   title('Hs Output of Swan');

   figure;
   %    Ub_swan(find(Ub_swan == 0.001)) = NaN;
   pcolor(xp,yp,squeeze(Ub_swan(:,:))); shading interp; colorbar;
   title('Ub Output of Swan');

   figure;
   Pwave_bot(find(Pwave_bot == 9999)) = NaN;
   pcolor(xp,yp,squeeze(Pwave_bot(:,:))); shading interp;
   title('Pwave_bot Output of Swan');
   caxis([2 15]);colorbar
   
   if write_wdir
   figure;
   %    Dwave(find(Dwave == 0)) = NaN;
   pcolor(xp,yp,squeeze(Dwave(:,:))); shading interp; colorbar;
   title('Dwave Output of Swan');
   end

   if write_setup 
   figure;
   %    Dwave(find(Dwave == 0)) = NaN;
   pcolor(xp,yp,squeeze(Wave_setup(:,:))); shading interp; colorbar;
   title('setup Output of Swan');
   end

   if write_wlen 
   figure;
   %    Dwave(find(Dwave == 0)) = NaN;
   pcolor(xp,yp,squeeze(Lwave(:,:))); shading interp;
   title('Lwave Output of Swan');
   caxis([0 75]);colorbar
   end
   
   if write_dissip
   figure
   %    Dwave(find(Dwave == 0)) = NaN;
   pcolor(xp,yp,squeeze(Wave_dissip(:,:))); shading interp; colorbar;
   end

   if write_force 
   figure
   subplot(121)
   %    Dwave(find(Dwave == 0)) = NaN;
   pcolor(xp,yp,squeeze(Wave_forcex(:,:))); shading interp; colorbar;
   title('forcex Output of Swan');

   subplot(122)
   %    Dwave(find(Dwave == 0)) = NaN;
   pcolor(xp,yp,squeeze(Wave_forcey(:,:))); shading interp; colorbar;
   title('forcey Output of Swan');
   end

   if write_force 
   figure
   subplot(121)
%    vel(find(vel == 0)) = NaN;
   pcolor(xp,yp,squeeze(velx(:,:))); shading interp; colorbar;
   title('velx Output of Swan');
   
   subplot(122)
   %    vel(find(vel == 0)) = NaN;
   pcolor(xp,yp,squeeze(vely(:,:))); shading interp; colorbar;
   title('vely Output of Swan');
   end
   if write_wind
   figure
   subplot(121)
%    vel(find(vel == 0)) = NaN;
   pcolor(xp,yp,squeeze(windx(:,:))); shading interp; colorbar;
   title('wind x Output of Swan');
   
   subplot(122)
   %    vel(find(vel == 0)) = NaN;
   pcolor(xp,yp,squeeze(windy(:,:))); shading interp; colorbar;
   title('wind y Output of Swan');
   end
end
end