Merge pull request #684 from CliMA/ts/0M-sat-excess

Expose 0M microphysics tendency based on saturation excess; bump patch release

Author: szy21

Project.toml

@@ -1,7 +1,7 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 authors = ["Climate Modeling Alliance"]
-version = "0.31.6"
+version = "0.31.7"
 
 [deps]
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"

src/BulkMicrophysicsTendencies.jl

@@ -296,14 +296,21 @@ end
 # --- 0-Moment Microphysics ---
 
 """
-    bulk_microphysics_tendencies(
-        ::Microphysics0Moment,
-        mp,
-        tps,
-        T,
-        q_lcl,
-        q_icl,
-    )
+    _precip_energy(tps, T, q_lcl, q_icl)
+
+Internal energy of removed condensate, weighted by the liquid fraction.
+Shared helper for the two 0-moment `bulk_microphysics_tendencies` methods.
+"""
+@inline function _precip_energy(tps, T, q_lcl, q_icl)
+    λ = TDI.liquid_fraction(tps, T, q_lcl, q_icl)
+    I_liq = TDI.internal_energy_liquid(tps, T)
+    I_ice = TDI.internal_energy_ice(tps, T)
+    return λ * I_liq + (1 - λ) * I_ice
+end
+
+"""
+    bulk_microphysics_tendencies(::Microphysics0Moment, mp, tps, T, q_lcl, q_icl)
+    bulk_microphysics_tendencies(::Microphysics0Moment, mp, tps, T, q_lcl, q_icl, q_vap_sat)
 
 Compute 0-moment microphysics tendencies in one fused call.
 
@@ -313,13 +320,16 @@ Returns a NamedTuple with:
 
 Caller adds geopotential Φ for energy tendency: `e_tot = dq_tot_dt * (e_int_precip + Φ)`
 
+The first form uses the fixed condensate threshold `qc_0`;
+the second form uses the supersaturation threshold `S_0 * q_vap_sat`.
+
 # Arguments
-- `mp`: NamedTuple with microphysics parameters:
-  - `params_0M`: Parameters0M struct (contains τ_precip, qc_0 or S_0)
+- `mp`: Microphysics0MParams (contains τ_precip, qc_0, S_0)
 - `tps`: Thermodynamics parameters
 - `T`: Temperature [K]
 - `q_lcl`: Cloud liquid specific content [kg/kg]
 - `q_icl`: Cloud ice specific content [kg/kg]
+- `q_vap_sat`: (second method only) Saturation specific humidity [kg/kg]
 
 # Notes
 - Does NOT apply limiters (caller applies based on timestep)
@@ -333,22 +343,26 @@ Caller adds geopotential Φ for energy tendency: `e_tot = dq_tot_dt * (e_int_pre
     q_lcl,
     q_icl,
 )
-    # Clamp negative specific contents to zero (robustness against numerical errors)
     q_lcl = UT.clamp_to_nonneg(q_lcl)
     q_icl = UT.clamp_to_nonneg(q_icl)
+    dq_tot_dt = CM0.remove_precipitation(mp.precip, q_lcl, q_icl)
+    e_int_precip = _precip_energy(tps, T, q_lcl, q_icl)
+    return (; dq_tot_dt, e_int_precip)
+end
 
-    # Unpack microphysics parameters
-    params_0M = mp.precip
-
-    # Precipitation removal rate
-    dq_tot_dt = CM0.remove_precipitation(params_0M, q_lcl, q_icl)
-
-    # Internal energy of removed condensate (liquid fraction weighted)
-    λ = TDI.liquid_fraction(tps, T, q_lcl, q_icl)
-    I_liq = TDI.internal_energy_liquid(tps, T)
-    I_ice = TDI.internal_energy_ice(tps, T)
-    e_int_precip = λ * I_liq + (1 - λ) * I_ice
-
+@inline function bulk_microphysics_tendencies(
+    ::Microphysics0Moment,
+    mp::CMP.Microphysics0MParams,
+    tps,
+    T,
+    q_lcl,
+    q_icl,
+    q_vap_sat,
+)
+    q_lcl = UT.clamp_to_nonneg(q_lcl)
+    q_icl = UT.clamp_to_nonneg(q_icl)
+    dq_tot_dt = CM0.remove_precipitation(mp.precip, q_lcl, q_icl, q_vap_sat)
+    e_int_precip = _precip_energy(tps, T, q_lcl, q_icl)
     return (; dq_tot_dt, e_int_precip)
 end
 

test/bulk_tendencies_tests.jl

@@ -548,6 +548,53 @@ function test_bulk_microphysics_0m_tendencies(FT)
         )
         @test tendencies isa NamedTuple{(:dq_tot_dt, :e_int_precip), NTuple{2, FT}}
     end
+
+    @testset "BulkMicrophysicsTendencies 0M - S_0 precipitation removal" begin
+        ρ = FT(1.2)
+        T = T_freeze + FT(10)
+        q_lcl = FT(2e-3)
+        q_icl = FT(0)
+        q_vap_sat = TDI.saturation_vapor_specific_content_over_liquid(tps, T, ρ)
+
+        tendencies = BMT.bulk_microphysics_tendencies(
+            BMT.Microphysics0Moment(),
+            mp, tps, T, q_lcl, q_icl, q_vap_sat,
+        )
+
+        # Should be negative (removing condensate above S_0 * q_vap_sat)
+        @test tendencies.dq_tot_dt < FT(0)
+        @test isfinite(tendencies.e_int_precip)
+    end
+
+    @testset "BulkMicrophysicsTendencies 0M - S_0 below threshold" begin
+        ρ = FT(1.2)
+        T = T_freeze + FT(10)
+        q_vap_sat = TDI.saturation_vapor_specific_content_over_liquid(tps, T, ρ)
+        # Condensate below S_0 * q_vap_sat
+        q_lcl = FT(1e-8)
+        q_icl = FT(0)
+
+        tendencies = BMT.bulk_microphysics_tendencies(
+            BMT.Microphysics0Moment(),
+            mp, tps, T, q_lcl, q_icl, q_vap_sat,
+        )
+
+        @test tendencies.dq_tot_dt == FT(0)
+    end
+
+    @testset "BulkMicrophysicsTendencies 0M - S_0 type stability" begin
+        ρ = FT(1.2)
+        T = T_freeze + FT(5)
+        q_lcl = FT(1e-3)
+        q_icl = FT(5e-4)
+        q_vap_sat = TDI.saturation_vapor_specific_content_over_liquid(tps, T, ρ)
+
+        tendencies = @inferred BMT.bulk_microphysics_tendencies(
+            BMT.Microphysics0Moment(),
+            mp, tps, T, q_lcl, q_icl, q_vap_sat,
+        )
+        @test tendencies isa NamedTuple{(:dq_tot_dt, :e_int_precip), NTuple{2, FT}}
+    end
 end
 
 function test_bulk_microphysics_2m_tendencies(FT)

test/gpu_tests.jl

@@ -314,6 +314,16 @@ end
     output[i] = BMT.bulk_microphysics_tendencies(CM0M, mp, tps, T[i], ql, qi)
 end
 
+@kernel inbounds = true function test_bulk_tendencies_0m_S0_kernel!(
+    mp, tps, output, liquid_frac, qc, T, q_vap_sat,
+)
+    i = @index(Global, Linear)
+    ql = qc[i] * liquid_frac[i]
+    qi = (1 - liquid_frac[i]) * qc[i]
+    CM0M = BMT.Microphysics0Moment()
+    output[i] = BMT.bulk_microphysics_tendencies(CM0M, mp, tps, T[i], ql, qi, q_vap_sat[i])
+end
+
 @kernel inbounds = true function test_bulk_tendencies_1m_kernel!(
     mp, tps, output, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno,
 )
@@ -934,6 +944,18 @@ function test_gpu(FT)
             TT.@test isfinite(tendencies.e_int_precip)
         end
 
+        # 0M S_0 tests (with q_vap_sat)
+        (; output) = setup_output(ndrange, DT)
+        q_vap_sat = constant_data(FT(0.01); ndrange)
+        kernel_s0! = test_bulk_tendencies_0m_S0_kernel!(backend, work_groups)
+        TT.@testset "0M S_0" begin
+            kernel_s0!(mp_0m, tps, output, liquid_frac, qc, T, q_vap_sat; ndrange)
+            TT.@test allequal(Array(output))
+            tendencies = Array(output)[1]
+            TT.@test tendencies.dq_tot_dt ≤ 0
+            TT.@test isfinite(tendencies.e_int_precip)
+        end
+
 
         # 1M tests
         DT = @NamedTuple{dq_lcl_dt::FT, dq_icl_dt::FT, dq_rai_dt::FT, dq_sno_dt::FT}

test/type_stability_tests.jl

@@ -40,6 +40,27 @@ function run_type_stability_tests()
                 @test getproperty(val0, k) isa FT
             end
 
+            # --- 0-Moment (S_0 mode) ---
+            ρ_0M = fill(FT(1.2), N)
+            q_vap_sat_0M = fill(FT(0.01), N)
+            tendencies_0M_S0 =
+                BMT.bulk_microphysics_tendencies.(
+                    Ref(BMT.Microphysics0Moment()),
+                    Ref(mp0),
+                    Ref(tps),
+                    T_0M,
+                    q_lcl_0M,
+                    q_icl_0M,
+                    q_vap_sat_0M,
+                )
+
+            @test tendencies_0M_S0 isa Vector
+            @test eltype(tendencies_0M_S0) <: NamedTuple
+            val0s = tendencies_0M_S0[1]
+            for k in keys(val0s)
+                @test getproperty(val0s, k) isa FT
+            end
+
             # --- 1-Moment ---
             mp1 = CMP.Microphysics1MParams(FT)
             ρ = fill(FT(1.2), N)