subroutine ndbox ( func, ndim, norder, xtab, weight, result )
!
!***********************************************************************
!
!! NDBOX computes the integral of F(X) over an ND product box.
!
!
!  Discussion:
!
!    The routine creates an N-dimensional product rule from a 1D rule
!    supplied by the user.  The routine is fairly inflexible.  If
!    you supply a rule for integration from -1 to 1, then your product
!    box must be a product of NDIM copies of the interval [-1,1].
!
!  Modified:
!
!    30 October 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real, external FUNC, the name of a routine which evaluates
!    the function to be integrated, of the form:
!      function func ( ndim, x )
!      integer ndim
!      real func
!      real x(ndim)
!      func = ...
!      return
!      end
!
!    Input, integer NDIM, the dimension of the product region.
!
!    Input, integer NORDER, the number of points used in the 1D rule.
!
!    Input, real XTAB(NORDER), the abscissas of the 1D rule.
!
!    Input, real WEIGHT(NORDER), the weights of the 1D rule.
!
!    Output, real RESULT, the approximate value of the integral.
!
  implicit none
!
  integer ndim
  integer norder
!
  real, external :: func
  integer indx(ndim)
  integer k
  real result
  real w
  real weight(norder)
  real x(ndim)
  real xtab(norder)
!
  if ( ndim < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NDBOX - Fatal error!'
    write ( *, '(a)' ) '  NDIM < 1'
    write ( *, '(a,i6)' ) '  NDIM = ', ndim
    stop
  end if

  if ( norder < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NDBOX - Fatal error!'
    write ( *, '(a)' ) '  NORDER < 1'
    write ( *, '(a,i6)' ) '  NORDER = ', norder
    stop
  end if

  k = 0
  result = 0.0E+00

  do

    call tuple_next ( norder, ndim, k, indx )

    if ( k == 0  ) then
      exit
    end if

    w = product ( weight(indx(1:ndim)) )

    x(1:ndim) = xtab(indx(1:ndim))

    result = result + w * func ( ndim, x )

  end do

  return
end
subroutine ndp5 ( func, ndim, a, b, work, result )
!
!***********************************************************************
!
!! NDP5 approximates the integral of F(X) over a product region.
!
!
!  Discussion:
!
!    The routine uses a method which is exact for polynomials of total 
!    degree 5 or less.
!
!  Reference:
!
!    Philip Davis and Philip Rabinowitz,
!    Numerical Integration,
!    Blaisdell Publishing, 1967.
!
!  Modified:
!
!    30 October 2000
!
!  Parameters:
!
!    Input, real, external FUNC, the name of a routine which evaluates
!    the function to be integrated, of the form:
!      function func ( ndim, x )
!      integer ndim
!      real func
!      real x(ndim)
!      func = ...
!      return
!      end
!
!    Input, integer NDIM, the dimension of the region over
!    which integration is carried out.
!
!    Input, real A(NDIM), contains the lower limits of integration.
!
!    Input, real B(NDIM), contains the upper limits of integration.
!
!    Workspace, real WORK(NDIM).
!
!    Output, real RESULT, the approximate value of the integral.
!
  implicit none
!
  integer ndim
!
  real a(ndim)
  real a0
  real a1
  real a2
  real a3
  real a4
  real a5
  real b(ndim)
  real en
  real, external :: func
  integer i
  integer j
  real result
  real sum1
  real sum2
  real sum3
  real volume
  real work(ndim)
!
  if ( ndim < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NDP5 - Fatal error!'
    write ( *, '(a,i6)' ) '  NDIM < 1, NDIM = ', ndim
    stop
  end if
 
  a2 = 25.0E+00 / 324.0E+00
  a3 = sqrt ( 0.6E+00 )
  en = real ( ndim )
  a0 = ( 25.0E+00 * en * en - 115.0 * en + 162.0E+00 ) / 162.0E+00
  a1 = ( 70.0E+00 - 25.0E+00 * en ) / 162.0E+00
 
  volume = 1.0E+00
  do i = 1, ndim
    volume = volume * ( b(i) - a(i) )
    work(i) = 0.5E+00 * ( a(i) + b(i) )
  end do
 
  result = 0.0E+00
  if ( volume == 0.0E+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NDP5 - Warning!'
    write ( *, '(a)' ) '  Volume = 0, integral=0.'
    return
  end if
 
  sum1 = a0 * func(ndim,work)
  sum2 = 0.0E+00
  sum3 = 0.0E+00
 
  do i = 1, ndim
    work(i) = 0.5E+00 * ((a(i)+b(i))+a3*(b(i)-a(i)))
    sum2 = sum2 + func(ndim,work)
    work(i) = 0.5 * ((a(i)+b(i))-a3*(b(i)-a(i)))
    sum2 = sum2 + func(ndim,work)
    work(i) = 0.5 * (a(i)+b(i))
  end do
 
  if ( ndim <= 1 ) then
    go to 80
  end if
 
  a4 = a3
 
10    continue
 
  do i = 1, ndim-1
 
    work(i) = 0.5 * ((a(i)+b(i))+a4*(b(i)-a(i)))
    a5 = a3
 
50  continue
 
    do j = i+1, ndim
      work(j) = 0.5E+00 * ((a(j)+b(j))+a5*(b(j)-a(j)))
      sum3 = sum3 + func(ndim,work)
      work(j) = 0.5E+00 * (a(j)+b(j))
    end do
 
    a5 = -a5
    if ( a5 < 0.0E+00 ) then
      go to 50
    end if
 
    work(i) = 0.5 * ( a(i) + b(i) )
 
  end do
 
  a4 = -a4
  if (a4 < 0.0E+00 ) go to 10
 
80    continue
 
  result = volume * (sum1+a1*sum2+a2*sum3)
 
  return
end
subroutine ndromb ( func, a, b, ndim, nsub, maxit, eps, iwork, work, &
  nwork, result, ind )
!
!***********************************************************************
!
!! NDROMB approximates the integral of a F(X) over an ND product region.
!
!
!  Discussion:
!
!    The routine uses a Romberg method based on the midpoint rule.
!
!  Reference:
!
!    Philip Davis and Philip Rabinowitz,
!    Numerical Integration,
!    Blaisdell Publishing, 1967.
!
!  Modified:
!
!    30 October 2000
!
!  Parameters:
!
!    Input, real, external FUNC, the name of a routine which evaluates
!    the function to be integrated, of the form:
!      function func ( ndim, x )
!      integer ndim
!      real func
!      real x(ndim)
!      func = ...
!      return
!      end
!
!    Input, real A(NDIM), contains the lower limits of integration.
!
!    Input, real B(NDIM), contains the upper limits of integration.
!
!    Input, integer NDIM, the dimension of the region which is
!    to be integrated.
!
!    Input, integer NSUB(NDIM), the number of subintervals into
!    which the I-th integration interval (A(I), B(I)) is
!    initially subdivided.  NSUB(I) must be greater than 0.
!    The program will increase all of the NSUB(I)'s uniformly
!    as the iteration proceeds.
!
!    Input, integer MAXIT, the maximum number of iterations to
!    be performed.  The number of function evaluations on
!    iteration J is at least J**NDIM, which grows very rapidly.
!    MAXIT should be small!
!
!    Input, real EPS, an error tolerance for the approximation
!    of the integral.
!
!    Workspace, integer IWORK(MAXIT+NDIM).
!
!    Workspace, real WORK(MAXIT+NDIM).
!
!    Input, integer NWORK, the dimension of vectors IWORK and
!    WORK.  NWORK must be at least MAXIT+NDIM.
!
!    Output, real RESULT, the approximate value of the integral.
!
!    Output, integer IND, error return flag.
!    IND = -1 if the error tolerance could not be achieved.
!    IND = 1 if the error tolerance was achieved.
!
  implicit none
!
  integer ndim
  integer nwork
!
  real a(ndim)
  real b(ndim)
  real en
  real eps
  real factor
  real, external :: func
  integer i
  integer ind
  integer istep
  integer iwork(nwork)
  integer kdim
  integer ll
  integer maxit
  integer nsub(ndim)
  real resold
  real result
  real rnderr
  real submid
  real sum1
  real weight
  real work(nwork)
!
!  First NDIM entries of WORK contain point X
!  next MAXIT entries contain difference table
!
!  First NDIM entries of IWORK contain a pointer IPOINT
!  next MAXIT entries contain a counter of number of points, NPOINT
!
  if(ndim<1)then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NDROMB - Fatal error!'
    write ( *, '(a,i6)' ) '  NDIM is less than 1.  NDIM = ',ndim
    stop
  end if
 
  if ( maxit < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NDROMB - Fatal error!'
    write ( *, '(a,i6)' ) '  MAXIT is less than 1.  MAXIT = ',maxit
    stop
  end if
 
  do i = 1, ndim
    if ( nsub(i) <= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NDROMB - Fatal error!'
      write ( *, '(a)' ) '  NSUB(I) is less than 1.'
      write ( *, '(a,i6)' ) '  for I = ',i
      write ( *, '(a,i6)' ) '  NSUB(I)=',nsub(i)
      stop
    end if
  end do
 
  ind = 0
  rnderr = epsilon ( 1.0E+00 )
  iwork(ndim+1) = 1
 
  if ( maxit > 1 ) then
    iwork(ndim+2) = 2
  end if
 
  istep = 1
 
10    continue
 
  sum1 = 0.0E+00
 
  weight = 1.0E+00
  do i = 1, ndim
    weight = (b(i)-a(i))*weight / real(nsub(i))
    iwork(i) = 1
  end do
 
   50 continue
 
  kdim = 1
  do i = 1, ndim
    submid = real(iwork(i)) - 0.5E+00
    work(i) = a(i)+(b(i)-a(i))*submid/real(nsub(i))
  end do
  sum1 = sum1 + func(ndim,work)
!
!  See if more subintervals to treat in KDIM component
!
   70 continue
 
  if ( iwork(kdim) < nsub(kdim) ) then
    iwork(kdim) = iwork(kdim) + 1
    go to 50
  end if
!
!  Done with KDIM component.  Reset pointer and go to next component.
!
  iwork(kdim) = 1
  kdim = kdim+1

  if ( kdim <= ndim ) then
    go to 70
  end if
!
!  Done with summing.
!
  work(istep+ndim) = weight * sum1
 
  if ( istep <= 1 ) then
    result = work(ndim+1)
    resold = result
    ind = 1
    if(istep>=maxit)return
    istep = istep+1
    go to 110
  end if
!
!  Compute difference table for Richardson extrapolation
!
  en = real(iwork(ndim+istep))
 
  do ll = 2, istep
    i = istep+1-ll
    factor = real(iwork(ndim+i)**2-1)/en
    work(ndim+i) = work(ndim+i+1)+(work(ndim+i+1)-work(ndim+i))*factor
  end do
 
  result = work(ndim+1)
!
!  Error acceptance check
!
  ind = 1

  if ( abs(result-resold)<abs(result*(eps+rnderr)) ) then
    return
  end if

  ind = -1

  if ( istep >= maxit ) then
    return
  end if

  resold = result
  istep = istep+1
!
!  Determine new number of points.
!
  iwork(ndim+istep) = int(1.5*real(iwork(ndim+istep-1)))
 
  110 continue
 
  do i = 1, ndim
    nsub(i) = iwork(ndim+istep)*nsub(i)
  end do
 
  go to 10
end
subroutine ndsamp ( func, k1, k2, ndim, est1, err1, dev1, est2, &
  err2, dev2 )
!
!***********************************************************************
!
!! NDSAMP estimates the integral of F(X) over the N-dimensional unit box.
!
!
!  Discussion:
!
!    It computes two sequences of integral estimates, EST1 and EST2, for
!    indices K going from K1 to K2.  These estimates are produced by
!    the generation of 'random' abscissas in the region.  The process
!    can become very expensive if high accuracy is needed.
!
!    The total number of function evaluations is
!    4*(K1**NDIM+(K1+1)**NDIM+...+(K2-1)**NDIM+K2**NDIM), and K2
!    should be chosen so as to make this quantity reasonable.
!    In most situations, EST2(K) are much better estimates than EST1(K).
!
!  Reference:
!
!    Philip Davis and Philip Rabinowitz,
!    Numerical Integration,
!    Blaisdell Publishing, 1967.
!
!  Modified:
!
!    30 October 2000
!
!  Parameters:
!
!    Input, real, external FUNC, the name of a routine which evaluates
!    the function to be integrated, of the form:
!      function func ( ndim, x )
!      integer ndim
!      real func
!      real x(ndim)
!      func = ...
!      return
!      end
!
!    Input, integer K1, the beginning index for the iteration.
!    1 <= K1 <= K2.
!
!    Input, integer K2, the final index for the iteration.  K1 <= K2.
!    Increasing K2 increases the accuracy of the calculation,
!    but vastly increases the work and running time of the code.
!
!    Input, integer NDIM, the dimension of the region.  1 <= NDIM <= 10.
!
!    Output, real EST1(K2).  Entries K1 through K2 contain
!    successively better estimates of the integral.
!
!    Output, real ERR1(K2).  Entries K1 through K2 contain
!    the corresponding estimates of the integration errors.
!
!    Output, real DEV1(K2).  Entries K1 through K2 contain
!    estimates of the reliability of the the integration.
!    If consecutive values DEV1(K) and DEV1(K+1) do not differ
!    by more than 10 percent, then ERR1(K) can be taken as
!    a reliable upper bound on the difference between EST1(K)
!    and the true value of the integral.
!
!    Output, real EST2(K2).  Entries K2 through K2 contain
!    successively better estimates of the integral.
!
!    Output, real ERR2(K2).  Entries K2 through K2 contain
!    the corresponding estimates of the integration errors.
!
!    Output, real DEV2(K2).  Entries K2 through K2 contain
!    estimates of the reliability of the the integration.
!    If consecutive values DEV2(K) and DEV2(K+2) do not differ
!    by more than 10 percent, then ERR2(K) can be taken as
!    a reliable upper bound on the difference between EST2(K)
!    and the true value of the integral.
!
  implicit none
!
  integer k2
  integer, parameter :: maxdim = 10
  integer ndim
!
  real ak
  real ak1
  real akn
  real al(10)
  real b
  real be(ndim)
  real bk
  real d1
  real d2
  real dev1(k2)
  real dev2(k2)
  real dex(ndim)
  real err1(k2)
  real err2(k2)
  real est1(k2)
  real est2(k2)
  real, external :: func
  real g
  real ga(ndim)
  integer i
  integer j
  integer k
  integer k1
  integer key
  real p1(ndim)
  real p2(ndim)
  real p3(ndim)
  real p4(ndim)
  real s1
  real s2
  real t
  real y1
  real y2
  real y3
  real y4
!
  save al
!
  data al/ &
       0.4142135623730950, &
       0.7320508075688773, &
       0.2360679774997897, &
       0.6457513110645906, &
       0.3166247903553998, &
       0.6055512754639893, &
       0.1231056256176605, &
       0.3589989435406736, &
       0.7958315233127195, &
       0.3851648071345040/
!
!  Check input
!
  if ( ndim < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NDSAMP - Fatal error!'
    write ( *, '(a,i6)' ) '  NDIM must be at least 1, but NDIM = ', ndim
    stop
  end if
 
  if ( ndim > maxdim ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NDSAMP - Fatal error!'
    write ( *, '(a,i6)' ) '  NDIM must be no more than MAXDIM = ', maxdim
    write ( *, '(a,i6)' ) '  but NDIM = ', ndim
    stop
  end if
 
  if ( k1 < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NDSAMP - Fatal error!'
    write ( *, '(a,i6)' ) '  K1 must be at least 1, but K1 = ', k1
    stop
  end if
 
  if(k1>k2)then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NDSAMP - Fatal error!'
    write ( *, '(a)' ) '  K1 may not be greater than K2, but '
    write ( *, '(a,i6)' ) '  K1 = ',k1
    write ( *, '(a,i6)' ) '  K2 = ', k2
    stop
  end if
 
  do i = 1, ndim
    be(i) = al(i)
    ga(i) = al(i)
    dex(i) = 0.0E+00
  end do
 
  do k = k1, k2
 
    ak = real(k)
    key = 0
    ak1 = ak-1.1
    s1 = 0.0E+00
    d1 = 0.0E+00
    s2 = 0.0E+00
    d2 = 0.0E+00
    akn = ak**ndim
    t = sqrt(ak**ndim)*ak
    bk = 1.0E+00 / ak
 
20  continue
 
    key = key+1
    if (key==1) go to 50
    key = key-1
    j = 1
 
30  continue
 
    if(dex(j)<=ak1)then
      dex(j) = dex(j)+1.0E+00
      go to 50
    end if
 
    dex(j) = 0.0E+00
    j = j+1
    if (j<=ndim) go to 30
    go to 70
 
50  continue
 
    do i = 1, ndim
      b = be(i)+al(i)
      if ( b > 1.0 ) b = b-1.0E+00
      g = ga(i)+b
      if ( g > 1.0 ) g = g-1.0E+00
      be(i) = b+al(i)
      if ( be(i) > 1.0 ) be(i) = be(i)-1.0E+00
      ga(i) = be(i) + g
      if ( ga(i) > 1.0 ) ga(i) = ga(i)-1.0E+00
      p1(i) = ( dex(i) + g ) * bk
      p2(i) = ( dex(i) + 1.0 - g ) *bk
      p3(i) = ( dex(i) + ga(i) ) * bk
      p4(i) = ( dex(i) + 1.0 - ga(i) ) * bk
    end do
 
    y1 = func(ndim,p1)
!
!  There may be an error in the next two lines,
!  but oddly enough, that is how the original reads
!
    y3 = func(ndim,p2)
    y2 = func(ndim,p3)
    y4 = func(ndim,p4)
    s1 = s1 + y1 + y2
    d1 = d1 + ( y1 - y2 )**2
    s2 = s2 + y3 + y4
    d2 = d2 + ( y1 + y3 - y2 - y4 )**2
    go to 20
 
70  continue
 
    est1(k) = 0.5E+00 * s1 / akn
    err1(k) = 1.5E+00 * sqrt ( d1 ) / akn
    dev1(k) = err1(k) * t
    est2(k) = 0.25E+00 * ( s1 + s2 ) / akn
    err2(k) = 0.75E+00 * sqrt ( d2 ) / akn
    dev2(k) = err2(k) * t * ak
 
  end do
 
  return
end
subroutine ndsum2 ( func, xtab, weight, ntab, ndim, nmax, result )
!
!***********************************************************************
!
!! NDSUM2 computes the integral of F(X) over an N-dimensional product region.
!
!
!  Discussion:
!
!    The routine uses a product rule supplied by the user.
!
!    The region may be a product of any combination of finite,
!    semi-infinite, or infinite intervals.
!
!    For each factor in the region, it is assumed that an integration
!    rule is given, and hence, the region is defined implicitly by
!    the integration rule chosen.
!
!  Reference:
!
!    Philip Davis and Philip Rabinowitz,
!    Numerical Integration,
!    Blaisdell Publishing, 1967.
!
!  Modified:
!
!    30 October 2000
!
!  Parameters:
!
!    Input, real, external FUNC, the name of a routine which evaluates
!    the function to be integrated, of the form:
!      function func ( ndim, x )
!      integer ndim
!      real func
!      real x(ndim)
!      func = ...
!      return
!      end
!
!    Input, real XTAB(NMAX,NDIM).  XTAB(I,J) is the I-th abscissa
!    for the J-th rule.
!
!    Input, real WEIGHT(NMAX,NDIM).  WEIGHT(I,J) is the I-th
!    weight for the J-th rule.
!
!    Input, integer NTAB(NDIM).  NTAB(I) is the number of
!    abscissas to be used in the J-th rule.  NTAB(I) must be
!    greater than 0 and less than or equal to NMAX.
!
!    Input, integer NDIM, the dimension of the product region,
!    the second dimension of the vectors XTAB and WEIGHT, the
!    dimension of WORK, IWORK and NTAB, and the number of integration rules.
!
!    Input, integer NMAX, the declared first dimension of WEIGHT and XTAB
!    in the calling program.
!
!    Output, real RESULT, the approximate value of the integral.
!
  implicit none
!
  integer ndim
  integer nmax
!
  real, external :: func
  integer i
  integer iwork(ndim)
  integer k
  integer m1
  integer ntab(ndim)
  real result
  real w1
  real weight(nmax,ndim)
  real work(ndim)
  real xtab(nmax,ndim)
!
  if ( ndim < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NDSUM2 - Fatal error!'
    write ( *, '(a)' ) '  NDIM < 1'
    write ( *, '(a,i6)' ) '  NDIM = ', ndim
    stop
  end if
 
  do i = 1, ndim
    if ( ntab(i) < 1 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NDSUM2 - Fatal error!'
      write ( *, '(a)' ) '  NTAB(I) < 1.'
      write ( *, '(a,i6)' ) '  For I = ', i
      write ( *, '(a,i6)' ) '  NTAB(I) = ', ntab(i)
      stop
    end if

    if ( ntab(i) > nmax )then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NDSUM2 - Fatal error!'
      write ( *, '(a)' ) '  NTAB(I) > NMAX.'
      write ( *, '(a,i6)' ) '  For I = ', i
      write ( *, '(a,i6)' ) '  NTAB(I) = ', ntab(i)
      write ( *, '(a,i6)' ) '  NMAX = ', nmax
      stop
    end if
  end do
 
  iwork(1:ndim) = 1
 
  result = 0.0E+00
 
  do
 
    k = 1
 
    w1 = 1.0E+00
    do i = 1, ndim
      m1 = iwork(i)
      work(i) = xtab(m1,i)
      w1 = w1 * weight(m1,i)
    end do

    result = result + w1 * func(ndim,work)
 
    do while ( iwork(k) == ntab(k) )

      iwork(k) = 1
      k = k + 1
      if ( k > ndim ) then
        return
      end if

    end do
 
    iwork(k) = iwork(k) + 1

  end do

  return 
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end
subroutine tuple_next ( m, n, k, x )
!
!*******************************************************************************
!
!! TUPLE_NEXT computes the next element of a tuple space.
!
!
!  Discussion:
!
!    The elements are N vectors.  Each entry is constrained to lie
!    between 1 and M.  The elements are produced one at a time.
!    The first element is
!      (1,1,...,1)
!    and the last element is
!      (M,M,...,M)
!    Intermediate elements are produced in lexicographic order.
!
!  Examples:
!
!    N = 2, M = 3
!
!    K   X
!    1   1 1
!    2   1 2
!    3   1 3
!    4   2 1
!    5   2 2
!    6   2 3
!    7   3 1
!    8   3 2
!    9   3 3
!
!  Modified:
!
!    30 October 2000
!
!  Parameters:
!
!    Input, integer M, the maximum entry.
!
!    Input, integer N, the number of components.
!
!    Input/output, integer K, counts the elements.
!    On first call, set K to 0.  Thereafter, K will indicate the
!    order of the element returned.  When there are no more elements,
!    K will be returned as 0.
!
!    Input/output, integer X(N), on input the previous tuple.
!    On output, the next tuple.
!
  implicit none
!
  integer n
!
  integer i
  integer k
  integer m
  integer x(n)
!
  if ( m < 1 ) then
    return
  end if

  if ( k <= 0 ) then

    x(1:n) = 1
    k = 1

  else

    k = k + 1
    i = n

    do

      if ( x(i) < m ) then
        x(i) = x(i) + 1
        exit
      end if

      x(i) = 1

      if ( i == 1 ) then
        k = 0
        exit
      end if

      i = i - 1

    end do

  end if

  return
end