Listing 3. Optimized Code

      poiint = 0.0
      do it=1, ntetra
      ch = tetvol(it)*(p(itetra(1,it)) +
     &     p(itetra(2,it)) + p(itetra(3,it))
     &     + p(itetra(4,it)))
      poiint = poiint +
     &     ch/sqrt((tcentr(1,it)-x(1))**2
     &     + (tcentr(2,it)-x(2))**2 +
     &     (tcentr(3,it)-x(3))**2)
      enddo
      poiint = 0.25*poiint