diff --git a/documentation/source/physics-models/profiles/plasma_density_profile.md b/documentation/source/physics-models/profiles/plasma_density_profile.md
index df5525e059..c801827681 100644
--- a/documentation/source/physics-models/profiles/plasma_density_profile.md
+++ b/documentation/source/physics-models/profiles/plasma_density_profile.md
@@ -148,4 +148,31 @@ $$\begin{aligned}
 
 5. Profile is then integrated with `integrate_profile_y()` using Simpsons integration from the profile abstract base class
 
+-----------------------
+
+### Setting pedestal and separatrix values | `set_pedestal_and_separatrix_values()`
+
+The switch `i_nd_plasma_pedestal_separatrix` controls how the values of the density pedestal and separatrix are set.
+
+#### User input
+
+If `i_nd_plasma_pedestal_separatrix == 0` then the values of `nd_plasma_pedestal_electron` and `nd_plasma_separatrix_electron` are taken directly from the input file
+
+#### Fraction of Greenwald Limit
+
+If `i_nd_plasma_pedestal_separatrix == 1`, the values of $n_{\text{ped}}$ and $n_{\text{sep}}$ are set as fractions of the [Greenwald](https://wiki.fusion.ciemat.es/wiki/Greenwald_limit) limit such as:
+
+$$
+n_{\text{ped}} = \overbrace{f_{\text{GW,ped}}}^{\texttt{f_nd_plasma_pedestal_greenwald}} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14}
+$$
+
+$$
+n_{\text{sep}} = \overbrace{f_{\text{GW,sep}}}^{\texttt{f_nd_plasma_separatrix_greenwald}} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14}
+$$
+
+
+`f_nd_plasma_pedestal_greenwald` and `f_nd_plasma_separatrix_greenwald` can be set as iteration variables respectively by using `ixc = 145`
+and `ixc = 152` respectively
+
+
 [^1]: Jean, J. (2011). *HELIOS: A Zero-Dimensional Tool for Next Step and Reactor Studies*. Fusion Science and Technology, 59(2), 308–349. <https://doi.org/10.13182/FST11-A11650>
diff --git a/documentation/source/physics-models/profiles/plasma_profiles.md b/documentation/source/physics-models/profiles/plasma_profiles.md
index c6514b5d51..1befc3407b 100644
--- a/documentation/source/physics-models/profiles/plasma_profiles.md
+++ b/documentation/source/physics-models/profiles/plasma_profiles.md
@@ -636,24 +636,6 @@ The same function is run from the `i_plasma_pedestal == 0 ` profile case, found
 
 --------
 
-### Setting pedestal values as fractions of the Greenwald limit
-
-By default, the values of $n_{\text{ped}}$ and $n_{\text{sep}}$ are set as fractions of the [Greenwald](https://wiki.fusion.ciemat.es/wiki/Greenwald_limit) limit such as:
-
-$$
-n_{\text{ped}} = \overbrace{f_{\text{GW,ped}}}^{\texttt{f_nd_plasma_pedestal_greenwald}} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14}
-$$
-
-$$
-n_{\text{sep}} = \overbrace{f_{\text{GW,sep}}}^{\texttt{f_nd_plasma_separatrix_greenwald}} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14}
-$$
-
-To set the values of $n_{\text{ped}}$ and $n_{\text{sep}}$ directly, the user can input the value of $\texttt{f_nd_plasma_pedestal_greenwald}$ or $\texttt{f_nd_plasma_separatrix_greenwald}$ to be less than 0.0 (i.e negative) to prevent the Greenwald fraction value being set.
-
-$\texttt{f_nd_plasma_pedestal_greenwald}$ and $\texttt{f_nd_plasma_separatrix_greenwald}$ can be set as iteration variables respectively by using `ixc = 45`
-and `ixc = 152` respectively
-
-------
 
 ### Pedestal Density Upper limit
 
diff --git a/process/core/init.py b/process/core/init.py
index 5d12d87f90..a7ed435b93 100644
--- a/process/core/init.py
+++ b/process/core/init.py
@@ -30,6 +30,7 @@
     init_superconducting_tf_coil_variables,
 )
 from process.data_structure.tfcoil_variables import init_tfcoil_variables
+from process.models.physics.profiles import DensityProfilePedestalType
 from process.models.stellarator.initialization import st_init
 from process.models.superconductors import (
     SuperconductorMaterial,
@@ -443,45 +444,53 @@ def check_process(inputs, data):  # noqa: ARG001
             )
 
         # Density checks
-        # Case where pedestal density is set manually
+        # Issue #589: Pedestal density is lower than separatrix density
+        pedestal_type = DensityProfilePedestalType(
+            data_structure.physics_variables.i_nd_plasma_pedestal_separatrix
+        )
         if (
-            data_structure.physics_variables.f_nd_plasma_pedestal_greenwald < 0
-            or not (
-                data_structure.numerics.ixc[: data_structure.numerics.nvar] == 145
-            ).any()
+            pedestal_type == DensityProfilePedestalType.USER_INPUT
+            and data_structure.physics_variables.nd_plasma_pedestal_electron
+            < data_structure.physics_variables.nd_plasma_separatrix_electron
+        ) or (
+            pedestal_type == DensityProfilePedestalType.GREENWALD_FRACTION
+            and data_structure.physics_variables.f_nd_plasma_pedestal_greenwald
+            < data_structure.physics_variables.f_nd_plasma_separatrix_greenwald
         ):
-            # Issue #589 Pedestal density is set manually using nd_plasma_pedestal_electron but it is less than nd_plasma_separatrix_electron.
-            if (
-                data_structure.physics_variables.nd_plasma_pedestal_electron
-                < data_structure.physics_variables.nd_plasma_separatrix_electron
-            ):
-                raise ProcessValidationError(
-                    "Density pedestal is lower than separatrix density",
-                    nd_plasma_pedestal_electron=data_structure.physics_variables.nd_plasma_pedestal_electron,
-                    nd_plasma_separatrix_electron=data_structure.physics_variables.nd_plasma_separatrix_electron,
-                )
+            raise ProcessValidationError(
+                "Density pedestal is lower than separatrix density",
+                **(
+                    {
+                        "nd_plasma_pedestal_electron": data_structure.physics_variables.nd_plasma_pedestal_electron,
+                        "nd_plasma_separatrix_electron": data_structure.physics_variables.nd_plasma_separatrix_electron,
+                    }
+                    if pedestal_type == DensityProfilePedestalType.USER_INPUT
+                    else {
+                        "f_nd_plasma_pedestal_greenwald": data_structure.physics_variables.f_nd_plasma_pedestal_greenwald,
+                        "f_nd_plasma_separatrix_greenwald": data_structure.physics_variables.f_nd_plasma_separatrix_greenwald,
+                    }
+                ),
+            )
 
-            # Issue #589 Pedestal density is set manually using nd_plasma_pedestal_electron,
-            # but pedestal width = 0.
-            if (
-                abs(
-                    data_structure.physics_variables.radius_plasma_pedestal_density_norm
-                    - 1.0
-                )
-                <= 1e-7
-                and (
-                    data_structure.physics_variables.nd_plasma_pedestal_electron
-                    - data_structure.physics_variables.nd_plasma_separatrix_electron
-                )
-                >= 1e-7
-            ):
-                warn(
-                    "Density pedestal is at plasma edge "
-                    f"({data_structure.physics_variables.radius_plasma_pedestal_density_norm = }), but nd_plasma_pedestal_electron "
-                    f"({data_structure.physics_variables.nd_plasma_pedestal_electron}) differs from "
-                    f"nd_plasma_separatrix_electron ({data_structure.physics_variables.nd_plasma_separatrix_electron})",
-                    stacklevel=2,
-                )
+        if (
+            abs(
+                data_structure.physics_variables.radius_plasma_pedestal_density_norm
+                - 1.0
+            )
+            <= 1e-7
+            and (
+                data_structure.physics_variables.nd_plasma_pedestal_electron
+                - data_structure.physics_variables.nd_plasma_separatrix_electron
+            )
+            >= 1e-7
+        ):
+            warn(
+                "Density pedestal is at plasma edge "
+                f"({data_structure.physics_variables.radius_plasma_pedestal_density_norm = }), but nd_plasma_pedestal_electron "
+                f"({data_structure.physics_variables.nd_plasma_pedestal_electron}) differs from "
+                f"nd_plasma_separatrix_electron ({data_structure.physics_variables.nd_plasma_separatrix_electron})",
+                stacklevel=2,
+            )
 
         # Issue #862 : Variable nd_plasma_electron_on_axis/nd_plasma_pedestal_electron ratio without constraint eq 81 (nd_plasma_electron_on_axis>nd_plasma_pedestal_electron)
         #  -> Potential hollowed density profile
diff --git a/process/core/input.py b/process/core/input.py
index 565a130876..2ffa9b73b7 100644
--- a/process/core/input.py
+++ b/process/core/input.py
@@ -176,10 +176,10 @@ def __post_init__(self):
     ),
     "ffwal": InputVariable(data_structure.physics_variables, float, range=(0.0, 10.0)),
     "f_nd_plasma_pedestal_greenwald": InputVariable(
-        data_structure.physics_variables, float, range=(-1.0, 5.0)
+        data_structure.physics_variables, float, range=(0.1, 1.5)
     ),
     "f_nd_plasma_separatrix_greenwald": InputVariable(
-        data_structure.physics_variables, float, range=(-1.0, 5.0)
+        data_structure.physics_variables, float, range=(0.001, 0.9)
     ),
     "f_plasma_fuel_helium3": InputVariable(
         data_structure.physics_variables, float, range=(-1.0, 5.0)
@@ -692,6 +692,9 @@ def __post_init__(self):
     "nd_plasma_separatrix_electron": InputVariable(
         data_structure.physics_variables, float, range=(0.0, 1e21)
     ),
+    "i_nd_plasma_pedestal_separatrix": InputVariable(
+        data_structure.physics_variables, int, choices=[0, 1]
+    ),
     "nflutfmax": InputVariable("constraints", float, range=(0.0, 1e24)),
     "oacdcp": InputVariable(
         data_structure.tfcoil_variables, float, range=(10000.0, 1000000000.0)
diff --git a/process/core/io/mfile/comparison.py b/process/core/io/mfile/comparison.py
index c55e95d80f..9ed37ecae4 100644
--- a/process/core/io/mfile/comparison.py
+++ b/process/core/io/mfile/comparison.py
@@ -117,7 +117,7 @@
     "temp_plasma_electron_on_axis_kev",
     "nd_plasma_electrons_vol_avg",
     "nd_plasma_electron_on_axis",
-    "dnla_gw",
+    "f_nd_plasma_greenwald",
     "temp_plasma_separatrix_kev",
     "nd_plasma_separatrix_electron",
     "temp_plasma_pedestal_kev",
diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index ac1e0951d5..ff26e4bd06 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -3163,7 +3163,7 @@ def plot_main_plasma_information(
         f"$\\mathbf{{Density \\ limit:}}$\n"
         f"({DensityLimitModel(int(mfile.get('i_density_limit', scan=scan))).full_name})\n"
         f"$n_{{\\text{{e,limit}}}}: {mfile.get('nd_plasma_electrons_max', scan=scan):.3e} \\ m^{{-3}}$\n"
-        f"$f_{{\\text{{GW}}}}$: {mfile.get('dnla_gw', scan=scan):.4f}"
+        f"$f_{{\\text{{GW}}}}$: {mfile.get('f_nd_plasma_greenwald', scan=scan):.4f}"
     )
 
     axis.text(
@@ -3857,8 +3857,12 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
     nd_ions_total = mfile.get("nd_plasma_ions_total_vol_avg", scan=scan)
     nd_fuel_ions = mfile.get("nd_plasma_fuel_ions_vol_avg", scan=scan)
     alphan = mfile.get("alphan", scan=scan)
-    fgwped_out = mfile.get("fgwped_out", scan=scan)
-    fgwsep_out = mfile.get("fgwsep_out", scan=scan)
+    f_nd_plasma_pedestal_greenwald = mfile.get(
+        "f_nd_plasma_pedestal_greenwald", scan=scan
+    )
+    f_nd_plasma_separatrix_greenwald = mfile.get(
+        "f_nd_plasma_separatrix_greenwald", scan=scan
+    )
     nd_plasma_electrons_vol_avg = mfile.get("nd_plasma_electrons_vol_avg", scan=scan)
     # find impurity densities
     imp_frac = np.array([
@@ -4048,7 +4052,10 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
 
     # Add text box with density profile parameters
     textstr_density = "\n".join((
-        rf"$\langle n_{{\text{{e}}}} \rangle$: {nd_plasma_electrons_vol_avg:.3e} m$^{{-3}}$",
+        (
+            rf"$\langle n_{{\text{{e}}}} \rangle$: {nd_plasma_electrons_vol_avg:.3e} m$^{{-3}}$"
+            rf"$\hspace{{4}} \overline{{n_{{e}}}}$: {mfile.get('nd_plasma_electron_line', scan=scan):.3e} m$^{{-3}}$"
+        ),
         (
             rf"$n_{{\text{{e,0}}}}$: {ne0:.3e} m$^{{-3}}$"
             rf"$\hspace{{4}} \alpha_{{\text{{n}}}}$: {alphan:.3f}"
@@ -4059,7 +4066,7 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
             f"{nd_fuel_ions / nd_plasma_electrons_vol_avg:.3f}"
         ),
         (
-            rf"$f_{{\text{{GW e,ped}}}}$: {fgwped_out:.3f}"
+            rf"$f_{{\text{{GW e,ped}}}}$: {f_nd_plasma_pedestal_greenwald:.3f}"
             r"$ \hspace{7} \frac{n_{e,0}}{\langle n_e \rangle}$: "
             f"{ne0 / nd_plasma_electrons_vol_avg:.3f}"
         ),
@@ -4069,12 +4076,12 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
             f"{mfile.get('nd_plasma_electron_line', scan=scan) / mfile.get('nd_plasma_electron_max_array(7)', scan=scan):.3f}"
         ),
         rf"$n_{{\text{{e,sep}}}}$: {nd_plasma_separatrix_electron:.3e} m$^{{-3}}$",
-        rf"$f_{{\text{{GW e,sep}}}}$: {fgwsep_out:.3f}",
+        rf"$f_{{\text{{GW e,sep}}}}$: {f_nd_plasma_separatrix_greenwald:.3f}",
     ))
 
     props_density = {"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5}
     ax_main.text(
-        0.0,
+        -0.05,
         -0.175,
         textstr_density,
         transform=ax_impurity.transAxes,
@@ -4293,33 +4300,34 @@ def plot_t_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
     # Add text box with temperature profile parameters
     textstr_temperature = "\n".join((
         (
-            rf"$\langle T_{{\text{{e}}}} \rangle_\text{{V}}$: {mfile.get('temp_plasma_electron_vol_avg_kev', scan=scan):.3f} keV"
-            rf"$\hspace{{3}} \langle T_{{\text{{e}}}} \rangle_\text{{n}}$: {mfile.get('temp_plasma_electron_density_weighted_kev', scan=scan):.3f} keV"
+            rf"$\langle T_{{\text{{e}}}} \rangle_\text{{V}}$:  {mfile.get('temp_plasma_electron_vol_avg_kev', scan=scan):.3f} keV"
+            rf"$\hspace{{2}} \langle T_{{\text{{e}}}} \rangle_\text{{n}}$: {mfile.get('temp_plasma_electron_density_weighted_kev', scan=scan):.3f} keV"
+            rf"$\hspace{{2}} \overline{{T_{{e}}}}$: {mfile.get('temp_plasma_electron_line_avg_kev', scan=scan):.3f} keV"
         ),
         (
-            rf"$T_{{\text{{e,0}}}}$: {te0:.3f} keV"
-            rf"$\hspace{{4}} \alpha_{{\text{{T}}}}$: {alphat:.3f}"
+            rf"$T_{{\text{{e,0}}}}$:    {te0:.3f} keV"
+            rf"$\hspace{{2}} \alpha_{{\text{{T}}}}$:   {alphat:.3f}"
         ),
         (
             rf"$T_{{\text{{e,ped}}}}$: {temp_plasma_pedestal_kev:.3f} keV"
-            r"$ \hspace{4} \frac{\langle T_i \rangle}{\langle T_e \rangle}$: "
+            r"$ \hspace{3} \frac{\langle T_i \rangle}{\langle T_e \rangle}$: "
             f"{f_temp_plasma_ion_electron:.3f}"
         ),
         (
             rf"$\rho_{{\text{{ped,T}}}}$: {radius_plasma_pedestal_temp_norm:.3f}"
-            r"$ \hspace{6} \frac{T_{e,0}}{\langle T_e \rangle}$: "
+            r"$ \hspace{5} \frac{T_{e,0}}{\langle T_e \rangle}$: "
             f"{te0 / te:.3f}"
         ),
         (
             rf"$T_{{\text{{e,sep}}}}$: {temp_plasma_separatrix_kev:.3f} keV"
-            r"$ \hspace{4} \frac{{{\langle T_e \rangle_n}}}{{{\langle T_e \rangle_V}}}$: "
+            r"$ \hspace{3} \frac{{{\langle T_e \rangle_n}}}{{{\langle T_e \rangle_V}}}$: "
             f"{mfile.get('f_temp_plasma_electron_density_vol_avg', scan=scan):.3f}"
         ),
     ))
 
     props_temperature = {"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5}
     prof.text(
-        0.0,
+        -0.1,
         -0.125,
         textstr_temperature,
         transform=prof.transAxes,
diff --git a/process/core/io/variable_metadata.py b/process/core/io/variable_metadata.py
index b0bf93117d..aed8636100 100644
--- a/process/core/io/variable_metadata.py
+++ b/process/core/io/variable_metadata.py
@@ -187,7 +187,7 @@ class VariableMetadata:
         description="Average electron density",
         units="m^{-3}",
     ),
-    "dnla_gw": VariableMetadata(
+    "f_nd_plasma_greenwald": VariableMetadata(
         latex=r"$f_{\mathrm{GW}}$", description="Greenwald fraction", units=""
     ),
     "normalised_toroidal_beta": VariableMetadata(
diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py
index 17392b4390..db9d70ad4b 100644
--- a/process/core/solver/iteration_variables.py
+++ b/process/core/solver/iteration_variables.py
@@ -235,10 +235,10 @@ class IterationVariable:
         1.0e20,
     ),
     145: IterationVariable(
-        "f_nd_plasma_pedestal_greenwald", data_structure.physics_variables, 0.1, 0.9
+        "f_nd_plasma_pedestal_greenwald", data_structure.physics_variables, 0.1, 1.5
     ),
     152: IterationVariable(
-        "f_nd_plasma_separatrix_greenwald", data_structure.physics_variables, 0.001, 0.5
+        "f_nd_plasma_separatrix_greenwald", data_structure.physics_variables, 0.001, 0.9
     ),
     155: IterationVariable("pfusife", "ife", 5.0e2, 3.0e3),
     156: IterationVariable("rrin", "ife", 1.0, 1.0e1),
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 6e12a4941e..4e3f146383 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -430,18 +430,17 @@ class DivertorNumberModels(IntEnum):
 load calculation (`i_pflux_fw_neutron=1`)
 """
 
+f_nd_plasma_greenwald: float = None
+"""Greenwald fraction of the line averaged electron density. The classic Greenwald
+limit value"""
 
 f_nd_plasma_pedestal_greenwald: float = None
-"""fraction of Greenwald density to set as pedestal-top density. If `<0`, pedestal-top
-density set manually using nd_plasma_pedestal_electron (`i_plasma_pedestal==1`).
-(`iteration variable 145`)
+"""Greenwald fraction of the pedestal density
 """
 
 
 f_nd_plasma_separatrix_greenwald: float = None
-"""fraction of Greenwald density to set as separatrix density. If `<0`, separatrix
-density set manually using nd_plasma_separatrix_electron (`i_plasma_pedestal==1`).
-(`iteration variable 152`)
+"""Greenwald fraction of the separatrix density
 """
 
 
@@ -611,6 +610,12 @@ class DivertorNumberModels(IntEnum):
 nd_plasma_separatrix_electron: float = None
 """electron density at separatrix [m-3] (`i_plasma_pedestal==1)"""
 
+i_nd_plasma_pedestal_separatrix: int = None
+"""switch for pedestal and separatrix density calculation:
+- =0 User input pedestal and separatrix density
+- =1 Calculate pedestal and separatrix density as fraction of Greenwald limit (see `f_nd_plasma_pedestal_greenwald` and `f_nd_plasma_separatrix_greenwald`)
+"""
+
 
 alpha_crit: float = None
 """critical ballooning parameter value"""
@@ -1217,6 +1222,8 @@ class DivertorNumberModels(IntEnum):
 temp_plasma_electron_density_weighted_kev: float = None
 """density weighted average electron temperature (keV)"""
 
+temp_plasma_electron_line_avg_kev: float = None
+"""Line averaged electron temperature (keV)"""
 
 temp_plasma_ion_vol_avg_kev: float = None
 """volume averaged ion temperature (keV). N.B. calculated from temp_plasma_electron_vol_avg_kev if `f_temp_plasma_ion_electron > 0.0`"""
@@ -1628,8 +1635,10 @@ def init_physics_variables():
         f_plasma_fuel_deuterium, \
         f_p_div_lower, \
         ffwal, \
+        f_nd_plasma_greenwald, \
         f_nd_plasma_pedestal_greenwald, \
         f_nd_plasma_separatrix_greenwald, \
+        i_nd_plasma_pedestal_separatrix, \
         f_plasma_fuel_helium3, \
         figmer, \
         fkzohm, \
@@ -1796,6 +1805,7 @@ def init_physics_variables():
         temp_plasma_electron_vol_avg_kev, \
         temp_plasma_electron_on_axis_kev, \
         temp_plasma_electron_density_weighted_kev, \
+        temp_plasma_electron_line_avg_kev, \
         temp_plasma_ion_vol_avg_kev, \
         temp_plasma_ion_on_axis_kev, \
         temp_plasma_ion_density_weighted_kev, \
@@ -1959,8 +1969,10 @@ def init_physics_variables():
     f_plasma_fuel_deuterium = 0.5
     f_p_div_lower = 1.0
     ffwal = 0.92
+    f_nd_plasma_greenwald = 1.0
     f_nd_plasma_pedestal_greenwald = 0.85
     f_nd_plasma_separatrix_greenwald = 0.50
+    i_nd_plasma_pedestal_separatrix = 1
     f_plasma_fuel_helium3 = 0.0
     figmer = 0.0
     fkzohm = 1.0
@@ -2129,6 +2141,7 @@ def init_physics_variables():
     temp_plasma_electron_vol_avg_kev = 12.9
     temp_plasma_electron_on_axis_kev = 0.0
     temp_plasma_electron_density_weighted_kev = 0.0
+    temp_plasma_electron_line_avg_kev = 0.0
     temp_plasma_ion_vol_avg_kev = 12.9
     temp_plasma_ion_on_axis_kev = 0.0
     temp_plasma_ion_density_weighted_kev = 0.0
diff --git a/process/models/physics/density_limit.py b/process/models/physics/density_limit.py
index 8ef3c46772..6b416e2918 100644
--- a/process/models/physics/density_limit.py
+++ b/process/models/physics/density_limit.py
@@ -90,7 +90,9 @@ def run(self):
             zeff=physics_variables.n_charge_plasma_effective_vol_avg,
         )
 
-        # Calculate beta_norm_max based on i_beta_norm_max
+        # Convert the chosen density limit model to the actual density limit value and
+        # store in physics variables. This is the value that will be used for all
+        # comparisons to the plasma density in the rest of the code.
         try:
             model = DensityLimitModel(int(physics_variables.i_density_limit))
             physics_variables.nd_plasma_electrons_max = self.get_density_limit_value(
@@ -102,6 +104,12 @@ def run(self):
                 i_density_limit=physics_variables.i_density_limit,
             ) from None
 
+        # Assign the Greenwald fraction for the rest of the code
+        physics_variables.f_nd_plasma_greenwald = (
+            physics_variables.nd_plasma_electron_line
+            / physics_variables.nd_plasma_electron_max_array[6]
+        )
+
     @staticmethod
     def get_density_limit_value(model: DensityLimitModel) -> float:
         """
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 34ce980ec6..0d6b24fce6 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -24,7 +24,10 @@
 )
 from process.models.physics import impurity_radiation
 from process.models.physics.plasma_geometry import PlasmaGeom
-from process.models.physics.profiles import PlasmaProfileShapeType
+from process.models.physics.profiles import (
+    DensityProfilePedestalType,
+    PlasmaProfileShapeType,
+)
 
 if TYPE_CHECKING:
     from process.models.physics.bootstrap_current import PlasmaBootstrapCurrent
@@ -358,24 +361,8 @@ def run(self):
         if (
             PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
             == PlasmaProfileShapeType.PEDESTAL_PROFILE
-        ) and (physics_variables.f_nd_plasma_pedestal_greenwald >= 0e0):
-            physics_variables.nd_plasma_pedestal_electron = (
-                physics_variables.f_nd_plasma_pedestal_greenwald
-                * 1.0e14
-                * physics_variables.plasma_current
-                / (np.pi * physics_variables.rminor * physics_variables.rminor)
-            )
-
-        if (
-            PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
-            == PlasmaProfileShapeType.PEDESTAL_PROFILE
-        ) and (physics_variables.f_nd_plasma_separatrix_greenwald >= 0e0):
-            physics_variables.nd_plasma_separatrix_electron = (
-                physics_variables.f_nd_plasma_separatrix_greenwald
-                * 1.0e14
-                * physics_variables.plasma_current
-                / (np.pi * physics_variables.rminor * physics_variables.rminor)
-            )
+        ):
+            self.plasma_profile.neprofile.set_pedestal_and_separatrix_values()
 
         self.plasma_profile.run()
 
@@ -2438,6 +2425,12 @@ def output_temperature_density_profile_info(self) -> None:
             physics_variables.temp_plasma_electron_on_axis_kev,
             "OP ",
         )
+        po.ovarrf(
+            self.outfile,
+            "Line averaged electron temperature (keV)",
+            "(temp_plasma_electron_line_avg_kev)",
+            physics_variables.temp_plasma_electron_line_avg_kev,
+        )
         po.ovarrf(
             self.outfile,
             "Volume averaged density weighted electron temperature (⟨Tₑ⟩ₙ) (keV)",
@@ -2486,68 +2479,90 @@ def output_temperature_density_profile_info(self) -> None:
         )
         po.oblnkl(self.outfile)
         if (
-            PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
+            physics_variables.i_plasma_pedestal
             == PlasmaProfileShapeType.PEDESTAL_PROFILE
         ):
+            po.ovarin(
+                self.outfile,
+                "Pedestal and separatrix density model selected",
+                "(i_nd_plasma_pedestal_separatrix)",
+                physics_variables.i_nd_plasma_pedestal_separatrix,
+            )
+            po.ocmmnt(
+                self.outfile,
+                f"Pedestal and separatrix values set by: "
+                f"{DensityProfilePedestalType(physics_variables.i_nd_plasma_pedestal_separatrix).description}",
+            )
+            po.oblnkl(self.outfile)
+
             po.ovarrf(
                 self.outfile,
                 "Density pedestal r/a location (ρₙ,pedestal)",
                 "(radius_plasma_pedestal_density_norm)",
                 physics_variables.radius_plasma_pedestal_density_norm,
             )
-            if physics_variables.f_nd_plasma_pedestal_greenwald >= 0e0:
-                po.ovarre(
+            if (
+                physics_variables.i_nd_plasma_pedestal_separatrix
+                == DensityProfilePedestalType.USER_INPUT
+            ):
+                po.ovarin(
                     self.outfile,
                     "Electron density pedestal height (nₑ_pedestal) (/m³)",
                     "(nd_plasma_pedestal_electron)",
                     physics_variables.nd_plasma_pedestal_electron,
+                )
+                po.ovarin(
+                    self.outfile,
+                    "Electron separatrix density (nₑ,ₛₑₚ) (/m³)",
+                    "(nd_plasma_separatrix_electron)",
+                    physics_variables.nd_plasma_separatrix_electron,
+                )
+                po.ovarre(
+                    self.outfile,
+                    "Pedestal Greenwald fraction",
+                    "(f_nd_plasma_pedestal_greenwald)",
+                    physics_variables.f_nd_plasma_pedestal_greenwald,
                     "OP ",
                 )
-            else:
+                po.ovarre(
+                    self.outfile,
+                    "Separatrix Greenwald fraction",
+                    "(f_nd_plasma_separatrix_greenwald)",
+                    physics_variables.f_nd_plasma_separatrix_greenwald,
+                    "OP ",
+                )
+
+            elif (
+                physics_variables.i_nd_plasma_pedestal_separatrix
+                == DensityProfilePedestalType.GREENWALD_FRACTION
+            ):
                 po.ovarre(
                     self.outfile,
                     "Electron density pedestal height (nₑ_pedestal) (/m³)",
                     "(nd_plasma_pedestal_electron)",
                     physics_variables.nd_plasma_pedestal_electron,
+                    "OP ",
                 )
-            # must be assigned to their exisiting values#
-            fgwped_out = (
-                physics_variables.nd_plasma_pedestal_electron
-                / physics_variables.nd_plasma_electron_max_array[6]
-            )
-            fgwsep_out = (
-                physics_variables.nd_plasma_separatrix_electron
-                / physics_variables.nd_plasma_electron_max_array[6]
-            )
-            if physics_variables.f_nd_plasma_pedestal_greenwald >= 0e0:
-                physics_variables.f_nd_plasma_pedestal_greenwald = (
-                    physics_variables.nd_plasma_pedestal_electron
-                    / physics_variables.nd_plasma_electron_max_array[6]
+                po.ovarre(
+                    self.outfile,
+                    "Electron separatrix density (nₑ,ₛₑₚ) (/m³)",
+                    "(nd_plasma_separatrix_electron)",
+                    physics_variables.nd_plasma_separatrix_electron,
+                    "OP ",
                 )
-            if physics_variables.f_nd_plasma_separatrix_greenwald >= 0e0:
-                physics_variables.f_nd_plasma_separatrix_greenwald = (
-                    physics_variables.nd_plasma_separatrix_electron
-                    / physics_variables.nd_plasma_electron_max_array[6]
+                po.ovarin(
+                    self.outfile,
+                    "Pedestal Greenwald fraction",
+                    "(f_nd_plasma_pedestal_greenwald)",
+                    physics_variables.f_nd_plasma_pedestal_greenwald,
+                )
+                po.ovarin(
+                    self.outfile,
+                    "Separatrix Greenwald fraction",
+                    "(f_nd_plasma_separatrix_greenwald)",
+                    physics_variables.f_nd_plasma_separatrix_greenwald,
                 )
 
-            po.ovarre(
-                self.outfile,
-                "Pedestal Greenwald fraction",
-                "(fgwped_out)",
-                fgwped_out,
-            )
-            po.ovarre(
-                self.outfile,
-                "Electron density at separatrix (nₑ,ₛₑₚ) (/m³)",
-                "(nd_plasma_separatrix_electron)",
-                physics_variables.nd_plasma_separatrix_electron,
-            )
-            po.ovarre(
-                self.outfile,
-                "Separatrix Greenwald fraction",
-                "(fgwsep_out)",
-                fgwsep_out,
-            )
             po.oblnkl(self.outfile)
 
         po.ovarre(
@@ -2574,9 +2589,8 @@ def output_temperature_density_profile_info(self) -> None:
             po.ovarre(
                 self.outfile,
                 "Greenwald fraction (f_GW)",
-                "(dnla_gw)",
-                physics_variables.nd_plasma_electron_line
-                / physics_variables.nd_plasma_electron_max_array[6],
+                "(f_nd_plasma_greenwald)",
+                physics_variables.f_nd_plasma_greenwald,
                 "OP ",
             )
         po.oblnkl(self.outfile)
diff --git a/process/models/physics/plasma_profiles.py b/process/models/physics/plasma_profiles.py
index 1956ed38af..f890fe1e86 100644
--- a/process/models/physics/plasma_profiles.py
+++ b/process/models/physics/plasma_profiles.py
@@ -143,6 +143,14 @@ def parabolic_paramterisation(self):
             / sp.special.gamma(physics_variables.alphan + 1.5)
         )
 
+        physics_variables.temp_plasma_electron_line_avg_kev = (
+            physics_variables.temp_plasma_electron_vol_avg_kev
+            * (1.0 + physics_variables.alphat)
+            * (sp.special.gamma(0.5) / 2.0)
+            * sp.special.gamma(physics_variables.alphat + 1.0)
+            / sp.special.gamma(physics_variables.alphat + 1.5)
+        )
+
         #  Density-weighted temperatures
 
         physics_variables.temp_plasma_electron_density_weighted_kev = (
@@ -224,11 +232,15 @@ def pedestal_parameterisation(self):
             / physics_variables.temp_plasma_electron_vol_avg_kev
         )
 
-        #  Line-averaged electron density
+        #  Line-averaged electron density and temperature
         #  = integral(n(rho).drho)
 
         physics_variables.nd_plasma_electron_line = self.neprofile.profile_integ
 
+        physics_variables.temp_plasma_electron_line_avg_kev = (
+            self.teprofile.profile_integ
+        )
+
         #  Scrape-off density / volume averaged density
         #  (Input value is used if i_plasma_pedestal = 0)
 
diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py
index 6344511098..1f93e300f0 100644
--- a/process/models/physics/profiles.py
+++ b/process/models/physics/profiles.py
@@ -13,6 +13,7 @@
 import scipy as sp
 
 from process.data_structure import physics_variables
+from process.models.physics.density_limit import PlasmaDensityLimit
 
 logger = logging.getLogger(__name__)
 
@@ -105,6 +106,31 @@ def integrate_profile_y(self):
         )
 
 
+class DensityProfilePedestalType(IntEnum):
+    """Enum for i_nd_plasma_pedestal_separatrix types"""
+
+    USER_INPUT = (0, "User input direct values")
+    GREENWALD_FRACTION = (1, "Fractions of the Greenwald limit")
+
+    def __new__(cls, value: int, description: str):
+        """Create a new DensityProfilePedestalType instance.
+
+        Parameters
+        ----------
+            value: Integer value for the enum member.
+            description: Human-readable description for the enum member.
+        """
+        obj = int.__new__(cls, value)
+        obj._value_ = value
+        obj._description_ = description
+        return obj
+
+    @DynamicClassAttribute
+    def description(self):
+        """Get the description of the plasma profile shape."""
+        return self._description_
+
+
 class NeProfile(Profile):
     """Electron density profile class. Contains a function to calculate the electron
     density profile and store the data.
@@ -250,6 +276,51 @@ def ncore(
             ncore = 1.0e-6
         return ncore
 
+    def set_pedestal_and_separatrix_values(self):
+        """Sets the pedestal and separatrix density values based on the user input or greenwald fraction method."""
+        i_nd_plasma_pedestal_separatrix = DensityProfilePedestalType(
+            physics_variables.i_nd_plasma_pedestal_separatrix
+        )
+
+        if i_nd_plasma_pedestal_separatrix == DensityProfilePedestalType.USER_INPUT:
+            physics_variables.f_nd_plasma_pedestal_greenwald = (
+                physics_variables.nd_plasma_pedestal_electron
+                / (
+                    PlasmaDensityLimit.calculate_greenwald_density_limit(
+                        c_plasma=physics_variables.plasma_current,
+                        rminor=physics_variables.rminor,
+                    )
+                )
+            )
+
+            physics_variables.f_nd_plasma_separatrix_greenwald = (
+                physics_variables.nd_plasma_separatrix_electron
+                / (
+                    PlasmaDensityLimit.calculate_greenwald_density_limit(
+                        c_plasma=physics_variables.plasma_current,
+                        rminor=physics_variables.rminor,
+                    )
+                )
+            )
+        elif (
+            i_nd_plasma_pedestal_separatrix
+            == DensityProfilePedestalType.GREENWALD_FRACTION
+        ):
+            physics_variables.nd_plasma_pedestal_electron = (
+                physics_variables.f_nd_plasma_pedestal_greenwald
+                * PlasmaDensityLimit.calculate_greenwald_density_limit(
+                    c_plasma=physics_variables.plasma_current,
+                    rminor=physics_variables.rminor,
+                )
+            )
+            physics_variables.nd_plasma_separatrix_electron = (
+                physics_variables.f_nd_plasma_separatrix_greenwald
+                * PlasmaDensityLimit.calculate_greenwald_density_limit(
+                    c_plasma=physics_variables.plasma_current,
+                    rminor=physics_variables.rminor,
+                )
+            )
+
     def set_physics_variables(self):
         """Calculates and sets physics variables required for the profile."""
         if (
diff --git a/tests/unit/models/physics/test_plasma_profiles.py b/tests/unit/models/physics/test_plasma_profiles.py
index 562792b275..31dcad54f1 100644
--- a/tests/unit/models/physics/test_plasma_profiles.py
+++ b/tests/unit/models/physics/test_plasma_profiles.py
@@ -216,6 +216,8 @@ class PlasmaProfilesParam(NamedTuple):
 
     expected_nd_electron_line: float = 0.0
 
+    expected_temp_plasma_electron_line_avg_kev: float = 0.0
+
     expected_ti: float = 0.0
 
 
@@ -265,6 +267,7 @@ class PlasmaProfilesParam(NamedTuple):
             expected_ne0=1.0585658890823703e20,
             expected_ti0=27.370104119511087,
             expected_nd_electron_line=8.8687354645836431e19,
+            expected_temp_plasma_electron_line_avg_kev=17.582796426026842,
             expected_ti=13.07,
         ),
         PlasmaProfilesParam(
@@ -310,6 +313,7 @@ class PlasmaProfilesParam(NamedTuple):
             expected_ne0=1.0585658890823703e20,
             expected_ti0=27.370104119511087,
             expected_nd_electron_line=8.8687354645836431e19,
+            expected_temp_plasma_electron_line_avg_kev=17.582796426026842,
             expected_ti=13.07,
         ),
     ],
@@ -520,6 +524,10 @@ def test_plasma_profiles(plasmaprofilesparam, monkeypatch, plasmaprofile):
         plasmaprofilesparam.expected_nd_electron_line
     )
 
+    assert physics_variables.temp_plasma_electron_line_avg_kev == pytest.approx(
+        plasmaprofilesparam.expected_temp_plasma_electron_line_avg_kev
+    )
+
     assert physics_variables.temp_plasma_ion_vol_avg_kev == pytest.approx(
         plasmaprofilesparam.expected_ti
     )