diff --git a/process/caller.py b/process/caller.py
index 0a93243a15..3f5b0a0eff 100644
--- a/process/caller.py
+++ b/process/caller.py
@@ -276,6 +276,8 @@ def _call_models_once(self, xc: np.ndarray) -> None:
         # First wall model
         self.models.fw.run()
 
+        self.models.vacuum_vessel.run()
+
         # Blanket model
         """Blanket switch values
         No.  |  model
diff --git a/process/main.py b/process/main.py
index 2b0c1028bd..393990d1b0 100644
--- a/process/main.py
+++ b/process/main.py
@@ -106,7 +106,7 @@
 from process.structure import Structure
 from process.superconducting_tf_coil import SuperconductingTFCoil
 from process.tf_coil import TFCoil
-from process.vacuum import Vacuum
+from process.vacuum import Vacuum, VacuumVessel
 from process.water_use import WaterUse
 
 os.environ["PYTHON_PROCESS_ROOT"] = os.path.join(os.path.dirname(__file__))
@@ -660,6 +660,7 @@ def __init__(self):
         self.availability = Availability()
         self.buildings = Buildings()
         self.vacuum = Vacuum()
+        self.vacuum_vessel = VacuumVessel()
         self.water_use = WaterUse()
         self.pulse = Pulse()
         self.ife = IFE(availability=self.availability, costs=self.costs)
diff --git a/process/vacuum.py b/process/vacuum.py
index 2e884980e7..c48463d07a 100644
--- a/process/vacuum.py
+++ b/process/vacuum.py
@@ -5,7 +5,10 @@
 
 from process import constants
 from process import process_output as po
+from process.blanket_library import dshellvol, eshellvol
+from process.data_structure import blanket_library as blanket_library
 from process.data_structure import build_variables as buv
+from process.data_structure import fwbs_variables as fwbs_variables
 from process.data_structure import physics_variables as pv
 from process.data_structure import tfcoil_variables as tfv
 from process.data_structure import times_variables as tv
@@ -689,3 +692,141 @@ def vacuum(
             )
 
         return pumpn, nduct, dlscalc, mvdsh, dimax
+
+
+class VacuumVessel:
+    """Class containing vacuum vessel routines"""
+
+    def __init__(self) -> None:
+        pass
+
+    def run(self) -> None:
+        blanket_library.dz_vv_half = self.calculate_vessel_half_height(
+            z_tf_inside_half=buv.z_tf_inside_half,
+            dz_shld_vv_gap=buv.dz_shld_vv_gap,
+            dz_vv_lower=buv.dz_vv_lower,
+            n_divertors=pv.n_divertors,
+            dz_blkt_upper=buv.dz_blkt_upper,
+            dz_shld_upper=buv.dz_shld_upper,
+            z_plasma_xpoint_upper=buv.z_plasma_xpoint_upper,
+            dr_fw_plasma_gap_inboard=buv.dr_fw_plasma_gap_inboard,
+            dr_fw_plasma_gap_outboard=buv.dr_fw_plasma_gap_outboard,
+            dr_fw_inboard=buv.dr_fw_inboard,
+            dr_fw_outboard=buv.dr_fw_outboard,
+        )
+
+        (
+            blanket_library.vol_vv_inboard,
+            blanket_library.vol_vv_outboard,
+            fwbs_variables.vol_vv,
+        ) = self.calculate_dshaped_vessel_volumes(
+            rsldi=buv.rsldi,
+            rsldo=buv.rsldo,
+            dz_vv_half=blanket_library.dz_vv_half,
+            dr_vv_inboard=buv.dr_vv_inboard,
+            dr_vv_outboard=buv.dr_vv_outboard,
+            dz_vv_upper=buv.dz_vv_upper,
+            dz_vv_lower=buv.dz_vv_lower,
+        )
+
+    def calculate_vessel_half_height(
+        self,
+        z_tf_inside_half: float,
+        dz_shld_vv_gap: float,
+        dz_vv_lower: float,
+        n_divertors: int,
+        dz_blkt_upper: float,
+        dz_shld_upper: float,
+        z_plasma_xpoint_upper: float,
+        dr_fw_plasma_gap_inboard: float,
+        dr_fw_plasma_gap_outboard: float,
+        dr_fw_inboard: float,
+        dr_fw_outboard: float,
+    ) -> float:
+        """Calculate vacuum vessel internal half-height (m)"""
+
+        z_bottom = z_tf_inside_half - dz_shld_vv_gap - dz_vv_lower
+
+        # Calculate component internal upper half-height (m)
+        # If a double null machine then symmetric
+        if n_divertors == 2:
+            z_top = z_bottom
+        else:
+            z_top = z_plasma_xpoint_upper + 0.5 * (
+                dr_fw_plasma_gap_inboard
+                + dr_fw_plasma_gap_outboard
+                + dr_fw_inboard
+                + dr_fw_outboard
+                + dz_blkt_upper
+                + dz_shld_upper
+            )
+
+        # Average of top and bottom (m)
+        return 0.5 * (z_top + z_bottom)
+
+    def calculate_dshaped_vessel_volumes(
+        self,
+        rsldi: float,
+        rsldo: float,
+        dz_vv_half: float,
+        dr_vv_inboard: float,
+        dr_vv_outboard: float,
+        dz_vv_upper: float,
+        dz_vv_lower: float,
+    ) -> tuple[float, float, float]:
+        """Calculate volumes of D-shaped vacuum vessel segments"""
+
+        r_1 = rsldi
+        r_2 = rsldo - r_1
+
+        (
+            vol_vv_inboard,
+            vol_vv_outboard,
+            vol_vv,
+        ) = dshellvol(
+            rmajor=r_1,
+            rminor=r_2,
+            zminor=dz_vv_half,
+            drin=dr_vv_inboard,
+            drout=dr_vv_outboard,
+            dz=(dz_vv_upper + dz_vv_lower) / 2,
+        )
+
+        return vol_vv_inboard, vol_vv_outboard, vol_vv
+
+    def calculate_elliptical_vessel_volumes(
+        self,
+        rmajor: float,
+        rminor: float,
+        triang: float,
+        rsldi: float,
+        rsldo: float,
+        dz_vv_half: float,
+        dr_vv_inboard: float,
+        dr_vv_outboard: float,
+        dz_vv_upper: float,
+        dz_vv_lower: float,
+    ) -> tuple[float, float, float]:
+        # Major radius to centre of inboard and outboard ellipses (m)
+        # (coincident in radius with top of plasma)
+        r_1 = rmajor - rminor * triang
+
+        # Calculate distance between r1 and outer edge of inboard ...
+        # ... section (m)
+        r_2 = r_1 - rsldi
+        r_3 = rsldo - r_1
+
+        (
+            vol_vv_inboard,
+            vol_vv_outboard,
+            vol_vv,
+        ) = eshellvol(
+            r_1,
+            r_2,
+            r_3,
+            dz_vv_half,
+            dr_vv_inboard,
+            dr_vv_outboard,
+            (dz_vv_upper + dz_vv_lower) / 2,
+        )
+        return vol_vv_inboard, vol_vv_outboard, vol_vv