RIT MHD dev:

Fix to GRHydro_RegisterVarsM.cc and allows MHD to work without Y_E bing evolved.

Add C2P polytype interface.

Change schedule.ccl to represent a more detailed version of choosing what 
gridfucntions get sync'd in Boundaries depending on what is being evolved.

Un-hardwire fixed values of gamma for all MHD Con2Prim routines.




git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@199 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

6 files changed, 411 insertions, 193 deletions

diff --git a/schedule.ccl b/schedule.ccl
index 1f79f13..5fec310 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -385,21 +385,39 @@ if (EoS_Change)
   }
   if (CCTK_Equals(EoS_Change_type,"GammaKS"))
   {
-    schedule GRHydro_EoSChangeGammaK_Shibata AT CCTK_Initial AFTER HydroBase_Initial BEFORE GRHydro_IVP
-    { 
-      LANG: Fortran
-      SYNC: dens
-      SYNC: tau
-      SYNC: scon
-      SYNC: w_lorentz
-      SYNC: hydrobase::rho
-      SYNC: hydrobase::press
-      SYNC: hydrobase::eps
-      SYNC: hydrobase::vel
-      SYNC: hydrobase::temperature
-      SYNC: hydrobase::Y_e
-      SYNC: Y_e_con
-    } "Reset the hydro variables if the EoS Gamma and K change between ID and evolution"
+
+    if(CCTK_Equals(Y_e_evolution_method,"GRHydro")) {
+      schedule GRHydro_EoSChangeGammaK_Shibata AT CCTK_Initial AFTER HydroBase_Initial BEFORE GRHydro_IVP
+      { 
+        LANG: Fortran
+        SYNC: dens
+        SYNC: tau
+        SYNC: scon
+        SYNC: w_lorentz
+        SYNC: hydrobase::rho
+        SYNC: hydrobase::press
+        SYNC: hydrobase::eps
+        SYNC: hydrobase::vel
+        SYNC: hydrobase::temperature
+        SYNC: hydrobase::Y_e
+        SYNC: Y_e_con
+      } "Reset the hydro variables if the EoS Gamma and K change between ID and evolution"
+    } else {
+      schedule GRHydro_EoSChangeGammaK_Shibata AT CCTK_Initial AFTER HydroBase_Initial BEFORE GRHydro_IVP
+      { 
+        LANG: Fortran
+        SYNC: dens
+        SYNC: tau
+        SYNC: scon
+        SYNC: w_lorentz
+        SYNC: hydrobase::rho
+        SYNC: hydrobase::press
+        SYNC: hydrobase::eps
+        SYNC: hydrobase::vel
+        SYNC: hydrobase::temperature
+        SYNC: Y_e_con
+      } "Reset the hydro variables if the EoS Gamma and K change between ID and evolution"
+    }
   }
 }
 
@@ -1072,91 +1090,167 @@ schedule group Do_GRHydro_Boundaries IN HydroBase_Boundaries
 {
 } "GRHydro Boundary conditions group"
 
-
-if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
-{
-  if (evolve_tracer)
+if(CCTK_Equals(Y_e_evolution_method,"GRHydro")) {
+  if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
   {
-    schedule GRHydro_BoundariesM IN HydroBase_Select_Boundaries AS GRHydro_Bound
+    if (evolve_tracer)
     {
-      LANG: Fortran
-      OPTIONS: LEVEL
-      SYNC: dens
-      SYNC: tau
-      SYNC: scon
-      SYNC: w_lorentz
-      SYNC: HydroBase::rho
-      SYNC: HydroBase::press
-      SYNC: HydroBase::eps
-      SYNC: HydroBase::vel
-      SYNC: HydroBase::Bvec
-      SYNC: GRHydro_cons_tracers
-      SYNC: GRHydro_tracers
-      SYNC: hydrobase::temperature
-      SYNC: hydrobase::Y_e
-      SYNC: Y_e_con
-    } "Select GRHydro boundary conditions - MHD version"
+      schedule GRHydro_BoundariesM IN HydroBase_Select_Boundaries AS GRHydro_Bound
+      {
+        LANG: Fortran
+        OPTIONS: LEVEL
+        SYNC: dens
+        SYNC: tau
+        SYNC: scon
+        SYNC: w_lorentz
+        SYNC: HydroBase::rho
+        SYNC: HydroBase::press
+        SYNC: HydroBase::eps
+        SYNC: HydroBase::vel
+        SYNC: HydroBase::Bvec
+        SYNC: GRHydro_cons_tracers
+        SYNC: GRHydro_tracers
+        SYNC: hydrobase::temperature
+        SYNC: hydrobase::Y_e
+        SYNC: Y_e_con
+      } "Select GRHydro boundary conditions - MHD version"
+    } else {
+      schedule GRHydro_BoundariesM IN HydroBase_Select_Boundaries AS GRHydro_Bound
+      {
+        LANG: Fortran
+        OPTIONS: LEVEL
+        SYNC: dens
+        SYNC: tau
+        SYNC: scon
+        SYNC: w_lorentz
+        SYNC: HydroBase::rho
+        SYNC: HydroBase::press
+        SYNC: HydroBase::eps
+        SYNC: HydroBase::vel
+        SYNC: HydroBase::Bvec
+        SYNC: hydrobase::temperature
+        SYNC: hydrobase::Y_e
+        SYNC: Y_e_con
+      } "Select GRHydro boundary conditions - MHD version"
+    }
   } else {
-    schedule GRHydro_BoundariesM IN HydroBase_Select_Boundaries AS GRHydro_Bound
+    if (evolve_tracer)
     {
-      LANG: Fortran
-      OPTIONS: LEVEL
-      SYNC: dens
-      SYNC: tau
-      SYNC: scon
-      SYNC: w_lorentz
-      SYNC: HydroBase::rho
-      SYNC: HydroBase::press
-      SYNC: HydroBase::eps
-      SYNC: HydroBase::vel
-      SYNC: HydroBase::Bvec
-      SYNC: hydrobase::temperature
-      SYNC: hydrobase::Y_e
-      SYNC: Y_e_con
-    } "Select GRHydro boundary conditions - MHD version"
+      schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound
+      {
+        LANG: Fortran
+        OPTIONS: LEVEL
+        SYNC: dens
+        SYNC: tau
+        SYNC: scon
+        SYNC: w_lorentz
+        SYNC: HydroBase::rho
+        SYNC: HydroBase::press
+        SYNC: HydroBase::eps
+        SYNC: HydroBase::vel
+        SYNC: GRHydro_cons_tracers
+        SYNC: GRHydro_tracers
+        SYNC: hydrobase::temperature
+        SYNC: hydrobase::Y_e
+        SYNC: Y_e_con
+      } "Select GRHydro boundary conditions"
+    } else
+    {
+      schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound
+      {
+        LANG: Fortran
+        OPTIONS: LEVEL
+        SYNC: dens
+        SYNC: tau
+        SYNC: scon
+        SYNC: w_lorentz
+        SYNC: HydroBase::rho
+        SYNC: HydroBase::press
+        SYNC: HydroBase::eps
+        SYNC: HydroBase::vel
+        SYNC: hydrobase::temperature
+        SYNC: hydrobase::Y_e
+        SYNC: Y_e_con
+      } "Select GRHydro boundary conditions"
+    }
   }
 } else {
-  if (evolve_tracer)
+  if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
   {
-    schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound
+    if (evolve_tracer)
     {
-      LANG: Fortran
-      OPTIONS: LEVEL
-      SYNC: dens
-      SYNC: tau
-      SYNC: scon
-      SYNC: w_lorentz
-      SYNC: HydroBase::rho
-      SYNC: HydroBase::press
-      SYNC: HydroBase::eps
-      SYNC: HydroBase::vel
-      SYNC: GRHydro_cons_tracers
-      SYNC: GRHydro_tracers
-      SYNC: hydrobase::temperature
-      SYNC: hydrobase::Y_e
-      SYNC: Y_e_con
-    } "Select GRHydro boundary conditions"
-  } else
-  {
-    schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound
+      schedule GRHydro_BoundariesM IN HydroBase_Select_Boundaries AS GRHydro_Bound
+      {
+        LANG: Fortran
+        OPTIONS: LEVEL
+        SYNC: dens
+        SYNC: tau
+        SYNC: scon
+        SYNC: w_lorentz
+        SYNC: HydroBase::rho
+        SYNC: HydroBase::press
+        SYNC: HydroBase::eps
+        SYNC: HydroBase::vel
+        SYNC: HydroBase::Bvec
+        SYNC: GRHydro_cons_tracers
+        SYNC: GRHydro_tracers
+        SYNC: hydrobase::temperature
+      } "Select GRHydro boundary conditions - MHD version"
+    } else {
+      schedule GRHydro_BoundariesM IN HydroBase_Select_Boundaries AS GRHydro_Bound
+      {
+        LANG: Fortran
+        OPTIONS: LEVEL
+        SYNC: dens
+        SYNC: tau
+        SYNC: scon
+        SYNC: w_lorentz
+        SYNC: HydroBase::rho
+        SYNC: HydroBase::press
+        SYNC: HydroBase::eps
+        SYNC: HydroBase::vel
+        SYNC: HydroBase::Bvec
+        SYNC: hydrobase::temperature
+      } "Select GRHydro boundary conditions - MHD version"
+    }
+  } else {
+    if (evolve_tracer)
     {
-      LANG: Fortran
-      OPTIONS: LEVEL
-      SYNC: dens
-      SYNC: tau
-      SYNC: scon
-      SYNC: w_lorentz
-      SYNC: HydroBase::rho
-      SYNC: HydroBase::press
-      SYNC: HydroBase::eps
-      SYNC: HydroBase::vel
-      SYNC: hydrobase::temperature
-      SYNC: hydrobase::Y_e
-      SYNC: Y_e_con
-    } "Select GRHydro boundary conditions"
+      schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound
+      {
+        LANG: Fortran
+        OPTIONS: LEVEL
+        SYNC: dens
+        SYNC: tau
+        SYNC: scon
+        SYNC: w_lorentz
+        SYNC: HydroBase::rho
+        SYNC: HydroBase::press
+        SYNC: HydroBase::eps
+        SYNC: HydroBase::vel
+        SYNC: GRHydro_cons_tracers
+        SYNC: GRHydro_tracers
+        SYNC: hydrobase::temperature
+      } "Select GRHydro boundary conditions"
+    } else
+    {
+      schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound
+      {
+        LANG: Fortran
+        OPTIONS: LEVEL
+        SYNC: dens
+        SYNC: tau
+        SYNC: scon
+        SYNC: w_lorentz
+        SYNC: HydroBase::rho
+        SYNC: HydroBase::press
+        SYNC: HydroBase::eps
+        SYNC: HydroBase::vel
+        SYNC: hydrobase::temperature
+      } "Select GRHydro boundary conditions"
+    }
   }
-}
-
+}   
 ############################################################
 ### Compute first differences of rho for mesh refinement ###
 ############################################################
diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90
index bbb803d..984569a 100644
--- a/src/GRHydro_Con2PrimM.F90
+++ b/src/GRHydro_Con2PrimM.F90
@@ -58,11 +58,11 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
   CCTK_INT :: type_bits, atmosphere
   CCTK_INT :: type2_bits
 
-  CCTK_REAL :: local_min_tracer
+  CCTK_REAL :: local_min_tracer, local_gam, local_pgam,local_K,sc
 
 ! begin EOS Omni vars                                                                                       
   integer :: n,keytemp,anyerr,keyerr(1)
-  real*8  :: xpress,xtemp,xye,xeps
+  real*8  :: xpress,xtemp,xye,xeps,xrho
   n=1;keytemp=0;anyerr=0;keyerr(1)=0
   xpress=0.0d0;xtemp=0.0d0;xye=0.0d0;xeps=0.0d0
 ! end EOS Omni vars                 
@@ -87,6 +87,11 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
   call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
 
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+         1.d0,1.d0,xtemp,xye,local_gam,keyerr,anyerr)
+  local_gam = local_gam+1.d0
+  
+
   !$OMP PARALLEL DO PRIVATE(i,j,itracer,&
   !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative)
   do k = 1, nz 
@@ -154,7 +159,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
            end if
            
            if(evolve_temper.eq.0) then
-              call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, dens(i,j,k), &
+              call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam, dens(i,j,k), &
                    scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), &
                    rho(i,j,k),&
                    vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
@@ -219,14 +224,24 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
               call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype.')
               !$OMP END CRITICAL
               
-!!$              call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, dens(i,j,k), &
-!!$                   scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), &
-!!$                   rho(i,j,k),&
-!!$                   vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
-!!$                   gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
-!!$                   uxx,uxy,uxz,uyy,uyz,uzz,det, &
-!!$                   Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
-!!$                   epsnegative,GRHydro_C2P_failed(i,j,k))
+              xrho=10.0d0
+              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   1.d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+              local_K = xpress
+
+              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+              local_pgam=log(xpress/local_K)/log(xrho)
+              sc = local_K*dens(i,j,k)
+
+              call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), &
+                   scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, &
+                   rho(i,j,k),&
+                   vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
+                   gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
+                   uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                   Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
+                   epsnegative,GRHydro_C2P_failed(i,j,k))
               
 #endif          
               
@@ -276,7 +291,7 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
        uxxr, uxyr, uxzr, uyyr, uyzr, uzzr, pmin, epsmin
   CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
        gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
-  CCTK_REAL :: b2minus, b2plus
+  CCTK_REAL :: b2minus, b2plus, local_gam, local_pgam,local_K,sc
   CCTK_INT :: epsnegative
   character(len=100) warnline
  
@@ -285,9 +300,9 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
 
   CCTK_REAL :: local_min_tracer
 
-! begin EOS Omni vars                                                                                       
+! begin EOS Omni vars                                            
   integer :: n,keytemp,anyerr,keyerr(1)
-  real*8  :: xpress,xtemp,xye,xeps
+  real*8  :: xpress,xtemp,xye,xeps,xrho
   n=1;keytemp=0;anyerr=0;keyerr(1)=0
   xpress=0.0d0;xeps=0.0d0;xtemp=0.0d0;xye=0.0d0
 ! end EOS Omni vars                       
@@ -298,7 +313,11 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
 
   call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,epsmin,xtemp,xye,pmin,epsmin,keyerr,anyerr)
-                                                                              
+
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+         1.d0,1.0d0,xtemp,xye,local_gam,keyerr,anyerr)
+  local_gam=local_gam+1.0
+
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
   type2_bits = -1
@@ -363,7 +382,7 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
            endif
            
            
-           call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, densminus(i,j,k), &
+           call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam,densminus(i,j,k), &
                 sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k), &
                 rhominus(i,j,k),&
                 velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),&
@@ -376,14 +395,25 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
               !$OMP CRITICAL
               call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype!')
               !$OMP END CRITICAL
-!!$              call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, densminus(i,j,k), &
-!!$                   sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k), &
-!!$                   rhominus(i,j,k),&
-!!$                   velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),&
-!!$                   gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
-!!$                   uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, &
-!!$                   Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,&
-!!$                   epsnegative,GRHydro_C2P_failed(i,j,k))
+              
+              xrho=10.0d0
+              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   1.d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+              local_K = xpress
+
+              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+              local_pgam=log(xpress/local_K)/log(xrho)
+              sc = local_K*densminus(i,j,k)
+
+              call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, densminus(i,j,k), &
+                   sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), sc, &
+                   rhominus(i,j,k),&
+                   velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),&
+                   gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
+                   uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, &
+                   Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,&
+                   epsnegative,GRHydro_C2P_failed(i,j,k))
            end if
            
            if (epsminus(i,j,k) .lt. 0.0d0) then
@@ -408,7 +438,7 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
            endif
            
            epsnegative = 0
-           call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, densplus(i,j,k), &
+           call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam, densplus(i,j,k), &
                 sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), tauplus(i,j,k),&
                 rhoplus(i,j,k),&
                 velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),&
@@ -421,14 +451,25 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
               !$OMP CRITICAL
               call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype!!')
               !$OMP END CRITICAL
-!!$              call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, densplus(i,j,k), &
-!!$                   sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), tauplus(i,j,k),&
-!!$                   rhoplus(i,j,k),&
-!!$                   velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),&
-!!$                   gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, &
-!!$                   uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, &
-!!$                   Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,&
-!!$                   epsnegative,GRHydro_C2P_failed(i,j,k))
+
+              xrho=10.0d0
+              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   1.d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+              local_K = xpress
+
+              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+              local_pgam=log(xpress/local_K)/log(xrho)
+              sc = local_K*densplus(i,j,k)
+
+              call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, densplus(i,j,k), &
+                   sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), sc,&
+                   rhoplus(i,j,k),&
+                   velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),&
+                   gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, &
+                   uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, &
+                   Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,&
+                   epsnegative,GRHydro_C2P_failed(i,j,k))
            end if
            
            if (epsplus(i,j,k) .lt. 0.0d0) then
@@ -492,9 +533,15 @@ subroutine Conservative2PrimitivePolytypeM(CCTK_ARGUMENTS)
 
   CCTK_INT :: epsnegative
   
-  CCTK_REAL :: local_min_tracer
+  CCTK_REAL :: local_min_tracer, local_gam, local_pgam,local_K, sc
   !  character(len=400) :: warnline
   
+! begin EOS Omni vars                                                                                       
+  integer :: n,keytemp,anyerr,keyerr(1)
+  real*8  :: xpress,xtemp,xye,xeps,xrho
+  n=1;keytemp=0;anyerr=0;keyerr(1)=0
+  xpress=0.0d0;xtemp=0.0d0;xye=0.0d0;xeps=0.0d0
+! end EOS Omni vars                 
   
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
@@ -547,15 +594,24 @@ subroutine Conservative2PrimitivePolytypeM(CCTK_ARGUMENTS)
               Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
            endif
            
-
-!!$           call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, dens(i,j,k), &
-!!$                scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k),&
-!!$                rho(i,j,k),&
-!!$                vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
-!!$                gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
-!!$                uxx,uxy,uxz,uyy,uyz,uzz,det, &
-!!$                Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
-!!$                epsnegative,GRHydro_C2P_failed(i,j,k))
+           xrho=10.0d0
+           call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                1.d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+           local_K = xpress
+           
+           call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+           local_pgam=log(xpress/local_K)/log(xrho)
+           sc = local_K*dens(i,j,k)
+           
+           call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), &
+                scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc,&
+                rho(i,j,k),&
+                vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
+                gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
+                uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
+                epsnegative,GRHydro_C2P_failed(i,j,k))
 
         end do
      end do
@@ -603,6 +659,15 @@ subroutine Con2PrimBoundsPolytypeM(CCTK_ARGUMENTS)
 
   CCTK_INT :: type_bits, atmosphere
   CCTK_INT :: type2_bits, epsnegative
+
+  CCTK_REAL :: local_gam, local_pgam,local_K,scplus,scminus
+
+! begin EOS Omni vars                                                                                       
+  integer :: n,keytemp,anyerr,keyerr(1)
+  real*8  :: xpress,xtemp,xye,xeps,xrho
+  n=1;keytemp=0;anyerr=0;keyerr(1)=0
+  xpress=0.0d0;xtemp=0.0d0;xye=0.0d0;xeps=0.0d0
+! end EOS Omni vars                 
   
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
@@ -640,22 +705,34 @@ subroutine Con2PrimBoundsPolytypeM(CCTK_ARGUMENTS)
            call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,&
                 gxxr,gxyr,gxzr,gyyr,gyzr,gzzr)        
 
-!!$           call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, densminus(i,j,k), &
-!!$                sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k),&
-!!$                rhominus(i,j,k),&
-!!$                velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),&
-!!$                gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
-!!$                uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, &
-!!$                Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,&
-!!$                epsnegative,GRHydro_C2P_failed(i,j,k))
-!!$           call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, densplus(i,j,k), &
-!!$                sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), tauplus(i,j,k),&
-!!$                rhoplus(i,j,k),&
-!!$                velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),&
-!!$                gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, &
-!!$                uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, &
-!!$                Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,&
-!!$                epsnegative,GRHydro_C2P_failed(i,j,k))
+           xrho=10.0d0
+           call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                1.d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+           local_K = xpress
+           
+           call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+           local_pgam=log(xpress/local_K)/log(xrho)
+           scminus = local_K*densminus(i,j,k)
+           scplus = local_K*densplus(i,j,k)
+           
+           call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densminus(i,j,k), &
+                sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), scminus,&
+                rhominus(i,j,k),&
+                velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),&
+                gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
+                uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, &
+                Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,&
+                epsnegative,GRHydro_C2P_failed(i,j,k))
+
+           call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densplus(i,j,k), &
+                sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), scplus,&
+                rhoplus(i,j,k),&
+                velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),&
+                gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, &
+                uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, &
+                Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,&
+                epsnegative,GRHydro_C2P_failed(i,j,k))
         end do
      end do
   end do
@@ -682,6 +759,7 @@ end subroutine Con2PrimBoundsPolytypeM
 subroutine Conservative2PrimitiveGeneralM(CCTK_ARGUMENTS)
   
   USE GRHydro_Scalars
+  USE Con2PrimM_fortran_interfaces
   
   implicit none
   
@@ -699,7 +777,7 @@ subroutine Conservative2PrimitiveGeneralM(CCTK_ARGUMENTS)
 
   CCTK_INT :: count
   CCTK_REAL :: b2,det
-  CCTK_REAL :: uxx,uxy,uxz,uyy,uyz,uzz
+  CCTK_REAL :: uxx,uxy,uxz,uyy,uyz,uzz, local_gam, local_pgam
 
   CCTK_INT :: type_bits, atmosphere
   CCTK_INT :: type2_bits,epsnegative
@@ -707,6 +785,24 @@ subroutine Conservative2PrimitiveGeneralM(CCTK_ARGUMENTS)
   CCTK_REAL :: local_min_tracer
 
   CCTK_INT, dimension(3) :: loc, loc2
+  
+!!$ begin EOS Omni vars
+  integer :: n,keytemp,anyerr,keyerr(1)
+  real*8  :: xpress(1),xtemp(1),xye(1),xdens(1),xeps(1)
+  n=1 
+  keytemp=0
+  anyerr=0
+  keyerr(1)=0
+  xpress(1)=0.0d0
+  xtemp(1)=0.0d0
+  xye(1)=0.0d0
+!!$ end EOS Omni vars                 
+
+  xdens(1)=1.d0
+  xeps(1)=1.d0
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+         xdens,xeps,xtemp,xye,xpress,keyerr,anyerr)
+  local_gam=xpress(1)+1.d0
 
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
@@ -776,7 +872,7 @@ subroutine Conservative2PrimitiveGeneralM(CCTK_ARGUMENTS)
                 gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
                 gyz(i,j,k),gzz(i,j,k))        
 
-           call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, dens(i,j,k), &
+           call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam,dens(i,j,k), &
                 scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k),&
                 rho(i,j,k),&
                 vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
@@ -828,7 +924,14 @@ subroutine Con2PrimPolytypeGeneralM(CCTK_ARGUMENTS)
   CCTK_INT :: atmosphere
   CCTK_INT :: type2_bits,epsnegative
 
-  CCTK_REAL :: local_min_tracer
+  CCTK_REAL :: local_min_tracer, local_pgam, local_k, sc
+
+! begin EOS Omni vars                                                                                       
+  integer :: n,keytemp,anyerr,keyerr(1)
+  real*8  :: xpress(1),xtemp(1),xye(1),xeps(1),xrho(1)
+  n=1;keytemp=0;anyerr=0;keyerr(1)=0
+  xpress(1)=0.0d0;xtemp(1)=0.0d0;xye(1)=0.0d0;xeps(1)=0.0d0
+! end EOS Omni vars                 
                                                                                
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
@@ -894,15 +997,26 @@ subroutine Con2PrimPolytypeGeneralM(CCTK_ARGUMENTS)
         call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,&
              gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
              gyz(i,j,k),gzz(i,j,k))        
+
+        xrho(1)=1.d0
+        call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+             xrho,xeps,xtemp,xye,xpress,keyerr,anyerr)
+        local_K = xpress(1)
+        
+        xrho(1)=10.d0
+        call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+             xrho,xeps,xtemp,xye,xpress,keyerr,anyerr)
+        local_pgam=log(xpress(1)/local_K)/log(xrho(1))
+        sc = local_K*dens(i,j,k)
         
-!!$        call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, dens(i,j,k), &
-!!$             scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k),&
-!!$             rho(i,j,k),&
-!!$             vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
-!!$             gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
-!!$             uxx,uxy,uxz,uyy,uyz,uzz,det, &
-!!$             Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
-!!$             epsnegative,GRHydro_C2P_failed(i,j,k))
+        call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,dens(i,j,k), &
+             scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc,&
+             rho(i,j,k),&
+             vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
+             gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
+             uxx,uxy,uxz,uyy,uyz,uzz,det, &
+             Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
+             epsnegative,GRHydro_C2P_failed(i,j,k))
 
      end do
   end do
diff --git a/src/GRHydro_Con2PrimM_pt.c b/src/GRHydro_Con2PrimM_pt.c
index d9ba2a6..bf888ff 100644
--- a/src/GRHydro_Con2PrimM_pt.c
+++ b/src/GRHydro_Con2PrimM_pt.c
@@ -76,36 +76,33 @@ static CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq ;
 // Declarations: 
 static CCTK_REAL vsq_calc(CCTK_REAL W);
 
-static CCTK_INT twod_newton_raphson( CCTK_REAL x[], 
+static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos,
 				     void (*funcd) (CCTK_REAL [], CCTK_REAL [], 
 						    CCTK_REAL [], CCTK_REAL [][2], 
-						    CCTK_REAL *, CCTK_REAL *) );
+						    CCTK_REAL *, CCTK_REAL *, CCTK_REAL) );
 
 static  void func_vsq( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][2], 
-		       CCTK_REAL *f, CCTK_REAL *df);
+		       CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos);
 static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) ;
 
 
 // EOS STUFF:
-#define gam (5./3.)
-#define GAMMA	        (gam)  
-static CCTK_REAL gam_m1_o_gam = ((GAMMA-1.)/GAMMA);
-static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq);
+static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos);
 /* pressure as a function of rho0 and u */
-static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u)
+static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u, CCTK_REAL gammaeos)
 {
-  return((GAMMA - 1.)*u) ;
+  return((gammaeos - 1.)*u) ;
 }
 
 /* Pressure as a function of rho0 and w = rho0 + u + p */
-static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w)
+static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w,CCTK_REAL gammaeos)
 {
-  return((GAMMA-1.)*(w - rho0)/GAMMA) ;
+  return((gammaeos-1.)*(w - rho0)/gammaeos) ;
 }
 
 
 void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) (
-    CCTK_INT  *handle, 
+    CCTK_INT  *handle, CCTK_REAL *gamma_eos,
     CCTK_REAL *dens_in, 
     CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, 
     CCTK_REAL *tau_in, 
@@ -186,7 +183,7 @@ RETURN VALUE: of retval = (i*100 + j)  where
 
 **********************************************************************************/
 void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  ( 
-    CCTK_INT  *handle, 
+    CCTK_INT  *handle, CCTK_REAL *gamma_eos,
     CCTK_REAL *dens_in, 
     CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, 
     CCTK_REAL *tau_in, 
@@ -208,7 +205,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   CCTK_REAL x_2d[2];
   CCTK_REAL sx, sy, sz;
   CCTK_REAL usx, usy, usz;
-  CCTK_REAL tau, dens;
+  CCTK_REAL tau, dens, gammaeos;
   CCTK_REAL QdotB;
   CCTK_REAL rho0,u,p,w,gammasq,gamma,gtmp,W_last,W,vsq;
   CCTK_REAL g_o_WBsq, QdB_o_W;
@@ -217,6 +214,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; 
   CCTK_INT i,j, i_increase ;
 
+  gammaeos = *gamma_eos;
+
   /* Assume ok initially: */
   *retval = 0.;
   *epsnegative = 0; 
@@ -314,7 +313,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   //   i.e. you don't get positive values for dP/d(vsq) . 
   rho0 = D / gamma ;
   u = (*epsilon) * rho0;
-  p = pressure_rho0_u(rho0,u) ;  // EOS
+  p = pressure_rho0_u(rho0,u,gammaeos) ;  // EOS
   w = rho0 + u + p ;
 
   W_last = w*gammasq ;
@@ -332,7 +331,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   // Calculate W and vsq: 
   x_2d[0] =  fabs( W_last );
   x_2d[1] = x1_of_x0( W_last ) ;
-  *retval = 1.0*twod_newton_raphson( x_2d, func_vsq ) ;  
+  *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, func_vsq ) ;  
 
   W = x_2d[0];
   vsq = x_2d[1];
@@ -361,7 +360,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   rho0 = D * gtmp;
 
   w = W * (1. - vsq) ;
-  p = pressure_rho0_w(rho0,w) ;  // EOS
+  p = pressure_rho0_w(rho0,w,gammaeos) ;  // EOS
   u = w - (rho0 + p) ;
 
   // User may want to handle this case differently, e.g. do NOT return upon 
@@ -479,10 +478,10 @@ static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] )
     -- inspired in part by Num. Rec.'s routine newt();
 
 *****************************************************************/
-static CCTK_INT twod_newton_raphson( CCTK_REAL x[], 
+static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, 
 			    void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], 
 					   CCTK_REAL [][2], CCTK_REAL *, 
-					   CCTK_REAL *) )
+					   CCTK_REAL *, CCTK_REAL) )
 {
   CCTK_REAL f, df, dx[2], x_old[2];
   CCTK_REAL resid[2], jac[2][2];
@@ -508,7 +507,7 @@ static CCTK_INT twod_newton_raphson( CCTK_REAL x[],
   keep_iterating = 1;
   while( keep_iterating ) { 
 
-    (*funcd) (x, dx, resid, jac, &f, &df);  /* returns with new dx, f, df */
+    (*funcd) (x, dx, resid, jac, &f, &df, gammaeos);  /* returns with new dx, f, df */
       
 
     /* Save old values before calculating the new: */
@@ -600,7 +599,7 @@ static CCTK_INT twod_newton_raphson( CCTK_REAL x[],
  *********************************************************************************/
 
 static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], 
-		     CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df)
+		     CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos)
 {
 
   
@@ -613,7 +612,7 @@ static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
   
   Wsq = W*W;
   
-  p_tmp = eos_info(W, vsq, &dPdW, &dPdvsq);
+  p_tmp = eos_info(W, vsq, &dPdW, &dPdvsq, gammaeos);
 
   // These expressions were calculated using Mathematica, but made into efficient 
   // code using Maple.  Since we know the analytic form of the equations, we can 
@@ -667,13 +666,15 @@ static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
  
       -- returns with all the EOS-related values needed;
  **********************************************************************/
-static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq)
+static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos)
 {
   register double ftmp,gtmp;
 
   ftmp = 1. - vsq;
   gtmp = sqrt(ftmp);
 
+  double gam_m1_o_gam = ((gammaeos-1.)/gammaeos);
+
   *dpdw =  gam_m1_o_gam * ftmp ;
   *dpdvsq =  gam_m1_o_gam * ( 0.5 * D/gtmp  -  W ) ;
 
diff --git a/src/GRHydro_InterfacesM.h b/src/GRHydro_InterfacesM.h
index 7479278..c6b49af 100644
--- a/src/GRHydro_InterfacesM.h
+++ b/src/GRHydro_InterfacesM.h
@@ -4,7 +4,7 @@ module Con2PrimM_fortran_interfaces
   interface
 
      subroutine GRHydro_Con2PrimM_pt(  handle, &
-                        dens, &
+				       local_gam, dens, &
                         sx,  sy,  sz, &
                         tau, &
                         rho, &
@@ -23,6 +23,7 @@ module Con2PrimM_fortran_interfaces
 
        implicit none       
        CCTK_INT handle
+       CCTK_REAL local_gam
        CCTK_REAL dens 
        CCTK_REAL sx, sy, sz
        CCTK_REAL tau 
@@ -41,10 +42,10 @@ module Con2PrimM_fortran_interfaces
        CCTK_REAL retval
      end subroutine  GRHydro_Con2PrimM_pt
 
-     subroutine GRHydro_Con2PrimM_Polytype_pt(  handle, &
+  subroutine GRHydro_Con2PrimM_Polytype_pt(  handle, local_gam,&
                         dens, &
                         sx,  sy,  sz, &
-                        tau, &
+                        sc, &
                         rho, &
                         velx,  vely,  velz,&
                         epsilon,  pressure,&
@@ -61,9 +62,10 @@ module Con2PrimM_fortran_interfaces
 
        implicit none       
        CCTK_INT handle
+       CCTK_REAL local_gam
        CCTK_REAL dens 
        CCTK_REAL sx, sy, sz
-       CCTK_REAL tau 
+       CCTK_REAL sc 
        CCTK_REAL rho
        CCTK_REAL velx, vely, velz
        CCTK_REAL epsilon, pressure
diff --git a/src/GRHydro_RegisterVarsM.cc b/src/GRHydro_RegisterVarsM.cc
index e65fae2..efe6c6a 100644
--- a/src/GRHydro_RegisterVarsM.cc
+++ b/src/GRHydro_RegisterVarsM.cc
@@ -94,16 +94,22 @@ extern "C"void GRHydro_RegisterM(CCTK_ARGUMENTS)
 
     // tracer
     if (evolve_tracer != 0)
-      register_evolved("GRHydro::GRHydro_cons_tracers", "GRHydro::GRHydro_tracer_rhs");
-
-    if (evolve_Y_e!=0) 
+      register_evolved("GRHydro::GRHydro_cons_tracers", 
+		       "GRHydro::GRHydro_tracer_rhs");
+    
+    // note that we set in pararamcheck
+    // evolve_y_e and evolve_temper, but MoLRegister
+    // happens before paramcheck...
+    if (CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) {
       register_constrained("HydroBase::Y_e");
       register_evolved("GRHydro::Y_e_con", 
 		       "GRHydro::Y_e_con_rhs");
+    }
 
-    if (evolve_temper!=0)
-     register_constrained("HydroBase::temperature");
-
+    if (CCTK_EQUALS(temperature_evolution_method,"GRHydro")) {
+      register_constrained("HydroBase::temperature");
+      register_constrained("HydroBase::entropy");
+    }
     // particles
     if (number_of_particles > 0)
       register_evolved("GRHydro::particles", "GRHydro::particle_rhs");
diff --git a/src/make.code.defn b/src/make.code.defn
index b28c974..30e5805 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -51,6 +51,7 @@ SRCS = 	Utils.F90 \
 	GRHydro_CalcUpdateM.F90 \
 	GRHydro_Con2PrimM.F90 \
 	GRHydro_Con2PrimM_pt.c \
+	GRHydro_Con2PrimM_polytype_pt.c \
 	GRHydro_EigenproblemM.F90 \
 	GRHydro_FluxM.F90 \
 	GRHydro_HLLE_EOSGM.F90 \