C     ********************************************
C     ** TWO-DIMENSIONAL BOUNDARY VALUE PROBLEM **
C     **                                        **
C     ** JACOBI METHOD                          **
C     **                                        **
C     ** 2000/01/29              Pyon           **
C     ********************************************
C
C
      PARAMETER (MX=301,MY=301)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION T(MX,MY),WK(MX,MY)
C
C     ***** NUMBER OF FIELD DIVIDE *****
C
        NX  =51
        NY  =51
        NXM1=NX-1
        NYM1=NY-1
C
C     ***** DIVIDED LENGTH *****
C
        DX  =1.0D0/DFLOAT(NXM1)
        DY  =1.0D0/DFLOAT(NYM1)
        DX2I=1.0D0/DX/DX
        DY2I=1.0D0/DY/DY
C
C     **** ITERATION CONSTANT FOR LAPLACE EQ.*****
C
        ITER1 =5000
        EPS1  =1.0D-6
C
C   新しいFORTRANならばいらないが，古いものを使っている人のため
C
      DO 100 J=1,NY
      DO 100 I=1,NX
C
  100    T(I,J) =0.0D0
C
C     *** BOUNDARY CONDITION ON TOP WALL ***
C
      DO 110 I=1,NX
C
  110    T(I,NY)=1.0D0
C
C     ***** SUCCESSIVE OVER RELAXATION LOOP *****
C
         COEF=1.D0/(2.D0*(DX2I+DY2I))
C
      DO 300 LOOP=1,ITER1
C
         ERR=0.0D0
C
C     ***** CAL. INTERNAL POINT *****
C
      DO 120 J=2,NYM1
      DO 120 I=2,NXM1
C
         TO=(T(I+1,J)+T(I-1,J))*DX2I
     &     +(T(I,J+1)+T(I,J-1))*DY2I
C
         DR=TO*COEF-T(I,J)
C
         WK(I,J)=T(I,J)+DR
C
         ERR=ERR+DR*DR
C
  120 CONTINUE
C
      DO 130 J=2,NYM1
      DO 130 I=2,NXM1
C
  130   T(I,J)=WK(I,J)
C
        ERR=DSQRT(ERR/(DFLOAT((NX-2)*(NY-2))))
        IF(ERR.LT.EPS1) GOTO 90
C
  300 CONTINUE
      WRITE(*,*)(' NOT CONVERGED!')
   90 CONTINUE
C
      WRITE(*,600) LOOP,ERR
  600 FORMAT(' LOOP & ERROR=',I5,E12.5)
C
C     *** OUTPUT FILE ***
C
      WRITE(3,610) NX,NY
  610 FORMAT(2I5)
C
      WRITE(3,620)((T(I,J),I=1,NX),J=1,NY)
  620 FORMAT(F12.7)
C
C
      STOP
      END