! fortran90 interface for Opendx
MODULE dxdump

implicit none

   INTERFACE dxPolar
     module procedure dxPol_scalar_dp
     module procedure dxPol_scalar_dp_pos
     module procedure dxPol_scalar_sp
     module procedure dxPol_scalar_sp_pos
     module procedure dxPol_vector_dp
     module procedure dxPol_vector_dp_pos
     module procedure dxPol_vector_sp
     module procedure dxPol_vector_sp_pos
   end INTERFACE

   INTERFACE dxCartesian
     module procedure dxCart_scalar_dp
     module procedure dxCart_scalar_dp_pos
     module procedure dxCart_scalar_sp
     module procedure dxCart_scalar_sp_pos
     module procedure dxCart_vector_dp
     module procedure dxCart_vector_dp_pos
     module procedure dxCart_vector_sp
     module procedure dxCart_vector_sp_pos
     module procedure dxCart_2d_scalar_dp
     module procedure dxCart_2d_scalar_dp_pos
     module procedure dxCart_2d_scalar_sp
     module procedure dxCart_2d_scalar_sp_pos
   end INTERFACE

   character(len=3),private :: endian='unk'
CONTAINS
!*******************************************************************************
   SUBROUTINE dxout(array, x, y, z, filenm)

     real, dimension(:,:,:):: array
     real, dimension(:)::  x, y, z
     character (len=*) :: filenm 

     integer :: nx, ny, nz, nitems
     integer :: iu

!     print*,('entering dxout') 

     nx=size(x, dim=1)
     ny=size(y, dim=1)
     nz=size(z, dim=1)

     nitems=nx*ny*nz

     iu=findiu()
     open(iu, file=trim(filenm)//'.dx')

!    create header file for OpenDX
	
     call getendian()
!  First object
     write(iu,*)'object 1 class array type float rank 0 items ', nitems,& 
                &endian,' ieee'
     write(iu,*)'data file ',trim(filenm)
     write(iu,*)' attribute "dep" string "positions"'

!  Second object
     write(iu,*)'object 2 class gridpositions counts ', nz, ny, nx
     write(iu,*)'origin ', z(1), y(1), x(1)
     write(iu,*)'delta ', 0, 0, z(2)-z(1)
     write(iu,*)'delta ', 0, y(2)-y(1), 0
     write(iu,*)'delta ', x(2)-x(1), 0, 0 

!  Third object
     write(iu,*)'object 3 class gridconnections counts ', nz, ny, nx
     write(iu,*)'attribute "element type" string "cubes"'
     write(iu,*)'attribute "ref" string "positions"'

!  connect objects
     write(iu,*)'object "', trim(filenm), '" class field'
     write(iu,*)' component "data" 1'
     write(iu,*)'component "positions" 2'
     write(iu,*)'component "connections" 3'

     close(iu)

     open(iu, file=trim(filenm),form='unformatted')
     write(iu) array
     close(iu)

   end SUBROUTINE dxout

!*******************************************************************************
SUBROUTINE dxPol_scalar_dp(data, r, t, z, filenm)

 character(len=*), intent(in) :: filenm
 real*8, dimension(:,:,:), intent(in) :: data
 real*8, dimension(:), intent(in) :: r, t, z
 real, allocatable, dimension(:,:,:):: datasp
 real, allocatable, dimension(:) :: rsp, tsp, zsp

! print*,('entering dxPol_scalar_dp')

 allocate(datasp(size(data,dim=1), size(data,dim=2), size(data,dim=3)))
 allocate(rsp(size(r)), tsp(size(t)), zsp(size(z)))

 datasp=data
 rsp=r
 tsp=t
 zsp=z

 call dxPol_scalar_sp(datasp, rsp, tsp, zsp, filenm)

 deallocate(datasp, rsp, tsp, zsp)
end SUBROUTINE dxPol_scalar_dp

!*******************************************************************************
SUBROUTINE dxPol_scalar_dp_pos(data, r, t, z, filenm, posfile)

 character(len=*), intent(in) :: filenm, posfile
 real*8, dimension(:,:,:), intent(in) :: data
 real*8, dimension(:), intent(in) :: r, t, z
 real, allocatable, dimension(:,:,:):: datasp
 real, allocatable, dimension(:) :: rsp, tsp, zsp

! print*,('entering dxPol_scalar_dp_pos')

 allocate(datasp(size(data,dim=1), size(data,dim=2), size(data,dim=3)))
 allocate(rsp(size(r)), tsp(size(t)), zsp(size(z)))

 datasp=data
 rsp=r
 tsp=t
 zsp=z

 call dxPol_scalar_sp_pos(datasp, rsp, tsp, zsp, filenm, posfile)

 deallocate(datasp, rsp, tsp, zsp)
end SUBROUTINE dxPol_scalar_dp_pos

!*******************************************************************************
SUBROUTINE dxPol_scalar_sp(data, r, t, z, filenm)

 real,allocatable :: pos(:,:)
 character(len=*), intent(in) :: filenm
 real, dimension(:,:,:), intent(in) :: data
 integer :: count, nr, nt, nz, nrc, ntc, nzc, i, j, k, iu
 real, dimension(:), intent(in) :: r, t, z
 real :: x, y

! print*,('entering dxPol_scalar_sp')

  nr=size(r)
  nt=size(t)
  nz=size(z)


  nrc=size(data, dim=1)
  ntc=size(data, dim=2)
  nzc=size(data, dim=3)

  if(.not.(nr==nrc.and.nt==ntc.and.nz==nzc)) then
    print*,'arguments do not match'
    stop 1
  endif

     iu=findiu()

! print out positions to file.pos
  open(iu, file=trim(filenm)//'.pos', form='unformatted')
!print *, ('ALLOCATING POSITIONS')
    allocate(pos(3,nr*nt*nz))
    count=0
    do i=1,nr
    do j=1,nt
    do k=1,nz
     count=count+1
! calculate x and y values
     x=r(i)*cos(t(j))
     y=r(i)*sin(t(j))
     pos(1,count)=x
     pos(2,count)=y
     pos(3,count)=z(k)
   enddo
   enddo
   enddo
!print *, ('ALLOCATING POSITIONS DONE')
   write(iu)pos
  close(iu)

! print out data to file.dat
  open(iu, file=trim(filenm)//'.dat', form='unformatted')
   write(iu)(((data(i,j,k),k=1,nz),j=1,nt),i=1,nr)
  close(iu)
  
  call getendian()
!print first object to header and positions
  open(iu, file=trim(filenm)//'.dx', form='formatted')
  write(iu,*)'object 1 class array type float rank 1 shape 3 items ', nr*nt*nz,&
                &endian,' ieee'
  write(iu,*)'data file '//trim(filenm)//'.pos, 4'
!print second object
   write(iu,*)'object 2 class gridconnections counts ', nr, nt, nz

!print third object and data to header
   write(iu,*)'object 3 class array type float rank 0 items ', nr*nt*nz,&
                &endian,' ieee'
   write(iu,*)'data file '//trim(filenm)//'.dat, 4'
  write(iu,*)'attribute "dep" string "positions"'
!print class field connecting all components to header
  write(iu,*)'object "irreg positions regular connections" class field'
  write(iu,*)'component "positions" value 1'
  write(iu,*)'component "connections" value 2'
  write(iu,*)'component "data" value 3'
  write(iu,*)'end'
  close(iu)

  deallocate(pos)
end SUBROUTINE dxPol_scalar_sp

!*******************************************************************************
SUBROUTINE dxPol_scalar_sp_pos(data, r, t, z, filenm, posfile)

 character(len=*), intent(in) :: filenm, posfile
 real, dimension(:,:,:), intent(in) :: data
 integer :: count, nr, nt, nz, nrc, ntc, nzc, i, j, k, iu
 real, dimension(:), intent(in) :: r, t, z
 real :: x, y
 logical :: ex

! print*,('entering dxPol_scalar_sp')

  nr=size(r)
  nt=size(t)
  nz=size(z)


  nrc=size(data, dim=1)
  ntc=size(data, dim=2)
  nzc=size(data, dim=3)

  if(.not.(nr==nrc.and.nt==ntc.and.nz==nzc)) then
    print*,'arguments do not match'
    stop 1
  endif

  inquire(file=posfile, exist=ex)
    if(.not.ex)then
      print*,'positions file does not exist'
      stop 2
    endif

     iu=findiu()

! print out data to file.dat
  open(iu, file=trim(filenm)//'.dat', form='unformatted')
   write(iu)(((data(i,j,k),k=1,nz),j=1,nt),i=1,nr)
  close(iu)

!print first object to header and positions
  call getendian()
  open(iu, file=trim(filenm)//'.dx', form='formatted')
  write(iu,*)'object 1 class array type float rank 1 shape 3 items ', nr*nt*nz,&
                &endian,' ieee'
  write(iu,*)'data file ',posfile,', 4'
!print second object
   write(iu,*)'object 2 class gridconnections counts ', nr, nt, nz

!print third object and data to header
   write(iu,*)'object 3 class array type float rank 0 items ', nr*nt*nz,&
                &endian,' ieee'
   write(iu,*)'data file '//trim(filenm)//'.dat, 4'
  write(iu,*)'attribute "dep" string "positions"'
!print class field connecting all components to header
  write(iu,*)'object "irreg positions regular connections" class field'
  write(iu,*)'component "positions" value 1'
  write(iu,*)'component "connections" value 2'
  write(iu,*)'component "data" value 3'
  write(iu,*)'end'
  close(iu)

end SUBROUTINE dxPol_scalar_sp_pos

!*******************************************************************************
SUBROUTINE dxCart_scalar_dp(data, x, y, z, filenm)

 character(len=*), intent(in) :: filenm
 real*8, dimension(:,:,:), intent(in) :: data
 real*8, dimension(:), intent(in) :: x, y, z
 real, allocatable, dimension(:,:,:):: datasp
 real, allocatable, dimension(:) :: xsp, ysp, zsp

! print*,('entering dxCart_scalar_dp')

 allocate(datasp(size(data,dim=1), size(data,dim=2), size(data,dim=3)))
 allocate(xsp(size(x)), ysp(size(y)), zsp(size(z)))

 datasp=data
 xsp=x
 ysp=y
 zsp=z

 call dxCart_scalar_sp(datasp, xsp, ysp, zsp, filenm)

 deallocate(datasp, xsp, ysp, zsp)
end SUBROUTINE dxCart_scalar_dp

!*******************************************************************************
SUBROUTINE dxCart_scalar_dp_pos(data, x, y, z, filenm, posfile)

 character(len=*), intent(in) :: filenm, posfile
 real*8, dimension(:,:,:), intent(in) :: data
 real*8, dimension(:), intent(in) :: x, y, z
 real, allocatable, dimension(:,:,:):: datasp
 real, allocatable, dimension(:) :: xsp, ysp, zsp

! print*,('entering dxCart_scalar_dp_pos')

 allocate(datasp(size(data,dim=1), size(data,dim=2), size(data,dim=3)))
 allocate(xsp(size(x)), ysp(size(y)), zsp(size(z)))

 datasp=data
 xsp=x
 ysp=y
 zsp=z

 call dxCart_scalar_sp_pos(datasp, xsp, ysp, zsp, filenm, posfile)

 deallocate(datasp, xsp, ysp, zsp)
end SUBROUTINE dxCart_scalar_dp_pos

!*******************************************************************************
SUBROUTINE dxCart_2d_scalar_dp(data, x, y, filenm)

 character(len=*), intent(in) :: filenm
 real*8, dimension(:,:), intent(in) :: data
 real*8, dimension(:), intent(in) :: x, y
 real, allocatable, dimension(:,:):: datasp
 real, allocatable, dimension(:) :: xsp, ysp

! print*,('entering dxCart_2d_scalar_dp')

 allocate(datasp(size(data,dim=1), size(data,dim=2)))
 allocate(xsp(size(x)), ysp(size(y)))

 datasp=data
 xsp=x
 ysp=y

 call dxCart_2d_scalar_sp(datasp, xsp, ysp, filenm)

 deallocate(datasp, xsp, ysp)
end SUBROUTINE dxCart_2d_scalar_dp

!*******************************************************************************
SUBROUTINE dxCart_2d_scalar_dp_pos(data, x, y, filenm, posfile) 

 character(len=*), intent(in) :: filenm, posfile
 real*8, dimension(:,:), intent(in) :: data
 real*8, dimension(:), intent(in) :: x, y
 real, allocatable, dimension(:,:):: datasp
 real, allocatable, dimension(:) :: xsp, ysp

! print*,('entering dxCart_2d_scalar_dp_pos')

 allocate(datasp(size(data,dim=1), size(data,dim=2)))
 allocate(xsp(size(x)), ysp(size(y)))

 datasp=data
 xsp=x
 ysp=y

 call dxCart_2d_scalar_sp_pos(datasp, xsp, ysp, filenm, posfile)

 deallocate(datasp, xsp, ysp)
end SUBROUTINE dxCart_2d_scalar_dp_pos

!*******************************************************************************
SUBROUTINE dxCart_2d_scalar_sp(data, x, y, filenm)

 real,allocatable :: pos(:,:)
 character(len=*), intent(in) :: filenm
 real, dimension(:,:), intent(in) :: data
 integer :: count, nx, ny, nxc, nyc, i, j, iu
 real, dimension(:), intent(in) :: x, y

! print*,('entering dxCart_2d_scalar_sp')
  nx=size(x)
  ny=size(y)

  nxc=size(data, dim=1)
  nyc=size(data, dim=2)

  if(.not.(nx==nxc.and.ny==nyc)) then
    print*,'arguments do not match'
    stop 1
  endif
     iu=findiu()

! print out positions to file.pos
  open(iu, file=trim(filenm)//'.pos', form='unformatted')
!print *, ('ALLOCATING POSITIONS')
    allocate(pos(2,nx*ny))
    count=0
    do i=1,nx
    do j=1,ny
     count=count+1
     pos(1,count)=x(i)
     pos(2,count)=y(j)
   enddo
   enddo
!print *, ('ALLOCATING POSITIONS DONE')
   write(iu)pos
  close(iu)

! print out data to file.dat
  open(iu, file=trim(filenm)//'.dat', form='unformatted')
   write(iu)((data(i,j),j=1,ny),i=1,nx)
  close(iu)

  call getendian()
!print first object to header and positions
  open(iu, file=trim(filenm)//'.dx', form='formatted')
  write(iu,*)'object 1 class array type float rank 1 shape 2 items ', nx*ny,endian,' ieee'
  write(iu,*)'data file '//trim(filenm)//'.pos, 4'
!print second object
   write(iu,*)'object 2 class gridconnections counts ', nx, ny

!print third object and data to header
   write(iu,*)'object 3 class array type float rank 0 items ', nx*ny,&
                &endian,' ieee'
   write(iu,*)'data file '//trim(filenm)//'.dat, 4'
  write(iu,*)'attribute "dep" string "positions"'
!print class field connecting all components to header
  write(iu,*)'object "irreg positions regular connections" class field'
  write(iu,*)'component "positions" value 1'
  write(iu,*)'component "connections" value 2'
  write(iu,*)'component "data" value 3'
  write(iu,*)'end'
  close(iu)

  deallocate(pos)
end SUBROUTINE dxCart_2d_scalar_sp  

!*******************************************************************************
SUBROUTINE dxCart_2d_scalar_sp_pos(data, x, y, filenm, posfile)

 character(len=*), intent(in) :: filenm, posfile
 real, dimension(:,:), intent(in) :: data
 integer :: count, nx, ny, nxc, nyc, i, j, iu
 real, dimension(:), intent(in) :: x, y
 logical :: ex

! print*,('entering dxCart_2d_scalar_sp_pos')
  nx=size(x)
  ny=size(y)

  nxc=size(data, dim=1)
  nyc=size(data, dim=2)

  if(.not.(nx==nxc.and.ny==nyc)) then
    print*,'arguments do not match'
    stop 1
  endif

  inquire(file=posfile, exist=ex)
    if(.not.ex)then
      print*,'positions file does not exist'
      stop 2
    endif

     iu=findiu()

! print out data to file.dat
  open(iu, file=trim(filenm)//'.dat', form='unformatted')
   write(iu)((data(i,j),j=1,ny),i=1,nx)
  close(iu)

  call getendian()
!print first object to header and positions
  open(iu, file=trim(filenm)//'.dx', form='formatted')
  write(iu,*)'object 1 class array type float rank 1 shape 2 items ', nx*ny,endian,' ieee'
  write(iu,*)'data file ', posfile,', 4'
!print second object
   write(iu,*)'object 2 class gridconnections counts ', nx, ny

!print third object and data to header
   write(iu,*)'object 3 class array type float rank 0 items ', nx*ny,&
                &endian,' ieee'
   write(iu,*)'data file '//trim(filenm)//'.dat, 4'
  write(iu,*)'attribute "dep" string "positions"'
!print class field connecting all components to header
  write(iu,*)'object "irreg positions regular connections" class field'
  write(iu,*)'component "positions" value 1'
  write(iu,*)'component "connections" value 2'
  write(iu,*)'component "data" value 3'
  write(iu,*)'end'
  close(iu)

end SUBROUTINE dxCart_2d_scalar_sp_pos
  
!*******************************************************************************
SUBROUTINE dxCart_scalar_sp(data, x, y, z, filenm)

 real,allocatable :: pos(:,:)
 character(len=*), intent(in) :: filenm
 real, dimension(:,:,:), intent(in) :: data
 integer :: count, nx, ny, nz, nxc, nyc, nzc, i, j, k, iu
 real, dimension(:), intent(in) :: x, y, z

! print*,('entering dxCart_scalar_sp')
  nx=size(x)
  ny=size(y)
  nz=size(z)


  nxc=size(data, dim=1)
  nyc=size(data, dim=2)
  nzc=size(data, dim=3)

  if(.not.(nx==nxc.and.ny==nyc.and.nz==nzc)) then
    print*,'arguments do not match'
    stop 1
  endif

     iu=findiu()

! print out positions to file.pos
  open(iu, file=trim(filenm)//'.pos', form='unformatted')
!print *, ('ALLOCATING POSITIONS')
    allocate(pos(3,nx*ny*nz))
    count=0
    do i=1,nx
    do j=1,ny
    do k=1,nz
     count=count+1
     pos(1,count)=x(i)
     pos(2,count)=y(j)
     pos(3,count)=z(k)
   enddo
   enddo
   enddo
!print *, ('ALLOCATING POSITIONS DONE')
   write(iu)pos
  close(iu)

! print out data to file.dat
  open(iu, file=trim(filenm)//'.dat', form='unformatted')
   write(iu)(((data(i,j,k),k=1,nz),j=1,ny),i=1,nx)
  close(iu)

  call getendian()
!print first object to header and positions
  open(iu, file=trim(filenm)//'.dx', form='formatted')
  write(iu,*)'object 1 class array type float rank 1 shape 3 items ', nx*ny*nz,&
                &endian,' ieee'
  write(iu,*)'data file '//trim(filenm)//'.pos, 4'
!print second object
   write(iu,*)'object 2 class gridconnections counts ', nx, ny, nz

!print third object and data to header
   write(iu,*)'object 3 class array type float rank 0 items ', nx*ny*nz,&
                &endian,' ieee'
   write(iu,*)'data file '//trim(filenm)//'.dat, 4'
  write(iu,*)'attribute "dep" string "positions"'
!print class field connecting all components to header
  write(iu,*)'object "irreg positions regular connections" class field'
  write(iu,*)'component "positions" value 1'
  write(iu,*)'component "connections" value 2'
  write(iu,*)'component "data" value 3'
  write(iu,*)'end'
  close(iu)

  deallocate(pos)
end SUBROUTINE dxCart_scalar_sp  

!*******************************************************************************
SUBROUTINE dxCart_scalar_sp_pos(data, x, y, z, filenm, posfile)

 character(len=*), intent(in) :: filenm, posfile
 real, dimension(:,:,:), intent(in) :: data
 integer :: count, nx, ny, nz, nxc, nyc, nzc, i, j, k, iu
 real, dimension(:), intent(in) :: x, y, z
 logical :: ex

! print*,('entering dxCart_scalar_sp_pos')
  nx=size(x)
  ny=size(y)
  nz=size(z)


  nxc=size(data, dim=1)
  nyc=size(data, dim=2)
  nzc=size(data, dim=3)

  if(.not.(nx==nxc.and.ny==nyc.and.nz==nzc)) then
    print*,'arguments do not match'
    stop 1
  endif

  inquire(file=posfile, exist=ex)
  if(.not.ex)then
    print*,'positions file does not exist'
    stop 1
  endif

     iu=findiu()

! print out data to file.dat
  open(iu, file=trim(filenm)//'.dat', form='unformatted')
   write(iu)(((data(i,j,k),k=1,nz),j=1,ny),i=1,nx)
  close(iu)

  call getendian()
!print first object to header and positions
  open(iu, file=trim(filenm)//'.dx', form='formatted')
  write(iu,*)'object 1 class array type float rank 1 shape 3 items ', nx*ny*nz,&
                &endian,' ieee'
  write(iu,*)'data file ', posfile,', 4'
!print second object
   write(iu,*)'object 2 class gridconnections counts ', nx, ny, nz

!print third object and data to header
   write(iu,*)'object 3 class array type float rank 0 items ', nx*ny*nz,&
                &endian,' ieee'
   write(iu,*)'data file '//trim(filenm)//'.dat, 4'
  write(iu,*)'attribute "dep" string "positions"'
!print class field connecting all components to header
  write(iu,*)'object "irreg positions regular connections" class field'
  write(iu,*)'component "positions" value 1'
  write(iu,*)'component "connections" value 2'
  write(iu,*)'component "data" value 3'
  write(iu,*)'end'
  close(iu)

end SUBROUTINE dxCart_scalar_sp_pos
  
!*******************************************************************************
SUBROUTINE dxPol_vector_dp(vr, vt, vz, r, t, z, filenm)

 character(len=*), intent(in) :: filenm
 real*8, dimension(:), intent(in) :: r, t, z
 real*8, dimension(:,:,:), intent(in) :: vr, vt, vz

 real, allocatable, dimension(:,:,:) :: vrsp, vtsp, vzsp
 real, allocatable, dimension(:) :: rsp, tsp, zsp

! print*,('entering dxPol_vector_dp')

 allocate(vrsp(size(vr,dim=1), size(vr,dim=2), size(vr,dim=3)))
 allocate(vtsp(size(vt,dim=1), size(vt,dim=2), size(vt,dim=3)))
 allocate(vzsp(size(vz,dim=1), size(vz,dim=2), size(vz,dim=3)))
 allocate(rsp(size(r)), tsp(size(t)), zsp(size(z)))
 
 vrsp=vr
 vtsp=vt
 vzsp=vz
 rsp=r
 tsp=t
 zsp=z

 call dxPol_vector_sp(vrsp, vtsp, vzsp, rsp, tsp, zsp, filenm)

 deallocate(vrsp, vtsp, vzsp, rsp, tsp, zsp)

end SUBROUTINE dxPol_vector_dp

!*******************************************************************************
SUBROUTINE dxPol_vector_dp_pos(vr, vt, vz, r, t, z, filenm, posfile)

 character(len=*), intent(in) :: filenm, posfile
 real*8, dimension(:), intent(in) :: r, t, z
 real*8, dimension(:,:,:), intent(in) :: vr, vt, vz

 real, allocatable, dimension(:,:,:) :: vrsp, vtsp, vzsp
 real, allocatable, dimension(:) :: rsp, tsp, zsp

! print*,('entering dxPol_vector_dp_pos')

 allocate(vrsp(size(vr,dim=1), size(vr,dim=2), size(vr,dim=3)))
 allocate(vtsp(size(vt,dim=1), size(vt,dim=2), size(vt,dim=3)))
 allocate(vzsp(size(vz,dim=1), size(vz,dim=2), size(vz,dim=3)))
 allocate(rsp(size(r)), tsp(size(t)), zsp(size(z)))
 
 vrsp=vr
 vtsp=vt
 vzsp=vz
 rsp=r
 tsp=t
 zsp=z

 call dxPol_vector_sp_pos(vrsp, vtsp, vzsp, rsp, tsp, zsp, filenm, posfile)

 deallocate(vrsp, vtsp, vzsp, rsp, tsp, zsp)

end SUBROUTINE dxPol_vector_dp_pos

!*******************************************************************************
SUBROUTINE dxPol_vector_sp(vr, vt, vz, r, t, z, filenm)

 real,allocatable :: pos(:,:), dat(:,:)
 character(len=*), intent(in) :: filenm
 integer :: count, nr, nt, nz, nvr, nvt, nvz, i, j, k, iu
 real, dimension(:), intent(in) :: r, t, z
 real, dimension(:,:,:), intent(in) :: vr, vt, vz
 real :: x, y, vx, vy 

! print*,('enteringdxPol_vector_sp')

  nr=size(r)
  nt=size(t)
  nz=size(z)

  nvr=size(vr,dim=1)
  nvt=size(vt,dim=2)
  nvz=size(vz,dim=3)


     iu=findiu()

! print out positions to file.pos
  open(iu, file=trim(filenm)//'.pos', form='unformatted')
!    print *,'allocating pos'
    allocate(pos(3,nr*nt*nz))
!    print *,'allocated pos'
    count=0
    do i=1,nr
    do j=1,nt
    do k=1,nz
     count=count+1
! calculate x and y values
     x=r(i)*cos(t(j))
     y=r(i)*sin(t(j))
     pos(1,count)=x
     pos(2,count)=y
     pos(3,count)=z(k)
   enddo
   enddo
   enddo
   write(iu)pos
  close(iu)

! print out data to file.dat
  open(iu, file=trim(filenm)//'.dat', form='unformatted')
!    print *,'allocating dat',nvr*nvt*nvz
    allocate(dat(3,nvr*nvt*nvz))
!    print *,'allocated dat'
    count=0
    do i=1,nvr
    do j=1,nvt
    do k=1,nvz
     count=count+1
! calculate vx and vy values
     vx=vr(i,j,k)*cos(t(j))-vt(i,j,k)*sin(t(j))
     vy=vr(i,j,k)*sin(t(j))+vt(i,j,k)*cos(t(j))  !KATE sin->cos
     dat(1,count)=vx
     dat(2,count)=vy
     dat(3,count)=vz(i,j,k)
   enddo
   enddo
   enddo
   write(iu)dat
  close(iu)

  call getendian()
!print first object to header and positions
  open(iu, file=trim(filenm)//'.dx', form='formatted')
  write(iu,*)'object 1 class array type float rank 1 shape 3 items ', nr*nt*nz,&
                &endian,' ieee'
  write(iu,*)'data file '//trim(filenm)//'.pos, 4'

!print second object
   write(iu,*)'object 2 class gridconnections counts ', nr, nt, nz

!print third object and data to header
   write(iu,*)'object 3 class array type float rank 1 shape 3 items ',& 
               & nr*nt*nz, endian,' ieee'
   write(iu,*)'data file '//trim(filenm)//'.dat, 4'
  write(iu,*)'attribute "dep" string "positions"'

!print class field connecting all components to header
  write(iu,*)'object "irreg positions regular connections" class field'
  write(iu,*)'component "positions" value 1'
  write(iu,*)'component "connections" value 2'
  write(iu,*)'component "data" value 3'
  write(iu,*)'end'
  close(iu)

  deallocate(pos, dat)

end SUBROUTINE dxPol_vector_sp

!*******************************************************************************
SUBROUTINE dxPol_vector_sp_pos(vr, vt, vz, r, t, z, filenm, posfile)

 real,allocatable :: dat(:,:)
 character(len=*), intent(in) :: filenm, posfile
 integer :: count, nr, nt, nz, nvr, nvt, nvz, i, j, k, iu
 real, dimension(:), intent(in) :: r, t, z
 real, dimension(:,:,:), intent(in) :: vr, vt, vz
 real :: x, y, vx, vy 
 logical :: ex

! print*,('enteringdxPol_vector_sp_pos')

  nr=size(r)
  nt=size(t)
  nz=size(z)

  nvr=size(vr,dim=1)
  nvt=size(vt,dim=2)
  nvz=size(vz,dim=3)

  inquire(file=posfile, exist=ex)
   if(.not.ex)then
     print*,'positions file does not exist'
     stop 1 
   endif


     iu=findiu()

! print out data to file.dat
  open(iu, file=trim(filenm)//'.dat', form='unformatted')
!    print *,'allocating dat',nvr*nvt*nvz
    allocate(dat(3,nvr*nvt*nvz))
!    print *,'allocated dat'
    count=0
    do i=1,nvr
    do j=1,nvt
    do k=1,nvz
     count=count+1
! calculate vx and vy values
     vx=vr(i,j,k)*cos(t(j))-vt(i,j,k)*sin(t(j))
     vy=vr(i,j,k)*sin(t(j))+vt(i,j,k)*cos(t(j)) ! changed sin->cos KATE
     dat(1,count)=vx
     dat(2,count)=vy
     dat(3,count)=vz(i,j,k)
   enddo
   enddo
   enddo
   write(iu)dat
  close(iu)

  call getendian()
!print first object to header and positions
  open(iu, file=trim(filenm)//'.dx', form='formatted')
  write(iu,*)'object 1 class array type float rank 1 shape 3 items ', nr*nt*nz,&
                &endian,' ieee'
  write(iu,*)'data file ',posfile,', 4'

!print second object
   write(iu,*)'object 2 class gridconnections counts ', nr, nt, nz

!print third object and data to header
   write(iu,*)'object 3 class array type float rank 1 shape 3 items ',& 
               & nr*nt*nz, endian,' ieee'
   write(iu,*)'data file '//trim(filenm)//'.dat, 4'
  write(iu,*)'attribute "dep" string "positions"'

!print class field connecting all components to header
  write(iu,*)'object "irreg positions regular connections" class field'
  write(iu,*)'component "positions" value 1'
  write(iu,*)'component "connections" value 2'
  write(iu,*)'component "data" value 3'
  write(iu,*)'end'
  close(iu)

  deallocate(dat)

end SUBROUTINE dxPol_vector_sp_pos

!*******************************************************************************
SUBROUTINE dxCart_vector_dp(vx, vy, vz, x, y, z, filenm)

 character(len=*), intent(in) :: filenm
 real*8, dimension(:), intent(in) :: x, y, z
 real*8, dimension(:,:,:), intent(in) :: vx, vy, vz

 real, allocatable, dimension(:,:,:) :: vxsp, vysp, vzsp
 real, allocatable, dimension(:) :: xsp, ysp, zsp

! print*,('entering dxCart_vector_dp')

 allocate(vxsp(size(vx,dim=1), size(vx,dim=2), size(vx,dim=3)))
 allocate(vysp(size(vy,dim=1), size(vy,dim=2), size(vy,dim=3)))
 allocate(vzsp(size(vz,dim=1), size(vz,dim=2), size(vz,dim=3)))
 allocate(xsp(size(x)), ysp(size(y)), zsp(size(z)))
 
 vxsp=vx
 vysp=vy
 vzsp=vz
 xsp=x
 ysp=y
 zsp=z

 call dxCart_vector_sp(vxsp, vysp, vzsp, xsp, ysp, zsp, filenm)

 deallocate(vxsp, vysp, vzsp, xsp, ysp, zsp)

end SUBROUTINE dxCart_vector_dp

!*******************************************************************************
SUBROUTINE dxCart_vector_dp_pos(vx, vy, vz, x, y, z, filenm, posfile)

 character(len=*), intent(in) :: filenm, posfile
 real*8, dimension(:), intent(in) :: x, y, z
 real*8, dimension(:,:,:), intent(in) :: vx, vy, vz

 real, allocatable, dimension(:,:,:) :: vxsp, vysp, vzsp
 real, allocatable, dimension(:) :: xsp, ysp, zsp

! print*,('entering dxCart_vector_dp_pos')

 allocate(vxsp(size(vx,dim=1), size(vx,dim=2), size(vx,dim=3)))
 allocate(vysp(size(vy,dim=1), size(vy,dim=2), size(vy,dim=3)))
 allocate(vzsp(size(vz,dim=1), size(vz,dim=2), size(vz,dim=3)))
 allocate(xsp(size(x)), ysp(size(y)), zsp(size(z)))
 
 vxsp=vx
 vysp=vy
 vzsp=vz
 xsp=x
 ysp=y
 zsp=z

 call dxCart_vector_sp_pos(vxsp, vysp, vzsp, xsp, ysp, zsp, filenm, posfile)

 deallocate(vxsp, vysp, vzsp, xsp, ysp, zsp)

end SUBROUTINE dxCart_vector_dp_pos

!*******************************************************************************
SUBROUTINE dxCart_vector_sp(vx, vy, vz, x, y, z, filenm)

 real,allocatable :: pos(:,:), dat(:,:)
 character(len=*), intent(in) :: filenm
 integer :: count, nx, ny, nz, nvx, nvy, nvz, i, j, k, iu
 real, dimension(:), intent(in) :: x, y, z
 real, dimension(:,:,:), intent(in) :: vx, vy, vz

! print *,'entering dxCart_vector_sp'

  nx=size(x)
  ny=size(y)
  nz=size(z)

  nvx=size(vx,dim=1)
  nvy=size(vy,dim=2)
  nvz=size(vz,dim=3)

     iu=findiu()

! print out positions to file.pos
  open(iu, file=trim(filenm)//'.pos', form='unformatted')
!    print *,'allocating pos'
    allocate(pos(3,nx*ny*nz))
!    print *,'allocated pos'
    count=0
    do i=1,nx
    do j=1,ny
    do k=1,nz
     count=count+1
     pos(1,count)=x(i)
     pos(2,count)=y(j)
     pos(3,count)=z(k)
   enddo
   enddo
   enddo
   write(iu)pos
  close(iu)

! print out data to file.dat
  open(iu, file=trim(filenm)//'.dat', form='unformatted')
!    print *,'allocating dat',nvx*nvy*nvz
    allocate(dat(3,nvx*nvy*nvz))
!    print *,'allocated dat'
    count=0
    do i=1,nvx
    do j=1,nvy
    do k=1,nvz
     count=count+1
     dat(1,count)=vx(i,j,k)
     dat(2,count)=vy(i,j,k)
     dat(3,count)=vz(i,j,k)
   enddo
   enddo
   enddo
   write(iu)dat
  close(iu)

  call getendian()
!print first object to header and positions
  open(iu, file=trim(filenm)//'.dx', form='formatted')
  write(iu,*)'object 1 class array type float rank 1 shape 3 items ', nx*ny*nz,&
                &endian,' ieee'
  write(iu,*)'data file '//trim(filenm)//'.pos, 4'

!print second object
   write(iu,*)'object 2 class gridconnections counts ', nx, ny, nz

!print third object and data to header
   write(iu,*)'object 3 class array type float rank 1 shape 3 items ',& 
               & nx*ny*nz, endian,' ieee'
   write(iu,*)'data file '//trim(filenm)//'.dat, 4'
  write(iu,*)'attribute "dep" string "positions"'

!print class field connecting all components to header
  write(iu,*)'object "irreg positions regular connections" class field'
  write(iu,*)'component "positions" value 1'
  write(iu,*)'component "connections" value 2'
  write(iu,*)'component "data" value 3'
  write(iu,*)'end'
  close(iu)

  deallocate(pos, dat)

end SUBROUTINE dxCart_vector_sp

!*******************************************************************************
SUBROUTINE dxCart_vector_sp_pos(vx, vy, vz, x, y, z, filenm, posfile)

 real,allocatable :: dat(:,:)
 character(len=*), intent(in) :: filenm, posfile
 integer :: count, nx, ny, nz, nvx, nvy, nvz, i, j, k, iu
 real, dimension(:), intent(in) :: x, y, z
 real, dimension(:,:,:), intent(in) :: vx, vy, vz
 logical :: ex

! print *,'entering dxCart_vector_sp_pos'

  nx=size(x)
  ny=size(y)
  nz=size(z)

  nvx=size(vx,dim=1)
  nvy=size(vy,dim=2)
  nvz=size(vz,dim=3)

  inquire(file=posfile, exist=ex)
    if(.not.ex)then
      print*,'positions file does not exist'
      stop 1
    endif
     iu=findiu()

! print out data to file.dat
  open(iu, file=trim(filenm)//'.dat', form='unformatted')
!    print *,'allocating dat',nvx*nvy*nvz
    allocate(dat(3,nvx*nvy*nvz))
!    print *,'allocated dat'
    count=0
    do i=1,nvx
    do j=1,nvy
    do k=1,nvz
     count=count+1
     dat(1,count)=vx(i,j,k)
     dat(2,count)=vy(i,j,k)
     dat(3,count)=vz(i,j,k)
   enddo
   enddo
   enddo
   write(iu)dat
  close(iu)

  call getendian()
!print first object to header and positions
  open(iu, file=trim(filenm)//'.dx', form='formatted')
  write(iu,*)'object 1 class array type float rank 1 shape 3 items ', nx*ny*nz,&
                &endian,' ieee'
  write(iu,*)'data file ', posfile,', 4'

!print second object
   write(iu,*)'object 2 class gridconnections counts ', nx, ny, nz

!print third object and data to header
   write(iu,*)'object 3 class array type float rank 1 shape 3 items ',& 
               & nx*ny*nz, endian,' ieee'
   write(iu,*)'data file '//trim(filenm)//'.dat, 4'
  write(iu,*)'attribute "dep" string "positions"'

!print class field connecting all components to header
  write(iu,*)'object "irreg positions regular connections" class field'
  write(iu,*)'component "positions" value 1'
  write(iu,*)'component "connections" value 2'
  write(iu,*)'component "data" value 3'
  write(iu,*)'end'
  close(iu)

  deallocate(dat)

end SUBROUTINE dxCart_vector_sp_pos

!*******************************************************************************
INTEGER FUNCTION findiu()
     logical :: od
	do findiu=1,1000
           inquire(findiu, opened=od)
           if(.not.od) exit
	enddo
        return
   end FUNCTION findiu

 subroutine getendian()
 integer*4 il
 integer*2 is
 integer :: iu
  if(endian.eq.'unk')then
    iu=findiu()
    open(iu,status='scratch',form='unformatted')
    il=7  ! anything with asymetric bit pattern
!  write out 4 byte integer and read back in the first 2 bytes
    write(iu)il
    rewind(iu)
    read(iu)is
    close(iu)
!  test for equality
    if(is.eq.il)then
      endian='lsb'
    else
      endian='msb'
    endif
  endif
  return
 end subroutine getendian

end MODULE dxdump