anchor.diff:      real (kind=x_precision), dimension(n_cell) :: f_exp
   
       real (kind=x_precision)                    :: overdx, overdx2
  -    real (kind=x_precision), dimension(n_cell) :: dS_over_dt, dS_over_x_over_dx, S_over_x, dT_over_dx, nuS
  +    real (kind=x_precision), dimension(n_cell) :: dS_over_x_over_dx, S_over_x, dT_over_dx, nuS
   
       overdx  = 1._x_precision / params%dx
       overdx2 = overdx**2
  
positive.code:     !using only on the superior branch of the S curve NOT VERIFIED YET
      implicit none
  
      type (state), intent(inout)                :: states
      real (kind=x_precision), intent(in)        :: dt
  
      real (kind=x_precision), dimension(0:n_cell+1) :: nuS ! two more for beginning / end
  
      nuS(0) = 0
      nuS(1:n_cell) = states%nu * states%S
      nuS(n_cell+1) = params%dx + nuS(n_cell)
  
      states%S = states%S + dt * (1._x_precision / x_state%x**2 * &
                 (nuS(2:n_cell+1) - 2._x_precision * nuS(1:n_cell) + nuS(0:n_cell-1)) / &
                 params%dx**2)
  
    end subroutine do_timestep_S_exp
  
  
negative.code:   public :: do_timestep_T, do_timestep_S_imp, do_timestep_S_exp
  
  contains
  
    function dS_over_dt(s)
      implicit none
  
      type(state), intent(in)                    :: s
      real(kind=x_precision), dimension(n_cell) :: dS_over_dt
      real(kind=x_precision), dimension(0:n_cell+1) :: nuS ! two more for beginning / end
  
      nuS(0) = 0
      nuS(1:n_cell) = s%nu * s%S
      nuS(n_cell+1) = params%dx + nuS(n_cell)
  
      dS_over_dt = 1._x_precision / x_state%x**2 * &
                   (nuS(2:n_cell+1) - 2._x_precision * nuS(1:n_cell) + nuS(0:n_cell-1)) / &
                   params%dx**2
  
    end function dS_over_dt
  
    function f(s)
    !process the right term of \partial T* / \partial t* = (see Recapitulatif des adimensionenemts in report)
  