diff --git a/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl b/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl index b8e59278c52..548fa8aefae 100644 --- a/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl +++ b/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl @@ -2,144 +2,152 @@ ##### Non-orographic gravity wave parameterization ##### using UnrolledUtilities +import ClimaCore.Domains +import ClimaCore.Meshes import ClimaCore.Spaces as Spaces import ClimaCore.Fields as Fields import ClimaCore.Geometry as Geometry import ClimaCore.Operators as Operator non_orographic_gravity_wave_cache(Y, atmos::AtmosModel) = - non_orographic_gravity_wave_cache( - Y, - atmos.non_orographic_gravity_wave, - atmos.model_config, - ) + non_orographic_gravity_wave_cache(Y, atmos.non_orographic_gravity_wave) -non_orographic_gravity_wave_cache(Y, ::Nothing, ::AbstractModelConfig) = (;) +non_orographic_gravity_wave_cache(Y, ::Nothing) = (;) non_orographic_gravity_wave_tendency!(Yₜ, Y, p, t, ::Nothing) = nothing -function non_orographic_gravity_wave_cache( - Y, - gw::NonOrographyGravityWave, - ::SingleColumnModel, -) - FT = Spaces.undertype(axes(Y.c)) - (; source_height, Bw, Bn, Bt_0, dc, cmax, c0, nk, cw, cn) = gw - - nc = Int(floor(FT(2 * cmax / dc + 1))) - c = ntuple(n -> FT((n - 1) * dc - cmax), Val(nc)) - source_ρ_z_u_v_level = - similar(Fields.level(Y.c.ρ, 1), Tuple{FT, FT, FT, FT, FT}) - ᶜlevel = similar(Y.c.ρ, FT) - for i in 1:Spaces.nlevels(axes(Y.c.ρ)) - fill!(Fields.level(ᶜlevel, i), i) - end +function iscolumn(space) + # TODO: Our columns are 2+1D boxes with one element at the base. Fix this + isbox = + Meshes.domain(Spaces.topology(Spaces.horizontal_space(space))) isa + Domains.RectangleDomain + isbox || return false + has_one_element = + Meshes.nelements(Spaces.topology(Spaces.horizontal_space(Y)).mesh) == 1 + has_one_element && return true +end - return (; - gw_source_height = source_height, - gw_source_ampl = Bt_0 .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), - gw_Bw = Bw .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), - gw_Bn = Bn .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), - gw_c = c, - gw_dc = dc, - gw_cmax = cmax, - gw_cw = cw .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), - gw_cn = cn .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), - gw_c0 = c0, - gw_flag = ones(FT, axes(Fields.level(Y.c.ρ, 1))), - gw_nk = Int(nk), - ᶜbuoyancy_frequency = similar(Y.c.ρ), - ᶜdTdz = similar(Y.c.ρ), - source_ρ_z_u_v_level, - source_level = similar(Fields.level(Y.c.ρ, 1)), - damp_level = similar(Fields.level(Y.c.ρ, 1)), - ᶜlevel, - u_waveforcing = similar(Y.c.ρ), - v_waveforcing = similar(Y.c.ρ), - uforcing = similar(Y.c.ρ), - vforcing = similar(Y.c.ρ), - gw_ncval = Val(nc), - ) +function issphere(space) + return Meshes.domain(Spaces.topology(Spaces.horizontal_space(space))) isa + Domains.SphereDomain end -function non_orographic_gravity_wave_cache( - Y, - gw::NonOrographyGravityWave, - ::SphericalModel, -) +function non_orographic_gravity_wave_cache(Y, gw::NonOrographyGravityWave) + if iscolumn(axes(Y.c)) + FT = Spaces.undertype(axes(Y.c)) + (; source_height, Bw, Bn, Bt_0, dc, cmax, c0, nk, cw, cn) = gw + + nc = Int(floor(FT(2 * cmax / dc + 1))) + c = ntuple(n -> FT((n - 1) * dc - cmax), Val(nc)) + source_ρ_z_u_v_level = + similar(Fields.level(Y.c.ρ, 1), Tuple{FT, FT, FT, FT, FT}) + ᶜlevel = similar(Y.c.ρ, FT) + for i in 1:Spaces.nlevels(axes(Y.c.ρ)) + fill!(Fields.level(ᶜlevel, i), i) + end - FT = Spaces.undertype(axes(Y.c)) - (; source_pressure, damp_pressure, Bw, Bn, Bt_0, Bt_n, Bt_s, Bt_eq) = gw - (; ϕ0_s, ϕ0_n, dϕ_n, dϕ_s, dc, cmax, c0, nk, cw, cw_tropics, cn) = gw - - nc = Int(floor(FT(2 * cmax / dc + 1))) - c = ntuple(n -> FT((n - 1) * dc - cmax), Val(nc)) - - ᶜlocal_geometry = Fields.local_geometry_field(Fields.level(Y.c, 1)) - lat = ᶜlocal_geometry.coordinates.lat - - gw_Bn = @. ifelse(dϕ_s <= lat <= dϕ_n, FT(0), Bn) - gw_cw = @. ifelse(dϕ_s <= lat <= dϕ_n, cw_tropics, cw) - gw_flag = @. ifelse(dϕ_s <= lat <= dϕ_n, FT(0), FT(1)) - gw_Bw = ones(FT, axes(lat)) .* Bw - gw_cn = ones(FT, axes(lat)) .* cn - - source_p_ρ_z_u_v_level = - similar(Fields.level(Y.c.ρ, 1), Tuple{FT, FT, FT, FT, FT, FT}) - ᶜlevel = similar(Y.c.ρ, FT) - for i in 1:Spaces.nlevels(axes(Y.c.ρ)) - fill!(Fields.level(ᶜlevel, i), i) - end + return (; + gw_source_height = source_height, + gw_source_ampl = Bt_0 .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), + gw_Bw = Bw .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), + gw_Bn = Bn .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), + gw_c = c, + gw_dc = dc, + gw_cmax = cmax, + gw_cw = cw .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), + gw_cn = cn .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), + gw_c0 = c0, + gw_flag = ones(FT, axes(Fields.level(Y.c.ρ, 1))), + gw_nk = Int(nk), + ᶜbuoyancy_frequency = similar(Y.c.ρ), + ᶜdTdz = similar(Y.c.ρ), + source_ρ_z_u_v_level, + source_level = similar(Fields.level(Y.c.ρ, 1)), + damp_level = similar(Fields.level(Y.c.ρ, 1)), + ᶜlevel, + u_waveforcing = similar(Y.c.ρ), + v_waveforcing = similar(Y.c.ρ), + uforcing = similar(Y.c.ρ), + vforcing = similar(Y.c.ρ), + gw_ncval = Val(nc), + ) + elseif issphere(axes(Y.c)) + + FT = Spaces.undertype(axes(Y.c)) + (; source_pressure, damp_pressure, Bw, Bn, Bt_0, Bt_n, Bt_s, Bt_eq) = gw + (; ϕ0_s, ϕ0_n, dϕ_n, dϕ_s, dc, cmax, c0, nk, cw, cw_tropics, cn) = gw + nc = Int(floor(FT(2 * cmax / dc + 1))) + c = ntuple(n -> FT((n - 1) * dc - cmax), Val(nc)) - # This is GFDL source specs -> a smooth function - # source_ampl = @. Bt_0 + - # Bt_n * FT(0.5) * (FT(1) + tanh((lat - ϕ0_n) / dϕ_n)) + - # Bt_s * FT(0.5) * (FT(1) + tanh((lat - ϕ0_s) / dϕ_s)) - - # This latitude depend source follows MiMA specs - source_ampl = @. ifelse( - (lat > ϕ0_n) | (lat < ϕ0_s), - Bt_0 + - Bt_n * FT(0.5) * (FT(1) + tanh((lat - ϕ0_n) / dϕ_n)) + - Bt_s * FT(0.5) * (FT(1) + tanh((lat - ϕ0_s) / dϕ_s)), - ifelse( - dϕ_s <= lat <= dϕ_n, - Bt_eq, + ᶜlocal_geometry = Fields.local_geometry_field(Fields.level(Y.c, 1)) + lat = ᶜlocal_geometry.coordinates.lat + + gw_Bn = @. ifelse(dϕ_s <= lat <= dϕ_n, FT(0), Bn) + gw_cw = @. ifelse(dϕ_s <= lat <= dϕ_n, cw_tropics, cw) + gw_flag = @. ifelse(dϕ_s <= lat <= dϕ_n, FT(0), FT(1)) + gw_Bw = ones(FT, axes(lat)) .* Bw + gw_cn = ones(FT, axes(lat)) .* cn + + source_p_ρ_z_u_v_level = + similar(Fields.level(Y.c.ρ, 1), Tuple{FT, FT, FT, FT, FT, FT}) + ᶜlevel = similar(Y.c.ρ, FT) + for i in 1:Spaces.nlevels(axes(Y.c.ρ)) + fill!(Fields.level(ᶜlevel, i), i) + end + + + # This is GFDL source specs -> a smooth function + # source_ampl = @. Bt_0 + + # Bt_n * FT(0.5) * (FT(1) + tanh((lat - ϕ0_n) / dϕ_n)) + + # Bt_s * FT(0.5) * (FT(1) + tanh((lat - ϕ0_s) / dϕ_s)) + + # This latitude depend source follows MiMA specs + source_ampl = @. ifelse( + (lat > ϕ0_n) | (lat < ϕ0_s), + Bt_0 + + Bt_n * FT(0.5) * (FT(1) + tanh((lat - ϕ0_n) / dϕ_n)) + + Bt_s * FT(0.5) * (FT(1) + tanh((lat - ϕ0_s) / dϕ_s)), ifelse( - dϕ_n <= lat <= ϕ0_n, - Bt_0 + (Bt_eq - Bt_0) / (ϕ0_n - dϕ_n) * (ϕ0_n - lat), - Bt_0 + (Bt_eq - Bt_0) / (ϕ0_s - dϕ_s) * (ϕ0_s - lat), + dϕ_s <= lat <= dϕ_n, + Bt_eq, + ifelse( + dϕ_n <= lat <= ϕ0_n, + Bt_0 + (Bt_eq - Bt_0) / (ϕ0_n - dϕ_n) * (ϕ0_n - lat), + Bt_0 + (Bt_eq - Bt_0) / (ϕ0_s - dϕ_s) * (ϕ0_s - lat), + ), ), - ), - ) + ) - return (; - gw_source_pressure = source_pressure, - gw_damp_pressure = damp_pressure, - gw_source_ampl = source_ampl, - gw_Bw = gw_Bw, - gw_Bn = gw_Bn, - gw_c = c, - gw_cw = gw_cw, - gw_cn = gw_cn, - gw_dc = dc, - gw_cmax = cmax, - gw_c0 = c0, - gw_flag = gw_flag, - gw_nk = Int(nk), - ᶜbuoyancy_frequency = similar(Y.c.ρ), - ᶜdTdz = similar(Y.c.ρ), - source_p_ρ_z_u_v_level, - source_level = similar(Fields.level(Y.c.ρ, 1)), - damp_level = similar(Fields.level(Y.c.ρ, 1)), - ᶜlevel, - u_waveforcing = similar(Y.c.ρ), - v_waveforcing = similar(Y.c.ρ), - uforcing = similar(Y.c.ρ), - vforcing = similar(Y.c.ρ), - gw_ncval = Val(nc), - ) + return (; + gw_source_pressure = source_pressure, + gw_damp_pressure = damp_pressure, + gw_source_ampl = source_ampl, + gw_Bw = gw_Bw, + gw_Bn = gw_Bn, + gw_c = c, + gw_cw = gw_cw, + gw_cn = gw_cn, + gw_dc = dc, + gw_cmax = cmax, + gw_c0 = c0, + gw_flag = gw_flag, + gw_nk = Int(nk), + ᶜbuoyancy_frequency = similar(Y.c.ρ), + ᶜdTdz = similar(Y.c.ρ), + source_p_ρ_z_u_v_level, + source_level = similar(Fields.level(Y.c.ρ, 1)), + damp_level = similar(Fields.level(Y.c.ρ, 1)), + ᶜlevel, + u_waveforcing = similar(Y.c.ρ), + v_waveforcing = similar(Y.c.ρ), + uforcing = similar(Y.c.ρ), + vforcing = similar(Y.c.ρ), + gw_ncval = Val(nc), + ) + else + error("Only sphere and columns are supported") + end end function non_orographic_gravity_wave_tendency!( @@ -165,7 +173,6 @@ function non_orographic_gravity_wave_tendency!( ᶜlevel, gw_ncval, ) = p.non_orographic_gravity_wave - (; model_config) = p.atmos ᶜρ = Y.c.ρ ᶜz = Fields.coordinate_field(Y.c).z @@ -191,7 +198,7 @@ function non_orographic_gravity_wave_tendency!( ᶜu = Geometry.UVVector.(Y.c.uₕ).components.data.:1 ᶜv = Geometry.UVVector.(Y.c.uₕ).components.data.:2 - if model_config isa SingleColumnModel + if iscolumn(axes(Y.c)) # source level: the index of the level that is closest to the source height (; gw_source_height, source_ρ_z_u_v_level) = p.non_orographic_gravity_wave @@ -216,7 +223,7 @@ function non_orographic_gravity_wave_tendency!( fill!(damp_level, Spaces.nlevels(axes(ᶜz))) - elseif model_config isa SphericalModel + elseif issphere(axes(Y.c)) (; ᶜp) = p.precomputed (; gw_source_pressure, gw_damp_pressure, source_p_ρ_z_u_v_level) = p.non_orographic_gravity_wave @@ -255,7 +262,8 @@ function non_orographic_gravity_wave_tendency!( return (level_prev, p_prev) end end - + else + error("Only sphere and columns are supported") end ᶜu = Geometry.UVVector.(Y.c.uₕ).components.data.:1 diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index 29ef6a11603..0bf6a2d9a56 100644 --- a/src/solver/model_getters.jl +++ b/src/solver/model_getters.jl @@ -10,30 +10,6 @@ function get_moisture_model(parsed_args) end end -function get_model_config(parsed_args) - config = parsed_args["config"] - - valid_configurations = ("sphere", "column", "box", "plane") - - if !(config ∈ valid_configurations) - error_message = string( - "config = $config is not one of the ", - "valid configurations $valid_configurations", - ) - throw(ArgumentError(error_message)) - end - - return if config == "sphere" - SphericalModel() - elseif config == "column" - SingleColumnModel() - elseif config == "box" - BoxModel() - elseif config == "plane" - PlaneModel() - end -end - function get_sfc_temperature_form(parsed_args) surface_temperature = parsed_args["surface_temperature"] @assert surface_temperature in @@ -187,9 +163,9 @@ function get_non_orographic_gravity_wave_model( nogw_name = parsed_args["non_orographic_gravity_wave"] @assert nogw_name in (true, false) return if nogw_name == true - if model_config isa SingleColumnModel + if parsed_args["model_config"] == "column" NonOrographyGravityWave{FT}(; Bw = 1.2, Bn = 0.0, Bt_0 = 4e-3) - elseif model_config isa SphericalModel + elseif parsed_args["model_config"] == "sphere" NonOrographyGravityWave{FT}(; Bw = 0.4, Bn = 0.0, diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index 952f51a03e9..51cde77ff1d 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -50,13 +50,11 @@ function get_atmos(config::AtmosConfig, params) filter = Val(parsed_args["edmfx_filter"]), ) - model_config = get_model_config(parsed_args) vert_diff = get_vertical_diffusion_model(diffuse_momentum, parsed_args, params, FT) atmos = AtmosModel(; moisture_model, - model_config, ozone, radiation_mode, subsidence = get_subsidence_model(parsed_args, radiation_mode, FT), @@ -73,7 +71,6 @@ function get_atmos(config::AtmosConfig, params) turbconv_model = get_turbconv_model(FT, parsed_args, turbconv_params), non_orographic_gravity_wave = get_non_orographic_gravity_wave_model( parsed_args, - model_config, FT, ), orographic_gravity_wave = get_orographic_gravity_wave_model( diff --git a/src/solver/types.jl b/src/solver/types.jl index 63943892233..7bcf0af04b6 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -21,12 +21,6 @@ struct GridScaleCloud <: AbstractCloudModel end struct QuadratureCloud <: AbstractCloudModel end struct SGSQuadratureCloud <: AbstractCloudModel end -abstract type AbstractModelConfig end -struct SingleColumnModel <: AbstractModelConfig end -struct SphericalModel <: AbstractModelConfig end -struct BoxModel <: AbstractModelConfig end -struct PlaneModel <: AbstractModelConfig end - abstract type AbstractSST end struct ZonallySymmetricSST <: AbstractSST end struct ZonallyAsymmetricSST <: AbstractSST end @@ -420,7 +414,6 @@ Base.@kwdef struct EDMFXModel{ end Base.@kwdef struct AtmosModel{ - MC, MM, PM, CM, @@ -451,7 +444,6 @@ Base.@kwdef struct AtmosModel{ SA, NUM, } - model_config::MC = nothing moisture_model::MM = nothing precip_model::PM = nothing cloud_model::CM = nothing diff --git a/test/solver/model_getters.jl b/test/solver/model_getters.jl index fb9b0e37da7..66ef4db1b77 100644 --- a/test/solver/model_getters.jl +++ b/test/solver/model_getters.jl @@ -2,16 +2,6 @@ using ClimaComms ClimaComms.@import_required_backends import ClimaAtmos as CA -@testset "Model config" begin - config = CA.AtmosConfig(Dict("config" => "box"), job_id = "model1") - parsed_args = config.parsed_args - @test CA.get_model_config(parsed_args) isa CA.BoxModel - - config = CA.AtmosConfig(Dict("config" => "plane"), job_id = "model2") - parsed_args = config.parsed_args - @test CA.get_model_config(parsed_args) isa CA.PlaneModel -end - @testset "Hyperdiffusion config" begin @info "CAM_SE (Special case of ClimaHyperdiffusion)" config = CA.AtmosConfig(