anchor.diff:  
       j = 0
   
  -    call variables(T, Smin, Omega, f_min, optical_depth, tau_eff)
  -
  -    call variables(T, Smax, Omega, f_max, optical_depth, tau_eff)
  +    call variables(k, T, Smin, f_min, optical_depth, tau_eff)
  +
  +    call variables(k, T, Smax, f_max, optical_depth, tau_eff)
   
       !write(*,*)'fmin = ',f_min
       !write(*,*)'fmax = ',f_max
  
positive.code:         Smax          = S_max
  
          ! S found with the dichotomy approach
          sigma_t_thick(i) = dichotomy(Smin, Smax, eps, temperature(i), omega, optical_depth)
  
          ! Optical thin case (tau < 1)
          optical_depth = 0
  
negative.code:     real(kind = x_precision)             :: Q_minus
      !------------------------------------------------------------------------
  
      Omega   = x_state%Omega(k)
      coeff_a = (Omega * state_0%Omega_0)**2 * Sigma * state_0%S_0 / 2._x_precision
      coeff_b = (-1._x_precision/3._x_precision) * cst_rad * (T * state_0%T_0)**4 / state_0%H_0
      coeff_c = - params%RTM * T * Sigma * state_0%S_0 / (2._x_precision * state_0%H_0**2)
  