anchor.diff:                prev_T = s%T
   
                 ! Integrate T
  -              call do_timestep_T(s)
  +              call do_timestep_T(s, dt_T)
  +              ! Increment time
  +              t = t + dt_V
                 ! 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
        
  
positive.code:            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)
             
             ! Recompute variables
             call compute_variables(s)
  
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
  