!-----------------------------------------------
! program to numerically integrate a function
! using the trapezoidal rule, reading parameter 
! from a input file
! 
! serial version
! 
! compile on bgq:
! module load profile/advanced front-end-xl/1.0
! xlf90 Trapez.f90 -o Trapez.x
!
! run on the front-end node
!-----------------------------------------------


  program integral
  implicit none

  real*8   :: f, a, b, delx, x, y, tot
  integer  :: n, i

  open(10,file='input')
  read(10,*)a ! starting point
  read(10,*)b ! eneding point
  read(10,*)n ! number of bins
  close(10)

  delx = (b-a)/real(n)
  x = a
  tot = 0.d0

  DO i = 1,n-1
    x = a + i * delx
    y= f(x)
    tot = tot + y
  END DO

  tot = delx * ((f(a)+f(b))/2.0 + tot)

  print*,"approximate integral using ",n," intervals =",tot

  end program

  function f(x)
  implicit none
  real*8 :: f,x
  f = sin(x)!!x**2+1
  end function