Wip on adding cloud condensate sedimentation

config/model_configs/sphere_aquaplanet_rhoe_nonequilmoist_allsky.yml

@@ -9,7 +9,7 @@ implicit_diffusion: true
 approximate_linear_solve_iters: 2
 cloud_model: "grid_scale"
 moist: "nonequil"
-precip_model: "nothing"
+precip_model: "1M"
 rad: "allskywithclear"
 insolation: "timevarying"
 rayleigh_sponge: true

examples/Manifest.toml

@@ -440,9 +440,9 @@ version = "0.1.13"
 
 [[deps.CloudMicrophysics]]
 deps = ["ClimaParams", "DocStringExtensions", "ForwardDiff", "HCubature", "LazyArtifacts", "QuadGK", "RootSolvers", "SpecialFunctions", "Thermodynamics"]
-git-tree-sha1 = "8f93d0c730dd4f41fade18a1954abf1eb5ac69ce"
+git-tree-sha1 = "e3b2ae212b68aea23d11c03300abb1268a56af87"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
-version = "0.22.4"
+version = "0.22.5"
 
     [deps.CloudMicrophysics.extensions]
     EmulatorModelsExt = ["DataFrames", "MLJ"]

src/ClimaAtmos.jl

@@ -29,6 +29,7 @@ include(joinpath("parameterized_tendencies", "radiation", "radiation.jl"))
 
 include(joinpath("cache", "prognostic_edmf_precomputed_quantities.jl"))
 include(joinpath("cache", "diagnostic_edmf_precomputed_quantities.jl"))
+include(joinpath("cache", "sedimentation_precomputed_quantities.jl"))
 include(joinpath("cache", "precipitation_precomputed_quantities.jl"))
 include(joinpath("cache", "precomputed_quantities.jl"))
 include(joinpath("cache", "cloud_fraction.jl"))

src/cache/precomputed_quantities.jl

@@ -148,6 +148,12 @@ function precomputed_quantities(Y, atmos)
     else
         (;)
     end
+    sedimentation_quantities =
+        atmos.moisture_model isa NonEquilMoistModel ?
+        (;
+            ᶜwₗ = similar(Y.c, FT),
+            ᶜwᵢ = similar(Y.c, FT),
+        ) : (;)
     precipitation_quantities =
         atmos.precip_model isa Microphysics1Moment ?
         (;
@@ -162,6 +168,7 @@ function precomputed_quantities(Y, atmos)
         advective_sgs_quantities...,
         diagnostic_sgs_quantities...,
         vert_diff_quantities...,
+        sedimentation_quantities...,
         precipitation_quantities...,
         cloud_diagnostics_tuple,
     )
@@ -502,6 +509,10 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
     #     compute_gm_mixing_length!(ᶜmixing_length, Y, p)
     # end
 
+    if moisture_model isa NonEquilMoistModel
+        set_sedimentation_precomputed_quantities!(Y, p, t)
+    end
+
     if precip_model isa Microphysics1Moment
         set_precipitation_precomputed_quantities!(Y, p, t)
     end

src/cache/sedimentation_precomputed_quantities.jl

@@ -0,0 +1,32 @@
+#####
+##### Precomputed quantities
+#####
+import CloudMicrophysics.MicrophysicsNonEq as CMNe
+
+"""
+    set_sedimentation_precomputed_quantities!(Y, p, t)
+
+Updates the sedimentation terminal velocity stored in `p`
+for the non-equilibrium microphysics scheme
+"""
+function set_sedimentation_precomputed_quantities!(Y, p, t)
+    @assert (p.atmos.moisture_model isa NonEquilMoistModel)
+
+    (; ᶜwₗ, ᶜwᵢ) = p.precomputed
+    cmc = CAP.microphysics_cloud_params(p.params)
+
+    # compute the precipitation terminal velocity [m/s]
+    @. ᶜwₗ = CMNe.terminal_velocity(
+        cmc.liquid,
+        cmc.Ch2022.rain,
+        Y.c.ρ,
+        max(0, Y.c.ρq_liq / Y.c.ρ),
+    )
+    @. ᶜwᵢ = CMNe.terminal_velocity(
+        cmc.ice,
+        cmc.Ch2022.small_ice,
+        Y.c.ρ,
+        max(0, Y.c.ρq_ice / Y.c.ρ),
+    )
+    return nothing
+end

src/parameters/create_parameters.jl

@@ -96,6 +96,7 @@ function create_parameter_set(config::AtmosConfig)
         (;
             liquid = CM.Parameters.CloudLiquid(toml_dict),
             ice = CM.Parameters.CloudIce(toml_dict),
+            Ch2022 = CM.Parameters.Chen2022VelType(toml_dict),
         )
     else
         nothing

src/prognostic_equations/implicit/implicit_solver.jl

@@ -402,6 +402,10 @@ NVTX.@annotate function Wfact!(A, Y, p, dtγ, t)
         p.precomputed.ᶜK,
         p.precomputed.ᶜts,
         p.precomputed.ᶜp,
+        (
+            p.atmos.moisture_model isa NonEquilMoistModel ?
+            (; p.precomputed.ᶜwₗ, p.precomputed.ᶜwᵢ) : (;)
+        )...,
         (
             p.atmos.precip_model isa Microphysics1Moment ?
             (; p.precomputed.ᶜwᵣ, p.precomputed.ᶜwₛ) : (;)
@@ -679,8 +683,12 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
         end
 
         ᶠlg = Fields.local_geometry_field(Y.f)
-        precip_info =
-            ((@name(c.ρq_rai), @name(ᶜwᵣ)), (@name(c.ρq_sno), @name(ᶜwₛ)))
+        precip_info = (
+            (@name(c.ρq_liq), @name(ᶜwₗ)),
+            (@name(c.ρq_ice), @name(ᶜwᵢ)),
+            (@name(c.ρq_rai), @name(ᶜwᵣ)),
+            (@name(c.ρq_sno), @name(ᶜwₛ)),
+        )
         MatrixFields.unrolled_foreach(precip_info) do (ρqₚ_name, wₚ_name)
             MatrixFields.has_field(Y, ρqₚ_name) || return
             ∂ᶜρqₚ_err_∂ᶜρqₚ = matrix[ρqₚ_name, ρqₚ_name]

src/prognostic_equations/implicit/implicit_tendency.jl

@@ -146,11 +146,23 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
         )
     end
 
+    # Advection of cloud condensate and precipitation with the mean flow
+    # is done with other passive tracers in the explicit tendency.
+    # Here we add the advection with cloud condensate and precipitation
+    # terminal velocity using downward biasing
+    # and free outflow bottom boundary condition
+    if moisture_model isa NonEquilMoistModel
+        (; ᶜwₗ, ᶜwᵢ) = p.precomputed
+        @. Yₜ.c.ρq_liq -= ᶜprecipdivᵥ(
+            ᶠwinterp(ᶜJ, Y.c.ρ) *
+            ᶠright_bias(Geometry.WVector(-(ᶜwₗ)) * ᶜspecific.q_liq),
+        )
+        @. Yₜ.c.ρq_ice -= ᶜprecipdivᵥ(
+            ᶠwinterp(ᶜJ, Y.c.ρ) *
+            ᶠright_bias(Geometry.WVector(-(ᶜwᵢ)) * ᶜspecific.q_ice),
+        )
+    end
     if precip_model isa Microphysics1Moment
-        # Advection of precipitation with the mean flow
-        # is done with other passive tracers in the explicit tendency.
-        # Here we add the advection with precipitation terminal velocity
-        # using downward biasing and free outflow bottom boundary condition
         (; ᶜwᵣ, ᶜwₛ) = p.precomputed
         @. Yₜ.c.ρq_rai -= ᶜprecipdivᵥ(
             ᶠwinterp(ᶜJ, Y.c.ρ) *