subroutine doughnut_deriv ( t, y, dydt )

!*****************************************************************************80
!
!! doughnut_deriv() evaluates the right hand side of the doughnut ODE.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    10 October 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk8 ) T, Y(3): the time and solution value.
!
!  Output:
!
!    real ( kind = rk8 ) DYDT(3): the derivative value.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) dydt(3)
  real ( kind = rk8 ) m
  real ( kind = rk8 ) n
  real ( kind = rk8 ) t
  real ( kind = rk8 ) y(3)

  interface
    subroutine doughnut_parameters ( &
      m_in,  n_in,  t0_in,  y0_in,  tstop_in, &
      m_out, n_out, t0_out, y0_out, tstop_out )
      integer, parameter :: rk8 = kind ( 1.0D+00 )
      real ( kind = rk8 ), optional :: m_in
      real ( kind = rk8 ), optional :: m_out
      real ( kind = rk8 ), optional :: n_in
      real ( kind = rk8 ), optional :: n_out
      real ( kind = rk8 ), optional :: t0_in
      real ( kind = rk8 ), optional :: t0_out
      real ( kind = rk8 ), optional :: y0_in(3)
      real ( kind = rk8 ), optional :: y0_out(3)
      real ( kind = rk8 ), optional :: tstop_in
      real ( kind = rk8 ), optional :: tstop_out
    end subroutine
  end interface

  call r8_fake_use ( t )

  call doughnut_parameters ( m_out = m, n_out = n )

  dydt(1) = - m * y(2) + n * y(1) * y(3)
  dydt(2) =   m * y(1) + n * y(2) * y(3)
  dydt(3) =   0.5 * n * ( 1.0 - y(1)**2 - y(2)**2 + y(3)**2 )

  return
end
subroutine doughnut_parameters ( &
  m_in,  n_in,  t0_in,  y0_in,  tstop_in, &
  m_out, n_out, t0_out, y0_out, tstop_out )

!*****************************************************************************80
!
!! doughnut_parameters() returns parameters of the doughnut ODE.
!
!  Discussion:
!
!    If input values are specified, this resets the default parameters.
!    Otherwise, the output will be the current defaults.
!
!    Suggested choices for m, n, (y0) include:
!
!      3,  5, (1, 1, 3.0)
!      4,  5, (1, 1, 0.3)
!      pi, 5, (1, 1, 0.3)
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    10 October 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk8 ) m_in, n_in: parameters.
!
!    real ( kind = rk8 ) t0_in, y0_in: the initial time and value.
!
!    real ( kind = rk8 ) tstop_in: the final time.
!
!  Persistent:
!
!    real ( kind = rk8 ) m_default, n_default.
!
!    real ( kind = rk8 ) t0_default.
!
!    real ( kind = rk8 ) y0_default.
!
!    real ( kind = rk8 ) tstop_default.
!
!  Output:
!
!    real ( kind = rk8 ) m_out, n_out: parameters.
!
!    real ( kind = rk8 ) t0_out, y0_out: the initial time and value.
!
!    real ( kind = rk8 ) tstop_out: the final time.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ), save     :: m_default = 3.0D+00
  real ( kind = rk8 ), optional :: m_in
  real ( kind = rk8 ), optional :: m_out
  real ( kind = rk8 ), save     :: n_default = 5.0D+00
  real ( kind = rk8 ), optional :: n_in
  real ( kind = rk8 ), optional :: n_out
  real ( kind = rk8 ), save     :: t0_default = 0.0D+00
  real ( kind = rk8 ), optional :: t0_in
  real ( kind = rk8 ), optional :: t0_out
  real ( kind = rk8 ), save     :: y0_default(3) = (/ 1.0D+00, 1.0D+00, 3.0D+00 /)
  real ( kind = rk8 ), optional :: y0_in(3)
  real ( kind = rk8 ), optional :: y0_out(3)
  real ( kind = rk8 ), save     :: tstop_default = 10.0D+00
  real ( kind = rk8 ), optional :: tstop_in
  real ( kind = rk8 ), optional :: tstop_out
!
!  New values, if supplied on input, overwrite the current values.
!
  if ( present ( m_in ) ) then
    m_default = m_in
  end if

  if ( present ( n_in ) ) then
    n_default = n_in
  end if

  if ( present ( t0_in ) ) then
    t0_default = t0_in
  end if

  if ( present ( y0_in ) ) then
    y0_default = y0_in
  end if

  if ( present ( tstop_in ) ) then
    tstop_default = tstop_in
  end if
!
!  The current values are copied to the output.
!
  if ( present ( m_out ) ) then
    m_out = m_default
  end if

  if ( present ( n_out ) ) then
    n_out = n_default
  end if

  if ( present ( t0_out ) ) then
    t0_out = t0_default
  end if

  if ( present ( y0_out ) ) then
    y0_out = y0_default
  end if

  if ( present ( tstop_out ) ) then
    tstop_out = tstop_default
  end if

  return
end
subroutine r8_fake_use ( x )

!*****************************************************************************80
!
!! r8_fake_use() pretends to use an R8 variable.
!
!  Discussion:
!
!    Some compilers will issue a warning if a variable is unused.
!    Sometimes there's a good reason to include a variable in a program,
!    but not to use it.  Calling this function with that variable as
!    the argument will shut the compiler up.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    21 April 2020
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk8 ) X, the variable to be "used".
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) x

  if ( x /= x ) then
    write ( *, '(a)' ) '  r8_fake_use(): variable is NAN.'
  end if

  return
end