subroutine cgrid ( xmin, xmax, nx, ymin, ymax, ny )

!*****************************************************************************80
!
!! CGRID draws evenly spaced Cartesian X and Y grid lines.
!
!  Discussion:
!
!    The color used to draw the grid lines may be changed by calling
!    LINCLR before calling this routine.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real XMIN, XMAX, are the horizontal limits of the location
!    of the grid.
!    XMIN and XMAX should be "within" your picture.  If you are using
!    the simple DRAWCGM coordinate system, which uses 0 as the minimum
!    X value and 1 as the maximum X value, then XMIN should be 0 or
!    greater, and XMAX should be 1 or less.
!    However, if you have used SETWCD or SETSCL to allow a different
!    range of X values, then XMIN and XMAX may be any values in that
!    range.
!    Similar remarks apply to the YMIN and YMAX values.
!
!    Input, integer NX, the number of grid lines to draw along the
!    X direction.
!
!    Input, real YMIN, YMAX, are the vertical limits of the location
!    of the grid.
!
!    Input, integer NY, the number of grid lines to draw along the
!    Y direction.
!
  implicit none

  integer i
  integer nx
  integer ny
  real xmax
  real xmin
  real xval
  real ymax
  real ymin
  real yval

  do i = 1, nx

    if ( nx /= 1 ) then
      xval = ( real ( nx - i ) * xmin + real ( i - 1 ) * xmax ) / &
        real ( nx - 1 )
    else
      xval = 0.5E+00 * ( xmin + xmax )
    end if

    call movcgm ( xval, ymin )
    call drwcgm ( xval, ymax )

  end do

  do i = 1, ny

    if ( ny /= 1 ) then
      yval = ( real ( ny - i ) * ymin + real ( i - 1 ) * ymax ) / &
        real ( ny - 1 )
    else
      yval = 0.5E+00 * ( ymin + ymax )
    end if

    call movcgm ( xmin, yval )
    call drwcgm ( xmax, yval )

  end do

  return
end
subroutine circle ( xcentr, ycentr, radius, filled )

!*****************************************************************************80
!
!! CIRCLE draws an open or filled circle of a given radius.
!
!  Discussion:
!
!    CIRCLE calls PLYLIN to draw an open circle, or PLYGON to draw a
!    filled circle.
!
!    To control the color of the circle, simply call the DRAWCGM
!    routine LINCLR before drawing open circles or FILCLR before
!    drawing closed circles.
!
!  Note:
!
!    If the X and Y coordinate systems do not have the same scale,
!    the circle drawn will be "flattened".  This can happen if the
!    routine SETWCD has been used to set the X and Y dimensions to
!    different extents.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real XCENTR, YCENTR, the coordinates of the circle center.
!
!    Input, real RADIUS, the radius of the circle.
!
!    Input, logical FILLED,
!    .TRUE. if the circle is to be filled.
!    .FALSE. if the circle is to be open.
!
  implicit none
!
  integer, parameter :: npts = 64
!
  real angle
  integer i
  logical filled
  real r_pi
  real radius
  real xcentr
  real xpoint(npts)
  real ycentr
  real ypoint(npts)
!
!  Set up the data defining the circle, using NPTS equally spaced
!  points along the circumference.  The first and the last points
!  are the same.
!
  do i = 1, npts
    angle = 2.0E+00 * r_pi ( ) * real ( i - 1 ) / real ( npts - 1 )
    xpoint(i) = xcentr + radius * sin ( angle )
    ypoint(i) = ycentr + radius * cos ( angle )
  end do
!
!  Draw the circle.
!
  if ( filled ) then
    call plygon ( npts-1, xpoint, ypoint )
  else
    call plylin ( npts, xpoint, ypoint )
  end if

  return
end
subroutine clrrct ( xmin, ymin, xmax, ymax )

!*****************************************************************************80
!
!! CLRRCT clears a rectangular area of the display.
!
!  Discussion:
!
!    The area is cleared by filling it with the current background color.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real XMIN, YMIN, XMAX, YMAX, the coordinates of the lower
!    left and upper right corners of the area to be cleared.
!
  implicit none
!
  integer icolor
  integer npts
  real x(4)
  real xmax
  real xmin
  real y(4)
  real ymax
  real ymin
!
  x(1) = xmin
  x(2) = xmax
  x(3) = xmax
  x(4) = xmin
  y(1) = ymin
  y(2) = ymin
  y(3) = ymax
  y(4) = ymax
  icolor = 0
  call filclr ( icolor )
  npts = 4
  call plygon ( npts, x, y )

  return
end
subroutine device ( devnam )

!*****************************************************************************80
!
!! DEVICE sets the output device.
!
!  Discussion:
!
!    The setting is made immediately;  an appropriate output file name
!    is saved for later use by GRFINI.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, character ( len = * ) DEVNAM, the name of the device.  Legal
!    values include:
!
!    'cgmb', CGM binary file;
!    'cgmc', CGM text file;
!    'gks', GKS;
!    'igl', ???;
!    'peritek', Peritek;
!    'ps', PostScript;
!    'tek', Tektronix terminal;
!    'tek4207', Tektronix 4207 color;
!    'uis', ?;
!    'x11', X windows;
!    'xws', X windows.
!
  implicit none

  logical devflg
  character ( len = * ) devnam
  character ( len = 128 ) filnam
  logical fnmflg
  integer ierr
  logical iniflg
  logical psflg

  common /drwchr/ filnam
  common /drwstt/ iniflg, psflg, devflg, fnmflg

  save /drwchr/
  save /drwstt/
!
!  If the device has already been set, ignore this call.
!
  if ( devflg ) then

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DEVICE - Warning!'
    write ( *, '(a)' ) '  The graphical device has already been set.'
    write ( *, '(a)' ) '  This call will be ignored.'
!
!  Otherwise, set the device.
!
  else

    ierr = 0
    call setdev ( devnam, ierr )

    if ( ierr /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DEVICE - Error!'
      write ( *, '(a,i6)' ) '  SETDEV returned IERR = ', ierr
      return 
    end if

    devflg = .true.

    if ( .not. fnmflg ) then
      if ( devnam(1:3) == 'tek' .or. &
           devnam(1:1) == 'x' .or. &
           devnam(1:3) == 'uis' .or. &
           devnam(1:3) == 'gks' .or. &
           devnam(1:3) == 'igl' .or. &
           devnam(1:7) == 'peritek' ) then
        filnam = '-'
      else
        filnam = 'cgmout'
      end if
    end if
  end if

  return
end
subroutine drawit ( ipixel, nxdim, nydim, xmin, ymin, xmax, ymax )

!*****************************************************************************80
!
!! DRAWIT displays a cell array.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer IPIXEL(NXDIM,NYDIM), the cell array.
!
!    Input, integer NXDIM, NYDIM, the dimensions of IPIXEL.
!
!    Input, real XMIN, YMIN, the device coordinates where the lower
!    left corner of the cell array should appear.
!
!    Input, real XMAX, YMAX, the device coordinates where the upper
!    right corner of the cell array should appear.
!
  implicit none
!
  integer nxdim
  integer nydim
!
  integer ierr
  integer ipixel(nxdim,nydim)
  real xmax
  real xmin
  real ymax
  real ymin
!
  ierr = 0

  call wrtcla ( ipixel, nxdim, nydim, xmin, ymin, xmax, ymax, &
    xmax, ymin, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DRAWIT - Error!'
    write ( *, '(a,i6)' ) '  WRTCLA returned IERR = ', ierr
  end if

  return
end
subroutine drwpix ( istat, xorg, yorg, xmax, ymax, ipixel, nxdim, nydim, &
  itable )

!*****************************************************************************80
!
!! DRWPIX displays a cell array, and does housekeeping.
!
!  Discussion:
!
!    On the first call, DRWPIX opens the graphics package.
!    On intermediate calls, DRWPIX clears the screen.
!    On the final call, DRWPIX closes the graphics package.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer ISTAT:
!    0 for first call;
!    1 for subsequent;
!    2 for last call.
!
!    Input, real XORG, YORG, the device coordinates where the lower
!    left corner of the cell array should appear.
!
!    Input, real XMAX, YMAX, the device coordinates where the upper
!    right corner of the cell array should appear.
!
!    Input, integer IPIXEL(NXDIM,NYDIM), the cell array.
!
!    Input, integer NXDIM, NYDIM, the dimensions of IPIXEL.
!
!    Input, integer ITABLE, the index of a color table.
!
  implicit none

  integer nxdim
  integer nydim

  logical devflg
  logical fnmflg
  logical iniflg
  integer ipixel(nxdim,nydim)
  integer istat
  integer itable
  logical psflg
  real xmax
  real xorg
  real ymax
  real yorg

  common /drwstt/ iniflg, psflg, devflg, fnmflg

  save /drwstt/
!
!  Set the color table.
!
  call setctb ( itable )

  if ( istat == 0 ) then

    call grfini
    call drawit ( ipixel, nxdim, nydim, xorg, yorg, xmax, ymax )

  else if ( istat == 1 ) then

    if ( iniflg ) then
      call newfrm
    else
      call grfini
    end if

    call drawit ( ipixel, nxdim, nydim, xorg, yorg, xmax, ymax )

  else if ( istat == 2 ) then

    if ( iniflg ) then
      call newfrm
    else
      call grfini
    end if

    call drawit ( ipixel, nxdim, nydim, xorg, yorg, xmax, ymax )
    call grfcls

  end if

  return
end
subroutine filclr ( icolor )

!*****************************************************************************80
!
!! FILCLR sets the color of filled polygons.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer ICOLOR, the index of the color to use for filling
!    polygons.
!
  implicit none

  integer icolor
  integer ierr
!
!  Set the polygon fill color.
!
  ierr = 0
  call wrpgnc ( icolor, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FILCLR - Error!'
    write ( *, '(a,i6)' ) '  WRPGNC returned IERR = ', ierr
  end if

  return
end
subroutine getclr ( icolor, rval, gval, bval, setflg )

!*****************************************************************************80
!
!! GETCLR gets the red, green, and blue values of a given color index.
!
!  Discussion:
!
!    GETCLR is included to insulate the user level
!    from the base level subroutine calls.  setflg is .true. if the
!    particular color index has been set.  if the color index has not
!    been set, setflg returns .false. and rval, gval, and bval are all
!    zero.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer ICOLOR, the index of the color.
!
!    Output, real RVAL, GVAL, BVAL, the RGB values of the color.
!
!    Output, logical SETFLG, is TRUE if the color was set,
!    and FALSE otherwise.
!
  implicit none

  integer, parameter :: mtblsz = 256

  real barray(mtblsz)
  real bval
  logical cstflg(mtblsz)
  real garray(mtblsz)
  real gval
  integer icolor
  real rarray(mtblsz)
  real rval
  logical setflg

  common /clrcom/ rarray, garray, barray, cstflg

  save /clrcom/

  if ( 1 <= icolor+1 .and. icolor+1 <= mtblsz .and. cstflg(icolor+1) ) then
    rval = rarray(icolor+1)
    gval = garray(icolor+1)
    bval = barray(icolor+1)
    setflg = .TRUE.
  else
    rval = 0.0E+00
    gval = 0.0E+00
    bval = 0.0E+00
    setflg = .FALSE.
  end if

  return
end
subroutine getctb ( minc, maxc, fname, ierror )

!*****************************************************************************80
!
!! GETCTB reads a color table from a file FNAME.
!
!  Discussion:
!
!    GETCTB maps the color table by linear interpolation into the
!    interval between color indices min and max inclusive.  The color
!    table in the file is assumed to contain all indices between some
!    minimum and maximum value of the color index.  The maximum allowed
!    value of the color indices in the table is mtblsz, set below.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer MINC, MAXC, the minimum and maximum color indices
!    described in the file.
!
!    Input, character ( len = * ) FNAME, the name of the color table file.
!
!    Output, integer IERROR, error flag.
!    0, no error.
!    1, the file could not be opened.
!    2, an error occurred while reading the file.
!
  implicit none

  integer, parameter :: mtblsz = 256

  real barray(mtblsz)
  logical cstflg(mtblsz)
  character ( len = * ) fname
  real garray(mtblsz)
  integer i
  integer ierror
  integer maxc
  integer maxdup
  integer minc
  integer mindup
  real rarray(mtblsz)

  common / clrcom / rarray, garray, barray, cstflg

  save / clrcom /
!
!  Check that the given minimum and maximum are acceptable, and correct if
!  necessary.
!
  mindup = min ( maxc, minc )
  mindup = max ( mindup, 0 )
  mindup = min ( mindup, mtblsz-1 )

  maxdup = max ( maxc, minc )
  maxdup = max ( maxdup, 0 )
  maxdup = min ( maxdup, mtblsz-1 )
!
!  Read in the color table.
!
  call rdclis ( rarray, garray, barray, mtblsz, mindup, maxdup, fname, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GETCTB - Warning!'
    write ( *, '(a)' ) '  An error occurred in RDCLIS.'
    return
  end if
!
!  Set the 'color set' flags.
!
  do i = mindup, maxdup
    cstflg(i+1) = .true.
  end do

  ierror = 0

  return
end
subroutine grfcls

!*****************************************************************************80
!
!! GRFCLS is the graphics shutdown routine.
!
!  Discussion:
!
!    GRFCLS should be called after all graphics activity is completed.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    None
!
  implicit none

  logical devflg
  logical fnmflg
  integer ierr
  logical iniflg
  logical psflg

  common /drwstt/ iniflg, psflg, devflg, fnmflg

  save /drwstt/

  ierr = 0
  call wrendp ( ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GRFCLS - Error!'
    write ( *, '(a,i6)' ) '  WRENDP returned IERR = ', ierr
  end if

  if ( psflg ) then
    write ( *, '(a,$)' ) char(7)
    read ( *, * )
  end if

  call wrtend ( ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GRFCLS - Error!'
    write ( *, '(a,i6)' ) '  WRTEND returned IERR = =', ierr
  end if

  iniflg = .false.

  return
end
subroutine grfini

!*****************************************************************************80
!
!! GRFINI is the graphics initializer.
!
!  Discussion:
!
!    GRFINI knows if it has been called before in this run (via the
!    variable iniflg), and will do no harm if called twice.  It will
!    use a device name if it has been set (via device), or the default
!    name cgmout if none has been set.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    None
!
  implicit none

  integer, parameter :: mtblsz = 256

  real barray(mtblsz)
  logical cstflg(mtblsz)
  logical devflg
  character ( len = 128 ) filnam
  logical fnmflg
  real garray(mtblsz)
  integer i
  integer ierr
  logical iniflg
  integer minclr
  integer mxclr
  integer nclrs
  integer nctot
  logical psflg
  real rarray(mtblsz)

  common /clrcom/ rarray, garray, barray, cstflg
  common /drwchr/ filnam
  common /drwstt/ iniflg, psflg, devflg, fnmflg

  save /clrcom/
  save /drwchr/
  save /drwstt/
!
!     if ( iniflg ) return
!
!     psflg = .false.
!     devflg = .false.
!     fnmflg = .false.
!
!     do i = 1, mtblsz
!       cstflg(i) = .false.
!     end do
!
  if ( .not. fnmflg ) then
    filnam = 'cgmout'
    devflg = .true.
  end if

  mxclr = mtblsz-1
  ierr = 0

  call wrtopn ( filnam, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GRFINI - Error!'
    write ( *, '(a,i6)' ) '  WRTOPN returned IERR = ', ierr
  end if

  call wrmxci ( mxclr, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GRFINI - Error!'
    write ( *, '(a,i6)' ) ' WRMXCI returned IERR = ', ierr
  end if

  call wrbegp ( ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GRFINI - Error!'
    write ( *, '(a,i6)' ) '  WRBEGP returned IERR = ', ierr
  end if

  if ( cstflg(1) ) then
    call wrbgdc ( rarray(1), garray(1), barray(1), ierr )
    if ( ierr /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'GRFINI - Error!'
      write ( *, '(a,i6)' ) '  WRBGDC returned IERR = ', ierr
    end if
  end if

  call wrbgpb ( ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GRFINI - Error!'
    write ( *, '(a,i6)' ) '  WRBGPB returned IERR = ', ierr
  end if

  call wristl (1, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GRFINI - Error!'
    write ( *, '(a,i6)' ) '  WRISTL returned IERR = ', ierr
  end if

  minclr = 1
  nclrs = 0
  nctot = 0

  do i = 1, mtblsz

    if ( cstflg(i) ) then

      if ( nclrs  ==  0 ) then
        minclr = i
      end if

      nclrs = nclrs + 1
      nctot = nctot + 1

    else if ( nclrs > 0 ) then

      call wrctbl ( rarray(minclr), garray(minclr), barray(minclr), minclr-1, &
        nclrs, ierr )

      if ( ierr /= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'GRFINI - Error!'
        write ( *, '(a,i6)' ) '  WRCTBL returned IERR = ', ierr
      end if

      nclrs = 0

    end if

  end do

  if ( nclrs > 0 ) then

    call wrctbl ( rarray(minclr), garray(minclr), barray(minclr), minclr-1, &
      nclrs, ierr )

    if ( ierr /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'GRFINI - Error!'
      write ( *, '(a,i6)' ) '  WRCTBL returned IERR = ', ierr
    end if

  end if

  if ( nctot == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GRFINI - Warning'
    write ( *, '(a)' ) '  GRFINI called with no color table set!'
  end if

  iniflg = .true.

  return
end
subroutine horcbr ( xmin, ymin, xmax, ymax, ncmin, ncmax, lstr, rstr, ilbclr, &
  size )

!*****************************************************************************80
!
!! HORCBR draws a horizontal color bar.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real XMIN, YMIN, XMAX, YMAX, the X and Y coordinates of the
!    lower left and upper right coordinates of the color bar.
!
!    Input, integer NCMIN, NCMAX, the minimum and maximum color indices
!    to be used in drawing the color bar.
!
!    Input, character ( len = * ) LSTR, RSTR, strings to be placed at the left
!    and right of the color bar.
!
!    Input, integer ILBCLR, the color index to use for the labels.
!
!    Input, real SIZE, the character size to use for the labels.
!
  implicit none

  integer, parameter :: mtblsz = 256

  integer icolor
  integer ilbclr
  integer j
  character ( len = * ) lstr
  integer maxcdp
  integer mincdp
  integer ncmax
  integer ncmin
  integer ncolrs
  character ( len = * ) rstr
  real size
  real x(5)
  real xleft
  real xmax
  real xmin
  real xright
  real xstep
  real y(5)
  real ymax
  real ymin
!
!  Set up the coordinates.
!
  mincdp = min ( ncmin, ncmax )
  maxcdp = max ( ncmin, ncmax )
  mincdp = max ( mincdp, 0 )
  maxcdp = max ( maxcdp, 0 )
  maxcdp = min ( maxcdp, mtblsz-1 )
  mincdp = min ( mincdp, mtblsz-1 )

  ncolrs = maxcdp - mincdp + 1
  xstep = ( xmax - xmin ) / real ( ncolrs )
!
!  Draw the color bar.
!
  y(1) = ymin
  y(2) = ymin
  y(3) = ymax
  y(4) = ymax
  xleft = xmin
  do j = 1, ncolrs
    icolor = j + mincdp - 1
    call filclr ( icolor )
    xright = xleft + xstep
    x(1) = xleft
    x(2) = xright
    x(3) = xright
    x(4) = xleft
    call plygon ( 4, x, y )
    xleft = xright
  end do
!
!  Draw a box around the color bar.
!
  x(1) = xmin
  x(2) = xmax
  x(3) = xmax
  x(4) = xmin
  x(5) = xmin
  y(1) = ymin
  y(2) = ymin
  y(3) = ymax
  y(4) = ymax
  y(5) = ymin
  call linclr ( 1 )
  call plylin ( 5, x, y )
!
!  Draw labels.
!
  call label ( xmin-0.5*size*len(lstr), ymin-1.5*size, lstr, ilbclr, size )

  call label ( xmax-0.5*size*len(rstr), ymin-1.5*size, rstr, ilbclr, size )

  return
end
subroutine horflp ( ipixel, nxdim, nydim )

!*****************************************************************************80
!
!! HORFLP inverts the image left-to-right, in place.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input/output, integer IPIXEL(NXDIM,NYDIM), the cell array.
!
!    Input, integer NXDIM, NYDIM, the dimensions of IPIXEL.
!
  implicit none

  integer nxdim
  integer nydim

  integer ileft
  integer ipixel(nxdim,nydim)
  integer iright
  integer j

  iright = nxdim + 1

  do ileft = 1, nxdim/2
    iright = iright - 1
    do j = 1, nydim 
      call i4_swap ( ipixel(ileft,j), ipixel(iright,j) )
    end do
  end do

  return
end
subroutine i4_swap ( i, j )

!*****************************************************************************80
!
!! I4_SWAP switches two integer values.
!
!  Modified:
!
!    30 November 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer I, J.  On output, the values of I and
!    J have been interchanged.
!
  integer i
  integer j
  integer k

  k = i
  i = j
  j = k

  return
end
subroutine imgmsk ( mfield, nxsmll, nysmll, ipixel, nxbig, nybig )

!*****************************************************************************80
!
!! IMGMSK masks out regions of a pixel array.
!
!  Discussion:
!
!    Regions of the mask array which are nonzero are set to the same 
!    value in the larger array.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer MFIELD(NXSMLL,NYSMLL), the mask array.
!
!    Input, integer NXSMLL, NYSMLL, the dimensions of the mask array.
!
!    Input/output, integer IPIXEL(NXBIG,NYBIG), the cell array.
!
!    Input, integer NXBIG, NYBIG, the dimensions of IPIXEL.
!
  implicit none

  integer nxsmll
  integer nysmll
  integer nxbig
  integer nybig

  integer i
  integer ipixel(nxbig,nybig)
  integer ismall
  integer j
  integer jsmall
  integer mfield(nxsmll,nysmll)
  integer nxbgst
  integer nybgst

  nxbgst = nxbig-1
  nybgst = nybig-1

  do i = 1, nxbig-1

    ismall = 1 + ( ( i - 1 ) * nxsmll ) / nxbgst

    do j = 1, nybig-1
      jsmall = 1 + ((j-1)*nysmll)/nybgst
      if ( mfield(ismall,jsmall) /= 0 ) then
        ipixel(i,j) = mfield(ismall,jsmall)
      end if
    end do

    if ( mfield(ismall,nysmll) /= 0 ) then
      ipixel(i,nybig) = mfield(ismall,nysmll)
    end if

  end do

  do j = 1, nybig-1

    jsmall = 1 + ((j-1)*nysmll)/nybgst

    if ( mfield(nxsmll,jsmall) /= 0 ) then
      ipixel(nxbig,j) = mfield(nxsmll,jsmall)
    end if

  end do

  if ( mfield(nxsmll,nysmll) /= 0 ) then
    ipixel(nxbig,nybig) = mfield(nxsmll,nysmll)
  end if

  return
end
subroutine interp_pixel_array ( imin, imax, ipixel, nxdim, nydim, iframe, nfrms )

!*****************************************************************************80
!
!! INTERP computes a pixel array by interpolation.
!
!  Discussion:
!
!    The interpolated location is iframe-1/nfrms of the way from 
!    imin to imax, so that if iframe = 1 ipixel takes the values of imin,
!    while if iframe=nfrms+1 ipixel is equal to imax.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer IMIN(NXDIM,NYDIM), ?
!
!    Input, integer IMAX(NXDIM,NYDIM), ?
!
!    Output, integer IPIXEL(NXDIM,NYDIM), the cell array.
!
!    Input, integer NXDIM, NYDIM, the dimensions of IPIXEL.
!
!    Input, integer IFRAME, ?
!
!    Input, integer NFRMS, ?
!
  implicit none

  integer nxdim
  integer nydim

  integer i
  integer iframe
  integer imax(nxdim,nydim)
  integer imin(nxdim,nydim)
  integer ipixel(nxdim,nydim)
  integer istep
  integer j
  integer nfrms

  istep = iframe - 1

  do i = 1, nxdim
    do j = 1, nydim
      ipixel(i,j) = imin(i,j) + ( istep * ( imax(i,j) - imin(i,j) ) ) / nfrms
    end do
  end do

  return
end
subroutine label ( xorg, yorg, string, icolor, size )

!*****************************************************************************80
!
!! LABEL draws simple labels.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real XORG, YORG, the lower left hand corner of the beginning
!    of the label.
!
!    Input, character ( len = * ) STRING, the label.
!
!    Input, integer ICOLOR, the index of the color to use for the label.
!
!    Input, real SIZE, the character size to use for the label.
!
  implicit none

  integer icolor
  integer ierr
  real size
  character ( len = * ) string
  real xorg
  real yorg

  ierr = 0
  call wrtxts ( size, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'LABEL - Error!'
    write ( *, '(a,i6)' ) '  WRTXTS returned IERR = ', ierr
  end if
!
!  Set the text color.
!
  ierr = 0
  call wrtxtc ( icolor, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'LABEL - Error!'
    write ( *, '(a,i6)' ) '  WRTXTC returned IERR = ', ierr
  end if

  ierr = 0
  call wrftxt ( string, xorg, yorg, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'LABEL - Error!'
    write ( *, '(a,i6)' ) '  WRFTXT returned IERR = ', ierr
  end if

  return
end
subroutine linclr ( icolor )

!*****************************************************************************80
!
!! LINCLR sets the color of polylines.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer ICOLOR, the index of the polyline color.
!
  implicit none

  integer icolor
  integer ierr

  ierr = 0
  call wrplnc ( icolor, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'LINCLR - Error!'
    write ( *, '(a,i6)' ) '  WRPLNC returned IERR = ', ierr
  end if

  return
end
subroutine linwid ( width )

!*****************************************************************************80
!
!! LINWID sets the width of polylines.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real WIDTH, the line width to use for polylines.
!
  implicit none

  integer ierr
  real width

  ierr = 0

  call wrplnw ( width, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'LINWID - Error!'
    write ( *, '(a,i6)' ) '  WRPLNC returned IERR = ', ierr
  end if

  return
end
subroutine movcgm ( x, y )

!*****************************************************************************80
!
!! MOVCGM moves to a user-specified point with the "plotting pen" up.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real X, Y, the coordinates of the point.
!
  implicit none

  real x
  real, save, dimension ( 2 ) :: xd = (/ 0, 0 /)
  real y
  real, save, dimension ( 2 ) :: yd = (/ 0, 0 /)

  xd(1) = x
  yd(1) = y

  return
entry drwcgm ( x, y )

!*****************************************************************************80
!
!! DRWCGM resumes drawing from the last point defined by MOVCGM or DRWCGM.
!
!  Discussion:
!
!    second entry point: draw from previous point set by movcgm to point
!    specified in drwcgm with "plotting pen" down:
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real X, Y, the coordinates of the new point.
!
  xd(2) = x
  yd(2) = y
  call plylin ( 2, xd, yd )
  xd(1) = x
  yd(1) = y

  return
end
subroutine mrkclr ( icolor )

!*****************************************************************************80
!
!! MRKCLR sets the color of polymarkers.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer ICOLOR, the index of the color to use for polymarkers.
!
  implicit none

  integer icolor
  integer ierr

  ierr = 0
  call wrpmkc ( icolor, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MRKCLR - Error!'
    write ( *, '(a,i6)' ) '  WRPMKC returned IERR = ', ierr
  end if

  return
end
subroutine mrksiz ( size )

!*****************************************************************************80
!
!! MRKSIZ sets the size of polymarkers.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real SIZE, the size, relative to the default size.
!    Thus SIZE = 2.0 means twice as big.  Note, however, that
!    '.' markers are unaffected by changes in SIZE; they are
!    always drawn as small as possible.
!
  implicit none

  integer ierr
  real size

  ierr = 0
  call wrpmks ( size, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MRKSIZ - Error!'
    write ( *, '(a,i6)' ) '  WRPMKS returned IERR = ', ierr
  end if

  return
end
subroutine mrktyp ( itype )

!*****************************************************************************80
!
!! MRKTYP sets the marker type of polymarkers.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer ITYPE, specifies the marker type:
!    1:'.',
!    2:'+', 
!    3:'*', 
!    4:'o', 
!    5:'x'
!
  implicit none
!
  integer ierr
  integer itype
!
  ierr = 0
  call wrpmkt ( itype, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MRKTYP - Error!'
    write ( *, '(a,i6)' ) '  WRPMKT returned IERR = ', ierr
  end if

  return
end
subroutine newfrm

!*****************************************************************************80
!
!! NEWFRM ends a frame, and begins the next.
!
!  Discussion:
!
!    Color table entries which have been set are written to the new frame
!    automatically.  NEWFRM calls GRFINI if necessary.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    None
!
  implicit none

  integer, parameter :: mtblsz = 256

  real barray(mtblsz)
  logical cstflg(mtblsz)
  logical devflg
  logical fnmflg
  real garray(mtblsz)
  integer i
  integer ierr
  logical iniflg
  integer minclr
  integer nclrs
  logical psflg
  real rarray(mtblsz)

  common /clrcom/ rarray, garray, barray, cstflg
  common /drwstt/ iniflg, psflg, devflg, fnmflg

  save /clrcom/
  save /drwstt/

  if ( .not. iniflg ) then

    call grfini

  else

    ierr = 0
    call wrendp ( ierr )

    if ( ierr /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NEWFRM - Error!'
      write ( *, '(a,i6)' ) '  WRENDP returned IERR = ', ierr
    end if
!
!  If we are pausing in between frames, ring the bell and
!  wait for the user to hit RETURN.
!
    if ( psflg ) then
      write ( *, '(a,$)' ) char(7)
      read ( *, * )
    end if

    call wrbegp ( ierr )

    if ( ierr /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NEWFRM - Error!'
      write ( *, '(a,i6)' ) '  WRBEGP returned IERR = ', ierr
    end if

    if ( cstflg(1) ) then
      call wrbgdc ( rarray(1), garray(1), barray(1), ierr )
      if ( ierr /= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'NEWFRM - Error!'
        write ( *, '(a,i6)' ) '  WRBGDC returned IERR = ', ierr
      end if
    end if

    call wrbgpb ( ierr )

    if ( ierr /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NEWFRM - Error!'
      write ( *, '(a,i6)' ) '  WRBGPB returned IERR = ', ierr
    end if

    call wristl ( 1, ierr )

    if ( ierr /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'NEWFRM - Error!'
      write ( *, '(a,i6)' ) '  WRISTL returned IERR = ', ierr
    end if

    minclr = 1
    nclrs = 0

    do i = 1, mtblsz

      if ( cstflg(i) ) then

        if ( nclrs == 0 ) then
          minclr = i
        end if

        nclrs = nclrs + 1

      else if ( nclrs > 0 ) then

        call wrctbl ( rarray(minclr), garray(minclr), barray(minclr), &
          minclr-1, nclrs, ierr )

        if ( ierr /= 0 ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'NEWFRM - Error!'
          write ( *, '(a,i6)' ) '  WRCTBL returned IERR = ', ierr
        end if

        nclrs = 0

      end if

    end do

    if ( nclrs > 0 ) then

      call wrctbl ( rarray(minclr), garray(minclr), barray(minclr), &
        minclr-1, nclrs, ierr )

      if ( ierr /= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'NEWFRM - Error!'
        write ( *, '(a,i6)' ) '  WRCTBL returned IERR = ', ierr
      end if

    end if

  end if

  return
end
subroutine outfil ( fname )

!*****************************************************************************80
!
!! OUTFIL sets the output graphics file name.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, character ( len = * ) FNAME, the name to use for the output file.
!
  implicit none

  logical devflg
  character ( len = 128 ) filnam
  character ( len = * ) fname
  logical fnmflg
  logical iniflg
  logical psflg

  common /drwchr/ filnam
  common /drwstt/ iniflg, psflg, devflg, fnmflg

  save /drwchr/
  save /drwstt/

  filnam = fname
  fnmflg = .true.

  return
end
function r_pi ( )

!*****************************************************************************80
!
!! R_PI returns the value of pi.
!
!  Modified:
!
!    04 December 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real R_PI, the value of pi.
!
  implicit none

  real r_pi

  r_pi = 3.14159265358979323846264338327950288419716939937510E+00

  return
end
subroutine pltbar ( xmin, ymin, xmax, ymax, minclr, maxclr, irclr, igclr, &
  ibclr )

!*****************************************************************************80
!
!! PLTBAR draws a graph of the color intensities of the color bar.
!
!  Discussion:
!
!    irclr, igclr, and ibclr are used to trace out the color curves for
!    red, green, and blue respectively.  
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real XMIN, YMIN, XMAX, YMAX, ?
!
!    Input, integer MINCLR, MAXCLR, the minimum and maximum indices
!    of colors to be used in the color bar.
!
!    Input, integer IRCLR, IGCLR, IBCLR, ?
!
  implicit none

  integer, parameter :: mtblsz = 256

  real barray(mtblsz)
  real bval
  real garray(mtblsz)
  real gval
  integer i
  integer ibclr
  integer igclr
  integer irclr
  integer maxcdp
  integer maxclr
  integer mincdp
  integer minclr
  real rarray(mtblsz)
  real rval
  logical setflg
  real x(mtblsz)
  real xcmin
  real xmax
  real xmin
  real xsize
  real xstep
  real y(5)
  real ybmax
  real ybmin
  real ygmax
  real ygmin
  real ymax
  real ymin
  real yrmax
  real yrmin
  real ysize
!
!  Set up coordinates
!
  mincdp = min ( minclr, maxclr )
  maxcdp = max ( minclr, maxclr )
  mincdp = max ( mincdp, 0 )
  maxcdp = max ( maxcdp, 0 )
  maxcdp = min ( maxcdp, mtblsz-1 )
  mincdp = min ( mincdp, mtblsz-1 )

  xsize = xmax - xmin
  ysize = ( ymax - ymin ) / 3.0
  ybmin = ymin
  ybmax = ybmin + ysize
  ygmin = ybmax
  ygmax = ygmin + ysize
  yrmin = ygmax
  yrmax = yrmin + ysize
  xcmin = xmin

  if ( maxcdp > mincdp ) then
    xstep = xsize / real ( maxcdp - mincdp )
  else
    xstep = 0.0E+00
  end if
!
!  Get the color table, translating values into y coordinates.
!
  do i = 0, maxcdp - mincdp
    call getclr ( i+mincdp, rval, gval, bval, setflg )
    rarray(i+1) = yrmin + rval * ysize
    garray(i+1) = ygmin + gval * ysize
    barray(i+1) = ybmin + bval * ysize

    x(i+1) = xcmin + real ( i ) * xstep

  end do
!
!  Clear the drawing area.
!
  call clrrct ( xmin, ymin, xmax, ymax )
!
!  Draw the curves for the case of maxclr >minclr, or alternately
!  the degenerate case (maxclr = minclr).
!
  if ( maxcdp > mincdp ) then

    call linclr ( irclr )
    call plylin ( maxcdp-mincdp+1, x, rarray )
    call linclr ( igclr )
    call plylin ( maxcdp-mincdp+1, x, garray )
    call linclr ( ibclr )
    call plylin ( maxcdp-mincdp+1, x, barray )

  else

    rarray(2) = rarray(1)
    garray(2) = garray(1)
    barray(2) = barray(1)
    x(2) = xmax
    call linclr ( irclr )
    call plylin ( 2, x, rarray )
    call linclr ( igclr )
    call plylin ( 2, x, garray )
    call linclr ( ibclr )
    call plylin ( 2, x, barray )

  end if
!
!  Draw a rectangle around the drawing area, and individual rectangles
!  around the curves.
!
  x(1) = xmin
  x(2) = xmax
  x(3) = xmax
  x(4) = xmin
  x(5) = xmin
  y(1) = yrmin
  y(2) = yrmin
  y(3) = yrmax
  y(4) = yrmax
  y(5) = yrmin
  call linclr ( irclr )
  call plylin ( 5, x, y )

  x(1) = xmin
  x(2) = xmax
  x(3) = xmax
  x(4) = xmin
  x(5) = xmin
  y(1) = ygmin
  y(2) = ygmin
  y(3) = ygmax
  y(4) = ygmax
  y(5) = ygmin
  call linclr ( igclr )
  call plylin ( 5, x, y )

  x(1) = xmin
  x(2) = xmax
  x(3) = xmax
  x(4) = xmin
  x(5) = xmin
  y(1) = ybmin
  y(2) = ybmin
  y(3) = ybmax
  y(4) = ybmax
  y(5) = ybmin
  call linclr ( ibclr )
  call plylin ( 5, x, y )

  x(1) = xmin
  x(2) = xmax
  x(3) = xmax
  x(4) = xmin
  x(5) = xmin
  y(1) = ymin
  y(2) = ymin
  y(3) = yrmax
  y(4) = yrmax
  y(5) = ymin
  call linclr ( 1 )
  call plylin ( 5, x, y )

  return
end
subroutine plygon ( npts, x, y )

!*****************************************************************************80
!
!! PLYGON draws polygons.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer NPTS, the number of points defining the polygon.
!
!    Input, real X(NPTS), Y(NPTS), the X and Y coordinates of the
!    points defining the polygon.
!
  implicit none

  integer npts

  integer ierr
  real x(npts)
  real y(npts)

  ierr = 0

  call wrtpgn ( x, y, npts, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PLYGON - Error!'
    write ( *, '(a,i6)' ) '  WRTPGN returned IERR = ', ierr
  end if

  return
end
subroutine plylin ( npts, x, y )

!*****************************************************************************80
!
!! PLYLIN draws polylines.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer NPTS, the number of points.
!
!    Input, real X(NPTS), Y(NPTS), the X and Y coordinates of the
!    points defining the polyline.
!
  implicit none

  integer npts

  integer ierr
  real x(npts)
  real y(npts)

  ierr = 0
  call wrplin ( x, y, npts, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PLYLIN - Error!'
    write ( *, '(a,i6)' ) '  WRPLIN returned IERR = ', ierr
  end if

  return
end
subroutine plymrk ( npts, x, y )

!*****************************************************************************80
!
!! PLYMRK draws polymarkers.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer NPTS, the number of points.
!
!    Input, real X(NPTS), Y(NPTS), the X and Y coordinates where the
!    markers are to be placed.
!
  implicit none

  integer npts

  integer ierr
  real x(npts)
  real y(npts)

  ierr = 0
  call wrtpmk ( x, y, npts, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PLYMRK - Error!'
    write ( *, '(a,i6)' ) '  WRTPMK returned IERR = ', ierr
  end if

  return
end
subroutine putctb ( imin, imax, fname, ierror )

!*****************************************************************************80
!
!! PUTCTB dumps the current working color table.
!
!  Modified:
!
!    07 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer IMIN, IMAX, the minimum and maximum color indices
!    to be written to the file.
!
!    Input, character ( len = * ) FNAME, the name of the file to be created
!    by this call.
!
!    Output, integer IERROR, error flag.
!    0, no error.
!    1, the file could not be opened.
!    2, an error occurred while writing the file.
!
  implicit none

  integer, parameter :: mtblsz = 256

  real barray(mtblsz)
  logical cstflg(mtblsz)
  character ( len = * ) fname
  real garray(mtblsz)
  integer ierror
  integer imax
  integer imin
  real rarray(mtblsz)

  common /clrcom/ rarray, garray, barray, cstflg

  save /clrcom/

  call wrclis ( rarray, garray, barray, mtblsz, imin, imax, fname, ierror )

  return
end
subroutine rdclis ( rarray, garray, barray, nclrs, imin, imax, fname, ierror )

!*****************************************************************************80
!
!! RDCLIS reads a color list in a file.
!
!  Discussion:
!
!    The routine maps the color list by linear interpolation into the 
!    interval between color indices IMIN+1 and IMAX+1 inclusive in the r, g,
!    and b color arrays.  The color table in the file is assumed to contain 
!    all indices between some minimum and maximum value of the color index.
!    The maximum allowed value of the color indices in the table
!    is mtblsz, set below.  
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Output, real RARRAY(NCLRS), GARRAY(NCLRS), BARRAY(NCLRS), the interpolated
!    color arrays.
!
!    Input, integer NCLRS, ?
!
!    Input, integer IMIN, IMAX, ?
!
!    Input, character ( len = * ) FNAME, the name of the file containing the
!    color list information.
!
!    Output, integer IERROR, error flag.
!    0 if everything went ok, 
!    1 if it was impossible to open the file, 
!    2 if it was impossible to read from the file, 
!    3 if IMIN or IMAX would result in writing outside the bounds of 
!      the color array.
!
  implicit none

  integer, parameter :: mtblsz = 256

  integer nclrs

  real barray(nclrs)
  real braw(mtblsz)
  real bval
  character ( len = * ) fname
  real gap
  real garray(nclrs)
  real graw(mtblsz)
  real gval
  integer i
  integer ierror
  integer imax
  integer imin
  real interp
  integer itbl
  integer maxdup
  integer maxtbl
  integer mindup
  integer mintbl
  real rarray(nclrs)
  real rloc
  real rraw(mtblsz)
  real rval
  real stride
  real v1
  real v2
  real x
!
!  Check array bounds
!
  if ( imin < 0 .or. imax >= nclrs .or. imin > imax ) then
    ierror = 3
    return
  end if
!
!  Read in the color table
!
  maxtbl = 0
  mintbl = mtblsz - 1

  open ( unit = 99, file = fname, status = 'old', err = 40 )

10    continue

  read ( 99, *, end = 20, err=50 ) itbl, rval, gval, bval

  if ( 0 <= itbl .and. itbl <= mtblsz-1 ) then

    rraw(itbl+1) = rval
    graw(itbl+1) = gval
    braw(itbl+1) = bval

    maxtbl = max ( maxtbl, itbl )
    mintbl = min ( mintbl, itbl )

  end if

  go to 10
!
!  Map it into the requested interval.
!
20    continue

  close ( unit = 99 )

  mindup = min ( imin, imax )
  mindup = max ( mindup, 0 )
  mindup = min ( mindup, mtblsz-1 )

  maxdup = max ( imin, imax )
  maxdup = max ( maxdup, 0 )
  maxdup = min ( maxdup, mtblsz-1 )

  if ( maxdup > mindup ) then
    stride = real ( maxtbl - mintbl ) / real ( maxdup - mindup )
  else
    stride = 0.0
  end if

  rloc = mintbl
!
!  There is something wrong with this code.
!  In its original form, the indexing of the RAW arrays goes out
!  of bounds at the end.
!  And why does the first indexing begin at location 2 rather than 1?
!
  do i = mindup+1, maxdup+1
    itbl = min ( int ( rloc ), maxtbl - 2 )
    gap = rloc - itbl
    rarray(i) = interp ( rraw(itbl+1), rraw(itbl+2), gap )
    garray(i) = interp ( graw(itbl+1), graw(itbl+2), gap )
    barray(i) = interp ( braw(itbl+1), braw(itbl+2), gap )
    rloc = rloc + stride
  end do

  ierror = 0
  return
!
!  Error handling:  open error, then read error.
!
40 continue
  ierror = 1
  return

50 continue
  ierror = 2
  return
end
function interp ( v1, v2, x ) 

!*****************************************************************************80
!
!! INTERP linearly interpolates between two values.
!
!  Modified:
!
!    09 June 2000
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real V1, V2, the two values to interpolate.
!
!    Input, real X, the weight assigned to V2.  V1 will be assigned
!    the weight 1-X.  Normally, X will be between 0 and 1.
!
!    Output, real INTERP, the interpolated value.
!
  implicit none

  real interp
  real v1
  real v2
  real x

  interp = ( 1.0E+00 - x ) * v1 + x * v2

return
end
subroutine rtoint ( rarray, iarray, nxdim, nydim, rmn, rmx, imn, imx )

!*****************************************************************************80
!
!! RTOINT converts a real array to integers, given ranges.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real RARRAY(NXDIM,NYDIM), ?
!
!    Output, integer IARRAY(NXDIM,NYDIM), ?
!
!    Input, integer NXDIM, NYDIM, the dimensions of RARRAY and IARRAY.
!
!    Input, real RMN, RMX, ?
!
!    Input, integer IMN, IMX, ?
!
  implicit none

  integer nxdim
  integer nydim

  real bndhi
  real bndlow
  logical errflg
  integer i
  integer iarray(nxdim,nydim)
  integer imn
  integer imx
  integer irange
  integer j
  real rarray(nxdim,nydim)
  real rmn
  real rmx
  real rrange

  errflg = .false.

  bndlow = min ( rmn, rmx )
  bndhi = max ( rmn, rmx )

  rrange = rmx - rmn

  if ( rrange == 0.0E+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'RTOINT - Error!'
    write ( *, '(a)' ) '  RMX = RMN'
    return
  end if

  irange = imx - imn

  do i = 1, nxdim

    do j = 1, nydim

      if ( rarray(i,j) < bndlow .or. rarray(i,j) > bndhi ) then

        iarray(i,j) = 1

        if ( .not. errflg ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'RTOINT - Warning!'
          write ( *, '(a)' ) '  Found a value outside the range '
          errflg = .true.
        end if

      else

        iarray(i,j) = imn + int ( irange * ( rarray(i,j) - rmn ) / rrange )

      end if

    end do

  end do

  return
end
subroutine setclr ( icolor, rval, gval, bval )

!*****************************************************************************80
!
!! SETCLR sets the red, green and blue values of a given color.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer ICOLOR, the index of the color, between 0 and 255.
!
!    Input, real RVAL, GVAL, BVAL, the red, green, and blue values
!    of the color, between 0.0 and 1.0.
!
  implicit none
!
  integer, parameter :: mtblsz = 256

  real barray(mtblsz)
  real bval
  logical cstflg(mtblsz)
  real garray(mtblsz)
  real gval
  integer icolor
  real rarray(mtblsz)
  real rval

  common /clrcom/ rarray, garray, barray, cstflg

  save /clrcom/

  if ( icolor >= 0 .and. icolor < mtblsz ) then
    rarray(icolor+1) = min ( max ( rval, 0.0E+00 ), 1.0E+00 )
    garray(icolor+1) = min ( max ( gval, 0.0E+00 ), 1.0E+00 )
    barray(icolor+1) = min ( max ( bval, 0.0E+00 ), 1.0E+00 )
    cstflg(icolor+1) = .TRUE.
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SETCLR - Warning!'
    write ( *, '(a,i6)' ) '  Color index out of range: ', icolor
  end if

  return
end
subroutine setctb ( index )

!*****************************************************************************80
!
!! SETCTB calculates a full set of color table entries.
!
!  Discussion:
!
!    An out-of-range index causes no color table to be generated;  i.e.
!    the previous colors will remain in effect.  color indices 0 and 1
!    are set to white and black respectively, if the index is in the
!    above range.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, integer INDEX, selects the color table.
!    1, gray scale.
!    2, blue fading into yellow.
!    3, waves of green in red fading to blue.
!    4, pseudospectral (red to blue).
!    5, inverted pseudospectral (blue to red).
!
  implicit none

  integer, parameter :: mtblsz = 256

  real bval
  real gauss
  real gval
  integer i
  integer index
  real r_pi
  real ratio
  real rval
  real theta
  real tval

  if ( index >= 1 .and. index <= 5 ) then

    do i = 0, mtblsz - 1
!
!  gray scale
!
      if ( index == 1 ) then
        rval = real ( i ) / real ( mtblsz - 1 )
        gval = rval
        bval = rval
      end if
!
!  blue into yellow
!
      if ( index == 2 ) then
        tval = real ( mtblsz - 1 )
        rval = real ( i ) / tval
        gval = rval
        bval = 1.0E+00 - rval + 1.0E+00 / tval
      end if
!
!  waves
!
      if ( index == 3 ) then
        theta = 0.5E+00 * r_pi ( ) * real ( i ) / real ( mtblsz - 1 )
        rval = cos ( theta )**2
        bval = sin ( theta )**2
        gval = 0.8E+00 * sin ( 10.0E+00 * theta )**6
      end if
!
!  pseudospectral
!
      if ( index == 4 ) then

        theta = 4.0E+00 * real ( i ) / real ( mtblsz - 1 )
        tval = gauss ( theta - 4.0E+00 )
        rval = gauss ( theta - 1.0E+00 ) + tval
        gval = gauss ( theta - 2.0E+00 ) + tval
        bval = gauss ( theta - 3.0E+00 ) + tval

        rval = min ( rval, 1.0E+00 )
        gval = min ( gval, 1.0E+00 )
        bval = min ( bval, 1.0E+00 )

      end if
!
!  inverted pseudospectral
!
      if ( index == 5 ) then
        theta = 4.0E+00 * real ( mtblsz - i - 1 ) / real ( mtblsz - 1 )
        tval = gauss ( theta - 4.0E+00 )
        rval = gauss ( theta - 1.0E+00 ) + tval
        gval = gauss ( theta - 2.0E+00 ) + tval
        bval = gauss ( theta - 3.0E+00 ) + tval
        rval = min ( rval, 1.0E+00 )
        gval = min ( gval, 1.0E+00 )
        bval = min ( bval, 1.0E+00 )
      end if
!
!  Background is to be white
!
      if ( i == 0 ) then
        rval = 1.0E+00
        gval = 1.0E+00
        bval = 1.0E+00
      end if
!
!  color 1 (foreground) is black
!
      if ( i == 1 ) then
        rval = 0.0E+00
        gval = 0.0E+00
        bval = 0.0E+00
      end if

      call setclr ( i, rval, gval, bval )

    end do

  end if

  return
end
function gauss ( ratio ) 

!*****************************************************************************80
!
!! GAUSS ...
!
!  Modified:
!
!    09 June 2000
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real RATIO, ...
!
!    Output, real GAUSS, ...
!
  implicit none

  real gauss
  real ratio

  gauss = exp ( - ratio**2 )

return
end
subroutine setscl ( xcoord, ycoord, npts )

!*****************************************************************************80
!
!! SETSCL sets the coordinate system for drawing.
!
!  Discussion:
!
!    The given input arrays specify the x and y locations of a set of points;
!    the drawing coordinate system is set so that the aspect ratio
!    of the coordinates is preserved but the set of points is framed
!    within the display boundaries.  the most widely separated
!    points in the most widely separated direction will be placed
!    5 percent of the device width from the window sides.  no
!    rotation is used.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real XCOORD(NPTS), YCOORD(NPTS), ?
!
!    Input, integer NPTS, the number of points.
!
  implicit none

  integer npts

  integer i
  integer ierr
  real, parameter :: margin = 0.525E+00
  real xave
  real xcoord(npts)
  real xmax
  real xmin
  real xrange
  real yave
  real ycoord(npts)
  real ymax
  real ymin
  real yrange
!
!  Find the extremes in each direction.
!
  xmin = xcoord(1)
  xmax = xcoord(1)
  ymin = ycoord(1)
  ymax = ycoord(1)
  do i = 2, npts
    xmin = min ( xmin, xcoord(i) )
    xmax = max ( xmax, xcoord(i) )
    ymin = min ( ymin, ycoord(i) )
    ymax = max ( ymax, ycoord(i) )
  end do

  xrange = xmax - xmin
  yrange = ymax - ymin
  xave = xmin + 0.5E+00 * xrange
  yave = ymin + 0.5E+00 * yrange
!
!  If the coordinate ranges are zero in both directions, exit
!  with an error message.
!
  if ( xrange == 0.0E+00 .and. yrange == 0.0E+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SETSCL - Error!'
    write ( *, '(a)' ) '  Degenerate coordinate ranges.'
    write ( *, '(a,g14.6)' ) '  XRANGE = ',xrange
    write ( *, '(a,g14.6)' ) '  YRANGE = ',yrange
    write ( *, '(a)' ) '  The coordinates were not rescaled.'
    return
  end if
!
!  Two cases:  if the span in the X direction is greater, set
!  coordinates so that direction is optimally framed.   The
!  factor of margin provides a 5 percent margin around the given
!  points.
!
  ierr = 0
  if ( xrange >= yrange ) then

    call setwcd ( xave-margin*xrange, yave-margin*xrange, &
      xave+margin*xrange, yave+margin*xrange, ierr )
!
!  Otherwise, set coordinates so that the Y direction is optimally framed.
!
  else

    call setwcd ( xave-margin*yrange, yave-margin*yrange, &
      xave+margin*yrange, yave+margin*yrange, ierr )

  end if

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SETSCL - Error!'
    write ( *, '(a,i6)' ) '  SETWCD returned IERR = ', ierr
  end if
 
  return
end
subroutine splie2 ( x2a, ya, m, n, y2a )

!*****************************************************************************80
!
!! SPLIE2 is a utility routine for spline calculations.
!
!
!  Modified:
!
!    06 March 1999
!
!  Parameters:
!
  implicit none

  integer, parameter :: nn = 1200

  integer n

  integer j
  integer k
  integer m
  real x2a(n)
  real y2a(m,n)
  real y2tmp(nn)
  real ya(m,n)
  real yp1
  real ypn
  real ytmp(nn)

  do j = 1, m

    ytmp(1:n) = ya(j,1:n)

    yp1 = 1.0E+30
    ypn = 1.0E+30
    call spline ( x2a, ytmp, n, yp1, ypn, y2tmp )

    y2a(j,1:n) = y2tmp(1:n)

  end do

  return
end
subroutine spline ( x, y, n, yp1, ypn, y2 )

!*****************************************************************************80
!
!! SPLINE computes a cubic spline through given data.
!
!  comes from numerical recipes by press, flannery,
!       teukolsky, and vetterling, chapter 3.
!
!  Modified:
!
!    06 March 1999
!
!  Parameters:
!
  implicit none
!
  integer, parameter :: nmax = 1200

  integer n

  integer i
  integer k
  real p
  real qn
  real sig
  real u(nmax)
  real un
  real x(n)
  real y(n)
  real y2(n)
  real yp1
  real ypn

  if ( yp1 > 0.99E+30 ) then
    y2(1) = 0.0E+00
    u(1) = 0.0E+00
  else
    y2(1) = - 0.5E+00
    u(1) = ( 3.0E+00 / ( x(2) - x(1) ) ) * &
      ( ( y(2) - y(1) ) / ( x(2) - x(1) ) - yp1 )
  end if

  do i = 2, n-1
    sig = ( x(i) - x(i-1) ) / ( x(i+1) - x(i-1) )
    p = sig * y2(i-1) + 2.0E+00
    y2(i) = ( sig - 1.0E+00 ) / p
    u(i) = ( 6.0E+00 * ( ( y(i+1) - y(i) ) / ( x(i+1) - x(i) ) &
      - ( y(i) - y(i-1) ) / ( x(i) - x(i-1) ) ) & 
      / ( x(i+1) - x(i-1) ) - sig * u(i-1) ) / p
  end do

  if ( ypn > 0.99E+30 ) then
    qn = 0.0E+00
    un = 0.0E+00
  else
    qn = 0.5E+00
    un = ( 3.0E+00 / ( x(n) - x(n-1) ) ) * ( ypn - ( y(n) - y(n-1) ) &
      / ( x(n) - x(n-1) ) )
  end if

  y2(n) = ( un - qn * u(n-1) ) / ( qn * y2(n-1) + 1.0E+00 )

  do k = n-1, 1, -1
    y2(k) = y2(k) * y2(k+1) + u(k)
  end do

  return
end
subroutine stpaus

!*****************************************************************************80
!
!! STPAUS sets the 'pause' flag.
!
!  Discussion:
!
!    STPAUS is called so that the user will be asked to hit return at 
!    each end of frame.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    None
!
  implicit none

  logical devflg
  logical fnmflg
  logical iniflg
  logical psflg

  common /drwstt/ iniflg, psflg, devflg, fnmflg

  save /drwstt/

  psflg = .true.

  return
end
subroutine stpspl ( xa, ya, y2a, n, x, y, indx )

!*****************************************************************************80
!
!! STPSPL does spline interpolation for ordered points.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    ?,  real XA(N), YA(N), ?
!
!    ?, real Y2A(N), ?
!
!    Input, integer N, the dimension of XA, YA and Y2A.
!
!    Input, real X, the abscissa at which interpolation is needed.
!
!    Output, real Y, the value of the spline at X.
!
!    Input/output, integer INDX, the index in XA of the left hand
!    limit of the interval containing X.  On input, this may be set
!    to a suggested value.
!
  implicit none

  integer n

  real a
  real b
  real h
  integer indx
  integer khi
  integer klo
  real x
  real xa(n)
  real y
  real y2a(n)
  real ya(n)
!
!  Make sure INDX is legal.
!
  if ( indx < 1 .or. indx > n-1 ) then
    indx = 1
  end if

5     continue

  if ( xa(indx) > x .or. xa(indx+1) < x ) then
    indx = indx+1
    if ( indx == n ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'STPSPL - Error!'
      write ( *, '(a)' ) '  Bad XA or INDX input.'
      call exit(2)
    end if
    go to 5
  end if

  klo = indx
  khi = indx+1
  h = xa(khi) - xa(klo)

  if ( h == 0.0E+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'STPSPL - Error!'
    write ( *, '(a)' ) '  Bad XA input.'
    call exit(2)
  end if

  a = ( xa(khi) - x ) / h
  b = ( x - xa(klo) ) / h

  y = a * ya(klo) + b * ya(khi) &
    + ( ( a**3 - 1.0E+00 ) * y2a(klo) &
    +   ( b**3 - b ) * y2a(khi) ) * ( h**2 ) / 6.0E+00

  return
end
subroutine rmat_expand_linear ( small, nsx, nsy, big, nx, ny )

!*****************************************************************************80
!
!! RMAT_EXPAND_LINEAR expands a real array by linear interpolation.
!
!  Modified:
!
!    02 May 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real SMALL(NSX,NSY), a "small" array.
!
!    Input, integer NSX, NSY, the number of rows and columns in SMALL.
!
!    Output, real BIG(NX,NY), the big array, which contains an 
!    interpolated version of the data in SMALL.
!
!    Input, integer NX, NY, the number of rows and columns in BIG.
!
  implicit none

  integer nsx
  integer nsy
  integer nx
  integer ny

  real big(nx,ny)
  integer i
  integer i1
  integer i2
  integer j
  integer j1
  integer j2
  real r
  real r1
  real r2
  real s
  real s1
  real s2
  real small(nsx,nsy)

  do i = 1, nx

    if ( nx == 1 ) then
      r = 0.5E+00
    else
      r = real ( i - 1 ) / real ( nx - 1 )
    end if

    i1 = 1 + int ( r * ( nsx - 1 ) ) 
    i2 = i1 + 1
    if ( i2 > nsx ) then
      i1 = nsx - 1
      i2 = nsx
    end if

    r1 = real ( i1 - 1 ) / real ( nsx - 1 )
    r2 = real ( i2 - 1 ) / real ( nsx - 1 )

    do j = 1, ny

      if ( ny == 1 ) then
        s = 0.5E+00
      else
        s = real ( j - 1 ) / real ( ny - 1 )
      end if

      j1 = 1 + int ( s * real ( nsy - 1 ) ) 
      j2 = j1 + 1
      if ( j2 > nsy ) then
        j1 = nsy - 1
        j2 = nsy
      end if

      s1 = real ( j1 - 1 ) / real ( nsy - 1 )
      s2 = real ( j2 - 1 ) / real ( nsy - 1 )

      big(i,j) = &
        ( ( r2 - r ) * ( s2 - s ) * small(i1,j1) &
        + ( r - r1 ) * ( s2 - s ) * small(i2,j1) &
        + ( r2 - r ) * ( s - s1 ) * small(i1,j2) &
        + ( r - r1 ) * ( s - s1 ) * small(i2,j2) ) &
        / ( ( r2 - r1 ) * ( s2 - s1 ) )

    end do

  end do

  return
end
subroutine strspl ( small, y2a, nsx, nsy, big, nx, ny )

!*****************************************************************************80
!
!! STRSPL expands a real array using bicubic natural spline interpolation.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real SMALL(NSX,NSY), a "small" array.
!
!    Input, integer NSX, NSY, the number of rows and columns in SMALL.
!
!    Output, real BIG(NX,NY), the big array, which contains an 
!    interpolated version of the data in SMALL.
!
!    Input, integer NX, NY, the number of rows and columns in BIG.
!
  implicit none

  integer, parameter :: nmax = 1200

  integer nsx
  integer nsy
  integer nx
  integer ny

  real big(nx,ny)
  real f2temp(nmax)
  real fftemp(nmax)
  real ftemp(nmax)
  integer i
  integer indx1
  integer indx2
  integer j
  integer k
  integer l
  real small(nsx,nsy)
  real val
  real x
  real xscale
  real xtemp(nmax)
  real y
  real yp1
  real ypn
  real yscale
  real ytemp(nmax)
  real y2a(nsx,nsy)
!
  xscale = real ( nx - 1 ) / real ( nsx - 1 )
  yscale = real ( ny - 1 ) / real ( nsy - 1 )

  do i = 1, nsx
    xtemp(i) = real ( i - 1 ) * xscale
  end do

  do j = 1, nsy
    ytemp(j) = real ( j - 1 ) * yscale
  end do

  call splie2 ( ytemp, small, nsx, nsy, y2a )

  do j = 1, ny

    y = real ( j - 1 )
    indx1 = 1

    do k = 1, nsx

      ftemp(1:nsy) = small(k,1:nsy)
      f2temp(1:nsy) = y2a(k,1:nsy)

      call stpspl ( ytemp, ftemp, f2temp, nsy, y, fftemp(k), indx1 )

    end do

    yp1 = 1.0E+30
    ypn = 1.0E+30
    call spline ( xtemp, fftemp, nsx, yp1, ypn, f2temp )
    indx2 = 1

    do i = 1, nx
      x = real ( i - 1 )
      call stpspl ( xtemp, fftemp, f2temp, nsx, x, val, indx2 )
      big(i,j) = val
    end do

  end do

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! 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 unpaus

!*****************************************************************************80
!
!! UNPAUS clears the 'pause' flag.
!
!  Discussion:
!
!    Once UNPAUS is called, the user will be not asked to hit return 
!    at the end of each graphics frame.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    None
!
  implicit none

  logical devflg
  logical fnmflg
  logical iniflg
  logical psflg

  common /drwstt/ iniflg, psflg, devflg, fnmflg

  save /drwstt/

  psflg = .false.

  return
end
subroutine vrtcbr ( xmin, ymin, xmax, ymax, ncmin, ncmax, bstr, tstr, ilbclr, &
  size )

!*****************************************************************************80
!
!! VRTCBR draws a vertical color bar.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real XMIN, YMIN, XMAX, YMAX, the X and Y coordinates of the
!    lower left and upper right coordinates of the color bar.
!
!    Input, integer NCMIN, NCMAX, the minimum and maximum color indices
!    to be used in drawing the color bar.
!
!    Input, character ( len = * ) BSTR, TSTR, strings to be placed at the left
!    and right of the color bar.
!
!    Input, integer ILBCLR, the color index to use for the labels.
!
!    Input, real SIZE, the character size to use for the labels.
!
  implicit none

  integer, parameter :: mtblsz = 256

  character ( len = * ) bstr
  integer icolor
  integer ilbclr
  integer j
  integer maxcdp
  integer mincdp
  integer ncmax
  integer ncmin
  integer ncolrs
  real size
  character ( len = * ) tstr
  real x(5)
  real xmax
  real xmin
  real y(5)
  real ybot
  real ymax
  real ymin
  real ystep
  real ytop
!
!  Set up coordinates.
!
  mincdp = min ( ncmin, ncmax )
  maxcdp = max ( ncmin, ncmax )
  mincdp = max ( mincdp, 0 )
  maxcdp = max ( maxcdp, 0 )
  maxcdp = min ( maxcdp, mtblsz-1 )
  mincdp = min ( mincdp, mtblsz-1 )

  ncolrs = maxcdp - mincdp + 1
  ystep = ( ymax - ymin ) / ncolrs
!
!  Draw color bar
!
  x(1) = xmin
  x(2) = xmax
  x(3) = xmax
  x(4) = xmin
  ybot = ymin

  do j = 1, ncolrs
    icolor = j + mincdp - 1
    call filclr ( icolor )
    ytop = ybot + ystep
    y(1) = ybot
    y(2) = ybot
    y(3) = ytop
    y(4) = ytop
    call plygon ( 4, x, y )
    ybot = ytop
  end do
!
!  Draw box around color bar
!
  x(1) = xmin
  x(2) = xmax
  x(3) = xmax
  x(4) = xmin
  x(5) = xmin
  y(1) = ymin
  y(2) = ymin
  y(3) = ymax
  y(4) = ymax
  y(5) = ymin
  icolor = 1
  call linclr ( icolor )
  call plylin ( 5, x, y )
!
!  Draw labels
!
  call label ( xmax+(xmax-xmin)/5, ymin-1.0/120.0, bstr, ilbclr, size )

  call label ( xmax+(xmax-xmin)/5, ymax-1.0/120.0, tstr, ilbclr, size )

  return
end
subroutine vrtflp ( ipixel, nxdim, nydim )

!*****************************************************************************80
!
!! VRTFLP inverts a cell array image top-to-bottom, in place.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input/output, integer IPIXEL(NXDIM,NYDIM), the cell array.
!
!    Input, integer NXDIM, NYDIM, the dimensions of IPIXEL.
!
  implicit none

  integer nxdim
  integer nydim

  integer i
  integer ipixel(nxdim,nydim)
  integer jbot
  integer jtop

  jbot = nydim + 1

  do jtop = 1, nydim / 2
    jbot = jbot - 1
    do i = 1, nxdim
      call i4_swap ( ipixel(i,jtop), ipixel(i,jbot) )
    end do
  end do

  return
end
subroutine winfrm

!*****************************************************************************80
!
!! WINFRM draws a window frame in the (full-screen) display window.
!
!  Discussion:
!
!    The window matches the window of our animation system.
!
!  Modified:
!
!    06 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    None
!
  implicit none

  integer, parameter :: nmax = 5

  integer n
  real x(nmax)
  real y(nmax)
  real ybot
  real ytop

  ybot = 0.1E+00
  ytop = 1.0E+00
  x(1) = 0.0E+00
  x(2) = 1.0E+00
  x(3) = 1.0E+00
  x(4) = 0.0E+00
  x(5) = 0.0E+00
  y(1) = ybot
  y(2) = ybot
  y(3) = ytop
  y(4) = ytop
  y(5) = ybot
  n = 5
  call plylin ( n, x, y )

  x(1) = 0.0E+00
  x(2) = 1.0E+00
  y(1) = ( ytop + ybot ) / 2.0E+00
  y(2) = y(1)
  n = 2
  call plylin ( n, x, y )

  x(1) = 0.5E+00
  x(2) = 0.5E+00
  y(1) = ytop
  y(2) = ybot
  n = 2
  call plylin ( n, x, y )

  return
end
subroutine wrclis ( rarray, garray, barray, nclrs, minc, maxc, fname, ierror )

!*****************************************************************************80
!
!! WRCLIS writes color data to a file.
!
!  Modified:
!
!    07 March 1999
!
!  Author:
!
!    Joel Welling, Pittsburgh Supercomputing Center
!
!  Parameters:
!
!    Input, real RARRAY(NCLRS), GARRAY(NCLRS), BARRAY(NCLRS), the
!    R, G and B values for a set of colors.
!
!    Input, integer NCLRS, the number of colors.
!
!    Input, integer MINC, MAXC, the first and last indices of the
!    colors to be written out.
!
!    Input, character ( len = * ) FNAME, the name of the file.
!
!    Output, integer IERROR, error flag.
!    0, no error.
!    1, the file could not be opened.
!    2, an error occurred while writing the file.
!
  implicit none

  integer, parameter :: mtblsz = 256

  integer nclrs

  real barray(nclrs)
  character ( len = * ) fname
  real garray(nclrs)
  integer i
  integer ierror
  integer maxc
  integer maxdup
  integer minc
  integer mindup
  real rarray(nclrs)
!
!  Open a new copy of the file.
!
  open ( unit = 99, file = fname, status = 'replace', err = 20 )
!
!  Set the array bounds.
!
  maxdup = max ( maxc, minc )
  mindup = min ( maxc, minc )
  mindup = max ( mindup, 0 )
  maxdup = max ( maxdup, 0 )
  maxdup = min ( maxdup, mtblsz-1 )
  mindup = min ( mindup, mtblsz-1 )
!
!  Write out the color table
!
  do i = mindup+1, maxdup+1
    write ( 99, *, err = 30 ) i-1, rarray(i), garray(i), barray(i)
  end do

  close ( unit = 99 )
  ierror = 0
  return
!
!  Could not open the file.
!
20    continue

  ierror = 1
  return
!
!  Error while writing to the file.
!
30    continue

  ierror = 2
 
  return
end