anchor.diff:                call snapshot(s, iteration, t, 13)
                 print*,'snapshot', iteration, t
              end if
  +        end do
   
  -           ! Do T integrations
  -           T_converged = .false.
  -           do while (.not. T_converged)
  -              ! Do a single T integration
  -              call do_timestep_T(s, dt_T, T_converged, delta_T_max)
  -
  -              ! Increment time, number of iterations
  -              t = t + dt_T
  -              iteration = iteration + 1
  -
  -              ! Do a snapshot
  -              if (mod(iteration, output_freq) == 0) then
  -                 call snapshot(s, iteration, t, 13)
  -                 print*,'snapshot', iteration, t
  -              end if
  -
  -           end do
  -
  -        end if
  -     else
  +     end if
  +     if (maxval(abs((prev_S - s%S)/s%S)) < delta_S_max) then
           ! Mdot kick
           
           iteration = n_iterations
  
positive.code:   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_nu, t_T
    real(kind = x_precision), dimension(n_cell) :: dt_nu, dt_T
    logical                                     :: T_converged
    real (kind = x_precision), dimension(n_cell):: dist, dist_crit
    
  
    ! FIXME
    S_crit = 1.e99_x_precision
    ! FIXME
    delta_S_max = 1e-4
    ! FIXME I love it
    delta_T_max = 1e-3
  
    ! Read the parameters, generate state_0 and create adim state
    call get_parameters()
  
negative.code: 
    ! Start iterations
    do while (iteration < n_iterations)
       ! Check that we are not close to S_critical
       if (maxval(s%S - S_crit) > 0) then
          print*, 'Exiting because of S'
          ! Switch to explicit scheme
       else
          ! Implicit schem
          prev_S = s%S
  
          ! Do a single S integration
          call do_timestep_S(s, dt_nu)
  
          ! Increment time, number of iterations
          t = t + dt_nu
          iteration = iteration + 1
  
          ! Do a snapshot
          if (mod(iteration, output_freq) == 0) then
             call snapshot(s, iteration, t, 13)
             print*,'snapshot', iteration, t
          end if
  
          ! Do T integrations
          T_converged = .false.
          do while (.not. T_converged)
             dt_tmp = dt_T
             ! Do a single T integration
             call do_timestep_T(s, dt_tmp, T_converged, delta_T_max)
  
             ! Increment time, number of iterations
             t = t + dt_tmp
             iteration = iteration + 1
  
             ! Do a snapshot
  