anchor.diff:          end if
           sol_2 = (coeff_c / (coeff_a * sol_1))
           sol = max(sol_1,sol_2)
  +        if (sol /= sol_1) then
  +          write(*,*)'The solution was the other one !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
  +        end if
         end if
   
       end if
  
positive.code:   real(kind=x_precision) function dichotomy(Smin, Smax, eps, T, omega, optical_depth)
      implicit none
  
      real(kind=x_precision),intent(inout)                     :: Smin,Smax
      real(kind=x_precision),intent(in)                        :: eps
      real(kind=x_precision),intent(in)                        :: T
      real(kind=x_precision),intent(in)                        :: omega
      integer,intent(in)                                       :: optical_depth
  
      real(kind=x_precision)                                   :: f_min
      real(kind=x_precision)                                   :: f_max
      real(kind=x_precision)                                   :: f_center
      real(kind=x_precision)                                   :: S_center = 0._x_precision
      integer                                                  :: j = 0
  
      !-------------------------------------------------------------------------
      ! N-> Number of iterations for the dichotomy
  
negative.code:   subroutine quadratic(coeff_a, coeff_b, coeff_c, sol)
      implicit none
  
      real(kind=x_precision), intent(in)  :: coeff_a,coeff_b,coeff_c
      real(kind=x_precision), intent(out) :: sol
      real(kind=x_precision)              :: sol_1 = 0.0d0
      real(kind=x_precision)              :: sol_2 = 0.0d0
      real(kind=x_precision)              :: delta = 0.0d0
      !------------------------------------------------------------------------
  
      delta          = coeff_b**2 - 4_x_precision * coeff_a * coeff_c
  