program readhmxb

! read in a 3D binary dataset from an HMXB simulation
!---------------------------------------------------------------------------

! LOCALS
CHARACTER(LEN=50) :: varfile, geofile, prefix

! dimensions of 3D dataset (radius, theta, phi)
INTEGER :: imax, jmax, kmax

! nsurf(j,k) gives the index of the first real radial zone along that radial ray
INTEGER, DIMENSION(:,:), ALLOCATABLE :: nsurf

! coordinate arrays (radius, theta, phi)
REAL, DIMENSION(:), ALLOCATABLE :: zxc, zyc, zzc

! data array - add more as necessary
REAL, DIMENSION(:,:,:), ALLOCATABLE :: density

!------------------------------------------------------------------------------
write(6,fmt="('Input prefix of data files: ')")
read (5,*) prefix

! read raw binary geometry information

geofile = trim(prefix) // '.geo'
open(unit=25,file=geofile,form='unformatted')
  read(25) imax, jmax, kmax
    ! allocate coordinate arrays
    ALLOCATE( zxc(imax) )
    ALLOCATE( zyc(jmax) )
    ALLOCATE( zzc(kmax) )
    ALLOCATE( nsurf(jmax,kmax) )
  read(25) zxc
  read(25) zyc
  read(25) zzc
  read(25) nsurf
close(25)


! ***************   DENSITY   ********************

  ALLOCATE(density(imax,jmax,kmax) )
  varfile = trim(prefix) // '.den'
  open(unit=15,file=varfile,form='unformatted')
   read(15) density
  close(15)


end