diff --git a/process/caller.py b/process/caller.py
index cb8469ca41..6f04225dee 100644
--- a/process/caller.py
+++ b/process/caller.py
@@ -1,10 +1,22 @@
+from __future__ import annotations
 from process import fortran as ft
+import numpy as np
+import logging
+from process.final import finalise
+from process.io.mfile import MFile
+from process.utilities.f2py_string_patch import f2py_compatible_to_string
+from typing import Union, Tuple, TYPE_CHECKING
+
+if TYPE_CHECKING:
+    from process.main import Models
+
+logger = logging.getLogger(__name__)
 
 
 class Caller:
     """Calls physics and engineering models."""
 
-    def __init__(self, models, x):
+    def __init__(self, models: Models) -> None:
         """Initialise all physics and engineering models.
 
         To ensure that, at the start of a run, all physics/engineering
@@ -12,14 +24,155 @@ def __init__(self, models, x):
         called with the initial optimisation paramters, x.
 
         :param models: physics and engineering model objects
-        :type models: process.main.Models
-        :param x: optimisation parameters
-        :type x: np.ndarray
+        :type models: Models
         """
         self.models = models
-        self.call_models(x)
 
-    def call_models(self, xc):
+    @staticmethod
+    def check_agreement(
+        previous: Union[float, np.ndarray], current: Union[float, np.ndarray]
+    ) -> bool:
+        """Compare previous and current arrays for agreement within a tolerance.
+
+        :param previous: value(s) from previous models evaluation
+        :type previous: float | np.ndarray
+        :param current: value(s) from current models evaluation
+        :type current: float | np.ndarray
+        :return: whether values agree or not
+        :rtype: bool
+        """
+        return np.allclose(previous, current, rtol=1.0e-6)
+
+    def call_models(self, xc: np.ndarray, m: int) -> Tuple[float, np.ndarray]:
+        """Evalutate models until results are idempotent.
+
+        Ensure objective function and constraints are idempotent before returning.
+
+        :param xc: optimisation parameters
+        :type xc: np.ndarray
+        :param m: number of constraints
+        :type m: int
+        :raises RuntimeError: if values are non-idempotent after successive
+        evaluations
+        :return: objective function and constraints
+        :rtype: Tuple[float, np.ndarray]
+        """
+        objf_prev = None
+        conf_prev = None
+
+        # Evaluate models up to 10 times; any more implies non-converging values
+        for _ in range(10):
+            self._call_models_once(xc)
+            # Evaluate objective function and constraints
+            objf = ft.function_evaluator.funfom()
+            conf, _, _, _, _ = ft.constraints.constraint_eqns(m, -1)
+
+            if objf_prev is None and conf_prev is None:
+                # First run: run again to check idempotence
+                logger.debug("New optimisation parameter vector being evaluated")
+                objf_prev = objf
+                conf_prev = conf
+                continue
+
+            # Check for idempotence
+            if self.check_agreement(objf_prev, objf) and self.check_agreement(
+                conf_prev, conf
+            ):
+                # Idempotent: no longer changing, so return
+                logger.debug(
+                    "Model evaluations idempotent, returning objective "
+                    "function and constraints"
+                )
+                return objf, conf
+
+            # Not idempotent: still changing, so evaluate models again
+            logger.debug("Model evaluations not idempotent: evaluating again")
+            objf_prev = objf
+            conf_prev = conf
+
+        raise RuntimeError(
+            "After 10 model evaluations at the current optimisation parameter "
+            "vector, values for the objective function and constraints haven't "
+            "converged (don't produce idempotent values)."
+        )
+
+    def call_models_and_write_output(self, xc: np.ndarray, ifail: int) -> None:
+        """Evaluate models until results are idempotent, then write output files.
+
+        Ensure all outputs in mfile are idempotent before returning, by
+        evaluating models multiple times. Typically used at the end of an
+        optimisation, or in a non-optimising evaluation. Writes OUT.DAT and
+        MFILE.DAT with final results.
+
+        :param xc: optimisation parameter
+        :type xc: np.ndarray
+        :param ifail: return code of solver
+        :type ifail: int
+        :raises RuntimeError: if values are non-idempotent after successive
+        evaluations
+        """
+        # TODO The only way to ensure idempotence in all outputs is by comparing
+        # mfiles at this stage
+        previous_mfile_arr = None
+
+        try:
+            # Evaluate models up to 10 times; any more implies non-converging values
+            for _ in range(10):
+                # Divert OUT.DAT and MFILE.DAT output to scratch files for
+                # idempotence checking
+                ft.init_module.open_idempotence_files()
+                self._call_models_once(xc)
+                # Write mfile
+                finalise(self.models, ifail)
+
+                # Extract data from intermediate idempotence-checking mfile
+                mfile_path = (
+                    f2py_compatible_to_string(ft.global_variables.output_prefix)
+                    + "IDEM_MFILE.DAT"
+                )
+                mfile = MFile(mfile_path)
+                mfile_data = {}
+                for var in mfile.data.keys():
+                    mfile_data[var] = mfile.data[var].get_scan(-1)
+
+                # Extract floats from mfile dict into array for straightforward
+                # comparison: only compare floats
+                current_mfile_arr = np.array(
+                    [val for val in mfile_data.values() if isinstance(val, float)]
+                )
+                if previous_mfile_arr is None:
+                    # First run: need another run to compare with
+                    logger.debug(
+                        "New mfile created: evaluating models again to check idempotence"
+                    )
+                    previous_mfile_arr = np.copy(current_mfile_arr)
+                    continue
+
+                if self.check_agreement(previous_mfile_arr, current_mfile_arr):
+                    # Previous and current mfiles agree (idempotent)
+                    logger.debug("Mfiles idempotent, returning")
+                    # Divert OUT.DAT and MFILE.DAT output back to original files
+                    # now idempotence checking complete
+                    ft.init_module.close_idempotence_files()
+                    # Write final output file and mfile
+                    finalise(self.models, ifail)
+                    return
+
+                # Mfiles not yet idempotent: re-evaluate models
+                logger.debug("Mfiles not idempotent, evaluating models again")
+                previous_mfile_arr = np.copy(current_mfile_arr)
+
+            raise RuntimeError(
+                "Model evaluations at the current optimisation parameter vector "
+                "don't produce idempotent values in the final output."
+            )
+        except Exception:
+            # If exception in model evaluations or idempotence can't be
+            # achieved, delete intermediate idempotence files to clean up
+            ft.init_module.close_idempotence_files()
+            raise
+
+    def _call_models_once(self, xc: np.ndarray) -> None:
         """Call the physics and engineering models.
 
         This method is the principal caller of all the physics and
@@ -155,3 +308,21 @@ def call_models(self, xc):
         self.models.costs.run()
 
         # FISPACT and LOCA model (not used)- removed
+
+
+def write_output_files(models: Models, ifail: int) -> None:
+    """Evaluate models and write output files (OUT.DAT and MFILE.DAT).
+
+    :param models: physics and engineering models
+    :type models: Models
+    :param ifail: solver return code
+    :type ifail: int
+    """
+    n = ft.numerics.nvar
+    x = ft.numerics.xcm[:n]
+    # Call models, ensuring output mfiles are fully idempotent
+    caller = Caller(models)
+    caller.call_models_and_write_output(
+        xc=x,
+        ifail=ifail,
+    )
diff --git a/process/evaluators.py b/process/evaluators.py
index a4d2db18b5..3d3ad3f9fe 100644
--- a/process/evaluators.py
+++ b/process/evaluators.py
@@ -1,12 +1,9 @@
 from process.caller import Caller
 from process.fortran import global_variables as gv
-from process.fortran import constraints
 from process.fortran import cost_variables as cv
 from process.fortran import numerics
 from process.fortran import physics_variables as pv
-from process.fortran import stellarator_variables as sv
 from process.fortran import times_variables as tv
-from process.fortran import function_evaluator
 import numpy as np
 import math
 import logging
@@ -25,7 +22,7 @@ def __init__(self, models, x):
         :param x: optimisation parameters
         :type x: np.ndarray
         """
-        self.caller = Caller(models, x)
+        self.caller = Caller(models)
 
     def fcnvmc1(self, n, m, xv, ifail):
         """Function evaluator for VMCON.
@@ -53,29 +50,7 @@ def fcnvmc1(self, n, m, xv, ifail):
         conf = np.zeros(m, dtype=np.float64, order="F")
 
         # Evaluate machine parameters at xv
-        self.caller.call_models(xv)
-
-        # Convergence loop to ensure burn time consistency
-        if sv.istell == 0:
-            for _ in range(10):
-                if abs((tv.tburn - tv.tburn0) / max(tv.tburn, 0.01)) <= 0.001:
-                    break
-
-                self.caller.call_models(xv)
-                if gv.verbose == 1:
-                    print("Internal tburn consistency check: ", tv.tburn, tv.tburn0)
-            else:
-                print(
-                    "Burn time values are not consistent in iteration: ",
-                    numerics.nviter,
-                )
-                print("tburn, tburn0: ", tv.tburn, tv.tburn0)
-
-        # Evaluate figure of merit (objective function)
-        objf = function_evaluator.funfom()
-
-        # Evaluate constraint equations
-        conf, _, _, _, _ = constraints.constraint_eqns(m, -1)
+        objf, conf = self.caller.call_models(xv, m)
 
         # Verbose diagnostics
         if gv.verbose == 1:
@@ -141,14 +116,10 @@ def fcnvmc2(self, n, m, xv, lcnorm):
                     xbac[i] = xv[j] * (1.0 - numerics.epsfcn)
 
             # Evaluate at (x+dx)
-            self.caller.call_models(xfor)
-            ffor = function_evaluator.funfom()
-            cfor, _, _, _, _ = constraints.constraint_eqns(m, -1)
+            ffor, cfor = self.caller.call_models(xfor, m)
 
             # Evaluate at (x-dx)
-            self.caller.call_models(xbac)
-            fbac = function_evaluator.funfom()
-            cbac, _, _, _, _ = constraints.constraint_eqns(m, -1)
+            fbac, cbac = self.caller.call_models(xbac, m)
 
             # Calculate finite difference gradients
             fgrd[i] = (ffor - fbac) / (xfor[i] - xbac[i])
@@ -162,6 +133,6 @@ def fcnvmc2(self, n, m, xv, lcnorm):
         # variable in the solution vector is inconsistent with its value
         # shown elsewhere in the output file, which is a factor (1-epsfcn)
         # smaller (i.e. its xbac value above).
-        self.caller.call_models(xv)
+        self.caller.call_models(xv, m)
 
         return fgrd, cnorm
diff --git a/process/final.py b/process/final.py
index 94a61e6502..06a31142fb 100644
--- a/process/final.py
+++ b/process/final.py
@@ -1,4 +1,5 @@
 """Final output at the end of a scan."""
+
 from process import fortran as ft
 from process.fortran import final_module as fm
 from process import output as op
@@ -7,6 +8,8 @@
 def finalise(models, ifail):
     """Routine to print out the final point in the scan.
 
+    Writes to OUT.DAT and MFILE.DAT.
+
     :param models: physics and engineering model objects
     :type models: process.main.Models
     :param ifail: error flag
@@ -18,7 +21,7 @@ def finalise(models, ifail):
     if ft.numerics.ioptimz == -2:
         fm.no_optimisation()
 
-    # Write output to OUT.DAT
+    # Write output to OUT.DAT and MFILE.DAT
     op.write(models, ft.constants.nout)
 
     fm.final_output()
diff --git a/process/main.py b/process/main.py
index 3d6859988d..ef91d5e740 100644
--- a/process/main.py
+++ b/process/main.py
@@ -49,7 +49,6 @@
 from process.plasma_geometry import PlasmaGeom
 from process.pulse import Pulse
 from process.scan import Scan
-from process import final
 from process.stellarator import Stellarator
 from process.structure import Structure
 from process.build import Build
@@ -61,7 +60,6 @@
 from process.availability import Availability
 from process.ife import IFE
 from process.costs_2015 import Costs2015
-from process.caller import Caller
 from process.power import Power
 from process.cs_fatigue import CsFatigue
 from process.physics import Physics
@@ -73,6 +71,8 @@
 from process.fw import Fw
 from process.current_drive import CurrentDrive
 from process.impurity_radiation import initialise_imprad
+from process.caller import write_output_files
+
 
 from pathlib import Path
 import os
@@ -459,25 +459,21 @@ def run_scan(self, solver):
         :param solver: which solver to use, as specified in solver.py
         :type solver: str
         """
-        if fortran.numerics.ioptimz >= 0:
+        if fortran.numerics.ioptimz == 1:
+            # VMCON optimisation
             self.scan = Scan(self.models, solver)
-        else:
-            # If no optimisation will be done, compute the OP variables now
-            if fortran.numerics.ioptimz == -2:
-                # Get optimisation parameters x, perform first evaluation of models
-                fortran.define_iteration_variables.loadxc()
-                n = fortran.numerics.nvar
-                x = fortran.numerics.xcm[:n]
-                caller = Caller(self.models, x)
-
-                # To ensure that, at the start of a run, all physics/engineering
-                # variables are fully initialised with consistent values, we perform
-                # a second evaluation call here
-                caller.call_models(x)
-                self.ifail = 6
-
-            final.finalise(self.models, self.ifail)
+        elif fortran.numerics.ioptimz == -2:
+            # No optimisation: compute the output variables now
+            # Get optimisation parameters x, evaluate models
+            fortran.define_iteration_variables.loadxc()
+            self.ifail = 6
+            write_output_files(models=self.models, ifail=self.ifail)
             self.show_errors()
+        else:
+            raise ValueError(
+                f"Invalid ioptimz value: {fortran.numerics.ioptimz}. Please "
+                "select either 1 (optimise) or -2 (no optimisation)."
+            )
 
     def show_errors(self):
         """Report all informational/error messages encountered."""
diff --git a/process/scan.py b/process/scan.py
index 664c4eff48..56689fd737 100644
--- a/process/scan.py
+++ b/process/scan.py
@@ -2,8 +2,8 @@
 from process.fortran import scan_module
 from process.fortran import numerics
 from process.optimiser import Optimiser
-from process import final
 import numpy as np
+from process.caller import write_output_files
 
 
 class Scan:
@@ -34,7 +34,7 @@ def run_scan(self):
 
         if scan_module.isweep == 0:
             ifail = self.doopt()
-            final.finalise(self.models, ifail)
+            write_output_files(models=self.models, ifail=ifail)
             error_handling.show_errors()
             return
 
@@ -76,7 +76,7 @@ def scan_1d(self):
             scan_module.scan_1d_write_point_header(iscan)
             ifail = self.doopt()
             scan_1d_ifail_dict[iscan] = ifail
-            final.finalise(self.models, ifail)
+            write_output_files(models=self.models, ifail=ifail)
 
             # outvar is an intent(out) of scan_1d_store_output()
             outvar = scan_module.scan_1d_store_output(
@@ -146,7 +146,7 @@ def scan_2d(self):
                 )
                 ifail = self.doopt()
 
-                final.finalise(self.models, ifail)
+                write_output_files(models=self.models, ifail=ifail)
 
                 outvar, sweep_1_vals, sweep_2_vals = scan_module.scan_2d_store_output(
                     ifail,
diff --git a/source/fortran/final_module.f90 b/source/fortran/final_module.f90
index b67fc728e4..06ae0c2680 100644
--- a/source/fortran/final_module.f90
+++ b/source/fortran/final_module.f90
@@ -73,7 +73,7 @@ subroutine no_optimisation()
   end subroutine no_optimisation
 
   subroutine final_output()
-    use constants, only: iotty
+    use constants, only: iotty, mfile
     use numerics, only: nfev1, ncalls, xcm, ioptimz, nviter, &
       nvar
     use define_iteration_variables, only: loadxc
@@ -117,5 +117,9 @@ subroutine final_output()
           t2,'The HYBRD point required ',i5,' iterations',/, &
           t2,'The optimisation required ',i5,' iterations',/, &
           t2,'There were ',i6,' function calls')
+
+    ! Required to ensure mfile fully written before any model evaluation
+    ! idempotence checks
+    flush(mfile)
   end subroutine final_output
 end module final_module
diff --git a/source/fortran/init_module.f90 b/source/fortran/init_module.f90
index 008332c9ca..18d80364be 100644
--- a/source/fortran/init_module.f90
+++ b/source/fortran/init_module.f90
@@ -168,6 +168,39 @@ subroutine init
 
    end subroutine init
 
+   subroutine open_idempotence_files
+      ! Open new output file and mfile to write output to
+      ! This is used when checking model evaluation idempotence, to avoid
+      ! polluting the final output file and mfile with intermediate result checks
+      use global_variables, only: output_prefix
+      use constants, only: nout, mfile
+      implicit none
+
+      ! Close existing output file and mfile (could be original out and mfiles
+      ! or idem scratch files)
+      close(unit = nout)
+      close(unit = mfile)
+      ! Open scratch files with same units
+      open(unit=nout, file=trim(output_prefix)//'IDEM_OUT.DAT', action='write', status='replace')
+      open(unit=mfile, file=trim(output_prefix)//'IDEM_MFILE.DAT', action='write', status='replace')
+   end subroutine open_idempotence_files
+
+   subroutine close_idempotence_files
+      ! Close the intermediate idempotence-check files, deleting them in the process
+      ! Re-open the original OUT.DAT and MFILE.DAT output files, ready to write
+      ! the final data, now model evaluation idempotence has been checked
+      use global_variables, only: output_prefix
+      use constants, only: nout, mfile
+      implicit none
+
+      ! Close idempotence files, deleting them in the process
+      close(unit = nout, status="delete")
+      close(unit = mfile, status="delete")
+      ! Re-open original output file and mfile, appending future output to them
+      open(unit=nout, file=trim(output_prefix)//'OUT.DAT', action='write', position='append')
+      open(unit=mfile, file=trim(output_prefix)//'MFILE.DAT', action='write', position='append')
+   end subroutine close_idempotence_files
+
    subroutine finish
       ! Originally at the end of the "program", this subroutine writes some final
       ! lines via the output module and then closes any open files. This is