module random_int_module
public :: random_int_array
contains
subroutine random_int_array(VALUES, low, high)
   integer, dimension (:), intent (out) :: VALUES
   integer, intent (in) :: low, high
   real, dimension (size (VALUES)) :: B
   call random_number(B)
   VALUES = int ((high - low + 1) * B + low)
end subroutine random_int_array
end module random_int_module
 
program probability_with_arrays
use random_int_module 
integer, parameter :: no_of_simulations = 100000
integer :: i ,wins
integer, dimension(10, no_of_simulations) :: RESULTS
! Initialize RESULTS
RESULTS=0
! Fill in the outcome of 5 coin tosses.
call random_int_array(RESULTS(1,:), 0, 1)
call random_int_array(RESULTS(2,:), 0, 1)
call random_int_array(RESULTS(3,:), 0, 1)
call random_int_array(RESULTS(4,:), 0, 1)
call random_int_array(RESULTS(5,:), 0, 1)
! Find the number of heads in the these 5 throws.
RESULTS(6,:) = count(RESULTS(1:5,:)==1,dim=1)
! Fill in the outcome of three dice throws
call random_int_array(RESULTS(7,:), 1, 6)
call random_int_array(RESULTS(8,:), 1, 6)
call random_int_array(RESULTS(9,:), 1, 6) 
! Delete outcome of the third die where less than 4 heads
! showed up.
where (RESULTS(6,:) < 4)
     RESULTS(9,:)=0
end where  
RESULTS(10,:) = sum(RESULTS(7:9,:),dim=1)
wins= count(RESULTS(10,:) > 10) 
print*,"Probability = ",real(wins)/no_of_simulations*100,"%" 
end program probability_with_arrays