anchor.diff:      ! Compute variables related to pressure
       state_out%Pgaz = state_out%rho * state_out%T
       state_out%Prad = state_out%T**4
  -    !state_out%beta = state_out%P_gaz / (state_out%P_gaz + state_0%beta_0 * state_out%P_rad)
       state_out%beta = 1._x_precision / (1._x_precision + state_0%beta_0 * state_out%Prad / state_out%Pgaz)
       state_out%Cv   = params%RTM * ((12._x_precision * (gammag - 1._x_precision) * &
                        (1._x_precision - state_out%beta)) + state_out%beta) / &
  
positive.code:     type(state), intent(inout)                  :: state_out
  
      ! Compute trinomial coefficients for H
      !a1 = r_state%Omega_r**2 * (state_out%S * state_0%S_0)
      !b1 = - 2._x_precision * cst_rad * (state_out%T * state_0%T_0)**4 * x_state%x / 3._x_precision
      !c1 = - (params%RTM * state_out%T * state_out%S * state_0%S_0)
      a1 = (x_state%Omega * state_0%Omega_0)**2 * (state_out%S * state_0%S_0)
      b1 = - 2._x_precision * cst_rad * (state_out%T * state_0%T_0)**4 * x_state%x / (3._x_precision * state_0%H_0)
      c1 = - params%RTM * state_out%T * state_out%S * state_0%S_0 / state_0%H_0**2
      Delta = b1**2 - 4._x_precision * a1 * c1
  
      ! Compute variable depending on S, H
      !state_out%H    = - 0.5_x_precision * (b1 + sign(sqrt(Delta),b1)) / a1 / state_0%H_0
      state_out%H    = - 0.5_x_precision * (b1 + sign(sqrt(Delta),b1)) / a1
      state_out%cs   = x_state%Omega * state_out%H
      state_out%nu   = params%alpha * state_out%cs * state_out%H
      state_out%rho  = state_out%S / (state_out%H * x_state%x)
  
negative.code:   subroutine compute_variables(state_out)
      implicit none
  
      real(kind = x_precision), dimension(n_cell) :: a, b, c
      real(kind = x_precision), dimension(n_cell) :: Delta
      real(kind = x_precision), dimension(n_cell) :: kappa_ff, tau, epsilo
      real(kind=x_precision) :: truc
  