diff --git a/Build/Scripts/HYPRE/build_hypre.bat b/Build/Scripts/HYPRE/build_hypre.bat
index 434543107ab..c4fd6442263 100644
--- a/Build/Scripts/HYPRE/build_hypre.bat
+++ b/Build/Scripts/HYPRE/build_hypre.bat
@@ -1,5 +1,5 @@
 @echo off
-set LIB_TAG=v2.33.0
+set LIB_TAG=v3.0.0
 
 ::*** library and tag name are the same
 
diff --git a/Build/Scripts/HYPRE/build_hypre.sh b/Build/Scripts/HYPRE/build_hypre.sh
index 767a3bf6e20..b6dfc5ec91a 100755
--- a/Build/Scripts/HYPRE/build_hypre.sh
+++ b/Build/Scripts/HYPRE/build_hypre.sh
@@ -1,5 +1,5 @@
 #!/bin/bash
-HYPRE_LIB_TAG=v2.33.0
+HYPRE_LIB_TAG=v3.0.0
 
 CONFMAKE=$1
 CLEAN_HYPRE=$2
@@ -14,10 +14,10 @@ fi
 
 echo "Checking for hypre library..."
 
-if [ -d "$FIREMODELS/libs/hypre" ]; then
+if [ -d "$FIREMODELS/libs/hypre/$HYPRE_LIB_TAG" ]; then
   echo "Hypre library exists.  Skipping hypre build."
   # List all directories under $FIREMODELS/libs/hypre
-  hypre_lib_dir=$(ls -d $FIREMODELS/libs/hypre/*/)
+  hypre_lib_dir=$(ls -d $FIREMODELS/libs/hypre/$HYPRE_LIB_TAG/)
   # Extract the version part (removes the leading path)
   HYPRE_VERSION=$(basename $hypre_lib_dir)
   export HYPRE_HOME=$FIREMODELS/libs/hypre/$HYPRE_VERSION
diff --git a/Build/Scripts/SUNDIALS/build_sundials.bat b/Build/Scripts/SUNDIALS/build_sundials.bat
index 0cb3c630320..2db4f153877 100644
--- a/Build/Scripts/SUNDIALS/build_sundials.bat
+++ b/Build/Scripts/SUNDIALS/build_sundials.bat
@@ -1,5 +1,5 @@
 @echo off
-set LIB_TAG=v6.7.0
+set LIB_TAG=v7.5.0
 
 ::*** library and tag name are the same
 
@@ -123,8 +123,8 @@ cmake ..\  ^
 -DBUILD_SHARED_LIBS=OFF ^
 -DCMAKE_INSTALL_LIBDIR="lib" ^
 -DCMAKE_MAKE_PROGRAM="%CMAKE_MAKE_PROGRAM%" ^
--DCMAKE_C_FLAGS_RELEASE="${CMAKE_C_FLAGS_RELEASE} /MT" ^
--DCMAKE_C_FLAGS_DEBUG="${CMAKE_C_FLAGS_DEBUG} /MTd"
+-DCMAKE_MSVC_RUNTIME_LIBRARY="MultiThreaded" ^
+-DSUNDIALS_LOGGING_LEVEL=0
 
 ::*** build and install sundials
 
diff --git a/Build/Scripts/SUNDIALS/build_sundials.sh b/Build/Scripts/SUNDIALS/build_sundials.sh
index 0b0478e186a..12808d212a6 100755
--- a/Build/Scripts/SUNDIALS/build_sundials.sh
+++ b/Build/Scripts/SUNDIALS/build_sundials.sh
@@ -1,5 +1,5 @@
 #!/bin/bash
-SUNDIALS_LIB_TAG=v6.7.0
+SUNDIALS_LIB_TAG=v7.5.0
 
 CONFMAKE=$1
 CLEAN_SUNDIALS=$2
@@ -15,10 +15,10 @@ fi
 
 echo "Checking for sundials library..."
 
-if [ -d "$FIREMODELS/libs/sundials" ]; then
+if [ -d "$FIREMODELS/libs/sundials/$SUNDIALS_LIB_TAG" ]; then
   echo "Sundials library exists.  Skipping sundials build."
   # List all directories under $FIREMODELS/libs/sundials
-  sundials_lib_dir=$(ls -d $FIREMODELS/libs/sundials/*/)
+  sundials_lib_dir=$(ls -d $FIREMODELS/libs/sundials/$SUNDIALS_LIB_TAG/)
   # Extract the version part (removes the leading path)
   SUNDIALS_VERSION=$(basename $sundials_lib_dir)
   export SUNDIALS_HOME=$FIREMODELS/libs/sundials/$SUNDIALS_VERSION
diff --git a/Build/Scripts/SUNDIALS/confmake.sh b/Build/Scripts/SUNDIALS/confmake.sh
index b6f0caf98c0..ff35dc023f6 100755
--- a/Build/Scripts/SUNDIALS/confmake.sh
+++ b/Build/Scripts/SUNDIALS/confmake.sh
@@ -12,6 +12,7 @@ cmake_args=(
   -DEXAMPLES_ENABLE_F2003=OFF
   -DENABLE_OPENMP=OFF
   -DCMAKE_INSTALL_LIBDIR="lib"
+  -DSUNDIALS_LOGGING_LEVEL=0
 )
 
 # Add OSX deployment target if building for macOS
diff --git a/Build/Scripts/set_compilers.sh b/Build/Scripts/set_compilers.sh
index 11bc912ae40..ce295e30527 100755
--- a/Build/Scripts/set_compilers.sh
+++ b/Build/Scripts/set_compilers.sh
@@ -72,9 +72,9 @@ else  # Default to GNU compilers
     select_compiler_from_system COMP_FC mpifort
 fi
 
-echo "Thirdparty libs C Compiler COMP_CC=$COMP_CC"
-echo "Thirdparty libs C++ compiler COMP_CXX=$COMP_CXX"
-echo "Firemodels and Thirdparty libs Fortran compiler COMP_FC=$COMP_FC"
+echo "Third party libs C Compiler COMP_CC=$COMP_CC"
+echo "Third party libs C++ compiler COMP_CXX=$COMP_CXX"
+echo "Firemodels and third party libs Fortran compiler COMP_FC=$COMP_FC"
 
 export COMP_CC=$COMP_CC
 export COMP_CXX=$COMP_CXX
diff --git a/Build/makefile b/Build/makefile
index 9cc9b95db3e..ac28083a584 100644
--- a/Build/makefile
+++ b/Build/makefile
@@ -11,9 +11,8 @@
 
 VPATH = ../Source
 
-# 3rd part library versions
-HYPRE_VERSION=v2.33.0
-SUNDIALS_VERSION=v6.7.0
+HYPRE_VERSION=v3.0.0
+SUNDIALS_VERSION=v7.5.0
 HDF5_VERSION=v1.14.5
 
 ifeq ($(shell echo "check_quotes"),"check_quotes")
@@ -115,8 +114,8 @@ endif # MKLROOT test
 ifdef SUNDIALS_HOME # This assumes the SUNDIALS library is available.
    FFLAGS_SUNDIALS = -DWITH_SUNDIALS ${SUNDIALS_INFO} -I"${SUNDIALS_HOME}/fortran"
    FFLAGS_SUNDIALS_WIN = -DWITH_SUNDIALS ${SUNDIALS_INFO} /I"${SUNDIALS_HOME}\fortran"
-   LFLAGS_SUNDIALS = ${SUNDIALS_HOME}/lib/libsundials_fcvode_mod.a ${SUNDIALS_HOME}/lib/libsundials_fnvecserial_mod.a ${SUNDIALS_HOME}/lib/libsundials_cvode.a
-   LFLAGS_SUNDIALS_WIN = ${SUNDIALS_HOME}\lib\sundials_fcvode_mod.lib ${SUNDIALS_HOME}\lib\sundials_fnvecserial_mod.lib ${SUNDIALS_HOME}\lib\sundials_cvode.lib
+   LFLAGS_SUNDIALS = ${SUNDIALS_HOME}/lib/libsundials_fcvode_mod.a ${SUNDIALS_HOME}/lib/libsundials_fnvecserial_mod.a ${SUNDIALS_HOME}/lib/libsundials_cvode.a ${SUNDIALS_HOME}/lib/libsundials_fcore_mod.a ${SUNDIALS_HOME}/lib/libsundials_core.a
+   LFLAGS_SUNDIALS_WIN = ${SUNDIALS_HOME}\lib\sundials_fcvode_mod_static.lib ${SUNDIALS_HOME}\lib\sundials_fnvecserial_mod_static.lib ${SUNDIALS_HOME}\lib\sundials_cvode_static.lib ${SUNDIALS_HOME}\lib\sundials_fcore_mod.lib ${SUNDIALS_HOME}\lib\sundials_core_static.lib /link  /NODEFAULTLIB:MSVCRTD /NODEFAULTLIB:libcmtd.lib
 endif
 
 ifdef HYPRE_HOME # This assumes the HYPRE library is available.
diff --git a/CMakeLists.txt b/CMakeLists.txt
index bae03c25a7e..63050ba1577 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -157,7 +157,7 @@ if(USE_HYPRE OR USE_HYPRE_NVIDIA OR USE_HYPRE_AMDGPU OR USE_HYPRE_INTELGPU)
         include(FetchContent)
         # As we are not using the system hypre, we need to choose the version we
         # want
-        set(HYPRE_GIT_VERSION "v2.32.0-24-63331f19c" )
+        set(HYPRE_GIT_VERSION "3.0.0" )
         FetchContent_Declare(
             HYPRE
             GIT_REPOSITORY https://github.com/hypre-space/hypre.git
@@ -216,7 +216,7 @@ if(USE_SUNDIALS)
         include(FetchContent)
         # As we are not using the system sundials, we need to choose the version
         # we want
-        set(SUNDIALS_GIT_VERSION "6.7.0")
+        set(SUNDIALS_GIT_VERSION "7.5.0")
         FetchContent_Declare(
             SUNDIALS
             GIT_REPOSITORY https://github.com/LLNL/sundials.git
diff --git a/Manuals/FDS_Technical_Reference_Guide/Equation_Chapter.tex b/Manuals/FDS_Technical_Reference_Guide/Equation_Chapter.tex
index a5e799c7687..651a5d3b630 100644
--- a/Manuals/FDS_Technical_Reference_Guide/Equation_Chapter.tex
+++ b/Manuals/FDS_Technical_Reference_Guide/Equation_Chapter.tex
@@ -149,9 +149,6 @@ \section{Mass and Species Transport}
 \be \dod{\rho}{t} + \nabla\!\cdot (\rho \bu)  =  \dm_{\rm b}'''  \label{mass} \ee
 because $\sum Z_\alpha=1$ and $\sum \dm_\alpha''' = 0$ and $\sum \dm_{\rm b,\alpha}'''=\dm_{\rm b}'''$, by definition, and because it is assumed that $\sum \rho D_\alpha \nabla Z_\alpha = 0$. This last assertion is not true, in general.  The diffusive flux for the most abundant local species is corrected to enforce the constraint.
 
-
-
-
 \paragraph{Enforcing Realizability} ~\\
 
 \noindent Realizability of species mass fractions requires $Y_\alpha \ge 0$ for all $\alpha$ and $\sum Y_\alpha = 1$.  Note that this is more restrictive than the boundedness constraint, which simply requires $0 \le Y_\alpha \le 1$.
@@ -160,6 +157,10 @@ \section{Mass and Species Transport}
 
 With this approach we must take care to ensure $\sum \rho D_\alpha \nabla Z_\alpha = 0$.  Our strategy is to absorb any errors in diffusive transport into the most abundant species \emph{locally}.  That is, for a given cell face we set $\rho D_m \nabla Z_m = -\sum_{\alpha\ne m} \rho D_\alpha \nabla Z_\alpha$, where $m$ is the most abundant species adjacent to that face.  Note that since FDS is typically used as an LES code mass transport by molecular diffusion may be two or three orders of magnitude less than turbulent transport, which uses the same turbulent diffusion coefficient for all species. Therefore, the errors in summation of the diffusive fluxes tend to be small.
 
+\paragraph{Enforcing the Equation of State} ~\\
+
+\noindent Given the strategy of solving $N_s$ transport equations and obtaining $\rho = \sum_{\alpha=1}^{N_s} (\rho Y)_\alpha$ to enforce realizability, we must take care to ensure that we maintain the equation of state at cell faces after flux limited interpolation.  In Sec.~\ref{sec_flux_limiters}, we describe a correction scheme that enforces the equation of state at cell faces and thus maintains isothermal flows for multi-component mixtures with variable molecular weights.
+
 \section{Low Mach Number Approximation}
 
 For low speed applications like fire, Rehm and Baum~\cite{Rehm:1} observed that the spatially and temporally resolved pressure, $p$, can be decomposed into a ``background'' pressure, $\bp(z,t)$, plus a perturbation, $\tp(x,y,z,t)$, with only the background pressure retained in the equation of state (ideal gas law):
diff --git a/Manuals/FDS_Technical_Reference_Guide/Mass_Chapter.tex b/Manuals/FDS_Technical_Reference_Guide/Mass_Chapter.tex
index 88092426b92..4924ef55476 100644
--- a/Manuals/FDS_Technical_Reference_Guide/Mass_Chapter.tex
+++ b/Manuals/FDS_Technical_Reference_Guide/Mass_Chapter.tex
@@ -78,7 +78,7 @@ \subsection{Flux Limiters}
 
 For uniform flow velocity, a fundamental property of the exact solution to the equations governing scalar transport is that the total variation of the scalar field (the sum of the absolute values of the scalar differences between neighboring cells) is either preserved or diminished (never increased).  In other words, no new extrema are created.  Numerical schemes which preserve this property are referred to as total variation diminishing (TVD) schemes.  The practical importance of using a TVD scheme for fire modeling is that such a scheme is able to accurately track coherent vortex structure in turbulent flames and does not develop spurious reaction zones.
 
-FDS employs two second-order TVD schemes as options for scalar transport: Superbee and CHARM.  Superbee \cite{Roe:1986} is recommended for LES because it more accurately preserves the scalar variance for coarse grid solutions that are not expected to be smooth.  Due to the gradient steepening applied in Superbee, however, the convergence degrades at small grid spacing for smooth solutions (the method will revert to a stair-step pattern instead of the exact solution).  CHARM \cite{Zhou:1995}, though slightly more dissipative than Superbee, is convergent, and is therefore the better choice for DNS calculations where the flame front is well resolved.
+FDS employs two second-order TVD schemes as options for scalar transport: Superbee and CHARM.  Superbee \cite{Roe:1986} is recommended for VLES because it more accurately preserves the scalar variance for coarse grid solutions that are not expected to be smooth.  Due to the gradient steepening applied in Superbee, however, the convergence degrades at small grid spacing for smooth solutions (the method will revert to a stair-step pattern instead of the exact solution).  CHARM \cite{Zhou:1995}, though slightly more dissipative than Superbee, is convergent, and is therefore the better choice for LES and DNS calculations.
 
 To illustrate how flux limiters are applied to the scalar transport equations, below we discretize the advection terms in Eq.~(\ref{species}) in one dimension:
 \be  \frac{(\rho Z)_{i}^* - (\rho Z)_{i}^n}{\dt}
@@ -94,12 +94,12 @@ \subsection{Flux Limiters}
 \begin{table}[H]
 \begin{center}
 \begin{tabular}{lc}
-Flux Limiter                           & $B(r)$                         \\
+Flux Limiter                                    & $B(r)$                         \\
 \hline
-Central Difference                     & 1                              \\
-Godunov                                & 0                              \\
-Superbee \cite{Roe:1986} (LES default) & $\max(0,\min(2r,1),\min(r,2))$ \\
-CHARM \cite{Zhou:1995} (DNS default)   & $s(3s+1)/(s+1)^2$; $s=1/r$     \\
+Central Difference                              & 1                              \\
+Godunov                                         & 0                              \\
+Superbee \cite{Roe:1986} (VLES default)         & $\max(0,\min(2r,1),\min(r,2))$ \\
+CHARM \cite{Zhou:1995} (LES and DNS default)    & $s(3s+1)/(s+1)^2$; $s=1/r$     \\
 \end{tabular}
 \end{center}
 \end{table}
@@ -123,21 +123,38 @@ \subsubsection{Notes on Implementation}
 
 The Central Difference and Godunov limiters are included for completeness, debugging, and educational purposes.  These schemes have little utility for typical FDS applications.
 
-% \subsubsection{Dealing with Variable Molecular Weights}
+\subsubsection{Dealing with Variable Molecular Weights}
 
-%% leave this here for a moment as a reminder to write up the constant limiter coefficient method we are now using.
+Handling multi-component mixtures with variable molecular weights requires special care when constructing fluxes, otherwise we may violate the equation of state.  The problem is easiest to illustrate by considering an isothermal (not necessarily constant density) flow.  Maintaining isothermal flow requires
+\begin{equation}
+T = \frac{\overline{W} \bar{p}}{\rho R} = \frac{\bar{p}}{R \rho \sum_\alpha \frac{Z_\alpha}{W_\alpha}} = \frac{\bar{p}}{R \sum_\alpha \frac{\rho Z_\alpha}{W_\alpha}}
+\end{equation}
+to be constant and uniform at all cells and faces. Therefore, with $\bar{p}$ and $R$ constant and uniform, we must maintain
+\begin{equation}
+\label{eq:rho_mw}
+\sum_\alpha \frac{(\rho Z_\alpha)}{W_\alpha} = \frac{\rho}{\overline{W}}
+\end{equation}
+The above condition is automatically satisfied in the cases of using Godunov or Central differencing or in the case of binary flow (two species).  However, if we apply a second-order flux limiter, such as Superbee or CHARM, independently to each species in a multi-component (three or more species) flow with variable molecular weights, then this condition is easily violated at a cell face.
+
+For example, consider the case of a multi-component, constant density flow.  A flux limited face value for species $\alpha$ is a sum of the values in adjacent cells times a coefficient $c_{i,\alpha}$ for each cell $i$.  At a face we require have $\sum_\alpha \sum_i c_{i,\alpha} (\rho Z_\alpha)_i = \rho = \mbox{constant}$.  The only way to guarantee $\rho$ is constant is if the $c_{i,\alpha}$ are the same for each $\alpha$, that is, $c_{i,\alpha} = c_i$.  Rearranging then gives $\sum_i c_i \sum_\alpha (\rho Z_\alpha)_i = \sum_i c_i \rho_i = \rho$, which just says that we need to obey the constraint $\sum_i c_i = 1$, which is true by construction of the limiter.
 
-% Maintaining isothermal flow requires
-% \begin{equation}
-% T = \frac{\overline{W} \bar{p}}{\rho R} = \frac{\bar{p}}{R \rho \sum_\alpha \frac{Z_\alpha}{W_\alpha}} = \frac{\bar{p}}{R \sum_\alpha \frac{\rho Z_\alpha}{W_\alpha}}
-% \end{equation}
-% to be constant and uniform at all cells and faces. Therefore, with $\bar{p}$ and $R$ constant and uniform, we must maintain
-% \begin{equation}
-% \label{eq:rho_mw}
-% \sum_\alpha \frac{(\rho Z_\alpha)}{W_\alpha} = \frac{\rho}{\overline{W}}
-% \end{equation}
+As mentioned above, constant limiter coefficients are only possible with simple schemes like Godunov (low order) or central (not TVD).  In practice, applying a second-order limiter to all species independently and then taking the worst case also reduces to a low order scheme.  Therefore, it is important to maintain the independent second-order limiters for the individual species.  An algorithm to achieve this goal and satisfy the equation of state is described next.
 
-% The above condition is automatically satisfied in the cases of using Godunov or Central differencing or in the case of binary flow (two species).  However, if we apply a second-order flux limiter, such as Superbee or CHARM, independently to each species in a multi-component (three or more species) flow with variable molecular weights, then this condition is easily violated.
+\subsubsection{Molecular Weight Correction Algorithm}
+
+Whether the flow is isothermal or not, satisfying Eq.~(\ref{eq:rho_mw}) is required to obey the equation of state.  The scheme introduced here uses a flux limiter to find a face value of $\rho/\overline{W}$ and then corrects the most abundant species face density to satisfy Eq.~(\ref{eq:rho_mw}).
+
+The following algorithm computes flux-limited species face densities for multi-component mixtures with variable molecular weights and preserves the equation of state, Eq.~(\ref{eq:rho_mw}), at cell faces.  For simplicity, the following considers a 1-D domain with cell face between cells $i$ and $i+1$.
+\begin{enumerate}
+\item Compute flux limited face values as usual, $\overline{(\rho Z_\alpha)}^{FL}_{i+1/2}$, for $\alpha=1..N_s$.
+\item Also, for each face, determine the {\tt maxloc} of the flux-limited $\overline{(\rho Z_\alpha)}^{FL}_{i+1/2}$ (the most abundant species on the face) and call it $\gamma$.  We will absorb the error into this species.
+\item For each cell, compute $\displaystyle \frac{\rho}{\overline{W}}$ from Eq.~(\ref{eq:rho_mw}).
+\item For each cell face, compute $\displaystyle \overline{\left(\frac{\rho}{\overline{W}}\right)}^{FL}_{i+1/2}$.
+\item Last, determine the face component density for $\gamma$ as
+\begin{equation}
+\overline{(\rho Z_\gamma)}^{FL}_{i+1/2} = W_\gamma \left[\overline{\left(\frac{\rho}{\overline{W}}\right)}^{FL}_{i+1/2} - \sum_{\alpha\ne\gamma} \frac{\overline{(\rho Z_{\alpha})}^{FL}_{i+1/2}}{W_\alpha} \right]
+\end{equation}
+\end{enumerate}
 
 
 \subsection{Time Splitting for Mass Source Terms}
diff --git a/Source/chem.f90 b/Source/chem.f90
index 8998f5d56c8..d7d2cc2d327 100644
--- a/Source/chem.f90
+++ b/Source/chem.f90
@@ -15,6 +15,7 @@ MODULE CVODE_INTERFACE
 USE TYPES
 USE CHEMCONS
 USE, INTRINSIC :: ISO_C_BINDING
+USE FSUNDIALS_CORE_MOD
 
 IMPLICIT NONE
 
@@ -44,7 +45,6 @@ MODULE CVODE_INTERFACE
 INTEGER(C_INT) FUNCTION RHSFN(TN_C, SUNVEC_Y, SUNVEC_F, C_USER_DATA) &
     RESULT(IERR) BIND(C,NAME='RHSFN')
 
-USE FSUNDIALS_NVECTOR_MOD
 
 ! CALLING VARIABLES
 REAL(C_DOUBLE), VALUE :: TN_C      ! CURRENT TIME
@@ -268,9 +268,7 @@ INTEGER(C_INT) FUNCTION JACFN(TN_C, SUNVEC_Y, SUNVEC_F, SUNMAT_J, &
  C_USER_DATA, TMP1, TMP2, TMP3) &
  RESULT(IERR) BIND(C,NAME='JACFN')
 
-USE FSUNDIALS_NVECTOR_MOD
 USE FSUNMATRIX_DENSE_MOD
-USE FSUNDIALS_MATRIX_MOD
 
 ! CALLING VARIABLES
 REAL(C_DOUBLE), VALUE :: TN_C             ! CURRENT TIME
@@ -810,13 +808,9 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, TMP_UNMIX, PR_IN, ZETA0, TAU_MIX, CELL_
 USE CHEMCONS, ONLY: CVODE_WARNING_CELLS,CVODE_ERR_CODE_MIN,CVODE_ERR_CODE_MAX
 USE GLOBAL_CONSTANTS
 USE FCVODE_MOD                 ! FORTRAN INTERFACE TO CVODE
-USE FSUNDIALS_CONTEXT_MOD      ! FORTRAN INTERFACE TO SUNCONTEXT
 USE FNVECTOR_SERIAL_MOD        ! FORTRAN INTERFACE TO SERIAL N_VECTOR
 USE FSUNMATRIX_DENSE_MOD       ! FORTRAN INTERFACE TO DENSE SUNMATRIX
 USE FSUNLINSOL_DENSE_MOD       ! FORTRAN INTERFACE TO DENSE SUNLINEARSOLVER
-USE FSUNDIALS_LINEARSOLVER_MOD ! FORTRAN INTERFACE TO GENERIC SUNLINEARSOLVER
-USE FSUNDIALS_MATRIX_MOD       ! FORTRAN INTERFACE TO GENERIC SUNMATRIX
-USE FSUNDIALS_NVECTOR_MOD      ! FORTRAN INTERFACE TO GENERIC N_VECTOR
 
 REAL(EB), INTENT(INOUT) :: CC(N_TRACKED_SPECIES)
 REAL(EB), INTENT(IN)    :: ZZ_0(N_TRACKED_SPECIES),TMP_IN,TMP_UNMIX,PR_IN,ZETA0,TAU_MIX,CELL_MASS,TCUR,TEND
@@ -864,7 +858,7 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, TMP_UNMIX, PR_IN, ZETA0, TAU_MIX, CELL_
 ATOLVEC_C(N_TRACKED_SPECIES+2) = 0.001_EB
 
 ! CREATE SUNDIALS CONTEXT
-IERR_C = FSUNCONTEXT_CREATE(C_NULL_PTR, SUNCTX)
+IERR_C = FSUNCONTEXT_CREATE(SUN_COMM_NULL, SUNCTX)
 
 ! CREATE SUNDIALS N_VECTOR
 SUNVEC_Y => FN_VMAKE_SERIAL(NEQ, CVEC_C, SUNCTX)
@@ -943,7 +937,7 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, TMP_UNMIX, PR_IN, ZETA0, TAU_MIX, CELL_
 ENDIF
 
 ! SET ERROR HANDLER
-IERR_C = FCVODESETERRHANDLERFN(CVODE_MEM, C_FUNLOC(FDS_CVODE_ERR_HANDLER), C_NULL_PTR)
+IERR_C = FSUNCONTEXT_PUSHERRHANDLER(SUNCTX, C_FUNLOC(FDS_CVODE_ERR_HANDLER), C_NULL_PTR)
 IF (IERR_C /= 0) THEN
    WRITE(LU_ERR,*) 'ERROR IN FCVODESETMAXNUMSTEPS, IERR = ', IERR_C, '; HALTING'
    STOP 1
@@ -1069,21 +1063,25 @@ END SUBROUTINE CVODE_SERIAL
 !> \param FUNC_NAME The functio name where error occured send by CVODE
 !> \param MESSAGE The error message
 !> \param USER_DATA User data, not used in FDS.
-SUBROUTINE FDS_CVODE_ERR_HANDLER( ERR_CODE, MOD_NAME, FUNC_NAME, MESSAGE, USER_DATA) &
+SUBROUTINE FDS_CVODE_ERR_HANDLER(LINE, FUNC_NAME, FILE_NAME, MESSAGE, ERR_CODE, USER_DATA, CTX) &
    BIND(C,NAME='FDS_CVODE_ERR_HANDLER')
-INTEGER(C_INT), VALUE :: ERR_CODE
-CHARACTER(KIND=C_CHAR),DIMENSION(*),INTENT(IN) :: MOD_NAME
-CHARACTER(KIND=C_CHAR),DIMENSION(*),INTENT(IN) :: FUNC_NAME
-CHARACTER(KIND=C_CHAR),DIMENSION(*),INTENT(IN) :: MESSAGE
-TYPE(C_PTR), VALUE    :: USER_DATA        ! USER-DEFINED DATA
+
+INTEGER(C_INT),     VALUE :: LINE
+CHARACTER(KIND=C_CHAR),DIMENSION(*),INTENT(IN):: FUNC_NAME
+CHARACTER(KIND=C_CHAR),DIMENSION(*),INTENT(IN):: FILE_NAME
+CHARACTER(KIND=C_CHAR),DIMENSION(*),INTENT(IN):: MESSAGE
+INTEGER(C_INT),     VALUE :: ERR_CODE
+TYPE(C_PTR),        VALUE :: USER_DATA
+TYPE(C_PTR),        VALUE :: CTX       ! SUNCONTEXT POINTER
+
 CHARACTER(LEN=200) :: TEMP_STRING
 LOGICAL :: FOUND_NULL
 INTEGER :: J
 
-
 IF (DEBUG) THEN
    WRITE(LU_ERR,'(A, E18.8)')" WARN: CVODE message at CFD time. CUR_CFD_TIME=", CUR_CFD_TIME
    WRITE(LU_ERR,*) '  CVODE CODE   : ', ERR_CODE
+   WRITE(LU_ERR,*) '  LINE NUMBER  : ', LINE
 
    ! Print Message
    FOUND_NULL = .FALSE.
@@ -1104,14 +1102,14 @@ SUBROUTINE FDS_CVODE_ERR_HANDLER( ERR_CODE, MOD_NAME, FUNC_NAME, MESSAGE, USER_D
    TEMP_STRING = ''
    J = 1
    DO WHILE (.NOT. FOUND_NULL)
-       IF (MOD_NAME(J) == C_NULL_CHAR) THEN
+       IF (FILE_NAME(J) == C_NULL_CHAR) THEN
            FOUND_NULL = .TRUE.
        ELSE
-           TEMP_STRING(J:J) = MOD_NAME(J)
+           TEMP_STRING(J:J) = FILE_NAME(J)
            J = J + 1
        END IF
    END DO
-   WRITE(LU_ERR,*) '  MODULE       : ', TRIM(TEMP_STRING)
+   WRITE(LU_ERR,*) '  FILE NAME       : ', TRIM(TEMP_STRING)
 
    ! Print func name
    FOUND_NULL = .FALSE.
@@ -1128,6 +1126,7 @@ SUBROUTINE FDS_CVODE_ERR_HANDLER( ERR_CODE, MOD_NAME, FUNC_NAME, MESSAGE, USER_D
    WRITE(LU_ERR,*) '  FUNCTION     : ', TRIM(TEMP_STRING)
 
    IF (.NOT. C_ASSOCIATED(USER_DATA)) WRITE(LU_ERR,*)" NO USER_DATA IS PROVIDED"
+   IF (.NOT. C_ASSOCIATED(CTX)) WRITE(LU_ERR,*)" SUNCONTEXT IS NOT VALID"
 ENDIF
 
 END SUBROUTINE FDS_CVODE_ERR_HANDLER
diff --git a/Utilities/Python/FDS_verification_script.py b/Utilities/Python/FDS_verification_script.py
index 1a4aa81bfb7..05b8444baf6 100644
--- a/Utilities/Python/FDS_verification_script.py
+++ b/Utilities/Python/FDS_verification_script.py
@@ -10,9 +10,6 @@
 print("Using:", fdsplotlib.__file__)
 # Scripts to run prior to dataplot
 
-# print("ignition_delay...");   runpy.run_path("./scripts/cantera_ignition_delay.py", run_name="__main__")
-# print("reaction_rates...");   runpy.run_path("./scripts/cantera_reaction_rates.py", run_name="__main__")
-# print("turbulent_batch_reactor...");   runpy.run_path("./scripts/cantera_turbulent_batch_reactor.py", run_name="__main__")
 print("ashrae_7...");                            runpy.run_path("./scripts/ashrae_7.py", run_name="__main__")
 print('burke_schumann...');                      runpy.run_path("./scripts/burke_schumann.py", run_name="__main__")
 print('cat_propane_depo...');                    runpy.run_path("./scripts/cat_propane_depo.py", run_name="__main__")
diff --git a/Utilities/Python/scripts/atmospheric_boundary_layer.py b/Utilities/Python/scripts/atmospheric_boundary_layer.py
index eb7977e220a..ce81c01dfc9 100644
--- a/Utilities/Python/scripts/atmospheric_boundary_layer.py
+++ b/Utilities/Python/scripts/atmospheric_boundary_layer.py
@@ -13,14 +13,14 @@
 outdir = '../../Verification/Atmospheric_Effects/'
 pltdir = '../../Manuals/FDS_User_Guide/SCRIPT_FIGURES/'
 
-git_file = outdir + 'atmospheric_boundary_layer_1_git.txt'
-version_string = fdsplotlib.get_version_string(git_file)
-
 basein  = outdir + 'atmospheric_boundary_layer'
 baseout = pltdir + 'atmospheric_boundary_layer'
 
 for i in range(1, 5):  # Loop from 1 to 4 (inclusive)
 
+    git_file = outdir + f'atmospheric_boundary_layer_{i}_git.txt'
+    version_string = fdsplotlib.get_version_string(git_file)
+
     datafile = basein + '_' + str(i)
     outfile = baseout + '_' + str(i)
 
@@ -52,6 +52,7 @@
     theta_star = u_star**2 * theta_0 / (g * kappa * L)
 
     z = np.array([z_0[i], 10*z_0[i], 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25, 30, 50, 100])
+    z = np.sort(z)
 
     # Create figure 1 for velocity
 
@@ -63,7 +64,7 @@
         psi_m = -5*z/L
         psi_h = psi_m
 
-    u = (u_star/kappa) * (np.log(z/z_0[i]) - psi_m)
+    u = (u_star/kappa) * np.maximum(np.log(z/z_0[i]) - psi_m, 0.)
     theta = theta_0 + (theta_star/kappa) * (np.log(z/z_0[i]) - psi_h)
     T = theta * (p_0 / (p_0 - rho_0*g*(z - z_0[i])))**(-0.285)
 
@@ -75,7 +76,9 @@
 
     fig = fdsplotlib.plot_to_fig(x_data=M2.iloc[:, 1].values, y_data=M2.iloc[:, 0].values, marker_style='k-', data_label='FDS',
                                  x_label='Velocity (m/s)', y_label='Height (m)',
-                                 x_min=0, x_max=u_high[i], y_min=0, y_max=100)
+                                 x_min=0, x_max=u_high[i], y_min=0, y_max=100,
+                                 revision_label=version_string,
+                                 legend_location='lower right')
 
     fdsplotlib.plot_to_fig(x_data=u, y_data=z, figure_handle=fig, marker_style='ko', data_label='M-O Theory')
 
@@ -100,7 +103,9 @@
 
     fig = fdsplotlib.plot_to_fig(x_data=M2.iloc[:,2].values, y_data=M2.iloc[:,0].values, marker_style='k-', data_label='FDS',
                                  x_label='Temperature (°C)', y_label='Height (m)', 
-                                 x_min=T_low[i], x_max=T_high[i], y_min=0, y_max=100)
+                                 x_min=T_low[i], x_max=T_high[i], y_min=0, y_max=100,
+                                 revision_label=version_string,
+                                 legend_location='lower left')
 
     fdsplotlib.plot_to_fig(x_data=T - 273.15, y_data=z, figure_handle=fig, marker_style='ko', data_label='M-O Theory')
 
diff --git a/Utilities/Python/scripts/cantera_ignition_delay.py b/Utilities/Python/scripts/cantera_ignition_delay.py
index 40004134ecb..12b9319a338 100644
--- a/Utilities/Python/scripts/cantera_ignition_delay.py
+++ b/Utilities/Python/scripts/cantera_ignition_delay.py
@@ -1,3 +1,7 @@
+"""
+This script calculates ignition delays of different chamical mechanisms using Cantera.
+This script creates the cantera_ignition_delay.csv file under Verification/Chemistry folder.
+"""
 import numpy as np
 import time
 import os
diff --git a/Utilities/Python/scripts/cantera_reaction_rates.py b/Utilities/Python/scripts/cantera_reaction_rates.py
index 0d2e5bc9e76..7f7616840f9 100644
--- a/Utilities/Python/scripts/cantera_reaction_rates.py
+++ b/Utilities/Python/scripts/cantera_reaction_rates.py
@@ -1,3 +1,7 @@
+"""
+This script calculates reaction rates using Cantera to validate the FDS Chemistry implementation for two-three steps mechanisms.
+This script creates the <caseName>_soln.csv file under Verification/Species folder.
+"""
 import numpy as np
 import os
 import cantera as ct
@@ -61,7 +65,7 @@ def calc_reaction_rate(mechanismFile, T0, P0, Y0, tstart, tend, dt, caseName, co
     stateArrReduced = stateArr[::writeInterval]
     columnNames = ['Time'] + gas.species_names  + ['Temperature'] + ['Pressure'] 
     csvdata = pd.DataFrame(stateArrReduced, columns=columnNames)
-    # save_csv_file(Species_DIR+caseName+"_soln.csv",csvdata)
+    save_csv_file(Species_DIR+caseName+"_soln.csv",csvdata)
 
 
 
diff --git a/Utilities/Python/scripts/cantera_turbulent_batch_reactor.py b/Utilities/Python/scripts/cantera_turbulent_batch_reactor.py
index 3eade7e7205..9c153b1f17b 100644
--- a/Utilities/Python/scripts/cantera_turbulent_batch_reactor.py
+++ b/Utilities/Python/scripts/cantera_turbulent_batch_reactor.py
@@ -1,9 +1,6 @@
 """
-McDermott
-21 May 2024
-
-Turbulent Batch Reactor (TBR) model
-Constant-pressure, adiabatic kinetics simulation.
+Turbulent batch reactor script using Cantera to validate the FDS turbulence batch reactor model implementation.
+This script creates the <caseName>_soln.csv file under Verification/Chemistry folder.
 """
 
 import sys
diff --git a/Utilities/Python/setup_python_env.sh b/Utilities/Python/setup_python_env.sh
index 76654cd604c..51ba198f25b 100755
--- a/Utilities/Python/setup_python_env.sh
+++ b/Utilities/Python/setup_python_env.sh
@@ -24,16 +24,18 @@ function error_exit {
 
 
 # Check Python version > 3.7
-PYTHON_VERSION=$(python3 -c 'import sys; print("{}.{}".format(sys.version_info.major, sys.version_info.minor))')
-REQUIRED_VERSION=3.7
-if (( $(echo "$PYTHON_VERSION < $REQUIRED_VERSION" | bc -l) )); then
-    echo "Python 3.7 or higher is required. Found $PYTHON_VERSION."
+major=$(python3 -c 'import sys; print(sys.version_info.major)')
+minor=$(python3 -c 'import sys; print(sys.version_info.minor)')
+
+if (( major < 3 || (major == 3 && minor < 7) )); then
+    echo "Python 3.7 or higher is required. Found $major.$minor."
     echo "Stopping python environment setup."
     return
 else
-    echo "Python version is OK: $PYTHON_VERSION"
+    echo "Python version is OK: $major.$minor"
 fi
 
+
 # Save current directory to return later
 curdir=$(pwd)