anchor.diff:             prev_S = s%S
   
              ! Integrate S
  -           call do_timestep_S(s)
  -
  -           print *, "After S timestep", maxval(abs(prev_S - s%S)/s%S)
  -           call snapshot(s, iteration, t, 13)
  +           call do_timestep_S(s, dt_V)
  +           ! Increment time
  +           t = t + dt_V
              
              ! Recompute variables
              call compute_variables(s)
  
positive.code:               prev_T = s%T
  
                ! Integrate T
                call do_timestep_T(s)
                ! Recompute variables
                call compute_variables(s)
                print *, "After integration of T", j, maxval(abs(prev_T - s%T)/s%T)
                call snapshot(s, iteration, t, 13)
                j = j+1
  
                ! FIXME : increment time
             end do
             ! Output things here
             write(13, *) "YOLO DUDE"
             call snapshot(s, iteration, t, 13)
          end if
       end do
       
  
negative.code:   end if
  
    ! Initialize prev_S and prev_T
    ! we add initially 2*delta to prevent the code from thinking it converged
    prev_S = s%S + 2*delta_S_max
    prev_T = s%T + 2*delta_T_max
  
    call snapshot(s, iteration, t, 13)
    
    ! Start iterating
    do iteration = 1, 1!n_iterations
       ! Check that S is at a fixed point
       do while (maxval(abs((prev_S - s%S)/s%S)) > delta_S_max)
          print *, "While prev_S/s%S", iteration
  