!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 6.3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! !
! Please Note: !
! !
! (1) This computer program is written by Tao Pang in conjunction with !
! his book, "An Introduction to Computational Physics," published !
! by Cambridge University Press in 1997. !
! !
! (2) No warranties, express or implied, are made for this program. !
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
PROGRAM G_WATER
!
! This program solves the groundwater dynamics problem in the
! rectangular geometry through the relaxation method.
! Copyright (c) Tao Pang 1997.
!
IMPLICIT NONE
INTEGER, PARAMETER :: NX=101,NY=51,ISKX=10,ISKY=5,ITMX=5
INTEGER :: I,J,ISTP
REAL :: PI,A0,B0,H0,CH,SX,SY,HX,HY,P,X,Y
REAL, DIMENSION (NX,NY) :: PHI,CK,SN
!
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
!
! Set up boundary conditions and initial guess of the solution
!
DO I = 1, NX
X = (I-1)*HX
DO J = 1, NY
Y = (J-1)*HY
SN(I,J) = 0.0
CK(I,J) = A0+B0*Y
PHI(I,J) = H0+CH*COS(PI*X/SX)*Y/SY
END DO
END DO
!
DO ISTP = 1, ITMX
!
! Ensure the boundary conditions by the 4-point formula
!
DO 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
END DO
!
CALL RX2D (PHI,CK,SN,NX,NY,P,HX,HY)
END DO
!
DO I = 1, NX, ISKX
X = (I-1)*HX
DO J = 1, NY, ISKY
Y = (J-1)*HY
WRITE (6,"(3F16.8)") X,Y,PHI(I,J)
END DO
END DO
END PROGRAM G_WATER
!
SUBROUTINE RX2D (FN,DN,S,NX,NY,P,HX,HY)
!
! Subroutine performing one iteration of the relaxation for
! the two-dimensional Poisson equation.
! Copyright (c) Tao Pang 1997.
!
IMPLICIT NONE
INTEGER, INTENT (IN) :: NX,NY
INTEGER :: I,J
REAL, INTENT (IN) :: HX,HY,P
REAL :: HX2,A,B,Q,CIP,CIM,CJP,CJM
REAL, INTENT (IN), DIMENSION (NX,NY) :: DN,S
REAL, INTENT (INOUT), DIMENSION (NX,NY) :: FN
!
HX2 = HX*HX
A = HX2/(HY*HY)
B = 1.0/(4.0*(1.0+A))
Q = 1.0-P
!
DO I = 2, NX-1
DO 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))
END DO
END DO
END SUBROUTINE RX2D