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 COORDINA
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 MOR
C   DETAIL.  CONTACT THOMAS H. PULLIAM, MS 202A-1 NASA AMES RESEARCH CEN
C   MOFFETT FIELD, CA.,  (415-694-6417) (EMAIL: PULLIAM@PRANDTL.NAS.NASA
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       CPUSEC:  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.6
C       DIS4Y:  FOURTH ORDER Y COEFFICIENT FOR FILTER       (DEFAULT=0.6
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 SCHEME
C      CAN BE EMPLOYED.  FOR STEADY STATE EULER IMPLICIT PERFORMS THE BE
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, E
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                                                                       
      PROGRAM ARC2D                                                     
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
                                                                        
      PARAMETER (MAXJ=289, MAXK=89, 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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       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 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                                                                       
      REAL T1,T2,T3,T4, T5, T6                                          
                                                                        
      COMMON/IO/IREAD,IWRIT,IVER                                        
CX00001 CALL MAINIX (1)                                                 
CY00001 CALL MAINIY (1)                                                 
CX60000 CALL CNTINX ()                                                  
CY60000 CALL CNTINY ()                                                  
C                                                                       
C     EXTERNAL SECOND                                                   
C     SECOND REPLACED BY CPUTIME AT CSRD JULY, 1989 BY L. POINTER       
C                                                                       
C                                                                       
      IREAD = 5                                                         
      IWRIT = 6                                                         
      IVER  = 7                                                         
C                                                                       
      OPEN(UNIT=IREAD,FILE='SRI5',FORM='FORMATTED',STATUS='OLD')        
      OPEN(UNIT=IWRIT,FILE='SRO6',FORM='FORMATTED',STATUS='NEW')        
      OPEN(UNIT=IVER, FILE='SRV', FORM='FORMATTED',STATUS='NEW')        
                                                                        
      CALL CPUTIM(T1)                                                   
      CALL ELAPSE(T2)                                                   
                                                                        
C  INPUT DATA                                                           
C                                                                       
       CALL INPUT(MAXJ,MAXJK)                                           
C                                                                       
C      INITIALIZE CPU TIME COUNTER                                      
C                                                                       
C         CPUSC0 = CPUTIM(DUMMY)                                        
          CALL CPUTIM(T5)                                               
C                                                                       
          WRITE(IWRIT,*) '--------------------------------------------' 
C         WRITE(IWRIT,*)'        CPUSEC AT START = ',CPUSC0             
          WRITE(IWRIT,*)'        CPUSEC AT START = ',T5                 
          WRITE(IWRIT,*) '--------------------------------------------' 
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                                                
CX10501 CALL DOINX( 501 )                                               
CY10501 CALL DOINY( 501 )                                               
         DO 10 NUMITE  = ISTART,IEND                                    
                                                                        
         CALL INTEGR(JDIMEN,KDIMEN,Q,S,PRESS,SNDSP,SLOPE,               
     >                 X,Y,XY,XYJ,XIT,ETT,DS,                           
     >                 UU,VV,CCX,CCY,COEF2X,COEF2Y,COEF4X,COEF4Y)       
                                                                        
C                                                                       
C  COMPUTE CPUSEC                                                       
C                                                                       
C        CPUSEC = CPUTIM(DUMMY) - CPUSC0                                
         CALL CPUTIM(T6)                                                
         CPUSEC = T6 -T5                                                
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*********************************************************************  
CX20502 CALL DOOUTX( 502 )                                              
CY20502 CALL DOOUTY( 502 )                                              
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 CPUTIM(T3)                                                
         CALL ELAPSE(T4)                                                
                                                                        
         WRITE(IVER,*) 'ELAPSED CPU TIME IN SECONDS: ',T3-T1            
         WRITE(IVER,*) 'MLOPS RATE: ', 1842.113398/(T3-T1)              
         WRITE(IVER,*) 'ELAPSED WALLCLOCK TIME IN SECONDS: ',           
     *                  MOD(T4-T2+1000000.0,1000000.0)                  
                                                                        
         CLOSE(IREAD)                                                   
         CLOSE(IWRIT)                                                   
         CLOSE(IVER)                                                    
CX70000 CALL CNTOTX ()                                                  
CY70000 CALL CNTOTY ()                                                  
CX50000 CALL MAINOX ()                                                  
CY50000 CALL MAINOY ()                                                  
                                                                        
         STOP                                                           
         END                                                            
C***********************************************************************
C****************** BICONVEX AIRFOIL BC             ********************
C***********************************************************************
      SUBROUTINE BC(JDIM,KDIM,Q,PRESS,SNDSP,XY,XIT,ETT,XYJ,X,Y,SLOPE)   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      PARAMETER (MAXJ=289, MAXK=89, 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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       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                          
      COMMON/WORKSP/F(MAXJ),WORK0(MAXJ,91)

CX30003 CALL SBINX (3)                                                  
CY30003 CALL SBINY (3)                                                  
C                                                                       
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                                
CX10503 CALL DOINX( 503 )                                               
CY10503 CALL DOINY( 503 )                                               
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                                                          
CX20504 CALL DOOUTX( 504 )                                              
CY20504 CALL DOOUTY( 504 )                                              
C                                                                       
C  SATISFY MOMENTUM RELATION FOR NORMAL PRESSURE DERIVATION             
C                                                                       
      K = 1                                                             
CX10505 CALL DOINX( 505 )                                               
CY10505 CALL DOINY( 505 )                                               
      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                                           
CX20506 CALL DOOUTX( 506 )                                              
CY20506 CALL DOOUTY( 506 )                                              
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                                                             
CX10507 CALL DOINX( 507 )                                               
CY10507 CALL DOINY( 507 )                                               
      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                                                          
CX20508 CALL DOOUTX( 508 )                                              
CY20508 CALL DOOUTY( 508 )                                              
C                                                                       
C                                                                       
C  OUTFLOW BOUNDARY                                                     
C                                                                       
C                                                                       
      J = JMAX                                                          
CX10509 CALL DOINX( 509 )                                               
CY10509 CALL DOINY( 509 )                                               
      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                                                          
CX20510 CALL DOOUTX( 510 )                                              
CY20510 CALL DOOUTY( 510 )                                              
                                                                        
      GO TO 100                                                         
C                                                                       
75    CONTINUE                                                          
C                                                                       
C  FOR SUPERSONIC INFLOW  ALL VARIABLES FIXED TO FREE STREAM            
C                                                                       
       J = 1                                                            
CX10511 CALL DOINX( 511 )                                               
CY10511 CALL DOINY( 511 )                                               
       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                                                         
CX20512 CALL DOOUTX( 512 )                                              
CY20512 CALL DOOUTY( 512 )                                              
C                                                                       
C   FOR SUPERSONIC OUTFLOW EXTRAPOLATE ALL VARIABLES                    
C                                                                       
       J = JMAX                                                         
CX10513 CALL DOINX( 513 )                                               
CY10513 CALL DOINY( 513 )                                               
       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                                                         
CX20514 CALL DOOUTX( 514 )                                              
CY20514 CALL DOOUTY( 514 )                                              
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                                          
CX10515 CALL DOINX( 515 )                                               
CY10515 CALL DOINY( 515 )                                               
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  (~ENTROP
C                     WITH ONE VAR. COMPUTED.  R2                       
C             FOR OUTFLOW : ONE VARIABLE IS FIXED R1  WITH THREE COMPUTE
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                                                          
CX20516 CALL DOOUTX( 516 )                                              
CY20516 CALL DOOUTY( 516 )                                              
      ELSE                                                              
C                                                                       
C                         SUPERSONIC FREESTREAM                         
C                                                                       
      K = KMAX                                                          
CX10517 CALL DOINX( 517 )                                               
CY10517 CALL DOINY( 517 )                                               
      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                                          
      QNEXT = QN1EXT                                                    
C                                                                       
           IF(QNEXT.GT.0.)THEN                                          
           RMET1 = XYJ(J,K-1)/XYJ(J,K)                                  
           DO 65 N = 1,4                                                
           Q(J,K,N) = Q(J,K-1,N)*RMET1                                  
65         CONTINUE                                                     
           ENDIF                                                        
70    CONTINUE                                                          
CX20518 CALL DOOUTX( 518 )                                              
CY20518 CALL DOOUTY( 518 )                                              
      ENDIF                                                             
CX40004 CALL SBOUTX (4)                                                 
CY40004 CALL SBOUTY (4)                                                 
                                                                        
C                                                                       
       RETURN                                                           
       END                                                              
                                                                        
C***********************************************************************
C****************** READ IN DATA AND SET PARAMETERS ********************
C***********************************************************************
      SUBROUTINE INPUT(JLARGE,JKLARG)                                   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
      COMMON/IO/IREAD,IWRIT,IVER                                        
CX30005 CALL SBINX (5)                                                  
CY30005 CALL SBINY (5)                                                  
C                                                                       
C                                                                       
C         SET PARAMETER DEFAULTS                                        
C                                                                       
      JMAX       = 65                                                   
      KMAX       = 21                                                   
      NP         = 10000                                                
      NPCP       = 1                                                    
      IPRINT     = 1                                                    
      JTAIL1     = 15                                                   
      JTAIL2     = 51                                                   
      JACDT      = 1                                                    
      THETAD     = 1.                                                   
      PHIDT      = 0.                                                   
      PERIDC     = .FALSE.                                              
      FSMACH     = 0.8                                                  
      ALPHA      = 0.                                                   
      DIS2X      = 0.0                                                  
      DIS4X      = 0.64                                                 
      DIS2Y      = 0.0                                                  
      DIS4Y      = 0.64                                                 
      DYM        = 0.02                                                 
      XMIN       = 2.0                                                  
      XMAX       = 3.0                                                  
      YMAX       = 3.0                                                  
      THICK      = 0.1                                                  
      NSTEPS     = 100                                                  
      DT         = 1.0                                                  
                                                                        
C                                                                       
      WRITE(IWRIT,*)'             ',                                    
     1   '       *         ******        ****     *****     *****  '    
      WRITE(IWRIT,*)'             ',                                    
     1   '      * *        *     *      *        *     *    *    * '    
      WRITE(IWRIT,*)'             ',                                    
     1   '     *   *       *      *    *              *     *     *'    
      WRITE(IWRIT,*)'             ',                                    
     1   '    *******      ******      *             *      *     *'    
      WRITE(IWRIT,*)'             ',                                    
     1   '   *       *     *     *      *          *        *    * '    
      WRITE(IWRIT,*)'             ',                                    
     1   '  *         *    *      *      ****     ******    *****  '    
      WRITE(IWRIT,*) ' '                                                
C                                                                       
      WRITE(IWRIT,*)'ARC2D : AN IMPLICIT FINITE ',                      
     1               'DIFFERENCE CODE FOR AIRFOILS'                     
      WRITE(IWRIT,*)'FOR FURTHER DETAILS SEE SELF ',                    
     1               'DOCUMENTATION OR CONTACT '                        
      WRITE(IWRIT,*)'THOMAS H. PULLIAM'                                 
      WRITE(IWRIT,*)'MS 202A-1, NASA AMES RESEARCH CENTER'              
      WRITE(IWRIT,*)'MOFFETT FIELD, CA., 94035'                         
      WRITE(IWRIT,*)'PHONE : 415 - 694 - 6417 '                         
C                                                                       
C               READ IN THE DATA FOR THE RUN FROM STANDARD INPUT        
C  *1*                                                                  
         WRITE(IWRIT,*) 'ENTER NSTEPS, DT DEFAULTS ARE ', NSTEPS,DT     
         READ(IREAD,*) NSTEPS,DT                                        
C  *2*                                                                  
         WRITE(IWRIT,*) 'ENTER JMAX, KMAX  DEFAULTS ARE ', JMAX, KMAX   
         READ(IREAD,*) JMAX, KMAX                                       
C  *3*                                                                  
         WRITE(IWRIT,*) 'ENTER FSMACH, ALPHA  DEFAULTS ARE ',           
     1                   FSMACH, ALPHA                                  
         READ(IREAD,*) FSMACH, ALPHA                                    
C  *4*                                                                  
         WRITE(IWRIT,*) 'ENTER DIS2X, DIS2X  DEFAULTS ARE ',            
     *                   DIS2X, DIS4X                                   
         READ(IREAD,*) DIS2X, DIS4X                                     
         WRITE(IWRIT,*) 'ENTER DIS2Y, DIS4Y  DEFAULTS ARE ',            
     *                   DIS2Y, DIS4Y                                   
         READ(IREAD,*) DIS2Y, DIS4Y                                     
C  *5*                                                                  
         WRITE(IWRIT,*) 'ENTER PERIDC    DEFAULT IS ',  PERIDC          
         READ(IREAD,*) PERIDC                                           
             IF(.NOT.PERIDC  )THEN                                      
                 WRITE(IWRIT,*) 'ENTER JTAIL1, JTAIL2  DEFAULTS ARE ',  
     >                           JTAIL1, JTAIL2                         
                 READ(IREAD,*) JTAIL1, JTAIL2                           
             ELSE                                                       
                 JTAIL1 = 1                                             
                 JTAIL2 = JMAX                                          
             ENDIF                                                      
C  *6*                                                                  
         WRITE(IWRIT,*) 'ENTER GRID SPECIFICATIONS: '                   
         WRITE(IWRIT,*) 'ENTER XMIN, XMAX, YMAX  DEFAULTS ARE ',        
     *                   XMIN,XMAX,YMAX                                 
         READ(IREAD,*) XMIN, XMAX, YMAX                                 
         WRITE(IWRIT,*)'ENTER DYM, THICK  DEFAULTS ARE ',DYM, THICK     
         READ(IREAD,*) DYM, THICK                                       
C  *7*                                                                  
         WRITE(IWRIT,*) 'ENTER PHIDT, THETAD   DEFAULTS ARE ',          
     *                   PHIDT, THETAD                                  
         READ(IREAD,*) PHIDT, THETAD                                    
C  *8*                                                                  
         WRITE(IWRIT,*) 'ENTER JACDT  DEFAULT  ', JACDT                 
         READ(IREAD,*) JACDT                                            
C  *9*                                                                  
         WRITE(IWRIT,*) 'ENTER IPRINT  DEFAULT  ', IPRINT               
         READ(IREAD,*) IPRINT                                           
C  *10*                                                                 
         WRITE(IWRIT,*) 'ENTER NP, NPCP  DEFAULTS ARE ',NP,'  ',NPCP    
         READ(IREAD,*) NP, NPCP                                         
C  *11*                                                                 
C                                                                       
C                                                                       
C  CHECK MESH SIZES                                                     
C                                                                       
      IF ( JMAX .GT. JLARGE ) THEN                                      
         WRITE(IWRIT, 36) JMAX, JLARGE                                  
36       FORMAT( 1X, 'JMAX IS ', I5, 5X, 'ONLY DIMENSIONED FOR ', I5)   
CX40006 CALL SBOUTX (6)                                                 
CY40006 CALL SBOUTY (6)                                                 
         STOP 'ERROR 1 IN IOALL'                                        
      ENDIF                                                             
      IF ( JMAX*KMAX .GT. JKLARG ) THEN                                 
         WRITE(IWRIT, 37) JMAX*KMAX, JKLARG                             
37       FORMAT( 1X,'JMAX*KMAX IS ',I5,5X,'ONLY DIMENSIONED FOR ',I5)   
CX40006 CALL SBOUTX (6)                                                 
CY40006 CALL SBOUTY (6)                                                 
         STOP 'ERROR 1 IN IOALL'                                        
      ENDIF                                                             
      IF ( JMAX .GT. 999 .OR. KMAX .GT. 999 ) THEN                      
         WRITE(IWRIT, 38) JMAX,KMAX                                     
38       FORMAT( 1X, 'JMAX,KMAX ARE ', 2I5, 5X,                         
     *    'BIGGER THAN 999 ', I5)                                       
CX40006 CALL SBOUTX (6)                                                 
CY40006 CALL SBOUTY (6)                                                 
         STOP 'ERROR 1 IN IOALL'                                        
      ENDIF                                                             
C                                                                       
      WRITE(IWRIT,*) ' '                                                
C                                                                       
      WRITE(IWRIT,*)'    '                                              
      WRITE(IWRIT,*)'...........FLOW FIELD CONDITIONS...............'   
      WRITE(IWRIT,*)' FSMACH = ',FSMACH,' ALPHA = ',ALPHA               
C                                                                       
      WRITE(IWRIT,*)'    '                                              
      WRITE(IWRIT,*)'.............GRID SPECIFICATIONS...............'   
      WRITE(IWRIT,*)' BICONVEX AIRFOIL GRID ',                          
     1               'IS BEING INTERNALLY GENERATED'                    
      WRITE(IWRIT,*)'...............................................'   
      WRITE(IWRIT,*)'  MINIMUM SPACING AT BODY'                         
      WRITE(IWRIT,*)'                DYM = ',DYM                        
      WRITE(IWRIT,*)' UPPER DIMENSION IN CHORDS (C=1) '                 
      WRITE(IWRIT,*)'                 YMAX = ',YMAX                     
      WRITE(IWRIT,*)' LEFT XMIN = ',XMIN,'  RIGHT  XMAX = ',XMAX        
      WRITE(IWRIT,*)' THICKNESS OF BICONVEX THICK = ',THICK             
      WRITE(IWRIT,*)'...............................................'   
                                                                        
      IF(PERIDC  )THEN                                                  
        WRITE(IWRIT,*)' THIS IS A PERIDC   GRID',                       
     1   ' JMAX = ',JMAX,' KMAX = ',KMAX                                
      ELSE                                                              
        WRITE(IWRIT,*)' THIS IS A NONPERIDC   GRID',                    
     1   ' JMAX = ',JMAX,' KMAX = ',KMAX,                               
     1   ' JTAIL1 = ',JTAIL1,' JTAIL2 = ',JTAIL2                        
      ENDIF                                                             
C                                                                       
C                                                                       
      WRITE(IWRIT,*)'    '                                              
      WRITE(IWRIT,*)'.............I/O SPECIFICATIONS.............'      
      WRITE(IWRIT,*)' NP = ',NP,'  IPRINT = ',IPRINT                    
      WRITE(IWRIT,*)'.............ALGORITHM PARAMETERS.............'    
      WRITE(IWRIT,*)' THETAD  = ',THETAD ,' PHIDT = ',PHIDT             
      WRITE(IWRIT,*)' JACDT = ',JACDT                                   
C                                                                       
C                                                                       
      IF(JACDT.EQ.0)WRITE(IWRIT,*)                                      
     1     ' <<<<<<<<<<<<  TIME ACCURATE  >>>>>>>>>>>>>>>>>>'           
      IF(JACDT.EQ.1)WRITE(IWRIT,*)                                      
     1     ' <<<<  VARIABLE DT = 1./(1+SQRT(J)) >>>>>>>>>>>>>'          
      IF(JACDT.EQ.2)WRITE(IWRIT,*)                                      
     1       ' <<<<  VARIABLE DT : MIN. EIGEN. >>>>>>>>>>>>>>>'         
      IF(JACDT.EQ.3)WRITE(IWRIT,*)                                      
     1      ' <<<<  VARIABLE DT : CONST. CFL  >>>>>>>>>>>>>>>'          
      WRITE(IWRIT,*)'   ********  NUMBER OF TIME STEPS == ',NSTEPS      
      WRITE(IWRIT,*)'   ********  TIME STEP            == ',DT          
C                                                                       
      WRITE(IWRIT,*)'*************************************************' 
      WRITE(IWRIT,*)'......   PENTADIAGONAL ALGORITHM      ...........' 
      WRITE(IWRIT,*)'*************************************************' 
      WRITE(IWRIT,*)' DIS2X = ',DIS2X,' DIS4X = ',DIS4X                 
      WRITE(IWRIT,*)' DIS2Y = ',DIS2Y,' DIS4Y = ',DIS4Y                 
CX40006 CALL SBOUTX (6)                                                 
CY40006 CALL SBOUTY (6)                                                 
                                                                        
C                                                                       
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C****************** INITIALIZE VARIABLES AND DATA **********************
C***********************************************************************
      SUBROUTINE INITIA(JDIM,KDIM,Q,PRESS,SNDSP,S,XY,XYJ,               
     *                  XIT,ETT,DS,X,Y,SLOPE)                           
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       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 S(JDIM,KDIM,4),XY(JDIM,KDIM,4),XYJ(JDIM,KDIM)           
      DIMENSION XIT(JDIM,KDIM),ETT(JDIM,KDIM),DS(JDIM,KDIM)             
      DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM),SLOPE(JDIM)                   
      COMMON/INFNTY/RHOINF,UINF,VINF,EINF,PINF                          
CX30007 CALL SBINX (7)                                                  
CY30007 CALL SBINY (7)                                                  
                                                                        
C                                                                       
C                                                                       
C               COMPUTE VARIOUS CONSTANTS                               
C                                                                       
      GAMMA = 1.4                                                       
      GAMI = GAMMA -1.                                                  
      PI = 4.*ATAN(1.0)                                                 
C                                                                       
      COSANG = COS( PI*ALPHA/180.)                                      
      SINANG = SIN(PI*ALPHA/180.)                                       
      RHOINF = 1.                                                       
      UINF   = FSMACH*COSANG                                            
      VINF   = FSMACH*SINANG                                            
      EINF   = 1./(GAMMA*GAMI) + .5*FSMACH**2                           
      PINF   = 1./GAMMA                                                 
C                                                                       
C              GRID CALCULATION                                         
C                                                                       
      CALL BICONG(JDIM,KDIM,JMAX,KMAX,JTAIL1,JTAIL2,DYM,THICK,          
     1 YMAX,XMAX,XMIN,X,Y)                                              
CX10519 CALL DOINX( 519 )                                               
CY10519 CALL DOINY( 519 )                                               
C                                                                       
C                INITIALIZE TO FREE STREAM                              
C                                                                       
C                LOAD FREE STREAM AT ANGLE OF ATTACK ALPHA              
C                                                                       
      DO 32 K=1,KMAX                                                    
      DO 32 J=1,JMAX                                                    
         Q(J,K,1) = RHOINF                                              
         Q(J,K,2) = RHOINF*UINF                                         
         Q(J,K,3) = RHOINF*VINF                                         
         Q(J,K,4) = EINF                                                
   32 CONTINUE                                                          
CX20520 CALL DOOUTX( 520 )                                              
CY20520 CALL DOOUTY( 520 )                                              
C                                                                       
C                                                                       
C  SET UP VARIABLES, INDECIES, AND COMPUTE METRICS                      
C                                                                       
C                                                                       
C  JTAIL2 = JMAX  FOR PERIDC   MESHES                                   
C                                                                       
      IF(PERIDC  )JTAIL2 = JMAX                                         
C                                                                       
C            SET PERIDC   INDEX                                         
C                                                                       
      JM = JMAX-1                                                       
      KM = KMAX-1                                                       
      JBEGIN = 1                                                        
      JEND = JMAX                                                       
      KBEGIN = 1                                                        
      KEND = KMAX                                                       
      KLOW = 2                                                          
      KUP = KM                                                          
CX10521 CALL DOINX( 521 )                                               
CY10521 CALL DOINY( 521 )                                               
      DO 10 J = 1, JMAX                                                 
         JPLUS(J) = J+1                                                 
         JMINU (J) = J-1                                                
10    CONTINUE                                                          
CX20522 CALL DOOUTX( 522 )                                              
CY20522 CALL DOOUTY( 522 )                                              
C                                                                       
      IF ( .NOT.PERIDC    ) THEN                                        
         JPLUS(JMAX) = JMAX                                             
         JMINU (1) = 1                                                  
         JLOW = 2                                                       
         JUP = JMAX-1                                                   
      ELSE                                                              
         JPLUS(JMAX) = 1                                                
         JMINU (1) = JMAX                                               
         JLOW = 1                                                       
         JUP = JMAX                                                     
         JTAIL1 = 1                                                     
         JTAIL2 = JMAX                                                  
      ENDIF                                                             
CX10523 CALL DOINX( 523 )                                               
CY10523 CALL DOINY( 523 )                                               
C                                                                       
      DO 8 K=1,KMAX                                                     
      DO 8 J=1,JMAX                                                     
         XIT(J,K) = 0.                                                  
         ETT(J,K) = 0.                                                  
8     CONTINUE                                                          
CX20524 CALL DOOUTX( 524 )                                              
CY20524 CALL DOOUTY( 524 )                                              
C                                                                       
C                                                                       
C                METRIC CALCULATION                                     
C   S USED FOR WORK SPACE HERE                                          
C                                                                       
      CALL XYMETS(JDIM,KDIM,X,Y,XY,XIT,ETT,XYJ,S)                       
CX10525 CALL DOINX( 525 )                                               
CY10525 CALL DOINY( 525 )                                               
C                                                                       
      DO 9 N=1,4                                                        
      DO 9 K=1,KMAX                                                     
      DO 9 J=1,JMAX                                                     
         S(J,K,N) = 0.                                                  
9     CONTINUE                                                          
CX20526 CALL DOOUTX( 526 )                                              
CY20526 CALL DOOUTY( 526 )                                              
C                                                                       
C               INITIALIZE VARIABLE TIME STEP IF JACDT=1                
C                                                                       
      FACTOR=1.0                                                        
      IF(JACDT.EQ.0) FACTOR=0.0                                         
CX10527 CALL DOINX( 527 )                                               
CY10527 CALL DOINY( 527 )                                               
      DO 65 K=1,KMAX                                                    
      DO 65 J=1,JMAX                                                    
         DS(J,K)=1.0/(1.0+FACTOR*SQRT(XYJ(J,K)))                        
   65 CONTINUE                                                          
CX20528 CALL DOOUTX( 528 )                                              
CY20528 CALL DOOUTY( 528 )                                              
CX10529 CALL DOINX( 529 )                                               
CY10529 CALL DOINY( 529 )                                               
C                                                                       
C                                                                       
C               SCALE VARIABLES BY JACOBIAN STORED BY XYMETS IN XYJ(J,K)
C                                                                       
      DO 50 N=1,4                                                       
      DO 50 K=1,KMAX                                                    
      DO 50 J=1,JMAX                                                    
       Q(J,K,N) = Q(J,K,N)/XYJ(J,K)                                     
   50 CONTINUE                                                          
CX20530 CALL DOOUTX( 530 )                                              
CY20530 CALL DOOUTY( 530 )                                              
CX10531 CALL DOINX( 531 )                                               
CY10531 CALL DOINY( 531 )                                               
C                                                                       
C                                                                       
C  AIRFOIL BODY SLOPE                                                   
C                                                                       
C    USED ONLY IF RUNNING THE BICONVEX AIRFOIL                          
C                                                                       
      DO 22 J = 1,JMAX                                                  
         SLOPE(J) = 0.0                                                 
22    CONTINUE                                                          
CX20532 CALL DOOUTX( 532 )                                              
CY20532 CALL DOOUTY( 532 )                                              
CX10533 CALL DOINX( 533 )                                               
CY10533 CALL DOINY( 533 )                                               
C                                                                       
      DO 24 J = JTAIL1,JTAIL2                                           
         XPT = X( J, 1)                                                 
         DGDX = 2. * THICK * ( 1. - 2. * XPT )                          
         SLOPE( J) = ATAN( DGDX)                                        
24    CONTINUE                                                          
CX20534 CALL DOOUTX( 534 )                                              
CY20534 CALL DOOUTY( 534 )                                              
CX10535 CALL DOINX( 535 )                                               
CY10535 CALL DOINY( 535 )                                               
      DO 26 J = 1,JMAX                                                  
                                                                        
  25     FORMAT(' ','****J,SLOPE',I5,2X,E20.8)                          
  26  CONTINUE                                                          
CX20536 CALL DOOUTX( 536 )                                              
CY20536 CALL DOOUTY( 536 )                                              
C                                                                       
C   INITIALIZE PRESS,SNDSP                                              
C                                                                       
      CALL CALCPS(JDIM,KDIM,Q,PRESS,SNDSP)                              
CX40008 CALL SBOUTX (8)                                                 
CY40008 CALL SBOUTY (8)                                                 
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C***************** INPUT/OUTPUT ROUTINE ********************************
C***********************************************************************
      SUBROUTINE IOALL(IWHAT,JDIM,KDIM,Q,PRESS,SNDSP,XY,XYJ,X,Y)        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
C                                                                       
      PARAMETER (MAXJ=289, MAXK=89, 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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       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 X(JDIM,KDIM),Y(JDIM,KDIM)                               
C                                                                       
      COMMON/WORKSP/CPBOD(MAXJ),XBOD(MAXJ), WORK2(MAXJ,90)              
C                                                                       
      COMMON/IO/IREAD,IWRIT,IVER                                        
CX30009 CALL SBINX (9)                                                  
CY30009 CALL SBINY (9)                                                  
C                                                                       
      GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17) IWHAT           
C                                                                       
    1 CONTINUE                                                          
C                                                                       
      GO TO 10000                                                       
C                                                                       
C                                                                       
C***********************************************************************
C                                                                       
    2 CONTINUE                                                          
      GO TO 10000                                                       
C                                                                       
C********************************************************************** 
C                                                                       
C                                                                       
    3 CONTINUE                                                          
C                                                                       
C                                                                       
      GO TO 10000                                                       
C                                                                       
C                                                                       
C***********************************************************************
C          READ THE RESTART                                             
C                                                                       
    4 CONTINUE                                                          
C                                                                       
C      READ IN A UNFORMATTED RESTART FILE FROM TAPE3.                   
C      THE Q READ IN IS NOT YET SCALED BY THE METRIC JACOBIAN XYJ,      
C      THE RESCALING IS DONE IN INITIA.                                 
C                                                                       
C  CHANGED TO BE COMPATIBLE WITH  NEW PLOT3D  4/2/86                    
C                                                                       
      OPEN(UNIT=3, FILE='SRI3',STATUS='OLD', FORM='UNFORMATTED')        
      READ(3) JMAX2,KMAX2                                               
      READ(3) FSMTEM,ALPTEM,RETEM,TOTIME                                
C                                                                       
      WRITE(IWRIT, 433) TOTIME, ISTART                                  
  433 FORMAT('0','TOTAL TIME =',F12.6,                                  
     1       10X, 'STARTING ITERATION =', I5)                           
C                                                                       
      READ(3) ( ( ( Q(J,K,N),J=1,JMAX2) ,K=1,KMAX2) ,N=1,4)             
      READ(3) JTITEM, ISTART                                            
C                                                                       
      GO TO 10000                                                       
C                                                                       
C                                                                       
C***********************************************************************
C          ITERATION NUMBER AND RESIDUAL WRITTEN TO OUTPUT              
C                                                                       
    5 CONTINUE                                                          
C                                                                       
C                                                                       
C  COMPUTE NUMBER OF SUPERSONIC POINTS                                  
C                                                                       
      NSUPER = 0                                                        
CX10537 CALL DOINX( 537 )                                               
CY10537 CALL DOINY( 537 )                                               
      DO 501 K = 1,KMAX                                                 
          DO 502 J = JLOW,JUP                                           
          RI = 1./Q(J,K,1)                                              
          UV2 = (Q(J,K,2)**2 + Q(J,K,3)**2)*RI**2                       
C            CHECK MACH**2  ( SAME AS MACH)                             
          CSPED2 = SNDSP(J,K)**2                                        
          ZS = UV2/CSPED2                                               
             IF( ZS .GE. 1.)NSUPER = NSUPER + 1                         
502       CONTINUE                                                      
501   CONTINUE                                                          
CX20538 CALL DOOUTX( 538 )                                              
CY20538 CALL DOOUTY( 538 )                                              
C                                                                       
C          LIFT AND DRAG CALCULATIONS AND OUTPUT ON TAPE6               
C                                                                       
      RSD = RESID                                                       
C                                                                       
C  OUTPUT                                                               
C                                                                       
C                                                                       
      IF (NUMITE.EQ.100) THEN                                           
         CALL VERIFY(NUMITE,RSD)                                        
      ENDIF                                                             
                                                                        
      WRITE(IWRIT, 505)NUMITE , RSD, NSUPER, CPUSEC                     
505   FORMAT(' ',' N = ',I6, '  RESID = ',E12.6,                        
     >           '  # OF SS PTS = ',I5,'  CPUSEC = ',F10.2)             
C                                                                       
      GO TO 10000                                                       
C                                                                       
C                                                                       
C***********************************************************************
C          COEFFICIENT OF PRESSURE DATA WRITTEN TO OUTPUT               
C                                                                       
    6 CONTINUE                                                          
C                                                                       
C                                                                       
      WRITE(IWRIT, 666)                                                 
666   FORMAT( 1X, 'COEFFICIENT OF PRESSURE ')                           
      CPC = 2./(GAMMA*FSMACH**2)                                        
C   COMPUTE CP AT GRID POINTS AND STORE IN Z ARRAY                      
CX10539 CALL DOINX( 539 )                                               
CY10539 CALL DOINY( 539 )                                               
      DO 665 J=JTAIL1,JTAIL2                                            
         PP = PRESS(J,1)*XYJ(J,1)                                       
         CP = (PP*GAMMA -1.)*CPC                                        
         WRITE(IWRIT, 667) J, X(J,1), CP                                
665   CONTINUE                                                          
CX20540 CALL DOOUTX( 540 )                                              
CY20540 CALL DOOUTY( 540 )                                              
C                                                                       
667   FORMAT( 1X, I5, 5X, F14.7, 5X, F14.7)                             
C                                                                       
C                                                                       
      GO TO 10000                                                       
C                                                                       
C***********************************************************************
   7    CONTINUE                                                        
      GO TO 10000                                                       
C                                                                       
C***********************************************************************
C                                                                       
    8 CONTINUE                                                          
C                                                                       
      GO TO 10000                                                       
C                                                                       
C                                                                       
C***********************************************************************
C                                                                       
    9 CONTINUE                                                          
C                                                                       
      GO TO 10000                                                       
C                                                                       
C                                                                       
C***********************************************************************
C          STORE A UNFORMATTED RESTART FILE                             
C                                                                       
   10 CONTINUE                                                          
C                                                                       
C     NOTE THAT NUMITE  IS ACTUALLY ONE GREATER THAN THE LAST USED.     
C     IT GETS STORED, AND USED AS THE FIRST ITERATION NEXT TIME.        
C                                                                       
C       REWIND BECAUSE THIS MAY BE CALLED PERIDC  ALLY                  
      OPEN(UNIT=4, FILE='SRO4',STATUS='NEW', FORM='UNFORMATTED')        
      WRITE(4) JMAX,KMAX                                                
      WRITE(4) FSMACH,ALPHA,FSMACH,ALPHA                                
      WRITE(4) ( ( ( Q(J,K,N)*XYJ(J,K),J=1,JMAX) ,K=1,KMAX) ,N=1,4)     
      WRITE(4) JTAIL1, NUMITE                                           
      REWIND 4                                                          
C                                                                       
      GO TO 10000                                                       
C                                                                       
C***********************************************************************
C                                                                       
   11 CONTINUE                                                          
      GOTO 10000                                                        
C***********************************************************************
C                                                                       
   12 CONTINUE                                                          
      GOTO 10000                                                        
C***********************************************************************
   13 CONTINUE                                                          
      GOTO 10000                                                        
C***********************************************************************
C***********************************************************************
   14 CONTINUE                                                          
      GOTO 10000                                                        
C***********************************************************************
C***********************************************************************
C                                                                       
   15 CONTINUE                                                          
C                                                                       
C                                                                       
C    OUTPUT ON TAPE6 9 VARIABLES ALONG THE BODY                         
C                                                                       
      K = 1                                                             
      WRITE(IWRIT,  1510) K                                             
1510  FORMAT( 1X, 'K = , I4 ')                                          
      WRITE(IWRIT, 1515)                                                
1515  FORMAT( 1X,                                                       
     1   '  J   K     RHO            U             V           ',       
     2   'ENERGY      PRESSURE     ',                                   
     3   ' COEFF P    LOCAL MACH     X             Y ')                 
CX10541 CALL DOINX( 541 )                                               
CY10541 CALL DOINY( 541 )                                               
C                                                                       
      DO 1590 J = 1, JMAX                                               
         RHO = Q( J, K, 1) * XYJ( J, K)                                 
         U = Q( J, K, 2) / Q( J, K, 1)                                  
         V = Q( J, K, 3) / Q( J, K, 1)                                  
         UVSQ = U**2 + V**2                                             
         ENERGY = Q( J, K, 4) * XYJ( J, K)                              
         PPP = PRESS(J,K)*XYJ(J,K)                                      
         CPC = 2./(GAMMA*FSMACH**2)                                     
         CP = (PPP * GAMMA - 1. ) * CPC                                 
         CSQ = SNDSP(J,K)                                               
         EMACH = SQRT( UVSQ )/ CSQ                                      
         XCOORD = X( J, K)                                              
         YCOORD = Y( J, K)                                              
         WRITE(IWRIT, 1520) J, K, RHO, U, V, ENERGY, PPP,               
     1      CP, EMACH, XCOORD, YCOORD                                   
1520     FORMAT( 1X, 2( I3, 1X), 9( F11.6, 2X ) )                       
1590   CONTINUE                                                         
CX20542 CALL DOOUTX( 542 )                                              
CY20542 CALL DOOUTY( 542 )                                              
C                                                                       
       GOTO 10000                                                       
C                                                                       
C***********************************************************************
C***********************************************************************
C          GRID OUTPUT                                                  
C                                                                       
   16 CONTINUE                                                          
C                                                                       
      OPEN(UNIT=11, FILE='SRO11',STATUS='NEW', FORM='UNFORMATTED')      
      WRITE( 11) JMAX,KMAX                                              
      WRITE( 11)  ( ( X(J,K),J=1,JMAX) ,K=1,KMAX) ,                     
     &             ( ( Y(J,K),J=1,JMAX) ,K=1,KMAX)                      
C                                                                       
C                                                                       
      GO TO 10000                                                       
C                                                                       
C                                                                       
C***********************************************************************
C***********************************************************************
C          COMPUTE CIRCULATION, AND PRINT IT OUT ( NOT IN COMMON )      
C                                                                       
   17 CONTINUE                                                          
C                                                                       
C                                                                       
      GO TO 10000                                                       
C                                                                       
C                                                                       
C***********************************************************************
C                                                                       
C                                                                       
10000 CONTINUE                                                          
CX40010 CALL SBOUTX (10)                                                
CY40010 CALL SBOUTY (10)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
C*************  CALCULATE PRESSURE AND SOUND SPEED *********************
C***********************************************************************
      SUBROUTINE CALCPS(JDIM,KDIM,Q,PRESS,SNDSP)                        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION Q(JDIM,KDIM,4),PRESS(JDIM,KDIM),SNDSP(JDIM,KDIM)        
CX30011 CALL SBINX (11)                                                 
CY30011 CALL SBINY (11)                                                 
CX10543 CALL DOINX( 543 )                                               
CY10543 CALL DOINY( 543 )                                               
C                                                                       
C  CALCULATE PRESSURE AND SOUND SPEED AND STORE                         
C                                                                       
C  NOTE PRESSURE HAS JACOBIAN                                           
C       SOUND SPEED DOES NOT HAVE JACOBIAN                              
C       CALL AFTER BC ROUTINE SO THAT BOUNDARY VALUES ARE CORRECT       
C                                                                       
        DO 1 K = KBEGIN,KEND                                            
        DO 1 J = JBEGIN,JEND                                            
C                                                                       
        PRESS(J,K) = GAMI*(Q(J,K,4) -                                   
     *               0.5*(Q(J,K,2)**2 + Q(J,K,3)**2)/Q(J,K,1))          
        SNDSP(J,K) = SQRT(GAMMA*PRESS(J,K)/Q(J,K,1))                    
C                                                                       
1       CONTINUE                                                        
CX20544 CALL DOOUTX( 544 )                                              
CY20544 CALL DOOUTY( 544 )                                              
CX40012 CALL SBOUTX (12)                                                
CY40012 CALL SBOUTY (12)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C***********************************************************************
C********************** METRIC CALCULATIONS ****************************
C***********************************************************************
      SUBROUTINE XYMETS(JDIM,KDIM,X,Y,XY,XIT,ETT,XYJ,WORK)              
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION XY(JDIM*KDIM*4),XYJ(JDIM*KDIM)                          
      DIMENSION XIT(JDIM*KDIM),ETT(JDIM*KDIM)                           
      DIMENSION X(JDIM*KDIM),Y(JDIM*KDIM)                               
      DIMENSION WORK(JDIM*KDIM*4)                                       
CX30013 CALL SBINX (13)                                                 
CY30013 CALL SBINY (13)                                                 
C                                                                       
C                                                                       
C    COMPUTE XI DERIVATIVES OF X, Y                                     
C                                                                       
       CALL XIDIF(JDIM,KDIM,X,Y,XY)                                     
C                                                                       
C   COMPUTE  ETA DERIVATIVES OF X, Y                                    
C                                                                       
       CALL ETADIF(JDIM,KDIM,X,Y,XY)                                    
C                                                                       
C   FORM METRICS AND JACOBIAN                                           
C                                                                       
       CALL CALCME(JDIM,KDIM,XY,XIT,ETT,XYJ)                            
CX40014 CALL SBOUTX (14)                                                
CY40014 CALL SBOUTY (14)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C********************** METRIC CALCULATIONS ****************************
C***********************************************************************
      SUBROUTINE CALCME(JDIM,KDIM,XY,XIT,ETT,XYJ)                       
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
      COMMON/IO/IREAD,IWRIT,IVER                                        
C                                                                       
      DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM)                          
      DIMENSION XIT(JDIM,KDIM),ETT(JDIM,KDIM)                           
CX30015 CALL SBINX (15)                                                 
CY30015 CALL SBINY (15)                                                 
CX10545 CALL DOINX( 545 )                                               
CY10545 CALL DOINY( 545 )                                               
C                                                                       
C  FORM JACOBIAN                                                        
C                                                                       
      DO 30 K=KBEGIN,KEND                                               
      DO 30 J=JBEGIN,JEND                                               
         DINV = 1./ ( XY(J,K,4) * XY(J,K,1) - XY(J,K,3) * XY(J,K,2) )   
         IF(DINV .LE. 0.0)WRITE(IWRIT,*)'JACOBIAN( ',J,K,' )= ',DINV    
         XYJ(J,K) = DINV                                                
   30 CONTINUE                                                          
CX20546 CALL DOOUTX( 546 )                                              
CY20546 CALL DOOUTY( 546 )                                              
CX10547 CALL DOINX( 547 )                                               
CY10547 CALL DOINY( 547 )                                               
C                                                                       
C  CALCULATE METRICS                                                    
C                                                                       
      DO 32 K=KBEGIN,KEND                                               
      DO 32 J=JBEGIN,JEND                                               
         DINV = XYJ(J,K)                                                
C                                                                       
C XIX = XY1, XIY=XY2,ETAX=XY3, ETAY = XY4                               
C                                                                       
         XY(J,K,1) =   XY(J,K,1)*DINV                                   
         XY(J,K,2) = - XY(J,K,2)*DINV                                   
         XY(J,K,3) = - XY(J,K,3)*DINV                                   
         XY(J,K,4) =   XY(J,K,4)*DINV                                   
32    CONTINUE                                                          
CX20548 CALL DOOUTX( 548 )                                              
CY20548 CALL DOOUTY( 548 )                                              
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
C                                                                       
      RETURN                                                            
      END                                                               
                                                                        
C********************************************************************** 
C********************** DIFFERENCE X, Y  IN XI *************************
C***********************************************************************
      SUBROUTINE XIDIF(JDIM,KDIM,X,Y,XY)                                
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
C                                                                       
      COMMON/BASE/                                                      
     1  JMAX,       KMAX,     JM,          KM,          JBEGIN,         
     1  JEND,                                                           
     1  KBEGIN,     KEND,     JPLUS(999),  JMINU (999), JLOW,           
     1  JUP,                                                            
     1  KLOW,       KUP,      PERIDC  ,    NP,          DT,             
     1  CPUSEC,                                                         
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION XY(JDIM,KDIM,4)                                         
      DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM)                               
CX30017 CALL SBINX (17)                                                 
CY30017 CALL SBINY (17)                                                 
CX10549 CALL DOINX( 549 )                                               
CY10549 CALL DOINY( 549 )                                               
C                                                                       
C                                                                       
C  XY4 =XXI,XY3 = YXI                                                   
C                                                                       
      DO 11 K=KBEGIN,KEND                                               
         DO 10 J=JLOW,JUP                                               
            JP1 = JPLUS(J)                                              
            JM1 = JMINU (J)                                             
            XY(J,K,4) = ( X(JP1,K) - X(JM1,K))*.5                       
            XY(J,K,3) = ( Y(JP1,K) - Y(JM1,K))*.5                       
10       CONTINUE                                                       
C                                                                       
          IF(.NOT.PERIDC  )THEN                                         
            J = JBEGIN                                                  
            XY(J,K,4) = .5*( -3.*X(J,K) +4.*X(J+1,K) - X(J+2,K))        
CX20550 CALL DOOUTX( 550 )                                              
CY20550 CALL DOOUTY( 550 )                                              
            XY(J,K,3) = .5*( -3.*Y(J,K) +4.*Y(J+1,K) - Y(J+2,K))        
            J = JEND                                                    
            XY(J,K,4) = ( 3.*X(J,K) -4.*X(J-1,K) + X(J-2,K))*.5         
            XY(J,K,3) = ( 3.*Y(J,K) -4.*Y(J-1,K) + Y(J-2,K))*.5         
          ENDIF                                                         
C                                                                       
   11 CONTINUE                                                          
CX40018 CALL SBOUTX (18)                                                
CY40018 CALL SBOUTY (18)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C***********************************************************************
C********************** DIFFERENCE X, Y  IN XI *************************
C***********************************************************************
      SUBROUTINE ETADIF(JDIM,KDIM,X,Y,XY)                               
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION XY(JDIM,KDIM,4)                                         
      DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM)                               
CX30019 CALL SBINX (19)                                                 
CY30019 CALL SBINY (19)                                                 
CX10551 CALL DOINX( 551 )                                               
CY10551 CALL DOINY( 551 )                                               
C                                                                       
C  XY2 = XETA, XY1 = YETA                                               
C                                                                       
      DO 21 J=JBEGIN,JEND                                               
         DO 20 K=KLOW,KUP                                               
            XY(J,K,2) = ( X(J,K+1) - X(J,K-1))*.5                       
            XY(J,K,1) = ( Y(J,K+1) - Y(J,K-1))*.5                       
20       CONTINUE                                                       
         K = KBEGIN                                                     
         XY(J,K,2) = ( -3.*X(J,K) +4.*X(J,K+1) - X(J,K+2))*.5           
         XY(J,K,1) = ( -3.*Y(J,K) +4.*Y(J,K+1) -Y(J,K+2))*.5            
         K = KEND                                                       
         XY(J,K,2) = ( 3.*X(J,K) -4.*X(J,K-1) + X(J,K-2))*.5            
         XY(J,K,1) = ( 3.*Y(J,K) -4.*Y(J,K-1) + Y(J,K-2))*.5            
   21 CONTINUE                                                          
CX20552 CALL DOOUTX( 552 )                                              
CY20552 CALL DOOUTY( 552 )                                              
CX40020 CALL SBOUTX (20)                                                
CY40020 CALL SBOUTY (20)                                                
C                                                                       
C                                                                       
         RETURN                                                         
         END                                                            
C***********************************************************************
C******************** DEFAULT GRID FOR *********************************
C**************** BICONVEX AIFOIL CALCULATIONS *************************
C***********************************************************************
      SUBROUTINE BICONG(JDIM,KDIM,NX,NY,J1,J2,DY,THICK1,TOP,RIGHT,ALEFT 
     >                   ,X,Y)                                          
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      PARAMETER (MAXJ=289, MAXK=89, MAXJK=MAXJ*MAXK  )                  
      DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM)                               
      COMMON/WORKSP/CLUS(MAXJ), WORK2(MAXJ,91)                          
CX30021 CALL SBINX (21)                                                 
CY30021 CALL SBINY (21)                                                 
C                                                                       
C                                                                       
C      J      REFERS TO THE X ( STREAMWISE ) DIRECTION                  
C      K      REFERS TO THE Y DIRECTION                                 
C                                                                       
C                                                                       
C POINTS IN X DIRECTION                                                 
      NJ = NX                                                           
C POINTS IN Y DIRECTION                                                 
      NK = NY                                                           
C                                                                       
C THICKNESS OF AIRFOIL ( 10% PARABOLIC BICONVEX )                       
C      THICK = THICK1                                                   
C                                                                       
C      BELOW, REGION IS [ -X1, X2 ] X [ 0, Y1 ]                         
C      AND AIRFOIL IS [ 0, 1] ON THE X AXIS                             
C -X1 IS LEFT BOUNDARY OF REGION                                        
C      X1 = -ALEFT                                                      
C X2 IS LENGTH OF REGION TO THE RIGHT OF THE AIRFOIL                    
C      X2 = RIGHT - 1.                                                  
C TOP OF GRID, NUMBER OF CHORD LENGTHS PERP. TO                         
C STREAMWISE DIRECTION                                                  
C      Y1 = TOP                                                         
C                                                                       
C        DEFINES THE FINE MESH                                          
C                                                                       
C      NJM1=NJ-1                                                        
C      NKM1=NK-1                                                        
C      NJM2=NJ-2                                                        
C      NKM2=NK-2                                                        
C                                                                       
C     NUMBER OF POINTS ON AIRFOIL                                       
C                                                                       
      JBODY=J2-J1+1                                                     
      PI=3.141592653589793                                              
C                                                                       
C     PUT THE POINTS ON THE AIRFOIL                                     
C                                                                       
      DANG=PI/JBODY                                                     
CX10553 CALL DOINX( 553 )                                               
CY10553 CALL DOINY( 553 )                                               
C                                                                       
C      TWOTH=2.*THICK                                                   
      DO 1 J=J1,J2                                                      
      XI=J-J1+.5                                                        
C THIS USES COS ON [ 0, PI ]                                            
      DEV=COS(XI*DANG)                                                  
C      X(J,1)=.5*(1.-DEV)                                               
CJB                                                                     
CJB      THIS ONE GIVES A LOT OF STRETCHING                             
CJB            X(1,J)=.5*(1.-DEV)                                       
CJB                                                                     
C        THIS GIVES LESS STRETCHING                                     
              X( J, 1) = .5 * ( .5*(1.-DEV) + XI/JBODY )                
C                                                                       
    1 CONTINUE                                                          
CX20554 CALL DOOUTX( 554 )                                              
CY20554 CALL DOOUTY( 554 )                                              
CJB                                                                     
CJB    PUT THE POINTS BEFORE THE AIRFOIL                                
CJB                                                                     
      DX = X(J1,1)                                                      
      CALL CLUSTR(JDIM,ALEFT,0.D0,DX,J1,CLUS)                           
CX10555 CALL DOINX( 555 )                                               
CY10555 CALL DOINY( 555 )                                               
       DO 2 I = 2,J1                                                    
       JJ = J1 - I + 1                                                  
       X(JJ,1) = -CLUS(I)                                               
2      CONTINUE                                                         
CX20556 CALL DOOUTX( 556 )                                              
CY20556 CALL DOOUTY( 556 )                                              
CJB                                                                     
CJB   PUT THE POINTS AFTER THE AIRFOIL                                  
CJB                                                                     
      DX = 1. - X(J2,1)                                                 
      NI = NJ-J2+1                                                      
      XMM = RIGHT - 1.0                                                 
      CALL CLUSTR(JDIM,XMM,0.D0,DX,NI,CLUS)                             
CX10557 CALL DOINX( 557 )                                               
CY10557 CALL DOINY( 557 )                                               
       DO 3 J = J2+1,NJ                                                 
       I = J - J2 + 1                                                   
       X(J,1) = 1.0 + CLUS(I)                                           
3      CONTINUE                                                         
CX20558 CALL DOOUTX( 558 )                                              
CY20558 CALL DOOUTY( 558 )                                              
CX10559 CALL DOINX( 559 )                                               
CY10559 CALL DOINY( 559 )                                               
C                                                                       
      DO 7 J=1,NJ                                                       
      XJ=X(J,1)                                                         
      DO 7 K=2,NK                                                       
      X(J,K)=XJ                                                         
    7 CONTINUE                                                          
CX20560 CALL DOOUTX( 560 )                                              
CY20560 CALL DOOUTY( 560 )                                              
CJB                                                                     
CJB   PUT THE POINTS ABOVE THE AIRFOIL                                  
CJB                                                                     
      CALL CLUSTR(JDIM,TOP,0.D0,DY,NY,CLUS)                             
CX10561 CALL DOINX( 561 )                                               
CY10561 CALL DOINY( 561 )                                               
       DO 5 K = 1,NY                                                    
       Y(1,K) = CLUS(K)                                                 
5      CONTINUE                                                         
CX20562 CALL DOOUTX( 562 )                                              
CY20562 CALL DOOUTY( 562 )                                              
CX10563 CALL DOINX( 563 )                                               
CY10563 CALL DOINY( 563 )                                               
       DO 6 K = 1,NY                                                    
       DO 6 J = 2,NX                                                    
       Y(J,K) = Y(1,K)                                                  
6      CONTINUE                                                         
CX20564 CALL DOOUTX( 564 )                                              
CY20564 CALL DOOUTY( 564 )                                              
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C   !********** SUBROUTINE CLUSTR                                       
      SUBROUTINE CLUSTR(JDIM,YMAX,YMIN,DYM,LMAX,R)                      
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION R(JDIM)                                                 
CX30023 CALL SBINX (23)                                                 
CY30023 CALL SBINY (23)                                                 
C                                                                       
C                                                                       
C                                                                       
      EPS = EPSIL(YMAX,YMIN,DYM,LMAX,0.1D-2,100,1,2)                    
      K = -1                                                            
      R(1) = YMIN                                                       
CX10565 CALL DOINX( 565 )                                               
CY10565 CALL DOINY( 565 )                                               
      DO 5 L = 2,LMAX                                                   
      K = K+1                                                           
      GN = (1.0 +EPS/SQRT(FLOAT(K+4)))**K                               
      R(L) = R(L-1)+DYM*GN                                              
    5 CONTINUE                                                          
CX20566 CALL DOOUTX( 566 )                                              
CY20566 CALL DOOUTY( 566 )                                              
CX40024 CALL SBOUTX (24)                                                
CY40024 CALL SBOUTY (24)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      DOUBLE PRECISION FUNCTION EPSIL                                   
     *  (FMX,FMIN,DFM,NPT,FPCC,ICC,NCALL,KEY)                           
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
                                                                        
      COMMON/IO/IREAD,IWRIT,IVER                                        
CX30025 CALL SBINX (25)                                                 
CY30025 CALL SBINY (25)                                                 
C                                                                       
C                                                                       
      FMXL = FMX                                                        
      FMINL = FMIN                                                      
      DFML = DFM                                                        
      FPCCL = FPCC                                                      
      ICCL = ICC                                                        
C                                                                       
C  SLOWER STRETCHING IF KEY = 2                                         
C                                                                       
      IF(KEY.EQ.2) GO TO 50                                             
C                                                                       
      FNPTM2 = NPT-2                                                    
      IF(NCALL.EQ.1) EPS = (FMXL/DFML)**(1./FNPTM2)-1.0                 
CX10567 CALL DOINX( 567 )                                               
CY10567 CALL DOINY( 567 )                                               
C                                                                       
      DO 3 NIT = 1,ICCL                                                 
      EP1 = EPS+1.                                                      
      EP1TN = EP1**FNPTM2                                               
      DFMOE = DFML/EPS                                                  
      F = FMXL-FMINL-DFMOE*(EP1TN*EP1-1.0)                              
      IF(ABS(F).LT.FPCCL) GO TO 4                                       
      FPN = DFMOE*FNPTM2*EP1TN                                          
      EPS = EPS+F/FPN                                                   
    3 CONTINUE                                                          
CX20568 CALL DOOUTX( 568 )                                              
CY20568 CALL DOOUTY( 568 )                                              
      GO TO 5                                                           
   50 CONTINUE                                                          
      IF(NCALL.EQ.1) EPS=20.                                            
CX10569 CALL DOINX( 569 )                                               
CY10569 CALL DOINY( 569 )                                               
      DO 11 NIT=1,ICCL                                                  
      SUMG = 0.0                                                        
      SUMGP = 0.0                                                       
      DO 12 M = 2,NPT                                                   
      N = M-2                                                           
      GN = (1.0+EPS/SQRT(FLOAT(N+4)))**N                                
      GPN = N*(1.0+EPS/SQRT(FLOAT(N+4)))**(N-1)*(1.0/SQRT(FLOAT(N+4)))  
      SUMG = SUMG+GN                                                    
   12 SUMGP = SUMGP+GPN                                                 
      F = FMXL-FMINL-DFML*SUMG                                          
      IF(ABS(F).LT.FPCCL) GO TO 4                                       
      FP=-DFML*SUMGP                                                    
   11 EPS = EPS-F/FP                                                    
CX20570 CALL DOOUTX( 570 )                                              
CY20570 CALL DOOUTY( 570 )                                              
C                                                                       
    5 CONTINUE                                                          
      EPSIL = EPS                                                       
      WRITE(IWRIT,100)                                                  
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
    4 EPSIL = EPS                                                       
      WRITE(IWRIT,101) EPSIL,F,NIT                                      
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
C                                                                       
  100 FORMAT(/'EXCEEDED MAX NO. OF ITERATIONS IN EPSIL.')               
  101 FORMAT(/'EPSIL=',F12.5,5X,' AND F=',F12.5,5X,'AFTER ',I3,         
     * ' ITERATIONS.')                                                  
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
C                                                                       
      END                                                               
C*********************************************************************  
C********************** MAIN INTEGRATION ROUTINE *********************  
C*********************************************************************  
C                                                                       
      SUBROUTINE INTEGR(JDIM,KDIM,Q,S,PRESS,SNDSP,SLOPE,X,Y,XY,XYJ,     
     >                    XIT,ETT,DS,                                   
     >                    UU,VV,CCX,CCY,COEF2X,COEF2Y,COEF4X,COEF4Y)    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
C                                                                       
C  MAIN INTEGRATION CALLS                                               
C                                                                       
      PARAMETER (MAXJ=289, MAXK=89, 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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       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 S(JDIM*KDIM*4),XY(JDIM*KDIM*4),XYJ(JDIM*KDIM)           
      DIMENSION XIT(JDIM*KDIM),ETT(JDIM*KDIM),DS(JDIM*KDIM)             
      DIMENSION X(JDIM*KDIM),Y(JDIM*KDIM),SLOPE(JDIM)                   
C                                                                       
      DIMENSION UU(JDIM*KDIM),CCX(JDIM*KDIM)                            
      DIMENSION VV(JDIM*KDIM),CCY(JDIM*KDIM)                            
      DIMENSION COEF2X(JDIM*KDIM),COEF2Y(JDIM*KDIM)                     
      DIMENSION COEF4X(JDIM*KDIM),COEF4Y(JDIM*KDIM)                     
C                                                                       
      COMMON/WORK/WORK1(MAXJK,10)                                       
CX30027 CALL SBINX (27)                                                 
CY30027 CALL SBINY (27)                                                 
C                                                                       
C                                                                       
C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
C <<<                    EXPLICIT OPERATORS                        >>>  
C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
C                                                                       
C                                                                       
C*********************************************************************  
C***    0. UPDATE PRESSURE AND SOUND SPEED                         ***  
C*********************************************************************  
C                                                                       
        CALL CALCPS(JDIM,KDIM,Q,PRESS,SNDSP)                            
C                                                                       
C*********************************************************************  
C***   1. EXPLICIT BOUNDARY CONDITIONS                             ***  
C*********************************************************************  
C                                                                       
            CALL BC(JDIM,KDIM,Q,PRESS,SNDSP,XY,XIT,ETT,XYJ,X,Y,SLOPE)   
C                                                                       
C*********************************************************************  
C***   2.  EXPLICIT RHS                                            ***  
C*********************************************************************  
C---------------------------------------------------------------------  
C---      EXPLICIT OPERATOR IN XI                                  ---  
C---------------------------------------------------------------------  
       CALL RHSX(JDIM,KDIM,Q,S,PRESS,XY,XYJ,XIT,WORK1)                  
C                                                                       
C---------------------------------------------------------------------  
C        EXPLICIT OPERATORS IN ETA                                      
C---------------------------------------------------------------------  
       CALL RHSY(JDIM,KDIM,Q,S,PRESS,XY,XYJ,ETT,WORK1)                  
C                                                                       
C                                                                       
C---------------------------------------------------------------------  
C---            3-A-1. CALCULATE EIGENVALUES FOR XI AND ETA          ---
C---------------------------------------------------------------------  
C                                                                       
                IX = 1                                                  
                IY = 2                                                  
                CALL EIGVAL(IX,IY,JDIM,KDIM,Q,PRESS,SNDSP,XYJ,          
     *                                              XY,XIT,UU,CCX)      
                IX = 3                                                  
                IY = 4                                                  
                CALL EIGVAL(IX,IY,JDIM,KDIM,Q,PRESS,SNDSP,XYJ,          
     *                                               XY,ETT,VV,CCY)     
C                                                                       
C                                                                       
C---------------------------------------------------------------------  
C---            3-A-2. FORM SPECTRAL RADII FOR DISSIPATION MODEL ---    
C---------------------------------------------------------------------  
C                                                                       
              CALL SPECT(JDIM,KDIM,UU,VV,CCX,CCY,XYJ,XY,COEF4X,COEF4Y)  
C                                                                       
C*********************************************************************  
C***   3-B. DISSIPATION TERMS                                      ***  
C*********************************************************************  
C                                                                       
C---------------------------------------------------------------------  
C---           3-B-1-A. FORM PRESSURE COEFFICIENTS                ---   
C---------------------------------------------------------------------  
                IXY = 1                                                 
                CALL GRADCO(IXY,JDIM,KDIM,PRESS,XYJ,COEF2X,WORK1)       
                IXY = 2                                                 
                CALL GRADCO(IXY,JDIM,KDIM,PRESS,XYJ,COEF2Y,WORK1)       
C                                                                       
C---------------------------------------------------------------------  
C---            3-B-1-B. ARTIFICIAL DISSIPATION COEFFICIENTS       ---  
C---------------------------------------------------------------------  
C                                                                       
                CALL COEF24(JDIM,KDIM,COEF2X,COEF4X,COEF2Y,COEF4Y,      
     >                      WORK1(1,1),WORK1(1,2))                      
C                                                                       
C---------------------------------------------------------------------  
C---      3-B-2. ARTIFICIAL DISSIPATION                            ---  
C---------------------------------------------------------------------  
C                                                                       
          CALL FILERX(JDIM,KDIM,Q,S,XYJ,COEF2X,COEF4X,WORK1)            
          CALL FILERY(JDIM,KDIM,Q,S,XYJ,COEF2Y,COEF4Y,WORK1)            
C                                                                       
C---------------------------------------------------------------------  
C                                                                       
C*********************************************************************  
C       4.  COMPUTE RESIDUAL AND SCALE BY VARIABLE DT                   
C*********************************************************************  
C---------------------------------------------------------------------  
C---      4-A. COMPUTE L2 AND MAX NORMS OF RESIDUAL                ---  
C---------------------------------------------------------------------  
C                                                                       
C    DENSITY RESIDUAL                                                   
C                                                                       
C                                                                       
C                                                                       
          NRES = 1                                                      
          CALL RESIL2(JDIM,KDIM,S,NRES)                                 
C*********************************************************************  
C*********************************************************************  
C                                                                       
C---------------------------------------------------------------------  
C---      4-B. SCALE WITH VARIABLE DT                              ---  
C---------------------------------------------------------------------  
C                                                                       
          CALL SCALDT(JDIM,KDIM,S,DS)                                   
C                                                                       
C                                                                       
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
C<<<                     IMPLICIT OPERATORS                        >>>  
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
C                                                                       
C---------------------------------------------------------------------  
C    5. SCALAR PENTADIAGONAL ALGORITHM                                  
C---------------------------------------------------------------------  
C                                                                       
C---------------------------------------------------------------------  
C---      5-A. INTEGRATION   XI - ETA ORDER                        ---  
C---------------------------------------------------------------------  
C                                                                       
C.....................................................................  
C...                      -1                                       ...  
C...           5-A-1.    T   S    EIGENVECTOR INVERSE MULTIPLY BY RHS   
C...                      XI       CALL WITH 1,2 FOR XI METRICS    ...  
C.....................................................................  
C                                                                       
              IX = 1                                                    
              IY = 2                                                    
              CALL TKINV(JDIM,KDIM,Q,SNDSP,S,XYJ,XY,IX,IY)              
C                                                                       
C.....................................................................  
C...         5-A-2.    XI DIAGONAL INVERSIONS                      ...  
C.....................................................................  
C                                                                       
             CALL STEPFX(JDIM,KDIM,Q,S,PRESS,SNDSP,XY,XYJ,DS,           
     *                                    COEF2X,COEF4X,UU,CCX,WORK1)   
C                                                                       
C.....................................................................  
C...              N INVERSE MULTIPLY BETWEEN DIAGONAL OPERATORS    ...  
C.....................................................................  
C                                                                       
                  IXX = 1                                               
                  IXY = 2                                               
                  IYX = 3                                               
                  IYY = 4                                               
                  CALL NINVER(JDIM,KDIM,S,XY,IXX,IXY,IYX,IYY)           
C.....................................................................  
C...              ETA DIAGONAL INTEGRATION                         ...  
C.....................................................................  
C                                                                       
                  CALL STEPFY(JDIM,KDIM,Q,S,PRESS,SNDSP,                
     *                            XY,XYJ,DS,COEF2Y,COEF4Y,VV,CCY,WORK1) 
C.....................................................................  
C...              T    S  MATRIX MULTIPLY TO CLOSE DIAGONALIZATION ...  
C...               ETA     3, 4  FOR ETA METRICS                   ...  
C.....................................................................  
C                                                                       
                  IX = 3                                                
                  IY = 4                                                
                  CALL TK(JDIM,KDIM,Q,SNDSP,S,XYJ,XY,IX,IY)             
C                                                                       
C.....................................................................  
C---------------------------------------------------------------------  
C---    END 5-A.                                                   ---  
C---------------------------------------------------------------------  
C*********************************************************************  
C***   END 5.                                                      ***  
C*********************************************************************  
C                                                                       
C*********************************************************************  
C                                                                       
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
C<<<                    FINISH INTEGRATION                         >>>  
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
C                                                                       
C*********************************************************************  
C***    7. UPDATE                                                  ***  
C*********************************************************************  
C                                                                       
        CALL UPDATE(JDIM,KDIM,Q,S)                                      
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
C                                                                       
C                                                                       
                     RETURN                                             
                     END                                                
C***********************************************************************
C*** SECOND -  FOURTH :  ORDER EXPLICIT FILTER *************************
C***********************************************************************
      SUBROUTINE FILERX(JDIM,KDIM,Q,S,XYJ,COEF2,COEF4,WORK)             
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
                                                                        
C                                                                       
      DIMENSION Q(JDIM,KDIM,4),S(JDIM,KDIM,4),XYJ(JDIM,KDIM)            
      DIMENSION COEF2(JDIM,KDIM),COEF4(JDIM,KDIM)                       
C                                                                       
      DIMENSION WORK(JDIM,KDIM,3)                                       
CX30029 CALL SBINX (29)                                                 
CY30029 CALL SBINY (29)                                                 
CX10571 CALL DOINX( 571 )                                               
CY10571 CALL DOINY( 571 )                                               
C                                                                       
C   FOURTH ORDER SMOOTHING, ADDED EXPLICITLY TO RHS                     
C   SECOND ORDER NEAR SHOCKS WITH PRESSURE GRD COEFF.                   
C                                                                       
C  XI DIRECTION                                                         
C                                                                       
C   START DIFFERENCES EACH VARIABLE SEPARATELY                          
C                                                                       
      DO 15 N = 1,4                                                     
C                                                                       
      DO 16 J = JLOW,JUP                                                
      JPL = JPLUS(J)                                                    
      DO 16 K = KLOW,KUP                                                
      WORK(J,K,1) = Q(JPL,K,N)*XYJ(JPL,K) - Q(J,K,N)*XYJ(J,K)           
16    CONTINUE                                                          
C                                                                       
C   C MESH BC                                                           
C                                                                       
         IF(.NOT.PERIDC  )THEN                                          
         J1 = JBEGIN                                                    
         J2 = JEND                                                      
         DO 76 K = KLOW,KUP                                             
         WORK(J1,K,1) = Q(J1+1,K,N)*XYJ(J1+1,K) -                       
     *                    Q(J1,K,N)*XYJ(J1,K)                           
         WORK(J2,K,1) = WORK(J2-1,K,1)                                  
76       CONTINUE                                                       
         ENDIF                                                          
C                                                                       
C   FOR FOURTH ORDER                                                    
C                                                                       
      DO 17 J = JLOW,JUP                                                
      JPL = JPLUS(J)                                                    
      JMI = JMINU (J)                                                   
      DO 17 K = KLOW,KUP                                                
      WORK(J,K,2) = WORK(JPL,K,1) - 2.* WORK(J,K,1) + WORK(JMI,K,1)     
17    CONTINUE                                                          
C                                                                       
C   C MESH BC                                                           
C                                                                       
         IF(.NOT.PERIDC  )THEN                                          
         J1 = JBEGIN                                                    
         J2 = JEND                                                      
         DO 77 K = KLOW,KUP                                             
         WORK(J1,K,2) = Q(J1+2,K,N)*XYJ(J1+2,K) -                       
     *         2.*Q(J1+1,K,N)*XYJ(J1+1,K) + Q(J1,K,N)*XYJ(J1,K)         
         WORK(J2,K,2) = 0.                                              
77       CONTINUE                                                       
         ENDIF                                                          
C                                                                       
C                                                                       
C   FORM DISSIPATION TERM                                               
C                                                                       
      DO 19 J = JLOW,JUP                                                
      DO 19 K = KLOW,KUP                                                
      WORK(J,K,3) = (COEF2(J,K)*WORK(J,K,1) - COEF4(J,K)*WORK(J,K,2))   
19    CONTINUE                                                          
C                                                                       
C   C MESH BC                                                           
C                                                                       
         IF(.NOT.PERIDC  )THEN                                          
         J1 = JBEGIN                                                    
         J2 = JEND                                                      
         DO 79 K = KLOW,KUP                                             
         WORK(J1,K,3) = (COEF2(J1,K)*WORK(J1,K,1) -                     
     >                   COEF4(J1,K)*WORK(J1,K,2))                      
         WORK(J2,K,3) = 0.                                              
79       CONTINUE                                                       
         ENDIF                                                          
C                                                                       
C  ADD IN DISSIPATION                                                   
C                                                                       
C                                                                       
      DTD = DT / (1. + PHIDT)                                           
C                                                                       
      DO 20 J = JLOW,JUP                                                
      JMI = JMINU (J)                                                   
      DO 20 K = KLOW,KUP                                                
      S(J,K,N) = S(J,K,N) + (WORK(J,K,3) - WORK(JMI,K,3))*DTD           
20    CONTINUE                                                          
C                                                                       
15    CONTINUE                                                          
CX20572 CALL DOOUTX( 572 )                                              
CY20572 CALL DOOUTY( 572 )                                              
CX40030 CALL SBOUTX (30)                                                
CY40030 CALL SBOUTY (30)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C*** SECOND -  FOURTH :  ORDER EXPLICIT FILTER *************************
C***********************************************************************
      SUBROUTINE FILERY(JDIM,KDIM,Q,S,XYJ,COEF2,COEF4,WORK)             
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      COMMON/BASE/                                                      
     1  JMAX,       KMAX,     JM,          KM,          JBEGIN,   JEND, 
     1  KBEGIN,     KEND,     JPLUS(999),  JMINU (999), JLOW,     JUP,  
     1  KLOW,       KUP,      PERIDC  ,    NP,          DT,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION Q(JDIM,KDIM,4),S(JDIM,KDIM,4),XYJ(JDIM,KDIM)            
      DIMENSION COEF2(JDIM,KDIM),COEF4(JDIM,KDIM)                       
C                                                                       
      DIMENSION WORK(JDIM,KDIM,3)                                       
CX30031 CALL SBINX (31)                                                 
CY30031 CALL SBINY (31)                                                 
CX10573 CALL DOINX( 573 )                                               
CY10573 CALL DOINY( 573 )                                               
C                                                                       
C   FOURTH ORDER SMOOTHING, ADDED EXPLICITLY TO RHS                     
C   SECOND ORDER NEAR SHOCKS WITH PRESSURE GRD COEFF.                   
C                                                                       
C  ETA DIRECTION                                                        
C                                                                       
C   START DIFFERENCES EACH VARIABLE SEPARATLY                           
C                                                                       
      DO 39 N = 1,4                                                     
C                                                                       
      DO 35 K = KBEGIN,KUP                                              
      KPL = K+1                                                         
      DO 35 J = JLOW,JUP                                                
      WORK(J,K,1) = Q(J,KPL,N)*XYJ(J,KPL) - Q(J,K,N)*XYJ(J,K)           
35    CONTINUE                                                          
C                                                                       
C   FOR FOURTH ORDER                                                    
C                                                                       
      DO 36 K = KLOW,KUP-1                                              
      KPL = K+1                                                         
      KMI = K-1                                                         
      DO 36 J = JLOW,JUP                                                
      WORK(J,K,2) = WORK(J,KPL,1) - 2.* WORK(J,K,1) + WORK(J,KMI,1)     
36    CONTINUE                                                          
C                                                                       
C  BOUNDARY                                                             
C                                                                       
      DO 37 J = JLOW,JUP                                                
      WORK(J,KBEGIN,2) = 0.                                             
      WORK(J,KUP,2) = WORK(J,KUP-1,1) - WORK(J,KUP,1)                   
37    CONTINUE                                                          
C                                                                       
C                                                                       
C   FORM DISSIPATION TERM                                               
C                                                                       
      DO 38 K = KBEGIN,KUP                                              
      DO 38 J = JLOW,JUP                                                
       WORK(J,K,3) = (COEF2(J,K)*WORK(J,K,1) - COEF4(J,K)*WORK(J,K,2))  
38    CONTINUE                                                          
C                                                                       
C  ADD IN DISSIPATION                                                   
C                                                                       
C                                                                       
      DTD = DT / (1. + PHIDT)                                           
C                                                                       
      DO 40 K = KLOW,KUP                                                
      DO 40 J = JLOW,JUP                                                
      S(J,K,N) = S(J,K,N) + (WORK(J,K,3) - WORK(J,K-1,3))*DTD           
40    CONTINUE                                                          
C                                                                       
39    CONTINUE                                                          
CX20574 CALL DOOUTX( 574 )                                              
CY20574 CALL DOOUTY( 574 )                                              
CX40032 CALL SBOUTX (32)                                                
CY40032 CALL SBOUTY (32)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C***************** MAIN INTEGRATION ROUTINE ****************************
C********************* DIAGONAL VERSION ********************************
C********  WITH IMPLICIT FOURTH ORDER DISSIPATION   ********************
C    ROUTINE VECTORIZED  THP  10/18/83                                  
C***********************************************************************
      SUBROUTINE STEPFX(JDIM,KDIM,Q,S,PRESS,SNDSP,XY,XYJ,DS,            
     *                    COEF2,COEF4,UU,CC,WORK)                       
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION Q(JDIM,KDIM,4),PRESS(JDIM,KDIM),SNDSP(JDIM,KDIM)        
      DIMENSION S(JDIM,KDIM,4),XY(JDIM,KDIM,4),XYJ(JDIM,KDIM)           
      DIMENSION DS(JDIM,KDIM)                                           
C                                                                       
      DIMENSION COEF2(JDIM,KDIM),COEF4(JDIM,KDIM)                       
      DIMENSION UU(JDIM,KDIM),CC(JDIM,KDIM)                             
C                                                                       
      DIMENSION WORK(JDIM,KDIM,10)                                      
CX30033 CALL SBINX (33)                                                 
CY30033 CALL SBINY (33)                                                 
C                                                                       
C  NOTE : D USED FOR VECTORIZATION                                      
C                                                                       
      DTD   = DT * THETAD /(1. + PHIDT)                                 
      HD    = 0.5*DTD                                                   
CX10575 CALL DOINX( 575 )                                               
CY10575 CALL DOINY( 575 )                                               
C                                                                       
C  XI DIRECTION                                                         
C                                                                       
C  FILL MATRIX FOR INVERSION FOR DIAGONAL IN XI                         
C                                                                       
C  SET SWITCHES FOR EIGENVALUES                                         
C                                                                       
C  DO FIRST TWO AT THE SAME TIME SINCE THEY HAVE THE SAME COEFFICIENTS  
C                                                                       
         DO 300 N = 2,4                                                 
            IF ( N .EQ. 2) SN = 0.                                      
            IF ( N .EQ. 3) SN = 1.                                      
            IF ( N .EQ. 4) SN = -1.                                     
C                                                                       
         DO 210 J = JLOW,JUP                                            
         JP1 = JPLUS(J)                                                 
         JM1 = JMINU (J)                                                
         JP2 = JPLUS(JP1)                                               
         JM2 = JMINU (JM1)                                              
            DO 212 K = KLOW,KUP                                         
            C2M = COEF2(JM1,K)*DTD                                      
            C4M = COEF4(JM1,K)*DTD                                      
              C2 = COEF2(J,K)*DTD                                       
              C4 = COEF4(J,K)*DTD                                       
            WORK(J,K,1) =  XYJ(JM2,K)*C4M                               
            WORK(J,K,2) = -(C2M + 3.*C4M + C4)*XYJ(JM1,K)               
            WORK(J,K,3) = XYJ(J,K)*(C2M+3.*C4M+C2+3.*C4)                
            WORK(J,K,4) = -(C2 + 3.*C4 + C4M)*XYJ(JP1,K)                
            WORK(J,K,5) = XYJ(JP2,K)*C4                                 
212         CONTINUE                                                    
210       CONTINUE                                                      
C                                                                       
            IF(.NOT.PERIDC  )THEN                                       
         J = JLOW                                                       
         JP1 = JPLUS(J)                                                 
         JM1 = JMINU (J)                                                
         JP2 = JPLUS(JP1)                                               
            DO 220 K = KLOW,KUP                                         
            C2M = COEF2(JM1,K)*DTD                                      
            C4M = COEF4(JM1,K)*DTD                                      
              C2 = COEF2(J,K)*DTD                                       
              C4 = COEF4(J,K)*DTD                                       
            WORK(J,K,1) =  0.                                           
            WORK(J,K,2) = -(C2M + C4M + C4)*XYJ(JM1,K)                  
            WORK(J,K,3) = XYJ(J,K)*(C2M+2.*C4M+C2+3.*C4)                
            WORK(J,K,4) = -(C2 + 3.*C4 + C4M)*XYJ(JP1,K)                
            WORK(J,K,5) = XYJ(JP2,K)*C4                                 
220         CONTINUE                                                    
C                                                                       
         J = JUP                                                        
         JP1 = JPLUS(J)                                                 
         JM1 = JMINU (J)                                                
         JM2 = JMINU (JM1)                                              
            DO 222 K = KLOW,KUP                                         
            C2M = COEF2(JM1,K)*DTD                                      
            C4M = COEF4(JM1,K)*DTD                                      
              C2 = COEF2(J,K)*DTD                                       
              C4 = COEF4(J,K)*DTD                                       
            WORK(J,K,1) =  XYJ(JM2,K)*C4M                               
            WORK(J,K,2) = -(C2M + 3.*C4M + C4)*XYJ(JM1,K)               
            WORK(J,K,3) = XYJ(J,K)*(C2M+3.*C4M+C2+2.*C4)                
            WORK(J,K,4) = -(C2 + C4 + C4M)*XYJ(JP1,K)                   
            WORK(J,K,5) = 0.                                            
222         CONTINUE                                                    
C                                                                       
            ENDIF                                                       
C                                                                       
C   ADD IN FLUX JACOBIANS AND VARIABLE DT AND IDENTITY                  
C                                                                       
         DO 230 J = JLOW,JUP                                            
         JP1 = JPLUS(J)                                                 
         JM1 = JMINU (J)                                                
            DO 232 K = KLOW,KUP                                         
            WORK(J,K,2) = WORK(J,K,2) - (UU(JM1,K)+SN*CC(JM1,K))*HD     
            WORK(J,K,4) = WORK(J,K,4) + (UU(JP1,K)+SN*CC(JP1,K))*HD     
            WORK(J,K,1) = WORK(J,K,1)*DS(J,K)                           
            WORK(J,K,2) = WORK(J,K,2)*DS(J,K)                           
            WORK(J,K,3) = WORK(J,K,3)*DS(J,K)                           
            WORK(J,K,4) = WORK(J,K,4)*DS(J,K)                           
            WORK(J,K,5) = WORK(J,K,5)*DS(J,K)                           
            WORK(J,K,3) = WORK(J,K,3) + 1.                              
232         CONTINUE                                                    
230       CONTINUE                                                      
C                                                                       
C  FOR N = 1 AND 2 DO BOTH TOGETHER                                     
C   INVERT ON FIRST TWO S ELEMENTS                                      
C                                                                       
          IF(N.EQ.2)THEN                                                
C                                                                       
        IF(.NOT.PERIDC  )THEN                                           
        CALL XPENT2(JDIM,KDIM,WORK(1,1,1),WORK(1,1,2),WORK(1,1,3),      
     *       WORK(1,1,4),WORK(1,1,5),WORK(1,1,7),WORK(1,1,8),           
     *       S(1,1,1),JLOW,JUP,KLOW,KUP)                                
C                                                                       
               ELSE                                                     
C              PERIDC                                                   
        CALL XPENP2(JDIM,KDIM,WORK(1,1,1),WORK(1,1,2),WORK(1,1,3),      
     *       WORK(1,1,4),WORK(1,1,5),WORK(1,1,7),WORK(1,1,8),           
     *       WORK(1,1,9),WORK(1,1,10),S(1,1,1),JLOW,JUP,KLOW,KUP)       
               ENDIF                                                    
C                                                                       
            ELSE                                                        
C                                                                       
C  N = 3 OR 4                                                           
C                                                                       
               IF(.NOT.PERIDC  )THEN                                    
        CALL XPENTA(JDIM,KDIM,WORK(1,1,1),WORK(1,1,2),WORK(1,1,3),      
     *       WORK(1,1,4),WORK(1,1,5),WORK(1,1,7),WORK(1,1,8),           
     *       S(1,1,N),JLOW,JUP,KLOW,KUP)                                
               ELSE                                                     
C              PERIDC                                                   
        CALL XPENTP(JDIM,KDIM,WORK(1,1,1),WORK(1,1,2),WORK(1,1,3),      
     *       WORK(1,1,4),WORK(1,1,5),WORK(1,1,7),WORK(1,1,8),           
     *       WORK(1,1,9),WORK(1,1,10),S(1,1,N),JLOW,JUP,KLOW,KUP)       
               ENDIF                                                    
C                                                                       
        ENDIF                                                           
C                                                                       
300       CONTINUE                                                      
CX20576 CALL DOOUTX( 576 )                                              
CY20576 CALL DOOUTY( 576 )                                              
CX40034 CALL SBOUTX (34)                                                
CY40034 CALL SBOUTY (34)                                                
C                                                                       
C  END XI DIRECTION                                                     
C                                                                       
        RETURN                                                          
        END                                                             
C***********************************************************************
C***************** MAIN INTEGRATION ROUTINE ****************************
C********************* DIAGONAL VERSION ********************************
C********  WITH IMPLICIT FOURTH ORDER DISSIPATION   ********************
C    ROUTINE VECTORIZED  THP  10/18/83                                  
C***********************************************************************
      SUBROUTINE STEPFY(JDIM,KDIM,Q,S,PRESS,SNDSP,XY,XYJ,DS,            
     *                    COEF2,COEF4,UU,CC,WORK)                       
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION Q(JDIM,KDIM,4),PRESS(JDIM,KDIM),SNDSP(JDIM,KDIM)        
      DIMENSION S(JDIM,KDIM,4),XY(JDIM,KDIM,4),XYJ(JDIM,KDIM)           
      DIMENSION DS(JDIM,KDIM)                                           
C                                                                       
      DIMENSION COEF2(JDIM,KDIM),COEF4(JDIM,KDIM)                       
      DIMENSION UU(JDIM,KDIM),CC(JDIM,KDIM)                             
C                                                                       
      DIMENSION WORK(JDIM,KDIM,10)                                      
CX30035 CALL SBINX (35)                                                 
CY30035 CALL SBINY (35)                                                 
C                                                                       
C                                                                       
      DTD   = DT * THETAD /(1. + PHIDT)                                 
      HD    = 0.5*DTD                                                   
CX10577 CALL DOINX( 577 )                                               
CY10577 CALL DOINY( 577 )                                               
C                                                                       
C  ETA INVERSION PART OF DIAGONAL SCHEME                                
C                                                                       
C  SET EIGENVALUE SWTICHES                                              
C                                                                       
C  DO N = 1 AND 2 TOGETHER                                              
C                                                                       
            DO 420 N = 2,4                                              
               IF ( N .EQ. 2) SN = 0.                                   
               IF ( N .EQ. 3) SN = 1.                                   
               IF ( N .EQ. 4) SN = -1.                                  
C                                                                       
C                                                                       
C  FILL TRIDIAGONALS                                                    
C                                                                       
          K   = KLOW                                                    
          KM1 = K-1                                                     
          DO 422 J = JLOW,JUP                                           
            C2M = COEF2(J,KM1)*DTD                                      
            C4M = COEF4(J,KM1)*DTD                                      
              C2 = COEF2(J,K)*DTD                                       
              C4 = COEF4(J,K)*DTD                                       
            WORK(J,K,1) =  0.                                           
            WORK(J,K,2) = -(C2M + C4)*XYJ(J,K-1)                        
            WORK(J,K,3) = XYJ(J,K)*(C2M+C2+3.*C4)                       
            WORK(J,K,4) = -(C2 + 3.*C4 )*XYJ(J,K+1)                     
            WORK(J,K,5) = XYJ(J,K+2)*C4                                 
422       CONTINUE                                                      
C                                                                       
          K   = KUP                                                     
          KM1 = K-1                                                     
          DO 426 J = JLOW,JUP                                           
            C2M = COEF2(J,KM1)*DTD                                      
            C4M = COEF4(J,KM1)*DTD                                      
              C2 = COEF2(J,K)*DTD                                       
              C4 = COEF4(J,K)*DTD                                       
            WORK(J,K,1) =  XYJ(J,K-2)*C4M                               
            WORK(J,K,2)  = -(C2M + 3.*C4M + C4)*XYJ(J,K-1)              
            WORK(J,K,3)  = XYJ(J,K)*(C2M+3.*C4M+C2+2.*C4)               
            WORK(J,K,4)  = -(C2 + C4 + C4M)*XYJ(J,K+1)                  
            WORK(J,K,5)  = 0.                                           
426       CONTINUE                                                      
C                                                                       
C   FOURTH ORDER IN INTERIOR                                            
C                                                                       
        DO 428 K = KLOW+1,KUP-1                                         
          KM1 = K-1                                                     
          DO 430 J = JLOW,JUP                                           
            C2M = COEF2(J,KM1)*DTD                                      
            C4M = COEF4(J,KM1)*DTD                                      
              C2 = COEF2(J,K)*DTD                                       
              C4 = COEF4(J,K)*DTD                                       
            WORK(J,K,1)  =  XYJ(J,K-2)*C4M                              
            WORK(J,K,2)  = -(C2M + 3.*C4M + C4)*XYJ(J,K-1)              
            WORK(J,K,3)  = XYJ(J,K)*(C2M+3.*C4M+C2+3.*C4)               
            WORK(J,K,4)  = -(C2 + 3.*C4 + C4M)*XYJ(J,K+1)               
            WORK(J,K,5)  = XYJ(J,K+2)*C4                                
430       CONTINUE                                                      
428     CONTINUE                                                        
C                                                                       
C                                                                       
C  ADD IN FLUX JACOBIANS, VARIABLE DT AND IDENTITY                      
C                                                                       
        DO 435 K = KLOW,KUP                                             
          DO 435 J = JLOW,JUP                                           
          WORK(J,K,2) = WORK(J,K,2) - (UU(J,K-1)+SN*CC(J,K-1))*HD       
          WORK(J,K,4) = WORK(J,K,4) + (UU(J,K+1)+SN*CC(J,K+1))*HD       
          WORK(J,K,1) = WORK(J,K,1)*DS(J,K)                             
          WORK(J,K,2) = WORK(J,K,2)*DS(J,K)                             
          WORK(J,K,3) = WORK(J,K,3)*DS(J,K)                             
          WORK(J,K,4) = WORK(J,K,4)*DS(J,K)                             
          WORK(J,K,5) = WORK(J,K,5)*DS(J,K)                             
          WORK(J,K,3) = WORK(J,K,3) + 1.0                               
435     CONTINUE                                                        
C                                                                       
C   COMBINE N = 1 AND 2                                                 
C    INVERT ON FIRST TWO S ELEMENTS                                     
C                                                                       
         IF(N.EQ.2)THEN                                                 
C                                                                       
        CALL YPENT2(JDIM,KDIM,WORK(1,1,1),WORK(1,1,2),WORK(1,1,3),      
     *       WORK(1,1,4),WORK(1,1,5),WORK(1,1,7),WORK(1,1,8),           
     *       S(1,1,1),KLOW,KUP,JLOW,JUP)                                
C                                                                       
        ELSE                                                            
C  DO N = 3 AND 4                                                       
C                                                                       
        CALL YPENTA(JDIM,KDIM,WORK(1,1,1),WORK(1,1,2),WORK(1,1,3),      
     *       WORK(1,1,4),WORK(1,1,5),WORK(1,1,7),WORK(1,1,8),           
     *       S(1,1,N),KLOW,KUP,JLOW,JUP)                                
C                                                                       
         ENDIF                                                          
C                                                                       
420         CONTINUE                                                    
CX20578 CALL DOOUTX( 578 )                                              
CY20578 CALL DOOUTY( 578 )                                              
CX40036 CALL SBOUTX (36)                                                
CY40036 CALL SBOUTY (36)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
C***********************************************************************
C**************  COMPUTE EIGENVALUES ***********************************
C***********************************************************************
      SUBROUTINE SPECT(JDIM,KDIM,UU,VV,CCX,CCY,XYJ,XY,SPECTX,SPECTY)    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
C                                                                       
C                                                                       
C  FORM SPECTX AND SPECTY INDEPENDENTLY                                 
C                                                                       
C     SPECTX = (ABS(UU) + CCX)/XYJ                                      
C     SPECTY = (ABS(VV) + CCY)/XYJ                                      
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION XYJ(JDIM,KDIM),XY(JDIM,KDIM,4)                          
      DIMENSION SPECTX(JDIM,KDIM),UU(JDIM,KDIM),CCX(JDIM,KDIM)          
      DIMENSION SPECTY(JDIM,KDIM),VV(JDIM,KDIM),CCY(JDIM,KDIM)          
CX30037 CALL SBINX (37)                                                 
CY30037 CALL SBINY (37)                                                 
CX10579 CALL DOINX( 579 )                                               
CY10579 CALL DOINY( 579 )                                               
C                                                                       
C   COMPUTE DISSIPATION SCALING                                         
C                                                                       
C                                                                       
C  SPECTRAL RADIUS                                                      
C                                                                       
           DO 100 K=KBEGIN,KEND                                         
           DO 100 J=JBEGIN,JEND                                         
C                                                                       
           SIGX = ABS(UU(J,K)) + CCX(J,K)                               
           SIGY = ABS(VV(J,K)) + CCY(J,K)                               
C                                                                       
           RJ = 1./XYJ(J,K)                                             
           SPECTX(J,K) = SIGX*RJ                                        
           SPECTY(J,K) = SIGY*RJ                                        
C                                                                       
100        CONTINUE                                                     
CX20580 CALL DOOUTX( 580 )                                              
CY20580 CALL DOOUTY( 580 )                                              
CX40038 CALL SBOUTX (38)                                                
CY40038 CALL SBOUTY (38)                                                
C                                                                       
                RETURN                                                  
                END                                                     
C***********************************************************************
C***********************************************************************
C*** GRADIENT COEFFICENT FOR ARTIFICIAL DISSIPATION ********************
C***********************************************************************
      SUBROUTINE GRADCO(IXY,JDIM,KDIM,PRESS,XYJ,COEF,PRSS)              
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
C                                                                       
C  IXY = 1  FOR XI  DIRECTION                                           
C  IXY = 2  FOR ETA DIRECTION                                           
      COMMON/BASE/                                                      
     1  JMAX,       KMAX,     JM,          KM,          JBEGIN,   JEND, 
     1  KBEGIN,     KEND,     JPLUS(999),  JMINU (999), JLOW,     JUP,  
     1  KLOW,       KUP,      PERIDC  ,    NP,          DT,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION PRESS(JDIM,KDIM),XYJ(JDIM,KDIM)                         
      DIMENSION COEF(JDIM,KDIM),PRSS(JDIM,KDIM)                         
CX30039 CALL SBINX (39)                                                 
CY30039 CALL SBINY (39)                                                 
C                                                                       
       IF(IXY.EQ.1)THEN                                                 
C                                                                       
            JEDGE1 = JLOW                                               
            JEDGE2 = JUP                                                
CX10581 CALL DOINX( 581 )                                               
CY10581 CALL DOINY( 581 )                                               
C                                                                       
C   COMPUTE PRESSURE GRAD                                               
C                                                                       
      DO 1 K = KBEGIN,KEND                                              
      DO 1 J = JBEGIN,JEND                                              
      COEF(J,K) = 0.                                                    
      PRSS(J,K) = PRESS(J,K)*XYJ(J,K)                                   
1     CONTINUE                                                          
CX20582 CALL DOOUTX( 582 )                                              
CY20582 CALL DOOUTY( 582 )                                              
C                                                                       
C   BYPASS LOGIC IF DIS2X = 0.0                                         
C                                                                       
      IF(DIS2X.EQ.0.0) THEN                                             
CX40040 CALL SBOUTX (40)                                                
CY40040 CALL SBOUTY (40)                                                
      RETURN                                                            
      ENDIF                                                             
CX10583 CALL DOINX( 583 )                                               
CY10583 CALL DOINY( 583 )                                               
C                                                                       
      DO 2 J = JLOW,JUP                                                 
      JPL = JPLUS(J)                                                    
      JMI = JMINU (J)                                                   
      DO 2 K = KLOW,KUP                                                 
      R1 = PRSS(JPL,K) -2.* PRSS(J,K) + PRSS(JMI,K)                     
      R2 = PRSS(JPL,K) +2.* PRSS(J,K) + PRSS(JMI,K)                     
      COEF(J,K) = ABS( R1/R2 )                                          
2     CONTINUE                                                          
CX20584 CALL DOOUTX( 584 )                                              
CY20584 CALL DOOUTY( 584 )                                              
C                                                                       
C   C MESH BC                                                           
C                                                                       
         IF(.NOT.PERIDC  )THEN                                          
CX10585 CALL DOINX( 585 )                                               
CY10585 CALL DOINY( 585 )                                               
         DO 3 K = KLOW,KUP                                              
         COEF(JBEGIN,K) = COEF(JLOW,K)                                  
         COEF(JEND,K) = COEF(JUP,K)                                     
         PRSS(JBEGIN,K) = COEF(JLOW,K)                                  
         PRSS(JEND,K) = COEF(JUP,K)                                     
3        CONTINUE                                                       
CX20586 CALL DOOUTX( 586 )                                              
CY20586 CALL DOOUTY( 586 )                                              
         ENDIF                                                          
CX10587 CALL DOINX( 587 )                                               
CY10587 CALL DOINY( 587 )                                               
C                                                                       
C  SMOOTH OUT PRESSURE COEFFICIENT                                      
C                                                                       
      DO 4 J = JEDGE1,JEDGE2                                            
      JPL = JPLUS(J)                                                    
      JPLPL = JPLUS(JPL)                                                
      JMI = JMINU (J)                                                   
      DO 4 K = KLOW,KUP                                                 
      PRSS(J,K) = MAX( COEF(JMI,K),COEF(J,K),COEF(JPL,K)                
     *                    ,COEF(JPLPL,K) )                              
4     CONTINUE                                                          
CX20588 CALL DOOUTX( 588 )                                              
CY20588 CALL DOOUTY( 588 )                                              
CX10589 CALL DOINX( 589 )                                               
CY10589 CALL DOINY( 589 )                                               
C                                                                       
      DO 5 J = JLOW,JUP                                                 
      JPL = JPLUS(J)                                                    
      JMI = JMINU (J)                                                   
      DO 5 K = KLOW,KUP                                                 
      COEF(J,K) = 0.5*PRSS(J,K) + 0.25*(PRSS(JPL,K)+PRSS(JMI,K))        
5     CONTINUE                                                          
CX20590 CALL DOOUTX( 590 )                                              
CY20590 CALL DOOUTY( 590 )                                              
C                                                                       
      IF(.NOT.PERIDC  )THEN                                             
C  TURN OFF SECOND ORDER NEAR TRAILING EDGE                             
CX10591 CALL DOINX( 591 )                                               
CY10591 CALL DOINY( 591 )                                               
C                                                                       
C   C MESH BC                                                           
C                                                                       
         DO 12 K = KLOW,KUP                                             
         COEF(JBEGIN,K) = COEF(JLOW,K)                                  
         COEF(JEND,K) = COEF(JUP,K)                                     
12       CONTINUE                                                       
CX20592 CALL DOOUTX( 592 )                                              
CY20592 CALL DOOUTY( 592 )                                              
C                                                                       
         ENDIF                                                          
C                                                                       
       ELSEIF(IXY .EQ. 2)THEN                                           
CX10593 CALL DOINX( 593 )                                               
CY10593 CALL DOINY( 593 )                                               
C                                                                       
C                                                                       
C  GRADIANT COEFFICENT IN ETA DIRECTION                                 
C                                                                       
      DO 21 K = KBEGIN,KEND                                             
      DO 21 J = JBEGIN,JEND                                             
      COEF(J,K) = 0.                                                    
      PRSS(J,K) = PRESS(J,K)*XYJ(J,K)                                   
21    CONTINUE                                                          
CX20594 CALL DOOUTX( 594 )                                              
CY20594 CALL DOOUTY( 594 )                                              
C                                                                       
C   BYPASS LOGIC IF DIS2Y = 0.0                                         
C                                                                       
      IF(DIS2Y.EQ.0.0) THEN                                             
CX40040 CALL SBOUTX (40)                                                
CY40040 CALL SBOUTY (40)                                                
      RETURN                                                            
      ENDIF                                                             
CX10595 CALL DOINX( 595 )                                               
CY10595 CALL DOINY( 595 )                                               
C                                                                       
      DO 22 K = KLOW,KUP                                                
      KPL = K+1                                                         
      KMI = K-1                                                         
      DO 22 J = JLOW,JUP                                                
      R1 = PRSS(J,KPL) -2.* PRSS(J,K) + PRSS(J,KMI)                     
      R2 = PRSS(J,KPL) +2.* PRSS(J,K) + PRSS(J,KMI)                     
      COEF(J,K) = ABS( R1/R2 )                                          
22    CONTINUE                                                          
CX20596 CALL DOOUTX( 596 )                                              
CY20596 CALL DOOUTY( 596 )                                              
CX10597 CALL DOINX( 597 )                                               
CY10597 CALL DOINY( 597 )                                               
C                                                                       
C   BC                                                                  
C                                                                       
         DO 24 J = JLOW,JUP                                             
         COEF(J,KBEGIN) = COEF(J,KLOW)                                  
         COEF(J,KEND) = COEF(J,KUP)                                     
         PRSS(J,KBEGIN) = COEF(J,KLOW)                                  
         PRSS(J,KEND) = COEF(J,KUP)                                     
24       CONTINUE                                                       
CX20598 CALL DOOUTX( 598 )                                              
CY20598 CALL DOOUTY( 598 )                                              
CX10599 CALL DOINX( 599 )                                               
CY10599 CALL DOINY( 599 )                                               
C                                                                       
C  SMOOTH OUT PRESSURE COEFFICIENT                                      
C                                                                       
      DO 26 K = KLOW,KUP                                                
      KPL = K+1                                                         
      KMI = K-1                                                         
      DO 26 J = JLOW,JUP                                                
      PRSS(J,K) = MAX( COEF(J,KMI),COEF(J,K),COEF(J,KPL))               
26    CONTINUE                                                          
CX20600 CALL DOOUTX( 600 )                                              
CY20600 CALL DOOUTY( 600 )                                              
CX10601 CALL DOINX( 601 )                                               
CY10601 CALL DOINY( 601 )                                               
C                                                                       
      DO 28 K = KLOW,KUP                                                
      KPL = K+1                                                         
      KMI = K-1                                                         
      DO 28 J = JLOW,JUP                                                
      COEF(J,K) = 0.5*PRSS(J,K) + 0.25*(PRSS(J,KPL)+PRSS(J,KMI))        
28    CONTINUE                                                          
CX20602 CALL DOOUTX( 602 )                                              
CY20602 CALL DOOUTY( 602 )                                              
C                                                                       
       ENDIF                                                            
CX40040 CALL SBOUTX (40)                                                
CY40040 CALL SBOUTY (40)                                                
C                                                                       
         RETURN                                                         
         END                                                            
C***********************************************************************
C**************  COMPUTE NONLINEAR 2ND AND 4TH ORDER COEFFICIENTS  *****
C***********************************************************************
      SUBROUTINE COEF24(JDIM,KDIM,COEF2X,COEF4X,COEF2Y,COEF4Y,          
     >           PCOEF,SPECT)                                           
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION COEF2X(JDIM,KDIM),COEF4X(JDIM,KDIM)                     
      DIMENSION COEF2Y(JDIM,KDIM),COEF4Y(JDIM,KDIM)                     
      DIMENSION PCOEF(JDIM,KDIM),SPECT(JDIM,KDIM)                       
CX30041 CALL SBINX (41)                                                 
CY30041 CALL SBINY (41)                                                 
CX10603 CALL DOINX( 603 )                                               
CY10603 CALL DOINY( 603 )                                               
C                                                                       
C   COEF2X,Y COMES IN AS COEFX,Y (PRESSURE GRADIANT COEFFICIENT)        
C   COEF4X,Y COMES IN AS SPECX,Y (SPECTRAL RADIUS)                      
C                                                                       
      DO 1 K = KBEGIN,KEND                                              
      DO 1 J = JBEGIN,JEND                                              
         PCOEF(J,K) = COEF2Y(J,K)                                       
         SPECT(J,K) = COEF4Y(J,K)                                       
1     CONTINUE                                                          
CX20604 CALL DOOUTX( 604 )                                              
CY20604 CALL DOOUTY( 604 )                                              
C                                                                       
       EPS4Y = DIS4Y/64.                                                
       EPS2Y = DIS2Y                                                    
CX10605 CALL DOINX( 605 )                                               
CY10605 CALL DOINY( 605 )                                               
C                                                                       
C  FORM 2ND AND 4TH ORDER DISSIPATION COEFFICIENTS IN Y                 
C                                                                       
       DO 10 K = KBEGIN,KUP                                             
       DO 10 J = JLOW,JUP                                               
C                                                                       
         FIL = SPECT(J,K+1) + SPECT(J,K)                                
         C2 = PCOEF(J,K)*FIL*EPS2Y                                      
         C4 = EPS4Y*FIL                                                 
         C4 = C4 - MIN(C4,C2)                                           
         COEF2Y(J,K) = C2                                               
         COEF4Y(J,K) = C4                                               
10       CONTINUE                                                       
CX20606 CALL DOOUTX( 606 )                                              
CY20606 CALL DOOUTY( 606 )                                              
CX10607 CALL DOINX( 607 )                                               
CY10607 CALL DOINY( 607 )                                               
                                                                        
C                                                                       
C                                                                       
      DO 2 K = KBEGIN,KEND                                              
      DO 2 J = JBEGIN,JEND                                              
         PCOEF(J,K) = COEF2X(J,K)                                       
         SPECT(J,K) = COEF4X(J,K)                                       
2     CONTINUE                                                          
CX20608 CALL DOOUTX( 608 )                                              
CY20608 CALL DOOUTY( 608 )                                              
C                                                                       
       EPS4X = DIS4X/64.                                                
       EPS2X = DIS2X                                                    
CX10609 CALL DOINX( 609 )                                               
CY10609 CALL DOINY( 609 )                                               
C                                                                       
C  FORM 2ND AND 4TH ORDER DISSIPATION COEFFICIENTS IN Y                 
C                                                                       
       DO 20 J = JBEGIN,JUP                                             
       JPL = JPLUS(J)                                                   
       DO 20 K = KLOW,KUP                                               
C                                                                       
         FIL = SPECT(JPL,K) + SPECT(J,K)                                
         C2 = PCOEF(J,K)*FIL*EPS2X                                      
         C4 = EPS4X*FIL                                                 
         C4 = C4 - MIN(C4,C2)                                           
         COEF2X(J,K) = C2                                               
         COEF4X(J,K) = C4                                               
20       CONTINUE                                                       
CX20610 CALL DOOUTX( 610 )                                              
CY20610 CALL DOOUTY( 610 )                                              
CX40042 CALL SBOUTX (42)                                                
CY40042 CALL SBOUTY (42)                                                
                                                                        
         RETURN                                                         
         END                                                            
C***********************************************************************
C*************  SCALE RHS WITH VARIABLE DT     *************************
C***********************************************************************
      SUBROUTINE SCALDT(JDIM,KDIM,S,DS)                                 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION S(JDIM,KDIM,4),DS(JDIM,KDIM)                            
CX30043 CALL SBINX (43)                                                 
CY30043 CALL SBINY (43)                                                 
CX10611 CALL DOINX( 611 )                                               
CY10611 CALL DOINY( 611 )                                               
C                                                                       
C                                                                       
C  SCALE WITH SPATIALLY VARIABLE TIME STEP                              
C                                                                       
      DO 10 M = 1,4                                                     
      DO 10 K = KBEGIN,KEND                                             
      DO 10 J = JBEGIN,JEND                                             
         S(J,K,M) = S(J,K,M)*DS(J,K)                                    
10    CONTINUE                                                          
CX20612 CALL DOOUTX( 612 )                                              
CY20612 CALL DOOUTY( 612 )                                              
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C***************** TK MULTIPLY ON S ************************************
C***********************************************************************
      SUBROUTINE TK(JDIM,KDIM,Q,SNDSP,S,XYJ,XY,IX,IY)                   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
C                                                                       
C  IX = 1, IY = 2  XI  METRICS                                          
C  IX = 3, IY = 4  ETA METRICS                                          
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION Q(JDIM,KDIM,4),S(JDIM,KDIM,4),SNDSP(JDIM,KDIM)          
      DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM)                          
CX30045 CALL SBINX (45)                                                 
CY30045 CALL SBINY (45)                                                 
C                                                                       
      BT = 1./SQRT(2.)                                                  
      RGAMI = 1./GAMI                                                   
CX10613 CALL DOINX( 613 )                                               
CY10613 CALL DOINY( 613 )                                               
C                                                                       
C  T  MATRIX MULTIPLY FOR DIAGONAL ALGORITHM                            
C   K                                                                   
C                                                                       
      DO 1 K = KLOW, KUP                                                
       DO 1 J = JLOW,JUP                                                
          RX      = XY(J,K,IX)                                          
          RY      = XY(J,K,IY)                                          
          RRXY    = 1./SQRT(RX**2+RY**2)                                
          RX      = RX*RRXY                                             
          RY      = RY*RRXY                                             
             RHO    = Q(J,K,1)*XYJ(J,K)                                 
             RHOINV = 1./Q(J,K,1)                                       
             U = Q(J,K,2)*RHOINV                                        
             V = Q(J,K,3)*RHOINV                                        
             ALP0 = RHO*BT/SNDSP(J,K)                                   
                 S1 = S(J,K,1)                                          
                 S2 = S(J,K,2)                                          
                 SA = (S(J,K,3) + S(J,K,4))*ALP0                        
                 SM = (S(J,K,3) - S(J,K,4))*RHO*BT                      
                 UV2 = 0.5*(U**2 + V**2)                                
                 RRX = RHO*RX                                           
                 RRY = RHO*RY                                           
C                                                                       
                    S(J,K,1) = S1 + SA                                  
                    S(J,K,2) = U*S1 + RRY*S2 + U*SA + RX*SM             
                    S(J,K,3) = V*S1 - RRX*S2 + V*SA + RY*SM             
                    S(J,K,4) = UV2*S1 + (RRY*U - RRX*V)*S2 +            
     >                         (UV2+SNDSP(J,K)**2*RGAMI)*SA +           
     >                         (RX*U + RY*V)*SM                         
C                                                                       
1       CONTINUE                                                        
CX20614 CALL DOOUTX( 614 )                                              
CY20614 CALL DOOUTY( 614 )                                              
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
C                                                                       
          RETURN                                                        
          END                                                           
C***********************************************************************
C***************** TK INVERSE MULTIPLY ON S ****************************
C***********************************************************************
      SUBROUTINE TKINV(JDIM,KDIM,Q,SNDSP,S,XYJ,XY,IX,IY)                
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
C                                                                       
C  IX = 1, IY = 2  XI  METRICS                                          
C  IX = 3, IY = 4  ETA METRICS                                          
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION Q(JDIM,KDIM,4),S(JDIM,KDIM,4),SNDSP(JDIM,KDIM)          
      DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM)                          
CX30047 CALL SBINX (47)                                                 
CY30047 CALL SBINY (47)                                                 
C                                                                       
      BT = 1./SQRT(2.)                                                  
CX10615 CALL DOINX( 615 )                                               
CY10615 CALL DOINY( 615 )                                               
C                                                                       
C  T INVERSE MATRIX MULTIPLY FOR DIAGONAL ALGORITHM                     
C   K                                                                   
C                                                                       
      DO 1 K = KLOW, KUP                                                
       DO 1 J = JLOW,JUP                                                
         RX        = XY(J,K,IX)                                         
         RY        = XY(J,K,IY)                                         
         RRXY      = 1./SQRT(RX**2+RY**2)                               
         RX        = RX*RRXY                                            
         RY        = RY*RRXY                                            
            RHO       = Q(J,K,1)*XYJ(J,K)                               
            RHOINV    = 1./RHO                                          
            RHOINR   = RHOINV*XYJ(J,K)                                  
            U = Q(J,K,2)*RHOINR                                         
            V = Q(J,K,3)*RHOINR                                         
            SNR   = 1./SNDSP(J, K)                                      
            BET0  = BT*RHOINV*SNR                                       
C                                                                       
               TERM1 = GAMI*( U*S(J,K,2) + V*S(J,K,3) - S(J,K,4)        
     >                 - 0.5*(U**2+V**2)*S(J,K,1) )                     
C                                                                       
               S1    =  S(J,K,1) + TERM1*SNR**2                         
               S2    = ( (- RY*U + RX*V)*S(J,K,1)                       
     >                 + RY*S(J,K,2) - RX*S(J,K,3) )*RHOINV             
C                                                                       
               TERM1 =  TERM1*BET0                                      
               TERM2 =  ((RX*U + RY*V)*S(J,K,1)                         
     >                 - RX*S(J,K,2) - RY*S(J,K,3) )*BT*RHOINV          
C                                                                       
               S(J,K,1) = S1                                            
               S(J,K,2) = S2                                            
               S(J,K,3) = -(TERM1 + TERM2)                              
               S(J,K,4) = -(TERM1 - TERM2)                              
C                                                                       
1       CONTINUE                                                        
CX20616 CALL DOOUTX( 616 )                                              
CY20616 CALL DOOUTY( 616 )                                              
CX40048 CALL SBOUTX (48)                                                
CY40048 CALL SBOUTY (48)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C***********************************************************************
C************* EVALUATION OF RIGHT HAND SIDE ( EXPLICIT PART ) *********
C***********************************************************************
      SUBROUTINE RHSY(JDIM,KDIM,Q,S,PRESS,XY,XYJ,ETT,FLUX)              
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION Q(JDIM,KDIM,4),PRESS(JDIM,KDIM),S(JDIM,KDIM,4)          
      DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM),ETT(JDIM,KDIM)           
      DIMENSION FLUX(JDIM,KDIM,4)                                       
CX30049 CALL SBINX (49)                                                 
CY30049 CALL SBINY (49)                                                 
CX10617 CALL DOINX( 617 )                                               
CY10617 CALL DOINY( 617 )                                               
C                                                                       
C      CENTRAL DIFFERENCING USED IN ETA DIRECTION                       
      DO 20 K=KBEGIN,KEND                                               
      DO 20 J=JLOW,JUP                                                  
         R1 = XY(J,K,3)                                                 
         R2 = XY(J,K,4)                                                 
         R0 = ETT(J,K)                                                  
C       RR IS J/RHO                                                     
         RR = 1./Q(J,K,1)                                               
         U = Q(J,K,2)*RR                                                
         V = Q(J,K,3)*RR                                                
C                                                                       
C       R1, R2 ARE EITHER D XI /DX, D ETA/DX, OR D XI /DY, D ETA/DY     
C       QS IS EITHER CAP-U OR CAP-V ( CONTRAVARIANT VELOCITIES )        
         QS = R0 +R1*U + R2*V                                           
C                                                                       
C       PRESSJ IS PRESSURE / J                                          
         PRESSJ= PRESS(J,K)                                             
C                                                                       
         FLUX(J,K,1)= Q(J,K,1)*QS                                       
         FLUX(J,K,2)= Q(J,K,2)*QS + R1*PRESSJ                           
         FLUX(J,K,3)= Q(J,K,3)*QS + R2*PRESSJ                           
         FLUX(J,K,4) = QS*( Q(J,K,4) + PRESSJ) - PRESSJ*R0              
C                                                                       
20    CONTINUE                                                          
CX20618 CALL DOOUTX( 618 )                                              
CY20618 CALL DOOUTY( 618 )                                              
C                                                                       
      R0 = - .5*DT / (1. + PHIDT)                                       
CX10619 CALL DOINX( 619 )                                               
CY10619 CALL DOINY( 619 )                                               
C          CENTRAL DIFFERENCE F                                         
      DO 30 N=1,4                                                       
      DO 30 K=KLOW,KUP                                                  
      DO 30 J = JLOW,JUP                                                
           S(J,K,N) = R0*( FLUX(J,K+1,N) - FLUX(J,K-1,N))               
     1                 + S(J,K,N)                                       
30    CONTINUE                                                          
CX20620 CALL DOOUTX( 620 )                                              
CY20620 CALL DOOUTY( 620 )                                              
CX40050 CALL SBOUTX (50)                                                
CY40050 CALL SBOUTY (50)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C************* EVALUATION OF RIGHT HAND SIDE ( EXPLICIT PART ) *********
C***********************************************************************
      SUBROUTINE RHSX(JDIM,KDIM,Q,S,PRESS,XY,XYJ,XIT,FLUX)              
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION Q(JDIM,KDIM,4),PRESS(JDIM,KDIM),S(JDIM,KDIM,4)          
      DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM),XIT(JDIM,KDIM)           
      DIMENSION FLUX(JDIM,KDIM,4)                                       
CX30051 CALL SBINX (51)                                                 
CY30051 CALL SBINY (51)                                                 
CX10621 CALL DOINX( 621 )                                               
CY10621 CALL DOINY( 621 )                                               
C                                                                       
C      CENTRAL DIFFERENCING USED IN XI                                  
C     COMPUTE FLUX VECTORS                                              
C                                                                       
      DO 200 K=KLOW,KUP                                                 
      DO 200 J=JBEGIN,JEND                                              
         R1 = XY(J,K,1)                                                 
         R2 = XY(J,K,2)                                                 
         R0 = XIT(J,K)                                                  
C       RR IS J/RHO                                                     
         RR = 1./Q(J,K,1)                                               
         U = Q(J,K,2)*RR                                                
         V = Q(J,K,3)*RR                                                
C                                                                       
C      R1, R2 ARE EITHER D XI /DX, D ETA/DX, OR D XI /DY, D ETA/DY      
C      QS IS EITHER CAP-U OR CAP-V ( CONTRAVARIANT VELOCITIES )         
C                                                                       
         QS = R0 +R1*U + R2*V                                           
C                                                                       
C      PRESSJ IS PRESSURE / J                                           
C                                                                       
         PRESSJ= PRESS(J,K)                                             
C                                                                       
         FLUX(J,K,1)= Q(J,K,1)*QS                                       
         FLUX(J,K,2)= Q(J,K,2)*QS + R1*PRESSJ                           
         FLUX(J,K,3)= Q(J,K,3)*QS + R2*PRESSJ                           
         FLUX(J,K,4) = QS*( Q(J,K,4) + PRESSJ) - PRESSJ*R0              
C                                                                       
200   CONTINUE                                                          
CX20622 CALL DOOUTX( 622 )                                              
CY20622 CALL DOOUTY( 622 )                                              
C                                                                       
C                                                                       
      RT0 = - .5*DT / (1. + PHIDT)                                      
CX10623 CALL DOINX( 623 )                                               
CY10623 CALL DOINY( 623 )                                               
      DO 400 N=1,4                                                      
      DO 400 J=JLOW,JUP                                                 
         JP1 = JPLUS(J)                                                 
         JM1 = JMINU (J)                                                
         DO 400 K = KLOW,KUP                                            
            S(J,K,N) = S(J,K,N) +                                       
     *          RT0*( FLUX(JP1,K,N) - FLUX(JM1,K,N))                    
400   CONTINUE                                                          
CX20624 CALL DOOUTX( 624 )                                              
CY20624 CALL DOOUTY( 624 )                                              
CX40052 CALL SBOUTX (52)                                                
CY40052 CALL SBOUTY (52)                                                
C                                                                       
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C**************  COMPUTE EIGENVALUES ***********************************
C***********************************************************************
      SUBROUTINE EIGVAL(IX,IY,JDIM,KDIM,Q,PRESS,SNDSP,XYJ,XY,TMET,      
     *                   UU,CC)                                         
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
C                                                                       
C  IX = 1  IY = 2 FOR XI   DIRECTION  METRICS                           
C  IX = 3  IY = 4 FOR ETA  DIRECTION  METRICS                           
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION Q(JDIM,KDIM,4),PRESS(JDIM,KDIM),SNDSP(JDIM,KDIM)        
      DIMENSION XY(JDIM,KDIM,4)                                         
      DIMENSION TMET(JDIM,KDIM),XYJ(JDIM,KDIM)                          
C                                                                       
      DIMENSION UU(JDIM,KDIM),CC(JDIM,KDIM)                             
CX30053 CALL SBINX (53)                                                 
CY30053 CALL SBINY (53)                                                 
CX10625 CALL DOINX( 625 )                                               
CY10625 CALL DOINY( 625 )                                               
C                                                                       
C   COMPUTE SPECTRAL RADIUS SCALING                                     
C                                                                       
      DO 100 K=KBEGIN,KEND                                              
      DO 100 J=JBEGIN,JEND                                              
C        RI IS J/RHO                                                    
         RI = 1./Q(J,K,1)                                               
         U = Q(J,K,2)*RI                                                
         V = Q(J,K,3)*RI                                                
C                                                                       
C  NOTE: TIME METRICS PUT IN (SLOPPY BUT THEY WORK)                     
C                                                                       
         UU(J,K) = TMET(J,K) + U*XY(J,K,IX)+V*XY(J,K,IY)                
         CC(J,K) = SNDSP(J,K)*SQRT(XY(J,K,IX)**2+XY(J,K,IY)**2)         
C                                                                       
  100 CONTINUE                                                          
CX20626 CALL DOOUTX( 626 )                                              
CY20626 CALL DOOUTY( 626 )                                              
CX40054 CALL SBOUTX (54)                                                
CY40054 CALL SBOUTY (54)                                                
C                                                                       
                RETURN                                                  
                END                                                     
C***********************************************************************
C***************** SOLUTION UPDATE ROUTINE  ****************************
C***********************************************************************
      SUBROUTINE UPDATE(JDIM,KDIM,Q,S)                                  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
C                                                                       
      DIMENSION Q(JDIM,KDIM,4),S(JDIM,KDIM,4)                           
CX30055 CALL SBINX (55)                                                 
CY30055 CALL SBINY (55)                                                 
C                                                                       
C                                                                       
         PHIC = PHIDT/(1.+PHIDT)                                        
CX10627 CALL DOINX( 627 )                                               
CY10627 CALL DOINY( 627 )                                               
C                                                                       
C  UPDATE SOLUTION                                                      
C                                                                       
         DO 600 N = 1,4                                                 
         DO 600 K = KLOW,KUP                                            
         DO 600 J = JLOW,JUP                                            
C  DELTA Q N                                                            
C  UPDATE                                                               
            Q(J,K,N) = Q(J,K,N)+S(J,K,N)                                
C  START NEXT                                                           
            S(J,K,N) = PHIC * S(J,K,N)                                  
600      CONTINUE                                                       
CX20628 CALL DOOUTX( 628 )                                              
CY20628 CALL DOOUTY( 628 )                                              
CX40056 CALL SBOUTX (56)                                                
CY40056 CALL SBOUTY (56)                                                
C                                                                       
C                                                                       
        RETURN                                                          
        END                                                             
C***********************************************************************
C***************** N INVERSE MULTIPLY ON S **************************** 
C***********************************************************************
      SUBROUTINE NINVER(JDIM,KDIM,S,XY,IXX,IXY,IYX,IYY)                 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
C                                                                       
C   IXX = 1, IXY = 2, IYX = 3, IYY = 4                                  
C               -1                                                      
C   N  =  T    T                                                        
C          XI   ETA                                                     
C                                                                       
C   IXX = 3, IXY = 4, IYX = 1, IYY = 2                                  
C               -1                                                      
C   N  =  T    T                                                        
C          ETA  XI                                                      
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION S(JDIM,KDIM,4),XY(JDIM,KDIM,4)                          
C                                                                       
      DIMENSION RG(4,4),SR(4)                                           
CX30057 CALL SBINX (57)                                                 
CY30057 CALL SBINY (57)                                                 
C                                                                       
C  PARAMETERS FOR DIAGONALIZATION                                       
C                                                                       
      BT = 1./SQRT(2.)                                                  
      BT2 = 0.5                                                         
CX10629 CALL DOINX( 629 )                                               
CY10629 CALL DOINY( 629 )                                               
C                                                                       
        DO 1 K = KLOW,KUP                                               
          DO 1 J = JLOW,JUP                                             
           RR1   = XY(J,K,IXX)**2+XY(J,K,IXY)**2                        
           RR2   = XY(J,K,IYX)**2+XY(J,K,IYY)**2                        
           RR1   = 1./SQRT(RR1)                                         
           RR2   = 1./SQRT(RR2)                                         
           RM1   = (XY(J,K,IXX)*XY(J,K,IYX)+XY(J,K,IXY)*XY(J,K,IYY))    
     *                                          *RR1*RR2                
           RM2   = (XY(J,K,IXX)*XY(J,K,IYY)-XY(J,K,IXY)*XY(J,K,IYX))    
     *                                          *RR1*RR2                
               S34M  = S(J,K,3) - S(J,K,4)                              
               S34P  = S(J,K,3) + S(J,K,4)                              
               S2 =  RM1*S(J,K,2)    + BT*S34M*RM2                      
               S3 = -BT*S(J,K,2)*RM2 + BT2*(S34P+RM1*S34M)              
               S4 =  BT*RM2*S(J,K,2) + BT2*(S34P-RM1*S34M)              
               S(J,K,2) = S2                                            
               S(J,K,3) = S3                                            
               S(J,K,4) = S4                                            
1      CONTINUE                                                         
CX20630 CALL DOOUTX( 630 )                                              
CY20630 CALL DOOUTY( 630 )                                              
CX40058 CALL SBOUTX (58)                                                
CY40058 CALL SBOUTY (58)                                                
C                                                                       
         RETURN                                                         
         END                                                            
C***********************************************************************
C************   COMPUTE L 2 NORM OF S  , THE RESIDUAL ******************
C***********************************************************************
      SUBROUTINE RESIL2(JDIM,KDIM,S,NRES)                               
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
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,   CPUSEC,   
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,             
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,          
     1  THETAD ,    JACDT,    RESID,       IPRINT,      NPCP,           
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS          
      LOGICAL PERIDC                                                    
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK                              
C                                                                       
      DIMENSION S(JDIM,KDIM,4)                                          
CX30059 CALL SBINX (59)                                                 
CY30059 CALL SBINY (59)                                                 
C                                                                       
C       COMPUTE L2  NORM OF DENSITY  RESIDUAL ( RIGHT HAND SIDE )       
C                                                                       
         RESID = 0.                                                     
CX10631 CALL DOINX( 631 )                                               
CY10631 CALL DOINY( 631 )                                               
         DO 5 K=KLOW,KUP                                                
         DO 5 J=JLOW,JUP                                                
            RESID = RESID + S(J,K,NRES)**2                              
5        CONTINUE                                                       
CX20632 CALL DOOUTX( 632 )                                              
CY20632 CALL DOOUTY( 632 )                                              
         RESID = RESID/( (JUP-JLOW+1)*(KUP-KLOW+1))                     
         RESID = SQRT( RESID ) / DT                                     
6     CONTINUE                                                          
CX40060 CALL SBOUTX (60)                                                
CY40060 CALL SBOUTY (60)                                                
                                                                        
         RETURN                                                         
         END                                                            
C*****************************************************************      
C*****************************************************************      
C  ********  VECTORIZED IN K SCALAR PENTA SOLVER   ***************      
C*****************************************************************      
C*****************************************************************      
       SUBROUTINE XPENTA(JDIM,KDIM,A,B,C,D,E,X,Y,F,JL,JU,KL,KU)         
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
       DOUBLE PRECISION LD,LD1,LD2,LDI                                  
       DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM),   
     *           E(JDIM,KDIM),X(JDIM,KDIM),Y(JDIM,KDIM),F(JDIM,KDIM)    
CX30061 CALL SBINX (61)                                                 
CY30061 CALL SBINY (61)                                                 
C                                                                       
C      ! START FORWARD GENERATION PROCESS AND SWEEP                     
C                                                                       
        J = JL                                                          
CX10633 CALL DOINX( 633 )                                               
CY10633 CALL DOINY( 633 )                                               
        DO 1 K = KL,KU                                                  
        LD = C(J,K)                                                     
        LDI = 1./LD                                                     
        F(J,K) = F(J,K)*LDI                                             
        X(J,K) = D(J,K)*LDI                                             
        Y(J,K) = E(J,K)*LDI                                             
1       CONTINUE                                                        
CX20634 CALL DOOUTX( 634 )                                              
CY20634 CALL DOOUTY( 634 )                                              
C                                                                       
        J = JL+1                                                        
CX10635 CALL DOINX( 635 )                                               
CY10635 CALL DOINY( 635 )                                               
        DO 2 K = KL,KU                                                  
        LD1 = B(J,K)                                                    
        LD = C(J,K) - LD1*X(J-1,K)                                      
        LDI = 1./LD                                                     
        F(J,K) = (F(J,K) - LD1*F(J-1,K))*LDI                            
        X(J,K) = (D(J,K) - LD1*Y(J-1,K))*LDI                            
        Y(J,K) = E(J,K)*LDI                                             
2       CONTINUE                                                        
CX20636 CALL DOOUTX( 636 )                                              
CY20636 CALL DOOUTY( 636 )                                              
CX10637 CALL DOINX( 637 )                                               
CY10637 CALL DOINY( 637 )                                               
C                                                                       
            DO 3 J = JL+2,JU-2                                          
                DO 11 K = KL,KU                                         
                LD2 = A(J,K)                                            
                LD1 = B(J,K) - LD2*X(J-2,K)                             
                LD = C(J,K) - (LD2*Y(J-2,K) + LD1*X(J-1,K))             
                LDI = 1./LD                                             
                F(J,K) = (F(J,K) - LD2*F(J-2,K) - LD1*F(J-1,K))*LDI     
                X(J,K) = (D(J,K) - LD1*Y(J-1,K))*LDI                    
                Y(J,K) = E(J,K)*LDI                                     
11              CONTINUE                                                
3            CONTINUE                                                   
CX20638 CALL DOOUTX( 638 )                                              
CY20638 CALL DOOUTY( 638 )                                              
C                                                                       
        J = JU-1                                                        
CX10639 CALL DOINX( 639 )                                               
CY10639 CALL DOINY( 639 )                                               
        DO 12 K = KL,KU                                                 
        LD2 = A(J,K)                                                    
        LD1 = B(J,K) - LD2*X(J-2,K)                                     
        LD = C(J,K) - (LD2*Y(J-2,K) + LD1*X(J-1,K))                     
        LDI = 1./LD                                                     
        F(J,K) = (F(J,K) - LD2*F(J-2,K) - LD1*F(J-1,K))*LDI             
        X(J,K) = (D(J,K) - LD1*Y(J-1,K))*LDI                            
12      CONTINUE                                                        
CX20640 CALL DOOUTX( 640 )                                              
CY20640 CALL DOOUTY( 640 )                                              
C                                                                       
                J = JU                                                  
CX10641 CALL DOINX( 641 )                                               
CY10641 CALL DOINY( 641 )                                               
                DO 13 K = KL,KU                                         
                LD2 = A(J,K)                                            
                LD1 = B(J,K) - LD2*X(J-2,K)                             
                LD = C(J,K) - (LD2*Y(J-2,K) + LD1*X(J-1,K))             
                LDI = 1./LD                                             
                F(J,K) = (F(J,K) - LD2*F(J-2,K) - LD1*F(J-1,K))*LDI     
13              CONTINUE                                                
CX20642 CALL DOOUTX( 642 )                                              
CY20642 CALL DOOUTY( 642 )                                              
CX10643 CALL DOINX( 643 )                                               
CY10643 CALL DOINY( 643 )                                               
C                                                                       
C        !  BACK SWEEP SOLUTION                                         
C                                                                       
        DO 14 K = KL,KU                                                 
        F(JU,K) = F(JU,K)                                               
        F(JU-1,K) = F(JU-1,K) - X(JU-1,K)*F(JU,K)                       
14      CONTINUE                                                        
CX20644 CALL DOOUTX( 644 )                                              
CY20644 CALL DOOUTY( 644 )                                              
CX10645 CALL DOINX( 645 )                                               
CY10645 CALL DOINY( 645 )                                               
C                                                                       
        DO 4 J = 2,JU-JL                                                
           JX = JU-J                                                    
           DO 15 K = KL,KU                                              
           F(JX,K) = F(JX,K) - X(JX,K)*F(JX+1,K) - Y(JX,K)*F(JX+2,K)    
15         CONTINUE                                                     
4       CONTINUE                                                        
CX20646 CALL DOOUTX( 646 )                                              
CY20646 CALL DOOUTY( 646 )                                              
CX40062 CALL SBOUTX (62)                                                
CY40062 CALL SBOUTY (62)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C*****************************************************************      
C*****************************************************************      
C  ********  VECTORIZED IN K SCALAR PENTA SOLVER   ***************      
C    SOLVES FOR TWO F'S AT A TIME                  ***************      
C*****************************************************************      
C*****************************************************************      
       SUBROUTINE XPENT2(JDIM,KDIM,A,B,C,D,E,X,Y,F,JL,JU,KL,KU)         
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
       DOUBLE PRECISION LD,LD1,LD2,LDI                                  
       DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM),   
     *           E(JDIM,KDIM),X(JDIM,KDIM),Y(JDIM,KDIM),F(JDIM,KDIM,2)  
CX30063 CALL SBINX (63)                                                 
CY30063 CALL SBINY (63)                                                 
C                                                                       
C      ! START FORWARD GENERATION PROCESS AND SWEEP                     
C                                                                       
        J = JL                                                          
CX10647 CALL DOINX( 647 )                                               
CY10647 CALL DOINY( 647 )                                               
        DO 1 K = KL,KU                                                  
        LD = C(J,K)                                                     
        LDI = 1./LD                                                     
        F(J,K,1) = F(J,K,1)*LDI                                         
        F(J,K,2) = F(J,K,2)*LDI                                         
        X(J,K) = D(J,K)*LDI                                             
        Y(J,K) = E(J,K)*LDI                                             
1       CONTINUE                                                        
CX20648 CALL DOOUTX( 648 )                                              
CY20648 CALL DOOUTY( 648 )                                              
C                                                                       
        J = JL+1                                                        
CX10649 CALL DOINX( 649 )                                               
CY10649 CALL DOINY( 649 )                                               
        DO 2 K = KL,KU                                                  
        LD1 = B(J,K)                                                    
        LD = C(J,K) - LD1*X(J-1,K)                                      
        LDI = 1./LD                                                     
        F(J,K,1) = (F(J,K,1) - LD1*F(J-1,K,1))*LDI                      
        F(J,K,2) = (F(J,K,2) - LD1*F(J-1,K,2))*LDI                      
        X(J,K) = (D(J,K) - LD1*Y(J-1,K))*LDI                            
        Y(J,K) = E(J,K)*LDI                                             
2       CONTINUE                                                        
CX20650 CALL DOOUTX( 650 )                                              
CY20650 CALL DOOUTY( 650 )                                              
CX10651 CALL DOINX( 651 )                                               
CY10651 CALL DOINY( 651 )                                               
C                                                                       
            DO 3 J = JL+2,JU-2                                          
                DO 11 K = KL,KU                                         
                LD2 = A(J,K)                                            
                LD1 = B(J,K) - LD2*X(J-2,K)                             
                LD = C(J,K) - (LD2*Y(J-2,K) + LD1*X(J-1,K))             
                LDI = 1./LD                                             
                F(J,K,1) = (F(J,K,1) - LD2*F(J-2,K,1)                   
     *                      - LD1*F(J-1,K,1))*LDI                       
                F(J,K,2) = (F(J,K,2) - LD2*F(J-2,K,2)                   
     *                      - LD1*F(J-1,K,2))*LDI                       
                X(J,K) = (D(J,K) - LD1*Y(J-1,K))*LDI                    
                Y(J,K) = E(J,K)*LDI                                     
11              CONTINUE                                                
3            CONTINUE                                                   
CX20652 CALL DOOUTX( 652 )                                              
CY20652 CALL DOOUTY( 652 )                                              
C                                                                       
        J = JU-1                                                        
CX10653 CALL DOINX( 653 )                                               
CY10653 CALL DOINY( 653 )                                               
        DO 12 K = KL,KU                                                 
        LD2 = A(J,K)                                                    
        LD1 = B(J,K) - LD2*X(J-2,K)                                     
        LD = C(J,K) - (LD2*Y(J-2,K) + LD1*X(J-1,K))                     
        LDI = 1./LD                                                     
        F(J,K,1) = (F(J,K,1) - LD2*F(J-2,K,1) - LD1*F(J-1,K,1))*LDI     
        F(J,K,2) = (F(J,K,2) - LD2*F(J-2,K,2) - LD1*F(J-1,K,2))*LDI     
        X(J,K) = (D(J,K) - LD1*Y(J-1,K))*LDI                            
12      CONTINUE                                                        
CX20654 CALL DOOUTX( 654 )                                              
CY20654 CALL DOOUTY( 654 )                                              
C                                                                       
                J = JU                                                  
CX10655 CALL DOINX( 655 )                                               
CY10655 CALL DOINY( 655 )                                               
                DO 13 K = KL,KU                                         
                LD2 = A(J,K)                                            
                LD1 = B(J,K) - LD2*X(J-2,K)                             
                LD = C(J,K) - (LD2*Y(J-2,K) + LD1*X(J-1,K))             
                LDI = 1./LD                                             
                F(J,K,1) = (F(J,K,1) - LD2*F(J-2,K,1) -                 
     *                     LD1*F(J-1,K,1))*LDI                          
                F(J,K,2) = (F(J,K,2) - LD2*F(J-2,K,2) -                 
     *                      LD1*F(J-1,K,2))*LDI                         
13              CONTINUE                                                
CX20656 CALL DOOUTX( 656 )                                              
CY20656 CALL DOOUTY( 656 )                                              
CX10657 CALL DOINX( 657 )                                               
CY10657 CALL DOINY( 657 )                                               
C                                                                       
C        !  BACK SWEEP SOLUTION                                         
C                                                                       
        DO 14 K = KL,KU                                                 
        F(JU,K,1) = F(JU,K,1)                                           
        F(JU,K,2) = F(JU,K,2)                                           
        F(JU-1,K,1) = F(JU-1,K,1) - X(JU-1,K)*F(JU,K,1)                 
        F(JU-1,K,2) = F(JU-1,K,2) - X(JU-1,K)*F(JU,K,2)                 
14      CONTINUE                                                        
CX20658 CALL DOOUTX( 658 )                                              
CY20658 CALL DOOUTY( 658 )                                              
CX10659 CALL DOINX( 659 )                                               
CY10659 CALL DOINY( 659 )                                               
C                                                                       
        DO 4 J = 2,JU-JL                                                
           JX = JU-J                                                    
           DO 15 K = KL,KU                                              
           F(JX,K,1) = F(JX,K,1) - X(JX,K)*F(JX+1,K,1) -                
     *                 Y(JX,K)*F(JX+2,K,1)                              
           F(JX,K,2) = F(JX,K,2) - X(JX,K)*F(JX+1,K,2) -                
     *                 Y(JX,K)*F(JX+2,K,2)                              
15         CONTINUE                                                     
4       CONTINUE                                                        
CX20660 CALL DOOUTX( 660 )                                              
CY20660 CALL DOOUTY( 660 )                                              
CX40064 CALL SBOUTX (64)                                                
CY40064 CALL SBOUTY (64)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C*********************************************************************  
C*********************************************************************  
C************  VECTORIZED IN K SCALAR PERIDC   PENTA SOLVER  *********  
C*********************************************************************  
C*********************************************************************  
        SUBROUTINE XPENTP(JDIM,KDIM,A,B,C,D,E,U1,U2,U3,U4,F,            
     *                     JL,JU,KL,KU)                                 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      PARAMETER (MAXJ=289, MAXK=89, MAXJK=MAXJ*MAXK  )                  
        DOUBLE PRECISION LD,LD1,LD2,LDI,LDD,LU                          
        DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM)   
        DIMENSION E(JDIM,KDIM),F(JDIM,KDIM),U1(JDIM,KDIM)               
        DIMENSION U2(JDIM,KDIM),U3(JDIM,KDIM),U4(JDIM,KDIM)             
C                                                                       
        COMMON/WORKSP/LDA(MAXJ),LDB(MAXJ),LUS(MAXJ),                    
     *         LDS(MAXJ),WORK(MAXJ,88)                                  
        DOUBLE PRECISION LDA,LDB,LUS,LDS                                
CX30065 CALL SBINX (65)                                                 
CY30065 CALL SBINY (65)                                                 
C                                                                       
C        ! START FORWARD GENERATION PROCESS AND SWEEP                   
C                                                                       
        I = JL                                                          
CX10661 CALL DOINX( 661 )                                               
CY10661 CALL DOINY( 661 )                                               
        DO 10 K = KL,KU                                                 
        LD = C(I,K)                                                     
        LDI = 1./LD                                                     
        U1(I,K) = D(I,K)*LDI                                            
        U2(I,K) = E(I,K)*LDI                                            
        U3(I,K) = A(I,K)*LDI                                            
        U4(I,K) = B(I,K)*LDI                                            
        F(I,K) = F(I,K)*LDI                                             
10      CONTINUE                                                        
CX20662 CALL DOOUTX( 662 )                                              
CY20662 CALL DOOUTY( 662 )                                              
C                                                                       
        I = JL+1                                                        
CX10663 CALL DOINX( 663 )                                               
CY10663 CALL DOINY( 663 )                                               
        DO 20 K = KL,KU                                                 
        LD1 = B(I,K)                                                    
        LD = C(I,K) - LD1*U1(I-1,K)                                     
        LDI = 1./LD                                                     
        U1(I,K) = (D(I,K) - LD1*U2(I-1,K))*LDI                          
        U2(I,K) = E(I,K)*LDI                                            
        U3(I,K) = -LD1*U3(I-1,K)*LDI                                    
        U4(I,K) = (A(I,K) -LD1*U4(I-1,K))*LDI                           
        F(I,K) = (F(I,K) - LD1*F(I-1,K))*LDI                            
20      CONTINUE                                                        
CX20664 CALL DOOUTX( 664 )                                              
CY20664 CALL DOOUTY( 664 )                                              
CX10665 CALL DOINX( 665 )                                               
CY10665 CALL DOINY( 665 )                                               
C                                                                       
          DO 1 I = JL+2,JU-4                                            
                DO 11 K = KL,KU                                         
                LD2 = A(I,K)                                            
                LD1 = B(I,K) - LD2*U1(I-2,K)                            
                LD = C(I,K) - (LD2*U2(I-2,K) + LD1*U1(I-1,K))           
                LDI = 1./LD                                             
                U1(I,K) = (D(I,K) - LD1*U2(I-1,K))*LDI                  
                U2(I,K) = E(I,K)*LDI                                    
                U3(I,K) = (-LD2*U3(I-2,K) - LD1*U3(I-1,K))*LDI          
                U4(I,K) = (-LD2*U4(I-2,K) - LD1*U4(I-1,K))*LDI          
                F(I,K) = (F(I,K) - LD2*F(I-2,K) - LD1*F(I-1,K))*LDI     
11              CONTINUE                                                
1           CONTINUE                                                    
CX20666 CALL DOOUTX( 666 )                                              
CY20666 CALL DOOUTY( 666 )                                              
C                                                                       
                I = JU-3                                                
CX10667 CALL DOINX( 667 )                                               
CY10667 CALL DOINY( 667 )                                               
                DO 12 K = KL,KU                                         
                LD2 = A(I,K)                                            
                LD1 = B(I,K) - LD2*U1(I-2,K)                            
                LD = C(I,K) - (LD2*U2(I-2,K) + LD1*U1(I-1,K))           
                LDI = 1./LD                                             
                U1(I,K) = (D(I,K) - LD1*U2(I-1,K))*LDI                  
                U3(I,K) = (E(I,K) - LD2*U3(I-2,K) - LD1*U3(I-1,K))*LDI  
                U4(I,K) = (-LD2*U4(I-2,K) - LD1*U4(I-1,K))*LDI          
                F(I,K) = (F(I,K) - LD2*F(I-2,K) - LD1*F(I-1,K))*LDI     
12              CONTINUE                                                
CX20668 CALL DOOUTX( 668 )                                              
CY20668 CALL DOOUTY( 668 )                                              
C                                                                       
        I = JU-2                                                        
CX10669 CALL DOINX( 669 )                                               
CY10669 CALL DOINY( 669 )                                               
        DO 22 K = KL,KU                                                 
        LD2 = A(I,K)                                                    
        LD1 = B(I,K) - LD2*U1(I-2,K)                                    
        LD = C(I,K) - (LD2*U2(I-2,K) + LD1*U1(I-1,K))                   
        LDI = 1./LD                                                     
        U3(I,K) = (D(I,K) - LD2*U3(I-2,K) - LD1*U3(I-1,K))*LDI          
        U4(I,K) = (E(I,K) - LD2*U4(I-2,K) - LD1*U4(I-1,K))*LDI          
        F(I,K) = (F(I,K) - LD2*F(I-2,K) - LD1*F(I-1,K))*LDI             
22      CONTINUE                                                        
CX20670 CALL DOOUTX( 670 )                                              
CY20670 CALL DOOUTY( 670 )                                              
C                                                                       
        I = JU-1                                                        
CX10671 CALL DOINX( 671 )                                               
CY10671 CALL DOINY( 671 )                                               
        DO 32 K = KL,KU                                                 
        LDA(K) = E(I,K)                                                 
        LDB(K) = -LDA(K)*U1(JL,K)                                       
        F(I,K) = F(I,K) - LDA(K)*F(JL,K) - LDB(K)*F(JL+1,K)             
        LDS(K) = LDA(K)*U3(JL,K) + LDB(K)*U3(JL+1,K)                    
        U4(I,K) = LDA(K)*U4(JL,K) + LDB(K)*U4(JL+1,K)                   
32      CONTINUE                                                        
CX20672 CALL DOOUTX( 672 )                                              
CY20672 CALL DOOUTY( 672 )                                              
CX10673 CALL DOINX( 673 )                                               
CY10673 CALL DOINY( 673 )                                               
                                                                        
        DO 3 M = JL+2,JU-4                                              
           DO 14 K = KL,KU                                              
           LDD = -LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)                   
           LDA(K) = LDB(K)                                              
           LDB(K) = LDD                                                 
           LDS(K) = LDS(K) + LDD*U3(M,K)                                
           U4(I,K) = U4(I,K) + LDD*U4(M,K)                              
           F(I,K) = F(I,K) - LDD*F(M,K)                                 
14         CONTINUE                                                     
3       CONTINUE                                                        
CX20674 CALL DOOUTX( 674 )                                              
CY20674 CALL DOOUTY( 674 )                                              
C                                                                       
        M = JU-3                                                        
CX10675 CALL DOINX( 675 )                                               
CY10675 CALL DOINY( 675 )                                               
        DO 15 K = KL,KU                                                 
        LDD = A(I,K) - LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)              
        LDA(K) = LDB(K)                                                 
        LDB(K) = LDD                                                    
        LDS(K) = LDS(K) + LDD*U3(M,K)                                   
        U4(I,K) = U4(I,K) + LDD*U4(M,K)                                 
        F(I,K) = F(I,K) - LDD*F(M,K)                                    
15      CONTINUE                                                        
CX20676 CALL DOOUTX( 676 )                                              
CY20676 CALL DOOUTY( 676 )                                              
        M = JU-2                                                        
CX10677 CALL DOINX( 677 )                                               
CY10677 CALL DOINY( 677 )                                               
        DO 25 K = KL,KU                                                 
        LDD = B(I,K) - LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)              
        LDA(K) = LDB(K)                                                 
        LDB(K) = LDD                                                    
        LDS(K) = LDS(K)+ LDD*U3(M,K)                                    
        U4(I,K) = U4(I,K) + LDD*U4(M,K)                                 
        F(I,K) = F(I,K) - LDD*F(M,K)                                    
25      CONTINUE                                                        
CX20678 CALL DOOUTX( 678 )                                              
CY20678 CALL DOOUTY( 678 )                                              
CX10679 CALL DOINX( 679 )                                               
CY10679 CALL DOINY( 679 )                                               
C                                                                       
        DO 35 K = KL,KU                                                 
        LDS(K) = C(I,K) - LDS(K)                                        
        LDI = 1./LDS(K)                                                 
        U4(I,K) = (D(I,K) - U4(I,K))*LDI                                
        F(I,K) = F(I,K)*LDI                                             
35      CONTINUE                                                        
CX20680 CALL DOOUTX( 680 )                                              
CY20680 CALL DOOUTY( 680 )                                              
C                                                                       
C    THIS BREAKUP OF K LOOPS IS VERY IMPORTANT                          
C                                                                       
        I = JU                                                          
CX10681 CALL DOINX( 681 )                                               
CY10681 CALL DOINY( 681 )                                               
        DO 16 K = KL,KU                                                 
        LDA(K) = D(I,K)                                                 
        LDB(K) = E(I,K) - LDA(K)*U1(JL,K)                               
        LUS(K) = LDA(K)*U3(JL,K) + LDB(K)*U3(JL+1,K)                    
        LDS(K) = LDA(K)*U4(JL,K) + LDB(K)*U4(JL+1,K)                    
        F(I,K) = F(I,K) - LDA(K)*F(JL,K) - LDB(K)*F(JL+1,K)             
C                                                                       
16      CONTINUE                                                        
CX20682 CALL DOOUTX( 682 )                                              
CY20682 CALL DOOUTY( 682 )                                              
CX10683 CALL DOINX( 683 )                                               
CY10683 CALL DOINY( 683 )                                               
C                                                                       
        DO 4 M = JL+2,JU-3                                              
           DO 17 K = KL,KU                                              
           LDD = -LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)                   
           LDA(K) = LDB(K)                                              
           LDB(K) = LDD                                                 
           LUS(K) = LUS(K) + LDD*U3(M,K)                                
           LDS(K) = LDS(K) + LDD*U4(M,K)                                
           F(I,K) = F(I,K) - LDD*F(M,K)                                 
17         CONTINUE                                                     
4       CONTINUE                                                        
CX20684 CALL DOOUTX( 684 )                                              
CY20684 CALL DOOUTY( 684 )                                              
C                                                                       
        M = JU-2                                                        
CX10685 CALL DOINX( 685 )                                               
CY10685 CALL DOINY( 685 )                                               
        DO 18 K = KL,KU                                                 
        LDD = A(I,K) - LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)              
        LDA(K) = LDB(K)                                                 
        LDB(K) = LDD                                                    
        LUS(K) = LUS(K) + LDD*U3(M,K)                                   
        LDS(K) = LDS(K) + LDD*U4(M,K)                                   
        F(I,K) = F(I,K) - LDD*F(M,K)                                    
18      CONTINUE                                                        
CX20686 CALL DOOUTX( 686 )                                              
CY20686 CALL DOOUTY( 686 )                                              
        M = JU-1                                                        
CX10687 CALL DOINX( 687 )                                               
CY10687 CALL DOINY( 687 )                                               
        DO 28 K = KL,KU                                                 
        LDD = B(I,K) - LUS(K)                                           
        F(I,K) = F(I,K) - LDD*F(M,K)                                    
        LDS(K) = C(I,K) - LDD*U4(M,K) - LDS(K)                          
        LDI = 1./LDS(K)                                                 
        F(I,K) = F(I,K)*LDI                                             
28      CONTINUE                                                        
CX20688 CALL DOOUTX( 688 )                                              
CY20688 CALL DOOUTY( 688 )                                              
CX10689 CALL DOINX( 689 )                                               
CY10689 CALL DOINY( 689 )                                               
C                                                                       
C        !  BACK SWEEP SOLUTION                                         
C                                                                       
      DO 19 K = KL,KU                                                   
      F(JU,K) = F(JU,K)                                                 
      F(JU-1,K) = F(JU-1,K) - U4(JU-1,K)*F(JU,K)                        
      F(JU-2,K) = F(JU-2,K) - U4(JU-2,K)*F(JU,K) - U3(JU-2,K)*F(JU-1,K) 
      F(JU-3,K) = F(JU-3,K) - U4(JU-3,K)*F(JU,K) - U3(JU-3,K)*F(JU-1,K) 
     *                - U1(JU-3,K)*F(JU-2,K)                            
19      CONTINUE                                                        
CX20690 CALL DOOUTX( 690 )                                              
CY20690 CALL DOOUTY( 690 )                                              
CX10691 CALL DOINX( 691 )                                               
CY10691 CALL DOINY( 691 )                                               
C                                                                       
        DO 2 I = 4,JU-JL                                                
        IX = JU-I                                                       
          DO 50 K = KL,KU                                               
          F(IX,K) = F(IX,K) - U4(IX,K)*F(JU,K) - U3(IX,K)*F(JU-1,K)     
     *              - U1(IX,K)*F(IX+1,K) - U2(IX,K)*F(IX+2,K)           
50        CONTINUE                                                      
2       CONTINUE                                                        
CX20692 CALL DOOUTX( 692 )                                              
CY20692 CALL DOOUTY( 692 )                                              
CX40066 CALL SBOUTX (66)                                                
CY40066 CALL SBOUTY (66)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C*********************************************************************  
C*********************************************************************  
C************  VECTORIZED IN K SCALAR PERIDC   PENTA SOLVER  *********  
C   DOES TWO F'S AT A TIME                                   *********  
C*********************************************************************  
C*********************************************************************  
        SUBROUTINE XPENP2(JDIM,KDIM,A,B,C,D,E,U1,U2,U3,U4,F,            
     *                      JL,JU,KL,KU)                                
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      PARAMETER (MAXJ=289, MAXK=89, MAXJK=MAXJ*MAXK  )                  
        DOUBLE PRECISION LD,LD1,LD2,LDI,LDD,LU                          
        DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM)   
        DIMENSION E(JDIM,KDIM),F(JDIM,KDIM,2),U4(JDIM,KDIM)             
        DIMENSION U1(JDIM,KDIM),U2(JDIM,KDIM),U3(JDIM,KDIM)             
C                                                                       
        COMMON/WORKSP/LDA(MAXJ),LDB(MAXJ),LUS(MAXJ),                    
     *         LDS(MAXJ),WORK(MAXJ,88)                                  
        DOUBLE PRECISION LDA,LDB,LUS,LDS                                
CX30067 CALL SBINX (67)                                                 
CY30067 CALL SBINY (67)                                                 
C                                                                       
C        ! START FORWARD GENERATION PROCESS AND SWEEP                   
C                                                                       
        I = JL                                                          
CX10693 CALL DOINX( 693 )                                               
CY10693 CALL DOINY( 693 )                                               
        DO 10 K = KL,KU                                                 
        LD = C(I,K)                                                     
        LDI = 1./LD                                                     
        U1(I,K) = D(I,K)*LDI                                            
        U2(I,K) = E(I,K)*LDI                                            
        U3(I,K) = A(I,K)*LDI                                            
        U4(I,K) = B(I,K)*LDI                                            
        F(I,K,1) = F(I,K,1)*LDI                                         
        F(I,K,2) = F(I,K,2)*LDI                                         
10      CONTINUE                                                        
CX20694 CALL DOOUTX( 694 )                                              
CY20694 CALL DOOUTY( 694 )                                              
C                                                                       
        I = JL+1                                                        
CX10695 CALL DOINX( 695 )                                               
CY10695 CALL DOINY( 695 )                                               
        DO 20 K = KL,KU                                                 
        LD1 = B(I,K)                                                    
        LD = C(I,K) - LD1*U1(I-1,K)                                     
        LDI = 1./LD                                                     
        U1(I,K) = (D(I,K) - LD1*U2(I-1,K))*LDI                          
        U2(I,K) = E(I,K)*LDI                                            
        U3(I,K) = -LD1*U3(I-1,K)*LDI                                    
        U4(I,K) = (A(I,K) -LD1*U4(I-1,K))*LDI                           
        F(I,K,1) = (F(I,K,1) - LD1*F(I-1,K,1))*LDI                      
        F(I,K,2) = (F(I,K,2) - LD1*F(I-1,K,2))*LDI                      
20      CONTINUE                                                        
CX20696 CALL DOOUTX( 696 )                                              
CY20696 CALL DOOUTY( 696 )                                              
CX10697 CALL DOINX( 697 )                                               
CY10697 CALL DOINY( 697 )                                               
C                                                                       
          DO 1 I = JL+2,JU-4                                            
                DO 11 K = KL,KU                                         
                LD2 = A(I,K)                                            
                LD1 = B(I,K) - LD2*U1(I-2,K)                            
                LD = C(I,K) - (LD2*U2(I-2,K) + LD1*U1(I-1,K))           
                LDI = 1./LD                                             
                U1(I,K) = (D(I,K) - LD1*U2(I-1,K))*LDI                  
                U2(I,K) = E(I,K)*LDI                                    
                U3(I,K) = (-LD2*U3(I-2,K) - LD1*U3(I-1,K))*LDI          
                U4(I,K) = (-LD2*U4(I-2,K) - LD1*U4(I-1,K))*LDI          
                F(I,K,1) = (F(I,K,1) - LD2*F(I-2,K,1) -                 
     *                     LD1*F(I-1,K,1))*LDI                          
                F(I,K,2) = (F(I,K,2) - LD2*F(I-2,K,2) -                 
     *                      LD1*F(I-1,K,2))*LDI                         
11              CONTINUE                                                
1           CONTINUE                                                    
CX20698 CALL DOOUTX( 698 )                                              
CY20698 CALL DOOUTY( 698 )                                              
C                                                                       
                I = JU-3                                                
CX10699 CALL DOINX( 699 )                                               
CY10699 CALL DOINY( 699 )                                               
                DO 12 K = KL,KU                                         
                LD2 = A(I,K)                                            
                LD1 = B(I,K) - LD2*U1(I-2,K)                            
                LD = C(I,K) - (LD2*U2(I-2,K) + LD1*U1(I-1,K))           
                LDI = 1./LD                                             
                U1(I,K) = (D(I,K) - LD1*U2(I-1,K))*LDI                  
                U3(I,K) = (E(I,K) - LD2*U3(I-2,K) - LD1*U3(I-1,K))*LDI  
                U4(I,K) = (-LD2*U4(I-2,K) - LD1*U4(I-1,K))*LDI          
                F(I,K,1) = (F(I,K,1) - LD2*F(I-2,K,1) -                 
     *                     LD1*F(I-1,K,1))*LDI                          
                F(I,K,2) = (F(I,K,2) - LD2*F(I-2,K,2) -                 
     *                     LD1*F(I-1,K,2))*LDI                          
12              CONTINUE                                                
CX20700 CALL DOOUTX( 700 )                                              
CY20700 CALL DOOUTY( 700 )                                              
C                                                                       
        I = JU-2                                                        
CX10701 CALL DOINX( 701 )                                               
CY10701 CALL DOINY( 701 )                                               
        DO 22 K = KL,KU                                                 
        LD2 = A(I,K)                                                    
        LD1 = B(I,K) - LD2*U1(I-2,K)                                    
        LD = C(I,K) - (LD2*U2(I-2,K) + LD1*U1(I-1,K))                   
        LDI = 1./LD                                                     
        U3(I,K) = (D(I,K) - LD2*U3(I-2,K) - LD1*U3(I-1,K))*LDI          
        U4(I,K) = (E(I,K) - LD2*U4(I-2,K) - LD1*U4(I-1,K))*LDI          
        F(I,K,1) = (F(I,K,1) - LD2*F(I-2,K,1) - LD1*F(I-1,K,1))*LDI     
        F(I,K,2) = (F(I,K,2) - LD2*F(I-2,K,2) - LD1*F(I-1,K,2))*LDI     
22      CONTINUE                                                        
CX20702 CALL DOOUTX( 702 )                                              
CY20702 CALL DOOUTY( 702 )                                              
C                                                                       
        I = JU-1                                                        
CX10703 CALL DOINX( 703 )                                               
CY10703 CALL DOINY( 703 )                                               
        DO 32 K = KL,KU                                                 
        LDA(K) = E(I,K)                                                 
        LDB(K) = -LDA(K)*U1(JL,K)                                       
        F(I,K,1) = F(I,K,1) - LDA(K)*F(JL,K,1) - LDB(K)*F(JL+1,K,1)     
        F(I,K,2) = F(I,K,2) - LDA(K)*F(JL,K,2) - LDB(K)*F(JL+1,K,2)     
        LDS(K) = LDA(K)*U3(JL,K) + LDB(K)*U3(JL+1,K)                    
        U4(I,K) = LDA(K)*U4(JL,K) + LDB(K)*U4(JL+1,K)                   
32      CONTINUE                                                        
CX20704 CALL DOOUTX( 704 )                                              
CY20704 CALL DOOUTY( 704 )                                              
CX10705 CALL DOINX( 705 )                                               
CY10705 CALL DOINY( 705 )                                               
                                                                        
        DO 3 M = JL+2,JU-4                                              
           DO 14 K = KL,KU                                              
           LDD = -LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)                   
           LDA(K) = LDB(K)                                              
           LDB(K) = LDD                                                 
           LDS(K) = LDS(K) + LDD*U3(M,K)                                
           U4(I,K) = U4(I,K) + LDD*U4(M,K)                              
           F(I,K,1) = F(I,K,1) - LDD*F(M,K,1)                           
           F(I,K,2) = F(I,K,2) - LDD*F(M,K,2)                           
14         CONTINUE                                                     
3       CONTINUE                                                        
CX20706 CALL DOOUTX( 706 )                                              
CY20706 CALL DOOUTY( 706 )                                              
C                                                                       
        M = JU-3                                                        
CX10707 CALL DOINX( 707 )                                               
CY10707 CALL DOINY( 707 )                                               
        DO 15 K = KL,KU                                                 
        LDD = A(I,K) - LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)              
        LDA(K) = LDB(K)                                                 
        LDB(K) = LDD                                                    
        LDS(K) = LDS(K) + LDD*U3(M,K)                                   
        U4(I,K) = U4(I,K) + LDD*U4(M,K)                                 
        F(I,K,1) = F(I,K,1) - LDD*F(M,K,1)                              
        F(I,K,2) = F(I,K,2) - LDD*F(M,K,2)                              
15      CONTINUE                                                        
CX20708 CALL DOOUTX( 708 )                                              
CY20708 CALL DOOUTY( 708 )                                              
        M = JU-2                                                        
CX10709 CALL DOINX( 709 )                                               
CY10709 CALL DOINY( 709 )                                               
        DO 25 K = KL,KU                                                 
        LDD = B(I,K) - LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)              
        LDA(K) = LDB(K)                                                 
        LDB(K) = LDD                                                    
        LDS(K) = LDS(K)+ LDD*U3(M,K)                                    
        U4(I,K) = U4(I,K) + LDD*U4(M,K)                                 
        F(I,K,1) = F(I,K,1) - LDD*F(M,K,1)                              
        F(I,K,2) = F(I,K,2) - LDD*F(M,K,2)                              
25      CONTINUE                                                        
CX20710 CALL DOOUTX( 710 )                                              
CY20710 CALL DOOUTY( 710 )                                              
CX10711 CALL DOINX( 711 )                                               
CY10711 CALL DOINY( 711 )                                               
C                                                                       
        DO 35 K = KL,KU                                                 
        LDS(K) = C(I,K) - LDS(K)                                        
        LDI = 1./LDS(K)                                                 
        U4(I,K) = (D(I,K) - U4(I,K))*LDI                                
        F(I,K,1) = F(I,K,1)*LDI                                         
        F(I,K,2) = F(I,K,2)*LDI                                         
35      CONTINUE                                                        
CX20712 CALL DOOUTX( 712 )                                              
CY20712 CALL DOOUTY( 712 )                                              
C                                                                       
C    THIS BREAKUP OF K LOOPS IS VERY IMPORTANT                          
C                                                                       
        I = JU                                                          
CX10713 CALL DOINX( 713 )                                               
CY10713 CALL DOINY( 713 )                                               
        DO 16 K = KL,KU                                                 
        LDA(K) = D(I,K)                                                 
        LDB(K) = E(I,K) - LDA(K)*U1(JL,K)                               
        LUS(K) = LDA(K)*U3(JL,K) + LDB(K)*U3(JL+1,K)                    
        LDS(K) = LDA(K)*U4(JL,K) + LDB(K)*U4(JL+1,K)                    
        F(I,K,1) = F(I,K,1) - LDA(K)*F(JL,K,1) - LDB(K)*F(JL+1,K,1)     
        F(I,K,2) = F(I,K,2) - LDA(K)*F(JL,K,2) - LDB(K)*F(JL+1,K,2)     
C                                                                       
16      CONTINUE                                                        
CX20714 CALL DOOUTX( 714 )                                              
CY20714 CALL DOOUTY( 714 )                                              
CX10715 CALL DOINX( 715 )                                               
CY10715 CALL DOINY( 715 )                                               
C                                                                       
        DO 4 M = JL+2,JU-3                                              
           DO 17 K = KL,KU                                              
           LDD = -LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)                   
           LDA(K) = LDB(K)                                              
           LDB(K) = LDD                                                 
           LUS(K) = LUS(K) + LDD*U3(M,K)                                
           LDS(K) = LDS(K) + LDD*U4(M,K)                                
           F(I,K,1) = F(I,K,1) - LDD*F(M,K,1)                           
           F(I,K,2) = F(I,K,2) - LDD*F(M,K,2)                           
17         CONTINUE                                                     
4       CONTINUE                                                        
CX20716 CALL DOOUTX( 716 )                                              
CY20716 CALL DOOUTY( 716 )                                              
C                                                                       
        M = JU-2                                                        
CX10717 CALL DOINX( 717 )                                               
CY10717 CALL DOINY( 717 )                                               
        DO 18 K = KL,KU                                                 
        LDD = A(I,K) - LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)              
        LDA(K) = LDB(K)                                                 
        LDB(K) = LDD                                                    
        LUS(K) = LUS(K) + LDD*U3(M,K)                                   
        LDS(K) = LDS(K) + LDD*U4(M,K)                                   
        F(I,K,1) = F(I,K,1) - LDD*F(M,K,1)                              
        F(I,K,2) = F(I,K,2) - LDD*F(M,K,2)                              
18      CONTINUE                                                        
CX20718 CALL DOOUTX( 718 )                                              
CY20718 CALL DOOUTY( 718 )                                              
        M = JU-1                                                        
CX10719 CALL DOINX( 719 )                                               
CY10719 CALL DOINY( 719 )                                               
        DO 28 K = KL,KU                                                 
        LDD = B(I,K) - LUS(K)                                           
        F(I,K,1) = F(I,K,1) - LDD*F(M,K,1)                              
        F(I,K,2) = F(I,K,2) - LDD*F(M,K,2)                              
        LDS(K) = C(I,K) - LDD*U4(M,K) - LDS(K)                          
        LDI = 1./LDS(K)                                                 
        F(I,K,1) = F(I,K,1)*LDI                                         
        F(I,K,2) = F(I,K,2)*LDI                                         
28      CONTINUE                                                        
CX20720 CALL DOOUTX( 720 )                                              
CY20720 CALL DOOUTY( 720 )                                              
CX10721 CALL DOINX( 721 )                                               
CY10721 CALL DOINY( 721 )                                               
C                                                                       
C        !  BACK SWEEP SOLUTION                                         
C                                                                       
      DO 19 K = KL,KU                                                   
      F(JU,K,1) = F(JU,K,1)                                             
      F(JU,K,2) = F(JU,K,2)                                             
      F(JU-1,K,1) = F(JU-1,K,1) - U4(JU-1,K)*F(JU,K,1)                  
      F(JU-1,K,2) = F(JU-1,K,2) - U4(JU-1,K)*F(JU,K,2)                  
      F(JU-2,K,1) = F(JU-2,K,1) - U4(JU-2,K)*F(JU,K,1) -                
     *              U3(JU-2,K)*F(JU-1,K,1)                              
      F(JU-2,K,2) = F(JU-2,K,2) - U4(JU-2,K)*F(JU,K,2) -                
     *              U3(JU-2,K)*F(JU-1,K,2)                              
      F(JU-3,K,1) = F(JU-3,K,1) - U4(JU-3,K)*F(JU,K,1) -                
     *              U3(JU-3,K)*F(JU-1,K,1)- U1(JU-3,K)*F(JU-2,K,1)      
      F(JU-3,K,2) = F(JU-3,K,2) - U4(JU-3,K)*F(JU,K,2) -                
     *              U3(JU-3,K)*F(JU-1,K,2)- U1(JU-3,K)*F(JU-2,K,2)      
19      CONTINUE                                                        
CX20722 CALL DOOUTX( 722 )                                              
CY20722 CALL DOOUTY( 722 )                                              
CX10723 CALL DOINX( 723 )                                               
CY10723 CALL DOINY( 723 )                                               
C                                                                       
        DO 2 I = 4,JU-JL                                                
        IX = JU-I                                                       
          DO 50 K = KL,KU                                               
          F(IX,K,1) = F(IX,K,1) - U4(IX,K)*F(JU,K,1) -                  
     *                U3(IX,K)*F(JU-1,K,1)- U1(IX,K)*F(IX+1,K,1) -      
     *                U2(IX,K)*F(IX+2,K,1)                              
          F(IX,K,2) = F(IX,K,2) - U4(IX,K)*F(JU,K,2) -                  
     *                U3(IX,K)*F(JU-1,K,2)- U1(IX,K)*F(IX+1,K,2) -      
     *                U2(IX,K)*F(IX+2,K,2)                              
50        CONTINUE                                                      
2       CONTINUE                                                        
CX20724 CALL DOOUTX( 724 )                                              
CY20724 CALL DOOUTY( 724 )                                              
CX40068 CALL SBOUTX (68)                                                
CY40068 CALL SBOUTY (68)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C*****************************************************************      
C*****************************************************************      
C  ********  VECTORIZED IN J SCALAR PENTA SOLVER   ***************      
C*****************************************************************      
C*****************************************************************      
       SUBROUTINE YPENTA(JDIM,KDIM,A,B,C,D,E,X,Y,F,KL,KU,JL,JU)         
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
       DOUBLE PRECISION LD,LD1,LD2,LDI                                  
       DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM),   
     *           E(JDIM,KDIM),X(JDIM,KDIM),Y(JDIM,KDIM),F(JDIM,KDIM)    
CX30069 CALL SBINX (69)                                                 
CY30069 CALL SBINY (69)                                                 
C                                                                       
C      ! START FORWARD GENERATION PROCESS AND SWEEP                     
C                                                                       
        K = KL                                                          
CX10725 CALL DOINX( 725 )                                               
CY10725 CALL DOINY( 725 )                                               
        DO 10 J = JL,JU                                                 
        LD = C(J,K)                                                     
        LDI = 1./LD                                                     
        F(J,K) = F(J,K)*LDI                                             
        X(J,K) = D(J,K)*LDI                                             
        Y(J,K) = E(J,K)*LDI                                             
10      CONTINUE                                                        
CX20726 CALL DOOUTX( 726 )                                              
CY20726 CALL DOOUTY( 726 )                                              
C                                                                       
        K = KL+1                                                        
CX10727 CALL DOINX( 727 )                                               
CY10727 CALL DOINY( 727 )                                               
        DO 20 J = JL,JU                                                 
        LD1 = B(J,K)                                                    
        LD = C(J,K) - LD1*X(J,K-1)                                      
        LDI = 1./LD                                                     
        F(J,K) = (F(J,K) - LD1*F(J,K-1))*LDI                            
        X(J,K) = (D(J,K) - LD1*Y(J,K-1))*LDI                            
        Y(J,K) = E(J,K)*LDI                                             
20      CONTINUE                                                        
CX20728 CALL DOOUTX( 728 )                                              
CY20728 CALL DOOUTY( 728 )                                              
CX10729 CALL DOINX( 729 )                                               
CY10729 CALL DOINY( 729 )                                               
C                                                                       
            DO 1 K = KL+2,KU-2                                          
                DO 11 J = JL,JU                                         
                LD2 = A(J,K)                                            
                LD1 = B(J,K) - LD2*X(J,K-2)                             
                LD = C(J,K) - (LD2*Y(J,K-2) + LD1*X(J,K-1))             
                LDI = 1./LD                                             
                F(J,K) = (F(J,K) - LD2*F(J,K-2) - LD1*F(J,K-1))*LDI     
                X(J,K) = (D(J,K) - LD1*Y(J,K-1))*LDI                    
                Y(J,K) = E(J,K)*LDI                                     
11              CONTINUE                                                
1            CONTINUE                                                   
CX20730 CALL DOOUTX( 730 )                                              
CY20730 CALL DOOUTY( 730 )                                              
C                                                                       
        K = KU-1                                                        
CX10731 CALL DOINX( 731 )                                               
CY10731 CALL DOINY( 731 )                                               
        DO 12 J = JL,JU                                                 
        LD2 = A(J,K)                                                    
        LD1 = B(J,K) - LD2*X(J,K-2)                                     
        LD = C(J,K) - (LD2*Y(J,K-2) + LD1*X(J,K-1))                     
        LDI = 1./LD                                                     
        F(J,K) = (F(J,K) - LD2*F(J,K-2) - LD1*F(J,K-1))*LDI             
        X(J,K) = (D(J,K) - LD1*Y(J,K-1))*LDI                            
12      CONTINUE                                                        
CX20732 CALL DOOUTX( 732 )                                              
CY20732 CALL DOOUTY( 732 )                                              
C                                                                       
                K = KU                                                  
CX10733 CALL DOINX( 733 )                                               
CY10733 CALL DOINY( 733 )                                               
                DO 22 J = JL,JU                                         
                LD2 = A(J,K)                                            
                LD1 = B(J,K) - LD2*X(J,K-2)                             
                LD = C(J,K) - (LD2*Y(J,K-2) + LD1*X(J,K-1))             
                LDI = 1./LD                                             
                F(J,K) = (F(J,K) - LD2*F(J,K-2) - LD1*F(J,K-1))*LDI     
22              CONTINUE                                                
CX20734 CALL DOOUTX( 734 )                                              
CY20734 CALL DOOUTY( 734 )                                              
CX10735 CALL DOINX( 735 )                                               
CY10735 CALL DOINY( 735 )                                               
C                                                                       
C        !  BACK SWEEP SOLUTION                                         
C                                                                       
        DO 32 J = JL,JU                                                 
        F(J,KU) = F(J,KU)                                               
        F(J,KU-1) = F(J,KU-1) - X(J,KU-1)*F(J,KU)                       
32      CONTINUE                                                        
CX20736 CALL DOOUTX( 736 )                                              
CY20736 CALL DOOUTY( 736 )                                              
CX10737 CALL DOINX( 737 )                                               
CY10737 CALL DOINY( 737 )                                               
C                                                                       
        DO 2 K = 2,KU-KL                                                
           KX = KU-K                                                    
           DO 13 J = JL,JU                                              
           F(J,KX) = F(J,KX) - X(J,KX)*F(J,KX+1) - Y(J,KX)*F(J,KX+2)    
13         CONTINUE                                                     
2       CONTINUE                                                        
CX20738 CALL DOOUTX( 738 )                                              
CY20738 CALL DOOUTY( 738 )                                              
CX40070 CALL SBOUTX (70)                                                
CY40070 CALL SBOUTY (70)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C*****************************************************************      
C*****************************************************************      
C  ********  VECTORIZED IN J SCALAR PENTA SOLVER   ***************      
C  TWO F'S DONE AT A TIME                          ***************      
C*****************************************************************      
C*****************************************************************      
       SUBROUTINE YPENT2(JDIM,KDIM,A,B,C,D,E,X,Y,F,KL,KU,JL,JU)         
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
       DOUBLE PRECISION LD,LD1,LD2,LDI                                  
       DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM),   
     *           E(JDIM,KDIM),X(JDIM,KDIM),Y(JDIM,KDIM),F(JDIM,KDIM,2)  
CX30071 CALL SBINX (71)                                                 
CY30071 CALL SBINY (71)                                                 
C                                                                       
C      ! START FORWARD GENERATION PROCESS AND SWEEP                     
C                                                                       
        K = KL                                                          
CX10739 CALL DOINX( 739 )                                               
CY10739 CALL DOINY( 739 )                                               
        DO 10 J = JL,JU                                                 
        LD = C(J,K)                                                     
        LDI = 1./LD                                                     
        F(J,K,1) = F(J,K,1)*LDI                                         
        F(J,K,2) = F(J,K,2)*LDI                                         
        X(J,K) = D(J,K)*LDI                                             
        Y(J,K) = E(J,K)*LDI                                             
10      CONTINUE                                                        
CX20740 CALL DOOUTX( 740 )                                              
CY20740 CALL DOOUTY( 740 )                                              
C                                                                       
        K = KL+1                                                        
CX10741 CALL DOINX( 741 )                                               
CY10741 CALL DOINY( 741 )                                               
        DO 20 J = JL,JU                                                 
        LD1 = B(J,K)                                                    
        LD = C(J,K) - LD1*X(J,K-1)                                      
        LDI = 1./LD                                                     
        F(J,K,1) = (F(J,K,1) - LD1*F(J,K-1,1))*LDI                      
        F(J,K,2) = (F(J,K,2) - LD1*F(J,K-1,2))*LDI                      
        X(J,K) = (D(J,K) - LD1*Y(J,K-1))*LDI                            
        Y(J,K) = E(J,K)*LDI                                             
20      CONTINUE                                                        
CX20742 CALL DOOUTX( 742 )                                              
CY20742 CALL DOOUTY( 742 )                                              
CX10743 CALL DOINX( 743 )                                               
CY10743 CALL DOINY( 743 )                                               
C                                                                       
            DO 1 K = KL+2,KU-2                                          
                DO 11 J = JL,JU                                         
                LD2 = A(J,K)                                            
                LD1 = B(J,K) - LD2*X(J,K-2)                             
                LD = C(J,K) - (LD2*Y(J,K-2) + LD1*X(J,K-1))             
                LDI = 1./LD                                             
                F(J,K,1) = (F(J,K,1) - LD2*F(J,K-2,1) -                 
     *                     LD1*F(J,K-1,1))*LDI                          
                F(J,K,2) = (F(J,K,2) - LD2*F(J,K-2,2) -                 
     *                     LD1*F(J,K-1,2))*LDI                          
                X(J,K) = (D(J,K) - LD1*Y(J,K-1))*LDI                    
                Y(J,K) = E(J,K)*LDI                                     
11              CONTINUE                                                
1            CONTINUE                                                   
CX20744 CALL DOOUTX( 744 )                                              
CY20744 CALL DOOUTY( 744 )                                              
C                                                                       
        K = KU-1                                                        
CX10745 CALL DOINX( 745 )                                               
CY10745 CALL DOINY( 745 )                                               
        DO 12 J = JL,JU                                                 
        LD2 = A(J,K)                                                    
        LD1 = B(J,K) - LD2*X(J,K-2)                                     
        LD = C(J,K) - (LD2*Y(J,K-2) + LD1*X(J,K-1))                     
        LDI = 1./LD                                                     
        F(J,K,1) = (F(J,K,1) - LD2*F(J,K-2,1) - LD1*F(J,K-1,1))*LDI     
        F(J,K,2) = (F(J,K,2) - LD2*F(J,K-2,2) - LD1*F(J,K-1,2))*LDI     
        X(J,K) = (D(J,K) - LD1*Y(J,K-1))*LDI                            
12      CONTINUE                                                        
CX20746 CALL DOOUTX( 746 )                                              
CY20746 CALL DOOUTY( 746 )                                              
C                                                                       
                K = KU                                                  
CX10747 CALL DOINX( 747 )                                               
CY10747 CALL DOINY( 747 )                                               
                DO 22 J = JL,JU                                         
                LD2 = A(J,K)                                            
                LD1 = B(J,K) - LD2*X(J,K-2)                             
                LD = C(J,K) - (LD2*Y(J,K-2) + LD1*X(J,K-1))             
                LDI = 1./LD                                             
                F(J,K,1) = (F(J,K,1) - LD2*F(J,K-2,1) -                 
     *                     LD1*F(J,K-1,1))*LDI                          
                F(J,K,2) = (F(J,K,2) - LD2*F(J,K-2,2) -                 
     *                      LD1*F(J,K-1,2))*LDI                         
22              CONTINUE                                                
CX20748 CALL DOOUTX( 748 )                                              
CY20748 CALL DOOUTY( 748 )                                              
CX10749 CALL DOINX( 749 )                                               
CY10749 CALL DOINY( 749 )                                               
C                                                                       
C        !  BACK SWEEP SOLUTION                                         
C                                                                       
        DO 32 J = JL,JU                                                 
        F(J,KU,1) = F(J,KU,1)                                           
        F(J,KU,2) = F(J,KU,2)                                           
        F(J,KU-1,1) = F(J,KU-1,1) - X(J,KU-1)*F(J,KU,1)                 
        F(J,KU-1,2) = F(J,KU-1,2) - X(J,KU-1)*F(J,KU,2)                 
32      CONTINUE                                                        
CX20750 CALL DOOUTX( 750 )                                              
CY20750 CALL DOOUTY( 750 )                                              
CX10751 CALL DOINX( 751 )                                               
CY10751 CALL DOINY( 751 )                                               
C                                                                       
        DO 2 K = 2,KU-KL                                                
           KX = KU-K                                                    
           DO 13 J = JL,JU                                              
           F(J,KX,1) = F(J,KX,1) - X(J,KX)*F(J,KX+1,1) -                
     *                 Y(J,KX)*F(J,KX+2,1)                              
           F(J,KX,2) = F(J,KX,2) - X(J,KX)*F(J,KX+1,2) -                
     *                 Y(J,KX)*F(J,KX+2,2)                              
13         CONTINUE                                                     
2       CONTINUE                                                        
CX20752 CALL DOOUTX( 752 )                                              
CY20752 CALL DOOUTY( 752 )                                              
CX40072 CALL SBOUTX (72)                                                
CY40072 CALL SBOUTY (72)                                                
C                                                                       
        RETURN                                                          
        END                                                             
C***********************************************************************
C*************** SCALAR TRIDIAGONAL ************************************
C***********************************************************************
      SUBROUTINE TRIB(JDIM,A,B,C,X,F,NL,NU)                             
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
C                                                                       
      DIMENSION A(JDIM),B(JDIM),C(JDIM),X(JDIM),F(JDIM)                 
CX30073 CALL SBINX (73)                                                 
CY30073 CALL SBINY (73)                                                 
C                                                                       
      X(NL)=C(NL)/B(NL)                                                 
      F(NL)=F(NL)/B(NL)                                                 
      NLP1 = NL +1                                                      
CX10753 CALL DOINX( 753 )                                               
CY10753 CALL DOINY( 753 )                                               
      DO 1 J=NLP1,NU                                                    
         Z=1./(B(J)-A(J)*X(J-1))                                        
         X(J)=C(J)*Z                                                    
         F(J)=(F(J)-A(J)*F(J-1))*Z                                      
1     CONTINUE                                                          
CX20754 CALL DOOUTX( 754 )                                              
CY20754 CALL DOOUTY( 754 )                                              
C                                                                       
      NUPNL=NU+NL                                                       
CX10755 CALL DOINX( 755 )                                               
CY10755 CALL DOINY( 755 )                                               
      DO 2 J1=NLP1,NU                                                   
         J=NUPNL-J1                                                     
         F(J)=F(J)-X(J)*F(J+1)                                          
2     CONTINUE                                                          
CX20756 CALL DOOUTX( 756 )                                              
CY20756 CALL DOOUTY( 756 )                                              
CX40074 CALL SBOUTX (74)                                                
CY40074 CALL SBOUTY (74)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C*************** SCALAR PERIDC   TRIDIAGONAL ***************************
C***********************************************************************
      SUBROUTINE TRIP(JDIM,A,B,C,F,Q,S,J1,J2)                           
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION A(JDIM),B(JDIM),C(JDIM),F(JDIM),Q(JDIM),S(JDIM)         
CX30075 CALL SBINX (75)                                                 
CY30075 CALL SBINY (75)                                                 
      JA = J1 + 1                                                       
      FN = F(J2)                                                        
C                                                                       
C         FORWARD ELIMINATION SWEEP                                     
C                                                                       
      Q(J1) = -C(J1)/B(J1)                                              
      F(J1) = F(J1)/B(J1)                                               
      S(J1) = - A(J1)/B(J1)                                             
CX10757 CALL DOINX( 757 )                                               
CY10757 CALL DOINY( 757 )                                               
      DO 10 J=JA,J2                                                     
         P =1./( B(J) + A(J)*Q(J-1))                                    
         Q(J) = - C(J)*P                                                
         F(J) = ( F(J) - A(J)*F(J-1))*P                                 
         S(J) = - A(J)*S(J-1)*P                                         
   10 CONTINUE                                                          
C                                                                       
C         BACKWARD PASS                                                 
C                                                                       
      JJ = J1 + J2                                                      
      Q(J2) = 0.                                                        
      S(J2) = 1.                                                        
      DO 11 I=JA,J2                                                     
         J = JJ - I                                                     
         S(J) = S(J) + Q(J)*S(J+1)                                      
         Q(J) = F(J) + Q(J)*Q(J+1)                                      
11    CONTINUE                                                          
      F(J2) = ( FN - C(J2)*Q(J1) - A(J2)*Q(J2-1))/(C(J2)*S(J1) +        
     1     A(J2)*S(J2-1)  +B(J2))                                       
C                                                                       
C         BACKWARD ELIMINATION PASS                                     
C                                                                       
      DO 12 I=JA,J2                                                     
         J = JJ -I                                                      
         F(J) = F(J2)*S(J) + Q(J)                                       
12    CONTINUE                                                          
      RETURN                                                            
      END                                                               
                                                                        
      SUBROUTINE VERIFY(NUMITE,RSD)                                     
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
                                                                        
C VERIFY THE CORRECTNESS OF THE OUTPUT USING THE 100TH ITERATION.       
C THE VERIFICATION PARAMETER IS THE 2-NORM OF THE RESDIUAL.             
                                                                        
      COMMON/IO/IREAD,IWRIT,IVER                                        
      IF (NUMITE.EQ.100) THEN                                           
C                                                                       
         WRITE(IVER,1)                                                  
1        FORMAT(1X,
     * 'OUTPUT FOR PERFECT CLUB BENCHMARK: ARC2D INSTRUMENTED VERSION'/
     *          1X,'$Revision: 1.A.1 $ $Author: psharma $'/)                
         WRITE(IVER,*)                                                  
         WRITE(IVER,2)                                                  
2        FORMAT(1X,'VERIFICATION PARAMETERS:'/)                         
C                                                                       
         WRITE(IVER,3) NUMITE,RSD                                       
3        FORMAT(1X,'ITERATION = ',I5,3X,'L2 RESIDUAL = ',E12.6)         
C                                                                       
         IF (ABS((RSD-0.142979E-04)/.142979E-04).LT.0.1E-03) THEN       
            WRITE(IVER,4)                                               
4           FORMAT(//1X,'RESULTS FOR THIS RUN ARE: VALID')              
         ELSE                                                           
            WRITE(IVER,5)                                               
5           FORMAT(//1X,'RESULTS FOR THIS RUN ARE: INVALID')            
         ENDIF                                                          
                                                                        
      ENDIF                                                             
                                                                        
      RETURN                                                            
      END                                                               
