anchor.diff:    subroutine variables(T, 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, Omega_0,rs, T_0, rho_0)
       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),intent(out)                     :: H
  -    real(kind = x_precision),intent(out)                     :: rho
  -    real(kind = x_precision),intent(out)                     :: cs
  -    real(kind = x_precision),intent(out)                     :: nu
  +    real(kind = x_precision),intent(out)                     :: H     
  +    real(kind = x_precision),intent(out)                     :: rho 
  +    real(kind = x_precision),intent(out)                     :: cs    
  +    real(kind = x_precision),intent(out)                     :: nu  
       real(kind = x_precision),intent(out)                     :: Q_plus
  -    real(kind = x_precision),intent(out)                     :: K_ff
  -    real(kind = x_precision),intent(out)                     :: K_e
  -    real(kind = x_precision),intent(out)                     :: E_ff
  +    real(kind = x_precision),intent(out)                     :: K_ff 
  +    real(kind = x_precision),intent(out)                     :: K_e  
  +    real(kind = x_precision),intent(out)                     :: E_ff  
       real(kind = x_precision),intent(out)                     :: tau_eff
  -    real(kind = x_precision),intent(out)                     :: Fz
  +    real(kind = x_precision),intent(out)                     :: Fz  
       real(kind = x_precision),intent(out)                     :: Q_minus
  -    real(kind = x_precision),intent(out)                     :: P_rad
  -    real(kind = x_precision),intent(out)                     :: P_gaz
  +    real(kind = x_precision),intent(out)                     :: P_rad 
  +    real(kind = x_precision),intent(out)                     :: P_gaz 
       real(kind = x_precision),intent(in)                      :: Sigma_0
       real(kind = x_precision),intent(in)                      :: Omega_0
       real(kind = x_precision),intent(in)                      :: T_0
       real(kind = x_precision),intent(in)                      :: rho_0
       real(kind = x_precision),intent(in)                      :: rs
  -    real(kind = x_precision)                                 :: nu_0
  -    real(kind = x_precision)                                 :: Fz_0
  -    real(kind = x_precision)                                 :: K_ff_0
  -    real(kind = x_precision)                                 :: Q_plus_0
  -    real(kind = x_precision)                                 :: Q_minus_0
  -
  -    real(kind = x_precision),intent(out)                     :: f
  +
  +    real(kind = x_precision),intent(out)                     :: f 
       integer                                                  :: optical_depth = 0
       !------------------------------------------------------------------------
  -
  -
  +   
       coeff_a              = (Omega**2 * Omega_0**2 * Sigma * Sigma_0)/2._x_precision
       coeff_b              = (-1._x_precision/3._x_precision) * cst_rad*T**4 * T_0**4 / rs
       coeff_c              = (-1._x_precision * params%RTM * T  *  Sigma * Sigma_0)/(2._x_precision * rs**2)
  
positive.code:     use mod_constants
      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_center
      real(kind=x_precision)                                   :: E_ff
  
negative.code:     real(kind = x_precision)                               :: Omega_0
      real(kind = x_precision)                               :: rho_0
      real(kind = x_precision)                               :: T_0
      real(kind = x_precision)                               :: temp_reel=0.0d0
      real(kind = x_precision)                               :: sigma_reel=0.0d0
      character(len = 8)                                     :: number_of_cell
      character(len = 64)                                    :: fname
      integer                                                :: fid
      !------------------------------------------------------------------------
      call display_parameters()
      call initial_variables(rs, rmin, Mdot_0, Sigma_0, Omega_0, T_0, rho_0)
      call display_initial_variables(rs, rmin, Mdot_0, Sigma_0, Omega_0, T_0, rho_0)
      !-------------------------------------------------------------------------
  