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 ch3003
  use precision_module
  use timing_module
  implicit none
  real (long) :: fortran_internal_pi
  real (long) :: partial_pi
  real (long) :: coarray_pi
  real (long) :: width
  real (long) :: total_sum
  real (long) :: x
  real (long) , codimension[*] :: partial_sum
  integer :: n_intervals
  integer :: i
  integer :: j
  integer :: current_image
  integer :: n_images

  fortran_internal_pi = 4.0_long*atan(1.0_long)
  n_images = num_images()
  current_image = this_image()
  if (current_image==1) then
    print *, ' Number of images = ', n_images
  end if
  n_intervals = 100000
  do j = 1, 5
    if (current_image==1) then
      call start_timing()
    end if
    width = 1.0_long/real(n_intervals,long)
    total_sum = 0.0_long
    partial_sum = 0.0_long
    do i = current_image, n_intervals, n_images
      x = (real(i,long)-0.5_long)*width
      partial_sum = partial_sum + f(x)
    end do
    partial_sum = partial_sum*width
    sync all
    if (current_image==1) then
      do i = 1, n_images
        total_sum = total_sum + partial_sum[ i]
      end do
      coarray_pi = total_sum
      print 20, n_intervals, time_difference()
20    format (' n intervals = ',i12,' time =',f8.3)
      print 30, coarray_pi, abs(coarray_pi-fortran_internal_pi)
30    format (' pi = ',f20.16,/,' difference = ',f20.16)
    end if
    n_intervals = n_intervals*10
    sync all
  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 ch3003