program write_rect
!
! This program sets up a 3-D rectilinear grid,
! initializes data on the grid, and uses the Rocketeer
! API to write an HDF file (rect.hdf).
!
! Written by Robert Fiedler, 12/13/02.
!
      USE mod_hdf_writer_API

      implicit NONE
!
      character*(8) filename     ! Name of the HDF file
      character*(32) long_name    ! Long name of variable
      integer nx,ny,nz           ! Dimensions of the arrays
      integer lghost             ! Number of layers of ghost cells
      integer mesh                ! mesh type
      LOGICAL :: append           ! Append or overwrite HDF file
      parameter(nx=21,ny=21,nz=21)
      REAL(MYREAL) x(nx),y(ny),z(nz)     ! Cartesian coordinate arrays
      REAL(MYREAL) dens(nx-1,ny-1,nz-1)  ! Cell centered data
      REAL(MYREAL) xdsp(nx,ny,nz)        ! Vertex centered data
      CHARACTER*(8) time                 ! Physical problem time
!
      integer i,j,k
      REAL(MYREAL) ratio,ax,ay,az,dens0,densm,r
      data ax,ay,az,dens0,densm /1.0,1.0,1.0,1.0,100.0/
!
      filename = 'rect.hdf'
      write(time,'(F4.2)') 1.0
      lghost = 0
      mesh = 1  ! Rectilinear
      append = .false.
!
! Set up a non-uniform rectilinear grid
!
      ratio = 1.20
      x(1) = 0.0
      x(2) = 0.01
      do i=3,nx
        x(i) = x(i-1) + (x(i-1) - x(i-2))*ratio
      end do
!
      ratio = 1.10
      y(1) = 0.0
      y(2) = 0.01
      do i=3,ny
        y(i) = y(i-1) + (y(i-1) - y(i-2))*ratio
      end do
!
      ratio = 1.05
      z(1) = 0.0
      z(2) = 0.01
      do i=3,nz
        z(i) = z(i-1) + (z(i-1) - z(i-2))*ratio
      end do
!
! Make up some data
!
      do k=1,nz-1
        do j=1,ny-1
          do i=1,nx-1
            r = (x(i)/ax)**2 + (y(j)/ay)**2 + (z(k)/az)**2
            dens(i,j,k) = densm * exp(-r) + dens0
          end do
        end do
      end do
!
      do k=1,nz
        do j=1,ny
          do i=1,nx
            r = (x(i)/ax)**3 + (y(j)/ay)**2 + (z(k)/az)**2
            xdsp(i,j,k) = exp(-r)
          end do
        end do
      end do
!
! Call routines in the Rocketeer API to create a valid HDF
! file for Rocketeer.  First, we'll write the volumetric data.

! Write block header

      call w_block_header  &
       (filename, filename, 'block_001', time, 'gas', mesh, lghost,  &
        0, append)
!
! Write geometry data
!
      call w_geometry_rect  &
       (filename, filename, 'block_001', 'gas', 'm', 'x', 'y', 'z',  &
        nx, ny, nz, lghost, x, y, z)
!
! Write scalar variable (cell centered)
!
      long_name = 'Density at time '//time

      call w_scalar_str  &
       (filename, long_name, 'd', 'kg/m^3',  &
        nx-1, ny-1, nz-1, lghost, dens)

!
! Write scalar variable (vertex centered)
!
      long_name = 'Stretch at time '//time

      call w_scalar_str  &
       (filename, long_name, 'xd', 'kg/m^3',  &
        nx, ny, nz, lghost, xdsp)

      end