anchor.diff:  
     ! Start iterations
     do while (iteration < n_iterations)
  -     
  -     ! Check that S is at a fixed point
  -     ! by computing the maximum difference between two consecutive S step
  -     if (maxval(abs((prev_S - s%S)/s%S)) > delta_S_max) then
  +     ! 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
   
  -        ! 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)
   
  -           ! Do a single S integration
  -           call do_timestep_S(s, dt_nu)
  +        ! Increment time, number of iterations
  +        t = t + dt_nu
  +        iteration = iteration + 1
  +
  +        ! Do a snapshot
  +        if (mod(iteration, output_freq) == 0) then
  +           call snapshot(s, iteration, t, 13)
  +           print*,'snapshot', iteration, t
  +        end if
  +
  +        ! Do T integrations
  +        T_converged = .false.
  +        do while (.not. T_converged)
  +           dt_tmp = dt_T
  +           ! Do a single T integration
  +           call do_timestep_T(s, dt_tmp, T_converged, delta_T_max)
   
              ! Increment time, number of iterations
  -           t = t + dt_nu
  +           t = t + dt_tmp
              iteration = iteration + 1
   
              ! Do a snapshot
  
positive.code:               call snapshot(s, iteration, t, 13)
                print*,'snapshot', iteration, t
             end if
  
             ! Do T integrations
             T_converged = .false.
             do while (.not. T_converged)
                ! Do a single T integration
                call do_timestep_T(s, dt_T, T_converged, delta_T_max)
  
                ! Increment time, number of iterations
                t = t + dt_T
                iteration = iteration + 1
  
                ! Do a snapshot
                if (mod(iteration, output_freq) == 0) then
                   call snapshot(s, iteration, t, 13)
                   print*,'snapshot', iteration, t
                end if
  
             end do
  
          end if
       else
          ! Mdot kick
          
          iteration = n_iterations
  
negative.code:   type(state)                                 :: s
    real(kind = x_precision)                    :: delta_S_max, delta_T_max, t
    real(kind = x_precision), dimension(n_cell) :: prev_S, S_crit
    real(kind = x_precision)                    :: t_nu, t_T
    real(kind = x_precision)                    :: dt_nu, dt_T, dt_tmp
    logical                                     :: T_converged
    real (kind = x_precision), dimension(n_cell):: dist, dist_crit
    
  
    ! FIXME
    S_crit = 1.e99_x_precision
  
    !----------------------------------------------
    ! Convergence criterion for S and T
    !----------------------------------------------
    delta_S_max = 1e-6
    delta_T_max = 1e-4
  
    ! Read the parameters, generate state_0 and create adim state
    call get_parameters()
  