!********************************************************************************
!* This program generates a mixed grid in a unit cubic domain, and generates
!* grid and boundary condition files in an unstructured format (.ugrid).
!*
!*        written by Dr. Katate Masatsuka (info[at]cfdbooks.com),
!*
!* the author of useful CFD books, "I do like CFD" (http://www.cfdbooks.com).
!*
!* The domain is divided into two parts, upper and lower.
!* A tetrahedral grid is generated in the upper part and 
!* a prismatic grid is generated in the lower part.
!* It resembles a typical high-Reynolds-number viscous grid.
!*
!* This is Version 3 (2011).
!* 
!* This F90 program is written and made available for an educational purpose.
!*
!* This file may be updated in future.
!*
!* Katate Masatsuka, March 2011. http://www.cfdbooks.com
!********************************************************************************

!********************************************************************************
!* Input : nx, ny, nz ! Number of nodes in each coordinate direction (nz>>nx=ny)
!*                 zm ! Z-coordinate separating tetra and prismatic layers. 
!*
!* Output: mixgrid.ugrid     ! Grid file (in .ugrid format)
!*         mixgrid.mapbc     ! Boundary condition file
!*         mixgrid.lines_fmt ! Line information in the prismatic layer
!*         mixgrid_tecplot.dat          !  Volume grid file for viewing by Tecplot
!*         mixgrid_boundary_tecplot.dat !Boundary grid file for viewing by Tecplot
!*
!* Note: Information on UGRID format can be found at
!*       http://www.simcenter.msstate.edu/docs/solidmesh/ugridformat.html
!*
!* Note: mixgrid.ugrid and mixgrid.mapbc are the files to be read by a flow solver
!*       for running a simulation.
!*
!* Note: BC numbers in mixgrid.mapbc must be chosen from available numbers in the
!*       flow solver.
!*
!********************************************************************************
 program mixgrid_cube

  implicit none
  integer , parameter :: dp = selected_real_kind(15) !Double precision

! Structured data
  real(dp), allocatable, dimension(:,:,:) :: x, y, z

! Mapping from (i,j,k) to the node number
  integer , allocatable, dimension(:,:,:) :: ijk2n

! Unstructured data: these data will be used to write output files.
  real(dp), allocatable, dimension(:,:) :: xyz       ! Coordinates
  integer , allocatable, dimension(:,:) :: prz, tet, tria, quad ! Element connectivities

! Number of nodes in each coordinate direction
  integer :: nx, ny, nz
! Number of nodes in z-direction: lower and upper.
  integer :: nm, nz_lower, nz_upper
! Number of prisms, tetrahedra, triangles, quadrilaterals
  integer :: nprz, ntet, ntria, nquad
! Number of nodes
  integer :: nnodes
! Local numbering of a hexahedron
  integer :: hex_loc(8)
! Grid spacing in each coordinate direction
  real(dp) :: dx, dy
! Z-coordinate of the intersection between tetrahedra and prismatic parts.
  real(dp) :: zm
! Spacing in z for the upper and lower regions (uniform in each region).
  real(dp) :: dz_lower, dz_upper
! Local variables
  integer  :: i, j, k, os
  real(dp) :: zs
! Line info
  integer  :: n_lines, n_total_points, max_points, min_points, i_count, inode

!*******************************************************************************
! Unit cube to be gridded.
!        ______________  
!       /             /|      Boundary group numbers are assigned as below:
!      /             / |       1:  Boundary face on x = 0.0
! z=1 /____________ /  |       2:  Boundary face on x = 1.0
!     |  |         |   |       3:  Boundary face on y = 0.0
!     |  |         |   |       4:  Boundary face on y = 1.0
!     |  |         |   |       5:  Boundary face on z = 0.0
!     |  |y=1______|___|       6:  Boundary face on z = 1.0
!     |  /         |  /
!     | /          | /
!     |/___________|/
! x=y=z=0        x=1
!
!
!  Here is the side-view of the grid.
!   Isotropic tetrahedra grid in the region: z > zm (node index k =  1, nm)
!              Prismatic grid in the region: z < zm (node index k = nm, nz)
!
!                                                        Index in z:
!                    z ^
!                      |  1 ____________________________ k=nz=(nz_upper-1)+nz_lower 
!                      |   !                            ! 
!                      |   !                            !
!                      |   !                            !
!                      |   |                            |
!                      |   | Isotropic Tetrahedral grid |
!                      |   |                            |
!                      |   |                            |
!                      |   |                            |
!                      |   |                            |
!                      | zm|____________________________|k = nm = nz_lower 
!                      |   |                            |
!                      |   |                            |
!                      |   |          Prismatic         |
!                      |   |                            |
!                      |   |____________________________|
!                      |  0                              1
!                      -------------------------------------> x
!
!
!*******************************************************************************
! Number of points in each direction.
! Set these values such that nz >> nx = ny.
! Also specify the vertical limit of the prismatic region by zm.

  nx = 17
  ny = 17
  nz = 33
  zm = 0.10_dp

! Allocate structured arrays

  allocate( x(nx,ny,nz), y(nx,ny,nz), z(nx,ny,nz) )

! Grid spacing (uniform)

  dx = 1.0_dp / real(nx-1)
  dy = 1.0_dp / real(ny-1)

! First, find the number of nodes in z-direction in the upper part that makes
! the upper part a uniform grid, and store it in "nm".
! NB: we assume that nz > nx = ny.

   nz_upper = 0

  do k = 1, nz
   zs = 1.0_dp - dx*real(k-1)
   if (zs > zm) then
    nz_upper = nz_upper + 1
   else
    exit
   endif
  end do

! Add the intersection point to the number of nodes in the upper part.

   nz_upper = nz_upper + 1

! Now set "nm" as the number of nodes in the lower part.

   nz_lower = (nz - nz_upper) + 1
   write(*,*) " nz_upper = ", nz_upper
   write(*,*) " nz_lower = ", nz_lower

! If nz_lower<=0, nz is not large enough. Stop

  if (nz_lower <= 0) then
   write(*,*) " Negative nz_lower = ", nz_lower
   write(*,*) " Not enough nz. Increase nz and try again."
   stop
  endif

! Spacings in the lower and upper parts.

  dz_lower =          zm / real(nz_lower-1)
  dz_upper = (1.0_dp -zm)/ real(nz_upper-1)

  nm = nz_lower

! Generate points.
! Uniform in each region: upper and lower domains.

  do i = 1, nx
   do j = 1, ny
    do k = 1, nz

     x(i,j,k) = dx*real(i-1)
     y(i,j,k) = dy*real(j-1)

     if (k <= nm) then
       z(i,j,k) = dz_lower*real(k-1)
     else
       z(i,j,k) = dz_upper*real(k-1 -nm+1) + zm
     endif

    end do
   end do
  end do

! For smooth transition from prismatic layer to the tetrahedral region,
! you may want to apply some stretching in z-direction in the prismatic
! layer. I'll leave it to you.

! Construct a node array

! (1)Construct a map first.

    allocate(ijk2n(nx,ny,nz)) ! (i,j,k) to node map

     nnodes = 0

    do i = 1, nx
     do j = 1, ny
      do k = 1, nz
              nnodes = nnodes + 1
        ijk2n(i,j,k) = nnodes
      end do
     end do
    end do

! (2)Construct a node array using ijk2n(:,:,:)
!     xyz(i,1) = x at node i.
!     xyz(i,2) = y at node i.
!     xyz(i,3) = z at node i.

    allocate(xyz(nnodes,3))

    do i = 1, nx
     do j = 1, ny
      do k = 1, nz
       xyz(ijk2n(i,j,k), 1) = x(i,j,k)
       xyz(ijk2n(i,j,k), 2) = y(i,j,k)
       xyz(ijk2n(i,j,k), 3) = z(i,j,k)
      end do
     end do
    end do

! Nodes are now ordered from 1 to nnodes(=total # of nodes).
! The nodal coordinates are stored in a one-dimensional fashion:
! x_i = xyz(i,1), y_i = xyz(i,2), z_i = xyz(i,3), i=1,2,3,...,nnodes.

! So, we don't need x, y, z, in the (i,j,k) data format.
! Goodbye, x, y, z!

   deallocate(x,y,z)

! Now, construct tet element and tria boundary connectivity data:
!  Prz: 6 nodes that define the prism.
! Quad: 4 nodes that define the quad face.
!  Tet: 4 nodes that define the tetrahedron.
! Tria: 3 nodes that define the tria face.

! Allocate the arrays with known dimensions:
!  - Hexahedron is divided into 6 tetrahedra
!  - Quad faces are divided into 2 triangles.

  nprz  = 2 * (nx-1)*(ny-1)*(nm-1)
  ntet  = 6 * (nx-1)*(ny-1)*(nz-1 -(nm-1))
  ntria = 2 * ( 2*(nx-1)*(ny-1) + 2*(ny-1)*(nz_upper-1) + 2*(nz_upper-1)*(nx-1) )
  nquad = 2 * (ny-1)*(nz_lower-1) + 2*(nz_lower-1)*(nx-1)

  allocate( prz(nprz,8), tet(ntet,4), tria(ntria,5), quad(nquad,5) )

! Reset ntet and ntria, nhex and nquad.
   nprz = 0
   ntet = 0
  ntria = 0
  nquad = 0

! Loop over (i,j,k) and construct connectivity data.

   i_loop : do i = 1, nx-1
    j_loop : do j = 1, ny-1
     k_loop : do k = 1, nz-1

! At each step, we look at a hexahedron as shown below.
! We work with local numbering because it is convenient.
!
!  Local node numbering, {1,2,3,4,5,6,7,8}, is defined as follows.
!
!
!                         (i,j+1,k+1) 8----------------------7 (i+1,j+1,k+1)
!                                    /.                     /|
!                                   / .                    / |
!                                  /  .                   /  |
!                                 /   .                  /   |
!                                /    .                 /    |
!                               /     .                /     |
!                              /      .               /      |
!                             /       .              /       |
!                            /        .             /        |
!                           /         .            /         |
!                          /          .           /          |
!               (i,j,k+1) 5----------------------6 (i+1,j,k+1)
!                         !           .          !           !
!                         !           .          !           !
!                         !           .          !           !
!                         |           4..........|...........3
!                         | (i,j+1,k).           |          /(i+1,j+1,k)
!                         |         .            |         /
!                         |        .             |        /
!                         |       .              |       /
!                         |      .               |      /
!                         |     .                |     /
!                         |    .                 |    /
!                         |   .                  |   /
!                         |  .                   |  /
!                         | .                    | /
!                         |.                     |/
!                         1----------------------2
!                   (i,j,k)                     (i+1,j,k)
!

! Store the node numbers into hex_loc(1:8), so that
! we can look at the hex as defined by local numbers {1,2,3,4,5,6,7,8}.
!
!    Lower 4 vertices
     hex_loc(1) = ijk2n(i  , j  , k  )
     hex_loc(2) = ijk2n(i+1, j  , k  )
     hex_loc(3) = ijk2n(i+1, j+1, k  )
     hex_loc(4) = ijk2n(i  , j+1, k  )

!    Upper 4 vertices
     hex_loc(5) = ijk2n(i  , j  , k+1)
     hex_loc(6) = ijk2n(i+1, j  , k+1)
     hex_loc(7) = ijk2n(i+1, j+1, k+1)
     hex_loc(8) = ijk2n(i  , j+1, k+1)

! From here, we can think of the hex defined by the 
! numbers {1,2,3,4,5,6,7,8}.

!--- Hexahedron defined by the nodes 1265-4378

!--------------------------------------------------------------------------
! Lower Part: Divide the hexahedra into 2 prisms as follows.
!--------------------------------------------------------------------------

     if (k < nm) then

!     Prism #1
      nprz = nprz + 1
      prz(nprz,1) = hex_loc(1)
      prz(nprz,2) = hex_loc(2)
      prz(nprz,3) = hex_loc(4)
      prz(nprz,4) = hex_loc(5)
      prz(nprz,5) = hex_loc(6)
      prz(nprz,6) = hex_loc(8)

!     Prism #2
      nprz = nprz + 1
      prz(nprz,1) = hex_loc(4)
      prz(nprz,2) = hex_loc(2)
      prz(nprz,3) = hex_loc(3)
      prz(nprz,4) = hex_loc(8)
      prz(nprz,5) = hex_loc(6)
      prz(nprz,6) = hex_loc(7)

!--- Boundary faces if the hexa is adjacent to the boundary:
!    Note: Nodes are ordered so that the boundary face points INWARD.
!
!    1. Left boundary face (x=0) : Quad face
      if     (i == 1) then

              nquad = nquad + 1
      quad(nquad,1) = hex_loc(1)
      quad(nquad,2) = hex_loc(4)
      quad(nquad,3) = hex_loc(8)
      quad(nquad,4) = hex_loc(5)
      quad(nquad,5) = 1  !<------ Boundary group number

!    2. Right boundary face (x=1) : Quad face
      elseif (i == nx-1) then

              nquad = nquad + 1
      quad(nquad,1) = hex_loc(3)
      quad(nquad,2) = hex_loc(2)
      quad(nquad,3) = hex_loc(6)
      quad(nquad,4) = hex_loc(7)
      quad(nquad,5) = 2  !<------ Boundary group number

      endif

!    3. Front boundary face (y=0) : Quad face
      if     (j == 1) then

              nquad = nquad + 1
      quad(nquad,1) = hex_loc(2)
      quad(nquad,2) = hex_loc(1)
      quad(nquad,3) = hex_loc(5)
      quad(nquad,4) = hex_loc(6)
      quad(nquad,5) = 3  !<------ Boundary group number

!    4. Back boundary face (y=1) : Quad face
      elseif (j == ny-1) then

              nquad = nquad + 1
      quad(nquad,1) = hex_loc(4)
      quad(nquad,2) = hex_loc(3)
      quad(nquad,3) = hex_loc(7)
      quad(nquad,4) = hex_loc(8)
      quad(nquad,5) = 4  !<------ Boundary group number

      endif

!    5. Bottom boundary face (z=0) : Tria face
      if     (k == 1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(1)
      tria(ntria,3) = hex_loc(2)
      tria(ntria,4) = 5  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(2)
      tria(ntria,3) = hex_loc(3)
      tria(ntria,4) = 5  !<------ Boundary group number

!    6. Top boundary face (z=1) : Tria face
      elseif (k == nz-1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(8)
      tria(ntria,2) = hex_loc(7)
      tria(ntria,3) = hex_loc(6)
      tria(ntria,4) = 6  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(8)
      tria(ntria,2) = hex_loc(6)
      tria(ntria,3) = hex_loc(5)
      tria(ntria,4) = 6  !<------ Boundary group number

      endif

     else

!--------------------------------------------------------------------------
! Upper Part: Divide the hexahedra into 6 tetrahedra as follows.
!--------------------------------------------------------------------------

!    Tet #1
      ntet = ntet + 1
      tet(ntet,1) = hex_loc(1)
      tet(ntet,2) = hex_loc(4)
      tet(ntet,3) = hex_loc(5)
      tet(ntet,4) = hex_loc(6)

!    Tet #2
      ntet = ntet + 1
      tet(ntet,1) = hex_loc(4)
      tet(ntet,2) = hex_loc(8)
      tet(ntet,3) = hex_loc(5)
      tet(ntet,4) = hex_loc(6)


!    Tet #3
      ntet = ntet + 1
      tet(ntet,1) = hex_loc(3)
      tet(ntet,2) = hex_loc(2)
      tet(ntet,3) = hex_loc(6)
      tet(ntet,4) = hex_loc(4)
 
!    Tet #4
      ntet = ntet + 1
      tet(ntet,1) = hex_loc(3)
      tet(ntet,2) = hex_loc(6)
      tet(ntet,3) = hex_loc(7)
      tet(ntet,4) = hex_loc(4)


!    Tet #5
      ntet = ntet + 1
      tet(ntet,1) = hex_loc(1)
      tet(ntet,2) = hex_loc(2)
      tet(ntet,3) = hex_loc(4)
      tet(ntet,4) = hex_loc(6)

!    Tet #6
      ntet = ntet + 1
      tet(ntet,1) = hex_loc(4)
      tet(ntet,2) = hex_loc(7)
      tet(ntet,3) = hex_loc(8)
      tet(ntet,4) = hex_loc(6)

!--- Boundary faces if the hexa is adjacent to the boundary:
!    Note: Nodes are ordered so that the boundary face points INWARD.
!
!    1. Left boundary face (x=0)
      if     (i == 1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(8)
      tria(ntria,3) = hex_loc(5)
      tria(ntria,4) = 1  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(5)
      tria(ntria,3) = hex_loc(1)
      tria(ntria,4) = 1  !<------ Boundary group number

!    2. Right boundary face (x=1)
      elseif (i == nx-1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(3)
      tria(ntria,2) = hex_loc(2)
      tria(ntria,3) = hex_loc(6)
      tria(ntria,4) = 2  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(3)
      tria(ntria,2) = hex_loc(6)
      tria(ntria,3) = hex_loc(7)
      tria(ntria,4) = 2  !<------ Boundary group number

      endif

!    3. Front boundary face (y=0)
      if     (j == 1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(1)
      tria(ntria,2) = hex_loc(5)
      tria(ntria,3) = hex_loc(6)
      tria(ntria,4) = 3  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(1)
      tria(ntria,2) = hex_loc(6)
      tria(ntria,3) = hex_loc(2)
      tria(ntria,4) = 3  !<------ Boundary group number

!    4. Back boundary face (y=1)
      elseif (j == ny-1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(3)
      tria(ntria,3) = hex_loc(7)
      tria(ntria,4) = 4  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(7)
      tria(ntria,3) = hex_loc(8)
      tria(ntria,4) = 4  !<------ Boundary group number

      endif

!    5. Bottom boundary face (z=0)
      if     (k == 1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(1)
      tria(ntria,3) = hex_loc(2)
      tria(ntria,4) = 5  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(4)
      tria(ntria,2) = hex_loc(2)
      tria(ntria,3) = hex_loc(3)
      tria(ntria,4) = 5  !<------ Boundary group number

!    6. Top boundary face (z=1)
      elseif (k == nz-1) then

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(8)
      tria(ntria,2) = hex_loc(7)
      tria(ntria,3) = hex_loc(6)
      tria(ntria,4) = 6  !<------ Boundary group number

              ntria = ntria + 1
      tria(ntria,1) = hex_loc(8)
      tria(ntria,2) = hex_loc(6)
      tria(ntria,3) = hex_loc(5)
      tria(ntria,4) = 6  !<------ Boundary group number

      endif

     endif

    end do k_loop
   end do j_loop
  end do i_loop

  write(*,*)
  write(*,*) "            Nodes = ", nnodes
  write(*,*) "           Prisms = ", nprz
  write(*,*) "       Tetrahedra = ", ntet
  write(*,*) "   Quadrilaterals = ", nquad
  write(*,*) "        Triangles = ", ntria

!*******************************************************************************
! Write a UGRID file
!*******************************************************************************
  open(unit=2, file="mixgrid.ugrid", status="unknown", iostat=os)

!                    #nodes, #tri_faces, #quad_faces, #tetra, #pyr, #prz, #hex
   write(2,'(7i10)') nnodes,      ntria,       nquad,   ntet,    0, nprz,    0
   write(2, *)

! Nodal coordinates
  do i = 1, nnodes
   write(2,'(3es27.15)') (xyz(i,j), j=1,3)
  end do
 
! Triangular boundary faces
   do i = 1, ntria
    write(2,'(3i10)') tria(i,1), tria(i,2), tria(i,3)
   end do

! Quad boundary faces
   do i = 1, nquad
    write(2,'(4i10)') quad(i,1), quad(i,2), quad(i,3), quad(i,4)
   end do

! Face tag: Boundary group number
   do i = 1, ntria
    write(2,'(i10)') tria(i,4)
   end do

   do i = 1, nquad
    write(2,'(i10)') quad(i,5)
   end do

! Tetrahedra
   do i = 1, ntet
   write(2,'(4i10)') tet(i,1), tet(i,2), tet(i,3), tet(i,4)
   end do

! Prisms
   do i = 1, nprz
   write(2,'(6i10)') prz(i,1), prz(i,2), prz(i,3), prz(i,4), prz(i,5), prz(i,6)
   end do

  close(2)

!*******************************************************************************
! Write a boundary condition map file: boundary marks
! Note: Set appropriate boundary condition numbers in this file.
!*******************************************************************************
 open(unit=3, file="mixgrid.mapbc", status="unknown", iostat=os)

  write(3,*) "Boundary Group", "   BC Index"
  write(3,*) 1, 5000 !<--- These numbers must be chosen from available boundary
  write(3,*) 2, 5000 !     condition numbers in a flow solver that read this
  write(3,*) 3, 5000 !     file. Here, I just set 5000; I don't know what it is.
  write(3,*) 4, 5000 !
  write(3,*) 5, 5000 !
  write(3,*) 6, 5000 !

 close(3)

!*******************************************************************************
! Write a line information file which may be used for performing
! line-relaxations in a CFD solver. Line information is created for the prismatic
! region only (up to k=nm): list of nodes for each line.
!*******************************************************************************
  open(unit=5, file="mixgrid.lines_fmt", status="unknown", iostat=os)

  max_points = nm
  min_points = nm
  n_lines = nx*ny
  n_total_points = n_lines * nm

  write(5,*) n_lines, n_total_points, "Total lines and points"
  write(5,*) min_points, max_points, "Min and max points in line"

  i_count = 0

  do j = 1, ny
   do i = 1, nx

     i_count = i_count + 1
     write(5,*) nm, " Points in line for line = ", i_count
     inode = ijk2n(i,j,1)
     write(5,'(i10,a10,3es30.20)') inode, " x/y/z= ", xyz(inode,1), xyz(inode,2), xyz(inode,3)

    do k = 2, nm-1
     inode = ijk2n(i,j,k)
     write(5,'(i10)') inode
    end do

     inode = ijk2n(i,j,nm)
     write(5,'(i10,a10,3es30.20)') inode, " x/y/z= ", xyz(inode,1), xyz(inode,2), xyz(inode,3)

   end do

  end do

  if (i_count /= n_lines) write(*,*) "Error: i_count /= n_lines"

 close(5)

!*******************************************************************************
! Write a Tecplot volume grid file for viewing. Just for viewing.
!*******************************************************************************
 open(unit=1, file="mixgrid_tecplot.dat", status="unknown", iostat=os)

  write(1,*) 'title = "tet grid"'
  write(1,*) 'variables = "x","y","z"'

! Tetra Zone
  write(1,*) 'zone  n=', nnodes,',e=', ntet,' , et=tetrahedron, f=fepoint'
  do i = 1, nnodes
    write(1,'(3ES20.10)') (xyz(i,j),j=1,3)
  end do

  do i = 1, ntet
   write(1,'(4i10)') tet(i,1), tet(i,2), tet(i,3), tet(i,4)
  end do

! Prism zone
  write(1,*) 'zone  n=', nnodes,',e=', nprz,' , et=brick, f=fepoint'
  do i = 1, nnodes
    write(1,'(3es20.10)') (xyz(i,j),j=1,3)
  end do

  do i = 1, nprz
   write(1,'(8i10)') prz(i,1), prz(i,2), prz(i,3), prz(i,3), &
                     prz(i,4), prz(i,5), prz(i,6), prz(i,6)
  end do

 close(1)

!*******************************************************************************
! Write a Tecplot boundary grid file for viewing. Just for viewing.
!******************************************************************************
 open(unit=4, file="mixgrid_boundary_tecplot.dat", status="unknown", iostat=os)

  write(4,*) 'title = "tet grid boundary"'
  write(4,*) 'variables = "x","y","z"'
  write(4,*) 'zone  n=', nnodes,',e=', ntria+nquad,' , et=quadrilateral, f=fepoint'

  do i = 1, nnodes
    write(4,'(3es27.15)') (xyz(i,j),j=1,3)
  end do

  do i = 1, ntria
   write(4,'(4i10)') tria(i,1), tria(i,2), tria(i,3), tria(i,1)
  end do

  do i = 1, nquad
   write(4,'(4i10)') quad(i,1), quad(i,2), quad(i,3), quad(i,4)
  end do

 close(4)


!*******************************************************************************

  write(*,*)
  write(*,*) " Grid files have been successfully generated:"
  write(*,*) " --- mixgrid.ugrid (boundary faces are ordered pointing inward.)"
  write(*,*) " --- mixgrid.mapbc"
  write(*,*)
  write(*,*) " Tecplot files have been successfully generated:"
  write(*,*) " --- mixgrid_tecplot.dat"
  write(*,*) " --- mixgrid_boundary_tecplot.dat"
  write(*,*)
  write(*,*) " Line info file have been successfully generated:"
  write(*,*) " --- mixgrid.lines_fmt"

 stop

 end program mixgrid_cube