subroutine twodinit( a, b, f, nx, sx, ex, sy, ey )
      integer nx, sx, ex, sy, ey
      double precision a(sx-1:ex+1, sy-1:ey+1), b(sx-1:ex+1, sy-1:ey+1), &
                       f(sx-1:ex+1, sy-1:ey+1)
!
      integer i, j
!
      do j=sy-1,ey+1
         do i=sx-1,ex+1
            a(i,j) = 0.0d0
            b(i,j) = 0.0d0
            f(i,j) = 0.0d0
         enddo
      enddo
!
!    Handle boundary conditions
!
      if (sx .eq. 1) then
         do j=sy,ey
            a(0,j) = 1.0d0
            b(0,j) = 1.0d0
         enddo
      endif
      if (ex .eq. nx) then
         do j=sy,ey
            a(nx+1,j) = 0.0d0
            b(nx+1,j) = 0.0d0
         enddo
      endif
      if (sy .eq. 1) then
         do i=sx,ex
            a(i,0) = 1.0d0
            b(i,0) = 1.0d0
         enddo
      endif
!
      return
      end