SUBROUTINE CUBIC(DELTA)
C     =======================
C
C     Solution of the cubic equation for determination of
C     the true gradient DELTA
C
C     Input:   A,B   - coefficients; transmitted by COMMON/CUBCON
C              DEL   - DELTA(RAD) - DELTA(ADIAB); also transm. by CUBCON
C     Output:  DELTA - true gradient
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      COMMON/CUBCON/A,B,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD
      PARAMETER (THIRD = 0.333333333333333D0)
       data ipri /0/
C
C     first solve the cubic equation
C         A*X**3 + X**2 + B*X = DEL
C     where X = (DELTA - DELTA(ELEM))**(1/2)
C
      AA=THIRD/A
      BB=B/A
      CC=-DEL/A
      P=BB*THIRD-AA*AA
      Q=AA**3-(BB*AA-CC)/2.D0
      D=Q*Q+P*P*P
      IF(D.GT.0.) THEN
         D=SQRT(D)
         if(d-abs(q).lt.1.e-14*d) then
            SOL=(2.D0*D)**THIRD-AA
          ELSE
            D1=ABS(D-Q)
            D2=ABS(D+Q)
            SOL=D1/(D-Q)*D1**THIRD-D2/(D+Q)*D2**THIRD-AA
          END IF
       ELSE
         COSF=-Q/SQRT(ABS(P*P*P))
         TANF=SQRT(UN-COSF*COSF)/COSF
         FI=ATAN(TANF)*THIRD
         SOL=2.D0*SQRT(ABS(P))*COS(FI)-AA
      END IF
C
C     if the previous formalism gives an unphysical solution
C     x > DEL, then find the physical solution in the range (0, DEL)
C     by a Newton-Raphson solution of the cubic equation
C
      DELDA=SOL*(B+SOL)
      IF(DELDA.GT.DEL.OR.DELDA.LT.0.) THEN
         X0=sol
         J=0
   10    DELX=(DEL-X0*(B+X0+A*X0*X0))/(3.D0*A*X0*X0+2.D0*X0+B)
         X0=X0+DELX
         J=J+1
         IF(ABS(DELX/X0).GT.1D-6.AND.J.LT.50) GO TO 10
         SOL=X0
      END IF
C
C     finally, the actual gradient delta
C
      DELTA=GRDADB+B*SOL+SOL*SOL
      RETURN
      END