anchor.diff:      rmin               = 3._x_precision * rs
       Mdot_0             = param%Mdot 
       Sigma_0            = Mdot_0 /(state_0%Omega_0 * rs**2 * 2 * pi)
  -    Omega_0            = sqrt( G*param%M/rs**3 )
  +    Omega_0            = state_0%Omega_0
   
     end subroutine initial_variables
   
  
positive.code:     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
  
negative.code:     integer                                                :: j     = 0
      real(kind = x_precision)                               :: temp  = 0.0d0
      real(kind = x_precision), parameter                    :: t_min = 8.5d-3
      real(kind = x_precision), parameter                    :: t_max = 2.3d0
      integer,                  parameter                    :: nb_it = 1000
      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-9
      real(kind = x_precision)                               :: omega = 0.0d0
      real(kind = x_precision)                               :: r     = 0.0d0
      real(kind = x_precision)                               :: H     = 0.0d0    
      real(kind = x_precision)                               :: rho   = 0.0d0       
      real(kind = x_precision)                               :: cs    = 0.0d0
  