anchor.diff:             prev_S = s%S
   
              ! Integrate S
  -           print *, '#S integration'
  -           call snapshot(s, iteration, t, 13)
              call do_timestep_S(s, dt_V)
              ! Increment time
              t = t + dt_V
  
positive.code:            do while (.not. T_converged)
                ! Integrate T
                ! print *, '# T integration'
                call do_timestep_T(s, dt_T, T_converged, epst)
                ! Increment time
                t = t + dt_T
                
                ! Recompute variables
                call compute_variables(s)
                if (mod(j, 1000) == 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
          end if
       else
         ! print*, 'Mdot kick'
  
negative.code: 
    iteration = 0
    ! Start iterating
    do while (iteration < n_iterations)
       ! Check that S is at a fixed point
       if (maxval(abs((prev_S - s%S)/s%S)) > delta_S_max) then
          ! Check here that S < S_crit
  