#include "dims.h"
      SUBROUTINE VQUART(A,ROOT,W1,W2,W3,LEN,I1,I2,j)
      implicit none
C     ****
C     ****     VECTOR PROCEDURE DETERMINES +IVE ROOTS OF EQUATIONS
C     ****     A(I,5)*X**4+4.*A(I,4)*X**3+6.*A(I,3)*X**2+4.*A(I,2)*X
C     ****     +A(I,1)=0., WHERE I=1,LEN
C     ****
C     ****     PROCEDURE IS SPECIFICALLY DESIGNED FOR REAL QUARTICS WITH
C     ****     REAL ROOTS ONLY ONE OF WHICH IS POSITIVE
C     ****
#include "params.h"
!
! Args:
      integer,intent(in) :: i1,i2,len,j
      real,intent(in) :: A(LEN,5)
      real,intent(out) :: ROOT(1),W1(1),W2(1),W3(1)
!
! Local:
      real :: e
      integer :: i
C     ****
C     ****     COEFFICIENTS OF QUARTICS TRANSMITTED IN ARRAY A(LEN,5)
C     ****     POSITIVE ROOTS RETURNED IN ROOT(LEN)
C     ****     W1(LEN),W2(LEN),W3(LEN) ARE WORKING ARRAYS
C     ****     LEN IS FIRST DIMENSION OF A
C     ****     I1,I2 GIVE RANGE OVER WHICH QUARTICS ARE TO BE SOLVED
C     ****
!     E=1.E-1200    ! this was used on Cray
      E=1.E-300     ! largest exponent allowed on sgi is ~308
!
      DO 1 I=I1,I2
C     ****
C     ****     W1=CH
C     ****
      W1(I)=-(A(I,5)*A(I,1)-4.*A(I,4)*A(I,2)+3.*A(I,3)**2)/12.
C     ****
C     ****     W2=CG
C     ****
      W2(I)=(A(I,5)*(A(I,3)*A(I,1)-A(I,2)**2)-A(I,4)*(A(I,4)*A(I,1)-
     1A(I,2)*A(I,3))+A(I,3)*(A(I,4)*A(I,2)-A(I,3)**2))/4.
C     ****
C     ****     ROOT=RLAM=-2.*REAL((.5*(CMPLX(CG,0.)+CSQRT(CMPLX(CG**2+4.
C     ****     *CH**3+E,0.)))+CMPLX(E,0.))**(1./3.))
C     ****
      ROOT(I)=-2.*REAL((.5*(CMPLX(W2(I),0.)+CSQRT(CMPLX(W2(I)**2+4.*
     1W1(I)**3+E,0.)))+CMPLX(E,0.))**(1./3.))
C     ****
C     ****     W1=P=SQRT(A(5)*RLAM+A(4)**2-A(5)*A(3)+E)
C     ****
      W1(I)=A(I,5)*ROOT(I)+A(I,4)**2-A(I,5)*A(I,3)+E
      if (w1(i) < 0.) w1(i) = 0.
      W1(I)=SQRT(W1(I))
C     ****
C     ****     W2=Q=SQRT((2.*RLAM+A(3))**2-A(5)*A(1)+E)
C     ****
      W2(I)=SQRT((2.*ROOT(I)+A(I,3))**2-A(I,5)*A(I,1)+E)
C     ****
C     ****     W3=PQ=2.*A(4)*RLAM+A(4)*A(3)-A(5)*A(2)+E
C     ****
      W3(I)=2.*A(I,4)*ROOT(I)+A(I,4)*A(I,3)-A(I,5)*A(I,2)+E
C     ****
C     ****     W1=P=SIGN(P,Q*PQ)
C     ****
      W1(I)=SIGN(W1(I),W2(I)*W3(I))
C     ****
C     ****     W3=P-A4
C     ****
      W3(I)=W1(I)-A(I,4)
C     ****
C     ****     FINAL EVALUATION OF REQUIRED ROOT
C     ****
      ROOT(I)=(W3(I)+SQRT(W3(I)**2-A(I,5)*(A(I,3)+2.*ROOT(I)-W2(I))))/
     1A(I,5)
    1 CONTINUE

!     call addfsech('W1'    ,' ',' ',w1,zimxp,zkmxp,zkmx,j)
!     call addfsech('W2'    ,' ',' ',w1(len+1),zimxp,zkmxp,zkmx,j)
!     call addfsech('W3'    ,' ',' ',w1(len*2+1),zimxp,zkmxp,zkmx,j)
!     call addfsech('VQROOT'    ,' ',' ',root,zimxp,zkmxp,zkmx,j)

      RETURN
      END
C