anchor.diff:             prev_S = s%S
   
              ! Integrate S
  +           print *, 'S integration'
              call do_timestep_S(s, dt_V)
              ! Increment time
              t = t + dt_V
  
positive.code:   real(kind = x_precision)                    :: t_V, t_T  
    real(kind = x_precision)                    :: dt_V, dt_T
  
    t_V = params%t_nu / state_0%temps_0
    t_T = params%t_T / state_0%temps_0
    dt_V = t_V / 100._x_precision
    dt_T = t_T / 100._x_precision
  
    ! FIXME
    S_crit = 1.e99_x_precision
    ! FIXME
  
negative.code:            ! Recompute variables
             call compute_variables(s)
  
             ! Iterate while T hasn't converged
             j = 0
             do while (maxval(abs(prev_T - s%T)/s%T) > delta_T_max)
  