%% Analyzing output of the model
% TO 100415
 
%% Bootstrapping to get the model data back in a generic way
load('name.mat') % loads name.mat, which only contains the variable "basename"
load(basename);  % having retrieved the baasename value, we can load

%% Size of the model
[NROW,NCOL,NLAY]=size(IBOUND); % get the size of the model

%% load the unformatted head file
[mfparnams,mfparvals]=getExcelData(basename,'MFLOW','vertical');
HNOFLO=mfparvals(strmatchi('HNOFLO',mfparnams),1);
HDRY  =mfparvals(strmatchi('HDRY',mfparnams),1);

H=readDat([basename,'','.hds']);  % use readDAT to read the heads  output file
H=maskHC(H,IBOUND);               % throws out NaN's
H=maskHC(H,HNOFLO);
H=maskHC(H,HDRY);

%% plot heads as surfaces
[xGr,yGr,xm,ym,Dx,Dy,Ny,Nx]=modelsize(xGr,yGr); % xGr and yGr were saved in mf_adapt
zm=0.5*(Z(:,:,1)+Z(:,:,2));

figure; hold on
plot(xm,Z(:,:,end),'color','k','linewidth',3);
plot(xm,STRTHD(:,:,end),'m');

for it=1:length(H)
    for iL=1:size(H(1).values,3)
        plot(xm,H(it).values(:,:,iL),'g'); hold on
    end
end
plot(xm,H(end).values(:,:,end),'r');

%% Read unformatted budget file and mask noflow cells if they exist

B=readBud([basename,'','.bgt']);  B=maskHC(B,IBOUND);

totCHD=sum( B(1).term{ strmatchi('CONSTANTHEAD',B(1).label) }(:) );

fprintf('total constant head discharge = %12.2f\n',totCHD);

%% In and outflow through constant head cells
CH=B(1).term{strmatch('CONSTANTHEAD',B(1).label)};
QCHin =sum(CH(CH(:)>0));
QCHout=sum(CH(CH(:)<0));

fprintf(['Total inflow  across constant head cells    =  %12.2f m^3/d\n',...
         'Total outflow across constant head cells    =  %12.2f m^3/d\n',...
         '                                               -----------------\n',...
         'Total net inflow across constant head cells =  %12.2f m^3/d\n'],...
         QCHin,QCHout,QCHin+QCHout);

%% Zonebudget
% see ex1 for usage of Zonebudget
