replace copy_stack_field_down! with boradcasts of reshaped arrays

src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl

@@ -11,7 +11,7 @@ using ...MPIStateArrays
 using ...Mesh.Filters: apply!
 using ...Mesh.Grids: VerticalDirection
 using ...Mesh.Geometry
-using ...DGMethods: DGModel, copy_stack_field_down!
+using ...DGMethods
 using ...DGMethods.NumericalFluxes
 using ...BalanceLaws
 using ..Ocean
@@ -678,13 +678,14 @@ function update_auxiliary_state_gradient!(
     indefinite_stack_integral!(dg, m, Q, A, t, elems) # bottom -> top
     reverse_indefinite_stack_integral!(dg, m, Q, A, t, elems) # top -> bottom
 
+    Nq, Nqk, _, _, nelemv, _, nelemh, _ = basic_grid_info(dg)
+
     # project w(z=0) down the stack
-    # [1] to convert from range to integer
-    # copy_stack_field_down! doesn't like ranges
-    # eventually replace this with a reshape and broadcast
-    index_w = varsindex(vars_state_auxiliary(m, FT), :w)[1]
-    index_wz0 = varsindex(vars_state_auxiliary(m, FT), :wz0)[1]
-    copy_stack_field_down!(dg, m, A, index_w, index_wz0, elems)
+    boxy_w = reshape(A.w, Nq^2, Nqk, 1, nelemv, nelemh)
+    flat_w = @view boxy_w[:, end, :, end, :]
+    flat_wz0 = reshape(flat_w, Nq^2, 1, 1, 1, nelemh)
+    boxy_wz0 = reshape(A.wz0, Nq^2, Nqk, 1, nelemv, nelemh)
+    boxy_wz0 .= flat_wz0
 
     return true
 end