use initial_atmosphere_factor

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

1 files changed, 18 insertions, 12 deletions

diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90
index 7620ff1..e13c165 100644
--- a/src/GRHydro_UpdateMask.F90
+++ b/src/GRHydro_UpdateMask.F90
@@ -249,6 +249,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
 
   CCTK_INT :: i, j, k
   CCTK_REAL :: det, psi4pt
+  CCTK_REAL :: rho_min
 
   CCTK_INT :: type_bits, atmosphere, not_atmosphere
   CCTK_INT :: eos_handle
@@ -281,12 +282,17 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
   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. GRHydro_rho_min) then
-          rho(i,j,k) = GRHydro_rho_min
+        if (rho(i,j,k) .le. rho_min) then
+          rho(i,j,k) = rho_min
           velx(i,j,k) = 0.0d0
           vely(i,j,k) = 0.0d0
           velz(i,j,k) = 0.0d0
@@ -296,8 +302,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
           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)
 #else
-          press(i,j,k) = EOS_Pressure(eos_handle, GRHydro_rho_min, eps(i,j,k))
-          eps(i,j,k)   = EOS_SpecificIntEnergy(eos_handle, GRHydro_rho_min, press(i,j,k))
+          press(i,j,k) = EOS_Pressure(eos_handle, rho_min, eps(i,j,k))
+          eps(i,j,k)   = EOS_SpecificIntEnergy(eos_handle, rho_min, press(i,j,k))
 #endif
           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))
@@ -309,8 +315,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
                velz(i,j,k), eps(i,j,k), press(i,j,k), 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
+          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
@@ -320,8 +326,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
           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)
 #else
-            press_p(i,j,k) = EOS_Pressure(eos_handle, GRHydro_rho_min, eps_p(i,j,k))
-            eps_p(i,j,k)   = EOS_SpecificIntEnergy(eos_handle, GRHydro_rho_min, press_p(i,j,k))
+            press_p(i,j,k) = EOS_Pressure(eos_handle, rho_min, eps_p(i,j,k))
+            eps_p(i,j,k)   = EOS_SpecificIntEnergy(eos_handle, rho_min, press_p(i,j,k))
 #endif
             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))
@@ -334,8 +340,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
           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
+          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
@@ -345,8 +351,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
           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)
 #else
-            press_p_p(i,j,k) = EOS_Pressure(eos_handle, GRHydro_rho_min, eps_p_p(i,j,k))
-            eps_p_p(i,j,k)   = EOS_SpecificIntEnergy(eos_handle, GRHydro_rho_min, press_p_p(i,j,k))
+            press_p_p(i,j,k) = EOS_Pressure(eos_handle, rho_min, eps_p_p(i,j,k))
+            eps_p_p(i,j,k)   = EOS_SpecificIntEnergy(eos_handle, rho_min, press_p_p(i,j,k))
 #endif
             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))