diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/nuc_eos/nuc_eos.F90 | 45 |
1 files changed, 32 insertions, 13 deletions
diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90 index b576ce0..64a78a5 100644 --- a/src/nuc_eos/nuc_eos.F90 +++ b/src/nuc_eos/nuc_eos.F90 @@ -160,6 +160,25 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& end subroutine nuc_eos_full +subroutine nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + + implicit none + real*8 xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe + real*8,parameter:: polyK1 = 5.0d14 + real*8,parameter :: polyGamma = 1.33333333333333d0 + integer keytemp + + xprs=polyK1*xrho**(polyGamma) + xenr=xprs/xrho/(polyGamma-1.d0) + xcs2=polyGamma*xprs/xrho + xdpderho=0.0d0 + xdpdrhoe=0.0d0 + + xenr=max(8.0d18,xenr) + xenr=min(3.0d20,xenr) + +end subroutine nuc_poly_eos + subroutine nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) implicit none @@ -212,9 +231,9 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& call CCTK_WARN (0, "nuc_eos: rho > rhomax") endif - if(xrho.lt.eos_rhomin*1.2d0) then - call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) - xent = 4.0d0 + if(xrho.lt.eos_rhomin*1.126d0) then + call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + xent = 0.0d0 return endif @@ -234,8 +253,8 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& endif if(xtemp.lt.eos_tempmin) then - call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) - xent = 4.0d0 + call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + xent = 0.0d0 return endif endif @@ -335,8 +354,8 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& call CCTK_WARN (0, "nuc_eos: rho > rhomax") endif - if(xrho.lt.eos_rhomin*1.2d0) then - call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + if(xrho.lt.eos_rhomin*1.126d0) then + call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) return endif @@ -354,9 +373,9 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& if(xtemp.gt.eos_tempmax) then call CCTK_WARN (0, "nuc_eos: temp > tempmax") endif - + if(xtemp.lt.eos_tempmin) then - call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) return endif endif @@ -414,8 +433,8 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& call CCTK_WARN (0, "nuc_eos: rho > rhomax") endif - if(xrho.lt.eos_rhomin*1.2d0) then - call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + if(xrho.lt.eos_rhomin*1.126d0) then + call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) return endif @@ -433,9 +452,9 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& if(xtemp.gt.eos_tempmax) then call CCTK_WARN (0, "nuc_eos: temp > tempmax") endif - + if(xtemp.lt.eos_tempmin) then - call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) return endif endif |