/*@@
   @file      primitive2conservative
   @date      Thu Jan  11 11:03:32 2002
   @author    Pedro Montero, Ian Hawke
   @desc 
   Primitive to conservative routine
   @enddesc 
 @@*/

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

#define velx(i,j,k) vel(i,j,k,1)
#define vely(i,j,k) vel(i,j,k,2)
#define velz(i,j,k) vel(i,j,k,3)
#define sx(i,j,k) scon(i,j,k,1)
#define sy(i,j,k) scon(i,j,k,2)
#define sz(i,j,k) scon(i,j,k,3)


 /*@@
   @routine   primitive2conservative.f90 
   @date       Thu Jan 11 11:03:32 2002
   @author     Pedro Montero, Ian Hawke
   @desc 
   Converts primitive to conserved variables for the boundary extended data.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine primitive2conservative(CCTK_ARGUMENTS)

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  
  integer :: i, j, k
  CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
       gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
  CCTK_REAL :: xtemp(1)

  if(evolve_temper.ne.1) then
     !$OMP PARALLEL DO PRIVATE(i, j, avg_detl, avg_detr,&
     !$OMP                      gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
     !$OMP                      gxxr,gxyr,gxzr,gyyr,gyzr,gzzr)
     do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
        do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
           do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
              
              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 prim2con(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,& 
                   gyzl,gzzl, &
                   avg_detl,densminus(i,j,k),sxminus(i,j,k),&
                   syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),&
                   rhominus(i,j,k), &
                   velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),&
                   epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k))
              
              call prim2con(GRHydro_eos_handle, gxxr,gxyr,gxzr,&
                   gyyr,gyzr,gzzr, &
                   avg_detr, densplus(i,j,k),sxplus(i,j,k),&
                   syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),&
                   rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
                   velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
                   w_lorentzplus(i,j,k)) 
              
           end do
        end do
     end do
     !$OMP END PARALLEL DO
  else
     !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,&
     !$OMP                      gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
     !$OMP                      gxxr,gxyr,gxzr,gyyr,gyzr,gzzr)
     do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
        do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
           do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
              
              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)

              ! we do not have plus/minus vars for temperature since
              ! eps is reconstructed. Hence, we do not update the
              ! variable 'temperature' in prim2con at the interfaces 
              ! We will instead use an average temperature as an initial
              ! guess.
              xtemp = 0.5d0*(temperature(i,j,k) + &
                   temperature(i-xoffset,j-yoffset,k-zoffset))
              call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,&
                   i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxxl,gxyl,gxzl,gyyl,& 
                   gyzl,gzzl, &
                   avg_detl,densminus(i,j,k),sxminus(i,j,k),&
                   syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),&
                   rhominus(i,j,k), &
                   velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),&
                   epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k),&
                   xtemp,y_e_minus(i,j,k))

              ! we do not have plus/minus vars for temperature since
              ! eps is reconstructed. Hence, we do not update the
              ! variable 'temperature' in prim2con at the interfaces
              ! We will instead use an average temperature as an initial
              ! guess.
              xtemp = 0.5d0*(temperature(i,j,k) + &
                   temperature(i+xoffset,j+yoffset,k+zoffset))
              call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel, &
                   i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxxr,gxyr,gxzr,&
                   gyyr,gyzr,gzzr, &
                   avg_detr, densplus(i,j,k),sxplus(i,j,k),&
                   syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),&
                   rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
                   velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
                   w_lorentzplus(i,j,k),xtemp, &
                   y_e_plus(i,j,k)) 
              
           end do
        end do
     end do
     !$OMP END PARALLEL DO
  endif

end subroutine primitive2conservative

 /*@@
   @routine    prim2con
   @date       Sat Jan 26 01:52:18 2002
   @author     Pedro Montero, Ian Hawke
   @desc 
   Converts from primitive to conservative at a single point
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine prim2con(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, &
     dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w) 
  
  implicit none

  DECLARE_CCTK_PARAMETERS

  CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det
  CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,&
       deps, dpress, w, vlowx, vlowy, vlowz   
  CCTK_INT :: handle

! begin EOS Omni vars
  integer :: n,keytemp,anyerr,keyerr(1) 
  real*8  :: xye,xtemp
  n = 1;keytemp = 0;anyerr = 0;keyerr(1) = 0
  xtemp = 0.0d0; xye = 0.0d0
! end EOS Omni vars

  w = 1.d0 / sqrt(1.d0 - (gxx*dvelx*dvelx + gyy*dvely*dvely + gzz &
       *dvelz*dvelz + 2*gxy*dvelx*dvely + 2*gxz*dvelx *dvelz + 2*gyz&
       *dvely*dvelz))  

  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
       drho,deps,xtemp,xye,dpress,keyerr,anyerr)  

  vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz
  vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
  vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz

  ddens = sqrt(det) * drho * w 
  dsx = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowx
  dsy = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowy
  dsz = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowz
  dtau = sqrt(det) * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens 

end subroutine prim2con


subroutine prim2con_hot(handle, GRHydro_reflevel, ii, jj, kk, &
     x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, &
     dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w, &
     temp,ye) 
  
  implicit none

  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_FUNCTIONS

  CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det
  CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho(1), dvelx, dvely, dvelz,&
       deps(1), dpress(1), w, vlowx, vlowy, vlowz   
  CCTK_REAL :: temp(1),ye(1), x, y, z
  CCTK_INT :: handle, GRHydro_reflevel, ii, jj, kk
  character(len=512) warnline

! begin EOS Omni vars
  integer :: n,keytemp,anyerr,keyerr(1)
  real*8  :: temp0(1)
  n = 1;keytemp = 0;anyerr = 0;keyerr(1) = 0
! end EOS Omni vars

  temp0 = temp
  w = 1.d0 / sqrt(1.d0 - (gxx*dvelx*dvelx + gyy*dvely*dvely + gzz &
       *dvelz*dvelz + 2*gxy*dvelx*dvely + 2*gxz*dvelx *dvelz + 2*gyz&
       *dvely*dvelz))  

  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
       drho,deps,temp,ye,dpress,keyerr,anyerr)  
  ! error handling
  if(anyerr.ne.0) then
     if(GRHydro_reflevel.lt.GRHydro_c2p_warn_from_reflevel) then
        ! in this case (coarse grid error that is hopefully restricted
        ! away), we use the average temperature between cells and call
        ! the EOS with keytemp=1
        keytemp=1
        call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
             drho,deps,temp,ye,dpress,keyerr,anyerr)  
        keytemp=0
        if(anyerr.ne.0) then
           call CCTK_WARN(1,"EOS error in prim2con_hot: lev 2")
           write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
           call CCTK_WARN(1,warnline)
           write(warnline,"(1P10E15.6)") drho,deps,temp,ye
           call CCTK_WARN(1,warnline)
           write(warnline,"(A7,i8)") "code: ",keyerr(1)
           call CCTK_WARN(1,warnline)
           write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
           call CCTK_WARN(1,warnline)
        endif
     else
        ! This is a way of recovering even on finer refinement levels:
        ! Use the average temperature at the interface instead of the
        ! reconstructed specific internal energy.
        call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!")
        write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
        call CCTK_WARN(1,warnline)
        write(warnline,"(1P10E15.6)") drho,deps,temp0,ye
        call CCTK_WARN(1,warnline)
        write(warnline,"(A7,i8)") "code: ",keyerr(1)
        call CCTK_WARN(1,warnline)
        write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
        call CCTK_WARN(1,warnline)
        keytemp=1
        temp = temp0
        call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
             drho,deps,temp,ye,dpress,keyerr,anyerr)  
        keytemp=0
        if(anyerr.ne.0) then
           call CCTK_WARN(1,"EOS error in prim2con_hot")
           write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
           call CCTK_WARN(1,warnline)
           write(warnline,"(1P10E15.6)") drho,deps,temp,ye
           call CCTK_WARN(1,warnline)
           write(warnline,"(A7,i8)") "code: ",keyerr(1)
           call CCTK_WARN(1,warnline)
           write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
           call CCTK_WARN(1,warnline)
           call CCTK_WARN(0,"Aborting!!!")
        endif
     endif
  endif


  vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz
  vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
  vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz

  ddens = sqrt(det) * drho(1) * w 
  dsx = sqrt(det) * (drho(1)*(1+deps(1))+dpress(1))*w*w * vlowx
  dsy = sqrt(det) * (drho(1)*(1+deps(1))+dpress(1))*w*w * vlowy
  dsz = sqrt(det) * (drho(1)*(1+deps(1))+dpress(1))*w*w * vlowz
  dtau = sqrt(det) * ((drho(1)*(1+deps(1))+dpress(1))*w*w - dpress(1)) - ddens 


end subroutine prim2con_hot


 /*@@
   @routine    Primitive2ConservativeCells
   @date       Sun Mar 10 21:16:20 2002
   @author     
   @desc 
   Wrapper function that converts primitive to conservative at the 
     cell centres. 
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/


subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS)

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  
  CCTK_INT :: i, j, k
  CCTK_REAL :: det
  character(len=512) :: warnline

  if(evolve_temper.ne.1) then
     !$OMP PARALLEL DO PRIVATE(i, j, det)
     do k = 1,cctk_lsh(3)
        do j = 1,cctk_lsh(2)
           do i = 1,cctk_lsh(1)
          
              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 prim2con(GRHydro_eos_handle,gxx(i,j,k),&
                   gxy(i,j,k),gxz(i,j,k),&
                   gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
                   det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
                   tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
                   eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
           end do
        end do
     end do
     !$OMP END PARALLEL DO
  else
     !$OMP PARALLEL DO PRIVATE(i, j, det)
     do k = 1,cctk_lsh(3)
        do j = 1,cctk_lsh(2)
           do i = 1,cctk_lsh(1)
          
              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 prim2con_hot(GRHydro_eos_handle,GRHydro_reflevel,i,j,k,&
                   x(i,j,k),y(i,j,k),z(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),&
                   det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
                   tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
                   eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
                   temperature(i,j,k),y_e(i,j,k))
           end do
        end do
     end do
     !$OMP END PARALLEL DO
  endif

  if(evolve_Y_e.ne.0) then
     !$OMP PARALLEL DO PRIVATE(i, j)
     do k = 1,cctk_lsh(3)
        do j = 1,cctk_lsh(2)
           do i = 1,cctk_lsh(1)
              Y_e_con(i,j,k) = Y_e(i,j,k) * dens(i,j,k)
           enddo
        enddo
     enddo
     !$OMP END PARALLEL DO
  endif


end subroutine Primitive2ConservativeCells


 /*@@
   @routine    Prim2ConservativePolytype
   @date       Tue Mar 19 22:52:21 2002
   @author     Ian Hawke
   @desc 
   Same as first routine, only for polytropes.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/


subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS)

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  
  integer :: i, j, k
  CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
       gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
  
  !$OMP PARALLEL DO PRIVATE(i, j, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
  !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr)
  do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
    do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
      do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
        
        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 prim2conpolytype(GRHydro_eos_handle, gxxl,gxyl,gxzl,&
             gyyl,gyzl,gzzl, &
             avg_detl,densminus(i,j,k),sxminus(i,j,k),&
             syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),rhominus(i,j,k), &
             velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),&
             epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k))
 
        call prim2conpolytype(GRHydro_eos_handle, gxxr,gxyr,gxzr,&
             gyyr,gyzr,gzzr, &
             avg_detr, densplus(i,j,k),sxplus(i,j,k),&
             syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),&
             rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
             velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
             w_lorentzplus(i,j,k)) 
      end do
    end do
  end do
  !$OMP END PARALLEL DO


  ! Note on Y_e: We use Y_e_plus and Y_e_minus directly
  ! in the Riemann solver. That's why it is not necessary
  ! to do a prim2con for Y_e

end subroutine Prim2ConservativePolytype

 /*@@
   @routine    prim2conpolytype
   @date       Sat Jan 26 01:52:18 2002
   @author     Pedro Montero, Ian Hawke
   @desc 
   Converts from primitive to conservative at a single point
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, &
     gzz, det, ddens, &
     dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w) 
  
  implicit none
  
  DECLARE_CCTK_PARAMETERS

  CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det
  CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,&
       deps, dpress, w_tmp, w, vlowx, vlowy, vlowz, sqrtdet   
  CCTK_INT :: handle
  character(len=256) NaN_WarnLine

! begin EOS Omni vars
  integer :: n, keytemp, anyerr, keyerr(1)
  real*8  :: xpress,xeps,xtemp,xye
  n = 1; keytemp = 0; anyerr = 0; keyerr(1) = 0
  xpress = 0.0d0; xeps = 0.0d0; xtemp = 0.0d0; xye = 0.0d0
! end EOS Omni vars
  
  w_tmp = gxx*dvelx*dvelx + gyy*dvely*dvely + gzz *dvelz*dvelz + &
          2*gxy*dvelx*dvely + 2*gxz*dvelx*dvelz + 2*gyz*dvely*dvelz
  if (w_tmp .ge. 1.d0) then
    ! In theory this should not happen, and even when accepting the fact
    ! that numerically it can, one might be tempted to set w to some large
    ! value in that case. However, this would lead to completely bogus
    ! and hard to trace wrong values below. There is no good value to
    ! choose in this case, but something small is probably the best of
    ! all bad choices.
    !$OMP CRITICAL
    write(NaN_WarnLine,'(a80,2g15.6)') 'Infinite Lorentz factor reset. rho, w_tmp: ', drho, w_tmp
    call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine)
    !$OMP END CRITICAL
    w = 1.d-20
  else
    w = 1.d0 / sqrt(1.d0 - w_tmp)
  endif

  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
       drho,xeps,xtemp,xye,dpress,keyerr,anyerr)

  call EOS_Omni_EpsFromPress(handle,keytemp,GRHydro_eos_rf_prec,n,&
       drho,xeps,xtemp,xye,dpress,deps,keyerr,anyerr)

  vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz
  vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
  vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz

  sqrtdet = sqrt(det)
  ddens = sqrtdet * drho * w 
  dsx = sqrtdet * (drho*(1+deps)+dpress)*w*w * vlowx
  dsy = sqrtdet * (drho*(1+deps)+dpress)*w*w * vlowy
  dsz = sqrtdet * (drho*(1+deps)+dpress)*w*w * vlowz
  dtau = sqrtdet * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens 

end subroutine prim2conpolytype


 /*@@
   @routine    Primitive2ConservativePolyCells
   @date       Sun Mar 10 21:16:20 2002
   @author     
   @desc 
   Wrapper function that converts primitive to conservative at the 
     cell centres. 
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/


subroutine Primitive2ConservativePolyCells(CCTK_ARGUMENTS)

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  
  CCTK_INT :: i, j, k
  CCTK_REAL :: det
  
  !$OMP PARALLEL DO PRIVATE(i, j, det)
     do k = 1,cctk_lsh(3)
        do j = 1,cctk_lsh(2)
           do i = 1,cctk_lsh(1)
        
        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 prim2conpolytype(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
             gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
             det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
             tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
             eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))

      end do
    end do
  end do
  !$OMP END PARALLEL DO

  if(evolve_Y_e.ne.0) then
     !$OMP PARALLEL DO PRIVATE(i, j)
     do k = 1,cctk_lsh(3)
        do j = 1,cctk_lsh(2)
           do i = 1,cctk_lsh(1)
              Y_e_con(i,j,k) = Y_e(i,j,k) * dens(i,j,k)
           enddo
        enddo
     enddo
     !$OMP END PARALLEL DO
  endif

end subroutine Primitive2ConservativePolyCells

 /*@@
   @routine    Prim2ConservativeTracer
   @date       Mon Mar  8 13:32:32 2004
   @author     Ian Hawke
   @desc 
   Gets the conserved tracer variable from the primitive.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS)

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  
  integer :: i, j, k
  CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
       gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
  
  do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
    do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
      do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
        
        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)
        cons_tracerplus(i,j,k,:) = tracerplus(i,j,k,:) * &
             sqrt(avg_detr) * rhoplus(i,j,k) / &
             sqrt(1.d0 - &
                   (gxxr * velxplus(i,j,k)**2 + &
                    gyyr * velyplus(i,j,k)**2 + &
                    gzzr * velzplus(i,j,k)**2 + &
                    2.d0 * (gxyr * velxplus(i,j,k) * velyplus(i,j,k) + &
                            gxzr * velxplus(i,j,k) * velzplus(i,j,k) + &
                            gyzr * velyplus(i,j,k) * velzplus(i,j,k) ) ) )
        cons_tracerminus(i,j,k,:) = tracerminus(i,j,k,:) * &
             sqrt(avg_detl) * rhominus(i,j,k) / &
             sqrt(1.d0 - &
                   (gxxl * velxminus(i,j,k)**2 + &
                    gyyl * velyminus(i,j,k)**2 + &
                    gzzl * velzminus(i,j,k)**2 + &
                    2.d0 * (gxyl * velxminus(i,j,k) * velyminus(i,j,k) + &
                            gxzl * velxminus(i,j,k) * velzminus(i,j,k) + &
                            gyzl * velyminus(i,j,k) * velzminus(i,j,k) ) ) )

      end do
    end do
  end do

end subroutine Prim2ConservativeTracer

 /*@@
   @routine   primitive2conservativegeneral
   @date       Thu Jan 11 11:03:32 2002
   @author     Pedro Montero, Ian Hawke
   @desc 
   Converts primitive to conserved variables for the boundary extended data.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine primitive2conservativegeneral(CCTK_ARGUMENTS)

  USE GRHydro_Scalars

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_FUNCTIONS
  
  CCTK_INT :: i, j, k, ierr
  CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl, &
       gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr, &
       vlowx, vlowy, vlowz
  
  ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConBndCallPlus)
  ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConBndCallMinus)

  !$OMP PARALLEL DO PRIVATE(i, j, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
  !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr)
  do k = GRHydro_stencil , cctk_lsh(3) - GRHydro_stencil + 1
    do j = GRHydro_stencil , cctk_lsh(2) - GRHydro_stencil + 1
      do i = GRHydro_stencil , cctk_lsh(1) - GRHydro_stencil + 1
        
        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 prim2conpolytype(-1, gxxl,gxyl,gxzl,&
             gyyl,gyzl,gzzl, &
             avg_detl,densminus(i,j,k),sxminus(i,j,k),&
             syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),rhominus(i,j,k), &
             velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),&
             epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k))

        call prim2conpolytype(-1, gxxr,gxyr,gxzr,&
             gyyr,gyzr,gzzr, &
             avg_detr, densplus(i,j,k),sxplus(i,j,k),&
             syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),&
             rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
             velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
             w_lorentzplus(i,j,k))
      end do
    end do
  end do
  !$OMP END PARALLEL DO

  ! Note on Y_e: We use Y_e_plus and Y_e_minus directly
  ! in the Riemann solver. That's why it is not necessary
  ! to do a prim2con for Y_e

end subroutine primitive2conservativegeneral

 /*@@
   @routine   Primitive2ConservativeCellsGeneral
   @date       Thu Jan 11 11:03:32 2002
   @author     Pedro Montero, Ian Hawke
   @desc 
   Converts primitive to conserved variables for the boundary extended data.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine Primitive2ConservativeCellsGeneral(CCTK_ARGUMENTS)

  USE GRHydro_Scalars

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_FUNCTIONS
  
  CCTK_INT :: i, j, k, ierr
  CCTK_REAL :: det,gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt, &
       vlowx, vlowy, vlowz
  
  ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConCellsCall)

  do k = 1,cctk_lsh(3)
    do j = 1,cctk_lsh(2)
      do i = 1,cctk_lsh(1)
        
        gxxpt = gxx(i,j,k)
        gxypt = gxy(i,j,k)
        gxzpt = gxz(i,j,k)
        gyypt = gyy(i,j,k)
        gyzpt = gyz(i,j,k)
        gzzpt = gzz(i,j,k)

        det = SPATIAL_DETERMINANT(gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt)

        w_lorentz(i,j,k) = 1.d0 / sqrt(1.d0 - &
             (gxxpt * velx(i,j,k)**2 + &
              gyypt * vely(i,j,k)**2 + &
              gzzpt * velz(i,j,k)**2 + &
              2.d0 * gxypt * velx(i,j,k) * vely(i,j,k) + &
              2.d0 * gxzpt * velx(i,j,k) * velz(i,j,k) + &
              2.d0 * gyzpt * vely(i,j,k) * velz(i,j,k) ))
        vlowx = gxxpt * velx(i,j,k) + &
                gxypt * vely(i,j,k) + &
                gxzpt * velz(i,j,k)
        vlowy = gxypt * velx(i,j,k) + &
                gyypt * vely(i,j,k) + &
                gyzpt * velz(i,j,k)
        vlowz = gxzpt * velx(i,j,k) + &
                gyzpt * vely(i,j,k) + &
                gzzpt * velz(i,j,k)
        dens(i,j,k) = sqrt(det) * rho(i,j,k) * &
             w_lorentz(i,j,k)
        sx(i,j,k) =  sqrt(det) * &
             ( rho(i,j,k) * &
               (1.d0 + eps(i,j,k)) + press(i,j,k) ) * &
               w_lorentz(i,j,k)**2 * vlowx
        sy(i,j,k) =  sqrt(det) * &
             ( rho(i,j,k) * &
               (1.d0 + eps(i,j,k)) + press(i,j,k) ) * &
               w_lorentz(i,j,k)**2 * vlowy
        sz(i,j,k) =  sqrt(det) * &
             ( rho(i,j,k) * &
               (1.d0 + eps(i,j,k)) + press(i,j,k) ) * &
               w_lorentz(i,j,k)**2 * vlowz
        tau(i,j,k) = sqrt(det) * &
             ( ( rho(i,j,k) * &
                 (1.d0 + eps(i,j,k)) + press(i,j,k) ) * &
                 w_lorentz(i,j,k)**2 - press(i,j,k) ) - &
             dens(i,j,k)

      end do
    end do
  end do

  if(evolve_Y_e.ne.0) then
     !$OMP PARALLEL DO PRIVATE(i, j)
     do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
        do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
           do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
              Y_e_con(i,j,k) = Y_e(i,j,k) * dens(i,j,k)
           enddo
        enddo
     enddo
     !$OMP END PARALLEL DO
  endif


end subroutine Primitive2ConservativeCellsGeneral