%% Analyzing model output

load('name'); load(basename); 

zGr  =Z(1,1,:);

[xGr,yGr,zGr,xm,ym,zm,dx,dy,dz,Nx,Ny,Nz]=modelsize3(xGr,yGr,zGr);

H=readDat([basename,'.hds']);
C=readMT3D('MT3D001.UCN');
B=readBud([basename,'.bgt'],'FLOWRIGHTFACE');

Psi=XS(B(end).term{1}); Psi=Psi(:,1:end-1);
Psi=[flipud(cumsum(Psi(end:-1:1,:),1)); zeros(size(Psi(1,:)))];

figure
contourf(xm,squeeze(zm),XS(C(end).values)); hold on
contour (xm,squeeze(zm),XS(H(end).values),'r');
contour (xGr(2:end-1),squeeze(zGr),Psi,'y');
title('Henry problem Seawat, [conc, heads and stream function]');
xlabel('x [m]');
ylabel('z [m]');
