C	  Integration par la methode de Romberg
C	  Dominique Lefebvre  fevrier 2007
C     Inspire du code de A.Garcia dans Numerical Methods for Physics
C
C     Modifiee le 30/08/09 suite ŕ une erreur d'estimation de la precision
C     signalee par Maelle Nodet (Universite de Grenoble)
C  
C     a = borne inferieure d'integration
C     b = borne superieure d'integration
C     eps = précision souhaitee
C     aire = surface retournee

		SUBROUTINE IntRomberg(fn,a,b,eps,aire)

      	REAL fn,a,b,eps,aire
		EXTERNAL fn

        INTEGER out, pieces, i, nx(16)
      	INTEGER nt, ii, n, nn, l, ntra, k, m, j, l2
      	REAL h,xi, delta, sum
      	REAL c, fotom, t(136)
        LOGICAL fini
       
C     Calcul du premier terme
      	h = b-a
		pieces  = 1
      	nx(1) = 1
      	delta = h / pieces
      	c = (fn(a) + fn(b)) / 2.0
      	sum = c
      	t(1) = delta * c
      	n = 1
      	nn = 2

     	fini =.false.
        DO WHILE (.not. fini)
		  n = n + 1
       	  fotom = 4.0
       	  nx(n) = nn
          pieces = pieces * 2
       	  l = pieces - 1
      	  l2 = (l+1) / 2
       	  delta = h / pieces

C	   Evaluation par la methode des trapezes	
          DO ii = 1, l2
            i = ii * 2 - 1
            xi = a + delta * i
            sum = sum + fn(xi)
          ENDDO
      
          t(nn) = sum * delta
          ntra = nx(n-1)
          k = n-1

C	  Calcul de la nieme colonne du tableau	
          DO m = 1, k
            j = nn + m
            nt = nx(n-1) + m - 1
            t(j) = (fotom * t(j-1) - t(nt)) / (fotom - 1)
            fotom = fotom * 4
          ENDDO

C     Estimation de la precision
          IF (n .lt. 5) THEN
            nn = j+1
          ELSE IF (t(nn+1) .eq. 0.0) THEN
            nn = j+1
          ELSE
            IF (abs(t(ntra+1)-t(nn+1)) .le. abs(t(nn+1)*eps)) THEN 
              fini = .true.
            ELSEIF (abs(t(nn-1) - t(j)) .le. abs(t(j)*eps)) THEN 
               fini = .true.
            ELSE
              nn = j+1
            ENDIF
		  ENDIF	

       ENDDO

C     On retourne la valeur finale de l'aire
       aire = t(j)

       RETURN
	   END