diff --git a/documentation/eng-models/pf-coil.md b/documentation/eng-models/pf-coil.md
index 493a823256..8160d2b04b 100644
--- a/documentation/eng-models/pf-coil.md
+++ b/documentation/eng-models/pf-coil.md
@@ -9,37 +9,101 @@ position and shape during the flat-top period.
 The positions and sizes of te PF coils are partly input, and partly calculated after consideration 
 of the required currents and allowable current density.
 
-The PF coil locations are controlled using a set of switched stored in array `i_pf_location` (see 
-Figure 1), and are calculated in routine `PFCOIL`. The coils are (usually) organised into groups 
+The PF coil locations are controlled using a set of switches stored in array `i_pf_location[]`. The coils are (usually) organised into groups 
 containing two PF coils placed symmetrically above and below the midplane, and each group `j` has 
 an element `i_pf_location(j)` assigned to it. Input parameter `n_pf_coil_groups` should be set to the number of groups, 
 and `n_pf_coils_in_group(j)` should be assigned the number of coils in each group - which should be 2 in each case.
 
-<figure markdown>
-![Machine build](../images/vertical-build.png){ width="100%"}
-<figcaption>Figure 1: Machine build for D-shaped major components</figcaption>
-</figure>
-
-In the following, all variables are defined in the variable descriptor file `vardes.html`. The 
-values for `rpf1`, `rpf2`, `zref(j)` and `routr` should be adjusted by the user to locate the PF 
+The values for `dr_pf_cs_middle_offset`, `rpf2`, `zref(j)` and `dr_pf_tf_outboard_out_offset` should be adjusted by the user to locate the PF 
 coils accurately.
 
-The three possible values of `i_pf_location(j)` correspond to the following PF coil positions: (Redo taking 
-into account `i_single_null` and other recent changes e.g. rclsnorm)
+For PF coils going to be placed outside the TF coils the key radius is defined as:
+
+$$
+\overbrace{R_{\text{PF, outside-TF}}}^{\texttt{r_pf_outside_tf_midplane}} = \overbrace{R_{\text{TF,outboard-out}}}^{\texttt{r_tf_outboard_out}} + \overbrace{dR_{\text{PF,TF-offset}}}^{\texttt{dr_tf_outboard_out_offset}}
+$$
+
+The four possible values of `i_pf_location(j)` correspond to the following PF coil positions: 
+
+### Above the central solenoid (one group only) | `place_pf_above_cs()`
+
+- `i_pf_location(j) = 1`;
+
+    $$
+    R = \overbrace{R_{\text{CS,middle}}}^{\texttt{r_cs_middle}} + \overbrace{dR_{\text{offset}}}^{\texttt{dr_pf_cs_middle_offset}}
+    $$
+
+    $$
+    Z = \pm \ \overbrace{Z_{\text{CS,top}}}^{\texttt{z_cs_coil_upper}} + 0.1 
+    \\ + \frac{1}{2}\left(\left(\underbrace{Z_{\text{TF,inside}}}_{\texttt{z_tf_inside_half}}-Z_{\text{CS,top}}\right)+\underbrace{dR_{\text{TF,inboard}}}_{\texttt{dr_tf_inboard}}+0.1\right)
+    $$
+
+-----------------------------
+
+### Above the TF coils (one group only) | `place_pf_above_tf()`
+
+- `i_pf_location(j) = 2`;
+
+    $$
+    R = R_0 + \texttt{rpf2} \times \delta \times a
+    $$
+
+    Due to the nautre of the TF coils not being top-down symmetric for single null cases the positioning slightly differs if the coil is above or below the midplane.
+    
+    For above the midplane:
+
+    $$
+    Z = Z_{\text{TF,top}} + 0.86
+    $$    
+
+    For below the midplane:
+
+    $$
+    Z = - \left(Z_{\text{TF,top}} - \overbrace{dZ_{\text{TF,upper-lower-midplane}}}^{\texttt{dz_tf_upper_lower_midplane}} + 0.86  \right)
+    $$
+
+-------------------
+
+
+### Outside the TF coils | `place_pf_outside_tf()`
+
+- `i_pf_location(j) = 3`;
+
+    The PF coils can either be stacked vertically outside the TF (ideal for a picture frame coil) or follow the TF coil curve (D-shaped):
+
+    If the chosen value for `i_tf_shape` is that of a picture frame or `i_r_pf_outside_tf_placement == 1` then:
+
+    $$
+    R = \overbrace{R_{\text{PF, outside-TF}}}^{\texttt{r_pf_outside_tf_midplane}}
+    $$
+
+    Else, the coils follow a D-shape curve
+
+    $$
+    R =  \sqrt{\left(\overbrace{R_{\text{PF, outside-TF}}}^{\texttt{r_pf_outside_tf_midplane}}\right)^2 - Z^2}
+    $$
+
+    $$
+    Z = \pm \ a \times \texttt{zref}
+    $$
+
+---------------------
+
+### General placement | `place_pf_generally()`
+
+- `i_pf_location(j) = 4`;
 
-`i_pf_location(j)` = 1: PF coils are placed above the central solenoid (one group only);
-*R* = `r_cs_middle` + `rpf1`<br>
-*Z* = $\pm$(`z_tf_inside_half` * `f_z_cs_tf_internal` + 0.1 + 0.5 * (`z_tf_inside_half` * (1 - `f_z_cs_tf_internal`) + `dr_tf_inboard` + 0.1))
+    The PF coils are placed generally in units of minor radius ($a$) relative to the mid-plane and plasma major radius ($R_0$);
 
-`i_pf_location(j)` = 2: PF coils are placed above the TF coils (one group only);<br>
-*R* = `rmajor` + `rpf2`<br>
-*Z* = $\pm$(`z_tf_inside_half` * `dr_tf_inboard` + 0.86)
+    $$
+    R = \pm \ a \times \texttt{rref[]} + R_0
+    $$
 
-`i_pf_location(j)` = 3: PF coils are placed radially outside the TF coils (any number of groups);<br>
-*R* = `rtot` + `dr_tf_outboard`/2 + `routr`<br>
-*Z* = $\pm$(`rminor` * `zref(j)`
+    $$
+    Z = \pm \ a \times \texttt{zref[]}
+    $$    
 
-The void fraction (for coolant) in each coil `i`'s winding pack is given by `f_a_pf_coil_void(i)`.
+-------------------------
 
 ## Coil currents
 
diff --git a/process/build.py b/process/build.py
index 3e97ede4f1..18ff3c962e 100644
--- a/process/build.py
+++ b/process/build.py
@@ -826,7 +826,7 @@ def calculate_vertical_build(self, output: bool) -> None:
             build_variables.z_tf_top = (
                 build_variables.z_tf_inside_half + build_variables.dr_tf_inboard
             )
-            build_variables.hpfdif = 0.0e0
+            build_variables.dz_tf_upper_lower_midplane = 0.0e0
         else:
             build_variables.z_tf_top = (
                 build_variables.dr_tf_inboard
@@ -842,10 +842,9 @@ def calculate_vertical_build(self, output: bool) -> None:
                 + build_variables.dz_fw_plasma_gap
                 + build_variables.z_plasma_xpoint_upper
             )
-            build_variables.hpfdif = (
-                build_variables.z_tf_top
-                - (build_variables.z_tf_inside_half + build_variables.dr_tf_inboard)
-            ) / 2.0e0
+            build_variables.dz_tf_upper_lower_midplane = build_variables.z_tf_top - (
+                build_variables.z_tf_inside_half + build_variables.dr_tf_inboard
+            )
 
     def divgeom(self, output: bool):
         """
diff --git a/process/data_structure/build_variables.py b/process/data_structure/build_variables.py
index 405908604f..6828986121 100644
--- a/process/data_structure/build_variables.py
+++ b/process/data_structure/build_variables.py
@@ -153,7 +153,7 @@
 """maximum (half-)height of TF coil (inside edge) (m)"""
 
 
-hpfdif: float = None
+dz_tf_upper_lower_midplane: float = None
 """difference in distance from midplane of upper and lower portions of TF
 legs (non-zero for single-null devices) (m)
 """
@@ -443,7 +443,7 @@ def init_build_variables():
     global gapomin
     global dr_shld_vv_gap_outboard
     global z_tf_inside_half
-    global hpfdif
+    global dz_tf_upper_lower_midplane
     global z_tf_top
     global hr1
     global iohcl
@@ -538,7 +538,7 @@ def init_build_variables():
     gapomin = 0.234
     dr_shld_vv_gap_outboard = 0.0
     z_tf_inside_half = 0.0
-    hpfdif = 0.0
+    dz_tf_upper_lower_midplane = 0.0
     z_tf_top = 0.0
     hr1 = 0.0
     iohcl = 1
diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
index 63692064e9..0289ec4563 100644
--- a/process/data_structure/pfcoil_variables.py
+++ b/process/data_structure/pfcoil_variables.py
@@ -43,9 +43,11 @@
 
 xind: list[float] = None
 
-rcls: list[float] = None
+r_pf_coil_middle_group_array: list[float] = None
+"""2D array of PF coil middle radii, indexed by group and coil in group"""
 
-zcls: list[float] = None
+z_pf_coil_middle_group_array: list[float] = None
+"""2D array of PF coil middle heights, indexed by group and coil in group"""
 
 ccls: list[float] = None
 
@@ -245,9 +247,8 @@
 """
 
 
-i_sup_pf_shape: int = None
+i_r_pf_outside_tf_placement: int = None
 """Switch for the placement of Location 3 (outboard) PF coils
-when the TF coils are superconducting (i_tf_sup = 1)
 - =0 (Default) Outboard PF coils follow TF shape
 in an ellipsoidal winding surface
 - =1 Outboard PF coils all have same radius, cylindrical
@@ -379,7 +380,7 @@
 """Full height of the central solenoid (m)"""
 
 
-routr: float = None
+dr_pf_tf_outboard_out_offset: float = None
 """radial distance (m) from outboard TF coil leg to centre of `i_pf_location=3` PF coils"""
 
 
@@ -387,7 +388,7 @@
 """radius of PF coil i (m)"""
 
 
-rpf1: float = None
+dr_pf_cs_middle_offset: float = None
 """offset (m) of radial position of `i_pf_location=1` PF coils from being directly above
 the central solenoid
 """
@@ -571,8 +572,8 @@ def init_pfcoil_module():
     global zfxf
     global cfxf
     global xind
-    global rcls
-    global zcls
+    global r_pf_coil_middle_group_array
+    global z_pf_coil_middle_group_array
     global ccls
     global ccl0
     global bpf2
@@ -592,8 +593,8 @@ def init_pfcoil_module():
     zfxf = np.zeros(NFIXMX)
     cfxf = np.zeros(NFIXMX)
     xind = np.zeros(NFIXMX)
-    rcls = np.zeros((N_PF_GROUPS_MAX, N_PF_COILS_IN_GROUP_MAX))
-    zcls = np.zeros((N_PF_GROUPS_MAX, N_PF_COILS_IN_GROUP_MAX))
+    r_pf_coil_middle_group_array = np.zeros((N_PF_GROUPS_MAX, N_PF_COILS_IN_GROUP_MAX))
+    z_pf_coil_middle_group_array = np.zeros((N_PF_GROUPS_MAX, N_PF_COILS_IN_GROUP_MAX))
     ccls = np.zeros(N_PF_GROUPS_MAX)
     ccl0 = np.zeros(N_PF_GROUPS_MAX)
     bpf2 = np.zeros(NGC2)
@@ -634,7 +635,7 @@ def init_pfcoil_variables():
     global j_crit_str_cs
     global j_crit_str_pf
     global i_pf_current
-    global i_sup_pf_shape
+    global i_r_pf_outside_tf_placement
     global j_cs_conductor_critical_pulse_start
     global j_cs_conductor_critical_flat_top_end
     global jcableoh_bop
@@ -664,9 +665,9 @@ def init_pfcoil_variables():
     global j_pf_wp_critical
     global r_cs_middle
     global dz_cs_full
-    global routr
+    global dr_pf_tf_outboard_out_offset
     global r_pf_coil_middle
-    global rpf1
+    global dr_pf_cs_middle_offset
     global rpf2
     global rref
     global s_shear_cs_peak
@@ -736,7 +737,7 @@ def init_pfcoil_variables():
     j_crit_str_cs = 0.0
     j_crit_str_pf = 0.0
     i_pf_current = 1
-    i_sup_pf_shape = 0
+    i_r_pf_outside_tf_placement = 0
     j_cs_conductor_critical_pulse_start = 0.0
     j_cs_conductor_critical_flat_top_end = 0.0
     jcableoh_bop = 0.0
@@ -766,9 +767,9 @@ def init_pfcoil_variables():
     j_pf_wp_critical = np.zeros(NGC2)
     r_cs_middle = 0.0
     dz_cs_full = 0.0
-    routr = 1.5
+    dr_pf_tf_outboard_out_offset = 1.5
     r_pf_coil_middle = np.zeros(NGC2)
-    rpf1 = 0.0
+    dr_pf_cs_middle_offset = 0.0
     rpf2 = -1.63
     rref = np.full(N_PF_GROUPS_MAX, 7.0)
     s_shear_cs_peak = 0.0
diff --git a/process/input.py b/process/input.py
index e38ec51f85..c16ec00578 100644
--- a/process/input.py
+++ b/process/input.py
@@ -1490,9 +1490,13 @@ def __post_init__(self):
     "roughness_fw_channel": InputVariable(
         data_structure.fwbs_variables, float, range=(0.0, 0.01)
     ),
-    "routr": InputVariable(data_structure.pfcoil_variables, float, range=(-3.0, 3.0)),
+    "dr_pf_tf_outboard_out_offset": InputVariable(
+        data_structure.pfcoil_variables, float, range=(-3.0, 3.0)
+    ),
     "row": InputVariable(data_structure.buildings_variables, float, range=(0.0, 10.0)),
-    "rpf1": InputVariable(data_structure.pfcoil_variables, float, range=(0.0, 3.0)),
+    "dr_pf_cs_middle_offset": InputVariable(
+        data_structure.pfcoil_variables, float, range=(0.0, 3.0)
+    ),
     "rpf2": InputVariable(data_structure.pfcoil_variables, float, range=(-3.0, 3.0)),
     "rrin": InputVariable(data_structure.ife_variables, float, range=(0.1, 50.0)),
     "rrmax": InputVariable(data_structure.ife_variables, float, range=(1.0, 50.0)),
@@ -2057,7 +2061,7 @@ def __post_init__(self):
         data_structure.physics_variables, int, choices=[0, 1]
     ),
     "i_str_wp": InputVariable(data_structure.tfcoil_variables, int, choices=[0, 1]),
-    "i_sup_pf_shape": InputVariable(
+    "i_r_pf_outside_tf_placement": InputVariable(
         data_structure.pfcoil_variables, int, choices=[0, 1]
     ),
     "i_tf_bucking": InputVariable(data_structure.tfcoil_variables, int, range=(0, 3)),
diff --git a/process/io/obsolete_vars.py b/process/io/obsolete_vars.py
index f597796750..d661fd9738 100644
--- a/process/io/obsolete_vars.py
+++ b/process/io/obsolete_vars.py
@@ -425,6 +425,9 @@
     "tdmptf": "t_tf_superconductor_quench",
     "tmaxpro": "temp_tf_conductor_quench_max",
     "coreradiationfraction": "f_p_plasma_core_rad_reduction",
+    "routr": "dr_pf_tf_outboard_out_offset",
+    "rpf1": "dr_pf_cs_middle_offset",
+    "i_sup_pf_shape": "i_r_pf_outside_tf_placement",
 }
 
 OBS_VARS_HELP = {
diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
index f1fdb474b0..c5baaaad60 100644
--- a/process/io/plot_proc.py
+++ b/process/io/plot_proc.py
@@ -8764,6 +8764,9 @@ def plot_tf_coil_structure(axis, mfile_data, scan, colour_scheme=1):
     r_tf_inboard_in = mfile_data.data["r_tf_inboard_in"].get_scan(scan)
     dr_tf_outboard = mfile_data.data["dr_tf_outboard"].get_scan(scan)
     len_tf_coil = mfile_data.data["len_tf_coil"].get_scan(scan)
+    dz_tf_upper_lower_midplane = mfile_data.data["dz_tf_upper_lower_midplane"].get_scan(
+        scan
+    )
 
     # Plot the points as black dots, number them, and connect them with lines
     xs = [x1, x2, x3, x4, x5]
@@ -9036,7 +9039,7 @@ def plot_tf_coil_structure(axis, mfile_data, scan, colour_scheme=1):
 
     # ==============================================================
 
-    # Add a label for the inboard thickness
+    # Add a label for the length of the coil
     axis.text(
         (r_tf_outboard_in + 2 * dr_tf_outboard),
         0.0,
@@ -9049,6 +9052,19 @@ def plot_tf_coil_structure(axis, mfile_data, scan, colour_scheme=1):
 
     # ==============================================================
 
+    # Add a label for the length of the coil
+    axis.text(
+        (r_tf_outboard_in + 2 * dr_tf_outboard),
+        -1.0,
+        f"$\\Delta Z$ upper and lower to midplane = {dz_tf_upper_lower_midplane:.3f} m",
+        fontsize=7,
+        color="black",
+        verticalalignment="center",
+        bbox={"boxstyle": "round", "facecolor": "pink", "alpha": 1.0},
+    )
+
+    # ==============================================================
+
     # Add arow for inboard coil radius
     axis.annotate(
         "",
diff --git a/process/pfcoil.py b/process/pfcoil.py
index 03627f0b7f..935f7aa5c5 100644
--- a/process/pfcoil.py
+++ b/process/pfcoil.py
@@ -14,7 +14,11 @@
 from process.data_structure import constraint_variables as ctv
 from process.data_structure import cs_fatigue_variables as csfv
 from process.data_structure import fwbs_variables as fwbsv
-from process.data_structure import numerics, pfcoil_variables
+from process.data_structure import (
+    numerics,
+    pfcoil_variables,
+    superconducting_tf_coil_variables,
+)
 from process.data_structure import physics_variables as pv
 from process.data_structure import rebco_variables as rcv
 from process.data_structure import tfcoil_variables as tfv
@@ -212,97 +216,125 @@ def pfcoil(self):
         # Scale PF coil locations
         signn[0] = 1.0e0
         signn[1] = -1.0e0
-        pfcoil_variables.rclsnorm = (
-            bv.r_tf_outboard_mid + 0.5e0 * bv.dr_tf_outboard + pfcoil_variables.routr
+        pfcoil_variables.r_pf_outside_tf_midplane = (
+            superconducting_tf_coil_variables.r_tf_outboard_out
+            + pfcoil_variables.dr_pf_tf_outboard_out_offset
         )
 
         # Place the PF coils:
 
-        # N.B. Problems here if k=n_pf_coils_in_group(group) is greater than 2.
-        for j in range(pfcoil_variables.n_pf_coil_groups):
-            if pfcoil_variables.i_pf_location[j] == 1:
+        # N.B. Problems here if coil=n_pf_coils_in_group(group) is greater than 2.
+        for group in range(pfcoil_variables.n_pf_coil_groups):
+            if pfcoil_variables.i_pf_location[group] == 1:
                 # PF coil is stacked on top of the Central Solenoid
-                for k in range(pfcoil_variables.n_pf_coils_in_group[j]):
-                    pfcoil_variables.rcls[j, k] = (
-                        pfcoil_variables.r_cs_middle + pfcoil_variables.rpf1
+                # Use a helper function to compute r_pf_coil_middle_group_array and
+                # z_pf_coil_middle_group_array arrays for this group
+
+                r_pf_coil_middle_group_array, z_pf_coil_middle_group_array = (
+                    self.place_pf_above_cs(
+                        n_pf_coils_in_group=pfcoil_variables.n_pf_coils_in_group,
+                        n_pf_group=group,
+                        r_cs_middle=pfcoil_variables.r_cs_middle,
+                        dr_pf_cs_middle_offset=pfcoil_variables.dr_pf_cs_middle_offset,
+                        z_tf_inside_half=bv.z_tf_inside_half,
+                        dr_tf_inboard=bv.dr_tf_inboard,
+                        z_cs_coil_upper=pfcoil_variables.dz_cs_full / 2,
                     )
-
-                    # Z coordinate of coil enforced so as not
-                    # to occupy the same space as the Central Solenoid
-                    pfcoil_variables.zcls[j, k] = signn[k] * (
-                        bv.z_tf_inside_half * pfcoil_variables.f_z_cs_tf_internal
-                        + 0.1e0
-                        + 0.5e0
-                        * (
-                            bv.z_tf_inside_half
-                            * (1.0e0 - pfcoil_variables.f_z_cs_tf_internal)
-                            + bv.dr_tf_inboard
-                            + 0.1e0
-                        )
+                )
+                for coil in range(pfcoil_variables.n_pf_coils_in_group[group]):
+                    pfcoil_variables.r_pf_coil_middle_group_array[group, coil] = (
+                        r_pf_coil_middle_group_array[group, coil]
                     )
+                    pfcoil_variables.z_pf_coil_middle_group_array[group, coil] = (
+                        z_pf_coil_middle_group_array[group, coil]
+                    )
+
+            # =========================================================================
 
-            elif pfcoil_variables.i_pf_location[j] == 2:
+            elif pfcoil_variables.i_pf_location[group] == 2:
                 # PF coil is on top of the TF coil
-                for k in range(pfcoil_variables.n_pf_coils_in_group[j]):
-                    pfcoil_variables.rcls[j, k] = (
-                        pv.rmajor + pfcoil_variables.rpf2 * pv.triang * pv.rminor
-                    )
-                    if pv.itart == 1 and pv.itartpf == 0:
-                        pfcoil_variables.zcls[j, k] = (
-                            bv.z_tf_inside_half - pfcoil_variables.zref[j]
-                        ) * signn[k]
-                    else:
-                        # pfcoil_variables.zcls(j,k) = (bv.z_tf_inside_half + bv.dr_tf_inboard + 0.86e0) * signn(k)
-                        if top_bottom == 1:  # this coil is above midplane
-                            pfcoil_variables.zcls[j, k] = bv.z_tf_top + 0.86e0
-                            top_bottom = -1
-                        else:  # this coil is below midplane
-                            pfcoil_variables.zcls[j, k] = -1.0e0 * (
-                                bv.z_tf_top - 2.0e0 * bv.hpfdif + 0.86e0
-                            )
-                            top_bottom = 1
+                (
+                    r_pf_coil_middle_group_array,
+                    z_pf_coil_middle_group_array,
+                    top_bottom,
+                ) = self.place_pf_above_tf(
+                    n_pf_coils_in_group=pfcoil_variables.n_pf_coils_in_group,
+                    n_pf_group=group,
+                    rmajor=pv.rmajor,
+                    triang=pv.triang,
+                    rminor=pv.rminor,
+                    itart=pv.itart,
+                    itartpf=pv.itartpf,
+                    z_tf_inside_half=bv.z_tf_inside_half,
+                    dz_tf_upper_lower_midplane=bv.dz_tf_upper_lower_midplane,
+                    z_tf_top=bv.z_tf_top,
+                    top_bottom=top_bottom,
+                    rpf2=pfcoil_variables.rpf2,
+                    zref=pfcoil_variables.zref,
+                )
 
-            elif pfcoil_variables.i_pf_location[j] == 3:
+                for coil in range(pfcoil_variables.n_pf_coils_in_group[group]):
+                    pfcoil_variables.r_pf_coil_middle_group_array[group, coil] = (
+                        r_pf_coil_middle_group_array[group, coil]
+                    )
+                    pfcoil_variables.z_pf_coil_middle_group_array[group, coil] = (
+                        z_pf_coil_middle_group_array[group, coil]
+                    )
+            # =========================================================================
+
+            elif pfcoil_variables.i_pf_location[group] == 3:
                 # PF coil is radially outside the TF coil
-                for k in range(pfcoil_variables.n_pf_coils_in_group[j]):
-                    pfcoil_variables.zcls[j, k] = (
-                        pv.rminor * pfcoil_variables.zref[j] * signn[k]
-                    )
-                    # Coil radius follows TF coil curve for SC TF (D-shape)
-                    # otherwise stacked for resistive TF (rectangle-shape)
-                    if tfv.i_tf_sup != 1 or pfcoil_variables.i_sup_pf_shape == 1:
-                        pfcoil_variables.rcls[j, k] = pfcoil_variables.rclsnorm
-                    else:
-                        pfcoil_variables.rcls[j, k] = math.sqrt(
-                            pfcoil_variables.rclsnorm**2
-                            - pfcoil_variables.zcls[j, k] ** 2
-                        )
-                        try:
-                            assert pfcoil_variables.rcls[j, k] < np.inf
-                        except AssertionError:
-                            logger.exception(
-                                "Element of pfcoil_variables.rcls is inf. Kludging to 1e10."
-                            )
-                            pfcoil_variables.rcls[j, k] = 1e10
+                (
+                    r_pf_coil_middle_group_array,
+                    z_pf_coil_middle_group_array,
+                ) = self.place_pf_outside_tf(
+                    n_pf_coils_in_group=pfcoil_variables.n_pf_coils_in_group,
+                    n_pf_group=group,
+                    rminor=pv.rminor,
+                    zref=pfcoil_variables.zref,
+                    i_tf_shape=tfv.i_tf_shape,
+                    i_r_pf_outside_tf_placement=pfcoil_variables.i_r_pf_outside_tf_placement,
+                    r_pf_outside_tf_midplane=pfcoil_variables.r_pf_outside_tf_midplane,
+                )
+
+                for coil in range(pfcoil_variables.n_pf_coils_in_group[group]):
+                    pfcoil_variables.r_pf_coil_middle_group_array[group, coil] = (
+                        r_pf_coil_middle_group_array[group, coil]
+                    )
+                    pfcoil_variables.z_pf_coil_middle_group_array[group, coil] = (
+                        z_pf_coil_middle_group_array[group, coil]
+                    )
+
+            # =========================================================================
 
-            elif pfcoil_variables.i_pf_location[j] == 4:
-                # PF coil is in general location
-                # See issue 1418
-                # https://git.ccfe.ac.uk/process/process/-/issues/1418
-                for k in range(pfcoil_variables.n_pf_coils_in_group[j]):
-                    pfcoil_variables.zcls[j, k] = (
-                        pv.rminor * pfcoil_variables.zref[j] * signn[k]
+            elif pfcoil_variables.i_pf_location[group] == 4:
+                (
+                    r_pf_coil_middle_group_array,
+                    z_pf_coil_middle_group_array,
+                ) = self.place_pf_generally(
+                    n_pf_coils_in_group=pfcoil_variables.n_pf_coils_in_group,
+                    n_pf_group=group,
+                    rminor=pv.rminor,
+                    rmajor=pv.rmajor,
+                    zref=pfcoil_variables.zref,
+                    rref=pfcoil_variables.rref,
+                )
+
+                for coil in range(pfcoil_variables.n_pf_coils_in_group[group]):
+                    pfcoil_variables.r_pf_coil_middle_group_array[group, coil] = (
+                        r_pf_coil_middle_group_array[group, coil]
                     )
-                    pfcoil_variables.rcls[j, k] = (
-                        pv.rminor * pfcoil_variables.rref[j] + pv.rmajor
+                    pfcoil_variables.z_pf_coil_middle_group_array[group, coil] = (
+                        z_pf_coil_middle_group_array[group, coil]
                     )
 
             else:
                 raise ProcessValueError(
                     "Illegal i_pf_location value",
-                    j=j,
-                    i_pf_location=pfcoil_variables.i_pf_location[j],
+                    group=group,
+                    i_pf_location=pfcoil_variables.i_pf_location[group],
                 )
+            # =========================================================================
 
         # Allocate current to the PF coils:
         # "Flux swing coils" participate in cancellation of the CS
@@ -347,8 +379,8 @@ def pfcoil(self):
                 pfcoil_variables.cfxf,
                 pfcoil_variables.n_pf_coil_groups,
                 pfcoil_variables.n_pf_coils_in_group,
-                pfcoil_variables.rcls,
-                pfcoil_variables.zcls,
+                pfcoil_variables.r_pf_coil_middle_group_array,
+                pfcoil_variables.z_pf_coil_middle_group_array,
                 pfcoil_variables.alfapf,
                 bfix,
                 gmat,
@@ -416,12 +448,12 @@ def pfcoil(self):
                         pfcoil_variables.ccls[i] = 0.0e0
                         nfxf0 = nfxf0 + pfcoil_variables.n_pf_coils_in_group[i]
                         for ccount in range(pfcoil_variables.n_pf_coils_in_group[i]):
-                            pfcoil_variables.rfxf[nocoil] = pfcoil_variables.rcls[
-                                i, ccount
-                            ]
-                            pfcoil_variables.zfxf[nocoil] = pfcoil_variables.zcls[
-                                i, ccount
-                            ]
+                            pfcoil_variables.rfxf[nocoil] = (
+                                pfcoil_variables.r_pf_coil_middle_group_array[i, ccount]
+                            )
+                            pfcoil_variables.zfxf[nocoil] = (
+                                pfcoil_variables.z_pf_coil_middle_group_array[i, ccount]
+                            )
                             pfcoil_variables.cfxf[nocoil] = pfcoil_variables.ccls[i]
                             nocoil = nocoil + 1
 
@@ -435,17 +467,19 @@ def pfcoil(self):
                             * (
                                 1.0e0
                                 - (pv.kappa * pv.rminor)
-                                / abs(pfcoil_variables.zcls[i, 0])
+                                / abs(
+                                    pfcoil_variables.z_pf_coil_middle_group_array[i, 0]
+                                )
                             )
                         )
                         nfxf0 = nfxf0 + pfcoil_variables.n_pf_coils_in_group[i]
                         for ccount in range(pfcoil_variables.n_pf_coils_in_group[i]):
-                            pfcoil_variables.rfxf[nocoil] = pfcoil_variables.rcls[
-                                i, ccount
-                            ]
-                            pfcoil_variables.zfxf[nocoil] = pfcoil_variables.zcls[
-                                i, ccount
-                            ]
+                            pfcoil_variables.rfxf[nocoil] = (
+                                pfcoil_variables.r_pf_coil_middle_group_array[i, ccount]
+                            )
+                            pfcoil_variables.zfxf[nocoil] = (
+                                pfcoil_variables.z_pf_coil_middle_group_array[i, ccount]
+                            )
                             pfcoil_variables.cfxf[nocoil] = pfcoil_variables.ccls[i]
                             nocoil = nocoil + 1
 
@@ -474,18 +508,26 @@ def pfcoil(self):
 
                 for ccount in range(ngrp0):
                     ncls0[ccount] = 2
-                    pfcoil_variables.rcls0[ccount, 0] = pfcoil_variables.rcls[
-                        pcls0[ccount] - 1, 0
-                    ]
-                    pfcoil_variables.rcls0[ccount, 1] = pfcoil_variables.rcls[
-                        pcls0[ccount] - 1, 1
-                    ]
-                    pfcoil_variables.zcls0[ccount, 0] = pfcoil_variables.zcls[
-                        pcls0[ccount] - 1, 0
-                    ]
-                    pfcoil_variables.zcls0[ccount, 1] = pfcoil_variables.zcls[
-                        pcls0[ccount] - 1, 1
-                    ]
+                    pfcoil_variables.rcls0[ccount, 0] = (
+                        pfcoil_variables.r_pf_coil_middle_group_array[
+                            pcls0[ccount] - 1, 0
+                        ]
+                    )
+                    pfcoil_variables.rcls0[ccount, 1] = (
+                        pfcoil_variables.r_pf_coil_middle_group_array[
+                            pcls0[ccount] - 1, 1
+                        ]
+                    )
+                    pfcoil_variables.zcls0[ccount, 0] = (
+                        pfcoil_variables.z_pf_coil_middle_group_array[
+                            pcls0[ccount] - 1, 0
+                        ]
+                    )
+                    pfcoil_variables.zcls0[ccount, 1] = (
+                        pfcoil_variables.z_pf_coil_middle_group_array[
+                            pcls0[ccount] - 1, 1
+                        ]
+                    )
 
                 npts0 = 1
                 rpts[0] = pv.rmajor
@@ -596,8 +638,12 @@ def pfcoil(self):
         ncl = 0
         for nng in range(pfcoil_variables.n_pf_coil_groups):
             for ng2 in range(pfcoil_variables.n_pf_coils_in_group[nng]):
-                pfcoil_variables.r_pf_coil_middle[ncl] = pfcoil_variables.rcls[nng, ng2]
-                pfcoil_variables.z_pf_coil_middle[ncl] = pfcoil_variables.zcls[nng, ng2]
+                pfcoil_variables.r_pf_coil_middle[ncl] = (
+                    pfcoil_variables.r_pf_coil_middle_group_array[nng, ng2]
+                )
+                pfcoil_variables.z_pf_coil_middle[ncl] = (
+                    pfcoil_variables.z_pf_coil_middle_group_array[nng, ng2]
+                )
 
                 # Currents at different times:
 
@@ -1023,6 +1069,274 @@ def pfcoil(self):
             pfcoil_variables.n_pf_cs_plasma_circuits - 1, 5
         ] = 0.0e0
 
+    def place_pf_above_cs(
+        self,
+        n_pf_coils_in_group: np.ndarray,
+        n_pf_group: int,
+        r_cs_middle: float,
+        dr_pf_cs_middle_offset: float,
+        z_tf_inside_half: float,
+        dr_tf_inboard: float,
+        z_cs_coil_upper: float,
+    ) -> tuple[np.ndarray, np.ndarray]:
+        """
+        Calculate the placement of PF coils stacked above the Central Solenoid.
+
+        :param n_pf_coils_in_group: Array containing the number of coils in each PF group.
+        :type n_pf_coils_in_group: np.ndarray
+        :param n_pf_group: Index of the PF coil group.
+        :type n_pf_group: int
+        :param r_cs_middle: Radial coordinate of CS coil centre (m).
+        :type r_cs_middle: float
+        :param dr_pf_cs_middle_offset: Radial offset for PF coil placement (m).
+        :type dr_pf_cs_middle_offset: float
+        :param z_tf_inside_half: Half-height of the TF bore (m).
+        :type z_tf_inside_half: float
+        :param dr_tf_inboard: Thickness of the TF inboard leg (m).
+        :type dr_tf_inboard: float
+        :param z_cs_coil_upper: Upper z coordinate of the CS coil (m).
+        :type z_cs_coil_upper: float
+        :return: Tuple of arrays containing the radial and vertical coordinates of PF coils in the group.
+        :rtype: tuple[np.ndarray, np.ndarray]
+        """
+
+        # Initialise as empty arrays; will be resized in the loop
+        r_pf_coil_middle_group_array = np.zeros((
+            pfcoil_variables.n_pf_coil_groups,
+            pfcoil_variables.N_PF_COILS_IN_GROUP_MAX,
+        ))
+        z_pf_coil_middle_group_array = np.zeros((
+            pfcoil_variables.n_pf_coil_groups,
+            pfcoil_variables.N_PF_COILS_IN_GROUP_MAX,
+        ))
+
+        for coil in range(n_pf_coils_in_group[n_pf_group]):
+            # Positions PF coil directly above centre of CS with offset from dr_pf_cs_middle_offset
+            r_pf_coil_middle_group_array[n_pf_group, coil] = (
+                r_cs_middle + dr_pf_cs_middle_offset
+            )
+
+            # Z coordinate of coil enforced so as not
+            # to occupy the same space as the Central Solenoid
+            # Set sign: +1 for coil 0, -1 for coil 1
+            sign = 1.0 if coil == 0 else -1.0
+            z_pf_coil_middle_group_array[n_pf_group, coil] = sign * (
+                z_cs_coil_upper
+                + 0.1e0
+                + 0.5e0 * ((z_tf_inside_half - z_cs_coil_upper) + dr_tf_inboard + 0.1e0)
+            )
+        return r_pf_coil_middle_group_array, z_pf_coil_middle_group_array
+
+    def place_pf_above_tf(
+        self,
+        n_pf_coils_in_group: np.ndarray,
+        n_pf_group: int,
+        rmajor: float,
+        triang: float,
+        rminor: float,
+        itart: int,
+        itartpf: int,
+        z_tf_inside_half: float,
+        dz_tf_upper_lower_midplane: float,
+        z_tf_top: float,
+        top_bottom: int,
+        rpf2: float,
+        zref: np.ndarray,
+    ) -> tuple[np.ndarray, np.ndarray, int]:
+        """
+        Calculates and places poloidal field (PF) coils above the toroidal field (TF) coils for a given group.
+
+        :param n_pf_coils_in_group: Array containing the number of PF coils in each group.
+        :type n_pf_coils_in_group: np.ndarray
+        :param n_pf_group: Index of the PF coil group to process.
+        :type n_pf_group: int
+        :param rmajor: Major radius of the device.
+        :type rmajor: float
+        :param triang: Triangularity parameter for coil placement.
+        :type triang: float
+        :param rminor: Minor radius of the device.
+        :type rminor: float
+        :param itart: Flag indicating ST configuration.
+        :type itart: int
+        :param itartpf: Flag indicating PF coil configuration for ST.
+        :type itartpf: int
+        :param z_tf_inside_half: Half-height of the TF coil inside region.
+        :type z_tf_inside_half: float
+        :param dz_tf_upper_lower_midplane: Height difference parameter for PF coil placement.
+        :type dz_tf_upper_lower_midplane: float
+        :param z_tf_top: Top z-coordinate of the TF coil.
+        :type z_tf_top: float
+        :param top_bottom: Indicator for coil placement above (+1) or below (-1) the midplane.
+        :type top_bottom: int
+        :param rpf2: Radial offset parameter for PF coil placement.
+        :type rpf2: float
+        :param zref: Array of reference z-coordinates for PF coil placement.
+        :type zref: np.ndarray
+
+        :returns: Tuple containing arrays of radial and vertical positions of PF coil middles for the specified group,
+                  and the updated top_bottom indicator.
+        :rtype: tuple[np.ndarray, np.ndarray, int]
+        """
+        # Initialise as empty arrays; will be resized in the loop
+        r_pf_coil_middle_group_array = np.zeros((
+            pfcoil_variables.n_pf_coil_groups,
+            pfcoil_variables.N_PF_COILS_IN_GROUP_MAX,
+        ))
+        z_pf_coil_middle_group_array = np.zeros((
+            pfcoil_variables.n_pf_coil_groups,
+            pfcoil_variables.N_PF_COILS_IN_GROUP_MAX,
+        ))
+
+        for coil in range(n_pf_coils_in_group[n_pf_group]):
+            # Place PF coils at radius determined by rmajor, triang and rminor
+            r_pf_coil_middle_group_array[n_pf_group, coil] = (
+                rmajor + rpf2 * triang * rminor
+            )
+
+            # Set sign: +1 for coil 0, -1 for coil 1
+            sign = 1.0 if coil == 0 else -1.0
+            if itart == 1 and itartpf == 0:
+                z_pf_coil_middle_group_array[n_pf_group, coil] = (
+                    z_tf_inside_half - zref[n_pf_group]
+                ) * sign
+            else:
+                if top_bottom == 1:  # this coil is above midplane
+                    z_pf_coil_middle_group_array[n_pf_group, coil] = z_tf_top + 0.86e0
+                    top_bottom = -1
+                else:  # this coil is below midplane
+                    z_pf_coil_middle_group_array[n_pf_group, coil] = -1.0e0 * (
+                        z_tf_top - dz_tf_upper_lower_midplane + 0.86e0
+                    )
+                    top_bottom = 1
+
+        return r_pf_coil_middle_group_array, z_pf_coil_middle_group_array, top_bottom
+
+    def place_pf_outside_tf(
+        self,
+        n_pf_coils_in_group: np.ndarray,
+        n_pf_group: int,
+        rminor: float,
+        zref: np.ndarray,
+        i_tf_shape: int,
+        i_r_pf_outside_tf_placement: int,
+        r_pf_outside_tf_midplane: float,
+    ) -> tuple[np.ndarray, np.ndarray]:
+        """
+        Calculates the radial and vertical positions of poloidal field (PF) coils placed outside the toroidal field (TF) coil.
+
+        :param n_pf_coils_in_group: Array containing the number of PF coils in each group.
+        :type n_pf_coils_in_group: np.ndarray
+        :param n_pf_group: Index of the PF coil group to process.
+        :type n_pf_group: int
+        :param rminor: Minor radius of the device.
+        :type rminor: float
+        :param zref: Reference vertical positions for each PF coil group.
+        :type zref: np.ndarray
+        :param i_tf_shape: Integer flag indicating TF coil shape (2 for picture frame, others for D-shape).
+        :type i_tf_shape: int
+        :param i_r_pf_outside_tf_placement: Placement switch for PF coil radius (1 for constant/stacked, 0 for following TF curve).
+        :type i_r_pf_outside_tf_placement: int
+        :param r_pf_outside_tf_midplane: Radial position of PF coil at the midplane.
+        :type r_pf_outside_tf_midplane: float
+
+        :returns: Tuple containing arrays of radial and vertical positions of PF coil centers for the specified group.
+        :rtype: tuple[np.ndarray, np.ndarray]
+        """
+
+        # Initialise as empty arrays; will be resized in the loop
+        r_pf_coil_middle_group_array = np.zeros((
+            pfcoil_variables.n_pf_coil_groups,
+            pfcoil_variables.N_PF_COILS_IN_GROUP_MAX,
+        ))
+        z_pf_coil_middle_group_array = np.zeros((
+            pfcoil_variables.n_pf_coil_groups,
+            pfcoil_variables.N_PF_COILS_IN_GROUP_MAX,
+        ))
+
+        # PF coil is radially outside the TF coil
+
+        for coil in range(n_pf_coils_in_group[n_pf_group]):
+            sign = 1.0 if coil == 0 else -1.0
+
+            z_pf_coil_middle_group_array[n_pf_group, coil] = (
+                rminor * zref[n_pf_group] * sign
+            )
+            # Coil radius is constant / stacked for picture frame TF or if placement switch is set
+            if i_tf_shape == 2 or i_r_pf_outside_tf_placement == 1:
+                r_pf_coil_middle_group_array[n_pf_group, coil] = (
+                    r_pf_outside_tf_midplane
+                )
+            else:
+                # Coil radius follows TF coil curve for TF (D-shape)
+                r_pf_coil_middle_group_array[n_pf_group, coil] = math.sqrt(
+                    r_pf_outside_tf_midplane**2
+                    - z_pf_coil_middle_group_array[n_pf_group, coil] ** 2
+                )
+                try:
+                    assert r_pf_coil_middle_group_array[n_pf_group, coil] < np.inf
+                except AssertionError:
+                    logger.exception(
+                        "Element of pfcoil_variables.r_pf_coil_middle_group_array is inf. Kludging to 1e10."
+                    )
+                    r_pf_coil_middle_group_array[n_pf_group, coil] = 1e10
+        return (
+            r_pf_coil_middle_group_array,
+            z_pf_coil_middle_group_array,
+        )
+
+    def place_pf_generally(
+        self,
+        n_pf_coils_in_group: np.ndarray,
+        n_pf_group: int,
+        rminor: float,
+        rmajor: float,
+        zref: np.ndarray,
+        rref: np.ndarray,
+    ) -> tuple[np.ndarray, np.ndarray]:
+        """
+        Calculates the radial and vertical positions of poloidal field (PF) coils placed in a general location.
+
+        :param n_pf_coils_in_group: Array containing the number of PF coils in each group.
+        :type n_pf_coils_in_group: numpy.ndarray
+        :param n_pf_group: Index of the PF coil group to process.
+        :type n_pf_group: int
+        :param rminor: Minor radius of the device.
+        :type rminor: float
+        :param rmajor: Major radius of the device.
+        :type rmajor: float
+        :param zref: Reference vertical positions for each PF coil group.
+        :type zref: numpy.ndarray
+        :param rref: Reference radial positions for each PF coil group.
+        :type rref: numpy.ndarray
+
+        :returns: Tuple containing arrays of radial and vertical positions of PF coil centers for the specified group.
+        :rtype: tuple[numpy.ndarray, numpy.ndarray]
+        """
+        r_pf_coil_middle_group_array: np.ndarray = np.zeros((
+            pfcoil_variables.n_pf_coil_groups,
+            pfcoil_variables.N_PF_COILS_IN_GROUP_MAX,
+        ))
+        z_pf_coil_middle_group_array: np.ndarray = np.zeros((
+            pfcoil_variables.n_pf_coil_groups,
+            pfcoil_variables.N_PF_COILS_IN_GROUP_MAX,
+        ))
+
+        for coil in range(n_pf_coils_in_group[n_pf_group]):
+            sign: float = 1.0 if coil == 0 else -1.0
+
+            # Place as mutiples of minor radius from the midplane
+            z_pf_coil_middle_group_array[n_pf_group, coil] = (
+                rminor * zref[n_pf_group] * sign
+            )
+            # Place as multiples of minor radius from the plasma centre
+            r_pf_coil_middle_group_array[n_pf_group, coil] = (
+                rminor * rref[n_pf_group] + rmajor
+            )
+        return (
+            r_pf_coil_middle_group_array,
+            z_pf_coil_middle_group_array,
+        )
+
     def efc(
         self,
         npts,
@@ -1036,8 +1350,8 @@ def efc(
         cfix,
         n_pf_coil_groups,
         n_pf_coils_in_group,
-        rcls,
-        zcls,
+        r_pf_coil_middle_group_array,
+        z_pf_coil_middle_group_array,
         alfa,
         bfix,
         gmat,
@@ -1078,10 +1392,10 @@ def efc(
         :type n_pf_coil_groups: int
         :param n_pf_coils_in_group: number of coils in each group, each value <= n_pf_coils_in_group_max
         :type n_pf_coils_in_group: np.ndarray
-        :param rcls: coords R(i,j), Z(i,j) of coil j in group i (m)
-        :type rcls: np.ndarray
-        :param zcls: coords R(i,j), Z(i,j) of coil j in group i (m)
-        :type zcls: np.ndarray
+        :param r_pf_coil_middle_group_array: coords R(i,j), Z(i,j) of coil j in group i (m)
+        :type r_pf_coil_middle_group_array: np.ndarray
+        :param z_pf_coil_middle_group_array: coords R(i,j), Z(i,j) of coil j in group i (m)
+        :type z_pf_coil_middle_group_array: np.ndarray
         :param alfa: smoothing parameter (0 = no smoothing, 1.0D-9 = large
         smoothing)
         :type alfa: float
@@ -1126,8 +1440,8 @@ def efc(
             bzin,
             int(n_pf_coil_groups),
             n_pf_coils_in_group,
-            rcls,
-            zcls,
+            r_pf_coil_middle_group_array,
+            z_pf_coil_middle_group_array,
             alfa,
             bfix,
             int(pfcoil_variables.N_PF_COILS_IN_GROUP_MAX),
@@ -1156,19 +1470,19 @@ def tf_pf_collision_detector(self):
             for i in range(pfcoil_variables.n_pf_coil_groups):
                 for ii in range(pfcoil_variables.n_pf_coil_groups):
                     for ij in range(pfcoil_variables.n_pf_coils_in_group[ii]):
-                        if pfcoil_variables.rcls[
+                        if pfcoil_variables.r_pf_coil_middle_group_array[
                             ii, ij
                         ] <= (  # Outboard TF coil collision
-                            pfcoil_variables.rclsnorm
-                            - pfcoil_variables.routr
+                            pfcoil_variables.r_pf_outside_tf_midplane
+                            - pfcoil_variables.dr_pf_tf_outboard_out_offset
                             + pfcoil_variables.r_pf_coil_middle[i]
-                        ) and pfcoil_variables.rcls[ii, ij] >= (
+                        ) and pfcoil_variables.r_pf_coil_middle_group_array[ii, ij] >= (
                             bv.r_tf_outboard_mid
                             - (0.5 * bv.dr_tf_outboard)
                             - pfcoil_variables.r_pf_coil_middle[i]
                         ):
                             pf_tf_collision += 1
-                        if pfcoil_variables.rcls[
+                        if pfcoil_variables.r_pf_coil_middle_group_array[
                             ii, ij
                         ] <= (  # Inboard TF coil collision
                             bv.dr_bore
@@ -1177,7 +1491,7 @@ def tf_pf_collision_detector(self):
                             + bv.dr_cs_tf_gap
                             + bv.dr_tf_inboard
                             + pfcoil_variables.r_pf_coil_middle[i]
-                        ) and pfcoil_variables.rcls[ii, ij] >= (
+                        ) and pfcoil_variables.r_pf_coil_middle_group_array[ii, ij] >= (
                             bv.dr_bore
                             + bv.dr_cs
                             + bv.dr_cs_precomp
@@ -1186,9 +1500,11 @@ def tf_pf_collision_detector(self):
                         ):
                             pf_tf_collision += 1
                         if (  # Vertical TF coil collision
-                            abs(pfcoil_variables.zcls[ii, ij])
+                            abs(pfcoil_variables.z_pf_coil_middle_group_array[ii, ij])
                             <= bv.z_tf_top + pfcoil_variables.r_pf_coil_middle[i]
-                            and abs(pfcoil_variables.zcls[ii, ij])
+                            and abs(
+                                pfcoil_variables.z_pf_coil_middle_group_array[ii, ij]
+                            )
                             >= bv.z_tf_top
                             - (0.5 * bv.dr_tf_outboard)
                             - pfcoil_variables.r_pf_coil_middle[i]
@@ -3947,8 +4263,8 @@ def mtrx(
     bzin,
     n_pf_coil_groups,
     n_pf_coils_in_group,
-    rcls,
-    zcls,
+    r_pf_coil_middle_group_array,
+    z_pf_coil_middle_group_array,
     alfa,
     bfix,
     n_pf_coils_in_group_max,
@@ -3983,10 +4299,10 @@ def mtrx(
     :type n_pf_coil_groups: int
     :param n_pf_coils_in_group: number of coils in each group, each value <= n_pf_coils_in_group_max
     :type n_pf_coils_in_group: numpy.ndarray
-    :param rcls: coords R(i,j), Z(i,j) of coil j in group i (m)
-    :type rcls: numpy.ndarray
-    :param zcls: coords R(i,j), Z(i,j) of coil j in group i (m)
-    :type zcls: numpy.ndarray
+    :param r_pf_coil_middle_group_array: coords R(i,j), Z(i,j) of coil j in group i (m)
+    :type r_pf_coil_middle_group_array: numpy.ndarray
+    :param z_pf_coil_middle_group_array: coords R(i,j), Z(i,j) of coil j in group i (m)
+    :type z_pf_coil_middle_group_array: numpy.ndarray
     :param alfa: smoothing parameter (0 = no smoothing, 1.0D-9 = large
     smoothing)
     :type alfa: float
@@ -4010,7 +4326,11 @@ def mtrx(
             nc = n_pf_coils_in_group[j]
 
             _, gmat[i, j], gmat[i + npts, j], _ = bfield(
-                rcls[j, :nc], zcls[j, :nc], cc[:nc], rpts[i], zpts[i]
+                r_pf_coil_middle_group_array[j, :nc],
+                z_pf_coil_middle_group_array[j, :nc],
+                cc[:nc],
+                rpts[i],
+                zpts[i],
             )
 
     # Add constraint equations
diff --git a/process/tf_coil.py b/process/tf_coil.py
index 92cca87391..eaa684d5db 100644
--- a/process/tf_coil.py
+++ b/process/tf_coil.py
@@ -762,6 +762,13 @@ def outtf(self):
             build_variables.z_tf_top,
             "OP ",
         )
+        po.ovarre(
+            self.outfile,
+            "Height difference in upper and lower TF from midplane (m)",
+            "(dz_tf_upper_lower_midplane)",
+            build_variables.dz_tf_upper_lower_midplane,
+            "OP ",
+        )
         if physics_variables.itart == 1:
             po.ovarre(
                 self.outfile,
diff --git a/tests/integration/test_pfcoil_int.py b/tests/integration/test_pfcoil_int.py
index ca4cc5b883..755cc5d16c 100644
--- a/tests/integration/test_pfcoil_int.py
+++ b/tests/integration/test_pfcoil_int.py
@@ -17,7 +17,7 @@
 from process.cs_fatigue import CsFatigue
 from process.data_structure import build_variables as bv
 from process.data_structure import fwbs_variables as fwbsv
-from process.data_structure import pfcoil_variables
+from process.data_structure import pfcoil_variables, superconducting_tf_coil_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
@@ -46,7 +46,7 @@ def test_pfcoil(monkeypatch, pfcoil):
     """
 
     monkeypatch.setattr(bv, "iohcl", 1)
-    monkeypatch.setattr(bv, "hpfdif", 0.0)
+    monkeypatch.setattr(bv, "dz_tf_upper_lower_midplane", 0.0)
     monkeypatch.setattr(bv, "z_tf_top", 4.0)  # guess
     monkeypatch.setattr(bv, "z_tf_inside_half", 8.8)
     monkeypatch.setattr(bv, "dr_cs", 0.65)
@@ -55,7 +55,7 @@ def test_pfcoil(monkeypatch, pfcoil):
     monkeypatch.setattr(bv, "r_tf_outboard_mid", 1.66e1)
     monkeypatch.setattr(bv, "dr_bore", 2.15)
     monkeypatch.setattr(fwbsv, "den_steel", 7.8e3)
-    monkeypatch.setattr(pfcoil_variables, "rpf1", 0.0)
+    monkeypatch.setattr(pfcoil_variables, "dr_pf_cs_middle_offset", 0.0)
     monkeypatch.setattr(pfcoil_variables, "m_pf_coil_structure_total", 0.0)
     monkeypatch.setattr(pfcoil_variables, "c_pf_cs_coil_flat_top_ma", np.full(22, 0.0))
     monkeypatch.setattr(pfcoil_variables, "n_cs_pf_coils", 0)
@@ -85,7 +85,8 @@ def test_pfcoil(monkeypatch, pfcoil):
     monkeypatch.setattr(
         pfcoil_variables, "c_pf_cs_coil_pulse_start_ma", np.full(22, 0.0)
     )
-    monkeypatch.setattr(pfcoil_variables, "routr", 1.5)
+    monkeypatch.setattr(pfcoil_variables, "dr_pf_tf_outboard_out_offset", 1.5)
+    monkeypatch.setattr(superconducting_tf_coil_variables, "r_tf_outboard_out", 10.0)
     monkeypatch.setattr(pfcoil_variables, "c_pf_cs_coils_peak_ma", np.full(22, 0.0))
     monkeypatch.setattr(pfcoil_variables, "f_j_cs_start_end_flat_top", 2.654e-1)
     monkeypatch.setattr(pfcoil_variables, "rpf2", -1.825)
@@ -122,7 +123,7 @@ def test_pfcoil(monkeypatch, pfcoil):
     monkeypatch.setattr(pfcoil_variables, "fcupfsu", 6.900e-1)
     monkeypatch.setattr(pfcoil_variables, "j_cs_pulse_start", 1.693e7)
     monkeypatch.setattr(pfcoil_variables, "j_pf_wp_critical", np.full(22, 0.0))
-    monkeypatch.setattr(pfcoil_variables, "i_sup_pf_shape", 0)
+    monkeypatch.setattr(pfcoil_variables, "i_r_pf_outside_tf_placement", 0)
     monkeypatch.setattr(pfcoil_variables, "rref", np.full(10, 7.0))
     monkeypatch.setattr(pfcoil_variables, "i_pf_current", 1)
     monkeypatch.setattr(pfcoil_variables, "ccl0_ma", np.full(10, 0.0))
@@ -378,7 +379,7 @@ def test_efc(pfcoil: PFCoil, monkeypatch: pytest.MonkeyPatch):
     # This 2D array argument discovered via gdb prints as a 1D array, therefore
     # needs to be reshaped into its original 2D. Fortran ordering is essential
     # when passing greater-than-1D arrays from Python to Fortran
-    rcls = np.reshape(
+    r_pf_coil_middle_group_array = np.reshape(
         [
             6.7651653417201345,
             6.7651653417201345,
@@ -404,7 +405,7 @@ def test_efc(pfcoil: PFCoil, monkeypatch: pytest.MonkeyPatch):
         (10, 2),
         order="F",
     )
-    zcls = np.reshape(
+    z_pf_coil_middle_group_array = np.reshape(
         [
             9.8904697261474404,
             -11.124884737289973,
@@ -447,8 +448,8 @@ def test_efc(pfcoil: PFCoil, monkeypatch: pytest.MonkeyPatch):
         cfix,
         n_pf_coil_groups,
         n_pf_coils_in_group,
-        rcls,
-        zcls,
+        r_pf_coil_middle_group_array,
+        z_pf_coil_middle_group_array,
         alfa,
         bfix,
         gmat,
@@ -519,7 +520,7 @@ def test_mtrx(pfcoil: PFCoil):
     bzin = np.zeros(nptsmx)
     n_pf_coil_groups = 4
     n_pf_coils_in_group = np.array([1, 1, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0])
-    rcls = np.reshape(
+    r_pf_coil_middle_group_array = np.reshape(
         [
             0,
             0,
@@ -545,7 +546,7 @@ def test_mtrx(pfcoil: PFCoil):
         (10, 2),
         order="F",
     )
-    zcls = np.reshape(
+    z_pf_coil_middle_group_array = np.reshape(
         [
             0,
             0,
@@ -659,8 +660,8 @@ def test_mtrx(pfcoil: PFCoil):
         bzin,
         n_pf_coil_groups,
         n_pf_coils_in_group,
-        rcls,
-        zcls,
+        r_pf_coil_middle_group_array,
+        z_pf_coil_middle_group_array,
         alfa,
         bfix,
         int(pfcoil_variables.N_PF_COILS_IN_GROUP_MAX),
diff --git a/tests/regression/input_files/spherical_tokamak_eval.IN.DAT b/tests/regression/input_files/spherical_tokamak_eval.IN.DAT
index 02178642ba..3d3087bc32 100644
--- a/tests/regression/input_files/spherical_tokamak_eval.IN.DAT
+++ b/tests/regression/input_files/spherical_tokamak_eval.IN.DAT
@@ -242,7 +242,7 @@ f_nd_impurity_electrons(14) = 5e-05
 i_pf_location   = 2,3,3,4 * Switch for location of PF coil group i;
 i_pf_conductor   = 0 * switch for PF &amp; CS coil conductor type;
 i_pf_superconductor = 9 * switch for superconductor material in PF coils;
-i_sup_pf_shape = 1 * Switch for the placement of Location 3 (outboard) PF coils
+i_r_pf_outside_tf_placement = 1 * Switch for the placement of Location 3 (outboard) PF coils
 n_pf_coils_in_group     = 2,2,2,2 * number of PF coils in group j
 n_pf_coil_groups     = 4 * number of groups of PF coils; Symmetric coil pairs should all be in the same group
 rref     = 7.0D0, 7.0D0, 7.0D0, 2.0, 7.0D0, 7.0D0, 7.0D0, 7.0D0, 7.0D0, 7.0D0 * PF coil radial positioning adjuster;
diff --git a/tests/regression/input_files/st_regression.IN.DAT b/tests/regression/input_files/st_regression.IN.DAT
index a3f7183e4c..82be122222 100644
--- a/tests/regression/input_files/st_regression.IN.DAT
+++ b/tests/regression/input_files/st_regression.IN.DAT
@@ -1865,7 +1865,7 @@ i_pf_location = 2,3,3,4
 * JUSTIFICATION: Design choice, one above x-point and two outboard
 
 
-i_sup_pf_shape = 1
+i_r_pf_outside_tf_placement = 1
 * DESCRIPTION:   Switch for the placement of i_pf_location = 3 (outboard) PF coils
 *  when the TF coils are superconducting (i_tf_sup = 1)
 *  =0 (Default) Outboard PF coils follow TF shape
@@ -1874,12 +1874,12 @@ i_sup_pf_shape = 1
 *  winding surface
 * JUSTIFICATION: Not used, no i_pf_location = 3 coils
 
-*routr =  
+*dr_pf_tf_outboard_out_offset =  
 * DESCRIPTION:  Radial distance (m) from outboard TF coil leg to centre of i_pf_location=3 PF coils
 * (default = 1.5) 
 * JUSTIFICATION: Not used, no i_pf_location = 3 coils
 
-*rpf1 =  
+*dr_pf_cs_middle_offset =  
 * DESCRIPTION:  offset (m) of radial position of i_pf_location=1 PF coils from being directly above
 * the central solenoid
 * JUSTIFICATION: Not used, no i_pf_location = 1 coils
diff --git a/tests/unit/test_pfcoil.py b/tests/unit/test_pfcoil.py
index 7bb602a56a..b68e43ffcb 100644
--- a/tests/unit/test_pfcoil.py
+++ b/tests/unit/test_pfcoil.py
@@ -2103,3 +2103,127 @@ def test_calculate_cs_geometry(
         dr_bore=dr_bore,
     )
     assert pytest.approx(result) == expected
+
+
+@pytest.mark.parametrize(
+    "n_pf_coils_in_group, n_pf_group, r_cs_middle, dr_pf_cs_middle_offset, z_tf_inside_half, dr_tf_inboard, z_cs_coil_upper, expected_r, expected_z",
+    [
+        # Single coil, above CS
+        (np.array([1, 0]), 0, 2.0, 0.1, 5.0, 0.2, 2.5, [2.1], [4.0]),
+        # Two coils, above CS
+        (np.array([2, 0]), 0, 3.0, 0.2, 6.0, 0.3, 3.0, [3.2, 3.2], [4.8, -4.8]),
+        # Single coil, different offset
+        (np.array([1, 0]), 0, 1.5, 0.05, 4.0, 0.1, 2.0, [1.55], [3.2]),
+        # Two coils, different offset and heights
+        (np.array([2, 0]), 0, 2.5, 0.15, 7.0, 0.25, 3.5, [2.65, 2.65], [5.525, -5.525]),
+    ],
+)
+def test_place_pf_above_cs(
+    pfcoil,
+    n_pf_coils_in_group,
+    n_pf_group,
+    r_cs_middle,
+    dr_pf_cs_middle_offset,
+    z_tf_inside_half,
+    dr_tf_inboard,
+    z_cs_coil_upper,
+    expected_r,
+    expected_z,
+):
+    r_array, z_array = pfcoil.place_pf_above_cs(
+        n_pf_coils_in_group=n_pf_coils_in_group,
+        n_pf_group=n_pf_group,
+        r_cs_middle=r_cs_middle,
+        dr_pf_cs_middle_offset=dr_pf_cs_middle_offset,
+        z_tf_inside_half=z_tf_inside_half,
+        dr_tf_inboard=dr_tf_inboard,
+        z_cs_coil_upper=z_cs_coil_upper,
+    )
+    # Only check the relevant group and number of coils
+    n_coils = n_pf_coils_in_group[n_pf_group]
+    assert_array_almost_equal(r_array[n_pf_group, :n_coils], expected_r)
+    assert_array_almost_equal(z_array[n_pf_group, :n_coils], expected_z)
+
+
+@pytest.mark.parametrize(
+    "n_pf_coils_in_group, n_pf_group, rmajor, triang, rminor, itart, itartpf, z_tf_inside_half, dz_tf_upper_lower_midplane, z_tf_top, top_bottom, rpf2, zref, expected_r, expected_z, expected_top_bottom",
+    [
+        # Test case 1: ST configuration, coil above midplane
+        (
+            np.array([2, 2]),
+            0,
+            3.0,
+            0.5,
+            1.0,
+            1,
+            0,
+            2.0,
+            0.5,
+            2.5,
+            1,
+            1.2,
+            np.array([0.3, 0.4]),
+            [3.6, 3.6],
+            [1.7, -1.7],
+            1,
+        ),
+        # Test case 2: Conventional configuration, coil above and below midplane
+        (
+            np.array([2, 2]),
+            1,
+            4.0,
+            0.3,
+            1.5,
+            0,
+            1,
+            2.5,
+            0.7,
+            3.0,
+            1,
+            1.0,
+            np.array([0.2, 0.5]),
+            [4.45, 4.45],
+            [3.86, -3.16],
+            1,
+        ),
+    ],
+)
+def test_place_pf_above_tf(
+    pfcoil,
+    n_pf_coils_in_group,
+    n_pf_group,
+    rmajor,
+    triang,
+    rminor,
+    itart,
+    itartpf,
+    z_tf_inside_half,
+    dz_tf_upper_lower_midplane,
+    z_tf_top,
+    top_bottom,
+    rpf2,
+    zref,
+    expected_r,
+    expected_z,
+    expected_top_bottom,
+):
+    r_arr, z_arr, top_bottom_out = pfcoil.place_pf_above_tf(
+        n_pf_coils_in_group=n_pf_coils_in_group,
+        n_pf_group=n_pf_group,
+        rmajor=rmajor,
+        triang=triang,
+        rminor=rminor,
+        itart=itart,
+        itartpf=itartpf,
+        z_tf_inside_half=z_tf_inside_half,
+        dz_tf_upper_lower_midplane=dz_tf_upper_lower_midplane,
+        z_tf_top=z_tf_top,
+        top_bottom=top_bottom,
+        rpf2=rpf2,
+        zref=zref,
+    )
+    # Only check the relevant group and number of coils
+    n_coils = n_pf_coils_in_group[n_pf_group]
+    assert np.allclose(r_arr[n_pf_group, :n_coils], expected_r)
+    assert np.allclose(z_arr[n_pf_group, :n_coils], expected_z)
+    assert top_bottom_out == expected_top_bottom