anchor.diff:    !process the temporal evolution of T
       implicit none
   
  -    type (state), intent(inout)       :: s
  -    logical, intent(out)              :: converge
  -    real(kind=x_precision),intent(in) :: epst
  +    real(x_precision), intent(in)    :: epst
  +    real(x_precision), intent(in)    :: dt
  +    type(state),       intent(inout) :: s
  +    logical,           intent(out)   :: converge
   
  -    type (state)                                :: s_deriv
  -    real(kind = x_precision), intent(in)        :: dt
   
  -    real(kind = x_precision), dimension(n_cell) :: dtemp, rhs, newT
  -    real(kind = x_precision), dimension(n_cell) :: f0, fT
  -    real(kind = x_precision) :: maxi
  +    real(x_precision), dimension(n_cell) :: dtemp, rhs, newT
  +    real(x_precision), dimension(n_cell) :: f0, fT
  +    real(x_precision)                    :: maxi
  +    type(state)                          :: s_deriv
   
       dtemp     = s%T * 1.e-3_x_precision
       s_deriv   = s
  
positive.code:     diag_low = -s%nu(1:n_cell-1) / x2(2:  n_cell) * dtoverdx2
  
      diag(n_cell)       = 1 + s%nu(n_cell) / x2(n_cell) * dtoverdx2
      diag_low(n_cell-1) = -s%nu(n_cell-1) / x2(n_cell)  * dtoverdx2
  
      s%S(n_cell) = s%S(n_cell) + dt * overdx / x2(n_cell)
  
  
negative.code:     end if
    end subroutine do_timestep_S_imp
  
    subroutine do_timestep_S_exp(states, dt)
      !process the temporal evolution of S with an explicit algorithm
      !using only on the superior branch of the S curve NOT VERIFIED YET
      implicit none
  
      real(x_precision), intent(in)    :: dt
      type(state),       intent(inout) :: states
  
      states%S = states%S + dt * dS_dt(states)
  
    end subroutine do_timestep_S_exp
  
  