%%% Clear screen and ignore warnings

fclose all;
clc;
warning off all;


%%% STANDARD GUI CHECK ON DIRECTORIES

convertGuiDirectoriesCheck;


%%% GUI CHECKS

% Check if the mdf file name has been specified (Delft3D)
filemdf     = get(handles.edit3,'String');
filemdf     = deblank2(filemdf);
if isempty(filemdf);
    if exist('wb'); close(wb); end;
    errordlg('The mdf file name has not been specified.','Error');
    break;
end

% Check if the mdu file name has been specified (D-Flow FM)
mdufile     = get(handles.edit4,'String');
if isempty(mdufile);
    if exist('wb'); close(wb); end;
    errordlg('The mdu file name has not been specified.','Error');
    break;
end

% Check if the net file name has been specified (D-Flow FM)
netfile     = get(handles.edit8,'String');
if isempty(netfile);
    if exist('wb'); close(wb); end;
    errordlg('The net file name has not been specified.','Error');
    break;
end

% Check if the ext file name has been specified (D-Flow FM)
extfile     = get(handles.edit10,'String');
if isempty(extfile);
    if exist('wb'); close(wb); end;
    errordlg('The external forcings file name has not been specified.','Error');
    break;
end

% Check if the obs file name has been specified (D-Flow FM)
obsfile     = get(handles.edit22,'String');
if isempty(obsfile);
    obsfile = [];
end

% Check if the crs file name has been specified (D-Flow FM)
crsfile     = get(handles.edit31,'String');
if isempty(crsfile);
    crsfile = [];
end

% Check if the roughness file name has been specified (D-Flow FM)
rghfile     = get(handles.edit24,'String');
jargh       = 1;
if isempty(rghfile);
    rghfile = [];
    jargh   = 0;
end

% Check if the viscosity file name has been specified (D-Flow FM)
visfile     = get(handles.edit25,'String');
javis       = 1;
if isempty(visfile);
    visfile = [];
    javis   = 0;
end

% Check if the thin dam file name has been specified (D-Flow FM)
thdfile     = get(handles.edit27,'String');
jathd       = 1;
if isempty(thdfile);
    thdfile = [];
    jathd   = 0;
end

% Put the output directory name in the filenames
filemdf     = [pathin ,'\',filemdf];
mdufile     = [pathout,'\',mdufile];


%%% WRITE THE MDU FILE

% Read the mdf file
mdfcontents = delft3d_io_mdf('read',filemdf);
mdfkeywds   = mdfcontents.keywords;

% Assign some variables to be written in mdu-file
itdate      = mdfkeywds.itdate;
itdate(itdate=='-') = [];

% Write the mdu file
fid2        = fopen(mdufile,'wt'); 
klok        = clock;
fprintf(fid2,'%s\n'  ,['# Generated by the Matlab program dflowfmConverter on ',num2str(      klok(4) ,'%0.2d'),':' ,...
                                                                                num2str(      klok(5) ,'%0.2d'),':' ,...
                                                                                num2str(floor(klok(6)),'%0.2d'),', ',...
                                                                                num2str(      klok(3) ,'%0.2d'),'-' ,...
                                                                                num2str(      klok(2) ,'%0.2d'),'-' ,...
                                                                                num2str(      klok(1) ,'%0.4d')        ]);
fprintf(fid2,'%s\n\n',['# Deltares, Matlab-generated mdu-file for D-Flow FM']);  

fprintf(fid2,'%s\n'  ,'[model]');
fprintf(fid2,'%s\n'  ,['Program                             = D-Flow FM']);
fprintf(fid2,'%s\n'  ,['Version                             = 1.1.81']);
fprintf(fid2,'%s\n\n',['AutoStart                           = 0                                  # Autostart simulation after loading MDU or not (0=no, 1=autostart, 2=autostartstop).']);

fprintf(fid2,'%s\n'  ,'[geometry]');
fprintf(fid2,'%s\n'  ,['NetFile                             = ',netfile,'                        # *_net.nc']);
fprintf(fid2,'%s\n'  ,['BathymetryFile                      =                                    # *.xyb']);
fprintf(fid2,'%s\n'  ,['WaterLevIniFile                     =                                    # Initial water levels sample file *.xyz']);
fprintf(fid2,'%s\n'  ,['LandBoundaryFile                    =                                    # Only for plotting']);
if jathd == 1;
    fprintf(fid2,'%s\n'  ,['ThinDamFile                         = ',thdfile,'                       # *_thd.pli, Polyline(s) for tracing thin dams.']);
else
    fprintf(fid2,'%s\n'  ,['ThinDamFile                         =                                    # *_thd.pli, Polyline(s) for tracing thin dams.']);
end
fprintf(fid2,'%s\n'  ,['ThindykeFile                        =                                    # *._tdk.pli, Polyline(s) x,y,z, z = thin dyke top levels']);
fprintf(fid2,'%s\n'  ,['ProflocFile                         =                                    # *_proflocation.xyz)    x,y,z, z = profile refnumber']);
fprintf(fid2,'%s\n'  ,['ProfdefFile                         =                                    # *_profdefinition.def) definition for all profile nrs']);
fprintf(fid2,'%s\n'  ,['ManholeFile                         =                                    # *...']);
if isfield(mdfkeywds,'zeta0' ) == 0; 
    mdfkeywds.zeta0  =  0.0;
end
jasal = 1;
if isfield(mdfkeywds,'s0' ) == 0;
    if isfield(mdfkeywds,'salw') == 0;
        mdfkeywds.s0  =  0;
        jasal         =  0;
    else
        mdfkeywds.s0  =  mdfkeywds.salw;
        jasal         =  0;
    end
end
if isfield(mdfkeywds,'depuni') == 0;
    mdfkeywds.depuni = -5.0;
end
fprintf(fid2,'%s\n'  ,['WaterLevIni                         = ',num2str(mdfkeywds.zeta0,'%5.7f')                 ,'                          # Initial water level']);
fprintf(fid2,'%s\n'  ,['Bedlevuni                           = ',num2str(mdfkeywds.depuni,'%5.7f')                ,'                         # Uniform bottom level, (only if bedlevtype>=3, used at missing z values in netfile']);
fprintf(fid2,'%s\n'  ,['Bedslope                            =                                    # bedslopeination, sets zk = bedlevuni + x*bedslope ans sets zbndz = xbndz*bedslope']);
fprintf(fid2,'%s\n'  ,['BedLevType                          = 3                                  # 1 : Bottom levels at waterlevel cells (=flow nodes), like tiles xz, yz, bl , bob = max(bl left, bl right)']);
fprintf(fid2,'%s\n'  ,['                                                                         # 2 : Bottom levels at velocity points  (=flow links),            xu, yu, blu, bob = blu,    bl = lowest connected link']);
fprintf(fid2,'%s\n'  ,['                                                                         # 3 : Bottom levels at velocity points  (=flow links), using mean network levels xk, yk, zk  bl = lowest connected link']);
fprintf(fid2,'%s\n'  ,['                                                                         # 4 : Bottom levels at velocity points  (=flow links), using min  network levels xk, yk, zk  bl = lowest connected link']);
fprintf(fid2,'%s\n'  ,['                                                                         # 5 : Bottom levels at velocity points  (=flow links), using max  network levels xk, yk, zk  bl = lowest connected link']);
if isfield(mdfkeywds,'anglat') == 0; 
    mdfkeywds.anglat = 0;
end
fprintf(fid2,'%s\n'  ,['AngLat                              = ',num2str(mdfkeywds.anglat,'%2.2f')                ,'                              # Angle of latitude (deg), 0=no Coriolis']);
fprintf(fid2,'%s\n'  ,['Conveyance2D                        = 3                                  # -1:R=HU,0:R=H, 1:R=A/P, 2:K=analytic-1D conv, 3:K=analytic-2D conv']);
fprintf(fid2,'%s\n'  ,['Sillheightmin                       =                                    # Weir treatment only if both sills larger than this value (m) ']);
fprintf(fid2,'%s\n'  ,['Kmx                                 =                                    # Max nr of vertical layers']);
fprintf(fid2,'%s\n'  ,['Makeorthocenters                    = 0                                  # 1=yes, 0=no switch from circumcentres to orthocentres in geominit']);
fprintf(fid2,'%s\n'  ,['Dcenterinside                       = 1.                                 # limit cell center; 1.0:in cell <-> 0.0:on c/g']);
fprintf(fid2,'%s\n\n',['PartitionFile                       =                                    # *_part.pol, polyline(s) x,y']);

fprintf(fid2,'%s\n'  ,'[numerics]');
fprintf(fid2,'%s\n'  ,['CFLMax                              = 0.7                                # Max. Courant nr.']);
fprintf(fid2,'%s\n'  ,['CFLWaveFrac                         = 0.1                                # Wave velocity fraction, total courant vel = u + cflw*wavevelocity']);
fprintf(fid2,'%s\n'  ,['AdvecType                           = 3                                  # Adv type, 0=no, 1= Wenneker, qu-udzt, 2=1, q(uio-u), 3=Perot q(uio-u), 4=Perot q(ui-u), 5=Perot q(ui-u) without itself']);
fprintf(fid2,'%s\n'  ,['Icoriolistype                       =                                    # Coriolis weigthing, icorio 0, 4 or 5 (4 : v = tang. comp of total velocity at u-point, continuity weighted ucnxq, ucnyq = default so far, 5 : 4, scaling down coriolis force below trshcorio)']);
fprintf(fid2,'%s\n'  ,['Limtyphu                            =                                    # Limiter type for waterdepth in continuity eq., 0=no, 1=minmod,2=vanLeer,3=Kooren,4=Monotone Central']);
fprintf(fid2,'%s\n'  ,['Limtypmom                           = 4                                  # Limiter type for cell center advection velocity, 0=no, 1=minmod,2=vanLeer,3=Kooren,4=Monotone Central']);
fprintf(fid2,'%s\n'  ,['Limtypsa                            = 4                                  # Limiter type for salinity transport,           0=no, 1=minmod,2=vanLeer,3=Kooren,4=Monotone Central']);
fprintf(fid2,'%s\n'  ,['Icgsolver                           = 4                                  # Solver type , 1 = sobekGS_OMP, 2 = sobekGS_OMPthreadsafe, 3 = sobekGS, 4 = sobekGS + Saadilud, 5 = parallel/global Saad, 6 = parallel/Petsc, 7 = parallel/GS']);
fprintf(fid2,'%s\n'  ,['Maxdegree                           = 6                                  # Maximum degree in Gauss elimination']);
fprintf(fid2,'%s\n'  ,['Thindykescheme                      =                                    # 0 = no, 1 = compact stencil, 2 = whole tile lifted, full subgrid weir + factor']);
fprintf(fid2,'%s\n'  ,['Thindykecontraction                 =                                    # flow width = flow width*thindykecontraction']);
fprintf(fid2,'%s\n'  ,['Izbndpos                            = 0                                  # Position of z boundary, 0=D3Dflow, 1=on net boundary, 2 = on specified polyline']);
fprintf(fid2,'%s\n'  ,['Jbasqbnddownwindhs                  = 0                                  # 0 : original hu on qbnd, 1 = downwind hs on qbnd']);
fprintf(fid2,'%s\n'  ,['Tlfsmo                              = 0                                  # Fourier smoothing time on waterlevel boundaries (s) ']);
fprintf(fid2,'%s\n'  ,['Slopedrop2D                         = 0.3                                # Apply droplosses only if local bottom slope > Slopedrop2D ']);
fprintf(fid2,'%s\n'  ,['Lincontin                           =                                    # 0 = no, 1 = yes linear continuity']);
fprintf(fid2,'%s\n'  ,['Jaembed1D                           =                                    # 1 : use embedded 1d channels, first run: Operations: Embed 1D channels']);
fprintf(fid2,'%s\n'  ,['Chkadvd                             =                                    # Check advection terms if depth < chkadvdp, => less setbacks ']);
fprintf(fid2,'%s\n'  ,['Teta0                               =                                    # Teta of time integration, 0.5 < Teta < 1d0']);
fprintf(fid2,'%s\n\n',['Maxitverticalforester               =                                    # 0 : no vertical filter, > 0 = Max nr of iterations']);

if strcmp(mdfkeywds.roumet,'C'); roumet = 0; end;
if strcmp(mdfkeywds.roumet,'M'); roumet = 1; end;
if strcmp(mdfkeywds.roumet,'W'); roumet = 2; end;
if isfield(mdfkeywds,'ccofu') == 0; 
    mdfkeywds.ccofu = 0;
    mdfkeywds.ccofv = 0;
    rouval = 0;
else
    rouval = (mdfkeywds.ccofu + mdfkeywds.ccofv)/2;
end
if javis == 0;
    if isfield(mdfkeywds,'vicouv') == 0;
        mdfkeywds.vicouv = 0;
        mdfkeywds.dicouv = 0;
    elseif isempty(mdfkeywds.vicouv) == 1;
        mdfkeywds.vicouv = 0;
        mdfkeywds.dicouv = 0;
    end
end
if isfield(mdfkeywds,'vicoww') == 0;
    mdfkeywds.vicoww = 0;
elseif isempty(mdfkeywds.vicoww) == 1;
    mdfkeywds.vicoww = 0;
end
if isfield(mdfkeywds,'dicoww') == 0; 
    mdfkeywds.dicoww = 0;
elseif isempty(mdfkeywds.dicoww) == 1;
    mdfkeywds.dicoww = 0;
end
if isfield(mdfkeywds,'rhow') == 0; 
    mdfkeywds.rhow = 0;
end
if isfield(mdfkeywds,'salw') == 0; 
    mdfkeywds.salw = 0;
    jasal          = 0;
else
    jasal          = 1;
end
jasal              = 0; % temp!

fprintf(fid2,'%s\n'  ,'[physics]');
if jargh == 1;
    fprintf(fid2,'%s\n'  ,['UnifFrictCoef                       =                                    # Uniform friction coefficient, 0=no friction']);
else
    fprintf(fid2,'%s\n'  ,['UnifFrictCoef                       = ',num2str(rouval,'%5.7f')                         ,'                          # Uniform friction coefficient, 0=no friction']);
end
fprintf(fid2,'%s\n'  ,['UnifFrictType                       = ',num2str(roumet,'%1.0f')                         ,'                                  # 0=Chezy, 1=Manning, 2=White Colebrook, 3=z0 etc']);
fprintf(fid2,'%s\n'  ,['Umodlin                             =                                    # linear friction umod, ifrctyp 4,5,6']);
if javis == 1;
    fprintf(fid2,'%s\n'  ,['Vicouv                              =                                    # Uniform horizontal eddy viscosity (m2/s)']);
    fprintf(fid2,'%s\n'  ,['Dicouv                              =                                    # Uniform horizontal eddy diffusivity (m2/s)']);
else
    fprintf(fid2,'%s\n'  ,['Vicouv                              = ',num2str(mdfkeywds.vicouv,'%5.7f')               ,'                          # Uniform horizontal eddy viscosity (m2/s)']);
    fprintf(fid2,'%s\n'  ,['Dicouv                              = ',num2str(mdfkeywds.dicouv,'%5.7f')               ,'                          # Uniform horizontal eddy diffusivity (m2/s)']);
end
fprintf(fid2,'%s\n'  ,['Vicoww                              = ',num2str(mdfkeywds.vicoww,'%5.7f')               ,'                          # Uniform vertical eddy viscosity (m2/s)']);
fprintf(fid2,'%s\n'  ,['Dicoww                              = ',num2str(mdfkeywds.dicoww,'%5.7f')               ,'                          # Uniform vertical eddy diffusivity (m2/s)']);
fprintf(fid2,'%s\n'  ,['Smagorinsky                         = 0.                                 # Add Smagorinsky horizontal turbulence : vicu = vicu + ( (Smagorinsky*dx)**2)*S, e.g. 0.1']);
fprintf(fid2,'%s\n'  ,['Elder                               = 0.                                 # Add Elder contribution                : vicu = vicu + Elder*kappa*ustar*H/6),   e.g. 1.0']);
fprintf(fid2,'%s\n'  ,['irov                                = ',num2str(mdfkeywds.irov,'%1.0f')                 ,'                                  # 0=free slip, 1 = partial slip using wall_ks']);
fprintf(fid2,'%s\n'  ,['wall_ks                             = 0.                                 # Nikuradse roughness for side walls, wall_z0=wall_ks/30']);
fprintf(fid2,'%s\n'  ,['TidalForcing                        = 1                                  # Tidal forcing (0=no, 1=yes) (only for jsferic == 1)']);
fprintf(fid2,'%s\n'  ,['Salinity                            = ',num2str(jasal,'%1.1d')                          ,'                                  # Include salinity, (0=no, 1=yes)']);
fprintf(fid2,'%s\n'  ,['InitialSalinity                     = ',num2str(mdfkeywds.s0  ,'%5.7f')                 ,'                          # Initial salinity concentration (ppt)']);
fprintf(fid2,'%s\n\n',['Rhomean                             = ',num2str(mdfkeywds.rhow,'%5.7f')                 ,'                       # Average water density (kg/m3)']);

% fprintf(fid2,'%s\n'  ,'[sediment]');
% fprintf(fid2,'%s\n'  ,['Sedimentmodelnr                     =                                    # Sediment model nr, (0=no, 1=Krone, 2=SvR2007)']);
% fprintf(fid2,'%s\n'  ,['Nr_of_sedfractions                  =                                    # Nr of sediment fractions, (specify the next parameters for each fraction) ']);
% fprintf(fid2,'%s\n'  ,['MxgrKrone                           =                                    # Highest fraction index treated by Krone']);
% fprintf(fid2,'%s\n'  ,['D50                                 =                                    # Mean Sandgrain diameter (m), e.g. 0.0001']);
% fprintf(fid2,'%s\n'  ,['Rhosed                              =                                    # Mean Sandgrain rho (kg/m3) , e.g. 2650']);
% fprintf(fid2,'%s\n'  ,['Ws                                  =                                    # Fall velocity (m/s), e.g. 0.0005 m/s']);
% fprintf(fid2,'%s\n'  ,['Erosionpar                          =                                    # Krone Partheniades erosion parameter, e.g. 0.0001  (kg/(m2s)']);
% fprintf(fid2,'%s\n'  ,['Taucre                              =                                    # Critical shear stress for erosion    (N/m2), e.g. 0.3']);
% fprintf(fid2,'%s\n'  ,['InitialSedimentConcentration        =                                    # Initial sediment concentration (kg /m3)']);
% fprintf(fid2,'%s\n'  ,['Uniformerodablethickness            =                                    # Uniform erodable layer thickness (m)']);
% fprintf(fid2,'%s\n'  ,['Numintverticaleinstein              =                                    # Number of vertical intervals in Einstein integrals ( )']);
% fprintf(fid2,'%s\n'  ,['Jaceneqtr                           =                                    # 1=equilibriumtransport at cell centre, 2= at netnode (default) ( ) ']);
% fprintf(fid2,'%s\n'  ,['Morfac                              =                                    # Morphological acceleration factor (), bottom updates active for morfac > 0, 1d0=realtime, etc']);
% fprintf(fid2,'%s\n\n',['TMorfspinup                         =                                    # Spin up time for morphological adaptations (s)']);

if isfield(mdfkeywds,'rhoa') == 0; 
    mdfkeywds.rhoa = 0;
end
if isfield(mdfkeywds,'pavbnd') == 0; 
    mdfkeywds.pavbnd = 101325.0;
end
if isfield(mdfkeywds,'gapres') == 0; 
    mdfkeywds.gapres = 101325.0;
end

fprintf(fid2,'%s\n'  ,'[wind]');
fprintf(fid2,'%s\n'  ,['ICdtyp                              = ',num2str(length(mdfkeywds.wstres)/2,'%1.0f')     ,'                                  # ( ), Cd = const, 2=S&B 2 breakpoints, 3= S&B 3 breakpoints']);
fprintf(fid2,'%s\n'  ,['Cdbreakpoints                       = ',num2str(mdfkeywds.wstres(1:2:end),'%1.7f\t')    ,'                # ( ),   e.g. 0.00063  0.00723']);
fprintf(fid2,'%s\n'  ,['Windspeedbreakpoints                = ',num2str(mdfkeywds.wstres(2:2:end),'%1.7f\t')    ,'              # (m/s), e.g. 0.0      100.0']);
fprintf(fid2,'%s\n'  ,['Rhoair                              = ',num2str(mdfkeywds.rhoa,'%5.7f')                 ,'                          # Air density (kg/m3)']);
fprintf(fid2,'%s\n'  ,['PavBnd                              = ',num2str(mdfkeywds.pavbnd,'%5.7f')               ,'                     # Average air Pressure on open boundaries, (N/m2), only applied if value > 0']);
fprintf(fid2,'%s\n\n',['Gapres                              = ',num2str(mdfkeywds.gapres,'%5.7f')               ,'                     # Global Atmospheric Pressure, (N/m2)']);

% fprintf(fid2,'%s\n'  ,'[waves]');
% fprintf(fid2,'%s\n'  ,['Wavemodelnr                         =                                    # Wave model nr, 0=no, 1=fetch/depth limited hurdlestive, 2=youngverhagen, 3 = SWAN']);
% fprintf(fid2,'%s\n'  ,['Wavenikuradse                       =                                    # Wave friction Nikuradse ks coefficient, used in Krone Swart, e.g. 0.01 m ']);
% fprintf(fid2,'%s\n\n',['Rouwav                              =                                    # Friction model for wave induced shear stress']);

fprintf(fid2,'%s\n'  , '[time]');
fprintf(fid2,'%s\n'  ,['RefDate                             = ',itdate            ,'                           # Reference date (yyyymmdd)']);
fprintf(fid2,'%s\n'  ,['Tunit                               = ',mdfkeywds.tunit   ,'                                  # Time units in MDU (H, M or S)']);
fprintf(fid2,'%s\n'  ,['DtUser                              = ',num2str(mdfkeywds.dt*60,'%5.7f'),'                         # User timestep in seconds (interval for external forcing update & his/map output)']);
fprintf(fid2,'%s\n'  ,['DtMax                               = ',num2str(mdfkeywds.dt*60,'%5.7f'),'                         # Max timestep in seconds']);
fprintf(fid2,'%s\n'  ,['DtInit                              = 1.                                 # Initial timestep in seconds']);
fprintf(fid2,'%s\n'  ,['AutoTimestep                        = 1                                  # Use CFL timestep limit or not (1/0)']);
fprintf(fid2,'%s\n'  ,['TStart                              = ',num2str(mdfkeywds.tstart,'%1.7e'),'                      # Start time w.r.t. RefDate (in TUnit)']);
fprintf(fid2,'%s\n\n',['TStop                               = ',num2str(mdfkeywds.tstop ,'%1.7e'),'                      # Stop  time w.r.t. RefDate (in TUnit)']);

fprintf(fid2,'%s\n'  ,['[restart]']);
fprintf(fid2,'%s\n'  ,['RestartFile                         =                                    # Restart file, only map, hence: *_map.nc']);
fprintf(fid2,'%s\n\n',['RestartDateTime                     = ',itdate,'000000','                     # Restart time (yyyymmddhhmmss)']);

fprintf(fid2,'%s\n'  ,['[external forcing]']);
fprintf(fid2,'%s\n\n',['ExtForceFile                        = ',extfile,'                         # *.ext']);

fprintf(fid2,'%s\n'  ,['[output]']);
fprintf(fid2,'%s\n'  ,['OutputDir                           =                                    # Output directory of map-, his-, rst-, dat- and timings-files, default: DFM_OUTPUT_<modelname>. Set to . for no dir/current dir.']);
fprintf(fid2,'%s\n'  ,['ObsFile                             = ',obsfile,'                                   # *.xyn Coords+name of observation stations.']);
fprintf(fid2,'%s\n'  ,['CrsFile                             = ',crsfile,'                                   # *_crs.pli Polyline(s) definining cross section(s).']);
fprintf(fid2,'%s\n'  ,['HisFile                             = ',netfile(1:end-7),'_his.nc','                        # *_his.nc History file in NetCDF format.']);
fprintf(fid2,'%s\n'  ,['HisInterval                         = ',num2str(mdfkeywds.flhis(2)*60,'%5.7f'),'                         # History output, given as "interval" "start period" "end period" (s)']);
fprintf(fid2,'%s\n'  ,['XLSInterval                         =                                    # Interval (s) between XLS history']);
fprintf(fid2,'%s\n'  ,['FlowGeomFile                        =                                    # *_flowgeom.nc Flow geometry file in NetCDF format.']);
fprintf(fid2,'%s\n'  ,['MapFile                             = ',netfile(1:end-7),'_map.nc','                        # *_map.nc Map file in NetCDF format.']);
fprintf(fid2,'%s\n'  ,['MapInterval                         = ',num2str(mdfkeywds.flmap(2)*60,'%5.7f'),'                         # Map file output, given as "interval" "start period" "end period" (s)']);
fprintf(fid2,'%s\n'  ,['MapFormat                           = 1                                  # Map file format, 1: netCDF, 2: Tecplot, 3: netCFD and Tecplot']);
fprintf(fid2,'%s\n'  ,['RstInterval                         =                                    # Interval (s) between map file outputs']);
fprintf(fid2,'%s\n'  ,['S1incinterval                       =                                    # Interval (m) in incremental file for waterlevels S1']);
fprintf(fid2,'%s\n'  ,['WaqFileBase                         =                                    # Basename (without extension) for all Delwaq files to be written.']);
fprintf(fid2,'%s\n'  ,['WaqInterval                         =                                    # Interval (in s) between Delwaq file outputs']);
fprintf(fid2,'%s\n'  ,['StatsInterval                       =                                    # Interval (in s) between simulation statistics output.']);
fprintf(fid2,'%s\n'  ,['WriteBalancefile                    = 0                                  # Write Balancefile, 1=yes, 0=no']);
fprintf(fid2,'%s\n'  ,['SnapshotDir                         =                                    # Directory where snapshots/screendumps are saved.']);
fprintf(fid2,'%s\n'  ,['TimingsInterval                     =                                    # Timings output interval']);
fprintf(fid2,'%s\n'  ,['TimeSplitInterval                   = 0X                                 # Time splitting interval, after which a new output file is started. value+unit, e.g. ''1 M'', valid units: Y,M,D,h,m,s.']);
fprintf(fid2,'%s\n'  ,['MapOutputTimeVector                 =                                    # File (.mpt) containing fixed map output times (s) w.r.t. RefDate']);
fprintf(fid2,'%s\n\n',['FullGridOutput                      = 0                                  # 0:compact, 1:full time-varying grid data']);

% fprintf(fid2,'%s\n'  ,['[equatorial]']);
% fprintf(fid2,'%s\n'  ,['Ampfreeleft                         =                                    # ampl of free wave left']);
% fprintf(fid2,'%s\n'  ,['Ampfreeright                        =                                    # ampl of free wave right']);
% fprintf(fid2,'%s\n'  ,['Ampforcedzerofr                     =                                    # ampl of forced wave at zero friction']);
% fprintf(fid2,'%s\n'  ,['Nmode                               =                                    # wavelength=nmode/Earthradius']);
% fprintf(fid2,'%s\n\n',['Nfreq                               =                                    # wavefreq  =nfreq/(24*3600)']);

fclose all;