anchor.diff:      ! 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
  
positive.code:   ! FIXME I love it
    delta_T_max = 1e-5
  
    ! Initial time = 0
    t = 0._x_precision
  
    ! Read the parameters, generate state_0 and create adim state
    call get_parameters()
    call snapshot(s, 0, 0._x_precision, 13)
  
    ! Initiate the S_curve
    ! call s_curve(foo, bar)
  
negative.code: 
          ! Integrate S
          call do_timestep_S(s, dt_V)
  
          ! Recompute variables
          call compute_variables(s)
  
          ! Increment time, number of iterations
          t = t + dt_V
          iteration = iteration + 1
  
          ! 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'
  