anchor.diff:      !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
  +    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)
  +    states%S = states%S + dt * dS_over_dt(states)
   
     end subroutine do_timestep_S_exp
   
  
positive.code:     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 ) &
              - 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
  
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)
  