program thrasher
  ! remember that size on file is not same as size in memory (it's about 4.5 times bigger on file)
  use omp_lib
  implicit none

  integer, parameter:: numFiles=32, mbVarSize=40
  real, parameter:: bytesMB = 1024.0**2
  character*8:: fileName                               ! fort.xxx
  integer:: i, myStatus
  double precision:: start

! vars for writing files
  integer, parameter:: realSize=4                                  ! 4 bytes for a real
  integer, parameter:: numReals=mbVarSize*(bytesMB/float(realSize))
  real,allocatable:: random(:), inputs(:)

! vars for par
  integer:: stream

! creat data
  allocate(random(numReals), stat=myStatus)
  if (myStatus /= 0) stop 'error allocating random'
  call random_number(random)
  write(*,*) 'have created',numReals,' reals. Size is',size(random)

!   ! write files
!   do i=1, numFiles
!      filename(1:5)="fort."
!      filename(6:8)=val2str(i)
!      write(*,*) 'opening file: ', filename
!      open(55, file=filename, status="new", form="formatted", action="write", iostat=myStatus)
!      if (myStatus /= 0) stop 'error opening new file'
!      write(55,*) random
!      write(*,*) '   file written'
!      close(55, iostat=myStatus)
!      if (myStatus /= 0) stop 'error closing file'
!      write(*,*) '   file closed successfully'
! end do

  deallocate(random)


!* HEC - there are various options: 
!* (a) single alloc pre-PAR - since this fails we do this last
!* (b) alloc per thread instantiation within PAR
!* (c) alloc with PAR but only once per thread

  
  ! read files SEQ
  allocate(inputs(numReals), stat=myStatus)
  if (myStatus /= 0) stop 'error allocating seq inputs'

  start=omp_get_wtime()
  do i=1, numFiles
     filename(1:5)="fort."
     filename(6:8)=val2str(i)
     write(*,*) 'opening file: ', filename
     open(55, file=filename, status="old", form="formatted", action="read", iostat=myStatus)
     if (myStatus /= 0) stop 'error opening existing file'
     read(55, *) inputs(1:numReals)
     write(*,*) '   file read: chkSum:', sum(inputs)
     close(55, iostat=myStatus)
     if (myStatus /= 0) stop 'error closing file'
     write(*,*) '   file closed successfully'
  end do
  write(*,*) 'seq reads took:', omp_get_wtime()-start
  deallocate(inputs)


!* HEC --------------------------------------------------------------------------------------
  ! read files in PARALLEL (b) alloc per thread instantiation within PAR
  start=omp_get_wtime()
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(filename, myStatus, stream, inputs)
  do i=1, numFiles
     stream=50+i
     filename(1:5)="fort."
     filename(6:8)=val2str(i)
     write(*,*) 'opening file: ', filename
     open(stream, file=filename, status="old", form="formatted", action="read", iostat=myStatus)
     if (myStatus /= 0) stop 'error opening existing file'

     allocate(inputs(numReals), stat=myStatus)
     if (myStatus /= 0) stop 'error allocating seq inputs'

     read(stream, *) inputs(1:numReals)
     write(*,*) '   file read: chkSum:', sum(inputs)

     deallocate(inputs)

     close(stream, iostat=myStatus)
     if (myStatus /= 0) stop 'error closing file'
     write(*,*) '   file closed successfully'
  end do
  write(*,*) 'PAR-(b) reads took:', omp_get_wtime()-start

!* HEC --------------------------------------------------------------------------------------
  ! read files in PARALLEL (c) alloc with PAR but only once per thread
  start=omp_get_wtime()
!$OMP PARALLEL DEFAULT(NONE) PRIVATE(filename, myStatus, stream, inputs) 
  allocate(inputs(numReals), stat=myStatus)
  if (myStatus /= 0) stop 'error allocating seq inputs'

!$OMP DO
  do i=1, numFiles
     stream=50+i
     filename(1:5)="fort."
     filename(6:8)=val2str(i)
     write(*,*) 'opening file: ', filename
     open(stream, file=filename, status="old", form="formatted", action="read", iostat=myStatus)
     if (myStatus /= 0) stop 'error opening existing file'

     read(stream, *) inputs(1:numReals)
     write(*,*) '   file read: chkSum:', sum(inputs)

     close(stream, iostat=myStatus)
     if (myStatus /= 0) stop 'error closing file'
     write(*,*) '   file closed successfully'
  end do
!$OMP END DO NOWAIT
  deallocate(inputs)
!$OMP END PARALLEL 
  write(*,*) 'PAR-(c) reads took:', omp_get_wtime()-start

!  deallocate(inputs)

!* HEC --------------------------------------------------------------------------------------
  ! read files in PARALLEL (a) single alloc pre-PAR
  start=omp_get_wtime()
  allocate(inputs(numReals), stat=myStatus)
  if (myStatus /= 0) stop 'error allocating par-(a) inputs'

!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(filename, myStatus, stream, inputs)
  do i=1, numFiles
     stream=50+i
     filename(1:5)="fort."
     filename(6:8)=val2str(i)
     write(*,*) 'opening file: ', filename
     open(stream, file=filename, status="old", form="formatted", action="read", iostat=myStatus)
     if (myStatus /= 0) stop 'error opening existing file'
     read(stream, *) inputs(1:numReals)
     write(*,*) '   file read: chkSum:', sum(inputs)
     close(stream, iostat=myStatus)
     if (myStatus /= 0) stop 'error closing file'
     write(*,*) '   file closed successfully'
  end do
  write(*,*) 'PAR-(a) reads took:', omp_get_wtime()-start

  deallocate(inputs)

  write(*,*) 'PROGRAM COMPLETE'

contains
  character*(3) function val2str(n)
    implicit none
    integer, intent(in):: n
    ! pad with zeros if required
    write(val2str,'(I3.3)') n
  end function val2str

end program
