program omp_sum
  use omp_lib
  implicit none

  integer, parameter :: dp = selected_real_kind(15,300)
  real(kind=dp):: start_time, stop_time
  real(kind=dp), allocatable, dimension(:) :: array
  real(kind=dp):: array_sum 
  integer :: i, N, ierror

  !set up program size
  N = 1000000
  array_sum = 0.0_dp

  !allocate dynamic storage
  allocate(array(1:N),stat=ierror)
  if (ierror/=0) stop 'Error in allocating array'

  !initialise array values
  do i=1,N
     array(i)=i
  end do

  !Now do the calculation - calculate the sum of the first N integers
  !$OMP parallel private(start_time, stop_time, i) shared(N,array, array_sum) default(none)
  start_time=omp_get_wtime()

  !$OMP do reduction(+:array_sum)
  do i=1,N
     array_sum=array_sum+array(i)
  end do
  !$OMP end do

  stop_time=omp_get_wtime()
  print *,'Elapsed time on thread num',omp_get_thread_num(),'=',stop_time-start_time
  !$OMP end parallel

  print *,'The sum of the first ',N,' integers is ',array_sum

  deallocate(array,stat=ierror)
  if (ierror/=0) stop 'Error in deallocating array'

end program omp_sum
