diff --git a/process/blanket_library.py b/process/blanket_library.py
index 694c9211db..453c09c318 100644
--- a/process/blanket_library.py
+++ b/process/blanket_library.py
@@ -115,7 +115,7 @@ def component_half_height(self, icomponent: int):
 
         # Calculate component internal upper half-height (m)
         # If a double null machine then symmetric
-        if physics_variables.n_divertors == 2:
+        if divertor_variables.n_divertors == 2:
             htop = hbot
         else:
             # Blanket
@@ -324,7 +324,7 @@ def apply_coverage_factors(self):
         Apply coverage factors to volumes
         """
         # Apply blanket coverage factors
-        if physics_variables.n_divertors == 2:
+        if divertor_variables.n_divertors == 2:
             # double null configuration
             build_variables.a_blkt_outboard_surface = (
                 build_variables.a_blkt_total_surface
@@ -1192,7 +1192,7 @@ def blanket_module_poloidal_height(self):
 
             # Calculate blanket poloidal length and segment, subtracting divertor length (m)
             # kit hcll version only had the single null option
-            if physics_variables.n_divertors == 2:
+            if divertor_variables.n_divertors == 2:
                 # Double null configuration
                 blanket_library.len_blkt_outboard_segment_poloidal = (
                     0.5
@@ -1234,7 +1234,7 @@ def blanket_module_poloidal_height(self):
             # Assume divertor lies between the two ellipses, so fraction f_ster_div_single still applies
 
             # kit hcll version only had the single null option
-            if physics_variables.n_divertors == 2:
+            if divertor_variables.n_divertors == 2:
                 # Double null configuration
                 blanket_library.len_blkt_inboard_segment_poloidal = (
                     0.5
@@ -1264,7 +1264,7 @@ def blanket_module_poloidal_height(self):
 
             # kit hcll version only had the single null option
             # Calculate outboard blanket poloidal length and segment, subtracting divertor length (m)
-            if physics_variables.n_divertors == 2:
+            if divertor_variables.n_divertors == 2:
                 # Double null configuration
                 blanket_library.len_blkt_outboard_segment_poloidal = (
                     0.5
diff --git a/process/build.py b/process/build.py
index 0317324e3c..f8adcf3f55 100644
--- a/process/build.py
+++ b/process/build.py
@@ -987,7 +987,7 @@ def divgeom(self, output: bool):
         divht = max(zplti, zplto) - min(zplbo, zplbi)
 
         if output:
-            if physics_variables.n_divertors == 1:
+            if divertor_variables.n_divertors == 1:
                 po.oheadr(self.outfile, "Divertor build and plasma position")
                 po.ocmmnt(self.outfile, "Divertor Configuration = Single Null Divertor")
                 po.oblnkl(self.outfile)
@@ -1200,7 +1200,7 @@ def divgeom(self, output: bool):
                     "OP ",
                 )
 
-            elif physics_variables.n_divertors == 2:
+            elif divertor_variables.n_divertors == 2:
                 po.oheadr(self.outfile, "Divertor build and plasma position")
                 po.ocmmnt(self.outfile, "Divertor Configuration = Double Null Divertor")
                 po.oblnkl(self.outfile)
@@ -2018,7 +2018,7 @@ def calculate_radial_build(self, output: bool) -> None:
             - 0.5e0 * (build_variables.dr_fw_inboard + build_variables.dr_fw_outboard)
         )
         if (
-            physics_variables.n_divertors == 2
+            divertor_variables.n_divertors == 2
         ):  # (i.e. physics_variables.i_single_null=0)
             htop = hbot
         else:
@@ -2089,7 +2089,7 @@ def calculate_radial_build(self, output: bool) -> None:
 
         #  Apply area coverage factor
 
-        if physics_variables.n_divertors == 2:
+        if divertor_variables.n_divertors == 2:
             # Double null configuration
             build_variables.a_fw_outboard = (
                 build_variables.a_fw_outboard_full_coverage
diff --git a/process/data_structure/divertor_variables.py b/process/data_structure/divertor_variables.py
index 091932a04b..0102992aea 100644
--- a/process/data_structure/divertor_variables.py
+++ b/process/data_structure/divertor_variables.py
@@ -56,6 +56,21 @@
 xpertin: float = None
 """perpendicular heat transport coefficient (m2/s)"""
 
+p_div_lower_nuclear_heat_mw: float = None
+"""Lower divertor neutron nuclear heat load on (MW)"""
+
+p_div_upper_nuclear_heat_mw: float = None
+"""Upper divertor neutron nuclear heat load on (MW)"""
+
+p_div_upper_rad_mw: float = None
+"""Upper divertor incident radiation power radiation power (MW)"""
+
+p_div_lower_rad_mw: float = None
+"""Lower divertor incident radiation power radiation power (MW)"""
+
+n_divertors: int = None
+"""Number of divertors (calculated from `i_single_null`)"""
+
 
 def init_divertor_variables():
     global anginc
@@ -76,6 +91,11 @@ def init_divertor_variables():
     global prn1
     global tdiv
     global xpertin
+    global p_div_lower_nuclear_heat_mw
+    global p_div_upper_nuclear_heat_mw
+    global p_div_upper_rad_mw
+    global p_div_lower_rad_mw
+    global n_divertors
 
     anginc = 0.262
     deg_div_field_plate = 1.0
@@ -95,3 +115,8 @@ def init_divertor_variables():
     prn1 = 0.285
     tdiv = 2.0
     xpertin = 2.0
+    p_div_lower_nuclear_heat_mw = 0.0
+    p_div_upper_nuclear_heat_mw = 0.0
+    p_div_upper_rad_mw = 0.0
+    p_div_lower_rad_mw = 0.0
+    n_divertors = 2
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index ef365893bc..8e07783e76 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -629,10 +629,6 @@
 """
 
 
-n_divertors: int = None
-"""number of divertors (calculated from `i_single_null`)"""
-
-
 i_beta_fast_alpha: int = None
 """switch for fast alpha pressure calculation
 - =0 ITER physics rules (Uckan) fit
@@ -1486,7 +1482,6 @@ def init_physics_variables():
     global i_plasma_current
     global i_diamagnetic_current
     global i_density_limit
-    global n_divertors
     global i_beta_fast_alpha
     global i_plasma_ignited
     global i_plasma_pedestal
@@ -1746,7 +1741,6 @@ def init_physics_variables():
     i_plasma_current = 4
     i_diamagnetic_current = 0
     i_density_limit = 8
-    n_divertors = 2
     i_beta_fast_alpha = 1
     i_plasma_ignited = 0
     i_plasma_pedestal = 1
diff --git a/process/dcll.py b/process/dcll.py
index e1a255dadb..514d4f15c3 100644
--- a/process/dcll.py
+++ b/process/dcll.py
@@ -8,6 +8,7 @@
     build_variables,
     current_drive_variables,
     dcll_variables,
+    divertor_variables,
     fwbs_variables,
     heat_transport_variables,
     physics_variables,
@@ -134,7 +135,7 @@ def dcll_neutronics_and_power(self, output: bool):
              - f_nuc_pow_bz_liq
         """
 
-        if physics_variables.n_divertors == 2:
+        if divertor_variables.n_divertors == 2:
             # Double null configuration
             covf = (
                 1
diff --git a/process/divertor.py b/process/divertor.py
index 738afae3ff..f833b2fc6e 100644
--- a/process/divertor.py
+++ b/process/divertor.py
@@ -39,13 +39,13 @@ def run(self, output: bool) -> None:
         fwbs.p_div_nuclear_heat_total_mw = self.incident_neutron_power(
             p_plasma_neutron_mw=pv.p_plasma_neutron_mw,
             f_ster_div_single=fwbs.f_ster_div_single,
-            n_divertors=pv.n_divertors,
+            n_divertors=dv.n_divertors,
         )
 
         fwbs.p_div_rad_total_mw = self.incident_radiation_power(
             p_plasma_rad_mw=pv.p_plasma_rad_mw,
             f_ster_div_single=fwbs.f_ster_div_single,
-            n_divertors=pv.n_divertors,
+            n_divertors=dv.n_divertors,
         )
 
         if dv.i_div_heat_load == 0 and output:
@@ -346,7 +346,7 @@ def divwade(
         hldiv_base = p_plasma_separatrix_mw * (1 - rad_fraction_sol) / area_wetted
 
         # For double null, calculate heat loads to upper and lower divertors and use the highest
-        if pv.n_divertors == 2:
+        if dv.n_divertors == 2:
             hldiv_lower = f_p_div_lower * hldiv_base
             hldiv_upper = (1.0 - f_p_div_lower) * hldiv_base
             dv.pflux_div_heat_load_mw = max(hldiv_lower, hldiv_upper)
@@ -423,6 +423,36 @@ def incident_neutron_power(
 class LowerDivertor(Divertor):
     """Module containing lower divertor routines"""
 
+    def run(self, output: bool) -> None:
+        super().run(output=output)
+
+        dv.p_div_lower_nuclear_heat_mw = self.incident_neutron_power(
+            p_plasma_neutron_mw=pv.p_plasma_neutron_mw,
+            f_ster_div_single=fwbs.f_ster_div_single,
+            n_divertors=1,
+        )
+
+        dv.p_div_lower_rad_mw = self.incident_radiation_power(
+            p_plasma_rad_mw=pv.p_plasma_rad_mw,
+            f_ster_div_single=fwbs.f_ster_div_single,
+            n_divertors=1,
+        )
+
 
 class UpperDivertor(Divertor):
     """Module containing upper divertor routines"""
+
+    def run(self, output: bool) -> None:
+        super().run(output=output)
+
+        dv.p_div_upper_nuclear_heat_mw = self.incident_neutron_power(
+            p_plasma_neutron_mw=pv.p_plasma_neutron_mw,
+            f_ster_div_single=fwbs.f_ster_div_single,
+            n_divertors=1,
+        )
+
+        dv.p_div_upper_rad_mw = self.incident_radiation_power(
+            p_plasma_rad_mw=pv.p_plasma_rad_mw,
+            f_ster_div_single=fwbs.f_ster_div_single,
+            n_divertors=1,
+        )
diff --git a/process/hcpb.py b/process/hcpb.py
index 7c8d29c312..288beb7a51 100644
--- a/process/hcpb.py
+++ b/process/hcpb.py
@@ -177,7 +177,7 @@ def run(self, output: bool):
         # Solid angle fraction taken by the breeding blankets/shields
         f_geom_blanket = (
             1
-            - physics_variables.n_divertors * fwbs_variables.f_ster_div_single
+            - divertor_variables.n_divertors * fwbs_variables.f_ster_div_single
             - f_geom_cp
         )
 
@@ -323,7 +323,7 @@ def component_masses(self):
             * physics_variables.rmajor
             * physics_variables.rminor
         )
-        if physics_variables.n_divertors == 2:
+        if divertor_variables.n_divertors == 2:
             divertor_variables.a_div_surface_total = (
                 divertor_variables.a_div_surface_total * 2.0
             )
diff --git a/process/init.py b/process/init.py
index 15b7598cc0..9bdd1c4b4c 100644
--- a/process/init.py
+++ b/process/init.py
@@ -600,7 +600,7 @@ def check_process(inputs):  # noqa: ARG001
             )
 
     if data_structure.physics_variables.i_single_null == 0:
-        data_structure.physics_variables.n_divertors = 2
+        data_structure.divertor_variables.n_divertors = 2
         data_structure.build_variables.dz_fw_plasma_gap = (
             data_structure.build_variables.dz_xpoint_divertor
         )
@@ -612,7 +612,7 @@ def check_process(inputs):  # noqa: ARG001
         )
         warn("Double-null: Upper vertical build forced to match lower", stacklevel=2)
     else:  # i_single_null == 1
-        data_structure.physics_variables.n_divertors = 1
+        data_structure.divertor_variables.n_divertors = 1
 
     #  Tight aspect ratio options (ST)
     if data_structure.physics_variables.itart == 1:
diff --git a/process/physics.py b/process/physics.py
index 3227317fa1..acdfaeda6f 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -2373,7 +2373,7 @@ def physics(self):
                 / physics_variables.a_plasma_surface
             )
         else:
-            if physics_variables.n_divertors == 2:
+            if divertor_variables.n_divertors == 2:
                 # Double null configuration
                 physics_variables.pflux_fw_neutron_mw = (
                     (
@@ -2507,7 +2507,7 @@ def physics(self):
         # if double null configuration share the power
         # over the upper and lower divertor, where physics_variables.f_p_div_lower gives
         # the factor of power conducted to the lower divertor
-        if physics_variables.n_divertors == 2:
+        if divertor_variables.n_divertors == 2:
             physics_variables.p_div_lower_separatrix_mw = (
                 physics_variables.f_p_div_lower
                 * physics_variables.p_plasma_separatrix_mw
@@ -2738,7 +2738,7 @@ def physics(self):
                 / physics_variables.a_plasma_surface
             )
         else:
-            if physics_variables.n_divertors == 2:
+            if divertor_variables.n_divertors == 2:
                 # Double Null configuration in - including SoL radiation
                 physics_variables.pflux_fw_rad_mw = (
                     (
@@ -2820,7 +2820,7 @@ def physics(self):
                 )
             )
         )
-        if physics_variables.n_divertors == 2:
+        if divertor_variables.n_divertors == 2:
             # Double Null configuration
             # Find all the power fractions accross the targets
             # Taken from D3-D conventional divertor design
@@ -4120,16 +4120,16 @@ def outplas(self):
         po.oheadr(self.outfile, "Plasma")
 
         if stellarator_variables.istell == 0:
-            if physics_variables.n_divertors == 0:
+            if divertor_variables.n_divertors == 0:
                 po.ocmmnt(self.outfile, "Plasma configuration = limiter")
-            elif physics_variables.n_divertors == 1:
+            elif divertor_variables.n_divertors == 1:
                 po.ocmmnt(self.outfile, "Plasma configuration = single null divertor")
-            elif physics_variables.n_divertors == 2:
+            elif divertor_variables.n_divertors == 2:
                 po.ocmmnt(self.outfile, "Plasma configuration = double null divertor")
             else:
                 raise ProcessValueError(
                     "Illegal value of n_divertors",
-                    n_divertors=physics_variables.n_divertors,
+                    n_divertors=divertor_variables.n_divertors,
                 )
         else:
             po.ocmmnt(self.outfile, "Plasma configuration = stellarator")
@@ -5719,7 +5719,7 @@ def outplas(self):
                 physics_variables.lambdaio,
                 "OP ",
             )
-            if physics_variables.n_divertors == 2:
+            if divertor_variables.n_divertors == 2:
                 po.ovarre(
                     self.outfile,
                     "Midplane seperation of the two magnetic closed flux surfaces (m)",
@@ -5749,7 +5749,7 @@ def outplas(self):
                 physics_variables.flo,
                 "OP ",
             )
-            if physics_variables.n_divertors == 2:
+            if divertor_variables.n_divertors == 2:
                 po.ovarre(
                     self.outfile,
                     "Fraction of power incident on the upper inner target",
@@ -5779,7 +5779,7 @@ def outplas(self):
                 physics_variables.plomw,
                 "OP ",
             )
-            if physics_variables.n_divertors == 2:
+            if divertor_variables.n_divertors == 2:
                 po.ovarre(
                     self.outfile,
                     "Power incident on the upper innner target (MW)",
@@ -5889,7 +5889,7 @@ def outplas(self):
             )
             po.oblnkl(self.outfile)
 
-        if physics_variables.n_divertors == 2:
+        if divertor_variables.n_divertors == 2:
             # Double null divertor configuration
             po.ovarre(
                 self.outfile,
diff --git a/process/vacuum.py b/process/vacuum.py
index 2e884980e7..85cb21c27a 100644
--- a/process/vacuum.py
+++ b/process/vacuum.py
@@ -6,6 +6,7 @@
 from process import constants
 from process import process_output as po
 from process.data_structure import build_variables as buv
+from process.data_structure import divertor_variables as dv
 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
@@ -74,7 +75,7 @@ def run(self, output: bool) -> None:
                 tfv.n_tf_coils,
                 tv.t_plant_pulse_dwell,
                 pv.nd_plasma_electrons_vol_avg,
-                pv.n_divertors,
+                dv.n_divertors,
                 qtorus,
                 gasld,
                 output=output,
diff --git a/tests/unit/test_blanket_library.py b/tests/unit/test_blanket_library.py
index 24b8b34910..6cce8aa5b1 100644
--- a/tests/unit/test_blanket_library.py
+++ b/tests/unit/test_blanket_library.py
@@ -488,7 +488,7 @@ def test_component_half_height(
         componenthalfheightparam.z_plasma_xpoint_upper,
     )
     monkeypatch.setattr(
-        physics_variables, "n_divertors", componenthalfheightparam.n_divertors
+        divertor_variables, "n_divertors", componenthalfheightparam.n_divertors
     )
     monkeypatch.setattr(
         divertor_variables, "dz_divertor", componenthalfheightparam.dz_divertor
@@ -1374,7 +1374,7 @@ def test_apply_coverage_factors(
     monkeypatch.setattr(fwbs_variables, "vol_vv", applycoveragefactorsparam.vol_vv)
     monkeypatch.setattr(fwbs_variables, "fvoldw", applycoveragefactorsparam.fvoldw)
     monkeypatch.setattr(
-        physics_variables, "n_divertors", applycoveragefactorsparam.n_divertors
+        divertor_variables, "n_divertors", applycoveragefactorsparam.n_divertors
     )
     monkeypatch.setattr(
         blanket_library, "vol_shld_inboard", applycoveragefactorsparam.vol_shld_inboard
@@ -1523,7 +1523,7 @@ def test_blanket_mod_pol_height(
     monkeypatch.setattr(physics_variables, "itart", blanketmodpolheightparam.itart)
     monkeypatch.setattr(physics_variables, "rminor", blanketmodpolheightparam.rminor)
     monkeypatch.setattr(
-        physics_variables, "n_divertors", blanketmodpolheightparam.n_divertors
+        divertor_variables, "n_divertors", blanketmodpolheightparam.n_divertors
     )
     monkeypatch.setattr(physics_variables, "rmajor", blanketmodpolheightparam.rmajor)
     monkeypatch.setattr(physics_variables, "triang", blanketmodpolheightparam.triang)
diff --git a/tests/unit/test_build.py b/tests/unit/test_build.py
index b301f4f96e..4ec96d2069 100644
--- a/tests/unit/test_build.py
+++ b/tests/unit/test_build.py
@@ -141,7 +141,7 @@ def test_divgeom(divgeomparam, monkeypatch, build):
 
     monkeypatch.setattr(physics_variables, "rminor", divgeomparam.rminor)
 
-    monkeypatch.setattr(physics_variables, "n_divertors", divgeomparam.n_divertors)
+    monkeypatch.setattr(divertor_variables, "n_divertors", divgeomparam.n_divertors)
 
     monkeypatch.setattr(physics_variables, "kappa", divgeomparam.kappa)
 
diff --git a/tests/unit/test_ccfe_hcpb.py b/tests/unit/test_ccfe_hcpb.py
index 50cc429bce..9983c2468a 100644
--- a/tests/unit/test_ccfe_hcpb.py
+++ b/tests/unit/test_ccfe_hcpb.py
@@ -998,7 +998,7 @@ def test_powerflow_calc(powerflowcalcparam, monkeypatch, ccfe_hcpb):
     )
 
     monkeypatch.setattr(
-        physics_variables, "n_divertors", powerflowcalcparam.n_divertors
+        divertor_variables, "n_divertors", powerflowcalcparam.n_divertors
     )
 
     monkeypatch.setattr(
@@ -1469,7 +1469,7 @@ def test_component_masses(componentmassesparam, monkeypatch, ccfe_hcpb):
     monkeypatch.setattr(physics_variables, "rminor", componentmassesparam.rminor)
     monkeypatch.setattr(physics_variables, "rmajor", componentmassesparam.rmajor)
     monkeypatch.setattr(
-        physics_variables, "n_divertors", componentmassesparam.n_divertors
+        divertor_variables, "n_divertors", componentmassesparam.n_divertors
     )
     monkeypatch.setattr(
         physics_variables, "a_plasma_surface", componentmassesparam.a_plasma_surface
diff --git a/tests/unit/test_dcll.py b/tests/unit/test_dcll.py
index a1a33e8b87..b479f0eabd 100644
--- a/tests/unit/test_dcll.py
+++ b/tests/unit/test_dcll.py
@@ -6,6 +6,7 @@
     build_variables,
     current_drive_variables,
     dcll_variables,
+    divertor_variables,
     fwbs_variables,
     physics_variables,
 )
@@ -272,7 +273,7 @@ def test_dcll_neutronics_and_power(dcllneutronicsandpowerparam, monkeypatch, dcl
     )
 
     monkeypatch.setattr(
-        physics_variables, "n_divertors", dcllneutronicsandpowerparam.n_divertors
+        divertor_variables, "n_divertors", dcllneutronicsandpowerparam.n_divertors
     )
 
     monkeypatch.setattr(
diff --git a/tests/unit/test_power.py b/tests/unit/test_power.py
index 506c6105e4..72415be514 100644
--- a/tests/unit/test_power.py
+++ b/tests/unit/test_power.py
@@ -9,6 +9,7 @@
     constraint_variables,
     cost_variables,
     current_drive_variables,
+    divertor_variables,
     fwbs_variables,
     heat_transport_variables,
     numerics,
@@ -2805,7 +2806,7 @@ def test_power2(power2param, monkeypatch, power):
 
     monkeypatch.setattr(physics_variables, "p_fw_alpha_mw", power2param.p_fw_alpha_mw)
 
-    monkeypatch.setattr(physics_variables, "n_divertors", power2param.n_divertors)
+    monkeypatch.setattr(divertor_variables, "n_divertors", power2param.n_divertors)
 
     monkeypatch.setattr(
         physics_variables, "p_plasma_ohmic_mw", power2param.p_plasma_ohmic_mw