* extrapolate to negative (!) temperatures, so that we can use this

  with the new Con2PrimHot treatment in GRHydro.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Omni/trunk@78 8e189c6b-2ab8-4400-aa02-70a9cfce18b9

3 files changed, 482 insertions, 63 deletions

diff --git a/src/nuc_eos/findtemp.F90 b/src/nuc_eos/findtemp.F90
index 646994d..125b0cd 100644
--- a/src/nuc_eos/findtemp.F90
+++ b/src/nuc_eos/findtemp.F90
@@ -16,10 +16,10 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,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
+  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 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps)
   eps0=epsin
   eps1=eps0
   
-
   ltmax=logtemp(ntemp)
   ltmin=logtemp(1)
 
@@ -42,25 +41,17 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps)
   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 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps)
      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 @@ subroutine findtemp_entropy(lr,lt0,y,sin,keyerrt,rfeps)
  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 @@ subroutine findtemp_entropy(lr,lt0,y,sin,keyerrt,rfeps)
           write(*,*) "Tried calling bisection... didn't help... :-/"
           write(*,*) "Bisection error: ",keyerrt
     endif
+#endif
     
     lt0=lt
     return
@@ -207,3 +194,89 @@ subroutine findtemp_entropy(lr,lt0,y,sin,keyerrt,rfeps)
 
 
 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
diff --git a/src/nuc_eos/linterp_many.F90 b/src/nuc_eos/linterp_many.F90
index d0ddb08..42b819f 100644
--- a/src/nuc_eos/linterp_many.F90
+++ b/src/nuc_eos/linterp_many.F90
@@ -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
diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90
index 3420ee0..1dfb9e8 100644
--- a/src/nuc_eos/nuc_eos.F90
+++ b/src/nuc_eos/nuc_eos.F90
@@ -69,16 +69,18 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
 
   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 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
   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 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
      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 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
      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 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
   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 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,&
 
   ! 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 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,&
      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 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,&
      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 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,&
 
   ! 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 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,&
      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 @@ subroutine findthis(lr,lt,y,value,array,d1,d2,d3)
 ! 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 @@ subroutine findall_short(lr,lt,y,ff)
 
 end subroutine findall_short
 
-
-subroutine findall_press_eps(lr,lt,y,ff)
+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(2,1)
+  real*8 ff(2)
   real*8 lr,lt,y
   integer i
   integer :: nvarsx = 2
@@ -577,13 +639,32 @@ subroutine findall_press_eps(lr,lt,y,ff)
 
 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 @@ subroutine findall_dpdr_dpde(lr,lt,y,ff)
   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