ncfile = 'p:\1002266-xbeach\coastal_curvature\simulations\00000000\xbeach\run20\radius01900\xboutput.nc';

info = nc_info(ncfile);
tdim = info.Dimension(strcmp({info.Dimension.Name}, 'globaltime')).Length;

x = nc_varget(ncfile, 'globalx');
y = nc_varget(ncfile, 'globaly');
t = nc_varget(ncfile, 'globaltime');
g = xb_stagger(x, y);

fprintf('tstart: %g\ntstop: %g\nduration: %g\n', t(1), t(tdim), diff(t([1 tdim])))

%% check cell area
dsz = g.dsz;
dsz(1,:) = 0;
sz = cumsum(dsz,1);
zb0 = nc_varget(ncfile, 'zb', [0 0 0], [1 -1 -1]);
%%
for i = 1:size(sz,2)
    sref(i) = findCrossings(sz(:,i), zb0(i,:)', [min(sz(:)) max(sz(:))], [5 5]);
end
rz = repmat(1900 + sref, 160, 1) - sz;
dalfaz = diff(g.alfaz,1,2);
dalfaz(:,end+1) = dalfaz(:,end);

dsdnz = (rz(1:end-1,:).^2 - rz(2:end,:).^2) .* dalfaz(1:end-1,:) / 2;

%% check integrated sedero
sedero = nc_varget(ncfile, 'sedero', [tdim-1 0 0], [1 -1 -1]);
% sedero =  nc_varget(ncfile, 'dzS', [tdim-1 0 0], [1 -1 -1]);
S_sedero = g.dsdnz .* sedero';
fprintf('Total erosion: %g\n', -sum(S_sedero(S_sedero<0)))
fprintf('Total sedimentation: %g\n', sum(S_sedero(S_sedero>0)))
fprintf('Integrated sedero: %g\n', sum(S_sedero(:)));

%% plot integrated sedero cross-shore
figure;
plot(cumsum(sum(S_sedero,2)), '-o')

%% check integrated dzav
dzav = nc_varget(ncfile, 'dzav', [tdim-1 0 0], [1 -1 -1]);
S_dzav = g.dsdnz .* dzav';
fprintf('Integrated dzav: %g\n', sum(S_dzav(:)));

%% plot integrated dzav cross-shore
figure;
plot(cumsum(sum(S_dzav,2)), '-gs')

%% check water volume
zs = nc_varget(ncfile, 'zs', [0 0 0], [2 -1 -1], [tdim-1 1 1]);
zb = nc_varget(ncfile, 'zb', [0 0 0], [2 -1 -1], [tdim-1 1 1]);
hh = zs - zb;
S_hh0 = squeeze(hh(1,:,:))' .* g.dsdnz;
fprintf('Integrated hh @ t=0: %g\n', sum(S_hh0(:)));
S_hh1 = squeeze(hh(end,:,:))' .* g.dsdnz;
fprintf('Integrated hh @ t=end: %g\n', sum(S_hh1(:)));

fprintf('Difference integrated hh: %g\n', sum(S_hh1(:)) - sum(S_hh0(:)))

%% check sediment transport at lateral boundaries
Svsg = nc_varget(ncfile, 'Svsg_mean');
S_Svsg = Svsg' .* g.dsu;
fprintf('Integrated transport over lateral boundary (n=1): %g\n', sum(S_Svsg(:,1)))
fprintf('Integrated transport over lateral boundary (n=end): %g\n', sum(S_Svsg(:,end)))

%% check sediment transport over offshore boundary
Susg = nc_varget(ncfile, 'Susg_mean');
S_Susg = Susg' .* g.dnv;
fprintf('Integrated transport over offshore boundary (s=1): %g\n', sum(S_Susg(1,:)))

%%
morfac = 10;
por = .4;
dsdnzi = 1./g.dsdnz;
sedero = nc_varget(ncfile, 'sedero');
%%
Svsg = nc_varget(ncfile, 'Svsg');
Susg = nc_varget(ncfile, 'Susg');
Svbg = nc_varget(ncfile, 'Svbg');
Subg = nc_varget(ncfile, 'Subg');
%%
dt = .2;
it = 1;
i = 2:size(x,2);
j = 2:size(x,1);
imin1 = i-1;
jmin1 = j-1;
for it = 1:size(Susg,1)-1
dzg(it+1,:,:) = morfac * dt / por * ...
(...
(...
squeeze(Susg(it,j,i))' .* g.dnu(i,j) - squeeze(Susg(it,jmin1,imin1))' .* g.dnu(imin1,jmin1) +...
squeeze(Svsg(it,j,i))' .* g.dnv(i,j) - squeeze(Svsg(it,jmin1,imin1))' .* g.dnv(imin1,jmin1) +...
squeeze(Subg(it,j,i))' .* g.dnu(i,j) - squeeze(Subg(it,jmin1,imin1))' .* g.dnu(imin1,jmin1) +...
squeeze(Svbg(it,j,i))' .* g.dnv(i,j) - squeeze(Svbg(it,jmin1,imin1))' .* g.dnv(imin1,jmin1)  ...
) .* dsdnzi(i,j) ...
)';
end
dzg(1,:,:) = 0;
sedero2 = cumsum(dzg,1);
%%
temp = sedero(it,j,i) - dzg(it,:,:);
figure; pcolor(x(j,i), y(j,i), squeeze(temp))
colorbar