subroutine evapot(x,y,th,vf,bl,br,bb,bt,
     +                  padd,uadd,vadd,nx,ny,np)
*     ****************************************************
*     Add the potential flow to satisfy the no-flux condition
*     on the boundary. The normal velocity of the front is
*     modified.
*     ****************************************************
      integer nx,ny,np,i,j,k,ii,jj
      double precision x(0:8192),y(0:8192),th(0:8192),
     +                 vf(0:8192)
      double precision padd(0:nx+1,0:ny+1),uadd(0:nx,ny),
     +                 vadd(nx,0:ny)
      double precision bl(ny),br(ny),bb(nx),bt(nx)
      double precision hm,hinv,xp,yp,dx,dy,up,vp,coef1,sum
      common/darcyc/coef1
      common/meshsz/hm,hinv
*     determine the velocity correction:
      do 10 j=1,ny
      do 10 i=1,nx-1
         uadd(i,j) = (padd(i,j) - padd(i+1,j))*hinv
  10  continue
      do 15 j=1,ny
         uadd(0,j) = bl(j)
         uadd(nx,j) = -br(j)
  15  continue
      do 20 j=1,ny-1
      do 20 i=1,nx
         vadd(i,j) = (padd(i,j) - padd(i,j+1))*hinv
  20  continue
      do 25 i=1,nx
         vadd(i,0) = bb(i)
         vadd(i,ny) = -bt(i)
  25  continue

*     bilinear interpolation to obtain the velocity at
*     front positions:

      do 30 k=1,np
         xp = x(k)
         yp = y(k)
*     for u:
         ii = int(xp*hinv)
         jj = int(yp*hinv+0.5D0)
         dx = xp*hinv - dfloat(ii)
         if(jj .eq. 0) then
            up = uadd(ii,1)*(1.D0-dx) +
     +           uadd(ii+1,1)*dx
         elseif(jj .eq. ny) then
            up = uadd(ii,ny)*(1.D0-dx) +
     +           uadd(ii+1,ny)*dx
         else
            dy = yp*hinv - dfloat(jj) + 0.5D0
            up = uadd(ii,jj)*(1.D0-dx)*(1.D0-dy) +
     +           uadd(ii+1,jj)*dx*(1.D0-dy) +
     +           uadd(ii,jj+1)*(1.D0-dx)*dy +
     +           uadd(ii+1,jj+1)*dx*dy
         endif
*     for v:
         ii = int(xp*hinv+0.5D0)
         ii = int(yp*hinv)
         dy = yp*hinv - dfloat(jj)
         if(ii .eq. 0) then
            vp = vadd(1,jj)*(1.D0-dy) +
     +           vadd(1,jj+1)*dy
         elseif(ii .eq. nx) then
            vp = vadd(nx,jj)*(1.D0-dy) +
     +           vadd(nx,jj+1)*dy
         else
            dx = xp*hinv - dfloat(ii) + 0.5D0
            vp = vadd(ii,jj)*(1.D0-dx)*(1.D0-dy) +
     +           vadd(ii+1,jj)*dx*(1.D0-dy) +
     +           vadd(ii,jj+1)*(1.D0-dx)*dy +
     +           vadd(ii+1,jj+1)*dx*dy
         endif
*     add to the single layer source velocity
         vf(k) = vf(k) - coef1*up*dsin(th(k)) + 
     +                   coef1*vp*dcos(th(k))
  30  continue
      vf(0) = vf(np)

      sum = 0.D0
      do 40 k=1,np
         sum = sum + vf(k)
 40   continue
      do 50 k=0,np
         vf(k) = vf(k) - sum/dfloat(np)
 50   continue

      end