anchor.diff:      real(kind = x_precision)                      :: dt, overdx, overdx2, overdt, dx
       integer                                       :: info
   
  -    real(kind = x_precision), dimension(n_cell)   :: diag
  +    real(kind = x_precision), dimension(n_cell)   :: dtemp, diag
       real(kind = x_precision), dimension(n_cell-1) :: diag_low, diag_up
       
       ! Get the timestep and the spacestep
  
positive.code:     dt=1d-6
  
      dtemp = 0.01_x_precision * states%T
      T_deriv = states%T + dtemp
      !call compute_variables(states%x, states%Omega, T_deriv, states%S, states_deriv)
      !modification of compute_varaible--just enter a state now -- FIXME
  
      f0 = f(states, dt, dx)
      f1 = (f(states_deriv, dt, dx) - f(states, dt, dx)) / dtemp
      
      ! Solve for T
      states%T = states%T + f0 / (1_x_precision / dt - f1)
  
      ! Copy the result
    end subroutine do_timestep_T
  
negative.code:   public :: do_timestep_T, do_timestep_S
    
  contains
    ! (T(t+dt) - T(t))/dt = f(s, dt, dx)
    function f (s, dt, dx)
  !process the right term of \partial T* / \partial t* = (see Recapitulatif des adimensionenemts in report)
      use mod_constants
  