Merge pull request #3567 from CliMA/aj/working_fluid_part2_RHS

Working fluid part 2 - RHS updates

config/longrun_configs/longrun_moist_baroclinic_wave.yml

@@ -10,4 +10,4 @@ dt_save_state_to_disk: "10days"
 toml: [toml/longrun_baroclinic_wave.toml]
 diagnostics:
   - short_name: [pfull, wa, va, ua, ta, rhoa, rv, hus, ke, clw, cli]
-    period: 1days
+    period: 10days

config/longrun_configs/longrun_moist_baroclinic_wave_1M.yml

@@ -0,0 +1,16 @@
+h_elem: 30
+z_elem: 43
+dz_bottom: 30.0
+dt: "60secs"
+approximate_linear_solve_iters: 3
+max_newton_iters_ode: 1
+rayleigh_sponge: true
+t_end: "15days"
+initial_condition: "MoistBaroclinicWave"
+moist: "nonequil"
+precip_model: "1M"
+dt_save_state_to_disk: "15days"
+toml: [toml/longrun_baroclinic_wave_1M.toml]
+diagnostics:
+  - short_name: [pfull, wa, va, ua, ta, rhoa, rv, hus, ke, clw, cli, husra, hussn]
+    period: 5days

src/cache/cloud_fraction.jl

@@ -2,6 +2,13 @@ import NVTX
 import StaticArrays as SA
 import ClimaCore.RecursiveApply: rzero, ⊞, ⊠
 
+"""
+    Helper function to populate the cloud diagnostics named tuple
+"""
+function make_named_tuple(t1, t2, t3)
+    return NamedTuple{(:cf, :q_liq, :q_ice)}(tuple(t1, t2, t3))
+end
+
 # TODO: write a test with scalars that are linear with z
 """
     Diagnose horizontal covariances based on vertical gradients
@@ -27,7 +34,7 @@ end
 NVTX.@annotate function set_cloud_fraction!(
     Y,
     p,
-    ::Union{EquilMoistModel, NonEquilMoistModel},
+    moist_model::Union{EquilMoistModel, NonEquilMoistModel},
     ::GridScaleCloud,
 )
     (; params) = p
@@ -48,13 +55,23 @@ NVTX.@annotate function set_cloud_fraction!(
         end
         compute_gm_mixing_length!(ᶜmixing_length, Y, p)
     end
-    @. cloud_diagnostics_tuple = NamedTuple{(:cf, :q_liq, :q_ice)}(
-        tuple(
+    if moist_model isa EquilMoistModel
+        @. cloud_diagnostics_tuple = make_named_tuple(
             ifelse(TD.has_condensate(thermo_params, ᶜts), 1, 0),
             TD.PhasePartition(thermo_params, ᶜts).liq,
             TD.PhasePartition(thermo_params, ᶜts).ice,
-        ),
-    )
+        )
+    else
+        @. cloud_diagnostics_tuple = make_named_tuple(
+            ifelse(
+                p.precomputed.ᶜspecific.q_liq + p.precomputed.ᶜspecific.q_ice > 0,
+                1,
+                0,
+            ),
+            p.precomputed.ᶜspecific.q_liq,
+            p.precomputed.ᶜspecific.q_ice,
+        )
+    end
 end
 NVTX.@annotate function set_cloud_fraction!(
     Y,

src/cache/precomputed_quantities.jl

@@ -319,30 +319,36 @@ function thermo_state(
     return get_ts(ρ, p, θ, e_int, q_tot, q_pt)
 end
 
-function thermo_vars(moisture_model, specific, K, Φ)
+function thermo_vars(moisture_model, precip_model, specific, K, Φ)
     energy_var = (; e_int = specific.e_tot - K - Φ)
     moisture_var = if moisture_model isa DryModel
         (;)
     elseif moisture_model isa EquilMoistModel
         (; specific.q_tot)
     elseif moisture_model isa NonEquilMoistModel
-        q_pt_args = (specific.q_tot, specific.q_liq, specific.q_ice)
+        q_pt_args = (
+            specific.q_tot,
+            specific.q_liq + specific.q_rai,
+            specific.q_ice + specific.q_sno,
+        )
         (; q_pt = TD.PhasePartition(q_pt_args...))
     end
     return (; energy_var..., moisture_var...)
 end
 
-ts_gs(thermo_params, moisture_model, specific, K, Φ, ρ) = thermo_state(
-    thermo_params;
-    thermo_vars(moisture_model, specific, K, Φ)...,
-    ρ,
-)
+ts_gs(thermo_params, moisture_model, precip_model, specific, K, Φ, ρ) =
+    thermo_state(
+        thermo_params;
+        thermo_vars(moisture_model, precip_model, specific, K, Φ)...,
+        ρ,
+    )
 
-ts_sgs(thermo_params, moisture_model, specific, K, Φ, p) = thermo_state(
-    thermo_params;
-    thermo_vars(moisture_model, specific, K, Φ)...,
-    p,
-)
+ts_sgs(thermo_params, moisture_model, precip_model, specific, K, Φ, p) =
+    thermo_state(
+        thermo_params;
+        thermo_vars(moisture_model, precip_model, specific, K, Φ)...,
+        p,
+    )
 
 function eddy_diffusivity_coefficient_H(D₀, H, z_sfc, z)
     return D₀ * exp(-(z - z_sfc) / H)
@@ -476,7 +482,7 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
     (; call_cloud_diagnostics_per_stage) = p.atmos
     thermo_params = CAP.thermodynamics_params(p.params)
     n = n_mass_flux_subdomains(turbconv_model)
-    thermo_args = (thermo_params, moisture_model)
+    thermo_args = (thermo_params, moisture_model, precip_model)
     (; ᶜΦ) = p.core
     (; ᶜspecific, ᶜu, ᶠu³, ᶜK, ᶜts, ᶜp) = p.precomputed
     ᶠuₕ³ = p.scratch.ᶠtemp_CT3

src/initial_conditions/initial_conditions.jl

@@ -381,9 +381,9 @@ end
 """
     overwrite_initial_conditions!(initial_condition, args...)
 
-Do-nothing fallback method for the operation overwriting initial conditions 
-(this functionality required in instances where we interpolate initial conditions from NetCDF files). 
-Future work may revisit this design choice. 
+Do-nothing fallback method for the operation overwriting initial conditions
+(this functionality required in instances where we interpolate initial conditions from NetCDF files).
+Future work may revisit this design choice.
 """
 function overwrite_initial_conditions!(
     initial_condition::InitialCondition,
@@ -419,26 +419,26 @@ function overwrite_initial_conditions!(
     @info "Overwriting initial conditions with data from file $(file_path)"
     center_space = Fields.axes(Y.c)
     face_space = Fields.axes(Y.f)
-    # Using surface pressure, air temperature and specific humidity 
-    # from the dataset, compute air pressure. 
+    # Using surface pressure, air temperature and specific humidity
+    # from the dataset, compute air pressure.
     p_sfc = Fields.level(
         SpaceVaryingInputs.SpaceVaryingInput(file_path, "p", face_space),
         Fields.half,
     )
     ᶜT = SpaceVaryingInputs.SpaceVaryingInput(file_path, "t", center_space)
     ᶜq_tot = SpaceVaryingInputs.SpaceVaryingInput(file_path, "q", center_space)
 
-    # With the known temperature (ᶜT) and moisture (ᶜq_tot) profile, 
+    # With the known temperature (ᶜT) and moisture (ᶜq_tot) profile,
     # recompute the pressure levels assuming hydrostatic balance is maintained.
-    # Uses the ClimaCore `column_integral_indefinite!` function to solve 
+    # Uses the ClimaCore `column_integral_indefinite!` function to solve
     # ∂(ln𝑝)/∂z = -g/(Rₘ(q)T), where
     # p is the local pressure
     # g is the gravitational constant
     # q is the specific humidity
     # Rₘ is the gas constant for moist air
     # T is the air temperature
     # p is then updated with the integral result, given p_sfc,
-    # following which the thermodynamic state is constructed. 
+    # following which the thermodynamic state is constructed.
     ᶜ∂lnp∂z = @. -thermo_params.grav /
        (TD.gas_constant_air(thermo_params, TD.PhasePartition(ᶜq_tot)) * ᶜT)
     ᶠlnp_over_psfc = zeros(face_space)
@@ -1261,7 +1261,7 @@ function (initial_condition::PrecipitatingColumn)(params)
             thermo_params,
             p(z),
             θ(z),
-            TD.PhasePartition(q_tot(z), qₗ(z), qᵢ(z)),
+            TD.PhasePartition(q_tot(z), qₗ(z) + qᵣ(z), qᵢ(z) + qₛ(z)),
         )
         return LocalState(;
             params,

src/parameterized_tendencies/microphysics/cloud_condensate.jl

@@ -27,9 +27,11 @@ function cloud_condensate_tendency!(
 )
     (; ᶜts) = p.precomputed
     (; params, dt) = p
+    (; q_rai, q_sno) = p.precomputed.ᶜspecific
+    FT = eltype(params)
     thp = CAP.thermodynamics_params(params)
     cmc = CAP.microphysics_cloud_params(params)
 
-    @. Yₜ.c.ρq_liq += cloud_sources(cmc.liquid, thp, ᶜts, dt)
-    @. Yₜ.c.ρq_ice += cloud_sources(cmc.ice, thp, ᶜts, dt)
+    @. Yₜ.c.ρq_liq += cloud_sources(cmc.liquid, thp, ᶜts, q_rai, dt)
+    @. Yₜ.c.ρq_ice += cloud_sources(cmc.ice, thp, ᶜts, q_sno, dt)
 end