program tiler_3d
!
!*******************************************************************************
!
!! TILER_3D illustrates the use of 3D blending.
!
!
!  Discussion:
!
!    This main program works in a 3D rectangular space we can think
!    of as "UVW" space.  We are interested in the space of data values
!    bounded by (umin,vmin,wmin) and (umax,vmax,wmax), and we plan to
!    divide this up, using indices (I,J,K), into ni by nj by nk sub-boxes.
!
!    The code below considers each sub-box indexed by (I,J,K) and determines
!    the values (u0,v0,w0) and (u1,v1,w1) that characterize its corners.
!    The coordinates of this sub-box and the coordinates of the big box
!    are then passed to SUB_BOX_TILER_3D.
!
!    The picture would be a LOT more interesting if the boundary were
!    a bit more wiggly, there were more sub-boxes, and the object in
!    each sub-box had more parts.
!
!  Reference:
!
!    W N Gordon and Charles A Hall,
!    Construction of Curvilinear Coordinate Systems and Application to
!      Mesh Generation,
!    International Journal of Numerical Methods in Engineering,
!    Volume 7, pages 461-477, 1973.
!
!    Joe Thompson, Bharat Soni, Nigel Weatherill,
!    Handbook of Grid Generation,
!    CRC Press,
!    1999.
!
!  Modified:
!
!    01 July 1999
!
!  Author:
!
!    John Burkardt
!
  implicit none
!
  integer i
  integer j
  integer k
  integer, parameter :: ni = 3
  integer, parameter :: nj = 3
  integer, parameter :: nk = 3
  real u0
  real u1
  real, parameter :: umax =   30.0E+00
  real, parameter :: umin =  150.0E+00
  real v0
  real v1
  real, parameter :: vmax =    5.0E+00
  real, parameter :: vmin =    1.0E+00
  real w0
  real w1
  real, parameter :: wmax = 30.0E+00
  real, parameter :: wmin = -30.0E+00
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TILER_3D'
  write ( *, '(a)' ) '  A simple example of transfinite'
  write ( *, '(a)') '  interpolation in 3D.'
!
!  Write the first line of the output file, which is the number
!  of triangles in the output 3D shape.  (In our simple example,
!  each sub-box will contain a tetrahedron.)
!
  open ( unit = 1, file = 'tiler_3d.tri', status = 'replace' )

  write ( 1, * ) 4 * ni * nj * nk
!
!  Consider items with index (I,*,*):
!
  do i = 1, ni

    u0 = ( real ( ni + 1 - i ) * umin + real ( i - 1 ) * umax ) / real ( ni )

    u1 = ( real ( ni - i ) * umin + real ( i ) * umax ) / real ( ni )
!
!  Consider items with index (I,J,*):
!
    do j = 1, nj

      v0 = ( real ( nj + 1 - j ) * vmin + real ( j - 1 ) * vmax ) / real ( nj )

      v1 = ( real ( nj - j ) * vmin + real ( j ) * vmax ) / real ( nj )
!
!  Consider items with index (I,J,K):
!
      do k = 1, nk

        w0 = ( real ( nk + 1 - k ) * wmin + real ( k - 1 ) * wmax ) / real ( nk )

        w1 = ( real ( nk - k ) * wmin + real ( k ) * wmax ) / real ( nk )
!
!  Fill sub-box (I,J,K) with the 3-D "tile".
!         
        call sub_box_tiler_3d ( umin, vmin, wmin, umax, vmax, wmax, &
          u0, v0, w0, u1, v1, w1 )

      end do

    end do

  end do

  close ( unit = 1 )

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

  stop
end
subroutine sub_box_tiler_3d ( umin, vmin, wmin, umax, vmax, wmax, &
  u0, v0, w0, u1, v1, w1 )
!
!*******************************************************************************
!
!! SUB_BOX_TILER_3D "tiles" a 3D sub-box with a given pattern.
!
!
!  Discussion:
!
!    This routine knows the (U,V,W) coordinates of the big box and the
!    sub box, and knows the shape of the object to be place in the sub-box.
!    It uses transfinite interpolation to put the shape in the box.
!    This requires that, for each point of the shape to be mapped, the 
!    (X,Y,Z) coordinates be evaluated for points on the surface of
!    the big box, namely, at 8 corners, 12 edges, and 6 faces.  These
!    values are then blended to give a sensible (X,Y,Z) coordinate for
!    the point belonging to the shape.
!
!  Reference:
!
!    W N Gordon and Charles A Hall,
!    Construction of Curvilinear Coordinate Systems and Application to
!      Mesh Generation,
!    International Journal of Numerical Methods in Engineering,
!    Volume 7, pages 461-477, 1973.
!
!    Joe Thompson, Bharat Soni, Nigel Weatherill,
!    Handbook of Grid Generation,
!    CRC Press,
!    1999.
!
!  Modified:
!
!    01 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real UMIN, VMIN, WMIN, UMAX, VMAX, WMAX, the (U,V,W) coordinates
!    of the two opposite corners of the big box.
!
!    Input, real U0, V0, W0, U1, V1, W1, the (U,V,W) coordinates of the
!    two oppositie corners of the sub-box.
!
  implicit none
!
  integer, parameter :: npoint = 4
!
  integer i
  real r
  real r_tab(npoint)
  real, parameter :: r0 = 0.0E+00
  real, parameter :: r1 = 1.0E+00
  real s
  real s_tab(npoint)
  real, parameter :: s0 = 0.0E+00
  real, parameter :: s1 = 1.0E+00
  real t
  real t_tab(npoint)
  real, parameter :: t0 = 0.0E+00
  real, parameter :: t1 = 1.0E+00
  real u
  real u0
  real u1
  real umax
  real umin
  real v
  real v0
  real v1
  real vmax
  real vmin
  real w
  real w0
  real w1
  real wmax
  real wmin
  real x(npoint)
  real x000
  real x00t
  real x001
  real x010
  real x011
  real x01t
  real x0s0
  real x0s1
  real x0st
  real x100
  real x101
  real x10t
  real x110
  real x111
  real x11t
  real x1s0
  real x1s1
  real x1st
  real xr00
  real xr01
  real xr0t
  real xr10
  real xr11
  real xr1t
  real xrs0
  real xrs1
  real y(npoint)
  real y000
  real y00t
  real y001
  real y010
  real y011
  real y01t
  real y0s0
  real y0s1
  real y0st
  real y100
  real y101
  real y10t
  real y110
  real y111
  real y11t
  real y1s0
  real y1s1
  real y1st
  real yr00
  real yr01
  real yr0t
  real yr10
  real yr11
  real yr1t
  real yrs0
  real yrs1
  real z(npoint)
  real z000
  real z00t
  real z001
  real z010
  real z011
  real z01t
  real z0s0
  real z0s1
  real z0st
  real z100
  real z101
  real z10t
  real z110
  real z111
  real z11t
  real z1s0
  real z1s1
  real z1st
  real zr00
  real zr01
  real zr0t
  real zr10
  real zr11
  real zr1t
  real zrs0
  real zrs1
!
!  Here are the (R,S,T) coordinates of the tetrahedron that we can think
!  of as the "tile" we need to place in each sub-box.
!
!  The (R,S,T) coordinate system is assumed to range from 0 to 1.
!
  r_tab(1) =  0.2E+00
  s_tab(1) =  0.2E+00
  t_tab(1) =  0.2E+00

  r_tab(2) =  0.8E+00
  s_tab(2) =  0.2E+00
  t_tab(2) =  0.2E+00

  r_tab(3) =  0.5E+00
  s_tab(3) =  0.8E+00
  t_tab(3) =  0.2E+00

  r_tab(4) =  0.5E+00
  s_tab(4) =  0.5E+00
  t_tab(4) =  0.8E+00
!
!  Compute the (X,Y,Z) coordinates of the corners of the (U,V,W) box.
!  These really only need to be computed once ever.
!
  call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umin, vmin, wmin, &
    x000, y000, z000 )

  call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umax, vmin, wmin, &
    x100, y100, z100 )

  call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umin, vmax, wmin, &
    x010, y010, z010 )

  call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umax, vmax, wmin, &
    x110, y110, z110 )

  call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umin, vmin, wmax, &
    x001, y001, z001 )

  call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umax, vmin, wmax, &
    x101, y101, z101 )

  call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umin, vmax, wmax, &
    x011, y011, z011 )

  call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umax, vmax, wmax, &
    x111, y111, z111 )
!
!  Now figure out the (X,Y,Z) coordinates of the tile point with
!  given (R,S,T) coordinates.  This depends on the positions of 
!  all sorts of points on the corners, edges, and faces of the big box.
!
  do i = 1, npoint
!
!  Get the (R,S,T) coordinates of point I.
!
    r = r_tab(i)
    s = s_tab(i)
    t = t_tab(i)
!
!  Get the corresponding point (U,V,W) in the rectangular space.
!
    u = ( ( r1 - r ) * u0 + ( r - r0 ) * u1 ) / ( r1 - r0 )
    v = ( ( s1 - s ) * v0 + ( s - s0 ) * v1 ) / ( s1 - s0 )
    w = ( ( t1 - t ) * w0 + ( t - t0 ) * w1 ) / ( t1 - t0 )
!
!  Evaluate (X,Y,Z) on the 12 edges "near" (U,V,W).
!
    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umin, vmin, w, &
      x00t, y00t, z00t )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umin, vmax, w, &
      x01t, y01t, z01t )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umax, vmin, w, &
      x10t, y10t, z10t )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umax, vmax, w, &
      x11t, y11t, z11t )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umin, v, wmin, &
      x0s0, y0s0, z0s0 )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umin, v, wmax, &
      x0s1, y0s1, z0s1 )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umax, v, wmin, &
      x1s0, y1s0, z1s0 )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umax, v, wmax, &
      x1s1, y1s1, z1s1 )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, u, vmin, wmin, &
      xr00, yr00, zr00 )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, u, vmin, wmax, &
      xr01, yr01, zr01 )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, u, vmax, wmin, &
      xr10, yr10, zr10 )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, u, vmax, wmax, &
      xr11, yr11, zr11 )
!
!  Evaluate (X,Y,Z) on the six faces near (U,V,W).
!
    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umin, v, w, &
      x0st, y0st, z0st )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, umax, v, w, &
      x1st, y1st, z1st )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, u, vmin, w, &
      xr0t, yr0t, zr0t )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, u, vmax, w, &
      xr1t, yr1t, zr1t )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, u, v, wmin, &
      xrs0, yrs0, zrs0 )

    call boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, u, v, wmax, &
      xrs1, yrs1, zrs1 )
!
!  Now figure out the location of point 
!    I => (R,S,T) => (U,V,W) => (X(I),Y(I),Z(I)).
!
    x(i) = ( ( umax - u    ) * ( vmax - v    ) * ( wmax - w    ) * x000 &
           - ( umax - u    ) * ( vmax - v    ) * ( wmax - wmin ) * x00t &
           + ( umax - u    ) * ( vmax - v    ) * ( w    - wmin ) * x001 &
           - ( umax - u    ) * ( vmax - vmin ) * ( wmax - w    ) * x0s0 &
           + ( umax - u    ) * ( vmax - vmin ) * ( wmax - wmin ) * x0st &
           - ( umax - u    ) * ( vmax - vmin ) * ( w    - wmin ) * x0s1 &
           + ( umax - u    ) * ( v    - vmin ) * ( wmax - w    ) * x010 &
           - ( umax - u    ) * ( v    - vmin ) * ( wmax - wmin ) * x01t &
           + ( umax - u    ) * ( v    - vmin ) * ( w    - wmin ) * x011 &
           - ( umax - umin ) * ( vmax - v    ) * ( wmax - w    ) * xr00 &
           + ( umax - umin ) * ( vmax - v    ) * ( wmax - wmin ) * xr0t &
           - ( umax - umin ) * ( vmax - v    ) * ( w    - wmin ) * xr01 &
           + ( umax - umin ) * ( vmax - vmin ) * ( wmax - w    ) * xrs0 &
           + ( umax - umin ) * ( vmax - vmin ) * ( w    - wmin ) * xrs1 &
           - ( umax - umin ) * ( v    - vmin ) * ( wmax - w    ) * xr10 &
           + ( umax - umin ) * ( v    - vmin ) * ( wmax - wmin ) * xr1t &
           - ( umax - umin ) * ( v    - vmin ) * ( w    - wmin ) * xr11 &
           + ( u    - umin ) * ( vmax - v    ) * ( wmax - w    ) * x100 &
           - ( u    - umin ) * ( vmax - v    ) * ( wmax - wmin ) * x10t &
           + ( u    - umin ) * ( vmax - v    ) * ( w    - wmin ) * x101 &
           - ( u    - umin ) * ( vmax - vmin ) * ( wmax - w    ) * x1s0 &
           + ( u    - umin ) * ( vmax - vmin ) * ( wmax - wmin ) * x1st &
           - ( u    - umin ) * ( vmax - vmin ) * ( w    - wmin ) * x1s1 &
           + ( u    - umin ) * ( v    - vmin ) * ( wmax - w    ) * x110 &
           - ( u    - umin ) * ( v    - vmin ) * ( wmax - wmin ) * x11t &
           + ( u    - umin ) * ( v    - vmin ) * ( w    - wmin ) * x111 ) &
         / ( ( umax - umin ) * ( vmax - vmin ) * ( wmax - wmin ) )

    y(i) = ( ( umax - u    ) * ( vmax - v    ) * ( wmax - w    ) * y000 &
           - ( umax - u    ) * ( vmax - v    ) * ( wmax - wmin ) * y00t &
           + ( umax - u    ) * ( vmax - v    ) * ( w    - wmin ) * y001 &
           - ( umax - u    ) * ( vmax - vmin ) * ( wmax - w    ) * y0s0 &
           + ( umax - u    ) * ( vmax - vmin ) * ( wmax - wmin ) * y0st &
           - ( umax - u    ) * ( vmax - vmin ) * ( w    - wmin ) * y0s1 &
           + ( umax - u    ) * ( v    - vmin ) * ( wmax - w    ) * y010 &
           - ( umax - u    ) * ( v    - vmin ) * ( wmax - wmin ) * y01t &
           + ( umax - u    ) * ( v    - vmin ) * ( w    - wmin ) * y011 &
           - ( umax - umin ) * ( vmax - v    ) * ( wmax - w    ) * yr00 &
           + ( umax - umin ) * ( vmax - v    ) * ( wmax - wmin ) * yr0t &
           - ( umax - umin ) * ( vmax - v    ) * ( w    - wmin ) * yr01 &
           + ( umax - umin ) * ( vmax - vmin ) * ( wmax - w    ) * yrs0 &
           + ( umax - umin ) * ( vmax - vmin ) * ( w    - wmin ) * yrs1 &
           - ( umax - umin ) * ( v    - vmin ) * ( wmax - w    ) * yr10 &
           + ( umax - umin ) * ( v    - vmin ) * ( wmax - wmin ) * yr1t &
           - ( umax - umin ) * ( v    - vmin ) * ( w    - wmin ) * yr11 &
           + ( u    - umin ) * ( vmax - v    ) * ( wmax - w    ) * y100 &
           - ( u    - umin ) * ( vmax - v    ) * ( wmax - wmin ) * y10t &
           + ( u    - umin ) * ( vmax - v    ) * ( w    - wmin ) * y101 &
           - ( u    - umin ) * ( vmax - vmin ) * ( wmax - w    ) * y1s0 &
           + ( u    - umin ) * ( vmax - vmin ) * ( wmax - wmin ) * y1st &
           - ( u    - umin ) * ( vmax - vmin ) * ( w    - wmin ) * y1s1 &
           + ( u    - umin ) * ( v    - vmin ) * ( wmax - w    ) * y110 &
           - ( u    - umin ) * ( v    - vmin ) * ( wmax - wmin ) * y11t &
           + ( u    - umin ) * ( v    - vmin ) * ( w    - wmin ) * y111 ) &
         / ( ( umax - umin ) * ( vmax - vmin ) * ( wmax - wmin ) )

    z(i) = ( ( umax - u    ) * ( vmax - v    ) * ( wmax - w    ) * z000 &
           - ( umax - u    ) * ( vmax - v    ) * ( wmax - wmin ) * z00t &
           + ( umax - u    ) * ( vmax - v    ) * ( w    - wmin ) * z001 &
           - ( umax - u    ) * ( vmax - vmin ) * ( wmax - w    ) * z0s0 &
           + ( umax - u    ) * ( vmax - vmin ) * ( wmax - wmin ) * z0st &
           - ( umax - u    ) * ( vmax - vmin ) * ( w    - wmin ) * z0s1 &
           + ( umax - u    ) * ( v    - vmin ) * ( wmax - w    ) * z010 &
           - ( umax - u    ) * ( v    - vmin ) * ( wmax - wmin ) * z01t &
           + ( umax - u    ) * ( v    - vmin ) * ( w    - wmin ) * z011 &
           - ( umax - umin ) * ( vmax - v    ) * ( wmax - w    ) * zr00 &
           + ( umax - umin ) * ( vmax - v    ) * ( wmax - wmin ) * zr0t &
           - ( umax - umin ) * ( vmax - v    ) * ( w    - wmin ) * zr01 &
           + ( umax - umin ) * ( vmax - vmin ) * ( wmax - w    ) * zrs0 &
           + ( umax - umin ) * ( vmax - vmin ) * ( w    - wmin ) * zrs1 &
           - ( umax - umin ) * ( v    - vmin ) * ( wmax - w    ) * zr10 &
           + ( umax - umin ) * ( v    - vmin ) * ( wmax - wmin ) * zr1t &
           - ( umax - umin ) * ( v    - vmin ) * ( w    - wmin ) * zr11 &
           + ( u    - umin ) * ( vmax - v    ) * ( wmax - w    ) * z100 &
           - ( u    - umin ) * ( vmax - v    ) * ( wmax - wmin ) * z10t &
           + ( u    - umin ) * ( vmax - v    ) * ( w    - wmin ) * z101 &
           - ( u    - umin ) * ( vmax - vmin ) * ( wmax - w    ) * z1s0 &
           + ( u    - umin ) * ( vmax - vmin ) * ( wmax - wmin ) * z1st &
           - ( u    - umin ) * ( vmax - vmin ) * ( w    - wmin ) * z1s1 &
           + ( u    - umin ) * ( v    - vmin ) * ( wmax - w    ) * z110 &
           - ( u    - umin ) * ( v    - vmin ) * ( wmax - wmin ) * z11t &
           + ( u    - umin ) * ( v    - vmin ) * ( w    - wmin ) * z111 ) &
         / ( ( umax - umin ) * ( vmax - vmin ) * ( wmax - wmin ) )

  end do

  write ( 1, '(6f12.6)' ) x(1), y(1), z(1), 0.0, 0.0, 0.0E+00
  write ( 1, '(6f12.6)' ) x(2), y(2), z(2), 0.0, 0.0, 0.0E+00
  write ( 1, '(6f12.6)' ) x(3), y(3), z(3), 0.0, 0.0, 0.0E+00

  write ( 1, '(6f12.6)' ) x(2), y(2), z(2), 0.0, 0.0, 0.0E+00
  write ( 1, '(6f12.6)' ) x(1), y(1), z(1), 0.0, 0.0, 0.0E+00
  write ( 1, '(6f12.6)' ) x(4), y(4), z(4), 0.0, 0.0, 0.0E+00

  write ( 1, '(6f12.6)' ) x(3), y(3), z(3), 0.0, 0.0, 0.0E+00
  write ( 1, '(6f12.6)' ) x(4), y(4), z(4), 0.0, 0.0, 0.0E+00
  write ( 1, '(6f12.6)' ) x(1), y(1), z(1), 0.0, 0.0, 0.0E+00

  write ( 1, '(6f12.6)' ) x(4), y(4), z(4), 0.0, 0.0, 0.0E+00
  write ( 1, '(6f12.6)' ) x(3), y(3), z(3), 0.0, 0.0, 0.0E+00
  write ( 1, '(6f12.6)' ) x(2), y(2), z(2), 0.0, 0.0, 0.0E+00

  return
end
subroutine boundary_3d ( umin, vmin, wmin, umax, vmax, wmax, u, v, w, x, y, z )
!
!*******************************************************************************
!
!! BOUNDARY_3D returns the (X,Y,Z) coordinates of a point (U,V,W).
!
!
!  Discussion:
!
!    In this example, a single formula describes the mapping
!
!      (U,V,W) => (X,Y,Z).
!
!    It is more common (and more interesting) for the formula
!    to depend on which face of the boundary is being considered.
!
!    This routine is only called for points (U,V,W) where one of the 
!    values U, V and W is "extreme", so there are generally six cases
!    to consider, for general boundaries.  The coding has been set
!    up with this in mind.
!
!  Reference:
!
!    W N Gordon and Charles A Hall,
!    Construction of Curvilinear Coordinate Systems and Application to
!      Mesh Generation,
!    International Journal of Numerical Methods in Engineering,
!    Volume 7, pages 461-477, 1973.
!
!    Joe Thompson, Bharat Soni, Nigel Weatherill,
!    Handbook of Grid Generation,
!    CRC Press,
!    1999.
!
!  Modified:
!
!    01 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real UMIN, VMIN, WMIN, UMAX, VMAX, WMAX, the (U,V,W) coordinates
!    of the two opposite corners of the big box.
!
!    Input, real U, V, W, the (U,V,W) coordinates of a point in the big box.
!
!    Output, real X, Y, Z, the (X,Y,Z) coordinates of the point.
!
  implicit none
!
  real, parameter :: deg2rad = 3.14159265E+00 / 180.0E+00
  real psi
  real theta
  real u
  real umax
  real umin
  real v
  real vmax
  real vmin
  real w
  real wmax
  real wmin
  real x
  real y
  real z

  theta = u * deg2rad
  psi = w * deg2rad

       if ( u == umin ) then
    x = v * cos ( theta ) * cos ( psi )
    y = v * sin ( theta ) * cos ( psi )
    z = v                 * sin ( psi )
  else if ( u == umax ) then
    x = v * cos ( theta ) * cos ( psi )
    y = v * sin ( theta ) * cos ( psi )
    z = v                 * sin ( psi )
  else if ( v == vmin ) then
    x = v * cos ( theta ) * cos ( psi )
    y = v * sin ( theta ) * cos ( psi )
    z = v                 * sin ( psi )
  else if ( v == vmax ) then
    x = v * cos ( theta ) * cos ( psi )
    y = v * sin ( theta ) * cos ( psi )
    z = v                 * sin ( psi )
  else if ( w == wmin ) then
    x = v * cos ( theta ) * cos ( psi )
    y = v * sin ( theta ) * cos ( psi )
    z = v                 * sin ( psi )
  else if ( w == wmax ) then
    x = v * cos ( theta ) * cos ( psi )
    y = v * sin ( theta ) * cos ( psi )
    z = v                 * sin ( psi )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'BOUNDARY_3D - Fatal error!'
    stop
  end if

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

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

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

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

  return
end