      SUBROUTINE INITAL
C...
C...  TWO DIMENSIONAL ADVECTION EQUATION
C...
C...     U  = -VX*U  - VY*U                               (1)
C...      T        X       Y
C...
C...     U(0,X,Y) = F1(-X/VX - Y/VY)                      (2)
C...
C...     U(T,0,Y) = F1(2*T - Y/VY)                        (3)
C...
C...     U(T,X,0) = F1(2*T - X/VX)                        (4)
C...
C...  ANALYTICAL SOLUTION
C...
C...     U(T,X,Y) = F1(2*T -X/VX - Y/VY)                  (5)
C...
C...  F1 IS AN ARBITRARILY DEFINED FUNCTION.  FOR THE FOLLOWING
C...  CODE, IT IS AN EXPONENTIAL FUNCTION CORRESPONDING TO A
C...  SMOOTH SOLUTION.  ALSO, IN THE FOLLOWING CODE, AN 11 X 11
C...  POINT GRID IS USED.  
C...  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (M=11, N=11)
      COMMON     /T/         T,      NFIN,      NRUN
     1           /Y/    U(M,N)
     2           /F/   UT(M,N)
     3           /S/   UX(M,N),   UY(M,N)
     4           /C/        VX,        VY,      X(M),      Y(N)
     5           /B/        XL,        XU,        YL,        YU
C...
C...  VELOCITIES
      VX=1.0D0
      VY=1.0D0
C...
C...  BOUNDARY VALUES OF X AND Y
      XL=0.D+00
      XU=1.D+00
      YL=0.D+00
      YU=1.D+00
C...
C...  GRID IN X
      DO 1 I=1,M
         X(I)=DFLOAT(I-1)/DFLOAT(M-1)*(XU-XL)
1     CONTINUE
C...
C...  GRID IN Y
      DO 2 J=1,N
         Y(J)=DFLOAT(J-1)/DFLOAT(N-1)*(YU-YL)
2     CONTINUE
C...
C...  INITIAL CONDITION
      DO 3 I=1,M
      DO 3 J=1,N
         ARG=-X(I)/VX-Y(J)/VY
         U(I,J)=F1(ARG)
3     CONTINUE
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (M=11, N=11)
      COMMON     /T/         T,      NFIN,      NRUN
     1           /Y/    U(M,N)
     2           /F/   UT(M,N)
     3           /S/   UX(M,N),   UY(M,N)
     4           /C/        VX,        VY,      X(M),      Y(N)
     5           /B/        XL,        XU,        YL,        YU
C...
C...  BOUNDARY CONDITION AT X = XL
      DO 1 J=1,N
         ARG=2.0D0*T-Y(J)/VY
         U(1,J)=F1(ARG)
1     CONTINUE
C...
C...  BOUNDARY CONDITION AT Y = YL
      DO 2 I=1,M
         ARG=2.0D0*T-X(I)/VX
         U(I,1)=F1(ARG)
2     CONTINUE
C...
C...  UX
      CALL DSS034(XL,XU,M,N,1,U,UX,1.0D+00)
C...
C...  UY
      CALL DSS034(YL,YU,M,N,2,U,UY,1.0D+00)
C...
C...  PDE
      DO 3 I=1,M
      DO 3 J=1,N
         UT(I,J)=-VX*UX(I,J)-VY*UY(I,J)
3     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (M=11, N=11)
      COMMON     /T/         T,      NFIN,      NRUN
     1           /Y/    U(M,N)
     2           /F/   UT(M,N)
     3           /S/   UX(M,N)
     4           /C/        VX,        VY,      X(M),      Y(N)
     5           /B/        XL,        XU,        YL,        YU
C...
C...  ARRAY FOR THE EXACT SOLUTION
      DIMENSION UE(M,N)
C...
C...  MONITOR SOLUTION ON SCREEN
      WRITE(*,*)T
C...
C...  COMPUTE THE EXACT SOLUTION
      DO 3 J=1,N,2
      DO 3 I=1,M,2
         ARG=2.0D0*T-X(I)/VX-Y(J)/VY
         UE(I,J)=F1(ARG)
3     CONTINUE
C...
C...  PRINT THE NUMERICAL AND EXACT SOLUTIONS
      WRITE(NO,1)T
1     FORMAT(//,5H T = ,F9.2,//,16X,
     1 10H    X=0   ,10H   X=0.2  ,10H   X=0.4  ,
     2 10H   X=0.6  ,10H   X=0.8  ,10H   X=1.0   )
      WRITE(NO,2)((U(I,J),I=1,M,2),(UE(I,J),I=1,M,2),J=1,N,2)
2     FORMAT(6H  Y=0 ,2X,5H  NUM,6F10.4,/,6X,2X,5HEXACT,6F10.4,//,
     1       6H Y=0.2,2X,5H  NUM,6F10.4,/,6X,2X,5HEXACT,6F10.4,//,
     1       6H Y=0.4,2X,5H  NUM,6F10.4,/,6X,2X,5HEXACT,6F10.4,//,
     1       6H Y=0.6,2X,5H  NUM,6F10.4,/,6X,2X,5HEXACT,6F10.4,//,
     1       6H Y=0.8,2X,5H  NUM,6F10.4,/,6X,2X,5HEXACT,6F10.4,//,
     1       6H Y=1.0,2X,5H  NUM,6F10.4,/,6X,2X,5HEXACT,6F10.4,//)
      RETURN
      END
      DOUBLE PRECISION FUNCTION F1(ARG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      F1=DEXP(ARG)
      RETURN
      END


