anchor.diff:             prev_S = s%S
   
              ! Integrate S
  -           print *, '#S integration'
  -           call snapshot(s, iteration, t, 13)
              call do_timestep_S(s, dt_V)
              ! Increment time
              t = t + dt_V
  
positive.code: 
    implicit none
  
    integer                                     :: iteration, ios, j, i
    type(state)                                 :: s
    real(kind = x_precision)                    :: delta_S_max, delta_T_max, t
    real(kind = x_precision), dimension(n_cell) :: prev_S, S_crit
    real(kind = x_precision)                    :: t_V, t_T
    real(kind = x_precision)                    :: dt_V, dt_T
    real(kind=x_precision)                      :: epst = 1.e-5_x_precision
    logical                                     :: T_converged
  
    ! FIXME
  
negative.code: 
    iteration = 0
    ! Start iterating
    do while (iteration < n_iterations)
       ! Check that S is at a fixed point
       if (maxval(abs((prev_S - s%S)/s%S)) > delta_S_max) then
          ! Check here that S < S_crit
  