diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 2f4b200b4b..a62cff029f 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -322,16 +322,16 @@ steps: - group: "Sphere Examples (Aquaplanet)" steps: - - label: ":umbrella: aquaplanet (ρe_tot) equil allsky monin_obukhov varying insol gravity wave (gfdl_restart) high top 1-moment" + - label: ":umbrella: aquaplanet (ρe_tot) nonequil allsky monin_obukhov varying insol gravity wave (gfdl_restart) high top 1-moment" command: > julia --color=yes --project=examples examples/hybrid/driver.jl - --config_file $CONFIG_PATH/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml - --job_id sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res + --config_file $CONFIG_PATH/sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res.yml + --job_id sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res julia --color=yes --project=examples reproducibility_tests/test_mse.jl - --job_id sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res - --out_dir sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res/output_active - artifact_paths: "sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res/output_active/*" + --job_id sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res + --out_dir sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res/output_active + artifact_paths: "sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res/output_active/*" agents: slurm_mem: 20GB @@ -349,7 +349,7 @@ steps: slurm_mem: 20GB slurm_constraint: icelake|cascadelake|skylake|epyc - - label: ":computer: aquaplanet (ρe_tot) nonequilmoist allsky radiation monin_obukhov varying insolation high top 1-moment" + - label: ":umbrella: aquaplanet (ρe_tot) nonequilmoist allsky radiation monin_obukhov varying insolation high top 1-moment" command: > julia --color=yes --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/sphere_aquaplanet_rhoe_nonequilmoist_allsky.yml diff --git a/Project.toml b/Project.toml index 09cc5802cd..1bb1396387 100644 --- a/Project.toml +++ b/Project.toml @@ -45,7 +45,7 @@ ClimaComms = "0.6.4" ClimaCore = "0.14.19" ClimaDiagnostics = "0.2.4" ClimaParams = "0.10.17" -ClimaTimeSteppers = "0.7.39" +ClimaTimeSteppers = "0.7, 0.8" ClimaUtilities = "0.1.14" CloudMicrophysics = "0.22.3" Dates = "1" diff --git a/config/longrun_configs/longrun_moist_baroclinic_wave.yml b/config/longrun_configs/longrun_moist_baroclinic_wave.yml index 2033944d58..2aa405122f 100644 --- a/config/longrun_configs/longrun_moist_baroclinic_wave.yml +++ b/config/longrun_configs/longrun_moist_baroclinic_wave.yml @@ -1,14 +1,17 @@ -h_elem: 30 -z_elem: 43 -dz_bottom: 30.0 -dt: "90secs" -t_end: "120days" +h_elem: 6 #30 +z_elem: 10 #43 +dz_bottom: 300.0 +dt: "200secs" +approximate_linear_solve_iters: 3 +max_newton_iters_ode: 1 +rayleigh_sponge: true +t_end: "10days" ode_algo: ARS343 initial_condition: "MoistBaroclinicWave" -moist: "equil" -precip_model: "0M" +moist: "nonequil" +precip_model: "1M" dt_save_state_to_disk: "10days" toml: [toml/longrun_baroclinic_wave.toml] diagnostics: - - short_name: [pfull, wa, va, rv, hus, ke] - period: 1days + - short_name: [pfull, wa, va, ua, ta, rhoa, rv, hus, ke, clw, cli, husra, hussn, ] + period: 600secs diff --git a/config/model_configs/diagnostic_edmfx_aquaplanet.yml b/config/model_configs/diagnostic_edmfx_aquaplanet.yml index bcef701417..baa41b5bf2 100644 --- a/config/model_configs/diagnostic_edmfx_aquaplanet.yml +++ b/config/model_configs/diagnostic_edmfx_aquaplanet.yml @@ -17,7 +17,7 @@ edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true moist: equil cloud_model: "quadrature_sgs" -precip_model: 1M +precip_model: 0M dt: 120secs t_end: 3hours reproducibility_test: true diff --git a/config/model_configs/diagnostic_edmfx_dycoms_rf02_box.yml b/config/model_configs/diagnostic_edmfx_dycoms_rf02_box.yml index e613149bd6..520c20528d 100644 --- a/config/model_configs/diagnostic_edmfx_dycoms_rf02_box.yml +++ b/config/model_configs/diagnostic_edmfx_dycoms_rf02_box.yml @@ -15,7 +15,7 @@ edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true moist: equil cloud_model: "quadrature_sgs" -precip_model: "1M" +precip_model: "0M" call_cloud_diagnostics_per_stage: true config: box x_max: 1e8 @@ -35,6 +35,4 @@ diagnostics: period: 10mins - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix] period: 10mins - - short_name: [husra, hussn] - period: 10mins ode_algo: ARS343 diff --git a/config/model_configs/diagnostic_edmfx_rico_box.yml b/config/model_configs/diagnostic_edmfx_rico_box.yml index 3efdffd675..c93def193a 100644 --- a/config/model_configs/diagnostic_edmfx_rico_box.yml +++ b/config/model_configs/diagnostic_edmfx_rico_box.yml @@ -14,7 +14,7 @@ edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true moist: equil cloud_model: "quadrature_sgs" -precip_model: "1M" +precip_model: "0M" call_cloud_diagnostics_per_stage: true config: box x_max: 1e8 @@ -34,6 +34,4 @@ diagnostics: period: 10mins - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix] period: 10mins - - short_name: [husra, hussn] - period: 10mins ode_algo: ARS343 diff --git a/config/model_configs/diagnostic_edmfx_trmm_box.yml b/config/model_configs/diagnostic_edmfx_trmm_box.yml index 54369f1321..dbe6663bb9 100644 --- a/config/model_configs/diagnostic_edmfx_trmm_box.yml +++ b/config/model_configs/diagnostic_edmfx_trmm_box.yml @@ -15,7 +15,7 @@ moist: equil cloud_model: "quadrature_sgs" call_cloud_diagnostics_per_stage: true apply_limiter: false -precip_model: "1M" +precip_model: "0M" config: box x_max: 1e8 y_max: 1e8 @@ -35,6 +35,4 @@ diagnostics: period: 10mins - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix] period: 10mins - - short_name: [husra, hussn] - period: 10mins ode_algo: ARS343 diff --git a/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml b/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml index fa098bae39..161eb7958e 100644 --- a/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml +++ b/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml @@ -15,7 +15,7 @@ edmfx_sgs_diffusive_flux: true moist: equil cloud_model: "quadrature_sgs" call_cloud_diagnostics_per_stage: true -precip_model: "1M" +precip_model: "0M" config: box x_max: 1e8 y_max: 1e8 @@ -35,6 +35,4 @@ diagnostics: period: 10mins - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix] period: 10mins - - short_name: [husra, hussn] - period: 10mins ode_algo: ARS343 diff --git a/config/model_configs/sphere_aquaplanet_rhoe_nonequilmoist_allsky.yml b/config/model_configs/sphere_aquaplanet_rhoe_nonequilmoist_allsky.yml index 26a58f0a51..725f2ec5ac 100644 --- a/config/model_configs/sphere_aquaplanet_rhoe_nonequilmoist_allsky.yml +++ b/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 diff --git a/config/model_configs/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml b/config/model_configs/sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res.yml similarity index 97% rename from config/model_configs/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml rename to config/model_configs/sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res.yml index 2a7fc1e5bc..f164ac4a24 100644 --- a/config/model_configs/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml +++ b/config/model_configs/sphere_aquaplanet_rhoe_nonequilmoist_allsky_gw_res.yml @@ -8,7 +8,7 @@ dt_save_state_to_disk: "24hours" vert_diff: "FriersonDiffusion" implicit_diffusion: true approximate_linear_solve_iters: 2 -moist: "equil" +moist: "nonequil" precip_model: "1M" rad: "allskywithclear" aerosol_radiation: true diff --git a/examples/Manifest-v1.11.toml b/examples/Manifest-v1.11.toml index 7a20391539..0eefa59643 100644 --- a/examples/Manifest-v1.11.toml +++ b/examples/Manifest-v1.11.toml @@ -359,7 +359,9 @@ weakdeps = ["CUDA", "MPI"] [[deps.ClimaCore]] deps = ["Adapt", "BandedMatrices", "BlockArrays", "ClimaComms", "CubedSphere", "DataStructures", "ForwardDiff", "GaussQuadrature", "GilbertCurves", "HDF5", "InteractiveUtils", "IntervalSets", "KrylovKit", "LinearAlgebra", "MultiBroadcastFusion", "NVTX", "PkgVersion", "RecursiveArrayTools", "RootSolvers", "SparseArrays", "StaticArrays", "Statistics", "Unrolled"] -git-tree-sha1 = "61d071cfe584d99a1400a3bf96a122c45c282121" +git-tree-sha1 = "ce9de8f8c74e523a560b66b067c220438b287974" +repo-rev = "main" +repo-url = "https://github.com/CliMA/ClimaCore.jl.git" uuid = "d414da3d-4745-48bb-8d80-42e94e092884" version = "0.14.24" weakdeps = ["CUDA", "Krylov"] @@ -418,9 +420,11 @@ version = "0.1.1" [[deps.ClimaTimeSteppers]] deps = ["ClimaComms", "Colors", "DataStructures", "DiffEqBase", "KernelAbstractions", "Krylov", "LinearAlgebra", "LinearOperators", "NVTX", "SciMLBase", "StaticArrays"] -git-tree-sha1 = "ff967e27a56a938d3e719bb8d2e025b086c3f808" +git-tree-sha1 = "facfe93769c3343bbe512fb61b19cfdf85804118" +repo-rev = "main" +repo-url = "https://github.com/CliMA/ClimaTimeSteppers.jl.git" uuid = "595c0a79-7f3d-439a-bc5a-b232dc3bde79" -version = "0.7.39" +version = "0.8.1" [deps.ClimaTimeSteppers.extensions] ClimaTimeSteppersBenchmarkToolsExt = ["CUDA", "BenchmarkTools", "OrderedCollections", "StatsBase", "PrettyTables"] diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl index 75d1c4d143..6c83828181 100644 --- a/post_processing/ci_plots.jl +++ b/post_processing/ci_plots.jl @@ -1190,17 +1190,16 @@ EDMFBoxPlots = Union{ Val{:prognostic_edmfx_bomex_box}, Val{:rcemipii_box_diagnostic_edmfx}, Val{:prognostic_edmfx_soares_column}, -} - -EDMFBoxPlotsWithPrecip = Union{ Val{:diagnostic_edmfx_dycoms_rf02_box}, - Val{:prognostic_edmfx_rico_column}, - Val{:prognostic_edmfx_trmm_column}, Val{:diagnostic_edmfx_rico_box}, Val{:diagnostic_edmfx_trmm_box}, Val{:diagnostic_edmfx_trmm_stretched_box}, } +EDMFBoxPlotsWithPrecip = Union{ + Val{:prognostic_edmfx_rico_column}, + Val{:prognostic_edmfx_trmm_column}, +} """ plot_edmf_vert_profile!(grid_loc, var_group) diff --git a/src/cache/cloud_fraction.jl b/src/cache/cloud_fraction.jl index f05f78d39f..6fd4ac0a71 100644 --- a/src/cache/cloud_fraction.jl +++ b/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 @@ -24,10 +31,11 @@ NVTX.@annotate function set_cloud_fraction!(Y, p, ::DryModel, _) p.precomputed.cloud_diagnostics_tuple .= ((; cf = FT(0), q_liq = FT(0), q_ice = FT(0)),) end + NVTX.@annotate function set_cloud_fraction!( Y, p, - ::Union{EquilMoistModel, NonEquilMoistModel}, + moist_model::Union{EquilMoistModel, NonEquilMoistModel}, ::GridScaleCloud, ) (; params) = p @@ -48,13 +56,24 @@ 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, @@ -64,7 +83,7 @@ NVTX.@annotate function set_cloud_fraction!( ) SG_quad = qc.SG_quad (; params) = p - + #TODO - do we want to substract precipitation for the noneq option? FT = eltype(params) thermo_params = CAP.thermodynamics_params(params) (; ᶜts, ᶜmixing_length, cloud_diagnostics_tuple) = p.precomputed @@ -103,6 +122,7 @@ NVTX.@annotate function set_cloud_fraction!( moisture_model::Union{EquilMoistModel, NonEquilMoistModel}, cloud_model::SGSQuadratureCloud, ) + #TODO - do we want to substract precipitation for the noneq option? (; turbconv_model) = p.atmos set_cloud_fraction!(Y, p, moisture_model, cloud_model, turbconv_model) @@ -119,6 +139,8 @@ NVTX.@annotate function set_cloud_fraction!( SG_quad = qc.SG_quad (; params) = p + #TODO - do we want to substract precipitation for the noneq option? + FT = eltype(params) thermo_params = CAP.thermodynamics_params(params) (; ᶜts, ᶜmixing_length, cloud_diagnostics_tuple) = p.precomputed @@ -143,21 +165,18 @@ NVTX.@annotate function set_cloud_fraction!( n = n_mass_flux_subdomains(turbconv_model) for j in 1:n - @. cloud_diagnostics_tuple += NamedTuple{(:cf, :q_liq, :q_ice)}( - tuple( - ifelse( - TD.has_condensate(thermo_params, ᶜtsʲs.:($$j)), - draft_area(ᶜρaʲs.:($$j), ᶜρʲs.:($$j)), - 0, - ), - draft_area(ᶜρaʲs.:($$j), ᶜρʲs.:($$j)) * - TD.PhasePartition(thermo_params, ᶜtsʲs.:($$j)).liq, - draft_area(ᶜρaʲs.:($$j), ᶜρʲs.:($$j)) * - TD.PhasePartition(thermo_params, ᶜtsʲs.:($$j)).ice, + @. cloud_diagnostics_tuple += make_named_tuple( + ifelse( + TD.has_condensate(thermo_params, ᶜtsʲs.:($$j)), + draft_area(ᶜρaʲs.:($$j), ᶜρʲs.:($$j)), + 0, ), + draft_area(ᶜρaʲs.:($$j), ᶜρʲs.:($$j)) * + TD.PhasePartition(thermo_params, ᶜtsʲs.:($$j)).liq, + draft_area(ᶜρaʲs.:($$j), ᶜρʲs.:($$j)) * + TD.PhasePartition(thermo_params, ᶜtsʲs.:($$j)).ice, ) end - end NVTX.@annotate function set_cloud_fraction!( @@ -176,6 +195,7 @@ NVTX.@annotate function set_cloud_fraction!( (; ᶜρʲs, ᶜtsʲs, ᶜρa⁰, ᶜρ⁰) = p.precomputed (; turbconv_model) = p.atmos + # TODO - do we want to substract precipitation for the noneq option? # TODO - we should make this default when using diagnostic edmf # environment diagnostic_covariance_coeff = CAP.diagnostic_covariance_coeff(params) @@ -191,29 +211,25 @@ NVTX.@annotate function set_cloud_fraction!( ) # weight cloud diagnostics by environmental area - @. cloud_diagnostics_tuple *= NamedTuple{(:cf, :q_liq, :q_ice)}( - tuple( - draft_area(ᶜρa⁰, ᶜρ⁰), - draft_area(ᶜρa⁰, ᶜρ⁰), - draft_area(ᶜρa⁰, ᶜρ⁰), - ), + @. cloud_diagnostics_tuple *= make_named_tuple( + draft_area(ᶜρa⁰, ᶜρ⁰), + draft_area(ᶜρa⁰, ᶜρ⁰), + draft_area(ᶜρa⁰, ᶜρ⁰), ) # updrafts n = n_mass_flux_subdomains(turbconv_model) for j in 1:n - @. cloud_diagnostics_tuple += NamedTuple{(:cf, :q_liq, :q_ice)}( - tuple( - ifelse( - TD.has_condensate(thermo_params, ᶜtsʲs.:($$j)), - draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)), - 0, - ), - draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)) * - TD.PhasePartition(thermo_params, ᶜtsʲs.:($$j)).liq, - draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)) * - TD.PhasePartition(thermo_params, ᶜtsʲs.:($$j)).ice, + @. cloud_diagnostics_tuple += make_named_tuple( + ifelse( + TD.has_condensate(thermo_params, ᶜtsʲs.:($$j)), + draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)), + 0, ), + draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)) * + TD.PhasePartition(thermo_params, ᶜtsʲs.:($$j)).liq, + draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)) * + TD.PhasePartition(thermo_params, ᶜtsʲs.:($$j)).ice, ) end end diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl index 59d88ba0c3..9631799275 100644 --- a/src/cache/diagnostic_edmf_precomputed_quantities.jl +++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl @@ -535,8 +535,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( # To be applied in updraft continuity, moisture and energy # for updrafts and grid mean if precip_model isa Microphysics0Moment - @. S_q_totʲ_prev_level = q_tot_precipitation_sources( - precip_model, + @. S_q_totʲ_prev_level = q_tot_0M_precipitation_sources( thermo_params, microphys_0m_params, dt, @@ -1034,8 +1033,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita (; q_tot) = p.precomputed.ᶜspecific # Environment precipitation sources (to be applied to grid mean) - @. ᶜSqₜᵖ⁰ = q_tot_precipitation_sources( - precip_model, + @. ᶜSqₜᵖ⁰ = q_tot_0M_precipitation_sources( thermo_params, microphys_0m_params, dt, diff --git a/src/cache/precipitation_precomputed_quantities.jl b/src/cache/precipitation_precomputed_quantities.jl index b25ff7168f..f19e888fcc 100644 --- a/src/cache/precipitation_precomputed_quantities.jl +++ b/src/cache/precipitation_precomputed_quantities.jl @@ -2,11 +2,7 @@ ##### Precomputed quantities ##### import CloudMicrophysics.Microphysics1M as CM1 - -# helper function to safely get precipitation from state -function qₚ(ρqₚ, ρ) - return max(zero(ρ), ρqₚ / ρ) -end +import CloudMicrophysics.MicrophysicsNonEq as CMNe """ set_precipitation_precomputed_quantities!(Y, p, t) @@ -17,33 +13,28 @@ for the 1-moment microphysics scheme function set_precipitation_precomputed_quantities!(Y, p, t) @assert (p.atmos.precip_model isa Microphysics1Moment) - (; ᶜwᵣ, ᶜwₛ, ᶜqᵣ, ᶜqₛ) = p.precomputed - + (; ᶜwᵣ, ᶜwₛ) = p.precomputed cmp = CAP.microphysics_1m_params(p.params) - # compute the precipitation specific humidities - @. ᶜqᵣ = qₚ(Y.c.ρq_rai, Y.c.ρ) - @. ᶜqₛ = qₚ(Y.c.ρq_sno, Y.c.ρ) - # compute the precipitation terminal velocity [m/s] @. ᶜwᵣ = CM1.terminal_velocity( cmp.pr, cmp.tv.rain, Y.c.ρ, - abs(Y.c.ρq_rai / Y.c.ρ), + max(zero(Y.c.ρ), Y.c.ρq_rai / Y.c.ρ), ) @. ᶜwₛ = CM1.terminal_velocity( cmp.ps, cmp.tv.snow, Y.c.ρ, - abs(Y.c.ρq_sno / Y.c.ρ), + max(zero(Y.c.ρ), Y.c.ρq_sno / Y.c.ρ), ) return nothing end - """ set_sedimentation_precomputed_quantities!(Y, p, t) + Updates the sedimentation terminal velocity stored in `p` for the non-equilibrium microphysics scheme """ @@ -51,12 +42,21 @@ function set_sedimentation_precomputed_quantities!(Y, p, t) @assert (p.atmos.moisture_model isa NonEquilMoistModel) (; ᶜwₗ, ᶜwᵢ) = p.precomputed - - FT = eltype(Y) + cmc = CAP.microphysics_cloud_params(p.params) + FT = eltype(p.params) # compute the precipitation terminal velocity [m/s] - # TODO - the actual parameterization will be added in the next PR - @. ᶜwₗ = FT(0) - @. ᶜwᵢ = FT(0) + @. ᶜwₗ = FT(0) #CMNe.terminal_velocity( + # cmc.liquid, + # cmc.Ch2022.rain, + # Y.c.ρ, + # max(zero(Y.c.ρ), Y.c.ρq_liq / Y.c.ρ), + #) + @. ᶜwᵢ = FT(0) #CMNe.terminal_velocity( + # cmc.ice, + # cmc.Ch2022.small_ice, + # Y.c.ρ, + # max(zero(Y.c.ρ), Y.c.ρq_ice / Y.c.ρ), + #) return nothing end diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index 4859fe0388..aab8caae5a 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -64,19 +64,26 @@ function precomputed_quantities(Y, atmos) sedimentation_quantities = atmos.moisture_model isa NonEquilMoistModel ? (; ᶜwₗ = similar(Y.c, FT), ᶜwᵢ = similar(Y.c, FT)) : (;) + precipitation_quantities = + atmos.precip_model isa Microphysics1Moment ? + (; ᶜwᵣ = similar(Y.c, FT), ᶜwₛ = similar(Y.c, FT)) : (;) precipitation_sgs_quantities = atmos.precip_model isa Microphysics0Moment ? (; ᶜSqₜᵖʲs = similar(Y.c, NTuple{n, FT}), ᶜSqₜᵖ⁰ = similar(Y.c, FT)) : atmos.precip_model isa Microphysics1Moment ? (; - ᶜSeₜᵖʲs = similar(Y.c, NTuple{n, FT}), - ᶜSqₜᵖʲs = similar(Y.c, NTuple{n, FT}), + ᶜSqₗᵖʲs = similar(Y.c, NTuple{n, FT}), + ᶜSqᵢᵖʲs = similar(Y.c, NTuple{n, FT}), ᶜSqᵣᵖʲs = similar(Y.c, NTuple{n, FT}), ᶜSqₛᵖʲs = similar(Y.c, NTuple{n, FT}), - ᶜSeₜᵖ⁰ = similar(Y.c, FT), - ᶜSqₜᵖ⁰ = similar(Y.c, FT), + ᶜSqₗᵖ⁰ = similar(Y.c, FT), + ᶜSqᵢᵖ⁰ = similar(Y.c, FT), ᶜSqᵣᵖ⁰ = similar(Y.c, FT), ᶜSqₛᵖ⁰ = similar(Y.c, FT), + ᶜq_liq⁰ = similar(Y.c, FT), + ᶜq_ice⁰ = similar(Y.c, FT), + ᶜq_rai⁰ = similar(Y.c, FT), + ᶜq_sno⁰ = similar(Y.c, FT), ) : (;) advective_sgs_quantities = atmos.turbconv_model isa PrognosticEDMFX ? @@ -111,6 +118,8 @@ function precomputed_quantities(Y, atmos) ᶜgradᵥ_θ_virt⁰ = Fields.Field(C3{FT}, cspace), ᶜgradᵥ_q_tot⁰ = Fields.Field(C3{FT}, cspace), ᶜgradᵥ_θ_liq_ice⁰ = Fields.Field(C3{FT}, cspace), + ᶜwₜʲs = similar(Y.c, NTuple{n, Geometry.WVector{FT}}), + ᶜwₕʲs = similar(Y.c, NTuple{n, Geometry.WVector{FT}}), precipitation_sgs_quantities..., ) : (;) sgs_quantities = (; @@ -152,14 +161,6 @@ function precomputed_quantities(Y, atmos) else (;) end - precipitation_quantities = - 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), - ) : (;) smagorinsky_lilly_quantities = if atmos.smagorinsky_lilly isa SmagorinskyLilly uvw_vec = UVW(FT(0), FT(0), FT(0)) @@ -313,30 +314,43 @@ 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 + elseif moisture_model isa EquilMoistModel && + precip_model isa NoPrecipitation (; specific.q_tot) - elseif moisture_model isa NonEquilMoistModel - q_pt_args = (specific.q_tot, specific.q_liq, specific.q_ice) + elseif moisture_model isa EquilMoistModel && + precip_model isa Microphysics0Moment + (; specific.q_tot) + elseif moisture_model isa NonEquilMoistModel && + precip_model isa Microphysics1Moment + 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...)) + else + error("Unsupported moisture and precipitation combination") 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) @@ -470,7 +484,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 @@ -526,27 +540,23 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) (; ᶜwₜqₜ, ᶜwₕhₜ) = p.precomputed @. ᶜwₜqₜ = Geometry.WVector(0) @. ᶜwₕhₜ = Geometry.WVector(0) - # - # TODO - uncomment in the next PR. Right now for the purpose of testing - # we want to merge with 0 sedimentation and precipitation - # if moisture_model isa NonEquilMoistModel set_sedimentation_precomputed_quantities!(Y, p, t) - # (; ᶜwₗ, ᶜwᵢ) = p.precomputed - # @. ᶜwₜqₜ += Geometry.WVector(ᶜwₗ * Y.c.ρq_liq + ᶜwᵢ * Y.c.ρq_ice) / Y.c.ρ - # @. ᶜwₕhₜ += Geometry.WVector( - # ᶜwₗ * Y.c.ρq_liq * (TD.internal_energy_liquid(thermo_params, ᶜts) + ᶜΦ + norm_sqr(Geometry.UVWVector(0, 0, -ᶜwₗ) + Geometry.UVWVector(ᶜu))/2) + - # ᶜwᵢ * Y.c.ρq_ice * (TD.internal_energy_ice(thermo_params, ᶜts) + ᶜΦ + norm_sqr(Geometry.UVWVector(0, 0, -ᶜwᵢ) + Geometry.UVWVector(ᶜu))/2) - # ) / Y.c.ρ + #(; ᶜwₗ, ᶜwᵢ) = p.precomputed + #@. ᶜwₜqₜ += Geometry.WVector(ᶜwₗ * Y.c.ρq_liq + ᶜwᵢ * Y.c.ρq_ice) / Y.c.ρ + #@. ᶜwₕhₜ += Geometry.WVector( + # ᶜwₗ * Y.c.ρq_liq * (TD.internal_energy_liquid(thermo_params, ᶜts) + ᶜΦ - norm_sqr(Geometry.UVWVector(0, 0, ᶜwₗ) + Geometry.UVWVector(ᶜu))/2) + + # ᶜwᵢ * Y.c.ρq_ice * (TD.internal_energy_ice(thermo_params, ᶜts) + ᶜΦ - norm_sqr(Geometry.UVWVector(0, 0, ᶜwᵢ) + Geometry.UVWVector(ᶜu))/2) + #) / Y.c.ρ end if precip_model isa Microphysics1Moment set_precipitation_precomputed_quantities!(Y, p, t) - # (; ᶜwᵣ, ᶜwₛ) = p.precomputed - # @. ᶜwₜqₜ += Geometry.WVector(ᶜwᵣ * Y.c.ρq_rai + ᶜwₛ * Y.c.ρq_sno) / Y.c.ρ - # @. ᶜwₕhₜ += Geometry.WVector( - # ᶜwᵣ * Y.c.ρq_rai * (TD.internal_energy_liquid(thermo_params, ᶜts) + ᶜΦ + norm_sqr(Geometry.UVWVector(0, 0, -ᶜwᵣ) + Geometry.UVWVector(ᶜu))/2) + - # ᶜwₛ * Y.c.ρq_sno * (TD.internal_energy_ice(thermo_params, ᶜts) + ᶜΦ + norm_sqr(Geometry.UVWVector(0, 0, -ᶜwₛ) + Geometry.UVWVector(ᶜu))/2) - # ) / Y.c.ρ + (; ᶜwᵣ, ᶜwₛ) = p.precomputed + #@. ᶜwₜqₜ += Geometry.WVector(ᶜwᵣ * Y.c.ρq_rai + ᶜwₛ * Y.c.ρq_sno) / Y.c.ρ + #@. ᶜwₕhₜ += Geometry.WVector( + # ᶜwᵣ * Y.c.ρq_rai * (TD.internal_energy_liquid(thermo_params, ᶜts) + ᶜΦ - norm_sqr(Geometry.UVWVector(0, 0, ᶜwᵣ) + Geometry.UVWVector(ᶜu))/2) + + # ᶜwₛ * Y.c.ρq_sno * (TD.internal_energy_ice(thermo_params, ᶜts) + ᶜΦ - norm_sqr(Geometry.UVWVector(0, 0, ᶜwₛ) + Geometry.UVWVector(ᶜu))/2) + #) / Y.c.ρ end if turbconv_model isa PrognosticEDMFX diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl index cf6b889f41..ee26b504a8 100644 --- a/src/cache/prognostic_edmf_precomputed_quantities.jl +++ b/src/cache/prognostic_edmf_precomputed_quantities.jl @@ -5,6 +5,9 @@ import NVTX import Thermodynamics as TD import ClimaCore: Spaces, Fields +#TODO - what to do about prognostic and diagnostic edmf precipitation +#TODO - we only want to support dry + nothing, equil + 0M and nonequil + 1M + """ set_prognostic_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t) @@ -19,11 +22,14 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_environment!( @assert !(p.atmos.moisture_model isa DryModel) thermo_params = CAP.thermodynamics_params(p.params) - (; turbconv_model) = p.atmos + (; turbconv_model, moisture_model) = p.atmos (; ᶜΦ,) = p.core (; ᶜp, ᶜh_tot, ᶜK) = p.precomputed (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜmse⁰, ᶜq_tot⁰) = p.precomputed + if moisture_model isa NonEquilMoistModel + (; ᶜq_liq⁰, ᶜq_ice⁰, ᶜq_rai⁰, ᶜq_sno⁰) = p.precomputed + end @. ᶜρa⁰ = ρa⁰(Y.c) @. ᶜtke⁰ = divide_by_ρa(Y.c.sgs⁰.ρatke, ᶜρa⁰, 0, Y.c.ρ, turbconv_model) @@ -41,10 +47,50 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_environment!( Y.c.ρ, turbconv_model, ) + if moisture_model isa NonEquilMoistModel + @. ᶜq_liq⁰ = divide_by_ρa( + Y.c.ρq_liq - ρaq_liq⁺(Y.c.sgsʲs), + ᶜρa⁰, + Y.c.ρq_liq, + Y.c.ρ, + turbconv_model, + ) + @. ᶜq_ice⁰ = divide_by_ρa( + Y.c.ρq_ice - ρaq_ice⁺(Y.c.sgsʲs), + ᶜρa⁰, + Y.c.ρq_ice, + Y.c.ρ, + turbconv_model, + ) + @. ᶜq_rai⁰ = divide_by_ρa( + Y.c.ρq_rai - ρaq_rai⁺(Y.c.sgsʲs), + ᶜρa⁰, + Y.c.ρq_rai, + Y.c.ρ, + turbconv_model, + ) + @. ᶜq_sno⁰ = divide_by_ρa( + Y.c.ρq_sno - ρaq_sno⁺(Y.c.sgsʲs), + ᶜρa⁰, + Y.c.ρq_sno, + Y.c.ρ, + turbconv_model, + ) + end + set_sgs_ᶠu₃!(u₃⁰, ᶠu₃⁰, Y, turbconv_model) set_velocity_quantities!(ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶠu₃⁰, Y.c.uₕ, ᶠuₕ³) # @. ᶜK⁰ += ᶜtke⁰ - @. ᶜts⁰ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmse⁰ - ᶜΦ, ᶜq_tot⁰) + if moisture_model isa EquilMoistModel + @. ᶜts⁰ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmse⁰ - ᶜΦ, ᶜq_tot⁰) + else + @. ᶜts⁰ = TD.PhaseNonEquil_phq( + thermo_params, + ᶜp, + ᶜmse⁰ - ᶜΦ, + TD.PhasePartition(ᶜq_tot⁰, ᶜq_liq⁰ + ᶜq_rai⁰, ᶜq_ice⁰ + ᶜq_sno⁰), + ) + end @. ᶜρ⁰ = TD.air_density(thermo_params, ᶜts⁰) return nothing end @@ -87,10 +133,29 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft_and_bc! ᶜρʲ = ᶜρʲs.:($j) ᶜmseʲ = Y.c.sgsʲs.:($j).mse ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot + if moisture_model isa NonEquilMoistModel + ᶜq_liqʲ = Y.c.sgsʲs.:($j).q_liq + ᶜq_iceʲ = Y.c.sgsʲs.:($j).q_ice + ᶜq_raiʲ = Y.c.sgsʲs.:($j).q_rai + ᶜq_snoʲ = Y.c.sgsʲs.:($j).q_sno + end set_velocity_quantities!(ᶜuʲ, ᶠu³ʲ, ᶜKʲ, ᶠu₃ʲ, Y.c.uₕ, ᶠuₕ³) @. ᶠKᵥʲ = (adjoint(CT3(ᶠu₃ʲ)) * ᶠu₃ʲ) / 2 - @. ᶜtsʲ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmseʲ - ᶜΦ, ᶜq_totʲ) + if moisture_model isa EquilMoistModel + @. ᶜtsʲ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜmseʲ - ᶜΦ, ᶜq_totʲ) + else + @. ᶜtsʲ = TD.PhaseNonEquil_phq( + thermo_params, + ᶜp, + ᶜmseʲ - ᶜΦ, + TD.PhasePartition( + ᶜq_totʲ, + ᶜq_liqʲ + ᶜq_raiʲ, + ᶜq_iceʲ + ᶜq_snoʲ, + ), + ) + end @. ᶜρʲ = TD.air_density(thermo_params, ᶜtsʲ) # EDMFX boundary condition: @@ -396,8 +461,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation # Sources from the updrafts n = n_mass_flux_subdomains(p.atmos.turbconv_model) for j in 1:n - @. ᶜSqₜᵖʲs.:($$j) = q_tot_precipitation_sources( - Microphysics0Moment(), + @. ᶜSqₜᵖʲs.:($$j) = q_tot_0M_precipitation_sources( thp, cmp, dt, @@ -406,14 +470,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation ) end # sources from the environment - @. ᶜSqₜᵖ⁰ = q_tot_precipitation_sources( - Microphysics0Moment(), - thp, - cmp, - dt, - ᶜq_tot⁰, - ᶜts⁰, - ) + @. ᶜSqₜᵖ⁰ = q_tot_0M_precipitation_sources(thp, cmp, dt, ᶜq_tot⁰, ᶜts⁰) return nothing end NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation!( @@ -428,9 +485,14 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation thp = CAP.thermodynamics_params(params) cmp = CAP.microphysics_1m_params(params) - (; ᶜSeₜᵖʲs, ᶜSqₜᵖʲs, ᶜSqᵣᵖʲs, ᶜSqₛᵖʲs, ᶜρʲs, ᶜtsʲs) = p.precomputed - (; ᶜSeₜᵖ⁰, ᶜSqₜᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰, ᶜρ⁰, ᶜts⁰) = p.precomputed - (; ᶜqᵣ, ᶜqₛ) = p.precomputed + (; ᶜSqₗᵖʲs, ᶜSqᵢᵖʲs, ᶜSqᵣᵖʲs, ᶜSqₛᵖʲs, ᶜρʲs, ᶜtsʲs) = p.precomputed + (; ᶜSqₗᵖ⁰, ᶜSqᵢᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰, ᶜρ⁰, ᶜts⁰) = p.precomputed + (; ᶜq_rai⁰, ᶜq_sno⁰) = p.precomputed + + (; ᶜuʲs) = p.precomputed + + + (; ᶜwₗʲs, ᶜwᵢʲs, ᶜwᵣʲs, ᶜwₛʲs) = p.precomputed # TODO - can I re-use them between js and env? ᶜSᵖ = p.scratch.ᶜtemp_scalar @@ -443,34 +505,105 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation compute_precipitation_sources!( ᶜSᵖ, ᶜSᵖ_snow, - ᶜSqₜᵖʲs.:($j), + ᶜSqₗᵖʲs.:($j), + ᶜSqᵢᵖʲs.:($j), ᶜSqᵣᵖʲs.:($j), ᶜSqₛᵖʲs.:($j), - ᶜSeₜᵖʲs.:($j), ᶜρʲs.:($j), - ᶜqᵣ, - ᶜqₛ, + Y.c.sgsʲs.:($j).q_rai, + Y.c.sgsʲs.:($j).q_sno, ᶜtsʲs.:($j), - ᶜΦ, dt, cmp, thp, ) + + @. ᶜwₗʲs.:($$j) = Geometry.WVector(0) # TODO + @. ᶜwᵢʲs.:($$j) = Geometry.WVector(0) # TODO + @. ᶜwᵣʲs.:($$j) = Geometry.WVector(0) # TODO + @. ᶜwₛʲs.:($$j) = Geometry.WVector(0) # TODO + + + @. ᶜwₜʲs.:($$j) += ifelse( + Y.c.sgsʲs.q_tot.:($$j) == 0, + 0, + ( + ᶜwₗʲs.:($$j) * Y.c.sgsʲs.q_liq.:($$j) + + ᶜwᵢʲs.:($$j) * Y.c.sgsʲs.q_ice.:($$j) + ) / Y.c.sgsʲs.q_tot.:($$j), + ) + @. ᶜwₜʲs.:($$j) += ifelse( + Y.c.sgsʲs.q_tot.:($$j) == 0, + 0, + ( + ᶜwᵣʲs.:($$j) * Y.c.sgsʲs.q_rai.:($$j) + + ᶜwₛʲs.:($$j) * Y.c.sgsʲs.q_sno.:($$j) + ) / Y.c.sgsʲs.q_tot.:($$j), + ) + + #TODO - what if denominator is zero? + @. ᶜwₕʲs += + Geometry.WVector( + ᶜwₗʲs.:($$j) * + Y.c.sgsʲs.q_liq.:($$j) * + ( + TD.internal_energy_liquid(thermo_params, ᶜtsʲs.:($$j)) + + p.core.ᶜΦ + + norm_sqr( + Geometry.UVWVector(0, 0, -1 * ᶜwₗʲs.:($$j)) + + Geometry.UVWVector(ᶜuʲs), + ) / 2 + ) + + ᶜwᵢʲs.:($$j) * + Y.c.sgsʲs.q_ice.:($$j) * + ( + TD.internal_energy_ice(thermo_params, ᶜtsʲs.:($$j)) + + p.core.ᶜΦ + + norm_sqr( + Geometry.UVWVector(0, 0, -1 * ᶜwᵢʲs.:($$j)) + + Geometry.UVWVector(ᶜuʲs), + ) / 2 + ), + ) / (Y.c.sgsʲs.:($j).mse - p.coreᶜΦ) + @. ᶜwₕʲs += + Geometry.WVector( + ᶜwᵣ * + Y.c.ρq_rai * + ( + TD.internal_energy_liquid(thermo_params, ᶜts) + + ᶜΦ + + norm_sqr( + Geometry.UVWVector(0, 0, -1 * ᶜwᵣ) + + Geometry.UVWVector(ᶜuʲs), + ) / 2 + ) + + ᶜwₛ * + Y.c.ρq_sno * + ( + TD.internal_energy_ice(thermo_params, ᶜts) + + ᶜΦ + + norm_sqr( + Geometry.UVWVector(0, 0, -1 * ᶜwₛ) + + Geometry.UVWVector(ᶜuʲs), + ) / 2 + ), + ) / (Y.c.sgsʲs.:($j).mse - p.coreᶜΦ) + + end # Sources from the environment compute_precipitation_sources!( ᶜSᵖ, ᶜSᵖ_snow, - ᶜSqₜᵖ⁰, + ᶜSqₗᵖ⁰, + ᶜSqᵢᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰, - ᶜSeₜᵖ⁰, ᶜρ⁰, - ᶜqᵣ, - ᶜqₛ, + ᶜq_rai⁰, + ᶜq_sno⁰, ᶜts⁰, - ᶜΦ, dt, cmp, thp, diff --git a/src/initial_conditions/initial_conditions.jl b/src/initial_conditions/initial_conditions.jl index 8e61c1f81c..a83b20db0f 100644 --- a/src/initial_conditions/initial_conditions.jl +++ b/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,8 +419,8 @@ 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, @@ -428,9 +428,9 @@ function overwrite_initial_conditions!( ᶜ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 @@ -438,7 +438,7 @@ function overwrite_initial_conditions!( # 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) @@ -1260,7 +1260,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, diff --git a/src/parameterized_tendencies/microphysics/cloud_condensate.jl b/src/parameterized_tendencies/microphysics/cloud_condensate.jl index f4739fb0e6..6fe3379499 100644 --- a/src/parameterized_tendencies/microphysics/cloud_condensate.jl +++ b/src/parameterized_tendencies/microphysics/cloud_condensate.jl @@ -2,19 +2,34 @@ ##### DryModel, EquilMoistModel ##### -cloud_condensate_tendency!(Yₜ, p, _) = nothing +cloud_condensate_tendency!(Yₜ, p, _, _) = nothing ##### ##### NonEquilMoistModel ##### -function cloud_condensate_tendency!(Yₜ, p, ::NonEquilMoistModel) - +function cloud_condensate_tendency!( + Yₜ, + p, + ::NonEquilMoistModel, + ::Union{NoPrecipitation, Microphysics0Moment}, +) + error( + "NonEquilMoistModel can only be run with Microphysics1Moment precipitation", + ) +end +function cloud_condensate_tendency!( + Yₜ, + p, + ::NonEquilMoistModel, + ::Microphysics1Moment, +) (; ᶜts) = p.precomputed + (; q_rai, q_sno) = p.precomputed.ᶜspecific (; params, dt) = p 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 diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl index 341a3a26bb..9eaf425656 100644 --- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl +++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl @@ -8,17 +8,13 @@ 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 -const Iᵢ = TD.internal_energy_ice -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 -cᵥₗ(thp) = TD.Parameters.cv_l(thp) -cᵥᵢ(thp) = TD.Parameters.cv_i(thp) + +qₗ(thp, ts, qᵣ) = max(0, TD.PhasePartition(thp, ts).liq - qᵣ) +qᵢ(thp, ts, qₛ) = max(0, TD.PhasePartition(thp, ts).ice - qₛ) # helper function to limit the tendency function limit(q, dt, n::Int) @@ -36,7 +32,13 @@ end Returns the condensation/evaporation or deposition/sublimation rate for non-equilibrium Morrison and Milbrandt 2015 cloud formation. """ -function cloud_sources(cm_params::CMP.CloudLiquid{FT}, thp, ts, dt) where {FT} +function cloud_sources( + cm_params::CMP.CloudLiquid{FT}, + thp, + ts, + qᵣ, + dt, +) where {FT} q = TD.PhasePartition(thp, ts) ρ = TD.air_density(thp, ts) @@ -47,10 +49,10 @@ function cloud_sources(cm_params::CMP.CloudLiquid{FT}, thp, ts, dt) where {FT} return ifelse( S > FT(0), min(S, limit(qᵥ(thp, ts), dt, 2)), - -min(abs(S), limit(qₗ(thp, ts), dt, 2)), + -min(abs(S), limit(qₗ(thp, ts, qᵣ), dt, 2)), ) end -function cloud_sources(cm_params::CMP.CloudIce{FT}, thp, ts, dt) where {FT} +function cloud_sources(cm_params::CMP.CloudIce{FT}, thp, ts, qₛ, dt) where {FT} q = TD.PhasePartition(thp, ts) ρ = TD.air_density(thp, ts) @@ -61,33 +63,22 @@ function cloud_sources(cm_params::CMP.CloudIce{FT}, thp, ts, dt) where {FT} return ifelse( S > FT(0), min(S, limit(qᵥ(thp, ts), dt, 2)), - -min(abs(S), limit(qᵢ(thp, ts), dt, 2)), + -min(abs(S), limit(qᵢ(thp, ts, qₛ), dt, 2)), ) end """ - q_tot_precipitation_sources(precip_model, thp, cmp, dt, qₜ, ts) + q_tot_0M_precipitation_sources(thp, cmp, dt, qₜ, ts) - - precip_model - a type for precipitation scheme choice - thp, cmp - structs with thermodynamic and microphysics parameters - dt - model time step - qₜ - total water specific humidity - ts - thermodynamic state (see Thermodynamics.jl package for details) Returns the qₜ source term due to precipitation formation -defined as Δm_tot / (m_dry + m_tot) +defined as Δm_tot / (m_dry + m_tot) for the 0-moment scheme """ -function q_tot_precipitation_sources(::NoPrecipitation, thp, cmp, dt, qₜ, ts) - return zero(qₜ) -end -function q_tot_precipitation_sources( - ::Microphysics0Moment, - thp, - cmp::CMP.Parameters0M, - dt, - qₜ, - ts, -) +function q_tot_0M_precipitation_sources(thp, cmp::CMP.Parameters0M, dt, qₜ, ts) return -min(max(qₜ, 0) / dt, -CM0.remove_precipitation(cmp, PP(thp, ts))) end @@ -104,39 +95,38 @@ for the 0-moment scheme function e_tot_0M_precipitation_sources_helper(thp, ts, Φ) λ = TD.liquid_fraction(thp, ts) + Iₗ = TD.internal_energy_liquid(thp, ts) + Iᵢ = TD.internal_energy_ice(thp, ts) - return λ * Iₗ(thp, ts) + (1 - λ) * Iᵢ(thp, ts) + Φ + return λ * Iₗ + (1 - λ) * Iᵢ + Φ end """ - compute_precipitation_sources!(Sᵖ, Sᵖ_snow, Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ, ρ, qᵣ, qₛ, ts, Φ, dt, mp, thp) + compute_precipitation_sources!(Sᵖ, Sᵖ_snow, Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ, ρ, qᵣ, qₛ, ts, dt, mp, thp) - Sᵖ, Sᵖ_snow - temporary containters to help compute precipitation source terms - - Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ - cached storage for precipitation source terms + - Sqₗᵖ, Sqᵢᵖ, Sqᵣᵖ, Sqₛᵖ - cached storage for precipitation source terms - ρ - air density - qᵣ, qₛ - precipitation (rain and snow) specific humidity - ts - thermodynamic state (see td package for details) - - Φ - geopotential - dt - model time step - thp, cmp - structs with thermodynamic and microphysics parameters Returns the q source terms due to precipitation formation from the 1-moment scheme. The specific humidity source terms are defined as defined as Δmᵢ / (m_dry + m_tot) where i stands for total, rain or snow. -Also returns the total energy source term due to the microphysics processes. """ function compute_precipitation_sources!( Sᵖ, Sᵖ_snow, - Sqₜᵖ, + Sqₗᵖ, + Sqᵢᵖ, Sqᵣᵖ, Sqₛᵖ, - Seₜᵖ, ρ, qᵣ, qₛ, ts, - Φ, dt, mp, thp, @@ -144,155 +134,104 @@ function compute_precipitation_sources!( FT = eltype(thp) # @. Sqₜᵖ = FT(0) should work after fixing # https://github.com/CliMA/ClimaCore.jl/issues/1786 - @. Sqₜᵖ = ρ * FT(0) + @. Sqₗᵖ = ρ * FT(0) + @. Sqᵢᵖ = ρ * FT(0) @. Sqᵣᵖ = ρ * FT(0) @. Sqₛᵖ = ρ * FT(0) - @. Seₜᵖ = ρ * FT(0) #! format: off # rain autoconversion: q_liq -> q_rain @. Sᵖ = ifelse( mp.Ndp <= 0, - CM1.conv_q_liq_to_q_rai(mp.pr.acnv1M, qₗ(thp, ts), true), - CM2.conv_q_liq_to_q_rai(mp.var, qₗ(thp, ts), ρ, mp.Ndp), + CM1.conv_q_liq_to_q_rai(mp.pr.acnv1M, qₗ(thp, ts, qᵣ), true), + CM2.conv_q_liq_to_q_rai(mp.var, qₗ(thp, ts, qᵣ), ρ, mp.Ndp), ) - @. Sᵖ = min(limit(qₗ(thp, ts), dt, 5), Sᵖ) - @. Sqₜᵖ -= Sᵖ + @. Sᵖ = min(limit(qₗ(thp, ts, qᵣ), dt, 5), Sᵖ) + @. Sqₗᵖ -= Sᵖ @. Sqᵣᵖ += Sᵖ - @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ) # snow autoconversion assuming no supersaturation: q_ice -> q_snow @. Sᵖ = min( - limit(qᵢ(thp, ts), dt, 5), - CM1.conv_q_ice_to_q_sno_no_supersat(mp.ps.acnv1M, qᵢ(thp, ts), true), + limit(qᵢ(thp, ts, qₛ), dt, 5), + CM1.conv_q_ice_to_q_sno_no_supersat(mp.ps.acnv1M, qᵢ(thp, ts, qₛ), true), ) - @. Sqₜᵖ -= Sᵖ + @. Sqᵢᵖ -= Sᵖ @. Sqₛᵖ += Sᵖ - @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ) # accretion: q_liq + q_rain -> q_rain @. Sᵖ = min( - limit(qₗ(thp, ts), dt, 5), - CM1.accretion(mp.cl, mp.pr, mp.tv.rain, mp.ce, qₗ(thp, ts), qᵣ, ρ), + limit(qₗ(thp, ts, qᵣ), dt, 5), + CM1.accretion(mp.cl, mp.pr, mp.tv.rain, mp.ce, qₗ(thp, ts, qᵣ), qᵣ, ρ), ) - @. Sqₜᵖ -= Sᵖ + @. Sqₗᵖ -= Sᵖ @. Sqᵣᵖ += Sᵖ - @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ) # accretion: q_ice + q_snow -> q_snow @. Sᵖ = min( - limit(qᵢ(thp, ts), dt, 5), - CM1.accretion(mp.ci, mp.ps, mp.tv.snow, mp.ce, qᵢ(thp, ts), qₛ, ρ), - ) - @. Sqₜᵖ -= Sᵖ - @. Sqₛᵖ += Sᵖ - @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ) - - # 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, 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 - # else the snow melts and both cloud water and snow become rain - α(thp, ts) = cᵥₗ(thp) / Lf(thp, ts) * (Tₐ(thp, ts) - mp.ps.T_freeze) - @. Sᵖ_snow = ifelse( - Tₐ(thp, ts) < mp.ps.T_freeze, - Sᵖ, - FT(-1) * min(Sᵖ * α(thp, ts), limit(qₛ, dt, 5)), - ) - @. Sqₛᵖ += Sᵖ_snow - @. Sqₜᵖ -= Sᵖ - @. Sqᵣᵖ += ifelse(Tₐ(thp, ts) < mp.ps.T_freeze, FT(0), Sᵖ - Sᵖ_snow) - @. Seₜᵖ -= ifelse( - Tₐ(thp, ts) < mp.ps.T_freeze, - Sᵖ * (Iᵢ(thp, ts) + Φ), - Sᵖ * (Iₗ(thp, ts) + Φ) - Sᵖ_snow * (Iₗ(thp, ts) - Iᵢ(thp, ts)), - ) - - # accretion: q_ice + q_rai -> q_sno - @. Sᵖ = min( - limit(qᵢ(thp, ts), dt, 5), - CM1.accretion(mp.ci, mp.pr, mp.tv.rain, mp.ce, qᵢ(thp, ts), qᵣ, ρ), - ) - @. Sqₜᵖ -= Sᵖ - @. Sqₛᵖ += Sᵖ - @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ) - # sink of rain via accretion cloud ice - rain - @. Sᵖ = min( - limit(qᵣ, dt, 5), - CM1.accretion_rain_sink(mp.pr, mp.ci, mp.tv.rain, mp.ce, qᵢ(thp, ts), qᵣ, ρ), + limit(qᵢ(thp, ts, qₛ), dt, 5), + CM1.accretion(mp.ci, mp.ps, mp.tv.snow, mp.ce, qᵢ(thp, ts, qₛ), qₛ, ρ), ) - @. Sqᵣᵖ -= Sᵖ + @. Sqᵢᵖ -= Sᵖ @. Sqₛᵖ += Sᵖ - @. Seₜᵖ += Sᵖ * Lf(thp, ts) - # accretion: q_rai + q_sno -> q_rai or q_sno - @. Sᵖ = ifelse( - Tₐ(thp, ts) < mp.ps.T_freeze, - min( - 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, 5), - CM1.accretion_snow_rain(mp.pr, mp.ps, mp.tv.snow, mp.tv.rain, mp.ce, qᵣ, qₛ, ρ), - ), - ) - @. Sqₛᵖ += Sᵖ - @. Sqᵣᵖ -= Sᵖ - @. Seₜᵖ += Sᵖ * Lf(thp, ts) - #! format: on + ## 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, qᵣ), dt, 5), + # CM1.accretion(mp.cl, mp.ps, mp.tv.snow, mp.ce, qₗ(thp, ts, qᵣ), qₛ, ρ), + #) + ## if T < T_freeze cloud droplets freeze to become snow + ## else the snow melts and both cloud water and snow become rain + #α(thp, ts) = TD.Parameters.cv_l(thp) / TD.latent_heat_fusion(thp, ts) * (Tₐ(thp, ts) - mp.ps.T_freeze) + #@. Sᵖ_snow = ifelse( + # Tₐ(thp, ts) < mp.ps.T_freeze, + # Sᵖ, + # FT(-1) * min(Sᵖ * α(thp, ts), limit(qₛ, dt, 5)), + #) + #@. Sqₛᵖ += Sᵖ_snow + #@. Sqₗᵖ -= Sᵖ + #@. Sqᵣᵖ += ifelse(Tₐ(thp, ts) < mp.ps.T_freeze, FT(0), Sᵖ - Sᵖ_snow) + + ## accretion: q_ice + q_rai -> q_sno + #@. Sᵖ = min( + # limit(qᵢ(thp, ts, qₛ), dt, 5), + # CM1.accretion(mp.ci, mp.pr, mp.tv.rain, mp.ce, qᵢ(thp, ts, qₛ), qᵣ, ρ), + #) + #@. Sqᵢᵖ -= Sᵖ + #@. Sqₛᵖ += Sᵖ + ## sink of rain via accretion cloud ice - rain + #@. Sᵖ = min( + # limit(qᵣ, dt, 5), + # CM1.accretion_rain_sink(mp.pr, mp.ci, mp.tv.rain, mp.ce, qᵢ(thp, ts, qₛ), qᵣ, ρ), + #) + #@. Sqᵣᵖ -= Sᵖ + #@. Sqₛᵖ += Sᵖ + + ## accretion: q_rai + q_sno -> q_rai or q_sno + #@. Sᵖ = ifelse( + # Tₐ(thp, ts) < mp.ps.T_freeze, + # min( + # 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, 5), + # CM1.accretion_snow_rain(mp.pr, mp.ps, mp.tv.snow, mp.tv.rain, mp.ce, qᵣ, qₛ, ρ), + # ), + #) + #@. Sqₛᵖ += Sᵖ + #@. Sqᵣᵖ -= Sᵖ + ##! format: on end """ - compute_precipitation_heating(Seₜᵖ, ᶜwᵣ, ᶜwₛ, ᶜu, qᵣ, qₛ, ᶜts, thp) - - - Seₜᵖ - cached storage for precipitation energy source terms - - ᶜwᵣ, ᶜwₛ - rain and snow terminal velocities - - ᶜu - air velocity - - qᵣ, qₛ - precipitation (rain and snow) specific humidities - - ᶜts - thermodynamic state (see td package for details) - - ᶜ∇T - cached temporary variable to store temperature gradient - - thp - structs with thermodynamic and microphysics parameters - - Augments the energy source terms with heat exchange between air - and precipitating species, following eq. 36 from Raymond 2013 - doi:10.1002/jame.20050 and assuming that precipitation has the same - temperature as the surrounding air. -""" -function compute_precipitation_heating!( - ᶜSeₜᵖ, - ᶜwᵣ, - ᶜwₛ, - ᶜu, - ᶜqᵣ, - ᶜqₛ, - ᶜts, - ᶜ∇T, - thp, -) - # TODO - at some point we want to switch to assuming that precipitation - # is at wet bulb temperature - - # 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ᵣ - @. ᶜ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) + compute_precipitation_sinks!(Sᵖ, Sqᵣᵖ, Sqₛᵖ, ρ, qᵣ, qₛ, ts, dt, mp, thp) - Sᵖ - a temporary containter to help compute precipitation source terms - - Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ - cached storage for precipitation source terms + - Sqᵣᵖ, Sqₛᵖ - cached storage for precipitation source terms - ρ - air density - qᵣ, qₛ - precipitation (rain and snow) specific humidities - ts - thermodynamic state (see td package for details) - - Φ - geopotential - dt - model time step - thp, cmp - structs with thermodynamic and microphysics parameters @@ -303,51 +242,43 @@ Also returns the total energy source term due to the microphysics processes. """ function compute_precipitation_sinks!( Sᵖ, - Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, - Seₜᵖ, ρ, qᵣ, qₛ, ts, - Φ, dt, mp, thp, ) - FT = eltype(Sqₜᵖ) + FT = eltype(Sᵖ) sps = (mp.ps, mp.tv.snow, mp.aps, thp) rps = (mp.pr, mp.tv.rain, mp.aps, thp) - #! format: off - # evaporation: q_rai -> q_vap - @. Sᵖ = -min( - limit(qᵣ, dt, 5), - -CM1.evaporation_sublimation(rps..., PP(thp, ts), qᵣ, ρ, Tₐ(thp, ts)), - ) - @. Sqₜᵖ -= Sᵖ - @. Sqᵣᵖ += Sᵖ - @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ) - - # melting: q_sno -> q_rai - @. Sᵖ = min( - limit(qₛ, dt, 5), - CM1.snow_melt(sps..., qₛ, ρ, Tₐ(thp, ts)), - ) - @. Sqᵣᵖ += Sᵖ - @. Sqₛᵖ -= Sᵖ - @. Seₜᵖ -= Sᵖ * Lf(thp, ts) - - # deposition/sublimation: q_vap <-> q_sno - @. Sᵖ = CM1.evaporation_sublimation(sps..., PP(thp, ts), qₛ, ρ, Tₐ(thp, ts)) - @. Sᵖ = ifelse( - Sᵖ > FT(0), - min(limit(qᵥ(thp, ts), dt, 5), Sᵖ), - -min(limit(qₛ, dt, 5), FT(-1) * Sᵖ), - ) - @. Sqₜᵖ -= Sᵖ - @. Sqₛᵖ += Sᵖ - @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ) - #! format: on + ##! format: off + ## evaporation: q_rai -> q_vap + #@. Sᵖ = -min( + # limit(qᵣ, dt, 5), + # -CM1.evaporation_sublimation(rps..., PP(thp, ts), qᵣ, ρ, Tₐ(thp, ts)), + #) + #@. Sqᵣᵖ += Sᵖ + + ## melting: q_sno -> q_rai + #@. Sᵖ = min( + # limit(qₛ, dt, 5), + # CM1.snow_melt(sps..., qₛ, ρ, Tₐ(thp, ts)), + #) + #@. Sqᵣᵖ += Sᵖ + #@. Sqₛᵖ -= Sᵖ + + ## deposition/sublimation: q_vap <-> q_sno + #@. Sᵖ = CM1.evaporation_sublimation(sps..., PP(thp, ts), qₛ, ρ, Tₐ(thp, ts)) + #@. Sᵖ = ifelse( + # Sᵖ > FT(0), + # min(limit(qᵥ(thp, ts), dt, 5), Sᵖ), + # -min(limit(qₛ, dt, 5), FT(-1) * Sᵖ), + #) + #@. Sqₛᵖ += Sᵖ + ##! format: on end diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl index 86606a5253..89aed8218e 100644 --- a/src/parameterized_tendencies/microphysics/precipitation.jl +++ b/src/parameterized_tendencies/microphysics/precipitation.jl @@ -24,7 +24,7 @@ function precipitation_cache(Y, precip_model::NoPrecipitation) surface_snow_flux = zeros(axes(Fields.level(Y.f, half))), ) end -precipitation_tendency!(Yₜ, Y, p, t, ::NoPrecipitation, _) = nothing +precipitation_tendency!(Yₜ, Y, p, t, _, ::NoPrecipitation, _) = nothing ##### ##### 0-Moment without sgs scheme or with diagnostic/prognostic edmf @@ -50,8 +50,7 @@ function compute_precipitation_cache!(Y, p, ::Microphysics0Moment, _) cm_params = CAP.microphysics_0m_params(params) thermo_params = CAP.thermodynamics_params(params) @. ᶜS_ρq_tot = - Y.c.ρ * q_tot_precipitation_sources( - Microphysics0Moment(), + Y.c.ρ * q_tot_0M_precipitation_sources( thermo_params, cm_params, dt, @@ -158,10 +157,21 @@ function precipitation_tendency!( Y, p, t, - precip_model::Microphysics0Moment, + ::DryModel, + ::Microphysics0Moment, _, ) - (; turbconv_model) = p.atmos + error("Microphysics0Moment precipitation should not be run with DryModel.") +end +function precipitation_tendency!( + Yₜ, + Y, + p, + t, + ::EquilMoistModel, + precip_model::Microphysics0Moment, + turbconv_model, +) (; ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precipitation # Compute the ρq_tot and ρe_tot precipitation source terms @@ -176,7 +186,19 @@ function precipitation_tendency!( return nothing end - +function precipitation_tendency!( + Yₜ, + Y, + p, + t, + ::NonEquilMoistModel, + ::Microphysics0Moment, + _, +) + error( + "Microphysics0Moment precipitation and NonEquilibriumMost model precipitation_tendency has not been implemented.", + ) +end ##### ##### 1-Moment without sgs scheme ##### @@ -184,10 +206,10 @@ end function precipitation_cache(Y, precip_model::Microphysics1Moment) FT = Spaces.undertype(axes(Y.c)) return (; - ᶜSqₜᵖ = similar(Y.c, FT), + ᶜSqₗᵖ = similar(Y.c, FT), + ᶜSqᵢᵖ = similar(Y.c, FT), ᶜSqᵣᵖ = similar(Y.c, FT), ᶜSqₛᵖ = similar(Y.c, FT), - ᶜSeₜᵖ = similar(Y.c, FT), surface_rain_flux = zeros(axes(Fields.level(Y.f, half))), surface_snow_flux = zeros(axes(Fields.level(Y.f, half))), ) @@ -196,13 +218,12 @@ end function compute_precipitation_cache!(Y, p, ::Microphysics1Moment, _) FT = Spaces.undertype(axes(Y.c)) (; dt) = p - (; ᶜts, ᶜqᵣ, ᶜqₛ, ᶜwᵣ, ᶜwₛ, ᶜu) = p.precomputed - (; ᶜΦ) = p.core - (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation + (; ᶜts) = p.precomputed + (; ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precipitation + (; q_rai, q_sno) = p.precomputed.ᶜspecific ᶜ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 @@ -213,15 +234,14 @@ function compute_precipitation_cache!(Y, p, ::Microphysics1Moment, _) compute_precipitation_sources!( ᶜSᵖ, ᶜSᵖ_snow, - ᶜSqₜᵖ, + ᶜSqₗᵖ, + ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, - ᶜSeₜᵖ, Y.c.ρ, - ᶜqᵣ, - ᶜqₛ, + q_rai, + q_sno, ᶜts, - ᶜΦ, dt, cmp, thp, @@ -231,21 +251,16 @@ function compute_precipitation_cache!(Y, p, ::Microphysics1Moment, _) # (For now only done on the grid mean) compute_precipitation_sinks!( ᶜSᵖ, - ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, - ᶜSeₜᵖ, Y.c.ρ, - ᶜqᵣ, - ᶜqₛ, + q_rai, + q_sno, ᶜts, - ᶜΦ, dt, cmp, thp, ) - # first term of eq 36 from Raymond 2013 - compute_precipitation_heating!(ᶜSeₜᵖ, ᶜwᵣ, ᶜwₛ, ᶜu, ᶜqᵣ, ᶜqₛ, ᶜts, ᶜ∇T, thp) end function compute_precipitation_cache!( Y, @@ -255,13 +270,12 @@ function compute_precipitation_cache!( ) FT = Spaces.undertype(axes(Y.c)) (; dt) = p - (; ᶜts, ᶜqᵣ, ᶜqₛ, ᶜwᵣ, ᶜwₛ, ᶜu) = p.precomputed - (; ᶜΦ) = p.core + (; ᶜts) = p.precomputed + (; q_rai, q_sno) = p.precomputed.ᶜspecific # Grid mean precipitation sinks - (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation + (; ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precipitation # additional scratch storage ᶜSᵖ = p.scratch.ᶜtemp_scalar - ᶜ∇T = p.scratch.ᶜtemp_CT123 # get thermodynamics and 1-moment microphysics params (; params) = p @@ -269,29 +283,22 @@ function compute_precipitation_cache!( thp = CAP.thermodynamics_params(params) # zero out the helper source terms - @. ᶜSqₜᵖ = FT(0) @. ᶜSqᵣᵖ = FT(0) @. ᶜSqₛᵖ = FT(0) - @. ᶜSeₜᵖ = FT(0) # compute precipitation sinks # (For now only done on the grid mean) compute_precipitation_sinks!( ᶜSᵖ, - ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, - ᶜSeₜᵖ, Y.c.ρ, - ᶜqᵣ, - ᶜqₛ, + q_rai, + q_sno, ᶜts, - ᶜΦ, dt, cmp, thp, ) - # first term of eq 36 from Raymond 2013 - compute_precipitation_heating!(ᶜSeₜᵖ, ᶜwᵣ, ᶜwₛ, ᶜu, ᶜqᵣ, ᶜqₛ, ᶜts, ᶜ∇T, thp) end function compute_precipitation_surface_fluxes!( @@ -333,23 +340,46 @@ function precipitation_tendency!( Y, p, t, + ::DryModel, + precip_model::Microphysics1Moment, + _, +) + error("Microphysics1Moment precipitation should not be used with DryModel") +end +function precipitation_tendency!( + Yₜ, + Y, + p, + t, + ::EquilMoistModel, + precip_model::Microphysics1Moment, + _, +) + error( + "Microphysics1Moment precipitation and EquilMoistModel precipitation_tendency is not implemented", + ) +end +function precipitation_tendency!( + Yₜ, + Y, + p, + t, + ::NonEquilMoistModel, precip_model::Microphysics1Moment, _, ) (; turbconv_model) = p.atmos - (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation + (; ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = 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_liq += Y.c.ρ * ᶜSqₗᵖ + @. Yₜ.c.ρq_ice += Y.c.ρ * ᶜSqᵢᵖ @. Yₜ.c.ρq_rai += Y.c.ρ * ᶜSqᵣᵖ @. Yₜ.c.ρq_sno += Y.c.ρ * ᶜSqₛᵖ - return nothing end function precipitation_tendency!( @@ -357,15 +387,16 @@ function precipitation_tendency!( Y, p, t, + ::NonEquilMoistModel, precip_model::Microphysics1Moment, turbconv_model::DiagnosticEDMFX, ) # Source terms from EDMFX environment - (; ᶜSeₜᵖ⁰, ᶜSqₜᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰) = p.precomputed + (; ᶜSqₗᵖ⁰, ᶜSqᵢᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰) = p.precomputed # Source terms from EDMFX updrafts - (; ᶜSeₜᵖʲs, ᶜSqₜᵖʲs, ᶜSqᵣᵖʲs, ᶜSqₛᵖʲs) = p.precomputed + (; ᶜSqₗᵖʲs, ᶜSqᵢᵖʲs, ᶜSqᵣᵖʲs, ᶜSqₛᵖʲs) = p.precomputed # Grid mean precipitation sinks - (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation + (; ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precipitation (; ᶜρaʲs) = p.precomputed @@ -375,18 +406,16 @@ function precipitation_tendency!( # Update from environment precipitation sources # and the grid mean precipitation sinks - @. Yₜ.c.ρ += Y.c.ρ * (ᶜSqₜᵖ⁰ + ᶜSqₜᵖ) - @. Yₜ.c.ρq_tot += Y.c.ρ * (ᶜSqₜᵖ⁰ + ᶜSqₜᵖ) - @. Yₜ.c.ρe_tot += Y.c.ρ * (ᶜSeₜᵖ⁰ + ᶜSeₜᵖ) + @. Yₜ.c.ρq_liq += Y.c.ρ * (ᶜSqₗᵖ⁰ + ᶜSqₗᵖ) + @. Yₜ.c.ρq_ice += Y.c.ρ * (ᶜSqᵢᵖ⁰ + ᶜSqᵢᵖ) @. Yₜ.c.ρq_rai += Y.c.ρ * (ᶜSqᵣᵖ⁰ + ᶜSqᵣᵖ) @. Yₜ.c.ρq_sno += Y.c.ρ * (ᶜSqₛᵖ⁰ + ᶜSqₛᵖ) # Update from the updraft precipitation sources n = n_mass_flux_subdomains(p.atmos.turbconv_model) for j in 1:n - @. Yₜ.c.ρ += ᶜρaʲs.:($$j) * ᶜSqₜᵖʲs.:($$j) - @. Yₜ.c.ρq_tot += ᶜρaʲs.:($$j) * ᶜSqₜᵖʲs.:($$j) - @. Yₜ.c.ρe_tot += ᶜρaʲs.:($$j) * ᶜSeₜᵖʲs.:($$j) + @. Yₜ.c.ρq_liq += ᶜρaʲs.:($$j) * ᶜSqₗᵖʲs.:($$j) + @. Yₜ.c.ρq_ice += ᶜρaʲs.:($$j) * ᶜSqᵢᵖʲs.:($$j) @. Yₜ.c.ρq_rai += ᶜρaʲs.:($$j) * ᶜSqᵣᵖʲs.:($$j) @. Yₜ.c.ρq_sno += ᶜρaʲs.:($$j) * ᶜSqₛᵖʲs.:($$j) end @@ -396,15 +425,16 @@ function precipitation_tendency!( Y, p, t, + ::NonEquilMoistModel, precip_model::Microphysics1Moment, turbconv_model::PrognosticEDMFX, ) # Source terms from EDMFX environment - (; ᶜSeₜᵖ⁰, ᶜSqₜᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰, ᶜρa⁰) = p.precomputed + (; ᶜSqₗᵖ⁰, ᶜSqᵢᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰, ᶜρa⁰) = p.precomputed # Source terms from EDMFX updrafts - (; ᶜSeₜᵖʲs, ᶜSqₜᵖʲs, ᶜSqᵣᵖʲs, ᶜSqₛᵖʲs) = p.precomputed + (; ᶜSqₗᵖʲs, ᶜSqᵢᵖʲs, ᶜSqᵣᵖʲs, ᶜSqₛᵖʲs) = p.precomputed # Grid mean precipitation sinks - (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation + (; ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precipitation # Populate the cache and precipitation surface fluxes compute_precipitation_cache!(Y, p, precip_model, turbconv_model) @@ -412,18 +442,16 @@ function precipitation_tendency!( # Update from environment precipitation sources # and the grid mean precipitation sinks - @. Yₜ.c.ρ += ᶜρa⁰ * ᶜSqₜᵖ⁰ + Y.c.ρ * ᶜSqₜᵖ - @. Yₜ.c.ρq_tot += ᶜρa⁰ * ᶜSqₜᵖ⁰ + Y.c.ρ * ᶜSqₜᵖ - @. Yₜ.c.ρe_tot += ᶜρa⁰ * ᶜSeₜᵖ⁰ + Y.c.ρ * ᶜSeₜᵖ + @. Yₜ.c.ρq_liq += ᶜρa⁰ * ᶜSqₗᵖ⁰ + Y.c.ρ * ᶜSqₗᵖ + @. Yₜ.c.ρq_ice += ᶜρa⁰ * ᶜSqᵢᵖ⁰ + Y.c.ρ * ᶜSqᵢᵖ @. Yₜ.c.ρq_rai += ᶜρa⁰ * ᶜSqᵣᵖ⁰ + Y.c.ρ * ᶜSqᵣᵖ @. Yₜ.c.ρq_sno += ᶜρa⁰ * ᶜSqₛᵖ⁰ + Y.c.ρ * ᶜSqₛᵖ # Update from the updraft precipitation sources n = n_mass_flux_subdomains(p.atmos.turbconv_model) for j in 1:n - @. Yₜ.c.ρ += Y.c.sgsʲs.:($$j).ρa * ᶜSqₜᵖʲs.:($$j) - @. Yₜ.c.ρq_tot += Y.c.sgsʲs.:($$j).ρa * ᶜSqₜᵖʲs.:($$j) - @. Yₜ.c.ρe_tot += Y.c.sgsʲs.:($$j).ρa * ᶜSeₜᵖʲs.:($$j) + @. Yₜ.c.ρq_liq += Y.c.sgsʲs.:($$j).ρa * ᶜSqₗᵖʲs.:($$j) + @. Yₜ.c.ρq_ice += Y.c.sgsʲs.:($$j).ρa * ᶜSqᵢᵖʲs.:($$j) @. Yₜ.c.ρq_rai += Y.c.sgsʲs.:($$j).ρa * ᶜSqᵣᵖʲs.:($$j) @. Yₜ.c.ρq_sno += Y.c.sgsʲs.:($$j).ρa * ᶜSqₛᵖʲs.:($$j) end diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl index 35e056b0b7..c36fd6b070 100644 --- a/src/parameters/create_parameters.jl +++ b/src/parameters/create_parameters.jl @@ -110,6 +110,7 @@ atmos_name_map = (; cloud_parameters(FT_or_toml) = (; liquid = CM.Parameters.CloudLiquid(FT_or_toml), ice = CM.Parameters.CloudIce(FT_or_toml), + Ch2022 = CM.Parameters.Chen2022VelType(FT_or_toml), ) microphys_1m_parameters(::Type{FT}) where {FT <: AbstractFloat} = diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index 44dcbf6ba5..51a68583d1 100644 --- a/src/prognostic_equations/advection.jl +++ b/src/prognostic_equations/advection.jl @@ -180,8 +180,9 @@ function edmfx_sgs_vertical_advection_tendency!( n = n_prognostic_mass_flux_subdomains(turbconv_model) (; dt) = p ᶜJ = Fields.local_geometry_field(Y.c).J + ᶜlg = Fields.local_geometry_field(Y.c) (; edmfx_upwinding) = p.atmos.numerics - (; ᶠu³ʲs, ᶠKᵥʲs, ᶜρʲs) = p.precomputed + (; ᶠu³ʲs, ᶠKᵥʲs, ᶜρʲs, ᶜwₕʲs, ᶜwₜʲs) = p.precomputed (; ᶠgradᵥ_ᶜΦ) = p.core ᶠz = Fields.coordinate_field(Y.f).z @@ -220,19 +221,44 @@ function edmfx_sgs_vertical_advection_tendency!( dt, edmfx_upwinding, ) - vertical_advection!( Yₜ.c.sgsʲs.:($j).mse, - ᶠu³ʲs.:($j), + ᶠu³ʲs.:($j) .- ᶠright_bias.(CT3.(ᶜwₕʲs.:($j))), Y.c.sgsʲs.:($j).mse, edmfx_upwinding, ) - vertical_advection!( Yₜ.c.sgsʲs.:($j).q_tot, - ᶠu³ʲs.:($j), + ᶠu³ʲs.:($j) .- ᶠright_bias.(CT3.(ᶜwₜʲs.:($j))), Y.c.sgsʲs.:($j).q_tot, edmfx_upwinding, ) + # TODO boundary condition for right bias + if p.atmos.moisture_model isa NonEquilMoistModel + vertical_advection!( + Yₜ.c.sgsʲs.:($j).q_liq, + ᶠu³ʲs.:($j) .- ᶠright_bias.(CT3.(ᶜwₗʲs.:($j))), + Y.c.sgsʲs.:($j).q_liq, + edmfx_upwinding, + ) + vertical_advection!( + Yₜ.c.sgsʲs.:($j).q_ice, + ᶠu³ʲs.:($j) .- ᶠright_bias.(CT3.(ᶜwᵢʲs.:($j))), + Y.c.sgsʲs.:($j).q_ice, + edmfx_upwinding, + ) + vertical_advection!( + Yₜ.c.sgsʲs.:($j).q_rai, + ᶠu³ʲs.:($j) .- ᶠright_bias.(CT3.(ᶜwᵣʲs.:($j))), + Y.c.sgsʲs.:($j).q_rai, + edmfx_upwinding, + ) + vertical_advection!( + Yₜ.c.sgsʲs.:($j).q_sno, + ᶠu³ʲs.:($j) .- ᶠright_bias.(CT3.(ᶜwₛʲs.:($j))), + Y.c.sgsʲs.:($j).q_sno, + edmfx_upwinding, + ) + end end end diff --git a/src/prognostic_equations/edmfx_precipitation.jl b/src/prognostic_equations/edmfx_precipitation.jl index 52a697d161..ca6852c5dd 100644 --- a/src/prognostic_equations/edmfx_precipitation.jl +++ b/src/prognostic_equations/edmfx_precipitation.jl @@ -36,30 +36,3 @@ function edmfx_precipitation_tendency!( end return nothing end - -function edmfx_precipitation_tendency!( - Yₜ, - Y, - p, - t, - turbconv_model::PrognosticEDMFX, - precip_model::Microphysics1Moment, -) - n = n_mass_flux_subdomains(turbconv_model) - (; ᶜSeₜᵖʲs, ᶜSqₜᵖʲs, ᶜtsʲs) = p.precomputed - thp = CAP.thermodynamics_params(p.params) - (; ᶜΦ) = p.core - - for j in 1:n - - @. Yₜ.c.sgsʲs.:($$j).ρa += Y.c.sgsʲs.:($$j).ρa * ᶜSqₜᵖʲs.:($$j) - - @. Yₜ.c.sgsʲs.:($$j).mse += - ᶜSeₜᵖʲs.:($$j) - - ᶜSqₜᵖʲs.:($$j) * TD.internal_energy(thp, ᶜtsʲs.:($$j)) - - @. Yₜ.c.sgsʲs.:($$j).q_tot += - ᶜSqₜᵖʲs.:($$j) * (1 - Y.c.sgsʲs.:($$j).q_tot) - end - return nothing -end diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl index 908671a42b..207c0d87eb 100644 --- a/src/prognostic_equations/remaining_tendency.jl +++ b/src/prognostic_equations/remaining_tendency.jl @@ -61,7 +61,12 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t) edmfx_filter_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model) edmfx_tke_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model) # Non-equilibrium cloud formation - cloud_condensate_tendency!(Yₜ, p, p.atmos.moisture_model) + cloud_condensate_tendency!( + Yₜ, + p, + p.atmos.moisture_model, + p.atmos.precip_model, + ) edmfx_precipitation_tendency!( Yₜ, Y, @@ -75,6 +80,7 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t) Y, p, t, + p.atmos.moisture_model, p.atmos.precip_model, p.atmos.turbconv_model, ) diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index 675738a696..42398041aa 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -576,8 +576,8 @@ function args_integrator(parsed_args, Y, p, tspan, ode_algo, callback) # Can we just pass implicit_tendency! and jac_prototype etc.? lim! = limiters_func!, dss!, - post_explicit! = set_precomputed_quantities!, - post_implicit! = set_precomputed_quantities!, + cache! = set_precomputed_quantities!, + cache_imp! = set_precomputed_quantities!, ) else SciMLBase.SplitFunction(implicit_func, remaining_tendency!) diff --git a/src/utils/variable_manipulations.jl b/src/utils/variable_manipulations.jl index 721bc78b0e..f354f922fb 100644 --- a/src/utils/variable_manipulations.jl +++ b/src/utils/variable_manipulations.jl @@ -211,6 +211,38 @@ mass-flux subdomain states are stored in `sgsʲs`. """ ρaq_tot⁺(sgsʲs) = mapreduce_with_init(sgsʲ -> sgsʲ.ρa * sgsʲ.q_tot, +, sgsʲs) +""" + ρaq_liq⁺(sgsʲs) + +Computes the total mass-flux subdomain area-weighted ρq_liq, assuming that the +mass-flux subdomain states are stored in `sgsʲs`. +""" +ρaq_liq⁺(sgsʲs) = mapreduce_with_init(sgsʲ -> sgsʲ.ρa * sgsʲ.q_liq, +, sgsʲs) + +""" + ρaq_ice⁺(sgsʲs) + +Computes the total mass-flux subdomain area-weighted ρq_ice, assuming that the +mass-flux subdomain states are stored in `sgsʲs`. +""" +ρaq_ice⁺(sgsʲs) = mapreduce_with_init(sgsʲ -> sgsʲ.ρa * sgsʲ.q_ice, +, sgsʲs) + +""" + ρaq_rai⁺(sgsʲs) + +Computes the total mass-flux subdomain area-weighted ρq_rai, assuming that the +mass-flux subdomain states are stored in `sgsʲs`. +""" +ρaq_rai⁺(sgsʲs) = mapreduce_with_init(sgsʲ -> sgsʲ.ρa * sgsʲ.q_rai, +, sgsʲs) + +""" + ρaq_sno⁺(sgsʲs) + +Computes the total mass-flux subdomain area-weighted ρq_sno, assuming that the +mass-flux subdomain states are stored in `sgsʲs`. +""" +ρaq_sno⁺(sgsʲs) = mapreduce_with_init(sgsʲ -> sgsʲ.ρa * sgsʲ.q_sno, +, sgsʲs) + """ ρa⁰(gs) @@ -220,7 +252,7 @@ subdomain states are stored in `gs.sgsʲs`. ρa⁰(gs) = gs.ρ - mapreduce_with_init(sgsʲ -> sgsʲ.ρa, +, gs.sgsʲs) """ - u₃⁺(ρaʲs, u₃ʲs, ρ, u₃, turbconv_model) + u₃⁺(ρaʲs, u₃ʲs, ρ, u₃, turbconv_model) Computes the average mass-flux subdomain vertical velocity `u₃⁺` by dividing the total momentum `ρaw⁺` by the total area-weighted density `ρa⁺`, both of which diff --git a/test/parameterized_tendencies/microphysics/precipitation.jl b/test/parameterized_tendencies/microphysics/precipitation.jl index 2d6a6fc161..57b25baaf9 100644 --- a/test/parameterized_tendencies/microphysics/precipitation.jl +++ b/test/parameterized_tendencies/microphysics/precipitation.jl @@ -9,33 +9,34 @@ import CloudMicrophysics as CM include("../../test_helpers.jl") -### Common Objects ### -@testset begin - "Precipitation tendency functions" +import Test + +@testset "Equilibrium Moisture + 0-moment precipitation RHS terms" begin + ### Boilerplate default integrator objects config = CA.AtmosConfig( Dict( - "initial_condition" => "PrecipitatingColumn", - "moist" => "nonequil", + "initial_condition" => "DYCOMS_RF02", + "moist" => "equil", "precip_model" => "0M", "config" => "column", "output_default_diagnostics" => false, ), - job_id = "precipitation1", + job_id = "equil_0M", ) (; Y, p, params) = generate_test_simulation(config) FT = eltype(Y) ᶜYₜ = zero(Y) - ### Component test begins here - @info "0M Scheme" - ### 0-Moment Scheme - precip_model = CA.Microphysics0Moment() - precip_cache = CA.precipitation_cache(Y, precip_model) + # Set all model choices + (; turbconv_model, moisture_model, precip_model) = p.atmos + # Test cache to verify expected variables exist in tendency function + precip_cache = CA.precipitation_cache(Y, precip_model) test_varnames = ( :ᶜS_ρq_tot, + :ᶜS_ρe_tot, :ᶜ3d_rain, :ᶜ3d_snow, :surface_rain_flux, @@ -44,25 +45,38 @@ include("../../test_helpers.jl") for var_name in test_varnames @test var_name ∈ propertynames(precip_cache) end - turbconv_model = nothing # Extend to other turbulence convection models + + # No NaNs in cache CA.compute_precipitation_cache!(Y, p, precip_model, turbconv_model) @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. - CA.precipitation_tendency!(ᶜYₜ, Y, p, FT(0), precip_model, turbconv_model) + CA.precipitation_tendency!( + ᶜYₜ, + Y, + p, + FT(0), + moisture_model, + precip_model, + turbconv_model, + ) @test ᶜYₜ.c.ρ == ᶜYₜ.c.ρq_tot @test ᶜYₜ.c.ρ == p.precipitation.ᶜS_ρq_tot - # test nonequilibrium cloud condensate - moist_model = CA.NonEquilMoistModel() - CA.cloud_condensate_tendency!(ᶜYₜ, p, moist_model) - @assert !any(isnan, ᶜYₜ.c.ρq_liq) - @assert !any(isnan, ᶜYₜ.c.ρq_ice) + # No cloud condensation tendency for the equilibrium model + @test CA.cloud_condensate_tendency!( + ᶜYₜ, + p, + moisture_model, + precip_model, + ) isa Nothing +end - ### 1-Moment Scheme - @info "1M Scheme" +@testset "NonEquilibrium Moisture + 1-moment precipitation RHS terms" begin + + ### Boilerplate default integrator objects config = CA.AtmosConfig( Dict( "initial_condition" => "PrecipitatingColumn", @@ -74,51 +88,60 @@ include("../../test_helpers.jl") job_id = "precipitation2", ) (; Y, p, params) = generate_test_simulation(config) - precip_model = CA.Microphysics1Moment() - (; turbconv_model) = p.atmos + + FT = eltype(Y) + ᶜYₜ = zero(Y) + + # Set all model choices + (; turbconv_model, moisture_model, precip_model) = p.atmos + + # Test cache to verify expected variables exist in tendency function precip_cache = CA.precipitation_cache(Y, precip_model) - ᶜYₜ = Y .* FT(0) test_varnames = - (:ᶜSqₜᵖ, :ᶜSqᵣᵖ, :ᶜSqₛᵖ, :ᶜSeₜᵖ, :surface_rain_flux, :surface_snow_flux) + (:ᶜSqₗᵖ, :ᶜSqᵢᵖ, :ᶜSqᵣᵖ, :ᶜSqₛᵖ, :surface_rain_flux, :surface_snow_flux) for var_name in test_varnames @test var_name ∈ propertynames(precip_cache) end # test helper functions - @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), 5) == FT(1) # compute source terms based on the last model state - CA.precipitation_tendency!(ᶜYₜ, Y, p, FT(0), precip_model, turbconv_model) + CA.precipitation_tendency!( + ᶜYₜ, + Y, + p, + FT(0), + moisture_model, + precip_model, + turbconv_model, + ) # check for nans @assert !any(isnan, ᶜYₜ.c.ρ) @assert !any(isnan, ᶜYₜ.c.ρq_tot) @assert !any(isnan, ᶜYₜ.c.ρe_tot) + @assert !any(isnan, ᶜYₜ.c.ρq_liq) + @assert !any(isnan, ᶜYₜ.c.ρq_ice) @assert !any(isnan, ᶜYₜ.c.ρq_rai) @assert !any(isnan, ᶜYₜ.c.ρq_sno) + @assert !any(isnan, p.precomputed.ᶜwₗ) + @assert !any(isnan, p.precomputed.ᶜwᵢ) @assert !any(isnan, p.precomputed.ᶜwᵣ) @assert !any(isnan, p.precomputed.ᶜwₛ) # test water budget @test ᶜYₜ.c.ρ == ᶜYₜ.c.ρq_tot - @test ᶜYₜ.c.ρ == Y.c.ρ .* p.precipitation.ᶜSqₜᵖ - @test all( - isapprox( - .-p.precipitation.ᶜSqₛᵖ .- p.precipitation.ᶜSqᵣᵖ, - p.precipitation.ᶜSqᵣᵖ, - atol = eps(FT), - ), - ) + @assert iszero(ᶜYₜ.c.ρ) # test nonequilibrium cloud condensate - moist_model = CA.NonEquilMoistModel() - CA.cloud_condensate_tendency!(ᶜYₜ, p, moist_model) + CA.cloud_condensate_tendency!(ᶜYₜ, p, moisture_model, precip_model) @assert !any(isnan, ᶜYₜ.c.ρq_liq) @assert !any(isnan, ᶜYₜ.c.ρq_ice) # test if terminal velocity is positive + @test minimum(p.precomputed.ᶜwₗ) >= FT(0) + @test minimum(p.precomputed.ᶜwᵢ) >= FT(0) @test minimum(p.precomputed.ᶜwᵣ) >= FT(0) @test minimum(p.precomputed.ᶜwₛ) >= FT(0) diff --git a/test/restart.jl b/test/restart.jl index 04b500a80f..4d8c0adca9 100644 --- a/test/restart.jl +++ b/test/restart.jl @@ -351,7 +351,7 @@ if MANYTESTS for configuration in configurations if configuration == "sphere" - moistures = ["equil", "nonequil"] + moistures = ["nonequil"] precips = ["1M"] topography = "Earth" turbconv_models = [nothing, "diagnostic_edmfx"]