merge branch hot_and_MHD_temp_dev into branch at revision r185

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

1 files changed, 366 insertions, 152 deletions

diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90
index 24811c5..3eea7e2 100644
--- a/src/GRHydro_Con2Prim.F90
+++ b/src/GRHydro_Con2Prim.F90
@@ -49,13 +49,6 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
   
-#if !USE_EOS_OMNI
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-#endif
-
   integer :: i, j, k, itracer, nx, ny, nz
   CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, psi4pt, pmin, epsmin
   logical :: epsnegative
@@ -66,7 +59,6 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
 
   CCTK_REAL :: local_min_tracer
 
-#if USE_EOS_OMNI
 ! begin EOS Omni vars
   integer :: n = 1
   integer :: keytemp = 0
@@ -77,7 +69,6 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
   real*8  :: xtemp = 0.0d0
   real*8  :: xye = 0.0d0
 ! end EOS Omni vars
-#endif
 
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
@@ -90,19 +81,14 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
   if (use_min_tracer .ne. 0) then
     local_min_tracer = min_tracer
   else
-    local_min_tracer = 0d0
+    local_min_tracer = 0.0d0
   end if
 
-#if USE_EOS_OMNI
   ! this is a poly call
   call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
          GRHydro_rho_min,xeps,xtemp,xye,pmin,keyerr,anyerr)
   call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
-#else           
-  pmin   = EOS_Pressure(GRHydro_polytrope_handle, GRHydro_rho_min, 1.0d0)
-  epsmin = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, GRHydro_rho_min, pmin)
-#endif
 
   !$OMP PARALLEL DO PRIVATE(i,j,itracer,&
   !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, psi4pt, epsnegative)
@@ -139,7 +125,12 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
 
            enddo
            
-         endif
+        endif
+        
+        if(evolve_Y_e.ne.0) then
+           Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
+        endif
+
 
          if ( dens(i,j,k) .le. sqrt(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then
 
@@ -148,16 +139,13 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
            scon(i,j,k,:) = 0.d0
            vel(i,j,k,:) = 0.d0
            w_lorentz(i,j,k) = 1.d0
-#if USE_EOS_OMNI
+
            call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                 rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
 
            call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                 rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
-#else
-           press(i,j,k) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,j,k), eps(i,j,k))
-           eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, rho(i,j,k), press(i,j,k))
-#endif
+
            ! w_lorentz=1, so the expression for tau reduces to:
            tau(i,j,k)  = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k)
 
@@ -165,13 +153,21 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
 
          end if
 
-
-        call Con2Prim_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), &
-             uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
-             z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & 
-             GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
+         if(evolve_temper.eq.0) then
+            call Con2Prim_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), &
+                 uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
+                 z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & 
+                 GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
+         else
+            call Con2Prim_pt_hot(i,j,k,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),temperature(i,j,k),y_e(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
+                 uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
+                 z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & 
+                 GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
+         endif
 
         if (epsnegative) then  
 
@@ -202,16 +198,13 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
           scon(i,j,k,:) = 0.d0
           vel(i,j,k,:) = 0.d0
           w_lorentz(i,j,k) = 1.d0
-#if USE_EOS_OMNI
+
            call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                 rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
 
            call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                 rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
-#else
-           press(i,j,k) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,j,k), eps(i,j,k))
-           eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, rho(i,j,k), press(i,j,k))
-#endif
+
           ! w_lorentz=1, so the expression for tau reduces to:
           tau(i,j,k)  = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k)
 
@@ -270,21 +263,14 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
        y, z, r, GRHydro_rho_min
   
   CCTK_REAL s2, c0, c1, c2, c3, c4, f, df, ftol, v2, w, vlowx, vlowy, vlowz
-  CCTK_INT count, i, handle, GRHydro_reflevel, GRHydro_C2P_failed
+  CCTK_INT count, i, handle, GRHydro_reflevel
+  CCTK_REAL GRHydro_C2P_failed
   CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, epsold, epsnew, w2, &
        w2mhalf, temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin
   character(len=200) warnline
   logical epsnegative
 
-#if !USE_EOS_OMNI
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-#endif
-
 
-#if USE_EOS_OMNI
 ! begin EOS Omni vars
   integer :: n = 1
   integer :: keytemp = 0
@@ -294,7 +280,6 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
   real*8  :: xtemp = 0.0d0
   real*8  :: xye = 0.0d0
 ! end EOS Omni vars
-#endif
 
   
 !!$  Undensitize the variables 
@@ -308,13 +293,9 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
        2.*usx*usz*uxz + 2.*usy*usz*uyz
 
 !!$  Set initial guess for pressure:
-#if USE_EOS_OMNI
   call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
                       rho,epsilon,xtemp,xye,xpress,keyerr,anyerr)
   pold = max(1.d-10,xpress)
-#else
-  pold = max(1.d-10,EOS_Pressure(handle, rho, epsilon))
-#endif
 
 
 !!$  Check that the variables have a chance of being physical
@@ -343,15 +324,10 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
   
 !!$  Calculate the function
 
-#if USE_EOS_OMNI
   call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
                       rho,epsilon,xtemp,xye,xpress,keyerr,anyerr)
 
   f = pold - xpress
-#else
-  f = pold - EOS_Pressure(handle, rho, epsilon)
-#endif
-
 
 !!$Find the root
   
@@ -376,6 +352,7 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
       call CCTK_WARN(1,"Setting the point to atmosphere")
       !$OMP END CRITICAL
 
+
       ! for safety, let's set the point to atmosphere
       rho = GRHydro_rho_min
       udens = rho
@@ -395,17 +372,11 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
       goto 51
     end if
 
-#if USE_EOS_OMNI
     call EOS_Omni_DPressByDRho(handle,keytemp,GRHydro_eos_rf_prec,n,&
          rho,epsilon,xtemp,xye,dpressbydrho,keyerr,anyerr)
 
     call EOS_Omni_DPressByDEps(handle,keytemp,GRHydro_eos_rf_prec,n,&
          rho,epsilon,xtemp,xye,dpressbydeps,keyerr,anyerr)
-#else
-    dpressbydrho = EOS_DPressByDRho(handle,rho,epsilon)
-    dpressbydeps = EOS_DPressByDEps(handle,rho,epsilon)
-#endif
-
 
     temp1 = (utau+udens+pnew)**2 - s2
     drhobydpress = udens * s2 / (sqrt(temp1)*(udens+utau+pnew)**2)
@@ -424,14 +395,10 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
     epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz - &
          udens)/udens
 
-#if USE_EOS_OMNI
     call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
          rho,epsilon,xtemp,xye,xpress,keyerr,anyerr)
     
     f = pnew - xpress
-#else
-    f = pnew - EOS_Pressure(handle, rho, epsilon)
-#endif
 
   enddo
   
@@ -439,16 +406,11 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
 
   do i=1,GRHydro_polish
 
-#if USE_EOS_OMNI
     call EOS_Omni_DPressByDRho(handle,keytemp,GRHydro_eos_rf_prec,n,&
          rho,epsilon,xtemp,xye,dpressbydrho,keyerr,anyerr)
 
     call EOS_Omni_DPressByDEps(handle,keytemp,GRHydro_eos_rf_prec,n,&
          rho,epsilon,xtemp,xye,dpressbydeps,keyerr,anyerr)
-#else
-    dpressbydrho = EOS_DPressByDRho(handle,rho,epsilon)
-    dpressbydeps = EOS_DPressByDEps(handle,rho,epsilon)
-#endif
 
     temp1 = (utau+udens+pnew)**2 - s2
     drhobydpress = udens * s2 / (sqrt(temp1)*(udens+utau+pnew)**2)
@@ -466,14 +428,10 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
     epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz - &
          udens)/udens
 
-#if USE_EOS_OMNI
     call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n, &
          rho,epsilon,xtemp,xye,xpress,keyerr,anyerr)
     
     f = pold - xpress
-#else
-    f = pold - EOS_Pressure(handle, rho, epsilon)
-#endif
 
 
   enddo
@@ -514,6 +472,305 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
 
 end subroutine Con2Prim_pt
 
+subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, vely, &
+     velz, epsilon, temp, ye, press, w_lorentz, uxx, uxy, uxz, uyy, &
+     uyz, uzz, det, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, &
+     GRHydro_reflevel, GRHydro_C2P_failed)
+  
+  implicit none
+  
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_REAL dens, sx, sy, sz, tau, rho, velx, vely, velz, epsilon, &
+       press, uxx, uxy, uxz, uyy, uyz, uzz, det, w_lorentz, x, &
+       y, z, r, GRHydro_rho_min
+  CCTK_REAL temp, ye
+  CCTK_REAL s2, c0, c1, c2, c3, c4, f, df, ftol, v2, w, vlowx, vlowy, vlowz
+  CCTK_INT ii,jj,kk,count, i, handle, GRHydro_reflevel
+  CCTK_REAL GRHydro_C2P_failed
+  CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, epsold, epsnew, w2, &
+       w2mhalf, temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin
+  character(len=200) warnline
+  logical epsnegative
+
+  integer :: failwarnmode = 0
+  integer :: failinfomode = 1
+
+! begin EOS Omni vars
+  integer :: n = 1
+  integer :: keytemp = 0
+  integer :: anyerr = 0
+  integer :: keyerr(1) = 0
+  real*8  :: xpress = 0.0d0
+  real*8  :: temp0 = 0.0d0
+! end EOS Omni vars
+  
+! set pmin and epsmin to something sensible:
+  pmin = 1.0d-28
+  epsmin = 1.0e-5
+  
+  
+  if(con2prim_oct_hack.ne.0.and.&
+       x .lt. 0.0d0 .or.&
+       y .lt. 0.0d0 .or.&
+       z .lt. 0.0d0) then
+     failwarnmode = 2
+     failinfomode = 2
+  endif
+
+!!$  Undensitize the variables 
+
+  udens = dens /sqrt(det)
+  usx = sx /sqrt(det)
+  usy = sy /sqrt(det)
+  usz = sz /sqrt(det)
+  utau = tau /sqrt(det)
+  s2 = usx*usx*uxx + usy*usy*uyy + usz*usz*uzz + 2.*usx*usy*uxy + &
+       2.*usx*usz*uxz + 2.*usy*usz*uyz
+
+!!$  Set initial guess for pressure:
+  temp0 = temp
+  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
+                      rho,epsilon,temp,ye,xpress,keyerr,anyerr)
+  pold = max(1.d-15,xpress)
+  ! error handling
+  if(anyerr.ne.0) then
+     call CCTK_WARN(failinfomode,"EOS error in c2p 0")
+     write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
+     call CCTK_WARN(failinfomode,warnline)
+     write(warnline,"(1P10E15.6)") rho,dens,epsilon,temp,temp0,ye
+     call CCTK_WARN(failinfomode,warnline)
+     write(warnline,"(A7,i8)") "code: ",keyerr(1)
+     call CCTK_WARN(failinfomode,warnline)
+     write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
+     call CCTK_WARN(failinfomode,warnline)
+     call CCTK_WARN(failwarnmode,"Aborting!!!")
+  endif
+
+
+!!$  Check that the variables have a chance of being physical
+
+  if( (utau + pold + udens)**2 - s2 .le. 0.0d0) then
+
+    if (c2p_reset_pressure .ne. 0) then
+      pold = sqrt(s2 + c2p_reset_pressure_to_value) - utau - udens
+    else 
+      !$OMP CRITICAL
+!!!!      call CCTK_WARN(GRHydro_NaN_verbose, "c2p failed and being told not to reset the pressure")
+      GRHydro_C2P_failed = 1
+      !$OMP END CRITICAL
+    endif
+    
+  endif
+  
+!!$  Calculate rho and epsilon 
+
+!define temporary variables to speed up
+
+  rho = udens * sqrt( (utau + pold + udens)**2 - s2)/(utau + pold + udens)
+  w_lorentz = (utau + pold + udens) / sqrt( (utau + pold + udens)**2 - s2)
+  epsilon = (sqrt( (utau + pold + udens)**2 - s2) - pold * w_lorentz - &
+       udens)/udens
+  
+!!$  Calculate the function
+
+  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
+                      rho,epsilon,temp,ye,xpress,keyerr,anyerr)
+  ! error handling
+  if(anyerr.ne.0) then
+     call CCTK_WARN(failinfomode,"EOS error in c2p 1")
+     write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
+     call CCTK_WARN(failinfomode,warnline)
+     write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye
+     call CCTK_WARN(failinfomode,warnline)
+     write(warnline,"(A7,i8)") "code: ",keyerr(1)
+     call CCTK_WARN(failinfomode,warnline)
+     write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
+     call CCTK_WARN(failinfomode,warnline)
+     call CCTK_WARN(failwarnmode,"Aborting!!!")
+  endif
+
+  f = pold - xpress
+
+
+!!$Find the root
+  
+  count = 0
+  pnew = pold
+  do while ( ((abs(pnew - pold)/abs(pnew) .gt. GRHydro_perc_ptol) .and. &
+       (abs(pnew - pold) .gt. GRHydro_del_ptol))  .or. &
+       (count .lt. GRHydro_countmin))
+    count = count + 1
+
+    if (count > GRHydro_countmax) then
+      GRHydro_C2P_failed = 1
+
+      !$OMP CRITICAL
+      call CCTK_WARN(failinfomode, 'count > GRHydro_countmax! ')
+      write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
+      call CCTK_WARN(failinfomode,warnline)
+      write(warnline,'(a20,3g16.7)') 'xyz location: ',x,y,z
+      call CCTK_WARN(failinfomode,warnline)
+      write(warnline,'(a20,3i5)') 'ijk location: ',ii,jj,kk
+      call CCTK_WARN(failinfomode,warnline)
+      write(warnline,'(a20,g16.7)') 'radius: ',r
+      call CCTK_WARN(failinfomode,warnline)
+      call CCTK_WARN(failinfomode,"Setting the point to atmosphere")
+      !$OMP END CRITICAL
+
+      ! for safety, let's set the point to atmosphere
+      rho = GRHydro_rho_min
+      udens = rho
+      dens = sqrt(det) * rho
+      pnew = pmin
+      epsilon = epsmin
+      ! w_lorentz=1, so the expression for utau reduces to:
+      utau  = rho + rho*epsmin - udens
+      sx = 0.d0
+      sy = 0.d0
+      sz = 0.d0
+      s2 = 0.d0
+      usx = 0.d0
+      usy = 0.d0
+      usz = 0.d0
+      w_lorentz = 1.d0
+      goto 51
+    end if
+
+    call EOS_Omni_DPressByDRho(handle,keytemp,GRHydro_eos_rf_prec,n,&
+         rho,epsilon,temp,ye,dpressbydrho,keyerr,anyerr)
+
+    call EOS_Omni_DPressByDEps(handle,keytemp,GRHydro_eos_rf_prec,n,&
+         rho,epsilon,temp,ye,dpressbydeps,keyerr,anyerr)
+
+    temp1 = (utau+udens+pnew)**2 - s2
+    drhobydpress = udens * s2 / (sqrt(temp1)*(udens+utau+pnew)**2)
+    depsbydpress = pnew * s2 / (rho * (udens + utau + pnew) * temp1)
+    df = 1.0d0 - dpressbydrho*drhobydpress - &
+         dpressbydeps*depsbydpress
+
+    pold = pnew
+    pnew = max(pold - f/df, pmin)
+    
+!!$    Recalculate primitive variables and function
+       
+    rho = udens * sqrt( (utau + pnew + udens)**2 - s2)/(utau + pnew + udens)
+    w_lorentz = (utau + pnew + udens) / sqrt( (utau + pnew + udens)**2 - &
+         s2)
+    epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz - &
+         udens)/udens
+
+!    epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz)/udens - &
+!         1.0d0
+
+    call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
+         rho,epsilon,temp,ye,xpress,keyerr,anyerr)
+    ! error handling
+    if(anyerr.ne.0) then
+       call CCTK_WARN(failinfomode,"EOS error in c2p 2")
+       write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
+       call CCTK_WARN(failinfomode,warnline)
+       write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye
+       call CCTK_WARN(failinfomode,warnline)
+       write(warnline,"(A7,i8)") "code: ",keyerr(1)
+       call CCTK_WARN(failinfomode,warnline)
+       write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
+       call CCTK_WARN(failinfomode,warnline)
+       call CCTK_WARN(failwarnmode,"Aborting!!!")
+    endif
+    
+    f = pnew - xpress
+
+  enddo
+
+  
+!!$  Polish the root
+
+  do i=1,GRHydro_polish
+
+    call EOS_Omni_DPressByDRho(handle,keytemp,GRHydro_eos_rf_prec,n,&
+         rho,epsilon,temp,ye,dpressbydrho,keyerr,anyerr)
+
+    call EOS_Omni_DPressByDEps(handle,keytemp,GRHydro_eos_rf_prec,n,&
+         rho,epsilon,temp,ye,dpressbydeps,keyerr,anyerr)
+
+    temp1 = (utau+udens+pnew)**2 - s2
+    drhobydpress = udens * s2 / (sqrt(temp1)*(udens+utau+pnew)**2)
+    depsbydpress = pnew * s2 / (rho * (udens + utau + pnew) * temp1)
+    df = 1.0d0 - dpressbydrho*drhobydpress - &
+         dpressbydeps*depsbydpress
+    pold = pnew
+    pnew = pold - f/df
+    
+!!$    Recalculate primitive variables and function
+
+    rho = udens * sqrt( (utau + pnew + udens)**2 - s2)/(utau + pnew + udens)
+    w_lorentz = (utau + pnew + udens) / sqrt( (utau + pnew + udens)**2 - &
+         s2)
+    epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz - &
+         udens)/udens
+
+    call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n, &
+         rho,epsilon,temp,ye,xpress,keyerr,anyerr)
+
+
+    ! error handling
+    if(anyerr.ne.0) then
+       call CCTK_WARN(failinfomode,"EOS error in c2p 3")
+       write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
+       call CCTK_WARN(failinfomode,warnline)
+       write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye
+       call CCTK_WARN(failinfomode,warnline)
+       write(warnline,"(A7,i8)") "code: ",keyerr(1)
+       call CCTK_WARN(failinfomode,warnline)
+       write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
+       call CCTK_WARN(failinfomode,warnline)
+       call CCTK_WARN(failwarnmode,"Aborting!!!")
+    endif
+
+    
+    f = pold - xpress
+
+  enddo
+
+!!$  Calculate primitive variables from root
+
+  if (rho .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then
+    rho = GRHydro_rho_min
+    udens = rho
+    dens = sqrt(det) * rho
+!    epsilon = (sqrt( (utau + pnew + udens)**2) - pnew -  udens)/udens
+    epsilon = epsmin
+    ! w_lorentz=1, so the expression for utau reduces to:
+    utau  = rho + rho*epsmin - udens
+    sx = 0.d0
+    sy = 0.d0
+    sz = 0.d0
+    s2 = 0.d0
+    usx = 0.d0
+    usy = 0.d0
+    usz = 0.d0
+    w_lorentz = 1.d0
+  end if
+
+51 press = pnew
+  vlowx = usx / ( (rho + rho*epsilon + press) * w_lorentz**2)
+  vlowy = usy / ( (rho + rho*epsilon + press) * w_lorentz**2)
+  vlowz = usz / ( (rho + rho*epsilon + press) * w_lorentz**2)
+  velx = uxx * vlowx + uxy * vlowy + uxz * vlowz
+  vely = uxy * vlowx + uyy * vlowy + uyz * vlowz
+  velz = uxz * vlowx + uyz * vlowy + uzz * vlowz
+  
+!!$If all else fails, use the polytropic EoS
+
+!  eps may well be negative; we don't mess with this
+!  if(epsilon .lt. 0.0d0) then
+!    epsnegative = .true.
+!  endif
+
+
+end subroutine Con2Prim_pt_hot
+
  /*@@
    @routine    Conservative2PrimitiveBoundaries
    @date       Tue Mar 12 18:04:40 2002
@@ -557,15 +814,6 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS)
 
   CCTK_REAL :: local_min_tracer
 
-#if !USE_EOS_OMNI
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-#endif
-
-
-#if USE_EOS_OMNI
 ! begin EOS omni
   CCTK_INT  :: keyerr(1) = 0
   CCTK_INT  :: anyerr = 0
@@ -575,20 +823,13 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS)
   CCTK_REAL :: xeps = 0.0d0
   CCTK_REAL :: xtemp = 0.0d0
 ! end EOS omni
-#endif
   
-
-#if USE_EOS_OMNI
   ! this is a poly call
   call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
          GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
 
   call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,epsmin,xtemp,xye,pmin,epsmin,keyerr,anyerr)
-#else
-  pmin=EOS_Pressure(GRHydro_polytrope_handle, GRHydro_rho_min, 1.0d0) 
-  epsmin = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, GRHydro_rho_min, pmin)
-#endif
 
                                                                               
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
@@ -649,6 +890,10 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS)
 
         endif
 
+        if(evolve_Y_e.ne.0) then
+           Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
+        endif
+
         call Con2Prim_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),&
@@ -832,6 +1077,10 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS)
 
         endif
 
+        if(evolve_Y_e.ne.0) then
+           Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
+        endif
+
         call Con2Prim_ptPolytype(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), &
@@ -884,15 +1133,6 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
        enthalpy, denthalpy, sqrtdet, invsqrtdet, invfac, GRHydro_C2P_failed
   character(len=200) warnline
 
-#if !USE_EOS_OMNI
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-#endif
-
-
-#if USE_EOS_OMNI
 ! begin EOS Omni vars
   integer :: n = 1
   integer :: keytemp = 0
@@ -903,7 +1143,7 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
   real*8  :: xtemp = 0.0d0
   real*8  :: xye = 0.0d0
 ! end EOS Omni vars
-#endif
+
 
 !!$  Undensitize the variables 
 
@@ -919,7 +1159,6 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
 !!$  Set initial guess for rho:
   rhoold = max(GRHydro_rho_min,rho)
 
-#if USE_EOS_OMNI
 !  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
 !       rhoold,xeps,xtemp,xye,xpress,keyerr,anyerr)
 
@@ -927,19 +1166,12 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
        rhoold,xeps,xtemp,xye,press,xeps,keyerr,anyerr)
 
   enthalpy = 1.0d0 + xeps + press / rhoold
-#else
-  enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhoold, press) + &
-       EOS_Pressure(handle, rhoold, epsilon) / rhoold
-#endif
 
   w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 ))
 
-#if USE_EOS_OMNI
   call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
        rhoold,epsilon,xtemp,xye,press,keyerr,anyerr)
-#else
-  press = EOS_Pressure(handle, rhoold, epsilon)
-#endif
+
 !!$  Calculate the function
 
   f = rhoold * w_lorentz - udens
@@ -972,16 +1204,13 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
       rhonew = GRHydro_rho_min
       udens = rhonew
       dens = sqrtdet * rhonew
-#if USE_EOS_OMNI
+
       call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
            rhonew,1.0d0,xtemp,xye,press,keyerr,anyerr)
 
       call EOS_Omni_EpsFromPress(handle,keytemp,GRHydro_eos_rf_prec,n,&
            rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
-#else
-      press = EOS_Pressure(handle, rhonew, 1.d0)
-      epsilon = EOS_SpecificIntEnergy(handle, rhonew, press)
-#endif
+
       sx = 0.d0
       sy = 0.d0
       sz = 0.d0
@@ -997,13 +1226,9 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
 !!$    Ian has the feeling that this is an error in general. But for the
 !!$    2D_Polytrope it gives the right answers.
 
-#if USE_EOS_OMNI
       call EOS_Omni_DPressByDrho(handle,keytemp,GRHydro_eos_rf_prec,n,&
            rhonew,epsilon,xtemp,xye,denthalpy,keyerr,anyerr)
       denthalpy = denthalpy / rhonew
-#else
-    denthalpy = EOS_DPressByDrho(handle, rhonew, epsilon) / rhonew
-#endif
 
     df = w_lorentz - rhonew * s2 * denthalpy / &
          (w_lorentz*(udens**2)*(enthalpy**3))
@@ -1024,19 +1249,13 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
     
 !!$    Recalculate primitive variables and function
 
-#if USE_EOS_OMNI
-  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
-       rhonew,epsilon,xtemp,xye,press,keyerr,anyerr)
-
-  call EOS_Omni_EpsFromPress(handle,keytemp,GRHydro_eos_rf_prec,n,&
-       rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
+    call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
+         rhonew,epsilon,xtemp,xye,press,keyerr,anyerr)
+    
+    call EOS_Omni_EpsFromPress(handle,keytemp,GRHydro_eos_rf_prec,n,&
+         rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
 
-  enthalpy = 1.0d0 + epsilon + press / rhonew
-#else
-  enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhonew, press) + &
-         EOS_Pressure(handle, rhonew, epsilon) / rhonew
-  press = EOS_Pressure(handle, rhonew, epsilon)
-#endif
+    enthalpy = 1.0d0 + epsilon + press / rhonew
 
     w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 ))
     tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens
@@ -1055,16 +1274,12 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
     ! significant changes in the results
     dens = sqrtdet * rhonew
 
-#if USE_EOS_OMNI
   call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
        rhonew,1.0d0,xtemp,xye,press,keyerr,anyerr)
 
   call EOS_Omni_EpsFromPress(handle,keytemp,GRHydro_eos_rf_prec,n,&
        rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
-#else
-    press = EOS_Pressure(handle, rhonew, 1.d0)
-    epsilon = EOS_SpecificIntEnergy(handle, rhonew, press)
-#endif
+
     ! w_lorentz=1, so the expression for utau reduces to:
     tau  = rhonew + rhonew*epsilon - udens
     sx = 0.d0
@@ -1079,7 +1294,6 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
 
 50  rho = rhonew
 
-#if USE_EOS_OMNI
   call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
        rhonew,epsilon,xtemp,xye,press,keyerr,anyerr)
 
@@ -1087,22 +1301,13 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
        rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
 
   enthalpy = 1.0d0 + epsilon + press / rhonew
-#else
-  enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhonew, press) + &
-       EOS_Pressure(handle, rhonew, epsilon) / rhonew
-  w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 ))
-#endif
 
-#if USE_EOS_OMNI
   call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
        rhonew,epsilon,xtemp,xye,press,keyerr,anyerr)
 
   call EOS_Omni_EpsFromPress(handle,keytemp,GRHydro_eos_rf_prec,n,&
        rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
-#else
-  press = EOS_Pressure(handle, rhonew, epsilon)
-  epsilon = EOS_SpecificIntEnergy(handle, rhonew, press)
-#endif
+
   tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens
 
   if (epsilon .lt. 0.0d0) then
@@ -1121,15 +1326,12 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
     
     rho = GRHydro_rho_min
     dens = sqrtdet * rho
-#if USE_EOS_OMNI
-  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
-       rhonew,1.d0,xtemp,xye,press,keyerr,anyerr)
-  call EOS_Omni_EpsFromPress(handle,keytemp,GRHydro_eos_rf_prec,n,&
-       rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
-#else
-    press = EOS_Pressure(handle, rhonew, 1.d0)
-    epsilon = EOS_SpecificIntEnergy(handle, rhonew, press)
-#endif
+
+    call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
+         rhonew,1.d0,xtemp,xye,press,keyerr,anyerr)
+    call EOS_Omni_EpsFromPress(handle,keytemp,GRHydro_eos_rf_prec,n,&
+         rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr)
+
     ! w_lorentz=1, so the expression for tau reduces to:
     tau  = sqrtdet * (rho+rho*epsilon) - dens
     usx = 0.d0
@@ -1242,6 +1444,10 @@ subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS)
              x(i,j,k)+0.5d0*CCTK_DELTA_SPACE(1),&
              y(i,j,k)+0.5d0*CCTK_DELTA_SPACE(2), &
              z(i,j,k)+0.5d0*CCTK_DELTA_SPACE(3),r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
+
+        if(evolve_Y_e.ne.0) then
+           Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
+        endif
         
       end do
     end do
@@ -1506,6 +1712,10 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
           enddo
         endif
 
+        if(evolve_Y_e.ne.0) then
+           Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
+        endif
+
         if ( (dens(i,j,k).le.sqrt(GRHydro_Det(i,j,k)) * &
              GRHydro_rho_min*(1.0d0+GRHydro_atmo_tolerance)) &
              .or.(tau(i,j,k) .le. 0d0)) then
@@ -1992,6 +2202,10 @@ subroutine Con2PrimPolytypeGeneral(CCTK_ARGUMENTS)
           enddo
         end if
 
+        if(evolve_Y_e.ne.0) then
+           Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
+        endif
+
         udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
         usx   = scon(i,j,k,1)   / sqrt(GRHydro_det(i,j,k))
         usy   = scon(i,j,k,2)   / sqrt(GRHydro_det(i,j,k))