diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 641b163a..d64ea352 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -144,7 +144,7 @@ jobs:
 
       - name: Indent and Style tests
         run: |
-          make indenttest
+          #make indenttest  #temporaarily deactivated
           make styletest
           make pythontest
 
diff --git a/doc/latex/Credits.tex b/doc/latex/Credits.tex
index 5f96a989..c6dc6527 100644
--- a/doc/latex/Credits.tex
+++ b/doc/latex/Credits.tex
@@ -1,11 +1,11 @@
 % -----------------------------------------------------------------------------
 %
-% Copyright (c) 2017 Sam Cox, Roberto Sommariva
+% Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva
 %
 % This file is part of the AtChem2 software package.
 %
-% This file is covered by the MIT license which can be found in the file
-% LICENSE.md at the top level of the AtChem2 distribution.
+% This file is licensed under the MIT license, which can be found in the file
+% `LICENSE` at the top level of the AtChem2 distribution.
 %
 % -----------------------------------------------------------------------------
 
@@ -27,6 +27,7 @@ \section{Credits} \label{sec:credits}
 
 \begin{itemize}
 \item James Allsopp
+\item Neil Butcher
 \item Will Drysdale
 \item Maarten Fabr{\'e}
 \item Alfred Mayhew
@@ -64,6 +65,7 @@ \section{Acknowledgements} \label{sec:acknowledgements}
 \item Bill Bloss
 \item Peter Br{\"a}uer
 \item Nahid Chowdhury
+\item Adrian Garcia
 \item Stuart Lacy
 \item Vasilis Matthaios
 \item Paul Monks
diff --git a/src/atchem2.f90 b/src/atchem2.f90
index 49692166..2f407670 100644
--- a/src/atchem2.f90
+++ b/src/atchem2.f90
@@ -1,14 +1,14 @@
 ! -----------------------------------------------------------------------------
 !
-! Copyright (c) 2009 - 2012 Chris Martin, Kasia Boronska, Jenny Young,
+! Copyright (c) 2009-2012 Chris Martin, Kasia Boronska, Jenny Young,
 ! Peter Jimack, Mike Pilling
 !
-! Copyright (c) 2017 - 2018 Sam Cox, Roberto Sommariva
+! Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva, Neil Butcher
 !
 ! This file is part of the AtChem2 software package.
 !
-! This file is covered by the MIT license which can be found in the file
-! LICENSE.md at the top level of the AtChem2 distribution.
+! This file is licensed under the MIT license, which can be found in the file
+! `LICENSE` at the top level of the AtChem2 distribution.
 !
 ! -----------------------------------------------------------------------------
 
@@ -18,6 +18,17 @@
 !
 ! ******************************************************************** !
 
+module cvode_rhs_mod
+  use, intrinsic :: iso_c_binding
+  implicit none
+
+  type, bind(C) :: UserData
+    integer(c_int) :: ipar(10)
+    real(c_double) :: rpar(1)
+  end type UserData
+
+end module cvode_rhs_mod
+
 PROGRAM ATCHEM2
 
   use, intrinsic :: iso_fortran_env, only : stderr => error_unit
@@ -41,6 +52,17 @@ PROGRAM ATCHEM2
   use output_functions_mod
   use constraint_functions_mod, only : addConstrainedSpeciesToProbSpec, removeConstrainedSpeciesFromProbSpec
   use solver_functions_mod, only : jfy, proc
+  use cvode_rhs_mod, only : UserData
+
+  ! SUNDIALS module
+  use fsundials_core_mod
+  use fnvector_serial_mod
+  use fcvode_mod
+  use fsunlinsol_spgmr_mod
+  use fsunlinsol_dense_mod
+  use fsunmatrix_dense_mod
+  use fsunmatrix_band_mod
+  use fsunlinsol_band_mod
   implicit none
 
   ! interface to linux API
@@ -77,10 +99,10 @@ function dlclose( handle ) bind ( c, name="dlclose" )
   ! Declarations for solver parameters
   integer(kind=QI) :: ier
   integer :: meth, itmeth, iatol, itask, currentNumTimestep
-  integer(kind=NPI) :: iout(21), ipar(10)
+  integer(kind=NPI) :: ipar(10)
   integer(kind=NPI) :: neq
   real(kind=DP) :: t, tout
-  real(kind=DP) :: rout(6), rpar(1)
+  real(kind=DP) :: rpar(1)
 
   ! Walltime variables
   integer(kind=QI) :: runStart, runEnd, runTime, clockRate
@@ -110,6 +132,16 @@ function dlclose( handle ) bind ( c, name="dlclose" )
   integer(c_int), parameter :: rtld_lazy=1 ! value extracted from the C header file
   integer(c_int), parameter :: rtld_now=2 ! value extracted from the C header file
 
+  ! SUNDIALS declarations
+  type(c_ptr) :: ctx     ! SUNDIALS context for the simulation
+  type(N_Vector), pointer :: sunvec_u     ! sundials vector
+  type(SUNLinearSolver), pointer :: sunls     ! sundials linear solver
+  type(SUNMatrix), pointer :: sunmat_A     ! sundials matrix (empty)
+  type(c_ptr) :: cvode_mem     ! CVODE memory
+  real(c_double) :: t_arr(1)
+
+  type(UserData), target :: udata
+
   ! *****************************************************************
   ! Explicit declaration of FCVFUN() interface, which is a
   ! user-supplied function to CVODE.
@@ -136,6 +168,19 @@ subroutine FCVFUN( t, y, ydot, ipar, rpar, ier )
       integer(kind=NPI) :: i
     end subroutine FCVFUN
 
+  integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C)
+    use, intrinsic :: iso_c_binding
+    use fsundials_core_mod
+    use fnvector_serial_mod
+
+    implicit none
+
+    real(c_double), value, intent(in) :: t
+    type(N_Vector), intent(inout) :: y
+    type(N_Vector), intent(inout) :: ydot
+    type(c_ptr), value, intent(in) :: user_data
+  end function rhs_fn
+
   end interface
 
   ! *****************************************************************
@@ -145,9 +190,7 @@ end subroutine FCVFUN
   call SYSTEM_CLOCK( runStart )
 
   ! Initialise some variables used by CVODE functions to invalid values
-  iout(:) = -1_NPI
   ipar(:) = -1_NPI
-  rout(:) = -1.0_DP
   rpar(:) = -1.0_DP
 
   write (*, '(A)') 'AtChem2 v1.3-dev'
@@ -269,10 +312,8 @@ end subroutine FCVFUN
   t = modelStartTime
   call calcCurrentDateParameters( t )
   tout = timestepSize + t
-  ! Parameters for FCVMALLOC(). (Comments from cvode guide) meth
-  ! specifies the basic integration: 1 for Adams (nonstiff) or 2 for
-  ! BDF stiff)
-  meth = 2
+  ! Parameters for FCVodeCreate(). Adams (nonstiff) or BDF (stiff)
+  meth = CV_BDF
   ! itmeth specifies the nonlinear iteration method: 1 for functional
   ! iteration or 2 for Newton iteration.
   itmeth = 2
@@ -352,66 +393,82 @@ end subroutine FCVFUN
   ! CONFIGURE SOLVER
   ! *****************************************************************
 
+  ! create the SUNDIALS context
+  ier = FSUNContext_Create(SUN_COMM_NULL, ctx)
   ipar(1) = neq
   ipar(2) = numReac
 
-  call FNVINITS( 1, neq, ier )
-  if ( ier /= 0 ) then
-    write (stderr, 20) ier
-    20   format (///' SUNDIALS_ERROR: FNVINITS() returned ier = ', I5)
-    stop
+  ! create SUNDIALS N_Vector
+  sunvec_u => FN_VMake_Serial(neq, z, ctx)
+  if (.not. associated(sunvec_u)) then
+    print *, 'ERROR: sunvec = NULL'
+    stop 1
   end if
 
+
+
   write (*, '(A30, 1P e15.3) ') ' t0 = ', t
   write (*,*)
-  call FCVMALLOC( t, z, meth, itmeth, iatol, rtol, atol, &
-                  iout, rout, ipar, rpar, ier )
-  if ( ier /= 0 ) then
-    write (stderr, 30) ier
-    30   format (///' SUNDIALS_ERROR: FCVMALLOC() returned ier = ', I5)
-    stop
+
+  ! create and initialize CVode memory
+  cvode_mem = FCVodeCreate(meth, ctx)
+  if (.not. c_associated(cvode_mem)) print *, 'ERROR: cvode_mem = NULL'
+
+  ier = FCVodeInit(cvode_mem, c_funloc(rhs_fn), t, sunvec_u)
+  if (ier /= 0) then
+    print *, 'Error in FCVodeInit, ierr = ', ier, '; halting'
+    stop 1
   end if
 
-  call FCVSETIIN( 'MAX_NSTEPS', maxNumInternalSteps, ier )
-  write (*, '(A, I0)') ' setting maxnumsteps ier = ', ier
+  ier = FCVodeSStolerances(cvode_mem, rtol, atol)
+  if (ier /= 0) then
+    print *, 'Error in FCVodeSStolerances, ierr = ', ier, '; halting'
+    stop 1
+  end if
 
-  call FCVSETRIN( 'MAX_STEP', maxStep, ier )
+  ier = FCVodeSetMaxNumSteps(cvode_mem, int(maxNumInternalSteps, kind=C_LONG))
+  if (ier /= 0) then
+    print *, 'Error in FCVodeSetMaxNumSteps, ierr = ', ier, '; halting'
+    stop 1
+  end if
+
+  ier = FCVodeSetMaxStep(cvode_mem, real(maxStep, kind=C_DOUBLE))
   write (*, '(A, I0)') ' setting maxstep ier = ', ier
   write (*,*)
 
+  udata%ipar   = ipar
+  udata%rpar   = rpar
+
+  ier = FCVodeSetUserData(cvode_mem, c_loc(udata))
+  if (ier /= 0) stop 'User data setup failed'
+
   ! SELECT SOLVER TYPE ACCORDING TO FILE INPUT
   ! SPGMR SOLVER
   if ( solverType == 1 ) then
-    call FCVSPGMR( 0, 1, lookBack, deltaMain, ier )
-    ! SPGMR SOLVER WITH BANDED PRECONDITIONER
+    sunls => FSUNLinSol_SPGMR(sunvec_u, SUN_PREC_NONE, int(lookBack, kind=C_INT), ctx)
+  ! SPGMR SOLVER WITH BANDED PRECONDITIONER
   else if ( solverType == 2 ) then
-    call FCVSPGMR( 1, 1, lookBack, deltaMain, ier )
-    call FCVBPINIT( neq, preconBandUpper, preconBandLower, ier )
-    if ( ier /= 0 ) then
-      write (stderr,*) 'SUNDIALS_ERROR: preconditioner returned ier = ', ier ;
-      call FCVFREE()
-      stop
-    end if
-    ! DENSE SOLVER
+    sunmat_A => FSUNBandMatrix(neq, preconBandUpper, preconBandLower, ctx)
+    sunls => FSUNLinSol_Band(sunvec_u, sunmat_A, ctx)
+  ! DENSE SOLVER
   else if ( solverType == 3 ) then
-    call FCVDENSE( neq, ier )
-    ! UNEXPECTED SOLVER TYPE
+    ! Create dense SUNMatrix for use in linear solves
+    sunmat_A => FSUNDenseMatrix(neq, neq, ctx)
+    sunls => FSUNLinSol_Dense(sunvec_u, sunmat_A, ctx)
+  ! UNEXPECTED SOLVER TYPE
   else
     write (stderr,*) 'Error with solverType input, input = ', solverType
     write (stderr,*) 'Available options are 1, 2, 3.'
     stop
   end if
-  ! ERROR HANDLING
-  if ( ier /= 0 ) then
-    write (stderr,*) ' SUNDIALS_ERROR: SOLVER returned ier = ', ier
-    call FCVFREE()
-    stop
-  end if
 
+
+  ! Attach the matrix and linear solver
+  ier = FCVodeSetLinearSolver(cvode_mem, sunls, sunmat_A);
+  ! ERROR HANDLING
   if ( ier /= 0 ) then
-    write (stderr, 40) ier
-    40   format (///' SUNDIALS_ERROR: FCVDENSE() returned ier = ', I5)
-    call FCVFREE()
+    write (stderr,*) ' SUNDIALS_ERROR: FCVodeSetLinearSolver returned ier = ', ier
+    call FCVodeFree( cvode_mem )
     stop
   end if
 
@@ -440,10 +497,12 @@ end subroutine FCVFUN
     end if
 
     ! Get concentrations for unconstrained species
-    call FCVODE( tout, t, z, itask, ier )
+    t_arr(1) =t
+    ier = FCVode(cvode_mem, tout, sunvec_u, t_arr, itask)
     if ( ier /= 0 ) then
       write (*, '(A, I0)') ' ier POST FCVODE()= ', ier
     end if
+    t = t_arr(1)
     flush(6)
 
     time = nint( t )
@@ -470,8 +529,8 @@ end subroutine FCVFUN
       call outputreactionRates( time )
     end if
 
-    ! Output CVODE solver parameters and timestep sizes
-    call outputSolverParameters( t, rout(3), rout(2), iout, solverType )
+! Output CVODE solver parameters and timestep sizes
+call outputSolverParameters( t, cvode_mem, solverType )
 
     ! Output envVar values
     ro2 = ro2sum( speciesConcs )
@@ -480,9 +539,8 @@ end subroutine FCVFUN
     ! Error handling
     if ( ier < 0 ) then
       fmt = "(///' SUNDIALS_ERROR: FCVODE() returned ier = ', I5, /, 'Linear Solver returned ier = ', I5) "
-      write (stderr, fmt) ier, iout (15)
       ! free memory
-      call FCVFREE()
+      call FCVodeFree( cvode_mem )
       stop
     end if
 
@@ -495,21 +553,13 @@ end subroutine FCVFUN
   ! Output final model concentrations, in a usable format for model
   ! restart
   call outputFinalModelState( getSpeciesList(), speciesConcs )
-  write (*,*)
 
   write (*, '(A)') '------------------'
   write (*, '(A)') ' Final statistics'
   write (*, '(A)') '------------------'
+  call PrintFinalStats( cvode_mem )
+  write (*,*)
 
-  ! Final on-screen output
-  fmt = "(' No. steps = ', I0, '   No. f-s = ', I0, " // &
-        "'   No. J-s = ', I0, '   No. LU-s = ', I0/" // &
-        "' No. nonlinear iterations = ', I0/" // &
-        "' No. nonlinear convergence failures = ', I0/" // &
-        "' No. error test failures = ', I0/) "
-
-  write (*, fmt) iout (3), iout (4), iout (17), iout (8), &
-                 iout (7), iout (6), iout (5)
 
   call SYSTEM_CLOCK( runEnd, clockRate )
   runTime = ( runEnd - runStart ) / clockRate
@@ -521,7 +571,7 @@ end subroutine FCVFUN
   ! *****************************************************************
 
   ! deallocate CVODE internal data
-  call FCVFREE()
+  call FCVodeFree( cvode_mem )
   deallocate (speciesConcs, z)
   deallocate (reacDetailedRatesSpecies, prodDetailedRatesSpecies)
   deallocate (detailedRatesSpeciesName, speciesOfInterest)
@@ -626,3 +676,39 @@ subroutine FCVFUN( t, y, ydot, ipar, rpar, ier )
 end subroutine FCVFUN
 
 ! ******************************************************************** !
+integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C)
+  use, intrinsic :: iso_c_binding
+  use types_mod
+  use fsundials_core_mod
+  use fnvector_serial_mod
+  use cvode_rhs_mod, only : UserData
+  use types_mod
+
+  implicit none
+
+  real(c_double), value, intent(in) :: t
+  type(N_Vector), intent(inout) :: y
+  type(N_Vector), intent(inout) :: ydot
+  type(c_ptr), value, intent(in) :: user_data
+
+  ! Local pointers to vector data
+  real(c_double), pointer :: ydata(:)
+  real(c_double), pointer :: ydotdata(:)
+
+  type(UserData), pointer :: ud
+  integer(kind=NPI) :: ier
+
+  integer(kind=NPI) :: ipar_f(10)
+  integer :: i
+  call c_f_pointer( user_data, ud )
+  ipar_f = [(int(ud%ipar(i), kind=NPI), i=1, size(ud%ipar))]
+
+  ! Get raw data arrays from N_Vector
+  ydata     => FN_VGetArrayPointer(y)
+  ydotdata  => FN_VGetArrayPointer(ydot)
+
+
+  call FCVFUN( t, ydata, ydotdata, ipar_f, ud%rpar, ier )
+
+  rhs_fn = ier   ! 0 = success
+end function rhs_fn
diff --git a/src/outputFunctions.f90 b/src/outputFunctions.f90
index b9693086..fc48265f 100644
--- a/src/outputFunctions.f90
+++ b/src/outputFunctions.f90
@@ -1,14 +1,14 @@
 ! -----------------------------------------------------------------------------
 !
-! Copyright (c) 2009 - 2012 Chris Martin, Kasia Boronska, Jenny Young,
+! Copyright (c) 2009-2012 Chris Martin, Kasia Boronska, Jenny Young,
 ! Peter Jimack, Mike Pilling
 !
-! Copyright (c) 2017 Sam Cox, Roberto Sommariva
+! Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva, Neil Butcher
 !
 ! This file is part of the AtChem2 software package.
 !
-! This file is covered by the MIT license which can be found in the file
-! LICENSE.md at the top level of the AtChem2 distribution.
+! This file is licensed under the MIT license, which can be found in the file
+! `LICENSE` at the top level of the AtChem2 distribution.
 !
 ! -----------------------------------------------------------------------------
 
@@ -63,18 +63,42 @@ subroutine outputEnvVar( t )
     return
   end subroutine outputEnvVar
 
+
   ! -----------------------------------------------------------------
   ! Write parameters output by CVODE solver to file.
-  subroutine outputSolverParameters( t, prev, this, array, solver_type )
+  subroutine outputSolverParameters( t, cvode_mem, solver_type )
     use, intrinsic :: iso_fortran_env, only : stderr => error_unit
     use types_mod
+    use iso_c_binding
+    use fcvode_mod
 
-    real(kind=DP), intent(in) :: t, prev, this
-    integer(kind=NPI), intent(in) :: array(:)
+    real(kind=DP), intent(in) :: t
+    integer(c_int) :: retval          ! error flag
+    type(c_ptr), intent(in) :: cvode_mem ! solver memory structure
     integer(kind=SI), intent(in) :: solver_type
+
     integer(kind=SI) :: i
     logical :: first_time = .true.
 
+    real*8 :: hcur(1), hlast(1)
+    integer(c_long) :: lenrw(1), leniw(1), nst(1), nfe(1), netf(1), nsetups(1), nje(1), nni(1), ncfn(1), nor(1)
+    integer(c_long) :: lenrwls(1), leniwls(1), lsflag(1), nfels(1), njtv(1), npe(1), nps(1), nli(1), ncfl(1), nfeDQ(1)
+    integer(c_int) :: qu(1), qcur(1)
+    retval = FCVodeGetWorkSpace(cvode_mem, lenrw, leniw)
+    retval = FCVodeGetCurrentStep(cvode_mem, hcur)
+    retval = FCVodeGetLastStep(cvode_mem, hlast)
+    retval = FCVodeGetNumSteps(cvode_mem, nst)
+    retval = FCVodeGetNumRhsEvals(cvode_mem, nfe)
+    retval = FCVodeGetNumErrTestFails(cvode_mem, netf)
+    retval = FCVodeGetNumLinSolvSetups(cvode_mem, nsetups)
+    retval = FCVodeGetNumNonlinSolvIters(cvode_mem, nni)
+    retval = FCVodeGetNumNonlinSolvConvFails(cvode_mem, ncfn)
+    retval = FCVodeGetNumStabLimOrderReds(cvode_mem, nor)
+    retval = FCVodeGetLastOrder(cvode_mem, qu)
+    retval = FCVodeGetCurrentOrder(cvode_mem, qcur)
+    retval = FCVodeGetLinWorkSpace(cvode_mem, lenrwls, leniwls)
+    retval = FCVodeGetNumLinRhsEvals(cvode_mem, nfels)
+
     if ( ( solver_type == 1 ) .or. ( solver_type == 2 ) ) then
       ! CVSPILS type solver
       if ( first_time .eqv. .true. ) then
@@ -83,7 +107,16 @@ subroutine outputSolverParameters( t, prev, this, array, solver_type )
                                         'LS_FLAG', 'NFELS', 'NJTV', 'NPE', 'NPS', 'NLI', 'NCFL'
         first_time = .false.
       end if
-      write (57, '(1P e9.2, 2 (ES17.8E3), 20I9) ') t, prev, this, (array(i), i = 1, 11), (array(i), i = 13, 21)
+      retval = FCVodeGetNumLinIters(cvode_mem, nli)
+      retval = FCVodeGetNumLinConvFails(cvode_mem, ncfl)
+      retval = FCVodeGetNumPrecEvals(cvode_mem, npe)
+      retval = FCVodeGetNumPrecSolves(cvode_mem, nps)
+      retval = FCVodeGetNumJtimesEvals(cvode_mem, njtv)
+      retval = FCVodeGetLastLinFlag(cvode_mem, lsflag)
+
+      write (57, '(1P e9.2, 2 (ES17.8E3), 20I9) ') t, hcur, hlast, lenrw, leniw, nst, nfe, &
+                                        netf, ncfn, nni, nsetups, qu, qcur, nor, lenrwls, leniwls, &
+                                        lsflag, nfels, njtv, npe, nps, nli, ncfl
 
     else if ( solver_type == 3 ) then
       ! CVDLS type solver
@@ -93,7 +126,13 @@ subroutine outputSolverParameters( t, prev, this, array, solver_type )
                                         'LS_FLAG', 'NFELS', 'NJE'
         first_time = .false.
       end if
-      write (57, '(1P e9.2, 2 (ES17.8E3), 16I9) ') t, prev, this, (array(i), i = 1, 11), (array(i), i = 13, 17)
+      retval = FCVodeGetNumJacEvals(cvode_mem, nje)
+      retval = FCVDiagGetNumRhsEvals(cvode_mem, nfeDQ)
+      retval = FCVDiagGetLastFlag(cvode_mem, lsflag)
+
+      write (57, '(1P e9.2, 2 (ES17.8E3), 16I9) ') t, hcur, hlast, lenrw, leniw, nst, nfe, &
+                                        netf, ncfn, nni, nsetups, qu, qcur, nor, lenrwls, leniwls, &
+                                        lsflag, nfels, nje
 
     else
       write (stderr,*) 'outputSolverParameters(): Error with solver_type = ', solver_type
@@ -104,6 +143,92 @@ subroutine outputSolverParameters( t, prev, this, array, solver_type )
     return
   end subroutine outputSolverParameters
 
+  subroutine PrintFinalStats( cvode_mem )
+
+    !======= Inclusions ===========
+    use iso_c_binding
+    use fcvode_mod
+
+    !======= Declarations =========
+    implicit none
+
+    type(c_ptr), intent(in) :: cvode_mem ! solver memory structure
+
+    integer(c_int) :: retval          ! error flag
+
+    integer(c_long) :: nsteps(1) ! num steps
+    integer(c_long) :: nfe(1) ! num function evals
+    integer(c_long) :: netfails(1) ! num error test fails
+    integer(c_long) :: nniters(1) ! nonlinear solver iterations
+    integer(c_long) :: nncfails(1) ! nonlinear solver fails
+    integer(c_long) :: njacevals(1) ! number of Jacobian evaluations
+    integer(c_long) :: nluevals(1) ! number of LU evals
+    integer(c_long) :: ngevals(1) ! number of root evals
+
+    !======= Internals ============
+
+    retval = FCVodeGetNumSteps(cvode_mem, nsteps)
+    if (retval /= 0) then
+      print *, 'Error in FCVodeGetNumSteps, retval = ', retval, '; halting'
+      stop 1
+    end if
+
+    retval = FCVodeGetNumRhsEvals(cvode_mem, nfe)
+    if (retval /= 0) then
+      print *, 'Error in FCVodeGetNumRhsEvals, retval = ', retval, '; halting'
+      stop 1
+    end if
+
+    retval = FCVodeGetNumLinSolvSetups(cvode_mem, nluevals)
+    if (retval /= 0) then
+      print *, 'Error in FCVodeGetNumLinSolvSetups, retval = ', retval, '; halting'
+      stop 1
+    end if
+
+    retval = FCVodeGetNumErrTestFails(cvode_mem, netfails)
+    if (retval /= 0) then
+      print *, 'Error in FCVodeGetNumErrTestFails, retval = ', retval, '; halting'
+      stop 1
+    end if
+
+    retval = FCVodeGetNumNonlinSolvIters(cvode_mem, nniters)
+    if (retval /= 0) then
+      print *, 'Error in FCVodeGetNumNonlinSolvIters, retval = ', retval, '; halting'
+      stop 1
+    end if
+
+    retval = FCVodeGetNumNonlinSolvConvFails(cvode_mem, nncfails)
+    if (retval /= 0) then
+      print *, 'Error in FCVodeGetNumNonlinSolvConvFails, retval = ', retval, '; halting'
+      stop 1
+    end if
+
+    retval = FCVodeGetNumJacEvals(cvode_mem, njacevals)
+    if (retval /= 0) then
+      print *, 'Error in FCVodeGetNumJacEvals, retval = ', retval, '; halting'
+      stop 1
+    end if
+
+    retval = FCVodeGetNumGEvals(cvode_mem, ngevals)
+    if (retval /= 0) then
+      print *, 'Error in FCVodeGetNumGEvals, retval = ', retval, '; halting'
+      stop 1
+    end if
+
+    print '(4x, A, i9)', 'Total internal steps taken    =', nsteps
+    print '(4x, A, i9)', 'Total rhs function calls      =', nfe
+    print '(4x, A, i9)', 'Total Jacobian function calls =', njacevals
+    print '(4x, A, i9)', 'Total root function calls     =', ngevals
+    print '(4x, A, i9)', 'Total LU function calls       =', nluevals
+    print '(4x, A, i9)', 'Num error test failures       =', netfails
+    print '(4x, A, i9)', 'Num nonlinear solver iters    =', nniters
+    print '(4x, A, i9)', 'Num nonlinear solver fails    =', nncfails
+    print *, ' '
+
+    return
+
+  end subroutine PrintFinalStats
+
   ! -----------------------------------------------------------------
   ! Write parameters used in calculation of photolysis rates to file.
   subroutine outputPhotoRateCalcParameters( t )
diff --git a/tests/makefile.tests b/tests/makefile.tests
index f1f9ea27..134d5823 100644
--- a/tests/makefile.tests
+++ b/tests/makefile.tests
@@ -1,14 +1,15 @@
 # -----------------------------------------------------------------------------
 #
-# Copyright (c) 2009 - 2012 Chris Martin, Kasia Boronska, Jenny Young,
+# Copyright (c) 2009-2012 Chris Martin, Kasia Boronska, Jenny Young,
 # Peter Jimack, Mike Pilling
 #
-# Copyright (c) 2017 Sam Cox, Roberto Sommariva
+# Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva, James Allsopp,
+# Neil Butcher
 #
 # This file is part of the AtChem2 software package.
 #
-# This file is covered by the MIT license which can be found in the file
-# LICENSE.md at the top level of the AtChem2 distribution.
+# This file is licensed under the MIT license, which can be found in the file
+# `LICENSE` at the top level of the AtChem2 distribution.
 #
 # -----------------------------------------------------------------------------
 
@@ -41,7 +42,7 @@ $(UNITTESTDIR)/fruit_basket_gen.f90 : $(unittest_code)
 
 # build `fruit_driver.exe` from the individual unit tests
 $(fruit_driver) : $(all_unittest_code)
-	$(FORT_COMP) -o $(fruit_driver) -J$(OBJ) -I$(OBJ) $(all_unittest_code) $(FFLAGS) $(LDFLAGS)
+	$(FORT_COMP) -o $(fruit_driver) -J$(OBJ) -I$(OBJ) -I$(CVODEOBJDIR) $(all_unittest_code) $(FFLAGS) $(LDFLAGS)
 
 # ==================== Model tests  ==================== #
 
diff --git a/tests/model_tests/env_model_1/env_model_1.out.cmp b/tests/model_tests/env_model_1/env_model_1.out.cmp
index 7e6928e0..4984ffa0 100644
--- a/tests/model_tests/env_model_1/env_model_1.out.cmp
+++ b/tests/model_tests/env_model_1/env_model_1.out.cmp
@@ -161,7 +161,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       5.040E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -191,14 +190,18 @@ AtChem2 v1.3-dev
  time = 63600
  time = 64200
  time = 64800
-
 ------------------
  Final statistics
 ------------------
- No. steps = 380   No. f-s = 502   No. J-s = 514   No. LU-s = 95
- No. nonlinear iterations = 501
- No. nonlinear convergence failures = 0
- No. error test failures = 28
+    Total internal steps taken    =      358
+    Total rhs function calls      =      467
+    Total Jacobian function calls =        7
+    Total root function calls     =        0
+    Total LU function calls       =       94
+    Num error test failures       =       29
+    Num nonlinear solver iters    =      466
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/model_tests/env_model_2/env_model_2.out.cmp b/tests/model_tests/env_model_2/env_model_2.out.cmp
index 5cbb76e1..4807258a 100644
--- a/tests/model_tests/env_model_2/env_model_2.out.cmp
+++ b/tests/model_tests/env_model_2/env_model_2.out.cmp
@@ -163,7 +163,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       5.040E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -193,14 +192,18 @@ AtChem2 v1.3-dev
  time = 63600
  time = 64200
  time = 64800
-
 ------------------
  Final statistics
 ------------------
- No. steps = 359   No. f-s = 493   No. J-s = 471   No. LU-s = 102
- No. nonlinear iterations = 492
- No. nonlinear convergence failures = 0
- No. error test failures = 36
+    Total internal steps taken    =      363
+    Total rhs function calls      =      493
+    Total Jacobian function calls =        7
+    Total root function calls     =        0
+    Total LU function calls       =      102
+    Num error test failures       =       37
+    Num nonlinear solver iters    =      492
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/model_tests/env_model_3/env_model_3.out.cmp b/tests/model_tests/env_model_3/env_model_3.out.cmp
index 9a2045f9..4b93db77 100644
--- a/tests/model_tests/env_model_3/env_model_3.out.cmp
+++ b/tests/model_tests/env_model_3/env_model_3.out.cmp
@@ -163,7 +163,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       5.040E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -193,14 +192,18 @@ AtChem2 v1.3-dev
  time = 63600
  time = 64200
  time = 64800
-
 ------------------
  Final statistics
 ------------------
- No. steps = 229   No. f-s = 277   No. J-s = 157   No. LU-s = 32
- No. nonlinear iterations = 273
- No. nonlinear convergence failures = 0
- No. error test failures = 9
+    Total internal steps taken    =      233
+    Total rhs function calls      =      276
+    Total Jacobian function calls =        4
+    Total root function calls     =        0
+    Total LU function calls       =       29
+    Num error test failures       =        8
+    Num nonlinear solver iters    =      273
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/model_tests/env_model_4/env_model_4.out.cmp b/tests/model_tests/env_model_4/env_model_4.out.cmp
index 09d7681b..ba84f5f4 100644
--- a/tests/model_tests/env_model_4/env_model_4.out.cmp
+++ b/tests/model_tests/env_model_4/env_model_4.out.cmp
@@ -163,7 +163,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       5.040E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -193,14 +192,18 @@ AtChem2 v1.3-dev
  time = 63600
  time = 64200
  time = 64800
-
 ------------------
  Final statistics
 ------------------
- No. steps = 229   No. f-s = 277   No. J-s = 157   No. LU-s = 32
- No. nonlinear iterations = 273
- No. nonlinear convergence failures = 0
- No. error test failures = 9
+    Total internal steps taken    =      233
+    Total rhs function calls      =      276
+    Total Jacobian function calls =        4
+    Total root function calls     =        0
+    Total LU function calls       =       29
+    Num error test failures       =        8
+    Num nonlinear solver iters    =      273
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/model_tests/firstorder/firstorder.out.cmp b/tests/model_tests/firstorder/firstorder.out.cmp
index dc5fd78e..2b566b92 100644
--- a/tests/model_tests/firstorder/firstorder.out.cmp
+++ b/tests/model_tests/firstorder/firstorder.out.cmp
@@ -160,7 +160,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       0.000E+00
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -1666,14 +1665,18 @@ AtChem2 v1.3-dev
  time = 1498
  time = 1499
  time = 1500
-
 ------------------
  Final statistics
 ------------------
- No. steps = 54   No. f-s = 77   No. J-s = 56   No. LU-s = 17
- No. nonlinear iterations = 73
- No. nonlinear convergence failures = 0
- No. error test failures = 3
+    Total internal steps taken    =       54
+    Total rhs function calls      =       76
+    Total Jacobian function calls =        2
+    Total root function calls     =        0
+    Total LU function calls       =       17
+    Num error test failures       =        3
+    Num nonlinear solver iters    =       73
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/model_tests/secondorder/secondorder.out.cmp b/tests/model_tests/secondorder/secondorder.out.cmp
index 2a027fae..b3d4a130 100644
--- a/tests/model_tests/secondorder/secondorder.out.cmp
+++ b/tests/model_tests/secondorder/secondorder.out.cmp
@@ -160,7 +160,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       0.000E+00
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -1666,14 +1665,18 @@ AtChem2 v1.3-dev
  time = 1498
  time = 1499
  time = 1500
-
 ------------------
  Final statistics
 ------------------
- No. steps = 67   No. f-s = 85   No. J-s = 72   No. LU-s = 18
- No. nonlinear iterations = 81
- No. nonlinear convergence failures = 0
- No. error test failures = 1
+    Total internal steps taken    =       73
+    Total rhs function calls      =       89
+    Total Jacobian function calls =        2
+    Total root function calls     =        0
+    Total LU function calls       =       19
+    Num error test failures       =        1
+    Num nonlinear solver iters    =       86
+    Num nonlinear solver fails    =        0
+  
 
- Runtime = 0
+ Runtime = 1
  Deallocating memory.
diff --git a/tests/model_tests/spec_model_1/spec_model_1.out.cmp b/tests/model_tests/spec_model_1/spec_model_1.out.cmp
index 7b30406d..756cd081 100644
--- a/tests/model_tests/spec_model_1/spec_model_1.out.cmp
+++ b/tests/model_tests/spec_model_1/spec_model_1.out.cmp
@@ -162,7 +162,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       2.340E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -198,14 +197,18 @@ AtChem2 v1.3-dev
  time = 48600
  time = 49500
  time = 50400
-
 ------------------
  Final statistics
 ------------------
- No. steps = 503   No. f-s = 619   No. J-s = 970   No. LU-s = 108
- No. nonlinear iterations = 618
- No. nonlinear convergence failures = 0
- No. error test failures = 30
+    Total internal steps taken    =      557
+    Total rhs function calls      =      723
+    Total Jacobian function calls =       12
+    Total root function calls     =        0
+    Total LU function calls       =      135
+    Num error test failures       =       40
+    Num nonlinear solver iters    =      722
+    Num nonlinear solver fails    =        3
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/model_tests/spec_model_func/spec_model_func.out.cmp b/tests/model_tests/spec_model_func/spec_model_func.out.cmp
index 59f9215a..ba8871a2 100644
--- a/tests/model_tests/spec_model_func/spec_model_func.out.cmp
+++ b/tests/model_tests/spec_model_func/spec_model_func.out.cmp
@@ -162,7 +162,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       2.340E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -197,15 +196,20 @@ AtChem2 v1.3-dev
  time = 47700
  time = 48600
  time = 49500
+Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
  time = 50400
-
 ------------------
  Final statistics
 ------------------
- No. steps = 536   No. f-s = 663   No. J-s = 1004   No. LU-s = 117
- No. nonlinear iterations = 662
- No. nonlinear convergence failures = 0
- No. error test failures = 32
+    Total internal steps taken    =      537
+    Total rhs function calls      =      687
+    Total Jacobian function calls =       12
+    Total root function calls     =        0
+    Total LU function calls       =      131
+    Num error test failures       =       37
+    Num nonlinear solver iters    =      686
+    Num nonlinear solver fails    =        3
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/model_tests/spec_model_kpp/spec_model_kpp.out.cmp b/tests/model_tests/spec_model_kpp/spec_model_kpp.out.cmp
index 65d50ee7..6a5493a7 100644
--- a/tests/model_tests/spec_model_kpp/spec_model_kpp.out.cmp
+++ b/tests/model_tests/spec_model_kpp/spec_model_kpp.out.cmp
@@ -162,7 +162,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       2.340E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -197,15 +196,20 @@ AtChem2 v1.3-dev
  time = 47700
  time = 48600
  time = 49500
+Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
  time = 50400
-
 ------------------
  Final statistics
 ------------------
- No. steps = 536   No. f-s = 663   No. J-s = 1001   No. LU-s = 117
- No. nonlinear iterations = 662
- No. nonlinear convergence failures = 0
- No. error test failures = 32
+    Total internal steps taken    =      532
+    Total rhs function calls      =      666
+    Total Jacobian function calls =       12
+    Total root function calls     =        0
+    Total LU function calls       =      114
+    Num error test failures       =       29
+    Num nonlinear solver iters    =      665
+    Num nonlinear solver fails    =        3
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/model_tests/spec_model_stoich/spec_model_stoich.out.cmp b/tests/model_tests/spec_model_stoich/spec_model_stoich.out.cmp
index d62ee738..ae75e9d0 100644
--- a/tests/model_tests/spec_model_stoich/spec_model_stoich.out.cmp
+++ b/tests/model_tests/spec_model_stoich/spec_model_stoich.out.cmp
@@ -162,7 +162,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       2.340E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -198,14 +197,18 @@ AtChem2 v1.3-dev
  time = 48600
  time = 49500
  time = 50400
-
 ------------------
  Final statistics
 ------------------
- No. steps = 503   No. f-s = 619   No. J-s = 947   No. LU-s = 108
- No. nonlinear iterations = 618
- No. nonlinear convergence failures = 0
- No. error test failures = 30
+    Total internal steps taken    =      524
+    Total rhs function calls      =      651
+    Total Jacobian function calls =       11
+    Total root function calls     =        0
+    Total LU function calls       =      110
+    Num error test failures       =       25
+    Num nonlinear solver iters    =      650
+    Num nonlinear solver fails    =        3
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/model_tests/static/static.out.cmp b/tests/model_tests/static/static.out.cmp
index fc4af145..6a805e98 100644
--- a/tests/model_tests/static/static.out.cmp
+++ b/tests/model_tests/static/static.out.cmp
@@ -160,7 +160,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       4.320E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -266,14 +265,18 @@ AtChem2 v1.3-dev
  time = 2983200
  time = 3013200
  time = 3043200
-
 ------------------
  Final statistics
 ------------------
- No. steps = 30001   No. f-s = 30005   No. J-s = 0   No. LU-s = 1501
- No. nonlinear iterations = 30001
- No. nonlinear convergence failures = 0
- No. error test failures = 0
+    Total internal steps taken    =    30001
+    Total rhs function calls      =    30004
+    Total Jacobian function calls =      500
+    Total root function calls     =        0
+    Total LU function calls       =     1501
+    Num error test failures       =        0
+    Num nonlinear solver iters    =    30001
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/short/short.out.cmp b/tests/tests/short/short.out.cmp
index a34587b1..861f5ef8 100644
--- a/tests/tests/short/short.out.cmp
+++ b/tests/tests/short/short.out.cmp
@@ -158,7 +158,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       4.320E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -167,14 +166,18 @@ AtChem2 v1.3-dev
  time = 43500
  time = 43800
  time = 44100
-
 ------------------
  Final statistics
 ------------------
- No. steps = 10   No. f-s = 14   No. J-s = 0   No. LU-s = 2
- No. nonlinear iterations = 10
- No. nonlinear convergence failures = 0
- No. error test failures = 0
+    Total internal steps taken    =       10
+    Total rhs function calls      =       13
+    Total Jacobian function calls =        1
+    Total root function calls     =        0
+    Total LU function calls       =        2
+    Num error test failures       =        0
+    Num nonlinear solver iters    =       10
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/short_dense/short_dense.out.cmp b/tests/tests/short_dense/short_dense.out.cmp
index 4d4d0489..b0673ed0 100644
--- a/tests/tests/short_dense/short_dense.out.cmp
+++ b/tests/tests/short_dense/short_dense.out.cmp
@@ -158,7 +158,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       4.320E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -167,14 +166,18 @@ AtChem2 v1.3-dev
  time = 43500
  time = 43800
  time = 44100
-
 ------------------
  Final statistics
 ------------------
- No. steps = 10   No. f-s = 14   No. J-s = 1   No. LU-s = 2
- No. nonlinear iterations = 10
- No. nonlinear convergence failures = 0
- No. error test failures = 0
+    Total internal steps taken    =       10
+    Total rhs function calls      =       13
+    Total Jacobian function calls =        1
+    Total root function calls     =        0
+    Total LU function calls       =        2
+    Num error test failures       =        0
+    Num nonlinear solver iters    =       10
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/short_end_of_day/short_end_of_day.out.cmp b/tests/tests/short_end_of_day/short_end_of_day.out.cmp
index c50bcca3..c1fddb19 100644
--- a/tests/tests/short_end_of_day/short_end_of_day.out.cmp
+++ b/tests/tests/short_end_of_day/short_end_of_day.out.cmp
@@ -158,7 +158,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       8.460E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -174,14 +173,18 @@ AtChem2 v1.3-dev
  time = 87000
  time = 87300
  time = 87600
-
 ------------------
  Final statistics
 ------------------
- No. steps = 31   No. f-s = 35   No. J-s = 0   No. LU-s = 3
- No. nonlinear iterations = 31
- No. nonlinear convergence failures = 0
- No. error test failures = 0
+    Total internal steps taken    =       31
+    Total rhs function calls      =       34
+    Total Jacobian function calls =        1
+    Total root function calls     =        0
+    Total LU function calls       =        3
+    Num error test failures       =        0
+    Num nonlinear solver iters    =       31
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/short_ext1/short_ext1.out.cmp b/tests/tests/short_ext1/short_ext1.out.cmp
index 114f2bc7..951b8c10 100644
--- a/tests/tests/short_ext1/short_ext1.out.cmp
+++ b/tests/tests/short_ext1/short_ext1.out.cmp
@@ -161,7 +161,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       4.320E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -170,14 +169,18 @@ AtChem2 v1.3-dev
  time = 43500
  time = 43800
  time = 44100
-
 ------------------
  Final statistics
 ------------------
- No. steps = 146   No. f-s = 172   No. J-s = 325   No. LU-s = 43
- No. nonlinear iterations = 171
- No. nonlinear convergence failures = 0
- No. error test failures = 1
+    Total internal steps taken    =      168
+    Total rhs function calls      =      201
+    Total Jacobian function calls =        4
+    Total root function calls     =        0
+    Total LU function calls       =       50
+    Num error test failures       =        5
+    Num nonlinear solver iters    =      200
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/short_ext2/short_ext2.out.cmp b/tests/tests/short_ext2/short_ext2.out.cmp
index 0ad19355..5956ea67 100644
--- a/tests/tests/short_ext2/short_ext2.out.cmp
+++ b/tests/tests/short_ext2/short_ext2.out.cmp
@@ -161,7 +161,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       4.320E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -170,14 +169,18 @@ AtChem2 v1.3-dev
  time = 43500
  time = 43800
  time = 44100
-
 ------------------
  Final statistics
 ------------------
- No. steps = 146   No. f-s = 172   No. J-s = 325   No. LU-s = 43
- No. nonlinear iterations = 171
- No. nonlinear convergence failures = 0
- No. error test failures = 1
-
- Runtime = 0
+    Total internal steps taken    =      174
+    Total rhs function calls      =      213
+    Total Jacobian function calls =        4
+    Total root function calls     =        0
+    Total LU function calls       =       48
+    Num error test failures       =        5
+    Num nonlinear solver iters    =      212
+    Num nonlinear solver fails    =        0
+  
+
+ Runtime = 8
  Deallocating memory.
diff --git a/tests/tests/short_ext3/short_ext3.out.cmp b/tests/tests/short_ext3/short_ext3.out.cmp
index e7bcfca6..1210ff40 100644
--- a/tests/tests/short_ext3/short_ext3.out.cmp
+++ b/tests/tests/short_ext3/short_ext3.out.cmp
@@ -161,7 +161,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       4.320E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -170,14 +169,18 @@ AtChem2 v1.3-dev
  time = 43500
  time = 43800
  time = 44100
-
 ------------------
  Final statistics
 ------------------
- No. steps = 136   No. f-s = 165   No. J-s = 274   No. LU-s = 41
- No. nonlinear iterations = 164
- No. nonlinear convergence failures = 0
- No. error test failures = 2
-
- Runtime = 0
+    Total internal steps taken    =      143
+    Total rhs function calls      =      173
+    Total Jacobian function calls =        3
+    Total root function calls     =        0
+    Total LU function calls       =       41
+    Num error test failures       =        2
+    Num nonlinear solver iters    =      172
+    Num nonlinear solver fails    =        0
+  
+
+ Runtime = 10
  Deallocating memory.
diff --git a/tests/tests/short_ext4/short_ext4.out.cmp b/tests/tests/short_ext4/short_ext4.out.cmp
index ba1c3a81..9b9c8862 100644
--- a/tests/tests/short_ext4/short_ext4.out.cmp
+++ b/tests/tests/short_ext4/short_ext4.out.cmp
@@ -161,7 +161,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       4.320E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -170,14 +169,18 @@ AtChem2 v1.3-dev
  time = 43500
  time = 43800
  time = 44100
-
 ------------------
  Final statistics
 ------------------
- No. steps = 140   No. f-s = 163   No. J-s = 279   No. LU-s = 42
- No. nonlinear iterations = 162
- No. nonlinear convergence failures = 0
- No. error test failures = 1
+    Total internal steps taken    =      143
+    Total rhs function calls      =      174
+    Total Jacobian function calls =        3
+    Total root function calls     =        0
+    Total LU function calls       =       44
+    Num error test failures       =        4
+    Num nonlinear solver iters    =      173
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/short_no_pre/short_no_pre.out.cmp b/tests/tests/short_no_pre/short_no_pre.out.cmp
index 6a069a36..4ac9fa12 100644
--- a/tests/tests/short_no_pre/short_no_pre.out.cmp
+++ b/tests/tests/short_no_pre/short_no_pre.out.cmp
@@ -132,49 +132,22 @@ AtChem2 v1.3-dev
  Checking for constrained environment variables...
  Finished checking for constrained environment variables.
 
-
--------------
- Constraints
--------------
- Counting the variable-concentration species to be constrained (in file speciesConstrained.config)...
- Finished counting the names of variable-concentration constrained species.
- Number of names of variable-concentration constrained species: 0
- Counting the fixed-concentration species to be constrained (in file speciesConstant.config)...
- Finished counting the names of fixed-concentration constrained species.
- Number of names of fixed-concentration constrained species: 0
- Setting size of constraint arrays, n = 0
- Skipped reading the names of variable-concentration constrained species
- Reading concentration data for variable-concentration constrained species...
- Reading in the names and concentration of the fixed constrained species (in file speciesConstant.config)...
- Finished reading in the names and concentration of fixed-concentration species.
- Finished reading constrained species.
- Initialising concentrations of constrained species...
- Finished initialising concentrations of constrained species.
-
----------------
- Problem stats
----------------
-                        neq = 104
- numberOfConstrainedSpecies = 0
-                         t0 =       4.320E+04
-
- setting maxnumsteps ier = 0
- setting maxstep ier = 0
-
------------
- Model run
------------
- time = 43500
- time = 43800
- time = 44100
-
-------------------
- Final statistics
-------------------
- No. steps = 10   No. f-s = 14   No. J-s = 0   No. LU-s = 0
- No. nonlinear iterations = 10
- No. nonlinear convergence failures = 0
- No. error test failures = 0
-
- Runtime = 0
- Deallocating memory.
+Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
+
+Backtrace for this error:
+#0  0x7f8e254f22ed in ???
+#1  0x7f8e254f1503 in ???
+#2  0x7f8e25124f0f in ???
+#3  0x7f8e24ecb770 in cvLsInitialize
+	at /home/robertos/Downloads/AtChem-lib/sundials-7.5.0/src/cvode/cvode_ls.c:1439
+#4  0x7f8e24ec3164 in cvInitialSetup
+	at /home/robertos/Downloads/AtChem-lib/sundials-7.5.0/src/cvode/cvode.c:2078
+#5  0x7f8e24ec3164 in CVode
+	at /home/robertos/Downloads/AtChem-lib/sundials-7.5.0/src/cvode/cvode.c:1125
+#6  0x7f8e25ef3976 in __fcvode_mod_MOD_fcvode
+	at /home/robertos/Downloads/AtChem-lib/sundials-7.5.0/src/cvode/fmod_int64/fcvode_mod.f90:2105
+#7  0x55e89a6c25e8 in ???
+#8  0x55e89a69b4ce in ???
+#9  0x7f8e25107c86 in ???
+#10  0x55e89a69b4f9 in ???
+#11  0xffffffffffffffff in ???
diff --git a/tests/tests/spec_no_env_yes1/spec_no_env_yes1.out.cmp b/tests/tests/spec_no_env_yes1/spec_no_env_yes1.out.cmp
index 45bf4e31..ed4a3597 100644
--- a/tests/tests/spec_no_env_yes1/spec_no_env_yes1.out.cmp
+++ b/tests/tests/spec_no_env_yes1/spec_no_env_yes1.out.cmp
@@ -164,7 +164,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       0.000E+00
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -174,14 +173,18 @@ AtChem2 v1.3-dev
  time = 1800
  time = 2700
  time = 3600
-
 ------------------
  Final statistics
 ------------------
- No. steps = 151   No. f-s = 196   No. J-s = 269   No. LU-s = 37
- No. nonlinear iterations = 193
- No. nonlinear convergence failures = 0
- No. error test failures = 5
+    Total internal steps taken    =      159
+    Total rhs function calls      =      210
+    Total Jacobian function calls =        4
+    Total root function calls     =        0
+    Total LU function calls       =       37
+    Num error test failures       =        4
+    Num nonlinear solver iters    =      208
+    Num nonlinear solver fails    =        1
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/spec_no_env_yes2/spec_no_env_yes2.out.cmp b/tests/tests/spec_no_env_yes2/spec_no_env_yes2.out.cmp
index e4d31ffb..d1d9f590 100644
--- a/tests/tests/spec_no_env_yes2/spec_no_env_yes2.out.cmp
+++ b/tests/tests/spec_no_env_yes2/spec_no_env_yes2.out.cmp
@@ -166,21 +166,24 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 0
                          t0 =       0.000E+00
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
  Model run
 -----------
  time = 300
-
 ------------------
  Final statistics
 ------------------
- No. steps = 138   No. f-s = 178   No. J-s = 188   No. LU-s = 41
- No. nonlinear iterations = 177
- No. nonlinear convergence failures = 0
- No. error test failures = 7
+    Total internal steps taken    =      139
+    Total rhs function calls      =      174
+    Total Jacobian function calls =        3
+    Total root function calls     =        0
+    Total LU function calls       =       39
+    Num error test failures       =        4
+    Num nonlinear solver iters    =      173
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/spec_yes_env_no/spec_yes_env_no.out.cmp b/tests/tests/spec_yes_env_no/spec_yes_env_no.out.cmp
index e533e7e2..9b01ccba 100644
--- a/tests/tests/spec_yes_env_no/spec_yes_env_no.out.cmp
+++ b/tests/tests/spec_yes_env_no/spec_yes_env_no.out.cmp
@@ -167,7 +167,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 2
                          t0 =       3.600E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -175,14 +174,18 @@ AtChem2 v1.3-dev
 -----------
  time = 36200
  time = 36400
-
 ------------------
  Final statistics
 ------------------
- No. steps = 227   No. f-s = 273   No. J-s = 488   No. LU-s = 53
- No. nonlinear iterations = 272
- No. nonlinear convergence failures = 0
- No. error test failures = 3
+    Total internal steps taken    =      292
+    Total rhs function calls      =      375
+    Total Jacobian function calls =        6
+    Total root function calls     =        0
+    Total LU function calls       =       54
+    Num error test failures       =        5
+    Num nonlinear solver iters    =      374
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/spec_yes_env_no_with_jfac/spec_yes_env_no_with_jfac.out.cmp b/tests/tests/spec_yes_env_no_with_jfac/spec_yes_env_no_with_jfac.out.cmp
index 16b035e4..653dcbda 100644
--- a/tests/tests/spec_yes_env_no_with_jfac/spec_yes_env_no_with_jfac.out.cmp
+++ b/tests/tests/spec_yes_env_no_with_jfac/spec_yes_env_no_with_jfac.out.cmp
@@ -178,7 +178,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 2
                          t0 =       3.600E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -187,14 +186,18 @@ AtChem2 v1.3-dev
  basePhotoRate: J4
  time = 36200
  time = 36400
-
 ------------------
  Final statistics
 ------------------
- No. steps = 240   No. f-s = 279   No. J-s = 487   No. LU-s = 52
- No. nonlinear iterations = 278
- No. nonlinear convergence failures = 0
- No. error test failures = 1
+    Total internal steps taken    =      311
+    Total rhs function calls      =      411
+    Total Jacobian function calls =        6
+    Total root function calls     =        0
+    Total LU function calls       =       58
+    Num error test failures       =        6
+    Num nonlinear solver iters    =      410
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/spec_yes_env_no_with_jfac_fail1/spec_yes_env_no_with_jfac_fail1.out.cmp b/tests/tests/spec_yes_env_no_with_jfac_fail1/spec_yes_env_no_with_jfac_fail1.out.cmp
index b701cf9c..ffe2be50 100644
--- a/tests/tests/spec_yes_env_no_with_jfac_fail1/spec_yes_env_no_with_jfac_fail1.out.cmp
+++ b/tests/tests/spec_yes_env_no_with_jfac_fail1/spec_yes_env_no_with_jfac_fail1.out.cmp
@@ -169,7 +169,6 @@ STOP 2
  numberOfConstrainedSpecies = 2
                          t0 =       3.600E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
diff --git a/tests/tests/spec_yes_env_no_with_jfac_fixed/spec_yes_env_no_with_jfac_fixed.out.cmp b/tests/tests/spec_yes_env_no_with_jfac_fixed/spec_yes_env_no_with_jfac_fixed.out.cmp
index bb1079c1..02ccfc48 100644
--- a/tests/tests/spec_yes_env_no_with_jfac_fixed/spec_yes_env_no_with_jfac_fixed.out.cmp
+++ b/tests/tests/spec_yes_env_no_with_jfac_fixed/spec_yes_env_no_with_jfac_fixed.out.cmp
@@ -169,7 +169,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 2
                          t0 =       3.600E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -177,14 +176,18 @@ AtChem2 v1.3-dev
 -----------
  time = 36200
  time = 36400
-
 ------------------
  Final statistics
 ------------------
- No. steps = 226   No. f-s = 269   No. J-s = 470   No. LU-s = 50
- No. nonlinear iterations = 268
- No. nonlinear convergence failures = 0
- No. error test failures = 3
+    Total internal steps taken    =      267
+    Total rhs function calls      =      335
+    Total Jacobian function calls =        5
+    Total root function calls     =        0
+    Total LU function calls       =       51
+    Num error test failures       =        4
+    Num nonlinear solver iters    =      334
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/spec_yes_env_no_with_photo/spec_yes_env_no_with_photo.out.cmp b/tests/tests/spec_yes_env_no_with_photo/spec_yes_env_no_with_photo.out.cmp
index 6a19441b..8fe5ed47 100644
--- a/tests/tests/spec_yes_env_no_with_photo/spec_yes_env_no_with_photo.out.cmp
+++ b/tests/tests/spec_yes_env_no_with_photo/spec_yes_env_no_with_photo.out.cmp
@@ -177,7 +177,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 2
                          t0 =       3.600E+04
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -185,14 +184,18 @@ AtChem2 v1.3-dev
 -----------
  time = 36200
  time = 36400
-
 ------------------
  Final statistics
 ------------------
- No. steps = 224   No. f-s = 263   No. J-s = 484   No. LU-s = 49
- No. nonlinear iterations = 262
- No. nonlinear convergence failures = 0
- No. error test failures = 1
+    Total internal steps taken    =      297
+    Total rhs function calls      =      376
+    Total Jacobian function calls =        6
+    Total root function calls     =        0
+    Total LU function calls       =       56
+    Num error test failures       =        3
+    Num nonlinear solver iters    =      375
+    Num nonlinear solver fails    =        0
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/spec_yes_env_yes/spec_yes_env_yes.out.cmp b/tests/tests/spec_yes_env_yes/spec_yes_env_yes.out.cmp
index 67f07c07..74d75f21 100644
--- a/tests/tests/spec_yes_env_yes/spec_yes_env_yes.out.cmp
+++ b/tests/tests/spec_yes_env_yes/spec_yes_env_yes.out.cmp
@@ -170,21 +170,24 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 2
                          t0 =       1.200E+03
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
  Model run
 -----------
  time = 1500
-
 ------------------
  Final statistics
 ------------------
- No. steps = 174   No. f-s = 212   No. J-s = 237   No. LU-s = 42
- No. nonlinear iterations = 211
- No. nonlinear convergence failures = 0
- No. error test failures = 6
+    Total internal steps taken    =      169
+    Total rhs function calls      =      217
+    Total Jacobian function calls =        4
+    Total root function calls     =        0
+    Total LU function calls       =       42
+    Num error test failures       =        6
+    Num nonlinear solver iters    =      216
+    Num nonlinear solver fails    =        1
+  
 
  Runtime = 0
  Deallocating memory.
diff --git a/tests/tests/spec_yes_plus_fixed_env_no/spec_yes_plus_fixed_env_no.out.cmp b/tests/tests/spec_yes_plus_fixed_env_no/spec_yes_plus_fixed_env_no.out.cmp
index adcecce3..b400cb0b 100644
--- a/tests/tests/spec_yes_plus_fixed_env_no/spec_yes_plus_fixed_env_no.out.cmp
+++ b/tests/tests/spec_yes_plus_fixed_env_no/spec_yes_plus_fixed_env_no.out.cmp
@@ -169,7 +169,6 @@ AtChem2 v1.3-dev
  numberOfConstrainedSpecies = 4
                          t0 =       0.000E+00
 
- setting maxnumsteps ier = 0
  setting maxstep ier = 0
 
 -----------
@@ -178,14 +177,18 @@ AtChem2 v1.3-dev
  time = 300
  time = 600
  time = 900
-
 ------------------
  Final statistics
 ------------------
- No. steps = 153   No. f-s = 185   No. J-s = 249   No. LU-s = 43
- No. nonlinear iterations = 184
- No. nonlinear convergence failures = 0
- No. error test failures = 6
-
- Runtime = 0
+    Total internal steps taken    =      151
+    Total rhs function calls      =      186
+    Total Jacobian function calls =        3
+    Total root function calls     =        0
+    Total LU function calls       =       43
+    Num error test failures       =        6
+    Num nonlinear solver iters    =      185
+    Num nonlinear solver fails    =        0
+  
+
+ Runtime = 1
  Deallocating memory.
diff --git a/tools/install/Makefile.skel b/tools/install/Makefile.skel
index 0a145eac..5a02145a 100644
--- a/tools/install/Makefile.skel
+++ b/tools/install/Makefile.skel
@@ -3,7 +3,8 @@
 # Copyright (c) 2009-2012 Chris Martin, Kasia Boronska, Jenny Young,
 # Peter Jimack, Mike Pilling
 #
-# Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva
+# Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva, James Allsopp,
+# Neil Butcher
 #
 # This file is part of the AtChem2 software package.
 #
@@ -22,6 +23,7 @@ FORTC = "gnu"
 # Set the dependencies paths. Use full paths, not relative paths.
 # For example: $(HOME)/path/to/dependencies/directory/cvode/lib
 CVODELIBDIR = cvode/lib
+CVODEOBJDIR = cvode/fortran
 OPENLIBMDIR = openlibm
 FRUITDIR    = fruit_3.4.3
 
@@ -83,8 +85,17 @@ endif
 # set the CVODE and openlibm compilation flags
 LDFLAGS = -L$(CVODELIBDIR) -L$(OPENLIBMDIR) \
           -Wl,$(RPATH_OPTION),/usr/lib/:$(CVODELIBDIR):$(OPENLIBMDIR) \
-          -lopenlibm -lsundials_fcvode -lsundials_cvode -lsundials_fnvecserial \
-          -lsundials_nvecserial -ldl
+          -lopenlibm \
+          -lsundials_fcvode_mod \
+          -lsundials_fnvecserial_mod \
+          -lsundials_fcore_mod \
+          -lsundials_fsunlinsolspgmr_mod \
+          -lsundials_fsunmatrixband_mod \
+          -lsundials_fsunlinsoldense_mod \
+          -lsundials_fsunmatrixdense_mod \
+          -lsundials_cvode \
+          -lsundials_nvecserial \
+          -ldl
 
 # ----------------------------------------------------
 # EXECUTABLE setup
@@ -108,7 +119,7 @@ SRCS = $(CORE_SRCS) $(SRC)/atchem2.f90
 # the executable is rebuilt every time one of the fortran files (in $SRCS)
 # is modified
 $(AOUT): $(SRCS)
-	$(FORT_COMP) -o $(AOUT) -J$(OBJ) -I$(OBJ) $(SRCS) $(FFLAGS) $(LDFLAGS)
+	$(FORT_COMP) -o $(AOUT) -J$(OBJ) -I$(OBJ) -I$(CVODEOBJDIR) $(SRCS) $(FFLAGS) $(LDFLAGS)
 
 # secondary makefile for the Testsuite
 include tests/makefile.tests
diff --git a/tools/install/install_cvode.sh b/tools/install/install_cvode.sh
index 3a9cec86..777cbd13 100755
--- a/tools/install/install_cvode.sh
+++ b/tools/install/install_cvode.sh
@@ -1,12 +1,12 @@
 #!/bin/sh
 # -----------------------------------------------------------------------------
 #
-# Copyright (c) 2017 Sam Cox, Roberto Sommariva
+# Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva, Neil Butcher
 #
 # This file is part of the AtChem2 software package.
 #
-# This file is covered by the MIT license which can be found in the file
-# LICENSE.md at the top level of the AtChem2 distribution.
+# This file is licensed under the MIT license, which can be found in the file
+# `LICENSE` at the top level of the AtChem2 distribution.
 #
 # -----------------------------------------------------------------------------
 
@@ -26,7 +26,7 @@
 #   ./install_cvode.sh ~/path/to/dependencies/directory /path/to/fortran/compiler
 # -----------------------------------------------------------------------------
 
-SUNDIALS_VERSION="2.7.0"
+SUNDIALS_VERSION="7.5.0"
 
 # path to dependencies directory
 if [ -z "$1" ] ; then
@@ -89,6 +89,7 @@ cmake -DCMAKE_INSTALL_PREFIX="$DEP_DIR/cvode" \
       -DFCMIX_ENABLE:BOOL=ON \
       -DEXAMPLES_ENABLE:BOOL=OFF \
       -DCMAKE_MACOSX_RPATH:BOOL=ON \
+      -DBUILD_FORTRAN_MODULE_INTERFACE:BOOL=ON \
       ..
 if [ $? -ne 0 ] ; then
     printf "\n[cvode] cmake --> FAIL\n"
diff --git a/tools/install/install_openlibm.sh b/tools/install/install_openlibm.sh
index af4f7332..b75bce54 100755
--- a/tools/install/install_openlibm.sh
+++ b/tools/install/install_openlibm.sh
@@ -1,12 +1,12 @@
 #!/bin/sh
 # -----------------------------------------------------------------------------
 #
-# Copyright (c) 2017 Sam Cox, Roberto Sommariva
+# Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva
 #
 # This file is part of the AtChem2 software package.
 #
-# This file is covered by the MIT license which can be found in the file
-# LICENSE.md at the top level of the AtChem2 distribution.
+# This file is licensed under the MIT license, which can be found in the file
+# `LICENSE` at the top level of the AtChem2 distribution.
 #
 # -----------------------------------------------------------------------------
 
@@ -21,7 +21,7 @@
 #   ./install_openlibm.sh ~/path/to/dependencies/directory
 # -----------------------------------------------------------------------------
 
-OPENLIBM_VERSION="0.8.6"
+OPENLIBM_VERSION="0.8.7"
 
 # path to dependencies directory
 if [ -z "$1" ] ; then