Merge pull request #2361 from CliMA/zs/remove_theta

remove theta as an energy variable

.buildkite/pipeline.yml

@@ -279,19 +279,6 @@ steps:
         command: "julia --color=yes --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/sphere_baroclinic_wave_rhoe_equilmoist_conservation_ft64.yml"
         artifact_paths: "sphere_baroclinic_wave_rhoe_equilmoist_conservation_ft64/*"
 
-      - label: ":computer: held suarez (ρθ)"
-        command: >
-          julia --color=yes --project=examples examples/hybrid/driver.jl
-          --config_file $CONFIG_PATH/sphere_held_suarez_rhotheta.yml
-
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_held_suarez_rhotheta --out_dir sphere_held_suarez_rhotheta
-
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
-          --nc_dir sphere_held_suarez_rhotheta --fig_dir sphere_held_suarez_rhotheta
-          --case_name dry_held_suarez
-        artifact_paths: "sphere_held_suarez_rhotheta/*"
-
       - label: ":computer: held suarez (ρe) hightop"
         key: sphere_held_suarez_rhoe_hightop
         command: >

config/default_configs/default_config.yml

@@ -138,9 +138,6 @@ forcing:
 test_dycore_consistency:
   help: "Test dycore consistency [`false` (default), `true`]"
   value: false
-energy_name:
-  help: "Energy variable name [`rhoe` (default), `rhotheta`]"
-  value: "rhoe"
 dt_save_to_disk:
   help: "Time between saving to disk. Examples: [`10secs`, `1hours`, `Inf` (do not save)]"
   value: "Inf"

post_processing/post_processing_funcs.jl

@@ -7,12 +7,7 @@ import ClimaCoreTempestRemap: def_space_coord
 import ClimaCoreSpectra: power_spectrum_1d, power_spectrum_2d
 import ClimaCore: Geometry, Fields, Spaces
 using LinearAlgebra: norm
-import ClimaAtmos:
-    DryModel,
-    EquilMoistModel,
-    NonEquilMoistModel,
-    PotentialTemperature,
-    TotalEnergy
+import ClimaAtmos: DryModel, EquilMoistModel, NonEquilMoistModel
 
 function process_name(s::AbstractString)
     # "c_ρ", "c_ρe", "c_uₕ_1", "c_uₕ_2", "f_w_1"

post_processing/remap/remap_helpers.jl

@@ -110,14 +110,7 @@ function remap2latlon(filein, data_dir, remap_tmpdir, weightfile, nlat, nlon)
     nc_time = def_time_coord(nc)
     # define variables for the prognostic states
     nc_rho = defVar(nc, "rho", FT, cspace, ("time",))
-    thermo_var = if :ρe_tot in propertynames(Y.c)
-        "e_tot"
-    elseif :ρθ in propertynames(Y.c)
-        "theta"
-    else
-        error("Unfound thermodynamic variable")
-    end
-    nc_thermo = defVar(nc, thermo_var, FT, cspace, ("time",))
+    nc_thermo = defVar(nc, "e_tot", FT, cspace, ("time",))
     nc_u = defVar(nc, "u", FT, cspace, ("time",))
     nc_v = defVar(nc, "v", FT, cspace, ("time",))
     nc_w = defVar(nc, "w", FT, cspace, ("time",))
@@ -199,11 +192,7 @@ function remap2latlon(filein, data_dir, remap_tmpdir, weightfile, nlat, nlon)
     # density
     nc_rho[:, 1] = Y.c.ρ
     # thermodynamics
-    if :ρe_tot in propertynames(Y.c)
-        nc_thermo[:, 1] = Y.c.ρe_tot ./ Y.c.ρ
-    elseif :ρθ in propertynames(Y.c)
-        nc_thermo[:, 1] = Y.c.ρθ ./ Y.c.ρ
-    end
+    nc_thermo[:, 1] = Y.c.ρe_tot ./ Y.c.ρ
     # physical horizontal velocity
     uh_phy = Geometry.transform.(tuple(Geometry.UVAxis()), Y.c.uₕ)
     nc_u[:, 1] = uh_phy.components.data.:1
@@ -279,7 +268,7 @@ function remap2latlon(filein, data_dir, remap_tmpdir, weightfile, nlat, nlon)
         joinpath(out_dir, first(splitext(basename(filein))) * ".nc")
     dry_variables = [
         "rho",
-        thermo_var,
+        "e_tot",
         "u",
         "v",
         "w",

src/cache/precomputed_quantities.jl

@@ -16,8 +16,6 @@ Allocates and returns the precomputed quantities:
     - `ᶜts`: the thermodynamic state on cell centers
     - `ᶜp`: the air pressure on cell centers
     - `sfc_conditions`: the conditions at the surface (at the bottom cell faces)
-
-If the `energy_form` is TotalEnergy, there is an additional quantity:
     - `ᶜh_tot`: the total enthalpy on cell centers
 
 If the `turbconv_model` is EDMFX, there also two SGS versions of every quantity
@@ -36,14 +34,10 @@ TODO: Rename `ᶜK` to `ᶜκ`.
 """
 function precomputed_quantities(Y, atmos)
     FT = eltype(Y)
-    @assert (
-        !(atmos.moisture_model isa DryModel) &&
-        atmos.energy_form isa TotalEnergy
-    ) || !(atmos.turbconv_model isa DiagnosticEDMFX)
-    @assert (
-        !(atmos.moisture_model isa DryModel) &&
-        atmos.energy_form isa TotalEnergy
-    ) || !(atmos.turbconv_model isa PrognosticEDMFX)
+    @assert !(atmos.moisture_model isa DryModel) ||
+            !(atmos.turbconv_model isa DiagnosticEDMFX)
+    @assert !(atmos.moisture_model isa DryModel) ||
+            !(atmos.turbconv_model isa PrognosticEDMFX)
     @assert !(atmos.edmfx_detr_model isa ConstantAreaDetrainment) ||
             !(atmos.turbconv_model isa DiagnosticEDMFX)
     TST = thermo_state_type(atmos.moisture_model, FT)
@@ -56,10 +50,7 @@ function precomputed_quantities(Y, atmos)
         ᶜK = similar(Y.c, FT),
         ᶜts = similar(Y.c, TST),
         ᶜp = similar(Y.c, FT),
-        (
-            atmos.energy_form isa TotalEnergy ?
-            (; ᶜh_tot = similar(Y.c, FT)) : (;)
-        )...,
+        ᶜh_tot = similar(Y.c, FT),
         sfc_conditions = Fields.Field(SCT, Spaces.level(axes(Y.f), half)),
     )
     advective_sgs_quantities =
@@ -252,12 +243,8 @@ function thermo_state(
     return get_ts(ρ, p, θ, e_int, q_tot, q_pt)
 end
 
-function thermo_vars(energy_form, moisture_model, specific, K, Φ)
-    energy_var = if energy_form isa PotentialTemperature
-        (; specific.θ)
-    elseif energy_form isa TotalEnergy
-        (; e_int = specific.e_tot - K - Φ)
-    end
+function thermo_vars(moisture_model, specific, K, Φ)
+    energy_var = (; e_int = specific.e_tot - K - Φ)
     moisture_var = if moisture_model isa DryModel
         (;)
     elseif moisture_model isa EquilMoistModel
@@ -269,19 +256,17 @@ function thermo_vars(energy_form, moisture_model, specific, K, Φ)
     return (; energy_var..., moisture_var...)
 end
 
-ts_gs(thermo_params, energy_form, moisture_model, specific, K, Φ, ρ) =
-    thermo_state(
-        thermo_params;
-        thermo_vars(energy_form, moisture_model, specific, K, Φ)...,
-        ρ,
-    )
+ts_gs(thermo_params, moisture_model, specific, K, Φ, ρ) = thermo_state(
+    thermo_params;
+    thermo_vars(moisture_model, specific, K, Φ)...,
+    ρ,
+)
 
-ts_sgs(thermo_params, energy_form, moisture_model, specific, K, Φ, p) =
-    thermo_state(
-        thermo_params;
-        thermo_vars(energy_form, moisture_model, specific, K, Φ)...,
-        p,
-    )
+ts_sgs(thermo_params, moisture_model, specific, K, Φ, p) = thermo_state(
+    thermo_params;
+    thermo_vars(moisture_model, specific, K, Φ)...,
+    p,
+)
 
 """
     set_precomputed_quantities!(Y, p, t)
@@ -300,10 +285,10 @@ function instead of recomputing the value yourself. Otherwise, it will be
 difficult to ensure that the duplicated computations are consistent.
 """
 NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
-    (; energy_form, moisture_model, turbconv_model) = p.atmos
+    (; moisture_model, turbconv_model) = p.atmos
     thermo_params = CAP.thermodynamics_params(p.params)
     n = n_mass_flux_subdomains(turbconv_model)
-    thermo_args = (thermo_params, energy_form, moisture_model)
+    thermo_args = (thermo_params, moisture_model)
     (; ᶜΦ) = p.core
     (; ᶜspecific, ᶜu, ᶠu³, ᶜK, ᶜts, ᶜp) = p.precomputed
     ᶠuₕ³ = p.scratch.ᶠtemp_CT3
@@ -323,7 +308,7 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
         # quantities of the form ᶜρaχ⁰ / ᶜρ and ᶜρaχʲ / ᶜρ. However, we cannot
         # compute ᶜρ⁰ and ᶜρʲ without first computing ᶜts⁰ and ᶜtsʲ, both of
         # which depend on the value of ᶜp, which in turn depends on ᶜts. Since
-        # ᶜts depends on ᶜK (at least when the energy_form is TotalEnergy), this
+        # ᶜts depends on ᶜK, this
         # means that the amount by which ᶜK needs to be incremented is a
         # function of ᶜK itself. So, unless we run a nonlinear solver here, this
         # circular dependency will prevent us from computing the exact value of
@@ -335,11 +320,8 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
     @. ᶜts = ts_gs(thermo_args..., ᶜspecific, ᶜK, ᶜΦ, Y.c.ρ)
     @. ᶜp = TD.air_pressure(thermo_params, ᶜts)
 
-    if energy_form isa TotalEnergy
-        (; ᶜh_tot) = p.precomputed
-        @. ᶜh_tot =
-            TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜspecific.e_tot)
-    end
+    (; ᶜh_tot) = p.precomputed
+    @. ᶜh_tot = TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜspecific.e_tot)
 
     if !isnothing(p.sfc_setup)
         SurfaceConditions.update_surface_conditions!(Y, p, t)

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -50,9 +50,8 @@ Updates the draft thermo state and boundary conditions
 precomputed quantities stored in `p` for edmfx.
 """
 function set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ³, t)
-    (; energy_form, moisture_model, turbconv_model) = p.atmos
+    (; moisture_model, turbconv_model) = p.atmos
     #EDMFX BCs only support total energy as state variable
-    @assert energy_form isa TotalEnergy
     @assert !(moisture_model isa DryModel)
 
     FT = Spaces.undertype(axes(Y.c))

src/callbacks/callbacks.jl

@@ -326,31 +326,26 @@ NVTX.@annotate function compute_diagnostics(integrator)
         turbulence_convection_diagnostic = NamedTuple()
     end
 
-    if p.atmos.energy_form isa TotalEnergy
-        sfc_local_geometry =
-            Fields.level(Fields.local_geometry_field(Y.f), Fields.half)
-        surface_ct3_unit =
-            CT3.(unit_basis_vector_data.(CT3, sfc_local_geometry))
-        (; ρ_flux_uₕ, ρ_flux_h_tot) = p.precomputed.sfc_conditions
-        sfc_flux_momentum =
-            Geometry.UVVector.(
-                adjoint.(ρ_flux_uₕ ./ Spaces.level(ᶠinterp.(Y.c.ρ), half)) .*
-                surface_ct3_unit
-            )
+    sfc_local_geometry =
+        Fields.level(Fields.local_geometry_field(Y.f), Fields.half)
+    surface_ct3_unit = CT3.(unit_basis_vector_data.(CT3, sfc_local_geometry))
+    (; ρ_flux_uₕ, ρ_flux_h_tot) = p.precomputed.sfc_conditions
+    sfc_flux_momentum =
+        Geometry.UVVector.(
+            adjoint.(ρ_flux_uₕ ./ Spaces.level(ᶠinterp.(Y.c.ρ), half)) .*
+            surface_ct3_unit
+        )
+    vert_diff_diagnostic = (;
+        sfc_flux_u = sfc_flux_momentum.components.data.:1,
+        sfc_flux_v = sfc_flux_momentum.components.data.:2,
+        sfc_flux_energy = dot.(ρ_flux_h_tot, surface_ct3_unit),
+    )
+    if :ρq_tot in propertynames(Y.c)
+        (; ρ_flux_q_tot) = p.precomputed.sfc_conditions
         vert_diff_diagnostic = (;
-            sfc_flux_u = sfc_flux_momentum.components.data.:1,
-            sfc_flux_v = sfc_flux_momentum.components.data.:2,
-            sfc_flux_energy = dot.(ρ_flux_h_tot, surface_ct3_unit),
+            vert_diff_diagnostic...,
+            sfc_evaporation = dot.(ρ_flux_q_tot, surface_ct3_unit),
         )
-        if :ρq_tot in propertynames(Y.c)
-            (; ρ_flux_q_tot) = p.precomputed.sfc_conditions
-            vert_diff_diagnostic = (;
-                vert_diff_diagnostic...,
-                sfc_evaporation = dot.(ρ_flux_q_tot, surface_ct3_unit),
-            )
-        end
-    else
-        vert_diff_diagnostic = NamedTuple()
     end
 
     if p.atmos.radiation_mode isa RRTMGPI.AbstractRRTMGPMode

src/diagnostics/Diagnostics.jl

@@ -19,9 +19,6 @@ import ..DryModel
 import ..EquilMoistModel
 import ..NonEquilMoistModel
 
-# energy_form
-import ..TotalEnergy
-
 # precip_model
 import ..Microphysics0Moment
 

src/diagnostics/core_diagnostics.jl

@@ -418,10 +418,7 @@ add_diagnostic_variable!(
 ###
 # Eastward and northward surface drag component (2d)
 ###
-compute_tau!(_, _, _, _, energy_form::T) where {T} =
-    error_diagnostic_variable("tau", energy_form)
-
-function compute_tau!(out, state, cache, component, energy_form::TotalEnergy)
+function compute_tau!(out, state, cache, component)
     sfc_local_geometry =
         Fields.level(Fields.local_geometry_field(state.f), Fields.half)
     surface_ct3_unit = CT3.(unit_basis_vector_data.(CT3, sfc_local_geometry))
@@ -446,10 +443,8 @@ function compute_tau!(out, state, cache, component, energy_form::TotalEnergy)
     return
 end
 
-compute_tauu!(out, state, cache, time) =
-    compute_tau!(out, state, cache, :1, cache.atmos.energy_form)
-compute_tauv!(out, state, cache, time) =
-    compute_tau!(out, state, cache, :2, cache.atmos.energy_form)
+compute_tauu!(out, state, cache, time) = compute_tau!(out, state, cache, :1)
+compute_tauv!(out, state, cache, time) = compute_tau!(out, state, cache, :2)
 
 add_diagnostic_variable!(
     short_name = "tauu",
@@ -472,12 +467,7 @@ add_diagnostic_variable!(
 ###
 # Surface energy flux (2d) - TODO: this may need to be split into sensible and latent heat fluxes
 ###
-compute_hfes!(out, state, cache, time) =
-    compute_hfes!(out, state, cache, time, cache.atmos.energy_form)
-compute_hfes!(_, _, _, _, energy_form::T) where {T} =
-    error_diagnostic_variable("hfes", energy_form)
-
-function compute_hfes!(out, state, cache, time, energy_form::TotalEnergy)
+function compute_hfes!(out, state, cache, time)
     (; ρ_flux_h_tot) = cache.precomputed.sfc_conditions
     sfc_local_geometry =
         Fields.level(Fields.local_geometry_field(state.f), Fields.half)
@@ -500,32 +490,17 @@ add_diagnostic_variable!(
 ###
 # Surface evaporation (2d)
 ###
-compute_evspsbl!(out, state, cache, time) = compute_evspsbl!(
-    out,
-    state,
-    cache,
-    time,
-    cache.atmos.moisture_model,
-    cache.atmos.energy_form,
-)
-compute_evspsbl!(
-    _,
-    _,
-    _,
-    _,
-    moisture_model::T1,
-    energy_form::T2,
-) where {T1, T2} = error_diagnostic_variable(
-    "Can only compute surface_evaporation with energy_form = TotalEnergy() and with a moist model",
-)
+compute_evspsbl!(out, state, cache, time) =
+    compute_evspsbl!(out, state, cache, time, cache.atmos.moisture_model)
+compute_evspsbl!(_, _, _, _, model::T) where {T} =
+    error_diagnostic_variable("evspsbl", model)
 
 function compute_evspsbl!(
     out,
     state,
     cache,
     time,
     moisture_model::T,
-    energy_form::TotalEnergy,
 ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}}
     (; ρ_flux_q_tot) = cache.precomputed.sfc_conditions
     sfc_local_geometry =

src/diagnostics/default_diagnostics.jl

@@ -193,7 +193,7 @@ hourly_average(short_names; output_writer = HDF5Writer()) =
 ########
 function core_default_diagnostics()
     core_diagnostics =
-        ["ts", "ta", "thetaa", "ha", "pfull", "rhoa", "ua", "va", "wa"]
+        ["ts", "ta", "thetaa", "ha", "pfull", "rhoa", "ua", "va", "wa", "hfes"]
 
     return [
         daily_averages(core_diagnostics...)...,
@@ -202,16 +202,6 @@ function core_default_diagnostics()
     ]
 end
 
-###############
-# Energy form #
-###############
-function default_diagnostics(::TotalEnergy)
-    total_energy_diagnostics = ["hfes"]
-
-    return [daily_averages(total_energy_diagnostics...)...]
-end
-
-
 ##################
 # Moisture model #
 ##################

src/initial_conditions/InitialConditions.jl

@@ -1,8 +1,6 @@
 module InitialConditions
 
 import ..AtmosModel
-import ..PotentialTemperature
-import ..TotalEnergy
 import ..DryModel
 import ..EquilMoistModel
 import ..NonEquilMoistModel