anchor.diff:      real(kind = x_precision)                               :: rmin
       real(kind = x_precision)                               :: Mdot_0
       real(kind = x_precision)                               :: Sigma_0
  +    character(len = 8)                                     :: number_of_cell
  +    character(len = 64)                                    :: fname
  +    integer                                                :: fid
       type(parameters)                                       :: param
  -    !----------------------------------------------------------------------
  -    
  +    !------------------------------------------------------------------------
       call get_parameters(param)
       call display_parameters()
       call initial_variables(rs, rmin, Mdot_0, Sigma_0)
  
positive.code:     !-------------------------------------------------------------------------
      ! Test for 1 value of r
      !-------------------------------------------------------------------------
      r              = 6*G*param%M/(c**2)     
      omega          = G*param%M/r**3
  
      open(18,file   = 'Temperature_Sigma.dat',status='unknown',action='readwrite') 
      write(18,*)"T Sigma"
      do i           = 1 , nb_it
         Smin        = 1.0d-5
         Smax        = 1.0d-0
         temp        = (t_max-t_min)/(nb_it-1)*(i-1) + t_min
         sigma       = dichotomy(Smin, Smax, eps, temp, omega)
         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)
         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)
         
         write(18,*)temp," ",sigma
      enddo
      close(18)
    end subroutine curve
  
  
  
negative.code:     use mod_variables
      implicit none
      
      integer                                                  :: j = 0   
      real(kind=x_precision),intent(inout)                     :: Smin,Smax   
      real(kind=x_precision),intent(in)                        :: eps         
      real(kind=x_precision),intent(in)                        :: T           
      real(kind=x_precision),intent(in)                        :: omega
      real(kind=x_precision)                                   :: H         
      real(kind=x_precision)                                   :: rho           
      real(kind=x_precision)                                   :: cs
      real(kind=x_precision)                                   :: nu        
      real(kind=x_precision)                                   :: Q_plus           
      real(kind=x_precision)                                   :: Q_minus
      real(kind=x_precision)                                   :: K_ff        
      real(kind=x_precision)                                   :: K_e           
      real(kind=x_precision)                                   :: tau_eff
      real(kind=x_precision)                                   :: P_rad        
      real(kind=x_precision)                                   :: P_gaz
      real(kind=x_precision)                                   :: f_min          
      real(kind=x_precision)                                   :: f_max
      real(kind=x_precision)                                   :: f_dichotomy
      real(kind=x_precision)                                   :: E_ff
      real(kind=x_precision)                                   :: Fz
      real(kind=x_precision)                                   :: Sigma_0
      
      !-------------------------------------------------------------------------
      ! N-> Number of iterations for the dichotomy
      ! Smin, Smax -> Starting range
      ! 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)
      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)
      call variables(T, dichotomy, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
           K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f_dichotomy,Sigma_0)
      
      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)
            if(f_min * f_dichotomy .lt. 0.) then
               Smax         = dichotomy
            else
               Smin         = dichotomy
            endif
            dichotomy       = (Smin+Smax)/2.0d0
            j               = j + 1
         end do iteration
      endif
    end function dichotomy
    
    
    !-------------------------------------------------------------------------
    !Subroutine in order to display parameters
    !-------------------------------------------------------------------------
    subroutine display_parameters()
      type(parameters)                            :: param
      call get_parameters(param)
      write(*,*)'           Input Parameters             '
      write(*,*)'****************************************'
      write(*,"(' BH_mass     =',1p,E12.4)") param%M
      write(*,"(' Mdot        =',1p,E12.4)") param%Mdot
      write(*,"(' rmax        =',1p,E12.4)") param%rmax
      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)
      implicit none
      real(kind = x_precision), intent(in)                     :: rs
      real(kind = x_precision), intent(in)                     :: rmin
      real(kind = x_precision), intent(in)                     :: Mdot_0
      real(kind = x_precision), intent(in)                     :: Sigma_0
      !-----------------------------------------------------------------------
  
      write(*,*)'           Initial Variables            '
      write(*,*)'****************************************'
      write(*,"(' Temp_0      =',1p,E12.4)") state_0%T_0
      write(*,"(' Simga_0     =',1p,E12.4)") Sigma_0
      write(*,"(' Omega_0     =',1p,E12.4)") state_0%Omega_0
      write(*,"(' H_0         =',1p,E12.4)") state_0%H_0
      write(*,"(' Mdot_0      =',1p,E12.4)") Mdot_0
      write(*,"(' rho_0       =',1p,E12.4)") state_0%rho_0
      write(*,"(' rmin        =',1p,E12.4)") rmin
      write(*,"(' rs          =',1p,E12.4)") rs
      write(*,*)'****************************************'
      read(*,*)
    end subroutine display_initial_variables
  
  
    !-------------------------------------------------------------------------
    !Subroutine in order to display variables
    !-------------------------------------------------------------------------
    subroutine 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)
      implicit none
      real(kind = x_precision), intent(in)                     :: temp
      real(kind = x_precision), intent(in)                     :: Omega
      real(kind = x_precision), intent(in)                     :: r
      real(kind = x_precision), intent(in)                     :: sigma
      real(kind = x_precision), intent(in)                     :: H     
      real(kind = x_precision), intent(in)                     :: rho        
      real(kind = x_precision), intent(in)                     :: cs  
      real(kind = x_precision), intent(in)                     :: nu    
      real(kind = x_precision), intent(in)                     :: Q_plus          
      real(kind = x_precision), intent(in)                     :: Q_minus  
      real(kind = x_precision), intent(in)                     :: K_ff         
      real(kind = x_precision), intent(in)                     :: K_e            
      real(kind = x_precision), intent(in)                     :: tau_eff  
      real(kind = x_precision), intent(in)                     :: P_rad        
      real(kind = x_precision), intent(in)                     :: P_gaz    
      real(kind = x_precision), intent(in)                     :: f        
      real(kind = x_precision), intent(in)                     :: E_ff     
      real(kind = x_precision), intent(in)                     :: Fz  
      !------------------------------------------------------------------------
      
      write(*,*)'               Variables                '
      write(*,*)'****************************************'
      write(*,"(' r           =',1p,E12.4)") r
      write(*,"(' Temp        =',1p,E12.4)") temp
      write(*,"(' Sigma       =',1p,E12.4)") sigma
      write(*,"(' Omega       =',1p,E12.4)") Omega
      write(*,"(' rho         =',1p,E12.4)") rho
      write(*,"(' sigma       =',1p,E12.4)") sigma
      write(*,"(' H           =',1p,E12.4)") H
      write(*,"(' rho         =',1p,E12.4)") rho
      write(*,"(' cs          =',1p,E12.4)") cs
      write(*,"(' nu          =',1p,E12.4)") nu
      write(*,"(' Q plus      =',1p,E12.4)") Q_plus
      write(*,"(' Q minus     =',1p,E12.4)") Q_minus
      write(*,"(' K_ff        =',1p,E12.4)") K_ff
      write(*,"(' K_e         =',1p,E12.4)") K_e
      write(*,"(' tau_eff     =',1p,E12.4)") tau_eff
      write(*,"(' P_gaz       =',1p,E12.4)") P_gaz
      write(*,"(' P_rad       =',1p,E12.4)") P_rad
      write(*,"(' E_ff        =',1p,E12.4)") E_ff
      write(*,"(' delta Q     =',1p,E12.4)") f
      write(*,"(' Fz          =',1p,E12.4)") Fz
      write(*,*)'****************************************'
      
    end subroutine display_variables
    
  
  end module mod_S_curve
  