module integrate_rc_mod
  implicit none

  type:: internal_state
    integer:: position
    integer:: i, n
    real::    h, sum
  end type internal_state
  

contains

  subroutine integrate_rc(state,iflag,x,fx,a,b,n,val)
    type(internal_state), intent(inout):: state
    integer, intent(inout):: iflag
    integer, intent(in)   :: n
    real, intent(inout)   :: x, val
    real, intent(in)      :: a, b, fx
  
    if ( iflag == 0 ) then
       ! Initialization
       state%position = 1
       state%n = n
       state%h = (b-a)/n
       x = a        ! evaluate function at a
       iflag = 1
    else if ( iflag == 1 ) then
       select case ( state%position )
       case ( 1 ) ! position 1
          state%sum = 0.5*fx
          x = b      ! evaluate function at b
          iflag = 1
          state%position = 2
       case ( 2 ) ! position 2
          state%sum = state%sum + 0.5*fx
          state%i = 1      ! start of loop
          x = a + state%i*state%h
          iflag = 1
          state%position = 3
       case ( 3 ) ! position 3
          state%sum = state%sum + fx
          state%i = state%i + 1
          x = a + state%i*state%h
          iflag = 1
          if ( state%i == n-1 ) then
             state%position = 4
          end if
       case ( 4 ) ! position 4
          state%sum = state%sum + fx
          val = state%sum*state%h
          iflag = 0    ! That's all folks!
       case default
          iflag = -1   ! Error !
       end select
    else
       iflag = -1     ! Error !
    end if
  end subroutine integrate_rc
  

  subroutine integrate(f,a,b,n,val)
    real, intent(in) :: a, b
    real, intent(out):: val
    real, external   :: f     ! function
    integer          :: i, n
    real             :: h, sum
    h = (b-a)/n
    ! position 1
    sum = 0.5*f(a)
    ! position 2
    sum = sum + 0.5*f(b)
    do i = 1, n-1
      ! position 3
      sum = sum + f(a+i*h)
    end do
    ! position 4
    val = sum*h;
  end subroutine integrate
  

  subroutine integrate2(f,a,b,n,val)
    real :: a, b, val
    real :: myparam = 37.0
    integer :: n
    real, external :: f
    integer:: iflag
    real   :: x, fx
    type(internal_state):: state
    
    iflag = 0
    call integrate_rc(state,iflag,x,fx,a,b,n,val)
    do while ( iflag /= 0 )
      if ( iflag == 1 ) then
        ! code for evaluating fx = f(x)
        fx = f(myparam,x) ! or something horribly complicated...
        call integrate_rc(state,iflag,x,fx,a,b,n,val)
      else ! if ( iflag /= 0 ) then
        exit      ! exit do while loop
      end if
    end do
    ! val is now the computed integral
    
  end subroutine integrate2

  real function myfunc1(x)
    real :: x
    myfunc1 = x**2 * cos(x)
  end function myfunc1

  real function myfunc2(param,x)
    real :: param, x
    myfunc2 = myfunc1(x)
  end function myfunc2

  subroutine test_rc
    real :: pi, val1, val2
    integer :: n
    n = 10
    pi = 4*atan(1.0)
    val1 = 0.0
    call integrate(myfunc1,0.0,pi/2,n,val1)
    print *, 'Result of trapezoidal rule =', val1
    val2 = 0.0
    call integrate2(myfunc2,0.0,pi/2,n,val2)
    print *, 'Result of trapezoidal rule with RC =', val2
    print *, 'Difference = ', val1-val2
  end subroutine test_rc

end module integrate_rc_mod