program write_unsurf
!
! This program sets up an unstructured quad surface mesh and 
! an unstructured hex element volumetric (full) mesh, 
! initializes cell-centered and vertex-centered arrays on
! the meshes, and then calls routines in the Rocketeer API to write
! out an HDF file (unsurf.hdf).
!
! The surface mesh covers only part of the surface of the
! volumetric mesh.
!
! Written by Robert Fiedler, 12/13/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) long_name

      integer ne,nn                    ! Number of elements, nodes
      parameter(ne=9,nn=16)
      REAL(MYREAL) xyz(nn,3)                   ! Cartesian coordinates
      integer conn(ne,4)                       ! Connectivity
      REAL(MYREAL) dens(ne)                    ! Cell centered data
      REAL(MYREAL) dsp(nn,3)                   ! nodal displacement
      CHARACTER*(8) time                       ! Physical problem time
!
      integer nef,nnf                  ! Elements, node of full mesh
      parameter(nef=ne*3,nnf=nn*4)     ! Simply related to surface mesh
      REAL(MYREAL) xyzf(nnf,3),densf(nef),dspf(nnf,3)  ! Full mesh quantities
      integer connf(nef,8)             ! Full mesh connectivity
!
      integer i,l, mesh
      LOGICAL append
      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 flat windowpaine surface mesh
!
      data xyz /0.,1.,2.,3.,0.,1.,2.,3.,0.,1.,2.,3.,0.,1.,2.,3.,  &
                0.,0.,0.,0.,1.,1.,1.,1.,2.,2.,2.,2.,3.,3.,3.,3.,  &
                0.,.2,.2,0.,0.,.2,.2,0.,0.,.2,.2,0.,0.,.2,.2,0./
!
          data conn/ 1, 2, 3, 5, 6, 7, 9,10,11,  &
                     2, 3, 4, 6, 7, 8,10,11,12,  &
                     6, 7, 8,10,11,12,14,15,16,  &
                     5, 6, 7, 9,10,11,13,14,15 /
!
      filename = 'unsurf.hdf'   ! Name of HDF file
      write(time,'(F4.2)') 0.
!
! Make up some data on the surface mesh
!
      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)**2 + (z/az)**2
        dens(i) = densm * exp(-r) + dens0
      end do
!
! Surface displacement
!
        do i=1,nn
          dsp(i,1) = xyz(i,1)*(-0.3)
          dsp(i,2) = xyz(i,2)*( 0.1)
          dsp(i,3) = xyz(i,3)*( 0.2)
        end do
!
! Now set up the full mesh coordinates and displacement
! The full mesh is obtained by sweeping the surface mesh in z.
!
      do l=1,4
        do i=1,nn
          xyzf(i+(l-1)*nn,1) = xyz(i,1)
          xyzf(i+(l-1)*nn,2) = xyz(i,2)
          xyzf(i+(l-1)*nn,3) = xyz(i,3) + real(l-1)
          dspf(i+(l-1)*nn,1) = dsp(i,1)
          dspf(i+(l-1)*nn,2) = dsp(i,2)
          dspf(i+(l-1)*nn,3) = dsp(i,3)
        enddo
      enddo
!
! The connectivity of the full mesh can be derived simply
! from the surface mesh in this case.
!
      do l=1,3
        do i=1,ne
          connf(i+(l-1)*ne,1) = conn(i,1) + (l-1)*nn
          connf(i+(l-1)*ne,2) = conn(i,2) + (l-1)*nn
          connf(i+(l-1)*ne,3) = conn(i,3) + (l-1)*nn
          connf(i+(l-1)*ne,4) = conn(i,4) + (l-1)*nn

          connf(i+(l-1)*ne,5) = conn(i,1) + (l-0)*nn
          connf(i+(l-1)*ne,6) = conn(i,2) + (l-0)*nn
          connf(i+(l-1)*ne,7) = conn(i,3) + (l-0)*nn
          connf(i+(l-1)*ne,8) = conn(i,4) + (l-0)*nn

          densf(i+(l-1)*ne) = dens(i)
        enddo
      enddo
!
! Call routines in the Rocketeer API to create a valid HDF
! file for Rocketeer.

! Begin with the quad surface mesh

        append = .false.
        mesh = 5

! Write block header

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

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

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

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


! Add the volume mesh (hexahedra)

        append = .true.
        mesh = 3

! Write block header

        call w_block_header  &
         (filename, filename, 'block_001', time, 'full', mesh, 0,  &
          0, append)
!
! Write geometry data
!
        call w_geometry_unstr  &
         (filename, filename, 'block_001', 'full', 'm',  &
          'x', 'y', 'z', nnf, 0, nef, 8, connf, xyzf(:,1), xyzf(:,2),  &
                                                           xyzf(:,3))
!
! Write scalar variable (cell centered)
!
        long_name = 'Density at time '//time

        call w_scalar_unstr  &
         (filename, long_name, 'df', 'kg/m^3',  &
          nef, 0, densf)
!
! Write vector variable (node centered)
!
        long_name = 'Displacement at time '//time

        call w_vector_unstr  &
         (filename, long_name, 'dispf', 'm',  &
          nnf, 0, dspf(:,1), dspf(:,2), dspf(:,3))

      end