function H=readDat(fname,userPeriod,userTstp,userLayers,userRows,userCols)
% H=readDat([-]fname [,userPeriod [,tsteps [,userLayers [,userRows [,userCols]]]]]])
% ---- read binary MODFLOW output into a Data(nP,nT) array of structs
% use
%      H=readDat(fname,-1)
% to get info on the contents of the file
%
% make subselections using the optoinal arguments
% they may be omittted as shown or used empty '' or [], for instance
%
%   H=readDat(fname,'',[1,3,6],'',3)
%
% fetches all userPeriod, time setps 1,3 and 6 all layers, row 3 and all
% columns
%
%   H=readDat(fname,[1 3 5:8],'',5,3:10,12:21}
%
% fetches userPeriod 1,3 and 5..8, all time steps, layer 3, userRows 3..10 and
% columns 12..21
%
% prepend '-' to fname to get less verbose output
%
% H is a cell array whose length is the number of records in the file
% it has the following fields:
%    .values is a 3D array containing the requested concentration values
%    .label        label telling the contents
%    .period       stress period number
%    .tstp         time step number within stress period
%    .pertim       time since begin of stress period
%    .totim        time since begin of simulation
%    .NCOL         number of original columns in file fname
%    .NROW         number of original rows    in file fname
%    .NLAY         number of original layers  in file fname
%    .userCols         list of fetched columns
%    .userRows         list of fetched rows
%    .userLayers         list of fetched layers
%
% TO 090104 091214
% TO 100818 logic redone to clean up necessary to fix number of ouput recs
%
% Copyright 2009 Theo Olsthoorn, TU-Delft and Waternet, without any warranty
% under free software foundation GNU license version 3 or later

if fname(1)=='-', fname=fname(2:end); verbose=0; else verbose=1; end

fp=fopen(fname); if fp<1, error('READDAT: Can''t find or open file %s!',fname); end

fprintf('Reading MODFLOW binary output file <<%s>> verbose= %d\n',fname,verbose);

%% Advance a single cell-by-cell record & compute the number or records in file

bytes=contrRec(fp);

nbyt=ftell(fp); fseek(fp,0,'eof'); Nbyt=ftell(fp); nRec=Nbyt/nbyt; % number of records in file

if rem(Nbyt,nbyt)~=0,
    error('BINARY file %s not of standard type. Sorry, can''t read it, use a better compiler!',fname);
end

%% Get the information about the contents of the file from record headings

periodsInFile = zeros(nRec,1); % stress period number of each record
tstepsInFile  = zeros(nRec,1); % time step number of each record in file 
labelInFile   = cell( nRec,1); % label of each record in file
pertim        = zeros(nRec,1); % time since beginning of stress period
totim         = zeros(nRec,1); % time since beginning of simulation
layersInFile  = zeros(nRec,1); % layer number of each record in file
NROW          = zeros(nRec,1); % NROW specified in each record in file
NCOL          = zeros(nRec,1); % NCOL specified in each record in file

fprintf('Scanning headers\n');
for i=1:nRec  % get meta data for each record in file
    fseek(fp,(i-1)*nbyt,'bof');  % move to start of last layer in file
    
    [tstepsInFile(i),...
        periodsInFile(i),...
        pertim(i),...
        totim(i),...
        labelInFile{i},...
        layersInFile(i),...
        NROW(i),...
        NCOL(i)...
        ]=contrRec(fp,bytes);
    
    if rem(i, 100)==0, fprintf('.'); end
    if rem(i,2500)==0, fprintf('%d records read\n',i); end
end
fprintf('finished, %d records read\n',i);

%% Evaluate what we've got in terms of the userPeriod, time steps unique lables etcetera before we continue to select the actual data

NPER=max(periodsInFile);
NSTP=max(tstepsInFile);
NLAY=max(layersInFile);
NROW=max(NROW);
NCOL=max(NCOL);

fprintf('File contains the following:\n');
fprintf('Number of records in file  : %10d\n',nRec);
fprintf('Number of stress userPeriod: %10d\n',NPER);
fprintf('Number of time steps       : %10d\n',NSTP);
fprintf('Number of layers           : %10d\n',NLAY);
fprintf('Number of rows             : %10d\n',NROW);
fprintf('Number of columns          : %10d\n',NCOL);

if exist('userPeriod','var') && ~isempty(userPeriod) && userPeriod(1)<=0, return; end

%% select the speciied zones or all if nothing specified
%periodsInFile  the period numbers in the InFile
%tstepsInFile   the time step numbers in the InFile
%layersInFile   the layer numbefs in the InFile

%% the period numbers requested by the user or all of them
if  exist('userPeriod','var') && ~isempty(userPeriod)   % if the user specified which stress periods
    if ~isempty(userPeriod(userPeriod<1 | userPeriod>NPER))
        error('Requested userPeriod beyond range(1..%d) in dat file %s!',NPER,fname);
    end
    userPeriod=unique(userPeriod);  % this is set of user requrested stress userPeriod
else
    userPeriod=unique(periodsInFile);
end

%% the time step numbers requested by the user or all of them
if exist('userTstp'  ,'var') && ~isempty(userTstp)  % if the user specified which time steps he wants
    if ~isempty(userTstp(userTstp<1 | userTstp>NSTP))
        error('Requested userPeriod beyond range(1..%d) in buget file %s!',NSTP,fname);
    end
    userTstep=unique(userTstp);
else
    userTstep=unique(tstepsInFile);
end

%% the layer numbers requested by the user of all of them
if exist('userLayers'   ,'var') && ~isempty(userLayers)  % if the user specifies which layers he/she wants
    if ~isempty(userLayers(userLayers<1 | userLayers>NLAY))
        error('Requested layers are beyond range(1..%d) in dat file %s!',NLAY,fname);
    end
    userLayers=unique(userLayers);
else
    userLayers=unique(layersInFile);
end

%% Selectors have values for all records of the inptut file. There
%  contents are userPeriod index, the userTstep index, the userLayers index
%  the combined userPeriod and userTstep index plus
%  the output record index and the output layer index
UPCOL=1; UTCOL=2; ULCOL=3; UPTCOL=4; OLCOL=6;

Select=zeros(nRec,OLCOL);  
for i=1:length(userPeriod), Select(periodsInFile==userPeriod(i),UPCOL)=i; end
for i=1:length(userTstep),  Select(tstepsInFile ==userTstep(i) ,UTCOL)=i; end
for i=1:length(userLayers), Select(layersInFile ==userLayers(i),ULCOL)=i; end

% input recoreds to deal with
IPT=find(Select(:,UPCOL)>0 & Select(:,UTCOL)>0);

% output records
userPT = unique([periodsInFile(IPT) tstepsInFile(IPT)],'rows');

for i=1:size(userPT,1), Select(periodsInFile ==userPT(i,1) & tstepsInFile==userPT(i,2),UPTCOL)=i; end

%% rows and cols do not need an indicater array to pick out the correct
%  data, the lists are used directly

if exist('userRows'   ,'var') && ~isempty(userRows)
    if ~isempty(userRows(userRows<1 | userRows>NROW))
        error('Requested rows are beyond the range(1..%d) in dat file %s!',NROW,fname);
    end
else
    userRows=1:NROW;
end

if exist('userCols'   ,'var') && ~isempty(userCols)
    if ~isempty(userCols(userCols>NCOL | userCols<1))
        error('Requested columns are beyond the range(1..%d) in dat file %s!',NCOL,fname);
    end
else
    userCols=1:NCOL;
end

%% What are the unique periods and tstep combinations? Because the number of

nRecOut = size(userPT,1);

B.values=NaN(length(userRows),length(userCols),length(userLayers)); % dimension of single outrec

H=repmat(B,[nRecOut,1]); % size of output sttruct array is allocated here

%% Get the actual data values

for i=1:length(IPT)  % indices of required records in infile
    iRecIn=IPT(i);
    fseek(fp,nbyt*(iRecIn-1),'bof');
    iROut=Select(iRecIn,UPTCOL);   % output record Nr
    iL=Select(iRecIn,ULCOL);    % output layer  Nr
    
   [H(iROut).tstp,...
    H(iROut).period,...
    H(iROut).pertim,...
    H(iROut).totim,...
    H(iROut).label,...
    iLay,...         % the layer number of this record
    H(iROut).NROW,...
    H(iROut).NCOL,...
    values...        % entire layer array
    ]=contrRec(fp,bytes);
        
    H(iROut).values(:,:,iL)=values(userRows,userCols);
    H(iROut).lays=userLayers;
    H(iROut).rows=userRows;
    H(iROut).cols=userCols;
    if verbose
        if iL==1,  % everty time iL corresponds to the first user layer
            fprintf('iRecIn=%3d iRecOut=%3d periodsInFile=%3d, tstepsInFile%3d, pertim=%12g, totim=%12g Layers=%3d',...
                iRecIn,...
                iROut,...
                H(iROut).period,...
                H(iROut).tstp,...
                H(iROut).pertim,...
                H(iROut).totim,...
                length(H(iROut).lays));
        else
            if rem(10,iROut)==0, fprintf('.'); end
        end
        if iL==length(userLayers)
            fprintf('%d\n',iROut);
        end
    else
        fprintf('.');
        if rem(iROut,50)==0, fprintf('%d\n',iROut); end
    end
end

if ~verbose && rem(iROut,50)~=0, fprintf('\n'); end

fclose(fp);

end

function [kstp,kper,pertim,totim,label,nlay,nrow,ncol,values]=contrRec(fp,bytes)
% [kstp,kper,text,nlay,nrow,ncol,data]=contrRec(fp,bytes)
% --- reads a complete layer from a MODFLOW binary file
% --- if noread exists, only the info  record is read to safe time
% TO 070703 091214

if nargin==1, % just get the byte offset in the Binary file
    fread(fp, 2,'int32');     % kstp, kper
    fread(fp, 2,'float32');   % pertim, totim
    
    
    % FORTRAN bytes at beginning and end of each record, are recognized by
    % assuming they are not in the file and checking wheather the first 4
    % bytes contain any non readable charachters if so, we assume that each
    % record is preceded and ended with 4 bytes
    % If compiled with a Fortran compiler that adds a different number of
    % bytes than 0 or 4 this has to be adapted
    % TO 100502
    
    p=ftell(fp);
    label =char(fread(fp, 4,'uchar')');
    if any(label<'0' & label~=' ') || any(label>'9' & label<'A') || any(label>'Z' & label<'a') || any(label>'z')
        bytes=length(label);
    else
        bytes=0;
    end
    fseek(fp,p,'bof');

    fread(fp,bytes,'int8');

    label =char(fread(fp,16,'uchar')');
    
    ncol  =fread(fp, 1,'int32');
    nrow  =fread(fp, 1,'int32');
    nlay  =fread(fp, 1,'int32');

    fread(fp,bytes,'int8');
    fread(fp,bytes,'int8');

    n=ftell(fp); fread(fp,1,'float'); floatlen=ftell(fp)-n;
    fseek(fp,floatlen*(ncol*nrow-1),0);
    
    fread(fp,bytes,'int8');

    kstp=bytes;
    return
end


% if bytes are known, nargin>1

fread(fp,bytes,'int8');

kstp  =fread(fp, 1,'int32');
kper  =fread(fp, 1,'int32');

pertim=fread(fp, 1,'float32');
totim =fread(fp, 1,'float32');

label =char(fread(fp,16,'uchar')'); label=label(label~=' ');

ncol  =fread(fp, 1,'int32');
nrow  =fread(fp, 1,'int32');
nlay  =fread(fp, 1,'int32');

fread(fp,bytes,'int8');
fread(fp,bytes,'int8');

if nargout<9 % don't fetch
    n=ftell(fp); fread(fp,1,'float32'); floatlen=ftell(fp)-n;
    fseek(fp,floatlen*(ncol*nrow-1),0);
else
    values=permute(reshape(fread(fp,ncol*nrow,'float'),[ncol,nrow]),[2,1]);    
end

fread(fp,bytes,'int8');

end
