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