change large scale tendencies to energy

CliMA · May 3, 2022 · c159d17 · c159d17
1 parent 1e4bcb1
commit c159d17
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 12 deletions.
diff --git a/driver/dycore.jl b/driver/dycore.jl
@@ -405,13 +405,13 @@ 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))
-    RBq = CCO.RightBiasedC2F(; top = CCO.SetValue(q_tot_gm_toa))
+    ρh_tot_gm_toa = aux_gm.subsidence[kc_toa] * (prog_gm.ρe_tot[kc_toa] + p0_c[kc_toa])
+    ρq_tot_gm_toa = aux_gm.subsidence[kc_toa] * prog_gm.ρq_tot[kc_toa]
+    RBe = CCO.RightBiasedC2F(; top = CCO.SetValue(ρh_tot_gm_toa))
+    RBq = CCO.RightBiasedC2F(; top = CCO.SetValue(ρq_tot_gm_toa))
     wvec = CC.Geometry.WVector
     ∇c = CCO.DivergenceF2C()
-    @. ∇h_tot_subsidence_flux = ∇c(wvec(RBθ(aux_gm.subsidence * (prog_gm.ρe_tot + p0_c))))
+    @. ∇h_tot_subsidence_flux = ∇c(wvec(RBe(aux_gm.subsidence * (prog_gm.ρe_tot + p0_c))))
     @. ∇q_tot_subsidence_flux = ∇c(wvec(RBq(aux_gm.subsidence * prog_gm.ρq_tot)))
 
     if edmf.moisture_model isa TC.NonEquilibriumMoisture
@@ -428,6 +428,8 @@ function compute_gm_tendencies!(
     @inbounds for k in TC.real_center_indices(grid)
         # Apply large-scale horizontal advection tendencies
         c_pm = TD.cp_m(param_set, ts_gm[k])
+        Lv = TD.latent_heat_vapor(param_set, ts_gm[k])
+        Π = TD.exner(param_set, ts_gm[k])
 
         if force.apply_coriolis
             tendencies_gm.u[k] -= force.coriolis_param * (aux_gm.vg[k] - prog_gm.v[k])
@@ -438,6 +440,7 @@ function compute_gm_tendencies!(
         end
         if TC.force_type(force) <: TC.ForcingDYCOMS_RF01
             tendencies_gm.ρq_tot[k] += ρ0_c[k] * aux_gm.dqtdt[k]
+            tendencies_gm.ρe_tot[k] += ρ0_c[k] * Lv * aux_gm.dqtdt[k]
             # Apply large-scale subsidence tendencies
             tendencies_gm.ρe_tot[k] -= ∇h_tot_subsidence_flux[k]
             tendencies_gm.ρq_tot[k] -= ∇q_tot_subsidence_flux[k]
@@ -454,8 +457,8 @@ function compute_gm_tendencies!(
                 tendencies_gm.ρe_tot[k] -= ∇h_tot_subsidence_flux[k]
                 tendencies_gm.ρq_tot[k] -= ∇q_tot_subsidence_flux[k]
             end
-            tendencies_gm.ρe_tot[k] += ρ0_c[k] * c_pm * aux_gm.dTdt[k]
             tendencies_gm.ρq_tot[k] += ρ0_c[k] * aux_gm.dqtdt[k]
+            tendencies_gm.ρe_tot[k] += ρ0_c[k] * (c_pm * aux_gm.dTdt[k] + Lv * aux_gm.dqtdt[k])
             if edmf.moisture_model isa TC.NonEquilibriumMoisture
                 if force.apply_subsidence
                     tendencies_gm.q_liq[k] -= ∇q_liq_gm[k] * aux_gm.subsidence[k]
@@ -467,16 +470,16 @@ function compute_gm_tendencies!(
         end
 
         if TC.force_type(force) <: TC.ForcingLES
-            H_horz_adv = aux_gm.dTdt_hadv[k]
-            H_fluc = aux_gm.dTdt_fluc[k]
+            H_horz_adv = ρ0_c[k] * c_pm * aux_gm.dTdt_hadv[k]
+            H_fluc = ρ0_c[k] * c_pm * aux_gm.dTdt_fluc[k]
 
             gm_U_nudge_k = (aux_gm.u_nudge[k] - prog_gm.u[k]) / force.nudge_tau
             gm_V_nudge_k = (aux_gm.v_nudge[k] - prog_gm.v[k]) / force.nudge_tau
 
             Γᵣ = TC.compute_les_Γᵣ(grid.zc[k])
             if Γᵣ != 0
                 tau_k = 1 / Γᵣ
-                gm_H_nudge_k = (aux_gm.H_nudge[k] - prog_gm.ρe_tot[k] / ρ0_c[k]) / tau_k
+                gm_H_nudge_k = ρ0_c[k] * c_pm * Π * (aux_gm.H_nudge[k] - aux_gm.θ_liq_ice[k]) / tau_k
                 gm_q_tot_nudge_k = (aux_gm.qt_nudge[k] - prog_gm.ρq_tot[k] / ρ0_c[k]) / tau_k
             else
                 gm_H_nudge_k = 0.0
@@ -512,9 +515,13 @@ function compute_gm_tendencies!(
                 end
             end
 
-            tendencies_gm.ρe_tot[k] += ρ0_c[k] * (H_horz_adv + gm_H_nudge_k + H_fluc + gm_H_subsidence_k)
             tendencies_gm.ρq_tot[k] +=
                 ρ0_c[k] * (aux_gm.dqtdt_hadv[k] + gm_q_tot_nudge_k + aux_gm.dqtdt_fluc[k] + gm_QT_subsidence_k)
+            tendencies_gm.ρe_tot[k] +=
+                ρ0_c[k] * (
+                    cp_m * Π * (H_horz_adv + gm_H_nudge_k + H_fluc + gm_H_subsidence_k) +
+                    Lv * (aux_gm.dqtdt_hadv[k] + gm_q_tot_nudge_k + aux_gm.dqtdt_fluc[k] + gm_QT_subsidence_k)
+                )
             tendencies_gm.u[k] += gm_U_nudge_k
             tendencies_gm.v[k] += gm_V_nudge_k
             if edmf.moisture_model isa TC.NonEquilibriumMoisture

diff --git a/src/EDMF_functions.jl b/src/EDMF_functions.jl
@@ -106,7 +106,7 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
         q_tot_up = aux_up_i.q_tot
         @. aux_up_f[i].massflux = ρ0_f * Ifau(a_up) * (w_up - w_gm)
         # We know that, since W = 0 at z = 0, m = 0 also, and
-        # therefore θ_liq_ice / q_tot values do not matter
+        # therefore h_tot / q_tot values do not matter
         ts_up_i = copy(ts_en)
         h_tot_up_i = copy(q_tot_up)
         @. ts_up_i = thermo_state_pθq(param_set, p0_c, θ_liq_ice_up, q_tot_up)
@@ -160,7 +160,7 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
     sgs_flux_v = aux_gm_f.sgs_flux_v
 
     @. sgs_flux_h_tot = diffusive_flux_h + massflux_h
-    @. sgs_flux_q_tot = diffusive_flux_qt# + massflux_qt
+    @. sgs_flux_q_tot = diffusive_flux_qt + massflux_qt
     @. sgs_flux_u = diffusive_flux_u # + massflux_u
     @. sgs_flux_v = diffusive_flux_v # + massflux_v