diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 07c5dde1ab..8a232cfee7 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -194,7 +194,7 @@ steps: - label: ":computer: hydrostatic balance (ρe) float64" command: > - julia --color=yes --project=examples examples/hybrid/driver.jl + julia --color=yes --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/sphere_hydrostatic_balance_rhoe_ft64.yml julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl @@ -449,7 +449,7 @@ steps: - label: ":computer: held suarez (ρe) equilmoist topography (dcmip)" command: > - julia --color=yes --project=examples examples/hybrid/driver.jl + julia --color=yes --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/sphere_held_suarez_rhoe_equilmoist_topography_dcmip.yml julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl @@ -464,8 +464,8 @@ steps: julia --color=yes --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200.yml - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl - --data_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200 + julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl + --data_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200 --out_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200 julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl @@ -478,14 +478,14 @@ steps: command: > julia --color=yes --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth.yml - + julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl - --data_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth + --data_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth --out_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth - - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl - --nc_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth - --fig_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth + + julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl + --nc_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth + --fig_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth --case_name elevation_spectrum artifact_paths: "sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth/*" @@ -616,7 +616,7 @@ steps: artifact_paths: "diagnostic_edmfx_test_box/*" agents: slurm_mem: 20GB - + - label: ":genie: Diagnostic EDMFX GABLS in a box" command: > julia --color=yes --project=examples examples/hybrid/driver.jl @@ -632,7 +632,7 @@ steps: artifact_paths: "diagnostic_edmfx_bomex_box/*" agents: slurm_mem: 20GB - + - label: ":genie: Diagnostic EDMFX Bomex stretched grid in a box" command: > julia --color=yes --project=examples examples/hybrid/driver.jl @@ -656,7 +656,7 @@ steps: artifact_paths: "diagnostic_edmfx_rico_box/*" agents: slurm_mem: 20GB - + - label: ":genie: Diagnostic EDMFX TRMM in a box" command: > julia --color=yes --project=examples examples/hybrid/driver.jl @@ -744,7 +744,7 @@ steps: command: > julia --color=yes --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/edmf_nieuwstadt_jfnk_imex.yml - + julia --color=yes --project=examples regression_tests/test_mse.jl --job_id edmf_nieuwstadt_jfnk_imex --out_dir edmf_nieuwstadt_jfnk_imex artifact_paths: "edmf_nieuwstadt_jfnk_imex/*" @@ -863,7 +863,7 @@ steps: key: "gpu_baroclinic_wave_rhoe" command: > julia --project -e 'using CUDA; CUDA.versioninfo()' - + julia --color=yes --project=examples examples/hybrid/driver.jl --config_file $PERF_CONFIG_PATH/gpu_baroclinic_wave_rhoe.yml artifact_paths: "gpu_implicit_barowave_ref/*" @@ -884,7 +884,7 @@ steps: key: "gpu_held_suarez_rhoe_hightop" command: > julia --project -e 'using CUDA; CUDA.versioninfo()' - + julia --color=yes --project=examples examples/hybrid/driver.jl --config_file $PERF_CONFIG_PATH/gpu_held_suarez_rhoe_hightop.yml artifact_paths: "gpu_held_suarez_rhoe_hightop/*" diff --git a/config/model_configs/diagnostic_edmfx_trmm_box.yml b/config/model_configs/diagnostic_edmfx_trmm_box.yml index 57cb3f5427..e3f7497153 100644 --- a/config/model_configs/diagnostic_edmfx_trmm_box.yml +++ b/config/model_configs/diagnostic_edmfx_trmm_box.yml @@ -1,28 +1,28 @@ job_id: diagnostic_edmfx_trmm_box -initial_condition: TRMM_LBA -rad: TRMM_LBA +initial_condition: TRMM_LBA +rad: TRMM_LBA surface_setup: TRMM_LBA -turbconv: diagnostic_edmfx +turbconv: diagnostic_edmfx prognostic_tke: true -edmfx_upwinding: first_order -edmfx_entr_detr: true -edmfx_nh_pressure: true +edmfx_upwinding: first_order +edmfx_entr_detr: true +edmfx_nh_pressure: true edmfx_sgs_flux: true moist: equil apply_limiter: false precip_model: "0M" -config: box +config: box hyperdiff: "true" kappa_4: 1e12 -x_max: 1e5 -y_max: 1e5 -x_elem: 2 -y_elem: 2 -z_elem: 82 -z_max: 16400 +x_max: 1e5 +y_max: 1e5 +x_elem: 2 +y_elem: 2 +z_elem: 82 +z_max: 16400 z_stretch: false -dt: 10secs -t_end: 6hours -dt_save_to_disk: 10mins -FLOAT_TYPE: "Float64" \ No newline at end of file +dt: 10secs +t_end: 6hours +dt_save_to_disk: 10mins +FLOAT_TYPE: "Float64" diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl index 58095ab0de..d717afa6a4 100644 --- a/src/ClimaAtmos.jl +++ b/src/ClimaAtmos.jl @@ -88,6 +88,13 @@ include(joinpath("prognostic_equations", "cloud_fraction.jl")) include( joinpath("parameterized_tendencies", "microphysics", "precipitation.jl"), ) +include( + joinpath( + "parameterized_tendencies", + "microphysics", + "microphysics_wrappers.jl", + ), +) include( joinpath("prognostic_equations", "vertical_diffusion_boundary_layer.jl"), ) diff --git a/src/TurbulenceConvection_deprecated/closures/microphysics.jl b/src/TurbulenceConvection_deprecated/closures/microphysics.jl index 63d4b9c28e..702aed6537 100644 --- a/src/TurbulenceConvection_deprecated/closures/microphysics.jl +++ b/src/TurbulenceConvection_deprecated/closures/microphysics.jl @@ -54,7 +54,9 @@ function precipitation_formation( qi_tendency += S_qt * (1 - λ) θ_liq_ice_tendency -= S_qt / Π_m / c_pm * (L_v0 * λ + L_s0 * (1 - λ)) - e_tot_tendency += (λ * I_l + (1 - λ) * I_i + Φ) * S_qt + e_tot_tendency += + S_qt * + e_tot_0M_precipitation_sources_helper(thermo_params, ts, Φ) end if precip_model isa Microphysics1Moment diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl index b174d814b8..a636ac906a 100644 --- a/src/cache/diagnostic_edmf_precomputed_quantities.jl +++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl @@ -83,7 +83,7 @@ end Updates the precomputed quantities stored in `p` for diagnostic edmfx. """ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) - (; moisture_model, turbconv_model) = p.atmos + (; moisture_model, turbconv_model, precip_model) = p.atmos FT = eltype(Y) n = n_mass_flux_subdomains(turbconv_model) ᶜz = Fields.coordinate_field(Y.c).z @@ -105,10 +105,13 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) ᶜρʲs, ᶜentr_detrʲs, ᶜnh_pressureʲs, + ᶜS_q_totʲs, ) = p - (; ᶠu³⁰, ᶜu⁰, ᶜtke⁰, ᶜlinear_buoygrad, ᶜshear², ᶜmixing_length) = p + (; ᶠu³⁰, ᶜu⁰, ᶜtke⁰, ᶜlinear_buoygrad, ᶜshear², ᶜmixing_length, ᶜS_q_tot⁰) = + p (; ᶜK_h, ᶜK_u, ρatke_flux) = p thermo_params = CAP.thermodynamics_params(params) + microphys_params = CAP.microphysics_params(params) ᶜlg = Fields.local_geometry_field(Y.c) @. ᶜtke⁰ = Y.c.sgs⁰.ρatke / Y.c.ρ @@ -154,6 +157,7 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) ᶜtsʲ = ᶜtsʲs.:($j) ᶜρʲ = ᶜρʲs.:($j) ᶜq_totʲ = ᶜq_totʲs.:($j) + ᶜS_q_totʲ = ᶜS_q_totʲs.:($j) ρaʲ_int_level = Fields.field_values(Fields.level(ᶜρaʲ, 1)) u³ʲ_int_halflevel = Fields.field_values(Fields.level(ᶠu³ʲ, half)) @@ -237,6 +241,7 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) end_index = fieldcount(eltype(∂x∂ξ_level)) # This will be 4 in 2D and 9 in 3D. ∂x³∂ξ³_level = ∂x∂ξ_level.:($end_index) + Φ_prev_level = Fields.field_values(Fields.level(ᶜΦ, i - 1)) ρ_ref_prev_level = Fields.field_values(Fields.level(ᶜρ_ref, i - 1)) ∇Φ³_prev_level = Fields.field_values(Fields.level(ᶜ∇Φ³, i - 1)) ∇Φ³_data_prev_level = ∇Φ³_prev_level.components.data.:1 @@ -249,6 +254,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) ts_prev_level = Fields.field_values(Fields.level(ᶜts, i - 1)) p_prev_level = Fields.field_values(Fields.level(ᶜp, i - 1)) z_prev_level = Fields.field_values(Fields.level(ᶜz, i - 1)) + S_q_tot⁰_prev_level = + Fields.field_values(Fields.level(ᶜS_q_tot⁰, i - 1)) local_geometry_prev_level = Fields.field_values( Fields.level(Fields.local_geometry_field(Y.c), i - 1), @@ -267,6 +274,7 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) ᶜq_totʲ = ᶜq_totʲs.:($j) ᶜentr_detrʲ = ᶜentr_detrʲs.:($j) ᶜnh_pressureʲ = ᶜnh_pressureʲs.:($j) + ᶜS_q_totʲ = ᶜS_q_totʲs.:($j) ρaʲ_level = Fields.field_values(Fields.level(ᶜρaʲ, i)) u³ʲ_halflevel = Fields.field_values(Fields.level(ᶠu³ʲ, i - half)) @@ -292,6 +300,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) Fields.field_values(Fields.level(ᶜnh_pressureʲ, i - 1)) scale_height = CAP.R_d(params) * CAP.T_surf_ref(params) / CAP.grav(params) + S_q_totʲ_prev_level = + Fields.field_values(Fields.level(ᶜS_q_totʲ, i - 1)) @. entr_detrʲ_prev_level = pi_groups_entr_detr( params, @@ -330,6 +340,18 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) nh_pressureʲ_data_prev_level = nh_pressureʲ_prev_level.components.data.:1 + # Updraft q_tot sources from precipitation formation + # To be applied in updraft continuity, moisture and energy + # for updrafts and grid mean + @. S_q_totʲ_prev_level = q_tot_precipitation_sources( + precip_model, + thermo_params, + microphys_params, + dt, + q_totʲ_prev_level, + tsʲ_prev_level, + ) + @. ρaʲu³ʲ_data = (1 / local_geometry_halflevel.J) * ( local_geometry_prev_halflevel.J * @@ -341,7 +363,10 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) (1 / local_geometry_halflevel.J) * ( local_geometry_prev_level.J * ρaʲ_prev_level * - (entr_detrʲ_prev_level.entr - entr_detrʲ_prev_level.detr) + ( + entr_detrʲ_prev_level.entr - + entr_detrʲ_prev_level.detr + S_q_totʲ_prev_level + ) ) # Using constant exponents in broadcasts allocate, so we use @@ -432,7 +457,13 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) ρaʲ_prev_level * ( entr_detrʲ_prev_level.entr * h_tot_prev_level - - entr_detrʲ_prev_level.detr * h_totʲ_prev_level + entr_detrʲ_prev_level.detr * h_totʲ_prev_level + + S_q_totʲ_prev_level * + e_tot_0M_precipitation_sources_helper( + thermo_params, + tsʲ_prev_level, + Φ_prev_level, + ) ) ) @. h_totʲ_level = ifelse( @@ -459,7 +490,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) ρaʲ_prev_level * ( entr_detrʲ_prev_level.entr * q_tot_prev_level - - entr_detrʲ_prev_level.detr * q_totʲ_prev_level + entr_detrʲ_prev_level.detr * q_totʲ_prev_level + + S_q_totʲ_prev_level ) ) @. q_totʲ_level = ifelse( @@ -489,6 +521,16 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) ) end + # Environment precipitation sources (to be applied to grid mean) + @. S_q_tot⁰_prev_level = q_tot_precipitation_sources( + precip_model, + thermo_params, + microphys_params, + dt, + q_tot_prev_level, + ts_prev_level, + ) + ρaʲs_level = Fields.field_values(Fields.level(ᶜρaʲs, i)) u³ʲs_halflevel = Fields.field_values(Fields.level(ᶠu³ʲs, i - half)) u³⁰_halflevel = Fields.field_values(Fields.level(ᶠu³⁰, i - half)) @@ -506,18 +548,48 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t) i_top = Spaces.nlevels(axes(Y.c)) u³⁰_halflevel = Fields.field_values(Fields.level(ᶠu³⁰, i_top + half)) @. u³⁰_halflevel = CT3(0) + for j in 1:n ᶠu³ʲ = ᶠu³ʲs.:($j) ᶜuʲ = ᶜuʲs.:($j) ᶠu³ʲ = ᶠu³ʲs.:($j) ᶜentr_detrʲ = ᶜentr_detrʲs.:($j) + ᶜS_q_totʲ = ᶜS_q_totʲs.:($j) + ᶜq_totʲ = ᶜq_totʲs.:($j) + ᶜtsʲ = ᶜtsʲs.:($j) + u³ʲ_halflevel = Fields.field_values(Fields.level(ᶠu³ʲ, i_top + half)) @. u³ʲ_halflevel = CT3(0) + entr_detrʲ_level = Fields.field_values(Fields.level(ᶜentr_detrʲ, i_top)) fill!(entr_detrʲ_level, RecursiveApply.rzero(eltype(entr_detrʲ_level))) @. ᶜuʲ = C123(Y.c.uₕ) + ᶜinterp(C123(ᶠu³ʲ)) + + S_q_totʲ_level = Fields.field_values(Fields.level(ᶜS_q_totʲ, i_top)) + q_totʲ_level = Fields.field_values(Fields.level(ᶜq_totʲ, i_top)) + tsʲ_level = Fields.field_values(Fields.level(ᶜtsʲ, i_top)) + @. S_q_totʲ_level = q_tot_precipitation_sources( + precip_model, + thermo_params, + microphys_params, + dt, + q_totʲ_level, + tsʲ_level, + ) end + q_tot_level = Fields.field_values(Fields.level(q_tot, i_top)) + ts_level = Fields.field_values(Fields.level(ᶜts, i_top)) + S_q_tot⁰_level = Fields.field_values(Fields.level(ᶜS_q_tot⁰, i_top)) + @. S_q_tot⁰_level = q_tot_precipitation_sources( + precip_model, + thermo_params, + microphys_params, + dt, + q_tot_level, + ts_level, + ) + @. ᶜu⁰ = C123(Y.c.uₕ) + ᶜinterp(C123(ᶠu³⁰)) @. ᶜlinear_buoygrad = buoyancy_gradients( diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index 9a4d2668d0..86d1981a89 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -106,6 +106,8 @@ function precomputed_quantities(Y, atmos) NTuple{n, NamedTuple{(:entr, :detr), NTuple{2, FT}}}, ), ᶜnh_pressureʲs = similar(Y.c, NTuple{n, CT3{FT}}), + ᶜS_q_totʲs = similar(Y.c, NTuple{n, FT}), + ᶜS_q_tot⁰ = similar(Y.c, FT), ᶠu³⁰ = similar(Y.f, CT3{FT}), ᶜu⁰ = similar(Y.c, C123{FT}), ᶜtke⁰ = similar(Y.c, FT), diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl new file mode 100644 index 0000000000..db1c381ebb --- /dev/null +++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl @@ -0,0 +1,66 @@ +# A set of wrappers for using CloudMicrophysics.jl functions inside EDMFX loops + +import Thermodynamics as TD +import CloudMicrophysics.Microphysics0M as CM0 + +""" + q_tot_precipitation_sources(precip_model, t_prs, m_prs, dt, q_tot, ts) + + - precip_model - a type for precipitation scheme choice + - t_prs, m_prs - stricts with thermodynamic and microphysics parameters + - dt - model time step + - q_tot - total water specific humidity + - ts - thermodynamic state (see td package for details) + +Returns the q_tot source term due to precipitation formation +defined as Δm_tot / (m_dry + m_tot) +""" +function q_tot_precipitation_sources( + ::NoPrecipitation, + t_prs, + m_prs, + dt, + q_tot::FT, + ts, +) where {FT <: Real} + return FT(0) +end +function q_tot_precipitation_sources( + ::Microphysics0Moment, + t_prs, + m_prs, + dt, + q_tot::FT, + ts, +) where {FT <: Real} + return -min( + max(q_tot, FT(0)) / dt, + -CM0.remove_precipitation( + m_prs, + TD.PhasePartition(t_prs, ts), + TD.q_vap_saturation(t_prs, ts), + ), + ) +end + +""" + e_tot_0M_precipitation_sources_helper(t_prs, ts, Φ) + + - t_prs - set with thermodynamics parameters + - ts - thermodynamic state (see td package for details) + - Φ - geopotential + +Returns the total energy source term multiplier from precipitation formation +""" +function e_tot_0M_precipitation_sources_helper( + t_prs, + ts, + Φ::FT, +) where {FT <: Real} + + λ::FT = TD.liquid_fraction(t_prs, ts) + I_l::FT = TD.internal_energy_liquid(t_prs, ts) + I_i::FT = TD.internal_energy_ice(t_prs, ts) + + return λ * I_l + (1 - λ) * I_i + Φ +end diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl index f37d70cd70..f8e2579545 100644 --- a/src/parameterized_tendencies/microphysics/precipitation.jl +++ b/src/parameterized_tendencies/microphysics/precipitation.jl @@ -25,7 +25,6 @@ function precipitation_cache(Y, precip_model::Microphysics0Moment) return (; precip_model, ᶜS_ρq_tot = similar(Y.c, FT), - ᶜλ = similar(Y.c, FT), ᶜ3d_rain = similar(Y.c, FT), ᶜ3d_snow = similar(Y.c, FT), col_integrated_rain = similar(Fields.level(Y.c.ρ, 1), FT), @@ -63,6 +62,24 @@ function compute_precipitation_cache!( (qt_tendency_precip_formation_bulk + qt_tendency_precip_formation_en) end +function compute_precipitation_cache!( + Y, + p, + colidx, + ::Microphysics0Moment, + ::DiagnosticEDMFX, +) + # For environment we multiply by grid mean ρ and not byᶜρa⁰ + # I.e. assuming a⁰=1 + + (; ᶜS_ρq_tot, ᶜS_q_tot⁰, ᶜS_q_totʲs, ᶜρaʲs) = p + ρ = Y.c.ρ + FT = Spaces.undertype(axes(Y.c)) + + @. ᶜS_ρq_tot[colidx] = + sum(ᶜS_q_totʲs[colidx] * ᶜρaʲs[colidx]) + ᶜS_q_tot⁰[colidx] * ρ[colidx] +end + function precipitation_tendency!( Yₜ, Y, @@ -79,7 +96,6 @@ function precipitation_tendency!( ᶜ3d_rain, ᶜ3d_snow, ᶜS_ρq_tot, - ᶜλ, col_integrated_rain, col_integrated_snow, params, @@ -111,17 +127,12 @@ function precipitation_tendency!( @. col_integrated_snow[colidx] = col_integrated_snow[colidx] / CAP.ρ_cloud_liq(params) - # liquid fraction - @. ᶜλ[colidx] = TD.liquid_fraction(thermo_params, ᶜts[colidx]) - if :ρe_tot in propertynames(Y.c) @. Yₜ.c.ρe_tot[colidx] += - ᶜS_ρq_tot[colidx] * ( - ᶜλ[colidx] * - TD.internal_energy_liquid(thermo_params, ᶜts[colidx]) + - (1 - ᶜλ[colidx]) * - TD.internal_energy_ice(thermo_params, ᶜts[colidx]) + - ᶜΦ[colidx] + ᶜS_ρq_tot[colidx] * e_tot_0M_precipitation_sources_helper( + thermo_params, + ᶜts[colidx], + ᶜΦ[colidx], ) end return nothing