anchor.diff:      do j              = 1 , 1
          !do j              = 1 , n_cell
          !  r = (param%rmax-rmin)/(n_cell-1)*(j-1) + rmin
  -       r              = 10._x_precision*G*param%M/(c**2)     
  -       omega          = sqrt(G*param%M/r**3) / Omega_0 
  +       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') 
  +       open(fid,file  = fname, status='unknown',action='readwrite')
          !  write(fid,*)"T     &         Sigma"
           do i          = 1, nb_it
             Smin        = 1d-3
             Smax        = 1d51
  -          
  +
             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)
  -          
  +
             call variables(temp, sigma, omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
                 K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f,Sigma_0, Omega_0,rs, T_0, rho_0)
   
  -          
  +
             call display_variables(temp,Omega,r, sigma, H, rho, cs, nu, Q_plus, Q_minus,&
                  K_ff, K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f)
  -          
  -          temp_reel   = temp * T_0
  -          sigma_reel  = sigma * sigma_0
  -          write(fid,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sigma_reel,temp_reel, tau_eff
  +
  +          temp_real   = temp * T_0
  +          sigma_real  = sigma * sigma_0
  +          write(fid,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sigma_real,temp_real, tau_eff
          enddo
          close(fid)
       enddo
  -    
  +
     end subroutine curve
   
   
  
positive.code:     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
      nu_0                 = 2._x_precision * rs**2 * Omega_0/3._x_precision
      
      
      !-------------------------------------------------------------------------
      !Select the case for the opticaly depth to compute Fz
      !-------------------------------------------------------------------------
  
      optical_depth        = 1
      
     ! if (tau_eff .ge. 1.)  then
     !    optical_depth     = 1
     ! else
  
negative.code:     write(*,"(' alpha       =',1p,E12.4)") param%alpha
      write(*,"(' X           =',1p,E12.4)") param%X
      write(*,*)'****************************************'
  
      read(*,*)
    end subroutine display_parameters
  
  
    !-------------------------------------------------------------------------
    !Subroutine in order to display initial variables
    !-------------------------------------------------------------------------
    subroutine display_initial_variables(rs, rmin, Mdot_0, Sigma_0, Omega_0, T_0, rho_0)
      implicit none
  
      real(kind = x_precision), intent(in)                     :: rs
      real(kind = x_precision), intent(in)                     :: rmin
      real(kind = x_precision), intent(in)                     :: Mdot_0
  