subroutine boundf(x,y,th,bl,br,bb,bt,np,nx,ny)
*     *******************************************************
*     obtain the normal velocity of the fluid at the boundaries.
*
*     trapezoidal rule is used to discretize the integral.
*     *******************************************************
      integer nx,ny,np,i,j,k
      double precision x(0:8192),y(0:8192),th(0:8192)
      double precision bl(ny),br(ny),bb(nx),bt(nx)
      double precision sa,hm,hinv,h,hp,xb,yb,dlg,rsq
      common/salpha/sa
      common/gridsz/h,hp
      common/meshsz/hm,hinv
      do 10 j=1,ny
         bl(j) = 0.D0
         br(j) = 0.D0
         yb = (dfloat(j)-0.5D0)*hm
         do 15 k=1,np
            rsq = x(k)**2 + (yb-y(k))**2
            dlg = x(k)
            bl(j) = bl(j) - dlg*dcos(th(k))/rsq
            rsq = (1.D0-x(k))**2 + (yb-y(k))**2
            dlg = 1.D0 - x(k)
            br(j) = br(j) - dlg*dcos(th(k))/rsq
 15      continue
         bl(j) = sa*bl(j)*h
         br(j) = sa*br(j)*h
 10   continue

      do 20 i=1,nx
         bb(i) = 0.D0
         bt(i) = 0.D0
         xb = (dfloat(i)-0.5D0)*hm
         do 25 k=1,np
            rsq = (xb-x(k))**2 + y(k)**2
            dlg = y(k)
            bb(i) = bb(i) - dlg*dcos(th(k))/rsq
            rsq = (xb-x(k))**2 + (1.D0-y(k))**2
            dlg = 1.D0 - y(k)
            bt(i) = bt(i) - dlg*dcos(th(k))/rsq
  25     continue
         bb(i) = sa*bb(i)*h
         bt(i) = sa*bt(i)*h
  20  continue

      end