mf_settings;

load('name.mat') % loads name.mat, which only contains the variable "basename"                                
load(basename);  % having retrieved the baasename value, we can load
                 
load('adapt_settings.mat') % Load settings added by Yannick

H=readDat([basename,'','.hds']);  % use readDAT to read the heads  output file
B=readBud([basename,'','.bgt']);  B=maskHC(B,IBOUND);


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

[xGr,yGr,xm,ym,DELX,DELY,NX,NY]=modelsize(xGr,yGr); % xGr and yGr were saved in mf_adapt
% Read unformatted budget file and mask noflow cells if they exist

if strcmpi( FlowDirection, 'horizontal'),    
    A = (Z(1,1,1)-Z(1,1,2)) * (yGr(1)-yGr(end));
    dl = abs( xGr(end)-xGr(1) );
elseif strcmpi( FlowDirection, 'vertical')   
    A = (Z(1,1,1)-Z(1,1,2)) * (xGr(end)-xGr(1));
    dl = abs( yGr(end)-yGr(1) );
end
% hcellsize = dl / size(MAP, 2);

%% In and outflow through constant head cells
dh = STRTHD(1, 1) - STRTHD(end, end);

CH=B.term{strmatch('CONSTANTHEAD',B.label)};
QCHin =sum(CH(CH(:)>0));
QCHout=sum(CH(CH(:)<0));

bulkconductivity = QCHin / ( A * dh/dl) ;
bulkpermeability = bulkconductivity / (2.75 .* 0.001 .* (1/3.28)) ;

mass_balance_error = abs( 100*( abs(QCHin)-abs(QCHout))/abs(QCHin) );

fprintf( '\nBulk K: %e mD    MBE: %e %% \n', [ bulkpermeability, mass_balance_error]);