anchor.diff:    Q_minus             = 2._x_precision * Fz /Sigma
     Q_plus              = 9._x_precision /4._x_precision * nu * Omega**2
     f                   = Q_plus - Q_minus
  -  end function f
  +  end subroutine variables
   
   
     !-------------------------------------------------------------------------
  
positive.code:   real(kind=x_precision),intent(in)                        :: eps         
    real(kind=x_precision),intent(in)                        :: T           
    real(kind=x_precision),intent(in)                        :: omega
    !-------------------------------------------------------------------------
    ! N-> Number of iterations for the dichotomy
    ! Smin, Smax -> Starting range
  
negative.code:   !-------------------------------------------------------------------------
    subroutine curve()  
      implicit none
      integer                                                :: i     = 0
      real(kind = x_precision)                               :: temp  = 0.0d0
      real(kind = x_precision), parameter                    :: t_min = 1.2d-2
      real(kind = x_precision), parameter                    :: t_max = 1.0d0
      integer,                  parameter                    :: nb_it = 10
      real(kind = x_precision)                               :: sigma = 0.0d0
      real(kind = x_precision)                               :: Smin  = 0.0d0
      real(kind = x_precision)                               :: Smax  = 0.0d0
      real(kind = x_precision)                               :: eps   = 1.0d-4
      real(kind = x_precision)                               :: omega = 0.0d0
      real(kind = x_precision)                               :: r     = 5.0d0
      real(kind = x_precision)                               :: H     = 0.0d0    
      real(kind = x_precision)                               :: rho   = 0.0d0       
      real(kind = x_precision)                               :: cs    = 0.0d0
      real(kind = x_precision)                               :: nu    = 0.0d0  
      real(kind = x_precision)                               :: Q_plus = 0.0d0        
      real(kind = x_precision)                               :: Q_minus= 0.0d0
      real(kind = x_precision)                               :: K_ff   = 0.0d0    
      real(kind = x_precision)                               :: K_e    = 0.0d0      
      real(kind = x_precision)                               :: tau_eff= 0.0d0
      real(kind = x_precision)                               :: P_rad  = 0.0d0    
      real(kind = x_precision)                               :: P_gaz  = 0.0d0
      real(kind = x_precision)                               :: f      = 0.0d0
      real(kind = x_precision)                               :: E_ff   = 0.0d0
      real(kind = x_precision)                               :: Fz     = 0.0d0
      real(kind = x_precision)                               :: rs
      real(kind = x_precision)                               :: rmin
      real(kind = x_precision)                               :: Mdot_0
      real(kind = x_precision)                               :: Sigma_0
      type(parameters)                                       :: param
      !----------------------------------------------------------------------
      
      call get_parameters(param)
      call display_parameters()
      call initial_variables(rs, rmin, Mdot_0, Sigma_0)
      call display_initial_variables(rs, rmin, Mdot_0, Sigma_0)
      !-------------------------------------------------------------------------
      ! Test for 1 value of r
      !-------------------------------------------------------------------------
  