aboutsummaryrefslogtreecommitdiff
path: root/src/xoper.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/xoper.f')
-rw-r--r--src/xoper.f2780
1 files changed, 2780 insertions, 0 deletions
diff --git a/src/xoper.f b/src/xoper.f
new file mode 100644
index 0000000..3ffcfdb
--- /dev/null
+++ b/src/xoper.f
@@ -0,0 +1,2780 @@
+C***********************************************************************
+C Module: xoper.f
+C
+C Copyright (C) 2000 Mark Drela
+C
+C This program is free software; you can redistribute it and/or modify
+C it under the terms of the GNU General Public License as published by
+C the Free Software Foundation; either version 2 of the License, or
+C (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program; if not, write to the Free Software
+C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+C***********************************************************************
+C
+ SUBROUTINE OPER
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*1 ANS
+ CHARACTER*4 COMAND, COMOLD
+ LOGICAL LRECALC, LCPX
+C
+ CHARACTER*128 COMARG, ARGOLD, LINE
+C
+ PARAMETER (NPRX = 101)
+ DIMENSION XPR(NPRX), YPR(NPRX), FPR(NPRX)
+C
+ DIMENSION NBLP(NPX)
+ DIMENSION IPPAI(NPX), NAPOLT(NPX)
+C
+ DIMENSION IINPUT(20)
+ DIMENSION RINPUT(20)
+ LOGICAL ERROR
+C
+C---- retain last-command info if OPER is exited and then re-entered
+ SAVE COMOLD, ARGOLD
+C
+C---- logical units for polar save file, polar dump file
+ LUPLR = 9
+ LUPLX = 11
+C
+ COMAND = '****'
+ COMARG = ' '
+ LRECALC = .FALSE.
+ LCPX = .FALSE.
+C
+ IF(N.EQ.0) THEN
+ WRITE(*,*)
+ WRITE(*,*) '*** No airfoil available ***'
+ RETURN
+ ENDIF
+C
+ IF(IPACT.NE.0) THEN
+ WRITE(*,5000) IPACT
+ 5000 FORMAT(/' Polar', I3,' is active')
+ ENDIF
+C
+ccc 500 CONTINUE
+ COMOLD = COMAND
+ ARGOLD = COMARG
+C
+C====================================================
+C---- start of menu loop
+ 500 CONTINUE
+C
+ IF(LVISC) THEN
+ IF(LPACC) THEN
+ CALL ASKC('.OPERva^',COMAND,COMARG)
+ ELSE
+ CALL ASKC('.OPERv^',COMAND,COMARG)
+ ENDIF
+ ELSE
+ IF(LPACC) THEN
+ CALL ASKC('.OPERia^',COMAND,COMARG)
+ ELSE
+ CALL ASKC('.OPERi^',COMAND,COMARG)
+ ENDIF
+ ENDIF
+C
+C---- process previous command ?
+ IF(COMAND(1:1).EQ.'!') THEN
+ IF(COMOLD.EQ.'****') THEN
+ WRITE(*,*) 'Previous .OPER command not valid'
+ GO TO 500
+ ELSE
+ COMAND = COMOLD
+ COMARG = ARGOLD
+ LRECALC = .TRUE.
+ ENDIF
+ ELSE
+ LRECALC = .FALSE.
+ ENDIF
+C
+ IF(COMAND.EQ.' ') THEN
+C----- just <return> was typed... clean up plotting and exit OPER
+ IF(LPLOT) CALL PLEND
+ LPLOT = .FALSE.
+ CALL CLRZOOM
+ RETURN
+ ENDIF
+C
+C---- extract command line numeric arguments
+ DO I=1, 20
+ IINPUT(I) = 0
+ RINPUT(I) = 0.0
+ ENDDO
+ NINPUT = 20
+ CALL GETFLT(COMARG,RINPUT,NINPUT,ERROR)
+C
+C---- don't try to read integers, since might get integer overflow
+ DO I=1, NINPUT
+ IF(ABS(RINPUT(I)) .GT. 2.1E9) THEN
+ IINPUT(I) = 2**31
+ ELSE
+ IINPUT(I) = INT(RINPUT(I))
+ ENDIF
+ ENDDO
+C
+ccc NINPUT = 20
+ccc CALL GETINT(COMARG,IINPUT,NINPUT,ERROR)
+C
+C--------------------------------------------------------
+ IF(COMAND.EQ.'? ') THEN
+ WRITE(*,1050)
+ 1050 FORMAT(
+ & /' <cr> Return to Top Level'
+ & /' ! Redo last ALFA,CLI,CL,ASEQ,CSEQ,VELS'
+ &//' Visc r Toggle Inviscid/Viscous mode'
+ & /' .VPAR Change BL parameter(s)'
+ & /' Re r Change Reynolds number'
+ & /' Mach r Change Mach number'
+ & /' Type i Change type of Mach,Re variation with CL'
+ & /' ITER Change viscous-solution iteration limit'
+ & /' INIT Toggle BL initialization flag'
+ &//' Alfa r Prescribe alpha'
+ & /' CLI r Prescribe inviscid CL'
+ & /' Cl r Prescribe CL'
+ & /' ASeq rrr Prescribe a sequence of alphas'
+ & /' CSeq rrr Prescribe a sequence of CLs'
+ &//' SEQP Toggle polar/Cp(x) sequence plot display'
+ & /' CINC Toggle minimum Cp inclusion in polar'
+ & /' HINC Toggle hinge moment inclusion in polar'
+ & /' Pacc i Toggle auto point accumulation to active polar'
+ & /' PGET f Read new polar from save file'
+ & /' PWRT i Write polar to save file'
+ & /' PSUM Show summary of stored polars'
+ & /' PLIS i List stored polar(s)'
+ & /' PDEL i Delete stored polar'
+ & /' PSOR i Sort stored polar'
+ & /' PPlo ii. Plot stored polar(s)'
+ & /' APlo ii. Plot stored airfoil(s) for each polar'
+ & /' ASET i Copy stored airfoil into current airfoil'
+ & /' PREM ir. Remove point(s) from stored polar'
+ & /' PNAM i Change airfoil name of stored polar'
+ & /' PPAX Change polar plot axis limits'
+ &//' RGET f Read new reference polar from file'
+ & /' RDEL i Delete stored reference polar'
+ &//' GRID Toggle Cp vs x grid overlay'
+ & /' CREF Toggle reference Cp data overlay'
+ & /' FREF Toggle reference CL,CD.. data display'
+ &//' CPx Plot Cp vs x'
+ & /' CPV Plot airfoil with pressure vectors (gee wiz)'
+ & /' .VPlo BL variable plots'
+ & /' .ANNO Annotate current plot'
+ & /' HARD Hardcopy current plot'
+ & /' SIZE r Change plot-object size'
+ & /' CPMI r Change minimum Cp axis annotation'
+ &//' BL i Plot boundary layer velocity profiles'
+ & /' BLC Plot boundary layer velocity profiles at cursor'
+ & /' BLWT r Change velocity profile scale weight'
+ &//' FMOM Calculate flap hinge moment and forces'
+ & /' FNEW rr Set new flap hinge point'
+ & /' VELS rr Calculate velocity components at a point'
+ & /' DUMP f Output Ue,Dstar,Theta,Cf vs s,x,y to file'
+ & /' CPWR f Output x vs Cp to file'
+ & /' CPMN Report minimum surface Cp'
+ & /' NAME s Specify new airfoil name'
+ & /' NINC Increment name version number')
+c &//' IMAG Toggle image-airfoil'
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'Z ') THEN
+ CALL USETZOOM(.TRUE.,.TRUE.)
+ CALL REPLOT(IDEV)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'U ') THEN
+ CALL CLRZOOM
+ CALL REPLOT(IDEV)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'VISC' .OR.
+ & COMAND.EQ.'V ' ) THEN
+ IF(LPACC) THEN
+ WRITE(*,2100)
+ GO TO 500
+ ENDIF
+C
+ LVISC = .NOT. LVISC
+C
+ IF(LVISC) THEN
+ IF(NINPUT.GE.1) THEN
+ REINF1 = RINPUT(1)
+ ELSE IF(REINF1 .EQ. 0.0) THEN
+ CALL ASKR('Enter Reynolds number^',REINF1)
+ ENDIF
+C
+ CALL MRSHOW(.TRUE.,.TRUE.)
+ ENDIF
+ LVCONV = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'HARD') THEN
+ IF(LPLOT) CALL PLEND
+ LPLOT = .FALSE.
+ CALL REPLOT(IDEVRP)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'SIZE') THEN
+ IF(NINPUT.GE.1) THEN
+ SIZE = RINPUT(1)
+ ELSE
+ WRITE(*,*) 'Current plot-object size =', SIZE
+ CALL ASKR('Enter new plot-object size^',SIZE)
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CPMI') THEN
+ IF(NINPUT.GE.1) THEN
+ CPMIN = RINPUT(1)
+ ELSE
+ WRITE(*,*) 'Current CPmin =', CPMIN
+ CALL ASKR('Enter new CPmin^',CPMIN)
+ ENDIF
+C
+ PFAC = PLOTAR/(CPMAX-CPMIN)
+ CPDEL = -0.5
+ IF(CPMIN .LT. -4.01) CPDEL = -1.0
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'UEMA') THEN
+ IF(NINPUT.GE.1) THEN
+ UEMAX = RINPUT(1)
+ ELSE
+ WRITE(*,*) 'Current UEmax =', UEMAX
+ CALL ASKR('Enter new UEMAX^',UEMAX)
+ ENDIF
+C
+ UFAC = PLOTAR/(UEMAX-UEMIN)
+ UEDEL = 0.2
+ IF((UEMAX-UEMIN) .GT. 2.51) UEDEL = 0.5
+ IF((UEMAX-UEMIN) .GT. 5.01) UEDEL = 1.0
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'UEMI') THEN
+ IF(NINPUT.GE.1) THEN
+ UEMIN = RINPUT(1)
+ ELSE
+ WRITE(*,*) 'Current UEmin =', UEMIN
+ CALL ASKR('Enter new UEMIN^',UEMIN)
+ ENDIF
+C
+ UFAC = PLOTAR/(UEMAX-UEMIN)
+ UEDEL = 0.2
+ IF((UEMAX-UEMIN) .GT. 2.51) UEDEL = 0.5
+ IF((UEMAX-UEMIN) .GT. 5.01) UEDEL = 1.0
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'VPAR') THEN
+ CALL VPAR
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'RE ' .OR.
+ & COMAND.EQ.'R ' ) THEN
+ IF(LPACC .AND. LVISC) THEN
+ WRITE(*,2100)
+ GO TO 500
+ ENDIF
+C
+ IF(NINPUT.GE.1) THEN
+ REINF1 = RINPUT(1)
+ ELSE
+ WRITE(*,*)
+ WRITE(*,*) 'Currently...'
+ CALL MRSHOW(.FALSE.,.TRUE.)
+ CALL ASKR('Enter new Reynolds number^',REINF1)
+ ENDIF
+C
+ccc CALL MRSHOW(.FALSE.,.TRUE.)
+ CALL MRCL(1.0,MINF_CL,REINF_CL)
+ LVCONV = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'MACH' .OR.
+ & COMAND.EQ.'M ' ) THEN
+ IF(LPACC) THEN
+ WRITE(*,2100)
+ GO TO 500
+ ENDIF
+C
+ 15 CONTINUE
+ IF(NINPUT.GE.1) THEN
+ MINF1 = RINPUT(1)
+ ELSE
+ WRITE(*,*)
+ WRITE(*,*) 'Currently...'
+ CALL MRSHOW(.TRUE.,.FALSE.)
+ CALL ASKR('Enter Mach number^',MINF1)
+ ENDIF
+C
+ IF(MINF1.GE.1.0) THEN
+ WRITE(*,*) 'Supersonic freestream not allowed'
+ NINPUT = 0
+ GO TO 15
+ ENDIF
+ccc CALL MRSHOW(.TRUE.,.FALSE.)
+ CALL MRCL(1.0,MINF_CL,REINF_CL)
+ CALL COMSET
+C
+ IF(MINF.GT.0.0) WRITE(*,1300) CPSTAR, QSTAR/QINF
+ 1300 FORMAT(/' Sonic Cp =', F10.2, ' Sonic Q/Qinf =', F10.3/)
+C
+ CALL CPCALC(N,QINV,QINF,MINF,CPI)
+ IF(LVISC) CALL CPCALC(N+NW,QVIS,QINF,MINF,CPV)
+ CALL CLCALC(N,X,Y,GAM,GAM_A,ALFA,MINF,QINF, XCMREF,YCMREF,
+ & CL,CM,CDP, CL_ALF,CL_MSQ)
+ CALL CDCALC
+ LVCONV = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'TYPE' .OR.
+ & COMAND.EQ.'T' ) THEN
+ IF(LPACC) THEN
+ WRITE(*,2100)
+ GO TO 500
+ ENDIF
+C
+ 17 CONTINUE
+ IF(NINPUT.GE.1) THEN
+ ITYP = IINPUT(1)
+ ELSE
+ WRITE(*,1105)
+ 1105 FORMAT(
+ & /' Type parameters held constant varying fixed '
+ & /' ---- ------------------------ ------- -----------'
+ & /' 1 M , Re .. lift chord, vel.'
+ & /' 2 M sqrt(CL) , Re sqrt(CL) .. vel. chord, lift'
+ & /' 3 M , Re CL .. chord lift , vel.')
+ CALL ASKI('Enter type of Mach,Re variation with CL^',ITYP)
+ ENDIF
+C
+ IF(ITYP.EQ.1) THEN
+ MATYP = 1
+ RETYP = 1
+ ELSE IF(ITYP.EQ.2) THEN
+ MATYP = 2
+ RETYP = 2
+ ELSE IF(ITYP.EQ.3) THEN
+ MATYP = 1
+ RETYP = 3
+ ENDIF
+C
+ IF(ITYP.LT.1 .OR. ITYP.GT.3) THEN
+ NINPUT = 0
+ GO TO 17
+ ENDIF
+C
+ CALL MRSHOW(.TRUE.,.TRUE.)
+ CALL MRCL(1.0,MINF_CL,REINF_CL)
+ CALL COMSET
+ LVCONV = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'ITER') THEN
+ 18 CONTINUE
+ IF(NINPUT.GE.1) THEN
+ ITMAX = IINPUT(1)
+ ELSE
+ WRITE(*,*) 'Current iteration limit:', ITMAX
+ CALL ASKI('Enter new iteration limit^',ITMAX)
+ ENDIF
+C
+ IF(ITMAX.LT.1) THEN
+ NINPUT = 0
+ GO TO 18
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'INIT') THEN
+ LBLINI = .NOT.LBLINI
+ IF(LBLINI) THEN
+ WRITE(*,*) 'BLs are assumed to be initialized'
+ ELSE
+ WRITE(*,*) 'BLs will be initialized on next point'
+ LIPAN = .FALSE.
+ ENDIF
+C
+C--------------------------------------------------------
+c ELSEIF(COMAND.EQ.'IMAG') THEN
+c LIMAGE = .NOT.LIMAGE
+c IF(LIMAGE) THEN
+c CALL ASKR('Enter y-position of image plane^',YIMAGE)
+c CALL ASKI('Specify image type (1=wall -1=free jet)^',KIMAGE)
+c ELSE
+c WRITE(*,*) 'Image airfoil removed'
+c ENDIF
+c LGAMU = .FALSE.
+c LQAIJ = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'ALFA' .OR.
+ & COMAND.EQ.'A ' ) THEN
+ IF(.NOT.LRECALC) THEN
+C------- set inviscid solution only if point is not being recalculated
+ IF(NINPUT.GE.1) THEN
+ ADEG = RINPUT(1)
+ ELSE
+ ADEG = ALFA/DTOR
+ CALL ASKR('Enter angle of attack (deg)^',ADEG)
+ ENDIF
+ LALFA = .TRUE.
+ ALFA = DTOR*ADEG
+ QINF = 1.0
+ CALL SPECAL
+ IF(ABS(ALFA-AWAKE) .GT. 1.0E-5) LWAKE = .FALSE.
+ IF(ABS(ALFA-AVISC) .GT. 1.0E-5) LVCONV = .FALSE.
+ IF(ABS(MINF-MVISC) .GT. 1.0E-5) LVCONV = .FALSE.
+ ENDIF
+C
+ IF(LVISC) CALL VISCAL(ITMAX)
+ CALL CPX
+ CALL FCPMIN
+C
+ccc IF( LVISC .AND. LPACC .AND. LVCONV ) THEN
+ IF( LPACC .AND. (LVCONV .OR. .NOT.LVISC)) THEN
+ CALL PLRADD(LUPLR,IPACT)
+ CALL PLXADD(LUPLX,IPACT)
+ ENDIF
+C
+ IF(LVISC .AND. .NOT.LPACC .AND. .NOT.LVCONV) THEN
+ WRITE(*,*) 'Type "!" to continue iterating'
+ ENDIF
+C
+C WRITE(*,*) 'N NW =', N, NW
+C call aski('Enter i^',ioff)
+C call askr('Enter dmass^',dms)
+Cc
+C do 43 is=1, 2
+C do 430 ibl=2, nbl(is)
+C i = ipan(ibl,is)
+C mass(ibl,is) = 0.
+C if(i.eq.ioff) mass(ibl,is) = dms
+C 430 continue
+C 43 continue
+Cc
+C call ueset
+C call qvfue
+C call gamqv
+C call cpcalc(N+NW,QVIS,QINF,MINF,CPV)
+C call cdcalc
+c CALL CLCALC(N,X,Y,GAM,GAM_A,ALFA,MINF,QINF,XCMREF,YCMREF,
+c & CL,CM,CDP, CL_ALF,CL_MSQ)
+C call cpx
+Cc
+C
+ COMOLD = COMAND
+ ARGOLD = COMARG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CLI ') THEN
+ IF(.NOT.LRECALC) THEN
+ IF(NINPUT.GE.1) THEN
+ CLSPEC = RINPUT(1)
+ ELSE
+ CALL ASKR('Enter inviscid lift coefficient^',CLSPEC)
+ ENDIF
+ LALFA = .TRUE.
+ ALFA = 0.0
+ QINF = 1.0
+ CALL SPECCL
+ ADEG = ALFA/DTOR
+ IF(ABS(ALFA-AWAKE) .GT. 1.0E-5) LWAKE = .FALSE.
+ IF(ABS(ALFA-AVISC) .GT. 1.0E-5) LVCONV = .FALSE.
+ IF(ABS(MINF-MVISC) .GT. 1.0E-5) LVCONV = .FALSE.
+ ENDIF
+C
+ IF(LVISC) CALL VISCAL(ITMAX)
+ CALL CPX
+ CALL FCPMIN
+C
+ccc IF( LVISC .AND. LPACC .AND. LVCONV ) THEN
+ IF( LPACC .AND. (LVCONV .OR. .NOT.LVISC)) THEN
+ CALL PLRADD(LUPLR,IPACT)
+ CALL PLXADD(LUPLX,IPACT)
+ ENDIF
+C
+ IF(LVISC .AND. .NOT.LPACC .AND. .NOT.LVCONV) THEN
+ WRITE(*,*) 'Type "!" to continue iterating'
+ ENDIF
+C
+ COMOLD = COMAND
+ ARGOLD = COMARG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CL ' .OR.
+ & COMAND.EQ.'C ' ) THEN
+ IF(.NOT.LRECALC) THEN
+ IF(NINPUT.GE.1) THEN
+ CLSPEC = RINPUT(1)
+ ELSE
+ CALL ASKR('Enter lift coefficient^',CLSPEC)
+ ENDIF
+ LALFA = .FALSE.
+ ALFA = 0.0
+ QINF = 1.0
+ CALL SPECCL
+ ADEG = ALFA/DTOR
+ IF(ABS(ALFA-AWAKE) .GT. 1.0E-5) LWAKE = .FALSE.
+ IF(ABS(ALFA-AVISC) .GT. 1.0E-5) LVCONV = .FALSE.
+ IF(ABS(MINF-MVISC) .GT. 1.0E-5) LVCONV = .FALSE.
+ ENDIF
+ IF(LVISC) CALL VISCAL(ITMAX)
+ CALL FCPMIN
+C
+ CALL CPX
+ccc IF( LVISC .AND. LPACC .AND. LVCONV ) THEN
+ IF( LPACC .AND. (LVCONV .OR. .NOT.LVISC)) THEN
+ CALL PLRADD(LUPLR,IPACT)
+ CALL PLXADD(LUPLX,IPACT)
+ ENDIF
+C
+ IF(LVISC .AND. .NOT.LPACC .AND. .NOT.LVCONV) THEN
+ WRITE(*,*) 'Type "!" to continue iterating'
+ ENDIF
+C
+ COMOLD = COMAND
+ ARGOLD = COMARG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'ASEQ' .OR.
+ & COMAND.EQ.'AS ' .OR.
+ & COMAND.EQ.'CSEQ' .OR.
+ & COMAND.EQ.'CS ' ) THEN
+ LALFA = COMAND.EQ.'ASEQ' .OR.
+ & COMAND.EQ.'AS '
+C
+ IF(LALFA) THEN
+ IF (NINPUT.GE.3) THEN
+ AA1 = RINPUT(1)
+ AA2 = RINPUT(2)
+ DAA = RINPUT(3)
+ ELSEIF(NINPUT.GE.2) THEN
+ AA1 = RINPUT(1)
+ AA2 = RINPUT(2)
+ DAA = DAA/DTOR
+ CALL ASKR('Enter alfa increment (deg)^',DAA)
+ ELSEIF(NINPUT.GE.1) THEN
+ AA1 = RINPUT(1)
+ AA2 = AA2/DTOR
+ CALL ASKR('Enter last alfa value (deg)^',AA2)
+ DAA = DAA/DTOR
+ CALL ASKR('Enter alfa increment (deg)^',DAA)
+ ELSE
+ AA1 = AA1/DTOR
+ CALL ASKR('Enter first alfa value (deg)^',AA1)
+ AA2 = AA2/DTOR
+ CALL ASKR('Enter last alfa value (deg)^',AA2)
+ DAA = DAA/DTOR
+ CALL ASKR('Enter alfa increment (deg)^',DAA)
+ ENDIF
+ IF(AA2.LT.AA1) THEN
+ DAA = -ABS(DAA)
+ ELSE
+ DAA = ABS(DAA)
+ ENDIF
+ AA1 = AA1*DTOR
+ AA2 = AA2*DTOR
+ DAA = DAA*DTOR
+ NPOINT = 1
+ IF(DAA .NE. 0.0) NPOINT = INT((AA2-AA1)/DAA + 0.5) + 1
+C
+ ELSE
+ IF (NINPUT.GE.3) THEN
+ CL1 = RINPUT(1)
+ CL2 = RINPUT(2)
+ DCL = RINPUT(3)
+ ELSEIF(NINPUT.GE.2) THEN
+ CL1 = RINPUT(1)
+ CL2 = RINPUT(2)
+ CALL ASKR('Enter CL increment ^',DCL)
+ ELSEIF(NINPUT.GE.1) THEN
+ CL1 = RINPUT(1)
+ CALL ASKR('Enter last CL value^',CL2)
+ CALL ASKR('Enter CL increment ^',DCL)
+ ELSE
+ CALL ASKR('Enter first CL value^',CL1)
+ CALL ASKR('Enter last CL value^',CL2)
+ CALL ASKR('Enter CL increment ^',DCL)
+ ENDIF
+ IF(CL2.LT.CL1) THEN
+ DCL = -ABS(DCL)
+ ELSE
+ DCL = ABS(DCL)
+ ENDIF
+ NPOINT = 1
+ IF(DCL .NE. 0.0) NPOINT = INT((CL2-CL1)/DCL + 0.5) + 1
+ ENDIF
+C
+C- - - - - - - - - - - - - - - - - -
+C
+C----- initialize plot
+ CALL PLTINI
+C
+ IF(LPPSHO) THEN
+C------ set up for polar plot
+C
+ ELSE
+C------ set up for Cp(x) plot
+C
+C------ Cp scaling factor
+ PFAC = PLOTAR/(CPMAX-CPMIN)
+C
+C------ determine airfoil box size and location
+ CALL AIRLIM(N,X,Y,XMIN,XMAX,YMIN,YMAX)
+C
+C------ y-offset for airfoil in Cp vs x plot
+ FACA = FACAIR/(XMAX-XMIN)
+ XOFA = XOFAIR*(XMAX-XMIN) - XMIN
+ YOFA = YOFAIR*(XMAX-XMIN) - YMAX - CPMAX*PFAC/FACA
+C
+C------ re-origin for Cp vs x plot
+ CALL PLOT(0.09 , 0.04 + CPMAX*PFAC + (YMAX-YMIN)*FACA, -3)
+C
+C------ draw axes and airfoil picture for Cp vs x plot
+ CALL CPAXES(LCPGRD,
+ & N,X,Y,XOFA,YOFA,FACA,
+ & CPMIN,CPMAX,CPDEL,PFAC,CH,
+ & 'XFOIL',VERSION)
+C
+C------ set initial x,y-positions of sequence plot label top
+ XL = 0.65
+ IF(LVISC) XL = 0.48
+ YL = -CPMIN*PFAC
+C
+C------ draw sequence plot label
+ CALL SEQLAB(XL,YL,XL1,XL2,XL3,XL4,XL5,XL6,CHSEQ,1,LVISC)
+C
+ CALL PLFLUSH
+C
+C------ set label y position
+ YL = YL - 0.2*CH
+ ENDIF
+C
+C----- initialize unconverged-point counter
+ ISEQEX = 0
+ ALAST = ADEG
+ CLAST = CL
+C
+C----- calculate each point, add Cp distribution to plot, and save to polar
+ DO 115 IPOINT=1, NPOINT
+C
+C------- set proper alpha for this point
+ IF(LALFA) THEN
+ ALFA = AA1 + DAA*FLOAT(IPOINT-1)
+ ELSE
+ CLSPEC = CL1 + DCL*FLOAT(IPOINT-1)
+ CALL SPECCL
+ ENDIF
+C
+ IF(ABS(ALFA-AWAKE) .GT. 1.0E-5) LWAKE = .FALSE.
+ IF(ABS(ALFA-AVISC) .GT. 1.0E-5) LVCONV = .FALSE.
+ IF(ABS(MINF-MVISC) .GT. 1.0E-5) LVCONV = .FALSE.
+ CALL SPECAL
+ ITMAXS = ITMAX + 5
+ IF(LVISC) CALL VISCAL(ITMAXS)
+C
+ ADEG = ALFA/DTOR
+C
+ CALL FCPMIN
+C
+C------- add point to buffer polar and/or disk files
+ccc IF( LVISC .AND. LPACC .AND. LVCONV ) THEN
+ IF( LPACC .AND. (LVCONV .OR. .NOT.LVISC)) THEN
+ CALL PLRADD(LUPLR,IPACT)
+ CALL PLXADD(LUPLX,IPACT)
+ ENDIF
+C
+ IF(LPPSHO) THEN
+ CALL PLTINI
+ccc CALL PLOTABS(0.5,0.5,-3)
+ PSIZE = 1.0*SIZE
+ CALL NEWFACTOR(PSIZE)
+ CALL PLOT(5.0*CH,7.0*CH,-3)
+C
+ CH1 = CH*0.90
+ CH2 = CH*0.75
+ CLEXP = 1.0
+C
+ DO IP=1, NPOL
+ NBLP(IP) = 1
+ ENDDO
+C
+ CALL POLPLT(NAX,NPOL,NAPOL,CPOL,
+ & REYNP1,MACHP1,ACRITP,PTRATP,ETAPP,
+ & NAMEPOL,ICOLP,ILINP,
+ & NFX,NPOLREF,NDREF,CPOLREF,NAMEREF,ICOLR,ISYMR,
+ & ISX,NBLP,CPOLSD ,IMATYP,IRETYP,
+ & ' ','XFOIL',VERSION,
+ & PLOTAR,XCDWID,XALWID,XOCWID,CH1,CH2,CLEXP,
+ & LPGRID,LPCDW,LPLIST,LPLEGN,LAECEN,LPCDH,LPCMDOT,
+ & CPOLPLF,' ',0)
+ ELSE
+C-------- add alpha, CL, etc. to plot
+ CALL SEQPLT(YL,XL1,XL2,XL3,XL4,XL5,XL6,CHSEQ,ADEG,CL,CM,LVISC)
+C
+C-------- add sonic Cp dashed line if within plot
+ IF(CPSTAR.GE.CPMIN) CALL DASH(0.0,XL-CH,-CPSTAR*PFAC)
+C
+ CALL NEWPEN(2)
+ IF(LVISC) THEN
+C--------- Plot viscous -Cp distribution on airfoil
+ CALL XYLINE(N+NW,X,CPV,-XOFA,FACA,0.0,-PFAC,1)
+ ELSE
+C--------- Plot inviscid -Cp distribution on airfoil
+ CALL XYLINE(N,X,CPI,-XOFA,FACA,0.0,-PFAC,1)
+ ENDIF
+ ENDIF
+C
+ CALL PLFLUSH
+c###
+ccc call dcpout
+C
+ IF(LVISC .AND. .NOT.LVCONV) THEN
+C-------- increment unconverged-point counter
+ ISEQEX = ISEQEX + 1
+ IF(ISEQEX .GE. NSEQEX) THEN
+ WRITE(*,1150) ISEQEX, ALAST, CLAST
+ 1150 FORMAT(
+ & /' Sequence halted since previous',I3,' points did not converge'
+ & /' Last-converged alpha =', F8.3, ' CL =', F10.5)
+ GO TO 116
+ ENDIF
+ ELSE
+C-------- converged OK... reset unconverged-point counter
+ ISEQEX = 0
+ ALAST = ADEG
+ CLAST = CL
+ ENDIF
+C
+ 115 CONTINUE
+ 116 CONTINUE
+ccc CALL ASKC('hit <cr>^',DUMMY,COMARG)
+C
+ COMOLD = COMAND
+ ARGOLD = COMARG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'SEQP') THEN
+ LPPSHO = .NOT.LPPSHO
+ IF(LPPSHO) THEN
+ WRITE(*,*) 'Polar will be plotted during point sequence'
+ ELSE
+ WRITE(*,*) 'Cp(x) will be plotted during point sequence'
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'PACC' .OR.
+ & COMAND.EQ.'P ' ) THEN
+ LPACC = .NOT.LPACC
+C
+ IF(LPACC) THEN
+ IF(NINPUT.GE.1) THEN
+C------- slot into which accumulated polar will go
+ IP = MIN( MAX( IINPUT(1) , 0 ) , NPOL+1 )
+ ELSE
+C------- no command argument was given... just use next available slot
+ IP = NPOL+1
+ PFNAME(IP) = ' '
+ PFNAMX(IP) = ' '
+ ENDIF
+C
+ IF(IP.GT.NPOL) THEN
+ IF(NPOL.EQ.NPX) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'Number of polars is at array limit'
+ WRITE(*,*) 'New polar will not be stored'
+ IPACT = 0
+ ELSE
+ IPACT = NPOL + 1
+ PFNAME(IPACT) = ' '
+ PFNAMX(IPACT) = ' '
+ ENDIF
+C
+ ELSE
+ IPACT = IP
+C
+ ENDIF
+C
+C------ set up for appending to new or existing polar (if IPACT > 0)
+ CALL PLRSET(IPACT)
+C
+C------ jump out if decision was made to abort polar accumulation
+ IF(IPACT.LE.0) THEN
+ LPACC = .FALSE.
+ GO TO 500
+ ENDIF
+C
+ CALL PLRINI(LUPLR,IPACT)
+ CALL PLXINI(LUPLX,IPACT)
+ WRITE(*,*)
+ WRITE(*,*) 'Polar accumulation enabled'
+C
+ ELSE
+ WRITE(*,*)
+ WRITE(*,*) 'Polar accumulation disabled'
+ IPACT = 0
+C
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'PGET') THEN
+ IF(NPOL.GE.NPX) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'Number of polars is at array limit'
+ WRITE(*,*) 'Delete with PDEL if necessary'
+ GO TO 500
+ ENDIF
+C
+ IP = NPOL+1
+C
+ IF(COMARG.EQ.' ') THEN
+ CALL ASKS('Enter polar filename^',FNAME)
+ ELSE
+ FNAME = COMARG
+ ENDIF
+C
+ LU = 17
+ CALL POLREAD(LU,FNAME,ERROR,
+ & NAX,NAPOL(IP),CPOL(1,1,IP),
+ & REYNP1(IP),MACHP1(IP),ACRITP(IP),XSTRIPP(1,IP),
+ & PTRATP(IP),ETAPP(IP),
+ & NAMEPOL(IP),IRETYP(IP),IMATYP(IP),
+ & ISX,NBLP(IP),CPOLSD(1,1,1,IP),
+ & CODEPOL(IP),VERSPOL(IP) )
+ IF(ERROR) THEN
+ WRITE(*,*) 'Polar file READ error'
+ ELSE
+ NPOL = IP
+ NXYPOL(IP) = 0
+ CALL STRIP(NAMEPOL(IP),NNAMEP)
+ NEL = 1
+ CALL POLWRIT(6,' ',ERROR, .TRUE.,
+ & NAX, 1,NAPOL(IP), CPOL(1,1,IP), IPOL,NIPOL,
+ & REYNP1(IP),MACHP1(IP),ACRITP(IP),XSTRIPP(1,IP),
+ & PTRATP(IP),ETAPP(IP),
+ & NAMEPOL(IP),IRETYP(IP),IMATYP(IP),
+ & ISX,NEL,CPOLSD(1,1,1,IP), JPOL,NJPOL,
+ & CODEPOL(IP),VERSPOL(IP), .FALSE. )
+ PFNAME(IP) = FNAME
+ WRITE(*,5500) IP
+ 5500 FORMAT(/' Stored as Polar', I4)
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'PWRT') THEN
+ 75 CONTINUE
+ IF(NPOL.EQ.1) THEN
+ IP = 1
+ ELSEIF(NINPUT.EQ.0) THEN
+ CALL PLRSUM(1,NPOL,IPACT)
+ CALL ASKI(
+ & 'Enter index of polar to write (0=all, -1=abort)^',IP)
+ IF(IP.EQ.-1) GO TO 500
+ ELSE
+ IP = IINPUT(1)
+ ENDIF
+C
+ IF(IP.EQ.0) THEN
+ IP1 = 1
+ IP2 = NPOL
+ ELSEIF(IP.GE.1 .AND. IP.LE.NPOL) THEN
+ IP1 = IP
+ IP2 = IP
+ ELSE
+ NINPUT = 0
+ GO TO 75
+ ENDIF
+C
+ NEL = 1
+ DO IP = IP1, IP2
+ LU = 19
+ CALL PLRSUM(IP,IP,IPACT)
+ CALL STRIP(PFNAME(IP),NPF)
+ IF(NPF.EQ.0) THEN
+ LINE = 'Enter polar output filename^'
+ ELSE
+ LINE = 'Enter polar output filename ['
+ & // PFNAME(IP)(1:NPF) // ']^'
+ ENDIF
+ CALL ASKS(LINE,FNAME)
+ IF(NPF.NE.0 .AND. FNAME.EQ.' ') FNAME = PFNAME(IP)
+C
+ NIPOL = NIPOL0
+ IF(LCMINP) THEN
+ NIPOL = NIPOL + 1
+ IPOL(IMC) = NIPOL
+ ENDIF
+ IF(LHMOMP) THEN
+ NIPOL = NIPOL + 1
+ IPOL(ICH) = NIPOL
+ ENDIF
+C
+ CALL POLWRIT(LU,FNAME,ERROR, .TRUE.,
+ & NAX, 1,NAPOL(IP),CPOL(1,1,IP), IPOL,NIPOL,
+ & REYNP1(IP),MACHP1(IP),ACRITP(IP),XSTRIPP(1,IP),
+ & PTRATP(IP),ETAPP(IP),
+ & NAMEPOL(IP),IRETYP(IP),IMATYP(IP),
+ & ISX,NEL,CPOLSD(1,1,1,IP), JPOL,NJPOL,
+ & 'XFOIL',VERSION, .TRUE. )
+ IF(ERROR) THEN
+ WRITE(*,1075) IP
+ 1075 FORMAT(' Polar', I3,' not written')
+ ELSE
+ PFNAME(IP) = FNAME
+ ENDIF
+ ENDDO
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'RGET') THEN
+ IF(NPOLREF.GE.NPX) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'Number of reference polars is at array limit'
+ WRITE(*,*) 'Delete with RDEL if necessary'
+ GO TO 500
+ ENDIF
+C
+ IR = NPOLREF+1
+C
+ IF(COMARG.EQ.' ') THEN
+ CALL ASKS('Enter reference polar filename^',FNAME)
+ ELSE
+ FNAME = COMARG
+ ENDIF
+C
+ LU = 9
+ OPEN(LU,FILE=FNAME,STATUS='OLD',ERR=27)
+ CALL POLREF(LU, FNAME, ERROR,
+ & NFX, NDREF(1,IR), CPOLREF(1,1,1,IR), NAMEREF(IR))
+ CLOSE(LU)
+ IF(ERROR) GO TO 27
+C
+ NPOLREF = IR
+C
+ CALL STRIP(NAMEREF(IR),NNREF)
+ IF(NNREF.EQ.0) THEN
+ CALL ASKS('Enter label for reference polar^',NAMEREF(IR))
+ CALL STRIP(NAMEREF(IR),NNREF)
+ ELSE
+ WRITE(*,*)
+ WRITE(*,*) NAMEREF(IR)
+ ENDIF
+C
+ccc ICOLR(IR) = NCOLOR - IR + 1
+ ICOLR(IR) = 2 + IR
+ ISYMR(IR) = MOD(IR,10)
+ 25 CONTINUE
+ GO TO 500
+C
+ 27 CONTINUE
+ WRITE(*,*) 'File OPEN error'
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'RDEL') THEN
+ IF(NPOLREF.EQ.0) THEN
+ WRITE(*,*) 'No reference polars are stored'
+ GO TO 500
+ ENDIF
+C
+ IF(NINPUT.GE.1) THEN
+ IR = IINPUT(1)
+ ELSE
+ IR = NPOLREF+1
+ ENDIF
+C
+ 35 CONTINUE
+C
+ IF(IR.EQ.0) THEN
+C------- delete all polars
+ NPOLREF = 0
+C
+ ELSEIF(IR.EQ.-1) THEN
+C------- abort
+ GO TO 500
+C
+ ELSEIF(IR.LT.-1 .OR. IR.GT.NPOLREF) THEN
+ CALL PRFSUM(1,NPOLREF)
+ CALL ASKI(
+ & 'Specify ref. polar to delete (0 = all, -1 = abort)^',IR)
+ GO TO 35
+C
+ ELSE
+C------- delete ref. polar IR
+ DO JR = IR+1, NPOLREF
+ CALL PRFCOP(JR,JR-1)
+ WRITE(*,1310) JR, JR-1
+ 1410 FORMAT(' Ref.polar',I3,' moved into ref.polar',I3)
+ ENDDO
+ NPOLREF = NPOLREF-1
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'PSUM') THEN
+ IF(NPOL.EQ.0) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'No polars are stored'
+ GO TO 500
+ ENDIF
+C
+ CALL PLRSUM(1,NPOL,IPACT)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'PLIS') THEN
+ IF(NPOL.EQ.0) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'No polars are stored'
+ GO TO 500
+ ENDIF
+C
+ IF(NINPUT.EQ.0) THEN
+ IP1 = 1
+ IP2 = NPOL
+ ELSE
+ IP = IINPUT(1)
+ IF(IP.EQ.0) THEN
+ IP1 = 1
+ IP2 = NPOL
+ ELSEIF(IP.GE.1 .AND. IP.LE.NPOL) THEN
+ IP1 = IP
+ IP2 = IP
+ ELSE
+ WRITE(*,*)
+ WRITE(*,*) 'Specified stored polar does not exist'
+ GO TO 500
+ ENDIF
+ ENDIF
+C
+ NIPOL = NIPOL0
+ IF(LCMINP) THEN
+ NIPOL = NIPOL + 1
+ IPOL(IMC) = NIPOL
+ ENDIF
+ IF(LHMOMP) THEN
+ NIPOL = NIPOL + 1
+ IPOL(ICH) = NIPOL
+ ENDIF
+C
+ NEL = 1
+ DO IP = IP1, IP2
+ WRITE(*,3100) IP
+ 3100 FORMAT(
+ &/' =============================================================='
+ &/' Polar', I3)
+ IA1 = 1
+ IA2 = NAPOL(IP)
+ CALL POLWRIT(6,' ',ERROR, .TRUE.,
+ & NAX, IA1,IA2, CPOL(1,1,IP), IPOL,NIPOL,
+ & REYNP1(IP),MACHP1(IP),ACRITP(IP),XSTRIPP(1,IP),
+ & PTRATP(IP),ETAPP(IP),
+ & NAMEPOL(IP), IRETYP(IP),IMATYP(IP),
+ & ISX,NEL,CPOLSD(1,1,1,IP), JPOL,NJPOL,
+ & 'XFOIL',VERSION, .FALSE.)
+ ENDDO
+ NIPOL = NIPOL0
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'PDEL') THEN
+ IF(NPOL.EQ.0) THEN
+ WRITE(*,*) 'No polars are stored'
+ GO TO 500
+ ENDIF
+C
+ IF(NINPUT.GE.1) THEN
+C------- use command argument
+ IP = IINPUT(1)
+ ELSE
+C------- no argument given... set up for user query test below
+ IP = NPOL+1
+ ENDIF
+C
+ 40 CONTINUE
+ IF(IP.EQ.0) THEN
+C------- delete all polars
+ NPOL = 0
+ IPACT = 0
+ LPACC = .FALSE.
+C
+ ELSEIF(IP.EQ.-1) THEN
+C------- abort
+ GO TO 500
+C
+ ELSEIF(IP.LT.-1 .OR. IP.GT.NPOL) THEN
+ CALL PLRSUM(1,NPOL,IPACT)
+ CALL ASKI(
+ & 'Specify polar to delete (0 = all, -1 = abort)^',IP)
+ GO TO 40
+C
+ ELSE
+C------- delete polar IP
+ IF(IPACT.EQ.IP) THEN
+ WRITE(*,*) 'Active polar deleted. Accumulation turned off'
+ IPACT = 0
+ LPACC = .FALSE.
+ ENDIF
+C
+ DO JP = IP+1, NPOL
+ CALL PLRCOP(JP,JP-1)
+ WRITE(*,1310) JP, JP-1
+ 1310 FORMAT(' Polar',I3,' moved into polar',I3)
+ IF(IPACT.EQ.JP) THEN
+ IPACT = JP-1
+ ENDIF
+ ENDDO
+ NPOL = NPOL-1
+C
+ ENDIF
+C
+ IF(IPACT.GT.0) THEN
+ WRITE(*,1320) IPACT
+ 1320 FORMAT(' Polar',I3,' is now active')
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'PSOR') THEN
+ IF(NPOL.EQ.0) THEN
+ WRITE(*,*) 'No polars are stored'
+ GO TO 500
+ ENDIF
+C
+ IF(NINPUT.GE.1) THEN
+C------- use command argument
+ IP = IINPUT(1)
+ ELSE
+C------- no argument given... set up for user query test below
+ IP = NPOL+1
+ ENDIF
+C
+C------ sort polars in increasing alpha
+ IDSORT = IAL
+C
+ 42 CONTINUE
+ IF (IP.EQ.-1) THEN
+C------- abort
+ GO TO 500
+C
+ ELSEIF(IP.LT.-1 .OR. IP.GT.NPOL) THEN
+ CALL PLRSUM(1,NPOL,IPACT)
+ CALL ASKI(
+ & 'Specify polar to sort (0 = all, -1 = abort)^',IP)
+ GO TO 42
+C
+ ELSE
+C------- sort polar(s)
+ IF(IP.EQ.0) THEN
+ IP1 = 1
+ IP2 = NPOL
+ ELSE
+ IP1 = IP
+ IP2 = IP
+ ENDIF
+ DO JP = IP1, IP2
+ CALL PLRSRT(JP,IDSORT)
+ ENDDO
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'PPLO' .OR.
+ & COMAND.EQ.'PP ' ) THEN
+C------ set temporary polar-size array to plot only selected polars
+ IF(NINPUT.EQ.0) THEN
+C------- no polars specified... plot all of them
+ DO IP=1, NPOL
+ NAPOLT(IP) = NAPOL(IP)
+ ENDDO
+ ELSE
+C------- set up to plot only specified polars
+ DO IP=1, NPOL
+ NAPOLT(IP) = 0
+ ENDDO
+ DO K=1, NINPUT
+ IP = IINPUT(K)
+ IF(IP.GE.1 .AND. IP.LE.NPOL) NAPOLT(IP) = NAPOL(IP)
+ ENDDO
+ ENDIF
+C
+ CALL PLTINI
+ccc CALL PLOTABS(0.5,0.5,-3)
+ PSIZE = 1.0*SIZE
+ CALL NEWFACTOR(PSIZE)
+ CALL PLOT(5.0*CH,7.0*CH,-3)
+C
+ CH1 = CH*0.90
+ CH2 = CH*0.75
+ CLEXP = 1.0
+ DO IP=1, NPOL
+ NBLP(IP) = 1
+ ENDDO
+C
+ CALL POLPLT(NAX,NPOL,NAPOLT,CPOL,
+ & REYNP1,MACHP1,ACRITP,PTRATP,ETAPP,
+ & NAMEPOL,ICOLP,ILINP,
+ & NFX,NPOLREF,NDREF,CPOLREF,NAMEREF,ICOLR,ISYMR,
+ & ISX,NBLP,CPOLSD ,IMATYP,IRETYP,
+ & ' ','XFOIL',VERSION,
+ & PLOTAR,XCDWID,XALWID,XOCWID,CH1,CH2,CLEXP,
+ & LPGRID,LPCDW,LPLIST,LPLEGN,LAECEN,LPCDH,LPCMDOT,
+ & CPOLPLF,' ',0)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'APLO' .OR.
+ & COMAND.EQ.'AP ' ) THEN
+ IF(NPOL.EQ.0) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'No polars are stored'
+ GO TO 500
+ ENDIF
+C
+ IF(NINPUT.EQ.0) THEN
+ NPPAI = NPOL
+ DO K=1, NPPAI
+ IPPAI(K) = K
+ ENDDO
+ ELSE
+ NPPAI = MIN( NINPUT , NPX )
+ DO K=1, NPPAI
+ IINP = IINPUT(K)
+ IF(IINP.GE.1 .AND. IINP.LE.NPOL) THEN
+ IPPAI(K) = IINP
+ ELSE
+ IPPAI(K) = 0
+ ENDIF
+ ENDDO
+ ENDIF
+C
+ CALL PPAPLT(NPPAI,IPPAI)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'ASET') THEN
+ IF(NPOL.EQ.0) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'No polar airfoils are stored'
+ GO TO 500
+ ENDIF
+C
+ 50 CONTINUE
+ IF(NINPUT.EQ.0) THEN
+ IF(NPOL.EQ.1) THEN
+ IP = 1
+ ELSE
+ CALL PLRSUM(1,NPOL,IPACT)
+ CALL ASKI('Enter index of polar airfoil to set^',IP)
+ ENDIF
+ ELSE
+ IP = IINPUT(1)
+ ENDIF
+C
+ IF(IP.EQ.0) THEN
+ GO TO 500
+ ELSEIF(IP.LT.1 .OR. IP.GT.NPOL) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'Specified polar airfoil does not exist'
+ NINPUT = 0
+ GO TO 50
+ ENDIF
+C
+ WRITE(*,*)
+ WRITE(*,*) 'Current airfoil will be overwritten. Proceed? Y'
+ READ(*,1000) ANS
+ 1000 FORMAT(A)
+C
+ IF(INDEX('Nn',ANS) .NE. 0) THEN
+ WRITE(*,*) 'No action taken'
+ GO TO 500
+ ELSE
+ CALL APCOPY(IP)
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'PREM') THEN
+ IF(NPOL.EQ.0) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'No polars are stored'
+ GO TO 500
+ ENDIF
+C
+ 52 CONTINUE
+ IF(NINPUT.EQ.0) THEN
+ IF(NPOL.EQ.1) THEN
+ IP = 1
+ ELSE
+ CALL PLRSUM(1,NPOL,IPACT)
+ CALL ASKI('Enter index of polar to modify^',IP)
+ ENDIF
+ ELSE
+ IP = IINPUT(1)
+ ENDIF
+C
+ IF(IP.EQ.0) THEN
+ GO TO 500
+ ELSEIF(IP.LT.1 .OR. IP.GT.NPOL) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'Specified polar airfoil does not exist'
+ NINPUT = 0
+ GO TO 52
+ ENDIF
+C
+ IF(NINPUT.GE.2) THEN
+ NREM = NINPUT - 1
+ ELSE
+ NIPOL = NIPOL0
+ IF(LCMINP) THEN
+ NIPOL = NIPOL + 1
+ IPOL(IMC) = NIPOL
+ ENDIF
+ IF(LHMOMP) THEN
+ NIPOL = NIPOL + 1
+ IPOL(ICH) = NIPOL
+ ENDIF
+C
+ WRITE(*,3100) IP
+ IA1 = 1
+ IA2 = NAPOL(IP)
+ CALL POLWRIT(6,' ',ERROR, .TRUE.,
+ & NAX, IA1,IA2, CPOL(1,1,IP), IPOL,NIPOL,
+ & REYNP1(IP),MACHP1(IP),ACRITP(IP),XSTRIPP(1,IP),
+ & PTRATP(IP),ETAPP(IP),
+ & NAMEPOL(IP), IRETYP(IP),IMATYP(IP),
+ & ISX,1,CPOLSD(1,1,1,IP), JPOL,NJPOL,
+ & 'XFOIL',VERSION, .FALSE. )
+ 53 WRITE(*,3220)
+ 3220 FORMAT(/' Enter alpha(s) of points to be removed: ', $)
+ READ(*,1000) LINE
+ NREM = 19
+ CALL GETFLT(LINE,RINPUT(2),NREM,ERROR)
+ IF(ERROR) GO TO 53
+ ENDIF
+C
+C----- go over specified alphas to be removed
+ DO 55 IREM = 1, NREM
+C------- check all alpha points in polar IP
+ DO IA = 1, NAPOL(IP)
+ ADIF = CPOL(IA,IAL,IP) - RINPUT(IREM+1)
+ IF(ABS(ADIF) .LT. 0.0005) THEN
+C---------- alphas match within 3-digit print tolerance...
+C- remove point by pulling down all points above it
+ DO JA = IA, NAPOL(IP)-1
+ DO K = 1, IPTOT
+ CPOL(JA,K,IP) = CPOL(JA+1,K,IP)
+ ENDDO
+ DO K = 1, JPTOT
+ CPOLSD(JA,1,K,IP) = CPOLSD(JA+1,1,K,IP)
+ CPOLSD(JA,2,K,IP) = CPOLSD(JA+1,2,K,IP)
+ ENDDO
+ ENDDO
+C---------- shrink polar by 1
+ NAPOL(IP) = NAPOL(IP) - 1
+C
+ IF(NAPOL(IP).LE.0) THEN
+C----------- last point has been removed... eliminate this polar IP
+ DO JP = IP+1, NPOL
+ CALL PLRCOP(JP,JP-1)
+ IF(IPACT.EQ.JP) IPACT = JP-1
+ WRITE(*,1310) JP, JP-1
+ ENDDO
+ NPOL = NPOL-1
+C
+ IF(IPACT.GT.0) THEN
+ WRITE(*,1320) IPACT
+ ENDIF
+C
+ GO TO 500
+ ENDIF
+C
+C---------- go to next specified alpha to be removed
+ GO TO 55
+ ENDIF
+ ENDDO
+ 55 CONTINUE
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'PNAM') THEN
+ IF(NPOL.EQ.0) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'No polars are stored'
+ GO TO 500
+ ENDIF
+C
+ 58 CONTINUE
+ IF(NINPUT.EQ.0) THEN
+ IF(NPOL.EQ.1) THEN
+ IP = 1
+ ELSE
+ CALL PLRSUM(1,NPOL,IPACT)
+ CALL ASKI('Enter index of polar to modify^',IP)
+ ENDIF
+ ELSE
+ IP = IINPUT(1)
+ ENDIF
+C
+ IF(IP.EQ.0) THEN
+ GO TO 500
+ ELSEIF(IP.LT.1 .OR. IP.GT.NPOL) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'Specified polar airfoil does not exist'
+ NINPUT = 0
+ GO TO 58
+ ENDIF
+C
+ NIPOL = NIPOL0
+ IF(LCMINP) THEN
+ NIPOL = NIPOL + 1
+ IPOL(IMC) = NIPOL
+ ENDIF
+ IF(LHMOMP) THEN
+ NIPOL = NIPOL + 1
+ IPOL(ICH) = NIPOL
+ ENDIF
+C
+ WRITE(*,3100) IP
+ IA1 = 0
+ IA2 = -1
+ CALL POLWRIT(6,' ',ERROR, .TRUE.,
+ & NAX, IA1,IA2, CPOL(1,1,IP), IPOL,NIPOL,
+ & REYNP1(IP),MACHP1(IP),ACRITP(IP),XSTRIPP(1,IP),
+ & PTRATP(IP),ETAPP(IP),
+ & NAMEPOL(IP), IRETYP(IP),IMATYP(IP),
+ & ISX,1,CPOLSD(1,1,1,IP), JPOL,NJPOL,
+ & 'XFOIL',VERSION, .FALSE. )
+ NIPOL = NIPOL0
+ WRITE(*,3320)
+ 3320 FORMAT(/' Enter new airfoil name of polar: ', $)
+ READ(*,1000) NAMEPOL(IP)
+ CALL STRIP(NAMEPOL(IP),NNP)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'PPAX') THEN
+ CALL POLAXI(CPOLPLF,XCDWID,XALWI,XOCWID)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CREF') THEN
+ LCPREF = .NOT. LCPREF
+ IF(LCPREF) THEN
+ WRITE(*,*) 'Reference Cp plotting enabled'
+ ELSE
+ WRITE(*,*) 'Reference Cp plotting disabled'
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'FREF') THEN
+ LFOREF = .NOT. LFOREF
+ IF(LFOREF) THEN
+ WRITE(*,*) 'Reference force plotting enabled'
+ ELSE
+ WRITE(*,*) 'Reference force plotting disabled'
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CPX ' .OR.
+ & COMAND.EQ.'CP ' ) THEN
+ CALL CPX
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'UEX ' .OR.
+ & COMAND.EQ.'UE ' ) THEN
+ CALL UEX
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'GRID') THEN
+ LCPGRD = .NOT.LCPGRD
+ IF(LCPGRD) THEN
+ WRITE(*,*) 'Cp grid overlay enabled'
+ ELSE
+ WRITE(*,*) 'Cp grid overlay disabled'
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CPV ') THEN
+ CALL CPVEC
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'BL ') THEN
+ IF(.NOT.LVCONV) THEN
+ WRITE(*,*) 'Compute valid viscous solution first'
+ GO TO 500
+ ENDIF
+C
+ IF(NINPUT.GE.1) THEN
+ NPR = MIN( IINPUT(1) , NPRX )
+ ELSE
+ NPR = 21
+ WRITE(*,*) 'Using default number of profiles:', NPR
+ ENDIF
+C
+ IF(NPR.GT.1) THEN
+C------ set NPR points along surface, offset slightly for the locating logic
+ DOFF = 0.00001*(S(N)-S(1))
+ DO IPR = 1, NPR
+ FRAC = FLOAT(IPR-1)/FLOAT(NPR-1)
+ SPR = S(1) + (S(N)-S(1))*FRAC
+ XPR(IPR) = SEVAL(SPR,X,XP,S,N) + DOFF*DEVAL(SPR,Y,YP,S,N)
+ YPR(IPR) = SEVAL(SPR,Y,YP,S,N) - DOFF*DEVAL(SPR,X,XP,S,N)
+ ENDDO
+ ENDIF
+C
+ CALL CPX
+ CALL DPLOT(NPR,XPR,YPR)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'BLC ') THEN
+ IF(.NOT.LVCONV) THEN
+ WRITE(*,*) 'Compute valid viscous solution first'
+ GO TO 500
+ ENDIF
+C
+ NPR = 0
+ CALL CPX
+ CALL DPLOT(NPR,XPR,YPR)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'BLF ') THEN
+C
+ NPR = 2
+ DO IPR = 1, NPR
+C
+c WRITE(*,'(1X,A,$)') 'Enter x/c_BL, delta : '
+c READ(*,1000) LINE
+c NINP = 2
+c CALL GETFLT(LINE,RINPUT,NINP,ERROR)
+c IF(ERROR .OR. NINP.EQ.0) THEN
+c GO TO 500
+c ELSE
+c SGN = SIGN(1.0,RINPUT(1))
+c XOC = ABS(RINPUT(1))
+c DPR = RINPUT(2)
+c ENDIF
+
+ if (ipr.eq.1) then
+ xoc = 0.4
+ sgn = 1.0
+ elseif(ipr.eq.2) then
+ xoc = 0.4
+ sgn = -1.0
+ endif
+
+ if(ninput .gt. 0) then
+ dpr = rinput(1)
+ else
+ dpr = 0.01
+ endif
+
+C
+ IF(SGN .GT. 0.0) THEN
+ SPR = SLE + (S(1)-SLE)*XOC
+ ELSE
+ SPR = SLE + (S(N)-SLE)*XOC
+ ENDIF
+C
+ XPRI = XLE + (XTE-XLE)*XOC
+ CALL SINVRT(SPR,XPRI,X,XP,S,N)
+C
+ DOFF = 0.00001*(S(N)-S(1))
+ XPR(IPR) = SEVAL(SPR,X,XP,S,N) + DOFF*DEVAL(SPR,Y,YP,S,N)
+ YPR(IPR) = SEVAL(SPR,Y,YP,S,N) - DOFF*DEVAL(SPR,X,XP,S,N)
+C
+ CALL FBLGET(XPR(IPR),YPR(IPR), DPR,FPR(IPR) )
+C
+ enddo
+
+ WRITE(*,*)
+ DO IPR = 1, NPR
+ WRITE(*,7720) 'xBL, Fint =', XPR(IPR), FPR(IPR)*1.0E4
+ ENDDO
+ 7720 FORMAT(1X,A,F7.3,F12.6)
+
+ccc GO TO 770
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'BLWT') THEN
+ IF(NINPUT.GE.1) THEN
+ UPRWT = RINPUT(1)
+ ELSE
+ WRITE(*,*) 'Current u/Qinf profile plot weight =', UPRWT
+ CALL ASKR('Enter new plot weight^',UPRWT)
+ ENDIF
+C
+ CALL CPX
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'FMOM') THEN
+ CALL MHINGE
+ WRITE(*,1500) XOF,YOF,HMOM,HFX,HFY
+ 1500 FORMAT(/' Flap hinge x,y :', 2F8.4/
+ & ' 2 2'/
+ & ' Hinge moment/span = ',F8.6,' x 1/2 rho V c '/
+ & ' 2 '/
+ & ' x-Force /span = ',F8.6,' x 1/2 rho V c '/
+ & ' 2 '/
+ & ' y-Force /span = ',F8.6,' x 1/2 rho V c '/)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'FNEW') THEN
+ IF (NINPUT.GE.2) THEN
+ XOF = RINPUT(1)
+ YOF = RINPUT(2)
+ ELSEIF(NINPUT.GE.1) THEN
+ XOF = RINPUT(1)
+ YOF = -999.0
+ ELSE
+ XOF = -999.0
+ YOF = -999.0
+ ENDIF
+ LFLAP = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'VELS') THEN
+ IF (NINPUT.GE.2) THEN
+ XXX = RINPUT(1)
+ YYY = RINPUT(2)
+ ELSEIF(NINPUT.GE.1) THEN
+ XXX = RINPUT(1)
+ CALL ASKR('Enter y^',YYY)
+ ELSE
+ CALL ASKR('Enter x^',XXX)
+ CALL ASKR('Enter y^',YYY)
+ ENDIF
+ CALL PSILIN(0,XXX,YYY,-1.0,0.0,PSI,VVV,.FALSE.,.TRUE.)
+ CALL PSILIN(0,XXX,YYY, 0.0,1.0,PSI,UUU,.FALSE.,.TRUE.)
+ QQQ = SQRT(UUU**2 + VVV**2)
+ CPP = 1.0 - (UUU**2 + VVV**2)
+ WRITE(*,1800) UUU,VVV,QQQ,CPP
+ 1800 FORMAT(/' u/Uinf = ', F8.4, ' v/Uinf = ', F8.4
+ & /' q/Uinf = ', F8.4, ' Cp = ', F8.4 / )
+C
+ COMOLD = COMAND
+ ARGOLD = COMARG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'DUMP') THEN
+ CALL BLDUMP(COMARG)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CPWR') THEN
+ CALL CPDUMP(COMARG)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CPMN') THEN
+ IF(LVISC)THEN
+ WRITE(*,1769) CPMNI, XCPMNI, CPMNV, XCPMNV
+ 1769 FORMAT(' Minimum Inviscid Cp =',F8.4,' at x =',F8.4
+ & / ' Minimum Viscous Cp =',F8.4,' at x =',F8.4 )
+ ELSE
+ WRITE(*,1779) CPMNI, XCPMNI
+ 1779 FORMAT(' Minimum Inviscid Cp =',F8.4,' at x =',F8.4)
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CINC') THEN
+ LCMINP = .NOT.LCMINP
+ IF(LCMINP) THEN
+ WRITE(*,*) 'Min Cp will be written to polar save file'
+ ELSE
+ WRITE(*,*) 'Min Cp won''t be written to polar save file'
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'HINC') THEN
+ LHMOMP = .NOT.LHMOMP
+ IF(LHMOMP) THEN
+ WRITE(*,*) 'Hinge moment will be written to polar save file'
+ IF(.NOT.LFLAP) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'Note: Flap hinge location not defined'
+ WRITE(*,*) ' Set it with FNEW command'
+ ENDIF
+ ELSE
+ WRITE(*,*) 'Hinge moment won''t be written to polar save file'
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'ANNO') THEN
+ IF(LPLOT) THEN
+ CALL ANNOT(CH)
+ ELSE
+ WRITE(*,*) 'No active plot to annotate'
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'VPLO' .OR.
+ & COMAND.EQ.'VP ' ) THEN
+ CALL BLPLOT
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'NAME') THEN
+ IF(COMARG.EQ.' ') THEN
+ CALL NAMMOD(NAME,0,-1)
+ ELSE
+ NAME = COMARG
+ ENDIF
+ CALL STRIP(NAME,NNAME)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'NINC') THEN
+ CALL NAMMOD(NAME,1,1)
+ CALL STRIP(NAME,NNAME)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'NDEC') THEN
+ CALL NAMMOD(NAME,-1,1)
+ CALL STRIP(NAME,NNAME)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'DAMP') THEN
+ IF(IDAMP.EQ.0) THEN
+ IDAMP = 1
+ WRITE(*,*) 'Modified amplification used'
+ ELSE
+ IDAMP = 0
+ WRITE(*,*) 'Original amplification used'
+ ENDIF
+C--------------------------------------------------------
+ ELSE
+ WRITE(*,8000) COMAND
+ 8000 FORMAT(1X,A4,' command not recognized. Type a "?" for list')
+
+ ENDIF
+C
+C---- go back to top of menu loop
+ GO TO 500
+C
+C--------------------------------------------
+ 2100 FORMAT(/' * Polar is being accumulated.'
+ & /' * Cannot change its parameters in midstream.')
+ END ! OPER
+
+
+ SUBROUTINE FCPMIN
+C------------------------------------------------
+C Finds minimum Cp on dist for cavitation work
+C------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+C
+ XCPMNI = X(1)
+ XCPMNV = X(1)
+ CPMNI = CPI(1)
+ CPMNV = CPV(1)
+C
+ DO I = 2, N + NW
+ IF(CPI(I) .LT. CPMNI) THEN
+ XCPMNI = X(I)
+ CPMNI = CPI(I)
+ ENDIF
+ IF(CPV(I) .LT. CPMNV) THEN
+ XCPMNV = X(I)
+ CPMNV = CPV(I)
+ ENDIF
+ ENDDO
+C
+
+ IF(LVISC)THEN
+ CPMN = CPMNV
+ ELSE
+ CPMN = CPMNI
+C
+ CPMNV = CPMNI
+ XCPMNV = XCPMNI
+ ENDIF
+C
+ RETURN
+ END ! FCPMIN
+
+
+
+ SUBROUTINE MRSHOW(LM,LR)
+ INCLUDE 'XFOIL.INC'
+ LOGICAL LM, LR
+C
+ IF(LM .OR. LR) WRITE(*,*)
+C
+ IF(LM) THEN
+ IF(MATYP.EQ.1) WRITE(*,1100) MINF1
+ IF(MATYP.EQ.2) WRITE(*,1100) MINF1, ' / sqrt(CL)'
+ IF(MATYP.EQ.3) WRITE(*,1100) MINF1, ' / CL'
+ ENDIF
+C
+ IF(LR) THEN
+ IF(RETYP.EQ.1) WRITE(*,1200) REINF1
+ IF(RETYP.EQ.2) WRITE(*,1200) REINF1, ' / sqrt(CL)'
+ IF(RETYP.EQ.3) WRITE(*,1200) REINF1, ' / CL'
+ ENDIF
+C
+ RETURN
+C
+ 1100 FORMAT(1X,'M =' , F10.4, A)
+ 1200 FORMAT(1X,'Re =' , G12.4, A)
+ END ! MRSHOW
+
+
+
+ SUBROUTINE NAMMOD(NAME,KDEL,KMOD0)
+ CHARACTER*(*) NAME
+C-------------------------------------------
+C Requests new modified NAME with
+C version number in brackets, e.g.
+C NACA 0012 [5]
+C
+C If bracketed index exists in NAME,
+C it is incremented by KDEL.
+C If no bracketed index exists, it
+C is added with initial value KMOD0,
+C unless KMOD0 is negative in which
+C case nothing is added.
+C-------------------------------------------
+ CHARACTER*48 NAMDEF
+C
+ CALL STRIP(NAME,NNAME)
+ KBRACK1 = INDEX(NAME,'[')
+ KBRACK2 = INDEX(NAME,']')
+C
+ NAMDEF = NAME(1:NNAME)
+C
+ IF(KBRACK1.NE.0 .AND.
+ & KBRACK2.NE.0 .AND. KBRACK2-KBRACK1.GT.1) THEN
+C----- brackets exist... get number, (go get user's input on READ error)
+ READ(NAME(KBRACK1+1:KBRACK2-1),*,ERR=40) KMOD
+ KMOD = IABS(KMOD)
+ KMODP = MOD( KMOD+KDEL , 100 )
+ IF(KBRACK1.GE.2) THEN
+ NAME = NAME(1:KBRACK1-1)
+ ELSE
+ NAME = ' '
+ ENDIF
+ CALL STRIP(NAME,NNAME)
+ ELSEIF(KMOD0.GT.0) THEN
+ KMODP = MOD( KMOD0 , 100 )
+ ELSE
+ KMODP = 0
+ ENDIF
+C
+ IF (KMODP.GE.10) THEN
+ NAMDEF = NAME(1:NNAME) // ' [ ]'
+ WRITE(NAMDEF(NNAME+3:NNAME+4),1020) KMODP
+ 1020 FORMAT(I2)
+ ELSEIF(KMODP.GE. 1) THEN
+ NAMDEF = NAME(1:NNAME) // ' [ ]'
+ WRITE(NAMDEF(NNAME+3:NNAME+3),1025) KMODP
+ 1025 FORMAT(I1)
+ ENDIF
+C
+ 40 WRITE(*,1040) NAMDEF
+ 1040 FORMAT(/' Enter airfoil name or <return> for default: ',A)
+ READ(*,1000) NAME
+ 1000 FORMAT(A)
+ IF(NAME .EQ. ' ') NAME = NAMDEF
+C
+ RETURN
+ END ! NAMMOD
+
+
+
+ SUBROUTINE BLDUMP(FNAME1)
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*(*) FNAME1
+C
+ CHARACTER*80 FILDEF
+C
+ CHARACTER*1 DELIM
+ CHARACTER*256 LINE
+C
+ IF (KDELIM.EQ.0) THEN
+ DELIM = ' '
+ ELSEIF(KDELIM.EQ.1) THEN
+ DELIM = ','
+ ELSEIF(KDELIM.EQ.2) THEN
+ DELIM = CHAR(9)
+ ELSE
+ WRITE(*,*) '? Illegal delimiter. Using blank.'
+ DELIM = ' '
+ ENDIF
+C
+ 1000 FORMAT(50A)
+C
+ IF(FNAME1(1:1).NE.' ') THEN
+ FNAME = FNAME1
+ ELSE
+C----- no argument... get it somehow
+ IF(NPREFIX.GT.0) THEN
+C------ offer default using existing prefix
+ FILDEF = PREFIX(1:NPREFIX) // '.bl'
+ WRITE(*,1100) FILDEF
+ 1100 FORMAT(/' Enter filename: ', A)
+ READ(*,1000) FNAME
+ CALL STRIP(FNAME,NFN)
+ IF(NFN.EQ.0) FNAME = FILDEF
+ ELSE
+C------ nothing available... just ask for filename
+ CALL ASKS('Enter filename^',FNAME)
+ ENDIF
+ ENDIF
+C
+ LU = 19
+ OPEN(LU,FILE=FNAME,STATUS='UNKNOWN')
+ REWIND(LU)
+C
+ IF(KDELIM.EQ.0) THEN
+ WRITE(LU,1000)
+ & '# s x y Ue/Vinf Dstar Theta ',
+ & ' Cf H'
+C 1.23456 0.23451 0.23451 0.23451 0.012345 0.001234 0.004123 10.512
+ ELSE
+ WRITE(LU,1000)
+ & '#s' ,DELIM,
+ & 'x' ,DELIM,
+ & 'y' ,DELIM,
+ & 'Ue/Vinf',DELIM,
+ & 'Dstar' ,DELIM,
+ & 'Theta' ,DELIM,
+ & 'Cf' ,DELIM,
+ & 'H'
+ ENDIF
+C
+ CALL COMSET
+ HSTINV = GAMM1*(MINF/QINF)**2 / (1.0 + 0.5*GAMM1*MINF**2)
+C
+ DO 10 I=1, N
+ IS = 1
+ IF(GAM(I) .LT. 0.0) IS = 2
+C
+ IF(LIPAN .AND. LVISC) THEN
+ IF(IS.EQ.1) THEN
+ IBL = IBLTE(IS) - I + 1
+ ELSE
+ IBL = IBLTE(IS) + I - N
+ ENDIF
+ DS = DSTR(IBL,IS)
+ TH = THET(IBL,IS)
+ CF = TAU(IBL,IS)/(0.5*QINF**2)
+ IF(TH.EQ.0.0) THEN
+ H = 1.0
+ ELSE
+ H = DS/TH
+ ENDIF
+ ELSE
+ DS = 0.
+ TH = 0.
+ CF = 0.
+ H = 1.0
+ ENDIF
+ UE = (GAM(I)/QINF)*(1.0-TKLAM) / (1.0 - TKLAM*(GAM(I)/QINF)**2)
+ AMSQ = UE*UE*HSTINV / (GAMM1*(1.0 - 0.5*UE*UE*HSTINV))
+ CALL HKIN( H, AMSQ, HK, DUMMY, DUMMY)
+C
+ IF(KDELIM.EQ.0) THEN
+ WRITE(LU,8500) S(I), X(I), Y(I), UE, DS, TH, CF, HK
+ 8500 FORMAT(1X, 4F9.5, 3F10.6, F10.3)
+C
+ ELSE
+ WRITE(LINE,8510)
+ & S(I),DELIM,
+ & X(I),DELIM,
+ & Y(I),DELIM,
+ & UE ,DELIM,
+ & DS ,DELIM,
+ & TH ,DELIM,
+ & CF ,DELIM,
+ & HK
+ 8510 FORMAT(1X, 4(F9.5,A), 3(F10.6,A), F10.3)
+ CALL BSTRIP(LINE,NLINE)
+ WRITE(LU,1000) LINE(1:NLINE)
+ ENDIF
+C
+ 10 CONTINUE
+C
+ IF(LWAKE) THEN
+ IS = 2
+ DO 20 I=N+1, N+NW
+ IBL = IBLTE(IS) + I - N
+ DS = DSTR(IBL,IS)
+ TH = THET(IBL,IS)
+ H = DS/TH
+ CF = 0.
+ UI = UEDG(IBL,IS)
+ UE = (UI/QINF)*(1.0-TKLAM) / (1.0 - TKLAM*(UI/QINF)**2)
+ AMSQ = UE*UE*HSTINV / (GAMM1*(1.0 - 0.5*UE*UE*HSTINV))
+ CALL HKIN( H, AMSQ, HK, DUMMY, DUMMY)
+C
+ IF(KDELIM.EQ.0) THEN
+ WRITE(LU,8500) S(I), X(I), Y(I), UE, DS, TH, CF, HK
+C
+ ELSE
+ WRITE(LINE,8510)
+ & S(I),DELIM,
+ & X(I),DELIM,
+ & Y(I),DELIM,
+ & UE ,DELIM,
+ & DS ,DELIM,
+ & TH ,DELIM,
+ & CF ,DELIM,
+ & HK
+ CALL BSTRIP(LINE,NLINE)
+ WRITE(LU,1000) LINE(1:NLINE)
+ ENDIF
+ 20 CONTINUE
+ ENDIF
+C
+ CLOSE(LU)
+ RETURN
+ END ! BLDUMP
+
+
+
+ SUBROUTINE CPDUMP(FNAME1)
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*(*) FNAME1
+C
+ CHARACTER*80 FILDEF
+C
+ CHARACTER*1 DELIM
+ CHARACTER*128 LINE
+C
+ IF (KDELIM.EQ.0) THEN
+ DELIM = ' '
+ ELSEIF(KDELIM.EQ.1) THEN
+ DELIM = ','
+ ELSEIF(KDELIM.EQ.2) THEN
+ DELIM = CHAR(9)
+ ELSE
+ WRITE(*,*) '? Illegal delimiter. Using blank.'
+ DELIM = ' '
+ ENDIF
+C
+ 1000 FORMAT(8A)
+C
+ IF(FNAME1(1:1).NE.' ') THEN
+ FNAME = FNAME1
+ ELSE
+C----- no argument... get it somehow
+ IF(NPREFIX.GT.0) THEN
+C------ offer default using existing prefix
+ FILDEF = PREFIX(1:NPREFIX) // '.cp'
+ WRITE(*,1100) FILDEF
+ 1100 FORMAT(/' Enter filename: ', A)
+ READ(*,1000) FNAME
+ CALL STRIP(FNAME,NFN)
+ IF(NFN.EQ.0) FNAME = FILDEF
+ ELSE
+C------ nothing available... just ask for filename
+ CALL ASKS('Enter filename^',FNAME)
+ ENDIF
+ ENDIF
+C
+C
+ LU = 19
+ OPEN(LU,FILE=FNAME,STATUS='UNKNOWN')
+ REWIND(LU)
+C
+ IF(KDELIM.EQ.0) THEN
+ WRITE(LU,1000)
+ & '# x Cp '
+C 0.23451 0.23451
+ ELSE
+ WRITE(LU,1000)
+ & '#x', DELIM,
+ & 'Cp'
+C
+ ENDIF
+C
+ CALL COMSET
+C
+ BETA = SQRT(1.0 - MINF**2)
+ BFAC = 0.5*MINF**2 / (1.0 + BETA)
+C
+ DO 10 I=1, N
+ CPINC = 1.0 - (GAM(I)/QINF)**2
+ DEN = BETA + BFAC*CPINC
+ CPCOM = CPINC / DEN
+C
+ IF(KDELIM.EQ.0) THEN
+ WRITE(LU,8500) X(I), CPCOM
+ 8500 FORMAT(1X,2F11.5)
+ ELSE
+ WRITE(LINE,8510)
+ & X(I) , DELIM,
+ & CPCOM
+ 8510 FORMAT(1X,2(F11.5,A))
+ CALL BSTRIP(LINE,NLINE)
+ WRITE(LU,1000) LINE(1:NLINE)
+ ENDIF
+ 10 CONTINUE
+C
+ CLOSE(LU)
+ RETURN
+ END ! CPDUMP
+
+
+
+ SUBROUTINE MHINGE
+C----------------------------------------------------
+C Calculates the hinge moment of the flap about
+C (XOF,YOF) by integrating surface pressures.
+C----------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+C
+ IF(.NOT.LFLAP) THEN
+C
+ CALL GETXYF(X,XP,Y,YP,S,N, TOPS,BOTS,XOF,YOF)
+ LFLAP = .TRUE.
+C
+ ELSE
+C
+C------ find top and bottom y at hinge x location
+ TOPS = XOF
+ BOTS = S(N) - XOF
+ CALL SINVRT(TOPS,XOF,X,XP,S,N)
+ CALL SINVRT(BOTS,XOF,X,XP,S,N)
+C
+ ENDIF
+C
+ TOPX = SEVAL(TOPS,X,XP,S,N)
+ TOPY = SEVAL(TOPS,Y,YP,S,N)
+ BOTX = SEVAL(BOTS,X,XP,S,N)
+ BOTY = SEVAL(BOTS,Y,YP,S,N)
+C
+C
+ HMOM = 0.
+ HFX = 0.
+ HFY = 0.
+C
+C---- integrate pressures on top and bottom sides of flap
+ DO 20 I=2, N
+ IF(S(I-1).GE.TOPS .AND. S(I).LE.BOTS) GO TO 20
+C
+ DX = X(I) - X(I-1)
+ DY = Y(I) - Y(I-1)
+ XMID = 0.5*(X(I)+X(I-1)) - XOF
+ YMID = 0.5*(Y(I)+Y(I-1)) - YOF
+ IF(LVISC) THEN
+ PMID = 0.5*(CPV(I) + CPV(I-1))
+ ELSE
+ PMID = 0.5*(CPI(I) + CPI(I-1))
+ ENDIF
+ HMOM = HMOM + PMID*(XMID*DX + YMID*DY)
+ HFX = HFX - PMID* DY
+ HFY = HFY + PMID* DX
+ 20 CONTINUE
+C
+C---- find S(I)..S(I-1) interval containing s=TOPS
+ DO I=2, N
+ IF(S(I).GT.TOPS) GO TO 31
+ ENDDO
+C
+ 31 CONTINUE
+C---- add on top surface chunk TOPS..S(I-1), missed in the DO 20 loop.
+ DX = TOPX - X(I-1)
+ DY = TOPY - Y(I-1)
+ XMID = 0.5*(TOPX+X(I-1)) - XOF
+ YMID = 0.5*(TOPY+Y(I-1)) - YOF
+ IF(S(I) .NE. S(I-1)) THEN
+ FRAC = (TOPS-S(I-1))/(S(I)-S(I-1))
+ ELSE
+ FRAC = 0.
+ ENDIF
+ IF(LVISC) THEN
+ TOPP = CPV(I)*FRAC + CPV(I-1)*(1.0-FRAC)
+ PMID = 0.5*(TOPP+CPV(I-1))
+ ELSE
+ TOPP = CPI(I)*FRAC + CPI(I-1)*(1.0-FRAC)
+ PMID = 0.5*(TOPP+CPI(I-1))
+ ENDIF
+ HMOM = HMOM + PMID*(XMID*DX + YMID*DY)
+ HFX = HFX - PMID* DY
+ HFY = HFY + PMID* DX
+C
+C---- add on inside flap surface contribution from hinge to top surface
+ DX = XOF - TOPX
+ DY = YOF - TOPY
+ XMID = 0.5*(TOPX+XOF) - XOF
+ YMID = 0.5*(TOPY+YOF) - YOF
+ HMOM = HMOM + PMID*(XMID*DX + YMID*DY)
+ HFX = HFX - PMID* DY
+ HFY = HFY + PMID* DX
+C
+C---- find S(I)..S(I-1) interval containing s=BOTS
+ DO I=N, 2, -1
+ IF(S(I-1).LT.BOTS) GO TO 41
+ ENDDO
+C
+ 41 CONTINUE
+C---- add on bottom surface chunk BOTS..S(I), missed in the DO 20 loop.
+ DX = X(I) - BOTX
+ DY = Y(I) - BOTY
+ XMID = 0.5*(BOTX+X(I)) - XOF
+ YMID = 0.5*(BOTY+Y(I)) - YOF
+ IF(S(I) .NE. S(I-1)) THEN
+ FRAC = (BOTS-S(I-1))/(S(I)-S(I-1))
+ ELSE
+ FRAC = 0.
+ ENDIF
+ IF(LVISC) THEN
+ BOTP = CPV(I)*FRAC + CPV(I-1)*(1.0-FRAC)
+ PMID = 0.5*(BOTP+CPV(I))
+ ELSE
+ BOTP = CPI(I)*FRAC + CPI(I-1)*(1.0-FRAC)
+ PMID = 0.5*(BOTP+CPI(I))
+ ENDIF
+ HMOM = HMOM + PMID*(XMID*DX + YMID*DY)
+ HFX = HFX - PMID* DY
+ HFY = HFY + PMID* DX
+C
+C---- add on inside flap surface contribution from hinge to bottom surface
+ DX = BOTX - XOF
+ DY = BOTY - YOF
+ XMID = 0.5*(BOTX+XOF) - XOF
+ YMID = 0.5*(BOTY+YOF) - YOF
+ HMOM = HMOM + PMID*(XMID*DX + YMID*DY)
+ HFX = HFX - PMID* DY
+ HFY = HFY + PMID* DX
+C
+C---- add on TE base thickness contribution
+ DX = X(1) - X(N)
+ DY = Y(1) - Y(N)
+ XMID = 0.5*(X(1)+X(N)) - XOF
+ YMID = 0.5*(Y(1)+Y(N)) - YOF
+ IF(LVISC) THEN
+ PMID = 0.5*(CPV(1)+CPV(N))
+ ELSE
+ PMID = 0.5*(CPI(1)+CPI(N))
+ ENDIF
+ HMOM = HMOM + PMID*(XMID*DX + YMID*DY)
+ HFX = HFX - PMID* DY
+ HFY = HFY + PMID* DX
+C
+ RETURN
+ END ! MHINGE
+
+
+ SUBROUTINE VPAR
+C---------------------------------------------
+C Viscous parameter change menu routine.
+C---------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ INCLUDE 'BLPAR.INC'
+ CHARACTER*4 COMAND
+ CHARACTER*128 COMARG
+C
+ DIMENSION IINPUT(20)
+ DIMENSION RINPUT(20)
+ LOGICAL ERROR
+C
+C
+ 10 TURB = 100.0 * EXP( -(ACRIT + 8.43)/2.4 )
+ WRITE(*,1200) XSTRIP(1), XSTRIP(2), ACRIT, TURB, VACCEL,
+ & SCCON, DUXCON, GACON, GBCON, CTCON, CTRCON, CTRCEX
+ 1200 FORMAT(/' Xtr/c =', F8.4, ' top side'
+ & /' Xtr/c =', F8.4, ' bottom side'
+ & /' Ncrit =', F8.2, ' (', F6.3, ' % turb. level )'
+ & /' Vacc =', F8.4,
+ & //' Klag =', F8.4,' Uxwt =', F8.2
+ & /' A =', F8.4,' B =', F8.4,' KCt =', F8.5
+ & /' CtiniK=', F8.4,' CtiniX=', F8.4 )
+C
+C======================================================================
+C---- start of user interaction loop
+ 500 CONTINUE
+ CALL ASKC('..VPAR^',COMAND,COMARG)
+C
+ DO I=1, 20
+ IINPUT(I) = 0
+ RINPUT(I) = 0.0
+ ENDDO
+ NINPUT = 20
+ CALL GETINT(COMARG,IINPUT,NINPUT,ERROR)
+ NINPUT = 20
+ CALL GETFLT(COMARG,RINPUT,NINPUT,ERROR)
+C
+C--------------------------------------------------------------
+ IF(COMAND.EQ.' ') THEN
+ RETURN
+C
+C--------------------------------------------------------------
+ ELSEIF(COMAND.EQ.'? ') THEN
+ WRITE(*,1050)
+ 1050 FORMAT(
+ & /' <cr> Return to OPER menu'
+ & /' SHOW Display viscous parameters'
+ & /' XTR rr Change trip positions Xtr/c'
+ & /' N r Change critical amplification exponent Ncrit'
+ & /' VACC r Change Newton solution acceleration parameter'
+ & /' INIT BL initialization flag toggle'
+ & //' LAG change lag equation constants'
+ & /' GB change G-beta constants'
+ & /' CTR change initial transition-Ctau constants'
+ & /' REST restore BL calibration to baseline')
+C
+C--------------------------------------------------------------
+ ELSEIF(COMAND.EQ.'SHOW') THEN
+ GO TO 10
+C
+C--------------------------------------------------------------
+ ELSEIF(COMAND.EQ.'XTR ') THEN
+ IF(LPACC .AND. LVISC) THEN
+ WRITE(*,2100)
+ GO TO 500
+ ENDIF
+ IF(NINPUT.GE.2) THEN
+ XSTRIP(1) = RINPUT(1)
+ XSTRIP(2) = RINPUT(2)
+ ELSE
+ CALL ASKR('Enter top side Xtrip/c^',XSTRIP(1))
+ CALL ASKR('Enter bottom side Xtrip/c^',XSTRIP(2))
+ ENDIF
+ LVCONV = .FALSE.
+C
+C--------------------------------------------------------------
+ ELSEIF(COMAND.EQ.'N ') THEN
+ IF(LPACC .AND. LVISC) THEN
+ WRITE(*,2100)
+ GO TO 500
+ ENDIF
+ IF(NINPUT.GE.1) THEN
+ ACRIT = RINPUT(1)
+ ELSE
+ CALL ASKR('Enter critical amplification ratio^',ACRIT)
+ ENDIF
+ LVCONV = .FALSE.
+C
+C--------------------------------------------------------------
+ ELSEIF(COMAND.EQ.'VACC') THEN
+ IF(NINPUT.GE.1) THEN
+ VACCEL = RINPUT(1)
+ ELSE
+ CALL ASKR('Enter viscous acceleration parameter^',VACCEL)
+ ENDIF
+C
+C--------------------------------------------------------------
+ ELSEIF(COMAND.EQ.'INIT') THEN
+ LBLINI = .NOT.LBLINI
+ IF(.NOT.LBLINI) WRITE(*,*)'BLs will be initialized on next point'
+ IF( LBLINI) WRITE(*,*)'BLs are assumed to be initialized'
+ IF(.NOT.LBLINI) LIPAN = .FALSE.
+C
+C--------------------------------------------------------------
+ ELSEIF(COMAND.EQ.'LAG ') THEN
+ IF(LPACC .AND. LVISC) THEN
+ WRITE(*,2100)
+ GO TO 500
+ ENDIF
+ IF(NINPUT.GE.2) THEN
+ SCCON = RINPUT(1)
+ DUXCON = RINPUT(2)
+ ELSE
+ CALL ASKR('Enter shear lag constant^',SCCON)
+ CALL ASKR('Enter shear lag UxEQ weight^',DUXCON)
+ ENDIF
+C
+C--------------------------------------------------------------
+ ELSEIF(COMAND.EQ.'GB ') THEN
+ IF(LPACC .AND. LVISC) THEN
+ WRITE(*,2100)
+ GO TO 500
+ ENDIF
+ IF(NINPUT.GE.2) THEN
+ GACON = RINPUT(1)
+ GBCON = RINPUT(2)
+ ELSE
+ CALL ASKR('Enter G-beta constant A^',GACON)
+ CALL ASKR('Enter G-beta constant B^',GBCON)
+ ENDIF
+ CTCON = 0.5/(GACON**2 * GBCON)
+C
+C--------------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CTR ') THEN
+ IF(LPACC .AND. LVISC) THEN
+ WRITE(*,2100)
+ GO TO 500
+ ENDIF
+ IF(NINPUT.GE.2) THEN
+ CTRCON = RINPUT(1)
+ CTRCEX = RINPUT(2)
+ ELSE
+ CALL ASKR('Enter initial-Ctau constant^',CTRCON)
+ CALL ASKR('Enter initial-Ctau exponent^',CTRCEX)
+ ENDIF
+C
+C--------------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CFAC') THEN
+ IF(NINPUT.GE.1) THEN
+ CFFAC = RINPUT(1)
+ ELSE
+ CALL ASKR('Enter Cf scaling factor^',CFFAC)
+ ENDIF
+C
+C--------------------------------------------------------------
+ ELSEIF(COMAND.EQ.'REST') THEN
+ IF(LPACC .AND. LVISC) THEN
+ WRITE(*,2100)
+ GO TO 500
+ ENDIF
+ CALL BLPINI
+C
+C--------------------------------------------------------------
+ ELSE
+ WRITE(*,1000) COMAND
+ 1000 FORMAT(1X,A4,' command not recognized. Type a "?" for list')
+C
+ ENDIF
+C
+ GO TO 500
+C--------------------------------------------
+ 2100 FORMAT(/' * Polar is being accumulated.'
+ & /' * Cannot change its parameters in midstream.')
+ END ! VPAR
+
+
+
+
+ SUBROUTINE SPECAL
+C-----------------------------------
+C Converges to specified alpha.
+C-----------------------------------
+ INCLUDE 'XFOIL.INC'
+ REAL MINF_CLM, MSQ_CLM
+C
+C---- calculate surface vorticity distributions for alpha = 0, 90 degrees
+ IF(.NOT.LGAMU .OR. .NOT.LQAIJ) CALL GGCALC
+C
+ COSA = COS(ALFA)
+ SINA = SIN(ALFA)
+C
+C---- superimpose suitably weighted alpha = 0, 90 distributions
+ DO 50 I=1, N
+ GAM(I) = COSA*GAMU(I,1) + SINA*GAMU(I,2)
+ GAM_A(I) = -SINA*GAMU(I,1) + COSA*GAMU(I,2)
+ 50 CONTINUE
+ PSIO = COSA*GAMU(N+1,1) + SINA*GAMU(N+1,2)
+C
+ CALL TECALC
+ CALL QISET
+C
+C---- set initial guess for the Newton variable CLM
+ CLM = 1.0
+C
+C---- set corresponding M(CLM), Re(CLM)
+ CALL MRCL(CLM,MINF_CLM,REINF_CLM)
+ CALL COMSET
+C
+C---- set corresponding CL(M)
+ CALL CLCALC(N,X,Y,GAM,GAM_A,ALFA,MINF,QINF, XCMREF,YCMREF,
+ & CL,CM,CDP, CL_ALF,CL_MSQ)
+C
+C---- iterate on CLM
+ DO 100 ITCL=1, 20
+C
+ MSQ_CLM = 2.0*MINF*MINF_CLM
+ DCLM = (CL - CLM)/(1.0 - CL_MSQ*MSQ_CLM)
+C
+ CLM1 = CLM
+ RLX = 1.0
+C
+C------ under-relaxation loop to avoid driving M(CL) above 1
+ DO 90 IRLX=1, 12
+C
+ CLM = CLM1 + RLX*DCLM
+C
+C-------- set new freestream Mach M(CLM)
+ CALL MRCL(CLM,MINF_CLM,REINF_CLM)
+C
+C-------- if Mach is OK, go do next Newton iteration
+ IF(MATYP.EQ.1 .OR. MINF.EQ.0.0 .OR. MINF_CLM.NE.0.0) GO TO 91
+C
+ RLX = 0.5*RLX
+ 90 CONTINUE
+ 91 CONTINUE
+C
+C------ set new CL(M)
+ CALL COMSET
+ CALL CLCALC(N,X,Y,GAM,GAM_A,ALFA,MINF,QINF, XCMREF,YCMREF,
+ & CL,CM,CDP,CL_ALF,CL_MSQ)
+C
+ IF(ABS(DCLM).LE.1.0E-6) GO TO 110
+C
+ 100 CONTINUE
+ WRITE(*,*) 'SPECAL: Minf convergence failed'
+ 110 CONTINUE
+C
+C---- set final Mach, CL, Cp distributions, and hinge moment
+ CALL MRCL(CL,MINF_CL,REINF_CL)
+ CALL COMSET
+ CALL CLCALC(N,X,Y,GAM,GAM_A,ALFA,MINF,QINF, XCMREF,YCMREF,
+ & CL,CM,CDP, CL_ALF,CL_MSQ)
+ CALL CPCALC(N,QINV,QINF,MINF,CPI)
+ IF(LVISC) THEN
+ CALL CPCALC(N+NW,QVIS,QINF,MINF,CPV)
+ CALL CPCALC(N+NW,QINV,QINF,MINF,CPI)
+ ELSE
+ CALL CPCALC(N,QINV,QINF,MINF,CPI)
+ ENDIF
+ IF(LFLAP) CALL MHINGE
+C
+ RETURN
+ END ! SPECAL
+
+
+ SUBROUTINE SPECCL
+C-----------------------------------------
+C Converges to specified inviscid CL.
+C-----------------------------------------
+ INCLUDE 'XFOIL.INC'
+C
+C---- calculate surface vorticity distributions for alpha = 0, 90 degrees
+ IF(.NOT.LGAMU .OR. .NOT.LQAIJ) CALL GGCALC
+C
+C---- set freestream Mach from specified CL -- Mach will be held fixed
+ CALL MRCL(CLSPEC,MINF_CL,REINF_CL)
+ CALL COMSET
+C
+C---- current alpha is the initial guess for Newton variable ALFA
+ COSA = COS(ALFA)
+ SINA = SIN(ALFA)
+ DO 10 I=1, N
+ GAM(I) = COSA*GAMU(I,1) + SINA*GAMU(I,2)
+ GAM_A(I) = -SINA*GAMU(I,1) + COSA*GAMU(I,2)
+ 10 CONTINUE
+ PSIO = COSA*GAMU(N+1,1) + SINA*GAMU(N+1,2)
+C
+C---- get corresponding CL, CL_alpha, CL_Mach
+ CALL CLCALC(N,X,Y,GAM,GAM_A,ALFA,MINF,QINF, XCMREF,YCMREF,
+ & CL,CM,CDP, CL_ALF,CL_MSQ)
+C
+C---- Newton loop for alpha to get specified inviscid CL
+ DO 100 ITAL=1, 20
+C
+ DALFA = (CLSPEC - CL) / CL_ALF
+ RLX = 1.0
+C
+ ALFA = ALFA + RLX*DALFA
+C
+C------ set new surface speed distribution
+ COSA = COS(ALFA)
+ SINA = SIN(ALFA)
+ DO 40 I=1, N
+ GAM(I) = COSA*GAMU(I,1) + SINA*GAMU(I,2)
+ GAM_A(I) = -SINA*GAMU(I,1) + COSA*GAMU(I,2)
+ 40 CONTINUE
+ PSIO = COSA*GAMU(N+1,1) + SINA*GAMU(N+1,2)
+C
+C------ set new CL(alpha)
+ CALL CLCALC(N,X,Y,GAM,GAM_A,ALFA,MINF,QINF, XCMREF,YCMREF,
+ & CL,CM,CDP,CL_ALF,CL_MSQ)
+C
+ IF(ABS(DALFA).LE.1.0E-6) GO TO 110
+ 100 CONTINUE
+ WRITE(*,*) 'SPECCL: CL convergence failed'
+ 110 CONTINUE
+C
+C---- set final surface speed and Cp distributions
+ CALL TECALC
+ CALL QISET
+ IF(LVISC) THEN
+ CALL CPCALC(N+NW,QVIS,QINF,MINF,CPV)
+ CALL CPCALC(N+NW,QINV,QINF,MINF,CPI)
+ ELSE
+ CALL CPCALC(N,QINV,QINF,MINF,CPI)
+ ENDIF
+ IF(LFLAP) CALL MHINGE
+C
+ RETURN
+ END ! SPECCL
+
+
+ SUBROUTINE VISCAL(NITER1)
+C----------------------------------------
+C Converges viscous operating point
+C----------------------------------------
+ INCLUDE 'XFOIL.INC'
+C
+C---- convergence tolerance
+ DATA EPS1 / 1.0E-4 /
+C
+ NITER = NITER1
+C
+C---- calculate wake trajectory from current inviscid solution if necessary
+ IF(.NOT.LWAKE) THEN
+ CALL XYWAKE
+ ENDIF
+C
+C---- set velocities on wake from airfoil vorticity for alpha=0, 90
+ CALL QWCALC
+C
+C---- set velocities on airfoil and wake for initial alpha
+ CALL QISET
+C
+ IF(.NOT.LIPAN) THEN
+C
+ IF(LBLINI) CALL GAMQV
+C
+C----- locate stagnation point arc length position and panel index
+ CALL STFIND
+C
+C----- set BL position -> panel position pointers
+ CALL IBLPAN
+C
+C----- calculate surface arc length array for current stagnation point location
+ CALL XICALC
+C
+C----- set BL position -> system line pointers
+ CALL IBLSYS
+C
+ ENDIF
+C
+C---- set inviscid BL edge velocity UINV from QINV
+ CALL UICALC
+C
+ IF(.NOT.LBLINI) THEN
+C
+C----- set initial Ue from inviscid Ue
+ DO IBL=1, NBL(1)
+ UEDG(IBL,1) = UINV(IBL,1)
+ ENDDO
+C
+ DO IBL=1, NBL(2)
+ UEDG(IBL,2) = UINV(IBL,2)
+ ENDDO
+C
+ ENDIF
+C
+ IF(LVCONV) THEN
+C----- set correct CL if converged point exists
+ CALL QVFUE
+ IF(LVISC) THEN
+ CALL CPCALC(N+NW,QVIS,QINF,MINF,CPV)
+ CALL CPCALC(N+NW,QINV,QINF,MINF,CPI)
+ ELSE
+ CALL CPCALC(N,QINV,QINF,MINF,CPI)
+ ENDIF
+ CALL GAMQV
+ CALL CLCALC(N,X,Y,GAM,GAM_A,ALFA,MINF,QINF, XCMREF,YCMREF,
+ & CL,CM,CDP, CL_ALF,CL_MSQ)
+ CALL CDCALC
+ ENDIF
+C
+C---- set up source influence matrix if it doesn't exist
+ IF(.NOT.LWDIJ .OR. .NOT.LADIJ) CALL QDCALC
+C
+C---- Newton iteration for entire BL solution
+ IF(NITER.EQ.0) CALL ASKI('Enter number of iterations^',NITER)
+ WRITE(*,*)
+ WRITE(*,*) 'Solving BL system ...'
+ DO 1000 ITER=1, NITER
+C
+C------ fill Newton system for BL variables
+ CALL SETBL
+C
+C------ solve Newton system with custom solver
+ CALL BLSOLV
+C
+C------ update BL variables
+ CALL UPDATE
+C
+ IF(LALFA) THEN
+C------- set new freestream Mach, Re from new CL
+ CALL MRCL(CL,MINF_CL,REINF_CL)
+ CALL COMSET
+ ELSE
+C------- set new inviscid speeds QINV and UINV for new alpha
+ CALL QISET
+ CALL UICALC
+ ENDIF
+C
+C------ calculate edge velocities QVIS(.) from UEDG(..)
+ CALL QVFUE
+C
+C------ set GAM distribution from QVIS
+ CALL GAMQV
+C
+C------ relocate stagnation point
+ CALL STMOVE
+C
+C------ set updated CL,CD
+ CALL CLCALC(N,X,Y,GAM,GAM_A,ALFA,MINF,QINF, XCMREF,YCMREF,
+ & CL,CM,CDP,CL_ALF,CL_MSQ)
+ CALL CDCALC
+C
+C------ display changes and test for convergence
+ IF(RLX.LT.1.0)
+ & WRITE(*,2000) ITER, RMSBL, RMXBL, VMXBL,IMXBL,ISMXBL,RLX
+ IF(RLX.EQ.1.0)
+ & WRITE(*,2010) ITER, RMSBL, RMXBL, VMXBL,IMXBL,ISMXBL
+ CDP = CD - CDF
+ WRITE(*,2020) ALFA/DTOR, CL, CM, CD, CDF, CDP
+C
+ IF(RMSBL .LT. EPS1) THEN
+ LVCONV = .TRUE.
+ AVISC = ALFA
+ MVISC = MINF
+ GO TO 90
+ ENDIF
+C
+ 1000 CONTINUE
+ WRITE(*,*) 'VISCAL: Convergence failed'
+C
+ 90 CONTINUE
+ CALL CPCALC(N+NW,QINV,QINF,MINF,CPI)
+ CALL CPCALC(N+NW,QVIS,QINF,MINF,CPV)
+ IF(LFLAP) CALL MHINGE
+ RETURN
+C....................................................................
+ 2000 FORMAT
+ & (/1X,I3,' rms: ',E10.4,' max: ',E10.4,3X,A1,' at ',I4,I3,
+ & ' RLX:',F6.3)
+ 2010 FORMAT
+ & (/1X,I3,' rms: ',E10.4,' max: ',E10.4,3X,A1,' at ',I4,I3)
+ 2020 FORMAT
+ & ( 1X,3X,' a =', F7.3,' CL =',F8.4 /
+ & 1X,3X,' Cm =', F8.4, ' CD =',F9.5,
+ & ' => CDf =',F9.5,' CDp =',F9.5)
+ END ! VISCAL
+
+
+ subroutine dcpout
+ include 'XFOIL.INC'
+c
+c Computes and writes upper and lower-surface
+c Cp values at two specified x locations
+c
+c
+ x1 = 0.05
+ x2 = 0.15
+c
+ lu = 60
+ open(lu,file='dcp.out',status='old',access='append',err=10)
+ go to 20
+c
+ 10 continue
+ open(lu,file='dcp.out',status='new')
+ write(lu,*) '# ', name
+ write(lu,*) '# alpha CL ',
+ & ' Cpl05 Cpu05 dCp05 ',
+ & ' Cpl15 Cpu15 dCp15 '
+ 20 continue
+c
+ call spline(cpv,w1,s,n)
+c
+ su1 = sle + x1*(s(1)-sle)
+ sl1 = sle + x1*(s(n)-sle)
+ su2 = sle + x2*(s(1)-sle)
+ sl2 = sle + x2*(s(n)-sle)
+c
+ call sinvrt(sl1,x1,x,xp,s,n)
+ call sinvrt(su1,x1,x,xp,s,n)
+ call sinvrt(sl2,x2,x,xp,s,n)
+ call sinvrt(su2,x2,x,xp,s,n)
+c
+ cpl1 = seval(sl1,cpv,w1,s,n)
+ cpu1 = seval(su1,cpv,w1,s,n)
+ cpl2 = seval(sl2,cpv,w1,s,n)
+ cpu2 = seval(su2,cpv,w1,s,n)
+c
+ write(lu,1200) alfa/dtor, cl,
+ & cpl1, cpu1, cpl1-cpu1,
+ & cpl2, cpu2, cpl2-cpu2
+
+ 1200 format(1x, f7.3, f9.4, 8f10.5)
+c
+ close(lu)
+c
+ return
+ end