RIT GRMHD dev: split reconstruct routine.

Comments:

1) The reconstruct routine is the most expensive routine
in GRHydro and may require some serious thought in the
future regarding optimization.

2) This routine was merged by the end of the last year.
With too many loops, the compiler wasn't able to fully
optimize it (and it would warn about that).

3) With this split, the code became easier to read,
the compiler was able to fully optimize it, but we only
have a marginal gain in the routine performance (3.2%). 
Further optimization and clean-up is still desirable.



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

1 files changed, 403 insertions, 0 deletions

diff --git a/src/GRHydro_ENOReconstruct_drv.F90 b/src/GRHydro_ENOReconstruct_drv.F90
new file mode 100644
index 0000000..656fd6e
--- /dev/null
+++ b/src/GRHydro_ENOReconstruct_drv.F90
@@ -0,0 +1,403 @@
+ /*@@
+   @file      GRHydro_ENOReconstruct_drv.F90
+   @date      Tue Jul 19 13:22:03 EDT 2011
+   @author    Bruno C. Mundim, Joshua Faber
+   @desc 
+   Driver routine to perform the ENO reconstruction.
+   @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    GRHydro_ENOReconstruct_drv
+   @date       Tue Jul 19 13:24:34 EDT 2011
+   @author     Luca Baiotti, Ian Hawke, Bruno C. Mundim, Joshua Faber
+   @desc 
+   A driver routine to do ENO reconstruction. 
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_ENOReconstruct_drv(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
+
+  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
+  
+!!$ The dum variables are used as dummies if MHD on but divergence cleaning off
+  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
+  
+  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.d0
+  !$OMP END WORKSHARE NOWAIT
+  if (shift_state .ne. 0) then
+    !$OMP WORKSHARE
+    lbetax = betax
+    lbetay = betay
+    lbetaz = betaz
+    !$OMP END WORKSHARE NOWAIT
+  else
+    !$OMP WORKSHARE
+    lbetax = 0.d0
+    lbetay = 0.d0
+    lbetaz = 0.d0
+    !$OMP END WORKSHARE NOWAIT
+  end if
+  !$OMP END PARALLEL
+
+!!$ ENO starts:
+     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))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                 eps(:,j,k),epsminus(:,j,k),epsplus(:,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))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                 tau(:,j,k),tauminus(:,j,k),tauplus(:,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))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 eps(j,:,k),epsminus(j,:,k),epsplus(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))
+            endif
+         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))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 tau(j,:,k),tauminus(j,:,k),tauplus(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))
+            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(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,:))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                 eps(j,k,:),epsminus(j,k,:),epsplus(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,:))
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                 tau(j,k,:),tauminus(j,k,:),tauplus(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
+!!$ ENO ends.
+
+  deallocate(trivial_rp)
+  deallocate(psi4, lbetax, lbetay, lbetaz)
+  deallocate(dum,dump,dumm)
+
+end subroutine GRHydro_ENOReconstruct_drv
+