anchor.diff:      cs                   = Omega * H
       nu                   = params%alpha * cs * H
       K_ff                 = 6.13d22 * rho_0 * rho * (T_0 * T)**(-3.5_x_precision)
  -    K_e                  =  params%kappa_e
  +    K_e                  = params%kappa_e
       E_ff                 = 6.22d20 * (rho_0 * rho)**2 * sqrt(T_0 * T)
       tau_eff              = 0.5_x_precision * sqrt(K_e * K_ff) * Sigma * Sigma_0
   
  
positive.code:     real(kind = x_precision),dimension(nb_it)              :: sigma_real_1=0.0d0
      real(kind = x_precision),dimension(nb_it)              :: temp_real_0 =0.0d0
      real(kind = x_precision),dimension(nb_it)              :: sigma_real_0=0.0d0
  
      character(len = 8)                                     :: number_of_cell
      character(len = 64)                                    :: fname
      integer                                                :: fid
      integer                                                :: n_cell = 1
      !------------------------------------------------------------------------
  
      call display_parameters()
  
negative.code:     call variables(T, Smax, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e, &
           tau_eff, P_rad, P_gaz,E_ff,Fz,f_max,Sigma_0, Omega_0,rs,T_0, rho_0, optical_depth)
  
    !   write(*,*)'fmin = ',f_min
    !   write(*,*)'fmax = ',f_max
  
      if ( f_max * f_min .gt. 0.) then
     !    write(*,*)'This function image does not switch its sign in this particular interval.'
         dichotomy = 0
  
      else if( f_max * f_min .lt. 0.) then
         iteration:do while ( dabs( Smax - Smin ) .ge. eps .and. j .lt. 10000)
  
  
  