diff --git a/Project.toml b/Project.toml
index c7a7459..186dc73 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,9 +4,9 @@ authors = ["Climate Modeling Alliance"]
 version = "0.1.5"
 
 [deps]
-Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
-Dierckx = "0.5"
+Interpolations = "0.14"
 julia = "1"
diff --git a/docs/Project.toml b/docs/Project.toml
index 3c2fcdf..f060804 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,10 +1,10 @@
 [deps]
 AtmosphericProfilesLibrary = "86bc3604-9858-485a-bdbe-831ec50de11d"
-Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [compat]
-Dierckx = "0.5"
+Interpolations = "0.14"
diff --git a/docs/src/define_save_plots.jl b/docs/src/define_save_plots.jl
index 733060e..63f473b 100644
--- a/docs/src/define_save_plots.jl
+++ b/docs/src/define_save_plots.jl
@@ -1,5 +1,4 @@
 import Plots
-import Dierckx
 
 # https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
 ENV["GKSwstype"] = "nul"
diff --git a/docs/src/index.md b/docs/src/index.md
index f4e029f..e897fa0 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -14,7 +14,7 @@ This repository is meant to serve as a library of static functions in efforts to
    - `prof = z->2*z` callable by `prof(z)`
    - `prof = t->2*t` callable by `prof(t)`
    - `prof = (t,z) -> z*t` callable by `prof(t, z)`
-   - `prof = Dierckx.Spline1D(FT[0, 1, 2], FT[0, 2, 4]; k = 1)` callable by `prof(z)`
+   - `prof = Interpolations.interpolate(FT[0, 1, 2], FT[0, 2, 4])` callable by `prof(z)`
 
 ## How to contribute
 
diff --git a/src/AtmosphericProfilesLibrary.jl b/src/AtmosphericProfilesLibrary.jl
index 1a8427a..60d3625 100644
--- a/src/AtmosphericProfilesLibrary.jl
+++ b/src/AtmosphericProfilesLibrary.jl
@@ -1,6 +1,6 @@
 module AtmosphericProfilesLibrary
 
-import Dierckx
+import Interpolations as Intp
 
 # Large data-based profiles
 include("profiles/Soares.jl")
diff --git a/src/profiles/ARM_SGP.jl b/src/profiles/ARM_SGP.jl
index fabfbb9..75fb725 100644
--- a/src/profiles/ARM_SGP.jl
+++ b/src/profiles/ARM_SGP.jl
@@ -7,7 +7,7 @@ end
 function ARM_SGP_θ_liq_ice(::Type{FT}) where {FT}
     z_in = ARM_SGP_z(FT)
     θ_liq_ice_in = FT[299.0, 301.5, 302.5, 303.53, 303.7, 307.13, 314.0, 343.2] # K
-    profile = Dierckx.Spline1D(z_in, θ_liq_ice_in; k = 1)
+    profile = Intp.interpolate((z_in, ), θ_liq_ice_in, Intp.Gridded(Intp.Linear()))
     return profile
 end
 """ [Brown2002](@cite) """
@@ -18,7 +18,7 @@ function ARM_SGP_q_tot(::Type{FT}) where {FT}
     z_in = ARM_SGP_z(FT)
     r_in = FT[15.2,15.17,14.98,14.8,14.7,13.5,3.0,3.0] ./ 1000 # qt should be in kg/kg
     q_tot_in = r_in ./ (1 .+ r_in)
-    profile = Dierckx.Spline1D(z_in, q_tot_in; k = 1)
+    profile = Intp.interpolate((z_in, ), q_tot_in, Intp.Gridded(Intp.Linear()))
     return profile
 end
 """ [Brown2002](@cite) """
@@ -49,7 +49,7 @@ function ARM_SGP_tke_prescribed(::Type{FT}) where {FT}
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-    not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
+    not_type_stable_spline = Intp.interpolate((z_in, ), tke_in, Intp.Gridded(Intp.Linear()))
     return x -> FT(not_type_stable_spline(x))
 end
 
@@ -66,8 +66,8 @@ function ARM_SGP_dTdt(::Type{FT}) where {FT}
     AT_in = FT[0.0, 0.0, 0.0, -0.08, -0.016, -0.016] ./ 3600
     # Radiative forcing for theta [K/h] converted to [K/sec]
     RT_in = FT[-0.125, 0.0, 0.0, 0.0, 0.0, -0.1] ./ 3600
-    dTdt_A = Dierckx.Spline1D(t_in, AT_in; k = 1)
-    dTdt_R = Dierckx.Spline1D(t_in, RT_in; k = 1)
+    dTdt_A = Intp.interpolate((t_in, ), AT_in, Intp.Gridded(Intp.Linear()))
+    dTdt_R = Intp.interpolate((t_in, ), RT_in, Intp.Gridded(Intp.Linear()))
     return (t, z) -> if z <= 1000.0
         dTdt_A(t)+dTdt_R(t)
     elseif z > 1000.0 && z <= 2000.0
@@ -82,7 +82,7 @@ function ARM_SGP_dqtdt(::Type{FT}) where {FT}
     t_in = ARM_SGP_time(FT)
     # Radiative forcing for qt converted to [kg/kg/sec]
     Rqt_in = FT[0.08, 0.02, 0.04, -0.1, -0.16, -0.3] ./ 1000 ./ 3600
-    dqtdt = Dierckx.Spline1D(t_in, Rqt_in; k = 1)
+    dqtdt = Intp.interpolate((t_in, ), Rqt_in, Intp.Gridded(Intp.Linear()))
     return (Π, t, z) -> if z <= 1000.0
         dqtdt(t) * Π
     elseif z > 1000.0 && z <= 2000.0
@@ -96,7 +96,7 @@ end
 function ARM_SGP_shf(::Type{FT}) where {FT}
     t_Sur_in = FT[0.0, 4.0, 6.5, 7.5, 10.0, 12.5, 14.5] .* 3600 #LES time is in sec
     shf = FT[-30.0, 90.0, 140.0, 140.0, 100.0, -10, -10] # W/m^2
-    profile = Dierckx.Spline1D(t_Sur_in, shf; k = 1)
+    profile = Intp.interpolate((t_Sur_in, ), shf, Intp.Gridded(Intp.Linear()))
     return profile
 end
 
@@ -104,6 +104,6 @@ end
 function ARM_SGP_lhf(::Type{FT}) where {FT}
     t_Sur_in = FT[0.0, 4.0, 6.5, 7.5, 10.0, 12.5, 14.5] .* 3600 #LES time is in sec
     lhf = FT[5.0, 250.0, 450.0, 500.0, 420.0, 180.0, 0.0] # W/m^2
-    profile = Dierckx.Spline1D(t_Sur_in, lhf; k = 1)
+    profile = Intp.interpolate((t_Sur_in, ), lhf, Intp.Gridded(Intp.Linear()))
     return profile
 end
diff --git a/src/profiles/Bomex.jl b/src/profiles/Bomex.jl
index 9764bb5..466c17f 100644
--- a/src/profiles/Bomex.jl
+++ b/src/profiles/Bomex.jl
@@ -45,7 +45,7 @@ function Bomex_tke_prescribed(::Type{FT}) where {FT}
                 0.1468, 0.1470, 0.1458, 0.1423, 0.1351, 0.1224, 0.1033, 0.0793, 0.0542, 0.0322,
                 0.0162, 0.0068, 0.0024, 0.0007, 0.0001, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-    not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
+    not_type_stable_spline = Intp.interpolate((z_in, ), tke_in, Intp.Gridded(Intp.Linear()))
     return x -> FT(not_type_stable_spline(x))
 end
 
diff --git a/src/profiles/DryBubble.jl b/src/profiles/DryBubble.jl
index 54219fb..fe8d301 100644
--- a/src/profiles/DryBubble.jl
+++ b/src/profiles/DryBubble.jl
@@ -61,7 +61,7 @@ function DryBubble_θ_liq_ice(::Type{FT}) where {FT}
         299.9837, 299.9837, 299.9837, 299.9837, 299.9837, 299.9837,
         299.9837, 299.9837
     ]
-    profile = Dierckx.Spline1D(z_in, θ_liq_ice_in; k = 1)
+    profile = Intp.interpolate((z_in, ), θ_liq_ice_in, Intp.Gridded(Intp.Linear()))
     return profile
 end
 
@@ -97,7 +97,7 @@ function DryBubble_updrafts_θ_liq_ice(::Type{FT}) where {FT}
         300.5212, 300.482 , 300.4272, 300.3875, 300.3354, 300.2968,
         300.2587, 300.2216, 300.1782, 300.1452, 300.1143, 300.0859,
         300.0603, 300.0408, 300.0211, 300.0067, 299.9963, 299.9884]
-    profile = Dierckx.Spline1D(z_in, θ_liq_in; k = 1)
+    profile = Intp.interpolate((z_in, ), θ_liq_in, Intp.Gridded(Intp.Linear()))
     return profile
 end
 
@@ -114,7 +114,7 @@ function DryBubble_updrafts_area(::Type{FT}) where {FT}
         0.15 , 0.145, 0.145, 0.145, 0.14 , 0.14 , 0.14 , 0.135, 0.135,
         0.13 , 0.13 , 0.125, 0.12 , 0.115, 0.115, 0.11 , 0.105, 0.1  ,
         0.095, 0.085, 0.08 , 0.07 , 0.055, 0.04]
-    profile = Dierckx.Spline1D(z_in, Area_in; k = 1)
+    profile = Intp.interpolate((z_in, ), Area_in, Intp.Gridded(Intp.Linear()))
     return profile
 end
 
@@ -132,7 +132,7 @@ function DryBubble_updrafts_w(::Type{FT}) where {FT}
         0.1601, 0.1559, 0.1589, 0.1543, 0.1496, 0.1514, 0.1464, 0.1475,
         0.1422, 0.1425, 0.1424, 0.1419, 0.1361, 0.135 , 0.1335, 0.1316,
         0.1294, 0.1302, 0.1271, 0.1264, 0.1269, 0.1256]
-    profile = Dierckx.Spline1D(z_in, W_in; k = 1)
+    profile = Intp.interpolate((z_in, ), W_in, Intp.Gridded(Intp.Linear()))
     return profile
 end
 
@@ -153,6 +153,6 @@ function DryBubble_updrafts_T(::Type{FT}) where {FT}
         270.4339, 269.9094, 269.3711, 268.8465, 268.311 , 267.7877,
         267.2649, 266.7432, 266.2159, 265.698 , 265.1821, 264.6685,
         264.1574, 263.6518, 263.1461, 262.6451, 262.1476, 261.6524]
-    profile = Dierckx.Spline1D(z_in, T_in; k = 1)
+    profile = Intp.interpolate((z_in, ), T_in, Intp.Gridded(Intp.Linear()))
     return profile
 end
diff --git a/src/profiles/Dycoms_RF01.jl b/src/profiles/Dycoms_RF01.jl
index 6024d76..18a9593 100644
--- a/src/profiles/Dycoms_RF01.jl
+++ b/src/profiles/Dycoms_RF01.jl
@@ -30,7 +30,7 @@ function Dycoms_RF01_tke_prescribed(::Type{FT}) where {FT}
               775.0, 825.0, 875.0, 925.0, 975.0, 1025.0, 1075.0, 1125.0, 1175.0, 1225.0, 1275.0, 1325.0, 1375.0, 1425.0, 1475.0]
     tke_in = FT[0.2726, 0.5479, 0.6597, 0.7079, 0.7285, 0.7343, 0.7319, 0.7252, 0.7166, 0.7064, 0.6887, 0.6317,
                 0.6362, 0.6266, 0.5832, 0.4633, 0.0504, 0.0001, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-    not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
+    not_type_stable_spline = Intp.interpolate((z_in, ), tke_in, Intp.Gridded(Intp.Linear()))
     return x -> FT(not_type_stable_spline(x))
 end
 
diff --git a/src/profiles/GABLS.jl b/src/profiles/GABLS.jl
index 81b54d5..3e6e4d5 100644
--- a/src/profiles/GABLS.jl
+++ b/src/profiles/GABLS.jl
@@ -26,7 +26,7 @@ GABLS_tke(::Type{FT}) where {FT} =
 function GABLS_tke_prescribed(::Type{FT}) where {FT}
     z_in = FT[25.0, 75.0, 125.0, 175.0, 225.0, 275.0, 325.0, 375.0]
     tke_in = FT[0.4662, 0.3873, 0.2777, 0.0277, 0.0003, 5.89e-8, 0.0, 0.0]
-    not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
+    not_type_stable_spline = Intp.interpolate((z_in, ), tke_in, Intp.Gridded(Intp.Linear()))
     return x -> FT(not_type_stable_spline(x))
 end
 
diff --git a/src/profiles/GATE_III.jl b/src/profiles/GATE_III.jl
index 59a64a9..7708d5e 100644
--- a/src/profiles/GATE_III.jl
+++ b/src/profiles/GATE_III.jl
@@ -15,7 +15,7 @@ function GATE_III_q_tot(::Type{FT}) where {FT}
                0.9, 0.5, 0.25, 0.125, 0.065, 0.003, 0.0015, 0.0007,  0.0003,  0.0001,  0.0001,  0.0001,  0.0001,
                0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001, 0.0001] ./ 1000 # mixing ratio should be in kg/kg
     q_tot = r_in ./ (1 .+ r_in) # convert mixing ratio to specific humidity
-    profile = Dierckx.Spline1D(z_in, q_tot; k = 1)
+    profile = Intp.interpolate((z_in, ), q_tot, Intp.Gridded(Intp.Linear()))
     return profile
 end
 
@@ -26,7 +26,7 @@ function GATE_III_u(::Type{FT}) where {FT}
     U_in = FT[  -1, -1.75, -2.5, -3.6, -6.0, -8.75, -11.75, -13.0, -13.1, -12.1, -11.0, -8.5, -5.0, -2.6, 0.0,
                0.5, 0.4,  0.3,   0.0,  -1.0, -2.5,   -3.5,   -4.5, -4.8, -5.0, -3.5, -2.0, -1.0, -1.0, -1.0,
                -1.5, -2.0, -2.5, -2.6, -2.7, -3.0, -3.0, -3.0] # [m/s]
-    profile = Dierckx.Spline1D(z_in, U_in; k = 1)
+    profile = Intp.interpolate((z_in, ), U_in, Intp.Gridded(Intp.Linear()))
     return profile
 end
 
@@ -36,7 +36,7 @@ function GATE_III_T(::Type{FT}) where {FT}
     z_in = GATE_III_z(FT)
     T_in = FT[299.184, 294.836, 294.261, 288.773, 276.698, 265.004, 253.930, 243.662, 227.674, 214.266, 207.757, 201.973, 198.278, 197.414, 198.110, 198.110]
     z_T_in = FT[0.0, 0.492, 0.700, 1.698, 3.928, 6.039, 7.795, 9.137, 11.055, 12.645, 13.521, 14.486, 15.448, 16.436, 17.293, 22.0] .* 1000 # for km
-    profile = Dierckx.Spline1D(z_T_in, T_in; k = 1)
+    profile = Intp.interpolate((z_T_in, ), T_in, Intp.Gridded(Intp.Linear()))
     return profile
 end
 
@@ -77,7 +77,7 @@ function GATE_III_dqtdt(::Type{FT}) where {FT}
                    0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0] ./ (24 * 3600) ./ 1000
 
     Qtend_in = r_tend_in ./ (1 .+ r_tend_in) # convert mixing ratio to specific humidity
-    profile = Dierckx.Spline1D(z_in, Qtend_in; k = 1)
+    profile = Intp.interpolate((z_in, ), Qtend_in, Intp.Gridded(Intp.Linear()))
     return profile
 end
 
@@ -97,7 +97,7 @@ function GATE_III_dTdt(::Type{FT}) where {FT}
                   -3.7, -3.5, -3.25, -3.0, -2.8,  -2.5, -2.1,  -1.7,  -1.3,   -1.0,   -0.7,   -0.5, -0.4,
                   -0.3, -0.2, -0.1,-0.05,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0] ./ (24*3600)
 
-    profile_T = Dierckx.Spline1D(z_in, Ttend_in; k = 1)
-    profile_R = Dierckx.Spline1D(z_in, RAD_in; k = 1)
+    profile_T = Intp.interpolate((z_in, ), Ttend_in, Intp.Gridded(Intp.Linear()))
+    profile_R = Intp.interpolate((z_in, ), RAD_in, Intp.Gridded(Intp.Linear()))
     return z -> profile_T(z) + profile_R(z)
 end
diff --git a/src/profiles/LifeCycleTan2018.jl b/src/profiles/LifeCycleTan2018.jl
index d1e8381..c9e490b 100644
--- a/src/profiles/LifeCycleTan2018.jl
+++ b/src/profiles/LifeCycleTan2018.jl
@@ -57,7 +57,7 @@ function LifeCycleTan2018_tke_prescribed(::Type{FT}) where {FT}
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0]
-    not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
+    not_type_stable_spline = Intp.interpolate((z_in, ), tke_in, Intp.Gridded(Intp.Linear()))
     return x -> FT(not_type_stable_spline(x))
 end
 
diff --git a/src/profiles/Nieuwstadt.jl b/src/profiles/Nieuwstadt.jl
index 8d2fa53..b450aef 100644
--- a/src/profiles/Nieuwstadt.jl
+++ b/src/profiles/Nieuwstadt.jl
@@ -51,6 +51,6 @@ function Nieuwstadt_tke_prescribed(::Type{FT}) where {FT}
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-    not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
+    not_type_stable_spline = Intp.interpolate((z_in, ), tke_in, Intp.Gridded(Intp.Linear()))
     return x -> FT(not_type_stable_spline(x))
 end
diff --git a/src/profiles/Rico.jl b/src/profiles/Rico.jl
index 216338e..8f96150 100644
--- a/src/profiles/Rico.jl
+++ b/src/profiles/Rico.jl
@@ -78,6 +78,6 @@ function Rico_tke_prescribed(::Type{FT}) where {FT}
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0]
-    not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
+    not_type_stable_spline = Intp.interpolate((z_in, ), tke_in, Intp.Gridded(Intp.Linear()))
     return x -> FT(not_type_stable_spline(x))
 end
diff --git a/src/profiles/Soares.jl b/src/profiles/Soares.jl
index afbfb07..5520492 100644
--- a/src/profiles/Soares.jl
+++ b/src/profiles/Soares.jl
@@ -46,6 +46,6 @@ function Soares_tke_prescribed(::Type{FT}) where {FT}
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-    not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
+    not_type_stable_spline = Intp.interpolate((z_in, ), tke_in, Intp.Gridded(Intp.Linear()))
     return x -> FT(not_type_stable_spline(x))
 end
diff --git a/src/profiles/TRMM_LBA.jl b/src/profiles/TRMM_LBA.jl
index c3fa9c8..8b4cb39 100644
--- a/src/profiles/TRMM_LBA.jl
+++ b/src/profiles/TRMM_LBA.jl
@@ -21,7 +21,7 @@ function TRMM_LBA_p(::Type{FT}) where {FT}
               182.3, 167.9, 154.9, 143.0, 131.1, 119.7, 108.9,
               100.1,  92.1,  84.6,  77.5,  71.4,  65.9,  60.7,
                55.9,  51.3,  47.2,  43.3,  10.3] .* 100 # LES pres is in pasc
-    return Dierckx.Spline1D(z_in, p_in; k = 1)
+    return Intp.interpolate((z_in, ), p_in, Intp.Gridded(Intp.Linear()))
 end
 """ [Grabowski2006](@cite) """
 function TRMM_LBA_T(::Type{FT}) where {FT}
@@ -33,7 +33,7 @@ function TRMM_LBA_T(::Type{FT}) where {FT}
              -59.16, -63.60, -67.68, -70.77, -74.41, -77.51, -80.64,
              -80.69, -80.00, -81.38, -81.17, -78.32, -74.77, -74.52,
              -72.62, -70.87, -69.19, -66.90, -66.90] .+ FT(273.15) # LES T is in deg K
-    return Dierckx.Spline1D(z_in, T_in; k = 1)
+    return Intp.interpolate((z_in, ), T_in, Intp.Gridded(Intp.Linear()))
 end
 """ [Grabowski2006](@cite) """
 function TRMM_LBA_RH(::Type{FT}) where {FT}
@@ -45,7 +45,7 @@ function TRMM_LBA_RH(::Type{FT}) where {FT}
                45.33,  39.78,  33.78,  28.78,  24.67,  20.67,  17.67,
                17.11,  16.22,  14.22,  13.00,  13.00,  12.22,   9.56,
                 7.78,   5.89,   4.33,   3.00,   3.00]
-    return Dierckx.Spline1D(z_in, RH_in; k = 1)
+    return Intp.interpolate((z_in, ), RH_in, Intp.Gridded(Intp.Linear()))
 end
 """ [Grabowski2006](@cite) """
 function TRMM_LBA_u(::Type{FT}) where {FT}
@@ -57,7 +57,7 @@ function TRMM_LBA_u(::Type{FT}) where {FT}
              -9.00,  -7.77,  -5.37,  -3.88,  -1.15,  -2.36,  -9.20,
              -8.01,  -5.68,  -8.83, -14.51, -15.55, -15.36, -17.67,
             -17.82, -18.94, -15.92, -15.32, -15.32]
-    return Dierckx.Spline1D(z_in, u_in; k = 1)
+    return Intp.interpolate((z_in, ), u_in, Intp.Gridded(Intp.Linear()))
 end
 """ [Grabowski2006](@cite) """
 function TRMM_LBA_v(::Type{FT}) where {FT}
@@ -69,7 +69,7 @@ function TRMM_LBA_v(::Type{FT}) where {FT}
                3.14,   3.93,   7.57,   2.58,   2.50,   6.44,   6.84,
                0.19,  -2.20,  -3.60,   0.56,   6.68,   9.41,   7.03,
                5.32,   1.14,  -0.65,   5.27,   5.27]
-    return Dierckx.Spline1D(z_in, v_in; k = 1)
+    return Intp.interpolate((z_in, ), v_in, Intp.Gridded(Intp.Linear()))
 end
 
 """ [Grabowski2006](@cite) """
@@ -105,7 +105,7 @@ function TRMM_LBA_tke_prescribed(::Type{FT}) where {FT}
        1.7938 , 1.56451, 1.37531, 1.17515, 0.96797, 0.61262, 0.26423,
        0.14929, 0.07465, 0.00635, 0.     , 0.     , 0.     , 0.     ,
        0.     , 0.     , 0.     , 0.     , 0.     ]
-    not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
+    not_type_stable_spline = Intp.interpolate((z_in, ), tke_in, Intp.Gridded(Intp.Linear()))
     return x -> FT(not_type_stable_spline(x))
 end
 
@@ -229,6 +229,6 @@ function TRMM_LBA_radiation(::Type{FT}) where {FT}
     rad_in = reduce(hcat, rad_in)::Matrix{FT}
     rad_in .= rad_in ./ 86400
     rad_in = (rad_in')::LinearAlgebra.Adjoint{FT, Matrix{FT}}
-    profile = Dierckx.Spline2D(rad_time, z_in, rad_in; kx = 1, ky = 1)
+    profile = Intp.interpolate((rad_time, z_in), rad_in, (Intp.Gridded(Intp.Linear()), Intp.Gridded(Intp.Linear())))
     return profile
 end