/*@@
   @file      GRHydro_PPMReconstruct_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 PPM reconstructions.
   @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)


 /*@@
   @routine    GRHydro_PPMReconstruct_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 PPM reconstructions. 
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine GRHydro_PPMReconstruct_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 gaa, gab, gac, gbb, gbc, gcc
  TARGET gxx, gxy, gxz, gyy, gyz, gzz
  TARGET betaa, betab, betac
  TARGET betax, betay, betaz
  TARGET lvel, vel

  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_FUNCTIONS

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

  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 :: g11, g12, g13, g22, g23, g33
  CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3
  CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup

  if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
    g11 => gaa
    g12 => gab
    g13 => gac
    g22 => gbb
    g23 => gbc
    g33 => gcc
    beta1 => betaa
    beta2 => betab
    beta3 => betac
    vup => lvel
  else
    g11 => gxx
    g12 => gxy
    g13 => gxz
    g22 => gyy
    g23 => gyz
    g33 => gzz
    beta1 => betax
    beta2 => betay
    beta3 => betaz
    vup => vel
  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_tracer .ne. 0) then
              tracerplus(i,j,k,:) = 0.0d0
              tracerminus(i,j,k,:) = 0.0d0
           endif

           if (evolve_Y_e .ne. 0) then
              ! set this to the cell center values
              ! to make sure we have good Y_e even in
              ! the boundary region (for full GF EOS calls)
              Y_e_plus(i,j,k) = Y_e(i,j,k)
              Y_e_minus(i,j,k) = Y_e(i,j,k)
           endif

           if (evolve_temper .ne. 0) then
              ! set this to temperature (piecewise constant), because we won't
              ! set it if reconstruct_temper .eq. 0
              tempplus(i,j,k) = temperature(i,j,k)
              tempminus(i,j,k) = temperature(i,j,k)
           endif

        enddo
     enddo
  enddo
  !$OMP END PARALLEL DO

!!$ PPM starts:
    if (flux_direction == 1) 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_1d(GRHydro_eos_handle,0,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),&
                g11(:,j,k), g12(:,j,k), g13(:,j,k), g22(:,j,k), g23(:,j,k), &
                g33(:,j,k), beta1(:,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
      if(evolve_temper.ne.0.and.reconstruct_temper.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_temperature_1d(nx,CCTK_DELTA_SPACE(1),velx(:,j,k), & 
                    temperature(:,j,k),press(:,j,k),&
                    tempminus(:,j,k),tempplus(:,j,k), hydro_excision_mask)
            end do
         end do
         !$OMP END PARALLEL DO
      endif
      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_ye_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
        !$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,0,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),&
                 g22(j,:,k), g23(j,:,k), g12(j,:,k), g33(j,:,k), g13(j,:,k), &
                 g11(j,:,k), beta2(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
      if(evolve_temper.ne.0.and.reconstruct_temper.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_temperature_1d(ny,CCTK_DELTA_SPACE(2),vely(j,:,k), & 
                    temperature(j,:,k),press(j,:,k),&
                    tempminus(j,:,k),tempplus(j,:,k), hydro_excision_mask)
            end do
         end do
         !$OMP END PARALLEL DO
      endif
      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_ye_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), &
                    vely(j,:,k),velz(j,:,k),velx(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
        !$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,0,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,:),&
                 g33(j,k,:), g13(j,k,:), g23(j,k,:), g11(j,k,:), g12(j,k,:), &
                 g22(j,k,:), beta3(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
      if(evolve_temper.ne.0.and.reconstruct_temper.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_temperature_1d(nz,CCTK_DELTA_SPACE(3),velz(j,k,:), & 
                    temperature(j,k,:),press(j,k,:),&
                    tempminus(j,k,:),tempplus(j,k,:), hydro_excision_mask)
            end do
         end do
         !$OMP END PARALLEL DO
      endif
      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_ye_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), &
                    velz(j,k,:),velx(j,k,:),vely(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
!!$ PPM ends.

  deallocate(trivial_rp)

end subroutine GRHydro_PPMReconstruct_drv