      SUBROUTINE BC(JDIM,KDIM,Q,PRESS,SNDSP,XY,XIT,ETT,XYJ,X,Y,SLOPE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAXJ=165, MAXK=48, MAXJK=MAXJ*MAXK  )
C
      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(JDIM,KDIM,4)
      DIMENSION PRESS(JDIM,KDIM),SNDSP(JDIM,KDIM)
      DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM)
      DIMENSION XIT(JDIM,KDIM),ETT(JDIM,KDIM)
      DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM),SLOPE(JDIM)
      COMMON/INFNTY/RHOINF,UINF,VINF,EINF,PINF
C
      COMMON/WORKSP/F(MAXJ),WORK0(MAXJ,91)
C
C    BICONVEX AIFOIL BC
C
      STRTIT = 12.
      T = ( NUMITE  -1.)/ STRTIT
      IF( T .GT.1.) T = 1.
      SCAL = (10. -15.*T + 6.*T*T) *T**3
C
C  U AND V ON LOWER SURFACE
C    U IS EXTRAPOLATED 
C    V IS SET TO SATISFY TANGENCY
C    A RECTILINEAR MESH IS ASSUMED
C
      DO 25 J = 2,JM
C
C         FIRST ORDER EXTRAPOLATION OF U ( DIVIDE BY RHO/J )
C          U2 = 2. * Q(J,2,2)/Q(J,2,1) - Q(J,3,2)/Q(J,3,1)
C
C         ZERO'TH ORDER EXTRAPOLATION OF U ( DIVIDE BY RHO/J )
          U2 = Q(J,2,2) / Q(J,2,1)
C
C         ZERO'TH ORDER EXTRAPOLATION OF U
          Q(J,1,2)  = U2
C
C         TANGENCY
          Q(J,1,3) = U2 * SLOPE(J) * SCAL
C
C
25    CONTINUE
C
C  SATISFY MOMENTUM RELATION FOR NORMAL PRESSURE DERIVATION
C
      K = 1
      DO 30 J=2,JM
C
C         PRESSURE IS EXTRAPOLATED TO FIRST ORDER INTO F
C
         P2 = GAMI*(Q(J,2,4)-0.5*(Q(J,2,2)**2+Q(J,2,3)**2)/Q(J,2,1))
         P3 = GAMI*(Q(J,3,4)-0.5*(Q(J,3,2)**2+Q(J,3,3)**2)/Q(J,3,1))
         F(J) = (2.*P2*XYJ(J,2)-P3*XYJ(J,3))/XYJ(J,1)
C
C RESCALE Q2 AND Q3 TO NEW DENSITY SUCH THAT WHEN REFORMING PRESSURE
C      FROM Q1 TO Q4, THE NORMAL DERIVATIVE OF PRESSURE IS SATIFIED
C
C         ZERO'TH ORDER EXTRAPOLATION OF RHO
          Q( J, K, 1)  = Q( J, K+1, 1) * XYJ( J, K+1) / XYJ( J, K)
C
C          U EXTRAP ABOVE, THIS PUTS IN RHO ( WITH METRIC JACOBIAN )
          Q(J,K,2) = Q(J,K,2)*Q(J,K,1)
C
C          V SET ABOVE, THIS PUTS IN RHO ( WITH METRIC JACOBIAN )
          Q(J,K,3) = Q(J,K,3)*Q(J,K,1)
C
C          ENERGY
C           RECALCULATE THE ENERGY USING THE NEW RHO, U, V 
C           AND EXTRAPOLATED PRESSURE.
           Q(J,K,4) = F(J)/GAMI +.5*( Q(J,K,2)**2+Q(J,K,3)**2) /Q(J,K,1)
C
   30 CONTINUE
C.....................................................................
C   INFLOW/OUTFLOW BOUNDARIES
C               CHECK FOR SUPERSONIC FREE STREAM
C.....................................................................
C
      IF(FSMACH.GT.1.)GO TO 75
C
C
C  INFLOW BOUNDARY
C
          GI = 1./GAMMA
          GM1I = 1./GAMI
C
         AINF = SQRT(GAMMA*PINF/RHOINF)
         HSTFS = 1./GAMI + 0.5*FSMACH**2
C
      J = 1
      DO 50 K = 2,KM
      SNORM = 1./SQRT(XY(J,K,1)**2+XY(J,K,2)**2)
      XY1H = XY(J,K,1)*SNORM
      XY2H = XY(J,K,2)*SNORM
C
C  EXTRAPOLATE TEMPORARY VARIABLES
C 
      XYJM1 = XYJ(J+1,K)
      RHOEXT = Q(J+1,K,1)*XYJM1 
      RJM1 = 1./Q(J+1,K,1)
      UEXT = Q(J+1,K,2)*RJM1 
      VEXT = Q(J+1,K,3)*RJM1 
      EEXT = Q(J+1,K,4)*XYJM1 
      PEXT = GAMI*(EEXT - 0.5*RHOEXT*(UEXT**2+VEXT**2))
C
C  FIX INCOMING RIEMANN INVARIANT
      RFIX = XY1H*UINF + XY2H*VINF + 2.*AINF*GM1I
      REXT = XY1H*UEXT + XY2H*VEXT - 2.*SQRT(GAMMA*PEXT/RHOEXT)*GM1I
C
C              COMPUTE FLOW VARIABLES BASED ON ABOVE CALC
C
      QN = (RFIX + REXT)*0.5
      CSPE = (RFIX - REXT)*GAMI*0.25
      C2 = CSPE**2
C
           QT = (-XY2H*UINF + XY1H*VINF)
           ENTRO = GAMMA
C
C             COMPUTE FLOW VARIABLES
C
      U = XY1H*QN - XY2H*QT
      V = XY2H*QN + XY1H*QT
      Q(J,K,1) = (C2*ENTRO*GI)**GM1I
      PRES = C2*Q(J,K,1)*GI
C
C              ADD JACOBIAN
C
          RJ = 1./XYJ(J,K)
          Q(J,K,1) = Q(J,K,1)*RJ
          Q(J,K,2) = Q(J,K,1)*U
          Q(J,K,3) = Q(J,K,1)*V
          Q(J,K,4) = PRES*GM1I*RJ + 0.5*Q(J,K,1)*(U**2+V**2)
50    CONTINUE
C
C
C  OUTFLOW BOUNDARY
C
C
      J = JMAX
      DO 55 K = 2,KM
      SNORM = 1./SQRT(XY(J,K,1)**2+XY(J,K,2)**2)
      XY1H = XY(J,K,1)*SNORM
      XY2H = XY(J,K,2)*SNORM
C
C  EXTRAPOLATE TEMPORARY VARIABLES
C 
      XYJM1 = XYJ(J-1,K)
      RHOEXT = Q(J-1,K,1)*XYJM1 
      RJM1 = 1./Q(J-1,K,1)
      UEXT = Q(J-1,K,2)*RJM1 
      VEXT = Q(J-1,K,3)*RJM1 
      EEXT = Q(J-1,K,4)*XYJM1 
      PEXT = GAMI*(EEXT - 0.5*RHOEXT*(UEXT**2+VEXT**2))
C
C  FIX INCOMING RIEMANN INVARIANT
      RFIX = XY1H*UINF + XY2H*VINF - 2.*AINF*GM1I
      REXT = XY1H*UEXT + XY2H*VEXT + 2.*SQRT(GAMMA*PEXT/RHOEXT)*GM1I
C
C              COMPUTE FLOW VARIABLES BASED ON ABOVE CALC
C
      QN = (RFIX + REXT)*0.5
      CSPE = (REXT - RFIX)*GAMI*0.25
      C2 = CSPE**2
C
           QT = (-XY2H*UEXT + XY1H*VEXT)
           ENTRO = RHOEXT**GAMMA/PEXT
C
C             COMPUTE FLOW VARIABLES
C
      U = XY1H*QN - XY2H*QT
      V = XY2H*QN + XY1H*QT
      Q(J,K,1) = (C2*ENTRO*GI)**GM1I
      PRES = C2*Q(J,K,1)*GI
C
C              ADD JACOBIAN
C
          RJ = 1./XYJ(J,K)
          Q(J,K,1) = Q(J,K,1)*RJ
          Q(J,K,2) = Q(J,K,1)*U
          Q(J,K,3) = Q(J,K,1)*V
          Q(J,K,4) = PRES*GM1I*RJ + 0.5*Q(J,K,1)*(U**2+V**2)
55    CONTINUE

      GO TO 100
C
75    CONTINUE
C
C  FOR SUPERSONIC INFLOW  ALL VARIABLES FIXED TO FREE STREAM
C
       J = 1
       DO 80 K = 1,KMAX
       RJJ = 1./XYJ(J,K)
       Q(J,K,1) = RHOINF*RJJ
       Q(J,K,2) = UINF*Q(J,K,1)      
       Q(J,K,3) = VINF*Q(J,K,1)      
       Q(J,K,4) = EINF*RJJ
80     CONTINUE
C
C   FOR SUPERSONIC OUTFLOW EXTRAPOLATE ALL VARIABLES
C  
       J = JMAX
       DO 85 K = 1,KMAX
       RJJ = XYJ(J-1,K)/XYJ(J,K)
       Q(J,K,1) = Q(J-1,K,1)*RJJ
       Q(J,K,2) = Q(J-1,K,2)*RJJ
       Q(J,K,3) = Q(J-1,K,3)*RJJ
       Q(J,K,4) = Q(J-1,K,4)*RJJ
85     CONTINUE
C
100    CONTINUE
C
C  TOP BOUNDARY
C
C        *****   FAR-FIELD STUFF  *****
C
C             TOP BOUNDARY
C
C
         AINF = SQRT(GAMMA*PINF/RHOINF)
         HSTFS = 1./GAMI + 0.5*FSMACH**2
         GM1I  = 1./GAMI
         GI    = 1./GAMMA
         K = KMAX
         IF(FSMACH.LT.1.0)THEN
C
C                         SUBSONIC FREESTREAM
C
      DO 60 J = 2,JM
         UF = UINF
         VF = VINF
         AF2 = GAMI*(HSTFS - .5*(UF**2+VF**2))
         AF = SQRT(AF2)
C
C              CHOOSE A REFERENCE FRAME IN TERMS OF NORMAL AND 
C              TANGENTIAL COMPONENTS
C
C                   METRIC TERMS
C  
      SNORM = 1./SQRT(XY(J,K,3)**2+XY(J,K,4)**2)
      XY3H = XY(J,K,3)*SNORM
      XY4H = XY(J,K,4)*SNORM
C
C                  CHECK FOR INFLOW OR OUTFLOW
C             FOR INFLOW : THREE VARIABLES ARE SPECIFIED  R1 AND  
C                     QT  ( Q_TANGENTIAL)  AND   S~ R**GAMMA/P  (~ENTROPY)
C                     WITH ONE VAR. COMPUTED.  R2
C             FOR OUTFLOW : ONE VARIABLE IS FIXED R1  WITH THREE COMPUTED
C                     R2, QT AND S
C
C             COMPUTE RIEMANN INVARIANTS
C             FIX          R1 = QN - 2.*A/(GAMMA-1)  AT FREE STREAM
C             COMPUTE      R2 = QN + 2.*A/(GAMMA-1)  FROM INTERIOR
C                          EXTRAPOLATION OF FLOW VARIABLES
C
C             GET EXTRAPOLATED VARIABLES
C
      RHOEXT = Q(J,K-1,1)*XYJ(J,K-1)
      RJM1   = 1./Q(J,K-1,1)
      UEXT   = Q(J,K-1,2)*RJM1 
      VEXT   = Q(J,K-1,3)*RJM1 
      EEXT   = Q(J,K-1,4)*XYJ(J,K-1)
      PEXT   = GAMI*(EEXT - 0.5*RHOEXT*(UEXT**2+VEXT**2))
C
C               SET RIEMANN INVARIANTS
C
      R1 = XY3H*UF + XY4H*VF - 2.*AF*GM1I
      R2 = XY3H*UEXT + XY4H*VEXT + 2.*SQRT(GAMMA*PEXT/RHOEXT)*GM1I
C
      QN = (R1 + R2)*0.5
      CSPE = (R2 - R1)*GAMI*0.25
      C2 = CSPE**2
C
C               SET OTHER FIXED OR EXTRAPOLATED VARIABLES
C
           IF(QN .LE. 0.0)THEN
           QT = XY4H*UF - XY3H*VF 
           ENTRO = GAMMA
           ELSE
           QT = XY4H*UEXT - XY3H*VEXT 
           ENTRO = RHOEXT**GAMMA/PEXT
           ENDIF
C
C              COMPUTE FLOW VARIABLES
C
      U = XY3H*QN + XY4H*QT
      V = XY4H*QN - XY3H*QT
C
      Q(J,K,1) = (C2*ENTRO*GI)**GM1I
      PRES = C2*Q(J,K,1)*GI
C
C              ADD JACOBIAN
C
          RJJ = 1./XYJ(J,K)
          Q(J,K,1) = Q(J,K,1)*RJJ
          Q(J,K,2) = Q(J,K,1)*U
          Q(J,K,3) = Q(J,K,1)*V
          Q(J,K,4) = PRES*GM1I*RJJ + 0.5*Q(J,K,1)*(U**2+V**2)
60    CONTINUE
      ELSE
C
C                         SUPERSONIC FREESTREAM
C
      K = KMAX
      DO 70 J = 2,JM
C
C               CHOOSE A REFERENCE FRAME IN TERMS OF NORMAL AND 
C               TANGENTIAL COMPONENTS
C
C
C               METRIC TERMS
C  
C          USE EXTRAPOLATED QN TO DETERMINE INFLOW/OUTFLOW
C
      RHOINV = 1./Q(J,K-1,1)
      U = Q(J,K-1,2)*RHOINV
      V = Q(J,K-1,3)*RHOINV
      SNORM = 1./SQRT(XY(J,K-1,3)**2+XY(J,K-1,4)**2)
      XY3H = XY(J,K-1,3)*SNORM
      XY4H = XY(J,K-1,4)*SNORM
      QN1EXT = XY3H*U + XY4H*V
      RHOINV = 1./Q(J,K-2,1)
      U = Q(J,K-2,2)*RHOINV
      V = Q(J,K-2,3)*RHOINV
      SNORM = 1./SQRT(XY(J,K-2,3)**2+XY(J,K-2,4)**2)
      XY3H = XY(J,K-2,3)*SNORM
      XY4H = XY(J,K-2,4)*SNORM
      QN2EXT = XY3H*U + XY4H*V
      QNEXT = QN1EXT
C
           IF(QNEXT.GT.0.)THEN
           RMET1 = XYJ(J,K-1)/XYJ(J,K)
           RMET2 = XYJ(J,K-2)/XYJ(J,K)
           DO 65 N = 1,4
           Q(J,K,N) = Q(J,K-1,N)*RMET1
65         CONTINUE
           ENDIF
70    CONTINUE
      ENDIF

C
       RETURN
       END
