PROGRAM Definite_Integral
!-----------------------------------------------------------------------
! Program to approximate the integral of a function over the interval
! [A,B] using the trapezoidal method.  This approximation is calculated
! by the function subprogram Integral; the integrand, the interval of
! integration, and the # of subintervals are passed as arguments to
! Integral.  Identifiers used are:
!   A, B      : endpoints of interval of integration
!   Integral  : approximates integral of F on [A, B]
!   Integrand : the integrand
!   Number_of_subintervals : # of subintervals into which [A, B] is cut
!
! Input:  A, B, and Number_of_subintervals
! Output: Approximation to integral of F on [A, B]
!-----------------------------------------------------------------------

   IMPLICIT NONE
   REAL :: A, B, Integral

   INTERFACE
      REAL FUNCTION Integrand(X)
         REAL, INTENT(IN) :: X
      END FUNCTION Integrand
   END INTERFACE

   INTEGER :: Number_of_subintervals

   PRINT *, "Enter the interval endpoints and the # of subintervals:"
   READ *, A, B, Number_of_subintervals

   PRINT 10, Number_of_subintervals, &
             Integral(Integrand, A, B, Number_of_subintervals)
10 FORMAT (1X, "Trapezoidal approximate value using", I4, &
               " subintervals is", F10.5)

END PROGRAM Definite_Integral


! Integral ------------------------------------------------------------
! Function to calculate the trapezoidal approximation of the integral
! of the function F over the interval [A,B] using N subintervals.
! Local variables used are:
!     I       : counter
!     Delta_X : the length of the subintervals
!     X       : a point of subdivision
!     Y       : the value of the function  at X
!
!Accepts: Function F, endpoints A and B, and number N of subintervals
!Returns: Approximate value of integral of F over [A, B]
!-----------------------------------------------------------------------

REAL FUNCTION Integral(F, A, B, N)

   IMPLICIT NONE

   INTERFACE
      REAL FUNCTION F(X)
         REAL, INTENT(IN) :: X
      END FUNCTION F
   END INTERFACE

   REAL, INTENT(IN) :: A, B
   INTEGER, INTENT(IN) :: N
   REAL :: Delta_X,  X, Y
   INTEGER ::  I

!  Calculate subinterval length
!  and initialize the approximating sum and X

   Delta_X = (B - A) / REAL(N)
   X = A
   Integral = 0.0

!  Now calculate the approximating sum
   DO I = 1, N - 1
      X = X + Delta_X
      Y = F(X)
      Integral = Integral + Y
   END DO
   Integral = Delta_X * ((F(A) + F(B)) / 2.0 + Integral)

END FUNCTION Integral


! Integrand ------------------------------------------------------------
!                         The integrand
!-----------------------------------------------------------------------

REAL FUNCTION Integrand(X)

   REAL, INTENT(IN) :: X
   Integrand = X ** 2 + 1

END FUNCTION Integrand
