anchor.diff:      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
  -  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 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 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 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(*,*)'****************************************'
  -  read(*,*)
  -
  -end subroutine display_variables
  -
  +    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
  
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:   subroutine curve()  
      implicit none
      integer                                                :: i     = 0
      integer                                                :: j     = 0
      real(kind = x_precision)                               :: temp  = 0.0d0
      real(kind = x_precision), parameter                    :: t_min = 1.2d-2
      real(kind = x_precision), parameter                    :: t_max = 1.0d0
  