anchor.diff:    integer                                     :: iteration, ios, j, i
     type(state)                                 :: s
     real(kind = x_precision)                    :: delta_S_max, delta_T_max, t
  -  real(kind = x_precision), dimension(n_cell) :: prev_S, prev_T, S_crit
  +  real(kind = x_precision), dimension(n_cell) :: prev_S, S_crit
     real(kind = x_precision)                    :: t_V, t_T
     real(kind = x_precision)                    :: dt_V, dt_T
  -  real(kind=x_precision)                      :: epst=1.e-3
  +  real(kind=x_precision)                      :: epst=1.e-5_x_precision
     logical                                     :: T_converged
   
     ! FIXME
  
positive.code:      stop "Error while opening output file."
    end if
  
    ! Initialize prev_S and prev_T
    ! we multiply delta to prevent the code from thinking it converged
    prev_S = 1.2*s%S
    prev_T = 1.2*s%T
    
    ! FIXME
    ! write(*,*) state_0%temps_0
  
negative.code:            j = 0
             T_converged = .false.
             do while (.not. T_converged)
                ! Integrate T
                ! print *, '# T integration'
                call do_timestep_T(s, dt_T, T_converged,epst)
  