diff --git a/examples/single_model_evaluation.ex.py b/examples/single_model_evaluation.ex.py
index 80f17024e7..4f67a420b2 100644
--- a/examples/single_model_evaluation.ex.py
+++ b/examples/single_model_evaluation.ex.py
@@ -70,7 +70,7 @@ def print_values():
 process.data_structure.impurity_radiation_module.f_nd_impurity_electron_array[13] = (
     5.0e-5
 )
-single_run.models.physics.physics()
+single_run.models.physics.run()
 print_values()
 
 # %% [markdown]
@@ -97,7 +97,7 @@ def run_impurities(w_imp_fracs):
         process.data_structure.impurity_radiation_module.f_nd_impurity_electron_array[
             13
         ] = imp_frac
-        single_run.models.physics.physics()
+        single_run.models.physics.run()
 
         # Evaluate constraint equation 15 (L-H threshold constraint)
         con15_value = ConstraintManager.evaluate_constraint(15).normalised_residual
diff --git a/process/core/caller.py b/process/core/caller.py
index a16d55a3c7..317b353ad9 100644
--- a/process/core/caller.py
+++ b/process/core/caller.py
@@ -259,43 +259,43 @@ def _call_models_once(self, xc: np.ndarray):
         # Perform the various function calls
         # Stellarator caller
         if data_structure.stellarator_variables.istell != 0:
-            self.models.stellarator.run(output=False)
+            self.models.stellarator.run()
             # TODO Is this return safe?
             return
 
         # Inertial Fusion Energy calls
         if data_structure.ife_variables.ife != 0:
-            self.models.ife.run(output=False)
+            self.models.ife.run()
             return
 
         # Tokamak calls
         # Plasma geometry model
-        self.models.plasma_geom.plasma_geometry()
+        self.models.plasma_geom.run()
 
         # Machine Build Model
         # Radial build
         self.models.build.run()
 
-        self.models.physics.physics()
+        self.models.physics.run()
 
         # Toroidal field coil model
 
         # Toroidal field coil resistive model
         if data_structure.tfcoil_variables.i_tf_sup == 0:
-            self.models.copper_tf_coil.run(output=False)
+            self.models.copper_tf_coil.run()
 
         # Toroidal field coil superconductor model
         if data_structure.tfcoil_variables.i_tf_sup == 1:
-            self.models.sctfcoil.run(output=False)
+            self.models.sctfcoil.run()
 
         if data_structure.tfcoil_variables.i_tf_sup == 2:
-            self.models.aluminium_tf_coil.run(output=False)
+            self.models.aluminium_tf_coil.run()
 
         # Poloidal field and central solenoid model
         self.models.pfcoil.run()
 
         # Pulsed reactor model
-        self.models.pulse.run(output=False)
+        self.models.pulse.run()
 
         # First wall model
         self.models.fw.run()
@@ -316,52 +316,37 @@ def _call_models_once(self, xc: np.ndarray):
         """
         if data_structure.fwbs_variables.i_blanket_type == 1:
             # CCFE HCPB model
-            self.models.ccfe_hcpb.run(output=False)
+            self.models.ccfe_hcpb.run()
 
         elif data_structure.fwbs_variables.i_blanket_type == 5:
             # DCLL model
-            self.models.dcll.run(output=False)
+            self.models.dcll.run()
 
-        self.models.divertor.run(output=False)
+        self.models.divertor.run()
 
         self.models.cryostat.run()
 
         # Structure Model
-        self.models.structure.run(output=False)
+        self.models.structure.run()
 
         # Tight aspect ratio machine model
         if (
             data_structure.physics_variables.itart == 1
             and data_structure.tfcoil_variables.i_tf_sup != 1
         ):
-            self.models.tfcoil.cntrpst()
+            self.models.tfcoil.run()
 
-        # Toroidal field coil power model
-        self.models.power.tfpwr(output=False)
-
-        # Poloidal field coil power model
-        self.models.power.pfpwr(output=False)
-
-        # Plant heat transport part 1
-        self.models.power.component_thermal_powers()
-
-        # Cryoplant loads
-        self.models.power.calculate_cryo_loads()
+        # Power model
+        self.models.power.run()
 
         # Vacuum model
-        self.models.vacuum.run(output=False)
+        self.models.vacuum.run()
 
         # Buildings model
-        self.models.buildings.run(output=False)
-
-        # Plant AC power requirements
-        self.models.power.acpow(output=False)
-
-        # Plant heat transport pt 2 & 3
-        self.models.power.plant_electric_production()
+        self.models.buildings.run()
 
         # Availability model
-        self.models.availability.run(output=False)
+        self.models.availability.run()
 
         # Water usage in secondary cooling system
         self.models.water_use.run()
diff --git a/process/core/output.py b/process/core/output.py
index a69d2303f8..86421f4f23 100644
--- a/process/core/output.py
+++ b/process/core/output.py
@@ -22,12 +22,12 @@ def write(models, _outfile):
 
     # Call stellarator output routine instead if relevant
     if data_structure.stellarator_variables.istell != 0:
-        models.stellarator.run(output=True)
+        models.stellarator.output()
         return
 
     #  Call IFE output routine instead if relevant
     if data_structure.ife_variables.ife != 0:
-        models.ife.run(output=True)
+        models.ife.output()
         return
 
     # Costs model
@@ -37,70 +37,56 @@ def write(models, _outfile):
     # 0    |  1990 costs model
     # 1    |  2015 Kovari model
     # 2    |  Custom model
-    models.costs.run()
     models.costs.output()
 
     # Availability model
-    models.availability.run(output=True)
+    models.availability.output()
 
-    # Writing the output from physics.f90 into OUT.DAT + MFILE.DAT
-    models.physics.calculate_effective_charge_ionisation_profiles()
-    models.physics.outplas()
+    # Physics model
+    models.physics.output()
 
     # Detailed physics, currently only done at final point as values are not used
     # by any other functions
-    models.physics_detailed.run()
-    models.physics_detailed.output_detailed_physics()
+    models.physics_detailed.output()
 
-    # TODO what is this? Not in caller.f90?
-    models.current_drive.output_current_drive()
+    # TODO what is this? Not in caller.py?
+    models.current_drive.output()
 
     # Pulsed reactor model
-    models.pulse.run(output=True)
-    models.physics.outtim()
+    models.pulse.output()
 
-    models.divertor.run(output=True)
+    models.divertor.output()
 
     # Machine Build Model
-    # Radial build
-    models.build.calculate_radial_build(output=True)
-
-    # Vertical build
-    models.build.calculate_vertical_build(output=True)
+    models.build.output()
 
     # Cryostat build
-    models.cryostat.cryostat_output()
+    models.cryostat.output()
 
     # Toroidal field coil copper model
     if data_structure.tfcoil_variables.i_tf_sup == 0:
-        models.copper_tf_coil.run(output=True)
+        models.copper_tf_coil.output()
 
     # Toroidal field coil superconductor model
     if data_structure.tfcoil_variables.i_tf_sup == 1:
-        models.sctfcoil.run(output=True)
-        models.sctfcoil.output_tf_superconductor_info()
+        models.sctfcoil.output()
 
     # Toroidal field coil aluminium model
     if data_structure.tfcoil_variables.i_tf_sup == 2:
-        models.aluminium_tf_coil.run(output=True)
+        models.aluminium_tf_coil.output()
 
     # Tight aspect ratio machine model
     if (
         data_structure.physics_variables.itart == 1
         and data_structure.tfcoil_variables.i_tf_sup != 1
     ):
-        models.tfcoil.iprint = 1
-        models.tfcoil.cntrpst()
-        models.tfcoil.iprint = 0
+        models.tfcoil.output()
 
     # Poloidal field coil model
     models.pfcoil.output()
 
     # Structure Model
-    models.structure.run(output=True)
-
-    # Poloidal field coil inductance calculation
-    models.pfcoil.output_induct()
+    models.structure.output()
 
     # Blanket model
     # Blanket switch values
@@ -112,48 +98,30 @@ def write(models, _outfile):
     # 4    |  KIT HCLL model
     # 5    |  DCLL model
 
-    models.shield.output_shld_areas_and_volumes()
-    models.vacuum_vessel.output_vv_areas_and_volumes()
+    models.shield.output()
+    models.vacuum_vessel.output()
 
     # First wall geometry
-    models.fw.output_fw_geometry()
-
-    # First wall surface loads
-    models.fw.output_fw_surface_loads()
-
-    # First wall pumping
-    models.fw.output_fw_pumping()
+    models.fw.output()
 
     if data_structure.fwbs_variables.i_blanket_type == 1:
         # CCFE HCPB model
-        models.ccfe_hcpb.run(output=True)
+        models.ccfe_hcpb.output()
 
     elif data_structure.fwbs_variables.i_blanket_type == 5:
         # DCLL model
-        models.dcll.run(output=True)
+        models.dcll.output()
 
     # FISPACT and LOCA model (not used)- removed
 
-    # Toroidal field coil power model
-    models.power.tfpwr(output=True)
-
-    # Poloidal field coil power model !
-    models.power.pfpwr(output=True)
+    # Power model
+    models.power.output()
 
     # Vacuum model
-    models.vacuum.run(output=True)
+    models.vacuum.output()
 
     # Buildings model
-    models.buildings.run(output=True)
-
-    # Plant AC power requirements
-    models.power.acpow(output=True)
-
-    # Plant heat transport pt 2 & 3
-    models.power.output_cryogenics()
-    models.power.output_plant_thermal_powers()
-    models.power.output_plant_electric_powers()
-    models.power.output_power_profiles_over_time()
+    models.buildings.output()
 
     # Water usage in secondary cooling system
     models.water_use.output()
diff --git a/process/models/availability.py b/process/models/availability.py
index e920976c23..77bcdee860 100644
--- a/process/models/availability.py
+++ b/process/models/availability.py
@@ -6,6 +6,7 @@
 from process.core import constants
 from process.core import process_output as po
 from process.core.exceptions import ProcessValueError
+from process.core.model import Model
 from process.data_structure import constraint_variables as ctv
 from process.data_structure import cost_variables as cv
 from process.data_structure import divertor_variables as dv
@@ -35,7 +36,7 @@
 }
 
 
-class Availability:
+class Availability(Model):
     """Module containing plant availability routines
 
 
@@ -46,6 +47,9 @@ class Availability:
     def __init__(self):
         self.outfile = constants.NOUT  # output file unit
 
+    def output(self):
+        self.run(output=True)
+
     def run(self, output: bool = False):
         """Run appropriate availability model
 
diff --git a/process/models/blankets/blanket_library.py b/process/models/blankets/blanket_library.py
index 197bd434d9..9a2d24f976 100644
--- a/process/models/blankets/blanket_library.py
+++ b/process/models/blankets/blanket_library.py
@@ -8,6 +8,7 @@
 from process.core import process_output as po
 from process.core.coolprop_interface import FluidProperties
 from process.core.exceptions import ProcessValueError
+from process.core.model import Model
 from process.data_structure import (
     blanket_library,
     build_variables,
@@ -36,12 +37,18 @@
 # FCI         Flow Channel Insert
 
 
-class BlanketLibrary:
+class BlanketLibrary(Model):
     def __init__(self, fw):
         self.outfile = constants.NOUT
 
         self.fw = fw
 
+    def output(self):
+        """This model doesn't have any output"""
+
+    def run(self):
+        """This model doesn't need to be run"""
+
     def component_volumes(self):
         """Calculate the blanket, shield, vacuum vessel and cryostat volumes
 
diff --git a/process/models/blankets/dcll.py b/process/models/blankets/dcll.py
index aaf26d5f94..52f364e595 100644
--- a/process/models/blankets/dcll.py
+++ b/process/models/blankets/dcll.py
@@ -86,10 +86,13 @@ class DCLL(InboardBlanket, OutboardBlanket):
                      Design 167, 112380
 
 
-    Note: request for when CCFE Bluemira nutronics work is added: output maximum values, as well as average values, for wall neutronics calculation if possible.
+    Note: request for when CCFE Bluemira neutronics work is added: output maximum values, as well as average values, for wall neutronics calculation if possible.
     """
 
-    def run(self, output: bool):
+    def output(self):
+        self.run(output=True)
+
+    def run(self, output: bool = False):
         self.component_volumes()
         dia_blkt_channel = self.pipe_hydraulic_diameter(i_channel_shape=1)
         fwbs_variables.radius_blkt_channel = dia_blkt_channel / 2
diff --git a/process/models/blankets/hcpb.py b/process/models/blankets/hcpb.py
index 73bb970e0f..f8dd375c72 100644
--- a/process/models/blankets/hcpb.py
+++ b/process/models/blankets/hcpb.py
@@ -40,7 +40,10 @@ class CCFE_HCPB(OutboardBlanket, InboardBlanket):
         doi: https://doi.org/10.1016/j.fusengdes.2016.01.007.
     """
 
-    def run(self, output: bool):
+    def output(self):
+        self.run(output=True)
+
+    def run(self, output: bool = False):
         # Coolant type
         fwbs_variables.i_blkt_coolant_type = 1
         # Note that the first wall coolant is now input separately.
diff --git a/process/models/build.py b/process/models/build.py
index 9f88bb11fb..c4e0656f92 100644
--- a/process/models/build.py
+++ b/process/models/build.py
@@ -4,6 +4,7 @@
 
 from process.core import constants
 from process.core import process_output as po
+from process.core.model import Model
 from process.data_structure import (
     build_variables,
     buildings_variables,
@@ -21,11 +22,18 @@
 logger = logging.getLogger(__name__)
 
 
-class Build:
+class Build(Model):
     def __init__(self):
         self.outfile = constants.NOUT
         self.mfile = constants.MFILE
 
+    def output(self):
+        # Radial build
+        self.calculate_radial_build(output=True)
+
+        # Vertical build
+        self.calculate_vertical_build(output=True)
+
     def run(self):
         self.calculate_radial_build(output=False)
         self.calculate_vertical_build(output=False)
diff --git a/process/models/buildings.py b/process/models/buildings.py
index dde467b438..6e5a661532 100644
--- a/process/models/buildings.py
+++ b/process/models/buildings.py
@@ -4,6 +4,7 @@
 
 from process.core import constants
 from process.core import process_output as po
+from process.core.model import Model
 from process.data_structure import (
     build_variables,
     buildings_variables,
@@ -20,7 +21,7 @@
 logger = logging.getLogger(__name__)
 
 
-class Buildings:
+class Buildings(Model):
     """
 
     This module contains routines for calculating the
@@ -33,6 +34,9 @@ def __init__(self):
         """
         self.outfile = constants.NOUT  # output file unit
 
+    def output(self):
+        self.run(output=True)
+
     def run(self, output: bool = False):
         # Find TF coil radial positions
         # outboard edge: outboard mid-leg radial position + half-thickness of outboard leg
diff --git a/process/models/costs/costs.py b/process/models/costs/costs.py
index 258d18093f..a700da5633 100644
--- a/process/models/costs/costs.py
+++ b/process/models/costs/costs.py
@@ -94,6 +94,7 @@ def run(self):
             self.coelc()
 
     def output(self):
+        self.run()
         if cost_variables.output_costs == 0:
             return
 
diff --git a/process/models/costs/costs_2015.py b/process/models/costs/costs_2015.py
index f9dfb3b6a9..f832d933ce 100644
--- a/process/models/costs/costs_2015.py
+++ b/process/models/costs/costs_2015.py
@@ -4,6 +4,7 @@
 
 from process.core import constants
 from process.core import process_output as po
+from process.core.model import Model
 from process.data_structure import (
     build_variables,
     cost_2015_variables,
@@ -22,7 +23,7 @@
 logger = logging.getLogger(__name__)
 
 
-class Costs2015:
+class Costs2015(Model):
     def __init__(self):
         self.outfile = constants.NOUT
 
diff --git a/process/models/cryostat.py b/process/models/cryostat.py
index 3fd475b02a..edfdaa35a1 100644
--- a/process/models/cryostat.py
+++ b/process/models/cryostat.py
@@ -2,6 +2,7 @@
 
 from process.core import constants
 from process.core import process_output as po
+from process.core.model import Model
 from process.data_structure import (
     blanket_library,
     build_variables,
@@ -11,7 +12,7 @@
 )
 
 
-class Cryostat:
+class Cryostat(Model):
     def __init__(self):
         self.outfile = constants.NOUT
 
@@ -83,7 +84,7 @@ def external_cryo_geometry():
             fwbs_variables.vol_vv + fwbs_variables.vol_cryostat
         ) * fwbs_variables.den_steel
 
-    def cryostat_output(self):
+    def output(self):
         """Outputs the cryostat geometry details to the output file."""
         po.oheadr(self.outfile, "Cryostat build")
 
diff --git a/process/models/divertor.py b/process/models/divertor.py
index cf0a76d988..e044ac238a 100644
--- a/process/models/divertor.py
+++ b/process/models/divertor.py
@@ -5,6 +5,7 @@
 from process.core import constants
 from process.core import process_output as po
 from process.core.exceptions import ProcessValueError
+from process.core.model import Model
 from process.data_structure import build_variables as bv
 from process.data_structure import divertor_variables as dv
 from process.data_structure import fwbs_variables as fwbs
@@ -12,7 +13,7 @@
 from process.data_structure import tfcoil_variables as tfv
 
 
-class Divertor:
+class Divertor(Model):
     """Module containing divertor routines
 
     This module contains routines relevant for calculating the
@@ -22,7 +23,10 @@ class Divertor:
     def __init__(self):
         self.outfile = constants.NOUT  # output file unit
 
-    def run(self, output: bool):
+    def output(self):
+        self.run(output=True)
+
+    def run(self, output: bool = False):
         """Routine to call the divertor model
 
         This subroutine calls the divertor routine. This routine scales
diff --git a/process/models/fw.py b/process/models/fw.py
index d59342e24d..8a59dd6e06 100644
--- a/process/models/fw.py
+++ b/process/models/fw.py
@@ -6,6 +6,7 @@
 from process.core import process_output as po
 from process.core.coolprop_interface import FluidProperties
 from process.core.exceptions import ProcessValueError
+from process.core.model import Model
 from process.data_structure import (
     blanket_library,
     build_variables,
@@ -24,11 +25,21 @@
 logger = logging.getLogger(__name__)
 
 
-class FirstWall:
+class FirstWall(Model):
     def __init__(self):
         self.outfile = constants.NOUT
         self.blanket_library = BlanketLibrary(fw=self)
 
+    def output(self):
+        # First wall geometry
+        self.output_fw_geometry()
+
+        # First wall surface loads
+        self.output_fw_surface_loads()
+
+        # First wall pumping
+        self.output_fw_pumping()
+
     def run(self):
         fwbs_variables.dz_fw_half = self.calculate_first_wall_half_height(
             z_plasma_xpoint_lower=build_variables.z_plasma_xpoint_lower,
diff --git a/process/models/ife.py b/process/models/ife.py
index d87083a39b..825a921dfd 100644
--- a/process/models/ife.py
+++ b/process/models/ife.py
@@ -9,6 +9,7 @@
 
 from process.core import constants, process_output
 from process.core.exceptions import ProcessValueError
+from process.core.model import Model
 from process.data_structure import (
     buildings_variables,
     cost_variables,
@@ -34,7 +35,7 @@
 ]
 
 
-class IFE:
+class IFE(Model):
     """Module containing Inertial Fusion Energy device routines
 
     N/A
@@ -55,7 +56,10 @@ def __init__(self, availability, costs):
         self.availability = availability
         self.costs = costs
 
-    def run(self, output: bool):
+    def output(self):
+        self.run(output=True)
+
+    def run(self, output: bool = False):
         """Routine to output the physics and engineering information
         relevant to inertial fusion energy power plants
 
diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py
index eefcbd1df2..6ab4a15c1e 100644
--- a/process/models/pfcoil.py
+++ b/process/models/pfcoil.py
@@ -11,6 +11,7 @@
 from process.core import constants
 from process.core import process_output as op
 from process.core.exceptions import ProcessValueError
+from process.core.model import Model
 from process.data_structure import build_variables as bv
 from process.data_structure import constraint_variables as ctv
 from process.data_structure import cs_fatigue_variables as csfv
@@ -28,7 +29,7 @@
 logger = logging.getLogger(__name__)
 
 
-class PFCoil:
+class PFCoil(Model):
     """Calculate poloidal field coil system parameters."""
 
     def __init__(self, cs_fatigue):
@@ -54,6 +55,7 @@ def output(self):
         self.cs_coil.output_cs_structure()
         self.outpf()
         self.outvolt()
+        self.output_induct()
 
     def output_induct(self):
         """Output poloidal field coil inductance calculation."""
diff --git a/process/models/physics/bootstrap_current.py b/process/models/physics/bootstrap_current.py
index 6a100c11c0..aae2bb302f 100644
--- a/process/models/physics/bootstrap_current.py
+++ b/process/models/physics/bootstrap_current.py
@@ -9,6 +9,7 @@
 from process.core import constants
 from process.core import process_output as po
 from process.core.exceptions import ProcessValueError
+from process.core.model import Model
 from process.data_structure import (
     current_drive_variables,
     physics_variables,
@@ -37,7 +38,7 @@ class BootstrapCurrentFractionModel(IntEnum):
     SUGIYAMA_H_MODE = 13
 
 
-class PlasmaBootstrapCurrent:
+class PlasmaBootstrapCurrent(Model):
     """Class to hold plasma bootstrap current for plasma processing."""
 
     def __init__(self, plasma_profile: PlasmaProfile) -> None:
@@ -1195,7 +1196,7 @@ def bootstrap_fraction_sugiyama_h_mode(
             * temp_plasma_pedestal_kev**0.0552
         )
 
-    def output_bootstrap_current_information(self):
+    def output(self):
         """Output the calculated bootstrap current information to the output file."""
 
         po.ovarrf(
diff --git a/process/models/physics/current_drive.py b/process/models/physics/current_drive.py
index 765c55a566..ab8d18bce2 100644
--- a/process/models/physics/current_drive.py
+++ b/process/models/physics/current_drive.py
@@ -7,6 +7,7 @@
     process_output as po,
 )
 from process.core.exceptions import ProcessError, ProcessValueError
+from process.core.model import Model
 from process.data_structure import (
     current_drive_variables,
     heat_transport_variables,
@@ -1433,7 +1434,7 @@ def lower_hybrid_ehst(
         )
 
 
-class CurrentDrive:
+class CurrentDrive(Model):
     def __init__(
         self,
         plasma_profile: PlasmaProfile,
@@ -1452,6 +1453,9 @@ def __init__(
         self.neutral_beam = neutral_beam
         self.electron_bernstein = electron_bernstein
 
+    def run(self):
+        """This model doesn't need to be run"""
+
     def cudriv(self):
         """Calculate the current drive power requirements.
 
@@ -2154,7 +2158,7 @@ def calculate_dimensionless_current_drive_efficiency(
             * (c_hcd_driven / p_hcd_injected)
         )
 
-    def output_current_drive(self):
+    def output(self):
         """Output the current drive information to the output file.
         This method writes the current drive information to the output file.
         """
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index d313c02eca..34f987e425 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -11,6 +11,7 @@
 from process.core import constants
 from process.core import process_output as po
 from process.core.exceptions import ProcessValueError
+from process.core.model import Model
 from process.data_structure import (
     constraint_variables,
     current_drive_variables,
@@ -200,7 +201,7 @@ def ps_fraction_scene(beta: float) -> float:
     return -9e-2 * beta
 
 
-class Physics:
+class Physics(Model):
     def __init__(
         self,
         plasma_profile,
@@ -227,7 +228,12 @@ def __init__(
         self.plasma_transition = plasma_transition
         self.current = plasma_current
 
-    def physics(self):
+    def output(self):
+        self.calculate_effective_charge_ionisation_profiles()
+        self.outplas()
+        self.outtim()
+
+    def run(self):
         """Routine to calculate tokamak plasma physics information
         This routine calculates all the primary plasma physics parameters for a tokamak fusion reactor.
 
@@ -3234,7 +3240,7 @@ def outplas(self):
 
             self.inductance.output_volt_second_information()
         if stellarator_variables.istell == 0:
-            self.plasma_bootstrap_current.output_bootstrap_current_information()
+            self.plasma_bootstrap_current.output()
 
         po.osubhd(self.outfile, "Fuelling :")
         po.ovarre(
@@ -4735,7 +4741,7 @@ def output_volt_second_information(self):
         po.oblnkl(self.outfile)
 
 
-class DetailedPhysics:
+class DetailedPhysics(Model):
     """Class to hold detailed physics models for plasma processing."""
 
     def __init__(self, plasma_profile):
@@ -4743,6 +4749,10 @@ def __init__(self, plasma_profile):
         self.mfile = constants.MFILE
         self.plasma_profile = plasma_profile
 
+    def output(self):
+        self.run()
+        self.output_detailed_physics()
+
     def run(self):
         # ---------------------------
         #  Debye length calculation
diff --git a/process/models/physics/plasma_geometry.py b/process/models/physics/plasma_geometry.py
index 89edde767a..ef0ff60bc1 100644
--- a/process/models/physics/plasma_geometry.py
+++ b/process/models/physics/plasma_geometry.py
@@ -12,7 +12,7 @@ class PlasmaGeom:
     def __init__(self):
         self.outfile = constants.NOUT
 
-    def plasma_geometry(self):
+    def run(self):
         """Plasma geometry parameters
 
         This method calculates the plasma geometry parameters based on various shaping terms and input values.
diff --git a/process/models/power.py b/process/models/power.py
index 5fbac9962e..9a85ed9465 100644
--- a/process/models/power.py
+++ b/process/models/power.py
@@ -7,6 +7,7 @@
 from process.core import constants
 from process.core import process_output as po
 from process.core.exceptions import ProcessValueError
+from process.core.model import Model
 from process.data_structure import (
     build_variables,
     buildings_variables,
@@ -28,11 +29,46 @@
 logger = logging.getLogger(__name__)
 
 
-class Power:
+class Power(Model):
     def __init__(self):
         self.outfile = constants.NOUT
         self.mfile = constants.MFILE
 
+    def output(self):
+        # Toroidal field coil power model
+        self.tfpwr(output=True)
+
+        # Poloidal field coil power model !
+        self.pfpwr(output=True)
+
+        # Plant AC power requirements
+        self.acpow(output=True)
+
+        # Plant heat transport pt 2 & 3
+        self.output_cryogenics()
+        self.output_plant_thermal_powers()
+        self.output_plant_electric_powers()
+        self.output_power_profiles_over_time()
+
+    def run(self):
+        # Toroidal field coil power model
+        self.tfpwr(output=False)
+
+        # Poloidal field coil power model
+        self.pfpwr(output=False)
+
+        # Plant heat transport part 1
+        self.component_thermal_powers()
+
+        # Cryoplant loads
+        self.calculate_cryo_loads()
+
+        # Plant AC power requirements
+        self.acpow(output=False)
+
+        # Plant heat transport pt 2 & 3
+        self.plant_electric_production()
+
     def pfpwr(self, output: bool):
         """PF coil power supply requirements
 
diff --git a/process/models/pulse.py b/process/models/pulse.py
index 5c6e031586..ffe533442e 100644
--- a/process/models/pulse.py
+++ b/process/models/pulse.py
@@ -2,6 +2,7 @@
 
 from process.core import constants
 from process.core import process_output as po
+from process.core.model import Model
 from process.data_structure import (
     constraint_variables,
     numerics,
@@ -15,11 +16,14 @@
 logger = logging.getLogger(__name__)
 
 
-class Pulse:
+class Pulse(Model):
     def __init__(self):
         self.outfile = constants.NOUT
 
-    def run(self, output: bool):
+    def output(self):
+        self.run(output=True)
+
+    def run(self, output: bool = False):
         """Caller for the pulsed reactor model
 
         This calls the routines relevant to a pulsed reactor scenario.
diff --git a/process/models/shield.py b/process/models/shield.py
index 2d9a67ba04..d13a95be8f 100644
--- a/process/models/shield.py
+++ b/process/models/shield.py
@@ -2,6 +2,7 @@
 
 from process.core import constants
 from process.core import process_output as po
+from process.core.model import Model
 from process.data_structure import blanket_library as blanket_library
 from process.data_structure import build_variables as build_variables
 from process.data_structure import ccfe_hcpb_module as ccfe_hcpb_module
@@ -18,10 +19,13 @@
 logger = logging.getLogger(__name__)
 
 
-class Shield:
+class Shield(Model):
     def __init__(self):
         self.outfile = constants.NOUT
 
+    def output(self):
+        self.output_shld_areas_and_volumes()
+
     def run(self):
         blanket_library.dz_shld_half = self.calculate_shield_half_height(
             z_plasma_xpoint_lower=build_variables.z_plasma_xpoint_lower,
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index c7e904f169..74865ffeaa 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -12,6 +12,7 @@
 from process.core import process_output as po
 from process.core.coolprop_interface import FluidProperties
 from process.core.exceptions import ProcessValueError
+from process.core.model import Model
 from process.data_structure import (
     build_variables,
     constraint_variables,
@@ -58,7 +59,7 @@
 KEV = 1e3 * constants.ELECTRON_CHARGE  # Kiloelectron-volt (keV)
 
 
-class Stellarator:
+class Stellarator(Model):
     """Module containing stellarator routines
 
     This module contains routines for calculating the
@@ -118,7 +119,10 @@ def __init__(
         self.beta = plasma_beta
         self.bootstrap = plasma_bootstrap
 
-    def run(self, output: bool):
+    def output(self):
+        self.run(output=True)
+
+    def run(self, output: bool = False):
         """Routine to call the physics and engineering modules
         relevant to stellarators
 
diff --git a/process/models/structure.py b/process/models/structure.py
index af782d114c..bfc536978e 100644
--- a/process/models/structure.py
+++ b/process/models/structure.py
@@ -5,6 +5,7 @@
 
 from process.core import constants
 from process.core import process_output as po
+from process.core.model import Model
 from process.data_structure import build_variables as bv
 from process.data_structure import divertor_variables as divv
 from process.data_structure import fwbs_variables as fwbsv
@@ -16,7 +17,7 @@
 logger = logging.getLogger(__name__)
 
 
-class Structure:
+class Structure(Model):
     """Class containing support structure calculations
 
     This class contains routines for calculating the
@@ -27,6 +28,9 @@ class Structure:
     def __init__(self):
         self.outfile = constants.NOUT  # output file unit
 
+    def output(self):
+        self.run(output=True)
+
     def run(self, output: bool = False):
         """Structure calculation caller
 
diff --git a/process/models/tfcoil/base.py b/process/models/tfcoil/base.py
index 258f243d47..52edec0b20 100644
--- a/process/models/tfcoil/base.py
+++ b/process/models/tfcoil/base.py
@@ -12,6 +12,7 @@
 from process.core import constants
 from process.core import process_output as po
 from process.core.exceptions import ProcessValueError
+from process.core.model import Model
 from process.data_structure import (
     build_variables,
     fwbs_variables,
@@ -43,16 +44,18 @@ class TFCoilShapeModel(IntEnum):
     PICTURE_FRAME = 2
 
 
-class TFCoil:
+class TFCoil(Model):
     """Calculates the parameters of a resistive TF coil system for a fusion power plant"""
 
     def __init__(self, build: Build):
         """Initialise Fortran module variables."""
         self.outfile = constants.NOUT  # output file unit
-        self.iprint = 0  # switch for writing to output file (1=yes)
         self.build = build
         self.a_tf_inboard_total = tfcoil_variables.a_tf_inboard_total
 
+    def run(self):
+        self.cntrpst()
+
     def run_base_tf(self):
         """Run main tfcoil subroutine without outputting."""
 
@@ -126,8 +129,7 @@ def run_base_tf(self):
 
     def output(self):
         """Run main tfcoil subroutine and write output."""
-        self.iprint = 1
-        self.tfcoil(output=bool(self.iprint))
+        self.cntrpst(output=True)
 
     def tf_global_geometry(
         self,
@@ -2176,11 +2178,9 @@ def circumference(aaa, bbb):
             * (1.0e0 + (3.0e0 * hh) / (10.0e0 + np.sqrt(4.0e0 - 3.0e0 * hh)))
         )
 
-    def cntrpst(self):
+    def cntrpst(self, output: bool = False):
         """Evaluates the properties of a TART centrepost
 
-        outfile : input integer : output file unit
-        iprint : input integer : switch for writing to output file (1=yes)
         This subroutine evaluates the parameters of the centrepost for a
         tight aspect ratio tokamak. The centrepost is assumed to be tapered,
         i.e. narrowest on the midplane (z=0).
@@ -2468,7 +2468,7 @@ def cntrpst(self):
         presin = psat + dpres
 
         # Output section
-        if self.iprint == 1:
+        if output:
             po.oheadr(self.outfile, "Centrepost Coolant Parameters")
             po.ovarre(
                 self.outfile,
diff --git a/process/models/tfcoil/resistive.py b/process/models/tfcoil/resistive.py
index 43bb7dd4de..6cec6d4887 100644
--- a/process/models/tfcoil/resistive.py
+++ b/process/models/tfcoil/resistive.py
@@ -23,7 +23,10 @@ class ResistiveTFCoil(TFCoil):
     def __init__(self):
         self.outfile = constants.NOUT
 
-    def run(self, output: bool):
+    def output(self):
+        self.run(output=True)
+
+    def run(self, output: bool = False):
         """Run main tfcoil subroutine without outputting.
 
         Parameters
@@ -31,8 +34,6 @@ def run(self, output: bool):
         output: bool
 
         """
-        self.iprint = 0
-
         # Set up TF values share by all coil types
         self.run_base_tf()
 
diff --git a/process/models/tfcoil/superconducting.py b/process/models/tfcoil/superconducting.py
index 9732174f13..30c3887fb8 100644
--- a/process/models/tfcoil/superconducting.py
+++ b/process/models/tfcoil/superconducting.py
@@ -42,7 +42,11 @@ class SuperconductingTFCoil(TFCoil):
     def __init__(self):
         self.outfile = constants.NOUT
 
-    def run(self, output: bool):
+    def output(self):
+        self.run(output=True)
+        self.output_tf_superconductor_info()
+
+    def run(self, output: bool = False):
         """Routine to call the superconductor module for the TF coils
 
         Parameters
diff --git a/process/models/vacuum.py b/process/models/vacuum.py
index 73ba2ba25f..345153a9a7 100644
--- a/process/models/vacuum.py
+++ b/process/models/vacuum.py
@@ -5,6 +5,7 @@
 
 from process.core import constants, process_output
 from process.core import process_output as po
+from process.core.model import Model
 from process.data_structure import (
     blanket_library,
     build_variables,
@@ -20,7 +21,7 @@
 logger = logging.getLogger(__name__)
 
 
-class Vacuum:
+class Vacuum(Model):
     """Module containing vacuum system routines
 
     This module contains routines for calculating the
@@ -30,7 +31,10 @@ class Vacuum:
     def __init__(self):
         self.outfile: int = constants.NOUT
 
-    def run(self, output: bool):
+    def output(self):
+        self.run(output=True)
+
+    def run(self, output: bool = False):
         """Routine to call the vacuum module
         This routine calls the main vacuum package.
 
@@ -715,7 +719,7 @@ def vacuum(
         return pumpn, nduct, dlscalc, mvdsh, dimax
 
 
-class VacuumVessel:
+class VacuumVessel(Model):
     """Class containing vacuum vessel routines"""
 
     def __init__(self):
@@ -947,7 +951,7 @@ def calculate_elliptical_vessel_volumes(
         )
         return vol_vv_inboard, vol_vv_outboard, vol_vv
 
-    def output_vv_areas_and_volumes(self):
+    def output(self):
         """Output shield areas and volumes to log."""
 
         po.oheadr(self.outfile, "Vacuum Vessel Areas and Volumes")