anchor.diff:      dS_over_x_over_dx(n_cell)       = ( params%dx - nuS(n_cell) + nuS(n_cell-1) ) / &
                                         ( x_state%x(n_cell)**2 * s%v(n_cell) * params%dx**2 )
   
  -    ! Compute dS/dt using the corresponding equation, with d²(nuS)/dx² as (d(left)+d(right))/2
  -    dS_over_dt(1:n_cell-1) = ( nuS(2:n_cell) &
  -                             - 2._x_precision*nuS(1:n_cell-1) &
  -                             + nuS(1:n_cell-1) ) &
  -                             * overdx2 / x_state%x(1:n_cell-1)**2
  -    dS_over_dt(n_cell)     = ( params%dx - nuS(n_cell) + nuS(n_cell-1) ) / &
  -                             ( x_state%x(n_cell)**2 * params%dx**2 )
  -    ! FIX THIS : dS_over_dt is wrongly defined !
  -
  -
       !right term
       f_exp = ( 3._x_precision * state_0%v_0**2 * s%nu * x_state%Omega**2 &
               - s%Fz * x_state%x / s%s &
               + params%RTM * ( 4._x_precision - 3._x_precision * s%beta ) / s%beta * s%T / s%S &
  -            * ( dS_over_dt + s%v * dS_over_x_over_dx ) &
  +            * (dS_over_dt(s) + s%v * dS_over_x_over_dx ) &
               - s%Cv * s%v / x_state%x * dT_over_dx ) &
               / s%Cv
   
     end function f_exp
  -
   
     subroutine do_timestep_S_imp (s, dt)
     !process the temporal evolution of S
  
positive.code:   public :: do_timestep_T, do_timestep_S_imp, do_timestep_S_exp
  
  contains
  
    function f(s)
    !process the right term of \partial T* / \partial t* = (see Recapitulatif des adimensionenemts in report)
  
negative.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
  
      states%S = states%S + dt * dS_over_dt(states)
  
    end subroutine do_timestep_S_exp
  
  