anchor.diff:      real(kind=x_precision)                                 :: delta = 0.0d0
       !------------------------------------------------------------------------
       delta          = coeff_b**2 - 4_x_precision * coeff_a * coeff_c
  -    if (delta < 0.) then
  -      write(*,*)'No solutions in the R field.'
  -   else
  -      sol_1         = -0.5_x_precision *(coeff_b + sign(sqrt(delta),coeff_b))/coeff_a
  -      if (sol_1 < 0.d-12) then
  -         write(*,*)'There is a problem with the solution of the quadratic equation.'
  -         stop
  -      endif
  -      sol_2 = (coeff_c / (coeff_a * sol_1))
  -      sol = max(sol_1,sol_2)
  -
  -    end if
  -  end subroutine quadratic
  -
  -
  -
  +    if (delta .lt. 0.) then
  +       write(*,*)'No solutions in the R field.'
  +    else
  +       sol_1         = -0.5_x_precision *(coeff_b + sign(sqrt(delta),coeff_b))/coeff_a
  +       !      if (sol_1 .lt. 0.d-12) then
  +       !        write(*,*)'Problem'
  +       !        stop
  +       !endif
  +    sol_2 = (coeff_c / (coeff_a * sol_1))
  +    sol = max(sol_1,sol_2)
  +      
  + end if
  +end subroutine quadratic
  +
  +
  +  
     !-------------------------------------------------------------------------
     !Subroutine in order to compute initial variables rs, rmin, Mdot_0,
     !Sigma_0
  
positive.code:     real(kind=x_precision),intent(in)                        :: T_0
      real(kind=x_precision),intent(in)                        :: rho_0
      real(kind=x_precision),intent(in)                        :: rs
      real(kind=x_precision)                                   :: Scenter = 0._x_precision
  
  
      !-------------------------------------------------------------------------
  
negative.code:   ! Range of surface density [Smin, Smax] and temperature [T_min, T_max]
    ! values should be changed in order to find more coherent results
    !-------------------------------------------------------------------------
    subroutine curve()  
      implicit none
      integer                                                :: i     = 0
      integer                                                :: j     = 0
      real(kind = x_precision)                               :: temp  = 0.0d0
      real(kind = x_precision), parameter                    :: t_min = 3.d-1
      real(kind = x_precision), parameter                    :: t_max = 2.49d0
      !real(kind = x_precision), parameter                    :: t_min = 2.d0
      !real(kind = x_precision), parameter                    :: t_max = 5.49d0
      
      integer,                  parameter                    :: nb_it = 200
      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     = 0.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
  