Merge #174

174: Add bcs to interpc2f r=charleskawczynski a=charleskawczynski



Co-authored-by: Charles Kawczynski <kawczynski.charles@gmail.com>

src/Operators.jl

@@ -56,13 +56,26 @@ end
 # top of the domain.
 ∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value
 
-function interpc2f(f, grid::Grid, k::Int)
-    return 0.5 * (f[k + 1] + f[k])
-end
+# Used when traversing cell faces
+
+interpc2f(f, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError()) =
+    interpc2f(dual_centers(f, grid, k), grid, k; bottom, top)
+
+interpc2f(f, grid::Grid, k::Int, i_up::Int; bottom = NoBCGivenError(), top = NoBCGivenError()) =
+    interpc2f(dual_centers(f, grid, k, i_up), grid, k; bottom, top)
 
-function interpc2f(f, grid::Grid, k::Int, i_up::Int)
-    return 0.5 * (f[i_up, k + 1] + f[i_up, k])
+function interpc2f(f_dual::SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError())
+    if is_surface_face(grid, k)
+        return interpc2f(f_dual, grid, BottomBCTag(), bottom)
+    elseif is_toa_face(grid, k)
+        return interpc2f(f_dual, grid, TopBCTag(), top)
+    else
+        return interpc2f(f_dual, grid, InteriorTag())
+    end
 end
+interpc2f(f::SVector, grid::Grid, ::InteriorTag) = (f[1] + f[2]) / 2
+interpc2f(f::SVector, grid::Grid, ::TopBCTag, bc::SetValue) = bc.value
+interpc2f(f::SVector, grid::Grid, ::BottomBCTag, bc::SetValue) = bc.value
 
 function interpf2c(f, grid::Grid, k::Int)
     return 0.5 * (f[k] + f[k - 1])
@@ -261,3 +274,22 @@ end
 # A 2-point field stencil for ordinary and updraft variables
 dual_faces(f::AbstractVector, grid, k::Int) = SVector(f[k - 1], f[k])
 dual_faces(f::AbstractMatrix, grid, k::Int, i_up::Int) = SVector(f[i_up, k - 1], f[i_up, k])
+
+function dual_centers(f::AbstractVector, grid, k::Int)
+    if is_surface_face(grid, k)
+        return SVector(f[k + 1])
+    elseif is_toa_face(grid, k)
+        return SVector(f[k])
+    else
+        return SVector(f[k], f[k + 1])
+    end
+end
+function dual_centers(f::AbstractMatrix, grid, k::Int, i_up::Int)
+    if is_surface_face(grid, k)
+        return SVector(f[i_up, k + 1])
+    elseif is_toa_face(grid, k)
+        return SVector(f[i_up, k])
+    else
+        return SVector(f[i_up, k], f[i_up, k + 1])
+    end
+end

src/Turbulence_PrognosticTKE.jl

@@ -710,13 +710,18 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl
         end
     end
 
+    kf_surf = kf_surface(grid)
+    kf_toa = kf_top_of_atmos(grid)
     @inbounds for k in real_face_indicies(grid)
         @inbounds for i in xrange(self.n_updrafts)
 
             # pressure
-            a_kfull = interpc2f(self.UpdVar.Area.values, grid, k, i)
+            a_bcs = (; bottom = SetValue(self.area_surface_bc[i]), top = SetValue(0))
+            a_kfull = interpc2f(self.UpdVar.Area.values, grid, k, i; a_bcs...)
             if a_kfull > 0.0
-                b_kfull = interpc2f(self.UpdVar.B.values, grid, k, i)
+                B = self.UpdVar.B.values
+                b_bcs = (; bottom = SetValue(B[i, kf_surf]), top = SetValue(B[i, kf_toa]))
+                b_kfull = interpc2f(self.UpdVar.B.values, grid, k, i; b_bcs...)
                 w_cut = fcut(self.UpdVar.W.values, grid, k, i)
                 ∇w_up = f∇(w_cut, grid, k; bottom = SetValue(0), top = SetGradient(0))
                 asp_ratio = 1.0