* updates needed to make TOVSolverHot work

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

1 files changed, 67 insertions, 36 deletions

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