ccccccccccccccccccccccccc     Program 6.3     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      PROGRAM G_WATER
C
C This program solves the groundwater dynamics problem in the
C rectangular geometry through the relaxation method.
C
      PARAMETER (NX=1001,NY=501,ITMX=5)
      DIMENSION PHI(NX,NY),CK(NX,NY),SN(NX,NY)
C
      PI = 4.0*ATAN(1.0)
      A0 = 1.0
      B0 = -0.04
      H0 = 200.0
      CH = -20.0
      SX = 1000.0
      SY = 500.0
      HX = SX/(NX-1)
      HY = SY/(NY-1)
      P  = 0.5
C
C Set up boundary conditions and initial guess of the solution
C
      DO    150  I = 1, NX
        X = (I-1)*HX
        DO  100  J = 1, NY
          Y = (J-1)*HY
          CK(I,J)  = A0+B0*Y
          PHI(I,J) = H0+CH*COS(PI*X/SX)*Y/SY 
          SN(I,J)  = 0.0
  100   CONTINUE
  150 CONTINUE
C
      DO 200 ISTP = 1, ITMX
C
C Ensure the boundary conditions by the 4-point formula
C
        DO  180 J = 1, NY
          PHI(1,J)  =(4.0*PHI(2,J)-PHI(3,J))/3.0
          PHI(NX,J) =(4.0*PHI(NX-1,J)-PHI(NX-2,J))/3.0
  180   CONTINUE
C
        CALL RX2D (PHI,CK,SN,NX,NY,P,HX,HY)
  200 CONTINUE
C
      DO   400  I = 1, NX, 100
        X = (I-1)*HX
        DO 390  J = 1, NY, 50
          Y = (J-1)*HY
          WRITE (6,9999) X, Y, PHI(I,J)
  390   CONTINUE
  400 CONTINUE
 9999 FORMAT (3F14.6)
      STOP
      END
C
      SUBROUTINE RX2D (FN,DN,S,NX,NY,P,HX,HY)
C
C Subroutine performing one iteration of the relaxation for
C the two-dimensional Poisson equation.
C
      DIMENSION FN(NX,NY),DN(NX,NY),S(NX,NY)
      HX2 = HX*HX
      A   = HX2/(HY*HY)
      B   = 1.0/(4.0*(1.0+A))
      Q   = 1.0-P
      DO    200 I = 2, NX-1
        DO  100 J = 2, NY-1
          CIP     = B*(DN(I+1,J)/DN(I,J)+1.0)
          CIM     = B*(DN(I-1,J)/DN(I,J)+1.0)
          CJP     = A*B*(DN(I,J+1)/DN(I,J)+1.0)
          CJM     = A*B*(DN(I,J-1)/DN(I,J)+1.0)
          FN(I,J) = Q*FN(I,J)+P*(CIP*FN(I+1,J)+CIM*FN(I-1,J)
     *             +CJP*FN(I,J+1)+CJM*FN(I,J-1)+HX2*S(I,J))
  100   CONTINUE
  200 CONTINUE
      RETURN
      END