anchor.diff:        kappa_ff(i) = 6.13e22_x_precision * state_0%rho_0 * state_out%rho(i) * &
                       (state_0%T_0 * state_out%T(i))**(-7._x_precision/2._x_precision)
   
  -      tau(i)      = 0.5_x_precision * sqrt(0.2_x_precision * (1._x_precision * &
  -                    para%X) * 6.13e22_x_precision * state_0%rho_0 * state_out%rho(i) * &
  -                    (state_0%T_0 * state_out%T(i))**(-7._x_precision/2._x_precision)) * &
  -                    state_0%S_0 * state_out%S(i) / state_out%x(i)
  +     ! tau(i)      = 0.5_x_precision * sqrt(0.2_x_precision * (1._x_precision * &
  +     !               para%X) * 6.13e22_x_precision * state_0%rho_0 * state_out%rho(i) * &
  +     !               (state_0%T_0 * state_out%T(i))**(-7._x_precision/2._x_precision)) * &
  +     !               state_0%S_0 * state_out%S(i) / state_out%x(i)
  +      tau(i)      = 0.5_x_precision * sqrt(kappa_e * kappa_ff(i)) * (state_out%S(i)/state_out%x(i) * state_0%S_0)
   
         epsil(i)    = 6.22e20_x_precision * (state_0%rho_0 * state_out%rho(i))**2 * &
                       sqrt((state_0%T_0 * state_out%T(i)))
  
positive.code:     real (kind = x_precision),              dimension(n_cell) :: kappa_ff, tau, epsil
      real (kind = x_precision)                                 :: kappa_e !cgs
  
      type(state),               intent(inout)                    :: state_out
  
      call get_parameters(para)
      kappa_e = 0.2_x_precision * (1._x_precision + para%X)
  
negative.code:     Delta = b1**2 - (4._x_precision * a1 * c1)
      !--------------------------------------------------------
      do i=1,n_cell
        state_out%H(i)     = - 0.5_x_precision * (b1(i) + sign(sqrt(Delta(i)),b1(i))) / a1(i) / state_0%H_0
        state_out%P_rad(i) = state_out%T(i)**4
        state_out%cs(i)    = state_out%Omega(i) * state_out%H(i)
        state_out%rho(i)   = state_out%S(i) / (state_out%H(i) * state_out%x(i))
        state_out%nu(i)    = para%alpha * state_out%cs(i) * state_out%H(i)
        state_out%P_gaz(i) = state_out%rho(i) * state_out%T(i)
        state_out%beta(i)  = state_out%P_gaz(i) / (state_out%P_gaz(i) + state_out%P_rad(i))
        state_out%Cv(i)    = para%RTM * ((12._x_precision * gammag -1._x_precision ) * &
                             (1._x_precision - state_out%beta(i)) + state_out%beta(i)) / &
  