program tsfoil2 C PROGRAM TSFOIL(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT,TAPE3,TAPE7) C*********************************************************************** C C MAIN PROGRAM FOR TSFOIL C PROGRAM COMPUTES TRANSONIC FLOW PAST A TWODIMENSIONAL C LIFTING AIRFOIL USING TRANSONIC SMALL DISTURBANCE THE C C PROGRAM WRITTEN BY C EARLL M. MURMAN AND FRANK R. BAILEY C NASA-AMES RESEARCH CENTER C AND C MARGARET L. JOHNSON C COMPUTER SCIENCES CORPORATION C version for VT Aerospace Design Space software C C REFERENCES FOR PROGRAM USAGE ARE C ( TO BE FILLED IN LATER) C c from Joh, VPI, with modifications c by w.h. mason for mac operation c and with diamond airfoil option c c mod to I/O to fix bug, April 2, 1998 c mod to output file for plotting, April 3, 2000 c mod to run on PC C - Added directory selection C - Added Pause statement at end of program C - Added more comments in the header C - Program outputs to a file now C Andy Ko 4/2/03. C Fixed bug in Jameson input section 4/9/03 - Andy Ko c C*********************************************************************** c USE DFPORT COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK LOGICAL ABORT1 COMMON / COM3/ IREF , ABORT1 , ICUT , KSTEP LOGICAL AMESH COMMON / COM4/ XIN(100) , YIN(100), AMESH COMMON / COM5/ XDIFF(100),YDIFF(100) COMMON / COM6/ FL(100) , FXL(100), FU(100) , FXU(100), 1 CAMBER(100), THICK(100),VOL , XFOIL(100), IFOIL COMMON / COM7/ CJUP , CJUP1 , CJLOW , CJLOW1 COMMON / COM8/ CVERGE , DVERGE , IPRTER , MAXIT , 1 WE(3) , EPS INTEGER BCFOIL COMMON / COM9/ BCFOIL , NL , NU , XL(100) , XU(100) , 1 YL(100) , YU(100) , RIGF,IFLAP,DELFLP,FLPLOC COMMON /COM10/ YFREE(100),YTUN(100),GAM , JMXF , JMXT INTEGER PSTART LOGICAL PSAVE COMMON /COM11/ ALPHAO , CLOLD , DELTAO , DUBO , EMACHO , 1 IMINO , IMAXO , IMAXI , JMINO , JMAXO , 2 JMAXI , PSAVE , PSTART ,TITLE(20),TITLEO(20), 3 VOLO , XOLD(100),YOLD(100) COMMON /COM12/ F , H , HALFPI , PI , RTKPOR , 1 TWOPI COMMON /COM13/ CDFACT , CLFACT , CMFACT , CPFACT , CPSTAR LOGICAL FCR , KUTTA COMMON /COM14/ CLSET , FCR , KUTTA , WCIRC COMMON /COM15/ B , BETA0 , BETA1 , BETA2 , PSI0 , 1 PSI1 , PSI2 REAL JET COMMON /COM16/ ALPHA0 , ALPHA1 , ALPHA2 , XSING , 1 OMEGA0 , OMEGA1 , OMEGA2 , JET COMMON /COM17/ CYYBLC , CYYBLD , CYYBLU , CYYBUC , CYYBUD , 1 CYYBUU ,FXLBC(100),FXUBC(100), ITEMP1, ITEMP2 LOGICAL OUTERR COMMON /COM18/ ERROR , I1 , I2 , IERROR , JERROR , 1 OUTERR , EMU(100,2) , VC(100) , 2 WI , DCIRC , POLD(100,2) COMMON /COM19/ DIAG(100), RHS(100), SUB(100), SUP(100) COMMON /COM20/ XMID(100), YMID(100) COMMON /COM21/ CYYLC , CYYLD , CYYLU , CYYUC , CYYUD , 1 CYYUU COMMON /COM22/ CXC(100) , CXL(100), CXR(100), CXXC(100),CXXL(100), 1 CXXR(100), C1(100) COMMON /COM23/ CYYC(100), CYYD(100),CYYU(100), IVAL COMMON /COM24/ DTOP(100), DBOT(100),DUP(100), DDOWN(100), 1 VTOP(100), VBOT(100),VUP(100), VDOWN(100) COMMON /COM25/ CPL(100) , CPU(100) , IDLA COMMON /COM26/ PJUMP(100) LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT INTEGER BCTYPE COMMON /COM28/ BCTYPE , CIRCFF , FHINV , POR , CIRCTE COMMON /COM32/ BIGRL , IRL , JRL COMMON /COM33/ THETA(100,100) COMMON /COM34/ NWDGE , WSLP(100,2) , XSHK(2,3) , 1 THAMAX(2,3) , AM1(2,3), ZETA(2,3) , 2 NVWPRT(2), WCONST , REYNLD , NISHK CHARACTER*50 newdir CHARACTER*50 output INTEGER*4 istatus C WRITE HEADER PAGE INFORMATION C C WRITE (6,1000) WRITE (6,1001) WRITE (6,1002) WRITE (6,1003) WRITE (6,1002) WRITE (6,1004) WRITE (6,1002) WRITE (6,1005) WRITE (6,1002) WRITE (6,1006) WRITE (6,1002) WRITE (6,1001) 1000 FORMAT(1H1/////////////////) 1001 FORMAT(4X,69(1H*)) 1002 FORMAT(4X,1H*,67X,1H*) 1003 FORMAT(4X,1H*,25X,15H PROGRAM TSFOIL,27X,1H*/ 1 4X,1H*,29X, 7H SOLVES,31X,1H*/ 2 4X,1H*, 5X,'INVISCID FLOW PAST THIN TWO DIMENSIONAL LIFTING', 3 ' AIRFOIL',7X,1H*/ 4 4X,1H*,29X, 6H USING,32X,1H*/ 5 4X,1H*,15X,35H TRANSONIC SMALL DISTURBANCE THEORY,17X,1H*/ 6 4X,1H*, 9X,'FULLY CONSERVATIVE FINITE DIFFERENCE EQUATIONS', 7 12X,1H*/ 8 4X,1H*,17X,31H SUCCESSIVE LINE OVERRELAXATION,19X,1H*) 1004 FORMAT(4X,1H*,27X,11H WRITTEN BY,29X,1H*) 1005 FORMAT(4X,1H*,15X,36H EARLL M. MURMAN AND FRANK R. BAILEY,16X,1H* 1 / 4X,1H*,19X,26H NASA-AMES RESEARCH CENTER,22X,1H*/ 2 4X,1H*,19X,26H MOFFETT FIELD, CALIFORNIA,22X,1H*/ 3 4X,1H*,31X, 4H AND,32X,1H*/ 4 4X,1H*,23X,20H MARGARET L. JOHNSON,24X,1H*/ 5 4X,1H*,18X,30H COMPUTER SCIENCES CORPORATION,19X,1H*/ 6 4X,1H*,20X,26H MOUNTAIN VIEW, CALIFORNIA,21X,1H*/ 7 4X,1H*,67X,1H*/ 8 4X,1H*,12X,43H Documented in NASA SP-347 and NASA CR-3064 9 ,12X,1H*) 1006 FORMAT(4X,1H*,27X,12H Version for, 28X,1H*/ 1 4X,1H*,15X,36H VT Aerospace Design Software Series,16X,1H*/ 2 4X,1H*,67X,1H*/ 2 4X,1H*,15X,14H Contact info:,38X,1H*/ 3 4X,1H*,15X,21H Dr. William H. Mason,31X,1H*/ 3 4X,1H*,15X,45H Aerospace & Ocean Engineering, Virginia Tech 4 ,7X,1H*/ 5 4X,1H*,15X, 20H Blacksburg, VA24061, 32X,1H*/ 6 4X,1H*,15X, 22H Email: whmason@vt.edu,30X,1H*) C Set the input/output directory. Default is C:\ - Andy Ko 4/2/03 istatus = CHDIR('C:\') IF (istatus .eq. enoent) THEN WRITE(*,*) 'Cannot change directory to C:\' WRITE(*,*) 'Please log into the network' Pause 'Program ended.Press the ENTER key to exit' STOP ENDIF WRITE (*,*) 'Default directory is C:\' WRITE (*,*) 'Please enter a new directory name. ' WRITE (*,*) 'Type "default" for the default directory' READ (*,*) newdir IF (newdir .eq. 'default') THEN WRITE (*,*) 'Default directory selected' ELSE istatus = CHDIR(newdir) IF (istatus .eq. enoent) THEN WRITE(*,*) 'The directory does not exist' WRITE(*,*) 'The directory is set to the default' ELSEIF (istatus .eq. enotdir) THEN WRITE(*,*) 'Directory could not be changed' WRITE(*,*) 'The directory is set to the default' ELSE WRITE(*,*) 'directory has been changed to ',newdir ENDIF ENDIF C WRITE OUTPUT TO FILE WRITE (*,*) 'Enter the name of the output file' READ (*,*) output OPEN (UNIT=15,FILE=output) WRITE (15,1001) WRITE (15,1002) WRITE (15,1003) WRITE (15,1002) WRITE (15,1004) WRITE (15,1002) WRITE (15,1005) WRITE (15,1002) WRITE (15,1006) WRITE (15,1002) WRITE (15,1001) C C THE MAIN PROGRAM DOES NO COMPUTATIONS C ALL COMPUTATIONS ARE DONE IN SUBROUTINES CALLED BY C TSFOIL. C SUBROUTINE ECHINP PROVIDES A LISTING OF ALL DATA C CARDS FOR ENTIRE JOB. CAN BE DELETED, IF DESIRED. C CALL ECHINP C SUBROUTINE READIN READS ALL INPUT AND CHECKS IT 1 CONTINUE CALL READIN C C SUBROUTNIE SCALE RESCALES ALL PHYSICAL VARIABLES TO C TRANSONIC SIMILARITY FORM CALL SCALE C C SUBROUTINE FARFLD SETS FAR FIELD BOUNDARY CONDITIONS. CALL FARFLD C C SUBROUTINE BODY COMPUTES AIRFOIL GEOMETRY AND PRINTS C OUT GEOMETRICAL INFORMATION CALL BODY C C SUBROUTINE CUTOUT REMOVES MESH POINTS FROM THE INPUT C MESH. CALCULATIONS ARE DONE FIRST ON COARSE MESH, C AND THEN ON PROGRESSIVELY REFINED MESHES UNTIL C INPUT MESH IS ACHIEVED. CALL CUTOUT C C SUBROUTINE GUESSP INITIALIZES P ARRAY CALL GUESSP C C SUBROUTINE DIFCOE CALCULATES FINITE DIFFERENCE C EQUATION COEFFICIENTS WHICH DEPEND ON MESH SIZE. CALL DIFCOE C C SUBROUTINE SETBC ADJUSTS THE AIRFOIL SLOPE BOUNDARY C CONDITION FOR THE CURRENT X Y MESH. ALSO, THE LIMITS C ON I AND J INDICIES FOR SOLVING THE DIFFERENCE C EQUATIONS ARE SET. CALL SETBC(0) C C SUBROUTINE SOLVE EXECUTES THE MAIN RELAXATION C SOLUTION OF THE DIFFERENCE EQUATIONS CALL SOLVE C C IF FINAL MESH HAS BEEN REACHED, RESULTS ARE PRINTED C OUT IN FINAL FORM. IF NOT, INTERMEDIATE RESULTS C ARE PRINTED OUT AND THE MESH IS REFINED. THE ABOVE C SEQUENCE OF CALCULATIONS ARE THEN REPEATED IF(IREF .LE. 0) GO TO 5 IF (ABORT1) GO TO 5 C C SUBROUTINE PRINT1 PRINTS OUT BODY PRESSURE C DISTRIBUTION. CALL PRINT1 IF (ABORT1) GO TO 1 C C SUBROUTINE REFINE ADDS MESH POINTS CALL REFINE C C REPEAT SEQUENCE OF RELAXATION CALCULATIONS CALL DIFCOE CALL SETBC(0) CALL SOLVE C C IF FINAL MESH HAS BEEN REACHED, RESULTS ARE PRINTED C OUT IN FINAL FORM. IF NOT, INTERMEDIATE RESULTS C ARE PRINTED OUT AND THE MESH IS REFINED. THE ABOVE C SEQUENCE OF CALCULATIONS ARE THEN REPEATED IF(IREF .LE. 0) GO TO 5 IF (ABORT1) GO TO 5 CALL PRINT1 CALL REFINE CALL DIFCOE CALL SETBC(0) CALL SOLVE C C RELAXATION SOLUTION IS COMPLETED 5 CONTINUE C C PRINTOUT FINAL INFORMATION CALL PRINT IF (IREF .GT. 0 ) CALL REFINE IF (IREF .GT. 0 ) CALL REFINE IF (IDLA.NE.0) CALL DLAOUT(ILE,ITE,ALPHA,DELFLP,EMACH, & VFACT,REYNLD) C C STORE SOLUTION FOR NEXT CASE OR ON TAPE 3 CALL SAVEP C C RETURN TO READIN TO COMPUTE NEXT CASE OR TERMINATE C CALCULATIONS. cwhm GO TO 1 CLOSE(UNIT=15) Pause 'Press the ENTER key to exit' stop END SUBROUTINE ANGLE C COMPUTES THETA AT EVERY MESH POINT. C CALLED BY - FARFLD. COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK LOGICAL AMESH COMMON / COM4/ XIN(100) , YIN(100), AMESH COMMON /COM12/ F , H , HALFPI , PI , RTKPOR , 1 TWOPI REAL JET COMMON /COM16/ ALPHA0 , ALPHA1 , ALPHA2 , XSING , 1 OMEGA0 , OMEGA1 , OMEGA2 , JET COMMON /COM33/ THETA(100,100) C C SUBROUTINE TO COMPUTE THE ANGLE THETA AT C EACH MESH POINT. C R2PI = 1.0 / TWOPI DO 20 I=IMIN,IMAX XX = XIN(I) - XSING DO 10 J=JMIN,JMAX YY = YIN(J) * RTK R = SQRT(YIN(J)**2 + XX*XX) ATN = ATAN2(YY,XX) Q = PI - SIGN(PI,YY) THETA(J,I) = -(ATN + Q) * R2PI IF ( R .LE. 1.0 ) THETA(J,I) = THETA(J,I) * R 10 CONTINUE 20 CONTINUE RETURN END FUNCTION ARF (X) C C EVALUATES ERF WITH AN ERROR .LT. 1.5E-7 BY RATIONAL C APPROXIMATION 7.1.26 FO HANDBOOK OF MATH. FUNCTIONS C U. S. DEPT. OF COMMERCE. NBS APPL MATH SER 55. C CALLED BY - AYMESH. C DIMENSION C(5) DATA C /1.061405429,-1.453152027,1.421413741, -.284496736, 1 .254829592/ C Y = X IF ( X .LT. 0.0 ) Y = -Y IF ( Y .LT. 10. ) GO TO 10 ARF = 1.0 GO TO 30 10 CONTINUE T = 1.0 / (1.0 + .3275911*Y) POLY = 0.0 DO 20 I=1,5 POLY = (POLY + C(I)) * T 20 CONTINUE ARF = 1.0 - POLY * EXP(-Y*Y) 30 CONTINUE IF (X .LT. 0.0) ARF = -ARF RETURN END SUBROUTINE AYMESH C COMPUTES ANALYTICAL X AND Y MESH POINTS. C CALLED BY - READIN. LOGICAL AMESH COMMON / COM4/ XIN(100) , YIN(100), AMESH INTEGER PSTART LOGICAL PSAVE COMMON /COM11/ ALPHAO , CLOLD , DELTAO , DUBO , EMACHO , 1 IMINO , IMAXO , IMAXI , JMINO , JMAXO , 2 JMAXI , PSAVE , PSTART ,TITLE(20),TITLEO(20), 3 VOLO , XOLD(100),YOLD(100) COMMON /COM12/ F , H , HALFPI , PI , RTKPOR , 1 TWOPI INTEGER BCTYPE COMMON /COM28/ BCTYPE , CIRCFF , FHINV , POR , CIRCTE COMMON /COM30/ BXH(401) , REST(3) C DIMENSION XH(401) , BX(100) DATA A0/.225/, A1/1.4/, A2/1.6/, A3/.6188/, A4/.75/, A5/30./, 1 A6/.603/, A7/2.0/ DATA CT1/2.0/, CT2/2.0/, CT3/1.0/, CF1/1.0/, CF2/1.0/, CF3/5.2/ C IF (IMAXI .EQ. 77) IMAXI=81 IMA = (IMAXI - 1) / 2 A22 = A2 * A2 A72 = A7 * A7 FACT = HALFPI * .005 DO 10 I=201,400 EX = TAN(FACT * (I-201) ) XH(I) = EX EX2 = EX * EX EX22 = A22 * EX2 EX72 = A72 * EX2 IF (EX72 .LT. 173.0) GO TO 4 TEX7 = 0.0 IF (EX22 .LT. 173.0) GO TO 6 TEX2 = 0.0 GO TO 8 4 TEX7 = 1.0/EXP(EX72) 6 TEX2 = 1.0/EXP(EX22) 8 BXH(I) = A1*EX*TEX2 + (1.0-TEX7)*ARF(A4*EX) 10 CONTINUE XH(401) = 1.E30 BXH(401) = 1.0 DO 20 I=1,200 XH(I) = - XH(402-I) BXH(I) = -BXH(402-I) 20 CONTINUE TOOPI = 2.0 / PI C C ADD 2/PI ATAN(AK(X+A6)) TO BXH TO GIVE MORE C POINTS NEAR THE LEADING EDGE. C DO 30 I=1,401 BXH(I) = BXH(I) * (1.-A0) + TOOPI*A0*ATAN(A5*(XH(I)+A6)) 30 CONTINUE C DX = 1.0 / IMA IM2 = IMA * 2 IM2P1= IM2 + 1 C DO 42 I=1,IM2P1 IF (I .EQ. 1 .OR. I .EQ. IM2P1) GO TO 31 BX(I) = (I-1) * DX - 1.0 GO TO 32 31 CONTINUE IF (I .EQ. 1) BX(I) = -.999 IF (I .EQ. IM2P1) BX(I) = .999 32 CONTINUE J = 0 33 CONTINUE J = J + 1 IF (BX(I) .GT. BXH(J)) GO TO 33 IF (BX(I) .LT. BXH(J)) GO TO 34 XI = XH(J) GO TO 40 C 34 CONTINUE BT1 = BX(I) - BXH(J-1) BHT1 = BXH(J) - BXH(J-1) XHT1 = XH(J) - XH(J-1) T1 = XHT1 / BHT1 XI = XH(J-1) + BT1 * T1 IF (J .EQ. 2) GO TO 40 C BT2 = BX(I) - BXH(J) BHT2 = BXH(J) - BXH(J-2) BHT3 = BXH(J-1) - BXH(J-2) XHT2 = XH(J-1) - XH(J-2) T2 = XHT2 / BHT3 T12 = T1 - T2 BT12 = BT1 * BT2 TBT2 = T12 / BHT2 XI = XI + BT12 * TBT2 IF (J .GE. 400) GO TO 40 C BT3 = BX(I) - BXH(J-2) BHT4 = BXH(J+1) - BXH(J-2) BHT5 = BXH(J+1) - BXH(J) BHT6 = BXH(J+1) - BXH(J-1) XHT3 = XH(J+1) - XH(J) T3 = XHT3 / BHT5 XI = XI + BT12 * BT3/BHT4 * ((T3-T1)/BHT6 - TBT2) 40 CONTINUE XIN(I) = XI + A3 42 CONTINUE IMAXI = IM2P1 JM2 = JMAXI JMA = JM2 / 2 FJ = JMA IF (BCTYPE .NE. 1) GO TO 44 C1 = CF1 C2 = CF2 C3 = CF3 GO TO 45 44 CONTINUE C1 = CT1 C2 = CT2 C3 = CT3 45 CONTINUE DETA = 1.0 / (FJ*C1) IF (BCTYPE .EQ. 1) DETA = 1.0 / ((FJ+1.0)*C1) C = C3 / (TAN(HALFPI*DETA*FJ))**C2 DO 50 I=1,JMA J = JMA + I ETAI = I * DETA YIN(J) = C * (TAN(HALFPI*ETAI))**C2 YIN(J-2*I+1) = -YIN(J) 50 CONTINUE RETURN END BLOCK DATA C C AND ALL COMMONS WHETHER THEY ARE NEEDED OR NOT. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK LOGICAL ABORT1 COMMON / COM3/ IREF , ABORT1 , ICUT , KSTEP LOGICAL AMESH COMMON / COM4/ XIN(100) , YIN(100), AMESH COMMON / COM5/ XDIFF(100),YDIFF(100) COMMON / COM6/ FL(100) , FXL(100), FU(100) , FXU(100), 1 CAMBER(100), THICK(100),VOL , XFOIL(100), IFOIL COMMON / COM7/ CJUP , CJUP1 , CJLOW , CJLOW1 COMMON / COM8/ CVERGE , DVERGE , IPRTER , MAXIT , 1 WE(3) , EPS INTEGER BCFOIL COMMON / COM9/ BCFOIL , NL , NU , XL(100) , XU(100) , 1 YL(100) , YU(100) , RIGF,IFLAP,DELFLP,FLPLOC COMMON /COM10/ YFREE(100),YTUN(100),GAM , JMXF , JMXT INTEGER PSTART LOGICAL PSAVE COMMON /COM11/ ALPHAO , CLOLD , DELTAO , DUBO , EMACHO , 1 IMINO , IMAXO , IMAXI , JMINO , JMAXO , 2 JMAXI , PSAVE , PSTART ,TITLE(20),TITLEO(20), 3 VOLO , XOLD(100),YOLD(100) COMMON /COM12/ F , H , HALFPI , PI , RTKPOR , 1 TWOPI COMMON /COM13/ CDFACT , CLFACT , CMFACT , CPFACT , CPSTAR LOGICAL FCR , KUTTA COMMON /COM14/ CLSET , FCR , KUTTA , WCIRC COMMON /COM15/ B , BETA0 , BETA1 , BETA2 , PSI0 , 1 PSI1 , PSI2 REAL JET COMMON /COM16/ ALPHA0 , ALPHA1 , ALPHA2 , XSING , 1 OMEGA0 , OMEGA1 , OMEGA2 , JET COMMON /COM17/ CYYBLC , CYYBLD , CYYBLU , CYYBUC , CYYBUD , 1 CYYBUU ,FXLBC(100),FXUBC(100), ITEMP1, ITEMP2 LOGICAL OUTERR COMMON /COM18/ ERROR , I1 , I2 , IERROR , JERROR , 1 OUTERR , EMU(100,2) , VC(100) , 2 WI , DCIRC , POLD(100,2) COMMON /COM19/ DIAG(100), RHS(100), SUB(100), SUP(100) COMMON /COM20/ XMID(100), YMID(100) COMMON /COM21/ CYYLC , CYYLD , CYYLU , CYYUC , CYYUD , 1 CYYUU COMMON /COM22/ CXC(100) , CXL(100), CXR(100), CXXC(100),CXXL(100), 1 CXXR(100), C1(100) COMMON /COM23/ CYYC(100), CYYD(100),CYYU(100), IVAL COMMON /COM24/ DTOP(100), DBOT(100),DUP(100), DDOWN(100), 1 VTOP(100), VBOT(100),VUP(100), VDOWN(100) COMMON /COM25/ CPL(100) , CPU(100) , IDLA COMMON /COM26/ PJUMP(100) LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT INTEGER BCTYPE COMMON /COM28/ BCTYPE , CIRCFF , FHINV , POR , CIRCTE COMMON /COM32/ BIGRL , IRL , JRL COMMON /COM33/ THETA(100,100) COMMON /COM34/ NWDGE , WSLP(100,2) , XSHK(2,3) , 1 THAMAX(2,3) , AM1(2,3), ZETA(2,3) , 2 NVWPRT(2), WCONST , REYNLD , NISHK C DATA WSLP /200*0.0/ , NWDGE / 0 / , WCONST / 4.0 / , 1 REYNLD / 4.0E+06 / , ZETA /6*0.0/ , IDLA / 0 / DATA IFLAP / 0 /, DELFLP / 5.0 /, FLPLOC / 0.77 / DATA BCFOIL / 3 / , BCTYPE / 1 / , PSTART / 1 / , 1 PRTFLO / 1 / , SIMDEF / 3 / DATA PHYS /.TRUE. / , PSAVE / .FALSE. / , FCR / .TRUE. / , 1 KUTTA / .TRUE. /, ABORT1 / .TRUE. /, AMESH / .FALSE. / DATA EMACH /.75/ , DELTA /.115/ , ALPHA /.12/ , 1 AK /0.0/, GAM/1.4/, RIGF/0.0/, EPS/.2/ DATA CLSET /0.0/, CVERGE/.00001/, DVERGE/10./, 1 F / 0.0 /, H / 0.0 / , POR / 0.0/, 1 WCIRC/1.0/, WE/1.8,1.9,1.95/ DATA IMAXI/77/, JMXF/56/, MAXIT/500/, 1 NL / 75 /, NU / 100 / , IPRTER / 10 / DATA PI /3.14159265/ , HALFPI /1.570796325/ , 1 TWOPI/ 6.28318531/ DATA IMIN /1/ , JMIN /1/ , ICUT /2/ DATA YIN /100*0.0 /, JMAXI / 64 / , JMXT /48 / DATA XIN / -1.075 , -.950 , 1 -.825 , -.7 , -.575 , -.45 , -.35 , 2 -.25 , -.175 , -.125 , -.075 , -.0525 , 3 -.035 , -.0225 , -.015 , -.0075 , -.0025 , 4 .0025 , .0075 , .0125 , .0175 , .0225 , 5 .0275 , .0325 , .0375 , .045 , .055 , 6 .065 , .075 , .085 , .0975 , .115 , 7 .140625 , .171875, .203125, .234375, .265625, 8 .296875 , .328125, .359375, .390625, .421875, 9 .453125 , .484375, .515625, .546875, .578125, 1 .609375 , .640625, .671875, .703125, .734375, 2 .765625 , .796875, .828125, .859375, .885 , 3 .9 , .915 , .93 , .945 , .96 , 4 .975 , .99 , 1.0 , 1.01 , 1.025 , 5 1.05 , 1.09 , 1.15 , 1.225 , 1.3 , 6 1.4 , 1.5 , 1.625 , 1.75 , 1.875 , 7 23*0.0 / DATA YFREE / -5.2 , -4.4 , -3.6 , -3.0 , -2.4 , 1 -1.95 , -1.6 , -1.35 , -1.15 , -.95 , 2 -.80 , -.65 , -.55 , -.45 , -.39 , 3 -.34 , -.30 , -.27 , -.24 , -.21 , 4 -.18 , -.15 , -.125 , -.1 , -.075 , 5 - .05 , - .03 , - .01 , .01 , .03 , 6 .05 , .075 , .1 , .125 , .15 , 7 .18 , .21 , .24 , .27 , .30 , 8 .34 , .39 , .45 , .55 , .65 , 9 .8 , .95 , 1.15 , 1.35 , 1.60 , 1 1.95 , 2.4 , 3.0 , 3.6 , 4.4 , 2 5.2 , 44*0.0 / DATA YTUN / -2.0 , -1.8 , -1.6 , -1.4 , -1.2 , 1 -1.0 , -.8 , -.65 , -.55 , -.45 , 2 -.39 , -.34 , -.30 , -.27 , -.24 , 3 -.21 , -.18 , -.15 , -.125 , -.1 , 4 - .075 , - .05 , - .03 , - .01 , .01 , 5 .03 , .05 , .075 , .1 , .125 , 6 .15 , .18 , .21 , .24 , .27 , 7 .3 , .34 , .39 , .45 , .55 , 8 .65 , .8 , 1.0 , 1.2 , 1.4 , 9 1.6 , 1.8 , 2.0 , 52*0.0 / DATA XU / 0.000008, 0.000167, 0.000391, 0.000799, 0.001407, 1 0.002153, 0.003331, 0.005336, 0.008648, 0.014583, 2 0.023481, 0.033891, 0.040887, 0.053973, 0.056921, 3 0.058456, 0.059966, 0.061445, 0.062909, 0.065925, 4 0.068785, 0.071482, 0.074007, 0.075322, 0.076603, 5 0.077862, 0.079112, 0.080445, 0.081819, 0.083269, 6 0.084841, 0.086702, 0.088848, 0.091378, 0.094413, 7 0.098308, 0.103104, 0.109010, 0.116244, 0.125452, 8 0.136635, 0.150037, 0.165853, 0.184699, 0.195177, 9 0.206361, 0.218244, 0.230813, 0.244047, 0.257917, 1 0.272371, 0.287410, 0.302990, 0.319057, 0.335555, 2 0.352421, 0.369591, 0.386995, 0.404133, 0.421391, 3 0.438708, 0.456013, 0.473246, 0.490343, 0.507242, 4 0.523881, 0.539536, 0.554867, 0.569823, 0.584351, 5 0.598405, 0.611936, 0.624904, 0.637273, 0.648435, 6 0.659016, 0.668987, 0.678321, 0.687012, 0.695090, 7 0.706936, 0.728406, 0.738649, 0.761390, 0.777010, 8 0.792241, 0.809068, 0.824992, 0.836953, 0.857188, 9 0.875621, 0.898268, 0.913686, 0.927686, 0.939804, 1 0.952002, 0.971789, 0.989100, 0.997860, 1.000000/ DATA YU / 0.000787, 0.003092, 0.004538, 0.006137, 0.007683, 1 0.009056, 0.010675, 0.012803, 0.015607, 0.019624, 2 0.024441, 0.029035, 0.031698, 0.035966, 0.036837, 3 0.037277, 0.037700, 0.038103, 0.038497, 0.039276, 4 0.039986, 0.040625, 0.041195, 0.041483, 0.041756, 5 0.042019, 0.042274, 0.042539, 0.042804, 0.043079, 6 0.043368, 0.043700, 0.044072, 0.044497, 0.044989, 7 0.045595, 0.046312, 0.047154, 0.048132, 0.049301, 8 0.050626, 0.052089, 0.053663, 0.055351, 0.056210, 9 0.057068, 0.057918, 0.058751, 0.059559, 0.060335, 1 0.061068, 0.061751, 0.062381, 0.062947, 0.063445, 2 0.063867, 0.064213, 0.064473, 0.064646, 0.064733, 3 0.064735, 0.064651, 0.064477, 0.064218, 0.063871, 4 0.063438, 0.062945, 0.062376, 0.061731, 0.061014, 5 0.060232, 0.059389, 0.058496, 0.057562, 0.056650, 6 0.055721, 0.054791, 0.053867, 0.052965, 0.052086, 7 0.050722, 0.048045, 0.046680, 0.043441, 0.041053, 8 0.038606, 0.035768, 0.032958, 0.030775, 0.026954, 9 0.023361, 0.018848, 0.015750, 0.012954, 0.010567, 1 0.008213, 0.004559, 0.001620, 0.000293, 0.000000/ DATA XL / 0.000000, 0.000012, 0.000043, 0.000183, 0.000249, 1 0.000348, 0.000455, 0.000680, 0.001011, 0.001481, 2 0.001875, 0.002316, 0.003055, 0.004201, 0.004747, 3 0.005779, 0.007035, 0.008265, 0.009969, 0.012286, 4 0.015346, 0.019276, 0.025335, 0.029379, 0.039095, 5 0.052516, 0.062469, 0.073329, 0.085290, 0.099822, 6 0.118563, 0.140987, 0.167184, 0.202933, 0.228511, 7 0.247895, 0.263995, 0.282047, 0.297045, 0.310147, 8 0.324075, 0.344872, 0.363502, 0.387644, 0.404492, 9 0.426308, 0.450016, 0.475378, 0.521837, 0.549843, 1 0.578612, 0.605305, 0.623479, 0.642152, 0.657543, 2 0.671212, 0.690340, 0.708891, 0.726684, 0.746683, 3 0.768502, 0.784892, 0.801149, 0.819187, 0.838548, 4 0.858817, 0.879431, 0.903723, 0.926504, 0.943652, 5 0.958668, 0.973623, 0.986187, 0.996582, 1.000000, 6 25*0.0/ DATA YL / 0.000000,-0.000700,-0.001385,-0.002868,-0.003330, 1 -0.003880,-0.004379,-0.005199,-0.006133,-0.007183, 2 -0.007933,-0.008676,-0.009776,-0.011204,-0.011815, 3 -0.012861,-0.013983,-0.014962,-0.016175,-0.017636, 4 -0.019336,-0.021258,-0.023836,-0.025373,-0.028634, 5 -0.032423,-0.034840,-0.037182,-0.039456,-0.041862, 6 -0.044483,-0.047017,-0.049298,-0.051443,-0.052406, 7 -0.052859,-0.053062,-0.053117,-0.053027,-0.052849, 8 -0.052562,-0.051951,-0.051218,-0.050013,-0.049004, 9 -0.047495,-0.045601,-0.043288,-0.038336,-0.034916, 1 -0.031104,-0.027333,-0.024661,-0.021854,-0.019517, 2 -0.017429,-0.014527,-0.011771,-0.009228,-0.006537, 3 -0.003868,-0.002086,-0.000524, 0.000950, 0.002227, 4 0.003224, 0.003885, 0.004212, 0.004067, 0.003657, 5 0.003067, 0.002242, 0.001329, 0.000376, 0.000000, 6 25*0.0/ END SUBROUTINE BCEND C SUBROUTINE BCEND MODIFIES THE DIAG AND RHS VECTORS C ON EACH I LINE IN THE APPROPRIATE WAY TO INCLUDE THE C BOUNDARY CONDITIONS AT JBOT AND JTOP. C CALLED BY - SYOR. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK COMMON / COM5/ XDIFF(100),YDIFF(100) COMMON /COM19/ DIAG(100), RHS(100), SUB(100), SUP(100) COMMON /COM23/ CYYC(100), CYYD(100),CYYU(100), IVAL INTEGER BCTYPE COMMON /COM28/ BCTYPE , CIRCFF , FHINV , POR , CIRCTE C I = IVAL C BRANCH TO APPROPRIATE ADDRESS FOR BCTYPE GO TO (10,20,30,40,50,60) , BCTYPE C BCTYPE = 1, FREE AIR 10 CONTINUE C DIRCHLET BOUNDARY CONDITION FOR SUBSONIC FREESTREAM IF (AK .GT. 0.0) RETURN C NEUMAN BOUNDARY CONDITION FOR SUPERSONIC FREESTREAM DFACL = -CYYD(JBOT) * RTK * XDIFF(I) DFACU = -CYYU(JTOP) * RTK * XDIFF(I) RFACL = DFACL * (P(JMIN,I) - P(JMIN,I-1)) RFACU = DFACU * (P(JMAX,I) - P(JMAX,I-1)) GO TO 95 C BCTYPE = 2, SOLID WALL 20 CONTINUE C NEUMAN BOUNDARY CONDITION = 0. C NO MODIFICATION NECESSARY TO DIAG OR RHS RETURN C BCTYPE = 3, FREE JET 30 CONTINUE C DIRCHLET BOUNDARY CONDITION IF (AK .LT. 0.0) GO TO 31 PJMIN = -.75 * CIRCFF PJMAX = -.25 * CIRCFF GO TO 32 31 CONTINUE PJMIN = 0.0 PJMAX = 0.0 32 CONTINUE GO TO 90 C BCTYPE = 4, IDEAL SLOTTED WALL 40 CONTINUE C NEUMAN BOUNDARY CONDITION DFACL =-FHINV * CYYD(JBOT) DFACU =-FHINV * CYYU(JTOP) IF (AK .LT. 0.0) GO TO 41 RFACL = DFACL * ( .75 * CIRCFF + P(JBOT,I)) RFACU = DFACU * ( .25 * CIRCFF + P(JTOP,I)) GO TO 42 41 CONTINUE RFACL = DFACL * P(JBOT,I) RFACU = DFACU * P(JTOP,I) 42 CONTINUE GO TO 95 C BCTYPE = 5, POROUS/PERFORATED WALL 50 CONTINUE IF(POR .GT. 1.5) GO TO 55 C NEUMAN BOUNDARY CONDITION FOR POR .LT. 1.5 DFACL = -CYYD(JBOT) * POR * XDIFF(I) DFACU = -CYYU(JTOP) * POR * XDIFF(I) RFACL = DFACL * (P(JMIN,I) - P(JMIN,I-1)) RFACU = DFACU * (P(JMAX,I) - P(JMAX,I-1)) GO TO 95 55 CONTINUE C DIRCHLET BOUNDARY CONDITION FOR POR .GT. 1.5 IF (I .NE. IUP) RETURN C SET VALUES OF P ON BOUNDARY BY INTEGRATING PX USING C OLD VALUES OF POTENTIAL PJMIN = P(JMIN,IUP) TERM = -.5 / (POR * (Y(JMIN) - Y(JMIN+1))) DO 57 II=IUP,IDOWN P(JMIN,II) = P(JMIN,II-1) - TERM * (X(II)-X(II-1)) * 1 (P(JMIN,II)+P(JMIN,II-1)-P(JMIN+1,II)-P(JMIN+1,II-1)) 57 CONTINUE PJMAX = P(JMAX,IUP) TERM = .5 / (POR * (Y(JMAX) - Y(JMAX-1))) DO 58 II=IUP,IDOWN P(JMAX,II) = P(JMAX,II-1) - TERM*(X(II) - X(II-1)) * 1 (P(JMAX,II)+P(JMAX,II-1)-P(JMAX-1,II)-P(JMAX-1,II-1)) 58 CONTINUE RHS(JBOT) = RHS(JBOT) - (CYYD(JBOT)*(P(JBOT-1,I)-PJMIN)) RHS(JTOP) = RHS(JTOP) - (CYYU(JTOP)*(P(JTOP+1,I)-PJMAX)) RETURN C BCTYPE = 6, GENERAL WALL BOUNDARY CONDITION C DIFFERENCE EQUATIONS FOR THIS BOUNDARY CONDITION C HAVE NOT YET BEEN WORKED OUT. USER MUST INSERT C INFORMATION NEEDED FOR CALCULATION 60 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,1000) WRITE(15,1000) 1000 FORMAT(34H1ABNORMAL STOP IN SUBROUTINE BCEND/ 1 24H BCTYPE=6 IS NOT USEABLE) STOP C DIRCHLET BOUNDARY CONDITIONS 90 CONTINUE RHS(JBOT) = RHS(JBOT) - (CYYD(JBOT)*(PJMIN-P(JBOT-1,I))) RHS(JTOP) = RHS(JTOP) - (CYYU(JTOP)*(PJMAX-P(JTOP+1,I))) RETURN 95 CONTINUE C NEUMAN BOUNDARY CONDITIONS DIAG(JBOT) = DIAG(JBOT) + DFACL DIAG(JTOP) = DIAG(JTOP) + DFACU RHS(JBOT) = RHS(JBOT) - RFACL + CYYD(JBOT)*P(JBOT-1,I) RHS(JTOP) = RHS(JTOP) - RFACU + CYYU(JTOP)*P(JTOP+1,I) RETURN END SUBROUTINE BODY C COMPUTES BODY GEOMETRY INFORMATION FOR BOUNDARY C CONDITIONS AND OUTPUT INFORMATION. Five CHOICES OF C BODY DESCRIPTION ARE AVAILABLE AS FOLLOWS C BCFOIL = 1 NACA 00XX AIRFOIL C BCFOIL = 2 PARABOLIC ARC AIRFOIL C BCFOIL = 3 AIRFOIL ORDINATES READ IN C BCFOIL = 4 JAMESON'S AIRFOIL INPUT FORMAT c BCFOIL = 5 diamond airfoil C BODY ORDINATES AND SLOPES ARE COMPUTED AT THE INPUT C X MESH LOCATIONS AND ARE DIVIDED BY THE THICKNESS C RATIO DELTA. THE BODY VOLUME, CAMBER, AND THICKNESS C ARE ALSO COMPUTED. C ACTUAL BOUNDARY CONDITION IS SET IN SUBROUTINE SETBC C CALLED BY - TSFOIL. C COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK LOGICAL AMESH COMMON / COM4/ XIN(100) , YIN(100), AMESH COMMON / COM6/ FL(100) , FXL(100), FU(100) , FXU(100), 1 CAMBER(100), THICK(100),VOL , XFOIL(100), IFOIL INTEGER BCFOIL COMMON / COM9/ BCFOIL , NL , NU , XL(100) , XU(100) , 1 YL(100) , YU(100) , RIGF,IFLAP,DELFLP,FLPLOC LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT common /com99/ iread COMMON /SPLN/ A(200) , B(200) , DY1 , DY2 , K1 , 1 K2 , XP , YP ,DYP C C SET NUMBER OF POINTS ON AIRFOIL IFOIL = ITE - ILE + 1 C ZERO ALL THICKNESSES AND SLOPES DO 10 I=IMIN,IMAX FU(I) = 0. FL(I) = 0. FXU(I) = 0. FXL(I) = 0. 10 CONTINUE C BRANCH TO APPROPRIATE AIRFOIL SPECIFICATION GO TO (100,200,300,400,250), BCFOIL 100 CONTINUE C BCFOIL = 1 C FORMULA FOR NACA 00XX SHAPE IC = 0 DO 125 I = ILE,ITE IC = IC + 1 Z = XIN(I) XFOIL(IC) = Z RTZ = SQRT(Z) Z2 = Z*Z Z3 = Z*Z2 Z4 = Z*Z3 FU(IC) = 1.4845*RTZ - .63*Z - 1.758*Z2 + 1.4215*Z3 - .5075*Z4 FL(IC) = -FU(IC) FXU(IC) = .74225/RTZ - .63 - 3.516*Z + 4.2645*Z2 - 2.03*Z3 FXL(IC) = -FXU(IC) 125 CONTINUE GO TO 500 200 CONTINUE C BCFOIL = 2 C PARABOLIC ARC AIRFOIL*******BI-CONVEX. IC = 0 DO 225 I = ILE,ITE IC = IC + 1 Z = XIN(I) XFOIL(IC) = Z Z2 = Z*Z FU(IC) = 2.*(Z - Z2) FL(IC) = -FU(IC) FXU(IC) = 2. - 4.*Z FXL(IC) = -FXU(IC) 225 CONTINUE go to 500 c mod to treat diamond airfoil 250 continue c bcfoil = 5 c diamond airfoil ic = 0 do 265 i = ile,ite ic = ic + 1 z = xin(i) xfoil(ic) = z if (z .le. 0.5) fu(ic) = .5*z if (z .gt. 0.5) fu(ic) = (1. - z) fl(ic) = -fu(ic) if (z .le. 0.5) fxu(ic) = 1.0 if (z .gt. 0.5) fxu(ic) = -1.0 fxl(ic) = -fxu(ic) 265 continue GO TO 500 300 CONTINUE C BCFOIL = 3 C BODY ORDINATES READ IN NAMELIST DELINV = 1. IF(PHYS) DELINV = 1./DELTA C COMPUTE ORDINATES AND SLOPES AT X MESH LOCATION ON C AIRFOIL BY CUBIC SPLINE INTERPOLATION. C DERIVATIVE END CONDITIONS ARE SPECIFIED AT X=0 C AND X=1 K1 = 1 K2 = 1 C UPPER SURFACE C CALCULATE DY/DX AT END POINTS BY FINITE DIFFERENCE C FORMULA. DY1 = (YU(2) - YU(1)) / (XU(2) - XU(1)) DY2 = (YU(NU) - YU(NU-1)) / (XU(NU) - XU(NU-1)) C INITIALIZE CUBIC SPLINE INTERPOLATION CALL SPLN1(XU,YU,NU) C CALCULATE ORDINATES AND SLOPES IC = 0 DO 320 I=ILE,ITE IC = IC + 1 XP = XIN(I) XFOIL(IC) = XP CALL SPLN1X(XU,YU,NU) C SET ORDINATE AND SLOPE OF AIRFOIL TO INTERPOLATED C VALUE DIVIDED BY THICKNESS RATIO FU(IC) = YP*DELINV FXU(IC) = DYP*DELINV 320 CONTINUE C LOWER SURFACE C CALCULATE DY/DX AT END POINTS BY FINITE DIFFERENCE C FORMULA. DY1 = (YL(2) - YL(1)) / (XL(2) - XL(1)) DY2 = (YL(NL) - YL(NL-1)) / (XL(NL) - XL(NL-1)) C INITIALIZE CUBIC SPLINE INTERPOLATION CALL SPLN1(XL,YL,NL) C CALCULATE ORDINATES AND SLOPES IC = 0 DO 330 I=ILE,ITE IC = IC + 1 XP = XIN(I) CALL SPLN1X(XL,YL,NL) C SET ORDINATE AND SLOPE OF AIRFOIL TO INTERPOLATED C VALUE DIVIDED BY THICKNESS RATIO FL(IC) = YP* DELINV FXL(IC) = DYP*DELINV 330 CONTINUE GO TO 500 400 CONTINUE C BCFOIL = 4 C AIRFOIL INPUT IN JAMESON'S FORMAT K1 = 1 K2 = 1 DELINV = 1. IF (PHYS) DELINV = 1./DELTA READ(iread,470) READ(iread,480) FSYM,FNU,FNL ISYM = FSYM NU = FNU NL = FNL READ(iread,470) READ(iread,490) (XU(N),YU(N),N=1,NU) DY1 = (YU(2) - YU(1)) / (XU(2) - XU(1)) C Bug: N = NU + 1, and therefore, the derivative is wrong - Andy Ko 4/9/03 C DY2 = (YU(N) - YU(N-1)) / (XU(N) - XU(N-1)) DY2 = (YU(NU) - YU(NU-1)) / (XU(NU) - XU(NU-1)) CALL SPLN1(XU,YU,NU) IC = 0 DO 410 I=ILE,ITE IC = IC + 1 XP = XIN(I) XFOIL(IC) = XP CALL SPLN1X(XU,YU,NU) FU(IC) = YP * DELINV FXU(IC) = DYP * DELINV 410 CONTINUE IF (ISYM.EQ.0) GO TO 440 NL = NU IC = 0 DO 420 I=ILE,ITE IC = IC + 1 FL(IC) = -FU(IC) 420 FXL(IC) = -FXU(IC) DO 430 N=1,NL XL(N) = XU(N) 430 YL(N) = -YU(N) GO TO 500 440 CONTINUE READ(iread,470) READ(iread,490) (XL(N),YL(N),N=1,NL) DY1 = (YL(2) - YL(1)) / (XL(2) - XL(1)) C Bug: N = NL + 1 and therefore derivative is calculated wrong - Andy Ko 4/9/03 C DY2 = (YL(N) - YL(N-1)) / (XL(N) - XL(N-1)) DY2 = (YL(NL) - YL(NL-1)) / (XL(NL) - XL(NL-1)) CALL SPLN1(XL,YL,NL) IC = 0 DO 450 I=ILE,ITE IC = IC + 1 XP = XIN(I) CALL SPLN1X(XL,YL,NL) FL(IC) = YP * DELINV FXL(IC) = DYP * DELINV 450 CONTINUE 470 FORMAT(1X) 480 FORMAT(3F10.0) 490 FORMAT(2F10.0) 500 CONTINUE C EXECUTE CALCULATIONS COMMON TO ALL AIRFOILS C C COMPUTE AIRFOIL VOLUMNE CALL SIMP(VOLU,XFOIL,FU,IFOIL,IERR) CALL SIMP(VOLL,XFOIL,FL,IFOIL,IERR) VOL = VOLU - VOLL C ADD IN FLAP DEFLECTION IF (IFLAP .EQ. 0) GO TO 524 DFLAP = DELFLP/57.29578 SDFLAP = SIN(DFLAP) DO 505 I=1,IFOIL IF (XFOIL(I) .GE. FLPLOC) GO TO 510 505 CONTINUE GO TO 524 510 IFP = I DO 520 I = IFP,IFOIL DELY = (XFOIL(I)-FLPLOC)*SDFLAP*DELINV FU(I) = FU(I) - DELY FL(I) = FL(I) - DELY FXU(I) = FXU(I) - DFLAP*DELINV FXL(I) = FXL(I) - DFLAP*DELINV 520 CONTINUE C COMPUTE CAMBER AND THICKNESS 524 DO 525 I = 1,IFOIL CAMBER(I) = .5*(FU(I) + FL(I)) THICK(I) = .5*(FU(I) - FL(I)) 525 CONTINUE DO 530 I = 1,IFOIL FXU(I) = FXU(I) / SQRT(1.0 + RIGF *(DELTA*FXU(I))**2 ) FXL(I) = FXL(I) / SQRT(1.0 + RIGF *(DELTA*FXL(I))**2 ) 530 CONTINUE CALL PRBODY RETURN END SUBROUTINE CDCOLE C COMPUTES THE DRAG BY MOMENTUM INTEGRAL METHOD. C INTEGRATES AROUND A CONTOUR ENCLOSING THE BODY AND C ALONG ALL SHOCKS INSIDE THE CONTOUR C CALLED BY - PRINT. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK COMMON / COM6/ FL(100) , FXL(100), FU(100) , FXU(100), 1 CAMBER(100), THICK(100),VOL , XFOIL(100), IFOIL COMMON / COM7/ CJUP , CJUP1 , CJLOW , CJLOW1 COMMON /COM13/ CDFACT , CLFACT , CMFACT , CPFACT , CPSTAR LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT COMMON/COM30/ Z(100), ARG(100), REST(204) GAM123 = GAM1*2./3. C SET LOCATIONS OF CONTOUR BOUNDARIES C C UPSTREAM BOUNDARY C IF AK = 0.0 CDCOLE WILL NOT BE CALLED. C AMACH MAY NOT BE = 1.0 . IF(AK .GT. 0.) IU = (ILE + IMIN)*.5 IF(AK .LT. 0.) IU = IUP C TOP AND BOTTOM BOUNDARIES C SUBSONIC FREESTREAM C SET JB,JT TO INCLUDE AS MUCH OF SHOCKS AS POSSIBLE JT = JMAX - 1 JB = JMIN + 1 IF (AK .GT. 0.) GO TO 30 C SUPERSONIC FREESTREAM C SET JB,JT TO INCLUDE ONLY SUBSONIC PART OF C DETACHED BOW WAVE C C FIND BOW SHOCK WAVE ISTOP = ILE - 3 CALL FINDSK(IUP,ISTOP,JUP,IBOW) IF(IBOW .LT. 0) GO TO 325 C SEARCH UP SHOCK TO FIND TIP OF SUBSONIC REGION ISK = IBOW JSTART = JUP + 1 JT = JUP - 1 DO 10 J = JSTART,JMAX JT = JT + 1 ISKOLD = ISK CALL NEWISK(ISKOLD,J, ISK) IF(ISK .LT. 0) GO TO 15 10 CONTINUE 15 CONTINUE C SEARCH DOWN SHOCK TO FIND TIP OF SUBSONIC REGION ISK = IBOW JB = JLOW + 2 DO 20 J = JMIN,JLOW JJ = JLOW - J + JMIN JB = JB - 1 ISKOLD = ISK CALL NEWISK(ISKOLD,JJ,ISK) IF(ISK .LT. 0) GO TO 25 20 CONTINUE 25 CONTINUE C SAVE I LOCATION OF BOW SHOCK WAVE ON LOWER BOUNDARY IBOW = ISKOLD 30 CONTINUE C DOWNSTREAM BOUNDARY ID = (ITE + IMAX) * .5 IF(PX(ITE+1,JUP) .LT. SONVEL) GO TO 40 C TRAILING EDGE IS SUPERSONIC. PLACE DOWNSTREAM C BOUNDARY AHEAD OF TRAILING EDGE TO AVOID TAIL SHOCK I = ITE 35 CONTINUE I = I-1 IF(X(I) .GT. .75) GO TO 35 ID = I 40 CONTINUE C C ALL BOUNDARIES ARE FIXED C C COMPUTE INTEGRALS ALONG BOUNDARIES C C INTEGRAL ON UPSTREAM BOUNDARY CDUP = 0. IF(AK .LT. 0.) GO TO 120 L = 0 DO 110 J = JB,JT L = L+1 Z(L) = Y(J) U = PX(IU,J) V = PY(IU,J) ARG(L) =((AK - GAM123*U)*U*U - V*V)*.5 110 CONTINUE CALL TRAP(Z,ARG,L,SUM) CDUP = 2.*CDFACT*SUM 120 CONTINUE C INTEGRAL ON TOP BOUNDARY L = 0 DO 130 I = IU,ID L = L + 1 Z(L) = X(I) ARG(L) = -PX(I,JT)*PY(I,JT) 130 CONTINUE CALL TRAP(Z,ARG,L,SUM) CDTOP = 2.*CDFACT*SUM C INTEGRAL ON BOTTOM BOUNDARY L = 0 DO 140 I = IU,ID L = L + 1 ARG(L) = PX(I,JB)*PY(I,JB) 140 CONTINUE CALL TRAP(Z,ARG,L,SUM) CDBOT = 2.*CDFACT*SUM C INTEGRAL ON DOWNSTREAM BOUNDARY L = 0 DO 150 J = JB,JT L = L + 1 Z(L) = Y(J) U = PX(ID,J) C IF FLOW SUPERSONIC, USE BACKWARD DIFFERENCE FORMULA IF(U .GT. SONVEL) U = PX(ID-1,J) V = PY(ID,J) ARG(L) = ((GAM123*U - AK)*U*U + V*V)*.5 150 CONTINUE CALL TRAP(Z,ARG,L,SUM) CDDOWN = 2.*CDFACT*SUM C INTEGRAL ON BODY BOUNDARY CDBODY = 0. IF(ID .GT. ITE) GO TO 200 ILIM = ITE + 1 L = 0 DO 160 I=ID,ILIM IB = I-ILE + 1 L = L + 1 Z(L) = X(I) UU = CJUP*PX(I,JUP) - CJUP1*PX(I,JUP+1) UL = CJLOW*PX(I,JLOW) - CJLOW1*PX(I,JLOW-1) ARG(L) = -UU*FXU(IB) + UL*FXL(IB) 160 CONTINUE CALL TRAP(Z,ARG,L,SUM) CDBODY = 2.*CDFACT*SUM 200 CONTINUE C C INTEGRATION ALONG SHOCK WAVES CDWAVE = 0. LPRT1 = 0 LPRT2 = 0 NSHOCK = 0 IF(AK .GT. 0.) GO TO 220 C INTEGRATE ALONG DETACHED BOW WAVE NSHOCK = NSHOCK + 1 LPRT1 = 1 LPRT2 = 1 L = 0 ISK = IBOW DO 210 J = JB,JT L = L + 1 ISKOLD = ISK CALL NEWISK(ISKOLD,J,ISK) Z(L) = Y(J) ARG(L) = (PX(ISK+1,J) - PX(ISK-2,J))**3 210 CONTINUE CALL TRAP(Z,ARG,L,SUM) CDSK = -GAM1/6.*CDFACT*SUM CDWAVE = CDWAVE + CDSK CALL PRTSK(Z,ARG,L,NSHOCK, CDSK,LPRT1) 220 CONTINUE C INTEGRATE ALONG SHOCKS ABOVE AIRFOIL ISTART = ILE 225 CONTINUE CALL FINDSK(ISTART,ITE,JUP,ISK) IF(ISK .LT. 0) GO TO 250 C SHOCK WAVE FOUND ISTART = ISK + 1 NSHOCK = NSHOCK + 1 LPRT1 = 0 L = 1 Z(L) = 0. ARG(L) = (CJUP*(PX(ISK+1,JUP) - PX(ISK-2,JUP)) + - CJUP1*(PX(ISK+1,JUP+1) - PX(ISK-2,JUP+1)))**3 DO 230 J = JUP , JT L = L+1 Z(L) = Y(J) ARG(L) = (PX(ISK+1,J) - PX(ISK-2,J))**3 ISKOLD = ISK JSK = J + 1 CALL NEWISK(ISKOLD,JSK,ISK) IF(ISK .LT. 0) GO TO 240 IF(ISK .GT. ID) GO TO 235 230 CONTINUE 235 CONTINUE LPRT1 = 1 240 CONTINUE CALL TRAP(Z,ARG,L,SUM) CDSK = -GAM1/6.*CDFACT*SUM CDWAVE = CDWAVE + CDSK CALL PRTSK(Z,ARG,L,NSHOCK,CDSK,LPRT1) IF(LPRT1 .EQ. 1) LPRT2 = 1 C RETURN TO FIND NEXT SHOCK GO TO 225 C INTEGRATE ALONG SHOCKS BELOW AIRFOIL 250 CONTINUE ISTART = ILE 260 CONTINUE CALL FINDSK(ISTART,ITE,JLOW,ISK) IF(ISK .LT. 0) GO TO 300 C SHOCK WAVE FOUND ISTART = ISK + 1 NSHOCK = NSHOCK + 1 LPRT1 = 0 L = 1 Z(L) = 0. ARG(L) = (CJLOW*(PX(ISK+1,JLOW) - PX(ISK-2,JLOW)) + - CJLOW1*(PX(ISK+1,JLOW-1) - PX(ISK-2,JLOW-1)))**3 DO 270 JJ = JB,JLOW J = JLOW + JB - JJ L = L+1 Z(L) = Y(J) ARG(L) = (PX(ISK+1,J) - PX(ISK-2,J))**3 ISKOLD = ISK JSK = J - 1 CALL NEWISK(ISKOLD,JSK,ISK) IF(ISK .LT. 0) GO TO 280 IF(ISK .GT. ID) GO TO 275 270 CONTINUE 275 CONTINUE LPRT1 = 1 280 CONTINUE CALL TRAP(Z,ARG,L,SUM) CDSK = -GAM1/6.*CDFACT*(-SUM) CDWAVE = CDWAVE + CDSK CALL PRTSK(Z,ARG,L,NSHOCK,CDSK,LPRT1) IF(LPRT1 .EQ. 1) LPRT2 = 1 C RETURN TO FIND NEXT SHOCK GO TO 260 300 CONTINUE C C INTEGRATION ALONG SHOCKS IS COMPLETE C C PRINTOUT CD INFORMATION XU = X(IU) XD = X(ID) YT = Y(JT)*YFACT YB = Y(JB)*YFACT CDC = CDUP + CDTOP + CDBOT + CDDOWN + CDBODY CD = CDC + CDWAVE C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,1001) WRITE(15,1001) C WRITE(6,1002) XU,CDUP, XD,CDDOWN, YT,CDTOP, YB,CDBOT C IF(XD .LT. 1.) WRITE(6,1003) XD,CDBODY C WRITE(6,1004) CDC C WRITE(6,1005) NSHOCK,CDWAVE C IF(NSHOCK .GT. 0 .AND. LPRT2 .EQ. 0) WRITE(6,1007) C IF(NSHOCK .GT. 0 .AND. LPRT2 .EQ. 1) WRITE(6,1008) C WRITE(6,1006) CD WRITE(15,1002) XU,CDUP, XD,CDDOWN, YT,CDTOP, YB,CDBOT IF(XD .LT. 1.) WRITE(15,1003) XD,CDBODY WRITE(15,1004) CDC WRITE(15,1005) NSHOCK,CDWAVE IF(NSHOCK .GT. 0 .AND. LPRT2 .EQ. 0) WRITE(15,1007) IF(NSHOCK .GT. 0 .AND. LPRT2 .EQ. 1) WRITE(15,1008) WRITE(15,1006) CD RETURN 325 CONTINUE C SHOCK IS TOO CLOSE TO BODY TO DO CONTOUR INTEGRAL. C WRITE MESSAGE AND RETURN ULE = PX(ILE,JUP) C Changed to write to a file - Andy Ko 4/3/03 C IF(ULE .GT. SONVEL) WRITE(6,1011) C IF(ULE .LE. SONVEL) WRITE(6,1012) IF(ULE .GT. SONVEL) WRITE(15,1011) IF(ULE .LE. SONVEL) WRITE(15,1012) CD = DRAG(CDFACT) C WRITE(6,1013) CD WRITE(15,1013) CD RETURN 1001 FORMAT(60H1CALCULATION OF DRAG COEFFICIENT BY MOMENTUM INTEGRAL ME 1THOD) 1002 FORMAT(27H0BOUNDARIES OF CONTOUR USED,15X,18HCONTRIBUTION TO CD/ & 16H UPSTREAM X =,F12.6,15X,8HCDUP =,F12.6/ & 16H DOWNSTREAM X =,F12.6,15X,8HCDDOWN =,F12.6/ & 16H TOP Y =,F12.6,15X,8HCDTOP =,F12.6/ & 16H BOTTOM Y =,F12.6,15X,8HCDBOT =,F12.6) 1003 FORMAT(16H BODY AFT OF X =,F12.6,15X,8HCDBODY =,F12.6) 1004 FORMAT(15X,36HTOTAL CONTRIBUTIONS AROUND CONTOUR =,F12.6) 1005 FORMAT(10H0THERE ARE,I3,38H SHOCKS INSIDE CONTOUR. TOTAL CDWAVE =, & F12.6) 1006 FORMAT(51H0DRAG CALCULATED FROM MOMENTUM INTEGRAL CD =, & F12.6) 1007 FORMAT(43H0NOTE - ALL SHOCKS CONTAINED WITHIN CONTOUR/ & 30H CDWAVE EQUALS TOTAL WAVE DRAG) 1008 FORMAT(52H0NOTE - ONE OR MORE SHOCKS EXTEND OUTSIDE OF CONTOUR/ & 38H CDWAVE DOES NOT EQUAL TOTAL WAVE DRAG) 1011 FORMAT(31H1SHOCK WAVE IS ATTACHED TO BODY/ & 33H MOMENTUM INTEGRAL CANNOT BE DONE/ & 45H DRAG OBTAINED FROM SURFACE PRESSURE INTEGRAL/) 1012 FORMAT(41H1DETACHED SHOCK WAVE IS TOO CLOSE TO BODY/ & 33H MOMENTUM INTEGRAL CANNOT BE DONE/ & 45H DRAG OBTAINED FROM SURFACE PRESSURE INTEGRAL/) 1013 FORMAT(4H0CD=,F12.6) END SUBROUTINE CKMESH C C CHECK X MESH AND ADJUST TO CONTAIN ODD NO. OF C POINTS BEFORE TAIL AND ODD NO. AFTER TAIL. C ITE IS INCLUDED IN BOTH COUNTS. C CKMESH IS CALLED ONCE BY READIN. C COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 LOGICAL ABORT1 COMMON / COM3/ IREF , ABORT1 , ICUT , KSTEP LOGICAL AMESH COMMON / COM4/ XIN(100) , YIN(100), AMESH C IF (ICUT .GT. 0) GO TO 4 IREF = -1 RETURN 4 CONTINUE C C TEST TO BE SURE THAT ADJUSTING THE NO. OF C POINTS WONT MAKE IMAX OR JMAX LARGER THAN 100. IF (IMAX .LE. 98 .AND. JMAX .LE. 98) GO TO 8 C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,100) WRITE (15,100) IREF = -1 RETURN 8 CONTINUE IF (MOD(ITE-IMIN+1, 2) .NE. 0) GO TO 20 C ADD EXTRA MESH POINT AHEAD OF AIRFOIL. LP = IMAX + IMIN + 1 DO 10 I=IMIN,IMAX L = LP - I XIN(L) = XIN(L-1) 10 CONTINUE IMAX = IMAX + 1 XIN(IMIN) = 2. * XIN(IMIN+1) - XIN(IMIN+2) CALL ISLIT ( XIN ) 20 CONTINUE C ADD EXTRA MESH POINT AFTER AIRFOIL. IF (MOD(IMAX-ITE+1, 2) .NE. 0) GO TO 30 IMAX = IMAX + 1 XIN(IMAX) = 2. * XIN(IMAX-1) - XIN(IMAX-2) 30 CONTINUE C CHECK Y MESH AND ADJUST TO CONTAIN EVEN NO. OF C POINTS ABOVE AND BELOW SLIT. IF (MOD(JLOW-JMIN ,2) .NE. 0) GO TO 50 C ADD EXTRA MESH POINT BELOW SLIT. LP = JMAX + JMIN + 1 DO 40 J=JMIN,JMAX L = LP - J YIN(L) = YIN(L-1) 40 CONTINUE JMAX = JMAX + 1 YIN(JMIN) = 2. * YIN(JMIN+1) - YIN(JMIN+2) CALL JSLIT ( YIN ) 50 CONTINUE C ADD EXTRA MESH POINT ABOVE SLIT. IF (MOD(JMAX-JUP ,2) .NE. 0) GO TO 60 JMAX = JMAX + 1 YIN(JMAX) = 2.0 * YIN(JMAX-1) - YIN(JMAX-2) 60 CONTINUE RETURN C 100 FORMAT(96H0 THE MESH CANNOT BE ADJUSTED FOR CUTOUT, BECAUSE IMAX O 1R JMAX IS TOO CLOSE TO THE LIMIT OF 100./ 25X,19HIREF WAS SET TO 0 ) C END SUBROUTINE CPPLOT (X, Y, Z, W, NP) C C SUBROUTINE CPPLOT PRODUCES A PRINTER PLOT C OF CRITICAL PRESSURE VS X . C CALLED BY - FIXPLT. LOGICAL AMESH COMMON / COM4/ XIN(100) , YIN(100), AMESH C DIMENSION X(101) , Y(101) , Z(101) , W(101) , M(120), ISYM(8) DIMENSION A(3) , IC(3) DATA IC/ 1, 1024, 1048576/ DATA ISYM/1H ,1HU,1HL,1HB,1H-,1HU,1HL,1HB/ C NC IS THE NUMBER OF COLUMNS. C NR IS THE NUMBER OF ROWS. DATA NC /120/ , NR /50/ C INITIALIZE RANGES IF ( AMESH ) GO TO 3 NPL = 1 NPR = NP - 1 NL5 = 2 GO TO 4 3 CONTINUE NPL = 2 NPR = NP - 2 NL5 = 3 4 CONTINUE HL = X(NPL) HR = X(NPL) VB=AMIN1(Y(1),Z(1),W(1)) VT=AMAX1(Y(1),Z(1),W(1)) C DETERMINE RANGES DO 5 I = NL5,NPR HL=AMIN1(HL,X(I)) HR=AMAX1(HR,X(I)) VB=AMIN1(VB,Y(I),Z(I),W(I)) VT=AMAX1(VT,Y(I),Z(I),W(I)) 5 CONTINUE C C SKIP TO NEW PAGE AND WRITE PLOT HEADING. C C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,900) WRITE(15,900) VDEL=(VT-VB)/FLOAT(NR) HDEL=(HR-HL)/FLOAT(NC) HDELM = 1.0 / HDEL VL=VT C C DO 100 IROW = 1,NR VH=VL VL=FLOAT(NR-IROW)*VDEL+VB DO 15 I=1,NC 15 M(I)=0 C DO 30 I =NPL,NPR J = MAX0(1,MIN0(NC,1+INT((X(I)-HL) * HDELM))) A(1) = Y(I) A(2) = Z(I) A(3) = W(I) DO 20 K=1,3 IF (A(K) .GT. VH) GO TO 20 IF (A(K) .GT. VL .OR. (A(K) .LE. VB .AND. IROW .EQ. NR)) 1 M(J) = M(J) + IC(K) 20 CONTINUE 30 CONTINUE C DO 90 I=1,NC J = 1 IF (M(I) .LT. IC(3)) GO TO 70 J = J + 4 M(I) = MOD(M(I),IC(3)) 70 CONTINUE IF (M(I) .LT. IC(2)) GO TO 75 J = J + 2 M(I) = MOD(M(I),IC(2)) 75 CONTINUE IF (M(I) .LE. 0) GO TO 80 J = J + 1 80 CONTINUE M(I) = ISYM(J) 90 CONTINUE C C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,901) (M(I),I=1,NC) WRITE(15,901) (M(I),I=1,NC) 100 CONTINUE C RETURN C 900 FORMAT(1H1,34X,50HPRINTER PLOT OF CP ON BODY AND DIVIDING STREAMLI 1NE//9X,17H U FOR CP(UPPER),5X,17H L FOR CP(LOWER), 2 5X,27H B FOR CP(UPPER)=CP(LOWER),5X,18H --- FOR CP SONIC) 901 FORMAT(1X,120A1) C END SUBROUTINE CUTOUT C C SUBROUTINE TO REDUCE THE NUMBER OF MESH POINTS C FOR THE FIRST CUT AT SOLUTION. THE X-MESH AND C Y-MESH WILL BE HALVED, AND IF POSSIBLE BE C HALVED AGAIN. C IREF =-1, CUTOUT SHOULDNT HAVE BEEN CALLED C IREF = 0, IF NOT HALVED AT ALL. C IREF = 1, IF HALVED ONCE. C IREF = 2, IF HALVED TWICE. C CALLED BY - TSFOIL. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 LOGICAL ABORT1 COMMON / COM3/ IREF , ABORT1 , ICUT , KSTEP LOGICAL AMESH COMMON / COM4/ XIN(100) , YIN(100), AMESH COMMON /COM20/ XMID(100), YMID(100) C IF(IREF .NE. -1) GO TO 5 C MESH CANNOT BE REFINED, LOAD XIN,YIN INTO X,Y DO 1 I = IMIN,IMAX X(I) = XIN(I) 1 CONTINUE DO 2 J = JMIN,JMAX Y(J) = YIN(J) 2 CONTINUE IREF = 0 RETURN 5 CONTINUE K = IMIN - 1 DO 10 I=IMIN,IMAX,2 K = K + 1 XMID(K) = XIN(I) 10 CONTINUE IMAX = (IMAX-IMIN) * .5 + IMIN CALL ISLIT ( XMID ) K = JMIN - 1 JE = JLOW - 1 DO 15 J=JMIN,JE,2 K = K + 1 YMID(K) = YIN(J) 15 CONTINUE JST = JUP + 1 DO 20 J=JST,JMAX,2 K = K + 1 YMID(K) = YIN(J) 20 CONTINUE JMAX = (JMAX-JMIN) * .5 + JMIN CALL JSLIT ( YMID ) IREF = 1 C FIRST HALVING COMPLETE. CHECK IF NO. OF POINTS C IS ODD. IF (ICUT .EQ. 1) GO TO 30 IF (MOD(ITE-IMIN+1, 2) .EQ. 0) GO TO 30 IF (MOD(IMAX-ITE+1,2) .EQ. 0) GO TO 30 IF (MOD(JLOW-JMIN ,2) .EQ. 0) GO TO 30 IF (MOD(JMAX-JUP ,2) .NE. 0) GO TO 60 30 CONTINUE C ONLY ONE MESH REFINEMENT POSSIBLE. DO 40 I=IMIN,IMAX X(I) = XMID(I) 40 CONTINUE DO 50 J=JMIN,JMAX Y(J) = YMID(J) 50 CONTINUE RETURN C ALL POINTS ARE ODD SO CUT AGAIN. 60 CONTINUE K = IMIN - 1 DO 70 I=IMIN,IMAX,2 K = K + 1 X(K) = XMID(I) 70 CONTINUE IMAX = (IMAX-IMIN) * .5 + IMIN CALL ISLIT ( X ) K = JMIN - 1 JE = JLOW - 1 DO 75 J=JMIN,JE,2 K = K + 1 Y(K) = YMID(J) 75 CONTINUE JST = JUP + 1 DO 80 J=JST,JMAX,2 K = K + 1 Y(K) = YMID(J) 80 CONTINUE JMAX = (JMAX-JMIN) * .5 + JMIN CALL JSLIT ( Y ) IREF = 2 RETURN C END SUBROUTINE DIFCOE C C COMPUTE DIFFERENCE COEFFICIENTS IN FIELD C CALLED BY - TSFOIL. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK COMMON / COM5/ XDIFF(100),YDIFF(100) COMMON / COM7/ CJUP , CJUP1 , CJLOW , CJLOW1 COMMON /COM17/ CYYBLC , CYYBLD , CYYBLU , CYYBUC , CYYBUD , 1 CYYBUU ,FXLBC(100),FXUBC(100), ITEMP1, ITEMP2 COMMON /COM22/ CXC(100) , CXL(100), CXR(100), CXXC(100),CXXL(100), 1 CXXR(100), C1(100) COMMON /COM23/ CYYC(100), CYYD(100),CYYU(100), IVAL C C COEFFICIENTS FOR (P)X AND (P)XX AT IMIN CXXL(IMIN) = 0.0 CXXR(IMIN) = 0.0 CXXC(IMIN) = 0.0 CXL (IMIN) = 0.0 CXR (IMIN) = 0.0 CXC (IMIN) = 0.0 C COEFFICIENTS FOR (P)X AND (P)XX FROM I=IMIN+1 TO I=IMAX-1 C2 = GAM1*.5 ISTART = IMIN+1 IEND=IMAX-1 DO 1 I=ISTART,IEND DXL=X(I)-X(I-1) DXR=X(I+1)-X(I) DXC = .5 * (X(I+1) - X(I-1)) C FOR VC C1(I)= AK/DXC C FOR (P)X CXL(I) = -C2 / (DXL*DXC) CXR(I) = C2 / (DXR*DXC) CXC(I) = -CXL(I) - CXR(I) C FOR (P)XX CXXL(I) = 1.0 / DXL CXXR(I) = 1.0 / DXR CXXC(I) = CXXL(I) + CXXR(I) 1 CONTINUE C COEFFICIENTS FOR (P)X AND (P)XX AT IMAX DX = X(IMAX) - X(IMAX-1) Q = 1.0 / (DX*DX) C1(IMAX) = AK / DX CXL(IMAX) = -C2 * Q CXR(IMAX) = C2 * Q CXC(IMAX) = 0.0 CXXL(IMAX) = 1.0 / DX CXXR(IMAX) = 1.0 / DX CXXC(IMAX) = CXXL(IMAX) + CXXR(IMAX) C COEFFICIENTS FOR (P)YY AT JMIN. C DYU = Y(JMIN+1) - Y(JMIN) CYYD(JMIN) = 2.0 / DYU CYYU(JMIN) = 2./(DYU*DYU) CYYC(JMIN) = CYYU(JMIN) C COEFICIENTS FOR (P)YY FROM J=JMIN+1 TO J=JMAX-1 JSTART = JMIN + 1 JEND=JMAX-1 DO 2 J=JSTART,JEND DYD=Y(J)-Y(J-1) DYU=Y(J+1)-Y(J) DYC=Y(J+1)-Y(J-1) CYYD(J)=2./(DYD*DYC) CYYU(J)=2./(DYU*DYC) CYYC(J) = CYYD(J) + CYYU(J) 2 CONTINUE C COEFFICIENTS FOR (P)YY AT JMAX. DYD = Y(JMAX) - Y(JMAX-1) CYYD(JMAX) = 2./(DYD*DYD) CYYU(JMAX) = 2./DYD CYYC(JMAX) = CYYD(JMAX) C C COEFFICIENTS FOR VELOCITY FORMULAS C ISTART = IMIN + 1 DO 3 I = ISTART,IMAX XDIFF(I) = 1./(X(I) - X(I-1)) 3 CONTINUE JSTART = JMIN +1 DO 5 J = JSTART , JMAX YDIFF(J) = 1./(Y(J) - Y(J-1)) 5 CONTINUE C C COEFFICIENTS FOR EXTRAPOLATION FORMULAS FOR AIRFOIL SURFACE PROPERTIES C CJLOW = -Y(JLOW-1) / (Y(JLOW) - Y(JLOW-1)) CJLOW1 = -Y(JLOW) / (Y(JLOW) - Y(JLOW-1)) CJUP = Y(JUP+1) / (Y(JUP+1) - Y(JUP)) CJUP1 = Y(JUP) / (Y(JUP+1) - Y(JUP)) C COMPUTE SPECIAL DIFFERENCE COEFFICIENTS FOR PYY C TO USE FOR AIRFOIL BOUNDARY CONDITION C DIFFERENCE COEFFICIENTS FOR UPPER SURFACE CYYBUD = -2./(Y(JUP+1) + Y(JUP)) CYYBUC = -CYYBUD / (Y(JUP+1) - Y(JUP)) CYYBUU = CYYBUC C DIFFERENCE COEFFICIENTS FOR LOWER SURFACE CYYBLU=-2.0/(Y(JLOW) + Y(JLOW-1)) CYYBLC = CYYBLU / (Y(JLOW) - Y(JLOW-1)) CYYBLD = CYYBLC RETURN END C DECK DRAG FUNCTION DRAG(CDFACT) C COMPUTES PRESSURE DRAG BY INTEGRATING U*V AROUND C AIRFOIL USING TRAPEZODIAL RULE. C CALLED BY - CDCOLE. C CALL BLANK COMMON P(102,101),X(100) , Y(100) CALL COM1 COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 CALL COM6 COMMON / COM6/ FL(100) , FXL(100), FU(100) , FXU(100), 1 CAMBER(100), THICK(100),VOL , XFOIL(100), IFOIL CALL COM7 COMMON / COM7/ CJUP , CJUP1 , CJLOW , CJLOW1 COMMON /COM30/ XI(100) , ARG(100) , REST(204) C K = 1 ARG(1) = 0. XI(1) = X(ILE-1) DO 10 I = ILE,ITE K = K+1 PXUP = CJUP*PX(I,JUP) - CJUP1*PX(I,JUP+1) PXLOW = CJLOW*PX(I,JLOW) - CJLOW1*PX(I,JLOW-1) ARG(K) = FXU(K-1)*PXUP - FXL(K-1)*PXLOW XI(K) = X(I) 10 CONTINUE K = K + 1 ARG(K) = 0. XI(K) = X(ITE+1) CALL TRAP(XI,ARG,K,SUM) DRAG = - SUM*CDFACT*2. RETURN END SUBROUTINE DROOTS C COMPUTE CONSTANTS ALPHA0,APLHA1,ALPHA2,OMEGA0,OMEGA1 C OMEGA2 USED IN FORMULA FOR DOUBLET IN SLOTTED WIND C TUNNEL WITH SUBSONIC FREESTREAM C CALLED BY - FARFLD. C COMMON /COM12/ F , H , HALFPI , PI , RTKPOR , 1 TWOPI REAL JET COMMON /COM16/ ALPHA0 , ALPHA1 , ALPHA2 , XSING , 1 OMEGA0 , OMEGA1 , OMEGA2 , JET C ERROR = .00001 C COMPUTE ALPHA0 ALPHA0 = 0. DO 10 I = 1,100 TEMP = ALPHA0 Q = F*TEMP - RTKPOR ALPHA0 = HALFPI - ATAN(Q) DALPHA = ABS(ALPHA0 - TEMP) IF(DALPHA .LT. ERROR) GO TO 15 10 CONTINUE N = 0 GO TO 9999 15 CONTINUE C COMPUTE ALPHA1 ALPHA1 = 0. DO 20 I = 1,100 TEMP = ALPHA1 Q = F*(TEMP - PI) - RTKPOR ALPHA1 = HALFPI - ATAN(Q) DALPHA = ABS(ALPHA1 - TEMP) IF(DALPHA .LT. ERROR) GO TO 25 20 CONTINUE N = 1 GO TO 9999 25 CONTINUE C COMPUTE ALPHA2 ALPHA2 = 0. DO 30 I=1,100 TEMP = ALPHA2 Q = F*(TEMP - TWOPI) - RTKPOR ALPHA2 = HALFPI - ATAN(Q) DALPHA = ABS(ALPHA2 - TEMP) IF(DALPHA .LT. ERROR) GO TO 35 30 CONTINUE N = 2 GO TO 9999 35 CONTINUE C COMPUTE OMEGA0,OMEGA1,OMEGA2 TEMP = 1.0 / TAN(ALPHA0) OMEGA0 = 1./(1. + F/(1.+ TEMP*TEMP)) TEMP = 1.0 / TAN(ALPHA1) OMEGA1 = 1./(1. + F/(1. + TEMP*TEMP)) TEMP = 1.0 / TAN(ALPHA2) OMEGA2 = 1./(1. + F/(1. + TEMP*TEMP)) RETURN C ABNORMAL STOP IF ITERATION FOR ALPHAS NOT CONVERGED 9999 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,1000) N WRITE(15,1000) N 1000 FORMAT(35H1ABNORMAL STOP IN SUBROUTINE DROOTS/ 1 38H0NONCONVERGENCE OF ITERATION FOR ALPHA, I1) STOP END SUBROUTINE ECHINP C PRINTS INPUT CARDS USED FOR RUN. C CALLED BY - TSFOIL. DIMENSION CRD(20) C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,110) WRITE (15,110) 10 CONTINUE READ (5,120,END=30) CRD C YES, NO 20 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,100) CRD WRITE (15,100) CRD GO TO 10 30 CONTINUE REWIND 5 RETURN C 100 FORMAT(1H ,20A4) 110 FORMAT(1H1) 120 FORMAT(20A4) C END FUNCTION EMACH1(U) C FUNCTION EMACH1 COMPUTES LOCAL SIMILARTY PARAMETER C OR LOCAL MACH NUMBER C CALLED BY - MACHMP, PRINT1, PRTFLD. C COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT C C COMPUTE SIMILARITY PARAMETER BASED ON LOCAL VELOCITY AK1 = AK - GAM1*U IF(PHYS) GO TO 5 C RETURN VALUE OF LOCAL SIMILARITY PARAMETER EMACH1 = AK1 RETURN 5 CONTINUE C COMPUTE VALUE OF LOCAL MACH NUMBER AND RETURN C COLE SCALING ARG = DELRT2*AK1 C SPREITER SCALING IF(SIMDEF .EQ. 2) ARG = ARG*EMROOT*EMROOT C KRUPP SCALING IF(SIMDEF .EQ. 3) ARG = ARG*EMACH ARG = 1. - ARG EMACH1 = 0. IF(ARG .GT. 0.)EMACH1 = SQRT(ARG) RETURN END SUBROUTINE EXTRAP(XP,YP,PNEW) C COMPUTE P AT X, YP USING FAR FIELD SOLUTION C FOR SUBSONIC FLOW C CALLED BY - GUESSP. C COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK COMMON /COM12/ F , H , HALFPI , PI , RTKPOR , 1 TWOPI COMMON /COM15/ B , BETA0 , BETA1 , BETA2 , PSI0 , 1 PSI1 , PSI2 REAL JET COMMON /COM16/ ALPHA0 , ALPHA1 , ALPHA2 , XSING , 1 OMEGA0 , OMEGA1 , OMEGA2 , JET INTEGER BCTYPE COMMON /COM28/ BCTYPE , CIRCFF , FHINV , POR , CIRCTE C IF(BCTYPE .NE. 1) GO TO 100 C FREE AIR BOUNDARY CONDITION IF(ABS(YP) .LT. 1.E-6) YP = -1.E-6 XI = XP - XSING ETA = YP*RTK PNEW = - CIRCFF/TWOPI*(ATAN2(ETA,XI) + PI - SIGN(PI,ETA)) 1 + DUB /TWOPI/RTK*(XI/(XI*XI + ETA*ETA)) RETURN 100 CONTINUE C TUNNEL WALL BOUNDARY CONDITION ETA = YP/H XI = (XP - XSING)/(H*RTK) IF(XI .LT. 0.) GO TO 200 C XP IS DOWNSTREAM OF AIRFOIL TERM = ETA IF(BCTYPE .NE. 3) TERM = SIN(ETA*BETA0)/BETA0 PNEW = - .5*CIRCFF*(1. - SIGN(1.,ETA) 1 + (1.-JET)*PSI0*TERM*EXP(-BETA0*XI)) 2 + DUB *.5/(AK*H)*(B + OMEGA0*COS(ETA*ALPHA0) 3 *EXP(-ALPHA0*XI)) RETURN 200 CONTINUE C XP IS UPSTREAM OF AIRFOIL TERM = 0. IF(JET .NE. 0.) TERM = JET*ETA/(1.+F) ARG1 = PI - ALPHA1 ARG2 = PI - BETA2 PNEW = -.5*CIRCFF*(1.-TERM-PSI2*SIN(ETA*ARG2)/ARG2*EXP(ARG2*XI)) 1 -.5*DUB /(AK*H)*((1.-B)*OMEGA1*COS(ETA*ARG1)*EXP(XI*ARG1)) RETURN END SUBROUTINE FARFLD C SUBROUTINE COMPUTES BOUNDARY DATA FOR OUTER C BOUNDARIES. C CALLED BY - TSFOIL. C COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK LOGICAL AMESH COMMON / COM4/ XIN(100) , YIN(100), AMESH COMMON /COM12/ F , H , HALFPI , PI , RTKPOR , 1 TWOPI COMMON /COM15/ B , BETA0 , BETA1 , BETA2 , PSI0 , 1 PSI1 , PSI2 REAL JET COMMON /COM16/ ALPHA0 , ALPHA1 , ALPHA2 , XSING , 1 OMEGA0 , OMEGA1 , OMEGA2 , JET COMMON /COM24/ DTOP(100), DBOT(100),DUP(100), DDOWN(100), 1 VTOP(100), VBOT(100),VUP(100), VDOWN(100) INTEGER BCTYPE COMMON /COM28/ BCTYPE , CIRCFF , FHINV , POR , CIRCTE C C TEST FOR SUPERSONIC OR SUBSONIC FREESTREAM IF(AK .GT. 0.) GO TO 99 IF (F .NE. 0.0 .AND. H .NE. 0.0) GO TO 10 FHINV = 1.0 GO TO 99 10 CONTINUE FHINV = 1.0 / (F * H) C SUPERSONIC FREESTREAM C UPSTREAM BOUNDARY CONDITIONS CORRESPOND TO UNIFORM C UNDISTURBED FLOW C DOWNSTREAM BOUNDARY REQUIRED TO BE SUPERSONIC C TOP AND BOTTOM BOUNDARIES USE SIMPLE WAVE SOLUTION. RETURN C 99 CONTINUE C C SUBSONIC FREESTREAM C C FUNCTIONAL FORM OF THE POTENTIAL ON OUTER BOUNDARIES C IS PRESCRIBED. EQUATIONS REPRESENT ASYMPTOTIC FORM C FOR DOUBLET AND VORTEX IN FREE AIR AND WIND TUNNEL C ENVIRONMENT. DOUBLET AND VORTEX ARE LOCATED AT C X=XSING, Y=0. C ACTUAL BOUNDARY VALUES ARE SET IN SUBROUTINES RECIRC C AND REDUB WHERE THE FUNCTIONAL FORMS ARE MULTIPLIED C BY THE VORTEX AND DOUBLET STRENGTHS. THE BOUNDARY C CONDITIONS ARE CALCULATED HERIN FOR THE INPUT X AND C Y MESH AND VALUES ARE DELETED FOR THE COARSE MESH IN C SUBROUTINE SETBC. C C SET LOCATION OF SINGLULAR VORTEX AND DOUBLET XSING = .5 C SET DEFAULT VALUES FOR TUNNEL WALL PARAMETERS B = 0. OMEGA0 = 1. OMEGA1 = 1. OMEGA2 = 1. JET = 0. PSI0 = 1. PSI1 = 1. PSI2 = 1. C BRANCH TO APPROPRIATE FORMULAS DEPENDING ON BCTYPE GO TO (100,200,300,400,500,600),BCTYPE 100 CONTINUE C BCTYPE = 1 C FREE AIR BOUNDARY CONDITION C SET BOUNDARY ORDINATES YT = YIN(JMAX)*RTK YB = YIN(JMIN)*RTK XU = XIN(IMIN) - XSING XD = XIN(IMAX) - XSING YT2 = YT*YT YB2 = YB*YB XU2 = XU*XU XD2 = XD*XD COEF1 = 1./TWOPI COEF2 = 1.0/(TWOPI*RTK) C COMPUTE DOUBLET AND VORTEX TERMS ON TOP AND BOTTOM C BOUNDARIES. DO 110 I=IMIN,IMAX XP = XIN(I) - XSING XP2 = XP*XP DTOP(I) = XP/(XP2 + YT2)*COEF2 DBOT(I) = XP/(XP2 + YB2)*COEF2 VTOP(I) = -ATAN2(YT,XP)*COEF1 VBOT(I) = -(ATAN2(YB,XP) + TWOPI)*COEF1 110 CONTINUE C COMPUTE DOUBLET AND VORTEX TERMS ON UPSTREAM C AND DOWNSTREAM BOUNDARIES DO 120 J=JMIN,JMAX YJ = YIN(J)*RTK YJ2 = YJ*YJ DUP(J) = XU/(XU2 + YJ2)*COEF2 DDOWN(J) = XD/(XD2 + YJ2)*COEF2 Q = PI - SIGN(PI,YJ) VUP(J) = -(ATAN2(YJ,XU) + Q)*COEF1 VDOWN(J) = -(ATAN2(YJ,XD) + Q)*COEF1 120 CONTINUE IF (AK .GT. 0.0) CALL ANGLE RETURN C COMPUTE WALL CONSTANTS FOR VARIOUS WIND TUNNEL C CONDITIONS. 200 CONTINUE C BCTYPE = 2 C SOLID WALL TUNNEL POR = 0. C SET CONSTANTS FOR DOUBLET SOLUTION B = .5 ALPHA0 = PI ALPHA1 = PI ALPHA2 = PI C SET CONSTANTS FOR VORTEX SOLUTION BETA0 = HALFPI BETA1 = HALFPI BETA2 = HALFPI GO TO 700 300 CONTINUE C BCTYPE = 3 C FREE JET F = 0. RTKPOR = 0. C SET CONSTANTS FOR DOUBLET SOLUTION ALPHA0 = HALFPI ALPHA1 = HALFPI ALPHA2 = HALFPI C SET CONSTANTS FOR VORTEX SOLUTION JET = .5 BETA0 = 0. BETA1 = 0. BETA2 = 0. GO TO 700 400 CONTINUE C BCTYPE = 4 C IDEAL SLOTTED WALL RTKPOR = 0. FHINV = 1.0 / ( F * H ) C SET CONSTANTS FOR DOUBLET SOLUTION CALL DROOTS C SET CONSTANTS FOR VORTEX SOLUTION JET = .5 CALL VROOTS GO TO 700 500 CONTINUE C BCTYPE = 5 C IDEAL PERFORATED/POROUS WALL F = 0. RTKPOR = RTK/POR C SET CONSTANTS FOR DOUBLET SOLUTION ALPHA0 = HALFPI - ATAN(-RTKPOR) ALPHA1 = ALPHA0 ALPHA2 = ALPHA0 C SET CONSTANTS FOR VORTEX SOLUTION BETA0 = ATAN(RTKPOR) BETA1 = BETA0 BETA2 = BETA1 GO TO 700 600 CONTINUE C BCTYPE = 6 C GENERAL HOMOGENEOUS WALL BOUNDARY CONDITION C BOUNDARY CONDITION IS NOT OPERABLE YET IN FINITE C DIFFERENCE SUBROUTINES. C FAR FIELD SOLUTION HAS BEEN DERIVED AND IS INCLUDED C HERE FOR FUTURE USE RTKPOR = RTK/POR CALL DROOTS CALL VROOTS C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,1000) WRITE(15,1000) 1000 FORMAT(35H1ABNORMAL STOP IN SUBROUTINE FARFLD/ 1 24H BCTYPE=6 IS NOT USEABLE) STOP 700 CONTINUE C COMPUTE FUNCTIONAL FORMS FOR UPSTREAM AND DOWNSTREAM C BOUNDARY CONDITIONS FOR DOUBLET AND VORTEX XU = (XIN(IMIN) - XSING)/(RTK*H) XD = (XIN(IMAX) - XSING)/(RTK*H) C DOUBLET TERMS COEF1 = .5/AK/H ARG0 = ALPHA0 ARG1 = PI - ALPHA1 ARG2 = TWOPI - ALPHA2 EXARG0 = EXP(-ARG0*XD) EXARG1 = EXP(ARG1*XU) EXARG2 = EXP(ARG2*XU) DO 720 J = JMIN,JMAX YJ = YIN(J) / H DDOWN(J) = COEF1*(B + OMEGA0*COS(YJ*ARG0)*EXARG0) DUP(J) = -COEF1*((1.-B)*OMEGA1*COS(YJ*ARG1)*EXARG1 + 1 OMEGA2*COS(YJ*ARG2)*EXARG2) 720 CONTINUE C VORTEX TERMS ARG0 = BETA0 ARG1 = PI + BETA1 ARG2 = PI - BETA2 EXARG0 = EXP(-ARG0*XD) EXARG1 = EXP(-ARG1*XD) EXARG2 = EXP(ARG2*XU) DO 740 J = JMIN,JMAX YJ = YIN(J) / H TERM = YJ IF (JET .EQ. 0.0) TERM = SIN(YJ*ARG0) / ARG0 VDOWN(J) = -.5*(1. - SIGN(1.,YJ) + (1. - JET)*PSI0*TERM*EXARG0 + 1 PSI1*SIN(YJ*ARG1)*EXARG1/ARG1) TERM = 0. IF (JET .NE. 0.0) TERM = JET * YJ / (1.0 + F) VUP(J) = -.5*(1. - TERM - PSI2*SIN(YJ*ARG2)*EXARG2/ARG2) 740 CONTINUE RETURN END SUBROUTINE FINDSK(ISTART,IEND,J, ISK) C SUBROUTINE LOCATES SHOCK WAVE ON LINE J BETWEEN C ISTART AND IEND. SHOCK IS LOCATED AT SHOCK POINT. C IF NO SHOCK FOUND, ISK IS SET NEGATIVE. C CALLED BY - CDCOLE. LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT ISK = ISTART - 1 U2 = PX(ISK,J) 5 ISK = ISK + 1 U1 = U2 U2 = PX(ISK,J) IF(U1 .GT. SONVEL .AND. U2 .LE. SONVEL) GO TO 10 IF(ISK .LT. IEND) GO TO 5 ISK = - IEND 10 CONTINUE RETURN END SUBROUTINE FIXPLT C C SETS UP ARRAYS FOR SUBROUTINE CPPLOT. C CALLED BY - PRINT. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON /COM13/ CDFACT , CLFACT , CMFACT , CPFACT , CPSTAR COMMON /COM25/ CPL(100) , CPU(100) , IDLA COMMON /COM30/ CPUP(101), CPLO(101) , CPS(101) , XP(101) C YMX = 5.*CPFACT YMN = -5.*CPFACT K = 0 DO 150 I = IMIN,IMAX K = K + 1 QCP = -CPU(I) QCP = AMAX1(QCP,YMN) QCP = AMIN1(QCP,YMX) CPUP(K) = QCP QC1 = -CPL(I) QC1 = AMAX1(QC1,YMN) QC1 = AMIN1(QC1,YMX) CPLO(K) = QC1 QC2 = -CPSTAR QC2 = AMAX1(QC2,YMN) QC2 = AMIN1(QC2,YMX) CPS(K) = QC2 XP(K) = X(I) 150 CONTINUE IMP = K + 1 CPUP(IMP) = YMX CPLO(IMP) = YMN CPS(IMP) = 0.0 XP(IMP) = X(IMAX) + .001 CALL CPPLOT (XP, CPUP, CPLO, CPS, IMP) RETURN END SUBROUTINE GUESSP C SUBROUTINE INITIALIZES P ARRAY AS FOLLOWS C PSTART = 1 P SET TO ZERO C PSTART = 2 P READ FROM UNIT 7 C PSTART = 3 P SET TO VALUES IN CORE C IF PSTART = 2 OR 3, SOLUTION IS INTERPOLATED FROM C XOLD, YOLD TO X,Y C BOUNDARY CONDITIONS FOR P ON OUTER BOUNDARIES ARE C AUTOMATICALLY SET DURING INITIALIZATION C CALLED BY - TSFOIL. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK LOGICAL ABORT1 COMMON / COM3/ IREF , ABORT1 , ICUT , KSTEP LOGICAL AMESH COMMON / COM4/ XIN(100) , YIN(100), AMESH INTEGER PSTART LOGICAL PSAVE COMMON /COM11/ ALPHAO , CLOLD , DELTAO , DUBO , EMACHO , 1 IMINO , IMAXO , IMAXI , JMINO , JMAXO , 2 JMAXI , PSAVE , PSTART ,TITLE(20),TITLEO(20), 3 VOLO , XOLD(100),YOLD(100) COMMON /COM13/ CDFACT , CLFACT , CMFACT , CPFACT , CPSTAR INTEGER BCTYPE COMMON /COM28/ BCTYPE , CIRCFF , FHINV , POR , CIRCTE COMMON /COM30/ PT(100) , REST(304) C C BRANCH TO APPROPRIATE LOCATION GO TO ( 200, 300, 400) , PSTART 200 CONTINUE C PSTART = 1 C P SET TO ZERO DO 205 I=1,101 DO 205 J=1,102 P(J,I) = 0. 205 CONTINUE DUB = 0. CIRCFF = 0. CIRCTE=0.0 RETURN 300 CONTINUE C PSTART = 2 C P, X, Y ARRAYS READ FROM UNIT 7 IN SUBROUTINE READIN C TOGETHER WITH INFORMATION ABOUT OLD SOLUTION 400 CONTINUE C PSTART = 3 C ARRAYS FROM PREVIOUS CASE ARE ALREADY IN P,XOLD,YOLD DUB = DUBO CIRCFF = CLOLD/CLFACT CIRCTE=CIRCFF C C FOR PSTART = 2 OR 3, OLD P ARRAY ON XOLD, YOLD MESH C MUST BE INTERPOLATED ONTO NEW X Y MESH. C C INTERPOLATE P FROM XOLD,YOLD TO X,YOLD C CHECK TO SEE IF XOLD AND XIN ARE THE SAME MESH IF(IMAXI .NE. IMAXO) GO TO 450 DO 410 I=IMIN, IMAXI TEST = ABS(XIN(I) - XOLD(I)) IF(TEST .GT. .0001) GO TO 450 410 CONTINUE C XIN AND XOLD ARE SAME MESH. C P ARRAY MAY BE INTERPOLATED BY SIMPLE DELETION OF C VALUES AT MESH POINTS DELETED IN SUBROUTINE CUTOUT C IF IREF .LE. ZERO, NO INTERPOLATION IS NEEDED IF(IREF .LE. 0) GO TO 500 ISTEP = 2*IREF DO 430 J = JMINO,JMAXO INEW = 0 DO 420 I = IMINO,IMAXO,ISTEP INEW = INEW + 1 P(J,INEW) = P(J,I) 420 CONTINUE 430 CONTINUE C INTERPOLATION IN X DIRECTION COMPLETE IF XIN AND C XOLD ARE THE SAME. GO TO 500 450 CONTINUE C INTERPOLATE FROM XOLD TO X FOR ARBITRARY CASE DO 490 J = JMINO, JMAXO YP = YOLD(J) DO 475 I = IMIN, IMAX XP = X(I) IF(XP .LT. XOLD(IMINO)) GO TO 470 IF(XP .GT. XOLD(IMAXO)) GO TO 470 C NEW X MESH POINT WITHIN RANGE OF OLD X MESH C FIND VALUE OF XOLD .GT. XP X2 = XOLD(1) K = 0 455 K = K +1 X1 = X2 X2 = XOLD(K) IF (X2 .LT. XP) GO TO 455 IF (X2 .GT. XP) GO TO 460 PT(I) = P(J,K) GO TO 475 460 CONTINUE P1 = P(J,K-1) P2 = P(J,K) PT(I) = P1 + (P2 - P1)/(X2 - X1) * (XP - X1) GO TO 475 470 CONTINUE C NEW X MESH POINT IS OUTSIDE RANGE OF OLD X MESH C FOR SUPERSONIC FREESTREAM SET P=0, FOR SUBSONIC C FREESTREAM, EXTRAPOLATE USING FAR FIELD SOLUTION PT(I) = 0. IF(AK .GT. 0.) CALL EXTRAP(XP,YP,PT(I)) 475 CONTINUE C WRITE NEW VALUES FOR P INTO P ARRAY DO 480 I = IMIN, IMAX P(J,I) = PT(I) 480 CONTINUE 490 CONTINUE 500 CONTINUE C INTERPOLATE FROM X,YOLD TO X,Y C CHECK TO SEE IF YIN AND YOLD ARE THE SAME MESH C IF(JMAXI .NE. JMAXO) GO TO 550 DO 510 J = JMIN, JMAXI TEST = ABS(YIN(J) - YOLD(J)) IF(TEST .GT. .0001) GO TO 550 510 CONTINUE C YIN AND YOLD ARE THE SAME MESH C P ARRAY MAY BE INTERPOLATED BY SIMPLE DELETION OF C VALUES AT MESH POINTS DELETED IN SUBROUTINE CUTOUT C IF IREF .LE. ZERO, NO INTERPOLATION IS NEEDED IF(IREF .LE. 0) GO TO 600 JSTEP = 2*IREF DO 530 I= IMIN, IMAX JNEW = 0 DO 520 J = JMINO, JMAXO, JSTEP JNEW = JNEW + 1 P(JNEW,I) = P(J,I) 520 CONTINUE 530 CONTINUE C INTERPOLATION IN Y DIRECTION COMPLETE IF YIN AND C YOLD ARE THE SAME. GO TO 600 550 CONTINUE C INTERPOLATE YOLD TO Y FOR ARBITRARY CASE DO 590 I = IMIN,IMAX XP = X(I) K = 2 Y1 = YOLD(1) DO 575 J = JMIN,JMAX YP = Y(J) IF(YP .LT. YOLD(JMINO)) GO TO 570 IF(YP .GT. YOLD(JMAXO)) GO TO 571 C NEW Y MESH POINT WITHIN RANGE OF OLD Y MESH C FIND VALUE OF YOLD .GT. YP Y2 = Y1 K = K - 1 555 K = K + 1 Y1 = Y2 Y2 = YOLD(K) IF(Y2 .LE. YP) GO TO 555 P1 = P(K-1,I) P2 = P(K,I) PT(J) = P1 + (P2 - P1) / (Y2 - Y1) * (YP - Y1) GO TO 575 C NEW Y MESH POINT OUTSIDE RANGE OF OLD Y MESH. C SET P(J,I) = P(JMINO,I) OR, IF SUBSONIC FREESTREAM C FREE AIRFLOW, EXTRAPOLATE FOR P(J,I) USING FARFLD C FORMULA. 570 CONTINUE PT(J) = P(JMINO,I) GO TO 572 571 PT(J) = P(JMAXO,I) 572 CONTINUE IF(AK .GT. 0. .AND. BCTYPE .EQ. 1) CALL EXTRAP(XP,YP,PT(J)) 575 CONTINUE C PUT NEW P VALUES INTO P ARRAY DO 580 J = JMIN,JMAX P(J,I) = PT(J) 580 CONTINUE 590 CONTINUE 600 CONTINUE RETURN END SUBROUTINE INPERR (I) C C FATAL ERROR IN INPUT. WRITE MESSAGE AND STOP. C CALLED BY - READIN, SCALE. C C Changed to write to a file - Andy Ko 4/3/03 C IF (I .EQ. 1) WRITE (6,901) C IF (I .EQ. 2) WRITE (6,902) C IF (I .EQ. 3) WRITE (6,903) C IF (I .EQ. 4) WRITE (6,904) C IF (I .EQ. 5) WRITE (6,905) C IF (I .EQ. 6) WRITE (6,906) C IF (I .EQ. 7) WRITE (6,907) C IF (I .EQ. 8) WRITE(6,908) IF (I .EQ. 1) WRITE (15,901) IF (I .EQ. 2) WRITE (15,902) IF (I .EQ. 3) WRITE (15,903) IF (I .EQ. 4) WRITE (15,904) IF (I .EQ. 5) WRITE (15,905) IF (I .EQ. 6) WRITE (15,906) IF (I .EQ. 7) WRITE (15,907) IF (I .EQ. 8) WRITE (15,908) STOP 901 FORMAT(1H0,5X,45HIMAX OR JMAX IS GREATER THAN 100,NOT ALLOWED.) 902 FORMAT(1H0,5X,39HX MESH POINTS NOT MONOTONIC INCREASING.) 903 FORMAT(1H0,5X,39HY MESH POINTS NOT MONOTONIC INCREASING.) 904 FORMAT(1H0,5X,44HMACH NUMBER NOT IN PERMITTED RANGE. (.5,2.0)) 905 FORMAT(1H0,5X,41HALPHA NOT IN PERMITTED RANGE. (-9.0, 9.0)) 906 FORMAT(1H0,5X,41HDELTA NOT IN PERMITTED RANGE. ( 0.0, 1.0)) 907 FORMAT(1H0,5X,45HAK=0. VALUE OF AK MUST BE INPUT SINCE PHYS=F.) 908 FORMAT(1H0,5X,60HMACH NUMBER IS NOT LESS THAN 1.0 FOR VISCOUS WEDG 1E INCLUSION ) END SUBROUTINE ISLIT ( X ) C C ISLIT COMPUTES ILE AND ITE FOR ANY XMESH ARRAY. C CALLED BY - CKMESH, CUTOUT, READIN, REFINE. C COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 C DIMENSION X(100) I = IMIN - 1 10 CONTINUE I = I + 1 IF (X(I) .LT. 0.0) GO TO 10 ILE = I 20 CONTINUE I = I + 1 IF (X(I) .LE. 1.0) GO TO 20 ITE = I - 1 RETURN END SUBROUTINE JSLIT ( Y ) C C JSLIT COMPUTES JLOW AND JUP FROM Y ARRAY. C CALLED BY - CKMESH, CUTOUT, READIN, REFINE. C COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 C DIMENSION Y(100) J = JMIN - 1 10 CONTINUE J = J + 1 IF (Y(J) .LT. 0.0) GO TO 10 JLOW = J - 1 JUP = J RETURN END REAL FUNCTION LIFT (CLFACT) C COMPUTE LIFT COEFFICIENT FROM JUMP IN P AT TRAILING C EDGE OF AIRFOIL C CALLED BY - PRINT1, SOLVE. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM7/ CJUP , CJUP1 , CJLOW , CJLOW1 C PTOP = CJUP*P(JUP,ITE) - CJUP1*P(JUP + 1 , ITE) PBOT = CJLOW*P(JLOW,ITE) - CJLOW1*P(JLOW-1,ITE) LIFT = 2.*CLFACT*(PTOP-PBOT) RETURN END SUBROUTINE MACHMP C C SUBROUTINE TO PRINT MAP OF MACH NO. ROUNDED C TO NEAREST .1. CALLED BY - PRINT. C COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON /COM30/ MM(100) , REST(304) C DATA IB/1H / , IP/1H+/ , IM/1H-/ , IL/1HL/ , IT/1HT/ C C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,400) WRITE (15,400) DO 60 K=2,JMAX J = JMAX - K + 2 IJC = IB IF ( J .EQ. JUP) IJC = IP IF ( J .EQ. JLOW ) IJC = IM DO 10 I=1,IMAX MM(I) = IB 10 CONTINUE DO 50 I=2,IMAX U = PX(I,J) EM = EMACH1 (U) IF (EM .GT. 0.0) GO TO 20 MM(I) = 0 GO TO 40 20 CONTINUE IF (EM .LE. 1.05) GO TO 30 EM = EM - 1.0 GO TO 20 30 CONTINUE MM(I) = INT(10. * EM + .5) 40 CONTINUE 50 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,410) IJC, (MM(I),I=2,IMAX) WRITE (15,410) IJC, (MM(I),I=2,IMAX) 60 CONTINUE DO 70 I=1,IMAX MM(I) = IB IF (I .EQ. ILE) MM(I) = IL IF (I .EQ. ITE) MM(I) = IT 70 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,420) (MM(I),I=2,IMAX) WRITE (15,420) (MM(I),I=2,IMAX) RETURN C 400 FORMAT(40H1 MACH NO. MAP. ROUNDED TO NEAREST .1///) 410 FORMAT (11X,A1,99I1) 420 FORMAT (12X,99A1) C END SUBROUTINE M1LINE C PRINTS COORDINATES WHERE SONIC VELOCITY IS COMPUTED C LINEAR INTERPOLATION BETWEEN MESH POINTS IS USED C CALLED BY - PRINT. C ********************************************************************** C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT INTEGER BCTYPE COMMON /COM28/ BCTYPE , CIRCFF , FHINV , POR , CIRCTE COMMON /COM30/ XSLPRT(200) , YSLPRT(200) , REST(4) C DIMENSION XSONIC(10) NPTS = 0 KMIN = JMIN KMAX = JMAX JP = JMAX + JMIN DO 200 K=KMIN,KMAX J = JP - K YPR = YFACT*Y(J) PX2 = PX(IMIN,J) M = 0 IF (J .NE. JLOW) GO TO 170 C Changed to write to a file - Andy Ko 4/3/03 C IF (NPTS .NE. 0) WRITE (6,2070) IF (NPTS .NE. 0) WRITE (15,2070) 170 CONTINUE IMM = IMIN + 1 DO 180 I=IMM,IMAX PX1 = PX2 PX2 = PX(I,J) IF(PX1 .GT. SONVEL .AND. PX2 .GT. SONVEL) GO TO 180 IF(PX1 .LT. SONVEL .AND. PX2 .LT. SONVEL) GO TO 180 C Changed to write to a file - Andy Ko 4/3/03 C IF (NPTS .EQ. 0) WRITE (6,2000) IF (NPTS .EQ. 0) WRITE (15,2000) M = M+1 RATIO = (SONVEL-PX1)/(PX2-PX1) XSONIC(M) = X(I-1) + (X(I)-X(I-1))*RATIO NPTS = NPTS + 1 XSLPRT(NPTS) = XSONIC(M) YSLPRT(NPTS) = YPR IF(NPTS .GE. 200) GO TO 250 180 CONTINUE IF(M .EQ. 0) GO TO 200 C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,2040) YPR,(XSONIC(L),L=1,M) WRITE(15,2040) YPR,(XSONIC(L),L=1,M) 200 CONTINUE GO TO 300 250 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,2060) WRITE(15,2060) GO TO 400 300 CONTINUE IF(NPTS .EQ.0) GO TO 400 YM = Y(JMIN) YX = Y(JMAX) DO 320 N=1,NPTS IF (YSLPRT(N) .NE. YM .AND. YSLPRT(N) .NE. YX) GO TO 320 C Changed to write to a file - Andy Ko 4/3/03 C IF ( AK .GT. 0.0 ) WRITE (6,2050) C IF ( AK .LT. 0.0 .AND. BCTYPE .EQ. 1 ) WRITE (6,2050) IF ( AK .GT. 0.0 ) WRITE (15,2050) IF ( AK .LT. 0.0 .AND. BCTYPE .EQ. 1 ) WRITE (15,2050) 320 CONTINUE XMIN = -.75 XMAX = 1.75 XINCR = .25 YMIN = -1.0 YMAX = 1.5 YINCR = .5 CALL PLTSON(XSLPRT,YSLPRT,XMIN,XMAX,XINCR,YMIN,YMAX,YINCR,NPTS) 400 CONTINUE RETURN C 2000 FORMAT(23H1SONIC LINE COORDINATES/6X,1HY,10X,6HXSONIC//) 2040 FORMAT(11F10.5) 2050 FORMAT(20H0***** CAUTION *****/ 1 34H SONIC LINE HAS REACHED A BOUNDARY/ 2 61H THIS VIOLATES ASSUMPTIONS USED TO DERIVE BOUNDARY CONDITIONS/ 3 29H SOLUTION IS PROBABLY INVALID) 2060 FORMAT(20H0***** CAUTION *****/36H NUMBER OF SONIC POINTS EXDEEDED 1 200/25H ARRAY DIMENSION EXCEEDED/42H EXECUTION OF SUBROUTINE M1LI 2NE TERMINATED) 2070 FORMAT(2X,13HBODY LOCATION) C END SUBROUTINE NEWISK(ISKOLD,J,ISKNEW) C FIND NEW LOCATION OF SHOCKWAVE (ISKNEW) ON LINE J C GIVEN AN INITIAL GUESS FOR LOCATION (ISKOLD). C SHOCK LOCATION IS DEFINED AS LOCATION OF SHOCK POINT. C IF NO SHOCK IS FOUND, ISKNEW IS SET NEGATIVE. C CALLED BY - CDCOLE. C LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT I2 = ISKOLD + 2 ISKNEW = ISKOLD - 3 U2 = PX(ISKNEW,J) 5 ISKNEW = ISKNEW + 1 U1 = U2 U2 = PX(ISKNEW,J) IF(U1 .GT. SONVEL .AND. U2 .LE. SONVEL) GO TO 10 IF(ISKNEW .LT. I2) GO TO 5 C NO SHOCK POINT FOUND, TIP OF SHOCK REACHED. ISKNEW = - ISKNEW 10 RETURN END FUNCTION PITCH(CMFACT) C COMPUTE AIRFOIL PITCHING MOMENT ABOUT X = XM, Y=0 C THE INTEGRAL OF(X-XM)(CPUPPER-CPLOWER) IS INTEGRATED C BY PARTS TO GIVE AN INTEGRAL IN THE JUMP IN P PLUS A C CONSTANT EQUAL TO THE JUMP IN P AT X = 1. TIMES THE C MOMENT ARM AT THE TRAILING EDGE C CALLED BY - PRINT1, SOLVE. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM7/ CJUP , CJUP1 , CJLOW , CJLOW1 COMMON /COM30/ XI(100) , ARG(100) , REST(204) C SET XM TO QUARTER CHORD C XM = .25 K = 0 DO 10 I=ILE,ITE K = K + 1 PTOP = CJUP*P(JUP,I) - CJUP1*P(JUP+1,I) PBOT = CJLOW*P(JLOW,I) - CJLOW1*P(JLOW-1,I) ARG(K) = PTOP - PBOT XI(K) = X(I) 10 CONTINUE CALL TRAP(XI,ARG,K,SUM) PITCH = CMFACT*((1.-XM)*ARG(K) - SUM) * (-2.) RETURN END SUBROUTINE PLTSON (X,Y,XAXMIN,XMAX,XINCR,YAXMIN,YMAX,YINCR,NPTS) C C THIS IS A MODIFIED VERSION OF PRNPLT AND C XMAX, XINCR, YMAX, YINCR MUST BE SUPPLIED. C XMIN, AND YMIN MUST ALSO BE SUPPLIED. C PRINTER PLOT ROUTINE M. S. ITZKOWITZ MAY,1967 C PLOTS THE @NPTS@ POINTS GIVEN BY, X(I),Y(I) C ON A 51 X 101 GRID USING A TOTAL OF 56 LINES ON C THE PRINTER. C IF EITHER INCREMENTAL STEP SIZE IS ZERO, THE C PROGRAM EXITS. C NEITHER OF THE INPUT ARRAYS ARE DESTROYED. C CALLED BY - M1LINE. C DIMENSION X(NPTS), Y(NPTS), IGRID(105), XAXIS(11) DIMENSION IBS(9) C INTEGER BLANK , DOT , STAR , DASH DATA BLANK , DOT , STAR , DASH / 1H , 1H., 1H*, 1H-/ DATA IBS/1HB,1HO,1HD,1HY,1H ,1HS,1HL,1HI,1HT/ C C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,900) C WRITE(6,910) WRITE(15,900) WRITE(15,910) IF (XINCR .EQ. 0. .OR. YINCR .EQ. 0.) GO TO 800 YRNGE = YMAX - YAXMIN XRNGE = XMAX - XAXMIN YVAL = YRNGE * .02 XVAL = XRNGE * .01 JSLT = IFIX( 51.*(YMAX/YRNGE) ) + 1 OXV = 1.0 / XVAL OYV = 1.0 / YVAL IZERO = YMAX * OYV + 1.5 JZERO = 103.5 - XMAX * OXV IF (JZERO .GT. 103 .OR. JZERO .LT. 4) JZERO = 2 FIZERO = IZERO FJZERO = JZERO IGRID(1) = BLANK IGRID(2) = DOT IGRID(104) = DOT IGRID(105) = BLANK FIZER5 = FIZERO + .5 FJZER5 = FJZERO + .5 C NOPW IS THE NUMBER OF PRINT WHEELS USED TO C SPAN THE BODY LENGTH. NOPW = OXV C JS - THE NUMBER OF PRT WHEEL POSITIONS TO BE C FILLED AT EACH END OF THE WORD ?BODY SLIT? JS = 0 IF (NOPW .GT. 9) JS = (NOPW-9) / 2 C C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,905) C WRITE (6,903) WRITE (15,905) WRITE (15,903) C LOOP TO SET UP ONE LINE EACH TIME THRU. DO 30 I=1,51 C BLANK THE LINE TO BE PRINTED. DO 2 J=3,103 IGRID(J) = BLANK 2 CONTINUE C C SEARCH ARRAY FOR POINTS ON THIS LINE. C DO 10 K=1,NPTS ITEST = FIZER5 - Y(K) * OYV IF (ITEST .NE. I) GO TO 7 J = FJZER5 + X(K) * OXV IF (J .GT. 103) J = 103 IF (J .LT. 3) J = 3 IGRID(J) = STAR 7 CONTINUE IF (JSLT .EQ. I) GO TO 9 GO TO 10 9 CONTINUE C WRITE ?BODY SLIT? IF THIS IS THE J0 LINE. J = JZERO + 1 IF (JS .EQ. 0) GO TO 4 DO 3 L=1,JS IF (IGRID(J) .NE. STAR) IGRID(J) = DASH J = J + 1 3 CONTINUE 4 CONTINUE DO 5 L=1,9 IF (IGRID(J) .NE. STAR) IGRID(J) = IBS(L) J = J + 1 5 CONTINUE IF (JS .EQ. 0) GO TO 10 DO 6 L=1,JS IF (IGRID(J) .NE. STAR) IGRID(J) = DASH J = J + 1 6 CONTINUE 10 CONTINUE IF (MOD(I,10) .EQ. 1) GO TO 13 C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,901) IGRID WRITE (15,901) IGRID GO TO 30 13 CONTINUE FI = I - 1 YAXIS = YMAX - FI * YINCR * .1 IF (ABS(YAXIS) .LT. YAXMIN) YAXIS = 0. C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,902) YAXIS, (IGRID(J),J=1,105) WRITE (15,902) YAXIS, (IGRID(J),J=1,105) 30 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,903) C WRITE (6,905) WRITE (15,903) WRITE (15,905) DO 40 M=1,11 FM = 11 - M XAXIS(M) = XMAX - XINCR * FM IF (XAXIS(M) .LT. XAXMIN) XAXIS(M) = XAXMIN 40 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,904) XAXIS, NPTS C WRITE(6,900) WRITE (15,904) XAXIS, NPTS WRITE(15,900) RETURN C C Changed to write to a file - Andy Ko 4/3/03 C 800 WRITE (6,9800) 800 WRITE (15,9800) STOP C 900 FORMAT(1H1) 901 FORMAT(14X,105A1) 902 FORMAT(1X,F10.1,2X,1H+,105A1,1H+) 903 FORMAT(15X,103(1H.)) 904 FORMAT(7X,11(F10.2),2H (,I4,5H PTS) ) 905 FORMAT(16X,11(1H+,9X)) 910 FORMAT(35X,15HSONIC LINE PLOT,10X,6HY VS X,10X,20H * FOR SONIC PO 1INTS//) 9800 FORMAT(46H1SCALING ERROR IN PRNPLT, EXECUTION TERMINATED ) C END SUBROUTINE PRBODY C PRINTS OUT BODY GEOMETRY C IF PHYS = .TRUE. ALL DIMENSIONS ARE NORMALIZED BY C AIRFOIL CHORD C IF PHYS = .FALSE. ALL DIMENSIONS EXCEPT X ARE C NORMALIZED BY CHORD LENGTH AND C THICKNESS RATIO. C CALLED BY - BODY. C COMMON / COM6/ FL(100) , FXL(100), FU(100) , FXU(100), 1 CAMBER(100), THICK(100),VOL , XFOIL(100), IFOIL LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT C C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,1000) WRITE(15,1000) C FIND MAXIMUM THICKNESS AND CAMBER THMAX = 0. CAMAX = 0. DO 10 I = 1,IFOIL THMAX = AMAX1(THMAX,THICK(I)) CAMAX = AMAX1(CAMAX,CAMBER(I)) 10 CONTINUE THMAX = 2.*THMAX IF(PHYS) GO TO 100 C PRINTOUT IN SIMILARITY VARIABLES C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,1001) THMAX C WRITE(6,1002) VOL,CAMAX C WRITE(6,1003) C DO 40 I=1,IFOIL C WRITE (6,1005) XFOIL(I),FU(I),FXU(I),FL(I),FXL(I),THICK(I), C 1 CAMBER(I) WRITE(15,1001) THMAX WRITE(15,1002) VOL,CAMAX WRITE(15,1003) DO 40 I=1,IFOIL WRITE (15,1005) XFOIL(I),FU(I),FXU(I),FL(I),FXL(I),THICK(I), 1 CAMBER(I) 40 CONTINUE RETURN 100 CONTINUE C PRINTOUT IN PHYSICAL VARIALBES THMAX = DELTA*THMAX CAMAX = DELTA*CAMAX C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,1006) THMAX C VOLUME = VOL*DELTA C WRITE(6,1002) VOLUME,CAMAX C WRITE(6,1004) WRITE(15,1006) THMAX VOLUME = VOL*DELTA WRITE(15,1002) VOLUME,CAMAX WRITE(15,1004) DO 110 I = 1,IFOIL YUP = DELTA*FU(I) YXUP = DELTA*FXU(I) YLO = DELTA*FL(I) YXLO = DELTA*FXL(I) TH = DELTA*THICK(I) CA = DELTA*CAMBER(I) C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,1005) XFOIL(I),YUP,YXUP,YLO,YXLO,TH,CA WRITE(15,1005) XFOIL(I),YUP,YXUP,YLO,YXLO,TH,CA 110 CONTINUE RETURN 1000 FORMAT(1x,'AIRFOIL GEOMETRY INFORMATION') 1001 FORMAT(1x,'PRINTOUT IN SIMILARITY VARIABLES', 1 28X, 15HMAX THICKNESS =,F12.8) 1002 FORMAT(1x,'AIRFOIL VOLUME =',F12.8,32X,15HMAX CAMBER =,F10.6// 1 20X,14HUPPER SURFACE,14X,14HLOWER SURFACE) 1003 FORMAT(8X,1HX,1X,2(12X,1HF,9X,5HDF/DX),8X,9HTHICKNESS,4X, 1 6HCAMBER) 1004 FORMAT(8X,1HX,1X,2(12X,1HY,9X,5HDY/DX),8X,9HTHICKNESS,4X, 1 6HCAMBER) 1005 FORMAT(1X,F12.8,2X,2F12.8,2(3X,2F12.8)) 1006 FORMAT(1x,'PRINTOUT IN PHYSICAL VARIABLES NORMALIZED BY CHORD ', 1 'LENGTH', 3X, 15HMAX THICKNESS =,F10.6) END SUBROUTINE PRINT C C SUBROUTINE FOR MAIN OUTPUT PRINT CONTROL. C PRINTS RELATIVE PARAMETERS AND CALLS C SPECIALIZED PRINT/PLOT SUBROUTINES AS REQUIRED. C C PRINTOUT DESCRIPTION C IF PHYS IS .TRUE. THE FOLLOWING INFORMATION IS C PRINTED IN PHYSICALLY SCALED VARIABLES (I.E. USUAL C AERODYNAMIC TERMS) C X,Y ARE COORDINATES NON-DIMENSIONALIZED BY C AIRFOIL CHORD. C CP = (P - P(INF))/Q(INF) = -2.*(DP/DX) C CD = DRAG COEFF = D/Q(INF)*AIRFOIL CHORD C CPSTAR = CRITICAL PRESSURE COEFFICIENT C IF PHYS IS .FALSE. THE ABOVE INFORMATION IS PRINTED C IN TRANSONIC SIMILARITY VARIABLES. C C PRINTOUT OF SONIC LINE ORDINATES IN SUBR. M1LINE C INCLUDING PRINTPLOT OF SONIC LINE ORDINATES C CALLED BY - TSFOIL. C ********************************************************************** C COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK LOGICAL ABORT1 COMMON / COM3/ IREF , ABORT1 , ICUT , KSTEP INTEGER PSTART LOGICAL PSAVE COMMON /COM11/ ALPHAO , CLOLD , DELTAO , DUBO , EMACHO , 1 IMINO , IMAXO , IMAXI , JMINO , JMAXO , 2 JMAXI , PSAVE , PSTART ,TITLE(20),TITLEO(20), 3 VOLO , XOLD(100),YOLD(100) COMMON /COM13/ CDFACT , CLFACT , CMFACT , CPFACT , CPSTAR LOGICAL FCR , KUTTA COMMON /COM14/ CLSET , FCR , KUTTA , WCIRC LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT INTEGER BCTYPE COMMON /COM28/ BCTYPE , CIRCFF , FHINV , POR , CIRCTE C DIMENSION TPH(6) , SIM(8) , BCT(15) , FCP(14) C DATA TPH /4HSIMI,4HLARI,2HTY,4H PHY,4HSICA,2HL / DATA SIM /4HCOLE,4H ,4HSPRE,4HITER,4HKRUP,4HP , 1 4HUSER,4H / DATA BCT /4HFREE,4H AIR,4H ,4HSOLI,4HD WA,4HLL ,4HFREE, 1 4H JET,4H ,4HSLOT,4HTED ,4HWALL,4HPORO,4HUS W,4HALL / DATA FCP /4HFULL,4HY CO,4HNSER,4HVATI,4HVE. ,4H ,4H , 1 4HNOT ,4HCONS,4HERVA,4HTIVE,4H AT ,4HSHOC,4HK. / C C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,900) WRITE (15,900) IS = 1 IF (PHYS) IS = 4 IE = IS + 2 C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,901) (TPH(I),I=IS,IE) WRITE (15,901) (TPH(I),I=IS,IE) IE = 2 * SIMDEF IS = IE - 1 C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,902) (SIM(I),I=IS,IE) WRITE (15,902) (SIM(I),I=IS,IE) IE = 3 * BCTYPE IS = IE - 2 C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,903) (BCT(I),I=IS,IE) WRITE (15,903) (BCT(I),I=IS,IE) IS = 8 IF (FCR) IS = 1 IE = IS + 6 C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,904) (FCP(I),I=IS,IE) WRITE (15,904) (FCP(I),I=IS,IE) IF (KUTTA) GO TO 10 C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,905) WRITE (15,905) GO TO 20 10 CONTINUE C C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,906) WRITE (15,906) 20 CONTINUE ALPHVF = ALPHA * VFACT C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,911) C IF ( PHYS ) WRITE (6,909) EMACH, DELTA C WRITE (6,907) ALPHVF, AK C IF (AK .GT. 0.0) WRITE (6,908) DUB C IF ( PHYS ) WRITE(6,910)CPFACT,CDFACT,CMFACT,CLFACT,YFACT,VFACT WRITE (15,911) IF ( PHYS ) WRITE (15,909) EMACH, DELTA WRITE (15,907) ALPHVF, AK IF (AK .GT. 0.0) WRITE (15,908) DUB IF ( PHYS ) WRITE(15,910)CPFACT,CDFACT,CMFACT,CLFACT,YFACT,VFACT CALL PRINT1 CALL PRTMC CALL MACHMP IF ( ABORT1 ) GO TO 60 CALL FIXPLT IF (BCTYPE .EQ. 1) GO TO 40 IF (BCTYPE .EQ. 3) GO TO 40 CALL PRTWAL 40 CONTINUE CALL M1LINE IF (PRTFLO .EQ. 1) GO TO 50 CALL PRTFLD 50 CONTINUE CALL CDCOLE 60 CONTINUE RETURN C 900 FORMAT(1H1) 901 FORMAT(14H0 PRINTOUT IN ,2A4,A2,11H VARIABLES.) 902 FORMAT(41H0 DEFINITION OF SIMILARITY PARAMETERS BY ,2A4) 903 FORMAT(25H0 BOUNDARY CONDITION FOR ,3A4) 904 FORMAT(27H0 DIFFERENCE EQUATIONS ARE ,7A4) 905 FORMAT(37H0 LIFT COEFFICIENT SPECIFIED BY USER.) 906 FORMAT(30H0 KUTTA CONDITION IS ENFORCED.) 907 FORMAT(13X,7HALPHA =,F12.7/17X,3HK =,F12.7) 908 FORMAT(2X,18HDOUBLET STRENGTH =,F12.7) 909 FORMAT(14X,6HMACH =,F12.7/13X,7HDELTA =,F12.7) 910 FORMAT(1H0,8X,38HPARAMETERS USED TO TRANSFORM VARIABLES/ 1 18X,20HTO TRANSONIC SCALING//12X,8HCPFACT =,F12.7/ 2 12X,8HCDFACT =,F12.7/12X,8HCMFACT =,F12.7/12X,8HCLFACT =,F12.7 3 /13X,7HYFACT =,F12.7/ 13X,7HVFACT =,F12.7) 911 FORMAT(1H0) C END SUBROUTINE PRINT1 C C PRINTS PRESSURE COEFFICIENT AND MACH NUMBER C ON Y=0 LINE, AND PLOTS CP ALONG SIDE OF PRINT. C CALLED BY - TSFOIL, PRINT. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 LOGICAL ABORT1 COMMON / COM3/ IREF , ABORT1 , ICUT , KSTEP COMMON / COM7/ CJUP , CJUP1 , CJLOW , CJLOW1 COMMON /COM13/ CDFACT , CLFACT , CMFACT , CPFACT , CPSTAR COMMON /COM25/ CPL(100) , CPU(100) , IDLA LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT COMMON /COM30/ EM1L(100), EM1U(100), YM(100), REST(104) C REAL LIFT DIMENSION LINE1(60) , TMAC(2) C DATA TMAC / 2HM1, 2HK1/ , IB /1H /, IL / 1HL/,IU /1HU/, IS/1H*/ DATA IBB/ 1HB/ C CL = LIFT (CLFACT) CM = PITCH (CMFACT) CPMIN = 1.0E37 CPMAX = -CPMIN IEM = 0 CJ01 = -Y(JLOW)/(Y(JUP)-Y(JLOW)) CJ02 = Y(JUP)/(Y(JUP)-Y(JLOW)) C DO 10 I=IMIN,IMAX UL = CJLOW*PX(I,JLOW) - CJLOW1*PX(I,JLOW-1) IF(I .GT. ITE) UL=CJ01*PX(I,JUP) + CJ02*PX(I,JLOW) IF(I .LT. ILE) UL=CJ01*PX(I,JUP) + CJ02*PX(I,JLOW) CPL(I) = -2.0 * UL * CPFACT EM1L(I) = EMACH1(UL) IF (EM1L(I) .GT. 1.3) IEM = 1 UU = CJUP *PX(I,JUP) - CJUP1 *PX(I,JUP+1) IF (I .GT. ITE) UU = UL IF (I .LT. ILE) UU = UL CPU(I) = -2.0 * UU * CPFACT EM1U(I) = EMACH1(UU) IF (EM1U(I) .GT. 1.3) IEM = 1 CPMAX = AMAX1(CPMAX,CPU(I), CPL(I)) CPMIN = AMIN1(CPMIN,CPU(I), CPL(I)) 10 CONTINUE C CPLARG = AMAX1(CPMAX,ABS(CPMIN)) UNPCOL = CPLARG / 29. C C LOCATE CP* FOR PRINTER PLOT C COL = -CPSTAR / UNPCOL NCOL = SIGN((ABS(COL)+.5), COL) NCOLS = NCOL + 30 C C PRINT SINGLE VARIABLES C C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,260) C IF (IREF .EQ. 2) WRITE (6,270) C IF (IREF .EQ. 1) WRITE (6,280) C IF (IREF .EQ. 0) WRITE (6,290) C WRITE (6,200) CL, CM, CPSTAR WRITE (15,260) IF (IREF .EQ. 2) WRITE (15,270) IF (IREF .EQ. 1) WRITE (15,280) IF (IREF .EQ. 0) WRITE (15,290) WRITE (15,200) CL, CM, CPSTAR IF (CPL(IMIN) .LT. CPSTAR .AND. CPL(IMIN+1) .GT. CPSTAR) GO TO 70 15 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,210) WRITE (15,210) KT = 2 IF ( PHYS ) KT = 1 C C PRINT COLUMN HEADERS. C C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,220) TMAC(KT) , TMAC(KT) WRITE (15,220) TMAC(KT) , TMAC(KT) C iout = 1 iplot = 0 if(iref .eq. 0) write(iout,340) EMACH,CL DO 40 I=IMIN,IMAX C DO 30 K=1,60 LINE1(K) = IB 30 CONTINUE C COL = -CPU(I) / UNPCOL NCOL = SIGN((ABS(COL) + .5) , COL) NCOLU = NCOL + 30 LINE1(NCOLU) = IU C COL = -CPL(I) / UNPCOL NCOL = SIGN((ABS(COL) + .5) , COL) NCOLL = NCOL + 30 LINE1(NCOLL) = IL IF (NCOLL .EQ. NCOLU) LINE1(NCOLL) = IBB IF (IABS(NCOLS) .LT. 61) LINE1(NCOLS) = IS C C Changed to write to a file - Andy Ko 4/3/03 C IF (I .EQ. ILE) WRITE (6,230) C WRITE (6,250) I, X(I), CPL(I), EM1L(I), CPU(I), EM1U(I), LINE1 C IF (I .EQ. ITE) WRITE (6,240) IF (I .EQ. ILE) WRITE (15,230) WRITE (15,250) I, X(I), CPL(I), EM1L(I), CPU(I), EM1U(I), LINE1 IF (I .EQ. ITE) WRITE (15,240) C C Save DATA FOR PLOTTING C IF (IREF.NE.0 .OR. I.LT.ILE. OR. I.GT.ITE) GO TO 40 iplot = iplot + 1 WRITE (1,330) iplot, X(I), CPU(i), EM1U(i),CPL(i),EM1L(i) 40 CONTINUE c IF ( IEM .EQ. 0 ) GO TO 50 IF ( PHYS ) WRITE (6,300) 50 CONTINUE DO 60 J=JMIN,JMAX YM(J) = Y(J) * YFACT 60 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,310) JMIN,JMAX,(YM(J),J=JMIN,JMAX) WRITE (15,310) JMIN,JMAX,(YM(J),J=JMIN,JMAX) RETURN 70 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,320) WRITE (15,320) IF (IREF .EQ. 2) GO TO 15 ABORT1 = .TRUE. RETURN C 200 FORMAT(1H0,9X,4HCL =,F10.6/10X,4HCM =,F10.6/9X,5HCP* =,F10.6) 210 FORMAT(1H0,27X,5HLOWER,23X,5HUPPER/28X,4HY=0-,24X,4HY=0+) 220 FORMAT(3X,1HI,8X,1HX,10X,2HCP,10X,A2,14X,2HCP,10X,A2/) 230 FORMAT(25X,20HAIRFOIL LEADING EDGE,45X,20HAIRFOIL LEADING EDGE) 240 FORMAT(25X,21HAIRFOIL TRAILING EDGE,44X,21HAIRFOIL TRAILING EDGE) 250 FORMAT(1H ,I3,3F12.6,4X,2F12.6,2X,60A1) 260 FORMAT(60H1 FORCE COEFFICIENTS, PRESSURE COEFFICIENT, AND MACH NUM 1BER /2X,59H(OR SIMILARITY PARAMETER) ON BODY AND DIVIDING STREAM L 2INE.) 270 FORMAT(20X,11HCOARSE MESH) 280 FORMAT(20X,11HMEDIUM MESH) 290 FORMAT(20X,11H FINAL MESH) 300 FORMAT(20H0***** CAUTION *****/ 1 32H MAXIMUM MACH NUMBER EXCEEDS 1.3/ 2 69H SHOCK JUMPS IN ERROR IF UPSTREAM NORMAL MACH NUMBER GREATER T 3HAN 1.3) 310 FORMAT(1H0//9X,7HY(J) J=,I3,3H TO,I3/(6X,6F12.6)) 320 FORMAT(1H0,// 1 ' DETACHED SHOCK WAVE UPSTREAM OF X-MESH,SOLUTION TERMINATED.') 330 format(2x,i3,2x,f7.4,2x,f10.5,2x,f7.4,2x,f10.5,2x,f7.4) 340 format(2x,'TSFOIL2',3x,'Mach = ',f7.3,3x,'CL = ',f7.3/ 1 4x,'i',5x,'X/C',8x,'Cp-up',5x,'M-up',6x,'Cp-low',4x,'M-low') C END SUBROUTINE PRTFLD C C PRINTS PRESSURE COEFFICIENT, FLOW ANGLE AND C MACH NUMBER IN FLOW FIELD. NUMBER OF J LINES C PRINTED IS DETERMINED FROM THE INPUT VALUE OF C PRTFLO. C PRTFLO = 1 , NONE. C PRTFLO = 2 , ALL J LINES EXCEPT J0. C PRTFLO = 3 , THREE J LINES AROUND JERROR. C CALLED BY - PRINT. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON /COM13/ CDFACT , CLFACT , CMFACT , CPFACT , CPSTAR LOGICAL OUTERR COMMON /COM18/ ERROR , I1 , I2 , IERROR , JERROR , 1 OUTERR , EMU(100,2) , VC(100) , 2 WI , DCIRC , POLD(100,2) LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT COMMON /COM30/ JLIN(100), REST(304) C DIMENSION TMAC(2), CPPR(3), PYPR(3), YPRINT(3), PRT(10), EM1(3) DATA TMAC / 2HM1, 2HK1/ DATA PRT / 4HMACH,4H NUM,4HBERS,4H ,4H , 1 4HSIMI,4HLARI,4HTY P,4HARAM,4HETER/ IF (PRTFLO.EQ. 2) GO TO 25 C LOCATE LINES AROUND JERROR. JL = 3 IF (JERROR .EQ. JMIN) GO TO 10 IF (JERROR .EQ. JLOW) GO TO 20 IF (JERROR .EQ. JUP ) GO TO 10 IF (JERROR .EQ. JMAX) GO TO 20 JLIN(1) = JERROR - 1 JLIN(2) = JERROR JLIN(3) = JERROR + 1 GO TO 30 10 CONTINUE JLIN(1) = JERROR JLIN(2) = JERROR + 1 JLIN(3) = JERROR + 2 GO TO 30 20 CONTINUE JLIN(1) = JERROR - 2 JLIN(2) = JERROR - 1 JLIN(3) = JERROR GO TO 30 25 CONTINUE JL = JMAX - JMIN + 1 C FILL JLIN ARRAY WITH VALUES OF J. K = 1 DO 28 J=JMIN,JMAX JLIN(K) = J K = K + 1 28 CONTINUE 30 CONTINUE C C PRINT FLOW FIELD IN 3 J LINES PER PAGE. DO 100 MPR = 1,JL,3 MPREND = MIN0(MPR+2,JL) DO 40 M=MPR,MPREND MQ = M-MPR+1 JQ = JLIN(M) YPRINT(MQ) = Y(JQ)*YFACT 40 CONTINUE C WRITE PAGE HEADER. IS = 1 IF ( PHYS ) IS = 6 IE = IS + 4 C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,900) (PRT(I),I=IS,IE) , CPSTAR C WRITE(6, 901) (JLIN(M),M=MPR,MPREND) C WRITE(6, 902)(YPRINT(M),M=1,MQ) WRITE (15,900) (PRT(I),I=IS,IE) , CPSTAR WRITE(15, 901) (JLIN(M),M=MPR,MPREND) WRITE(15, 902)(YPRINT(M),M=1,MQ) KT = 2 IF ( PHYS ) KT = 1 C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,903) TMAC(KT), TMAC(KT), TMAC(KT) WRITE (15,903) TMAC(KT), TMAC(KT), TMAC(KT) DO 60 I=IMIN,IMAX DO 50 M=MPR,MPREND MQ = M-MPR+1 J = JLIN(M) U = PX(I,J) CPPR(MQ) = -2.0 * CPFACT * U PYPR(MQ) = VFACT*PY(I,J) EM1(MQ) = EMACH1(U) 50 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,904) I,X(I),(CPPR(M),PYPR(M),EM1(M),M=1,MQ) WRITE (15,904) I,X(I),(CPPR(M),PYPR(M),EM1(M),M=1,MQ) 60 CONTINUE 100 CONTINUE RETURN C 900 FORMAT(47H1PRESSURE COEFFICIENTS, FLOW ANGLES, AND LOCAL ,5A4/ 1 20H ON Y=CONSTANT LINES/9H CPSTAR =,F12.7//) 901 FORMAT(13X,3(15X,2HJ=,I4,15X)) 902 FORMAT(13X,3(12X,2HY=,F10.6,12X)) 903 FORMAT(4H0 I,8X,1HX,5X,3(6X,2HCP,8X,5HTHETA,7X,A2,6X)//) 904 FORMAT(1X,I3,2X,F10.6,1X,3(2X,3F11.6,1X)) C END SUBROUTINE PRTMC C SUBROUTINE TO PRINT A CHARACTER FOR EACH C POINT IN THE GRID DESCRIBING THE TYPE OF C FLOW. S , SHOCK POINT C H , HYPERBOLIC POINT C P , PARABOLIC POINT C - , ELLIPTIC POINT C CALLED BY - PRINT. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK COMMON /COM22/ CXC(100) , CXL(100), CXR(100), CXXC(100),CXXL(100), 1 CXXR(100), C1(100) C COMMON /COM30/ IPC(100) , VT(100,2), REST(104) DATA IHP/1HP/, IHH/1HH/, IHS/1HS/, IHD/1H-/, IB/1H / C C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,100) WRITE (15,100) DO 5 I=1,50,2 IPC(I) = IB IPC(I+1) = IB 5 CONTINUE DO 10 J=JMIN,JMAX VT(J,1) = C1(2) 10 CONTINUE DO 60 K=JMIN,JMAX J = JMAX - K + 1 DO 50 I=IUP,IDOWN VT(J,2) = VT(J,1) VT(J,1) = C1(I)-(CXL(I)*P(J,I-1)+CXC(I)*P(J,I)+CXR(I)*P(J,I+1)) IF (VT(J,1) .GT. 0.0) GO TO 30 IF (VT(J,2) .LT. 0.0) GO TO 20 C PARABOLIC POINT (SONIC) IPC(I) = IHP GO TO 50 20 CONTINUE C HYPERBOLIC POINT (SUPERSONIC) IPC(I) = IHH GO TO 50 30 CONTINUE IF (VT(J,2) .LT. 0.0) GO TO 40 C ELLIPTIC POINT (SUBSONIC) IPC(I) = IHD GO TO 50 40 CONTINUE C SHOCK POINT IPC(I) = IHS 50 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,110) J, (IPC(I),I=IUP,IDOWN) WRITE (15,110) J, (IPC(I),I=IUP,IDOWN) 60 CONTINUE RETURN C 100 FORMAT(39H1 FLOW AT EACH GRID POINT. P PARABOLIC/ 1 28X,12HH HYPERBOLIC/28X,7HS SHOCK/28X,10H- ELLIPTIC//) 110 FORMAT(10X,I3,5X,100A1) C END SUBROUTINE PRTSK(Z,ARG,L,NSHOCK,CDSK,LPRT1) C PRINTOUT WAVE DRAG CONTRIBUTION AND TOTAL PRESSURE C LOSS ALONG SHOCK WAVE C CALLED BY - CDCOLE. C COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK COMMON /COM13/ CDFACT , CLFACT , CMFACT , CPFACT , CPSTAR LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT DIMENSION Z(1), ARG(1) CDYCOF = -CDFACT*GAM1/(6.*YFACT) POYCOF = DELTA*DELTA*GAM1*(GAM1-1.)/12. C Changed to write to a file - Andy Ko 4/3/03 C IF(NSHOCK .EQ.1) WRITE(6,1001) C WRITE(6,1002) NSHOCK, CDSK IF(NSHOCK .EQ.1) WRITE(15,1001) WRITE(15,1002) NSHOCK, CDSK DO 10 K = 1,L YY = Z(K)*YFACT CDY = CDYCOF*ARG(K) POY = 1. + POYCOF*ARG(K) C Changed to write to a file - Andy Ko 4/3/03 C WRITE(6,1003) YY,CDY,POY WRITE(15,1003) YY,CDY,POY 10 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C IF(LPRT1 .EQ. 1) WRITE(6,1004) IF(LPRT1 .EQ. 1) WRITE(15,1004) RETURN 1001 FORMAT(74H1INVISCID WAKE PROFILES FOR INDIVIDUAL SHOCK WAVES WITHI 1N MOMENTUM CONTOUR) 1002 FORMAT(6H0SHOCK,I3/26H WAVE DRAG FOR THIS SHOCK=,F12.6/ 1 6X,1HY,9X,5HCD(Y),8X,8HPO/POINF/) 1003 FORMAT(1X,3F12.8) 1004 FORMAT(35H0SHOCK WAVE EXTENDS OUTSIDE CONTOUR/ 1 61H PRINTOUT OF SHOCK LOSSES ARE NOT AVAILABLE FOR REST OF SHOCK) END SUBROUTINE PRTWAL C C PRINTS PRESSURE COEFFICIENT AND FLOW ANGLE C ON Y=-H AND Y=+H, AND PLOTS CP ALONG SIDE OF C TABULATION. C CALLED BY - PRINT. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM5/ XDIFF(100),YDIFF(100) COMMON /COM12/ F , H , HALFPI , PI , RTKPOR , 1 TWOPI COMMON /COM13/ CDFACT , CLFACT , CMFACT , CPFACT , CPSTAR LOGICAL PHYS INTEGER PRTFLO , SIMDEF COMMON /COM27/ CL , DELTA , DELRT2 , EMACH , EMROOT , 1 PHYS , PRTFLO , SIMDEF , SONVEL , VFACT , 2 YFACT INTEGER BCTYPE COMMON /COM28/ BCTYPE , CIRCFF , FHINV , POR , CIRCTE COMMON /COM30/ CPLW(100), CPUW(100), VLW(100), VUW(100), REST(4) C DIMENSION BCT(15) , LINE1(60) C DATA IB /1H /, IL / 1HL/, IU /1HU/, IS / 1H*/, IBB /1HB/ DATA BCT /4H ,4HFREE,4H AIR,4H SO,4HLID ,4HWALL,4H , 1 4HFREE,4H JET,4HSLOT,4HTED ,4HWALL,4H POR,4HOUS ,4HWALL/ C C PRINT SINGLE VARIABLES C I2 = 3 * BCTYPE I1 = I2 - 2 C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,903) (BCT(I),I=I1,I2) WRITE (15,903) (BCT(I),I=I1,I2) THH = H * YFACT C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,904) THH WRITE (15,904) THH IF( BCTYPE .LT. 5 ) GO TO 4 PORF = POR / YFACT C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,905) PORF WRITE (15,905) PORF 4 CONTINUE IF( BCTYPE .NE. 4 .AND. BCTYPE .NE. 6 ) GO TO 6 C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,906) F WRITE (15,906) F 6 CONTINUE C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,907) CPSTAR WRITE (15,907) CPSTAR C CPMIN = 1.0E37 CPMAX = -CPMIN C CPT = -2.0 * CPFACT DO 10 I=IUP,IDOWN CPLW(I) = CPT * PX(I,JMIN) CPUW(I) = CPT * PX(I,JMAX) CPMAX = AMAX1(CPMAX,CPUW(I), CPLW(I)) CPMIN = AMIN1(CPMIN,CPUW(I), CPLW(I)) 10 CONTINUE DO 20 I=IUP,IDOWN IF (BCTYPE .NE. 2) GO TO 11 C SOLID WALL VLW(I) = 0.0 VUW(I) = 0.0 GO TO 20 11 CONTINUE IF (BCTYPE .NE. 3) GO TO 12 C FREE JET VLW(I) = VFACT * PY(I,JMIN) VUW(I) = VFACT * PY(I,JMAX) GO TO 20 12 CONTINUE IF (BCTYPE .NE. 4) GO TO 13 C SLOTTED WALL VLW(I) = VFACT * FHINV * (P(JBOT,I) + .75 * CIRCFF) VUW(I) = -VFACT * FHINV * (P(JTOP,I) - .25 * CIRCFF) GO TO 20 13 CONTINUE C POROUS WALL IF (POR .GT. 1.5) GO TO 14 VLW(I) = VFACT * POR * XDIFF(I)*(P(JMIN,I)-P(JMIN,I-1)) VUW(I) = -VFACT * POR * XDIFF(I)*(P(JMAX,I)-P(JMAX,I-1)) GO TO 20 14 CONTINUE VLW(I) = VFACT *.25*(P(JMIN+1,I+1)+2.*P(JMIN+1,I)+P(JMIN+1,I-1) 1 - P(JMIN ,I+1)-2.*P(JMIN ,I)-P(JMIN ,I-1)) 2 / (Y(JMIN+1)-Y(JMIN)) VUW(I) = VFACT *.25*(P(JMAX,I+1) +2.*P(JMAX,I) +P(JMAX,I-1) 1 - P(JMAX-1,I+1)-2.*P(JMAX-1,I)-P(JMAX-1,I-1)) 2 / (Y(JMAX)-Y(JMAX-1)) 20 CONTINUE C CPLARG = AMAX1(CPMAX,ABS(CPMIN)) UNPCOL = CPLARG / 29. C C LOCATE CP* FOR PRINTER PLOT C COL = -CPSTAR / UNPCOL NCOL = SIGN((ABS(COL)+.5), COL) NCOLS = NCOL + 30 C C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,210) WRITE (15,210) C C PRINT COLUMN HEADERS. C C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,220) WRITE (15,220) C DO 40 I=IUP,IDOWN C DO 30 K=1,60 LINE1(K) = IB 30 CONTINUE C COL = -CPUW(I) / UNPCOL NCOL = SIGN((ABS(COL) + .5) , COL) NCOLU = NCOL + 30 LINE1(NCOLU) = IU C COL = -CPLW(I) / UNPCOL NCOL = SIGN((ABS(COL) + .5) , COL) NCOLL = NCOL + 30 LINE1(NCOLL) = IL IF ( NCOLL .EQ. NCOLU ) LINE1(NCOLL) = IBB IF (IABS(NCOLS) .LT. 61) LINE1(NCOLS) = IS C C Changed to write to a file - Andy Ko 4/3/03 C WRITE (6,250) I, X(I), CPLW(I), VLW(I), CPUW(I), VUW(I), LINE1 WRITE (15,250) I, X(I), CPLW(I), VLW(I), CPUW(I), VUW(I), LINE1 40 CONTINUE RETURN C 210 FORMAT(1H0,27X,5HLOWER,23X,5HUPPER/28X,4HY=-H,24X,4HY=+H) 220 FORMAT(3X,1HI,8X,1HX,10X,2HCP,9X,5HTHETA,12X,2HCP,9X,5HTHETA/) 250 FORMAT(1H ,I3,3F12.6,4X,2F12.6,2X,60A1) 903 FORMAT(2H1 ,3A4,20H BOUNDARY CONDITION.) 904 FORMAT(1H0,10X,24HH (TUNNEL HALF HEIGHT) =,F9.6) 905 FORMAT(1H0,11X,23HPOR (POROSITY FACTOR) =,F9.6) 906 FORMAT(1H0,14X,20HF (SLOT PARAMETER) =,F9.6) 907 FORMAT(1H0,29X,5HCP* =,F9.6) C END FUNCTION PX(I,J) C FUNCTION PX COMPUTES U = DP/DX AT POINT I,J C CALLED BY - CDCOLE, DRAG, FINDSK, MACHMP, M1LINE, C NEWSK, PRINT1, PRTFLD, PRTWAL. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM5/ XDIFF(100),YDIFF(100) C C TEST TO LOCATE END POINTS IF(I .EQ. IMIN) GO TO 10 IF(I .EQ. IMAX) GO TO 20 C INTERIOR MESH POINT PJI = P(J,I) PX = .5*(XDIFF(I+1)*(P(J,I+1)-PJI) + XDIFF(I)*(PJI-P(J,I-1))) RETURN 10 CONTINUE C UPSTREAM BOUNDARY PX = 1.5*XDIFF(I+1)*(P(J,I+1)-P(J,I)) - 1 0.5*XDIFF(I+2)*(P(J,I+2) - P(J,I+1)) RETURN 20 CONTINUE C DOWNSTREAM BOUNDARY PX = 1.5*XDIFF(I)*(P(J,I) - P(J,I-1)) 1 -0.5*XDIFF(I-1)*(P(J,I-1) - P(J,I-2)) RETURN END FUNCTION PY(I,J) C FUNCTION PY COMPUTES V = DP/DY AT POINT I,J C CALLED BY - CDCOLE, PRTFLD, PRTWAL. C COMMON P(102,101),X(100) , Y(100) COMMON / COM1/ IMIN , IMAX , IUP , IDOWN , ILE , 1 ITE , JMIN , JMAX , JUP , JLOW , 2 JTOP , JBOT , J1 , J2 COMMON / COM2/ AK , ALPHA , DUB , GAM1 , RTK COMMON / COM5/ XDIFF(100),YDIFF(100) COMMON / COM6/ FL(100) , FXL(100), FU(100) , FXU(100), 1 CAMBER(100), THICK(100),VOL , XFOIL(100), IFOIL COMMON /COM26/ PJUMP(100) INTEGER BCTYPE COMMON /COM28/ BCTYPE , CIRCFF , FHINV , POR , CIRCTE C C TEST FOR END POINTS OR POINTS NEAR AIRFOIL SLIT IF(J .EQ. JMIN) GO TO 10 IF(J .EQ. JLOW) GO TO 20 IF(J .EQ. JUP ) GO TO 40 IF(J .EQ. JMAX) GO TO 50 C I,J IS AN INTERIOR POINT PJI = P(J,I) PY = .5*(YDIFF(J+1)*(P(J+1,I)-PJI) + YDIFF(J)*(PJI-P(J-1,I))) RETURN 10 CONTINUE C I,J IS ON LOWER BOUNDARY. USE ONE SIDED DERIVATIVE PY = 1.5* YDIFF(J+1)*(P(J+1,I) - P(J,I)) - 1 0.5* YDIFF(J+2)*(P(J+2,I) - P(J+1,I)) RETURN 20 CONTINUE C I,J IS ON ROW OF MESH PIINTS BELOW AIRFOIL. C VMINUS = YDIFF(J)*(P(J,I) - P(J-1,I)) C C TEST TO SEE IF I,J IS AHEAD, UNDER, OR BEHIND SLIT. IF(I .LT. ILE) GO TO 21 IF(I .GT. ITE) GO TO 22 C I,J IS UNDER AIRFOIL. USE DERIVATIVE BOUNDARY C CONDITION. IC = I - ILE + 1 PY = .5 * (FXL(IC) - ALPHA + VMINUS) RETURN 21 CONTINUE C I,J IS AHEAD OF AIRFOIL. PY = .5*((P(JUP,I) - P(JLOW,I)) * YDIFF(JUP) + VMINUS) RETURN 22 CONTINUE C I,J IS BEHIND AIRFOIL PY = .5*((P(JUP,I) - PJUMP(I) - P(JLOW,I)) * YDIFF(JUP) + VMINUS) RETURN 40 CONTINUE C I,J IS ON ROW OF MESH POINTS ABOVE AIRFOIL VPLUS = YDIFF(J+1)*(P(J+1,I) - P(J,I)) C TEST TO SEE IF I IS AHEAD OF, OVER, OR BEHIND C AIRFOIL SLIT. IF(I .LT. ILE) GO TO 41 IF(I .GT. ITE) GO TO 42 IC = I - ILE + 1 PY = .5 * (VPLUS + FXU(IC) - ALPHA) RETURN 41 CONTINUE C I,J IS AHEAD OF AIRFOIL. PY = .5*((P(JUP,I) - P(JLOW,I)) * YDIFF(JUP) + VPLUS) RETURN 42 CONTINUE C I,J IS BEHIND AIRFOIL PY = .5*((P(JUP,I) - PJUMP(I) - P(JLOW,I)) * YDIFF(JUP) + VPLUS) RETURN 50 CONTINUE C I,J IS ON TOP ROW OF MESH POINTS. USE ONE SIDED C FO