      PROGRAM laplace
*     implementation of Jacobi relaxation method to solve
*     Laplace's equation in rectangular geometry
      IMPLICIT NONE
      REAL*8 V(0:100,0:100), min_change
      INTEGER iterations, nx, ny
      CALL assign(V,nx,ny,min_change)
      CALL iterate(V,nx,ny,min_change,iterations)
      END
C
      SUBROUTINE assign(V,nx,ny,min_change)
*     linear dimension of rectangular region
      IMPLICIT NONE
      REAL*8 V0, V(0:100,0:100), min_change
      INTEGER x, y, nx, ny
      nx = 9                                ! number of interior points in x-direction
      ny = 9                                ! number of interior points in y-direction
      V0 = 10                               ! boundary potential of rectangle
      WRITE(*,'(A,$)') 'percentage change = '
      READ(*,*) min_change
      min_change = min_change/100
*     fix potential on boundary of rectangle
      DO x = 0, nx + 1
         V(x,0) = V0
         V(x,ny+1) = V0
      END DO
      DO y = 0, ny + 1
         V(0,y) = V0
         V(nx+1,y) = V0
      END DO
*     guess initial values of potentials of interior sites
      DO y = 1, ny
         DO x = 1, nx
            V(x,y) = 0.9D0*V0
         END DO
      END DO
      CALL show_output(V,nx,ny,0)
      END
C
      SUBROUTINE iterate(V,nx,ny,min_change,iterations)
      IMPLICIT NONE
      REAL*8 Vave(100,100), change, dV, V(0:100,0:100), min_change
      INTEGER x, y, nx, ny, iterations
      LOGICAL Ltmp1_
      iterations = 0
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
*        maximum difference of iterated potential at any site
         change = 0
         DO y = 1, ny
            DO x = 1, nx
*              average of potential of neighboring cells
               Vave(x,y) = V(x+1,y) + V(x-1,y)
               Vave(x,y) = Vave(x,y) + V(x,y+1) + V(x,y-1)
               Vave(x,y) = 0.25D0*Vave(x,y)
*              compute percentage change in potential
               IF(Vave(x,y) .NE. 0) THEN
                  dV = abs((V(x,y) - Vave(x,y))/Vave(x,y))
                  IF(dV .GT. change) change = dV
               END IF
            END DO
         END DO
         DO y = 1, ny
            DO x = 1, nx
               V(x,y) = Vave(x,y)
            END DO
         END DO
         iterations = iterations + 1
         CALL show_output(V,nx,ny,iterations)
         Ltmp1_ = change .GT. min_change
      END DO
      END
C
      SUBROUTINE show_output(V,nx,ny,iterations)
      IMPLICIT NONE
      REAL*8 V(0:100,0:100)
      INTEGER x, y, nx, ny, iterations
      WRITE(*,*)
      WRITE(*,*) 'iteration =', iterations
      DO y = 0, ny + 1
         DO x = nx + 1, 0, - 1
            WRITE(*, '(F7.2,$)') V(x,y)
         END DO
         WRITE(*,*)
      END DO
      END
