User: cott Date: 2013/03/23 03:10 PM
Modified: /trunk/src/nuc_eos/ findtemp.F90, linterp_many.F90, nuc_eos.F90
Log: * extrapolate to negative (!) temperatures, so that we can use this with the new Con2PrimHot treatment in GRHydro.
File Changes:
Directory: /trunk/src/nuc_eos/ ==============================
File [modified]: findtemp.F90 Delta lines: +90 -17 =================================================================== --- trunk/src/nuc_eos/findtemp.F90 2013-03-23 20:08:12 UTC (rev 77) +++ trunk/src/nuc_eos/findtemp.F90 2013-03-23 20:10:20 UTC (rev 78) @@ -16,10 +16,10 @@ integer :: rl = 0 integer itmax,i,keyerrt integer ii,jj,kk - + keyerrt=0
- tol=rfeps ! need to find energy to less than 1 in 10^10 + tol=rfeps ! precision to which we need to find temp itmax=20 ! use at most 20 iterations, then bomb
lt=lt0 @@ -28,7 +28,6 @@ eps0=epsin eps1=eps0
- ltmax=logtemp(ntemp) ltmin=logtemp(1)
@@ -42,25 +41,17 @@ lt1=lt eps1=eps
-! write(*,"(i4,1P12E19.10)") 0,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0,d2 -! write(*,"(i4,1P12E19.10)") 0,lr,lt0,y,ltmin,ltmax - do i=1,itmax !d2 is the derivative deps/dlogtemp; ldt = -(eps - eps0)/d2 -! if(ldt.gt.0.d0) ltmin=lt -! if(ldt.le.0.d0) ltmax=lt ltn = lt+ldt ltn = min(ltn,ltmax) ltn = max(ltn,ltmin) lt1=lt lt=ltn eps1=eps -! write(*,"(i4,1P12E19.10)") i,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0,d2,ldt call findthis(lr,lt,y,eps,alltables(:,:,:,2),d1,d2,d3) -! call findthis(lr,lt,y,eps,epst,d1,d2,d3) if (abs(eps - eps0).lt.tol*abs(eps0)) then -! write(*,"(1P12E19.10)") tol,abs(eps-eps0)/eps0 exit endif !setup new d2 @@ -72,12 +63,6 @@ if(abs(eps-eps0).lt.1.0d-3*abs(eps0)) then d2 = (eps-eps1)/(lt-lt1) endif -! if(i.ge.10) then -! write(*,*) "EOS: Did not converge in findtemp!" -! write(*,*) "rl,logrho,logtemp0,ye,lt,eps,eps0,abs(eps-eps0)/eps0" -! if(i.gt.5) then -! write(*,"(i4,1P10E22.14)") i,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0 -! endif enddo
@@ -190,6 +175,7 @@ if(i.ge.itmax) then keyerrt=667 call bisection(lr,lt0,y,s0,lt,alltables(:,:,:,3),keyerrt,2) +#if 0 if(keyerrt.eq.667) then write(*,*) "EOS: Did not converge in findtemp_entropy!" write(*,*) "rl,logrho,logtemp0,ye,lt,s,s0,abs(s-s0)/s0" @@ -197,6 +183,7 @@ write(*,*) "Tried calling bisection... didn't help... :-/" write(*,*) "Bisection error: ",keyerrt endif +#endif
lt0=lt return @@ -207,3 +194,89 @@
end subroutine findtemp_entropy + + +subroutine findtemp_low(lr,t0,y,epsin,keyerrt,rfeps) + + ! this routine is for finding fake temperatures + ! outside of the table range; we use linear + ! extrapolation in temperature + + use eosmodule + + implicit none + + real*8 lr,t0,y,epsin + real*8 eps,t,dt + real*8 tol + real*8 dlepsdt + real*8 eps0,eps1,t1 + + real*8 tn,tmin,tmax + real*8 tinput,rfeps + + integer :: rl = 0 + integer itmax,i,keyerrt + integer ii,jj,kk + + keyerrt=0 + + tol=rfeps ! need to find energy to less than 1 in 10^-10 + itmax=100 ! use at most 20 iterations, then bomb + + t=t0 + t1=t + + eps0=epsin + eps1=eps0 + + tmax=10.0d0**logtemp(2) + tmin=-20.0d0 + + !preconditioning 1: do we already have the right temperature? + call intp3d_linearTlow(lr,t,y,eps,1,alltables(:,:,:,2),nrho,& + ntemp,nye,logrho,logtemp,ye,dlepsdt) + + if (abs(eps-eps0).lt.tol*abs(eps0)) then + return + endif + t1=t + eps1=eps + + do i=1,itmax + dt = -(eps - eps0)/dlepsdt + tn = t+dt + if(tn >= tmax) then + tn = tmax*0.9d0 + endif + tn = min(tn,tmax) + tn = max(tn,tmin) + t1=t + t=tn + eps1=eps + + call intp3d_linearTlow(lr,t,y,eps,1,alltables(:,:,:,2),nrho,& + ntemp,nye,logrho,logtemp,ye,dlepsdt) + + if (abs(eps - eps0).lt.tol*abs(eps0)) then + exit + endif + + ! if we are closer than 10^-2 to the + ! root (eps-eps0)=0, we are switching to + ! the secant method, since the table is rather coarse and the + ! derivatives may be garbage. + if(abs(eps-eps0).lt.1.0d-3*abs(eps0)) then + dlepsdt = (eps-eps1)/(t-t1) + endif + enddo + + t0 = t + + if(i.ge.itmax) then + keyerrt=667 + return + endif + + +end subroutine findtemp_low
File [modified]: linterp_many.F90 Delta lines: +250 -4 =================================================================== --- trunk/src/nuc_eos/linterp_many.F90 2013-03-23 20:08:12 UTC (rev 77) +++ trunk/src/nuc_eos/linterp_many.F90 2013-03-23 20:10:20 UTC (rev 78) @@ -1,5 +1,3 @@ -#include "cctk.h" - SUBROUTINE intp3d_many ( x, y, z, f, kt, ft, nx, ny, nz, nvars, xt, yt, zt) ! implicit none @@ -45,7 +43,7 @@ real*8 dx,dy,dz,dxi,dyi,dzi,dxyi,dxzi,dyzi,dxyzi integer n,ix,iy,iz
- IF (kt .GT. ktx) call CCTK_WARN(0, '***KTX**') + IF (kt .GT. ktx) STOP'***KTX**' ! ! !------ determine spacing parameters of (equidistant!!!) table @@ -67,7 +65,7 @@ ! !------- loop over all points to be interpolated ! - do n = 1, kt + dO n = 1, kt ! !------- determine location in (equidistant!!!) table ! @@ -125,3 +123,251 @@
end SUBROUTINE intp3d_many
+ + SUBROUTINE intp3d_linearTlow ( x, y, z, f, kt, ft, nx, ny, nz, xt, yt, zt, & + d2) +! + implicit none +! +!--------------------------------------------------------------------- +! +! purpose: interpolation of a function of three variables in an +! equidistant(!!!) table. +! +! method: 8-point Lagrange linear interpolation formula +! +! x input vector of first variable +! y temperature +! z input vector of third variable +! +! f output vector of interpolated function values +! +! kt vector length of input and output vectors +! +! ft 3d array of tabulated function values +! nx x-dimension of table +! ny y-dimension of table +! nz z-dimension of table +! xt vector of x-coordinates of table +! yt vector of y-coordinates of table +! zt vector of z-coordinates of table +! +!--------------------------------------------------------------------- + + integer kt,nx,ny,nz + real*8 :: ft(nx,ny,nz) + + real*8 x(kt),y(kt),z(kt),f(kt) + real*8 xt(nx),yt(ny),zt(nz) + real*8 d1,d2,d3 +! +! + integer,parameter :: ktx = 1 + real*8 fh(ktx,8), delx(ktx), dely(ktx), delz(ktx), & + a1(ktx), a2(ktx), a3(ktx), a4(ktx), & + a5(ktx), a6(ktx), a7(ktx), a8(ktx) + + real*8 dx,dy,dz,dxi,dyi,dzi,dxyi,dxzi,dyzi,dxyzi + integer n,ix,iy,iz + + IF (kt .GT. ktx) STOP'***KTX**' +! +! +!------ determine spacing parameters of (equidistant!!!) table +! + dx = (xt(nx) - xt(1)) / FLOAT(nx-1) + dy = 10.0d0**yt(2) - 10.0d0**yt(1) + dz = (zt(nz) - zt(1)) / FLOAT(nz-1) +! + dxi = 1. / dx + dyi = 1. / dy + dzi = 1. / dz +! + dxyi = dxi * dyi + dxzi = dxi * dzi + dyzi = dyi * dzi +! + dxyzi = dxi * dyi * dzi +! +! +!------- loop over all points to be interpolated +! + dO n = 1, kt +! +!------- determine location in (equidistant!!!) table +! + ix = 2 + INT( (x(n) - xt(1) - 1.e-10) * dxi ) + iy = 2 + iz = 2 + INT( (z(n) - zt(1) - 1.e-10) * dzi ) +! + ix = MAX( 2, MIN( ix, nx ) ) + iz = MAX( 2, MIN( iz, nz ) ) +! +! write(*,*) iy-1,iy,iy+1 +! +!------- set-up auxiliary arrays for Lagrange interpolation +! + delx(n) = xt(ix) - x(n) + dely(n) = 10.0d0**yt(iy) - y(n) + delz(n) = zt(iz) - z(n) +! + fh(n,1) = ft(ix , iy , iz) + fh(n,2) = ft(ix-1, iy , iz) + fh(n,3) = ft(ix , iy-1, iz) + fh(n,4) = ft(ix , iy , iz-1) + fh(n,5) = ft(ix-1, iy-1, iz) + fh(n,6) = ft(ix-1, iy , iz-1) + fh(n,7) = ft(ix , iy-1, iz-1) + fh(n,8) = ft(ix-1, iy-1, iz-1) +! +!------ set up coefficients of the interpolation polynomial and +! evaluate function values + ! + a1(n) = fh(n,1) + a2(n) = dxi * ( fh(n,2) - fh(n,1) ) + a3(n) = dyi * ( fh(n,3) - fh(n,1) ) + a4(n) = dzi * ( fh(n,4) - fh(n,1) ) + a5(n) = dxyi * ( fh(n,5) - fh(n,2) - fh(n,3) + fh(n,1) ) + a6(n) = dxzi * ( fh(n,6) - fh(n,2) - fh(n,4) + fh(n,1) ) + a7(n) = dyzi * ( fh(n,7) - fh(n,3) - fh(n,4) + fh(n,1) ) + a8(n) = dxyzi * ( fh(n,8) - fh(n,1) + fh(n,2) + fh(n,3) + & + fh(n,4) - fh(n,5) - fh(n,6) - fh(n,7) ) +! + d2 = -a3(n) + f(n) = a1(n) + a2(n) * delx(n) & + + a3(n) * dely(n) & + + a4(n) * delz(n) & + + a5(n) * delx(n) * dely(n) & + + a6(n) * delx(n) * delz(n) & + + a7(n) * dely(n) * delz(n) & + + a8(n) * delx(n) * dely(n) * delz(n) +! + + enddo +! + + end SUBROUTINE intp3d_linearTlow + + SUBROUTINE intp3d_many_linearTLow ( x, y, z, f, kt, ft, nx, ny, nz, nvars, xt, yt, zt) +! + implicit none +! +!--------------------------------------------------------------------- +! +! purpose: interpolation of a function of three variables in an +! equidistant(!!!) table. +! +! method: 8-point Lagrange linear interpolation formula +! +! x input vector of first variable +! y input vector of second variable +! z input vector of third variable +! +! f output vector of interpolated function values +! +! kt vector length of input and output vectors +! +! ft 3d array of tabulated function values +! nx x-dimension of table +! ny y-dimension of table +! nz z-dimension of table +! xt vector of x-coordinates of table +! yt vector of y-coordinates of table +! zt vector of z-coordinates of table +! +!--------------------------------------------------------------------- + + integer kt,nx,ny,nz,iv,nvars + real*8 :: ft(nx,ny,nz,nvars) + + real*8 x(kt),y(kt),z(kt),f(kt,nvars) + real*8 xt(nx),yt(ny),zt(nz) + real*8 d1,d2,d3 +! +! + integer,parameter :: ktx = 1 + real*8 fh(ktx,8,nvars), delx(ktx), dely(ktx), delz(ktx), & + a1(ktx,nvars), a2(ktx,nvars), a3(ktx,nvars), a4(ktx,nvars), & + a5(ktx,nvars), a6(ktx,nvars), a7(ktx,nvars), a8(ktx,nvars) + + real*8 dx,dy,dz,dxi,dyi,dzi,dxyi,dxzi,dyzi,dxyzi + integer n,ix,iy,iz + + IF (kt .GT. ktx) STOP'***KTX**' +! +! +!------ determine spacing parameters of (equidistant!!!) table +! + dx = (xt(nx) - xt(1)) / FLOAT(nx-1) + dy = 10.0d0**yt(2) - 10.0d0**yt(1) + dz = (zt(nz) - zt(1)) / FLOAT(nz-1) +! + dxi = 1. / dx + dyi = 1. / dy + dzi = 1. / dz +! + dxyi = dxi * dyi + dxzi = dxi * dzi + dyzi = dyi * dzi +! + dxyzi = dxi * dyi * dzi +! +! +!------- loop over all points to be interpolated +! + dO n = 1, kt +! +!------- determine location in (equidistant!!!) table +! + ix = 2 + INT( (x(n) - xt(1) - 1.e-10) * dxi ) + iy = 2 + iz = 2 + INT( (z(n) - zt(1) - 1.e-10) * dzi ) +! + ix = MAX( 2, MIN( ix, nx ) ) + iz = MAX( 2, MIN( iz, nz ) ) +! +! write(*,*) iy-1,iy,iy+1 +! +!------- set-up auxiliary arrays for Lagrange interpolation +! + delx(n) = xt(ix) - x(n) + dely(n) = 10.0d0**yt(iy) - y(n) + delz(n) = zt(iz) - z(n) +! + do iv = 1, nvars + fh(n,1,iv) = ft(ix , iy , iz, iv ) + fh(n,2,iv) = ft(ix-1, iy , iz, iv ) + fh(n,3,iv) = ft(ix , iy-1, iz, iv ) + fh(n,4,iv) = ft(ix , iy , iz-1, iv) + fh(n,5,iv) = ft(ix-1, iy-1, iz, iv ) + fh(n,6,iv) = ft(ix-1, iy , iz-1, iv) + fh(n,7,iv) = ft(ix , iy-1, iz-1, iv) + fh(n,8,iv) = ft(ix-1, iy-1, iz-1, iv) +! +!------ set up coefficients of the interpolation polynomial and +! evaluate function values + ! + a1(n,iv) = fh(n,1,iv) + a2(n,iv) = dxi * ( fh(n,2,iv) - fh(n,1,iv) ) + a3(n,iv) = dyi * ( fh(n,3,iv) - fh(n,1,iv) ) + a4(n,iv) = dzi * ( fh(n,4,iv) - fh(n,1,iv) ) + a5(n,iv) = dxyi * ( fh(n,5,iv) - fh(n,2,iv) - fh(n,3,iv) + fh(n,1,iv) ) + a6(n,iv) = dxzi * ( fh(n,6,iv) - fh(n,2,iv) - fh(n,4,iv) + fh(n,1,iv) ) + a7(n,iv) = dyzi * ( fh(n,7,iv) - fh(n,3,iv) - fh(n,4,iv) + fh(n,1,iv) ) + a8(n,iv) = dxyzi * ( fh(n,8,iv) - fh(n,1,iv) + fh(n,2,iv) + fh(n,3,iv) + & + fh(n,4,iv) - fh(n,5,iv) - fh(n,6,iv) - fh(n,7,iv) ) +! + f(n,iv) = a1(n,iv) + a2(n,iv) * delx(n) & + + a3(n,iv) * dely(n) & + + a4(n,iv) * delz(n) & + + a5(n,iv) * delx(n) * dely(n) & + + a6(n,iv) * delx(n) * delz(n) & + + a7(n,iv) * dely(n) * delz(n) & + + a8(n,iv) * delx(n) * dely(n) * delz(n) +! + enddo + + enddo +! + + end SUBROUTINE intp3d_many_linearTlow
File [modified]: nuc_eos.F90 Delta lines: +142 -42 =================================================================== --- trunk/src/nuc_eos/nuc_eos.F90 2013-03-23 20:08:12 UTC (rev 77) +++ trunk/src/nuc_eos/nuc_eos.F90 2013-03-23 20:10:20 UTC (rev 78) @@ -69,16 +69,18 @@
if(keytemp.eq.1) then if(xtemp.gt.eos_tempmax) then - call CCTK_WARN (0, "nuc_eos: temp > tempmax") + keyerr = 105 + return endif
if(xtemp.lt.eos_tempmin) then - call CCTK_WARN (0, "nuc_eos: temp < tempmin") + keyerr = 106 + return endif endif
lr = log10(xrho) - lt = log10(xtemp) + lt = log10(max(xtemp,eos_tempmin)) y = xye xeps = xenr + energy_shift leps = log10(max(xeps,1.0d0)) @@ -221,7 +223,7 @@ integer, intent(out) :: keyerr
! local variables - real*8 :: lr,lt,y,xx,xeps,leps,xs,xpressure + real*8 :: lr,lt,y,xt,xx,xeps,leps,xs,xpressure real*8 :: d1,d2,d3,ff(8) integer :: keyerrt integer :: keyerrr @@ -256,14 +258,13 @@ endif
if(xtemp.lt.eos_tempmin) then - call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) - xent = 0.0d0 + keyerr = 106 return endif endif
lr = log10(xrho) - lt = log10(xtemp) + lt = log10(max(xtemp,eos_tempmin)) y = xye
keyerr = 0 @@ -273,12 +274,22 @@ xeps = xenr + energy_shift leps = log10(max(xeps,1.0d0)) call findtemp(lr,lt,y,leps,keyerrt,rfeps) - if(keyerrt.ne.0) then + if(keyerrt.eq.667 .and. lt .lt. 10.0d0) then + ! turns out, we can still converge, but + ! too a bogus, possibly negative temperature + xt = xtemp + call findtemp_low(lr,xt,y,leps,keyerrt,rfeps) + if(keyerrt.eq.0) then + keyerrt = 668 + xtemp = xt + else + keyerrt = 667 + return + endif keyerr = keyerrt - return + else + xtemp = 10.0d0**lt endif - xtemp = 10.0d0**lt - elseif(keytemp.eq.2) then !need to find temperature based on xent xs = xent @@ -302,7 +313,11 @@ endif
! have rho,temp,ye; proceed: - call findall_short(lr,lt,y,ff) + if(keyerr.ne.668) then + call findall_short(lr,lt,y,ff) + else + call findall_short_lowT(lr,xt,y,ff) + endif
!unless we want xprs to be constant (keytemp==3), reset xprs if(.not.keytemp.eq.3) then @@ -347,7 +362,7 @@
! local variables real*8 :: xcs2,xdpderho,xdpdrhoe - real*8 :: lr,lt,y,xx,xeps,leps,xs + real*8 :: lr,lt,y,xt,xx,xeps,leps,xs real*8 :: d1,d2,d3,ff(8) integer :: keyerrt
@@ -379,7 +394,7 @@ endif
if(xtemp.lt.eos_tempmin) then - call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + keyerr = 106 return endif endif @@ -390,29 +405,44 @@ call CCTK_WARN (0, "eos_nuc_press does not support keytemp other than 0 and 1") endif
- lr = log10(xrho) - lt = log10(xtemp) + lr = log10(xrho) + lt = log10(max(xtemp,eos_tempmin)) y = xye if(keytemp.eq.0) then !need to find temperature based on xeps xeps = xenr + energy_shift leps = log10(max(xeps,1.0d0)) call findtemp(lr,lt,y,leps,keyerrt,rfeps) - if(keyerrt.ne.0) then + if(keyerrt.eq.667 .and. lt .lt. 10.0d0) then + ! turns out, we can still converge, but + ! too a bogus, possibly negative temperature + xt = xtemp + call findtemp_low(lr,xt,y,leps,keyerrt,rfeps) + if(keyerrt.eq.0) then + keyerrt = 668 + xtemp = xt + else + keyerrt = 667 + return + endif keyerr = keyerrt - return + else + xtemp = 10.0d0**lt endif - xtemp = 10.0d0**lt endif ! have temperature, proceed: - call findall_press_eps(lr,lt,y,ff) + if(keyerr.ne.668) then + call findall_press_eps(lr,lt,y,ff) + else + call findall_press_eps_lowT(lr,xt,y,ff) + endif xprs = 10.0d0**ff(1) xenr = 10.0d0**ff(2) - energy_shift
end subroutine nuc_eos_press_eps
-subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& - xdpderho,& +subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpderho,& + xdpdrhoe,& keytemp,keyerr,rfeps)
use eosmodule @@ -427,8 +457,8 @@
! local variables real*8 :: xcs2,xprs - real*8 :: lr,lt,y,xx,xeps,leps,xs - real*8 :: d1,d2,d3,ff(8) + real*8 :: lr,lt,y,xt,xx,xeps,leps,xs + real*8 :: d1,d2,d3,ff(2) integer :: keyerrt
keyerrt = 0 @@ -459,38 +489,52 @@ endif
if(xtemp.lt.eos_tempmin) then - call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + keyerr = 106 return endif endif
lr = log10(xrho) - lt = log10(xtemp) + lt = log10(max(xtemp,eos_tempmin)) y = xye - xeps = xenr + energy_shift - leps = log10(max(xeps,1.0d0))
keyerr = 0 - if(keytemp.eq.0) then !need to find temperature based on xeps + xeps = xenr + energy_shift + leps = log10(max(xeps,1.0d0)) call findtemp(lr,lt,y,leps,keyerrt,rfeps) - if(keyerrt.ne.0) then + if(keyerrt.eq.667 .and. lt .lt. 10.0d0) then + ! turns out, we can still converge, but + ! too a bogus, possibly negative temperature + xt = xtemp + call findtemp_low(lr,xt,y,leps,keyerrt,rfeps) + if(keyerrt.eq.0) then + keyerrt = 668 + xtemp = xt + else + keyerrt = 667 + return + endif keyerr = keyerrt - return + else + xtemp = 10.0d0**lt endif - xtemp = 10.0d0**lt - elseif(keytemp.gt.1) then call CCTK_WARN (0, "eos_nuc_press does not support keytemp > 1") endif
! have temperature, proceed: - call findall_dpdr_dpde(lr,lt,y,ff) - xdpdrhoe = 10.0d0**ff(1) - xdpderho = 10.0d0**ff(2) + if(keyerr.ne.668) then + call findall_dpdr_dpde(lr,lt,y,ff) + else + call findall_dpdr_dpde_lowT(lr,xtemp,y,ff) + endif
- if(keytemp.eq.1) then + xdpdrhoe = ff(1) + xdpderho = ff(2) + + if(keytemp.eq.1.and.keyerr.ne.668) then call findthis(lr,lt,y,xeps,alltables(:,:,:,2),d1,d2,d3) xeps = 10.0d0**xeps - energy_shift xenr = xeps @@ -515,7 +559,6 @@ ! Ewald's interpolator call intp3d(lr,lt,y,value,1,array,nrho,ntemp,nye,logrho,logtemp,ye,d1,d2,d3)
- end subroutine findthis
@@ -557,14 +600,33 @@
end subroutine findall_short
+subroutine findall_short_lowT(lr,t,y,ff)
+ use eosmodule + implicit none + + real*8 ffx(8,1) + real*8 ff(8) + real*8 lr,t,y + integer i + integer :: nvarsx = 8 + + +! Ewald's interpolator + call intp3d_many_linearTlow(lr,t,y,ffx,1,alltables(:,:,:,1:8), & + nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) + ff(:) = ffx(:,1) + +end subroutine findall_short_lowT + + subroutine findall_press_eps(lr,lt,y,ff)
use eosmodule implicit none
- real*8 ffx(8,1) - real*8 ff(8) + real*8 ffx(2,1) + real*8 ff(2) real*8 lr,lt,y integer i integer :: nvarsx = 2 @@ -577,13 +639,32 @@
end subroutine findall_press_eps
+subroutine findall_press_eps_lowT(lr,t,y,ff) + + use eosmodule + implicit none + + real*8 ffx(2,1) + real*8 ff(2) + real*8 lr,t,y + integer i + integer :: nvarsx = 2 + +! Ewald's interpolator + call intp3d_many_linearTlow(lr,t,y,ffx,1,alltables(:,:,:,1:2), & + nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) + ff(:) = ffx(:,1) + +end subroutine findall_press_eps_lowT + + subroutine findall_dpdr_dpde(lr,lt,y,ff)
use eosmodule implicit none
- real*8 ffx(8,1) - real*8 ff(8) + real*8 ffx(2,1) + real*8 ff(2) real*8 lr,lt,y integer i integer :: nvarsx = 2 @@ -595,3 +676,22 @@ ff(:) = ffx(:,1)
end subroutine findall_dpdr_dpde + +subroutine findall_dpdr_dpde_lowT(lr,t,y,ff) + + use eosmodule + implicit none + + real*8 ffx(2,1) + real*8 ff(2) + real*8 lr,t,y + integer i + integer :: nvarsx = 2 + + +! Ewald's interpolator + call intp3d_many_linearTlow(lr,t,y,ffx,1,alltables(:,:,:,7:8), & + nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) + ff(:) = ffx(:,1) + +end subroutine findall_dpdr_dpde_lowT