!#######################################################################
!   SUBROUTINE   ---- sub_eps_5phase ----                            ###
!                                                                    ###
!   Subroutine for calculating the pressure,                         ###
!   sound velocity, and temperature from the Saumon EOS              ###
!                                                                    ###
!                                                         2015/05/28 ###
!#######################################################################
!=======================================================================
!   Input Data  : d    : Density [kg/m^3]
!                 e    : Energy [J/kg]
!   Output Data : p    : Pressure [Pa]
!                 sv   : Sound Velocity [m/s]
!                 t    : Temperature [K]
!=======================================================================

subroutine sub_eos_5phase(d,e,p,sv,flag_out,t,s)
  use mod_parameter_5phase_eos
  implicit none

!---- Arguments in This Subroutine ----------------------------------
  real(8), intent(INOUT)  :: d
  real(8), intent(IN)  :: e
  real(8), intent(OUT) :: p
  real(8), intent(OUT) :: sv
  integer(4), intent(IN) :: flag_out
  real(8), intent(OUT) :: t
  real(8), intent(OUT) :: s

!---- Variables Used in This Subroutine -----------------------------
  integer(4) :: i,j
  integer(4) :: ja,jb
  integer(4) :: ia,ib,ic,id
  real(8) :: da,db
  real(8) :: ea,eb,ec,ed
  real(8) :: pa,pb,pc,pd,pac,pbd
  real(8) :: sva,svb,svc,svd,svac,svbd
  real(8) :: ta,tb,tc,td,tac,tbd
  real(8) :: sa,sb,sc,sd,sac,sbd
!--------------------------------------------------------------------

  !--- Check for Density is in Table --------------------------------
  if (d < d_5ph(1) .or. d > d_5ph(nd_5ph)) then
     write(60,'(a39,e16.8)') 'Density is out of table for 5-Phase EOS', d
     write(*,'(a39,e16.8)')  'Density is out of table for 5-Phase EOS', d
     if (d < d_5ph(1)) d = d_5ph(1)
     if (d > d_5ph(nd_5ph)) d = d_5ph(nd_5ph)
  endif
  !------------------------------------------------------------------

  do j=1,nd_5ph-1
     if (d >= d_5ph(j) .and. d <= d_5ph(j+1)) then
        ja = j
        jb = j+1
        exit
     endif
  enddo

  !--- Check for Energy is in Table ---------------------------------
  if (e < e_5ph(0,ja) .or. e < e_5ph(0,jb) .or.              &
       &   e > e_5ph(nt_5ph,ja) .or. e > e_5ph(nt_5ph,jb)) then
     write(60,'(a38,e16.8)') 'Energy is out of table for 5-Phase EOS', e
     write(*,'(a38,e16.8)')  'Energy is out of table for 5-Phase EOS', e
     write(*,*) ja,jb,d_5ph(ja),d_5ph(jb)
     write(*,*) e_5ph(1,ja),e_5ph(1,jb)
     write(*,*) e_5ph(nt_5ph,ja),e_5ph(nt_5ph,jb)
     stop
  endif
  !------------------------------------------------------------------

  do i=0,nt_5ph-1
     if (e >= e_5ph(i,ja) .and. e < e_5ph(i+1,ja)) then
        ia = i
        ic = i+1
        exit
     endif
  enddo

  do i=0,nt_5ph-1
     if (e >= e_5ph(i,jb) .and. e < e_5ph(i+1,jb)) then
        ib = i
        id = i+1
        exit
     endif
  enddo

  da = d_5ph(ja)
  db = d_5ph(jb)

  ea = e_5ph(ia,ja)
  eb = e_5ph(ib,jb)
  ec = e_5ph(ic,ja)
  ed = e_5ph(id,jb)

  !--- Calculate Pressure ------------------------
  pa = p_5ph(ia,ja)
  pb = p_5ph(ib,jb)
  pc = p_5ph(ic,ja)
  pd = p_5ph(id,jb)

  pac = pa+(pc-pa)/(ec-ea)*(e-ea)
  pbd = pb+(pd-pb)/(ed-eb)*(e-eb)

  p = pac+(pbd-pac)/(db-da)*(d-da)

  !====
  if (p <= 0.0e0_8) p=0.0e0_8
  !====
  
!  write(88,'(3i5,3e16.8)') ja,ia,ib,d,da,db

  !--- Calculate Sound Velocity ------------------
  sva = sv_5ph(ia,ja)
  svb = sv_5ph(ib,jb)
  svc = sv_5ph(ic,ja)
  svd = sv_5ph(id,jb)

  svac = sva+(svc-sva)/(ec-ea)*(e-ea)
  svbd = svb+(svd-svb)/(ed-eb)*(e-eb)

  sv = svac+(svbd-svac)/(db-da)*(d-da)

  if (sv < SV_MIN_5PHASE) sv = SV_MIN_5PHASE

  t=0.0e0_8
  s=0.0e0_8
  if(flag_out == 1) then
     !--- Calculate Temperature ---------------------
     ta = t_5ph(ia)
     tb = t_5ph(ib)
     tc = t_5ph(ic)
     td = t_5ph(id)

     tac = ta+(tc-ta)/(ec-ea)*(e-ea)
     tbd = tb+(td-tb)/(ed-eb)*(e-eb)

     t = tac+(tbd-tac)/(db-da)*(d-da)

     !--- Calculate Entropy -------------------------
     sa = s_5ph(ia,ja)
     sb = s_5ph(ib,jb)
     sc = s_5ph(ic,ja)
     sd = s_5ph(id,jb)

     sac = sa+(sc-sa)/(ec-ea)*(e-ea)
     sbd = sb+(sd-sb)/(ed-eb)*(e-eb)

     s = sac+(sbd-sac)/(db-da)*(d-da)
  endif

  return
end subroutine sub_eos_5phase