RECURSIVE FUNCTION ADAPTIVE_QUAD (F, A, B, FA, FB, TOL, ABS_ERROR) &
                RESULT (RESULT)
IMPLICIT NONE

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

       REAL, INTENT(IN)        :: A, B, FA, FB, TOL
       REAL, INTENT(OUT)       :: ABS_ERROR
       REAL                    :: RESULT

       REAL                    :: STEP, MIDDLE_POINT, FMIDDLE
       REAL                    :: ONE_TRAPEZOIDAL_AREA, TWO_TRAPEZOIDAL_AREAS
       REAL                    :: LEFT_AREA, RIGHT_AREA
       REAL                    :: DIFF, ABS_ERROR_L, ABS_ERROR_R

       STEP = B-A
       MIDDLE_POINT= 0.5 * (A+B)
       FMIDDLE = F(MIDDLE_POINT)
       ONE_TRAPEZOIDAL_AREA = STEP * 0.5 * (FA + FB)
       TWO_TRAPEZOIDAL_AREAS = STEP * 0.25 * (FA + 2.0*FMIDDLE + FB)
       DIFF = TWO_TRAPEZOIDAL_AREAS - ONE_TRAPEZOIDAL_AREA

       IF ( ABS (DIFF) < TOL ) THEN
              RESULT = TWO_TRAPEZOIDAL_AREAS + DIFF/3.0
              ABS_ERROR = ABS(DIFF)
       ELSE
              LEFT_AREA = ADAPTIVE_QUAD (F, A, MIDDLE_POINT, FA, FMIDDLE, &
                      0.5*TOL, ABS_ERROR_L)
              RIGHT_AREA = ADAPTIVE_QUAD (F, MIDDLE_POINT, B, FMIDDLE, FB, &
                       0.5*TOL, ABS_ERROR_R)
              RESULT = LEFT_AREA + RIGHT_AREA
              ABS_ERROR = ABS_ERROR_L + ABS_ERROR_R
       END IF
END FUNCTION ADAPTIVE_QUAD