diff options
Diffstat (limited to 'src/nuc_eos/nuc_eos.F90')
| -rw-r--r-- | src/nuc_eos/nuc_eos.F90 | 63 |
1 files changed, 48 insertions, 15 deletions
diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90 index 1dfb9e8..65b0562 100644 --- a/src/nuc_eos/nuc_eos.F90 +++ b/src/nuc_eos/nuc_eos.F90 @@ -83,13 +83,17 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& 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 - call findtemp(lr,lt,y,leps,keyerrt,rfeps) + if(xeps .gt. 0.0d0) then + leps = log10(xeps) + call findtemp(lr,lt,y,leps,keyerrt,rfeps) + else + keyerrt = 667 + endif if(keyerrt.ne.0) then call CCTK_WARN (0, "Did not find temperature") endif @@ -272,18 +276,23 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& 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.eq.667 .and. lt .lt. 10.0d0) then + if(xeps .gt. 0.0d0) then + leps = log10(xeps) + call findtemp(lr,lt,y,leps,keyerrt,rfeps) + else + keyerrt = 667 + endif + if(keyerrt.eq.667) 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) + call findtemp_low(lr,xt,y,xeps,keyerrt,rfeps) if(keyerrt.eq.0) then keyerrt = 668 xtemp = xt else keyerrt = 667 + keyerr = keyerrt return endif keyerr = keyerrt @@ -365,6 +374,7 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& real*8 :: lr,lt,y,xt,xx,xeps,leps,xs real*8 :: d1,d2,d3,ff(8) integer :: keyerrt + character(len=256) :: warnline keyerrt = 0 @@ -411,20 +421,38 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& 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.eq.667 .and. lt .lt. 10.0d0) then + if(xeps .gt. 0.0d0) then + leps = log10(xeps) + call findtemp(lr,lt,y,leps,keyerrt,rfeps) + else + keyerrt = 667 + endif +! if(xenr .lt. -9.00725347549139d20) then +! if(xenr.lt.-2.28659899d20) then +! write(warnline,"(A4,i5,1P10E15.6)") "UGGx",keyerrt,lr,lt,xt,xtemp,y,xeps,xenr +! call CCTK_WARN(1,warnline) +! endif + if(keyerrt.eq.667) 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) + call findtemp_low(lr,xt,y,xeps,keyerrt,rfeps) if(keyerrt.eq.0) then keyerrt = 668 + if(xeps.lt.0.0d0) then + keyerrt = 668 ! impossible to do anything + endif xtemp = xt else keyerrt = 667 + keyerr = keyerrt return endif +! if(xenr .lt. -9.00725347549139d20) then +! if(xenr.lt.-2.28659899d20) then +! write(warnline,"(A4,i5,1P10E15.6)") "UGG",keyerrt,lr,lt,xt,xtemp,y,xeps,xenr +! call CCTK_WARN(0,warnline) +! endif keyerr = keyerrt else xtemp = 10.0d0**lt @@ -502,20 +530,25 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpderho,& 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.eq.667 .and. lt .lt. 10.0d0) then + if(xeps .gt. 0.0d0) then + leps = log10(xeps) + call findtemp(lr,lt,y,leps,keyerrt,rfeps) + else + keyerrt = 667 + endif + if(keyerrt.eq.667) 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) + call findtemp_low(lr,xt,y,xeps,keyerrt,rfeps) if(keyerrt.eq.0) then keyerrt = 668 xtemp = xt else keyerrt = 667 + keyerr = keyerrt return - endif + endif keyerr = keyerrt else xtemp = 10.0d0**lt |