anchor.diff:  
       implicit none
       real(kind=x_precision), intent(in)                     :: coeff_a,coeff_b,coeff_c
  -    real(kind=x_precision), intent(inout)                  :: sol
  -    real(kind=x_precision)                                 :: delta=0.0d0
  +    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
       if (delta .lt. 0.) then
         write(*,*)'No solutions in the R field.'
  -    else
  -      sol          = -0.5d0 *(coeff_b + sign(sqrt(delta),coeff_b))/coeff_a
  -    end if
  +   else
  +      sol_1=(-1.*coeff_b+sqrt(delta))/(2.*coeff_a)
  +      sol_2=(-1.*coeff_b-sqrt(delta))/(2.*coeff_a)
  +   end if
  +   if (sol_1 .gt. 0. ) then
  +      sol=sol_1
  +   else
  +      sol=sol_2
  +   end if
  +      !sol          = -0.5_x_precision *(coeff_b + sign(sqrt(delta),coeff_b))/coeff_a
  +    !end if
     end subroutine quadratic
   
   
  
positive.code:     !Select the case for the opticaly depth to compute Fz
      !-------------------------------------------------------------------------
  
      optical_depth        = 1
  
      !if (tau_eff .ge. 1.)  then
      !   optical_depth     = 1
      !else
      !   optical_depth     = 0
      !end if
      
       select case(optical_depth)
        case(1)
           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 * rs * E_ff * H / (Omega_0 * Sigma_0)
           Fz = 6.22d20 * 2._x_precision / Omega_0 * rho_0 * H * rho**2 * sqrt(T * T_0)
       end select
      Q_minus             = 2._x_precision * Fz /Sigma
      Q_plus              = 9._x_precision /4._x_precision * nu * Omega**2
      f                   = Q_plus - Q_minus
      
    end subroutine variables
    
    
    !-------------------------------------------------------------------------
    ! Function in order to determine the change of sign in a given interval
    ! [Smin,Smax] with an epsilon precision
    !-------------------------------------------------------------------------
    real(kind=x_precision) function dichotomy(Smin, Smax, eps, T, omega, sigma_0, Omega_0,rs, T_0, rho_0)
      use mod_read_parameters
  
negative.code:     write(*,"(' rmin        =',1p,E12.4)") rmin
      write(*,"(' rs          =',1p,E12.4)") rs
      write(*,*)'****************************************'
  
      read(*,*)
    end subroutine display_initial_variables
  
  