      PROGRAM ARC2D
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C*************************************************************************
C*************************************************************************
C*************************************************************************
C                                                        *****************
C   ARC2D :  EULER VERSION FOR BENCHMARK  USE.           *****************
C            11/05/86  THOMAS H. PULLIAM                 *****************
C            NASA AMES RESEARCH CENTER                   *****************
C                                                        *****************
C*************************************************************************
C  ARC2D SOLVES THE  EULER EQUATIONS IN GENERALIZED CURVILINEAR COORDINATES
C  USING AN IMPLICIT FINITE DIFFERENCE ALGORITHM WITH APPROXIMATED 
C  FACTORIZATION AND A DIAGONALIZATION OF THE IMPLICIT OPERATORS.
C
C
C  EULER EQUATIONS :   Q_T + E_XI + F_ETA = 0
C
C      -1 | RHO |       -1|     RHO U     |      -1|    RHO V        |
C Q = J   |RHO U|, E = J  |RHO UU + XI_X P|, F = J |RHO U V + ETA_X P|
C         |RHO V|         |RHO VU + XI_Y P|        |RHO V V + ETA_Y P|
C         |  E  |         |U(E+P) - XI_T P|        |V(E+P)  - ETA_T P|
C
C    U = XI_T + XI_X U + XI_Y V,    V = ETA_T + ETA_X U + ETA_Y V
C    
C     -1
C    J   = (X_XI Y_ETA - X_ETA Y_XI)
C
C    XI_X  =  J Y_ETA, XI_Y  = -J X_ETA,  XI_T  = -X_T  XI_X - Y_T XI_Y
C   ETA_X = -J Y_XI,  ETA_Y =  J X_XI,   ETA_T = -X_T ETA_X - Y_T ETA_Y
C
C    P = (GAMMA-1) (E - 0.5 RHO (U^2 + V^2))  PRESSURE
C    U IS X COMPONENT OF VELOCITY,  V IS Y COMPONENT OF VELOCITY
C    E IS TOTAL ENERGY,  RHO IS DENSITY
C    GAMMA IS THE RATIO OF SPECIFIC HEATS (USUALLY 1.4)
C
C    NONDIMENSIONALIZATION :   RHO_INFINITY FOR RHO
C                              A_INFINITY (SPEED OF SOUND) FOR U AND V
C                              RHO_INFINITY A_INFINITY^2 FOR E
C                              REFERENCE LENGTH IS L 
C                                      (USUALLY CHORD FOR AN AIRFOIL)
C   NOTE THIS GIVES P AT INFINITY = 1./GAMMA
C
C
C   THE NUMERICAL ALGORITHM IS THE CLASS OF TWO TIME LEVEL 1ST OR 2ND
C   ORDER ACCURATE TIME IMPLICIT SCHEMES WITH SECOND ORDER CENTRAL
C   DIFFERENCES IN COMPUTATIONAL SPACE.
C   THE IMPLICIT SIDE IS APPROXIMATELY FACTORED AND DIAGONALIZED.
C   SEE: PULLIAM T. H. AND CHAUSSEE D. S. , 
C        "A DIAGONAL FORM OF AN IMPLICIT APPROXIMATE FACTORIZATION 
C         ALGORITHM", JCP 39, 1981, P347.
C
C   THE BENCHMARK TEST CASE SET UP HERE IS FOR THE UPPER HALF OF A
C   BICONVEX AIRFOIL SITTING IN A STRETCH RECTANGULAR DOMAIN.
C   THE GRID IS GENERATED INTERNALLY, SEE BICONG.
C   THE BOUNDARY CONDITIONS ARE THIN AIRFOIL TANGENCY CONDITIONS ON THE
C   AIRFOIL SURFACE,  QUASI-ONE-DIMENSIONAL RIEMANN INVARIANT CONDITIONS
C   AT INFLOW, OUTFLOW AND THE TOP OF THE DOMAIN, SEE BC.
C
C   INPUT IS HANDLED IN INPUT AND INITIA.
C   OUPUT IS HANDLED IN IOALL.
C   CONDITIONS FOR THE SUPPLIED TEST CASES ARE VARIOUS GRID SIZES AT
C   A MACH NUMBER = 0.8.  THIS GIVES TRANSONIC FLOW AND A SHOCK FOR 
C   A THICKNESS OF THE BICONVEX = 0.1
C
C   THE CODE IS SELF DOCUMENTED, SEE THE DEFINITIONS GIVEN BELOW FOR MORE 
C   DETAIL.  CONTACT THOMAS H. PULLIAM, MS 202A-1 NASA AMES RESEARCH CENTER,
C   MOFFETT FIELD, CA.,  (415-694-6417) (EMAIL: PULLIAM@PRANDTL.NAS.NASA.GOV)
C   TO REPORT PROBLEMS OR FOR QUESTIONS.
C
C*************************************************************************
C*************************************************************************
C
C     COMMON/BASE/
C
C------------------------------------------------------------------------
C                       INDICES  J,K
C------------------------------------------------------------------------
C       JMAX:   TOTAL NUMBER OF POINTS IN XI DIRECTION   (DEFAULT=193)
C       KMAX:   TOTAL NUMBER OF POINTS IN ETA DIRECTION  (DEFAULT=33 )
C       JM:     JMAX - 1
C       KM:     KMAX - 1
C       JBEGIN, JEND : ABSOLUTE LOWER AND UPPER LIMITS ON J
C                      TYPICALLY 1 AND JMAX
C       KBEGIN, KEND : ABSOLUTE LOWER AND UPPER LIMITS ON K
C                      TYPICALLY 1 AND KMAX
C------------------------------------------------------------------------
C     AUTOMATIC PERIODIC AND NONPERIODIC LOGIC IS AVAILABLE.
C
C       PERIODIC    SWITCH FOR NONPERIODIC VRS PERIODIC CODING
C               PERIODIC = F  NONPERIODIC OPTION
C               PERIODIC = T   PERIODIC OPTION           (DEFAULT=FALSE)
C       JPLUS   INDEX ARRAY FOR J DIFFERENCING
C               JPLUS(J) = J+1   
C                 JPLUS(JMAX) = JMAX FOR PERIODIC = F
C                 JPLUS(JMAX) = 1    FOR PERIODIC = T
C       JMINU   INDEX ARRAY FOR J DIFFERENCING
C               JMINU (J) = J+1
C                 JMINU (1) = 1    FOR PERIODIC = F
C                 JMINU (1) = JMAX FOR PERIODIC = T
C       JLOW    LOWER DO LOOP LIMIT IN J DIRECTION
C                  JLOW = 2  FOR PERIODIC = F
C                  JLOW = 1  FOR PERIODIC = T
C       JUP     UPPER DO LOOP LIMIT IN J DIRECTION
C                  JUP = JMAX-1  FOR PERIODIC = F
C                  JUP = JMAX    FOR PERIODIC = T
C       KLOW    LOWER DO LOOP LIMIT IN K DIRECTION
C                  KLOW = 2  TYPICALLY
C       KUP     UPPER DO LOOP LIMIT IN K DIRECTION
C                  KUP = KMAX-1  TYPICALLY
C------------------------------------------------------------------------
C
C      NP:      CALLS IOALL(15) TO SHOW RHO, U, V, E, P, CP, MACH, X, Y
C               (DEFAULT=10000)
C
C       DT:     TIME STEP            (DEFAULT=5.0)
C               EITHER A TRUE DT FOR TIME ACCURATE COMPUTATION
C               OR A SCALING FACTOR FOR THE SPATIALLY VARIABLE TIME STEP
C               [ IN STEADY STATE CALCULATION DT IS ORDER(1.0) ]
C
C       CP2TIM:  ELAPSED CPU TIME OF THE RUN.  
C                 NOTE CHECK FUNCTION SECOND FOR INSTRUCTIONS FOR 
C                 MACHINE DEPENDENT CHANGES.
C------------------------------------------------------------------------
C       FSMACH: FREE STREAM MACH NUMBER                (DEFAULT=0.397)
C               READ IN IN IOALL(1)
C       ALPHA:  ANGLE OF ATTACK OF AIRFOIL             (DEFAULT=0.0 )
C               READ IN IN IOALL(1)
C       GAMMA:  GAS CONSTANT ( 1.4 )
C       GAMI:   GAMMA - 1
C       PI  :   4.*ATAN(1.0)
C------------------------------------------------------------------------
C                 DISSIPATION PARAMETERS
C------------------------------------------------------------------------
C
C       DIS2X:  SECOND ORDER X COEFFICIENT FOR FILTER       (DEFAULT=0.0)
C       DIS2Y:  SECOND ORDER Y COEFFICIENT FOR FILTER       (DEFAULT=0.0)
C               READ IN IN IOALL(1) ONLY USED FOR METH = 2,3  
C               USED IN FILERX, FILERY, STEPFX, STEPFY 
C       DIS4X:  FOURTH ORDER X COEFFICIENT FOR FILTER       (DEFAULT=0.64)
C       DIS4Y:  FOURTH ORDER Y COEFFICIENT FOR FILTER       (DEFAULT=0.64)
C               READ IN IN IOALL(1) ONLY USED FOR METH = 2,3  
C               USED IN FILERX, FILERY,  STEPFX, STEPFY 
C------------------------------------------------------------------------
C       JACDT:  SWITCH TO TURN ON VARIABLE TIME STEP.    (DEFAULT = 1)
C                      FOR UNSTEADY PROBLEMS USE CONSTANT TIME STEP.
C                      FOR STEADY - STATE USE VARIABLE TIME STEP.
C               READ IN IN IOALL(1)
C                       =0 : CONSTANT TIME STEP
C                       =1 : TIME STEP SCALED BY JACOBIAN
C                       =2 : TIME STEP BASED ON MINIMUM 1/EIGENVALUES
C                       =3 : TIME STEP BASED CONSTANT CFL AT EACH POINT
C
C------------------------------------------------------------------------
C------------------------------------------------------------------------
C       RESID:  ASSIGNED TO IN RESIL2 OUTPUT IN IOALL(8) ONTO TAPE7
C                       OUTPUT IN IOALL(5) TO TAPE6
C                       ZEROED IN INITIA
C                       2-NORM OF RESIDUAL
C
C       IPRINT: OUTPUT SWITCH FOR FREQUENCY OF RESIDUAL HISTORY OUTPUT
C               ON UNIT 6 ( MAIN )                      (DEFAULT=1)
C
C------------------------------------------------------------------------
C                     GRID PARAMETERS
C------------------------------------------------------------------------
C      DYM:    MINIMUM SPACING AT BODY FOR BICONVEX AIRFOIL
C      YMAX:      UPPER DIMENSION IN CORD LENGTHS FOR BICONVEX GRID
C      XMIN:      LEFT DIMENSION IN CORD LENGTHS FOR BICONVEX GRID
C      XMAX:      RIGHT DIMENSION IN CORD LENGTHS FOR BICONVEX GRID
C      THICK:      THICKNESS RATIO FOR BICONVEX AIRFOIL 
C      JTAIL1:      FIRST SOLID BODY POINT 
C      JTAIL2: LAST SOLID BODY POINT 
C       SLOPE:  SLOPE OF THE AIRFOIL, USED FOR  BICONVEX CASE 
C               SET IN INITIA, USED IN BC
C
C------------------------------------------------------------------------
C                    TIME ACCURACY PARAMETERS
C------------------------------------------------------------------------
C      A CLASS OF THREE STEP FIRST OR SECOND ORDER TIME ACCURACTE SCHEMES 
C      CAN BE EMPLOYED.  FOR STEADY STATE EULER IMPLICIT PERFORMS THE BEST.
C      FOR UNSTEADY, 3PT IMPLICIT IS A GOOD CHOICE.
C
C       PHIDT         TIME INTEGRATION PARAMETERS : 
C       THETAD        PHIDT = 0   THETAD  = 1    EULER IMPLICIT (DEFAULT)
C                     PHIDT = 0   THETAD  = 1/2  TRAP. IMPLICIT
C                     PHIDT = 1/2 THETAD  = 1    3PT   IMPLICIT
C
C------------------------------------------------------------------------
C      NUMITE :      ITERATION COUNTER FOR MAIN STEPS
C       ISTART:         THE STARTING ITERATION,
C                               READ IN FROM THE RESTART
C       NSTEPS:          TOTAL NUMBER OF ITERATION STEPS
C----------------------------------------------------------------------
C-----------------------   VARIABLES   --------------------------------
C----------------------------------------------------------------------
C      Q:      CONSERVATIVE VARIABLES
C            Q(I,J, 1)     DENSITY / DETERMINANT OF METRIC JACOBIAN
C            Q(I,J, 2)     DENSITY * U / DETERMINANT OF METRIC JACOBIAN
C            Q(I,J, 3)     DENSITY * V / DETERMINANT OF METRIC JACOBIAN
C            Q(I,J, 4)     ENERGY / DETERMINANT OF METRIC JACOBIAN
C
C       PRESS :  PRESSURE = GAMI*(Q(4) - 0.5*(Q(2)**2+Q(3)**2)/Q(1))
C                NOTE:  HAS JACOBIAN IN IT
C       SNDSP :  A = SQRT(GAMMA*PRESS/Q(1))
C                NOTE:  DOESN'T HAVE JACOBIAN IN IT
C      S:      
C          USED FOR EXPLICIT RIGHT HAND SIDE IN RHS,
C           ALTERED IN VISRHS,
C          SMOOTHED IN SMOOTH,FILTER 
C           USED FOR INTERMEDIATE QUANTITIES IN THE AF SCHEME IN STEP, ETC. 
C           AND ULTIMATELY HOLDS THE CHANGE IN Q.
C      XY:      METRIC TRANSFORMATIONS
C            XY(I,J, 1)     D XI / DX
C            XY(I,J, 2)     D XI / DY
C            XY(I,J, 3)     D ETA / DX
C            XY(I,J, 4)     D ETA / DY
C      XYJ:      JACOBIANS OF METRIC TRANSFORMATIONS
C            = ( DXI/DX * DETA/DY ) - ( DXI/DY * DETA/DX )
C      XIT:      D XI / DT
C      ETT:      D ETA / DT
C      DS:      SPATIALLY VARIABLE TIME STEP, SCALED ON METRIC JACOBIAN,
C                  SET IN INITIA, USED IN STEP, STEPF, 
C                                         EIGEN, FLMATX2, FLMATY2
C      X:      CARTESIAN X COORDINATES OF THE GRID.
C      Y:      CARTESIAN Y COORDINATES OF THE GRID.
C
C
C   TO CHANGE DIMENSIONS:::  MODIFY ALL THE PARAMETRER STATEMENTS
C                            MAXJ = MAXIMUM J DIMENSION
C                            MAXK = MAXIMUM K DIMENSION
C                            MAXJK = PRODUCT MAXJ*MAXK
C                            
      PARAMETER (MAXJ=165, MAXK=48, MAXJK=MAXJ*MAXK  )
      COMMON/BASE/
     1  JMAX,       KMAX,     JM,          KM,          JBEGIN,   JEND,
     1  KBEGIN,     KEND,     JPLUS(999),  JMINU (999), JLOW,     JUP,
     1  KLOW,       KUP,      PERIDC  ,    NP,          DT,   CP2TIM,
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,    
     1  THETAD ,    RESID,    JACDT,       IPRINT,      NPCP,
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS
      LOGICAL PERIDC  
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK
C
      DIMENSION Q(MAXJK*4)
      DIMENSION PRESS(MAXJK),SNDSP(MAXJK)
      DIMENSION S(MAXJK*4),XY(MAXJK*4),XYJ(MAXJK)
      DIMENSION XIT(MAXJK),ETT(MAXJK),DS(MAXJK)
      DIMENSION X(MAXJK),Y(MAXJK)
      DIMENSION DELTQN(MAXJK*4),SLOPE(MAXJ)
      DIMENSION UU(MAXJK),VV(MAXJK),CCX(MAXJK),CCY(MAXJK)
      DIMENSION COEF4X(MAXJK),COEF4Y(MAXJK)
      DIMENSION COEF2X(MAXJK),COEF2Y(MAXJK)
C
C
      COMMON/WORKSP/WORK2(MAXJ,92)
      COMMON/WORK/WORK1(MAXJK*10)
C
      EXTERNAL ZECOND
C      
      CALL INITIM()
C
C  INPUT DATA 
C
       CALL INPUT(MAXJ,MAXJK)
C
C           INITIALIZE CPU TIME COUNTER
C
       CPTIM0 = ZECOND(DUMMY)
C
          PRINT *, '--------------------------------------------'
          PRINT *,'        CPUTIM AT START = ',CPTIM0
          PRINT *, '--------------------------------------------'
C
C
C  INITIALIZE VARIABLES
C
C   SET JDIMEN,KDIMEN
C
      JDIMEN  = JMAX
      KDIMEN  = KMAX
C
C
      CALL INITIA(JDIMEN,KDIMEN,Q,PRESS,SNDSP,S,XY,XYJ,
     *                  XIT,ETT,DS,X,Y,SLOPE)
C
C
C  PERFORM INTEGRATION
C
         ISTART = 1
         IEND   = NSTEPS
         DO 10 NUMITE  = ISTART,IEND
       
         CALL INTEGR(JDIMEN,KDIMEN,Q,S,PRESS,SNDSP,SLOPE,
     +           X,Y,XY,XYJ,XIT,ETT,DS,DELTQN,
     +           UU,VV,CCX,CCY,COEF2X,COEF2Y,COEF4X,COEF4Y)

C
C  COMPUTE CPUTIM
C
         CP2TIM = ZECOND(DUMMY) - CPTIM0
C
C*********************************************************************
C***  RESIDUAL, SS PTS,                                            ***
C*********************************************************************
      IP = 5
         IF( MOD( NUMITE  - ISTART + 1, IPRINT) .EQ. 0 ) 
     * CALL IOALL(IP,JDIMEN,KDIMEN,Q,PRESS,SNDSP,XY,XYJ,X,Y)
C*********************************************************************
C***            WRITE OUT RHO, U, V, E, P, CP, MACH, X, Y          ***
C***                     ONTO TAPE6 EVERY NP ITERATIONS            ***
C*********************************************************************
      IP = 15
         IF ( MOD( NUMITE  - ISTART + 1, NP) .EQ. 0 ) 
     * CALL IOALL(IP,JDIMEN,KDIMEN,Q,PRESS,SNDSP,XY,XYJ,X,Y)
C
C*********************************************************************
10       CONTINUE
C*********************************************************************
C
C
C  FINAL OUTPUT
C
C          Q SOLUTION FILE FOR PLOTTING
      CALL  IOALL(10,JDIMEN,KDIMEN,Q,PRESS,SNDSP,XY,XYJ,X,Y)
C          XY GRID FILE FOR PLOTTING
      CALL  IOALL(16,JDIMEN,KDIMEN,Q,PRESS,SNDSP,XY,XYJ,X,Y)
C          CP VRS X
      CALL  IOALL(6,JDIMEN,KDIMEN,Q,PRESS,SNDSP,XY,XYJ,X,Y)
C***********************************************************************
      CALL SHOTIM()
      STOP
      END
