subroutine reflec(pf,nx,ny)
*     ***********************************************
*     impose the reflection boundary condition for pf
*     ***********************************************
      integer i,j,nx,ny
      double precision pf(0:nx+1,0:ny+1)
      do 10 i=1,nx
         pf(i,0) = pf(i,1)
         pf(i,ny+1) = pf(i,ny)
 10   continue
      do 20 j=1,ny
         pf(0,j) = pf(1,j)
         pf(nx+1,j) = pf(nx,j)
 20   continue
      end