diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index c30155c8a4..df27334cf0 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -292,6 +292,9 @@
 b_plasma_total: float = None
 """Sum of plasma total toroidal + poloidal field (T)"""
 
+e_plasma_magnetic_stored: float = None
+"""Plasma stored magnetic energy (J)"""
+
 
 burnup: float = None
 """fractional plasma burnup"""
@@ -1392,6 +1395,7 @@ def init_physics_variables():
     global b_plasma_poloidal_average
     global b_plasma_toroidal_on_axis
     global b_plasma_total
+    global e_plasma_magnetic_stored
     global burnup
     global burnup_in
     global b_plasma_vertical_required
@@ -1644,6 +1648,7 @@ def init_physics_variables():
     b_plasma_poloidal_average = 0.0
     b_plasma_toroidal_on_axis = 5.68
     b_plasma_total = 0.0
+    e_plasma_magnetic_stored = 0.0
     burnup = 0.0
     burnup_in = 0.0
     b_plasma_vertical_required = 0.0
diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index b2a418eae5..adcd0b016c 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -2681,7 +2681,8 @@ def plot_main_plasma_information(
         f"$\\Omega_{{\\text{{p}}}}$: {mfile_data.data['res_plasma'].get_scan(scan):.4e} $\\Omega$\n"
         f"Plasma resistive diffusion time: {mfile_data.data['t_plasma_res_diffusion'].get_scan(scan):.4f} s\n"
         f"Plasma inductance: {mfile_data.data['ind_plasma'].get_scan(scan):.4e} H\n"
-        f"Plasma normalised internal inductance: {mfile_data.data['ind_plasma_internal_norm'].get_scan(scan):.4f}\n"
+        f"Plasma stored magnetic energy: {mfile_data.data['e_plasma_magnetic_stored'].get_scan(scan) / 1e9:.4f} GJ\n"
+        f"Plasma normalised internal inductance: {mfile_data.data['ind_plasma_internal_norm'].get_scan(scan):.4f}"
     )
 
     axis.text(
diff --git a/process/physics.py b/process/physics.py
index b3c8867600..9d9ab8ec83 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -2544,6 +2544,10 @@ def physics(self):
             physics_variables.ind_plasma_internal_norm,
         )
 
+        physics_variables.e_plasma_magnetic_stored = (
+            0.5e0 * physics_variables.ind_plasma * physics_variables.plasma_current**2
+        )
+
         # Calculate auxiliary physics related information
         sbar = 1.0e0
         (
@@ -6234,6 +6238,13 @@ def outplas(self):
                 physics_variables.ind_plasma,
                 "OP ",
             )
+            po.ovarre(
+                self.outfile,
+                "Plasma magnetic energy stored (J)",
+                "(e_plasma_magnetic_stored)",
+                physics_variables.e_plasma_magnetic_stored,
+                "OP ",
+            )
             po.ovarrf(
                 self.outfile,
                 "Plasma normalised internal inductance",