Move hyperdif into numerics, clean up

Author: nefrathenrici

src/diagnostics/default_diagnostics.jl

@@ -215,15 +215,15 @@ function default_diagnostics(
 end
 
 function default_diagnostics(
-    atmos_hydrology::AtmosWater,
+    atmos_water::AtmosWater,
     duration,
     start_date;
     output_writer,
 )
     diagnostics = []
 
     # Add moisture and precipitation model diagnostics
-    for model in (atmos_hydrology.moisture_model, atmos_hydrology.precip_model)
+    for model in (atmos_water.moisture_model, atmos_water.precip_model)
         !isnothing(model) && append!(
             diagnostics,
             default_diagnostics(model, duration, start_date; output_writer),

src/parameterized_tendencies/radiation/radiation.jl

@@ -28,8 +28,12 @@ radiation_model_cache(Y, atmos::AtmosModel, args...) =
 ##### No Radiation
 #####
 
-radiation_model_cache(Y, radiation_mode::Nothing; args...) = (;)
-radiation_tendency!(Yₜ, Y, p, t, ::Nothing) = nothing
+radiation_model_cache(
+    Y,
+    radiation_mode::Union{Nothing, HeldSuarezForcing};
+    args...,
+) = (;)
+radiation_tendency!(Yₜ, Y, p, t, ::Union{Nothing, HeldSuarezForcing}) = nothing
 
 #####
 ##### RRTMGP Radiation

src/prognostic_equations/remaining_tendency.jl

@@ -141,7 +141,8 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
     ᶜuₕ = Y.c.uₕ
     ᶠu₃ = Y.f.u₃
     ᶜρ = Y.c.ρ
-    (; radiation_mode, moisture_model, rayleigh_sponge, viscous_sponge) = p.atmos
+    (; radiation_mode, moisture_model, rayleigh_sponge, viscous_sponge) =
+        p.atmos
     (; ls_adv, scm_coriolis) = p.atmos
     (; params) = p
     thermo_params = CAP.thermodynamics_params(params)

src/solver/model_getters.jl

@@ -246,6 +246,7 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
         "gray",
         "allsky",
         "allskywithclear",
+        "held_suarez",
         "DYCOMS",
         "TRMM_LBA",
         "ISDAC",
@@ -286,6 +287,8 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
             reset_rng_seed,
             deep_atmosphere,
         )
+    elseif radiation_name == "held_suarez"
+        HeldSuarezForcing()
     elseif radiation_name == "DYCOMS"
         RadiationDYCOMS{FT}()
     elseif radiation_name == "TRMM_LBA"
@@ -352,11 +355,11 @@ end
 function get_forcing_type(parsed_args)
     forcing = parsed_args["forcing"]
     @assert forcing in (nothing, "held_suarez")
-    return if forcing == nothing
-        nothing
-    elseif forcing == "held_suarez"
-        HeldSuarezForcing()
+    if forcing == "held_suarez"
+        @warn "The 'held_suarez' forcing option is deprecated. Use rad='held_suarez' instead to set HeldSuarezForcing as a radiation mode."
+        return HeldSuarezForcing()  # Still return the object for backward compatibility
     end
+    return nothing
 end
 
 struct CallCloudDiagnosticsPerStage end
@@ -372,7 +375,7 @@ end
 
 function get_subsidence_model(parsed_args, radiation_mode, FT)
     subsidence = parsed_args["subsidence"]
-    subsidence == nothing && return nothing
+    isnothing(subsidence) && return nothing
 
     prof = if subsidence == "Bomex"
         APL.Bomex_subsidence(FT)
@@ -382,6 +385,8 @@ function get_subsidence_model(parsed_args, radiation_mode, FT)
         APL.Rico_subsidence(FT)
     elseif subsidence == "DYCOMS"
         @assert radiation_mode isa RadiationDYCOMS
+        # For DYCOMS case, subsidence is linearly proportional to height
+        # with slope equal to the divergence rate specified in radiation mode
         z -> -z * radiation_mode.divergence
     elseif subsidence == "ISDAC"
         APL.ISDAC_subsidence(FT)

src/solver/type_getters.jl

@@ -48,7 +48,11 @@ function get_atmos(config::AtmosConfig, params)
     (!isnothing(co2) && !with_rrtmgp) &&
         @warn ("$(co2) does nothing if RRTMGP is not used")
 
-    disable_momentum_vertical_diffusion = forcing_type isa HeldSuarezForcing
+    # HeldSuarezForcing can be set via radiation_mode or legacy forcing option for now
+    final_radiation_mode =
+        forcing_type isa HeldSuarezForcing ? forcing_type : radiation_mode
+    disable_momentum_vertical_diffusion =
+        final_radiation_mode isa HeldSuarezForcing
 
     advection_test = parsed_args["advection_test"]
     @assert advection_test in (false, true)
@@ -102,11 +106,10 @@ function get_atmos(config::AtmosConfig, params)
         scm_coriolis = get_scm_coriolis(parsed_args, FT),
 
         # AtmosRadiation
-        radiation_mode,
+        radiation_mode = final_radiation_mode,
         ozone,
         co2,
         insolation = get_insolation_form(parsed_args),
-        held_suarez_forcing = forcing_type,
 
         # AtmosTurbconv - Turbulence & Convection
         edmfx_model,
@@ -139,8 +142,7 @@ function get_atmos(config::AtmosConfig, params)
 
         # Top-level options (not grouped)
         vert_diff,
-        hyperdiff = get_hyperdiffusion_model(parsed_args, FT),
-        numerics = get_numerics(parsed_args),
+        numerics = get_numerics(parsed_args, FT),
         disable_surface_flux_tendency = parsed_args["disable_surface_flux_tendency"],
     )
     # TODO: Should this go in the AtmosModel constructor?
@@ -161,7 +163,7 @@ function get_scale_blending_method(parsed_args)
     end
 end
 
-function get_numerics(parsed_args)
+function get_numerics(parsed_args, FT)
     test_dycore =
         parsed_args["test_dycore_consistency"] ? TestDycoreConsistency() :
         nothing
@@ -189,6 +191,9 @@ function get_numerics(parsed_args)
 
     # wrap each upwinding mode in a Val for dispatch
     diff_mode = parsed_args["implicit_diffusion"] ? Implicit() : Explicit()
+
+    hyperdiff = get_hyperdiffusion_model(parsed_args, FT)
+
     numerics = AtmosNumerics(;
         energy_upwinding,
         tracer_upwinding,
@@ -197,6 +202,7 @@ function get_numerics(parsed_args)
         limiter,
         test_dycore_consistency = test_dycore,
         diff_mode,
+        hyperdiff,
     )
     @info "numerics $(summary(numerics))"