program write_str
!
! This program sets up a 3-D structured grid (cylindrical domain),
! initializes cell-centered and vertex-centered arrays on
! the grid, and then calls some fortran routines in the Rocketeer
! API to write out an HDF file (stsurf.hdf).
!
! Just for fun, surface meshes and data are also written
! to the HDF file. 
!
! Written by Robert Fiedler, 12/13/2002.
!
      USE mod_hdf_writer_API

      implicit NONE

      character*(10) filename     ! Name of the HDF file
      character*(32) long_name    ! Long name of variable
      integer ni,nj,nk            ! Dimensions of the arrays
      integer lghost              ! Number of layers of ghost cells
      integer mesh                ! mesh type
      LOGICAL :: append           ! Append or overwrite HDF file
      parameter(ni=21,nj=21,nk=21)
      REAL(MYREAL) x(ni,nj,nk),y(ni,nj,nk),z(ni,nj,nk) ! Coordinates
      REAL(MYREAL) dens(ni-1,nj-1,nk-1)  ! Cell centered data
      REAL(MYREAL) :: xdsp(ni,nj,nk)       ! Vertex centered vector comp
      REAL(MYREAL) :: ydsp(ni,nj,nk)       ! Vertex centered vector comp
      REAL(MYREAL) :: zdsp(ni,nj,nk)       ! Vertex centered vector comp
      REAL(MYREAL) :: xx(ni-1,nj-1,nk-1),xy(ni-1,nj-1,nk-1)
      REAL(MYREAL) :: xz(ni-1,nj-1,nk-1),yx(ni-1,nj-1,nk-1)
      REAL(MYREAL) :: yy(ni-1,nj-1,nk-1),yz(ni-1,nj-1,nk-1)
      REAL(MYREAL) :: zx(ni-1,nj-1,nk-1),zy(ni-1,nj-1,nk-1)
      REAL(MYREAL) :: zz(ni-1,nj-1,nk-1)   ! Cell centered tensor comps
      CHARACTER*(8) time                   ! Physical problem time
!
      integer i,j,k
      REAL(MYREAL) ratio,ax,ay,az,dens0,densm,r
      REAL(MYREAL) rj,rjm1,rjm2,phi,dphi,two_pi
      data ax,ay,az,dens0,densm /1.0,1.0,1.0,1.0,100.0/
!
      filename = 'stsurf.hdf'
      write(time,'(F4.2)') 1.0
      lghost = 0
      mesh = 2
      append = .false.
!
! Set up a non-uniform structured grid
!
      ratio = 1.05
      do k=1,nk
        do j=1,nj
          z(1,j,k) = 0.0
          z(2,j,k) = 0.01
          do i=3,ni
            z(i,j,k) = z(i-1,j,k) + (z(i-1,j,k) - z(i-2,j,k))*ratio
          end do
        end do
      end do
!
      ratio = 1.05
      two_pi = 4.0*asin(1.0)
      dphi = two_pi/(nk -1)
      do k=1,nk
        do i=1,ni
          do j=1,nj
            if (j.eq.1) then
              rj = 0.0
              rjm1 = 0.0
            else if (j.eq.2) then
              rj = 0.01
              rjm1 = 0.0
            else
              rjm2 = rjm1
              rjm1 = rj
              rj = rjm1 + (rjm1 - rjm2) * ratio
            endif
            phi = (k-1)*dphi
            x(i,j,k) = rj*cos(phi)
            y(i,j,k) = rj*sin(phi)
          end do
        end do
      end do
!
! Make up some scalar data
!
      do k=1,nk-1
        do j=1,nj-1
          do i=1,ni-1
            r = (x(i,j,k)/ax)**2 + (y(i,j,k)/ay)**2 + (z(i,j,k)/az)**2
            dens(i,j,k) = densm * exp(-r) + dens0
          end do
        end do
      end do
!
! Make up some vector data
!
      do k=1,nk
        do j=1,nj
          do i=1,ni
            r = (x(i,j,k)/ax)**3 + (y(i,j,k)/ay)**2 + (z(i,j,k)/az)**2
            xdsp(i,j,k) = 10.0*exp(-r)/densm
            ydsp(i,j,k) = -20.0*exp(-2.0*r)/densm
            zdsp(i,j,k) = 30.0*exp(-3.0*r)/densm
          end do
        end do
      end do
!
! Make up some tensor data
!
      do k=1,nk-1
        do j=1,nj-1
          do i=1,ni-1
            xx(i,j,k)=xdsp(i+1,j,k)-xdsp(i,j,k)
            xy(i,j,k)=xdsp(i,j+1,k)-xdsp(i,j,k)
            xz(i,j,k)=xdsp(i,j,k+1)-xdsp(i,j,k)
            yx(i,j,k)=ydsp(i+1,j,k)-ydsp(i,j,k)
            yy(i,j,k)=ydsp(i,j+1,k)-ydsp(i,j,k)
            yz(i,j,k)=ydsp(i,j,k+1)-ydsp(i,j,k)
            zx(i,j,k)=zdsp(i+1,j,k)-zdsp(i,j,k)
            zy(i,j,k)=zdsp(i,j+1,k)-zdsp(i,j,k)
            zz(i,j,k)=zdsp(i,j,k+1)-zdsp(i,j,k)
          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_str  &
       (filename, filename, 'block_001', 'gas', 'm', 'x', 'y', 'z',  &
        ni, nj, nk, 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',  &
        ni-1, nj-1, nk-1, lghost, dens)
!
! Write vector variable (node centered)
!
      long_name = 'Displacement at time '//time

      call w_vector_str  &
       (filename, long_name, 'disp', 'm',  &
        ni, nj, nk, lghost, xdsp, ydsp, zdsp)
!
! Write tensor variable (cell centered)
!
      long_name = 'Stress at time '//time

      call w_tensor_str  &
       (filename, long_name, 'stress', 'mmm',  &
        ni-1, nj-1, nk-1, lghost, xx, xy, xz, yx, yy, yz, zx, zy, zz)

!.......................................................................

! Add surface mesh for i = 1.

      append = .true.

! Write block header

      call w_block_header  &
       (filename, filename, 'block_001', time, 'sur1', mesh,  &
        lghost, 0, append)
!
! Write geometry data
!
      call w_geometry_str  &
       (filename, filename, 'block_001', 'sur1', 'm', 'x', 'y', 'z',  &
        1, nj, nk, lghost, x(1:1,:,:), y(1:1,:,:), z(1:1,:,:))
!
! Write scalar variable (cell centered)
!
      long_name = 'Density at time '//time

      call w_scalar_str  &
       (filename, long_name, 'ds1', 'kg/m^3',  &
        1, nj-1, nk-1, lghost, dens(1:1,:,:))
!
! Write vector variable (node centered)
!
      long_name = 'Displacement at time '//time

      call w_vector_str  &
       (filename, long_name, 'disps1', 'm',  &
        1, nj, nk, lghost, xdsp(1:1,:,:), ydsp(1:1,:,:), zdsp(1:1,:,:))


!.......................................................................

! Add surface mesh for j = nj.

! Write block header

      call w_block_header  &
       (filename, filename, 'block_001', time, 'sur2', mesh,  &
        lghost, 0, append)
!
! Write geometry data
!
      call w_geometry_str  &
       (filename, filename, 'block_001', 'sur2', 'm', 'x', 'y', 'z',  &
        ni, 1, nk, lghost, x(:,nj:nj,:), y(:,nj:nj,:), z(:,nj:nj,:))
!
! Write scalar variable (cell centered)
!
      long_name = 'Density at time '//time

      call w_scalar_str  &
       (filename, long_name, 'ds2', 'kg/m^3',  &
        ni-1, 1, nk-1, lghost, dens(:,nj-1:nj-1,:))
!
! Write vector variable (node centered)
!
      long_name = 'Displacement at time '//time

      call w_vector_str  &
       (filename, long_name, 'disps2', 'm',  &
        ni, 1, nk, lghost, xdsp(:,nj:nj,:), ydsp(:,nj:nj,:),  &
                           zdsp(:,nj:nj,:))

!.......................................................................

! Add surface mesh for k = 1.

! Write block header

      call w_block_header  &
       (filename, filename, 'block_001', time, 'sur3', mesh,  &
        lghost, 0, append)
!
! Write geometry data
!
      call w_geometry_str  &
       (filename, filename, 'block_001', 'sur3', 'm', 'x', 'y', 'z',  &
        ni, nj, 1, lghost, x(:,:,1:1), y(:,:,1:1), z(:,:,1:1))
!
! Write scalar variable (cell centered)
!
      long_name = 'Density at time '//time

      call w_scalar_str  &
       (filename, long_name, 'ds3', 'kg/m^3',  &
        ni-1, nj-1, 1, lghost, dens(:,:,1:1))
!
! Write vector variable (node centered)
!
      long_name = 'Displacement at time '//time

      call w_vector_str  &
       (filename, long_name, 'disps3', 'm',  &
        ni, nj, 1, lghost, xdsp(:,:,1:1), ydsp(:,:,1:1), zdsp(:,:,1:1))

      end