Read and process a Fortran90 binary file

Hi,

I have a Fortran90 code where a variable called 'stream' is defined such that it stores certain variable every nstep:

open(19,file='stream',form='UNFORMATTED',position='append')
if(mod(NSTEP,100).eq.0)then
    write(19)NSTEP,U,W,zeta
    endif

How can I properly read the file in Matlab for postprocessing?

FYI, the fortran code is built with double precision.

Thank you

[James Tursa]

Unformatted Fortran can contain beginning and end of line markers in the file. You need to give us a small sample file to work with to verify this, and also tell us the type and dimensions of the variables NSTEP, U, W, and zeta.

[Walter Roberson]

The FORTRAN runtime system embeds the record boundaries in the data by inserting an INTEGER*4 byte count at the beginning and end of each unformatted sequential record during an unformatted sequential WRITE. The trailing byte count enables BACKSPACE to operate on records.

https://docs.oracle.com/cd/E19957-01/805-4939/6j4m0vnc4/index.html

[Carola Forlini]

I am attaching an extract of the file.

NSTEP is an integer = 2000

U, W, zeta are matrices = 32x512

Thank you

Answer 1

[Les Beckham]

It would help to have a sample file to experiment with.

Nevertheless, assuming that the write call that you have shown is the only thing that writes to this file in your Fortran code, I would suggest something like this using fread and reshape:

fid = fopen('your_file_name.ext', 'rb');
A = fread(fid);
data = reshape(A, 4, numel(A)/4)';
nstep = A(:,1);
U = A(:,2);
W = A(:,3);
zeta = A(:,4);

Note that you might run into endian-ness issues if your Fortran code runs on a different architecture than your Matlab is running on. See the documentation for the machinefmt option to fread if necessary.

[Walter Roberson]

form='UNFORMATTED' has record length markers before and after each record.

[Walter Roberson]

If you start with fread(fid, 1, 'integer*4') and discard that, and thereafter fread() with a skip field of 8 then it should be properly positioned to fread() more.

Unfortunately we are not told the size of NSTEP, U, W, zeta

[Carola Forlini]

Hi, I provided a sample file and the dimensions of NSTEP, U, W, zeta in the post above.

Thank you.

[Les Beckham]

Thanks for adding the example.

Can you clarify what you mean by this?

NSTEP is an integer = 2000

Especially the "= 2000" part. Is that just an example of a possible value? What is the size of this integer (number of bits)?

Also, are U, W, and zeta double precision (64 bit) floating point 2-d arrays with the given size?

[Carola Forlini]

I apologize for the misunderstanding.The size of NSTEP is 32 bit.

U, W, zeta are single precision (32 bit) of size N1×N3, with N1=32 and N3=512.

Answer 2

[Les Beckham]

In light of the additional variable size information and sample file provided by the OP, and the information from Walter about Fortran inserting record boundary byte counts at the beginning and end of each record written with an unformatted sequential WRITE, here is a fresh answer.

Note that this makes an assumption about the order the data is written (column major order). If this assumption is not correct, the inner for loop will need to be modified accordingly.

Hopefully this is close enough to get you where you need to go.

filenames = unzip('stream.zip')

filenames = 1x1 cell array
    {'stream'}

fid = fopen(filenames{1}, 'rb');

% start/stop byte counts (4 bytes each) + NSTEP (4 bytes) + U, W, zeta (3*32*512*4 bytes)
recordSize = 2*4 + 4 + 3*32*512*4; 

fseek(fid, 0, 'eof'); % seek to the end of the file
nBytes = ftell(fid); % find number of bytes in the file
frewind(fid); % go back to the beginning of the file

nRecords = nBytes / recordSize;

% preallocate
nStep = zeros(nRecords, 1);      % column vector 
U     = zeros(3, 512, nRecords); % 3D arrays
W     = zeros(3, 512, nRecords);
zeta  = zeros(3, 512, nRecords);

for iRecord = 1:nRecords
    fseek(fid, 4, 'cof'); % skip begin record byte count field
    nStep(iRecord) = fread(fid, 1, 'integer*4');
    for iRow = 1:3
        U(iRow, :, iRecord)    = fread(fid, 512, 'real*4');
        W(iRow, :, iRecord)    = fread(fid, 512, 'real*4'); 
        zeta(iRow, :, iRecord) = fread(fid, 512, 'real*4');
    end
    fseek(fid, 4, 'cof'); % skip end record byte count field
end

fclose(fid); % close the file

whos nStep U W zeta % check variable sizes

  Name        Size                 Bytes  Class     Attributes

  U           3x512x21            258048  double              
  W           3x512x21            258048  double              
  nStep      21x1                    168  double              
  zeta        3x512x21            258048  double              

plot(nStep); % make a plot see if nStep makes sense
grid on

Well, I don't know what your NSTEP should look like, but I would have expected a monotonically increasing count. So, either I made a mistake, or there is some additional information needed. It does appear that the size information you gave is either correct, or just coincidentally came out to be an integer divisor of the file size.