anchor.diff:       iteration:do while ( dabs( Smax - Smin ) .ge. eps .and. j .lt. N)
           if( f(T,Smin,omega) * f(T,dichotomy,omega) .lt. 0.) then
              Smax         = dichotomy
  +           write(*,*)'smin=',smin,' smax=', smax
           else
              Smin         = dichotomy
  +           write(*,*)'smin=',smin,' smax=', smax
           endif
           dichotomy       = (Smin+Smax)/2.0d0
           j               = j + 1
  
positive.code:     use mod_variables
      implicit none
      
    integer                                                  ::j=0,N=100   
    real(kind=x_precision),intent(inout)                     ::Smin,Smax   
    real(kind=x_precision),intent(in)                        ::eps         
    real(kind=x_precision),intent(in)                        ::T           
  
negative.code:   end if
  
    select case(optical_depth)
    case(1)
       Fz = 4._x_precision * c**2 * T**4/(27. * sqrt(3._x_precision) * (K_ff + K_e) * Sigma * Sigma_0)
    case default
       Fz = 6.22d20 * 2._x_precision / state_0%Omega_0 * state_0%rho_0 *  H * rho**2 * sqrt(T*state_0%T_0)
  