Cleanup

Author: tapios

src/air_saturation_functions.jl

@@ -390,9 +390,19 @@ Derivative of saturation vapor specific humidity with respect to temperature.
 # 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.
+Computed via the Clausius-Clapeyron relation:
+```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
+```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/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