anchor.diff:      real(kind=x_precision),intent(in)                        :: T_0
       real(kind=x_precision),intent(in)                        :: rho_0
       real(kind=x_precision),intent(in)                        :: rs
  +    integer,intent(in)                                       :: optical_depth
   
       real(kind=x_precision)                                   :: H
       real(kind=x_precision)                                   :: rho
  
positive.code:           temp_real_0(i)   = log10(temp * T_0)
            sigma_real_0(i)  = log10(sigma * sigma_0)
  
         enddo
  
  
         temp = 0.0d0
         sigma = 0.0d0
  
  !********For the optically thick medium********
          do i          = 1, nb_it
            Smin        = 1d-2
            Smax        = 1d51
            t_min       = 3.d-1
            t_max       = 2.49d0
  
            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)
  
            optical_depth = 1
  
            call variables(temp, sigma, omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
  
negative.code:     dichotomy             = (Smin+Smax)/2.
      j = 0
      call variables(T, Smin, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e,&
           tau_eff, P_rad, P_gaz,E_ff,Fz,f_min,Sigma_0, Omega_0,rs,T_0,rho_0,optical_depth)
  
      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
  