anchor.diff:      !------------------------------------------------------------------------
   
       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)
  -
       delta   = coeff_b**2 - 4._x_precision * coeff_a * coeff_c
       H       = -0.5_x_precision * (coeff_b + sign(sqrt(delta),coeff_b)) / coeff_a
   
  
positive.code:       dichotomy = 0
  
      else if( f_max * f_min < 0.) then
        iteration:do while (dabs(Smax - Smin) >= eps .and. j < max_it)
  
          call variables(k, T, Smin, f_min, optical_depth, tau_eff)
  
  
negative.code:     integer, intent(in) :: optical_depth
      !-------------------------------------------------------------------------
      integer                  :: j
      real(kind = x_precision) :: T ! Temperature used for the dichotomy
      real(kind = x_precision) :: tau_eff
      real(kind = x_precision) :: Smin     ! Lowest point
      real(kind = x_precision) :: Smax     ! Highest point
      real(kind = x_precision) :: S_center ! Intermediate point
      real(kind = x_precision) :: f_min    ! Q+ − Q− at the lowest point
      real(kind = x_precision) :: f_max    ! Q+ − Q− at the highest point
      real(kind = x_precision) :: f_center ! Q+ − Q− at the intermediate point
      !-------------------------------------------------------------------------
  
      T = temperature(i)
  