anchor.diff:    !-------------------------------------------------------------------------
     subroutine curve()  
       implicit none
  -    integer::i = 0
  +    integer                                                :: i     = 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
  -    integer,                  parameter                    :: nb_it = 2
  +    integer,                  parameter                    :: nb_it = 10
       real(kind = x_precision)                               :: sigma = 0.0d0
       real(kind = x_precision)                               :: Smin  = 0.0d0
       real(kind = x_precision)                               :: Smax  = 0.0d0
       real(kind = x_precision)                               :: eps   = 1.0d-4
       real(kind = x_precision)                               :: omega = 0.0d0
       real(kind = x_precision)                               :: r     = 5.0d0
  +    real(kind = x_precision)                               :: H     = 0.0d0    
  +    real(kind = x_precision)                               :: rho   = 0.0d0       
  +    real(kind = x_precision)                               :: cs    = 0.0d0
  +    real(kind = x_precision)                               :: nu    = 0.0d0  
  +    real(kind = x_precision)                               :: Q_plus = 0.0d0        
  +    real(kind = x_precision)                               :: Q_minus= 0.0d0
  +    real(kind = x_precision)                               :: K_ff   = 0.0d0    
  +    real(kind = x_precision)                               :: K_e    = 0.0d0      
  +    real(kind = x_precision)                               :: tau_eff= 0.0d0
  +    real(kind = x_precision)                               :: P_rad  = 0.0d0    
  +    real(kind = x_precision)                               :: P_gaz  = 0.0d0
  +    real(kind = x_precision)                               :: f      = 0.0d0
  +    real(kind = x_precision)                               :: E_ff   = 0.0d0
  +    real(kind = x_precision)                               :: Fz     = 0.0d0
  +    real(kind = x_precision)                               :: rs
  +    real(kind = x_precision)                               :: rmin
  +    real(kind = x_precision)                               :: Mdot_0
  +    real(kind = x_precision)                               :: Sigma_0
       type(parameters)                                       :: param
       !----------------------------------------------------------------------
       
       call get_parameters(param)
       call display_parameters()
  -    call initial_variables()
  +    call initial_variables(rs, rmin, Mdot_0, Sigma_0)
  +    call display_initial_variables(rs, rmin, Mdot_0, Sigma_0)
       !-------------------------------------------------------------------------
       ! Test for 1 value of r
       !-------------------------------------------------------------------------
  
positive.code: 
  
  
  real(kind=x_precision) function f(T, Sigma, Omega)
    implicit none
   
    type(parameters)::param
    real(kind = x_precision),intent(in)                      :: T,Sigma,Omega
    real(kind = x_precision)                                 :: coeff_a=0.,coeff_b=0.,coeff_c=0.
    real(kind = x_precision)                                 :: H     = 0.
    real(kind = x_precision)                                 :: rho   = 0.
    real(kind = x_precision)                                 :: cs    = 0.
    real(kind = x_precision)                                 :: nu    = 0.
    real(kind = x_precision)                                 :: Q_plus= 0.
    real(kind = x_precision)                                 :: K_ff  = 0.
    real(kind = x_precision)                                 :: K_e   = 0.
    real(kind = x_precision)                                 :: E_ff  = 0.
    real(kind = x_precision)                                 :: tau_eff=0.
    real(kind = x_precision)                                 :: Fz     =0.
    real(kind = x_precision)                                 :: Q_minus=0.
    real(kind = x_precision)                                 :: P_rad  =0.
    real(kind = x_precision)                                 :: P_gaz = 0.
    real(kind = x_precision)                                 :: Mdot_0 =0.
    real(kind = x_precision)                                 :: Sigma_0=0.
    real(kind = x_precision)                                 :: rmin  = 0.
    real(kind = x_precision)                                 :: rs    = 0.
    integer                                                  :: optical_depth = 0
    call get_parameters(param)
    
    rs                   = 2._x_precision * G * param%M/(c**2)
    rmin                 = 3._x_precision * rs
    Mdot_0               = param%Mdot 
    Sigma_0              = Mdot_0 /(state_0%Omega_0 * rs**2 * 2 * pi)
    coeff_a              = (Omega**2 * state_0%Omega_0**2 * Sigma * Sigma_0)/2.
    coeff_b              = (-1._x_precision/3._x_precision) * cst_rad*T**4 * state_0%T_0**4
    coeff_c              = (-1._x_precision * kmp * T * state_0%T_0 * Sigma * Sigma_0)/(2._x_precision * mu)
    call quadratic(coeff_a,coeff_b,coeff_c,H)
    rho                  = Sigma/H
    cs                   = Omega * H
    nu                   = param%alpha * cs * H
  
negative.code:        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)
  