c*******************************************************************
c     p1d-seq.f - a solution to the Poisson problem using Jacobi 
c     interation 
c*******************************************************************
      program main 
c     
      integer maxn
      parameter (nx=8,ny=8)
      double precision  a(0:nx+1,0:ny+1),b(0:nx+1,0:ny+1)
      double precision  f(0:nx+1,0:ny+1),diff,df
      integer nx,ny
      do j=0,ny+1
         do i=0,nx+1
            a(i,j) = 0.0d0
            b(i,j) = 0.0d0
            f(i,j) = 0.0d0
         enddo
      enddo
c     
c     Handle boundary conditions
c     
      do j=1,ny
         a(0,j) = 1.0d0
         b(0,j) = 1.0d0
         a(nx+1,j) = 0.0d0
         b(nx+1,j) = 0.0d0
      enddo
      do i=1,nx
         a(i,0) = 1.0d0
         b(i,0) = 1.0d0
      enddo
      h = 1.0d0/dble(nx+1)
      diff=10.
      iter=0
c
c     Jacobi Iteration
c
      do while(dabs(diff).gt.1.e-4)
c      do k=1,2

         do j=1,ny
            do i=1,nx
               b(i,j)=0.25*(a(i-1,j)+a(i,j+1)+a(i,j-1)+a(i+1,j))- 
     *              h*h*f(i,j)
            enddo
         enddo
         iter=iter+1
         df=0.
         do i=1,nx
            do j=1,ny
               df=df+(a(i,j)-b(i,j))**2
            enddo
         enddo
         write(*,*) 'iter=',iter,'   df=',df
         do i=0,nx+1
            write(*,*) (b(i,j),j=0,ny+1)
         enddo
         do i=1,nx
            do j=1,ny
               df=df+(a(i,j)-b(i,j))**2
            enddo
         enddo
         do j=1,ny
            do i=1,nx
               a(i,j)=0.25*(b(i-1,j)+b(i,j+1)+b(i,j-1)+b(i+1,j)- 
     *              h*h*f(i,j))
            enddo
         enddo
         iter=iter+1
         df=0.
         do i=1,nx
            do j=1,ny
               df=df+(a(i,j)-b(i,j))**2
            enddo
         enddo
         write(*,*) 'iter=',iter,'   df=',df
         diff=df
c         ip=1
c         if(ip.eq.1) stop
      enddo
      do i=1,nx
         write(*,*) (a(i,j),j=1,ny)
      enddo
      stop
      end