/*@@
   @file      GRHydro_UpdateMaskM.F90
   @date      Sep 2, 2010
   @author    
   @desc 
   Alter the update terms if inside the atmosphere or excision region
   @enddesc 
 @@*/

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

#include "GRHydro_Macros.h"
#include "SpaceMask.h"

!!$ We don't need to adapt GRHydroUpdateAtmosphereMask, GRHydro_SetupMask, or  
!!$ since we need to evolve Bvec in the atmosphere

!!$ In GRHydro_AtmosphereResetM, we just need to switch the P2C calls to MHD

 /*@@
   @routine    GRHydro_AtmosphereResetM
   @date       Sep 2, 2010
   @author     Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
   @desc 
   After MoL has evolved, if a point is supposed to be reset then do so.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine GRHydro_AtmosphereResetM(CCTK_ARGUMENTS)

  implicit none

  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  CCTK_INT :: i, j, k
  CCTK_REAL :: det, psi4pt

  do k = 1, cctk_lsh(3)
    do j = 1, cctk_lsh(2)
      do i = 1, cctk_lsh(1)
        
        if ( atmosphere_mask(i, j, k) .eq. 1) then

          rho(i,j,k) = GRHydro_rho_min
          vel(i,j,k,1) = 0.0d0
          vel(i,j,k,2) = 0.0d0
          vel(i,j,k,3) = 0.0d0
          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 prim2conpolytypeM(GRHydro_polytrope_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), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
               tau(i,j,k), Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),&
               rho(i,j,k), vel(i,j,k,1), vel(i,j,k,2), &
               vel(i,j,k,3), eps(i,j,k), press(i,j,k), &
               Bvec(i,j,k,1),Bvec(i,j,k,2),Bvec(i,j,k,3),w_lorentz(i,j,k))
          if (wk_atmosphere .eq. 0) then
            atmosphere_mask(i, j, k) = 0
          end if

        end if

      end do
    end do
  end do

!!$  call GRHydro_BoundariesM(CCTK_PASS_FTOF)
  
end subroutine GRHydro_AtmosphereResetM

subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS)

  implicit none

  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  CCTK_INT :: i, j, k
  CCTK_REAL :: det, psi4pt

  CCTK_INT :: eos_handle

! 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;xye=0.0d0;xeps=0.0d0;xtemp=0.0d0
! end EOS Omni vars

  eos_handle = GRHydro_polytrope_handle
  
  do k = 1, cctk_lsh(3)
    do j = 1, cctk_lsh(2)
      do i = 1, cctk_lsh(1)
        
        if (rho(i,j,k) .le. GRHydro_rho_min) then
          rho(i,j,k) = GRHydro_rho_min
          vel(i,j,k,1) = 0.0d0
          vel(i,j,k,2) = 0.0d0
          vel(i,j,k,3) = 0.0d0

          call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                 rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
          call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                 rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)

          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 prim2conpolytypeM(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), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
               tau(i,j,k), Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),&
               rho(i,j,k), vel(i,j,k,1), vel(i,j,k,2), &
               vel(i,j,k,3), eps(i,j,k), press(i,j,k), &
               Bvec(i,j,k,1),Bvec(i,j,k,2),Bvec(i,j,k,3),w_lorentz(i,j,k))
        end if
        if (timelevels .gt. 1) then
          if (rho_p(i,j,k) .le. GRHydro_rho_min) then
            rho_p(i,j,k) = GRHydro_rho_min
            vel_p(i,j,k,1) = 0.0d0
            vel_p(i,j,k,2) = 0.0d0
            vel_p(i,j,k,3) = 0.0d0

          call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                 rho_p(i,j,k),eps_p(i,j,k),xtemp,xye,press_p(i,j,k),keyerr,anyerr)
          call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                 rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr)

            det = SPATIAL_DETERMINANT(gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), \
                  gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k))
            call prim2conpolytypeM(eos_handle, &
                 gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), &
                 gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k), det, &
                 dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), &
                 tau_p(i,j,k), Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),&
                 rho_p(i,j,k), vel_p(i,j,k,1), vel_p(i,j,k,2), &
                 vel_p(i,j,k,3), eps_p(i,j,k), press_p(i,j,k), &
                 Bvec_p(i,j,k,1),Bvec_p(i,j,k,2),Bvec_p(i,j,k,3),w_lorentz_p(i,j,k))
          endif
        end if
        if (timelevels .gt. 2) then
          if (rho_p_p(i,j,k) .le. GRHydro_rho_min) then
            rho_p_p(i,j,k) = GRHydro_rho_min
            vel_p_p(i,j,k,1) = 0.0d0
            vel_p_p(i,j,k,2) = 0.0d0
            vel_p_p(i,j,k,3) = 0.0d0
          call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                 rho_p_p(i,j,k),eps_p_p(i,j,k),xtemp,xye,press_p_p(i,j,k),keyerr,anyerr)
          call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                 rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr)

            det = SPATIAL_DETERMINANT(gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), \
                  gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k))
            call prim2conpolytypeM(eos_handle, &
                 gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), &
                 gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k), det, &
                 dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), &
                 tau_p_p(i,j,k), Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),&
                 rho_p_p(i,j,k), vel_p_p(i,j,k,1), vel_p_p(i,j,k,2), &
                 vel_p_p(i,j,k,3), eps_p_p(i,j,k), press_p_p(i,j,k), &
                 Bvec_p_p(i,j,k,1),Bvec_p_p(i,j,k,2),Bvec_p_p(i,j,k,3),w_lorentz_p_p(i,j,k))
          endif
        endif

      end do
    end do
  end do

  write(*,*) "     GRHydro_InitialAtmosphereReset"
!!$  call GRHydro_BoundariesM(CCTK_PASS_FTOF)

end subroutine GRHydro_InitialAtmosphereResetM