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: 
  end subroutine display_parameters
  
  subroutine initial_variables()
    write(*,*)'           Initial Variables            '
    write(*,*)'****************************************'
    write(*,"(' Temp_0      =',1p,E12.4)") state_0%T_0
    write(*,"(' Simga_0     =',1p,E12.4)") state_0%S_0/state_0%x_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)") state_0%M_dot_0
    write(*,"(' rho_0       =',1p,E12.4)") state_0%rho_0
  
    write(*,*)'****************************************'
    read(*,*)
  end subroutine initial_variables
  
  
  
  subroutine display_variables(temp,Omega,r,sigma)
    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)                     ::sol
  
    write(*,*)'               Variables                '
    write(*,*)'****************************************'
  
negative.code:   Q_minus             = 2._x_precision * Fz /Sigma
    Q_plus              = 9._x_precision /4._x_precision * nu * Omega**2
    f                   = Q_plus - Q_minus
    end subroutine variables
  
  
    !-------------------------------------------------------------------------
  