change prognostic variables from (ρθ_liq_ice,ρq_tot) to (ρe_tot,ρq_tot)

WIP

.buildkite/pipeline.yml

@@ -154,12 +154,12 @@ steps:
         command: "julia --color=yes --project=integration_tests integration_tests/driver.jl --case Rico --thermo_covariance_model prognostic --skip_tests true --suffix _prognostic_covarinaces"
         artifact_paths: "Output.Rico.01_prognostic_covarinaces/stats/comparison/*"
 
-  - group: "3D Experiments"
-    steps:
+  # - group: "3D Experiments"
+  #   steps:
 
-      - label: ":package: 3D Bomex"
-        command: "julia --color=yes --project=integration_tests integration_tests/3dBomex.jl"
-        artifact_paths: "integration_tests/output/3dBomex/*"
+  #     - label: ":package: 3D Bomex"
+  #       command: "julia --color=yes --project=integration_tests integration_tests/3dBomex.jl"
+  #       artifact_paths: "integration_tests/output/3dBomex/*"
 
   - group: "Performance monitoring"
     steps:

driver/Forcing.jl

@@ -9,19 +9,19 @@ function initialize(self::TC.ForcingBase{TC.ForcingLES}, grid, state, LESDat::TC
         zc_les = Array(TC.get_nc_data(data, "zc"))
 
         dTdt_hadv = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "dtdt_hadv", imin, imax))
-        H_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "thetali_mean", imin, imax))
+        T_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "temperature_mean", imin, imax))
         dTdt_fluc = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "dtdt_fluc", imin, imax))
         dqtdt_hadv = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "dqtdt_hadv", imin, imax))
         qt_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "qt_mean", imin, imax))
         dqtdt_fluc = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "dqtdt_fluc", imin, imax))
         subsidence = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "ls_subsidence", imin, imax))
         u_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "u_mean", imin, imax))
         v_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "v_mean", imin, imax))
-        (; dTdt_hadv, H_nudge, dTdt_fluc, dqtdt_hadv, qt_nudge, dqtdt_fluc, subsidence, u_nudge, v_nudge)
+        (; dTdt_hadv, T_nudge, dTdt_fluc, dqtdt_hadv, qt_nudge, dqtdt_fluc, subsidence, u_nudge, v_nudge)
     end
     for k in TC.real_center_indices(grid)
         aux_gm.dTdt_hadv[k] = nt.dTdt_hadv[k]
-        aux_gm.H_nudge[k] = nt.H_nudge[k]
+        aux_gm.T_nudge[k] = nt.T_nudge[k]
         aux_gm.dTdt_fluc[k] = nt.dTdt_fluc[k]
         aux_gm.dqtdt_hadv[k] = nt.dqtdt_hadv[k]
         aux_gm.qt_nudge[k] = nt.qt_nudge[k]

driver/Surface.jl

@@ -68,7 +68,7 @@ function get_surface(
         ch = result.Ch,
         ρu_flux = surf_params.zero_uv_fluxes ? FT(0) : result.ρτxz,
         ρv_flux = surf_params.zero_uv_fluxes ? FT(0) : result.ρτyz,
-        ρθ_liq_ice_flux = shf / TD.cp_m(param_set, ts_in),
+        ρe_tot_flux = shf + lhf,
         ρq_tot_flux = lhf / TD.latent_heat_vapor(param_set, ts_in),
         wstar = convective_vel,
         ρq_liq_flux = FT(0),
@@ -126,7 +126,7 @@ function get_surface(
         ustar = result.ustar,
         ρu_flux = result.ρτxz,
         ρv_flux = result.ρτyz,
-        ρθ_liq_ice_flux = shf / TD.cp_m(param_set, ts_in),
+        ρe_tot_flux = shf + lhf,
         ρq_tot_flux = lhf / TD.latent_heat_vapor(param_set, ts_in),
         bflux = result.buoy_flux,
         wstar = convective_vel,
@@ -183,7 +183,7 @@ function get_surface(
         ustar = result.ustar,
         ρu_flux = result.ρτxz,
         ρv_flux = result.ρτyz,
-        ρθ_liq_ice_flux = shf / TD.cp_m(param_set, ts_in),
+        ρe_tot_flux = shf + lhf,
         ρq_tot_flux = lhf / TD.latent_heat_vapor(param_set, ts_in),
         bflux = result.buoy_flux,
         wstar = convective_vel,
@@ -245,7 +245,7 @@ function get_surface(
         ustar = result.ustar,
         ρu_flux = result.ρτxz,
         ρv_flux = result.ρτyz,
-        ρθ_liq_ice_flux = shf / TD.cp_m(param_set, ts_in),
+        ρe_tot_flux = shf + lhf,
         ρq_tot_flux = lhf / TD.latent_heat_vapor(param_set, ts_in),
         bflux = result.buoy_flux,
         wstar = convective_vel,

driver/compute_diagnostics.jl

@@ -139,8 +139,6 @@ function compute_diagnostics!(
 
     @inbounds for k in TC.real_center_indices(grid)
         aux_gm.s[k] = TD.specific_entropy(param_set, ts_gm[k])
-        aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k]
-        aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k]
         aux_en.s[k] = TD.specific_entropy(param_set, ts_en[k])
     end
     @inbounds for k in TC.real_center_indices(grid)

driver/dycore.jl

@@ -62,15 +62,15 @@ cent_aux_vars_gm(FT, edmf) = (;
     subsidence = FT(0), #Large-scale subsidence
     dTdt_hadv = FT(0), #Horizontal advection of temperature
     dqtdt_hadv = FT(0), #Horizontal advection of moisture
-    H_nudge = FT(0), #Reference H profile for relaxation tendency
+    T_nudge = FT(0), #Reference T profile for relaxation tendency
     qt_nudge = FT(0), #Reference qt profile for relaxation tendency
     dTdt_fluc = FT(0), #Vertical turbulent advection of temperature
     dqtdt_fluc = FT(0), #Vertical turbulent advection of moisture
     u_nudge = FT(0), #Reference u profile for relaxation tendency
     v_nudge = FT(0), #Reference v profile for relaxation tendency
     ug = FT(0), #Geostrophic u velocity
     vg = FT(0), #Geostrophic v velocity
-    ∇θ_liq_ice_gm = FT(0),
+    ∇MSE_gm = FT(0),
     ∇q_tot_gm = FT(0),
     cent_aux_vars_gm_moisture(FT, edmf.moisture_model)...,
     θ_virt = FT(0),
@@ -89,7 +89,7 @@ face_aux_vars_gm(FT, edmf) = (;
     diffusive_flux_s = FT(0),
     total_flux_s = FT(0),
     f_rad = FT(0),
-    sgs_flux_θ_liq_ice = FT(0),
+    sgs_flux_h_tot = FT(0),
     sgs_flux_q_tot = FT(0),
     face_aux_vars_gm_moisture(FT, edmf.moisture_model)...,
     sgs_flux_u = FT(0),
@@ -139,7 +139,7 @@ cent_prognostic_vars_gm(::Type{FT}, local_geometry, edmf) where {FT} = (;
     ρ = FT(0),
     u = FT(0),
     v = FT(0),
-    ρθ_liq_ice = FT(0),
+    ρe_tot = FT(0),
     ρq_tot = FT(0),
     # TODO: Change to:
     # uₕ = CCG.Covariant12Vector(CCG.UVVector(FT(0), FT(0)), local_geometry),
@@ -221,7 +221,6 @@ function compute_ref_state!(
     @info "z_span = $z_span"
     prob = ODE.ODEProblem(rhs, logp, z_span)
     sol = ODE.solve(prob, ODE.Tsit5(), reltol = 1e-12, abstol = 1e-12)
-
     parent(p_f) .= sol.(vec(grid.zf))
     parent(p_c) .= sol.(vec(grid.zc))
 
@@ -241,23 +240,43 @@ function compute_ref_state!(
     return nothing
 end
 
-function set_prog_from_aux!(state, ::TC.ρθVar)
+
+function set_thermo_state_peq!(state, grid, moisture_model, param_set)
+    Ic = CCO.InterpolateF2C()
+    FT = eltype(grid)
+    ts_gm = TC.center_aux_grid_mean(state).ts
     prog_gm = TC.center_prog_grid_mean(state)
+    prog_gm_f = TC.face_prog_grid_mean(state)
     aux_gm = TC.center_aux_grid_mean(state)
-    @. prog_gm.ρθ_liq_ice = prog_gm.ρ * aux_gm.θ_liq_ice
-    @. prog_gm.ρq_tot .= prog_gm.ρ * aux_gm.q_tot
+    p_c = aux_gm.p
+    ρ_c = prog_gm.ρ
+    w_c = copy(prog_gm.ρe_tot)
+    @. w_c = Ic(FT(0) + prog_gm_f.w)
+    @inbounds for k in TC.real_center_indices(grid)
+        thermo_args = if moisture_model isa TC.EquilibriumMoisture
+            ()
+        elseif moisture_model isa TC.NonEquilibriumMoisture
+            (prog_gm.q_liq[k], prog_gm.q_ice[k])
+        else
+            error("Something went wrong. The moisture_model options are equilibrium or nonequilibrium")
+        end
+        e_kin = TC.kinetic_energy(prog_gm.u[k], prog_gm.v[k], w_c[k])
+        e_pot = TC.geopotential(param_set, grid.zc.z[k])
+        e_int = prog_gm.ρe_tot[k] / ρ_c[k] - e_kin - e_pot
+        ts_gm[k] = TC.thermo_state_peq(param_set, p_c[k], e_int, aux_gm.q_tot[k], thermo_args...)
+        aux_gm.θ_liq_ice[k] = TD.liquid_ice_pottemp(param_set, ts_gm[k])
+        aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k]
+    end
     return nothing
 end
 
-function set_thermo_state!(state, grid, moisture_model, param_set)
+function set_thermo_state_pθq!(state, grid, moisture_model, param_set)
+    Ic = CCO.InterpolateF2C()
     ts_gm = TC.center_aux_grid_mean(state).ts
     prog_gm = TC.center_prog_grid_mean(state)
     aux_gm = TC.center_aux_grid_mean(state)
     p_c = aux_gm.p
-    ρ_c = prog_gm.ρ
     @inbounds for k in TC.real_center_indices(grid)
-        aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k]
-        aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k]
         thermo_args = if moisture_model isa TC.EquilibriumMoisture
             ()
         elseif moisture_model isa TC.NonEquilibriumMoisture
@@ -270,14 +289,31 @@ function set_thermo_state!(state, grid, moisture_model, param_set)
     return nothing
 end
 
+function set_grid_mean_from_thermo_state!(param_set, state, grid)
+    Ic = CCO.InterpolateF2C()
+    ts_gm = TC.center_aux_grid_mean(state).ts
+    prog_gm = TC.center_prog_grid_mean(state)
+    prog_gm_f = TC.face_prog_grid_mean(state)
+    aux_gm = TC.center_aux_grid_mean(state)
+    ρ_c = prog_gm.ρ
+    FT = eltype(grid)
+    w_c = copy(prog_gm.ρe_tot)
+    @. w_c = Ic(FT(0) + prog_gm_f.w)
+    @inbounds for k in TC.real_center_indices(grid)
+        e_kin = TC.kinetic_energy(prog_gm.u[k], prog_gm.v[k], w_c[k])
+        e_pot = TC.geopotential(param_set, grid.zc.z[k])
+        prog_gm.ρe_tot[k] = ρ_c[k] * TD.total_energy(param_set, ts_gm[k], e_kin, e_pot)
+        prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
+    end
+    return nothing
+end
+
 function assign_thermo_aux!(state, grid, moisture_model, param_set)
     aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
     ts_gm = TC.center_aux_grid_mean(state).ts
     ρ_c = prog_gm.ρ
     @inbounds for k in TC.real_center_indices(grid)
-        aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k]
-        aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k]
         ts = ts_gm[k]
         aux_gm.q_liq[k] = TD.liquid_specific_humidity(param_set, ts)
         aux_gm.q_ice[k] = TD.ice_specific_humidity(param_set, ts)
@@ -311,7 +347,7 @@ function ∑tendencies!(tendencies::FV, prog::FV, params::NT, t::Real) where {NT
 
     state = TC.State(prog, aux, tendencies)
 
-    set_thermo_state!(state, grid, edmf.moisture_model, param_set)
+    set_thermo_state_peq!(state, grid, edmf.moisture_model, param_set)
 
     # TODO: where should this live?
     aux_gm = TC.center_aux_grid_mean(state)
@@ -365,14 +401,19 @@ function compute_gm_tendencies!(
     force::TC.ForcingBase,
     param_set::APS,
 )
+    Ic = CCO.InterpolateF2C()
+    R_d = CPP.R_d(param_set)
+    T_0 = CPP.T_0(param_set)
+    Lv_0 = CPP.LH_v0(param_set)
     tendencies_gm = TC.center_tendencies_grid_mean(state)
     kc_toa = TC.kc_top_of_atmos(grid)
     kf_surf = TC.kf_surface(grid)
     FT = eltype(grid)
     prog_gm = TC.center_prog_grid_mean(state)
+    prog_gm_f = TC.face_prog_grid_mean(state)
     aux_gm = TC.center_aux_grid_mean(state)
     aux_gm_f = TC.face_aux_grid_mean(state)
-    ∇θ_liq_ice_gm = TC.center_aux_grid_mean(state).∇θ_liq_ice_gm
+    ∇MSE_gm = TC.center_aux_grid_mean(state).∇MSE_gm
     ∇q_tot_gm = TC.center_aux_grid_mean(state).∇q_tot_gm
     aux_en = TC.center_aux_environment(state)
     aux_en_f = TC.face_aux_environment(state)
@@ -384,14 +425,20 @@ function compute_gm_tendencies!(
     aux_tc = TC.center_aux_turbconv(state)
     ts_gm = TC.center_aux_grid_mean(state).ts
 
-    θ_liq_ice_gm_toa = aux_gm.θ_liq_ice[kc_toa]
-    q_tot_gm_toa = aux_gm.q_tot[kc_toa]
-    RBθ = CCO.RightBiasedC2F(; top = CCO.SetValue(θ_liq_ice_gm_toa))
+    e_kin = copy(prog_gm.ρe_tot)
+    w_c = copy(prog_gm.ρe_tot)
+    h_tot_gm = copy(prog_gm.ρe_tot)
+    @. w_c = Ic(FT(0) + prog_gm_f.w)
+    @. h_tot_gm = TC.anelastic_total_enthalpy(param_set, prog_gm.ρe_tot / ρ_c, ts_gm)
+    @. e_kin = TC.kinetic_energy(prog_gm.u, prog_gm.v, w_c)
+    MSE_gm_toa = h_tot_gm[kc_toa] - e_kin[kc_toa]
+    q_tot_gm_toa = prog_gm.ρq_tot[kc_toa] / ρ_c[kc_toa]
+    RBe = CCO.RightBiasedC2F(; top = CCO.SetValue(MSE_gm_toa))
     RBq = CCO.RightBiasedC2F(; top = CCO.SetValue(q_tot_gm_toa))
     wvec = CC.Geometry.WVector
     ∇c = CCO.DivergenceF2C()
-    @. ∇θ_liq_ice_gm = ∇c(wvec(RBθ(aux_gm.θ_liq_ice)))
-    @. ∇q_tot_gm = ∇c(wvec(RBq(aux_gm.q_tot)))
+    @. ∇MSE_gm = ∇c(wvec(RBe(h_tot_gm - e_kin)))
+    @. ∇q_tot_gm = ∇c(wvec(RBq(prog_gm.ρq_tot / ρ_c)))
 
     if edmf.moisture_model isa TC.NonEquilibriumMoisture
         ∇q_liq_gm = TC.center_aux_grid_mean(state).∇q_liq_gm
@@ -406,28 +453,33 @@ function compute_gm_tendencies!(
 
     # Apply forcing and radiation
     @inbounds for k in TC.real_center_indices(grid)
-        Π = TD.exner(param_set, ts_gm[k])
+        cp_m = TD.cp_m(param_set, ts_gm[k])
+        cp_v = CPP.cp_v(param_set)
+        cv_v = CPP.cv_v(param_set)
+        R_v = CPP.R_v(param_set)
+        cv_m = TD.cv_m(param_set, ts_gm[k])
+        h_v = cv_v * (aux_gm.T[k] - T_0) + Lv_0 - R_v * T_0
 
         # Coriolis
         if force.apply_coriolis
             tendencies_gm.u[k] -= force.coriolis_param * (aux_gm.vg[k] - prog_gm.v[k])
             tendencies_gm.v[k] += force.coriolis_param * (aux_gm.ug[k] - prog_gm.u[k])
         end
         # LS Subsidence
-        tendencies_gm.ρθ_liq_ice[k] -= ρ_c[k] * ∇θ_liq_ice_gm[k] * aux_gm.subsidence[k]
-        tendencies_gm.ρq_tot[k] -= ρ_c[k] * ∇q_tot_gm[k] * aux_gm.subsidence[k]
+        tendencies_gm.ρe_tot[k] -= ρ_c[k] * aux_gm.subsidence[k] * ∇MSE_gm[k]
+        tendencies_gm.ρq_tot[k] -= ρ_c[k] * aux_gm.subsidence[k] * ∇q_tot_gm[k]
         if edmf.moisture_model isa TC.NonEquilibriumMoisture
             tendencies_gm.q_liq[k] -= ∇q_liq_gm[k] * aux_gm.subsidence[k]
             tendencies_gm.q_ice[k] -= ∇q_ice_gm[k] * aux_gm.subsidence[k]
         end
         # Radiation
         if TC.rad_type(radiation) <: Union{TC.RadiationDYCOMS_RF01, TC.RadiationLES, TC.RadiationTRMM_LBA}
-            tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * aux_gm.dTdt_rad[k] / Π
+            tendencies_gm.ρe_tot[k] += ρ_c[k] * cv_m * aux_gm.dTdt_rad[k]
         end
         # LS advection
         tendencies_gm.ρq_tot[k] += ρ_c[k] * aux_gm.dqtdt_hadv[k]
         if !(TC.force_type(force) <: TC.ForcingDYCOMS_RF01)
-            tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * aux_gm.dTdt_hadv[k] / Π
+            tendencies_gm.ρe_tot[k] += ρ_c[k] * (cp_m * aux_gm.dTdt_hadv[k] + h_v * aux_gm.dqtdt_hadv[k])
         end
         if edmf.moisture_model isa TC.NonEquilibriumMoisture
             tendencies_gm.q_liq[k] += aux_gm.dqldt[k]
@@ -436,21 +488,22 @@ function compute_gm_tendencies!(
 
         # LES specific forcings
         if TC.force_type(force) <: TC.ForcingLES
-            H_fluc = aux_gm.dTdt_fluc[k] / Π
+            T_fluc = aux_gm.dTdt_fluc[k]
 
             gm_U_nudge_k = (aux_gm.u_nudge[k] - prog_gm.u[k]) / force.wind_nudge_τᵣ
             gm_V_nudge_k = (aux_gm.v_nudge[k] - prog_gm.v[k]) / force.wind_nudge_τᵣ
 
             Γᵣ = TC.compute_les_Γᵣ(grid.zc[k], force.scalar_nudge_τᵣ, force.scalar_nudge_zᵢ, force.scalar_nudge_zᵣ)
-            gm_H_nudge_k = Γᵣ * (aux_gm.H_nudge[k] - aux_gm.θ_liq_ice[k])
+            gm_T_nudge_k = Γᵣ * (aux_gm.T_nudge[k] - aux_gm.T[k])
             gm_q_tot_nudge_k = Γᵣ * (aux_gm.qt_nudge[k] - aux_gm.q_tot[k])
             if edmf.moisture_model isa TC.NonEquilibriumMoisture
                 gm_q_liq_nudge_k = Γᵣ * (aux_gm.ql_nudge[k] - prog_gm.q_liq[k])
                 gm_q_ice_nudge_k = Γᵣ * (aux_gm.qi_nudge[k] - prog_gm.q_ice[k])
             end
 
-            tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * (gm_H_nudge_k + H_fluc)
             tendencies_gm.ρq_tot[k] += ρ_c[k] * (gm_q_tot_nudge_k + aux_gm.dqtdt_fluc[k])
+            tendencies_gm.ρe_tot[k] +=
+                ρ_c[k] * (cv_m * (gm_T_nudge_k + T_fluc) + h_v * (gm_q_tot_nudge_k + aux_gm.dqtdt_fluc[k]))
             tendencies_gm.u[k] += gm_U_nudge_k
             tendencies_gm.v[k] += gm_V_nudge_k
             if edmf.moisture_model isa TC.NonEquilibriumMoisture
@@ -466,38 +519,38 @@ function compute_gm_tendencies!(
                 aux_en.qt_tendency_precip_formation[k] +
                 aux_tc.qt_tendency_precip_sinks[k]
             )
-        tendencies_gm.ρθ_liq_ice[k] +=
+
+        tendencies_gm.ρe_tot[k] +=
             ρ_c[k] * (
-                aux_bulk.θ_liq_ice_tendency_precip_formation[k] +
-                aux_en.θ_liq_ice_tendency_precip_formation[k] +
-                aux_tc.θ_liq_ice_tendency_precip_sinks[k]
+                aux_bulk.e_tot_tendency_precip_formation[k] +
+                aux_en.e_tot_tendency_precip_formation[k] +
+                aux_tc.e_tot_tendency_precip_sinks[k]
             )
+
         if edmf.moisture_model isa TC.NonEquilibriumMoisture
             tendencies_gm.q_liq[k] += aux_bulk.ql_tendency_precip_formation[k] + aux_en.ql_tendency_precip_formation[k]
             tendencies_gm.q_ice[k] += aux_bulk.qi_tendency_precip_formation[k] + aux_en.qi_tendency_precip_formation[k]
         end
     end
 
-    TC.compute_sgs_flux!(edmf, grid, state, surf)
-
-    sgs_flux_θ_liq_ice = aux_gm_f.sgs_flux_θ_liq_ice
+    TC.compute_sgs_flux!(edmf, grid, state, surf, param_set)
+    sgs_flux_h_tot = aux_gm_f.sgs_flux_h_tot
     sgs_flux_q_tot = aux_gm_f.sgs_flux_q_tot
     sgs_flux_u = aux_gm_f.sgs_flux_u
     sgs_flux_v = aux_gm_f.sgs_flux_v
     # apply surface BC as SGS flux at lowest level
-    sgs_flux_θ_liq_ice[kf_surf] = surf.ρθ_liq_ice_flux
+    sgs_flux_h_tot[kf_surf] = surf.ρe_tot_flux
     sgs_flux_q_tot[kf_surf] = surf.ρq_tot_flux
     sgs_flux_u[kf_surf] = surf.ρu_flux
     sgs_flux_v[kf_surf] = surf.ρv_flux
 
-    tends_ρθ_liq_ice = tendencies_gm.ρθ_liq_ice
+    tends_ρe_tot = tendencies_gm.ρe_tot
     tends_ρq_tot = tendencies_gm.ρq_tot
     tends_u = tendencies_gm.u
     tends_v = tendencies_gm.v
 
     ∇sgs = CCO.DivergenceF2C()
-
-    @. tends_ρθ_liq_ice += -∇sgs(wvec(sgs_flux_θ_liq_ice))
+    @. tends_ρe_tot += -∇sgs(wvec(sgs_flux_h_tot))
     @. tends_ρq_tot += -∇sgs(wvec(sgs_flux_q_tot))
     @. tends_u += -∇sgs(wvec(sgs_flux_u)) / ρ_c
     @. tends_v += -∇sgs(wvec(sgs_flux_v)) / ρ_c