diff --git a/species/data/companion_data.json b/species/data/companion_data.json
index f79ad530..0b3eb665 100644
--- a/species/data/companion_data.json
+++ b/species/data/companion_data.json
@@ -63,7 +63,7 @@
             "Males et al. 2014",
             "Mirek Brandt et al. 2021",
             "Stolker et al. 2019",
-            "Stolker et al. 2020"
+            "Stolker et al. 2020a"
         ]
     },
     "beta Pic c": {
@@ -151,7 +151,7 @@
             "Cheetham et al. 2019",
             "Gaia Early Data Release 3",
             "Marleau et al. 2019",
-            "Stolker et al. 2020"
+            "Stolker et al. 2020a"
         ]
     },
     "51 Eri b": {
@@ -613,8 +613,7 @@
             "Gaia Early Data Release 3",
             "Haffert et al. 2019",
             "Hashimoto et al. 2020",
-            "Stolker et al. 2020",
-            "Stolker et al. 2020.",
+            "Stolker et al. 2020b",
             "Wang et al. 2020",
             "Wang et al. 2021"
         ]
@@ -651,7 +650,7 @@
         "references": [
             "Gaia Early Data Release 3",
             "Haffert et al. 2019",
-            "Stolker et al. 2020",
+            "Stolker et al. 2020b",
             "Wang et al. 2020",
             "Wang et al. 2021"
         ]
@@ -807,7 +806,7 @@
             "Gaia Early Data Release 3",
             "Grandjean et al. 2019",
             "Milli et al. 2017",
-            "Stolker et al. 2020",
+            "Stolker et al. 2020a",
             "Ward-Duong et al. 2020"
         ]
     },
@@ -1065,7 +1064,7 @@
             "Maire et al. 2015",
             "Maire et al. 2016",
             "Musso Barcucci et al. 2019",
-            "Stolker et al. 2020"
+            "Stolker et al. 2020a"
         ]
     },
     "kappa And b": {
@@ -2008,7 +2007,7 @@
             0.0
         ],
         "mass_companion": [
-            19.0,
+            17.0,
             5.0
         ],
         "accretion": true,
@@ -2016,7 +2015,8 @@
             "Gaia Early Data Release 3",
             "Hartmann et al. 1998",
             "Schmidt et al. 2008",
-            "Wu et al. 2015"
+            "Wu et al. 2015",
+            "Wu et al. 2020"
         ]
     },
     "SR 12 C": {
diff --git a/species/data/database.py b/species/data/database.py
index 878cd5e6..50448c95 100644
--- a/species/data/database.py
+++ b/species/data/database.py
@@ -6,6 +6,7 @@
 import json
 import os
 import pathlib
+import sys
 # import urllib.error
 import warnings
 
@@ -2632,10 +2633,10 @@ def get_pt_profiles(
         self, tag: str, random: Optional[int] = None, out_file: Optional[str] = None
     ) -> Tuple[np.ndarray, np.ndarray]:
         """
-        Function for returning the pressure-temperature profiles from
-        the posterior of the atmospheric retrieval with
-        ``petitRADTRANS``. The data can also optionally be written to
-        an output file.
+        Function for returning the pressure-temperature profiles
+        from the posterior of the atmospheric retrieval with
+        ``petitRADTRANS``. The data can also optionally be
+        written to an output file.
 
         Parameters
         ----------
@@ -2682,6 +2683,11 @@ def get_pt_profiles(
         elif "nparam" in dset.attrs:
             n_param = dset.attrs["nparam"]
 
+        if "temp_nodes" in dset.attrs:
+            temp_nodes = dset.attrs["temp_nodes"]
+        else:
+            temp_nodes = 15
+
         samples = np.asarray(dset)
 
         if random is None:
@@ -2752,15 +2758,15 @@ def get_pt_profiles(
                 else:
                     pt_smooth = 0.0
 
-                knot_press = np.logspace(np.log10(press[0]), np.log10(press[-1]), 15)
+                knot_press = np.logspace(np.log10(press[0]), np.log10(press[-1]), temp_nodes)
 
                 knot_temp = []
-                for j in range(15):
-                    knot_temp.append(item[param_index[f"t{i}"]])
+                for k in range(temp_nodes):
+                    knot_temp.append(item[param_index[f"t{k}"]])
 
                 knot_temp = np.asarray(knot_temp)
 
-                temp[:, j] = retrieval_util.pt_spline_interp(
+                temp[:, i] = retrieval_util.pt_spline_interp(
                     knot_press, knot_temp, press, pt_smooth
                 )
 
@@ -3132,8 +3138,10 @@ def add_retrieval(
 
             if "temp_nodes" not in radtrans or radtrans["temp_nodes"] is None:
                 dset.attrs["temp_nodes"] = "None"
+                temp_nodes = 15
             else:
                 dset.attrs["temp_nodes"] = radtrans["temp_nodes"]
+                temp_nodes = radtrans["temp_nodes"]
 
             if "pt_smooth" in radtrans:
                 dset.attrs["pt_smooth"] = radtrans["pt_smooth"]
@@ -3165,10 +3173,10 @@ def add_retrieval(
                     print(" [DONE]")
 
                     print("Importing chemistry module...", end="", flush=True)
-                    from poor_mans_nonequ_chem_FeH.poor_mans_nonequ_chem.poor_mans_nonequ_chem import (
-                        interpol_abundances,
-                    )
-
+                    if "poor_mans_nonequ_chem" in sys.modules:
+                        from poor_mans_nonequ_chem.poor_mans_nonequ_chem import interpol_abundances
+                    else:
+                        from petitRADTRANS.poor_mans_nonequ_chem.poor_mans_nonequ_chem import interpol_abundances
                     print(" [DONE]")
 
                     rt_object = Radtrans(
@@ -3221,11 +3229,11 @@ def add_retrieval(
                         or radtrans["pt_profile"] == "monotonic"
                     ):
                         knot_press = np.logspace(
-                            np.log10(pressure[0]), np.log10(pressure[-1]), 15
+                            np.log10(pressure[0]), np.log10(pressure[-1]), temp_nodes
                         )
 
                         knot_temp = []
-                        for k in range(15):
+                        for k in range(temp_nodes):
                             knot_temp.append(sample_dict[f"t{k}"])
 
                         knot_temp = np.asarray(knot_temp)
@@ -3325,11 +3333,11 @@ def add_retrieval(
                     or radtrans["pt_profile"] == "monotonic"
                 ):
                     knot_press = np.logspace(
-                        np.log10(pressure[0]), np.log10(pressure[-1]), 15
+                        np.log10(pressure[0]), np.log10(pressure[-1]), temp_nodes
                     )
 
                     knot_temp = []
-                    for k in range(15):
+                    for k in range(temp_nodes):
                         knot_temp.append(sample_dict[f"t{k}"])
 
                     knot_temp = np.asarray(knot_temp)
@@ -3862,6 +3870,16 @@ def petitcode_param(
         else:
             p_quench = None
 
+        if "temp_nodes" in sample_box.attributes:
+            temp_nodes = sample_box.attributes["temp_nodes"]
+        else:
+            temp_nodes = 15
+
+        if "pressure_grid" in sample_box.attributes:
+            pressure_grid = sample_box.attributes["pressure_grid"]
+        else:
+            pressure_grid = "smaller"
+
         pressure = np.logspace(-6.0, 3.0, 180)
 
         if sample_box.attributes["pt_profile"] == "molliere":
@@ -3876,10 +3894,10 @@ def petitcode_param(
             )
 
         else:
-            knot_press = np.logspace(np.log10(pressure[0]), np.log10(pressure[-1]), 15)
+            knot_press = np.logspace(np.log10(pressure[0]), np.log10(pressure[-1]), temp_nodes)
 
             knot_temp = []
-            for i in range(15):
+            for i in range(temp_nodes):
                 knot_temp.append(model_param[f"t{i}"])
 
             knot_temp = np.asarray(knot_temp)
@@ -3893,7 +3911,10 @@ def petitcode_param(
                 knot_press, knot_temp, pressure, pt_smooth=pt_smooth
             )
 
-        from poor_mans_nonequ_chem.poor_mans_nonequ_chem import interpol_abundances
+        if "poor_mans_nonequ_chem" in sys.modules:
+            from poor_mans_nonequ_chem.poor_mans_nonequ_chem import interpol_abundances
+        else:
+            from petitRADTRANS.poor_mans_nonequ_chem.poor_mans_nonequ_chem import interpol_abundances
 
         # Interpolate the abundances, following chemical equilibrium
         abund_in = interpol_abundances(
@@ -3976,7 +3997,7 @@ def petitcode_param(
             cloud_species=cloud_species,
             scattering=True,
             wavel_range=(0.5, 50.0),
-            pressure_grid=sample_box.attributes["pressure_grid"],
+            pressure_grid=pressure_grid,
             res_mode="c-k",
             cloud_wavel=cloud_wavel,
         )
diff --git a/species/plot/plot_retrieval.py b/species/plot/plot_retrieval.py
index acf13857..c8a96c8c 100644
--- a/species/plot/plot_retrieval.py
+++ b/species/plot/plot_retrieval.py
@@ -496,9 +496,7 @@ def plot_pt_profile(
         if "poor_mans_nonequ_chem" in sys.modules:
             from poor_mans_nonequ_chem.poor_mans_nonequ_chem import interpol_abundances
         else:
-            from petitRADTRANS.poor_mans_nonequ_chem.poor_mans_nonequ_chem import (
-                interpol_abundances,
-            )
+            from petitRADTRANS.poor_mans_nonequ_chem.poor_mans_nonequ_chem import interpol_abundances
 
         abund_in = interpol_abundances(
             np.full(pressure.shape[0], median["c_o_ratio"]),
@@ -526,17 +524,19 @@ def plot_pt_profile(
                     pressure_grid=radtrans.pressure_grid,
                 )
 
-        for item in cloud_species:
-            if item in radtrans.cloud_species:
-                sat_press, sat_temp = retrieval_util.return_T_cond_Fe_comb(
-                    median["metallicity"],
-                    median["c_o_ratio"],
-                    MMW=np.mean(abund_in["MMW"]),
-                )
+        for cloud_item in cloud_species:
+
+            if cloud_item in radtrans.cloud_species:
+                cond_temp = retrieval_util.get_condensation_curve(
+                    composition=cloud_item[:-3],
+                    press=pressure,
+                    metallicity=median["metallicity"],
+                    c_o_ratio=median["c_o_ratio"],
+                    mmw=np.mean(abund_in["MMW"]))
 
                 ax.plot(
-                    sat_temp,
-                    sat_press,
+                    cond_temp,
+                    pressure,
                     "--",
                     lw=0.8,
                     color=next(color_iter, "black"),
@@ -721,12 +721,12 @@ def plot_pt_profile(
 
             color_iter = iter(cloud_colors)
 
-            for item in cloud_species:
-                if item in radtrans.cloud_species:
-                    cloud_index = radtrans.rt_object.cloud_species.index(item)
+            for cloud_item in cloud_species:
+                if cloud_item in radtrans.cloud_species:
+                    cloud_index = radtrans.rt_object.cloud_species.index(cloud_item)
 
                     label = ""
-                    for char in item[:-3]:
+                    for char in cloud_item[:-3]:
                         if char.isnumeric():
                             label += f"$_{char}$"
                         else:
diff --git a/species/plot/plot_spectrum.py b/species/plot/plot_spectrum.py
index 3c05ad0b..61a276c3 100644
--- a/species/plot/plot_spectrum.py
+++ b/species/plot/plot_spectrum.py
@@ -48,6 +48,7 @@ def plot_spectrum(
     output: Optional[str] = "spectrum.pdf",
     leg_param: Optional[List[str]] = None,
     grid_hspace: float = 0.1,
+    inc_model_name: bool = False,
 ):
     """
     Function for plotting a spectral energy distribution and combining
@@ -147,6 +148,9 @@ def plot_spectrum(
         The relative height spacing between subplots, expressed
         as a fraction of the average axis height. The default
         value is set to 0.1.
+    inc_model_name : bool
+        Include the model name in the legend of any
+        :class:`~species.core.box.ModelBox`.
 
     Returns
     -------
@@ -467,6 +471,10 @@ def plot_spectrum(
                             if item not in leg_param:
                                 del param[item]
 
+                    if leg_param is not None:
+                        param_new = {k: param[k] for k in leg_param}
+                        param = param_new.copy()
+
                     par_key, par_unit, par_label = plot_util.quantity_unit(
                         param=list(param.keys()), object_type=object_type
                     )
@@ -475,6 +483,7 @@ def plot_spectrum(
                     # newline = False
 
                     for i, item in enumerate(par_key):
+
                         if item[:4] == "teff":
                             value = f"{param[item]:.0f}"
 
@@ -534,15 +543,21 @@ def plot_spectrum(
                         #     label += '\n'
                         #     newline = True
 
+                        model_name = plot_util.model_name(box_item.model)
+
                         if par_unit[i] is None:
                             if len(label) > 0:
                                 label += ", "
+                            elif inc_model_name:
+                                label += f"{model_name}: "
 
                             label += f"{par_label[i]} = {value}"
 
                         else:
                             if len(label) > 0:
                                 label += ", "
+                            elif inc_model_name:
+                                label += f"{model_name}: "
 
                             label += f"{par_label[i]} = {value} {par_unit[i]}"
 
diff --git a/species/util/plot_util.py b/species/util/plot_util.py
index c5244838..3f8b8378 100644
--- a/species/util/plot_util.py
+++ b/species/util/plot_util.py
@@ -693,13 +693,13 @@ def update_labels(param: List[str], object_type: str = "planet") -> List[str]:
 
 
 @typechecked
-def model_name(key: str) -> str:
+def model_name(in_name: str) -> str:
     """
     Function for updating a model name for use in plots.
 
     Parameters
     ----------
-    key : str
+    in_name : str
         Model name as used by species.
 
     Returns
@@ -708,64 +708,70 @@ def model_name(key: str) -> str:
         Updated model name for plots.
     """
 
-    if key == "drift-phoenix":
-        name = "DRIFT-PHOENIX"
+    if in_name == "drift-phoenix":
+        out_name = "DRIFT-PHOENIX"
 
-    elif key == "ames-cond":
-        name = "AMES-Cond"
+    elif in_name == "ames-cond":
+        out_name = "AMES-Cond"
 
-    elif key == "ames-dusty":
-        name = "AMES-Dusty"
+    elif in_name == "ames-dusty":
+        out_name = "AMES-Dusty"
 
-    elif key == "atmo":
-        name = "ATMO"
+    elif in_name == "atmo":
+        out_name = "ATMO"
 
-    elif key == "bt-cond":
-        name = "BT-Cond"
+    elif in_name == "bt-cond":
+        out_name = "BT-Cond"
 
-    elif key == "bt-cond-feh":
-        name = "BT-Cond"
+    elif in_name == "bt-cond-feh":
+        out_name = "BT-Cond"
 
-    elif key == "bt-settl":
-        name = "BT-Settl"
+    elif in_name == "bt-settl":
+        out_name = "BT-Settl"
 
-    elif key == "bt-settl-cifist":
-        name = "BT-Settl"
+    elif in_name == "bt-settl-cifist":
+        out_name = "BT-Settl"
 
-    elif key == "bt-nextgen":
-        name = "BT-NextGen"
+    elif in_name == "bt-nextgen":
+        out_name = "BT-NextGen"
 
-    elif key == "petitcode-cool-clear":
-        name = "petitCODE"
+    elif in_name == "petitcode-cool-clear":
+        out_name = "petitCODE"
 
-    elif key == "petitcode-cool-cloudy":
-        name = "petitCODE"
+    elif in_name == "petitcode-cool-cloudy":
+        out_name = "petitCODE"
 
-    elif key == "petitcode-hot-clear":
-        name = "petitCODE"
+    elif in_name == "petitcode-hot-clear":
+        out_name = "petitCODE"
 
-    elif key == "petitcode-hot-cloudy":
-        name = "petitCODE"
+    elif in_name == "petitcode-hot-cloudy":
+        out_name = "petitCODE"
 
-    elif key == "exo-rem":
-        name = "Exo-REM"
+    elif in_name == "exo-rem":
+        out_name = "Exo-REM"
 
-    elif key == "planck":
-        name = "Blackbody radiation"
+    elif in_name == "planck":
+        out_name = "Blackbody radiation"
 
-    elif key == "zhu2015":
-        name = "Zhu (2015)"
+    elif in_name == "zhu2015":
+        out_name = "Zhu (2015)"
 
-    elif key == "sonora-cholla":
-        name = "Sonora Cholla"
+    elif in_name == "sonora-cholla":
+        out_name = "Sonora Cholla"
 
-    elif key == "sonora-bobcat":
-        name = "Sonora Bobcat"
+    elif in_name == "sonora-bobcat":
+        out_name = "Sonora Bobcat"
 
-    elif key == "sonora-bobcat-co":
-        name = "Sonora Bobcat C/O"
+    elif in_name == "sonora-bobcat-co":
+        out_name = "Sonora Bobcat C/O"
 
-    return name
+    elif in_name == "petitradtrans":
+        out_name = "petitRADTRANS"
+
+    else:
+        raise ValueError(f"The model name '{in_name}' is not known.")
+
+    return out_name
 
 
 @typechecked
@@ -797,59 +803,39 @@ def quantity_unit(
     unit = []
     label = []
 
-    if "teff" in param:
-        quantity.append("teff")
-        unit.append("K")
-        label.append(r"$T_\mathrm{eff}$")
-
-    if "logg" in param:
-        quantity.append("logg")
-        unit.append(None)
-        label.append(r"$\log g$")
-
-    if "metallicity" in param:
-        quantity.append("metallicity")
-        unit.append(None)
-        label.append("[Fe/H]")
-
-    if "feh" in param:
-        quantity.append("feh")
-        unit.append(None)
-        label.append("[Fe/H]")
-
-    if "fsed" in param:
-        quantity.append("fsed")
-        unit.append(None)
-        label.append(r"$f_\mathrm{sed}$")
+    for item in param:
+        if item == "teff":
+            quantity.append("teff")
+            unit.append("K")
+            label.append(r"$T_\mathrm{eff}$")
 
-    if "c_o_ratio" in param:
-        quantity.append("c_o_ratio")
-        unit.append(None)
-        label.append("C/O")
+        if item == "logg":
+            quantity.append("logg")
+            unit.append(None)
+            label.append(r"$\log g$")
 
-    if "radius" in param:
-        quantity.append("radius")
+        if item == "metallicity":
+            quantity.append("metallicity")
+            unit.append(None)
+            label.append("[Fe/H]")
 
-        if object_type == "planet":
-            unit.append(r"$R_\mathrm{J}$")
+        if item == "feh":
+            quantity.append("feh")
+            unit.append(None)
+            label.append("[Fe/H]")
 
-        elif object_type == "star":
-            unit.append(r"$R_\mathrm{\odot}$")
+        if item == "fsed":
+            quantity.append("fsed")
+            unit.append(None)
+            label.append(r"$f_\mathrm{sed}$")
 
-        label.append(r"$R$")
+        if item == "c_o_ratio":
+            quantity.append("c_o_ratio")
+            unit.append(None)
+            label.append("C/O")
 
-    for i in range(100):
-        if f"teff_{i}" in param:
-            quantity.append(f"teff_{i}")
-            unit.append("K")
-            label.append(rf"$T_\mathrm{{{i+1}}}$")
-
-        else:
-            break
-
-    for i in range(100):
-        if f"radius_{i}" in param:
-            quantity.append(f"radius_{i}")
+        if item == "radius":
+            quantity.append("radius")
 
             if object_type == "planet":
                 unit.append(r"$R_\mathrm{J}$")
@@ -857,51 +843,72 @@ def quantity_unit(
             elif object_type == "star":
                 unit.append(r"$R_\mathrm{\odot}$")
 
-            label.append(rf"$R_\mathrm{{{i+1}}}$")
+            label.append(r"$R$")
 
-        else:
-            break
+        for i in range(100):
+            if item == f"teff_{i}":
+                quantity.append(f"teff_{i}")
+                unit.append("K")
+                label.append(rf"$T_\mathrm{{{i+1}}}$")
 
-    if "distance" in param:
-        quantity.append("distance")
-        unit.append("pc")
-        label.append(r"$d$")
+            else:
+                break
 
-    if "mass" in param:
-        quantity.append("mass")
+        for i in range(100):
+            if item == f"radius_{i}":
+                quantity.append(f"radius_{i}")
 
-        if object_type == "planet":
-            unit.append(r"$M_\mathrm{J}$")
+                if object_type == "planet":
+                    unit.append(r"$R_\mathrm{J}$")
 
-        elif object_type == "star":
-            unit.append(r"$M_\mathrm{\odot}$")
+                elif object_type == "star":
+                    unit.append(r"$R_\mathrm{\odot}$")
 
-        label.append("M")
+                label.append(rf"$R_\mathrm{{{i+1}}}$")
 
-    if "luminosity" in param:
-        quantity.append("luminosity")
-        unit.append(None)
-        label.append(r"$\log\,L/L_\mathrm{\odot}$")
+            else:
+                break
 
-    if "ism_ext" in param:
-        quantity.append("ism_ext")
-        unit.append(None)
-        label.append(r"$A_V$")
+        if item == "distance":
+            quantity.append("distance")
+            unit.append("pc")
+            label.append(r"$d$")
 
-    if "lognorm_ext" in param:
-        quantity.append("lognorm_ext")
-        unit.append(None)
-        label.append(r"$A_V$")
+        if item == "mass":
+            quantity.append("mass")
 
-    if "powerlaw_ext" in param:
-        quantity.append("powerlaw_ext")
-        unit.append(None)
-        label.append(r"$A_V$")
+            if object_type == "planet":
+                unit.append(r"$M_\mathrm{J}$")
 
-    if "pt_smooth" in param:
-        quantity.append("pt_smooth")
-        unit.append(None)
-        label.append(r"$\sigma_\mathrm{P-T}$")
+            elif object_type == "star":
+                unit.append(r"$M_\mathrm{\odot}$")
+
+            label.append("M")
+
+        if item == "luminosity":
+            quantity.append("luminosity")
+            unit.append(None)
+            label.append(r"$\log\,L/L_\mathrm{\odot}$")
+
+        if item == "ism_ext":
+            quantity.append("ism_ext")
+            unit.append(None)
+            label.append(r"$A_V$")
+
+        if item == "lognorm_ext":
+            quantity.append("lognorm_ext")
+            unit.append(None)
+            label.append(r"$A_V$")
+
+        if item == "powerlaw_ext":
+            quantity.append("powerlaw_ext")
+            unit.append(None)
+            label.append(r"$A_V$")
+
+        if item == "pt_smooth":
+            quantity.append("pt_smooth")
+            unit.append(None)
+            label.append(r"$\sigma_\mathrm{P-T}$")
 
     return quantity, unit, label
 
diff --git a/species/util/retrieval_util.py b/species/util/retrieval_util.py
index 7f50e3ae..91a5b062 100644
--- a/species/util/retrieval_util.py
+++ b/species/util/retrieval_util.py
@@ -1933,15 +1933,13 @@ def cloud_mass_fraction(
 
 
 @typechecked
-def find_cloud_deck(
+def get_condensation_curve(
     composition: str,
     press: np.ndarray,
-    temp: np.ndarray,
     metallicity: float,
     c_o_ratio: float,
     mmw: float = 2.33,
-    plotting: bool = False,
-) -> float:
+) -> np.ndarray:
     """
     Function to find the base of the cloud deck by intersecting the
     P-T profile with the saturation vapor pressure.
@@ -1953,21 +1951,17 @@ def find_cloud_deck(
         'Na2S', or 'KCL').
     press : np.ndarray
         Pressures (bar).
-    temp : np.ndarray
-        Temperatures (K).
     metallicity : float
         Metallicity [Fe/H].
     c_o_ratio : float
         Carbon-to-oxygen ratio.
     mmw : float
         Mean molecular weight.
-    plotting : bool
-        Create a plot.
 
     Returns
     -------
-    float
-        Pressure (bar) at the base of the cloud deck.
+    np.array
+        Condensation temperatures (K) for the provided input pressures.
     """
 
     if composition == "Fe":
@@ -1990,14 +1984,62 @@ def find_cloud_deck(
 
     else:
         raise ValueError(
-            f"The '{composition}' composition is not supported by find_cloud_deck."
+            f"The '{composition}' composition is not "
+            "supported by get_condensation_curve."
         )
 
     index = (Pc > 1e-8) & (Pc < 1e5)
     Pc, Tc = Pc[index], Tc[index]
 
     tcond_p = interp1d(Pc, Tc)
-    Tcond_on_input_grid = tcond_p(press)
+
+    return tcond_p(press)
+
+
+@typechecked
+def find_cloud_deck(
+    composition: str,
+    press: np.ndarray,
+    temp: np.ndarray,
+    metallicity: float,
+    c_o_ratio: float,
+    mmw: float = 2.33,
+    plotting: bool = False,
+) -> float:
+    """
+    Function to find the base of the cloud deck by intersecting the
+    P-T profile with the saturation vapor pressure.
+
+    Parameters
+    ----------
+    composition : str
+        Cloud composition ('Fe', 'MgSiO3', 'Mg2SiO4', 'Al2O3',
+        'Na2S', or 'KCL').
+    press : np.ndarray
+        Pressures (bar).
+    temp : np.ndarray
+        Temperatures (K).
+    metallicity : float
+        Metallicity [Fe/H].
+    c_o_ratio : float
+        Carbon-to-oxygen ratio.
+    mmw : float
+        Mean molecular weight.
+    plotting : bool
+        Create a plot.
+
+    Returns
+    -------
+    float
+        Pressure (bar) at the base of the cloud deck.
+    """
+
+    Tcond_on_input_grid = get_condensation_curve(
+        composition=composition,
+        press=press,
+        metallicity=metallicity,
+        c_o_ratio=c_o_ratio,
+        mmw=mmw)
 
     Tdiff = Tcond_on_input_grid - temp
     diff_vec = Tdiff[1:] * Tdiff[:-1]