joshuas
/
sbeach_toolbox
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Compare commits

base: joshuas:f25299af0bcaf2af2e6c8215f0df9bbe4093542b
Branches Tags
joshuas:master
...
compare: joshuas:cbc72a4da4867f4ed397acce88f9a09d5288290e
Branches Tags
joshuas:master

2 Commits
f25299af0b ... cbc72a4da4

Author SHA1 Message Date
Joshua Simmons cbc72a4da4 Merge branch 'master' of http://git.wrl.unsw.edu.au:3000/joshuas/sbeach_toolbox 7 years ago
Joshua Simmons af986a4297 init commit 7 years ago
45 changed files with 2795 additions and 0 deletions
Show Stats

120
RUN_sbeach_model.m
Unescape Escape View File

@ -0,0 +1,120 @@
% Written by Joshua Simmons 11/02/2015 and updated 17/11/2017
% to convert from fielddata storage struct to SBEACH toolbox.
clear
clc
%% Settings
%profiledata locations
basedatapath = '.\example\input_data\';
setupoutputpath = '.\temp\';
infooutpath = '.\example\setup_info\';
sbeachexepath = '.\sbeach_exe\final_exe\SBEACH.exe';
%meta
site = 'narrabeen';
stormid = '2015';
%write file to
writetofile = 1; %boolean 1 = write to file, 0 = just SBEACHout.
% other parameters;
sbts = 2; %time step length minutes
%% Code - format the data
%load profile data - account for different data names - lin/poly
pfdata = load([basedatapath 'morphology\processed\morphology_' site '_' stormid '.mat']);
fns1 = fieldnames(pfdata);
pfdata = pfdata.(fns1{1});
%load wavedata
wavedata = load([basedatapath 'fielddata\wave\processed\wavedata_' site '_' stormid '.mat']);
fns2 = fieldnames(wavedata);
wavedata = wavedata.(fns2{1});
%load tidedata
tidedata = load([basedatapath 'fielddata\tide\processed\tidedata_' site '_' stormid '.mat']);
fns3 = fieldnames(tidedata);
tidedata = tidedata.(fns3{1});
% load the profileinfo file
profileinfo = load([infooutpath 'profileinfo_' site '.mat']);
fns2 = fieldnames(profileinfo);
profileinfo = profileinfo.(fns2{1});
%check for dimension match
if length(pfdata.data) ~= length(wavedata.ns.data)
error('Wave data information does not match the morphology data')
end
%get the depth contour
sbdepth = wavedata.ns.contour(1:2);
pathrun = cell(length(pfdata.data),1);
%do the modelling
fprintf('Setting up %d SBEACH models...', length(pfdata.data))
for ii = 1:length(pfdata.data)
if ~isnan(pfdata.data(ii).initial.x)
fprintf('%d..',ii)
pfid = ['PF' num2str(ii)];
%get profile info
pfx = pfdata.data(ii).initial.x;
pfz = pfdata.data(ii).initial.z;
%angle of pf from true north
pfangle = profileinfo.data(ii).TN; %angle from true north
outpath = [setupoutputpath 'sb_' pfid '\'];
% outpath = [setupoutputpath 'sb_' site '_' stormid '_' pfid '\'];
sb_model = fsbeach_setupmodel(...
'bathy', { 'x', pfx, 'z', pfz}, ...
'waves', { 'time', wavedata.ns.data(ii).datenums, 'Hsig', wavedata.ns.data(ii).Hsig, 'Tp1', wavedata.ns.data(ii).Tp1, ...
'ang', wavedata.ns.data(ii).Wdir, 'pfang', pfangle}, ...
'tide', { 'time', tidedata.time,'elev', tidedata.tide}, ...
'settings', { 'DT', sbts, 'K', 2.5e-6, 'BMAX', 28.5, 'D50', 0.455 ...
}, ...
'path', outpath, ...
'name', ['sb_' pfid] ...
);
% variable grid can be specified in settings with e.g.:
% 'NGRID', 3, 'DXV', [1 5 10], 'NDXV', [115 30 28]
% hard bottom can be specified in bathy with e.g.:
% 'hbz', pfz-1
measured(ii).x = pfdata.data(ii).measured.x;
measured(ii).z = pfdata.data(ii).measured.z;
pathrun{ii} = sb_model.path;
end
end
fprintf('\n')
fprintf('Running %d SBEACH models...', length(pathrun))
%now run sbeach EXE on the files setup
for ii = 1:length(pathrun)
fprintf('%d..',ii)
[fail, cmdout] = fsbeach_runmodel(pathrun{ii},sbeachexepath);
if fail
error(cmdout);
end
end
fprintf('\n')
fprintf('Collecting SBEACH results...\n')
%now collect the sbeach results
sbout = fsbeach_collectresults(pathrun);
fprintf('Getting measured data...\n')
%attached the measured data
sbout.measured = measured;
fprintf('Processing SBEACH results...\n')
%evaulate results
sbout = fsbeach_processresults(sbout);
fprintf('Done!\n')

169
SBEACH_GLUE_setup.m
Unescape Escape View File

@ -0,0 +1,169 @@
% Written by Joshua Simmons 11/02/2015 and updated on 11/01/2016
% to convert from fielddata storage struct to SBEACH compatible format.
clear
clc
%% Settings
%all the initial files location
%specify folder with Mitch's bathymetry and field data setup for xbeach.
setupoutputpath = 'D:\7 Erosion Models\a_nb_model_comparison\sbeach\model_setup\setup_folders\';
setupbasepath = 'D:\7 Erosion Models\a_nb_model_comparison\sbeach\model_setup\';
basebatchrunpath = 'D:\7 Erosion Models\a_nb_model_comparison\sbeach\model_results\src\';
% meta
site = 'deewhy';
% stormid = '2007';
% stormid = '2011';
% stormid = '2014';
stormid = '2015';
%profile info - CHECK
pfno = 5;
% model parameters
%parameter vals - leave as '' if params to be generated
paramload = 'D:\7 Erosion Models\a_nb_model_comparison\sbeach\params\bigparams_K_EPS_BMAX.mat';
% paramload = '';
%or specify stats
paramsname = {'K','EPS','BMAX'}; %K divided by 3600 below to bring into seconds, Cs for R2;
% paramsrange = {[0.5e-6 2.5e-6],[0.001 0.003],[15 30]}; % give range as [min max] then fed into rand()
paramsrange = {[1.75e-6 1.75e-6],[0.002 0.002],[28.5 28.5]}; % give range as [min max] then fed into rand()
nruns = 10000; %number of parameter conmbinations
nrunst = 1;
nrunend = 10000;
%% Code - format the data
%generate param combos
if isempty(paramload)
for aa = 1:length(paramsname)
tmpmin = paramsrange{aa}(1);
tmpmax = paramsrange{aa}(2);
params(:,aa) = (tmpmax-tmpmin).*rand(nruns,1) + tmpmin;
end
bigparams.name = paramsname;
bigparams.value = params;
else
bigparams = load(paramload);
fnsparams = fieldnames(bigparams);
bigparams = bigparams.(fnsparams{1,1});
params = bigparams.value;
if sum(strcmp(bigparams.name,paramsname)) ~= length(paramsname)
error('Parameter names in loaded file do not match those specified')
end
end
%load pf info
profileinfo = load([setupbasepath 'setup_info\profileinfo_' site '.mat']);
fns2 = fieldnames(profileinfo);
profileinfo = profileinfo.(fns2{1});
%load cfg info
cfginfo = load([setupbasepath 'setup_info\cfginfo_' site '_' stormid '_PF' num2str(pfno) '.mat']);
fns3 = fieldnames(cfginfo);
cfginfo = cfginfo.(fns3{1});
%get the names of the files copying across
dns1 = dir([setupoutputpath 'sb_' site '_' stormid '_PF' num2str(pfno) '\']);
dns1dirs = {dns1.name}; % Get a list of the subdirectories
dns1ind2 = ~ismember(dns1dirs,{'.','..'}); % Find index of subdirectories
dns1 = dns1(dns1ind2);
allfilelocs = cell(nruns,1);
% not done VVVV
%do the modelling
for ii = nrunst:nrunend
CFGname = ['PF' num2str(pfno) '_Run_' num2str(ii) '.CFG'];
tmpprintdir = [basebatchrunpath site '_' stormid '_PF' num2str(pfno) '\PF' num2str(pfno) '_Run_' num2str(ii) '\'];
if ~exist(tmpprintdir,'dir')
mkdir(tmpprintdir)
end
%now copy the setup files across
copyfile([setupoutputpath 'sb_' site '_' stormid '_PF' num2str(pfno) '\*'],tmpprintdir);
% now change the contents of CFG file
fin = fopen([setupbasepath 'base.CFG']) ;
fout = fopen([tmpprintdir 'temp.CFG'], 'w');
tline = fgetl(fin);
while ischar(tline)
%replace:
%filename
if (~isempty(strfind(tline, '%RPLCNAME%')))
tline = strrep(tline, '%RPLCNAME%', ['_Run_' num2str(ii)]);
end
%profile ID
if (~isempty(strfind(tline, '%RPLCPFID%')))
tline = strrep(tline, '%RPLCPFID%', ['PF' num2str(pfno)]);
end
%run number
if (~isempty(strfind(tline, '%RPLCRNUM%')))
tline = strrep(tline, '%RPLCRNUM%', num2str(ii));
end
%calculation cells
%number
if (~isempty(strfind(tline, '%RPLCNDX%')))
tline = strrep(tline, '%RPLCNDX%', num2str(cfginfo.NDX));
end
%landward boundary
if (~isempty(strfind(tline, '%RPLCXSTART%')))
tline = strrep(tline, '%RPLCXSTART%', num2str(cfginfo.XSTART));
end
%spacing between
if (~isempty(strfind(tline, '%RPLCDXC%')))
tline = strrep(tline, '%RPLCDXC%', num2str(cfginfo.DXC));
end
%calculation time
%number
if (~isempty(strfind(tline, '%RPLCNDT%')))
tline = strrep(tline, '%RPLCNDT%', num2str(cfginfo.NDT));
end
%steps
if (~isempty(strfind(tline, '%RPLCDT%')))
tline = strrep(tline, '%RPLCDT%', num2str(cfginfo.DT));
end
%steps
if (~isempty(strfind(tline, '%RPLCD50%')))
tline = strrep(tline, '%RPLCD50%', num2str(profileinfo.data(pfno).d50*1000));
end
%check for params
for aa = 1:length(paramsname)
if (~isempty(strfind(tline, ['%RPLC' paramsname{aa} '%'])))
tline = strrep(tline, ['%RPLC' paramsname{aa} '%'], num2str(params(ii,aa)));
end
end
%now print the final string
fprintf(fout,'%s\r\n',tline);
tline = fgetl(fin);
end
fclose(fin);
fclose(fout);
%rename CFG to original name
try
movefile([tmpprintdir 'temp.CFG'],[tmpprintdir CFGname])
catch
fileattrib([tmpprintdir 'temp.CFG'],'+w')
movefile([tmpprintdir 'temp.CFG'],[tmpprintdir CFGname])
end
%change file name
for jj = 1:length(dns1)
movefile([tmpprintdir dns1(jj).name],[tmpprintdir 'PF' num2str(pfno) '_Run_' num2str(ii) dns1(jj).name(end-3:end)])
end
%save location of file
allfilelocs{ii,1} = ['PF' num2str(pfno) '_Run_' num2str(ii) '\PF' num2str(pfno) '_Run_' num2str(ii)];
end
%write allfilelocs to .txt
fid = fopen([basebatchrunpath site '_' stormid '_PF' num2str(pfno) '\SBEACHrun_loc.txt'],'w');
for kk = 1:length(allfilelocs)
fprintf(fid,'%s\r\n',allfilelocs{kk});
end
fclose(fid);

140
SBEACH_storm_setup.m
Unescape Escape View File

@ -0,0 +1,140 @@
% Written by Joshua Simmons 11/02/2015 and updated 20/01/2016
% to convert from fielddata storage struct to SBEACH compatible format.
clear
clc
%% Settings
%profiledata locations
basedatapath = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\';
setupoutputpath = 'D:\7 Erosion Models\a_nb_model_comparison\sbeach\model_setup\setup_folders\';
infooutpath = 'D:\7 Erosion Models\a_nb_model_comparison\sbeach\model_setup\setup_info\';
%meta
site = 'monavale';
% stormid = '2007';
% stormid = '2011';
% stormid = '2014';
stormid = '2015';
%write file to
writetofile = 1; %boolean 1 = write to file, 0 = just SBEACHout.
% other parameters;
sbts = 5; %time step length minutes
%% Code - format the data
%load profile data - account for different data names - lin/poly
pfdata = load([basedatapath 'morphology\processed\morphology_' site '_' stormid '.mat']);
fns1 = fieldnames(pfdata);
pfdata = pfdata.(fns1{1});
%load wavedata
wavedata = load([basedatapath 'fielddata\wave\processed\wavedata_' site '_' stormid '.mat']);
fns2 = fieldnames(wavedata);
wavedata = wavedata.(fns2{1});
%load tidedata
tidedata = load([basedatapath 'fielddata\tide\processed\tidedata_' site '_' stormid '.mat']);
fns3 = fieldnames(tidedata);
tidedata = tidedata.(fns3{1});
% load the profileinfo file
profileinfo = load([infooutpath 'profileinfo_' site '.mat']);
fns2 = fieldnames(profileinfo);
profileinfo = profileinfo.(fns2{1});
%check for dimension match
if length(pfdata.data) ~= length(wavedata.ns.data)
error('Wave data information does not match the morphology data')
end
%check start time
if tidedata.time(1) ~= wavedata.os.datenums(1)
error('Start times for wave and tide data dont match')
end
if tidedata.time(end) ~= wavedata.os.datenums(end)
error('End times for wave and tide data dont match. Make sure start/end times are hourly (for wave data).')
end
%check and allocate time parameters;
sbtidetime = round((tidedata.time - tidedata.time(1))*24*60);
%check time steps
diffsbtidetime = diff(sbtidetime);
if any(diffsbtidetime~=diffsbtidetime(1))
error('Inconsistent time steps in tide data');
end
%now check the wave data lines up (os only required the rest if copied)
sbwavetime = round((wavedata.os.datenums - wavedata.os.datenums(1))*24*60);
diffsbwavetime = diff(sbwavetime);
if any(diffsbwavetime~=diffsbwavetime(1))
error('Inconsistent time steps in wave data');
end
%get the depth contour
sbdepth = wavedata.ns.contour(1:2);
%do the modelling
for ii = 1:length(pfdata.data)
if ~isnan(pfdata.data(ii).initial.x)
pfid = ['PF' num2str(ii)];
%get profile info
prifile = [pfdata.data(ii).initial.x pfdata.data(ii).initial.z]; %initial pf
%angle of pf from true north
pfangle = profileinfo.data(ii).TN; %angle from true north
%create the other information required for files
% for .WAV and .ANG file - wave height/period and angle files
wavfile = [wavedata.ns.data(ii).Hsig wavedata.ns.data(ii).Tp1];
angfile = -1*(wavedata.ns.data(ii).Wdir-pfangle); %correct to SBEACH wave angle convention as per report 3 pg 19
wavets = sbwavetime(2)-sbwavetime(1); %wave time step minutes
% for .ELV file - water elevation file
elvfile = tidedata.tide;
tidets = sbtidetime(2)-sbtidetime(1); %tide time step
totaltimewav = sbwavetime(end);
totaltimetide = sbtidetime(end);
if totaltimewav - totaltimetide ~= 0
disp('Warning: difference between wave and tide total times')
end
totaltime = totaltimetide;
totalmodeltimesteps = totaltime/sbts;
disp([pfid ':'])
disp(['A.10 - Total Time steps: ' num2str(totalmodeltimesteps)])
disp(['B.6 - DTWAV: ' num2str(wavets)])
disp(['B.11 - DTANG: ' num2str(wavets)])
disp(['B.12 - DMEAS: ' sbdepth])
disp(['B.20 - DTELEV: ' num2str(tidets)])
outpath = [setupoutputpath 'sb_' site '_' stormid '_' pfid '\'];
cfginfo.NDT = totalmodeltimesteps;
cfginfo.DT = sbts;
cfginfo.DXC = pfdata.data(ii).initial.x(2)-pfdata.data(ii).initial.x(1);
cfginfo.NDX = length(pfdata.data(ii).initial.x)-1;
[cfginfo.XSTART,~] = min(pfdata.data(ii).initial.x);
cfginfo.DTWAV = wavets;
cfginfo.DTANG = wavets;
cfginfo.DMEAS = str2double(sbdepth);
cfginfo.DTELEV = tidets;
%write to files
if writetofile
if ~exist(outpath, 'dir')
mkdir(outpath)
else
delete([outpath '\*'])
end
writetofilepri(prifile, outpath, [site '_' stormid '_' pfid])
writetofilewav(wavfile, outpath, [site '_' stormid '_' pfid])
writetofileang(angfile, outpath, [site '_' stormid '_' pfid])
writetofileelv(elvfile, outpath, [site '_' stormid '_' pfid])
save([infooutpath 'cfginfo_' site '_' stormid '_' pfid '.mat'],'cfginfo')
end
end
end

38
calc_TNorientation.m
Unescape Escape View File

@ -0,0 +1,38 @@
% Written by Joshua Simmons 20/01/16 to calculate angle for TN for SBeach
% modelling purposes.
% NOTE: this code assumes that the furthest offshore point (-10m) is the
% x(end) and the furthest inshore point is in x(1) (in line with the xbeach
% co-ordinate system)
clear
clc
% load('D:\7 Erosion Models\a_nb_model_comparison\xbeach\other_data\narra_xyz_for_alfa.mat')
%
% for ii = 1:length(pfxy)
% x = [pfxy(ii).x(1) pfxy(ii).x(end)];
% y = [pfxy(ii).y(1) pfxy(ii).y(end)];
%
% if x(2)-x(1)>0
% intang = atand((y(2)-y(1))/(x(2)-x(1)));
% ang = 90 - intang;
% elseif x(2)-x(1)<0
% intang = atand((y(2)-y(1))/(x(2)-x(1)));
% ang = 270 - intang;
% end
% TN(ii,1) = ang;
% end
%or else get information from site specific file
addpath('J:\Coastal\')
coastalinit
% pfs = {'PF1','PF2','PF4','PF6','PF8'};
% pfs = {'PF1','PF2','PF3'};
pfs = {'PF1','PF2','PF3','PF4','PF5'};
for ii = 1:length(pfs)
[~,rot_angle] = get_rotation_specs('DEEWH',pfs{ii});
alfa(ii,1) = -rot_angle*180/pi;
end

BIN
example/input_data/fielddata/stormdates_narra.mat
View File

Binary file not shown.

BIN
example/input_data/fielddata/stormdates_other.mat
View File

Binary file not shown.

44
example/input_data/fielddata/tide/code/tidedataget.m
Unescape Escape View File

@ -0,0 +1,44 @@
% Written by Joshua Simmons 04/01/2016
% Description: used to tranform offshore to nearshore mike21 using the code
% OS2NS_MIKE21.m
clear
clc
%% Settings
%tide gauge data
guagedataloc = 'J:\Coastal\Data\Tide\Sydney\Processed\SydTideData.mat';
site = 'Narrabeen'; %embayment of interest
%storm details
stormdatesloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\stormdates_narra.mat';
storm = {'2007', '2011', '2014', '2015'};
% save?
savebool = 1;
% datasave location
saveloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\tide\processed\';
%% Code
%load buoy data,profile data and lookup table
guagedata = load(guagedataloc);
stormdates = load(stormdatesloc);
% narrow down to date range
for ii = 1:length(storm)
startdate = stormdates.(['storm_' storm{ii}]).startdate;
enddate = stormdates.(['storm_' storm{ii}]).enddate;
validentries = find(guagedata.SydTideData.time>=startdate & guagedata.SydTideData.time<=enddate);
tidedata.tide = guagedata.SydTideData.tide(validentries);
tidedata.time = guagedata.SydTideData.time(validentries);
if savebool
save([saveloc 'tidedata_' lower(site) '_' storm{ii} '.mat'],'tidedata')
end
clear tidedata startdate enddate validentries
end

44
example/input_data/fielddata/tide/code/tidedataget_other.m
Unescape Escape View File

@ -0,0 +1,44 @@
% Written by Joshua Simmons 04/01/2016
% Description: used to tranform offshore to nearshore mike21 using the code
% OS2NS_MIKE21.m
clear
clc
%% Settings
%tide gauge data
guagedataloc = 'J:\Coastal\Data\Tide\Sydney\Processed\SydTideData.mat';
site = 'Monavale'; %embayment of interest
%storm details
stormdatesloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\stormdates_other.mat';
storm = {'2014', '2015'};
% save?
savebool = 1;
% datasave location
saveloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\tide\processed\';
%% Code
%load buoy data,profile data and lookup table
guagedata = load(guagedataloc);
stormdates = load(stormdatesloc);
% narrow down to date range
for ii = 1:length(storm)
startdate = stormdates.(['storm_' storm{ii}]).startdate;
enddate = stormdates.(['storm_' storm{ii}]).enddate;
validentries = find(guagedata.SydTideData.time>=startdate & guagedata.SydTideData.time<=enddate);
tidedata.tide = guagedata.SydTideData.tide(validentries);
tidedata.time = guagedata.SydTideData.time(validentries);
if savebool
save([saveloc 'tidedata_' lower(site) '_' storm{ii} '.mat'],'tidedata')
end
clear tidedata startdate enddate validentries
end

BIN
example/input_data/fielddata/tide/processed/tidedata_narrabeen_2015.mat
View File

Binary file not shown.

105
example/input_data/fielddata/wave/code/fmikeos2ns.m
Unescape Escape View File

@ -0,0 +1,105 @@
function nswaves = fmikeos2ns(oswaves, LookUps, site, contour, location)
% Written by Joshua Simmons (WRL) 03/11/2014
% to transform waves from offshore (OS) to nearshore (NS) using MIKE21 model
% oswaves - OS wave data - oswaves.Hsig(Significant wave height)/
% oswaves.Tp1(Peak period)/oswaves.Wdir(Incident wave direction)
% LookUps - MIKE21 wave model output in file "OEH_ARC_PointOutput.mat"
% site - embayment of interest - i.e. 'Narrabeen'
% contour - 10 or 15 (m) contour for NS output
% location - either pure S co-ord (see Mitchell Harley's xyz2spz
% transformation) or location.E and location.N - refers to output
% location for NS waves (gets nearest MIKE21 point)
%
%V1.1 - updated by Joshua Simmons 21/01/2015 to patch Wamberal/Wanda
%error.
%find site and contour indices
siteidx = getnameidx(LookUps(:,2), site);
%apply fix for Wamberal and Wanda
% from TM on 21/01/2015: "I've found an error in the .mat file for the
% look-up tables; in the pdf I say that all output for all locations runs
% north to south - for Wamberal and Wanda it seems I've output them south to north.
% So if you are using output from these two sites, if you haven't realised already,
% alongshore gradients will look inverted."
for aa = 1:2;
LookUps{6,1}{aa,1} = flipud(LookUps{6,1}{aa,1}); LookUps{6,1}{aa,1}(:,1) = flipud(LookUps{6,1}{aa,1}(:,1));
LookUps{7,1}{aa,1} = flipud(LookUps{7,1}{aa,1}); LookUps{7,1}{aa,1}(:,1) = flipud(LookUps{7,1}{aa,1}(:,1));
end
if contour == 10
contouridx = 1;
elseif contour == 15;
contouridx = 2;
else
error('Contour information not available, select either 10 or 15m contour')
end
%extract scenarios for easy lookup
scenariosidx = getnameidx(LookUps(:,2), 'Offshore Wave Scenarios');
scenarios = LookUps{scenariosidx,1};
%get nearest location available for output
Eouts = LookUps{siteidx,1}{contouridx,1}(:,2);
Nouts = LookUps{siteidx,1}{contouridx,1}(:,3);
if strcmp(site,'Narrabeen')
%compares closes alongshore spacing (s) for log-spiral narrabeen
xyz_data.x = cell2mat(Eouts);
xyz_data.y = cell2mat(Nouts);
xyz_data.z = zeros(size(Eouts));
spz = xyz2spz(xyz_data,'NARRA');
[~,locationidx] = min(abs(location-spz.s));
else
%gets nearest output point
dists = sqrt((location.E-cell2mat(Eouts)).^2+(location.N-cell2mat(Nouts)).^2);
[~,locationidx] = min(abs(dists));
end
fcount = 1;
%loop through oswaves and lookup nswaves
for ii = 1:length(oswaves.Hsig)
tmpHs = oswaves.Hsig(ii);
tmpTp = oswaves.Tp1(ii);
tmpWdir = oswaves.Wdir(ii);
%find closest scenario match
[~,Hsidx] = min(abs(tmpHs-scenarios(:,1)));
tmpHsround = scenarios(Hsidx,1);
[~,Tpidx] = min(abs(tmpTp-scenarios(:,2)));
tmpTpround = scenarios(Tpidx,2);
[~,Wdiridx] = min(abs(tmpWdir-scenarios(:,3)));
tmpWdirround = scenarios(Wdiridx,3);
iter1ind = find(scenarios(:,3) == tmpWdirround);
iter2ind = find(scenarios(iter1ind,2) == tmpTpround);
iter3ind = find(scenarios(iter1ind(iter2ind),1) == tmpHsround);
iter4ind = iter1ind(iter2ind(iter3ind));
%now get ns lookup information if available
if isempty(iter4ind)
%get the closest wave height to scenario chosen
[~,iter3ind] = min(abs(tmpHs-scenarios(iter1ind(iter2ind),1)));
iter4ind = iter1ind(iter2ind(iter3ind));
diffHs = scenarios(iter4ind,1)-tmpHs;
disp(['No match found for oswaves(' num2str(ii) ') wave height off by ' num2str(diffHs) 'm']);
end
nstmp = LookUps{siteidx,1}{contouridx,1}{locationidx,4}(iter4ind,:);
%store the values from the lookup table
nswaves.Hsig(ii,1) = nstmp(1);
nswaves.Hmax(ii,1) = nstmp(2);
nswaves.Tp1(ii,1) = nstmp(3);
nswaves.T01(ii,1) = nstmp(4);
nswaves.Wdir(ii,1) = nstmp(5);
nswaves.power(ii,1) = nstmp(6);
end
try
nswaves.date = oswaves.date;
nswaves.time = oswaves.time;
nswaves.datenums = oswaves.datenums;
catch
nswaves.datenums = oswaves.datenums;
end
end

95
example/input_data/fielddata/wave/code/run_OS_to_NS_mike21.m
Unescape Escape View File

@ -0,0 +1,95 @@
% Written by Joshua Simmons 04/01/2016
% Description: used to tranform offshore to nearshore mike21 using the code
% OS2NS_MIKE21.m
clear
clc
%% Settings
%buoy data
buoydataloc = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed_CSIROGapFill.mat';
% the non CSIRO filtered data is the only data available after 2014 so fill
% from the recent data with gaps.
backupbuoydata = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed.mat';
%location of mike21 runcode
mike21loc = 'J:\Coastal\Tools\Nearshore Wave Transformations\MIKE21\';
%location of MALT
maltloc = 'J:\Coastal\Tools\MALT Logspiral Transformation\';
%location of profiles for which wave data is required.
profiledataloc = [437 850 1750 2606 3420]; % s co-ordinate
site = 'Narrabeen'; %embayment of interest
contour = 10; %m - 10m or 15m contour?
%storm details
stormdatesloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\stormdates_narra.mat';
storm = {'2007', '2011', '2014', '2015'};
% save?
savebool = 1;
% datasave location
saveloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\wave\processed\';
%% Code
% for MIKE21 code and MALT
addpath(mike21loc);
addpath(maltloc);
%load buoy data,profile data and lookup table
buoydata = load(buoydataloc);
backupbuoydata = load(backupbuoydata);
load('OEH_ARC_PointOutput.mat');
stormdates = load(stormdatesloc);
% narrow down to date range
for ii = 1:length(storm)
startdate = stormdates.(['storm_' storm{ii}]).startdate;
enddate = stormdates.(['storm_' storm{ii}]).enddate;
validentries = find(buoydata.dates>=startdate & buoydata.dates<=enddate);
% if nothing then revert to backup
if isempty(validentries) || validentries(end) == length(buoydata.dates)
validentries = find(backupbuoydata.dates>=startdate & backupbuoydata.dates<=enddate);
oswaves.Hsig = backupbuoydata.Hsig(validentries);
oswaves.Tp1 = backupbuoydata.Tp(validentries);
oswaves.Wdir = backupbuoydata.Wdir(validentries);
oswaves.datenums = backupbuoydata.dates(validentries);
else
oswaves.Hsig = buoydata.Hsig(validentries);
oswaves.Tp1 = buoydata.Tp(validentries);
oswaves.Wdir = buoydata.Wdir(validentries);
oswaves.datenums = buoydata.dates(validentries);
end
fns1 = fieldnames(oswaves);
for kk = 1:length(fns1)
if any(isnan(oswaves.(fns1{kk})))
error('NaN values in wave data')
end
end
tic
parfor jj = 1:length(profiledataloc)
location = profiledataloc(jj);
nswaves = OS2NS_MIKE21(oswaves, LookUps, site, contour, location);
% bug fix to ensure no inaccurate wave period estimates
nswaves.Tp1 = oswaves.Tp1;
tmp(jj) = nswaves;
end
wavedata.ns.data = tmp;
toc
wavedata.os = oswaves;
wavedata.ns.model = 'MIKE21';
wavedata.ns.contour = [num2str(contour) 'm'];
wavedata.ns.location = profiledataloc;
if savebool
save([saveloc 'wavedata_' lower(site) '_' storm{ii} '.mat'],'wavedata')
end
clear wavedata oswaves startdate enddate validentries tmp
end

99
example/input_data/fielddata/wave/code/run_OS_to_NS_mike21_deewhy.m
Unescape Escape View File

@ -0,0 +1,99 @@
% Written by Joshua Simmons 04/01/2016
% Description: used to tranform offshore to nearshore mike21 using the code
% OS2NS_MIKE21.m
clear
clc
%% Settings
%buoy data
buoydataloc = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed_CSIROGapFill.mat';
% the non CSIRO filtered data is the only data available after 2014 so fill
% from the recent data with gaps.
backupbuoydata = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed.mat';
%location of mike21 runcode
mike21loc = 'J:\Coastal\Tools\Nearshore Wave Transformations\MIKE21\';
%location of MALT
maltloc = 'J:\Coastal\Tools\MALT Logspiral Transformation\';
%location of profiles for which wave data is required.
% profiledataloc = [6275797.660794820 6275692.739810870 6275536.815570840]; % Northings co-ordinate - Bilgola (at 10m contour)
% profiledataloc = [6272148.185 6272338.391 6272389.496 6271924.185 6271823.724]; % Northings co-ordinate - Mona Vale (at 10m contour)
profiledataloc = [6264042.586 6264127.494 6263995.016 6263745.956 6263879.058]; % Northings co-ordinate - Dee Why (at 10m contour)
profiledataloc2 = [343633.9263 343038.1544 342940.1201 342757.2994 342848.2807]; % Eastings co-ordinate - Dee Why (at 10m contour)
site = 'DeeWhy'; %embayment of interest
contour = 10; %m - 10m or 15m contour?
%storm details
stormdatesloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\stormdates_other.mat';
storm = {'2014', '2015'};
% save?
savebool = 1;
% datasave location
saveloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\wave\processed\';
%% Code
% for MIKE21 code and MALT
addpath(mike21loc);
addpath(maltloc);
%load buoy data,profile data and lookup table
buoydata = load(buoydataloc);
backupbuoydata = load(backupbuoydata);
load('OEH_ARC_PointOutput.mat');
stormdates = load(stormdatesloc);
% narrow down to date range
for ii = 1:length(storm)
startdate = stormdates.(['storm_' storm{ii}]).startdate;
enddate = stormdates.(['storm_' storm{ii}]).enddate;
validentries = find(buoydata.dates>=startdate & buoydata.dates<=enddate);
% if nothing then revert to backup
if isempty(validentries) || validentries(end) == length(buoydata.dates)
validentries = find(backupbuoydata.dates>=startdate & backupbuoydata.dates<=enddate);
oswaves.Hsig = backupbuoydata.Hsig(validentries);
oswaves.Tp1 = backupbuoydata.Tp(validentries);
oswaves.Wdir = backupbuoydata.Wdir(validentries);
oswaves.datenums = backupbuoydata.dates(validentries);
else
oswaves.Hsig = buoydata.Hsig(validentries);
oswaves.Tp1 = buoydata.Tp(validentries);
oswaves.Wdir = buoydata.Wdir(validentries);
oswaves.datenums = buoydata.dates(validentries);
end
fns1 = fieldnames(oswaves);
for kk = 1:length(fns1)
if any(isnan(oswaves.(fns1{kk})))
error('NaN values in wave data')
end
end
tic
for jj = 1:length(profiledataloc)
location.N = profiledataloc(jj);
location.E = profiledataloc2(jj);
nswaves = fmikeos2ns(oswaves, LookUps, site, contour, location);
% bug fix to ensure no inaccurate wave period estimates
nswaves.Tp1 = oswaves.Tp1;
tmp(jj) = nswaves;
end
wavedata.ns.data = tmp;
toc
wavedata.os = oswaves;
wavedata.ns.model = 'MIKE21';
wavedata.ns.contour = [num2str(contour) 'm'];
wavedata.ns.location = profiledataloc;
if savebool
save([saveloc 'wavedata_' lower(site) '_' storm{ii} '.mat'],'wavedata')
end
clear wavedata oswaves startdate enddate validentries tmp
end

105
example/input_data/fielddata/wave/code/run_OS_to_NS_mike21_end_effects_check.m
Unescape Escape View File

@ -0,0 +1,105 @@
% Written by Joshua Simmons 04/01/2016
% Description: used to tranform offshore to nearshore mike21 using the code
% OS2NS_MIKE21.m
clear
clc
%% Settings
%buoy data
buoydataloc = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed_CSIROGapFill.mat';
% the non CSIRO filtered data is the only data available after 2014 so fill
% from the recent data with gaps.
backupbuoydata = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed.mat';
%location of mike21 runcode
mike21loc = 'J:\Coastal\Tools\Nearshore Wave Transformations\MIKE21\';
%location of MALT
maltloc = 'J:\Coastal\Tools\MALT Logspiral Transformation\';
%location of profiles for which wave data is required.
profiledataloc = [90:20:3870]; % s co-ordinate
site = 'Narrabeen'; %embayment of interest
contour = 10; %m - 10m or 15m contour?
%storm details
stormdatesloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\stormdates_narra.mat';
storm = {'2007', '2011', '2014', '2015'};
% save?
savebool = 0;
% datasave location
saveloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\wave\processed\';
%% Code
% for MIKE21 code and MALT
addpath(mike21loc);
addpath(maltloc);
%load buoy data,profile data and lookup table
buoydata = load(buoydataloc);
backupbuoydata = load(backupbuoydata);
load('OEH_ARC_PointOutput.mat');
stormdates = load(stormdatesloc);
% narrow down to date range
for ii = 1:length(storm)
startdate = stormdates.(['storm_' storm{ii}]).startdate;
enddate = stormdates.(['storm_' storm{ii}]).enddate;
validentries = find(buoydata.dates>=startdate & buoydata.dates<=enddate);
% if nothing then revert to backup
if isempty(validentries) || validentries(end) == length(buoydata.dates)
validentries = find(backupbuoydata.dates>=startdate & backupbuoydata.dates<=enddate);
oswaves.Hsig = backupbuoydata.Hsig(validentries);
oswaves.Tp1 = backupbuoydata.Tp(validentries);
oswaves.Wdir = backupbuoydata.Wdir(validentries);
oswaves.datenums = backupbuoydata.dates(validentries);
else
oswaves.Hsig = buoydata.Hsig(validentries);
oswaves.Tp1 = buoydata.Tp(validentries);
oswaves.Wdir = buoydata.Wdir(validentries);
oswaves.datenums = buoydata.dates(validentries);
end
fns1 = fieldnames(oswaves);
for kk = 1:length(fns1)
if any(isnan(oswaves.(fns1{kk})))
error('NaN values in wave data')
end
end
tic
parfor jj = 1:length(profiledataloc)
location = profiledataloc(jj);
nswaves = OS2NS_MIKE21(oswaves, LookUps, site, contour, location);
% bug fix to ensure no inaccurate wave period estimates
nswaves.Tp1 = oswaves.Tp1;
% tmp(jj) = nswaves;
Hsig(:,jj) = nswaves.Hsig(:,1);
Wdir(:,jj) = nswaves.Wdir(:,1);
power(:,jj) = nswaves.power(:,1);
end
dataout(ii).s = profiledataloc;
dataout(ii).Hsig = Hsig;
dataout(ii).Wdir = Wdir;
dataout(ii).power = power;
dataout(ii).osWdir = oswaves.Wdir;
dataout(ii).osHsig = oswaves.Hsig;
% wavedata.ns.data = tmp;
toc
wavedata.os = oswaves;
wavedata.ns.model = 'MIKE21';
wavedata.ns.contour = [num2str(contour) 'm'];
wavedata.ns.location = profiledataloc;
if savebool
save([saveloc 'wavedata_' lower(site) '_' storm{ii} '.mat'],'wavedata')
end
clear wavedata oswaves startdate enddate validentries
clear power s Hsig Wdir
end

99
example/input_data/fielddata/wave/code/run_OS_to_NS_mike21_other.m
Unescape Escape View File

@ -0,0 +1,99 @@
% Written by Joshua Simmons 04/01/2016
% Description: used to tranform offshore to nearshore mike21 using the code
% OS2NS_MIKE21.m
clear
clc
%% Settings
%buoy data
buoydataloc = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed_CSIROGapFill.mat';
% the non CSIRO filtered data is the only data available after 2014 so fill
% from the recent data with gaps.
backupbuoydata = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed.mat';
%location of mike21 runcode
mike21loc = 'J:\Coastal\Tools\Nearshore Wave Transformations\MIKE21\';
%location of MALT
maltloc = 'J:\Coastal\Tools\MALT Logspiral Transformation\';
%location of profiles for which wave data is required.
% profiledataloc = [6275797.660794820 6275692.739810870 6275536.815570840]; % Northings co-ordinate - Bilgola (at 10m contour)
% profiledataloc = [6272148.185 6272338.391 6272389.496 6271924.185 6271823.724]; % Northings co-ordinate - Mona Vale (at 10m contour)
profiledataloc = [6264042.586 6264127.494 6263995.016 6263745.956 6263879.058]; % Northings co-ordinate - Dee Why (at 10m contour)
site = 'DeeWhy'; %embayment of interest
contour = 10; %m - 10m or 15m contour?
%storm details
stormdatesloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\stormdates_other.mat';
storm = {'2014', '2015'};
% save?
savebool = 1;
% datasave location
saveloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\wave\processed\';
%% Code
% for MIKE21 code and MALT
addpath(mike21loc);
addpath(maltloc);
%load buoy data,profile data and lookup table
buoydata = load(buoydataloc);
backupbuoydata = load(backupbuoydata);
load('OEH_ARC_PointOutput.mat');
stormdates = load(stormdatesloc);
% narrow down to date range
for ii = 1:length(storm)
startdate = stormdates.(['storm_' storm{ii}]).startdate;
enddate = stormdates.(['storm_' storm{ii}]).enddate;
validentries = find(buoydata.dates>=startdate & buoydata.dates<=enddate);
% if nothing then revert to backup
if isempty(validentries) || validentries(end) == length(buoydata.dates)
validentries = find(backupbuoydata.dates>=startdate & backupbuoydata.dates<=enddate);
oswaves.Hsig = backupbuoydata.Hsig(validentries);
oswaves.Tp1 = backupbuoydata.Tp(validentries);
oswaves.Wdir = backupbuoydata.Wdir(validentries);
oswaves.datenums = backupbuoydata.dates(validentries);
else
oswaves.Hsig = buoydata.Hsig(validentries);
oswaves.Tp1 = buoydata.Tp(validentries);
oswaves.Wdir = buoydata.Wdir(validentries);
oswaves.datenums = buoydata.dates(validentries);
end
fns1 = fieldnames(oswaves);
for kk = 1:length(fns1)
if any(isnan(oswaves.(fns1{kk})))
error('NaN values in wave data')
end
end
tic
for jj = 1:length(profiledataloc)
location = profiledataloc(jj);
nswaves = OS2NS_MIKE21(oswaves, LookUps, site, contour, location);
% bug fix to ensure no inaccurate wave period estimates
nswaves.Tp1 = oswaves.Tp1;
tmp(jj) = nswaves;
end
wavedata.ns.data = tmp;
toc
wavedata.os = oswaves;
wavedata.ns.model = 'MIKE21';
wavedata.ns.contour = [num2str(contour) 'm'];
wavedata.ns.location = profiledataloc;
if savebool
save([saveloc 'wavedata_' lower(site) '_' storm{ii} '.mat'],'wavedata')
end
clear wavedata oswaves startdate enddate validentries tmp
end

BIN
example/input_data/fielddata/wave/processed/wavedata_narrabeen_2015.mat
View File

Binary file not shown.

BIN
example/input_data/fielddata/wave/~$2007_survey_to_survey_wave_data.xlsx
View File

Binary file not shown.

BIN
example/input_data/morphology/processed/morphology_narrabeen_2015.mat
View File

Binary file not shown.

BIN
example/setup_info/profileinfo_narrabeen.mat
View File

Binary file not shown.

21
formatpfdata.m
Unescape Escape View File

@ -0,0 +1,21 @@
function [pfdataout] = formatpfdata(xin,zin,interpspace)
% Written by Joshua Simmons 11/02/2015 - to format profile data as
% required (linear interploation etc etc
% Input:
% xin = cross shore
% zin = elevation
% interppspace = regular spacing of the profile to be linearly
% interpolated to. Leave as "NaN" if raw input points to
% be used.
% Output:
% structure with .x,.z formatted as required.
if isnan(interpspace)
pfdataout.x = xin;
pfdataout.z = zin;
else
minx = min(xin);
maxx = max(xin);
pfdataout.x = [minx:interpspace:maxx];
pfdataout.z = interp1(xin,zin,pfdataout.x);
end
end

36
fsbeach_checktsdata.m
Unescape Escape View File

@ -0,0 +1,36 @@
% SBEACH toolbox - fsbeach_checktsdata
% Written by Joshua Simmons 11/2017
% returns all the SBEACH default values
% Syntax:
% OPT = fsbeach_checktsdata()
% Input:
% waves = wave data structure - output by fsbeach_genwaves
% tides = tide data structure - output by fsbeach_gentides
function fsbeach_checktsdata(waves,tide)
% check the tide and wave data matches
if tide.time(1) ~= waves.time(1)
error('Start times for wave and tide data dont match')
end
if tide.time(end) ~= waves.time(end)
error('End times for wave and tide data dont match. Make sure start/end times are hourly (for wave data).')
end
%check and allocate time parameters;
tt = round((tide.time - tide.time(1))*24*60);
%check time steps
dtt = diff(tt);
if any(dtt~=dtt(1))
error('Inconsistent time steps in tide data');
end
%now check the wave data lines up (os only required the rest if copied)
wt = round((waves.time - waves.time(1))*24*60);
dwt = diff(wt);
if any(dwt~=dwt(1))
error('Inconsistent time steps in wave data');
end
end

34
fsbeach_collectresults.m
Unescape Escape View File

@ -0,0 +1,34 @@
% SBEACH toolbox - fsbeach_collectresults
% Written by Joshua Simmons 11/2017
% Syntax:
% sbout = fsbeach_collectresults(sblocs)
%
% Input:
% sblocs = Cell array of locations where SBEACH has been run
% Output:
% sbout = structure with sbeach data;
function sbout = fsbeach_collectresults(sblocs)
for ii = 1:length(sblocs)
fid = fopen([sblocs{ii} '.Excel.ASSUPPLIED'],'r');
formatSpec = '%f %f %f';
sizeA = [3 Inf];
tempresults = fscanf(fid,formatSpec,sizeA)';
fclose(fid);
clear fid
%get computed profile
sbout.computed(ii).z(:,1) = tempresults(:,3);
sbout.computed(ii).x = tempresults(:,1);
%input bathy
sbout.input(ii).x = tempresults(:,1);
sbout.input(ii).z = tempresults(:,2);
%now get wave elevations etc
fid = fopen([sblocs{ii} '.XVR'],'r');
sbout.xvr(ii) = fsbeach_readxvrfile(fid);
end
end

103
fsbeach_genbathy.m
Unescape Escape View File

@ -0,0 +1,103 @@
% SBEACH toolbox - fsbeach_genbathy
% Written by Joshua Simmons 11/2017
% Syntax:
% bathyout = fsbeach_genbathy(varargin)
%
% Input:
% varargin = x: x-coordinates of bathymetry
% z: z-coordinates of bathymetry
% hbz: z-coordinates of hard bottom
% ne: vector or matrix of the size of z containing
% either booleans indicating if a cell is
% non-erodable or a numeric value indicating the
% thickness of the erodable layer on top of a
% non-erodable layer
%
% Output:
% bathyout = stucture of x and z
function bathyout = fsbeach_genbathy(outpath,name,settings,varargin)
OPT = struct( ...
'x', [0 100 200 300 400 500]', ...
'z', [-10 -5 0 1 3 10]', ...
'hbz', NaN, ...
'ne', [] ...
);
OPT = setproperty(OPT, varargin{:});
%transpose them if needed
OPT = fsbeach_makecolumn(OPT,{'x','z','ne'});
%check they can be zipped
if length(OPT.x) ~= length(OPT.z)
error('x is of length %d, while z is of length %d...',length(OPT.x),length(OPT.z))
end
xin = OPT.x; %store x input
%check that constant grid is provided as required by SBEACH
diffx = diff(OPT.x);
dx = diffx(1);
if any(diffx~=dx)
dx = mode(diffx);
fprintf('\n')
warning('SBEACH requires constant spacing - interpolating to a regular grid with spacing of %d...',dx)
newx = [OPT.x(1):dx:OPT.x(end)]';
newz = interp1(OPT.x,OPT.z,newx);
OPT.x = newx;
OPT.z = newz;
clear newx newz
end
if settings.NGRID
tempdx = [0];
for ii = 1:settings.NGRID
tempdx = [tempdx; ones(settings.NDXV(ii),1).*settings.DXV(ii)];
end
newx = OPT.x(1) + cumsum(tempdx);
if newx(end) ~= OPT.x(end)
error(['Variable grid is wrong, please adjust DXV and NDXV so that'...
'xstart = %d and xend = %d. Currently xend = %d.'],newx(1),OPT.x(end),newx(end))
end
newz = interp1(OPT.x,OPT.z,newx);
OPT.x = newx;
OPT.z = newz;
clear newx newz
end
if ~isnan(OPT.hbz)
OPT.hbz = interp1(xin,OPT.hbz,OPT.x);
end
%check they are in the right direction
if mean(diff(OPT.x)) < 0
OPT.x = flip(OPT.x);
OPT.z = flip(OPT.z);
if ~isnan(OPT.hbz)
OPT.hbz = flip(OPT.hbz);
end
end
if mean(diff(OPT.z)) > 0
fprintf('\n')
error('Please ensure x coordinate system is increasing offshore')
end
%zip
prifile = [OPT.x OPT.z]; %initial pf
%and write
writetofilepri(prifile, outpath, name)
if ~isnan(OPT.hbz)
hdbfile = [OPT.x OPT.hbz]; %hard bottom pf
writetofilehdb(hdbfile, outpath, name);
end
bathyout.x = OPT.x;
bathyout.z = OPT.z;
bathyout.hbz = OPT.hbz;
end

33
fsbeach_gentide.m
Unescape Escape View File

@ -0,0 +1,33 @@
% SBEACH toolbox - fsbeach_gentide
% Written by Joshua Simmons 11/2017
% Syntax:
% tideout = fsbeach_gentide(varargin)
%
% Input:
% varargin = time: array of measurement times for tide (minutes)
% elev: array of waterlevels at seaward model border (m)
%
% Output:
% tideout = stucture of time and elevation
function tideout = fsbeach_gentide(outpath,name,varargin)
OPT = struct( ...
'time', [0 100 200 300 400 500]', ...
'elev', [1 1.2 1.1 1.5 1]' ...
);
OPT = setproperty(OPT, varargin{:});
%transpose them if needed
OPT = fsbeach_makecolumn(OPT,{'time','elev'});
%zip
elvfile = OPT.elev; %initial pf
%and write
writetofileelv(elvfile, outpath, name)
tideout.time = OPT.time;
tideout.elev = OPT.elev;
end

62
fsbeach_genwaves.m
Unescape Escape View File

@ -0,0 +1,62 @@
% SBEACH toolbox - fsbeach_genwaves
% Written by Joshua Simmons 11/2017
% Syntax:
% wavesout = fsbeach_genwaves(varargin)
%
% Input:
% varargin = time: array of measurement times for w (minutes)
%
% Hsig: significant wave height (m)
% Tp1: peak wave period (s)
% ang: wave direction (degrees from True North)
% pfang: Profile orientation - from the shore to the
% most seaward point (degrees from True North)
%
% Output:
% wavesout = structure array of wave inputs
function wavesout = fsbeach_genwaves(outpath,name,varargin)
OPT = struct( ...
'time', [0 100 200 300 400 500]', ...
'Hsig', [0 100 200 300 400 500]', ...
'Tp1', [-10 -5 0 1 3 10]', ...
'ang', [-10 -5 0 1 3 10]', ...
'pfang', NaN ...
);
OPT = setproperty(OPT, varargin{:});
%transpose them if needed
OPT = fsbeach_makecolumn(OPT,{'time','Hsig','Tp1','ang'});
%check they can be zipped
if length(OPT.Hsig) ~= length(OPT.Tp1)
error('Hsig is of length %d, while Tp1 is of length %d...',length(OPT.Hsig),length(OPT.Tp1))
end
%zip wave data
wavfile = [OPT.Hsig OPT.Tp1];
% correct for SBEACH wave angle convention
% See SBEACH Report 3 - User's manual (pg 19)
if isnan(OPT.pfang)
error('Please provide the orientation of your profile as pfang.')
end
angfile = OPT.pfang - OPT.ang;
%check
if any(abs(angfile)>90)
warning('Check the orientation of your profile relative to the waves.')
end
%and write
writetofilewav(wavfile, outpath, name)
writetofileang(angfile, outpath, name)
wavesout.time = OPT.time;
wavesout.Hsig = OPT.Hsig;
wavesout.Tp1 = OPT.Tp1;
wavesout.ang = angfile;
end

109
fsbeach_getdefaults.m
Unescape Escape View File

@ -0,0 +1,109 @@
% SBEACH toolbox - fsbeach_getdefaults
% Written by Joshua Simmons 11/2017
% returns all the SBEACH default values
% Syntax:
% OPT = fsbeach_getdefaults()
function OPT = fsbeach_getdefaults()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A - MODEL SETUP
OPT.TITLE = 'Test';
OPT.UNITS = 1;
OPT.NDX = NaN;
OPT.XSTART = NaN;
OPT.IDX = 0;
OPT.DXC = 1;
OPT.NGRID = 0;
OPT.DXV = 0;
OPT.NDXV = 0;
OPT.NDT = NaN;
OPT.DT = 5;
OPT.NWR = 0;
OPT.WRI = 0;
OPT.ICOMP = 0;
OPT.ELV1 = 0;
OPT.ELV2 = 0.25;
OPT.ELV3 = 0.5;
OPT.EDP1 = 0.5;
OPT.EDP2 = 1;
OPT.EDP3 = 1.5;
OPT.REFELV = 0;
% Calibration parameters
OPT.K = 1.75e-06;
OPT.EPS = 0.002;
OPT.LAMM = 0.5;
OPT.TEMPC = 20;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% B - WAVES/WATER ELEVATION/WIND
OPT.WVTYPE = 2;
OPT.IWAVE = 1;
OPT.HIN = 0;
OPT.T = 0;
OPT.DTWAV = NaN;
OPT.IANG = 1;
OPT.ZIN = 0;
OPT.DTANG = NaN;
OPT.DMEAS = NaN; %0 is deepwater
OPT.IRAND = 0;
OPT.ISEED = 4567;
OPT.RPERC = 5;
OPT.IELEV = 1;
OPT.TELEV = 0;
OPT.DTELV = NaN;
OPT.IWIND = 0;
OPT.W = 0;
OPT.ZWIND = 0;
OPT.DTWND = 60;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% C - BEACH
OPT.TPIN = 1;
OPT.SCHEM = 1;
OPT.XLAND = 0;
OPT.DLAND = [];
OPT.XLBDUNE = [];
OPT.DLBDUNE = [];
OPT.XLCDUNE = [];
OPT.DLCDUNE = [];
OPT.XSCDUNE = [];
OPT.DSCDUNE = [];
OPT.XBERMS = [];
OPT.DBERMS = [];
OPT.XBERME = [];
OPT.DBERME = [];
OPT.XFORS = [];
OPT.DFORS = [];
%Calibration params
OPT.DFS = 0.5;
OPT.D50 = 0.35;
OPT.BMAX = 17.2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% D - BEACH FILL
OPT.IBCHFILL = 0;
OPT.XBFS = 0;
OPT.XBFE = 0;
OPT.NFILL = 0;
OPT.XF = [];
OPT.EFILL = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% E - SEAWALL/REVETMENT
OPT.ISWALL = 0;
OPT.XSWALL = 0;
OPT.ISWFAIL = 0;
OPT.PEFAIL = 0;
OPT.WEFAIL = 0;
OPT.HFAIL = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% F - HARD BOTTOM
OPT.IHBOT = 0;
OPT.SCF = 1;
end

55
fsbeach_getstdsettings.m
Unescape Escape View File

@ -0,0 +1,55 @@
% SBEACH toolbox - fsbeach_getstdsettings
% Written by Joshua Simmons 11/2017
% returns all the SBEACH default values
% Syntax:
% OPT = fsbeach_getstdsettings()
% Input:
% bathy = bathy data structure - output by fsbeach_genbathy
% waves = wave data structure - output by fsbeach_genwaves
% tides = tide data structure - output by fsbeach_gentides
%
% Output:
% OPT = settings
function OPT = fsbeach_getstdsettings(bathy, waves, tide, OPT)
% preliminary work
diffx = diff(bathy.x);
dx = diffx(1);
wavestt = (waves.time(end)-waves.time(1))*24*60;
dwt = round((waves.time(2)-waves.time(1))*24*60);
if isstruct(tide)
tidett = (tide.time(end)-tide.time(1))*24*60; %total time
dtt = round((tide.time(2)-tide.time(1))*24*60); %time steps
end
% Section A
OPT.DXC = dx;
OPT.NDX = length(bathy.x); %number of cells
if isnan(OPT.XSTART)
OPT.XSTART = bathy.x(1); %landward calculation boundary
else
warning('XSTART is %d m from first profile point',OPT.XSTART - bathy.x(1))
end
OPT.NDT = round(wavestt/OPT.DT);
% Section B
OPT.DTWAV = dwt;
OPT.DTANG = dwt;
if isnan(OPT.DMEAS)
OPT.DMEAS = -min(bathy.z); %assume it is input at the boundary if not specified
end
if isstruct(tide)
OPT.DTELV = dtt;
end
% Section F - check for hard bottom
if ~isnan(bathy.hbz)
OPT.IHBOT = 1;
end
end

18
fsbeach_makecolumn.m
Unescape Escape View File

@ -0,0 +1,18 @@
% SBEACH toolbox - fsbeach_makecolumn
% Written by Joshua Simmons 11/2017
% convert to column vectors
% Syntax:
% OPT = fsbeach_makecolumn(OPT,vars)
%
% Input:
% OPT: toolbox OPT structure
% vars: variable names to be checked in structure
function OPT = fsbeach_makecolumn(OPT,vars)
for ii = 1:length(vars)
if ~iscolumn(OPT.(vars{ii}))
OPT.(vars{ii}) = OPT.(vars{ii})';
end
end
end

211
fsbeach_modelskill.m
Unescape Escape View File

@ -0,0 +1,211 @@
function [bss, rmse, dVc, dVm, ddV] = fsbeach_modelskill(measured, computed, initial)
%edited by Joshua Simmons 04/04/2016 to revise definitions of r2 etc
%XB_SKILL Computes a variety skill scores
%
% Computes a variety skill scores: R^2, Sci, Relative bias, Brier Skill
% Score. Special feature: within the XBeach testbed, the results are
% stored to be able to show the development of the different skill scores
% in time.
%
% Syntax:
% [r2 sci relbias bss] = xb_skill(measured, computed, varargin)
%
% Input:
% measured = Measured data where the first column contains independent
% values and the second column contains dependent values
% computed = Computed data where the first column contains independent
% values and the second column contains dependent values
% initial = Initial data where the first column contains independent
% values and the second column contains dependent values
% varargin = var: Name of the variable that is supplied
%
% Output:
% bss = Brier Skill Score
% rb20 = Regression through origin R2 value
%
% Example
% [r2 sci relbias bss] = xb_skill(measured, computed)
% [r2 sci relbias bss] = xb_skill(measured, computed, 'var', 'zb')
%
% See also xb_plot_skill
%% Copyright notice
% --------------------------------------------------------------------
% Copyright (C) 2011 Deltares
% Bas Hoonhout
%
% bas.hoonhout@deltares.nl
%
% Rotterdamseweg 185
% 2629HD Delft
%
% This library is free software: you can redistribute it and/or
% modify it under the terms of the GNU Lesser General Public
% License as published by the Free Software Foundation, either
% version 2.1 of the License, or (at your option) any later version.
%
% This library is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% Lesser General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public
% License along with this library. If not, see <http://www.gnu.org/licenses/>.
% --------------------------------------------------------------------
% This tool is part of <a href="http://OpenEarth.nl">OpenEarthTools</a>.
% OpenEarthTools is an online collaboration to share and manage data and
% programming tools in an open source, version controlled environment.
% Sign up to recieve regular updates of this function, and to contribute
% your own tools.
%% Version <http://svnbook.red-bean.com/en/1.5/svn.advanced.props.special.keywords.html>
% Created: 13 Apr 2011
% Created with Matlab version: 7.9.0.529 (R2009b)
% $Id: xb_skill.m 4725 2011-06-27 12:58:48Z hoonhout $
% $Date: 2011-06-27 14:58:48 +0200 (ma, 27 jun 2011) $
% $Author: hoonhout $
% $Revision: 4725 $
% $HeadURL: https://repos.deltares.nl/repos/OpenEarthTools/trunk/matlab/applications/xbeach/xb_analysis/xb_skill.m $
% $Keywords: $
%% read options
if ~exist('initial', 'var'); initial = []; end;
%% xbeach bugfix
% bug - sometime xbeach has suprious points on the end (a flat end to
% the profile) this is fake so remove it.
if length(computed)>=10
if nanmean(computed(1:10,2)) < nanmean(computed(end-9:end,2))
%xbeach direction
spr_len = length(diff(computed(:,2))) - find(diff(computed(:,2))~=0,1,'last');
computed(end-spr_len+1:end,2) = NaN;
else
spr_len = find(diff(computed(:,2))~=0,1,'first');
computed(1:spr_len-1,2) = NaN;
end
if spr_len > 2
disp(['Warning check computed points, ' num2str(spr_len) ' points found to be spurious'] )
end
end
%% remove nans
measured = measured(~any(isnan(measured), 2),:);
computed = computed(~any(isnan(computed), 2),:);
if ~isempty(initial)
initial = initial (~any(isnan(initial ), 2),:);
end
%% compute skills
[bss rmse dVc dVm ddV] = deal(nan);
if size(measured,1) > 1 && size(computed,1) > 1
x = unique([computed(:,1) ; measured(:,1)]);
x = x(~isnan(x));
zmt = interp1(measured(:,1), measured(:,2), x);
zct = interp1(computed(:,1), computed(:,2), x);
% determine active zone
if ~isempty(initial)
zit = interp1(initial(:,1), initial(:,2), x);
%average over 10 steps for small spikes
dzmi = smooth(abs(zmt-zit),10);
dzmi((isnan(zmt)|isnan(zit))) = 0;
dzcm = smooth(abs(zct-zmt),10);
dzcm((isnan(zct)|isnan(zmt))) = 0;
dz = max(dzmi,dzcm);
%has to be more change than a standard measurement error (10cm say)
zi1 = find(dz>0.1,1,'first');
zi2 = find(dz>0.1,1,'last');
x = x (zi1:zi2);
zmt = zmt(zi1:zi2);
zct = zct(zi1:zi2);
zit = zit(zi1:zi2);
end
%% clear all nan values in input (to avoid problems with interpolation)
%all equal at this point and assign values for BSS calc
x0 = x;
xc = x;
xm = x;
z0 = zit;
zc = zct;
zm = zmt;
x0(isnan(z0)) = [];
z0(isnan(z0)) = [];
xc(isnan(zc)) = [];
zc(isnan(zc)) = [];
xm(isnan(zm)) = [];
zm(isnan(zm)) = [];
%% determine new grid covered by all profiles
x_new = unique([x0' xc' xm']);
%% interpolate profiles onto new grid
[x0,id1] = unique(x0);
[xc,id2] = unique(xc);
[xm,id3] = unique(xm);
z0 = z0(id1);
zc = zc(id2);
zm = zm(id3);
z0_new = interp1(x0, z0, x_new);
zc_new = interp1(xc, zc, x_new);
zm_new = interp1(xm, zm, x_new);
% reduce vectors if any nans are left
id = ~(isnan(z0_new) | isnan(zc_new) | isnan(zm_new));
[x_new z0_new zc_new zm_new] = deal(x_new(id), z0_new(id), zc_new(id), zm_new(id));
%% Calculate volume change
try
Vc = volumefind(x_new,zc_new,0);
Vm = volumefind(x_new,zm_new,0);
V0 = volumefind(x_new,z0_new,0);
dVm = Vm-V0; %measured volume change
dVc = Vc-V0; %computed volume change
ddV = dVc-dVm; %predicted volume change - measured
nddV = ddV/abs(V0-Vm);
catch
dVm = NaN;
dVc = NaN;
ddV = NaN;
nddV = NaN;
end
%% derive weight
% weight is defined as half of diff left and half of diff right of each
% point; first and last point are considered to have only one side.
weight = diff([x_new(1) diff(x_new)/2 + x_new(1:end-1) x_new(end)]);
weight = weight / sum(weight);
%% calculate BSS only
% brier skill score
if ~isempty(initial)
mse_p = sum( (zm_new - zc_new).^2 .* weight );
mse_0 = sum( (zm_new - z0_new).^2 .* weight );
bss = 1. - (mse_p/mse_0);
else
bss = 1-(std(zc-zm))^2/var(zm);
end
%% Calculate others
rmse = sqrt(mean((zc_new-zm_new).^2));
rmsm = sqrt(mean(zm_new.^2));
sci = rmse/max(rmsm,abs(mean(zm_new)));
relbias = mean(zc_new-zm_new)/max(rmsm,abs(mean(zm_new)));
end

71
fsbeach_processresults.m
Unescape Escape View File

@ -0,0 +1,71 @@
% SBEACH toolbox - fsbeach_processresults
% Written by Joshua Simmons 11/2017
% Syntax:
% sbout = fsbeach_processresults(sblocs)
%
% Input:
% sbout = structure with input, computed and measured
% Output:
% sbout = structure with sbeach data and statistics;
function sbout = fsbeach_processresults(sbout)
for ii = 1:length(sbout.input)
xin = sbout.input(ii).x;
zin = sbout.input(ii).z;
xcomp = sbout.computed(ii).x;
zcomp = sbout.computed(ii).z;
xmeas = sbout.measured(ii).x;
zmeas = sbout.measured(ii).z;
try
[tmpbss, tmprmse, dVc, dVm, ddV] = fsbeach_modelskill([xmeas zmeas], [xcomp zcomp], [xin zin]);
%calculate BSS above MSL
%first get points above MSL - for measured and then apply that
%to all - cant take these independently as may only get say one
%point for modelled and then xb_skill interpolates.
%BSS deals well if based on measured as the -10 and x
%determinant - it only finds the interpolated values with
%info so will deal if initpf doesnt stretch to all of
%measured pf
keepme1 = find(zmeas>0);
if xmeas(keepme1(1))>xmeas(keepme1(end))
keepmo1 = find(xcomp>xmeas(keepme1(end)) & xcomp<xmeas(keepme1(1)));
keepin1 = find(xin>xmeas(keepme1(end)) & xin<xmeas(keepme1(1)));
elseif xmeas(keepme1(1))<xmeas(keepme1(end))
keepmo1 = find(xcomp<xmeas(keepme1(end)) & xcomp>xmeas(keepme1(1)));
keepin1 = find(xin<xmeas(keepme1(end)) & xin>xmeas(keepme1(1)));
end
[tmpabMSLbss, tmpabMSLrmse, ~, ~, ~] = fsbeach_modelskill([xmeas(keepme1) zmeas(keepme1)], [xcomp(keepmo1) zcomp(keepmo1)], [xin(keepin1) zin(keepin1)]);
%calculate BSS below MSL
keepme2 = find(zmeas<0);
if xmeas(keepme2(1))>xmeas(keepme2(end))
keepmo2 = find(xcomp>xmeas(keepme2(end)) & xcomp<xmeas(keepme2(1)));
keepin2 = find(xin>xmeas(keepme2(end)) & xin<xmeas(keepme2(1)));
elseif xmeas(keepme2(1))<xmeas(keepme2(end))
keepmo2 = find(xcomp<xmeas(keepme2(end)) & xcomp>xmeas(keepme2(1)));
keepin2 = find(xin<xmeas(keepme2(end)) & xin>xmeas(keepme2(1)));
end
[tmpbeMSLbss,tmpbeMSLrmse, ~, ~, ~] = fsbeach_modelskill([xmeas(keepme2) zmeas(keepme2)], [xcomp(keepmo2) zcomp(keepmo2)], [xin(keepin2) zin(keepin2)]);
%store results
sbout.stats(ii).bss = tmpbss;
sbout.stats(ii).rmse = tmprmse;
sbout.stats(ii).dVc = dVc;
sbout.stats(ii).dVm = dVm;
sbout.stats(ii).ddV = ddV;
sbout.stats(ii).aboveMSL.bss = tmpabMSLbss;
sbout.stats(ii).aboveMSL.rmse = tmpabMSLrmse;
sbout.stats(ii).belowMSL.bss = tmpbeMSLbss;
sbout.stats(ii).belowMSL.rmse = tmpbeMSLrmse;
catch
disp(['Run number ' num2str(ii) ' failed.'])
continue
end
end
end

43
fsbeach_readxvrfile.m
Unescape Escape View File

@ -0,0 +1,43 @@
% SBEACH toolbox - fsbeach_collectresults
% Written by Joshua Simmons 02/2018
% Syntax:
% dataout = fsbeach_readxvrfile(fid)
%
% Input:
% fid = file id of .XVR file from which data is to be extracted
% Output:
% dataout = structure with sbeach data;
function dataout = fsbeach_readxvrfile(fid)
ln = fgets(fid); % first line
%header material
while ln(1) ~= '_'
ln = fgets(fid);
end
for ii = 1:8
tmpd = [];
while any(isletter(ln)) || ln(1) == '_'
ln = fgets(fid);
end
while ln(1) ~= '_' && ~any(isletter(ln)) && ischar(ln)
tmpd = [tmpd; sscanf(ln,'%f')];
ln = fgets(fid);
end
tmpdataout(ii).data = tmpd;
end
padlen = length(tmpdataout(1).data)-length(tmpdataout(2).data);
%max value
dataout.Hmax = padarray(tmpdataout(2).data,padlen,0,'pre');
dataout.TWLmax = padarray(tmpdataout(4).data,padlen,0,'pre');
dataout.depmax = padarray(tmpdataout(6).data,padlen,0,'pre');
%time step of max
dataout.tsHmax = padarray(tmpdataout(3).data,padlen,0,'pre');
dataout.tsTWLmax = padarray(tmpdataout(5).data,padlen,0,'pre');
dataout.tsdepmax = padarray(tmpdataout(7).data,padlen,0,'pre');
end

18
fsbeach_runmodel.m
Unescape Escape View File

@ -0,0 +1,18 @@
% SBEACH toolbox - fsbeach_runmodel
% Written by Joshua Simmons 11/2017
% Syntax:
% success = fsbeach_runmodel(cfgpath,exepath)
%
% Input:
% cfgpath = path to SBEACH .CFG file (and other files setup by
% fsbeach_setupmodel()
% exepath = path to SBEACH exe
% Output:
% success = bool;
function [fail, cmdout] = fsbeach_runmodel(cfgpath,exepath)
cmdin = ['"' exepath '"' ' -cfgFN ' '"' cfgpath '"'];
[fail, cmdout] = system(cmdin);
end

6
fsbeach_runsetup.asv
Unescape Escape View File

@ -0,0 +1,6 @@
% SBEACH toolbox
% Written by Joshua Simmons
function out = fsbeach_runsetup()
end

6
fsbeach_runsetup.m
Unescape Escape View File

@ -0,0 +1,6 @@
% SBEACH toolbox
% Written by Joshua Simmons
function out = fsbeach_runsetup()
end

96
fsbeach_setupmodel.m
Unescape Escape View File

@ -0,0 +1,96 @@
% SBEACH toolbox - fsbeach_setupmodel
% Written by Joshua Simmons 11/2017
% Note this code is specifically written to align with (and borrows code
% from) the XBeach matlab toolbox so as to make the use of both models
% streamlined.
% Syntax:
% OPT = fsbeach_setupmodel(varargin)
%
% Input:
% varargin = bathy: cell array of name/value pairs of bathymetry
% settings supplied to fsbeach_genbathy
% waves: cell array of name/value pairs of waves
% settings supplied to fsbeach_genwaves
% tide: cell array of name/value pairs of tide
% settings supplied to fsbeach_gentide
% settings: cell array of name/value pairs of model
% settings supplied to fsbeach_writecfg
% path: destination directory of model setup, if
% written to disk
%
% Output:
% OPT = options
function out = fsbeach_setupmodel(varargin)
OPT = struct( ...
'bathy', {{}}, ...
'waves', {{}}, ...
'tide', {{}}, ...
'settings', {{}}, ...
'path', [pwd '\temp\'], ...
'name', 'sbeach_test' ...
);
OPT = setproperty(OPT, varargin{:});
%check we are ready to write!
if ~isdir(OPT.path)
mkdir(OPT.path)
else
%delete any old SBEACH files
delete(fullfile(OPT.path,'*.PRI'));
delete(fullfile(OPT.path,'*.WAV'));
delete(fullfile(OPT.path,'*.ANG'));
delete(fullfile(OPT.path,'*.ELV'));
delete(fullfile(OPT.path,'*.HDB'));
delete(fullfile(OPT.path,'*.CFG'));
delete(fullfile(OPT.path,'*.AsSupplied'));
delete(fullfile(OPT.path,'*.LOG'));
delete(fullfile(OPT.path,'*.PRC'));
delete(fullfile(OPT.path,'*.RPTNew'));
delete(fullfile(OPT.path,'*.XVR'));
end
%%% Initial settings
%get defaults
settingsOPT = fsbeach_getdefaults();
%apply any setttings as needed
settingsOPT = setproperty(settingsOPT, OPT.settings);
%%% Bathy
%print to file
bathy = fsbeach_genbathy(OPT.path,OPT.name,settingsOPT,OPT.bathy{:});
%%% Waves
%print to file
waves = fsbeach_genwaves(OPT.path,OPT.name,OPT.waves{:});
%%% Tide
%print to file
if isempty(OPT.tide) && ~any(strcmp(OPT.settings,'IELEV'))
error('Please supply tide or use the IELEV setting to supply constant WL.')
elseif ~isempty(OPT.tide)
tide = fsbeach_gentide(OPT.path,OPT.name,OPT.tide{:});
%check timeseries data
fsbeach_checktsdata(waves,tide);
else
tide = NaN; % no tide supplied
end
%%% Final settings
% get standard settings from bathy, waves and tides
settingsOPT = fsbeach_getstdsettings(bathy, waves, tide, settingsOPT);
%write the configuration file
fsbeach_writecfg(OPT.path,OPT.name,settingsOPT);
%%% Output
%structure out
out.bathy = bathy;
out.waves = waves;
out.tides = tide;
out.settings = settingsOPT;
out.path = [OPT.path OPT.name];
end

18
fsbeach_writecfg.m
Unescape Escape View File

@ -0,0 +1,18 @@
% SBEACH toolbox - fsbeach_writecfg
% Written by Joshua Simmons 11/2017
% Syntax:
% wavesout = fsbeach_writecfg(outpath,standard_settings,varargin)
%
% Input:
% outpath: output path for CFG file
% standard_settings: standard settings for bathy, waves, tide etc.
%
% varargin = sbeach settings
%
% Output:
% cfg = all config settings
function fsbeach_writecfg(outpath,name,OPT)
writetofilecfg(OPT, outpath, name);
end

BIN
sbeach_exe/final_exe/SBEACH.exe
View File

Binary file not shown.

471
setproperty.m
Unescape Escape View File

@ -0,0 +1,471 @@
function [OPT Set Default] = setproperty(OPT, inputCell, varargin)
% SETPROPERTY generic routine to set values in PropertyName-PropertyValue pairs
%
% Routine to set properties based on PropertyName-PropertyValue
% pairs (aka <keyword,value> pairs). Can be used in any function
% where PropertyName-PropertyValue pairs are used.
%
%
% [OPT Set Default] = setproperty(OPT, inputCell, varargin)
%
% input:
% OPT = structure in which fieldnames are the keywords and the values
% are the defaults
% inputCell = must be a cell array containing single struct, or a set of
% 'PropertyName', PropertyValue,... pairs. It is best practice
% to pass the varargin of the calling function as the 2nd argument.
% Property names must be strings. If they contain a dot (.) this is
% interpreted a field separator. This can be used to assign
% properties on a subfield level.
% varargin = series of 'PropertyName', PropertyValue,... pairs to set
% methods for setproperty itself.
%
% output:
% OPT = structure, similar to the input argument OPT, with possibly
% changed values in the fields
% Set = structure, similar to OPT, values are true where OPT has been
% set (and possibly changed)
% Default = structure, similar to OPT, values are true where the values of
% OPT are equal to the original OPT
%
% Example calls:
%
% [OPT Set Default] = setproperty(OPT, vararginOfCallingFunction)
% [OPT Set Default] = setproperty(OPT, {'PropertyName', PropertyValue,...})
% [OPT Set Default] = setproperty(OPT, {OPT2})
%
% Any number of leading structs with
% any number of trailing <keyword,value> pairs:
% [OPT Set Default] = setproperty(OPT1, OPT2, ..., OPTn)
% [OPT Set Default] = setproperty(OPT1, OPT2, ..., OPTn,'PropertyName', PropertyValue,...)
%
% Different methods for dealing with class changes of variables, or
% extra fields (properties that are not in the input structure) can
% be defined as property-value pairs. Valid properties are
% onExtraField and onClassChange:
%
% PROPERTY VALUE
% onClassChange: ignore ignore (default)
% warn throw warning
% error throw error
% onExtraField: silentIgnore silently ignore the field
% warnIgnore ignore the field and throw warning
% silentAppend silently append the field to OPT
% warnAppend append the field to OPT and throw warning
% error throw error (default)
%
% Example calls:
%
% [OPT Set Default] = setproperty(OPT, vararginOfCallingFunction,'onClassChange','warn')
% [OPT Set Default] = setproperty(OPT, {'PropertyName', PropertyValue},'onExtraField','silentIgnore')
%
%
% Example:
%
% +------------------------------------------->
% function y = dosomething(x,'debug',1)
% OPT.debug = 0;
% OPT = setproperty(OPT, varargin);
% y = x.^2;
% if OPT.debug; plot(x,y);pause; end
% +------------------------------------------->
%
% legacy syntax is also supported, but using legacy syntax prohibits the
% setting of onClassChange and onExtraField methods:
%
% [OPT Set Default] = setproperty(OPT, varargin{:})
% OPT = setproperty(OPT, 'PropertyName', PropertyValue,...)
% OPT = setproperty(OPT, OPT2)
%
% input:
% OPT = structure in which fieldnames are the keywords and the values are the defaults
% varargin = series of PropertyName-PropertyValue pairs to set
% OPT2 = is a structure with the same fields as OPT.
%
% Internally setproperty translates OPT2 into a set of
% PropertyName-PropertyValue pairs (see example below) as in:
% OPT2 = struct( 'propertyName1', 1,...
% 'propertyName2', 2);
% varcell = reshape([fieldnames(OPT2)'; struct2cell(OPT2)'], 1, 2*length(fieldnames(OPT2)));
% OPT = setproperty(OPT, varcell{:});
%
% Change log:
% 2011-09-30: full code rewrite to include:
% - setpropertyInDeeperStruct functionality
% - user defined handling of extra fields
% - class change warning/error message
%
% See also: VARARGIN, STRUCT, MERGESTRUCTS
%% Copyright notice
% --------------------------------------------------------------------
% Copyright (C) 2009 Delft University of Technology
% C.(Kees) den Heijer
%
% C.denHeijer@TUDelft.nl
%
% Faculty of Civil Engineering and Geosciences
% P.O. Box 5048
% 2600 GA Delft
% The Netherlands
%
% This library is free software; you can redistribute it and/or
% modify it under the terms of the GNU Lesser General Public
% License as published by the Free Software Foundation; either
% version 2.1 of the License, or (at your option) any later version.
%
% This library is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% Lesser General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public
% License along with this library; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
% USA
% or http://www.gnu.org/licenses/licenses.html, http://www.gnu.org/, http://www.fsf.org/
% --------------------------------------------------------------------
% This tools is part of <a href="http://OpenEarth.Deltares.nl">OpenEarthTools</a>.
% OpenEarthTools is an online collaboration to share and manage data and
% programming tools in an open source, version controlled environment.
% Sign up to recieve regular updates of this function, and to contribute
% your own tools.
%% Version <http://svnbook.red-bean.com/en/1.5/svn.advanced.props.special.keywords.html>
% Created: 26 Feb 2009
% Created with Matlab version: 7.4.0.287 (R2007a)
% $Id: setproperty.m 8036 2013-02-05 21:06:43Z boer_g $
% $Date: 2013-02-05 22:06:43 +0100 (Tue, 05 Feb 2013) $
% $Author: boer_g $
% $Revision: 8036 $
% $HeadURL: https://svn.oss.deltares.nl/repos/openearthtools/trunk/matlab/setproperty.m $
% $Keywords: $
%% shortcut function if there is nothing to set (1)
if nargin==1||isempty(inputCell)
if nargout > 1 % process Set
flds1 = fieldnames(OPT);
Set = [flds1,repmat({false},size(flds1))]';
Set = struct(Set{:});
end
if nargout > 2 % process Default
Default = [flds1,repmat({true},size(flds1))]';
Default = struct(Default{:});
end
return
end
%% check and parse inputCell (usually the varargin struct in the calling function)
% determine mode from class of the second argument
switch(class(inputCell))
case 'struct'
% recursively let setproperty peel all leading structs
% with optional trailing keyword,value> pairs
inputCell = setproperty(inputCell,varargin{:});
inputCell = struct2arg(inputCell);
varargin = {};
case 'char'
% legay syntax mode
inputCell = [{inputCell} varargin];
varargin = {};
case 'cell'
% recursively let setproperty peel all leading structs from cell
if isstruct(inputCell{1})
inputCell = {setproperty(inputCell{:})};
end
otherwise
error('SETPROPERTY:inputCell',...
'Second input must be a cell, (or a char or struct for legacy syntax)')
end
%% shortcut function if there is nothing to set (2)
if numel(inputCell) == 1 && isempty(inputCell{1})
if nargout > 1 % process Set
flds1 = fieldnames(OPT);
Set = [flds1,repmat({false},size(flds1))]';
Set = struct(Set{:});
end
if nargout > 2 % process Default
Default = [flds1,repmat({true},size(flds1))]';
Default = struct(Default{:});
end
return
end
if odd(length(inputCell))
%then length must be 1
if length(inputCell) == 1
%then the inputCell must be a structure
if ~isstruct(inputCell{1})
error('SETPROPERTY:inputCell',...
'Second argument inputCell must be a cell containing a single struct or property/value pairs')
end
flds2 = fieldnames(inputCell{1})';
newVal = struct2cell(inputCell{1})';
else
error('SETPROPERTY:inputCell',...
'Second argument inputCell must be a cell containing a single struct or property/value pairs')
end
else
flds2 = inputCell(1:2:end);
if ~all(cellfun(@ischar,flds2))
error('SETPROPERTY:inputCell',...
'Second argument inputCell, property names should be strings')
end
newVal = inputCell(2:2:end);
end
%% check and parse varargin
if length(varargin) > 4
error('Too many input arguments to setpropoerty.');
end
% error(nargchk(0, 4, length(varargin), 'struct'))
if odd(length(varargin))
error('SETPROPERTY:varargin',...
'Set onClassChange and onExtraField with the varargin, as keyword value pairs');
end
onExtraField = 'error';
onClassChange = 'ignore';
if ~isempty(varargin)
for ii = 1:2:length(varargin)
switch varargin{ii}
case 'onClassChange'
if ischar(varargin{ii+1}) && ismember(varargin(ii+1),...
{'ignore','warn','error'})
onClassChange = varargin{ii+1};
else
error('SETPROPERTY:InvalidMethod',...
'Valid methods for onClassChange are: ignore, warn and error');
end
case 'onExtraField'
if ischar(varargin{ii+1}) && ismember(varargin(ii+1),...
{'silentAppend','warnAppend','silentIgnore','warnIgnore','error'})
onExtraField = varargin{ii+1};
else
error('SETPROPERTY:InvalidMethod',...
'Valid methods for onExtraField are: append, silentIgnore, warnIgnore, error');
end
otherwise
error('SETPROPERTY:InvalidKeyword', ...
'Supported keywords are onClassChange and onExtraField');
end
end
end
switch onClassChange
case 'warn'
warnOnClassChange = true;
errorOnClassChange = false;
case 'error'
warnOnClassChange = false;
errorOnClassChange = true;
case 'ignore'
warnOnClassChange = false;
errorOnClassChange = false;
end
%% copy the original OPT only if Default has to be returned
if nargout > 2
OPToriginal = OPT;
end
%% check field names, and distinguish between field to set and extra fields
% flds1 contains field names of OPT
flds1 = fieldnames(OPT);
% fldsCell contains field names to assign, split per subfield
fldsCell = regexp(flds2,'[^\.]*','match');
fldsCellLen = cellfun(@length,fldsCell);
% check if all field names are valid to assign
allFieldNames = [fldsCell{:}];
validFieldName = cellfun(@(x) (isvarname(x) | iskeyword(x)),allFieldNames); % a field like 'case' is allowed also, but is not a valid varname
if ~all(validFieldName)
error('SETPROPERTY:invalidFieldName',...
'\nThe following field name is not valid: %s',allFieldNames{~validFieldName});
end
% identify fldsToSet for simple field names (without subfields)
% fldsToSet = ismember(flds2,flds1);
% As ismember is rather slow, this is replaced by more optimized code
fldsToSet = false(size(flds2));
for ii = find(fldsCellLen==1)
fldsToSet(ii) = any(strcmp(flds2(ii),flds1));
end
% identify fldsToSet for field names with subfields)
for ii = find(fldsCellLen>1)
fldsToSet(ii) = isfield2(OPT,fldsCell{ii});
end
% all other fields are either to be set, or extra
fldsExtra = ~fldsToSet;
%% process fldsToSet
if any(fldsToSet)
for ii = find(fldsToSet)
switch length(fldsCell{ii})
case 1
if warnOnClassChange || errorOnClassChange
class1 = class(OPT.(fldsCell{ii}{1}));
class2 = class(newVal{ii});
classChange(flds2{ii},class1,class2,warnOnClassChange,errorOnClassChange);
end
OPT.(fldsCell{ii}{1}) = newVal{ii};
case 2
if warnOnClassChange || errorOnClassChange
class1 = class(OPT.(fldsCell{ii}{1}).(fldsCell{ii}{2}));
class2 = class(newVal{ii});
classChange(flds2{ii},class1,class2,warnOnClassChange,errorOnClassChange);
end
OPT.(fldsCell{ii}{1}).(fldsCell{ii}{2}) = newVal{ii};
otherwise
if warnOnClassChange || errorOnClassChange
class1 = class(getfield(OPT,fldsCell{ii}{:}));
class2 = class(newVal{ii});
classChange(flds2{ii},class1,class2,warnOnClassChange,errorOnClassChange);
end
OPT = setfield(OPT,fldsCell{ii}{:},newVal{ii});
end
end
end
%% deal with fldsExtra
if any(fldsExtra)
switch lower(onExtraField)
case {'silentappend','warnappend'}
% append extra fields to OPT
for ii = find(fldsExtra)
switch length(fldsCell{ii})
case 1
OPT.(flds2{ii}) = newVal{ii};
otherwise
OPT = setfield(OPT,fldsCell{ii}{:},newVal{ii});
end
end
if strcmpi(onExtraField,'warnAppend')
% throw a warning
warning('SETPROPERTY:ExtraField',...
'\nThe following field is Extra and appended to OPT: %s',flds2{fldsExtra});
end
case 'silentignore'
% do nothing, silently ignore extra fields
case 'warnignore'
warning('SETPROPERTY:ExtraField',...
'\nThe following field is Extra and thus ignored: %s',flds2{fldsExtra});
case 'error'
error('SETPROPERTY:ExtraField',...
'\nThe following field is Extra: %s',flds2{fldsExtra});
end
end
%% assign Output
if nargout > 1 % process Set
Set = [flds1,repmat({false},size(flds1))]';
Set = struct(Set{:});
if any(fldsToSet)
for ii = find(fldsToSet)
switch length(fldsCell{ii})
case 1
Set.(fldsCell{ii}{1}) = true;
otherwise % this routine is rather slow, would be nice to add case 2 to cover the most common calls
for nn = 1:length(fldsCell{ii})-1
if ~isstruct(getfield(Set,fldsCell{ii}{1:nn})) %#ok<*GFLD>
substruct = [fieldnames(getfield(OPT,fldsCell{ii}{1:nn})),...
repmat({false},size(fieldnames(getfield(OPT,fldsCell{ii}{1:nn}))))]';
Set = setfield(Set,fldsCell{ii}{1:nn},[]); %#ok<*SFLD>
Set = setfield(Set,fldsCell{ii}{1:nn},struct(substruct{:}));
end
end
Set = setfield(Set,fldsCell{ii}{:},true);
end
end
end
if any(fldsExtra)
switch lower(onExtraField)
case {'silentappend','warnappend'}
for ii = find(fldsExtra)
switch length(fldsCell{ii})
case 1
Set.(flds2{ii}) = true;
otherwise
Set = setfield(Set,fldsCell{ii}{:},true);
end
end
end
end
end
if nargout > 2 % process Default
Default = [flds1,repmat({true},size(flds1))]';
Default = struct(Default{:});
if any(fldsToSet)
for ii = find(fldsToSet)
switch length(fldsCell{ii})
case 1
Default.(fldsCell{ii}{1}) = ...
isequalwithequalnans(OPToriginal.(fldsCell{ii}{1}),OPT.(fldsCell{ii}{1}));
otherwise % this routine is rather slow, would be nice to add case 2 to cover the most common calls
for nn = 1:length(fldsCell{ii})-1
if ~isstruct(getfield(Default,fldsCell{ii}{1:nn})) %#ok<*GFLD>
substruct = [fieldnames(getfield(OPT,fldsCell{ii}{1:nn})),...
repmat({true},size(fieldnames(getfield(OPT,fldsCell{ii}{1:nn}))))]';
Default = setfield(Default,fldsCell{ii}{1:nn},[]); %#ok<*SFLD>
Default = setfield(Default,fldsCell{ii}{1:nn},struct(substruct{:}));
end
end
Default = setfield(Default,fldsCell{ii}{:},...
isequalwithequalnans(getfield(OPToriginal,fldsCell{ii}{:}),getfield(OPT,fldsCell{ii}{:})));
end
end
end
if any(fldsExtra)
switch lower(onExtraField)
case {'silentappend','warnappend'}
for ii = find(fldsExtra)
switch length(fldsCell{ii})
case 1
Default.(flds2{ii}) = false;
otherwise
Default = setfield(Default,fldsCell{ii}{:},false);
end
end
end
end
end
function tf = isfield2(OPT,fldsCell)
if isfield(OPT,fldsCell{1})
if length(fldsCell) > 1
tf = isfield2(OPT.(fldsCell{1}),fldsCell(2:end));
else
tf = true;
end
else
tf = false;
end
function classChange(fld,class1,class2,warnOnClassChange,errorOnClassChange)
if ~strcmp(class1,class2)
if warnOnClassChange
warning('SETPROPERTY:ClassChange', ...
'Class change of field ''%s'' from %s to %s',...
fld,class1,class2);
elseif errorOnClassChange
error('SETPROPERTY:ClassChange', ...
'Class change of field ''%s'' from %s to %s not allowed',...
fld,class1,class2);
end
end
function out = odd(in)
%ODD test whether number if odd
out = mod(in,2)==1;
%% EOF

39
volumefind.m
Unescape Escape View File

@ -0,0 +1,39 @@
function [VolumeTotal, Beachwidth] = volumefind(x,z,shore)
% Calculates Volume and Beach Width
% Written by Matthew Phillips 2015
% Define Shoreline Contour and Fixed Landward Column Boundary
ShoreCont = shore; % m AHD
BackCol = 1; % column number
if ~iscolumn(z)
z = z';
end
if ~iscolumn(x)
x = x';
end
if z(end) > z(1)
z = flip(z);
x = flip(x);
end
if x(end) < x(1)
x = -x;
end
xstep = mean(diff(x));
shorecol = find(z>=ShoreCont,1,'last')+1;
% Width at MSL
Beachwidth = interp1(z(shorecol-1:shorecol,1),x(shorecol-1:shorecol,1),ShoreCont);
% Total Volume
Volx = [x(BackCol:shorecol-1); Beachwidth];
Volz = [z(BackCol:shorecol-1); ShoreCont];
VolumeTotal = trapz(Volx,Volz);
end

18
writetofileang.m
Unescape Escape View File

@ -0,0 +1,18 @@
function writetofileang(angfile, outpath, pfstring)
fid = fopen([outpath '\' pfstring '.ANG'],'w');
fprintf(fid,'C> *******************************************************\r\n');
fprintf(fid,'C> * SBEACH wave input file: %s.ANG *\r\n',pfstring);
fprintf(fid,'C> *******************************************************\r\n');
fprintf(fid,'C>\r\n');
fprintf(fid,'C>\r\n');
fprintf(fid,'C> DIRECTION (DEGREES)\r\n');
fprintf(fid,'E>----------------------------------------------------------------\r\n');
for ii = 1:length(angfile)
if ii == length(angfile)
fprintf(fid,'%.3f\r\n', angfile(ii));
else
fprintf(fid,'%.3f\r\n', angfile(ii));
end
end
fclose(fid);
end

173
writetofilecfg.m
Unescape Escape View File

@ -0,0 +1,173 @@
function writetofilecfg(OPT, outpath, pfstring)
fid = fopen([outpath '\' pfstring '.CFG'],'w');
fprintf(fid,' *******************************************************\r\n');
fprintf(fid,' * SBEACH model configuration file: %s.CFG *\r\n',pfstring);
fprintf(fid,' *******************************************************\r\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A - MODEL SETUP
fprintf(fid,'A---------------------------- MODEL SETUP -------------------------------A\r\n');
fprintf(fid,'A.1 RUN TITLE: TITLE\r\n');
fprintf(fid,' %s\r\n',pfstring);
fprintf(fid,'A.2 INPUT UNITS (SI=1, AMERICAN CUST.=2): UNITS\r\n');
fprintf(fid,' %d\r\n',OPT.UNITS);
fprintf(fid,'A.3 TOTAL NUMBER OF CALCULATION CELLS AND POSITION OF LANDWARD BOUNDARY\r\n');
fprintf(fid,' RELATIVE TO INITIAL PROFILE: NDX, XSTART\r\n');
fprintf(fid,' %d %d\r\n',OPT.NDX,OPT.XSTART);
fprintf(fid,'A.4 GRID TYPE (CONSTANT=0, VARIABLE=1): IDX\r\n');
fprintf(fid,' %d\r\n',OPT.IDX);
fprintf(fid,'A.5 COMMENT: IF GRID TYPE IS VARIABLE, CONTINUE TO A.8\r\n');
fprintf(fid,'A.6 CONSTANT GRID CELL WIDTH: DXC\r\n');
fprintf(fid,' %d\r\n',OPT.DXC);
fprintf(fid,'A.7 COMMENT: IF GRID TYPE IS CONSTANT CONTINUE TO A.10\r\n');
fprintf(fid,'A.8 NUMBER OF DIFFERENT GRID CELL REGIONS: NGRID\r\n');
fprintf(fid,' %d\r\n',OPT.NGRID);
fprintf(fid,'A.9 GRID CELL WIDTHS AND NUMBER OF CELLS IN EACH REGION FROM LANDWARD\r\n');
fprintf(fid,' TO SEAWARD BOUNDARY: (DXV(I), NDXV(I), I=1,NGRID)\r\n');
for ii = 1:length(OPT.NGRID)
fprintf(fid,' %d %d',OPT.DXV(ii),OPT.NDXV(ii));
end
fprintf(fid,'\r\n');
fprintf(fid,'A.10 NUMBER OF TIME STEPS AND VALUE OF TIME STEP IN MINUTES: NDT,DT\r\n');
fprintf(fid,' %d %d\r\n',OPT.NDT,OPT.DT);
fprintf(fid,'A.11 NUMBER OF TIME STEP(S) INTERMEDIATE OUTPUT IS WANTED: NWR\r\n');
fprintf(fid,' %d\r\n',OPT.NWR);
fprintf(fid,'A.12 TIME STEPS OF INTERMEDIATE OUTPUT: (WRI(I), I=1,NWR)\r\n');
for ii = 1:length(OPT.NWR)
fprintf(fid,' %d',OPT.WRI(ii));
end
fprintf(fid,'\r\n');
fprintf(fid,'A.13 IS A MEASURED PROFILE AVAILABLE FOR COMPARISON? (NO=0, YES=1): ICOMP\r\n');
fprintf(fid,' %d\r\n',OPT.ICOMP);
fprintf(fid,'A.14 THREE PROFILE ELEVATION CONTOURS (MAXIMUM HORIZONTAL RECESSION OF EACH\r\n');
fprintf(fid,' WILL BE DETERMINED): ELV1, ELV2, ELV3\r\n');
fprintf(fid,' %.2f %.2f %.2f\r\n',OPT.ELV1,OPT.ELV2,OPT.ELV3);
fprintf(fid,'A.15 THREE PROFILE EROSION DEPTHS AND REFERENCE ELEVATION (DISTANCE FROM\r\n');
fprintf(fid,' POSITION OF REFERENCE ELEVATION ON INITIAL PROFILE TO POSITION OF\r\n');
fprintf(fid,' LANDWARD MOST OCCURENCE OF EACH EROSION DEPTH WILL BE DETERMINED\r\n');
fprintf(fid,' EDP1, EDP2, EDP3, REFELV\r\n');
fprintf(fid,' %.2f %.2f %.2f %.2f\r\n',OPT.EDP1,OPT.EDP2,OPT.EDP3,OPT.REFELV);
fprintf(fid,'A.16 TRANSPORT RATE COEFFICIENT (m^4/N): K\r\n');
fprintf(fid,' %.3e\r\n',OPT.K);
fprintf(fid,'A.17 COEFFICIENT FOR SLOPE-DEPENDENT TERM (m^2/s): EPS\r\n');
fprintf(fid,' %.3f\r\n',OPT.EPS);
fprintf(fid,'A.18 TRANSPORT RATE DECAY COEFFICIENT MULTIPLIER: LAMM\r\n');
fprintf(fid,' %.3f\r\n',OPT.LAMM);
fprintf(fid,'A.19 WATER TEMPERATURE IN DEGREES C: TEMPC\r\n');
fprintf(fid,' %.2f\r\n',OPT.TEMPC);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% B - WAVES/WATER ELEVATION/WIND
fprintf(fid,'B---------------------- WAVES/WATER ELEVATION/WIND -----------------------B\r\n');
fprintf(fid,'B.1 WAVE TYPE (MONOCHROMATIC=1, IRREGULAR=2): WVTYPE\r\n');
fprintf(fid,' %d\r\n',OPT.WVTYPE);
fprintf(fid,'B.2 WAVE HEIGHT AND PERIOD INPUT (CONSTANT=0, VARIABLE=1): IWAVE\r\n');
fprintf(fid,' %d\r\n',OPT.IWAVE);
fprintf(fid,'B.3 COMMENT: IF WAVE HEIGHT AND PERIOD ARE VARIABLE, CONTINUE TO B.6\r\n');
fprintf(fid,'B.4 CONSTANT WAVE HEIGHT AND PERIOD: HIN, T\r\n');
fprintf(fid,' %.3f %.3f\r\n',OPT.HIN,OPT.T);
fprintf(fid,'B.5 COMMENT: IF WAVE HEIGHT AND PERIOD ARE CONSTANT, CONTINUE TO B.7\r\n');
fprintf(fid,'B.6 TIME STEP OF VARIABLE WAVE HEIGHT AND PERIOD INPUT IN MINUTES: DTWAV\r\n');
fprintf(fid,' %.2f\r\n',OPT.DTWAV);
fprintf(fid,'B.7 WAVE ANGLE INPUT (CONSTANT=0, VARIABLE=1): IANG\r\n');
fprintf(fid,' %d\r\n',OPT.IANG);
fprintf(fid,'B.8 COMMENT: IF WAVE ANGLE IS VARIABLE, CONTINUE TO B.11\r\n');
fprintf(fid,'B.9 CONSTANT WAVE ANGLE: ZIN\r\n');
fprintf(fid,' %.3f\r\n',OPT.ZIN);
fprintf(fid,'B.10 COMMENT: IF WAVE ANGLE IS CONSTANT, CONTINUE TO B.12\r\n');
fprintf(fid,'B.11 TIME STEP OF VARIABLE WAVE ANGLE INPUT IN MINUTES: DTANG\r\n');
fprintf(fid,' %.2f\r\n',OPT.DTANG);
fprintf(fid,'B.12 WATER DEPTH OF INPUT WAVES (DEEPWATER=0): DMEAS\r\n');
fprintf(fid,' %d\r\n',OPT.DMEAS);
fprintf(fid,'B.13 IS RANDOMIZATION OF WAVE HEIGHT DESIRED? (NO=0, YES=1): IRAND\r\n');
fprintf(fid,' %d\r\n',OPT.IRAND);
fprintf(fid,'B.14 COMMENT: IF RANDOMIZATION OF WAVE HEIGHT IS NOT DESIRED, CONTINUE TO B.16\r\n');
fprintf(fid,'B.15 SEED VALUE FOR RANDOMIZER AND PERCENT OF VARIABILITY: ISEED, RPERC\r\n');
fprintf(fid,' %d %.2f\r\n',OPT.ISEED,OPT.RPERC);
fprintf(fid,'B.16 TOTAL WATER ELEVATION INPUT (CONSTANT=0, VARIABLE=1): IELEV\r\n');
fprintf(fid,' %d\r\n',OPT.IELEV);
fprintf(fid,'B.17 COMMENT: IF WATER ELEVATION IS VARIABLE CONTINUE TO B.20\r\n');
fprintf(fid,'B.18 CONSTANT TOTAL WATER ELEVATION: TELEV\r\n');
fprintf(fid,' %.3f\r\n',OPT.TELEV);
fprintf(fid,'B.19 COMMENT: IF WATER ELEVATION IS CONSTANT, CONTINUE TO B.21\r\n');
fprintf(fid,'B.20 TIME STEP OF VARIABLE TOTAL WATER ELEVATION INPUT IN MINUTES: DTELV\r\n');
fprintf(fid,' %.2f\r\n',OPT.DTELV);
fprintf(fid,'B.21 WIND SPEED AND ANGLE INPUT (CONSTANT=0, VARIABLE=1): IWIND\r\n');
fprintf(fid,' %d\r\n',OPT.IWIND);
fprintf(fid,'B.22 COMMENT: IF WIND SPEED AND ANGLE ARE VARIABLE, CONTINUE TO B.25\r\n');
fprintf(fid,'B.23 CONSTANT WIND SPEED AND ANGLE: W,ZWIND\r\n');
fprintf(fid,' %.3f %.3f\r\n',OPT.W,OPT.ZWIND);
fprintf(fid,'B.24 COMMENT: IF WIND SPEED AND ANGLE ARE CONSTANT, CONTINUE TO C.\r\n');
fprintf(fid,'B.25 TIME STEP OF VARIABLE WIND SPEED AND ANGLE INPUT IN MINUTES: DTWND\r\n');
fprintf(fid,' %.2f\r\n',OPT.DTWND);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% C - BEACH
fprintf(fid,'C------------------------------- BEACH ----------------------------------C\r\n');
fprintf(fid,'C.1 TYPE OF INPUT PROFILE (ARBITRARY=1, SCHEMATIZED=2): TPIN\r\n');
fprintf(fid,' %d\r\n',OPT.TPIN);
fprintf(fid,'C.2 COMMENT: IF PROFILE TYPE IS ARBITRARY CONTINUE TO C.4\r\n');
fprintf(fid,'C.3 LOCATION AND ELEVATION OF LANDWARD BOUNDARY, LANDWARD BASE OF DUNE,\r\n');
fprintf(fid,' LANDWARD CREST OF DUNE, SEAWARD CREST OF DUNE, START OF BERM,\r\n');
fprintf(fid,' END OF BERM, AND FORESHORE: XLAND,DLAND,XLBDUNE,DLBDUNE,XLCDUNE,DLCDUNE,\r\n');
fprintf(fid,' XSCDUNE,DSCDUNE,XBERMS,DBERMS,XBERME,DBERME,XFORS,DFORS\r\n');
fprintf(fid,' %d %d %d %d %d %d\r\n',OPT.XLAND,OPT.DLAND,OPT.XLBDUNE,OPT.DLBDUNE,OPT.XLCDUNE,OPT.DLCDUNE);
fprintf(fid,' %d %d %d %d %d %d %d %d\r\n',OPT.XSCDUNE,OPT.DSCDUNE,OPT.XBERMS,OPT.DBERMS,OPT.XBERME,OPT.DBERME,OPT.XFORS,OPT.DFORS);
fprintf(fid,'C.4 DEPTH CORRESPONDING TO LANDWARD END OF SURF ZONE: DFS\r\n');
fprintf(fid,' %.3f\r\n',OPT.DFS);
fprintf(fid,'C.5 EFFECTIVE GRAIN SIZE DIAMETER IN MILLIMETERS: D50\r\n');
fprintf(fid,' %.3f\r\n',OPT.D50);
fprintf(fid,'C.6 MAXIMUM PROFILE SLOPE PRIOR TO AVALANCHING IN DEGREES: BMAX\r\n');
fprintf(fid,' %.3f\r\n',OPT.BMAX);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% D - BEACH FILL
fprintf(fid,'D------------------------------ BEACH FILL --------------------------------D\r\n');
fprintf(fid,'D.1 IS A BEACH FILL PRESENT? (NO=0, YES=1): IBCHFILL\r\n');
fprintf(fid,' %d\r\n',OPT.IBCHFILL);
fprintf(fid,'D.2 COMMENT: IF NO BEACH FILL, CONTINUE TO E.\r\n');
fprintf(fid,'D.3 POSITION OF START AND END OF BEACH FILL RELATIVE\r\n');
fprintf(fid,' TO INITIAL PROFILE: XBFS, XBFE\r\n');
fprintf(fid,' %.2f %.2f\r\n',OPT.XBFS,OPT.XBFE);
fprintf(fid,'D.4 NUMBER OF REPRESENTATIVE POINTS BETWEEN START\r\n');
fprintf(fid,'AND END OF BEACH FILL: NFILL\r\n');
fprintf(fid,' %d\r\n',OPT.NFILL);
fprintf(fid,'D.5 LOCATION AND ELEVATION OF REPRESENTATIVE POINTS RELATIVE TO THE\r\n');
fprintf(fid,'INITIAL PROFILE: (XF(I), EFILL(I), I=1,NFILL)\r\n');
for ii = 1:length(OPT.NFILL)
if isempty(OPT.XF)
fprintf(fid,'\r\n');
else
fprintf(fid,' %.2f %.2f',OPT.XF(ii),OPT.EFILL(ii));
end
end
fprintf(fid,'\r\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% E - SEAWALL/REVETMENT
fprintf(fid,'E--------------------------- SEAWALL/REVETMENT ----------------------------E\r\n');
fprintf(fid,'E.1 IS A SEAWALL PRESENT? (NO=0, YES=1): ISWALL\r\n');
fprintf(fid,' %d\r\n',OPT.ISWALL);
fprintf(fid,'E.2 COMMENT: IF NO SEAWALL, CONTINUE TO F.\r\n');
fprintf(fid,'E.3 LOCATION OF SEAWALL RELATIVE TO INITIAL PROFILE: XSWALL\r\n');
fprintf(fid,' %.2F\r\n',OPT.XSWALL);
fprintf(fid,'E.4 IS SEAWALL ALLOWED TO FAIL? (NO=0, YES =1): ISWFAIL\r\n');
fprintf(fid,' %d\r\n',OPT.ISWFAIL);
fprintf(fid,'E.5 COMMENT: IF NO SEAWALL FAILURE, CONTINUE TO F.\r\n');
fprintf(fid,'E.6 PROFILE ELEVATION AT SEAWALL WHICH CAUSES FAILURE, TOTAL WATER ELEVATION\r\n');
fprintf(fid,' AT SEAWALL WHICH CAUSES FAILURE, AND WAVE HEIGHT AT SEAWALL WHICH CAUSES\r\n');
fprintf(fid,' FAILURE: PEFAIL, WEFAIL,HFAIL\r\n');
fprintf(fid,' %.2f %.2f %.2f\r\n',OPT.PEFAIL,OPT.WEFAIL,OPT.HFAIL);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% F - HARD BOTTOM
fprintf(fid,'F------------------------------ HARD BOTTOM -----------------------------F\r\n');
fprintf(fid,'F.1 IS A HARD BOTTOM PRESENT? (NO=0, YES=1): IHBOT\r\n');
fprintf(fid,' %d\r\n',OPT.IHBOT);
fprintf(fid,'F.2 COEFFICIENT FOR TRANSPORT GROWTH DOWNDRIFT EXPOSED HARD BOTTOM: SCF\r\n');
fprintf(fid,' %.2f\r\n',OPT.SCF);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(fid,'----------------------------------- END ------------------------------------\r\n');
fclose(fid);
end

18
writetofileelv.m
Unescape Escape View File

@ -0,0 +1,18 @@
function writetofileelv(elvfile, outpath, pfstring)
fid = fopen([outpath '\' pfstring '.ELV'],'w');
fprintf(fid,'C> *******************************************************\r\n');
fprintf(fid,'C> * SBEACH water elevation input file: %s.ELV *\r\n',pfstring);
fprintf(fid,'C> *******************************************************\r\n');
fprintf(fid,'C>\r\n');
fprintf(fid,'C>\r\n');
fprintf(fid,'C> TOTAL WATER ELEVATION (m or ft)\r\n');
fprintf(fid,'E>----------------------------------------------------------------\r\n');
for ii = 1:length(elvfile)
if ii == length(elvfile)
fprintf(fid,'%.3f\r\n', elvfile(ii));
else
fprintf(fid,'%.3f\r\n', elvfile(ii));
end
end
fclose(fid);
end

20
writetofilehdb.m
Unescape Escape View File

@ -0,0 +1,20 @@
function writetofilehdb(hdbfile, outpath, namestr)
fid = fopen([outpath '\' namestr '.HDB'],'w');
fprintf(fid,'C> *******************************************************\r\n');
fprintf(fid,'C> * SBEACH hard bottom input file: %s.HDB *\r\n',namestr);
fprintf(fid,'C> *******************************************************\r\n');
fprintf(fid,'C>\r\n');
fprintf(fid,'C>\r\n');
fprintf(fid,'C> NUMBER OF POINTS\r\n');
fprintf(fid,'C> POSITION (m or ft), ELEVATION (m or ft)\r\n');
fprintf(fid,'E>----------------------------------------------------------------\r\n');
fprintf(fid,'%d\r\n',length(hdbfile));
for ii = 1:length(hdbfile)
if ii == length(hdbfile)
fprintf(fid,'%.3f %.3f\r\n', hdbfile(ii,1),hdbfile(ii,2));
else
fprintf(fid,'%.3f %.3f\r\n', hdbfile(ii,1),hdbfile(ii,2));
end
end
fclose(fid);
end

20
writetofilepri.m
Unescape Escape View File

@ -0,0 +1,20 @@
function writetofilepri(prifile, outpath, namestr)
fid = fopen([outpath '\' namestr '.PRI'],'w');
fprintf(fid,'C> *******************************************************\r\n');
fprintf(fid,'C> * SBEACH initial profile input file: %s.PRI *\r\n',namestr);
fprintf(fid,'C> *******************************************************\r\n');
fprintf(fid,'C>\r\n');
fprintf(fid,'C>\r\n');
fprintf(fid,'C> NUMBER OF POINTS\r\n');
fprintf(fid,'C> POSITION (m or ft), ELEVATION (m or ft)\r\n');
fprintf(fid,'E>----------------------------------------------------------------\r\n');
fprintf(fid,'%d\r\n',length(prifile));
for ii = 1:length(prifile)
if ii == length(prifile)
fprintf(fid,'%.3f %.3f\r\n', prifile(ii,1),prifile(ii,2));
else
fprintf(fid,'%.3f %.3f\r\n', prifile(ii,1),prifile(ii,2));
end
end
fclose(fid);
end

20
writetofileprm.m
Unescape Escape View File

@ -0,0 +1,20 @@
function writetofileprm(prmfile, outpath, pfstring, rundescription)
fid = fopen([outpath pfstring '\' pfstring '_' rundescription '.PRM'],'w');
fprintf(fid,'C> *******************************************************\n');
fprintf(fid,'C> * SBEACH initial profile input file: %s.PRM *\n',[pfstring '_' rundescription]);
fprintf(fid,'C> *******************************************************\n');
fprintf(fid,'C>\n');
fprintf(fid,'C>\n');
fprintf(fid,'C> NUMBER OF POINTS\n');
fprintf(fid,'C> POSITION (m or ft), ELEVATION (m or ft)\n');
fprintf(fid,'E>----------------------------------------------------------------\n');
fprintf(fid,'%d\n',length(prmfile));
for ii = 1:length(prmfile)
if ii == length(prmfile)
fprintf(fid,'%4.3f %3.3f', prmfile(ii,1),prmfile(ii,2));
else
fprintf(fid,'%4.3f %3.3f\n', prmfile(ii,1),prmfile(ii,2));
end
end
fclose(fid);
end

18
writetofilewav.m
Unescape Escape View File

@ -0,0 +1,18 @@
function writetofilewav(wavfile, outpath, pfstring)
fid = fopen([outpath '\' pfstring '.WAV'],'w');
fprintf(fid,'C> *******************************************************\r\n');
fprintf(fid,'C> * SBEACH wave input file: %s.WAV *\r\n',pfstring);
fprintf(fid,'C> *******************************************************\r\n');
fprintf(fid,'C>\r\n');
fprintf(fid,'C>\r\n');
fprintf(fid,'C> HEIGHT (m or ft), PERIOD (sec)\r\n');
fprintf(fid,'E>----------------------------------------------------------------\r\n');
for ii = 1:length(wavfile)
if ii == length(wavfile)
fprintf(fid,'%.3f %.3f\r\n', wavfile(ii,1),wavfile(ii,2));
else
fprintf(fid,'%.3f %.3f\r\n', wavfile(ii,1),wavfile(ii,2));
end
end
fclose(fid);
end
Write Preview
Loading…
Cancel
Save