Thursday, February 3, 2011

Read Modflow Head file using Matlab

When I first started working with Modflow, I used to use a utility program to convert the output data from binary to ascii. However, I have writen a function in Matlab that extracts the head field directly from the binary file to matlab workspace.

function H=readheadfile(filename,NR,NC,NL)
%NR, NC, NL: number of rows, columns, layers
fid=fopen(filename,'r');
A=fread(fid,'float'); %reads the entire file and stores it to the variable A
fclose(fid);
A=reshape(A,11+NC*NR,NL);
%{
The constant number 11 depends on the preprocessor platform. For Groundwater vistas, version 4 and 5, this number should be 11. For an old version of PmWin I had to change that to 14.
%}
A(1:11,:)=[];
for i=1:NL
H(:,:,i)=reshape(A(:,NL+1-i),NC,NR)';
end

Note that the output of the function readheadfile is a 3D array where the top layer corresponds to the bottom layer of the aquifer and vice versa.