!****************************************************************************
!   pi_omp.f - compute pi by integrating f(x) = 4/(1 + x**2) from 0 to 1
!   This is an OpenMP program
!****************************************************************************
      PROGRAM pi_omp

        IMPLICIT NONE
        DOUBLE PRECISION     PI25DT, pi, h, sum, x
        REAL                 t(2), t_start, t_end
        INTEGER              n, i, numprocs
        PARAMETER(PI25DT = 3.141592653589793238462643D0)

        numprocs = OMP_GET_THREAD_NUM()
        WRITE(*,*) "OpenMP initialized with ", numprocs, " threads"

        DO

           WRITE(*,*) "Enter the number of intervals (0 to quit)"
           READ(*,*) n
           IF ( n<=0 ) EXIT

           t_start = ETIME(t)

           h = 1.0d0/n
           sum  = 0.0d0

!$omp      PARALLEL
!$omp&       DEFAULT (SHARED), PRIVATE (V)
!$omp&       REDUCTION(+:GSUM)
!$omp        DO
             DO i = 1, n
               x = h * (DBLE(i) - 0.5d0)
               sum = sum + 4.d0/(1.d0 + x*x)
             END DO
!$omp      END PARALLEL

           pi = h * sum

           WRITE(*,*) "Approximation for pi is ", pi
           WRITE(*,*) "The error is ", abs(pi - PI25DT)

           t_end = ETIME(t)
           WRITE(*,*) "Elapsed time is", t_end-t_start

        END DO

      END PROGRAM pi_omp