Use more ClimaCore operators

CliMA · Dec 2, 2021 · d99c9d9 · d99c9d9
1 parent 970336f
commit d99c9d9
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 4 deletions.
diff --git a/integration_tests/utils/main.jl b/integration_tests/utils/main.jl
@@ -176,6 +176,8 @@ cent_aux_vars_edmf(FT, n_up) = (;
         KH = FT(0),
         mixing_length = FT(0),
         θ_virt = FT(0),
+        ∇θ_virt = FT(0),
+        ∇Ri_bulk = FT(0),
         massflux_tendency_h = FT(0),
         massflux_tendency_qt = FT(0),
         diffusive_tendency_h = FT(0),

diff --git a/src/turbulence_functions.jl b/src/turbulence_functions.jl
@@ -16,15 +16,20 @@ function get_inversion(grid::Grid, state, param_set, Ri_bulk_crit)
     u = center_prog_grid_mean(state).u
     v = center_prog_grid_mean(state).v
     Ri_bulk = center_aux_bulk(state).Ri
+    ∇θ_virt = center_aux_turbconv(state).∇θ_virt
+    ∇Ri_bulk = center_aux_turbconv(state).∇Ri_bulk
     θ_virt_b = θ_virt[kc_surf]
     z_c = grid.zc
+    ∇c = CCO.DivergenceF2C()
+    wvec = CC.Geometry.WVector
 
     # test if we need to look at the free convective limit
     if (u[kc_surf]^2 + v[kc_surf]^2) <= 0.01
         kmask = map(k -> (k, θ_virt[k] > θ_virt_b), real_center_indices(grid))
         k_star = first(kmask[findlast(km -> km[2], kmask)])
-        ∇θ_virt = c∇_upwind(θ_virt, grid, k_star; bottom = SetGradient(0), top = FreeBoundary())
-        h = (θ_virt_b - θ_virt[k_star - 1]) / ∇θ_virt + z_c[k_star - 1]
+        LB = CCO.LeftBiasedF2C(; bottom = CCO.SetValue(θ_virt[kc_surf]))
+        @. ∇θ_virt = ∇c(wvec(LB(θ_virt)))
+        h = (θ_virt_b - θ_virt[k_star - 1]) / ∇θ_virt[k_star] + z_c[k_star - 1]
     else
         Ri_bulk_fn(k) = g * (θ_virt[k] - θ_virt_b) * z_c[k] / θ_virt_b / (u[k] * u[k] + v[k] * v[k])
 
@@ -33,8 +38,9 @@ function get_inversion(grid::Grid, state, param_set, Ri_bulk_crit)
         end
         kmask = map(k -> (k, Ri_bulk_fn(k) > Ri_bulk_crit), real_center_indices(grid))
         k_star = first(kmask[findlast(km -> km[2], kmask)])
-        ∇Ri_bulk = c∇_upwind(Ri_bulk, grid, k_star; bottom = SetGradient(0), top = FreeBoundary())
-        h = (Ri_bulk_crit - Ri_bulk[k_star - 1]) / ∇Ri_bulk + z_c[k_star - 1]
+        LB = CCO.LeftBiasedF2C(; bottom = CCO.SetValue(Ri_bulk[kc_surf]))
+        @. ∇Ri_bulk = ∇c(wvec(LB(Ri_bulk)))
+        h = (Ri_bulk_crit - Ri_bulk[k_star - 1]) / ∇Ri_bulk[k_star] + z_c[k_star - 1]
     end
 
     return h