diff --git a/Project.toml b/Project.toml
index 21514968..98fcb8c5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,20 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
+# Still need this- not sure why https://pkgdocs.julialang.org/v1/creating-packages/#Using-an-extension-while-supporting-older-Julia-versions
+CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
+
+[weakdeps]
+CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
+
+[extras]
+CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
+
+[targets]
+test = ["CLIMAParameters"]
+
+[extensions]
+CreateParametersExt = "CLIMAParameters"
 
 [compat]
 DocStringExtensions = "0.8.1, 0.9"
@@ -15,3 +29,4 @@ KernelAbstractions = "0.7, 0.8, 0.9"
 Random = "1"
 RootSolvers = "0.4"
 julia = "1.6"
+CLIMAParameters = "0.7"
diff --git a/docs/src/Clausius_Clapeyron.jl b/docs/src/Clausius_Clapeyron.jl
index 79c56897..4c2c0909 100644
--- a/docs/src/Clausius_Clapeyron.jl
+++ b/docs/src/Clausius_Clapeyron.jl
@@ -4,23 +4,6 @@ import CLIMAParameters as CP
 using Thermodynamics.TestedProfiles
 import Thermodynamics.Parameters as TP
 
-function get_parameter_set(::Type{FT}) where {FT}
-    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
-    aliases = string.(fieldnames(TP.ThermodynamicsParameters))
-    param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
-    param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
-    logfilepath = joinpath(@__DIR__, "logfilepath_$FT.toml")
-    CP.log_parameter_information(toml_dict, logfilepath)
-    return param_set
-end
-const param_set_Float64 = get_parameter_set(Float64)
-const param_set_Float32 = get_parameter_set(Float32)
-parameter_set(::Type{Float64}) = param_set_Float64
-parameter_set(::Type{Float32}) = param_set_Float32
-ArrayType = Array{Float64}
-
-FT = eltype(ArrayType)
-param_set = parameter_set(FT)
 profiles = TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
 (; T, p, e_int, ρ, θ_liq_ice, phase_type) = profiles
 (; q_tot, q_liq, q_ice, q_pt, RH, e_kin, e_pot) = profiles
diff --git a/ext/CreateParametersExt.jl b/ext/CreateParametersExt.jl
new file mode 100644
index 00000000..32a47ede
--- /dev/null
+++ b/ext/CreateParametersExt.jl
@@ -0,0 +1,41 @@
+module CreateParametersExt
+
+import Thermodynamics.Parameters.ThermodynamicsParameters
+import CLIMAParameters as CP
+
+ThermodynamicsParameters(FT) = ThermodynamicsParameters(CP.create_toml_dict(FT))
+
+function ThermodynamicsParameters(toml_dict::CP.AbstractTOMLDict)
+    name_map = Dict{String, String}(
+        "temperature_min_at_reference" => "T_min_ref",
+        "entropy_water_vapor" => "entropy_water_vapor",
+        "entropy_dry_air" => "entropy_dry_air",
+        "potential_temperature_reference_pressure" => "p_ref_theta",
+        "entropy_reference_temperature" => "entropy_reference_temperature",
+        "temperature_saturation_adjustment_max" => "T_max",
+        "molar_mass_dry_air" => "molmass_dryair",
+        "pow_icenuc" => "pow_icenuc",
+        "temperature_triple_point" => "T_triple",
+        "adiabatic_exponent_dry_air" => "kappa_d",
+        "pressure_triple_point" => "press_triple",
+        "thermodynamics_temperature_reference" => "T_0",
+        "temperature_water_freeze" => "T_freeze",
+        "isobaric_specific_heat_ice" => "cp_i",
+        "latent_heat_sublimation_at_reference" => "LH_s0",
+        "isobaric_specific_heat_vapor" => "cp_v",
+        "molar_mass_water" => "molmass_water",
+        "mean_sea_level_pressure" => "MSLP",
+        "isobaric_specific_heat_liquid" => "cp_l",
+        "latent_heat_vaporization_at_reference" => "LH_v0",
+        "temperature_saturation_adjustment_min" => "T_min",
+        "gas_constant" => "gas_constant",
+        "temperature_mean_at_reference" => "T_surf_ref",
+        "gravitational_acceleration" => "grav",
+        "temperature_homogenous_nucleation" => "T_icenuc",
+    )
+    parameters = CP.get_parameter_values(toml_dict, name_map, "Thermodynamics")
+    FT = CP.float_type(toml_dict)
+    return ThermodynamicsParameters{FT}(; parameters...)
+end
+
+end
diff --git a/perf/common_micro_bm.jl b/perf/common_micro_bm.jl
index df9c3c19..064ff01c 100644
--- a/perf/common_micro_bm.jl
+++ b/perf/common_micro_bm.jl
@@ -10,17 +10,7 @@ using Test
 import UnPack
 import BenchmarkTools
 
-import CLIMAParameters
-const CP = CLIMAParameters
-function get_parameter_set(::Type{FT}) where {FT}
-    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
-    aliases = string.(fieldnames(TP.ThermodynamicsParameters))
-    param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
-    param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
-    # logfilepath = joinpath(@__DIR__, "logfilepath_$FT.toml")
-    # CP.log_parameter_information(toml_dict, logfilepath)
-    return param_set
-end
+import CLIMAParameters as CP
 
 #####
 ##### Finding indexes in profiles satisfying certain conditions
diff --git a/perf/jet.jl b/perf/jet.jl
index b47d7cf0..739a8024 100644
--- a/perf/jet.jl
+++ b/perf/jet.jl
@@ -5,7 +5,7 @@ include("common_micro_bm.jl")
 TD.print_warning() = false
 
 function jet_thermo_states(::Type{FT}) where {FT}
-    param_set = get_parameter_set(FT)
+    param_set = TP.ThermodynamicsParameters(FT)
     ArrayType = Array{FT}
     profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
 
diff --git a/perf/microbenchmarks.jl b/perf/microbenchmarks.jl
index 254078d6..157d9708 100644
--- a/perf/microbenchmarks.jl
+++ b/perf/microbenchmarks.jl
@@ -3,7 +3,7 @@ include("common_micro_bm.jl")
 function benchmark_thermo_states(::Type{FT}) where {FT}
     summary = OrderedCollections.OrderedDict()
     ArrayType = Array{FT}
-    param_set = get_parameter_set(FT)
+    param_set = TP.ThermodynamicsParameters(FT)
     profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
 
     for C in (
diff --git a/src/Thermodynamics.jl b/src/Thermodynamics.jl
index 63855122..cee3904e 100644
--- a/src/Thermodynamics.jl
+++ b/src/Thermodynamics.jl
@@ -88,4 +88,7 @@ include("TestedProfiles.jl")
 Base.broadcastable(dap::DryAdiabaticProcess) = tuple(dap)
 Base.broadcastable(phase::Phase) = tuple(phase)
 
+# For backwards compatibility with package extensions
+isdefined(Base, :get_extension) || include(joinpath("..", "ext", "CreateParametersExt.jl"))
+
 end #module Thermodynamics.jl
diff --git a/test/TemperatureProfiles.jl b/test/TemperatureProfiles.jl
index 96555fa6..bc845e09 100644
--- a/test/TemperatureProfiles.jl
+++ b/test/TemperatureProfiles.jl
@@ -1,31 +1,14 @@
 using Test
-import Thermodynamics
-const TD = Thermodynamics
-const TP = TD.Parameters
-const TDTP = TD.TemperatureProfiles
+import Thermodynamics as TD
+import Thermodynamics.Parameters as TP
+import Thermodynamics.TemperatureProfiles as TDTP
 using ForwardDiff
 
-import CLIMAParameters
-const CP = CLIMAParameters
-
-function get_parameter_set(::Type{FT}) where {FT}
-    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
-    aliases = string.(fieldnames(TP.ThermodynamicsParameters))
-    param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
-    param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
-    logfilepath = joinpath(@__DIR__, "logfilepath_$FT.toml")
-    CP.log_parameter_information(toml_dict, logfilepath)
-    return param_set
-end
-
-const param_set_Float64 = get_parameter_set(Float64)
-const param_set_Float32 = get_parameter_set(Float32)
-parameter_set(::Type{Float64}) = param_set_Float64
-parameter_set(::Type{Float32}) = param_set_Float32
+import CLIMAParameters as CP
 
 @testset "TemperatureProfiles - DecayingTemperatureProfile" begin
     for FT in [Float32, Float64]
-        param_set = parameter_set(FT)
+        param_set = TP.ThermodynamicsParameters(FT)
         _grav = FT(TP.grav(param_set))
         _R_d = FT(TP.R_d(param_set))
         _p_ref = FT(TP.MSLP(param_set))
diff --git a/test/data_tests.jl b/test/data_tests.jl
index 03391581..1b5e03bb 100644
--- a/test/data_tests.jl
+++ b/test/data_tests.jl
@@ -16,7 +16,7 @@ dycoms_dataset_path = get_data_folder(dycoms_dataset)
 
 @testset "Data tests" begin
     FT = Float64
-    param_set = parameter_set(FT)
+    param_set = TP.ThermodynamicsParameters(FT)
     data = joinpath(dycoms_dataset_path, "test_data_PhaseEquil.nc")
     ds_PhaseEquil = Dataset(data, "r")
     e_int = Array{FT}(ds_PhaseEquil["e_int"][:])
@@ -28,7 +28,7 @@ dycoms_dataset_path = get_data_folder(dycoms_dataset)
 end
 
 @testset "pθq data-driven tests" begin
-    param_set = parameter_set(Float64)
+    param_set = TP.ThermodynamicsParameters(Float64)
     #! format: off
     pθq_broken = [
         # (; p = , θ_liq_ice = , q_tot = ),
diff --git a/test/relations.jl b/test/relations.jl
index 6115e45a..10ebe92d 100644
--- a/test/relations.jl
+++ b/test/relations.jl
@@ -1,10 +1,10 @@
 using Test
 using Thermodynamics
+import Thermodynamics as TD
+import Thermodynamics.Parameters as TP
 using Thermodynamics.TemperatureProfiles
 using Thermodynamics.TestedProfiles
 
-const TD = Thermodynamics
-const TP = TD.Parameters
 using UnPack
 using NCDatasets
 using Random
@@ -15,28 +15,11 @@ const RS = RootSolvers
 using LinearAlgebra
 import ForwardDiff
 
-const TP = TD.Parameters
-
 import CLIMAParameters
 const CP = CLIMAParameters
 
-function get_parameter_set(::Type{FT}) where {FT}
-    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
-    aliases = string.(fieldnames(TP.ThermodynamicsParameters))
-    param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
-    param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
-    logfilepath = joinpath(@__DIR__, "logfilepath_$FT.toml")
-    # CP.log_parameter_information(toml_dict, logfilepath)
-    return param_set
-end
-
-const param_set_Float64 = get_parameter_set(Float64)
-const param_set_Float32 = get_parameter_set(Float32)
-parameter_set(::Type{Float64}) = param_set_Float64
-parameter_set(::Type{Float32}) = param_set_Float32
-
-
 # Tolerances for tested quantities:
+param_set_Float64 = TP.ThermodynamicsParameters(Float64)
 atol_temperature = 5e-1
 atol_energy = TP.cv_d(param_set_Float64) * atol_temperature
 rtol_temperature = 1e-1
@@ -63,7 +46,7 @@ compare_moisture(param_set, ts::PhaseNonEquil, q_pt::PhasePartition) = all((
 @testset "Thermodynamics - isentropic processes" begin
     for ArrayType in array_types
         FT = eltype(ArrayType)
-        param_set = parameter_set(FT)
+        param_set = TP.ThermodynamicsParameters(FT)
 
         _R_d = FT(TP.R_d(param_set))
         _molmass_ratio = FT(TP.molmass_ratio(param_set))
@@ -120,7 +103,7 @@ end
 
 @testset "Thermodynamics - correctness" begin
     FT = Float64
-    param_set = parameter_set(FT)
+    param_set = TP.ThermodynamicsParameters(FT)
     _R_d = FT(TP.R_d(param_set))
     _molmass_ratio = FT(TP.molmass_ratio(param_set))
     _cp_d = FT(TP.cp_d(param_set))
@@ -523,7 +506,7 @@ end
     or(a, b) = a || b
     for ArrayType in array_types
         FT = eltype(ArrayType)
-        param_set = parameter_set(FT)
+        param_set = TP.ThermodynamicsParameters(FT)
         profiles = TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
         @unpack T, p, e_int, ρ, θ_liq_ice, phase_type = profiles
         @unpack q_tot, q_liq, q_ice, q_pt, RH, e_kin, e_pot = profiles
@@ -899,7 +882,7 @@ end
 
     ArrayType = Array{Float64}
     FT = eltype(ArrayType)
-    param_set = parameter_set(FT)
+    param_set = TP.ThermodynamicsParameters(FT)
     profiles = TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
     @unpack T, p, e_int, ρ, θ_liq_ice, phase_type = profiles
     @unpack q_tot, q_liq, q_ice, q_pt, RH, e_kin, e_pot = profiles
@@ -1008,7 +991,7 @@ end
 
     for ArrayType in array_types
         FT = eltype(ArrayType)
-        param_set = parameter_set(FT)
+        param_set = TP.ThermodynamicsParameters(FT)
 
         profiles = TestedProfiles.PhaseDryProfiles(param_set, ArrayType)
         @unpack T, p, e_int, h, ρ, θ_liq_ice, phase_type = profiles
@@ -1368,7 +1351,7 @@ end
     # with converging to the same tolerances as `Float64`, so they're relaxed here.
     ArrayType = Array{Float32}
     FT = eltype(ArrayType)
-    param_set = parameter_set(FT)
+    param_set = TP.ThermodynamicsParameters(FT)
 
     profiles = TestedProfiles.PhaseDryProfiles(param_set, ArrayType)
     @unpack T, p, e_int, ρ, θ_liq_ice, phase_type = profiles
@@ -1509,7 +1492,7 @@ end
 
     ArrayType = Array{Float64}
     FT = eltype(ArrayType)
-    param_set = parameter_set(FT)
+    param_set = TP.ThermodynamicsParameters(FT)
     profiles = TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
     @unpack T, p, e_int, ρ, θ_liq_ice, phase_type = profiles
     @unpack q_tot, q_liq, q_ice, q_pt, RH, e_kin, e_pot = profiles
@@ -1677,7 +1660,7 @@ end
 @testset "Thermodynamics - ProfileSet Iterator" begin
     ArrayType = Array{Float64}
     FT = eltype(ArrayType)
-    param_set = parameter_set(FT)
+    param_set = TP.ThermodynamicsParameters(FT)
     profiles = TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
     @unpack T, q_pt, z, phase_type = profiles
     @test all(z .≈ (nt.z for nt in profiles))
@@ -1697,7 +1680,7 @@ end
 @testset "Thermodynamics - test T_guess" begin
     ArrayType = Array{Float64}
     FT = eltype(ArrayType)
-    param_set = parameter_set(FT)
+    param_set = TP.ThermodynamicsParameters(FT)
     profiles = TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
     @unpack p, ρ, e_int, h, θ_liq_ice, q_tot, T, phase_type = profiles
     T_guess = T .+ (FT(0.2) .* randn(FT, length(T)))
@@ -1734,7 +1717,7 @@ TD.solution_type() = RS.VerboseSolution()
 @testset "Test data collection" begin
     ArrayType = Array{Float64}
     FT = eltype(ArrayType)
-    param_set = parameter_set(FT)
+    param_set = TP.ThermodynamicsParameters(FT)
     profiles = TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
     @unpack ρ, e_int, q_tot = profiles
     ts = PhaseEquil_ρeq.(param_set, ρ, e_int, q_tot)
diff --git a/test/runtests_gpu.jl b/test/runtests_gpu.jl
index fe552dc4..b19158c4 100644
--- a/test/runtests_gpu.jl
+++ b/test/runtests_gpu.jl
@@ -22,22 +22,6 @@ const CP = CLIMAParameters
 
 const TP = TD.Parameters
 
-function get_parameter_set(::Type{FT}) where {FT}
-    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
-    aliases = string.(fieldnames(TP.ThermodynamicsParameters))
-    param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
-    param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)
-    logfilepath = joinpath(@__DIR__, "logfilepath_$FT.toml")
-    CP.log_parameter_information(toml_dict, logfilepath)
-    return param_set
-end
-
-const param_set_Float64 = get_parameter_set(Float64)
-const param_set_Float32 = get_parameter_set(Float32)
-parameter_set(::Type{Float64}) = param_set_Float64
-parameter_set(::Type{Float32}) = param_set_Float32
-
-
 if get(ARGS, 1, "Array") == "CuArray"
     import CUDA
     ArrayType = CUDA.CuArray
@@ -61,6 +45,7 @@ end
     i = @index(Group, Linear)
     @inbounds begin
 
+        param_set = TP.ThermodynamicsParameters(FT)
         ts = TD.PhaseEquil_ρeq(param_set, FT(ρ[i]), FT(e_int[i]), FT(q_tot[i]))
         dst[1, i] = TD.air_temperature(param_set, ts)
 
@@ -106,7 +91,7 @@ convert_profile_set(ps::TD.TestedProfiles.ProfileSet, ArrayType, slice) =
 
 @testset "Thermodynamics - kernels" begin
     FT = Float32
-    param_set = parameter_set(FT)
+    param_set = TP.ThermodynamicsParameters(FT)
     dev = device(ArrayType)
     profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, Array)
     slice = Colon()