Abstract away operator

CliMA · Oct 18, 2021 · 3e07852 · 3e07852
1 parent b7772ce
commit 3e07852
Show file tree

Hide file tree

Showing 2 changed files with 71 additions and 13 deletions.
diff --git a/src/Operators.jl b/src/Operators.jl
@@ -250,6 +250,49 @@ function c∇(f::SA.SVector, grid::Grid, ::BottomBCTag, ::Extrapolate)
     return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
 end
 
+# ∇(center data) for possibly vanishing subdomains.
+
+"""
+    c∇_vanishing_subdomain
+
+Used when the vertical gradient of a field must be computed over a conditional subdomain which may vanish
+below or above the current cell. If the subdomains vanishes, a default user-defined gradient is returned. 
+"""
+function c∇_vanishing_subdomain(
+    f_cut::SA.SVector,
+    sd_cut::SA.Svector,
+    default∇::FT,
+    grid::Grid,
+    k;
+    bottom = NoBCGivenError(),
+    top = NoBCGivenError(),
+) where {FT <: Real}
+    if is_surface_center(grid, k)
+        return c∇(f_cut, grid, BottomBCTag(), bottom)
+    elseif is_toa_center(grid, k)
+        return c∇(f_cut, grid, TopBCTag(), top)
+    else
+        return c∇_vanishing_subdomain(f_cut, sd_cut, default∇, grid, InteriorTag())
+    end
+end
+function c∇_vanishing_subdomain(
+    f::SA.SVector,
+    sd::SA.Svector,
+    default∇::FT,
+    grid::Grid,
+    ::InteriorTag,
+) where {FT <: Real}
+    @assert length(f) == 3
+    @assert length(sd) == 3
+    if sd[1] * sd[3] ≈ FT(0)
+        return default∇
+    else
+        f_dual⁺ = SA.SVector(f[2], f[3])
+        f_dual⁻ = SA.SVector(f[1], f[2])
+        return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
+    end
+end
+
 #####
 ##### ∇(face data)
 #####

diff --git a/src/update_aux.jl b/src/update_aux.jl
@@ -334,22 +334,37 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS)
 
         elseif edmf.bg_closure == BuoyGradQuadratures()
             # Second order approximation: Use dry and cloudy environmental fields.
-            cf = aux_en.cloud_fraction
+            cf_cut = ccut(aux_en.cloud_fraction, grid, k)
             QT_cloudy_cut = ccut(en_thermo.qt_cloudy, grid, k)
-            ∂qt∂z_cloudy = c∇(QT_cloudy_cut, grid, k; bottom = SetGradient(0), top = SetGradient(0))
+            ∂qt∂z_cloudy = c∇_vanishing_subdomain(
+                QT_cloudy_cut,
+                cf_cut,
+                ∂qt∂z,
+                grid,
+                k;
+                bottom = SetGradient(0),
+                top = SetGradient(0),
+            )
             THL_cloudy_cut = ccut(en_thermo.thl_cloudy, grid, k)
-            ∂θl∂z_cloudy = c∇(THL_cloudy_cut, grid, k; bottom = SetGradient(0), top = SetGradient(0))
+            ∂θl∂z_cloudy = c∇_vanishing_subdomain(
+                THL_cloudy_cut,
+                cf_cut,
+                ∂θl∂z,
+                grid,
+                k;
+                bottom = SetGradient(0),
+                top = SetGradient(0),
+            )
             THV_dry_cut = ccut(en_thermo.thv_dry, grid, k)
-            ∂θv∂z_dry = c∇(THV_dry_cut, grid, k; bottom = SetGradient(0), top = SetGradient(0))
-            # Gradient safeguards for dry/cloudy -- default to environmental mean
-            # Assumes no fog or domain-top cloudiness
-            if k != kc_surface(grid) && k != kc_top_of_atmos(grid)
-                if cf[k - 1] * cf[k + 1] ≈ 0.0
-                    ∂qt∂z_cloudy = ∂qt∂z
-                    ∂θl∂z_cloudy = ∂θl∂z
-                end
-                ∂θv∂z_dry = (1 - cf[k - 1]) * (1 - cf[k + 1]) ≈ 0.0 ? ∂θv∂z : ∂θv∂z_dry
-            end
+            ∂θv∂z_dry = c∇_vanishing_subdomain(
+                THV_dry_cut,
+                cf_cut,
+                ∂θv∂z,
+                grid,
+                k;
+                bottom = SetGradient(0),
+                top = SetGradient(0),
+            )
 
             bg_kwargs = (;
                 t_cloudy = en_thermo.t_cloudy[k],