C********************************************************************
C--- HOMEWORK #5
C--- Numericl Solution to Partial Differential Equations.
C---
C--- Solve Laplace's equation in 2 dimension using the 
C--- point-Jacobi iteration.
C--------------------------------------------------------------------
C--- Programmed by Sung-Hoon Ko(shko@npac.syr.edu)
C---
C********************************************************************
 
C--- Note Fortran77 version

      PROGRAM  HW5_LAPLACE_JACOBI
      IMPLICIT NONE
      REAL    GRID(1:27,1:27)
      REAL    eps
      INTEGER i, j, iter

C---  Install boundary values
      DO i=1, 27
         GRID(i,1)=1
         GRID(1,i)=1
         GRID(i,27)=1
         GRID(27,i)=1
      ENDDO

C--- Initialize interior points
      DO i=2,26
         DO j=2,26
            GRID(i,j)=0
         ENDDO
      ENDDO

      iter=0
      eps=1.0
      DO WHILE(eps .GT. 0.01)
         CALL JOCOBI(GRID, eps)
         iter=iter+1
      ENDDO

      WRITE(6,11) 'Fortran77 version :'
      WRITE(6,22) 'Error =',eps 
      WRITE(6,33) 'Number of iterations =',iter
 11   FORMAT(//,5X,A)
 22   FORMAT(//,8X,A,F9.6)
 33   FORMAT(/,8X,A,I3,//)

      STOP
      END



C=====================================================================
C---
C---  Perform Point-Jocobi on grid points G.
C---  Returns updated functions, writes error estimate in maxError.

      SUBROUTINE JOCOBI(G, maxError)
      REAL    G(1:27,1:27)   , Gnew(1:27,1:27)
      REAL    Diff(1:27,1:27), MaxError
      INTEGER i, j

      max = ABS(G(2,2))
      DO i=2, 26
         DO j=2, 26
            Gnew(i,j) = 0.25*(  G(i,j-1) 
     *                        + G(i,j+1)
     *                        + G(i-1,j) 
     *                        + G(i+1,j) )

            Diff(i,j) = ABS( Gnew(i,j) - G(i,j) )
         ENDDO
      ENDDO
      
      maxError=Diff(2,2)
      DO i=2,26
         DO j=2,26
            G(i,j)=Gnew(i,j)    ! Install updated Grid values
C---------- Find the maximum error 
            IF(maxError .LT. Diff(i,j)) THEN  
               maxError=Diff(i,j)
            ENDIF
         ENDDO
      ENDDO

      RETURN
      END
      



