Skip to content

Instantly share code, notes, and snippets.

@milancurcic
Last active June 12, 2020 09:58
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save milancurcic/6590d74fbbd3e765832b823474623a27 to your computer and use it in GitHub Desktop.
Save milancurcic/6590d74fbbd3e765832b823474623a27 to your computer and use it in GitHub Desktop.
Example Fortran program that reads a GrADS binary file
program read_grads_example
! Example program that reads a flat binary file
!
! With GNU Fortran, compile as:
!
! gfortran -fconvert=big-endian read_grads_example.f90
!
! With Intel Fortran, compile as:
!
! ifort -convert big_endian -assume byterecl read_grads_example.f90
! use portable 4-byte reals
use iso_fortran_env, only: real32
implicit none
character(len=*), parameter :: filename = 'taucw.grd' ! file to read
integer, parameter :: im = 2880, jm = 1440, nm = 60 ! grid dimensions
integer :: fileunit
integer :: record_length
real(kind=real32) :: field(im,jm,nm) ! 3-d array to store data
record_length = storage_size(field) / 8 * size(field)
open(newunit=fileunit, file=filename, access='direct', recl=record_length)
read(unit=fileunit, rec=1) field
close(fileunit)
write(*,*) 'Min / max: ', minval(field), maxval(field)
end program read_grads_example
@geoqueen
Copy link

geoqueen commented May 1, 2018

Worked great! Thank you :)

@robetatis
Copy link

Great contribution!
One short question on line 25: I don't see the declaration of storage_size. What is it? And why is record_length = storage_size(field) / 8 * size(field)?
Best
Roberto

@milancurcic
Copy link
Author

Thank you @robetatis. I'm happy you found it useful.

I don't see the declaration of storage_size. What is it?

storage_size is a built-in function that returns the size of its argument in bits. Pay special attention to "The result value is the size expressed in bits for an element of an array that has the dynamic type and type parameters of A", from the description on that page. Even if you pass an array to storage_size, the result is the number of bits occupied by one element of that array.

And why is record_length = storage_size(field) / 8 * size(field)?

Here the record length is expressed in bytes. As storage_size returns the number of bits and there are 8 bits in a byte, storage_size(field) / 8 is the storage size of one element in bytes, and then storage_size(field) / 8 * size(field) is the storage size of the whole array in bytes. We store one array per record, so that's our record length.

Note that by default gfortran assumes the value of recl in bytes, whereas ifort assumes it in number of elements. This is why we specify -assume byterecl flag with ifort--to have a portable solution between two compilers.

@robetatis
Copy link

Thanks @milancurcic !

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment