anchor.diff:    real(x_precision), dimension(n_cell) :: S_c
     !-------------------------------------------------------------------------
   
  -  integer                              :: iteration, ios, i, T_steps
  +  integer                              :: start_iteration, iteration, ios, i, T_steps, S_steps
     type(state)                          :: s
  -  real(x_precision), dimension(n_cell) :: prev_S, prev_T
  -  real(x_precision)                    :: delta_S_max, delta_T_max, S_rel_diff
  +  real(x_precision), dimension(n_cell) :: prev_S, prev_T, S_rel_diff
  +  real(x_precision)                    :: delta_S_max, delta_T_max, S_rel_diff_max
     real(x_precision)                    :: t
     real(x_precision)                    :: dt_nu, dt_T, pf
     integer :: p
  
positive.code:        unstable = .true.
       end if
  
    elseif (arg == 'start') then
       s%T = IC%T / state_0%T_0
       s%S = IC%Sigma / state_0%S_0 * x_state%x
       ! Initial time = 0
       t = 0._x_precision
       ! Start
       iteration = 0
       unstable = .false.
       s%Mdot(n_cell) = 1._x_precision
  
negative.code: 
          if (mod(iteration, output_freq) == 0) then
             call snapshot(s, iteration, t, 13)
             if ((arg == 'load' .or. arg == 'restart') &
                  .and. iteration-start_iteration > stp_value) then
                stop
             endif
             print*,'snapshot', iteration, t, 'exp', maxval(s%T - T_c), maxval(s%S - S_c), dt_T
  