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)
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 GETINT(COMARG,IINPUT,NINPUT,ERROR)
      NINPUT = 20
      CALL GETFLT(COMARG,RINPUT,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.'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, 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,
     &                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),
     &             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),
     &             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),
     &      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),
     &                 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, 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,
     &              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),
     &                 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
       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 = 0
        IA2 = 0
        CALL POLWRIT(6,' ',ERROR, .TRUE.,
     &                 NAX, IA1,IA2, CPOL(1,1,IP), IPOL,NIPOL,
     &                 REYNP1(IP),MACHP1(IP),ACRITP(IP),XSTRIPP(1,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)
       ENDIF
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.'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.'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) INT(REINF1)
       IF(RETYP.EQ.2) WRITE(*,1200) INT(REINF1), ' / sqrt(CL)'
       IF(RETYP.EQ.3) WRITE(*,1200) INT(REINF1), ' / CL'
      ENDIF
C
      RETURN
C
 1100 FORMAT(1X,'M  =' , F10.4, A)
 1200 FORMAT(1X,'Re =' , I10  , 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
 1000 FORMAT(A,A)
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
      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
C
      CALL COMSET
      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
        WRITE(LU,8500) S(I), X(I), Y(I), UE, DS, TH, CF, HK
 8500   FORMAT(1X,4F9.5,3F10.6,F10.3)
  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
          WRITE(LU,8500) S(I), X(I), Y(I), UE, DS, TH, CF, HK
 20     CONTINUE
      ENDIF
C
      CLOSE(LU)
      RETURN
      END ! BLDUMP



      SUBROUTINE CPDUMP(FNAME1)
      INCLUDE 'XFOIL.INC'
      CHARACTER*(*) FNAME1
C
      CHARACTER*80 FILDEF
C
 1000 FORMAT(A)
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
      WRITE(LU,1000)
     & '#    x        Cp  '
C         0.23451  0.23451
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
        WRITE(LU,8500) X(I), CPCOM
 8500   FORMAT(1X,2F9.5)
  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.'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
