anchor.diff:      if( f_max * f_min .lt. 0.) then
          iteration:do while ( dabs( Smax - Smin ) .ge. eps .and. j .lt. 1000)
             if(f_min * f_dichotomy .lt. 0.) then
  +
  + !   write(*,*)'fmin = ',f_min
  + !   write(*,*)'fmax = ',f_max
  + !   write(*,*)' '
  +
                Smax         = dichotomy
                
                call variables(T, Smax, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e, &
  
positive.code:     !-------------------------------------------------------------------------
      do j              = 1 , 1
         !do j              = 1 , n_cell
         !  r = param%dx*(j-1) + rmin
         r              = 10._x_precision*G*param%M/(c**2)     
         omega          = sqrt(G*param%M/r**3) / Omega_0 
         write(number_of_cell,'(I5.5)') j
         fid = 20 + j
         fname = 'Temperature_Sigma_'//trim(number_of_cell)//'.dat'
         open(fid,file  = fname, status='unknown',action='readwrite') 
         !  write(fid,*)"T              Sigma"
          do i          = 1 , nb_it
            Smin        = 1.0d-50
            Smax        = 1.0d50
            temp        = (t_max-t_min)/(nb_it-1)*(i-1) + t_min
            sigma       = dichotomy(Smin, Smax, eps, temp, omega, sigma_0, Omega_0,rs, T_0, rho_0)
            
  
negative.code:     integer                                                :: i     = 0
      integer                                                :: j     = 0
      real(kind = x_precision)                               :: temp  = 0.0d0
      real(kind = x_precision), parameter                    :: t_min = 5.715d-4
      real(kind = x_precision), parameter                    :: t_max = 1.815d-6
      integer,                  parameter                    :: nb_it = 200
      real(kind = x_precision)                               :: sigma = 0.0d0
      real(kind = x_precision)                               :: Smin  = 0.0d0
      real(kind = x_precision)                               :: Smax  = 0.0d0
      real(kind = x_precision)                               :: eps   = 1.0d-10
      real(kind = x_precision)                               :: omega = 0.0d0
      real(kind = x_precision)                               :: r     = 0.0d0
      real(kind = x_precision)                               :: H     = 0.0d0    
  