!!****f* source/materials/eos/helmholtz_threadsafe/multigamma/eos
!!
!! NAME 
!!
!!  eos
!!
!! SYNOPSIS
!!  
!!  subroutine eos(dens, temp, pres, ener, xn, abar, zbar,  
!!                 dpt, dpd, det, ded, c_v, c_p, gammac, pel, ne, eta,  
!!                 input)
!!
!! DESCRIPTION
!!
!!   This is a generic wrapper for the equation of state so the
!!   pressure, energy, etc. can be computed through an argument
!!   call, without requiring the calling routine to share a common
!!   block with the EOS.
!!
!!  This version wraps the Helmholtz equation of state.
!!
!!  ARGUMENTS
!!      dens   : the density of the material
!!      temp   : the temperature 
!!      pres   : pressure
!!      ener   : internal energy
!!      xn()   : mass fraction array
!!      abar   : average mass of the nuclei
!!      zbar   : average proton number
!!      dpt    : d(pres)/d(temp)
!!      dpd    : d(pres)/d(dens)
!!      det    : d(ener)/d(temp)
!!      ded    : d(ener)/d(dens)
!!      c_v    : specific heat at constant volume
!!      c_p    : specific heat at constant pressure
!!      gammac : d(log P)/d(log dens)
!!      pel    : electron pressure
!!      ne     : electron number density
!!      eta    : electron degeneracy parameter (chem pot / k_b*T)
!!      input  : what to do
!!
!!     If input=1, density and temperature are taken as inputs, and
!!     pressure and energy are generated as output.  If input=2,
!!     density and energy are taken as inputs, and pressure and
!!     temperature are generated as output.  Mass fractions are
!!     always taken as input; the remaining parameters are computed
!!     on output.
!!
!!
!! USED BY
!!
!!   init_block
!!   init_1d 
!!   tot_bnd
!!   user_bnd
!!   conductivity
!!   net
!!   perturbLib
!!
!!
!!***

module ModuleEos

  private
  public :: eos

contains

  subroutine eos (dens, temp, pres, ener, xn, abar, zbar, & 
     &                dpt, dpd, det, ded, c_v, c_p, gammac, pel, ne, eta, & 
     &                input)

!===============================================================================

       use logfile, ONLY: stamp_logfile
       use physical_constants, ONLY: get_constant_from_db
       use multifluid_database, ONLY: get_mfluid_property, & 
     &                                query_mfluid_suminv, & 
     &                                query_mfluid_sumfrc, & 
     &                                mfluid_status_fail,  & 
     &                                mf_prop_a, & 
     &                                mf_prop_z
       use dBase, ONLY: ionmax
       implicit none

!       Parameters

      real              dens, temp, pres, ener, dpt, dpd, det, ded, c_v, c_p,  & 
     &                  gammac, pel, ne, eta, abar, zbar
      real              xn(ionmax)
      integer           input

      integer           FAIL, ierr
      parameter         (FAIL = -1)
 
      real              gc(ionmax), gammam1j(ionmax), ggprodj(ionmax)
      real              ggprodinvj(ionmax), gam1invj(ionmax)
      real              weight(ionmax), rt
      real              gascon
      logical ::        first_call = .true.
      save              gascon, first_call
      save              gc, gammam1j,ggprodj,ggprodinvj,gam1invj

!===============================================================================

!       Initialization

      if (first_call) then

        first_call = .false.
        call get_constant_from_db ("ideal gas constant", gascon)

        call get_mfluid_property ("adiabatic index", gc, ierr)
        if (ierr == mfluid_status_fail) then
          call abort_flash ("eos1d:  trouble getting gammas")
        endif
 
        gammam1j   = 1. / (gc - 1.)
        ggprodj    = gammam1j * gascon
        ggprodinvj = 1. / ggprodj
        gam1invj   = 1. / gammam1j

      endif

!       Get avg atomic mass

      call query_mfluid_suminv (mf_prop_A, xn(:), abar)
      abar = 1.e0 / abar
      call query_mfluid_sumfrc (mf_prop_Z, xn(:), zbar)
      zbar = abar * zbar

!       Get weighted gamma

      weight = xn * gammam1j
      call query_mfluid_suminv (mf_prop_A, weight, rt)
      gammac = 1.0e0 + 1.0e0/(rt*abar)

!===============================================================================

!               input = 1:  temperature and density given

      if (input == 1) then


        pres   = gascon * dens * temp / abar
        ener   = gascon * temp / (abar*(gammac-1.))
        dpt    = gascon*dens/abar
        dpd    = gascon*temp/abar
        det    = gascon/(abar*(gammac-1.0))
        ded    = 0.
        pel    = 0.
        ne     = 0.
        eta    = 0.

        c_v = det
        c_p = gammac*c_v

!===============================================================================

!               input = 2:  internal energy and density given

      elseif (input == 2) then

        pres   = dens * ener * (gammac - 1.)
        temp   = ener * (gammac - 1.) * abar/gascon
        dpt    = gascon*dens/abar
        dpd    = gascon*temp/abar
        det    = gascon/(abar*(gammac-1.0))
        ded    = 0.
        pel    = 0.
        ne     = 0.
        eta    = 0.

        c_v = det
        c_p = gammac*c_v

!===============================================================================

!               input = 3:  pressure and density given

      elseif (input == 3) then

        ener   = pres /( (gammac - 1.) * dens )
        temp   = ener * (gammac - 1.) * abar/gascon
        dpt    = gascon*dens/abar
        dpd    = gascon*temp/abar
        det    = gascon/(abar*(gammac-1.0))
        ded    = 0.
        pel    = 0.
        ne     = 0.
        eta    = 0.

        c_v = det
        c_p = gammac*c_v

!===============================================================================

!               unrecognized value for input

      else

        call stamp_logfile ("eos:  unrecognized value for input")
        call abort_flash("eos: unrecognized value for input")
        stop

      endif

!===============================================================================

      return
    end subroutine eos

end module ModuleEos