*DECK PRO4P10
      SUBROUTINE INITAL
C...
C...  THE TEST PROBLEM FOR NEUMANN BOUNDARY CONDITIONS IS
C...
C...     U  = D*U                                                 (1)
C...      T      XX
C...
C...     U(X,0) = 0, 0 LE X LT 0.5
C...                                                              (2)
C...     U(X,0) = 1, 0.5 LT X LE 1
C...
C...     U (0,T) = 0, U (L,T) = 0                              (3)(4)
C...      X            X
C...
C...  WITH THE STEADY STATE SOLUTION U(X,INF) = 0.5
C...
C...  THE FOLLOWING CODE COMPUTES A NUMERICAL SOLUTION TO EQUATIONS
C...  (1) TO (4) BY
C...
C...     (1)  FINITE DIFFERENCES
C...
C...          NORUN = 1, THREE POINT DIFFERENCES
C...
C...          NORUN = 2, FIVE POINT DIFFERENCES
C...
C...     (2)  FINITE ELEMENTS
C...
C...          NORUN = 3, LINEAR FINITE ELEMENTS
C...
C...     (3)  FINITE VOLUMES
C...
C...          NORUN = 4, GRID WITH EQUAL VOLUMES
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/      T,  NSTOP,  NORUN
C...
C...  SOLUTION OF EQUATIONS (1) TO (4) BY FINITE DIFFERENCES
      IF((NORUN.EQ.1).OR.(NORUN.EQ.2))CALL INIT1
C...
C...  SOLUTION OF EQUATIONS (1) TO (4) BY LINEAR FINITE ELEMENTS
      IF(NORUN.EQ.3)CALL INIT2
C...
C...  SOLUTION OF EQUATIONS (1) TO (4) BY FINITE VOLUMES
      IF(NORUN.EQ.4)CALL INIT3
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      COMMON/T/      T,  NSTOP,  NORUN
C...
C...  SOLUTION OF EQUATIONS (1) TO (4) BY FINITE DIFFERENCES
      IF((NORUN.EQ.1).OR.(NORUN.EQ.2))CALL DERV1
C...
C...  SOLUTION OF EQUATIONS (1) TO (4) BY LINEAR FINITE ELEMENTS
      IF(NORUN.EQ.3)CALL DERV2
C...
C...  SOLUTION OF EQUATIONS (1) TO (4) BY FINITE VOLUMES
      IF(NORUN.EQ.4)CALL DERV3
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/      T,  NSTOP,  NORUN
C...
C...  SOLUTION OF EQUATIONS (1) TO (4) BY FINITE DIFFERENCES
      IF((NORUN.EQ.1).OR.(NORUN.EQ.2))CALL PRINT1(NI,NO)
C...
C...  SOLUTION OF EQUATIONS (1) TO (4) BY LINEAR FINITE ELEMENTS
      IF(NORUN.EQ.3)CALL PRINT2(NI,NO)
C...
C...  SOLUTION OF EQUATIONS (1) TO (4) BY FINITE VOLUMES
      IF(NORUN.EQ.4)CALL PRINT3(NI,NO)
      RETURN
      END
      SUBROUTINE INIT1
C...
C...  SOLUTION OF EQUATIONS (1) TO (4) BY FINITE DIFFERENCES
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=51)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  LENGTH
      L=1.0D0
C...
C...  GRID SPACING
      DX=L/DFLOAT(N-1)
C...
C...  INITIAL CONDITION (2)
      DO 1 I=1,N
         X(I)=DFLOAT(I-1)*DX
         IF(I.LT.(N+1)/2)U(I)=0.0D0
         IF(I.GT.(N+1)/2)U(I)=1.0D0
         IF(I.EQ.(N+1)/2)U(I)=0.5D0
1     CONTINUE
C...
C...  INITIALIZE THE TEMPORAL DERIVATIVES IN COMMON/F/
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE INIT2
C...
C...  SOLUTION OF EQUATIONS (1) TO (4) BY LINEAR FINITE ELEMENTS (HAT
C...  FUNCTIONS)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=51)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
     5     /FE/  AL(N),  BM(N),  CU(N), BRHS(N)
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  LENGTH
      L=1.0D0
C...
C...  GRID SPACING
      DX=L/DFLOAT(N-1)
C...
C...  SET UP THE COEFFICIENT MATRIX IN BAND STORAGE MODE. THIS IS
C...  DONE ONLY ONCE SINCE THE COEFFICIENT MATRIX IS CONSTANT
C...
C...  LOWER DIAGONAL
      DO 2 I=1,N
         IF(I.EQ.1)THEN
            AL(1)=0.0D0
         ELSE
     +   IF(I.EQ.N)THEN
            AL(N)=1.0D0/6.0D0
         ELSE
            AL(I)=1.0D0/6.0D0
         END IF
2     CONTINUE
C...
C...  MAIN DIAGONAL
      DO 3 I=1,N
         IF(I.EQ.1)THEN
            BM(1)=2.0D0/6.0D0
         ELSE
     +   IF(I.EQ.N)THEN
            BM(N)=2.0D0/6.0D0
         ELSE
            BM(I)=4.0D0/6.0D0
         END IF
3     CONTINUE
C...
C...  UPPER DIAGONAL
      DO 4 I=1,N
         IF(I.EQ.1)THEN
            CU(1)=1.0D0/6.0D0
         ELSE
     +   IF(I.EQ.N)THEN
            CU(N)=0.0D0
         ELSE
            CU(I)=1.0D0/6.0D0
         END IF
4     CONTINUE
C...
C...  INITIAL CONDITION (2)
      DO 1 I=1,N
         X(I)=DFLOAT(I-1)*DX
         IF(I.LT.(N+1)/2)U(I)=0.0D0
         IF(I.GT.(N+1)/2)U(I)=1.0D0
         IF(I.EQ.(N+1)/2)U(I)=0.5D0
1     CONTINUE
C...
C...  INITIALIZE THE TEMPORAL DERIVATIVES IN COMMON/F/
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE INIT3
C...
C...  SOLUTION OF EQUATIONS (1) TO (4) BY FINITE VOLUMES
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=50)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  LENGTH
      L=1.0D0
C...
C...  GRID SPACING
      DX=L/DFLOAT(N)
C...
C...  INITIAL CONDITION (2)
      DO 1 I=1,N
         X(I)=DFLOAT(I-1)*DX+0.5D0*DX
         IF(I.LE.(N/2))U(I)=0.0D0
         IF(I.GT.(N/2))U(I)=1.0D0
1     CONTINUE
C...
C...  INITIALIZE THE TEMPORAL DERIVATIVES IN COMMON/F/
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV1
C...
C...  SUBROUTINE DERV1 IS CALLED BY THE ODE INTEGRATOR TO DEFINE THE
C...  ODES IN THE FINITE DIFFERENCE SOLUTION OF EQUATIONS (1) TO (4)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=51)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  DEFINE ARRAYS FOR THE SPATIAL DERIVATIVES
      DIMENSION UX(N), UXX(N)
C...
C...  BOUNDARY CONDITION (3)
      UX(1)=0.0D0
      NL=2
C...
C...  BOUNDARY CONDITION (4)
      UX(N)=0.0D0
      NU=2
C...
C...  U  BY THREE POINT DIFFERENCES
C...   XX
      IF(NORUN.EQ.1)CALL DSS042(0.,L,N,U,UX,UXX,NL,NU)
C...
C...  U  BY FIVE POINT DIFFERENCES
C...   XX
      IF(NORUN.EQ.2)CALL DSS044(0.,L,N,U,UX,UXX,NL,NU)
C...
C...  EQUATION (1)
      DO 1 I=1,N
         UT(I)=UXX(I)
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV2
C...
C...  SUBROUTINE DERV2 IS CALLED BY THE ODE INTEGRATOR TO DEFINE THE
C...  ODES IN THE FINITE ELEMENT SOLUTION OF EQUATIONS (1) TO (4)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=51)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
     5     /FE/  AL(N),  BM(N),  CU(N), BRHS(N)
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  BOUNDARY CONDITION (3)
      H1=0.0D0
C...
C...  BOUNDARY CONDITION (4)
      HN=0.0D0
C...
C...  RIGHT HAND SIDE VECTOR (WITH D = 1)
      DO 1 I=1,N
         IF(I.EQ.1)THEN
            BRHS(I)=-H1/DX+(U(I+1)-U(I))/DX**2
         ELSE
     +   IF(I.EQ.N)THEN
            BRHS(I)= HN/DX+(U(I-1)-U(I))/DX**2
         ELSE
            BRHS(I)=(U(I+1)-2.0D0*U(I)+U(I-1))/DX**2
         END IF
1     CONTINUE
C...
C...  SOLVE THE LINEAR ALGEBRAIC EQUATIONS BY SUBROUTINE TRIDAG, WHICH
C...  RETURNS THE DERIVATIVE VECTOR UT IN COMMON/F/
      CALL TRIDAG(AL,BM,CU,BRHS,UT,N)
      RETURN
      END
      SUBROUTINE DERV3
C...
C...  SUBROUTINE DERV3 IS CALLED BY THE ODE INTEGRATOR TO DEFINE THE
C...  ODES IN THE FINITE VOLUME SOLUTION OF EQUATIONS (1) TO (4)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=50)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  DIMENSION ARRAY FOR THE FLUXES
      DIMENSION Q(0:N)
C...
C...  BOUNDARY CONDITION (3)
      Q(0)=0.0D0
C...
C...  BOUNDARY CONDITION (4)
      Q(N)=0.0D0
C...
C...  COMPUTE THE FLUXES (K/(RHO*CP) = 1)
      DO 2 I=1,N-1
         Q(I)=(U(I)-U(I+1))/DX
2     CONTINUE
C...
C...  EQUATION (1)
      DO 1 I=1,N
         UT(I)=(Q(I-1)-Q(I))/DX
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT1(NI,NO)
C...
C...  SUBROUTINE PRINT1 IS CALLED BY THE MAIN PROGRAM TO PRINT THE
C...  FINITE DIFFERENCE SOLUTION OF EQUATIONS (1) TO (4)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=51)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  DIMENSION ARRAYS FOR THE ANALYTICAL SOLUTION
      DIMENSION UE(N), DIFF(N)
C...
C...  COMPUTE THE TOTAL ENERGY (RHO = CP = 1)
      CALL DERV1
      ENERGY=0.0D0
      DO 1 I=2,N
         ENERGY=ENERGY+0.5D0*(U(I)+U(I-1))*DX
1     CONTINUE
C...
C...  CALCULATE THE EXACT SOLUTION, AND THE DIFFERENCE BETWEEN THE
C...  NUMERICAL AND EXACT SOLUTIONS
      DO 3 I=1,N,10
         UE(I)=SERIES(I)
         DIFF(I)=U(I)-UE(I)
3     CONTINUE
C...
C...  PRINT THE NUMERICAL AND EXACT SOLUTIONS TO EQUATIONS (1) TO (4)
      WRITE(NO,2)T,ENERGY,(U(I),I=1,N,10),
     1                   (UE(I),I=1,N,10),
     2                 (DIFF(I),I=1,N,10)
2     FORMAT('    T = ',F6.2,'  ENERGY = ',F9.5,
     1           /,14X,'   X=0',5X,'X=0.2',5X,'X=0.4',
     2               5X,'X=0.6',5X,'X=0.8',5X,'X=1.0D0'/,
     3                         '    U(X,T)',6F10.6,/,
     4                         '   UE(X,T)',6F10.6,/,
     5                         ' DIFF(X,T)',6F10.6,/)
      RETURN
      END
      SUBROUTINE PRINT2(NI,NO)
C...
C...  SUBROUTINE PRINT2 IS CALLED BY THE MAIN PROGRAM TO PRINT THE
C...  FINITE ELEMENT SOLUTION OF EQUATIONS (1) TO (4)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=51)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
     5     /FE/  AL(N),  BM(N),  CU(N), BRHS(N)
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  DIMENSION ARRAYS FOR THE ANALYTICAL SOLUTION
      DIMENSION UE(N), DIFF(N)
C...
C...  PRINT THE COEFFICIENT MATRIX FOR VERIFICATION
      IP=IP+1
      IF(IP.EQ.1)THEN
         WRITE(NO,13)
13       FORMAT(/,' COEFFICIENT MATRIX',/)
         DO 14 I=1,N
            WRITE(NO,15)I,AL(I),BM(I),CU(I)
15          FORMAT(I5,3F12.4)
14       CONTINUE
         WRITE(NO,16)
16       FORMAT(//)
      END IF
C...
C...  COMPUTE THE TOTAL ENERGY (RHO = CP = 1)
      CALL DERV2
      ENERGY=0.0D0
      DO 1 I=2,N
         ENERGY=ENERGY+0.5D0*(U(I)+U(I-1))*DX
1     CONTINUE
C...
C...  CALCULATE THE EXACT SOLUTION, AND THE DIFFERENCE BETWEEN THE
C...  NUMERICAL AND EXACT SOLUTIONS
      DO 3 I=1,N,10
         UE(I)=SERIES(I)
         DIFF(I)=U(I)-UE(I)
3     CONTINUE
C...
C...  PRINT THE NUMERICAL AND EXACT SOLUTIONS TO EQUATIONS (1) TO (4)
      WRITE(NO,2)T,ENERGY,(U(I),I=1,N,10),
     1                   (UE(I),I=1,N,10),
     2                 (DIFF(I),I=1,N,10)
2     FORMAT('    T = ',F6.2,'  ENERGY = ',F9.5,
     1           /,14X,'   X=0',5X,'X=0.2',5X,'X=0.4',
     2               5X,'X=0.6',5X,'X=0.8',5X,'X=1.0D0'/,
     3                         '    U(X,T)',6F10.6,/,
     4                         '   UE(X,T)',6F10.6,/,
     5                         ' DIFF(X,T)',6F10.6,/)
      RETURN
      END
      SUBROUTINE PRINT3(NI,NO)
C...
C...  SUBROUTINE PRINT3 IS CALLED BY THE MAIN PROGRAM TO PRINT THE
C...  FINITE VOLUME SOLUTION OF EQUATIONS (1) TO (4)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=50)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  DIMENSION ARRAYS FOR THE ANALYTICAL SOLUTION
      DIMENSION UE(N), DIFF(N)
C...
C...  COMPUTE THE TOTAL ENERGY (RHO = CP = 1)
      CALL DERV3
      ENERGY=0.0D0
      DO 1 I=1,N
         ENERGY=ENERGY+U(I)*DX
1     CONTINUE
C...
C...  CALCULATE THE EXACT SOLUTION, AND THE DIFFERENCE BETWEEN THE
C...  NUMERICAL AND EXACT SOLUTIONS
      DO 3 I=1,21,5
         UE(I)=SERIES(I)
         DIFF(I)=U(I)-UE(I)
3     CONTINUE
C...
C...  PRINT THE NUMERICAL AND EXACT SOLUTIONS TO EQUATIONS (1) TO (4)
      UAVG=(U(N/2)+U(N/2+1))/2.0D0
      WRITE(NO,2)T,ENERGY,(U(I),I=1,21,5),UAVG,
     1                   (UE(I),I=1,21,5),
     2                 (DIFF(I),I=1,21,5)
2     FORMAT('    T = ',F6.2,'  ENERGY = ',F9.5,
     1           /,14X,'X=0.01',4X,'X=0.11',4X,'X=0.21',
     2              4X,'X=0.31',4X,'X=0.41',4X,'X=0.50'/,
     3                         '    U(X,T)',6F10.6,/,
     4                         '   UE(X,T)',5F10.6,/,
     5                         ' DIFF(X,T)',5F10.6,/)
      RETURN
      END
      DOUBLE PRECISION FUNCTION SERIES(IG)
C...
C...  FUNCTION SERIES EVALUATES THE ANALYTICAL SOLUTION
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=51)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  DEFINE MAXIMUM NUMBER OF TERMS IN THE SERIES (THIS IS SET TO A
C...  HIGH VALUE TO ACHIEVE CONVERGENCE TO THREE FIGURES FOR T = 0;
C...  BEYOND T = 0, THE SERIES CONVERGES RAPIDLY)
      PARAMETER (IS=1000)
C...
C...  PI
      PI=4.0D0*DATAN(1.0D0)
C...
C...  EVALUATE SERIES SOLUTION
      SUM=0.25D0
      SIGN=1.0D0
      DO 1 I=1,IS,2
         EI=DFLOAT(I)*PI
         EXI=-(EI**2)*T
         SIGN=-1.0D0*SIGN
         CI=SIGN/EI
         SUM=SUM+CI*DEXP(EXI)*DCOS(EI*X(IG))
C...
C...     AVOID UNDERFLOW OF THE EXP FUNCTION
         IF(EXI.LT.-100.0D0)GO TO 2
1     CONTINUE
C...
C...  SUMMATION OF SERIES IS COMPLETE
2     SERIES=2.0*SUM
      END
