GRHydro: Add basic vector potential support

   Basic cell-centered, algebraic gauge vector potential
   method with place-holders for lorenz gauge.
   Initial Avec constrained to poloidal at the moment.

From: Tanja Bode <tanja.bode@physics.gatech.edu>

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

1 files changed, 1365 insertions, 0 deletions

diff --git a/src/GRHydro_Con2PrimAM.F90 b/src/GRHydro_Con2PrimAM.F90
new file mode 100644
index 0000000..50d701c
--- /dev/null
+++ b/src/GRHydro_Con2PrimAM.F90
@@ -0,0 +1,1365 @@
+/*@@
+   @file      GRHydro_Con2PrimAM.F90
+   @date      Sep 3, 2010
+   @author    Scott Noble, Joshua Faber, Bruno Mundim, Tanja Bode
+   @desc 
+   The routines for converting conservative to primitive variables.
+   Vector-potential MHD.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "SpaceMask.h"
+#include "GRHydro_InterfacesAM.h"
+#include "GRHydro_Macros.h"
+
+#define ITER_TOL (1.0e-8)
+#define MAXITER (50)
+
+ /*@@
+   @routine    Conservative2PrimitiveAM
+   @date       Sep 3, 2010
+   @author     Scott Noble, Joshua Faber, Bruno Mundim, Ian Hawke
+   @desc 
+   Wrapper routine that converts from conserved to primitive variables
+   at every grid cell centre.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   Trimmed and altered from the GR3D routines, original author Mark Miller.
+   2007?: Bruno excluded the points in the atmosphere and excision region from the computation.
+   Aug. 2008: Luca added a check on whether a failure at a given point may be disregarded, 
+   because that point will then be restricted from a finer level. This should be completely 
+   safe only if *regridding happens at times when all levels are evolved.*
+   Feb. 2009: The above procedure proved to be wrong, so Luca implemented another one. 
+   When a failure occurs, it is temporarily ignored, except for storing the location of where 
+   it occured in a mask. At the end, after a Carpet restriction, the mask is checked and if 
+   it still contains failures, the run is aborted with an error message. Only used routines 
+   have been updated to use this procedure.
+   @endhistory 
+
+@@*/
+
+subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS)
+
+  use Con2PrimM_fortran_interfaces
+
+  implicit none
+  
+  ! save memory when MP is not used
+  ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1
+  TARGET gaa, gab, gac, gbb, gbc, gcc
+  TARGET gxx, gxy, gxz, gyy, gyz, gzz
+  TARGET lvel, vel
+  TARGET lBvec, Bvec
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  integer :: i, j, k, itracer, nx, ny, nz
+  CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, sdet, pmin(1), epsmin(1)
+  CCTK_REAL :: oob, b2, d2, s2, bscon, bxhat, byhat, bzhat, bhatscon
+  CCTK_REAL :: Wm, Wm0, Wm_resid, Wmold
+  CCTK_REAL :: s2m, s2m0, s2m_resid, s2mold, s2max, taum
+  CCTK_INT :: niter
+  CCTK_INT :: epsnegative
+  character(len=100) warnline
+  
+  CCTK_REAL :: local_min_tracer, local_gam(1), local_pgam,local_K, sc
+  CCTK_REAL :: local_perc_ptol
+
+! begin EOS Omni vars                                                                                       
+  integer :: n,keytemp,anyerr,keyerr(1)
+  real*8  :: xpress(1),xtemp(1),xye(1),xeps(1),xrho(1),one(1)=1.0d0
+! end EOS Omni vars                 
+
+  ! save memory when MP is not used
+  CCTK_INT :: GRHydro_UseGeneralCoordinates
+  CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33
+  CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim
+
+  logical :: posdef
+
+  CCTK_REAL :: g11c, g12c, g13c, g22c, g23c, g33c
+  CCTK_REAL :: tmp1
+
+  ! Save the primitive variables to temporary functions before calling the
+  ! con2prim pointwise routines:
+  CCTK_REAL :: rho_tmp, press_tmp, eps_tmp
+  CCTK_REAL :: velx_tmp, vely_tmp, velz_tmp, w_lorentz_tmp
+  CCTK_REAL :: Bvecx_tmp, Bvecy_tmp, Bvecz_tmp
+  CCTK_REAL :: Bconsx_tmp, Bconsy_tmp, Bconsz_tmp
+
+  CCTK_REAL :: bdotv, magpress
+
+  ! Assume 3-metric is positive definite. Check deep inside the horizon 
+  ! if this is actually satisfied and if it is not then cast the metric 
+  !as conformally flat only for con2prim inversion purposes.
+  posdef = .true. 
+ 
+  !! Get Bvec from Avec before rest, so we can reuse pointwise MHD recovery
+  call GRHydro_BvecfromAvec(CCTK_PASS_FTOF)
+
+  if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
+    g11 => gaa
+    g12 => gab
+    g13 => gac
+    g22 => gbb
+    g23 => gbc
+    g33 => gcc
+    vup => lvel
+    Bprim => lBvec
+  else
+    g11 => gxx
+    g12 => gxy
+    g13 => gxz
+    g22 => gyy
+    g23 => gyz
+    g33 => gzz
+    vup => vel
+    Bprim => Bvec
+  end if
+#define gxx faulty_gxx
+#define gxy faulty_gxy
+#define gxz faulty_gxz
+#define gyy faulty_gyy
+#define gyz faulty_gyz
+#define gzz faulty_gzz
+#define betax faulty_betax
+#define betay faulty_betay
+#define betaz faulty_betaz
+#define vel faulty_vel
+#define Bvec faulty_Bvec
+
+! begin EOS Omni vars                                                                                       
+  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                 
+
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+  
+  if (use_min_tracer .ne. 0) then
+     local_min_tracer = min_tracer
+  else
+     local_min_tracer = 0d0
+  end if
+
+!  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)
+  ! this is a poly call
+  xrho(1)=GRHydro_rho_min
+  call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       xrho,xeps,xtemp,xye,pmin,keyerr,anyerr)
+  call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       xrho,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
+
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+         one,one,xtemp,xye,local_gam,keyerr,anyerr)
+  local_gam = local_gam+1.d0
+
+  !call CCTK_WARN(1,"In Con2PrimAM")
+
+  !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,&
+  !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, &
+  !$OMP b2,xrho,xeps,xpress,xtemp,local_K,local_pgam,sc,keyerr,anyerr,keytemp, &
+  !$OMP local_perc_ptol,posdef,g11c,g12c,g13c,g22c,g23c,g33c,tmp1, &
+  !$OMP sdet,d2,s2,oob,bscon,bxhat,byhat,bzhat, &
+  !$OMP bhatscon,Wm,Wm0,Wm_resid,Wmold,s2m,s2m0,s2m_resid,s2mold,s2max, &
+  !$OMP taum,niter,rho_tmp,press_tmp,eps_tmp,velx_tmp,vely_tmp,velz_tmp, &
+  !$OMP Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, &
+  !$OMP w_lorentz_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,bdotv,magpress)
+  do k = 1, nz 
+     do j = 1, ny 
+        do i = 1, nx
+           
+           !do not compute if in atmosphere or in excised region
+           if ((atmosphere_mask(i,j,k) .ne. 0) .or. &
+                (hydro_excision_mask(i,j,k) .ne. 0)) cycle
+           
+           epsnegative = 0
+           
+           det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k))
+           sdet = sqrt(det)
+
+
+           call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,&
+                g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),&
+                g23(i,j,k),g33(i,j,k))        
+           
+!!$ Tracers don't need an MHD treatment!
+           if (evolve_tracer .ne. 0) then
+              do itracer=1,number_of_tracers
+                 call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), tracer(i,j,k,itracer), &
+                      dens(i,j,k))
+                 
+                 if (use_min_tracer .ne. 0) then
+                    if (tracer(i,j,k,itracer) .le. local_min_tracer) then
+                       tracer(i,j,k,itracer) = local_min_tracer
+                    end if
+                 end if
+                 
+              enddo
+              
+           endif
+           
+           if(evolve_Y_e.ne.0) then
+              Y_e(i,j,k) = max(min(Y_e_con(i,j,k) / dens(i,j,k),GRHydro_Y_e_max),&
+                 GRHydro_Y_e_min)
+           endif
+
+           
+           b2=g11(i,j,k)*Bprim(i,j,k,1)**2+g22(i,j,k)*Bprim(i,j,k,2)**2+g33(i,j,k)*Bprim(i,j,k,3)**2+ &
+           2.0*(g12(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,2)+g13(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,3)+ &
+           g23(i,j,k)*Bprim(i,j,k,2)*Bprim(i,j,k,3))
+
+       
+           if ( dens(i,j,k) .le. sdet*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then
+              
+              !call CCTK_WARN(1,"Con2Prim: Resetting to atmosphere")
+              !write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k)
+              !call CCTK_WARN(1,warnline)
+              !write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k),&
+              !      temperature(i,j,k),y_e(i,j,k)
+              !call CCTK_WARN(1,warnline)
+
+              
+              dens(i,j,k) = sdet*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
+              rho(i,j,k) = GRHydro_rho_min
+              scon(i,j,k,:) = 0.d0
+              vup(i,j,k,:) = 0.d0
+              w_lorentz(i,j,k) = 1.d0
+              
+              if(evolve_temper.ne.0) then
+              ! set hot atmosphere values
+                temperature(i,j,k) = grhydro_hot_atmo_temp
+                y_e(i,j,k) = grhydro_hot_atmo_Y_e
+                y_e_con(i,j,k) = y_e(i,j,k) * dens(i,j,k)
+                keytemp = 1
+                call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
+                   press(i,j,k),keyerr,anyerr)
+              else
+                keytemp = 0
+                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)
+              endif
+
+              !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)
+
+              ! w_lorentz=1, so the expression for tau reduces to:
+              
+              !!$ tau does need to take into account the existing B-field
+              !!$ with w_lorentz=1, we find tau = sqrtdet*(rho (1+eps+b^2/2)) - dens  [Press drops out]
+              tau(i,j,k)  = sdet * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k)
+              
+              if(tau(i,j,k).le.sdet*b2*0.5d0)then
+                tau(i,j,k) = GRHydro_tau_min + sdet*b2*0.5d0
+              endif
+
+              cycle
+              
+           end if
+
+
+           if(evolve_temper.eq.0) then
+            
+
+              if(sqrtdet_thr.gt.0d0 .and. sdet.ge.sqrtdet_thr) then
+                    d2 = dens(i,j,k)**2
+                    s2 = uxx*scon(i,j,k,1)**2 + uyy*scon(i,j,k,2)**2 &
+                                              + uzz*scon(i,j,k,3)**2 &
+                       + 2.0d0*uxy*scon(i,j,k,1)*scon(i,j,k,2) &
+                       + 2.0d0*uxz*scon(i,j,k,1)*scon(i,j,k,3) &
+                       + 2.0d0*uyz*scon(i,j,k,2)*scon(i,j,k,3) 
+                   oob = 1.0d0/sqrt(b2)
+                 bxhat = oob*Bprim(i,j,k,1)
+                 byhat = oob*Bprim(i,j,k,2)
+                 bzhat = oob*Bprim(i,j,k,3)
+                 bhatscon = bxhat*scon(i,j,k,1)+byhat*scon(i,j,k,2) &
+                                               +bzhat*scon(i,j,k,3)
+                 bscon = Bprim(i,j,k,1)*scon(i,j,k,1) &
+                       + Bprim(i,j,k,2)*scon(i,j,k,2) &
+                       + Bprim(i,j,k,3)*scon(i,j,k,3)
+                 ! Initial guesses for iterative procedure to find Wm:
+                 Wm0 = sdet*sqrt(bhatscon**2+d2)
+                 s2m0 = (Wm0**2*s2+bhatscon**2*(b2+2.0d0*Wm0)) &
+                      / (Wm0+b2)**2 
+                 Wm = sdet*sqrt(s2m0+d2)
+                 s2m = (Wm**2*s2+bscon**2*(b2+2.0d0*Wm)) &
+                     / (Wm+b2)**2 
+                 s2m_resid = 1.0d60
+                 Wm_resid = 1.0d60
+                 niter = 0
+                 do while((s2m_resid.ge.ITER_TOL.and.Wm_resid.ge.ITER_TOL).and.&
+                          niter.le.MAXITER)
+                   Wmold = Wm
+                   s2mold = s2m
+                   Wm = sdet*sqrt(s2m+d2)
+                   s2m = (Wm**2*s2+bscon**2*(b2+2.0d0*Wm)) &
+                       / (Wm+b2)**2 
+                   Wm_resid = abs(Wmold-Wm)
+                   s2m_resid = abs(s2mold-s2m)
+                   niter = niter + 1
+                 end do
+                 !TODO: abort execution if niter .eq. MAXITER and warn user
+                 taum = tau(i,j,k) - 0.5d0*sdet*b2 -0.5d0*(b2*s2-bscon**2)/ &
+                                                         (sdet*(Wm+b2)**2) 
+                 s2max = taum*(taum+2.0d0*dens(i,j,k))
+                 if(taum.lt.GRHydro_tau_min)then
+                   tau(i,j,k) = GRHydro_tau_min + 0.5d0*sdet*b2 + 0.5d0* &
+                                (b2*s2-bscon**2)/(sdet*(Wm+b2)**2)
+                 end if
+                 if(s2.gt.s2max) then
+                   scon(i,j,k,1) = scon(i,j,k,1)*sqrt(s2max/s2)
+                   scon(i,j,k,2) = scon(i,j,k,2)*sqrt(s2max/s2)
+                   scon(i,j,k,3) = scon(i,j,k,3)*sqrt(s2max/s2)
+                 end if
+              endif 
+
+              rho_tmp = rho(i,j,k)
+              press_tmp = press(i,j,k)
+              eps_tmp = eps(i,j,k)
+              velx_tmp = vup(i,j,k,1)
+              vely_tmp = vup(i,j,k,2)
+              velz_tmp = vup(i,j,k,3)
+              w_lorentz_tmp = w_lorentz(i,j,k)
+              Bvecx_tmp = Bprim(i,j,k,1)
+              Bvecy_tmp = Bprim(i,j,k,2)
+              Bvecz_tmp = Bprim(i,j,k,3)
+
+              !! We've already create B from A, construct Bcons
+              Bconsx_tmp = sdet*Bvecx_tmp 
+              Bconsy_tmp = sdet*Bvecy_tmp 
+              Bconsz_tmp = sdet*Bvecz_tmp 
+
+              keytemp = 0
+
+              call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, &
+                   GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), &
+                   scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), &
+                   Bconsx_tmp, Bconsy_tmp, Bconsz_tmp,xye(1), &
+                   xtemp(1),rho_tmp,velx_tmp,vely_tmp,velz_tmp,&
+                   eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,&
+                   w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),&
+                   g22(i,j,k),g23(i,j,k),g33(i,j,k), &
+                   uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                   epsnegative,GRHydro_C2P_failed(i,j,k))
+
+              if(evolve_entropy.ne.0) then
+                 entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho(i,j,k)
+              endif
+
+              if(GRHydro_C2P_failed(i,j,k).ne.0) then
+                xrho=1.0d0; xtemp=0.0d0; xeps=1.0d0
+                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=10.0d0; xeps=1.0d0
+                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)
+
+                if(sdet.ge.sqrtdet_thr) then
+                  GRHydro_C2P_failed(i,j,k) = 0
+
+                  rho_tmp = rho(i,j,k)
+                  press_tmp = press(i,j,k)
+                  eps_tmp = eps(i,j,k)
+                  velx_tmp = vup(i,j,k,1)
+                  vely_tmp = vup(i,j,k,2)
+                  velz_tmp = vup(i,j,k,3)
+                  w_lorentz_tmp = w_lorentz(i,j,k)
+                  Bvecx_tmp = Bprim(i,j,k,1)
+                  Bvecy_tmp = Bprim(i,j,k,2)
+                  Bvecz_tmp = Bprim(i,j,k,3)
+                  Bconsx_tmp = sdet*Bvecx_tmp 
+                  Bconsy_tmp = sdet*Bvecy_tmp 
+                  Bconsz_tmp = sdet*Bvecz_tmp 
+
+                  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, &
+                       Bconsx_tmp, Bconsy_tmp, Bconsz_tmp,rho_tmp,&
+                       velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,&
+                       Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,&
+                       g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), &
+                       uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                       epsnegative,GRHydro_C2P_failed(i,j,k))
+
+                  if(GRHydro_C2P_failed(i,j,k).ne.0) then
+                    GRHydro_C2P_failed(i,j,k) = 0
+                    call prim2conAM(GRHydro_eos_handle,g11(i,j,k),g12(i,j,k), &
+                         g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k),det, &
+                       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),vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3), &
+                          eps(i,j,k),press(i,j,k), &
+                  Bprim(i,j,k,1),Bprim(i,j,k,2),Bprim(i,j,k,3),w_lorentz(i,j,k))
+                    cycle
+                  end if
+                end if
+
+                bdotv=g11(i,j,k)*Bprim(i,j,k,1)*vup(i,j,k,1)+ &
+                      g22(i,j,k)*Bprim(i,j,k,2)*vup(i,j,k,2)+ &
+                      g33(i,j,k)*Bprim(i,j,k,3)*vup(i,j,k,3)+ &
+                 2.0*(g12(i,j,k)*Bprim(i,j,k,1)*vup(i,j,k,2)+ &
+                      g13(i,j,k)*Bprim(i,j,k,1)*vup(i,j,k,3)+ &
+                      g23(i,j,k)*Bprim(i,j,k,2)*vup(i,j,k,3))
+
+                magpress = 0.5d0*(b2/w_lorentz(i,j,k)**2+bdotv**2)
+
+                if(rho(i,j,k)*eps(i,j,k)*max_magnetic_to_gas_pressure_ratio.le.magpress) then
+                  GRHydro_C2P_failed(i,j,k) = 0
+
+                  if(evolve_entropy.ne.0) then
+                    local_K = entropycons(i,j,k)/dens(i,j,k) 
+                    xrho=10.0d0; xeps=1.0d0
+                    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)
+                  end if
+
+                  rho_tmp = rho(i,j,k)
+                  press_tmp = press(i,j,k)
+                  eps_tmp = eps(i,j,k)
+                  velx_tmp = vup(i,j,k,1)
+                  vely_tmp = vup(i,j,k,2)
+                  velz_tmp = vup(i,j,k,3)
+                  w_lorentz_tmp = w_lorentz(i,j,k)
+                  Bvecx_tmp = Bprim(i,j,k,1)
+                  Bvecy_tmp = Bprim(i,j,k,2)
+                  Bvecz_tmp = Bprim(i,j,k,3)
+                  Bconsx_tmp = sdet*Bvecx_tmp
+                  Bconsy_tmp = sdet*Bvecy_tmp
+                  Bconsz_tmp = sdet*Bvecz_tmp
+
+                  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, &
+                       Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rho_tmp,&
+                       velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,&
+                       Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,&
+                       g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), &
+                       uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                       epsnegative,GRHydro_C2P_failed(i,j,k))
+
+                  rho(i,j,k) = rho_tmp 
+                  press(i,j,k) = press_tmp 
+                  eps(i,j,k) = eps_tmp 
+                  vup(i,j,k,1) = velx_tmp 
+                  vup(i,j,k,2) = vely_tmp 
+                  vup(i,j,k,3) = velz_tmp 
+                  w_lorentz(i,j,k) = w_lorentz_tmp 
+                  Bprim(i,j,k,1) = Bvecx_tmp 
+                  Bprim(i,j,k,2) = Bvecy_tmp 
+                  Bprim(i,j,k,3) = Bvecz_tmp 
+                  cycle
+                end if
+              end if
+
+           else    ! if(evolve_temper.eq.0) then
+
+              rho_tmp = rho(i,j,k)
+              press_tmp = press(i,j,k)
+              eps_tmp = eps(i,j,k)
+              velx_tmp = vup(i,j,k,1)
+              vely_tmp = vup(i,j,k,2)
+              velz_tmp = vup(i,j,k,3)
+              w_lorentz_tmp = w_lorentz(i,j,k)
+              Bvecx_tmp = Bprim(i,j,k,1)
+              Bvecy_tmp = Bprim(i,j,k,2)
+              Bvecz_tmp = Bprim(i,j,k,3)
+              Bconsx_tmp = sdet*Bvecx_tmp
+              Bconsy_tmp = sdet*Bvecy_tmp
+              Bconsz_tmp = sdet*Bvecz_tmp
+
+              keytemp = 0 
+
+              call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, &
+                   GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), &
+                   scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), &
+                   Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, Y_e(i,j,k), &
+                   temperature(i,j,k),rho_tmp,velx_tmp,vely_tmp,velz_tmp,&
+                   eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,&
+                   w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),&
+                   g22(i,j,k),g23(i,j,k),g33(i,j,k), &
+                   uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                   epsnegative,GRHydro_C2P_failed(i,j,k))
+              if(GRHydro_C2P_failed(i,j,k).ne.0) then
+                ! this means c2p did not converge.
+                ! In this case, we attempt to call c2p with a reduced
+                ! accuracy requirement; if it fails again, we abort
+                GRHydro_C2P_failed(i,j,k) = 0
+                local_perc_ptol = GRHydro_eos_rf_prec*100.0d0
+                ! Use the previous primitive values as initial guesses
+                rho_tmp = rho(i,j,k)
+                press_tmp = press(i,j,k)
+                eps_tmp = eps(i,j,k)
+                velx_tmp = vup(i,j,k,1)
+                vely_tmp = vup(i,j,k,2)
+                velz_tmp = vup(i,j,k,3)
+                w_lorentz_tmp = w_lorentz(i,j,k)
+                Bvecx_tmp = Bprim(i,j,k,1)
+                Bvecy_tmp = Bprim(i,j,k,2)
+                Bvecz_tmp = Bprim(i,j,k,3)
+                Bconsx_tmp = sdet*Bvecx_tmp
+                Bconsy_tmp = sdet*Bvecy_tmp
+                Bconsz_tmp = sdet*Bvecz_tmp
+                call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, &
+                     local_perc_ptol, local_gam(1), dens(i,j,k), &
+                     scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), &
+                     Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, Y_e(i,j,k), &
+                     temperature(i,j,k),rho_tmp,velx_tmp,vely_tmp,velz_tmp,&
+                     eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,&
+                     w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),&
+                     g22(i,j,k),g23(i,j,k),g33(i,j,k), &
+                     uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                     epsnegative,GRHydro_C2P_failed(i,j,k))
+                if(GRHydro_C2P_failed(i,j,k).ne.0) then
+                  !$OMP CRITICAL
+                  if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
+                     call CCTK_WARN(1,"Convergence problem in c2p")
+                     write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel
+                     call CCTK_WARN(1,warnline)
+                     write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k)
+                     call CCTK_WARN(1,warnline)
+                     write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k),&
+                          temperature(i,j,k),y_e(i,j,k)
+                     call CCTK_WARN(1,warnline)
+                     call CCTK_WARN(0,"Aborting!!!")
+                  endif
+                  !$OMP END CRITICAL
+                endif
+              endif
+           endif    ! if(evolve_temper.eq.0) then
+
+           
+           if (epsnegative .ne. 0) then  
+              
+#if 0
+              ! cott 2010/03/30:
+              ! Set point to atmosphere, but continue evolution -- this is better than using
+              ! the poly EOS -- it will lead the code to crash if this happens inside a (neutron) star,
+              ! but will allow the job to continue if it happens in the atmosphere or in a
+              ! zone that contains garbage (i.e. boundary, buffer zones)
+              ! Ultimately, we want this fixed via a new carpet mask presently under development
+              !           GRHydro_C2P_failed(i,j,k) = 1
+              
+              !$OMP CRITICAL
+              call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0! ')
+              write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
+              call CCTK_WARN(GRHydro_NaN_verbose+2,warnline)
+              write(warnline,'(a20,3g16.7)') 'xyz location: ',&
+                   x(i,j,k),y(i,j,k),z(i,j,k)
+              call CCTK_WARN(GRHydro_NaN_verbose+2,warnline)
+              write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k)
+              call CCTK_WARN(GRHydro_NaN_verbose+2,warnline)
+              call CCTK_WARN(GRHydro_NaN_verbose+2,"Setting the point to atmosphere")
+              !$OMP END CRITICAL
+              
+              ! for safety, let's set the point to atmosphere
+              dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
+              rho(i,j,k) = GRHydro_rho_min
+              scon(i,j,k,:) = 0.d0
+              vup(i,j,k,:) = 0.d0
+              w_lorentz(i,j,k) = 1.d0
+
+              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)
+
+              b2=g11(i,j,k)*Bprim(i,j,k,1)**2+g22(i,j,k)*Bprim(i,j,k,2)**2+g33(i,j,k)*Bprim(i,j,k,3)**2+ &
+                   2.0*(g12(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,2)+g13(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,3)+ &
+                   g23(i,j,k)*Bprim(i,j,k,2)*Bprim(i,j,k,3))
+              
+              
+              ! w_lorentz=1, so the expression for tau reduces to [see above]:
+              tau(i,j,k)  = sqrt(det) * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k)           
+#else
+              ! cott 2010/03/27:      
+              ! Honestly, this should never happen. We need to flag the point where
+              ! this happened as having led to failing con2prim.
+              
+              !$OMP CRITICAL
+              call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype.')
+              !$OMP END CRITICAL
+              
+              xrho=1.0d0; xtemp=0.0d0; xeps=1.0d0
+              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=10.0d0; xeps=1.0d0
+              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)
+
+              rho_tmp = rho(i,j,k)
+              press_tmp = press(i,j,k)
+              eps_tmp = eps(i,j,k)
+              velx_tmp = vup(i,j,k,1)
+              vely_tmp = vup(i,j,k,2)
+              velz_tmp = vup(i,j,k,3)
+              w_lorentz_tmp = w_lorentz(i,j,k)
+              Bvecx_tmp = Bprim(i,j,k,1)
+              Bvecy_tmp = Bprim(i,j,k,2)
+              Bvecz_tmp = Bprim(i,j,k,3)
+              Bconsx_tmp = sdet*Bvecx_tmp
+              Bconsy_tmp = sdet*Bvecy_tmp
+              Bconsz_tmp = sdet*Bvecz_tmp
+
+            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, &
+                   Bconsx_tmp, Bconsy_tmp, Bconsz_tmp,rho_tmp,&
+                   velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,&
+                   Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,&
+                   g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), &
+                   uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                   epsnegative,GRHydro_C2P_failed(i,j,k))
+
+#endif          
+              
+           end if    ! if (epsnegative .ne. 0) then  
+
+           rho(i,j,k) = rho_tmp 
+           press(i,j,k) = press_tmp 
+           eps(i,j,k) = eps_tmp 
+           vup(i,j,k,1) = velx_tmp 
+           vup(i,j,k,2) = vely_tmp 
+           vup(i,j,k,3) = velz_tmp 
+           w_lorentz(i,j,k) = w_lorentz_tmp 
+           Bprim(i,j,k,1) = Bvecx_tmp 
+           Bprim(i,j,k,2) = Bvecy_tmp 
+           Bprim(i,j,k,3) = Bvecz_tmp 
+          
+           if ( rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance)) then
+!           if ( rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) .or. GRHydro_C2P_failed(i,j,k) .ge. 1) then
+               
+              b2=g11(i,j,k)*Bprim(i,j,k,1)**2+g22(i,j,k)*Bprim(i,j,k,2)**2+g33(i,j,k)*Bprim(i,j,k,3)**2+ &
+                   2.0*(g12(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,2)+g13(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,3)+ &
+                   g23(i,j,k)*Bprim(i,j,k,2)*Bprim(i,j,k,3))
+              
+              dens(i,j,k) = sdet*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
+              rho(i,j,k) = GRHydro_rho_min
+              scon(i,j,k,:) = 0.d0
+              vup(i,j,k,:) = 0.d0
+              w_lorentz(i,j,k) = 1.d0
+
+              if(evolve_temper.ne.0) then
+                ! set hot atmosphere values
+                temperature(i,j,k) = grhydro_hot_atmo_temp
+                y_e(i,j,k) = grhydro_hot_atmo_Y_e
+                y_e_con(i,j,k) = y_e(i,j,k) * dens(i,j,k)
+                keytemp = 1
+                call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
+                   press(i,j,k),keyerr,anyerr)
+              else
+                keytemp = 0
+                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)
+              endif
+
+              ! w_lorentz=1, so the expression for tau reduces to:
+
+              !!$ tau does need to take into account the existing B-field
+              !!$ with w_lorentz=1, we find tau = sqrtdet*(rho (1+eps+b^2/2)) - dens  [Press drops out]
+              tau(i,j,k)  = sdet * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k)
+
+              if(tau(i,j,k).le.sdet*b2*0.5d0)then
+                tau(i,j,k) = GRHydro_tau_min + sdet*b2*0.5d0
+              endif
+
+           end if
+
+        end do
+     end do
+  end do
+  !$OMP END PARALLEL DO
+  
+  return
+#undef faulty_gxx
+#undef faulty_gxy
+#undef faulty_gxz
+#undef faulty_gyy
+#undef faulty_gyz
+#undef faulty_gzz
+#undef faulty_betax
+#undef faulty_betay
+#undef faulty_betaz
+#undef faulty_vel
+#undef faulty_Bvec
+  
+end subroutine Conservative2PrimitiveAM
+
+
+ /*@@
+   @routine    Conservative2PrimitiveBoundariesAM
+   @date       Sep 15, 2010
+   @author     Scott Noble, Joshua Faber, Bruno Mundim, The GRHydro Developers    
+   @desc 
+        This routine is used only if the reconstruction is performed on the conserved variables. 
+        It computes the primitive variables on cell boundaries.
+        Since reconstruction on conservative had not proved to be very successful, 
+        some of the improvements to the C2P routines (e.g. the check about 
+        whether a failure happens in a point that will be restriced anyway) 
+        are not implemented here yet.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+
+subroutine Conservative2PrimitiveBoundsAM(CCTK_ARGUMENTS)
+  
+  use Con2PrimM_fortran_interfaces
+  
+  implicit none
+  
+  ! save memory when MP is not used
+  ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1
+  TARGET gaa, gab, gac, gbb, gbc, gcc
+  TARGET gxx, gxy, gxz, gyy, gyz, gzz
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+  
+  integer :: i, j, k, itracer, nx, ny, nz
+  CCTK_REAL :: uxxl, uxyl, uxzl, uyyl, uyzl, uzzl,&
+       uxxr, uxyr, uxzr, uyyr, uyzr, uzzr, pmin(1), epsmin(1)
+  CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
+       gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
+  CCTK_REAL :: b2minus, b2plus, local_gam(1), local_pgam,local_K,scminus,scplus
+  CCTK_INT :: epsnegative
+  CCTK_REAL :: Bconsx_tmp, Bconsy_tmp, Bconsz_tmp
+  character(len=100) warnline
+ 
+  CCTK_REAL :: local_min_tracer
+
+! begin EOS Omni vars                                            
+  integer :: n,keytemp,anyerr,keyerr(1)
+  real*8  :: xpress(1),xtemp(1),xye(1),xeps(1),xrho(1),one(1)=1.0d0
+! end EOS Omni vars                       
+
+  ! save memory when MP is not used
+  CCTK_INT :: GRHydro_UseGeneralCoordinates
+  CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33
+
+  if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
+    g11 => gaa
+    g12 => gab
+    g13 => gac
+    g22 => gbb
+    g23 => gbc
+    g33 => gcc
+  else
+    g11 => gxx
+    g12 => gxy
+    g13 => gxz
+    g22 => gyy
+    g23 => gyz
+    g33 => gzz
+  end if
+#define gxx faulty_gxx
+#define gxy faulty_gxy
+#define gxz faulty_gxz
+#define gyy faulty_gyy
+#define gyz faulty_gyz
+#define gzz faulty_gzz
+#define betax faulty_betax
+#define betay faulty_betay
+#define betaz faulty_betaz
+#define vel faulty_vel
+#define Bvec faulty_Bvec
+
+! begin EOS Omni vars                                            
+  n=1;keytemp=0;anyerr=0;keyerr(1)=0
+  xpress(1)=0.0d0;xeps(1)=0.0d0;xtemp(1)=0.0d0;xye(1)=0.0d0
+! end EOS Omni vars                       
+
+  ! this is a poly call
+  xrho(1)=GRHydro_rho_min 
+  call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       xrho,one,xtemp,xye,pmin,keyerr,anyerr)
+
+  call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+       xrho,epsmin,xtemp,xye,pmin,epsmin,keyerr,anyerr)
+
+  call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+         one,one,xtemp,xye,local_gam,keyerr,anyerr)
+  local_gam=local_gam+1.0
+
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3) 
+  
+  if (use_min_tracer .ne. 0) then
+     local_min_tracer = min_tracer
+  else
+     local_min_tracer = 0d0
+  end if
+
+  do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+     do j = GRHydro_stencil, ny - GRHydro_stencil + 1
+        do i = GRHydro_stencil, nx - GRHydro_stencil + 1
+           
+           !do not compute if in atmosphere or in an excised region
+           if ((atmosphere_mask(i,j,k) .ne. 0) .or. &
+                GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
+           
+           gxxl = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset))
+           gxyl = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset))
+           gxzl = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset))
+           gyyl = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset))
+           gyzl = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset))
+           gzzl = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset))
+           gxxr = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset))
+           gxyr = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset))
+           gxzr = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset))
+           gyyr = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset))
+           gyzr = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset))
+           gzzr = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset))
+           
+           epsnegative = 0
+           
+           avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
+           avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
+           call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,&
+                gxxl,gxyl,gxzl,gyyl,gyzl,gzzl)        
+           call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,&
+                gxxr,gxyr,gxzr,gyyr,gyzr,gzzr)        
+           
+!!$ Tracers get no update for MHD!
+           if (evolve_tracer .ne. 0) then
+              do itracer=1,number_of_tracers
+                 call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), &
+                      tracer(i,j,k,itracer), dens(i,j,k))
+              enddo
+              
+              if (use_min_tracer .ne. 0) then
+                 if (tracer(i,j,k,itracer) .le. local_min_tracer) then
+                    tracer(i,j,k,itracer) = local_min_tracer
+                 end if
+              end if
+              
+           endif
+
+           if(evolve_Y_e.ne.0) then
+              Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
+           endif
+           Bconsx_tmp = sqrt(avg_detl)*Bvecxminus(i,j,k)
+           Bconsy_tmp = sqrt(avg_detl)*Bvecyminus(i,j,k)
+           Bconsz_tmp = sqrt(avg_detl)*Bveczminus(i,j,k)
+           call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1),densminus(i,j,k), &
+                sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k), &
+                Bconsx_tmp, Bconsy_tmp, Bconsz_tmp,  xye(1), xtemp(1), rhominus(i,j,k),&
+                velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),&
+                Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),&
+                gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
+                uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, &
+                epsnegative,GRHydro_C2P_failed(i,j,k))
+           
+           if (epsnegative .ne. 0) then
+              !$OMP CRITICAL
+              call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype!')
+              !$OMP END CRITICAL
+              
+              xrho=10.0d0
+              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   one,one,xtemp,xye,xpress,keyerr,anyerr)
+              local_K = xpress(1)
+
+              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   xrho,one,xtemp,xye,xpress,keyerr,anyerr)
+              local_pgam=log(xpress(1)/local_K)/log(xrho(1))
+              scminus = 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), scminus, &
+                   Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rhominus(i,j,k),&
+                   velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),&
+                   Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),&
+                   gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
+                   uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, &
+                   epsnegative,GRHydro_C2P_failed(i,j,k))
+          
+           end if
+           
+           if (epsminus(i,j,k) .lt. 0.0d0) then
+              if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
+                 !$OMP CRITICAL
+                 call CCTK_WARN(1,'Con2Prim: stopping the code.')
+                 call CCTK_WARN(1, '   specific internal energy just went below 0! ')
+                 write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
+                 call CCTK_WARN(1,warnline)
+                 write(warnline,'(a20,3g16.7)') 'xyz location: ',&
+                      x(i,j,k),y(i,j,k),z(i,j,k)
+                 call CCTK_WARN(1,warnline)
+                 write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k)
+                 call CCTK_WARN(1,warnline)
+                 write(warnline,'(a20,3g16.7)') 'velocities: ',&
+                      velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k)
+                 call CCTK_WARN(1,warnline)
+                 call CCTK_WARN(GRHydro_c2p_warnlevel, "Specific internal energy negative")
+                 !$OMP END CRITICAL
+                 exit
+              endif
+           endif
+           
+           epsnegative = 0
+
+           Bconsx_tmp = sqrt(avg_detr)*Bvecxplus(i,j,k)
+           Bconsy_tmp = sqrt(avg_detr)*Bvecyplus(i,j,k)
+           Bconsz_tmp = sqrt(avg_detr)*Bveczplus(i,j,k)
+           call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1),densplus(i,j,k), &
+                sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), tauplus(i,j,k), &
+                Bconsx_tmp, Bconsy_tmp, Bconsz_tmp,  xye(1), xtemp(1), rhoplus(i,j,k),&
+                velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
+                Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k),b2plus, w_lorentzplus(i,j,k),&
+                gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, &
+                uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, &
+                epsnegative,GRHydro_C2P_failed(i,j,k))
+ 
+           if (epsnegative .ne. 0) then
+              !$OMP CRITICAL
+              call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype!!')
+              !$OMP END CRITICAL
+
+              xrho=10.0d0
+              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   one,one,xtemp,xye,xpress,keyerr,anyerr)
+              local_K = xpress(1)
+
+              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   xrho,one,xtemp,xye,xpress,keyerr,anyerr)
+              local_pgam=log(xpress(1)/local_K)/log(xrho(1))
+              scplus = 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), scplus,&
+                   Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rhoplus(i,j,k),&
+                   velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
+                   Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k),b2plus, w_lorentzplus(i,j,k),&
+                   gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, &
+                   uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, &
+                   epsnegative,GRHydro_C2P_failed(i,j,k))
+           end if
+           
+           if (epsplus(i,j,k) .lt. 0.0d0) then
+              if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
+                 !$OMP CRITICAL
+                 call CCTK_WARN(1,'Con2Prim: stopping the code.')
+                 call CCTK_WARN(1, '   specific internal energy just went below 0! ')
+                 write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
+                 call CCTK_WARN(1,warnline)
+                 write(warnline,'(a20,3g16.7)') 'xyz location: ',&
+                      x(i,j,k),y(i,j,k),z(i,j,k)
+                 call CCTK_WARN(1,warnline)
+                 write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k)
+                 call CCTK_WARN(1,warnline)
+                 write(warnline,'(a20,3g16.7)') 'velocities: ',&
+                      velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k)
+                 call CCTK_WARN(1,warnline)
+                 call CCTK_WARN(GRHydro_c2p_warnlevel, "Specific internal energy negative")
+                 write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',&
+                      x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k)
+                 call CCTK_WARN(1,warnline)
+                 !$OMP END CRITICAL
+              endif
+           endif
+        end do
+     end do
+  end do
+  
+#undef faulty_gxx
+#undef faulty_gxy
+#undef faulty_gxz
+#undef faulty_gyy
+#undef faulty_gyz
+#undef faulty_gzz
+#undef faulty_betax
+#undef faulty_betay
+#undef faulty_betaz
+#undef faulty_vel
+#undef faulty_Bvec
+  
+end subroutine Conservative2PrimitiveBoundsAM
+
+/*@@
+@routine    Con2PrimPolytypeAM
+@date       Sep 16, 2010
+@author     SCott Noble, Joshua Faber, Bruno Mundim, Ian Hawke
+@desc 
+All routines below are identical to those above, just
+specialised from polytropic type EOS.
+@enddesc 
+@calls     
+@calledby   
+@history 
+
+@endhistory 
+
+@@*/
+
+subroutine Conservative2PrimitivePolytypeAM(CCTK_ARGUMENTS)
+  
+  implicit none
+  
+  ! save memory when MP is not used
+  ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1
+  TARGET gaa, gab, gac, gbb, gbc, gcc
+  TARGET gxx, gxy, gxz, gyy, gyz, gzz
+  TARGET lvel, vel
+  TARGET lBvec, Bvec
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  
+  integer :: i, j, k, itracer, nx, ny, nz
+  CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det,b2
+
+  CCTK_INT :: epsnegative
+  CCTK_REAL :: Bconsx_tmp, Bconsy_tmp, Bconsz_tmp
+  
+  CCTK_REAL :: local_min_tracer, 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
+! end EOS Omni vars                 
+
+  ! save memory when MP is not used
+  CCTK_INT :: GRHydro_UseGeneralCoordinates
+  CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33
+  CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim
+
+  call GRHydro_BvecfromAvec(CCTK_PASS_FTOF)
+
+  if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
+    g11 => gaa
+    g12 => gab
+    g13 => gac
+    g22 => gbb
+    g23 => gbc
+    g33 => gcc
+    vup => lvel
+    Bprim => lBvec
+  else
+    g11 => gxx
+    g12 => gxy
+    g13 => gxz
+    g22 => gyy
+    g23 => gyz
+    g33 => gzz
+    vup => vel
+    Bprim => Bvec
+  end if
+#define gxx faulty_gxx
+#define gxy faulty_gxy
+#define gxz faulty_gxz
+#define gyy faulty_gyy
+#define gyz faulty_gyz
+#define gzz faulty_gzz
+#define betax faulty_betax
+#define betay faulty_betay
+#define betaz faulty_betaz
+#define vel faulty_vel
+#define Bvec faulty_Bvec
+
+! begin EOS Omni vars                                                                                       
+  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                 
+  
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+  
+  if (use_min_tracer .ne. 0) then
+     local_min_tracer = min_tracer
+  else
+     local_min_tracer = 0d0
+  end if
+  
+!!$  do k = GRHydro_stencil + 1, nz - GRHydro_stencil
+!!$    do j = GRHydro_stencil + 1, ny - GRHydro_stencil
+!!$      do i = GRHydro_stencil + 1, nx - GRHydro_stencil
+
+
+  !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,&
+  !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, &
+  !$OMP Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, &
+  !$OMP b2, xrho, xpress, local_K, local_pgam, sc)
+  do k = 1, nz
+     do j = 1, ny
+        do i = 1, nx
+           
+           !do not compute if in atmosphere or in an excised region
+           if ((atmosphere_mask(i,j,k) .ne. 0) .or. &
+                GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
+           
+        det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k))
+           call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,&
+                g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),&
+                g23(i,j,k),g33(i,j,k))        
+
+!!$ No MHD changes to tracers
+           if (evolve_tracer .ne. 0) then
+              do itracer=1,number_of_tracers
+                 call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), & 
+                      tracer(i,j,k,itracer), dens(i,j,k))
+              enddo
+              
+              if (use_min_tracer .ne. 0) then
+                 if (tracer(i,j,k,itracer) .le. local_min_tracer) then
+                    tracer(i,j,k,itracer) = local_min_tracer
+                 end if
+              end if
+              
+           endif
+
+           if(evolve_Y_e.ne.0) then
+              Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
+           endif
+           
+           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)
+           Bconsx_tmp = sqrt(det)*Bprim(i,j,k,1)
+           Bconsy_tmp = sqrt(det)*Bprim(i,j,k,2)
+           Bconsz_tmp = sqrt(det)*Bprim(i,j,k,3)
+           
+           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, &
+                   Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rho(i,j,k),&
+                   vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3),eps(i,j,k),press(i,j,k),&
+                   Bprim(i,j,k,1), Bprim(i,j,k,2), Bprim(i,j,k,3),b2, w_lorentz(i,j,k),&
+                   g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), &
+                   uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                   epsnegative,GRHydro_C2P_failed(i,j,k))
+            
+        end do
+     end do
+  end do
+  
+  !$OMP END PARALLEL DO
+  
+  return
+
+#undef faulty_gxx
+#undef faulty_gxy
+#undef faulty_gxz
+#undef faulty_gyy
+#undef faulty_gyz
+#undef faulty_gzz
+#undef faulty_betax
+#undef faulty_betay
+#undef faulty_betaz
+#undef faulty_vel
+#undef faulty_Bvec
+  
+end subroutine Conservative2PrimitivePolytypeAM
+
+
+ /*@@
+   @routine    Cons2PrimBoundsPolytypeAM
+   @date       Sep 16, 2010
+   @author     Scott Noble, Joshua Faber, Bruno Mundim, The GRHydro Developers
+   @desc 
+        This routine is used only if the reconstruction is performed on the conserved variables. 
+        It computes the primitive variables on cell boundaries.
+        Since reconstruction on conservative had not proved to be very successful, 
+        some of the improvements to the C2P routines (e.g. the check about 
+        whether a failure happens in a point that will be restriced anyway) are not implemented here yet.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine Con2PrimBoundsPolytypeAM(CCTK_ARGUMENTS)
+  
+  implicit none
+  
+  ! save memory when MP is not used
+  ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1
+  TARGET gaa, gab, gac, gbb, gbc, gcc
+  TARGET gxx, gxy, gxz, gyy, gyz, gzz
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  
+  integer :: i, j, k, nx, ny, nz
+  CCTK_REAL :: uxxl, uxyl, uxzl, uyyl, uyzl, uzzl,&
+       uxxr, uxyr, uxzr, uyyr, uyzr, uzzr
+  CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
+       gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
+  CCTK_REAL :: b2minus, b2plus
+  CCTK_REAL :: Bconsx_tmp, Bconsy_tmp, Bconsz_tmp
+  CCTK_INT :: epsnegative
+
+  CCTK_REAL :: local_pgam,local_K,scplus,scminus
+
+! begin EOS Omni vars                                                                                       
+  integer :: n,keytemp,anyerr,keyerr(1)
+  real*8  :: xpress,xtemp,xye,xeps,xrho
+! end EOS Omni vars                       
+
+  ! save memory when MP is not used
+  CCTK_INT :: GRHydro_UseGeneralCoordinates
+  CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33
+
+  if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
+    g11 => gaa
+    g12 => gab
+    g13 => gac
+    g22 => gbb
+    g23 => gbc
+    g33 => gcc
+  else
+    g11 => gxx
+    g12 => gxy
+    g13 => gxz
+    g22 => gyy
+    g23 => gyz
+    g33 => gzz
+  end if
+#define gxx faulty_gxx
+#define gxy faulty_gxy
+#define gxz faulty_gxz
+#define gyy faulty_gyy
+#define gyz faulty_gyz
+#define gzz faulty_gzz
+#define betax faulty_betax
+#define betay faulty_betay
+#define betaz faulty_betaz
+#define vel faulty_vel
+#define Bvec faulty_Bvec
+
+! begin EOS Omni vars                                            
+  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                 
+  
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+  
+  do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+     do j = GRHydro_stencil, ny - GRHydro_stencil + 1
+        do i = GRHydro_stencil, nx - GRHydro_stencil + 1
+           
+           !do not compute if in atmosphere or in an excised region
+           if ((atmosphere_mask(i,j,k) .ne. 0) .or. &
+                GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
+           
+           gxxl = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset))
+           gxyl = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset))
+           gxzl = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset))
+           gyyl = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset))
+           gyzl = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset))
+           gzzl = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset))
+           gxxr = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset))
+           gxyr = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset))
+           gxzr = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset))
+           gyyr = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset))
+           gyzr = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset))
+           gzzr = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset))
+           
+           avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
+           avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
+           call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,&
+                gxxl,gxyl,gxzl,gyyl,gyzl,gzzl)        
+           call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,&
+                gxxr,gxyr,gxzr,gyyr,gyzr,gzzr)        
+
+           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)
+           Bconsx_tmp = sqrt(avg_detl)*Bvecxminus(i,j,k)
+           Bconsy_tmp = sqrt(avg_detl)*Bvecyminus(i,j,k)
+           Bconsz_tmp = sqrt(avg_detl)*Bveczminus(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,&
+                Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rhominus(i,j,k),&
+                velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),&
+                Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),&
+                gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
+                uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, &
+                epsnegative,GRHydro_C2P_failed(i,j,k))
+
+           Bconsx_tmp = sqrt(avg_detr)*Bvecxplus(i,j,k)
+           Bconsy_tmp = sqrt(avg_detr)*Bvecyplus(i,j,k)
+           Bconsz_tmp = sqrt(avg_detr)*Bveczplus(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,&
+                Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rhoplus(i,j,k),&
+                velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
+                Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,w_lorentzplus(i,j,k),&
+                gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, &
+                uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, &
+                epsnegative,GRHydro_C2P_failed(i,j,k))
+        end do
+     end do
+  end do
+  
+#undef faulty_gxx
+#undef faulty_gxy
+#undef faulty_gxz
+#undef faulty_gyy
+#undef faulty_gyz
+#undef faulty_gzz
+#undef faulty_betax
+#undef faulty_betay
+#undef faulty_betaz
+#undef faulty_vel
+#undef faulty_Bvec
+  
+end subroutine Con2PrimBoundsPolytypeAM
+
+!!$ Con2Prim_ptTracer, Con2Prim_BoundsTracer, and Con2Prim_ptBoundsTracer need not be rewritten!
+