*DECK DCFODE
SUBROUTINE DCFODE (METH, ELCO, TESCO)
C***BEGIN PROLOGUE DCFODE
C***SUBSIDIARY
C***PURPOSE Set ODE integrator coefficients.
C***TYPE DOUBLE PRECISION (SCFODE-S, DCFODE-D)
C***AUTHOR Hindmarsh, Alan C., (LLNL)
C***DESCRIPTION
C
C DCFODE is called by the integrator routine to set coefficients
C needed there. The coefficients for the current method, as
C given by the value of METH, are set for all orders and saved.
C The maximum order assumed here is 12 if METH = 1 and 5 if METH = 2.
C (A smaller value of the maximum order is also allowed.)
C DCFODE is called once at the beginning of the problem,
C and is not called again unless and until METH is changed.
C
C The ELCO array contains the basic method coefficients.
C The coefficients el(i), 1 .le. i .le. nq+1, for the method of
C order nq are stored in ELCO(i,nq). They are given by a genetrating
C polynomial, i.e.,
C l(x) = el(1) + el(2)*x + ... + el(nq+1)*x**nq.
C For the implicit Adams methods, l(x) is given by
C dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1), l(-1) = 0.
C For the BDF methods, l(x) is given by
C l(x) = (x+1)*(x+2)* ... *(x+nq)/K,
C where K = factorial(nq)*(1 + 1/2 + ... + 1/nq).
C
C The TESCO array contains test constants used for the
C local error test and the selection of step size and/or order.
C At order nq, TESCO(k,nq) is used for the selection of step
C size at order nq - 1 if k = 1, at order nq if k = 2, and at order
C nq + 1 if k = 3.
C
C***SEE ALSO DLSODE
C***ROUTINES CALLED (NONE)
C***REVISION HISTORY (YYMMDD)
C 791129 DATE WRITTEN
C 890501 Modified prologue to SLATEC/LDOC format. (FNF)
C 890503 Minor cosmetic changes. (FNF)
C 930809 Renamed to allow single/double precision versions. (ACH)
C***END PROLOGUE DCFODE
C**End
INTEGER METH
INTEGER I, IB, NQ, NQM1, NQP1
DOUBLE PRECISION ELCO, TESCO
DOUBLE PRECISION AGAMQ, FNQ, FNQM1, PC, PINT, RAGQ,
1 RQFAC, RQ1FAC, TSIGN, XPIN
DIMENSION ELCO(13,12), TESCO(3,12)
DIMENSION PC(12)
C
C***FIRST EXECUTABLE STATEMENT DCFODE
GO TO (100, 200), METH
C
100 ELCO(1,1) = 1.0D0
ELCO(2,1) = 1.0D0
TESCO(1,1) = 0.0D0
TESCO(2,1) = 2.0D0
TESCO(1,2) = 1.0D0
TESCO(3,12) = 0.0D0
PC(1) = 1.0D0
RQFAC = 1.0D0
DO 140 NQ = 2,12
C-----------------------------------------------------------------------
C The PC array will contain the coefficients of the polynomial
C p(x) = (x+1)*(x+2)*...*(x+nq-1).
C Initially, p(x) = 1.
C-----------------------------------------------------------------------
RQ1FAC = RQFAC
RQFAC = RQFAC/NQ
NQM1 = NQ - 1
FNQM1 = NQM1
NQP1 = NQ + 1
C Form coefficients of p(x)*(x+nq-1). ----------------------------------
PC(NQ) = 0.0D0
DO 110 IB = 1,NQM1
I = NQP1 - IB
110 PC(I) = PC(I-1) + FNQM1*PC(I)
PC(1) = FNQM1*PC(1)
C Compute integral, -1 to 0, of p(x) and x*p(x). -----------------------
PINT = PC(1)
XPIN = PC(1)/2.0D0
TSIGN = 1.0D0
DO 120 I = 2,NQ
TSIGN = -TSIGN
PINT = PINT + TSIGN*PC(I)/I
120 XPIN = XPIN + TSIGN*PC(I)/(I+1)
C Store coefficients in ELCO and TESCO. --------------------------------
ELCO(1,NQ) = PINT*RQ1FAC
ELCO(2,NQ) = 1.0D0
DO 130 I = 2,NQ
130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/I
AGAMQ = RQFAC*XPIN
RAGQ = 1.0D0/AGAMQ
TESCO(2,NQ) = RAGQ
IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/NQP1
TESCO(3,NQM1) = RAGQ
140 CONTINUE
RETURN
C
200 PC(1) = 1.0D0
RQ1FAC = 1.0D0
DO 230 NQ = 1,5
C-----------------------------------------------------------------------
C The PC array will contain the coefficients of the polynomial
C p(x) = (x+1)*(x+2)*...*(x+nq).
C Initially, p(x) = 1.
C-----------------------------------------------------------------------
FNQ = NQ
NQP1 = NQ + 1
C Form coefficients of p(x)*(x+nq). ------------------------------------
PC(NQP1) = 0.0D0
DO 210 IB = 1,NQ
I = NQ + 2 - IB
210 PC(I) = PC(I-1) + FNQ*PC(I)
PC(1) = FNQ*PC(1)
C Store coefficients in ELCO and TESCO. --------------------------------
DO 220 I = 1,NQP1
220 ELCO(I,NQ) = PC(I)/PC(2)
ELCO(2,NQ) = 1.0D0
TESCO(1,NQ) = RQ1FAC
TESCO(2,NQ) = NQP1/ELCO(1,NQ)
TESCO(3,NQ) = (NQ+2)/ELCO(1,NQ)
RQ1FAC = RQ1FAC/FNQ
230 CONTINUE
RETURN
C----------------------- END OF SUBROUTINE DCFODE ----------------------
END