function SWAN2ROMS(iname,fname,plot_key);
% SWAN2ROMS - Convert SWAN output files to ROMS input netCDF files
% Usage: SWAN2ROMS(iname,fname,plot_key);
% 
% Inputs: iname = SWAN input file name (e.g. 'INPUT')
%         fname = ROMS forcing file name (.e.g. 'swan_frc.nc')
%         plot_key = 1 to create plots, 0 otherwise

% Function to convert SWAN output data files of interest to ROMS (depth.dat, hsig.dat,
% rtp.dat, ubot.dat, wdir.dat, xp.dat, and yp.dat) to a netCDF file ('fname').
% If plot_key = 0, the program does not plot the variables; if plot_key = 1, the first
% time of each output variable is plotted in a matlab figure window.

% SWAN2ROMS expects to read 'INPUT' in Layout 4 format.  Snippet:

% ! Model specification
% CGRID CURVILINEAR 159 59 EXC 99999. CIRCLE 36 0.05 0.6 26
% ...
% ! Model output
% GROUP 'group1' SUBG 0 159 0 59
% BLOCK 'group1' NOHEAD 'xp.dat' LAY 4 XP
% BLOCK 'group1' NOHEAD 'yp.dat' LAY 4 YP
% BLOCK 'group1' NOHEAD 'depth.dat'  LAY 4 DEPTH
% BLOCK 'group1' NOHEAD 'hsig.dat'   LAY 4 HSIG   OUTPUT  20020912.000000 3 HR
% BLOCK 'group1' NOHEAD 'wdir.dat'   LAY 4 DIR    OUTPUT  20020912.000000 3 HR
% BLOCK 'group1' NOHEAD 'rtp.dat'    LAY 4 RTP    OUTPUT  20020912.000000 3 HR
% BLOCK 'group1' NOHEAD 'tmbot.dat'  LAY 4 TMBOT  OUTPUT  20020912.000000 3 HR
% BLOCK 'group1' NOHEAD 'ubot.dat'   LAY 4 UBOT   OUTPUT  20020912.000000 3 HR
% BLOCK 'group1' NOHEAD 'wind.dat'   LAY 4 WIND   OUTPUT  20020912.000000 3 HR
% BLOCK 'group1' NOHEAD 'wlen.dat'   LAY 4 WLEN   OUTPUT  20020912.000000 3 HR

% ! Wave-current interaction output
% BLOCK 'group1' NOHEAD 'dissip.dat' LAY 4 DISSIP OUTPUT  20020912.000000 3 HR
% BLOCK 'group1' NOHEAD 'force.dat'  LAY 4 FORCE  OUTPUT  20020912.000000 3 HR
% BLOCK 'group1' NOHEAD 'vel.dat'    LAY 4 VEL    OUTPUT  20020912.000000 3 HR
% BLOCK 'group1' NOHEAD 'setup.dat'  LAY 4 SETUP  OUTPUT  20020912.000000 3 HR

% ! Wave-current interaction output
% BLOCK 'COMPGRID' NOHEADER 'dissip.dat' LAY 4 DISSIP 1.
% BLOCK 'COMPGRID' NOHEADER 'force.dat' LAY 4 FORCE 1.
% BLOCK 'COMPGRID' NOHEADER 'setup.dat' LAY 4 SETUP 1.
% BLOCK 'COMPGRID' NOHEADER 'vel.dat' LAY 4 VEL 1.

%Soupy Alexander, 6/13/02
%jcwarner: added internal netcdf creation
%csherwood@usgs.gov: changed documentation a little, and made
%      wave-current interaction terms optional
% rsignell@usgs.gov: Rich Signell : rewrote to use "textread" instead of "fgetl" (25x faster)
%                                 : writes correct time vector in Modified Julian Day 
%                                 : based on info from INPUT 
%
%jcwarner 30Sep2005 - fix a few small plotting bug fixes in var names. Calc nt only once.


% 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;

if(exist('wind.dat','file')==2),write_wind=1;,end
if(exist('dissip.dat','file')==2),write_dissip=1;,end
if(exist('force.dat','file')==2),write_force=1;,end
if(exist('setup.dat','file')==2),write_setup=1;,end
if(exist('vel.dat','file')==2),write_vel=1;,end
if(exist('wlen.dat','file')==2),write_wlen=1;,end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read meta data from SWAN input file
s=swan_meta(iname)
xsize=s.xsize;
ysize=s.ysize;
jdmat0=s.jdmat0;
dt=s.dt;  % seconds
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%read in hsig because we need to find out how many time steps there are.
disp('hsig')
hsig=textread('hsig.dat','%f');
nt=length(hsig)/(xsize*ysize);
hsig=reshape(hsig,xsize,ysize,nt);
hsig=permute(hsig,[3 2 1]);
hsig(find(hsig == -9)) = 0.0;

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

%% Global attributes:

disp(' ## Defining Global Attributes...')
nc.history = ncchar(['Created by ROMS2SWAN on ' datestr(now)]);
nc.type = ncchar('forcing file from swan output via swan2roms.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'); %% 1 element.
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'}(:)  = squeeze(hsig);
clear hsig
end

if(write_depth)
   disp('depth')
   depth=textread('depth.dat','%f');
   depth=reshape(depth,xsize,ysize);
   depth=permute(depth,[2 1]);
   depth(find(depth == -99)) = 0.0;
   nc{'depth'}(:) = depth;
   clear depth
end

if(write_dissip)
   disp('dissip -> Wave_dissip')
   dissip=textread('dissip.dat','%f');
   dissip=reshape(dissip,xsize,ysize,nt);
   dissip=permute(dissip,[3 2 1]);
   dissip(find(dissip == -9)) = 0.0;
   nc{'Wave_dissip'}(:) = dissip;
   clear dissip
end

if(write_force)
   disp('force => Wave_forcex, Wave_forcey')
   force=textread('force.dat','%f');
   force=reshape(force,xsize,ysize,2,nt);
   force=permute(force,[4 2 1 3]);
   force(find(force == 0)) = 0.0;
   nc{'Wave_forcex'}(:) = force(:,:,:,1);
   nc{'Wave_forcey'}(:) = force(:,:,:,2);
   clear force 
end
if(write_setup)
   disp('setup -> Wave_setup')
   setup=textread('setup.dat','%f');
   setup=reshape(setup,xsize,ysize,nt);
   setup=permute(setup,[3 2 1]);
   setup(find(setup == -9)) = 0.0;
   nc{'Wave_setup'}(:) = setup;
   clear setup
end
if(write_rtp)
   disp('rtp - > Pwave_top')
   rtp=textread('rtp.dat','%f');
   rtp=reshape(rtp,xsize,ysize,nt);
   rtp=permute(rtp,[3 2 1]);
   rtp(find(rtp == -9)) = 10.0;
   nc{'Pwave_top'}(:) = rtp;
   clear rtp
end
if(write_tmbot)
   disp('tmbot - > Pwave_bot')
   tmbot=textread('tmbot.dat','%f');
   tmbot=reshape(tmbot,xsize,ysize,nt);
   tmbot=permute(tmbot,[3 2 1]);
   tmbot(find(tmbot == -9)) = 10.0;
   nc{'Pwave_bot'}(:) = tmbot;
   clear tmbot
end

if(write_ubot)
   disp('ubot -> Ub_swan')
   ubot=textread('ubot.dat','%f');
   ubot=reshape(ubot,xsize,ysize,nt);
   ubot=permute(ubot,[3 2 1]);
   ubot(find(ubot == -10)) = 0.0001;
   nc{'Ub_swan'}(:) = ubot;
   clear ubot
end

if(write_vel)
   disp('vel -> velx, vely')
   vel=textread('vel.dat','%f');
%   nt=length(vel)/(xsize*ysize*2)
   vel=reshape(vel,xsize,ysize,2,nt);
   vel=permute(vel,[4 2 1 3]);
   vel(find(vel == 0)) = 0.0;
   nc{'velx'}(:) = vel(:,:,:,1);
   nc{'vely'}(:) = vel(:,:,:,2);
   clear vel
end
if(write_wdir)
   disp('wdir -> Dwave')
   wdir=textread('wdir.dat','%f');
   wdir=reshape(wdir,xsize,ysize,nt);
   wdir=permute(wdir,[3 2 1]);
   wdir(find(wdir == -999)) = 0.0;
   nc{'Dwave'}(:) = wdir;
   clear wdir
end

if(write_wlen)
   disp('wlen -> Lwave')
   wlen=textread('wlen.dat','%f');
   wlen=reshape(wlen,xsize,ysize,nt);
   wlen=permute(wlen,[3 2 1]);
   wlen(find(wlen == -9)) = 10.0;
   nc{'Lwave'}(:) = wlen;
   clear wlen 
end

if(write_xp)
   disp('lon')
   xp=textread('xp.dat','%f');
   xp=reshape(xp,xsize,ysize);
   xp=permute(xp,[2 1]);
   nc{'lon'}(:) = xp;
   clear xp
end

if(write_yp)
   disp('lat')
   yp=textread('yp.dat','%f');
   yp=reshape(yp,xsize,ysize);
   yp=permute(yp,[2 1]);
   nc{'lat'}(:) = yp;
   clear yp
end

if(write_wind)
   disp('wind -> Windx, Windy')
   wind=textread('wind.dat','%f');
   wind=reshape(wind,xsize,ysize,2,nt);
   wind=permute(wind,[4 2 1 3]);
   wind(find(wind == -9)) = 0.0;
   nc{'Windx'}(:) = wind(:,:,:,1);
   nc{'Windy'}(:) = wind(:,:,:,2);
   clear wind
end


jdshift=datenum(1968,5,23,0,0,0);
mjd_days=jdmat0+[0:(nt-1)]*(dt/(24*3600))-jdshift;  %time in days
% write time
nc{'wave_time'}(:) = mjd_days;   % 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