anchor.diff:    subroutine compute_variables(state_out)
       implicit none
   
  -    real(kind = x_precision), dimension(n_cell) :: a, b, c
  +    real(kind = x_precision), dimension(n_cell) :: coeff_a, coeff_b, coeff_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
  
positive.code: 
    private
  
    public           :: compute_variables, dim_adim
  
  contains
  
  
negative.code:     type(state), intent(inout)                  :: state_out
  
      ! Compute trinomial coefficients for H
      coeff_a = (x_state%Omega * state_0%Omega_0)**2 * (state_out%S * state_0%S_0)
      coeff_b = - 2._x_precision * cst_rad * (state_out%T * state_0%T_0)**4 * x_state%x / (3._x_precision * state_0%H_0)
      coeff_c = - params%RTM * state_out%T * state_out%S * state_0%S_0 / state_0%H_0**2
      Delta = coeff_b**2 - 4._x_precision * coeff_a * coeff_c
  
      ! Compute variable depending on S, H
      state_out%H    = - 0.5_x_precision * (coeff_b + sign(sqrt(Delta), coeff_b)) / coeff_a
      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)
  