subroutine aptquar (c0, c1, c2, c3, tol, nroot, nrooti,
     &                    xroot, yroot, atype, itrun)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTQUAR
c
c     call aptquar (c0, c1, c2, c3, tol, nroot, nrooti,
c                   xroot, yroot, atype, itrun)
c
c     Version:  aptquar  Updated    1997 April 4 17:00.
c               aptquar  Originated 1997 March 14 18:40.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To find the real and any complex roots z = x + i*y of the
c               quartic equation:
c                    f(z) = c0 + c1 * z + c2 * z**2 + c3 * z**3 + z**4 = 0
c               and to return their number (nroot and nrooti), their values
c               (xroot and yroot), and their types (atype), and a trucation
c               indicator (itrun).
c               The method uses a tolerance limit (tol) to control and measure
c               the accuracy of the results.  Accuracy may be reduced if the
c               coefficients differ by many orders of magnitude. 
c               Any complex roots occur in conjugate pairs, so the quartic
c               equation may have only 0, 2 or 4 complex roots.
c
c               Note that for the four roots z1, z2, z3 and z4, the following
c               relationships may be used to check the accuracy of the roots:
c                    c0 =   z1*z2*z3*z4,
c                    c1 = -(z1*z2*z3 + z2*z3*z4 + z3*z4*z1 + z4*z1*z2),
c                    c2 =   z1*z2 + z2*z3 + z3*z4 + z4*z1 + z1*z3 + z2*z4,
c                    c3 = -(z1 + z2 + z3 + z4).
c
c               Also, if one real root x1 is known, the other three roots are
c               the roots of the cubic equation:
c                    q(z) = d0 + d1 * z + d2 * z**2 + z**3,
c               where d0 = c1 + c2 * x1 + c3 * x1**2 + x1**3
c                        = -c0 / x1  (x1 nonzero),
c                     d1 = c2 + c3 * x1 + x1**2
c                        = -(c0 + c1 * x1) / x1**2 (x1 nonzero),
c                     d2 = c3 + x1
c                        = -(c0 + c1 * x1 + c2 * x1**2) / x1**3 (x1 nonzero)
c               (for each coefficient, use the expression which minimizes the
c               numerical truncation error)
c               Either one or three of the roots of the cubic equation are real.
c               If three are real, one may be a double or triple root.
c
c               If a pair of roots z1 and z2 are known, the other two roots are
c               the roots of the quadratic equation:
c                    q(z) = d0 + d1 * z + z**2,
c               where d0 = c0 / (z1 * z2),
c               and   d1 = c3 + Z1 + z2.
c               There are either two real roots or a pair of conjugate complex
c               roots.  The two real roots may be the same.
c
c               If significant truncation error occurs (itrun = 2), try calling
c               aptquar with the coefficients c0' = c3 / c0, c1' = c2 / c0,
c               c2' = c1 / c0, c3' = 1 / c0, and use the reciprocals of the
c               resulting roots.  This may make very small roots more accurate.
c               It may not be possible to get accurate values of complex roots
c               for which the real and imaginary parts differ by many orders of
c               magnitude.
c
c     Method:   The method in the "Handbook of tables for Mathematics", 3rd Ed.,
c               1967, page 133 is used.
c
c               Any real root of maximum magnitude is then factored out of the
c               quartic equation, and the resulting cubic equation is solved.
c
c               Any zero roots are recalculated to eliminate as much of the
c               truncation error as possible.
c
c               Finally, all roots are reordered, starting with real roots in
c               ascending numerical order, followed by any pair of conjugate
c               complex roots.
c
c               The accuracy criterion tol is used throughout to replace
c               calculated quantities with zero when they have values less than
c               the estimated error in their calculation, based on tol.
c
c     Input:    c0, c1, c2, c3, tol.
c
c     Output:   nroot, nrooti, xroot, yroot, itrun.
c
c     Glossary:
c
c     atype     Output   Types of root:
c                          "single"   Single unique real root.
c                          "dbl max"  Double unique real root at maximum.
c                          "dbl min"  Double unique real root at minimum.
c                          "triple"   Triple unique real root.
c                          "quadrpl"  Quadruple unique real root.
c                          "complex"  One of two or four unique complex roots.
c                          "dbl com"  One of two double unique complex roots.
c                          "unknown"  Unevaluated.
c                          The array size must be at least 4.
c
c     c0        Input    Coefficient of 1    in the quartic equation.
c
c     c1        Input    Coefficient of z    in the quartic equation.
c
c     c2        Input    Coefficient of z**2 in the quartic equation.
c
c     c3        Input    Coefficient of z**3 in the quartic equation.
c
c     itrun     Output   Indicates the degree of truncation error that occurred
c                          in finding roots of the quartic function.
c                          0 if none (within machine accuracy).
c                          1 if the truncation error is less than tol.
c                          2 if the truncation error exceeds tol.  The residual
c                            value of the function at the roots may be
c                            excessive (try increasing tol).  The coefficients
c                            recalculated from the roots may differ excessively
c                            from the original coefficients.  This may mean the
c                            magnitudes of the coefficients differ too much to
c                            get an accurate result.
c
c     nroot     Output   Number of unique real roots. Possible values:
c                          0 (nrooti = 0):
c                            The method failed.  Try changing tol.
c                          0 (nrooti = 2)
c                            "dbl com", "dbl com"
c                          0 (nrooti = 4):
c                            "complex", "complex", "complex", "complex".
c                          1 (nrooti = 0):
c                            "quadrpl" (at the minimum point).
c                          1 (nrooti = 2):
c                            "dbl min", "complex", "complex".
c                          2 (nrooti = 0):
c                            "single", "triple" (at an inversion point), or
c                            "dbl min", "dbl min".
c                          2 (nrooti = 2):
c                            "single", "single", "complex", "complex".
c                          3 (nrooti = 0):
c                            "dbl max", "single", "single", or
c                            "dbl min", "single", "single".
c                          4 (nrooti = 0):
c                            "single", "single", "single", "single".
c
c     nrooti    Output   The number of complex roots, which must be either 0, 2
c                          or 4, as they only occur in conjugate complex pairs.
c                          See nroot.
c
c     tol       Input    Numerical tolerance limit and convergence control.
c                          If zero, itrun will almost certainly be 2.
c                          On 64-bit computers, or 32-bit computers with
c                          double precision, recommend 1.e-5 to 1.e-11.
c
c     xroot     Output   The real parts of the unique roots of the quartic
c                          function.  The array size must be at least 4.
c                          Real roots are first, in increasing numerical order,
c                          followed by any complex roots, in increasing
c                          numerical order of the imaginary part.
c                          May truncate to zero if yroot is many orders of
c                          magnitude greater than xroot.
c
c     yroot     Output   The imaginary parts of the complex roots of the quartic
c                          function.  The array size must be at least 4.
c                          Always zero for the first nroot roots.
c                          Nonzero for roots (nroot + 1) through
c                          (nroot + nrooti).
c                          May truncate to zero if xroot is many orders of
c                          magnitude greater than yroot.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

c.... Roots of the quartic equation z**4 + c3*z**3 + c2*z**2 + c1*z + c0 = 0.
      dimension xroot(4)              ! The real part of a root.
      dimension yroot(4)              ! The imaginary part of a root.
c.... Quadric root types:
c....   "single", "dbl max", "dbl min", "triple", "quadrpl", "complex".
      dimension atype(4)              ! Root type for quartic.
      character*8 atype               ! Root type for quartic.

c.... Cubic root types:
c....   "single", "dbl max", "dbl min", "triple", "complex".
      dimension atcub(3)              ! Root type for cubic.
      character*8 atcub               ! Root type for cubic.
      dimension xcubr(3)              ! The real part of a root of a cubic.
      dimension xcubi(3)              ! The imaginary part of a root of a cubic.


cbugc***DEBUG begins.
cbug 9901 format (/ 'aptquar finding roots of a quartic.  tol=',1pe10.2 /
cbug     &  '  Equation is c0 + c1 * z + c2 * z**2 + c3 * z**2 + z**4 = 0' /
cbug     &  '  c0,c1,c2=   ',1p3e22.14 /
cbug     &  '  c3=         ',1pe22.14)
cbug      write ( 3, 9901) tol, c0, c1, c2, c3
cbugc***DEBUG ends.

c.... Initialize.

      nroot    = 0
      nrooti   = 0
      itrun    = 0
      xroot(1) = -1.e+99
      xroot(2) = -1.e+99
      xroot(3) = -1.e+99
      xroot(4) = -1.e+99
      yroot(1) = 0.0
      yroot(2) = 0.0
      yroot(3) = 0.0
      yroot(4) = 0.0
      atype(1) = 'unknown'
      atype(2) = 'unknown'
      atype(3) = 'unknown'
      atype(4) = 'unknown'

c.... See if one of the roots is zero.

      if (c0 .eq. 0.0) then           ! At least one root is zero.
        xroot(1) = 0.0
        yroot(1) = 0.0

        call aptcubs (c1, c2, c3, tol, nroots, xroot(2), yroot(2),
     &                atype(2), itrun)

        if (nroots .eq. 0) then
          kbad  = 1
          itrun = 2
          go to 395
        elseif (atype(2) .eq. 'triple') then
          xroot(3) = xroot(2)
          yroot(3) = yroot(2)
          xroot(4) = xroot(2)
          yroot(4) = yroot(2)
        elseif ((atype(2) .eq. 'dbl min')  .or.
     &          (atype(2) .eq. 'dbl max')) then
          xroot(4) = xroot(2)
          yroot(4) = yroot(2)
        elseif ((atype(3) .eq. 'dbl min')  .or.
     &          (atype(3) .eq. 'dbl max')) then
          xroot(4) = xroot(3)
          yroot(4) = yroot(3)
        endif

        go to 200
      endif                           ! At least one root is zero.

c.... Find the coefficients of the resolvent cubic equation.

c....   f(u) = a0 + a1 * u + a2 * u**2 + u**3
c....        = (-c3**2 * c0 + 4 * c2 * c0 - c**2) +
c....          (c3 * c1 - 4 * c0) * u - c2 * u**2 + u**3 = 0

      a0    = -c3**2 * c0 + 4.0 * c2 * c0 - c1**2
      a0err = tol * (c3**2 * abs (c0) + 4.0 * abs (c2 * c0) + c1**2)
      if (abs (a0) .le. a0err) then
        if ((itrun .eq. 0) .and. (a0 .ne. 0.0)) itrun = 1
        a0 = 0.0
      endif
      a1    = c3 * c1 - 4.0 * c0
      a1err = tol * (abs (c3 * c1) + 4.0 * abs (c0))
      if (abs (a1) .le. a1err) then
        if ((itrun .eq. 0) .and. (a1 .ne. 0.0)) itrun = 1
        a1 = 0.0
      endif
      a2    = -c2

c.... Find the most positive real root of the resolvent cubic equation.

      call aptcubs (a0, a1, a2, tol,
     &              nroots, xcubr, xcubi, atcub, itrun)

      u1 = xcubr(nroots)
cbugc***DEBUG begins.
cbug 9701 format ('aptquar DEBUG.  Resolvent root u1 = ',1pe22.14)
cbug      write ( 3, 9701) u1
cbugc***DEBUG ends.

c.... Find R = sqrt (rsq) = rreal + i * rimag, where i = sqrt (-1),
c....   and rsq = c3**2 / 4 - c2 + u1.
c....   1 / R = rreal / abs (rsq) + i * rimag / abs (rsq).

      rsq = 0.25 * c3**2  - c2 + u1
      rsqerr = tol * (0.25 * c3**2 + abs (c2) + abs (u1))
      if (abs (rsq) .le. rsqerr) then
        if ((itrun .eq. 0) .and. (rsq .ne. 0.0)) itrun = 1
        rsq = 0.0
      endif
      if (rsq .lt. 0.0) then
        rreal = 0.0
        rimag = sqrt (-rsq)
      else
        rreal = sqrt (rsq)
        rimag = 0.0
      endif
cbugc***DEBUG begins.
cbug 9702 format ('aptquar DEBUG.  R   = ',1pe22.14,'  + i *',1pe22.14)
cbug      write ( 3, 9702) rreal, rimag
cbugc***DEBUG ends.

      rmag = sqrt (abs (rsq))

c.... Find dcap = sqrt (dsq) = sqrt (ff + gg) and
c....      ecap = sqrt (esq) = sqrt (ff - gg),
c....   where ff = 0.75 * c3**2 - 2.0 * c2 - rsq, and
c....   and   gg = 2 * sqrt (u1**2 - 4 * c0)              (R = 0),
c....   or    gg = (c3 * c2 - 2 * c1 - 0.25 * c3**3) / R  (R nonzero).
c....         gg = ggreal + i * ggimag

      ff    = 0.75 * c3**2 - 2.0 * c2 - rsq
      fferr = tol * (0.75 * c3**2 + 2.0 * abs (c2) + abs (rsq))
      if (abs (ff) .le. fferr) then
        if ((itrun .eq. 0) .and. (ff .ne. 0.0)) itrun = 1
        ff = 0.0
      endif
cbugc***DEBUG begins.
cbug 9703 format ('aptquar DEBUG.  ff  = ',1pe22.14)
cbug      write ( 3, 9703) ff
cbugc***DEBUG ends.

      if (rmag .eq. 0.0) then
        ff    = 0.75 * c3**2 - 2.0 * c2
        fferr = tol * (0.75 * c3**2 + 2.0 * abs (c2))
        if (abs (ff) .le. fferr) then
          if ((itrun .eq. 0) .and. (ff .ne. 0.0)) itrun = 1
          ff = 0.0
        endif
        ggarg = u1**2 - 4.0 * c0
        ggargerr = tol * (u1**2 + 4.0 * abs (c0))
        if (abs (ggarg) .le. ggargerr) then
          if ((itrun .eq. 0) .and. (ggarg .ne. 0.0)) itrun = 1
           ggarg = 0.0
        endif
        if (ggarg .lt. 0.0) then
          ggreal = 0.0
          ggimag = 2.0 * sqrt (-ggarg)
        else
          ggreal = 2.0 * sqrt (ggarg)
          ggimag = 0.0
        endif
      else
        ff    = 0.5 * c3**2 - c2 - u1
        fferr = tol * (0.5 * c3**2 + abs (c2) + abs (u1))
        if (abs (ff) .le. fferr) then
          if ((itrun .eq. 0) .and. (ff .ne. 0.0)) itrun = 1
          ff = 0.0
        endif
        fact    = c3 * c2 - 2.0 * c1 - 0.25 * c3**3
        facterr = tol * (abs (c3 * c2) + 2.0 * abs (c1) +
     &                   0.25 * abs (c3**3))
        if (abs (fact) .le. facterr) then
          if ((itrun .eq. 0) .and. (fact .ne. 0.0)) itrun = 1
          fact = 0.0
        endif
        ggreal =  fact * rreal / abs (rsq)
        ggimag = -fact * rimag / abs (rsq)
      endif
cbugc***DEBUG begins.
cbug 9704 format ('aptquar DEBUG.  gg  = ',1pe22.14,'  + i *',1pe22.14)
cbug      write ( 3, 9704) ggreal, ggimag
cbugc***DEBUG ends.

      dsqreal = ff + ggreal
      dsqerr = tol * (abs (ff) + abs (ggreal))
      if (abs (dsqreal) .le. dsqerr) then
        if ((itrun .eq. 0) .and. (dsqreal .ne. 0.0)) itrun = 1
        dsqreal = 0.0
      endif
      dsqimag = ggimag
      esqreal = ff - ggreal
      if (abs (esqreal) .le. dsqerr) then
        if ((itrun .eq. 0) .and. (esqreal .ne. 0.0)) itrun = 1
        esqreal = 0.0
      endif
      esqimag = -ggimag

cbugc***DEBUG begins.
cbug 9705 format ('aptquar DEBUG.  Dsq = ',1pe22.14,'  + i *',1pe22.14)
cbug 9706 format ('aptquar DEBUG.  Esq = ',1pe22.14,'  + i *',1pe22.14)
cbug      write ( 3, 9705) dsqreal, dsqimag
cbug      write ( 3, 9706) esqreal, esqimag
cbugc***DEBUG ends.

      if (dsqimag .eq. 0.0) then      ! dsq is real.
        if (dsqreal .gt. 0.0) then
          dcapreal = sqrt (dsqreal)
          dcapimag = 0.0
        elseif (dsqreal .eq. 0.0) then
          dcapreal = 0.0
          dcapimag = 0.0
        else
          dcapreal = 0.0
          dcapimag = sqrt (-dsqreal)
        endif
      else                            ! dsq is complex.
        dmag = sqrt (dsqreal**2 + dsqimag**2)
        if (abs (dsqreal) .gt. abs (dsqimag)) then
          if (dsqreal .gt. 0.0) then
            cosfun  = 1.0 + dsqreal / dmag
            sinth   = dsqimag / dmag
            cosdang = sqrt (0.5 * cosfun)
            sindang = sinth / sqrt (2.0 * cosfun)
          else
            cosfun  = 1.0 - dsqreal / dmag
            sinth   = dsqimag / dmag
            cosdang = sinth / sqrt (2.0 * cosfun)
            sindang = sqrt (0.5 * cosfun)
          endif
        else
          dang = 0.5 * atan2 (dsqimag, dsqreal)
          cosdang = cos (dang)
          sindang = sin (dang)
        endif
        dmagrt   = sqrt (dmag)
        dcapreal = dmagrt * cosdang
        dcapimag = dmagrt * sindang
      endif

      if (esqimag .eq. 0.0) then      ! esq is real.
        if (esqreal .gt. 0.0) then
          ecapreal = sqrt (esqreal)
          ecapimag = 0.0
        elseif (esqreal .eq. 0.0) then
          ecapreal = 0.0
          ecapimag = 0.0
        else
          ecapreal = 0.0
          ecapimag = sqrt (-esqreal)
        endif
      else                            ! esq is complex.
        emag = sqrt (esqreal**2 + esqimag**2)
        if (abs (esqreal) .gt. abs (esqimag)) then
          if (esqreal .gt. 0.0) then
            cosfun  = 1.0 + esqreal / emag
            sinth   = esqimag / emag
            coseang = sqrt (0.5 * cosfun)
            sineang = sinth / sqrt (2.0 * cosfun)
          else
            cosfun  = 1.0 - esqreal / emag
            sinth   = esqimag / emag
            coseang = sinth / sqrt (2.0 * cosfun)
            sineang = sqrt (0.5 * cosfun)
          endif
        else
          eang = 0.5 * atan2 (esqimag, esqreal)
          coseang = cos (eang)
          sineang = sin (eang)
        endif
        emagrt   = sqrt (emag)
        ecapreal = emagrt * coseang
        ecapimag = emagrt * sineang
      endif
cbugc***DEBUG begins.
cbug 9707 format ('aptquar DEBUG.  dcap= ',1pe22.14,'  + i *',1pe22.14)
cbug 9708 format ('aptquar DEBUG.  ecap= ',1pe22.14,'  + i *',1pe22.14)
cbug      write ( 3, 9707) dcapreal, dcapimag
cbug      write ( 3, 9708) ecapreal, ecapimag
cbugc***DEBUG ends.

c.... Find the four roots of the quartic equation.

      xroot(1) = -0.25 * c3 + 0.5 * (-rreal - ecapreal)
      xrerr    = tol * (0.25 * abs (c3) +
     &                  0.5 * (abs (rreal) + abs (ecapreal)))
      if (abs (xroot(1)) .le. xrerr) then
        if ((itrun .eq. 0) .and. (xroot(1) .ne. 0.0)) itrun = 1
        xroot(1) = 0.0
      endif

      yroot(1) =              0.5 * (-rimag - ecapimag)
      yrerr    = tol * (0.5 * (abs (rimag) + abs (ecapimag)))
      if (abs (yroot(1)) .le. yrerr) then
        if ((itrun .eq. 0) .and. (yroot(1) .ne. 0.0)) itrun = 1
        yroot(1) = 0.0
      endif

      xroot(2) = -0.25 * c3 + 0.5 * (-rreal + ecapreal)
      if (abs (xroot(2)) .le. xrerr) then
        if ((itrun .eq. 0) .and. (xroot(2) .ne. 0.0)) itrun = 1
        xroot(2) = 0.0
      endif

      yroot(2) =              0.5 * (-rimag + ecapimag)
      if (abs (yroot(2)) .le. yrerr) then
        if ((itrun .eq. 0) .and. (yroot(2) .ne. 0.0)) itrun = 1
        yroot(2) = 0.0
      endif

      xroot(3) = -0.25 * c3 + 0.5 * ( rreal - dcapreal)
      xrerr    = tol * (0.25 * abs (c3) +
     &                  0.5 * (abs (rreal) + abs (dcapreal)))
      if (abs (xroot(3)) .le. xrerr) then
        if ((itrun .eq. 0) .and. (xroot(3) .ne. 0.0)) itrun = 1
        xroot(3) = 0.0
      endif

      yroot(3) =              0.5 * ( rimag - dcapimag)
      yrerr    = tol * (0.5 * (abs (rimag) + abs (dcapimag)))
      if (abs (yroot(3)) .le. yrerr) then
        if ((itrun .eq. 0) .and. (yroot(3) .ne. 0.0)) itrun = 1
        yroot(3) = 0.0
      endif

      xroot(4) = -0.25 * c3 + 0.5 * ( rreal + dcapreal)
      if (abs (xroot(4)) .le. xrerr) then
        if ((itrun .eq. 0) .and. (xroot(4) .ne. 0.0)) itrun = 1
        xroot(4) = 0.0
      endif

      yroot(4) =              0.5 * ( rimag + dcapimag)
      if (abs (yroot(4)) .le. yrerr) then
        if ((itrun .eq. 0) .and. (yroot(4) .ne. 0.0)) itrun = 1
        yroot(4) = 0.0
      endif
cbugc***DEBUG begins.
cbug 9810 format ('  Root ',i1,'    r,i=',1p2e22.14)
cbug      write ( 3, 9810) (n, xroot(n), yroot(n), n = 1, 4)
cbugc***DEBUG ends.

c.... Find the root of maximum magnitude.

      nmax   = 0
      rsqmax = 0.0
      do 180 n = 1, 4
        rsqn = xroot(n)**2 + yroot(n)**2
        if (rsqn .gt. rsqmax) then
          nmax   = n
          rsqmax = rsqn
        endif
  180 continue
cbugc***DEBUG begins.
cbug 9812 format ('  Root ',i1,' (max) =',1p2e22.14,'  (real, imag)' )
cbug
cbug        write ( 3, 9812) nmax, xroot(nmax), yroot(nmax)
cbugc***DEBUG ends.

c.... Factor the maximum root(s) out of the quartic equation.

      if (yroot(nmax) .eq. 0.0) then  ! Maximum root is real.

c....   Factor out one real root.

        x1  = xroot(nmax)
        x1a = abs (x1)
        c0a = abs (c0)
        c1a = abs (c1)
        c2a = abs (c2)
        c3a = abs (c3)

        d0err  = tol * (c1a + x1a * (c2a + x1a * (c3a + x1a)))
        d0err2 = tol * (c0a / x1a)
        if (d0err2 .le. d0err) then
          d0 = -c0 / x1
        else
          d0 = c1 + x1 * (c2 + x1 * (c3 + x1))
        endif

        d1err  = tol * (c2a + x1a * (c3a + x1a))
        d1err2 = tol * (c0a + c1a * x1a) / x1**2
        if (d1err2 .le. d1err) then
          d1     = -(c0 + c1 * x1) / x1**2
        else
          d1     = c2 + x1 * (c3 + x1)
        endif

        d2err  = tol * abs (c3a + x1a)
        d2err2 = tol * (c0a + x1a * (c1 * x1a + c2 * x1) / x1a**3)
        if (d2err2 .le. d2err) then
          d2 = -(c0 + x1 * (c1 + c2 * x1)) / x1**3
        else
          d2 = c3 + x1
        endif
cbugc***DEBUG begins.
cbug 9808 format ('  d0,d1,d2= ',1p3e22.14)
cbug      write ( 3, 9808) d0, d1, d2
cbugc***DEBUG ends.

        xroot(1) = x1
        yroot(1) = 0.0

        call aptcubs (d0, d1, d2, tol, nroots, xroot(2), yroot(2),
     &                atype(2), itrun)

        if (nroots .eq. 0) then
          kbad  = 1
          itrun = 2
          go to 395
        elseif (atype(2) .eq. 'triple') then
          xroot(3) = xroot(2)
          yroot(3) = yroot(2)
          xroot(4) = xroot(2)
          yroot(4) = yroot(2)
        elseif ((atype(2) .eq. 'dbl min')  .or.
     &          (atype(2) .eq. 'dbl max')) then
          xroot(4) = xroot(2)
          yroot(4) = yroot(2)
        elseif ((atype(3) .eq. 'dbl min')  .or.
     &          (atype(3) .eq. 'dbl max')) then
          xroot(4) = xroot(3)
          yroot(4) = yroot(3)
        endif

cbugc***DEBUG begins.
cbug 9809 format ('aptquar DEBUG.  From aptcubs:  nroots=',i2)
cbug      write ( 3, 9809) nroots
cbug 9811 format ('  Root ',i1,' (',a7,')',1p2e22.14,'  (real, imag)' )
cbug      write ( 3, 9811) (n, atype(n), xroot(n), yroot(n), n = 1, 4)
cbugc***DEBUG ends.

cbugc***DEBUG begins.
cbug      write ( 3, 9810) (n, xroot(n), yroot(n), n = 1, 4)
cbugc***DEBUG ends.
      else                            ! Maximum roots are complex conjugates.

c....   Factor out two conjugate complex roots.

        x1    = xroot(nmax)
        y1    = yroot(nmax)

        d0    = c0 / (x1**2 + y1**2)
        d1    = c3 + 2.0 * x1
        d1err = tol * (abs (c3) + 2.0 * abs (x1))
        if (abs (d1) .le. d1err) d1 = 0.0
        d2    = 1.0
cbugc***DEBUG begins.
cbug 9910 format ('  d0,d1       = ',1p2e22.14)
cbug      write ( 3, 9910) d0, d1
cbugc***DEBUG ends.

        call aptqrts (0, d2, d1, d0, qq, tol,
     &                nrootx, xnew1, xnew2, itrun)

        xroot(1) = x1
        yroot(1) = -abs (y1)
        xroot(2) = x1
        yroot(2) =  abs (y1)

        if (nrootx .eq. 1) then
          xroot(3) = xnew1
          yroot(3) = 0.0
          xroot(4) = xnew1
          yroot(4) = 0.0
        elseif (nrootx .eq. 2) then
          xroot(3) = xnew1
          yroot(3) = 0.0
          xroot(4) = xnew2
          yroot(4) = 0.0
        else
          xroot(3) = xnew1
          yroot(3) = -abs (xnew2)
          xroot(4) = xnew1
          yroot(4) =  abs (xnew2)
        endif

cbugc***DEBUG begins.
cbug      write ( 3, 9810) (n, xroot(n), yroot(n), n = 1, 4)
cbugc***DEBUG ends.

      endif                           ! Factored out largest real root.

c.... Reorder the roots in ascending order of absolute magnitude.

      nzero = 0

      do 192 n = 1, 4
        xsqn = xroot(n)**2 + yroot(n)**2
        if (xsqn .eq. 0.0) nzero = nzero + 1
  192 continue

      do 195 n = 1, 3
        do 190 nn = n + 1, 4
          xsqn =  xroot(n)**2 + yroot(n)**2
          xsqnn = xroot(nn)**2 + yroot(nn)**2
          if (xsqnn .lt. xsqn) then
            xrootx    = xroot(n)
            yrootx    = yroot(n)
            xroot(n)  = xroot(nn)
            yroot(n)  = yroot(nn)
            xroot(nn) = xrootx
            yroot(nn) = yrootx
          endif
  190   continue
  195 continue
cbugc***DEBUG begins.
cbug 9815 format ('aptquar DEBUG.  nzero=',i1)
cbug      write ( 3, 9815) nzero
cbug      write ( 3, 9810) (n, xroot(n), yroot(n), n = 1, 4)
cbugc***DEBUG ends.

c.... See if there is a double root at zero.  If so, recalculate.

      if (nzero .eq. 2) then
        d0    = c0 / (xroot(3) * xroot(4) - yroot(3) * yroot(4))
        d1    = c3 + xroot(3) + xroot(4)
        d1err = tol * (abs (c3) + abs (xroot(3)) + abs (xroot(4)))
        if (abs (d1) .le. d1err) then
          if ((d1 .ne. 0.0) .and. (itrun .eq. 0)) itrun = 1
          d1 = 0.0
        endif
        d2    = 1.0

        call aptqrts (0, d2, d1, d0, qq, tol,
     &                nrootx, xnew1, xnew2, itrun)

        if (nrootx .eq. 2) then
          xroot(1) = xnew1
          yroot(1) = 0.0
          xroot(2) = xnew2
          yroot(2) = 0.0
        else
          xroot(1) =  xnew1
          yroot(1) = -xnew2
          xroot(2) =  xnew1
          yroot(2) =  xnew2
        endif
cbugc***DEBUG begins.
cbug 9818 format ('aptquar DEBUG.  Replaced double root at zero.')
cbug      write ( 3, 9818)
cbug      write ( 3, 9810) (n, xroot(n), yroot(n), n = 1, 4)
cbugc***DEBUG ends.

      endif

c.... Reorder the roots in ascending order of real, then complex.  Find types.

  200 do 220 n = 1, 3                 ! Loop over roots
        do 210 nn = n + 1, 4          ! Loop over roots.
c....     Put real roots in ascending order.
          if ((yroot(n)  .eq. 0.0)  .and.
     &        (yroot(nn) .eq. 0.0)) then
            atype(n)  = 'single'
            atype(nn) = 'single'
            if (xroot(n) .gt. xroot(nn)) then
              xrootx    = xroot(n)
              xroot(n)  = xroot(nn)
              xroot(nn) = xrootx
            endif
c....     Put complex roots in ascending order of real, then imaginary parts.
          elseif ((yroot(n)  .ne. 0.0)        .and.
     &            (yroot(nn) .ne. 0.0))       then
            atype(n)  = 'complex'
            atype(nn) = 'complex'
            if (xroot(n)  .gt. xroot(nn))  then
              xrootx    = xroot(n)
              xroot(n)  = xroot(nn)
              xroot(nn) = xrootx
              yrootx    = yroot(n)
              yroot(n)  = yroot(nn)
              yroot(nn) = yrootx
            elseif (xroot(n)  .eq. xroot(nn))  then
              if (abs (yroot(n)) .gt. abs (yroot(nn))) then
                xrootx    = xroot(n)
                xroot(n)  = xroot(nn)
                xroot(nn) = xrootx
                yrootx    = yroot(n)
                yroot(n)  = yroot(nn)
                yroot(nn) = yrootx
              elseif ((abs (yroot(n)) .eq. abs (yroot(nn))) .and.
     &                (yroot(nn) .lt. 0.0)) then
                xrootx    = xroot(n)
                xroot(n)  = xroot(nn)
                xroot(nn) = xrootx
                yrootx    = yroot(n)
                yroot(n)  = yroot(nn)
                yroot(nn) = yrootx
              endif
            endif
c....     Put real roots before complex roots.
          elseif ((yroot(n)  .ne. 0.0)  .and.
     &            (yroot(nn) .eq. 0.0)) then
            atype(n)  = 'single'
            atype(nn) = 'complex'
            xrootx    = xroot(n)
            xroot(n)  = xroot(nn)
            xroot(nn) = xrootx
            yrootx    = yroot(n)
            yroot(n)  = yroot(nn)
            yroot(nn) = yrootx
          elseif ((yroot(n)  .eq. 0.0)  .and.
     &            (yroot(nn) .ne. 0.0)) then
            atype(n)  = 'single'
            atype(nn) = 'complex'
          endif

  210   continue                      ! Loop over roots.
  220 continue                        ! Loop over roots

c.... Reconstruct the coefficients.

      kbad = 0
      c0p =   xroot(1) *  xroot(2) *  xroot(3) * xroot(4)   -
     &       (xroot(1) *  xroot(2) *  yroot(3) * yroot(4)   +
     &        xroot(2) *  xroot(3) *  yroot(4) * yroot(1)   +
     &        xroot(3) *  xroot(4) *  yroot(1) * yroot(2)   +
     &        xroot(4) *  xroot(1) *  yroot(2) * yroot(3)   +
     &        xroot(1) *  xroot(3) *  yroot(2) * yroot(4)   +
     &        xroot(2) *  xroot(4) *  yroot(1) * yroot(3))  +
     &        yroot(1) *  yroot(2) *  yroot(3) * yroot(4)

      c1p = -(xroot(1) *  xroot(2) * (xroot(3) + xroot(4))  +
     &        xroot(3) *  xroot(4) * (xroot(1) + xroot(2))) +
     &        xroot(1) * (yroot(2) * (yroot(3) + yroot(4))  +
     &                    yroot(3) *  yroot(4)           )  +
     &        xroot(2) * (yroot(3) * (yroot(4) + yroot(1))  +
     &                    yroot(4) *  yroot(1)           )  +
     &        xroot(3) * (yroot(4) * (yroot(1) + yroot(2))  +
     &                    yroot(1) *  yroot(2)           )  +
     &        xroot(4) * (yroot(1) * (yroot(2) + yroot(3))  +
     &                    yroot(2) *  yroot(3)              )

      c2p =   xroot(1) *  xroot(2) +  xroot(2) * xroot(3)   +
     &        xroot(3) *  xroot(4) +  xroot(4) * xroot(1)   +
     &        xroot(1) *  xroot(3) +  xroot(2) * xroot(4)   -
     &       (yroot(1) *  yroot(2) +  yroot(2) * yroot(3)   +
     &        yroot(3) *  yroot(4) +  yroot(4) * yroot(1)   +
     &        yroot(1) *  yroot(3) +  yroot(2) * yroot(4)   )

      c3p = -(xroot(1) +  xroot(2)  +  xroot(3)  + xroot(4))

      if (abs (c0 - c0p) .gt. tol * (abs (c0) + abs (c0p))) kbad = 1
      if (abs (c1 - c1p) .gt. tol * (abs (c1) + abs (c1p))) kbad = 1
      if (abs (c2 - c2p) .gt. tol * (abs (c2) + abs (c2p))) kbad = 1
      if (abs (c3 - c3p) .gt. tol * (abs (c3) + abs (c3p))) kbad = 1

cbugc***DEBUG begins.
cbug 9930 format ('aptquar DEBUG.  itrun =',i2, /
cbug     &        '  c0,c1(input) =',1p2e22.14 /
cbug     &        '  c0,c1(recalc)=',1p2e22.14 //
cbug     &        '  c2,c3(input) =',1p2e22.14 /
cbug     &        '  c2,c3(recalc)=',1p2e22.14 )
cbug      write ( 3, 9930) itrun, c0, c1, c0p, c1p, c2, c3, c2p, c3p
cbugC***DEBUG ends.

c.... Find and classify unique roots.

      nroots = 1
      do 230 n = 2, 4
        xdiff = abs (xroot(n) - xroot(n-1))
        ydiff = abs (yroot(n) - yroot(n-1))
        xderr = tol * (abs (xroot(n)) + abs (xroot(n-1)))
        yderr = tol * (abs (yroot(n)) + abs (yroot(n-1)))
        if ((xdiff .le. xderr)  .and.
     &      (ydiff .le. yderr)) then
          if (atype(nroots) .eq. 'single') then
            curve = c2 + 3.0 * c3 * xroot(n) + 6.0 * xroot(n)**2
            if (curve .lt. 0.0) then
              atype(nroots) = 'dbl max'
            else
              atype(nroots) = 'dbl min'
            endif
          elseif ((atype(nroots) .eq. 'dbl min')  .or.
     &            (atype(nroots) .eq. 'dbl max')) then
            atype(nroots) = 'triple'
          elseif (atype(nroots) .eq. 'triple') then
            atype(nroots) = 'quadrpl'
          elseif (atype(nroots) .eq. 'complex') then
            atype(nroots) = 'dbl com'
          endif
        else                          ! This root is unique.
          nroots = nroots + 1
          atype(nroots) = atype(n)
          xroot(nroots) = xroot(n)
          yroot(nroots) = yroot(n)
        endif
  230 continue

c.... Count real and imaginary roots.
      
      do 250 n = 1, nroots
        if (yroot(n) .eq. 0.0) then
          nroot  = nroot  + 1
        else
          nrooti = nrooti + 1
        endif
  250 continue

      if (nroots .lt. 4) then
        do 260 n = 4, nroots + 1, -1
          atype(n) = 'unknown'
  260   continue
      endif

c.... Store any duplicate roots in root locations nroots + 1, 4.

      next = nroots
      do 270 n = 1, nroots
        if ((atype(n) .eq. 'dbl min')  .or.
     &      (atype(n) .eq. 'dbl max')  .or.
     &      (atype(n) .eq. 'dbl com')) then
          next = next + 1
          xroot(next) = xroot(n)
          yroot(next) = yroot(n)
        elseif (atype(n) .eq. 'triple') then
          next = next + 1
          xroot(next) = xroot(n)
          yroot(next) = yroot(n)
          next = next + 1
          xroot(next) = xroot(n)
          yroot(next) = yroot(n)
        elseif (atype(n) .eq. 'quadrpl') then
          next = next + 1
          xroot(next) = xroot(n)
          yroot(next) = yroot(n)
          next = next + 1
          xroot(next) = xroot(n)
          yroot(next) = yroot(n)
          next = next + 1
          xroot(next) = xroot(n)
          yroot(next) = yroot(n)
        endif
  270 continue
cbugc***DEBUG begins.
cbug 9920 format (/ 'aptquar results.  nroot=',i2,'  nrooti=',i3)
cbug      write ( 3, 9920) nroot, nrooti
cbugc***DEBUG ends.

      do 390 n = 1, nroots

      frres = 0.0

        frres = xroot(n)**4 + c3 * xroot(n)**3 + c2 * xroot(n)**2 +
     &          c1 * xroot(n) + c0 +
     &          yroot(n)**2 * (yroot(n)**2 - 6.0 * xroot(n)**2 -
     &          3.0 * c3 * xroot(n) - c2)

        frerr = tol * (xroot(n)**4 + abs (c3 * xroot(n)**3) +
     &          abs (c2) * xroot(n)**2 + abs (c1 * xroot(n)) +
     &          abs (c0) +
     &          yroot(n)**2 * (yroot(n)**2 + 6.0 * xroot(n)**2 +
     &          3.0 * abs (c3 * xroot(n)) + abs (c2)))

        fires = yroot(n) * (4.0 * xroot(n) *
     &          (xroot(n)**2 - yroot(n)**2) +
     &          c3 * (3.0 * xroot(n)**2 - yroot(n)**2) +
     &          2.0 * c2 * xroot(n) + c1)

        fierr = tol * (abs (yroot(n)) * (4.0 * abs (xroot(n)) *
     &          (xroot(n)**2 + yroot(n)**2) +
     &          abs (c3) * (3.0 * xroot(n)**2 + yroot(n)**2) +
     &          2.0 * abs (c2 * xroot(n)) + abs (c1)))

        if (abs (frres) .gt. frerr) kbad = 1
        if (abs (fires) .gt. fierr) kbad = 1

cbugc***DEBUG begins.
cbug        if (kbad .eq. 1) then         ! Try to improve with Newton's method.
cbug
cbugc....     Find the partial derivatives of the real and imaginary parts.
cbug
cbug          xslop = c1 + 2.0 * c2 * xroot(n) + 3.0 * c3 * xroot(n)**2 +
cbug     &            4.0 * xroot(n)**3 -
cbug     &            3.0 * yroot(n)**2 * (c3 + 4.0 * xroot(n))
cbug
cbug          xserr = tol * (abs (c1) + 2.0 * abs (c2 * xroot(n)) +
cbug     &            3.0 * abs (c3) * xroot(n)**2 +
cbug     &            4.0 * xroot(n)**3 +
cbug     &            3.0 * yroot(n)**2 * (abs (c3) + 4.0 * abs (xroot(n))))
cbug
cbug          yslop = 4.0 * yroot(n)**3 - 2.0 * yroot(n) *
cbug     &            (c2 + 3.0 * c3 * xroot(n) + 6.0 * xroot(n)**2)
cbug
cbug          yserr = tol * (4.0 * abs (yroot(n)**3) +
cbug     &            2.0 * abs (yroot(n)) *
cbug     &            (abs (c2) + 3.0 * abs (c3 * xroot(n)) +
cbug     &            6.0 * xroot(n)**2))
cbug
cbug          if (abs (xslop) .le. xserr) xslop = 0.0
cbug          if (abs (yslop) .le. yserr) yslop = 0.0
cbug
cbug          slsq = xslop**2 + yslop**2
cbug          xdel = 0.0
cbug          ydel = 0.0
cbug          if (slsq .gt. 0.0) then
cbug            xdel = (fires * yslop - frres * xslop) / slsq
cbug            ydel = (frres * yslop + fires * xslop) / slsq
cbug          endif
cbug          xcorr = xroot(n) + xdel
cbug          ycorr = yroot(n) + ydel
cbug        endif                         ! Tried Newton's method.
cbugc***DEBUG ends.

cbugc***DEBUG begins.
cbug 9921 format ('  Root ',i1,'    r,i=',1p2e22.14, 2x,a8 /
cbug     &        '    P(z)    r,i=',1p2e22.14 /
cbug     &        '    allowed r,i=',1p2e22.14 )
cbug 9922 format ('    Slope   r,i=',1p2e22.14 /
cbug     &        '    Delta   r,i=',1p2e22.14 /
cbug     &        '    Iterate r,i=',1p2e22.14 )
cbug      write ( 3, 9921) n, xroot(n), yroot(n), atype(n), frres, fires,
cbug     &  frerr, fierr
cbug      if ((kbad .eq. 1)      .and.
cbug     &   ((xdel .ne. 0.0)   .or.
cbug     &    (ydel .ne. 0.0))) then
cbug        write ( 3, 9922) xslop, yslop, xdel, ydel, xcorr, ycorr
cbug      endif
cbugc***DEBUG ends.

  390 continue

  395 if ((kbad .eq. 0) .and. (itrun .eq. 2)) itrun = 1
      if ((kbad .eq. 1) .and. (itrun .eq. 0)) itrun = 1
cbugc***DEBUG begins.
cbug 9940 format (/ ' itrun=',i2)
cbug      write ( 3, 9940) itrun
cbugc***DEBUG ends.

      return

c.... End of subroutine aptquar.      (+1 line.)
      end

UCRL-WEB-209832