anchor.diff:             print*,'snapshot', iteration, t
           end if
   
  -        ! Do T integrations
  +        !----------------------------------------------
  +        ! T integrations
  +        ! iterate while t hasn't converged
  +        !----------------------------------------------
           T_converged = .false.
           do while (.not. T_converged)
              dt_tmp = dt_T
  
positive.code:   ! Start
    iteration = 0
  
    ! Write initial state
    call snapshot(s, iteration, t, 13)
  
    ! 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
  
  
negative.code:         end do
  
       end if
  
       !----------------------------------------------
       ! Mdot kick
       ! increase Mdot at the boundaries if S is stalled
       !----------------------------------------------
       if (maxval(abs((prev_S - s%S)/s%S)) < delta_S_max) then
          ! Mdot kick
          ! params%kick_factor = params%kick_factor * 1.5
          print*, 'Mdot kick!', params%kick_factor
       end if
  
       ! Recompute variables and loop
  