anchor.diff:  
       ! 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
  +    f0 = f(s)
  +    f1 = (f(s_deriv) - f(s)) / dtemp
   
       rhs = dt*f0 * (1._x_precision + dt*f1)
       
  
positive.code:                              * 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 )
      !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 &
  
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
  