anchor.diff:      Mdot_0             = param%Mdot
       Omega_0            = sqrt(G * param%M / rs**3 )
       Sigma_0            = Mdot_0 /(Omega_0 * rs**2 * 2_x_precision * pi)
  +    write(*,*)'Mdot_0',Mdot_0
       T_0                = (Mdot_0 * c**2 / (48._x_precision * pi * rs**2 * stefan * &
                            sqrt(27._x_precision) ) )**(1._x_precision/4._x_precision)
       rho_0              = Sigma_0 / (2._x_precision * rs)
  
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:     !-------------------------------------------------------------------------
      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)
      
  