%% XBeach advanced course: non-hydrostatic modelling
% Modeling of nature-based flood protection solutions with XBeach-NH
% v1.0  de Ridder     Nov-18
clear all
close all

%% 0. Input
destout     = 'p:\xbeach\course_material\DSD_2019\02_advanced_nonh\hands-on\vegetation\_run\';
length      = 75;   % width of the overtopping basin

%% 1. Read parameters
% Retrieve valeus from netcdf
fname            = [destout '\xboutput.nc'];
point_zs         = nc_varget(fname, 'point_zs'); 
var_zs           = nc_varget(fname, 'zs_var'); 
zb               = nc_varget(fname, 'zb'); zb = squeeze(zb(1,1,:));
globalx          = nc_varget(fname, 'globalx'); 
pointtime        = nc_varget(fname, 'pointtime');

% Compute wave height
ind  = find(zb<0.5);        % show results where the waterdepth is deeper than 0.5
Hrms = sqrt(8 * var_zs);

%% 2. Plot wave height
figure;
plot(globalx(ind),Hrms(ind))
hold on
plot(globalx,zb)
legend('H_{rms}','zb')
xlabel('Distance [m]');ylabel('Elevation [m]')

%% 3. Compute overtopping
% Water level variation over time
figure;
plot(pointtime,point_zs(:,2)); title('water level basin')
xlabel('time [s]');ylabel('Elevation [m]')

% Compute overtopping rate
p   = polyfit(pointtime,point_zs(:,2),1);

% See result
figure; hold on;
plot(pointtime,point_zs(:,2), '.'); title('water level basin')
plot(pointtime, polyval(p,pointtime));
xlabel('time [s]');ylabel('Elevation [m]')
legend('point','linear trend')

% Compute rate method 1
q   = p(1) * length * 1000;             % [m/s] * [m] * [m3 to l]
fprintf(['The mean discharge per m is: ' num2str(q) ' l/s/m \n'])

% compute rate - method 2
wl  = point_zs(end,2) - point_zs(1,2);
V   = wl * length * 1000; % volume in liter
q = V/(pointtime(end)-pointtime(1));
fprintf(['The mean discharge per m is: ' num2str(q) ' l/s/m \n'])