1 files changed, 114 insertions, 0 deletions
diff --git a/src/GRHydro_Only_Atmo.F90 b/src/GRHydro_Only_Atmo.F90
new file mode 100644
index 0000000..c2fb989
--- /dev/null
+++ b/src/GRHydro_Only_Atmo.F90
@@ -0,0 +1,114 @@
+ /*@@
+   @file      GRHydro_Only_Atmo.F90
+   @date      Sat Jan 29 04:44:25 2002
+   @author    Luca Baiotti
+   @desc 
+   Hydro initial data for vacuum.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.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 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_Only_Atmo
+   @date       Sat Jan 29 04:53:43 2002
+   @author     Luca Baiotti
+   @desc 
+   Hydro initial data for vacuum.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_Only_Atmo(CCTK_ARGUMENTS)
+
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  
+  CCTK_INT i,j,k,nx,ny,nz
+  CCTK_REAL det, atmo
+
+  call CCTK_INFO("Setting only atmosphere as initial data.")
+
+  
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+  
+  if (rho_abs_min > 0.0) then
+    atmo = rho_abs_min
+    GRHydro_rho_min = rho_abs_min
+  else 
+    call CCTK_WARN(0,"For initial data 'only_atmo' the parameter 'rho_abs_min' must be set")
+  end if
+    
+  if (initial_rho_abs_min > 0.0) then
+    atmo = initial_rho_abs_min
+  else if (initial_rho_rel_min > 0.0) then
+    call CCTK_WARN(3,"For initial data 'only_atmo' the parameter 'initial_rho_rel_min' is ignored")
+  end if
+
+  if (initial_atmosphere_factor > 0.0) then
+    atmo = atmo * initial_atmosphere_factor
+  end if
+
+  rho = atmo
+  if (rho_abs_min < initial_rho_abs_min) then
+    velx(:,:,:) = atmosphere_vel(1)
+    vely(:,:,:) = atmosphere_vel(2)
+    velz(:,:,:) = atmosphere_vel(3)
+  else
+    velx(:,:,:) = 0.d0
+    vely(:,:,:) = 0.d0
+    velz(:,:,:) = 0.d0
+  end if
+  if (attenuate_atmosphere .ne. 0) then
+    velx(:,:,:) = velx(:,:,:) * (1.d0 - exp(-2.d0 * (r - 1.d0)**2))
+    vely(:,:,:) = vely(:,:,:) * (1.d0 - exp(-2.d0 * (r - 1.d0)**2))
+    velz(:,:,:) = velz(:,:,:) * (1.d0 - exp(-2.d0 * (r - 1.d0)**2))
+  end if
+  eps = GRHydro_eps_min
+
+  do i=1,nx
+    do j=1,ny
+      do k=1,nz
+
+!!$        if (x(i,j,k) < 0.0) then
+!!$          velx(i,j,k) = -velx(i,j,k)
+!!$          vely(i,j,k) = -vely(i,j,k)
+!!$          velz(i,j,k) = -velz(i,j,k)
+!!$        end if
+        
+        call SpatialDeterminant(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)
+        call prim2con(GRHydro_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),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+             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))
+      enddo
+    enddo
+  enddo
+
+  densrhs = 0.d0
+  srhs = 0.d0
+  taurhs = 0.d0
+
+  return
+  
+end subroutine GRHydro_Only_Atmo