anchor.diff:  
     subroutine do_timestep_S_exp (states, dt)
       !process the temporal evolution of S with an explicit algortyhm
  -    !using only on the superior branch of the S curve
  +    !using only on the superior branch of the S curve NOT VERIFIED YET
       implicit none
       
       type (state), intent(inout)                   :: states 
  
positive.code: 
      ! Let dT/dt = f0 and d²T/dt² = f1 so T(t+1) = T(t) + dt * f0 * (1 + dt * f1)
      
      f0 = f(s, dt)
      f1 = (f(s_deriv, dt) - f(s, dt)) / dtemp
  
      rhs = dt*f0 * (1._x_precision + dt*f1)
      
  
negative.code:     real (kind=x_precision)                    :: overdx, overdx2, overdt
      real (kind=x_precision), dimension(n_cell) :: dS_over_dt, dS_over_x_over_dx, S_over_x, dT_over_dx, nuS
  
      overdx  = 1._x_precision / params%dx
      overdx2 = overdx**2
      overdt  = 1._x_precision /dt
      nuS     = s%nu * s%S
  
      ! Compute dT/dx as the right spatial derivative
      dT_over_dx(1:n_cell - 1) = (s%T(2:n_cell) - s%T(1:n_cell - 1)) * overdx
      dT_over_dx(n_cell)       = 0
  
      ! Compute d(S/x)/dx as the right spatial derivative 
      S_over_x                        = s%S / x_state%x
  