From b7d9443f2e81cfb358a3c7dfdf464cf927bbc667 Mon Sep 17 00:00:00 2001 From: Zhaoyi Shen <11598433+szy21@users.noreply.github.com> Date: Fri, 12 Apr 2024 17:06:25 -0700 Subject: [PATCH] add vertical fluctuation and nudging --- .buildkite/pipeline.yml | 8 +++ config/default_configs/default_config.yml | 6 ++ .../prognostic_edmfx_gcmdriven_column.yml | 40 ++++++++++++ post_processing/ci_plots.jl | 1 + src/ClimaAtmos.jl | 2 + src/cache/cache.jl | 8 +++ src/cache/temporary_quantities.jl | 1 + .../forcing/large_scale_advection.jl | 2 +- src/prognostic_equations/forcing/nudging.jl | 62 +++++++++++++++++++ .../forcing/vertical_fluctuation.jl | 51 +++++++++++++++ .../remaining_tendency.jl | 2 + src/solver/model_getters.jl | 38 +++++++++--- src/solver/type_getters.jl | 2 + src/solver/types.jl | 7 +++ test/coupler_compatibility.jl | 2 + 15 files changed, 224 insertions(+), 8 deletions(-) create mode 100644 config/model_configs/prognostic_edmfx_gcmdriven_column.yml create mode 100644 src/prognostic_equations/forcing/nudging.jl create mode 100644 src/prognostic_equations/forcing/vertical_fluctuation.jl diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 29ac02e5f6c..cddb32f68d7 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -629,6 +629,14 @@ steps: artifact_paths: "prognostic_edmfx_trmm_column_0M/*" agents: slurm_mem: 20GB + + - label: ":genie: Prognostic EDMFX GCM driven in a column" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/prognostic_edmfx_gcmdriven_column.yml + artifact_paths: "prognostic_edmfx_gcmdriven_column/output_active/*" + agents: + slurm_mem: 20GB - label: ":genie: Prognostic EDMFX Bomex in a box" command: > diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml index 176446dfb71..4c2981c0a49 100644 --- a/config/default_configs/default_config.yml +++ b/config/default_configs/default_config.yml @@ -229,6 +229,12 @@ check_precipitation: ls_adv: help: "Large-scale advection [`nothing` (default), `Bomex`, `LifeCycleTan2018`, `Rico`, `ARM_SGP`, `GATE_III`]" value: ~ +vert_fluc: + help: "Vertical fluctuation [`nothing` (default), `GCMDriven`]" + value: ~ +nudging: + help: "Nudging to a mean profile [`nothing` (default), `GCMDriven`]" + value: ~ fps: help: "Frames per second for animations" value: 5 diff --git a/config/model_configs/prognostic_edmfx_gcmdriven_column.yml b/config/model_configs/prognostic_edmfx_gcmdriven_column.yml new file mode 100644 index 00000000000..ca6c5ee3151 --- /dev/null +++ b/config/model_configs/prognostic_edmfx_gcmdriven_column.yml @@ -0,0 +1,40 @@ +job_id: "prognostic_edmfx_gcmdriven_column" +initial_condition: "Bomex" +subsidence: "Bomex" +edmf_coriolis: "Bomex" +ls_adv: "GCMDriven" +vert_fluc: "GCMDriven" +nudging: "GCMDriven" +surface_setup: "Bomex" +turbconv: "prognostic_edmfx" +implicit_diffusion: true +implicit_sgs_advection: true +approximate_linear_solve_iters: 2 +max_newton_iters_ode: 3 +edmfx_upwinding: first_order +edmfx_entr_model: "Generalized" +edmfx_detr_model: "Generalized" +edmfx_sgs_mass_flux: true +edmfx_sgs_diffusive_flux: true +edmfx_nh_pressure: true +edmfx_velocity_relaxation: true +prognostic_tke: true +moist: "equil" +config: "column" +z_max: 3e3 +z_elem: 60 +z_stretch: false +perturb_initstate: false +dt: "50secs" +t_end: "6hours" +dt_save_to_disk: "10mins" +toml: [toml/prognostic_edmfx_bomex.toml] +netcdf_output_at_levels: true +netcdf_interpolation_num_points: [2, 2, 60] +diagnostics: + - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl] + period: 10mins + - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke] + period: 10mins + - short_name: [entr, detr, lmix, bgrad, strain, edt, evu] + period: 10mins diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl index 2e03d5f0f17..322250793ab 100644 --- a/post_processing/ci_plots.jl +++ b/post_processing/ci_plots.jl @@ -945,6 +945,7 @@ EDMFBoxPlots = Union{ Val{:prognostic_edmfx_gabls_column}, Val{:prognostic_edmfx_bomex_fixtke_column}, Val{:prognostic_edmfx_bomex_column}, + Val{:prognostic_edmfx_gcmdriven_column}, Val{:prognostic_edmfx_bomex_stretched_column}, Val{:prognostic_edmfx_bomex_explicit_column}, Val{:prognostic_edmfx_dycoms_rf01_column}, diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl index 49a43c0b925..fda0d6f2cd7 100644 --- a/src/ClimaAtmos.jl +++ b/src/ClimaAtmos.jl @@ -50,6 +50,8 @@ include(joinpath("prognostic_equations", "implicit", "implicit_solver.jl")) include(joinpath("prognostic_equations", "remaining_tendency.jl")) include(joinpath("prognostic_equations", "forcing", "large_scale_advection.jl")) # TODO: should this be in tendencies/? +include(joinpath("prognostic_equations", "forcing", "vertical_fluctuation.jl")) +include(joinpath("prognostic_equations", "forcing", "nudging.jl")) include(joinpath("prognostic_equations", "forcing", "subsidence.jl")) include(joinpath("prognostic_equations", "surface_temp.jl")) diff --git a/src/cache/cache.jl b/src/cache/cache.jl index a31c75f0ef4..d7faac6a026 100644 --- a/src/cache/cache.jl +++ b/src/cache/cache.jl @@ -19,6 +19,8 @@ struct AtmosCache{ PR, SUB, LSAD, + VERTFLUC, + NUDGING, EDMFCOR, FOR, NONGW, @@ -81,6 +83,8 @@ struct AtmosCache{ precipitation::PR subsidence::SUB large_scale_advection::LSAD + vertical_fluctuation::VERTFLUC + nudging::NUDGING edmf_coriolis::EDMFCOR forcing::FOR non_orographic_gravity_wave::NONGW @@ -225,6 +229,8 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names) precipitation = precipitation_cache(Y, atmos) subsidence = subsidence_cache(Y, atmos) large_scale_advection = large_scale_advection_cache(Y, atmos) + vertical_fluctuation = vertical_fluctuation_cache(Y, atmos) + nudging = nudging_cache(Y, atmos) edmf_coriolis = edmf_coriolis_cache(Y, atmos) forcing = forcing_cache(Y, atmos) non_orographic_gravity_wave = non_orographic_gravity_wave_cache(Y, atmos) @@ -253,6 +259,8 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names) precipitation, subsidence, large_scale_advection, + vertical_fluctuation, + nudging, edmf_coriolis, forcing, non_orographic_gravity_wave, diff --git a/src/cache/temporary_quantities.jl b/src/cache/temporary_quantities.jl index 659caea84ed..ed4f1013369 100644 --- a/src/cache/temporary_quantities.jl +++ b/src/cache/temporary_quantities.jl @@ -32,6 +32,7 @@ function temporary_quantities(Y, atmos) temp_data_level_3 = Fields.field_values( Fields.level(Fields.Field(FT, center_space), 1), ), # ρaʲu³ʲ_datah_tot + ᶜtemp_C12 = Fields.Field(C12{FT}, center_space), # ᶜ∇Φ₃ ᶜtemp_C3 = Fields.Field(C3{FT}, center_space), # ᶜ∇Φ₃ ᶜtemp_CT3 = Fields.Field(CT3{FT}, center_space), # ᶜω³, ᶜ∇Φ³ ᶠtemp_CT3 = Fields.Field(CT3{FT}, face_space), # ᶠuₕ³ diff --git a/src/prognostic_equations/forcing/large_scale_advection.jl b/src/prognostic_equations/forcing/large_scale_advection.jl index beebd235a3c..15234f22f3f 100644 --- a/src/prognostic_equations/forcing/large_scale_advection.jl +++ b/src/prognostic_equations/forcing/large_scale_advection.jl @@ -40,7 +40,6 @@ function large_scale_advection_tendency!( cv_v = TD.Parameters.cv_v(thermo_params) R_v = TD.Parameters.R_v(thermo_params) - @. Yₜ.c.ρq_tot[colidx] += Y.c.ρ[colidx] * ᶜdqtdt_hadv # TODO: should `hv` be a thermo function? # (hv = cv_v * (ᶜT - T_0) + Lv_0 - R_v * T_0) @. Yₜ.c.ρe_tot[colidx] += @@ -51,5 +50,6 @@ function large_scale_advection_tendency!( R_v * T_0 ) * ᶜdqtdt_hadv ) + @. Yₜ.c.ρq_tot[colidx] += Y.c.ρ[colidx] * ᶜdqtdt_hadv return nothing end diff --git a/src/prognostic_equations/forcing/nudging.jl b/src/prognostic_equations/forcing/nudging.jl new file mode 100644 index 00000000000..a04a425cf15 --- /dev/null +++ b/src/prognostic_equations/forcing/nudging.jl @@ -0,0 +1,62 @@ +##### +##### Nudging/relaxation tendency for GCM-driven SCM +##### + +import ClimaCore.Fields as Fields + +nudging_cache(Y, atmos::AtmosModel) = nudging_cache(Y, atmos.nudging) + +nudging_cache(Y, ::Nothing) = (;) +nudging_tendency!(Yₜ, Y, p, t, colidx, ::Nothing) = nothing + +function nudging_cache(Y, ::Nudging) + FT = Spaces.undertype(axes(Y.c)) + ᶜT_mean = similar(Y.c, FT) + ᶜq_tot_mean = similar(Y.c, FT) + ᶜu_mean = similar(Y.c, FT) + ᶜv_mean = similar(Y.c, FT) + @. ᶜT_mean = 290 + @. ᶜq_tot_mean = 0.010 + @. ᶜu_mean = 5 + @. ᶜv_mean = 0 + hr = 3600 + ᶜτ_scalar = 24hr + ᶜτ_wind = 6hr + return (; ᶜτ_scalar, ᶜτ_wind, ᶜT_mean, ᶜq_tot_mean, ᶜu_mean, ᶜv_mean) +end + +function nudging_tendency!(Yₜ, Y, p, t, colidx, ::Nudging) + + thermo_params = CAP.thermodynamics_params(p.params) + T_0 = TD.Parameters.T_0(thermo_params) + Lv_0 = TD.Parameters.LH_v0(thermo_params) + cv_v = TD.Parameters.cv_v(thermo_params) + R_v = TD.Parameters.R_v(thermo_params) + + (; ᶜu_mean, ᶜv_mean, ᶜτ_wind, ᶜτ_scalar, ᶜT_mean, ᶜq_tot_mean) = p.nudging + (; ᶜspecific, ᶜts) = p.precomputed + ᶜlg = Fields.local_geometry_field(Y.c) + ᶜuₕ_mean = p.scratch.ᶜtemp_C12 + @. ᶜuₕ_mean[colidx] = + C12(Geometry.UVVector(ᶜu_mean[colidx], ᶜv_mean[colidx]), ᶜlg[colidx]) + @. Yₜ.c.uₕ[colidx] -= (Y.c.uₕ[colidx] - ᶜuₕ_mean[colidx]) / ᶜτ_wind + + ᶜdTdt_nudging = p.scratch.ᶜtemp_scalar + ᶜdqtdt_nudging = p.scratch.ᶜtemp_scalar_2 + @. ᶜdTdt_nudging[colidx] = + -(TD.air_temperature(thermo_params, ᶜts[colidx]) - ᶜT_mean[colidx]) / + ᶜτ_scalar + @. ᶜdqtdt_nudging[colidx] = + -(ᶜspecific.q_tot[colidx] - ᶜq_tot_mean[colidx]) / ᶜτ_scalar + + @. Yₜ.c.ρq_tot[colidx] += Y.c.ρ[colidx] * ᶜdqtdt_nudging[colidx] + + @. Yₜ.c.ρe_tot[colidx] += + Y.c.ρ[colidx] * ( + TD.cp_m(thermo_params, ᶜts[colidx]) * ᶜdTdt_nudging[colidx] + + ( + cv_v * (TD.air_temperature(thermo_params, ᶜts[colidx]) - T_0) + + Lv_0 - R_v * T_0 + ) * ᶜdqtdt_nudging[colidx] + ) +end diff --git a/src/prognostic_equations/forcing/vertical_fluctuation.jl b/src/prognostic_equations/forcing/vertical_fluctuation.jl new file mode 100644 index 00000000000..bef6c8babcc --- /dev/null +++ b/src/prognostic_equations/forcing/vertical_fluctuation.jl @@ -0,0 +1,51 @@ +##### +##### Vertical fluctuation (for single column experiments) +##### + +import Thermodynamics as TD +import ClimaCore.Spaces as Spaces +import ClimaCore.Fields as Fields + +vertical_fluctuation_cache(Y, atmos::AtmosModel) = + vertical_fluctuation_cache(Y, atmos.vert_fluc) + +vertical_fluctuation_cache(Y, ::Nothing) = (;) +function vertical_fluctuation_cache(Y, ::VerticalFluctuation) + FT = Spaces.undertype(axes(Y.c)) + ᶜdTdt_fluc = similar(Y.c, FT) + ᶜdqtdt_fluc = similar(Y.c, FT) + @. ᶜdTdt_fluc = 0 + @. ᶜdqtdt_fluc = 0 + return (; ᶜdTdt_fluc, ᶜdqtdt_fluc) +end + +vertical_fluctuation_tendency!(Yₜ, Y, p, t, colidx, ::Nothing) = nothing +function vertical_fluctuation_tendency!( + Yₜ, + Y, + p, + t, + colidx, + ::VerticalFluctuation, +) + (; params) = p + (; ᶜts) = p.precomputed + (; ᶜdTdt_fluc, ᶜdqtdt_fluc) = p.vertical_fluctuation + thermo_params = CAP.thermodynamics_params(params) + T_0 = TD.Parameters.T_0(thermo_params) + Lv_0 = TD.Parameters.LH_v0(thermo_params) + cv_v = TD.Parameters.cv_v(thermo_params) + R_v = TD.Parameters.R_v(thermo_params) + + @. Yₜ.c.ρe_tot[colidx] += + Y.c.ρ[colidx] * ( + TD.cp_m(thermo_params, ᶜts[colidx]) * ᶜdTdt_fluc[colidx] + + ( + cv_v * (TD.air_temperature(thermo_params, ᶜts[colidx]) - T_0) + Lv_0 - + R_v * T_0 + ) * ᶜdqtdt_fluc[colidx] + ) + @. Yₜ.c.ρq_tot[colidx] += Y.c.ρ[colidx] * ᶜdqtdt_fluc[colidx] + + return nothing +end diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl index c5f6c49e88a..ea7d958d4b4 100644 --- a/src/prognostic_equations/remaining_tendency.jl +++ b/src/prognostic_equations/remaining_tendency.jl @@ -33,6 +33,8 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t) subsidence_tendency!(Yₜ, Y, p, t, colidx, p.atmos.subsidence) edmf_coriolis_tendency!(Yₜ, Y, p, t, colidx, p.atmos.edmf_coriolis) large_scale_advection_tendency!(Yₜ, Y, p, t, colidx, p.atmos.ls_adv) + vertical_fluctuation_tendency!(Yₜ, Y, p, t, colidx, p.atmos.vert_fluc) + nudging_tendency!(Yₜ, Y, p, t, colidx, p.atmos.nudging) if p.atmos.sgs_adv_mode == Explicit() edmfx_sgs_vertical_advection_tendency!( diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index f0da55e4096..5f93dffbe65 100644 --- a/src/solver/model_getters.jl +++ b/src/solver/model_getters.jl @@ -333,18 +333,22 @@ function get_large_scale_advection_model(parsed_args, ::Type{FT}) where {FT} (APL.ARM_SGP_dTdt(FT), APL.ARM_SGP_dqtdt(FT)) elseif ls_adv == "GATE_III" (APL.GATE_III_dTdt(FT), APL.GATE_III_dqtdt(FT)) + elseif ls_adv == "GCMDriven" + (APL.Bomex_dTdt(FT), APL.Bomex_dqtdt(FT)) else error("Uncaught case") end # See https://clima.github.io/AtmosphericProfilesLibrary.jl/dev/ # for which functions accept which arguments. - prof_dqtdt = if ls_adv in ("Bomex", "LifeCycleTan2018", "Rico", "GATE_III") - (thermo_params, ᶜts, t, z) -> prof_dqtdt₀(z) - elseif ls_adv == "ARM_SGP" - (thermo_params, ᶜts, t, z) -> - prof_dqtdt₀(TD.exner(thermo_params, ᶜts), t, z) - end - prof_dTdt = if ls_adv in ("Bomex", "LifeCycleTan2018", "Rico") + prof_dqtdt = + if ls_adv in + ("Bomex", "LifeCycleTan2018", "Rico", "GATE_III", "GCMDriven") + (thermo_params, ᶜts, t, z) -> prof_dqtdt₀(z) + elseif ls_adv == "ARM_SGP" + (thermo_params, ᶜts, t, z) -> + prof_dqtdt₀(TD.exner(thermo_params, ᶜts), t, z) + end + prof_dTdt = if ls_adv in ("Bomex", "LifeCycleTan2018", "Rico", "GCMDriven") (thermo_params, ᶜts, t, z) -> prof_dTdt₀(TD.exner(thermo_params, ᶜts), z) elseif ls_adv == "ARM_SGP" @@ -356,6 +360,26 @@ function get_large_scale_advection_model(parsed_args, ::Type{FT}) where {FT} return LargeScaleAdvection(prof_dTdt, prof_dqtdt) end +function get_vertical_fluctuation_model(parsed_args, ::Type{FT}) where {FT} + vert_fluc_name = parsed_args["vert_fluc"] + @assert vert_fluc_name in (nothing, "GCMDriven") + return if isnothing(vert_fluc_name) + nothing + elseif vert_fluc_name == "GCMDriven" + VerticalFluctuation() + end +end + +function get_nudging_model(parsed_args) + nudging_name = parsed_args["nudging"] + @assert nudging_name in (nothing, "GCMDriven") + return if isnothing(nudging_name) + nothing + elseif nudging_name == "GCMDriven" + Nudging() + end +end + function get_edmf_coriolis(parsed_args, ::Type{FT}) where {FT} edmf_coriolis = parsed_args["edmf_coriolis"] edmf_coriolis == nothing && return nothing diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index 98e4d0de71a..238744fb1cb 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -60,6 +60,8 @@ function get_atmos(config::AtmosConfig, params) radiation_mode, subsidence = get_subsidence_model(parsed_args, radiation_mode, FT), ls_adv = get_large_scale_advection_model(parsed_args, FT), + vert_fluc = get_vertical_fluctuation_model(parsed_args, FT), + nudging = get_nudging_model(parsed_args), edmf_coriolis = get_edmf_coriolis(parsed_args, FT), advection_test, tendency_model = get_tendency_model(parsed_args), diff --git a/src/solver/types.jl b/src/solver/types.jl index 3ba4df6d10f..0dd44ed18e0 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -114,6 +114,9 @@ struct LargeScaleAdvection{PT, PQ} prof_dTdt::PT # Set large-scale cooling prof_dqtdt::PQ # Set large-scale drying end +# maybe need to <: AbstractForcing +struct VerticalFluctuation end +struct Nudging end struct EDMFCoriolis{U, V, FT} prof_ug::U @@ -332,6 +335,8 @@ Base.@kwdef struct AtmosModel{ S, RM, LA, + VF, + NUDGING, EC, AT, TM, @@ -365,6 +370,8 @@ Base.@kwdef struct AtmosModel{ subsidence::S = nothing radiation_mode::RM = nothing ls_adv::LA = nothing + vert_fluc::VF = nothing + nudging::NUDGING = nothing edmf_coriolis::EC = nothing advection_test::AT = nothing tendency_model::TM = nothing diff --git a/test/coupler_compatibility.jl b/test/coupler_compatibility.jl index 541f483f4b0..b8394d09501 100644 --- a/test/coupler_compatibility.jl +++ b/test/coupler_compatibility.jl @@ -77,6 +77,8 @@ const T2 = 290 p.precipitation, p.subsidence, p.large_scale_advection, + p.vertical_fluctuation, + p.nudging, p.edmf_coriolis, p.forcing, p.non_orographic_gravity_wave,