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/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/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,30 @@ 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, 
+    if moisture_model isa NonEquilMoistModel &&
+       microphysics_model isa Microphysics1Moment &&
+       p.atmos.water.moisture_fixer
+        moisture_species = (@name(c.ρq_liq), @name(c.ρq_ice),
+            @name(c.ρq_rai), @name(c.ρq_sno),
+        )
+        qᵥ = @. lazy(
+            specific(Y.c.ρq_tot - Y.c.ρq_liq - Y.c.ρq_ice - Y.c.ρq_rai - Y.c.ρq_sno,
+                Y.c.ρ),
+        )
+
+        MatrixFields.unrolled_foreach(
+            moisture_species,
+        ) do (ρq_name)
+            ᶜρq = MatrixFields.get_field(Y, ρq_name)
+            ᶜρqₜ = MatrixFields.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
+
+    # 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

@@ -604,12 +604,13 @@ end
 
 Groups moisture-related models and types.
 """
-Base.@kwdef struct AtmosWater{MM, PM, CM, NCFM, CCDPS}
+Base.@kwdef struct AtmosWater{MM, PM, CM, NCFM, CCDPS, MF}
     moisture_model::MM = nothing
     microphysics_model::PM = nothing
     cloud_model::CM = nothing
     noneq_cloud_formation_mode::NCFM = nothing
     call_cloud_diagnostics_per_stage::CCDPS = nothing
+    moisture_fixer::MF = false
 end
 
 """