subroutine bndry(ne,xc,yc,zc,xn,h,g,l,u,q,area)
      real*8 h(l,l),g(l,l),xc(l),yc(l),zc(l),xn(l,3)
      real*8 u(l),q(l),area(l)
      real*8 a11,a22,a33,ux(500),theta
      integer ipvt(500),info,job
      open(25,file='tmp.dat',form='unformatted')
c     
      do k=1,1
         do i=1,ne
            theta=dacos(xc(i))
            q(i)=xn(i,k)
            ux(i)=-0.5*dcos(theta)
         enddo
c     
c        compute the right-hand side of the integral equation
c     
         do i=1,ne
            u(i)=0.0
            do j=1,ne
               u(i)=u(i)+g(i,j)*q(j)
            enddo
         enddo
         write(25) ne
         do i=1,ne
            write(25) (h(i,j),j=1,ne)
         enddo
         write(25) (u(i),i=1,ne)
         rewind 25
c     
c        solve the simultaneous equation with netlib
c     
         iaa=l
         call dgefa(h,iaa,ne,ipvt,info)
         call dgesl(h,iaa,ne,ipvt,u,job)
         do i=1,ne
            q(i)=u(i)
            dff=(q(i)-ux(i))/ux(i)
            write(11,101) i,q(i),ux(i),dff    
            write(*,*) i,q(i),ux(i),dff    
         enddo
 101     format(i5,3(e13.6,2x))
         as=0.0
         pi=acos(-1.)
         do i=1,ne
            as=as+q(i)*xn(i,k)*area(i)
         enddo
c
c        a11, a22 and a33 are added masses
c
         as=as*3./(2.*pi)
         if(k.eq.1) a11=as
         if(k.eq.2) a22=as
         if(k.eq.3) a33=as
      enddo
      return
      end