Cleanup

Author: tapios

docs/src/API.md

@@ -112,12 +112,12 @@ air_pressure_given_θ
 ```@docs
 saturation_vapor_pressure
 q_vap_saturation
+∂q_vap_sat_∂T
 supersaturation
 saturation_excess
 q_vap_from_RH
 q_vap_from_p_vap
 q_vap_saturation_from_pressure
-∂q_vap_sat_∂T
 ```
 
 ### Saturation Adjustment

src/air_energies.jl

@@ -282,21 +282,32 @@ The total energy is `e_tot = e_int + e_kin + e_pot`.
 end
 
 """
-    total_enthalpy(e_tot, R_m, T)
+    total_enthalpy(param_set, e_tot, T, q_tot=0, q_liq=0, q_ice=0)
 
-The total specific enthalpy.
+The total enthalpy per unit mass, `h_tot = e_tot + R_m T`.
 
 # Arguments
- - `e_tot`: total specific energy [J/kg]
- - `R_m`: gas constant of moist air [J/(kg·K)], see [`gas_constant_air`](@ref)
- - `T`: air temperature [K]
+ - `param_set`: thermodynamics parameter set, see [`Thermodynamics`](@ref)
+ - `e_tot`: total specific energy [J/kg], see [`total_energy`](@ref)
+ - `T`: temperature [K]
+ - `q_tot`: total specific humidity [kg/kg]
+ - `q_liq`: liquid specific humidity [kg/kg]
+ - `q_ice`: ice specific humidity [kg/kg]
 
 # Returns
  - `h_tot`: total specific enthalpy [J/kg]
 
-The total enthalpy is computed as `h_tot = e_tot + R_m T`.
+In the dry limit (`q_tot = q_liq = q_ice = 0`, the default), this reduces to the dry-air expression.
 """
-@inline function total_enthalpy(e_tot, R_m, T)
+@inline function total_enthalpy(
+    param_set::APS,
+    e_tot,
+    T,
+    q_tot = 0,
+    q_liq = 0,
+    q_ice = 0,
+)
+    R_m = gas_constant_air(param_set, q_tot, q_liq, q_ice)
     return e_tot + R_m * T
 end
 

src/air_saturation_functions.jl

@@ -388,11 +388,21 @@ Derivative of saturation vapor specific humidity with respect to temperature.
 - `phase`: (Optional) Phase (Liquid() or Ice()).
 
 # Returns
-- `∂q_v^* / ∂T`: Derivative of saturation specific humidity with respect to temperature [kg/kg/K].
-
-Computed via the Clausius-Clapeyron relation: `∂q_sat/∂T = q_sat * (L / (Rv * T^2) - 1 / T)`, which 
-follows from `q_sat = p_sat / (ρ * R_v * T)` and the Clausius-Clapeyron relation 
-`∂ln(p_sat)/∂T = L / (Rv * T^2)` by differentiation with respect to `T` while holding `ρ` constant.
+- `∂q_vap_sat / ∂T`: Derivative of saturation specific humidity with respect to temperature [kg/kg/K].
+
+Computed as
+```math
+\\frac{\\partial q_v^*}{\\partial T} = q_v^* \\left( \\frac{L}{R_v T^2} - \\frac{1}{T} \\right),
+```
+which follows from the definition of saturation specific humidity and the ideal gas law,
+```math
+q_v^* = \\frac{p_v^*}{\\rho R_v T},
+```
+(where `` q_v^* `` is `q_vap_sat` and `` p_v^* `` is saturation vapor pressure) and the Clausius-Clapeyron relation
+```math
+\\frac{\\partial \\ln p_v^*}{\\partial T} = \\frac{L}{R_v T^2}
+```
+by differentiation with respect to `` T `` while holding `` \\rho `` constant.
 
 If `q_liq` and `q_ice` are provided, the saturation specific humidity derivative is computed
 assuming a phase mixture defined by the liquid fraction (see [`liquid_fraction`](@ref)).
@@ -407,16 +417,14 @@ parameterization (see [`liquid_fraction_ramp`](@ref)).
     q_vap_sat = q_vap_saturation(param_set, T, ρ)
     λ = liquid_fraction_ramp(param_set, T)
     L = latent_heat_mixed(param_set, T, λ)
-    R_v = TP.R_v(param_set)
-    return q_vap_sat * (L / (R_v * T^2) - 1 / T)
+    return ∂q_vap_sat_∂T_from_L(param_set, q_vap_sat, L, T)
 end
 
 @inline function ∂q_vap_sat_∂T(param_set::APS, T, ρ, q_liq, q_ice)
     q_vap_sat = q_vap_saturation(param_set, T, ρ, q_liq, q_ice)
     λ = liquid_fraction(param_set, T, q_liq, q_ice)
     L = latent_heat_mixed(param_set, T, λ)
-    R_v = TP.R_v(param_set)
-    return q_vap_sat * (L / (R_v * T^2) - 1 / T)
+    return ∂q_vap_sat_∂T_from_L(param_set, q_vap_sat, L, T)
 end
 
 @inline function ∂q_vap_sat_∂T(param_set::APS, T, ρ, phase::Phase)
@@ -426,6 +434,10 @@ end
         latent_heat_vapor(param_set, T),
         latent_heat_sublim(param_set, T),
     )
+    return ∂q_vap_sat_∂T_from_L(param_set, q_vap_sat, L, T)
+end
+
+@inline function ∂q_vap_sat_∂T_from_L(param_set::APS, q_vap_sat, L, T)
     R_v = TP.R_v(param_set)
     return q_vap_sat * (L / (R_v * T^2) - 1 / T)
 end
@@ -628,4 +640,3 @@ If the specific humidities are not given, the result is the saturation vapor pre
         vapor_pressure_deficit(param_set, T, p, q_tot, q_liq, q_ice, Ice()),
     )
 end
-

src/depr_phase_partition_methods.jl

@@ -583,28 +583,22 @@ and total specific humidity.
 end
 
 """
-    total_enthalpy(param_set, e_tot, T[, q::PhasePartition])
+    total_enthalpy(param_set, e_tot, T, q::PhasePartition)
 
 The total specific enthalpy, defined as `e_tot + R_m * T`.
 
  - `param_set` a thermodynamics parameter set, see the [`Thermodynamics`](@ref) for more details
  - `e_tot` total specific energy
  - `T` temperature
-
-and, optionally,
-
  - `q` [`PhasePartition`](@ref). 
- 
-When `q` is not provided, the results are for dry air.
 """
 @inline function total_enthalpy(
     param_set::APS,
     e_tot,
     T,
-    q::PhasePartition = q_pt_0(param_set),
+    q::PhasePartition,
 )
-    R_m = gas_constant_air(param_set, q)
-    return total_enthalpy(e_tot, R_m, T)
+    return total_enthalpy(param_set, e_tot, T, q.tot, q.liq, q.ice)
 end
 
 # -------------------------------------------------------------------------

src/depr_state_methods.jl

@@ -326,7 +326,7 @@ The total specific enthalpy.
 @inline function total_enthalpy(param_set::APS, ts::ThermodynamicState, e_tot)
     R_m = gas_constant_air(param_set, ts)
     T = air_temperature(param_set, ts)
-    return total_enthalpy(e_tot, R_m, T)
+    return e_tot + R_m * T
 end
 
 """

src/saturation_adjustment.jl

@@ -1215,4 +1215,3 @@ calculating necessary intermediate variables.
     e_int_sat = internal_energy_sat(param_set, T, ρ, q_tot)
     return ∂e_int_∂T(param_set, T, e_int_sat, ρ, q_tot, q_vap_sat)
 end
-

test/ad_tests.jl

@@ -295,4 +295,41 @@ using Thermodynamics: ρe, pe, ph
             end
         end
     end
+
+    @testset "∂q_vap_sat_∂T method consistency" begin
+        # Test that the different method signatures return consistent values
+        # Above freezing: Phase(Liquid()) should match (q_liq=0, q_ice=0) which defaults to liquid
+        # Below freezing: Phase(Ice()) should match (q_liq=0, q_ice=0) which defaults to ice
+
+        ρ = FT(1.0)
+
+        # Above freezing (T > T_freeze): expects liquid saturation
+        T_above = FT(290.0)
+        dq_dT_default = TD.∂q_vap_sat_∂T(param_set, T_above, ρ)
+        dq_dT_explicit = TD.∂q_vap_sat_∂T(param_set, T_above, ρ, FT(0), FT(0))
+        dq_dT_liquid = TD.∂q_vap_sat_∂T(param_set, T_above, ρ, TD.Liquid())
+
+        @test isapprox(dq_dT_default, dq_dT_explicit; rtol = FT(1e-10))
+        @test isapprox(dq_dT_default, dq_dT_liquid; rtol = FT(1e-10))
+
+        # Below freezing (T < T_icenuc): expects ice saturation  
+        T_below = FT(220.0)
+        dq_dT_default_cold = TD.∂q_vap_sat_∂T(param_set, T_below, ρ)
+        dq_dT_explicit_cold = TD.∂q_vap_sat_∂T(param_set, T_below, ρ, FT(0), FT(0))
+        dq_dT_ice = TD.∂q_vap_sat_∂T(param_set, T_below, ρ, TD.Ice())
+
+        @test isapprox(dq_dT_default_cold, dq_dT_explicit_cold; rtol = FT(1e-10))
+        @test isapprox(dq_dT_default_cold, dq_dT_ice; rtol = FT(1e-10))
+
+        # Mixed phase region (T_icenuc < T < T_freeze): 
+        # Default method uses liquid_fraction_ramp, should differ from pure phases
+        T_mixed = FT(260.0)
+        dq_dT_mixed_default = TD.∂q_vap_sat_∂T(param_set, T_mixed, ρ)
+        dq_dT_mixed_liquid = TD.∂q_vap_sat_∂T(param_set, T_mixed, ρ, TD.Liquid())
+        dq_dT_mixed_ice = TD.∂q_vap_sat_∂T(param_set, T_mixed, ρ, TD.Ice())
+
+        # In mixed phase, default should be between pure liquid and pure ice
+        @test min(dq_dT_mixed_liquid, dq_dT_mixed_ice) <= dq_dT_mixed_default <=
+              max(dq_dT_mixed_liquid, dq_dT_mixed_ice)
+    end
 end

test/correctness.jl

@@ -261,6 +261,13 @@ using the non-deprecated functional API (no `PhasePartition`/state types).
             s_vd =
                 TD.virtual_dry_static_energy(param_set, T, e_pot, q_tot, q_liq, q_ice)
             @test s_vd ≈ cp_d * (T_virt - T_0) + e_pot
+
+            # Total enthalpy: h_tot = e_tot + R_m * T
+            e_kin = FT(50)  # Kinetic energy [J/kg]
+            e_tot = TD.total_energy(param_set, e_kin, e_pot, T, q_tot, q_liq, q_ice)
+            h_tot = TD.total_enthalpy(param_set, e_tot, T, q_tot, q_liq, q_ice)
+            R_m = TD.gas_constant_air(param_set, q_tot, q_liq, q_ice)
+            @test h_tot ≈ e_tot + R_m * T
         end
 
         @testset "Humidity definitions (\$FT)" begin