diff options
Diffstat (limited to 'src/EOS_Hybrid_Analysis.F')
-rw-r--r-- | src/EOS_Hybrid_Analysis.F | 76 |
1 files changed, 76 insertions, 0 deletions
diff --git a/src/EOS_Hybrid_Analysis.F b/src/EOS_Hybrid_Analysis.F new file mode 100644 index 0000000..9f82a41 --- /dev/null +++ b/src/EOS_Hybrid_Analysis.F @@ -0,0 +1,76 @@ + /*@@ + @file EOS_Hybrid_Analysis.F + @date Fri Apr 26 16:14:47 2002 + @author Harry Dimmelmeier, Ian Hawke, Christian Ott + @desc + Calculates the polytropic and thermal contributions to the + total pressure. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + + /*@@ + @routine Check_Poly_Contrib + @date Fri Apr 26 16:25:35 2002 + @author Harry Dimmelmeier, Ian Hawke, Christian Ott + @desc + The routine that calculates the contributions. + @enddesc + @calls + @calledby + @history + + @endhistory +@@*/ + + subroutine Check_Poly_Contrib(CCTK_ARGUMENTS) + + USE EOS_Polytrope_Scalars + USE EOS_Hybrid_Scalars + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: i,j,k + + + CCTK_REAL local_eos_gamma, local_eos_k_cgs, d_p_poly, d_p_th_1, + . d_p_th_2, zero + + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + + if (rho(i,j,k) > rho_nuc) then + local_eos_gamma = eos_gamma_supernuclear + local_eos_k_cgs = eos_k_supernuclear_cgs + else + local_eos_gamma = eos_gamma + local_eos_k_cgs = eos_k_cgs + end if + + pressure_poly(i,j,k) = p_geom_factor * local_eos_k_cgs * + . (rho(i,j,k) * rho_geom_factor_inv)**local_eos_gamma + + pressure_th(i,j,k) = - p_geom_factor * local_eos_k_cgs * + . (eos_gamma_th - 1.d0) / (local_eos_gamma - 1.d0) * + . (rho(i,j,k) * rho_geom_factor_inv)**local_eos_gamma + + . (eos_gamma_th - 1.d0) * rho(i,j,k) * eps(i,j,k) - + . (eos_gamma_th - 1.d0) * (local_eos_gamma - eos_gamma) / + . (eos_gamma - 1.d0) / (eos_gamma_supernuclear - 1.d0) * + . p_geom_factor * eos_k_cgs * + . rho_geom_factor_inv**eos_gamma * + . rho_nuc**(eos_gamma - 1.d0) * rho(i,j,k) + +c$$$ pressure_th(i,j,k) = dmax1(0.d0, pressure_th(i,j,k)) + + end do + end do + end do + + end subroutine Check_Poly_Contrib |