From c1115d9d42eac4581b7a1721d21882f8b036e9a8 Mon Sep 17 00:00:00 2001
From: Simon Byrne <simonbyrne@gmail.com>
Date: Fri, 15 Sep 2023 10:28:28 -0700
Subject: [PATCH] Use weak divergence in scalar advection

---
 docs/src/equations.md                 |  7 +++----
 src/prognostic_equations/advection.jl | 18 +++++++++---------
 2 files changed, 12 insertions(+), 13 deletions(-)

diff --git a/docs/src/equations.md b/docs/src/equations.md
index 181eb2b8512..62f073a34fd 100644
--- a/docs/src/equations.md
+++ b/docs/src/equations.md
@@ -33,7 +33,6 @@ We make use of the following operators
 
 ### Differential operators
 
-- ``\mathcal{D}_h`` is the [discrete horizontal spectral divergence](https://clima.github.io/ClimaCore.jl/stable/operators/#ClimaCore.Operators.Divergence).
 - ``\hat{\mathcal{D}}_h`` is the [discrete horizontal spectral weak divergence](https://clima.github.io/ClimaCore.jl/stable/operators/#ClimaCore.Operators.WeakDivergence).
 - ``\mathcal{D}^c_v`` is the [face-to-center vertical divergence](https://clima.github.io/ClimaCore.jl/stable/operators/#ClimaCore.Operators.DivergenceF2C).
 !!! todo
@@ -107,7 +106,7 @@ Follows the continuity equation
 This is discretized using the following
 ```math
 \frac{\partial}{\partial t} \rho
-= - \mathcal{D}_h[ \rho \bar{\boldsymbol{u}}] - \mathcal{D}^c_v \left[WI^f( J, \rho) \tilde{\boldsymbol{u}} \right]
+= - \hat{\mathcal{D}}_h[ \rho \bar{\boldsymbol{u}}] - \mathcal{D}^c_v \left[WI^f( J, \rho) \tilde{\boldsymbol{u}} \right]
 ```
 
 with the
@@ -210,7 +209,7 @@ which is stabilized with the addition of a 4th-order hyperdiffusion term on tota
 This is discretized using
 ```math
 \frac{\partial}{\partial t} \rho e \approx
-- \mathcal{D}_h[ (\rho e + p) \bar{\boldsymbol{u}} ]
+- \hat{\mathcal{D}}_h[ (\rho e + p) \bar{\boldsymbol{u}} ]
 - \mathcal{D}^c_v \left[ WI^f(J,\rho) \,  \tilde{\boldsymbol{u}} \, I^f \left(\frac{\rho e + p}{\rho} \right)
   + \boldsymbol{F}_R \right] - \nu_h \hat{\mathcal{D}}_h( \rho \mathcal{G}_h(\psi) ).
 ```
@@ -247,7 +246,7 @@ This is stabilized with the addition of a 4th-order hyperdiffusion term
 This is discretized using
 ```math
 \frac{\partial}{\partial t} \rho \chi \approx
-- \mathcal{D}_h[ \rho \chi \bar{\boldsymbol{u}}]
+- \hat{\mathcal{D}}_h[ \rho \chi \bar{\boldsymbol{u}}]
 - \mathcal{D}^c_v \left[ WI^f(J,\rho) \, U^f\left( \tilde{\boldsymbol{u}},  \frac{\rho \chi}{\rho} \right) \right]
 - \nu_\chi \hat{\mathcal{D}}_h ( \rho \, \mathcal{G}_h (\psi) )
 ```
diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl
index aff688bd18e..9e73101eaa7 100644
--- a/src/prognostic_equations/advection.jl
+++ b/src/prognostic_equations/advection.jl
@@ -16,34 +16,34 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t)
         (; ᶜuʲs) = p
     end
 
-    @. Yₜ.c.ρ -= divₕ(Y.c.ρ * ᶜu)
+    @. Yₜ.c.ρ -= wdivₕ(Y.c.ρ * ᶜu)
     if p.atmos.turbconv_model isa EDMFX
         for j in 1:n
-            @. Yₜ.c.sgsʲs.:($$j).ρa -= divₕ(Y.c.sgsʲs.:($$j).ρa * ᶜuʲs.:($$j))
+            @. Yₜ.c.sgsʲs.:($$j).ρa -= wdivₕ(Y.c.sgsʲs.:($$j).ρa * ᶜuʲs.:($$j))
         end
     end
 
     if :ρθ in propertynames(Y.c)
-        @. Yₜ.c.ρθ -= divₕ(Y.c.ρθ * ᶜu)
+        @. Yₜ.c.ρθ -= wdivₕ(Y.c.ρθ * ᶜu)
     elseif :ρe_tot in propertynames(Y.c)
         (; ᶜh_tot) = p
-        @. Yₜ.c.ρe_tot -= divₕ(Y.c.ρ * ᶜh_tot * ᶜu)
+        @. Yₜ.c.ρe_tot -= wdivₕ(Y.c.ρ * ᶜh_tot * ᶜu)
     end
     if p.atmos.turbconv_model isa EDMFX
         for j in 1:n
             if :ρθ in propertynames(Y.c)
                 @. Yₜ.c.sgsʲs.:($$j).ρaθ -=
-                    divₕ(Y.c.sgsʲs.:($$j).ρaθ * ᶜuʲs.:($$j))
+                    wdivₕ(Y.c.sgsʲs.:($$j).ρaθ * ᶜuʲs.:($$j))
             elseif :ρe_tot in propertynames(Y.c)
                 (; ᶜh_totʲs) = p
                 @. Yₜ.c.sgsʲs.:($$j).ρae_tot -=
-                    divₕ(Y.c.sgsʲs.:($$j).ρa * ᶜh_totʲs.:($$j) * ᶜuʲs.:($$j))
+                    wdivₕ(Y.c.sgsʲs.:($$j).ρa * ᶜh_totʲs.:($$j) * ᶜuʲs.:($$j))
             end
         end
     end
 
     if use_prognostic_tke(p.atmos.turbconv_model)
-        @. Yₜ.c.sgs⁰.ρatke -= divₕ(Y.c.sgs⁰.ρatke * ᶜu⁰)
+        @. Yₜ.c.sgs⁰.ρatke -= wdivₕ(Y.c.sgs⁰.ρatke * ᶜu⁰)
     end
 
     @. Yₜ.c.uₕ -= C12(gradₕ(ᶜp - ᶜp_ref) / Y.c.ρ + gradₕ(ᶜK + ᶜΦ))
@@ -59,7 +59,7 @@ NVTX.@annotate function horizontal_tracer_advection_tendency!(Yₜ, Y, p, t)
     end
 
     for ρχ_name in filter(is_tracer_var, propertynames(Y.c))
-        @. Yₜ.c.:($$ρχ_name) -= divₕ(Y.c.:($$ρχ_name) * ᶜu)
+        @. Yₜ.c.:($$ρχ_name) -= wdivₕ(Y.c.:($$ρχ_name) * ᶜu)
     end
 
     if p.atmos.turbconv_model isa EDMFX
@@ -67,7 +67,7 @@ NVTX.@annotate function horizontal_tracer_advection_tendency!(Yₜ, Y, p, t)
             for ρaχ_name in
                 filter(is_tracer_var, propertynames(Y.c.sgsʲs.:($j)))
                 @. Yₜ.c.sgsʲs.:($$j).:($$ρaχ_name) -=
-                    divₕ(Y.c.sgsʲs.:($$j).:($$ρaχ_name) * ᶜuʲs.:($$j))
+                    wdivₕ(Y.c.sgsʲs.:($$j).:($$ρaχ_name) * ᶜuʲs.:($$j))
             end
         end
     end