anchor.diff:  
          rho = p/e/(gamma - 1.0_dp_t)
   
  +       if (present(T)) T = mu*m_amu*e/(1.5_dp_t*k)
  +       if (present(s)) s = c_v*log(p/rho**gamma) 
  +
  +    else if (input == eos_input_ps) then
  +
  +       ! entropy and pressure and inputs, return density and e 
  +       if (.not. present(s)) then
  +          print *, 'ERROR: entropy must be input'
  +          stop
  +       endif
  +       
  +       ! note: our definition of entropy neglects the constant since
  +       ! we really only care about changes in entropy.
  +       rho = (p*exp(-s/c_v))**(1.0_dp_t/gamma)
  +
  +       e = p / (rho * (gamma - 1.0_dp_t))
  +
  +       if (present(T)) T = mu*m_amu*e/(1.5_dp_t*k)
  +
       else
          print *, 'ERROR: invalid input in EOS'
          stop
       endif
   
  +
     end subroutine eos
   
   end module eos_module
  
positive.code: ! a simple gamma-law equation of state
  ! p = rho e (gamma - 1)
  module eos_module
  
    use datatypes_module
    use params_module, only: gamma
    implicit none
  
    integer, parameter :: eos_input_p = 1
    integer, parameter :: eos_input_e = 2
    integer, parameter :: eos_input_pe = 3
  
    private
  
    public :: eos_input_p, eos_input_e, eos_input_pe
    public :: eos
  
  contains
  
    subroutine eos(input, p, e, rho)
      
      integer,          intent(in   ) :: input
      real (kind=dp_t), intent(inout) :: p
      real (kind=dp_t), intent(inout) :: e
      real (kind=dp_t), intent(inout) :: rho
  
      if (input == eos_input_p) then
  
  
negative.code: 
         e = p / (rho * (gamma - 1.0_dp_t))
  
         if (present(T)) T = mu*m_amu*e/(1.5_dp_t*k)
         if (present(s)) s = c_v*log(p/rho**gamma) 
  
      else if (input == eos_input_e) then
  
         if (e < 0.0_dp_t) then
  