anchor.diff:      call quadratic(coeff_a , coeff_b , coeff_c , H)
   
   
  -
  +    
       rho                  = Sigma / H
       P_rad                = T**4
       P_gaz                = rho * T
       cs                   = Omega * H
  -    nu                   = params%alpha * cs * H
  +    nu                   = param%alpha * cs * H
       K_ff                 = 6.13d22 * rho_0 * rho * (T_0 * T)**(-3.5_x_precision)
  -    K_e                  = params%kappa_e
  +    K_e                  =  params%kappa_e
       E_ff                 = 6.22d20 * (rho_0 * rho)**2 * sqrt(T_0 * T)
       tau_eff              = 0.5_x_precision * sqrt(K_e * K_ff) * Sigma * Sigma_0
  -    nu_0                 = 2._x_precision * rs**2 * Omega_0/3._x_precision
  -    Fz_0                 = Sigma_0 * Omega_0 / 4._x_precision
  -
  -    !-------------------------------------------------------------------------
  -    !Select the case for the optical depth to compute Fz
  -    !-------------------------------------------------------------------------
  -
  -    optical_depth        = 1
  -
  -   ! if (tau_eff >= 1.)  then
  +    
  +    
  +    !-------------------------------------------------------------------------
  +    !Select the case for the opticaly depth to compute Fz
  +    !-------------------------------------------------------------------------
  +
  +    optical_depth        = 0
  +    
  +   ! if (tau_eff .ge. 1.)  then
      !    optical_depth     = 1
      ! else
      !    optical_depth     = 0
      ! end if
   
  -
  -
       select case(optical_depth)
  +
  +
       case(1)
   
          Fz = 4._x_precision * c**2 * T**4 /(27._x_precision * sqrt(3._x_precision) &
               * (K_ff + K_e) * Sigma * Sigma_0)
       case (0)
   
  -
  -       Fz = E_ff * rs * H / Fz_0
  -
  -
  +       Fz = 4._x_precision * rs * E_ff * H / (Omega_0 * Sigma_0)
  +
  +       
       end select
  -
  -    Q_plus_0 = Omega_0**3 * rs**2
  -
  -    Q_minus_0   = 2._x_precision * Fz_0 / sigma_0
  -
  -
  -
  +    
       Q_plus              = 3._x_precision  * rs**2 * nu * Omega**2 * Omega_0**2
  -    ! Q_plus =  nu * Omega**2
  -    !  Q_minus             = Fz  / Sigma * (1 / (3._x_precision * nu_0) * Omega_0 * rs**2)
       Q_minus             = Fz  / Sigma
   
       f                   = Q_plus - Q_minus
  -
  +    
     end subroutine variables
  -
  -
  +  
  +  
     !-------------------------------------------------------------------------
     ! Dichotomic function in order to determine the change of sign in a given
     ! interval [Smin,Smax] with an epsilon precision
  
positive.code:   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
  
      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)                     :: 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)                     :: tau_eff
      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(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
      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)
  
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)
      !-------------------------------------------------------------------------
  