anchor.diff:    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), dimension(n_cell) :: dt_nu, dt_T
  +  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
  -  ! FIXME
  -  delta_S_max = 1e-4
  -  ! FIXME I love it
  -  delta_T_max = 1e-3
  +
  +  !----------------------------------------------
  +  ! 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()
  
positive.code:   do i= 1, n_cell
      write(15, *)r_state%r(i), x_state%x(i), s%T(i), s%S(i), s%T(i)*state_0%T_0, s%S(i) * state_0%s_0 / x_state%x(i), CI%H_over_r(i)
    enddo
  
    close(15)
  
    call compute_variables(s)
    call timestep (s,dt_T, dt_nu)
  
    ! Open unit 13 to write
    ! do it more safely
    open(13, file="output.dat", status="replace", iostat=ios)
    if (ios /= 0) then
      stop "Error while opening output file."
    end if
  
    ! Initialize prev_S
    ! we multiply delta to prevent the code from thinking it converged
    prev_S = 1.2*s%S
  
   ! t_T  = params%t_T  !* state_0%temps_0
   ! t_nu = params%t_nu !*  state_0%temps_0
  
   ! dt_nu = t_nu / cst_dt
   ! dt_T  = t_T / cst_dt
  
    write(*,*) 'dt_T_min, dt_nu_min:', minval(dt_T), minval(dt_nu)
    write(*,*) 'dt_T_max, dt_nu_max:', maxval(dt_T), maxval(dt_nu)
  
    ! Initial time = 0
    t = 0._x_precision
  
  
negative.code: 
    ! 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
          prev_S = s%S
  
          ! 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_tmp
             iteration = iteration + 1
  
             ! Do a snapshot
  