cwharman
Mechanical
- Jul 10, 2005
- 1
Hi,
I am getting a run-time error when trying to compile my program and as I
am not very experienced in Fortran I am not sure what the error is or how
to fix it. The program has been lifted out of a book called "Interneal
Combustion Engines:Applied Thermosciences" and is a program to model the
heat release rate of a an engine. It requires a couplke of subroutines,
DVERK and LEQIF which are IMSL routines and I do not have. I have found
the code for DVERK but unfortunately have not found the code for LEQIF
which is why I assume I am receiving the run-time error. the error is a
zero or negative argument to a logarithm routine.
FARG - in file release2_good.for at line 664 [+1738]
main - in file release2_good.for at line 41 [+0190]
The main code is listed below and any help in resolving this matter would
be greatly appreciated.
Thanks
Chris
C GIVEN IN BLOCK DATA SUBPROGRAM
C COMPRESSION RATIO - R
C BORE - B (CM)
C STROKE - S (CM)
C HALF-STROKE TO ROD RATIO - EPS
C ENGINE SPEED - RPM
C HEAT TRANSFER COEFFICIENT - H (J/S/M**2/K)
C BLOWBY COEFFICIENT - C (L/S)
C EQUIVALENCE RATIO - PHI
C RESIDUAL FRACTION - F
C INITIAL PRESSURE - P1 (BAR)
C INITIAL TEMPERATURE - T1 (K)
C WALL TEMPERATURE - TW (K)
C BURN ANGLE - THETAB (DEG)
C START OF HEAT RELEASE - THETAS (DEG ATDC)
C
C FIND:
C INDICATED THERMAL EFFICIENCY - ETA
C INDICATED MEAN EFFECTIVE PRESSURE - IMEP (BAR)
C
C NOTES:
C 1. COSINE BURNING LAW IS EMPLOYED, EQN. 4.61
C ` 2. FUEL - GASOLINE (C7H17)
C
EXTERNAL RATES
REAL IMEP, MNOT, M
COMMON /ENGINE/ R, B, S, EPS, RPM, H, C, THETAB, THETAS,
& PHI, F, P1, T1, TW, MNOT, OMEGA, UNOT, UFINAL
DIMENSION Y1(6), Y2(10), Y(6), COM(24), WM(6,9)
DATA AO/47870./, FS/.06548/, Y/1., 10000., 350., 3*0.0/,
& N/6/, THETA/-180./, THETAE/-179./, TOL/.001/, IND/1/,
& NW/6/, PI/3.141593/, Y2/10*0.0/
C
C LABEL HEADER ON OUTPUT PAGE AND PRINT INITIAL CONDITIONS
C
OMEGA = RPM*PI/30
CALL AUXLRY(THETA, V, X, EM)
Y(1) = P1
Y(3) = T1
CALL FARG(Y(1), Y(3), PHI, F, HU, U, VU, SU, Y1, CP, DLVLT,
& DLVLP)
MNOT = V/VU
UNOT = MNOT*U
WRITE(6,10)
10 FORMAT(5X,'THETA',4X,'VOLUME',3X,'BURN FRAC',2X,'PRESSURE',3X,
& 'BURN TEMP',2X,'T UNBURN',3X,'WORK',3X,'HEAT LOSS',5X,
& 'MASS',3X,'H-LEAK')
WRITE(6,15)
15 FORMAT(4X,'DEG ATDC',2X,'CM**3',6X,'--',8X,'BARS',7X,'KELVIN',
& 5X,'KELVIN',3X,'JOULES',4X,'JOULES',6X,'GRAMS',2X,
& 'JOULES')
WRITE(6,20) THETA, V, X, (Y(I),I=1,5), MNOT, Y(6)
20 FORMAT(4X,F5.0,6X,F5.0,4X,F5.3,6X,F5.2,7X,F5.0,5X,F5.0,
& 5X,F5.0,3X,F5.0,7X,F5.3,3X,F5.2)
C
C INTEGRATE THE RATE EQUATIONS 4.80 - 4.85 AND PRINT ANSWERS
C EVERY 10 DEGREES
C
DO 30 I = 1,36
DO 25 J = 1,10
CALL DVERK (N, RATES, THETA, Y, THETAE, TOL, IND, COM, NW, WM)
CALL AUXLRY(THETAE, V, X, EM)
M = EM*MNOT
THETA = THETAE
THETAE = THETA + 1.
IF(THETAS .GE. THETA .AND. THETAS .LT. THETAE)
& CALL TINITL(Y(1), T(3), PHI, F, Y(2))
IF(THETA .GT. THETAS + THETAB) Y(3) = 10000.
25 CONTINUE
WRITE(6,20) THETA, V, X, (Y(J),J=1,5), M, Y(6)
30 CONTINUE
C
C COMPUTE AND WRITE THE EAT, IMEP AND ERRORS
C
CALL ECP(Y(1), Y(2), PHI, HB, U, VB, SB, Y2, CP, DLVLT, DLVLP,
& IER)
UFINAL = U*M
ERROR1 = 1.0 - VB*M/V
ERROR2 = 1.0 + Y(4)/(UFINAL - UNOT + Y(5) + Y(6))
ETA = Y(4)/MNOT*(1. + PHI*FS*(1. - F))/PHI/FS/(1. - F)/AO
IMEP = Y(4)/(PI/4.*B**2*S)*10.
WRITE(6,40) ETA, IMEP, ERROR1, ERROR2
40 FORMAT('ETA = ',1PE11.4,' IMEP = ',1PE11.4,'ERROR1 = ',1PE11.4,
& ' ERROR2 = ',1PE11.4)
STOP
END
C
SUBROUTINE RATES(N, THETA, Y, YPRIME)
DIMENSION Y1(6), Y2(10)
REAL Y(N), YPRIME(N), M, MNOT
COMMON /ENGINE/ R, BORE, STROKE, EPS, RPM, HEAT, CEE, THETAB,
& THETAS, PHI, F, P1, T1, TW, MNOT, OMEGA, UNOT,
& UFINAL
DATA PI/3.141593/, Y2/10*0.0/
CALL AUXLRY(THETA, V, X, EM)
M = EM*MNOT
RAD = THETA*PI/180.
DUMB = SQRT(1. - (EPS*SIN(RAD))**2)
DV = PI/8.*BORE**2*STROKE*SIN(RAD)*(1. + EPS*COS(RAD)/DUMB)
A = (DV + V*CEE/OMEGA)/M
C1 = HEAT*(PI*BORE**2/2. + 4.*V/BORE)/OMEGA/M/10000.
C0 = SQRT(X)
C
C THREE DIFFERENT COMPUTATIONS ARE REQUIRED DEPENDING
C UPON THE SIZE OF THE MASS FRACTION BURNED
C
IF( X .GT. .999 ) GO TO 20
IF( X .GT. .001 ) GO TO 10
C
C COMPRESSION
C
CALL FARG(Y(1), Y(3), PHI, F, HL, U, VU, S, Y1, CP, DLVLT,
& DLVLP)
B = C1*VU/CP*DLVLT*(1.0 - TW/Y(3))
C = 0.
D = 0.
E = VU**2/CP/Y(3)*DLVLT**2 + 10.*VU/Y(1)*DLVLP
YPRIME(1) = (A + B + C)/(D + E)*10.
YPRIME(2) = 0.
YPRIME(3) = -C1/CP*(Y(3) - TW) + VU/CP*DLVLT*YPRIME(1)/10.
GO TO 30
C
C COMBUSTION
C
10 CALL FARG(Y(1), Y(3), PHI, F, HU, U, VU, S, Y1, CPU, DLVLTU,
& DLVLP)
CALL ECP(Y(1), Y(2), PHI, HB, U, VB, S, Y2, CPB, DDLVLTB,
& DLVLPB, IER)
B = C1*(VB/CPB*DLVLTB*C0*(1. - TW/Y(2)) + VU/CPU*DLVLTU*
& (1. - C0)*(1.- TW/Y(3)))
DX = 0.5*SIN(PI*(THETA - THETAS)/THETAB)*180./THETAB
C = -(VB - VU)*DX - VB/Y(2)*DLVLTB*(HU - HB)/CPB*(DX - (X -
& X**2)*CEE/OMEGA)
D = X*(VB**2/CPB/Y(2)*DLVLTB**2 + 10.*VB/Y(1)*DLVLPB)
E = (1. - X)*(VU**2/CPU/Y(3)*DLVLTU**2 + 10.*VU/Y(1)*DLVLPU)
HL = (1.0 - X**2)*HU + X**2*HB
YPRIME(1) = (A + B + C)/(D + E)*10.
YPRIME(2) = -C1/CPB/C0*(Y(2) - TW) + VB/CPB*DLVLTB*
& YPRIME(1)/10. + (HU - HB)/CPB*(DX/X - (1. - X)*
& CEE/OMEGA)
YPRIME(3) = -C1/CPU/(1. + CO)*(Y(3) - TW) + VU/CPU*DLVLTU*
& YPRIME(1)/10.
GO TO 30
C
C EXPANSION
C
20 CALL ECP(Y(1), Y(2), PHI, HL, U, VB, S, Y2, CP, DLVLT, DLVLP,
& IER)
B = C1*VB/CP*DLVLT*(1. - TW/Y(2))
C = 0.
D = VB**2/CP/Y(2)*DLVLT**2 + 10.*VB/Y(1)*DLVLP
E = 0.
YPRIME(1) = (A + B + C)/(D + E)*10.
YPRIME(2) = -C1/CP*(Y(2) - TW) + VB/CP*DLVLT*YPRIME(1)/10.
YPRIME(3) = 0.
C
C ALL CASES
C
30 YPRIME(4) = Y(1)*DV/10.
YPRIME(5) = C1*M*(C0*(Y(2) - TW) + (1. - C0)*(Y(3) - TW))
YPRIME(6) = CEE*M/OMEGA*HL
DO 40 I = 1,N
40 YPRIME(I) = YPRIME(I)*PI/180.
RETURN
END
C
SUBROUTINE AUXLRY(THETA, V, X, EM)
COMMON /ENGINE/ R, B, S, EPS, RPM, H, C, THETAB, THETAS,
& PHI, F, P1, T1, TW, EMNOT, OMEGA, UNOT, UFINAL
DATA PI/3.141593/
RAD = THETA*PI/180.
VTDC = PI/4.*B**2*S/(R - 1.)
V = VTDC*(1. + (R - 1.)/2.*(1. - COS(RAD) + 1./EPS*(1. -
& SQRT(1. - (EPS*SIN(RAD))**2))))
X = 0.5*(1. - COS(PI*(THETA - THETAS)/THETAB))
IF(THETA .LE. THETAS) X = 0.
IF(THETA .GE. THETAS + THETAB) X = 1.
EM = EXP(-C*(RAD + PI)/OMEGA)
RETURN
END
C
BLOCK DATA
COMMON /ENGINE/ R, B, S, EPS, RPM, H, C, THETAB, THETAS,
& PHI, F, P1, T1, TW, EMNOT, OMEGA, UNOT, UFINAL
DATA R/10./, B/10./, S/8./, EPS/.25/, RPM/2000./, H/500./,
& C/0.8/, THETAB/60./, THETAS/-35./, PHI/0.8/, F/.1/,
& P1/1.0/, T1/350./, TW/420./
C
END
C
SUBROUTINE TINITL(P, TU, PHI, F, TB)
DIMENSION Y1(6), Y2(10)
DATA MAXITS/50/, TOL/.001/, Y2/10*0.0/
TB = 2000.
CALL FARG(P, TU, PHI, F, HU, U, V, S, Y1, CP, DLVLT, DLVLP)
DO 10 I = 1,MAXITS
CALL ECP(P, TB, PHI, HB, U, V, S, Y2, CP, DLVLT, DLVLP, IER)
DELT = (HU - HB)/CP
TB = TB + DELT
IF(ABS(DELT/TB) .LT. TOL) RETURN
10 CONTINUE
WRITE(6,20)
20 FORMAT('CONVERGENCE FAILURE IN TINITL')
RETURN
END
SUBROUTINE ECP (P, T, PHI, H, U, V, S, Y, CP, DLVLT, DLVLP, IER)
C
C PURPOSE:
C COMPUTE THE PROPERTIES OF EQUILIBRIUM COMBUSTION
C PRODUCTS
C
C GIVEN:
C P - PRESSURE (BARS)
C T - TEMPERATURE )K)
C PHI - FUEL AIR EQUIVALENCE RATIO
C
C RETURNS:
C H - ENTHALPY (J/G)
C U - INTERNAL ENERGY (J/G)
C V - SPECIFIC VOLUME (CM**3/G)
C S - ENTROPY (J/G/K)
C Y - A TEN DIMENSIONAL COMPOSITION VECTOR OF
C MOLE FRACTIONS
C 1=CO2. 2=H2O, 3=N2, 4=O2, 5=CO, 6=H2
C 7=H, 8=O, 9=OH, 10=NO
C CP - SPECIFIC HEAT AT CONSTANT PRESSURE (J/G/K)
C DLVLT - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C TEMPERATURE AT CONSTANT PRESSURE
C DLVLP - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C PRESSURE AT CONSTANT TEMPERATURE
C IER - AN ERROR MESSAGE
C 0 - NORMAL
C 1 - CONVERGENCE FAILURE: COMPOSITION LOOP
C 2 - CALLED WITH TOO HIGH AN EQUIVALENCE RATIO; C(S)
C AND OTHER SPECIES WOULD FORM
C
REAL M, MW, K, KP, MT, MP
LOGICAL CHECK
DIMENSION M(10), K(6), C(6), Y(10), D(3), B(4), A(4,4), Y0(6),
& KP(5,6), DFDT(4), DCDT(6), DKDT(6), DCDP(6),
& DYDT(10), DYDP(10), A0(7,10), DFDP(4), CP0(10),
& H0(10), S0(10), WK(12)
C
C FUEL DATA
C
DATA ALPHA/7.0/, BETA/17.0/, GAMMA/0./, DELTA/0./
C
C SPECIFIC HEAT DATA FROM GORDON AND MCBRDE (1971)
C
DATA A0/
1 .44608041E+01, .30981719E-02, -.12392571E-05, .22741325E-09,
1 -.15525954E-13, -.48961442E+05, -.98635982E+00,
2 .27167633E+01, .29451374E-02, -.80224374E-06, .10226682E-09,
2 -.48472145E-14, -.29905826E+05, .66305671E+01,
3 .28963194E+01, .15154866E-02, -.57235277E-06, .99807393E-10,
3 -.65223555E-14, -.90586184E+03, .61615148E+01,
4 .36219535E+01, .72518264E-03, -.19652228E-06, .36201558E-10,
4 -.28945627E-14, -.12019825E+04, .36150960E+01,
5 .29840696E+01, .14891390E-02, -.57899684E-06, .36201558E-10,
5 -.69353550E-14, -.14245228E+05, .63479156E+01,
6 .31001901E+01, .51119464E-03, .52644210E-07, -.34909973E-10,
6 .36945345E-14, -.87738042E+03, -.19629421E+01,
7 .25000000E+01, 0., 0., 0., 0., .25471627E+05, -.46011763E+00,
8 .55420596E+01, -.27550619E-04, -.31028033E-08, .45510674E-11,
8 -.43680515E-15, .29230803E+03, .49203080E+01,
9 .29106427E+01, .95931650E-03, -.19441702E-06, .13756646E-10,
9 .14224542E-15, .39353815E+04, .54423445E+01,
# .31890000E+01, .13382281E-02, -.52899318E-06, .95919332E-10,
# -.64847932E-14, .98283290E+04, .67458126E+01/
C
C EQUILIBRIUM CONSTANT DATA FROM OLIKARA AND BORMAN (1975)
C
DATA KP/
1 .432168E+00, -.112464E+05, .267269E+01, -.745744E-04,
1 .242484E-08,
2 .310805E+00, -.129540E+05, .321779E+01, -.738336E-04,
2 .344645E-08,
3 -.141784E+00, -.213308E+04, .853461E+00, .355015E-04,
3 -.310227E-08,
4 .150879E-01, -.470959E+04, .646096E+00, .272805E-05,
4 -.154444E-08,
5 -.752364E+00, .124210E+05, -.260284E+01, .259556E-03,
5 -.162687E-07,
6 -.415302E-02, .148627E+05, -.475746E+01, .124699E-03,
6 -.900227E-08/
C
C MOLECULAR WEIGHTS
C
DATA M/44.01, 18.02, 28.008, 32., 28.01, 2.018, 1.009, 16.,
& 17.009, 30.004/
C
C OTHER DATA
C
DATA MAXITS/100/, TOL/3.0E-05/, RU/8.31434/
C
C MAKE SURE SOLID CARBON WILL NOT FORM
C
IER = 2
EPS = .210/(ALPHA + 0.25*BETA - 0.5*GAMMA)
IF(PHI .GT. (.210/EPS/(0.5*ALPHA - 0.5*GAMMA)) ) RETURN
IER = 0
C
C DECIDE IF AN INITIAL ESTIMATE TO THE COMPOSITION IS
C REQUIRED OR IF T<1000 IN WHICH CASE FARG WILL SUFFICE
C
C
SUM = 0
DO 10 I = 1,10
10 SUM = SUM + Y(I)
IF( T .GT. 1000. .AND. SUM .GT. 0.998) GO TO 20
CALL FARG(P,T,PHI,1.,H,U,V,S,Y0,CP,DLVLT,DLVLP)
DO 30 I = 1,10
30 Y(I) = 0
DO 40 I = 1,6
40 Y(I) = Y0(I)
IF( T .LE. 1000. ) RETURN
20 CONTINUE
C
C
C EVALUATE THE REQUISITE CONSTANTS
C
PATM = .9869233*P
DO 50 I = 1,6
50 K(I) = 10.**(KP(1,I)*ALOG(T/1000.) + KP(2,I)/T + KP(3,I) +
& KP(4,I)*T + KP(5,I)*T*T)
C(1) = K(1)/SQRT(PATM)
C(2) = K(2)/SQRT(PATM)
C(3) = K(3)
C(4) = K(4)
C(5) = K(5)*SQRT(PATM)
C(6) = K(6)*SQRT(PATM)
D(1) = BETA/ALPHA
D(2) = (GAMMA + .42/EPS/PHI)/ALPHA
D(3) = (DELTA + 1.58/EPS/PHI)/ALPHA
C
C
C SET UP ITERATION LOOP FOR NEWTON-RAPHSON ITERATION AND
C INTRODUCE TRICK TO PREVENT INSTABILITIES NEAR PHI = 1.0
C AT LOW TEMPERATURES ON 24 BIT MACHINES (VAX). SET
C TRICK = 1.0 ON 48 BIT MACHINES (CDC)
C
TRICK = 10.
CHECK = ABS(PHI - 1.0) .LT. TOL
IF(CHECK) PHI = PHI*(1.0 + SIGN(TOL,PHI - 1.0))
ICHECK = 0
DO 5 J = 1,MAXITS
DO 4 JJ = 1,10
4 Y(JJ) = AMIN1(1.0,AMAX1(Y(JJ),1.0E-25))
D76 = 0.5*C(1)/SQRT(Y(6))
D84 = 0.5*C(1)/SQRT(Y(4))
D94 = 0.5*C(3)*SQRT(Y(6)/Y(4))
D96 = 0.5*C(3)*SQRT(Y(4)/Y(6))
D103 = 0.5*C(4)*SQRT(Y(4)/Y(3))
D104 = 0.5*C(4)*SQRT(Y(3)/Y(4))
D24 = 0.5*C(5)*Y(6)/SQRT(Y(4))
D26 = C(5)*SQRT(Y(4))
D14 = 0.5*C(6)*Y(5)/SQRT(Y(4))
D15 = C(6)*SQRT(Y(4))
A(1,1) = 1. + D103
A(1,2) = D14 + D24 + 1. +D84 + D104 + D96
A(1,3) = D15 + 1.0
A(1,4) = D26 + 1.0 + D76 +D96
A(2,1) = 0.0
A(2,2) = 2.0*D24 + D94 - D(1)*D14
A(2,3) = -D(1)*D15 - D(1)
A(2,4) = 2.0*D26 + 2.0 + D76 + D96
A(3,1) = D103
A(3,2) = 2.0*D14 + D24 + 2.0 + D84 + D04 + D104 - D(2)*D14
A(3,3) = 2.0*D15 + 1.0 - D(2)*D15 - D(2)
A(3,4) = D26 + D96
A(4,1) = 2.0 +D103
A(4,2) = D104 - D(3)*D14
A(4,3) = -D(3)*D15 - D(3)
A(4,4) = 0
C
C SOLVE THE MATRIX EQUATIONS 3.81 FOR COMPOSITION CORRECTIONS
C
SUM = 0
DO 55 I = 1,10
55 SUM = SUM + Y(I)
B(1) = -(SUM - 1.0)
B(2) = -(2.0*Y(2) + 2.0*Y(6) + Y(7) + Y(9) - D(1)*Y(1) -
& D(1)*Y(5))
B(3) = -(2.0*Y(1) + Y(2) + 2.0*Y(4) + Y(5) + Y(8) + Y(9)
& + Y(10) - D(2)*Y(1) - D(2)*Y(5))
B(4) = -(2.0*Y(3) + Y(10) - D(3)*Y(1) - D(3)*Y(5))
CALL LEQIF(A, 4, 4, 4, B, 4, 1, 0, WK, IERROR)
ERROR = 0.
DO 56 L = 3,6
LL = L - 2
Y(L) = Y(L) + B(LL)/TRICK
ERROR = AMAX1(ERROR,ABS(B(LL)))
Y(L) = AMIN(1.0,AMAX(Y(L),1.0E-25))
56 CONTINUE
Y(7) = C(1)*SQRT(Y(6))
Y(8) = C(2)*SQRT(Y(4))
Y(9) = C(3)*SQRT(Y(4)*Y(6))
Y(10) = C(4)*SQRT(Y(4)*Y(3))
Y(2) = C(5)*SQRT(Y(4))*Y(6)
Y(1) = C(6)*SQRT(Y(4))*Y(5)
IF(ERROR .LT. TOL) ICHECK = ICHECK +1
IF(ERROR .LT. TOL .AND. ICHECK .GE. 2) GO TO 57
5 CONTINUE
IER = 1
C
C COMPUTE THE REQUISITE CONSTANTS TO FIND PARTIAL DERIVATIVES
C
57 DO 60 I = 1,6
60 DKDT(I) = 2.302585*K(I)*(KP(1,I)/T - KP(2,I)/T/T + KP(4,I)
& + 2.0*KP(5,I)*T)
DCDT(1) = DKDT(1)/SQRT(PATM)
DCDT(2) = DKDT(2)/SQRT(PATM)
DCDT(3) = DKDT(3)
DCDT(4) = DKDT(4)
DCDT(5) = DKDT(5)*SQRT(PATM)
DCDT(6) = DKDT(6)*SQRT(PATM)
DCDP(1) = -0.5*C(1)/P
DCDP(2) = -0.5*C(2)/P
DCDP(5) = 0.5*C(5)/P
DCDP(6) = 0.5*C(6)/P
C
X1 = Y(1)/C(6)
X2 = Y(2)/C(5)
X7 = Y(7)/C(1)
X8 = Y(8)/C(2)
X9 = Y(9)/C(3)
X10 = Y(10)/C(4)
DFDT(1) = DCDT(6)*X1 + DCDT(5)*X2 + DCDT(1)*X7 +DCDT(2)*X8
& + DCDT(3)*X9 + DCDT(4)*X10
DFDT(2) = 2.0*DCDT(5)*X2 + DCDT(1)*X7 + DCDT(3)*X9 -
& D(1)*DCDT(6)*X1
DFDT(3) = 2.0*DCDT(6)*X1 + DCDT(5)*X2 + DCDT(2)*X8
& + DCDT(3)*X9 + DCDT(4)*X10 - D(2)*DCDT(6)*X1
DFDT(4) = DCDT(4)*X10 - D(3)*DCDT(6)*X1
DFDP(1) = DCDP(6)*X1 + DCDP(5)*X2 + DCDP(1)*X7 + DCDP(2)*X8
DFDP(2) = 2.0*DCDP(5)*X2 + DCDP(1)*X7 - D(1)*DCDP(6)*X1
DFDP(3) = 2.0*DCDP(6)*X1 + DCDP(5)*X2 + DCDP(2)*X8 -
& D(2)*DCDP(6)*X1
DFDP(4) = -D(3)*DCDP(6)*X1
C
C SOLVE THE MATRIX EQUATIONS 3.91 FOR THE INDEPENDENT DERIVATIVES
C AND THEN DETERMINE THE DEPENDENT DERIVATICES
C
DO 70 I = 1,4
70 B(I) = -DFDT(I)
CALL LEQIF(A, 4, 4, 4, B, 4, 1, 1, WK, IERROR)
DYDT(3) = B(1)
DYDT(4) = B(2)
DYDT(5) = B(3)
DYDT(6) = B(4)
DYDT(1) = SQRT(Y(4))*Y(5)*DCDT(6) + D14*DYDT(4) + D15*DYDT(5)
DYDT(2) = SQRT(Y(4))*Y(6)*DCDT(5) + D24*DYDT(4) + D26*DYDT(6)
DYDT(7) = SQRT(Y(6))*DCDT(1) + D76*DYDT(6)
DYDT(8) = SQRT(Y(4))*DCDT(2) + D84*DYDT(4)
DYDT(9) = SQRT(Y(4)*Y(6))*DCDT(3) + D94*DYDT(4) + D96*DYDT(6)
DYDT(10) = SQRT(Y(4)*Y(3))*DCDT(4) + D104*DYDT(4) + D103*DYDT(3)
C
DO 80 I = 1,4
80 B(I) = -DFDP(I)
CALL LEQIF(A, 4, 4, 4, B, 4, 1, 1, WK, IERROR)
DYDP(3) = B(1)
DYDP(4) = B(2)
DYDP(5) = B(3)
DYDP(6) = B(4)
DYDP(1) = SQRT(Y(4))*Y(5)*DCDP(6) + D14*DYDP(4) + D15*DYDP(5)
DYDP(2) = SQRT(Y(4))*Y(6)*DCDP(5) + D24*DYDP(4) + D26*DYDP(6)
DYDP(7) = SQRT(Y(6))*DCDP(1) + D76*DYDP(6)
DYDP(8) = SQRT(Y(4))*DCDP(2) + D84*DYDP(4)
DYDP(9) = D94*DYDP(4) + D96*DYDP(6)
DYDP(10) = D104*DYDP(4) + D103*DYDP(3)
C
DO 90 I = 1,10
CP0(I) = A0(1,I) + A0(2,I)*T + A0(3,I)*T*T + A0(4,I)*T**3
& + A0(5,I)*T**4
H0(I) = A0(1,I) + A0(2,I)/2.*T + A0(3,I)/3.*T*T + A0(4,I)/4.*
& T**3 + A0(5,1)/5.*T**4 + A0(6,I)/T
90 S0(I) = A0(1,I)*ALOG(T) + A0(2,I)*T + A0(3,I)/2.*T**2
& + A0(4,I)/3.*T**3 + A0(5,I)/4.*T**4 +A0(7,I)
C
C Y(1), Y(2) ARE REEVALUATED TO CLEAN UP ANY ROUND-OFF ERRORS WHICH
C CAN OCCUR FOR THE LOW TEMPERATURE STOICHIOMETRIC SOLUTIONS
C
Y(1) = (2.*Y(3) + Y(10))/D(3) - Y(5)
Y(2) = (D(1)/D(3)*(2.*Y(3) + Y(10)) - 2.*Y(6) - Y(7) - Y(9))/2
C
MW = 0.
CP = 0.
H = 0.
S = 0.
MT = 0.
MP = 0.
DO 100 I = 1,10
IF( Y(I) .LE. 1.0E-25 ) GO TO 100
H = H + H0(I)*Y(I)
MW = MW + M(I)*Y(I)
MT = MT + M(I)*DYDT(I)
MP = MP + M(I)*DYDP(I)
CP = CP + Y(I)*CP0(I) + H0(I)*T*DYDT(I)
S = S + Y(I)*(S0(I) - ALOG(Y(I)))
100 CONTINUE
C
R = RU/MW
V = 10.*R*T/P
CP = R*(CP - H*T*MT/MW)
DLVLT = 1.0 + AMAX1(-T*MT/MW,0.)
DLVLP = -1.0 - AMAX1(P*MP/MW,0.)
H = H*R*T
S = R*(-ALOG(PATM) + 5)
U = H - R*T
RETURN
END
C
C
SUBROUTINE FARG (P, T, PHI, F, H, U, V, S, Y, CP, DLVLT, DLVLP)
C
C PURPOSE:
C COMPUTE THE PROPERTIES OF A FUEL, AIR, RESIDUAL
C GAS MIXTURE
C
C GIVEN:
C P - PRESSURE (BARS)
C T - TEMPERATURE (K)
C PHI - FUEL AIR EQUIVALENCE RATIO
C F - RESIDUAL MASS FRACTION
C
C RETURNS:
C H - ENTHALPY (J/G)
C U - INTERNAL ENERGY (J/G)
C V - SPECIFIC VOLUME (CM**3/g)
C S - ENTROPY (J/G/K)
C Y - A SIX DIMENSIONAL COMPOSITION VECTOR OF
C MOLE FRACTIONS
C 1=CO2, 2=H2O, 3=N2, 4O2, 5=CO, 6=H2
C CP - SPECIFIC HEAT AT CONSTANT PRESSURE (J/G/K)
C DLVLT - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C TEMPERATURE AT CONSTANT PRESSURE
C DLVLP - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C PRESSURE AT CONSTANT TEMPERATURE
C
C REMARKS:
C 1. VALID FOR 300 < T < 1000
C 2. ASSUMES THE FUEL IS GASOLINE
C 3. ENTHALPIES OF O2, H2, N2 AND C(S) ARE SET TO ZERO
C AT 298K
C 4. THE FUEL MOLE FRACTION IS UNITY MINUS THE SUM OF Y(I)
C
REAL M, MW, MRES, MFA, MFUEL, K, NU, N2
LOGICAL RICH, LEAN
DIMENSION A(7,6), TABLE(6), M(6), NU(6), Y(6), CP0(6), H0(6),
& S0(6)
C
C FUEL DATA
C
DATA ALPHA/7.0/, BETA/17.0/, GAMMA/0./, DELTA/0./,
& A0/4.0652/, B0/6.0977E-02/, C0/-1.8801E-05/,
& D0/-3.5880E+04/, E0/15.45/
C
C TABLE 3.1, I = 1,6
C
DATA A/
1 .24007797E+01, .87350957E-02, -.6607878E-05, .20021861E-08,
1 .63274039E-15, -.48377527E+05, .96951457E+01,
2 .40701275E+01, -.11084499E-02, .41521180E-05, -.29637404E-08,
2 .80702103E-12, -.30279722E+05, -.32270046E+00,
3 .36748261E+01, -.12081500E-02, .23240102E-05, -.63217559E-09,
3 -.22577253E-12, -.10611588E+04, .23580424E+01,
4 .36255985E+01, -.18782184E-02, .70554544E-05, -.67635137E-08,
4 .21555993E-11, -.10475226E+04, .43052778E+01,
5 .37100928E+01, -.16190964E-02, .36923594E-05, -.20319674E-08,
5 .23953344E-12, -.14356310E+05, .29555355E+01,
6 .30574451E+01, .26765200E-02, -.58099162E-05, .55210391E-08,
6 -.18122739E-11, -.98890474E+03, -.22997056E+01/
C
C OTHER DATA
C
DATA RU/8.315/, TABLE/-1.,1.,0.,0.,1.,-1./, M/44.01, 18.02,
& 28.008, 32.000, 28.01, 2.018/
C
C COMPUTE RESIDUAL GAS COMPOSITION ACCORDING TO TABLE 3.3
C
RICH = PHI .GT. 1.0
LEAN = .NOT. RICH
DLVLT = 1.0
DLVLP = -1.0
EPS = .21/(ALPHA + 0.25*BETA - 0.5*GAMMA)
IF (RICH) GO TO 10
NU(1) = ALPHA*PHI*EPS
NU(2) = BETA*PHI*EPS/2
NU(3) = 0.79 + DELTA*PHI*EPS/2
NU(4) = 0.21*(1.0 - PHI)
NU(5) = 0.
NU(6) = 0.
DCDT = 0.
GO TO 20
10 Z = 1000./T
K = EXP(2.743 + Z*(-1.761 + Z*(-1.611 + Z*.2803)))
DKDT = -K*(-1.761 + Z*(-3.222 + Z*.8409))/1000.
A1 = 1.0 - K
B = .42 - PHI*EPS*(2.*ALPHA - GAMMA) + K*(.42*(PHI - 1.) +
& ALPHA*PHI*EPS)
C = -.42*ALPHA*PHI*EPS*(PHI - 1.0)*K
NU(5) = (-B + SQRT(B*B - 4.*A1*C))/2./A1
DCDT = DKDT*(NU(5)**2 - NU(5)*(.42*(PHI - 1.) + ALPHA*
& PHI*EPS) + .42*ALPHA*PHI*EPS*(PHI - 1.))/
& (2.*NU(5)*A1 + B)
NU(1) = ALPHA*PHI*EPS - NU(5)
NU(2) = .42 - PHI*EPS*(2.*ALPHA - GAMMA) + NU(5)
NU(3) = .79 + DELTA*PHI*EPS/2.
NU(4) = 0.
NU(6) = .42*(PHI - 1.) - NU(5)
C
C COMPUTE MOLE FRACTIONS AND MOLECULAR WEIGHT OF RESIDUAL
C
20 TMOLES = 0.
DO 30 I = 1,6
30 TMOLES = TMOLES + NU(I)
MRES = 0.
DO 40 I = 1,6
Y(I) = NU(I)/TMOLES
40 MRES = MRES + Y(I)*M(I)
C
C COMPUTE MOLE FRACTIONS AND MOLECULAR WEIGHT OF FURLE-AIR
C
FUEL = EPS*PHI/(1. + EPS*PHI)
O2 = .21/(1. + EPS*PHI)
N2 = .79/(1. + EPS*PHI)
MFA = FUEL*(12.01*ALPHA + 1.008*BETA + 16.*GAMMA + 14.01*
& DELTA) + 32.*O2 + 28.02*N2
C
C COMPUTE MOLE FRACTIONS OF FUEL-AIR-RESIDUAL GAS
C
YRES = F/(F + MRES/MFA*(1. - F))
DO 50 I = 1,6
50 Y(I) = Y(I)*YRES
YFUEL = FUEL*(1. - YRES)
Y(3) = Y(3) + N2*(1. - YRES)
Y(4) = Y(4) + O2*(1. - YRES)
C
C COMPUTE COMPONENT PROPERTIES
C
DO 60 I = 1,6
CP0(I) = A(1,I) + A(2,I)*T + A(3,I)*T**2 + A(4,I)*T**3 +
& A(5,I)*T**4
H0(I) = A(1,I) + A(2,I)/2.*T + A(3,I)/3.*T**2 + A(4,I)/4.*
& T**3 + A(5,I)/5.*T**4 + A(6,I)/T
60 S0(I) = A(1,I)*ALOG(T) + A(2,I)*T + A(3,I)/2.*T**2 +
& A(4,I)/3.*T**3 + A(5,I)/4.*T**4 + A(7,1)
C
MFUEL = 12.01*ALPHA + 1.008*BETA + 16.000*GAMMA + 14.01*DELTA
CPFUEL = A0 + B0*T + C0*T**2
HFUEL = A0 + B0/2.*T + C0/3.*T**2 + D0/T
S0FUEL = A0*ALOG(T) + B0*T + C0/2.*T**2 + E0
C
C COMPUTE PROPERTIES OF MIXTURE
C
H = HFUEL*YFUEL
S = (S0FUEL - ALOG(YFUEL))*YFUEL
CP = CPFUEL*YFUEL
MW = MFUEL*YFUEL
DO 70 I = 1,6
H = H + H0(I)*Y(I)
S = S + Y(I)*(S0(I) - ALOG(Y(I)))
CP = CP + CP0(I)*Y(I) + H0(I)*T*TABLE(I)*DCDT*YRES/TMOLES
70 MW = MW + Y(I)*M(I)
C
R = RU/MW
H = R*T*H
U = H - R*T
V = 10.*R*T/P
S = R*(-ALOG(.9869*P) + S)
CP = R*CP
C
RETURN
END
C
C
SUBROUTINE DVERK(N, FCN, X, Y, XEND, TOL, IND, C, NW, W)
INTEGER N, IND, NW, K
DOUBLE PRECISION X, Y(N), XEND, TOL, C(*), W(NW,9), TEMP
C Added external statement for fcn to avoid a warning message.
external fcn
C
C***********************************************************************
C *
C PURPOSE - THIS IS A RUNGE-KUTTA SUBROUTINE BASED ON VERNER'S *
C FIFTH AND SIXTH ORDER PAIR OF FORMULAS FOR FINDING APPROXIMATIONS TO *
C THE SOLUTION OF A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL *
C EQUATIONS WITH INITIAL CONDITIONS. IT ATTEMPTS TO KEEP THE GLOBAL *
C ERROR PROPORTIONAL TO A TOLERANCE SPECIFIED BY THE USER. (THE *
C PROPORTIONALITY DEPENDS ON THE KIND OF ERROR CONTROL THAT IS USED, *
C AS WELL AS THE DIFFERENTIAL EQUATION AND THE RANGE OF INTEGRATION.) *
C *
C VARIOUS OPTIONS ARE AVAILABLE TO THE USER, INCLUDING DIFFERENT *
C KINDS OF ERROR CONTROL, RESTRICTIONS ON STEP SIZES, AND INTERRUPTS *
C WHICH PERMIT THE USER TO EXAMINE THE STATE OF THE CALCULATION (AND *
C PERHAPS MAKE MODIFICATIONS) DURING INTERMEDIATE STAGES. *
C *
C THE PROGRAM IS EFFICIENT FOR NON-STIFF SYSTEMS. HOWEVER, A GOOD *
C VARIABLE-ORDER-ADAMS METHOD WILL PROBABLY BE MORE EFFICIENT IF THE *
C FUNCTION EVALUATIONS ARE VERY COSTLY. SUCH A METHOD WOULD ALSO BE *
C MORE SUITABLE IF ONE WANTED TO OBTAIN A LARGE NUMBER OF INTERMEDIATE *
C SOLUTION VALUES BY INTERPOLATION, AS MIGHT BE THE CASE FOR EXAMPLE *
C WITH GRAPHICAL OUTPUT. *
C *
C HULL-ENRIGHT-JACKSON 1/10/76 *
C *
C***********************************************************************
C *
C USE - THE USER MUST SPECIFY EACH OF THE FOLLOWING *
C *
C N NUMBER OF EQUATIONS *
C *
C FCN NAME OF SUBROUTINE FOR EVALUATING FUNCTIONS - THE SUBROUTINE *
C ITSELF MUST ALSO BE PROVIDED BY THE USER - IT SHOULD BE OF *
C THE FOLLOWING FORM *
C SUBROUTINE FCN(N, X, Y, YPRIME) *
C INTEGER N *
C DOUBLE PRECISION X, Y(N), YPRIME(N) *
C *** ETC *** *
C AND IT SHOULD EVALUATE YPRIME, GIVEN N, X AND Y *
C *
C X INDEPENDENT VARIABLE - INITIAL VALUE SUPPLIED BY USER *
C *
C Y DEPENDENT VARIABLE - INITIAL VALUES OF COMPONENTS Y(1), Y(2), *
C ..., Y(N) SUPPLIED BY USER *
C *
C XEND VALUE OF X TO WHICH INTEGRATION IS TO BE CARRIED OUT - IT MAY *
C BE LESS THAN THE INITIAL VALUE OF X *
C *
C TOL TOLERANCE - THE SUBROUTINE ATTEMPTS TO CONTROL A NORM OF THE *
C LOCAL ERROR IN SUCH A WAY THAT THE GLOBAL ERROR IS *
C PROPORTIONAL TO TOL. IN SOME PROBLEMS THERE WILL BE ENOUGH *
C DAMPING OF ERRORS, AS WELL AS SOME CANCELLATION, SO THAT *
C THE GLOBAL ERROR WILL BE LESS THAN TOL. ALTERNATIVELY, THE *
C CONTROL CAN BE VIEWED AS ATTEMPTING TO PROVIDE A *
C CALCULATED VALUE OF Y AT XEND WHICH IS THE EXACT SOLUTION *
C TO THE PROBLEM Y' = F(X,Y) + E(X) WHERE THE NORM OF E(X) *
C IS PROPORTIONAL TO TOL. (THE NORM IS A MAX NORM WITH *
C WEIGHTS THAT DEPEND ON THE ERROR CONTROL STRATEGY CHOSEN *
C BY THE USER. THE DEFAULT WEIGHT FOR THE K-TH COMPONENT IS *
C 1/MAX(1,ABS(Y(K))), WHICH THEREFORE PROVIDES A MIXTURE OF *
C ABSOLUTE AND RELATIVE ERROR CONTROL.) *
C *
C IND INDICATOR - ON INITIAL ENTRY IND MUST BE SET EQUAL TO EITHER *
C 1 OR 2. IF THE USER DOES NOT WISH TO USE ANY OPTIONS, HE *
C SHOULD SET IND TO 1 - ALL THAT REMAINS FOR THE USER TO DO *
C THEN IS TO DECLARE C AND W, AND TO SPECIFY NW. THE USER *
C MAY ALSO SELECT VARIOUS OPTIONS ON INITIAL ENTRY BY *
C SETTING IND = 2 AND INITIALIZING THE FIRST 9 COMPONENTS OF *
C C AS DESCRIBED IN THE NEXT SECTION. HE MAY ALSO RE-ENTER *
C THE SUBROUTINE WITH IND = 3 AS MENTIONED AGAIN BELOW. IN *
C ANY EVENT, THE SUBROUTINE RETURNS WITH IND EQUAL TO *
C 3 AFTER A NORMAL RETURN *
C 4, 5, OR 6 AFTER AN INTERRUPT (SEE OPTIONS C(8), C(9)) *
C -1, -2, OR -3 AFTER AN ERROR CONDITION (SEE BELOW) *
C *
C C COMMUNICATIONS VECTOR - THE DIMENSION MUST BE GREATER THAN OR *
C EQUAL TO 24, UNLESS OPTION C(1) = 4 OR 5 IS USED, IN WHICH *
C CASE THE DIMENSION MUST BE GREATER THAN OR EQUAL TO N+30 *
C *
C NW FIRST DIMENSION OF WORKSPACE W - MUST BE GREATER THAN OR *
C EQUAL TO N *
C *
C W WORKSPACE MATRIX - FIRST DIMENSION MUST BE NW AND SECOND MUST *
C BE GREATER THAN OR EQUAL TO 9 *
C *
C THE SUBROUTINE WILL NORMALLY RETURN WITH IND = 3, HAVING *
C REPLACED THE INITIAL VALUES OF X AND Y WITH, RESPECTIVELY, THE VALUE *
C OF XEND AND AN APPROXIMATION TO Y AT XEND. THE SUBROUTINE CAN BE *
C CALLED REPEATEDLY WITH NEW VALUES OF XEND WITHOUT HAVING TO CHANGE *
C ANY OTHER ARGUMENT. HOWEVER, CHANGES IN TOL, OR ANY OF THE OPTIONS *
C DESCRIBED BELOW, MAY ALSO BE MADE ON SUCH A RE-ENTRY IF DESIRED. *
C *
C THREE ERROR RETURNS ARE ALSO POSSIBLE, IN WHICH CASE X AND Y *
C WILL BE THE MOST RECENTLY ACCEPTED VALUES - *
C WITH IND = -3 THE SUBROUTINE WAS UNABLE TO SATISFY THE ERROR *
C REQUIREMENT WITH A PARTICULAR STEP-SIZE THAT IS LESS THAN OR *
C EQUAL TO HMIN, WHICH MAY MEAN THAT TOL IS TOO SMALL *
C WITH IND = -2 THE VALUE OF HMIN IS GREATER THAN HMAX, WHICH *
C PROBABLY MEANS THAT THE REQUESTED TOL (WHICH IS USED IN THE *
C CALCULATION OF HMIN) IS TOO SMALL *
C WITH IND = -1 THE ALLOWED MAXIMUM NUMBER OF FCN EVALUATIONS HAS *
C BEEN EXCEEDED, BUT THIS CAN ONLY OCCUR IF OPTION C(7), AS *
C DESCRIBED IN THE NEXT SECTION, HAS BEEN USED *
C *
C THERE ARE SEVERAL CIRCUMSTANCES THAT WILL CAUSE THE CALCULATIONS *
C TO BE TERMINATED, ALONG WITH OUTPUT OF INFORMATION THAT WILL HELP *
C THE USER DETERMINE THE CAUSE OF THE TROUBLE. THESE CIRCUMSTANCES *
C INVOLVE ENTRY WITH ILLEGAL OR INCONSISTENT VALUES OF THE ARGUMENTS, *
C SUCH AS ATTEMPTING A NORMAL RE-ENTRY WITHOUT FIRST CHANGING THE *
C VALUE OF XEND, OR ATTEMPTING TO RE-ENTER WITH IND LESS THAN ZERO. *
C *
C***********************************************************************
C *
C OPTIONS - IF THE SUBROUTINE IS ENTERED WITH IND = 1, THE FIRST 9 *
C COMPONENTS OF THE COMMUNICATIONS VECTOR ARE INITIALIZED TO ZERO, AND *
C THE SUBROUTINE USES ONLY DEFAULT VALUES FOR EACH OPTION. IF THE *
C SUBROUTINE IS ENTERED WITH IND = 2, THE USER MUST SPECIFY EACH OF *
C THESE 9 COMPONENTS - NORMALLY HE WOULD FIRST SET THEM ALL TO ZERO, *
C AND THEN MAKE NON-ZERO THOSE THAT CORRESPOND TO THE PARTICULAR *
C OPTIONS HE WISHES TO SELECT. IN ANY EVENT, OPTIONS MAY BE CHANGED ON *
C RE-ENTRY TO THE SUBROUTINE - BUT IF THE USER CHANGES ANY OF THE *
C OPTIONS, OR TOL, IN THE COURSE OF A CALCULATION HE SHOULD BE CAREFUL *
C ABOUT HOW SUCH CHANGES AFFECT THE SUBROUTINE - IT MAY BE BETTER TO *
C RESTART WITH IND = 1 OR 2. (COMPONENTS 10 TO 24 OF C ARE USED BY THE *
C PROGRAM - THE INFORMATION IS AVAILABLE TO THE USER, BUT SHOULD NOT *
C NORMALLY BE CHANGED BY HIM.) *
C *
C C(1) ERROR CONTROL INDICATOR - THE NORM OF THE LOCAL ERROR IS THE *
C MAX NORM OF THE WEIGHTED ERROR ESTIMATE VECTOR, THE *
C WEIGHTS BEING DETERMINED ACCORDING TO THE VALUE OF C(1) - *
C IF C(1)=1 THE WEIGHTS ARE 1 (ABSOLUTE ERROR CONTROL) *
C IF C(1)=2 THE WEIGHTS ARE 1/ABS(Y(K)) (RELATIVE ERROR *
C CONTROL) *
C IF C(1)=3 THE WEIGHTS ARE 1/MAX(ABS(C(2)),ABS(Y(K))) *
C (RELATIVE ERROR CONTROL, UNLESS ABS(Y(K)) IS LESS *
C THAN THE FLOOR VALUE, ABS(C(2)) ) *
C IF C(1)=4 THE WEIGHTS ARE 1/MAX(ABS(C(K+30)),ABS(Y(K))) *
C (HERE INDIVIDUAL FLOOR VALUES ARE USED) *
C IF C(1)=5 THE WEIGHTS ARE 1/ABS(C(K+30)) *
C FOR ALL OTHER VALUES OF C(1), INCLUDING C(1) = 0, THE *
C DEFAULT VALUES OF THE WEIGHTS ARE TAKEN TO BE *
C 1/MAX(1,ABS(Y(K))), AS MENTIONED EARLIER *
C (IN THE TWO CASES C(1) = 4 OR 5 THE USER MUST DECLARE THE *
C DIMENSION OF C TO BE AT LEAST N+30 AND MUST INITIALIZE THE *
C COMPONENTS C(31), C(32), ..., C(N+30).) *
C *
C C(2) FLOOR VALUE - USED WHEN THE INDICATOR C(1) HAS THE VALUE 3 *
C *
C C(3) HMIN SPECIFICATION - IF NOT ZERO, THE SUBROUTINE CHOOSES HMIN *
C TO BE ABS(C(3)) - OTHERWISE IT USES THE DEFAULT VALUE *
C 10*MAX(DWARF,RREB*MAX(WEIGHTED NORM Y/TOL,ABS(X))), *
C WHERE DWARF IS A VERY SMALL POSITIVE MACHINE NUMBER AND *
C RREB IS THE RELATIVE ROUNDOFF ERROR BOUND *
C *
C C(4) HSTART SPECIFICATION - IF NOT ZERO, THE SUBROUTINE WILL USE *
C AN INITIAL HMAG EQUAL TO ABS(C(4)), EXCEPT OF COURSE FOR *
C THE RESTRICTIONS IMPOSED BY HMIN AND HMAX - OTHERWISE IT *
C USES THE DEFAULT VALUE OF HMAX*(TOL)**(1/6) *
C *
C C(5) SCALE SPECIFICATION - THIS IS INTENDED TO BE A MEASURE OF THE *
C SCALE OF THE PROBLEM - LARGER VALUES OF SCALE TEND TO MAKE *
C THE METHOD MORE RELIABLE, FIRST BY POSSIBLY RESTRICTING *
C HMAX (AS DESCRIBED BELOW) AND SECOND, BY TIGHTENING THE *
C ACCEPTANCE REQUIREMENT - IF C(5) IS ZERO, A DEFAULT VALUE *
C OF 1 IS USED. FOR LINEAR HOMOGENEOUS PROBLEMS WITH *
C CONSTANT COEFFICIENTS, AN APPROPRIATE VALUE FOR SCALE IS A *
C NORM OF THE ASSOCIATED MATRIX. FOR OTHER PROBLEMS, AN *
C APPROXIMATION TO AN AVERAGE VALUE OF A NORM OF THE *
C JACOBIAN ALONG THE TRAJECTORY MAY BE APPROPRIATE *
C *
C C(6) HMAX SPECIFICATION - FOUR CASES ARE POSSIBLE *
C IF C(6).NE.0 AND C(5).NE.0, HMAX IS TAKEN TO BE *
C MIN(ABS(C(6)),2/ABS(C(5))) *
C IF C(6).NE.0 AND C(5).EQ.0, HMAX IS TAKEN TO BE ABS(C(6)) *
C IF C(6).EQ.0 AND C(5).NE.0, HMAX IS TAKEN TO BE *
C 2/ABS(C(5)) *
C IF C(6).EQ.0 AND C(5).EQ.0, HMAX IS GIVEN A DEFAULT VALUE *
C OF 2 *
C *
C C(7) MAXIMUM NUMBER OF FUNCTION EVALUATIONS - IF NOT ZERO, AN *
C ERROR RETURN WITH IND = -1 WILL BE CAUSED WHEN THE NUMBER *
C OF FUNCTION EVALUATIONS EXCEEDS ABS(C(7)) *
C *
C C(8) INTERRUPT NUMBER 1 - IF NOT ZERO, THE SUBROUTINE WILL *
C INTERRUPT THE CALCULATIONS AFTER IT HAS CHOSEN ITS *
C PRELIMINARY VALUE OF HMAG, AND JUST BEFORE CHOOSING HTRIAL *
C AND XTRIAL IN PREPARATION FOR TAKING A STEP (HTRIAL MAY *
C DIFFER FROM HMAG IN SIGN, AND MAY REQUIRE ADJUSTMENT IF *
C XEND IS NEAR) - THE SUBROUTINE RETURNS WITH IND = 4, AND *
C WILL RESUME CALCULATION AT THE POINT OF INTERRUPTION IF *
C RE-ENTERED WITH IND = 4 *
C *
C C(9) INTERRUPT NUMBER 2 - IF NOT ZERO, THE SUBROUTINE WILL *
C INTERRUPT THE CALCULATIONS IMMEDIATELY AFTER IT HAS *
C DECIDED WHETHER OR NOT TO ACCEPT THE RESULT OF THE MOST *
C RECENT TRIAL STEP, WITH IND = 5 IF IT PLANS TO ACCEPT, OR *
C IND = 6 IF IT PLANS TO REJECT - Y(*) IS THE PREVIOUSLY *
C ACCEPTED RESULT, WHILE W(*,9) IS THE NEWLY COMPUTED TRIAL *
C VALUE, AND W(*,2) IS THE UNWEIGHTED ERROR ESTIMATE VECTOR. *
C THE SUBROUTINE WILL RESUME CALCULATIONS AT THE POINT OF *
C INTERRUPTION ON RE-ENTRY WITH IND = 5 OR 6. (THE USER MAY *
C CHANGE IND IN THIS CASE IF HE WISHES, FOR EXAMPLE TO FORCE *
C ACCEPTANCE OF A STEP THAT WOULD OTHERWISE BE REJECTED, OR *
C VICE VERSA. HE CAN ALSO RESTART WITH IND = 1 OR 2.) *
C *
C***********************************************************************
C *
C SUMMARY OF THE COMPONENTS OF THE COMMUNICATIONS VECTOR *
C *
C PRESCRIBED AT THE OPTION DETERMINED BY THE PROGRAM *
C OF THE USER *
C *
C C(10) RREB(REL ROUNDOFF ERR BND) *
C C(1) ERROR CONTROL INDICATOR C(11) DWARF (VERY SMALL MACH NO) *
C C(2) FLOOR VALUE C(12) WEIGHTED NORM Y *
C C(3) HMIN SPECIFICATION C(13) HMIN *
C C(4) HSTART SPECIFICATION C(14) HMAG *
C C(5) SCALE SPECIFICATION C(15) SCALE *
C C(6) HMAX SPECIFICATION C(16) HMAX *
C C(7) MAX NO OF FCN EVALS C(17) XTRIAL *
C C(8) INTERRUPT NO 1 C(18) HTRIAL *
C C(9) INTERRUPT NO 2 C(19) EST *
C C(20) PREVIOUS XEND *
C C(21) FLAG FOR XEND *
C C(22) NO OF SUCCESSFUL STEPS *
C C(23) NO OF SUCCESSIVE FAILURES *
C C(24) NO OF FCN EVALS *
C *
C IF C(1) = 4 OR 5, C(31), C(32), ... C(N+30) ARE FLOOR VALUES *
C *
C RREB and DWARF are machine dependent constants currently set so *
C that they should be appropriate for most machines. However, it may *
C be appropriate to change them when this program is installed on a *
C new machine. K.R.J. 3 Oct 1991. *
C *
C***********************************************************************
C *
C AN OVERVIEW OF THE PROGRAM *
C *
C BEGIN INITIALIZATION, PARAMETER CHECKING, INTERRUPT RE-ENTRIES *
C ......ABORT IF IND OUT OF RANGE 1 TO 6 *
C . CASES - INITIAL ENTRY, NORMAL RE-ENTRY, INTERRUPT RE-ENTRIES *
C . CASE 1 - INITIAL ENTRY (IND .EQ. 1 OR 2) *
C V........ABORT IF N.GT.NW OR TOL.LE.0 *
C . IF INITIAL ENTRY WITHOUT OPTIONS (IND .EQ. 1) *
C . SET C(1) TO C(9) EQUAL TO ZERO *
C . ELSE INITIAL ENTRY WITH OPTIONS (IND .EQ. 2) *
C . MAKE C(1) TO C(9) NON-NEGATIVE *
C . MAKE FLOOR VALUES NON-NEGATIVE IF THEY ARE TO BE USED *
C . END IF *
C . INITIALIZE RREB, DWARF, PREV XEND, FLAG, COUNTS *
C . CASE 2 - NORMAL RE-ENTRY (IND .EQ. 3) *
C .........ABORT IF XEND REACHED, AND EITHER X CHANGED OR XEND NOT *
C . RE-INITIALIZE FLAG *
C . CASE 3 - RE-ENTRY FOLLOWING AN INTERRUPT (IND .EQ. 4 TO 6) *
C V TRANSFER CONTROL TO THE APPROPRIATE RE-ENTRY POINT....... *
C . END CASES . *
C . END INITIALIZATION, ETC. . *
C . V *
C . LOOP THROUGH THE FOLLOWING 4 STAGES, ONCE FOR EACH TRIAL STEP . *
C . STAGE 1 - PREPARE . *
C***********ERROR RETURN (WITH IND=-1) IF NO OF FCN EVALS TOO GREAT . *
C . CALC SLOPE (ADDING 1 TO NO OF FCN EVALS) IF IND .NE. 6 . *
C . CALC HMIN, SCALE, HMAX . *
C***********ERROR RETURN (WITH IND=-2) IF HMIN .GT. HMAX . *
C . CALC PRELIMINARY HMAG . *
C***********INTERRUPT NO 1 (WITH IND=4) IF REQUESTED.......RE-ENTRY.V *
C . CALC HMAG, XTRIAL AND HTRIAL . *
C . END STAGE 1 . *
C V STAGE 2 - CALC YTRIAL (ADDING 7 TO NO OF FCN EVALS) . *
C . STAGE 3 - CALC THE ERROR ESTIMATE . *
C . STAGE 4 - MAKE DECISIONS . *
C . SET IND=5 IF STEP ACCEPTABLE, ELSE SET IND=6 . *
C***********INTERRUPT NO 2 IF REQUESTED....................RE-ENTRY.V *
C . IF STEP ACCEPTED (IND .EQ. 5) *
C . UPDATE X, Y FROM XTRIAL, YTRIAL *
C . ADD 1 TO NO OF SUCCESSFUL STEPS *
C . SET NO OF SUCCESSIVE FAILURES TO ZERO *
C**************RETURN(WITH IND=3, XEND SAVED, FLAG SET) IF X .EQ. XEND *
C . ELSE STEP NOT ACCEPTED (IND .EQ. 6) *
C . ADD 1 TO NO OF SUCCESSIVE FAILURES *
C**************ERROR RETURN (WITH IND=-3) IF HMAG .LE. HMIN *
C . END IF *
C . END STAGE 4 *
C . END LOOP *
C . *
C BEGIN ABORT ACTION *
C OUTPUT APPROPRIATE MESSAGE ABOUT STOPPING THE CALCULATIONS, *
C ALONG WITH VALUES OF IND, N, NW, TOL, HMIN, HMAX, X, XEND, *
C PREVIOUS XEND, NO OF SUCCESSFUL STEPS, NO OF SUCCESSIVE *
C FAILURES, NO OF FCN EVALS, AND THE COMPONENTS OF Y *
C STOP *
C END ABORT ACTION *
C *
C***********************************************************************
C
C ******************************************************************
C * BEGIN INITIALIZATION, PARAMETER CHECKING, INTERRUPT RE-ENTRIES *
C ******************************************************************
C
C ......ABORT IF IND OUT OF RANGE 1 TO 6
IF (IND.LT.1 .OR. IND.GT.6) GO TO 500
C
C CASES - INITIAL ENTRY, NORMAL RE-ENTRY, INTERRUPT RE-ENTRIES
GO TO (5, 5, 45, 1111, 2222, 2222), IND
C CASE 1 - INITIAL ENTRY (IND .EQ. 1 OR 2)
C .........ABORT IF N.GT.NW OR TOL.LE.0
5 IF (N.GT.NW .OR. TOL.LE.0.D0) GO TO 500
IF (IND.EQ. 2) GO TO 15
C INITIAL ENTRY WITHOUT OPTIONS (IND .EQ. 1)
C SET C(1) TO C(9) EQUAL TO 0
DO 10 K = 1, 9
C(K) = 0.D0
10 CONTINUE
GO TO 35
15 CONTINUE
C INITIAL ENTRY WITH OPTIONS (IND .EQ. 2)
C MAKE C(1) TO C(9) NON-NEGATIVE
DO 20 K = 1, 9
C(K) = DABS(C(K))
20 CONTINUE
C MAKE FLOOR VALUES NON-NEGATIVE IF THEY ARE TO BE USED
IF (C(1).NE.4.D0 .AND. C(1).NE.5.D0) GO TO 30
DO 25 K = 1, N
C(K+30) = DABS(C(K+30))
25 CONTINUE
30 CONTINUE
35 CONTINUE
C INITIALIZE RREB, DWARF, PREV XEND, FLAG, COUNTS
C(10) = 2.D0**(-56)
C(11) = 1.D-35
C SET PREVIOUS XEND INITIALLY TO INITIAL VALUE OF X
C(20) = X
DO 40 K = 21, 24
C(K) = 0.D0
40 CONTINUE
GO TO 50
C CASE 2 - NORMAL RE-ENTRY (IND .EQ. 3)
C .........ABORT IF XEND REACHED, AND EITHER X CHANGED OR XEND NOT
45 IF (C(21).NE.0.D0 .AND.
+ (X.NE.C(20) .OR. XEND.EQ.C(20))) GO TO 500
C RE-INITIALIZE FLAG
C(21) = 0.D0
GO TO 50
C CASE 3 - RE-ENTRY FOLLOWING AN INTERRUPT (IND .EQ. 4 TO 6)
C TRANSFER CONTROL TO THE APPROPRIATE RE-ENTRY POINT..........
C THIS HAS ALREADY BEEN HANDLED BY THE COMPUTED GO TO .
C END CASES V
50 CONTINUE
C
C END INITIALIZATION, ETC.
C
C ******************************************************************
C * LOOP THROUGH THE FOLLOWING 4 STAGES, ONCE FOR EACH TRIAL STEP *
C * UNTIL THE OCCURRENCE OF ONE OF THE FOLLOWING *
C * (A) THE NORMAL RETURN (WITH IND .EQ. 3) ON REACHING XEND IN *
C * STAGE 4 *
C * (B) AN ERROR RETURN (WITH IND .LT. 0) IN STAGE 1 OR STAGE 4 *
C * (C) AN INTERRUPT RETURN (WITH IND .EQ. 4, 5 OR 6), IF *
C * REQUESTED, IN STAGE 1 OR STAGE 4 *
C ******************************************************************
C
99999 CONTINUE
C
C ***************************************************************
C * STAGE 1 - PREPARE - DO CALCULATIONS OF HMIN, HMAX, ETC., *
C * AND SOME PARAMETER CHECKING, AND END UP WITH SUITABLE *
C * VALUES OF HMAG, XTRIAL AND HTRIAL IN PREPARATION FOR TAKING *
C * AN INTEGRATION STEP. *
C ***************************************************************
C
C***********ERROR RETURN (WITH IND=-1) IF NO OF FCN EVALS TOO GREAT
IF (C(7).EQ.0.D0 .OR. C(24).LT.C(7)) GO TO 100
IND = -1
RETURN
100 CONTINUE
C
C CALCULATE SLOPE (ADDING 1 TO NO OF FCN EVALS) IF IND .NE. 6
IF (IND .EQ. 6) GO TO 105
CALL FCN(N, X, Y, W(1,1))
C(24) = C(24) + 1.D0
105 CONTINUE
C
C CALCULATE HMIN - USE DEFAULT UNLESS VALUE PRESCRIBED
C(13) = C(3)
IF (C(3) .NE. 0.D0) GO TO 165
C CALCULATE DEFAULT VALUE OF HMIN
C FIRST CALCULATE WEIGHTED NORM Y - C(12) - AS SPECIFIED
C BY THE ERROR CONTROL INDICATOR C(1)
TEMP = 0.D0
IF (C(1) .NE. 1.D0) GO TO 115
C ABSOLUTE ERROR CONTROL - WEIGHTS ARE 1
DO 110 K = 1, N
TEMP = DMAX1(TEMP, DABS(Y(K)))
110 CONTINUE
C(12) = TEMP
GO TO 160
115 IF (C(1) .NE. 2.D0) GO TO 120
C RELATIVE ERROR CONTROL - WEIGHTS ARE 1/DABS(Y(K)) SO
C WEIGHTED NORM Y IS 1
C(12) = 1.D0
GO TO 160
120 IF (C(1) .NE. 3.D0) GO TO 130
C WEIGHTS ARE 1/MAX(C(2),ABS(Y(K)))
DO 125 K = 1, N
TEMP = DMAX1(TEMP, DABS(Y(K))/C(2))
125 CONTINUE
C(12) = DMIN1(TEMP, 1.D0)
GO TO 160
130 IF (C(1) .NE. 4.D0) GO TO 140
C WEIGHTS ARE 1/MAX(C(K+30),ABS(Y(K)))
DO 135 K = 1, N
TEMP = DMAX1(TEMP, DABS(Y(K))/C(K+30))
135 CONTINUE
C(12) = DMIN1(TEMP, 1.D0)
GO TO 160
140 IF (C(1) .NE. 5.D0) GO TO 150
C WEIGHTS ARE 1/C(K+30)
DO 145 K = 1, N
TEMP = DMAX1(TEMP, DABS(Y(K))/C(K+30))
145 CONTINUE
C(12) = TEMP
GO TO 160
150 CONTINUE
C DEFAULT CASE - WEIGHTS ARE 1/MAX(1,ABS(Y(K)))
DO 155 K = 1, N
TEMP = DMAX1(TEMP, DABS(Y(K)))
155 CONTINUE
C(12) = DMIN1(TEMP, 1.D0)
160 CONTINUE
C(13) = 10.D0*DMAX1(C(11),C(10)*DMAX1(C(12)/TOL,DABS(X)))
165 CONTINUE
C
C CALCULATE SCALE - USE DEFAULT UNLESS VALUE PRESCRIBED
C(15) = C(5)
IF (C(5) .EQ. 0.D0) C(15) = 1.D0
C
C CALCULATE HMAX - CONSIDER 4 CASES
C CASE 1 BOTH HMAX AND SCALE PRESCRIBED
IF (C(6).NE.0.D0 .AND. C(5).NE.0.D0)
+ C(16) = DMIN1(C(6), 2.D0/C(5))
C CASE 2 - HMAX PRESCRIBED, BUT SCALE NOT
IF (C(6).NE.0.D0 .AND. C(5).EQ.0.D0) C(16) = C(6)
C CASE 3 - HMAX NOT PRESCRIBED, BUT SCALE IS
IF (C(6).EQ.0.D0 .AND. C(5).NE.0.D0) C(16) = 2.D0/C(5)
C CASE 4 - NEITHER HMAX NOR SCALE IS PROVIDED
IF (C(6).EQ.0.D0 .AND. C(5).EQ.0.D0) C(16) = 2.D0
C
C***********ERROR RETURN (WITH IND=-2) IF HMIN .GT. HMAX
IF (C(13) .LE. C(16)) GO TO 170
IND = -2
RETURN
170 CONTINUE
C
C CALCULATE PRELIMINARY HMAG - CONSIDER 3 CASES
IF (IND .GT. 2) GO TO 175
C CASE 1 - INITIAL ENTRY - USE PRESCRIBED VALUE OF HSTART, IF
C ANY, ELSE DEFAULT
C(14) = C(4)
IF (C(4) .EQ. 0.D0) C(14) = C(16)*TOL**(1./6.)
GO TO 185
175 IF (C(23) .GT. 1.D0) GO TO 180
C CASE 2 - AFTER A SUCCESSFUL STEP, OR AT MOST ONE FAILURE,
C USE MIN(2, .9*(TOL/EST)**(1/6))*HMAG, BUT AVOID POSSIBLE
C OVERFLOW. THEN AVOID REDUCTION BY MORE THAN HALF.
TEMP = 2.D0*C(14)
IF (TOL .LT. (2.D0/.9D0)**6*C(19))
+ TEMP = .9D0*(TOL/C(19))**(1./6.)*C(14)
C(14) = DMAX1(TEMP, .5D0*C(14))
GO TO 185
180 CONTINUE
C CASE 3 - AFTER TWO OR MORE SUCCESSIVE FAILURES
C(14) = .5D0*C(14)
185 CONTINUE
C
C CHECK AGAINST HMAX
C(14) = DMIN1(C(14), C(16))
C
C CHECK AGAINST HMIN
C(14) = DMAX1(C(14), C(13))
C
C***********INTERRUPT NO 1 (WITH IND=4) IF REQUESTED
IF (C(8) .EQ. 0.D0) GO TO 1111
IND = 4
RETURN
C RESUME HERE ON RE-ENTRY WITH IND .EQ. 4 ........RE-ENTRY..
1111 CONTINUE
C
C CALCULATE HMAG, XTRIAL - DEPENDING ON PRELIMINARY HMAG, XEND
IF (C(14) .GE. DABS(XEND - X)) GO TO 190
C DO NOT STEP MORE THAN HALF WAY TO XEND
C(14) = DMIN1(C(14), .5D0*DABS(XEND - X))
C(17) = X + DSIGN(C(14), XEND - X)
GO TO 195
190 CONTINUE
C HIT XEND EXACTLY
C(14) = DABS(XEND - X)
C(17) = XEND
195 CONTINUE
C
C CALCULATE HTRIAL
C(18) = C(17) - X
C
C END STAGE 1
C
C ***************************************************************
C * STAGE 2 - CALCULATE YTRIAL (ADDING 7 TO NO OF FCN EVALS). *
C * W(*,2), ... W(*,8) HOLD INTERMEDIATE RESULTS NEEDED IN *
C * STAGE 3. W(*,9) IS TEMPORARY STORAGE UNTIL FINALLY IT HOLDS *
C * YTRIAL. *
C ***************************************************************
C
TEMP = C(18)/1398169080000.D0
C
DO 200 K = 1, N
W(K,9) = Y(K) + TEMP*W(K,1)*233028180000.D0
200 CONTINUE
CALL FCN(N, X + C(18)/6.D0, W(1,9), W(1,2))
C
DO 205 K = 1, N
W(K,9) = Y(K) + TEMP*( W(K,1)*74569017600.D0
+ + W(K,2)*298276070400.D0 )
205 CONTINUE
CALL FCN(N, X + C(18)*(4.D0/15.D0), W(1,9), W(1,3))
C
DO 210 K = 1, N
W(K,9) = Y(K) + TEMP*( W(K,1)*1165140900000.D0
+ - W(K,2)*3728450880000.D0
+ + W(K,3)*3495422700000.D0 )
210 CONTINUE
CALL FCN(N, X + C(18)*(2.D0/3.D0), W(1,9), W(1,4))
C
DO 215 K = 1, N
W(K,9) = Y(K) + TEMP*( - W(K,1)*3604654659375.D0
+ + W(K,2)*12816549900000.D0
+ - W(K,3)*9284716546875.D0
+ + W(K,4)*1237962206250.D0 )
215 CONTINUE
CALL FCN(N, X + C(18)*(5.D0/6.D0), W(1,9), W(1,5))
C
DO 220 K = 1, N
W(K,9) = Y(K) + TEMP*( W(K,1)*3355605792000.D0
+ - W(K,2)*11185352640000.D0
+ + W(K,3)*9172628850000.D0
+ - W(K,4)*427218330000.D0
+ + W(K,5)*482505408000.D0 )
220 CONTINUE
CALL FCN(N, X + C(18), W(1,9), W(1,6))
C
DO 225 K = 1, N
W(K,9) = Y(K) + TEMP*( - W(K,1)*770204740536.D0
+ + W(K,2)*2311639545600.D0
+ - W(K,3)*1322092233000.D0
+ - W(K,4)*453006781920.D0
+ + W(K,5)*326875481856.D0 )
225 CONTINUE
CALL FCN(N, X + C(18)/15.D0, W(1,9), W(1,7))
C
DO 230 K = 1, N
W(K,9) = Y(K) + TEMP*( W(K,1)*2845924389000.D0
+ - W(K,2)*9754668000000.D0
+ + W(K,3)*7897110375000.D0
+ - W(K,4)*192082660000.D0
+ + W(K,5)*400298976000.D0
+ + W(K,7)*201586000000.D0 )
230 CONTINUE
CALL FCN(N, X + C(18), W(1,9), W(1,8))
C
C CALCULATE YTRIAL, THE EXTRAPOLATED APPROXIMATION AND STORE
C IN W(*,9)
DO 235 K = 1, N
W(K,9) = Y(K) + TEMP*( W(K,1)*104862681000.D0
+ + W(K,3)*545186250000.D0
+ + W(K,4)*446637345000.D0
+ + W(K,5)*188806464000.D0
+ + W(K,7)*15076875000.D0
+ + W(K,8)*97599465000.D0 )
235 CONTINUE
C
C ADD 7 TO THE NO OF FCN EVALS
C(24) = C(24) + 7.D0
C
C END STAGE 2
C
C ***************************************************************
C * STAGE 3 - CALCULATE THE ERROR ESTIMATE EST. FIRST CALCULATE *
C * THE UNWEIGHTED ABSOLUTE ERROR ESTIMATE VECTOR (PER UNIT *
C * STEP) FOR THE UNEXTRAPOLATED APPROXIMATION AND STORE IT IN *
C * W(*,2). THEN CALCULATE THE WEIGHTED MAX NORM OF W(*,2) AS *
C * SPECIFIED BY THE ERROR CONTROL INDICATOR C(1). FINALLY, *
C * MODIFY THIS RESULT TO PRODUCE EST, THE ERROR ESTIMATE (PER *
C * UNIT STEP) FOR THE EXTRAPOLATED APPROXIMATION YTRIAL. *
C ***************************************************************
C
C CALCULATE THE UNWEIGHTED ABSOLUTE ERROR ESTIMATE VECTOR
DO 300 K = 1, N
W(K,2) = ( W(K,1)*8738556750.D0
+ + W(K,3)*9735468750.D0
+ - W(K,4)*9709507500.D0
+ + W(K,5)*8582112000.D0
+ + W(K,6)*95329710000.D0
+ - W(K,7)*15076875000.D0
+ - W(K,8)*97599465000.D0)/1398169080000.D0
300 CONTINUE
C
C CALCULATE THE WEIGHTED MAX NORM OF W(*,2) AS SPECIFIED BY
C THE ERROR CONTROL INDICATOR C(1)
TEMP = 0.D0
IF (C(1) .NE. 1.D0) GO TO 310
C ABSOLUTE ERROR CONTROL
DO 305 K = 1, N
TEMP = DMAX1(TEMP,DABS(W(K,2)))
305 CONTINUE
GO TO 360
310 IF (C(1) .NE. 2.D0) GO TO 320
C RELATIVE ERROR CONTROL
DO 315 K = 1, N
TEMP = DMAX1(TEMP, DABS(W(K,2)/Y(K)))
315 CONTINUE
GO TO 360
320 IF (C(1) .NE. 3.D0) GO TO 330
C WEIGHTS ARE 1/MAX(C(2),ABS(Y(K)))
DO 325 K = 1, N
TEMP = DMAX1(TEMP, DABS(W(K,2))
+ / DMAX1(C(2), DABS(Y(K))) )
325 CONTINUE
GO TO 360
330 IF (C(1) .NE. 4.D0) GO TO 340
C WEIGHTS ARE 1/MAX(C(K+30),ABS(Y(K)))
DO 335 K = 1, N
TEMP = DMAX1(TEMP, DABS(W(K,2))
+ / DMAX1(C(K+30), DABS(Y(K))) )
335 CONTINUE
GO TO 360
340 IF (C(1) .NE. 5.D0) GO TO 350
C WEIGHTS ARE 1/C(K+30)
DO 345 K = 1, N
TEMP = DMAX1(TEMP, DABS(W(K,2)/C(K+30)))
345 CONTINUE
GO TO 360
350 CONTINUE
C DEFAULT CASE - WEIGHTS ARE 1/MAX(1,ABS(Y(K)))
DO 355 K = 1, N
TEMP = DMAX1(TEMP, DABS(W(K,2))
+ / DMAX1(1.D0, DABS(Y(K))) )
355 CONTINUE
360 CONTINUE
C
C CALCULATE EST - (THE WEIGHTED MAX NORM OF W(*,2))*HMAG*SCALE
C - EST IS INTENDED TO BE A MEASURE OF THE ERROR PER UNIT
C STEP IN YTRIAL
C(19) = TEMP*C(14)*C(15)
C
C END STAGE 3
C
C ***************************************************************
C * STAGE 4 - MAKE DECISIONS. *
C ***************************************************************
C
C SET IND=5 IF STEP ACCEPTABLE, ELSE SET IND=6
IND = 5
IF (C(19) .GT. TOL) IND = 6
C
C***********INTERRUPT NO 2 IF REQUESTED
IF (C(9) .EQ. 0.D0) GO TO 2222
RETURN
C RESUME HERE ON RE-ENTRY WITH IND .EQ. 5 OR 6 ...RE-ENTRY..
2222 CONTINUE
C
IF (IND .EQ. 6) GO TO 410
C STEP ACCEPTED (IND .EQ. 5), SO UPDATE X, Y FROM XTRIAL,
C YTRIAL, ADD 1 TO THE NO OF SUCCESSFUL STEPS, AND SET
C THE NO OF SUCCESSIVE FAILURES TO ZERO
X = C(17)
DO 400 K = 1, N
Y(K) = W(K,9)
400 CONTINUE
C(22) = C(22) + 1.D0
C(23) = 0.D0
C**************RETURN(WITH IND=3, XEND SAVED, FLAG SET) IF X .EQ. XEND
IF (X .NE. XEND) GO TO 405
IND = 3
C(20) = XEND
C(21) = 1.D0
RETURN
405 CONTINUE
GO TO 420
410 CONTINUE
C STEP NOT ACCEPTED (IND .EQ. 6), SO ADD 1 TO THE NO OF
C SUCCESSIVE FAILURES
C(23) = C(23) + 1.D0
C**************ERROR RETURN (WITH IND=-3) IF HMAG .LE. HMIN
IF (C(14) .GT. C(13)) GO TO 415
IND = -3
RETURN
415 CONTINUE
420 CONTINUE
C
C END STAGE 4
C
GO TO 99999
C END LOOP
C
C BEGIN ABORT ACTION
500 CONTINUE
C
WRITE(6,505) IND, TOL, X, N, C(13), XEND, NW, C(16), C(20),
+ C(22), C(23), C(24), (Y(K), K = 1, N)
505 FORMAT( /// 1H0, 58HCOMPUTATION STOPPED IN DVERK WITH THE FOLLOWIN
+ G VALUES-
+ / 1H0, 5HIND =, I4, 5X, 6HTOL =, 1PD13.6, 5X, 11HX =,
+ 1PD22.15
+ / 1H , 5HN =, I4, 5X, 6HHMIN =, 1PD13.6, 5X, 11HXEND =,
+ 1PD22.15
+ / 1H , 5HNW =, I4, 5X, 6HHMAX =, 1PD13.6, 5X, 11HPREV XEND =,
+ 1PD22.15
+ / 1H0, 14X, 27HNO OF SUCCESSFUL STEPS =, 0PF8.0
+ / 1H , 14X, 27HNO OF SUCCESSIVE FAILURES =, 0PF8.0
+ / 1H , 14X, 27HNO OF FUNCTION EVALS =, 0PF8.0
+ / 1H0, 23HTHE COMPONENTS OF Y ARE
+ // (1H , 1P5D24.15) )
C
STOP
C
C END ABORT ACTION
C
END
C
C
subroutine fcn(neq, t, y, yprime)
implicit none
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
include 'fcn.inc'
c-----------------------------------------------------------------------
c Arguments of subroutine and indices
c Note: dimension of y is nst*(nst+3) (nst = dims(4))
c dimension of yprime is neq = dims(1)
c-----------------------------------------------------------------------
integer ic, jc
integer minst, maxst, nst
real*8 t
integer neq
real*8 y, yprime,TOTIME,HBAR, V,W,TOCYC,THETA
dimension neq(4)
dimension y(1), yprime(1)
c-----------------------------------------------------------------------
c Temporaries for dimension of system
c-----------------------------------------------------------------------
minst = neq(2)
maxst = neq(3)
nst = neq(4)
c-----------------------------------------------------------------------
c Define derivatives
c Note: y(ic) = a(ic+maxst-nst)
c y(ic+nst) = b(ic+maxst-nst)
c y(ic+2*nst) = E(ic+maxst-nst)
c y(nst*ic+jc+2*nst) = x(jc+maxst-nst,ic+maxst-nst)
c
c Note: we will replace x(j,i) with x(i,j) in the calculations
c below for better vectorization. We can do this as long as
c the x(j,i) are real.
c-----------------------------------------------------------------------
if (t .lt. totime) then
do ic = 1, nst
yprime(ic) = y(ic+2*nst)*y(ic+nst)/hbar
yprime(ic+nst) = - y(ic+2*nst)*y(ic)/hbar
do jc = 1, nst
yprime(ic) = yprime(ic) -
& V*(sin(w*t/(4.0d0*tocyc))**2)*
& cos(w*t+theta)*
& y(nst*ic+jc+2*nst)*y(jc+nst)/hbar
yprime(ic+nst) = yprime(ic+nst)
& +V*(sin(w*t/(4.0d0*tocyc))**2)*
& cos(w*t+theta)*
& y(nst*ic+jc+2*nst)*y(jc)/hbar
end do
end do
else
do ic = 1, nst
yprime(ic) = y(ic+2*nst)*y(ic+nst)/hbar
yprime(ic+nst) = - y(ic+2*nst)*y(ic)/hbar
do jc = 1, nst
yprime(ic) = yprime(ic) -
& V*cos(w*t+theta)*y(nst*ic+jc+2*nst)
& *y(jc+nst)/hbar
yprime(ic+nst) = yprime(ic+nst) +
& V*cos(w*t+theta)*y(nst*ic+jc+2*nst)
& *y(jc)/hbar
end do
end do
end if
return
end
I am getting a run-time error when trying to compile my program and as I
am not very experienced in Fortran I am not sure what the error is or how
to fix it. The program has been lifted out of a book called "Interneal
Combustion Engines:Applied Thermosciences" and is a program to model the
heat release rate of a an engine. It requires a couplke of subroutines,
DVERK and LEQIF which are IMSL routines and I do not have. I have found
the code for DVERK but unfortunately have not found the code for LEQIF
which is why I assume I am receiving the run-time error. the error is a
zero or negative argument to a logarithm routine.
FARG - in file release2_good.for at line 664 [+1738]
main - in file release2_good.for at line 41 [+0190]
The main code is listed below and any help in resolving this matter would
be greatly appreciated.
Thanks
Chris
C GIVEN IN BLOCK DATA SUBPROGRAM
C COMPRESSION RATIO - R
C BORE - B (CM)
C STROKE - S (CM)
C HALF-STROKE TO ROD RATIO - EPS
C ENGINE SPEED - RPM
C HEAT TRANSFER COEFFICIENT - H (J/S/M**2/K)
C BLOWBY COEFFICIENT - C (L/S)
C EQUIVALENCE RATIO - PHI
C RESIDUAL FRACTION - F
C INITIAL PRESSURE - P1 (BAR)
C INITIAL TEMPERATURE - T1 (K)
C WALL TEMPERATURE - TW (K)
C BURN ANGLE - THETAB (DEG)
C START OF HEAT RELEASE - THETAS (DEG ATDC)
C
C FIND:
C INDICATED THERMAL EFFICIENCY - ETA
C INDICATED MEAN EFFECTIVE PRESSURE - IMEP (BAR)
C
C NOTES:
C 1. COSINE BURNING LAW IS EMPLOYED, EQN. 4.61
C ` 2. FUEL - GASOLINE (C7H17)
C
EXTERNAL RATES
REAL IMEP, MNOT, M
COMMON /ENGINE/ R, B, S, EPS, RPM, H, C, THETAB, THETAS,
& PHI, F, P1, T1, TW, MNOT, OMEGA, UNOT, UFINAL
DIMENSION Y1(6), Y2(10), Y(6), COM(24), WM(6,9)
DATA AO/47870./, FS/.06548/, Y/1., 10000., 350., 3*0.0/,
& N/6/, THETA/-180./, THETAE/-179./, TOL/.001/, IND/1/,
& NW/6/, PI/3.141593/, Y2/10*0.0/
C
C LABEL HEADER ON OUTPUT PAGE AND PRINT INITIAL CONDITIONS
C
OMEGA = RPM*PI/30
CALL AUXLRY(THETA, V, X, EM)
Y(1) = P1
Y(3) = T1
CALL FARG(Y(1), Y(3), PHI, F, HU, U, VU, SU, Y1, CP, DLVLT,
& DLVLP)
MNOT = V/VU
UNOT = MNOT*U
WRITE(6,10)
10 FORMAT(5X,'THETA',4X,'VOLUME',3X,'BURN FRAC',2X,'PRESSURE',3X,
& 'BURN TEMP',2X,'T UNBURN',3X,'WORK',3X,'HEAT LOSS',5X,
& 'MASS',3X,'H-LEAK')
WRITE(6,15)
15 FORMAT(4X,'DEG ATDC',2X,'CM**3',6X,'--',8X,'BARS',7X,'KELVIN',
& 5X,'KELVIN',3X,'JOULES',4X,'JOULES',6X,'GRAMS',2X,
& 'JOULES')
WRITE(6,20) THETA, V, X, (Y(I),I=1,5), MNOT, Y(6)
20 FORMAT(4X,F5.0,6X,F5.0,4X,F5.3,6X,F5.2,7X,F5.0,5X,F5.0,
& 5X,F5.0,3X,F5.0,7X,F5.3,3X,F5.2)
C
C INTEGRATE THE RATE EQUATIONS 4.80 - 4.85 AND PRINT ANSWERS
C EVERY 10 DEGREES
C
DO 30 I = 1,36
DO 25 J = 1,10
CALL DVERK (N, RATES, THETA, Y, THETAE, TOL, IND, COM, NW, WM)
CALL AUXLRY(THETAE, V, X, EM)
M = EM*MNOT
THETA = THETAE
THETAE = THETA + 1.
IF(THETAS .GE. THETA .AND. THETAS .LT. THETAE)
& CALL TINITL(Y(1), T(3), PHI, F, Y(2))
IF(THETA .GT. THETAS + THETAB) Y(3) = 10000.
25 CONTINUE
WRITE(6,20) THETA, V, X, (Y(J),J=1,5), M, Y(6)
30 CONTINUE
C
C COMPUTE AND WRITE THE EAT, IMEP AND ERRORS
C
CALL ECP(Y(1), Y(2), PHI, HB, U, VB, SB, Y2, CP, DLVLT, DLVLP,
& IER)
UFINAL = U*M
ERROR1 = 1.0 - VB*M/V
ERROR2 = 1.0 + Y(4)/(UFINAL - UNOT + Y(5) + Y(6))
ETA = Y(4)/MNOT*(1. + PHI*FS*(1. - F))/PHI/FS/(1. - F)/AO
IMEP = Y(4)/(PI/4.*B**2*S)*10.
WRITE(6,40) ETA, IMEP, ERROR1, ERROR2
40 FORMAT('ETA = ',1PE11.4,' IMEP = ',1PE11.4,'ERROR1 = ',1PE11.4,
& ' ERROR2 = ',1PE11.4)
STOP
END
C
SUBROUTINE RATES(N, THETA, Y, YPRIME)
DIMENSION Y1(6), Y2(10)
REAL Y(N), YPRIME(N), M, MNOT
COMMON /ENGINE/ R, BORE, STROKE, EPS, RPM, HEAT, CEE, THETAB,
& THETAS, PHI, F, P1, T1, TW, MNOT, OMEGA, UNOT,
& UFINAL
DATA PI/3.141593/, Y2/10*0.0/
CALL AUXLRY(THETA, V, X, EM)
M = EM*MNOT
RAD = THETA*PI/180.
DUMB = SQRT(1. - (EPS*SIN(RAD))**2)
DV = PI/8.*BORE**2*STROKE*SIN(RAD)*(1. + EPS*COS(RAD)/DUMB)
A = (DV + V*CEE/OMEGA)/M
C1 = HEAT*(PI*BORE**2/2. + 4.*V/BORE)/OMEGA/M/10000.
C0 = SQRT(X)
C
C THREE DIFFERENT COMPUTATIONS ARE REQUIRED DEPENDING
C UPON THE SIZE OF THE MASS FRACTION BURNED
C
IF( X .GT. .999 ) GO TO 20
IF( X .GT. .001 ) GO TO 10
C
C COMPRESSION
C
CALL FARG(Y(1), Y(3), PHI, F, HL, U, VU, S, Y1, CP, DLVLT,
& DLVLP)
B = C1*VU/CP*DLVLT*(1.0 - TW/Y(3))
C = 0.
D = 0.
E = VU**2/CP/Y(3)*DLVLT**2 + 10.*VU/Y(1)*DLVLP
YPRIME(1) = (A + B + C)/(D + E)*10.
YPRIME(2) = 0.
YPRIME(3) = -C1/CP*(Y(3) - TW) + VU/CP*DLVLT*YPRIME(1)/10.
GO TO 30
C
C COMBUSTION
C
10 CALL FARG(Y(1), Y(3), PHI, F, HU, U, VU, S, Y1, CPU, DLVLTU,
& DLVLP)
CALL ECP(Y(1), Y(2), PHI, HB, U, VB, S, Y2, CPB, DDLVLTB,
& DLVLPB, IER)
B = C1*(VB/CPB*DLVLTB*C0*(1. - TW/Y(2)) + VU/CPU*DLVLTU*
& (1. - C0)*(1.- TW/Y(3)))
DX = 0.5*SIN(PI*(THETA - THETAS)/THETAB)*180./THETAB
C = -(VB - VU)*DX - VB/Y(2)*DLVLTB*(HU - HB)/CPB*(DX - (X -
& X**2)*CEE/OMEGA)
D = X*(VB**2/CPB/Y(2)*DLVLTB**2 + 10.*VB/Y(1)*DLVLPB)
E = (1. - X)*(VU**2/CPU/Y(3)*DLVLTU**2 + 10.*VU/Y(1)*DLVLPU)
HL = (1.0 - X**2)*HU + X**2*HB
YPRIME(1) = (A + B + C)/(D + E)*10.
YPRIME(2) = -C1/CPB/C0*(Y(2) - TW) + VB/CPB*DLVLTB*
& YPRIME(1)/10. + (HU - HB)/CPB*(DX/X - (1. - X)*
& CEE/OMEGA)
YPRIME(3) = -C1/CPU/(1. + CO)*(Y(3) - TW) + VU/CPU*DLVLTU*
& YPRIME(1)/10.
GO TO 30
C
C EXPANSION
C
20 CALL ECP(Y(1), Y(2), PHI, HL, U, VB, S, Y2, CP, DLVLT, DLVLP,
& IER)
B = C1*VB/CP*DLVLT*(1. - TW/Y(2))
C = 0.
D = VB**2/CP/Y(2)*DLVLT**2 + 10.*VB/Y(1)*DLVLP
E = 0.
YPRIME(1) = (A + B + C)/(D + E)*10.
YPRIME(2) = -C1/CP*(Y(2) - TW) + VB/CP*DLVLT*YPRIME(1)/10.
YPRIME(3) = 0.
C
C ALL CASES
C
30 YPRIME(4) = Y(1)*DV/10.
YPRIME(5) = C1*M*(C0*(Y(2) - TW) + (1. - C0)*(Y(3) - TW))
YPRIME(6) = CEE*M/OMEGA*HL
DO 40 I = 1,N
40 YPRIME(I) = YPRIME(I)*PI/180.
RETURN
END
C
SUBROUTINE AUXLRY(THETA, V, X, EM)
COMMON /ENGINE/ R, B, S, EPS, RPM, H, C, THETAB, THETAS,
& PHI, F, P1, T1, TW, EMNOT, OMEGA, UNOT, UFINAL
DATA PI/3.141593/
RAD = THETA*PI/180.
VTDC = PI/4.*B**2*S/(R - 1.)
V = VTDC*(1. + (R - 1.)/2.*(1. - COS(RAD) + 1./EPS*(1. -
& SQRT(1. - (EPS*SIN(RAD))**2))))
X = 0.5*(1. - COS(PI*(THETA - THETAS)/THETAB))
IF(THETA .LE. THETAS) X = 0.
IF(THETA .GE. THETAS + THETAB) X = 1.
EM = EXP(-C*(RAD + PI)/OMEGA)
RETURN
END
C
BLOCK DATA
COMMON /ENGINE/ R, B, S, EPS, RPM, H, C, THETAB, THETAS,
& PHI, F, P1, T1, TW, EMNOT, OMEGA, UNOT, UFINAL
DATA R/10./, B/10./, S/8./, EPS/.25/, RPM/2000./, H/500./,
& C/0.8/, THETAB/60./, THETAS/-35./, PHI/0.8/, F/.1/,
& P1/1.0/, T1/350./, TW/420./
C
END
C
SUBROUTINE TINITL(P, TU, PHI, F, TB)
DIMENSION Y1(6), Y2(10)
DATA MAXITS/50/, TOL/.001/, Y2/10*0.0/
TB = 2000.
CALL FARG(P, TU, PHI, F, HU, U, V, S, Y1, CP, DLVLT, DLVLP)
DO 10 I = 1,MAXITS
CALL ECP(P, TB, PHI, HB, U, V, S, Y2, CP, DLVLT, DLVLP, IER)
DELT = (HU - HB)/CP
TB = TB + DELT
IF(ABS(DELT/TB) .LT. TOL) RETURN
10 CONTINUE
WRITE(6,20)
20 FORMAT('CONVERGENCE FAILURE IN TINITL')
RETURN
END
SUBROUTINE ECP (P, T, PHI, H, U, V, S, Y, CP, DLVLT, DLVLP, IER)
C
C PURPOSE:
C COMPUTE THE PROPERTIES OF EQUILIBRIUM COMBUSTION
C PRODUCTS
C
C GIVEN:
C P - PRESSURE (BARS)
C T - TEMPERATURE )K)
C PHI - FUEL AIR EQUIVALENCE RATIO
C
C RETURNS:
C H - ENTHALPY (J/G)
C U - INTERNAL ENERGY (J/G)
C V - SPECIFIC VOLUME (CM**3/G)
C S - ENTROPY (J/G/K)
C Y - A TEN DIMENSIONAL COMPOSITION VECTOR OF
C MOLE FRACTIONS
C 1=CO2. 2=H2O, 3=N2, 4=O2, 5=CO, 6=H2
C 7=H, 8=O, 9=OH, 10=NO
C CP - SPECIFIC HEAT AT CONSTANT PRESSURE (J/G/K)
C DLVLT - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C TEMPERATURE AT CONSTANT PRESSURE
C DLVLP - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C PRESSURE AT CONSTANT TEMPERATURE
C IER - AN ERROR MESSAGE
C 0 - NORMAL
C 1 - CONVERGENCE FAILURE: COMPOSITION LOOP
C 2 - CALLED WITH TOO HIGH AN EQUIVALENCE RATIO; C(S)
C AND OTHER SPECIES WOULD FORM
C
REAL M, MW, K, KP, MT, MP
LOGICAL CHECK
DIMENSION M(10), K(6), C(6), Y(10), D(3), B(4), A(4,4), Y0(6),
& KP(5,6), DFDT(4), DCDT(6), DKDT(6), DCDP(6),
& DYDT(10), DYDP(10), A0(7,10), DFDP(4), CP0(10),
& H0(10), S0(10), WK(12)
C
C FUEL DATA
C
DATA ALPHA/7.0/, BETA/17.0/, GAMMA/0./, DELTA/0./
C
C SPECIFIC HEAT DATA FROM GORDON AND MCBRDE (1971)
C
DATA A0/
1 .44608041E+01, .30981719E-02, -.12392571E-05, .22741325E-09,
1 -.15525954E-13, -.48961442E+05, -.98635982E+00,
2 .27167633E+01, .29451374E-02, -.80224374E-06, .10226682E-09,
2 -.48472145E-14, -.29905826E+05, .66305671E+01,
3 .28963194E+01, .15154866E-02, -.57235277E-06, .99807393E-10,
3 -.65223555E-14, -.90586184E+03, .61615148E+01,
4 .36219535E+01, .72518264E-03, -.19652228E-06, .36201558E-10,
4 -.28945627E-14, -.12019825E+04, .36150960E+01,
5 .29840696E+01, .14891390E-02, -.57899684E-06, .36201558E-10,
5 -.69353550E-14, -.14245228E+05, .63479156E+01,
6 .31001901E+01, .51119464E-03, .52644210E-07, -.34909973E-10,
6 .36945345E-14, -.87738042E+03, -.19629421E+01,
7 .25000000E+01, 0., 0., 0., 0., .25471627E+05, -.46011763E+00,
8 .55420596E+01, -.27550619E-04, -.31028033E-08, .45510674E-11,
8 -.43680515E-15, .29230803E+03, .49203080E+01,
9 .29106427E+01, .95931650E-03, -.19441702E-06, .13756646E-10,
9 .14224542E-15, .39353815E+04, .54423445E+01,
# .31890000E+01, .13382281E-02, -.52899318E-06, .95919332E-10,
# -.64847932E-14, .98283290E+04, .67458126E+01/
C
C EQUILIBRIUM CONSTANT DATA FROM OLIKARA AND BORMAN (1975)
C
DATA KP/
1 .432168E+00, -.112464E+05, .267269E+01, -.745744E-04,
1 .242484E-08,
2 .310805E+00, -.129540E+05, .321779E+01, -.738336E-04,
2 .344645E-08,
3 -.141784E+00, -.213308E+04, .853461E+00, .355015E-04,
3 -.310227E-08,
4 .150879E-01, -.470959E+04, .646096E+00, .272805E-05,
4 -.154444E-08,
5 -.752364E+00, .124210E+05, -.260284E+01, .259556E-03,
5 -.162687E-07,
6 -.415302E-02, .148627E+05, -.475746E+01, .124699E-03,
6 -.900227E-08/
C
C MOLECULAR WEIGHTS
C
DATA M/44.01, 18.02, 28.008, 32., 28.01, 2.018, 1.009, 16.,
& 17.009, 30.004/
C
C OTHER DATA
C
DATA MAXITS/100/, TOL/3.0E-05/, RU/8.31434/
C
C MAKE SURE SOLID CARBON WILL NOT FORM
C
IER = 2
EPS = .210/(ALPHA + 0.25*BETA - 0.5*GAMMA)
IF(PHI .GT. (.210/EPS/(0.5*ALPHA - 0.5*GAMMA)) ) RETURN
IER = 0
C
C DECIDE IF AN INITIAL ESTIMATE TO THE COMPOSITION IS
C REQUIRED OR IF T<1000 IN WHICH CASE FARG WILL SUFFICE
C
C
SUM = 0
DO 10 I = 1,10
10 SUM = SUM + Y(I)
IF( T .GT. 1000. .AND. SUM .GT. 0.998) GO TO 20
CALL FARG(P,T,PHI,1.,H,U,V,S,Y0,CP,DLVLT,DLVLP)
DO 30 I = 1,10
30 Y(I) = 0
DO 40 I = 1,6
40 Y(I) = Y0(I)
IF( T .LE. 1000. ) RETURN
20 CONTINUE
C
C
C EVALUATE THE REQUISITE CONSTANTS
C
PATM = .9869233*P
DO 50 I = 1,6
50 K(I) = 10.**(KP(1,I)*ALOG(T/1000.) + KP(2,I)/T + KP(3,I) +
& KP(4,I)*T + KP(5,I)*T*T)
C(1) = K(1)/SQRT(PATM)
C(2) = K(2)/SQRT(PATM)
C(3) = K(3)
C(4) = K(4)
C(5) = K(5)*SQRT(PATM)
C(6) = K(6)*SQRT(PATM)
D(1) = BETA/ALPHA
D(2) = (GAMMA + .42/EPS/PHI)/ALPHA
D(3) = (DELTA + 1.58/EPS/PHI)/ALPHA
C
C
C SET UP ITERATION LOOP FOR NEWTON-RAPHSON ITERATION AND
C INTRODUCE TRICK TO PREVENT INSTABILITIES NEAR PHI = 1.0
C AT LOW TEMPERATURES ON 24 BIT MACHINES (VAX). SET
C TRICK = 1.0 ON 48 BIT MACHINES (CDC)
C
TRICK = 10.
CHECK = ABS(PHI - 1.0) .LT. TOL
IF(CHECK) PHI = PHI*(1.0 + SIGN(TOL,PHI - 1.0))
ICHECK = 0
DO 5 J = 1,MAXITS
DO 4 JJ = 1,10
4 Y(JJ) = AMIN1(1.0,AMAX1(Y(JJ),1.0E-25))
D76 = 0.5*C(1)/SQRT(Y(6))
D84 = 0.5*C(1)/SQRT(Y(4))
D94 = 0.5*C(3)*SQRT(Y(6)/Y(4))
D96 = 0.5*C(3)*SQRT(Y(4)/Y(6))
D103 = 0.5*C(4)*SQRT(Y(4)/Y(3))
D104 = 0.5*C(4)*SQRT(Y(3)/Y(4))
D24 = 0.5*C(5)*Y(6)/SQRT(Y(4))
D26 = C(5)*SQRT(Y(4))
D14 = 0.5*C(6)*Y(5)/SQRT(Y(4))
D15 = C(6)*SQRT(Y(4))
A(1,1) = 1. + D103
A(1,2) = D14 + D24 + 1. +D84 + D104 + D96
A(1,3) = D15 + 1.0
A(1,4) = D26 + 1.0 + D76 +D96
A(2,1) = 0.0
A(2,2) = 2.0*D24 + D94 - D(1)*D14
A(2,3) = -D(1)*D15 - D(1)
A(2,4) = 2.0*D26 + 2.0 + D76 + D96
A(3,1) = D103
A(3,2) = 2.0*D14 + D24 + 2.0 + D84 + D04 + D104 - D(2)*D14
A(3,3) = 2.0*D15 + 1.0 - D(2)*D15 - D(2)
A(3,4) = D26 + D96
A(4,1) = 2.0 +D103
A(4,2) = D104 - D(3)*D14
A(4,3) = -D(3)*D15 - D(3)
A(4,4) = 0
C
C SOLVE THE MATRIX EQUATIONS 3.81 FOR COMPOSITION CORRECTIONS
C
SUM = 0
DO 55 I = 1,10
55 SUM = SUM + Y(I)
B(1) = -(SUM - 1.0)
B(2) = -(2.0*Y(2) + 2.0*Y(6) + Y(7) + Y(9) - D(1)*Y(1) -
& D(1)*Y(5))
B(3) = -(2.0*Y(1) + Y(2) + 2.0*Y(4) + Y(5) + Y(8) + Y(9)
& + Y(10) - D(2)*Y(1) - D(2)*Y(5))
B(4) = -(2.0*Y(3) + Y(10) - D(3)*Y(1) - D(3)*Y(5))
CALL LEQIF(A, 4, 4, 4, B, 4, 1, 0, WK, IERROR)
ERROR = 0.
DO 56 L = 3,6
LL = L - 2
Y(L) = Y(L) + B(LL)/TRICK
ERROR = AMAX1(ERROR,ABS(B(LL)))
Y(L) = AMIN(1.0,AMAX(Y(L),1.0E-25))
56 CONTINUE
Y(7) = C(1)*SQRT(Y(6))
Y(8) = C(2)*SQRT(Y(4))
Y(9) = C(3)*SQRT(Y(4)*Y(6))
Y(10) = C(4)*SQRT(Y(4)*Y(3))
Y(2) = C(5)*SQRT(Y(4))*Y(6)
Y(1) = C(6)*SQRT(Y(4))*Y(5)
IF(ERROR .LT. TOL) ICHECK = ICHECK +1
IF(ERROR .LT. TOL .AND. ICHECK .GE. 2) GO TO 57
5 CONTINUE
IER = 1
C
C COMPUTE THE REQUISITE CONSTANTS TO FIND PARTIAL DERIVATIVES
C
57 DO 60 I = 1,6
60 DKDT(I) = 2.302585*K(I)*(KP(1,I)/T - KP(2,I)/T/T + KP(4,I)
& + 2.0*KP(5,I)*T)
DCDT(1) = DKDT(1)/SQRT(PATM)
DCDT(2) = DKDT(2)/SQRT(PATM)
DCDT(3) = DKDT(3)
DCDT(4) = DKDT(4)
DCDT(5) = DKDT(5)*SQRT(PATM)
DCDT(6) = DKDT(6)*SQRT(PATM)
DCDP(1) = -0.5*C(1)/P
DCDP(2) = -0.5*C(2)/P
DCDP(5) = 0.5*C(5)/P
DCDP(6) = 0.5*C(6)/P
C
X1 = Y(1)/C(6)
X2 = Y(2)/C(5)
X7 = Y(7)/C(1)
X8 = Y(8)/C(2)
X9 = Y(9)/C(3)
X10 = Y(10)/C(4)
DFDT(1) = DCDT(6)*X1 + DCDT(5)*X2 + DCDT(1)*X7 +DCDT(2)*X8
& + DCDT(3)*X9 + DCDT(4)*X10
DFDT(2) = 2.0*DCDT(5)*X2 + DCDT(1)*X7 + DCDT(3)*X9 -
& D(1)*DCDT(6)*X1
DFDT(3) = 2.0*DCDT(6)*X1 + DCDT(5)*X2 + DCDT(2)*X8
& + DCDT(3)*X9 + DCDT(4)*X10 - D(2)*DCDT(6)*X1
DFDT(4) = DCDT(4)*X10 - D(3)*DCDT(6)*X1
DFDP(1) = DCDP(6)*X1 + DCDP(5)*X2 + DCDP(1)*X7 + DCDP(2)*X8
DFDP(2) = 2.0*DCDP(5)*X2 + DCDP(1)*X7 - D(1)*DCDP(6)*X1
DFDP(3) = 2.0*DCDP(6)*X1 + DCDP(5)*X2 + DCDP(2)*X8 -
& D(2)*DCDP(6)*X1
DFDP(4) = -D(3)*DCDP(6)*X1
C
C SOLVE THE MATRIX EQUATIONS 3.91 FOR THE INDEPENDENT DERIVATIVES
C AND THEN DETERMINE THE DEPENDENT DERIVATICES
C
DO 70 I = 1,4
70 B(I) = -DFDT(I)
CALL LEQIF(A, 4, 4, 4, B, 4, 1, 1, WK, IERROR)
DYDT(3) = B(1)
DYDT(4) = B(2)
DYDT(5) = B(3)
DYDT(6) = B(4)
DYDT(1) = SQRT(Y(4))*Y(5)*DCDT(6) + D14*DYDT(4) + D15*DYDT(5)
DYDT(2) = SQRT(Y(4))*Y(6)*DCDT(5) + D24*DYDT(4) + D26*DYDT(6)
DYDT(7) = SQRT(Y(6))*DCDT(1) + D76*DYDT(6)
DYDT(8) = SQRT(Y(4))*DCDT(2) + D84*DYDT(4)
DYDT(9) = SQRT(Y(4)*Y(6))*DCDT(3) + D94*DYDT(4) + D96*DYDT(6)
DYDT(10) = SQRT(Y(4)*Y(3))*DCDT(4) + D104*DYDT(4) + D103*DYDT(3)
C
DO 80 I = 1,4
80 B(I) = -DFDP(I)
CALL LEQIF(A, 4, 4, 4, B, 4, 1, 1, WK, IERROR)
DYDP(3) = B(1)
DYDP(4) = B(2)
DYDP(5) = B(3)
DYDP(6) = B(4)
DYDP(1) = SQRT(Y(4))*Y(5)*DCDP(6) + D14*DYDP(4) + D15*DYDP(5)
DYDP(2) = SQRT(Y(4))*Y(6)*DCDP(5) + D24*DYDP(4) + D26*DYDP(6)
DYDP(7) = SQRT(Y(6))*DCDP(1) + D76*DYDP(6)
DYDP(8) = SQRT(Y(4))*DCDP(2) + D84*DYDP(4)
DYDP(9) = D94*DYDP(4) + D96*DYDP(6)
DYDP(10) = D104*DYDP(4) + D103*DYDP(3)
C
DO 90 I = 1,10
CP0(I) = A0(1,I) + A0(2,I)*T + A0(3,I)*T*T + A0(4,I)*T**3
& + A0(5,I)*T**4
H0(I) = A0(1,I) + A0(2,I)/2.*T + A0(3,I)/3.*T*T + A0(4,I)/4.*
& T**3 + A0(5,1)/5.*T**4 + A0(6,I)/T
90 S0(I) = A0(1,I)*ALOG(T) + A0(2,I)*T + A0(3,I)/2.*T**2
& + A0(4,I)/3.*T**3 + A0(5,I)/4.*T**4 +A0(7,I)
C
C Y(1), Y(2) ARE REEVALUATED TO CLEAN UP ANY ROUND-OFF ERRORS WHICH
C CAN OCCUR FOR THE LOW TEMPERATURE STOICHIOMETRIC SOLUTIONS
C
Y(1) = (2.*Y(3) + Y(10))/D(3) - Y(5)
Y(2) = (D(1)/D(3)*(2.*Y(3) + Y(10)) - 2.*Y(6) - Y(7) - Y(9))/2
C
MW = 0.
CP = 0.
H = 0.
S = 0.
MT = 0.
MP = 0.
DO 100 I = 1,10
IF( Y(I) .LE. 1.0E-25 ) GO TO 100
H = H + H0(I)*Y(I)
MW = MW + M(I)*Y(I)
MT = MT + M(I)*DYDT(I)
MP = MP + M(I)*DYDP(I)
CP = CP + Y(I)*CP0(I) + H0(I)*T*DYDT(I)
S = S + Y(I)*(S0(I) - ALOG(Y(I)))
100 CONTINUE
C
R = RU/MW
V = 10.*R*T/P
CP = R*(CP - H*T*MT/MW)
DLVLT = 1.0 + AMAX1(-T*MT/MW,0.)
DLVLP = -1.0 - AMAX1(P*MP/MW,0.)
H = H*R*T
S = R*(-ALOG(PATM) + 5)
U = H - R*T
RETURN
END
C
C
SUBROUTINE FARG (P, T, PHI, F, H, U, V, S, Y, CP, DLVLT, DLVLP)
C
C PURPOSE:
C COMPUTE THE PROPERTIES OF A FUEL, AIR, RESIDUAL
C GAS MIXTURE
C
C GIVEN:
C P - PRESSURE (BARS)
C T - TEMPERATURE (K)
C PHI - FUEL AIR EQUIVALENCE RATIO
C F - RESIDUAL MASS FRACTION
C
C RETURNS:
C H - ENTHALPY (J/G)
C U - INTERNAL ENERGY (J/G)
C V - SPECIFIC VOLUME (CM**3/g)
C S - ENTROPY (J/G/K)
C Y - A SIX DIMENSIONAL COMPOSITION VECTOR OF
C MOLE FRACTIONS
C 1=CO2, 2=H2O, 3=N2, 4O2, 5=CO, 6=H2
C CP - SPECIFIC HEAT AT CONSTANT PRESSURE (J/G/K)
C DLVLT - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C TEMPERATURE AT CONSTANT PRESSURE
C DLVLP - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
C PRESSURE AT CONSTANT TEMPERATURE
C
C REMARKS:
C 1. VALID FOR 300 < T < 1000
C 2. ASSUMES THE FUEL IS GASOLINE
C 3. ENTHALPIES OF O2, H2, N2 AND C(S) ARE SET TO ZERO
C AT 298K
C 4. THE FUEL MOLE FRACTION IS UNITY MINUS THE SUM OF Y(I)
C
REAL M, MW, MRES, MFA, MFUEL, K, NU, N2
LOGICAL RICH, LEAN
DIMENSION A(7,6), TABLE(6), M(6), NU(6), Y(6), CP0(6), H0(6),
& S0(6)
C
C FUEL DATA
C
DATA ALPHA/7.0/, BETA/17.0/, GAMMA/0./, DELTA/0./,
& A0/4.0652/, B0/6.0977E-02/, C0/-1.8801E-05/,
& D0/-3.5880E+04/, E0/15.45/
C
C TABLE 3.1, I = 1,6
C
DATA A/
1 .24007797E+01, .87350957E-02, -.6607878E-05, .20021861E-08,
1 .63274039E-15, -.48377527E+05, .96951457E+01,
2 .40701275E+01, -.11084499E-02, .41521180E-05, -.29637404E-08,
2 .80702103E-12, -.30279722E+05, -.32270046E+00,
3 .36748261E+01, -.12081500E-02, .23240102E-05, -.63217559E-09,
3 -.22577253E-12, -.10611588E+04, .23580424E+01,
4 .36255985E+01, -.18782184E-02, .70554544E-05, -.67635137E-08,
4 .21555993E-11, -.10475226E+04, .43052778E+01,
5 .37100928E+01, -.16190964E-02, .36923594E-05, -.20319674E-08,
5 .23953344E-12, -.14356310E+05, .29555355E+01,
6 .30574451E+01, .26765200E-02, -.58099162E-05, .55210391E-08,
6 -.18122739E-11, -.98890474E+03, -.22997056E+01/
C
C OTHER DATA
C
DATA RU/8.315/, TABLE/-1.,1.,0.,0.,1.,-1./, M/44.01, 18.02,
& 28.008, 32.000, 28.01, 2.018/
C
C COMPUTE RESIDUAL GAS COMPOSITION ACCORDING TO TABLE 3.3
C
RICH = PHI .GT. 1.0
LEAN = .NOT. RICH
DLVLT = 1.0
DLVLP = -1.0
EPS = .21/(ALPHA + 0.25*BETA - 0.5*GAMMA)
IF (RICH) GO TO 10
NU(1) = ALPHA*PHI*EPS
NU(2) = BETA*PHI*EPS/2
NU(3) = 0.79 + DELTA*PHI*EPS/2
NU(4) = 0.21*(1.0 - PHI)
NU(5) = 0.
NU(6) = 0.
DCDT = 0.
GO TO 20
10 Z = 1000./T
K = EXP(2.743 + Z*(-1.761 + Z*(-1.611 + Z*.2803)))
DKDT = -K*(-1.761 + Z*(-3.222 + Z*.8409))/1000.
A1 = 1.0 - K
B = .42 - PHI*EPS*(2.*ALPHA - GAMMA) + K*(.42*(PHI - 1.) +
& ALPHA*PHI*EPS)
C = -.42*ALPHA*PHI*EPS*(PHI - 1.0)*K
NU(5) = (-B + SQRT(B*B - 4.*A1*C))/2./A1
DCDT = DKDT*(NU(5)**2 - NU(5)*(.42*(PHI - 1.) + ALPHA*
& PHI*EPS) + .42*ALPHA*PHI*EPS*(PHI - 1.))/
& (2.*NU(5)*A1 + B)
NU(1) = ALPHA*PHI*EPS - NU(5)
NU(2) = .42 - PHI*EPS*(2.*ALPHA - GAMMA) + NU(5)
NU(3) = .79 + DELTA*PHI*EPS/2.
NU(4) = 0.
NU(6) = .42*(PHI - 1.) - NU(5)
C
C COMPUTE MOLE FRACTIONS AND MOLECULAR WEIGHT OF RESIDUAL
C
20 TMOLES = 0.
DO 30 I = 1,6
30 TMOLES = TMOLES + NU(I)
MRES = 0.
DO 40 I = 1,6
Y(I) = NU(I)/TMOLES
40 MRES = MRES + Y(I)*M(I)
C
C COMPUTE MOLE FRACTIONS AND MOLECULAR WEIGHT OF FURLE-AIR
C
FUEL = EPS*PHI/(1. + EPS*PHI)
O2 = .21/(1. + EPS*PHI)
N2 = .79/(1. + EPS*PHI)
MFA = FUEL*(12.01*ALPHA + 1.008*BETA + 16.*GAMMA + 14.01*
& DELTA) + 32.*O2 + 28.02*N2
C
C COMPUTE MOLE FRACTIONS OF FUEL-AIR-RESIDUAL GAS
C
YRES = F/(F + MRES/MFA*(1. - F))
DO 50 I = 1,6
50 Y(I) = Y(I)*YRES
YFUEL = FUEL*(1. - YRES)
Y(3) = Y(3) + N2*(1. - YRES)
Y(4) = Y(4) + O2*(1. - YRES)
C
C COMPUTE COMPONENT PROPERTIES
C
DO 60 I = 1,6
CP0(I) = A(1,I) + A(2,I)*T + A(3,I)*T**2 + A(4,I)*T**3 +
& A(5,I)*T**4
H0(I) = A(1,I) + A(2,I)/2.*T + A(3,I)/3.*T**2 + A(4,I)/4.*
& T**3 + A(5,I)/5.*T**4 + A(6,I)/T
60 S0(I) = A(1,I)*ALOG(T) + A(2,I)*T + A(3,I)/2.*T**2 +
& A(4,I)/3.*T**3 + A(5,I)/4.*T**4 + A(7,1)
C
MFUEL = 12.01*ALPHA + 1.008*BETA + 16.000*GAMMA + 14.01*DELTA
CPFUEL = A0 + B0*T + C0*T**2
HFUEL = A0 + B0/2.*T + C0/3.*T**2 + D0/T
S0FUEL = A0*ALOG(T) + B0*T + C0/2.*T**2 + E0
C
C COMPUTE PROPERTIES OF MIXTURE
C
H = HFUEL*YFUEL
S = (S0FUEL - ALOG(YFUEL))*YFUEL
CP = CPFUEL*YFUEL
MW = MFUEL*YFUEL
DO 70 I = 1,6
H = H + H0(I)*Y(I)
S = S + Y(I)*(S0(I) - ALOG(Y(I)))
CP = CP + CP0(I)*Y(I) + H0(I)*T*TABLE(I)*DCDT*YRES/TMOLES
70 MW = MW + Y(I)*M(I)
C
R = RU/MW
H = R*T*H
U = H - R*T
V = 10.*R*T/P
S = R*(-ALOG(.9869*P) + S)
CP = R*CP
C
RETURN
END
C
C
SUBROUTINE DVERK(N, FCN, X, Y, XEND, TOL, IND, C, NW, W)
INTEGER N, IND, NW, K
DOUBLE PRECISION X, Y(N), XEND, TOL, C(*), W(NW,9), TEMP
C Added external statement for fcn to avoid a warning message.
external fcn
C
C***********************************************************************
C *
C PURPOSE - THIS IS A RUNGE-KUTTA SUBROUTINE BASED ON VERNER'S *
C FIFTH AND SIXTH ORDER PAIR OF FORMULAS FOR FINDING APPROXIMATIONS TO *
C THE SOLUTION OF A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL *
C EQUATIONS WITH INITIAL CONDITIONS. IT ATTEMPTS TO KEEP THE GLOBAL *
C ERROR PROPORTIONAL TO A TOLERANCE SPECIFIED BY THE USER. (THE *
C PROPORTIONALITY DEPENDS ON THE KIND OF ERROR CONTROL THAT IS USED, *
C AS WELL AS THE DIFFERENTIAL EQUATION AND THE RANGE OF INTEGRATION.) *
C *
C VARIOUS OPTIONS ARE AVAILABLE TO THE USER, INCLUDING DIFFERENT *
C KINDS OF ERROR CONTROL, RESTRICTIONS ON STEP SIZES, AND INTERRUPTS *
C WHICH PERMIT THE USER TO EXAMINE THE STATE OF THE CALCULATION (AND *
C PERHAPS MAKE MODIFICATIONS) DURING INTERMEDIATE STAGES. *
C *
C THE PROGRAM IS EFFICIENT FOR NON-STIFF SYSTEMS. HOWEVER, A GOOD *
C VARIABLE-ORDER-ADAMS METHOD WILL PROBABLY BE MORE EFFICIENT IF THE *
C FUNCTION EVALUATIONS ARE VERY COSTLY. SUCH A METHOD WOULD ALSO BE *
C MORE SUITABLE IF ONE WANTED TO OBTAIN A LARGE NUMBER OF INTERMEDIATE *
C SOLUTION VALUES BY INTERPOLATION, AS MIGHT BE THE CASE FOR EXAMPLE *
C WITH GRAPHICAL OUTPUT. *
C *
C HULL-ENRIGHT-JACKSON 1/10/76 *
C *
C***********************************************************************
C *
C USE - THE USER MUST SPECIFY EACH OF THE FOLLOWING *
C *
C N NUMBER OF EQUATIONS *
C *
C FCN NAME OF SUBROUTINE FOR EVALUATING FUNCTIONS - THE SUBROUTINE *
C ITSELF MUST ALSO BE PROVIDED BY THE USER - IT SHOULD BE OF *
C THE FOLLOWING FORM *
C SUBROUTINE FCN(N, X, Y, YPRIME) *
C INTEGER N *
C DOUBLE PRECISION X, Y(N), YPRIME(N) *
C *** ETC *** *
C AND IT SHOULD EVALUATE YPRIME, GIVEN N, X AND Y *
C *
C X INDEPENDENT VARIABLE - INITIAL VALUE SUPPLIED BY USER *
C *
C Y DEPENDENT VARIABLE - INITIAL VALUES OF COMPONENTS Y(1), Y(2), *
C ..., Y(N) SUPPLIED BY USER *
C *
C XEND VALUE OF X TO WHICH INTEGRATION IS TO BE CARRIED OUT - IT MAY *
C BE LESS THAN THE INITIAL VALUE OF X *
C *
C TOL TOLERANCE - THE SUBROUTINE ATTEMPTS TO CONTROL A NORM OF THE *
C LOCAL ERROR IN SUCH A WAY THAT THE GLOBAL ERROR IS *
C PROPORTIONAL TO TOL. IN SOME PROBLEMS THERE WILL BE ENOUGH *
C DAMPING OF ERRORS, AS WELL AS SOME CANCELLATION, SO THAT *
C THE GLOBAL ERROR WILL BE LESS THAN TOL. ALTERNATIVELY, THE *
C CONTROL CAN BE VIEWED AS ATTEMPTING TO PROVIDE A *
C CALCULATED VALUE OF Y AT XEND WHICH IS THE EXACT SOLUTION *
C TO THE PROBLEM Y' = F(X,Y) + E(X) WHERE THE NORM OF E(X) *
C IS PROPORTIONAL TO TOL. (THE NORM IS A MAX NORM WITH *
C WEIGHTS THAT DEPEND ON THE ERROR CONTROL STRATEGY CHOSEN *
C BY THE USER. THE DEFAULT WEIGHT FOR THE K-TH COMPONENT IS *
C 1/MAX(1,ABS(Y(K))), WHICH THEREFORE PROVIDES A MIXTURE OF *
C ABSOLUTE AND RELATIVE ERROR CONTROL.) *
C *
C IND INDICATOR - ON INITIAL ENTRY IND MUST BE SET EQUAL TO EITHER *
C 1 OR 2. IF THE USER DOES NOT WISH TO USE ANY OPTIONS, HE *
C SHOULD SET IND TO 1 - ALL THAT REMAINS FOR THE USER TO DO *
C THEN IS TO DECLARE C AND W, AND TO SPECIFY NW. THE USER *
C MAY ALSO SELECT VARIOUS OPTIONS ON INITIAL ENTRY BY *
C SETTING IND = 2 AND INITIALIZING THE FIRST 9 COMPONENTS OF *
C C AS DESCRIBED IN THE NEXT SECTION. HE MAY ALSO RE-ENTER *
C THE SUBROUTINE WITH IND = 3 AS MENTIONED AGAIN BELOW. IN *
C ANY EVENT, THE SUBROUTINE RETURNS WITH IND EQUAL TO *
C 3 AFTER A NORMAL RETURN *
C 4, 5, OR 6 AFTER AN INTERRUPT (SEE OPTIONS C(8), C(9)) *
C -1, -2, OR -3 AFTER AN ERROR CONDITION (SEE BELOW) *
C *
C C COMMUNICATIONS VECTOR - THE DIMENSION MUST BE GREATER THAN OR *
C EQUAL TO 24, UNLESS OPTION C(1) = 4 OR 5 IS USED, IN WHICH *
C CASE THE DIMENSION MUST BE GREATER THAN OR EQUAL TO N+30 *
C *
C NW FIRST DIMENSION OF WORKSPACE W - MUST BE GREATER THAN OR *
C EQUAL TO N *
C *
C W WORKSPACE MATRIX - FIRST DIMENSION MUST BE NW AND SECOND MUST *
C BE GREATER THAN OR EQUAL TO 9 *
C *
C THE SUBROUTINE WILL NORMALLY RETURN WITH IND = 3, HAVING *
C REPLACED THE INITIAL VALUES OF X AND Y WITH, RESPECTIVELY, THE VALUE *
C OF XEND AND AN APPROXIMATION TO Y AT XEND. THE SUBROUTINE CAN BE *
C CALLED REPEATEDLY WITH NEW VALUES OF XEND WITHOUT HAVING TO CHANGE *
C ANY OTHER ARGUMENT. HOWEVER, CHANGES IN TOL, OR ANY OF THE OPTIONS *
C DESCRIBED BELOW, MAY ALSO BE MADE ON SUCH A RE-ENTRY IF DESIRED. *
C *
C THREE ERROR RETURNS ARE ALSO POSSIBLE, IN WHICH CASE X AND Y *
C WILL BE THE MOST RECENTLY ACCEPTED VALUES - *
C WITH IND = -3 THE SUBROUTINE WAS UNABLE TO SATISFY THE ERROR *
C REQUIREMENT WITH A PARTICULAR STEP-SIZE THAT IS LESS THAN OR *
C EQUAL TO HMIN, WHICH MAY MEAN THAT TOL IS TOO SMALL *
C WITH IND = -2 THE VALUE OF HMIN IS GREATER THAN HMAX, WHICH *
C PROBABLY MEANS THAT THE REQUESTED TOL (WHICH IS USED IN THE *
C CALCULATION OF HMIN) IS TOO SMALL *
C WITH IND = -1 THE ALLOWED MAXIMUM NUMBER OF FCN EVALUATIONS HAS *
C BEEN EXCEEDED, BUT THIS CAN ONLY OCCUR IF OPTION C(7), AS *
C DESCRIBED IN THE NEXT SECTION, HAS BEEN USED *
C *
C THERE ARE SEVERAL CIRCUMSTANCES THAT WILL CAUSE THE CALCULATIONS *
C TO BE TERMINATED, ALONG WITH OUTPUT OF INFORMATION THAT WILL HELP *
C THE USER DETERMINE THE CAUSE OF THE TROUBLE. THESE CIRCUMSTANCES *
C INVOLVE ENTRY WITH ILLEGAL OR INCONSISTENT VALUES OF THE ARGUMENTS, *
C SUCH AS ATTEMPTING A NORMAL RE-ENTRY WITHOUT FIRST CHANGING THE *
C VALUE OF XEND, OR ATTEMPTING TO RE-ENTER WITH IND LESS THAN ZERO. *
C *
C***********************************************************************
C *
C OPTIONS - IF THE SUBROUTINE IS ENTERED WITH IND = 1, THE FIRST 9 *
C COMPONENTS OF THE COMMUNICATIONS VECTOR ARE INITIALIZED TO ZERO, AND *
C THE SUBROUTINE USES ONLY DEFAULT VALUES FOR EACH OPTION. IF THE *
C SUBROUTINE IS ENTERED WITH IND = 2, THE USER MUST SPECIFY EACH OF *
C THESE 9 COMPONENTS - NORMALLY HE WOULD FIRST SET THEM ALL TO ZERO, *
C AND THEN MAKE NON-ZERO THOSE THAT CORRESPOND TO THE PARTICULAR *
C OPTIONS HE WISHES TO SELECT. IN ANY EVENT, OPTIONS MAY BE CHANGED ON *
C RE-ENTRY TO THE SUBROUTINE - BUT IF THE USER CHANGES ANY OF THE *
C OPTIONS, OR TOL, IN THE COURSE OF A CALCULATION HE SHOULD BE CAREFUL *
C ABOUT HOW SUCH CHANGES AFFECT THE SUBROUTINE - IT MAY BE BETTER TO *
C RESTART WITH IND = 1 OR 2. (COMPONENTS 10 TO 24 OF C ARE USED BY THE *
C PROGRAM - THE INFORMATION IS AVAILABLE TO THE USER, BUT SHOULD NOT *
C NORMALLY BE CHANGED BY HIM.) *
C *
C C(1) ERROR CONTROL INDICATOR - THE NORM OF THE LOCAL ERROR IS THE *
C MAX NORM OF THE WEIGHTED ERROR ESTIMATE VECTOR, THE *
C WEIGHTS BEING DETERMINED ACCORDING TO THE VALUE OF C(1) - *
C IF C(1)=1 THE WEIGHTS ARE 1 (ABSOLUTE ERROR CONTROL) *
C IF C(1)=2 THE WEIGHTS ARE 1/ABS(Y(K)) (RELATIVE ERROR *
C CONTROL) *
C IF C(1)=3 THE WEIGHTS ARE 1/MAX(ABS(C(2)),ABS(Y(K))) *
C (RELATIVE ERROR CONTROL, UNLESS ABS(Y(K)) IS LESS *
C THAN THE FLOOR VALUE, ABS(C(2)) ) *
C IF C(1)=4 THE WEIGHTS ARE 1/MAX(ABS(C(K+30)),ABS(Y(K))) *
C (HERE INDIVIDUAL FLOOR VALUES ARE USED) *
C IF C(1)=5 THE WEIGHTS ARE 1/ABS(C(K+30)) *
C FOR ALL OTHER VALUES OF C(1), INCLUDING C(1) = 0, THE *
C DEFAULT VALUES OF THE WEIGHTS ARE TAKEN TO BE *
C 1/MAX(1,ABS(Y(K))), AS MENTIONED EARLIER *
C (IN THE TWO CASES C(1) = 4 OR 5 THE USER MUST DECLARE THE *
C DIMENSION OF C TO BE AT LEAST N+30 AND MUST INITIALIZE THE *
C COMPONENTS C(31), C(32), ..., C(N+30).) *
C *
C C(2) FLOOR VALUE - USED WHEN THE INDICATOR C(1) HAS THE VALUE 3 *
C *
C C(3) HMIN SPECIFICATION - IF NOT ZERO, THE SUBROUTINE CHOOSES HMIN *
C TO BE ABS(C(3)) - OTHERWISE IT USES THE DEFAULT VALUE *
C 10*MAX(DWARF,RREB*MAX(WEIGHTED NORM Y/TOL,ABS(X))), *
C WHERE DWARF IS A VERY SMALL POSITIVE MACHINE NUMBER AND *
C RREB IS THE RELATIVE ROUNDOFF ERROR BOUND *
C *
C C(4) HSTART SPECIFICATION - IF NOT ZERO, THE SUBROUTINE WILL USE *
C AN INITIAL HMAG EQUAL TO ABS(C(4)), EXCEPT OF COURSE FOR *
C THE RESTRICTIONS IMPOSED BY HMIN AND HMAX - OTHERWISE IT *
C USES THE DEFAULT VALUE OF HMAX*(TOL)**(1/6) *
C *
C C(5) SCALE SPECIFICATION - THIS IS INTENDED TO BE A MEASURE OF THE *
C SCALE OF THE PROBLEM - LARGER VALUES OF SCALE TEND TO MAKE *
C THE METHOD MORE RELIABLE, FIRST BY POSSIBLY RESTRICTING *
C HMAX (AS DESCRIBED BELOW) AND SECOND, BY TIGHTENING THE *
C ACCEPTANCE REQUIREMENT - IF C(5) IS ZERO, A DEFAULT VALUE *
C OF 1 IS USED. FOR LINEAR HOMOGENEOUS PROBLEMS WITH *
C CONSTANT COEFFICIENTS, AN APPROPRIATE VALUE FOR SCALE IS A *
C NORM OF THE ASSOCIATED MATRIX. FOR OTHER PROBLEMS, AN *
C APPROXIMATION TO AN AVERAGE VALUE OF A NORM OF THE *
C JACOBIAN ALONG THE TRAJECTORY MAY BE APPROPRIATE *
C *
C C(6) HMAX SPECIFICATION - FOUR CASES ARE POSSIBLE *
C IF C(6).NE.0 AND C(5).NE.0, HMAX IS TAKEN TO BE *
C MIN(ABS(C(6)),2/ABS(C(5))) *
C IF C(6).NE.0 AND C(5).EQ.0, HMAX IS TAKEN TO BE ABS(C(6)) *
C IF C(6).EQ.0 AND C(5).NE.0, HMAX IS TAKEN TO BE *
C 2/ABS(C(5)) *
C IF C(6).EQ.0 AND C(5).EQ.0, HMAX IS GIVEN A DEFAULT VALUE *
C OF 2 *
C *
C C(7) MAXIMUM NUMBER OF FUNCTION EVALUATIONS - IF NOT ZERO, AN *
C ERROR RETURN WITH IND = -1 WILL BE CAUSED WHEN THE NUMBER *
C OF FUNCTION EVALUATIONS EXCEEDS ABS(C(7)) *
C *
C C(8) INTERRUPT NUMBER 1 - IF NOT ZERO, THE SUBROUTINE WILL *
C INTERRUPT THE CALCULATIONS AFTER IT HAS CHOSEN ITS *
C PRELIMINARY VALUE OF HMAG, AND JUST BEFORE CHOOSING HTRIAL *
C AND XTRIAL IN PREPARATION FOR TAKING A STEP (HTRIAL MAY *
C DIFFER FROM HMAG IN SIGN, AND MAY REQUIRE ADJUSTMENT IF *
C XEND IS NEAR) - THE SUBROUTINE RETURNS WITH IND = 4, AND *
C WILL RESUME CALCULATION AT THE POINT OF INTERRUPTION IF *
C RE-ENTERED WITH IND = 4 *
C *
C C(9) INTERRUPT NUMBER 2 - IF NOT ZERO, THE SUBROUTINE WILL *
C INTERRUPT THE CALCULATIONS IMMEDIATELY AFTER IT HAS *
C DECIDED WHETHER OR NOT TO ACCEPT THE RESULT OF THE MOST *
C RECENT TRIAL STEP, WITH IND = 5 IF IT PLANS TO ACCEPT, OR *
C IND = 6 IF IT PLANS TO REJECT - Y(*) IS THE PREVIOUSLY *
C ACCEPTED RESULT, WHILE W(*,9) IS THE NEWLY COMPUTED TRIAL *
C VALUE, AND W(*,2) IS THE UNWEIGHTED ERROR ESTIMATE VECTOR. *
C THE SUBROUTINE WILL RESUME CALCULATIONS AT THE POINT OF *
C INTERRUPTION ON RE-ENTRY WITH IND = 5 OR 6. (THE USER MAY *
C CHANGE IND IN THIS CASE IF HE WISHES, FOR EXAMPLE TO FORCE *
C ACCEPTANCE OF A STEP THAT WOULD OTHERWISE BE REJECTED, OR *
C VICE VERSA. HE CAN ALSO RESTART WITH IND = 1 OR 2.) *
C *
C***********************************************************************
C *
C SUMMARY OF THE COMPONENTS OF THE COMMUNICATIONS VECTOR *
C *
C PRESCRIBED AT THE OPTION DETERMINED BY THE PROGRAM *
C OF THE USER *
C *
C C(10) RREB(REL ROUNDOFF ERR BND) *
C C(1) ERROR CONTROL INDICATOR C(11) DWARF (VERY SMALL MACH NO) *
C C(2) FLOOR VALUE C(12) WEIGHTED NORM Y *
C C(3) HMIN SPECIFICATION C(13) HMIN *
C C(4) HSTART SPECIFICATION C(14) HMAG *
C C(5) SCALE SPECIFICATION C(15) SCALE *
C C(6) HMAX SPECIFICATION C(16) HMAX *
C C(7) MAX NO OF FCN EVALS C(17) XTRIAL *
C C(8) INTERRUPT NO 1 C(18) HTRIAL *
C C(9) INTERRUPT NO 2 C(19) EST *
C C(20) PREVIOUS XEND *
C C(21) FLAG FOR XEND *
C C(22) NO OF SUCCESSFUL STEPS *
C C(23) NO OF SUCCESSIVE FAILURES *
C C(24) NO OF FCN EVALS *
C *
C IF C(1) = 4 OR 5, C(31), C(32), ... C(N+30) ARE FLOOR VALUES *
C *
C RREB and DWARF are machine dependent constants currently set so *
C that they should be appropriate for most machines. However, it may *
C be appropriate to change them when this program is installed on a *
C new machine. K.R.J. 3 Oct 1991. *
C *
C***********************************************************************
C *
C AN OVERVIEW OF THE PROGRAM *
C *
C BEGIN INITIALIZATION, PARAMETER CHECKING, INTERRUPT RE-ENTRIES *
C ......ABORT IF IND OUT OF RANGE 1 TO 6 *
C . CASES - INITIAL ENTRY, NORMAL RE-ENTRY, INTERRUPT RE-ENTRIES *
C . CASE 1 - INITIAL ENTRY (IND .EQ. 1 OR 2) *
C V........ABORT IF N.GT.NW OR TOL.LE.0 *
C . IF INITIAL ENTRY WITHOUT OPTIONS (IND .EQ. 1) *
C . SET C(1) TO C(9) EQUAL TO ZERO *
C . ELSE INITIAL ENTRY WITH OPTIONS (IND .EQ. 2) *
C . MAKE C(1) TO C(9) NON-NEGATIVE *
C . MAKE FLOOR VALUES NON-NEGATIVE IF THEY ARE TO BE USED *
C . END IF *
C . INITIALIZE RREB, DWARF, PREV XEND, FLAG, COUNTS *
C . CASE 2 - NORMAL RE-ENTRY (IND .EQ. 3) *
C .........ABORT IF XEND REACHED, AND EITHER X CHANGED OR XEND NOT *
C . RE-INITIALIZE FLAG *
C . CASE 3 - RE-ENTRY FOLLOWING AN INTERRUPT (IND .EQ. 4 TO 6) *
C V TRANSFER CONTROL TO THE APPROPRIATE RE-ENTRY POINT....... *
C . END CASES . *
C . END INITIALIZATION, ETC. . *
C . V *
C . LOOP THROUGH THE FOLLOWING 4 STAGES, ONCE FOR EACH TRIAL STEP . *
C . STAGE 1 - PREPARE . *
C***********ERROR RETURN (WITH IND=-1) IF NO OF FCN EVALS TOO GREAT . *
C . CALC SLOPE (ADDING 1 TO NO OF FCN EVALS) IF IND .NE. 6 . *
C . CALC HMIN, SCALE, HMAX . *
C***********ERROR RETURN (WITH IND=-2) IF HMIN .GT. HMAX . *
C . CALC PRELIMINARY HMAG . *
C***********INTERRUPT NO 1 (WITH IND=4) IF REQUESTED.......RE-ENTRY.V *
C . CALC HMAG, XTRIAL AND HTRIAL . *
C . END STAGE 1 . *
C V STAGE 2 - CALC YTRIAL (ADDING 7 TO NO OF FCN EVALS) . *
C . STAGE 3 - CALC THE ERROR ESTIMATE . *
C . STAGE 4 - MAKE DECISIONS . *
C . SET IND=5 IF STEP ACCEPTABLE, ELSE SET IND=6 . *
C***********INTERRUPT NO 2 IF REQUESTED....................RE-ENTRY.V *
C . IF STEP ACCEPTED (IND .EQ. 5) *
C . UPDATE X, Y FROM XTRIAL, YTRIAL *
C . ADD 1 TO NO OF SUCCESSFUL STEPS *
C . SET NO OF SUCCESSIVE FAILURES TO ZERO *
C**************RETURN(WITH IND=3, XEND SAVED, FLAG SET) IF X .EQ. XEND *
C . ELSE STEP NOT ACCEPTED (IND .EQ. 6) *
C . ADD 1 TO NO OF SUCCESSIVE FAILURES *
C**************ERROR RETURN (WITH IND=-3) IF HMAG .LE. HMIN *
C . END IF *
C . END STAGE 4 *
C . END LOOP *
C . *
C BEGIN ABORT ACTION *
C OUTPUT APPROPRIATE MESSAGE ABOUT STOPPING THE CALCULATIONS, *
C ALONG WITH VALUES OF IND, N, NW, TOL, HMIN, HMAX, X, XEND, *
C PREVIOUS XEND, NO OF SUCCESSFUL STEPS, NO OF SUCCESSIVE *
C FAILURES, NO OF FCN EVALS, AND THE COMPONENTS OF Y *
C STOP *
C END ABORT ACTION *
C *
C***********************************************************************
C
C ******************************************************************
C * BEGIN INITIALIZATION, PARAMETER CHECKING, INTERRUPT RE-ENTRIES *
C ******************************************************************
C
C ......ABORT IF IND OUT OF RANGE 1 TO 6
IF (IND.LT.1 .OR. IND.GT.6) GO TO 500
C
C CASES - INITIAL ENTRY, NORMAL RE-ENTRY, INTERRUPT RE-ENTRIES
GO TO (5, 5, 45, 1111, 2222, 2222), IND
C CASE 1 - INITIAL ENTRY (IND .EQ. 1 OR 2)
C .........ABORT IF N.GT.NW OR TOL.LE.0
5 IF (N.GT.NW .OR. TOL.LE.0.D0) GO TO 500
IF (IND.EQ. 2) GO TO 15
C INITIAL ENTRY WITHOUT OPTIONS (IND .EQ. 1)
C SET C(1) TO C(9) EQUAL TO 0
DO 10 K = 1, 9
C(K) = 0.D0
10 CONTINUE
GO TO 35
15 CONTINUE
C INITIAL ENTRY WITH OPTIONS (IND .EQ. 2)
C MAKE C(1) TO C(9) NON-NEGATIVE
DO 20 K = 1, 9
C(K) = DABS(C(K))
20 CONTINUE
C MAKE FLOOR VALUES NON-NEGATIVE IF THEY ARE TO BE USED
IF (C(1).NE.4.D0 .AND. C(1).NE.5.D0) GO TO 30
DO 25 K = 1, N
C(K+30) = DABS(C(K+30))
25 CONTINUE
30 CONTINUE
35 CONTINUE
C INITIALIZE RREB, DWARF, PREV XEND, FLAG, COUNTS
C(10) = 2.D0**(-56)
C(11) = 1.D-35
C SET PREVIOUS XEND INITIALLY TO INITIAL VALUE OF X
C(20) = X
DO 40 K = 21, 24
C(K) = 0.D0
40 CONTINUE
GO TO 50
C CASE 2 - NORMAL RE-ENTRY (IND .EQ. 3)
C .........ABORT IF XEND REACHED, AND EITHER X CHANGED OR XEND NOT
45 IF (C(21).NE.0.D0 .AND.
+ (X.NE.C(20) .OR. XEND.EQ.C(20))) GO TO 500
C RE-INITIALIZE FLAG
C(21) = 0.D0
GO TO 50
C CASE 3 - RE-ENTRY FOLLOWING AN INTERRUPT (IND .EQ. 4 TO 6)
C TRANSFER CONTROL TO THE APPROPRIATE RE-ENTRY POINT..........
C THIS HAS ALREADY BEEN HANDLED BY THE COMPUTED GO TO .
C END CASES V
50 CONTINUE
C
C END INITIALIZATION, ETC.
C
C ******************************************************************
C * LOOP THROUGH THE FOLLOWING 4 STAGES, ONCE FOR EACH TRIAL STEP *
C * UNTIL THE OCCURRENCE OF ONE OF THE FOLLOWING *
C * (A) THE NORMAL RETURN (WITH IND .EQ. 3) ON REACHING XEND IN *
C * STAGE 4 *
C * (B) AN ERROR RETURN (WITH IND .LT. 0) IN STAGE 1 OR STAGE 4 *
C * (C) AN INTERRUPT RETURN (WITH IND .EQ. 4, 5 OR 6), IF *
C * REQUESTED, IN STAGE 1 OR STAGE 4 *
C ******************************************************************
C
99999 CONTINUE
C
C ***************************************************************
C * STAGE 1 - PREPARE - DO CALCULATIONS OF HMIN, HMAX, ETC., *
C * AND SOME PARAMETER CHECKING, AND END UP WITH SUITABLE *
C * VALUES OF HMAG, XTRIAL AND HTRIAL IN PREPARATION FOR TAKING *
C * AN INTEGRATION STEP. *
C ***************************************************************
C
C***********ERROR RETURN (WITH IND=-1) IF NO OF FCN EVALS TOO GREAT
IF (C(7).EQ.0.D0 .OR. C(24).LT.C(7)) GO TO 100
IND = -1
RETURN
100 CONTINUE
C
C CALCULATE SLOPE (ADDING 1 TO NO OF FCN EVALS) IF IND .NE. 6
IF (IND .EQ. 6) GO TO 105
CALL FCN(N, X, Y, W(1,1))
C(24) = C(24) + 1.D0
105 CONTINUE
C
C CALCULATE HMIN - USE DEFAULT UNLESS VALUE PRESCRIBED
C(13) = C(3)
IF (C(3) .NE. 0.D0) GO TO 165
C CALCULATE DEFAULT VALUE OF HMIN
C FIRST CALCULATE WEIGHTED NORM Y - C(12) - AS SPECIFIED
C BY THE ERROR CONTROL INDICATOR C(1)
TEMP = 0.D0
IF (C(1) .NE. 1.D0) GO TO 115
C ABSOLUTE ERROR CONTROL - WEIGHTS ARE 1
DO 110 K = 1, N
TEMP = DMAX1(TEMP, DABS(Y(K)))
110 CONTINUE
C(12) = TEMP
GO TO 160
115 IF (C(1) .NE. 2.D0) GO TO 120
C RELATIVE ERROR CONTROL - WEIGHTS ARE 1/DABS(Y(K)) SO
C WEIGHTED NORM Y IS 1
C(12) = 1.D0
GO TO 160
120 IF (C(1) .NE. 3.D0) GO TO 130
C WEIGHTS ARE 1/MAX(C(2),ABS(Y(K)))
DO 125 K = 1, N
TEMP = DMAX1(TEMP, DABS(Y(K))/C(2))
125 CONTINUE
C(12) = DMIN1(TEMP, 1.D0)
GO TO 160
130 IF (C(1) .NE. 4.D0) GO TO 140
C WEIGHTS ARE 1/MAX(C(K+30),ABS(Y(K)))
DO 135 K = 1, N
TEMP = DMAX1(TEMP, DABS(Y(K))/C(K+30))
135 CONTINUE
C(12) = DMIN1(TEMP, 1.D0)
GO TO 160
140 IF (C(1) .NE. 5.D0) GO TO 150
C WEIGHTS ARE 1/C(K+30)
DO 145 K = 1, N
TEMP = DMAX1(TEMP, DABS(Y(K))/C(K+30))
145 CONTINUE
C(12) = TEMP
GO TO 160
150 CONTINUE
C DEFAULT CASE - WEIGHTS ARE 1/MAX(1,ABS(Y(K)))
DO 155 K = 1, N
TEMP = DMAX1(TEMP, DABS(Y(K)))
155 CONTINUE
C(12) = DMIN1(TEMP, 1.D0)
160 CONTINUE
C(13) = 10.D0*DMAX1(C(11),C(10)*DMAX1(C(12)/TOL,DABS(X)))
165 CONTINUE
C
C CALCULATE SCALE - USE DEFAULT UNLESS VALUE PRESCRIBED
C(15) = C(5)
IF (C(5) .EQ. 0.D0) C(15) = 1.D0
C
C CALCULATE HMAX - CONSIDER 4 CASES
C CASE 1 BOTH HMAX AND SCALE PRESCRIBED
IF (C(6).NE.0.D0 .AND. C(5).NE.0.D0)
+ C(16) = DMIN1(C(6), 2.D0/C(5))
C CASE 2 - HMAX PRESCRIBED, BUT SCALE NOT
IF (C(6).NE.0.D0 .AND. C(5).EQ.0.D0) C(16) = C(6)
C CASE 3 - HMAX NOT PRESCRIBED, BUT SCALE IS
IF (C(6).EQ.0.D0 .AND. C(5).NE.0.D0) C(16) = 2.D0/C(5)
C CASE 4 - NEITHER HMAX NOR SCALE IS PROVIDED
IF (C(6).EQ.0.D0 .AND. C(5).EQ.0.D0) C(16) = 2.D0
C
C***********ERROR RETURN (WITH IND=-2) IF HMIN .GT. HMAX
IF (C(13) .LE. C(16)) GO TO 170
IND = -2
RETURN
170 CONTINUE
C
C CALCULATE PRELIMINARY HMAG - CONSIDER 3 CASES
IF (IND .GT. 2) GO TO 175
C CASE 1 - INITIAL ENTRY - USE PRESCRIBED VALUE OF HSTART, IF
C ANY, ELSE DEFAULT
C(14) = C(4)
IF (C(4) .EQ. 0.D0) C(14) = C(16)*TOL**(1./6.)
GO TO 185
175 IF (C(23) .GT. 1.D0) GO TO 180
C CASE 2 - AFTER A SUCCESSFUL STEP, OR AT MOST ONE FAILURE,
C USE MIN(2, .9*(TOL/EST)**(1/6))*HMAG, BUT AVOID POSSIBLE
C OVERFLOW. THEN AVOID REDUCTION BY MORE THAN HALF.
TEMP = 2.D0*C(14)
IF (TOL .LT. (2.D0/.9D0)**6*C(19))
+ TEMP = .9D0*(TOL/C(19))**(1./6.)*C(14)
C(14) = DMAX1(TEMP, .5D0*C(14))
GO TO 185
180 CONTINUE
C CASE 3 - AFTER TWO OR MORE SUCCESSIVE FAILURES
C(14) = .5D0*C(14)
185 CONTINUE
C
C CHECK AGAINST HMAX
C(14) = DMIN1(C(14), C(16))
C
C CHECK AGAINST HMIN
C(14) = DMAX1(C(14), C(13))
C
C***********INTERRUPT NO 1 (WITH IND=4) IF REQUESTED
IF (C(8) .EQ. 0.D0) GO TO 1111
IND = 4
RETURN
C RESUME HERE ON RE-ENTRY WITH IND .EQ. 4 ........RE-ENTRY..
1111 CONTINUE
C
C CALCULATE HMAG, XTRIAL - DEPENDING ON PRELIMINARY HMAG, XEND
IF (C(14) .GE. DABS(XEND - X)) GO TO 190
C DO NOT STEP MORE THAN HALF WAY TO XEND
C(14) = DMIN1(C(14), .5D0*DABS(XEND - X))
C(17) = X + DSIGN(C(14), XEND - X)
GO TO 195
190 CONTINUE
C HIT XEND EXACTLY
C(14) = DABS(XEND - X)
C(17) = XEND
195 CONTINUE
C
C CALCULATE HTRIAL
C(18) = C(17) - X
C
C END STAGE 1
C
C ***************************************************************
C * STAGE 2 - CALCULATE YTRIAL (ADDING 7 TO NO OF FCN EVALS). *
C * W(*,2), ... W(*,8) HOLD INTERMEDIATE RESULTS NEEDED IN *
C * STAGE 3. W(*,9) IS TEMPORARY STORAGE UNTIL FINALLY IT HOLDS *
C * YTRIAL. *
C ***************************************************************
C
TEMP = C(18)/1398169080000.D0
C
DO 200 K = 1, N
W(K,9) = Y(K) + TEMP*W(K,1)*233028180000.D0
200 CONTINUE
CALL FCN(N, X + C(18)/6.D0, W(1,9), W(1,2))
C
DO 205 K = 1, N
W(K,9) = Y(K) + TEMP*( W(K,1)*74569017600.D0
+ + W(K,2)*298276070400.D0 )
205 CONTINUE
CALL FCN(N, X + C(18)*(4.D0/15.D0), W(1,9), W(1,3))
C
DO 210 K = 1, N
W(K,9) = Y(K) + TEMP*( W(K,1)*1165140900000.D0
+ - W(K,2)*3728450880000.D0
+ + W(K,3)*3495422700000.D0 )
210 CONTINUE
CALL FCN(N, X + C(18)*(2.D0/3.D0), W(1,9), W(1,4))
C
DO 215 K = 1, N
W(K,9) = Y(K) + TEMP*( - W(K,1)*3604654659375.D0
+ + W(K,2)*12816549900000.D0
+ - W(K,3)*9284716546875.D0
+ + W(K,4)*1237962206250.D0 )
215 CONTINUE
CALL FCN(N, X + C(18)*(5.D0/6.D0), W(1,9), W(1,5))
C
DO 220 K = 1, N
W(K,9) = Y(K) + TEMP*( W(K,1)*3355605792000.D0
+ - W(K,2)*11185352640000.D0
+ + W(K,3)*9172628850000.D0
+ - W(K,4)*427218330000.D0
+ + W(K,5)*482505408000.D0 )
220 CONTINUE
CALL FCN(N, X + C(18), W(1,9), W(1,6))
C
DO 225 K = 1, N
W(K,9) = Y(K) + TEMP*( - W(K,1)*770204740536.D0
+ + W(K,2)*2311639545600.D0
+ - W(K,3)*1322092233000.D0
+ - W(K,4)*453006781920.D0
+ + W(K,5)*326875481856.D0 )
225 CONTINUE
CALL FCN(N, X + C(18)/15.D0, W(1,9), W(1,7))
C
DO 230 K = 1, N
W(K,9) = Y(K) + TEMP*( W(K,1)*2845924389000.D0
+ - W(K,2)*9754668000000.D0
+ + W(K,3)*7897110375000.D0
+ - W(K,4)*192082660000.D0
+ + W(K,5)*400298976000.D0
+ + W(K,7)*201586000000.D0 )
230 CONTINUE
CALL FCN(N, X + C(18), W(1,9), W(1,8))
C
C CALCULATE YTRIAL, THE EXTRAPOLATED APPROXIMATION AND STORE
C IN W(*,9)
DO 235 K = 1, N
W(K,9) = Y(K) + TEMP*( W(K,1)*104862681000.D0
+ + W(K,3)*545186250000.D0
+ + W(K,4)*446637345000.D0
+ + W(K,5)*188806464000.D0
+ + W(K,7)*15076875000.D0
+ + W(K,8)*97599465000.D0 )
235 CONTINUE
C
C ADD 7 TO THE NO OF FCN EVALS
C(24) = C(24) + 7.D0
C
C END STAGE 2
C
C ***************************************************************
C * STAGE 3 - CALCULATE THE ERROR ESTIMATE EST. FIRST CALCULATE *
C * THE UNWEIGHTED ABSOLUTE ERROR ESTIMATE VECTOR (PER UNIT *
C * STEP) FOR THE UNEXTRAPOLATED APPROXIMATION AND STORE IT IN *
C * W(*,2). THEN CALCULATE THE WEIGHTED MAX NORM OF W(*,2) AS *
C * SPECIFIED BY THE ERROR CONTROL INDICATOR C(1). FINALLY, *
C * MODIFY THIS RESULT TO PRODUCE EST, THE ERROR ESTIMATE (PER *
C * UNIT STEP) FOR THE EXTRAPOLATED APPROXIMATION YTRIAL. *
C ***************************************************************
C
C CALCULATE THE UNWEIGHTED ABSOLUTE ERROR ESTIMATE VECTOR
DO 300 K = 1, N
W(K,2) = ( W(K,1)*8738556750.D0
+ + W(K,3)*9735468750.D0
+ - W(K,4)*9709507500.D0
+ + W(K,5)*8582112000.D0
+ + W(K,6)*95329710000.D0
+ - W(K,7)*15076875000.D0
+ - W(K,8)*97599465000.D0)/1398169080000.D0
300 CONTINUE
C
C CALCULATE THE WEIGHTED MAX NORM OF W(*,2) AS SPECIFIED BY
C THE ERROR CONTROL INDICATOR C(1)
TEMP = 0.D0
IF (C(1) .NE. 1.D0) GO TO 310
C ABSOLUTE ERROR CONTROL
DO 305 K = 1, N
TEMP = DMAX1(TEMP,DABS(W(K,2)))
305 CONTINUE
GO TO 360
310 IF (C(1) .NE. 2.D0) GO TO 320
C RELATIVE ERROR CONTROL
DO 315 K = 1, N
TEMP = DMAX1(TEMP, DABS(W(K,2)/Y(K)))
315 CONTINUE
GO TO 360
320 IF (C(1) .NE. 3.D0) GO TO 330
C WEIGHTS ARE 1/MAX(C(2),ABS(Y(K)))
DO 325 K = 1, N
TEMP = DMAX1(TEMP, DABS(W(K,2))
+ / DMAX1(C(2), DABS(Y(K))) )
325 CONTINUE
GO TO 360
330 IF (C(1) .NE. 4.D0) GO TO 340
C WEIGHTS ARE 1/MAX(C(K+30),ABS(Y(K)))
DO 335 K = 1, N
TEMP = DMAX1(TEMP, DABS(W(K,2))
+ / DMAX1(C(K+30), DABS(Y(K))) )
335 CONTINUE
GO TO 360
340 IF (C(1) .NE. 5.D0) GO TO 350
C WEIGHTS ARE 1/C(K+30)
DO 345 K = 1, N
TEMP = DMAX1(TEMP, DABS(W(K,2)/C(K+30)))
345 CONTINUE
GO TO 360
350 CONTINUE
C DEFAULT CASE - WEIGHTS ARE 1/MAX(1,ABS(Y(K)))
DO 355 K = 1, N
TEMP = DMAX1(TEMP, DABS(W(K,2))
+ / DMAX1(1.D0, DABS(Y(K))) )
355 CONTINUE
360 CONTINUE
C
C CALCULATE EST - (THE WEIGHTED MAX NORM OF W(*,2))*HMAG*SCALE
C - EST IS INTENDED TO BE A MEASURE OF THE ERROR PER UNIT
C STEP IN YTRIAL
C(19) = TEMP*C(14)*C(15)
C
C END STAGE 3
C
C ***************************************************************
C * STAGE 4 - MAKE DECISIONS. *
C ***************************************************************
C
C SET IND=5 IF STEP ACCEPTABLE, ELSE SET IND=6
IND = 5
IF (C(19) .GT. TOL) IND = 6
C
C***********INTERRUPT NO 2 IF REQUESTED
IF (C(9) .EQ. 0.D0) GO TO 2222
RETURN
C RESUME HERE ON RE-ENTRY WITH IND .EQ. 5 OR 6 ...RE-ENTRY..
2222 CONTINUE
C
IF (IND .EQ. 6) GO TO 410
C STEP ACCEPTED (IND .EQ. 5), SO UPDATE X, Y FROM XTRIAL,
C YTRIAL, ADD 1 TO THE NO OF SUCCESSFUL STEPS, AND SET
C THE NO OF SUCCESSIVE FAILURES TO ZERO
X = C(17)
DO 400 K = 1, N
Y(K) = W(K,9)
400 CONTINUE
C(22) = C(22) + 1.D0
C(23) = 0.D0
C**************RETURN(WITH IND=3, XEND SAVED, FLAG SET) IF X .EQ. XEND
IF (X .NE. XEND) GO TO 405
IND = 3
C(20) = XEND
C(21) = 1.D0
RETURN
405 CONTINUE
GO TO 420
410 CONTINUE
C STEP NOT ACCEPTED (IND .EQ. 6), SO ADD 1 TO THE NO OF
C SUCCESSIVE FAILURES
C(23) = C(23) + 1.D0
C**************ERROR RETURN (WITH IND=-3) IF HMAG .LE. HMIN
IF (C(14) .GT. C(13)) GO TO 415
IND = -3
RETURN
415 CONTINUE
420 CONTINUE
C
C END STAGE 4
C
GO TO 99999
C END LOOP
C
C BEGIN ABORT ACTION
500 CONTINUE
C
WRITE(6,505) IND, TOL, X, N, C(13), XEND, NW, C(16), C(20),
+ C(22), C(23), C(24), (Y(K), K = 1, N)
505 FORMAT( /// 1H0, 58HCOMPUTATION STOPPED IN DVERK WITH THE FOLLOWIN
+ G VALUES-
+ / 1H0, 5HIND =, I4, 5X, 6HTOL =, 1PD13.6, 5X, 11HX =,
+ 1PD22.15
+ / 1H , 5HN =, I4, 5X, 6HHMIN =, 1PD13.6, 5X, 11HXEND =,
+ 1PD22.15
+ / 1H , 5HNW =, I4, 5X, 6HHMAX =, 1PD13.6, 5X, 11HPREV XEND =,
+ 1PD22.15
+ / 1H0, 14X, 27HNO OF SUCCESSFUL STEPS =, 0PF8.0
+ / 1H , 14X, 27HNO OF SUCCESSIVE FAILURES =, 0PF8.0
+ / 1H , 14X, 27HNO OF FUNCTION EVALS =, 0PF8.0
+ / 1H0, 23HTHE COMPONENTS OF Y ARE
+ // (1H , 1P5D24.15) )
C
STOP
C
C END ABORT ACTION
C
END
C
C
subroutine fcn(neq, t, y, yprime)
implicit none
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
include 'fcn.inc'
c-----------------------------------------------------------------------
c Arguments of subroutine and indices
c Note: dimension of y is nst*(nst+3) (nst = dims(4))
c dimension of yprime is neq = dims(1)
c-----------------------------------------------------------------------
integer ic, jc
integer minst, maxst, nst
real*8 t
integer neq
real*8 y, yprime,TOTIME,HBAR, V,W,TOCYC,THETA
dimension neq(4)
dimension y(1), yprime(1)
c-----------------------------------------------------------------------
c Temporaries for dimension of system
c-----------------------------------------------------------------------
minst = neq(2)
maxst = neq(3)
nst = neq(4)
c-----------------------------------------------------------------------
c Define derivatives
c Note: y(ic) = a(ic+maxst-nst)
c y(ic+nst) = b(ic+maxst-nst)
c y(ic+2*nst) = E(ic+maxst-nst)
c y(nst*ic+jc+2*nst) = x(jc+maxst-nst,ic+maxst-nst)
c
c Note: we will replace x(j,i) with x(i,j) in the calculations
c below for better vectorization. We can do this as long as
c the x(j,i) are real.
c-----------------------------------------------------------------------
if (t .lt. totime) then
do ic = 1, nst
yprime(ic) = y(ic+2*nst)*y(ic+nst)/hbar
yprime(ic+nst) = - y(ic+2*nst)*y(ic)/hbar
do jc = 1, nst
yprime(ic) = yprime(ic) -
& V*(sin(w*t/(4.0d0*tocyc))**2)*
& cos(w*t+theta)*
& y(nst*ic+jc+2*nst)*y(jc+nst)/hbar
yprime(ic+nst) = yprime(ic+nst)
& +V*(sin(w*t/(4.0d0*tocyc))**2)*
& cos(w*t+theta)*
& y(nst*ic+jc+2*nst)*y(jc)/hbar
end do
end do
else
do ic = 1, nst
yprime(ic) = y(ic+2*nst)*y(ic+nst)/hbar
yprime(ic+nst) = - y(ic+2*nst)*y(ic)/hbar
do jc = 1, nst
yprime(ic) = yprime(ic) -
& V*cos(w*t+theta)*y(nst*ic+jc+2*nst)
& *y(jc+nst)/hbar
yprime(ic+nst) = yprime(ic+nst) +
& V*cos(w*t+theta)*y(nst*ic+jc+2*nst)
& *y(jc)/hbar
end do
end do
end if
return
end