anchor.diff:      !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
  +    optical_depth        = 0
  +
  +  !  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 =  E_ff * H 
  +       ! 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
  +  ! Dichotomic 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
  
positive.code: 
      call get_parameters(param)
     
      coeff_a              = (Omega**2 * Omega_0**2 * Sigma * Sigma_0)/2.
      coeff_b              = (-1._x_precision/3._x_precision) * cst_rad*T**4 * T_0**4
      coeff_c              = (-1._x_precision * param%RTM * T * T_0* Sigma * Sigma_0)/2._x_precision
  
      call quadratic(coeff_a , coeff_b , coeff_c , H)
  
      H                    = H /rs
      rho                  = Sigma / H
      P_rad                = T**4
      P_gaz                = rho * T
      cs                   = Omega * H
      nu                   = param%alpha * cs * H
      K_ff                 = 6.13d22 * rho_0 * rho * (T_0 * T)**(-3.5_x_precision)
      K_e                  = 0.2_x_precision * (1_x_precision + param%X)
      E_ff                 = 6.22d20 * (rho_0 * rho)**2 * sqrt(T_0 * T)
      tau_eff              = 0.5_x_precision * sqrt(K_e * K_ff) * Sigma * Sigma_0    
      
  
negative.code:     real(kind=x_precision)                                   :: rho           
      real(kind=x_precision)                                   :: cs
      real(kind=x_precision)                                   :: nu        
      real(kind=x_precision)                                   :: Q_plus          
      real(kind=x_precision)                                   :: Q_minus
      real(kind=x_precision)                                   :: K_ff        
      real(kind=x_precision)                                   :: K_e           
  