1
2      SUBROUTINE CFOD(METH,ELCO,TESCO)
3C***BEGIN PROLOGUE  CFOD
4C***REFER TO  DEBDF
5C
6C   CFOD defines coefficients needed in the integrator package DEBDF
7C***ROUTINES CALLED  (NONE)
8C***END PROLOGUE  CFOD
9C
10C
11CLLL. OPTIMIZE
12      INTEGER METH, I, IB, NQ, NQM1, NQP1
13      REAL ELCO, TESCO, AGAMQ, FNQ, FNQM1, PC, PINT, RAGQ,
14     1   RQFAC, RQ1FAC, TSIGN, XPIN
15      DIMENSION ELCO(13,12), TESCO(3,12)
16C-----------------------------------------------------------------------
17C CFOD  IS CALLED BY THE INTEGRATOR ROUTINE TO SET COEFFICIENTS
18C NEEDED THERE.  THE COEFFICIENTS FOR THE CURRENT METHOD, AS
19C GIVEN BY THE VALUE OF METH, ARE SET FOR ALL ORDERS AND SAVED.
20C THE MAXIMUM ORDER ASSUMED HERE IS 12 IF METH = 1 AND 5 IF METH = 2.
21C (A SMALLER VALUE OF THE MAXIMUM ORDER IS ALSO ALLOWED.)
22C CFOD  IS CALLED ONCE AT THE BEGINNING OF THE PROBLEM,
23C AND IS NOT CALLED AGAIN UNLESS AND UNTIL METH IS CHANGED.
24C
25C THE ELCO ARRAY CONTAINS THE BASIC METHOD COEFFICIENTS.
26C THE COEFFICIENTS EL(I), 1 .LE. I .LE. NQ+1, FOR THE METHOD OF
27C ORDER NQ ARE STORED IN ELCO(I,NQ).  THEY ARE GIVEN BY A GENETRATING
28C POLYNOMIAL, I.E.,
29C     L(X) = EL(1) + EL(2)*X + ... + EL(NQ+1)*X**NQ.
30C FOR THE IMPLICIT ADAMS METHODS, L(X) IS GIVEN BY
31C     DL/DX = (X+1)*(X+2)*...*(X+NQ-1)/FACTORIAL(NQ-1),    L(-1) = 0.
32C FOR THE BDF METHODS, L(X) IS GIVEN BY
33C     L(X) = (X+1)*(X+2)* ... *(X+NQ)/K,
34C WHERE         K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ).
35C
36C THE TESCO ARRAY CONTAINS TEST CONSTANTS USED FOR THE
37C LOCAL ERROR TEST AND THE SELECTION OF STEP SIZE AND/OR ORDER.
38C AT ORDER NQ, TESCO(K,NQ) IS USED FOR THE SELECTION OF STEP
39C SIZE AT ORDER NQ - 1 IF K = 1, AT ORDER NQ IF K = 2, AND AT ORDER
40C NQ + 1 IF K = 3.
41C-----------------------------------------------------------------------
42      DIMENSION PC(12)
43C
44C***FIRST EXECUTABLE STATEMENT  CFOD
45      GO TO (100, 200), METH
46C
47 100  ELCO(1,1) = 1.0E0
48      ELCO(2,1) = 1.0E0
49      TESCO(1,1) = 0.0E0
50      TESCO(2,1) = 2.0E0
51      TESCO(1,2) = 1.0E0
52      TESCO(3,12) = 0.0E0
53      PC(1) = 1.0E0
54      RQFAC = 1.0E0
55      DO 140 NQ = 2,12
56C-----------------------------------------------------------------------
57C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL
58C     P(X) = (X+1)*(X+2)*...*(X+NQ-1).
59C INITIALLY, P(X) = 1.
60C-----------------------------------------------------------------------
61        RQ1FAC = RQFAC
62        RQFAC = RQFAC/FLOAT(NQ)
63        NQM1 = NQ - 1
64        FNQM1 = FLOAT(NQM1)
65        NQP1 = NQ + 1
66C FORM COEFFICIENTS OF P(X)*(X+NQ-1). ----------------------------------
67        PC(NQ) = 0.0E0
68        DO 110 IB = 1,NQM1
69          I = NQP1 - IB
70 110      PC(I) = PC(I-1) + FNQM1*PC(I)
71        PC(1) = FNQM1*PC(1)
72C COMPUTE INTEGRAL, -1 TO 0, OF P(X) AND X*P(X). -----------------------
73        PINT = PC(1)
74        XPIN = PC(1)/2.0E0
75        TSIGN = 1.0E0
76        DO 120 I = 2,NQ
77          TSIGN = -TSIGN
78          PINT = PINT + TSIGN*PC(I)/FLOAT(I)
79 120      XPIN = XPIN + TSIGN*PC(I)/FLOAT(I+1)
80C STORE COEFFICIENTS IN ELCO AND TESCO. --------------------------------
81        ELCO(1,NQ) = PINT*RQ1FAC
82        ELCO(2,NQ) = 1.0E0
83        DO 130 I = 2,NQ
84 130      ELCO(I+1,NQ) = RQ1FAC*PC(I)/FLOAT(I)
85        AGAMQ = RQFAC*XPIN
86        RAGQ = 1.0E0/AGAMQ
87        TESCO(2,NQ) = RAGQ
88      IF(NQ.LT.12)TESCO(1,NQP1)=RAGQ*RQFAC/FLOAT(NQP1)
89        TESCO(3,NQM1) = RAGQ
90 140    CONTINUE
91      RETURN
92C
93 200  PC(1) = 1.0E0
94      RQ1FAC = 1.0E0
95      DO 230 NQ = 1,5
96C-----------------------------------------------------------------------
97C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL
98C     P(X) = (X+1)*(X+2)*...*(X+NQ).
99C INITIALLY, P(X) = 1.
100C-----------------------------------------------------------------------
101        FNQ = FLOAT(NQ)
102        NQP1 = NQ + 1
103C FORM COEFFICIENTS OF P(X)*(X+NQ). ------------------------------------
104        PC(NQP1) = 0.0E0
105        DO 210 IB = 1,NQ
106          I = NQ + 2 - IB
107 210      PC(I) = PC(I-1) + FNQ*PC(I)
108        PC(1) = FNQ*PC(1)
109C STORE COEFFICIENTS IN ELCO AND TESCO. --------------------------------
110        DO 220 I = 1,NQP1
111 220      ELCO(I,NQ) = PC(I)/PC(2)
112        ELCO(2,NQ) = 1.0E0
113        TESCO(1,NQ) = RQ1FAC
114        TESCO(2,NQ) = FLOAT(NQP1)/ELCO(1,NQ)
115        TESCO(3,NQ) = FLOAT(NQ+2)/ELCO(1,NQ)
116        RQ1FAC = RQ1FAC/FNQ
117 230    CONTINUE
118      RETURN
119C----------------------- END OF SUBROUTINE CFOD  -----------------------
120      END
121