Merge #1220

1220: change grid mean variable to theta r=yairchn a=yairchn

Reverting back to theta ...

Co-authored-by: yairchn <yairchn@caltech.edu>

README.md

@@ -128,7 +128,7 @@ Upon running a particular experiment (described above), comparison plots (agains
 | `u`           | [✓]           |  [x]            |  [x]         |
 | `v`           | [✓]           |  [x]            |  [x]         |
 | `w`           | [-]           |  [x]            |  [✓]         |
-| `e_tot`       | [✓]           |  [x]            |  [✓]         |
+| `θ_liq_ice`   | [✓]           |  [x]            |  [✓]         |
 | `q_tot`       | [✓]           |  [x]            |  [✓]         |
 | `a`           | [-]           |  [x]            |  [✓]         |
 | `tke`         | [x]           |  [✓]            |  [-]         |

driver/Forcing.jl

@@ -10,19 +10,19 @@ function initialize(::ForcingBase{ForcingLES}, grid, state, LESDat::LESData)
         getvar(var) = TC.pyinterp(vec(grid.zc.z), zc_les, TC.mean_nc_data(data, "profiles", var, imin, imax))
 
         dTdt_hadv = getvar("dtdt_hadv")
-        T_nudge = getvar("temperature_mean")
+        H_nudge = getvar("thetali_mean")
         dTdt_fluc = getvar("dtdt_fluc")
         dqtdt_hadv = getvar("dqtdt_hadv")
         qt_nudge = getvar("qt_mean")
         dqtdt_fluc = getvar("dqtdt_fluc")
         subsidence = getvar("ls_subsidence")
         u_nudge = getvar("u_mean")
         v_nudge = getvar("v_mean")
-        (; dTdt_hadv, T_nudge, dTdt_fluc, dqtdt_hadv, qt_nudge, dqtdt_fluc, subsidence, u_nudge, v_nudge)
+        (; dTdt_hadv, H_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.T_nudge[k] = nt.T_nudge[k]
+        aux_gm.H_nudge[k] = nt.H_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

@@ -66,7 +66,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,
-        ρe_tot_flux = shf + lhf,
+        ρθ_liq_ice_flux = shf / TD.cp_m(thermo_params, ts_in),
         ρq_tot_flux = lhf / TD.latent_heat_vapor(thermo_params, ts_in),
         wstar = convective_vel,
         ρq_liq_flux = FT(0),
@@ -125,7 +125,7 @@ function get_surface(
         ustar = result.ustar,
         ρu_flux = result.ρτxz,
         ρv_flux = result.ρτyz,
-        ρe_tot_flux = shf + lhf,
+        ρθ_liq_ice_flux = shf / TD.cp_m(thermo_params, ts_in),
         ρq_tot_flux = lhf / TD.latent_heat_vapor(thermo_params, 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,
-        ρe_tot_flux = shf + lhf,
+        ρθ_liq_ice_flux = shf / TD.cp_m(thermo_params, ts_in),
         ρq_tot_flux = lhf / TD.latent_heat_vapor(thermo_params, ts_in),
         bflux = result.buoy_flux,
         wstar = convective_vel,

driver/compute_diagnostics.jl

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

driver/dycore.jl

@@ -110,18 +110,15 @@ function compute_ref_state!(
 end
 
 
-function set_thermo_state_peq!(state, grid, moisture_model, param_set)
+function set_thermo_state_from_prog!(state, grid, moisture_model, param_set)
     Ic = CCO.InterpolateF2C()
     thermo_params = TCP.thermodynamics_params(param_set)
     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_uₕ = TC.grid_mean_uₕ(state)
     p_c = aux_gm.p
     ρ_c = prog_gm.ρ
-    C123 = CCG.Covariant123Vector
-    @. aux_gm.e_kin = LA.norm_sqr(C123(prog_gm_uₕ) + C123(Ic(prog_gm_f.w))) / 2
 
     @inbounds for k in TC.real_center_indices(grid)
         thermo_args = if moisture_model isa TC.EquilibriumMoisture
@@ -131,14 +128,18 @@ function set_thermo_state_peq!(state, grid, moisture_model, param_set)
         else
             error("Something went wrong. The moisture_model options are equilibrium or nonequilibrium")
         end
-        e_pot = TC.geopotential(param_set, grid.zc.z[k])
-        e_int = prog_gm.ρe_tot[k] / ρ_c[k] - aux_gm.e_kin[k] - e_pot
-        ts_gm[k] = TC.thermo_state_peq(param_set, p_c[k], e_int, prog_gm.ρq_tot[k] / ρ_c[k], thermo_args...)
+        ts_gm[k] = TC.thermo_state_pθq(
+            param_set,
+            p_c[k],
+            prog_gm.ρθ_liq_ice[k] / ρ_c[k],
+            prog_gm.ρq_tot[k] / ρ_c[k],
+            thermo_args...,
+        )
     end
     return nothing
 end
 
-function set_thermo_state_pθq!(state, grid, moisture_model, param_set)
+function set_thermo_state_from_aux!(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)
@@ -159,23 +160,13 @@ end
 
 function set_grid_mean_from_thermo_state!(param_set, state, grid)
     thermo_params = TCP.thermodynamics_params(param_set)
-    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)
-    prog_gm_uₕ = TC.grid_mean_uₕ(state)
+    ts_gm = TC.center_aux_grid_mean(state).ts
     ρ_c = prog_gm.ρ
-    C123 = CCG.Covariant123Vector
-    @. prog_gm.ρe_tot =
-        ρ_c * TD.total_energy(
-            thermo_params,
-            ts_gm,
-            LA.norm_sqr(C123(prog_gm_uₕ) + C123(Ic(prog_gm_f.w))) / 2,
-            TC.geopotential(param_set, grid.zc.z),
-        )
-    @. prog_gm.ρq_tot = ρ_c * aux_gm.q_tot
-
+    @. prog_gm.ρθ_liq_ice = ρ_c * TD.liquid_ice_pottemp(thermo_params, ts_gm)
+    @. prog_gm.ρq_tot = ρ_c * TD.total_specific_humidity(thermo_params, ts_gm)
     return nothing
 end
 
@@ -188,12 +179,11 @@ function assign_thermo_aux!(state, grid, moisture_model, param_set)
     @inbounds for k in TC.real_center_indices(grid)
         ts = ts_gm[k]
         aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k]
+        aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k]
         aux_gm.q_liq[k] = TD.liquid_specific_humidity(thermo_params, ts)
         aux_gm.q_ice[k] = TD.ice_specific_humidity(thermo_params, ts)
         aux_gm.T[k] = TD.air_temperature(thermo_params, ts)
         aux_gm.RH[k] = TD.relative_humidity(thermo_params, ts)
-        aux_gm.θ_liq_ice[k] = TD.liquid_ice_pottemp(thermo_params, ts)
-        aux_gm.h_tot[k] = TC.total_enthalpy(param_set, prog_gm.ρe_tot[k] / ρ_c[k], ts)
     end
     return
 end
@@ -233,7 +223,7 @@ function ∑tendencies!(tendencies::FV, prog::FV, params::NT, t::Real) where {NT
         state = TC.column_state(prog, aux, tendencies, inds...)
         grid = TC.Grid(state)
 
-        set_thermo_state_peq!(state, grid, edmf.moisture_model, param_set)
+        set_thermo_state_from_prog!(state, grid, edmf.moisture_model, param_set)
         assign_thermo_aux!(state, grid, edmf.moisture_model, param_set)
 
         aux_gm = TC.center_aux_grid_mean(state)
@@ -288,9 +278,6 @@ function compute_gm_tendencies!(
 )
     thermo_params = TCP.thermodynamics_params(param_set)
     Ic = CCO.InterpolateF2C()
-    R_d = TCP.R_d(param_set)
-    T_0 = TCP.T_0(param_set)
-    Lv_0 = TCP.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)
@@ -299,7 +286,7 @@ function compute_gm_tendencies!(
     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)
-    ∇MSE_gm = TC.center_aux_grid_mean(state).∇MSE_gm
+    ∇θ_liq_ice_gm = TC.center_aux_grid_mean(state).∇θ_liq_ice_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)
@@ -311,13 +298,13 @@ function compute_gm_tendencies!(
     aux_tc = TC.center_aux_turbconv(state)
     ts_gm = TC.center_aux_grid_mean(state).ts
 
-    MSE_gm_toa = aux_gm.h_tot[kc_toa] - aux_gm.e_kin[kc_toa]
+    θ_liq_ice_gm_toa = prog_gm.ρθ_liq_ice[kc_toa] / ρ_c[kc_toa]
     q_tot_gm_toa = prog_gm.ρq_tot[kc_toa] / ρ_c[kc_toa]
-    RBe = CCO.RightBiasedC2F(; top = CCO.SetValue(MSE_gm_toa))
+    RBθ = CCO.RightBiasedC2F(; top = CCO.SetValue(θ_liq_ice_gm_toa))
     RBq = CCO.RightBiasedC2F(; top = CCO.SetValue(q_tot_gm_toa))
     wvec = CC.Geometry.WVector
     ∇c = CCO.DivergenceF2C()
-    @. ∇MSE_gm = ∇c(wvec(RBe(aux_gm.h_tot - aux_gm.e_kin)))
+    @. ∇θ_liq_ice_gm = ∇c(wvec(RBθ(prog_gm.ρθ_liq_ice / ρ_c)))
     @. ∇q_tot_gm = ∇c(wvec(RBq(prog_gm.ρq_tot / ρ_c)))
 
     if edmf.moisture_model isa TC.NonEquilibriumMoisture
@@ -334,9 +321,7 @@ function compute_gm_tendencies!(
     # Apply forcing and radiation
     prog_gm_uₕ = TC.grid_mean_uₕ(state)
     aux_gm_uₕ_g = TC.grid_mean_uₕ_g(state)
-    # prog_gm_v = TC.grid_mean_v(state)
     tendencies_gm_uₕ = TC.tendencies_grid_mean_uₕ(state)
-    # tendencies_gm_v = TC.tendencies_grid_mean_v(state)
     prog_gm_u = TC.physical_grid_mean_u(state)
     prog_gm_v = TC.physical_grid_mean_v(state)
 
@@ -355,28 +340,22 @@ function compute_gm_tendencies!(
 
 
     @inbounds for k in TC.real_center_indices(grid)
-        cp_m = TD.cp_m(thermo_params, ts_gm[k])
-        cp_v = TCP.cp_v(param_set)
-        cv_v = TCP.cv_v(param_set)
-        R_v = TCP.R_v(param_set)
-        cv_m = TD.cv_m(thermo_params, ts_gm[k])
-        h_v = cv_v * (aux_gm.T[k] - T_0) + Lv_0 - R_v * T_0
-
+        Π = TD.exner(thermo_params, ts_gm[k])
         # LS Subsidence
-        tendencies_gm.ρe_tot[k] -= ρ_c[k] * aux_gm.subsidence[k] * ∇MSE_gm[k]
+        tendencies_gm.ρθ_liq_ice[k] -= ρ_c[k] * aux_gm.subsidence[k] * ∇θ_liq_ice_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 Cases.rad_type(radiation) <: Union{Cases.RadiationDYCOMS_RF01, Cases.RadiationLES, Cases.RadiationTRMM_LBA}
-            tendencies_gm.ρe_tot[k] += ρ_c[k] * cv_m * aux_gm.dTdt_rad[k]
+            tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * aux_gm.dTdt_rad[k] / Π
         end
         # LS advection
         tendencies_gm.ρq_tot[k] += ρ_c[k] * aux_gm.dqtdt_hadv[k]
         if !(Cases.force_type(force) <: Cases.ForcingDYCOMS_RF01)
-            tendencies_gm.ρe_tot[k] += ρ_c[k] * (cp_m * aux_gm.dTdt_hadv[k] + h_v * aux_gm.dqtdt_hadv[k])
+            tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * aux_gm.dTdt_hadv[k] / Π
         end
         if edmf.moisture_model isa TC.NonEquilibriumMoisture
             tendencies_gm.q_liq[k] += aux_gm.dqldt[k]
@@ -385,22 +364,21 @@ function compute_gm_tendencies!(
 
         # LES specific forcings
         if Cases.force_type(force) <: Cases.ForcingLES
-            T_fluc = aux_gm.dTdt_fluc[k]
+            H_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_τᵣ
 
             Γᵣ = compute_les_Γᵣ(grid.zc[k].z, force.scalar_nudge_τᵣ, force.scalar_nudge_zᵢ, force.scalar_nudge_zᵣ)
-            gm_T_nudge_k = Γᵣ * (aux_gm.T_nudge[k] - aux_gm.T[k])
+            gm_H_nudge_k = Γᵣ * (aux_gm.H_nudge[k] - aux_gm.θ_liq_ice[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.ρ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.ρθ_liq_ice[k] += ρ_c[k] * (gm_H_nudge_k + H_fluc)
             tendencies_gm_uₕ[k] += CCG.Covariant12Vector(CCG.UVVector(gm_U_nudge_k, gm_V_nudge_k), lg[k])
             if edmf.moisture_model isa TC.NonEquilibriumMoisture
                 tendencies_gm.q_liq[k] += aux_gm.dqldt_hadv[k] + gm_q_liq_nudge_k + aux_gm.dqldt_fluc[k]
@@ -416,28 +394,28 @@ function compute_gm_tendencies!(
                 aux_tc.qt_tendency_precip_sinks[k]
             )
 
-        tendencies_gm.ρe_tot[k] +=
+        tendencies_gm.ρθ_liq_ice[k] +=
             ρ_c[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]
+                aux_bulk.θ_liq_ice_tendency_precip_formation[k] +
+                aux_en.θ_liq_ice_tendency_precip_formation[k] +
+                aux_tc.θ_liq_ice_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, param_set)
-    sgs_flux_h_tot = aux_gm_f.sgs_flux_h_tot
+    TC.compute_sgs_flux!(edmf, grid, state, surf)
+    sgs_flux_θ_liq_ice = aux_gm_f.sgs_flux_θ_liq_ice
     sgs_flux_q_tot = aux_gm_f.sgs_flux_q_tot
     sgs_flux_uₕ = aux_gm_f.sgs_flux_uₕ
-    tends_ρe_tot = tendencies_gm.ρe_tot
+    tends_ρθ_liq_ice = tendencies_gm.ρθ_liq_ice
     tends_ρq_tot = tendencies_gm.ρq_tot
     tends_uₕ = TC.tendencies_grid_mean_uₕ(state)
 
     ∇sgs = CCO.DivergenceF2C()
-    @. tends_ρe_tot += -∇sgs(wvec(sgs_flux_h_tot))
+    @. tends_ρθ_liq_ice += -∇sgs(wvec(sgs_flux_θ_liq_ice))
     @. tends_ρq_tot += -∇sgs(wvec(sgs_flux_q_tot))
     @. tends_uₕ += -∇sgs(sgs_flux_uₕ) / ρ_c
 

driver/dycore_variables.jl

@@ -48,23 +48,21 @@ cent_aux_vars_gm(FT, local_geometry, edmf) = (;
     subsidence = FT(0), #Large-scale subsidence
     dTdt_hadv = FT(0), #Horizontal advection of temperature
     dqtdt_hadv = FT(0), #Horizontal advection of moisture
-    T_nudge = FT(0), #Reference T profile for relaxation tendency
+    H_nudge = FT(0), #Reference thetali 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
     uₕ_g = CCG.Covariant12Vector(CCG.UVVector(FT(0), FT(0)), local_geometry), #Geostrophic u velocity
-    ∇MSE_gm = FT(0),
+    ∇θ_liq_ice_gm = FT(0),
     ∇q_tot_gm = FT(0),
     cent_aux_vars_gm_moisture(FT, edmf.moisture_model)...,
     θ_virt = FT(0),
     Ri = FT(0),
     θ_liq_ice = FT(0),
     q_tot = FT(0),
     p = FT(0),
-    e_kin = FT(0),
-    h_tot = FT(0),
 )
 cent_aux_vars(FT, local_geometry, edmf) =
     (; cent_aux_vars_gm(FT, local_geometry, edmf)..., TC.cent_aux_vars_edmf(FT, local_geometry, edmf)...)
@@ -77,7 +75,7 @@ face_aux_vars_gm(FT, local_geometry, edmf) = (;
     diffusive_flux_s = FT(0),
     total_flux_s = FT(0),
     f_rad = FT(0),
-    sgs_flux_h_tot = FT(0),
+    sgs_flux_θ_liq_ice = FT(0),
     sgs_flux_q_tot = FT(0),
     face_aux_vars_gm_moisture(FT, edmf.moisture_model)...,
     sgs_flux_uₕ = CCG.Covariant3Vector(FT(0)) ⊗ CCG.Covariant12Vector(FT(0), FT(0)),
@@ -128,7 +126,7 @@ cent_prognostic_vars_gm_moisture(::Type{FT}, ::TC.EquilibriumMoisture) where {FT
 cent_prognostic_vars_gm(::Type{FT}, local_geometry, edmf) where {FT} = (;
     ρ = FT(0),
     uₕ = CCG.Covariant12Vector(CCG.UVVector(FT(0), FT(0)), local_geometry),
-    ρe_tot = FT(0),
+    ρθ_liq_ice = FT(0),
     ρq_tot = FT(0),
     cent_prognostic_vars_gm_moisture(FT, edmf.moisture_model)...,
 )

driver/main.jl

@@ -266,7 +266,7 @@ function initialize(sim::Simulation1d)
             end
         end
         Cases.initialize_profiles(case, grid, param_set, state; les_data_kwarg...)
-        set_thermo_state_pθq!(state, grid, edmf.moisture_model, param_set)
+        set_thermo_state_from_aux!(state, grid, edmf.moisture_model, param_set)
         set_grid_mean_from_thermo_state!(param_set, state, grid)
         assign_thermo_aux!(state, grid, edmf.moisture_model, param_set)
         Cases.initialize_forcing(case, forcing, grid, state, param_set; les_data_kwarg...)

integration_tests/sphere_utils.jl

@@ -79,7 +79,7 @@ function test_zero_horizontal_variance(Y; filter_prop_chain = pc -> true)
         (:cent, :ρ) => false,
         (:cent, :uₕ, :components, :data, 1) => false,
         (:cent, :uₕ, :components, :data, 2) => false,
-        (:cent, :ρe_tot) => false,
+        (:cent, :ρθ_liq_ice) => false,
         (:cent, :ρq_tot) => false,
         (:cent, :turbconv, :en, :ρatke) => false,
         (:cent, :turbconv, :up, 1, :ρarea) => false,