

subroutine EOS_Omni_HelmholtzReadTable(eosfilename)

  use EOS_Omni_Module
  implicit none

  character(len=512) :: warnline
  character(len=256) :: eosfilename
  integer :: slength
  LOGICAL :: tablethere
  real*8 :: tstp, dstp
  integer :: i, j, iostatus, AllocateStatus

  ! onle allocate memory if we are actually using the Helmholtz EOS
  allocate( heos_dt(HEOSJMAX), heos_dtSqr(HEOSJMAX), heos_dtInv(HEOSJMAX), &
            heos_dtSqrInv(HEOSJMAX), heos_t(HEOSJMAX), STAT = AllocateStatus )
  if (AllocateStatus /= 0) then
    write(warnline,*) "Not enough memory to read EOS file" 
    write(*,*) warnline
    STOP
  endif     
            
  allocate( heos_dd(HEOSIMAX), heos_ddSqr(HEOSIMAX), heos_ddInv(HEOSIMAX), &
            heos_ddSqrInv(HEOSIMAX), heos_d(HEOSIMAX), STAT = AllocateStatus )
  if (AllocateStatus /= 0) then
    write(warnline,*) "Not enough memory to read EOS file" 
    write(*,*) warnline
    STOP
  endif
  
  allocate( heos_f(HEOSIMAX,HEOSJMAX), heos_ft(HEOSIMAX,HEOSJMAX), heos_ftt(HEOSIMAX,HEOSJMAX), &
    heos_fd(HEOSIMAX,HEOSJMAX), heos_fdd(HEOSIMAX,HEOSJMAX), heos_fdt(HEOSIMAX,HEOSJMAX), &
    heos_fddt(HEOSIMAX,HEOSJMAX), heos_fdtt(HEOSIMAX,HEOSJMAX), heos_fddtt(HEOSIMAX,HEOSJMAX), &
    heos_dpdf(HEOSIMAX,HEOSJMAX), heos_dpdft(HEOSIMAX,HEOSJMAX), heos_dpdfd(HEOSIMAX,HEOSJMAX), &
    heos_dpdfdt(HEOSIMAX,HEOSJMAX), heos_ef(HEOSIMAX,HEOSJMAX), heos_eft(HEOSIMAX,HEOSJMAX), &
    heos_efd(HEOSIMAX,HEOSJMAX), heos_efdt(HEOSIMAX,HEOSJMAX), heos_xf(HEOSIMAX,HEOSJMAX), &
    heos_xft(HEOSIMAX,HEOSJMAX), heos_xfd(HEOSIMAX,HEOSJMAX), heos_xfdt(HEOSIMAX,HEOSJMAX), &
    STAT = AllocateStatus )
  if (AllocateStatus /= 0) then
    write(warnline,*) "Not enough memory to read EOS file" 
    write(*,*) warnline
    STOP
  endif

  

  ! convert file name to fortran string
  !call CCTK_FortranString ( slength, helmeos_table_name, &
  !     eosfilename )

  ! read EOS table file
  inquire(file=trim(adjustl(eosfilename)),exist=tablethere)

  if(.not.tablethere) then
     write(warnline,"(A10,A,A15)") "EOS file ", trim(adjustl(eosfilename)), " does not exist!" 
     write(*,*) warnline
     STOP
  endif

  write(*,*) "Reading Helmholtz EOS table"
  write(*,*) 'filename: ', trim(adjustl(eosfilename))

  open(unit=19,file=trim(adjustl(eosfilename)),status='old')
      
  !..read the helmholtz free energy table
  do j=1,HEOSJMAX
    do i=1,HEOSIMAX
      read(19,*) heos_f(i,j),heos_fd(i,j),heos_ft(i,j), &
            heos_fdd(i,j),heos_ftt(i,j), & 
            heos_fdt(i,j), heos_fddt(i,j),heos_fdtt(i,j),heos_fddtt(i,j)
    enddo
  enddo

  !..read the pressure derivative with density table
  do j=1,HEOSJMAX
    do i=1,HEOSIMAX
      read(19,*) heos_dpdf(i,j),heos_dpdfd(i,j), &
            heos_dpdft(i,j),heos_dpdfdt(i,j)
    enddo
  enddo

  !..read the electron chemical potential table
  do j=1,HEOSJMAX
    do i=1,HEOSIMAX
      read(19,*) heos_ef(i,j),heos_efd(i,j),&
            heos_eft(i,j),heos_efdt(i,j)
    enddo
  enddo

  !..read the number density table
  do j=1,HEOSJMAX
    do i=1,HEOSIMAX
      read(19,*) heos_xf(i,j),heos_xfd(i,j),&
            heos_xft(i,j),heos_xfdt(i,j)
    enddo
  enddo

  close(unit=19)
       
  heos_tlo   = 4.0d0
  tstp  = (11.0d0 - heos_tlo)/float(HEOSJMAX-1)
  heos_tstpi = 1.0d0/tstp
  heos_dlo   = -10.0d0
  dstp  = (11.0d0 - heos_dlo)/float(HEOSIMAX-1)
  heos_dstpi = 1.0d0/dstp
  do j=1,HEOSJMAX
    heos_t(j) = 10.0d0**(heos_tlo + (j-1)*tstp)
    do i=1,HEOSIMAX
      heos_d(i) = 10.0d0**(heos_dlo + (i-1)*dstp)
    enddo
  enddo


  !..store the temperature and density differences and their inverses 
  do j=1,HEOSJMAX-1
    heos_dt(j)   = heos_t(j+1) - heos_t(j)
    heos_dtSqr(j)  = heos_dt(j)*heos_dt(j)
    heos_dtInv(j)  = 1.0d0/heos_dt(j)
    heos_dtSqrInv(j) = 1.0d0/heos_dtSqr(j)
  enddo
  do i=1,HEOSIMAX-1
    heos_dd(i)   = heos_d(i+1) - heos_d(i)
    heos_ddSqr(i)  = heos_dd(i)*heos_dd(i)
    heos_ddInv(i)  = 1.0d0/heos_dd(i)
    heos_ddSqrInv(i) = 1.0d0/heos_ddSqr(i)
  enddo 
       
  ! set some EOS limits
  eos_tempmin = heos_t(1) * temp_gf
  eos_epsmin = -0.007594153911d0 !HelmEOS: -0.007594153911 = - Qalpha/(4*m_u)
  eos_yemax = 0.99 ! set arbitrarily
  eos_yemin = 0.01 ! set arbitrarily
  
  write(*,*) 'imax=',HEOSIMAX,' jmax=',HEOSJMAX
  write(*,*) 'temp_min =',heos_t(1),' temp_max =',heos_t(HEOSJMAX)
  write(*,*) 'ye*den(1) =',heos_d(1),' ye*den(imax) =',heos_d(HEOSIMAX)
  write(*,*) "Done reading Helmholtz EOS table"


end subroutine EOS_Omni_HelmholtzReadTable