/*@@
   @file      GRHydro_UpdateMask.F90
   @date      Wed Mar 13 14:18:38 2002
   @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"

#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 velx_p(i,j,k) vel_p(i,j,k,1)
#define vely_p(i,j,k) vel_p(i,j,k,2)
#define velz_p(i,j,k) vel_p(i,j,k,3)
#define velx_p_p(i,j,k) vel_p_p(i,j,k,1)
#define vely_p_p(i,j,k) vel_p_p(i,j,k,2)
#define velz_p_p(i,j,k) vel_p_p(i,j,k,3)

subroutine GRHydroUpdateAtmosphereMask(CCTK_ARGUMENTS)

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  CCTK_INT :: i,j,k
  CCTK_REAL :: frac

  CCTK_INT :: type_bits, atmosphere

  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
  call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")

  frac = CCTK_DELTA_TIME

  if(evolve_temper.ne.1) then
     !$OMP PARALLEL DO PRIVATE(j,i)
     do k = 1, cctk_lsh(3)
        do j = 1, cctk_lsh(2)
           do i = 1, cctk_lsh(1)
              if ( GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0) .or. &
                   (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \
                 type_bits, atmosphere)) .or. &
                      (tau(i,j,k) + frac * taurhs(i,j,k) .le. 0.d0) .or. &
                      (dens(i,j,k) + frac * densrhs(i,j,k) .le. 0.d0) ) then
                 densrhs(i,j,k) = 0.0d0
                 srhs(i,j,k,:)   = 0.0d0
                 taurhs(i,j,k)  = 0.0d0
                 atmosphere_mask(i,j,k) = 1
                 
                 SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
              end if
           end do
        end do
     end do
     !$OMP END PARALLEL DO
  else
     ! allow negative total energy
     !$OMP PARALLEL DO PRIVATE(j,i)
     do k = 1, cctk_lsh(3)
        do j = 1, cctk_lsh(2)
           do i = 1, cctk_lsh(1)
              if ( GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0) .or. &
                   (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \
                 type_bits, atmosphere)) .or. &
                      (dens(i,j,k) + frac * densrhs(i,j,k) .le. 0.d0) ) then
                 densrhs(i,j,k) = 0.0d0
                 srhs(i,j,k,:)   = 0.0d0
                 taurhs(i,j,k)  = 0.0d0
                 atmosphere_mask(i,j,k) = 1
                 
                 SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
              end if
           end do
        end do
     end do
     !$OMP END PARALLEL DO
  endif

end subroutine GRHydroUpdateAtmosphereMask

 /*@@
   @routine    GRHydro_SetupMask
   @date       Thu Jun 20 13:27:28 2002
   @author     Ian Hawke
   @desc 
   Initialize the mask to be zero.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine GRHydro_SetupMask(CCTK_ARGUMENTS)

  implicit none

  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  CCTK_INT :: type_bits, not_atmosphere, i, j, k

  ! Initialize all rhs variables and the mask.
  ! The former vars need to be initialized since there is
  ! no rhs computation in CCTK_INITIAL or POSTINITIAL.

  densrhs = 0.0d0
  taurhs  = 0.0d0 
  srhs    = 0.0d0

  if (evolve_tracer .ne. 0) then
     cons_tracerrhs = 0.0d0
  endif

  if (evolve_y_e .ne. 0) then
     y_e_con_rhs = 0.0d0
  endif

  atmosphere_mask = 0
  
  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
  call SpaceMask_GetStateBits(not_atmosphere, &
       &"Hydro_Atmosphere", "not_in_atmosphere")

  do k = 1, cctk_lsh(3)
    do j = 1, cctk_lsh(2)
      do i = 1, cctk_lsh(1)
  
        SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, not_atmosphere)

      end do
    end do
  end do

  call CCTK_INFO("Setting up the atmosphere mask: all points are not_atmosphere")
  
end subroutine GRHydro_SetupMask

 /*@@
   @routine    GRHydro_InitAtmosMask
   @date       Thu Jun 20 13:27:28 2002
   @author     Ian Hawke
   @desc 
   Initialize the mask based on rho_min. This is used only if wk_atmosphere=yes.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine GRHydro_InitAtmosMask(CCTK_ARGUMENTS)

  implicit none

  DECLARE_CCTK_ARGUMENTS

  CCTK_INT :: type_bits, atmosphere, i, j, k

  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
  call SpaceMask_GetStateBits(atmosphere, &
       &"Hydro_Atmosphere", "in_atmosphere")

  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
  
          SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
          atmosphere_mask(i,j,k) = 1

        end if
        
      end do
    end do
  end do

  call CCTK_INFO("Setting up the atmosphere mask: points with rho<rho_min are set to be atmosphere")
  
end subroutine GRHydro_InitAtmosMask

 /*@@
   @routine    GRHydro_AtmosphereReset
   @date       Thu Jun 20 13:30:51 2002
   @author     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_AtmosphereReset(CCTK_ARGUMENTS)

  implicit none

  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

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

  CCTK_INT :: type_bits, atmosphere, not_atmosphere

  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
  call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",&
                              "in_atmosphere")
  call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",&
                              "not_in_atmosphere")

  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) &
             &.or. (SpaceMask_CheckStateBitsF90(space_mask,i, j, k, type_bits,\
                    atmosphere)) &
             &) then

!!$          write(*,*) 'Resetting at ',i,j,k, atmosphere_mask(i, j, k), &
!!$             & (SpaceMask_CheckStateBitsF90(space_mask,i, j, k, type_bits,\
!!$                    atmosphere)) 

          rho(i,j,k) = GRHydro_rho_min
          velx(i,j,k) = 0.0d0
          vely(i,j,k) = 0.0d0
          velz(i,j,k) = 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 prim2conpolytype(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), rho(i,j,k), velx(i,j,k), vely(i,j,k), &
               velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k))
          if (wk_atmosphere .eq. 0) then
            atmosphere_mask(i, j, k) = 0
            SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\
                                      not_atmosphere)
          end if

        end if

      end do
    end do
  end do


!!$  call GRHydro_Boundaries(CCTK_PASS_FTOF)

end subroutine GRHydro_AtmosphereReset

subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)

  implicit none

  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

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

  CCTK_INT :: type_bits, atmosphere, not_atmosphere
  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
  
  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
  call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",&
                              "in_atmosphere")
  call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",&
                              "not_in_atmosphere")

  rho_min = GRHydro_rho_min
  if (initial_atmosphere_factor .gt. 0) then
    rho_min = rho_min * initial_atmosphere_factor
  endif

  do k = 1, cctk_lsh(3)
    do j = 1, cctk_lsh(2)
      do i = 1, cctk_lsh(1)
        
        if (rho(i,j,k) .le. rho_min) then
          rho(i,j,k) = rho_min
          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))

          velx(i,j,k) = 0.0d0
          vely(i,j,k) = 0.0d0
          velz(i,j,k) = 0.0d0

          if(evolve_temper.ne.0) then
             ! set the temperature to be relatively low
             temperature(i,j,k) = 0.1d0
             keytemp = 1
             call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                  rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
                  press(i,j,k),keyerr,anyerr)

             call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,&
                  i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),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), rho(i,j,k), velx(i,j,k), vely(i,j,k), &
                  velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k),&
                  temperature(i,j,k),y_e(i,j,k))
          else
             call EOS_Omni_press(GRHydro_polytrope_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(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                  rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
             call prim2conpolytype(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), rho(i,j,k), velx(i,j,k), vely(i,j,k), &
                  velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k))
          endif
        end if
        if (timelevels .gt. 1) then
          if (rho_p(i,j,k) .le. rho_min) then
            rho_p(i,j,k) = rho_min
            velx_p(i,j,k) = 0.0d0
            vely_p(i,j,k) = 0.0d0
            velz_p(i,j,k) = 0.0d0

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

            if(evolve_temper.ne.0) then
             ! set the temperature to be relatively low
             temperature_p(i,j,k) = 0.1d0
             keytemp = 1
             call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                  rho_p(i,j,k),eps_p(i,j,k),temperature_p(i,j,k),y_e_p(i,j,k),&
                  press_p(i,j,k),keyerr,anyerr)

             call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,&
                  i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),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), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), &
                  velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k),&
                  temperature_p(i,j,k),y_e_p(i,j,k))
            else
               call EOS_Omni_press(GRHydro_polytrope_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(GRHydro_polytrope_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)
               call prim2conpolytype(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), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), &
                    velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k))
            endif

          endif
        end if
        if (timelevels .gt. 2) then
          if (rho_p_p(i,j,k) .le. rho_min) then
            rho_p_p(i,j,k) = rho_min
            velx_p_p(i,j,k) = 0.0d0
            vely_p_p(i,j,k) = 0.0d0
            velz_p_p(i,j,k) = 0.0d0
            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))

            if(evolve_temper.ne.0) then
             ! set the temperature to be relatively low
             temperature_p_p(i,j,k) = 0.1d0
             keytemp = 1
             call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                  rho_p_p(i,j,k),eps_p_p(i,j,k),temperature_p_p(i,j,k),y_e_p_p(i,j,k),&
                  press_p_p(i,j,k),keyerr,anyerr)
             call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,&
                  i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),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), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), &
                  velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k),&
                  temperature_p_p(i,j,k),y_e_p_p(i,j,k))
            else
               call EOS_Omni_press(GRHydro_polytrope_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(GRHydro_polytrope_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)
               call prim2conpolytype(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), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), &
                    velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k))
            endif
          endif
        endif

      end do
    end do
  end do

end subroutine GRHydro_InitialAtmosphereReset