SUBROUTINE DEFENG IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) C INPUT ENGINE RELATED DATA - NAMELIST ENGDIN AND THE ENGINE DECK LOGICAL ENGOPT COMMON /ENGPRF/ CDT, CDP, VJET, CDTMAX, CDPMAX, VJMAX, ENGFN, 1 ENGSFC, OVEFF, SPECT, STMIN, 2 ARATIO, ARMAX, XACRS, XMCRS, IGENEN COMMON /CONTRL/ NICE, IANAL COMMON /UNITS / IU5, IU6, IU7, IU8, IU9, IU16, IU17, IU18 COMMON /REGEN / EDV(6), IREGEN COMMON /THRPLT/ PCODE(16), NPCODE, IPLTTH COMMON /NOXDAT/ OX(16,15,20), FFUEL(6), FNOX(6), RTNOX, TNOX, 1 GNOX, OFNOX, NOX COMMON /ENDTA / EXTFAC,FFFAC ,DFFAC ,TAKOFF,EMACH(20) , 1 ALT(15,20) ,FF(16,15,20) ,THR(16,15,20), 2 TXFUFL,NM ,NA ,NT ,IXTRAP,MAXCR C COMMON BLOCK ENDTA IS USED TO TRANSFER THE ENGINE DATA BETWEEN C SUBROUTINES. THE VARIABLES USED ARE C REAL VARIABLES C DFFAC = (0.0) FUEL FLOW SCALING CONSTANT TERM C EXTFAC = (1.0) SLOPE FACTOR FOR EXTRAPOLATING FUEL FLOWS C FFFAC = (0.0) FUEL FLOW SCALING LINEAR TERM C TAKOFF = TAKEOFF FUEL FLOW C TXFUFL = TAXI FUEL FLOW C REAL ARRAYS C ALT = ALTITUDES - MAX NUMBER IS 15 FOR EACH MACH NO. C EMACH = MACH NUMBERS - MAX NUMBER IS 20 C FF = FUEL FLOWS - ONE FOR EACH THRUST LEVEL C THR = THRUST LEVELS - MAX NUMBER IS 16 FOR EACH ALT FOR C EACH MACH NO. C INTEGER VARIABLES C IXTRAP = (1) PREVENTS IMPROVEMENT IN EXTRAPOLATED SFC'S C MAXCR = (2) MAX POWER SETTING USED FOR CRUISE C NA = NUMBER OF ALTITUDES C NM = NUMBER OF MACH NUMBERS C NT = NUMBER OF THRUST LEVELS COMMON /SCRTCH/ EM(1600), AL(1600), TH(1600), FL(1600), XO(1600), 1 K(1600), IPC(1600) COMMON / SYNT / ACCUX ,EPS ,RK ,GLM ,PEN ,XBJ ,AMULT , 1 EF ,DEP ,NDD ,JECT ,JVKC COMMON /CDFIL / CDFILE, A10, A9REF, A10REF, XNOZ, XNREF, RCRV, 1 NAB, NABREF COMMON / FFFF / FFFSUB, FFFSUP DIMENSION SFC(16) CHARACTER*80 EIFILE, TITLE, CDFILE SAVE ICHECK, IST NAMELIST /ENGDIN/ EXTFAC, FFFAC, DFFAC, IDLE, IXTRAP, IFILL, ALT, 1 MAXCR, BOOST, NGPRT, IGEO, IGENEN, EMACH, NONEG, 2 NPCODE, PCODE, EIFILE, NOX, FIDMIN, FIDMAX, 3 INSDRG, CDFILE, A10, A9REF, A10REF, XNOZ, XNREF, 4 NAB, NABREF, FFFSUB, FFFSUP, RCRV DATA BOOST /0./, IFILL /2/, NIP /1600/, IDLE, ISTOP, IGEO /3*0/, 1 REARTH /6367533./, GR /9.80665/, GNS /9.823693/, CM1 /.9985/, 2 OC2 /26.76566E-10/, NGPRT /1/, EIFILE /'ENGDEK'/, 3 FIDMIN /0.08/, FIDMAX /1.00/ DATA ICHECK, NONEG, INSDRG /3*0/, IST /1/ DATA FFPEN / 1.15 / C ZERO OUT ARRAYS DO 90 M = 1,20 DO 90 N = 1,15 DO 90 L = 1,16 THR(L,N,M) = 0.0 FF (L,N,M) = 0.0 OX (L,N,M) = 0.0 90 CONTINUE IFLAG = 0 C READ NAMELIST $ENGDIN UNLESS PARAMETRIC C VARIATION ON ENGINE DESIGN ONLY. IF ( IANAL .EQ. 4 ) THEN IGENEN = 2 ELSE IF ( IREGEN .GT. 0 ) THEN IGENEN = IGENEN + 10 GO TO 92 ENDIF IGENEN = 0 READ (IU5, NML=ENGDIN, ERR=9000, END=9000) IF ( MAXCR .EQ. 0 ) MAXCR = 2 WRITE(IU6, 340) NGPRT, IGENEN, EXTFAC, FFFSUB, FFFSUP, IDLE, 1 NONEG IF ( IDLE .GT. 0 ) WRITE(IU6, 342) FIDMIN, FIDMAX WRITE(IU6, 343) IXTRAP, IFILL, MAXCR, BOOST, DFFAC, FFFAC, NOX, 1 INSDRG IF ( INSDRG .GT. 0 ) WRITE(IU6, 344) CDFILE, NAB, NABREF, A10, 1 A9REF, A10REF, XNOZ, XNREF IF ( IGENEN .LE. 0 ) THEN IF ( IGENEN .EQ. 0 ) EIFILE = 'INPUT ' WRITE(IU6, 345) EIFILE GO TO 6 ENDIF ENDIF 92 CONTINUE C----------------------------------------------------------------------C C THIS BLOCK OF CODE AND A BLOCK IN ANALYS MAY BE COMMENTED OUT TO C DISCONNECT THE ENGINE CYCLE ANALYSIS MODULE C GENERATE ENGINE DECK cjmg IF ( ICHECK .EQ. 1 ) GO TO 5 cjmg ICHECK = 1 cjmg NA = 0 cjmg IEM = 0 cjmg IEA = 0 cjmg IF ( EMACH(1) .GE. 0. ) WRITE(IU6,2500) cjmg DO 3 M = 1,20 cjmg IF ( EMACH(M) .LT. 0. ) GO TO 4 cjmg IF ( M .GT. 1 .AND. EMACH(M) .GE. EMACH(M-1) ) IEM = 1 cjmg DO 1 N = 2,15 cjmg IF ( ALT(N,M) .LT. 0. ) GO TO 2 cjmg IF ( ALT(N,M) .GE. ALT(N-1,M) ) IEA = 1 cjmg 1 CONTINUE cjmg N = 16 cjmg 2 IF ( N - 1 .GT. NA ) NA = N - 1 cjmg WRITE(IU6,2510) EMACH(M), (ALT(I,M), I = 1,N-1) cjmg 3 CONTINUE cjmg M = 21 cjmg 4 NM = M - 1 cjmg IF ( IEM .EQ. 1 ) WRITE(IU6,2520) cjmg IF ( IEA .EQ. 1 ) WRITE(IU6,2530) cjmg IF ( IEM .EQ. 1 .OR. IEA .EQ. 1 ) STOP cjmg 5 ENGOPT = .FALSE. cjmg CALL GENENG ( DESSFC, ENGOPT ) cjmg IF ( JECT .NE. 0 ) RETURN cjmg IF ( IGENEN .EQ. 2 ) GO TO 335 cjmg GO TO 145 C----------------------------------------------------------------------C C READ ENGINE DECK 6 J = 1 I = 0 IU11 = IU5 IF ( IGENEN .EQ. 0 ) GO TO 7 IU11 = 11 OPEN (UNIT=IU11, FILE=EIFILE, STATUS='UNKNOWN', ERR=3140) 7 K(J) = J I = I + 1 READ (IU11, 350, ERR=3000, END=20) EM(J), AL(J), PC, TG, RD, 1 FL(J), XO(J), A9 IF ( EM(J) .GT. 5.0 ) GO TO 20 IF ( FL(J) .EQ. 0.0 ) GO TO 7 TH(J) = TG - RD C IF POWER CODES ARE TO BE USED, SET UP ARRAY IF ( NPCODE .LT. 2 ) GO TO 12 DO 8 NP = 1,NPCODE IF ( ABS(PC-PCODE(NP)) .LT. .001 ) GO TO 10 8 CONTINUE C IF THERE IS NO MATCH, IGNORE POINT GO TO 7 10 IPC(J) = NP C ESTIMATION OF NOZZLE INSTALLATION DRAG VIA TABLE LOOK-UP 12 IF ( INSDRG .GT. 0 ) 1 CALL ENGCD ( EM(J), AL(J), TH(J), A9, INSDRG ) J = J + 1 IF ( J .LE. NIP ) GO TO 7 WRITE(IU6,2000) J-1, I 20 IF ( IFILL .EQ. 0 ) NPCODE = 0 NIT = J - 2 IF ( NIT .LE. 0 ) GO TO 3120 C SORT ENGINE DECK BY MACH NUMBER DO 40 J = 1,NIT NIP = NIT - J + 1 DO 30 I = 1,NIP K1 = K(I) K2 = K(I+1) IF ( EM(K1) .GE. EM(K2) ) GO TO 30 K(I) = K2 K(I+1) = K1 30 CONTINUE 40 CONTINUE C BY ALTITUDE NAT = NA * NT - 1 DO 60 J = 1,NAT NIP = NIT - J + 1 DO 50 I = 1,NIP K1 = K(I) K2 = K(I+1) IF ( EM(K1) .NE. EM(K2) .OR. AL(K1) .GE. AL(K2) ) GO TO 50 K(I) = K2 K(I+1) = K1 50 CONTINUE 60 CONTINUE C BY THRUST LEVEL NT1 = NT - 1 DO 80 J = 1,NT1 NIP = NIT - J + 1 DO 70 I = 1,NIP K1 = K(I) K2 = K(I+1) IF ( EM(K1) .NE. EM(K2) .OR. AL(K1) .NE. AL(K2) .OR. 1 TH(K1) .GE. TH(K2) ) GO TO 70 K(I) = K2 K(I+1) = K1 70 CONTINUE 80 CONTINUE C PACK ENGINE DECK INTO ARRAYS NI = NIT + 1 M = 0 EMT = -3.0 NA = 1 NT = 1 IF ( NGPRT .GT. 1 ) WRITE(IU6, 360) DO 140 I = 1,NI J = K(I) IF ( ABS( EM(J) - EMT ) .LT. 0.009 ) GO TO 100 C NEW MACH NUMBER IF ( NGPRT .GT. 1 ) WRITE(IU6, 370) M = M + 1 EMT = EM(J) N = 1 ALS = AL(J) L = 0 IF ( M .LE. 20 ) THEN EMACH(M) = EMT ALT(N,M) = ALS ENDIF GO TO 120 100 IF ( ABS( AL(J) - ALS ) .LT. 100.0 ) GO TO 110 C NEW ALTITUDE FOR SAME MACH NUMBER IF ( NGPRT .GT. 1 ) WRITE(IU6, 370) N = N + 1 IF ( N .GT. NA ) NA = N L = 0 ALS = AL(J) IF ( N .LE. 15 .AND. M .LE. 20 ) ALT(N,M) = ALS GO TO 120 C CHECK FOR DUPLICATE THRUST 110 IF ( ABS( TH(J) - THR(L,N,M) ) .LT. 1.0 ) GO TO 130 120 L = L + 1 IF ( NPCODE .GT. 1 ) L = IPC(J) IF ( L .GT. NT ) NT = L IF ( L .GT. 16 .OR. N .GT. 15 .OR. M .GT. 20 ) GO TO 130 FF(L,N,M) = FL(J) IF ( NOX .EQ. 1 .OR. NOX .EQ. 3 ) OX(L,N,M) = XO(J) IF ( NOX .EQ. 2 .AND. FL(J) .GT. 0. ) 1 OX(L,N,M) = XO(J) * 1000. / FL(J) THR(L,N,M) = TH(J) 130 IF ( NGPRT .GT. 1 ) WRITE(IU6, 380) EM(J), AL(J), TH(J), FL(J), 1 XO(J) 140 CONTINUE IF ( NOX .EQ. 2 ) NOX = 1 NM = M C CHECK ARRAY DIMENSIONS 145 IF ( NM .LE. 1 ) THEN WRITE(IU6, 440) ISTOP = 1 ENDIF IF ( NM .GT. 20 ) THEN WRITE(IU6,2100) NM ISTOP = 1 NM = 20 ENDIF IF ( NA .GT. 15 ) THEN WRITE(IU6,2200) NA ISTOP = 1 NA = 15 ENDIF IF ( NT .GT. 16 ) THEN WRITE(IU6,2300) NT ISTOP = 1 NT = 16 ENDIF C FILL DATA TABLES IF POWER CODES ARE BEING USED IF ( NPCODE .LT. 2 ) GO TO 600 C FIND FIRST COMPLETE SET OF PART POWER DATA DO 470 MF = 1,NM DO 460 NF = 1,NA DO 450 I = 1,NT IF ( FF(I,NF,MF) .LE. 0. ) GO TO 460 450 CONTINUE GO TO 490 460 CONTINUE 470 CONTINUE C NO COMPLETE SETS - STOP WRITE(IU6, 480) 480 FORMAT(/' * * * NO COMPLETE SETS OF PART POWER DATA - STOP * * *') STOP C LOOP ON MACH NUMBER - CHECK FOR COMPLETE SETS AT THAT MACH NUMBER 490 DO 590 M = 1,NM IF ( MF .GE. M ) GO TO 520 DO 510 N = 1,NA DO 500 I = 1,NT IF ( FF(I,N,M) .LE. 0. ) GO TO 510 500 CONTINUE MF = M NF = N GO TO 520 510 CONTINUE C CHECK FOR UNDEFINED POINTS 520 DO 580 N = 1,NA IF ( N .EQ. NF .AND. M .EQ. MF ) GO TO 580 DO 570 L = 1,NT IF ( FF(L,N,M) .GT. 0. ) GO TO 570 I = L - 1 IF ( L .GT. 1 ) GO TO 540 C IF FIRST POINT IS UNDEFINED, FIND FIRST DEFINED POINT DO 530 I = 2,NT IF ( FF(I,N,M) .GT. 0. ) GO TO 540 530 CONTINUE GO TO 590 C INTERPOLATE/EXTRAPOLATE TO FILL DATA AT UNDEFINED POINT 540 THR(L,N,M) = THR(I,N,M) * THR(L,NF,MF) / THR(I,NF,MF) FF(L,N,M) = FF(I,N,M) * FF(L,NF,MF) / FF(I,NF,MF) IF ( NOX .GT. 0 .AND. OX(I,NF,MF) .GT. 0. ) 1 OX(L,N,M) = OX(I,N,M) * OX(L,NF,MF) / OX(I,NF,MF) 570 CONTINUE 580 CONTINUE 590 CONTINUE GO TO 250 C FILL IN PART POWER DATA IF POWER CODES ARE NOT USED 600 IF ( IFILL .EQ. 0 ) GO TO 250 DO 240 M = 1,NM DO 235 N = 1,NA C WHEN USING THE ENGINE CYCLE GENERATION CAPABILITY IN C OPTIMIZATION OR PARAMETRIC VARIATIONS, AN INCOMPLETE C ENGINE DECK MAY BE ACCEPTABLE. IF ( NICE .NE. 0 .AND. ALT(N,M) .GE. 0.0 .AND. * FF(1,N,M) .LE. 0.0 ) THEN C FIND THE NEXT ALTITUDE WITH GOOD DATA IF ( N .EQ. 1 ) THEN DO 560 KA = 2,NA 560 IF ( FF(1,KA,M) .GT. 0.0 ) GO TO 561 JECT = 1 RETURN 561 CONTINUE ELSE KA = N - 1 ENDIF C IF POSSIBLE, INTERPOLATE BETWEEN TWO ALTITUDES OR MACH NUMBERS IF (( N .NE. 1 .AND. N .NE. 15 ) .AND. * FF(1,N-1,M) .GT. 0.0 .AND. FF(1,N+1,M) .GT. 0.0 ) THEN THR(1,N,M) = ( THR(1,N-1,M) + THR(1,N+1,M) ) / 2.0 FF (1,N,M) = ( FF (1,N-1,M) + FF (1,N+1,M) ) / 2.0 IF ( NOX .GT. 0 ) * OX (1,N,M) = ( OX (1,N-1,M) + OX (1,N+1,M) ) / 2.0 ELSE IF (( M .NE. 1 .AND. M .NE. 20 ) .AND. * FF(1,N,M-1) .GT. 0.0 .AND. FF(1,N,M+1) .GT. 0.0 ) THEN THR(1,N,M) = ( THR(1,N,M-1) + THR(1,N,M+1) ) / 2.0 FF (1,N,M) = ( FF (1,N,M-1) + FF (1,N,M+1) ) / 2.0 IF ( NOX .GT. 0 ) * OX (1,N,M) = ( OX (1,N,M-1) + OX (1,N,M+1) ) / 2.0 ELSE c ohad 16/7/08 - 0.0d0 rather than 0. c CALL ATMO(ALT(KA,M),0.,DELT1,THET1,ASTAR,TM,RE,HFT) CALL ATMO(ALT(KA,M),0.0d0,DELT1,THET1,ASTAR,TM,RE, & HFT) SQTH1 = SQRT ( THET1 ) c ohad 16/7/08 - 0.0d0 rather than 0. c CALL ATMO(ALT(N,M),0.,DELT2,THET2,ASTAR,TM,RE,HFT) CALL ATMO(ALT(N,M),0.0d0,DELT2,THET2,ASTAR,TM,RE,HFT) SQTH2 = SQRT ( THET2 ) CTHR = THR(1,KA,M) / DELT1 CFF = FF (1,KA,M) / ( DELT1 * SQTH1 ) THR(1,N,M) = CTHR * DELT2 FF (1,N,M) = CFF * DELT2 * SQTH2 * FFPEN IF ( NOX .GT. 0 ) THEN COX = OX(1,KA,M) / ( DELT1 * SQTH1 ) OX (1,N,M) = COX * DELT2 * SQTH2 ENDIF ENDIF ISTOP = 0 IF ( IFLAG .EQ. 0 ) WRITE(IU6,2405) IFLAG = 1 ENDIF IF ( FF(1,N,M) .LE. 0.0 ) GO TO 235 DO 150 I = 1,IFILL IP1 = I + 1 IF ( FF(IP1,N,M) .LE. 0.0 ) GO TO 160 150 CONTINUE GO TO 235 160 M1 = M - 1 N1 = 0 IF ( N .GT. 1 ) N1 = N - 1 170 M1 = M1 + 1 IF ( M1 .GT. NM ) GO TO 250 DO 180 N2 = N,NA IF ( FF(IFILL + 1,N2,M1) .GT. 0.0 ) GO TO 190 180 CONTINUE IF ( N1 .GT. 0 ) GO TO 220 IF ( M .EQ. 1 ) GO TO 170 M1 = M - 1 N1 = 1 GO TO 220 190 IF ( N1 .GT. 0 ) GO TO 200 N1 = N2 GO TO 220 200 ALTF = (ALT(N1,M) - ALT(N,M)) / (ALT(N1,M) - ALT(N2,M)) DO 210 L = IP1,NT THR(L,N,M) = THR(I,N,M) * ( (1.0 - ALTF) * THR(L,N1,M) / 1 THR(I,N1,M) + ALTF * THR(L,N2,M) / 2 THR(I,N2,M) ) FF(L,N,M) = FF(I,N,M) * ( (1.-ALTF) * FF(L,N1,M) / 1 FF(I,N1,M) + ALTF * FF(L,N2,M) / FF(I,N2,M)) IF ( NOX .GT. 0 ) THEN IF ( OX(I,N1,M) .GT. 0. ) OX(L,N,M) = OX(I,N,M) * 1 (1.-ALTF) * OX(L,N1,M) / OX(I,N1,M) IF ( OX(I,N2,M) .GT. 0. ) OX(L,N,M) = OX(L,N,M) + 1 OX(I,N,M) * ALTF * OX(L,N2,M) / OX(I,N2,M) ENDIF 210 CONTINUE GO TO 235 220 DO 230 L = IP1,NT THR(L,N,M) = THR(I,N,M) * THR(L,N1,M1) / THR(I,N1,M1) FF(L,N,M) = FF(I,N,M) * FF(L,N1,M1) / FF(I,N1,M1) IF ( NOX .GT. 0 .AND. OX(I,N1,M1) .GT. 0. ) 1 OX(L,N,M) = OX(I,N,M) * OX(L,N1,M1) / OX(I,N1,M1) 230 CONTINUE 235 CONTINUE 240 CONTINUE C IDLE, BOOST AND ENGINE DATA OUTPUT 250 IF ( IGEO .GT. 0 ) WRITE(IU6, 430) IF ( NGPRT .GT. 0 ) WRITE(IU6, 385) IF ( IDLE .GT. 0 ) NT = NT + 1 IF ( NT .GT. 16 ) THEN WRITE(IU6,2300) NT ISTOP = 1 NT = 16 ENDIF DO 330 M = 1,NM DO 325 N = 1,NA IF ( FF(1,N,M) .LE. 0.0 ) GO TO 323 C CONVERT FROM GEOPOTENTIAL TO GEOMETRIC ALTITUDES IF ( IGEO .EQ. 0 ) GO TO 258 HFT = ALT(N,M) HO = HFT * .3048 Z = ( HFT + 4.37E-8 * HFT**2.0085 ) * .3048 257 R = REARTH + Z GN = GNS * (REARTH/R)**(CM1+1.) H = (R * GN * ((R/REARTH)**CM1 - 1.) / CM1 1 - Z * (R - Z/2.) * OC2) / GR DH = HO - H Z = Z + DH IF ( ABS(DH) .GT. 1. ) GO TO 257 ALT(N,M) = Z / .3048 258 IF ( NGPRT .GT. 0 ) WRITE(IU6, 390) EMACH(M),ALT(N,M) C FILL IN FLIGHT IDLE IF ( IDLE .EQ. 0 ) GO TO 290 IDROP = 0 IF ( FF(2,N,M) .GT. 0.0 ) GO TO 260 FF(2,N,M) = FIDMIN * FF(1,N,M) IF ( NOX .GT. 0 ) OX(2,N,M) = 0.50 * OX(1,N,M) GO TO 290 260 DO 270 L = 3,NT IF ( FF(L,N,M) .LE. 0.0 ) GO TO 280 270 CONTINUE GO TO 300 280 IF ( NONEG .EQ. 1 .AND. THR(L-1,N,M) .LT. 0. ) THEN L = L - 1 IDROP = 1 GO TO 280 ENDIF IF ( THR(L-1,N,M) .LE. 0. ) GO TO 300 LL = L + IDROP FF(L,N,M) = FF(LL-1,N,M) - (FF(LL-2,N,M) - FF(LL-1,N,M)) 1 * THR(LL-1,N,M) / (THR(LL-2,N,M) - THR(LL-1,N,M)) I = IDLE IF ( I .GE. L ) I = L - 1 IF ( FF(L,N,M) .LT. FIDMIN * FF(I,N,M) ) 1 FF(L,N,M) = FIDMIN * FF(I,N,M) IF ( FF(L,N,M) .GT. FIDMAX * FF(I,N,M) ) 1 FF(L,N,M) = FIDMAX * FF(I,N,M) IF ( FF(L,N,M) .GT. FF(L-1,N,M) ) FF(L,N,M) = FF(L-1,N,M) IF ( NOX .GT. 0 ) THEN OX(L,N,M) = OX(LL-1,N,M) - (OX(LL-2,N,M) - OX(LL-1,N,M)) 1 * THR(LL-1,N,M) / (THR(LL-2,N,M) - THR(LL-1,N,M)) IF ( OX(L,N,M) .LT. FIDMIN * OX(I,N,M) ) 1 OX(L,N,M) = FIDMIN * OX(I,N,M) ENDIF THR(L,N,M) = 0. C CHECK FOR ZERO FUEL FLOW AT SECOND POWER SETTING 290 IF ( FF(2,N,M) .GT. 0.0 ) GO TO 300 WRITE(IU6, 400) EMACH(M), ALT(N,M) ISTOP = 1 C FUDGE SUBSONIC AND SUPERSONIC FUEL FLOWS 300 DO 302 L = 1,NT IF ( EMACH(M) .LT. 1. ) THEN FF(L,N,M) = FFFSUB * FF(L,N,M) ELSE FF(L,N,M) = FFFSUP * FF(L,N,M) ENDIF 302 CONTINUE C REMOVE NEGATIVE THRUST POINTS IF ( NONEG .EQ. 1 ) THEN DO 305 L = 3,NT IF ( THR(L,N,M) .LT. 0. ) THEN THR(L,N,M) = 0. FF (L,N,M) = 0. OX (L,N,M) = 0. ENDIF 305 CONTINUE ENDIF C SPECIAL OPTION FOR BOOSTER ENGINES IF ( BOOST .LE. 0.0 .OR. 1 THR(1,N,M) .LE. 100000.0 ) GO TO 310 THR(1,N,M) = BOOST * ( THR(1,N,M) - 100000.0 ) + 1 THR(2,N,M) FF (1,N,M) = BOOST * FF(1,N,M) + FF(2,N,M) C PRINT ENGINE DECK 310 IF ( NGPRT .LE. 0 ) GO TO 325 WRITE(IU6, 410) (THR(L,N,M), L = 1,NT) WRITE(IU6, 410) (FF (L,N,M), L = 1,NT) DO 320 L = 1,NT IF ( THR(L,N,M) .EQ. 0.0 ) GO TO 315 SFC(L) = FF(L,N,M) / THR(L,N,M) GO TO 320 315 SFC(L) = 0.0 320 CONTINUE WRITE(IU6, 420) (SFC(L), L = 1,NT) IF ( NOX .EQ. 1 ) WRITE(IU6, 420) (OX(L,N,M), L = 1,NT) IF ( NOX .EQ. 3 ) WRITE(IU6, 410) (OX(L,N,M), L = 1,NT) GO TO 325 C CHECK FOR AT LEAST TWO ALTITUDES PER MACH NUMBER 323 IF ( N .GT. 2 ) GO TO 325 WRITE(IU6,2400) EMACH(M) ISTOP = 1 325 CONTINUE 330 CONTINUE C IF THERE IS AN ERROR IN THE ENGINE DECK, QUIT IF ( ISTOP .EQ. 1 ) STOP 335 CONTINUE CALL PLTTHR ( IST ) RETURN 340 FORMAT(///'# NAMELIST $ENGDIN',/ 1 ' ENGINE DECK CONTROL, SCALING AND USAGE DATA',// 2 14H DESCRIPTION,16X,4HNAME,9X,17HVALUE DIMENSIONS// 3 37H ENGINE DECK PRINT CONTROL NGPRT ,I11/ 4 37H ENGINE DECK SOURCE SWITCH IGENEN ,I11/ 5 37H SLOPE FACTOR FOR / 6 37H EXTRAPOLATING FUEL FLOWS EXTFAC ,F11.4/ 7 37H SUBSONIC FUEL FLOW FACTOR FFFSUB ,F11.4/ 8 37H SUPERSONIC FUEL FLOW FACT FFFSUP ,F11.4/ 9 37H FLIGHT IDLE SWITCH IDLE ,I11/ * 37H IGNORE NEGATIVE THRUSTS NONEG ,I11) 342 FORMAT( 37H MIN IDLE FUEL FLOW FRACT FIDMIN ,F11.4/ 1 37H MAX IDLE FUEL FLOW FRACT FIDMAX ,F11.4) 343 FORMAT( 37H SFC EXTRAPOLATION SWITCH IXTRAP ,I11/ 1 37H PART POWER DATA SWITCH IFILL ,I11/ 2 37H MAX. CRUISE POWER SETTING MAXCR ,I11/ 3 37H BOOST ENGINE SWITCH BOOST ,F11.4/ 4 37H FUEL FLOW SCALING / 5 37H CONSTANT TERM DFFAC ,F11.4/ 6 37H FUEL FLOW SCALING / 7 37H LINEAR TERM FFFAC ,F11.4/ 8 37H NITROGEN OXIDES SWITCH NOX ,I11/ 9 37H INSTALLATION DRAG SWITCH INSDRG ,I11) 344 FORMAT (37H DRAG COEF TABLE FILE NAME CDFILE ,5X,A / 1 37H AFTERBODY DRAG TABLE NO. NAB ,I11/ 2 37H REF AFTERBODY DRAG TABLE NABREF ,I11/ 3 37H MAXIMUM NOZZLE AREA A10 ,F11.2,6H SQ IN/ 4 37H REF NOZZLE EXIT AREA A9REF ,F11.2,6H SQ IN/ 5 37H REF MAXIMUM NOZZLE AREA A10REF ,F11.2,6H SQ IN/ 6 37H NOZZLE LENGTH XNOZ ,F11.2,3H IN/ 7 37H REFERENCE NOZZLE LENGTH XNREF ,F11.2,3H IN) 345 FORMAT (37H ENGINE DECK FILE NAME EIFILE ,5X,A ) 350 FORMAT ( F5.2,F10.0,F5.0,3F10.0,10X,2F10.0 ) 360 FORMAT (/'# SORTED ENGINE DECK',//' MACH NUMBER ALTITUDE', 1 ' NET THRUST FUEL FLOW NOX RATE') 370 FORMAT ( 1X ) 380 FORMAT ( F13.4,3F13.1,F13.3 ) 385 FORMAT ( /'# ALL POINT ENGINE DECK SUMMARY' ) 390 FORMAT ( /8H MACH =,F7.3,13H, ALTITUDE =,F8.0, 1 44H, THRUSTS/FUEL FLOWS/SFCS/NOX RATIOS FOLLOW) 400 FORMAT ( /44H * * * ONLY ONE THRUST LEVEL FOR MACH NUMBER,F6.3, 1 13H FOR ALTITUDE,F8.0,6H * * */ ) 410 FORMAT ( 1X,16F8.1 ) 420 FORMAT ( 1X,16F8.3 ) 430 FORMAT ( /42H * ALTITUDES CHANGED FROM GEOPOTENTIAL TO , 1 11HGEOMETRIC *) 440 FORMAT ( /40H * * * ONLY ONE MACH NUMBER - STOP * * * ) 2000 FORMAT ( /11H * * * ONLY,I5,13H OF THE FIRST,I6, 1 30H ENGINE POINTS WERE USED * * */ ) 2100 FORMAT ( /6H * * *,I3,31H IS TOO MANY MACH NUMBERS * * */ ) 2200 FORMAT ( /6H * * *,I3,28H IS TOO MANY ALTITUDES * * */ ) 2300 FORMAT ( /6H * * *,I3,32H IS TOO MANY THRUST LEVELS * * */ ) 2400 FORMAT ( /40H * * * ONLY ONE ALTITUDE FOR MACH NUMBER, 1 F6.3,6H * * */ ) 2405 FORMAT ( /, * ' * * * DATA FOR MISSING ALTITUDES ARE APPROXIMATED. * * *') 2500 FORMAT ( /' MACH NUMBERS AND ALTITUDES FOR ENGINE DECK GENERATION' 1 ,//' MACH ALTITUDES',/) 2510 FORMAT ( F7.2,3X,15F8.0) 2520 FORMAT ( /' * * * INPUT MACH NUMBERS ARE NOT IN DESCENDING ORDER', 1 ' * * *') 2530 FORMAT ( /' * * * INPUT ALTITUDES ARE NOT IN DESCENDING ORDER', 1 ' * * *') 3000 IF ( I .GT. 1 ) GO TO 3100 IF ( IU11 .EQ. IU5 ) GO TO 7 REWIND ( IU11 ) READ (IU11,3010) TITLE 3010 FORMAT ( A80 ) WRITE(IU6,3020) TITLE 3020 FORMAT (/' ENGINE DECK TITLE',//,' ',A80) GO TO 7 3100 WRITE(IU6,3110) 3110 FORMAT ( /' * * * ILLEGAL DATA IN ENGINE DECK * * *',/ ) STOP 3120 WRITE(IU6,3130) 3130 FORMAT ( /' * * * ENGINE DECK MISSING * * *',/ ) STOP 3140 WRITE(IU6,3150) EIFILE, IU11 3150 FORMAT(/' * * * ERROR OPENING FILE ',A,/ 1 ' AS UNIT',I3,' * * *') STOP 9000 WRITE(IU6,9001) IU5 STOP 9001 FORMAT (//' ERROR READING NAMELIST $ENGDIN FROM UNIT',I3,/, 2 ' PROGRAM ABORTED IN SUBROUTINE DEFENG.',/) END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEPARATOR CCCCCCCCCCCCCCCCCCCCCCCCC BLOCK DATA EDATIN IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) C DEFAULT VALUES FOR ENGINE RELATED INPUT DATA NAMELIST $ENGDIN COMMON /ENDTA / EXTFAC,FFFAC ,DFFAC ,TAKOFF,EMACH(20) , 1 ALT(15,20) ,FF(16,15,20) ,THR(16,15,20), 2 TXFUFL,NM ,NA ,NT ,IXTRAP,MAXCR COMMON /NOXDAT/ OX(16,15,20), FFUEL(6), FNOX(6), RTNOX, TNOX, 1 GNOX, OFNOX, NOX COMMON /CDFIL / CDFILE, A10, A9REF, A10REF, XNOZ, XNREF, RCRV, 1 NAB, NABREF COMMON /THRPLT/ PCODE(16), NPCODE, IPLTTH COMMON / FFFF / FFFSUB, FFFSUP CHARACTER*80 CDFILE DATA DFFAC, FFFAC, TAKOFF, TXFUFL /4*0./, EXTFAC/1./, MAXCR /2/, 1 IXTRAP /1/, NM/20/, NA/15/, NT/16/, EMACH, ALT /320*-1./ DATA FFUEL/6*1./, FNOX, GNOX, OFNOX /8*0./, NOX/0/ DATA FFFSUB, FFFSUP /2*1./ DATA CDFILE/'ENDRAG'/, NAB/6969/, NABREF/6969/ DATA A10, A9REF, A10REF, XNOZ, XNREF/5*0./, RCRV/-1./ DATA IPLTTH, NPCODE /2*0/, PCODE /1., 2., 3., 4., 5., 6., 7., 8., 1 9., 10., 11., 12., 13., 14., 15., 16./ END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEPARATOR CCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE THTOL ( THRUST ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) C PREPARE THRUST TABLES FOR TAKEOFF AND LANDING MODULE C THRUST CURRENT ENGINE THRUST REFERENCE VALUE COMMON /TOLTH / APA, ANS1, DTCT, SPDSND, THFACT, THREF, THRTO(10), 1 THRLD(10), VELTO(10), VELLD(10), INTHTO, INTHLD COMMON /TREVRS/ RVFACT, TIRVRS, TIRVA, REVCUT, CLREV, CDREV, 1 VELRV(10), THRRV(10), INTHRV CALL ATMO ( APA,DTCT,DELTA,THETA,ASTAR,TM,RE,HFT ) SPDSND = ASTAR * 6076.1155 / 3600.0 ANS1 = 0.0023769 * DELTA / THETA ESCALE = THRUST / THREF C TAKEOFF THRUST DATA IF ( INTHTO- 1 ) 50,10,30 C SCALE INPUT VALUES 10 DO 20 I = 1,10 THRTO(I) = THRTO(I) * ESCALE 20 CONTINUE GO TO 50 C USE MAIN ENGINE DECK SPECIFIED SETTING 30 DO 40 I = 1,10 XMACH = VELTO(I) / SPDSND CALL ENINT ( XMACH,APA,THU,FDUM,1,TMX,ITD,INTHTO- 1 ) THRTO(I) = THU * THFACT 40 CONTINUE C LANDING THRUST DATA 50 IF ( INTHLD- 1 ) 100,60,80 C SCALE INPUT VALUES 60 DO 70 I = 1,10 THRLD(I) = THRLD(I) * ESCALE 70 CONTINUE GO TO 100 C USE MAIN ENGINE DECK IDLE SETTINGS 80 DO 90 I = 1,10 XMACH = VELLD(I) / SPDSND CALL ENINT ( XMACH,APA,THU,FDUM,3,TMX,ITD,1 ) THRLD(I) = THU 90 CONTINUE C REVERSE THRUST DATA 100 IF ( INTHRV .GT. 1 ) GO TO 130 IF ( INTHRV ) 160,180,110 C SCALE INPUT VALUES 110 DO 120 I = 1,10 THRRV(I) = THRRV(I) * ESCALE 120 CONTINUE GO TO 180 C USE MAIN ENGINE DECK SPECIFIED POWER SETTING 130 DO 140 I = 1,10 XMACH = VELRV(I) / SPDSND CALL ENINT ( XMACH,APA,THU,FDUM,1,TMX,ITD,INTHRV - 1 ) THRRV(I) = THU 140 CONTINUE GO TO 180 C USE TAKEOFF THRUST DATA 160 DO 170 I = 1,10 THRRV(I) = THRTO(I) VELRV(I) = VELTO(I) 170 CONTINUE 180 THREF = THRUST RETURN END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEPARATOR CCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SKLENG ( ENGSKL,IFLAG ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) C SCALES THRUST TABLES USING ENGSKL AND FUEL FLOW TABLES USING ENGSKL C AND NONLINEAR FACTORS FFFAC AND DFFAC C ENGSKL = ENGINE SCALE FACTOR C IFLAG = PRINT FLAG COMMON /ENDTA / EXTFAC,FFFAC ,DFFAC ,TAKOFF,EMACH(20) , 1 ALT(15,20) ,FF(16,15,20) ,THR(16,15,20), 2 TXFUFL,NM ,NA ,NT ,IXTRAP,MAXCR COMMON /UNITS / IU5, IU6, IU7, IU8, IU9, IU16, IU17, IU18 C----------------------------------------------------------------------C C MODIFICATIONS FOR POINT MODE OPERATION C LOGICAL ENGOPT C COMMON /ENGPRF/ CDT, CDP, VJET, CDTMAX, CDPMAX, VJMAX, ENGFN, C 1 ENGSFC, OVEFF, SPECT, STMIN, C 2 ARATIO, ARMAX, XACRS, XMCRS, IGENEN C COMMON /RCOMENG/ TABMAX, C * EFFAB, XMDES, XADES, DESFN, C * COSTBL, HPEXT, FHV, XIDLE, YEAR C----------------------------------------------------------------------C DIMENSION SFC(16) SAVE TOTSKL, RSCAL DATA TOTSKL, RSCAL / 2*1. / C IF THE INPUT SCALING FACTOR IS 1.0, NO SCALING IF PERFORMED IF ( ENGSKL .EQ. 1.0 ) GO TO 40 C CALCULATE NEW TOTAL SCALE FACTOR AND NONLINEAR FUEL FLOW SCALE C FACTOR AND SCALE TAXI AND TAKEOFF FUEL FLOWS TOTSKL = TOTSKL * ENGSKL SCALE = TOTSKL * ( FFFAC * ( 1.0 - TOTSKL ) + DFFAC + 1.0 ) C----------------------------------------------------------------------C C MODIFICATIONS FOR POINT MODE OPERATION C IF ( IGENEN .EQ. 2 ) THEN C DESFN = DESFN * ENGSKL C IGENEN = 12 C ENGOPT = .FALSE. C CALL GENENG ( DESSFC, ENGOPT ) C RETURN C ENDIF C----------------------------------------------------------------------C C SCALE MAIN ENGINE DECK DO 10 L = 1,NM DO 5 N = 1,NA DO 3 M = 1,NT FF (M,N,L) = FF(M,N,L) * SCALE / RSCAL THR(M,N,L) = THR(M,N,L) * ENGSKL 3 CONTINUE 5 CONTINUE 10 CONTINUE RSCAL = SCALE C PRINT SCALED ENGINE TABLES IF ( IFLAG .LT. 3 ) GO TO 40 WRITE(IU6, 50) TOTSKL,FFFAC ,DFFAC DO 30 M = 1,NM DO 25 N = 1,NA IF ( FF(1,N,M) .LE. 0.0 ) GO TO 30 WRITE(IU6, 70) EMACH(M),ALT(N,M) WRITE(IU6, 80) ( THR(L,N,M),L = 1,NT ) WRITE(IU6, 80) ( FF (L,N,M),L = 1,NT ) DO 20 L = 1,NT SFC(L) = 0.0 IF (THR(L,N,M) .GT. 0.) SFC(L) = FF(L,N,M) / THR(L,N,M) 20 CONTINUE WRITE(IU6, 90) ( SFC(L),L = 1,NT ) 25 CONTINUE 30 CONTINUE 40 ENGSKL = TOTSKL RETURN 50 FORMAT (/' SCALED ENGINE DATA - ENGSKAL = ',F7.4,5X,'FFFAC=', 1 F7.4,5X,'DFFAC=',F7.4 ) 70 FORMAT ( /8H MACH =,F7.3,13H, ALTITUDE =,F8.0, 1 33H, THRUSTS/FUEL FLOWS/SFCS FOLLOW) 80 FORMAT ( 1X,16F8.1 ) 90 FORMAT ( 1X,16F8.5 ) END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEPARATOR CCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ENINT ( E,A,T,F,K,TMAX,IT,IPP ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) C ENGINE DECK INTERPOLATION - FOR A GIVEN MACH NUMBER AND ALTITUDE, C FIND THE MAX THRUST OR FLIGHT IDLE THRUST AND CORRESPONDING C FUEL FLOW OR FOR A GIVEN THRUST FIND THE FUEL FLOW C E = INPUT MACH NUMBER C A = INPUT ALTITUDE C T = THRUST - EITHER INPUT OR OUTPUT C F = FUEL FLOW - OUTPUT C K = 1, FIND THRUST AND FUEL FLOW FOR POWER SETTING NUMBER IPP C K = 2, FIND FUEL FLOW FOR THRUST T C K = 3, FIND FLIGHT IDLE THRUST AND FUEL FLOW C TMAX = THRUST FOR POWER SETTING IPP (MAX THRUST IF IPP = 1) C IT = 1 INDICATES THAT REQUESTED THRUST EXCEEDS TMAX C IPP = DESIRED POWER SETTING FOR K = 1 COMMON /ENDTA / EXTFAC,FFFAC ,DFFAC ,TAKOFF,EMACH(20) , 1 ALT(15,20) ,FF(16,15,20) ,THR(16,15,20), 2 TXFUFL,NM ,NA ,NT ,IXTRAP,MAXCR COMMON /NOXDAT/ OX(16,15,20), FFUEL(6), FNOX(6), RTNOX, TNOX, 1 GNOX, OFNOX, NOX DIMENSION NB(2), FA(2), FT(4), TT(4), OT(4) CP COMMON /ENGPRF/ CDT, CDP, VJET, CDTMAX, CDPMAX, VJMAX, ENGFN, CP 1 ENGSFC, OVEFF, SPECT, STMIN, CP 2 ARATIO, ARMAX, XACRS, XMCRS, IGENEN CPC SAVE THE MACH AND ALTITUDE FROM THE LAST CALL CP SAVE EL, AL CPC SAVE THE CORRECTED (THRUST AND FUEL FLOW) FROM THE LAST CALL CP SAVE TCORL, FCORL CPC SAVE THE CORRECTED MAXIMUM THRUST FROM THE LAST CALL CP SAVE TMAXCL CP SAVE KLAST, IPLAST IA = 0 IT = 0 CPC--------------------------------------------------------------------C CPC MODIFICATIONS FOR POINT MODE OPERATION ARE INDICATED WITH A 'CP' CPC INSERTED AT THE START OF EACH LINE TO DISABLE IT. POINT MODE SEEMS CPC TO BE VERY EXPENSIVE WITH LITTLE OR NO ADVANTAGE. THERE IS ALSO CPC POINT MODE CODE IN SUBROUTINE SKLENG. CPC USE IPP=-2 FOR MAX A/B, IPP=-3 FOR MAX DRY, AND IPP=-4 FOR IDLE CP IF ( IGENEN .EQ. 2 ) THEN CP TMAX = 0.0 CP CPC--------------------------------------------------------------------C CPC USE APPROXIMATIONS FOR INCREASED SPEED CPC IF THE MACH AND CORRECTED THRUST ARE ESSENTIALLY UNCHANGED, CPC USE NON-DIMENSIONAL ANALYSIS TO ESTIMATE THE FUEL FLOW FOR CPC SMALL CHANGES IN ALTITUDE. CP DELT = DELTA(A) CP TCOR = T / DELT CP TDIFF = ABS((TCOR-TCORL)/MAX(1.0,TCORL)) CP IF ( K.EQ.1 .AND. KLAST.EQ.1 .AND. IPP.EQ.IPLAST ) TDIFF = 0.0 CPC CALL POINT WITH TMAX < -100 TO USE LAST TEMP AS INITIAL GUESS. CP IF ( TDIFF .LE. 0.05000 .AND. TCORL .GT. 1.0 ) TMAX = -1000.0 CP IF ( ABS((E-EL)/MAX(1.0,EL)) .LE. 0.00010 .AND. CP * ABS((A-AL)/MAX(1.0,AL)) .LE. 0.05000 .AND. CP * TDIFF .LE. 0.00500 ) THEN CP STHT = SQRT(THETA(A)) CP IF ( K .EQ. 3 ) T = TCORL * DELT CP F = FCORL * DELT * STHT CP TMAX = TMAXCL * DELT CP RETURN CP ENDIF CPC--------------------------------------------------------------------C CP CP IF ( K .EQ. 1 .OR. K .EQ. 3 ) THEN CP F = -2 CP TMAX = 0.0 CP IF ( IPP .LE. -3 .OR. IPP .EQ. NT .OR. K .EQ. 3 ) F = -1.0 CP CALL POINT ( E, A, T, F, TMAX , 0 ) CP IF ( IPP .EQ. -4 .OR. IPP .EQ. NT .OR. K .EQ. 3 ) THEN CP T = T * 0.050 CP F = 0.00 CP CALL POINT ( E, A, T, F, TMAX , 0 ) CP ENDIF CP ELSE IF ( K .EQ. 2 ) THEN CP FN = T CP CALL POINT ( E, A, FN, F, TMAX , 0 ) CP IF ( FN .LT. T*1.000001) THEN CP IT = 1 CP T = FN CP ENDIF CP ENDIF CP EL = E CP AL = A CP DELT = DELTA(A) CP TCORL = T / DELT CP TMAXCL = TMAX / DELT CP FCORL = F / (DELT * SQRT(THETA(A))) CP KLAST = K CP IPLAST = IPP CP RETURN CP ENDIF CPC--------------------------------------------------------------------C C FIND MACH NUMBER BRACKET DO 10 M = 2,NM IF ( EMACH(M) .LE. E ) GO TO 20 10 CONTINUE M = NM 20 MA = M - 1 FM = 1.0 - ( E - EMACH(MA) ) / ( EMACH(MA + 1) - EMACH(MA) ) C FIND ALTITUDE BRACKET FOR EACH MACH NUMBER DO 50 M = 1,2 I = MA - 1 + M DO 30 N = 2,NA IF ( ALT(N,I) .LE. A .OR. 1 FF(1,N + 1,I) .LE. 0.0 ) GO TO 40 30 CONTINUE N = NA 40 NB(M) = N - 1 FA(M) = 1.0 - ( A - ALT(NB(M),I) ) / ( ALT(NB(M) + 1,I) - 1 ALT(NB(M),I) ) 50 CONTINUE IF ( A .GT. ( FM * ALT(1,MA) + (1.0 - FM) * ALT(1,MA + 1) ) ) 1 IA=1 IF ( A .LT. ( FM * ALT(NB(1) + 1,MA) + (1.0 - FM) * 1 ALT(NB(2) + 1,MA + 1) ) ) IA = -1 J = 0 IF ( K .EQ. 3 ) GO TO 130 C FIND MAXIMUM THRUST (K = 1 FOR CLIMB OR CRUISE) C ALSO REQUIRED FOR GIVEN THRUST (K = 2) J1 = IPP J2 = IPP J3 = IPP J4 = IPP IF ( FF(3,NB(1),MA) .EQ. 0.0 ) J1 = 1 IF ( FF(3,NB(1)+1,MA) .EQ. 0.0 ) J2 = 1 IF ( FF(3,NB(2),MA+1) .EQ. 0.0 ) J3 = 1 IF ( FF(3,NB(2)+1,MA+1) .EQ. 0.0 ) J4 = 1 TMAX = FM * ( FA(1) * THR(J1,NB(1),MA) + (1.0 - FA(1)) * 1 THR(J2,NB(1)+1,MA) ) + (1.0 - FM) * ( FA(2) * 2 THR(J3,NB(2),MA+1) + (1.0 - FA(2)) * 3 THR(J4,NB(2)+1,MA+1) ) IF ( K .EQ. 2 ) GO TO 70 C FIND FUEL FLOW FOR MAXIMUM THRUST T = TMAX F = FM * ( FA(1) * FF(J1,NB(1),MA) + (1.0 - FA(1)) * 1 FF(J2,NB(1)+1,MA) ) + (1.0 - FM) * ( FA(2) * 2 FF(J3,NB(2),MA+1) + (1.0 - FA(2)) * FF(J4,NB(2)+1,MA+1)) IF ( NOX .GT. 0 ) RTNOX = FM * ( FA(1) * OX(J1,NB(1),MA) + (1.0 - 1 FA(1)) * OX(J2,NB(1)+1,MA) ) + (1.0 - FM) * ( FA(2) * 2 OX(J3,NB(2),MA+1) + (1.0 - FA(2)) * OX(J4,NB(2)+1,MA+1)) IF ( IA * IXTRAP .EQ. 0 ) RETURN C LIMIT SFC FOR ALTITUDE EXTRAPOLATION IF(IA.LT.0) GO TO 60 C INPUT ALTITUDE HIGHER THAN TABLES SD = ( FM * THR(J1,1,MA) + (1.0 - FM) * THR(J3,1,MA+1) ) IF ( SD .LE. 0. .OR. T .LE. 0. ) RETURN SFCM = ( FM * FF(J1,1,MA) + (1.0 - FM) * FF(J3,1,MA+1) ) / SD IF ( F / T .LT. SFCM ) F = T * SFCM RETURN C INPUT ALTITUDE LOWER THAN TABLES 60 SD = FM * THR(J2,NB(1)+1,MA) + (1.0 - FM) * THR(J4,NB(2)+1,MA+1) IF ( SD .LE. 0. .OR. T .LE. 0. ) RETURN SFCM = ( FM * FF(J2,NB(1)+1,MA) + (1.0 - FM) * 1 FF(J4,NB(2)+1,MA+1) ) / SD IF ( F / T .LT. SFCM ) F = T * SFCM RETURN C FIND FUEL FLOW FOR GIVEN THRUST (K = 2) 70 IF ( T .LE. 2. ) THEN P = T T = P * TMAX ELSE P = T / TMAX ENDIF IF ( P .GT. 1.0 ) IT = 1 EF = 1.0 IF ( P .GT. 1.0 ) EF = EXTFAC DO 120 M = 1,2 I = MA - 1 + M DO 110 N = 1,2 II = NB(M) - 1 + N TM = THR(IPP,II,I) J = J + 1 TT(J) = TM DO 80 L = 2,NT IF ( THR(L,II,I) / TM .LE. P ) GO TO 90 80 CONTINUE L = NT 90 LA = L - 1 IF ( FF(L,II,I) .GT. 0.0 ) GO TO 100 IF ( LA .EQ. 1 ) GO TO 100 L = L - 1 GO TO 90 100 FFF = ( P * TM - THR(LA,II,I)) / 1 ( THR(LA+1,II,I) - THR(LA,II,I) ) FT(J) = FF(LA,II,I) + FFF*EF*(FF(LA+1,II,I) - FF(LA,II,I)) IF ( NOX .GT. 0 ) OT(J) = OX(LA,II,I) 1 + FFF*EF*(OX(LA+1,II,I) - OX(LA,II,I)) 110 CONTINUE 120 CONTINUE GO TO 180 C FIND FLIGHT IDLE THRUST AND FUEL FLOW (K = 3) 130 DO 170 M = 1,2 I = MA - 1 + M DO 160 N = 1,2 II = NB(M) - 1 + N J = J + 1 DO 140 L = 2,NT LP = NT - L + 2 IF ( FF(LP,II,I) .GT. 0.0 ) GO TO 150 140 CONTINUE 150 TT(J) = THR(LP,II,I) FT(J) = FF(LP,II,I) OT(J) = OX(LP,II,I) 160 CONTINUE 170 CONTINUE T = ( FA(1) * TT(1) + (1.0 - FA(1)) * TT(2) ) * FM 1 + ( FA(2) * TT(3) + (1.0 - FA(2)) * TT(4) ) * (1.0 - FM) 180 F = ( FA(1)*FT(1) + (1.0 - FA(1))*FT(2) ) * FM 1 +( FA(2)*FT(3) + (1.0 - FA(2))*FT(4) ) * (1.0 - FM) IF ( NOX .GT. 0 ) RTNOX = ( FA(1)*OT(1) + (1.0 - FA(1))*OT(2) ) 1 * FM + ( FA(2)*OT(3) + (1.0 - FA(2))*OT(4) ) * (1.0 - FM) C LIMIT SFC FOR ALTITUDE EXTRAPOLATION IF ( IA * IXTRAP .EQ. 0 .AND. F .GE. 0.0 ) RETURN IF ( K .EQ. 3 ) GO TO 200 IF ( IA .LT. 0 ) GO TO 190 SFCM = ( FM * FT(1) + (1.0 - FM) * FT(3) ) / 1 ( FM * TT(1) + (1.0 - FM) * TT(3) ) / P IF ( F / T .LT. SFCM ) F = SFCM * T RETURN 190 SFCM = ( FM * FT(2) + (1.0 - FM) * FT(4) ) / 1 ( FM * TT(2) + (1.0 - FM) * TT(4) ) / P IF ( F / T .LT. SFCM ) F = SFCM * T RETURN C LIMIT FUEL FLOW FOR FLIGHT IDLE 200 IF ( IA .LT. 0 ) GO TO 210 SFCM = FM * FT(1) + (1.0 - FM) * FT(3) IF ( F .LT. SFCM ) F = SFCM RETURN 210 SFCM = FM * FT(2) + (1.0 - FM) * FT(4) IF ( F .LT. SFCM ) F = SFCM RETURN END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEPARATOR CCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ATMO ( ZFT, DTC, DELTA, THETA, ASTAR, TM, RE, HFT ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) C 1962 STANDARD ATMOSPHERIC PROPERTIES GOOD UP TO 88743 METERS C GEOPOTENTIAL ALTITUDE (90 KM GEOMETRIC ALTITUDE OR 291152 FEET) C ALSO SAME AS 1976 STD ATMOSPHERE TO 51 KM (167323 FEET) C INPUT/OUTPUT IN ENGLISH UNITS, CALCULATIONS IN SI UNITS C BASE PRESSURES AND EXPONENTS FOR EACH LAYER ARE RECOMPUTED TO ASSURE C CONTINUITY AT THE CORNERS REGARDLESS OF THE COMPUTER USED C ZFT INPUT ALTITUDE - FEET C DTC DELTA TEMPERATURE FROM STD - DEG C C DELTA PRESSURE RATIO C THETA TEMPERATURE RATIO C ASTAR SPEED OF SOUND - KNOTS C TM MOLECULAR-SCALE TEMPERATURE - DEG KELVIN C RE REYNOLDS NUMBER PER FOOT AT MACH 1. C HFT GEOPOTENTIAL ALTITUDE - FEET DOUBLE PRECISION ZFT, DTC, DELTA, THETA, ASTAR, TM, RE, HFT, 1 SFT, STC DIMENSION P(9), E(9) SAVE P, E, DDLTA, DHETA, DSTAR, DTM, DRE, DFT, SFT, STC, IFIR DATA REARTH/6367533./,GR/9.80665/,GNS/9.823693/,CM1/.9985/, 1 OC2/26.76566D-10/,IFIR/1/,SFT,STC/2*-1000./ C PRECALCULATE EXPONENTS AND BASE PRESSURE RATIOS IF ( IFIR .NE. 1 ) GO TO 5 P(1) = 1. GMOR = 9.80665 * 1.225 * 288.15 / 101.325 E(1) = GMOR / 6.5 P(2) = (216.65/288.15)**E(1) E(2) = -GMOR / 216.65 P(3) = P(2) * EXP(E(2) * 9.) E(3) = GMOR P(4) = P(3) * (216.65/228.65)**GMOR E(4) = GMOR / 2.8 P(5) = P(4) * (228.65/270.65)**E(4) E(5) = -GMOR / 270.65 P(6) = P(5) * EXP(E(5) * 5.) E(6) = GMOR / 2. P(7) = P(6) * (252.65/270.65)**E(6) E(7) = GMOR / 4. P(8) = P(7) * (180.65/252.65)**E(7) E(8) = -GMOR / 180.65 Z90 = 90000. R = REARTH + Z90 GN = GNS * (REARTH / R)**(CM1 + 1.) H90 = (R * GN * ( (R/REARTH)**CM1 - 1.) / CM1 1 - Z90 * (R - Z90/2.) * OC2) / GR DH = H90/1000. - 79. P(9) = P(8) * EXP(E(8) * DH) E(9) = 11.056 IFIR = 0 C CONVERT INPUT FEET TO METERS 5 IF ( ZFT .EQ. SFT .AND. DTC .EQ. STC ) GO TO 110 SFT = ZFT STC = DTC Z = ZFT * .3048 C CALCULATE GEOPOTENTIAL ALTITUDE R = REARTH + Z RPOW = (REARTH / R)**CM1 GN = GNS * (REARTH / R) * RPOW H = (R * GN * ( 1./RPOW - 1.) / CM1 1 - Z * (R - Z/2.) * OC2) / GR DFT = H / .3048 C CONVERT H TO KILOMETERS H = H/1000. C SEA LEVEL TO 11 KM IF ( H .GT. 11. ) GO TO 11 DTM = 288.15 - 6.5 * H DDLTA = ((DTM)/288.15)**E(1) GO TO 100 C 11 KM TO 20 KM 11 IF ( H .GT. 20. ) GO TO 20 DH = H - 11. DTM = 216.65 DDLTA = P(2) * EXP(E(2) * DH) GO TO 100 C 20 KM TO 32 KM 20 IF ( H .GT. 32. ) GO TO 32 DH = H - 20. DTM = 216.65 + DH DDLTA = P(3) * (216.65/DTM)**E(3) GO TO 100 C 32 KM TO 47 KM 32 IF ( H .GT. 47. ) GO TO 47 DH = H - 32. DTM = 228.65 + 2.8 * DH DDLTA = P(4) * (228.65/DTM)**E(4) GO TO 100 C 47 KM TO 52 KM 47 IF ( H .GT. 52. ) GO TO 52 DH = H - 47. DTM = 270.65 DDLTA = P(5) * EXP(E(5) * DH) GO TO 100 C 52 KM TO 61 KM 52 IF ( H .GT. 61. ) GO TO 61 DH = H - 52. DTM = 270.65 - 2.0 * DH DDLTA = P(6) * (DTM/270.65)**E(6) GO TO 100 C 61 KM TO 79 KM 61 IF ( H .GT. 79. ) GO TO 79 DH = H - 61. DTM = 252.65 - 4.0 * DH DDLTA = P(7) * (DTM/252.65)**E(7) GO TO 100 C 79 KM TO 88743 METERS 79 IF ( Z .GT. 90000. ) GO TO 90 DH = H - 79. DTM = 180.65 DDLTA = P(8) * EXP(E(8) * DH) GO TO 100 C ABOVE 88743 M, 1962 STD ATMOSPHERE SWITCHES TO GEOMETRIC ALTITUDE C THE EQUATIONS BELOW ARE CLOSE UP TO 100 KM AND DIVERGE AFTER THAT 90 DH = Z/1000. - 90. DTM = 180.65 + 3.0 * DH DDLTA = P(9) * (180.65/DTM)**E(9) C CALCULATE TEMPERATURE RATIO, SPEED OF SOUND, AND REYNOLDS NUMBER 100 DHETA = (DTM + DTC) / 288.15 DSTAR = 661.479 * SQRT(DHETA) DRE = 1.479301E+9 * DDLTA * (DTM + 110.4) / DTM**2 110 DELTA = DDLTA THETA = DHETA ASTAR = DSTAR TM = DTM RE = DRE HFT = DFT RETURN END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEPARATOR CCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ORIDE ( VAL, OVAL ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) C PROVIDES A CAPABILITY TO OVERRIDE CERTAIN CALCULATED VALUES C VAL INPUT CALCULATED VALUE, OUTPUT OVERRIDDEN VALUE C OVAL OVERRIDE PARAMETER C > 5. REPLACEMENT VALUE C < 0. NEGATIVE OF REPLACEMENT VALUE THAT IS SCALED IN C SUBSEQUENT ITERATIONS C 0. < OVAL < 5. MULTIPLIER IF ( OVAL .GT. 5.0 ) GO TO 10 IF ( VAL .EQ. 0.0 ) RETURN IF ( OVAL .LT. 0.0 ) OVAL = -OVAL / VAL VAL = VAL * OVAL RETURN 10 VAL = OVAL RETURN END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEPARATOR CCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION FUNCTION XINT1 ( X, XST, NX, FST ) C LINEAR INTERPOLATION ROUTINE C X INDEPENDENT VARIABLE C XST ARRAY OF VALUES OF X C NX NUMBER OF ARRAY ELEMENTS C FST DEPENDENT VARIABLE ARRAY C XINT1 RESULT OF INTERPOLATION IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) DIMENSION XST(NX), FST(NX) DELX = SEARCH ( X, XST, IX, NX ) XINT1 = DELX * FST(IX+1) + (1.-DELX) * FST(IX) RETURN END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEPARATOR CCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION FUNCTION SEARCH ( X, XST, IX, NX ) C FIND LOCATION OF INDEPENDENT VARIABLE X IN XST ARRAY C IX ELEMENT NUMBER PRECEDING X POSITION C NX NUMBER OF ELEMENTS IN ARRAY C SEARCH FRACTION OF DISTANCE FROM X(IX) TO X(IX+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) DIMENSION XST(NX) NXM1 = NX - 1 IF ( XST(1) .GT. XST(2) ) GO TO 12 C ASCENDING ARRAY DO 13 I = 1,NXM1 IF ( X .LE. XST(I+1) ) GO TO 14 13 CONTINUE 16 I = NXM1 14 SEARCH = (X-XST(I)) / (XST(I+1)-XST(I)) IX = I RETURN C DESCENDING ARRAY 12 DO 15 I = 1,NXM1 IF ( X .GE. XST(I+1) ) GO TO 14 15 CONTINUE GO TO 16 END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEPARATOR CCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ENGCD ( XM, ALT, THRUST, A9, INSDRG ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) C----------------------------------------------------------------------C C THIS SUBROUTINE ESTIMATES INSTALLATION DRAG VIA TABLE LOOK-UP. C C C C CURRENTLY USED TO ACCOUNT FOR CHANGES IN DRAG DUE TO A CHANGE C C IN NOZZLE DESIGN (CHANGES IN A10 RELATIVE TO THE BASELINE). C C C C INPUT (NAMELIST $ENGDIN): C C CDFILE NAME OF THE FILE CONTAINING THE TABLE OF DRAG C C COEFFICIENTS. C C NAB TABLE NUMBER TO BE USED FOR AFTERBODY DRAG. C C NABREF TABLE NUMBER TO BE USED FOR REFERENCE AFTERBODY DRAG. C C A10 MAXIMUM NOZZLE AREA, IN*IN C C A10REF REFERENCE MAXIMUM NOZZLE AREA, IN*IN C C XNOZ NOZZLE LENGTH, IN C C XNREF REFERENCE NOZZLE LENGTH, IN C C A9REF REFERENCE NOZZLE EXIT AREA, IN*IN C C C C INPUT (ENGINE DECK): COL C C XM MACH NUMBER 1- 5 C C ALT ALTITUDE, FT 6-15 C C THRUST NET THRUST, LB 21-30 MINUS 31-40 C C A9 NOZZLE EXIT AREA, IN*IN 71-80 C C INSDRG INSTALLATION DRAG OPTION C C =0, NO INSTALLATION DRAG C C =1, SCALE INSTALLATION DRAG FOR CHANGES IN A10 C C =2, CALCULATE INSTALLATION DRAG; BASED ON A10 C C =3, CALCULATE INSTALLATION DRAG; CD=0 @ A9=A9REF C C C C----------------------------------------------------------------------C CHARACTER*80 CDFILE COMMON /UNITS / IU5, IU6, IU7, IU8, IU9, IU16, IU17, IU18 COMMON /ICDTAB/ ICDDAT(10), NTBL COMMON /CDFIL / CDFILE, A10, A9REF, A10REF, XNOZ, XNREF, RCRV, 1 NAB, NABREF COMMON /TRNSF / ASTAR ,CDT ,CL ,D ,DCD ,DIST , 1 DELTA ,ENG ,FUEL ,G ,HP ,HPL ,DYNPR , 2 RCI ,RCM ,RSM ,RSR ,SR ,SREF ,THETA , 3 TIME ,VH ,VT ,WF ,WFD ,WT ,DTC , 4 FACT ,ENGO ,VMIN , 5 ELD ,AMCH ,NW ,IFLAG ,IDOQ ,II SAVE INIT, LOCAB, LOCREF PARAMETER ( PI = 3.1415927 ) DATA INIT, LOCAB, LOCREF /3*0/ C----------------------------------------------------------------------C IF ( INSDRG .EQ. 0 ) RETURN C ... READ THE TABLE OF DRAG COEFFICIENTS AND INITIALIZE ARRAY LOCS. NPRINT = 0 IF ( INIT .EQ. 0 ) THEN IF ( A10*A10REF*XNOZ*XNREF .LE. 0. ) GO TO 2000 CALL CDRD ( NPRINT ) DO 10 N = 1, NTBL IF ( ICDDAT(N) .EQ. NAB ) LOCAB = N IF ( ICDDAT(N) .EQ. NABREF ) LOCREF = N 10 CONTINUE INIT = 1 ENDIF C ... GET THE DYNAMIC PRESSURE CALL ATMO ( ALT, DTC, DELTA, THETA, ASTAR, TM, RE, HFT ) Q = 2116.22 * 0.7 * XM * XM * DELTA C ... COMPUTE THE AFTERBODY DRAG GO TO ( 100, 200, 300 ) INSDRG GO TO 1000 C ... SCALE INITIAL BOATTAIL DRAG FOR A CHANGE IN A10 C----------------------------------------------------------------------C 100 CONTINUE IF ( LOCAB .NE. 0 ) THEN C ...... COMPUTE REFERENCE NOZZLE BOATTAIL ANGLE. R9REF = SQRT ( A9 / PI ) R10REF = SQRT ( A10REF / PI ) XLREF = 0.78 * XNREF BREF = ATAN ( ( R10REF - R9REF ) / XLREF ) * 180. / PI C ...... COMPUTE NOZZLE BOATTAIL ANGLE. R9 = R9REF R10 = SQRT ( A10 / PI ) XL = 0.78 * XNOZ BETA = ATAN ( ( R10 - R9 ) / XL ) * 180. / PI C ...... GET THE REFERENCE DRAG COEFFICIENT A9A10R = A9 / A10REF CALL CDLU ( LOCAB, XM, BREF, A9A10R, CDREF ) C ...... GET THE NEW DRAG COEFFICIENT A9A10 = A9 / A10 CALL CDLU ( LOCAB, XM, BETA, A9A10, CD ) C ...... UPDATE THE THRUST THRUST = THRUST + Q * ( CDREF * A10REF - CD * A10 ) / 144.0 ENDIF GO TO 1000 C ... COMPUTE THROTTLE SENSATIVE BOATTAIL DRAG; BASED ON A10 C----------------------------------------------------------------------C 200 CONTINUE IF ( LOCAB .NE. 0 ) THEN C ...... COMPUTE ACTUAL NOZZLE BOATTAIL ANGLE. R9 = SQRT ( A9 / PI ) R10 = SQRT ( A10 / PI ) XL = 0.78 * XNOZ BETA = ATAN ( ( R10 - R9 ) / XL ) * 180. / PI C ...... GET THE DRAG COEFFICIENT A9A10 = A9 / A10 CALL CDLU ( LOCAB, XM, BETA, A9A10, CD ) C ...... UPDATE THE THRUST THRUST = THRUST - Q * CD * A10 / 144.0 ENDIF GO TO 1000 C ... COMPUTE THROTTLE SENSATIVE BOATTAIL DRAG; CD=0 @ A9=A9REF C----------------------------------------------------------------------C 300 CONTINUE IF ( LOCAB .NE. 0 .AND. LOCREF .NE. 0 ) THEN C ...... COMPUTE REFERENCE NOZZLE BOATTAIL ANGLE. THIS IS FOR C THE NACELLE USED TO COMPUTE THE TOTAL AIRCRAFT DRAG. R9REF = SQRT ( A9REF / PI ) R10REF = SQRT ( A10REF / PI ) XLREF = 0.78 * XNREF BREF = ATAN ( ( R10REF - R9REF ) / XLREF ) * 180. / PI C ...... GET THE REFERENCE DRAG COEFFICIENT. THIS PORTION IS C ACCOUNTED FOR IN THE CONFIGURATION AERODYNAMICS. A9A10R = A9REF / A10REF CALL CDLU ( LOCREF, XM, BREF, A9A10R, CDAERO ) C ...... COMPUTE ACTUAL NOZZLE BOATTAIL ANGLE. R9 = SQRT ( A9 / PI ) R10 = SQRT ( A10 / PI ) XL = 0.78 * XNOZ BETA = ATAN ( ( R10 - R9 ) / XL ) * 180. / PI C ...... GET THE ACTUAL TOTAL DRAG COEFFICIENT A9A10 = A9 / A10 CALL CDLU ( LOCAB, XM, BETA, A9A10, CDACT ) C ...... GET THE ACTUAL THROTTLE SENSATIVE DRAG COEFFICIENT CDTS = CDACT * A10 - CDAERO * A10REF C ...... UPDATE THE THRUST THRUST = THRUST - Q * CDTS / 144.0 ENDIF GO TO 1000 C ... GO HERE TO RETURN 1000 RETURN C ... INPUT ERROR 2000 WRITE(IU6, 3000) 3000 FORMAT ( //'ILLEGAL VALUE FOR NOZZLE AREA OR LENGTH') STOP END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEPARATOR CCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE CDLU ( II, X, Y, Z, FXYZ ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) C*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=C C THIS SUBROUTINE DOES TABLE LOOK UP USING LINEAR INTERPOLATION C C OF DRAG COEFFICIENTS. C C C C ARGUMENTS C C II TABLE NUMBER FOR INTERPOLATION C C X,Y,Z INDEPENDENT VARIABLES C C FXYZ INTERPOLATED FUNCTION VALUE C C C C TABLE DEFINITION C C CDNOZ VECTOR CONTAINING ALL TABLES C C LOCCD(II) LOCATION OF TABLE II FIRST ELEMENT IN CDNOZ VECTOR C C C C ELEMENT DEFINITION FOR TABLE II IN CDNOZ VECTOR C C INZ = NUMBER OF Z'S (FIRST ELEMENT) C C INZ Z'S C C INZ+5 UNUSED ELEMENTS C C ............ C C INY = NUMBER OF Y'S . C C INY Y'S . C C INY+5 UNUSED ELEMENTS . REPEAT C C ......... . INZ TIMES C C INX = NUMBER OF X'S . . C C INX X'S . REPEAT . C C INX FUNCTION VALUES . INY TIMES . C C 5 UNUSED ELEMENTS ......... ............ C C C C THE UNUSED ELEMENTS WERE USED IN THE CUBIC SPLINE INTERPOLATION C C AND HAVE BEEN RETAINED FOR COMPATIBILITY C C*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=C COMMON /CDSIZE/ CDNOZ(5000) COMMON /CDLOC / LOCCD(5) DIMENSION FY(2), FX(2) C======================================================================C IF (II .LE. 0) RETURN C FIND START OF TABLE, NUMBER OF Z'S, AND START OF Y'S IZ = LOCCD(II) INZ = CDNOZ(IZ) IY = IZ + 2*INZ + 6 C IF INCOMING Z = FIRST Z VALUE, PRETEND IT'S THE ONLY ONE IF ( ABS(Z-CDNOZ(IZ+1)) .LT. .00001 ) INZ = 1 MZ = 0 C LOOP ON NUMBER OF Z'S IN TABLE DO 70 IZS = 1,INZ C FIND NUMBER OF Y'S AND START OF FIRST SET OF X'S INY = CDNOZ(IY) IX = IY + 2*INY + 6 C FIND IF Z NEEDS EXTRAPOLATION OR HAS BEEN BOUNDED IF ( IZS .GE. INZ-1 .OR. CDNOZ(IZ+IZS+1) .GT. Z ) MZ = MZ + 1 IF ( MZ .EQ. 1 ) JZ = IZS MY = 0 C LOOP ON NUMBER OF Y'S IN TABLE FOR THIS Z DO 50 IYS = 1,INY C FIND NUMBER OF X'S INX = CDNOZ(IX) C DETERMINE IF CURRENT Y IS A BOUNDING VALUE IF ( MY .EQ. 2 .OR. MZ .EQ. 0 ) GO TO 40 IF ( MY .EQ. 0 .AND. CDNOZ(IY+IYS+1) .LT. Y 1 .AND. IYS .LT. INY-1 ) GO TO 40 MY = MY + 1 IF ( MY .EQ. 1 ) JY = IYS IF ( INX .GT. 1 ) GO TO 10 C ONLY ONE X USED FX(MY) = CDNOZ(IX+2) GO TO 40 C LOOP ON NUMBER OF X'S TO FIND BOUNDING VALUES 10 DO 20 IXS = 2,INX IF ( CDNOZ(IX+IXS) .GT. X ) GO TO 30 20 CONTINUE C INTERPOLATE ON X IXS = INX 30 DX = (X - CDNOZ(IX+IXS-1)) / (CDNOZ(IX+IXS) - CDNOZ(IX+IXS-1)) FX(MY) = (1.-DX)*CDNOZ(IX+IXS-1+INX) + DX*CDNOZ(IX+IXS+INX) C UPDATE START OF NEXT SET OF X'S 40 IX = IX + 2*INX + 6 50 CONTINUE IF ( MZ .EQ. 0 ) GO TO 60 C ONLY ONE Y USED FY(MZ) = FX(1) IF ( MY .EQ. 1 ) GO TO 60 C INTERPOLATE ON Y DY = (Y - CDNOZ(JY+IY)) / (CDNOZ(JY+IY+1) - CDNOZ(JY+IY)) FY(MZ) = (1. - DY) * FX(1) + DY * FX(2) C UPDATE START OF NEXT SET OF Y'S 60 IY = IX C CHECK TO SEE IF FINISHED IF ( MZ .GE. 2 ) GO TO 80 70 CONTINUE C ONLY ONE Z USED FXYZ = FY(1) RETURN C INTERPOLATE ON Z 80 DZ = (Z - CDNOZ(JZ+IZ)) / (CDNOZ(JZ+IZ+1) - CDNOZ(JZ+IZ)) FXYZ = (1. - DZ) * FY(1) + DZ * FY(2) RETURN END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEPARATOR CCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE CDRD ( II ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) C*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=C C THIS SUBROUTINE READS PROPULSION SYSTEM DRAG DATA. C C C C USE NON-ZERO "II" TO WRITE TABLE DATA. C C*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=C CHARACTER*80 TITLE, CDFILE DIMENSION A(100),IREP(30),IDD(4) COMMON /UNITS / IU5, IU6, IU7, IU8, IU9, IU16, IU17, IU18 COMMON /CYUNIT/ IU3, IU4, IU10, IU11, IU12, IU13, IU14, IU15 COMMON /CYICOM/ IENG, IPRINT, NPRINT, NPDRY, NPAB, NITMAX, 1 IVAT, LIMCD COMMON /CDSIZE/ CDNOZ(5000) COMMON /ICDTAB/ ICDDAT(10), NTBL COMMON /CDLOC/ LOCCD(5) COMMON /CDFIL / CDFILE, A10, A9REF, A10REF, XNOZ, XNREF, RCRV, 1 NAB, NABREF C======================================================================C IU12 = 12 IUNIT = IU12 OPEN (UNIT=IUNIT, FILE=CDFILE, STATUS='OLD', ERR=9992) REWIND (IUNIT) NTBL = 0 NMAX = 5000 MR = 0 IF ( NTBL .EQ. 0 ) THEN LOC = 1 DO 10 M = 1,NMAX 10 CDNOZ(M) = 0. ENDIF 20 READ (IU12, 5000, END=400, ERR=6666 ) IP, ITABNO, TITLE 30 IF(ITABNO .EQ. 0) THEN NREM = NMAX-LOC IF ( II .GT. 0 ) THEN WRITE(IU6,5010) NTBL WRITE(IU6,5020) (N, ICDDAT(N), LOCCD(N), N = 1,NTBL) WRITE(IU6,5030) NMAX, NREM ENDIF IF( MR .GT. 0 ) WRITE(IU6,5040) (IREP(IQ), IQ = 1,MR) RETURN ELSE NTBL = NTBL+1 ICDDAT(NTBL) = ITABNO LOCCD(NTBL) = LOC IF ( NTBL .GT. 1 ) THEN DO 40 I = 2,NTBL IF ( ITABNO .EQ. ICDDAT(I-1) ) GO TO 50 40 CONTINUE GO TO 150 50 MLOC = LOCCD(I-1) ILOC = MLOC LOCCD(I-1) = LOC NTBL = NTBL - 1 MR = MR + 1 IREP(MR) = ITABNO LIBLOC = 100000 DO 60 IRK = 1,NTBL LIBRA = LOCCD(IRK) IF ( (LIBRA .GE. MLOC) .AND. (LIBRA .LT. LIBLOC) ) * LIBLOC = LIBRA 60 CONTINUE LEO = LIBLOC - MLOC IF ( LEO .NE. 0 ) THEN IZAR = LOC - 1 DO 70 I = LIBLOC,IZAR CDNOZ(MLOC) = CDNOZ(I) MLOC = MLOC + 1 70 CONTINUE DO 80 I = 1,NTBL LLOC = LOCCD(I) IF ( LLOC .GT. ILOC ) LOCCD(I) = LLOC - LEO 80 CONTINUE ENDIF DO 90 I = MLOC,LOC 90 CDNOZ(I) = 0. LOC = MLOC ENDIF ENDIF 150 IF ( II .GT. 0 ) THEN WRITE ( IU6, 5060 ) ITABNO, TITLE IP = 1 ENDIF ICC = 0 LZ = LOC DO 180 ICL = 1,4 180 IDD(ICL) = 4H 190 IDL = ID READ (IU12, 5070, END=200, ERR=6666 ) ID, N, (A(I), I = 1,N) 200 IF ( ID .EQ. 4HEOT ) GO TO 20 ICC = ICC + 1 IF ( ICC .LE. 4 ) IDD(ICC) = ID IF ( ID .EQ. IDL ) GO TO 280 IF ( ID .NE. IDD(4) ) GO TO 210 IF ( IDL .EQ. IDD(2) ) GO TO 280 LOC = LE L = LX GO TO 220 210 CDNOZ(LOC) = N L = LOC IF ( ID .EQ. IDD(3) ) LX = L IF ( ID .NE. IDD(2) ) GO TO 260 LY = L LZ = LZ + 1 GO TO 260 C ... WRITE THE TABULAR DATA. 220 IF ( IP .NE. 0 ) THEN LY = LY+1 WRITE(IU6,5080) IDD(1), CDNOZ(LZ), IDD(2), CDNOZ(LY) JF = 0 M = N LF = LX 240 NP = M IF ( M .GT. 8 ) NP = 8 M = M - NP LE = LF + NP LF = LF + 1 WRITE(IU6,5090) IDD(3), (CDNOZ(I), I = LF,LE) LF = LE JE = JF + NP JF = JF + 1 WRITE(IU6,5090) IDD(4), (A(I), I = JF,JE) JF = JE IF( M .GT. 0 ) GO TO 240 ENDIF 260 DO 270 I = 1,N LOC = LOC + 1 CDNOZ(LOC) = A(I) 270 CONTINUE LE = LOC IF ( ID .EQ. IDD(4) ) CDNOZ(LOC+2) = 1. LOC = L + 2*N + 6 IF ( LOC .GT. NMAX ) GO TO 300 GO TO 190 280 CDNOZ(LOC) = CDNOZ(LX) L = LOC DO 290 I = 1,N LOC = LOC + 1 CDNOZ(LOC) = CDNOZ(LX+I) 290 CONTINUE GO TO 220 300 WRITE(IU6,5100) ITABNO ITABNO = 0 GO TO 30 400 CONTINUE ITABNO = 0 GO TO 30 6666 CONTINUE WRITE(IU6,6669) IU12 STOP 5000 FORMAT (I1,I4,A) 5010 FORMAT (//,33X,'TABLE DATA INPUT SUMMARY, ',I3,' TABLES',//, 1 28X,'TABLE NUMBER REFERENCE NUMBER ARRAY LOCATION') 5020 FORMAT (33X,I2,12X,I5,14X,I5) 5030 FORMAT (/,36X,'DATA STORAGE ALLOCATION ',I7,/, 1 36X,'DATA STORAGE NOT USED ',I7,/) 5040 FORMAT(10X,'THE FOLLOWING TABLES HAVE BEEN REPLACED',10(I5,',')) 5060 FORMAT (//,3X,I4,A,/) 5070 FORMAT(A4,I3,3X,(T11,7F10.0)) 5080 FORMAT(1X,A4,' = ',E13.5,1X,A4,' = ',E13.5) 5090 FORMAT(20X,A4,1X,8E13.5) 5100 FORMAT (' ********* TABLE OVER FLOW, TABLE ',I5,' NOT LOADED') 6669 FORMAT(/,' ERROR READING ENGINE TABULAR INPUT DATA FROM UNIT', * I3,'.',/,' PROGRAM ABORTED IN SUBROUTINE TABRD.') 9992 CONTINUE WRITE(IU6,9999) 'ERROR OPENING THE INPUT FILE ', CDFILE, * ' AS UNIT ', IUNIT, '.', * 'PROGRAM ABORTED IN SUBROUTINE CDRD.' 9999 FORMAT(//,1X,2A,/,A,I2,A,/,1X,A,/) STOP END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEPARATOR CCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE PLTTHR ( IST ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c ohad 15/7/08 c IMPLICIT INTEGER*8 (I-N) IMPLICIT INTEGER*4 (I-N) C PREPARES FULL PLOT FILE FOR ENTIRE ENGINE DECK C IST CALLING ROUTINE INDICATOR C = 1, FROM DEFENG C = 2, FROM MAIN PROGRAM C = 3, PLOT FILE PREVIOUSLY PREPARED COMMON /UNITS / IU5, IU6, IU7, IU8, IU9, IU16, IU17, IU18 COMMON /THRPLT/ PCODE(16), NPCODE, IPLTTH COMMON /NOXDAT/ OX(16,15,20), FFUEL(6), FNOX(6), RTNOX, TNOX, 1 GNOX, OFNOX, NOX COMMON /ENDTA / EXTFAC,FFFAC ,DFFAC ,TAKOFF,EMACH(20) , 1 ALT(15,20) ,FF(16,15,20) ,THR(16,15,20), 2 TXFUFL,NM ,NA ,NT ,IXTRAP,MAXCR IF ( IPLTTH .NE. IST ) RETURN WRITE (IU17, 10) 10 FORMAT (' MACH ALTITUDE PC THRUST RAM DRAG FUEL FLOW', 1 ' SFC NOX EM') DO 100 I = 1,NM DO 90 J = 1,NA IF ( FF(1,J,I) .LE. 0. ) GO TO 100 DO 80 K = 1,NT FFF = FF(K,J,I) IF ( FFF .LE. 0. ) GO TO 90 SFC = 0. IF (THR(K,J,I) .GT. THR(1,J,I)/100.) SFC = FFF/THR(K,J,I) WRITE (IU17, 50) EMACH(I), ALT(J,I), PCODE(K), 1 THR(K,J,I), FFF, SFC, OX(K,J,I) 50 FORMAT(F5.2,F10.0,F5.0,F10.1,10X,F10.1,1X,F9.5,F10.3) 80 CONTINUE 90 CONTINUE 100 CONTINUE CLOSE(IU17) IST = 3 RETURN END