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:        converge = .false.
      end if
  
      ! Copy the result
    end subroutine do_timestep_T
  
  end module mod_integrator
  
negative.code:     call compute_variables(s)
      ! Let dT/dt = f0 and d²T/dt² = f1 so T(t+1) = T(t) + dt * f0 * (1 + dt * f1)
  
      f0 = dT_dt_imp(s)
      fT = (dT_dt_imp(s_deriv) - f0) / dtemp
  
      rhs = f0 / fT * (exp(fT*dt) - 1._x_precision)
      ! rhs = dt*f0 * (1._x_precision + 0.5_x_precision*dt*fT)
  