* fix a number of small issues

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

4 files changed, 140 insertions, 14 deletions

diff --git a/src/EOS_Omni_Handles.c b/src/EOS_Omni_Handles.c
index 41dbced..b744cac 100644
--- a/src/EOS_Omni_Handles.c
+++ b/src/EOS_Omni_Handles.c
@@ -10,6 +10,8 @@ CCTK_INT EOS_Omni_GetHandle_(CCTK_STRING name)
         return 2;
     if (CCTK_EQUALS(name, "Hybrid"))
         return 3;
+    if (CCTK_EQUALS(name, "nuc_eos"))
+        return 4;
     return 0;
 }
 
diff --git a/src/EOS_Omni_Module.F90 b/src/EOS_Omni_Module.F90
index 75c6c2d..f98930d 100644
--- a/src/EOS_Omni_Module.F90
+++ b/src/EOS_Omni_Module.F90
@@ -3,23 +3,22 @@ module EOS_Omni_Module
   implicit none
 
   real*8,parameter :: rho_gf = 1.61620075314614d-18
-!  real*8,parameter :: rho_gf = 1.61930347d-18
   real*8,parameter :: press_gf = 1.7982953469278d-39
-!  real*8,parameter :: press_gf = 1.80171810d-39
   real*8,parameter :: eps_gf = 1.11265006d-21
   real*8,parameter :: time_gf = 2.03001708d+05
   real*8,parameter :: mass_gf = 5.02765209d-34
   real*8,parameter :: length_gf = 6.77140812d-06
   
   real*8,parameter :: inv_rho_gf = 6.18735016707159d17
-!  real*8,parameter :: inv_rho_gf = 6.17549470205236d17
   real*8,parameter :: inv_press_gf = 5.56082181777535d38
-!  real*8,parameter :: inv_press_gf = 5.55025783445257d38
   real*8,parameter :: inv_eps_gf = 8.98755175549085d20
   real*8,parameter :: inv_time_gf = 4.92606692747629d-6
   real*8,parameter :: inv_mass_gf = 1.98899999860571d33
   real*8,parameter :: inv_length_gf = 147679.77092481d0
 
+  real*8,parameter :: clite = 2.99792458d10
+  real*8,parameter :: cliteinv2 = 1.11265005605362d-21
+
   real*8 :: poly_k_cgs = 0.0d0
   real*8 :: gl_k_cgs   = 0.0d0
 
diff --git a/src/EOS_Omni_SingleVarCalls.F90 b/src/EOS_Omni_SingleVarCalls.F90
index 217c8f3..ce8a22e 100644
--- a/src/EOS_Omni_SingleVarCalls.F90
+++ b/src/EOS_Omni_SingleVarCalls.F90
@@ -30,6 +30,10 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
   real*8           :: hybrid_local_gamma, hybrid_local_k_cgs, &
                       hybrid_p_poly, hybrid_p_th
   real*8,parameter :: zero = 0.0d0
+  ! temporary vars for nuc_eos
+  real*8           :: xrho,xye,xtemp,xenr,xent
+  real*8           :: xprs,xmunu,xcs2
+  real*8           :: xdedt,xdpderho,xdpdrhoe
 
   anyerr    = 0
   keyerr(:) = 0
@@ -84,12 +88,31 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
            press(i) = hybrid_p_poly + hybrid_p_th
         enddo
 
-!!!     case (4)
-!!!
-!!!        do i=1,npoints
-!!!
-!!!        enddo
+     case (4)
 
+        do i=1,npoints
+
+           xrho = rho(i) * inv_rho_gf
+           xtemp = temp(i)
+           xye = ye(i)
+           xenr = eps(i) * inv_eps_gf
+           call nuc_eos_short(xrho,xtemp,xye,xenr,xprs,& 
+                xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,&
+                keytemp,keyerr(i),rf_precision)
+           
+           if(keyerr(i).ne.0) then
+              anyerr = 1
+           endif
+
+           if(keytemp.eq.1) then
+              eps(i) = xenr * eps_gf
+           else
+              temp(i) = xtemp
+           endif
+
+           press(i) = xprs * press_gf
+
+        enddo
 
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
@@ -121,6 +144,10 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
   real*8           :: hybrid_local_gamma, hybrid_local_k_cgs, &
                       hybrid_dp_poly, hybrid_dp_th1, hybrid_dp_th2
   real*8,parameter :: zero = 0.0d0
+  ! temporary vars for nuc_eos
+  real*8           :: xrho,xye,xtemp,xenr,xent
+  real*8           :: xprs,xmunu,xcs2
+  real*8           :: xdedt,xdpderho,xdpdrhoe
 
   anyerr    = 0
   keyerr(:) = 0
@@ -181,6 +208,28 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
 
         enddo
 
+     case (4)
+        do i=1,npoints
+           xrho = rho(i) * inv_rho_gf
+           xtemp = temp(i)
+           xye = ye(i)
+           xenr = eps(i) * inv_eps_gf
+           call nuc_eos_short(xrho,xtemp,xye,xenr,xprs,& 
+                xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,&
+                keytemp,keyerr(i),rf_precision)
+           
+           if(keyerr(i).ne.0) then
+              anyerr = 1
+           endif
+
+           if(keytemp.eq.1) then
+              eps(i) = xenr * eps_gf
+           else
+              temp(i) = xtemp
+           endif
+
+           dpdrhoe(i) = xdpdrhoe * press_gf * inv_rho_gf
+        enddo
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
         call CCTK_WARN(0,warnstring)
@@ -206,6 +255,11 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,rf_precision,npoints,&
   ! local vars
   integer          :: i
   character(256)   :: warnstring
+  ! temporary vars for nuc_eos
+  real*8           :: xrho,xye,xtemp,xenr,xent
+  real*8           :: xprs,xmunu,xcs2
+  real*8           :: xdedt,xdpderho,xdpdrhoe
+
   anyerr    = 0
   keyerr(:) = 0
 
@@ -241,6 +295,32 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,rf_precision,npoints,&
            dpdepsrho(i) = (hybrid_gamma_th - 1.0d0) * rho(i)
         enddo
 
+     case (4)
+        ! nuc_eos
+        do i=1,npoints
+
+           xrho = rho(i) * inv_rho_gf
+           xtemp = temp(i)
+           xye = ye(i)
+           xenr = eps(i) * inv_eps_gf
+           call nuc_eos_short(xrho,xtemp,xye,xenr,xprs,& 
+                xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,&
+                keytemp,keyerr(i),rf_precision)
+           
+           if(keyerr(i).ne.0) then
+              anyerr = 1
+           endif
+
+           if(keytemp.eq.1) then
+              eps(i) = xenr * eps_gf
+           else
+              temp(i) = xtemp
+           endif
+
+           dpdepsrho(i) = xdpderho * press_gf * inv_eps_gf 
+
+        enddo
+
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
         call CCTK_WARN(0,warnstring)
@@ -270,6 +350,10 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
   real*8           :: hybrid_local_gamma, hybrid_local_k_cgs, &
                       hybrid_p_poly, hybrid_p_th
   real*8,parameter :: zero = 0.0d0
+  ! temporary vars for nuc_eos
+  real*8           :: xrho,xye,xtemp,xenr,xent
+  real*8           :: xprs,xmunu,xcs2
+  real*8           :: xdedt
 
   anyerr    = 0
   keyerr(:) = 0
@@ -331,6 +415,35 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
            cs2(i) = (hybrid_local_gamma * hybrid_p_poly + hybrid_gamma_th * hybrid_p_th) / &
                 rho(i) / (1.0d0 + eps(i) + xpress/rho(i))
         enddo
+     case(4)
+        ! nuc_eos
+
+        do i=1,npoints
+           
+           xrho = rho(i) * inv_rho_gf
+           xtemp = temp(i)
+           xye = ye(i)
+           xenr = eps(i) * inv_eps_gf
+           call nuc_eos_short(xrho,xtemp,xye,xenr,xprs,& 
+                xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,&
+                keytemp,keyerr(i),rf_precision)
+           
+           if(keyerr(i).ne.0) then
+              anyerr = 1
+           endif
+
+           if(keytemp.eq.1) then
+              eps(i) = xenr * eps_gf
+           else
+              temp(i) = xtemp
+           endif
+
+           cs2(i) = xcs2 * cliteinv2 / &
+                (1.0d0 + eps(i) + xprs * press_gf / rho(i))
+
+        enddo
+
+
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
         call CCTK_WARN(0,warnstring)
@@ -357,6 +470,10 @@ subroutine EOS_Omni_EOS_eps_from_press(eoskey,keytemp,rf_precision,npoints,&
   ! local vars
   integer        :: i
   character(256) :: warnstring
+  ! temporary vars for nuc_eos
+  real*8           :: xrho,xye,xtemp,xenr,xent
+  real*8           :: xprs,xmunu,xcs2
+  real*8           :: xdedt,xdpderho,xdpdrhoe
 
   if(keytemp.eq.1) then
      anyerr = 1
@@ -381,6 +498,10 @@ subroutine EOS_Omni_EOS_eps_from_press(eoskey,keytemp,rf_precision,npoints,&
         ! hybrid EOS
         write(warnstring,*) "EOS_Omni_EpsFromPress call not supported for hybrid EOS"
         call CCTK_WARN(0,warnstring)
+     case (4)
+        ! nuc EOS
+        write(warnstring,*) "EOS_Omni_EpsFromPress call not supported for nuc_eos yet"
+        call CCTK_WARN(0,warnstring)
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
         call CCTK_WARN(0,warnstring)
@@ -431,6 +552,10 @@ subroutine EOS_Omni_EOS_RestMassDensityFromEpsPress(eoskey,keytemp,rf_precision,
         ! hybrid EOS
         write(warnstring,*) "EOS_Omni_RestMassDensityFromEpsPress not supported for hybrid EOS"
         call CCTK_WARN(0,warnstring)
+     case (4)
+        ! nuc EOS
+        write(warnstring,*) "EOS_Omni_RestMassDensityFromEpsPress call not supported for nuc_eos yet"
+        call CCTK_WARN(0,warnstring)
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
         call CCTK_WARN(0,warnstring)
diff --git a/src/nuc_eos/readtable.F90 b/src/nuc_eos/readtable.F90
index d8ae2bc..1c3e7e2 100644
--- a/src/nuc_eos/readtable.F90
+++ b/src/nuc_eos/readtable.F90
@@ -21,13 +21,13 @@ subroutine nuc_eos_readtable(eos_filename)
   real*8 buffer1,buffer2,buffer3,buffer4
   accerr=0
 
-  write(*,*) "Reading Nuclear EOS Table"
+!  write(*,*) "Reading Nuclear EOS Table"
 
   call h5open_f(error)
 
   call h5fopen_f (trim(adjustl(eos_filename)), H5F_ACC_RDONLY_F, file_id, error)
 
-  write(6,*) trim(adjustl(eos_filename))
+!  write(6,*) trim(adjustl(eos_filename))
 
 ! read scalars
   dims1(1)=1
@@ -57,8 +57,8 @@ subroutine nuc_eos_readtable(eos_filename)
      stop "Could not read EOS table file"
   endif
 
-  write(message,"(a25,1P3i5)") "We have nrho ntemp nye: ", nrho,ntemp,nye
-  write(*,*) message
+!  write(message,"(a25,1P3i5)") "We have nrho ntemp nye: ", nrho,ntemp,nye
+!  write(*,*) message
 
   allocate(alltables(nrho,ntemp,nye,nvars))
 
@@ -226,7 +226,7 @@ subroutine nuc_eos_readtable(eos_filename)
   eos_tempmin = 10.0d0**logtemp(1)
   eos_tempmax = 10.0d0**logtemp(ntemp)
 
-  write(6,*) "Done reading eos tables"
+!  write(6,*) "Done reading eos tables"
 
 
 end subroutine nuc_eos_readtable