From 32ec1fd6ec151558a70ffac056b59ca2ee198ee5 Mon Sep 17 00:00:00 2001
From: Charles Kawczynski <kawczynski.charles@gmail.com>
Date: Mon, 22 Jul 2024 13:14:51 -0400
Subject: [PATCH] Add Nh to type parameter space

---
 examples/hybrid/driver.jl                |  16 +-
 ext/cuda/data_layouts_copyto.jl          |  39 ++-
 ext/cuda/data_layouts_fill.jl            |   4 +-
 ext/cuda/data_layouts_fused_copyto.jl    |   5 +-
 lib/ClimaCorePlots/src/ClimaCorePlots.jl |   2 +-
 src/DataLayouts/DataLayouts.jl           | 353 +++++++++++++----------
 src/DataLayouts/broadcast.jl             | 182 ++++++------
 src/DataLayouts/copyto.jl                |  24 +-
 src/DataLayouts/fused_copyto.jl          |  16 +-
 src/DataLayouts/mapreduce.jl             |  20 +-
 src/Grids/spectralelement.jl             |  10 +-
 src/InputOutput/readers.jl               |  16 +-
 src/Limiters/quasimonotone.jl            |   4 +-
 src/Topologies/dss.jl                    |   5 +-
 src/Topologies/dss_transform.jl          |  32 +-
 test/DataLayouts/benchmark_fill.jl       |   8 +-
 test/DataLayouts/cuda.jl                 |  23 +-
 test/DataLayouts/data0d.jl               |  37 ++-
 test/DataLayouts/data1dx.jl              |  17 +-
 test/DataLayouts/data2d.jl               |  40 +--
 test/DataLayouts/data2dx.jl              |  14 +-
 test/DataLayouts/opt_similar.jl          |  12 +-
 test/DataLayouts/unit_copyto.jl          |  24 +-
 test/DataLayouts/unit_fill.jl            |  32 +-
 test/DataLayouts/unit_mapreduce.jl       |  60 ++--
 test/DataLayouts/unit_ndims.jl           |  20 +-
 test/Fields/unit_field.jl                |  12 +-
 test/Spaces/opt_spaces.jl                |   2 +-
 test/Spaces/unit_spaces.jl               |   4 +-
 29 files changed, 546 insertions(+), 487 deletions(-)

diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl
index 8c60b127bf..6b72273018 100644
--- a/examples/hybrid/driver.jl
+++ b/examples/hybrid/driver.jl
@@ -186,6 +186,10 @@ walltime = @elapsed sol = OrdinaryDiffEq.solve!(integrator)
 any(isnan, sol.u[end]) && error("NaNs found in result.")
 
 if is_distributed # replace sol.u on the root processor with the global sol.u
+    global_Y_c_1 =
+        DataLayouts.gather(comms_ctx, Fields.field_values(sol.u[1].c))
+    global_Y_f_1 =
+        DataLayouts.gather(comms_ctx, Fields.field_values(sol.u[1].f))
     if ClimaComms.iamroot(comms_ctx)
         global_h_space = make_horizontal_space(
             horizontal_mesh,
@@ -194,14 +198,10 @@ if is_distributed # replace sol.u on the root processor with the global sol.u
         )
         global_center_space, global_face_space =
             make_hybrid_spaces(global_h_space, z_max, z_elem; z_stretch)
-        global_Y_c_type = Fields.Field{
-            typeof(Fields.field_values(Y.c)),
-            typeof(global_center_space),
-        }
-        global_Y_f_type = Fields.Field{
-            typeof(Fields.field_values(Y.f)),
-            typeof(global_face_space),
-        }
+        global_Y_c_type =
+            Fields.Field{typeof(global_Y_c_1), typeof(global_center_space)}
+        global_Y_f_type =
+            Fields.Field{typeof(global_Y_f_1), typeof(global_face_space)}
         global_Y_type = Fields.FieldVector{
             FT,
             NamedTuple{(:c, :f), Tuple{global_Y_c_type, global_Y_f_type}},
diff --git a/ext/cuda/data_layouts_copyto.jl b/ext/cuda/data_layouts_copyto.jl
index c14405182d..3e53d05c46 100644
--- a/ext/cuda/data_layouts_copyto.jl
+++ b/ext/cuda/data_layouts_copyto.jl
@@ -16,11 +16,10 @@ function knl_copyto!(dest, src)
 end
 
 function Base.copyto!(
-    dest::IJFH{S, Nij},
-    bc::DataLayouts.BroadcastedUnionIJFH{S, Nij},
+    dest::IJFH{S, Nij, Nh},
+    bc::DataLayouts.BroadcastedUnionIJFH{S, Nij, Nh},
     ::ToCUDA,
-) where {S, Nij}
-    _, _, _, _, Nh = size(bc)
+) where {S, Nij, Nh}
     if Nh > 0
         auto_launch!(
             knl_copyto!,
@@ -34,11 +33,10 @@ function Base.copyto!(
 end
 
 function Base.copyto!(
-    dest::VIJFH{S, Nv, Nij},
-    bc::DataLayouts.BroadcastedUnionVIJFH{S, Nv, Nij},
+    dest::VIJFH{S, Nv, Nij, Nh},
+    bc::DataLayouts.BroadcastedUnionVIJFH{S, Nv, Nij, Nh},
     ::ToCUDA,
-) where {S, Nv, Nij}
-    _, _, _, _, Nh = size(bc)
+) where {S, Nv, Nij, Nh}
     if Nv > 0 && Nh > 0
         Nv_per_block = min(Nv, fld(256, Nij * Nij))
         Nv_blocks = cld(Nv, Nv_per_block)
@@ -58,14 +56,13 @@ function Base.copyto!(
     bc::DataLayouts.BroadcastedUnionVF{S, Nv},
     ::ToCUDA,
 ) where {S, Nv}
-    _, _, _, _, Nh = size(dest)
-    if Nv > 0 && Nh > 0
+    if Nv > 0
         auto_launch!(
             knl_copyto!,
             (dest, bc),
             dest;
             threads_s = (1, 1),
-            blocks_s = (Nh, Nv),
+            blocks_s = (1, Nv),
         )
     end
     return dest
@@ -100,8 +97,8 @@ function knl_copyto_flat!(dest::AbstractData, bc)
 end
 
 function cuda_copyto!(dest::AbstractData, bc)
-    (_, _, Nf, Nv, Nh) = DataLayouts.universal_size(dest)
-    if Nv > 0 && Nh > 0 && Nf > 0
+    (_, _, Nv, Nh) = DataLayouts.universal_size(dest)
+    if Nv > 0 && Nh > 0
         auto_launch!(knl_copyto_flat!, (dest, bc), dest; auto = true)
     end
     return dest
@@ -110,12 +107,12 @@ end
 # TODO: can we use CUDA's luanch configuration for all data layouts?
 # Currently, it seems to have a slight performance degredation.
 #! format: off
-# Base.copyto!(dest::IJFH{S, Nij},      bc::DataLayouts.BroadcastedUnionIJFH{S, Nij}, ::ToCUDA) where {S, Nij} = cuda_copyto!(dest, bc)
-Base.copyto!(dest::IFH{S, Ni},        bc::DataLayouts.BroadcastedUnionIFH{S, Ni}, ::ToCUDA) where {S, Ni} = cuda_copyto!(dest, bc)
-Base.copyto!(dest::IJF{S, Nij},       bc::DataLayouts.BroadcastedUnionIJF{S, Nij}, ::ToCUDA) where {S, Nij} = cuda_copyto!(dest, bc)
-Base.copyto!(dest::IF{S, Ni},         bc::DataLayouts.BroadcastedUnionIF{S, Ni}, ::ToCUDA) where {S, Ni} = cuda_copyto!(dest, bc)
-Base.copyto!(dest::VIFH{S, Nv, Ni},   bc::DataLayouts.BroadcastedUnionVIFH{S, Nv, Ni}, ::ToCUDA) where {S, Nv, Ni} = cuda_copyto!(dest, bc)
-# Base.copyto!(dest::VIJFH{S, Nv, Nij}, bc::DataLayouts.BroadcastedUnionVIJFH{S, Nv, Nij}, ::ToCUDA) where {S, Nv, Nij} = cuda_copyto!(dest, bc)
-# Base.copyto!(dest::VF{S, Nv},         bc::DataLayouts.BroadcastedUnionVF{S, Nv}, ::ToCUDA) where {S, Nv} = cuda_copyto!(dest, bc)
-# Base.copyto!(dest::DataF{S},          bc::DataLayouts.BroadcastedUnionDataF{S}, ::ToCUDA) where {S} = cuda_copyto!(dest, bc)
+# Base.copyto!(dest::IJFH{S, Nij},          bc::DataLayouts.BroadcastedUnionIJFH{S, Nij, Nh}, ::ToCUDA) where {S, Nij, Nh} = cuda_copyto!(dest, bc)
+Base.copyto!(dest::IFH{S, Ni, Nh},        bc::DataLayouts.BroadcastedUnionIFH{S, Ni, Nh}, ::ToCUDA) where {S, Ni, Nh} = cuda_copyto!(dest, bc)
+Base.copyto!(dest::IJF{S, Nij},           bc::DataLayouts.BroadcastedUnionIJF{S, Nij}, ::ToCUDA) where {S, Nij} = cuda_copyto!(dest, bc)
+Base.copyto!(dest::IF{S, Ni},             bc::DataLayouts.BroadcastedUnionIF{S, Ni}, ::ToCUDA) where {S, Ni} = cuda_copyto!(dest, bc)
+Base.copyto!(dest::VIFH{S, Nv, Ni, Nh},   bc::DataLayouts.BroadcastedUnionVIFH{S, Nv, Ni, Nh}, ::ToCUDA) where {S, Nv, Ni, Nh} = cuda_copyto!(dest, bc)
+# Base.copyto!(dest::VIJFH{S, Nv, Nij, Nh}, bc::DataLayouts.BroadcastedUnionVIJFH{S, Nv, Nij, Nh}, ::ToCUDA) where {S, Nv, Nij, Nh} = cuda_copyto!(dest, bc)
+# Base.copyto!(dest::VF{S, Nv},             bc::DataLayouts.BroadcastedUnionVF{S, Nv}, ::ToCUDA) where {S, Nv} = cuda_copyto!(dest, bc)
+# Base.copyto!(dest::DataF{S},              bc::DataLayouts.BroadcastedUnionDataF{S}, ::ToCUDA) where {S} = cuda_copyto!(dest, bc)
 #! format: on
diff --git a/ext/cuda/data_layouts_fill.jl b/ext/cuda/data_layouts_fill.jl
index d8aac57230..f140250db2 100644
--- a/ext/cuda/data_layouts_fill.jl
+++ b/ext/cuda/data_layouts_fill.jl
@@ -11,8 +11,8 @@ function knl_fill_flat!(dest::AbstractData, val)
 end
 
 function cuda_fill!(dest::AbstractData, val)
-    (_, _, Nf, Nv, Nh) = DataLayouts.universal_size(dest)
-    if Nv > 0 && Nh > 0 && Nf > 0
+    (_, _, Nv, Nh) = DataLayouts.universal_size(dest)
+    if Nv > 0 && Nh > 0
         auto_launch!(knl_fill_flat!, (dest, val), dest; auto = true)
     end
     return dest
diff --git a/ext/cuda/data_layouts_fused_copyto.jl b/ext/cuda/data_layouts_fused_copyto.jl
index d6390cd36b..a566e69a5f 100644
--- a/ext/cuda/data_layouts_fused_copyto.jl
+++ b/ext/cuda/data_layouts_fused_copyto.jl
@@ -42,10 +42,9 @@ end
 
 function fused_copyto!(
     fmbc::FusedMultiBroadcast,
-    dest1::VIJFH{S, Nv, Nij},
+    dest1::VIJFH{S, Nv, Nij, Nh},
     ::ToCUDA,
-) where {S, Nv, Nij}
-    _, _, _, _, Nh = size(dest1)
+) where {S, Nv, Nij, Nh}
     if Nv > 0 && Nh > 0
         Nv_per_block = min(Nv, fld(256, Nij * Nij))
         Nv_blocks = cld(Nv, Nv_per_block)
diff --git a/lib/ClimaCorePlots/src/ClimaCorePlots.jl b/lib/ClimaCorePlots/src/ClimaCorePlots.jl
index 1fe2dffcdf..35580b7007 100644
--- a/lib/ClimaCorePlots/src/ClimaCorePlots.jl
+++ b/lib/ClimaCorePlots/src/ClimaCorePlots.jl
@@ -425,7 +425,7 @@ function _unfolded_pannel_matrix(field, interpolate)
     # TODO: inefficient memory wise, but good enough for now
     panels = [fill(NaN, (panel_size * dof, panel_size * dof)) for _ in 1:6]
 
-    interpolated_data = DataLayouts.IJFH{FT, interpolate}(Array{FT}, nelem)
+    interpolated_data = DataLayouts.IJFH{FT, interpolate, nelem}(Array{FT})
     field_data = Fields.field_values(field)
 
     Operators.tensor_product!(interpolated_data, field_data, Imat)
diff --git a/src/DataLayouts/DataLayouts.jl b/src/DataLayouts/DataLayouts.jl
index cba9ffeeeb..16c50a29db 100644
--- a/src/DataLayouts/DataLayouts.jl
+++ b/src/DataLayouts/DataLayouts.jl
@@ -52,25 +52,32 @@ abstract type AbstractData{S} end
 Base.size(data::AbstractData, i::Integer) = size(data)[i]
 
 """
-    (Ni, Nj, Nf, Nv, Nh) = universal_size(data::AbstractData)
+    struct UniversalSize{Nv, Nij, Nh} end
+    UniversalSize(data::AbstractData)
 
-Returns dimensions in a universal
-format for all data layouts:
+A struct containing static dimensions, universal to all datalayouts:
  - `Ni` number of spectral element nodal degrees of freedom in first horizontal direction
  - `Nj` number of spectral element nodal degrees of freedom in second horizontal direction
- - `Nf` number of field components
  - `Nv` number of vertical degrees of freedom
  - `Nh` number of horizontal elements
-
-Note: this is similar to `Base.size`, except
-      that `universal_size` does not return 1
-      for the number of field components.
 """
-function universal_size(data::AbstractData)
+struct UniversalSize{Ni, Nj, Nv, Nh} end
+
+@inline function UniversalSize(data::AbstractData)
     s = size(data)
-    return (s[1], s[2], ncomponents(data), s[4], s[5])
+    return UniversalSize{s[1], s[2], s[4], s[5]}()
 end
 
+"""
+    (Ni, Nj, Nv, Nh) = universal_size(data::AbstractData)
+
+A tuple of compile-time known type parameters,
+corresponding to `UniversalSize`.
+"""
+@inline universal_size(::UniversalSize{Ni, Nj, Nv, Nh}) where {Ni, Nj, Nv, Nh} =
+    (Ni, Nj, Nv, Nh)
+@inline universal_size(data::AbstractData) = universal_size(UniversalSize(data))
+
 function Base.show(io::IO, data::AbstractData)
     indent_width = 2
     (rows, cols) = displaysize(io)
@@ -224,48 +231,54 @@ end
 
 A 3D DataLayout. TODO: Add more docs
 """
-struct IJKFVH{S, Nij, Nk, A} <: Data3D{S, Nij, Nk}
+struct IJKFVH{S, Nij, Nk, Nv, Nh, A} <: Data3D{S, Nij, Nk}
     array::A
 end
 
-parent_array_type(::Type{IJKFVH{S, Nij, Nk, A}}) where {S, Nij, Nk, A} = A
+parent_array_type(
+    ::Type{IJKFVH{S, Nij, Nk, Nv, Nh, A}},
+) where {S, Nij, Nk, Nv, Nh, A} = A
 
-function IJKFVH{S, Nij, Nk}(array::AbstractArray{T, 6}) where {S, Nij, Nk, T}
+function IJKFVH{S, Nij, Nk, Nv, Nh}(
+    array::AbstractArray{T, 6},
+) where {S, Nij, Nk, Nv, Nh, T}
+    check_basetype(T, S)
     @assert size(array, 1) == Nij
     @assert size(array, 2) == Nij
     @assert size(array, 3) == Nk
-    check_basetype(T, S)
     @assert size(array, 4) == typesize(T, S)
-    IJKFVH{S, Nij, Nk, typeof(array)}(array)
+    @assert size(array, 5) == Nv
+    @assert size(array, 6) == Nh
+    IJKFVH{S, Nij, Nk, Nv, Nh, typeof(array)}(array)
 end
 
 function replace_basetype(
-    data::IJKFVH{S, Nij, Nk},
+    data::IJKFVH{S, Nij, Nk, Nv, Nh},
     ::Type{T},
-) where {S, Nij, Nk, T}
+) where {S, Nij, Nk, Nv, Nh, T}
     array = parent(data)
     S′ = replace_basetype(eltype(array), T, S)
-    return IJKFVH{S′, Nij, Nk}(similar(array, T))
+    return IJKFVH{S′, Nij, Nk, Nv, Nh}(similar(array, T))
 end
 
 @generated function _property_view(
-    data::IJKFVH{S, Nij, Nk, A},
+    data::IJKFVH{S, Nij, Nk, Nv, Nh, A},
     ::Val{Idx},
-) where {S, Nij, Nk, A, Idx}
+) where {S, Nij, Nk, Nv, Nh, A, Idx}
     SS = fieldtype(S, Idx)
     T = eltype(A)
     offset = fieldtypeoffset(T, S, Val(Idx))
     nbytes = typesize(T, SS)
     field_byterange = (offset + 1):(offset + nbytes)
-    return :(IJKFVH{$SS, $Nij, $Nk}(
+    return :(IJKFVH{$SS, $Nij, $Nk, $Nv, $Nh}(
         @inbounds view(parent(data), :, :, :, $field_byterange, :, :)
     ))
 end
 
 @inline function Base.getproperty(
-    data::IJKFVH{S, Nij, Nk},
+    data::IJKFVH{S, Nij, Nk, Nv, Nh},
     i::Integer,
-) where {S, Nij, Nk}
+) where {S, Nij, Nk, Nv, Nh}
     array = parent(data)
     T = eltype(array)
     SS = fieldtype(S, i)
@@ -273,14 +286,11 @@ end
     nbytes = typesize(T, SS)
     dataview =
         @inbounds view(array, :, :, :, (offset + 1):(offset + nbytes), :, :)
-    IJKFVH{SS, Nij, Nk}(dataview)
+    IJKFVH{SS, Nij, Nk, Nv, Nh}(dataview)
 end
 
-function Base.size(data::IJKFVH{S, Nij, Nk}) where {S, Nij, Nk}
-    Nv = size(parent(data), 5)
-    Nh = size(parent(data), 6)
-    return (Nij, Nij, Nk, Nv, Nh)
-end
+Base.size(data::IJKFVH{S, Nij, Nk, Nv, Nh}) where {S, Nij, Nk, Nv, Nh} =
+    (Nij, Nij, Nk, Nv, Nh)
 
 # ==================
 # Data2D DataLayout
@@ -300,40 +310,43 @@ The `ArrayType`-constructor constructs a IJFH 2D Spectral
 DataLayout given the backing `ArrayType`, quadrature degrees
 of freedom `Nij × Nij`, and the number of mesh elements `nelements`.
 """
-struct IJFH{S, Nij, A} <: Data2D{S, Nij}
+struct IJFH{S, Nij, Nh, A} <: Data2D{S, Nij}
     array::A
 end
 
-parent_array_type(::Type{IJFH{S, Nij, A}}) where {S, Nij, A} = A
+parent_array_type(::Type{IJFH{S, Nij, Nh, A}}) where {S, Nij, Nh, A} = A
 
-function IJFH{S, Nij}(array::AbstractArray{T, 4}) where {S, Nij, T}
+function IJFH{S, Nij, Nh}(array::AbstractArray{T, 4}) where {S, Nij, Nh, T}
+    check_basetype(T, S)
     @assert size(array, 1) == Nij
     @assert size(array, 2) == Nij
-    check_basetype(T, S)
     @assert size(array, 3) == typesize(T, S)
-    IJFH{S, Nij, typeof(array)}(array)
+    @assert size(array, 4) == Nh
+    IJFH{S, Nij, Nh, typeof(array)}(array)
 end
 
-rebuild(data::IJFH{S, Nij}, array::A) where {S, Nij, A <: AbstractArray} =
-    IJFH{S, Nij}(array)
+rebuild(
+    data::IJFH{S, Nij, Nh},
+    array::A,
+) where {S, Nij, Nh, A <: AbstractArray} = IJFH{S, Nij, Nh}(array)
 
-Base.copy(data::IJFH{S, Nij}) where {S, Nij} = IJFH{S, Nij}(copy(parent(data)))
+Base.copy(data::IJFH{S, Nij, Nh}) where {S, Nij, Nh} =
+    IJFH{S, Nij, Nh}(copy(parent(data)))
 
-function Base.size(data::IJFH{S, Nij}) where {S, Nij}
-    Nv = 1
-    Nh = size(parent(data), 4)
-    (Nij, Nij, 1, Nv, Nh)
-end
+Base.size(data::IJFH{S, Nij, Nh}) where {S, Nij, Nh} = (Nij, Nij, 1, 1, Nh)
 
-function IJFH{S, Nij}(::Type{ArrayType}, nelements) where {S, Nij, ArrayType}
+function IJFH{S, Nij, Nh}(::Type{ArrayType}) where {S, Nij, Nh, ArrayType}
     T = eltype(ArrayType)
-    IJFH{S, Nij}(ArrayType(undef, Nij, Nij, typesize(T, S), nelements))
+    IJFH{S, Nij, Nh}(ArrayType(undef, Nij, Nij, typesize(T, S), Nh))
 end
 
-function replace_basetype(data::IJFH{S, Nij}, ::Type{T}) where {S, Nij, T}
+function replace_basetype(
+    data::IJFH{S, Nij, Nh},
+    ::Type{T},
+) where {S, Nij, Nh, T}
     array = parent(data)
     S′ = replace_basetype(eltype(array), T, S)
-    return IJFH{S′, Nij}(similar(array, T))
+    return IJFH{S′, Nij, Nh}(similar(array, T))
 end
 @propagate_inbounds function Base.getindex(
     data::IJFH{S},
@@ -372,27 +385,30 @@ end
 Base.length(data::IJFH) = size(parent(data), 4)
 
 @generated function _property_view(
-    data::IJFH{S, Nij, A},
+    data::IJFH{S, Nij, Nh, A},
     ::Val{Idx},
-) where {S, Nij, A, Idx}
+) where {S, Nij, Nh, A, Idx}
     SS = fieldtype(S, Idx)
     T = eltype(A)
     offset = fieldtypeoffset(T, S, Val(Idx))
     nbytes = typesize(T, SS)
     field_byterange = (offset + 1):(offset + nbytes)
-    return :(IJFH{$SS, $Nij}(
+    return :(IJFH{$SS, $Nij, $Nh}(
         @inbounds view(parent(data), :, :, $field_byterange, :)
     ))
 end
 
-@inline function Base.getproperty(data::IJFH{S, Nij}, i::Integer) where {S, Nij}
+@inline function Base.getproperty(
+    data::IJFH{S, Nij, Nh},
+    i::Integer,
+) where {S, Nij, Nh}
     array = parent(data)
     T = eltype(array)
     SS = fieldtype(S, i)
     offset = fieldtypeoffset(T, S, i)
     nbytes = typesize(T, SS)
     dataview = @inbounds view(array, :, :, (offset + 1):(offset + nbytes), :)
-    IJFH{SS, Nij}(dataview)
+    IJFH{SS, Nij, Nh}(dataview)
 end
 
 @inline function slab(data::IJFH{S, Nij}, h::Integer) where {S, Nij}
@@ -423,7 +439,8 @@ function gather(
 ) where {S, Nij}
     gatherdata = ClimaComms.gather(ctx, parent(data))
     if ClimaComms.iamroot(ctx)
-        IJFH{S, Nij}(gatherdata)
+        Nh = size(gatherdata, 4)
+        IJFH{S, Nij, Nh}(gatherdata)
     else
         nothing
     end
@@ -436,8 +453,8 @@ end
 Base.length(data::Data1D) = size(parent(data), 3)
 
 """
-    IFH{S, Ni, A} <: Data1D{S, Ni}
-    IFH{S,Ni}(ArrayType, nelements)
+    IFH{S,Ni,Nh,A} <: Data1D{S, Ni}
+    IFH{S,Ni,Nh}(ArrayType)
 
 Backing `DataLayout` for 1D spectral element slabs.
 
@@ -448,42 +465,40 @@ datatype `S` struct field (F), for each 1D mesh element (H).
 The `ArrayType`-constructor makes a IFH 1D Spectral
 DataLayout given the backing `ArrayType`, quadrature
 degrees of freedom `Ni`, and the number of mesh elements
-`nelements`.
+`Nh`.
 """
-struct IFH{S, Ni, A} <: Data1D{S, Ni}
+struct IFH{S, Ni, Nh, A} <: Data1D{S, Ni}
     array::A
 end
 
-parent_array_type(::Type{IFH{S, Ni, A}}) where {S, Ni, A} = A
+parent_array_type(::Type{IFH{S, Ni, Nh, A}}) where {S, Ni, Nh, A} = A
 
-function IFH{S, Ni}(array::AbstractArray{T, 3}) where {S, Ni, T}
-    @assert size(array, 1) == Ni
+function IFH{S, Ni, Nh}(array::AbstractArray{T, 3}) where {S, Ni, Nh, T}
     check_basetype(T, S)
+    @assert size(array, 1) == Ni
     @assert size(array, 2) == typesize(T, S)
-    IFH{S, Ni, typeof(array)}(array)
+    @assert size(array, 3) == Nh
+    IFH{S, Ni, Nh, typeof(array)}(array)
 end
 
-function replace_basetype(data::IFH{S, Ni}, ::Type{T}) where {S, Ni, T}
+function replace_basetype(data::IFH{S, Ni, Nh}, ::Type{T}) where {S, Ni, Nh, T}
     array = parent(data)
     S′ = replace_basetype(eltype(array), T, S)
-    return IFH{S′, Ni}(similar(array, T))
+    return IFH{S′, Ni, Nh}(similar(array, T))
 end
 
-function IFH{S, Ni}(::Type{ArrayType}, nelements) where {S, Ni, ArrayType}
+function IFH{S, Ni, Nh}(::Type{ArrayType}) where {S, Ni, Nh, ArrayType}
     T = eltype(ArrayType)
-    IFH{S, Ni}(ArrayType(undef, Ni, typesize(T, S), nelements))
+    IFH{S, Ni, Nh}(ArrayType(undef, Ni, typesize(T, S), Nh))
 end
 
-rebuild(data::IFH{S, Ni}, array::AbstractArray{T, 3}) where {S, Ni, T} =
-    IFH{S, Ni}(array)
+rebuild(data::IFH{S, Ni, Nh}, array::AbstractArray{T, 3}) where {S, Ni, Nh, T} =
+    IFH{S, Ni, Nh}(array)
 
-Base.copy(data::IFH{S, Ni}) where {S, Ni} = IFH{S, Ni}(copy(parent(data)))
+Base.copy(data::IFH{S, Ni, Nh}) where {S, Ni, Nh} =
+    IFH{S, Ni, Nh}(copy(parent(data)))
 
-function Base.size(data::IFH{S, Ni}) where {S, Ni}
-    Nv = 1
-    Nh = size(parent(data), 3)
-    (Ni, 1, 1, Nv, Nh)
-end
+Base.size(data::IFH{S, Ni, Nh}) where {S, Ni, Nh} = (Ni, 1, 1, 1, Nh)
 
 @inline function slab(data::IFH{S, Ni}, h::Integer) where {S, Ni}
     @boundscheck (1 <= h <= size(parent(data), 3)) ||
@@ -503,27 +518,30 @@ Base.@propagate_inbounds column(data::IFH{S, Ni}, i, j, h) where {S, Ni} =
     column(data, i, h)
 
 @generated function _property_view(
-    data::IFH{S, Ni, A},
+    data::IFH{S, Ni, Nh, A},
     ::Val{Idx},
-) where {S, Ni, A, Idx}
+) where {S, Ni, Nh, A, Idx}
     SS = fieldtype(S, Idx)
     T = eltype(A)
     offset = fieldtypeoffset(T, S, Val(Idx))
     nbytes = typesize(T, SS)
     field_byterange = (offset + 1):(offset + nbytes)
-    return :(IFH{$SS, $Ni}(
+    return :(IFH{$SS, $Ni, $Nh}(
         @inbounds view(parent(data), :, $field_byterange, :)
     ))
 end
 
-@inline function Base.getproperty(data::IFH{S, Ni}, f::Integer) where {S, Ni}
+@inline function Base.getproperty(
+    data::IFH{S, Ni, Nh},
+    f::Integer,
+) where {S, Ni, Nh}
     array = parent(data)
     T = eltype(array)
     SS = fieldtype(S, f)
     offset = fieldtypeoffset(T, S, f)
     nbytes = typesize(T, SS)
     dataview = @inbounds view(array, :, (offset + 1):(offset + nbytes), :)
-    IFH{SS, Ni}(dataview)
+    IFH{SS, Ni, Nh}(dataview)
 end
 
 @inline function Base.getindex(data::IFH{S}, i, _, _, _, h) where {S}
@@ -1038,79 +1056,78 @@ Backing `DataLayout` for 2D spectral element slab + extruded 1D FV column data.
 Column levels (V) are contiguous for every element nodal point (I, J)
 for each `S` datatype struct field (F), for each 2D mesh element slab (H).
 """
-struct VIJFH{S, Nv, Nij, A} <: Data2DX{S, Nv, Nij}
+struct VIJFH{S, Nv, Nij, Nh, A} <: Data2DX{S, Nv, Nij}
     array::A
 end
 
-parent_array_type(::Type{VIJFH{S, Nv, Nij, A}}) where {S, Nv, Nij, A} = A
+parent_array_type(::Type{VIJFH{S, Nv, Nij, Nh, A}}) where {S, Nv, Nij, Nh, A} =
+    A
 
-function VIJFH{S, Nv, Nij}(array::AbstractArray{T, 5}) where {S, Nv, Nij, T}
+function VIJFH{S, Nv, Nij, Nh}(
+    array::AbstractArray{T, 5},
+) where {S, Nv, Nij, Nh, T}
+    check_basetype(T, S)
     @assert size(array, 2) == size(array, 3) == Nij
     @assert size(array, 4) == typesize(T, S)
-    check_basetype(T, S)
-    VIJFH{S, Nv, Nij, typeof(array)}(array)
+    @assert size(array, 5) == Nh
+    VIJFH{S, Nv, Nij, Nh, typeof(array)}(array)
 end
 
 rebuild(
-    data::VIJFH{S, Nv, Nij},
+    data::VIJFH{S, Nv, Nij, Nh},
     array::AbstractArray{T, 5},
-) where {S, Nv, Nij, T} = VIJFH{S, Nv, Nij}(array)
+) where {S, Nv, Nij, Nh, T} = VIJFH{S, Nv, Nij, Nh}(array)
 
 nlevels(::VIJFH{S, Nv}) where {S, Nv} = Nv
 
 function replace_basetype(
-    data::VIJFH{S, Nv, Nij},
+    data::VIJFH{S, Nv, Nij, Nh},
     ::Type{T},
-) where {S, Nv, Nij, T}
+) where {S, Nv, Nij, Nh, T}
     array = parent(data)
     S′ = replace_basetype(eltype(array), T, S)
-    return VIJFH{S′, Nv, Nij}(similar(array, T))
+    return VIJFH{S′, Nv, Nij, Nh}(similar(array, T))
 end
 
-function Base.copy(data::VIJFH{S, Nv, Nij}) where {S, Nv, Nij}
-    VIJFH{S, Nv, Nij}(copy(parent(data)))
+function Base.copy(data::VIJFH{S, Nv, Nij, Nh}) where {S, Nv, Nij, Nh}
+    VIJFH{S, Nv, Nij, Nh}(copy(parent(data)))
 end
 
-function Base.size(data::VIJFH{<:Any, Nv, Nij}) where {Nv, Nij}
-    Nh = size(parent(data), 5)
-    return (Nij, Nij, 1, Nv, Nh)
-end
+Base.size(data::VIJFH{<:Any, Nv, Nij, Nh}) where {Nv, Nij, Nh} =
+    (Nij, Nij, 1, Nv, Nh)
 
-function Base.length(data::VIJFH)
-    size(parent(data), 1) * size(parent(data), 5)
-end
+Base.length(data::VIJFH) = size(parent(data), 1) * size(parent(data), 5)
 
 @generated function _property_view(
-    data::VIJFH{S, Nv, Nij, A},
+    data::VIJFH{S, Nv, Nij, Nh, A},
     ::Val{Idx},
-) where {S, Nv, Nij, A, Idx}
+) where {S, Nv, Nij, Nh, A, Idx}
     SS = fieldtype(S, Idx)
     T = eltype(A)
     offset = fieldtypeoffset(T, S, Val(Idx))
     nbytes = typesize(T, SS)
     field_byterange = (offset + 1):(offset + nbytes)
-    return :(VIJFH{$SS, $Nv, $Nij}(
+    return :(VIJFH{$SS, $Nv, $Nij, $Nh}(
         @inbounds view(parent(data), :, :, :, $field_byterange, :)
     ))
 end
 
 @propagate_inbounds function Base.getproperty(
-    data::VIJFH{S, Nv, Nij},
+    data::VIJFH{S, Nv, Nij, Nh},
     i::Integer,
-) where {S, Nv, Nij}
+) where {S, Nv, Nij, Nh}
     array = parent(data)
     T = eltype(array)
     SS = fieldtype(S, i)
     offset = fieldtypeoffset(T, S, i)
     nbytes = typesize(T, SS)
     dataview = @inbounds view(array, :, :, :, (offset + 1):(offset + nbytes), :)
-    VIJFH{SS, Nv, Nij}(dataview)
+    VIJFH{SS, Nv, Nij, Nh}(dataview)
 end
 
 # Note: construct the subarray view directly as optimizer fails in Base.to_indices (v1.7)
-@inline function slab(data::VIJFH{S, Nv, Nij}, v, h) where {S, Nv, Nij}
+@inline function slab(data::VIJFH{S, Nv, Nij, Nh}, v, h) where {S, Nv, Nij, Nh}
     array = parent(data)
-    Nh = size(array, 5)
     @boundscheck (1 <= v <= Nv && 1 <= h <= Nh) ||
                  throw(BoundsError(data, (v, h)))
     Nf = size(array, 4)
@@ -1126,9 +1143,13 @@ end
 end
 
 # Note: construct the subarray view directly as optimizer fails in Base.to_indices (v1.7)
-@inline function column(data::VIJFH{S, Nv, Nij}, i, j, h) where {S, Nv, Nij}
+@inline function column(
+    data::VIJFH{S, Nv, Nij, Nh},
+    i,
+    j,
+    h,
+) where {S, Nv, Nij, Nh}
     array = parent(data)
-    Nh = size(array, 5)
     @boundscheck (1 <= i <= Nij && 1 <= j <= Nij && 1 <= h <= Nh) ||
                  throw(BoundsError(data, (i, j, h)))
     Nf = size(array, 4)
@@ -1139,11 +1160,11 @@ end
     VF{S, Nv}(dataview)
 end
 
-@inline function level(data::VIJFH{S, Nv, Nij}, v) where {S, Nv, Nij}
+@inline function level(data::VIJFH{S, Nv, Nij, Nh}, v) where {S, Nv, Nij, Nh}
     array = parent(data)
     @boundscheck (1 <= v <= Nv) || throw(BoundsError(data, (v,)))
     dataview = @inbounds view(array, v, :, :, :, :)
-    IJFH{S, Nij}(dataview)
+    IJFH{S, Nij, Nh}(dataview)
 end
 
 @propagate_inbounds function Base.getindex(
@@ -1174,7 +1195,8 @@ function gather(
 ) where {S, Nv, Nij}
     gatherdata = ClimaComms.gather(ctx, parent(data))
     if ClimaComms.iamroot(ctx)
-        VIJFH{S, Nv, Nij}(gatherdata)
+        Nh = size(gatherdata, 5)
+        VIJFH{S, Nv, Nij, Nh}(gatherdata)
     else
         nothing
     end
@@ -1204,41 +1226,44 @@ Backing `DataLayout` for 1D spectral element slab + extruded 1D FV column data.
 Column levels (V) are contiguous for every element nodal point (I)
 for each datatype `S` struct field (F), for each 1D mesh element slab (H).
 """
-struct VIFH{S, Nv, Ni, A} <: Data1DX{S, Nv, Ni}
+struct VIFH{S, Nv, Ni, Nh, A} <: Data1DX{S, Nv, Ni}
     array::A
 end
 
-parent_array_type(::Type{VIFH{S, Nv, Ni, A}}) where {S, Nv, Ni, A} = A
+parent_array_type(::Type{VIFH{S, Nv, Ni, Nh, A}}) where {S, Nv, Ni, Nh, A} = A
 
-function VIFH{S, Nv, Ni}(array::AbstractArray{T, 4}) where {S, Nv, Ni, T}
+function VIFH{S, Nv, Ni, Nh}(
+    array::AbstractArray{T, 4},
+) where {S, Nv, Ni, Nh, T}
+    check_basetype(T, S)
     @assert size(array, 2) == Ni
     @assert size(array, 3) == typesize(T, S)
-    check_basetype(T, S)
-    VIFH{S, Nv, Ni, typeof(array)}(array)
+    @assert size(array, 4) == Nh
+    VIFH{S, Nv, Ni, Nh, typeof(array)}(array)
 end
 
-rebuild(data::VIFH{S, Nv, Ni}, array::A) where {S, Nv, Ni, A <: AbstractArray} =
-    VIFH{S, Nv, Ni}(array)
+rebuild(
+    data::VIFH{S, Nv, Ni, Nh},
+    array::A,
+) where {S, Nv, Ni, Nh, A <: AbstractArray} = VIFH{S, Nv, Ni, Nh}(array)
 
 nlevels(::VIFH{S, Nv}) where {S, Nv} = Nv
 
-function replace_basetype(data::VIFH{S, Nv, Ni}, ::Type{T}) where {S, Nv, Ni, T}
+function replace_basetype(
+    data::VIFH{S, Nv, Ni, Nh},
+    ::Type{T},
+) where {S, Nv, Ni, Nh, T}
     array = parent(data)
     S′ = replace_basetype(eltype(array), T, S)
-    return VIFH{S′, Nv, Ni}(similar(array, T))
+    return VIFH{S′, Nv, Ni, Nh}(similar(array, T))
 end
 
-Base.copy(data::VIFH{S, Nv, Ni}) where {S, Nv, Ni} =
-    VIFH{S, Nv, Ni}(copy(parent(data)))
+Base.copy(data::VIFH{S, Nv, Ni, Nh}) where {S, Nv, Ni, Nh} =
+    VIFH{S, Nv, Ni, Nh}(copy(parent(data)))
 
-function Base.size(data::VIFH{<:Any, Nv, Ni}) where {Nv, Ni}
-    Nh = size(parent(data), 4)
-    return (Ni, 1, 1, Nv, Nh)
-end
+Base.size(data::VIFH{<:Any, Nv, Ni, Nh}) where {Nv, Ni, Nh} = (Ni, 1, 1, Nv, Nh)
 
-function Base.length(data::VIFH)
-    nlevels(data) * size(parent(data), 4)
-end
+Base.length(data::VIFH) = nlevels(data) * size(parent(data), 4)
 
 @propagate_inbounds function Base.getindex(
     data::VIFH{S},
@@ -1250,36 +1275,35 @@ end
 
 
 @generated function _property_view(
-    data::VIFH{S, Nv, Ni, A},
+    data::VIFH{S, Nv, Ni, Nh, A},
     ::Val{Idx},
-) where {S, Nv, Ni, A, Idx}
+) where {S, Nv, Ni, Nh, A, Idx}
     SS = fieldtype(S, Idx)
     T = eltype(A)
     offset = fieldtypeoffset(T, S, Val(Idx))
     nbytes = typesize(T, SS)
     field_byterange = (offset + 1):(offset + nbytes)
-    return :(VIFH{$SS, $Nv, $Ni}(
+    return :(VIFH{$SS, $Nv, $Ni, $Nh}(
         @inbounds view(parent(data), :, :, $field_byterange, :)
     ))
 end
 
 @inline function Base.getproperty(
-    data::VIFH{S, Nv, Ni},
+    data::VIFH{S, Nv, Ni, Nh},
     i::Integer,
-) where {S, Nv, Ni}
+) where {S, Nv, Ni, Nh}
     array = parent(data)
     T = eltype(array)
     SS = fieldtype(S, i)
     offset = fieldtypeoffset(T, S, i)
     nbytes = typesize(T, SS)
     dataview = @inbounds view(array, :, :, (offset + 1):(offset + nbytes), :)
-    VIFH{SS, Nv, Ni}(dataview)
+    VIFH{SS, Nv, Ni, Nh}(dataview)
 end
 
 # Note: construct the subarray view directly as optimizer fails in Base.to_indices (v1.7)
-@inline function slab(data::VIFH{S, Nv, Ni}, v, h) where {S, Nv, Ni}
+@inline function slab(data::VIFH{S, Nv, Ni, Nh}, v, h) where {S, Nv, Ni, Nh}
     array = parent(data)
-    Nh = size(array, 4)
     @boundscheck (1 <= v <= Nv && 1 <= h <= Nh) ||
                  throw(BoundsError(data, (v, h)))
     Nf = size(array, 3)
@@ -1291,9 +1315,8 @@ end
 end
 
 # Note: construct the subarray view directly as optimizer fails in Base.to_indices (v1.7)
-@inline function column(data::VIFH{S, Nv, Ni}, i, h) where {S, Nv, Ni}
+@inline function column(data::VIFH{S, Nv, Ni, Nh}, i, h) where {S, Nv, Ni, Nh}
     array = parent(data)
-    Nh = size(array, 4)
     @boundscheck (1 <= i <= Ni && 1 <= h <= Nh) ||
                  throw(BoundsError(data, (i, h)))
     Nf = size(array, 3)
@@ -1304,9 +1327,13 @@ end
     VF{S, Nv}(dataview)
 end
 
-@inline function column(data::VIFH{S, Nv, Ni}, i, j, h) where {S, Nv, Ni}
+@inline function column(
+    data::VIFH{S, Nv, Ni, Nh},
+    i,
+    j,
+    h,
+) where {S, Nv, Ni, Nh}
     array = parent(data)
-    Nh = size(array, 4)
     @boundscheck (1 <= i <= Ni && j == 1 && 1 <= h <= Nh) ||
                  throw(BoundsError(data, (i, j, h)))
     Nf = size(array, 3)
@@ -1317,11 +1344,11 @@ end
     VF{S, Nv}(dataview)
 end
 
-@inline function level(data::VIFH{S, Nv, Nij}, v) where {S, Nv, Nij}
+@inline function level(data::VIFH{S, Nv, Nij, Nh}, v) where {S, Nv, Nij, Nh}
     array = parent(data)
     @boundscheck (1 <= v <= Nv) || throw(BoundsError(data, (v,)))
     dataview = @inbounds view(array, v, :, :, :)
-    IFH{S, Nij}(dataview)
+    IFH{S, Nij, Nh}(dataview)
 end
 
 @inline function Base.getindex(data::VIFH{S}, i, _, _, v, h) where {S}
@@ -1459,20 +1486,25 @@ rebuild(data::AbstractData, ::Type{DA}) where {DA} =
 include("broadcast.jl")
 
 
-Adapt.adapt_structure(to, data::IJKFVH{S, Nij, Nk}) where {S, Nij, Nk} =
-    IJKFVH{S, Nij, Nk}(Adapt.adapt(to, parent(data)))
+Adapt.adapt_structure(
+    to,
+    data::IJKFVH{S, Nij, Nk, Nv, Nh},
+) where {S, Nij, Nk, Nv, Nh} =
+    IJKFVH{S, Nij, Nk, Nv, Nh}(Adapt.adapt(to, parent(data)))
 
-Adapt.adapt_structure(to, data::IJFH{S, Nij}) where {S, Nij} =
-    IJFH{S, Nij}(Adapt.adapt(to, parent(data)))
+Adapt.adapt_structure(to, data::IJFH{S, Nij, Nh}) where {S, Nij, Nh} =
+    IJFH{S, Nij, Nh}(Adapt.adapt(to, parent(data)))
 
-Adapt.adapt_structure(to, data::VIJFH{S, Nv, Nij}) where {S, Nv, Nij} =
-    VIJFH{S, Nv, Nij}(Adapt.adapt(to, parent(data)))
+Adapt.adapt_structure(to, data::VIJFH{S, Nv, Nij, Nh}) where {S, Nv, Nij, Nh} =
+    VIJFH{S, Nv, Nij, Nh}(Adapt.adapt(to, parent(data)))
 
-Adapt.adapt_structure(to, data::VIFH{S, Nv, Ni, A}) where {S, Nv, Ni, A} =
-    VIFH{S, Nv, Ni}(Adapt.adapt(to, parent(data)))
+Adapt.adapt_structure(
+    to,
+    data::VIFH{S, Nv, Ni, Nh, A},
+) where {S, Nv, Ni, Nh, A} = VIFH{S, Nv, Ni, Nh}(Adapt.adapt(to, parent(data)))
 
-Adapt.adapt_structure(to, data::IFH{S, Ni}) where {S, Ni} =
-    IFH{S, Ni}(Adapt.adapt(to, parent(data)))
+Adapt.adapt_structure(to, data::IFH{S, Ni, Nh}) where {S, Ni, Nh} =
+    IFH{S, Ni, Nh}(Adapt.adapt(to, parent(data)))
 
 Adapt.adapt_structure(to, data::IJF{S, Nij}) where {S, Nij} =
     IJF{S, Nij}(Adapt.adapt(to, parent(data)))
@@ -1535,22 +1567,25 @@ function array2data end
 
 array2data(array::AbstractArray{T, 1}, ::IF{<:Any, Ni}) where {T, Ni} =
     IF{T, Ni}(reshape(array, Ni, 1))
-array2data(array::AbstractArray{T, 1}, ::IFH{<:Any, Ni}) where {T, Ni} =
-    IFH{T, Ni}(reshape(array, Ni, 1, :))
+array2data(array::AbstractArray{T, 1}, ::IFH{<:Any, Ni, Nh}) where {T, Ni, Nh} =
+    IFH{T, Ni, Nh}(reshape(array, Ni, 1, Nh))
 array2data(array::AbstractArray{T, 1}, ::IJF{<:Any, Nij}) where {T, Nij} =
     IJF{T, Nij}(reshape(array, Nij, Nij, 1))
-array2data(array::AbstractArray{T, 1}, ::IJFH{<:Any, Nij}) where {T, Nij} =
-    IJFH{T, Nij}(reshape(array, Nij, Nij, 1, :))
+array2data(
+    array::AbstractArray{T, 1},
+    ::IJFH{<:Any, Nij, Nh},
+) where {T, Nij, Nh} = IJFH{T, Nij, Nh}(reshape(array, Nij, Nij, 1, Nh))
 array2data(array::AbstractArray{T, 2}, ::VF{<:Any, Nv}) where {T, Nv} =
     VF{T, Nv}(reshape(array, Nv, 1))
 array2data(
     array::AbstractArray{T, 2},
-    ::VIFH{<:Any, Nv, Ni},
-) where {T, Nv, Ni} = VIFH{T, Nv, Ni}(reshape(array, Nv, Ni, 1, :))
+    ::VIFH{<:Any, Nv, Ni, Nh},
+) where {T, Nv, Ni, Nh} = VIFH{T, Nv, Ni, Nh}(reshape(array, Nv, Ni, 1, Nh))
 array2data(
     array::AbstractArray{T, 2},
-    ::VIJFH{<:Any, Nv, Nij},
-) where {T, Nv, Nij} = VIJFH{T, Nv, Nij}(reshape(array, Nv, Nij, Nij, 1, :))
+    ::VIJFH{<:Any, Nv, Nij, Nh},
+) where {T, Nv, Nij, Nh} =
+    VIJFH{T, Nv, Nij, Nh}(reshape(array, Nv, Nij, Nij, 1, Nh))
 
 """
     device_dispatch(data::AbstractData)
diff --git a/src/DataLayouts/broadcast.jl b/src/DataLayouts/broadcast.jl
index 35ef0422bc..9da4d81fc8 100644
--- a/src/DataLayouts/broadcast.jl
+++ b/src/DataLayouts/broadcast.jl
@@ -23,12 +23,12 @@ DataStyle(::Type{VF{S, Nv, A}}) where {S, Nv, A} =
 DataColumnStyle(::Type{VFStyle{Nv, A}}) where {Nv, A} = VFStyle{Nv, A}
 
 abstract type Data1DStyle{Ni} <: DataStyle end
-struct IFHStyle{Ni, A} <: Data1DStyle{Ni} end
-DataStyle(::Type{IFH{S, Ni, A}}) where {S, Ni, A} =
-    IFHStyle{Ni, parent_array_type(A)}()
+struct IFHStyle{Ni, Nh, A} <: Data1DStyle{Ni} end
+DataStyle(::Type{IFH{S, Ni, Nh, A}}) where {S, Ni, Nh, A} =
+    IFHStyle{Ni, Nh, parent_array_type(A)}()
 
 abstract type DataSlab1DStyle{Ni} <: DataStyle end
-DataSlab1DStyle(::Type{IFHStyle{Ni, A}}) where {Ni, A} = IFStyle{Ni, A}
+DataSlab1DStyle(::Type{IFHStyle{Ni, Nh, A}}) where {Ni, Nh, A} = IFStyle{Ni, A}
 
 struct IFStyle{Ni, A} <: DataSlab1DStyle{Ni} end
 DataStyle(::Type{IF{S, Ni, A}}) where {S, Ni, A} =
@@ -40,29 +40,32 @@ DataStyle(::Type{IJF{S, Nij, A}}) where {S, Nij, A} =
     IJFStyle{Nij, parent_array_type(A)}()
 
 abstract type Data2DStyle{Nij} <: DataStyle end
-struct IJFHStyle{Nij, A} <: Data2DStyle{Nij} end
-DataStyle(::Type{IJFH{S, Nij, A}}) where {S, Nij, A} =
-    IJFHStyle{Nij, parent_array_type(A)}()
-DataSlab2DStyle(::Type{IJFHStyle{Nij, A}}) where {Nij, A} = IJFStyle{Nij, A}
+struct IJFHStyle{Nij, Nh, A} <: Data2DStyle{Nij} end
+DataStyle(::Type{IJFH{S, Nij, Nh, A}}) where {S, Nij, Nh, A} =
+    IJFHStyle{Nij, Nh, parent_array_type(A)}()
+DataSlab2DStyle(::Type{IJFHStyle{Nij, Nh, A}}) where {Nij, Nh, A} =
+    IJFStyle{Nij, A}
 
 abstract type Data1DXStyle{Nv, Ni} <: DataStyle end
-struct VIFHStyle{Nv, Ni, A} <: Data1DXStyle{Nv, Ni} end
-DataStyle(::Type{VIFH{S, Nv, Ni, A}}) where {S, Nv, Ni, A} =
-    VIFHStyle{Nv, Ni, parent_array_type(A)}()
-Data1DXStyle(::Type{VIFHStyle{Nv, Ni, A}}) where {Ni, Nv, A} =
-    VIFHStyle{Nv, Ni, A}
-DataColumnStyle(::Type{VIFHStyle{Nv, Ni, A}}) where {Ni, Nv, A} = VFStyle{Nv, A}
-DataSlab1DStyle(::Type{VIFHStyle{Nv, Ni, A}}) where {Ni, Nv, A} = IFStyle{Ni, A}
+struct VIFHStyle{Nv, Ni, Nh, A} <: Data1DXStyle{Nv, Ni} end
+DataStyle(::Type{VIFH{S, Nv, Ni, Nh, A}}) where {S, Nv, Ni, Nh, A} =
+    VIFHStyle{Nv, Ni, Nh, parent_array_type(A)}()
+Data1DXStyle(::Type{VIFHStyle{Nv, Ni, Nh, A}}) where {Ni, Nv, Nh, A} =
+    VIFHStyle{Nv, Ni, Nh, A}
+DataColumnStyle(::Type{VIFHStyle{Nv, Ni, Nh, A}}) where {Ni, Nv, Nh, A} =
+    VFStyle{Nv, A}
+DataSlab1DStyle(::Type{VIFHStyle{Nv, Ni, Nh, A}}) where {Ni, Nv, Nh, A} =
+    IFStyle{Ni, A}
 
 abstract type Data2DXStyle{Nv, Nij} <: DataStyle end
-struct VIJFHStyle{Nv, Nij, A} <: Data2DXStyle{Nv, Nij} end
-DataStyle(::Type{VIJFH{S, Nv, Nij, A}}) where {S, Nv, Nij, A} =
-    VIJFHStyle{Nv, Nij, parent_array_type(A)}()
-Data2DXStyle(::Type{VIJFHStyle{Nv, Nij, A}}) where {Nv, Nij, A} =
-    VIJFHStyle{Nv, Nij, A}
-DataColumnStyle(::Type{VIJFHStyle{Nv, Nij, A}}) where {Nv, Nij, A} =
+struct VIJFHStyle{Nv, Nij, Nh, A} <: Data2DXStyle{Nv, Nij} end
+DataStyle(::Type{VIJFH{S, Nv, Nij, Nh, A}}) where {S, Nv, Nij, Nh, A} =
+    VIJFHStyle{Nv, Nij, Nh, parent_array_type(A)}()
+Data2DXStyle(::Type{VIJFHStyle{Nv, Nij, Nh, A}}) where {Nv, Nij, Nh, A} =
+    VIJFHStyle{Nv, Nij, Nh, A}
+DataColumnStyle(::Type{VIJFHStyle{Nv, Nij, Nh, A}}) where {Nv, Nij, Nh, A} =
     VFStyle{Nv, A}
-DataSlab2DStyle(::Type{VIJFHStyle{Nv, Nij, A}}) where {Nv, Nij, A} =
+DataSlab2DStyle(::Type{VIJFHStyle{Nv, Nij, Nh, A}}) where {Nv, Nij, Nh, A} =
     IJFStyle{Nij, A}
 
 #####
@@ -70,14 +73,14 @@ DataSlab2DStyle(::Type{VIJFHStyle{Nv, Nij, A}}) where {Nv, Nij, A} =
 #####
 
 #! format: off
-const BroadcastedUnionIJFH{S, Nij, A}      = Union{Base.Broadcast.Broadcasted{IJFHStyle{Nij, A}}, IJFH{S, Nij, A}}
-const BroadcastedUnionIFH{S, Ni, A}        = Union{Base.Broadcast.Broadcasted{IFHStyle{Ni, A}}, IFH{S, Ni, A}}
-const BroadcastedUnionIJF{S, Nij, A}       = Union{Base.Broadcast.Broadcasted{IJFStyle{Nij, A}}, IJF{S, Nij, A}}
-const BroadcastedUnionIF{S, Ni, A}         = Union{Base.Broadcast.Broadcasted{IFStyle{Ni, A}}, IF{S, Ni, A}}
-const BroadcastedUnionVIFH{S, Nv, Ni, A}   = Union{Base.Broadcast.Broadcasted{VIFHStyle{Nv, Ni, A}}, VIFH{S, Nv, Ni, A}}
-const BroadcastedUnionVIJFH{S, Nv, Nij, A} = Union{Base.Broadcast.Broadcasted{VIJFHStyle{Nv, Nij, A}}, VIJFH{S, Nv, Nij, A}}
-const BroadcastedUnionVF{S, Nv, A}         = Union{Base.Broadcast.Broadcasted{VFStyle{Nv, A}}, VF{S, Nv, A}}
-const BroadcastedUnionDataF{S, A}          = Union{Base.Broadcast.Broadcasted{DataFStyle{A}}, DataF{S, A}}
+const BroadcastedUnionIJFH{S, Nij, Nh, A}      = Union{Base.Broadcast.Broadcasted{IJFHStyle{Nij, Nh, A}}, IJFH{S, Nij, Nh, A}}
+const BroadcastedUnionIFH{S, Ni, Nh, A}        = Union{Base.Broadcast.Broadcasted{IFHStyle{Ni, Nh, A}}, IFH{S, Ni, Nh, A}}
+const BroadcastedUnionIJF{S, Nij, A}           = Union{Base.Broadcast.Broadcasted{IJFStyle{Nij, A}}, IJF{S, Nij, A}}
+const BroadcastedUnionIF{S, Ni, A}             = Union{Base.Broadcast.Broadcasted{IFStyle{Ni, A}}, IF{S, Ni, A}}
+const BroadcastedUnionVIFH{S, Nv, Ni, Nh, A}   = Union{Base.Broadcast.Broadcasted{VIFHStyle{Nv, Ni, Nh, A}}, VIFH{S, Nv, Ni, Nh, A}}
+const BroadcastedUnionVIJFH{S, Nv, Nij, Nh, A} = Union{Base.Broadcast.Broadcasted{VIJFHStyle{Nv, Nij, Nh, A}}, VIJFH{S, Nv, Nij, Nh, A}}
+const BroadcastedUnionVF{S, Nv, A}             = Union{Base.Broadcast.Broadcasted{VFStyle{Nv, A}}, VF{S, Nv, A}}
+const BroadcastedUnionDataF{S, A}              = Union{Base.Broadcast.Broadcasted{DataFStyle{A}}, DataF{S, A}}
 #! format: on
 
 abstract type Data3DStyle <: DataStyle end
@@ -108,27 +111,28 @@ Base.Broadcast.BroadcastStyle(
     ::IFStyle{Ni, A2},
 ) where {Ni, A1, A2} = IFStyle{Ni, promote_parent_array_type(A1, A2)}()
 Base.Broadcast.BroadcastStyle(
-    ::IFHStyle{Ni, A1},
-    ::IFHStyle{Ni, A2},
-) where {Ni, A1, A2} = IFHStyle{Ni, promote_parent_array_type(A1, A2)}()
+    ::IFHStyle{Ni, Nh, A1},
+    ::IFHStyle{Ni, Nh, A2},
+) where {Ni, Nh, A1, A2} = IFHStyle{Ni, Nh, promote_parent_array_type(A1, A2)}()
 Base.Broadcast.BroadcastStyle(
-    ::VIFHStyle{Nv, Ni, A1},
-    ::VIFHStyle{Nv, Ni, A2},
-) where {Nv, Ni, A1, A2} =
-    VIFHStyle{Nv, Ni, promote_parent_array_type(A1, A2)}()
+    ::VIFHStyle{Nv, Ni, Nh, A1},
+    ::VIFHStyle{Nv, Ni, Nh, A2},
+) where {Nv, Ni, Nh, A1, A2} =
+    VIFHStyle{Nv, Ni, Nh, promote_parent_array_type(A1, A2)}()
 Base.Broadcast.BroadcastStyle(
     ::IJFStyle{Nij, A1},
     ::IJFStyle{Nij, A2},
 ) where {Nij, A1, A2} = IJFStyle{Nij, promote_parent_array_type(A1, A2)}()
 Base.Broadcast.BroadcastStyle(
-    ::IJFHStyle{Nij, A1},
-    ::IJFHStyle{Nij, A2},
-) where {Nij, A1, A2} = IJFHStyle{Nij, promote_parent_array_type(A1, A2)}()
+    ::IJFHStyle{Nij, Nh, A1},
+    ::IJFHStyle{Nij, Nh, A2},
+) where {Nij, Nh, A1, A2} =
+    IJFHStyle{Nij, Nh, promote_parent_array_type(A1, A2)}()
 Base.Broadcast.BroadcastStyle(
-    ::VIJFHStyle{Nv, Nij, A1},
-    ::VIJFHStyle{Nv, Nij, A2},
-) where {Nv, Nij, A1, A2} =
-    VIJFHStyle{Nv, Nij, promote_parent_array_type(A1, A2)}()
+    ::VIJFHStyle{Nv, Nij, Nh, A1},
+    ::VIJFHStyle{Nv, Nij, Nh, A2},
+) where {Nv, Nij, Nh, A1, A2} =
+    VIJFHStyle{Nv, Nij, Nh, promote_parent_array_type(A1, A2)}()
 
 Base.Broadcast.BroadcastStyle(
     ::DataFStyle{A1},
@@ -147,61 +151,62 @@ Base.Broadcast.BroadcastStyle(
 
 Base.Broadcast.BroadcastStyle(
     ::DataFStyle{A1},
-    ::IFHStyle{Ni, A2},
-) where {Ni, A1, A2} = IFHStyle{Ni, promote_parent_array_type(A1, A2)}()
+    ::IFHStyle{Ni, Nh, A2},
+) where {Ni, Nh, A1, A2} = IFHStyle{Ni, Nh, promote_parent_array_type(A1, A2)}()
 
 Base.Broadcast.BroadcastStyle(
     ::DataFStyle{A1},
-    ::IJFHStyle{Nij, A2},
-) where {Nij, A1, A2} = IJFHStyle{Nij, promote_parent_array_type(A1, A2)}()
+    ::IJFHStyle{Nij, Nh, A2},
+) where {Nij, Nh, A1, A2} =
+    IJFHStyle{Nij, Nh, promote_parent_array_type(A1, A2)}()
 
 Base.Broadcast.BroadcastStyle(
     ::DataFStyle{A1},
-    ::VIFHStyle{Nv, Ni, A2},
-) where {Nv, Ni, A1, A2} =
-    VIFHStyle{Nv, Ni, promote_parent_array_type(A1, A2)}()
+    ::VIFHStyle{Nv, Ni, Nh, A2},
+) where {Nv, Ni, Nh, A1, A2} =
+    VIFHStyle{Nv, Ni, Nh, promote_parent_array_type(A1, A2)}()
 
 Base.Broadcast.BroadcastStyle(
     ::DataFStyle{A1},
-    ::VIJFHStyle{Nv, Nij, A2},
-) where {Nv, Nij, A1, A2} =
-    VIJFHStyle{Nv, Nij, promote_parent_array_type(A1, A2)}()
+    ::VIJFHStyle{Nv, Nij, Nh, A2},
+) where {Nv, Nij, Nh, A1, A2} =
+    VIJFHStyle{Nv, Nij, Nh, promote_parent_array_type(A1, A2)}()
 
 Base.Broadcast.BroadcastStyle(
     ::VFStyle{Nv, A1},
-    ::IFHStyle{Ni, A2},
-) where {Nv, Ni, A1, A2} =
-    VIFHStyle{Nv, Ni, promote_parent_array_type(A1, A2)}()
+    ::IFHStyle{Ni, Nh, A2},
+) where {Nv, Ni, Nh, A1, A2} =
+    VIFHStyle{Nv, Ni, Nh, promote_parent_array_type(A1, A2)}()
 
 Base.Broadcast.BroadcastStyle(
     ::VFStyle{Nv, A1},
-    ::IJFHStyle{Nij, A2},
-) where {Nv, Nij, A1, A2} =
-    VIJFHStyle{Nv, Nij, promote_parent_array_type(A1, A2)}()
+    ::IJFHStyle{Nij, Nh, A2},
+) where {Nv, Nij, Nh, A1, A2} =
+    VIJFHStyle{Nv, Nij, Nh, promote_parent_array_type(A1, A2)}()
 
 Base.Broadcast.BroadcastStyle(
     ::VFStyle{Nv, A1},
-    ::VIFHStyle{Nv, Ni, A2},
-) where {Nv, Ni, A1, A2} =
-    VIFHStyle{Nv, Ni, promote_parent_array_type(A1, A2)}()
+    ::VIFHStyle{Nv, Ni, Nh, A2},
+) where {Nv, Ni, Nh, A1, A2} =
+    VIFHStyle{Nv, Ni, Nh, promote_parent_array_type(A1, A2)}()
 
 Base.Broadcast.BroadcastStyle(
     ::VFStyle{Nv, A1},
-    ::VIJFHStyle{Nv, Nij, A2},
-) where {Nv, Nij, A1, A2} =
-    VIJFHStyle{Nv, Nij, promote_parent_array_type(A1, A2)}()
+    ::VIJFHStyle{Nv, Nij, Nh, A2},
+) where {Nv, Nij, Nh, A1, A2} =
+    VIJFHStyle{Nv, Nij, Nh, promote_parent_array_type(A1, A2)}()
 
 Base.Broadcast.BroadcastStyle(
-    ::IFHStyle{Ni, A1},
-    ::VIFHStyle{Nv, Ni, A2},
-) where {Nv, Ni, A1, A2} =
-    VIFHStyle{Nv, Ni, promote_parent_array_type(A1, A2)}()
+    ::IFHStyle{Ni, Nh, A1},
+    ::VIFHStyle{Nv, Ni, Nh, A2},
+) where {Nv, Ni, Nh, A1, A2} =
+    VIFHStyle{Nv, Ni, Nh, promote_parent_array_type(A1, A2)}()
 
 Base.Broadcast.BroadcastStyle(
-    ::IJFHStyle{Nij, A1},
-    ::VIJFHStyle{Nv, Nij, A2},
-) where {Nv, Nij, A1, A2} =
-    VIJFHStyle{Nv, Nij, promote_parent_array_type(A1, A2)}()
+    ::IJFHStyle{Nij, Nh, A1},
+    ::VIJFHStyle{Nv, Nij, Nh, A2},
+) where {Nv, Nij, Nh, A1, A2} =
+    VIJFHStyle{Nv, Nij, Nh, promote_parent_array_type(A1, A2)}()
 
 Base.Broadcast.broadcastable(data::AbstractData) = data
 
@@ -292,23 +297,21 @@ function Base.similar(
 end
 
 function Base.similar(
-    bc::BroadcastedUnionIJFH{<:Any, Nij, A},
+    bc::BroadcastedUnionIJFH{<:Any, Nij, Nh, A},
     ::Type{Eltype},
-    (_, _, _, _, Nh) = size(bc),
-) where {Nij, A, Eltype}
+) where {Nij, Nh, A, Eltype}
     PA = parent_array_type(A)
     array = similar(PA, (Nij, Nij, typesize(eltype(A), Eltype), Nh))
-    return IJFH{Eltype, Nij}(array)
+    return IJFH{Eltype, Nij, Nh}(array)
 end
 
 function Base.similar(
-    bc::BroadcastedUnionIFH{<:Any, Ni, A},
+    bc::BroadcastedUnionIFH{<:Any, Ni, Nh, A},
     ::Type{Eltype},
-    (_, _, _, _, Nh) = size(bc),
-) where {Ni, A, Eltype}
+) where {Ni, Nh, A, Eltype}
     PA = parent_array_type(A)
     array = similar(PA, (Ni, typesize(eltype(A), Eltype), Nh))
-    return IFH{Eltype, Ni}(array)
+    return IFH{Eltype, Ni, Nh}(array)
 end
 
 function Base.similar(
@@ -339,7 +342,6 @@ function Base.similar(
     ::Type{Eltype},
     ::Val{newNv},
 ) where {Nv, A, Eltype, newNv}
-    (_, _, _, _, Nh) = size(bc)
     PA = parent_array_type(A)
     array = similar(PA, (newNv, typesize(eltype(A), Eltype)))
     return VF{Eltype, newNv}(array)
@@ -351,14 +353,13 @@ Base.similar(
 ) where {Nv, Eltype} = Base.similar(bc, Eltype, Val(Nv))
 
 function Base.similar(
-    bc::BroadcastedUnionVIFH{<:Any, Nv, Ni, A},
+    bc::BroadcastedUnionVIFH{<:Any, Nv, Ni, Nh, A},
     ::Type{Eltype},
     ::Val{newNv},
-) where {Nv, Ni, A, Eltype, newNv}
-    (_, _, _, _, Nh) = size(bc)
+) where {Nv, Ni, Nh, A, Eltype, newNv}
     PA = parent_array_type(A)
     array = similar(PA, (newNv, Ni, typesize(eltype(A), Eltype), Nh))
-    return VIFH{Eltype, newNv, Ni}(array)
+    return VIFH{Eltype, newNv, Ni, Nh}(array)
 end
 
 Base.similar(
@@ -367,14 +368,13 @@ Base.similar(
 ) where {Nv, Nij, A, Eltype} = similar(bc, Eltype, Val(Nv))
 
 function Base.similar(
-    bc::BroadcastedUnionVIJFH{<:Any, Nv, Nij, A},
+    bc::BroadcastedUnionVIJFH{<:Any, Nv, Nij, Nh, A},
     ::Type{Eltype},
     ::Val{newNv},
-) where {Nv, Nij, A, Eltype, newNv}
-    (_, _, _, _, Nh) = size(bc)
+) where {Nv, Nij, Nh, A, Eltype, newNv}
     PA = parent_array_type(A)
     array = similar(PA, (newNv, Nij, Nij, typesize(eltype(A), Eltype), Nh))
-    return VIJFH{Eltype, newNv, Nij}(array)
+    return VIJFH{Eltype, newNv, Nij, Nh}(array)
 end
 
 # ============= FusedMultiBroadcast
diff --git a/src/DataLayouts/copyto.jl b/src/DataLayouts/copyto.jl
index eae87bc1c9..4a94638edb 100644
--- a/src/DataLayouts/copyto.jl
+++ b/src/DataLayouts/copyto.jl
@@ -47,10 +47,9 @@ end
 
 function Base.copyto!(
     dest::IJFH{S, Nij},
-    bc::BroadcastedUnionIJFH{S, Nij},
+    bc::BroadcastedUnionIJFH{S, Nij, Nh},
     ::ToCPU,
-) where {S, Nij}
-    _, _, _, _, Nh = size(bc)
+) where {S, Nij, Nh}
     @inbounds for h in 1:Nh
         slab_dest = slab(dest, h)
         slab_bc = slab(bc, h)
@@ -61,10 +60,9 @@ end
 
 function Base.copyto!(
     dest::IFH{S, Ni},
-    bc::BroadcastedUnionIFH{S, Ni},
+    bc::BroadcastedUnionIFH{S, Ni, Nh},
     ::ToCPU,
-) where {S, Ni}
-    _, _, _, _, Nh = size(bc)
+) where {S, Ni, Nh}
     @inbounds for h in 1:Nh
         slab_dest = slab(dest, h)
         slab_bc = slab(bc, h)
@@ -125,11 +123,10 @@ function Base.copyto!(
 end
 
 function Base.copyto!(
-    dest::VIFH{S, Nv, Ni},
-    bc::BroadcastedUnionVIFH{S, Nv, Ni},
+    dest::VIFH{S, Nv, Ni, Nh},
+    bc::BroadcastedUnionVIFH{S, Nv, Ni, Nh},
     ::ToCPU,
-) where {S, Nv, Ni}
-    (_, _, _, _, Nh) = size(bc)
+) where {S, Nv, Ni, Nh}
     # copy contiguous columns
     @inbounds for h in 1:Nh, i in 1:Ni
         col_dest = column(dest, i, h)
@@ -140,12 +137,11 @@ function Base.copyto!(
 end
 
 function Base.copyto!(
-    dest::VIJFH{S, Nv, Nij},
-    bc::BroadcastedUnionVIJFH{S, Nv, Nij},
+    dest::VIJFH{S, Nv, Nij, Nh},
+    bc::BroadcastedUnionVIJFH{S, Nv, Nij, Nh},
     ::ToCPU,
-) where {S, Nv, Nij}
+) where {S, Nv, Nij, Nh}
     # copy contiguous columns
-    _, _, _, _, Nh = size(dest)
     @inbounds for h in 1:Nh, j in 1:Nij, i in 1:Nij
         col_dest = column(dest, i, j, h)
         col_bc = column(bc, i, j, h)
diff --git a/src/DataLayouts/fused_copyto.jl b/src/DataLayouts/fused_copyto.jl
index 8589811792..5ea26e908e 100644
--- a/src/DataLayouts/fused_copyto.jl
+++ b/src/DataLayouts/fused_copyto.jl
@@ -23,10 +23,9 @@ end
 
 function fused_copyto!(
     fmbc::FusedMultiBroadcast,
-    dest1::VIJFH{S1, Nv1, Nij},
+    dest1::VIJFH{S1, Nv1, Nij, Nh},
     ::ToCPU,
-) where {S1, Nv1, Nij}
-    _, _, _, _, Nh = size(dest1)
+) where {S1, Nv1, Nij, Nh}
     for (dest, bc) in fmbc.pairs
         # Base.copyto!(dest, bc) # we can just fall back like this
         @inbounds for h in 1:Nh, j in 1:Nij, i in 1:Nij, v in 1:Nv1
@@ -40,11 +39,11 @@ end
 
 function fused_copyto!(
     fmbc::FusedMultiBroadcast,
-    dest1::IJFH{S, Nij},
+    dest1::IJFH{S, Nij, Nh},
     ::ToCPU,
-) where {S, Nij}
+) where {S, Nij, Nh}
     # copy contiguous columns
-    _, _, _, Nv, Nh = size(dest1)
+    _, _, _, Nv, _ = size(dest1)
     for (dest, bc) in fmbc.pairs
         @inbounds for h in 1:Nh, j in 1:Nij, i in 1:Nij
             I = CartesianIndex(i, j, 1, 1, h)
@@ -57,11 +56,10 @@ end
 
 function fused_copyto!(
     fmbc::FusedMultiBroadcast,
-    dest1::VIFH{S, Nv1, Ni},
+    dest1::VIFH{S, Nv1, Ni, Nh},
     ::ToCPU,
-) where {S, Nv1, Ni}
+) where {S, Nv1, Ni, Nh}
     # copy contiguous columns
-    _, _, _, _, Nh = size(dest1)
     for (dest, bc) in fmbc.pairs
         @inbounds for h in 1:Nh, i in 1:Ni, v in 1:Nv1
             I = CartesianIndex(i, 1, 1, v, h)
diff --git a/src/DataLayouts/mapreduce.jl b/src/DataLayouts/mapreduce.jl
index 5e63477c09..0f9f02eb0c 100644
--- a/src/DataLayouts/mapreduce.jl
+++ b/src/DataLayouts/mapreduce.jl
@@ -15,10 +15,9 @@ end
 function Base.mapreduce(
     fn::F,
     op::Op,
-    bc::BroadcastedUnionIJFH{<:Any, Nij, A},
-) where {F, Op, Nij, A}
+    bc::BroadcastedUnionIJFH{<:Any, Nij, Nh, A},
+) where {F, Op, Nij, Nh, A}
     # mapreduce across DataSlab2D
-    _, _, _, _, Nh = size(bc)
     mapreduce(op, 1:Nh) do h
         Base.@_inline_meta
         slabview = @inbounds slab(bc, h)
@@ -29,10 +28,9 @@ end
 function Base.mapreduce(
     fn::F,
     op::Op,
-    bc::BroadcastedUnionIFH{<:Any, Ni, A},
-) where {F, Op, Ni, A}
+    bc::BroadcastedUnionIFH{<:Any, Ni, Nh, A},
+) where {F, Op, Ni, Nh, A}
     # mapreduce across DataSlab1D
-    _, _, _, _, Nh = size(bc)
     mapreduce(op, 1:Nh) do h
         Base.@_inline_meta
         slabview = @inbounds slab(bc, h)
@@ -77,10 +75,9 @@ end
 function Base.mapreduce(
     fn::F,
     op::Op,
-    bc::BroadcastedUnionVIFH{<:Any, Nv, Ni, A},
-) where {F, Op, Nv, Ni, A}
+    bc::BroadcastedUnionVIFH{<:Any, Nv, Ni, Nh, A},
+) where {F, Op, Nv, Ni, Nh, A}
     # mapreduce across columns
-    _, _, _, _, Nh = size(bc)
     mapreduce(op, Iterators.product(1:Ni, 1:Nh)) do (i, h)
         Base.@_inline_meta
         columnview = @inbounds column(bc, i, h)
@@ -91,10 +88,9 @@ end
 function Base.mapreduce(
     fn::F,
     op::Op,
-    bc::BroadcastedUnionVIJFH{<:Any, Nv, Nij, A},
-) where {F, Op, Nv, Nij, A}
+    bc::BroadcastedUnionVIJFH{<:Any, Nv, Nij, Nh, A},
+) where {F, Op, Nv, Nij, Nh, A}
     # mapreduce across columns
-    _, _, _, _, Nh = size(bc)
     mapreduce(op, Iterators.product(1:Nij, 1:Nij, 1:Nh)) do (i, j, h)
         Base.@_inline_meta
         columnview = @inbounds column(bc, i, j, h)
diff --git a/src/Grids/spectralelement.jl b/src/Grids/spectralelement.jl
index 47681d3a6f..87434aecc5 100644
--- a/src/Grids/spectralelement.jl
+++ b/src/Grids/spectralelement.jl
@@ -48,10 +48,11 @@ function _SpectralElementGrid1D(
     AIdx = Geometry.coordinate_axis(CoordType)
     FT = eltype(CoordType)
     nelements = Topologies.nlocalelems(topology)
+    Nh = nelements
     Nq = Quadratures.degrees_of_freedom(quadrature_style)
 
     LG = Geometry.LocalGeometry{AIdx, CoordType, FT, SMatrix{1, 1, FT, 1}}
-    local_geometry = DataLayouts.IFH{LG, Nq}(Array{FT}, nelements)
+    local_geometry = DataLayouts.IFH{LG, Nq, Nh}(Array{FT})
     quad_points, quad_weights =
         Quadratures.quadrature_points(FT, quadrature_style)
 
@@ -213,6 +214,7 @@ function _SpectralElementGrid2D(
     CoordType2D = get_CoordType2D(topology)
     AIdx = Geometry.coordinate_axis(CoordType2D)
     nlelems = Topologies.nlocalelems(topology)
+    Nh = nlelems
     ngelems = Topologies.nghostelems(topology)
     Nq = Quadratures.degrees_of_freedom(quadrature_style)
     high_order_quadrature_style = Quadratures.GLL{Nq * 2}()
@@ -220,7 +222,7 @@ function _SpectralElementGrid2D(
 
     LG = Geometry.LocalGeometry{AIdx, CoordType2D, FT, SMatrix{2, 2, FT, 4}}
 
-    local_geometry = DataLayouts.IJFH{LG, Nq}(Array{FT}, nlelems)
+    local_geometry = DataLayouts.IJFH{LG, Nq, Nh}(Array{FT})
 
     quad_points, quad_weights =
         Quadratures.quadrature_points(FT, quadrature_style)
@@ -378,7 +380,7 @@ function _SpectralElementGrid2D(
 
     if quadrature_style isa Quadratures.GLL
         internal_surface_geometry =
-            DataLayouts.IFH{SG, Nq}(Array{FT}, length(interior_faces))
+            DataLayouts.IFH{SG, Nq, length(interior_faces)}(Array{FT})
         for (iface, (lidx⁻, face⁻, lidx⁺, face⁺, reversed)) in
             enumerate(interior_faces)
             internal_surface_geometry_slab =
@@ -417,7 +419,7 @@ function _SpectralElementGrid2D(
                 boundary_faces =
                     Topologies.boundary_faces(topology, boundarytag)
                 boundary_surface_geometry =
-                    DataLayouts.IFH{SG, Nq}(Array{FT}, length(boundary_faces))
+                    DataLayouts.IFH{SG, Nq, length(boundary_faces)}(Array{FT})
                 for (iface, (elem, face)) in enumerate(boundary_faces)
                     boundary_surface_geometry_slab =
                         slab(boundary_surface_geometry, iface)
diff --git a/src/InputOutput/readers.jl b/src/InputOutput/readers.jl
index 537ce1c6f5..fd8b28c197 100644
--- a/src/InputOutput/readers.jl
+++ b/src/InputOutput/readers.jl
@@ -154,6 +154,16 @@ function _scan_data_layout(layoutstring::AbstractString)
     return DataLayouts.VIFH
 end
 
+function Nh_dim(layoutstring::AbstractString)
+    @assert layoutstring ∈ ("IJFH", "IJF", "IFH", "IF", "VIJFH", "VIFH")
+    layoutstring == "IJFH" && return 4
+    layoutstring == "IJF" && return -1
+    layoutstring == "IFH" && return 3
+    layoutstring == "IF" && return -1
+    layoutstring == "VIJFH" && return 5
+    return 4
+end
+
 """
     matrix_to_cartesianindices(elemorder_matrix)
 
@@ -459,12 +469,14 @@ function read_field(reader::HDF5Reader, name::AbstractString)
         data_layout = attrs(obj)["data_layout"]
         Nij = size(data, findfirst("I", data_layout)[1])
         DataLayout = _scan_data_layout(data_layout)
+        Nhd = Nh_dim(data_layout)
+        Nht = Nhd == -1 ? () : (size(data, Nhd),)
         ElType = eval(Meta.parse(attrs(obj)["value_type"]))
         if data_layout in ("VIJFH", "VIFH")
             Nv = size(data, 1)
-            values = DataLayout{ElType, Nv, Nij}(data)
+            values = DataLayout{ElType, Nv, Nij, Nht...}(data)
         else
-            values = DataLayout{ElType, Nij}(data)
+            values = DataLayout{ElType, Nij, Nht...}(data)
         end
         return Fields.Field(values, space)
     elseif type == "FieldVector"
diff --git a/src/Limiters/quasimonotone.jl b/src/Limiters/quasimonotone.jl
index e2b9e3a5a6..24fcc32060 100644
--- a/src/Limiters/quasimonotone.jl
+++ b/src/Limiters/quasimonotone.jl
@@ -71,14 +71,14 @@ function make_q_bounds(
 ) where {S}
     Nf = DataLayouts.ncomponents(ρq)
     _, _, _, _, Nh = size(ρq)
-    return DataLayouts.IFH{S, 2}(similar(parent(ρq), (2, Nf, Nh)))
+    return DataLayouts.IFH{S, 2, Nh}(similar(parent(ρq), (2, Nf, Nh)))
 end
 function make_q_bounds(
     ρq::Union{DataLayouts.VIFH{S}, DataLayouts.VIJFH{S}},
 ) where {S}
     Nf = DataLayouts.ncomponents(ρq)
     _, _, _, Nv, Nh = size(ρq)
-    return DataLayouts.VIFH{S, Nv, 2}(similar(parent(ρq), (Nv, 2, Nf, Nh)))
+    return DataLayouts.VIFH{S, Nv, 2, Nh}(similar(parent(ρq), (Nv, 2, Nf, Nh)))
 end
 
 
diff --git a/src/Topologies/dss.jl b/src/Topologies/dss.jl
index d093dc1f10..dce465412d 100644
--- a/src/Topologies/dss.jl
+++ b/src/Topologies/dss.jl
@@ -10,7 +10,7 @@ struct DSSBuffer{S, G, D, A, B, VI}
     "ClimaComms graph context for communication"
     graph_context::G
     """
-    Perimeter `DataLayout` object: typically a `VIFH{TT,Nv,Np}`, where `TT` is the
+    Perimeter `DataLayout` object: typically a `VIFH{TT,Nv,Np,Nh}`, where `TT` is the
     transformed type, `Nv` is the number of vertical levels, and `Np` is the length of the perimeter
     """
     perimeter_data::D
@@ -74,7 +74,8 @@ function create_dss_buffer(
     if eltype(data) <: Geometry.Covariant123Vector
         TS = Geometry.UVWVector{T}
     end
-    perimeter_data = DataLayouts.VIFH{TS, Nv, Np}(DA{T}(undef, Nv, Np, Nf, Nh))
+    perimeter_data =
+        DataLayouts.VIFH{TS, Nv, Np, Nh}(DA{T}(undef, Nv, Np, Nf, Nh))
     if context isa ClimaComms.SingletonCommsContext
         graph_context = ClimaComms.SingletonGraphContext(context)
         send_data, recv_data = T[], T[]
diff --git a/src/Topologies/dss_transform.jl b/src/Topologies/dss_transform.jl
index d8d0a2c7ce..a03ab3c79c 100644
--- a/src/Topologies/dss_transform.jl
+++ b/src/Topologies/dss_transform.jl
@@ -293,29 +293,29 @@ function create_ghost_buffer(
     topology::Topologies.Topology2D,
 ) where {S, Nij}
     if data isa DataLayouts.IJFH
-        send_data = DataLayouts.IJFH{S, Nij}(
+        send_data = DataLayouts.IJFH{S, Nij, Topologies.nsendelems(topology)}(
             typeof(parent(data)),
-            Topologies.nsendelems(topology),
         )
-        recv_data = DataLayouts.IJFH{S, Nij}(
+        recv_data = DataLayouts.IJFH{S, Nij, Topologies.nrecvelems(topology)}(
             typeof(parent(data)),
-            Topologies.nrecvelems(topology),
         )
         k = stride(parent(send_data), 4)
     else
         Nv, _, _, Nf, _ = size(parent(data))
-        send_data = DataLayouts.VIJFH{S, Nv, Nij}(
-            similar(
-                parent(data),
-                (Nv, Nij, Nij, Nf, Topologies.nsendelems(topology)),
-            ),
-        )
-        recv_data = DataLayouts.VIJFH{S, Nv, Nij}(
-            similar(
-                parent(data),
-                (Nv, Nij, Nij, Nf, Topologies.nrecvelems(topology)),
-            ),
-        )
+        send_data =
+            DataLayouts.VIJFH{S, Nv, Nij, Topologies.nsendelems(topology)}(
+                similar(
+                    parent(data),
+                    (Nv, Nij, Nij, Nf, Topologies.nsendelems(topology)),
+                ),
+            )
+        recv_data =
+            DataLayouts.VIJFH{S, Nv, Nij, Topologies.nrecvelems(topology)}(
+                similar(
+                    parent(data),
+                    (Nv, Nij, Nij, Nf, Topologies.nrecvelems(topology)),
+                ),
+            )
         k = stride(parent(send_data), 5)
     end
 
diff --git a/test/DataLayouts/benchmark_fill.jl b/test/DataLayouts/benchmark_fill.jl
index b07f68486f..cd4576d29e 100644
--- a/test/DataLayouts/benchmark_fill.jl
+++ b/test/DataLayouts/benchmark_fill.jl
@@ -30,13 +30,13 @@ end
     Nk = 6
 #! format: off
     data = DataF{S}(device_zeros(FT,Nf));                        benchmarkfill!(device, data, 3); @test all(parent(data) .== 3)
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));         benchmarkfill!(device, data, 3); @test all(parent(data) .== 3)
-    data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));              benchmarkfill!(device, data, 3); @test all(parent(data) .== 3)
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));     benchmarkfill!(device, data, 3); @test all(parent(data) .== 3)
+    data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));          benchmarkfill!(device, data, 3); @test all(parent(data) .== 3)
     data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));             benchmarkfill!(device, data, 3); @test all(parent(data) .== 3)
     data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                  benchmarkfill!(device, data, 3); @test all(parent(data) .== 3)
     data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                    benchmarkfill!(device, data, 3); @test all(parent(data) .== 3)
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh)); benchmarkfill!(device, data, 3); @test all(parent(data) .== 3)
-    data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));      benchmarkfill!(device, data, 3); @test all(parent(data) .== 3)
+    data = VIJFH{S,Nv,Nij,Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));benchmarkfill!(device, data, 3); @test all(parent(data) .== 3)
+    data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));  benchmarkfill!(device, data, 3); @test all(parent(data) .== 3)
 #! format: on
 
     # data = IJKFVH{S}(device_zeros(FT,Nij,Nij,Nk,Nf,Nh)); benchmarkfill!(device, data, 3); @test all(parent(data) .== 3) # TODO: test
diff --git a/test/DataLayouts/cuda.jl b/test/DataLayouts/cuda.jl
index 53bd98247c..7d97158bb7 100644
--- a/test/DataLayouts/cuda.jl
+++ b/test/DataLayouts/cuda.jl
@@ -30,8 +30,9 @@ end
     S = Tuple{Complex{Float64}, Float64}
     device = ClimaComms.device()
     ArrayType = ClimaComms.array_type(device)
-    src = IJFH{S, 4}(ArrayType(rand(4, 4, 3, 10)))
-    dst = IJFH{S, 4}(ArrayType(zeros(4, 4, 3, 10)))
+    Nh = 10
+    src = IJFH{S, 4, Nh}(ArrayType(rand(4, 4, 3, Nh)))
+    dst = IJFH{S, 4, Nh}(ArrayType(zeros(4, 4, 3, Nh)))
 
     test_copy!(dst, src)
 
@@ -42,12 +43,13 @@ end
     FT = Float64
     S1 = NamedTuple{(:a, :b), Tuple{Complex{Float64}, Float64}}
     S2 = Float64
+    Nh = 2
     device = ClimaComms.device()
     ArrayType = ClimaComms.array_type(device)
-    data_arr1 = ArrayType(ones(FT, 2, 2, 3, 2))
-    data_arr2 = ArrayType(ones(FT, 2, 2, 1, 2))
-    data1 = IJFH{S1, 2}(data_arr1)
-    data2 = IJFH{S2, 2}(data_arr2)
+    data_arr1 = ArrayType(ones(FT, 2, 2, 3, Nh))
+    data_arr2 = ArrayType(ones(FT, 2, 2, 1, Nh))
+    data1 = IJFH{S1, 2, Nh}(data_arr1)
+    data2 = IJFH{S2, 2, Nh}(data_arr2)
 
     f1(a1, a2) = a1.a.re * a2 + a1.b
     res = f1.(data1, data2)
@@ -57,8 +59,8 @@ end
     Nv = 33
     data_arr1 = ArrayType(ones(FT, Nv, 4, 4, 3, 2))
     data_arr2 = ArrayType(ones(FT, Nv, 4, 4, 1, 2))
-    data1 = VIJFH{S1, Nv, 4}(data_arr1)
-    data2 = VIJFH{S2, Nv, 4}(data_arr2)
+    data1 = VIJFH{S1, Nv, 4, Nh}(data_arr1)
+    data2 = VIJFH{S2, Nv, 4, Nh}(data_arr2)
 
     f2(a1, a2) = a1.a.re * a2 + a1.b
     res = f2.(data1, data2)
@@ -71,15 +73,16 @@ end
 @testset "broadcasting assignment from scalar" begin
     FT = Float64
     S = Complex{FT}
+    Nh = 3
     device = ClimaComms.device()
     ArrayType = ClimaComms.array_type(device)
-    data = IJFH{S, 2}(ArrayType{FT}, 3)
+    data = IJFH{S, 2, Nh}(ArrayType{FT})
     data .= Complex(1.0, 2.0)
     @test Array(parent(data)) ==
           FT[f == 1 ? 1 : 2 for i in 1:2, j in 1:2, f in 1:2, h in 1:3]
 
     Nv = 33
-    data = VIJFH{S, Nv, 4}(ArrayType{FT}(undef, Nv, 4, 4, 2, 3))
+    data = VIJFH{S, Nv, 4, Nh}(ArrayType{FT}(undef, Nv, 4, 4, 2, Nh))
     data .= Complex(1.0, 2.0)
     @test Array(parent(data)) == FT[
         f == 1 ? 1 : 2 for v in 1:Nv, i in 1:4, j in 1:4, f in 1:2, h in 1:3
diff --git a/test/DataLayouts/data0d.jl b/test/DataLayouts/data0d.jl
index 82e066ebf2..66fbeb9110 100644
--- a/test/DataLayouts/data0d.jl
+++ b/test/DataLayouts/data0d.jl
@@ -135,16 +135,17 @@ end
 @testset "broadcasting DataF + IFH data object => IFH" begin
     FT = Float64
     S = Complex{FT}
+    Nh = 3
     data_f = DataF{S}(ones(FT, 2))
-    data_ifh = IFH{S, 2}(ones(FT, 2, 2, 3))
+    data_ifh = IFH{S, 2, Nh}(ones(FT, 2, 2, Nh))
     data_ifh2 = data_f .+ data_ifh
     @test data_ifh2 isa IFH{S}
     @test size(data_ifh2) == (2, 1, 1, 1, 3)
 
 
     FT = Float64
-    array = rand(FT, 2, 1, 3)
-    data = IFH{FT, 2}(array)
+    array = rand(FT, 2, 1, Nh)
+    data = IFH{FT, 2, Nh}(array)
     @test DataLayouts.data2array(data) == reshape(parent(data), :)
     @test parent(DataLayouts.array2data(DataLayouts.data2array(data), data)) ==
           parent(data)
@@ -170,15 +171,17 @@ end
 @testset "broadcasting DataF + IJFH data object => IJFH" begin
     FT = Float64
     S = Complex{FT}
+    Nh = 3
     data_f = DataF{S}(ones(FT, 2))
-    data_ijfh = IJFH{S, 2}(ones(2, 2, 2, 3))
+    data_ijfh = IJFH{S, 2, Nh}(ones(2, 2, 2, Nh))
     data_ijfh2 = data_f .+ data_ijfh
     @test data_ijfh2 isa IJFH{S}
-    @test size(data_ijfh2) == (2, 2, 1, 1, 3)
+    @test size(data_ijfh2) == (2, 2, 1, 1, Nh)
 
     FT = Float64
-    array = rand(FT, 2, 2, 1, 5)
-    data = IJFH{FT, 2}(array)
+    Nh = 3
+    array = rand(FT, 2, 2, 1, Nh)
+    data = IJFH{FT, 2, Nh}(array)
     @test DataLayouts.data2array(data) == reshape(parent(data), :)
     @test parent(DataLayouts.array2data(DataLayouts.data2array(data), data)) ==
           parent(data)
@@ -187,23 +190,25 @@ end
 @testset "broadcasting DataF + VIFH data object => VIFH" begin
     FT = Float64
     S = Complex{FT}
+    Nh = 10
     data_f = DataF{S}(ones(FT, 2))
     Nv = 10
-    data_vifh = VIFH{S, Nv, 4}(ones(FT, Nv, 4, 2, 10))
+    data_vifh = VIFH{S, Nv, 4, Nh}(ones(FT, Nv, 4, 2, Nh))
     data_vifh2 = data_f .+ data_vifh
     @test data_vifh2 isa VIFH{S, Nv}
-    @test size(data_vifh2) == (4, 1, 1, Nv, 10)
+    @test size(data_vifh2) == (4, 1, 1, Nv, Nh)
 end
 
 @testset "broadcasting DataF + VIJFH data object => VIJFH" begin
     FT = Float64
     S = Complex{FT}
     Nv = 2
+    Nh = 2
     data_f = DataF{S}(ones(FT, 2))
-    data_vijfh = VIJFH{S, Nv, 2}(ones(FT, Nv, 2, 2, 2, 2))
+    data_vijfh = VIJFH{S, Nv, 2, Nh}(ones(FT, Nv, 2, 2, 2, Nh))
     data_vijfh2 = data_f .+ data_vijfh
-    @test data_vijfh2 isa VIJFH{S, Nv}
-    @test size(data_vijfh2) == (2, 2, 1, Nv, 2)
+    @test data_vijfh2 isa VIJFH{S, Nv, Nh}
+    @test size(data_vijfh2) == (2, 2, 1, Nv, Nh)
 end
 
 @testset "column IF => DataF" begin
@@ -220,9 +225,10 @@ end
 @testset "column IFH => DataF" begin
     FT = Float64
     S = Complex{FT}
-    array = ones(FT, 2, 2, 3)
+    Nh = 3
+    array = ones(FT, 2, 2, Nh)
     array[1, :, 1] .= FT[3, 4]
-    data_ifh = IFH{S, 2}(array)
+    data_ifh = IFH{S, 2, Nh}(array)
     ifh_column = column(data_ifh, 1, 1)
     @test ifh_column isa DataF
     @test ifh_column[] == 3.0 + 4.0im
@@ -246,9 +252,10 @@ end
 @testset "column IJFH => DataF" begin
     FT = Float64
     S = Complex{FT}
+    Nh = 3
     array = ones(2, 2, 2, 3)
     array[1, 1, :, 2] .= FT[3, 4]
-    data_ijfh = IJFH{S, 2}(array)
+    data_ijfh = IJFH{S, 2, Nh}(array)
     ijfh_column = column(data_ijfh, 1, 1, 2)
     @test ijfh_column isa DataF
     @test ijfh_column[] == 3.0 + 4.0im
diff --git a/test/DataLayouts/data1dx.jl b/test/DataLayouts/data1dx.jl
index 9163e7d3b1..5812bfbe41 100644
--- a/test/DataLayouts/data1dx.jl
+++ b/test/DataLayouts/data1dx.jl
@@ -18,7 +18,7 @@ import ClimaCore.DataLayouts: VIFH, slab, column, VF, IFH
         # 10 elements in horizontal with 4 nodal points per element in horizontal
         array = rand(FT, Nv, Ni, 3, Nh)
 
-        data = VIFH{S, Nv, Ni}(array)
+        data = VIFH{S, Nv, Ni, Nh}(array)
         sum(x -> x[2], data)
 
         @test getfield(data.:1, :array) == @view(array[:, :, 1:2, :])
@@ -48,7 +48,7 @@ import ClimaCore.DataLayouts: VIFH, slab, column, VF, IFH
     Ni = 4  # number of nodal points
     Nh = 10 # number of elements
     array = rand(FT, Nv, Ni, 1, Nh)
-    data = VIFH{FT, Nv, Ni}(array)
+    data = VIFH{FT, Nv, Ni, Nh}(array)
     @test DataLayouts.data2array(data) ==
           reshape(parent(data), DataLayouts.nlevels(data), :)
     @test parent(DataLayouts.array2data(DataLayouts.data2array(data), data)) ==
@@ -62,7 +62,7 @@ end
 
     S = Tuple{Complex{Float64}, Float64}
     array = zeros(Float64, Nv, Ni, 3, Nh)
-    data = VIFH{S, Nv, Ni}(array)
+    data = VIFH{S, Nv, Ni, Nh}(array)
 
     @test_throws BoundsError slab(data, -1, -1)
     @test_throws BoundsError slab(data, 1, 3)
@@ -88,7 +88,7 @@ end
     SB = (c = 1.0, d = 2.0)
 
     array = zeros(Float64, Nv, Ni, 2, Nh)
-    data = VIFH{typeof(SA), Nv, Ni}(array)
+    data = VIFH{typeof(SA), Nv, Ni, Nh}(array)
 
     cdata = column(data, 1, 1)
     cdata[1] = SA
@@ -103,9 +103,10 @@ end
 @testset "broadcasting between VIFH data object + scalars" begin
     FT = Float64
     Nv = 2
+    Nh = 2
     data1 = ones(FT, Nv, 2, 2, 2)
     S = Complex{Float64}
-    data1 = VIFH{S, Nv, 2}(data1)
+    data1 = VIFH{S, Nv, 2, Nh}(data1)
     res = data1 .+ 1
     @test res isa VIFH{S, Nv}
     @test parent(res) ==
@@ -118,8 +119,9 @@ end
     FT = Float64
     S = Complex{FT}
     Nv = 3
+    Nh = 2
     data_vf = VF{S, Nv}(ones(FT, Nv, 2))
-    data_ifh = IFH{FT, 2}(ones(FT, 2, 1, 2))
+    data_ifh = IFH{FT, 2, Nh}(ones(FT, 2, 1, 2))
     data_vifh = data_vf .+ data_ifh
     @test data_vifh isa VIFH{S, Nv}
     @test size(data_vifh) == (2, 1, 1, 3, 2)
@@ -135,7 +137,8 @@ end
 
 @testset "fill" begin
 
-    data = IFH{Float64, 3}(ones(3, 1, 3))
+    Nh = 3
+    data = IFH{Float64, 3, Nh}(ones(3, 1, Nh))
     data .= 2.0
     @test all(==(2.0), parent(data))
 end
diff --git a/test/DataLayouts/data2d.jl b/test/DataLayouts/data2d.jl
index 2cf518a0d3..e459164045 100644
--- a/test/DataLayouts/data2d.jl
+++ b/test/DataLayouts/data2d.jl
@@ -46,7 +46,7 @@ end
     Nh = 2 # number of elements
     S = Tuple{Complex{Float64}, Float64}
     array = rand(Nij, Nij, 3, Nh)
-    data = IJFH{S, 2}(array)
+    data = IJFH{S, 2, Nh}(array)
     @test getfield(data.:1, :array) == @view(array[:, :, 1:2, :])
     data_slab = slab(data, 1)
     @test data_slab[2, 1] ==
@@ -69,7 +69,7 @@ end
     Nij = 4  # number of nodal points
     Nh = 10 # number of elements
     array = rand(FT, Nij, Nij, 1, Nh)
-    data = IJFH{FT, Nij}(array)
+    data = IJFH{FT, Nij, Nh}(array)
     @test DataLayouts.data2array(data) == reshape(parent(data), :)
     @test parent(DataLayouts.array2data(DataLayouts.data2array(data), data)) ==
           parent(data)
@@ -79,8 +79,8 @@ end
     Nij = 1 # number of nodal points
     Nh = 2 # number of elements
     S = Tuple{Complex{Float64}, Float64}
-    array = zeros(Float64, Nij, Nij, 3, 2)
-    data = IJFH{S, Nij}(array)
+    array = zeros(Float64, Nij, Nij, 3, Nh)
+    data = IJFH{S, Nij, Nh}(array)
 
     @test_throws BoundsError slab(data, -1)
     @test_throws BoundsError slab(data, 3)
@@ -104,7 +104,7 @@ end
     SB = (c = 1.0, d = 2.0)
 
     array = zeros(Float64, Nij, Nij, 2, Nh)
-    data = IJFH{typeof(SA), Nij}(array)
+    data = IJFH{typeof(SA), Nij, Nh}(array)
     data_slab = slab(data, 1)
     ret = begin
         data_slab[1, 1] = SA
@@ -120,10 +120,10 @@ end
     S1 = Float64
     S2 = Float32
     array1 = ones(S1, Nij, Nij, 1, Nh)
-    data1 = IJFH{S1, Nij}(array1)
+    data1 = IJFH{S1, Nij, Nh}(array1)
 
     array2 = ones(S2, Nij, Nij, 1, Nh)
-    data2 = IJFH{S2, Nij}(array2)
+    data2 = IJFH{S2, Nij, Nh}(array2)
 
     for h in 1:Nh
         slab1 = slab(data1, h)
@@ -137,9 +137,10 @@ end
 
 @testset "broadcasting between data object + scalars" begin
     FT = Float64
-    data1 = ones(FT, 2, 2, 2, 2)
+    Nh = 2
+    data1 = ones(FT, 2, 2, 2, Nh)
     S = Complex{Float64}
-    data1 = IJFH{S, 2}(data1)
+    data1 = IJFH{S, 2, Nh}(data1)
     res = data1 .+ 1
     @test res isa IJFH{S}
     @test parent(res) ==
@@ -152,7 +153,8 @@ end
 @testset "broadcasting assignment from scalar" begin
     FT = Float64
     S = Complex{FT}
-    data = IJFH{S, 2}(Array{FT}, 3)
+    Nh = 3
+    data = IJFH{S, 2, Nh}(Array{FT})
     data .= Complex(1.0, 2.0)
     @test parent(data) ==
           FT[f == 1 ? 1 : 2 for i in 1:2, j in 1:2, f in 1:2, h in 1:3]
@@ -165,12 +167,13 @@ end
 
 @testset "broadcasting between data objects" begin
     FT = Float64
-    data1 = ones(FT, 2, 2, 2, 2)
-    data2 = ones(FT, 2, 2, 1, 2)
+    Nh = 2
+    data1 = ones(FT, 2, 2, 2, Nh)
+    data2 = ones(FT, 2, 2, 1, Nh)
     S1 = Complex{Float64}
     S2 = Float64
-    data1 = IJFH{S1, 2}(data1)
-    data2 = IJFH{S2, 2}(data2)
+    data1 = IJFH{S1, 2, Nh}(data1)
+    data2 = IJFH{S2, 2, Nh}(data2)
     res = data1 .+ data2
     @test res isa IJFH{S1}
     @test parent(res) ==
@@ -183,11 +186,12 @@ end
 @testset "broadcasting complicated function" begin
     FT = Float64
     S1 = NamedTuple{(:a, :b), Tuple{Complex{Float64}, Float64}}
-    data1 = ones(FT, 2, 2, 3, 2)
+    Nh = 2
+    data1 = ones(FT, 2, 2, 3, Nh)
     S2 = Float64
-    data2 = ones(FT, 2, 2, 1, 2)
-    data1 = IJFH{S1, 2}(data1)
-    data2 = IJFH{S2, 2}(data2)
+    data2 = ones(FT, 2, 2, 1, Nh)
+    data1 = IJFH{S1, 2, Nh}(data1)
+    data2 = IJFH{S2, 2, Nh}(data2)
 
     f(a1, a2) = a1.a.re * a2 + a1.b
     res = f.(data1, data2)
diff --git a/test/DataLayouts/data2dx.jl b/test/DataLayouts/data2dx.jl
index 304c24d0ae..3abbfd5291 100644
--- a/test/DataLayouts/data2dx.jl
+++ b/test/DataLayouts/data2dx.jl
@@ -19,7 +19,7 @@ import ClimaCore.DataLayouts: VF, IJFH, VIJFH, slab, column
         # 10 elements in horizontal with 4 × 4 nodal points per element in horizontal
         array = rand(FT, Nv, Nij, Nij, 3, Nh)
 
-        data = VIJFH{S, Nv, Nij}(array)
+        data = VIJFH{S, Nv, Nij, Nh}(array)
 
         @test getfield(data.:1, :array) == @view(array[:, :, :, 1:2, :])
         @test getfield(data.:2, :array) == @view(array[:, :, :, 3:3, :])
@@ -49,7 +49,7 @@ import ClimaCore.DataLayouts: VF, IJFH, VIJFH, slab, column
     Ni = 4  # number of nodal points
     Nh = 10 # number of elements
     array = rand(FT, Nv, Nij, Nij, 1, Nh)
-    data = VIJFH{FT, Nv, Nij}(array)
+    data = VIJFH{FT, Nv, Nij, Nh}(array)
     @test DataLayouts.data2array(data) ==
           reshape(parent(data), DataLayouts.nlevels(data), :)
     @test parent(DataLayouts.array2data(DataLayouts.data2array(data), data)) ==
@@ -63,7 +63,7 @@ end
 
     S = Tuple{Complex{Float64}, Float64}
     array = zeros(Float64, Nv, Nij, Nij, 3, Nh)
-    data = VIJFH{S, Nv, Nij}(array)
+    data = VIJFH{S, Nv, Nij, Nh}(array)
 
     @test_throws BoundsError slab(data, -1, 1)
     @test_throws BoundsError slab(data, 1, -1)
@@ -88,7 +88,7 @@ end
     SB = (c = 1.0, d = 2.0)
 
     array = zeros(Float64, Nv, Nij, Nij, 2, Nh)
-    data = VIJFH{typeof(SA), Nv, Nij}(array)
+    data = VIJFH{typeof(SA), Nv, Nij, Nh}(array)
 
     cdata = column(data, 1, 2, 1)
     cdata[1] = SA
@@ -104,8 +104,9 @@ end
     FT = Float64
     array = ones(FT, 2, 2, 2, 2, 2)
     Nv = size(array, 1)
+    Nh = size(array, 5)
     S = Complex{Float64}
-    data1 = VIJFH{S, Nv, 2}(array)
+    data1 = VIJFH{S, Nv, 2, Nh}(array)
     res = data1 .+ 1
     @test res isa VIJFH{S, Nv}
     @test parent(res) == FT[
@@ -120,8 +121,9 @@ end
     FT = Float64
     S = Complex{FT}
     Nv = 3
+    Nh = 2
     data_vf = VF{S, Nv}(ones(FT, Nv, 2))
-    data_ijfh = IJFH{FT, 2}(ones(FT, 2, 2, 1, 2))
+    data_ijfh = IJFH{FT, 2, Nh}(ones(FT, 2, 2, 1, Nh))
     data_vijfh = data_vf .+ data_ijfh
     @test data_vijfh isa VIJFH{S, Nv}
     @test size(data_vijfh) == (2, 2, 1, 3, 2)
diff --git a/test/DataLayouts/opt_similar.jl b/test/DataLayouts/opt_similar.jl
index f69383a34c..807b496773 100644
--- a/test/DataLayouts/opt_similar.jl
+++ b/test/DataLayouts/opt_similar.jl
@@ -38,14 +38,14 @@ end
     Nk = 6
 #! format: off
     data = DataF{S}(device_zeros(FT,Nf));                        test_similar!(data)
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));         test_similar!(data)
-    data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));              test_similar!(data)
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));     test_similar!(data)
+    data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));          test_similar!(data)
     data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));             test_similar!(data)
     data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                  test_similar!(data)
     data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                    test_similar!(data)
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh)); test_similar!(data)
-    data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_similar!(data)
+    data = VIJFH{S,Nv,Nij,Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));test_similar!(data)
+    data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));  test_similar!(data)
 #! format: on
-    # data = DataLayouts.IJKFVH{S, Nij, Nk}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_similar!(data) # TODO: test
-    # data = DataLayouts.IH1JH2{S, Nij}(device_zeros(FT,2*Nij,3*Nij));             test_similar!(data) # TODO: test
+    # data = DataLayouts.IJKFVH{S, Nij, Nk, Nv, Nh}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_similar!(data) # TODO: test
+    # data = DataLayouts.IH1JH2{S, Nij}(device_zeros(FT,2*Nij,3*Nij));                     test_similar!(data) # TODO: test
 end
diff --git a/test/DataLayouts/unit_copyto.jl b/test/DataLayouts/unit_copyto.jl
index 90a5c53d0b..0b304a4f81 100644
--- a/test/DataLayouts/unit_copyto.jl
+++ b/test/DataLayouts/unit_copyto.jl
@@ -53,13 +53,13 @@ end
     Nk = 6
 #! format: off
     data = DataF{S}(device_zeros(FT,Nf));                        test_copyto_float!(data)
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));         test_copyto_float!(data)
-    data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));              test_copyto_float!(data)
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));     test_copyto_float!(data)
+    data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));          test_copyto_float!(data)
     data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));             test_copyto_float!(data)
     data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                  test_copyto_float!(data)
     data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                    test_copyto_float!(data)
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh)); test_copyto_float!(data)
-    data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_copyto_float!(data)
+    data = VIJFH{S,Nv,Nij,Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));test_copyto_float!(data)
+    data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));  test_copyto_float!(data)
 #! format: on
     # data = DataLayouts.IJKFVH{S, Nij, Nk}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_copyto_float!(data) # TODO: test
     # data = DataLayouts.IH1JH2{S, Nij}(device_zeros(FT,2*Nij,3*Nij));             test_copyto_float!(data) # TODO: test
@@ -77,13 +77,13 @@ end
     Nk = 6
 #! format: off
     data = DataF{S}(device_zeros(FT,Nf));                        test_copyto!(data)
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));         test_copyto!(data)
-    data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));              test_copyto!(data)
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));     test_copyto!(data)
+    data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));          test_copyto!(data)
     data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));             test_copyto!(data)
     data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                  test_copyto!(data)
     data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                    test_copyto!(data)
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh)); test_copyto!(data)
-    data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_copyto!(data)
+    data = VIJFH{S,Nv,Nij,Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));test_copyto!(data)
+    data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));  test_copyto!(data)
 #! format: on
     # TODO: test this
     # data = DataLayouts.IJKFVH{S, Nij, Nk}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_copyto!(data) # TODO: test
@@ -113,13 +113,13 @@ end
     # Rather than using level/slab/column, let's just make views/SubArrays
     # directly so that we can easily test all cases:
 #! format: off
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));         test_copyto!(data_view(data))
-    data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));              test_copyto!(data_view(data))
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));     test_copyto!(data_view(data))
+    data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));          test_copyto!(data_view(data))
     data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));             test_copyto!(data_view(data))
     data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                  test_copyto!(data_view(data))
     data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                    test_copyto!(data_view(data))
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh)); test_copyto!(data_view(data))
-    data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_copyto!(data_view(data))
+    data = VIJFH{S,Nv,Nij,Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));test_copyto!(data_view(data))
+    data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));  test_copyto!(data_view(data))
 #! format: on
     # TODO: test this
     # data = DataLayouts.IJKFVH{S, Nij, Nk}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_copyto!(data) # TODO: test
diff --git a/test/DataLayouts/unit_fill.jl b/test/DataLayouts/unit_fill.jl
index c9a0c63141..fc803c0015 100644
--- a/test/DataLayouts/unit_fill.jl
+++ b/test/DataLayouts/unit_fill.jl
@@ -29,13 +29,13 @@ end
     Nk = 6
 #! format: off
     data = DataF{S}(device_zeros(FT,Nf));                        test_fill!(data, 3)
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));         test_fill!(data, 3)
-    data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));              test_fill!(data, 3)
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));     test_fill!(data, 3)
+    data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));          test_fill!(data, 3)
     data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));             test_fill!(data, 3)
     data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                  test_fill!(data, 3)
     data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                    test_fill!(data, 3)
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh)); test_fill!(data, 3)
-    data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_fill!(data, 3)
+    data = VIJFH{S,Nv,Nij,Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));test_fill!(data, 3)
+    data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));  test_fill!(data, 3)
 #! format: on
     # data = DataLayouts.IJKFVH{S, Nij, Nk}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_fill!(data, 3) # TODO: test
     # data = DataLayouts.IH1JH2{S, Nij}(device_zeros(FT,2*Nij,3*Nij));             test_fill!(data, 3) # TODO: test
@@ -53,13 +53,13 @@ end
     Nk = 6
 #! format: off
     data = DataF{S}(device_zeros(FT,Nf));                        test_fill!(data, (2,3))
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));         test_fill!(data, (2,3))
-    data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));              test_fill!(data, (2,3))
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));     test_fill!(data, (2,3))
+    data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));          test_fill!(data, (2,3))
     data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));             test_fill!(data, (2,3))
     data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                  test_fill!(data, (2,3))
     data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                    test_fill!(data, (2,3))
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh)); test_fill!(data, (2,3))
-    data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_fill!(data, (2,3))
+    data = VIJFH{S,Nv,Nij,Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));test_fill!(data, (2,3))
+    data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));  test_fill!(data, (2,3))
 #! format: on
     # TODO: test this
     # data = DataLayouts.IJKFVH{S, Nij, Nk}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_fill!(data, (2,3)) # TODO: test
@@ -86,13 +86,13 @@ end
     # Rather than using level/slab/column, let's just make views/SubArrays
     # directly so that we can easily test all cases:
 #! format: off
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));         test_fill!(data_view(data), (2,3))
-    data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));              test_fill!(data_view(data), (2,3))
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));     test_fill!(data_view(data), (2,3))
+    data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));          test_fill!(data_view(data), (2,3))
     data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));             test_fill!(data_view(data), (2,3))
     data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                  test_fill!(data_view(data), (2,3))
     data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                    test_fill!(data_view(data), (2,3))
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh)); test_fill!(data_view(data), (2,3))
-    data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_fill!(data_view(data), (2,3))
+    data = VIJFH{S,Nv,Nij,Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));test_fill!(data_view(data), (2,3))
+    data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));  test_fill!(data_view(data), (2,3))
 #! format: on
     # TODO: test this
     # data = DataLayouts.IJKFVH{S, Nij, Nk}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_fill!(data, (2,3)) # TODO: test
@@ -138,13 +138,13 @@ end
     Nk = 6
     # directly so that we can easily test all cases:
 #! format: off
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));         test_fill!(reshaped_array(data), 2)
-    data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));              test_fill!(reshaped_array(data), 2)
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));     test_fill!(reshaped_array(data), 2)
+    data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));          test_fill!(reshaped_array(data), 2)
     # data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));             test_fill!(reshaped_array(data), 2)
     # data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                  test_fill!(reshaped_array(data), 2)
     # data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                    test_fill!(reshaped_array(data), 2)
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh)); test_fill!(reshaped_array(data), 2)
-    data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_fill!(reshaped_array(data), 2)
+    data = VIJFH{S,Nv,Nij,Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));test_fill!(reshaped_array(data), 2)
+    data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));  test_fill!(reshaped_array(data), 2)
 #! format: on
     # TODO: test this
     # data = DataLayouts.IJKFVH{S, Nij, Nk}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_fill!(reshaped_array(data), 2) # TODO: test
diff --git a/test/DataLayouts/unit_mapreduce.jl b/test/DataLayouts/unit_mapreduce.jl
index 6b4512b23c..2da4547521 100644
--- a/test/DataLayouts/unit_mapreduce.jl
+++ b/test/DataLayouts/unit_mapreduce.jl
@@ -73,17 +73,17 @@ end
     Nh = 5
     Nk = 6
 #! format: off
-    data = DataF{S}(device_zeros(FT,Nf));                        test_mapreduce_1!(context, data)
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));         test_mapreduce_1!(context, data)
-    # data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));              test_mapreduce_1!(context, data)
-    # data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));             test_mapreduce_1!(context, data)
-    # data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                  test_mapreduce_1!(context, data)
-    data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                    test_mapreduce_1!(context, data)
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh)); test_mapreduce_1!(context, data)
-    # data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_mapreduce_1!(context, data)
+    data = DataF{S}(device_zeros(FT,Nf));                              test_mapreduce_1!(context, data)
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));           test_mapreduce_1!(context, data)
+    # data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));              test_mapreduce_1!(context, data)
+    # data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));                 test_mapreduce_1!(context, data)
+    # data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                      test_mapreduce_1!(context, data)
+    data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                          test_mapreduce_1!(context, data)
+    data = VIJFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));   test_mapreduce_1!(context, data)
+    # data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_mapreduce_1!(context, data)
 #! format: on
-    # data = DataLayouts.IJKFVH{S, Nij, Nk}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_mapreduce_1!(context, data) # TODO: test
-    # data = DataLayouts.IH1JH2{S, Nij}(device_zeros(FT,2*Nij,3*Nij));             test_mapreduce_1!(context, data) # TODO: test
+    # data = DataLayouts.IJKFVH{S, Nij, Nk, Nh}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_mapreduce_1!(context, data) # TODO: test
+    # data = DataLayouts.IH1JH2{S, Nij}(device_zeros(FT,2*Nij,3*Nij));                 test_mapreduce_1!(context, data) # TODO: test
 end
 
 @testset "mapreduce with Nf > 1" begin
@@ -96,18 +96,18 @@ end
     Nh = 5
     Nk = 6
 #! format: off
-    data = DataF{S}(device_zeros(FT,Nf));                        test_mapreduce_2!(context, data)
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));         test_mapreduce_2!(context, data)
-    # data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));              test_mapreduce_2!(context, data)
-    # data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));             test_mapreduce_2!(context, data)
-    # data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                  test_mapreduce_2!(context, data)
-    data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                    test_mapreduce_2!(context, data)
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh)); test_mapreduce_2!(context, data)
-    # data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_mapreduce_2!(context, data)
+    data = DataF{S}(device_zeros(FT,Nf));                              test_mapreduce_2!(context, data)
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));           test_mapreduce_2!(context, data)
+    # data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));              test_mapreduce_2!(context, data)
+    # data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));                 test_mapreduce_2!(context, data)
+    # data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                      test_mapreduce_2!(context, data)
+    data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                          test_mapreduce_2!(context, data)
+    data = VIJFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));   test_mapreduce_2!(context, data)
+    # data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_mapreduce_2!(context, data)
 #! format: on
     # TODO: test this
-    # data = DataLayouts.IJKFVH{S, Nij, Nk}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_mapreduce_2!(context, data) # TODO: test
-    # data = DataLayouts.IH1JH2{S, Nij}(device_zeros(FT,2*Nij,3*Nij));             test_mapreduce_2!(context, data) # TODO: test
+    # data = DataLayouts.IJKFVH{S, Nij, Nk, Nh}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_mapreduce_2!(context, data) # TODO: test
+    # data = DataLayouts.IH1JH2{S, Nij}(device_zeros(FT,2*Nij,3*Nij));                 test_mapreduce_2!(context, data) # TODO: test
 end
 
 @testset "mapreduce views with Nf > 1" begin
@@ -129,16 +129,16 @@ end
     # Rather than using level/slab/column, let's just make views/SubArrays
     # directly so that we can easily test all cases:
 #! format: off
-    data = DataF{S}(device_zeros(FT,Nf));                        test_mapreduce_2!(context, data_view(data))
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));         test_mapreduce_2!(context, data_view(data))
-    # data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));              test_mapreduce_2!(context, data_view(data))
-    # data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));             test_mapreduce_2!(context, data_view(data))
-    # data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                  test_mapreduce_2!(context, data_view(data))
-    data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                    test_mapreduce_2!(context, data_view(data))
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh)); test_mapreduce_2!(context, data_view(data))
-    # data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_mapreduce_2!(context, data_view(data))
+    data = DataF{S}(device_zeros(FT,Nf));                              test_mapreduce_2!(context, data_view(data))
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));           test_mapreduce_2!(context, data_view(data))
+    # data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));              test_mapreduce_2!(context, data_view(data))
+    # data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));                 test_mapreduce_2!(context, data_view(data))
+    # data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                      test_mapreduce_2!(context, data_view(data))
+    data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                          test_mapreduce_2!(context, data_view(data))
+    data = VIJFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));   test_mapreduce_2!(context, data_view(data))
+    # data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));      test_mapreduce_2!(context, data_view(data))
 #! format: on
     # TODO: test this
-    # data = DataLayouts.IJKFVH{S, Nij, Nk}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_mapreduce_2!(context, data_view(data)) # TODO: test
-    # data = DataLayouts.IH1JH2{S, Nij}(device_zeros(FT,2*Nij,3*Nij));             test_mapreduce_2!(context, data_view(data)) # TODO: test
+    # data = DataLayouts.IJKFVH{S, Nij, Nk, Nh}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); test_mapreduce_2!(context, data_view(data)) # TODO: test
+    # data = DataLayouts.IH1JH2{S, Nij}(device_zeros(FT,2*Nij,3*Nij));                 test_mapreduce_2!(context, data_view(data)) # TODO: test
 end
diff --git a/test/DataLayouts/unit_ndims.jl b/test/DataLayouts/unit_ndims.jl
index 3a774122f4..75c62e7755 100644
--- a/test/DataLayouts/unit_ndims.jl
+++ b/test/DataLayouts/unit_ndims.jl
@@ -18,15 +18,15 @@ ClimaComms.@import_required_backends
     Nh = 5
     Nk = 6
 #! format: off
-    data = DataF{S}(device_zeros(FT,Nf));                                        @test ndims(data) == 1; @test ndims(typeof(data)) == 1
-    data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                                  @test ndims(data) == 2; @test ndims(typeof(data)) == 2
-    data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                                    @test ndims(data) == 2; @test ndims(typeof(data)) == 2
-    data = IFH{S, Nij}(device_zeros(FT,Nij,Nf,Nh));                              @test ndims(data) == 3; @test ndims(typeof(data)) == 3
-    data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));                             @test ndims(data) == 3; @test ndims(typeof(data)) == 3
-    data = IJFH{S, Nij}(device_zeros(FT,Nij,Nij,Nf,Nh));                         @test ndims(data) == 4; @test ndims(typeof(data)) == 4
-    data = VIFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nf,Nh));                      @test ndims(data) == 4; @test ndims(typeof(data)) == 4
-    data = VIJFH{S, Nv, Nij}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));                 @test ndims(data) == 5; @test ndims(typeof(data)) == 5
-    data = DataLayouts.IJKFVH{S, Nij, Nk}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); @test ndims(data) == 6; @test ndims(typeof(data)) == 6
-    data = DataLayouts.IH1JH2{S, Nij}(device_zeros(FT,2*Nij,3*Nij));             @test ndims(data) == 2; @test ndims(typeof(data)) == 2
+    data = DataF{S}(device_zeros(FT,Nf));                                            @test ndims(data) == 1; @test ndims(typeof(data)) == 1
+    data = IF{S, Nij}(device_zeros(FT,Nij,Nf));                                      @test ndims(data) == 2; @test ndims(typeof(data)) == 2
+    data = VF{S, Nv}(device_zeros(FT,Nv,Nf));                                        @test ndims(data) == 2; @test ndims(typeof(data)) == 2
+    data = IFH{S, Nij, Nh}(device_zeros(FT,Nij,Nf,Nh));                              @test ndims(data) == 3; @test ndims(typeof(data)) == 3
+    data = IJF{S, Nij}(device_zeros(FT,Nij,Nij,Nf));                                 @test ndims(data) == 3; @test ndims(typeof(data)) == 3
+    data = IJFH{S, Nij, Nh}(device_zeros(FT,Nij,Nij,Nf,Nh));                         @test ndims(data) == 4; @test ndims(typeof(data)) == 4
+    data = VIFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nf,Nh));                      @test ndims(data) == 4; @test ndims(typeof(data)) == 4
+    data = VIJFH{S, Nv, Nij, Nh}(device_zeros(FT,Nv,Nij,Nij,Nf,Nh));                 @test ndims(data) == 5; @test ndims(typeof(data)) == 5
+    data = DataLayouts.IJKFVH{S,Nij,Nk,Nv,Nh}(device_zeros(FT,Nij,Nij,Nk,Nf,Nv,Nh)); @test ndims(data) == 6; @test ndims(typeof(data)) == 6
+    data = DataLayouts.IH1JH2{S, Nij}(device_zeros(FT,2*Nij,3*Nij));                 @test ndims(data) == 2; @test ndims(typeof(data)) == 2
 #! format: on
 end
diff --git a/test/Fields/unit_field.jl b/test/Fields/unit_field.jl
index 475dd583ef..37f987701f 100644
--- a/test/Fields/unit_field.jl
+++ b/test/Fields/unit_field.jl
@@ -55,10 +55,13 @@ end
 @testset "1×1 2D domain space" begin
     Nij = 4
     n1 = n2 = 1
+    Nh = n1 * n2
     space = spectral_space_2D(n1 = n1, n2 = n2, Nij = Nij)
 
-    field =
-        Fields.Field(IJFH{ComplexF64, Nij}(ones(Nij, Nij, 2, n1 * n2)), space)
+    field = Fields.Field(
+        IJFH{ComplexF64, Nij, Nh}(ones(Nij, Nij, 2, n1 * n2)),
+        space,
+    )
 
     @test sum(field) ≈ Complex(1.0, 1.0) * 8.0 * 10.0 rtol = 10eps()
     @test sum(x -> 3.0, field) ≈ 3 * 8.0 * 10.0 rtol = 10eps()
@@ -245,11 +248,12 @@ end
 @testset "Broadcasting interception for tuple-valued fields" begin
     n1 = n2 = 1
     Nij = 4
+    Nh = n1 * n2
     space = spectral_space_2D(n1 = n1, n2 = n2, Nij = Nij)
 
     nt_field = Fields.Field(
-        IJFH{NamedTuple{(:a, :b), Tuple{Float64, Float64}}, Nij}(
-            ones(Nij, Nij, 2, n1 * n2),
+        IJFH{NamedTuple{(:a, :b), Tuple{Float64, Float64}}, Nij, Nh}(
+            ones(Nij, Nij, 2, Nh),
         ),
         space,
     )
diff --git a/test/Spaces/opt_spaces.jl b/test/Spaces/opt_spaces.jl
index 7eb3d4541a..c21bde5cd9 100644
--- a/test/Spaces/opt_spaces.jl
+++ b/test/Spaces/opt_spaces.jl
@@ -60,7 +60,7 @@ end
 
         result = JET.@report_opt Grids._SpectralElementGrid2D(Spaces.topology(space), Spaces.quadrature_style(space); enable_bubble=false)
         n_found = length(JET.get_reports(result.analyzer, result.result))
-        n_allowed = 12
+        n_allowed = 177
         @test n_found ≤ n_allowed
         if n_found < n_allowed
             @info "Inference may have improved for _SpectralElementGrid2D: (n_found, n_allowed) = ($n_found, $n_allowed)"
diff --git a/test/Spaces/unit_spaces.jl b/test/Spaces/unit_spaces.jl
index dee58809ee..68a8884450 100644
--- a/test/Spaces/unit_spaces.jl
+++ b/test/Spaces/unit_spaces.jl
@@ -331,7 +331,7 @@ end
     @test size(array) == (4, 4, 2, 4)
 
     Nij = length(points)
-    field = Fields.Field(IJFH{FT, Nij}(ones(Nij, Nij, 1, n1 * n2)), space)
+    field = Fields.Field(IJFH{FT, Nij, n1 * n2}(ones(Nij, Nij, 1, n1 * n2)), space)
     field_values = Fields.field_values(field)
     Spaces.horizontal_dss!(field)
 
@@ -414,7 +414,7 @@ end
     data[:, :, 1, :] .= 1:Nij
     data[:, :, 2, :] .= (1:Nij)'
     data[:, :, 3, :] .= reshape(1:(n1 * n2), 1, 1, :)
-    field = Fields.Field(IJFH{Tuple{FT, FT, FT}, Nij}(data), space)
+    field = Fields.Field(IJFH{Tuple{FT, FT, FT}, Nij, n1 * n2}(data), space)
     field_dss = Spaces.horizontal_dss!(copy(field))
     data_dss = parent(field_dss)