c************************************************************************
c                                                                       *
c                          PROGRAM THINFOIL                             *
c                                                                       *
c         FINITE-DIFFERENCE SOLUTION OF INCOMPRESSIBLE POTENTIAL FLOW   *
c         PAST A THIN SYMMETRIC AIRFOIL AT ZERO ANGLE OF ATTACK         *
c                                                                       *
c         SOR, SLOR and ADI Iterative schemes                           *
c                                                                       *
c         by Valery Razgonyaev for W.H. Mason                           *
c                                                                       *
c         inspired by MORAN's THINAIR code,   PAGE 326-330              *
c                                                                       *
c         code emphasizes clarity over efficiency (hopefully)           *
c                                                                       *
c         Last modification: April 9, 1992                              *
c************************************************************************

      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*1 ICASE
      COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91),
     1             X(91),Y(41),NX,NY,ILE,ITE,
     2             AX(91),BX(91),CX(91),AY(41),BY(41),CY(41),
     3             RELAX,TCMAX
      DATA IREAD /5/ IWRIT /6/ IFILE /7/
      OPEN(UNIT = IFILE, FILE = 'THINAIR.OUT', STATUS = 'UNKNOWN')

      WRITE(IWRIT,1000)
      WRITE(IFILE,1000)

c     define the mesh, initialize the potential, define bc's
      
      call setup

c     solve the equations with iteration

   10 call solve

c     output the solution on the surface

      call output

c     ask if more iterations are desired
  
      write(IWRIT,1040)
      write(IFILE,1040)
      read(IREAD,'(a)') icase
      write(IFILE,1050) icase
              KOUNT = KOUNT+maxit
      IF (icase .eq. 'Y' .or. ICASE .eq. 'y') go to 10

      close(ifile)
      
      stop

c     formats   

 1000 format(/5X,'Program THINFOIL'//5x,'Finite difference solution of'/
     1        5x,'the Laplace Eqn. over a biconvex airfoil'/)
 1040 format(/5x,'Additional iteration? (Y/N):',$)
 1050 format(5x,A1/)

      END


      subroutine solve

c     select an iteration scheme and solve the equations

      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*1 ICASE
      CHARACTER*4 SCHEME(4)
      COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91),
     1             X(91),Y(41),NX,NY,ILE,ITE,
     2             AX(91),BX(91),CX(91),AY(41),BY(41),CY(41),
     3             RELAX,TCMAX
      DATA IREAD /5/ IWRIT /6/ IFILE /7/

*************************************************************************
*                INPUT ITERATIVE SCHEME TYPE (ITYPE),                   *
*               MAXIMUM NUMBER OF ITERATIONS (maxit)                    *
*                 AND OVERRELAXATION FACTOR (RELAX)                     *
*************************************************************************
 
           SCHEME(1) = 'SOR '
           SCHEME(2) = 'SLOR'
           SCHEME(3) = 'AF-1'
           SCHEME(4) = 'AF-2'
               KOUNT = 0
 030  WRITE(IWRIT,1000)
      READ(IREAD,*)     ITYPE
      WRITE(IWRIT,1005) SCHEME(ITYPE)      
      WRITE(IFILE,1005) SCHEME(ITYPE)      
      WRITE(IWRIT,1010)
      READ(IREAD,*)     RELAX
      write(iwrit,1015) 
      read(iread,*)     maxit
      WRITE(IFILE,1020) maxit, RELAX
      WRITE(IWRIT,1030)      
      WRITE(IFILE,1030)      

*************************************************************************
*                  COMPUTATION OF THE POTENTIAL                         *
*************************************************************************
      
      DO  80    iter = 1,maxit
      
*************************************************************************
*              SET  UPDATE OUTSIDE THE COMPUTATIONAL SPACE              *
*************************************************************************
     
      DO  40      I = 2,NX
             C(I,1) =   C(I,3)
          C(I,NY+1) =   C(I,NY-1)
  40  continue
      DO  50      J = 2,NY
             C(1,J) =   C(3,J)
          C(NX+1,J) =   C(NX-1,J)
  50  continue

*************************************************************************
*                      UPDATE THE POTENTIAL                             *
*************************************************************************

      DO  60      I = 1,NX+1
      DO  60      J = 1,NY+1
           PHI(I,J) = PHI(I,J)+C(I,J)
  60  continue
     

*************************************************************************
*                      UPDATE THE RESIDUAL                              *
*************************************************************************

      do  70      I = 2,NX
      do  70      J = 2,NY
             R(I,J) = AX(I)*PHI(I-1,J)+BX(I)*PHI(I,J)+CX(I)*PHI(I+1,J)+
     1                AY(J)*PHI(I,J-1)+BY(J)*PHI(I,J)+CY(J)*PHI(I,J+1)
  70  continue
  
*************************************************************************
*                      OUTPUT ITERATION HISTORY                         *
*************************************************************************
 
      CALL STAT(KOUNT+iter,itype)

*************************************************************************
*              COMPUTATION OF THE POTENTIAL UPDATE                      *
*************************************************************************
      
      IF      (ITYPE .EQ. 1) THEN 
                                 CALL SOR
      ELSE IF (ITYPE .EQ. 2) THEN
                                 CALL SLOR
      ELSE IF (ITYPE .EQ. 3) THEN
                                 CALL AF1(iter)
      ELSE
                                 CALL AF2(iter)
      END IF
  80  continue

      return


 1000 format(//5X,'Iterative Scheme Options:'//
     1         8x,'1 - SOR '/
     2         8x,'2 - SLOR'/
     3         8x,'3 - AF-1'/
     4         8x,'4 - AF-2'/
     5        20x,'Select option:',$)
 1005 format(/5X,'Iterative Scheme Selected: - ',A4/)
 1010 format(/5X,'INPUT relaxation factor:'/
     1        8x,'Suggested values for the ADI Relaxation Factors:'/
     2'                 AF-1 - 2.0'/
     3'                 AF-2 - 1.333              :',$)
 1015 format(//5X,'INPUT number of iterations, maxit:'$)
 1020 format(/4X,I6,F20.3)
 1030 format(//33x,'Convergence History'/
     1        20x,'Residual', 27x,'Correction'/
     2        9x,'-------------------------------',
     3        6x,'-------------------------------'/
     4        3X,'ITER IMAX JMAX     MAX         AVRG',
     5        7x,     'IMAX JMAX     MAX         AVRG')

      end



      subroutine setup

c     setup mesh, initialize th epotential and define bc's

      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*1 ICASE
      CHARACTER*4 SCHEME(4)
      COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91),
     1             X(91),Y(41),NX,NY,ILE,ITE,
     2             AX(91),BX(91),CX(91),AY(41),BY(41),CY(41),
     3             RELAX,TCMAX
      DATA IREAD /5/ IWRIT /6/ IFILE /7/


c     SET UP MESH 

   
      CALL MESH
      

c     MAKE THE INITIAL GUESS FOR POTENTIAL & UPDATE 

      
      DO  10       I = 1,NX+1
      DO  10       J = 1,NY+1
            PHI(I,J) = 0.D0
              C(I,J) = 0.D0
  10  continue
 

c     DEFINE THE BOUNDARY CONDITION ON THE SURFACE                 *

 
              tcmax = 0.05D0

      do 20       i = 1,nx

              fp(i) = 0.0

      if(i .ge. ile .and. i .le. ite) then
                     fp(i) = 2.D0*tcmax*(1.D0-2.D0*X(i))
                     endif

            PHI(i,1) =-(Y(3)-Y(1))*fp(i)

  20  continue


c     output the grid and surface slope values

      WRITE(IFILE,4020) 
      WRITE(Iwrit,4020)
      i       = 1
      write(iwrit,3031) i,x(i),fp(i)
      WRITE(IFILE,3030) I,X(I),fp(i)
      DO 80 I = 2,NX
      dxgrid  = x(i) - x(i-1)
      write(iwrit,3030) i,x(i),dxgrid,fp(i)
   80 WRITE(IFILE,3030) I,X(I),dxgrid,fp(i)

      WRITE(IFILE,4030)
      WRITE(iwrit,4030)
      j       = 1
      WRITE(Iwrit,3030) J,Y(J)
      WRITE(IFILE,3030) J,Y(J)
      DO 90 J = 2,NY
      dygrid  = y(j) - y(j-1)
      WRITE(Iwrit,3030) J,Y(J),dygrid
   90 WRITE(IFILE,3030) J,Y(J),dygrid

      return

 3030 format(4x,I3,5E14.5)
 3031 format(4x,i3,e14.5,14x,3e14.5)
 4020 format(/8X,'X GRID and slope BC, fp'//
     1        6x,'i',6x,'x(i)',7x,'x(i)-x(i-1)',6x,'fp(i)')
 4030 format(/8X,'Y GRID'//
     1        6x,'j',6x,'y(j)',7x,'y(j)-y(j-1)')
 4040 format(/6X,I6,f20.3)

      end

      SUBROUTINE SOR

**************************************************************
*    POINT GAUSS-SEIDEL (SOR) ITERATIVE SCHEME               *
**************************************************************

      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91),
     1             X(91),Y(41),NX,NY,ILE,ITE,
     2             AX(91),BX(91),CX(91),AY(41),BY(41),CY(41),
     3             RELAX,TCMAX

      DO 100      I = 2,NX     
      DO 100      J = 2,NY
                AXI = AX(I)
                AYJ = AY(J)
      IF (I. EQ. 2)
     1          AXI = 0.D0
      IF (I. EQ. NX)
     1          AXI =-BX(NX)
      IF (J .EQ. 2)
     1          AYJ = 0.D0
      IF (J .EQ. NY)
     1          AYJ =-BY(NY)
             C(I,J) =-RELAX*(R(I,J)+AXI*C(I-1,J)+AYJ*C(I,J-1))/
     1                (BX(I)+BY(J))
 100  continue

      RETURN
      END
           

      
      SUBROUTINE SLOR

**************************************************************
*    LINE GAUSS-SEIDEL (SLOR) ITERATIVE SCHEME               *
**************************************************************

      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 MAIND(41),LOWD(41),UPPD(41),RHS(41)
      COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91),
     1             X(91),Y(41),NX,NY,ILE,ITE,
     2             AX(91),BX(91),CX(91),AY(41),BY(41),CY(41),
     3             RELAX,TCMAX
      
*************************************************************************
*                   SWEEP ALONG Y DIRECTION                             *
*************************************************************************
           
      DO 210      I = 2,NX
         
*************************************************************************
*     SET TRIDIAGONAL MATRIX & RHS VECTOR FOR UPDATE IN I-th COLUMN     *
*************************************************************************

      DO 200      J = 2,NY
                AXI = AX(I)
                AYJ = AY(J)
                CYJ = CY(J)
      IF (I. EQ. 2)
     1          AXI = 0.D0
      IF (I. EQ. NX)
     1          AXI =-BX(NX)
      IF (J .EQ. 2)  
     1          CYJ =-BY(2)
      IF (J .EQ. NY)  
     1          AYJ =-BY(NY)
      
         MAIND(J-1) = BX(I)+BY(J) 
          UPPD(J-1) = CYJ
          LOWD(J-1) = AYJ
           RHS(J-1) =-RELAX*(R(I,J)+AXI*C(I-1,J))
 200  continue
   
*************************************************************************
*                      CALL FOR THOMAS ALGORITHM                        *
*************************************************************************

      CALL THOMAS(LOWD,MAIND,UPPD,RHS,NY-1)
      DO 210      J = 2,NY
             C(I,J) = RHS(J-1)
 210  continue

      RETURN
      END


      
      SUBROUTINE AF1(NIT)

************************************************************************
*     APPROXIMATE FACTORISATION AF-1 (ADI) ITERATIVE SCHEME            *
************************************************************************

      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 F(91,41),MAIND(91),LOWD(91),UPPD(91),RHS(91),M
      COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91),
     1             X(91),Y(41),NX,NY,ILE,ITE,
     2             AX(91),BX(91),CX(91),AY(41),BY(41),CY(41),
     3             RELAX,TCMAX

                  K = NIT-1
                  M = 11.D0
                  P = 2.D0*((K-M*INT(K/M))/(M-1.D0)-1.D0)
               ALFA = (0.5D0*(Y(3)-Y(2)))**P

*************************************************************************
*                    INVERT OUTER OPERATOR                              *
*                             OR                                        *
*                   SWEEP ALONG I DIRECTION                             *
*************************************************************************

      DO 320      J = 2,NY
      DO 310      I = 2,NX
                AXI = AX(I)
                CXI = CX(I)
      IF (I .EQ. 2)  
     1          CXI =-BX(2)
      IF (I .EQ. NX)  
     1          AXI =-BX(NX)
 
         MAIND(I-1) = BX(I)-ALFA
          UPPD(I-1) = CXI
          LOWD(I-1) = AXI     
           RHS(I-1) = ALFA*RELAX*R(I,J)
 310  continue
    
c************************************************************************
c                      CALL FOR THOMAS ALGORITHM                        *
c************************************************************************

      CALL THOMAS(LOWD,MAIND,UPPD,RHS,NX-1)
      DO 320      I = 2,NX
             F(I,J) = RHS(I-1)
 320  continue

c************************************************************************
c                    INVERT INNER OPERATOR                              *
c                             OR                                        *
c                   SWEEP ALONG J DIRECTION                             *
c************************************************************************
           
      DO 340      I = 2,NX
      DO 330      J = 2,NY
                AYJ = AY(J)
                CYJ = CY(J)
      IF (J .EQ. 2)  
     1          CYJ =-BY(2)
      IF (J .EQ. NY)  
     1          AYJ =-BY(NY)
 
         MAIND(J-1) = BY(J)-ALFA
          UPPD(J-1) = CYJ
          LOWD(J-1) = AYJ     
           RHS(J-1) = F(I,J)
 330  continue
    
c************************************************************************
c                      CALL FOR THOMAS ALGORITHM                        *
c************************************************************************

      CALL THOMAS(LOWD,MAIND,UPPD,RHS,NY-1)
      DO 340      J = 2,NY
             C(I,J) = RHS(J-1)
 340  continue

      RETURN
      END

      
      SUBROUTINE AF2(NIT)

c************************************************************************
c       APPROXIMATE FACTORISATION AF-2 (ADI) ITERATIVE SCHEME            *
c************************************************************************

      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 F(91,41),MAIND(91),LOWD(91),UPPD(91),RHS(91),M
      COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91),
     1             X(91),Y(41),NX,NY,ILE,ITE,
     2             AX(91),BX(91),CX(91),AY(41),BY(41),CY(41),
     3             RELAX,TCMAX

                  K = NIT-1
                  M = 11.D0
                  P = (K-M*INT(K/M))/(M-1.D0)-1.D0
               ALFA = ((Y(3)-Y(2))/2.D0)**P
*************************************************************************
*                    INVERT OUTER OPERATOR                              *
*                             OR                                        *
*                   SWEEP ALONG I DIRECTION                             *
*************************************************************************

      DO 410      J = 2,NY
      DO 400      I = 2,NX
               DLOW = 0.D0
      IF (I .EQ. NX)  
     1         DLOW =-1.D0
 
         MAIND(I-1) = ALFA*(X(I+1)-X(I))+1.D0
          UPPD(I-1) =-1.D0
          LOWD(I-1) = DLOW
           RHS(I-1) =-ALFA*RELAX*R(I,J)*(X(I+1)-X(I))
 400  continue
    
*************************************************************************
*                      CALL FOR THOMAS ALGORITHM                        *
*************************************************************************

      CALL THOMAS(LOWD,MAIND,UPPD,RHS,NX-1)
      DO 410      I = 2,NX
             F(I,J) = RHS(I-1)
 410  continue

*************************************************************************
*                    INVERT INNER OPERATOR                              *
*                             OR                                        *
*                   SWEEP ALONG J DIRECTION                             *
*************************************************************************
           
      DO 430      I = 2,NX
      DO 420      J = 2,NY
                AYJ = AY(J)
                CYJ = CY(J)
                CIJ = C(I-1,J)
      IF (J .EQ. 2)  
     1          CYJ =-BY(2)
      IF (J .EQ. NY)  
     1          AYJ =-BY(NY)
      IF (I .EQ. 2)  
     1          CIJ = 0.D0
 
         MAIND(J-1) = BY(J)-ALFA/(X(I)-X(I-1))
          UPPD(J-1) = CYJ
          LOWD(J-1) = AYJ
           RHS(J-1) = F(I,J)-ALFA*CIJ/(X(I)-X(I-1))
 420  continue
    
*************************************************************************
*                      CALL FOR THOMAS ALGORITHM                        *
*************************************************************************

      CALL THOMAS(LOWD,MAIND,UPPD,RHS,NY-1)
      DO 430      J = 2,NY
             C(I,J) = RHS(J-1)
 430  continue

      RETURN
      END



      SUBROUTINE THOMAS(A,B,C,D,N)

c************************************************************************
c           INVERT TRIDIAGONAL MATRIX USING THOMAS ALGORITHM            *
c************************************************************************
      REAL*8 A(91), B(91), C(91), D(91), M
      DO 400      K = 2,N 
                  M = A(K)/B(K-1)
               B(K) = B(K)-M*C(K-1)
               D(K) = D(K)-M*D(K-1)
 400  continue
               D(N) = D(N)/B(N)
      DO 410      K = N-1,1,-1
               D(K) = (D(K)-C(K)*D(K+1))/B(K)
 410  continue
      RETURN 
      END


      
      SUBROUTINE MESH

c***************************************************
c    MESH GENERATION                               *
c***************************************************

      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91),
     1             X(91),Y(41),NX,NY,ILE,ITE,
     2             AX(91),BX(91),CX(91),AY(41),BY(41),CY(41),
     3             RELAX,TCMAX
      DATA PI /3.141592685/ IREAD /5/ IWRIT /6/ IFILE /7/

c************************************************************************
c              INPUT FIELD DIMENSIONS (XW < X < XE, 0 < Y < YN)         *
c              AND NUMBERS OF INTERVALS UPSTREAM (NUP)                  *
c              DOWNSTREAM (NDOWN), ON(NON), AND ABOVE (NABOVE)          *
c                            THE AIRFOIL                                *
c************************************************************************

      WRITE(IWRIT,2010)
      READ(IREAD,*) XW,XE,YN
      WRITE(IWRIT,2020)
      READ(IREAD,*) NUP

      WRITE(IWRIT,2022)
      READ(IREAD,*) NDOWN

      WRITE(IWRIT,2024)
      READ(IREAD,*) NON

      WRITE(IWRIT,2026)
      READ(IREAD,*) NABOVE

                ILE = NUP + 2
                ITE = ILE + NON
                 NX = ITE + NDOWN
                 NY = NABOVE + 2

      IF (NX .GT. 86 ) THEN
                WRITE(IWRIT,2030) NX
                STOP
      END IF
      IF (NY .GT. 36) THEN
                WRITE(IWRIT,2040) NY
                STOP
      END IF

      WRITE(IFILE,2010)
      WRITE(IFILE,2050) XW,XE,YN
      WRITE(IFILE,2020)
      WRITE(IFILE,2060) NUP,NDOWN,NON,NABOVE

c************************************************************************
c                         GENERATE MESH                                 *
c************************************************************************

      DO 510     I = 2,ILE
 510          X(I) = XW*(1.0D0-DCOS(PI/2.D0*(ILE-I)/DBLE(ILE-2)))
              X(1) = 2.D0*X(2)-X(3)
      DO 520     I = ILE,ITE
 520          X(I) = 0.5D0*(1.D0-COS(PI*(I-ILE)/DBLE(ITE-ILE)))
      DO 530     I = ITE,NX
 530          X(I) = XE-(XE-1.D0)*COS(PI/2.D0*(I-ITE)/DBLE(NX-ITE))
           X(NX+1) = X(NX)+(X(NX)-X(NX-1))
      DO 540     J = 2,NY
 540          Y(J) = YN*(1.D0-COS(PI/2.D0*(J-2)/DBLE(NY-2)))
              Y(1) = 2.D0*Y(2) - Y(3)
           Y(NY+1) = Y(NY) + (Y(NY) - Y(NY-1))


  
c************************************************************************
c               COMPUTE COEFFICIENTS FOR LAPLACIAN                      * 
c                                                                       *
c      L(PHI(I,J)) = AX(I)*PHI(I-1,J)+BX(I)*PHI(I,J)+CX(I)*PHI(I+1,J)+  *
c                                                                       *
c                    AY(I)*PHI(I,J-1)+BY(J)*PHI(I,J)+CY(J)*PHI(I,J+1)   *
c                                                                       *
c************************************************************************

      DO 550     I = 2,NX
             AX(I) = 2.D0/(X(I+1)-X(I-1))/(X(I)-X(I-1))
             CX(I) = 2.D0/(X(I+1)-X(I-1))/(X(I+1)-X(I))
             BX(I) =-(AX(I)+CX(I))
 550  continue
      DO 560     J = 2,NY
             AY(J) = 2.D0/(Y(J+1)-Y(J-1))/(Y(J)-Y(J-1))
             CY(J) = 2.D0/(Y(J+1)-Y(J-1))/(Y(J+1)-Y(J))
             BY(J) =-(AY(J)+CY(J))
 560  continue

      RETURN

 2010 format(/5X,
     1 'INPUT grid boundary locations; Xmin, Xmax, Ymax:',$)
 2020 format(/5X,
     1 'INPUT number of points upstream of airfoil,     Nup:',$)
 2022 format(5x,
     1 '                       downstream of airfoil, Ndown:',$)
 2024 format(5X,
     1 '                       on the airfoil,          NON:',$)
 2026 format(5X,
     1 '                       above the airfoil,    NABOVE:',$)
 2030 format(/5x,'Too many x grid points. Max is 86, your nx is',i3/)
 2040 format(/5x,'Too many y grid points. Max is 36, your ny is',i3/)
 2050 format(/4X,3F12.4)
 2060 format(/4X,4I6)

      END


      
      SUBROUTINE STAT(KOUNT,itype)

c*****************************************************
c    Solution HISTORY                                *
c*****************************************************

      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91),
     1             X(91),Y(41),NX,NY,ILE,ITE,
     2             AX(91),BX(91),CX(91),AY(41),BY(41),CY(41),
     3             RELAX,TCMAX
      DATA IREAD /5/ IWRIT /6/ IFILE /7/

c************************************************************************
c             FIND MAXIMAL RESIDUAL AND UPDATE                          *
c************************************************************************

             IMAXR = 2
             JMAXR = 2
             IMAXC = 2
             JMAXC = 2
             RAVRG = 0.D0
             CAVRG = 0.D0
      DO 600     I = 2,NX   
      DO 600     J = 2,NY
      IF (DABS(R(I,J)) .GT. DABS(R(IMAXR,JMAXR))) THEN
             IMAXR = I
             JMAXR = J
      END IF
      IF (DABS(C(I,J)) .GT. DABS(C(IMAXC,JMAXC))) THEN
             IMAXC = I
             JMAXC = J
      END IF
             RAVRG = RAVRG+R(I,J)**2
             CAVRG = CAVRG+C(I,J)**2
 600  continue
              RMAX = DABS(R(IMAXR,JMAXR))
              CMAX = DABS(C(IMAXC,JMAXC))
             RAVRG = DSQRT(RAVRG/DBLE((NX-1)*(NY-1)))
             CAVRG = DSQRT(CAVRG/DBLE((NX-1)*(NY-1)))

c      if(itype .lt. 3 .and. kount*(kount/10) .ne. kount) return 

c************************************************************************
c                    OUTPUT ITERATION HISTORY                           *
c************************************************************************

      WRITE(IFILE,3010) KOUNT,IMAXR,JMAXR,RMAX,
     1                  RAVRG,IMAXC,JMAXC,CMAX,CAVRG       
      WRITE(IWRIT,3010) KOUNT,IMAXR,JMAXR,RMAX,
     1                  RAVRG,IMAXC,JMAXC,CMAX,CAVRG      


      RETURN 

3000  format(/3X,A50)
3010  format(1X,3I5,2E12.4,3X,2I5,2E12.4)

      END



      SUBROUTINE output

c************************************************************************
c          COMPUTE AND PRINT PRESSURE DISTRIBUTION ON THE AIRFOIL       *
c************************************************************************

      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91),
     1             X(91),Y(41),NX,NY,ILE,ITE,
     2             AX(91),BX(91),CX(91),AY(41),BY(41),CY(41),
     3             RELAX,TCMAX
      DATA IREAD /5/ IWRIT /6/ IFILE /7/

      WRITE(IWRIT,4000)
      WRITE(IFILE,4000)
      ILEP      = ILE + 1
      DO 700  I = ILEP,ITE
      XP        = 0.5D0*(X(I) + X(I-1))
      U         = (PHI(I,2) - PHI(I-1,2))/(X(I) - X(I-1))
      CP        = - 2.D0*U
      WRITE(IWRIT,4010) I,XP,U,CP
      WRITE(IFILE,4010) I,XP,U,CP
  700 continue

 4000 format(//3X,' PRESSURE DISTRIBUTION ON THE AIRFOIL'//5X,
     1 'I',10X,'x/c',13X,'U',14X,'Cp')
 4010 format(3X,I3,3F15.5)

      RETURN
      END























