      PROGRAM LAPLACE
      INTEGER MAXX, MAXY
      PRINT *,'MAXX = (z.B. 64) '
      READ *,MAXX
      PRINT *,'MAXY = (z.B. 64) '
      READ *,MAXY
      CALL RELAX (MAXX,MAXY)
      STOP
      END

      SUBROUTINE RELAX (MAXX,MAXY)
      INTEGER MAXX, MAXY
      REAL F(MAXX,MAXY), DF(MAXX,MAXY), HF(MAXX,MAXY)
      LOGICAL CMASK(MAXX,MAXY)
      REAL FMAX
      INTEGER ITER
C
      CMASK = .FALSE.
      CMASK (2:MAXX-1,2:MAXY-1) = .TRUE.
      F = 2.
      F(:,MAXY) = 1.
      WHERE (CMASK) F = 0.
      ITER = 0
      FMAX = 1
      DF = 0.0
      DO WHILE (FMAX > 0.001)
         ITER = ITER + 1
         HF = 0.25 * (CSHIFT(F,1,1) + CSHIFT(F,1,-1)
     $              + CSHIFT(F,2,1) + CSHIFT(F,2,-1)) - F
         WHERE (CMASK)
            DF = HF
            F = F + DF
         ENDWHERE 
         DF = ABS(DF)
         FMAX = MAXVAL (DF)
         PRINT *,'Iteration ',ITER,'  Max = ',FMAX
      END DO
      END SUBROUTINE
 
