From d1cea8dd27af5d380e2376c7ba2f64c3257da60e 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 --- 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/prognostic_equations/forcing/nudging.jl | 64 +++++++++++++++++++ .../forcing/vertical_fluctuation.jl | 55 ++++++++++++++++ .../remaining_tendency.jl | 2 + src/solver/model_getters.jl | 51 +++++++++++++-- src/solver/type_getters.jl | 2 + src/solver/types.jl | 9 +++ 11 files changed, 233 insertions(+), 7 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/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/prognostic_equations/forcing/nudging.jl b/src/prognostic_equations/forcing/nudging.jl new file mode 100644 index 00000000000..f087a3d94ae --- /dev/null +++ b/src/prognostic_equations/forcing/nudging.jl @@ -0,0 +1,64 @@ +##### +##### 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 = 6hr + ᶜτ_wind = 1hr + 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_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) * ᶜdTdt_nudging[colidx] + ( + cv_v * (TD.air_temperature(thermo_params, ᶜts[colidx]) - T_0) + Lv_0 - + R_v * T_0 + ) * ᶜdqtdt_nudging[colidx] + ) +end \ No newline at end of file diff --git a/src/prognostic_equations/forcing/vertical_fluctuation.jl b/src/prognostic_equations/forcing/vertical_fluctuation.jl new file mode 100644 index 00000000000..f253ed271b2 --- /dev/null +++ b/src/prognostic_equations/forcing/vertical_fluctuation.jl @@ -0,0 +1,55 @@ +##### +##### 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, vert_fluc::Nothing) = (;) +function vertical_fluctuation_cache(Y, vert_fluc::VerticalFluctuation) + FT = Spaces.undertype(axes(Y.c)) + return (; ᶜdqtdt_fluc = similar(Y.c, FT), ᶜdTdt_fluc = similar(Y.c, FT)) +end + +vertical_fluctuation_tendency!(Yₜ, Y, p, t, colidx, ::Nothing) = nothing +function vertical_fluctuation_tendency!( + Yₜ, + Y, + p, + t, + colidx, + vert_fluc::VerticalFluctuation, +) + (; prof_dTdt, prof_dqtdt) = vert_fluc + + thermo_params = CAP.thermodynamics_params(p.params) + ᶜts = p.precomputed.ᶜts[colidx] + ᶜdqtdt_fluc = p.vertical_fluctuation.ᶜdqtdt_fluc[colidx] + ᶜdTdt_fluc = p.vertical_fluctuation.ᶜdTdt_fluc[colidx] + z = Fields.coordinate_field(axes(ᶜdqtdt_fluc)).z + + @. ᶜdTdt_fluc = prof_dTdt(thermo_params, ᶜts, t, z) + @. ᶜdqtdt_fluc = prof_dqtdt(thermo_params, ᶜts, t, z) + + 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.ρq_tot[colidx] += Y.c.ρ[colidx] * ᶜdqtdt_fluc + # TODO: should `hv` be a thermo function? + # (hv = cv_v * (ᶜT - T_0) + Lv_0 - R_v * T_0) + @. Yₜ.c.ρe_tot[colidx] += + Y.c.ρ[colidx] * ( + TD.cp_m(thermo_params, ᶜts) * ᶜdTdt_fluc + + ( + cv_v * (TD.air_temperature(thermo_params, ᶜts) - T_0) + Lv_0 - + R_v * T_0 + ) * ᶜdqtdt_fluc + ) + 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..67d234e2efd 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,39 @@ 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 = parsed_args["vert_fluc"] + vert_fluc == nothing && return nothing + + (prof_dTdt₀, prof_dqtdt₀) = if vert_fluc == "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 vert_fluc == "GCMDriven" + (thermo_params, ᶜts, t, z) -> prof_dqtdt₀(z) + end + prof_dTdt = if vert_fluc == "GCMDriven" + (thermo_params, ᶜts, t, z) -> + prof_dTdt₀(TD.exner(thermo_params, ᶜts), z) + end + + return VerticalFluctuation(prof_dTdt, prof_dqtdt) +end + +function get_nudging_model(parsed_args) + nudging_name = parsed_args["nudging"] + return if isnothing(nudging_name) + nothing + elseif nudging_name == "GCMDriven" + Nudging() + else + error("Uncaught nudging model `$nudging_name`.") + 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..73404ad7765 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -114,6 +114,11 @@ struct LargeScaleAdvection{PT, PQ} prof_dTdt::PT # Set large-scale cooling prof_dqtdt::PQ # Set large-scale drying end +struct VerticalFluctuation{PT, PQ} + prof_dTdt::PT # Set large-scale cooling + prof_dqtdt::PQ # Set large-scale drying +end +struct Nudging end # maybe need to <: AbstractForcing struct EDMFCoriolis{U, V, FT} prof_ug::U @@ -332,6 +337,8 @@ Base.@kwdef struct AtmosModel{ S, RM, LA, + VF, + NUDGING, EC, AT, TM, @@ -365,6 +372,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