/*@@
   @file      GRHydro_TVDReconstruct_drv.F90
   @date      Tue Jul 19 13:22:03 EDT 2011
   @author    Bruno C. Mundim, Joshua Faber, Christian D. Ott
   @desc 
   Driver routine to perform the TVD 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) vup(i,j,k,1)
#define vely(i,j,k) vup(i,j,k,2)
#define velz(i,j,k) vup(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) Bprim(i,j,k,1)
#define Bvecy(i,j,k) Bprim(i,j,k,2)
#define Bvecz(i,j,k) Bprim(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_TVDReconstruct_drv
   @date       Tue Jul 19 13:24:34 EDT 2011
   @author     Luca Baiotti, Ian Hawke, Bruno C. Mundim, Joshua Faber, Christian D. Ott
   @desc 
   A driver routine to do TVD reconstruction. Currently just does
   TVD on the primitive variables.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine GRHydro_TVDReconstruct_drv(CCTK_ARGUMENTS)
  
  implicit none
  
  ! save memory when MP is not used
  ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1
  TARGET vel, lvel
  TARGET Bvec, lBvec

  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_FUNCTIONS

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

  logical, dimension(:,:,:), allocatable :: trivial_rp

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

  CCTK_INT :: ierr

  ! save memory when MP is not used
  CCTK_INT :: GRHydro_UseGeneralCoordinates
  CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim

  if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
    vup => lvel
    Bprim => lBvec
  else
    vup => vel
    Bprim => Bvec
  end if

  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
  
  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)

!!$ Initialize variables that store reconstructed quantities:

  !$OMP PARALLEL DO PRIVATE(i,j,k)
  do k=1,cctk_lsh(3)
     do j=1,cctk_lsh(2)
        do i=1,cctk_lsh(1)
           trivial_rp(i,j,k) = .false.
           rhoplus(i,j,k) = 0.0d0
           rhominus(i,j,k)= 0.0d0
           epsplus(i,j,k) = 0.0d0
           epsminus(i,j,k) = 0.0d0
           velxplus(i,j,k) = 0.0d0
           velxminus(i,j,k) = 0.0d0
           velyplus(i,j,k) = 0.0d0
           velyminus(i,j,k) = 0.0d0
           velzplus(i,j,k) = 0.0d0
           velzminus(i,j,k) = 0.0d0

           if(evolve_mhd.ne.0) then
             Bvecxplus(i,j,k) = 0.0d0
             Bvecxminus(i,j,k) = 0.0d0
             Bvecyplus(i,j,k) = 0.0d0
             Bvecyminus(i,j,k) = 0.0d0
             Bveczplus(i,j,k) = 0.0d0
             Bveczminus(i,j,k) = 0.0d0
             if(clean_divergence.ne.0) then
                psidcplus(i,j,k) = 0.0d0
                psidcminus(i,j,k) = 0.0d0
             endif
           endif

           if (evolve_entropy .ne. 0) then
              entropyplus(i,j,k) = 0.0d0
              entropyminus(i,j,k) = 0.0d0
           endif

           if (evolve_tracer .ne. 0) then
              tracerplus(i,j,k,:) = 0.0d0
              tracerminus(i,j,k,:) = 0.0d0
           endif

           if (evolve_Y_e .ne. 0) then
              Y_e_plus(i,j,k) = 0.0d0
              Y_e_minus(i,j,k) = 0.0d0
           endif

        enddo
     enddo
  enddo
  !$OMP END PARALLEL DO

!!$ TVD 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 (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, &
           vup(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask)
      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
           vup(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask)
      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
           vup(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask)
      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
           eps, epsplus, epsminus, trivial_rp, hydro_excision_mask)
      if(evolve_mhd.ne.0) then
         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
              Bprim(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask)
         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
              Bprim(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask)
         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
              Bprim(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask)
      endif
      if (evolve_entropy .ne. 0) then
         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
              entropy, entropyplus, entropyminus, 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)
      call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
           tau, tauplus, tauminus, 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
      if (evolve_entropy .ne. 0) then
         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
              entropycons, entropyconsplus, entropyconsminus, 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
!!$ TVD ends.

  deallocate(trivial_rp)

end subroutine GRHydro_TVDReconstruct_drv