      SUBROUTINE INITAL
C...
C...  RIVER POLLUTION MODEL
C...
C...  THE FOLLOWING MODEL HAS BEEN PROPOSED FOR POLLUTION TRANSPORT
C...  IN A RIVER (1):
C...
C...  U  = A*U   - B*U  - C*U + PS(X,T)                         (1)
C...   T      XX      X
C...
C...  WHERE
C...
C...  U       BIOCHEMICAL OXYGEN DEMAND (BOD) (M/L**3)
C...
C...  X       DISTANCE ALONG THE RIVER (L)
C...
C...  T       TIME (T)
C...
C...  A       DISPERSION COEFFICIENT (L**2/T)
C...
C...  B       VELOCITY OF FLOW (L/T)
C...
C...  C       BOD DECAY COEFFICIENT (1/T)
C...
C...  PS(X,T) POLLUTION SOURCES (M/L**3-T)
C...
C...  (1)  SHIBATA, Y., AND W. J. KARPLUS, DISTRIBUTED SOURCE IDENTIFI-
C...       CATION IN RIVER POLLUTION SYSTEMS BY PATTERN RECOGNITION,
C...       TRANS. SOCIETY FOR COMPUTER SIMULATION, VOL. 3, NO. 1, PP.
C...       1-29, JAN., 1986
C...
C...  EQAUTION (1) COULD, FOR EXAMPLE, BE USED TO SIMULATE THE POLLUTION
C...  OF THE RHINE RIVER ON NOVEMBER 1, 1986 RESULTING FROM A LARGE
C...  DISCHARGE OF PESTICIDE FROM A CHEMICAL PLANT IN SWITZERLAND.  THE
C...  SOURCE FUNCTION, SP(X,T), WOULD BE THE SPATIAL-TIME DESCRIPTION OF
C...  THE DISCHARGE.
C...
C...  THE FOLLOWING REFERENCES PERTAIN TO THE SPILL:
C...
C...  (2)  LEWIS, PAUL, 'HUGE CHEMICAL SPILL IN THE RHINE CREATES HAVOC
C...       IN FOUR COUNTRIES', THE NEW YORK TIMES, NOVEMBER 11, 1986
C...
C...  (3)  ROTH, TERRENCE, 'SWISS ACCIDENT THAT POLLUTED THE RHINE MAY
C...       HAVE HEAVIEST IMPACT ON GERMANY', THE WALL STREET JOURNAL,
C...       NOVEMBER 11, 1986
C...
C...  (4)  NETTER, THOMAS W., 'ANGER ALONG THE RHINE GROWS AFTER CHEMI-
C...       CAL SPILL', THE NEW YORK TIMES, NOVEMBER 12, 1986
C...
C...  (5)  STUDER, MARGARET A., 'SWISS NEVER THOUGHT SANDOZ DISASTER
C...       COULD HAPPEN IN THEIR ORDERLY WORLD', THE WALL STREET JOUR-
C...       NAL, NOVEMBER 13, 1986
C...
C...  (6)  'EUROPEAN OFFICIALS AGREE ON MEASURES TO PROTECT THE RHINE',
C...       THE WALL STREET JOURNAL, NOVEMBER 13, 1986
C...
C...  (7)  TAGLIABUE, JOHN, 'RHINE POISONING STRETCHING 185 MILES', THE
C...       NEW YORK TIMES, NOVEMBER 13, 1986
C...
C...  (8)  LEWIS, PAUL, 'SWISS WILL PAY COMPENSATION FOR SPILL', THE
C...       NEW YORK TIMES, NOVEMBER 13, 1986
C...
C...  (9)  MARCOM, JOHN, JR., 'EUROPE'S ENVIRONMENTAL POLICY FOUNDERS',
C...       THE WALL STREET JOURNAL, NOVEMBER 14, 1986
C...
C... (10)  'SWISS ANNOUNCE PROCEDURE FOR A CLEANUP OF THE RHINE', THE
C...       STREET JOURNAL, NOVEMBER 14, 1986
C...
C... (11)  NETTER, THOMAS W., 'MERCURY A KEY CONCERN IN RHINE SPILL',
C...       THE NEW YORK TIMES, NOVEMBER 15, 1986
C...
C...  THE INITIAL CONDITION FOR EQUATION (1) IS
C...
C...  U(X,0) = 0, X0 LE X LE XF                                 (2)
C...
C...  AND THE BOUNDARY CONDITIONS ARE
C...
C...  U(X0,T) = U(XF,T) = 0, T GT 0                          (3)(4)
C...
C...  WHERE
C...
C...  X0      X AT THE BEGINNING OF THE RIVER
C...
C...  XF      X AT THE END OF THE RIVER
C...
C...  COMPUTE A SOLUTION TO EQUATIONS (1) TO (4) BY THE NUMERICAL
C...  METHOD OF LINES ON A 51-POINT GRID (49 INTERIOR POINTS).  USE
C...  THE FOLLOWING NUMERICAL VALUES OF THE PARAMETERS
C...
C...  A      25 KM**2/HR
C...
C...  B      5 KM/HR
C...
C...  C      0 1/HR
C...
C...  X0     0 KM
C...
C...  XF     500 KM
C...
C...  PS(X,T)      0, X NE 250, T GE 0
C...
C...             100, X = 250, 0 LE T LE 10 HR
C...
C...               0, X = 250, T GT 10 HR
C...
C...  COMPUTE U(X,T) (U(X,T) AT THE 51 GRID POINTS) AND PLOT U(X,T)
C...  VS X FOR SIX VALUES OF T.
C...
C...  THE FOLLOWING CODING IS FOR THREE VARIATIONS OF THE MODEL.  IN
C...  PARTICULAR, THE EXIT BOUNDARY CONDITION (4) IS CHANGED, AS DES-
C...  CRIBED IN SUBROUTINE DERV.
C...
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(51)
     2      /F/ UT(51)
     3      /C/      A,      B,      C,     X0,     XF,     DX,    DXS,
     4           X(51),      N,     IP
C...
C...  MODEL PARAMETERS
      A=25.
      B=5.
      C=0.
      X0=0.
      XF=500.
      N=51
C...
C...  SPATIAL GRID
      DX=(XF-X0)/FLOAT(N-1)
      DXS=DX*DX
      DO 1 I=1,N
      X(I)=X0+DX*FLOAT(I-1)
1     CONTINUE
      DX=2.0*DX
C...
C...  MODEL INITIAL CONDITION (2)
      DO 2 I=1,N
      U(I)=0.
2     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
C...
C...  SUBROUTINE DERV CALLS SUBROUTINES DERV1 TO DERV5 FOR ALTERNATE
C...  NUMOL FORMULATIONS
C...
      COMMON/T/      T,  NSTOP,  NORUN
C...
C...  *****************************************************************
C...
C...  DERV1 IMPLEMENTS THE DIRICHLET BOUNDARY CONDITION U(XF,T) = 0.
C...  THIS CONDITION REPRESENTS A CONDITION OF NO POLLUTION AT THE END
C...  OF THE RIVER.
      IF(NORUN.EQ.1)CALL DERV1
C...
C...  *****************************************************************
C...
C...  DERV2 IMPLEMENTS THE NEUMANN BOUNDARY CONDITION U (XF,T) = 0.
C...                                                   X
C...  THIS CONDITION REPRESENTS PURE CONVECTION (FLOW) FROM THE END OF
C...  THE RIVER.
      IF(NORUN.EQ.2)CALL DERV2
C...
C...  *****************************************************************
C...
C...  DERV3 IMPLEMENTS THE TIME-DEPENDENT BOUNDARY CONDITION U (XF,T)
C...                                                          T
C...  = -B*U (XF,T).  THIS CONDITION REPRESENTS CONVECTION (FLOW) PLUS
C...        X
C...  ACCUMULATION AT THE END OF THE RIVER.
      IF(NORUN.EQ.3)CALL DERV3
C...
C...  *****************************************************************
C...
C...  DERV4 IMPLEMENTS THE NEUMANN BOUNDARY CONDITION U (XF,T) = 0.
C...                                                   X
C...  THIS CONDITION REPRESENTS PURE CONVECTION (FLOW) FROM THE END OF
C...  THE RIVER.  SPATIAL DERIVATIVES COMPUTED BY DSS020 AND DSS044.
      IF(NORUN.EQ.4)CALL DERV4
C...
C...  *****************************************************************
C...
C...  DERV5 IMPLEMENTS THE TIME-DEPENDENT BOUNDARY CONDITION U (XF,T)
C...                                                          T
C...  = -B*U (XF,T).  THIS CONDITION REPRESENTS CONVECTION (FLOW) PLUS
C...        X
C...  ACCUMULATION AT THE END OF THE RIVER. SPATIAL DERIVATIVES COM-
C...  PUTED BY DSS020 AND DSS044.
      IF(NORUN.EQ.5)CALL DERV5
C...
C...  *****************************************************************
C...
      RETURN
      END
      SUBROUTINE DERV1
C...
C...  SUBROUTINE DERV1 IMPLEMENTS THE NUMERICAL METHOD OF LINES
C...  SOLUTION OF EQUATIONS (1) TO (4)
C...
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(51)
     2      /F/ UT(51)
     3      /C/      A,      B,      C,     X0,     XF,     DX,    DXS,
     4           X(51),      N,     IP
C...
C...  BOUNDARY CONDITION (3)
      U(1)=0.
      UT(1)=0.
C...
C...  BOUNDARY CONDITION (4)
      U(N)=0.
      UT(N)=0.
C...
C...  EQUATION (1)
      DO 1 I=2,N-1
      UT(I)=A*(U(I+1)-2.0*U(I)+U(I-1))/DXS
     1     -B*(U(I+1)-U(I-1))/DX
     2     -C*U(I)+PS(X(I),T)
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV2
C...
C...  SUBROUTINE DERV2 IMPLEMENTS THE NUMERICAL METHOD OF LINES
C...  SOLUTION OF EQUATIONS (1) TO (3) WITH BOUNDARY CONDITION (4)
C...  REPLACED BY THE NEUMANN BOUNDARY CONDITION U (XF,T) = 0
C...                                              X
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(51)
     2      /F/ UT(51)
     3      /C/      A,      B,      C,     X0,     XF,     DX,    DXS,
     4           X(51),      N,     IP
C...
C...  BOUNDARY CONDITION (3)
      U(1)=0.
      UT(1)=0.
C...
C...  NEUMANN BOUNDARY CONDITION
      U(N)=(1.0/3.0)*(4.0*U(N-1)-U(N-2))
      UT(N)=0.
C...
C...  EQUATION (1)
      DO 1 I=2,N-1
      UT(I)=A*(U(I+1)-2.0*U(I)+U(I-1))/DXS
     1     -B*(U(I+1)-U(I-1))/DX
     2     -C*U(I)+PS(X(I),T)
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV3
C...
C...  SUBROUTINE DERV3 IMPLEMENTS THE NUMERICAL METHOD OF LINES
C...  SOLUTION OF EQUATIONS (1) TO (3) WITH BOUNDARY CONDITION (4)
C...  REPLACED BY THE TIME-DEPENDENT BOUNDARY CONDITION U (XF,T) =
C...  = -B*U (XF,T)                                      T
C...        X
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(51)
     2      /F/ UT(51)
     3      /C/      A,      B,      C,     X0,     XF,     DX,    DXS,
     4           X(51),      N,     IP
C...
C...  BOUNDARY CONDITION (3)
      U(1)=0.
      UT(1)=0.
C...
C...  TIME-DEPENDENT BOUNDARY CONDITION
      UT(N)=-B*(1.0/(2.0*DX))*(3.0*U(N)-4.0*U(N-1)+U(N-2))
C...
C...  EQUATION (1)
      DO 1 I=2,N-1
      UT(I)=A*(U(I+1)-2.0*U(I)+U(I-1))/DXS
     1     -B*(U(I+1)-U(I-1))/DX
     2     -C*U(I)+PS(X(I),T)
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV4
C...
C...  SUBROUTINE DERV4 IMPLEMENTS THE NUMERICAL METHOD OF LINES
C...  SOLUTION OF EQUATIONS (1) TO (3) WITH BOUNDARY CONDITION (4)
C...  REPLACED BY THE NEUMANN BOUNDARY CONDITION U (XF,T) = 0
C...                                              X
C...  AND THE SPATIAL DERIVATIVES COMPUTED BY DSS020 AND DSS044
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(51)
     2      /F/ UT(51)
     3      /C/      A,      B,      C,     X0,     XF,     DX,    DXS,
     4           X(51),      N,     IP
      REAL      UX(51),UXX(51)
C...
C...  BOUNDARY CONDITION (3)
      U(1)=0.
      UT(1)=0.
C...
C...  UX
      CALL DSS020(X0,XF,N,U,UX,1.0)
C...
C...  NEUMANN BOUNDARY CONDITION (INPUT TO DSS044)
      UX(N)=0.
C...
C...  UXX
      NL=1
      NU=2
      CALL DSS044(X0,XF,N,U,UX,UXX,NL,NU)
C...
C...  EQUATION (1)
      DO 1 I=2,N
      UT(I)=A*UXX(I)
     1     -B*UX(I)
     2     -C*U(I)+PS(X(I),T)
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV5
C...
C...  SUBROUTINE DERV5 IMPLEMENTS THE NUMERICAL METHOD OF LINES
C...  SOLUTION OF EQUATIONS (1) TO (3) WITH BOUNDARY CONDITION (4)
C...  REPLACED BY THE TIME-DEPENDENT BOUNDARY CONDITION U (XF,T) =
C...  = -B*U (XF,T)                                      T
C...        X
C...  AND THE SPATIAL DERIVATIVES COMPUTED BY DSS020 AND DSS044
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(51)
     2      /F/ UT(51)
     3      /C/      A,      B,      C,     X0,     XF,     DX,    DXS,
     4           X(51),      N,     IP
      REAL      UX(51),UXX(51)
C...
C...  BOUNDARY CONDITION (3)
      U(1)=0.
      UT(1)=0.
C...
C...  UX
      CALL DSS020(X0,XF,N,U,UX,1.0)
C...
C...  UXX
      NL=1
      NU=2
      CALL DSS044(X0,XF,N,U,UX,UXX,NL,NU)
C...
C...  EQUATION (1)
      DO 1 I=2,N-1
      UT(I)=A*UXX(I)
     1     -B*UX(I)
     2     -C*U(I)+PS(X(I),T)
1     CONTINUE
C...
C...  TIME-DEPENDENT BOUNDARY CONDITION IS THE NTH ODE
      UT(N)=-B*UX(N)
      RETURN
      END
      FUNCTION PS(X,T)
C...
C...  FUNCTION PS(X,T) DEFINES THE POLLUTION SOURCE FUNCTION
C...
C...  POLLUTION ENTERS THE RIVER AT X = 250
      IF((X.GT.249.5).AND.(X.LT.250.5))THEN
C...
C...  THE POLLUTION ENTERS THE RIVER FOR 10 TIME UNITS
         IF(T.LE.10.0)THEN
            PS=100.
         ELSE
C...
C...  AFTER TEN TIME UNITS, THE POLLUTION IS ENDED
            PS=0.
         END IF
      ELSE
C...
C...  THE RIVER IS NOT POLLUTED AT OTHER LOCATIONS (THAN X = 250)
         PS=0.
      END IF
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(51)
     2      /F/ UT(51)
     3      /C/      A,      B,      C,     X0,     XF,     DX,    DXS,
     4           X(51),      N,     IP
C...
C...  ADDITIONAL ARRAY FOR PLOTTING
      COMMON/CAL/ UP(6,51)
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IP=IP+1
      IF(IP.EQ.1)WRITE(NO,1)
1     FORMAT(9X,'T',4X,'U(0,T)',2X,'U(100,T)',2X,'U(200,T)',
     1            2X,'U(300,T)',2X,'U(400,T)',2X,'U(500,T)')
C...
C...  PRINT THE NUMERICAL SOLUTION
      WRITE(NO,2)T,U(1),U(11),U(21),U(31),U(41),U(51)
2     FORMAT(F10.2,6F10.4)
C...
C...  STORE THE SOLUTION FOR PLOTTING
      DO 3 I=1,N
      UP(IP,I)=U(I)
3     CONTINUE
C...
C...  PLOT THE SOLUTION
      IF(IP.LT.6)RETURN
      CALL TPLOTS(IP,N,X,UP)
      WRITE(NO,4)
4     FORMAT(1H ,/,' BIOCHEMICAL OXYGEN DEMAND (BOD) VS X',/,
     1             ' WITH T AS A PARAMETER, T = 0, 20,..., 100')
      RETURN
      END

