      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,   CP2TIM,
     1  FSMACH,     ALPHA,    GAMMA,       GAMI,        PI,
     1  DIS2X,      DIS2Y,    DIS4X,       DIS4Y,       PHIDT,    
     1  THETAD ,    RESID,    JACDT,       IPRINT,      NPCP,
     1  JTAIL1,     JTAIL2,   NUMITE ,     ISTART,      NSTEPS
      LOGICAL PERIDC  
      COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK
C
      DIMENSION Q(JDIM,KDIM,4),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)
C
C
      DTD   = DT * THETAD /(1. + PHIDT)
      HD    = 0.5*DTD
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
C
      RETURN
      END
