subroutine Simpsons(points, ref, PM_S)
            REAL, intent(in), dimension(1:800,1:2) :: points
            REAL, intent(in):: ref
            REAL, intent(out):: PM_S
            INTEGER :: I
            PM_S=0
            DO I=1,800-1
                PM_S= PM_S + ( ref - points(I,1) ) * ( ( points(I+1,1) &
                  & - points(I,1) ) / 6 ) * ( points(I,2) + 2 * &
                  & ( points(I,2) + points(I+1,2) ) + points(I+1,2) )
            END DO
            ! Assume linear interpolation for the f(A+B/2) since we don’t
            ! have a point between A and B
        END subroutine Simpsons
   END PROGRAM Pitching_Moment_Cp