anchor.diff:  
     ! Conversion into S*_crit
     S_c = sigma_c * x_state%x
  -  dist_crit = 500. * x_state%x
  +  dist_crit = 400. * x_state%x
   
     ! Initiate the S_curve
     ! call s_curve(foo, bar)
  
positive.code:      
       ! 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 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
          end if
          iteration = iteration + 1
  
  
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*
  