SUBROUTINE D9RPLY(SSS, NZ, IFLAG, P, QP, K, QK)
C VARIABLE-SHIFT H POLYNOMIAL ITERATION FOR A REAL
C ZERO.
C SSS   - STARTING ITERATE
C NZ    - NUMBER OF ZERO FOUND
C IFLAG - FLAG TO INDICATE A PAIR OF ZEROS NEAR REAL
C         AXIS.
C
      COMMON /P77PLY/ SR, SI, U,
     1 V, A, B, C, D, A1, A2, A3, A6, A7, E, F, G,
     2 H, SZR, SZI, LZR, LZI, ETA, ARE, MRE, N, NN
C
      INTEGER N, NN
      INTEGER NZ, IFLAG, I, J, NM1
C
      DOUBLE PRECISION ETA, ARE, MRE
      DOUBLE PRECISION MS, MP, OMP, EE
C
      DOUBLE PRECISION P(1), QP(1), K(1),
     1   QK(1), SR, SI, U, V, A, B, C, D,
     2   A1, A2, A3, A6, A7, E, F, G, H, SZR, SZI,
     3   LZR, LZI
      DOUBLE PRECISION PV, KV, T, S, SSS
      NM1 = N - 1
      NZ = 0
      S = SSS
      IFLAG = 0
      J = 0
C MAIN LOOP
   10 PV = P(1)
C EVALUATE P AT S
      QP(1) = PV
      DO 20 I=2,NN
        PV = PV*S + P(I)
        QP(I) = PV
   20 CONTINUE
      MP = DABS(PV)
C COMPUTE A RIGOROUS BOUND ON THE ERROR IN EVALUATING
C P
      MS = DABS(S)
      EE = (MRE/(ARE+MRE))*ABS(QP(1))
      DO 30 I=2,NN
        EE = EE*MS + ABS(QP(I))
   30 CONTINUE
C ITERATION HAS CONVERGED SUFFICIENTLY IF THE
C POLYNOMIAL VALUE IS LESS THAN 20 TIMES THIS BOUND
      IF (MP.GT.20.*((ARE+MRE)*EE-MRE*MP)) GO TO 40
      NZ = 1
      SZR = S
      SZI = 0.D0
      RETURN
   40 J = J + 1
C STOP ITERATION AFTER 10 STEPS
      IF (J.GT.10) RETURN
      IF (J.LT.2) GO TO 50
      IF (DABS(T).GT..001*DABS(S-T) .OR. MP.LE.OMP)
     1   GO TO 50
C A CLUSTER OF ZEROS NEAR THE REAL AXIS HAS BEEN
C ENCOUNTERED RETURN WITH IFLAG SET TO INITIATE A
C QUADRATIC ITERATION
      IFLAG = 1
      SSS = S
      RETURN
C RETURN IF THE POLYNOMIAL VALUE HAS INCREASED
C SIGNIFICANTLY
   50 OMP = MP
C COMPUTE T, THE NEXT POLYNOMIAL, AND THE NEW ITERATE
      KV = K(1)
      QK(1) = KV
      DO 60 I=2,N
        KV = KV*S + K(I)
        QK(I) = KV
   60 CONTINUE
      IF (DABS(KV).LE.DABS(K(N))*10.*ETA) GO TO 80
C USE THE SCALED FORM OF THE RECURRENCE IF THE VALUE
C OF K AT S IS NONZERO
      T = -PV/KV
      K(1) = QP(1)
      DO 70 I=2,N
        K(I) = T*QK(I-1) + QP(I)
   70 CONTINUE
      GO TO 100
C USE UNSCALED FORM
   80 K(1) = 0.0D0
      DO 90 I=2,N
        K(I) = QK(I-1)
   90 CONTINUE
  100 KV = K(1)
      DO 110 I=2,N
        KV = KV*S + K(I)
  110 CONTINUE
      T = 0.D0
      IF (DABS(KV).GT.DABS(K(N))*10.*ETA) T = -PV/KV
      S = S + T
      GO TO 10
      END

.