anchor.diff:      real(kind=x_precision),intent(out)                      :: Omega_0
       real(kind=x_precision),intent(out)                      :: T_0
       real(kind=x_precision),intent(out)                      :: rho_0
  -    !real(kind=x_precision)                                  :: nu_0
  -
       !------------------------------------------------------------------------
   
       rs                 = 2._x_precision * G * params%M/(c**2)
  
positive.code:     real(kind = x_precision)                               :: Omega_0
      real(kind = x_precision)                               :: rho_0
      real(kind = x_precision)                               :: T_0
      real(kind = x_precision)                               :: temp_real=0.0d0
      real(kind = x_precision)                               :: sigma_real=0.0d0
      character(len = 8)                                     :: number_of_cell
      character(len = 64)                                    :: fname
      integer                                                :: fid
      !------------------------------------------------------------------------
      call initial_variables(rs, rmin, Mdot_0, Sigma_0, Omega_0, T_0, rho_0)
      call display_initial_variables(rs, rmin, Mdot_0, Sigma_0, Omega_0, T_0, rho_0)
      !-------------------------------------------------------------------------
  
negative.code:     ! eps -> Precision
      ! T-> Fixed variable
      !-------------------------------------------------------------------------
      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)
      
      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)
      
      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.'
      endif
  
      if( f_max * f_min .lt. 0.) then
         iteration:do while ( dabs( Smax - Smin ) .ge. eps .and. j .lt. 10000)
  
            
      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)
      
      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)
      
      call variables(T, Scenter, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
           K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f_center,Sigma_0, Omega_0,rs,T_0, rho_0)
  
            
            if(f_min * f_center .gt. 0.) then
               
               Smin         = Scenter
               
            else if (f_max * f_center .gt. 0.) then
               Smax         = Scenter
               
            endif
            
            Scenter       = (Smin + Smax) * 1._x_precision / 2._x_precision
  
            j               = j + 1
      
         end do iteration
      endif
    dichotomy = Scenter
    end function dichotomy
    
    
    !-------------------------------------------------------------------------
    !Subroutine in order to display parameters
    !-------------------------------------------------------------------------
    subroutine display_parameters()
  
      write(*,*)'           Input Parameters             '
      write(*,*)'****************************************'
      write(*,"(' BH_mass     =',1p,E12.4)") params%M
      write(*,"(' Mdot        =',1p,E12.4)") params%Mdot
      write(*,"(' alpha       =',1p,E12.4)") params%alpha
      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
  