program write_spheres
!
! This program reads in point data from standard input
! and calls a routine to write out an HDF
! file (spher1.hdf) for use with ROCKETEER.
!
! Compile on the CSAR Suns with:
!
! f90 hdf_spheres.f90 -L/usr/commercial/HDF-4.1r5_32bit/lib -ldf -lz -ljpeg
!
! Compile on turing with:
!
! f90 hdf_spheres.f90 -L/home/jiao/lib -ldf -lz -ljpeg
!
! Compile on modi4 with:
!
! f90 -64 hdf_spheres.f90 -L/afs/ncsa/packages/hdf/IRIX64_6.5 -ldf -lz -ljpeg
!
! Written by Robert Fiedler, 12/13/2002.
!
      USE mod_hdf_writer_API

      implicit NONE
!
      character*(10) filename          ! Name of the HDF file
      character*(13) infile            ! Input data file
      character*(32) long_name
      integer nn                       ! Number of spheres
      REAL(MYREAL), ALLOCATABLE :: xyz(:,:)   ! Cartesian coordinates
      REAL(MYREAL), ALLOCATABLE :: radius(:)  ! Sphere radii
      REAL(MYREAL), ALLOCATABLE :: v(:,:)     ! Shpere velocity vector
      CHARACTER*(8) time                      ! Physical problem time
      REAL(MYREAL) r1,r2,r3,r4         ! Dummy
!
      integer i,n,indx,itype,iblock,ios
      integer mesh, elems, nconn, conn(1,1)  ! For API
      logical append

      data indx /1/    ! A time index to name the HDF file 
      data iblock /1/  ! A block index for the HDF file
!
! Generate HDF file name
!
      write(filename,'("spher",i1,".hdf")') indx
      write(time,'(F4.2)') real(indx)  ! A fake time to label the data
!
! Read in the coords of the centers of the spheres and the radii
!
      infile = 'spheres_input'
      open(8,file=infile,iostat=ios)
      if (ios .ne. 0) then
        print *,'hdf_spheres requires data file ', infile
        stop
      endif

! Count the lines to find array dimensions

      do n=1,2000000
        read(8,*,end=10) i,itype,r1,r2,r3,r4
      end do
10    continue
      nn = n-1
      print *,'Number of spheres = ',nn

! Allocate arrays

      allocate(xyz(1:nn,1:3),radius(1:nn),v(1:nn,1:3))
  
! Save the data this time

      rewind 8
      do n=1,nn
        read(8,*) i,itype,xyz(n,1),xyz(n,2),xyz(n,3),radius(n)
!     1                   ,v(n,1),v(n,2),v(n,3)
!
! Do not read velocity (for now)
 
! Make the radii vary (they were constant in the Rocpart file)
!
        radius(n) = radius(n) + i/1000.

      end do
!
! Make up a velocity
!
      do i=1,nn
        v(i,1) = i
        v(i,2) = abs(i - 1000)
        v(i,3) = 1000000/i
      end do  ! i

! Call routines in the Rocketeer API to create a valid HDF
! file for Rocketeer.

      mesh = 3  ! Unstructured
      nconn = 1  ! Implies point data
      elems = 1  ! Will be ignored
      append = .false.

! Write block header

      call w_block_header  &
       (filename, filename, 'block_001', time, 'spheres', mesh,  &
        0, 0, append)
!
! Write geometry data
!
      call w_geometry_unstr  &
       (filename, filename, 'block_001', 'spheres', 'm',  &
        'x', 'y', 'z', nn, 0, elems, nconn, conn,  &
         xyz(:,1), xyz(:,2), xyz(:,3))
!
! Write scalar variable (node centered)
!
      long_name = 'Radius at time '//time

      call w_scalar_unstr  &
       (filename, long_name, 'radius', 'm',  &
        nn, 0, radius)
!
! Write vector variable (node centered)
!
      long_name = 'Velocity at time '//time

      call w_vector_unstr  &
       (filename, long_name, 'v', 'm',  &
        nn, 0, v(:,1), v(:,2), v(:,3))

      end