program write_tets
!
! This program sets up a 3-D unstructured grid,
! initializes cell-centered and vertex-centered arrays on
! the grid, and then calls the Rocketeer API to write out an HDF
! file (ustets.hdf) for use with ROCKETEER.
!
! The grid is tetrahedral.
!
! Written by Robert Fiedler, 12/12/2002.
!
      USE mod_hdf_writer_API

      implicit NONE
!
      character*(10) filename     ! Name of the HDF file
      character*(10) filegeom          ! HDF geometry data file
      character*(32) string
      integer ne,nn               ! Number of elements, nodes
      parameter(ne=4,nn=6)
      REAL(MYREAL) xyz(nn,3)              ! Cartesian coordinates
      integer conn(ne,4)          ! Connectivity
      REAL(MYREAL) dens(ne)               ! Cell centered scalar
      REAL(MYREAL) dsp(nn,3)              ! Vertex centered vector
      CHARACTER*(8) time                  ! Physical problem time
!
      integer i,j,k,i1,i2,i3,i4
      REAL(MYREAL) x,y,z,r
      REAL(MYREAL) ax,ay,az,dens0,densm
      data ax,ay,az,dens0,densm /1.0,2.0,3.0,1.0,100.0/
!
! Set up tetrahedral mesh.
!
      xyz(1,1) = 0.
      xyz(1,2) = 0.
      xyz(1,3) = 0.
      xyz(2,1) = 1.
      xyz(2,2) = 0.
      xyz(2,3) = 0.
      xyz(3,1) = 1.
      xyz(3,2) = 1.
      xyz(3,3) = 0.
      xyz(4,1) = 0.
      xyz(4,2) = 1.
      xyz(4,3) = 0.
      xyz(5,1) = 0.5
      xyz(5,2) = 0.5
      xyz(5,3) = 0.5
      xyz(6,1) = 0.5
      xyz(6,2) = 0.5
      xyz(6,3) =-0.5
!
      conn(1,1) = 5
      conn(1,2) = 2
      conn(1,3) = 1
      conn(1,4) = 3
      conn(2,1) = 5
      conn(2,2) = 4
      conn(2,3) = 3
      conn(2,4) = 1
      conn(3,1) = 6
      conn(3,2) = 2
      conn(3,3) = 3
      conn(3,4) = 1
      conn(4,1) = 6
      conn(4,2) = 4
      conn(4,3) = 1
      conn(4,4) = 3
!
      filename = 'ustets.hdf'
      filegeom = 'ustets.hdf'
      write(time,'(F4.2)') 1.0
!
! Make up some data
!
      do i=1,ne
        x = 0.25*(xyz(conn(i,1),1) + xyz(conn(i,2),1) +  &
                  xyz(conn(i,3),1) + xyz(conn(i,4),1) )
        y = 0.25*(xyz(conn(i,1),2) + xyz(conn(i,2),2) +  &
                  xyz(conn(i,3),2) + xyz(conn(i,4),2) )
        z = 0.25*(xyz(conn(i,1),3) + xyz(conn(i,2),3) +  &
                  xyz(conn(i,3),3) + xyz(conn(i,4),3) )
        r = (x/ax)**2 + (y/ay)**4 + (z/az)
        dens(i) = densm * exp(-r) + dens0
        print *,'i = ',i,'  x,y,z = ',x,y,z,'  dens = ',dens(i)
      end do
!
      do i=1,nn
        r = xyz(i,1)**2 + xyz(i,2)**4 + xyz(i,3)
        dsp(i,1) = ax*exp(-r/ax)
        dsp(i,2) = ay*exp(-r/ay)
        dsp(i,3) = az*exp(-r/az)
      end do
!
! Call routines in the Rocketeer API to create a valid HDF
! file for Rocketeer.

! Write block header

      call w_block_header  &
       (filename, filegeom, 'block_001', time, 'solid', 3, 0,  &
        0, .false.)
!
! Write geometry data
!
      call w_geometry_unstr  &
       (filename, filegeom, 'block_001', 'solid', 'm',  &
        'x', 'y', 'z', nn, 0, ne, 4, conn, xyz(:,1), xyz(:,2), xyz(:,3))
!
! Write scalar variable (cell centered)
!
      string = 'Density at time '//time

      call w_scalar_unstr  &
       (filename, string, 'd', 'kg/m^3',  &
        ne, 0, dens)
!
! Write vector variable (node centered)
!
      string = 'Displacement at time '//time

      call w_vector_unstr  &
       (filename, string, 'disp', 'm',  &
        nn, 0, dsp(:,1), dsp(:,2), dsp(:,3))

      end