!====================================================================
!   SUBROUTINE   ---- sub_eps_saumon ----                         ===
!                                                                 ===
!   Subroutine for calculating the pressure                       ===
!   and sound velocity from the Saumon EOS                        ===
!                                                      2006/06/18 ===
!                                                  mod 2012/09/12 ===
!====================================================================
!====================================================================
!   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_saumon(d,e,p,sv,t)
  use mod_parameter_saumon_eos
  implicit none

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

!---- Variables Used in This Subroutine -----------------------------
  real(8) :: sv2     ! Square of sv
  real(8) :: d_log
  real(8) :: e_log
  integer(4) :: id,ie
  real(8) :: dx,dy
  real(8) :: p_log_Saumon_a,t_log_Saumon_a
  real(8) :: p_log_Saumon_b,t_log_Saumon_b
  real(8) :: p_log_Saumon_c,t_log_Saumon_c
  real(8) :: p_log_Saumon_d,t_log_Saumon_d
  real(8) :: p_log_ab,t_log_ab
  real(8) :: p_log_cd,t_log_cd
  real(8) :: p_log_ac,t_log_ac
  real(8) :: p_log_bd,t_log_bd
  real(8) :: p_log,t_log
  real(8) :: dp_log_dd_log
  real(8) :: dp_log_de_log
!--------------------------------------------------------------------

  if (e .le. 0.0d0) then
     p=0.0d0
     t=0.0d0
     sv=1.0d-10
     write(60,'(a7)') '   e=<0'
     go to 555
  end if

  d_log=log10(d)
  e_log=log10(e)

  id=int((d_log-d_log_min_Saumon)/dd_log_Saumon)
  ie=int((e_log-e_log_min_Saumon)/de_log_Saumon)

  if (id .lt. 0 .or. id .ge. id_Saumon) then
     p=0.0d0
     t=0.0d0
     sv=1.0d-10
     write(60,*) id,d,d_log,d_log_min_Saumon,dd_log_Saumon,'   d:out'
     go to 555
  end if

  if (ie .lt. 0 .or. ie .ge. ie_Saumon) then
     p=0.0d0
     t=0.0d0
     sv=1.0d-10
     write(60,*) e,'   e:out'
     go to 555
  end if

  dx=(d_log-d_log_Saumon(id))/dd_log_Saumon
  dy=(e_log-e_log_Saumon(ie))/de_log_Saumon

  p_log_Saumon_a=p_log_Saumon(id,ie)
  p_log_Saumon_b=p_log_Saumon(id+1,ie)
  p_log_Saumon_c=p_log_Saumon(id,ie+1)
  p_log_Saumon_d=p_log_Saumon(id+1,ie+1)

  t_log_Saumon_a=t_log_Saumon(id,ie)
  t_log_Saumon_b=t_log_Saumon(id+1,ie)
  t_log_Saumon_c=t_log_Saumon(id,ie+1)
  t_log_Saumon_d=t_log_Saumon(id+1,ie+1)

  if (p_log_Saumon_a .eq. p_log_out_Saumon .or.          &
  &    p_log_Saumon_b .eq. p_log_out_Saumon .or.         &
  &    p_log_Saumon_c .eq. p_log_out_Saumon .or.         &
  &    p_log_Saumon_d .eq. p_log_out_Saumon) then
     p=0.0d0
     t=0.0d0
     sv=1.0d-10
     write(60,'(a9)') '   NoDATA'

!cccc
!         write(55,'(2i7)') id,ie
!         write(55,'(2d15.5)') d_log,e_log
!         write(55,'(4d15.5)') p_log_Saumon_a,p_log_Saumon_b,
!     $                p_log_Saumon_c,p_log_Saumon_d
!cccc
     go to 555
  end if

  p_log_ab=p_log_Saumon_a+(p_log_Saumon_b-p_log_Saumon_a)*dx
  p_log_cd=p_log_Saumon_c+(p_log_Saumon_d-p_log_Saumon_c)*dx
  p_log_ac=p_log_Saumon_a+(p_log_Saumon_c-p_log_Saumon_a)*dy
  p_log_bd=p_log_Saumon_b+(p_log_Saumon_d-p_log_Saumon_b)*dy

  t_log_ab=t_log_Saumon_a+(t_log_Saumon_b-t_log_Saumon_a)*dx
  t_log_cd=t_log_Saumon_c+(t_log_Saumon_d-t_log_Saumon_c)*dx
  t_log_ac=t_log_Saumon_a+(t_log_Saumon_c-t_log_Saumon_a)*dy
  t_log_bd=t_log_Saumon_b+(t_log_Saumon_d-t_log_Saumon_b)*dy

  p_log=p_log_ab+(p_log_cd-p_log_ab)*dy
  t_log=t_log_ab+(t_log_cd-t_log_ab)*dy

  p=10.0d0**p_log
  t=10.0d0**t_log

  dp_log_dd_log=(p_log_bd-p_log_ac)/dd_log_Saumon
  dp_log_de_log=(p_log_cd-p_log_ab)/de_log_Saumon

  sv2=p/d*dp_log_dd_log+p*p/d/d/e*dp_log_de_log
  sv=sqrt(sv2)

!  write(55,'(d15.5)') sv

555 continue

  return
end subroutine sub_eos_saumon