      SUBROUTINE PDEDIF(FUNC,BNBDY,NPDE,NZ,TEND,EPS,ZPTS,ITYPE,MBLFT,
     1 MBRHT,NBUG,T,DT,INDEX,USOLN,WKCM,IWCM,NWKCM,NIWCM)
C     .
C     VERSION OF PDEDIF PACKAGE DATED   02/27/80
C     SOLVES A PARABOLIC SYSTEM OF EQUATIONS IN ONE SPACE DIMENSION
C     AND TIME.  TWO TYPES OF EQUATIONS ARE ALLOWED, FOR ITYPE=1
C     .
C     DU/DT = D(A(U)*DU/DZ)/DZ + D(B(U)*U)/DZ + C(U,UZ)
C     .
C     AND FOR ITYPE=2
C     .
C     DU/DT = A(U,UZ)*D2U/D2Z + B(U,UZ)*DU/DZ +C(U,UZ)
C     .
C     THE COEFFICIENTS A AND B ARE VECTOR RATHER THAN MATRIX, THAT IS
C     THERE IS NO COUPLING OF THE EQUATIONS THRU THE Z DERIVATIVES.
C     ACTUALLY THE COEFFICIENTS DEPEND ON Z AND T AS WELL AS U,
C     THAT IS  A(M)=A(M,Z,T,U1,...UN)  WHERE M IS THE EQUATION NUMBER
C     THE COEFFICIENTS A, B, AND C ARE COMPUTED BY THE USER SUPPLIED
C     ROUTINE FUNC.  THE BOUNDARY CONDITIONS CAN BE EITHER
C     DIRICHLET OR MIXED, THE DIRICHLET MUST BE LINEAR IN U.  THE
C     TYPE OF THE BOUNDARY CONDITION CAN BE DIFFERENT FOR EACH
C     COMPONENT OF THE SOLUTION VECTOR.
C     THE BOUNDARY CONDITIONS HAVE THE FORM
C     U = BL(T) AT THE LEFT OR LOWER BOUNDARY AND
C     U = BR(T) AT THE RIGHT OR UPPER BOUNDARY
C     THE FOLLOWING FORM IS ALSO ALLOWED
C     DU/DZ = BL(M,UL,T)  OR  DU/DZ = BR(M,UR,T)
C     A SECOND ORDER FINITE DIFFERENCE SCHEME IS USED IN SPACE
C     AND THE METHOD OF LINES IS USED FOR TIME,  THE GEARB
C     INTEGRATOR OF HINDMARSH IS USED FOR TIME INTEGRATION.
C
C     THE FUNC SUBROUTINE COMPUTES THE COEFFICIENTS A,B, AND C AND
C     RETURNS THE VALUES IN THE ARRAYS AN,BN, AND CN .
C     .
C     SUBROUTINE FUNC(NPDE,NZP2,ZN,T,UN,UZN,AN,BN,CN)
C     .
C     HERE AN,BN, AND CN ARE OUTPUT PARAMETERS AND THE REMAINING
C     ARGUMENTS ARE INPUT PARAMETERS.
C     THE DIMENSIONS OF THE ARGUMENT ARRAYS ARE
C     REAL ZN(NZP2),UN(NPDE,NZP2),UZN(NPDE,NZP2),AN(NPDE,NZP2),
C    1 BN(NPDE,NZP2),CN(NPDE,NZP2)
C     WHERE NZP2=NZ+2.
C     NOTE THAT THE ORIGINAL MESH ARRAY ZPTS IS SHIFTED AND TWO
C     FICTITIOUS POINTS AT I=1 AND I=NZ+2 ARE ADDED TO FORM THE MESH
C     ARRAY ZN WHICH IS USED INTERNALLY IN THE PDEDIF PACKAGE.
C     BEWARE... THE ZPTS AND ZN ARRAYS ARE NOT THE SAME.
C     THE ARRAYS AN AND BN CAN DEPEND ON THE SOLUTION UN(M,I)
C     THE SOLUTION UN IS EXTRAPOLATED TO THE FICTITIOUS POINTS.
C     THE VALUES OF AN AND BN MUST BE COMPUTED WITHIN THIS ROUTINE AT
C     THESE FICTITIOUS POINTS.  THEREFORE AN AND BN CAN NOT DEPEND ON
C     THE DERIVATIVE UZN OF UN.  HOWEVER CN NEED NOT BE COMPUTED AT
C     THE FICTITIOUS POINTS AND THEREFORE CAN DEPEND ON UZN.  NOTE
C     THAT UZN IS AN INPUT PARAMTETER OF FUNC.  UZN IS NOT
C     DEFINED AT THE FICTITIOUS POINTS I=1 AND I=NZP2=NZ+2.
C     IF THE BOUNDARY TYPE IS MBLFT=1 OR MBLFT=3, THEN ALTHOUGH THE
C     VALUE OF AN AND BN AT THE FICTITIOUS POINTS MUST BE SET
C     THEY HAVE NO EFFECT ON THE SOLUTION.
C     .
C     THE SECOND USER SUPPLIED ROUTINE BNBDY IS USED TO DEFINE THE
C     BOUNDARY CONDITIONS.  THIS ROUTINE HAS THE ARGUMENTS
C     .
C     SUBROUTINE BNBDY(NPDE,T,ZL,ZR,UL,UR,BL,BR)
C     .
C     WHERE THE ARGUMENT ARRAYS HAVE THE DIMENSIONS
C     .
C     REAL UL(NPDE),UR(NPDE),BL(NPDE),BR(NPDE)
C     .
C     .
C     HERE T IS THE TIME.  THE LOWER OR LEFT BOUNDARY IS LOCATED AT
C     Z=ZL AND THE RIGHT OR UPPER BOUNDARY IS AT Z=ZR.  ZL AND ZR
C     ARE INPUT PARAMETERS.  THE BOUNDARY CONDITIONS HAVE THE FORM
C     .
C     U(M,ZL,T)=BL(M,T)     IF MBLFT(M)=1
C     OR
C     DU/DZ = BL(M,UL,T)      IF MBLFT(M)=2
C     OR
C     DU/DT(M,ZL,T) = BL(M,UL,T)      IF MBLFT(M)=3
C     .
C      THE UPPER OR RIGHT BOUNDARY IS SIMILAR.
C     .
C     U(M,ZR,T) = BR(M,T)
C     OR
C     DU/DZ(M,ZR,T) = BR(M,UR,T)
C     OR
C     DU/DT(M,ZR,T) = BR(M,UR,T)
C     .
C     .
C     .
C     THE INPUT/OUTPUT PARAMETERS ARE THE FOLLOWING
C     .
C     FUNC    - THE USER SUPPLIED SUBROUTINE TO DEFINE THE RIGHT SIDE.
C               SEE ABOVE FOR MORE DESCRIPTION
C     BNBDY   - THE USER SUPPLIED SUBROUTINE TO SUPPLY THE BOUNDARY
C               CONDITIONS.  SEE ABOVE.
C     NPDE    - THE NUMBER OF EQUATIONS OR SOLUTION VECTOR COMPONENTS
C     NZ      - THE NUMBER OF MESH POINTS IN Z
C     TEND    - THE STOPPING TIME FOR THE INTEGRATION
C     EPS     - THE ERROR TOLERANCE FOR THE GEARB INTEGRATOR
C     ZPTS    - THE ARRAY OF MESH POINTS  ZPTS(NZ)
C     ITYPE   - SELECTS THE FORMAT OF THE SYSTEM OF EQUATIONS
C             =1 FOR CONSERVATION FORM D(A*DU/DZ)/DZ
C             =2 FOR NON CONSERVATION FORM   A*D2U/D2Z + B*DU/DZ + C
C     MBLFT   - ARRAY GIVING TYPE OF BOUNDARY CONDITION FOR M-TH
C     SOLUTION COMPONENT AT THE LOWER BOUNDARY.
C                MBLFT(M) = 1  FOR   U(M,ZL,T) = BL(M,T)
C                MBLFT(M) = 2  FOR   DU/DZ(M,ZL,T) = BL(M,UL,T)
C                MBLFT(M) = 3  FOR   DU/DT(M,ZL,T) = BL(M,UL,T)
C     MBRHT   - ARRAY GIVING TYPE OF BOUNDARY CONDITION AT UPPER
C               OR RIGHT BOUNDARY.
C     NBUG    - SELECTS DEBUGGING PRINTING
C             =1 TO PRINT MESSAGE WHEN GEARB REJECTS A TIME STEP BECAUSE
C                OF A LARGE ERROR ESTIMATE.
C             =2 TO PRINT CERTAIN GEARB PARAMETERS ON EACH TIME
C                STEP
C             =4 TO PRINT DIAGNOSTIC INFORMATION IN THE DIFFN ROUTINE
C             =8 TO PRINT DIAGNOSTIC INFOMMATION IN THE DIFFUN ROUTINE
C     T       - ON INPUT THIS IS THE INITIAL TIME FOR THE INTEGRATION.
C               ON OUTPUT THIS IS THE TIME REACHED, WHICH IS TEND FOR
C               A SUCCESSFUL INTEGRATION.
C     DT       - ON INPUT THIS GIVES THE INITIAL GUESS FOR THE TIME STEP
C                IF IN DOUBT, UNDERESTIMATE THIS VALUE.  ON OUTPUT
C                THIS GIVES THE VALUE OF THE FINAL TIME STEP USED.
C     INDEX   - DETERMINES HOW THE TIME INTEGRATION WILL BE STARTED
C               AND ALSO IS USED AS THE ERROR FLAG ON OUTPUT.
C             =1 TO START A NEW INTEGRATION
C             =0 TO CONTINUE AN INTEGRATION
C             =0 ON OUTPUT IF THE INTEGRATION WAS SUCCESSFUL.
C     USOLN   - AN ARRAY HOLDING THE INITIAL VALUE OF THE SOLUTION ON
C               INPUT AND THE FINAL VALUE ON OUTPUT.  THE DIMENSION IS
C               USOLN(NPDE,NZ)
C     WKCM    - REAL WORK ARRAY.
C     NWKCM - DIMENSION OF WKCM.  MUST EXCEED (24+6*NPDE)*NPDE*NZP2
C              +5*NZP2 + 4*NPDE   WHERE NZP2=NZ+2
C     IWCM  - INTEGER WORK ARRAY
C     NIWCM    DIMENSION OF IWCM.  MUST EXCEED 2*NPDE
C     .
C     .
      REAL ZPTS(NZ),USOLN(NPDE,NZ), WKCM(NWKCM)
      INTEGER IWCM(NIWCM)
      INTEGER MBLFT(NPDE),MBRHT(NPDE)
      COMMON /PDVAR/ NPDEC,NZC,NBUGC,ITYPEC
      COMMON /SYSPAR/ UROUND,MIN,MOUT
      EXTERNAL FUNC,BNBDY
C     .
C     CHECK THE INPUT DATA
C     .
      CALL CKINPT(ZPTS,NZ,NPDE,T,TEND,DT,MBLFT,MBRHT,NWKCM,NIWCM,
     1 INDEX)
      IF(INDEX .LT. -31)RETURN
C     .
C     SET VARIABLES IN COMMON
C     .
      NPDEC=NPDE
      NZC=NZ
      ITYPEC=ITYPE
      NBUGC=NBUG
C     .
C     MOVE INPUT VARIABLES TO WORK ARRAY FOR COMMUNICATION WITH
C     OTHER ROUTINES
C     .
      NZP2=NZ+2
      LML=1
      LMR=LML+NPDE
      LZN=1
      LUN=LZN+NZP2
      LU=LUN+NPDE*NZP2
      LAL=5*NZP2+7*NPDE*NZP2+1
      LAR=LAL+NPDE
      LBL=LAR+NPDE
      LBR=LBL+NPDE
      LY=LBR+NPDE
      LYMAX=LY+13*NPDE*NZP2
      LERROR=LYMAX+NPDE*NZP2
      LSAV1=LERROR+NPDE*NZP2
      LSAV2=LSAV1+NPDE*NZP2
      LPW=LSAV2+NPDE*NZP2
      N0PW=(6*NPDE-2)*NPDE*NZP2
      LIPIV=LPW+N0PW
      DO 40 M=1,NPDE
      IWCM(M+LML-1)=MBLFT(M)
40    IWCM(M+LMR-1)=MBRHT(M)
C     .
C     ADD FICTITIOUS POINTS TO THE MESH
C     .
      DO 50 I=1,NZ
      WKCM(LZN+I)=ZPTS(I)
      DO 45 M=1,NPDE
      IM=LUN+M-1+I*NPDE
45    WKCM(IM)=USOLN(M,I)
50    CONTINUE
      WKCM(LZN)=ZPTS(1)-(ZPTS(2)-ZPTS(1))
      WKCM(LZN+NZP2-1)=ZPTS(NZ)+(ZPTS(NZ)-ZPTS(NZ-1))
C     .
C     COUNT THE NUMBER OF UNKNOWNS IN THE DIFFERENTIAL SYSTEM, THAT
C     IS SUBSTRACT THE DIRICHLET BOUNDARY VALUES AND THE FICTITIOUS
C     POINTS
C     .
      NK=0
      DO 60 M=1,NPDE
      IF(MBLFT(M) .NE. 1)NK=NK+1
      IF(MBRHT(M) .NE. 1)NK=NK+1
60    CONTINUE
      NEQ=NK + NPDE*(NZ-2)
C     .
C     MOVE INITIAL VALUES TO THE ARRAY U FOR INTEGRATION BY GEARB
C     NOTE THAT DIRICHLET BOUNDARY VALUES ARE DROPPED FROM U.
C     .
      CALL MVTOU(NPDE,NZ,NEQ,USOLN,MBLFT,MBRHT,WKCM(LU))
C     .
C     CALL THE GEARB ODE SOLVER
C     .
      MF=22
      ML=2*NPDE-1
      MU=ML
      TBGN=T
      TLAST=TEND

      CALL DRIVEB(NEQ,TBGN,DT,WKCM(LU),TLAST,EPS,MF,INDEX,ML,MU,
     1 N0PW,WKCM(LY),WKCM(LYMAX),WKCM(LERROR),WKCM(LSAV1),
     2 WKCM(LSAV2),WKCM(LPW),WKCM(LIPIV),WKCM,IWCM,NWKCM,NIWCM,NBUG,
     1 FUNC,BNBDY)
C     .
C     MOVE THE FINAL SOLUTION TO THE ARRAY USOLN
C     .
      CALL MVTOUN(TLAST,NPDE,NZP2,NEQ,WKCM(LU),MBLFT,MBRHT,
     1 WKCM(LBL),WKCM(LBR),WKCM(LUN),WKCM(LZN),BNBDY)
      CALL MVUNUS(NPDE,NZP2,NZ,WKCM(LUN),USOLN)
      T=TLAST
      RETURN
      END
      SUBROUTINE DIFFUN(NEQ,T,U,UT,WKCM,IWCM,NWKCM,NIWCM,
     1 FUNC,BNBDY)
C     .
C     THE ROUTINE REQUIRED BY GEARB TO EVALUTE THE TIME DERIVATIVE
C     .
      REAL U(NEQ),UT(NEQ),WKCM(NWKCM)
      INTEGER IWCM(NIWCM)
      COMMON /PDVAR/ NPDE,NZ,NBUG,ITYPE
      COMMON /SYSPAR/ UROUND,MIN,MOUT
      EXTERNAL FUNC,BNBDY
C     .
      NZP2=NZ+2
      LML=1
      LMR=LML+NPDE
      LZN=1
      LUN=LZN+NZP2
      LU=LUN+NPDE*NZP2
      LUT=LU+NPDE*NZP2
      LUZN=LUT+NPDE*NZP2
      LUT1=LUZN+NPDE*NZP2
      LDM=LUT1+NZP2
      LAN=LDM+3*NZP2
      LBN=LAN+NPDE*NZP2
      LCN=LBN+NPDE*NZP2
      LAL=LCN+NPDE*NZP2
      LAR=LAL+NPDE
      LBL=LAR+NPDE
      LBR=LBL+NPDE

      CALL DIFFN(NPDE,NZP2,NEQ,T,WKCM(LZN),U,UT,WKCM(LUT1),WKCM(LDM),
     1 WKCM(LUN),WKCM(LAN),WKCM(LBN),WKCM(LCN),WKCM(LAL),WKCM(LAR),
     2 WKCM(LBL),WKCM(LBR),IWCM(LML),IWCM(LMR),WKCM(LUZN),FUNC,BNBDY)
C     .
C     DEBUG PRINTOUT WHEN NBUG=4
C     .
      IF(NBUG/4-2*(NBUG/8) .NE. 1)GOTO 75
      WRITE (MOUT,65) T
65    FORMAT(//3X,25HNBUG=1 PRINTING IN DIFFUN,6H AT T=,
     1 E12.5/5X,1HK,8X,4HU(I),7X,5HUT(I))
      DO 70 K=1,NEQ
      WRITE (MOUT,67) K,U(K),UT(K)
67    FORMAT(3X,I3,3X,2E13.5)
70    CONTINUE
75    CONTINUE
      RETURN
      END
      SUBROUTINE DIFFN(NPDE,NZP2,NEQ,T,ZN,U,UT,UT1,DM,UN,AN,BN,CN,AL,
     1 AR,BL,BR,MBLFT,MBRHT,UZN,FUNC,BNBDY)
C     .
C     THIS ROUTINE COMPUTES THE TIME DERIVATIVE GIVEN THE SPATIAL
C     DERIVATIVES.
C     .
      REAL ZN(NZP2),UN(NPDE,NZP2),U(NEQ),UT(NEQ),UT1(NZP2),
     1 BN(NPDE,NZP2),CN(NPDE,NZP2),AL(NPDE),AR(NPDE),DM(NZP2,3),
     2 AN(NPDE,NZP2),BL(NPDE),BR(NPDE),UZN(NPDE,NZP2)
      INTEGER MBLFT(NPDE),MBRHT(NPDE)
      COMMON /PDVAR/ NPDEC,NZC,NBUG,ITYPE
      COMMON /SYSPAR/ UROUND,MIN,MOUT
      EXTERNAL FUNC,BNBDY
C     .
      NZ=NZP2-2
      NZP1=NZP2-1

C     .
C     MOVE THE SOLUTION IN THE U ARRAY TO THE ARRAY UN AND SET THE
C     BOUNDARY VALUES FROM THE BOUNDARY CONDITIONS
C     .
      CALL MVTOUN(T,NPDE,NZP2,NEQ,U,MBLFT,MBRHT,BL,BR,UN,ZN,BNBDY)
C     .
C     COMPUTE THE TIME DERIVATIVE AND STORE TEMPORARILY IN THE UT1 ARRAY
C     .
C     .
C     DETERMINE THE NUMBER OF MIXED BOUNDARY VARIABLES AT THE LEFT
C     (NQLFT) AND AT THE RIGHT (NQRHT)
C     .
      NQLFT=0
      NQRHT=0
      DO 50 M=1,NPDE
      IF(MBLFT(M) .NE. 1)NQLFT=NQLFT+1
      IF(MBRHT(M) .NE. 1)NQRHT=NQRHT+1
50    CONTINUE
C     .
C     COMPUTE THE DERIVATIVE UZN BY 2-ND ORDER DIFFERENCES
C     .
      DO 110 I=2,NZP1
      DM(I,1)=1./((ZN(I)-ZN(I-1))*(ZN(I+1)-ZN(I-1)))
      DM(I,2)=1./((ZN(I)-ZN(I-1))*(ZN(I+1)-ZN(I)))
      DM(I,3)=1./((ZN(I+1)-ZN(I))*(ZN(I+1)-ZN(I-1)))
110   CONTINUE
      DO 120 M=1,NPDE
      DO 120 I=2,NZP1
      UZN(M,I)= -(ZN(I+1)-ZN(I))*UN(M,I-1)*DM(I,1)
     1    +(ZN(I+1)-2.*ZN(I)+ZN(I-1))*UN(M,I)*DM(I,2)
     2    +(ZN(I)-ZN(I-1))*UN(M,I+1)*DM(I,3)
120   CONTINUE
C
C     COMPUTE THE TIME DERIVATVE AT THE BOUNDARY IF THERE IS A TYPE
C     THREE BOUNDARY CONDITION (MBLFT=3 OR MBRHT=3)
C
      DO 130 M=1,NPDE
      IF(MBLFT(M) .EQ. 3)GOTO 140
      IF(MBRHT(M) .EQ. 3)GOTO 140
130   CONTINUE
      GOTO 150
140   CALL BNBDY(NPDE,T,ZN(2),ZN(NZP1),UN(1,2),UN(1,NZP1),
     1 BL,BR)
150   CONTINUE
C     .
C     EVALUATE COEFFICIENT FUNCTIONS
C     .
cjpe      write(6,"(' CN=',/,(6e12.4))") CN(m,:)
      CALL FUNC(NPDE,NZP2,ZN,T,UN,UZN,AN,BN,CN)
cjpe      write(6,"(' CN=',/,(6e12.4))") CN(m,:)
C     .
C     DEBUG PRINTOUT NBUG=8
C     .
      IF(NBUG/8-2*(NBUG/16) .NE. 1)GOTO 190
      WRITE (MOUT,180) T
180   FORMAT(//3X,24HNBUG=2 PRINTING IN DIFFN,6H AT T=,
     1 E12.5)
      DO 185 M=1,NPDE
      WRITE (MOUT,181) M
181   FORMAT(/3X,2HM=,I3/5X,1HI,9X,2HZN,9X,2HUN,8X,3HUZN,
     1 9X,2HAN,9X,2HBN,9X,2HCN)
      DO 185 I=1,NZP2
      WRITE (MOUT,182) I,ZN(I),UN(M,I),UZN(M,I),AN(M,I),BN(M,I),
     1 CN(M,I)
182   FORMAT(1X,I3,2X,6E12.4)
185   CONTINUE
190   CONTINUE
C     .
      KTL=0
      KTR=NPDE*(NZ-2) + NQLFT
      DO 280 M=1,NPDE
      IF(ITYPE .EQ. 2)GOTO 222
C     .
C     CONSERVATION FORM OF THE SYSTEM
C     .
      DO 220 I=2,NZP1
      T1=.5*(AN(M,I+1)+AN(M,I))*(UN(M,I+1)-UN(M,I))
      T2=T1/(ZN(I+1)-ZN(I))
      T3=.5*(AN(M,I)+AN(M,I-1))*(UN(M,I)-UN(M,I-1))
      T4=T3/(ZN(I)-ZN(I-1))
      T5=.25*(BN(M,I+1)+BN(M,I))*(UN(M,I+1)+UN(M,I))
      T6=.25*(BN(M,I)+BN(M,I-1))*(UN(M,I)+UN(M,I-1))
      T7=(T2-T4+T5-T6)/(.5*(ZN(I+1)-ZN(I-1)))
      T8=T7+CN(M,I)
      UT1(I)=T8
220   CONTINUE
      GOTO 228
C     .
C     NON CONSERVATION FORM OF THE EQUATIONS
C     .
222   DO 224 I=2,NZP1
      UT1(I)=AN(M,I)*2.*(UN(M,I-1)*DM(I,1) - UN(M,I)*DM(I,2)
     1                 +UN(M,I+1)*DM(I,3))
     2 +BN(M,I)*(-(ZN(I+1)-ZN(I))*UN(M,I-1)*DM(I,1)
     3         +(ZN(I+1)-2.*ZN(I)+ZN(I-1))*UN(M,I)*DM(I,2)
     4         +(ZN(I)-ZN(I-1))*UN(M,I+1)*DM(I,3) )
     5   + CN(M,I)
224   CONTINUE
cc      write(6,"(' UN=',/,(6e12.4))") UN(m,:)
cc      write(6,"(' BN=',/,(6e12.4))") BN(m,:)

228   CONTINUE
C
C     CORRECT TIME DERIVATIVE AT THE BOUNDARY
C
      IF(MBLFT(M) .EQ. 3)UT1(2)=BL(M)
      IF(MBRHT(M) .EQ. 3)UT1(NZP1)=BR(M)
C     .
C     STORE THE TIME DERIVATIVE IN THE ARRAY UT FOR USE BY GEARB
C     .
      IF(MBLFT(M) .EQ. 1)GOTO 230
      KTL=KTL+1
      UT(KTL)=UT1(2)
230   IF(MBRHT(M) .EQ. 1)GOTO 240
      KTR=KTR+1
      UT(KTR)=UT1(NZP1)
240   KT=NQLFT+M
      DO 245 I=3,NZ
      UT(KT)=UT1(I)
245   KT=KT+NPDE
C     .
C     DEBUG PRINTOUT NBUG=8
C     .
      IF(NBUG/8-2*(NBUG/16) .NE. 1)GOTO 280
      WRITE (MOUT,270) T,M
270   FORMAT(/3X,24HNBUG=4 PRINTING IN DIFFN,6H AT T=,
     1 E12.5,3X,2HM=,I3/10X,12HUT1(2..NZP1))
      WRITE (MOUT,272) (UT1(I),I=2,NZP1)
272   FORMAT(1X,6E12.4)
280   CONTINUE
C     .
      RETURN
      END
      SUBROUTINE MVTOUN(T,NPDE,NZP2,NEQ,U,MBLFT,MBRHT,
     1 BL,BR,UN,ZN,BNBDY)
C     .
C     MOVE THE SOLUTION FROM THE U ARRAY TO THE UN ARRAY AND ADD
C     THE BOUNDARY VALUES USING THE BOUNDARY VALUES FROM BNBDY
C     .
      REAL U(NEQ),UN(NPDE,NZP2),BL(NPDE),BR(NPDE)
     1 ,ZN(NZP2)
      INTEGER MBLFT(NPDE),MBRHT(NPDE)
C     .
      NZ=NZP2-2
      NZP1=NZP2-1
      KT=0
C     .
      DO 20 M=1,NPDE
      IF(MBLFT(M) .EQ. 1)GOTO 20
      KT=KT+1
      UN(M,2)=U(KT)
20    CONTINUE
      NQLFT=KT
      DO 40 M=1,NPDE
      KT=NQLFT+M
      DO 35 I=3,NZ
      UN(M,I)=U(KT)
35    KT=KT+NPDE
40    CONTINUE
      KT=NQLFT+NPDE*(NZ-2)
      DO 50 M=1,NPDE
      IF(MBRHT(M) .EQ. 1)GOTO 50
      KT=KT+1
      UN(M,NZP2-1)=U(KT)
50    CONTINUE
C     .
C     SET DIRICHLET POINTS AND ALSO SET FICTITIOUS POINTS
C     .

      CALL BNBDY(NPDE,T,ZN(2),ZN(NZP2-1),UN(1,2),UN(1,NZP1),BL,BR)
      DO 125 M=1,NPDE
C     .
C     LOWER BOUNARY
C     .
      IF(MBLFT(M) .GT. 1)GOTO 110
C     .
C     DIRICHLET BOUNDARY
C     .
      UN(M,2)=BL(M)
      UN(M,1)=UN(M,2)+((ZN(1)-ZN(2))/(ZN(3)-ZN(2)))*(UN(M,3)-UN(M,2))
      GOTO 115
C     .
C     NEUMANN BOUNDARY
C     .
110   UN(M,1)=UN(M,3)-(ZN(3)-ZN(1))*BL(M)
C     .
C     UPPER BOUNDARY
C     .
115   IF(MBRHT(M) .GT. 1)GOTO 120
C     .
C     DIRICHLET BOUNDARY
C     .
      UN(M,NZP2-1)=BR(M)
      UN(M,NZP2)=UN(M,NZP2-1)+((ZN(NZP2)-ZN(NZP1))/(ZN(NZP1)-ZN(NZ)))*
     1                         (UN(M,NZP1)-UN(M,NZ))
      GOTO 125
C     .
C     NEUMANN BOUNDARY
C     .
120   UN(M,NZP2)=UN(M,NZ)+(ZN(NZP2)-ZN(NZ))*BR(M)
125   CONTINUE
      RETURN
      END
      SUBROUTINE MVTOU(NPDE,NZ,NEQ,USOLN,MBLFT,MBRHT,U)
C     .
C     MOVE SOLUTION FROM ARRAY USOLN TO ARRAY U.  DROPS DIRICHLET
C     BOUNDARY VALUES
C     .
      REAL U(NEQ),USOLN(NPDE,NZ)
      INTEGER MBLFT(NPDE),MBRHT(NPDE)
C     .
      NZM1=NZ-1
C     .
      KT=0
      DO 20 M=1,NPDE
      IF(MBLFT(M) .EQ. 1)GOTO 20
      KT=KT+1
      U(KT)=USOLN(M,1)
20    CONTINUE
      NQLFT=KT
      DO 40 M=1,NPDE
      KT=NQLFT+M
      DO 35 I=2,NZM1
      U(KT)=USOLN(M,I)
35    KT=KT+NPDE
40    CONTINUE
      KT=NQLFT+NPDE*(NZ-2)
      DO 50 M=1,NPDE
      IF(MBRHT(M) .EQ. 1)GOTO 50
      KT=KT+1
      U(KT)=USOLN(M,NZ)
50    CONTINUE
      RETURN
      END
      SUBROUTINE MVUNUS(NPDE,NZP2,NZ,UN,USOLN)
C     .
C     MOVE SOLUTION IN UN ARRAY TO USOLN ARRAY
C     .
      REAL UN(NPDE,NZP2),USOLN(NPDE,NZ)
C     .
      DO 20 M=1,NPDE
      DO 10 I=1,NZ
10    USOLN(M,I)=UN(M,I+1)
20    CONTINUE
      RETURN
      END
      SUBROUTINE CKINPT(ZPTS,NZ,NPDE,TBGN,TEND,DT,MBLFT,MBRHT,
     1 NWKCM,NIWCM,INDEX)
C     .
C     CHECK THE INPUT PARAMETERS FOR ERRORS
C     .
      REAL ZPTS(NZ)
      INTEGER MBLFT(NPDE),MBRHT(NPDE)
      COMMON /SYSPAR/ UROUND,MIN,MOUT
      DATA NZSAV,NPDSAV,ZPSAV/-1,-1,-1./
C     .
      IDX=0
      IF(INDEX .GE. -1 .AND. INDEX .LE. 3)GOTO 40
      WRITE (MOUT,35) INDEX
35    FORMAT(5H ****,7H INDEX=,I6,16H IS OUT OF RANGE)
      IDX=IDX-32
      GOTO 230
C     .
40    IF(NPDE .GE. 1)GOTO 50
      IDX=IDX-32
      WRITE (MOUT,45)
45    FORMAT(//5H ****,12H NPDE .LT. 1)
      GOTO 230
C     .
50    IF(NZ .GE. 3)GOTO 60
      IDX=IDX-64
      WRITE (MOUT,55)
55    FORMAT(//5H ****,10H NZ .LT. 3)
      GOTO 230
C     .
C     .
C     CHECK THAT PARAMETERS WERE NOT CHANGED ON CONTINUATION OF
C     INTEGRATION
C     .
60    IF(INDEX .EQ. 1)GOTO 138
      IF(NZ .NE. NZSAV)GOTO 134
      IF(NPDE .NE. NPDSAV)GOTO 134
      NZD2=NZ/2 + 1
      ZP=ZPTS(1)+ZPTS(NZD2)+ZPTS(NZ)
      IF(ABS(ZP-ZPSAV) .LE. 100.*UROUND)GOTO 140
134   IDX=IDX-128
      WRITE (MOUT,135)
135   FORMAT(5H ****,17H NZ NPDE OR ZPTS ,
     1 38HCHANGED ON CONTINUATION OF INTEGRATION /
     2 5X,31HSET INDEX=1 FOR NEW INTEGRATION )
      GOTO 230
138   NZSAV=NZ
      NPDSAV=NPDE
      NZD2=NZ/2 + 1
      ZPSAV=ZPTS(1)+ZPTS(NZD2)+ZPTS(NZ)
C     .
C     CHECK FOR MONOTONE ZPTS
C     .
140   I1=0
      NZP2=NZ+2
      DO 145 K=2,NZ
      IF(ZPTS(K) .LE. ZPTS(K-1))I1=1
145   CONTINUE
      IF(I1 .EQ. 0)GOTO 160
      WRITE (MOUT,150)
150   FORMAT(//5H ****,30H ZPTS(I) NOT INCREASING WITH I)
      IDX=IDX-256
C     .
160   IF(TBGN .LT. TEND)GOTO 170
      IDX=IDX-512
      WRITE (MOUT,165)
165   FORMAT(//5H ****,15H TBGN .GE. TEND)
C     .
170   IF(DT .GT. 0.)GOTO 180
      IDX=IDX-1024
      WRITE (MOUT,175)
175   FORMAT(//5H ****,11H DT .LE. 0.)
C     .
180   I1=0
      DO 195 M=1,NPDE
      IF(MBLFT(M) .LT. 1 .OR. MBLFT(M) .GT.  3)I1=1
      IF(MBRHT(M) .LT. 1 .OR. MBRHT(M) .GT. 3)I1=1
195   CONTINUE
      IF(I1 .EQ. 0)GOTO 210
      IDX=IDX-2048
      WRITE (MOUT,205)
205   FORMAT(//5H ****,28H MBLFT OR MBRHT OUT OF RANGE)
C     .
210   NLIM=(24+6*NPDE)*NPDE*NZP2 + 5*NZP2 + 4*NPDE
      IF(NWKCM .GE. NLIM)GOTO 220
      IDX=IDX-4096
      WRITE (MOUT,215) NLIM
215   FORMAT(//5H ****,18H NWKCM MUST EXCEED,I6)
C     .
220   IF(NIWCM .GE. 2*NPDE)GOTO 230
      WRITE (MOUT,225)
225   FORMAT(5H ****,23HNIWCM MUST BE GE 2*NPDE)
      IDX=IDX-8192
C     .
230   CONTINUE
      IF(IDX .LT. 0)INDEX=IDX
      RETURN
      END
      SUBROUTINE PDB
C     .
C     DUMMY ROUTINE FOR COMPUTATION OF JACOBIAN. SHOULD NOT BE USED
C     .
      WRITE (MOUT,10)
10    FORMAT(5H ****,11HSTOP IN PDB)
      STOP
      END
      SUBROUTINE DRIVEB (N, T0, H0, Y0, TOUT, EPS, MF, INDEX, ML, MU,
     1 N0PW,Y,YMAX,ERROR,SAVE1,SAVE2,PW,IPIV,WKCM,IWCM,NWKCM,
     2 NIWCM,NBUG,FUNC,BNBDY)
C        THIS IS MODIFIED FROM THE MARCH 19, 1975 VERSION OF
C        GEARB, A PACKAGE FOR THE SOLUTION OF THE INITIAL VALUE
C        PROBLEM FOR SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS,
C                  DY/DT = F(Y,T),    Y = (Y(1),Y(2),...,Y(N)).
C        GEARB IS A VARIANT OF THE GEAR PACKAGE TO BE USED WHEN
C        THE JACOBIAN MATRIX DF/DY HAS BANDED OR NEARLY BANDED FORM.
C        SUBROUTINE DRIVEB IS A DRIVER ROUTINE FOR THE GEARB PACKAGE.
C     .
C                       REFERENCES
C       1.  A. C. HINDMARSH, GEAR.. ORDINARY DIFFERENTIAL EQUATION
C             SYSTEM SOLVER, UCID-30001 REV. 3, LAWRENCE LIVERMORE
C             LABORATORY, P.O.BOX 808, LIVERMORE, CA 94550, DEC. 1974.
C     .
C       2.  A. C. HINDMARSH, GEARB.. SOLUTION OF ORDINARY
C             DIFFERENTIAL EQUATIONS HAVING BANDED JACOBIAN,
C             UCID-30059 REV. 1, L.L.L., MARCH 1975.
C     .
C-----------------------------------------------------------------------
C  DRIVEB IS TO BE CALLED ONCE FOR EACH OUTPUT VALUE OF T, AND
C  IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR, STIFFB.
C     .
C  THE INPUT PARAMETERS ARE..
C    N     =  THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.
C               N CAN BE REDUCED, BUT NEVER INCREASED, DURING A PROBLEM.
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    Y0    =  A VECTOR OF LENGTH N CONTAINING THE INITIAL VALUES OF
C               Y (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 INTRRPOLATE 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 YMAX(I) WILL BE KEPT LESS THAN
C               EPS IN ROOT-MEAN-SQUARE NORM (I.E. EUCLIDEAN NORM
C               DIVIDED BY SQRT(N) ).  THE VECTOR YMAX OF WEIGHTS
C               IS COMPUTED IN DRIVEB.  INITIALLY YMAX(I) IS
C               ABS(Y(I)), WITH A DEFAULT VALUE OF 1 IF Y(I) = 0
C               INITIALLY.  THEREAFTER, YMAX(I) IS THE LARGEST VALUE
C               OF ABS(Y(I)) SEEN SO FAR, OR THE INITIAL YMAX(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 INTRRPOLATION 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    ML,MU =  THE WIDTHS OF THE LOWER AND UPPER PARTS, RESPECTIVELY,
C               OF THE BAND IN THE JACOBIAN MATRIX, NOT COUNTING THE
C               MAIN DIAGONAL.  THE FULL BANDWIDTH IS ML + MU + 1.
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    Y0    =  THE COMPUTED VALUES OF Y 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.  INTRRPOLATION 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     .
C   IN ADDITION TO DRIVEB, THE FOLLOWING ROUTINES ARE PROVIDED IN
C   THE PACKAGE..
C     INTRRP(TOUT,Y,N0,Y0)  INTRRPOLATES TO GET THE OUTPUT VALUES
C                           AT T = TOUT, FROM THE DATA IN THE Y ARRAY.
C     STIFFB(Y,N0)  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     PSETB(Y,N0,CON,MITER,IER)  COMPUTES AND PROCESSES THE JACOBIAN
C                           MATRIX J = DF/DY.
C     DECB(N0,N,ML,MU,B,IP,IER)  AND  SOLB(N0,N,ML,MU,B,C,X,IP)  ARE
C                           USED TO SOLVE BANDED LINEAR SYSTEMS.
C   NOTE.. PSETB, DECB, AND SOLB ARE CALLED ONLY IF MITER = 1 OR 2.
C     .
C   THE FOLLOWING ROUTINES ARE TO BE SUPPLIED BY THE USER..
C     DIFFUN(N,T,Y,YDOT)  COMPUTES THE FUNCTION YDOT = F(Y,T), THE
C                           RIGHT-HAND SIDE OF THE O.D.E.
C                           HERE Y AND YDOT ARE VECTORS OF LENGTH N.
C     PDB(N,T,Y,PD,N0,ML,MU)  COMPUTES THE N BY N JACOBIAN MATRIX OF
C                           PARTIAL DERIVATIVES, AND STORES IT IN PD
C                           AS AN N0 BY ML+MU+1 ARRAY.  PD(I,J-I+ML+1)
C                           IS TO BE SET TO THE PARTIAL DERIVATIVE OF
C                           YDOT(I) WITH RESPECT TO Y(J).  CALLED ONLY
C                           IF MITER = 1.  OTHERWISE A DUMMY ROUTINE
C                           CAN BE SUBSTITUTED.
C     .
C  THE DIMENSIONS IN THE FOLLOWING DECLARATIONS ARE SET FOR A
C  MAXIMUM OF 100 EQUATIONS AND FOR 2*ML+MU .LE. 30.  IF THESE LIMITS
C  ARE TO BE EXCEEDED, THE DIMENSIONS SHOULD BE INCREASED ACCORDINGLY.
C  THE DIMENSION OF PW BELOW MUST BE AT LEAST N*(2*ML+MU+1) IF MITER =
C  1 OR 2, BUT CAN BE REDUCED TO N IF MITER = 3, OR TO 1 IF MITER = 0.
C  THE DIMENSIONS OF YMAX, ERROR, SAVE1, SAVE2, IPIV, AND THE FIRST
C  DIMENSION OF Y SHOULD ALL BE AT LEAST N.  THE COLUMN LENGTH OF
C  THE Y ARRAY AS USED ELSEWHERE IS N0, NOT 100.  THE ROW LENGTH OF Y
C  CAN BE REDUCED FROM 13 TO 6 IF METH = 2.
C  THE IPIV ARRAY IS USED ONLY IF MITER IS 1 OR 2.
C     .
C  THE COMMON BLOCK GEAR9 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 (DIFFUN CALLS) SO FAR,
C  AND THE NUMBER OF JACOBIAN EVALUATIONS SO FAR.
C     .
C  IN THE FOLLOWING DATA STATEMENT, SET..
C    UROUND =  THE UNIT ROUNDOFF OF THE MACHINE, I.E. THE SMALLEST
C              POSITIVE U SUCH THAT 1. + U .NE. 1. ON THE MACHINE.
C    LOUT   =  THE LOGICAL UNIT NUMBER FOR THE OUTPUT OF MESSAGES
C              DURING THE INTEGRATION.
C-----------------------------------------------------------------------
      DIMENSION Y0(N)
      DIMENSION Y(N,13),WKCM(NWKCM),IWCM(NIWCM)
      COMMON /GEAR1/ T,H,HMIN,HMAX,EPSC,UROUND,NC,MFC,KFLAG,JSTART
      DIMENSION YMAX(N),ERROR(N),SAVE1(N),SAVE2(N),PW(N0PW),IPIV(N)
      COMMON /GEAR8/ EPSJ,MLC,MUC,MW,NM1,N0ML,N0W
      COMMON /GEAR9/ HUSED,NQUSED,NSTEP,NFE,NJE
C     .
      COMMON /SYSPAR/ URUND,LDUMMY,LOUT
      EXTERNAL FUNC,BNBDY
C     .
      UROUND=URUND
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.) GO TO 400
      IF (N .LE. 0) GO TO 410
      IF ((T0-TOUT)*H0 .GE. 0.) GO TO 420
C-----------------------------------------------------------------------
C  IF INITIAL VALUES OF YMAX OTHER THAN THOSE SET BELOW ARE DESIRED,
C  THEY SHOULD BE SET HERE.  ALL YMAX(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-----------------------------------------------------------------------
      DO 10 I = 1,N
        YMAX(I) = ABS(Y0(I))
        IF (YMAX(I) .EQ. 0.) YMAX(I) = 1.
 10     Y(I,1) = Y0(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.
      EPSC = EPS
      MFC = MF
      JSTART = 0
      N0 = N
      EPSJ = SQRT(UROUND)
      MLC = ML
      MUC = MU
      MW = ML + MU + 1
      NM1 = N0 - 1
      N0ML = N0*ML
      N0W = N0*MW
      NHCUT = 0
      GO TO 50
C     .
C  TOUTP IS THE PREVIOUS VALUE OF TOUT FOR USE IN HMAX.-----------------
 20   HMAX = ABS(TOUT-TOUTP)*10.
      GO TO 80
C     .
 25   HMAX = ABS(TOUT-TOUTP)*10.
      IF ((T-TOUT)*H .GE. 0.) GO TO 500
      GO TO 85
C     .
 30   IF ((T-TOUT)*H .GE. 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 STIFFB (Y, N0,N0PW,YMAX,ERROR,SAVE1,SAVE2,PW,IPIV,
     1 WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY)
C     .
      IF(NBUG/2-2*(NBUG/4) .NE. 1)GOTO 1055
      WRITE (LOUT,1050) T,NSTEP,H,NFE,NJE,NQUSED
1050  FORMAT(2X,2HT=,E11.4,7H NSTEP=,I4,3H H=,E10.3,8H F EVAL=,
     1 I4,10H JAC EVAL=,I4,7H ORDER=,I2)
1055  CONTINUE
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 YMAX(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, Y0 IS SET TO THE CURRENT Y VALUES ON RETURN.
C  IF INDEX = 2, H IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF
C  ERROR), AND THEN THE CURRENT Y VALUES ARE PUT IN Y0 ON RETURN.
C  FOR ANY OTHER VALUE OF INDEX, CONTROL RETURNS TO THE INTEGRATOR
C  UNLESS TOUT HAS BEEN REACHED.  THEN INTRRPOLATED VALUES OF Y ARE
C  COMPUTED AND STORED IN Y0 ON RETURN.
C  IF INTRRPOLATION IS NOT DESIRED, THE CALL TO INTRRP SHOULD BE
C  REMOVED AND CONTROL TRANSFERRED TO STATEMENT 500 INSTEAD OF 520.
C-----------------------------------------------------------------------
      D = 0.
      DO 70 I = 1,N
        AYI = ABS(Y(I,1))
        YMAX(I) = max(YMAX(I), AYI)
 70     D = D + (AYI/YMAX(I))**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.) GO TO 40
      CALL INTRRP (TOUT, Y, N0, Y0)
      GO TO 520
 85   IF (((T+H)-TOUT)*H .LE. 0.) GO TO 40
      IF (ABS(T-TOUT) .LE. 100.*UROUND*HMAX) GO TO 500
      IF ((T-TOUT)*H .GE. 0.) GO TO 500
      H = (TOUT - T)*(1. - 4.*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   IF(NBUG - 2*(NBUG/2) .NE. 1)GOTO 110
      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 = .1*HMIN
      H = .1*H
      IF(NBUG-2*(NBUG/2) .NE. 1)GOTO 1115
      WRITE (LOUT,115) H
 115  FORMAT(24H  H HAS BEEN REDUCED TO ,E16.8,
     1   26H  AND STEP WILL BE RETRIED//)
1115  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 INTRRPOLATION WAS DONE AS ON NORMAL RETURN./
     2   41H DESIRED PARAMETER CHANGES WERE NOT MADE.)
      CALL INTRRP (TOUT, Y, N0, Y0)
      INDEX = -5
      RETURN
C     .
 500  TOUT = T
      DO 510 I = 1,N
 510    Y0(I) = Y(I,1)
 520  INDEX = KFLAG
      TOUTP = TOUT
      H0 = HUSED
      IF (KFLAG .NE. 0) H0 = H
      RETURN
C----------------------- END OF SUBROUTINE DRIVEB ----------------------
      END
      SUBROUTINE INTRRP (TOUT, Y, N0, Y0)
      DIMENSION Y0(N0),Y(N0,1)
      COMMON /GEAR1/ T,H,DUMMY(4),N,IDUMMY(2),JSTART
C-----------------------------------------------------------------------
C  SUBROUTINE INTRRP COMPUTES INTRRPOLATED VALUES OF THE DEPENDENT
C  VARIABLE Y AND STORES THEM IN Y0.  THE INTRRPOLATION IS TO THE
C  POINT T = TOUT, AND USES THE NORDSIECK HISTORY ARRAY Y, AS FOLLOWS..
C                              NQ
C                   Y0(I)  =  SUM  Y(I,J+1)*S**J ,
C                             J=0
C  WHERE S = -(T-TOUT)/H.
C-----------------------------------------------------------------------
      DO 10 I = 1,N
 10     Y0(I) = Y(I,1)
      L = JSTART + 1
      S = (TOUT - T)/H
      S1 = 1.
      DO 30 J = 2,L
        S1 = S1*S
        DO 20 I = 1,N
 20       Y0(I) = Y0(I) + S1*Y(I,J)
 30     CONTINUE
      RETURN
C----------------------- END OF SUBROUTINE INTRRP ----------------------
      END
      SUBROUTINE STIFFB (Y, N0,N0PW,YMAX,ERROR,SAVE1,SAVE2,PW,IPIV,
     1 WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY)
      DIMENSION Y(N0,13), WKCM(NWKCM),IWCM(NIWCM)
      COMMON /GEAR1/ T,H,HMIN,HMAX,EPS,UROUND,N,MF,KFLAG,JSTART
      DIMENSION YMAX(N0),ERROR(N0),SAVE1(N0),SAVE2(N0),PW(N0PW),
     1 IPIV(N0)
      COMMON /GEAR8/ EPSJ,ML,MU,MW,NM1,N0ML,N0W
      COMMON /GEAR9/ HUSED,NQUSED,NSTEP,NFE,NJE

      EXTERNAL FUNC,BNBDY
C-----------------------------------------------------------------------
C  STIFFB PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE
C  PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.
C  STIFFB IS A VERSION FOR BANDED FORM OF THE JACOBIAN MATRIX.
C  COMMUNICATION WITH STIFFB IS DONE WITH THE FOLLOWING VARIABLES..
C     .
C    Y       AN N0 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              Y(I,J+1) CONTAINS THE J-TH DERIVATIVE OF Y(I), SCALED BY
C              H**J/FACTORIAL(J)  (J = 0,1,...,NQ).
C    N0      A CONSTANT INTEGER .GE. N, USED FOR DIMENSIONING PURPOSES.
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    N       THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.
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 Y 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    YMAX    AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL
C              ERRORS IN Y ARE COMPARED.
C    ERROR   AN ARRAY OF N ELEMENTS.  ERROR(I)/TQ(2) IS THE ESTIMATED
C              ONE-STEP ERROR IN Y(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    IPIV    AN INTEGER ARRAY OF LENGTH N USED FOR PIVOT
C              INFORMATION IF MITER = 1 OR 2.
C    ML,MU   THE LOWER AND UPPER HALF BANDWIDTHS, RESPECTIVELY, OF
C              THE JACOBIAN.  SEE DESCRIPTION IN DRIVER.
C-----------------------------------------------------------------------
      DIMENSION EL(13),TQ(4)
      DATA EL(2)/1./, OLDL0/1./
      KFLAG = 0
      TOLD = T
c
c 6/97 btf: Init local vars mfold,meth,miter,l to avoid f90 warnings
c           re var used before defined:
      if (jstart.eq.0) then
        mfold = 0
        meth = 0
        miter = 0
        l = 0
      endif
      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 YDOT 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 DIFFUN (N, T, Y, SAVE1,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY)
      DO 110 I = 1,N
 110    Y(I,2) = H*SAVE1(I)
      METH = MF/10
      MITER = MF - 10*METH
      NQ = 1
      L = 2
      IDOUB = 3
      RMAX = 1.E4
      RC = 0.
      CRATE = 1.
      EPSOLD = EPS
      HOLD = H
      MFOLD = MF
      NOLD = N
      NSTEP = 0
      NSTEPJ = 0
      NFE = 1
      NJE = 0
      IRET = 3
      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, Y 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  continue
      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)
      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 = max(RH,HMIN/ABS(H))
 175  RH = min(RH,HMAX/ABS(H),RMAX)
      R1 = 1.
      DO 180 J = 2,L
        R1 = R1*RH
        DO 180 I = 1,N
 180      Y(I,J) = Y(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 Y 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  continue
      IF (ABS(RC-1.) .GT. 0.3) 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        Y(I,J) = Y(I,J) + Y(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 Y ARRAY IS NOT ALTERED IN THE CORRECTOR
C  LOOP.  THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1.
C-----------------------------------------------------------------------
 220  DO 230 I = 1,N
 230    ERROR(I) = 0.
      M = 0
      CALL DIFFUN (N, T, Y, SAVE2,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY)
      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.
      NJE = NJE + 1
      NSTEPJ = NSTEP
      GO TO (250, 240, 260), MITER
 240  NFE = NFE + MW
 250  CON = -H*EL(1)

      CALL PSETB (Y, N0, CON, MITER, IER,N0PW,YMAX,SAVE1,SAVE2,PW,IPIV,
     1 WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY)
      IF (IER .NE. 0) GO TO 420
      GO TO 350
 260  R = EL(1)*.1
      DO 270 I = 1,N
 270    PW(I) = Y(I,1) + R*(H*SAVE2(I) - Y(I,2))
      CALL DIFFUN (N, T, PW, SAVE1,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY)
      NFE = NFE + 1
      HL0 = H*EL(1)
      DO 280 I = 1,N
        R0 = H*SAVE2(I) - Y(I,2)
        PW(I) = 1.
        D = .1*R0 - H*(SAVE1(I) - SAVE2(I))
        SAVE1(I) = 0.
        IF (ABS(R0) .LT. UROUND*YMAX(I)) GO TO 280
        IF (ABS(D) .EQ. 0.) GO TO 420
        PW(I) = .1*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 Y DIRECTLY FROM
C  THE RESULT OF THE LAST DIFFUN CALL.
C-----------------------------------------------------------------------
      D = 0.
      DO 300 I = 1,N
        R = H*SAVE2(I) - Y(I,2)
        D = D + ( (R-ERROR(I))/YMAX(I) )**2
        SAVE1(I) = Y(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. - R*(1. - 1./PW(I))
        IF (ABS(D) .EQ. 0.) GO TO 440
 320    PW(I) = 1./D
 330  DO 340 I = 1,N
 340    SAVE1(I) = PW(I)*(H*SAVE2(I) - (Y(I,2) + ERROR(I)))
      GO TO 370
 350  DO 360 I = 1,N
 360    SAVE1(I) = H*SAVE2(I) - (Y(I,2) + ERROR(I))
      CALL SOLB (N0, N, ML, MU, PW, SAVE1, IPIV)
 370  D = 0.
      DO 380 I = 1,N
        ERROR(I) = ERROR(I) + SAVE1(I)
        D = D + (SAVE1(I)/YMAX(I))**2
 380    SAVE1(I) = Y(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 = max(.9*CRATE,D/D1)
      IF ((D*min(1.,2.*CRATE)) .LE. BND) GO TO 450
      D1 = D
      M = M + 1
      IF (M .EQ. 3) GO TO 410
      CALL DIFFUN (N, T, SAVE1, SAVE2,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY)
      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 Y 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.
      DO 430 J1 = 1,NQ
        DO 430 J2 = J1,NQ
          J = (NQ + J1) - J2
          DO 430 I = 1,N
 430        Y(I,J) = Y(I,J) - Y(I,J+1)
      IF (ABS(H) .LE. HMIN*1.00001) GO TO 680
      RH = .25
      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.
      DO 460 I = 1,N
 460    D = D + (ERROR(I)/YMAX(I))**2
      IF (D .GT. E) GO TO 500
C-----------------------------------------------------------------------
C  AFTER A SUCCESSFUL STEP, UPDATE THE Y 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      Y(I,J) = Y(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    Y(I,LMAX) = ERROR(I)
      GO TO 700
C-----------------------------------------------------------------------
C  THE ERROR TEST FAILED.  KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
C  RESTORE T AND THE Y 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        Y(I,J) = Y(I,J) - Y(I,J+1)
      RMAX = 2.
      IF (ABS(H) .LE. HMIN*1.00001) 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.
      DO 530 I = 1,N
 530    D1 = D1 + ((ERROR(I) - Y(I,LMAX))/YMAX(I))**2
      ENQ3 = .5/FLOAT(L+1)
      PR3 = ((D1/EUP)**ENQ3)*1.4 + 1.4E-6
 540  ENQ2 = .5/FLOAT(L)
      PR2 = ((D/E)**ENQ2)*1.2 + 1.2E-6
      PR1 = 1.E+20
      IF (NQ .EQ. 1) GO TO 560
      D = 0.
      DO 550 I = 1,N
 550    D = D + (Y(I,L)/YMAX(I))**2
      ENQ1 = .5/FLOAT(NQ)
      PR1 = ((D/EDN)**ENQ1)*1.3 + 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./PR2
      GO TO 620
 580  NEWQ = NQ - 1
      RH = 1./PR1
      GO TO 620
 590  NEWQ = L
      RH = 1./PR3
      IF (RH .LT. 1.1) GO TO 610
      DO 600 I = 1,N
 600    Y(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.1)) 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 Y 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  Y 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 = .1
      RH = max(HMIN/ABS(H),RH)
      H = H*RH
      CALL DIFFUN (N, T, Y, SAVE1,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY)
      NFE = NFE + 1
      DO 650 I = 1,N
 650    Y(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.
 700  HOLD = H
      JSTART = NQ
      RETURN
C----------------------- END OF SUBROUTINE STIFFB ----------------------
      END
      SUBROUTINE COSET (METH, NQ, EL, TQ, MAXDER)
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 Y 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.,2.,1.,.3158,.07407,.01391,.002182,
     1                 .0002945,.00003492,.000003692,.0000003524,
     2               1.,1.,.5,.1667,.04167,1.,1.,1.,1.,1.,1.,1.,
     3               2.,12.,24.,37.89,53.33,70.08,87.97,106.9,
     4                 126.7,147.4,168.8,191.0,
     5               2.0,4.5,7.333,10.42,13.7,1.,1.,1.,1.,1.,1.,1.,
     6               12.0,24.0,37.89,53.33,70.08,87.97,106.9,
     7                 126.7,147.4,168.8,191.0,1.,
     8               3.0,6.0,9.167,12.5,1.,1.,1.,1.,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.0
      GO TO 900
 102  EL(1) = 0.5
      EL(3) = 0.5
      GO TO 900
 103  EL(1) = 4.1666666666667E-01
      EL(3) = 0.75
      EL(4) = 1.6666666666667E-01
      GO TO 900
 104  EL(1) = 0.375
      EL(3) = 9.1666666666667E-01
      EL(4) = 3.3333333333333E-01
      EL(5) = 4.1666666666667E-02
      GO TO 900
 105  EL(1) = 3.4861111111111E-01
      EL(3) = 1.0416666666667
      EL(4) = 4.8611111111111E-01
      EL(5) = 1.0416666666667E-01
      EL(6) = 8.3333333333333E-03
      GO TO 900
 106  EL(1) = 3.2986111111111E-01
      EL(3) = 1.1416666666667
      EL(4) = 0.625
      EL(5) = 1.7708333333333E-01
      EL(6) = 0.025
      EL(7) = 1.3888888888889E-03
      GO TO 900
 107  EL(1) = 3.1559193121693E-01
      EL(3) = 1.225
      EL(4) = 7.5185185185185E-01
      EL(5) = 2.5520833333333E-01
      EL(6) = 4.8611111111111E-02
      EL(7) = 4.8611111111111E-03
      EL(8) = 1.9841269841270E-04
      GO TO 900
 108  EL(1) = 3.0422453703704E-01
      EL(3) = 1.2964285714286
      EL(4) = 8.6851851851852E-01
      EL(5) = 3.3576388888889E-01
      EL(6) = 7.7777777777778E-02
      EL(7) = 1.0648148148148E-02
      EL(8) = 7.9365079365079E-04
      EL(9) = 2.4801587301587E-05
      GO TO 900
 109  EL(1) = 2.9486800044092E-01
      EL(3) = 1.3589285714286
      EL(4) = 9.7655423280423E-01
      EL(5) = 0.4171875
      EL(6) = 1.1135416666667E-01
      EL(7) = 0.01875
      EL(8) = 1.9345238095238E-03
      EL(9) = 1.1160714285714E-04
      EL(10)= 2.7557319223986E-06
      GO TO 900
 110  EL(1) = 2.8697544642857E-01
      EL(3) = 1.4144841269841
      EL(4) = 1.0772156084656
      EL(5) = 4.9856701940035E-01
      EL(6) = 0.1484375
      EL(7) = 2.9060570987654E-02
      EL(8) = 3.7202380952381E-03
      EL(9) = 2.9968584656085E-04
      EL(10)= 1.3778659611993E-05
      EL(11)= 2.7557319223986E-07
      GO TO 900
 111  EL(1) = 2.8018959644394E-01
      EL(3) = 1.4644841269841
      EL(4) = 1.1715145502646
      EL(5) = 5.7935819003527E-01
      EL(6) = 1.8832286155203E-01
      EL(7) = 4.1430362654321E-02
      EL(8) = 6.2111441798942E-03
      EL(9) = 6.2520667989418E-04
      EL(10)= 4.0417401528513E-05
      EL(11)= 1.5156525573192E-06
      EL(12)= 2.5052108385442E-08
      GO TO 900
 112  EL(1) = 2.7426554003160E-01
      EL(3) = 1.5099386724387
      EL(4) = 1.2602711640212
      EL(5) = 6.5923418209877E-01
      EL(6) = 2.3045800264550E-01
      EL(7) = 5.5697246105232E-02
      EL(8) = 9.4394841269841E-03
      EL(9) = 1.1192749669312E-03
      EL(10)= 9.0939153439153E-05
      EL(11)= 4.8225308641975E-06
      EL(12)= 1.5031265031265E-07
      EL(13)= 2.0876756987868E-09
      GO TO 900
 201  EL(1) = 1.0
      GO TO 900
 202  EL(1) = 6.6666666666667E-01
      EL(3) = 3.3333333333333E-01
      GO TO 900
 203  EL(1) = 5.4545454545455E-01
      EL(3) = EL(1)
      EL(4) = 9.0909090909091E-02
      GO TO 900
 204  EL(1) = 0.48
      EL(3) = 0.7
      EL(4) = 0.2
      EL(5) = 0.02
      GO TO 900
 205  EL(1) = 4.3795620437956E-01
      EL(3) = 8.2116788321168E-01
      EL(4) = 3.1021897810219E-01
      EL(5) = 5.4744525547445E-02
      EL(6) = 3.6496350364964E-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 PSETB (Y, N0, CON, MITER, IER,N0PW,YMAX,SAVE1,SAVE2,
     1 PW,IPIV,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY)
      DIMENSION Y(N0,1),WKCM(NWKCM),IWCM(NIWCM)
      COMMON /GEAR1/ T,H,DUMMY(3),UROUND,N,IDUMMY(3)
      DIMENSION YMAX(N0),SAVE1(N0),SAVE2(N0),PW(N0PW),IPIV(N0)
      COMMON /GEAR8/ EPSJ,ML,MU,MW,NM1,N0ML,N0W
      EXTERNAL FUNC,BNBDY
C-----------------------------------------------------------------------
C  PSETB IS CALLED BY STIFFB TO COMPUTE AND PROCESS THE MATRIX
C  P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
C  THIS IS THE VERSION FOR BANDED FORM OF J.
C  J IS COMPUTED, EITHER BY THE USER-SUPPLIED ROUTINE PDB
C  IF MITER = 1, OR BY FINITE DIFFERENCING IF MITER = 2.
C  J IS STORED IN PW AND REPLACED BY P, USING CON = -H*EL(1).
C  THEN P IS SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR
C  LATER SOLUTION OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX.
C     .
C  IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
C  WITH PSETB USES THE FOLLOWING..
C    EPSJ    = SQRT(UROUND), USED IN THE NUMERICAL JACOBIAN INCREMENTS.
C    MW      = ML + MU + 1.
C    NM1     = N0 - 1.
C    N0ML    = N0*ML.
C    N0W     = N0*MW.
C-----------------------------------------------------------------------
      IF (MITER .EQ. 2) GO TO 20
C  IF MITER = 1, CALL PDB AND MULTIPLY BY SCALAR. ----------------------
! 6/14/05 btf: SGI compiler does not like this call because the PDB routine
!              defined above is a "dummy" sub, and has no dummy args.
!              PDB is not called here (MITER==2), so is commented out.
!     CALL PDB (N, T, Y, PW, N0, ML, MU)
      DO 10 I = 1,N0W
 10     PW(I) = PW(I)*CON
      GO TO 90
C  IF MITER = 2, MAKE MW CALLS TO DIFFUN TO APPROXIMATE J. -------------
 20   D = 0.
      DO 30 I = 1,N
 30     D = D + SAVE2(I)**2
      R0 = ABS(H)*SQRT(D)*1.E03*UROUND
      DO 80 J = 1,MW
        KMAX = (N - J)/MW + 1
        YM = 0.
        DO 40 K = 1,KMAX
          I = J + (K - 1)*MW
 40       YM = max(YM,YMAX(I))
        R = max(EPSJ*YM,R0)
        DO 50 K = 1,KMAX
          I = J + (K - 1)*MW
 50       Y(I,1) = Y(I,1) + R
        CALL DIFFUN (N, T, Y, SAVE1,WKCM,IWCM,NWKCM,NIWCM,FUNC,BNBDY)
        D = CON/R
        J1 = J*N0 + N0ML
        DO 70 K = 1,KMAX
          JJ = J + (K - 1)*MW
          I1 = MAX0(JJ-MU,1)
          I2 = MIN0(JJ+ML,N)
          II = J1 - NM1*I1
          DO 60 I = I1,I2
            PW(II) = (SAVE1(I) - SAVE2(I))*D
 60         II = II - NM1
          Y(JJ,1) = Y(JJ,1) - R
 70       J1 = J1 + N0W
 80     CONTINUE
C  ADD IDENTITY MATRIX. ------------------------------------------------
 90   DO 100 I = 1,N
 100    PW(N0ML+I) = PW(N0ML+I) + 1.
C  DO LU DECOMPOSITION ON P. -------------------------------------------
      CALL DECB (N0, N, ML, MU, PW, IPIV, IER)
      RETURN
C----------------------- END OF SUBROUTINE PSETB -----------------------
      END
      SUBROUTINE DECB (NDIM, N, ML, MU, B, IP, IER)
      DIMENSION B(NDIM,1),IP(N)
C-----------------------------------------------------------------------
C     LU DECOMPOSITION OF BAND MATRIX A..  L*U = P*A , WHERE P IS A
C        PERMUTATION MATRIX, L IS A UNIT LOWER TRIANGULAR MATRIX,
C        AND U IS AN UPPER TRIANGULAR MATRIX.
C     N     =  ORDER OF MATRIX.
C     B     =  N BY (2*ML+MU+1) ARRAY CONTAINING THE MATRIX A ON INPUT
C              AND ITS FACTORED FORM ON OUTPUT.
C              ON INPUT, B(I,K) (1.LE.I.LE.N) CONTAINS THE K-TH
C              DIAGONAL OF A, OR A(I,J) IS STORED IN B(I,J-I+ML+1).
C              ON OUTPUT, B CONTAINS THE L AND U FACTORS, WITH
C              U IN COLUMNS 1 TO ML+MU+1, AND L IN COLUMNS
C              ML+MU+2 TO 2*ML+MU+1.
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     NDIM  =  THE FIRST DIMENSION (COLUMN LENGTH) OF THE ARRAY B.
C              NDIM MUST BE .GE. N.
C     IP    =  ARRAY OF LENGTH N CONTAINING PIVOT INFORMATION.
C     IER   =  ERROR INDICATOR..
C           =  0  IF NO ERROR,
C           =  K  IF THE K-TH PIVOT CHOSEN WAS ZERO (A IS SINGULAR).
C     THE INPUT ARGUMENTS ARE  NDIM, N, ML, MU, B.
C     THE OUTPUT ARGUMENTS ARE  B, IP, IER.
C-----------------------------------------------------------------------
      IER = 0
      IF (N .EQ. 1) GO TO 92
      LL = ML + MU + 1
      N1 = N - 1
      IF (ML .EQ. 0) GO TO 32
      DO 30 I = 1,ML
        II = MU + I
        K = ML + 1 - I
        DO 10 J = 1,II
 10       B(I,J) = B(I,J+K)
        K = II + 1
        DO 20 J = K,LL
 20       B(I,J) = 0.
 30     CONTINUE
 32   LR = ML
      DO 90 NR = 1,N1
        NP = NR + 1
        IF (LR .NE. N) LR = LR + 1
        MX = NR
        XM = ABS(B(NR,1))
        IF (ML .EQ. 0) GO TO 42
        DO 40 I = NP,LR
          IF (ABS(B(I,1)) .LE. XM) GO TO 40
          MX = I
          XM = ABS(B(I,1))
 40       CONTINUE
 42     IP(NR) = MX
        IF (MX .EQ. NR) GO TO 60
        DO 50 I = 1,LL
          XX = B(NR,I)
          B(NR,I) = B(MX,I)
 50       B(MX,I) = XX
 60     XM = B(NR,1)
        IF (XM .EQ. 0.) GO TO 100
        B(NR,1) = 1./XM
        IF (ML .EQ. 0) GO TO 90
        XM = -B(NR,1)
        KK = MIN0(N-NR,LL-1)
        DO 80 I = NP,LR
          J = LL + I - NR
          XX = B(I,1)*XM
          B(NR,J) = XX
          DO 70 II = 1,KK
 70         B(I,II) = B(I,II+1) + XX*B(NR,II+1)
 80       B(I,LL) = 0.
 90     CONTINUE
 92   NR = N
      IF (B(N,1) .EQ. 0.) GO TO 100
      B(N,1) = 1./B(N,1)
      RETURN
 100  IER = NR
      RETURN
C----------------------- END OF SUBROUTINE DECB ------------------------
      END
      SUBROUTINE SOLB (NDIM, N, ML, MU, B, Y, IP)
      DIMENSION B(NDIM,1),Y(N),IP(N)
C-----------------------------------------------------------------------
C     SOLUTION OF  A*X = C  GIVEN LU DECOMPOSITION OF A FROM DECB.
C     Y  =  RIGHT-HAND VECTOR C, OF LENGTH N, ON INPUT,
C        =  SOLUTION VECTOR X ON OUTPUT.
C     ALL THE ARGUMENTS ARE INPUT ARGUMENTS.
C     THE OUTPUT ARGUMENT IS  Y.
C-----------------------------------------------------------------------
      IF (N .EQ. 1) GO TO 60
      N1 = N - 1
      LL = ML + MU + 1
      IF (ML .EQ. 0) GO TO 32
      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
 20       Y(NR+I) = Y(NR+I) + Y(NR)*B(NR,LL+I)
 30     CONTINUE
 32   LL = LL - 1
      Y(N) = Y(N)*B(N,1)
      KK = 0
      DO 50 NB = 1,N1
        NR = N - NB
        IF (KK .NE. LL) KK = KK + 1
        DP = 0.
        IF (LL .EQ. 0) GO TO 50
        DO 40 I = 1,KK
 40       DP = DP + B(NR,I+1)*Y(NR+I)
 50     Y(NR) = (Y(NR) - DP)*B(NR,1)
      RETURN
 60   Y(1) = Y(1)*B(1,1)
      RETURN
      END
c
C----------------------- END OF SUBROUTINE SOLB ------------------------

