* updates needed to make TOVSolverHot work

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

3 files changed, 151 insertions, 39 deletions

diff --git a/src/nuc_eos/bisection.F90 b/src/nuc_eos/bisection.F90
index f12b702..fb0eb43 100644
--- a/src/nuc_eos/bisection.F90
+++ b/src/nuc_eos/bisection.F90
@@ -7,7 +7,7 @@ subroutine bisection(lr,lt0,y,eps0,lt,bivar,keyerrt,keybisect)
   real*8 lr,lt0,y,eps0,lt
   integer keyerrt
 
-  integer keybisect
+  integer keybisect !this doesn't do anything...
   ! 1 -> operate on energy
   ! 2 -> operate on entropy
 
@@ -19,7 +19,7 @@ subroutine bisection(lr,lt0,y,eps0,lt,bivar,keyerrt,keybisect)
 
   integer bcount,i,itmax,maxbcount
 
-  real*8 bivar(*)
+  real*8 bivar
 
   bcount = 0
   maxbcount = 80
@@ -80,3 +80,84 @@ subroutine bisection(lr,lt0,y,eps0,lt,bivar,keyerrt,keybisect)
 
 
 end subroutine bisection
+
+subroutine bisection_rho(lr0,lt,y,lpressin,lr,bivar,keyerrr,keybisect)
+
+  use eosmodule
+
+  implicit none
+
+  real*8 lr,lr0,y,lpressin,lt
+  integer keyerrr
+
+  integer keybisect !this doesn't do anything
+  ! 3 -> operate on pressure
+
+  !temporary vars
+  real*8 lr1,lr2,lrmin,lrmax
+  real*8 f1,f2,fmid,dlr,lrmid
+  real*8 d1,d2,d3,tol
+  real*8 f1a,f2a
+
+  integer bcount,i,itmax,maxbcount
+
+  real*8 bivar
+
+  bcount = 0
+  maxbcount = 80
+
+  tol=1.d-9 ! need to find energy to less than 1 in 10^-9
+  itmax=50
+
+  lrmax=logrho(nrho)
+  lrmin=logrho(1)
+
+  lr = lr0
+  lr1 = dlog10(min(10.0d0**lrmax,1.10d0*(10.0d0**lr0)))
+  lr2 = dlog10(max(10.0d0**lrmin,0.90d0*(10.0d0**lr0)))
+
+  call findthis(lr1,lt,y,f1a,bivar,d1,d2,d3)
+  call findthis(lr2,lt,y,f2a,bivar,d1,d2,d3)
+
+  f1=f1a-lpressin
+  f2=f2a-lpressin
+
+  keyerrr=0
+  
+  do while(f1*f2.ge.0.0d0)
+     bcount=bcount+1
+     lr1=dlog10(min(10.0d0**lrmax,1.2d0*(10.0d0**lr1)))
+     lr2=dlog10(max(10.0d0**lrmin,0.8d0*(10.0d0**lr2)))
+     call findthis(lr1,lt,y,f1a,bivar,d1,d2,d3)
+     call findthis(lr2,lt,y,f2a,bivar,d1,d2,d3)
+     f1=f1a-lpressin
+     f2=f2a-lpressin
+     if(bcount.ge.maxbcount) then
+        keyerrr=667
+        return
+     endif
+  enddo
+
+  if(f1.lt.0.0d0) then
+     lr=lr1
+     dlr=dlog10((10.0D0**lr2)-(10.0d0**lr1))
+  else
+     lr=lr2
+     dlr=dlog10((10.0D0**lr1)-(10.0d0**lr2))
+  endif
+
+  do i=1,itmax
+     dlr=dlog10((10.0d0**dlr)*0.5d0)
+     lrmid=dlog10(10.0d0**lr+10.0d0**dlr)
+     call findthis(lrmid,lt,y,f2a,bivar,d1,d2,d3)
+     fmid=f2a-lpressin
+     if(fmid.le.0.0d0) lr=lrmid
+     if(abs(1.0d0-f2a/lpressin).lt.tol) then
+        lr=lrmid
+        return
+     endif
+  enddo
+
+
+
+end subroutine bisection_rho
diff --git a/src/nuc_eos/make.code.defn b/src/nuc_eos/make.code.defn
index 59c8118..1b5de60 100644
--- a/src/nuc_eos/make.code.defn
+++ b/src/nuc_eos/make.code.defn
@@ -1,5 +1,5 @@
 SRCS = eosmodule.F90 nuc_eos.F90 bisection.F90 \
-       findtemp.F90 linterp_many.F90 readtable.F90 \
+       findtemp.F90 findrho.F90 linterp_many.F90 readtable.F90 \
        linterp.f
 
 SUBDIRS = 
diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90
index 399d89b..51170d0 100644
--- a/src/nuc_eos/nuc_eos.F90
+++ b/src/nuc_eos/nuc_eos.F90
@@ -24,7 +24,7 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
   use eosmodule
   implicit none
 
-  real*8, intent(in)    :: xrho,xye
+  real*8, intent(inout) :: xrho,xye
   real*8, intent(inout) :: xtemp,xenr,xent
   real*8, intent(out)   :: xprs,xcs2,xdedt
   real*8, intent(out)   :: xdpderho,xdpdrhoe,xxa,xxh,xxn,xxp
@@ -35,10 +35,11 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
 
 
   ! local variables
-  real*8 :: lr,lt,y,xx,xeps,leps,xs
+  real*8 :: lr,lt,y,xx,xeps,leps,xs,xpressure
   real*8 :: d1,d2,d3
   real*8 :: ff(nvars)
   integer :: keyerrt = 0
+  integer :: keyerrr = 0
 
   if(xrho.gt.eos_rhomax) then
      stop "nuc_eos: rho > rhomax"
@@ -49,11 +50,13 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
   endif
 
   if(xye.gt.eos_yemax) then
-     stop "nuc_eos: ye > yemax"
+     keyerr = 101
+     return
   endif
 
   if(xye.lt.eos_yemin) then
-     stop "nuc_eos: ye < yemin"
+     keyerr = 102
+     return
   endif
 
   if(keytemp.eq.1) then
@@ -87,19 +90,33 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
      xs = xent
      call findtemp_entropy(lr,lt,y,xs,keyerrt,rfeps)
      xtemp = 10.0d0**lt
+
+  elseif(keytemp.eq.3) then
+     !need to find rho based on xprs
+     xpressure = log10(xprs)
+     call findrho_press(lr,lt,y,xpressure,keyerrr,rfeps)
+     if(keyerrr.ne.0) then
+        write(*,*) "Problem in findrho_press:", keyerr
+        keyerr = keyerrr
+	return
+     endif
+     xrho = 10.0d0**lr
+
   endif
 
   ! have temperature, proceed:
   call findall(lr,lt,y,ff)
 
-  xprs = 10.0d0**ff(1)
+  !unless we want xprs to be constant (keytemp==3), reset xprs
+  if(.not.keytemp.eq.3) then
+     xprs = 10.0d0**ff(1)
+  endif
 
-  xeps = 10.0d0**ff(2) - energy_shift
-  if(keytemp.eq.1.or.keytemp.eq.2) then
-     xenr = xeps
+  if(.not.keytemp.eq.0) then
+     xenr = 10.0d0**ff(2) - energy_shift
   endif
 
-  if(keytemp.eq.1.or.keytemp.eq.0) then
+  if(.not.keytemp.eq.2) then
      xent = ff(3)
   endif
 
@@ -167,7 +184,7 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
   use eosmodule
   implicit none
 
-  real*8, intent(in)    :: xrho,xye
+  real*8, intent(inout) :: xrho,xye
   real*8, intent(inout) :: xtemp,xenr,xent
   real*8, intent(in)    :: rfeps
   real*8, intent(out)   :: xprs,xmunu,xcs2,xdedt
@@ -176,9 +193,10 @@ 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
+  real*8 :: lr,lt,y,xx,xeps,leps,xs,xpressure
   real*8 :: d1,d2,d3,ff(8)
   integer :: keyerrt = 0
+  integer :: keyerrr = 0
 
   if(xrho.gt.eos_rhomax) then
      stop "nuc_eos: rho > rhomax"
@@ -191,13 +209,13 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
   endif
 
   if(xye.gt.eos_yemax) then
-     stop "nuc_eos: ye > yemax"
+     keyerr = 101
+     return
   endif
 
   if(xye.lt.eos_yemin) then
      keyerr = 102
-     write(6,*) "ye: ", xye
-     stop "nuc_eos ye < yemin"
+     return
   endif
 
   if(keytemp.eq.1) then
@@ -220,7 +238,7 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
 
   keyerr = 0
 
-  if(keytemp.eq.0) then
+if(keytemp.eq.0) then
      !need to find temperature based on xeps
      call findtemp(lr,lt,y,leps,keyerrt,rfeps)
      if(keyerrt.ne.0) then
@@ -239,18 +257,33 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
      endif
      xtemp = 10.0d0**lt
 
+  elseif(keytemp.eq.3) then
+     !need to find rho based on xprs
+     xpressure = log10(xprs)
+     call findrho_press(lr,lt,y,xpressure,keyerrr,rfeps)
+     if (keyerrr.ne.0) then
+        keyerr = keyerrr
+        write(*,*) "Problem in findrho_press:", keyerr
+        return
+     endif
+     xrho = 10.0d0**lr
   endif
 
-  ! have temperature, proceed:
+  ! have rho,temp,ye; proceed:
   call findall_short(lr,lt,y,ff)
-  xprs = 10.0d0**ff(1)
 
-  xeps = 10.0d0**ff(2) - energy_shift
-  if((keytemp.eq.1).or.(keytemp.eq.2)) then
-     xenr = xeps
+  !unless we want xprs to be constant (keytemp==3), reset xprs
+  if(.not.keytemp.eq.3) then
+     xprs = 10.0d0**ff(1)
   endif
 
-  if(keytemp.eq.1.or.keytemp.eq.0) then
+  !unless we want xenr to be constant (keytemp==0), reset xenr
+  if(.not.keytemp.eq.0) then
+     xenr = 10.0d0**ff(2) - energy_shift
+  endif
+
+  !unless we want xent to be constant (keytemp==2), reset xent
+  if(.not.keytemp.eq.2) then
      xent = ff(3)
   endif
 
@@ -296,13 +329,13 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,&
   endif
 
   if(xye.gt.eos_yemax) then
-     stop "nuc_eos: ye > yemax"
+     keyerr = 101
+     return
   endif
 
   if(xye.lt.eos_yemin) then
      keyerr = 102
-     write(6,*) "ye: ", xye
-     stop "nuc_eos ye < yemin"
+     return
   endif
 
   if(keytemp.eq.1) then
@@ -333,17 +366,16 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,&
      endif
      xtemp = 10.0d0**lt
 
-  elseif(keytemp.eq.2) then
-     stop "eos_nuc_press does not support keytemp.eq.2"
+  elseif(keytemp.gt.1) then
+     stop "eos_nuc_press does not support keytemp other than 0 and 1"
   endif
 
   ! have temperature, proceed:
   call findall_press_eps(lr,lt,y,ff)
   xprs = 10.0d0**ff(1)
 
-  xeps = 10.0d0**ff(2) - energy_shift
-  if((keytemp.eq.1).or.(keytemp.eq.2)) then
-     xenr = xeps
+  if(keytemp.eq.1) then
+     xenr = 10.0d0**ff(2) - energy_shift
   endif
 
 end subroutine nuc_eos_press_eps
@@ -378,13 +410,13 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,&
   endif
 
   if(xye.gt.eos_yemax) then
-     stop "nuc_eos: ye > yemax"
+     keyerr = 101
+     return
   endif
 
   if(xye.lt.eos_yemin) then
      keyerr = 102
-     write(6,*) "ye: ", xye
-     stop "nuc_eos ye < yemin"
+     return
   endif
 
   if(keytemp.eq.1) then
@@ -415,8 +447,8 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,&
      endif
      xtemp = 10.0d0**lt
 
-  elseif(keytemp.eq.2) then
-     stop "eos_nuc_press does not support keytemp.eq.2"
+  elseif(keytemp.gt.1) then
+     stop "eos_nuc_press does not support keytemp > 1"
   endif
 
   ! have temperature, proceed:
@@ -424,8 +456,7 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,&
   xdpdrhoe = 10.0d0**ff(1)
   xdpderho = 10.0d0**ff(2)
 
-
-  if((keytemp.eq.1).or.(keytemp.eq.2)) then
+  if(keytemp.eq.1) then
      call findthis(lr,lt,y,xeps,alltables(:,:,:,2),d1,d2,d3)
      xeps = 10.0d0**xeps - energy_shift
      xenr = xeps