%% Analyzing output of the model
% TO 091011 091129
 
%% Bootstrapping to get the model data back in a generic way
%  this bootstrapping allows clearing the workspace before to make sure
%  there are no other variables than those pertaining to the model

load('name.mat') % loads name.mat, which only contains the variable "basename"
                 % having as its value the basename of the model.
                 
load(basename);  % having retrieved the baasename value, we can load
                 % the <<basename>>.mat file which holds the
                 % model variables generated in mf_adapt

%% Size of the model
% IBOUND is always in the <<basename>>.mat file so we use it to regain the
% size of the model

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

%% load the unformatted head file
% To present output we will read the unformatted output files produced by
% running the model. In this case we have the unformatted heads file and
% the unformatted budget file

H=readDat([basename,'','.hds']);  % use readDAT to read the heads  output file
            % readDAT will read any file with the same format, like the
            % drawdown file. The result is a struct whose length is the
            % number of timesteps stored in the file. The struct has
            % various fields that can be used for processing and
            % visualisation.
            % In this case H is only one record long. Showing H(1) yields
            % this list of parameters held by the struct H
 %H(1)=
 %            values: [15x15x3 double]
 %             tstp: 1
 %           period: 1
 %           pertim: 86400
 %            totim: 86400
 %            label: 'HEAD'
 %   selectedLayers: [1 2 3]
 %             NROW: 15
 %             NCOL: 15
 %             lays: [3x1 double]
 %             rows: [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15]
 %             cols: [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15]
 %
 % Because readDAT has many options to select exactly what you desire from
 % the file, for instance just a couple of time steps and a couuple of
 % layers rows and columns, it also shows the row, column and layer numbers
 % of the model that were selected. The actual data are in H(1).values.
 % This is a 3D array of the heads, only the selected heads, which may be
 % just a couple of columns rows and layers instead of everything.
 % type help readDat to get more information on how to make precise
 % selections (which are useful if datasets get very large or in cases
 % where you know exactly which items and cells you want).

% Sometimes H contains nodata values, for instance to indicated inactive
% cells and when cell fell dry during the computation. Generally MODFLOW
% uses some very high value stored in HONOFLOW in the MFLOW worksheet. To
% present contours of the heads, we have to remove these values. The most
% straightforward method in MATLAB is replacing them with a NaN, as NaNs do
% not interfere with plotting and contouring; NaNs just leave gaps, which
% is eactly what we desire in such cases.

H.values(H.values>1000)=NaN;      % remove possible inactive cells

%H=maskHC(H,IBOUND);  % Dos the same but in a more general fashion.
% It uses the IBOUND array itself to mask the inactive cells with NaNs.
% However maskHC would not mask dry celle that have come about during the
% computations.

% To get a list of all time steps and times in struct H do this
%tstp= [H.tstp];
%totim=[H.totim];
% of course this is only useful if H has several of many time steps instead
% of just one as in this example

%% plot heads as surfaces

% Let's plot the computed heads, and let's use axis in units of 1000 ft or
% m to make easier to read distance labels
% Notice that xGr and yGr are the gridline coordinates but xm and ym are
% the coordinates of the cell centers. We need these because we only have
% values at the cell centers.

[xGr,yGr,xm,ym,DELX,DELY,NX,NY]=modelsize(xGr,yGr); % xGr and yGr were saved in mf_adapt

figure
for iLay=1:NLAY    % plot 3D surface of every layer of model
    subplot(NLAY,1,iLay)
    h=surf(xm/1000,ym/1000,H.values(:,:,iLay)); hold on;
    title(sprintf('Example 1 of mf2k Open file Report 00-92 p89 [Layer %d]',iLay));
    xlabel('x [1000ft]');  ylabel('y [1000ft]');
end

% This shows that the difference between te layers is small
%% Contouring heads

% the most general way of showing results of a model is by contouring the
% heads, drawdowns etc.
% we can just contour any layer
figure
iLay=3;
[c,hdl]=contour(xm,ym,H.values(:,:,iLay));  % default contouring
clabel(c,hdl);                          % default labels
title(sprintf('Head contours, Open file Report 00-92 p89 [Layer %d]',iLay));
xlabel('x [1000ft]');  ylabel('y [1000ft]');

%% Specific contouring
% sometimes we rather specify the head contour steps and the entire contour
% range ourselfs
dPhi=1;                              % contour step size
phiMin=floor(min(H(end).values(:))); % round down min head (first contour)
phiMax=ceil( max(H(end).values(:))); % round up max head (last contour)
phirange=phiMin:dPhi:phiMax;         % contours to be drawn

figure                               % new picture
for iLay=1:NLAY
    subplot(1,NLAY,iLay);
    [c,hdl]=contourf(xm/1000,ym/1000,H.values(:,:,iLay));
    clabel(c,hdl);
    axis('equal'); % make both axes same scale
    axis('tight'); % no white area around axes
    caxis(phirange([1 end]));  % fix color scale to these values
    colorbar;      % place colorbar but now with known scale
    xlabel('x [1000ft]'); if iLay==1, ylabel('y [1000ft]'); end
    title(sprintf('period=%d, tstp=%d, layer=%d',...
                   H(end).period,H(end).tstp,iLay));      
end

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

% The next analysis of model output may be a water balance of the amount of
% water discharged by drains etc.
% For this we need to read the unformatted budget file that modflow generated.

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

% Like it is the case with readDat, B is a struct array with as many
% members as there are time steps in the budget file. Type B to see what
% fields this struct has. Type B(5) to see the actual values of the fifths
% stored time step (works only if length(B)>5 of course).
% As is the case with readDat you can specify what exactly is to be read
% from the budget file, which flow terms (see B.labels for the flow terms
% that are in the file), which times, time steps, rows, columns and layers.
% So if you want the WELLS you can tell readBud to extract the
% corresponding cell by cell flow terms. Altenratively readBud can read out
% everything and you can get the wells from the struct B directly
% afterwards.
% The actual flow term data of the first time step are now in
% B(1).term
% however, term is a cell array with as many cells as there are unique
% labels in the file. So if WELLS is the third label in B.labels then
%B(1).term{3}
%is teh actual 3D array of the sought cell by cell WELL terms
% To get specific rows, columen and layers from it goes like this
%B{1}.term{3}(rows,columns,layers)
% Hence to get the total drain discharge from the model, we may do the
% following
totDRN=sum( B(1).term{ strmatch('DRAINS',      B(1).label) }(:) );
totRCH=sum( B(1).term{ strmatch('RECHARGE',    B(1).label) }(:) );
totCHD=sum( B(1).term{ strmatch('CONSTANTHEAD',B(1).label) }(:) );
totWEL=sum( B(1).term{ strmatch('WELLS'       ,B(1).label) }(:) );

fprintf(['total recharge is             = %12.2f ft^3/d\n',...
         'total nwt well discharge      = %12.2f ft^3/d\n',...    
         'total constant head discharge = %12.2f ft^3/d\n',...
         'total drain discharge         = %12.2f ft^3/d\n',...
         '                                ------------------\n',...
         'overall total (balance)       = %12.f  ft^3/d\n'],...
            totRCH,totWEL,totCHD,totDRN, totRCH+totWEL+totCHD+totDRN);

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

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

%% Zonebudget
iz=reshape(1:3,[1,1,3]);  % one zone per layer
Zonearray=iz(ones(NROW,1),ones(1,NCOL),:);
zonebudget(B,Zonearray,1);
zonebudget(B,Zonearray,2);
zonebudget(B,Zonearray,3);
zonebudget(B,Zonearray,[1 3]);
zonebudget(B,Zonearray,[1 2 3]);



