anchor.diff:      coeff_b = (-1._x_precision/3._x_precision) * cst_rad * (T * state_0%T_0)**4 / state_0%H_0
       coeff_c = - params%RTM * T * Sigma * state_0%S_0 / (2._x_precision * state_0%H_0**2)
   
  -    call quadratic(coeff_a , coeff_b , coeff_c , H)
  +    delta   = coeff_b**2 - 4._x_precision * coeff_a * coeff_c
  +    H       = -0.5_x_precision * (coeff_b + sign(sqrt(delta),coeff_b)) / coeff_a
   
       rho     = Sigma / H
       cs      = Omega * H
       nu      = params%alpha * cs * H
  -    K_ff    = 6.13d22 * state_0%rho_0 * rho * (state_0%T_0 * T)**(-3.5_x_precision)
  +    K_ff    = 6.13e22_x_precision * state_0%rho_0 * rho * (state_0%T_0 * T)**(-3.5_x_precision)
       K_e     = params%kappa_e
  -    E_ff    = 6.22d20 * (state_0%rho_0 * rho)**2 * sqrt(state_0%T_0 * T)
  +    E_ff    = 6.22e20_x_precision * (state_0%rho_0 * rho)**2 * sqrt(state_0%T_0 * T)
       tau_eff = 0.5_x_precision * sqrt(K_e * K_ff) * Sigma * state_0%S_0
   
       !-------------------------------------------------------------------------
  
positive.code:     integer                                                :: i     = 0
      integer                                                :: j     = 0
      integer                                                :: k     = 0
      integer                                                :: file_length = 1        ! Number of lines of the data saving file
  
   ! For the temperature and s
      real(kind = x_precision)                               :: temp  = 0.0d0
      real(kind = x_precision)                               :: Smin  = 0.0d0
      real(kind = x_precision)                               :: Smax  = 0.0d0
      real(kind = x_precision)                               :: Sigma = 0.0d0          ! S
  
      real(kind = x_precision)                               :: omega = 0.0d0          ! Angular velocity
      real(kind = x_precision)                               :: r     = 0.0d0          ! Radius
      real(kind = x_precision)                               :: f     = 0.0d0          ! Q+ - Q-
      real(kind=x_precision)                                 :: tau_eff = 0.0d0        ! Effective optical depth
      integer                                                :: optical_depth = 0      ! Indicator for the optical thickness
  
   ! For the separate curves given the thickness
      real(kind = x_precision),dimension(nb_it)              :: temp_t_thick
  
negative.code:   end subroutine write_critical_points
  
  
    !-------------------------------------------------------------------------
    !Subroutine in order to compute variables H, rho, cs, nu, Q_plus, Q_minus,
    !K_ff, K_e, tau_eff, E_ff,Fz given T, Sigma and Omega
  