anchor.diff:      real(kind=x_precision),intent(in)                        :: omega
       integer,intent(in)                                       :: optical_depth
   
  -    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)                                   :: f_min
       real(kind=x_precision)                                   :: f_max
       real(kind=x_precision)                                   :: f_center
  -    real(kind=x_precision)                                   :: E_ff
  -    real(kind=x_precision)                                   :: Fz
       real(kind=x_precision)                                   :: S_center = 0._x_precision
       integer                                                  :: j = 0
   
  
positive.code:   !Subroutine in order to compute variables H, rho, cs, nu, Q_plus, Q_minus,
    !K_ff, K_e, tau_eff, E_ff,Fz given T, Sigma and Omega
    !------------------------------------------------------------------------
    subroutine variables(T, Sigma, Omega, H, rho, cs, nu, Q_plus, Q_minus,&
         K_ff, K_e, tau_eff, E_ff,Fz,f, optical_depth)
      implicit none
  
      real(kind = x_precision),intent(in)                      :: T,Sigma,Omega
      integer, intent(in)                                      :: optical_depth
      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)                     :: f
      !------------------------------------------------------------------------
  
      coeff_a              = (Omega * state_0%Omega_0)**2 * Sigma * state_0%S_0 / 2._x_precision
  
negative.code:     real(kind = x_precision)                               :: eps   = 1.0d-4
      real(kind = x_precision)                               :: eps2   = 5.0d-1
  
      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)                               :: omega = 0.0d0
      real(kind = x_precision)                               :: r     = 0.0d0
      real(kind = x_precision)                               :: f     = 0.0d0
      integer                                                :: optical_depth =0
  
      real(kind = x_precision),dimension(nb_it)              :: temp_t_1  =0.0d0
  