×
INTELLIGENT WORK FORUMS
FOR ENGINEERING PROFESSIONALS

Log In

Come Join Us!

Are you an
Engineering professional?
Join Eng-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!
  • Students Click Here

*Eng-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Students Click Here

Run-time error in engine heat release rate program

Run-time error in engine heat release rate program

Run-time error in engine heat release rate program

(OP)
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
      

RE: Run-time error in engine heat release rate program

Chris

I have compiled your code using Salford FTN95 and get your same error message.

However, there is more wrong with this source than just missing routines, my compiler diagnostics are as follows:-

    NO ERRORS  [<main program>
0141) YPRIME(3) = -C1/CPU/(1. + CO)*(Y(3) - TW) + VU/CPU*DLVLTU*
0142) &                               YPRIME(1)/10.
WARNING - Variable 'CO' has been used without being given an initial value
0135) E = (1. - X)*(VU**2/CPU/Y(3)*DLVLTU**2 + 10.*VU/Y(1)*DLVLPU)
WARNING - Variable 'DLVLPU' has been used without being given an initial
    value
0129) B = C1*(VB/CPB*DLVLTB*C0*(1. - TW/Y(2)) + VU/CPU*DLVLTU*
0130) &               (1. - C0)*(1.- TW/Y(3)))
WARNING - Variable 'DLVLTB' has been used without being given an initial
    value
    NO ERRORS, 3 WARNINGS  [<RATES>
    NO ERRORS  [<AUXLRY>
    NO ERRORS  [<DATA@>
    NO ERRORS  [<TINITL>
0370) A(3,2) = 2.0*D14 + D24 + 2.0 + D84 + D04 + D104 - D(2)*D14
WARNING - Variable 'D04' has been used without being given an initial value
    NO ERRORS, 1 WARNING  [<ECP>
0550) LOGICAL RICH, LEAN
WARNING - Variable 'LEAN' has been given a value but never used
    NO ERRORS, 1 WARNING  [<FARG>
0681) DOUBLE PRECISION X, Y(N), XEND, TOL, C(*), W(NW,9), TEMP
WARNING - Variable 'TEMP' has been declared but not used
0680) INTEGER N, IND, NW, K
WARNING - Variable 'K' has been declared but not used
    NO ERRORS, 2 WARNINGS  [<DVERK>
0683) external fcn
*** Missing END SUBROUTINE
WARNING the following symbols are missing:
DVERK
T     
LEQIF
AMIN  
AMAX  



The use of variables that have not been assigned a value is a show stopper, until these are fixed it is not worth finding the missing routines!

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Eng-Tips Forums free from inappropriate posts.
The Eng-Tips staff will check this out and take appropriate action.

Reply To This Thread

Posting in the Eng-Tips forums is a member-only feature.

Click Here to join Eng-Tips and talk with other members! Already a Member? Login


Resources

White Paper - How ESI is Helping Move New Medical Device Product to Market Quicker & More Cost Effic
Early Supplier Involvement has long been a strategy employed by manufacturers to produce innovative products. Now, it almost seems like a necessity. Because decisions made in the design phase can positively affect product quality and costs, this can help add value to OEM bottom lines. This white paper will discuss many facets of ESI, including why it’s so valuable today, what challenges limit the benefits of ESI, how cost is impacted, and more. Download Now
White Paper - Moving to a Driverless Future
This white paper describes what we see as the best practices to support a sustainable engineering process for autonomous vehicle design. It exposes how to use simulation and testing in common frameworks to enable design exploration, verification and validation for the development of autonomous cars at a system, software and full-vehicle level to drive a mature product development process for automated driving. Download Now
Research Report - How Engineers are Using Remote Access
Remote access enables engineers to work from anywhere provided they have an internet connection. We surveyed our audience of engineers, designers and product managers to learn how they use remote access within their organizations. We wanted to know which industries have adopted remote access, which software they are using, and what features matter most. Download Now

Close Box

Join Eng-Tips® Today!

Join your peers on the Internet's largest technical engineering professional community.
It's easy to join and it's free.

Here's Why Members Love Eng-Tips Forums:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close