[CONTACT]

[ABOUT]

[POLICY]

[ADVERTISE]

C ALGORITHM CC SOLUTION

Found at: ftp.icm.edu.pl:70/packages/netlib/toms-2014-06-10/565

C ALGORITHM 565
C
C PDETWO/PSTEM/GEARB: SOLUTION OF SYSTEMS OF TWO DIMENSIONAL
C NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS
C
C BY D.K. MELGAARD AND R.F. SINCOVEC
C
C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE 7,1 (MARCH 1981)
C
C                                                                        CPDE
C     MAIN PROGRAM TO ILLUSTRATE THE USAGE OF THE PDETWO PACKAGE         CPDE
C     USING EXAMPLES 1, 2, AND 3.                                        CPDE
C                                                                        CPDE
      REAL  H,S,X,ERR1,DX,Y,TOUT,REALN1,DY,U1,HUSED,                     CPDE
     *      EXP,R,T0,EPS,ERMAX,ABS,WORK,DL,DLI,DEV,                      CPDE
     *      U2,U3,ERR2,ERR3,REALN2,REALN3                                CPDE
      INTEGER IX,NPDE,NSTEP,NX,NFE,NODE,MF,NY,                           CPDE
     *        NJE,NQUSED,IY,INDEX,I,IWORK,KODE,IK                        CPDE
      COMMON /GEAR3/ HUSED,NQUSED,NSTEP,NFE,NJE                          CPDE
      COMMON /PROB/ DL,DLI,KODE                                          CPDE
      DIMENSION U1(10,10),ERR1(10,10),REALN1(10,10)                      CPDE
      DIMENSION U2(5,31),ERR2(5,31),REALN2(5,31)                         CPDE
      DIMENSION U3(2,5,5),ERR3(2,5,5),REALN3(2,5,5)                      CPDE
      DIMENSION WORK(4186),IWORK(155),X(10),Y(31)                        CPDE
      EQUIVALENCE (U1(1,1),U3(1,1,1)),(ERR1(1,1),ERR3(1,1,1)),           CPDE
     *            (REALN1(1,1),REALN3(1,1,1))                            CPDE
      EQUIVALENCE (U2(1,1),U3(1,1,1)),(ERR2(1,1),ERR3(1,1,1)),           CPDE
     *            (REALN2(1,1),REALN3(1,1,1))                            CPDE
C                                                                        CPDE
      DO 1000 IK=1,3                                                     CPDE
      KODE = IK                                                          CPDE
      IF(KODE.EQ.1)  GO TO 100                                           CPDE
      IF(KODE.EQ.2)  GO TO 300                                           CPDE
      IF(KODE.EQ.3)  GO TO 500                                           CPDE
 100  CONTINUE                                                           CPDE
C***********************************************************************
C
C  EXAMPLE 1.  ELLIPTIC PDE
C
C***********************************************************************
      WRITE (6,110)
  110 FORMAT(1H1,3X,36HOUTPUT FROM EXAMPLE 1.  ELLIPTIC PDE//)
C
C  DEFINE THE PROBLEM PARAMETERS.
C
      NX=10
      NY=10
      NPDE=1
      NODE=NPDE*NX*NY
      MF=22
      INDEX=1
      T0=0.0
      H=0.1E-06
      EPS=0.1E-06
      DX = 1.0/(FLOAT(NX)-1.0)
      DY = 1.0/(FLOAT(NY)-1.0)
      DO 120 IX=1,NX
        Y(IX)=FLOAT(IX)*DY-DY
  120   X(IX)=FLOAT(IX)*DX-DX
      IWORK(1) = NPDE
      IWORK(2) = NX
      IWORK(3) = NY
      IWORK(4) = 5
      IWORK(5) = 4186
      IWORK(6) = 100
C
C  CALCULATE THE ACTUAL SOLUTION
C
      DO 130 IY = 1,NY
        DO 130 IX = 1,NX
          R = X(IX)
          S = Y(IY)
  130     REALN1(IX,IY)=3.0*EXP(R)*EXP(S)*(R-R*R)*(S-S*S)
      WRITE (6,140) ((REALN1(IX,IY),IX=1,NX),IY=1,NY)
  140 FORMAT (/23H STEADY STATE SOLUTION /10(/3H    ,10E11.4))
C
C  DEFINE THE INITIAL CONDITION
C
      DO 160 IY = 1,NY
        DO 160 IX = 1,NX
  160    U1(IX,IY)=0.0
      WRITE (6,170) NODE,T0,H,EPS,MF,U1
  170 FORMAT (//6H NODE ,I3,4H T0 ,F9.2,3H H ,E8.1,5H EPS ,E8.1,4H MF ,
     *        I2 //16H INITIAL VALUES / 10(/3H    ,10E11.2))
C
C  SET UP THE LOOP FOR CALLING THE INTEGRATOR AT DIFFERENT TOUT VALUES
C
      TOUT=.001
      DO 260 I=1,5
      WRITE (6,200) TOUT
  200 FORMAT (// 5H TOUT,E15.6)
C
C  CALL THE INTEGRATOR
C
      CALL DRIVEP (NODE,T0,H,U1,TOUT,EPS,MF,INDEX,WORK,IWORK,X,Y)
C
C  CHECK ERROR RETURN
C
      IF (INDEX .NE. 0) GO TO 1000
C
C  OUTPUT THE RESULTS
C
      WRITE (6,220) HUSED,NQUSED,NSTEP,NFE,NJE
  220 FORMAT (/7H HUSED ,E10.4,7H ORDER ,I3,7H NSTEP ,I6,5H NFE ,I5,
     *        5H NJE ,I5)
      WRITE (6,230) ((U1(IX,IY),IX=1,NX),IY=1,NY)
  230 FORMAT ( /9H U VALUES  //(3H    ,10E11.4))
C
C  DETERMINE THE ERROR
C
      ERMAX = 0.0
      DO 240 IY = 1,NY
         DO 240 IX = 1,NX
          ERR1(IX,IY) = U1(IX,IY) - REALN1(IX,IY)
          IF (ABS(ERMAX) .LT. ABS(ERR1(IX,IY))) ERMAX=ERR1(IX,IY)
  240     CONTINUE
      WRITE (6,250) ERMAX
  250 FORMAT (/51H MAXIMUM DIFFERENCE FROM THE STEADY STATE SOLUTION //
     *        (3H    ,E11.4))
        TOUT=TOUT*10.
  260 CONTINUE
C
C************************** END OF EXAMPLE 1. **************************
C
      GO TO 1000
 300  CONTINUE
C***********************************************************************
C
C  EXAMPLE 2.  BURGERS EQUATION
C
C***********************************************************************
      WRITE (6,310)
  310 FORMAT(1H1,3X,40HOUTPUT FROM EXAMPLE 2.  BURGERS EQUATION//)
C
C  DEFINE THE PROBLEM PARAMETERS.
C
      DL=.010
      DLI=1.0/(2.0*DL)
      NX=5
      NY=31
      NPDE=1
      NODE=NX*NY*NPDE
      MF=10
      INDEX=1
      T0=0.0
      H=0.1E-06
      EPS=0.1E-03
      DY = 1.0/(FLOAT(NY)-1.0)
      DX = DY
      DO 320 I=1,NX
  320   X(I)=FLOAT(I)*DX-DX
      DO 330 I=1,NY
  330   Y(I)=FLOAT(I)*DY-DY
      IWORK(1) = NPDE
      IWORK(2) = NX
      IWORK(3) = NY
      IWORK(4) = 12
      IWORK(5) = 2687
      IWORK(6) = 155
C
C  DEFINE THE INITIAL CONDITION
C
      DO 360 IY = 1,NY
        DO 360 IX = 1,NX
          DEV=EXP((-X(IX)-Y(IY)+T0)*DLI)
  360     U2(IX,IY)=DEV/(1.0+DEV)
      WRITE (6,370) NODE,T0,H,EPS,MF,((U2(IX,IY),IX=1,NX),IY=1,NY)
  370 FORMAT (6H NODE ,I3,4H T0 ,F9.2,3H H ,E8.1,5H EPS ,E8.1,4H MF ,I2
     *        //16H INITIAL VALUES / 31(/3H    , 5E11.4))
C
C  SET UP THE LOOP FOR CALLING THE INTEGRATOR AT DIFFERENT TOUT VALUES
C
      TOUT=0.0
      DO 460 I=1,4
      TOUT=TOUT+.10
      WRITE (6,400) TOUT
  400 FORMAT (// 5H TOUT,E15.6)
C
C  CALCULATE THE ACTUAL SOLUTION
C
      DO 410 IY = 1,NY
        DO 410 IX = 1,NX
          DEV=EXP((-X(IX)-Y(IY)+TOUT)*DLI)
  410     REALN2(IX,IY)=DEV/(1.0+DEV)
C
C  CALL THE INTEGRATOR
C
      CALL DRIVEP (NODE,T0,H,U2,TOUT,EPS,MF,INDEX,WORK,IWORK,X,Y)
C
C  CHECK ERROR RETURN
C
      IF (INDEX .NE. 0) GO TO 1000
C
C  OUTPUT THE RESULTS
C
      WRITE (6,420) HUSED,NQUSED,NSTEP,NFE,NJE
  420 FORMAT (/7H HUSED ,E10.4,7H ORDER ,I3,7H NSTEP ,I6,5H NFE ,I5,
     *        5H NJE ,I5)
      WRITE (6,430) ((U2(IX,IY),IX=1,NX),IY=1,NY)
  430 FORMAT ( /9H U VALUES  // (3H    ,5E11.4))
C
C  DETERMINE THE ERROR
C
      ERMAX=0.0
      DO 440 IY = 1,NY
         DO 440 IX = 1,NX
          ERR2(IX,IY)=U2(IX,IY)-REALN2(IX,IY)
          IF (ABS(ERMAX) .LT. ABS(ERR2(IX,IY))) ERMAX=ERR2(IX,IY)
  440     CONTINUE
      WRITE (6,450) ERMAX
  450 FORMAT (/13H MAX ERROR U ,E11.4)
      WRITE (6,455) ((ERR2(IX,IY),IX=1,NX),IY=1,NY)
  455 FORMAT (/6H ERROR // (3H    ,5E11.4))
  460 CONTINUE
C
C************************** END OF EXAMPLE 2. **************************
C
      GO TO 1000
 500  CONTINUE
C***********************************************************************
C
C  EXAMPLE 3.  COUPLED SYSTEM OF PDE*S
C
C***********************************************************************
      WRITE (6,510)
  510 FORMAT(1H1,3X,47HOUTPUT FROM EXAMPLE 3.  COUPLED SYSTEM OF PDE*S//
     *)
C
C  DEFINE THE PROBLEM PARAMETERS.  NPDE,NX, AND NY PRIMARILY DETERMINE
C  THE DIMENSIONS FOR THE ARRAYS IN PDETWO AND THE MODIFIED GEARB.
C
      NX=5
      NY=5
      NPDE=2
      NODE=NX*NY*NPDE
      MF=22
      INDEX=1
      T0=0.0
      H=0.1E-06
      EPS=0.1E-03
      DX=1.0/(FLOAT(NX)-1.0)
      DO 520 IX=1,NX
        X(IX)=FLOAT(IX)*DX-DX
  520   Y(IX)=X(IX)
      IWORK(1) = NPDE
      IWORK(2) = NX
      IWORK(3) = NY
      IWORK(4) = 5
      IWORK(5) = 2308
      IWORK(6) = 50
C
C  DEFINE THE INITIAL CONDITIONS
C
      DO 560 IY=1,NY
        DO 560 IX=1,NX
          U3(1,IX,IY)=X(IX)+Y(IY)
  560     U3(2,IX,IY)=2.0*X(IX)+3.0*Y(IY)
      WRITE (6,570) NODE,T0,H,EPS,MF,((U3(1,IX,IY),IX=1,NX),IY=1,NY),
     *             ((U3(2,IX,IY),IX=1,NX),IY=1,NY)
  570 FORMAT (6H NODE ,I3,4H T0 ,F9.2,3H H ,E8.1,5H EPS ,E8.1,4H MF ,I2
     *        //19H INITIAL U1 VALUES  / 5(/3H    ,5E11.4)
     *        //19H INITIAL U2 VALUES  / 5(/3H    ,5E11.4))
C
C  SET UP THE LOOP FOR CALLING THE INTEGRATOR AT DIFFERENT TOUT VALUES
C
      TOUT=.25
      DO 660 I=1,4
      WRITE (6,600) TOUT
  600 FORMAT (//5H TOUT ,E15.6)
C
C  CALCULATE THE ACTUAL SOLUTION
C
      DO 610 IY=1,NY
        DO 610 IX=1,NX
          REALN3(2,IX,IY)=TOUT+2.0*X(IX)+3.0*Y(IY)
  610     REALN3(1,IX,IY)=TOUT +X(IX)+Y(IY)
C
C  CALL THE INTEGRATOR
C
      CALL DRIVEP (NODE,T0,H,U3,TOUT,EPS,MF,INDEX,WORK,IWORK,X,Y)
C
C  CHECK ERROR RETURN
C
      IF (INDEX .NE. 0) GO TO 1000
C
C  OUTPUT THE RESULTS
C
      WRITE (6,620) HUSED,NQUSED,NSTEP,NFE,NJE
  620 FORMAT (/7H HUSED ,E10.4,7H ORDER ,I3,7H NSTEP ,I6,5H NFE ,I5,
     *        5H NJE ,I5)
      WRITE (6,630) ((U3(1,IX,IY),IX=1,NX),IY=1,NY)
  630 FORMAT (/11H U1 VALUES //(3H    ,5E11.4))
      WRITE (6,635) ((U3(2,IX,IY),IX=1,NX),IY=1,NY)
  635 FORMAT (/11H U2 VALUES //(3H    ,5E11.4))
C
C  DETERMINE THE ERROR
C
      DO 640 IY=1,NY
         DO 640 IX=1,NX
          ERR3(1,IX,IY)=U3(1,IX,IY)-REALN3(1,IX,IY)
          ERR3(2,IX,IY)=U3(2,IX,IY)-REALN3(2,IX,IY)
  640     CONTINUE
      WRITE (6,650) ((ERR3(1,IX,IY),IX=1,NX),IY=1,NY)
  650 FORMAT (/9H ERROR U1 //(3H    ,5E11.4))
      WRITE (6,655) ((ERR3(2,IX,IY),IX=1,NX),IY=1,NY)
  655 FORMAT (/9H ERROR U2 //(3H    ,5E11.4))
        TOUT=TOUT+.25
  660 CONTINUE
C
C************************** END OF EXAMPLE 3. **************************
C
 1000 CONTINUE
      STOP
      END
      SUBROUTINE BNDRYV (T,X,Y,U,AV,BV,CV,NPDE)                          BNDRYV
C
C  DEFINE THE VERTICAL BOUNDARY CONDITIONS
C
      REAL  T,U,X,Y,BV,AV,CV
      INTEGER NPDE
      COMMON /PROB/ DL,DLI,KODE
      DIMENSION U(NPDE),AV(NPDE),BV(NPDE),CV(NPDE)
      IF(KODE.EQ.1)  GO TO 100
      IF(KODE.EQ.2)  GO TO 200
      IF(KODE.EQ.3)  GO TO 300
C
C  EXAMPLE 1.  ELLIPTIC PDE
C
  100 CONTINUE
      AV(1) = 1.0
      BV(1) = 0.0
      CV(1) = 0.0
      GO TO 400
C
C  EXAMPLE 2.  BURGERS EQUATION
C
  200 CONTINUE
      AV(1) = 1.0
      BV(1) = 0.0
      DEV = EXP((-X-Y+T)*DLI)
      CV(1) = DEV/(1.0+DEV)
      GO TO 400
C
C  EXAMPLE 3.  COUPLED SYSTEM OF PDE*S
C
  300 CONTINUE
      IF (X .NE. 0.0) GO TO 310
      AV(1)=1.0
      BV(1)=0.0
      CV(1)=T+Y
      AV(2)=1.0
      BV(2)=0.0
      CV(2)=T+3.0*Y
      GO TO 400
  310 CONTINUE
      AV(1)=1.0
      BV(1)=0.0
      CV(1)=T+1.0+Y
      AV(2)=0.0
      BV(2)=1.0
      CV(2)=2.0
  400 CONTINUE
      RETURN
      END
      SUBROUTINE BNDRYH (T,X,Y,U,AH,BH,CH,NPDE)                          BNDRYH
C
C  DEFINE THE HORIZONTAL BOUNDARY CONDITIONS
C
      REAL  T,U,X,Y,BH,AH,CH
      INTEGER NPDE
      COMMON /PROB/ DL,DLI,KODE
      DIMENSION U(NPDE),AH(NPDE),BH(NPDE),CH(NPDE)
      IF(KODE.EQ.1)  GO TO 100
      IF(KODE.EQ.2)  GO TO 200
      IF(KODE.EQ.3)  GO TO 300
C
C  EXAMPLE 1.  ELLIPTIC PDE
C
  100 CONTINUE
      AH(1) = 1.0
      BH(1) = 0.0
      CH(1) = 0.0
      GO TO 400
C
C  EXAMPLE 2.  BURGERS EQUATION
C
  200 CONTINUE
      AH(1) = 1.0
      BH(1) = 0.0
      DEV = EXP((-X-Y+T)*DLI)
      CH(1) = DEV/(1.0+DEV)
      GO TO 400
C
C  EXAMPLE 3.  COUPLED SYSTEM OF PDE*S
C
  300 CONTINUE
      IF (Y .NE. 0.0) GO TO 310
      AH(1)=0.0
      BH(1)=1.0
      CH(1)=1.0
      AH(2)=1.0
      BH(2)=1.0
      CH(2)=T+2.0*X+3.0
      GO TO 400
  310 CONTINUE
      AH(1)=0.0
      BH(1)=U(2)
      CH(1)=T+2.0*X+3.0
      AH(2)=1.0
      BH(2)=0.0
      CH(2)=T+2.0*X+3.0
  400 CONTINUE
      RETURN
      END
      SUBROUTINE DIFFH (T,X,Y,U,DH,NPDE)                                 DIFFH
C
C  DEFINE THE HORIZONTAL DIFFUSION COEFFICIENTS
C
      REAL  T,U,X,Y,DH
      INTEGER NPDE
      COMMON /PROB/ DL,DLI,KODE
      DIMENSION U(NPDE),DH(NPDE,NPDE)
      IF(KODE.EQ.1)  GO TO 100
      IF(KODE.EQ.2)  GO TO 200
      IF(KODE.EQ.3)  GO TO 300
C
C  EXAMPLE 1.  ELLIPTIC PDE
C
  100 CONTINUE
      DH(1,1) = 1.0
      GO TO 400
C
C  EXAMPLE 2.  BURGERS EQUATION
C
  200 CONTINUE
      DH(1,1) = 1.0
      GO TO 400
C
C  EXAMPLE 3.  COUPLED SYSTEM OF PDE*S
C
  300 CONTINUE
      DH(1,1)=0.0
      DH(1,2)=0.0
      DH(2,1)=U(2)
      DH(2,2)=U(1)
  400 CONTINUE
      RETURN
      END
      SUBROUTINE DIFFV (T,X,Y,U,DV,NPDE)                                 DIFFV
C
C  DEFINE THE VERTICAL DIFFUSION COEFFICIENTS
C
      REAL  T,U,X,Y,DV
      INTEGER NPDE
      COMMON /PROB/ DL,DLI,KODE
      DIMENSION U(NPDE),DV(NPDE,NPDE)
      IF(KODE.EQ.1)  GO TO 100
      IF(KODE.EQ.2)  GO TO 200
      IF(KODE.EQ.3)  GO TO 300
C
C  EXAMPLE 1.  ELLIPTIC PDE
C
  100 CONTINUE
      DV(1,1) = 1.0
      GO TO 400
C
C  EXAMPLE 2.  BURGERS EQUATION
C
  200 CONTINUE
      DV(1,1) = 1.0
      GO TO 400
C
C  EXAMPLE 3.  COUPLED SYSTEM OF PDE*S
C
  300 CONTINUE
      DV(1,1)=U(2)
      DV(1,2)=U(1)
      DV(2,1)=0.0
      DV(2,2)=0.0
  400 CONTINUE
      RETURN
      END
      SUBROUTINE F(T,X,Y,U,UX,UY,DUXX,DUYY,DUDT,NPDE)                    F
C
C  DEFINE THE PDE
C
      REAL  T,U,X,Y,UX,UY,DUXX,DUYY,DUDT,EXP
      INTEGER NPDE
      COMMON /PROB/ DL,DLI,KODE
      DIMENSION U(NPDE),UX(NPDE),UY(NPDE),DUXX(NPDE,NPDE),
     * DUYY(NPDE,NPDE),DUDT(NPDE)
      IF(KODE.EQ.1)  GO TO 100
      IF(KODE.EQ.2)  GO TO 200
      IF(KODE.EQ.3)  GO TO 300
C
C  EXAMPLE 1.  ELLIPTIC PDE
C
  100 CONTINUE
      DUDT(1) = DUXX(1,1) + DUYY(1,1) - 6.0*X*Y*EXP(X)*EXP(Y)*
     * (X*Y+X+Y-3.0)
      GO TO 400
C
C  EXAMPLE 2.  BURGERS EQUATION
C
  200 CONTINUE
      DUDT(1) = -U(1)*(UX(1)+UY(1))+DL*(DUXX(1,1)+DUYY(1,1))
      GO TO 400
C
C  EXAMPLE 3.  COUPLED SYSTEM OF PDE*S
C
  300 CONTINUE
      DUDT(1)=U(2)*(DUYY(1,1)+DUYY(1,2))-3.0*U(2)*UX(2)+4.0*UY(1)-UY(2)
      DUDT(2)=DUXX(2,1)+DUXX(2,2)-UY(2)
  400 CONTINUE
      RETURN
      END
      SUBROUTINE PSETM (NPDE,NX,NY,X,Y,U,UMAX,USAVE,DUDTR,DUDT,PW,CON,   PSETM
     *                  MITER,IER,NEBAND,WORK,IWORK)
C
C***********************************************************************
C
C       PSETM IS INTENDED TO BE USED IN CONJUNCTION WITH THE
C  INTERFACE PDETWO FOR THE SOLUTION OF SECOND ORDER TIME DEP-
C  ENDENT PARTIAL DIFFERENTIAL EQUATIONS (PDE*S) DEFINED OVER A
C  RECTANGULAR REGION.  PSETM IS SPECIFICALLY DESIGNED TO REPLACE
C  THE ROUTINE PSETB USED IN THE ORDINARY DIFFERENTIAL EQUATION
C  INTEGRATOR GEARB (1) IN ORDER TO MINIMIZE THE COMPUTATIONS
C  REQUIRED TO GENERATE THE JACOBIAN MATRIX.
C
C       PSETM IS CALLED BY STIFFP, A ROUTINE IN THE MODIFIED GEARB,
C  WHEN THE INTEGRATION METHOD (MITER=1 OR 2) REQUIRES A JACOBIAN
C  MATRIX.  WHEN MITER=1, THE JACOBIAN IS DETERMINED BY A USER
C  DEFINED ROUTINE, PDB.  THIS METHOD SHOULD BE AVOIDED SINCE IT RE-
C  QUIRES THE USER TO COMPLETELY UNDERSTAND PDETWO.  IF THIS
C  METHOD IS USED, PSETM OFFERS NO ADVANTAGE OVER PSETB.  IF
C  MITER=2, THE JACOBIAN IS APPROXIMATED BY FINITE DIFFERENCES,
C  USING THE APPROXIMATION  J = (DUDT(U+R) - DUDT(U)) / R,
C  WHERE DUDT(U+R) AND DUDT(U) ARE THE VALUES OF THE RIGHT HAND
C  SIDE OF THE PDE*S EVALUATED AT U+R AND U RESPECTIVELY, AND R
C  IS SOME SMALL NUMBER.  PSETM TAKES ADVANTAGE OF THE STRUCTURE
C  OF THE JACOBIAN MATRIX REQUIRED BY PDETWO TO REDUCE THE COMPU-
C  TATIONS NEEDED TO GENERATE THE MATRIX (REQUIRING ONLY 5 * NPDE
C  CALLS TO PDETWO).  WITH EACH CALL TO PDETWO, PSETM DETERMINES
C  THE ENTRIES IN THE JACOBIAN FOR THE MESH POINTS IN THE FOLLOWING
C  PATTERN ..
C
C                  X O O O O X O O O O
C                  O O O X O O O O X O
C                  O X O O O O X O O O
C                  O O O O X O O O O X
C                  O O X O O O O X O O
C                  X O O O O X O O O O
C
C  WHERE X REPRESENTS THE POINTS FOR WHICH THE JACOBIAN IS
C  APPROXIMATED.  MITER=2 IS THE RECOMMENDED METHOD.
C
C       TO CREATE A DOUBLE PRECISION VERSION OF PSETM..
C  CHANGE THE REAL STATEMENT BELOW AND CHANGE THE SINGLE PRECISION
C  FUNCTIONS ABS, SQRT, AND AMAX1.
C
C***********************************************************************
C  THE  PARAMETERS
C***********************************************************************
C
C  NPDE         IS THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS.
C  NX           IS THE NUMBER OF MESH POINTS IN THE HORIZONTAL
C               DIRECTION.
C  NY           IS THE NUMBER OF MESH POINTS IN THE VERTICAL
C               DIRECTION.
C  X            ARE THE MESH POINTS IN THE HORIZONTAL DIRECTION.
C  Y            ARE THE MESH POINTS IN THE VERTICAL DIRECTION.
C  U            CONTAINS THE CURRENT SOLUTIONS OF THE PARTIAL
C               DIFFERENTIAL EQUATIONS.
C  UMAX         IS THE MAXIMUM U VALUES, WHICH ARE USED TO SCALE
C               THE VALUE IN R AND FOR ERROR CONTROL IN GEARB.
C  USAVE        IS A TEMPORY STORE FOR U VALUES.
C  DUDTR        IS THE VALUE OF DUDT(U+R) RETURNED FROM PDETWO.
C  DUDT         IS THE VALUE OF DUDT(U), WHICH IS PASSED TO PSETM
C               FROM STIFFP.
C  PW           STORES THE APPROXIMATION FOR THE JACOBIAN.
C  CON          IS THE CONSTANT (-H*EL(1)).
C  MITER        INDICATES THE TYPE OF ITERATION METHOD BEING USED
C               BY THE INTEGRATOR.
C               MITER = 0 MEANS FUNCTIONAL ITERATION.
C               MITER = 1 MEANS THE CHORD METHOD WITH AN ANALYTIC
C                         JACOBIAN SUPPLIED IN THE USER DEFINED
C                         ROUTINE PDB.  THIS METHOD IS IN GEARB,
C                         BUT IT SHOULD BE AVOIDED.
C               MITER = 2 MEANS THE CHORD METHOD WITH THE JACO-
C                         BIAN CALCULATED IN PSETM.  THIS IS THE
C                         ONLY METHOD WHERE PSETM IS USEFUL.
C               MITER = 3 MEANS THE CHORD METHOD WITH THE JACO-
C                         BIAN REPLACED BY A DIAGONAL APPROXI-
C                         MATION BASED ON A DIRECTIONAL DERIV-
C                         ATIVE.  THIS METHOD IS IN GEARB, BUT
C                         IS NOT GENERALLY USEFUL IN SOLVING PDE*S.
C  IER          IS AN ERROR INDICATOR USED IN THE ROUTINE DECBR.
C  NEBAND       = 3 * ML + 1
C  WORK         PROVIDES WORKING STORAGE FOR DYNAMIC DIMENSIONING
C               OF THE ARRAYS IN PDETWO AND THE MODIFIED GEARB.
C  IWORK        IS THE PIVOT VECTOR USED IN THE LU DECOMPOSITION.
C
C***********************************************************************
C  THE VARIABLES
C***********************************************************************
C
C  R            IS A SMALL INCREMENT USED IN EVALUATING THE
C               FUNCTION NEAR THE CURRENT SOLUTION U ( I.E.
C               DUDT (U+R)).
C  R0           IS A LOWER BOUND ON THE SIZE OF R.
C  D            IS A SMALL NUMBER USED TO APPROXIMATE THE
C               DERIVATIVE.
C  NXNPDE       = NX * NPDE
C
C  T            IS THE TIME BEING USED FOR THE INTEGRATION.
C  H            IS THE CURRENT TIME STEP SIZE BEING USED IN THE
C               INTEGRATION.
C  UROUND       IS THE UNIT ROUNDOFF OF THE MACHINE.
C  EPSJ         IS SQRT(UROUND).
C  ML           IS THE WIDTH OF THE LOWER (AND UPPER) HALF OF THE
C               BAND OF THE JACOBIAN.  ML = (NX + 1) * NPDE - 1.
C  IW1 - IW27   ARE THE SUBSCRIPT POINTERS USED TO INDICATE THE
C               STORAGE LOCATIONS IN WORK OF ARRAYS IN PDETWO AND
C               THE MODIFIED GEARB.
C
C***********************************************************************
C  THE ROUTINES CALLED BY PSETM
C***********************************************************************
C
C  PDB(NODE,T,U,PW,NODE,ML,ML) IS A USER DEFINED ROUTINE WHICH DEFINES
C               THE JACOBIAN IN PW IF MITER=1.  SINCE THIS METHOD IS
C               NOT RECOMMENDED, A DUMMY ROUTINE IS PROVIDED.
C  DECBR(NEBAND,NODE,ML,ML,PW,IWORK,IER) COMPUTES THE LU DECOMPOSITION
C               OF THE JACOBIAN FOR THE DIRECT SOLUTION OF THE
C               JACOBIAN.
C  PDETWO(NPDE,NX,NY,X,Y,T,U,...) IS THE INTERFACE WHICH DESCRETIZES
C               THE SPATIAL VARIABLES OF THE TWO DIMENSIONAL SYSTEM OF
C               PDE*S, THEREBY EVALUATING THE RIGHT HAND SIDE OF
C               THE SYSTEM OF ODE*S.
C
C***********************************************************************
C  REFERENCES
C***********************************************************************
C
C       1.  A. C. HINDMARSH, GEARB.. SOLUTION OF ORDINARY DIFFERENTIAL
C             EQUATIONS HAVING BANDED JACOBIAN, UCID-30059 REV. 2,
C             LAWRENCE LIVERMORE LABORATORY, P.O.BOX 808, LIVERMORE,
C             CA 94550, JUNE 1977.
C
C***********************************************************************
C
      REAL  U,UMAX,USAVE,DUDTR,DUDT,PW,CON,R,R0,D,T,H,UROUND
     *                 ,EPSJ,DUMMY,WORK
      COMMON /GEAR1/ T,H,DUMMY(3),UROUND,IDUMMY(4)
      COMMON /GEAR2/ EPSJ,ML,TDUMY(5)
      COMMON /ISTORE/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11,
     *                IW12,IW13,IW14,IW15,IW16,IW17,IW18,IW19,IW20,IW21,
     *                IW22,IW23,IW24,IW25,IW26,IW27
      DIMENSION UMAX(NPDE,NX,NY),USAVE(NPDE,NX,NY),DUDTR(NPDE,NX,NY),
     *          DUDT(NPDE,NX,NY),PW(NEBAND,1),U(NPDE,NX,NY),
     *          WORK(1),IWORK(1),X(1),Y(1)
C
C  SET THE NUMBER OF ORDINARY DIFFERENTIAL EQUATIONS
C
      NODE = NPDE*NX*NY
C
C  IF MITER = 1, CALL PDB TO GENERATE THE JACOBIAN
C
      IF (MITER .EQ. 2) GO TO 20
      CALL PDB (NODE,T,U,PW,NODE,ML,ML)
        DO 10 J=1,NODE
      DO 10 I=1,NEBAND
   10     PW(I,J)=PW(I,J)*CON
      GO TO 150
C
C  IF MITER = 2, APPROXIMATE THE JACOBIAN WITH FINITE DIFFERENCES
C
   20 D = 0.0
      DO 30 K=1,NY
        DO 30 J=1,NX
          DO 30 I=1,NPDE
            D=D+DUDT(I,J,K)**2
   30       USAVE(I,J,K)=U(I,J,K)
      R0 = ABS(H)*SQRT(D)*1.E03*UROUND
      ISIGN=-1
      NXNPDE=NX*NPDE
      DO 50 I=1,NODE
        DO 50 J=1,NEBAND
   50     PW(J,I)=0.0
      DO 140 M=1,5
C
C  APPROXIMATE THE JACOBIAN FOR EACH PDE
C
          DO 140 IPDE=1,NPDE
C
C  DETERMINE THE MESH POINTS FOR WHICH THE APPROXIMATIONS ARE TO
C  BE CALCULATED
C
            DO 70 IY=1,NY
              IXSTRT=1+MOD(M+2*IY,5)
              IF(IXSTRT.GT.NX) GO TO 70
              DO 60 IX=IXSTRT,NX,5
                R=AMAX1(EPSJ*UMAX(IPDE,IX,IY),R0)
   60           U(IPDE,IX,IY)=U(IPDE,IX,IY)+R
   70         CONTINUE
C
C  DETERMINE DUDT(U+R)
C
            CALL PDETWO (NPDE,NX,NY,X,Y,T,U,DUDTR,WORK,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))
            DO 140 IYY=1,NY
              IXSTRT=1+MOD(M+2*IYY,5)
              IF(IXSTRT.GT.NX) GO TO 140
              DO 130 IXX=IXSTRT,NX,5
                R=AMAX1(EPSJ*UMAX(IPDE,IXX,IYY),R0)
                D=CON/R
                IY=IYY
                IX=IXX
                IBEG1=ML+IPDE+1
                IBEG2=(IX-1)*NPDE+(IY-1)*NXNPDE
                I1=IBEG1
                I2=IBEG2
C
C  WITH FIVE POINT CENTERED DIFFERENCING, THE ONLY NONZERO
C  ENTRIES IN THE JACOBIAN ARE GENERATED FROM THE FIVE MESH
C  POINTS USED IN THE DIFFERENCING.
C
                DO 120 K=1,5
C
C  IF K=1, DETERMINE THE JACOBIAN APPROXIMATIONS FOR THE MESH
C  POINT (X(IXX),Y(IYY))
C
                  GO TO (100,80,80,90,90),K
C
C
C  IF K=2 OR K=3, DETERMINE THE JACOBIAN APPROXIMATIONS FOR THE
C  MESH POINTS (X(IXX-1),Y(IYY)) AND (X(IXX+1),Y(IYY))
C
   80             IX=IXX-ISIGN
                  IF (IX .EQ. 0 .OR. IX .GT. NX) GO TO 120
                  I1=IBEG1+NPDE*ISIGN
                  I2=IBEG2-NPDE*ISIGN
                  GO TO 100
C
C  IF K=4 OR K=5, DETERMINE THE JACOBIAN APPROXIMATIONS FOR THE
C  MESH POINTS (X(IXX),Y(IYY-1)) AND (X(IXX),Y(IYY+1))
C
   90             IX=IXX
                  IY=IYY-ISIGN
                  IF (IY .EQ. 0 .OR. IY .GT. NY) GO TO 120
                  I1=IBEG1+NXNPDE*ISIGN
                  I2=IBEG2-NXNPDE*ISIGN
C
C  CALCULATE THE FINITE DIFFERENCE APPROXIMATIONS
C   (DUDT(U+R) - DUDT(U)) / R
C
  100             DO 110 J=1,NPDE
                    K1=I1-J
                    K2=I2+J
  110               PW(K1,K2)=(DUDTR(J,IX,IY)-DUDT(J,IX,IY))*D
  120             ISIGN=-ISIGN
                U(IPDE,IXX,IYY)=USAVE(IPDE,IXX,IYY)
  130           USAVE(IPDE,IXX,IYY)=0.0
  140         CONTINUE
  150 DO 160 I=1,NODE
  160   PW(ML+1,I)=PW(ML+1,I)+1.0
C
C  DO THE LU DECOMPOSITION ON THE JACOBIAN
C
      CALL DECBR (NEBAND,NODE,ML,ML,PW,IWORK,IER)
      RETURN
      END
      SUBROUTINE PDETWO (NPDE,NX,NY,X,Y,T,U,DUDT,DVS,DVST,DV,DH,         PDETWO
     *                   AH,BH,CH,AV,BV,CV,UX,UY,DXI,DXIR,DXIC,
     *                   XAVG,UDVA,UDHR,UDVB,UDHL,UAVGH,UAVGV)
C
C***********************************************************************
C
C       PDETWO IS AN INTERFACE DESIGNED TO SOLVE A SYSTEM OF SECOND
C  ORDER TIME DEPENDENT PARTIAL DIFFERENTIAL EQUATIONS (PDE*S)
C  DEFINED OVER A RECTANGLE IN CONJUNCTION WITH AN ORDINARY
C  DIFFERENTIAL EQUATION (ODE) INTEGRATOR PACKAGE (PRIMARILY A MODI-
C  FIED VERSION OF GEARB(1)).  IT IS AN EXTENSION OF PDEONE (2,3), AN
C  INTERFACE DEVELOPED FOR THE SOLUTION OF ONE-DIMENSIONAL SYSTEMS OF
C  PDE*S.  THIS EXTENSION CONSISTS OF CHANGING FROM THREE POINT
C  DIFFERENCING TO FIVE POINT DIFFERENCING IN ORDER TO APPROXIMATE
C  THE SPATIAL DERIVATIVES IN THE TWO DIMENSIONAL PDE SYSTEM.  TO
C  SOLVE A PDE SYSTEM USING THIS SOFTWARE INTERFACE IN CONJUNCTION
C  WITH AN ODE INTEGRATOR, THE USER IS REQUIRED ONLY TO PROVIDE
C  THE SPATIAL MESH AND DEFINE THE PROBLEM TO BE EVALUATED.  THE
C  INTERFACE WILL THEN FORM AND EVALUATE A SEMIDISCRETE APPROXI-
C  MATION OF THIS SYSTEM OF PDE*S AND THUS PERMIT ONE TO SOLVE
C  THE ORIGINAL PDE SYSTEM AS A SYSTEM OF TIME DEPENDENT ODE*S
C  USING AN ODE INTEGRATOR.
C
C       TO CREATE A DOUBLE PRECISION VERSION OF PDETWO..
C  CHANGE THE REAL STATEMENT BELOW.
C
C***********************************************************************
C  PROBLEM DEFINITION
C***********************************************************************
C
C       PDETWO IS DESIGNED TO BE USED TO SOLVE A COUPLED SYSTEM OF
C  NPDE PDE*S.  THE L-TH PDE OF THIS SYSTEM IS OF THE FORM ..
C
C     DUDT(L) = F (T, X, Y, U, UX, UY, DUXX, DUYY )
C
C  WHERE
C
C       DUDT(L) REPRESENTS THE FIRST PARTIAL OF THE L-TH COMPONENT
C               OF U WITH RESPECT TO T,
C       T       IS THE CURRENT TIME,
C       X,Y     DEFINE THE POSITION IN THE HORIZONTAL AND VERTICAL
C               DIRECTIONS RESPECTIVELY,
C
C       U    =  (  U(1), ... ,  U(NPDE) ),
C       UX   =  ( UX(1), ... , UX(NPDE) ),
C       UY   =  ( UY(1), ... , UY(NPDE) ),
C
C       DUXX    IS AN NPDE BY NPDE ARRAY SUCH THAT
C                           D
C               DUXX(L,K) = -- (DH(L,K)*UX(K)),
C                           DX
C               THE L-TH ROW OF DUXX CORRESPONDS TO THE L-TH PDE,
C
C       DUYY    IS AN NPDE BY NPDE ARRAY SUCH THAT
C                           D
C               DUYY(L,K) = -- (DV(L,K)*UY(K)),
C                           DY
C               THE L-TH ROW OF DUYY CORRESPONDS TO THE L-TH PDE,
C
C  AND
C
C       U(K)    IS THE VALUE OF THE SOLUTION FOR THE K-TH PDE AT
C               (T, X, Y),
C       UX(K)   IS THE FIRST PARTIAL DERIVATIVE OF U(K) WITH RESPECT
C               TO X,
C       UY(K)   IS THE FIRST PARTIAL DERIVATIVE OF U(K) WITH RESPECT
C               TO Y,
C       DH(L,K) IS THE DIFFUSION COEFFICIENT IN THE HORIZONTAL
C               DIRECTION FOR THE L-TH PDE ASSOCIATED WITH U(K),
C       DV(L,K) IS THE DIFFUSION COEFFICIENT IN THE VERTICAL
C               DIRECTION FOR THE L-TH PDE ASSOCIATED WITH U(K).
C
C       DH(L,K) AND DV(L,K) MAY BE FUNCTIONS OF T, X, Y AND U.
C
C  HORIZONTAL BOUNDARY CONDITIONS
C
C     AH(L)*U(L) + BH(L)*UY(L) = CH(L)
C
C  VERTICAL BOUNDARY CONDITIONS
C
C     AV(L)*U(L) + BV(L)*UX(L) = CV(L)
C
C       AH(L), BH(L) AND CH(L) (AV(L), BV(L) AND CV(L)) ARE AT LEAST
C     PIECEWISE CONTINUOUS FUNCTIONS OF THEIR RESPECTIVE VARIABLES.
C     IF BH(L) .NE. 0 (BV(L) .NE. 0) THEN AH(L), BH(L) AND CH(L)
C     (AV(L), BV(L) AND CV(L)) MAY BE FUNCTIONS OF T, X, Y AND U
C     BUT OTHERWISE THEY MAY ONLY BE FUNCTIONS OF T, X AND Y. NULL
C     BOUNDARY CONDITIONS ARE NOT ALLOWED SINCE PDETWO WAS NOT
C     DESIGNED TO SOLVE HYPERBOLIS PDE*S. A ZERO DIVIDE WILL OCCUR
C     IF A NULL BOUNDARY CONDITION IS SPECIFIED.
C
C  INITIAL CONDITIONS
C
C       THE INITIAL SOLUTION IS DEFINED FOR EACH U(K) TO BE A
C     KNOWN FUNCTION OF X AND Y FOR THE INITIAL TIME T = T0.
C     THE INITIAL CONDITIONS NEED NOT BE CONSISTENT WITH THE
C     BOUNDARY CONDITIONS AND MAY CONTAIN DISCONTINUITIES.
C
C***********************************************************************
C  USER SUPPLIED ROUTINES
C***********************************************************************
C
C       THE USER MUST PROVIDE A MAIN PROGRAM AND FIVE SUBROUTINES
C  BNDRYH, BNDRYV, DIFFV, DIFFH AND F.  THESE ROUTINES ARE ALL THAT
C  IS REQUIRED TO DEFINE COMPLETELY THE PROBLEM GIVEN ABOVE.
C
C  1)     THE MAIN PROGRAM DEFINES THE SPATIAL MESH, THE INITIAL
C     CONDITIONS AND THE PARAMETERS FOR THE ODE INTEGRATOR, CALLS
C     THE INTEGRATOR AND PRINTS OR PLOTS THE RESULTS.  THE MAIN
C     PROGRAM SHOULD BE CONSTRUCTED AS FOLLOWS ..
C
C     DIMENSION U(***,*,**),X(*),Y(**)
C     DIMENSION WORK(*****),IWORK(****)
C
C         FOR THE DIMENSIONS ABOVE ENTER THE ACTUAL NUMERICAL
C       VALUES FOR * = NX, ** = NY, *** = NPDE , **** = NODE, AND
C       ***** = NPDE * ( NPDE*(NX*2+3) + 13 + NX ) + NX*4 + 4*NODE
C       + LPW + LU
C       WHERE
C          NODE = NPDE*NX*NY
C          LPW  = 1                       IF MITER = 0
C               = (3*(NX+1)*NPDE-2)*NODE  IF MITER = 1,2
C               = NODE                    IF MITER = 3
C          LU   = NODE*(MORDER+1)
C       SEE MF AND MORDER BELOW FOR THE DEFINITION OF MITER AND
C       MORDER.
C
C       DEFINE ..
C
C       1) THE NUMBER OF MESH POINTS IN THE X (HORIZONTAL) DIRECTION,
C          NX .GE. 3 AND THE NUMBER OF MESH POINTS IN THE Y (VERTICAL)
C          DIRECTION, NY .GE. 3  ( TO CONSERVE STORAGE ORIENT THE
C          PROBLEM SO THAT NX .LE. NY ),
C       2) THE MESH POINTS (I.E. X(1) .LT. X(2) .LT. ... .LT. X(NX)
C          AND Y(1) .LT. Y(2) .LT. ... .LT. Y(NY)),
C       3) NPDE, THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS,
C       4) INITIAL VALUES FOR ALL OF U,
C       5) MORDER, THE MAXIMUM ORDER OF THE METHOD USED IN THE
C          MODIFIED GEARB.  MORDER MUST BE LESS THAN OR EQUAL TO 12
C          IF METH = 1 OR IF METH = 2, IT MUST BE LESS THAN OR EQUAL
C          TO 5.  SEE MF BELOW FOR THE DEFINITION OF METH.
C       6) THE PARAMETERS FOR THE CALL TO THE INTEGRATOR, AND
C       7) THE FIRST SIX LOCATIONS OF THE ARRAY IWORK AS FOLLOWS..
C          IWORK(1) = NPDE
C          IWORK(2) = NX
C          IWORK(3) = NY
C          IWORK(4) = MORDER
C          IWORK(5) = NRWK
C          IWORK(6) = NRIWK
C          WHERE NRWK AND NRIWK ARE THE LENGTHS OF THE ARRAYS WORK
C          AND IWORK RESPECTIVELY.  THEY SHOULD BE AT LEAST EQUAL
C          TO ***** AND **** AS DEFINED IN THE DESCRIPTION OF THE
C          DIMENSIONS.
C
C         IN DETERMINING THE MESH SPACING, THE USER SHOULD BE AWARE
C       THAT THE DIFFERENCE APPROXIMATIONS AT THE INTERIOR POINTS
C       ARE SECOND ORDER ONLY FOR UNIFORM MESHES.  FOR NONUNIFORM
C       MESHES, ONLY FIRST ORDER APPROXIMATIONS SHOULD BE EXPECTED.
C       IN UNUSUAL SITUATIONS (SEE REMARK IN SECTION 3 OF THE
C       ACCOMPANYING PAPER) INVOLVING COUPLED SYSTEMS OF EQUATIONS,
C       CERTAIN APPROXIMATIONS AT THE BOUNDARY ARE ONLY FIRST ORDER.
C       IN THESE UNUSUAL SITUATIONS THE ERROR CAN BE MINIMIZED BY
C       CHOOSING SMALL MESH SPACINGS NEXT TO THE BOUNDARY.  IN ANY
C       CASE THE USER IS ADVISED TO CHOOSE A MESH WHICH IS LOCALLY
C       NEARLY UNIFORM.
C
C         FOR THE MODIFIED GEARB THE PARAMETERS INCLUDE THE DESIRED
C       OUTPUT TIME (TOUT), THE DESIRED LOCAL ACCURACY (EPS), THE
C       NUMBER OF ODE*S (NODE = NPDE*NX*NY), THE INITIAL TIME (T0),
C       THE INITIAL TIME STEP SIZE (H),THE TYPE OF CALL BEING MADE
C       (INDEX) AND THE TYPE OF INTEGRATION METHOD DESIRED (MF).
C       MF HAS TWO DECIMAL DIGITS, METH AND MITER (MF = 10*METH +
C       MITER).  METH IS THE BASIC METHOD INDICATOR..
C         METH  = 1 MEANS THE ADAMS METHODS.
C         METH  = 2 MEANS THE BACKWARD DIFFERENTIATION FORMULAS
C                   (BDF), OR STIFF METHODS OF GEAR.
C       MITER IS THE ITERATION METHOD INDICATOR..
C         MITER = 0 MEANS FUNCTIONAL ITERATION (NO PARTIAL
C                   DERIVATIVES NEEDED).
C         MITER = 1 MEANS THE CHORD METHOD WITH AN ANALYTIC
C                   JACOBIAN SUPPLIED IN THE USER DEFINED
C                   ROUTINE PDB.  THIS METHOD IS IN GEARB,
C                   BUT IT SHOULD BE AVOIDED.
C         MITER = 2 MEANS THE CHORD METHOD WITH THE JACO-
C                   BIAN CALCULATED IN PSETM.  THIS IS THE
C                   ONLY METHOD WHERE PSETM IS USEFUL.
C         MITER = 3 MEANS THE CHORD METHOD WITH THE JACO-
C                   BIAN REPLACED BY A DIAGONAL APPROXI-
C                   MATION BASED ON A DIRECTIONAL DERIV-
C                   ATIVE.  THIS METHOD IS IN GEARB, BUT
C                   IS NOT GENERALLY USEFUL IN SOLVING PDE*S.
C       SEE COMMENTS IN DRIVEP FOR ADDITIONAL INFORMATION ON THESE
C       PARAMETERS.  NATURALLY THESE PARAMETERS ARE DEPENDENT ON THE
C       INTEGRATOR BEING USED.  FINALLY THE USER SHOULD CALL THE
C       INTEGRATOR.
C
C     CALL DRIVEP (NODE,T0,H,U,TOUT,EPS,MF,INDEX,WORK,IWORK,X,Y)
C
C       THE INTEGRATOR WILL RETURN THE SOLUTIONS (U) AT T = TOUT
C       TO THE MAIN PROGRAM TO BE PRINTED OR PLOTTED.  IF A CONTINUA-
C       TION TO ANOTHER TOUT IS DESIRED, SIMPLY RESET TOUT AND CALL
C       DRIVEP AGAIN.
C
C     STOP
C     END
C
C  2)     THE BOUNDARY SUBROUTINES, BNDRYH AND BNDRYV PROVIDE
C     PDETWO WITH THE COEFFICIENTS FOR THE HORIZONTAL AND
C     VERTICAL BOUNDARY CONDITIONS RESPECTIVELY.  THE CONSTRUCTION
C     OF BNDRYV IS ANALOGOUS TO THE CONSTRUCTION OF BNDRYH.  THE
C     USER SHOULD CONSTRUCT BNDRYH AS FOLLOWS ..
C
C     SUBROUTINE BNDRYH (T,X,Y,U,AH,BH,CH,NPDE)
C     DIMENSION U(NPDE), AH(NPDE), BH(NPDE), CH(NPDE)
C
C         THE INCOMING PARAMETERS X AND Y REPRESENT ANY POINT
C       OF X(J) (J = 1,2,...,NX) AND EITHER Y(1) OR Y(NY) RESPECTIVELY.
C       DEFINE THE FUNCTIONS AH(K), BH(K) AND CH(K) (K = 1,2,...,NPDE)
C       FOR THE LOWER (Y=Y(1)) AND UPPER (Y=Y(NY)) BOUNDARIES. NULL
C       BOUNDARY CONDITIONS ARE NOT ALLOWED AND IF SPECIFIED WILL
C       RESULT IN A ZERO DIVIDE.
C
C         TO INSURE A COMPLETELY ACCURATE DEFINITION OF THE BOUNDARY
C       CONDITIONS AT THE CORNER MESH POINTS, THE USER SHOULD
C       CONSULT THE APPROPRIATE DIFFERENCE APPROXIMATIONS (SEE
C       SECTION 3 OF THE ACCOMPANYING PAPER).
C
C     RETURN
C     END
C
C  3)     DIFFH AND DIFFV DEFINE FOR PDETWO THE HORIZONTAL AND VERTICAL
C     DIFFUSION COEFFICIENTS RESPECTIVELY.  THE CONSTRUCTION OF
C     DIFFV IS ANALOGOUS TO THE CONSTRUCTION OF DIFFH.  THE
C     USER SHOULD CONSTRUCT DIFFH AS FOLLOWS ..
C
C     SUBROUTINE DIFFH (T,X,Y,U,DH,NPDE)
C     DIMENSION U(NPDE),DH(NPDE,NPDE)
C
C         IN THIS ROUTINE DEFINE THE DH(L,K) COEFFICIENTS (L,K =
C       1,2,...,NPDE).  THE INCOMING PARAMTERS X AND Y DENOTE EITHER
C       A BOUNDARY POINT OR A MESH MIDPOINT.
C
C     RETURN
C     END
C
C  4)     THE RIGHT HAND SIDE OF THE PDE IS DEFINED FOR PDETWO IN
C     THE SUBROUTINE F.  THIS ROUTINE IS CONSTRUCTED AS FOLLOWS..
C
C     SUBROUTINE F (T,X,Y,U,UX,UY,DUXX,DUYY,DUDT,NPDE)
C     DIMENSION U(NPDE), UX(NPDE), UY(NPDE), DUXX(NPDE,NPDE),
C    * DUYY(NPDE,NPDE), DUDT(NPDE)
C
C         IN THIS ROUTINE, THE INCOMING VALUES X AND Y REPRESENT
C       THE MESH POINT BEING EVALUATED AND UX, UY, DUXX AND DUYY
C       ARE THE VALUES DENOTED IN THE PROBLEM DEFINITION ABOVE.
C       USING THESE VALUES, DEFINE IN DUDT(L) (L = 1,2,...,NPDE)
C       THE RIGHT HAND SIDE OF THE PDE*S.
C
C     RETURN
C     END
C
C***********************************************************************
C  THE INPUT PARAMETERS
C***********************************************************************
C
C  NPDE         IS THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS.
C  NX           IS THE NUMBER OF MESH POINTS IN THE HORIZONTAL
C               DIRECTION.
C  NY           IS THE NUMBER OF MESH POINTS IN THE VERTICAL
C               DIRECTION.
C  X            ARE THE MESH POINTS IN THE HORIZONTAL DIRECTION.
C  Y            ARE THE MESH POINTS IN THE VERTICAL DIRECTION.
C  T            IS THE CURRENT TIME.
C  U            CONTAINS THE CURRENT SOLUTION VALUES FOR ALL THE
C               MESH POINTS.
C  DVS          SAVES THE VERTICAL DIFFUSION COEFFICIENTS FOR
C               FUTURE EVALUATIONS.
C  DVST         SAVES THE VERTICAL DIFFUSION COEFFICIENTS FOR THE
C               FUTURE EVALUATIONS OF THE TOP ROW OF MESH POINTS.
C  DH           RETURNS FROM DIFFH THE HORIZONTAL DIFFUSION COEFFI-
C               CIENTS AND CONTAINS DUXX, THE APPROXIMATIONS TO THE
C               HORIZONTAL DIVERGENCE TERMS, ON CALLS TO F AND STORES
C               THE PREVIOUS VALUE OF DUXX.
C  DV           RETURNS FROM DIFFV THE VERTICAL DIFFUSION COEFFI-
C               CIENTS AND CONTAINS DUYY, THE APPROXIMATIONS TO THE
C               VERTICAL DIVERGENCE TERMS, ON CALLS TO F.
C  AH, BH, CH   CONTAIN THE HORIZONTAL BOUNDARY COEFFICIENTS
C               PASSED TO PDETWO FROM BNDRYH.
C  AV, BV, CV   CONTAIN THE VERTICAL BOUNDARY COEFFICIENTS
C               PASSED TO PDETWO FROM BNDRYV.
C  UX,UY        STORE THE DIFFERENCE APPROXIMATIONS FOR THE FIRST
C               PARTIAL OF U WITH RESPECT TO X AND Y RESPECTIVELY.
C  DXI,DXIR,
C    DXIC       STORE INFORMATION ABOUT THE HORIZONTAL MESH SPACING.
C  XAVG         STORES INFORMATION ABOUT THE AVERAGE BETWEEN TWO MESH
C               POINTS ON THE HORIZONTAL AXIS.
C  UDVA,UDVB,
C    UDHR,UDHL  STORE INFORMATION FOR APPROXIMATING DUXX AND DUYY.
C  UAVGH,UAVGV  CONTAIN U AVERAGES FOR APPROXIMATING THE DIFFUSION
C               COEFFICIENTS AT THE MESH MID-POINTS.
C
C***********************************************************************
C  THE OUTPUT PARAMETERS
C***********************************************************************
C
C  U            CONTAINS SOLUTIONS UPDATED TO TIME T FOR THE
C               BOUNDARY MESH POINTS DEFINED BY DIRICHLET BOUNDARY
C               CONDITIONS.
C  DUDT         IS THE RIGHT HAND SIDE OF THE SYSTEM OF ODE*S PASSED
C               TO THE INTEGRATOR.  THESE VALUES ARE CALCULATED
C               FROM THE CENTERED DIFFERENCE APPROXIMATIONS OF
C               THE SPATIAL VARIABLES.
C
C***********************************************************************
C  REFERENCES
C***********************************************************************
C
C       1.  A. C. HINDMARSH, GEARB.. SOLUTION OF ORDINARY DIFFERENTIAL
C             EQUATIONS HAVING BANDED JACOBIAN, UCID-30059 REV. 2,
C             LAWRENCE LIVERMORE LABORATORY, P.O.BOX 808, LIVERMORE,
C             CA 94550, JUNE 1977.
C
C       2.  R.F. SINCOVEC AND N.K. MADSEN, SOFTWARE FOR NONLINEAR
C             PARTIAL DIFFERENTIAL EQUATIONS, ACM TRANSACTIONS ON
C             MATHEMATICAL SOFTWARE, SEPT. 1975, PP. 232-260.
C
C       3.  R.F. SINCOVEC AND N.K. MADSEN, ALGORITHM 494 PDEONE,
C             SOLUTIONS OF SYSTEMS OF PARTIAL DIFFERENTIAL EQUATIONS,
C             ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, SEPT.
C             1975, PP. 262-263.
C
C***********************************************************************
C
      REAL  T,U,X,Y,UX,UY,DXI,DYI,DXIC,DXIR,DYIA,DYIB,DYIC,
     * UDHL,UDHR,DUDT,UDVA,UDVB,XAVG,YAVG,BH,BV,DH,DV,
     * UAVGH,UAVGV,AH,AV,DVS,DVST,CH,CV,DX,DUXR,DUXL,DY,DUY
      DIMENSION DX(4),DUXR(2),DUXL(2),DY(4),DUY(2,2)
      DIMENSION DVS(NPDE,NPDE,NX), DVST(NPDE,NPDE,NX), DV(NPDE,NPDE),
     * DH(NPDE,NPDE,2), U(NPDE,NX,NY), DUDT(NPDE,NX,NY), AH(NPDE),
     * BH(NPDE), CH(NPDE), AV(NPDE), BV(NPDE), CV(NPDE), UX(NPDE),
     * UY(NPDE), DXI(NX), DXIR(NX), DXIC(NX), XAVG(NX), UDVA(NPDE,NX),
     * UDHR(NPDE), UDVB(NPDE), UDHL(NPDE), UAVGH(NPDE), UAVGV(NPDE),
     * X(NX), Y(NY)
C
C DETERMINE THE MESH SPACING ALONG THE HORIZONTAL AXIS
C
      DXI(1)=1./(X(2)-X(1))
      DXIR(1)=DXI(1)
      DXIC(1)=2.*DXI(1)
      XAVG(1)=.5*(X(2)+X(1))
      ILIM=NX-1
      ILIMY=NY-1
      DO 10 I=2,ILIM
        XAVG(I)=.5*(X(I+1)+X(I))
        DXI(I)=1./(X(I+1)-X(I-1))
        DXIR(I)=1./(X(I+1)-X(I))
   10   DXIC(I)=2.*DXI(I)
      DXI(NX)=DXIR(ILIM)
      DXIC(NX)=2.*DXI(NX)
C
C DETERMINE THE DIFFERENCING FOR THE BOUNDARIES
C
      DX(1)=X(2)-X(1)
      DX(2)=X(3)-X(2)
      DX(3)=X(NX-1)-X(NX-2)
      DX(4)=X(NX)-X(NX-1)
      DUXL(1)=(DX(1)+DX(2))/(DX(1)*DX(2))
      DUXL(2)=DX(1)/((DX(1)+DX(2))*DX(2))
      DUXR(1)=(DX(4)+DX(3))/(DX(4)*DX(3))
      DUXR(2)=DX(4)/((DX(4)+DX(3))*DX(3))
      DY(1)=Y(2)-Y(1)
      DY(2)=Y(3)-Y(2)
      DY(3)=Y(NY-1)-Y(NY-2)
      DY(4)=Y(NY)-Y(NY-1)
      DUY(2,1)=(DY(1)+DY(2))/(DY(1)*DY(2))
      DUY(2,2)=DY(1)/((DY(1)+DY(2))*DY(2))
      DUY(1,1)=(DY(4)+DY(3))/(DY(4)*DY(3))
      DUY(1,2)=DY(4)/((DY(4)+DY(3))*DY(3))
C
C THE FOLLOWING LOOP DETERMINES THE U APPROXIMATIONS FOR THE BOTTOM
C AND TOP BOUNDARIES.  IC HAS THE VALUE OF 1 FOR THE TOP AND 2 FOR THE
C BOTTOM
C
      DO 410 IC=1,2
        IF (IC .EQ. 1) GO TO 20
        M=1
        MN=2
        MNN=3
        SIGN=1.0
        GO TO 25
   20   M=NY
        MN=ILIMY
        MNN=MN-1
        SIGN=-1.0
   25   DYI=1./(Y(MN)-Y(M))*SIGN
        DYIA=DYI
        DYIC=2.*DYI
        YAVG=.5*(Y(M)+Y(MN))
C
C DETERMINE THE BOUNDARY CONDITIONS (LEFT CORNER)
C
        CALL BNDRYH (T,X(2),Y(M),U(1,2,M),AH,BH,CH,NPDE)
        CALL BNDRYV (T,X(1),Y(MN),U(1,1,MN),AV,BV,CV,NPDE)
        ITESTL=0
        ITESTH=0
        DO 40 K=1,NPDE
          IF (BH(K) .NE. 0.0) GO TO 30
          U(K,2,M) = CH(K)/AH(K)
          ITESTH=ITESTH+1
   30     IF (BV(K) .NE. 0.0) GO TO 40
          U(K,1,MN) = CV(K)/AV(K)
          ITESTL=ITESTL+1
   40     CONTINUE
        CALL BNDRYH (T,X(1),Y(M),U(1,1,M),AH,BH,CH,NPDE)
        CALL BNDRYV (T,X(1),Y(M),U(1,1,M),AV,BV,CV,NPDE)
        ITEST=0
        DO 70 K=1,NPDE
          IF (BV(K) .NE. 0.0 .AND. BH(K) .NE. 0.0) GO TO 70
          IF (BH(K) .EQ. 0.0 .AND. BV(K) .NE. 0.0) GO TO 50
          IF (BH(K) .NE. 0.0 .AND. BV(K) .EQ. 0.0) GO TO 60
          U(K,1,M)=(CH(K)/AH(K)+CV(K)/AV(K))*.5
          ITEST=ITEST+1
          GO TO 70
   50     U(K,1,M)=CH(K)/AH(K)
          GO TO 70
   60     U(K,1,M)=CV(K)/AV(K)
   70     CONTINUE
        IF (ITEST .EQ. NPDE) GO TO 130
        CALL BNDRYH (T,X(1),Y(M),U(1,1,M),AH,BH,CH,NPDE)
        CALL BNDRYV (T,X(1),Y(M),U(1,1,M),AV,BV,CV,NPDE)
C
C EVALUATE THE DIFFUSION COEFFICIENTS (LEFT CORNER)
C
        CALL DIFFH (T,X(1),Y(M),U(1,1,M),DH,NPDE)
        CALL DIFFV (T,X(1),Y(M),U(1,1,M),DV,NPDE)
C
C EVALUATE DU/DX AND DU/DY (LEFT CORNER)
C
        DO 120 K=1,NPDE
          IF (BV(K) .NE. 0.0) GO TO 90
          UX(K)=(U(K,2,M)-U(K,1,M))*DUXL(1)-(U(K,3,M)-U(K,1,M))*DUXL(2)
          GO TO 100
   90     UX(K)=(CV(K)-AV(K)*U(K,1,M))/BV(K)
  100     IF (BH(K) .NE. 0.0) GO TO 110
          UY(K)=((U(K,1,MN)-U(K,1,M))*DUY(IC,1)-(U(K,1,MNN)-U(K,1,M))
     *          *DUY(IC,2))*SIGN
          GO TO 120
  110     UY(K)=(CH(K)-AH(K)*U(K,1,M))/BH(K)
  120     CONTINUE
C
C CALCULATE THE U AVERAGES (LEFT CORNER)
C
  130   DO 140 K=1,NPDE
          UAVGH(K)=.5*(U(K,2,M)+U(K,1,M))
          UAVGV(K)=.5*(U(K,1,M)+U(K,1,MN))
          UDHR(K)=(U(K,2,M)-U(K,1,M))*DXIR(1)
          UDVA(K,1)=(U(K,1,MN)-U(K,1,M))*DYIA*SIGN
  140     CONTINUE
C
C CALCULATE THE DIFFUSION COEFFICIENTS AT THE MIDPOINTS OF THE
C INTERVALS BETWEEN THE LEFT CORNER AND THE NEIGHBORING POINTS
C
        CALL DIFFV (T,X(1),YAVG,UAVGV,DVS(1,1,1),NPDE)
        CALL DIFFH (T,XAVG(1),Y(M),UAVGH,DH(1,1,2),NPDE)
        IF (ITEST .EQ. NPDE) GO TO 160
C
C CALCULATE DUXX AND DUYY (LEFT CORNER)
C
        DO 150 L=1,NPDE
          DO 150 K=1,NPDE
            DH(K,L,1)=DXIC(1)*(DH(K,L,2)*UDHR(L)-DH(K,L,1)
     *                    *UX(L))
            DV(K,L)=DYIC*(DVS(K,L,1)*UDVA(L,1)-DV(K,L)
     *                    *UY(L))*SIGN
  150       CONTINUE
C
C EVALUATE THE RIGHT SIDE OF PDE (LEFT CORNER)
C
        CALL F (T,X(1),Y(M),U(1,1,M),UX,UY,DH,DV,DUDT(1,1,M),
     *          NPDE)
C
C SET DUDT = 0 FOR KNOWN BOUNDARY CONDITIONS
C
  160   DO 170 K=1,NPDE
          IF (BH(K) .EQ. 0.0 .OR. BV(K) .EQ. 0.0) DUDT(K,1,M)=0.0
  170     CONTINUE
C
C EVALUATE THE RIGHT CORNER BOUNDARY CONDITIONS
C
        CALL BNDRYH (T,X(NX),Y(M),U(1,NX,M),AH,BH,CH,NPDE)
        CALL BNDRYV (T,X(NX),Y(M),U(1,NX,M),AV,BV,CV,NPDE)
        ITEST = 0
        DO 200 K=1,NPDE
          IF (BV(K) .NE. 0.0 .AND. BH(K) .NE. 0.0) GO TO 200
          IF (BH(K) .EQ. 0.0 .AND. BV(K) .NE. 0.0) GO TO 180
          IF (BH(K) .NE. 0.0 .AND. BV(K) .EQ. 0.0) GO TO 190
          U(K,NX,M)=(CV(K)/AV(K)+CH(K)/AH(K))*.5
              ITEST=ITEST+1
          GO TO 200
  180     U(K,NX,M)=CH(K)/AH(K)
          GO TO 200
  190     U(K,NX,M)=CV(K)/AV(K)
  200       CONTINUE
        IBCK=1
        IFWD=2
C
C LOOP ON THE HORIZONTAL BOUNDARY MESH POINTS FROM X(2) TO X(NX-1)
C
        DO 300 I=2,ILIM
          K=IBCK
          IBCK=IFWD
          IFWD=K
C
C EVALUATE THE HORIZONTAL BOUNDARY CONDITIONS
C
          ITESTC=ITESTH
          IF (I .EQ. ILIM) GO TO 220
          CALL BNDRYH (T,X(I+1),Y(M),U(1,I+1,M),AH,BH,CH,
     *                 NPDE)
          ITESTH=0
          DO 210 K=1,NPDE
            IF (BH(K) .NE. 0.0) GO TO 210
            U(K,I+1,M) = CH(K)/AH(K)
            ITESTH=ITESTH+1
  210       CONTINUE
  220     CONTINUE
          CALL BNDRYH (T,X(I),Y(M),U(1,I,M),AH,BH,CH,NPDE)
          IF (ITESTC .EQ. NPDE) GO TO 260
C
C CALCULATE DU/DX AND DU/DY (HORIZONTAL BOUNDARY)
C
          DO 250 K=1,NPDE
            UX(K)=(U(K,I+1,M)-U(K,I-1,M))*DXI(I)
            IF (BH(K) .NE. 0.0) GO TO 240
            UY(K)=((U(K,I,MN)-U(K,I,M))*DUY(IC,1)-(U(K,I,MNN)-U(K,I,M))
     *            *DUY(IC,2))*SIGN
            GO TO 250
  240       UY(K)=(CH(K)-AH(K)*U(K,I,M))/BH(K)
  250       CONTINUE
C
C DETERMINE U AVERAGE (HORIZONTAL BOUNDARY)
C
  260     DO 270 K=1,NPDE
            UAVGV(K)=(U(K,I,M)+U(K,I,MN))*.5
            UAVGH(K)=(U(K,I+1,M)+U(K,I,M))*.5
            UDHL(K)=UDHR(K)
            UDHR(K)=(U(K,I+1,M)-U(K,I,M))*DXIR(I)
            UDVA(K,I)=(U(K,I,MN)-U(K,I,M))*DYIA*SIGN
  270       CONTINUE
C
C EVALUATE THE DIFFUSION COEFFICIENTS (HORIZONTAL BOUNDARY)
C
          CALL DIFFV (T,X(I),YAVG,UAVGV,DVS(1,1,I),NPDE)
          CALL DIFFH (T,XAVG(I),Y(M),UAVGH,DH(1,1,IFWD),NPDE)
          IF (ITESTC .EQ. NPDE) GO TO 290
          CALL DIFFV (T,X(I),Y(M),U(1,I,M),DV,NPDE)
C
C EVALUATE DUXX AND DUYY (HORIZONTAL BOUNDARY)
C
          DO 280 L=1,NPDE
            DO 280 K=1,NPDE
              DH(K,L,IBCK)=DXIC(I)*(DH(K,L,IFWD)*UDHR(L)-
     *                         DH(K,L,IBCK)*UDHL(L))
              DV(K,L)=DYIC*(DVS(K,L,I)*UDVA(L,I)-DV(K,L)
     *                         *UY(L))*SIGN
  280         CONTINUE
C
C EVALUATE THE RIGHT SIDE OF THE PDE (HORIZONTAL BOUNDARY)
C
          CALL F (T,X(I),Y(M),U(1,I,M),UX,UY,DH(1,1,IBCK),
     *            DV,DUDT(1,I,M),NPDE)
C
C SET DUDT = 0 FOR KNOWN BOUNDARY CONDITIONS
C
  290   DO 300 K=1,NPDE
          IF (BH(K) .EQ. 0.0) DUDT(K,I,M)=0.0
  300     CONTINUE
C
C COMPLETE EVALUATING THE RIGHT CORNER
C
        CALL BNDRYV (T,X(NX),Y(MN),U(1,NX,MN),AV,BV,CV,NPDE)
        ITESTR=0
        DO 305 K=1,NPDE
          IF (BV(K) .NE. 0.0) GO TO 305
          U(K,NX,MN) = CV(K)/AV(K)
          ITESTR=ITESTR+1
  305     CONTINUE
        CALL BNDRYH (T,X(NX),Y(M),U(1,NX,M),AH,BH,CH,NPDE)
        CALL BNDRYV (T,X(NX),Y(M),U(1,NX,M),AV,BV,CV,NPDE)
        IF (ITEST .EQ. NPDE) GO TO 360
C
C EVALUATE THE DIFFUSION COEFFICIENTS (RIGHT CORNER)
C
        CALL DIFFH (T,X(NX),Y(M),U(1,NX,M),DH(1,1,IBCK),NPDE)
        CALL DIFFV (T,X(NX),Y(M),U(1,NX,M),DV,NPDE)
C
C EVALUATE DU/DX AND DU/DY (RIGHT CORNER)
C
        DO 350 K=1,NPDE
          IF (BV(K) .NE. 0.0) GO TO 320
          UX(K)=(U(K,NX,M)-U(K,ILIM,M))*DUXR(1)-(U(K,NX,M)
     *          -U(K,ILIM-1,M))*DUXR(2)
          GO TO 330
  320     UX(K)=(CV(K)-AV(K)*U(K,NX,M))/BV(K)
  330     IF (BH(K) .NE. 0.0) GO TO 340
          UY(K)=((U(K,NX,MN)-U(K,NX,M))*DUY(IC,1)-(U(K,NX,MNN)
     *          -U(K,NX,M))*DUY(IC,2))*SIGN
          GO TO 350
  340     UY(K)=(CH(K)-AH(K)*U(K,NX,M))/BH(K)
  350     CONTINUE
C
C EVALUATE THE VERTICAL U AVERAGE (RIGHT CORNER)
C
  360   DO 370 K=1,NPDE
          UAVGV(K)=(U(K,NX,M)+U(K,NX,MN))*.5
          UDVA(K,NX)=(U(K,NX,MN)-U(K,NX,M))*DYI*SIGN
  370     CONTINUE
C
C EVALUATE THE VERTICAL DIFFUSION COEFFICIENT (ABOVE THE RIGHT CORNER)
C
        CALL DIFFV (T,X(NX),YAVG,UAVGV,DVS(1,1,NX),NPDE)
        IF (ITEST .EQ. NPDE) GO TO 390
C
C EVALUATE DUXX AND DUYY (RIGHT CORNER)
C
        DO 380 L=1,NPDE
          DO 380 K=1,NPDE
            DH(K,L,IBCK)=DXIC(NX)*(DH(K,L,IBCK)*UX(L)-
     *                       DH(K,L,IFWD)*UDHR(L))
            DV(K,L)=DYIC*(DVS(K,L,NX)*UDVA(L,NX)-DV(K,L)
     *                       *UY(L))*SIGN
  380       CONTINUE
C
C EVALUATE THE RIGHT SIDE OF THE PDE (RIGHT CORNER)
C
        CALL F (T,X(NX),Y(M),U(1,NX,M),UX,UY,DH(1,1,IBCK),DV,
     *          DUDT(1,NX,M),NPDE)
  390   DO 400 K=1,NPDE
          IF (BH(K) .EQ. 0.0 .OR. BV(K) .EQ. 0.0) DUDT(K,NX,M)=0.0
  400     CONTINUE
        IF (IC .NE. 1) GO TO 410
        DO 405 IIC=1,NX
          DO 405 L=1,NPDE
            DO 405 K=1,NPDE
  405         DVST(K,L,IIC)=DVS(K,L,IIC)
  410   CONTINUE
C
C DETERMINE THE U VALUES FOR THE J-TH ROW
C
      IC=2
      DO 780 J=2,ILIMY
        IF (J .NE. ILIMY) GO TO 500
        IC=1
  500   DYI=1./(Y(J+1)-Y(J-1))
        DYIC=2.*DYI
        DYIB=DYIA
        DYIA=1./(Y(J+1)-Y(J))
        YAVG=(Y(J+1)+Y(J))*.5
C
C DETERMINE THE LEFT BOUNDARY (J-TH ROW)
C
        ITESTC=ITESTL
        IF (IC .EQ. 1) GO TO 510
        CALL BNDRYV (T,X(1),Y(J+1),U(1,1,J+1),AV,BV,CV,NPDE)
        ITESTL=0
        DO 505 K=1,NPDE
          IF (BV(K) .NE. 0.0) GO TO 505
          U(K,1,J+1) = CV(K)/AV(K)
          ITESTL=ITESTL+1
  505     CONTINUE
  510   CONTINUE
        CALL BNDRYV (T,X(1),Y(J),U(1,1,J),AV,BV,CV,NPDE)
        IF (ITESTC .EQ. NPDE) GO TO 560
C
C EVALUATE DU/DX AND DU/DY (LEFT BOUNDARY,J-TH ROW)
C
        DO 550 K=1,NPDE
          IF (BV(K) .NE. 0.0) GO TO 530
          UX(K)=(U(K,2,J)-U(K,1,J))*DUXL(1)-(U(K,3,J)-U(K,1,J))*DUXL(2)
          GO TO 540
  530     UX(K)=(CV(K)-AV(K)*U(K,1,J))/BV(K)
  540     UY(K)=(U(K,1,J+1)-U(K,1,J-1))*DYI
  550     CONTINUE
C
C EVALUATE THE HORIZONTAL U AVERAGE (J-TH ROW)
C
  560   DO 570 K=1,NPDE
          UAVGH(K)=(U(K,2,J)+U(K,1,J))*.5
          UDHR(K)=(U(K,2,J)-U(K,1,J))*DXIR(1)
          UDVB(K)=UDVA(K,1)
          UDVA(K,1)=(U(K,1,J+1)-U(K,1,J))*DYIA
  570     CONTINUE
C
C EVALUATE THE DIFFUSION COEFFICIENTS (LEFT BOUNDARY,J-TH ROW)
C
        DO 580 L=1,NPDE
          DO 580 K=1,NPDE
            DV(K,L)=DVS(K,L,1)
  580       CONTINUE
        IF (IC .EQ. 2) GO TO 582
        DO 581 L=1,NPDE
          DO 581 K=1,NPDE
  581       DVS(K,L,1)=DVST(K,L,1)
  582   CONTINUE
        CALL DIFFH (T,X(1),Y(J),U(1,1,J),DH,NPDE)
        CALL DIFFH (T,XAVG(1),Y(J),UAVGH,DH(1,1,2),NPDE)
        IF (IC .EQ. 1) GO TO 590
C
C EVALUATE THE VERTICAL U AVERAGE (J-TH ROW)
C
        DO 585 K=1,NPDE
          UAVGV(K)=(U(K,1,J+1)+U(K,1,J))*.5
  585     CONTINUE
        CALL DIFFV (T,X(1),YAVG,UAVGV,DVS(1,1,1),NPDE)
C
C EVALUATE DUXX AND DUYY (LEFT BOUNDARY,J-TH ROW)
C
  590   IF (ITESTC .EQ. NPDE) GO TO 610
        DO 600 L=1,NPDE
          DO 600 K=1,NPDE
            DH(K,L,1)=DXIC(1)*(DH(K,L,2)*UDHR(L)-
     *                    DH(K,L,1)*UX(L))
            DV(K,L)=DYIC*(DVS(K,L,1)*UDVA(L,1)-DV(K,L)
     *                    *UDVB(L))
  600       CONTINUE
C
C EVALUATE THE RIGHT SIDE OF THE PDE (LEFT BOUNDARY,J-TH ROW)
C
        CALL F (T,X(1),Y(J),U(1,1,J),UX,UY,DH,DV,DUDT(1,1,J),NPDE)
C
C SET DUDT = 0 FOR KNOWN LEFT BOUNDARY CONDITIONS
C
  610   DO 620 K=1,NPDE
          IF (BV(K) .EQ. 0.0) DUDT(K,1,J)=0.0
  620     CONTINUE
C
C LOOP TO EVALUATE U IN THE CENTER OF THE GRID
C (2.LT.I.LT.NX-1, 2.LT.J.LT.NY-1)
C
        IBCK=1
        IFWD=2
        DO 680 I=2,ILIM
          K=IBCK
          IBCK=IFWD
          IFWD=K
C
C CALCULATE THE HORIZONTAL U AVERAGE, DU/DX AND DU/DY
C (I-TH POINT OF THE J-TH ROW)
C
          DO 640 K=1,NPDE
            UAVGH(K)=(U(K,I+1,J)+U(K,I,J))*.5
            UX(K)=(U(K,I+1,J)-U(K,I-1,J))*DXI(I)
            UY(K)=(U(K,I,J+1)-U(K,I,J-1))*DYI
  640       CONTINUE
C
C EVALUATE THE DIFFUSION COEFFICIENTS (I-TH POINT OF THE J-TH ROW)
C
          DO 650 L=1,NPDE
            DO 650 K=1,NPDE
              DV(K,L)=DVS(K,L,I)
  650         CONTINUE
          IF (IC .EQ. 2) GO TO 652
          DO 651 L=1,NPDE
            DO 651 K=1,NPDE
  651         DVS(K,L,I)=DVST(K,L,I)
  652     CONTINUE
          CALL DIFFH (T,XAVG(I),Y(J),UAVGH,DH(1,1,IFWD),NPDE)
          IF (IC .EQ. 1) GO TO 660
          DO 655 K=1,NPDE
            UAVGV(K)=(U(K,I,J+1)+U(K,I,J))*.5
  655       CONTINUE
          CALL DIFFV (T,X(I),YAVG,UAVGV,DVS(1,1,I),NPDE)
C
C EVALUATE DUXX AND DUYY (I-TH POINT OF THE J-TH ROW)
C
  660     DO 670 L=1,NPDE
            UDHL(L)=UDHR(L)
            UDHR(L)=(U(L,I+1,J)-U(L,I,J))*DXIR(I)
            UDVB(L)=UDVA(L,I)
            UDVA(L,I)=(U(L,I,J+1)-U(L,I,J))*DYIA
            DO 670 K=1,NPDE
              DH(K,L,IBCK)=DXIC(I)*(DH(K,L,IFWD)*UDHR(L)-
     *                         DH(K,L,IBCK)*UDHL(L))
              DV(K,L)=DYIC*(DVS(K,L,I)*UDVA(L,I)-DV(K,L)
     *                         *UDVB(L))
  670         CONTINUE
C
C EVALUATE THE RIGHT SIDE OF THE PDE (I-TH POINT OF THE J-TH ROW)
C
          CALL F (T,X(I),Y(J),U(1,I,J),UX,UY,DH(1,1,IBCK),DV,
     *            DUDT(1,I,J),NPDE)
  680     CONTINUE
C
C EVALUATE THE RIGHT BOUNDARY FOR THE J-TH ROW
C
        ITESTC=ITESTR
        IF (IC .EQ. 1) GO TO 690
        CALL BNDRYV (T,X(NX),Y(J+1),U(1,NX,J+1),AV,BV,CV,NPDE)
        ITESTR=0
        DO 685 K=1,NPDE
          IF (BV(K) .NE. 0.0) GO TO 685
          U(K,NX,J+1) = CV(K)/AV(K)
          ITESTR=ITESTR+1
  685     CONTINUE
  690   CONTINUE
        CALL BNDRYV (T,X(NX),Y(J),U(1,NX,J),AV,BV,CV,NPDE)
        IF (ITESTC .EQ. NPDE) GO TO 730
C
C EVALUATE DU/DX AND DU/DY (RIGHT BOUNDARY,J-TH ROW)
C
        DO 710 K=1,NPDE
          UY(K)=(U(K,NX,J+1)-U(K,NX,J-1))*DYI
          IF (BV(K) .NE. 0.0) GO TO 700
          UX(K)=(U(K,NX,J)-U(K,ILIM,J))*DUXR(1)-(U(K,NX,J)
     *          -U(K,ILIM-1,J))*DUXR(2)
          GO TO 710
  700     UX(K)=(CV(K)-AV(K)*U(K,NX,J))/BV(K)
  710     CONTINUE
C
C EVALUATE THE DIFFUSION COEFFICIENTS (RIGHT BOUNDARY,J-TH ROW)
C
        DO 720 L=1,NPDE
          DO 720 K=1,NPDE
            DV(K,L)=DVS(K,L,NX)
  720       CONTINUE
        IF (IC .EQ. 2) GO TO 725
        DO 721 L=1,NPDE
          DO 721 K=1,NPDE
  721       DVS(K,L,NX)=DVST(K,L,NX)
  725   CONTINUE
        CALL DIFFH (T,X(NX),Y(J),U(1,NX,J),DH(1,1,IBCK),NPDE)
  730   DO 740 K=1,NPDE
          UAVGV(K)=(U(K,NX,J+1)+U(K,NX,J))*.5
          UDVB(K)=UDVA(K,NX)
          UDVA(K,NX)=(U(K,NX,J+1)-U(K,NX,J))*DYIA
  740     CONTINUE
        IF (IC .EQ. 1) GO TO 750
        CALL DIFFV (T,X(NX),YAVG,UAVGV,DVS(1,1,NX),NPDE)
  750   IF (ITESTC .EQ. NPDE) GO TO 770
C
C EVALUATE DUXX AND DUYY (RIGHT BOUNDARY,J-TH ROW)
C
        DO 760 L=1,NPDE
          DO 760 K=1,NPDE
            DH(K,L,IBCK)=DXIC(NX)*(DH(K,L,IBCK)*UX(L)-
     *                       DH(K,L,IFWD)*UDHR(L))
            DV(K,L)=DYIC*(DVS(K,L,NX)*UDVA(L,NX)-DV(K,L)
     *                       *UDVB(L))
  760       CONTINUE
C
C EVALUATE THE RIGHT SIDE OF THE PDE (RIGHT BOUNDARY,J-TH ROW)
C
        CALL F (T,X(NX),Y(J),U(1,NX,J),UX,UY,DH(1,1,IBCK),DV,
     *          DUDT(1,NX,J),NPDE)
C
C SET DUDT = 0 FOR KNOWN BOUNDARY VALUES
C
  770 DO 780 K=1,NPDE
        IF (BV(K) .EQ. 0.0) DUDT(K,NX,J)=0.0
  780   CONTINUE
      RETURN
      END
      SUBROUTINE STRSET (N,MF,INDEX,IWORK,LOUT,ML,MU)                    STRSET
C
C***********************************************************************
C
C       STRSET IS DESIGNED TO EXTRACT AND CHECK THE INPUT INTEGER
C  PARAMETERS FROM THE FIRST SIX LOCATIONS OF THE ARRAY IWORK AND TO
C  PROVIDE DYNAMIC DIMENSIONING FOR THE ARRAYS IN THE ROUTINE PDETWO
C  AND THE MODIFIED GEARB PACKAGE.  STRSET ESTABLISHES THE DIMENSION
C  POINTERS, IW1 - IW27 FOR WORK, THE USER DIMENSIONED WORKING STOR-
C  AGE ARRAY.  THE CORRESPONDENCE BETWEEN THIS WORKING STORAGE AND THE
C  ARRAYS IN PDETWO AND THE MODIFIED GEARB IS LISTED BELOW.  STRSET
C  ALSO SETS THE LOWER AND THE UPPER BANDWIDTHS ML AND MU FOR THE
C  JACOBIAN MATRIX.
C
C       PARAMETERS ARE SELF-EXPLANATORY IN THE COMMENT STATEMENTS
C  BELOW.
C
C***********************************************************************
C  CORRESPONDENCE BETWEEN WORK AND THE ARRAYS IN PDETWO AND THE MODI-
C  FIED GEARB PACKAGE
C***********************************************************************
C
C                         CORRESPONDING ARRAY
C   LOCATION IN WORK           IN PDETWO           LENGTH
C  ------------------    ---------------------    ----------------------
C
C     WORK(1)                    DVS                NPDE * NPDE * NX
C     WORK(IW1)                  DVST               NPDE * NPDE * NX
C     WORK(IW2)                  DV                 NPDE * NPDE
C     WORK(IW3)                  DH                 NPDE * NPDE * 2
C     WORK(IW4)                  AH                 NPDE
C     WORK(IW5)                  BH                 NPDE
C     WORK(IW6)                  CH                 NPDE
C     WORK(IW7)                  AV                 NPDE
C     WORK(IW8)                  BV                 NPDE
C     WORK(IW9)                  CV                 NPDE
C     WORK(IW10)                 UX                 NPDE
C     WORK(IW11)                 UY                 NPDE
C     WORK(IW12)                 DXI                NX
C     WORK(IW13)                 DXIR               NX
C     WORK(IW14)                 DXIC               NX
C     WORK(IW15)                 XAVG               NX
C     WORK(IW16)                 UDVA               NPDE * NX
C     WORK(IW17)                 UDHR               NPDE
C     WORK(IW18)                 UDVB               NPDE
C     WORK(IW19)                 UDHL               NPDE
C     WORK(IW20)                 UAVGH              NPDE
C     WORK(IW21)                 UAVGV              NPDE
C
C                         CORRESPONDING ARRAY
C   LOCATION IN WORK     IN THE MODIFIED GEARB     LENGTH
C                               PACKAGE
C  ------------------    ---------------------    ----------------------
C
C     WORK(IW22)                 U                  NODE*(MORDER+1)
C     WORK(IW23)                 UMAX               NODE
C     WORK(IW24)                 ERROR              NODE
C     WORK(IW25)                 SAVE1              NODE
C     WORK(IW26)                 SAVE2              NODE
C     WORK(IW27)                 PW                 LPW
C
C  WHERE NODE = NPDE*NX*NY, NPDE IS THE NUMBER OF PARTIAL
C  DIFFERENTIAL EQUATIONS, NX IS THE NUMBER OF MESH POINTS IN
C  THE HORIZONTAL DIRECTION AND NY IS THE NUMBER OF MESH POINTS
C  IN THE VERTICAL DIRECTION.  SEE COMMENTS IN PDETWO FOR THE
C  DEFINITION OF MORDER AND LPW.
C
C***********************************************************************
C
      COMMON /IPARAM/ NPDE,NX,NY,MORDER,LWORK,LIWORK
      COMMON /ISTORE/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11,
     *                IW12,IW13,IW14,IW15,IW16,IW17,IW18,IW19,IW20,IW21,
     *                IW22,IW23,IW24,IW25,IW26,IW27
      DIMENSION IWORK(N)
C
C EXTRACT THE INPUT INTEGER PARAMETERS FROM THE ARRAY IWORK
C
      NPDE = IWORK(1)
      NX = IWORK(2)
      NY = IWORK(3)
      MORDER = IWORK(4)
      NRWK = IWORK(5)
      NRIWK = IWORK(6)
C
C CHECK IF NUMBER OF ODE, NX, AND NY ARE CORRECT
C
      NODE = NPDE*NX*NY
      IF (N .NE. NODE) GO TO 110
      IF (NX .LT. 3) GO TO 120
      IF (NY .LT. 3) GO TO 130
C
C CHECK IF METH AND MITER OBTAINED FROM THE METHOD FLAG MF ARE WITHIN
C THE LIMIT OF THE ALLOWED VALUES
C
      METH = MF/10
      IF (METH .LT. 1 .OR. METH .GT. 2) GO TO 140
      MITER = MF - 10*METH
      IF (MITER .LT. 0 .OR. MITER .GT. 3) GO TO 150
C
C CHECK IF THE MAXIMUM ORDER OF THE METHOD IS WITHIN THE LIMIT
C IF IT EXCEEDS THE UPPER LIMIT, THE MAXIMUM ORDER AVALIABLE IS USED
C
      IF (MORDER .LT. 1) GO TO 160
      IF (METH .EQ. 1 .AND. MORDER .GT. 12) MORDER = 12
      IF (METH .EQ. 2 .AND. MORDER .GT. 5) MORDER = 5
C
C CHECK IF THE LENGTHS OF THE ARRAYS WORK AND IWORK ARE LARGE ENOUGH
C
      LU = NODE*(MORDER+1)
      LPW = 1
      IF(MITER .EQ. 1 .OR. MITER .EQ. 2) LPW = (3*(NX+1)*NPDE - 2)*NODE
      IF(MITER .EQ. 3) LPW = NODE
      LWORK = NPDE*(NPDE*(2*NX+3) + 13 + NX) + 4*(NX + NODE) + LU + LPW
      IF(NRWK .LT. LWORK) GO TO 170
      LIWORK = NODE
      IF(NRIWK .LT. LIWORK) GO TO 180
      GO TO 200
C
C PRINT ERROR MESSAGE
C
 110  WRITE (LOUT,115) N,NODE
 115  FORMAT(//16H ILLEGAL INPUT..,7H NODE =,I9,2X,28HIS NOT EQUAL TO NP
     1DE.NX.NY =,I9//)
      INDEX = -4
      GO TO 300
C
 120  WRITE (LOUT,125) NX
 125  FORMAT(//16H ILLEGAL INPUT..,5H NX =,I5,2X,14HIS LESS THAN 3//)
      INDEX = -4
      GO TO 300
C
 130  WRITE (LOUT,135) NY
 135  FORMAT(//16H ILLEGAL INPUT..,5H NY =,I5,2X,14HIS LESS THAN 3//)
      INDEX = -4
      GO TO 300
C
 140  WRITE (LOUT,145) METH
 145  FORMAT(//16H ILLEGAL INPUT..,7H METH =,I3,2X,14HIS NOT ALLOWED//)
      INDEX = -4
      GO TO 300
C
 150  WRITE (LOUT,155) MITER
 155  FORMAT(//16H ILLEGAL INPUT..,8H MITER =,I3,2X,14HIS NOT ALLOWED//)
      INDEX = -4
      GO TO 300
C
 160  WRITE (LOUT,165) MORDER
 165  FORMAT(//16H ILLEGAL INPUT..,41H MAXIMUM ORDER OF THE METHODS AVAI
     1LABLE =,I3,2X,14HIS NOT ALLOWED//)
      INDEX = -4
      GO TO 300
C
 170  WRITE(LOUT,175) NRWK,LWORK
 175  FORMAT(//23H INSUFFICIENT STORAGE..,11H THE LENGTH,I9,2X,60HOF THE
     1 WORK ARRAY IS LESS THAN THE REQUIRED MINIMUM STORAGE ,I9//)
      INDEX = -6
      GO TO 300
C
 180  WRITE(LOUT,185) NRIWK,LIWORK
 185  FORMAT(//23H INSUFFICIENT STORAGE..,11H THE LENGTH,I9,2X,60HOF THE
     1 IWORK ARRAY IS LESS THAN THE REQUIRED MINIMUM STORAGE,I9//)
      INDEX = -6
      GO TO 300
C
C ASSIGN DIMENSIONS FOR THE ARRAYS
C
 200  CONTINUE
      NPDE2 = NPDE*NPDE
      IW1 = NPDE2*NX + 1
      IW2 = IW1 + NPDE2*NX
      IW3 = IW2 + NPDE2
      IW4 = IW3 + 2*NPDE2
      IW5 = IW4 + NPDE
      IW6 = IW5 + NPDE
      IW7 = IW6 + NPDE
      IW8 = IW7 + NPDE
      IW9 = IW8 + NPDE
      IW10 = IW9 + NPDE
      IW11 = IW10 + NPDE
      IW12 = IW11 + NPDE
      IW13 = IW12 + NX
      IW14 = IW13 + NX
      IW15 = IW14 + NX
      IW16 = IW15 + NX
      IW17 = IW16 + NPDE * NX
      IW18 = IW17 + NPDE
      IW19 = IW18 + NPDE
      IW20 = IW19 + NPDE
      IW21 = IW20 + NPDE
      IW22 = IW21 + NPDE
      IW23 = IW22 + NODE*(MORDER+1)
      IW24 = IW23 + NODE
      IW25 = IW24 + NODE
      IW26 = IW25 + NODE
      IW27 = IW26 + NODE
C
C SET BANDWIDTH MU AND ML FOR THE JACOBIAN
C
      MU = NPDE*(NX+1) - 1
      ML = MU
 300  RETURN
      END
      SUBROUTINE DRIVEP (N,T0,H0,U0,TOUT,EPS,MF,INDEX,WORK,IWORK,XM,YM)  DRIVEP
C        THIS IS A MODIFIED VERSION OF GEARB(1),
C        A PACKAGE FOR THE SOLUTION OF THE INITIAL VALUE
C        PROBLEM FOR SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS,
C                  DU/DT = F(U,T),    U = (U(1),U(2),...,U(N)).
C        GEARB IS A VARIANT OF THE GEAR PACKAGE TO BE USED WHEN
C        THE JACOBIAN MATRIX DF/DU HAS BANDED OR NEARLY BANDED FORM.
C        SUBROUTINE DRIVEP IS A DRIVER ROUTINE FOR THE MODIFIED
C        GEARB PACKAGE.  THE PACKAGE IS MODIFIED EXCLUSIVELY FOR USE
C        WITH PDETWO AND PSETM.
C
C        THIS PACKAGE WAS CONSTRUCTED SO AS TO CONFORM TO AS MANY
C        ANSI-FORTRAN RULES AS WAS CONVENIENTLY POSSIBLE.  THE FORTRAN
C        USED VIOLATES ONLY ONE ANSI STANDARD - AN ARRAY NAME APPEARS
C        IN A DATA STATEMENT IN THE SUBROUTINE COSET.
C
C
C                       REFERENCES
C       1.  A. C. HINDMARSH, GEARB.. SOLUTION OF ORDINARY DIFFERENTIAL
C             EQUATIONS HAVING BANDED JACOBIAN, UCID-30059 REV. 2,
C             LAWRENCE LIVERMORE LABORATORY, P.O.BOX 808, LIVERMORE,
C             CA 94550, JUNE 1977.
C
C-----------------------------------------------------------------------
C  DRIVEP IS TO BE CALLED ONCE FOR EACH OUTPUT VALUE OF T, AND
C  IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR, STIFFP.
C
C  THE INPUT PARAMETERS ARE..
C    N     =  THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.
C               IT SHOULD BE SET TO NODE = NPDE*NX*NY.  SEE IWORK
C               BELOW FOR THE DEFINITION OF NPDE, NX, NY.
C    T0    =  THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE
C               (USED ONLY ON FIRST CALL).
C    H0    =  THE NEXT STEP SIZE IN T (USED FOR INPUT ONLY ON THE
C               FIRST CALL).
C    U0    =  A VECTOR OF LENGTH N CONTAINING THE INITIAL VALUES OF
C               U (USED FOR INPUT ONLY ON FIRST CALL).
C    TOUT  =  THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT.
C               INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT
C               AND THE PACKAGE WILL INTERPOLATE TO T = TOUT.
C    EPS   =  THE RELATIVE ERROR BOUND  (USED ONLY ON THE
C               FIRST CALL, UNLESS INDEX = -1).  SINGLE STEP ERROR
C               ESTIMATES DIVIDED BY UMAX(I) WILL BE KEPT LESS THAN
C               EPS IN ROOT-MEAN-SQUARE NORM (I.E. EUCLIDEAN NORM
C               DIVIDED BY SQRT(N) ).  THE VECTOR UMAX OF WEIGHTS
C               IS COMPUTED IN DRIVEP.  INITIALLY UMAX(I) IS
C               ABS(U(I)), WITH A DEFAULT VALUE OF 1 IF Y(I) = 0
C               INITIALLY.  THEREAFTER, UMAX(I) IS THE LARGEST VALUE
C               OF ABS(U(I)) SEEN SO FAR, OR THE INITIAL UMAX(I) IF
C               THAT IS LARGER.  TO ALTER EITHER OF THESE, CHANGE THE
C               APPROPRIATE STATEMENTS IN THE DO-LOOPS ENDING AT
C               STATEMENTS 10 AND 70 BELOW.
C    MF    =  THE METHOD FLAG  (USED ONLY ON FIRST CALL, UNLESS
C               INDEX = -1).  ALLOWED VALUES ARE 10, 11, 12, 13,
C               20, 21, 22, 23.  MF HAS TWO DECIMAL DIGITS, METH
C               AND MITER  (MF = 10*METH + MITER).
C               METH IS THE BASIC METHOD INDICATOR..
C                 METH = 1  MEANS THE ADAMS METHODS.
C                 METH = 2  MEANS THE BACKWARD DIFFERENTIATION
C                           FORMULAS (BDF), OR STIFF METHODS OF GEAR.
C               MITER IS THE ITERATION METHOD INDICATOR..
C                 MITER = 0 MEANS FUNCTIONAL ITERATION (NO PARTIAL
C                           DERIVATIVES NEEDED).
C                 MITER = 1 MEANS CHORD METHOD WITH ANALYTIC JACOBIAN.
C                           FOR THIS USER SUPPLIES SUBROUTINE
C                           PDB  (SEE DESCRIPTION BELOW).
C                 MITER = 2 MEANS CHORD METHOD WITH JACOBIAN CALCULATED
C                           INTERNALLY BY FINITE DIFFERENCES.
C                 MITER = 3 MEANS CHORD METHOD WITH JACOBIAN REPLACED
C                           BY A DIAGONAL APPROXIMATION BASED ON A
C                           DIRECTIONAL DERIVATIVE.
C    INDEX =  INTEGER USED ON INPUT TO INDICATE TYPE OF CALL,
C               WITH THE FOLLOWING VALUES AND MEANINGS..
C                  1    THIS IS THE FIRST CALL FOR THIS PROBLEM.
C                  0    THIS IS NOT THE FIRST CALL FOR THIS PROBLEM,
C                       AND INTEGRATION IS TO CONTINUE.
C                 -1    THIS IS NOT THE FIRST CALL FOR THE PROBLEM,
C                       AND THE USER HAS RESET N, EPS, AND/OR MF.
C                  2    SAME AS 0 EXCEPT THAT TOUT IS TO BE HIT
C                       EXACTLY (NO INTERPOLATION IS DONE).
C                       ASSUMES TOUT .GE. THE CURRENT T.
C                  3    SAME AS 0 EXCEPT CONTROL RETURNS TO CALLING
C                       PROGRAM AFTER ONE STEP.  TOUT IS IGNORED.
C               SINCE THE NORMAL OUTPUT VALUE OF INDEX IS 0,
C               IT NEED NOT BE RESET FOR NORMAL CONTINUATION.
C    WORK  =  FLOATING POINT WORKING ARRAY.  SEE STRSET FOR ADDITIONAL
C               DETAILS.
C    IWORK =  INTEGER WORKING ARRAY.  THE FIRST SIX LOCATIONS MUST BE
C               DEFINED AS FOLLOWS..
C               IWORK(1) = NPDE
C               IWORK(2) = NX
C               IWORK(3) = NY
C               IWORK(4) = MORDER
C               IWORK(5) = NRWK
C               IWORK(6) = NRIWK
C               WHERE NPDE IS THE NUMBER OF PARTIAL DIFFERENTIAL EQUA-
C               TIONS, NX IS THE NUMBER OF MESH POINTS IN THE HORIZON-
C               TAL DIRECTION AND NY IS THE NUMBER OF MESH POINTS IN
C               THE VERTICAL DIRECTION, MORDER IS THE MAXIMUM ORDER OF
C               METHOD TO BE USED IN THE MODIFIED GEARB, NRWK IS THE
C               LENGTH OF THE ARRAY WORK AND NRIWK IS THE LENGTH OF THE
C               ARRAY IWORK.  MORDER MUST BE LESS THAN OR EQUAL TO 12
C               IF METH = 1 OR IF METH = 2, IT MUST BE LESS THAN OR
C               EQUAL TO 5.
C    XM    =  ARE THE MESH POINTS IN THE HORIZONTAL DIRECTION.
C    YM    =  ARE THE MESH POINTS IN THE VERTICAL DIRECTION.
C
C  AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURRED AND A NORMAL
C  CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN.
C  ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL.
C  A CHANGE OF PARAMETERS WITH INDEX = -1 CAN BE MADE AFTER
C  EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN.
C
C  THE OUTPUT PARAMETERS ARE..
C    H0    =  THE STEP SIZE H USED LAST, WHETHER SUCCESSFULLY OR NOT.
C    U0    =  THE COMPUTED VALUES OF U AT T = TOUT.
C    TOUT  =  THE OUTPUT VALUE OF T.  IF INTEGRATION WAS SUCCESSFUL,
C               AND THE INPUT VALUE OF INDEX WAS NOT 3, TOUT IS
C               UNCHANGED FROM ITS INPUT VALUE.  OTHERWISE, TOUT
C               IS THE CURRENT VALUE OF T TO WHICH INTEGRATION
C               HAS BEEN COMPLETED.
C    INDEX =  INTEGER USED ON OUTPUT TO INDICATE RESULTS,
C               WITH THE FOLLOWING VALUES AND MEANINGS..
C          0    INTEGRATION WAS COMPLETED TO TOUT OR BEYOND.
C         -1    THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE
C               ERROR TEST EVEN AFTER REDUCING H BY A FACTOR OF
C               1.E10 FROM ITS INITIAL VALUE.
C         -2    AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS
C               HALTED EITHER BY REPEATED ERROR TEST FAILURES OR BY
C               A TEST ON EPS.  TOO MUCH ACCURACY HAS BEEN REQUESTED.
C         -3    THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE
C               CORRECTOR CONVERGENCE EVEN AFTER REDUCING H BY A
C               FACTOR OF 1.E10 FROM ITS INITIAL VALUE.
C         -4    IMMEDIATE HALT BECAUSE OF ILLEGAL VALUES OF INPUT
C               PARAMETERS.  SEE PRINTED MESSAGE.
C         -5    INDEX WAS -1 ON INPUT, BUT THE DESIRED CHANGES OF
C               PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT
C               WAS NOT BEYOND T.  INTERPOLATION TO T = TOUT WAS
C               PERFORMED AS ON A NORMAL RETURN.  TO TRY AGAIN,
C               SIMPLY CALL AGAIN WITH INDEX = -1 AND A NEW TOUT.
C         -6    INSUFFICIENT STORAGE.  THE LENGTH OF THE ARRAY WORK
C               OR IWORK IS LESS THAN THE REQUIRED MINIMUM STORAGE.
C
C   IN ADDITION TO DRIVEP, THE FOLLOWING ROUTINES ARE PROVIDED IN
C   THE PACKAGE..
C     INTERP(TOUT,U,N,U0,WORK,XM,YM)  INTERPOLATES TO GET THE OUTPUT
C                           VALUES AT T = TOUT, FROM THE DATA IN THE U
C                           ARRAY.
C     STIFFP(U,UMAX,..)  IS THE CORE INTEGRATOR ROUTINE.  IT PERFORMS A
C                           SINGLE STEP AND ASSOCIATED ERROR CONTROL.
C     COSET(METH,NQ,EL,TQ,MAXDER)  SETS COEFFICIENTS FOR USE IN
C                           THE CORE INTEGRATOR.
C     DECBR(M0,N,ML,MU,B,IP,IER)  AND  SOLBR(M0,N,ML,MU,B,IP,X)
C                           ARE USED TO SOLVE BANDED LINEAR SYSTEMS,
C                           WITH A BAND MATRIX STORED BY ROWS IN B.
C   NOTE..  DECBR, AND SOLBR ARE CALLED ONLY IF MITER = 1 OR 2.
C
C  THE COMMON BLOCK GEAR3 CAN BE ACCESSED EXTERNALLY BY THE USER
C  IF DESIRED.  IT CONTAINS THE STEP SIZE LAST USED (SUCCESSFULLY),
C  THE ORDER LAST USED (SUCCESSFULLY), THE NUMBER OF STEPS TAKEN
C  SO FAR, THE NUMBER OF F EVALUATIONS SO FAR, AND THE NUMBER OF
C  JACOBIAN EVALUATIONS SO FAR.
C
C  IN THE DATA STATEMENT BELOW, SET..
C    LOUT   =  THE LOGICAL UNIT NUMBER FOR THE OUTPUT OF MESSAGES
C              DURING THE INTEGRATION.
C-----------------------------------------------------------------------
      INTEGER N,MF,INDEX,ML,MU
      INTEGER NC,MFC,KFLAG,JSTART,MLC,MUC,MW,NEBAND,
     1   NEB1,NNEB,NQUSED,NSTEP,NFE,NJE
      INTEGER LOUT,I,NHCUT,KGO
      INTEGER IWORK
      REAL  T0,H0,U0,TOUT,EPS
      REAL  T,H,HMIN,HMAX,EPSC,UROUND,EPSJ,HUSED
      REAL  TOUTP,AYI,D
      REAL  WORK,XM,YM
      COMMON /GEAR1/ T,H,HMIN,HMAX,EPSC,UROUND,NC,MFC,KFLAG,JSTART
      COMMON /GEAR2/ EPSJ,MLC,MUC,MW,NEBAND,NEB1,NNEB
      COMMON /GEAR3/ HUSED,NQUSED,NSTEP,NFE,NJE
      COMMON /IPARAM/ NPDE,NX,NY,MORDER,LWORK,LIWORK
      COMMON /ISTORE/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11,
     *                IW12,IW13,IW14,IW15,IW16,IW17,IW18,IW19,IW20,IW21,
     *                IW22,IW23,IW24,IW25,IW26,IW27
      DIMENSION U0(N),WORK(1),IWORK(1),XM(1),YM(1)
      DATA LOUT/6/
C
      IF (INDEX .EQ. 0) GO TO 20
      IF (INDEX .EQ. 2) GO TO 25
      IF (INDEX .EQ. -1) GO TO 30
      IF (INDEX .EQ. 3) GO TO 40
      IF (INDEX .NE. 1) GO TO 430
      IF (EPS .LE. 0.0) GO TO 400
      IF (N .LE. 0) GO TO 410
      IF ((T0-TOUT)*H0 .GE. 0.0) GO TO 420
C-----------------------------------------------------------------------
C  CALL STRSET TO PROVIDE DYNAMIC DIMENSIONING FOR THE ARRAYS IN THE
C  ROUTINE PDETWO AND THE MODIFIED GEARB PACKAGE.
C-----------------------------------------------------------------------
      CALL STRSET (N,MF,INDEX,IWORK,LOUT,ML,MU)
      IF (INDEX .EQ. -4 .OR. INDEX .EQ. -6) RETURN
C-----------------------------------------------------------------------
C  IF INITIAL VALUES OF UMAX OTHER THAN THOSE SET BELOW ARE DESIRED,
C  THEY SHOULD BE SET HERE.  ALL UMAX(I) MUST BE POSITIVE.
C  IF VALUES FOR HMIN OR HMAX, THE BOUNDS ON ABS(H), OTHER THAN
C  THOSE BELOW ARE DESIRED, THEY SHOULD BE SET BELOW.
C  IN THE FOLLOWING STATEMENT, SET
C    UROUND =  THE UNIT ROUNDOFF OF THE MACHINE, I.E. THE SMALLEST
C              POSITIVE UR SUCH THAT 1. + UR .NE. 1. ON THE MACHINE.
C              (.711E-14 IS THE UNIT ROUNDOFF FOR CDC 6000/7000 SERIES
C              MACHINE).
C-----------------------------------------------------------------------
C
      UROUND = 7.11E-14
      DO 10 I = 1,N
        IM1 = I - 1
        IWA = IW22 + IM1
        IWB = IW23 + IM1
        WORK(IWB) = ABS(U0(I))
        IF(WORK(IWB) .EQ. 0.0) WORK(IWB) = 1.0
 10     WORK(IWA) = U0(I)
      NC = N
      T = T0
      H = H0
      IF ((T+H) .EQ. T) WRITE(LOUT,15)
 15   FORMAT(35H WARNING..  T + H = T ON NEXT STEP.)
      HMIN = ABS(H0)
      HMAX = ABS(T0-TOUT)*10.0
      EPSC = EPS
      MFC = MF
      JSTART = 0
      EPSJ = SQRT(UROUND)
      MLC = ML
      MUC = MU
      MW = ML + MU + 1
      NEBAND = 2*ML + MU + 1
      NEB1 = NEBAND - 1
      NNEB = N*NEBAND
      NHCUT = 0
      GO TO 50
C
C  TOUTP IS THE PREVIOUS VALUE OF TOUT FOR USE IN HMAX.-----------------
 20   HMAX = ABS(TOUT-TOUTP)*10.0
      GO TO 80
C
 25   HMAX = ABS(TOUT-TOUTP)*10.0
      IF ((T-TOUT)*H .GE. 0.0) GO TO 500
      GO TO 85
C
 30   IF ((T-TOUT)*H .GE. 0.0) GO TO 440
      JSTART = -1
      NC = N
      EPSC = EPS
      MFC = MF
C
 40   IF ((T+H) .EQ. T) WRITE(LOUT,15)
C
 50   CALL STIFFP (WORK(IW22),WORK(IW23),WORK(IW24),WORK(IW25),
     *             WORK(IW26),WORK(IW27),N,WORK,IWORK,XM,YM)
C
      KGO = 1 - KFLAG
      GO TO (60, 100, 200, 300), KGO
C  KFLAG  =   0,  -1,  -2,  -3
C
 60   CONTINUE
C-----------------------------------------------------------------------
C  NORMAL RETURN FROM INTEGRATOR.
C
C  THE WEIGHTS UMAX(I) ARE UPDATED.  IF DIFFERENT VALUES ARE DESIRED,
C  THEY SHOULD BE SET HERE.  A TEST IS MADE FOR EPS BEING TOO SMALL
C  FOR THE MACHINE PRECISION.
C
C  ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EVERY
C  STEP SHOULD BE INSERTED HERE.
C
C  IF INDEX = 3, U0 IS SET TO THE CURRENT U VALUES ON RETURN.
C  IF INDEX = 2, H IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF
C  ERROR), AND THEN THE CURRENT U VALUES ARE PUT IN U0 ON RETURN.
C  FOR ANY OTHER VALUE OF INDEX, CONTROL RETURNS TO THE INTEGRATOR
C  UNLESS TOUT HAS BEEN REACHED.  THEN INTERPOLATED VALUES OF U ARE
C  COMPUTED AND STORED IN U0 ON RETURN.
C  IF INTERPOLATION IS NOT DESIRED, THE CALL TO INTERP SHOULD BE
C  REMOVED AND CONTROL TRANSFERRED TO STATEMENT 500 INSTEAD OF 520.
C-----------------------------------------------------------------------
      D = 0.0
      DO 70 I = 1,N
        IM1 = I - 1
        IWA = IW22 + IM1
        IWB = IW23 + IM1
        AYI = ABS(WORK(IWA))
        WORK(IWB) = AMAX1(WORK(IWB), AYI)
 70     D = D + (AYI/WORK(IWB))**2
      D = D*(UROUND/EPS)**2
      IF (D .GT. FLOAT(N)) GO TO 250
      IF (INDEX .EQ. 3) GO TO 500
      IF (INDEX .EQ. 2) GO TO 85
 80   IF ((T-TOUT)*H .LT. 0.0) GO TO 40
      CALL INTERP (TOUT,WORK(IW22),N,U0,WORK,XM,YM)
      GO TO 520
 85   IF (((T+H)-TOUT)*H .LE. 0.0) GO TO 40
      IF (ABS(T-TOUT) .LE. 100.0*UROUND*HMAX) GO TO 500
      IF ((T-TOUT)*H .GE. 0.0) GO TO 500
      H = (TOUT - T)*(1.0 - 4.0*UROUND)
      JSTART = -1
      GO TO 40
C-----------------------------------------------------------------------
C  ON AN ERROR RETURN FROM INTEGRATOR, AN IMMEDIATE RETURN OCCURS IF
C  KFLAG = -2, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE.
C  TO RECOVER, H AND HMIN ARE REDUCED BY A FACTOR OF .1 UP TO 10
C  TIMES BEFORE GIVING UP.
C-----------------------------------------------------------------------
 100  WRITE (LOUT,105) T
 105  FORMAT(//35H KFLAG = -1 FROM INTEGRATOR AT T = ,E16.8/
     1   38H  ERROR TEST FAILED WITH ABS(H) = HMIN/)
 110  IF (NHCUT .EQ. 10) GO TO 150
      NHCUT = NHCUT + 1
      HMIN = .10*HMIN
      H = .10*H
      WRITE (LOUT,115) H
 115  FORMAT(24H  H HAS BEEN REDUCED TO ,E16.8,
     1   26H  AND STEP WILL BE RETRIED//)
      JSTART = -1
      GO TO 40
C
 150  WRITE (LOUT,155)
 155  FORMAT(//44H PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT//)
      GO TO 500
C
 200  WRITE (LOUT,205) T,H
 205  FORMAT(//35H KFLAG = -2 FROM INTEGRATOR AT T = ,E16.8,5H  H =,
     1  E16.8/52H  THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED//)
      GO TO 500
C
 250  WRITE (LOUT,255) T
 255  FORMAT(//37H INTEGRATION HALTED BY DRIVER AT T = ,E16.8/
     1   56H  EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION/)
      KFLAG = -2
      GO TO 500
C
 300  WRITE (LOUT,305) T
 305  FORMAT(//35H KFLAG = -3 FROM INTEGRATOR AT T = ,E16.8/
     1   45H  CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED/)
      GO TO 110
C
 400  WRITE (LOUT,405)
 405  FORMAT(//28H ILLEGAL INPUT.. EPS .LE. 0.//)
      INDEX = -4
      RETURN
C
 410  WRITE (LOUT,415)
 415  FORMAT(//25H ILLEGAL INPUT.. N .LE. 0//)
      INDEX = -4
      RETURN
C
 420  WRITE (LOUT,425)
 425  FORMAT(//36H ILLEGAL INPUT.. (T0-TOUT)*H .GE. 0.//)
      INDEX = -4
      RETURN
C
 430  WRITE (LOUT,435) INDEX
 435  FORMAT(//24H ILLEGAL INPUT.. INDEX =,I5//)
      INDEX = -4
      RETURN
C
 440  WRITE(LOUT,445) T,TOUT,H
 445  FORMAT(//44H INDEX = -1 ON INPUT WITH (T-TOUT)*H .GE. 0./
     1   4H T =,E16.8,9H   TOUT =,E16.8,6H   H =,E16.8/
     1   44H INTERPOLATION WAS DONE AS ON NORMAL RETURN./
     2   41H DESIRED PARAMETER CHANGES WERE NOT MADE.)
      CALL INTERP (TOUT,WORK(IW22),N,U0,WORK,XM,YM)
      INDEX = -5
      RETURN
C
 500  TOUT = T
      DO 510 I = 1,N
        IM1 = I - 1
        IWA = IW22 + IM1
 510    U0(I) = WORK(IWA)
 520  INDEX = KFLAG
      TOUTP = TOUT
      H0 = HUSED
      IF (KFLAG .NE. 0) H0 = H
      RETURN
C----------------------- END OF SUBROUTINE DRIVEP ----------------------
      END
      SUBROUTINE INTERP (TOUT,U,N,U0,WORK,XM,YM)                         INTERP
      INTEGER N,IDUMMY,JSTART,I,L,J
      REAL  TOUT,U,U0,T,H,DUMMY,S,S1
      DIMENSION U0(N),U(N,1),WORK(1),XM(1),YM(1)
      COMMON /GEAR1/ T,H,DUMMY(4),IDUMMY(3),JSTART
      COMMON /IPARAM/ NPDE,NX,NY,MORDER,LWORK,LIWORK
      COMMON /ISTORE/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11,
     *                IW12,IW13,IW14,IW15,IW16,IW17,IW18,IW19,IW20,IW21,
     *                IW22,IW23,IW24,IW25,IW26,IW27
      LOGICAL  IX,JY
C-----------------------------------------------------------------------
C  SUBROUTINE INTERP COMPUTES INTERPOLATED VALUES OF THE DEPENDENT
C  VARIABLE U AND STORES THEM IN U0.  THE INTERPOLATION IS TO THE
C  POINT T = TOUT, AND USES THE NORDSIECK HISTORY ARRAY U, AS FOLLOWS..
C                              NQ
C                   U0(I)  =  SUM  U(I,J+1)*S**J ,
C                             J=0
C  WHERE S = -(T-TOUT)/H.  FOR THE PACKAGE PDETWO, IT IS MODIFIED TO
C  DETERMINE ALSO THE SOLUTIONS FOR CONSTANT, TIME DEPENDENT BOUNDARY
C  CONDITIONS.
C-----------------------------------------------------------------------
      DO 10 I = 1,N
 10     U0(I) = U(I,1)
      L = JSTART + 1
      S = (TOUT - T)/H
      S1 = 1.0
      DO 30 J = 2,L
        S1 = S1*S
        DO 20 I = 1,N
 20       U0(I) = U0(I) + S1*U(I,J)
 30     CONTINUE
C
C  DETERMINE THE SOLUTIONS FOR CONSTANT, TIME DEPENDENT BOUNDARY
C  CONDITIONS
C
      DO 60 J=1,NY
        JNX = (J-1)*NX
        JY = J.EQ.1.OR.J.EQ.NY
        DO 60 I=1,NX
          IJNPDE = (I+JNX-1)*NPDE
          IJ1 = 1 + IJNPDE
          IX = I.EQ.1.OR.I.EQ.NX
          IF (JY) CALL BNDRYH (TOUT,XM(I),YM(J),U0(IJ1),WORK(IW4),
     *                         WORK(IW5),WORK(IW6),NPDE)
          IF (IX) CALL BNDRYV (TOUT,XM(I),YM(J),U0(IJ1),WORK(IW7),
     *                         WORK(IW8),WORK(IW9),NPDE)
          DO 50 K=1,NPDE
            KIJ = K + IJNPDE
            KM1 = K - 1
            KIW4 = IW4 + KM1
            KIW5 = IW5 + KM1
            KIW6 = IW6 + KM1
            KIW7 = IW7 + KM1
            KIW8 = IW8 + KM1
            KIW9 = IW9 + KM1
            IF (IX.AND.JY) GO TO 40
            IF (JY.AND.WORK(KIW5).EQ.0.0) U0(KIJ) = WORK(KIW6)/
     *                                              WORK(KIW4)
            IF (IX.AND.WORK(KIW8).EQ.0.0) U0(KIJ) = WORK(KIW9)/
     *                                              WORK(KIW7)
            GO TO 50
 40         IF (WORK(KIW5).EQ.0.0.AND.WORK(KIW8).NE.0.0) U0(KIJ) =
     *        WORK(KIW6)/WORK(KIW4)
            IF (WORK(KIW5).NE.0.0.AND.WORK(KIW8).EQ.0.0) U0(KIJ) =
     *        WORK(KIW9)/WORK(KIW7)
            IF (WORK(KIW5).EQ.0.0.AND.WORK(KIW8).EQ.0.0) U0(KIJ) = 0.5*
     *        ((WORK(KIW6)/WORK(KIW4)) + (WORK(KIW9)/WORK(KIW7)))
 50       CONTINUE
 60   CONTINUE
      RETURN
C----------------------- END OF SUBROUTINE INTERP ----------------------
      END
      SUBROUTINE STIFFP (U,UMAX,ERROR,SAVE1,SAVE2,PW,N,WORK,IWORK,       STIFFP
     *                   XM,YM)
      INTEGER N,MF,KFLAG,JSTART,ML,MU,MW,NEBAND,
     1   NEB1,NNEB,NQUSED,NSTEP,NFE,NJE
      INTEGER I,METH,MITER,NQ,L,IDOUB,MFOLD,NOLD,IRET,MEO,MIO,
     1   IWEVAL,MAXDER,LMAX,IREDO,J,NSTEPJ,J1,J2,M,IER,NEWQ
      INTEGER IWORK
      REAL  U,T,H,HMIN,HMAX,EPS,UROUND,UMAX,ERROR,
     1   SAVE1,SAVE2,PW,EPSJ,HUSED
      REAL  EL,OLDL0,TOLD,RMAX,RC,CRATE,EPSOLD,
     1   HOLD,FN,EDN,E,EUP,BND,RH,R1,CON,R,HL0,R0,D,PHL0,
     1   PR3,D1,ENQ3,ENQ2,PR2,PR1,ENQ1
      REAL  WORK,XM,YM
      REAL  TQ
      COMMON /GEAR1/ T,H,HMIN,HMAX,EPS,UROUND,NC,MF,KFLAG,JSTART
      COMMON /GEAR2/ EPSJ,ML,MU,MW,NEBAND,NEB1,NNEB
      COMMON /GEAR3/ HUSED,NQUSED,NSTEP,NFE,NJE
      COMMON /IPARAM/ NPDE,NX,NY,MORDER,LWORK,LIWORK
      COMMON /ISTORE/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11,
     *                IW12,IW13,IW14,IW15,IW16,IW17,IW18,IW19,IW20,IW21,
     *                IW22,IW23,IW24,IW25,IW26,IW27
      DIMENSION U(N,1),UMAX(N),ERROR(N),SAVE1(N),SAVE2(N),PW(1),
     *          WORK(1),IWORK(1),XM(1),YM(1)
C-----------------------------------------------------------------------
C  STIFFP PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE
C  PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.
C  STIFFP IS A VERSION FOR BANDED FORM OF THE JACOBIAN MATRIX.
C  COMMUNICATION WITH STIFFP IS DONE WITH THE FOLLOWING VARIABLES..
C
C    U       AN N BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES
C              AND THEIR SCALED DERIVATIVES.  LMAX IS 13 FOR THE ADAMS
C              METHODS AND 6 FOR THE GEAR METHODS.  LMAX - 1 = MAXDER
C              IS THE MAXIMUM ORDER AVAILABLE.  SEE SUBROUTINE COSET.
C              U(I,J+1) CONTAINS THE J-TH DERIVATIVE OF U(I), SCALED BY
C              H**J/FACTORIAL(J)  (J = 0,1,...,NQ).
C    UMAX    AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL
C              ERRORS IN U ARE COMPARED.
C    ERROR   AN ARRAY OF N ELEMENTS.  ERROR(I)/TQ(2) IS THE ESTIMATED
C              ONE-STEP ERROR IN U(I).
C    SAVE1,  TWO ARRAYS OF WORKING STORAGE,
C     SAVE2    EACH OF LENGTH N.
C    PW      A BLOCK OF LOCATIONS USED FOR PARTIAL DERIVATIVES IF
C              MITER IS NOT 0.  SEE DESCRIPTION IN DRIVER.
C    N       THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.
C    WORK    FLOATING POINT WORKING ARRAY. SEE STRSET FOR ADDITIONAL
C              DETAILS.
C    IWORK   AN INTEGER ARRAY OF LENGTH N USED FOR PIVOT
C              INFORMATION IF MITER = 1 OR 2.
C    XM      ARE THE MESH POINTS IN THE HORIZONTAL DIRECTION.
C    YM      ARE THE MESH POINTS IN THE VERTICAL DIRECTION.
C    T       THE INDEPENDENT VARIABLE. T IS UPDATED ON EACH STEP TAKEN.
C    H       THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
C              H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE
C              PROBLEM.  H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS
C              SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM.
C    HMIN,   THE MINIMUM AND MAXIMUM ABSOLUTE VALUE OF THE STEP SIZE
C     HMAX     TO BE USED FOR THE STEP.  THESE MAY BE CHANGED AT ANY
C              TIME, BUT WILL NOT TAKE EFFECT UNTIL THE NEXT H CHANGE.
C    EPS     THE RELATIVE ERROR BOUND.  SEE DESCRIPTION IN DRIVER.
C    UROUND  THE UNIT ROUNDOFF OF THE MACHINE.
C    MF      THE METHOD FLAG.  SEE DESCRIPTION IN DRIVER.
C    KFLAG   A COMPLETION CODE WITH THE FOLLOWING MEANINGS..
C                      0  THE STEP WAS SUCCESFUL.
C                     -1  THE REQUESTED ERROR COULD NOT BE ACHIEVED
C                           WITH ABS(H) = HMIN.
C                     -2  THE REQUESTED ERROR IS SMALLER THAN CAN
C                           BE HANDLED FOR THIS PROBLEM.
C                     -3  CORRECTOR CONVERGENCE COULD NOT BE
C                           ACHIEVED FOR ABS(H) = HMIN.
C              ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND
C              THE U ARRAY ARE AS OF THE BEGINNING OF THE LAST
C              STEP, AND H IS THE LAST STEP SIZE ATTEMPTED.
C    JSTART  AN INTEGER USED ON INPUT AND OUTPUT.
C              ON INPUT, IT HAS THE FOLLOWING VALUES AND MEANINGS..
C                      0  PERFORM THE FIRST STEP.
C                  .GT.0  TAKE A NEW STEP CONTINUING FROM THE LAST.
C                  .LT.0  TAKE THE NEXT STEP WITH A NEW VALUE OF
C                           H, EPS, N, AND/OR MF.
C              ON EXIT, JSTART IS NQ, THE CURRENT ORDER OF THE METHOD.
C    ML,MU   THE LOWER AND UPPER HALF BANDWIDTHS, RESPECTIVELY, OF
C              THE JACOBIAN.  SEE DESCRIPTION IN DRIVER.
C    NEBAND  THE EXTENDED BANDWIDTH, 2*ML + MU + 1.
C    NNEB    THE MAXIMUM SIZE OF PW USED, N*NEBAND.
C-----------------------------------------------------------------------
      DIMENSION EL(13),TQ(4)
      DATA EL(2)/1.0/, OLDL0/1.0/
      KFLAG = 0
      TOLD = T
      IF (JSTART .GT. 0) GO TO 200
      IF (JSTART .NE. 0) GO TO 120
C-----------------------------------------------------------------------
C  ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL UDOT IS
C  CALCULATED.  RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED
C  IN A SINGLE STEP.  IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL
C  INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10.  IF A FAILURE
C  OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2
C  FOR THE NEXT INCREASE.
C-----------------------------------------------------------------------
      CALL PDETWO (NPDE,NX,NY,XM,YM, T, U, SAVE1, WORK,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))
      DO 110 I = 1,N
 110    U(I,2) = H*SAVE1(I)
      METH = MF/10
      MITER = MF - 10*METH
      NQ = 1
      L = 2
      IDOUB = 3
      RMAX = 1.E+4
      RC = 0.0
      CRATE = 1.0
      EPSOLD = EPS
      HOLD = H
      MFOLD = MF
      NOLD = N
      NSTEP = 0
      NSTEPJ = 0
      NFE = 1
      NJE = 0
      IRET = 3
      IF (MITER .NE. 1) GO TO 130
      DO 115 I = 1,NNEB
 115    PW(I) = 0.0
      GO TO 130
C-----------------------------------------------------------------------
C  IF THE CALLER HAS CHANGED METH, COSET IS CALLED TO SET
C  THE COEFFICIENTS OF THE METHOD.  IF THE CALLER HAS CHANGED
C  N, EPS, OR METH, THE CONSTANTS E, EDN, EUP, AND BND MUST BE RESET.
C  E IS A COMPARISON FOR ERRORS OF THE CURRENT ORDER NQ. EUP IS
C  TO TEST FOR INCREASING THE ORDER, EDN FOR DECREASING THE ORDER.
C  BND IS USED TO TEST FOR CONVERGENCE OF THE CORRECTOR ITERATES.
C  IF THE CALLER HAS CHANGED H, U MUST BE RESCALED.
C  IF H OR METH HAS BEEN CHANGED, IDOUB IS RESET TO L + 1 TO PREVENT
C  FURTHER CHANGES IN H FOR THAT MANY STEPS.
C-----------------------------------------------------------------------
 120  IF (MF .EQ. MFOLD) GO TO 150
      MEO = METH
      MIO = MITER
      METH = MF/10
      MITER = MF - 10*METH
      MFOLD = MF
      IF (MITER .NE. MIO) IWEVAL = MITER
      IF (METH .EQ. MEO) GO TO 150
      IDOUB = L + 1
      IRET = 1
 130  CALL COSET (METH, NQ, EL, TQ, MAXDER)
      MAXDER = MIN0(MAXDER,MORDER)
      LMAX = MAXDER + 1
      RC = RC*EL(1)/OLDL0
      OLDL0 = EL(1)
 140  FN = FLOAT(N)
      EDN = FN*(TQ(1)*EPS)**2
      E   = FN*(TQ(2)*EPS)**2
      EUP = FN*(TQ(3)*EPS)**2
      BND = FN*(TQ(4)*EPS)**2
      GO TO (160, 170, 200), IRET
 150  IF ((EPS .EQ. EPSOLD) .AND. (N .EQ. NOLD)) GO TO 160
      EPSOLD = EPS
      NOLD = N
      IRET = 1
      GO TO 140
 160  IF (H .EQ. HOLD) GO TO 200
      RH = H/HOLD
      H = HOLD
      IREDO = 3
      GO TO 175
 170  RH = AMAX1(RH,HMIN/ABS(H))
 175  RH = AMIN1(RH,HMAX/ABS(H),RMAX)
      R1 = 1.0
      DO 180 J = 2,L
        R1 = R1*RH
        DO 180 I = 1,N
 180      U(I,J) = U(I,J)*R1
      H = H*RH
      RC = RC*RH
      IDOUB = L + 1
      IF (IREDO .EQ. 0) GO TO 690
C-----------------------------------------------------------------------
C  THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY
C  MULTIPLYING THE U ARRAY BY THE PASCAL TRIANGLE MATRIX.
C  RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT  H*EL(1).
C  WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, OR THE CALLER HAS
C  CHANGED MITER, IWEVAL IS SET TO MITER TO FORCE THE PARTIALS TO BE
C  UPDATED, IF PARTIALS ARE USED.  IN ANY CASE, THE PARTIALS
C  ARE UPDATED AT LEAST EVERY 20-TH STEP.
C-----------------------------------------------------------------------
 200  IF (ABS(RC-1.0) .GT. 0.30) IWEVAL = MITER
      IF (NSTEP .GE. NSTEPJ+20) IWEVAL = MITER
      T = T + H
      DO 210 J1 = 1,NQ
        DO 210 J2 = J1,NQ
          J = (NQ + J1) - J2
          DO 210 I = 1,N
 210        U(I,J) = U(I,J) + U(I,J+1)
C-----------------------------------------------------------------------
C  UP TO 3 CORRECTOR ITERATIONS ARE TAKEN.  A CONVERGENCE TEST IS
C  MADE ON THE R.M.S. NORM OF EACH CORRECTION, USING BND, WHICH
C  IS DEPENDENT ON EPS.  THE SUM OF THE CORRECTIONS IS ACCUMULATED
C  IN THE VECTOR ERROR(I).  THE U ARRAY IS NOT ALTERED IN THE CORRECTOR
C  LOOP.  THE UPDATED U VECTOR IS STORED TEMPORARILY IN SAVE1.
C-----------------------------------------------------------------------
 220  DO 230 I = 1,N
 230    ERROR(I) = 0.0
      M = 0
      CALL PDETWO (NPDE,NX,NY,XM,YM, T, U, SAVE2, WORK,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))
      NFE = NFE + 1
      IF (IWEVAL .LE. 0) GO TO 290
C-----------------------------------------------------------------------
C  IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED BEFORE
C  STARTING THE CORRECTOR ITERATION.  IWEVAL IS SET TO 0 AS AN
C  INDICATOR THAT THIS HAS BEEN DONE.  IF MITER = 1 OR 2, P IS
C  COMPUTED AND PROCESSED IN PSETB.  IF MITER = 3, THE MATRIX USED
C  IS P = I - H*EL(1)*D, WHERE D IS A DIAGONAL MATRIX.
C-----------------------------------------------------------------------
      IWEVAL = 0
      RC = 1.0
      NJE = NJE + 1
      NSTEPJ = NSTEP
      GO TO (250, 240, 260), MITER
C240  NFE = NFE + MW
 240  NFE = NFE + 5 * NPDE
 250  CON = -H*EL(1)
      CALL PSETM (NPDE,NX,NY,XM,YM,U,UMAX,ERROR,SAVE1,SAVE2,PW,CON,
     *            MITER,IER,NEBAND,WORK,IWORK)
      IF (IER .NE. 0) GO TO 420
      GO TO 350
 260  R = EL(1)*.10
      DO 270 I = 1,N
 270    PW(I) = U(I,1) + R*(H*SAVE2(I) - U(I,2))
      CALL PDETWO (NPDE,NX,NY,XM,YM, T, PW, SAVE1, WORK,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))
      NFE = NFE + 1
      HL0 = H*EL(1)
      DO 280 I = 1,N
        R0 = H*SAVE2(I) - U(I,2)
        PW(I) = 1.0
        D = .10*R0 - H*(SAVE1(I) - SAVE2(I))
        SAVE1(I) = 0.0
        IF (ABS(R0) .LT. UROUND*UMAX(I)) GO TO 280
        IF (ABS(D) .EQ. 0.0) GO TO 420
        PW(I) = .10*R0/D
        SAVE1(I) = PW(I)*R0
 280    CONTINUE
      GO TO 370
 290  IF (MITER .NE. 0) GO TO (350, 350, 310), MITER
C-----------------------------------------------------------------------
C  IN THE CASE OF FUNCTIONAL ITERATION, UPDATE U DIRECTLY FROM
C  THE RESULT OF THE LAST DIFFUN CALL.
C-----------------------------------------------------------------------
      D = 0.0
      DO 300 I = 1,N
        R = H*SAVE2(I) - U(I,2)
        D = D + ( (R-ERROR(I))/UMAX(I) )**2
        SAVE1(I) = U(I,1) + EL(1)*R
 300    ERROR(I) = R
      GO TO 400
C-----------------------------------------------------------------------
C  IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR,
C  F SUB (M), AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND
C  SIDE AND P AS COEFFICIENT MATRIX, USING THE LU DECOMPOSITION
C  IF MITER = 1 OR 2.  IF MITER = 3, THE COEFFICIENT H*EL(1)
C  IN P IS UPDATED.
C-----------------------------------------------------------------------
 310  PHL0 = HL0
      HL0 = H*EL(1)
      IF (HL0 .EQ. PHL0) GO TO 330
      R = HL0/PHL0
      DO 320 I = 1,N
        D = 1.0 - R*(1.0 - 1.0/PW(I))
        IF (ABS(D) .EQ. 0.0) GO TO 440
 320    PW(I) = 1.0/D
 330  DO 340 I = 1,N
 340    SAVE1(I) = PW(I)*(H*SAVE2(I) - (U(I,2) + ERROR(I)))
      GO TO 370
 350  DO 360 I = 1,N
 360    SAVE1(I) = H*SAVE2(I) - (U(I,2) + ERROR(I))
      CALL SOLBR (NEBAND,N,ML,MU,PW,IWORK,SAVE1)
 370  D = 0.0
      DO 380 I = 1,N
        ERROR(I) = ERROR(I) + SAVE1(I)
        D = D + (SAVE1(I)/UMAX(I))**2
 380    SAVE1(I) = U(I,1) + EL(1)*ERROR(I)
C-----------------------------------------------------------------------
C  TEST FOR CONVERGENCE.  IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE
C  RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST.
C-----------------------------------------------------------------------
 400  IF (M .NE. 0) CRATE = AMAX1(.90*CRATE,D/D1)
      IF ((D*AMIN1(1.0,2.0*CRATE)) .LE. BND) GO TO 450
      D1 = D
      M = M + 1
      IF (M .EQ. 3) GO TO 410
      CALL PDETWO (NPDE,NX,NY,XM,YM, T, SAVE1, SAVE2, WORK,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))
      GO TO 290
C-----------------------------------------------------------------------
C  THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES.  IF PARTIALS
C  ARE INVOLVED BUT ARE NOT UP TO DATE, THEY ARE REEVALUATED FOR THE
C  NEXT TRY.  OTHERWISE THE U ARRAY IS RETRACTED TO ITS VALUES
C  BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE.  IF NOT, A
C  NO-CONVERGENCE EXIT IS TAKEN.
C-----------------------------------------------------------------------
 410  NFE = NFE + 2
      IF (IWEVAL .EQ. -1) GO TO 440
 420  T = TOLD
      RMAX = 2.0
      DO 430 J1 = 1,NQ
        DO 430 J2 = J1,NQ
          J = (NQ + J1) - J2
          DO 430 I = 1,N
 430        U(I,J) = U(I,J) - U(I,J+1)
      IF (ABS(H) .LE. HMIN*1.00001) GO TO 680
      RH = .250
      IREDO = 1
      GO TO 170
 440  IWEVAL = MITER
      GO TO 220
C-----------------------------------------------------------------------
C  THE CORRECTOR HAS CONVERGED.  IWEVAL IS SET TO -1 IF PARTIAL
C  DERIVATIVES WERE USED, TO SIGNAL THAT THEY MAY NEED UPDATING ON
C  SUBSEQUENT STEPS.  THE ERROR TEST IS MADE AND CONTROL PASSES TO
C  STATEMENT 500 IF IT FAILS.
C-----------------------------------------------------------------------
 450  IF (MITER .NE. 0) IWEVAL = -1
      NFE = NFE + M
      D = 0.0
      DO 460 I = 1,N
 460    D = D + (ERROR(I)/UMAX(I))**2
      IF (D .GT. E) GO TO 500
C-----------------------------------------------------------------------
C  AFTER A SUCCESSFUL STEP, UPDATE THE U ARRAY.
C  CONSIDER CHANGING H IF IDOUB = 1.  OTHERWISE DECREASE IDOUB BY 1.
C  IF IDOUB IS THEN 1 AND NQ .LT. MAXDER, THEN ERROR IS SAVED FOR
C  USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
C  IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER
C  BY ONE IS CONSIDERED ALSO.  A CHANGE IN H IS MADE ONLY IF IT IS BY A
C  FACTOR OF AT LEAST 1.1.  IF NOT, IDOUB IS SET TO 10 TO PREVENT
C  TESTING FOR THAT MANY STEPS.
C-----------------------------------------------------------------------
      KFLAG = 0
      IREDO = 0
      NSTEP = NSTEP + 1
      HUSED = H
      NQUSED = NQ
      DO 470 J = 1,L
        DO 470 I = 1,N
 470      U(I,J) = U(I,J) + EL(J)*ERROR(I)
      IF (IDOUB .EQ. 1) GO TO 520
      IDOUB = IDOUB - 1
      IF (IDOUB .GT. 1) GO TO 700
      IF (L .EQ. LMAX) GO TO 700
      DO 490 I = 1,N
 490    U(I,LMAX) = ERROR(I)
      GO TO 700
C-----------------------------------------------------------------------
C  THE ERROR TEST FAILED.  KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
C  RESTORE T AND THE U ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE
C  TO TRY THE STEP AGAIN.  COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR
C  ONE LOWER ORDER.
C-----------------------------------------------------------------------
 500  KFLAG = KFLAG - 1
      T = TOLD
      DO 510 J1 = 1,NQ
        DO 510 J2 = J1,NQ
          J = (NQ + J1) - J2
          DO 510 I = 1,N
 510        U(I,J) = U(I,J) - U(I,J+1)
      RMAX = 2.0
      IF (ABS(H) .LE. HMIN*1.000010) GO TO 660
      IF (KFLAG .LE. -3) GO TO 640
      IREDO = 2
      PR3 = 1.E+20
      GO TO 540
C-----------------------------------------------------------------------
C  REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS
C  PR1, PR2, AND PR3 ARE COMPUTED, BY WHICH H COULD BE DIVIDED
C  AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY.
C  IN THE CASE OF FAILURE, PR3 = 1.E20 TO AVOID AN ORDER INCREASE.
C  THE SMALLEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN
C  ACCORDINGLY.  IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE
C  ADDITIONAL SCALED DERIVATIVE.
C-----------------------------------------------------------------------
 520  PR3 = 1.E+20
      IF (L .EQ. LMAX) GO TO 540
      D1 = 0.0
      DO 530 I = 1,N
 530    D1 = D1 + ((ERROR(I) - U(I,LMAX))/UMAX(I))**2
      ENQ3 = .50/FLOAT(L+1)
      PR3 = ((D1/EUP)**ENQ3)*1.40 + 1.4E-6
 540  ENQ2 = .50/FLOAT(L)
      PR2 = ((D/E)**ENQ2)*1.20 + 1.2E-6
      PR1 = 1.E+20
      IF (NQ .EQ. 1) GO TO 560
      D = 0.0
      DO 550 I = 1,N
 550    D = D + (U(I,L)/UMAX(I))**2
      ENQ1 = .50/FLOAT(NQ)
      PR1 = ((D/EDN)**ENQ1)*1.30 + 1.3E-6
 560  IF (PR2 .LE. PR3) GO TO 570
      IF (PR3 .LT. PR1) GO TO 590
      GO TO 580
 570  IF (PR2 .GT. PR1) GO TO 580
      NEWQ = NQ
      RH = 1.0/PR2
      GO TO 620
 580  NEWQ = NQ - 1
      RH = 1.0/PR1
      GO TO 620
 590  NEWQ = L
      RH = 1.0/PR3
      IF (RH .LT. 1.10) GO TO 610
      DO 600 I = 1,N
 600    U(I,NEWQ+1) = ERROR(I)*EL(L)/FLOAT(L)
      GO TO 630
 610  IDOUB = 10
      GO TO 700
 620  IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.10)) GO TO 610
C-----------------------------------------------------------------------
C  IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS.
C  IN ANY CASE H IS RESET ACCORDING TO RH AND THE U ARRAY IS RESCALED.
C  THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE.
C-----------------------------------------------------------------------
      IF (NEWQ .EQ. NQ) GO TO 170
 630  NQ = NEWQ
      L = NQ + 1
      IRET = 2
      GO TO 130
C-----------------------------------------------------------------------
C  CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED.
C  IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE
C  U ARRAY HAVE ERRORS OF THE WRONG ORDER.  HENCE THE FIRST
C  DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1.  THEN
C  H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED.
C  AFTER A TOTAL OF 7 FAILURES, AN EXIT IS TAKEN WITH KFLAG = -2.
C-----------------------------------------------------------------------
 640  IF (KFLAG .EQ. -7) GO TO 670
      RH = .10
      RH = AMAX1(HMIN/ABS(H),RH)
      H = H*RH
      CALL PDETWO (NPDE,NX,NY,XM,YM, T, U, SAVE1, WORK,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))
      NFE = NFE + 1
      DO 650 I = 1,N
 650    U(I,2) = H*SAVE1(I)
      IWEVAL = MITER
      IDOUB = 10
      IF (NQ .EQ. 1) GO TO 200
      NQ = 1
      L = 2
      IRET = 3
      GO TO 130
C-----------------------------------------------------------------------
C  ALL RETURNS ARE MADE THROUGH THIS SECTION.  H IS SAVED IN HOLD
C  TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
C-----------------------------------------------------------------------
 660  KFLAG = -1
      GO TO 700
 670  KFLAG = -2
      GO TO 700
 680  KFLAG = -3
      GO TO 700
 690  RMAX = 10.0
 700  HOLD = H
      JSTART = NQ
      RETURN
C----------------------- END OF SUBROUTINE STIFFP ----------------------
      END
      SUBROUTINE COSET (METH, NQ, EL, TQ, MAXDER)                        COSET
C
      INTEGER METH,NQ,MAXDER,K
      REAL  EL
      REAL  TQ,PERTST
C-----------------------------------------------------------------------
C  COSET IS CALLED BY THE INTEGRATOR AND SETS COEFFICIENTS USED THERE.
C  THE VECTOR EL, OF LENGTH NQ + 1, DETERMINES THE BASIC METHOD.
C  THE VECTOR TQ, OF LENGTH 4, IS INVOLVED IN ADJUSTING THE STEP SIZE
C  IN RELATION TO TRUNCATION ERROR.  ITS VALUES ARE GIVEN BY THE
C  PERTST ARRAY.
C  THE VECTORS EL AND TQ DEPEND ON METH AND NQ.
C  COSET ALSO SETS MAXDER, THE MAXIMUM ORDER OF THE METHOD AVAILABLE.
C  CURRENTLY IT IS 12 FOR THE ADAMS METHODS AND 5 FOR THE BDF METHODS.
C  LMAX = MAXDER + 1 IS THE NUMBER OF COLUMNS IN THE U ARRAY.
C  THE MAXIMUM ORDER USED MAY BE REDUCED SIMPLY BY DECREASING THE
C  NUMBERS IN STATEMENTS 1 AND/OR 2 BELOW.
C
C  THE COEFFICIENTS IN PERTST NEED BE GIVEN TO ONLY ABOUT
C  ONE PERCENT ACCURACY.  THE ORDER IN WHICH THE GROUPS APPEAR BELOW
C  IS..  COEFFICIENTS FOR ORDER NQ - 1, COEFFICIENTS FOR ORDER NQ,
C  COEFFICIENTS FOR ORDER NQ + 1.  WITHIN EACH GROUP ARE THE
C  COEFFICIENTS FOR THE ADAMS METHODS, FOLLOWED BY THOSE FOR THE
C  BDF METHODS.
C-----------------------------------------------------------------------
      DIMENSION PERTST(12,2,3),EL(13),TQ(4)
      DATA PERTST(1,1,1),PERTST(2,1,1),PERTST(3,1,1),PERTST(4,1,1)
     1           /1., 1., 2., 1./,
     2     PERTST(5,1,1),PERTST(6,1,1),PERTST(7,1,1),PERTST(8,1,1)
     3           /.3158, .07407, .01391, .002182/,
     4     PERTST(9,1,1),PERTST(10,1,1),PERTST(11,1,1),PERTST(12,1,1)
     5           /.0002945, .00003492, .000003692, .0000003524/
      DATA PERTST(1,2,1),PERTST(2,2,1),PERTST(3,2,1),PERTST(4,2,1)
     1           /1., 1., .5, .1667/,
     2     PERTST(5,2,1),PERTST(6,2,1),PERTST(7,2,1),PERTST(8,2,1)
     3           /.04167, 1., 1., 1./,
     4     PERTST(9,2,1),PERTST(10,2,1),PERTST(11,2,1),PERTST(12,2,1)
     5           /1., 1., 1., 1./
      DATA PERTST(1,1,2),PERTST(2,1,2),PERTST(3,1,2),PERTST(4,1,2)
     1           /2., 12., 24., 37.89/,
     2     PERTST(5,1,2),PERTST(6,1,2),PERTST(7,1,2),PERTST(8,1,2)
     3           /53.33, 70.08, 87.97, 106.9/,
     4     PERTST(9,1,2),PERTST(10,1,2),PERTST(11,1,2),PERTST(12,1,2)
     5           /126.7, 147.4, 168.8, 191./
      DATA PERTST(1,2,2),PERTST(2,2,2),PERTST(3,2,2),PERTST(4,2,2)
     1           /2., 4.5, 7.333, 10.42/,
     2     PERTST(5,2,2),PERTST(6,2,2),PERTST(7,2,2),PERTST(8,2,2)
     3           /13.7, 1., 1., 1./,
     4     PERTST(9,2,2),PERTST(10,2,2),PERTST(11,2,2),PERTST(12,2,2)
     5           /1., 1., 1., 1./
      DATA PERTST(1,1,3),PERTST(2,1,3),PERTST(3,1,3),PERTST(4,1,3)
     1           /12., 24., 37.89, 53.33/,
     2     PERTST(5,1,3),PERTST(6,1,3),PERTST(7,1,3),PERTST(8,1,3)
     3           /70.08, 87.97, 106.9, 126.7/,
     4     PERTST(9,1,3),PERTST(10,1,3),PERTST(11,1,3),PERTST(12,1,3)
     5           /147.4, 168.8, 191., 1./
      DATA PERTST(1,2,3),PERTST(2,2,3),PERTST(3,2,3),PERTST(4,2,3)
     1           /3., 6., 9.167, 12.5/,
     2     PERTST(5,2,3),PERTST(6,2,3),PERTST(7,2,3),PERTST(8,2,3)
     3           /1., 1., 1., 1./,
     4     PERTST(9,2,3),PERTST(10,2,3),PERTST(11,2,3),PERTST(12,2,3)
     5           /1., 1., 1., 1./
C
      GO TO (1,2),METH
 1    MAXDER = 12
      GO TO (101,102,103,104,105,106,107,108,109,110,111,112),NQ
 2    MAXDER = 5
      GO TO (201,202,203,204,205),NQ
C-----------------------------------------------------------------------
C  THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO MACHINE ACCURACY.
C  FOR A GIVEN ORDER NQ, THEY CAN BE CALCULATED BY USE OF THE
C  GENERATING POLYNOMIAL L(T), WHOSE COEFFICIENTS ARE EL(I)..
C       L(T) = EL(1) + EL(2)*T + ... + EL(NQ+1)*T**NQ.
C  FOR THE IMPLICIT ADAMS METHODS, L(T) IS GIVEN BY
C       DL/DT = (T+1)*(T+2)* ... *(T+NQ-1)/K,    L(-1) = 0,
C  WHERE                 K = FACTORIAL(NQ-1).
C  FOR THE BDF METHODS,
C       L(T) = (T+1)*(T+2)* ... *(T+NQ)/K,
C  WHERE         K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ).
C
C  THE ORDER IN WHICH THE GROUPS APPEAR BELOW IS..
C  IMPLICIT ADAMS METHODS OF ORDERS 1 TO 12,
C  BDF METHODS OF ORDERS 1 TO 5.
C-----------------------------------------------------------------------
 101  EL(1) = 1.00
      GO TO 900
 102  EL(1) = 0.50
      EL(3) = EL(1)
      GO TO 900
 103  EL(1) = 4.1666666666666667E-01
      EL(3) = 0.750
      EL(4) = 1.6666666666666667E-01
      GO TO 900
 104  EL(1) = 0.3750
      EL(3) = 9.1666666666666667E-01
      EL(4) = 3.3333333333333333E-01
      EL(5) = 4.1666666666666667E-02
      GO TO 900
 105  EL(1) = 3.4861111111111111E-01
      EL(3) = 1.04166666666666670
      EL(4) = 4.8611111111111111E-01
      EL(5) = 1.0416666666666667E-01
      EL(6) = 8.3333333333333333E-03
      GO TO 900
 106  EL(1) = 3.2986111111111111E-01
      EL(3) = 1.1416666666666667E+00
      EL(4) = 0.625E+00
      EL(5) = 1.7708333333333333E-01
      EL(6) = 0.025E+00
      EL(7) = 1.3888888888888889E-03
      GO TO 900
 107  EL(1) = 3.1559193121693122E-01
      EL(3) = 1.225E+00
      EL(4) = 7.5185185185185185E-01
      EL(5) = 2.5520833333333333E-01
      EL(6) = 4.8611111111111111E-02
      EL(7) = 4.8611111111111111E-03
      EL(8) = 1.9841269841269841E-04
      GO TO 900
 108  EL(1) = 3.0422453703703704E-01
      EL(3) = 1.2964285714285714E+00
      EL(4) = 8.6851851851851852E-01
      EL(5) = 3.3576388888888889E-01
      EL(6) = 7.7777777777777778E-02
      EL(7) = 1.0648148148148148E-02
      EL(8) = 7.9365079365079365E-04
      EL(9) = 2.4801587301587302E-05
      GO TO 900
 109  EL(1) = 2.9486800044091711E-01
      EL(3) = 1.3589285714285714E+00
      EL(4) = 9.7655423280423280E-01
      EL(5) = 0.4171875E+00
      EL(6) = 1.1135416666666667E-01
      EL(7) = 0.01875E+00
      EL(8) = 1.9345238095238095E-03
      EL(9) = 1.1160714285714286E-04
      EL(10)= 2.7557319223985891E-06
      GO TO 900
 110  EL(1) = 2.8697544642857143E-01
      EL(3) = 1.4144841269841270E+00
      EL(4) = 1.0772156084656085E+00
      EL(5) = 4.9856701940035273E-01
      EL(6) = 0.1484375E+00
      EL(7) = 2.9060570987654321E-02
      EL(8) = 3.7202380952380952E-03
      EL(9) = 2.9968584656084656E-04
      EL(10)= 1.3778659611992945E-05
      EL(11)= 2.7557319223985891E-07
      GO TO 900
 111  EL(1) = 2.8018959644393672E-01
      EL(3) = 1.4644841269841270E+00
      EL(4) = 1.1715145502645503E+00
      EL(5) = 5.7935819003527337E-01
      EL(6) = 1.8832286155202822E-01
      EL(7) = 4.1430362654320988E-02
      EL(8) = 6.2111441798941799E-03
      EL(9) = 6.2520667989417989E-04
      EL(10)= 4.0417401528512640E-05
      EL(11)= 1.5156525573192240E-06
      EL(12)= 2.5052108385441719E-08
      GO TO 900
 112  EL(1) = 2.7426554003159906E-01
      EL(3) = 1.5099386724386724E+00
      EL(4) = 1.2602711640211640E+00
      EL(5) = 6.5923418209876543E-01
      EL(6) = 2.3045800264550265E-01
      EL(7) = 5.5697246105232216E-02
      EL(8) = 9.4394841269841270E-03
      EL(9) = 1.1192749669312169E-03
      EL(10)= 9.0939153439153439E-05
      EL(11)= 4.8225308641975309E-06
      EL(12)= 1.5031265031265031E-07
      EL(13)= 2.0876756987868099E-09
      GO TO 900
C
 201  EL(1) = 1.0E+00
      GO TO 900
 202  EL(1) = 6.6666666666666667E-01
      EL(3) = 3.3333333333333333E-01
      GO TO 900
 203  EL(1) = 5.4545454545454545E-01
      EL(3) = EL(1)
      EL(4) = 9.0909090909090909E-02
      GO TO 900
 204  EL(1) = 0.48E+00
      EL(3) = 0.7E+00
      EL(4) = 0.2E+00
      EL(5) = 0.02E+00
      GO TO 900
 205  EL(1) = 4.3795620437956204E-01
      EL(3) = 8.2116788321167883E-01
      EL(4) = 3.1021897810218978E-01
      EL(5) = 5.4744525547445255E-02
      EL(6) = 3.6496350364963504E-03
C
 900  DO 910 K = 1,3
 910    TQ(K) = PERTST(NQ,METH,K)
      TQ(4) = .5*TQ(2)/FLOAT(NQ+2)
      RETURN
C----------------------- END OF SUBROUTINE COSET -----------------------
      END
      SUBROUTINE DECBR (MDIM, N, ML, MU, B, IP, IER)                     DECBR
      INTEGER MDIM,N,ML,MU,IP,IER
      REAL  B
      DIMENSION B(MDIM,N), IP(N)
C-----------------------------------------------------------------------
C  THIS SUBROUTINE CONSTRUCTS THE LU DECOMPOSITION OF A BAND MATRIX A
C  IN THE FORM L*U = P*A , WHERE P IS A PERMUTATION MATRIX, L IS A
C  UNIT LOWER TRIANGULAR MATRIX, AND U IS AN UPPER TRIANGULAR MATRIX.
C  THE MATRIX A IS ASSUMED HERE TO BE STORED BY ROWS IN THE ARRAY B.
C  MDIM  =  THE FIRST DIMENSION (COLUMN LENGTH) OF THE ARRAY B.
C           MDIM MUST BE AT LEAST 2*ML+MU+1.
C  N     =  ORDER OF MATRIX.
C  ML,MU =  WIDTHS OF THE LOWER AND UPPER PARTS OF THE BAND, NOT
C           COUNTING THE MAIN DIAGONAL. TOTAL BANDWIDTH IS ML+MU+1.
C  B     =  MDIM BY N ARRAY CONTAINING THE MATRIX A ON INPUT
C           AND ITS FACTORED FORM ON OUTPUT.
C           ON INPUT, A IS STORED BY ROWS IN B AS AN (ML+MU+1) BY N
C           ARRAY (THE ROWS OF A ARE THE COLUMNS OF B), WITH A COLUMN
C           LENGTH OF MDIM FOR B.  THUS B(K,I) = B(K+(I-1)*MDIM)
C           (1 .LE. K .LE. ML+MU+1) CONTAINS THE PART OF THE I-TH ROW
C           OF A WITHIN THE BAND.  THAT IS, B(J-I+ML+1,I) CONTAINS
C           A(I,J) FOR -ML .LE. (J-I) .LE. MU.
C           ON OUTPUT, B CONTAINS THE L AND U FACTORS, WITH U IN ROWS
C           1 TO ML+MU+1, AND -L IN ROWS ML+MU+2 TO 2*ML+MU+1.
C           THE DIAGONAL OF L IS NOT STORED, AND THE RECIPROCALS
C           OF THE DIAGONAL ELEMENTS OF U ARE STORED.
C  IP    =  ARRAY OF LENGTH N CONTAINING PIVOT INFORMATION ON OUTPUT.
C  IER   =  ERROR INDICATOR..
C        =  0  IF NO ERROR,
C        = -1  IF MDIM, N, ML, AND/OR MU IS ILLEGAL,
C        =  K  IF THE K-TH PIVOT CHOSEN WAS ZERO (A IS SINGULAR).
C  THE INPUT ARGUMENTS ARE  MDIM, N, ML, MU, AND B.
C  THE OUTPUT ARGUMENTS ARE  B, IP, AND IER.
C-----------------------------------------------------------------------
      INTEGER I,II,J,K,KK,LL,LR,MX,NP,NR,N1
      REAL  XM,XX
      N1 = N - 1
      LL = ML + MU + 1
      IER = MIN0(N1,ML,MU,N1-ML,N1-MU,MDIM-LL-ML)
      IF (IER .LT. 0) GO TO 110
      IER = 0
      IF (N .EQ. 1) GO TO 92
      IF (ML .EQ. 0) GO TO 32
C  INITIALIZE STORAGE.  SHIFT AND ZERO OUT ELEMENTS OF B. --------------
      DO 30 I = 1,ML
        II = MU + I
        K = ML + 1 - I
        DO 10 J = 1,II
        JPK = J + K
 10     B(J,I) = B(JPK,I)
        K = II + 1
        DO 20 J = K,LL
 20       B(J,I) = 0.0
 30     CONTINUE
 32   LR = ML
      DO 90 NR = 1,N1
        NP = NR + 1
        IF (LR .NE. N) LR = LR + 1
C  FIND THE PIVOT IN COLUMN NR, WITHIN THE BAND, AT OR BELOW ROW NR. ---
        MX = NR
        XM = ABS(B(1,NR))
        IF (ML .EQ. 0) GO TO 42
        DO 40 I = NP,LR
          IF (ABS(B(1,I)) .LE. XM) GO TO 40
          MX = I
          XM = ABS(B(1,I))
 40       CONTINUE
 42     IP(NR) = MX
        IF (MX .EQ. NR) GO TO 60
C  INTERCHANGE ROWS NR AND MX. -----------------------------------------
        DO 50 I = 1,LL
          XX = B(I,NR)
          B(I,NR) = B(I,MX)
 50       B(I,MX) = XX
 60     XM = B(1,NR)
        IF (XM .EQ. 0.0) GO TO 100
C  CONSTRUCT ELEMENTS OF L AND U. --------------------------------------
        B(1,NR) = 1.0/XM
        IF (ML .EQ. 0) GO TO 90
        XM = -B(1,NR)
        KK = MIN0(N-NR,LL-1)
        DO 80 I = NP,LR
          J = LL + I - NR
          XX = B(1,I)*XM
          B(J,NR) = XX
          DO 70 II = 1,KK
 70         B(II,I) = B(II+1,I) + XX*B(II+1,NR)
 80       B(LL,I) = 0.0
 90     CONTINUE
 92   NR = N
      IF (B(1,N) .EQ. 0.0) GO TO 100
      B(1,N) = 1.0/B(1,N)
      RETURN
C
 100  IER = NR
      RETURN
C
 110  IER = -1
      RETURN
C-----------------------  END OF SUBROUTINE DECBR  ---------------------
      END
      SUBROUTINE SOLBR (MDIM, N, ML, MU, B, IP, Y)                       SOLBR
      INTEGER MDIM,N,ML,MU,IP
      REAL  B,Y
      DIMENSION B(MDIM,N), IP(N), Y(N)
C-----------------------------------------------------------------------
C  THIS SUBROUTINE COMPUTES THE SOLUTION OF THE BANDED LINEAR SYSTEM
C  A*X = C , GIVEN THE LU DECOMPOSITION OF THE MATRIX A FROM DECBR.
C  Y     =  RIGHT-HAND VECTOR C, OF LENGTH N, ON INPUT,
C        =  SOLUTION VECTOR X ON OUTPUT.
C  ALL OTHER ARGUMENTS ARE AS DESCRIBED IN DECBR.
C  ALL THE ARGUMENTS ARE INPUT ARGUMENTS.
C  THE OUTPUT ARGUMENT IS  Y.
C-----------------------------------------------------------------------
      INTEGER I,J,KK,LL,NB,NR,N1
      REAL  DP,XX
      IF (N .EQ. 1) GO TO 60
      N1 = N - 1
      LL = ML + MU + 1
C  IF ML = 0, SKIP THE FIRST BACKSOLUTION PHASE. -----------------------
      IF (ML .EQ. 0) GO TO 32
C  APPLY PERMUTATION AND L MULTIPLIERS TO Y. ---------------------------
      DO 30 NR = 1,N1
        IF (IP(NR) .EQ. NR) GO TO 10
        J = IP(NR)
        XX = Y(NR)
        Y(NR) = Y(J)
        Y(J) = XX
 10     KK = MIN0(N-NR,ML)
        DO 20 I = 1,KK
        NRPI = NR + I
        LLPI = LL + I
 20     Y(NRPI) = Y(NRPI) + Y(NR)*B(LLPI,NR)
 30     CONTINUE
 32   LL = LL - 1
C  BACKSOLVE FOR X USING U. --------------------------------------------
      Y(N) = Y(N)*B(1,N)
      KK = 0
      DO 50 NB = 1,N1
        NR = N - NB
        IF (KK .NE. LL) KK = KK + 1
        DP = 0.0
        IF (LL .EQ. 0) GO TO 50
        DO 40 I = 1,KK
        IP1 = I + 1
        NRPI = NR + I
 40     DP = DP + B(IP1,NR)*Y(NRPI)
 50     Y(NR) = (Y(NR) - DP)*B(1,NR)
      RETURN
 60   Y(1) = Y(1)*B(1,1)
      RETURN
C-----------------------  END OF SUBROUTINE SOLBR  ---------------------
      END
      SUBROUTINE PDB(NODE,T,U,PW,N0,ML,MU)                               PDB
      DIMENSION U(1),PW(1)
C-----------------------------------------------------------------------
C  A USER DEFINED ROUTINE WHICH DEFINES THE JACOBIAN IF MITER = 1.
C  SINCE THIS METHOD IS NOT RECOMMENDED, A DUMMY ROUTINE IS PROVIDED.
C  THE CALLING PARAMETERS ARE DEFINED IN THE COMMENTS OF SUBROUTINE
C  PSETM.  NOTE THAT NODE = NPDE*NX*NY, N0 = NODE, AND MU = ML.  THE
C  DETAILED COMMENTS ON HOW TO WRITE THIS ROUTINE ARE TO BE FOUND IN
C  REFERENCE 1 (SEE PSETM).
C-----------------------------------------------------------------------
      RETURN
      END

		
.

NEW PAGES:

[ODDNUGGET]

[GOPHER]