module precision_module
implicit none
integer, parameter :: long = selected_real_kind(15,307)
end module precision_module

module timing_module
implicit none
integer, dimension (8), private :: dt
real, private :: h, m, s, ms, tt
real, private :: last_tt

contains
  subroutine start_timing()
  implicit none
    call date_and_time(values=dt)
    print 100, dt(1:3), dt(5:8)
100 format (1x,i4,'/',i2,'/',i2,1x,i2,':',i2,':',i2,1x,i3)
    h = real(dt(5))
    m = real(dt(6))
    s = real(dt(7))
    ms = real(dt(8))
    last_tt = 60*(60*h+m) + s + ms/1000.0
  end subroutine start_timing

  subroutine print_date_and_time
  implicit none
    call date_and_time(values=dt)
    print 100, dt(1:3), dt(5:8)
100 format (1x,i4,'/',i2,'/',i2,1x,i2,':',i2,':',i2,1x,i3)
  end subroutine print_date_and_time

  subroutine print_hms
  implicit none
    call date_and_time(values=dt)
    print 100, dt(5:8)
100 format (1x,i2,':',i2,':',i2,1x,i3)
  end subroutine print_hms

  subroutine print_ms
  implicit none
    call date_and_time(values=dt)
    h = real(dt(5))
    m = real(dt(6))
    s = real(dt(7))
    ms = real(dt(8))
    tt = 60*(60*h+m) + s + ms/1000.0
    print 100, tt
100 format (1x,f14.3)
  end subroutine print_ms

  subroutine print_time_difference
  implicit none
    call date_and_time(values=dt)
    h = real(dt(5))
    m = real(dt(6))
    s = real(dt(7))
    ms = real(dt(8))
    tt = 60*(60*h+m) + s + ms/1000.0
    print 100, (tt-last_tt)
100 format (1x,f14.3)
    last_tt = tt
  end subroutine print_time_difference

  real function time_difference()
  implicit none
    tt = 0.0
    call date_and_time(values=dt)
    h = real(dt(5))
    m = real(dt(6))
    s = real(dt(7))
    ms = real(dt(8))
    tt = 60*(60*h+m) + s + ms/1000.0
    time_difference = tt - last_tt
  end function time_difference
end module timing_module

program ch2904
  use precision_module
  use timing_module
  use omp_lib
  implicit none
  real (long) :: fortran_internal_pi
  real (long) :: partial_pi
  real (long) :: openmp_pi
  real (long) :: width
  real (long) :: x
  integer :: nthreads
  integer :: i
  integer :: j
  integer :: k
  integer :: n

  nthreads = omp_get_max_threads()
  fortran_internal_pi = 4.0_long*atan(1.0_long)
  print *, ' Maximum number of threads is ', nthreads
  k = 1
  do
    call start_timing()
    n = 100000
    call omp_set_num_threads(k)
    print *, ' Number of threads = ', k
    do j = 1, 5
      width = 1.0_long/n
      partial_pi = 0.0_long
!$OMP parallel do private(x) shared(width)reduction(+:partial_pi)
      do i = 1, n
        x = width*(real(i,long)-0.5_long)
        partial_pi = partial_pi + f(x)
      end do
!$omp end parallel do
      openmp_pi = width*partial_pi
      print 20, n, time_difference()
20    format (' N intervals = ',i12,' time =',f8.3)
      print 30, openmp_pi, abs(openmp_pi-fortran_internal_pi)
30    format (' openmp_pi = ',f20.16,/,'difference = ',f20.16)
      n = n*10
    end do
    k = k*2
    if (k>nthreads) exit
  end do
contains
  real (long) function f(x)
    implicit none
    real (long), intent (in) :: x

    f = 4.0_long/(1.0_long+x*x)
  end function f
end program ch2904