!.....EOS module for the Hempel and Schaffner-Bielich EOS v1.0 in the
!.....Shen 98 format, which does not include electrons, positrons and
!.....photons
!.....RMF Parameterization: FSUgold
!.....Nuclear mass table: Roca-Maza
      module fsg_roca_shen98format_eos_module

      implicit none

!.....t_entries: number of temperature entries
!.....y_entries: number of electron fraction entries
!.....nb_entries: number of nb entries
!.....entries: number of EOS entries per grid-point
!.....eos(i,j,k,l): EOS array
!.....i=1, ..., t_entries: temperature index
!.....j=1, ..., y_entries: electron fraction index
!.....k=1, ..., nb_entries: baryon number density index
!.....l=1, ..., entries: EOS entry index

!.....brief list of the entries:
!.....==========================
!.....1	logarithm of baryon mass density: log(rho) [g/cm^3]
!.....2	baryon number density: nb [fm^-3]
!.....3	logarithm of total proton fraction: log(yp)
!.....4	total proton fraction: yp (=ye)
!.....5	free energy per baryon wrt 938 MeV: f/nb-938 [MeV]
!.....6	internal energy per baryon wrt amu=931.49432 MeV: e/nb-amu [MeV]
!.....7	entropy per baryon: s/nb [kB]
!.....8	averaged mass number of heavy nuclei: <A>
!.....9	averaged charge number of heavy nuclei: <Z>
!.....10	effective mass: m^* [MeV]
!.....11	unbound neutron mass fraction: X_n
!.....12	unbound proton mass fraction: X_p
!.....13	light nuclei mass fraction: X_s
!.....14	heavy nuclei mass fraction: X_A
!.....15	pressure p: [MeV/fm^3]
!.....16	neutron chemical potential wrt m_n: mu_n-m_n [MeV]
!.....17	proton chemical potential wrt m_n: mu_p-m_p [MeV]
!.....18	averaged mass number of light nuclei: <a>
!.....19	averaged charge number of light nuclei: <z>


!.....t: temperature array [MeV]
!.....y: nominal electron fraction array
!.....nb: nominal baryon number density array [fm^-3]
      integer, parameter :: t_entries = 81
      integer, parameter :: y_entries = 60
      integer, parameter :: nb_entries = 326
      integer, parameter :: entries = 19
      double precision, 
     & dimension(t_entries,y_entries,nb_entries,entries) :: eos
      double precision, dimension(t_entries) :: t
      double precision, dimension(y_entries) :: y
      double precision, dimension(nb_entries) :: nb
      
      contains

!.....This subroutine reads in the ASCII EOS table and calculates the
!.....nominal electron fraction and nominal baryon number density arrays
        subroutine eos_readin

        implicit none
        
!.....define the constants of the t, y, and nb array
        double precision, parameter :: nbin = 1.0d-12
        double precision, parameter :: tin = 0.1d0
        double precision, parameter :: yin = 0.0d0
        double precision, parameter :: dnb = 0.04d0
        double precision, parameter :: dt = 0.04d0
        double precision, parameter :: dy = 0.01d0  

        character(len=15) :: dummy
        integer :: i, j, k, l

        open(10,file='fsg_roca_eos_shen98format_v1.0.tab',status='old')
!.....t loop
        do i=1,t_entries
          read(10,*)
          read(10,21) dummy,t(i)
          read(10,*)
!.....ye loop
          do j=1,y_entries
!.....nb loop
            do k=1,nb_entries
              read(10,22) (eos(i,j,k,l),l=1,entries)
            end do
            read(10,*)
          end do
        end do
        close(10)

!.....Uncomment the following part, if you do not want to use the stored
!.....values of t (round-off error) which are read-in above, but the
!.....correctly calculated ones:
!!.....calculate the temperature array
!        do i=1,t_entries
!          t(i)=tin*10.0d0**(dble(i-1)*dt)
!        end do

!.....calculate the electron fraction array
        do j=0,y_entries
          y(j)=yin+dble(j)*dy
        end do
!.....calculate the baryon number density array
        do k=1,nb_entries
          nb(k)=nbin*10.0d0**(dble(k-1)*dnb)
        end do
        
  21  format(a15,1x,1pE12.6)
  22  format(1pE13.6,18(1x,1pE13.6))

        return
        end subroutine eos_readin

      end module fsg_roca_shen98format_eos_module