      SUBROUTINE INITAL
C...
C...  M3 - NONLINEAR, FRONT-SHARPENING, CONVECTIVE FLOW
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  SELECT CASE
      NCASE=3
C...
C...  SELECT INITIALIZATION ROUTINE
      IF(NCASE.EQ.1)CALL INIT1
      IF(NCASE.EQ.2)CALL INIT2
      IF(NCASE.EQ.3)CALL INIT3
      RETURN
      END
      SUBROUTINE INIT1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=101)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NX)
     +      /F/ UT(NX)
     +      /S/ UX(NX)
     +      /X/  X(NX)
     +      /C/      A,      B,      C,     XL,     DX
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  SELECT THE SOLUTION PARAMETERS
C...
C...     FRONT FLATTENING
         IF(NORUN.EQ.1)THEN
         A=1.0D0
         B=1.0D0
         C=0.0D0
         NT=11
C...
C...     FRONT SHARPENING
         ELSE
     +   IF(NORUN.EQ.2)THEN
         A= 1.01D0
         B=-1.0D0
         C= 1.0D0
         NT=11
         END IF
C...
C...  NUMBER OF GRID POINTS
      N=11
C...
C...  TOTAL LENGTH IN X, GRID SPACING
      XL=1.0D0
      DX=XL/DFLOAT(N-1)
C...
C...  INITIAL CONDITION
      DO 1 I=1,N
         X(I)=XL*DFLOAT(I-1)/DFLOAT(N-1)
         U(I)=ANAL(X(I),0.0D0)
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE INIT2
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=101)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NX)
     +      /F/ UT(NX)
     +      /S/ UX(NX)
     +      /X/  X(NX)
     +      /C/      A,      B,      C,     XL,     DX
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  SELECT THE SOLUTION PARAMETERS
C...
C...     FRONT FLATTENING
         IF(NORUN.EQ.1)THEN
         A=1.0D0
         B=1.0D0
         C=0.0D0
         NT=11
C...
C...     FRONT SHARPENING
         ELSE
     +   IF(NORUN.EQ.2)THEN
         A= 1.0D0
         B=-1.0D0
         C= 1.0D0
         NT=13
         END IF
C...
C...  NUMBER OF GRID POINTS
      N=101
C...
C...  TOTAL LENGTH IN X, GRID SPACING
      XL=1.0D0
      DX=XL/DFLOAT(N-1)
C...
C...  INITIAL CONDITION
      DO 1 I=1,N
         X(I)=XL*DFLOAT(I-1)/DFLOAT(N-1)
         U(I)=ANAL(X(I),0.0D0)
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE INIT3
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=101)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NX)
     +      /F/ UT(NX)
     +      /S/ UX(NX)
     +      /X/  X(NX)
     +      /C/      A,      B,      C,     XL,     DX
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  SELECT THE SOLUTION PARAMETERS
C...
C...     FRONT FLATTENING
         IF(NORUN.EQ.1)THEN
         A=1.0D0
         B=1.0D0
         C=0.0D0
         NT=11
C...
C...     FRONT SHARPENING
         ELSE
     +   IF(NORUN.EQ.2)THEN
         A= 1.0D0
         B=-1.0D0
         C= 1.0D0
         NT=11
         END IF
C...
C...  NUMBER OF GRID POINTS
      N=101
C...
C...  TOTAL LENGTH IN X, GRID SPACING
      XL=1.0D0
      DX=XL/DFLOAT(N-1)
C...
C...  INITIAL CONDITION
      DO 1 I=1,N
         X(I)=XL*DFLOAT(I-1)/DFLOAT(N-1)
         U(I)=ANAL(X(I),0.0D0)
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  SELECT DERIVATIVE SUROUTINE
      IF(NCASE.EQ.1)CALL DERV1
      IF(NCASE.EQ.2)CALL DERV2
      IF(NCASE.EQ.3)CALL DERV3
      RETURN
      END
      SUBROUTINE DERV1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=101)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NX)
     +      /F/ UT(NX)
     +      /S/ UX(NX)
     +      /X/  X(NX)
     +      /C/      A,      B,      C,     XL,     DX
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  BOUNDARY CONDITION AT X = 0
       U(1)=ANAL(0.0D0,T)
      UT(1)=0.0D0
C...
C...  PDE
      DO 1 I=2,N
         UT(I)=-U(I)*(U(I)-U(I-1))/DX
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV2
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=101)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NX)
     +      /F/ UT(NX)
     +      /S/ UX(NX)
     +      /X/  X(NX)
     +      /C/      A,      B,      C,     XL,     DX
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  BOUNDARY CONDITION AT X = 0
       U(1)=ANAL(0.0D0,0.0D0)
      UT(1)=0.0D0
C...
C...  PDE
      DO 1 I=2,N
         UT(I)=-U(I)*(U(I)-U(I-1))/DX
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV3
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=101)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NX)
     +      /F/ UT(NX)
     +      /S/ UX(NX)
     +      /X/  X(NX)
     +      /C/      A,      B,      C,     XL,     DX
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  BOUNDARY CONDITION AT X = 0
       U(1)=ANAL(0.0D0,0.0D0)
      UT(1)=0.0D0
C...
C...  PDE
C...  X NE 0 OR XL
      DO 1 I=2,N-1
         UT(I)=-U(I)*(U(I+1)-U(I-1))/(2.0D0*DX)
1     CONTINUE
C...
C...  X = XL
      UT(N)=-U(N)*(U(N)-U(N-1))/DX
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  SELECT PRINT ROUTINE
      IF(NCASE.EQ.1)CALL PRINT1(NI,NO)
      IF(NCASE.EQ.2)CALL PRINT2(NI,NO)
      IF(NCASE.EQ.3)CALL PRINT3(NI,NO)
      RETURN
      END
      SUBROUTINE PRINT1(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=101)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NX)
     +      /F/ UT(NX)
     +      /S/ UX(NX)
     +      /X/  X(NX)
     +      /C/      A,      B,      C,     XL,     DX
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  MONITOR OUTPUT ON SCREEN
      WRITE(*,*)NORUN,T
C...
C...  PRINT NUMERICAL AND ANALYTICAL SOLUTIONS
      WRITE(NO,3)T
3        FORMAT(/,' t = ',F5.2,/,9X,'x',4X,'u(x,t)',3X,'ua(x,t)')
         DO 1 I=1,N
         UA=ANAL(X(I),T)
         WRITE(NO,2)X(I),U(I),UA
2        FORMAT(F10.2,2F10.5)
1        CONTINUE
C...
C...  PLOT THE SOLUTION
      IP=IP+1
      CALL PLOTM
      RETURN
      END
      SUBROUTINE PRINT2(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=101)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NX)
     +      /F/ UT(NX)
     +      /S/ UX(NX)
     +      /X/  X(NX)
     +      /C/      A,      B,      C,     XL,     DX
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  MONITOR OUTPUT ON SCREEN
      WRITE(*,*)NORUN,T
C...
C...  PRINT NUMERICAL SOLUTION
      WRITE(NO,3)T
3        FORMAT(/,' t = ',F5.2,/,9X,'x',4X,'u(x,t)')
         DO 1 I=1,N
         WRITE(NO,2)X(I),U(I)
2        FORMAT(F10.2,F10.5)
1        CONTINUE
C...
C...  PLOT THE SOLUTION
      IP=IP+1
      CALL PLOTM
      RETURN
      END
      SUBROUTINE PRINT3(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=101)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NX)
     +      /F/ UT(NX)
     +      /S/ UX(NX)
     +      /X/  X(NX)
     +      /C/      A,      B,      C,     XL,     DX
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  MONITOR OUTPUT ON SCREEN
      WRITE(*,*)NORUN,T
C...
C...  PRINT NUMERICAL SOLUTION
      WRITE(NO,3)T
3        FORMAT(/,' t = ',F5.2,/,9X,'x',4X,'u(x,t)')
         DO 1 I=1,N
         WRITE(NO,2)X(I),U(I)
2        FORMAT(F10.2,F10.5)
1        CONTINUE
C...
C...  PLOT THE SOLUTION
      IP=IP+1
      CALL PLOTM
      RETURN
      END
      DOUBLE PRECISION FUNCTION ANAL(X,T)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/C/      A,      B,      C,     XL,     DX
      ANAL=(C+B*X)/(A+B*T)
      RETURN
      END
      SUBROUTINE PLOTM
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NX=101)
      COMMON/T/      T,  NSTOP,  NORUN
     +      /Y/  U(NX)
     +      /F/ UT(NX)
     +      /S/ UX(NX)
     +      /X/  X(NX)
     +      /C/      A,      B,      C,     XL,     DX
     +      /I/      N,  NCASE,     IP,     NT
C...
C...  OPEN A FILE FOR MATLAB PLOTTING
      IF((IP.EQ.1).AND.(NORUN.EQ.1))THEN
         OPEN(1,FILE='m3.out')
      END IF
C...
C...  WRITE THE NUMERICAL SOLUTION FOR SUBSEQUENT MATLAB PLOTTING
      DO 2 I=1,N
         WRITE(1,3)X(I),U(I)
3        FORMAT(F10.3,F10.5)
2     CONTINUE
      RETURN
      END
