/*@@
   @file      EOS_Omni_HelmholtzEOS_NSE.F90
   @date      Feb 12, 2016
   @author    Daniel Siegel
   @desc 
   implementation of the Helmholtz EOS
   @enddesc 
 @@*/

#include "HelmholtzEOS_NSE.h"

#if !STANDALONE
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "cctk_Functions.h"
#endif

#define inv_tol 1d-12


!--------------------------------------------------------------------------
!
! Wrapper routines for Helmholtz EOS
!
!--------------------------------------------------------------------------


subroutine helm_eos_kt1_press_eps(npoints, rho, temp, ye, &
       eps, press, NSEkeyerr, anyerr)

  implicit none

#if !STANDALONE
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, temp, ye
  CCTK_REAL, dimension(npoints), intent(OUT) :: press, eps
  CCTK_INT, intent(out) :: anyerr, NSEkeyerr(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, temp, ye
  real*8, dimension(npoints), intent(OUT) :: press, eps
  integer, intent(out) :: anyerr, NSEkeyerr(npoints)
#endif
  
  real*8, dimension(npoints) :: stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow

  call helm_eos(npoints, rho, temp, ye, &
       eps, press, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, NSEkeyerr, anyerr)

  return

end subroutine helm_eos_kt1_press_eps


subroutine helm_eos_kt1_press_eps_point(rho, temp, ye, eps, press, anyerr)

  implicit none

#if !STANDALONE
  CCTK_REAL, intent(IN)  :: rho, temp, ye
  CCTK_REAL, intent(OUT) :: press, eps
  CCTK_INT, intent(OUT)  :: anyerr
#else
  real*8, intent(IN)   :: rho, temp, ye
  real*8, intent(OUT)  :: press, eps
  integer, intent(OUT) :: anyerr
#endif
  
  real*8, dimension(1) :: rhoRow, tempRow, yeRow, epsRow, pressRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow
  integer, dimension(1) :: NSEkeyerr

  rhoRow(1) = rho
  tempRow(1) = temp
  yeRow(1) = ye

  call helm_eos(1, rhoRow, tempRow, yeRow, &
       epsRow, pressRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, NSEkeyerr, anyerr)

  eps = epsRow(1)
  press = pressRow(1)

  return

end subroutine helm_eos_kt1_press_eps_point


subroutine helm_eos_kt0_press(npoints, rho, temp, ye, &
             eps, press, inv_err, anyerr)

  implicit none

#if !STANDALONE
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, ye, eps
  CCTK_REAL, dimension(npoints), intent(OUT) :: press
  CCTK_REAL, dimension(npoints), intent(INOUT) :: temp
  CCTK_INT, intent(out) :: anyerr, inv_err(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, ye, eps
  real*8, dimension(npoints), intent(OUT) :: press
  real*8, dimension(npoints), intent(INOUT) :: temp
  integer, intent(out) :: anyerr, inv_err(npoints)
#endif
  
  ! local vars
  real*8, dimension(npoints) :: stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, tempguess
  integer, dimension(npoints) :: xnEOScalls, xnEOScalls_bisec     
  integer :: j   
  
  tempguess = temp  
  
  call HelmEOS_invert_et(npoints, rho, temp, ye, press, eps, &
          stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
          pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
          xxnRow, xxpRow, xxaRow, xxatmRow, inv_err, anyerr, xnEOScalls)

  ! if Newton-Raphson inversion fails or does not converge fast enough
  ! try bisection
  if (anyerr .eq. 1) then
    do j=1, npoints
      if (inv_err(j) .ne. 0) then
        call HelmEOS_invert_et_bisec(1, rho(j), tempguess(j), ye(j), press(j), eps(j), &
               stotRow(j), htotRow(j), AbarRow(j), dpdRow(j), dptRow(j), dedRow(j), &
               detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
               pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), &
               cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
               xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), inv_err(j), anyerr, xnEOScalls_bisec(j))   
        temp(j) = tempguess(j)
      end if
    end do
  end if
    
  return

end subroutine helm_eos_kt0_press


subroutine helm_eos_kt0_press_temp_point(rho, temp, ye, eps, press, anyerr, nEOScalls)

  implicit none

#if !STANDALONE
  CCTK_REAL, intent(IN)  :: rho, ye, eps
  CCTK_REAL, intent(OUT) :: press
  CCTK_REAL, intent(INOUT) :: temp
  CCTK_INT, intent(OUT) :: anyerr
  CCTK_INT, intent(INOUT) :: nEOScalls
#else
  real*8, intent(IN)  :: rho, ye, eps
  real*8, intent(OUT) :: press
  real*8, intent(INOUT) :: temp
  integer, intent(OUT) :: anyerr
  integer, intent(INOUT) :: nEOScalls
#endif
  
  ! local vars
  real*8, dimension(1) :: rhoRow, tempRow, yeRow, epsRow, pressRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, tempguess
  integer, dimension(1) :: key_err
  integer, dimension(1) :: xnEOScalls, xnEOScalls_bisec
  
  tempguess(1) = temp  
  
  rhoRow(1) = rho
  tempRow(1) = temp
  yeRow(1) = ye
  epsRow(1) = eps
  key_err(1) = 0
  
  call HelmEOS_invert_et(1, rhoRow, tempRow, yeRow, pressRow, epsRow, &
          stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
          pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
          xxnRow, xxpRow, xxaRow, xxatmRow, key_err, anyerr, xnEOScalls)

  nEOScalls = xnEOScalls(1)

  ! if Newton-Raphson inversion fails or does not converge fast enough
  ! try bisection
  if (anyerr .eq. 1) then
      if (key_err(1) .ne. 0) then
        call HelmEOS_invert_et_bisec(1, rhoRow, tempguess, yeRow, pressRow, epsRow, &
               stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow, detRow, &
               dsdRow, dstRow, dhdRow, dhtRow, &
               pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, &
               cvRow, cpRow, cs2Row, dpde_d_Row, &
               xxnRow, xxpRow, xxaRow, xxatmRow, key_err, anyerr, xnEOScalls_bisec)   
        tempRow(1) = tempguess(1)
        nEOScalls = nEOScalls + xnEOScalls_bisec(1)
      end if
  end if
  
  temp = tempRow(1)
  press = pressRow(1)
    
  return

end subroutine helm_eos_kt0_press_temp_point


subroutine helm_eos_kt1_press_eps_ent_abar(npoints, rho, temp, ye, &
       eps, press, ent, abar, NSEkeyerr, anyerr)

  implicit none

#if !STANDALONE
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, temp, ye
  CCTK_REAL, dimension(npoints), intent(OUT) :: press, eps, ent, abar
  CCTK_INT, intent(out) :: anyerr, NSEkeyerr(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, temp, ye
  real*8, dimension(npoints), intent(OUT) :: press, eps, ent, abar
  integer, intent(out) :: anyerr, NSEkeyerr(npoints)
#endif
  
  real*8, dimension(npoints) :: htotRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow

  call helm_eos(npoints, rho, temp, ye, &
       eps, press, ent, htotRow, abar, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, NSEkeyerr, anyerr)

  return

end subroutine helm_eos_kt1_press_eps_ent_abar


subroutine helm_eos_kt0_press_ent_abar(npoints, rho, temp, ye, &
             eps, press, ent, abar, inv_err, anyerr)

  implicit none

#if !STANDALONE
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, ye, eps
  CCTK_REAL, dimension(npoints), intent(OUT) :: press, ent, abar
  CCTK_REAL, dimension(npoints), intent(INOUT) :: temp
  CCTK_INT, intent(out) :: anyerr, inv_err(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, ye, eps
  real*8, dimension(npoints), intent(OUT) :: press, ent, abar
  real*8, dimension(npoints), intent(INOUT) :: temp
  integer, intent(out) :: anyerr, inv_err(npoints)
#endif
  
  ! local vars
  real*8, dimension(npoints) :: htotRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, tempguess
  integer, dimension(npoints) :: xnEOScalls, xnEOScalls_bisec     
  integer :: j   
  
  tempguess = temp  
  
  call HelmEOS_invert_et(npoints, rho, temp, ye, press, eps, &
          ent, htotRow, abar, dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
          pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
          xxnRow, xxpRow, xxaRow, xxatmRow, inv_err, anyerr, xnEOScalls)

  ! if Newton-Raphson inversion fails or does not converge fast enough
  ! try bisection
  if (anyerr .eq. 1) then
    do j=1, npoints
      if (inv_err(j) .ne. 0) then
        call HelmEOS_invert_et_bisec(1, rho(j), tempguess(j), ye(j), press(j), eps(j), &
               ent(j), htotRow(j), abar(j), dpdRow(j), dptRow(j), dedRow(j), detRow(j), &
               dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
               pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), &
               cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
               xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), inv_err(j), anyerr, xnEOScalls_bisec(j))   
        temp(j) = tempguess(j)
      end if
    end do
  end if
    
  return

end subroutine helm_eos_kt0_press_ent_abar


subroutine helm_eos_kt1_press_eps_ent_abar_point(rho, temp, ye, eps, press, ent, abar, anyerr)

  implicit none

#if !STANDALONE
  CCTK_REAL, intent(IN)  :: rho, temp, ye
  CCTK_REAL, intent(OUT) :: press, eps, ent, abar
  CCTK_INT, intent(OUT)  :: anyerr
#else
  real*8, intent(IN)   :: rho, temp, ye
  real*8, intent(OUT)  :: press, eps, ent, abar
  integer, intent(OUT) :: anyerr
#endif
  
  real*8, dimension(1) :: rhoRow, tempRow, yeRow, epsRow, pressRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow
  integer, dimension(1) :: NSEkeyerr

  rhoRow(1) = rho
  tempRow(1) = temp
  yeRow(1) = ye

  call helm_eos(1, rhoRow, tempRow, yeRow, &
       epsRow, pressRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, NSEkeyerr, anyerr)

  eps = epsRow(1)
  press = pressRow(1)
  ent = stotRow(1)
  abar = AbarRow(1)

  return

end subroutine helm_eos_kt1_press_eps_ent_abar_point


subroutine helm_eos_kt0_press_temp_ent_abar_point(rho, temp, ye, eps, press, ent, abar, anyerr, nEOScalls)

  implicit none

#if !STANDALONE
  CCTK_REAL, intent(IN)  :: rho, ye, eps
  CCTK_REAL, intent(OUT) :: press, ent, abar
  CCTK_REAL, intent(INOUT) :: temp
  CCTK_INT, intent(OUT) :: anyerr
  CCTK_INT, intent(INOUT) :: nEOScalls
#else
  real*8, intent(IN)  :: rho, ye, eps
  real*8, intent(OUT) :: press, ent, abar
  real*8, intent(INOUT) :: temp
  integer, intent(OUT) :: anyerr
  integer, intent(INOUT) :: nEOScalls
#endif
  
  ! local vars
  real*8, dimension(1) :: rhoRow, tempRow, yeRow, epsRow, pressRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, tempguess
  integer, dimension(1) :: key_err
  integer, dimension(1) :: xnEOScalls, xnEOScalls_bisec
  
  tempguess(1) = temp  
  
  rhoRow(1) = rho
  tempRow(1) = temp
  yeRow(1) = ye
  epsRow(1) = eps
  key_err(1) = 0
  
  call HelmEOS_invert_et(1, rhoRow, tempRow, yeRow, pressRow, epsRow, &
          stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
          pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
          xxnRow, xxpRow, xxaRow, xxatmRow, key_err, anyerr, xnEOScalls)

  nEOScalls = xnEOScalls(1)

  ! if Newton-Raphson inversion fails or does not converge fast enough
  ! try bisection
  if (anyerr .eq. 1) then
      if (key_err(1) .ne. 0) then
        call HelmEOS_invert_et_bisec(1, rhoRow, tempguess, yeRow, pressRow, epsRow, &
               stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow, &
               detRow, dsdRow, dstRow, dhdRow, dhtRow, &
               pelRow, neRow, etaeRow, etapRow, etanRow, &
               gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
               xxnRow, xxpRow, xxaRow, xxatmRow, key_err, anyerr, xnEOScalls_bisec)   
        tempRow(1) = tempguess(1)
        nEOScalls = nEOScalls + xnEOScalls_bisec(1)
      end if
  end if
  
  temp = tempRow(1)
  press = pressRow(1)
  ent = stotRow(1)
  abar = AbarRow(1)
    
  return

end subroutine helm_eos_kt0_press_temp_ent_abar_point


subroutine helm_eos_m_kt1_eps_cs2(npoints, rho, temp, ye, eps, cs2, NSEkeyerr, anyerr)

  implicit none

#if !STANDALONE	
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, temp, ye
  CCTK_REAL, dimension(npoints), intent(OUT) :: eps, cs2
  CCTK_INT, intent(out) :: anyerr, NSEkeyerr(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, temp, ye
  real*8, dimension(npoints), intent(OUT) :: eps, cs2
  integer, intent(out) :: anyerr, NSEkeyerr(npoints)
#endif

  ! local vars  
  real*8, dimension(npoints) :: ptotRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow

  call helm_eos(npoints, rho, temp, ye, &
       eps, ptotRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, NSEkeyerr, anyerr)   
       
  return     

end subroutine helm_eos_m_kt1_eps_cs2


subroutine helm_eos_m_kt1_eps_press_cs2(npoints, rho, temp, ye, eps, press, cs2, NSEkeyerr, anyerr)

  implicit none

#if !STANDALONE	
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, temp, ye
  CCTK_REAL, dimension(npoints), intent(OUT) :: eps, press, cs2
  CCTK_INT, intent(out) :: anyerr, NSEkeyerr(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, temp, ye
  real*8, dimension(npoints), intent(OUT) :: eps, press, cs2
  integer, intent(out) :: anyerr, NSEkeyerr(npoints)
#endif
    
  ! local vars  
  real*8, dimension(npoints) :: stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow

  call helm_eos(npoints, rho, temp, ye, &
       eps, press, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, NSEkeyerr, anyerr)   
       
  return     

end subroutine helm_eos_m_kt1_eps_press_cs2


subroutine helm_eos_m_kt0_cs2(npoints, rho, temp, ye, eps, cs2, inv_err, anyerr)

  implicit none

#if !STANDALONE	
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)    :: rho, ye, eps
  CCTK_REAL, dimension(npoints), intent(INOUT) :: temp
  CCTK_REAL, dimension(npoints), intent(OUT)   :: cs2 !eps, cs2
  CCTK_INT, intent(out) :: anyerr, inv_err(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)    :: rho, ye, eps
  real*8, dimension(npoints), intent(INOUT) :: temp
  real*8, dimension(npoints), intent(OUT)   :: cs2 !eps, cs2
  integer, intent(out) :: anyerr, inv_err(npoints)
#endif

  ! local vars  
  real*8, dimension(npoints) :: ptotRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, tempguess
  integer, dimension(npoints) :: xnEOScalls, xnEOScalls_bisec
  integer :: j  
  
  tempguess = temp
  
  call HelmEOS_invert_et(npoints, rho, temp, ye, ptotRow, eps, &
           stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
           pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2, dpde_d_Row, &
           xxnRow, xxpRow, xxaRow, xxatmRow, inv_err, anyerr, xnEOScalls)
  
  ! if Newton-Raphson inversion fails or does not converge fast enough
  ! try bisection
  if (anyerr .eq. 1) then
    do j=1, npoints
      if (inv_err(j) .ne. 0) then
        call HelmEOS_invert_et_bisec(1, rho(j), tempguess(j), ye(j), ptotRow(j), eps(j), &
               stotRow(j), htotRow(j), AbarRow(j), dpdRow(j), dptRow(j), dedRow(j), &
               detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
               pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), &
               cvRow(j), cpRow(j), cs2(j), dpde_d_Row(j), &
               xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), inv_err(j), anyerr, xnEOScalls_bisec(j))   
        temp(j) = tempguess(j)
      end if
    end do
  end if  
       
  return     

end subroutine helm_eos_m_kt0_cs2


subroutine helm_eos_m_kt0_press_cs2(npoints, rho, temp, ye, eps, press, cs2, inv_err, anyerr)

  implicit none

#if !STANDALONE	
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)    :: rho, ye, eps
  CCTK_REAL, dimension(npoints), intent(INOUT) :: temp
  CCTK_REAL, dimension(npoints), intent(OUT)   :: press, cs2 !eps, cs2
  CCTK_INT, intent(out) :: anyerr, inv_err(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)    :: rho, ye, eps
  real*8, dimension(npoints), intent(INOUT) :: temp
  real*8, dimension(npoints), intent(OUT)   :: press, cs2 !eps, cs2
  integer, intent(out) :: anyerr, inv_err(npoints)
#endif
    
  ! local vars  
  real*8, dimension(npoints) :: stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, tempguess
  integer, dimension(npoints) :: xnEOScalls, xnEOScalls_bisec
  integer :: j  
  
  tempguess = temp
  
  call HelmEOS_invert_et(npoints, rho, temp, ye, press, eps, &
           stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
           pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2, dpde_d_Row, &
           xxnRow, xxpRow, xxaRow, xxatmRow, inv_err, anyerr, xnEOScalls)
  
  ! if Newton-Raphson inversion fails or does not converge fast enough
  ! try bisection
  if (anyerr .eq. 1) then
    do j=1, npoints
      if (inv_err(j) .ne. 0) then
        call HelmEOS_invert_et_bisec(1, rho(j), tempguess(j), ye(j), press(j), eps(j), &
               stotRow(j), htotRow(j), AbarRow(j), dpdRow(j), dptRow(j), dedRow(j), &
               detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
               pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), &
               gam1Row(j), cvRow(j), cpRow(j), cs2(j), dpde_d_Row(j), &
               xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), inv_err(j), anyerr, xnEOScalls_bisec(j))   
        temp(j) = tempguess(j)
      end if
    end do
  end if  
       
  return     

end subroutine helm_eos_m_kt0_press_cs2


subroutine helm_eos_m_kt1_dpressdeps_eps(npoints, rho, temp, ye, eps, dpdeps_rho, NSEkeyerr, anyerr)

  implicit none

#if !STANDALONE	
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, temp, ye
  CCTK_REAL, dimension(npoints), intent(OUT) :: dpdeps_rho, eps
  CCTK_INT, intent(out) :: anyerr, NSEkeyerr(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, temp, ye
  real*8, dimension(npoints), intent(OUT) :: dpdeps_rho, eps
  integer, intent(out) :: anyerr, NSEkeyerr(npoints)
#endif
    
  ! local vars  
  real*8, dimension(npoints) :: ptotRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow

  call helm_eos(npoints, rho, temp, ye, &
       eps, ptotRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpdeps_rho, &
       xxnRow, xxpRow, xxaRow, xxatmRow, NSEkeyerr, anyerr)
       
  return     

end subroutine helm_eos_m_kt1_dpressdeps_eps


subroutine helm_eos_m_kt0_dpressdeps(npoints, rho, temp, ye, eps, dpdeps_rho, &
               inv_err, anyerr)

  implicit none

#if !STANDALONE	
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, ye, eps
  CCTK_REAL, dimension(npoints), intent(OUT) :: dpdeps_rho !, eps
  CCTK_REAL, dimension(npoints), intent(INOUT) :: temp
  CCTK_INT, intent(OUT) :: anyerr, inv_err(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, ye, eps
  real*8, dimension(npoints), intent(OUT) :: dpdeps_rho !, eps
  real*8, dimension(npoints), intent(INOUT) :: temp
  integer, intent(OUT) :: anyerr, inv_err(npoints)
#endif
    
  ! local vars  
  real*8, dimension(npoints) :: ptotRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, tempguess
  integer :: xnEOScalls(npoints), xnEOScalls_bisec(npoints)
  integer :: j  
  
  tempguess = temp
  
  call HelmEOS_invert_et(npoints, rho, temp, ye, ptotRow, eps, &
           stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
           pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpdeps_rho, &
           xxnRow, xxpRow, xxaRow, xxatmRow, inv_err, anyerr, xnEOScalls)

  ! if Newton-Raphson inversion fails or does not converge fast enough
  ! try bisection
  if (anyerr .eq. 1) then
    do j=1, npoints
      if (inv_err(j) .ne. 0) then
        call HelmEOS_invert_et_bisec(1, rho(j), tempguess(j), ye(j), ptotRow(j), eps(j), &
               stotRow(j), htotRow(j), AbarRow(j), dpdRow(j), dptRow(j), dedRow(j), &
               detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
               pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), &
               gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpdeps_rho(j), &
               xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), inv_err(j), anyerr, xnEOScalls_bisec(j))   
        temp(j) = tempguess(j)
      end if
    end do
  end if  

  return

end subroutine helm_eos_m_kt0_dpressdeps


subroutine helm_eos_m_kt1_dpressdrho_eps(npoints, rho, temp, ye, eps, dpdrho, &
               NSEkeyerr, anyerr)

  implicit none

#if !STANDALONE	
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, temp, ye
  CCTK_REAL, dimension(npoints), intent(OUT) :: dpdrho, eps
  CCTK_INT, intent(out) :: anyerr, NSEkeyerr(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, temp, ye
  real*8, dimension(npoints), intent(OUT) :: dpdrho, eps
  integer, intent(out) :: anyerr, NSEkeyerr(npoints)
#endif
    
  ! local vars  
  real*8, dimension(npoints) :: ptotRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow

  call helm_eos(npoints, rho, temp, ye, &
       eps, ptotRow, stotRow, htotRow, AbarRow, &
       dpdrho, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, NSEkeyerr, anyerr)
       
  return     

end subroutine helm_eos_m_kt1_dpressdrho_eps


subroutine helm_eos_m_kt0_dpressdrho(npoints, rho, temp, ye, eps, dpdrho, &
               inv_err, anyerr)

  implicit none
  
#if !STANDALONE	
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, ye, eps
  CCTK_REAL, dimension(npoints), intent(OUT) :: dpdrho !, eps
  CCTK_REAL, dimension(npoints), intent(INOUT) :: temp
  CCTK_INT, intent(OUT) :: anyerr, inv_err(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, ye, eps
  real*8, dimension(npoints), intent(OUT) :: dpdrho !, eps
  real*8, dimension(npoints), intent(INOUT) :: temp
  integer, intent(OUT) :: anyerr, inv_err(npoints)
#endif 
    
  ! local vars  
  real*8, dimension(npoints) :: ptotRow, stotRow, htotRow, AbarRow, &
       dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, tempguess
  integer, dimension(npoints) :: xnEOScalls, xnEOScalls_bisec
  integer :: j
  
  tempguess = temp  
  
  call HelmEOS_invert_et(npoints, rho, temp, ye, ptotRow, eps, &
           stotRow, htotRow, AbarRow, dpdrho, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
           pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
           xxnRow, xxpRow, xxaRow, xxatmRow, inv_err, anyerr, xnEOScalls)

  ! if Newton-Raphson inversion fails or does not converge fast enough
  ! try bisection
  if (anyerr .eq. 1) then
    do j=1, npoints
      if (inv_err(j) .ne. 0) then
        call HelmEOS_invert_et_bisec(1, rho(j), tempguess(j), ye(j), ptotRow(j), eps(j), &
               stotRow(j), htotRow(j), AbarRow(j), dpdrho(j), dptRow(j), dedRow(j), &
               detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
               pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), &
               gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
               xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), inv_err(j), anyerr, xnEOScalls_bisec(j))   
        temp(j) = tempguess(j)
      end if
    end do
  end if  

  return

end subroutine helm_eos_m_kt0_dpressdrho


subroutine helm_eos_m_kt0_dpressdrho_dpressdeps(npoints, rho, temp, ye, eps, dpdrho_eps, dpdeps_rho, &
inv_err, anyerr, nEOScalls)

implicit none

#if !STANDALONE
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, ye, eps
  CCTK_REAL, dimension(npoints), intent(OUT) :: dpdrho_eps, dpdeps_rho
  CCTK_REAL, dimension(npoints), intent(INOUT) :: temp
  CCTK_INT, intent(OUT) :: anyerr, inv_err(npoints)
  CCTK_INT, dimension(npoints), intent(INOUT) :: nEOScalls
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, ye, eps
  real*8, dimension(npoints), intent(OUT) :: dpdrho_eps, dpdeps_rho
  real*8, dimension(npoints), intent(INOUT) :: temp
  integer, intent(OUT) :: anyerr, inv_err(npoints)
  integer, dimension(npoints), intent(INOUT) :: nEOScalls
#endif

! local vars
real*8, dimension(npoints) :: ptotRow, stotRow, htotRow, AbarRow, &
  dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
  pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, &
  xxnRow, xxpRow, xxaRow, xxatmRow, tempguess
integer, dimension(npoints) :: xnEOScalls, xnEOScalls_bisec
integer :: j

tempguess = temp

call HelmEOS_invert_et(npoints, rho, temp, ye, ptotRow, eps, &
  stotRow, htotRow, AbarRow, dpdrho_eps, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
  pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpdeps_rho, &
  xxnRow, xxpRow, xxaRow, xxatmRow, inv_err, anyerr, xnEOScalls)

nEOScalls = xnEOScalls

! if Newton-Raphson inversion fails or does not converge fast enough
! try bisection
if (anyerr .eq. 1) then
  do j=1, npoints
    if (inv_err(j) .ne. 0) then
      call HelmEOS_invert_et_bisec(1, rho(j), tempguess(j), ye(j), ptotRow(j), eps(j), &
        stotRow(j), htotRow(j), AbarRow(j), dpdrho_eps(j), dptRow(j), dedRow(j), &
        detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
        pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), &
        gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpdeps_rho(j), &
        xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), inv_err(j), anyerr, xnEOScalls_bisec(j))
      temp(j) = tempguess(j)
      nEOScalls(j) = nEOScalls(j) + xnEOScalls_bisec(j)
    end if
  end do
end if

return

end subroutine helm_eos_m_kt0_dpressdrho_dpressdeps


subroutine helm_eos_m_kt0_dpressdrho_dpressdeps_hinv(npoints, rho, temp, enth, ye, eps, &
  press, dpdrho_eps, dpdeps_rho, entr, abar, inv_err, anyerr, nEOScalls)

implicit none
#if !STANDALONE
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, ye, enth
  CCTK_REAL, dimension(npoints), intent(OUT) :: dpdrho_eps, dpdeps_rho, press, eps, entr, abar
  CCTK_REAL, dimension(npoints), intent(INOUT) :: temp
  CCTK_INT, intent(OUT) :: anyerr, inv_err(npoints)
  CCTK_INT, dimension(npoints), intent(INOUT) :: nEOScalls
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, ye, enth
  real*8, dimension(npoints), intent(OUT) :: dpdrho_eps, dpdeps_rho, press, eps, entr, abar
  real*8, dimension(npoints), intent(INOUT) :: temp
  integer, intent(OUT) :: anyerr, inv_err(npoints)
  integer, dimension(npoints), intent(INOUT) :: nEOScalls
#endif

! local vars
real*8, dimension(npoints) :: dpdRow, &
  dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
  pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, &
  xxnRow, xxpRow, xxaRow, xxatmRow, tempguess
integer, dimension(npoints) :: xnEOScalls, xnEOScalls_bisec
integer :: j

tempguess = temp

call HelmEOS_invert_ht(npoints, rho, temp, ye, press, eps, &
  entr, enth, abar, dpdrho_eps, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
  pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpdeps_rho, &
  xxnRow, xxpRow, xxaRow, xxatmRow, inv_err, anyerr, xnEOScalls)

nEOScalls = xnEOScalls
!nEOScalls = 0
!xnEOScalls = 0
!anyerr = 1
!inv_err = 1
!if Newton-Raphson inversion fails or does not converge fast enough
!try bisection
!if (anyerr .eq. 1) then
!  do j=1, npoints
!    if (inv_err(j) .ne. 0) then
!!     write(*,*) 'BISECTON:',rho(j), tempguess(j), ye(j), enth(j)
!      call HelmEOS_invert_ht_bisec(1, rho(j), tempguess(j), ye(j), press(j), eps(j), &
!        entr(j), enth(j), abar(j), dpdrho_eps(j), dptRow(j), dedRow(j), &
!        detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
!        pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), &
!        gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpdeps_rho(j), &
!        xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), inv_err(j), anyerr, xnEOScalls_bisec(j))
!      temp(j) = tempguess(j)
!      nEOScalls(j) = xnEOScalls(j) + xnEOScalls_bisec(j)
!    end if
!  end do
!end if

return

end subroutine helm_eos_m_kt0_dpressdrho_dpressdeps_hinv


subroutine helm_eos_m_kt1_dpressdrho_dpressdeps_dpressdt_depsdt(npoints, rho, temp, ye, eps, press, &
  dpdrho_eps, dpdeps_rho, dpdt, depsdt, NSEkeyerr, anyerr)

implicit none

#if !STANDALONE
  CCTK_INT, intent(IN) :: npoints
  CCTK_REAL, dimension(npoints), intent(IN)  :: rho, ye, temp
  CCTK_REAL, dimension(npoints), intent(OUT) :: dpdrho_eps, dpdeps_rho, dpdt, depsdt, press
  CCTK_REAL, dimension(npoints), intent(INOUT) :: eps
  CCTK_INT, intent(OUT) :: anyerr, NSEkeyerr(npoints)
#else
  integer, intent(IN) :: npoints
  real*8, dimension(npoints), intent(IN)  :: rho, ye, temp
  real*8, dimension(npoints), intent(OUT) :: dpdrho_eps, dpdeps_rho, dpdt, depsdt, press
  real*8, dimension(npoints), intent(INOUT) :: eps
  integer, intent(OUT) :: anyerr, NSEkeyerr(npoints)
#endif

  ! local vars  
real*8, dimension(npoints) :: stotRow, htotRow, AbarRow, dedRow, dsdRow, dstRow, dhdRow, dhtRow, &
  pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, &
  xxnRow, xxpRow, xxaRow, xxatmRow



call helm_eos(npoints, rho, temp, ye, &
      eps, press, stotRow, htotRow, AbarRow, &
      dpdrho_eps, dpdt, dedRow, depsdt, dsdRow, dstRow, dhdRow, dhtRow, &
      pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpdeps_rho, &
      xxnRow, xxpRow, xxaRow, xxatmRow, NSEkeyerr, anyerr)

return

end subroutine helm_eos_m_kt1_dpressdrho_dpressdeps_dpressdt_depsdt


subroutine helm_eos_kt1_etas_xs_point(rho, temp, ye, &
       etae, etap, etan, xp, xn, xa)

  implicit none

#if !STANDALONE
  CCTK_REAL, intent(IN)  :: rho, temp, ye
  CCTK_REAL, intent(OUT) :: etae, etap, etan, xp, xn, xa
#else
  real*8, intent(IN)  :: rho, temp, ye
  real*8, intent(OUT) :: etae, etap, etan, xp, xn, xa
#endif  
  
  ! local vars
  integer :: anyerr, NSEkeyerr(1)
  real*8, dimension(1) :: rhoRow, tempRow, yeRow, epsRow, pressRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow

  rhoRow(1) = rho
  tempRow(1) = temp
  yeRow(1) = ye

  call helm_eos(1, rhoRow, tempRow, yeRow, &
       epsRow, pressRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, NSEkeyerr, anyerr)   
       
  etae = etaeRow(1)
  etap = etapRow(1)
  etan = etanRow(1)
  xp = xxpRow(1)
  xn = xxnRow(1)
  xa = xxaRow(1)

  return

end subroutine helm_eos_kt1_etas_xs_point


subroutine helm_eos_m_kt1_dpdrho_dpdt_dedrho_dedt_point(rho, temp, ye, eps, press, dpdrho, &
               dpdt, dedrho, dedt, anyerr)

  implicit none

#if !STANDALONE
  CCTK_REAL, intent(IN)  :: rho, temp, ye
  CCTK_REAL, intent(OUT) :: eps, press, dpdrho, dpdt, dedrho, dedt
  CCTK_INT, intent(OUT) :: anyerr
#else
  real*8, intent(IN)  :: rho, temp, ye
  real*8, intent(OUT) :: eps, press, dpdrho, dpdt, dedrho, dedt
  integer, intent(OUT) :: anyerr
#endif    
  
  ! local vars
  integer :: NSEkeyerr(1)
  real*8, dimension(1) :: rhoRow, tempRow, yeRow, epsRow, pressRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow

  rhoRow(1) = rho
  tempRow(1) = temp
  yeRow(1) = ye

  call helm_eos(1, rhoRow, tempRow, yeRow, &
       epsRow, pressRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, NSEkeyerr, anyerr)
       
  eps = epsRow(1)
  press = pressRow(1)
  dpdrho = dpdRow(1)
  dpdt = dptRow(1)
  dedrho = dedRow(1)
  dedt = detRow(1)

  return

end subroutine helm_eos_m_kt1_dpdrho_dpdt_dedrho_dedt_point


!--------------------------------------------------------------------------
!
! Inversion routines from specific internal energy and specific enthalpy
!
!--------------------------------------------------------------------------


subroutine HelmEOS_invert_et(num_inv, denRow, tempRow, YeRow, ptotRow, etotRow, &
               stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
               pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
               xxnRow, xxpRow, xxaRow, xxatmRow, conv_err, anyerr, numEOScalls)

  ! given the specific internal energy, density, Ye, 
  ! and an initial guess for the temperature, find everything else via
  ! Newton-Raphson root finding
      
  use EOS_Omni_Module    
      
  implicit none

  ! arguments
  integer, intent(IN)    :: num_inv
  real*8, dimension (num_inv), intent(IN) :: denRow, etotRow, YeRow
  real*8, dimension (num_inv), intent(INOUT) :: tempRow, xxnRow, xxpRow, xxaRow, xxatmRow
  real*8, dimension (num_inv), intent(OUT)   :: ptotRow, stotRow, htotRow, AbarRow, dpdRow, dptRow, &
        dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, &
        cvRow, cpRow, cs2Row, dpde_d_Row
  integer, intent(OUT)  :: conv_err(num_inv), anyerr
  integer, intent(INOUT) :: numEOScalls(num_inv)

  ! local variables
  integer :: i, j, NSEkeyerr(num_inv), anyerrloc
  real*8  :: f, df, temp, tempnew, tempmin, tempmax, yemin, yemax
  real*8, dimension(num_inv) :: error_T, error_eps, corr, eloc, qloc, initguess, YeRow_corr

  real*8, parameter  :: eostol = inv_tol
  real*8, parameter  :: eostol_eps = inv_tol
  integer, parameter :: itmax = 40       

  ! initialize
  error_T = 0.0d0
  error_eps = 0.0d0
  corr = 0.0d0
  qloc = etotRow
  initguess = tempRow
  
  conv_err = 0
  anyerr = 0
  anyerrloc = 0
  numEOScalls = 0
  
  tempmin = eos_tempmin
  tempmax = 2.0d2  ! set 200MeV arbitrarily as upper limit
  yemin = eos_yemin
  yemax = eos_yemax
    
  ! Newton-Raphson: loop over each element of the pipe individually
  do j = 1, num_inv
  
    ! limit to reasonable range
    tempRow(j)    = min( tempmax, max( tempRow(j), tempmin) )
    YeRow_corr(j) = min( yemax, max(YeRow(j), yemin) )
    
    do i=1,itmax

      call helm_eos(1, denRow(j), tempRow(j), YeRow_corr(j), &
         eloc(j), ptotRow(j), stotRow(j), htotRow(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
    
      ! increase EOS calls counter
      numEOScalls(j) = numEOScalls(j) + 1
      
      if (anyerrloc .eq. 1) then
        anyerr = 1
#if DEBUG        
        write(*,*) 'WARNING, HelmEOS_invert_et: NSE calculation did not converge somewhere...'
#endif
      end if 

      f       = eloc(j)/qloc(j) - 1.0d0 
      df      = detRow(j)/qloc(j) 
      corr(j) = f/df
    
      ! limit excursions to factor of two changes
      temp    = tempRow(j)
      tempnew = min( max( 0.5d0*temp, temp - corr(j) ), 2.0d0*temp )
      ! keep within limits
      tempnew = min( tempmax, max( tempnew, tempmin ) )

      ! compute the error
      error_T(j)  = abs( (tempnew - temp) / temp )
      error_eps(j) = abs( (eloc(j) - qloc(j)) / qloc(j) )

      ! store the new temp, keep it within limits
      tempRow(j)  = tempnew
      
      if ( error_T(j) .lt. eostol .and. error_eps(j) .lt. eostol_eps) then 
        exit
      else if (i .eq. int(0.5d0*itmax)) then
        if (tempRow(j) .eq. tempmin) then
          exit
        else
          if ( j.gt.1 ) then
            !write(*,*) 'WARNING: resetting initial guess for temperature...'
            tempRow(j) = tempRow(j-1)
          else
            cycle
          endif  
        end if  
      else if (i .eq. itmax) then 
        conv_err(j) = 1 
        anyerr = 1
        !write(*,*)
        !write(*,*) 'iteration nr.', i
        !write(*,*) 'newton-raphson failed in routine HelmEOS_invert_et'
        !write(*,*) 'pipeline element',j
        !write(*,*) 'eps  =',qloc(j)
        !write(*,*) 'rho =',denRow(j)
        !write(*,*) 'Ye =', YeRow(j)
        !write(*,*) 'eps_conv, eps =', etotRow(j), qloc(j)
        !write(*,*) 'error_T =',error_T(j), 'error_eps =',error_eps(j)
        !write(*,*) 'eostol =',eostol,'  eostol_eps =',eostol_eps
        !write(*,*) 'temp_loc  =',tempRow(j),'  temp_init=',initguess(j)
        !write(*,*) 'f/df  =',corr(j),' f   =',f,    ' df    =',df
        !write(*,*)
        !!stop 'could not find a temperature in routine HelmEOS_invert_et'
      end if
    end do
  end do

  return
  
end subroutine HelmEOS_invert_et


subroutine HelmEOS_invert_ht(num_inv, denRow, tempRow, YeRow, ptotRow, etotRow, &
               stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
               pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
               xxnRow, xxpRow, xxaRow, xxatmRow, conv_err, anyerr, numEOScalls)

  ! given the specific enthalpy h = rho + p/rho, density, Ye, 
  ! and an initial guess for the temperature, find everything else via
  ! Newton-Raphson root finding
  ! Note: non-relativistic specific enthalpy, without rest-mass contribution
      
  use EOS_Omni_Module    
      
  implicit none

  ! arguments
  integer, intent(IN)    :: num_inv
  real*8, dimension (num_inv), intent(IN) :: denRow, htotRow, YeRow
  real*8, dimension (num_inv), intent(INOUT) :: tempRow, xxnRow, xxpRow, xxaRow, xxatmRow
  real*8, dimension (num_inv), intent(OUT)   :: ptotRow, stotRow, etotRow, AbarRow, dpdRow, dptRow, &
        dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, &
        cvRow, cpRow, cs2Row, dpde_d_Row
  integer, intent(OUT)  :: conv_err(num_inv), anyerr
  integer, intent(INOUT) :: numEOScalls(num_inv)

  ! local variables
  integer :: i, j, NSEkeyerr(num_inv), anyerrloc
  real*8  :: f, df, temp, tempnew, tempmin, tempmax, yemin, yemax
  real*8, dimension(num_inv) :: error_T, error_h, corr, hloc, qloc, initguess, YeRow_corr

  real*8, parameter  :: eostol = inv_tol
  real*8, parameter  :: eostol_h = inv_tol
  integer, parameter :: itmax = 40       

  ! initialize
  error_T = 0.0d0
  error_h = 0.0d0
  corr = 0.0d0
  qloc = htotRow
  initguess = tempRow
  
  conv_err = 0
  anyerr = 0
  anyerrloc = 0
  numEOScalls = 0
  
  tempmin = eos_tempmin
  tempmax = 2.0d2  ! set 200MeV arbitrarily as upper limit
  yemin = eos_yemin
  yemax = eos_yemax
    
  ! Newton-Raphson: loop over each element of the pipe individually
  do j = 1, num_inv
  
    ! limit to reasonable range
    tempRow(j)    = min( tempmax, max( tempRow(j), tempmin) )
    YeRow_corr(j) = min( yemax, max(YeRow(j), yemin) )
  
    do i=1,itmax

      call helm_eos(1, denRow(j), tempRow(j), YeRow_corr(j), &
         etotRow(j), ptotRow(j), stotRow(j), hloc(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
    
      ! increase EOS calls counter
      numEOScalls(j) = numEOScalls(j) + 1
      
      if (anyerrloc .eq. 1) then
        anyerr = 1
#if DEBUG        
        write(*,*) 'WARNING, HelmEOS_invert_ht: NSE calculation did not converge somewhere...'
#endif      
      end if 

      f       = hloc(j)/qloc(j) - 1.0d0 
      df      = dhtRow(j)/qloc(j) 
      corr(j) = f/df
    
      ! limit excursions to factor of two changes
      temp    = tempRow(j)
      tempnew = min( max( 0.5d0*temp, temp - corr(j) ), 2.0d0*temp )
      ! keep within limits
      tempnew = min( tempmax, max( tempnew, tempmin ) )

      ! compute the error
      error_T(j)  = abs( (tempnew - temp) / temp )
      error_h(j) = abs( (hloc(j) - qloc(j)) / qloc(j) )

      ! store the new temp, keep it within limits
      tempRow(j)  = tempnew
      
      if ( error_T(j) .lt. eostol .and. error_h(j) .lt. eostol_h) then 
        exit
      else if (i .eq. int(0.5d0*itmax)) then
        if (tempRow(j) .eq. tempmin) then
          exit
        else
          if ( j.gt.1 ) then
            !write(*,*) 'WARNING: resetting initial guess for temperature...'
            tempRow(j) = tempRow(j-1)
          else
            cycle
          endif  
        end if  
      else if (i .eq. itmax) then 
        conv_err(j) = 1 
        anyerr = 1
        !write(*,*) 'h NR inversion failed!'
        !write(*,*)
        !write(*,*) 'iteration nr.', i
        !write(*,*) 'newton-raphson failed in routine HelmEOS_invert_ht'
        !write(*,*) 'pipeline element',j
        !write(*,*) 'h  =',qloc(j)
        !write(*,*) 'rho =',denRow(j)
        !write(*,*) 'Ye =', YeRow(j)
        !write(*,*) 'h_conv, h =', htotRow(j), qloc(j)
        !write(*,*) 'error_T =',error_T(j), 'error_h =',error_h(j)
        !write(*,*) 'eostol =',eostol,'  eostol_h =',eostol_h
        !write(*,*) 'temp_loc  =',tempRow(j),'  temp_init=',initguess(j)
        !write(*,*) 'f/df  =',corr(j),' f   =',f,    ' df    =',df
        !write(*,*)
        !!stop 'could not find a temperature in routine HelmEOS_invert_ht'
      end if
    end do
  end do

  return
  
end subroutine HelmEOS_invert_ht


subroutine HelmEOS_invert_et_bisec(num_inv, denRow, tempRow, YeRow, ptotRow, etotRow, &
               stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
               pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
               xxnRow, xxpRow, xxaRow, xxatmRow, conv_err, anyerr, numEOScalls)

  ! given the specific internal energy, density, Ye, 
  ! and an initial guess for the temperature, find everything else via bisection
      
  use EOS_Omni_Module    
      
  implicit none

  ! arguments
  integer, intent(IN)    :: num_inv
  real*8, dimension (num_inv), intent(IN) :: denRow, etotRow, YeRow
  real*8, dimension (num_inv), intent(INOUT) :: tempRow, xxnRow, xxpRow, xxaRow, xxatmRow
  real*8, dimension (num_inv), intent(OUT)   :: ptotRow, stotRow, htotRow, AbarRow, dpdRow, dptRow, &
        dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, &
        cvRow, cpRow, cs2Row, dpde_d_Row
  integer, intent(OUT)  :: conv_err(num_inv), anyerr
  integer, intent(INOUT) :: numEOScalls(num_inv)

  ! local variables
  integer :: i, j, NSEkeyerr(num_inv), anyerrloc, anyerr_bisec, count_bs
  real*8  :: temp_a, temp_b, temp_c, temp_c_new, tempmin, tempmax, yemin, yemax
  real*8  :: sign_f_a, sign_f_b, sign_f_c
  real*8  :: error_eps_a, error_eps_b, error_eps_c
  real*8, dimension(num_inv) :: error_T, error_eps, epsloc, initguess, temps_a, temps_b, temps_c, &
                                epss_a, epss_b, epss_c, YeRow_corr

  real*8, parameter  :: eostol_T = inv_tol
  real*8, parameter  :: eostol_eps = inv_tol
  integer, parameter :: itmax = 80
  integer, parameter :: count_bs_max = 50
    
  ! initialize
  error_T = 1.0d0
  error_eps = 1.0d0
  epsloc = etotRow
  initguess = tempRow
  
  conv_err = 0
  anyerr = 0
  anyerrloc = 0
  numEOScalls = 0
  
  ! make sure atmosphere level is included, so that it can be bracketed
  tempmin = eos_tempmin
  tempmax = 2.0d2  ! set 200MeV arbitrarily as upper limit
  yemin = eos_yemin
  yemax = eos_yemax
  
  ! bisection: loop over each element of the pipe individually
  do j = 1, num_inv
  
    anyerr_bisec = 0
    
    YeRow_corr(j)   = min( yemax, max(YeRow(j), yemin) )
  
    temp_a = max(0.98d0*initguess(j), tempmin)
    temp_b = min(1.02d0*initguess(j), tempmax)
    temps_a(j) = temp_a
    temps_b(j) = temp_b
  
    call helm_eos(1, denRow(j), temps_a(j), YeRow_corr(j), &
         epss_a(j), ptotRow(j), stotRow(j), htotRow(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
         
    call helm_eos(1, denRow(j), temps_b(j), YeRow_corr(j), &
         epss_b(j), ptotRow(j), stotRow(j), htotRow(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
    
    ! increase number of EOS calls counter
    numEOScalls(j) = numEOScalls(j) + 2
      
    sign_f_a = sign(1.0d0, (epss_a(j)/epsloc(j) - 1.0d0))
    sign_f_b = sign(1.0d0, (epss_b(j)/epsloc(j) - 1.0d0))
    
    ! first find bracket to initial guess
    count_bs = 0
    do while ((sign_f_a*sign_f_b .gt. 0.0d0) .and. (count_bs .lt. count_bs_max))
      temp_a = max(0.98d0*temp_a, tempmin)
      temp_b = min(1.02d0*temp_b, tempmax)
      temps_a(j) = temp_a
      temps_b(j) = temp_b
      
      call helm_eos(1, denRow(j), temps_a(j), YeRow_corr(j), &
         epss_a(j), ptotRow(j), stotRow(j), htotRow(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
         
      call helm_eos(1, denRow(j), temps_b(j), YeRow_corr(j), &
         epss_b(j), ptotRow(j), stotRow(j), htotRow(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
      
      ! increase number of EOS calls counter
      numEOScalls(j) = numEOScalls(j) + 2
      
      sign_f_a = sign(1.0d0, (epss_a(j)/epsloc(j) - 1.0d0))
      sign_f_b = sign(1.0d0, (epss_b(j)/epsloc(j) - 1.0d0))
      
      if ((temp_a .eq. tempmin .and. temp_b .eq. tempmax) .or. (count_bs .eq. (count_bs_max-1))) then
        !!call CCTK_ERROR("Initial guess for temperature cannot be bracketed, stopping the code...")
        !write(*,*) "It seems the initial guess for the temperature cannot be bracketed..."
        !write(*,*) 'eps =', epsloc(j)
        !write(*,*) 'den =', denRow(j)
        !write(*,*) 'ye =', YeRow_corr(j)
        !write(*,*) "pipeline element", j
        anyerr_bisec = 1
        exit
      end if
      count_bs = count_bs + 1  
    end do
  
    ! now start bisection or skip if there is no bracket
    if (anyerr_bisec .eq. 1) then
      conv_err(j) = 1
      anyerr = 1
    else
   
      temp_c = 0.5d0 * (temp_a + temp_b)
  
      temps_a(j) = temp_a
      temps_b(j) = temp_b
      temps_c(j) = temp_c
  
      call helm_eos(1, denRow(j), temps_a(j), YeRow_corr(j), &
         epss_a(j), ptotRow(j), stotRow(j), htotRow(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
         
      call helm_eos(1, denRow(j), temps_b(j), YeRow_corr(j), &
         epss_b(j), ptotRow(j), stotRow(j), htotRow(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
         
      call helm_eos(1, denRow(j), temps_c(j), YeRow_corr(j), &
         epss_c(j), ptotRow(j), stotRow(j), htotRow(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)                
    
      ! increase number of EOS calls counter
      numEOScalls(j) = numEOScalls(j) + 3
      
      sign_f_a = sign(1.0d0, (epss_a(j)/epsloc(j) - 1.0d0))
      sign_f_b = sign(1.0d0, (epss_b(j)/epsloc(j) - 1.0d0))
      sign_f_c = sign(1.0d0, (epss_c(j)/epsloc(j) - 1.0d0))
    
      error_T(j) = abs( (temp_c - temp_b)/ temp_b )
      error_eps_a = abs( (epss_a(j) - epsloc(j))/epsloc(j) )
      error_eps_b = abs( (epss_b(j) - epsloc(j))/epsloc(j) )
      error_eps_c = abs( (epss_c(j) - epsloc(j))/epsloc(j) )
    
      ! start bisection
      do i = 1,itmax
        if ( (epss_c(j)/epsloc(j) - 1.0d0) .eq. 0.0d0) then
          tempRow(j) = temp_c
          error_eps(j) = error_eps_c
          exit
        else if (sign_f_a * sign_f_c .lt. 0.0d0) then
          temp_b = temp_c
          temps_b(j) = temp_c
          epss_b(j) = epss_c(j)
          error_eps_b = error_eps_c
          sign_f_b = sign_f_c
        else
          temp_a = temp_c
          temps_a(j) = temp_c
          epss_a(j) = epss_c(j)
          error_eps_a = error_eps_c
          sign_f_a = sign_f_c
        end if 
       
        temp_c_new =  0.5d0 * (temp_a + temp_b)
      
        temps_c(j) = temp_c_new
      
        call helm_eos(1, denRow(j), temps_c(j), YeRow_corr(j), &
         epss_c(j), ptotRow(j), stotRow(j), htotRow(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), &
         cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)    
      
        ! increase number of EOS calls counter
        numEOScalls(j) = numEOScalls(j) + 1
      
        if (anyerrloc .eq. 1) then
          anyerr = 1
#if DEBUG          
          write(*,*) 'WARNING, HelmEOS_invert_et_bisec: NSE calculation did not converge...'
#endif        
        end if
      
        sign_f_c = sign(1.0d0, (epss_c(j)/epsloc(j) - 1.0d0))
        error_T(j) =  abs( (temp_c_new - temp_c)/ temp_c )
        error_eps_c = abs( (epss_c(j) - epsloc(j))/epsloc(j) )  
        temp_c = temp_c_new
      
        if (error_T(j) .lt. eostol_T .and. error_eps_c .lt. eostol_eps) then
          tempRow(j) = temp_c
          error_eps(j) = error_eps_c
          exit
        else if (i .eq. itmax) then
          tempRow(j) = temp_c
          error_eps(j) = error_eps_c
          conv_err(j) = 1 
          anyerr = 1
          !write(*,*)
          !write(*,*) 'iteration nr.', i
          !write(*,*) 'bisection failed in routine invert_helm_et_bisec'
          !write(*,*) 'pipeline element',j
          !write(*,*) 'eps  =',epsloc(j)
          !write(*,*) 'rho =',denRow(j)
          !write(*,*) 'Ye =', YeRow_corr(j)
          !write(*,*) 'eps_conv, eps =', epss_c(j), epsloc(j)
          !write(*,*) 'error_T =', error_T(j), 'error_eps =', error_eps(j)
          !write(*,*) 'eostol_T =', eostol_T, 'eostol_eps =', eostol_eps
          !write(*,*) 'temp_loc  =', tempRow(j),'  temp_init=', initguess(j)
          !write(*,*)
          !!write(*,*) "bisection failed, stopping the code..."
          !!call CCTK_ERROR("bisection failed, stopping the code...")
          
        end if
      end do
    end if
  end do

  return
  
end subroutine HelmEOS_invert_et_bisec


subroutine HelmEOS_invert_ht_bisec(num_inv, denRow, tempRow, YeRow, ptotRow, etotRow, &
               stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
               pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
               xxnRow, xxpRow, xxaRow, xxatmRow, conv_err, anyerr, numEOScalls)

  ! given the specific enthalpy h = eps + p/rho, density, Ye, 
  ! and an initial guess for the temperature, find everything else via bisection
  ! Note: non-relativistic enthalpy, rest mass contribution not included!     
      
  use EOS_Omni_Module    
      
  implicit none

  ! arguments
  integer, intent(IN)    :: num_inv
  real*8, dimension (num_inv), intent(IN) :: denRow, htotRow, YeRow
  real*8, dimension (num_inv), intent(INOUT) :: tempRow
  real*8, dimension (num_inv), intent(OUT)   :: ptotRow, stotRow, etotRow, AbarRow, dpdRow, dptRow, &
        dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, &
        cvRow, cpRow, cs2Row, dpde_d_Row, xxnRow, xxpRow, xxaRow, xxatmRow
  integer, intent(OUT)  :: conv_err(num_inv), anyerr
  integer, intent(INOUT) :: numEOScalls(num_inv)

  ! local variables
  integer :: i, j, NSEkeyerr(num_inv), anyerrloc, anyerr_bisec, count_bs
  real*8  :: temp_a, temp_b, temp_c, temp_c_new, tempmin, tempmax, yemin, yemax
  real*8  :: sign_f_a, sign_f_b, sign_f_c
  real*8  :: error_h_a, error_h_b, error_h_c
  real*8, dimension(num_inv) :: error_T, error_h, hloc, initguess, temps_a, temps_b, temps_c, &
                                hs_a, hs_b, hs_c, YeRow_corr

  real*8, parameter  :: eostol_T = inv_tol
  real*8, parameter  :: eostol_h = inv_tol
  integer, parameter :: itmax = 80
  integer, parameter :: count_bs_max = 80

#if DEBUG    
  write(*,*) "in bisec..."  
#endif
  ! initialize
  error_T = 1.0d0
  error_h = 1.0d0
  hloc = htotRow
  initguess = tempRow
  
  conv_err = 0
  anyerr = 0
  anyerrloc = 0
  numEOScalls = 0
  
  ! make sure atmosphere level is included, so that it can be bracketed
  tempmin = eos_tempmin
  tempmax = 2.0d2  ! set 200MeV arbitrarily as upper limit
  yemin = eos_yemin
  yemax = eos_yemax
  
  ! bisection: loop over each element of the pipe individually
  do j = 1, num_inv
  
    anyerr_bisec = 0
    
    YeRow_corr(j)   = min( yemax, max(YeRow(j), yemin) )
  
    temp_a = max(0.98d0*initguess(j), tempmin)
    temp_b = min(1.02d0*initguess(j), tempmax)
    temps_a(j) = temp_a
    temps_b(j) = temp_b
  
    call helm_eos(1, denRow(j), temps_a(j), YeRow_corr(j), &
         etotRow(j), ptotRow(j), stotRow(j), hs_a(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
         
    call helm_eos(1, denRow(j), temps_b(j), YeRow_corr(j), &
         etotRow(j), ptotRow(j), stotRow(j), hs_b(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
    
    ! increase number of EOS calls counter
    numEOScalls(j) = numEOScalls(j) + 2
      
    sign_f_a = sign(1.0d0, (hs_a(j)/hloc(j) - 1.0d0))
    sign_f_b = sign(1.0d0, (hs_b(j)/hloc(j) - 1.0d0))
    
    ! first find bracket to initial guess
    count_bs = 0
    do while ((sign_f_a*sign_f_b .gt. 0.0d0) .and. (count_bs .lt. count_bs_max))
    
      temp_a = max(0.98d0*temp_a, tempmin)
      temp_b = min(1.02d0*temp_b, tempmax)
      temps_a(j) = temp_a
      temps_b(j) = temp_b
#if DEBUG      
      write(*,*) count_bs, temp_a, temp_b, hs_a(j), hs_b(j), hloc(j)
#endif      
      call helm_eos(1, denRow(j), temps_a(j), YeRow_corr(j), &
         etotRow(j), ptotRow(j), stotRow(j), hs_a(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
         
      call helm_eos(1, denRow(j), temps_b(j), YeRow_corr(j), &
         etotRow(j), ptotRow(j), stotRow(j), hs_b(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
      
      ! increase number of EOS calls counter
      numEOScalls(j) = numEOScalls(j) + 2
      
      sign_f_a = sign(1.0d0, (hs_a(j)/hloc(j) - 1.0d0))
      sign_f_b = sign(1.0d0, (hs_b(j)/hloc(j) - 1.0d0))
      
      if ((temp_a .eq. tempmin .and. temp_b .eq. tempmax) .or. (count_bs .eq. (count_bs_max-1))) then
        !!call CCTK_ERROR("Initial guess for temperature cannot be bracketed, stopping the code...")
        !write(*,*) 'h =', hloc(j)
        !write(*,*) 'den =', denRow(j)
        !write(*,*) 'ye =', YeRow_corr(j)
        !write(*,*) "pipeline element", j
        anyerr_bisec = 1
        exit
      end if
      count_bs = count_bs + 1  
    end do
  
    ! now start bisection or skip if there is no bracket
    if (anyerr_bisec .eq. 1) then
      conv_err(j) = 1
      anyerr = 1
#if DEBUG
  write(*,*) 'bisection h: temperature cannot be bracketed...'
#endif      
    else
   
      temp_c = 0.5d0 * (temp_a + temp_b)
  
      temps_a(j) = temp_a
      temps_b(j) = temp_b
      temps_c(j) = temp_c
  
      call helm_eos(1, denRow(j), temps_a(j), YeRow_corr(j), &
         etotRow(j), ptotRow(j), stotRow(j), hs_a(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
         
      call helm_eos(1, denRow(j), temps_b(j), YeRow_corr(j), &
         etotRow(j), ptotRow(j), stotRow(j), hs_b(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)
         
      call helm_eos(1, denRow(j), temps_c(j), YeRow_corr(j), &
         etotRow(j), ptotRow(j), stotRow(j), hs_c(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)                
    
      ! increase number of EOS calls counter
      numEOScalls(j) = numEOScalls(j) + 3
      
      sign_f_a = sign(1.0d0, (hs_a(j)/hloc(j) - 1.0d0))
      sign_f_b = sign(1.0d0, (hs_b(j)/hloc(j) - 1.0d0))
      sign_f_c = sign(1.0d0, (hs_c(j)/hloc(j) - 1.0d0))
    
      error_T(j) = abs( (temp_c - temp_b)/ temp_b )
      error_h_a = abs( (hs_a(j) - hloc(j))/hloc(j) )
      error_h_b = abs( (hs_b(j) - hloc(j))/hloc(j) )
      error_h_c = abs( (hs_c(j) - hloc(j))/hloc(j) )
    
      ! start bisection
      do i = 1,itmax
        if ( (hs_c(j)/hloc(j) - 1.0d0) .eq. 0.0d0) then
          tempRow(j) = temp_c
          error_h(j) = error_h_c
          exit
        else if (sign_f_a * sign_f_c .lt. 0.0d0) then
          temp_b = temp_c
          temps_b(j) = temp_c
          hs_b(j) = hs_c(j)
          error_h_b = error_h_c
          sign_f_b = sign_f_c
        else
          temp_a = temp_c
          temps_a(j) = temp_c
          hs_a(j) = hs_c(j)
          error_h_a = error_h_c
          sign_f_a = sign_f_c
        end if 
       
        temp_c_new =  0.5d0 * (temp_a + temp_b)
      
        temps_c(j) = temp_c_new
      
        call helm_eos(1, denRow(j), temps_c(j), YeRow_corr(j), &
         etotRow(j), ptotRow(j), stotRow(j), hs_c(j), AbarRow(j), &
         dpdRow(j), dptRow(j), dedRow(j), detRow(j), dsdRow(j), dstRow(j), dhdRow(j), dhtRow(j), &
         pelRow(j), neRow(j), etaeRow(j), etapRow(j), etanRow(j), gam1Row(j), cvRow(j), &
         cpRow(j), cs2Row(j), dpde_d_Row(j), &
         xxnRow(j), xxpRow(j), xxaRow(j), xxatmRow(j), NSEkeyerr(j), anyerrloc)    
      
        ! increase number of EOS calls counter
        numEOScalls(j) = numEOScalls(j) + 1
      
        if (anyerrloc .eq. 1) then
          anyerr = 1
#if DEBUG          
          write(*,*) 'WARNING, HelmEOS_invert_et_bisec: NSE calculation did not converge...'
#endif        
        end if
      
        sign_f_c = sign(1.0d0, (hs_c(j)/hloc(j) - 1.0d0))
        error_T(j) =  abs( (temp_c_new - temp_c)/ temp_c )
        error_h_c = abs( (hs_c(j) - hloc(j))/hloc(j) )  
        temp_c = temp_c_new
        !write(*,*) temp_c, hs_c(j), ptotRow(j)
        if (error_T(j) .lt. eostol_T .and. error_h_c .lt. eostol_h) then
          tempRow(j) = temp_c
          error_h(j) = error_h_c
          exit
        else if (i .eq. itmax) then
          tempRow(j) = temp_c
          error_h(j) = error_h_c
          conv_err(j) = 1 
          anyerr = 1
#if DEBUG
          !write(*,*)
          !write(*,*) 'iteration nr.', i
          write(*,*) 'bisection failed in routine invert_helm_et_bisec'
          !write(*,*) 'pipeline element',j
          write(*,*) 'h  =',hloc(j)
          write(*,*) 'rho =',denRow(j)
          !write(*,*) 'Ye =', YeRow_corr(j)
          !write(*,*) 'h_conv, h =', hs_c(j), hloc(j)
          write(*,*) 'error_T =', error_T(j), 'error_h =', error_h(j)
          write(*,*) 'eostol_T =', eostol_T, 'eostol_h =', eostol_h
          write(*,*) 'temp_loc  =', tempRow(j),'  temp_init=', initguess(j)
          !write(*,*)
          !!write(*,*) "bisection failed, stopping the code..."
          !call CCTK_ERROR("bisection failed, stopping the code...")
#endif          
        end if
      end do
    end if
  end do

  return
  
end subroutine HelmEOS_invert_ht_bisec


!----------------------------------------------------------------
!
! Helmholtz EOS routine
!
!----------------------------------------------------------------


subroutine helm_eos(npoints, denRow, tempRow, YeRow, &
       etotRow, ptotRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row, &
       xxnRow, xxpRow, xxaRow, xxatmRow, NSEkeyerr, anyerr)

! This is a substantially modified implementation of the Helmholtz EOS, 
! based on the implementation by Timmes & Arnett 1999, ApJS, 125, 277 and Timmes & Swesty 2000,
! ApJS 126, 501, see http://cococubed.asu.edu/code_pages/eos.shtml,
! and based on modifications by R. Fernandez used in Fernandez & Metzger 2013, MNRAS 435, 502.

! This implementation includes abundance calculations for neutrons, 
! protons, and alpha particles in nuclear statistical equilibrium (NSE), as well as an
! implementation of analytic formulae for an extremely ultrarelativistic Fermi gas 
! when the limits of the EOS table are exceeded (T > 1e11K, rho > 1e11 g/cm^3).

! Note: all calculations are performed in cgs units, so conversion between Cactus code units
! and cgs units is required (see below)

! Note: this EOS assumes a baryon mass equal to the atomic mass unit,
! mb = mu = 1/N_A = 1.6605402e-24 g, where N_A is Avogadro's constant.

! Given temperature, density, and electron fraction Ye, this routine returns all
! other thermodynamic quantities.

! The fluid is assumed to be composed of planckian photons, 
! an ideal gas of ions with Coulomb corrections,
! and a fully ionized electron-positron gas with an arbitrary degree 
! of relativity and degeneracy. Interpolation in a table of the 
! Helmholtz free energy is used to return the electron-positron 
! thermodynamic quantities.
!  
!  references: cox & giuli chapter 24 ; timmes & swesty apj 1999
!
!  btemp    = temperature
!  den      = density
!  abar     = average number of nucleons per nuclei
!  zbar     = average number of protons per nuclei
!  z2bar    = square of zbar
!  ytot1    = total number of moles per gram
!  ye       = electron fraction
!

  use EOS_Omni_Module
  
  implicit none

  !DECLARE_CCTK_PARAMETERS

  ! arguments
#if !STANDALONE
   CCTK_INT, intent(IN) :: npoints
   CCTK_INT, intent(OUT) :: anyerr, NSEkeyerr(npoints)
   CCTK_REAL, dimension(npoints), intent(IN)     :: denRow, tempRow, YeRow
   CCTK_REAL, dimension(npoints), intent(OUT) :: xxnRow, xxpRow, xxaRow, xxatmRow
   CCTK_REAL, dimension(npoints), intent(OUT) :: ptotRow, etotRow, stotRow, htotRow, AbarRow, &
        dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
        pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row
#else
!f2py integer intent(in) :: npoints
!f2py integer intent(out) :: anyerr, NSEkeyerr
!f2py intent(in) :: denRow, tempRow, YeRow
!f2py intent(out) :: xxnRow, xxpRow, xxaRow, xxatmRow
!f2py intent(out) :: ptotRow, etotRow, stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow
!f2py intent(out) :: detRow, dsdRow, dstRow, dhdRow, dhtRow, pelRow, neRow, etaeRow
!f2py intent(out) :: etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row
!f2py depend(npoints) :: denRow, tempRow, YeRow
!f2py depend(npoints) :: NSEkeyerr
!f2py depend(npoints) :: xxnRow, xxpRow, xxaRow, xxatmRow
!f2py depend(npoints) :: ptotRow, etotRow, stotRow, htotRow, AbarRow, dpdRow, dptRow, dedRow
!f2py depend(npoints) :: detRow, dsdRow, dstRow, dhdRow, dhtRow, pelRow, neRow, etaeRow
!f2py depend(npoints) :: etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row
  integer, intent(IN) :: npoints
  integer, intent(OUT) :: anyerr
  integer, dimension(npoints),intent(OUT) :: NSEkeyerr
  real*8, dimension(npoints), intent(IN)  :: denRow, tempRow, YeRow
  real*8, dimension(npoints), intent(OUT) :: xxnRow, xxpRow, xxaRow, xxatmRow
  real*8, dimension(npoints), intent(OUT) :: ptotRow, etotRow, stotRow, htotRow, AbarRow, &
       dpdRow, dptRow, dedRow, detRow, dsdRow, dstRow, dhdRow, dhtRow, &
       pelRow, neRow, etaeRow, etapRow, etanRow, gam1Row, cvRow, cpRow, cs2Row, dpde_d_Row
#endif       

  
  ! local parameters
  
  ! constants
  double precision,parameter :: pi = acos(-1d0), &
       !Ggrav   = 6.67259d-8,          &
       hplanck = 6.6260755d-27,       &
       hbar    = 1.0545727d-27,       &
       avo     = 6.0221367d23,        &
       v_c     = 2.99792458d10,       &
       !kBerg   = 1.380658d-16,        &
       kBev    = 8.617385d-5,         &
       !amu     = 1.6605402d-24,       &
       msol    = 1.9892d33,           &
       ssol    = 5.67051d-5,          &
       asol    = 4.0d0 * ssol / v_c
       
  !.. derived constants
  double precision, parameter :: avoinv  = 1.0d0/avo, & 
       kergavo = kBerg * avo, &
       asoli3  = asol/3.0d0, &
       sioncon = (2.0d0 * pi * amu * kBerg)/(hplanck*hplanck)

  ! conversion factors
  !real*8, parameter :: conv_rho_CU2CGS = 1.0d0 / ( (Ggrav / v_c**2d0 )**3d0 * msol**2d0 )
  !real*8, parameter :: conv_rho_CGS2CU = ( Ggrav / v_c**2d0 )**3d0 * msol**2d0
  !real*8, parameter :: conv_t_MeV2K = 1.0d6 / kBev
  !real*8, parameter :: temp_gf = 1.0d-6 * kBev
  !real*8, parameter :: conv_press_CU2CGS = v_c**8d0 / ( Ggrav**3d0 * msol**2d0 )
  !real*8, parameter :: conv_press_CGS2CU = ( Ggrav / v_c**4d0 ) * ( Ggrav / v_c**2d0 )**2d0 * msol**2d0
  !real*8, parameter :: conv_eps_CU2CGS = v_c**2d0
  !real*8, parameter :: conv_eps_CGS2CU = 1.0d0 / v_c**2d0
  real*8, parameter :: conv_s_CU2CGS = kBerg / amu   ! Note: entropy in CU is in kB/baryon
  real*8, parameter :: conv_s_CGS2CU = amu / kBerg
  real*8, parameter :: conv_h_CU2CGS = inv_eps_gf
  real*8, parameter :: conv_h_CGS2CU = eps_gf
  
  integer ::  i, j, NSEconverge
  character(len=300) warnline

  real*8 ::  btemp,den,abar,zbar,z2bar,ytot1,ye
  real*8 ::  local_coulombMult

  real*8 ::  x1,x2,x3,x4,x5,x6,x7, &
                   y0,y1,y2,y3,y4, &  ! don't reuse variables -- it's confusing!
       deni,tempi,kt, & 
       prad,dpraddd,dpraddt,erad,deraddd,deraddt, &
       srad, dsraddd, dsraddt, & 
       xni,pion,dpiondd,dpiondt,eion,deiondd, deiondt,& 
       sion,dsiondd, dsiondt, & 
       pele,dpepdd,dpepdt,eele,deepdd,deepdt, &
       sele, dsepdd,dsepdt, & 
       pres,dpresdd,dpresdt,ener,denerdt, &
       entr, & 
       presi,chit,chid,gam1,gam2,gam3,nabad,kavoy

  real*8 :: cv, cp, zz, sound, soundCU, etaele, etap, etan, xnefer, &
         denerdd, dentrdd, dentrdt, dptotdetot_d


  !  for the interpolations
  integer :: iat,jat
  real*8 ::  free,df_d,df_t,df_dd,df_tt,df_dt
  real*8 ::  xt,xd,mxt,mxd, & 
       si0t,si1t,si2t,si0mt,si1mt,si2mt, & 
       si0d,si1d,si2d,si0md,si1md,si2md, & 
       dsi0t,dsi1t,dsi2t,dsi0mt,dsi1mt,dsi2mt, & 
       dsi0d,dsi1d,dsi2d,dsi0md,dsi1md,dsi2md, & 
       ddsi0t,ddsi1t,ddsi2t,ddsi0mt,ddsi1mt,ddsi2mt, & 
       zFunc,z0,z1,z2,z3,z4,z5,z6, & ! Split up confusing calculations
       psi0,dpsi0,ddpsi0,psi1,dpsi1,ddpsi1,psi2, & 
       dpsi2,ddpsi2,din,h5, & 
       xpsi0,xpsi1,h3dpd, & 
       w0t,w1t,w2t,w0mt,w1mt,w2mt, & 
       w0d,w1d,w2d,w0md,w1md,w2md
  real*8 :: h3e, h3x
  real*8 :: fi(36)


  !  for the coulomb corrections
  real*8 ::  ktinv,dxnidd,dsdd,lami,inv_lami,lamidd, & 
       s0,s1,s2,s3,s4, &   ! all temporary variables, reuse is confusing
       plasg,plasg_inv,plasgdd,plasgdt,a1,b1,c1,d1cc,e1cc,a2,b2,c2, & 
       ecoul,decouldd,decouldt,pcoul,dpcouldd,dpcouldt, & 
       scoul,dscouldd,dscouldt


  real*8, parameter :: third = 1.0d0/3.0d0, & 
       forth = 4.0d0/3.0d0, & 
       qe    = 4.8032068d-10,   & 
       esqu  = qe * qe


  !  for the uniform background coulomb correction
  data a1,b1,c1,d1cc,e1cc,a2,b2,c2 & 
       /-0.898004d0, 0.96786d0, 0.220703d0, -0.86097d0, & 
       2.5269d0  , 0.29561d0, 1.9885d0,    0.288675d0/


  ! variables for NSE calculations
  real*8            :: xxn, xxp, xxa, xxn_r, xxn_t, xxn_y, &
                                      xxp_r, xxp_t, xxp_y, &
                                      xxa_r, xxa_t, xxa_y
  real*8            :: xx_atm, xxa_guess
  real*8            :: dadd, dadt, dady
  !real*8, parameter :: Qalpha = 28.29566*1.60217648740d-6
  !logical, save   :: firstCall = .true.



! --------------------------------------------------------------------------

  !! ***********Beginning of statement function declarations **********

  !! The next few statements are effectively function declarations.
  !! Normally they should have been declared as, for instance 
  !! real function psi0(z)
  !!   real,intent(in) :: z
  !!   psi0=z**3 * ( z * (-6.0e0*z + 15.0e0) -10.0e0) + 1.0e0
  !! end function psi0

!!  quintic hermite polynomial statement functions

!!  psi0 and its derivatives
  psi0(zFunc)   = zFunc**3 * ( zFunc * (-6.0d0*zFunc + 15.0d0) -10.0d0) + 1.0d0
  dpsi0(zFunc)  = zFunc**2 * ( zFunc * (-30.0d0*zFunc + 60.0d0) - 30.0d0)
  ddpsi0(zFunc) = zFunc* ( zFunc*( -120.0d0*zFunc + 180.0d0) -60.0d0)


!!  psi1 and its derivatives
  psi1(zFunc)   = zFunc* ( zFunc**2 * ( zFunc * (-3.0d0*zFunc + 8.0d0) - 6.0d0) + 1.0d0)
  dpsi1(zFunc)  = zFunc*zFunc * ( zFunc * (-15.0d0*zFunc + 32.0d0) - 18.0d0) +1.0d0
  ddpsi1(zFunc) = zFunc * (zFunc * (-60.0d0*zFunc + 96.0d0) -36.0d0)


!!  psi2  and its derivatives
  psi2(zFunc)   = 0.5d0*zFunc*zFunc*( zFunc* ( zFunc * (-zFunc + 3.0d0) - 3.0d0) + 1.0d0)
  dpsi2(zFunc)  = 0.5d0*zFunc*( zFunc*(zFunc*(-5.0d0*zFunc + 12.0d0) - 9.0d0) + 2.0d0)
  ddpsi2(zFunc) = 0.5d0*(zFunc*( zFunc * (-20.0d0*zFunc + 36.0d0) - 18.0d0) + 2.0d0)


!! Define this choice if you want to decrease the operations.  It
!!  does change results slightly.

  h5(w0t,w1t,w2t,w0mt,w1mt,w2mt,w0d,w1d,w2d,w0md,w1md,w2md)= & 
       fi(1)  *w0d*w0t   + fi(2)  *w0md*w0t & 
       + fi(3)  *w0d*w0mt  + fi(4)  *w0md*w0mt & 
       + fi(5)  *w0d*w1t   + fi(6)  *w0md*w1t & 
       + fi(7)  *w0d*w1mt  + fi(8)  *w0md*w1mt & 
       + fi(9)  *w0d*w2t   + fi(10) *w0md*w2t & 
       + fi(11) *w0d*w2mt  + fi(12) *w0md*w2mt & 
       + fi(13) *w1d*w0t   + fi(14) *w1md*w0t & 
       + fi(15) *w1d*w0mt  + fi(16) *w1md*w0mt & 
       + fi(17) *w2d*w0t   + fi(18) *w2md*w0t & 
       + fi(19) *w2d*w0mt  + fi(20) *w2md*w0mt & 
       + fi(21) *w1d*w1t   + fi(22) *w1md*w1t & 
       + fi(23) *w1d*w1mt  + fi(24) *w1md*w1mt & 
       + fi(25) *w2d*w1t   + fi(26) *w2md*w1t & 
       + fi(27) *w2d*w1mt  + fi(28) *w2md*w1mt & 
       + fi(29) *w1d*w2t   + fi(30) *w1md*w2t & 
       + fi(31) *w1d*w2mt  + fi(32) *w1md*w2mt & 
       + fi(33) *w2d*w2t   + fi(34) *w2md*w2t & 
       + fi(35) *w2d*w2mt  + fi(36) *w2md*w2mt

!  !! Alternate attempt at reducing multiplications
!  h5(w0t,w1t,w2t,w0mt,w1mt,w2mt,w0d,w1d,w2d,w0md,w1md,w2md)= & 
!         (fi(1)  *w0d  + fi(2)  *w0md   & 
!       +  fi(17) *w2d  + fi(18) *w2md   & 
!       +  fi(13) *w1d  + fi(14) *w1md) *w0t & 
!       + (fi(3)  *w0d  + fi(4)  *w0md   & 
!       +  fi(15) *w1d  + fi(16) *w1md   & 
!       +  fi(19) *w2d  + fi(20) *w2md) *w0mt & 
!       + (fi(5)  *w0d  + fi(6)  *w0md   & 
!       +  fi(21) *w1d  + fi(22) *w1md   & 
!       +  fi(25) *w2d  + fi(26) *w2md) *w1t & 
!       + (fi(7)  *w0d  + fi(8)  *w0md   & 
!       +  fi(23) *w1d  + fi(24) *w1md   & 
!       +  fi(27) *w2d  + fi(28) *w2md) *w1mt & 
!       + (fi(9)  *w0d  + fi(10) *w0md   & 
!       +  fi(29) *w1d  + fi(30) *w1md   & 
!       +  fi(33) *w2d  + fi(34) *w2md) *w2t & 
!       + (fi(11) *w0d  + fi(12) *w0md   & 
!       +  fi(31) *w1d  + fi(32) *w1md   & 
!       +  fi(35) *w2d  + fi(36) *w2md) *w2mt

!!  cubic hermite polynomial statement functions
!!  psi0 & derivatives
  xpsi0(zFunc)  = zFunc * zFunc * (2.0d0*zFunc - 3.0d0) + 1.0
!UNUSED  xdpsi0(zFunc) = zFunc * (6.0e0*zFunc - 6.0e0)


!!  psi1 & derivatives
  xpsi1(zFunc)  = zFunc * ( zFunc * (zFunc - 2.0d0) + 1.0d0)
!UNUSED  xdpsi1(zFunc) = zFunc * (3.0e0*zFunc - 4.0e0) + 1.0e0



!!  bicubic hermite for the pressure derivative
  h3dpd(i,j,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) =  & 
       heos_dpdf(i,j)        *w0d*w0t   +   heos_dpdf(i+1,j)    *w0md*w0t  & 
       +   heos_dpdf(i,j+1)  *w0d*w0mt  +   heos_dpdf(i+1,j+1)  *w0md*w0mt & 
       +  heos_dpdft(i,j)    *w0d*w1t   +  heos_dpdft(i+1,j)    *w0md*w1t  & 
       +  heos_dpdft(i,j+1)  *w0d*w1mt  +  heos_dpdft(i+1,j+1)  *w0md*w1mt & 
       +  heos_dpdfd(i,j)    *w1d*w0t   +  heos_dpdfd(i+1,j)    *w1md*w0t  & 
       +  heos_dpdfd(i,j+1)  *w1d*w0mt  +  heos_dpdfd(i+1,j+1)  *w1md*w0mt & 
       + heos_dpdfdt(i,j)    *w1d*w1t   + heos_dpdfdt(i+1,j)    *w1md*w1t  & 
       + heos_dpdfdt(i,j+1)  *w1d*w1mt  + heos_dpdfdt(i+1,j+1)  *w1md*w1mt

!  h3dpd(i,j,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) =  & 
!         (heos_dpdf(i,j)     *w0d  +    heos_dpdf(i+1,j)   *w0md        & 
!       +  heos_dpdfd(i,j)    *w1d  +  heos_dpdfd(i+1,j)    *w1md) *w0t  & 
!       + (heos_dpdf(i,j+1)   *w0d  +   heos_dpdf(i+1,j+1)  *w0md        & 
!       +  heos_dpdfd(i,j+1)  *w1d  +  heos_dpdfd(i+1,j+1)  *w1md) *w0mt & 
!       + (heos_dpdft(i,j)    *w0d  +  heos_dpdft(i+1,j)    *w0md        & 
!       +  heos_dpdfdt(i,j)   *w1d  + heos_dpdfdt(i+1,j)    *w1md) *w1t  & 
!       + (heos_dpdft(i,j+1)  *w0d  +  heos_dpdft(i+1,j+1)  *w0md        & 
!       +  heos_dpdfdt(i,j+1) *w1d  + heos_dpdfdt(i+1,j+1)  *w1md) *w1mt


!!  bicubic hermite polynomial for the chemical potential
  h3e(i,j,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) =  & 
           heos_ef(i,j)    *w0d*w0t   +   heos_ef(i+1,j)    *w0md*w0t  & 
       +   heos_ef(i,j+1)  *w0d*w0mt  +   heos_ef(i+1,j+1)  *w0md*w0mt & 
       +  heos_eft(i,j)    *w0d*w1t   +  heos_eft(i+1,j)    *w0md*w1t  & 
       +  heos_eft(i,j+1)  *w0d*w1mt  +  heos_eft(i+1,j+1)  *w0md*w1mt & 
       +  heos_efd(i,j)    *w1d*w0t   +  heos_efd(i+1,j)    *w1md*w0t  & 
       +  heos_efd(i,j+1)  *w1d*w0mt  +  heos_efd(i+1,j+1)  *w1md*w0mt & 
       + heos_efdt(i,j)    *w1d*w1t   + heos_efdt(i+1,j)    *w1md*w1t  & 
       + heos_efdt(i,j+1)  *w1d*w1mt  + heos_efdt(i+1,j+1)  *w1md*w1mt

!  h3e(i,j,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) =  & 
!          (heos_ef(i,j)    *w0d   +   heos_ef(i+1,j)    *w0md        & 
!       +   heos_efd(i,j)   *w1d   +  heos_efd(i+1,j)    *w1md) *w0t  & 
!       +  (heos_ef(i,j+1)  *w0d   +   heos_ef(i+1,j+1)  *w0md        & 
!       +   heos_efd(i,j+1) *w1d   +  heos_efd(i+1,j+1)  *w1md) *w0mt & 
!       +  (heos_eft(i,j)   *w0d   +  heos_eft(i+1,j)    *w0md        & 
!      +   heos_efdt(i,j)  *w1d   + heos_efdt(i+1,j)    *w1md) *w1t  & 
!       +  (heos_eft(i,j+1) *w0d   +  heos_eft(i+1,j+1)  *w0md        & 
!       +   heos_efdt(i,j+1)*w1d   + heos_efdt(i+1,j+1)  *w1md) *w1mt



!!  bicubic hermite polynomial for electron positron number densities
  h3x(i,j,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) =  & 
           heos_xf(i,j)    *w0d*w0t   +   heos_xf(i+1,j)    *w0md*w0t  & 
       +   heos_xf(i,j+1)  *w0d*w0mt  +   heos_xf(i+1,j+1)  *w0md*w0mt & 
       +  heos_xft(i,j)    *w0d*w1t   +  heos_xft(i+1,j)    *w0md*w1t  & 
       +  heos_xft(i,j+1)  *w0d*w1mt  +  heos_xft(i+1,j+1)  *w0md*w1mt & 
       +  heos_xfd(i,j)    *w1d*w0t   +  heos_xfd(i+1,j)    *w1md*w0t  & 
       +  heos_xfd(i,j+1)  *w1d*w0mt  +  heos_xfd(i+1,j+1)  *w1md*w0mt & 
       + heos_xfdt(i,j)    *w1d*w1t   + heos_xfdt(i+1,j)    *w1md*w1t  & 
       + heos_xfdt(i,j+1)  *w1d*w1mt  + heos_xfdt(i+1,j+1)  *w1md*w1mt

!  h3x(i,j,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) =  & 
!         (heos_xf(i,j)    *w0d   +   heos_xf(i+1,j)    *w0md        & 
!       +  heos_xfd(i,j)   *w1d   +  heos_xfd(i+1,j)    *w1md) *w0t  & 
!       + (heos_xf(i,j+1)  *w0d   +   heos_xf(i+1,j+1)  *w0md        & 
!       +  heos_xfd(i,j+1) *w1d   +  heos_xfd(i+1,j+1)  *w1md) *w0mt & 
!       + (heos_xft(i,j)   *w0d   +  heos_xft(i+1,j)    *w0md        & 
!       +  heos_xfdt(i,j)  *w1d   + heos_xfdt(i+1,j)    *w1md) *w1t  & 
!       + (heos_xft(i,j+1) *w0d   +  heos_xft(i+1,j+1)  *w0md        & 
!       +  heos_xfdt(i,j+1)*w1d   + heos_xfdt(i+1,j+1)  *w1md) *w1mt

  !!************ This is the end statement function definitions *************
! -----------------------------------------------------------------------------

!!  popular format statements
!03 format(1x,4(a,1pe11.3))
!04 format(1x,4(a,i4))

! ------------------------------------------------------------------------------

  
  anyerr = 0

  do j=1,npoints

     NSEkeyerr(j) = 0

     ! convert density and temperature to cgs units
     
     btemp  = tempRow(j) * inv_temp_gf
     den    = denRow(j) * inv_rho_gf
     ye     = YeRow(j)
     !xx_atm = xxatmRow(j)
     
     ! ************************************
     !  NSE CALCULATION
     ! ************************************
     
     ! calculate abundances
     ! assume NSE is valid everywhere
     !if (tempRow(j) .gt. helmeos_temp_atmo .and. ) then
     !if (denRow(j) .gt. helmeos_rho_atmo) then
        ! NSE is valid
        
        if (ye .gt. 1.0d0 .or. ye .lt. 0.0d0) then
#if !STANDALONE           
        write(warnline,*) 'HelmEOS: Y_e out of range!! Y_e=:', ye
        anyerr = 1
        NSEkeyerr(j) = 2
        !call CCTK_ERROR(warnline)
#else           
        write(*,*) 'HelmEOS: Y_e out of range!! Y_e=:', ye
        anyerr = 1
        NSEkeyerr(j) = 2
        STOP
#endif    
        end if
        
        xx_atm = 0.0d0       
        xxa_guess = ye
        ! calculate NSE abundances
        call nse_abundance(den, btemp, ye, xxn, xxp, xxa, xxa_guess, NSEconverge)
        
        if (xxn.lt.0.0d0 .or. xxp .lt. 0.0d0 .or. xxa .lt. 0.0d0 .or. &
            xxn.gt.1.0d0 .or. xxp .gt. 1.0d0 .or. xxa .gt. 1.0d0 .or. NSEconverge .eq. 0) then
          xxa_guess = 1.0d0 - ye
          call nse_abundance(den, btemp, ye, xxn, xxp, xxa, xxa_guess, NSEconverge)
        end if
        
        call nse_derivatives(den,btemp,xxn,xxp,xxa, xxa_r,xxa_t,xxa_y,xxn_r,xxn_t,xxn_y,xxp_r,xxp_t,xxp_y)

        if (xxn.lt.0.0d0 .or. xxp .lt. 0.0d0 .or. xxa .lt. 0.0d0 .or. &
            xxn.gt.1.0d0 .or. xxp .gt. 1.0d0 .or. xxa .gt. 1.0d0 .or. NSEconverge .eq. 0) then
#if !STANDALONE           
            write(warnline,*) 'NSE calc failed: den, temp, ye, conv, xn, xp, xa:', &
              denRow(j), tempRow(j), ye, NSEconverge, xxn, xxp, xxa
            call CCTK_WARN(1,warnline)
#else
        
           write(*,*) 'HelmEOS: NSE calculations failed'
           write(*,*) 'den, temp, ye =', den, btemp, ye
           write(*,*) 'xn, xp, xa =', xxn, xxp, xxa
           write(*,*) 'conv =', NSEconverge
           write(*,*) 'ERROR: eos_helm: wrong abundances out of NSE!!!'

#endif           
        endif

        !Prevent anomalous abundances (for testing purposes)
        xxn = max(xxn,1d-10)
        !xxn = min(xxn,1d0)
        xxp = max(xxp,1d-10)
        !xxp = min(xxp,1d0)
        xxa = max(xxa,1d-10)
        !xxa = min(xxa,1d0)
        
        !abar  = 1d0/(xx_atm + xxn + xxp + xxa/4d0)
        abar  = 1d0/(xxn + xxp + xxa/4d0)
        zbar  = abar*ye
        ytot1 = 1d0/abar

     ! else
!         ! atmosphere
!         ! assume pure hydrogen
!         ! assume temperature > few in 1e4 K, so that hydrogen is completely ionized
!         
!         !ensure density and temperature are reset to atmosphere (in cgs units)
!         btemp  = helmeos_temp_atmo * inv_temp_gf
!         den    = helmeos_rho_atmo * inv_rho_gf
!         
!     
!         xxn = 0.0d0 !xxnRow(j)
!         xxp = 1.0d0 !xxpRow(j)
!         xxa = 0.0d0 !xxaRow(j)
!         ye  = 1.0d0
!         
!         xx_atm = 0.
!         NSEconverge = 1
! 
!         xxa_r=0.; xxa_t=0.; xxa_y=0.; xxn_r=0.; xxn_t=0.; xxn_y=0.; xxp_r=0.; xxp_t=0.; xxp_y=0.;
!  
!         abar  = 1d0
!         zbar  = 1d0
!         ytot1 = 1d0
!  
!      endif

     !  frequent combinations
     deni    = 1.0d0/den
     tempi   = 1.0d0/btemp 
     kt      = kBerg * btemp
     ktinv   = 1.0d0/kt
     kavoy   = kergavo * ytot1

     
     ! ************************************
     !  RADIATION SECTION:
     ! ************************************
     
     prad    = asoli3 * btemp * btemp * btemp * btemp
     dpraddt = 4.0d0 * prad * tempi
     dpraddd = 0.0d0

     x1      = prad * deni 
     erad    = 3.0d0 * x1
     deraddd = -erad*deni
     deraddt = 4.0d0 * erad * tempi

     srad    = (x1 + erad)*tempi
     dsraddd = (dpraddd*deni - x1*deni + deraddd)*tempi
     dsraddt = (dpraddt*deni + deraddt - srad)*tempi


     !  ION SECTION:
     
     !  RF: modified for NSE
     dxnidd  = avo * ytot1
     xni     = dxnidd * den

     pion    = xni * kt
     !RF:
     !dpiondd = avo * ytot1 * kt
     !dpiondt = xni * kBerg
     dpiondd = dxnidd * kt + den*kt*avo*(xxn_r + xxp_r + xxa_r/4d0)
     dpiondt = xni * kBerg  + den*kt*avo*(xxn_t + xxp_t + xxa_t/4d0)
     
   
     !RF:
     !eion    = 1.5e0 * pion * deni
     !deiondd = (1.5e0 * dpiondd - eion)*deni
     !deiondt = 1.5e0 * xni * kBerg *deni
     eion    = 1.5d0 * pion*deni - (Qalpha/4d0)*avo*xxa
     deiondd =                        avo*kt*(1.5d0*xxn_r + 1.5d0*xxp_r + xxa_r/4d0*(1.5d0-Qalpha/kt))
     deiondt = 1.5d0*avo*ytot1*kBerg + avo*kt*(1.5d0*xxn_t + 1.5d0*xxp_t + xxa_t/4d0*(1.5d0-Qalpha/kt))
     

     !!  sackur-tetrode equation for the ion entropy of 
     !!  a single ideal gas characterized by abar
     x2      = abar*abar*sqrt(abar) * deni*avoinv        !x
     y0      = sioncon * btemp                           !s
     z0      = x2 * y0 * sqrt(y0)                        !z
     !RF:
     !sion    = (pion*deni + eion)*tempi + kavoy*log(z0)
     !dsiondd = (dpiondd*deni - pion*deni*deni + deiondd)*tempi    &
     !     - kavoy * deni
     !dsiondt = (dpiondt*deni + deiondt)*tempi   &
     !     - (pion*deni + eion) * tempi*tempi  &
     !     + 1.5e0 * kavoy * tempi
     sion    = kergavo*ytot1*(2.5d0 + log(z0))
     dadd    = -abar*abar*(xxn_r + xxp_r + xxa_r/4d0)
     dadt    = -abar*abar*(xxn_t + xxp_t + xxa_t/4d0)
     dsiondd = -kergavo*ytot1*deni       + (2.5d0*kergavo*ytot1-sion)*ytot1*dadd
     dsiondt = 1.5d0*kergavo*ytot1*tempi + (2.5d0*kergavo*ytot1-sion)*ytot1*dadt

     ! normalized (dimensionless) chemical potentials of neutrons and protons: eta = mu / kT
     ! dissociation energies chi=0 for n and p, as the nucleon state is taken
     ! to be the zero point of energy (cf. definition of eion)
     etap = log( xxp * den * avo * hplanck**3d0 / ( 2.0d0 * pi * amu * kt )**1.5d0 )
     etan = log( xxn * den * avo * hplanck**3d0 / ( 2.0d0 * pi * amu * kt )**1.5d0 )

     
     ! ************************************
     !  ELECTRON-POSITRON SECTION:
     ! ************************************
     
     ! enter the table with ye*den, no checks of the input
     din = ye*den

     ! check whether we are in the validity range of the EOS table
     if (.not. ((din .lt. 1d+11) .and. (btemp .lt. 1d+11)) ) then
        
        ! temperature and/or density limits of the table are exceeded,
        ! use analytic formulas for ultr-relativistic Fermi gas
     
        !normalized fermi momentum
        x3     = 1.05457266d-27*2.9979d+10*(3*acos(-1d0)**2*din*avo)**(1d0/3d0)/kt
        !invert cubic equation
        !etaele = (0.5d0*x3**3d0 + sqrt(0.25d0*x3**6.0d0 + acos(-1d0)**6d0/27d0))**(1d0/3d0) &
        !       +   (0.5d0*x3**3d0 - sqrt(0.25d0*x3**6.0d0 + acos(-1d0)**6d0/27d0))**(1d0/3d0)
        !DS: corrected, as above expression leads to NaNs
        etaele = (0.5d0*x3**3d0 + sqrt(0.25d0*x3**6.0d0 + acos(-1d0)**6d0/27d0))**(1d0/3d0) &
          -   abs(0.5d0*x3**3d0 - sqrt(0.25d0*x3**6.0d0 + acos(-1d0)**6d0/27d0))**(1d0/3d0)       
        !----
        pele   = asoli3*btemp**4d0*(7d0/4d0 + 7.5d0*etaele**2/acos(-1d0)**2* &
                                  (1d0 + 0.5d0*etaele**2/acos(-1d0)**2))
                                 
        eele   = 3d0*pele*deni
        sele   = acos(-1d0)**2*ye*avo*kBerg*(7d0/15d0*acos(-1d0)**2+etaele**2)/etaele/(etaele**2+acos(-1d0)**2)

        !deta/dt
        x1     = -3d0/btemp/(1d0/etaele + 2d0*etaele/(etaele**2+acos(-1d0)**2))
        !deta/drho 
        x3     = deni/(1d0/etaele + 2d0*etaele/(etaele**2+acos(-1d0)**2))

        dpepdt = 4d0*pele/btemp + x1*15d0/acos(-1d0)**2*asoli3*btemp**4d0* &
                                   etaele*(1d0+etaele**2/acos(-1d0)**2)

        dpepdd =                x3*15d0/acos(-1d0)**2*asoli3*btemp**4d0* &
                                   etaele*(1d0+etaele**2/acos(-1d0)**2)

        deepdt =              3d0*dpepdt*deni
        deepdd = -eele*deni + 3d0*dpepdd*deni

        dsepdt = x1*sele*(2d0*etaele/(7d0/15d0*acos(-1d0)**2+etaele**2) - &
                          2d0*etaele/(acos(-1d0)**2+etaele**2) -1d0/etaele)
        dsepdd = x3*sele*(2d0*etaele/(7d0/15d0*acos(-1d0)**2+etaele**2) - &
                          2d0*etaele/(acos(-1d0)**2+etaele**2) -1d0/etaele)      

     else

       ! we are inside the validity range of the table, so use the EOS table

       !!  hash locate this temperature and density
       jat = int((log10(btemp) - heos_tlo)*heos_tstpi) + 1
       jat = max(1,min(jat,HEOSJMAX-1))
       iat = int((log10(din) - heos_dlo)*heos_dstpi) + 1
       iat = max(1,min(iat,HEOSIMAX-1))

       !!  access the table locations only once
       fi(1)  = heos_f(iat,jat)
       fi(2)  = heos_f(iat+1,jat)
       fi(3)  = heos_f(iat,jat+1)
       fi(4)  = heos_f(iat+1,jat+1)
       fi(5)  = heos_ft(iat,jat)
       fi(6)  = heos_ft(iat+1,jat)
       fi(7)  = heos_ft(iat,jat+1)
       fi(8)  = heos_ft(iat+1,jat+1)
       fi(9)  = heos_ftt(iat,jat)
       fi(10) = heos_ftt(iat+1,jat)
       fi(11) = heos_ftt(iat,jat+1)
       fi(12) = heos_ftt(iat+1,jat+1)
       fi(13) = heos_fd(iat,jat)
       fi(14) = heos_fd(iat+1,jat)
       fi(15) = heos_fd(iat,jat+1)
       fi(16) = heos_fd(iat+1,jat+1)
       fi(17) = heos_fdd(iat,jat)
       fi(18) = heos_fdd(iat+1,jat)
       fi(19) = heos_fdd(iat,jat+1)
       fi(20) = heos_fdd(iat+1,jat+1)
       fi(21) = heos_fdt(iat,jat)
       fi(22) = heos_fdt(iat+1,jat)
       fi(23) = heos_fdt(iat,jat+1)
       fi(24) = heos_fdt(iat+1,jat+1)
       fi(25) = heos_fddt(iat,jat)
       fi(26) = heos_fddt(iat+1,jat)
       fi(27) = heos_fddt(iat,jat+1)
       fi(28) = heos_fddt(iat+1,jat+1)
       fi(29) = heos_fdtt(iat,jat)
       fi(30) = heos_fdtt(iat+1,jat)
       fi(31) = heos_fdtt(iat,jat+1)
       fi(32) = heos_fdtt(iat+1,jat+1)
       fi(33) = heos_fddtt(iat,jat)
       fi(34) = heos_fddtt(iat+1,jat)
       fi(35) = heos_fddtt(iat,jat+1)
       fi(36) = heos_fddtt(iat+1,jat+1)


       !!  various differences
       xt  = max( (btemp - heos_t(jat))*heos_dtInv(jat), 0.0d0)
       xd  = max( (din - heos_d(iat))*heos_ddInv(iat), 0.0d0)
       mxt = 1.0d0 - xt
       mxd = 1.0d0 - xd


       !!  the density and temperature basis functions
       si0t =   psi0(xt)
       si1t =   psi1(xt)*heos_dt(jat)
       si2t =   psi2(xt)*heos_dtSqr(jat)

       si0mt =  psi0(mxt)
       si1mt = -psi1(mxt)*heos_dt(jat)
       si2mt =  psi2(mxt)*heos_dtSqr(jat)

       si0d =   psi0(xd)
       si1d =   psi1(xd)*heos_dd(iat)
       si2d =   psi2(xd)*heos_ddSqr(iat)

       si0md =  psi0(mxd)
       si1md = -psi1(mxd)*heos_dd(iat)
       si2md =  psi2(mxd)*heos_ddSqr(iat)


       !!  the first derivatives of the basis functions
       dsi0t =   dpsi0(xt)*heos_dtInv(jat)
       dsi1t =   dpsi1(xt)
       dsi2t =   dpsi2(xt)*heos_dt(jat)

       dsi0mt = -dpsi0(mxt)*heos_dtInv(jat)
       dsi1mt =  dpsi1(mxt)
       dsi2mt = -dpsi2(mxt)*heos_dt(jat)

       dsi0d =   dpsi0(xd)*heos_ddInv(iat)
       dsi1d =   dpsi1(xd)
       dsi2d =   dpsi2(xd)*heos_dd(iat)

       dsi0md = -dpsi0(mxd)*heos_ddInv(iat)
       dsi1md =  dpsi1(mxd)
       dsi2md = -dpsi2(mxd)*heos_dd(iat)


       !!  the second derivatives of the basis functions
       ddsi0t =   ddpsi0(xt)*heos_dtSqrInv(jat)
       ddsi1t =   ddpsi1(xt)*heos_dtInv(jat)
       ddsi2t =   ddpsi2(xt)

       ddsi0mt =  ddpsi0(mxt)*heos_dtSqrInv(jat)
       ddsi1mt = -ddpsi1(mxt)*heos_dtInv(jat)
       ddsi2mt =  ddpsi2(mxt)


       !!  the free energy
       free  = h5( & 
          si0t,   si1t,   si2t,   si0mt,   si1mt,   si2mt, & 
          si0d,   si1d,   si2d,   si0md,   si1md,   si2md)


       !!  derivative with respect to density
       df_d  = h5( & 
          si0t,   si1t,   si2t,   si0mt,   si1mt,   si2mt, & 
          dsi0d,  dsi1d,  dsi2d,  dsi0md,  dsi1md,  dsi2md)


       !!  derivative with respect to temperature
       df_t = h5( & 
          dsi0t,  dsi1t,  dsi2t,  dsi0mt,  dsi1mt,  dsi2mt, & 
          si0d,   si1d,   si2d,   si0md,   si1md,   si2md)


       !!  second derivative with respect to temperature
       df_tt = h5( & 
          ddsi0t, ddsi1t, ddsi2t, ddsi0mt, ddsi1mt, ddsi2mt, & 
          si0d,   si1d,   si2d,   si0md,   si1md,   si2md)


       !!  second derivative with respect to temperature and density
       df_dt = h5( & 
          dsi0t,  dsi1t,  dsi2t,  dsi0mt,  dsi1mt,  dsi2mt, & 
          dsi0d,  dsi1d,  dsi2d,  dsi0md,  dsi1md,  dsi2md)


       !!  now get the pressure derivative with density, chemical potential, and 
       !!  electron positron number densities
       !!  get the interpolation weight functions
       si0t   =  xpsi0(xt)
       si1t   =  xpsi1(xt)*heos_dt(jat)

       si0mt  =  xpsi0(mxt)
       si1mt  =  -xpsi1(mxt)*heos_dt(jat)

       si0d   =  xpsi0(xd)
       si1d   =  xpsi1(xd)*heos_dd(iat)

       si0md  =  xpsi0(mxd)
       si1md  =  -xpsi1(mxd)*heos_dd(iat)


       !!  pressure derivative with density
       dpepdd  = h3dpd(iat,jat, & 
          si0t,   si1t,   si0mt,   si1mt, & 
          si0d,   si1d,   si0md,   si1md)
       dpepdd  = max(ye * dpepdd,0.0d0)



       !!  electron chemical potential etaele
       etaele  = h3e(iat,jat, & 
                    si0t,   si1t,   si0mt,   si1mt, & 
                    si0d,   si1d,   si0md,   si1md)


       !!  electron + positron number densities
       xnefer   = h3x(iat,jat, & 
                  si0t,   si1t,   si0mt,   si1mt, & 
                  si0d,   si1d,   si0md,   si1md)


       !!  the desired electron-positron thermodynamic quantities
       x3      = din * din
       pele    = x3 * df_d
       dpepdt  = x3 * df_dt

       sele    = -df_t * ye
       dsepdt  = -df_tt * ye
       dsepdd  = -df_dt * ye * ye

 
       eele    = ye * free + btemp * sele
       deepdt  = btemp * dsepdt
       deepdd  = ye*ye*df_d + btemp*dsepdd

     endif
     
     ! ************************************
     !  COULOMB SECTION:
     ! ************************************
     
     !  initialize
     pcoul    = 0.0d0
     dpcouldd = 0.0d0
     dpcouldt = 0.0d0
     ecoul    = 0.0d0
     decouldd = 0.0d0
     decouldt = 0.0d0
     scoul    = 0.0d0
     dscouldd = 0.0d0
     dscouldt = 0.0d0

     !! Set the coulomb multiplier to a local value -- we might change it only within this call
     local_coulombMult = heos_coulombMult

     !!  uniform background corrections & only the needed parts for speed
     !!  plasg is the plasma coupling parameter
     !!  split up calculations below -- they all used to depend upon a redefined z
     z1        = forth * pi
     s1        = z1 * xni
     dsdd      = z1 * dxnidd
     lami      = 1.0d0/s1**third
     inv_lami  = 1.0d0/lami
     z2        = -third * lami/s1
     lamidd    = z2 * dsdd


     plasg     = zbar*zbar*esqu*ktinv*inv_lami
     plasg_inv = 1.0d0/plasg  
     z3        = -plasg * inv_lami
     plasgdd   = z3 * lamidd
     plasgdt   = -plasg * ktinv * kBerg


     !!  yakovlev & shalybkov 1989 equations 82, 85, 86, 87
     if (plasg .ge. 1.0) then
        x4       = plasg**(0.25d0)
        z4       = c1/x4
        ecoul    = dxnidd * kt * (a1*plasg + b1*x4 + z4 + d1cc)
        pcoul    = third * den * ecoul
        scoul    = -kavoy*(3.0d0*b1*x4 - 5.0d0*z4 &
             + d1cc*(log(plasg) - 1.0d0) - e1cc)

        y1       = dxnidd * kt * (a1 + 0.25d0*plasg_inv*(b1*x4 - z4))
        decouldd = y1 * plasgdd 
        decouldt = y1 * plasgdt + ecoul * tempi
        dpcouldd = third * (ecoul + den * decouldd)
        dpcouldt = third * den  * decouldt

        y2       = -kavoy*plasg_inv*(0.75d0*b1*x4 + 1.25d0*z4 + d1cc)
        dscouldd = y2 * plasgdd
        dscouldt = y2 * plasgdt



        !!  yakovlev & shalybkov 1989 equations 102, 103, 104
     else if (plasg .lt. 1.0) then
        x5       = plasg * sqrt(plasg)
        y3       = plasg**b2
        z5       = c2 * x5 - third * a2 * y3
        pcoul    = -pion * z5
        ecoul    = 3.0d0 * pcoul * deni
        scoul    = -kavoy*(c2*x5 - a2*(b2 - 1.0d0)/b2*y3)

        s2       = (1.5d0*c2*x5 - third*a2*b2*y3)*plasg_inv
        dpcouldd = -dpiondd*z5 - pion*s2*plasgdd
        dpcouldt = -dpiondt*z5 - pion*s2*plasgdt
        decouldd = 3.0d0*dpcouldd*deni - ecoul*deni
        decouldt = 3.0d0*dpcouldt*deni

        !! Equations for decoulda and decouldz can be added here....

        s3       = -kavoy*plasg_inv*(1.5d0*c2*x5 - a2*(b2 - 1.0d0)*y3)
        dscouldd = s3 * plasgdd
        dscouldt = s3 * plasgdt
     end if

     s4 = prad + pion + pele
     x6 = s4 + pcoul*heos_coulombMult

     if ( x6 <= 0.e0 ) then

       write(*,*) 'ERROR: Helmholtz EOS: Negative total pressure...'
       
       if ( s4 > 0.e0 ) then
       
         write(*,*) '...caused by Coulomb correction!'
         write(*,*) '...fixing problem by switching off Coulomb correction!'
         local_coulombMult = 0.e0
       
       else
         write(*,*) '...NOT caused by Coulomb correction!'
       
       endif
     
     endif  
    

     pcoul    = pcoul * local_coulombMult
     dpcouldd = dpcouldd * local_coulombMult
     dpcouldt = dpcouldt * local_coulombMult

     ecoul    = ecoul * local_coulombMult
     decouldd = decouldd * local_coulombMult
     decouldt = decouldt * local_coulombMult

     scoul    = scoul * local_coulombMult
     dscouldd = dscouldd * local_coulombMult 
     dscouldt = dscouldt * local_coulombMult

     
     ! ************************************
     !  TOTAL QUANTITIES (ALL SECTIONS ADDED)
     ! ************************************
     
     !  sum all the components
     pres    = prad    + pion    + pele   + pcoul
     ener    = erad    + eion    + eele   + ecoul
     entr    = srad    + sion    + sele   + scoul

     dpresdd = dpraddd + dpiondd + dpepdd + dpcouldd
     dpresdt = dpraddt + dpiondt + dpepdt + dpcouldt

     denerdd = deraddd + deiondd + deepdd + decouldd
     denerdt = deraddt + deiondt + deepdt + decouldt

     dentrdd = dsraddd + dsiondd + dsepdd + dscouldd
     dentrdt = dsraddt + dsiondt + dsepdt + dscouldt

     !  compute gamma_1, gamma_2, gamma_3, cv, cp, soundspeed
     presi = 1.0d0/pres
     chit  = btemp * presi * dpresdt
     chid  = dpresdd * den * presi
     x7    = pres * deni * chit/(btemp * denerdt)
     gam1  = chit * x7 + chid
     
     gam3  = x7 + 1.0d0
     nabad = x7/gam1
     gam2  = 1.0d0/(1.0d0 - nabad)
     cv    = denerdt
     cp    = cv * gam1 / chid
     zz    = 1.0d0 + (ener + v_c*v_c) * den * presi
     sound = v_c * sqrt( gam1 / zz )
     soundCU = sqrt( gam1 / zz )
     
     ! further derivatives
     dptotdetot_d = dpresdt / denerdt  ! dp/de at const. density
     
     ! define the output
     
     ! convert to Cactus units and store the output
     ptotRow(j)   = pres * press_gf   ! pressure
     etotRow(j)   = ener * eps_gf     ! specific internal energy
     stotRow(j)   = entr * conv_s_CGS2CU  ! specific entropy
     htotRow(j)   = (ener + pres*deni) * eps_gf ! specific enthalpy, this is the Newtonian definition without rest mass

     ! thermodynamic derivatives in CGS units 
     dpdRow(j)    = dpresdd * press_gf * inv_rho_gf
     dptRow(j)    = dpresdt * press_gf * inv_temp_gf

     dedRow(j)    = denerdd * eps_gf * inv_rho_gf
     detRow(j)    = denerdt * eps_gf * inv_temp_gf

     dsdRow(j)    = dentrdd * conv_s_CGS2CU * inv_rho_gf
     dstRow(j)    = dentrdt * conv_s_CGS2CU * inv_temp_gf
     
     dhdRow(j)    = (denerdd + dpresdd * deni - pres * deni*deni) * conv_h_CGS2CU * inv_rho_gf
     dhtRow(j)    = (denerdt + dpresdt * deni) * conv_h_CGS2CU * inv_temp_gf
     
     !UNUSED     pradRow(j)   = prad
     !UNUSED     eradRow(j)   = erad
     !UNUSED     sradRow(j)   = srad

     !UNUSED     pionRow(j)   = pion
     !UNUSED     eionRow(j)   = eion
     !UNUSED     sionRow(j)   = sion
 
     pelRow(j)   = pele * press_gf 
     !UNUSED     eeleRow(j)   = eele
     !UNUSED     seleRow(j)   = sele

     neRow(j)    = xnefer
     etaeRow(j)  = etaele
     etapRow(j)  = etap
     etanRow(j)  = etan

     gam1Row(j)   = gam1 

     cvRow(j)     = cv
     cpRow(j)     = cp
     cs2Row(j)    = soundCU * soundCU  ! sound speed squared in Cactus code units (in units of v_c)
     
     dpde_d_Row(j) = dptotdetot_d * press_gf * inv_eps_gf ! dp/de at const. density

     
     ! abundances
     AbarRow(j) = abar
     xxnRow(j)  = xxn
     xxpRow(j)  = xxp
     xxaRow(j)  = xxa  
     xxatmRow(j) = xx_atm
     
     if (NSEconverge .eq. 0) then
       NSEkeyerr(j) = 1 
       anyerr = 1
     !else
     !  NSEkeyerr(j) = 0
     end if
     
     ! end of main loop over npoints
  enddo

  return

end subroutine helm_eos


!----------------------------------------------------------------
!
! Some subroutines and functions for computing the NSE abundances
!
!----------------------------------------------------------------


subroutine nse_abundance(dens,temp,ye, xn,xp,xa, xa_guess, conv)
  
  ! written by Rodrigo Fernandez, modified by Daniel Siegel

  use EOS_Omni_Module
  
  implicit none

  double precision, intent(IN) :: dens, temp, ye, xa_guess
  double precision, intent(OUT):: xn, xp, xa
  integer, intent(OUT) :: conv

  double precision, parameter :: hbar = 1.0545727d-27
  integer, parameter :: itmax = 100

  integer :: i,j
  double precision :: x_n, x_p, n_Q, n
  double precision :: Fa, dFa, xa_new

  n   = dens/amu  ! the baryon mass is the atomic mass for this EOS
  n_Q = (amu*kBerg*temp/2d0/acos(-1d0)/hbar**2)**1.5d0

  xa = xa_guess
  do i=1,itmax
    Fa  = x_n(ye,xa)**2*x_p(ye,xa)**2 - 0.5d0*xa*(n_Q/n)**3*exp(-Qalpha/kBerg/temp)
    dFa = -x_p(ye,xa)*x_n(ye,xa)**2 - x_n(ye,xa)*x_p(ye,xa)**2 - 0.5d0*(n_Q/n)**3*exp(-Qalpha/kBerg/temp)
    xa_new = abs(xa - Fa/dFa)
    xa_new = min(xa_new, 1.0d0 / xa_new)
    
    xa = xa_new
    
    if (abs(Fa/dFa/xa).lt.(1d-10)) then
      conv = 1
      exit
    elseif (i .eq. itmax) then
      conv = 0
      !write(*,*)
      !write(*,*) 'WARNING: Helm EOS: NSE abundances did not converge!'
      !write(*,*) ' density =', dens
      !write(*,*) ' temp =', temp
      !write(*,*) ' Ye =', ye
      !write(*,*) ' current error =', Fa/dFa
      !write(*,*)
    end if
   
  enddo

  xn = x_n(ye,xa)
  xp = x_p(ye,xa)

  return
  
end subroutine nse_abundance



function x_n(ye,xa)
  
  implicit none

  double precision, intent(IN) :: ye, xa
  double precision ::x_n

  x_n = 1d0 - ye - 0.5d0*xa

  return

end function x_n
 
 

function x_p(ye,xa)

  implicit none
  
  double precision, intent(IN) :: ye, xa
  double precision :: x_p

  x_p = ye - 0.5d0*xa

  return

end function x_p



!---------------------------------------------------------------------------------
!Derivative of abundances in NSE.

subroutine nse_derivatives(dens,temp,xn,xp,xa, xa_r,xa_t,xa_y,xn_r,xn_t,xn_y,xp_r,xp_t,xp_y)

  ! written by Rodrigo Fernandez

  use EOS_Omni_Module

  implicit none

  double precision, intent(IN)  :: dens, temp, xn, xp, xa
  double precision, intent(OUT) :: xa_r,xa_t,xa_y,xn_r,xn_t,xn_y,xp_r,xp_t,xp_y

  xa_r = (3d0/dens)*xn*xp*xa/(xn*xp + xn*xa + xp*xa)
  xa_t = -1d0/temp*(4.5d0 + Qalpha/kBerg/temp)*xn*xp*xa/(xn*xp + xn*xa + xp*xa)
  xa_y = 2d0*(xn-xp)*xa/(xn*xp + xn*xa + xp*xa)

  xn_r = -0.5d0*xa_r
  xn_t = -0.5d0*xa_t
  xn_y = -1d0 -0.5d0*xa_y

  xp_r = xn_r
  xp_t = xn_t
  xp_y = 1d0 -0.5d0*xa_y

  return

end subroutine nse_derivatives
