PROGRAM Cubic !https://en.m.wikipedia.org/wiki/Cubic_equation#General_cubic_formula
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  INTEGER :: k
  REAL(dp) :: a, b, c, d
  COMPLEX(dp) :: Delta_0, Delta_1, radical, CC, CBRT, x
  COMPLEX(dp), PARAMETER :: xi = CMPLX(-1.0_dp / 2.0_dp, SQRT(3.0_dp) / 2.0_dp, dp) !primitive cube root of unity
  EXTERNAL CBRT !function
  PRINT *, " "
  WRITE(6, '(a6)', advance='no') "  a = "
  READ *, a
  IF (a > 0.0_dp .OR. a < 0.0_dp) THEN !a /= 0
    GO TO 1
  ELSE
    PRINT *, " "
    STOP
  END IF
1 WRITE(6, '(a6)', advance='no') "  b = "
  READ *, b
  WRITE(6, '(a6)', advance='no') "  c = "
  READ *, c
  WRITE(6, '(a6)', advance='no') "  d = "
  READ *, d
  PRINT *, " "
  Delta_0 = b**2 - 3.0_dp * a * c
  Delta_1 = 2.0_dp * b**3 - 9.0_dp * a * b * c + 27.0_dp * a**2 * d
  radical = SQRT(Delta_1**2 - 4.0_dp * Delta_0**3)
  CC = CBRT((Delta_1 + radical) / 2.0_dp)
  roots: DO k = 0, 2 !three roots
    x = -(b + xi**k * CC + Delta_0 / (xi**k * CC)) / (3.0_dp * a) !k = 0, always real
    IF (INT(x / CONJG(x)) == 1) THEN !are you for real?
      WRITE(6, '(f20.15)') REALPART(x)
    ELSE
      WRITE(6, '(2f20.15, a1)') REALPART(x), IMAGPART(x), "i"
    END IF
  END DO roots
  PRINT *, " "
END PROGRAM Cubic

COMPLEX(dp) FUNCTION CBRT(x)
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  COMPLEX(dp) :: x
  REAL(dp), PARAMETER :: ONE_THIRD = 1.0_dp / 3.0_dp
  CBRT = x**ONE_THIRD
END FUNCTION CBRT