function ph = duwa_check_volume_conservation(varargin)

OPT = struct(...
    'irun', 2,...
    'maindir', fullfile(fileparts(which(mfilename)), '..', 'simulations', '00000000', 'xbeach'));

OPT = setproperty(OPT, varargin{:});

ncfiles = dir(fullfile(OPT.maindir, sprintf('run%02i', OPT.irun), '*.nc'));
ncfile = fullfile(OPT.maindir, sprintf('run%02i', OPT.irun), ncfiles.name);

angledirs = dir(fullfile(fileparts(ncfile), 'angle*'));
anglestrings = {angledirs.name};

x = nc_varget(ncfile, 'globalx');
y = nc_varget(ncfile, 'globaly');
zb = nc_varget(ncfile, 'zb');

sedero = squeeze(zb(:,2,:,:) - zb(:,1,:,:));

xm = [x(:,1) mean(cat(3, x(:,1:end-1), x(:,2:end)),3) x(:,end)];
xd = diff(xm,1,2);

ym = [y(1,:); mean(cat(3, y(1:end-1,:), y(2:end,:)),3); y(end,:)];
yd = diff(ym,1,1);

[xx yy] = deal(NaN(size(sedero)));
for i = 1:size(yy,1);
    xx(i,:,:) = xd;
    yy(i,:,:) = yd;
end

Vnet = xx .* yy .* sedero;
Vero = -Vnet;
Vero(sedero>0) = 0;
Vsed = Vnet;
Vsed(sedero<0) = 0;

fprintf('Nett volumes:\n')
fprintf('%g\n',  sum(sum(Vnet,3),2));
fprintf('Erosion volumes:\n')
fprintf('%g\n',  sum(sum(Vero,3),2));
fprintf('Accretion volumes:\n')
fprintf('%g\n',  sum(sum(Vsed,3),2));

ph = plot(y(:,1:size(zb,1)), (sum(Vero,3) ./ squeeze(yy(:,:,1)))');
legend(anglestrings{:})
ylabel('Erosion volume [m^3/m]')
xlabel('Alongshore coordinate [m]')
