Replicable Monte Carlo simulation with OpenMP

Help
18

But it becomes reproducible if I insert the following at the start of the openmp region, before the loop.

if(omp_get_thread_num() /= 0) call random_number(x1)
!$omp barrier

Is it possible that the master thread’s seed is sometimes updated by random number generation, before all other threads generate a random number (and thus determine their seed)? The above code is designed to prevent that.

Full code below

program main
   use omp_lib
   implicit none
   logical verbose
   integer :: n1, n2, i, n, nt
   integer, allocatable :: seed(:)
   real, allocatable :: threadval(:)
   real(8) :: pi, x1, x2

   verbose = .true.
   if (command_argument_count() /= 0) verbose = .false.
 
   ! Initialize the master state to a repeatable sequence.
   call random_init(.true., .false.)

   n1 = 1000*1000

   !$omp parallel
   nt = omp_get_num_threads()
   !$omp end parallel

   allocate(threadval(nt))

   !$omp parallel private(x1, x2, i), shared(n1, threadval), reduction(+:n2)
   n2 = 0

   if(omp_get_thread_num() /= 0) call random_number(x1)
   !$omp barrier

   !$omp do
   do i = 1, n1
      call random_number(x1)
      call random_number(x2)
      if (x1**2 + x2**2 <= 1.0_8) n2 = n2 + 1
   end do
   !$omp end do

   threadval(omp_get_thread_num()+1) = x2  ! Store final RN on each thread

   !$omp end parallel

   pi = 4.0_8*dble(n2)/dble(n1)
   if (verbose) print '(A,F15.12,A,I0)', 'pi = ', pi, ', nt = ', nt

   !Check that threads are all distinct
   if (.not. verbose) print '(A,*(F9.6))', 'RN:', threadval
end program

EDIT – actually the above program is still not totally repeatable. I can get different values of threadval on repeated runs. But pi seems repeatable in any case (over hundreds of test runs).