anchor.diff:  
       case(1)
   
  -       Fz = 2._x_precision * c**2 * T**4 /(27._x_precision * sqrt(3._x_precision) &
  -            * (K_ff + K_e) * Sigma * state_0%S_0)
  +       Fz = 4._x_precision * c**2 * T**4 /(27._x_precision * sqrt(3._x_precision) &
  +            * (K_ff + K_e) * Sigma * Sigma_0)
       case (0)
   
  -       Fz = 4._x_precision * state_0%H_0 * E_ff * H / (state_0%Omega_0 * state_0%S_0)
  +       Fz = 4._x_precision * rs * E_ff * H / (Omega_0 * Sigma_0)
   
       end select
   
  -    Q_plus              = 3._x_precision  * state_0%H_0**2 * nu * Omega**2 * state_0%Omega_0**2
  +    Q_plus              = 3._x_precision  * rs**2 * nu * Omega**2 * Omega_0**2
       Q_minus             = Fz  / Sigma
   
       f                   = Q_plus - Q_minus
  
positive.code:     real(kind = x_precision),intent(out)                     :: f
      !------------------------------------------------------------------------
  
      coeff_a              = (Omega**2 * state_0%Omega_0**2 * Sigma * state_0%S_0)/2._x_precision
      coeff_b              = (-1._x_precision/3._x_precision) * cst_rad*T**4 * state_0%T_0**4 / state_0%H_0
      coeff_c              = (-1._x_precision * 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)
  
  
negative.code: 
         call build_s_curve(sigma_real_1, sigma_real_0, temp_real_1,nb_it, index_scp, sigma_real, temp_real)
  
          call display_critical_points(sigma_c_thin, temp_c_thin,sigma_c_thick, temp_c_thick, k)
  
  
         
            do l          = 1, nb_it
               write(fid_3,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sigma_real(l),temp_real(l)
            enddo
         
         
  