Merge pull request #690 from CliMA/aj/derivatives_for_2M

Add 2M derivatives and tests

Author: trontrytel

Project.toml

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

docs/src/API.md

@@ -88,6 +88,7 @@ Microphysics2M.rain_breakup
 Microphysics2M.rain_self_collection_and_breakup
 Microphysics2M.rain_terminal_velocity
 Microphysics2M.rain_evaporation
+Microphysics2M.∂rain_evaporation_∂N_rai_∂q_rai
 Microphysics2M.conv_q_lcl_to_q_rai
 Microphysics2M.number_increase_for_mass_limit
 Microphysics2M.number_decrease_for_mass_limit

src/BulkMicrophysicsTendencies.jl

@@ -369,6 +369,95 @@ as `bulk_microphysics_tendencies`.
     return (; ∂tendency_∂q_lcl, ∂tendency_∂q_icl, ∂tendency_∂q_rai, ∂tendency_∂q_sno)
 end
 
+# --- 2-Moment Microphysics derivatives ---
+
+"""
+    bulk_microphysics_derivatives(
+        ::Microphysics2Moment,
+        mp,
+        tps,
+        ρ,
+        T,
+        q_tot,
+        q_lcl,
+        q_icl,
+        q_rai,
+        q_sno,
+        n_lcl,
+        n_rai,
+    )
+
+Compute 2-moment microphysics tendency derivatives in one fused call.
+
+Returns leading-order derivatives of species tendencies w.r.t. their own
+specific content (q) and number (n). Rain uses the 2-moment rain evaporation
+derivatives; snow and cloud formation derivatives are zero for now.
+
+# Arguments
+- `mp`: Microphysics2MParams (warm rain; ice optional)
+- `n_lcl`, `n_rai`: Cloud and rain number per kg air (1/kg); N_rai = ρ * n_rai used for rain evaporation derivative
+
+# Returns
+`NamedTuple` with fields: `∂tendency_∂q_lcl`, `∂tendency_∂q_icl`, `∂tendency_∂q_rai`, `∂tendency_∂q_sno`, `∂tendency_∂n_lcl`, `∂tendency_∂n_rai`
+"""
+@inline function bulk_microphysics_derivatives(
+    ::Microphysics2Moment,
+    mp::CMP.Microphysics2MParams{FT, WR, ICE},
+    tps,
+    ρ,
+    T,
+    q_tot,
+    q_lcl,
+    q_icl,
+    q_rai,
+    q_sno,
+    n_lcl,
+    n_rai,
+) where {FT, WR, ICE}
+    ρ = UT.clamp_to_nonneg(ρ)
+    q_tot = UT.clamp_to_nonneg(q_tot)
+    q_lcl = UT.clamp_to_nonneg(q_lcl)
+    q_icl = UT.clamp_to_nonneg(q_icl)
+    q_rai = UT.clamp_to_nonneg(q_rai)
+    q_sno = UT.clamp_to_nonneg(q_sno)
+    n_lcl = UT.clamp_to_nonneg(n_lcl)
+    n_rai = UT.clamp_to_nonneg(n_rai)
+
+    sb = mp.warm_rain.seifert_beheng
+    aps = mp.warm_rain.air_properties
+    N_rai = ρ * n_rai
+
+    # TODO: Cloud formation — 2M bulk_microphysics_tendencies does not call cloud formation yet;
+    # once it does, set these from MM2015 (same as 1M) via ∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld.
+    ∂tendency_∂q_lcl = zero(ρ)
+    ∂tendency_∂q_icl = zero(ρ)
+
+    # 2-moment rain evaporation derivatives (return order: ∂n_rai, ∂q_rai, matching rain_evaporation)
+    (; ∂N_rai, ∂q_rai) =
+        CM2.∂rain_evaporation_∂N_rai_∂q_rai(sb, aps, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, N_rai, T)
+    ∂tendency_∂q_rai = ∂q_rai
+    ∂tendency_∂n_rai = ∂N_rai / ρ
+
+    # TODO: 2M scheme also applies number/mass adjustment terms (e.g. number_increase_for_mass_limit,
+    # number_decrease_for_mass_limit) to keep the assumed size distribution stable; consider adding
+    # their derivatives here later.
+
+    # Snow derivative zero for now
+    ∂tendency_∂q_sno = zero(ρ)
+
+    # Cloud number derivative zero for now
+    ∂tendency_∂n_lcl = zero(ρ)
+
+    return (;
+        ∂tendency_∂q_lcl,
+        ∂tendency_∂q_icl,
+        ∂tendency_∂q_rai,
+        ∂tendency_∂q_sno,
+        ∂tendency_∂n_lcl,
+        ∂tendency_∂n_rai,
+    )
+end
+
 # --- 0-Moment Microphysics ---
 
 """

src/Microphysics2M.jl

@@ -27,6 +27,7 @@ export autoconversion,
     rain_terminal_velocity,
     conv_q_lcl_to_q_rai,
     rain_evaporation,
+    ∂rain_evaporation_∂N_rai_∂q_rai,
     rain_self_collection,
     rain_breakup,
     rain_self_collection_and_breakup,
@@ -796,6 +797,39 @@ function rain_evaporation(
     return (; evap_rate_0, evap_rate_1)
 end
 
+"""
+    ∂rain_evaporation_∂N_rai_∂q_rai(sb, aps, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, N_rai, T)
+
+Returns the leading-order derivatives of the rain evaporation tendencies with
+respect to rain specific content `q_rai` and rain number concentration N_rai.
+
+Uses the same approximation pattern as
+`Microphysics1M.∂evaporation_sublimation_∂q_precip`:
+- ∂(evap_rate_0/ρ)/∂N_rai ≈ evap_rate_0 / N_rai  (number tendency, first)
+- ∂(evap_rate_1)/∂q_rai ≈ evap_rate_1 / q_rai  (mass tendency, second)
+
+# Returns
+`NamedTuple` with fields `(; ∂tendency_∂N_rai, ∂tendnecy_∂q_rai)`.
+"""
+@inline function ∂rain_evaporation_∂N_rai_∂q_rai(
+    sb::CMP.SB2006{FT},
+    aps::CMP.AirProperties,
+    tps::TDI.PS,
+    q_tot,
+    q_lcl,
+    q_icl,
+    q_rai,
+    q_sno,
+    ρ,
+    N_rai,
+    T,
+) where {FT}
+    result = rain_evaporation(sb, aps, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, N_rai, T)
+    ∂N_rai = N_rai > UT.ϵ_numerics_2M_N(FT) ? result.evap_rate_0 / N_rai : zero(result.evap_rate_0)
+    ∂q_rai = q_rai > UT.ϵ_numerics_2M_M(FT) ? result.evap_rate_1 / q_rai : zero(result.evap_rate_1)
+    return (; ∂N_rai, ∂q_rai)
+end
+
 """
     number_increase_for_mass_limit(numadj, x_max, q, ρ, N)