PROGRAM trap2
!
! Author:  Br. David Carlson
!
! Date:  January 1, 2000
!
! Revised:  April 22, 2004
!
! This program finds the approximate value of the integral of sin(x) from
! 0 to pi.  The trapezoidal rule is used, with double precision reals.

IMPLICIT NONE

INTEGER, PARAMETER::num = 100   ! Number of strips of area to add up.
REAL(KIND=8)::pi, delta, sum, area, x
INTEGER::k

! Function used:
REAL(KIND=8)::f

WRITE (*, *) 'Program to integrate sin(x) from 0 to pi'

pi = 2.0D0 * DASIN(1.0D0)  ! Double precision calculation to find pi.
delta = pi / num
x = delta
sum = (f(0.0D0) + f(pi)) / 2.0D0

DO k = 1, num - 1
   sum = sum + f(x)
   x = x + delta
END DO

area = delta * sum

WRITE (*, 100) area
100 FORMAT (1X, 'The trapezoidal rule gives an integral value of: ', E17.9)
WRITE (*, *) 'Of course, the exact value is 2.'

END PROGRAM


!  Given:   x   a real value
!  Task:    Compute f(x) and return it.
!  Return:  Computed value in function name.
REAL(KIND=8) FUNCTION f(x)

IMPLICIT NONE

REAL(KIND=8), INTENT(IN)::x

f = DSIN(x)   ! Double precision sine function
END FUNCTION