      SUBROUTINE CBCRT (AC, RZ)
C     ******************************************************************
C     FORTRAN SUBROUTINE FOR ROOTS OF CUBIC POLYNOMIAL
C     ******************************************************************
C     AC = POLYNOMIAL COEFFICIENTS (4-ARRAY)
C     RZ = POLYNOMIAL ROOTS        (6-ARRAY)
C
      DIMENSION AC(4), RZ(6)
  001 P=AC(3)/AC(4)
      Q=AC(2)/AC(4)
      R=AC(1)/AC(4)
      A=(3.0*Q-P*P)/3.0
      B=(2.0*P*P*P-9.0*P*Q+27.0*R)/27.0
      SQ=B*B/4.0+A*A*A/27.0
      RQ=SQRT(ABS(SQ))
      IF(SQ)010,002,004
  002 IF(B)004,003,004
  003 RA=0.0
      RB=0.0
      GO TO 008
  004 IF(B)005,005,006
  005 RA=CBRT(-0.5*B+RQ)
      GO TO 007
  006 RA=CBRT(-0.5*B-RQ)
  007 RB=-A/(3.0*RA)
  008 RZ(1)=RA+RB-P/3.0
      RZ(2)=0.0
      RZ(3)=-0.5*(RA+RB)-P/3.0
      RZ(4)=+0.8660254037844386*(RA-RB)
      RZ(5)=-0.5*(RA+RB)-P/3.0
      RZ(6)=-0.8660254037844386*(RA-RB)
      A=RZ(1)*RZ(1)
      B=RZ(3)*RZ(5)-RZ(4)*RZ(6)
      IF(A.GE.B)RETURN
  009 RZ(1)=-R/B
      RETURN
  010 F=ATAN2(RQ,-0.5*B)
      CF=COS(F/3.0)
      SF=SIN(F/3.0)
      RQ=SQRT(-A/3.0)
      RZ(1)=2.0*RQ*CF-P/3.0
      RZ(2)=0.0
      RZ(3)=-RQ*(CF+1.7320508075688773*SF)-P/3.0
      RZ(4)=0.0
      RZ(5)=-RQ*(CF-1.7320508075688773*SF)-P/3.0
      RZ(6)=0.0
      IF(ABS(RZ(1)).GE.ABS(RZ(3)))GO TO 011
      IF(ABS(RZ(1)).GE.ABS(RZ(5)))GO TO 013
      RZ(1)=-R/(RZ(3)*RZ(5))
      RETURN
  011 IF(ABS(RZ(3)).GE.ABS(RZ(5)))GO TO 014
  012 RZ(3)=-R/(RZ(5)*RZ(1))
      RETURN
  013 IF(ABS(RZ(5)).GE.ABS(RZ(3)))GO TO 012
  014 RZ(5)=-R/(RZ(1)*RZ(3))
      RETURN
      END
