anchor.diff:  
     integer, parameter :: eos_input_p = 1
     integer, parameter :: eos_input_e = 2
  +  integer, parameter :: eos_input_pe = 3
   
   contains
   
  
positive.code:     integer,          intent(in   ) :: input
      real (kind=dp_t), intent(inout) :: p
      real (kind=dp_t), intent(inout) :: e
      real (kind=dp_t), intent(in   ) :: rho
  
      if (input == eos_input_p) then
  
  
negative.code: 
         p = rho*e*(gamma - 1.0_dp_t)
  
      else if (input == eos_input_pe) then
  
         ! pressure and energy are inputs -- return density
         if (e < 0.0_dp_t .or. p < 0.0_dp_t) then
            print *, 'ERROR: input internal energy or pressure < 0 in EOS'
            stop
         endif
  
         rho = p/e/(gamma - 1.0_dp_t)
  
      else
         print *, 'ERROR: invalid input in EOS'
         stop
      endif
  
    end subroutine eos
  