C     ALGORITHM 621 COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 4,
C     DEC., 1984, P. 378.
C     PROGRAM EXAM1(OUTPUT,TAPE6=OUTPUT)
C-----------------------------------------------------------------
C   TEST PROGRAM 1 FOR THE ALGORITHM BDMG.
C   THIS TEST PROBLEM IS A TWO-DIMENSIONAL NONLINEAR PARABOLIC
C   EQUATION DISCUSSED IN SECTION 5 ("FIRST EXAMPLE") OF THE
C   COMPANION PAPER "SOFTWARE WITH LOW STORAGE REQUIREMENTS FOR
C   TWO-DIMENSIONAL NONLINEAR PARABOLIC DIFFERENTIAL EQUATIONS" BY
C   B.P. SOMMEIJER AND P.J. VAN DER HOUWEN, WHERE ALSO THE
C   NUMERICAL RESULTS OF THIS TEST CAN BE FOUND.
C-----------------------------------------------------------------
      DIMENSION V(25,25),X(49),Y(49),INFO(18)
      DIMENSION WORK(3850),IWORK(43)
      EXTERNAL G,SPRAD
C-----------------------------------------------------------
C   DEFINITION OF THE FINEST GRID AND OF THE NUMBER OF GRIDS
C-----------------------------------------------------------
      NX = 25
      NY = 25
      DX = 1.0/(NX - 1.0)
      DO 10 I = 1,NX
   10    X(I) = (I - 1)*DX
      DY = 1.0/(NY - 1.0)
      DO 20 I = 1,NY
   20    Y(I) = (I - 1)*DY
      M = 4
C------------------------------------------------------------
C   DEFINITION OF THE INTEGRATION INTERVAL AND INITIALIZATION
C------------------------------------------------------------
      T = 0.0
      TEND = 1.0
      DO 30 J = 1,NY
         DO 30 I = 1,NX
   30    V(I,J) = SOL (T, X(I), Y(J))
C---------------------------------------
C   DEFINITION OF THE PROBLEM PARAMETERS
C---------------------------------------
      TOL = 1.0E-4
      INFO(1) = 0
      INFO(2) = 2000
      INFO(3) = 3
      IWORK(43) = 3850
      METH = 1
C--------------------------
C   CALL FOR THE INTEGRATOR
C--------------------------
      CALL BDMG (NX, NY, M, X, Y,
     +           T, TEND, V, G, SPRAD, TOL, METH,
     +           WORK, IWORK, INFO, IFLAG)
C----------------------------
C   CHECK ERRORFLAG ON RETURN
C----------------------------
      IF (IFLAG .NE. 0) WRITE (6,1000) IFLAG
      IF (IFLAG .GE. 1 .AND. IFLAG .LE. 4) GO TO 60
C-------------------------
C   OUTPUT SOME STATISTICS
C-------------------------
      WRITE (6,1001) INFO(7), INFO(8), INFO(9)
      WRITE (6,1002) INFO(10)
C---------------------------------------------
C   PRINT THE SOLUTION AND DETERMINE THE ERROR
C---------------------------------------------
      WRITE (6,1003) T
      NXM1 = NX - 1
      NYM1 = NY - 1
      AE = -1.0
      DO 50 J = 2,NYM1
         WRITE (6,1004) J
         DO 40 I = 2,NXM1
            AE = AMAX1(AE,ABS(SOL(T, X(I), Y(J))-V(I,J)))
   40       CONTINUE
   50    WRITE (6,1005) (V(I,J),I = 2,NXM1)
      WRITE (6,1006) -ALOG10(AE)
 1000 FORMAT(27H ERRORFLAG HAS BEEN SET TO:,I5)
 1001 FORMAT(27H NUMBER OF STEPS PERFORMED:,I6/
     +       26H NUMBER OF REJECTED STEPS:,I6/
     +       31H NUMBER OF STEPS WITH OVERFLOW:,I6)
 1002 FORMAT(44H TOTAL NUMBER OF (FINEST GRID)G-EVALUATIONS:,I6)
 1003 FORMAT(//15H SOLUTION AT T=,F10.4/)
 1004 FORMAT(I4)
 1005 FORMAT(5(E20.10))
 1006 FORMAT(34H MINIMAL NUMBER OF CORRECT DIGITS:,F7.2)
C------------------------------------
C   END OF MAIN PROGRAM FOR EXAMPLE 1
C------------------------------------
   60 STOP
      END
      FUNCTION SOL (T, X, Y)
C-----------------------------------------------------
C   ANALYTIC FUNCTION USED FOR INITIALIZATION,BOUNDARY
C   CONDITIONS AND CALCULATION OF THE ERROR
C-----------------------------------------------------
      SOL = (0.8*(2.0*T + X+ Y ))**0.25
      RETURN
      END
      FUNCTION SPRAD (K, NXK, NYK, XK, YK, T, VK, TAU)
      DIMENSION XK(NXK),YK(NYK),VK(NXK,NYK)
C---------------------------------------------------------
C   DEFINE AN UPPER ESTIMATE OF THE SPECTRAL RADIUS ON THE
C   INTERVAL [T,T+TAU]
C---------------------------------------------------------
      SPRAD = 32.0*(1.0 + T + TAU)*((NXK - 1)**2 + (NYK - 1)**2)
      RETURN
      END
      SUBROUTINE G (K, NXK, NYK, XK, YK, T, VK, DVK)
      DIMENSION XK(NXK),YK(NYK),VK(NXK,NYK),DVK(NXK,NYK)
C----------------------------------------------
C   DEFINE THE DERIVATIVE FUNCTION G AT LEVEL K
C----------------------------------------------
      NXKM1 = NXK - 1
      NYKM1 = NYK - 1
C-----------------------------------
C   TREATMENT OF THE BOUNDARY POINTS
C-----------------------------------
      DO 10 J = 1,NYK
         DO 10 I = 1,NXK,NXKM1
         DVK(I,J) = 0.0
   10    VK(I,J) = SOL (T, XK(I), YK(J))
      DO 20 J = 1,NYK,NYKM1
         DO 20 I = 2,NXKM1
         DVK(I,J) = 0.0
   20    VK(I,J) = SOL (T, XK(I), YK(J))
C---------------------------------
C   DERIVATIVES AT INTERNAL POINTS
C---------------------------------
      DO 30 J = 2,NYKM1
         DO 30 I = 2,NXKM1
   30    DVK(I,J)=(VK(I+1,J)**5-2.0*VK(I,J)**5+VK(I-1,J)**5)*NXKM1**2 +
     +            (VK(I,J+1)**5-2.0*VK(I,J)**5+VK(I,J-1)**5)*NYKM1**2
      RETURN
      END
C     PROGRAM EXAM2 (OUTPUT,TAPE6=OUTPUT)
C-----------------------------------------------------------------
C   TEST PROGRAM 2 FOR THE ALGORITHM BDMG.
C   THIS TEST PROBLEM IS A TWO-DIMENSIONAL NONLINEAR PARABOLIC
C   EQUATION INCLUDING A MIXED-DERIVATIVE TERM, DISCUSSED IN
C   IN SECTION 5 ("SECOND EXAMPLE") OF THE COMPANION PAPER
C   "SOFTWARE WITH LOW STORAGE REQUIREMENTS FOR TWO-DIMENSIONAL
C   NONLINEAR PARABOLIC DIFFERENTIAL EQUATIONS" BY B.P. SOMMEIJER
C   AND P.J. VAN DER HOUWEN, WHERE ALSO THE NUMERICAL RESULTS OF
C   THIS TEST CAN BE FOUND.
C-----------------------------------------------------------------
      DIMENSION V(17,17),X(31),Y(31),INFO(16)
      DIMENSION WORK(1776),IWORK(33)
      EXTERNAL G,SPRAD
C-----------------------------------------------------------
C   DEFINITION OF THE FINEST GRID AND OF THE NUMBER OF GRIDS
C-----------------------------------------------------------
      NX = 17
      NY = 17
      DX = 1.0/(NX - 1.0)
      DO 10 I = 1,NX
   10    X(I) = (I - 1)*DX
      DY = 1.0/(NY - 1.0)
      DO 20 I = 1,NY
   20    Y(I) = (I - 1)*DY
      M = 3
C---------------------------------------------------------------
C   SUCCESSIVE INTEGRATION PROCESSES FOR DIFFERENT VALUES OF TOL
C---------------------------------------------------------------
      DO 60 ITOL = 1,6
C------------------------------------------------------------
C   DEFINITION OF THE INTEGRATION INTERVAL AND INITIALIZATION
C------------------------------------------------------------
         T = 0.0
         TEND = 1.0
         DO 30 J = 1,NY
            DO 30 I = 1,NX
   30       V(I,J) = SOL (T, X(I), Y(J))
C---------------------------------------
C   DEFINITION OF THE PROBLEM PARAMETERS
C---------------------------------------
         TOL = 10.0**(-ITOL)
         INFO(1) = 0
         INFO(2) = 1000
         INFO(3) = 3
         IWORK(33) = 1776
         METH = 1
C--------------------------
C   CALL FOR THE INTEGRATOR
C--------------------------
         CALL BDMG (NX, NY, M, X, Y,
     +              T, TEND, V, G, SPRAD, TOL, METH,
     +              WORK, IWORK, INFO, IFLAG)
C----------------------------
C   CHECK ERRORFLAG ON RETURN
C----------------------------
         IF (IFLAG .NE. 0) WRITE (6,1000) IFLAG
         IF (IFLAG .GE. 1 .AND. IFLAG .LE. 4) GO TO 60
C-------------------------
C   OUTPUT SOME STATISTICS
C-------------------------
         WRITE (6,1007) TOL
         WRITE (6,1001) INFO(7), INFO(8), INFO(9)
         WRITE (6,1002) INFO(10)
C---------------------------------------------
C   PRINT THE SOLUTION AND DETERMINE THE ERROR
C---------------------------------------------
         WRITE (6,1003) T
         NXM1 = NX - 1
         NYM1 = NY - 1
         AE = -1.0
         DO 50 J = 2,NYM1
            WRITE (6,1004) J
            DO 40 I = 2,NXM1
               AE = AMAX1(AE,ABS(SOL (T, X(I), Y(J))-V(I,J)))
   40          CONTINUE
   50       WRITE (6,1005) (V(I,J),I = 2,NXM1)
         WRITE (6,1006) -ALOG10(AE)
   60    CONTINUE
 1000 FORMAT(27H ERRORFLAG HAS BEEN SET TO:,I5)
 1001 FORMAT(27H NUMBER OF STEPS PERFORMED:,I6/
     +       26H NUMBER OF REJECTED STEPS:,I6/
     +       31H NUMBER OF STEPS WITH OVERFLOW:,I6)
 1002 FORMAT(44H TOTAL NUMBER OF (FINEST GRID)G-EVALUATIONS:,I6)
 1003 FORMAT(//15H SOLUTION AT T=,F10.4/)
 1004 FORMAT(I4)
 1005 FORMAT(5(E20.10))
 1006 FORMAT(34H MINIMAL NUMBER OF CORRECT DIGITS:,F7.2)
 1007 FORMAT(1H1,7H TOL = ,E10.2)
C------------------------------------
C   END OF MAIN PROGRAM FOR EXAMPLE 2
C------------------------------------
      STOP
      END
      FUNCTION SOL (T, X, Y)
C-----------------------------------------------------
C   ANALYTIC FUNCTION USED FOR INITIALIZATION,BOUNDARY
C   CONDITIONS AND CALCULATION OF THE ERROR
C-----------------------------------------------------
      SOL = EXP(-T)*X*Y*(X + Y)
      RETURN
      END
      FUNCTION SPRAD (K, NXK, NYK, XK, YK, T, VK, TAU)
      DIMENSION XK(NXK),YK(NYK),VK(NXK,NYK)
C---------------------------------------------------------
C   DEFINE AN UPPER ESTIMATE OF THE SPECTRAL RADIUS ON THE
C   INTERVAL [T,T+TAU]
C---------------------------------------------------------
      SPRAD = 6.0*((NXK-1)**2 + (NYK-1)**2) + 2.0*(NXK-1)*(NYK-1)
      RETURN
      END
      SUBROUTINE G (K, NXK, NYK, XK, YK, T, VK, DVK)
      DIMENSION XK(NXK),YK(NYK),VK(NXK,NYK),DVK(NXK,NYK)
C----------------------------------------------
C   DEFINE THE DERIVATIVE FUNCTION G AT LEVEL K
C----------------------------------------------
      NXKM1 = NXK - 1
      NYKM1 = NYK - 1
C-----------------------------------
C   TREATMENT OF THE BOUNDARY POINTS
C-----------------------------------
      DO 10 J = 1,NYK
         DO 10 I = 1,NXK,NXKM1
         DVK(I,J) = 0.0
   10    VK(I,J) = SOL (T, XK(I), YK(J))
      DO 20 J = 1,NYK,NYKM1
         DO 20 I = 2,NXKM1
         DVK(I,J) = 0.0
   20    VK(I,J) = SOL (T, XK(I), YK(J))
C---------------------------------
C   DERIVATIVES AT INTERNAL POINTS
C---------------------------------
      ET = EXP(-T)
      DO 30 J = 2,NYKM1
         Y2 = YK(J)**2
         DO 30 I = 2,NXKM1
         X2 = XK(I)**2
         A = 0.5*X2 + Y2
         B = -0.5*(X2 + Y2)
         C = X2 + 0.5*Y2
         D = 1.0/(1.0 + (X2*YK(J) + XK(I)*Y2)*ET)
         DVK(I,J) = A*(VK(I-1,J)-2.0*VK(I,J)+VK(I+1,J))*NXKM1*NXKM1 +
     +              C*(VK(I,J-1)-2.0*VK(I,J)+VK(I,J+1))*NYKM1*NYKM1 +
     +              (VK(I+1,J+1)-VK(I-1,J+1)+VK(I-1,J-1)-VK(I+1,J-1))*
     +              2.0*B*NXKM1*NYKM1/4.0
   30    DVK(I,J) = DVK(I,J)*(D*(1.0 + VK(I,J)))**10
      RETURN
      END
      SUBROUTINE BDMG (NX, NY, M, X, Y,
     *                 T, TEND, V, G, SPRAD, TOL, METH,
     *                 WORK, IWORK, INFO, IFLAG)
C
C-----------------------------------------------------------------------
C
C   PURPOSE AND METHOD
C
C-----------------------------------------------------------------------
C
C       BDMG IS A TIME-INTEGRATOR DESIGNED TO SOLVE A SYSTEM OF ORDINARY
C    DIFFERENTIAL EQUATIONS(ODE'S) ORIGINATING FROM SEMI-DISCRETIZATION
C    OF A (SCALAR) GENERAL PARABOLIC PARTIAL DIFFERENTIAL EQUATION(PDE)
C    IN TWO SPACE DIMENSIONS. THIS SEMI-DISCRETIZATION CAN BE PERFORMED
C    EITHER BY HAND OR BY PDETWO[3], AN INTERFACE TO FORM AND EVALUATE
C    A SEMI-DISCRETE APPROXIMATION TO THE ORIGINAL PDE.
C
C    THE METHOD ON WHICH BDMG IS BASED [2], CONSISTS OF THE APPLICATION
C    OF THE SECOND ORDER BACKWARD DIFFERENTIATION FORMULA TO THIS SEMI-
C    DISCRETE SYSTEM. THE RESULTING NON-LINEAR PROBLEM IS SOLVED USING A
C    MULTIGRID TECHNIQUE.
C
C    THE MAIN CHARACTERISTIC OF THE CODE BDMG IS ITS MINIMAL STORAGE
C    REQUIREMENTS.
C
C    THE APPLICABILITY OF BDMG IS RESTRICTED TO
C    (1)  TWO SPACE DIMENSIONS
C    (2)  RECTANGULAR DOMAINS
C    (3)  SCALAR PARTIAL DIFFERENTIAL EQUATIONS
C    (4)  PDE'S IN WHICH DIFFUSION GREATLY DOMINATES CONVECTION
C    (5)  PDE'S WITHOUT MIXED DERIVATIVES(THIS RESTRICTION ONLY APPLIES
C         IF PDETWO IS USED)
C
C-----------------------------------------------------------------------
C
C   PROBLEM DEFINITION
C
C-----------------------------------------------------------------------
C
C       BDMG APPLIES TO SEMI-DISCRETE SYSTEMS OF THE FORM
C
C
C             DV
C    (1)    [ -- = G (T, V) ]             ,I=1,...,NX
C             DT             I,J           J=1,...,NY   ,
C
C
C    WHERE THE (SCALAR) MATRIX ELEMENTS V(I,J) ARE ASSOCIATED TO THE
C    GRID POINTS (X(I),Y(J)) WHICH FORM A GRID WITH RECTANGULAR MESHES
C    IN THE (X,Y)-PLANE, INCLUDING THE BOUNDARY POINTS.
C    THIS SYSTEM IS SUPPOSED TO ORIGINATE FROM PARTIAL DIFFERENTIAL
C    EQUATIONS OF THE FORM
C
C
C    (2)     UT = F (T, X, Y, U, UX, UY, UXX, UYY, UXY) ,
C
C    WHERE UT     IS THE FIRST PARTIAL DERIVATIVE OF U WITH RESPECT TO T
C    AND   UX,UXX ARE THE FIRST AND SECOND PARTIAL DERIVATIVE OF U WITH
C                 RESPECT TO X
C          ETC.
C    THUS, ANY PDE OF THE FORM (2) DEFINED ON A RECTANGLE WITH INITIAL
C    CONDITION DEFINED ON THIS RECTANGLE AND BOUNDARY CONDITIONS OF
C    THE FORM
C
C    (3)    A(T,X,Y)*U + B(T,X,Y)*UN = C(T,X,Y),  UN IS DERIVATIVE OF U
C                                                 NORMAL TO THE BOUNDARY
C
C    CAN BE DEALT WITH BY BDMG AS SOON AS (2) AND (1) ARE SEMI-
C    DISCRETIZED ON THE GRID (X(I),Y(J)).
C
C    THE INITIAL CONDITION BEING DEFINED FOR EACH V(I,J) NEEDS NOT BE
C    CONSISTENT WITH THE BOUNDARY CONDITIONS.
C
C-----------------------------------------------------------------------
C
C   USER-REQUIRED INFORMATION
C
C-----------------------------------------------------------------------
C
C       IN DEFINING HIS PROBLEM, THE USER IS ASKED TO PROVIDE A MAIN
C    PROGRAM IN WHICH THE GRID IS DEFINED, THE INTEGRATION INTERVAL IS
C    SPECIFIED AND THE INITIAL CONDITION IS SET. BECAUSE BDMG IS A
C    TIME-INTEGRATOR FOR A SYSTEM OF ODE'S IT REQUIRES IN ADDITION A
C    SEMI-DISCRETE APPROXIMATION TO THE ORIGINAL PDE. THIS SEMI-
C    DISCRETIZATION CAN BE PERFORMED EITHER BY HAND OR BY PDETWO;
C    BY MEANS OF THE PARAMETER METH THE USER CAN INDICATE HIS CHOICE.
C    SELECTING THE FIRST POSSIBILITY, A SUBROUTINE IS REQUIRED
C    DELIVERING THE DERIVATIVE FUNCTION G IN (1). IN CASE OF THE USAGE
C    OF PDETWO IT SUFFICES TO SUPPLY ROUTINES DEFINING THE PARTIAL
C    DIFFERENTIAL EQUATION AND THE BOUNDARY CONDITIONS (SEE [3] FOR A
C    DETAILED DESCRIPTION OF THE USER-REQUIRED INFORMATION).
C
C-----------------------------------------------------------------------
C
C   THE PARAMETERS
C
C-----------------------------------------------------------------------
C
C       WE DISTINGUISH BETWEEN THREE TYPES OF PARAMETERS:
C
C    GRID-DEFINING PARAMETERS
C    ------------------------
C
C    NX,NY           ARE THE NUMBER OF MESH LINES IN X- AND Y-DIRECTION,
C                    RESPECTIVELY
C    M               IS THE NUMBER OF SUCCESSIVE GRIDS, USED IN THE
C                    MULTIGRID METHOD(FOR A DISCUSSION OF M SEE BELOW).
C    X(1),...,X(NX)  CONTAIN THE X-COORDINATES OF THE MESH POINTS
C                    ALONG EACH GRID LINE IN THE X-DIRECTION.
C                    (X(1) .LT. X(2) .LT. ... .LT. X(NX)).
C    Y(1),...,Y(NY)  CONTAIN THE ANALOGOUS INFORMATION IN THE
C                    Y-DIRECTION.
C                    (Y(1) .LT. Y(2) .LT. ... .LT. Y(NY)).
C
C    NX,NY,M,X,Y ARE INPUT PARAMETERS AND HAVE TO BE INITIALIZED BY
C    THE USER.
C    NOTE THAT X AND Y ARE NOT NECESSARILY EQUIDISTANT.
C
C       IT IS WORTH GOING INTO MORE DETAIL ABOUT THESE PARAMETERS.
C    BECAUSE THE METHOD IS BASED ON A MULTIGRID METHOD, A SEQUENCE OF M
C    SUCCESSIVE GRIDS IS NECESSARY; THESE GRIDS ARE GENERATED BY THE
C    CODE. THE LEVEL INDEX K RUNS FROM 1 (WHICH IS ASSOCIATED WITH THE
C    COARSEST GRID) UNTIL M (THE FINEST GRID, WHICH IS IDENTICAL TO THE
C    GRID CHOOSEN BY THE USER AND DEFINED BY MEANS OF THE PARAMETERS
C    NX,NY,X,Y). INTERNAL VALUES NXK, NYK ARE CALCULATED AS WELL AS
C    ARRAYS XK(I),I=1,...,NXK AND YK(J),J=1,...,NYK, DEFINING THE GRID
C    AT LEVEL K. THESE VALUES ARE USED, AMONG OTHER THINGS, TO PASS TO
C    THE SUBROUTINE G (SEE BELOW) IN WHICH THE USER HAS TO DEFINE THE
C    DERIVATIVE AT THAT PARTICULAR GRID (XK(I),YK(J)) (ONLY IF METH=1,
C    SEE BELOW).
C
C       A GRID AT LEVEL K IS OBTAINED FROM A GRID AT LEVEL K+1 (DENOTED
C    BY KP1) BY ALTERNATINGLY DELETING GRID LINES;HENCE, NXK=(NXKP1+1)/2
C    XK(1)=XKP1(1), XK(2)=XKP1(3), ETC. THIS MEANS THAT THE GRID LINES
C    ON LEVEL K COINCIDE WITH THOSE ON LEVEL K+1. BECAUSE THE VALUES OF
C    THE ARRAYS XK (K<M) ARE STORED IN THE X-ARRAY TOO, THE USER
C    SHOULD DIMENSION THIS ARRAY AS X(M+(NX-1)*(2**M-1)/(2**(M-1))).
C    SIMILARLY, THE Y-ARRAY MUST BE DECLARED OF LENGTH
C    M+(NY-1)*(2**M-1)/(2**(M-1)).
C
C       THE VALUES OF NX,NY AND M SHOULD BE CHOOSEN IN A SUITABLE
C    RELATIONSHIP TO EACH OTHER. THE NUMBER OF GRIDS SHOULD NOT BE
C    SPECIFIED SO LARGE THAT ONE CONFLICTS WITH ONE OF THE FOLLOWING
C    SITUATIONS:
C    (I)  THE ABOVE ALGORITHM OF DELETING GRID LINES DELIVERS A
C         NON-INTEGER NXK- OR NYK-VALUE FOR SOME K<M.
C    (II) THE GRID ON LEVEL 1 HAS LESS THAN FOUR GRID LINES IN EACH
C         DIRECTION, INCLUDING THE GRID LINES FORMING THE BOUNDARIES.
C
C    IF NONE OF THESE RESTRICTIONS IS APPLICABLE, THE USER IS ADVISED TO
C    CHOOSE M AS LARGE AS POSSIBLE, BECAUSE THE EFFICIENCY OF THE CODE
C    USUALLY INCREASES AS MORE GRIDS ARE AVAILABLE.
C
C    FINALLY, WE MENTION THAT SPECIFYING M=1 IS PROHIBITED.
C
C
C    TIME-INTEGRATION PARAMETERS
C    ---------------------------
C
C    T      IS THE CURRENT TIME; ON ENTRY, T CONTAINS THE INITIAL VALUE
C           OF THE INDEPENDENT VARIABLE IN (1). ON EXIT, T CONTAINS THE
C           VALUE TEND, UNLESS AN ERROR HAS OCCURRED (SEE THE PARAMETER
C           IFLAG BELOW).
C    TEND   SPECIFIES THE POINT AT WHICH THE SOLUTION IS DESIRED
C           (TEND MUST BE .GT. T ON ENTRY).
C    V      IS AN ARRAY TO BE DIMENSIONED IN THE CALLING PROGRAM AS
C           V(NX,NY); IT CONTAINS THE CURRENT SOLUTION VALUES FOR ALL
C           GRID POINTS (SEE ALSO THE REMARKS ON BOUNDARY POINTS IN THE
C           DESCRIPTION OF THE PARAMETER G). ON ENTRY, V SHOULD BE GIVEN
C           THE INITIAL VALUE OF THE DEPENDENT VARIABLE IN (1).
C           ON EXIT, V CONTAINS THE SOLUTION AT THE POINT T.
C    G      IS THE USER-SUPPLIED SUBROUTINE DEFINING THE RIGHT-HAND SIDE
C           FUNCTION G AT LEVEL K. THIS ROUTINE IS ONLY NECESSARY IN
C           CASE OF SEMI-DISCRETIZATION BY HAND (SEE THE PARAMETER METH
C           BELOW). ITS SPECIFICATION READS
C              G(K,NXK,NYK,XK,YK,T,VK,DVK)
C              DIMENSION XK(NXK),YK(NYK),VK(NXK,NYK),DVK(NXK,NYK)
C           THE GRID IS DETERMINED BY XK(I),I=1,...NXK AND YK(J),
C           J=1,...,NYK, CONTAINING THE COORDINATES OF THE GRID LINES AT
C           LEVEL K. GIVEN THE APPROXIMATIONS AT THE CURRENT TIME T IN
C           THE TWO-DIMENSIONAL ARRAY VK, THIS ROUTINE SHOULD DELIVER
C           THE DERIVATIVES AT ALL GRID POINTS IN THE TWO-DIMENSIONAL
C           ARRAY DVK, INCLUDING THE DERIVATIVES AT THE BOUNDARY POINTS.
C           G MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM.
C
C           IN ORDER TO PRESERVE THE STRUCTURE OF THE SYSTEM OF ODE'S,
C           IRRESPECTIVE OF THE TYPE OF BOUNDARY CONDITIONS, THE
C           BOUNDARY POINTS ARE INCLUDED IN THE GRID AND THUS IN THE
C           ARRAYS VK AND DVK.
C           HOWEVER, IN CASE OF DIRICHLET BOUNDARY CONDITIONS TIME-
C           INTEGRATION AT THESE POINTS SEEMS TO BE SUPERFLUOUS;
C           THEREFORE, THE DERIVATIVES AT THESE POINTS MAY BE GIVEN A
C           DUMMY VALUE. IF THIS POSSIBILITY IS USED, ONE SHOULD SPECIFY
C           THE VALUE ZERO FOR THESE DERIVATIVES.
C           AS A CONSEQUENCE OF THIS APPROACH, THE ELEMENTS IN V
C           CORRESPONDING TO BOUNDARY POINTS WHERE A DIRICHLET CONDITION
C           IS PRESCRIBED DO NOT CONTAIN THE APPROXIMATE SOLUTION AT
C           TIME T (EXCEPT WHEN THE DIRICHLET CONDITION IS A CONSTANT,
C           IN WHICH CASE THE VALUE ZERO EQUALS THE DERIVATIVE). HENCE,
C           IN CALCULATING THE DERIVATIVE IN A POINT ADJACENT TO A
C           "DIRICHLET BOUNDARY POINT" ONE SHOULD USE THE DIRICHLET
C           BOUNDARY CONDITION RATHER THAN THE VALUE OF VK.
C
C           IN CASE OF USING PDETWO TO PERFORM THE SEMI-DISCRETIZATION,
C           THE SAME APPROACH IS FOLLOWED IN "DIRICHLET BOUNDARY POINTS"
C    SPRAD  IS A FUNCTION THE USER MUST SUPPLY IN CASE OF THE USAGE OF
C           THE NON-AUTOMATIC OPTION CONCERNING THE SPECTRAL RADIUS (SEE
C           THE PARAMETER INFO BELOW).
C           SPRAD (K, NXK, NYK, XK, YK, T, VK, TAU) SHOULD DELIVER AN
C           ESTIMATE OF THE SPECTRAL RADIUS OF THE JACOBIAN MATRIX OF
C           THE RIGHT-HAND SIDE FUNCTION G AT LEVEL K. THE MEANING OF
C           THE PARAMETERS OF SPRAD IS THE SAME AS IN THE DESCRIPTION OF
C           G, TAU BEING THE CURRENT INTEGRATION STEP.
C           SPRAD MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM.
C    TOL    SPECIFIES THE LOCAL ERROR TOLERANCE PARAMETER
C           (TOL MUST BE .GT. 0).
C    METH   MUST BE GIVEN THE VALUE 1 OR 2.
C           IF METH=1 THE SEMI-DISCRETIZATION HAS TO BE PERFORMED BY
C                     HAND AND THE SUBROUTINE G DELIVERS THE DERIVATIVES
C                     REQUIRED BY THE INTEGRATOR.
C              METH=2 MEANS THE SEMI-DISCRETIZATION HAS TO BE PERFORMED
C                     BY PDETWO.
C
C
C    AUXILIARY PARAMETERS
C    --------------------
C
C
C    WORK   IS A REAL WORK ARRAY FOR INTERNAL STORAGE. IT MUST BE
C           DIMENSIONED AS WORK(**) IN THE CALLING PROGRAM, WHERE
C           ** .GE. (NX-1) * (NY-1) * (71 - 5/(4**(M-2)))/12 +
C                   (NX+NY-2) * (19 - 5/(2**(M-2)))/2 +
C                   7*M + 3 +
C                   (7*NX+16) * (METH-1)
C
C           THE USER SHOULD BE AWARE THAT THIS EXPRESSION MUST
C           DELIVER AN INTEGER VALUE. IF NOT, ONE SHOULD CAREFULLY
C           RECONSIDER THE VALUES OF NX,NY AND M.
C    IWORK  IS AN INTEGER WORK ARRAY CONTAINING POINTERS USED TO
C           PROVIDE DYNAMIC DIMENSIONING IN BDMG. IT MUST BE DIMENSIONED
C           AS IWORK(***) IN THE CALLING PROGRAM, WHERE
C           *** .GE. 10*M+3.
C           IWORK(10*M+3) MUST BE INITIALIZED WITH THE VALUE OF
C                         ** (SEE THE PARAMETER WORK).
C    INFO   IS AN INTEGER ARRAY CONTAINING INFORMATION ABOUT THE STATUS
C           OF THE INTEGRATION PROCESS; IT MUST BE DIMENSIONED AS
C           INFO(****), WHERE
C           **** .GE. 2*M+10
C
C           INFO(1),...,INFO(3)  ARE INPUT PARAMETERS AND MUST BE
C                                INITIALIZED IN THE CALLING PROGRAM
C           INFO(7),...,INFO(2*M+10)  ARE OUTPUT PARAMETERS.
C
C           THE ELEMENTS OF INFO HAVE THE FOLLOWING MEANINGS:
C           INFO(1)=0 TO INDICATE THAT THIS IS THE FIRST CALL TO BDMG.
C                     ON RETURN, BDMG HAS ASSIGNED INFO(1) THE APPROPRI-
C                     ATE VALUE WITH RESPECT TO SUBSEQUENT CALLS.
C           INFO(2)   MAXIMUM NUMBER OF EVALUATIONS OF THE RIGHT-HAND
C                     SIDE OF (1) TO BE SPENT DURING THE INTEGRATION
C                     PROCESS. THIS NUMBER IS THE CONVERTED EQUIVALENT
C                     OF EVALUATIONS ON THE FINEST (I.E. USER-DEFINED)
C                     GRID. TO GET AN IMPRESSION OF THE COSTS OF THE
C                     ALGORITHM WE REFER TO TABLE [2], 3.1 OF SECTION 3.
C           INFO(3)   CAN BE GIVEN THE VALUE 1,2 OR 3. IT IS USED TO
C                     SELECT THE OPTION WITH RESPECT TO THE SPECTRAL
C                     RADIUS OF THE JACOBIAN MATRIX OF THE RIGHT-HAND
C                     SIDE FUNCTION G.
C                     THE THREE OPTIONS AVAILABLE ARE:
C
C               INFO(3)=1 : AUTOMATIC OPTION.
C                           ----------------
C                     WITH THIS OPTION THE CODE ESTIMATES AND CONTROLS
C                     THE SPECTRAL RADIUS DURING THE ENTIRE INTEGRATION
C                     INTERVAL. A FIVE DIAGONAL STRUCTURE OF THE
C                     JACOBIAN MATRIX IS ASSUMED (THIS RESTRICTION
C                     IS ALSO MADE IN THE SEMI-DISCRETIZATION CODE
C                     PDETWO).
C
C               INFO(3)=2 : SEMI-AUTOMATIC OPTION.
C                           ---------------------
C                     IF THE PROBLEM HAS A CONSTANT SPECTRAL RADIUS
C                     (E.G. IN LINEAR PROBLEMS) OR IF IT IS AT LEAST
C                     NON-INCREASING WITH TIME, ONE MAY SELECT THE SEMI-
C                     AUTOMATIC OPTION WHICH ONLY CALCULATES THE
C                     SPECTRAL RADIUS AT THE BEGINNING OF THE INTEGRAT-
C                     ION PROCESS. NO CONTROL IS PERFORMED. IF THIS
C                     SITUATION APPLIES THIS OPTION MAY SAVE A LOT OF
C                     RIGHT-HAND SIDE EVALUATIONS. ONE SHOULD ONLY APPLY
C                     THIS OPTION WHEN THE RIGHT-HAND SIDE FUNCTION G
C                     HAS A FIVE-POINT COUPLING.
C
C               INFO(3)=3 : NON-AUTOMATIC OPTION.
C                           --------------------
C                     SOME PROBLEMS ALLOW THE USER THE CALCULATION
C                     BY HAND OF THE SPECTRAL RADIUS, SO THAT AN
C                     EXPLICIT EXPRESSION IN TERMS OF T, THE MESH SIZES
C                     AT LEVEL K, THE CURRENT SOLUTION V AND THE CUR-
C                     RENT TIME STEP TAU CAN BE GIVEN TO BE USED IN THE
C                     INTERVAL [T,T+TAU]. AGAIN, MANY RIGHT-HAND SIDE
C                     EVALUATIONS CAN BE AVOIDED AND, IF NO USE IS MADE
C                     OF PDETWO, THE RESTRICTION TO FIVE-POINT COUPLING
C                     NO LONGER APPLIES.
C
C                     IT SHOULD BE NOTED THAT AN UNDER ESTIMATE OF THE
C                     SPECTRAL RADIUS MAY GIVE RISE TO INTERNAL
C                     INSTABILITIES IN THE ITERATION PROCESS WITHIN
C                     ONE INTEGRATION STEP. BECAUSE THESE INSTABILITIES
C                     CAN QUICKLY GIVE RISE TO VERY LARGE NUMBERS, IT
C                     MIGHT HAPPEN THAT AN OVERFLOW WOULD OCCUR IN
C                     THE ITERATION PROCESS BEFORE CONTROL IS RETURNED
C                     TO THE DRIVER, WHERE THE LOCAL ERROR ESTIMATOR
C                     SURELY WOULD HAVE DETECTED THIS INSTABILITY.
C                     THEREFORE, WE HAVE PROTECTED THE CODE AGAINST
C                     OVERFLOW.
C                     WHATEVER OPTION IS CHOSEN, THE OVERFLOW CHECKING
C                     IS ALWAYS PERFORMED. IN CASE OF OPTIONS 2 AND 3,
C                     THE INTEGRATION PROCESS IS DISCONTINUED WHEN
C                     VIOLATING THE OVERFLOW CONDITION; CONTROL IS
C                     RETURNED TO THE CALLING PROGRAM AND THE USER CAN
C                     TAKE ACTION WITH RESPECT TO THE SPECTRAL RADIUS.
C
C           INFO(4),...,INFO(6)   ARE USED FOR INTERNAL CONTROL
C
C           INFO(7)   TOTAL NUMBER OF INTEGRATION STEPS PERFORMED,
C                     INCLUDING REJECTED ONES.
C           INFO(8)   TOTAL NUMBER OF REJECTED STEPS.
C           INFO(9)   NUMBER OF TIMES A STEP HAS BEEN ABANDONED BECAUSE
C                     OF VIOLATING THE OVERFLOW-TEST.
C           INFO(10)  TOTAL NUMBER OF DERIVATIVE-EVALUATIONS, EXPRESSED
C                     IN TERMS OF EVALUATIONS ON THE ORIGINAL GRID.
C           INFO(10+K) NUMBER OF TIMES THE DERIVATIVE IS EVALUATED AT
C                     LEVEL K (K=1,...,M) USED FOR THE INTEGRATION
C                     PROCESS ONLY.
C           INFO(10+M+K)  NUMBER OF TIMES THE DERIVATIVE IS EVALUATED AT
C                     LEVEL K (K=1,...,M) USED FOR THE EVALUATION AND
C                     CONTROL OF THE SPECTRAL RADIUS.
C
C    IFLAG  IS A FLAG USED TO INDICATE VARIOUS ERROR CONDITIONS.
C           ON RETURN, IT MAY HAVE THE FOLLOWING VALUES AND MEANINGS:
C
C           IFLAG=0  SUCCESSFUL INTEGRATION; THAT IS, THE END POINT TEND
C                    HAS BEEN REACHED. TO CONTINUE THE INTEGRATION
C                    IT SUFFICES TO DEFINE A NEW VALUE OF TEND AND
C                    RECALL BDMG.
C           IFLAG=1  THE VALUE OF M IS TOO LARGE IN CONNECTION WITH THE
C                    SPECIFIED VALUES OF NX AND NY.
C           IFLAG=2  X(I) .GE. X(I+1)   OR   Y(J) .GE. Y(J+1)  FOR SOME
C                    I OR J
C           IFLAG=3  THE LENGTH OF THE ARRAY WORK IS LESS THAN THE
C                    REQUIRED MINIMUM.
C           IFLAG=4  TOL .LE. 0  OR  TEND .LE. T
C
C           IN THE LAST FOUR CASES THE PROCESS IS NOT STARTED AND ONE
C           SHOULD RECONSIDER THE INPUT REQUIREMENTS GIVEN AT THE
C           DESCRIPTION OF THE APPROPRIATE PARAMETERS.
C
C           IFLAG=5  THE STEP LENGTH HAS BEEN REDUCED TO AN UNACCEPT-
C                    ABLY SMALL VALUE. THE PROBLEM SEEMS TO BE
C                    UNSOLVABLE TO BDMG. A POSSIBLE REASON MAY BE TOO
C                    STRINGENT ACCURACY REQUIREMENTS (PARAMETER TOL).
C           IFLAG=6  THE MAXIMUM NUMBER OF DERIVATIVE EVALUATIONS
C                    HAS BEEN SPENT; IF THE USER DECIDES TO CONTINUE,
C                    HE ONLY HAS TO INCREASE INFO(2) AND RECALL BDMG.
C           IFLAG=7  THE OVERFLOW-TEST WAS VIOLATED WHILE INFO(3) .NE. 1
C                    THE USER MUST TAKE ACTION WITH RESPECT TO THE
C                    SPECTRAL RADIUS OR CHANGE TO THE AUTOMATIC OPTION.
C
C-----------------------------------------------------------------------
C
C   BRIEF SUMMARY OF INPUT PARAMETERS
C
C-----------------------------------------------------------------------
C
C    NX,NY,M,X,Y   DEFINE THE GRID AND THE NUMBER OF GRIDS
C    T,TEND        DEFINE THE INTEGRATION INTERVAL
C    V             CONTAINS THE INITIAL VALUE
C    TOL           IS THE LOCAL TOLERANCE PARAMETER
C    METH          SELECTS THE OPTION CONCERNING THE SEMI-DISCRETIZATION
C    IWORK(10*M+3) SPECIFIES THE LENGTH OF THE ARRAY WORK.
C    INFO(1)       EQUALS ZERO, ONLY AT FIRST CALL
C    INFO(2)       IS THE MAXIMUM NUMBER OF DERIVATIVE-EVALUATIONS
C    INFO(3)       SELECTS THE OPTION CONCERNING THE SPECTRAL RADIUS
C
C    THE PARAMETERS TEND, TOL, INFO(2), AND INFO(3) MAY BE CHANGED FROM
C    CALL TO CALL.
C
C-----------------------------------------------------------------------
C
C   PROGRAMMING NOTES
C
C-----------------------------------------------------------------------
C
C     -   WE EMPHASIZE THAT, USING OPTION 3 WITH RESPECT TO THE
C         SPECTRAL RADIUS (CF. INFO(3)), THE CODE RELIES ON THE USER-
C         SPECIFIED ESTIMATE. THERE IS A CALL TO THE SUBROUTINE SPRAD
C         AT THE START OF EACH INTEGRATION STEP AND THE USER MUST BE
C         SURE TO DELIVER AN ESTIMATE WHICH HOLDS ON THE WHOLE
C         INTERVAL [T,T+TAU].
C     -   CONCERNING THE TOTAL AMOUNT OF STORAGE REQUIRED BY BDMG, WE
C         NOTE THAT THE CODE NEEDS AT MOST STORAGE EQUIVALENT WITH SEVEN
C         ARRAYS, THE SIZE OF WHICH CORRESPONDS WITH THE NUMBER OF
C         GRID POINTS IN THE ORIGINAL GRID.
C     -   IN SOLVING PARABOLIC EQUATIONS TWO TYPES OF ERRORS ARISE,VIZ.
C         SPACE-DISCRETIZATION ERRORS AND TIME-INTEGRATION ERRORS. WE
C         MENTION THAT ONLY THE LAST ONES ARE CONTROLLED BY BDMG (BY
C         MEANS OF THE PARAMETER TOL). HENCE, IN CHOOSING A VALUE OF
C         TOL THE USER IS ADVISED TO CONSIDER THE ACCURACY OF THE SPACE-
C         DISCRETIZATION.
C     -   THE MINIMUM STEP SIZE DEPENDS ON THE MACHINE ROUNDOFF U AND
C         THE UNDERFLOW NUMBER P, WHICH OBVIOUSLY ARE MACHINE-DEPENDENT.
C         ADAPT, IF NECESSARY, THE DATA STATEMENT BELOW, IN WHICH
C         FOURU=4*U AND TENP=10*P.
C     -   THE PROGRAM HAS BEEN WRITTEN IN FORTRAN IV AND MEETS,WITH THE
C         EXCEPTION OF THE RULES CONCERNING THE SUBSCRIPTING OF ARRAYS,
C         THE REQUIREMENTS OF THE PFORT VERIFYER [4]. MOREOVER, THE
C         PROGRAM CONFORMS THE ANSI X3.9-1978 FORTRAN 77 STANDARD [1].
C
C--------------------------------------------------------------------
C
C   REFERENCES
C
C----------------------------------------------------------------------
C
C   [1]   AMERICAN NATIONAL STANDARDS INSTITUTE. AMERICAN NATIONAL
C              STANDARD PROGRAMMING LANGUAGE FORTRAN, ANSI X3.9-1978,
C              AMERICAN NATIONAL STANDARDS INSTITUTE, NEW YORK, 1978.
C
C   [2]   HOUWEN,P.J. VAN DER AND B.P.SOMMEIJER, ANALYSIS OF CHEBYSHEV
C              RELAXATION IN MULTIGRID METHODS FOR NONLINEAR PARABOLIC
C              DIFFERENTIAL EQUATIONS, ZAMM 63, 1983, PP. 193-201.
C
C   [3]   MELGAARD,D.K. AND R.F.SINCOVEC, GENERAL SOFTWARE FOR TWO-
C              DIMENSIONAL NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS,
C              ACM TRANS. MATH. SOFTW., VOL 7, MARCH 1981, PP. 106-125.
C
C   [4]   RYDER,B.F. AND A.D.HALL, THE PFORT VERIFIER. BELL LABORATORIES
C              COMPUTER SCIENCE REP. 2 (MAY 1973-JANUARY 1981).
C
C-----------------------------------------------------------------------
C
C
C
      DIMENSION X(1),Y(1),V(NX,NY),WORK(1),IWORK(1),INFO(1)
      LOGICAL FIRST,LAST
      COMMON/COEF/A1,A2,B0
      COMMON/MACH/ZETA
      EXTERNAL G,SPRAD
C-----------------------------------------------------------------------
C   SET MACHINE-DEPENDENT CONSTANTS
C       FOURU = 4 * U,
C        TENP =10 * P,
C   WHERE U IS THE MACHINE ROUNDOFF AND P IS THE UNDERFLOW NUMBER
C-----------------------------------------------------------------------
      DATA FOURU/2.842E-14/,TENP/3.132E-293/
      ZETA=10.0**(AMAX1(10.0,-ALOG10(TENP)/10.0))
C-----------------------------------------------------------------------
C   CHECK IF PROBLEM PARAMETERS ARE CORRECT
C-----------------------------------------------------------------------
      IF (TOL .GT. 0.0 .AND. TEND .GT. T) GO TO 10
      IFLAG = 4
      RETURN
   10 M2 = M*2
      M3 = M*3
      M8 = M*8
      M210 = M2 + 10
      IFLAG = 0
C-----------------------------------------------------------------------
C   ON FIRST CALL, PRELIM IS CALLED TO CHECK INPUT PARAMETERS AND TO
C   PROVIDE DYNAMIC DIMENSIONING FOR THE ARRAYS USED IN THE INTEGRATION
C   PROCESS
C-----------------------------------------------------------------------
      IF (INFO(1) .EQ. 1) GO TO 50
      CALL PRELIM (NX, NY, M, X, Y, METH, IWORK, IFLAG)
      IF (IFLAG .NE. 0) RETURN
C-----------------------------------------------------------------------
C   IF PDETWO IS USED, CALL STRSET TO PROVIDE DYNAMIC DIMENSIONING FOR
C   THE ARRAYS USED IN PDETWO
C-----------------------------------------------------------------------
      IF (METH .EQ. 2) CALL STRSET (NX, NY, M, IWORK)
      JSE = IWORK(M8+1)
      JSU = IWORK(M8+2)
      N1 = IWORK(M2+1)
      JN = IWORK(M3+3)
      JNM1 = IWORK(M3+4)
      DO 20 L = 1,N1
         WORK(JN+L) = V(L,1)
         WORK(JNM1+L) = V(L,1)
   20    CONTINUE
      DO 30 I = 4,M210
         INFO(I) = 0
   30    CONTINUE
C-----------------------------------------------------------------------
C   CALCULATE AN INITIAL STEP SIZE AND ADAPT IT TO THE LENGTH OF THE
C   REMAINING INTEGRATION INTERVAL. THIS ADAPTATION IS PERFORMED TO
C   ARRIVE EXACTLY AT THE POINT TEND. HENCE, INTERPOLATING THE NUMERICAL
C   SOLUTION AT TEND CAN BE AVOIDED
C-----------------------------------------------------------------------
      TAUNEW = TAUST (NX, NY, M, X, Y, T, V, G, TOL, METH, WORK, IWORK,
     *                INFO)
      CALL FIT (T, TEND, TAUNEW, LAST)
      IF (TAUNEW .GT. FOURU*ABS(T) + TENP) GO TO 40
      IFLAG = 5
      RETURN
C-----------------------------------------------------------------------
C   SET COEFFICIENTS FOR THE ONE-STEP BACKWARD DIFFERENTIATION FORMULA
C-----------------------------------------------------------------------
   40 A1 = 1.0
      A2 = 0.0
      B0 = 1.0
      QNEW = 1.0
      FIRST = .TRUE.
      GO TO 70
C-----------------------------------------------------------------------
C   IN CASE OF A SUBSEQUENT CALL, RESET DATA FOR CONTINUING THE PROCESS
C-----------------------------------------------------------------------
   50 TAUNEW = WORK(IWORK(M8+2)+M)
      TAUOLD = WORK(IWORK(M8+2)+M+1)
      QOLD = WORK(IWORK(M8+2)+M+2)
      JNM1 = IWORK(M3+4)
      JSE = IWORK(M8+1)
      JSU = IWORK(M8+2)
      FIRST = .FALSE.
      CALL FIT (T, TEND, TAUNEW, LAST)
C-----------------------------------------------------------------------
C   SET COEFFICIENTS FOR THE TWO-STEP BACKWARD DIFFERENTIATION FORMULA
C-----------------------------------------------------------------------
   60 QNEW = TAUOLD/TAUNEW
      QAUX = QNEW*(QNEW + 2.0)
      A1 = (QAUX + 1.0)/QAUX
      A2 = 1.0 - A1
      B0 = (QAUX - QNEW)/QAUX
      IF (INFO(3).EQ.2 .OR. ERROR.GT.EPS .OR. INFO(4).EQ.1) GO TO 80
C-----------------------------------------------------------------------
C   CALCULATE THE VALUE OF THE SPECTRAL RADIUS AT EACH GRID
C-----------------------------------------------------------------------
   70 CALL SPRI (NX, NY, M, X, Y, T, TAUNEW, V, G, SPRAD, WORK(JSE),
     *           WORK(JSU), METH, WORK, IWORK, INFO, QNEW)
C-----------------------------------------------------------------------
C   PERFORM ONE INTEGRATION STEP
C-----------------------------------------------------------------------
   80 CALL STEP (NX, NY, M, X, Y, T, TAUNEW, V, G, WORK(JSU), METH,
     *           WORK, IWORK, INFO, QNEW)
C-----------------------------------------------------------------------
C   CHECK ON OVERFLOW CONDITION; VERIFY WHETHER THE VALUES OF THE
C   SPECTRAL RADII WERE UPDATED AT THE BEGINNING OF THE PARTICULAR STEP.
C   IF NOT, WE DO SO AND TRY AGAIN, OTHERWISE WE HALVE THE STEP SIZE.
C   IF, IN THE SAME INTERVAL, THE OVERFLOW CONDITION IS VIOLATED AGAIN,
C   WE SET TAU:=TAU/5, REPEATEDLY IF NECESSARY
C-----------------------------------------------------------------------
      IF (INFO(6) .EQ. 0) GO TO 130
      INFO(9) = INFO(9) + 1
      IF (INFO(3) .EQ. 1 .OR. (INFO(3) .EQ. 2 .AND. FIRST) ) GO TO 90
      IFLAG = 7
      GO TO 180
   90 IF (INFO(4) .EQ. 1) GO TO 100
      INFO(4) = 1
      IF (INFO(5) .EQ. 0) GO TO 70
      TAUNEW = TAUNEW/2.0
      GO TO 110
  100 TAUNEW = TAUNEW/5.0
  110 IF (TAUNEW .GT. FOURU*ABS(T) + TENP) GO TO 120
      IFLAG = 5
      GO TO 180
  120 LAST = .FALSE.
      INFO(6) = 0
      IF(FIRST) GO TO 80
      GO TO 60
C-----------------------------------------------------------------------
C   THE STEP HAS BEEN FINISHED WITHOUT DETECTION OF INSTABILITIES.
C   ESTIMATE THE LOCAL TRUNCATION ERROR
C-----------------------------------------------------------------------
  130 INFO(7) = INFO(7) + 1
      CALL LTE (NX, NY, M, X, Y, T, TAUNEW, V, G, TOL, METH, WORK,
     *          IWORK, INFO, QOLD, QNEW, EPS, ERROR)
      IF (ERROR. LT. EPS) GO TO 150
C-----------------------------------------------------------------------
C   THE ERROR TEST FAILED. CALULATE A NEW STEP SIZE, ADAPTED TO THE
C   LENGTH OF THE INTEGRATION INTERVAL
C-----------------------------------------------------------------------
      INFO(8) = INFO(8) + 1
      CALL NEWSTP (EPS, ERROR, TAUNEW, TAUNEW, FIRST)
      CALL FIT (T, TEND, TAUNEW, LAST)
      IF (TAUNEW .GT. FOURU*ABS(T) + TENP) GO TO 140
      IFLAG = 5
      GO TO 180
  140 IF (FIRST) GO TO 80
      GO TO 60
C-----------------------------------------------------------------------
C   THE STEP HAS BEEN SUCCESSFUL. SAVE, IF AVAILABLE, THE SOLUTION
C   AT THE POINT T-TAU TO BE USED IN THE ERROR-ESTIMATION IN THE NEXT
C   STEP
C-----------------------------------------------------------------------
  150 INFO(4) = 0
      IF (.NOT. FIRST) GO TO 160
      FIRST = .FALSE.
      INFO(1) = 1
      GO TO 170
  160 CALL INJECT (NX, NY, IWORK(2), IWORK(M+2), WORK(JNM1+1),
     *             WORK(IWORK(M8)+1))
C-----------------------------------------------------------------------
C   SHIFT THE SOLUTION VECTORS AND THE T-ARGUMENT
C-----------------------------------------------------------------------
  170 CALL SHIFT (NX, NY, M, T, TAUOLD, TAUNEW, V, WORK, IWORK)
      QOLD = QNEW
C-----------------------------------------------------------------------
C   CALCULATE A NEW STEP SIZE. STOP THE INTEGRATION IF TEND IS REACHED;
C   OTHERWISE, ADAPT THE STEP SIZE TO THE LENGTH OF THE INTEGRATION
C   INTERVAL
C-----------------------------------------------------------------------
      CALL NEWSTP (EPS, ERROR, TAUOLD, TAUNEW, FIRST)
      IF (LAST) GO TO 200
      CALL FIT (T, TEND, TAUNEW, LAST)
C-----------------------------------------------------------------------
C   TEST IF THE NUMBER OF G-EVALUATIONS SPENT SO FAR EXCEEDS THE
C   MAXIMALLY ALLOWED NUMBER OF G-EVALUATIONS
C-----------------------------------------------------------------------
      INFO(10) = INT(SUMG(M,INFO) + 0.5)
      IF (INFO(2) .GT. INFO(10)) GO TO 60
      IFLAG = 6
      GO TO 200
C-----------------------------------------------------------------------
C   CONTROL REACHES THIS POINT WHEN EITHER THE STEP SIZE BECOMES TOO
C   SMALL OR INSTABILITIES ARISED AND THE AUTOMATIC OPTION CONCERNING
C   THE SPECTRAL RADIUS WAS NOT SELECTED. IN BOTH CASES, THE INTEGRATION
C   PROCESS IS STOPPED
C-----------------------------------------------------------------------
  180 DO 190 L = 1,N1
         V(L,1) = WORK(JN+L)
  190    CONTINUE
C-----------------------------------------------------------------------
C   SAVE DATA FOR SUBSEQUENT CALLS
C-----------------------------------------------------------------------
  200 WORK(IWORK(M8+2)+M) = TAUNEW
      WORK(IWORK(M8+2)+M+1) = TAUOLD
      WORK(IWORK(M8+2)+M+2) = QOLD
      INFO(10) = INT(SUMG(M,INFO) + 0.5)
      RETURN
      END
      SUBROUTINE PRELIM (NX, NY, M, X, Y, METH, IWORK, IFLAG)
      DIMENSION X(1),Y(1),IWORK(1)
C-----------------------------------------------------------------------
C   PRELIM CHECKS THE INPUT PARAMETERS CONCERNING THE DEFINITION OF THE
C   GRID, DETERMINES THE COORDINATES OF THE ADDITIONAL GRIDS AND SETS
C   THE POINTERS IN THE ARRAY IWORK TO PROVIDE DYNAMIC DIMENSIONING FOR
C   THE ARRAYS IN THE INTEGRATION PROCESS. ACTUALLY, MOST OF THE
C   ELEMENTS OF IWORK ARE POINTERS TO THE ENTRY POINTS OF THE AUXILIARY
C   VECTORS, WHICH ARE COLLECTED IN THE ARRAY WORK, THE USER-DIMENSIONED
C   WORKING SPACE
C-----------------------------------------------------------------------
      WKL = FLOAT((NX-1)*(NY-1))*(71.0 - 5.0/4.0**(M-2))/12.0 +
     *      FLOAT(NX+NY-2)*(19.0 - 5.0/2.0**(M-2))/2.0 + 7.0*FLOAT(M) +
     *      3.0
      IF (METH .EQ. 2) WKL = WKL + 7.0*FLOAT(NX) + 16.0
      IF (IWORK(10*M+3) .GE. INT(WKL + 0.5)) GO TO 10
      IFLAG = 3
      RETURN
   10 M2 = M*2
      M3 = M*3
      M8 = M*8
      M9 = M*9
C-----------------------------------------------------------------------
C   CHECK THE GRID-DEFINING PARAMETERS
C-----------------------------------------------------------------------
      IF ((NX/2)*2 .EQ. NX .OR. (NY/2)*2 .EQ. NY) GO TO 50
      DO 20 I = 2,NX
         IF (X(I) .GT. X(I-1)) GO TO 20
         IFLAG = 2
         GO TO 120
   20    CONTINUE
      DO 30 I = 2,NY
         IF (Y(I) .GT. Y(I-1)) GO TO 30
         IFLAG = 2
         GO TO 120
   30    CONTINUE
C-----------------------------------------------------------------------
C   DEFINITION OF THE ADDITIONAL GRIDS
C-----------------------------------------------------------------------
      IWORK(1) = NX
      IWORK(M+1) = NY
      IWORK(M2+1) = NX*NY
      DO 40 I = 2,M
         NXI = (IWORK(I-1) + 1)/2
         IF ((NX/2)*2 .EQ. NX .AND. I .NE. M) GO TO 50
         IWORK(I) = NXI
         NYI = (IWORK(M+I-1) + 1)/2
         IF ((NY/2)*2 .EQ. NY .AND. I .NE. M) GO TO 50
         IWORK(M+I) = NYI
         IWORK(M2+I) = NXI*NYI
   40    CONTINUE
      GO TO 60
   50 IFLAG = 1
      RETURN
C-----------------------------------------------------------------------
C   CALCULATION OF THE ADDITIONAL GRIDS
C-----------------------------------------------------------------------
   60 ISXO = 0
      ISX = NX
      ISYO = 0
      ISY = NY
      IWORK(M8+3) = 1
      IWORK(M9+3) = 1
      DO 90 K = 2,M
         NXK = IWORK(K)
         DO 70 I = 1,NXK
            X(ISX+I) = X(ISXO+2*I-1)
   70       CONTINUE
         ISXO = ISX
         ISX = ISX + IWORK(K)
         IWORK(M8+2+K) = IWORK(M8+1+K) + IWORK(K-1)
         NYK = IWORK(M+K)
         DO 80 I = 1,NYK
            Y(ISY+I) = Y(ISYO+2*I-1)
   80       CONTINUE
         ISYO = ISY
         ISY = ISY + IWORK(M+K)
         IWORK(M9+2+K) = IWORK(M9+1+K) + IWORK(M+K-1)
   90    CONTINUE
C-----------------------------------------------------------------------
C   DEFINITION OF THE POINTERS TO PROVIDE DYNAMIC DIMENSIONING
C-----------------------------------------------------------------------
      DO 100 I = 1,4
         IWORK(M3+I) = (I - 1)*IWORK(M2+1)
  100    CONTINUE
      DO 110 I = 2,M
         NA = M3 + 5*I - 6
         DO 110 J = 1,5
         IWORK(NA+J) = IWORK(NA) + IWORK(M2+I-1) + (J - 1)*IWORK(M2+I)
  110    CONTINUE
      IWORK(M8) = IWORK(M8-1) + IWORK(M3)
      IWORK(M8+1) = IWORK(M8) + IWORK(M2+2) + 1
      IWORK(M8+2) = IWORK(M8+1) + M
  120 CONTINUE
      RETURN
      END
      SUBROUTINE STRSET (NX, NY, M, IWORK)
      DIMENSION IWORK(1)
      COMMON/ISTORE/IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,
     *              IW11,IW12,IW13,IW14,IW15,IW16,IW17,IW18,
     *              IW19,IW20,IW21,IW22
C-----------------------------------------------------------------------
C   STRSET IS DESIGNED TO PROVIDE DYNAMIC DIMENSIONING FOR THE ARRAYS
C   USED IN PDETWO
C-----------------------------------------------------------------------
      IWKL = 4*NX*NY + 2*M + 3 + IWORK(2*M+2)
      DO 10 I = 2,M
         IWKL = IWKL + 5*IWORK(2*M+I)
   10    CONTINUE
      IW1 = IWKL + 1
      IW2 = IW1 + NX
      IW3 = IW2 + NX
      IW4 = IW3 + 1
      IW5 = IW4 + 2
      IW6 = IW5 + 1
      IW7 = IW6 + 1
      IW8 = IW7 + 1
      IW9 = IW8 + 1
      IW10 = IW9 + 1
      IW11 = IW10 + 1
      IW12 = IW11 + 1
      IW13 = IW12 + 1
      IW14 = IW13 + NX
      IW15 = IW14 + NX
      IW16 = IW15 + NX
      IW17 = IW16 + NX
      IW18 = IW17 + NX
      IW19 = IW18 + 1
      IW20 = IW19 + 1
      IW21 = IW20 + 1
      IW22 = IW21 + 1
      RETURN
      END
      SUBROUTINE DERIV (K, NX, NY, X, Y, T, U, DU, G, METH, WORK)
      DIMENSION X(1),Y(1),U(NX,NY),DU(NX,NY),WORK(1)
      COMMON/ISTORE/IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,
     *              IW11,IW12,IW13,IW14,IW15,IW16,IW17,IW18,
     *              IW19,IW20,IW21,IW22
C-----------------------------------------------------------------------
C   DEPENDING ON THE PARAMETER METH, DERIV SELECTS THE DERIVATIVE
C   FUNCTION TO BE CALLED. IF METH = 2, PDETWO IS CALLED TO EVALUATE A
C   SEMI-DISCRETE APPROXIMATION TO THE PARTIAL DIFFERENTIAL EQUATION;
C   IF METH = 1, THE USER-SUPPLIED SUBROUTINE G IS CALLED
C-----------------------------------------------------------------------
      IF (METH .EQ. 2) GO TO 10
      CALL G (K, NX, NY, X, Y, T, U, DU)
      RETURN
   10 CALL PDETWO (1, NX, NY, X, Y, T, U, DU,
     *             WORK(IW1 ), WORK(IW2 ), WORK(IW3 ), WORK(IW4 ),
     *             WORK(IW5 ), WORK(IW6 ), WORK(IW7 ), WORK(IW8 ),
     *             WORK(IW9 ), WORK(IW10), WORK(IW11), WORK(IW12),
     *             WORK(IW13), WORK(IW14), WORK(IW15), WORK(IW16),
     *             WORK(IW17), WORK(IW18), WORK(IW19), WORK(IW20),
     *             WORK(IW21), WORK(IW22))
      RETURN
      END
      SUBROUTINE STEP (NX, NY, M, X, Y, T, TAU, V, G, SUSE, METH,
     *                 WORK, IWORK, INFO, Q)
      DIMENSION X(1),Y(1),V(NX,NY),SUSE(1),WORK(1),IWORK(1),INFO(1)
      COMMON/COEF/A1,A2,B0
      COMMON/OFLOW/VNNORM,FNNORM
      EXTERNAL G
C-----------------------------------------------------------------------
C   STEP ACTUALLY PERFORMES ONE STEP IN THE INTEGRATION OF THE INITIAL-
C   VALUE PROBLEM  DV/DT = G(T,V). GIVEN THE VALUE OF V AT TIME T, STEP
C   CALCULATES AN APPROXIMATION TO THE SOLUTION AT TIME T+TAU. THE
C   METHOD IS BASED ON THE TWO-STEP BACKWARD DIFFERENTIATION FORMULA.
C   THE OCCURING IMPLICIT RELATIONS ARE APPROXIMATELY SOLVED USING A
C   MULTIGRID TECHNIQUE. ALL ADDITIONAL INFORMATION IS STORED IN THE
C   WORKING ARRAY WORK
C-----------------------------------------------------------------------
      M2 = M*2
      M3 = M*3
      M8 = M*8
      M8P2 = M8 + 2
      M9P2 = M8P2 + M
      M10P2 = M9P2 + M
      N1 = NX*NY
      JB = IWORK(M3+1)
      JA1 = IWORK(M3+2)
      JN = IWORK(M3+3)
      JNM1 = IWORK(M3+4)
      DV = 2.222222222**(1.0/FLOAT(M - 1))
      ALFAV = 0.004**(1.0/FLOAT(M - 1))
      NX2 = IWORK(2)
      NY2 = IWORK(M+2)
C-----------------------------------------------------------------------
C   CONSTRUCT AN INITIAL APPROXIMATION ON THE FINEST GRID BY MEANS OF
C   EXTRAPOLATION
C-----------------------------------------------------------------------
      DO 10 L = 1,N1
         WORK(JB+L) = ((Q + 1.0)*WORK(JN+L) - WORK(JNM1+L))/Q
   10    CONTINUE
C-----------------------------------------------------------------------
C   PERFORM TWO SWEEPS; BOTH ARE OF SO-CALLED V-CYCLE TYPE
C-----------------------------------------------------------------------
      DO 160 ITER = 1,2
      CALL DERIV (M, NX, NY, X, Y, T+TAU, WORK(JB+1), WORK(JA1+1),
     *            G, METH, WORK)
      INFO(M+10) = INFO(M+10) + 1
      DO 20 L = 1,N1
         V(L,1) = WORK(JB+L) - B0*TAU*WORK(JA1+L) -
     *            (A1*WORK(JN+L) + A2*WORK(JNM1+L))
   20    CONTINUE
C-----------------------------------------------------------------------
C   CONSTRUCT THE SEQUENCE OF MODIFIED RIGHT-HAND SIDE FUNCTIONS, USED
C   IN THE MULTIGRID METHOD
C-----------------------------------------------------------------------
      DO 40 I = 2,M
         NXI = IWORK(I)
         NYI = IWORK(M+I)
         NI = NXI*NYI
         NXIM1 = IWORK(I-1)
         NYIM1 = IWORK(M+I-1)
         M3P5I = M3 + 5*I
         JI = IWORK(M3P5I-5)
         JIB = IWORK(M3P5I-4)
         JIM1B = IWORK(M3P5I-9)
         JIA1 = IWORK(M3P5I-3)
         IF (I .EQ .2) CALL RESTR (NXIM1, NYIM1, NXI, NYI, X, Y, V,
     *                             WORK(JI+1))
         IF (I .GT. 2) CALL RESTR (NXIM1, NYIM1, NXI, NYI,
     *                           X(IWORK(M8P2+I-1)), Y(IWORK(M9P2+I-1)),
     *                            WORK(IWORK(M3P5I-10)+1), WORK(JI+1))
         CALL RESTR (NXIM1, NYIM1, NXI, NYI, X(IWORK(M8P2+I-1)),
     *               Y(IWORK(M9P2+I-1)), WORK(JIM1B+1), WORK(JIB+1))
         CALL DERIV (M+1-I, NXI, NYI, X(IWORK(M8P2+I)),
     *               Y(IWORK(M9P2+I)), T+TAU, WORK(JIB+1),
     *               WORK(JIA1+1), G, METH, WORK)
         INFO(11+M-I) = INFO(11+M-I) + 1
         DO 30 L = 1,NI
            WORK(IWORK(M3P5I-1)+L) = WORK(JIB+L) - B0*TAU*WORK(JIA1+L) -
     *                               WORK(JI+L)
   30       CONTINUE
         IF (ITER .GT. 1) GO TO 40
         JIA2 = IWORK(M3P5I-2)
         JIM1A2 = IWORK(M3P5I-7)
         CALL RESTR (NXIM1, NYIM1, NXI, NYI, X(IWORK(M8P2+I-1)),
     *               Y(IWORK(M9P2+I-1)), WORK(JIM1A2+1), WORK(JIA2+1))
   40    CONTINUE
      NXM = IWORK(M)
      NYM = IWORK(M2)
      NM = NXM*NYM
      JM = IWORK(M8-5)
      JMB = IWORK(M8-4)
      JMA1 = IWORK(M8-3)
      JMA2 = IWORK(M8-2)
      JMS = IWORK(M8-1)
      IF (ITER .GT. 1) GO TO 60
C-----------------------------------------------------------------------
C   DURING THE FIRST SWEEP, WE SOLVE THE COARSEST-GRID PROBLEM USING A
C   STABILIZED RUNGE-KUTTA METHOD
C-----------------------------------------------------------------------
      DO 50 L = 1,NM
         WORK(JM+L) = WORK(JMA2+L)
   50    CONTINUE
      CALL RKR (NXM, NYM, X(IWORK(M9P2)), Y(IWORK(M10P2)), T, TAU,
     *          WORK(JM+1), WORK(JMB+1), WORK(JMA1+1), WORK(JMA2+1),
     *          WORK(JMS+1), G, SUSE(M), METH, WORK, INFO)
      GO TO 80
C-----------------------------------------------------------------------
C   IN THE SECOND SWEEP, WE PERFORM CHEBYSHEV ITERATION AS A RELAXATION
C   PROCESS FOR THE COARSEST-GRID PROBLEM
C-----------------------------------------------------------------------
   60 DO 70 L = 1,NM
         WORK(JM+L) = WORK(JMB+L)
   70    CONTINUE
      CALL CHEB (NXM, NYM, M, X(IWORK(M9P2)), Y(IWORK(M10P2)), T, TAU,
     *           WORK(JM+1), WORK(JMB+1), WORK(JMA1+1), WORK(JMA2+1),
     *           WORK(JMS+1), G, SUSE(M), METH, WORK, INFO, 10.0, 0.001,
     *           M)
   80 IF (INFO(6) .EQ. 1) RETURN
      DO 90 L = 1,NM
         WORK(JM+L) = WORK(JM+L) - WORK(JMB+L)
   90    CONTINUE
      IF (M .EQ. 2) GO TO 130
      MM1 = M - 1
C-----------------------------------------------------------------------
C   FOR THE PROBLEMS ON LEVEL K, K=M-1,...,2 THE CURRENT APPROXIMATION
C   VK IS IMPROVED USING (A PROLONGATION OF) THE APPROXIMATION ON THE
C   NEXT COARSER GRID. USE VK AS STARTING VALUE IN THE CHEBYSHEV PROCESS
C   ON LEVEL K
C-----------------------------------------------------------------------
      DO 120 IREV = 2,MM1
         I = M + 1 - IREV
         NXI = IWORK(I)
         NYI = IWORK(M+I)
         NI = NXI*NYI
         NXIP1 = IWORK(I+1)
         NYIP1 = IWORK(M+I+1)
         M3P5I = M3 + 5*I
         JI = IWORK(M3P5I-5)
         JIP1 = IWORK(M3P5I)
         JIB = IWORK(M3P5I-4)
         JIA1 = IWORK(M3P5I-3)
         JIA2 = IWORK(M3P5I-2)
         JIS = IWORK(M3P5I-1)
         CALL PROLON (NXIP1, NYIP1, NXI, NYI, X(IWORK(M8P2+I)),
     *                Y(IWORK(M9P2+I)), WORK(JIP1+1), WORK(JI+1))
         DO 100 L = 1,NI
            WORK(JI+L) = WORK(JI+L) + WORK(JIB+L)
  100       CONTINUE
         DI = 4.5*DV**(I - 1)
         ALFAI = .25*ALFAV**(I - 1)
         CALL CHEB (NXI, NYI, M, X(IWORK(M8P2+I)), Y(IWORK(M9P2+I)), T,
     *              TAU, WORK(JI+1), WORK(JIB+1), WORK(JIA1+1),
     *              WORK(JIA2+1), WORK(JIS+1), G, SUSE(I), METH, WORK,
     *              INFO, DI, ALFAI, I)
         IF (INFO(6) .EQ. 1) RETURN
         DO 110 L = 1,NI
            WORK(JI+L) = WORK(JI+L) - WORK(JIB+L)
  110       CONTINUE
  120    CONTINUE
C-----------------------------------------------------------------------
C   IMPROVEMENT OF THE FINEST-GRID APPROXIMATION
C-----------------------------------------------------------------------
  130 CALL PROLON (NX2, NY2, NX, NY, X, Y, WORK(IWORK(M3+5)+1), V)
      DO 140 L = 1,N1
         V(L,1) = V(L,1) + WORK(JB+L)
  140    CONTINUE
C-----------------------------------------------------------------------
C   CHEBYSHEV RELAXATION ON THE FINEST GRID LEVEL
C-----------------------------------------------------------------------
      CALL CHEB (NX, NY, M, X, Y, T, TAU, V, WORK(JB+1), WORK(JA1+1),
     *           WORK(JN+1), WORK(JNM1+1), G, SUSE(1), METH, WORK,
     *           INFO, 4.5, 0.25, 1)
      IF (INFO(6) .EQ. 1) RETURN
      IF (ITER .EQ. 2) GO TO 160
      DO 150 L = 1,N1
         WORK(JB+L) = V(L,1)
  150    CONTINUE
  160 CONTINUE
      RETURN
      END
      SUBROUTINE CHEB (NX, NY, M, X, Y, T, TAU, V, VB, VA1, VA2, VSIG,
     *                 G, SK, METH, WORK, INFO, D, ALFA, DEPTH)
      DIMENSION X(1),Y(1),V(NX,NY),VB(NX,NY),VA1(NX,NY),VA2(NX,NY),
     *          VSIG(NX,NY),WORK(1),INFO(1)
      COMMON/COEF/A1,A2,B0
      COMMON/OFLOW/VNNORM,FNNORM
      COMMON/MACH/ZETA
      INTEGER DEPTH
      REAL MU,LAMBDA
      EXTERNAL G
C-----------------------------------------------------------------------
C   CHEB PERFORMS A RELAXATION PROCESS BASED ON CHEBYSHEV ITERATION
C-----------------------------------------------------------------------
      ETA = 5.0*(VNNORM + TAU*FNNORM)
      DELTA = ALFA*SK
      W1 =  -2.0/(B0*TAU*(SK - DELTA))
      W0 = (SK + DELTA)/(SK - DELTA) - W1
      LAMBDA = W1/W0
      CH2 = W0
      CH1 = 2.0*W0*W0 - 1.0
      IF (DEPTH .EQ. 1) GO TO 30
      DO 10 J = 1,NY
         DO 10 I = 1,NX
         VA2(I,J) = V(I,J)
   10    CONTINUE
      CALL DERIV (M+1-DEPTH, NX, NY, X, Y, T+TAU, VA2, V, G, METH, WORK)
      INFO(11+M-DEPTH) = INFO(11+M-DEPTH) + 1
      VNORM = 0.0
      DO 20 J = 1,NY
         DO 20 I = 1,NX
         AUX = VA2(I,J) + LAMBDA*(VA2(I,J) - B0*TAU*V(I,J) - VSIG(I,J))
         VNORM = AMAX1(ABS(AUX),VNORM)
         VA1(I,J) = AUX
   20    CONTINUE
      IF (VNORM .LT. AMAX1(ETA,ZETA)) GO TO 60
      INFO(6) = 1
      RETURN
   30 DO 40 J = 1,NY
         DO 40 I = 1,NX
         VB(I,J) = V(I,J)
   40    CONTINUE
      CALL DERIV (M, NX, NY, X, Y, T+TAU, VB, V, G, METH, WORK)
      INFO(M+10) = INFO(M+10) + 1
      VNORM = 0.0
      DO 50 J = 1,NY
         DO 50 I = 1,NX
         AUX = VB(I,J) + LAMBDA*(VB(I,J) - B0*TAU*V(I,J) -
     *         (A1*VA2(I,J) + A2*VSIG(I,J)))
         VNORM = AMAX1(ABS(AUX),VNORM)
         V(I,J) = AUX
   50    CONTINUE
      IF (VNORM .LT. AMAX1(ETA,ZETA)) GO TO 60
      INFO(6) = 1
      RETURN
   60 MU = 2.0*W0*CH2/CH1
      IF (DEPTH .EQ. 1) GO TO 80
      CALL DERIV (M+1-DEPTH, NX, NY, X, Y, T+TAU, VA1, V, G, METH, WORK)
      INFO(11+M-DEPTH) = INFO(11+M-DEPTH) + 1
      VNORM = 0.0
      DO 70 J = 1,NY
         DO 70 I = 1,NX
         AUX = MU*VA1(I,J) + (1.0 - MU)*VA2(I,J) +
     *         MU*LAMBDA*(VA1(I,J) - B0*TAU*V(I,J) - VSIG(I,J))
         VNORM = AMAX1(ABS(AUX),VNORM)
         V(I,J) = AUX
   70    CONTINUE
      IF (VNORM .LT. AMAX1(ETA,ZETA)) GO TO 100
      INFO(6) = 1
      RETURN
   80 CALL DERIV (M, NX, NY, X, Y, T+TAU, V, VA1, G, METH, WORK)
      INFO(M+10) = INFO(M+10) + 1
      VNORM = 0.0
      DO 90 J = 1,NY
         DO 90 I = 1,NX
         AUX = MU*V(I,J) + (1.0 - MU)*VB(I,J) + MU*LAMBDA*(V(I,J) -
     *         B0*TAU*VA1(I,J) - (A1*VA2(I,J) + A2*VSIG(I,J)))
         VNORM = AMAX1(ABS(AUX),VNORM)
         V(I,J) = AUX
   90    CONTINUE
      IF (VNORM .LT. AMAX1(ETA,ZETA)) GO TO 120
      INFO(6) = 1
      RETURN
  100 IF (CH1 .GT. D) GO TO 120
      CH = 2.0*W0*CH1 - CH2
      CH2 = CH1
      CH1 = CH
      DO 110 J = 1,NY
         DO 110 I = 1,NX
         VA2(I,J) = VA1(I,J)
         VA1(I,J) = V(I,J)
  110    CONTINUE
      GO TO 60
  120 CONTINUE
      RETURN
      END
      SUBROUTINE RKR (NX, NY, X, Y, T, TAU, V, VB, VA1, VA2, VSIG, G,
     *                SM, METH, WORK, INFO)
      DIMENSION V(NX,NY),VB(NX,NY),VA1(NX,NY),VA2(NX,NY),VSIG(NX,NY),
     *          WORK(1),INFO(1),X(1),Y(1)
      COMMON/COEF/A1,A2,B0
      COMMON/OFLOW/VNNORM,FNNORM
      COMMON/MACH/ZETA
      REAL MU,LAMBDA
      EXTERNAL G
      DATA PI/3.1415926535898/
C-----------------------------------------------------------------------
C   RKR IS USED TO SOLVE THE COARSEST-GRID PROBLEM DURING THE FIRST
C   SWEEP. IT IS BASED ON AN EXPLICIT, STABILIZED RUNGE-KUTTA METHOD.
C-----------------------------------------------------------------------
      ITER = SQRT(B0*TAU*SM/1.50) + 1.0
      IF (ITER .EQ. 1) ITER = 2
      W0 = 1.0 + 0.05/FLOAT(ITER*ITER)
      W1 = COS(PI/FLOAT(2*ITER)) - W0
      LAMBDA = W1/W0
      CH2 = W0
      CH1 = 2.0*W0*W0 - 1.0
      TH2 = 0.0
      TH1 = -LAMBDA
      DO 10 J = 1,NY
         DO 10 I = 1,NX
         VA2(I,J) = V(I,J)
   10    CONTINUE
      CALL DERIV (1, NX, NY, X, Y, T, VA2, V, G, METH, WORK)
      INFO(11) = INFO(11) + 1
      VNNORM = 0.0
      FNNORM = 0.0
      DO 20 J = 1,NY
         DO 20 I = 1,NX
         VNNORM = AMAX1(ABS(VA2(I,J)),VNNORM)
         FNNORM = AMAX1(ABS(V(I,J)),FNNORM)
   20    CONTINUE
      VNORM = 0.0
      DO 30 J = 1,NY
         DO 30 I = 1,NX
         AUX = VA2(I,J) + LAMBDA*(VA2(I,J) - B0*TAU*V(I,J) - VSIG(I,J))
         VNORM = AMAX1(ABS(AUX),VNORM)
         VA1(I,J) = AUX
   30    CONTINUE
      IF (VNORM .LT. AMAX1(5.0*(VNNORM + TH1*TAU*FNNORM),ZETA)) GO TO 40
      INFO(6) = 1
      RETURN
   40 ITERM1 = ITER - 1
      DO 80 K = 1,ITERM1
         CALL DERIV (1, NX, NY, X, Y, T+TH1*TAU, VA1, V, G, METH, WORK)
         INFO(11) = INFO(11) + 1
         MU = 2.0*W0*CH2/CH1
         TH = MU*TH1 + (1.0 - MU)*TH2 + MU*LAMBDA*(TH1 - 1.0)
         VNORM = 0.0
         DO 50 J = 1,NY
            DO 50 I = 1,NX
            AUX = MU*VA1(I,J) + (1.0 - MU)*VA2(I,J) +
     *            MU*LAMBDA*(VA1(I,J) - B0*TAU*V(I,J) - VSIG(I,J))
            VNORM = AMAX1(ABS(AUX),VNORM)
            V(I,J) = AUX
   50       CONTINUE
         IF (VNORM .LT. AMAX1(5.0*(VNNORM + TH1*TAU*FNNORM),ZETA))
     *       GO TO 60
         INFO(6) = 1
         RETURN
   60    IF (K .EQ. ITERM1) GO TO 80
         CH = 2.0*W0*CH1 - CH2
         DO 70 J = 1,NY
            DO 70 I = 1,NX
            VA2(I,J) = VA1(I,J)
            VA1(I,J) = V(I,J)
   70       CONTINUE
         TH2 = TH1
         TH1 = TH
         CH2 = CH1
         CH1 = CH
   80    CONTINUE
      RETURN
      END
      SUBROUTINE RESTR (NXF, NYF, NXC, NYC, X, Y, VF, VC)
      DIMENSION X(1),Y(1),VF(NXF,NYF),VC(NXC,NYC)
C-----------------------------------------------------------------------
C   RESTR PERFORMS A RESTRICTION OPERATOR, NECESSARY IN THE MULTIGRID
C   TECHNIQUE. IT TRANSFORMS A GRID FUNCTION VF DEFINED ON SOME GRID
C   INTO A GRID FUNCTION VC, DEFINED ON THE NEXT COARSER GRID. THE
C   GRID LINES OF THE COARSE GRID COINCIDE WITH THOSE OF THE FINE GRID.
C-----------------------------------------------------------------------
      NXC1 = NXC - 1
      NYC1 = NYC - 1
C-----------------------------------------------------------------------
C   AT INTERNAL POINTS, RESTR USES A NINE-POINT MOLECULE
C-----------------------------------------------------------------------
      DO 10 JC = 2,NYC1
         JF = 2*JC - 1
         BETA = (Y(JF) - Y(JF-1))/(Y(JF+1) - Y(JF))
         DO 10 IC = 2,NXC1
         IF = 2*IC - 1
         ALPHA = (X(IF) - X(IF-1))/(X(IF+1) - X(IF))
         VC(IC,JC) = (VF(IF,JF) +
     *               (VF(IF-1,JF) + ALPHA*VF(IF+1,JF))/(1.0 + ALPHA) +
     *               (VF(IF,JF-1) + BETA*VF(IF,JF+1))/(1.0 + BETA) +
     *               (VF(IF-1,JF-1) + ALPHA*VF(IF+1,JF-1) +
     *                BETA*VF(IF-1,JF+1) + ALPHA*BETA*VF(IF+1,JF+1))/
     *               ((1.0 + ALPHA)*(1.0 + BETA)))/4.0
   10    CONTINUE
C-----------------------------------------------------------------------
C   THE POINTS AT THE BOUNDARY ARE TREATED SEPARATELY
C-----------------------------------------------------------------------
      DO 20 JC = 1,NYC,NYC1
         JF = 2*JC - 1
         DO 20 IC = 1,NXC
         IF = 2*IC - 1
         VC(IC,JC) = VF(IF,JF)
   20    CONTINUE
      DO 30 JC = 2,NYC1
         JF = 2*JC - 1
         DO 30 IC = 1,NXC,NXC1
         IF = 2*IC - 1
         VC(IC,JC) = VF(IF,JF)
   30    CONTINUE
      RETURN
      END
      SUBROUTINE INJECT (NXF, NYF, NXC, NYC, VF, VC)
      DIMENSION VF(NXF,NYF),VC(NXC,NYC)
C-----------------------------------------------------------------------
C   INJECT PERFORMS THE SIMPLEST RESTRICTION OPERATOR, CALLED INJECTION
C   BY WHICH THE ELEMENTS OF THE COARSE-GRID FUNCTION ARE OBTAINED FROM
C   COPYING THE CORRESPONDING ELEMENTS OF THE FINE-GRID FUNCTION
C-----------------------------------------------------------------------
      DO 10 JC = 1,NYC
         JF = 2*JC - 1
         DO 10 IC = 1,NXC
         IF = 2*IC - 1
         VC(IC,JC) = VF(IF,JF)
   10    CONTINUE
      RETURN
      END
      SUBROUTINE PROLON (NXC, NYC, NXF, NYF, X, Y, VC, VF)
      DIMENSION X(1),Y(1),VC(NXC,NYC),VF(NXF,NYF)
C-----------------------------------------------------------------------
C   PROLON PERFORMS A PROLONGATION OPERATOR, USED IN THE MULTIGRID
C   METHOD. STARTING WITH THE COARSE-GRID FUNCTION VC IT DELIVERS IN VF
C   THE GRID FUNCTION DEFINED ON THE NEXT FINER GRID. PROLON USES LINEAR
C   INTERPOLATION AT INTERNAL POINTS AND LINEAR EXTRAPOLATION IN POINTS
C   ADJACENT TO THE BOUNDARY. IN POINTS WHICH COINCIDE ON BOTH GRIDS,
C   THE VALUES OF VF ARE COPIED FROM THE CORRESPONDING VALUES OF THE
C   COARSE-GRID FUNCTION VC
C-----------------------------------------------------------------------
      NXFM1 = NXF - 1
      NXFM2 = NXF - 2
      NXFM3 = NXF - 3
      NYFM1 = NYF - 1
      NYFM2 = NYF - 2
      NYFM3 = NYF - 3
C-----------------------------------------------------------------------
C   PROLONGATION OF THE POINTS ON THE BOUNDARY
C-----------------------------------------------------------------------
      DO 30 JF = 1,NYF,NYFM1
         JC = (JF + 1)/2
         DO 10 IF = 1,NXF,2
            IC = (IF + 1)/2
            VF(IF,JF) = VC(IC,JC)
   10       CONTINUE
         DO 20 IF = 2,NXFM1,2
            IC = IF/2
            VF(IF,JF) = (VC(IC+1,JC)*(X(IF) - X(IF-1)) + VC(IC,JC)*
     *                  (X(IF+1) - X(IF)))/(X(IF+1) - X(IF-1))
   20       CONTINUE
   30    CONTINUE
      DO 60 IF = 1,NXF,NXFM1
         IC = (IF + 1)/2
         DO 40 JF = 3,NYFM2,2
            JC = (JF + 1)/2
            VF(IF,JF) = VC(IC,JC)
   40       CONTINUE
         DO 50 JF = 2,NYFM1,2
            JC = JF/2
            VF(IF,JF) = (VC(IC,JC+1)*(Y(JF) - Y(JF-1)) + VC(IC,JC)*
     *                  (Y(JF+1) - Y(JF)))/(Y(JF+1) - Y(JF-1))
   50       CONTINUE
   60    CONTINUE
C-----------------------------------------------------------------------
C   PROLONGATION OF THE INTERIOR POINTS, USING LINEAR INTERPOLATION
C-----------------------------------------------------------------------
      DO 90 JF = 3,NYFM2,2
         JC = (JF + 1)/2
         DO 70 IF = 3,NXFM2,2
            IC = (IF + 1)/2
            VF(IF,JF) = VC(IC,JC)
   70       CONTINUE
         DO 80 IF = 4,NXFM3,2
            IC = IF/2
            VF(IF,JF) = (VC(IC+1,JC)*(X(IF) - X(IF-1)) + VC(IC,JC)*
     *                  (X(IF+1) - X(IF)))/(X(IF+1) - X(IF-1))
   80       CONTINUE
   90    CONTINUE
      DO 120 JF = 4,NYFM3,2
         JC = JF/2
         DO 100 IF = 3,NXFM2,2
            IC = (IF + 1)/2
            VF(IF,JF) = (VC(IC,JC+1)*(Y(JF) - Y(JF-1)) + VC(IC,JC)*
     *                  (Y(JF+1) - Y(JF)))/(Y(JF+1) - Y(JF-1))
  100       CONTINUE
         DO 110 IF = 4,NXFM3,2
            VF(IF,JF) = (VF(IF+1,JF)*(X(IF) - X(IF-1)) +
     *                VF(IF-1,JF)*(X(IF+1) - X(IF)))/(X(IF+1) - X(IF-1))
  110       CONTINUE
  120    CONTINUE
C-----------------------------------------------------------------------
C   CALCULATION OF VF-VALUES CORRESPONDING TO POINTS LYING ON GRID LINES
C   ADJACENT TO THE BOUNDARY, USING EXTRAPOLATION
C-----------------------------------------------------------------------
      DY1 = Y(4) - Y(3)
      DY2 = Y(3) - Y(2)
      EY1 = Y(NYFM2) - Y(NYFM3)
      EY2 = Y(NYFM1) - Y(NYFM2)
      DO 130 IF = 3,NXFM2
        VF(IF,2) = (VF(IF,3)*(DY1 + DY2) - VF(IF,4)*DY2)/DY1
        VF(IF,NYFM1) = (VF(IF,NYFM2)*(EY1 + EY2) - VF(IF,NYFM3)*EY2)/EY1
  130   CONTINUE
      DX1 = X(4) - X(3)
      DX2 = X(3) - X(2)
      EX1 = X(NXFM2) - X(NXFM3)
      EX2 = X(NXFM1) - X(NXFM2)
      DO 140 JF = 3,NYFM2
        VF(2,JF) = (VF(3,JF)*(DX1 + DX2) - VF(4,JF)*DX2)/DX1
        VF(NXFM1,JF) = (VF(NXFM2,JF)*(EX1 + EX2) - VF(NXFM3,JF)*EX2)/EX1
  140   CONTINUE
      VF(2,2) = ((VF(3,2)*(DX1 + DX2) - VF(4,2)*DX2)/DX1 +
     *         (VF(2,3)*(DY1 + DY2) - VF(2,4)*DY2)/DY1)/2.0
      VF(NXFM1,2) = ((VF(NXFM2,2)*(EX1 + EX2) - VF(NXFM3,2)*EX2)/EX1 +
     *             (VF(NXFM1,3)*(DY1 + DY2) - VF(NXFM1,4)*DY2)/DY1)/2.0
      VF(2,NYFM1) = ((VF(3,NYFM1)*(DX1 + DX2) - VF(4,NYFM1)*DX2)/DX1 +
     *             (VF(2,NYFM2)*(EY1 + EY2) - VF(2,NYFM3)*EY2)/EY1)/2.0
      VF(NXFM1,NYFM1) = ((VF(NXFM2,NYFM1)*(EX1 + EX2) - VF(NXFM3,NYFM1)*
     *                 EX2)/EX1 + (VF(NXFM1,NYFM2)*(EY1 + EY2) -
     *                 VF(NXFM1,NYFM3)*EY2)/EY1)/2.0
      RETURN
      END
      FUNCTION TAUST (NX, NY, M, X, Y, T, V, G, TOL, METH, WORK, IWORK,
     *                INFO)
      DIMENSION X(1),Y(1),V(NX,NY),WORK(1),IWORK(1),INFO(1)
      EXTERNAL G
C-----------------------------------------------------------------------
C   TAUST CALCULATES AN INITIAL STEP SIZE
C-----------------------------------------------------------------------
      JA1 = IWORK(3*M+2)
      N = NX*NY
      CALL DERIV (M, NX, NY, X, Y, T, V, WORK(JA1+1), G, METH, WORK)
      INFO(M+10) = INFO(M+10) + 1
      VNORM = 0.0
      FNORM = 0.0
      DO 10 L = 1,N
         VNORM = VNORM + V(L,1)**2
         FNORM = FNORM + WORK(JA1+L)**2
   10    CONTINUE
      VNORM = SQRT(VNORM/FLOAT(N))
      FNORM = AMAX1(SQRT(FNORM/FLOAT(N)),1.0E-10)
      TAUST = 0.1*SQRT((TOL + TOL*VNORM)/FNORM)
      RETURN
      END
      SUBROUTINE SPRI (NX, NY, M, X, Y, T, TAU, V, G, SPRAD, SEST, SUSE,
     *                 METH, WORK, IWORK, INFO, Q)
      DIMENSION X(1),Y(1),V(NX,NY),SEST(1),SUSE(1),
     *          WORK(1),IWORK(1),INFO(1)
      EXTERNAL G
C-----------------------------------------------------------------------
C   SPRI DELIVERS THE SPECTRAL RADIUS SK OF THE JACOBIAN MATRIX OF THE
C   RIGHT-HAND SIDE FUNCTION G, DEFINED ON LEVEL K, K=1,...,M.
C   THERE ARE THREE DIFFERENT OPTIONS AVAILABLE WITH RESPECT TO THE
C   SPECTRAL RADIUS (SEE THE DOCUMENTATION IN SUBROUTINE BDMG).
C   SPRI DELIVERS TWO VALUES CONCERNING SK:
C   SEST(K), BEING THE CALCULATED ESTIMATE OF SK AT TIME T ON LEVEL K
C   SUSE(K), THE VALUE OF SK WHICH WILL BE ACTUALLY USED IN THE
C            RELAXATION PROCESS ON LEVEL K
C-----------------------------------------------------------------------
      M3 = M*3
      M2 = M*2
      M8 = M*8
      M8P2 = M8 + 2
      M9P2 = M8P2 + M
      M10P2 = M9P2 + M
      MM1 = M - 1
      IF (INFO(3) .LT. 3) GO TO 20
C-----------------------------------------------------------------------
C   IN CASE OF THE NON-AUTOMATIC OPTION, THE SPECTRAL RADIUS IS
C   DETERMINED BY MEANS OF THE USER-SUPPLIED FUNCTION SPRAD
C-----------------------------------------------------------------------
      SUSE(1) = SPRAD (M, NX, NY, X, Y, T, WORK(IWORK(M3+3)+1), TAU)
      DO 10 I = 2,M
         NXI = IWORK(I)
         NYI = IWORK(M+I)
         NXIM1 = IWORK(I-1)
         NYIM1 = IWORK(M+I-1)
         M3P5I = M3 + 5*I
         JIN = IWORK(M3P5I-2)
         JIM1N = IWORK(M3P5I-7)
         CALL RESTR (NXIM1, NYIM1, NXI, NYI, X(IWORK(M8P2+I-1)),
     *               Y(IWORK(M9P2+I-1)), WORK(JIM1N+1), WORK(JIN+1))
         SUSE(I) = SPRAD (M+1-I, NXI, NYI, X(IWORK(M8P2+I)),
     *                    Y(IWORK(M9P2+I)), T, WORK(JIN+1), TAU)
   10    CONTINUE
      RETURN
C-----------------------------------------------------------------------
C   AUTOMATIC OR SEMI-AUTOMATIC OPTION, WHICH MEANS THAT SPRI HAS TO
C   ESTIMATE THE SPECTRAL RADIUS
C-----------------------------------------------------------------------
   20 IEND = M
      IF (INFO(6) .EQ. 1) IEND = MM1
      DO 30 I = 2,IEND
         CALL RESTR (IWORK(I-1), IWORK(M+I-1), IWORK(I), IWORK(M+I),
     *               X(IWORK(M8P2+I-1)), Y(IWORK(M9P2+I-1)),
     *               WORK(IWORK(M3+5*I-7)+1), WORK(IWORK(M3+5*I-2)+1))
   30    CONTINUE
      IF (INFO(6) .EQ. 1) GO TO 50
C-----------------------------------------------------------------------
C   ESTIMATION OF THE SPECTRAL RADIUS ON THE COARSEST LEVEL
C-----------------------------------------------------------------------
      S = GERSG (IWORK(M), IWORK(M2), 1, X(IWORK(M9P2)),
     *           Y(IWORK(M10P2)), T, WORK(IWORK(M8-2)+1), G, METH,
     *           WORK, WORK(IWORK(M8-4)+1), WORK(IWORK(M8-3)+1))
      INFO(M+11) = INFO(M+11) + 6
      IF (INFO(7) .EQ. 0) SEST(M) = S
      FAC = S/SEST(M)
      IF (FAC.GT.1.15  .OR. FAC.LT.0.75 .OR. INFO(7).EQ.0) GO TO 50
C-----------------------------------------------------------------------
C   IF THE SPECTRAL RADIUS ON THE COARSEST LEVEL CHANGES SLOWLY WITH
C   TIME, THE SPECTRAL RADII ON ALL FINER LEVELS ARE UPDATED ACCORDING
C   TO THIS SMALL CHANGE
C-----------------------------------------------------------------------
      ZMULT = 1.0
      IF (FAC .GT. 1.0) ZMULT = (Q + 1.0 - 1.0/FAC)/Q
      DO 40 I = 1,M
         SEST(I) = SEST(I)*FAC
         SUSE(I) = SEST(I)*ZMULT
   40    CONTINUE
      INFO(5) = 0
      GO TO 170
C-----------------------------------------------------------------------
C   IF THE SPECTRAL RADIUS ON THE COARSEST LEVEL CHANGES RAPIDLY WITH
C   TIME, THE SPECTRAL RADII ON ALL FINER LEVELS ARE RECALCULATED
C-----------------------------------------------------------------------
   50 IIBG = 2
      IF (M .EQ. 2) IIBG = 1
      IF (INFO(6) .EQ. 1) S = SEST(M)
      DO 100 II = IIBG,M
         I = M + IIBG - II
         IF (I .EQ. M) GO TO 60
         NXI = IWORK(I)
         NYI = IWORK(M+I)
         M3P5I = M3 + 5*I
         JIN = IWORK(M3P5I-2)
         JIB = IWORK(M3P5I-4)
         JIA1 = IWORK(M3P5I-3)
         S = GERSG (NXI, NYI, M+1-I, X(IWORK(M8P2+I)), Y(IWORK(M9P2+I)),
     *              T, WORK(JIN+1), G, METH, WORK, WORK(JIB+1),
     *              WORK(JIA1+1))
         INFO(11+M2-I) = INFO(11+M2-I) + 6
   60    IF (INFO(7) .EQ. 0) GO TO 70
         IF (INFO(6) .EQ. 1) GO TO 80
         SUSE(I) = AMAX1(S,S + (S - SEST(I))/Q)
         GO TO 90
   70    SUSE(I) = S
         GO TO 90
   80    SUSE(I) = AMAX1(1.2*SUSE(I),SUSE(I) + (S - SEST(I))*
     *                                                (Q + 1.0)/Q)
   90    SEST(I) = S
  100    CONTINUE
      INFO(5) = 1
      IF (M .EQ. 2) GO TO 170
C-----------------------------------------------------------------------
C   IF WE ARE DEALING WITH MORE THAN TWO GRIDS, THE SPECTRAL RADIUS ON
C   THE FINEST LEVEL IS EXTRAPOLATED FROM THE VALUES OF THE SPECTRAL
C   RADII ON THE UNDERLYING COARSER LEVELS
C-----------------------------------------------------------------------
      DO 110 I = 2,M
         IF (SEST(I) .EQ. 0.0) SEST(I) = 1.0E-10
  110    CONTINUE
      S = 0.0
      DO 130 I = 2,MM1
         RLI = 1.0
         DO 120 J = 2,MM1
            IF (I .EQ. J) GO TO 120
            RLI = RLI*FLOAT(1 - J)/FLOAT(I - J)
  120       CONTINUE
         S = S + RLI*SEST(I)/SEST(I+1)
  130    CONTINUE
      SP = SEST(2)*AMAX1(S,4.0)
      IF (INFO(7) .EQ. 0) GO TO 140
      IF (INFO(6) .EQ. 1) GO TO 150
      SUSE(1) = AMAX1(SP,SP + (SP - SEST(1))/Q)
      GO TO 160
  140 SUSE(1) = SP
      GO TO 160
  150 SUSE(1) = AMAX1(1.2*SUSE(1),SUSE(1) + (SP - SEST(1))*
     *                                (Q + 1.0)/Q)
  160 SEST(1) = SP
C-----------------------------------------------------------------------
C   PROVISIONS IN ORDER TO AVOID ZERO VALUES FOR THE SPECTRAL RADII
C-----------------------------------------------------------------------
  170 DO 180 I = 1,M
         IF (SEST(I) .EQ. 0.0) SEST(I) = 1.0E-10
         IF (SUSE(I) .EQ. 0.0) SUSE(I) = 1.0E-10
  180    CONTINUE
      IF (INFO(6) .EQ. 1) INFO(6) = 0
      RETURN
      END
      FUNCTION GERSG (NX, NY, K, X, Y, T, U, G, METH, WORK, FU, FUR)
      DIMENSION X(1),Y(1),U(NX,NY),WORK(1),FU(NX,NY),FUR(NX,NY)
      EXTERNAL G
C-----------------------------------------------------------------------
C   GERSG ESTIMATES THE SPECTRAL RADIUS AT GRID LEVEL K. IT USES
C   GERSCHGORIN'S DISK THEOREM, ASSUMING A FIVE-DIAGONAL STRUCTURE IN
C   THE JACOBIAN MATRIX. GERSG USES SIX EVALUATIONS OF THE DERIVATIVE
C   FUNCTION AT LEVEL K
C-----------------------------------------------------------------------
      FAC = 1.0E-6
      GERSG = 0.0
      CALL DERIV (K, NX, NY, X, Y, T, U, FU, G, METH, WORK)
      DO 110 L = 1,5
         DO 20 J = 1,NY
            IST = 1 + MOD(L + 2*J,5)
            IF (IST .GT. NX) GO TO 20
            DO 10 I = IST,NX,5
               R = U(I,J)*FAC
               IF (R .EQ. 0.0) R = FAC
               U(I,J) = U(I,J) + R
   10          CONTINUE
   20       CONTINUE
         CALL DERIV (K, NX, NY, X, Y, T, U, FUR, G, METH, WORK)
         DO 100 J = 1,NY
            IST = 1 + MOD(L + 2*J,5)
            IF (IST .GT. NX) GO TO 100
            DO 90 I = IST,NX,5
               S = ABS(FUR(I,J) - FU(I,J))
               IF (I .EQ. 1) GO TO 30
               S = S + ABS(FUR(I-1,J) - FU(I-1,J))
   30          IF (I .EQ. NX) GO TO 40
               S = S + ABS(FUR(I+1,J) - FU(I+1,J))
   40          IF (J .EQ. 1) GO TO 50
               S = S + ABS(FUR(I,J-1) - FU(I,J-1))
   50          IF (J .EQ. NY) GO TO 60
               S = S + ABS(FUR(I,J+1) - FU(I,J+1))
   60          IF (U(I,J) .EQ. 0.0 .OR. U(I,J) .EQ. FAC) GO TO 70
               R = FAC*U(I,J)/(1.0 + FAC)
               GO TO 80
   70          R = FAC
   80          S = ABS(S/R)
               IF (S .GT. GERSG) GERSG = S
               U(I,J) = U(I,J) - R
   90          CONTINUE
  100       CONTINUE
  110    CONTINUE
      RETURN
      END
      SUBROUTINE LTE (NX, NY, M, X, Y, T, TAUNEW, V, G, TOL, METH, WORK,
     *                IWORK, INFO, QOLD, Q, EPS, ERROR)
      DIMENSION X(1),Y(1),V(NX,NY),WORK(1),IWORK(1),INFO(1)
      EXTERNAL G
C-----------------------------------------------------------------------
C   LTE ESTIMATES THE LOCAL TRUNCATION ERROR (PARAMETER ERROR) AND
C   CALCULATES THE LOCAL ERROR BOUND (PARAMETER EPS) TO BE USED IN THE
C   MIXED ERROR TEST
C-----------------------------------------------------------------------
      M3 = M*3
      J1A1 = IWORK(M3+2)
      J1N = IWORK(M3+3)
      J1NM1 = IWORK(M3+4)
      J2 = IWORK(M3+5)
      J2B = IWORK(M3+6)
      J2A1 = IWORK(M3+7)
      J2A2 = IWORK(M3+8)
      J2AUX = IWORK(8*M)
      NX2 = IWORK(2)
      NX2M1 = NX2 - 1
      NY2 = IWORK(M+2)
      NY2M1 = NY2 - 1
C-----------------------------------------------------------------------
C   DISTINGUISH BETWEEN THE FIRST, THE SECOND AND SUBSEQUENT STEPS
C-----------------------------------------------------------------------
      NSUCST = INFO(7) - INFO(8)
      IF (NSUCST .EQ. 1) GO TO 20
      IF (NSUCST .EQ. 2) GO TO 10
C-----------------------------------------------------------------------
C   SET COEFFICIENTS FOR A THIRD ORDER BDF, USED AS A REFERENCE SOLUTION
C-----------------------------------------------------------------------
      R = Q*(1.0 + QOLD)
      A3 = (Q*(Q + 2.0) + 1.0)/((Q + 2.0)*R*(R*R - 3.0) + R*(R + 2.0)*
     *     (3.0 - Q*Q))
      A2 = ( -1.0 - R*(R + 2.0)*A3)/(Q*(Q + 2.0))
      B0 = 1.0 + Q*A2 + R*A3
      A1 = 1.0 - A2 - A3
      COEF = A3
      GO TO 30
C-----------------------------------------------------------------------
C   SET COEFFICIENTS FOR A THIRD ORDER, TWO-STEP LINEAR MULTISTEP
C   FORMULA, USED AS A REFERENCE SOLUTION IN THE SECOND STEP
C-----------------------------------------------------------------------
   10 B0 = (Q + 1.0)/(2.0*Q + 3.0)
      A2 = 1.0/(Q*Q*(2.0*Q + 3.0))
      B1 = (Q + 1.0)*B0/Q
      A1 = 1.0 - A2
      COEF = B1*TAUNEW
      GO TO 30
C-----------------------------------------------------------------------
C   SET COEFFICIENTS FOR THE SECOND ORDER TRAPEZOIDAL METHOD, USED AS A
C   REFERENCE SOLUTION IN THE FIRST STEP
C-----------------------------------------------------------------------
   20 A1 = 1.0
      A2 = 0.0
      B0 = 0.5
      COEF = B0*TAUNEW
   30 CALL INJECT (NX, NY, NX2, NY2, V, WORK(J2+1))
      CALL INJECT (NX, NY, NX2, NY2, WORK(J1N+1), WORK(J2B+1))
      CALL INJECT (NX, NY, NX2, NY2, WORK(J1A1+1), WORK(J2A2+1))
      IF (NSUCST .EQ. 1) GO TO 40
      CALL INJECT (NX, NY, NX2, NY2, WORK(J1NM1+1), WORK(J2A1+1))
      IF (NSUCST .GT. 2) GO TO 50
   40 CALL DERIV (M-1, NX2, NY2, X(IWORK(8*M+4)), Y(IWORK(9*M+4)), T,
     *            WORK(J2B+1), WORK(J2AUX+1), G, METH, WORK)
      INFO(M+9) = INFO(M+9) + 1
C-----------------------------------------------------------------------
C   CALCULATE THE LOCAL ERROR BOUND AND THE LOCAL TRUNCATION ERROR
C-----------------------------------------------------------------------
   50 EPS = 0.0
      ERROR = 0.0
      DO 60 J = 2,NY2M1
         DO 60 I = 2,NX2M1
         L = (J - 1)*NX2 + I
         VNP1 = WORK(J2+L)
         EPS = EPS + VNP1*VNP1
         E = VNP1 - A1*WORK(J2B+L) - A2*WORK(J2A1+L) -
     *       COEF*WORK(J2AUX+L) - B0*TAUNEW*WORK(J2A2+L)
         ERROR = ERROR + E*E
   60    CONTINUE
      EPS = TOL + TOL*SQRT(EPS/FLOAT(NX2M1*NY2M1))
      ERROR = SQRT(ERROR/FLOAT(NX2M1*NY2M1))
      RETURN
      END
      SUBROUTINE NEWSTP (EPS, ERROR, TAUOLD, TAUNEW, FIRST)
      LOGICAL FIRST
C-----------------------------------------------------------------------
C   BASED ON THE VALUES OF EPS AND ERROR, NEWSTP CALCULATES A STEP SIZE
C   FOR THE NEXT INTEGRATION STEP. THE RATIO OF THE NEW STEP SIZE AND
C   THE OLD ONE IS RESTRICTED TO THE INTERVAL [0.1,2.4]
C-----------------------------------------------------------------------
      IF (ERROR .NE. 0.0) GO TO 10
      TAUNEW = TAUOLD*2.4
      RETURN
   10 IP = 3
      IF (FIRST) IP = 2
      ALPHA = (EPS/ERROR)**(1.0/FLOAT(IP))
      IF (                    ALPHA.LT.0.18) TAUNEW=TAUOLD*0.1
      IF (ALPHA.GE.0.18 .AND. ALPHA.LT.1.0) TAUNEW=TAUOLD*ALPHA/1.8
      IF (ALPHA.GE. 1.0 .AND. ALPHA.LT.3.6) TAUNEW=TAUOLD*(.6*ALPHA+.24)
      IF (ALPHA.GE. 3.6                   ) TAUNEW=TAUOLD*2.4
      RETURN
      END
      SUBROUTINE FIT (T, TEND, TAU, LAST)
C-----------------------------------------------------------------------
C   FIT ADJUSTS THE CURRENT STEP SIZE TAU IN SUCH A WAY THAT THE
C   REMAINING INTEGRATION INTERVAL [T,TEND] IS A MULTIPLE OF TAU.
C   BY APPLYING THIS ADAPTATION TO EVERY INTEGRATION STEP, THE PROCESS
C   ARRIVES EXACTLY AT TEND; HENCE, INTERPOLATING THE NUMERICAL SOLUTION
C   AT THIS POINT CAN BE AVOIDED
C-----------------------------------------------------------------------
      LOGICAL LAST
      RNS = (TEND - T)/TAU
      NS = INT(RNS)
      IF (RNS - FLOAT(NS) .GT. 0.05 .OR. NS .EQ. 0) NS = NS + 1
      TAU = (TEND - T)/FLOAT(NS)
      LAST = NS .EQ. 1
      RETURN
      END
      SUBROUTINE SHIFT (NX, NY, M, T, TAUOLD, TAUNEW, V, WORK, IWORK)
      DIMENSION V(NX,NY),WORK(1),IWORK(1)
C-----------------------------------------------------------------------
C   IN SHIFT, THE SOLUTION VECTORS AND THE T-ARGUMENT ARE SHIFTED
C-----------------------------------------------------------------------
      JN = IWORK(3*M+3)
      JNM1 = IWORK(3*M+4)
      N = NX*NY
      DO 10 L = 1,N
         WORK(JNM1+L) = WORK(JN+L)
         WORK(JN+L) = V(L,1)
   10    CONTINUE
      T = T + TAUNEW
      TAUOLD = TAUNEW
      RETURN
      END
      FUNCTION SUMG (M, INFO)
      DIMENSION INFO(1)
C-----------------------------------------------------------------------
C   SUMG CALCULATES THE TOTAL NUMBER OF G-EVALUATIONS SPENT SO FAR.
C   THE G-EVALUATIONS ON THE COARSER GRIDS ARE CONVERTED INTO THE
C   EQUIVALENTS OF FINE-GRID EVALUATIONS. THIS TOTAL NUMBER INCLUDES
C   THE EVALUATIONS PERFORMED IN THE INTEGRATION PROCESS AS WELL AS THE
C   EVALUATIONS REQUIRED TO ESTIMATE AND CONTROL THE SPECTRAL RADIUS
C-----------------------------------------------------------------------
      SUMG = INFO(M+10) + INFO(2*M+10)
      DO 10 K = 2,M
         SUMG = SUMG + FLOAT(INFO(M+11-K)+INFO(2*M+11-K)) / (4.0**(K-1))
   10    CONTINUE
      RETURN
      END