diff --git a/src/AtmosphericProfilesLibrary.jl b/src/AtmosphericProfilesLibrary.jl
index 1a8427a..507e98c 100644
--- a/src/AtmosphericProfilesLibrary.jl
+++ b/src/AtmosphericProfilesLibrary.jl
@@ -2,6 +2,32 @@ module AtmosphericProfilesLibrary
 
 import Dierckx
 
+abstract type AbstractProfile end
+struct TimeProfile{P} <: AbstractProfile
+  prof::P
+end
+@inline (prof::TimeProfile)(t) = prof.prof(t)
+
+struct TimeZProfile{P} <: AbstractProfile
+  prof::P
+end
+@inline (prof::TimeZProfile)(t, z) = prof.prof(t, z)
+
+struct ZProfile{P} <: AbstractProfile
+  prof::P
+end
+@inline (prof::ZProfile)(z) = prof.prof(z)
+
+struct ΠTimeZProfile{P} <: AbstractProfile
+  prof::P
+end
+@inline (prof::ΠTimeZProfile)(Π, t, z) = prof.prof(Π, t, z)
+
+struct ΠZProfile{P} <: AbstractProfile
+  prof::P
+end
+@inline (prof::ΠZProfile)(Π, z) = prof.prof(Π, z)
+
 # Large data-based profiles
 include("profiles/Soares.jl")
 include("profiles/Nieuwstadt.jl")
diff --git a/src/profiles/ARM_SGP.jl b/src/profiles/ARM_SGP.jl
index fabfbb9..251ec26 100644
--- a/src/profiles/ARM_SGP.jl
+++ b/src/profiles/ARM_SGP.jl
@@ -8,10 +8,10 @@ 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)
-    return profile
+    return ZProfile(profile)
 end
 """ [Brown2002](@cite) """
-ARM_SGP_u(::Type{FT}) where {FT} = z -> FT(10)
+ARM_SGP_u(::Type{FT}) where {FT} = ZProfile(z -> FT(10))
 
 """ [Brown2002](@cite) """
 function ARM_SGP_q_tot(::Type{FT}) where {FT}
@@ -19,15 +19,15 @@ function ARM_SGP_q_tot(::Type{FT}) where {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)
-    return profile
+    return ZProfile(profile)
 end
 """ [Brown2002](@cite) """
 function ARM_SGP_tke(::Type{FT}) where {FT}
-    return z -> if z <= 2500.0
-        FT(1) - z / 3000
-    else
-        FT(0)
-    end
+    return ZProfile(z -> if z <= 2500.0
+            FT(1) - z / 3000
+        else
+            FT(0)
+        end)
 end
 
 """ TMP TKE profile for testing """
@@ -50,7 +50,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]
     not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
-    return x -> FT(not_type_stable_spline(x))
+    return ZProfile(x -> FT(not_type_stable_spline(x)))
 end
 
 """ [Brown2002](@cite) """
@@ -68,13 +68,14 @@ function ARM_SGP_dTdt(::Type{FT}) where {FT}
     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)
-    return (t, z) -> if z <= 1000.0
-        dTdt_A(t)+dTdt_R(t)
+    prof = (t, z) -> if z <= 1000.0
+        FT(dTdt_A(t)+dTdt_R(t))
     elseif z > 1000.0 && z <= 2000.0
-        (dTdt_A(t)+dTdt_R(t)) * (1 - (z - 1000) / 1000)
+        FT((dTdt_A(t)+dTdt_R(t)) * (1 - (z - 1000) / 1000))
     else
         FT(0)
     end
+    return TimeZProfile(prof)
 end
 
 """ [Brown2002](@cite) """
@@ -83,13 +84,14 @@ function ARM_SGP_dqtdt(::Type{FT}) where {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)
-    return (Π, t, z) -> if z <= 1000.0
-        dqtdt(t) * Π
+    prof = (Π, t, z) -> if z <= 1000.0
+        FT(dqtdt(t) * Π)
     elseif z > 1000.0 && z <= 2000.0
-        dqtdt(t) * Π * (1 - (z - 1000) / 1000)
+        FT(dqtdt(t) * Π * (1 - (z - 1000) / 1000))
     else
         FT(0)
     end
+    return ΠTimeZProfile(prof)
 end
 
 """ [Brown2002](@cite) """
@@ -97,7 +99,7 @@ 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)
-    return profile
+    return TimeProfile(profile)
 end
 
 """ [Brown2002](@cite) """
@@ -105,5 +107,5 @@ 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)
-    return profile
+    return TimeProfile(profile)
 end
diff --git a/src/profiles/Bomex.jl b/src/profiles/Bomex.jl
index 9764bb5..1a9e06b 100644
--- a/src/profiles/Bomex.jl
+++ b/src/profiles/Bomex.jl
@@ -1,5 +1,5 @@
 """ [Siebesma2003](@cite) """
-Bomex_θ_liq_ice(::Type{FT}) where {FT} = z -> if z <= 520.0
+Bomex_θ_liq_ice(::Type{FT}) where {FT} = ZProfile(z -> if z <= 520.0
         FT(298.7)
     elseif z > 520.0 && z <= 1480.0
         FT(298.7) + (z - 520) * (FT(302.4) - FT(298.7)) / (1480 - 520)
@@ -9,29 +9,29 @@ Bomex_θ_liq_ice(::Type{FT}) where {FT} = z -> if z <= 520.0
         FT(308.2) + (z - 2000) * (FT(311.85) - FT(308.2)) / (3000 - 2000)
     else
         FT(0)
-    end
+    end)
 """ [Siebesma2003](@cite) """
-Bomex_q_tot(::Type{FT}) where {FT} = z -> if z <= 520
+Bomex_q_tot(::Type{FT}) where {FT} = ZProfile(z -> if z <= 520
         (FT(17.0) + z * (FT(16.3) - 17) / 520) / 1000
     elseif z > 520.0 && z <= 1480.0
         (FT(16.3) + (z - 520) * (FT(10.7) - FT(16.3)) / (1480 - 520)) / 1000
     elseif z > 1480.0 && z <= 2000.0
         (FT(10.7) + (z - 1480) * (FT(4.2) - FT(10.7)) / (2000 - 1480)) / 1000
-    elseif z > 2000.0
+    else
         (FT(4.2) + (z - 2000) * (3 - FT(4.2)) / (3000 - 2000)) / 1000
-    end
+    end)
 """ [Siebesma2003](@cite) """
-Bomex_u(::Type{FT}) where {FT} = z -> if z <= 700.0
+Bomex_u(::Type{FT}) where {FT} = ZProfile(z -> if z <= 700.0
         FT(-8.75)
-    elseif z > 700.0
+    else
         FT(-8.75) + (z - 700) * (FT(-4.61) - FT(-8.75)) / (3000 - 700)
-    end
+    end)
 """ [Siebesma2003](@cite) """
-Bomex_tke(::Type{FT}) where {FT} = z -> if (z <= 2500.0)
+Bomex_tke(::Type{FT}) where {FT} = ZProfile(z -> if (z <= 2500.0)
         FT(1) - z / 3000
     else
         FT(0)
-    end
+    end)
 """ TMP TKE profile for testing """
 function Bomex_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, 425.0, 475.0, 525.0,
@@ -46,38 +46,38 @@ function Bomex_tke_prescribed(::Type{FT}) where {FT}
                 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)
-    return x -> FT(not_type_stable_spline(x))
+    return ZProfile(x -> FT(not_type_stable_spline(x)))
 end
 
 # Geostrophic velocity profiles
 """ :( """
-Bomex_geostrophic_u(::Type{FT}) where {FT} = z -> -10 + FT(1.8e-3) * z
+Bomex_geostrophic_u(::Type{FT}) where {FT} = ZProfile(z -> -10 + FT(1.8e-3) * z)
 """ :( """
-Bomex_geostrophic_v(::Type{FT}) where {FT} = z -> FT(0)
+Bomex_geostrophic_v(::Type{FT}) where {FT} = ZProfile(z -> FT(0))
 
 """ dTdt(Π, z) """
-Bomex_dTdt(::Type{FT}) where {FT} = (Π, z) -> if z <= 1500.0
+Bomex_dTdt(::Type{FT}) where {FT} = ΠZProfile((Π, z) -> if z <= 1500.0
         (-2 / (3600 * 24)) * Π
     else
         (-2 / (3600 * 24) + (z - 1500) * (0 - -2 / (3600 * 24)) / (3000 - 1500)) * Π
-    end
+    end)
 
 # Large-scale drying
 """ :( """
-Bomex_dqtdt(::Type{FT}) where {FT} = z -> if z <= 300.0
+Bomex_dqtdt(::Type{FT}) where {FT} = ZProfile(z -> if z <= 300.0
         FT(-1.2e-8)   #kg/(kg * s)
     elseif z > 300.0 && z <= 500.0
         FT(-1.2e-8) + (z - 300) * (0 - FT(-1.2e-8)) / (500 - 300) #kg/(kg * s)
     else
         FT(0)
-    end
+    end)
 
 # Large scale subsidence
 """ :( """
-Bomex_subsidence(::Type{FT}) where {FT} = z -> if z <= 1500.0
+Bomex_subsidence(::Type{FT}) where {FT} = ZProfile(z -> if z <= 1500.0
         FT(0) + z * (FT(-0.65) / 100 - 0) / (1500 - 0)
     elseif z > 1500.0 && z <= 2100.0
         FT(-0.65) / 100 + (z - 1500) * (0 - FT(-0.65) / 100) / (2100 - 1500)
     else
         FT(0)
-    end
+    end)
diff --git a/src/profiles/DryBubble.jl b/src/profiles/DryBubble.jl
index 54219fb..445db6b 100644
--- a/src/profiles/DryBubble.jl
+++ b/src/profiles/DryBubble.jl
@@ -62,7 +62,7 @@ function DryBubble_θ_liq_ice(::Type{FT}) where {FT}
         299.9837, 299.9837
     ]
     profile = Dierckx.Spline1D(z_in, θ_liq_ice_in; k = 1)
-    return profile
+    return ZProfile(profile)
 end
 
 """ :( """
@@ -98,7 +98,7 @@ function DryBubble_updrafts_θ_liq_ice(::Type{FT}) where {FT}
         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)
-    return profile
+    return ZProfile(profile)
 end
 
 """ :( """
@@ -115,7 +115,7 @@ function DryBubble_updrafts_area(::Type{FT}) where {FT}
         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)
-    return profile
+    return ZProfile(profile)
 end
 
 """ :( """
@@ -133,7 +133,7 @@ function DryBubble_updrafts_w(::Type{FT}) where {FT}
         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)
-    return profile
+    return ZProfile(profile)
 end
 
 """ :( """
@@ -154,5 +154,5 @@ function DryBubble_updrafts_T(::Type{FT}) where {FT}
         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)
-    return profile
+    return ZProfile(profile)
 end
diff --git a/src/profiles/Dycoms_RF01.jl b/src/profiles/Dycoms_RF01.jl
index 6024d76..bfa2536 100644
--- a/src/profiles/Dycoms_RF01.jl
+++ b/src/profiles/Dycoms_RF01.jl
@@ -1,29 +1,29 @@
 """ [Stevens2005](@cite) """
-Dycoms_RF01_q_tot(::Type{FT}) where {FT} = z -> if z <= 840.0
+Dycoms_RF01_q_tot(::Type{FT}) where {FT} = ZProfile(z -> if z <= 840.0
         FT(9.0) / FT(1000.0)
     else
         FT(1.5) / FT(1000.0)
-    end
+    end)
 
 """ [Stevens2005](@cite) """
-Dycoms_RF01_θ_liq_ice(::Type{FT}) where {FT} = z -> if z <= 840.0
+Dycoms_RF01_θ_liq_ice(::Type{FT}) where {FT} = ZProfile(z -> if z <= 840.0
         FT(289)
     else
         FT(297.5) + (z - FT(840))^FT(1.0 / 3.0)
-    end
+    end)
 
 """ [Stevens2005](@cite) """
-Dycoms_RF01_u0(::Type{FT}) where {FT} = z -> FT(7)
+Dycoms_RF01_u0(::Type{FT}) where {FT} = ZProfile(z -> FT(7))
 
 """ [Stevens2005](@cite) """
-Dycoms_RF01_v0(::Type{FT}) where {FT} = z -> FT(-5.5)
+Dycoms_RF01_v0(::Type{FT}) where {FT} = ZProfile(z -> FT(-5.5))
 
 """ [Stevens2005](@cite) """
-Dycoms_RF01_tke(::Type{FT}) where {FT} = z -> if z <= 800.0
+Dycoms_RF01_tke(::Type{FT}) where {FT} = ZProfile(z -> if z <= 800.0
         FT(1) - z / FT(1000)
     else
         FT(0)
-    end
+    end)
 """ TMP TKE profile for testing """
 function Dycoms_RF01_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, 425.0, 475.0, 525.0, 575.0, 625.0, 675.0, 725.0,
@@ -31,6 +31,6 @@ function Dycoms_RF01_tke_prescribed(::Type{FT}) where {FT}
     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)
-    return x -> FT(not_type_stable_spline(x))
+    return ZProfile(x -> FT(not_type_stable_spline(x)))
 end
 
diff --git a/src/profiles/Dycoms_RF02.jl b/src/profiles/Dycoms_RF02.jl
index c449e74..3490fe9 100644
--- a/src/profiles/Dycoms_RF02.jl
+++ b/src/profiles/Dycoms_RF02.jl
@@ -1,33 +1,33 @@
 """ [Ackerman2009](@cite) """
 Dycoms_RF02_q_tot(::Type{FT}) where {FT} =
-    z -> if z <= 795.0
+    ZProfile(z -> if z <= 795.0
         FT(9.45) / FT(1000.0)
     else
         (FT(5) - FT(3) * (FT(1) - exp(-(z - FT(795)) / FT(500)))) / FT(1000)
-    end
+    end)
 
 """ [Ackerman2009](@cite) """
 Dycoms_RF02_θ_liq_ice(::Type{FT}) where {FT} =
-    z -> if z <= 795.0
+    ZProfile(z -> if z <= 795.0
         FT(288.3)
     else
         FT(295.0) + (z - FT(795))^FT(1.0 / 3.0)
-    end
+    end)
 
 """ [Ackerman2009](@cite) """
 function Dycoms_RF02_u(::Type{FT}) where {FT}
-    return z -> FT(3) + FT(4.3) * z / FT(1000) - FT(5)
+    return ZProfile(z -> FT(3) + FT(4.3) * z / FT(1000) - FT(5))
 end
 
 """ [Ackerman2009](@cite) """
 function Dycoms_RF02_v(::Type{FT}) where {FT}
-    return z -> FT(-9) + FT(5.6) * z / FT(1000) - FT(-5.5)
+    return ZProfile(z -> FT(-9) + FT(5.6) * z / FT(1000) - FT(-5.5))
 end
 
 """ [Ackerman2009](@cite) """
 Dycoms_RF02_tke(::Type{FT}) where {FT} =
-    z -> if z <= 795.0
+    ZProfile(z -> if z <= 795.0
         FT(1) - z / FT(1000)
     else
         FT(0)
-    end
+    end)
diff --git a/src/profiles/GABLS.jl b/src/profiles/GABLS.jl
index 81b54d5..4829936 100644
--- a/src/profiles/GABLS.jl
+++ b/src/profiles/GABLS.jl
@@ -1,36 +1,36 @@
 """ :( """
-GABLS_u(::Type{FT}) where {FT} = z -> FT(8.0)
+GABLS_u(::Type{FT}) where {FT} = ZProfile(z -> FT(8.0))
 """ :( """
-GABLS_v(::Type{FT}) where {FT} = z -> FT(0.0)
+GABLS_v(::Type{FT}) where {FT} = ZProfile(z -> FT(0.0))
 
 """ :( """
 GABLS_θ_liq_ice(::Type{FT}) where {FT} =
-    z -> if z <= 100.0
+    ZProfile(z -> if z <= 100.0
         FT(265)
     else
         FT(265) + (z - 100) * FT(0.01)
-    end
+    end)
 
 """ :( """
-GABLS_q_tot(::Type{FT}) where {FT} = z -> FT(0)
+GABLS_q_tot(::Type{FT}) where {FT} = ZProfile(z -> FT(0))
 
 """ :( """
 GABLS_tke(::Type{FT}) where {FT} =
-    z -> if z <= 250.0
+    ZProfile(z -> if z <= 250.0
         FT(0.4) * (1 - z / 250) * (1 - z / 250) * (1 - z / 250)
     else
         FT(0)
-    end
+    end)
 
 """ TMP TKE profile for testing """
 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)
-    return x -> FT(not_type_stable_spline(x))
+    return ZProfile(x -> FT(not_type_stable_spline(x)))
 end
 
 """ :( """
-GABLS_geostrophic_ug(::Type{FT}) where {FT} = z -> FT(8)
+GABLS_geostrophic_ug(::Type{FT}) where {FT} = ZProfile(z -> FT(8))
 """ :( """
-GABLS_geostrophic_vg(::Type{FT}) where {FT} = z -> FT(0)
+GABLS_geostrophic_vg(::Type{FT}) where {FT} = ZProfile(z -> FT(0))
diff --git a/src/profiles/GATE_III.jl b/src/profiles/GATE_III.jl
index 59a64a9..853fd4b 100644
--- a/src/profiles/GATE_III.jl
+++ b/src/profiles/GATE_III.jl
@@ -16,7 +16,7 @@ function GATE_III_q_tot(::Type{FT}) where {FT}
                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)
-    return profile
+    return ZProfile(profile)
 end
 
 
@@ -27,7 +27,7 @@ function GATE_III_u(::Type{FT}) where {FT}
                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)
-    return profile
+    return ZProfile(profile)
 end
 
 """ [Khairoutdinov2009](@cite) """
@@ -37,16 +37,16 @@ function GATE_III_T(::Type{FT}) where {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)
-    return profile
+    return ZProfile(profile)
 end
 
 """ [Khairoutdinov2009](@cite) """
 function GATE_III_tke(::Type{FT}) where {FT}
-    return z -> if z <= 2500.0
-            FT(1) - z / 3000.0
+    return ZProfile(z -> if z <= 2500.0
+            FT(1 - z / 3000)
         else
             FT(0)
-        end
+        end)
 end
 
 #=
@@ -78,7 +78,7 @@ function GATE_III_dqtdt(::Type{FT}) where {FT}
 
     Qtend_in = r_tend_in ./ (1 .+ r_tend_in) # convert mixing ratio to specific humidity
     profile = Dierckx.Spline1D(z_in, Qtend_in; k = 1)
-    return profile
+    return ZProfile(profile)
 end
 
 """ [Khairoutdinov2009](@cite) """
@@ -99,5 +99,5 @@ function GATE_III_dTdt(::Type{FT}) where {FT}
 
     profile_T = Dierckx.Spline1D(z_in, Ttend_in; k = 1)
     profile_R = Dierckx.Spline1D(z_in, RAD_in; k = 1)
-    return z -> profile_T(z) + profile_R(z)
+    return ZProfile(z -> profile_T(z) + profile_R(z))
 end
diff --git a/src/profiles/ISDAC.jl b/src/profiles/ISDAC.jl
index 47c0290..f0e5394 100644
--- a/src/profiles/ISDAC.jl
+++ b/src/profiles/ISDAC.jl
@@ -4,7 +4,7 @@
 
 """ [ovchinnikov_intercomparison_2014](@cite) """
 function ISDAC_q_tot(::Type{FT}) where {FT}
-    z -> FT(if z < 400
+    ZProfile(z -> FT(if z < 400
         1.5 - 0.00075 * (z - 400)
     elseif 400 ≤ z < 825
         1.5
@@ -12,12 +12,12 @@ function ISDAC_q_tot(::Type{FT}) where {FT}
         1.2
     else  #2045 ≤ z
         0.5 - 0.000075 * (z - 2045)
-    end / 1000)  # kg/kg
+    end / 1000))  # kg/kg
 end
 
 """ [ovchinnikov_intercomparison_2014](@cite) """
 function ISDAC_θ_liq_ice(::Type{FT}) where {FT}
-    z -> FT(if z < 400
+    ZProfile(z -> FT(if z < 400
         265 + 0.004 * (z - 400)
     elseif 400 ≤ z < 825
         265
@@ -25,43 +25,43 @@ function ISDAC_θ_liq_ice(::Type{FT}) where {FT}
         266 + (z - 825)^0.3
     else  #2045 ≤ z
         271 + (z - 2000)^0.33
-    end)  # K
+    end))  # K
 end
 
 """ [ovchinnikov_intercomparison_2014](@cite) """
-ISDAC_u(::Type{FT}) where {FT} = z -> FT(-7)  # m/s
+ISDAC_u(::Type{FT}) where {FT} = ZProfile(z -> FT(-7))  # m/s
 
 """ [ovchinnikov_intercomparison_2014](@cite) """
-ISDAC_v(::Type{FT}) where {FT} = z -> FT(-2 + 0.003z)  # m/s
+ISDAC_v(::Type{FT}) where {FT} = ZProfile(z -> FT(-2 + 0.003z))  # m/s
 
 """ [ovchinnikov_intercomparison_2014](@cite) """
 ISDAC_subsidence(::Type{FT}) where {FT} =
-    z -> FT(z < 825 ? -5e-6z : -0.4125e-2)  # m/s
+    ZProfile(z -> FT(z < 825 ? -5e-6z : -0.4125e-2))  # m/s
 
 """ [ovchinnikov_intercomparison_2014](@cite) """
-ISDAC_tke(::Type{FT}) where {FT} = z -> FT(0.1)  # m²/s²
+ISDAC_tke(::Type{FT}) where {FT} = ZProfile(z -> FT(0.1))  # m²/s²
 
 """ [ovchinnikov_intercomparison_2014](@cite) """
 function ISDAC_inv_τ_scalar(::Type{FT}) where {FT}
     z₁ = 1200  # m
     z₂ = 1500  # m
     hr = 3600  # s
-    z -> FT(if z < z₁
+    ZProfile(z -> FT(if z < z₁
         0
     elseif z₁ ≤ z ≤ z₂
         1 / hr * (1 - cos(π * (z - z₁) / (z₂ - z₁))) / 2
     else  #z > z₂
         1 / hr
-    end)  # s⁻¹
+    end))  # s⁻¹
 end
 
 """ [ovchinnikov_intercomparison_2014](@cite) """
 function ISDAC_inv_τ_wind(::Type{FT}) where {FT}
     zᵤᵥ = 825  # m
     hr = 3600  # s
-    z -> FT(if z ≤ zᵤᵥ
+    ZProfile(z -> FT(if z ≤ zᵤᵥ
         1 / 2hr * (1 - cos(π * z / zᵤᵥ)) / 2
     else  #z > zᵤᵥ
         1 / 2hr
-    end)  # s⁻¹
+    end))  # s⁻¹
 end
diff --git a/src/profiles/LifeCycleTan2018.jl b/src/profiles/LifeCycleTan2018.jl
index d1e8381..71305ae 100644
--- a/src/profiles/LifeCycleTan2018.jl
+++ b/src/profiles/LifeCycleTan2018.jl
@@ -1,6 +1,6 @@
 """ [Tan2018](@cite) """
 LifeCycleTan2018_θ_liq_ice(::Type{FT}) where {FT} =
-    z -> if z <= 520.0
+    ZProfile(z -> if z <= 520.0
         FT(298.7)
     elseif z > 520.0 && z <= 1480.0
         FT(298.7) + (z - 520) * (FT(302.4) - FT(298.7)) / (1480 - 520)
@@ -10,11 +10,11 @@ LifeCycleTan2018_θ_liq_ice(::Type{FT}) where {FT} =
         FT(308.2) + (z - 2000) * (FT(311.85) - FT(308.2)) / (3000 - 2000)
     else
         FT(0)
-    end
+    end)
 
 """ [Tan2018](@cite) """
 LifeCycleTan2018_q_tot(::Type{FT}) where {FT} =
-    z -> if z <= 520
+    ZProfile(z -> if z <= 520
         (FT(17) + z * (FT(16.3) - FT(17.0)) / 520) / 1000
     elseif z > 520.0 && z <= 1480.0
         (FT(16.3) + (z - 520) * (FT(10.7) - FT(16.3)) / (1480 - 520)) / 1000
@@ -24,22 +24,22 @@ LifeCycleTan2018_q_tot(::Type{FT}) where {FT} =
         (FT(4.2) + (z - 2000) * (3 - FT(4.2)) / (3000 - 2000)) / 1000
     else
         FT(0)
-    end
+    end)
 
 """ [Tan2018](@cite) """
 LifeCycleTan2018_u(::Type{FT}) where {FT} =
-    z -> if z <= 700.0
+    ZProfile(z -> if z <= 700.0
         FT(-8.75)
     else
         FT(-8.75) + (z - 700) * (FT(-4.61) - FT(-8.75)) / (3000 - 700)
-    end
+    end)
 """ [Tan2018](@cite) """
 LifeCycleTan2018_tke(::Type{FT}) where {FT} =
-    z -> if z <= 2500.0
+    ZProfile(z -> if z <= 2500.0
         FT(1) - z / 3000
     else
         FT(0)
-    end
+    end)
 
 """ TMP TKE profile for testing """
 function LifeCycleTan2018_tke_prescribed(::Type{FT}) where {FT}
@@ -58,45 +58,45 @@ 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]
     not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
-    return x -> FT(not_type_stable_spline(x))
+    return ZProfile(x -> FT(not_type_stable_spline(x)))
 end
 
 # Large-scale cooling
 """ [Tan2018](@cite) """
 LifeCycleTan2018_dTdt(::Type{FT}) where {FT} =
-    (Π, z) -> if z <= 1500.0
+    ΠZProfile((Π, z) -> if z <= 1500.0
         FT(-2 / (3600 * 24)) * Π
     else
         FT(-2 / (3600 * 24) + (z - 1500) * (0 - -2 / (3600 * 24)) / (3000 - 1500)) * Π
-    end
+    end)
 
 # geostrophic velocity profiles
 """ [Tan2018](@cite) """
 LifeCycleTan2018_geostrophic_u(::Type{FT}) where {FT} =
-    z -> -10 + FT(1.8e-3) * z
+    ZProfile(z -> -10 + FT(1.8e-3) * z)
 
 """ [Tan2018](@cite) """
 LifeCycleTan2018_geostrophic_v(::Type{FT}) where {FT} =
-    z -> FT(0)
+    ZProfile(z -> FT(0))
 
 # Large-scale drying
 """ [Tan2018](@cite) """
 LifeCycleTan2018_dqtdt(::Type{FT}) where {FT} =
-    z -> if z <= 300.0
+    ZProfile(z -> if z <= 300.0
         FT(-1.2e-8)   #kg/(kg * s)
     elseif z > 300.0 && z <= 500.0
         FT(-1.2e-8) + (z - 300) * (0 - FT(-1.2e-8)) / (500 - 300) #kg/(kg * s)
     else
         FT(0)
-    end
+    end)
 
 #Large scale subsidence
 """ [Tan2018](@cite) """
 LifeCycleTan2018_subsidence(::Type{FT}) where {FT} =
-    z -> if z <= 1500.0
+    ZProfile(z -> if z <= 1500.0
         FT(0) + z * (FT(-0.65) / 100 - 0) / (1500 - 0)
     elseif z > 1500.0 && z <= 2100.0
         FT(-0.65) / 100 + (z - 1500) * (0 - FT(-0.65) / 100) / (2100 - 1500)
     else
         FT(0)
-    end
+    end)
diff --git a/src/profiles/Nieuwstadt.jl b/src/profiles/Nieuwstadt.jl
index 8d2fa53..7a74c62 100644
--- a/src/profiles/Nieuwstadt.jl
+++ b/src/profiles/Nieuwstadt.jl
@@ -1,19 +1,19 @@
 """ [Nieuwstadt1993](@cite) """
-Nieuwstadt_θ_liq_ice(::Type{FT}) where {FT} = z -> if z <= 1350.0
+Nieuwstadt_θ_liq_ice(::Type{FT}) where {FT} = ZProfile(z -> if z <= 1350.0
         FT(300)
     else
         FT(300) + 3 * (z - 1350) / 1000
-    end
+    end)
 
 """ [Nieuwstadt1993](@cite) """
-Nieuwstadt_u(::Type{FT}) where {FT} = z -> FT(0.01)
+Nieuwstadt_u(::Type{FT}) where {FT} = ZProfile(z -> FT(0.01))
 
 """ [Nieuwstadt1993](@cite) """
-Nieuwstadt_tke(::Type{FT}) where {FT} = z -> if (z <= 1600.0)
+Nieuwstadt_tke(::Type{FT}) where {FT} = ZProfile(z -> if (z <= 1600.0)
         FT(0.1) * FT(1.46) * FT(1.46) * (1 - z / 1600)
     else
         FT(0)
-    end
+    end)
 
 """ TMP TKE profile for testing """
 function Nieuwstadt_tke_prescribed(::Type{FT}) where {FT}
@@ -52,5 +52,5 @@ 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]
     not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
-    return x -> FT(not_type_stable_spline(x))
+    return ZProfile(x -> FT(not_type_stable_spline(x)))
 end
diff --git a/src/profiles/Rico.jl b/src/profiles/Rico.jl
index 216338e..4d472ae 100644
--- a/src/profiles/Rico.jl
+++ b/src/profiles/Rico.jl
@@ -1,62 +1,62 @@
 """ [vanZanten2011](@cite) """
 function Rico_u(::Type{FT}) where {FT}
-    return z -> FT(-9.9) + FT(2.0e-3) * z
+    return ZProfile(z -> FT(-9.9) + FT(2.0e-3) * z)
 end
 
 """ [vanZanten2011](@cite) """
 function Rico_v(::Type{FT}) where {FT}
-    return z -> FT(-3.8)
+    return ZProfile(z -> FT(-3.8))
 end
 
 """ [vanZanten2011](@cite) """
 Rico_θ_liq_ice(::Type{FT}) where {FT} =
-    z -> if z <= 740.0
+    ZProfile(z -> if z <= 740.0
         FT(297.9)
     else
         FT(297.9) + (317 - FT(297.9)) / (4000 - 740) * (z - 740)
-    end
+    end)
 
 """ [vanZanten2011](@cite) """
 Rico_q_tot(::Type{FT}) where {FT} =
-    z -> if z <= 740.0
+    ZProfile(z -> if z <= 740.0
         (16 + (FT(13.8) - 16) / 740 * z) / 1000
     elseif z > 740.0 && z <= 3260.0
         (FT(13.8) + (FT(2.4) - FT(13.8)) / (3260 - 740) * (z - 740)) / 1000
     else
         (FT(2.4) + (FT(1.8) - FT(2.4)) / (4000 - 3260) * (z - 3260)) / 1000
-    end
+    end)
 
 """ [vanZanten2011](@cite) """
 function Rico_geostrophic_ug(::Type{FT}) where {FT}
-    return z -> FT(-9.9) + FT(2.0e-3) * z
+    return ZProfile(z -> FT(-9.9) + FT(2.0e-3) * z)
 end
 
 """ [vanZanten2011](@cite) """
 function Rico_geostrophic_vg(::Type{FT}) where {FT}
-    return z -> FT(-3.8)
+    return ZProfile(z -> FT(-3.8))
 end
 
 """ [vanZanten2011](@cite) """
 function Rico_dTdt(::Type{FT}) where {FT}
-    return (Π, z) -> (FT(-2.5) / (3600 * 24)) * Π
+    return ΠZProfile((Π, z) -> (FT(-2.5) / (3600 * 24)) * Π)
 end
 
 """ [vanZanten2011](@cite) """
 function Rico_dqtdt(::Type{FT}) where {FT}
-    return z -> if z <= 2980.0
+    return ZProfile(z -> if z <= 2980.0
         (-1 + FT(1.3456) / FT(2980.0) * z) / FT(86400.0) / 1000   #kg/(kg * s)
     else
         FT(0.3456) / 86400 / 1000
-    end
+    end)
 end
 
 """ [vanZanten2011](@cite) """
 function Rico_subsidence(::Type{FT}) where {FT}
-    return z -> if z <= 2260.0
+    return ZProfile(z -> if z <= 2260.0
             -(FT(0.005) / FT(2260.0)) * z
         else
             FT(-0.005)
-        end
+        end)
 end
 """ TMP TKE profile for testing """
 function Rico_tke_prescribed(::Type{FT}) where {FT}
@@ -79,5 +79,5 @@ 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]
     not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
-    return x -> FT(not_type_stable_spline(x))
+    return ZProfile(x -> FT(not_type_stable_spline(x)))
 end
diff --git a/src/profiles/SP.jl b/src/profiles/SP.jl
index c53e843..99bbade 100644
--- a/src/profiles/SP.jl
+++ b/src/profiles/SP.jl
@@ -1,24 +1,24 @@
 """ :( """
-SP_u(::Type{FT}) where {FT} = z -> FT(1)
+SP_u(::Type{FT}) where {FT} = ZProfile(z -> FT(1))
 
 """ :( """
-SP_v(::Type{FT}) where {FT} = z -> FT(0)
+SP_v(::Type{FT}) where {FT} = ZProfile(z -> FT(0))
 
 """ :( """
 SP_θ_liq_ice(::Type{FT}) where {FT} =
-    z -> if z <= 974.0
+    ZProfile(z -> if z <= 974.0
         FT(300)
     elseif z < 1074.0
         FT(300) + (z - 974) * FT(0.08)
     else
         FT(308) + (z - 1074) * FT(0.003)
-    end
+    end)
 
 """ :( """
-SP_q_tot(::Type{FT}) where {FT} = z -> FT(0)
+SP_q_tot(::Type{FT}) where {FT} = ZProfile(z -> FT(0))
 
 """ :( """
-SP_geostrophic_u(::Type{FT}) where {FT} = z -> FT(1)
+SP_geostrophic_u(::Type{FT}) where {FT} = ZProfile(z -> FT(1))
 
 """ :( """
-SP_geostrophic_v(::Type{FT}) where {FT} = z -> FT(0)
+SP_geostrophic_v(::Type{FT}) where {FT} = ZProfile(z -> FT(0))
diff --git a/src/profiles/Soares.jl b/src/profiles/Soares.jl
index afbfb07..276a986 100644
--- a/src/profiles/Soares.jl
+++ b/src/profiles/Soares.jl
@@ -1,24 +1,24 @@
 """ [Soares2004](@cite) """
-Soares_q_tot(::Type{FT}) where {FT} = z -> if z <= 1350.0
+Soares_q_tot(::Type{FT}) where {FT} = ZProfile(z -> if z <= 1350.0
         FT(5.0e-3) - FT(3.7e-4) * z / 1000
     else
         FT(5.0e-3) - FT(3.7e-4) * FT(1.35) - FT(9.4e-4) * (z - 1350) / 1000
-    end
+    end)
 """ [Soares2004](@cite) """
-Soares_θ_liq_ice(::Type{FT}) where {FT} = z -> if z <= 1350.0
+Soares_θ_liq_ice(::Type{FT}) where {FT} = ZProfile(z -> if z <= 1350.0
         FT(300.0)
     else
         FT(300) + 2 * (z - 1350) / 1000
-    end
+    end)
 """ [Soares2004](@cite) """
-Soares_u(::Type{FT}) where {FT} = z -> FT(0.01)
+Soares_u(::Type{FT}) where {FT} = ZProfile(z -> FT(0.01))
 
 """ [Soares2004](@cite) """
-Soares_tke(::Type{FT}) where {FT} = z -> if z <= 1600.0
+Soares_tke(::Type{FT}) where {FT} = ZProfile(z -> if z <= 1600.0
         FT(0.1) * FT(1.46) * FT(1.46) * (1 - z / 1600)
     else
         FT(0)
-    end
+    end)
 """ TMP TKE profile for testing """
 function Soares_tke_prescribed(::Type{FT}) where {FT}
     z_in = FT[12.5, 37.5, 62.5, 87.5, 112.5, 137.5, 162.5, 187.5, 212.5, 237.5, 262.5, 287.5, 312.5, 337.5,
@@ -47,5 +47,5 @@ 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]
     not_type_stable_spline = Dierckx.Spline1D(z_in, tke_in; k = 1)
-    return x -> FT(not_type_stable_spline(x))
+    return ZProfile(x -> FT(not_type_stable_spline(x)))
 end
diff --git a/src/profiles/TRMM_LBA.jl b/src/profiles/TRMM_LBA.jl
index c3fa9c8..a4872de 100644
--- a/src/profiles/TRMM_LBA.jl
+++ b/src/profiles/TRMM_LBA.jl
@@ -21,7 +21,8 @@ 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)
+    prof = Dierckx.Spline1D(z_in, p_in; k = 1)
+    return ZProfile(prof)
 end
 """ [Grabowski2006](@cite) """
 function TRMM_LBA_T(::Type{FT}) where {FT}
@@ -33,7 +34,8 @@ 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)
+    prof = Dierckx.Spline1D(z_in, T_in; k = 1)
+    return ZProfile(prof)
 end
 """ [Grabowski2006](@cite) """
 function TRMM_LBA_RH(::Type{FT}) where {FT}
@@ -45,7 +47,8 @@ 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)
+    prof = Dierckx.Spline1D(z_in, RH_in; k = 1)
+    return ZProfile(prof)
 end
 """ [Grabowski2006](@cite) """
 function TRMM_LBA_u(::Type{FT}) where {FT}
@@ -57,7 +60,8 @@ 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)
+    prof = Dierckx.Spline1D(z_in, u_in; k = 1)
+    return ZProfile(prof)
 end
 """ [Grabowski2006](@cite) """
 function TRMM_LBA_v(::Type{FT}) where {FT}
@@ -69,16 +73,17 @@ 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)
+    prof = Dierckx.Spline1D(z_in, v_in; k = 1)
+    return ZProfile(prof)
 end
 
 """ [Grabowski2006](@cite) """
 function TRMM_LBA_tke(::Type{FT}) where {FT}
-    z -> if z <= 2500.0
+    ZProfile(z -> if z <= 2500.0
         FT(1) - z / 3000
     else
         FT(0)
-    end
+    end)
 end
 """ TMP TKE profile for testing """
 function TRMM_LBA_tke_prescribed(::Type{FT}) where {FT}
@@ -106,7 +111,7 @@ function TRMM_LBA_tke_prescribed(::Type{FT}) where {FT}
        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)
-    return x -> FT(not_type_stable_spline(x))
+    return ZProfile(x -> FT(not_type_stable_spline(x)))
 end
 
 
@@ -230,5 +235,5 @@ function TRMM_LBA_radiation(::Type{FT}) where {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)
-    return profile
+    return TimeZProfile(profile)
 end
diff --git a/test/profiles.jl b/test/profiles.jl
new file mode 100644
index 0000000..218498c
--- /dev/null
+++ b/test/profiles.jl
@@ -0,0 +1,103 @@
+import AtmosphericProfilesLibrary as APL
+profiles(FT) = [
+(:ARM_SGP_dTdt, APL.ARM_SGP_dTdt(FT)),
+(:ARM_SGP_dqtdt, APL.ARM_SGP_dqtdt(FT)),
+(:ARM_SGP_lhf, APL.ARM_SGP_lhf(FT)),
+(:ARM_SGP_q_tot, APL.ARM_SGP_q_tot(FT)),
+(:ARM_SGP_shf, APL.ARM_SGP_shf(FT)),
+(:ARM_SGP_tke, APL.ARM_SGP_tke(FT)),
+(:ARM_SGP_tke_prescribed, APL.ARM_SGP_tke_prescribed(FT)),
+(:ARM_SGP_u, APL.ARM_SGP_u(FT)),
+(:ARM_SGP_θ_liq_ice, APL.ARM_SGP_θ_liq_ice(FT)),
+(:Bomex_dTdt, APL.Bomex_dTdt(FT)),
+(:Bomex_dqtdt, APL.Bomex_dqtdt(FT)),
+(:Bomex_geostrophic_u, APL.Bomex_geostrophic_u(FT)),
+(:Bomex_geostrophic_v, APL.Bomex_geostrophic_v(FT)),
+(:Bomex_q_tot, APL.Bomex_q_tot(FT)),
+(:Bomex_subsidence, APL.Bomex_subsidence(FT)),
+(:Bomex_tke, APL.Bomex_tke(FT)),
+(:Bomex_tke_prescribed, APL.Bomex_tke_prescribed(FT)),
+(:Bomex_u, APL.Bomex_u(FT)),
+(:Bomex_θ_liq_ice, APL.Bomex_θ_liq_ice(FT)),
+(:DryBubble_updrafts_T, APL.DryBubble_updrafts_T(FT)),
+(:DryBubble_updrafts_area, APL.DryBubble_updrafts_area(FT)),
+(:DryBubble_updrafts_w, APL.DryBubble_updrafts_w(FT)),
+(:DryBubble_updrafts_θ_liq_ice, APL.DryBubble_updrafts_θ_liq_ice(FT)),
+(:DryBubble_θ_liq_ice, APL.DryBubble_θ_liq_ice(FT)),
+(:Dycoms_RF01_q_tot, APL.Dycoms_RF01_q_tot(FT)),
+(:Dycoms_RF01_tke, APL.Dycoms_RF01_tke(FT)),
+(:Dycoms_RF01_tke_prescribed, APL.Dycoms_RF01_tke_prescribed(FT)),
+(:Dycoms_RF01_u0, APL.Dycoms_RF01_u0(FT)),
+(:Dycoms_RF01_v0, APL.Dycoms_RF01_v0(FT)),
+(:Dycoms_RF01_θ_liq_ice, APL.Dycoms_RF01_θ_liq_ice(FT)),
+(:Dycoms_RF02_q_tot, APL.Dycoms_RF02_q_tot(FT)),
+(:Dycoms_RF02_tke, APL.Dycoms_RF02_tke(FT)),
+(:Dycoms_RF02_u, APL.Dycoms_RF02_u(FT)),
+(:Dycoms_RF02_v, APL.Dycoms_RF02_v(FT)),
+(:Dycoms_RF02_θ_liq_ice, APL.Dycoms_RF02_θ_liq_ice(FT)),
+(:GABLS_geostrophic_ug, APL.GABLS_geostrophic_ug(FT)),
+(:GABLS_geostrophic_vg, APL.GABLS_geostrophic_vg(FT)),
+(:GABLS_q_tot, APL.GABLS_q_tot(FT)),
+(:GABLS_tke, APL.GABLS_tke(FT)),
+(:GABLS_tke_prescribed, APL.GABLS_tke_prescribed(FT)),
+(:GABLS_u, APL.GABLS_u(FT)),
+(:GABLS_v, APL.GABLS_v(FT)),
+(:GABLS_θ_liq_ice, APL.GABLS_θ_liq_ice(FT)),
+(:GATE_III_T, APL.GATE_III_T(FT)),
+(:GATE_III_dTdt, APL.GATE_III_dTdt(FT)),
+(:GATE_III_dqtdt, APL.GATE_III_dqtdt(FT)),
+(:GATE_III_q_tot, APL.GATE_III_q_tot(FT)),
+(:GATE_III_tke, APL.GATE_III_tke(FT)),
+(:GATE_III_u, APL.GATE_III_u(FT)),
+(:ISDAC_inv_τ_scalar, APL.ISDAC_inv_τ_scalar(FT)),
+(:ISDAC_inv_τ_wind, APL.ISDAC_inv_τ_wind(FT)),
+(:ISDAC_q_tot, APL.ISDAC_q_tot(FT)),
+(:ISDAC_subsidence, APL.ISDAC_subsidence(FT)),
+(:ISDAC_tke, APL.ISDAC_tke(FT)),
+(:ISDAC_u, APL.ISDAC_u(FT)),
+(:ISDAC_v, APL.ISDAC_v(FT)),
+(:ISDAC_θ_liq_ice, APL.ISDAC_θ_liq_ice(FT)),
+(:LifeCycleTan2018_dTdt, APL.LifeCycleTan2018_dTdt(FT)),
+(:LifeCycleTan2018_dqtdt, APL.LifeCycleTan2018_dqtdt(FT)),
+(:LifeCycleTan2018_geostrophic_u, APL.LifeCycleTan2018_geostrophic_u(FT)),
+(:LifeCycleTan2018_geostrophic_v, APL.LifeCycleTan2018_geostrophic_v(FT)),
+(:LifeCycleTan2018_q_tot, APL.LifeCycleTan2018_q_tot(FT)),
+(:LifeCycleTan2018_subsidence, APL.LifeCycleTan2018_subsidence(FT)),
+(:LifeCycleTan2018_tke, APL.LifeCycleTan2018_tke(FT)),
+(:LifeCycleTan2018_tke_prescribed, APL.LifeCycleTan2018_tke_prescribed(FT)),
+(:LifeCycleTan2018_u, APL.LifeCycleTan2018_u(FT)),
+(:LifeCycleTan2018_θ_liq_ice, APL.LifeCycleTan2018_θ_liq_ice(FT)),
+(:Nieuwstadt_tke, APL.Nieuwstadt_tke(FT)),
+(:Nieuwstadt_tke_prescribed, APL.Nieuwstadt_tke_prescribed(FT)),
+(:Nieuwstadt_u, APL.Nieuwstadt_u(FT)),
+(:Nieuwstadt_θ_liq_ice, APL.Nieuwstadt_θ_liq_ice(FT)),
+(:Rico_dTdt, APL.Rico_dTdt(FT)),
+(:Rico_dqtdt, APL.Rico_dqtdt(FT)),
+(:Rico_geostrophic_ug, APL.Rico_geostrophic_ug(FT)),
+(:Rico_geostrophic_vg, APL.Rico_geostrophic_vg(FT)),
+(:Rico_q_tot, APL.Rico_q_tot(FT)),
+(:Rico_subsidence, APL.Rico_subsidence(FT)),
+(:Rico_tke_prescribed, APL.Rico_tke_prescribed(FT)),
+(:Rico_u, APL.Rico_u(FT)),
+(:Rico_v, APL.Rico_v(FT)),
+(:Rico_θ_liq_ice, APL.Rico_θ_liq_ice(FT)),
+(:SP_geostrophic_u, APL.SP_geostrophic_u(FT)),
+(:SP_geostrophic_v, APL.SP_geostrophic_v(FT)),
+(:SP_q_tot, APL.SP_q_tot(FT)),
+(:SP_u, APL.SP_u(FT)),
+(:SP_v, APL.SP_v(FT)),
+(:SP_θ_liq_ice, APL.SP_θ_liq_ice(FT)),
+(:Soares_q_tot, APL.Soares_q_tot(FT)),
+(:Soares_tke, APL.Soares_tke(FT)),
+(:Soares_tke_prescribed, APL.Soares_tke_prescribed(FT)),
+(:Soares_u, APL.Soares_u(FT)),
+(:Soares_θ_liq_ice, APL.Soares_θ_liq_ice(FT)),
+(:TRMM_LBA_RH, APL.TRMM_LBA_RH(FT)),
+(:TRMM_LBA_T, APL.TRMM_LBA_T(FT)),
+(:TRMM_LBA_p, APL.TRMM_LBA_p(FT)),
+(:TRMM_LBA_radiation, APL.TRMM_LBA_radiation(FT)),
+(:TRMM_LBA_tke, APL.TRMM_LBA_tke(FT)),
+(:TRMM_LBA_tke_prescribed, APL.TRMM_LBA_tke_prescribed(FT)),
+(:TRMM_LBA_u, APL.TRMM_LBA_u(FT)),
+(:TRMM_LBA_v, APL.TRMM_LBA_v(FT)),
+]
diff --git a/test/runtests.jl b/test/runtests.jl
index 58dd102..699ef21 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,25 +1,40 @@
 using Test
 using AtmosphericProfilesLibrary
+import AtmosphericProfilesLibrary as APL
 
 iscallable(f) = !isempty(methods(f))
-@testset "AtmosphericProfilesLibrary" begin
-    for name in names(AtmosphericProfilesLibrary; all = true)
-        startswith(string(name), "#") && continue
-        name == :AtmosphericProfilesLibrary && continue
-        name == :TRMM_LBA_z && continue # returns an array for other profiles
-        name == :ARM_SGP_z && continue # returns an array for other profiles
-        name == :GATE_III_z && continue # returns an array for other profiles
-        name == :GATE_III_z_in && continue # returns an array for other profiles
-        name == :ARM_SGP_time && continue # returns an array for other profiles
-        name == :DryBubble_updrafts_z && continue # returns an array for other profiles
-        name == :eval && continue
-        name == :include && continue
-        prof = getproperty(AtmosphericProfilesLibrary, name)
-        iscallable(prof) || @show name
 
-        # Test that function is callable
-        @test iscallable(prof)
-        # Test for type-stability
-        @test @inferred Float32 iscallable(prof(Float32))
+function test_profile(::Type{FT}, prof, name) where {FT}
+    @debug "Testing $name"
+    if prof isa APL.TimeProfile
+        prof(1)
+        @inferred prof(Float32(1))
+    elseif prof isa APL.TimeZProfile
+        prof(1, 1)
+        @inferred prof(Float32(1), Float32(1))
+    elseif prof isa APL.ZProfile
+        prof(1)
+        @inferred prof(Float32(1))
+    elseif prof isa APL.ΠTimeZProfile
+        prof(1,1,1)
+        @inferred prof(Float32(1),Float32(1),Float32(1))
+    elseif prof isa APL.ΠZProfile
+        prof(1,1)
+        @inferred prof(Float32(1),Float32(1))
+    else
+        @show name
+        @show prof
+        error("uncaught case")
+    end
+end
+
+include("profiles.jl")
+
+@testset "AtmosphericProfilesLibrary" begin
+    for (name, prof) in profiles(Float64)
+        test_profile(Float64, prof, name)
+    end
+    for (name, prof) in profiles(Float32)
+        test_profile(Float32, prof, name)
     end
 end