/*@@
   @file      GRHydro_ReconstructPoly.F90
   @date      Sat Jan 26 02:13:25 2002
   @author    
   @desc 
   Wrapper routine to perform the reconstruction for polytropes.
   @enddesc 
 @@*/

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

#include "SpaceMask.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)
#define Bvecx(i,j,k) Bvec(i,j,k,1)
#define Bvecy(i,j,k) Bvec(i,j,k,2)
#define Bvecz(i,j,k) Bvec(i,j,k,3)
#define Bconsx(i,j,k) Bcons(i,j,k,1)
#define Bconsy(i,j,k) Bcons(i,j,k,2)
#define Bconsz(i,j,k) Bcons(i,j,k,3)


 /*@@
   @routine    ReconstructionPolytype
   @date       Tue Mar 19 11:36:55 2002
   @author     Ian Hawke  
   @desc 
   If using a polytropic type EOS, we do not have to reconstruct all the 
   variables.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
  
  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_FUNCTIONS

  integer :: nx, ny, nz, i, j, k, itracer

  logical, dimension(:,:,:), allocatable :: trivial_rp
!!$  logical, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: trivial_rp

  CCTK_INT :: type_bitsx, trivialx, not_trivialx, &
       &type_bitsy, trivialy, not_trivialy, &
       &type_bitsz, trivialz, not_trivialz

  CCTK_REAL, dimension(:,:,:),allocatable :: &
       &psi4, lbetax, lbetay, lbetaz, dum, dump, dumm
!!$  CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: &
!!$       &psi4, lbetax, lbetay, lbetaz

  CCTK_INT :: ierr

  CCTK_REAL :: local_min_tracer

  CCTK_INT :: GRHydro_UseGeneralCoordinates

  allocate(trivial_rp(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr)
  if (ierr .ne. 0) then
    call CCTK_WARN(0, "Allocation problems with trivial_rp")
  end if
  
  allocate(psi4(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
       lbetax(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
       lbetay(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
       lbetaz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
       dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
       dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
       dumm(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr)

  if (ierr .ne. 0) then
    call CCTK_WARN(0, "Allocation problems with lbeta")
  end if
  
  if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
    call CCTK_WARN(0,"MP not yet supported in GRHydro_ReconstructionPolytype")
  end if

  call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX")
  call SpaceMask_GetStateBits(trivialx, "Hydro_RiemannProblemX", &
       &"trivial")
  call SpaceMask_GetStateBits(not_trivialx, "Hydro_RiemannProblemX", &
       &"not_trivial")
  call SpaceMask_GetTypeBits(type_bitsy, "Hydro_RiemannProblemY")
  call SpaceMask_GetStateBits(trivialy, "Hydro_RiemannProblemY", &
       &"trivial")
  call SpaceMask_GetStateBits(not_trivialy, "Hydro_RiemannProblemY", &
       &"not_trivial")
  call SpaceMask_GetTypeBits(type_bitsz, "Hydro_RiemannProblemZ")
  call SpaceMask_GetStateBits(trivialz, "Hydro_RiemannProblemZ", &
       &"trivial")
  call SpaceMask_GetStateBits(not_trivialz, "Hydro_RiemannProblemZ", &
       &"not_trivial")

  nx = cctk_lsh(1)
  ny = cctk_lsh(2)
  nz = cctk_lsh(3)

  !$OMP PARALLEL
  !$OMP WORKSHARE
  trivial_rp = .false.
  !$OMP END WORKSHARE NOWAIT

!!$  Currently only option is reconstruction on primitive variables.
!!$  Should change this.

  !$OMP WORKSHARE
  psi4 = 1.0d0
  !$OMP END WORKSHARE NOWAIT

  !$OMP WORKSHARE
  lbetax = betax
  lbetay = betay
  lbetaz = betaz
  !$OMP END WORKSHARE NOWAIT

!!$ Initialize variables that store reconstructed quantities

  !$OMP WORKSHARE
  rhoplus = 0.0d0
  rhominus = 0.0d0
  epsplus = 0.0d0
  epsminus = 0.0d0
  velxplus = 0.0d0
  velxminus = 0.0d0
  velyplus = 0.0d0
  velyminus = 0.0d0
  velzplus = 0.0d0
  velzminus = 0.0d0
  !$OMP END WORKSHARE NOWAIT
  if(evolve_mhd.ne.0) then
     !$OMP WORKSHARE
     Bvecxplus = 0.0d0
     Bvecxminus = 0.0d0
     Bvecyplus = 0.0d0
     Bvecyminus = 0.0d0
     Bveczplus = 0.0d0
     Bveczminus = 0.0d0
     !$OMP END WORKSHARE NOWAIT
     if(clean_divergence.ne.0) then
        !$OMP WORKSHARE
        psidcplus = 0.0d0
        psidcminus=0.0d0
        !$OMP END WORKSHARE NOWAIT
     endif
  endif

  if (evolve_tracer .ne. 0) then
     !$OMP WORKSHARE
     tracerplus = 0.0d0
     tracerminus = 0.0d0
     !$OMP END WORKSHARE NOWAIT
  endif

  if (evolve_Y_e .ne. 0) then
     !$OMP WORKSHARE
     Y_e_plus = 0.0d0
     Y_e_minus = 0.0d0
     !$OMP END WORKSHARE NOWAIT
  endif
  !$OMP END PARALLEL

  if (CCTK_EQUALS(recon_method,"tvd")) then

     if (evolve_tracer .ne. 0) then
        do itracer=1,number_of_tracers
           call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
                tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), &
                tracerminus(:,:,:,itracer), &
                trivial_rp, hydro_excision_mask)
        enddo
     end if
     
     if (evolve_Y_e .ne. 0) then
        call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
             Y_e(:,:,:), Y_e_plus(:,:,:), &
             Y_e_minus(:,:,:), &
             trivial_rp, hydro_excision_mask)
     endif

    if (CCTK_EQUALS(recon_vars,"primitive")) then
      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
           rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask)
      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
           vel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask)
      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
           vel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask)
      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
           vel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask)
      if(evolve_mhd.ne.0) then
         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
              Bvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask)
         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
              Bvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask)
         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
              Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask)
      endif
    else if (CCTK_EQUALS(recon_vars,"conservative")) then
      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
           dens, densplus, densminus, trivial_rp, hydro_excision_mask)
      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
           scon(:,:,:,1), sxplus, sxminus, trivial_rp, hydro_excision_mask)
      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
           scon(:,:,:,2), syplus, syminus, trivial_rp, hydro_excision_mask)
      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
           scon(:,:,:,3), szplus, szminus, trivial_rp, hydro_excision_mask)
      if(evolve_mhd.ne.0) then
         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
              Bcons(:,:,:,1), Bconsxplus, Bconsxminus, trivial_rp, hydro_excision_mask)
         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
              Bcons(:,:,:,2), Bconsyplus, Bconsyminus, trivial_rp, hydro_excision_mask)
         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
              Bcons(:,:,:,3), Bconszplus, Bconszminus, trivial_rp, hydro_excision_mask)
      endif
    else
      call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
    end if

    if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
            psidc, psidcplus, psidcminus, trivial_rp, hydro_excision_mask)
    endif

    !$OMP PARALLEL DO PRIVATE(i, j, k)
    do k = 1, nz
      do j = 1, ny
        do i = 1, nx
          if (trivial_rp(i,j,k)) then
            if (flux_direction == 1) then
              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx)
            else if (flux_direction == 2) then
              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsy, trivialy)
            else if (flux_direction == 3) then
              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsz, trivialz)
            end if
          else
            if (flux_direction == 1) then
              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx)
            else if (flux_direction == 2) then
              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsy, not_trivialy)
            else if (flux_direction == 3) then
              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsz, not_trivialz)
            end if
          end if
        end do
      end do
    end do
    !$OMP END PARALLEL DO

  else if (CCTK_EQUALS(recon_method,"ppm")) then

    if (flux_direction == 1) then
      if(evolve_mhd.ne.0) then
        if(clean_divergence.ne.0) then
          !$OMP PARALLEL DO PRIVATE(i, j, k)
          do k = GRHydro_stencil, nz - GRHydro_stencil + 1
            do j = GRHydro_stencil, ny - GRHydro_stencil + 1
              call SimplePPM_1dM(GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),&
                   rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),&
                   Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),psidc(:,j,k),eps(:,j,k),press(:,j,k),&
                   rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),&
                   Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),psidcminus(:,j,k),epsminus(:,j,k),&
                   rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),&
                   Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),psidcplus(:,j,k),epsplus(:,j,k),&
                   1,trivial_rp(:,j,k), hydro_excision_mask(:,j,k),&
                   gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), &
                   gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),&
                   w_lorentz(:,j,k), &
                   1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, &
                   GRHydro_mppm_eigenvalue_x_right, &
                   GRHydro_mppm_xwind)
              do i = 1, nx
                if (trivial_rp(i,j,k)) then
                  SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx)
                else
                  SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx)
                end if
              end do
            end do
          end do
          !$OMP END PARALLEL DO
        else  !clean_divergence
          !$OMP PARALLEL DO PRIVATE(i, j, k)
          do k = GRHydro_stencil, nz - GRHydro_stencil + 1
            do j = GRHydro_stencil, ny - GRHydro_stencil + 1
              call SimplePPM_1dM(GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),&
                   rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),&
                   Bvecx(:,j,k),Bvecy(:,j,k),Bvecz(:,j,k),dum(:,j,k),eps(:,j,k),press(:,j,k),&
                   rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),&
                   Bvecxminus(:,j,k),Bvecyminus(:,j,k),Bveczminus(:,j,k),dumm(:,j,k),epsminus(:,j,k),&
                   rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),&
                   Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),dump(:,j,k),epsplus(:,j,k),&
                   0,trivial_rp(:,j,k), hydro_excision_mask(:,j,k),&
                   gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), &
                   gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),&
                   w_lorentz(:,j,k), &
                   1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, &
                   GRHydro_mppm_eigenvalue_x_right, &
                   GRHydro_mppm_xwind)
              do i = 1, nx
                if (trivial_rp(i,j,k)) then
                  SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx)
                else
                  SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx)
                end if
              end do
            end do
          end do
          !$OMP END PARALLEL DO
        endif  !clean_divergence
      else  !evolve_mhd
        !$OMP PARALLEL DO PRIVATE(i, j, k)
        do k = GRHydro_stencil, nz - GRHydro_stencil + 1
          do j = GRHydro_stencil, ny - GRHydro_stencil + 1
            call SimplePPM_1d(GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),&
                 rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),eps(:,j,k),&
                 press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),&
                 velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),&
                 velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),&
                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k),&
                 gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), &
                 gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),&
                 w_lorentz(:,j,k), &
                 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, &
                 GRHydro_mppm_eigenvalue_x_right, &
                 GRHydro_mppm_xwind)
            do i = 1, nx
              if (trivial_rp(i,j,k)) then
                SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx)
              else
                SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx)
              end if
            end do
          end do
        end do
        !$OMP END PARALLEL DO
      endif  !evolve_mhd
      if(evolve_tracer.ne.0) then
         !$OMP PARALLEL DO PRIVATE(j, k)
         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
            do j = GRHydro_stencil, ny - GRHydro_stencil + 1
               call SimplePPM_tracer_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), &
                    velx(:,j,k),vely(:,j,k),velz(:,j,k), &
                    tracer(:,j,k,:),tracerminus(:,j,k,:),tracerplus(:,j,k,:), &
                    press(:,j,k))
            end do
         end do
         !$OMP END PARALLEL DO
      end if
      if(evolve_Y_e.ne.0) then
         !$OMP PARALLEL DO PRIVATE(j, k)
         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
            do j = GRHydro_stencil, ny - GRHydro_stencil + 1
               call SimplePPM_tracer_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), &
                    velx(:,j,k),vely(:,j,k),velz(:,j,k), &
                    Y_e(:,j,k),Y_e_minus(:,j,k),Y_e_plus(:,j,k), &
                    press(:,j,k))
            end do
         end do
         !$OMP END PARALLEL DO
      end if

    else if (flux_direction == 2) then
      if(evolve_mhd.ne.0) then
        if(clean_divergence.ne.0) then
          !$OMP PARALLEL DO PRIVATE(i, j, k)
          do k = GRHydro_stencil, nz - GRHydro_stencil + 1
            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
              call SimplePPM_1dM(GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),&
                   rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),&
                   Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),psidc(j,:,k),eps(j,:,k),press(j,:,k),&
                   rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),&
                   Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),psidcminus(j,:,k),epsminus(j,:,k),&
                   rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),&
                   Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),psidcplus(j,:,k),epsplus(j,:,k),&
                   1,trivial_rp(j,:,k), hydro_excision_mask(j,:,k),&
                   gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), &
                   gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),&
                   w_lorentz(j,:,k), &
                   2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, &
                   GRHydro_mppm_eigenvalue_y_right, &
                   GRHydro_mppm_xwind)
              do i = 1, ny
                if (trivial_rp(j,i,k)) then
                  SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy)
                else
                  SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy)
                end if
              end do
            end do
          end do
          !$OMP END PARALLEL DO
        else  !clean_divergence
          !$OMP PARALLEL DO PRIVATE(i, j, k)
          do k = GRHydro_stencil, nz - GRHydro_stencil + 1
            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
              call SimplePPM_1dM(GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),&
                   rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),&
                   Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),dum(j,:,k),eps(j,:,k),press(j,:,k),&
                   rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),&
                   Bvecyminus(j,:,k),Bveczminus(j,:,k),Bvecxminus(j,:,k),dumm(j,:,k),epsminus(j,:,k),&
                   rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),&
                   Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),dump(j,:,k),epsplus(j,:,k),&
                   0,trivial_rp(j,:,k), hydro_excision_mask(j,:,k),&
                   gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), &
                   gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),&
                   w_lorentz(j,:,k), &
                   2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, &
                   GRHydro_mppm_eigenvalue_y_right, &
                   GRHydro_mppm_xwind)
              do i = 1, ny
                if (trivial_rp(j,i,k)) then
                  SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy)
                else
                  SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy)
                end if
              end do
            end do
          end do
          !$OMP END PARALLEL DO
        endif  !clean_divergence
      else  !evolve_mhd
        !$OMP PARALLEL DO PRIVATE(i, j, k)
        do k = GRHydro_stencil, nz - GRHydro_stencil + 1
          do j = GRHydro_stencil, nx - GRHydro_stencil + 1
            call SimplePPM_1d(GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),&
                 rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),&
                 press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),&
                 velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),&
                 velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),&
                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k),&
                 gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), &
                 gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),&
                 w_lorentz(j,:,k), &
                 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, &
                 GRHydro_mppm_eigenvalue_y_right, &
                 GRHydro_mppm_xwind)
            do i = 1, ny
              if (trivial_rp(j,i,k)) then
                SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy)
              else
                SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy)
              end if
            end do
          end do
        end do
        !$OMP END PARALLEL DO
      endif  !evolve_mhd
      if(evolve_tracer.ne.0) then
         !$OMP PARALLEL DO PRIVATE(j, k)
         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
              call SimplePPM_tracer_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), &
                   vely(j,:,k),velz(j,:,k),velx(j,:,k), &
                   tracer(j,:,k,:),tracerminus(j,:,k,:),tracerplus(j,:,k,:), &
                   press(j,:,k))
            end do
         end do
         !$OMP END PARALLEL DO
      end if
      if(evolve_Y_e.ne.0) then
         !$OMP PARALLEL DO PRIVATE(j, k)
         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
               call SimplePPM_tracer_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), &
                    velx(j,:,k),vely(j,:,k),velz(j,:,k), &
                    Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), &
                    press(j,:,k))
            end do
         end do
         !$OMP END PARALLEL DO
      end if
      
    else if (flux_direction == 3) then
      if(evolve_mhd.ne.0) then
        if(clean_divergence.ne.0) then
          !$OMP PARALLEL DO PRIVATE(i, j, k)
          do k = GRHydro_stencil, ny - GRHydro_stencil + 1
            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
              call SimplePPM_1dM(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),&
                   rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),&
                   Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),psidc(j,k,:),eps(j,k,:),press(j,k,:),&
                   rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),&
                   Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),psidcminus(j,k,:),epsminus(j,k,:),&
                   rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),&
                   Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),psidcplus(j,k,:),epsplus(j,k,:),&
                   1,trivial_rp(j,k,:), hydro_excision_mask(j,k,:),&
                   gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), &
                   gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),&
                   w_lorentz(j,k,:), &
                   3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, &
                   GRHydro_mppm_eigenvalue_z_right, &
                   GRHydro_mppm_xwind)
              do i = 1, nz
                if (trivial_rp(j,k,i)) then
                  SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz)
                else
                  SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz)
                end if
              end do
            end do
          end do
          !$OMP END PARALLEL DO
        else  !clean_divergence
          !$OMP PARALLEL DO PRIVATE(i, j, k)
          do k = GRHydro_stencil, ny - GRHydro_stencil + 1
            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
              call SimplePPM_1dM(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),&
                   rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),&
                   Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),dum(:,j,k),eps(j,k,:),press(j,k,:),&
                   rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),&
                   Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),dumm(j,k,:),epsminus(j,k,:),&
                   rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),&
                   Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),dump(j,k,:),epsplus(j,k,:),&
                   0,trivial_rp(j,k,:), hydro_excision_mask(j,k,:),&
                   gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), &
                   gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),&
                   w_lorentz(j,k,:), &
                   3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, &
                   GRHydro_mppm_eigenvalue_z_right, &
                   GRHydro_mppm_xwind)
              do i = 1, nz
                if (trivial_rp(j,k,i)) then
                  SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz)
                else
                  SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz)
                end if
              end do
            end do
          end do
          !$OMP END PARALLEL DO
        endif  !clean_divergence
      else  !evolve_mhd
        !$OMP PARALLEL DO PRIVATE(i, j, k)
        do k = GRHydro_stencil, ny - GRHydro_stencil + 1
          do j = GRHydro_stencil, nx - GRHydro_stencil + 1
            call SimplePPM_1d(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),&
                 rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),eps(j,k,:),&
                 press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),&
                 velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),&
                 velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),&
                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:),&
                 gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), &
                 gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),&
                 w_lorentz(j,k,:), &
                 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, &
                 GRHydro_mppm_eigenvalue_z_right, &
                 GRHydro_mppm_xwind)
            do i = 1, nz
              if (trivial_rp(j,k,i)) then
                SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz)
              else
                SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz)
              end if
            end do
          end do
        end do
        !$OMP END PARALLEL DO
      endif  !evolve_mhd
      if(evolve_tracer.ne.0) then
         !$OMP PARALLEL DO PRIVATE(j, k)
         do k = GRHydro_stencil, ny - GRHydro_stencil + 1
            do j = GRHydro_stencil, nx - GRHydro_stencil + 1

               call SimplePPM_tracer_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), &
                    velz(j,k,:),velx(j,k,:),vely(j,k,:), &
                    tracer(j,k,:,:),tracerminus(j,k,:,:),tracerplus(j,k,:,:), &
                    press(j,k,:))
            end do
         end do
         !$OMP END PARALLEL DO
      end if
      if(evolve_Y_e.ne.0) then
         !$OMP PARALLEL DO PRIVATE(j, k)
         do k = GRHydro_stencil, ny - GRHydro_stencil + 1
            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
               call SimplePPM_tracer_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), &
                    velx(j,k,:),vely(j,k,:),velz(j,k,:), &
                    Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(j,k,:), &
                    press(j,k,:))
            end do
         end do
         !$OMP END PARALLEL DO
      end if
    else
      call CCTK_WARN(0, "Flux direction not x,y,z")
    end if
  else if (CCTK_EQUALS(recon_method,"eno")) then

     if (evolve_tracer .ne. 0) then
        do itracer=1,number_of_tracers
           call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
                tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), &
                tracerminus(:,:,:,itracer), trivial_rp, &
                hydro_excision_mask)
        enddo
     end if

     if (evolve_Y_e .ne. 0) then
        call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
             Y_e(:,:,:), Y_e_plus(:,:,:), &
             Y_e_minus(:,:,:), &
             trivial_rp, hydro_excision_mask)
     endif

    if (flux_direction == 1) then
      !$OMP PARALLEL DO PRIVATE(i, j, k)
      do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1
        do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1
          if (CCTK_EQUALS(recon_vars,"primitive")) then
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                 rho(:,j,k),rhominus(:,j,k),rhoplus(:,j,k),&
                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                 velx(:,j,k),velxminus(:,j,k),velxplus(:,j,k),&
                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                 vely(:,j,k),velyminus(:,j,k),velyplus(:,j,k),&
                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                 velz(:,j,k),velzminus(:,j,k),velzplus(:,j,k),&
                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
            if(evolve_mhd.ne.0) then
               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                    Bvecx(:,j,k),Bvecxminus(:,j,k),Bvecxplus(:,j,k),&
                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                    Bvecy(:,j,k),Bvecyminus(:,j,k),Bvecyplus(:,j,k),&
                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                    Bvecz(:,j,k),Bveczminus(:,j,k),Bveczplus(:,j,k),&
                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
            endif
          else if (CCTK_EQUALS(recon_vars,"conservative")) then
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                 dens(:,j,k),densminus(:,j,k),densplus(:,j,k),&
                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                 sx(:,j,k),sxminus(:,j,k),sxplus(:,j,k),&
                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                 sy(:,j,k),syminus(:,j,k),syplus(:,j,k),&
                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                 sz(:,j,k),szminus(:,j,k),szplus(:,j,k),&
                 trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
            if(evolve_mhd.ne.0) then
                 call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                      Bconsx(:,j,k),Bconsxminus(:,j,k),Bconsxplus(:,j,k),&
                      trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
                 call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                      Bconsy(:,j,k),Bconsyminus(:,j,k),Bconsyplus(:,j,k),&
                      trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
                 call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                      Bconsz(:,j,k),Bconszminus(:,j,k),Bconszplus(:,j,k),&
                      trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
            endif
          else
            !$OMP CRITICAL
            call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
            !$OMP END CRITICAL
          end if

          if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  psidc(:,j,k),psidcminus(:,j,k),psidcplus(:,j,k),&
                  trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
          endif

          do i = 1, cctk_lsh(1)
            if (trivial_rp(i,j,k)) then
              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx)
            else
              SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx)
            end if
          end do
        end do
       end do
        !$OMP END PARALLEL DO
     else if (flux_direction == 2) then
        !$OMP PARALLEL DO PRIVATE(i, j, k)
        do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1
           do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1
              if (CCTK_EQUALS(recon_vars,"primitive")) then
                 call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                      rho(j,:,k),rhominus(j,:,k),rhoplus(j,:,k),&
                      trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                 call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                      velx(j,:,k),velxminus(j,:,k),velxplus(j,:,k),&
                      trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                 call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                      vely(j,:,k),velyminus(j,:,k),velyplus(j,:,k),&
                      trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                 call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                      velz(j,:,k),velzminus(j,:,k),velzplus(j,:,k),&
                      trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                 if(evolve_mhd.ne.0) then
                    call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                         Bvecx(j,:,k),Bvecxminus(j,:,k),Bvecxplus(j,:,k),&
                         trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                    call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                         Bvecy(j,:,k),Bvecyminus(j,:,k),Bvecyplus(j,:,k),&
                         trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                    call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                         Bvecz(j,:,k),Bveczminus(j,:,k),Bveczplus(j,:,k),&
                         trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                 end if

              else if (CCTK_EQUALS(recon_vars,"conservative")) then
                 call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                      dens(j,:,k),densminus(j,:,k),densplus(j,:,k),&
                      trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                 call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                      sx(j,:,k),sxminus(j,:,k),sxplus(j,:,k),&
                      trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                 call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                      sy(j,:,k),syminus(j,:,k),syplus(j,:,k),&
                      trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                 call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                      sz(j,:,k),szminus(j,:,k),szplus(j,:,k),&
                      trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                 if(evolve_mhd.ne.0) then
                    call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                         Bconsx(j,:,k),Bconsxminus(j,:,k),Bconsxplus(j,:,k),&
                         trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                    call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                         Bconsy(j,:,k),Bconsyminus(j,:,k),Bconsyplus(j,:,k),&
                         trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                    call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                         Bconsz(j,:,k),Bconszminus(j,:,k),Bconszplus(j,:,k),&
                         trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
                 end if
              else
                 !$OMP CRITICAL
                 call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
                 !$OMP END CRITICAL
              end if

              if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
                 call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                      psidc(j,:,k),psidcminus(j,:,k),psidcplus(j,:,k),&
                      trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
              endif

              do i = 1, cctk_lsh(2)
                 if (trivial_rp(j,i,k)) then
                    SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy)
                 else
                    SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy)
                 end if
              end do
           end do
        end do
      !$OMP END PARALLEL DO
    else if (flux_direction == 3) then
      !$OMP PARALLEL DO PRIVATE(i, j, k)
      do k = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1
        do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1
          if (CCTK_EQUALS(recon_vars,"primitive")) then
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                 rho(j,k,:),rhominus(j,k,:),rhoplus(j,k,:),&
                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                 velx(j,k,:),velxminus(j,k,:),velxplus(j,k,:),&
                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                 vely(j,k,:),velyminus(j,k,:),velyplus(j,k,:),&
                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                 velz(j,k,:),velzminus(j,k,:),velzplus(j,k,:),&
                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
            if(evolve_mhd.ne.0) then
               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                    Bvecx(j,k,:),Bvecxminus(j,k,:),Bvecxplus(j,k,:),&
                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                    Bvecy(j,k,:),Bvecyminus(j,k,:),Bvecyplus(j,k,:),&
                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                    Bvecz(j,k,:),Bveczminus(j,k,:),Bveczplus(j,k,:),&
                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
            endif
          else if (CCTK_EQUALS(recon_vars,"conservative")) then
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                 dens(j,k,:),densminus(j,k,:),densplus(j,k,:),&
                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                 sx(j,k,:),sxminus(j,k,:),sxplus(j,k,:),&
                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                 sy(j,k,:),syminus(j,k,:),syplus(j,k,:),&
                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                 sz(j,k,:),szminus(j,k,:),szplus(j,k,:),&
                 trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
            if(evolve_mhd.ne.0) then
               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                    Bconsx(j,k,:),Bconsxminus(j,k,:),Bconsxplus(j,k,:),&
                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                    Bconsy(j,k,:),Bconsyminus(j,k,:),Bconsyplus(j,k,:),&
                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                    Bconsz(j,k,:),Bconszminus(j,k,:),Bconszplus(j,k,:),&
                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
            endif
          else
            !$OMP CRITICAL
            call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
            !$OMP END CRITICAL
          end if

          if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  psidc(j,k,:),psidcminus(j,k,:),psidcplus(j,k,:),&
                  trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
          endif

          do i = 1, cctk_lsh(3)
            if (trivial_rp(j,k,i)) then
              SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz)
            else
              SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz)
            end if
          end do
        end do
      end do
      !$OMP END PARALLEL DO
    else
      call CCTK_WARN(0, "Flux direction not x,y,z")
    end if
  else
    call CCTK_WARN(0, "Reconstruction method not recognized!")
  end if

  deallocate(trivial_rp)
  deallocate(psi4, lbetax, lbetay, lbetaz)
  deallocate(dum, dump, dumm)

  !$OMP PARALLEL WORKSHARE
  where ( (rhoplus < GRHydro_rho_min).or.(rhominus < GRHydro_rho_min).or.&
          (epsplus < 0.d0).or.(epsminus < 0.d0) )
    rhoplus = rho
    rhominus = rho
    velxplus = vel(:,:,:,1)
    velxminus = vel(:,:,:,1)
    velyplus = vel(:,:,:,2)
    velyminus = vel(:,:,:,2)
    velzplus = vel(:,:,:,3)
    velzminus = vel(:,:,:,3)
    epsplus = eps
    epsminus = eps
  end where
  !$OMP END PARALLEL WORKSHARE

  if (evolve_tracer .ne. 0) then
    if (use_min_tracer .ne. 0) then
      local_min_tracer = min_tracer
    else
      local_min_tracer = 0d0
    end if
   
    !$OMP PARALLEL WORKSHARE
    where( (tracerplus  .le. local_min_tracer).or.&
           (tracerminus .le. local_min_tracer) )
      tracerplus = tracer
      tracerminus = tracer
    end where
    !$OMP END PARALLEL WORKSHARE
    ! Call the conserved tracer routine in any case because (accord. to
    ! Christian Ott) this is the only way this works
    call Prim2ConservativeTracer(CCTK_PASS_FTOF)
  endif

  if (CCTK_EQUALS(recon_vars,"primitive").or.&
       CCTK_EQUALS(recon_method,"ppm")) then
     if(evolve_mhd.ne.0) then
        call Prim2ConservativePolytypeM(CCTK_PASS_FTOF)
     else
        call Prim2ConservativePolytype(CCTK_PASS_FTOF)
     endif
  else if (CCTK_EQUALS(recon_vars,"conservative")) then
     if(evolve_mhd.ne.0) then
        call Con2PrimBoundsPolytype(CCTK_PASS_FTOF)
     else
        call Con2PrimBoundsPolytype(CCTK_PASS_FTOF)
     endif
  else
     call CCTK_WARN(0,"Variable type to reconstruct not recognized.")
  end if
  
  return
  
end subroutine ReconstructionPolytype