function swanmat2roms(input)
% swanmat2roms - Convert SWAN output mat files to ROMS input netCDF files
% Usage: swanmat2roms(swan);
%
% Inputs: input=INPUT file name i.e.
% 'Y:\barmstrong\data\Isabel\SWAN\run16\INPUT'

%swan_meta_v3 creates:
%           swan = structure containing:
%           swan.input = name of swan INPUT file
%           swan.var = cell array of variable names
%           swan.matfile = cell array of mat file names corresponding
%                                 to variable names
%   
%           Note: recognized variable names:
%           {'depth','dissip','hsig','qb','rtp','setup','tmbot',...
%            'ubot','vel','wdir','wind','wlen','xp','yp'}
%
% Outputs: A NetCDF file produced for each input variable
%
% Example:
%     swan.input='INPUT';
%     swan.var ={'depth',    'xp',    'yp',    'hsig',    'ubot'}
%     swan.matfile={'depth.mat','xp.mat','yp.mat','hsig.mat','ubot.mat'}
%     swanmat2roms(swan)
%       makes NetCDF files: depth.nc, xp.nc, yp.nc, hsig.nc and ubot.nc

% read meta data from SWAN input file
[swan]=swan_meta_v3(input);
s=swan;
%s=swan_meta_v2(swan.input)
xsize=s.xsize;
ysize=s.ysize;
jdmat0=s.jdmat0;
dt=s.dt;    % seconds
nt=s.nt;    % number of records each file

% Create netcdf files for each variable
for i=1%:length(swan.var);
    tic

    nc=netcdf([swan.var{i} '.nc'],'clobber');
    if isempty(nc), return, end


    %% Global attributes:

    disp(' ## Defining Global Attributes...')
    nc.history = ncchar(['Created by ROMS2SWAN_v3 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('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-5-23 00:00 UTC');
    nc{'wave_time'}.field = ncchar('wave_time, scalar, series');

    switch swan.var{i}
        case 'depth'
            nc{'swan_depth'} = ncfloat('wave_time','eta_rho', 'xi_rho');
            nc{'swan_depth'}.long_name = ncchar('swan_depth');
            nc{'swan_depth'}.units = ncchar('meter');
            nc{'swan_depth'}.coordinates = ncchar('xp yp');
            nc{'swan_depth'}.field = ncchar('depth, scalar, series');
            disp('   depth -> swan_depth')
            F=load('depth.mat');
            vname=fieldnames(F);
            Depth=[];
            for K=1:length(vname)
                Depth=getfield(F,char(vname(K)));
                Depth(Depth == -99) = 0.0;     % flag
                Depth(isnan(Depth))=0.0;       % land
                nc{'swan_depth'}(K,:,:)=squeeze(Depth);
                disp(sprintf('%g/%g swan_depth saved',K,length(vname)))
                clear Depth;
            end
            clear F vname K
        case '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'}.coordinates = ncchar('xp yp');
            nc{'Wave_dissip'}.field = ncchar('Wave_dissip, scalar, series');
            disp('   dissip -> Wave_dissip')
            F=load('dissip.mat');
            vname=fieldnames(F);
            Dissip=[];
            for K=1:length(vname)
                Dissip=getfield(F,char(vname(K)));
                Dissip(find(Dissip == -9)) = 0.0;    % flag
                Dissp(find(isnan(Dissip))) = 0.0;    % land
                nc{'Wave_dissip'}(K,:,:) = squeeze(Dissip);
                disp(sprintf('%g/%g Dissip saved',K,length(vname)))
                clear Dissip;
            end
            clear F vname K
        case 'force'
            nc{'Wave_forcex'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
            nc{'Wave_forcex'}.long_name = ncchar('Wave_forcex');
            nc{'Wave_forcex'}.coordinates = ncchar('xp yp');
            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'}.coordinates = ncchar('xp yp');
            nc{'Wave_forcey'}.units = ncchar('Newtons meter-2');
            nc{'Wave_forcey'}.field = ncchar('Wave_forcey, scalar, series');
            disp('   force => Wave_forcex, Wave_forcey')
            F=load('force.mat','WForce_x*');
            vname=fieldnames(F);
            for K=1:length(vname);
                Forcex=getfield(F,char(vname(K)));
                Forcex(Forcex == 0) = 0.0;     % flag
                Forcex(isnan(Forcex)) = 0.0;   % land
                nc{'Wave_forcex'}(K,:,:) = squeeze(Forcex);
                disp(sprintf('%g/%g WForce_x saved',K,length(vname)))
                clear Force*;
            end
            clear F vname K
            F=load('force.mat','WForce_y*');
            vname=fieldnames(F);
            for K=1:length(vname);
                Forcey=getfield(F,char(vname(K)));
                Forcey(Forcey == 0) = 0.0;     % flag
                Forcey(isnan(Forcey)) = 0.0;   % land
                nc{'Wave_forcey'}(K,:,:) = squeeze(Forcey);
                disp(sprintf('%g/%g WForce_y saved',K,length(vname)))
                clear Force*;
            end
            clear F vname K
        case 'hsig'
            nc{'Hwave'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
            nc{'Hwave'}.long_name = ncchar('Hwave');
            nc{'Hwave'}.units = ncchar('meter');
            nc{'Hwave'}.coordinates = ncchar('xp yp');
            nc{'Hwave'}.field = ncchar('Hwave, scalar, series');
            disp('   hsig -> Hwave')
            F=load('hsig.mat');
            vname=fieldnames(F);
            for K=1:length(vname)
                Hsig=getfield(F,char(vname(K)));
                Hsig(find(Hsig == -99)) = 0.0;   % flag
                Hsig(find(isnan(Hsig)))=0.0;     % land
                nc{'Hwave'}(K,:,:)=squeeze(Hsig);
                disp(sprintf('%g/%g Hwave saved',K,length(vname)))
                clear Hsig;
            end
            clear F vname K
        case 'qb'
            nc{'Wave_break'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
            nc{'Wave_break'}.long_name = ncchar('Wave_break');
            nc{'Wave_break'}.units = ncchar('meter');
            nc{'Wave_break'}.field = ncchar('Wave_break, scalar, series');
            disp('   qb -> Wave break')
            F=load('qb.mat');
            vname=fieldnames(F);
            for K=1:length(vname)
                Qb=getfield(F,char(vname(K)));
                Qb(find(Qb == -9)) = 10.0;   % flag
                Qb(find(isnan(Qb))) = 0.0;   % land
                nc{'Wave_break'}(K,:,:) = squeeze(Qb);
                disp(sprintf('%g/%g ''Qb'' saved',K,length(vname)))
                clear Qb;
            end
            clear F vname K
        case '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'}.coordinates = ncchar('xp yp');
            nc{'Pwave_top'}.field = ncchar('Pwave_top, scalar, series');
            disp('   rtp - > Pwave_top')
            F=load('rtp.mat');
            vname=fieldnames(F);
            for K=1:length(vname)
                Rtp=getfield(F,char(vname(K)));
                Rtp(find(Rtp == -9)) = 10.0;     % flag
                Rtp(find(isnan(Rtp))) = 0.0;     % land
                nc{'Pwave_top'}(K,:,:) = squeeze(Rtp);
                disp(sprintf('%g/%g ''Pwave_top'' saved',K,length(vname)))
                clear Rtp;
            end
            clear F vname K
        case 'setup'
            nc{'Wave_setup'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
            nc{'Wave_setup'}.long_name = ncchar('Wave_setup');
            nc{'Wave_setup'}.coordinates = ncchar('xp yp');
            nc{'Wave_setup'}.units = ncchar('meter');
            nc{'Wave_setup'}.field = ncchar('Wave_setup, scalar, series');
            disp('   setup -> Wave_setup')
            F=load('setup.mat');
            vname=fieldnames(F);
            Setup=[];
            for K=1:length(vname)
                Setup=getfield(F,char(vname(K)));
                Setup(find(Setup == -9)) = 0.0;      % flag
                Setup(find(isnan(Setup))) = 0.0;     % land
                nc{'Wave_setup'}(K,:,:) = squeeze(Setup);
                disp(sprintf('%g/%g ''Setup'' saved',K,length(vname)))
                clear Setup;
            end
            clear F vname K
        case 'tmbot'
            nc{'Pwave_bot'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
            nc{'Pwave_bot'}.long_name = ncchar('Pwave_bot');
            nc{'Pwave_bot'}.coordinates = ncchar('xp yp');
            nc{'Pwave_bot'}.units = ncchar('second');
            nc{'Pwave_bot'}.field = ncchar('Pwave_bot, scalar, series');
            disp('   tmbot - > Pwave_bot')
            F=load('tmbot.mat');
            vname=fieldnames(F);
            for K=1:length(vname)
                Tmbot=getfield(F,char(vname(K)));
                Tmbot(find(Tmbot == -9)) = 10.0;     % flag
                Tmbot(find(isnan(Tmbot))) = 0.0;     % land
                nc{'Pwave_bot'}(K,:,:) = squeeze(Tmbot);
                disp(sprintf('%g/%g ''Pwave_bottom'' saved',K,length(vname)))
                clear Tmbot;
            end
            clear F vname K

        case 'ubot'
            nc{'Ub_swan'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
            nc{'Ub_swan'}.long_name = ncchar('Ub_swan');
            nc{'Ub_swan'}.coordinates = ncchar('xp yp');
            nc{'Ub_swan'}.units = ncchar('meter second-1');
            nc{'Ub_swan'}.field = ncchar('Ub_swan, scalar, series');
            disp('   ubot -> Ub_swan')
            F=load('ubot.mat');
            vname=fieldnames(F);
            for K=1:length(vname)
                Ubot=getfield(F,char(vname(K)));
                Ubot(find(Ubot == -10)) = 0.0001;    % flag
                Ubot(find(isnan(Ubot))) = 0.0;       % land
                nc{'Ub_swan'}(K,:,:) = squeeze(Ubot);
                disp(sprintf('%g/%g ''Ub_swan'' saved',K,length(vname)))
                clear Ubot;
            end
            clear F vname K
        case 'vel'
            nc{'velx'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
            nc{'velx'}.long_name = ncchar('velx');
            nc{'velx'}.coordinates = ncchar('xp yp');
            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'}.coordinates = ncchar('xp yp');
            nc{'vely'}.units = ncchar('meter second-1');
            nc{'vely'}.field = ncchar('vely, scalar, series');
            disp('   vel -> velx, vely')


            F=load('vel.mat');
            vname=fieldnames(F);
            for K=1:2:length(vname)-1
                Velx=getfield(F,char(vname(K)));
                Vely=getfield(F,char(vname(K+1)));
                Velx(find(Velx == 0)) = 0.0;         % flag
                Vely(find(Vely == 0)) = 0.0;         % flag
                Velx(find(isnan(Velx))) = 0.0;       % land
                Vely(find(isnan(Vely))) = 0.0;       % land
                nc{'velx'}(K,:,:) = squeeze(Velx);
                nc{'vely'}(K,:,:) = squeeze(Vely);
                disp(sprintf('%g/%g ''Velocity'' saved',K,length(vname)))
                clear Vel*;
            end
            clear F vname K
        case 'wdir'
            nc{'Dwave'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
            nc{'Dwave'}.long_name = ncchar('Dwave');
            nc{'Dwave'}.coordinates = ncchar('xp yp');
            nc{'Dwave'}.units = ncchar('degrees');
            nc{'Dwave'}.field = ncchar('Dwave, scalar, series');
            disp('   wdir -> Dwave')
            F=load('wdir.mat');
            vname=fieldnames(F);
            for K=1:length(vname)
                Wdir=getfield(F,char(vname(K)));
                Wdir(find(Wdir == -999)) = 0.0;  % flag
                Wdir(find(isnan(Wdir))) = 0.0;   % land
                nc{'Dwave'}(K,:,:) = squeeze(Wdir);
                disp(sprintf('%g/%g ''Dwave'' saved',K,length(vname)))
                clear Wdir;
            end
            clear F vname K
        case 'wind'
            nc{'Windx'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
            nc{'Windx'}.long_name = ncchar('Windx');
            nc{'Windx'}.coordinates = ncchar('xp yp');
            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'}.coordinates = ncchar('xp yp');
            nc{'Windy'}.units = ncchar('meter second-1');
            nc{'Windy'}.field = ncchar('Windy, scalar, series');

            disp('   wind -> Windx, Windy')
            F=load('wind.mat','Windv_x*');
            vname=fieldnames(F);
            for K=1:length(vname);
                Windx=getfield(F,char(vname(K)));
                Windx(Windx == -9) = 0.0;      % flag
                Windx(isnan(Windx)) = 0.0;     % land
                nc{'Windx'}(K,:,:) = squeeze(Windx);
                disp(sprintf('%g/%g ''Wind_x'' saved',K,length(vname)))
                clear Wind*;
            end
            clear F vname K
            F=load('wind.mat','Windv_y*');
            vname=fieldnames(F);
            for K=1:length(vname);
                Windy=getfield(F,char(vname(K)));
                Windy(Windy == -9) = 0.0;      % flag
                Windy(isnan(Windy)) = 0.0;     % land
                nc{'Windy'}(K,:,:) = squeeze(Windy);
                disp(sprintf('%g/%g ''Wind_y'' saved',K,length(vname)))
                clear Wind*;
            end
            clear F vname K
        case 'wlen'
            nc{'Lwave'} = ncfloat('wave_time', 'eta_rho', 'xi_rho');
            nc{'Lwave'}.long_name = ncchar('Lwave');
            nc{'Lwave'}.units = ncchar('meter');
            nc{'Lwave'}.coordinates = ncchar('xp yp');
            nc{'Lwave'}.field = ncchar('Lwave, scalar, series');
            disp('   wlen -> Lwave')
            F=load('wlen.mat');
            vname=fieldnames(F);
            for K=1:length(vname)
                Wlen=getfield(F,char(vname(K)));
                Wlen(find(Wlen == -9)) = 10.0;   % flag
                Wlen(find(isnan(Wlen))) = 0.0;   % land
                nc{'Lwave'}(K,:,:) = squeeze(Wlen);
                disp(sprintf('%g/%g ''Lwave'' saved',K,length(vname)))
                clear Wlen;
            end
            clear F vname K

        case 'xp'
            nc{'xp'} = ncfloat('eta_rho', 'xi_rho');
            nc{'xp'}.long_name = ncchar('xp');
            nc{'xp'}.units = ncchar('degrees_east');
            nc{'xp'}.FillValue_ = ncfloat(100000.);
            nc{'xp'}.missing_value = ncfloat(100000.);
            nc{'xp'}.field = ncchar('xp, scalar, series');

            disp('   xp')
            F=load('xp.mat');
            vname=fieldnames(F);
            for K=1:1 %length(vname)
                Xp=getfield(F,char(vname(K)));
                nc{'xp'}(:,:) = squeeze(Xp);
                disp(sprintf('%g/%g ''xp'' saved',K,length(vname)))
                clear Xp;
            end
            clear F vname K
        case 'yp'
            nc{'yp'} = ncfloat('eta_rho', 'xi_rho');
            nc{'yp'}.long_name = ncchar('yp');
            nc{'yp'}.units = ncchar('degrees_north');
            nc{'yp'}.missing_value = ncfloat(100000.);
            nc{'yp'}.FillValue_ = ncfloat(100000.);
            nc{'yp'}.field = ncchar('yp, scalar, series');
            disp('   yp')
            F=load('yp.mat');
            vname=fieldnames(F);
            for K=1:1 %length(vname)
                Yp=getfield(F,char(vname(K)));
                nc{'yp'}(:,:) = squeeze(Yp);
                disp(sprintf('%g/%g ''yp'' saved',K,length(vname)))
                clear Yp;
            end
            clear F vname K

        otherwise
            disp(['Unrecognized variable ' swan.var{i} ' in swan.var']);
    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)

    toc


end