From ca4014ed03339c177ecc7e082938d2d2a56a33e5 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Thu, 20 Jun 2024 17:23:15 -0700 Subject: [PATCH] Wip on adding 2M microphysics to Atmos --- .../precipitation_precomputed_quantities.jl | 43 +++++- src/cache/precomputed_quantities.jl | 18 ++- src/diagnostics/core_diagnostics.jl | 67 ++++++++- src/initial_conditions/atmos_state.jl | 5 + src/initial_conditions/initial_conditions.jl | 40 ++++- src/initial_conditions/local_state.jl | 16 +- .../microphysics/microphysics_wrappers.jl | 139 ++++++++++++++++++ .../microphysics/precipitation.jl | 111 +++++++++++++- src/parameters/create_parameters.jl | 6 + src/prognostic_equations/hyperdiffusion.jl | 4 +- .../implicit/implicit_solver.jl | 11 +- .../implicit/implicit_tendency.jl | 11 ++ .../vertical_diffusion_boundary_layer.jl | 2 +- 13 files changed, 453 insertions(+), 20 deletions(-) diff --git a/src/cache/precipitation_precomputed_quantities.jl b/src/cache/precipitation_precomputed_quantities.jl index f4b25dbbe48..bf7aa3ebb19 100644 --- a/src/cache/precipitation_precomputed_quantities.jl +++ b/src/cache/precipitation_precomputed_quantities.jl @@ -2,6 +2,7 @@ ##### Precomputed quantities ##### import CloudMicrophysics.Microphysics1M as CM1 +import CloudMicrophysics.Microphysics2M as CM2 # helper function to safely get precipitation from state function qₚ(ρqₚ::FT, ρ::FT) where {FT} @@ -9,13 +10,14 @@ function qₚ(ρqₚ::FT, ρ::FT) where {FT} end """ - set_precipitation_precomputed_quantities!(Y, p, t) + set_precipitation_precomputed_quantities!(Y, p, t, precip_model) -Updates the precipitation terminal velocity stored in `p` -for the 1-moment microphysics scheme +Updates the precipitation terminal velocity and tracers stored in cache """ -function set_precipitation_precomputed_quantities!(Y, p, t) - @assert (p.atmos.precip_model isa Microphysics1Moment) +function set_precipitation_precomputed_quantities!(Y, p, t, _) + return nothing +end +function set_precipitation_precomputed_quantities!(Y, p, t, ::Microphysics1Moment) (; ᶜwᵣ, ᶜwₛ, ᶜqᵣ, ᶜqₛ) = p.precomputed @@ -40,3 +42,34 @@ function set_precipitation_precomputed_quantities!(Y, p, t) ) return nothing end + +function set_precipitation_precomputed_quantities!(Y, p, t, ::Microphysics2Moment) + + (; ᶜwᵣ, ᶜw_nᵣ, ᶜqᵣ) = p.precomputed + + cmp = CAP.microphysics_precipitation_params(p.params) + + # compute the precipitation specific humidities + @. ᶜqᵣ = qₚ(Y.c.ρq_rai, Y.c.ρ) + @. ᶜwᵣ = getindex( + CM2.rain_terminal_velocity( + cmp.SB2006, + cmp.SB2006Vel, + ᶜqᵣ, + Y.c.ρ, + Y.c.ρn_rai, + ), + 2, + ) + @. ᶜw_nᵣ = getindex( + CM2.rain_terminal_velocity( + cmp.SB2006, + cmp.SB2006Vel, + ᶜqᵣ, + Y.c.ρ, + Y.c.ρn_rai, + ), + 1, + ) + return nothing +end diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index b5f07893b8d..496ec12ae25 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -146,14 +146,22 @@ function precomputed_quantities(Y, atmos) else (;) end - precipitation_quantities = - atmos.precip_model isa Microphysics1Moment ? + precipitation_quantities if atmos.precip_model isa Microphysics1Moment (; ᶜwᵣ = similar(Y.c, FT), ᶜwₛ = similar(Y.c, FT), ᶜqᵣ = similar(Y.c, FT), ᶜqₛ = similar(Y.c, FT), - ) : (;) + ) + elseif atmos.precip_model isa Microphysics2Moment + (; + ᶜwᵣ = similar(Y.c, FT), + ᶜw_nᵣ = similar(Y.c, FT), + ᶜqᵣ = similar(Y.c, FT), + ) + else + (;) + end return (; gs_quantities..., sgs_quantities..., @@ -504,9 +512,7 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) # compute_gm_mixing_length!(ᶜmixing_length, Y, p) # end - if precip_model isa Microphysics1Moment - set_precipitation_precomputed_quantities!(Y, p, t) - end + set_precipitation_precomputed_quantities!(Y, p, t, p.atmos.precip_model) if turbconv_model isa PrognosticEDMFX set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ³, t) diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl index 520b6ac9424..38bbba2ceb6 100644 --- a/src/diagnostics/core_diagnostics.jl +++ b/src/diagnostics/core_diagnostics.jl @@ -648,7 +648,7 @@ function compute_husra!( state, cache, time, - precip_model::Microphysics1Moment, + precip_model::Union{Microphysics1Moment, Microphysics2Moment}, ) if isnothing(out) return state.c.ρq_rai ./ state.c.ρ @@ -700,6 +700,71 @@ add_diagnostic_variable!( compute! = compute_hussn!, ) +### +# Number concentrations for cloud and precipitation (3d) +### +compute_cdnc!(out, state, cache, time) = + compute_cdnc!(out, state, cache, time, cache.atmos.precip_model) +compute_cdnc!(_, _, _, _, model::T) where {T} = + error_diagnostic_variable("cdnc", model) + +function compute_cdnc!( + out, + state, + cache, + time, + precip_model::Microphysics2Moment, +) + if isnothing(out) + return state.c.ρn_liq + else + out .= state.c.ρn_liq + end +end + +add_diagnostic_variable!( + short_name = "cdnc", + long_name = "Cloud Droplet Number Concentration in liquid water clouds", + standard_name = "cloud_droplet_number_concentration", + units = "m^-3", + comments = """ + This is calculated as the number of cloud droplets in the grid cell per + volume of air. + """, + compute! = compute_cdnc!, +) + +compute_rdnc!(out, state, cache, time) = + compute_rdnc!(out, state, cache, time, cache.atmos.precip_model) +compute_rdnc!(_, _, _, _, model::T) where {T} = + error_diagnostic_variable("rdnc", model) + +function compute_rdnc!( + out, + state, + cache, + time, + precip_model::Microphysics2Moment, +) + if isnothing(out) + return state.c.ρn_rai + else + out .= state.c.ρn_rai + end +end + +add_diagnostic_variable!( + short_name = "rdnc", + long_name = "Rain Drop Number Concentration", + standard_name = "rain_drop_number_concentration", + units = "m^-3", + comments = """ + This is calculated as the number of rain drops in the grid cell per + volume of air. + """, + compute! = compute_rdnc!, +) + ### # Topography ### diff --git a/src/initial_conditions/atmos_state.jl b/src/initial_conditions/atmos_state.jl index 485b2eca078..2aa57ae69b6 100644 --- a/src/initial_conditions/atmos_state.jl +++ b/src/initial_conditions/atmos_state.jl @@ -122,6 +122,11 @@ precip_variables(ls, ::Microphysics1Moment) = (; ρq_rai = ls.ρ * ls.precip_state.q_rai, ρq_sno = ls.ρ * ls.precip_state.q_sno, ) +precip_variables(ls, ::Microphysics2Moment) = (; + ρn_liq = ls.ρ * ls.precip_state.n_liq, + ρq_rai = ls.ρ * ls.precip_state.q_rai, + ρn_rai = ls.ρ * ls.precip_state.n_rai, +) # We can use paper-based cases for LES type configurations (no TKE) # or SGS type configurations (initial TKE needed), so we do not need to assert diff --git a/src/initial_conditions/initial_conditions.jl b/src/initial_conditions/initial_conditions.jl index 87dced494bc..fd03be5ad34 100644 --- a/src/initial_conditions/initial_conditions.jl +++ b/src/initial_conditions/initial_conditions.jl @@ -426,7 +426,7 @@ function deep_atmos_baroclinic_wave_values(z, ϕ, λ, params, perturb) Ω = CAP.Omega(params) R = CAP.planet_radius(params) - # Constants from paper (See Table 1. in Ullrich et al (2014)) + # Constants from paper (See Table 1. in Ullrich et al (2014)) k = 3 # Power for temperature field T_e = FT(310) # Surface temperature at the equator T_p = FT(240) # Surface temperature at the pole @@ -1124,6 +1124,44 @@ function (initial_condition::PrecipitatingColumn)(params) end return local_state end +function (initial_condition::PrecipitatingColumn2M)(params) + FT = eltype(params) + thermo_params = CAP.thermodynamics_params(params) + p_0 = FT(101300.0) + qᵣ = prescribed_prof(FT, 2000, 5000, 1e-6) + qₗ = prescribed_prof(FT, 4000, 5500, 2e-5) + qᵢ = prescribed_prof(FT, 6000, 9000, 0) + nₗ = prescribed_prof(FT, 4000, 5500, 1e8) + nᵣ = prescribed_prof(FT, 2000, 5000, 1e6) + θ = APL.Rico_θ_liq_ice(FT) + q_tot = APL.Rico_q_tot(FT) + u = prescribed_prof(FT, 0, Inf, 0) + v = prescribed_prof(FT, 0, Inf, 0) + p = hydrostatic_pressure_profile(; thermo_params, p_0, θ, q_tot) + function local_state(local_geometry) + (; z) = local_geometry.coordinates + ts = TD.PhaseNonEquil_pθq( + thermo_params, + p(z), + θ(z), + TD.PhasePartition(q_tot(z), qₗ(z), qᵢ(z)), + ) + return LocalState(; + params, + geometry = local_geometry, + thermo_state = ts, + velocity = Geometry.UVVector(u(z), v(z)), + turbconv_state = nothing, + precip_state = PrecipState2M(; + n_liq = nₗ(z), + q_rai = qᵣ(z), + n_rai = nᵣ(z) + ), + ) + end + return local_state +end + """ GCMDriven <: InitialCondition diff --git a/src/initial_conditions/local_state.jl b/src/initial_conditions/local_state.jl index bc409240863..e0929aa605c 100644 --- a/src/initial_conditions/local_state.jl +++ b/src/initial_conditions/local_state.jl @@ -104,7 +104,7 @@ struct NoPrecipState{FT} <: PrecipState{FT} end """ PrecipState1M(; q_rai, q_sno) -Stores the values of `ρq_rai` and `ρq_sno` for the `precip_model`. +Stores the values of `q_rai` and `q_sno` for the `precip_model`. If no values are provided, they are set to zero. """ struct PrecipState1M{FT} <: PrecipState{FT} @@ -113,3 +113,17 @@ struct PrecipState1M{FT} <: PrecipState{FT} end PrecipState1M(; q_rai = 0, q_sno = 0) = PrecipState1M{typeof(q_rai)}(q_rai, q_sno) + +""" + PrecipState2M(; n_liq, q_rai, n_rai) + +Stores the values of `q_rai` and `n_rai` and `n_liq` for the `precip_model`. +If no values are provided, they are set to zero. +""" +struct PrecipState2M{FT} <: PrecipState{FT} + n_liq::FT + q_rai::FT + n_rai::FT +end +PrecipState2M(; n_liq = 0, q_rai = 0, n_rai = 0) = + PrecipState2M{typeof(q_rai)}(n_liq, q_rai, n_rai) diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl index 9c8a7642d09..ab05f2c9c4d 100644 --- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl +++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl @@ -264,6 +264,74 @@ function compute_precipitation_sources!( @. Seₜᵖ += Sᵖ * Lf(thp, ts) #! format: on end +function compute_precipitation_sources2M!( + Sᵖ, + Sqₜᵖ, + Sqᵣᵖ, + SNᵣᵖ, + SNₗᵖ, + Seₜᵖ, + ρ, + qᵣ, + Nᵣ, + Nₗ, + ts, + Φ, + dt, + mp, + thp, +) + FT = eltype(thp) + # @. Sqₜᵖ = FT(0) should work after fixing + # https://github.com/CliMA/ClimaCore.jl/issues/1786 + @. Sqₜᵖ = ρ * FT(0) + @. Sqᵣᵖ = ρ * FT(0) + @. SNᵣᵖ = ρ * FT(0) + @. SNₗᵖ = ρ * FT(0) + @. Seₜᵖ = ρ * FT(0) + + #! format: off + # rain autoconversion: q_liq -> q_rain + @. Sᵖ = min( + limit(qₗ(thp, ts), dt, 5), + CM1.conv_q_liq_to_q_rai(mp.pr.acnv1M, qₗ(thp, ts), true), + ) + @. Sqₜᵖ -= Sᵖ + @. Sqᵣᵖ += Sᵖ + @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ) + + #TODO - add sources + # # autoconversion liquid to rain (mass) + # @. S₁ = limit(q_liq, dt, CM2.autoconversion(sb2006.acnv, q_liq, q_rai, ρ, N_liq).dq_rai_dt) + # @. aux.precip_sources += to_sources(-S₁, -S₁, S₁, 0, 0, 0) + + # # autoconversion liquid to rain (number) + # @. S₂ = limit(N_liq, dt, CM2.autoconversion(sb2006.acnv, q_liq, q_rai, ρ, N_liq).dN_rai_dt, 2) + # @. aux.precip_sources += to_sources(0, 0, 0, 0, -2 * S₂, S₂) + # # liquid self_collection + # @. S₁ = -limit(N_liq, dt, -CM2.liquid_self_collection(sb2006.acnv, q_liq, ρ, -2 * S₂)) + # @. aux.precip_sources += to_sources(0, 0, 0, 0, S₁, 0) + + # # rain self_collection + # @. S₁ = CM2.rain_self_collection(sb2006.pdf_r, sb2006.self, q_rai, ρ, N_rai) + # @. aux.precip_sources += to_sources(0, 0, 0, 0, 0, -limit(N_rai, dt, -S₁)) + # # rain breakup + # @. aux.precip_sources += to_sources( + # 0, + # 0, + # 0, + # 0, + # 0, + # limit(N_rai, dt, CM2.rain_breakup(sb2006.pdf_r, sb2006.brek, q_rai, ρ, N_rai, S₁)), + # ) + + # # accretion cloud water + rain + # @. S₁ = limit(q_liq, dt, CM2.accretion(sb2006, q_liq, q_rai, ρ, N_liq).dq_rai_dt) + # @. S₂ = -limit(N_liq, dt, -CM2.accretion(sb2006, q_liq, q_rai, ρ, N_liq).dN_liq_dt) + # @. aux.precip_sources += to_sources(-S₁, -S₁, S₁, 0, S₂, 0) + + +end """ compute_precipitation_heating(Seₜᵖ, ᶜwᵣ, ᶜwₛ, ᶜu, qᵣ, qₛ, ᶜts, thp) @@ -303,6 +371,23 @@ function compute_precipitation_heating!( @. ᶜSeₜᵖ -= dot(ᶜ∇T, (ᶜu - C123(Geometry.WVector(ᶜwᵣ)))) * cᵥₗ(thp) * ᶜqᵣ @. ᶜSeₜᵖ -= dot(ᶜ∇T, (ᶜu - C123(Geometry.WVector(ᶜwₛ)))) * cᵥᵢ(thp) * ᶜqₛ end +function compute_precipitation_heating!( + ᶜSeₜᵖ, + ᶜwᵣ, + ᶜu, + ᶜqᵣ, + ᶜts, + ᶜ∇T, + thp, +) + # compute full temperature gradient + @. ᶜ∇T = CT123(ᶜgradᵥ(ᶠinterp(Tₐ(thp, ᶜts)))) + @. ᶜ∇T += CT123(gradₕ(Tₐ(thp, ᶜts))) + # dot product with effective velocity of precipitation + # (times q and specific heat) + @. ᶜSeₜᵖ -= dot(ᶜ∇T, (ᶜu - C123(Geometry.WVector(ᶜwᵣ)))) * cᵥₗ(thp) * ᶜqᵣ +end + """ compute_precipitation_sinks!(Sᵖ, Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ, ρ, qᵣ, qₛ, ts, Φ, dt, mp, thp) @@ -370,3 +455,57 @@ function compute_precipitation_sinks!( @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ) #! format: on end + +""" + compute_precipitation_sinks2M!(Sᵖ, Sqₜᵖ, Sqᵣᵖ, Snᵣᵖ, Seₜᵖ, ρ, qᵣ, Nᵣ, ts, Φ, dt, mp, thp) + + - Sᵖ - a temporary containter to help compute precipitation source terms + - Sqₜᵖ, Sqᵣᵖ, Snᵣᵖ, Seₜᵖ - cached storage for precipitation source terms + - ρ - air density + - qᵣ - rain specific humidity [kg/kg] + - Nᵣ - rain number concentration [1/m3] + - ts - thermodynamic state (see td package for details) + - Φ - geopotential + - dt - model time step + - thp, cmp - structs with thermodynamic and microphysics parameters + +Returns the q and N source terms due to precipitation sinks from the 2-moment scheme. +The specific humidity source terms are defined as defined as Δmᵣ / (m_dry + m_tot) +Also returns the total energy source term due to the microphysics processes. +""" +function compute_precipitation_sinks2M!( + Sᵖ, + Sqₜᵖ, + Sqᵣᵖ, + SNᵣᵖ, + Seₜᵖ, + ρ, + qᵣ, + Nᵣ + ts, + Φ, + dt, + mp, + thp, +) + FT = eltype(Sqₜᵖ) + rps = (mp.sb2006, mp.aps, thp) + + #! format: off + # evaporation: q_rai -> q_vap (mass) + @. Sᵖ = -min( + limit(qᵣ, dt, 5), + -CM2.rain_evaporation(rps..., PP(thp, ts), qᵣ, ρ, Nᵣ, Tₐ(thp.ts)).evap_rate_1, + ) + @. Sqₜᵖ -= Sᵖ + @. Sqᵣᵖ += Sᵖ + @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ) + + # evaporation: q_rai -> q_vap (mass) + @. Sᵖ = -min( + limit(Nᵣ, dt, 5), + -CM2.rain_evaporation(rps..., PP(thp, ts), qᵣ, ρ, Nᵣ, Tₐ(thp.ts)).evap_rate_0, + ) + @. SNᵣᵖ += Sᵖ + #! format: on +end diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl index c57d7302ccc..6920e105c8d 100644 --- a/src/parameterized_tendencies/microphysics/precipitation.jl +++ b/src/parameterized_tendencies/microphysics/precipitation.jl @@ -3,6 +3,7 @@ ##### import CloudMicrophysics.Microphysics1M as CM1 +import CloudMicrophysics.Microphysics2M as CM2 import CloudMicrophysics as CM import Thermodynamics as TD import ClimaCore.Spaces as Spaces @@ -301,7 +302,6 @@ function compute_precipitation_surface_fluxes!( precip_model::Microphysics1Moment, ) (; surface_rain_flux, surface_snow_flux) = p.precipitation - (; col_integrated_precip_energy_tendency,) = p.conservation_check (; ᶜwᵣ, ᶜwₛ, ᶜspecific) = p.precomputed (; ᶠtemp_scalar) = p.scratch @@ -429,3 +429,112 @@ function precipitation_tendency!( @. Yₜ.c.ρq_sno += Y.c.sgsʲs.:($$j).ρa * ᶜSqₛᵖʲs.:($$j) end end + +##### +##### 2-Moment without sgs scheme +##### + +function precipitation_cache(Y, precip_model::Microphysics1Moment) + FT = Spaces.undertype(axes(Y.c)) + return (; + ᶜSqₜᵖ = similar(Y.c, FT), + ᶜSqᵣᵖ = similar(Y.c, FT), + ᶜSnₗᵖ = similar(Y.c, FT), + ᶜSNᵣᵖ = similar(Y.c, FT), + ᶜSeₜᵖ = similar(Y.c, FT), + surface_rain_flux = zeros(axes(Fields.level(Y.f, half))), + ) +end + +function compute_precipitation_cache!(Y, p, ::Microphysics2Moment, _) + FT = Spaces.undertype(axes(Y.c)) + (; dt) = p + (; ᶜts, ᶜqᵣ, ᶜwᵣ, ᶜw_nᵣ, ᶜu) = p.precomputed + (; ᶜΦ) = p.core + (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSnₗᵖ, ᶜSNᵣᵖ, ᶜSeₜᵖ) = p.precipitation + + ᶜSᵖ = p.scratch.ᶜtemp_scalar + ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2 + ᶜ∇T = p.scratch.ᶜtemp_CT123 + + # get thermodynamics and 1-moment microphysics params + (; params) = p + cmp = CAP.microphysics_precipitation_params(params) + thp = CAP.thermodynamics_params(params) + + # compute precipitation source terms on the grid mean + compute_precipitation_sources2M!(...) + + # compute precipitation sinks + # (For now only done on the grid mean) + compute_precipitation_sinks2M!( + ᶜSᵖ, + ᶜSqₜᵖ, + ᶜSqᵣᵖ, + ᶜSNᵣᵖ, + ᶜSeₜᵖ, + Y.c.ρ, + ᶜqᵣ, + Y.c.ρn_rai, + ᶜts, + ᶜΦ, + dt, + cmp, + thp, + ) + + # first term of eq 36 from Raymond 2013 + compute_precipitation_heating!(ᶜSeₜᵖ, ᶜwᵣ, ᶜu, ᶜqᵣ, ᶜts, ᶜ∇T, thp) +end + +function compute_precipitation_surface_fluxes!( + Y, + p, + precip_model::Microphysics2Moment, +) + (; surface_rain_flux) = p.precipitation + (; ᶜwᵣ, ᶜspecific) = p.precomputed + + (; ᶠtemp_scalar) = p.scratch + slg = Fields.level(Fields.local_geometry_field(ᶠtemp_scalar), Fields.half) + + # Constant extrapolation: - put values from bottom cell center to bottom cell face + ˢρ = Fields.Field(Fields.field_values(Fields.level(Y.c.ρ, 1)), axes(slg)) + # For density this is equivalent with ᶠwinterp(ᶜJ, Y.c.ρ) and therefore + # consistent with the way we do vertical advection + ˢqᵣ = Fields.Field( + Fields.field_values(Fields.level(ᶜspecific.q_rai, 1)), + axes(slg), + ) + ˢwᵣ = Fields.Field(Fields.field_values(Fields.level(ᶜwᵣ, 1)), axes(slg)) + + # Project the flux to CT3 vector and convert to physical units. + @. surface_rain_flux = + -projected_vector_data(CT3, ˢρ * ˢqᵣ * Geometry.WVector(ˢwᵣ), slg) +end + +function precipitation_tendency!( + Yₜ, + Y, + p, + t, + precip_model::Microphysics2Moment, + _, +) + (; turbconv_model) = p.atmos + (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSNᵣᵖ, ᶜSnₗᵖ, ᶜSeₜᵖ) = p.precipitation + + # Populate the cache and precipitation surface fluxes + compute_precipitation_cache!(Y, p, precip_model, turbconv_model) + compute_precipitation_surface_fluxes!(Y, p, precip_model) + + # Update grid mean tendencies + @. Yₜ.c.ρ += Y.c.ρ * ᶜSqₜᵖ + @. Yₜ.c.ρq_tot += Y.c.ρ * ᶜSqₜᵖ + @. Yₜ.c.ρe_tot += Y.c.ρ * ᶜSeₜᵖ + @. Yₜ.c.ρq_rai += Y.c.ρ * ᶜSqᵣᵖ + @. Yₜ.c.ρn_rai += ᶜSNᵣᵖ + @. Yₜ.c.ρn_liq += Y.c.ρ * ᶜSnₗᵖ + + return nothing +end diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl index b3bff9be974..756820aaf06 100644 --- a/src/parameters/create_parameters.jl +++ b/src/parameters/create_parameters.jl @@ -97,6 +97,12 @@ function create_parameter_set(config::AtmosConfig) tv = CM.Parameters.Blk1MVelType(toml_dict), aps = CM.Parameters.AirProperties(toml_dict), ) + elseif precip_model == "2M" + (; + SB2006 = CM.Parameters.SB2006(toml_dict), + SB2006Vel = CM.Parameters.SB2006VelType(toml_dict), + aps = CM.Parameters.AirProperties(toml_dict), + ) else error("Invalid precip_model $(precip_model)") end diff --git a/src/prognostic_equations/hyperdiffusion.jl b/src/prognostic_equations/hyperdiffusion.jl index cf046c05d9c..d4e265b756c 100644 --- a/src/prognostic_equations/hyperdiffusion.jl +++ b/src/prognostic_equations/hyperdiffusion.jl @@ -245,9 +245,9 @@ NVTX.@annotate function apply_tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t) # TODO: Figure out why caching the duplicated tendencies in ᶜtemp_scalar # triggers allocations. for (ᶜρχₜ, ᶜ∇²χ, χ_name) in matching_subfields(Yₜ.c, ᶜ∇²specific_tracers) - ν₄_scalar = ifelse(χ_name in (:q_rai, :q_sno), 0 * ν₄_scalar, ν₄_scalar) + ν₄_scalar = ifelse(χ_name in (:q_rai, :q_sno, :n_rai), 0 * ν₄_scalar, ν₄_scalar) @. ᶜρχₜ -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²χ)) - if !(χ_name in (:q_rai, :q_sno)) + if !(χ_name in (:q_rai, :q_sno, :n_rai)) @. Yₜ.c.ρ -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²χ)) end end diff --git a/src/prognostic_equations/implicit/implicit_solver.jl b/src/prognostic_equations/implicit/implicit_solver.jl index d8f721e9679..277b4f9da81 100644 --- a/src/prognostic_equations/implicit/implicit_solver.jl +++ b/src/prognostic_equations/implicit/implicit_solver.jl @@ -10,7 +10,7 @@ use_derivative(::IgnoreDerivative) = false """ ImplicitEquationJacobian( - Y, atmos; + Y, atmos; approximate_solve_iters, diffusion_flag, topography_flag, sgs_advection_flag, transform_flag ) @@ -161,6 +161,8 @@ function ImplicitEquationJacobian( @name(c.ρq_ice), @name(c.ρq_rai), @name(c.ρq_sno), + @name(c.ρn_liq), + @name(c.ρn_rai), ) available_tracer_names = MatrixFields.unrolled_filter(is_in_Y, tracer_names) @@ -618,6 +620,8 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ) (@name(c.ρq_ice), @name(q_ice)), (@name(c.ρq_rai), @name(q_rai)), (@name(c.ρq_sno), @name(q_sno)), + (@name(c.ρn_liq), @name(n_liq)), + (@name(c.ρn_rai), @name(n_rai)), ) MatrixFields.unrolled_foreach(tracer_info) do (ρq_name, q_name) MatrixFields.has_field(Y, ρq_name) || return @@ -676,7 +680,10 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ) ᶠlg = Fields.local_geometry_field(Y.f) precip_info = - ((@name(c.ρq_rai), @name(ᶜwᵣ)), (@name(c.ρq_sno), @name(ᶜwₛ))) + ((@name(c.ρq_rai), @name(ᶜwᵣ)), + (@name(c.ρq_sno), @name(ᶜwₛ)), + (@name(c.ρn_rai), @name(ᶜw_nᵣ)), + ) 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] diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl index 645ec5d079f..4589218f686 100644 --- a/src/prognostic_equations/implicit/implicit_tendency.jl +++ b/src/prognostic_equations/implicit/implicit_tendency.jl @@ -161,6 +161,17 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t) ᶠright_bias(Geometry.WVector(-(ᶜwₛ)) * ᶜspecific.q_sno), ) end + if precip_model isa Microphysics2Moment + (; ᶜwᵣ, ᶜw_nᵣ) = p.precomputed + @. Yₜ.c.ρq_rai -= ᶜprecipdivᵥ( + ᶠwinterp(ᶜJ, Y.c.ρ) * + ᶠright_bias(Geometry.WVector(-(ᶜwᵣ)) * ᶜspecific.q_rai), + ) + @. Yₜ.c.ρn_rai -= ᶜprecipdivᵥ( + ᶠwinterp(ᶜJ, Y.c.ρ) * + ᶠright_bias(Geometry.WVector(-(ᶜw_nᵣ)) * ᶜspecific.n_rai), + ) + end @. Yₜ.f.u₃ -= ᶠgradᵥ(ᶜp) / ᶠinterp(Y.c.ρ) + ᶠgradᵥ_ᶜΦ diff --git a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl index c9bdf3a5814..fb089f10ae8 100644 --- a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl +++ b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl @@ -67,7 +67,7 @@ function vertical_diffusion_boundary_layer_tendency!( @. ᶜρχₜ_diffusion = ᶜdivᵥ_ρχ(-(ᶠinterp(Y.c.ρ) * ᶠinterp(ᶜK_h) * ᶠgradᵥ(ᶜχ))) @. ᶜρχₜ -= ᶜρχₜ_diffusion - if !(χ_name in (:q_rai, :q_sno)) + if !(χ_name in (:q_rai, :q_sno, :n_rai)) @. Yₜ.c.ρ -= ᶜρχₜ_diffusion end end