Adds a moisture fixer tendency for small moisture tracers

Author: trontrytel

config/default_configs/default_config.yml

@@ -172,6 +172,9 @@ c_amd:
 smagorinsky_lilly:
   help: "Smagorinsky-Lilly diffusive closure [`nothing` (default), `UVW`, `UV`, `W`, `UV_W_decoupled`]"
   value: ~
+moisture_fixer:
+  help: "A flag to switch on a grid mean tendency that moves mass from vapor to adjust small negative tracer values [false (default), true]"
+  value: false
 bubble:
   help: "Enable bubble correction for more accurate surface areas"
   value: false

config/model_configs/baroclinic_wave_nonequil_conservation_source.yml

@@ -2,6 +2,7 @@ FLOAT_TYPE: "Float32"
 initial_condition: "MoistBaroclinicWave"
 t_end: "5days"
 moist: "nonequil"
+moisture_fixer: true
 surface_setup: DefaultMoninObukhov
 prognostic_surface: "SlabOceanSST"
 rad: "clearsky"

docs/src/microphysics.md

@@ -140,6 +140,17 @@ The magnitude of diffusion acting on precipitation tracers can be scaled using t
   `tracer_vertical_diffusion_factor`.
 There is no such scaling applied when using the Smagorinsky-Lilly or AMD models.
 
+### Moisture Fixer
+
+The moisture fixer is an optional grid-mean tendency correction that can be enabled for
+  nonequilibrium moisture models with 1-moment microphysics.
+When enabled (via the `moisture_fixer` configuration option), the fixer adjusts small negative
+  tracer values by transferring mass from water vapor to the affected microphysics tracers
+  (cloud liquid, cloud ice, rain, and snow).
+This provides an additional safeguard against numerical errors that may produce negative tracer values,
+  complementing the limiters and diffusion schemes described above.
+The fixer operates on the grid-mean tendencies and conserves total water mass.
+
 ## Aerosol Activation for 2-Moment Microphysics
 
 Aerosol activation uses functions from the [CloudMicrophysics.jl](https://github.com/CliMA/CloudMicrophysics.jl) library, based on the Abdul-Razzak and Ghan (ARG) parameterization. ARG predicts the number of activated cloud droplets assuming a parcel of clear air rising adiabatically. This formulation is traditionally applied only at cloud base, where the maximum supersaturation typically occurs.

src/ClimaAtmos.jl

@@ -115,6 +115,10 @@ include(
         "microphysics_wrappers.jl",
     ),
 )
+include(
+    joinpath(
+        "parameterized_tendencies", "microphysics", "moisture_fixer.jl"),
+)
 include(
     joinpath("prognostic_equations", "vertical_diffusion_boundary_layer.jl"),
 )

src/parameterized_tendencies/microphysics/microphysics_wrappers.jl

@@ -29,6 +29,15 @@ function limit(q, dt, n::Int)
     return q / FT(dt) / n
 end
 
+function moisture_fixer(q, qᵥ, dt)
+    FT = eltype(q)
+    return triangle_inequality_limiter(
+        -min(FT(0), q / dt),
+        limit(qᵥ, dt, 5),
+        FT(0),
+    )
+end
+
 """
     cloud_sources(cm_params, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Tₐ, dt)
 

src/parameterized_tendencies/microphysics/moisture_fixer.jl

@@ -0,0 +1,35 @@
+import ClimaCore.MatrixFields as MF
+import CloudMicrophysics.ThermodynamicsInterface as TDI
+
+moisture_fixer_tendency!(Yₜ, Y, p, t, _, _) = nothing
+
+function moisture_fixer_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    ::NonEquilMoistModel,
+    ::Union{Microphysics1Moment, Microphysics2Moment},
+)
+    if p.atmos.water.moisture_fixer
+        moisture_species = (
+            MF.@name(c.ρq_liq), MF.@name(c.ρq_ice),
+            MF.@name(c.ρq_rai), MF.@name(c.ρq_sno),
+        )
+        qᵥ = @. lazy(
+            specific(
+                TDI.q_vap(Y.c.ρq_tot, Y.c.ρq_liq, Y.c.ρq_ice, Y.c.ρq_rai, Y.c.ρq_sno),
+                Y.c.ρ,
+            ),
+        )
+
+        MF.unrolled_foreach(moisture_species) do (ρq_name)
+            ᶜρq = MF.get_field(Y, ρq_name)
+            ᶜρqₜ = MF.get_field(Yₜ, ρq_name)
+            ᶜq = @. lazy(specific(ᶜρq, Y.c.ρ))
+            # Increase the grid mean small tracers if negative,
+            # using mass from grid mean vapor.
+            @. ᶜρqₜ += Y.c.ρ * moisture_fixer(ᶜq, qᵥ, p.dt)
+        end
+    end
+end

src/prognostic_equations/remaining_tendency.jl

@@ -7,7 +7,7 @@ state vector `Y`.
 This function follows a sequence:
 1. Prepares hyperdiffusion tendencies for tracers (stored in `Yₜ_lim`).
 2. Prepares hyperdiffusion tendencies for other state variables (e.g., momentum, energy, stored in `Yₜ`).
-3. If Direct Stiffness Summation (DSS) is required and hyperdiffusion is active, performs DSS on the 
+3. If Direct Stiffness Summation (DSS) is required and hyperdiffusion is active, performs DSS on the
    prepared hyperdiffusion tendencies.
 4. Applies the (potentially DSSed) hyperdiffusion tendencies to `Yₜ_lim` and `Yₜ`.
 
@@ -356,7 +356,11 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
     horizontal_amd_tendency!(Yₜ, Y, p, t, amd)
     vertical_amd_tendency!(Yₜ, Y, p, t, amd)
 
-    # NOTE: This will zero out all momentum tendencies in the EDMFX advection test, 
+    # Optional tendency to bring negative small tracers back from negative
+    # at the cost of water vapor.
+    moisture_fixer_tendency!(Yₜ, Y, p, t, moisture_model, microphysics_model)
+
+    # NOTE: This will zero out all momentum tendencies in the EDMFX advection test,
     # where velocities do not evolve
     # DO NOT add additional velocity tendencies after this function
     zero_velocity_tendency!(Yₜ, Y, p, t)

src/solver/type_getters.jl

@@ -122,6 +122,7 @@ function get_atmos(config::AtmosConfig, params)
         noneq_cloud_formation_mode = implicit_noneq_cloud_formation ?
                                      Implicit() : Explicit(),
         call_cloud_diagnostics_per_stage,
+        moisture_fixer = parsed_args["moisture_fixer"],
 
         # SCMSetup - Single-Column Model components
         subsidence = get_subsidence_model(parsed_args, radiation_mode, FT),

src/solver/types.jl

@@ -610,6 +610,7 @@ Base.@kwdef struct AtmosWater{MM, PM, CM, NCFM, CCDPS}
     cloud_model::CM = nothing
     noneq_cloud_formation_mode::NCFM = nothing
     call_cloud_diagnostics_per_stage::CCDPS = nothing
+    moisture_fixer::Bool = false
 end
 
 """