PROGRAM Definite_Integral_5
!-----------------------------------------------------------------------
! Program to approximate the integral of a function over the interval
! [A,B] using the trapezoidal method.  This approximation is calculated
! by the subroutine Integrate; the integrand, the interval of
! integration, and the # of subintervals are passed as arguments to
! Integrate.  Identifiers used are:
!   A, B      : endpoints of interval of integration
!   Integrate : subroutine to approximate 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

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

  INTEGER :: Number_of_Subintervals

  WRITE (*, '(1X, A)', ADVANCE = "NO") &
        "Enter the interval endpoints and the # of subintervals: "
  READ *, A, B, Number_of_Subintervals

  CALL Integrate(Integrand, A, B, Number_of_Subintervals)

CONTAINS

  ! Integrate ----------------------------------------------------------
  ! Subroutine 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
  !     DeltaX  : the length of the subintervals
  !     X       : a point of subdivision
  !     Y       : the value of the function  at X
  !     Sum     : the approximating sum
  !
  !Accepts: Function F, endpoints A and B, and number N of subintervals
  !Output:  Approximate value of integral of F over [A, B]
  !---------------------------------------------------------------------
  
  SUBROUTINE Integrate(F, A, B, N)
  
    IMPLICIT NONE
  
    REAL, INTENT(IN) :: A, B
    INTEGER, INTENT(IN) :: N
    REAL :: F, DeltaX, X, Y, Sum
    INTEGER ::  I
  
    ! Calculate subinterval length
    ! and initialize the approximating sum and X
  
    DeltaX = (B - A) / REAL(N)
    X = A
    Sum = 0.0
  
    ! Now calculate the approximating sum
    DO I = 1, N - 1
       X = X + DeltaX
       Y = F(X)
       Sum = Sum + Y
    END DO
    Sum = DeltaX * ((F(A) + F(B)) / 2.0 + Sum)
    PRINT 10, Number_of_Subintervals, Sum
    10 FORMAT (1X, "Trapezoidal approximate value using", I4, &
                   " subintervals is", F10.5)
  
  END SUBROUTINE Integrate

END PROGRAM Definite_Integral_5


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

REAL FUNCTION Integrand(X)

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

END FUNCTION Integrand