program refine
!
!*******************************************************************************
!
!! REFINE extends a computed flow solution to a refined grid.
!
!
!  Discussion:
!
!    To use this program, you start with the following output files
!    from a run of FLOW6:
!
!      a NODE file,
!      an ELEMENT file,
!      a COEFFICIENT file.
!
!    You specify refinement factors in the horizontal and vertical
!    dimensions.  A horizontal factor of 3, for instance, means that
!    you want to use three times as many elements in the horizontal
!    direction.
!
!    REFINE then evaluates the solution defined in the input files
!    at the appropriate nodes, and constructs the coefficient vector
!    for the refined grid, as though the original solution had been
!    computed there.
!
!
!    REFINE is not a perfect program.  It uses an accidental fact
!    as though it were guaranteed to happen, namely, that the ordering
!    of the nodes in the refined grid can be determined by sorting the
!    nodes by their X and Y coordinates.  This is really only true if
!    the nodes line up nicely in vertical rows.  If there's any significant
!    skew in the layout of the nodes, this idea will fail.  However,
!    it is better than assuming that the nodes are laid out in a
!    complete rectangular array, which is definitely not true for some 
!    of the regions we study!  Someday I'll get back to this.  I need
!    to add some data structures to keep track of where elements are,
!    or else change the way I define node numbers.
!
!  Modified:
!
!    28 June 2000
!
!  Author:
!
!    John Burkardt
!
  implicit none
!
  integer, parameter :: maxnpe = 6
  integer, parameter :: maxnx = 160
  integer, parameter :: maxny = 60

  integer, parameter :: maxelm = 2 * maxnx * maxny
  integer, parameter :: maxeqn = 2 * ( 2 * maxnx + 1 ) * ( 2 * maxny + 1 ) &
    + ( maxnx + 1 ) * ( maxny + 1 )
  integer, parameter :: maxnp = ( 2 * maxnx + 1 ) * ( 2 * maxny + 1 )
!
  integer col_factor
  real dist
  real distmin
  character ( len = 80 ) file_name
  real g1(maxeqn)
  real g2(maxeqn)
  integer i
  integer ieqn2
  integer ierror
  integer ihi
  integer indx1(3,maxnp)
  integer indx2(3,maxnp)
  integer iperm(maxnp)
  integer ip
  integer j
  integer kind(maxnp)
  integer mp2
  integer nelem1
  integer neqn1
  integer neqn2
  integer node1(maxnpe,maxelm)
  integer np1
  integer np2
  integer npe
  integer nx1
  integer nx2
  integer ny1
  integer ny2
  real p1(maxnp)
  real p2(maxnp)
  character ( len = 20 ) region
  real nu_inv
  integer row_factor
  real u1(maxnp)
  real u2(maxnp)
  real v1(maxnp)
  real v2(maxnp)
  real xc1(maxnp)
  real xc2(maxnp)
  real yc1(maxnp)
  real yc2(maxnp)
!
  call timestamp ( )
!
!  Say hello.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'REFINE'
  write ( *, '(a)' ) '  Extend a flow solution to a refined grid.'
  write ( *, '(a)' ) '  Last modified on 21 October 1999.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Main limits:'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '    MAXNX =  ', maxnx
  write ( *, '(a,i6)' ) '    MAXNY =  ', maxny
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Derived limits:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    MAXELM = ', maxelm
  write ( *, '(a,i6)' ) '    MAXEQN = ', maxeqn
  write ( *, '(a,i6)' ) '    MAXNP =  ', maxnp
!
!  Read the data.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Specify the input COEF file:'
  write ( *, '(a)' ) '  (Probably "G.TXT")'
  read ( *, '(a)' ) file_name
  if ( file_name == ' ' ) then
    file_name = 'g.txt'
  end if

  call coef_read ( file_name, g1, ierror, maxeqn, neqn1, nx1, ny1, region, &
    nu_inv )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Specify the input ELEMENT file:'
  write ( *, '(a)' ) '  (Probably "ELEMENTS.TXT").'
  read ( *, '(a)' ) file_name

  if ( file_name == ' ' ) then
    file_name = 'elements.txt'
  end if

  call element_read ( file_name, ierror, maxelm, maxnp, maxnpe, nelem1, &
    node1, np1, npe )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Specify the input NODE file:'
  write ( *, '(a)' ) '  (Probably "UVP.TXT").'
  read ( *, '(a)' ) file_name

  if ( file_name == ' ' ) then
    file_name = 'uvp.txt'
  end if

  call node_read ( file_name, ierror, np1, p1, u1, v1, xc1, yc1 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Specify the horizontal and vertical refinement factors.'
  read ( *, * ) col_factor, row_factor

  if ( col_factor < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'REFINE - Fatal error!'
    write ( *, '(a)' ) '  Illegal column refinement ratio.'
    stop
  end if

  if ( row_factor < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'REFINE - Fatal error!'
    write ( *, '(a)' ) '  Illegal row refinement ratio.'
    stop
  end if
!
!  Reconstruct the node-to-unknowns mapping vector INDX.
!
  call indx_set ( indx1, maxelm, maxnp, maxnpe, nelem1, neqn1, node1, np1 )
!
!  Evaluate (U,V,P) at the refinement points.
!
  call coef_extend ( col_factor, g1, ierror, indx1, kind, maxelm, maxeqn, &
    maxnp, maxnpe, nelem1, node1, np2, p2, row_factor, u2, v2, xc1, xc2, &
    yc1, yc2 )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'REFINE - Fatal error!'
    write ( *, '(a)' ) '  Error return from COEF_EXTEND.'
    stop
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) 'Initial number of refined nodes is NP2 = ', np2
!
!  Assume that the proper node ordering can be recovered by sorting
!  on (X,Y) coordinates.
!
  call xyvec_sort_heap_index ( np2, xc2, yc2, iperm )

  call ivec_permute ( np2, kind, iperm )
  call rvec_permute ( np2, xc2, iperm )
  call rvec_permute ( np2, yc2, iperm )
  call rvec_permute ( np2, u2, iperm )
  call rvec_permute ( np2, v2, iperm )
  call rvec_permute ( np2, p2, iperm )
!
!  Drop multiple copies of the same point.
!  Detect them by looking for close (X,Y) coordinates.
!
  mp2 = 1

  do i = 2, np2

    distmin = 0.0E+00

    do j = 1, mp2

      dist = abs ( xc2(i) - xc2(j) ) + abs ( yc2(i) - yc2(j) ) 
 
      if ( j == 1 ) then
        distmin = dist  
      else 
        distmin = min ( distmin, dist )
      end if

    end do
  
    if ( distmin > 0.00001E+00 ) then

      mp2 = mp2 + 1

      kind(mp2) = kind(i)
      p2(mp2) = p2(i)
      u2(mp2) = u2(i)
      v2(mp2) = v2(i)
      xc2(mp2) = xc2(i)
      yc2(mp2) = yc2(i)

    end if

  end do

  np2 = mp2
!
!  Construct INDX2.
!
  ieqn2 = 0

  do ip = 1, np2

    ieqn2 = ieqn2 + 1
    indx2(1,ip) = ieqn2
    ieqn2 = ieqn2 + 1
    indx2(2,ip) = ieqn2
    if ( kind(ip) == 3 ) then
      ieqn2 = ieqn2 + 1
      indx2(3,ip) = ieqn2
    else
      indx2(3,ip) = 0
    end if

  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) 'Final number of refined nodes is NP2 = ', np2
!
!  Construct the coefficient vector.
!
  neqn2 = 0
  do i = 1, np2
    neqn2 = neqn2 + 1
    g2(neqn2) = u2(i)
    neqn2 = neqn2 + 1
    g2(neqn2) = v2(i)
    if ( kind(i) == 3 ) then
      neqn2 = neqn2 + 1
      p2(neqn2) = p2(i)
    end if
  end do
!
!  Write out the coefficient vector.
!
  file_name = 'g2.txt'
  nx2 = nx1 * col_factor
  ny2 = ny1 * row_factor

  call coef_write ( file_name, g2, ierror, maxeqn, neqn2, nx2, ny2, region, &
    nu_inv )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Wrote refined coefficient file ' // trim ( file_name )
!
!  Write out the UVP file.
!
  file_name = 'uvp2.txt'

  call node_write ( file_name, g2, ierror, indx2, maxeqn, maxnp, np2, p2, &
    xc2, yc2 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Wrote refined node data file ' // trim ( file_name )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'REFINE'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine coef_extend ( col_factor, g1, ierror, indx1, kind, maxelm, maxeqn, &
  maxnp, maxnpe, nelem1, node1, np2, p2, row_factor, u2, v2, xc1, xc2, yc1, &
  yc2 )
!
!*******************************************************************************
!
!! COEF_EXTEND extends a vector of solution coefficients to a refined grid.
!
!
!  Discussion:
!
!    The refined grid is assumed to be related to the coarse grid in
!    a simple way.  The number of elements in the X direction must have
!    been increased by an integral factor, and each element subdivided
!    by that factor in the same way.  A similar process must have occurred
!    in the Y direction.  However, the X and Y refinement factors may be
!    different, and either one may be 1.
!
!  Modified:
!
!    21 October 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer COL_FACTOR, the X-direction refinement factor.  This
!    must be an integer.  A value of 2 means that twice as many elements
!    are to be used in the X direction as before.
!
!    Input, real G1(MAXEQN), contains the finite element coefficients of
!    the solution that was computed on the coarse grid.
!
!    Output, integer IERROR, error flag.
!    0, no error occurred.
!    nonzero, an error occurred.
!
!    Input, integer INDX1(3,MAXNP), mapping from nodes to degrees of freedom.
!
!    Output, integer KIND(MAXNP), is 2 for a node associated with U and V 
!    only, or 3 for a node associated with U, V and P.
!
!    Input, integer MAXELM, the maximum number of elements.
!
!    Input, integer MAXEQN, the maximum number of equations.
!
!    Input, integer MAXNP, the maximum number of nodes.
!
!    Input, integer NELEM1, the number of elements in the coarse grid.
!
!    Input, integer NODE1(MAXNPE,MAXELM), the nodes that make up each element.
!
!    Output, integer NP2, the number of nodes in the refined grid.  One
!    physical node may occur in two or more elements, and will be counted
!    several times in NP2.  These duplicate nodes must be weeded out later.
!
!    Output, P2(MAXNP2), the value of the pressure at the refined nodes.
!
!    Input, integer ROW_FACTOR, the Y-direction refinement factor.  This
!    must be an integer.  A value of 2 means that twice as many elements
!    are to be used in the Y direction as before.
!
!    Output, U2(MAXNP2), V2(MAXNP2), the value of the horizontal and vertical
!    velocities at the refined nodes.
!
!    Input, real XC1(MAXNP), the X coordinates of nodes in the coarse grid.
!
!    Input, real XC2(MAXNP), the X coordinates of nodes in the refined grid.
!
!    Input, real YC1(MAXNP), the Y coordinates of nodes in the coarse grid.
!
!    Input, real YC2(MAXNP), the Y coordinates of nodes in the refined grid.
!
  implicit none
!
  integer maxelm
  integer maxeqn
  integer maxnp
  integer maxnpe
!
  integer col_factor
  real eta
  real g1(maxeqn)
  integer i
  integer ielem1
  integer ierror
  integer indx1(3,maxnp)
  integer j
  integer kind(maxnp)
  integer nelem1
  integer node_local_col_max
  integer node_local_row_max
  integer node1(maxnpe,maxelm)
  integer np2
  real p
  real p2(maxnp)
  integer row_factor
  real u
  real u2(maxnp)
  real v
  real v2(maxnp)
  real x
  real xc1(maxnp)
  real xc2(maxnp)
  real xsi
  real y
  real yc1(maxnp)
  real yc2(maxnp)
!
  ierror = 0
  node_local_col_max = 2 * col_factor + 1
  node_local_row_max = 2 * row_factor + 1

  np2 = 0

  do ielem1 = 1, nelem1

    do i = 1, node_local_col_max
      xsi = real ( i - 1 ) / real ( node_local_col_max - 1 )
      do j = 1, node_local_row_max
        eta = real ( j - 1 ) / real ( node_local_row_max - 1 )

        if ( xsi + eta <= 1.0E+00 ) then

          np2 = np2 + 1

          if ( np2 > maxnp ) then
            ierror = 1
            write ( *, '(a)' ) ' '
            write ( *, '(a)' ) 'COEF_EXTEND - Fatal error!'
            write ( *, '(a)' ) '  MAXNP is too small.'
            return
          end if

          call uvp_value ( eta, g1, ielem1, indx1, maxelm, maxeqn, &
            maxnp, node1, p, u, v, xsi )

          call x_of_xsi ( eta, ielem1, maxelm, maxnp, node1, x, xc1, xsi, y, &
            yc1 )

          u2(np2) = u
          v2(np2) = v
          p2(np2) = p

          if ( mod ( i, 2 ) == 1 .and. mod ( j, 2 ) == 1 ) then
            kind(np2) = 3
          else
            kind(np2) = 2
          end if

          xc2(np2) = x
          yc2(np2) = y

        end if

      end do
    end do
    
  end do

  return
end
subroutine coef_read ( file_name, g, ierror, maxeqn, neqn, nx, ny, region, &
  nu_inv )
!
!*******************************************************************************
!
!! COEF_READ reads the coefficient data from a file.
!
!
!  Modified:
!
!    20 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real G(MAXEQN), the current solution vector.
!
!    Output, integer IERROR, error flag.
!    0, no error occurred.
!    nonzero, an error occurred.
!
!    Input, integer MAXEQN, the maximum number of equations.
!
!    Output, integer NEQN, the number of equations.
!
!    Output, integer NX, NY, the number of elements along a horizontal, or
!    vertical line.
!
!    Output, character ( len = 20 ) REGION, the flow problem.
!
!    Output, real NU_INV, the value of the inverse dynamic viscosity.
!
  implicit none
!
  integer maxeqn
!
  character ( len = * ) file_name
  real g(maxeqn)
  integer i
  integer ierror
  integer ios
  integer iunit
  integer neqn
  integer nx
  integer ny
  character ( len = * ) region
  real nu_inv
!
  ierror = 0

  neqn = 0
  nx = 0
  ny = 0
  region = ' '
  nu_inv = 0.0E+00

  call get_unit ( iunit )

  open ( unit = iunit, file = file_name, status = 'old', iostat = ios )

  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'COEF_READ - Warning!'
    write ( *, '(a)' ) '  Could not open the file:'
    write ( *, '(a)' ) trim ( file_name )
    return
  end if

  do i = 1, 10
    read ( iunit, '(1x)', iostat = ios )
    if ( ios /= 0 ) then
      ierror = 1
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'COEF_READ - Warning!'
      write ( *, '(a)' ) '  Error while reading the coefficient file.'
      return
    end if
  end do

  read ( iunit, *, iostat = ios ) neqn
  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'COEF_READ - Warning!'
    write ( *, '(a)' ) '  Error while reading the coefficient file.'
    return
  end if

  read ( iunit, *, iostat = ios ) nx
  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'COEF_READ - Warning!'
    write ( *, '(a)' ) '  Error while reading the coefficient file.'
    return
  end if

  read ( iunit, *, iostat = ios ) ny
  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'COEF_READ - Warning!'
    write ( *, '(a)' ) '  Error while reading the coefficient file.'
    return
  end if

  read ( iunit, '(a)', iostat = ios ) region
  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'COEF_READ - Warning!'
    write ( *, '(a)' ) '  Error while reading the coefficient file.'
    return
  end if

  read ( iunit, *, iostat = ios ) nu_inv
  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'COEF_READ - Warning!'
    write ( *, '(a)' ) '  Error while reading the coefficient file.'
    return
  end if

  read ( iunit, '(1x)', iostat = ios )
  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'COEF_READ - Warning!'
    write ( *, '(a)' ) '  Error while reading the coefficient file.'
    return
  end if

  if ( neqn > maxeqn ) then
    ierror = 3
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'COEF_READ - Warning!'
    write ( *, '(a)' ) '  Input NEQN = ', neqn
    write ( *, '(a)' ) '  Internal MAXEQN = ', maxeqn
    close ( unit = iunit )
    return
  end if

  do i = 1, neqn
    read ( iunit, *, iostat = ios ) g(i)
    if ( ios /= 0 ) then
      ierror = 1
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'COEF_READ - Warning!'
      write ( *, '(a)' ) '  Error while reading the coefficient file.'
      return
    end if
  end do

  close ( unit = iunit )

  return
end
subroutine coef_write ( file_name, g, ierror, maxeqn, neqn, nx, ny, region, &
  nu_inv )
!
!*******************************************************************************
!
!! COEF_WRITE writes the coefficient data to a file for possible restart.
!
!
!  Modified:
!
!    18 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real G(MAXEQN), the current solution vector.
!
!    Output, integer IERROR, error flag.
!    0, no error occurred.
!    nonzero, an error occurred.
!
!    Input, integer MAXEQN, the maximum number of equations.
!
!    Input, integer NEQN, the number of equations.
!
!    Input, integer NX, NY, the number of elements along a horizontal, or
!    vertical line.
!
!    Input, character ( len = 20 ) REGION, the flow problem.
!
!    Input, real NU_INV, the value of the inverse dynamic viscosity.
!
  implicit none
!
  integer maxeqn
!
  character ( len = * ) file_name
  real g(maxeqn)
  integer i
  integer ierror
  integer ios
  integer iunit
  integer lenc
  integer neqn
  integer nx
  integer ny
  character ( len = * ) region
  real nu_inv
!
  ierror = 0

  call get_unit ( iunit )

  open ( unit = iunit, file = file_name, status = 'replace', &
    iostat = ios )

  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'COEF_WRITE - Warning!'
    write ( *, '(a)' ) '  Could not open the file.'
    return
  end if

  write ( iunit, '(a)' ) '#  g.txt  FLOW6 solution vector'
  write ( iunit, '(a)' ) '#  '
  write ( iunit, '(a)' ) '#  NEQN'
  write ( iunit, '(a)' ) '#  NX'
  write ( iunit, '(a)' ) '#  NY'
  write ( iunit, '(a)' ) '#  REGION'
  write ( iunit, '(a)' ) '#  NU_INV'
  write ( iunit, '(a)' ) '#  '
  write ( iunit, '(a)' ) '#  (G(I),I=1,NEQN)'
  write ( iunit, '(a)' ) '#'

  write ( iunit, '(i6)' ) neqn
  write ( iunit, '(i6)' ) nx
  write ( iunit, '(i6)' ) ny
  write ( iunit, '(a)' ) trim ( region )
  write ( iunit, '(g14.6)' ) nu_inv
  write ( iunit, '(1x)' )

  do i = 1, neqn
    write ( iunit, '(g14.6)' ) g(i)
  end do

  close ( unit = iunit )

  return
end
subroutine element_read ( file_name, ierror, maxelm, maxnp, maxnpe, nelem, &
  node, np, npe )
!
!*******************************************************************************
!
!! ELEMENT_READ reads an element file.
!
!
!  Format:
!
!    The element file contains information about the organization of
!    the nodes into elements.  The format is as follows:
!
!      Line 1: NELEM, the number of elements
!      Line 2: NPE, the number of nodes per element.
!      Line I+2 through NELEM+2: the nodes in element I.
!
!  Modified:
!
!    28 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILELM, the name of the file
!    containing the element information.
!
!    Input, integer MAXELM, the maximum number of elements.
!
!    Input, integer MAXNP, the maximum number of nodes.
!
!    Input, integer MAXNPE, the maximum number of nodes per element.
!
!    Output, integer NELEM, the number of elements.
!
!    Input, integer NODE(MAXNPE,MAXELM), the global node numbers in each element
!
!    Output, integer NP, the number of nodes.
!
!    Output, integer NPE, the number of nodes per element.
!
  implicit none
!
  integer maxelm
  integer maxnpe
!
  character ( len = * ) file_name
  integer i
  integer ielem
  integer ierror
  integer input_unit
  integer ios
  integer itemp(6)
  integer j
  integer maxnp
  integer nelem
  integer node(maxnpe,maxelm)
  integer np
  integer npe
!
  ierror = 0

  call get_unit ( input_unit )
!
!  Open the data file.
!
  open ( unit = input_unit, file = file_name, status = 'old', iostat = ios )

  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'ELEMENT_READ - Fatal error!'
    write ( *, '(a)' ) '  Could not open the element file!'
    return
  end if
!
!  Initialize number of elements, number of nodes.
!
  ielem = 0
  np = 0
!
!  Read the number of elements.
!
  read ( input_unit, *, end = 30, err = 40 ) nelem

  if ( nelem <= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'ELEMENT_READ - Fatal error!'
    write ( *, '(a)' ) '  This problem has zero elements!'
    return
  end if

  if ( nelem > maxelm ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'ELEMENT_READ - Fatal error!'
    write ( *, '(a)' ) '  This problem has too many elements!'
    write ( *, '(a,i6)' ) '  Number of elements NELEM = ', nelem
    write ( *, '(a,i6)' ) '  DISPLAY can handle up to MAXELM = ', maxelm
    return
  end if
!
!  Read the number of nodes per element.
!
  read ( input_unit, *, end = 30, err = 40 ) npe

  if ( npe /= 3 .and. npe /= 4 .and. npe /= 6 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'ELEMENT_READ - Fatal error!'
    write ( *, '(a)' ) '  Legal values of NPE are 3, 4 and 6.'
    write ( *, '(a,i6)' ) '  User input value is NPE = ', npe
    return
  end if
!
!  Try to read the next line of the file, which should contain
!  the node numbers associated with the next element.
!
10    continue

  read ( input_unit, *, end = 30, err = 40 ) itemp(1:npe)

  ielem = ielem + 1

  if ( ielem > nelem ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'ELEMENT_READ - Fatal error!'
    write ( *, '(a)' ) '  This problem has too many elements!'
    write ( *, '(a,i6)' ) '  Current element number IELEM = ', ielem
    write ( *, '(a,i6)' ) '  Declared number of elements NELEM = ', nelem
    return
  end if

  do j = 1, npe

    if ( itemp(j) <= 0 ) then
      ierror = 1
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'ELEMENT_READ - Fatal error!'
      write ( *, '(a)' ) '  Element ', ielem, ' has an illegal node.'
      write ( *, '(a,i6,a,i6)' ) '  Local node number ', j , &
        ' has index ', itemp(j)
      return
    end if

  end do

  do i = 1, npe
    np = max ( np, itemp(i) )
    node(i,ielem) = itemp(i)
  end do

  go to 10
!
!  End of file.
!
30    continue

  close ( unit = input_unit )

  if ( nelem /= ielem ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'ELEMENT_READ - Fatal error!'
    write ( *, '(a,i6)' ) '  Number of elements declared was NELEM = ',nelem
    write ( *, '(a,i6)' ) '  Number of elements found was IELEM = ', ielem
    return
  end if

  if ( np > maxnp ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'ELEMENT_READ - Fatal error!'
    write ( *, '(a)' ) '  This problem has too many nodes!'
    write ( *, '(a,i6)' ) '  Number of nodes NP = ',np
    write ( *, '(a,i6)' ) '  The program can handle up to MAXNP = ', maxnp
    return
  end if

  if ( np <= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'ELEMENT_READ - Fatal error!'
    write ( *, '(a,i6)' ) '  This problem has NP = ',np
    write ( *, '(a)' ) '  but the number of nodes must be positive!'
    return
  end if
!
!  All seems well.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'ELEMENT_READ - Note:'
  write ( *, '(a,i6,a)' ) '  There are NELEM = ', nelem, ' elements,'
  write ( *, '(a,i6,a)' ) '  and NP = ', np, ' nodes.'

  return
!
!  Error while reading.
!
40    continue

  close ( unit = input_unit )
  ierror = 1
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'ELEMENT_READ - Fatal error!'
  write ( *, '(a)' ) '  A READ "ERR" condition occurred.'
  return
end
subroutine get_unit ( iunit )
!
!*******************************************************************************
!
!! GET_UNIT returns a free FORTRAN unit number.
!
!
!  Discussion:
!
!    A "free" FORTRAN unit number is an integer between 1 and 99 which
!    is not currently associated with an I/O device.  A free FORTRAN unit
!    number is needed in order to open a file with the OPEN command.
!
!  Modified:
!
!    02 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IUNIT.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5 and 6).
!
!    Otherwise, IUNIT is an integer between 1 and 99, representing a
!    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
  implicit none
!
  integer i
  integer ios
  integer iunit
  logical lopen

  iunit = 0

  do i = 1, 99

    if ( i /= 5 .and. i /= 6 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )

      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if

  end do

  return
end
subroutine indx_set ( indx, maxelm, maxnp, maxnpe, nelem, neqn, node, np )
!
!*******************************************************************************
!
!! INDX_SET determines indices and other data from the NODE array.
!
!
!  Modified:
!
!    16 September 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer INDX(3,MAXNP), mapping from nodes to degrees of freedom.
!
!    Input, integer MAXELM, the maximum number of elements.
!
!    Input, integer MAXNP, the maximum number of nodes.
!
!    Input, integer MAXNPE, the maximum number of nodes per element.
!
!    Input, integer NELEM, the number of elements.
!
!    Input, integer NEQN, the number of equations.
!
!    Input, integer NODE(MAXNPE,MAXELM), the nodes that make up each element.
!
!    Input, integer NP, the number of nodes.
!
  implicit none
!
  integer maxelm
  integer maxnp
  integer maxnpe
!
  integer i
  integer indx(3,maxnp)
  integer j
  integer k
  integer n
  integer nelem
  integer neqn
  integer node(maxnpe,maxelm)
  integer np
!
  do i = 1, 3
    do j = 1, maxnp
      indx(i,j) = 0
    end do
  end do
!
!  Now set INDX(1,N) to 1 if the node is a corner (UVP) node.
!
  do j = 1, nelem
    do i = 1, 3
      n = node(i,j)
      indx(1,n) = 1
    end do
  end do
!
!  Count the nodes.
!
  np = 0
  do j = 1, nelem
    do i = 1, 6
      n = node(i,j)
      np = max ( np, n )
    end do
  end do
!
!  Now we can march through the nodes, setting aside just enough
!  space for the unknowns.
!
  neqn = 0

  do i = 1, np

    k = indx(1,i)
    neqn = neqn + 1
    indx(1,i) = neqn
    neqn = neqn + 1
    indx(2,i) = neqn

    if ( k > 0 ) then
      neqn = neqn + 1
      indx(3,i) = neqn
    end if
    
  end do

  return
end
subroutine ivec_permute ( n, x, perm )
!
!*******************************************************************************
!
!! IVEC_PERMUTE permutes an integer vector in place.
!
!
!  Note:
!
!    This routine permutes an array of integer "objects", but the same
!    logic can be used to permute an array of objects of any arithmetic
!    type, or an array of objects of any complexity.  The only temporary
!    storage required is enough to store a single object.  The number
!    of data movements made is N + the number of cycles of order 2 or more,
!    which is never more than N + N/2.
!
!  Example:
!
!    Input:
!
!      N = 5
!      PERM = (   2,   4,   5,   1,   3 )
!      X    = (   1,   2,   3,   4,   5 )
!
!    Output:
!
!      X    = (   2,   4,   5,   1,   3 ).
!
!  Modified:
!
!    26 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of objects.
!
!    Input/output, integer X(N), the array to be permuted.
!
!    Input, integer PERM(N), the permutation.  PERM(I) = J means
!    that the I-th element of the output array should be the J-th
!    element of the input array.  PERM must be a legal permutation
!    of the integers from 1 to N, otherwise the algorithm will
!    fail catastrophically.
!
  implicit none
!
  integer n
!
  integer i
  integer iget
  integer iput
  integer istart
  integer perm(n)
  integer x(n)
  integer x_temp
!
  istart = 0
!
!  Search for the next element of PERM that has not been used.
!
10    continue

  istart = istart + 1

  if ( istart > n ) then

    do i = 1, n
      perm(i) = - perm(i)
    end do

    return

  end if

  if ( perm(istart) < 0 ) then
    go to 10
  end if

  if ( perm(istart) == istart ) then
    perm(istart) = - perm(istart)
    go to 10
  end if

  x_temp = x(istart)
  iget = istart
!
!  Copy the new value into the vacated entry.
!
20    continue

  iput = iget
  iget = perm(iget)

  perm(iput) = - perm(iput)

  if ( iget < 1 .or. iget > n ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'IVEC_PERMUTE - Fatal error!' 
    write ( *, '(a)' ) '  An illegal entry was encountered:'
    write ( *, '(a,i6,a,i6)' ) '  Entry ', iput, ' had the value ', iget
    write ( *, '(a,i6)' ) '  but values must be between 1 and ', n
    stop
  else if ( iget /= istart ) then
    x(iput) = x(iget)
    go to 20
  else
    x(iput) = x_temp
    go to 10
  end if

end
subroutine node_read ( file_name, ierror, np, p, u, v, xc, yc )
!
!*******************************************************************************
!
!! NODE_READ reads a file containing information about data at nodes.
!
!
!  Modified:
!
!    28 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILDAT, the file containing the data.
!    This file should be a "plain text" file.
!
!    Input, integer NP, the number of nodes.
!
!    Output, real P(NP), RHO(NP), U(NP), V(NP),
!    The pressure, density, horizontal and vertical velocities of each node.
!
!    Output, real XC(NP), the X coordinates of the nodes.
!
!    Output, real YC(NP), the Y coordinates of the nodes.
!
  implicit none
!
  integer np
!
  character ( len = * ) file_name
  integer i
  integer ierror
  integer input_unit
  integer ios
  real p(np)
  real u(np)
  real v(np)
  real xc(np)
  real yc(np)
!
  ierror = 0

  call get_unit ( input_unit )

  open ( unit = input_unit, file = file_name, status = 'old', iostat = ios )

  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NODE_READ - Fatal error!'
    write ( *, '(a)' ) '  Could not open the node file.'
    stop
  end if

  do i = 1, np

    read ( input_unit, *, iostat = ios ) xc(i), yc(i), u(i), v(i), p(i)

    if ( ios /= 0 ) then
      ierror = 2
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NODE_READ - Fatal error!'
      write ( *, '(a)' ) '  READ error while reading node data.'
      stop
    end if

  end do

  close ( unit = input_unit )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'NODE_READ - Note:'
  write ( *, '(a)' ) '  Read the node data.'

  return
end
subroutine node_write ( file_name, g, ierror, indx, maxeqn, maxnp, np, p, &
  xc, yc )
!
!*******************************************************************************
!
!! NODE_WRITE writes solution information to a file.
!
!
!  Modified:
!
!    28 June 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file.
!
!    Input, real G(MAXEQN), the current solution vector.
!
!    Output, integer IERROR, error flag.
!    0, no error occurred.
!    nonzero, an error occurred.
!
!    Input, integer MAXNP, the maximum number of nodes.
!
!    Input, integer NP, the number of nodes.
!
!    Input, real P(MAXNP), the pressure at each node.
!
!    Input, real XC(MAXNP), YC(MAXNP), the coordinates of the nodes.
!
  implicit none
!
  integer maxeqn
  integer maxnp
!
  character ( len = * ) file_name
  real g(maxeqn)
  integer i
  integer ierror
  integer indx(3,maxnp)
  integer ios
  integer iunit
  integer np
  real p(maxnp)
  real xc(maxnp)
  real yc(maxnp)
!
  ierror = 0

  call get_unit ( iunit )

  open ( unit = iunit, file = file_name, status = 'replace', &
    iostat = ios )

  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NODE_WRITE - Warning!'
    write ( *, '(a)' ) '  Could not open the node file.'
    return
  end if

  do i = 1, np
    write ( iunit, '(2g12.4,3g14.6)' ) xc(i), yc(i), g(indx(1,i)), &
      g(indx(2,i)), p(i)
  end do

  close ( unit = iunit )

  return
end
subroutine ref_bf_l3 ( q, dqdeta, dqdxsi, eta, iq, xsi )
!
!*******************************************************************************
!
!! REF_BF_L3 evaluates a reference element linear basis function.
!
!
!  Diagram:
!
!          ^
!          |
!        1 +    2
!          |    |\
!    ETA   |    | \
!          |    |  \
!        0 +    3---1
!          |
!          +----+---+--->
!               0   1
!
!              XSI
!
!  Discussion:
!
!    The basis function and its X and Y derivatives are evaluated at a
!    particular point (X,Y) in a particular element, by referring to the
!    corresponding points (XSI,ETA) in the reference triangle.
!
!  Modified:
!
!    28 September 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real Q, DQDETA, DQDXSI, the value of the basis
!    function, and its derivatives with respect to ETA and XSI, at
!    the point with reference coordinates (ETA,XSI).
!
!    Input, real ETA, XSI, the local coordinates of the
!    point at which the basis information is desired.
!
!    Input, integer IQ, the basis function, between 1 and 3.
!
  implicit none
!
  real dqdeta
  real dqdxsi
  real eta
  integer iq
  real q
  real xsi
!
  if ( iq == 1 ) then

    q = xsi
    dqdxsi = 1.0
    dqdeta =  0.0

  else if ( iq == 2 ) then

    q = eta
    dqdxsi = 0.0
    dqdeta = 1.0

  else if ( iq == 3 ) then

    q = 1.0 - xsi - eta
    dqdxsi = - 1.0
    dqdeta = - 1.0

  else

    q = 0.0
    dqdxsi = 0.0
    dqdeta = 0.0

  end if

  return
end
subroutine ref_bf_q6 ( w, dwdeta, dwdxsi, eta, iq, xsi )
!
!*******************************************************************************
!
!! REF_BF_Q6 evaluates shape functions for a 6 node triangle.
!
!
!  Diagram:
!
!    |
!    1  2
!    |  |\
!    |  | \
!    S  6  5
!    |  |   \
!    |  |    \
!    0  3--4--1
!    |
!    +--0--R--1-->
!
!  Modified:
!
!    04 October 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real W, DWDETA, DWDXSI, the value of the basis function, and
!    its derivatives with respect to ETA and XSI, at the point (XSI,ETA).
!
!    Input, real ETA, the ETA coordinate of the point.
!
!    Input, integer IQ, the basis function, between 1 and 6.
!
!    Input, real XSI, the XSI coordinate of the point.
!
  implicit none
!
  real dwdeta
  real dwdxsi
  real eta
  integer iq
  real w
  real xsi
!
  if ( iq == 1 ) then

    w =        2.0 * xsi * ( xsi - 0.5 )
    dwdxsi = - 1.0 + 4.0 * xsi
    dwdeta =   0.0

  else if ( iq == 2 ) then

    w =        2.0 * eta * ( eta - 0.5 )
    dwdxsi =   0.0
    dwdeta = - 1.0           + 4.0 * eta

  else if ( iq == 3 ) then

    w =        2.0 * ( 1.0 - xsi - eta ) * ( 0.5 - xsi - eta )
    dwdxsi = - 3.0 + 4.0 * xsi + 4.0 * eta
    dwdeta = - 3.0 + 4.0 * xsi + 4.0 * eta

  else if ( iq == 4 ) then

    w =        4.0 * xsi * ( 1.0 - xsi - eta )
    dwdxsi =   4.0 - 8.0 * xsi - 4.0 * eta
    dwdeta =       - 4.0 * xsi

  else if ( iq == 5 ) then

    w =      4.0 * xsi * eta
    dwdxsi = 4.0       * eta
    dwdeta = 4.0 * xsi

  else if ( iq == 6 ) then

    w =      4.0 * eta * ( 1.0 - xsi - eta )
    dwdxsi =                 - 4.0 * eta
    dwdeta = 4.0 - 4.0 * xsi - 8.0 * eta

  else

    w = 0.0
    dwdxsi = 0.0
    dwdeta = 0.0

  end if

  return
end
subroutine ref_map_q6 ( qdata, a, b, c, d, e, f )
!
!*******************************************************************************
!
!! REF_MAP_Q6 returns the interpolation map for data on a 6 node triangle.
!
!
!  Diagram:
!
!    |
!    1  2
!    |  |\
!    E  | \
!    T  6  5
!    A  |   \
!    |  |    \
!    0  3--4--1
!    |
!    +--0-XSI-1-->
!
!  Formula:
!
!    Q(R,S) = A * XSI**2 + B * XSI * ETA + C * ETA**2 + D * XSI + E * ETA + F
!
!  Note:
!
!    To find the polynomial form of, say, the third basis function,
!    set QDATA = ( 0, 0, 1, 0, 0, 0 ).
!
!  Modified:
!
!    14 October 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real QDATA(6), the values of a quantity at the basis nodes.
!
!    Output, real A, B, C, D, E, F, the polynomial coefficients of the
!    interpolation function over the element.
!
  implicit none
!
  real a
  real b
  real c
  real d
  real e
  real f
  real qdata(6)
!
  a =  2.0 * qdata(1) - 4.0 * qdata(4) + 2.0 * qdata(3)
  b =  4.0 * qdata(3) - 4.0 * qdata(4) + 4.0 * qdata(5) - 4.0 * qdata(6)
  c =  2.0 * qdata(2) - 4.0 * qdata(6) + 2.0 * qdata(3)
  d =      - qdata(1) + 4.0 * qdata(4) - 3.0 * qdata(3)
  e =      - qdata(2) + 4.0 * qdata(6) - 3.0 * qdata(3)
  f =        qdata(3)

  return
end
subroutine rvec_permute ( n, x, perm )
!
!*******************************************************************************
!
!! RVEC_PERMUTE permutes a real vector in place.
!
!
!  Note:
!
!    This routine permutes an array of real "objects", but the same
!    logic can be used to permute an array of objects of any arithmetic
!    type, or an array of objects of any complexity.  The only temporary
!    storage required is enough to store a single object.  The number
!    of data movements made is N + the number of cycles of order 2 or more,
!    which is never more than N + N/2.
!
!  Example:
!
!    Input:
!
!      N = 5
!      PERM = (   2,   4,   5,   1,   3 )
!      X    = ( 1.0, 2.0, 3.0, 4.0, 5.0 )
!
!    Output:
!
!      X    = ( 2.0, 4.0, 5.0, 1.0, 3.0 ).
!
!  Modified:
!
!    26 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of objects.
!
!    Input/output, real X(N), the array to be permuted.
!
!    Input, integer PERM(N), the permutation.  PERM(I) = J means
!    that the I-th element of the output array should be the J-th
!    element of the input array.  PERM must be a legal permutation
!    of the integers from 1 to N, otherwise the algorithm will
!    fail catastrophically.
!
  implicit none
!
  integer n
!
  integer i
  integer iget
  integer iput
  integer istart
  integer perm(n)
  real x(n)
  real x_temp
!
  istart = 0
!
!  Search for the next element of PERM that has not been used.
!
10    continue

  istart = istart + 1

  if ( istart > n ) then
    perm(1:n) = - perm(1:n)
    return
  end if

  if ( perm(istart) < 0 ) then
    go to 10
  end if

  if ( perm(istart) == istart ) then
    perm(istart) = - perm(istart)
    go to 10
  end if

  x_temp = x(istart)
  iget = istart
!
!  Copy the new value into the vacated entry.
!
20    continue

  iput = iget
  iget = perm(iget)

  perm(iput) = - perm(iput)

  if ( iget < 1 .or. iget > n ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'RVEC_PERMUTE - Fatal error!'
    stop
  else if ( iget /= istart ) then
    x(iput) = x(iget)
    go to 20
  else
    x(iput) = x_temp
    go to 10
  end if

end
subroutine uvp_value ( eta, g, ielem, indx, maxelm, maxeqn, maxnp, node, &
  p, u, v, xsi )
!
!*******************************************************************************
!
!! UVP_VALUE evaluates U, V and P at any point in a given element.
!
!
!  Modified:
!
!    19 February 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ETA, the ETA coordinate of the point.
!
!    Input, real G(MAXEQN), the current solution vector.
!
!    Input, integer IELEM, the index of the element.
!
!    Input, integer INDX(3,MAXNP), mapping from nodes to degrees of freedom.
!
!    Input, integer MAXELM, the maximum number of elements.
!
!    Input, integer MAXEQN, the maximum number of equations.
!
!    Input, integer MAXNP, the maximum number of nodes.
!
!    Input, integer NODE(6,MAXELM), the nodes that make up each element.
!
!    Output, real P, the value of the pressure.
!
!    Output, real U, the value of the horizontal velocity.
!
!    Output, real V, the value of the vertical velocity.
!
!    Input, real XSI, the XSI coordinate of the point.
!
  implicit none
!
  integer maxelm
  integer maxeqn
  integer maxnp
!
  real dqdeta
  real dqdxsi
  real dwdeta
  real dwdxsi
  real eta
  real g(maxeqn)
  integer icof
  integer ielem
  integer indx(3,maxnp)
  integer ip
  integer iq
  integer node(6,maxelm)
  real p
  real q
  real u
  real v
  real w
  real xsi
!
  p = 0.0

  do iq = 1, 3

    call ref_bf_l3 ( q, dqdeta, dqdxsi, eta, iq, xsi )

    ip = node(iq,ielem)

    icof = indx(3,ip)
    p = p + g(icof) * q

  end do

  u = 0.0
  v = 0.0
  do iq = 1, 6

    ip = node(iq,ielem)

    call ref_bf_q6 ( w, dwdeta, dwdxsi, eta, iq, xsi )

    ip = node(iq,ielem)

    icof = indx(1,ip)
    u = u + g(icof) * w

    icof = indx(2,ip)
    v = v + g(icof) * w

  end do

  return
end
subroutine x_of_xsi ( eta, ielem, maxelm, maxnp, node, x, xc, xsi, y, yc )
!
!*******************************************************************************
!
!! X_OF_XSI computes X and Y given XSI and ETA coordinates.
!
!
!  Modified:
!
!    16 September 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ETA, the ETA coordinate of the point.
!
!    Input, integer IELEM, the index of the element.
!
!    Input, integer MAXELM, the maximum number of elements.
!
!    Input, integer MAXNP, the maximum number of nodes.
!
!    Input, integer NODE(6,MAXELM), the nodes that make up each element.
!
!    Output, real X, the X coordinate of the point.
!
!    Input, real XC(MAXNP), the X coordinates of the nodes.
!
!    Input, real XSI, the XSI coordinate of the point.
!
!    Output, real Y, the Y coordinate of the point.
!
!    Input, real YC(MAXNP), the Y coordinates of the nodes.
!
  implicit none
!
  integer maxelm
  integer maxnp
!
  real a
  real b
  real c
  real d
  real e
  real eta
  real f
  integer i
  integer ielem
  integer node(6,maxelm)
  real x
  real xc(maxnp)
  real xsi
  real xx(6)
  real y
  real yc(maxnp)
!
!  Compute X(XSI,ETA).
!
  do i = 1, 6
    xx(i) = xc(node(i,ielem))
  end do

  call ref_map_q6 ( xx, a, b, c, d, e, f )

  x = a * xsi**2 + b * xsi * eta + c * eta**2 + d * xsi + e * eta + f
!
!  Compute Y(XSI,ETA).
!
  do i = 1, 6
    xx(i) = yc(node(i,ielem))
  end do

  call ref_map_q6 ( xx, a, b, c, d, e, f )

  y = a * xsi**2 + b * xsi * eta + c * eta**2 + d * xsi+ e * eta + f

  return
end
subroutine xyvec_sort_heap_index ( n, x, y, indx )
!
!*******************************************************************************
!
!! XYVEC_SORT_HEAP_INDEX does an indexed heap sort of (X,Y) data.
!
!
!  Discussion:
!
!    The sorting is not actually carried out.  Rather an index array is
!    created which defines the sorting.  This array may be used to sort
!    or index the array, or to sort or index related arrays keyed on the
!    original array.
!
!    ( X(I), Y(I) ) < ( X(J), Y(J) ) if:
!
!    * X(I) < X(J), or
!
!    * X(I) = X(J), and Y(I) < Y(J).
!
!    Once the index array is computed, the sorting can be carried out
!    "implicitly:
!
!      ( X(INDX(I)), Y(INDX(I) ), for I = 1 to N is sorted,
!
!    or explicitly, by the call
!
!      call RVEC_PERMUTE ( N, X, INDX )
!      call RVEC_PERMUTE ( N, Y, INDX )
!
!    after which ( X(I), Y(I) ), I = 1 to N is sorted.
!
!  Modified:
!
!    20 October 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in the array.
!
!    Input, real X(N),Y(N), pairs of X, Y coordinates of points.
!
!    Output, integer INDX(N), contains the sort index.  The
!    I-th element of the sorted array has coordinates ( X(INDX(I)), Y(INDX(I) ).
!
  implicit none
!
  integer n
!
  integer i
  integer indx(n)
  integer indxt
  integer ir
  integer j
  integer l
  real x(n)
  real xval
  real y(n)
  real yval
!
  do i = 1, n
    indx(i) = i
  end do

  l = n / 2 + 1
  ir = n

10    continue

  if ( l > 1 ) then

    l = l - 1
    indxt = indx(l)
    xval = x(indxt)
    yval = y(indxt)

  else

    indxt = indx(ir)
    xval = x(indxt)
    yval = y(indxt)
    indx(ir) = indx(1)
    ir = ir - 1

    if ( ir == 1 ) then
      indx(1) = indxt
      return
    end if

  end if

  i = l
  j = l + l

20    continue

  if ( j <= ir ) then

    if ( j < ir ) then

      if ( x(indx(j)) < x(indx(j+1)) .or. &
         ( x(indx(j)) == x(indx(j+1)) .and. y(indx(j)) < y(indx(j+1)) ) ) then
        j = j + 1
      end if

    end if

    if ( xval < x(indx(j)) .or. &
       ( xval == x(indx(j)) .and. yval < y(indx(j)) ) ) then
      indx(i) = indx(j)
      i = j
      j = j + j
    else
      j = ir + 1
    end if

    go to 20

  end if

  indx(i) = indxt
  go to 10
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