poiint = 0.0 do it=1, ntetra i1 = itetra(1,it) i2 = itetra(2,it) i3 = itetra(3,it) i4 = itetra(4,it) cx = 0.25*(vertex(1,i1)+vertex(1,i2)+ & vertex(1,i3)+vertex(1,i4)) cy = 0.25*(vertex(2,i1)+vertex(2,i2)+ & vertex(2,i3)+vertex(2,i4)) cz = 0.25*(vertex(3,i1)+vertex(3,i2)+ & vertex(3,i3)+vertex(3,i4)) ch = 0.25*tetvol(it)*(p(i1) + p(i2) + & p(i3) + p(i4)) poiint = poiint+ch/sqrt((cx-x(1))**2+ (cy-x(2))**2+(cz-x(3))**2) enddo