function postpr_xbeach (dirxbc,dirfig);

cd (dirxbc);

%% get model dimensions

fid=fopen('dims.dat','r');
nt=fread(fid,[1],'double'); % nr of timesteps
nx=fread(fid,[1],'double'); % nr of x-points
ny=fread(fid,[1],'double'); % nr of y-points
fclose(fid);

%% get x and y values of grid

fixy=fopen('xy.dat','r');
x=fread(fixy,[nx+1,ny+1],'double'); % x values
y=fread(fixy,[nx+1,ny+1],'double'); % y values
fclose(fixy);

%% get datatype to plot

fizb=fopen('zb.dat','r');
for i=1:nt;
    zb=fread(fizb,[nx+1,ny+1],'double');
    if i==1;  zb_0=zb;   end;
    if i==nt; zb_end=zb; end;
end
fclose(fizb);

figure(1);clf;pcolor(x,y,zb_0);
shading interp; axis equal; axis image;
title('Initial bathymetry Egmond');
caxis([-10 10]);
cb=colorbar;set(get(cb,'ylabel'),'string','Vertical elevation in m to NAP')
pngnam = [dirfig '\InitialBathymetryEgmond'];
print (gcf,'-dpng',pngnam);
close;


figure(2);clf;pcolor(x,y,zb_end);
shading interp; axis equal; axis image;
title('Final bathymetry Egmond');
caxis([-10 10]);
cb=colorbar;set(get(cb,'ylabel'),'string','Vertical elevation in m to NAP')
pngnam = [dirfig '\FinalBathymetryEgmond'];
print (gcf,'-dpng',pngnam);
close;

figure(3);clf;pcolor(x,y,zb_end-zb_0);
shading interp; axis equal; axis image;
title('\deltaz_b [m]');
caxis([-0.1 0.1]);
cb=colorbar;set(get(cb,'ylabel'),'string','Morphological changes in m')
pngnam = [dirfig '\SedimentationErosionEgmond'];
print (gcf,'-dpng',pngnam);
close;