anchor.diff:  
           ! Integrate S
           call do_timestep_S(s, dt_V)
  -        ! Increment time
  -        t = t + dt_V
  -
  -        ! print *, j, log10(s%T(50)*state_0%T_0), log10(s%S(50)*state_0%S_0)
   
           ! Recompute variables
           call compute_variables(s)
   
  +        ! Increment time, number of iterations
  +        t = t + dt_V
           iteration = iteration + 1
  -        ! Iterate while T hasn't converged
  -        T_converged = .false.
  -        do while (.not. T_converged)
  -          ! Integrate T
  -          ! print *, '# T integration'
  -          call do_timestep_T(s, dt_T, T_converged, delta_T_max)
  -          ! Increment time
  -          t = t + dt_T
  -
  -          ! Recompute variables
  -          call compute_variables(s)
  -          if (mod(iteration, output_freq) == 0) then
  -            call snapshot(s, iteration, t, 13)
  -            print*,'snapshot', iteration, t
  -          end if
  -
  -          ! print *, j, log10(s%T(50)*state_0%T_0), log10(s%S(50)*state_0%S_0), t
  -          iteration = iteration + 1
  -
  -        end do
   
           ! Output things here
           if (mod(iteration, output_freq) == 0) then
             call snapshot(s, iteration, t, 13)
             print*,'snapshot', iteration, t
           end if
  +
  +        ! Iterate while T hasn't converged
  +        T_converged = .false.
  +        do while (.not. T_converged)
  +          ! Integrate T
  +          call do_timestep_T(s, dt_T, T_converged, delta_T_max)
  +
  +          ! Recompute variables
  +          call compute_variables(s)
  +
  +          ! Increment time, number of iterations
  +          t = t + dt_T
  +          iteration = iteration + 1
  +
  +          if (mod(iteration, output_freq) == 0) then
  +            call snapshot(s, iteration, t, 13)
  +            print*,'snapshot', iteration, t
  +          end if
  +
  +        end do
  +
         end if
       else
         ! print*, 'Mdot kick'
  
positive.code:     ! state_0%M_dot_0 = state_0%M_dot_0 + 
    end do
  
    end do
         call system("rm "//"S_T.dat")
    close(13)
  
  end program black_hole_diffusion
  
negative.code:   ! FIXME I love it
    delta_T_max = 1e-5
  
    ! Read the parameters, generate state_0 and create adim state
    call get_parameters()
  
    ! Initiate the S_curve
    ! call s_curve(foo, bar)
  