anchor.diff:       
        ! Check that we are not close to S_critical
        if (1. - maxval(s%S/S_c) < 1.e-2) then
  -        call do_timestep_S_exp(s, min_dt_T*8.e-3)
  -        call do_timestep_T(s, min_dt_T*1.e-2, T_converged, delta_T_max)
  -        t = t + 2*min_dt_T*1e-2
  +        call do_timestep_S_exp(s, maxval(dt_T*1.e-1))
  +        call do_timestep_T(s, maxval(dt_T*1.e-1), T_converged, delta_T_max)
  +        t = t + 2*maxval(dt_T*1e-1)
           call compute_variables(s)
           call timestep (s, dt_T, dt_nu)
   
           if (mod(iteration, output_freq) == 0) then ! .or. &
              !  (iteration > 790000 .and. mod(iteration, 10) == 0)) then
              call snapshot(s, iteration, t, 13)
  -           print*,'snapshot', iteration, t, 1 - maxval(s%S/S_c), dt_pre_factor
  +           print*,'snapshot', iteration, t, 1 - maxval(s%S/S_c), 1 - minval(s%S/s_c), dt_pre_factor
           end if
           iteration = iteration + 1
   
  
positive.code:   implicit none
  
    !--------------------------- Parameters for s_curve-----------------------
    real(x_precision), parameter         :: eps_in = 1.0e-7_x_precision ! Precision required for the dichotomy
    real(x_precision), parameter         :: Tmin   = 2.5e-2_x_precision / 2.
    real(x_precision), parameter         :: Tmax   = 4.49e0_x_precision * 2.
    real(x_precision), parameter         :: Smin   = 2.36e1_x_precision / 5.
  
negative.code:   !----------------------------------------------
    ! Set the initial conditions for S, T
    !----------------------------------------------
    s%T = 1.8_x_precision*CI%T_ci / state_0%T_0
    s%S = CI%Sig_ci / state_0%S_0 * x_state%x
  
    ! H_over_r becoming H*
  