From e07142cccd9bc43a37c41861d418c58fb38506ea Mon Sep 17 00:00:00 2001
From: Austin Harris <jaharris87@gmail.com>
Date: Mon, 1 Jul 2024 14:57:23 -0400
Subject: [PATCH] scalar/vector interfaces for other EoS interface modules

---
 source/xnet_eos_bahcall.f90    | 319 ++++++++++++++++++++++++++++-----
 source/xnet_eos_helm.f90       | 297 +++++++++++++++++++++++++-----
 source/xnet_eos_starkiller.f90 |  10 +-
 3 files changed, 528 insertions(+), 98 deletions(-)

diff --git a/source/xnet_eos_bahcall.f90 b/source/xnet_eos_bahcall.f90
index 88a30e7a..b6331f96 100644
--- a/source/xnet_eos_bahcall.f90
+++ b/source/xnet_eos_bahcall.f90
@@ -6,7 +6,24 @@
 !***************************************************************************************************
 
 Module xnet_eos
+  Use xnet_types, Only: dp
   Implicit None
+  Real(dp), Allocatable :: ye(:), ytot(:), abar(:), zbar(:), z2bar(:), zibar(:), sratio(:)
+
+  Interface eos_interface
+    Module Procedure eos_interface_scalar
+    Module Procedure eos_interface_vector
+  End Interface
+
+  Interface eos_screen
+    Module Procedure eos_screen_scalar
+    Module Procedure eos_screen_vector
+  End Interface eos_screen
+
+  Interface salpeter_ratio
+    Module Procedure salpeter_ratio_scalar
+    Module Procedure salpeter_ratio_vector
+  End Interface
 
 Contains
 
@@ -14,49 +31,53 @@ Subroutine eos_initialize
     !-----------------------------------------------------------------------------------------------
     ! This routine initializes the EoS.
     !-----------------------------------------------------------------------------------------------
+    Use xnet_controls, Only: nzevolve
     Implicit None
+
+    Allocate (ye(nzevolve))
+    Allocate (ytot(nzevolve))
+    Allocate (abar(nzevolve))
+    Allocate (zbar(nzevolve))
+    Allocate (z2bar(nzevolve))
+    Allocate (zibar(nzevolve))
+    Allocate (sratio(nzevolve))
+    ye = 0.0
+    ytot = 0.0
+    abar = 0.0
+    zbar = 0.0
+    z2bar = 0.0
+    zibar = 0.0
+    sratio = 0.0
+
     Return
   End Subroutine eos_initialize
 
-  Subroutine eos_interface(rho,t9,y,ye,cv,efermkt,defermktdt9,xext,aext,zext)
+  Subroutine eosx(t9,rho,ye,ytot,abar,zbar,z53bar,cv,efermkt,defermktdt9)
     !-----------------------------------------------------------------------------------------------
-    ! This routine updates the Equation of State for changes in temperature and density.
-    !
-    ! Bahcall introduced a fit for the plasma parameters needed to calculate the screening factors
-    ! for nuclear reaction based on the easier to calculate efermkt (the ratio of the electron
-    ! Fermi Energy to the themal energy) instead of the more general expresions based on etae (the
-    ! ratio of the electron chemical potential to the thermal energy), which generally require a
-    ! complex Equation of State to calculate.   For basic cases, XNet uses Bahcall's approach
-    ! to simplify the
-    !
-    ! Note: the factors efermkt and defermktdt9 (the derivative of efermkt wrt to temperature) in the
-    ! subroutine interface are stored in XNet's conditions module as etae and detaedt9.  For more
-    ! general EoS, etae and detaedt9 are their proper selves.
+    ! This routine interfaces with and calls the underlying EoS.
     !-----------------------------------------------------------------------------------------------
-    Use nuclear_data, Only: ny, zz53
-    Use xnet_abundances, Only: y_moment
     Use xnet_constants, Only: asig, avn, bok, e2, pi, pi2, third, emass, clt, ele_en, hbar
-    Use xnet_controls, Only: idiag, iheat, iscrn, lun_diag
+    Use xnet_controls, Only: iheat, iscrn
     Use xnet_types, Only: dp
     Implicit None
 
     ! Input variables
-    Real(dp), Intent(in) :: t9, rho, y(ny), xext, aext, zext
+    Real(dp), Intent(in) :: t9, rho, ye, ytot, abar, zbar, z53bar
 
-    ! Output variables
-    Real(dp), Intent(out) :: ye, cv, efermkt, defermktdt9
+    ! Ouput variables
+    Real(dp), Intent(out) :: cv, efermkt, defermktdt9
 
     ! Local variables
     Real(dp), Parameter :: a1 = 0.898004, b1 = 0.96786, c1 = 0.220703, d1 = -0.86097
     Real(dp), Parameter :: a2 = -0.86602540378, b2 = 0.29561, c2 = 1.9885
-    Real(dp) :: ytot, bkt, abar, zbar, z2bar, zibar, z53bar
+    Real(dp) :: bkt
     Real(dp) :: etae_mb, ae, gam, gam14, gam32, gamc2, rel_ef
     Real(dp) :: eion, erad, ecoul
     Real(dp) :: deiondt9, deraddt9, decouldt9
 
-    ! Calculate Ye and other needed moments of the abundance distribution
-    Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext)
-
+    cv = 0.0
+    efermkt = 0.0
+    defermktdt9 = 0.0
     If ( iscrn > 0 .or. iheat > 0 ) Then
 
       bkt = bok*t9
@@ -73,7 +94,6 @@ Subroutine eos_interface(rho,t9,y,ye,cv,efermkt,defermktdt9,xext,aext,zext)
       deraddt9 = 4.0*erad / t9
 
       ae = (3.0 / (4.0*pi*avn*rho*ye))**third ! electron-sphere radius
-      z53bar = sum(zz53 * y) / ytot
       gam = z53bar*e2 / (ae*bkt) ! ionic Coulomb coupling parameter
       If ( gam >= 1.0 ) Then
         gam14 = gam**0.25
@@ -86,18 +106,122 @@ Subroutine eos_interface(rho,t9,y,ye,cv,efermkt,defermktdt9,xext,aext,zext)
         decouldt9 = ecoul/t9 - bok * ytot * (1.5*a2*gam32 + b2*c2*gamc2)
       EndIf
       cv = deiondt9 + deraddt9 + decouldt9
-    Else
-      efermkt = 0.0
-      defermktdt9 = 0.0
-      cv = 0.0
     EndIf
 
+  End Subroutine eosx
+
+  Subroutine eos_interface_scalar(rho,t9,y,ye,cv,efermkt,defermktdt9,xext,aext,zext)
+    !-----------------------------------------------------------------------------------------------
+    ! This routine updates the Equation of State for changes in temperature and density.
+    !
+    ! Bahcall introduced a fit for the plasma parameters needed to calculate the screening factors
+    ! for nuclear reaction based on the easier to calculate efermkt (the ratio of the electron
+    ! Fermi Energy to the themal energy) instead of the more general expresions based on etae (the
+    ! ratio of the electron chemical potential to the thermal energy), which generally require a
+    ! complex Equation of State to calculate.   For basic cases, XNet uses Bahcall's approach
+    ! to simplify the
+    !
+    ! Note: the factors efermkt and defermktdt9 (the derivative of efermkt wrt to temperature) in the
+    ! subroutine interface are stored in XNet's conditions module as etae and detaedt9.  For more
+    ! general EoS, etae and detaedt9 are their proper selves.
+    !-----------------------------------------------------------------------------------------------
+    Use nuclear_data, Only: ny, zz53
+    Use xnet_abundances, Only: y_moment
+    Use xnet_controls, Only: idiag, lun_diag
+    Use xnet_types, Only: dp
+    Implicit None
+
+    ! Input variables
+    Real(dp), Intent(in) :: t9, rho, y(ny), xext, aext, zext
+
+    ! Output variables
+    Real(dp), Intent(out) :: ye, cv, efermkt, defermktdt9
+
+    ! Local variables
+    Real(dp) :: ytot, abar, zbar, z2bar, zibar, z53bar
+
+    ! Calculate Ye
+    Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext)
+    z53bar = sum(zz53 * y) / ytot
+
+    ! Call the eos
+    Call eosx(t9,rho,ye,ytot,abar,zbar,z53bar,cv,efermkt,defermktdt9)
+
     If ( idiag >= 3 ) Write(lun_diag,"(a,6es23.15)") 'EOS',t9,rho,ye,cv,efermkt,defermktdt9
 
     Return
-  End Subroutine eos_interface
+  End Subroutine eos_interface_scalar
 
-  Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext)
+  Subroutine eos_interface_vector(t9,rho,y,ye,cv,efermkt,defermktdt9,xext,aext,zext,mask_in)
+    !-----------------------------------------------------------------------------------------------
+    ! This routine updates the Equation of State for changes in temperature and density.
+    !
+    ! Bahcall introduced a fit for the plasma parameters needed to calculate the screening factors
+    ! for nuclear reaction based on the easier to calculate efermkt (the ratio of the electron
+    ! Fermi Energy to the themal energy) instead of the more general expresions based on etae (the
+    ! ratio of the electron chemical potential to the thermal energy), which generally require a
+    ! complex Equation of State to calculate.   For basic cases, XNet uses Bahcall's approach
+    ! to simplify the
+    !
+    ! Note: the factors efermkt and defermktdt9 (the derivative of efermkt wrt to temperature) in the
+    ! subroutine interface are stored in XNet's conditions module as etae and detaedt9.  For more
+    ! general EoS, etae and detaedt9 are their proper selves.
+    !-----------------------------------------------------------------------------------------------
+    Use nuclear_data, Only: ny, zz53
+    Use xnet_abundances, Only: y_moment
+    Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive
+    Use xnet_types, Only: dp
+    Implicit None
+
+    ! Input variables
+    Real(dp), Intent(in) :: t9(zb_lo:zb_hi), rho(zb_lo:zb_hi), y(ny,zb_lo:zb_hi)
+    Real(dp), Intent(in) :: xext(zb_lo:zb_hi), aext(zb_lo:zb_hi), zext(zb_lo:zb_hi)
+
+    ! Ouput variables
+    Real(dp), Intent(out) :: ye(zb_lo:zb_hi), cv(zb_lo:zb_hi)
+    Real(dp), Intent(out) :: efermkt(zb_lo:zb_hi), defermktdt9(zb_lo:zb_hi)
+
+    ! Optional variables
+    Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi)
+
+    ! Local variables
+    Integer :: izb
+    Real(dp) :: z53bar
+    Logical, Pointer :: mask(:)
+
+    If ( present(mask_in) ) Then
+      mask(zb_lo:) => mask_in
+    Else
+      mask(zb_lo:) => lzactive(zb_lo:zb_hi)
+    EndIf
+    If ( .not. any(mask) ) Return
+
+    ! Calculate Ye
+    Call y_moment(y,ye,ytot(zb_lo:zb_hi), &
+      & abar(zb_lo:zb_hi),zbar(zb_lo:zb_hi),z2bar(zb_lo:zb_hi),zibar(zb_lo:zb_hi), &
+      & xext(zb_lo:zb_hi),aext(zb_lo:zb_hi),zext(zb_lo:zb_hi),mask_in = mask)
+
+    ! Call the eos
+    Do izb = zb_lo, zb_hi
+      If ( mask(izb) ) Then
+        z53bar = sum(zz53 * y(:,izb)) / ytot(izb)
+        Call eosx(t9(izb),rho(izb),ye(izb),ytot(izb),abar(izb),zbar(izb), &
+          & z53bar,cv(izb),efermkt(izb),defermktdt9(izb))
+      EndIf
+    EndDo
+
+    If ( idiag >= 3 ) Then
+      Do izb = zb_lo, zb_hi
+        If ( mask(izb) ) Then
+          Write(lun_diag,"(a,6es23.15)") 'EOS',t9(izb),rho(izb),ye(izb),cv(izb),efermkt(izb),defermktdt9(izb)
+        EndIf
+      EndDo
+    EndIf
+
+    Return
+  End Subroutine eos_interface_vector
+
+  Subroutine eos_screen_scalar(t9,rho,y,efermkt,defermktdt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext)
     !-----------------------------------------------------------------------------------------------
     ! This routine Bahcall's approach to calculate the factors needed for screening from the input
     ! temperature, density and composition.
@@ -107,10 +231,9 @@ Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztild
     !-----------------------------------------------------------------------------------------------
     Use nuclear_data, Only: ny
     Use xnet_abundances, Only: y_moment
-    Use xnet_constants, Only: avn, bok, clt, e2, ele_en, emass, hbar, pi, pi2, third, two3rd, &
-      & thbim2, twm2bi
-    Use xnet_controls, Only: idiag, iheat, lun_diag
+    Use xnet_controls, Only: idiag, lun_diag
     Use xnet_types, Only: dp
+    Use xnet_util, Only: plasma
     Implicit None
 
     ! Input variables
@@ -127,27 +250,89 @@ Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztild
     Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext)
 
     ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24)
-    Call salpeter_ratio(efermkt,sratio,dsratiodefermkt)
+    Call salpeter_ratio(efermkt,sratio,dztildedt9)
     ztilde = sqrt(z2bar + zbar*sratio)
-    If ( iheat > 0 ) Then
-      dztildedt9 = 0.5*zbar/ztilde * dsratiodefermkt*defermktdt9
-    Else
-      dztildedt9 = 0.0
-    EndIf
+    dztildedt9 = 0.5*zbar/ztilde * dztildedt9*defermktdt9
 
     ! Calculate plasma quantities
-    bkt = bok*t9
-    lambda0 = sqrt(4.0*pi*rho*avn*ytot) * (e2/bkt)**1.5 ! DGC, Eq. 3
-    ae = (3.0 / (4.0*pi*avn*rho*ye))**third ! electron-sphere radius
-    gammae = e2 / (ae*bkt) ! electron Coulomb coupling parameter
-    zinter = zibar / (ztilde**thbim2 * zbar**twm2bi)
-    If ( idiag >= 2 ) Write(lun_diag,"(a14,9es23.15)") 'EOS Screen', &
+    Call plasma(t9,rho,ytot,ye,zbar,zibar,ztilde,zinter,lambda0,gammae)
+    If ( idiag >= 3 ) Write(lun_diag,"(a14,9es23.15)") 'EOS Screen', &
       & t9,rho,ye,z2bar,zbar,sratio,ztilde,ztilde*lambda0,gammae
 
     Return
-  End Subroutine eos_screen
+  End Subroutine eos_screen_scalar
+
+  Subroutine eos_screen_vector(t9,rho,y,efermkt,defermktdt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext,mask_in)
+    !-----------------------------------------------------------------------------------------------
+    ! This routine Bahcall's approach to calculate the factors needed for screening from the input
+    ! temperature, density and composition.
+    !
+    ! Note: the efermkt and defermktdt9 in the subroutine interface are stored in XNet's conditions
+    ! module as etae and detaedt9.  For more general EoS, etae and detaedt9 are their proper selves.
+    !-----------------------------------------------------------------------------------------------
+    Use nuclear_data, Only: ny
+    Use xnet_abundances, Only: y_moment
+    Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive
+    Use xnet_types, Only: dp
+    Use xnet_util, Only: plasma
+    Implicit None
+
+    ! Input variables
+    Real(dp), Intent(in) :: t9(zb_lo:zb_hi), rho(zb_lo:zb_hi), y(ny,zb_lo:zb_hi)
+    Real(dp), Intent(in) :: efermkt(zb_lo:zb_hi), defermktdt9(zb_lo:zb_hi)
+    Real(dp), Intent(in) :: xext(zb_lo:zb_hi), aext(zb_lo:zb_hi), zext(zb_lo:zb_hi)
+
+    ! Output variables
+    Real(dp), Intent(out) :: ztilde(zb_lo:zb_hi), zinter(zb_lo:zb_hi)
+    Real(dp), Intent(out) :: lambda0(zb_lo:zb_hi), gammae(zb_lo:zb_hi)
+    Real(dp), Intent(out) :: dztildedt9(zb_lo:zb_hi)
+
+    ! Optional variables
+    Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi)
+
+    ! Local variables
+    Integer :: izb
+    Logical, Pointer :: mask(:)
+
+    If ( present(mask_in) ) Then
+      mask(zb_lo:) => mask_in
+    Else
+      mask(zb_lo:) => lzactive(zb_lo:zb_hi)
+    EndIf
+    If ( .not. any(mask) ) Return
+
+    ! Calculate Ye
+    Call y_moment(y,ye(zb_lo:zb_hi),ytot(zb_lo:zb_hi), &
+      & abar(zb_lo:zb_hi),zbar(zb_lo:zb_hi),z2bar(zb_lo:zb_hi),zibar(zb_lo:zb_hi), &
+      & xext(zb_lo:zb_hi),aext(zb_lo:zb_hi),zext(zb_lo:zb_hi),mask_in = mask)
+
+    ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24; DGC, Eq. 5)
+    Call salpeter_ratio(efermkt,sratio(zb_lo:zb_hi),dztildedt9,mask_in = mask)
+
+    Do izb = zb_lo, zb_hi
+      If ( mask(izb) ) Then
+        ztilde(izb) = sqrt(z2bar(izb) + sratio(izb)*zbar(izb)) ! DGC, Eq. 4
+        dztildedt9(izb) = 0.5*zbar(izb)/ztilde(izb) * dztildedt9(izb)*defermktdt9(izb)
+
+        ! Calculate plasma quantities
+        Call plasma(t9(izb),rho(izb),ytot(izb),ye(izb),zbar(izb), &
+          & zibar(izb),ztilde(izb),zinter(izb),lambda0(izb),gammae(izb))
+      EndIf
+    EndDo
+    If ( idiag >= 3 ) Then
+      Do izb = zb_lo, zb_hi
+        If ( mask(izb) ) Then
+          Write(lun_diag,"(a14,9es23.15)") 'EOS Screen', &
+            & t9(izb),rho(izb),ye(izb),z2bar(izb),zbar(izb),sratio(izb), &
+            & ztilde(izb),ztilde(izb)*lambda0(izb),gammae(izb)
+        EndIf
+      EndDo
+    EndIf
+
+    Return
+  End Subroutine eos_screen_vector
 
-  Subroutine salpeter_ratio(efmkt,ratio,dratiodefmkt)
+  Subroutine salpeter_ratio_scalar(efmkt,ratio,dratiodefmkt)
     !-----------------------------------------------------------------------------------------------
     ! This routine calculates the Salpeter (1954) ratio f'/f(eta) needed for electron screening, using
     ! a fit to Figure 24 in that paper. efmkt is the ratio of electron chemical potential to kT.
@@ -186,6 +371,44 @@ Subroutine salpeter_ratio(efmkt,ratio,dratiodefmkt)
     EndIf
 
     Return
-  End Subroutine salpeter_ratio
+  End Subroutine salpeter_ratio_scalar
+
+  Subroutine salpeter_ratio_vector(efmkt,ratio,dratiodefmkt,mask_in)
+    !-----------------------------------------------------------------------------------------------
+    ! This routine calculates the Salpeter (1954) ratio f'/f(eta) needed for electron screening, using
+    ! a fit to Figure 24 in that paper. efmkt is the ratio of electron chemical potential to kT.
+    !-----------------------------------------------------------------------------------------------
+    Use xnet_controls, Only: zb_lo, zb_hi, lzactive
+    Use xnet_types, Only: dp
+    Implicit None
+
+    ! Input variables
+    Real(dp), Intent(in) :: efmkt(zb_lo:zb_hi)
+
+    ! Output variables
+    Real(dp), Intent(out) :: ratio(zb_lo:zb_hi), dratiodefmkt(zb_lo:zb_hi)
+
+    ! Optional variables
+    Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi)
+
+    ! Local variables
+    Integer :: izb
+    Logical, Pointer :: mask(:)
+
+    If ( present(mask_in) ) Then
+      mask(zb_lo:) => mask_in
+    Else
+      mask(zb_lo:) => lzactive(zb_lo:zb_hi)
+    EndIf
+    If ( .not. any(mask) ) Return
+
+    Do izb = zb_lo, zb_hi
+      If ( mask(izb) ) Then
+        Call salpeter_ratio_scalar(efmkt(izb),ratio(izb),dratiodefmkt(izb))
+      EndIf
+    EndDo
+
+    Return
+  End Subroutine salpeter_ratio_vector
 
 End Module xnet_eos
diff --git a/source/xnet_eos_helm.f90 b/source/xnet_eos_helm.f90
index ebbbf408..aff85f84 100644
--- a/source/xnet_eos_helm.f90
+++ b/source/xnet_eos_helm.f90
@@ -6,8 +6,25 @@
 !***************************************************************************************************
 
 Module xnet_eos
+  Use xnet_types, Only: dp
   Implicit None
   include 'vector_eos.dek'
+  Real(dp), Allocatable :: ye(:), ytot(:), abar(:), zbar(:), z2bar(:), zibar(:), sratio(:)
+
+  Interface eos_interface
+    Module Procedure eos_interface_scalar
+    Module Procedure eos_interface_vector
+  End Interface
+
+  Interface eos_screen
+    Module Procedure eos_screen_scalar
+    Module Procedure eos_screen_vector
+  End Interface
+
+  Interface salpeter_ratio
+    Module Procedure salpeter_ratio_scalar
+    Module Procedure salpeter_ratio_vector
+  End Interface
 
 Contains
 
@@ -15,32 +32,47 @@ Subroutine eos_initialize
     !-----------------------------------------------------------------------------------------------
     ! This routine updates the Equation of State for changes in temperature and density.
     !-----------------------------------------------------------------------------------------------
+    Use xnet_controls, Only: nzevolve
     Implicit None
+
     Call read_helm_table
 
+    Allocate (ye(nzevolve))
+    Allocate (ytot(nzevolve))
+    Allocate (abar(nzevolve))
+    Allocate (zbar(nzevolve))
+    Allocate (z2bar(nzevolve))
+    Allocate (zibar(nzevolve))
+    Allocate (sratio(nzevolve))
+    ye = 0.0
+    ytot = 0.0
+    abar = 0.0
+    zbar = 0.0
+    z2bar = 0.0
+    zibar = 0.0
+    sratio = 0.0
+
     Return
   End Subroutine eos_initialize
 
-  Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext)
-    Use nuclear_data, Only: ny
-    Use xnet_abundances, Only: y_moment
-    Use xnet_constants, Only: avn, epmev
-    Use xnet_controls, Only: idiag, iheat, iscrn, lun_diag
+  Subroutine eosx(t9,rho,ye,abar,zbar,cv,etae,detaedt9)
+    !-----------------------------------------------------------------------------------------------
+    ! This routine interfaces with and calls the underlying EoS.
+    !-----------------------------------------------------------------------------------------------
+    Use xnet_constants, Only: amu
+    Use xnet_controls, Only: iheat, iscrn
     Use xnet_types, Only: dp
     Implicit None
 
     ! Input variables
-    Real(dp), Intent(in) :: t9, rho, y(ny), xext, aext, zext
+    Real(dp), Intent(in) :: t9, rho, ye, abar, zbar
 
     ! Ouput variables
-    Real(dp), Intent(out) :: ye, cv, etae, detaedt9
-
-    ! Local variables
-    Real(dp) :: ytot, abar, zbar, z2bar, zibar
-
-    ! Calculate Ye
-    Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext)
+    Real(dp), Intent(out) :: cv, etae, detaedt9
 
+    cv = 0.0
+    etae = 0.0
+    detaedt9 = 0.0
     If ( iscrn > 0 .or. iheat > 0 ) Then
 
       ! Load input variables for the eos
@@ -55,30 +87,108 @@ Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext)
       Call helmeos
 
       ! Convert units from ergs/g to MeV/nucleon and K to GK
+      cv = cv_row(1) * amu * 1e9
       etae = etaele_row(1)
-      detaedt9 = detat_row(1)*1e9
-      cv = cv_row(1)*1.0d9/epmev/avn
-    Else
-      etae = 0.0
-      detaedt9 = 0.0
-      cv = 0.0
+      detaedt9 = detat_row(1) * 1e9
     EndIf
+
+  End Subroutine eosx
+
+  Subroutine eos_interface_scalar(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext)
+    !-----------------------------------------------------------------------------------------------
+    ! This routine updates the equation of state for changes in temperature and density.
+    !-----------------------------------------------------------------------------------------------
+    Use nuclear_data, Only: ny
+    Use xnet_abundances, Only: y_moment
+    Use xnet_controls, Only: idiag, lun_diag
+    Use xnet_types, Only: dp
+    Implicit None
+
+    ! Input variables
+    Real(dp), Intent(in) :: t9, rho, y(ny), xext, aext, zext
+
+    ! Ouput variables
+    Real(dp), Intent(out) :: ye, cv, etae, detaedt9
+
+    ! Local variables
+    Real(dp) :: ytot, abar, zbar, z2bar, zibar
+
+    ! Calculate Ye
+    Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext)
+
+    ! Call the eos
+    Call eosx(t9,rho,ye,abar,zbar,cv,etae,detaedt9)
+
     If ( idiag >= 3 ) Write(lun_diag,"(a,6es23.15)") 'EOS',t9,rho,ye,cv,etae,detaedt9
 
     Return
-  End Subroutine eos_interface
+  End Subroutine eos_interface_scalar
+
+  Subroutine eos_interface_vector(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext,mask_in)
+    !-----------------------------------------------------------------------------------------------
+    ! This routine updates the equation of state for changes in temperature and density.
+    !-----------------------------------------------------------------------------------------------
+    Use nuclear_data, Only: ny
+    Use xnet_abundances, Only: y_moment
+    Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive
+    Use xnet_types, Only: dp
+    Implicit None
+
+    ! Input variables
+    Real(dp), Intent(in) :: t9(zb_lo:zb_hi), rho(zb_lo:zb_hi), y(ny,zb_lo:zb_hi)
+    Real(dp), Intent(in) :: xext(zb_lo:zb_hi), aext(zb_lo:zb_hi), zext(zb_lo:zb_hi)
+
+    ! Ouput variables
+    Real(dp), Intent(out) :: ye(zb_lo:zb_hi), cv(zb_lo:zb_hi)
+    Real(dp), Intent(out) :: etae(zb_lo:zb_hi), detaedt9(zb_lo:zb_hi)
+
+    ! Optional variables
+    Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi)
+
+    ! Local variables
+    Integer :: izb
+    Logical, Pointer :: mask(:)
+
+    If ( present(mask_in) ) Then
+      mask(zb_lo:) => mask_in
+    Else
+      mask(zb_lo:) => lzactive(zb_lo:zb_hi)
+    EndIf
+    If ( .not. any(mask) ) Return
+
+    ! Calculate Ye
+    Call y_moment(y,ye,ytot(zb_lo:zb_hi), &
+      & abar(zb_lo:zb_hi),zbar(zb_lo:zb_hi),z2bar(zb_lo:zb_hi),zibar(zb_lo:zb_hi), &
+      & xext(zb_lo:zb_hi),aext(zb_lo:zb_hi),zext(zb_lo:zb_hi),mask_in = mask)
+
+    ! Call the eos
+    Do izb = zb_lo, zb_hi
+      If ( mask(izb) ) Then
+        Call eosx(t9(izb),rho(izb),ye(izb),abar(izb),zbar(izb),cv(izb),etae(izb),detaedt9(izb))
+      EndIf
+    EndDo
+
+    If ( idiag >= 3 ) Then
+      Do izb = zb_lo, zb_hi
+        If ( mask(izb) ) Then
+          Write(lun_diag,"(a,6es23.15)") 'EOS',t9(izb),rho(izb),ye(izb),cv(izb),etae(izb),detaedt9(izb)
+        EndIf
+      EndDo
+    EndIf
+
+    Return
+  End Subroutine eos_interface_vector
 
-  Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext)
+  Subroutine eos_screen_scalar(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext)
     !-----------------------------------------------------------------------------------------------
     ! This routine uses the current composition and prior updates to the Equation of State to
     ! calculate the factors needed for screening.
     !-----------------------------------------------------------------------------------------------
     Use nuclear_data, Only: ny
     Use xnet_abundances, Only: y_moment
-    Use xnet_constants, Only: avn, bok, clt, e2, ele_en, emass, hbar, pi, pi2, third, two3rd, &
-      & thbim2, twm2bi
     Use xnet_controls, Only: idiag, iheat, lun_diag
     Use xnet_types, Only: dp
+    Use xnet_util, Only: plasma
     Implicit None
 
     ! Input variables
@@ -88,37 +198,97 @@ Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztild
     Real(dp), Intent(out) :: ztilde, zinter, lambda0, gammae, dztildedt9
 
     ! Local variables
-    Real(dp) :: ye, ytot, bkt, abar, zbar, z2bar, zibar
-    Real(dp) :: sratio, ae, dsratiodeta
+    Real(dp) :: ye, ytot, abar, zbar, z2bar, zibar
+    Real(dp) :: sratio
 
     ! Calculate Ye and other needed moments of the abundance distribution
     Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext)
 
-    ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24)
-    Call salpeter_ratio(etae,sratio,dsratiodeta)
-    ztilde = sqrt(z2bar + zbar*sratio)
-    If ( iheat > 0 ) Then
-      dztildedt9 = 0.5*zbar/ztilde * dsratiodeta*detaedt9
-    Else
-      dztildedt9 = 0.0
-    EndIf
+    ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24; DGC, Eq. 5)
+    Call salpeter_ratio(etae,sratio,dztildedt9)
+    ztilde = sqrt(z2bar + zbar*sratio) ! DGC, Eq. 4
+    dztildedt9 = 0.5*zbar/ztilde * dztildedt9*detaedt9
 
     ! Calculate plasma quantities
-    bkt = bok*t9
-    lambda0 = sqrt(4.0*pi*rho*avn*ytot) * (e2/bkt)**1.5 ! DGC, Eq. 3
-    ae = (3.0 / (4.0*pi*avn*rho*ye))**third ! electron-sphere radius
-    gammae = e2 / (ae*bkt) ! electron Coulomb coupling parameter
-    zinter = zibar / (ztilde**thbim2 * zbar**twm2bi)
+    Call plasma(t9,rho,ytot,ye,zbar,zibar,ztilde,zinter,lambda0,gammae)
     If ( idiag >= 3 ) Write(lun_diag,"(a14,9es23.15)") 'EOS Screen', &
       & t9,rho,ye,z2bar,zbar,sratio,ztilde,ztilde*lambda0,gammae
 
     Return
-  End Subroutine eos_screen
+  End Subroutine eos_screen_scalar
+
+  Subroutine eos_screen_vector(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext,mask_in)
+    !-----------------------------------------------------------------------------------------------
+    ! This routine uses the current composition and prior updates to the Equation of State to
+    ! calculate the factors needed for screening.
+    !-----------------------------------------------------------------------------------------------
+    Use nuclear_data, Only: ny
+    Use xnet_abundances, Only: y_moment
+    Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive
+    Use xnet_types, Only: dp
+    Use xnet_util, Only: plasma
+    Implicit None
+
+    ! Input variables
+    Real(dp), Intent(in) :: t9(zb_lo:zb_hi), rho(zb_lo:zb_hi), y(ny,zb_lo:zb_hi)
+    Real(dp), Intent(in) :: etae(zb_lo:zb_hi), detaedt9(zb_lo:zb_hi)
+    Real(dp), Intent(in) :: xext(zb_lo:zb_hi), aext(zb_lo:zb_hi), zext(zb_lo:zb_hi)
+
+    ! Output variables
+    Real(dp), Intent(out) :: ztilde(zb_lo:zb_hi), zinter(zb_lo:zb_hi)
+    Real(dp), Intent(out) :: lambda0(zb_lo:zb_hi), gammae(zb_lo:zb_hi)
+    Real(dp), Intent(out) :: dztildedt9(zb_lo:zb_hi)
+
+    ! Optional variables
+    Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi)
+
+    ! Local variables
+    Integer :: izb
+    Logical, Pointer :: mask(:)
+
+    If ( present(mask_in) ) Then
+      mask(zb_lo:) => mask_in
+    Else
+      mask(zb_lo:) => lzactive(zb_lo:zb_hi)
+    EndIf
+    If ( .not. any(mask) ) Return
+
+    ! Calculate Ye
+    Call y_moment(y,ye(zb_lo:zb_hi),ytot(zb_lo:zb_hi), &
+      & abar(zb_lo:zb_hi),zbar(zb_lo:zb_hi),z2bar(zb_lo:zb_hi),zibar(zb_lo:zb_hi), &
+      & xext(zb_lo:zb_hi),aext(zb_lo:zb_hi),zext(zb_lo:zb_hi),mask_in = mask)
+
+    ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24; DGC, Eq. 5)
+    Call salpeter_ratio(etae,sratio(zb_lo:zb_hi),dztildedt9,mask_in = mask)
+
+    Do izb = zb_lo, zb_hi
+      If ( mask(izb) ) Then
+        ztilde(izb) = sqrt(z2bar(izb) + sratio(izb)*zbar(izb)) ! DGC, Eq. 4
+        dztildedt9(izb) = 0.5*zbar(izb)/ztilde(izb) * dztildedt9(izb)*detaedt9(izb)
+
+        ! Calculate plasma quantities
+        Call plasma(t9(izb),rho(izb),ytot(izb),ye(izb),zbar(izb), &
+          & zibar(izb),ztilde(izb),zinter(izb),lambda0(izb),gammae(izb))
+      EndIf
+    EndDo
+    If ( idiag >= 3 ) Then
+      Do izb = zb_lo, zb_hi
+        If ( mask(izb) ) Then
+          Write(lun_diag,"(a14,9es23.15)") 'EOS Screen', &
+            & t9(izb),rho(izb),ye(izb),z2bar(izb),zbar(izb),sratio(izb), &
+            & ztilde(izb),ztilde(izb)*lambda0(izb),gammae(izb)
+        EndIf
+      EndDo
+    EndIf
+
+    Return
+  End Subroutine eos_screen_vector
 
-  Subroutine salpeter_ratio(eta,ratio,dratiodeta)
+  Subroutine salpeter_ratio_scalar(eta,ratio,dratiodeta)
     !-----------------------------------------------------------------------------------------------
     ! This routine calculates the Salpeter (1954) ratio f'/f(eta) needed for electron screening.
     ! eta is the ratio of electron chemical potential to kT.
+    ! f'/f is also defined as theta_e in DeWitt+ (1973).
     !
     ! Calculation uses Fermi function relation d/dx f_(k+1) = (k+1) f_k and the rational function
     ! expansions of Fukushima (2015; AMC 259 708) for the F-D integrals of order 1/2, -1/2, and -3/2.
@@ -142,7 +312,7 @@ Subroutine salpeter_ratio(eta,ratio,dratiodeta)
     fermim = fdm1h(eta)
     fermip = fd1h(eta)
 
-    ! Evalutate the Salpeter ratio
+    ! Evalutate the Salpeter ratio (extra factor of 1/2 from FD integral definition)
     ratio = 0.5 * fermim/fermip
     If ( iheat > 0 ) Then
       dfmdeta = -0.5 * fdm3h(eta)
@@ -153,7 +323,48 @@ Subroutine salpeter_ratio(eta,ratio,dratiodeta)
     EndIf
 
     Return
-  End Subroutine salpeter_ratio
+  End Subroutine salpeter_ratio_scalar
 
-End  Module xnet_eos
+  Subroutine salpeter_ratio_vector(eta,ratio,dratiodeta,mask_in)
+    !-----------------------------------------------------------------------------------------------
+    ! This routine calculates the Salpeter (1954) ratio f'/f(eta) needed for electron screening.
+    ! eta is the ratio of electron chemical potential to kT.
+    ! f'/f is also defined as theta_e in DeWitt+ (1973).
+    !
+    ! Calculation uses Fermi function relation d/dx f_(k+1) = (k+1) f_k and the rational function
+    ! expansions of Fukushima (2015; AMC 259 708) for the F-D integrals of order 1/2, -1/2, and -3/2.
+    !-----------------------------------------------------------------------------------------------
+    Use xnet_controls, Only: zb_lo, zb_hi, lzactive
+    Use xnet_types, Only: dp
+    Implicit None
+
+    ! Input variables
+    Real(dp), Intent(in) :: eta(zb_lo:zb_hi)
+
+    ! Output variables
+    Real(dp), Intent(out) :: ratio(zb_lo:zb_hi), dratiodeta(zb_lo:zb_hi)
+
+    ! Optional variables
+    Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi)
+
+    ! Local variables
+    Integer :: izb
+    Logical, Pointer :: mask(:)
+
+    If ( present(mask_in) ) Then
+      mask(zb_lo:) => mask_in
+    Else
+      mask(zb_lo:) => lzactive(zb_lo:zb_hi)
+    EndIf
+    If ( .not. any(mask) ) Return
+
+    Do izb = zb_lo, zb_hi
+      If ( mask(izb) ) Then
+        Call salpeter_ratio_scalar(eta(izb),ratio(izb),dratiodeta(izb))
+      EndIf
+    EndDo
+
+    Return
+  End Subroutine salpeter_ratio_vector
 
+End Module xnet_eos
diff --git a/source/xnet_eos_starkiller.f90 b/source/xnet_eos_starkiller.f90
index 0c223d91..321825b9 100644
--- a/source/xnet_eos_starkiller.f90
+++ b/source/xnet_eos_starkiller.f90
@@ -213,9 +213,7 @@ Subroutine eos_screen_scalar(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae
     ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24; DGC, Eq. 5)
     Call salpeter_ratio(etae,sratio,dztildedt9)
     ztilde = sqrt(z2bar + zbar*sratio) ! DGC, Eq. 4
-    If ( iheat > 0 ) Then
-      dztildedt9 = 0.5*zbar/ztilde * dztildedt9*detaedt9
-    EndIf
+    dztildedt9 = 0.5*zbar/ztilde * dztildedt9*detaedt9
 
     ! Calculate plasma quantities
     Call plasma(t9,rho,ytot,ye,zbar,zibar,ztilde,zinter,lambda0,gammae)
@@ -232,7 +230,7 @@ Subroutine eos_screen_vector(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae
     !-----------------------------------------------------------------------------------------------
     Use nuclear_data, Only: ny
     Use xnet_abundances, Only: y_moment
-    Use xnet_controls, Only: idiag, iheat, lun_diag, zb_lo, zb_hi, lzactive
+    Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive
     Use xnet_types, Only: dp
     Use xnet_util, Only: plasma
     Implicit None
@@ -272,9 +270,7 @@ Subroutine eos_screen_vector(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae
     Do izb = zb_lo, zb_hi
       If ( mask(izb) ) Then
         ztilde(izb) = sqrt(z2bar(izb) + sratio(izb)*zbar(izb)) ! DGC, Eq. 4
-        If ( iheat > 0 ) Then
-          dztildedt9(izb) = 0.5*zbar(izb)/ztilde(izb) * dztildedt9(izb)*detaedt9(izb)
-        EndIf
+        dztildedt9(izb) = 0.5*zbar(izb)/ztilde(izb) * dztildedt9(izb)*detaedt9(izb)
 
         ! Calculate plasma quantities
         Call plasma(t9(izb),rho(izb),ytot(izb),ye(izb),zbar(izb), &