!****************************************************************************
!   pi_f77.f - compute pi by integrating f(x) = 4/(1 + x**2) from 0 to 1
!   This is a serial g77 compatible program
!****************************************************************************
      PROGRAM pi_f77

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

        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
           DO i = 1, n
             x = h * (DBLE(i) - 0.5d0)
             sum = sum + 4.d0/(1.d0 + x*x)
           END DO
           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_f77