%% Analyzing model output

[xGr,yGr,zGr,xm,ym,zm,Dx,Dy,Dz,Nx,Ny,Nz]=modelsize3(xGr,yGr,zGr);


ZM=0.5*(Z(:,:,1:end-1)+Z(:,:,2:end));

%% Reading unformattted files
H=readDat([basename,'.hds']);
T=readMT3D('MT3D002.UCN');
C=readMT3D('MT3D001.UCN');
B=readBud([basename,'.cbb'],'FLOWRIGHTFACE');

%% Stream function ??
NT=length(B);
Psi=NaN(Nz+1,Nx-1,NT);
for it=1:NT
    FRF=XS(B(it).term{1}(:,1:end-1,:));
    Psi(:,:,it)=flipud(cumsum(flipud([FRF; zeros(1,Nx-1)])));
end
%% Contouring concentrations
Crange=0+[0.05 0.5 0.95]*(35-0);
Trange=5+[0.05 0.5 0.95]*(25-5);
Hrange=10;
Prange=10;
zm=squeeze(zm);
for it=1:NT
    figure
    contour(xm,zm,XS(C(it).values),Crange,'b'); hold on
    contour(xm,zm,XS(T(it).values),Trange,'k'); hold on
    contour(xm,zm,XS(H(it).values),Hrange,'r');  hold on
    contour(xGr(2:end-1),squeeze(zGr),Psi(:,:,it),Prange,'k'); % hold on
    title(sprintf('%s problem, t=%.0f',basename,C(it).time));
    xlabel('x [m]');
    ylabel('z [m]');
end
