anchor.diff:        dichotomy = 0
   
       else if( f_max * f_min < 0.) then
  -      iteration:do while (dabs(Smax - Smin) >= eps .and. j < max_it)
  +      do while (dabs(Smax - Smin) >= eps .and. j < max_it)
   
           call variables(k, T, Smin, f_min, optical_depth, tau_eff)
   
  
positive.code:         S_center = (Smin + Smax) / 2._x_precision
          j = j + 1
  
        end do iteration
  
        dichotomy = S_center
  
  
negative.code:   use mod_variables
    implicit none
  
    real(kind = x_precision), dimension(nb_it) :: temperature ! Temperatures used for the dichotomy
    real(kind = x_precision), dimension(nb_it) :: temp_real   ! Dimensioned version
    real(kind = x_precision)                   :: eps         ! Precision required for the dichotomy
    real(kind = x_precision)                   :: S_min       ! Minimum surface density limit
  