      SUBROUTINE INITAL
C...
C...  SOLUTION OF BURGERS' EQUATION BY SPLINE APPROXIMATION ON AN
C...  ADAPTIVE GRID
C...
C...  SINGLE PRECISION
C...
C...  SUBROUTINE INITAL
C...
C...     (1)  DEFINES THE INITIAL ADAPTIVE GRID AS THE BASIC UNIFORM
C...          GRID VIA A CALL TO SUBROUTINE GRIDE
C...
C...     (2)  DEFINES THE INITIAL CONDITION ON THE ADAPTIVE GRID OF
C...          (1) BY A CALL TO FUNCTION PHI WHICH IMPLEMENTS THE EXACT
C...          SOLUTION TO BURGERS EQUATION
C...
C...  COMMON/SYSTM1/ IS ACCESSED SO THE NUMBER OF ODES DEFINED BY THE
C...  ADAPTIVE GRID, N, CAN BE CHANGED AS THE GRID EVOLVES
      COMMON/SYSTM1/ T0,TF,TP,H,ERROR,N
C...
C...  THE FOLLOWING COMMON AREAS ARE BRIEFLY EXPLAINED BELOW
C...
C...     /T/ /Y/ /F/    STANDARD DSS/2 COMMON AREAS DIMENSIONED FOR UP
C...                    TO 150 ODES
C...
C...     /B/            THE SOLUTION, CB, AND ITS SECOND DERIVATIVE,
C...                    CBXX, DEFINED ON THE BASIC GRID
C...
C...     /SV/           INDEPENDENT VARIABLE DEFINED ON THE ADAPTIVE
C...                    GRID
C...
C...     /SC/           CUBIC SPLINE COEFFICIENTS (LINEAR, QUADRATIC
C...                    AND CUBIC COEFFICIENTS RESPECTIVELY)
C...
C...     /C/            CONSTANTS PERTAINING TO THE PROBLEM
C...
C...     /AB1/          PARAMETERS WHICH CONTROL THE ADAPTIVE GRID
C...
      COMMON  /T/         T
     1        /Y/   CA(150)
     2        /F/   CT(150)
     3        /B/   CB(120), CBXX(120)
     4       /SV/   XA(150)
     5       /SC/   CX(150),  CXX(150), CXXX(150)
     6        /C/        XL,        XU,         V,       VIS,
     7                   SN,        NB,        IP,        JP
      COMMON/AB1/        DL,       XFL,       XFR,        NF,
     1                  DXF,   XB(120),  LXB(120),  NOB(120)
C...
C...  SET THE VELOCITY AND VISCOSITY OF BURGERS EQUATION
      V=1.0E+00
      VIS=0.003E+00
C...
C...  INITIALIZE THE COUNTERS FOR PRINTED AND PLOTTED OUTPUT USED
C...  IN SUBROUTINE PRINT, THE RUNNING SUM OF THE NUMBER OF ODES
C...  USED IN THE ADAPTIVE GRID
      IP=0
      JP=0
      SN=0.0E+00
C...
C...  DEFINE THE LENGTH OF THE X AXIS
      XL=0.0E+00
      XU=1.0E+00
C...
C...  DEFINE THE NUMBER OF BASIC GRID POINTS, NUMBER OF SUBINTERVALS
C...  BETWEEN THE BASIC GRID POINTS WHEN ADAPTIVE POINTS ARE ADDED
      NB=41
      NF=2
C...
C...  DEFINE AN INITIAL UNIFORM GRID OF NB POINTS
      CALL GRIDE(NB,XL,XU,DXB,XB)
C...
C...  THE NUMBER OF ADAPTIVE GRID POINTS, N, INITIALLY EQUALS THE
C...  NUMBER OF BASIC GRID POINTS, NB
      N=NB
C...
C...  DEFINE THE ADAPTIVE GRID INTERVAL, DXF, FROM THE BASIC GRID
C...  INTERVAL, DXB
      DXF=DXB/FLOAT(NF)
C...
C...  DEFINE THE THRESHOLD FOR THE SECOND DERIVATIVE, ABOVE WHICH
C...  ADAPTIVE GRID POINTS ARE ADDED
      DL=20.0E+00
C...
C...  DEFINE THE MAXIMUM DISTANCES TO THE LEFT AND RIGHT OF A BASIC
C...  GRID POINT, XFL AND XFR, BEYOND WHICH ADAPTIVE GRID POINTS WILL
C...  NOT BE ADDED
      SFL=1.0E+00
      SFR=1.5E+00
      XFL=SFL*V*TP
      XFR=SFR*V*TP
C...
C...  DEFINE THE INITIAL BASIC AND ADAPTIVE GRIDS, EACH WITH NB POINTS
      DO 1 I=1,NB
C...
C...  INITIALIZE ARRAY LXB WHICH SPECIFIES THE EXISTENCE OF ADAPTIVE
C...  GRID POINTS
      LXB(I)=0
C...
C...  STORE THE SUBSCRIPTS OF THE BASIC GRID IN ARRAY NOB
      NOB(I)=I
C...
C...  DEFINE THE INITIAL CONDITIONS OF THE BASIC AND ADAPTIVE GRIDS
C...  FROM THE EXACT SOLUTION TO BURGERS EQUATION
      CB(I)=PHI(0.0E+00,XB(I))
      CA(I)=CB(I)
      XA(I)=XB(I)
1     CONTINUE
C...
C...  INITIALIZE THE DERIVATIVE CALCULATIONS IN SUBROUTINE DERV
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
C...
C...  SUBROUTINE DERV
C...
C...     (1)  COMPUTES THE FIRST AND SECOND DERIVATIVES, UX AND UXX,
C...          VIA A CALL TO SUBROUTINE NCSPLE
C...
C...     (2)  USES DERIVATIVES UX AND UXX TO ASSEMBLE BURGERS EQUATION
C...          IN THE USUAL METHOD OF LINES FORMAT
C...
      COMMON/SYSTM1/ T0,TF,TP,H,ERROR,N
      COMMON  /T/         T
     1        /Y/   CA(150)
     2        /F/   CT(150)
     3        /B/   CB(120), CBXX(120)
     4       /SV/   XA(150)
     5       /SC/   CX(150),  CXX(150), CXXX(150)
     6        /C/        XL,        XU,         V,       VIS,
     7                   SN,        NB,        IP,        JP
      COMMON/AB1/        DL,       XFL,       XFR,        NF,
     1                  DXF,   XB(120),  LXB(120),  NOB(120)
C...
C...  APPLY THE BOUNDARY CONDITIONS FOR BURGERS EQUATION AT X = 0
C...  AND X = 1
      CA(1)=PHI(T,0.0E+00)
      CA(N)=PHI(T,1.0E+00)
C...
C...  COMPUTE THE FIRST AND SECOND SPATIAL DERIVATIVES IN BURGERS
C...  EQUATION, UX AND UXX, BY THE NORMAL CUBIC SPLINE
      CALL NCSPLE(N,XA,CA,CX,CXX,CXXX)
C...
C...  ASSEMBLE BURGERS EQUATION OVER GRID POINTS 2 TO N-1.  COMMON
C...  /SC/ ACTUALLY CONTAINS THE COEFFICIENTS OF THE CUBIC SPLINE
C...  C(X) = A0 + A1*(X - XI) + A2*(X - XI)**2 + A3*(X - XI)***3,
C...  I.E., COMMON/SC/A1(150),A2(150),A3(150), AFTER NCSPLE IS
C...  CALLED.  THUS CX(XI) = A1, CXX(XI) = 2*A2, CXXX(XI) = 6*A3,
C...  AND THE SECOND ARRAY IN /SC/ MUST BE MULTIPLIED BY 2 TO OBTAIN
C...  THE SECOND DERIVATIVE, CXX, IN BURGERS EQUATION.  ALSO, FROM
C...  THE EQUATION FOR THE SPLINE, C(XI) = A0 = CA(XI) WHERE ARRAY
C...  CA IS IN COMMON/Y/.  SUBROUTINE NCSPLE THEREFORE NUMERICALLY
C...  DIFFERENTIATES CA TO CX, CXX AND CXXX AS WELL AS EVALUATING
C...  THE COEFFICIENTS OF THE CUBIC SPLINE (THESE TWO OPERATIONS
C...  ARE EQUIVALENT)
      NM1=N-1
      DO 1 I=2,NM1
      CT(I)=2.0E+00*VIS*CXX(I)-CA(I)*CX(I)
1     CONTINUE
C...
C...  ZERO THE TIME DERIVATIVES AT I = 1, N TO INSURE THE BOUNDARY
C...  CONDITIONS ARE OBSERVED
      CT(1)=0.0E+00
      CT(N)=0.0E+00
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
C...
C...  SUBROUTINE PRINT
C...
C...     (1)  STORES THE NUMERICAL AND EXACT SOLUTIONS TO BURGERS
C...          EQUATION DURING SUCCESSIVE CALLS FOR PLOTTING
C...
C...     (2)  REDEFINES THE ADAPTIVE GRID IN ACCORDANCE WITH CHANGES
C...          IN THE NUMERICAL SOLUTION VIA A CALL TO SUBROUTINE ANUGB1
C...
C...     (3)  PLOTS THE NUMBER OF POINTS IN THE ADAPTIVE GRID VS T AT
C...          THE END OF THE RUN AND PRINTS THE AVERAGE NUMBER OF GRID
C...          POINTS DURING THE RUN
C...
      COMMON/SYSTM1/ T0,TF,TP,H,ERROR,N
      COMMON  /T/         T
     1        /Y/   CA(150)
     2        /F/   CT(150)
     3        /B/   CB(120), CBXX(120)
     4       /SV/   XA(150)
     5       /SC/   CX(150),  CXX(150), CXXX(150)
     6        /C/        XL,        XU,         V,       VIS,
     7                   SN,        NB,        IP,        JP
      COMMON/AB1/        DL,       XFL,       XFR,        NF,
     1                  DXF,   XB(120),  LXB(120),  NOB(120)
C...
C...  DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION
      DIMENSION     GR(150),  CAA(150),    PN(51),    TN(51),
     1          CBXXP(2,51)
C...
C...  INCREMENT THE INDEX FOR THE PLOTTED SOLUTION
      IP=IP+1
C...
C...  COMPUTE THE RUNNING SUM OF THE NUMBER OF ODES
      SN=SN+FLOAT(N)
C...
C...  STORE THE NUMBER OF ODES FOR SUBSEQUENT PLOTTING
      PN(IP)=FLOAT(N)
      TN(IP)=T
C...
C...  CALL SUBROUTINE NCSPLE TO UPDATE THE SPLINE COEFFICIENTS USED
C...  IN THE NEW ADAPTIVE GRID
      CALL NCSPLE(N,XA,CA,CX,CXX,CXXX)
C...
C...  STORE THE DEPENDENT VARIABLE AND ITS SECOND DERIVATIVE ON THE
C...  ADAPTIVE GRID, CA(I) AND CXX(I), AT THE BASIC GRID POINTS WITH
C...  SUBSCRIPT K IN PREPARATION FOR UPDATING THE ADAPTIVE GRID.  THE
C...  FACTOR OF 2 MULTIPLYING CXX(K) IS EXPLAINED IN SUBROUTINE DERV
      DO 1 I=1,NB
      K=NOB(I)
      CB(I)=CA(K)
      CBXX(I)=2.0E+00*CXX(K)
1     CONTINUE
C...
C...  STORE THE MINIMUM AND MAXIMUM VALUES OF THE SECOND DERIVATIVE,
C...  CBXX, FOR SUBSEQUENT PLOTTING
      CBMAX=ABS(CBXX(1))
      CBMIN=ABS(CBXX(1))
      DO 10 I=2,NB
      CBMAX=AMAX1(ABS(CBXX(I)),CBMAX)
      CBMIN=AMIN1(ABS(CBXX(I)),CBMIN)
10    CONTINUE
      CBXXP(1,IP)=CBMAX
      CBXXP(2,IP)=CBMIN
C...
C...  STORE THE EXACT SOLUTION FOR SUBSEQUENT PLOTTING
      IF((IP-6)*(IP-16)*(IP-26)*(IP-36)*(IP-46).NE.0)GO TO 3
      JP=JP+1
      DO 2 I=1,N
      CAA(I)=PHI(T,XA(I))
      GR(I)=-0.05E+00*FLOAT(5-JP)
2     CONTINUE
C...
C...  PLOT THE NUMERICAL SOLUTION WITH THE ADAPTIVE GRID SUPERIMPOSED
C...  AT SUCCESSIVE VALUES OF THE INDEPENDENT VARIABLE T.  INTEGER
C...  JP SPECIFIES THAT ONLY ADDITIONAL CURVES ARE PLOTTED AFTER THE
C...  FIRST CALL TO QIKPLT
      IM=-N
      IF(JP.GT.1)GO TO 4
C...
C...  ******************************************************************
C...  THE FOLLOWING CODING CALLS A CONTINUOUS (CALCOMP) PLOTTER.  IT
C...  IS NOT STANDARD FORTRAN, AND THEREFORE WILL PROBABLY HAVE TO BE
C...  CHANGED FOR OTHER COMPUTERS AND PLOTTERS
C...  CALL QIKSET(4.,0.,0.25,7.,-0.3,0.2)
C...  CALL QIKSAX(3,3)
C...  CALL QIKFORM(2,2)
C...  CALL QIKPLT(XA,CAA,N,3H*X*,8H*U(X,T)*,18H*BURGERS EQUATION*,1)
C...  CALL PLOT(-5.,1.,-3)
C...  CALL QLINE(XA,CA,IM,4)
C...  CALL QLINE(XA,GR,IM,-13)
      GO TO 6
4     CONTINUE
C...  CALL QLINE(XA,CAA,N,M0)
C...  CALL QLINE(XA,CA,IM,4)
C...  CALL QLINE(XA,GR,IM,-13)
C...  ******************************************************************
C...
C...  PRINT THE NUMERICAL SOLUTION CORRESPONDING TO THE PLOTTED SOLUTION
6     WRITE(NO,7)T,N,(XA(I),I=1,N)
7     FORMAT(1H ,//,5H T = ,E11.3,3X,5H N = ,I3,/,3H  X,/,(10E11.3))
      WRITE(NO,8)(CAA(I),I=1,N)
8     FORMAT(1H ,7H U ANAL,/,(10E11.3))
      WRITE(NO,9)( CA(I),I=1,N)
9     FORMAT(1H ,7H U NUM ,/,(10E11.3))
C...
C...  CONTINUE WITH SUBROUTINE PRINT, PARTICULARLY THE REDEFINITION OF
C...  THE ADAPTIVE GRID
3     CONTINUE
C...
C...  CALL SUBROUTINE ANUGB1 TO REDEFINE THE ADAPTIVE GRID
      CALL ANUGB1(NB,CB,CBXX,N,XA,CA)
C...
C...  AT THE END OF THE RUN, PLOT THE NUMBER OF ODES VS T
      IF(IP.LT.51)RETURN
      CALL TPLOTS(1,IP,TN,PN)
C...
C...  COMPUTE AND PRINT THE AVERAGE NUMBER OF ODES
      AN=SN/FLOAT(IP)
      WRITE(NO,5)AN
5     FORMAT(1H ,//,
     1 7H N VS T,3X,8HAVG N = ,F6.2)
C...
C...  PLOT THE NUMBER OF ODES, THE MINIMUM AND MAXIMUM VALUES OF THE
C...  SECOND DERIVATIVE VS T
C...  ******************************************************************
C...  THE FOLLOWING CODING CALLS A CONTINUOUS (CALCOMP) PLOTTER.  IT
C...  IS NOT STANDARD FORTRAN, AND THEREFORE WILL PROBABLY HAVE TO BE
C...  CHANGED FOR OTHER COMPUTERS AND PLOTTERS
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(4.,0.,0.25,7.,40.,2.5)
C...  CALL QIKPLT(TN,PN,IP,3H*T*,3H*N*,16H*NUMBER OF ODES*,1)
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(4.,0.,0.25,7.,0.,400.)
C...  CALL QIKPLT(TN,CBXXP,IP,3H*T*,6H*CBXX*,19H*SECOND DERIVATIVE*,2)
C...  ******************************************************************
C...
      RETURN
      END
      REAL FUNCTION PHI(T,X)
C...
C...  FUNCTION PHI(T,X) COMPUTES THE EXACT SOLUTION TO BURGERS
C...  EQUATION.  IT IS USED IN THREE PLACES
C...
C...     (1)  TO PROVIDE THE EXACT SOLUTION IN SUBROUTINE PRINT FOR
C...          COMPARISON WITH THE NUMERICAL SOLUTION.
C...
C...     (2)  TO PROVIDE THE INITIAL CONDITION IN SUBROUTINE INITAL.
C...
C...     (3)  TO PROVIDE THE BOUNDARY CONDITIONS IN SUBROUTINE DERV.
C...
C...  ARGUMENT LIST
C...
C...      T   INITIAL-VALUE INDEPENDENT VARIABLE IN BURGERS EQUATION
C...          (INPUT)
C...
C...      X   BOUNDARY-VALUE INDEPENDENT VARIABLE IN BURGERS EQUATION
C...          (INPUT)
C...
C...  THE VALUE OF THE FUNCTION (PHI) IS THE EXACT SOLUTION TO BURGERS
C...  EQUATION AT X AND T.
C...
C...  TYPE REAL VARIABLES AS SINGLE PRECISION
C...
      COMMON  /C/        XL,        XU,         V,       VIS
      A=(0.05E+00/VIS)*(X-0.5E+00+4.95E+00*T)
      B=(0.25E+00/VIS)*(X-0.5E+00+0.75E+00*T)
      C=( 0.5E+00/VIS)*(X-0.375E+00)
      EA=EXP(-A)
      EB=EXP(-B)
      EC=EXP(-C)
C...
C...  THE FOLLOWING IF WAS ADDED FOR SHORT WORD LENGTH COMPUTERS, E.G.,
C...  32-BIT COMPUTERS SUCH AS THE VAX, WHICH CANNOT ACCOMMODATE WIDE
C...  VARIATIONS IN THE EXP FUNCTION.  THE INTENTION IS TO EFFECTIVELY
C...  AVOID A DIVISION BY ZERO.  IT IS BASED ON THE OBSERVATION THAT
C...  WHEN THE DIVISION BY ZERO OCCURS, EC MLT EB MLT EA MLT (WHERE
C...  MLT DENOTES "MUCH LESS THAN") SO THE EQUATION FOR PHI BECOMES
C...
C...     PHI = 0.1*EA/EA = 0.1
C...
      IF(ABS(EA+EB+EC).LT.1.0E-20)THEN
         PHI=0.1E+00
      ELSE
         PHI=(0.1E+00*EA+0.5E+00*EB+EC)/(EA+EB+EC)
      END IF
      RETURN
      END
