diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index af5b0b6e19..29ac02e5f6 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -288,6 +288,14 @@ steps: agents: slurm_mem: 20GB + - label: ":computer: aquaplanet (ρe_tot) nonequilmoist allsky radiation monin_obukhov varying insolation gravity wave (gfdl_restart) high top with 1-moment micro" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res.yml + artifact_paths: "sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res/output_active/*" + agents: + slurm_mem: 20GB + - label: ":computer: aquaplanet (ρe_tot) equilmoist allsky radiation monin_obukhov varying insolation gravity wave (raw_topo) high top zonally asymmetric" command: > julia --color=yes --project=examples examples/hybrid/driver.jl diff --git a/config/model_configs/sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res.yml b/config/model_configs/sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res.yml new file mode 100644 index 0000000000..3e7ed2a50e --- /dev/null +++ b/config/model_configs/sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res.yml @@ -0,0 +1,24 @@ +z_elem: 25 +z_max: 45000.0 +dz_bottom: 300.0 +dt: "400secs" +t_end: "18hours" +dt_save_state_to_disk: "18hours" +vert_diff: "FriersonDiffusion" +implicit_diffusion: true +approximate_linear_solve_iters: 2 +moist: "nonequil" +precip_model: "nothing" +rad: "allskywithclear" +idealized_insolation: false +rayleigh_sponge: true +non_orographic_gravity_wave: true +orographic_gravity_wave: "gfdl_restart" +surface_setup: "DefaultMoninObukhov" +prescribed_aerosols: ["CB1", "CB2"] +job_id: "sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res" +toml: [toml/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.toml] +diagnostics: + - short_name: [edt, evu] + reduction_time: average + period: "18hours" diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl index a661c0f66f..2e03d5f0f1 100644 --- a/post_processing/ci_plots.jl +++ b/post_processing/ci_plots.jl @@ -824,6 +824,7 @@ end AquaplanetPlots = Union{ Val{:mpi_sphere_aquaplanet_rhoe_equilmoist_clearsky}, + Val{:sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res}, Val{:longrun_aquaplanet_rhoe_equil_55km_nz63_gray_0M}, Val{:longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_0M}, Val{:longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_diagedmf_diffonly_0M}, diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl index 4294fe61b1..49a43c0b92 100644 --- a/src/ClimaAtmos.jl +++ b/src/ClimaAtmos.jl @@ -90,6 +90,9 @@ include(joinpath("prognostic_equations", "edmfx_precipitation.jl")) include( joinpath("parameterized_tendencies", "microphysics", "precipitation.jl"), ) +include( + joinpath("parameterized_tendencies", "microphysics", "cloud_condensate.jl"), +) include( joinpath( "parameterized_tendencies", diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl index c5f3971ccd..f64bb0d65a 100644 --- a/src/cache/diagnostic_edmf_precomputed_quantities.jl +++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl @@ -243,7 +243,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( (; ᶠu³⁰, ᶜK⁰, ᶜtke⁰) = p.precomputed thermo_params = CAP.thermodynamics_params(params) - microphys_params = CAP.microphysics_params(params) + microphys_params = CAP.microphysics_precipitation_params(params) ᶠΦ = p.scratch.ᶠtemp_scalar @. ᶠΦ = CAP.grav(params) * ᶠz @@ -874,7 +874,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita precip_model::Microphysics0Moment, ) thermo_params = CAP.thermodynamics_params(p.params) - microphys_params = CAP.microphysics_params(p.params) + microphys_params = CAP.microphysics_precipitation_params(p.params) (; dt) = p (; ᶜts, ᶜSqₜᵖ⁰) = p.precomputed (; q_tot) = p.precomputed.ᶜspecific diff --git a/src/cache/precipitation_precomputed_quantities.jl b/src/cache/precipitation_precomputed_quantities.jl index 21c0fc7f23..f4b25dbbe4 100644 --- a/src/cache/precipitation_precomputed_quantities.jl +++ b/src/cache/precipitation_precomputed_quantities.jl @@ -19,7 +19,7 @@ function set_precipitation_precomputed_quantities!(Y, p, t) (; ᶜwᵣ, ᶜwₛ, ᶜqᵣ, ᶜqₛ) = p.precomputed - cmp = CAP.microphysics_params(p.params) + cmp = CAP.microphysics_precipitation_params(p.params) # compute the precipitation specific humidities @. ᶜqᵣ = qₚ(Y.c.ρq_rai, Y.c.ρ) diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl index 8542205849..c63e764cd8 100644 --- a/src/cache/prognostic_edmf_precomputed_quantities.jl +++ b/src/cache/prognostic_edmf_precomputed_quantities.jl @@ -363,7 +363,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation (; params, dt) = p thp = CAP.thermodynamics_params(params) - cmp = CAP.microphysics_params(params) + cmp = CAP.microphysics_precipitation_params(params) (; ᶜts⁰, ᶜq_tot⁰, ᶜtsʲs, ᶜSqₜᵖʲs, ᶜSqₜᵖ⁰) = p.precomputed # Sources from the updrafts @@ -399,7 +399,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation (; params, dt) = p (; ᶜΦ,) = p.core thp = CAP.thermodynamics_params(params) - cmp = CAP.microphysics_params(params) + cmp = CAP.microphysics_precipitation_params(params) (; ᶜSeₜᵖʲs, ᶜSqₜᵖʲs, ᶜSqᵣᵖʲs, ᶜSqₛᵖʲs, ᶜρʲs, ᶜtsʲs) = p.precomputed (; ᶜSeₜᵖ⁰, ᶜSqₜᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰, ᶜρ⁰, ᶜts⁰) = p.precomputed diff --git a/src/parameterized_tendencies/microphysics/cloud_condensate.jl b/src/parameterized_tendencies/microphysics/cloud_condensate.jl new file mode 100644 index 0000000000..a8d9e6a40e --- /dev/null +++ b/src/parameterized_tendencies/microphysics/cloud_condensate.jl @@ -0,0 +1,20 @@ +##### +##### DryModel, EquilMoistModel +##### + +cloud_condensate_tendency!(Yₜ, p, colidx, _) = nothing + +##### +##### NonEquilMoistModel +##### + +function cloud_condensate_tendency!(Yₜ, p, colidx, ::NonEquilMoistModel) + + (; ᶜts) = p.precomputed + (; params, dt) = p + thp = CAP.thermodynamics_params(params) + cmc = CAP.microphysics_cloud_params(params) + + @. Yₜ.c.ρq_liq[colidx] += cloud_sources(cmc.liquid, thp, ᶜts[colidx], dt) + @. Yₜ.c.ρq_ice[colidx] += cloud_sources(cmc.ice, thp, ᶜts[colidx], dt) +end diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl index ca1559f243..e543bdaa0b 100644 --- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl +++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl @@ -2,6 +2,9 @@ import Thermodynamics as TD import CloudMicrophysics.Microphysics0M as CM0 +import CloudMicrophysics.Microphysics1M as CM1 +import CloudMicrophysics.MicrophysicsNonEq as CMNe +import CloudMicrophysics.Parameters as CMP # define some aliases and functions to make the code more readable const Iₗ = TD.internal_energy_liquid @@ -10,12 +13,64 @@ const Lf = TD.latent_heat_fusion const Tₐ = TD.air_temperature const PP = TD.PhasePartition const qᵥ = TD.vapor_specific_humidity +qₜ(thp, ts) = TD.PhasePartition(thp, ts).tot qₗ(thp, ts) = TD.PhasePartition(thp, ts).liq qᵢ(thp, ts) = TD.PhasePartition(thp, ts).ice # helper function to limit the tendency -function limit(q::FT, dt::FT) where {FT} - return q / dt / 5 +function limit(q::FT, dt::FT, n::Int) where {FT} + return q / dt / n +end + +""" + cloud_sources(cm_params, thp, ts, dt) + + - cm_params - CloudMicrophysics parameters struct for cloud water or ice condensate + - thp - Thermodynamics parameters struct + - ts - thermodynamics state + - dt - model time step + +Returns the condensation/evaporation or deposition/sublimation rate for +non-equilibrium cloud formation. +""" +function cloud_sources(cm_params::CMP.CloudLiquid{FT}, thp, ts, dt) where {FT} + + λ = TD.liquid_fraction(thp, Tₐ(thp, ts), TD.PhaseEquil) + qᵥ_ex_liquid = λ * (qₜ(thp, ts) - TD.q_vap_saturation_liquid(thp, ts)) + + # Ideally the logic whether to apply this source term, and what relaxation + # timescale to use, should be based on the availability of CCNs and INPs. + # We need a 2-moment microphysics scheme for that. + + # Additionally, to approximate the partitioning between the condensation/evaporation + # and deposition/sublimation, we are scaling down the excess q by the liquid fraction. + # This needs more thought. + S = CMNe.conv_q_vap_to_q_liq_ice( + cm_params, + PP(FT(0), qᵥ_ex_liquid, FT(0)), + PP(thp, ts), + ) + return ifelse( + S > FT(0), + min(S, limit(qᵥ(thp, ts), dt, 2)), + -min(abs(S), limit(qₗ(thp, ts), dt, 2)), + ) +end +function cloud_sources(cm_params::CMP.CloudIce{FT}, thp, ts, dt) where {FT} + + λ = TD.liquid_fraction(thp, Tₐ(thp, ts), TD.PhaseEquil) + qᵥ_ex_ice = (1 - λ) * (qₜ(thp, ts) - TD.q_vap_saturation_ice(thp, ts)) + + S = CMNe.conv_q_vap_to_q_liq_ice( + cm_params, + PP(FT(0), FT(0), qᵥ_ex_ice), + PP(thp, ts), + ) + return ifelse( + S > FT(0), + min(S, limit(qᵥ(thp, ts), dt, 2)), + -min(abs(S), limit(qᵢ(thp, ts), dt, 2)), + ) end """ @@ -114,7 +169,7 @@ function compute_precipitation_sources!( #! format: off # rain autoconversion: q_liq -> q_rain @. Sᵖ = min( - limit(qₗ(thp, ts), dt), + limit(qₗ(thp, ts), dt, 5), CM1.conv_q_liq_to_q_rai(mp.pr.acnv1M, qₗ(thp, ts), true), ) @. Sqₜᵖ -= Sᵖ @@ -123,7 +178,7 @@ function compute_precipitation_sources!( # snow autoconversion assuming no supersaturation: q_ice -> q_snow @. Sᵖ = min( - limit(qᵢ(thp, ts), dt), + limit(qᵢ(thp, ts), dt, 5), CM1.conv_q_ice_to_q_sno_no_supersat(mp.ps.acnv1M, qᵢ(thp, ts), true), ) @. Sqₜᵖ -= Sᵖ @@ -132,7 +187,7 @@ function compute_precipitation_sources!( # accretion: q_liq + q_rain -> q_rain @. Sᵖ = min( - limit(qₗ(thp, ts), dt), + limit(qₗ(thp, ts), dt, 5), CM1.accretion(mp.cl, mp.pr, mp.tv.rain, mp.ce, qₗ(thp, ts), qᵣ, ρ), ) @. Sqₜᵖ -= Sᵖ @@ -141,7 +196,7 @@ function compute_precipitation_sources!( # accretion: q_ice + q_snow -> q_snow @. Sᵖ = min( - limit(qᵢ(thp, ts), dt), + limit(qᵢ(thp, ts), dt, 5), CM1.accretion(mp.ci, mp.ps, mp.tv.snow, mp.ce, qᵢ(thp, ts), qₛ, ρ), ) @. Sqₜᵖ -= Sᵖ @@ -151,7 +206,7 @@ function compute_precipitation_sources!( # accretion: q_liq + q_sno -> q_sno or q_rai # sink of cloud water via accretion cloud water + snow @. Sᵖ = min( - limit(qₗ(thp, ts), dt), + limit(qₗ(thp, ts), dt, 5), CM1.accretion(mp.cl, mp.ps, mp.tv.snow, mp.ce, qₗ(thp, ts), qₛ, ρ), ) # if T < T_freeze cloud droplets freeze to become snow @@ -160,7 +215,7 @@ function compute_precipitation_sources!( @. Sᵖ_snow = ifelse( Tₐ(thp, ts) < mp.ps.T_freeze, Sᵖ, - FT(-1) * min(Sᵖ * α(thp, ts), limit(qₛ, dt)), + FT(-1) * min(Sᵖ * α(thp, ts), limit(qₛ, dt, 5)), ) @. Sqₛᵖ += Sᵖ_snow @. Sqₜᵖ -= Sᵖ @@ -173,7 +228,7 @@ function compute_precipitation_sources!( # accretion: q_ice + q_rai -> q_sno @. Sᵖ = min( - limit(qᵢ(thp, ts), dt), + limit(qᵢ(thp, ts), dt, 5), CM1.accretion(mp.ci, mp.pr, mp.tv.rain, mp.ce, qᵢ(thp, ts), qᵣ, ρ), ) @. Sqₜᵖ -= Sᵖ @@ -181,7 +236,7 @@ function compute_precipitation_sources!( @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ) # sink of rain via accretion cloud ice - rain @. Sᵖ = min( - limit(qᵣ, dt), + limit(qᵣ, dt, 5), CM1.accretion_rain_sink(mp.pr, mp.ci, mp.tv.rain, mp.ce, qᵢ(thp, ts), qᵣ, ρ), ) @. Sqᵣᵖ -= Sᵖ @@ -192,11 +247,11 @@ function compute_precipitation_sources!( @. Sᵖ = ifelse( Tₐ(thp, ts) < mp.ps.T_freeze, min( - limit(qᵣ, dt), + limit(qᵣ, dt, 5), CM1.accretion_snow_rain(mp.ps, mp.pr, mp.tv.rain, mp.tv.snow, mp.ce, qₛ, qᵣ, ρ), ), -min( - limit(qₛ, dt), + limit(qₛ, dt, 5), CM1.accretion_snow_rain(mp.pr, mp.ps, mp.tv.snow, mp.tv.rain, mp.ce, qᵣ, qₛ, ρ), ), ) @@ -245,7 +300,7 @@ function compute_precipitation_sinks!( #! format: off # evaporation: q_rai -> q_vap @. Sᵖ = -min( - limit(qᵣ, dt), + limit(qᵣ, dt, 5), -CM1.evaporation_sublimation(rps..., PP(thp, ts), qᵣ, ρ, Tₐ(thp, ts)), ) @. Sqₜᵖ -= Sᵖ @@ -254,7 +309,7 @@ function compute_precipitation_sinks!( # melting: q_sno -> q_rai @. Sᵖ = min( - limit(qₛ, dt), + limit(qₛ, dt, 5), CM1.snow_melt(sps..., qₛ, ρ, Tₐ(thp, ts)), ) @. Sqᵣᵖ += Sᵖ @@ -265,8 +320,8 @@ function compute_precipitation_sinks!( @. Sᵖ = CM1.evaporation_sublimation(sps..., PP(thp, ts), qₛ, ρ, Tₐ(thp, ts)) @. Sᵖ = ifelse( Sᵖ > FT(0), - min(limit(qᵥ(thp, ts), dt), Sᵖ), - -min(limit(qₛ, dt), FT(-1) * Sᵖ), + min(limit(qᵥ(thp, ts), dt, 5), Sᵖ), + -min(limit(qₛ, dt, 5), FT(-1) * Sᵖ), ) @. Sqₜᵖ -= Sᵖ @. Sqₛᵖ += Sᵖ diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl index 490220b44d..9abadeb21a 100644 --- a/src/parameterized_tendencies/microphysics/precipitation.jl +++ b/src/parameterized_tendencies/microphysics/precipitation.jl @@ -45,7 +45,7 @@ function compute_precipitation_cache!(Y, p, colidx, ::Microphysics0Moment, _) (; ᶜts) = p.precomputed (; ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precipitation (; ᶜΦ) = p.core - cm_params = CAP.microphysics_params(params) + cm_params = CAP.microphysics_precipitation_params(params) thermo_params = CAP.thermodynamics_params(params) @. ᶜS_ρq_tot[colidx] = Y.c.ρ[colidx] * q_tot_precipitation_sources( @@ -215,7 +215,7 @@ function compute_precipitation_cache!(Y, p, colidx, ::Microphysics1Moment, _) # get thermodynamics and 1-moment microphysics params (; params) = p - cmp = CAP.microphysics_params(params) + cmp = CAP.microphysics_precipitation_params(params) thp = CAP.thermodynamics_params(params) # compute precipitation source terms on the grid mean @@ -273,7 +273,7 @@ function compute_precipitation_cache!( # get thermodynamics and 1-moment microphysics params (; params) = p - cmp = CAP.microphysics_params(params) + cmp = CAP.microphysics_precipitation_params(params) thp = CAP.thermodynamics_params(params) # zero out the helper source terms diff --git a/src/parameters/Parameters.jl b/src/parameters/Parameters.jl index b403c635f4..f32c8b858c 100644 --- a/src/parameters/Parameters.jl +++ b/src/parameters/Parameters.jl @@ -44,12 +44,22 @@ Base.@kwdef struct TurbulenceConvectionParameters{FT} <: ATCP max_area_limiter_power::FT end -Base.@kwdef struct ClimaAtmosParameters{FT, TP, RP, IP, MPP, WP, SFP, TCP} <: - ACAP +Base.@kwdef struct ClimaAtmosParameters{ + FT, + TP, + RP, + IP, + MPC, + MPP, + WP, + SFP, + TCP, +} <: ACAP thermodynamics_params::TP rrtmgp_params::RP insolation_params::IP - microphysics_params::MPP + microphysics_cloud_params::MPC + microphysics_precipitation_params::MPP water_params::WP surface_fluxes_params::SFP turbconv_params::TCP diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl index 73ae5b1adc..b3bff9be97 100644 --- a/src/parameters/create_parameters.jl +++ b/src/parameters/create_parameters.jl @@ -64,6 +64,17 @@ function create_parameter_set(config::AtmosConfig) SF.Parameters.SurfaceFluxesParameters(toml_dict, UF.BusingerParams) SFP = typeof(surface_fluxes_params) + moisture_model = parsed_args["moist"] + microphysics_cloud_params = if moisture_model == "nonequil" + (; + liquid = CM.Parameters.CloudLiquid(toml_dict), + ice = CM.Parameters.CloudIce(toml_dict), + ) + else + nothing + end + MPC = typeof(microphysics_cloud_params) + # Microphysics scheme parameters (from CloudMicrophysics.jl) # TODO - repeating the logic from solver/model_getters.jl... if parsed_args["override_τ_precip"] @@ -71,7 +82,7 @@ function create_parameter_set(config::AtmosConfig) FT(CA.time_to_seconds(parsed_args["dt"])) end precip_model = parsed_args["precip_model"] - microphysics_params = + microphysics_precipitation_params = if precip_model == nothing || precip_model == "nothing" nothing elseif precip_model == "0M" @@ -89,7 +100,7 @@ function create_parameter_set(config::AtmosConfig) else error("Invalid precip_model $(precip_model)") end - MPP = typeof(microphysics_params) + MPP = typeof(microphysics_precipitation_params) name_map = (; :f_plane_coriolis_frequency => :f_plane_coriolis_frequency, @@ -115,12 +126,13 @@ function create_parameter_set(config::AtmosConfig) :water_refractive_index => :water_refractive_index, ) parameters = CP.get_parameter_values(toml_dict, name_map, "ClimaAtmos") - return CAP.ClimaAtmosParameters{FT, TP, RP, IP, MPP, WP, SFP, TCP}(; + return CAP.ClimaAtmosParameters{FT, TP, RP, IP, MPC, MPP, WP, SFP, TCP}(; parameters..., thermodynamics_params, rrtmgp_params, insolation_params, - microphysics_params, + microphysics_cloud_params, + microphysics_precipitation_params, water_params, surface_fluxes_params, turbconv_params, diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl index 5ba03a7872..c5f6c49e88 100644 --- a/src/prognostic_equations/remaining_tendency.jl +++ b/src/prognostic_equations/remaining_tendency.jl @@ -84,7 +84,8 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t) p.atmos.turbconv_model, ) edmfx_tke_tendency!(Yₜ, Y, p, t, colidx, p.atmos.turbconv_model) - # TODO - add the 1-moment precipitation microphysics here + # Non-equilibrium cloud formation + cloud_condensate_tendency!(Yₜ, p, colidx, p.atmos.moisture_model) edmfx_precipitation_tendency!( Yₜ, Y, diff --git a/test/parameterized_tendencies/microphysics/precipitation.jl b/test/parameterized_tendencies/microphysics/precipitation.jl index 2c50da307b..9dcee39c7c 100644 --- a/test/parameterized_tendencies/microphysics/precipitation.jl +++ b/test/parameterized_tendencies/microphysics/precipitation.jl @@ -55,8 +55,8 @@ include("../../test_helpers.jl") @test maximum(abs.(p.precipitation.ᶜS_ρq_tot)) <= sqrt(eps(FT)) # Test that tendencies result in correct water-mass budget, - # and that the tendency modification corresponds exactly to the - # cached source term. + # and that the tendency modification corresponds exactly to the + # cached source term. CC.Fields.bycolumn(axes(Y.c)) do colidx CA.precipitation_tendency!( ᶜYₜ, @@ -93,7 +93,7 @@ include("../../test_helpers.jl") end @test CA.qₚ(FT(10), FT(2)) == FT(5) @test CA.qₚ(FT(-10), FT(2)) == FT(0) - @test CA.limit(FT(10), FT(2)) == FT(1) + @test CA.limit(FT(10), FT(2), 5) == FT(1) CC.Fields.bycolumn(axes(Y.c)) do colidx CA.precipitation_tendency!( ᶜYₜ,