subroutine kdv_ift ( nx, nt, x, tplt, uplt )

!*****************************************************************************80
!
!! kdv_ift() solves the Korteweg-deVries equation using the IFT method.
!
!  Discussion:
!
!    The system being solved is:
!
!    ut + u ux + uxxx = 0 on -pi < x < pi
!    by FFT with integrating factor v = exp(-ik^3t)*uhat.
!
!    This code is related to p27.m in the Trefethen reference.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 April 2020
!
!  Author:
!
!    Original MATLAB version by Lloyd Trefethen.
!    This version by John Burkardt.
!
!  Reference:
!
!    Aly-Khan Kassam, Lloyd Trefethen,
!    Fourth-order time-stepping for stiff ODE's,
!    SIAM Journal on Scientific Computing,
!    Volume 26, Number 4, pages 1214-1233, 2005.
!
!    Lloyd Trefethen,
!    Spectral methods in MATLAB,
!    SIAM, 2000,
!    LC: QA377.T65
!    ISBN: 978-0-898714-65-4
!
!  Input:
!
!    integer NX: the number of nodes.
!    NX should be even.
!
!    integer NT: the number of time points.
!
!  Output:
!
!    real ( kind = rk ) X(NX): the spatial grid.
!
!    real ( kind = rk ) TPLT(NT): the time values.
!
!    real ( kind = rk ) UPLT(NX,NT): solution values.
!
  implicit none

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

  integer nt
  integer nx

  complex ( kind = ck ) a(nx)
  complex ( kind = ck ) b(nx)
  complex ( kind = ck ) c(nx)
  real ( kind = rk ) c1
  real ( kind = rk ) c2
  complex ( kind = ck ) d(nx)
  real ( kind = rk ) dt
  complex ( kind = ck ) e(nx)
  complex ( kind = ck ) e2(nx)
  complex ( kind = ck ) g(nx)
  complex ( kind = ck ), parameter :: i = cmplx ( 0.0D+00, 1.0D+00, kind = ck )
  integer ier
  integer, parameter :: inc = 1
  integer j
  integer jplt
  integer jstep
  integer k(nx)
  complex ( kind = ck ) l(nx)
  integer lensav
  integer lenwrk
  integer nmax
  real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00
  real ( kind = rk ) t
  real ( kind = rk ) tmax
  real ( kind = rk ) tplt(nt)
  real ( kind = rk ) u(nx)
  real ( kind = rk ) uplt(nx,nt)
  complex ( kind = ck ) v(nx)
  complex ( kind = ck ) w(nx)
  real ( kind = rk ), allocatable :: work ( : )
  real ( kind = rk ), allocatable :: wsave ( : )
  real ( kind = rk ) x(nx)
  real ( kind = rk ) xhi
  real ( kind = rk ) xlo
!
!  Initialize FFT.
!
  lenwrk = 2 * nx
  allocate ( work(1:lenwrk) )
  lensav = 2 * nx + int ( log ( real ( nx ) ) ) + 4
  allocate ( wsave(1:lensav) )
  call zfft1i ( nx, wsave, lensav, ier )
  if ( ier /= 0 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'kdv_ift(): Fatal error!'
    write ( *, '(a,i2)' ) '  zfft1i() returns ier = ', ier
    stop ( 1 )
  end if
!
!  Set up grid.
!
  dt = 0.4D+00 / nx**2
  xlo = - r8_pi
  xhi = + r8_pi
  do j = 1, nx
    x(j) = ( real ( nx - j + 1, kind = rk ) * xlo   &
           + real (      j - 1, kind = rk ) * xhi ) &
           / real ( nx,         kind = rk )
  end do
!
!  Set up two-soliton initial data.
!
  c1 = 25.0
  c2 = 16.0

  u = 3.0 * c1**2 / ( cosh ( 0.5 * ( c1 * ( x + 2.0 ) ) ) )**2 &
    + 3.0 * c2**2 / ( cosh ( 0.5 * ( c2 * ( x + 1.0 ) ) ) )**2
!
!  Compute V = FFT(U).
!
  v = u
  call zfft1f ( nx, inc, v, nx, wsave, lensav, work, lenwrk, ier )
  if ( ier /= 0 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'kdv_ift(): Fatal error!'
    write ( *, '(a,i2)' ) '  zfft1f returns ier = ', ier
    stop ( 1 )
  end if
!
!  Set the time step.
!
  dt = 0.4D+00 / nx**2
!
!  Set the wave numbers.
!
  do j = 1, ( nx / 2 )
    k(j) = j - 1
  end do

  k(nx/2+1) = 0
  do j = nx/2+2, nx
    k(j) = j - nx - 1
  end do
!
!  Fourier multipliers.
!
  l = cmplx ( 0.0D+00, k**3, kind = ck )
  e = exp ( dt * l )
  e2 = exp ( dt * l / 2.0 )
!
!  Solve the PDE.
!
  tmax = 0.006D+00
  nmax = nint ( tmax / dt )
  jstep = nmax / ( nt - 1 )
  
  jplt = 1
  tplt(jplt) = 0.0D+00
  uplt(1:nx,jplt) = u(1:nx)

  g = - 0.5D+00 * i * dt * k

  do j = 1, nmax

    t = j * dt
!
!   a = g * fft ( real ( ifft ( v ) ) ** 2 )
!
    a = v
    call zfft1b ( nx, inc, a, nx, wsave, lensav, work, lenwrk, ier )
    a = real ( a )
    a = a**2
    call zfft1f ( nx, inc, a, nx, wsave, lensav, work, lenwrk, ier )
    a = g * a
!
!   b = g * fft ( real ( ifft ( e  * ( v + a / 2.0 ) ) ) ** 2 )
!
    b = e * ( v + 0.5D+00 * a )
    call zfft1b ( nx, inc, b, nx, wsave, lensav, work, lenwrk, ier )
    b = real ( b )
    b = b**2
    call zfft1f ( nx, inc, b, nx, wsave, lensav, work, lenwrk, ier )
    b = g * b
!
!   c = g * fft ( real ( ifft ( e  *   v + b / 2.0   ) ) ** 2 )
!
    c = e * v + 0.5D+00 * b
    call zfft1b ( nx, inc, c, nx, wsave, lensav, work, lenwrk, ier )
    c = real ( c )
    c = c**2
    call zfft1f ( nx, inc, c, nx, wsave, lensav, work, lenwrk, ier )
    c = g * c
!
!   d = g * fft ( real ( ifft ( e2 *   v + e * c     ) ) ** 2 )
!
    d = e2 * v + e * c
    call zfft1b ( nx, inc, d, nx, wsave, lensav, work, lenwrk, ier )
    d = real ( d )
    d = d**2
    call zfft1f ( nx, inc, d, nx, wsave, lensav, work, lenwrk, ier )
    d = g * d

    v = e2 * v + ( e2 * a + 2.0 * e * ( b + c ) + d ) / 6.0

    if ( mod ( j, jstep ) == 0 ) then
      w = v
      call zfft1b ( nx, inc, w, nx, wsave, lensav, work, lenwrk, ier )
      if ( ier /= 0 ) then
        write ( *, '(a)' ) ''
        write ( *, '(a)' ) 'kdv_ift(): Fatal error!'
        write ( *, '(a,i2)' ) '  zfft1b returns ier = ', ier
        stop ( 1 )
      end if
      u = real ( w )
      jplt = jplt + 1
      tplt(jplt) = t
      uplt(1:nx,jplt) = u(1:nx)
    end if

  end do
!
!  Free FFT memory.
!
  deallocate ( work )
  deallocate ( wsave )

  return
end