*DECK FLOW24
      SUBROUTINE INITAL
C...
C...  AQUIFER SIMULATION IN CARTESIAN COORDINATES
C...
C...  THE AQUIFER TO BE ANALYZED IS DEPICTED BELOW*
C...
C                   .
C                   .
C                   .           H (X,YL,T) = 0
C                   .            Y
C              Y=YL . . . . . . . . . . . . . . . . . S
C                   .   +                             .
C                   .   .                             .
C                   .   .                             .
C                   .   .                             .
C                   .   .                             .
C                   .   .                             .
C                   .   .                             .
C     H (0,Y,T) = 0 . 2000 M        H(X,Y,T)          . H (XL,Y,T) = 0
C      X            .   .                             .  X
C                 + .   .                             .
C                 . .   .                             .
C                 . .   .                             .
C                 . .   .                             .
C                 Y .   .                             .
C                   .+ . . . . . . 2000 M . . . . . .+.
C               Y=0 .   +                             .
C         . . . . . W . . . . . . . . . . . . . . . . . . . .
C                  X=0          H (X,0,T) = 0        X=XL
C                                Y
C                       X ---+
C
C...
C...  *THIS PROBLEM WAS PROVIDED BY DR. SHIV O. PRASHER, DEPARTMENT OF
C...  AGRICULTURAL ENGINEERING, MACDONALD COLLGE, MCGILL UNIVERSITY,
C...  21,111 LAKESHORE ROAD, STE. ANNE DE BELLEVUE, P.Q., H9X 1C0,
C...  CANADA
C...
C...  THE AQUIFER IS MODELED AS A RECTANGULAR DOMAIN IN THE X-Y
C...  COORDINATE SYSTEM, MEASURING XL (2000 M) BY YL (2000 M). THE
C...  HEIGHT OF WATER AT AN POINT IN THE AQUIFER IS H(X,Y,T).  ZERO
C...  FLOW BOUNDARY CONDITIONS ARE ASSUMED AT X = 0, X = XL, Y = 0
C...  AND Y = YL.
C...
C...  TWO CASES ARE CONSIDERED
C...
C...     (1)  A SINK (WELL FOR WITHDRAWAL OF WATER) IS LOCATED AT X
C...          = 0, Y = 0 (DESIGNATED AS W)
C...
C...     (2)  A SOURCE IS ADDED TO THE SINK OF (1) AT X = XL, Y = YL
C...          (DESIGNATED AS S).  THE SOURCE AND SINK OF (1) AND (2)
C...          ARE OF EQUAL STRENGTH (FLOW RATE) SO THAT THE ENTIRE
C...          SYSTEM COMES TO A STEADY STATE.
C...
C...  THE FOLLOWING TWO-DIMENSIONAL PARTIAL DIFFERENTIAL EQUATION
C...  (PDE) RELATES THE WATER HEIGHT TO X, Y AND T (TIME)
C...
C...     SC*H  = TR*(H   + H  ) + RR                               (1)
C...         T        XX    YY
C...
C...  WHERE
C...
C...     H      HYDRAULIC HEAD
C...
C...     X      HORIZONTAL COORDINATE
C...
C...     Y      VERTICAL COORDINATE
C...
C...     T      TIME
C...
C...     TR     TRANSMISSIVITY
C...
C...     RR     RECHARGE RATE
C...
C...     SC     STORAGE COEFFICIENT
C...
C...  THE INITIAL AND BOUNDARY CONDITIONS FOR EQUATION (1) ARE
C...
C...     H (0,Y,T) = H (XL,Y,T) = H (X,0,T) = H (X,YL,T) = 0
C...      X           X            Y           Y
C...                                                   (2)(3)(4)(5)(6)
C...     H(X,Y,0) = H0
C...
C...     H0     INITIAL HYDRAULIC HEAD
C...
C...  THE FOLLOWING CODING IS FOR EQUATIONS (1) TO (6).
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=11)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    H(N,N)
     2      /F/   HT(N,N)
     3      /C/        QS,        TR,        SC,
     4                 XL,        YL,         A,
     5                 DX,        DY,        RR,
     6                 H0
     7      /I/        IP
C...
C...  MODEL PARAMETERS
C...
C...     SOURCE AND SINK FLOW RATES (1/DAY)
         QS=2000.0D0
C...
C...     TRANSMISSIVITY (M**3/DAY)
         TR=300.0D0
C...
C...     STORAGE COEFFICIENT (DIMENSIONLESS)
         SC=0.02D0
C...
C...     LENGTH OF THE AQUIFER IN THE X DIRECTION (M)
         XL=2000.0D0
C...
C...     LENGTH OF THE AQUIFER IN THE Y DIRECTION (M)
         YL=2000.0D0
C...
C...     GRID SPACING IN THE X DIRECTION
         DX=XL/DFLOAT(N-1)
C...
C...     GRID SPACING IN THE Y DIRECTION
         DY=YL/DFLOAT(N-1)
C...
C...     AREA OF A GRID UNIT (M**2)
         A=DX*DY
C...
C...     RECHARGE RATE (M/DAY)
         RR=QS/A
C...
C...     INITIAL HYDRAULIC HEAD (M)
         H0=10.0D0
C...
C...  INITIAL CONDITIONS
      DO 1 I=1,N
      DO 1 J=1,N
         H(I,J)=H0
1     CONTINUE
C...
C...  COMPUTE THE INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=11)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    H(N,N)
     2      /F/   HT(N,N)
     3      /C/        QS,        TR,        SC,
     4                 XL,        YL,         A,
     5                 DX,        DY,        RR,
     6                 H0
     7      /I/        IP
C...
C...  SET CONSTRAINTS ON THE SOLUTION
      DO 7 I=1,N
      DO 8 J=1,N
         IF(H(I,J).LT.0.0D0)H(I,J)=0.0D0
8     CONTINUE
7     CONTINUE
C...
C...  BOUNDARY POINTS FROM I = 2 TO N-1
      DO 1 I=2,N-1
C...
C...  Y = 0
      HT(I,1)=
     1  (TR/SC)*(      H(I+1,  1)-2.0D0*H(I  ,  1)+   H(I-1,  1))/DX**2
     2 +(TR/SC)*(2.0D0*H(I  ,  2)-2.0D0*H(I  ,  1)              )/DY**2
C...
C...  Y = YL
      HT(I,N)=
     1  (TR/SC)*(      H(I+1,N  )-2.0D0*H(I  ,N  )+   H(I-1,N  ))/DX**2
     2 +(TR/SC)*(2.0D0*H(I  ,N-1)-2.0D0*H(I  ,N  )              )/DY**2
1     CONTINUE
C...
C...  BOUNDARY POINTS FROM J = 2 TO N-1
      DO 2 J=2,N-1
C...
C...  X = 0
      HT(1,J)=
     1  (TR/SC)*(2.0D0*H(  2,J  )-2.0D0*H(  1,J  )              )/DX**2
     2 +(TR/SC)*(      H(  1,J+1)-2.0D0*H(  1,J  )+   H(  1,J-1))/DY**2
C...
C...  X = XL
      HT(N,J)=
     1  (TR/SC)*(2.0D0*H(N-1,J  )-2.0D0*H(N  ,J  )              )/DX**2
     2 +(TR/SC)*(      H(N  ,J+1)-2.0D0*H(N  ,J  )+   H(N  ,J-1))/DY**2
2     CONTINUE
C...
C...  CORNER POINTS
C...
C...  X = 0, Y = 0
      HT(1,1)=
     1  (TR/SC)*(2.0D0*H(  2,  1)-2.0D0*H(  1,  1)              )/DX**2
     2 +(TR/SC)*(2.0D0*H(  1,  2)-2.0D0*H(  1,  1)              )/DY**2
C...
C...  X = XL, Y = 0
      HT(N,1)=
     1  (TR/SC)*(2.0D0*H(N-1,  1)-2.0D0*H(N  ,  1)              )/DX**2
     2 +(TR/SC)*(2.0D0*H(N  ,  2)-2.0D0*H(N  ,  1)              )/DY**2
C...
C...  X = 0, Y = YL
      HT(1,N)=
     1  (TR/SC)*(2.0D0*H(  2,N  )-2.0D0*H(  1,N  )              )/DX**2
     2 +(TR/SC)*(2.0D0*H(  1,N-1)-2.0D0*H(  1,N  )              )/DY**2
C...
C...  X = XL, Y = YL
      HT(N,N)=
     1  (TR/SC)*(2.0D0*H(N-1,N  )-2.0D0*H(N  ,N  )              )/DX**2
     2 +(TR/SC)*(2.0D0*H(N  ,N-1)-2.0D0*H(N  ,N  )              )/DY**2
C...
C...  INTERIOR POINTS
      DO 4 I=2,N-1
      DO 4 J=2,N-1
      HT(I,J)=
     1  (TR/SC)*(      H(I+1,J  )-2.0D0*H(I  ,J  )+   H(I-1,J  ))/DX**2
     2 +(TR/SC)*(      H(I  ,J+1)-2.0D0*H(I  ,J  )+   H(I  ,J-1))/DY**2
4     CONTINUE
C...
C...  INCLUDE THE SOURCE AND SINK
      IF(NORUN.EQ.1)THEN
         HT(1,1)=HT(1,1)-RR/SC
      ELSE IF(NORUN.EQ.2)THEN
         HT(1,1)=HT(1,1)-RR/SC
         HT(N,N)=HT(N,N)+RR/SC
      END IF
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=11)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    H(N,N)
     2      /F/   HT(N,N)
     3      /C/        QS,        TR,        SC,
     4                 XL,        YL,         A,
     5                 DX,        DY,        RR,
     6                 H0
     7      /I/        IP
C...
C...  ARRAYS FOR PLOTTING
      DIMENSION TP(11), HP(3,11), HTP(3,11)
C...
C...  MONITOR THE CALCULATION
      WRITE(*,*)NORUN,T
C...
C...  PRINT THE SOLUTION
      NS=2
      WRITE(NO,1)T
1     FORMAT(1H ,//,5H T = ,F6.1,//,11H   H(X,Y,T))
      WRITE(NO,2)((  H(I,J),I=1,N,NS),J=1,N,NS)
      WRITE(NO,3)
3     FORMAT(11H  HT(X,Y,T))
      WRITE(NO,2)(( HT(I,J),I=1,N,NS),J=1,N,NS)
2     FORMAT(
     +           9X,10H     X = 0,10H   X = 400,10H   X = 800,
     +              10H  X = 1200,10H  X = 1600,10H  X = 2000,/,
     +               9H Y = 0   ,6F10.4,/,9H Y = 400 ,6F10.4,/,
     +               9H Y = 800 ,6F10.4,/,9H Y = 1200,6F10.4,/,
     +               9H Y = 1600,6F10.4,/,9H Y = 2000,6F10.4,//)
C...
C...  STORE THE SOLUTION AND THE TEMPORAL DERIVATIVE FOR PLOTTING
      IP=IP+1
      TP(IP)=T
      HP(1,IP)=H(1,1)
      HP(2,IP)=H(6,6)
      HP(3,IP)=H(N,N)
      HTP(1,IP)=HT(1,1)
      HTP(2,IP)=HT(6,6)
      HTP(3,IP)=HT(N,N)
C...
C...  PRINT THE SOLUTION AND THE TEMPORAL DERIVATIVE AT THE END OF THE
C...  RUN FOR THE POINTS (X=0,Y=0), (X=XL/2,Y=YL/2), (X=XL,Y=YL)
      IF(IP.EQ.11)THEN
         DO 5 J=1,IP
         WRITE(NO,4)TP(J),(HP(I,J),I=1,3),(HTP(I,J),I=1,3)
4        FORMAT(1X,'T = ',F6.1,/,3F12.4,/,3F12.4,/)
5        CONTINUE
C...
C...  PLOT THE SOLUTION AND THE TEMPORAL DERIVATIVE AT THE END OF THE
C...  RUN FOR THE POINTS (X=0,Y=0), (X=XL/2,Y=YL/2), (X=XL,Y=YL)
         CALL SPLOTS(3,IP,TP,HP)
         WRITE(NO,6),NORUN
6        FORMAT(/,' RUN ',I2,3X,
     1   ' 1 - H(0,0), 2 - H(XL/2,YL/2), 3 - H(XL,YL)',/)
         CALL SPLOTS(3,IP,TP,HTP)
         WRITE(NO,7),NORUN
7        FORMAT(/,' RUN ',I2,3X,
     1   ' 1 - HT(0,0), 2 - HT(XL/2,YL/2), 3 - HT(XL,YL)',/)
      END IF
      RETURN
      END
