/*@@
   @file      GRHydro_Con2PrimM.F90
   @date      Sep 3, 2010
   @author    Scott Noble, Joshua Faber, Bruno Mundim
   @desc 
   The routines for converting conservative to primitive variables.
   @enddesc 
 @@*/

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "SpaceMask.h"
#include "GRHydro_InterfacesM.h"
#include "GRHydro_Macros.h"

 /*@@
   @routine    Conservative2PrimitiveM
   @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 Conservative2PrimitiveM(CCTK_ARGUMENTS)

  use Con2PrimM_fortran_interfaces

  implicit none
  
  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, pmin(1), epsmin(1)
  CCTK_REAL :: b2
  CCTK_INT :: epsnegative
  character(len=100) warnline
  
  CCTK_REAL :: local_min_tracer, local_gam(1), 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),one(1)=1.0d0

  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

  ! 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,xpress,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

  !$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 )
  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(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 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))        
           
!!$ 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) = 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
              
              b2=gxx(i,j,k)*Bvec(i,j,k,1)**2+gyy(i,j,k)*Bvec(i,j,k,2)**2+gzz(i,j,k)*Bvec(i,j,k,3)**2+ &
                   2.0*(gxy(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,2)+gxz(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,3)+ &
                   gyz(i,j,k)*Bvec(i,j,k,2)*Bvec(i,j,k,3))
              
              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
              vel(i,j,k,:) = 0.d0
              w_lorentz(i,j,k) = 1.d0
              xtemp = 0.0d0

              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)  = sqrt(det) * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k)
              
              cycle
              
           end if

           if(evolve_temper.eq.0) then
              call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, 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), &
                   Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), 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),&
                   Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, 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, &
                   epsnegative,GRHydro_C2P_failed(i,j,k))
           else
              stop "Please implement finite T con2prim routine in MHD part!"
           endif

           
           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
              vel(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=gxx(i,j,k)*Bvec(i,j,k,1)**2+gyy(i,j,k)*Bvec(i,j,k,2)**2+gzz(i,j,k)*Bvec(i,j,k,3)**2+ &
                   2.0*(gxy(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,2)+gxz(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,3)+ &
                   gyz(i,j,k)*Bvec(i,j,k,2)*Bvec(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)

              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, &
                   Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), 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),&
                   Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, 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, &
                   epsnegative,GRHydro_C2P_failed(i,j,k))
             
#endif          
              
           end if
           
           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=gxx(i,j,k)*Bvec(i,j,k,1)**2+gyy(i,j,k)*Bvec(i,j,k,2)**2+gzz(i,j,k)*Bvec(i,j,k,3)**2+ &
                   2.0*(gxy(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,2)+gxz(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,3)+ &
                   gyz(i,j,k)*Bvec(i,j,k,2)*Bvec(i,j,k,3))
              
              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
              vel(i,j,k,:) = 0.d0
              w_lorentz(i,j,k) = 1.d0
              xtemp = 0.0d0
   
              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)  = sqrt(det) * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k)

           end if
            
        end do
     end do
  end do
  !$OMP END PARALLEL DO
  
  return
  
end subroutine Conservative2PrimitiveM


 /*@@
   @routine    Conservative2PrimitiveBoundariesM
   @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 Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
  
  use Con2PrimM_fortran_interfaces
  
  implicit none
  
  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
  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
  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 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
           gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
           gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
           gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
           gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
           gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
           gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
           gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
           gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
           gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
           gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
           gzzr = 0.5d0 * (gzz(i,j,k) + gzz(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
           
           
           call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam(1),densminus(i,j,k), &
                sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k), &
                Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(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),&
                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, &
                   Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(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),&
                   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
           call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam(1), densplus(i,j,k), &
                sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), tauplus(i,j,k),&
                Bconsxplus(i,j,k), Bconsyplus(i,j,k), Bconszplus(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),&
                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,&
                   Bconsxplus(i,j,k), Bconsyplus(i,j,k), Bconszplus(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),&
                   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
  
end subroutine Conservative2PrimitiveBoundsM

/*@@
@routine    Con2PrimPolytypeM
@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 Conservative2PrimitivePolytypeM(CCTK_ARGUMENTS)
  
  implicit none
  
  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 :: 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                 
  
  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 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(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 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))        

!!$ 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)
           
           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, &
                   Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), 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),&
                   Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, 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, &
                   epsnegative,GRHydro_C2P_failed(i,j,k))
            
        end do
     end do
  end do
  
  !$OMP END PARALLEL DO
  
  return

end subroutine Conservative2PrimitivePolytypeM


 /*@@
   @routine    Cons2PrimBoundsPolytypeM
   @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 Con2PrimBoundsPolytypeM(CCTK_ARGUMENTS)
  
  implicit none
  
  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_INT :: 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                 
  
  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 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
           gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
           gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
           gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
           gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
           gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
           gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
           gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
           gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
           gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
           gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
           gzzr = 0.5d0 * (gzz(i,j,k) + gzz(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)
           
           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,&
                Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(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),&
                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))

           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,&
                Bconsxplus(i,j,k), Bconsyplus(i,j,k), Bconszplus(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),&
                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
  
end subroutine Con2PrimBoundsPolytypeM

!!$ Con2Prim_ptTracer, Con2Prim_BoundsTracer, and Con2Prim_ptBoundsTracer need not be rewritten!