anchor.diff:    ! Start
     iteration = 0
   
  -  ! Write initial state
  +  !----------------------------------------------
  +  ! Save initial snapshot
  +  !----------------------------------------------
     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
  +        !----------------------------------------------
  +        ! S integration
  +        !----------------------------------------------
           prev_S = s%S
   
           ! Do a single S integration
           call do_timestep_S(s, dt_nu)
   
  -        ! Increment time, number of iterations
  +        ! Update time, number of iterations
           t = t + dt_nu
           iteration = iteration + 1
   
  
positive.code:         end do
  
       end if
       if (maxval(abs((prev_S - s%S)/s%S)) < delta_S_max) then
          ! Mdot kick
          
          iteration = n_iterations
          !s%Mdot(n_cell) = s%Mdot(n_cell) * 1.01_x_precision
       end if
  
       ! Recompute variables and loop
  
negative.code:            print*,'snapshot', iteration, t
          end if
  
          !----------------------------------------------
          ! T integrations
          ! iterate while t hasn't converged
          !----------------------------------------------
          T_converged = .false.
          do while (.not. T_converged)
             dt_tmp = dt_T
  