diff --git a/NEWS.md b/NEWS.md
index 5da18861b3..13d353d291 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,29 @@ ClimaCore.jl Release Notes
 main
 -------
 
+- ![][badge-🔥behavioralΔ] Unified strong/weak spectral element operator variants
+  via a `FormType` parameter. `Divergence`, `Gradient`, and `Curl` now carry a
+  second type parameter (`StrongForm` or `WeakForm`), and `WeakDivergence`,
+  `WeakGradient`, and `WeakCurl` are backwards-compatibility aliases for the
+  `WeakForm` variants (e.g. `WeakDivergence{I} = Divergence{I, WeakForm}`),
+  collected in `src/Operators/deprecated.jl` for eventual removal. All operator
+  names and their public constructors are unchanged. Results are bitwise
+  identical to the previous implementations, with one exception: the strong-form
+  `Divergence` on GPUs in 2D now keeps one accumulator per dimension (matching
+  the other unified operators) instead of summing both dimensions into a single
+  accumulator, which changes results by accumulation roundoff only (verified
+  identical in exact arithmetic; float differences are a few `eps` of the
+  accumulated terms). Two notes for downstream code:
+  `isa` checks against the strong operator types now also match the weak
+  variants (e.g. `WeakDivergence() isa Divergence` is `true`), so dispatch that
+  must distinguish forms should use the `FormType` parameter (e.g.
+  `Divergence{I, StrongForm}`); and the internal
+  `Divergence{()}(space)`-style constructors used to reset operator axes have
+  been replaced by `rebuild_operator(op, space)`. As a side effect of unifying
+  the CUDA `Gradient` kernels, `WeakGradient` now accepts rank-1 tensor fields
+  on GPUs, which previously raised an "Unsupported input type" error even though
+  both forms already supported them on CPUs.
+  [2556](https://github.com/CliMA/ClimaCore.jl/pull/2556)
 - Refactor DataLayouts module [2522](https://github.com/CliMA/ClimaCore.jl/pull/2522)
   - All data layout types are unified into a single `DataLayout` type, and all
     loops over data go through two communication primitives (`foreach_slice` and
diff --git a/docs/src/operators.md b/docs/src/operators.md
index b3b6a25d8b..022411b113 100644
--- a/docs/src/operators.md
+++ b/docs/src/operators.md
@@ -28,6 +28,18 @@ Curl
 WeakCurl
 ```
 
+### Strong and weak forms
+
+`Divergence`, `Gradient`, and `Curl` each have a strong and a weak variant,
+distinguished by a [`FormType`](@ref) type parameter. The weak variants are also
+available under the names `WeakDivergence`, `WeakGradient`, and `WeakCurl`.
+
+```@docs
+FormType
+StrongForm
+WeakForm
+```
+
 ### Interpolation Operators
 ```@docs
 Interpolate
diff --git a/ext/cuda/operators_sem_shmem.jl b/ext/cuda/operators_sem_shmem.jl
index 423a8b18a1..5a448716ba 100644
--- a/ext/cuda/operators_sem_shmem.jl
+++ b/ext/cuda/operators_sem_shmem.jl
@@ -2,25 +2,26 @@ import ClimaCore: DataLayouts, Spaces, Geometry, Operators, Quadratures
 import CUDA
 import ClimaCore.Operators:
     Divergence,
-    WeakDivergence,
     SplitDivergence,
     Gradient,
-    WeakGradient,
-    Curl,
-    WeakCurl
+    Curl
 import ClimaCore.Operators: operator_return_eltype, get_local_geometry
+import ClimaCore.Operators: form_jacobian, form_weighted_arg
 
+# Both forms of the divergence hold one scaled contravariant component per dimension in
+# shared memory, so they share these methods; they differ only in the Jacobian factor that
+# scales the components (see form_jacobian), so the shared arrays named Jv hold J uⁱ for the
+# strong form and WJ uⁱ for the weak form.
 Base.@propagate_inbounds function operator_shmem(
     space,
     ::Val{Nvt},
-    op::Operators.Divergence{(1,)},
+    op::Divergence{(1,)},
     arg,
 ) where {Nvt}
-    FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
     # allocate temp output
-    RT = Operators.operator_return_eltype(op, eltype(arg))
+    RT = operator_return_eltype(op, eltype(arg))
     Jv¹ = CUDA.CuStaticSharedArray(RT, (Nq, Nvt))
     return (Jv¹,)
 end
@@ -30,7 +31,6 @@ Base.@propagate_inbounds function operator_shmem(
     op::Divergence{(1, 2)},
     arg,
 ) where {Nvt}
-    FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
     # allocate temp output
@@ -41,93 +41,33 @@ Base.@propagate_inbounds function operator_shmem(
 end
 
 Base.@propagate_inbounds function operator_fill_shmem!(
-    op::Divergence{(1,)},
+    op::Divergence{(1,), F},
     (Jv¹,),
     space,
     ij,
     slabidx,
     arg,
-)
+) where {F}
     vt = threadIdx().z
     local_geometry = get_local_geometry(space, ij, slabidx)
     i, _ = ij.I
-    (; J) = local_geometry
-    Jv¹[i, vt] = J * Geometry.contravariant1(arg, local_geometry)
+    jacobian = form_jacobian(F(), local_geometry)
+    Jv¹[i, vt] = jacobian * Geometry.contravariant1(arg, local_geometry)
 end
 Base.@propagate_inbounds function operator_fill_shmem!(
-    op::Divergence{(1, 2)},
+    op::Divergence{(1, 2), F},
     (Jv¹, Jv²),
     space,
     ij,
     slabidx,
     arg,
-)
+) where {F}
     vt = threadIdx().z
     local_geometry = get_local_geometry(space, ij, slabidx)
     i, j = ij.I
-    (; J) = local_geometry
-    Jv¹[i, j, vt] = J * Geometry.contravariant1(arg, local_geometry)
-    Jv²[i, j, vt] = J * Geometry.contravariant2(arg, local_geometry)
-end
-
-Base.@propagate_inbounds function operator_shmem(
-    space,
-    ::Val{Nvt},
-    op::WeakDivergence{(1,)},
-    arg,
-) where {Nvt}
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    # allocate temp output
-    RT = operator_return_eltype(op, eltype(arg))
-    WJv¹ = CUDA.CuStaticSharedArray(RT, (Nq, Nvt))
-    return (WJv¹,)
-end
-Base.@propagate_inbounds function operator_shmem(
-    space,
-    ::Val{Nvt},
-    op::WeakDivergence{(1, 2)},
-    arg,
-) where {Nvt}
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    # allocate temp output
-    RT = operator_return_eltype(op, eltype(arg))
-    WJv¹ = CUDA.CuStaticSharedArray(RT, (Nq, Nq, Nvt))
-    WJv² = CUDA.CuStaticSharedArray(RT, (Nq, Nq, Nvt))
-    return (WJv¹, WJv²)
-end
-
-Base.@propagate_inbounds function operator_fill_shmem!(
-    op::WeakDivergence{(1,)},
-    (WJv¹,),
-    space,
-    ij,
-    slabidx,
-    arg,
-)
-    vt = threadIdx().z
-    local_geometry = get_local_geometry(space, ij, slabidx)
-    i, _ = ij.I
-    (; WJ) = local_geometry
-    WJv¹[i, vt] = WJ * Geometry.contravariant1(arg, local_geometry)
-end
-Base.@propagate_inbounds function operator_fill_shmem!(
-    op::WeakDivergence{(1, 2)},
-    (WJv¹, WJv²),
-    space,
-    ij,
-    slabidx,
-    arg,
-)
-    vt = threadIdx().z
-    local_geometry = get_local_geometry(space, ij, slabidx)
-    i, j = ij.I
-    (; WJ) = local_geometry
-    WJv¹[i, j, vt] = WJ * Geometry.contravariant1(arg, local_geometry)
-    WJv²[i, j, vt] = WJ * Geometry.contravariant2(arg, local_geometry)
+    jacobian = form_jacobian(F(), local_geometry)
+    Jv¹[i, j, vt] = jacobian * Geometry.contravariant1(arg, local_geometry)
+    Jv²[i, j, vt] = jacobian * Geometry.contravariant2(arg, local_geometry)
 end
 
 Base.@propagate_inbounds function operator_shmem(
@@ -199,6 +139,10 @@ Base.@propagate_inbounds function operator_fill_shmem!(
     psi[i, j, vt] = arg2
 end
 
+# Both forms of the gradient hold the argument in shared memory, so they share these
+# methods; they differ only in the quadrature weighting applied to it (see
+# form_weighted_arg), so the shared array holds f for the strong form and W f for the weak
+# form. It is wrapped in a tuple to match the other operators.
 Base.@propagate_inbounds function operator_shmem(
     space,
     ::Val{Nvt},
@@ -207,7 +151,8 @@ Base.@propagate_inbounds function operator_shmem(
 ) where {Nvt}
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
-    return CUDA.CuStaticSharedArray(eltype(arg), (Nq, Nvt))
+    f = CUDA.CuStaticSharedArray(eltype(arg), (Nq, Nvt))
+    return (f,)
 end
 Base.@propagate_inbounds function operator_shmem(
     space,
@@ -217,103 +162,51 @@ Base.@propagate_inbounds function operator_shmem(
 ) where {Nvt}
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
-    return CUDA.CuStaticSharedArray(eltype(arg), (Nq, Nq, Nvt))
-end
-
-Base.@propagate_inbounds function operator_fill_shmem!(
-    op::Gradient{(1,)},
-    input,
-    space,
-    ij,
-    slabidx,
-    arg,
-)
-    vt = threadIdx().z
-    i, _ = ij.I
-    input[i, vt] = arg
-end
-Base.@propagate_inbounds function operator_fill_shmem!(
-    op::Gradient{(1, 2)},
-    input,
-    space,
-    ij,
-    slabidx,
-    arg,
-)
-    vt = threadIdx().z
-    i, j = ij.I
-    input[i, j, vt] = arg
-end
-
-Base.@propagate_inbounds function operator_shmem(
-    space,
-    ::Val{Nvt},
-    op::WeakGradient{(1,)},
-    arg,
-) where {Nvt}
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    # allocate temp output
-    IT = eltype(arg)
-    Wf = CUDA.CuStaticSharedArray(IT, (Nq, Nvt))
-    return (Wf,)
-end
-Base.@propagate_inbounds function operator_shmem(
-    space,
-    ::Val{Nvt},
-    op::WeakGradient{(1, 2)},
-    arg,
-) where {Nvt}
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    # allocate temp output
-    IT = eltype(arg)
-    Wf = CUDA.CuStaticSharedArray(IT, (Nq, Nq, Nvt))
-    return (Wf,)
+    f = CUDA.CuStaticSharedArray(eltype(arg), (Nq, Nq, Nvt))
+    return (f,)
 end
 
 Base.@propagate_inbounds function operator_fill_shmem!(
-    op::WeakGradient{(1,)},
-    (Wf,),
+    op::Gradient{(1,), F},
+    (f,),
     space,
     ij,
     slabidx,
     arg,
-)
+) where {F}
     vt = threadIdx().z
     local_geometry = get_local_geometry(space, ij, slabidx)
-    W = local_geometry.WJ * local_geometry.invJ
     i, _ = ij.I
-    Wf[i, vt] = W * arg
+    f[i, vt] = form_weighted_arg(F(), local_geometry, arg)
 end
 Base.@propagate_inbounds function operator_fill_shmem!(
-    op::WeakGradient{(1, 2)},
-    (Wf,),
+    op::Gradient{(1, 2), F},
+    (f,),
     space,
     ij,
     slabidx,
     arg,
-)
+) where {F}
     vt = threadIdx().z
     local_geometry = get_local_geometry(space, ij, slabidx)
-    W = local_geometry.WJ * local_geometry.invJ
     i, j = ij.I
-    Wf[i, j, vt] = W * arg
+    f[i, j, vt] = form_weighted_arg(F(), local_geometry, arg)
 end
 
+# Both forms of the curl hold the covariant components of the argument in shared memory, so
+# they share these methods; they differ only in the quadrature weighting applied to those
+# components (see form_weighted_arg). `curl_result_type` always returns a
+# Contravariant123Vector, so all three components are allocated in 2D, while the first is
+# unused in 1D.
 Base.@propagate_inbounds function operator_shmem(
     space,
     ::Val{Nvt},
     op::Curl{(1,)},
     arg,
 ) where {Nvt}
-    FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
-    IT = eltype(arg)
-    ET = eltype(IT)
+    ET = eltype(eltype(arg))
     v₂ = CUDA.CuStaticSharedArray(ET, (Nq, Nvt))
     v₃ = CUDA.CuStaticSharedArray(ET, (Nq, Nvt))
     return (nothing, v₂, v₃)
@@ -324,11 +217,9 @@ Base.@propagate_inbounds function operator_shmem(
     op::Curl{(1, 2)},
     arg,
 ) where {Nvt}
-    FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
-    IT = eltype(arg)
-    ET = eltype(IT)
+    ET = eltype(eltype(arg))
     v₁ = CUDA.CuStaticSharedArray(ET, (Nq, Nq, Nvt))
     v₂ = CUDA.CuStaticSharedArray(ET, (Nq, Nq, Nvt))
     v₃ = CUDA.CuStaticSharedArray(ET, (Nq, Nq, Nvt))
@@ -336,102 +227,35 @@ Base.@propagate_inbounds function operator_shmem(
 end
 
 Base.@propagate_inbounds function operator_fill_shmem!(
-    op::Curl{(1,)},
+    op::Curl{(1,), F},
     work,
     space,
     ij,
     slabidx,
     arg,
-)
+) where {F}
     vt = threadIdx().z
     i, _ = ij.I
     local_geometry = get_local_geometry(space, ij, slabidx)
     _, v₂, v₃ = work
-    v₂[i, vt] = Geometry.covariant2(arg, local_geometry)
-    v₃[i, vt] = Geometry.covariant3(arg, local_geometry)
+    weighted(x) = form_weighted_arg(F(), local_geometry, x)
+    v₂[i, vt] = weighted(Geometry.covariant2(arg, local_geometry))
+    v₃[i, vt] = weighted(Geometry.covariant3(arg, local_geometry))
 end
 Base.@propagate_inbounds function operator_fill_shmem!(
-    op::Curl{(1, 2)},
+    op::Curl{(1, 2), F},
     work,
     space,
     ij,
     slabidx,
     arg,
-)
+) where {F}
     vt = threadIdx().z
     i, j = ij.I
     local_geometry = get_local_geometry(space, ij, slabidx)
     v₁, v₂, v₃ = work
-    v₁[i, j, vt] = Geometry.covariant1(arg, local_geometry)
-    v₂[i, j, vt] = Geometry.covariant2(arg, local_geometry)
-    v₃[i, j, vt] = Geometry.covariant3(arg, local_geometry)
-end
-
-Base.@propagate_inbounds function operator_shmem(
-    space,
-    ::Val{Nvt},
-    op::WeakCurl{(1,)},
-    arg,
-) where {Nvt}
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    IT = eltype(arg)
-    ET = eltype(IT)
-    # `curl_result_type` always returns Contravariant123Vector. Wv₁ is unused
-    # for WeakCurl{(1,)}.
-    Wv₂ = CUDA.CuStaticSharedArray(ET, (Nq, Nvt))
-    Wv₃ = CUDA.CuStaticSharedArray(ET, (Nq, Nvt))
-    return (nothing, Wv₂, Wv₃)
-end
-Base.@propagate_inbounds function operator_shmem(
-    space,
-    ::Val{Nvt},
-    op::WeakCurl{(1, 2)},
-    arg,
-) where {Nvt}
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    IT = eltype(arg)
-    ET = eltype(IT)
-    # `curl_result_type` always returns Contravariant123Vector; allocate all three.
-    Wv₁ = CUDA.CuStaticSharedArray(ET, (Nq, Nq, Nvt))
-    Wv₂ = CUDA.CuStaticSharedArray(ET, (Nq, Nq, Nvt))
-    Wv₃ = CUDA.CuStaticSharedArray(ET, (Nq, Nq, Nvt))
-    return (Wv₁, Wv₂, Wv₃)
-end
-
-Base.@propagate_inbounds function operator_fill_shmem!(
-    op::WeakCurl{(1,)},
-    work,
-    space,
-    ij,
-    slabidx,
-    arg,
-)
-    vt = threadIdx().z
-    i, _ = ij.I
-    local_geometry = get_local_geometry(space, ij, slabidx)
-    W = local_geometry.WJ * local_geometry.invJ
-    _, Wv₂, Wv₃ = work
-    Wv₂[i, vt] = W * Geometry.covariant2(arg, local_geometry)
-    Wv₃[i, vt] = W * Geometry.covariant3(arg, local_geometry)
-end
-Base.@propagate_inbounds function operator_fill_shmem!(
-    op::WeakCurl{(1, 2)},
-    work,
-    space,
-    ij,
-    slabidx,
-    arg,
-)
-    vt = threadIdx().z
-    i, j = ij.I
-    local_geometry = get_local_geometry(space, ij, slabidx)
-    W = local_geometry.WJ * local_geometry.invJ
-    Wv₁, Wv₂, Wv₃ = work
-    Wv₁[i, j, vt] = W * Geometry.covariant1(arg, local_geometry)
-    Wv₂[i, j, vt] = W * Geometry.covariant2(arg, local_geometry)
-    Wv₃[i, j, vt] = W * Geometry.covariant3(arg, local_geometry)
+    weighted(x) = form_weighted_arg(F(), local_geometry, x)
+    v₁[i, j, vt] = weighted(Geometry.covariant1(arg, local_geometry))
+    v₂[i, j, vt] = weighted(Geometry.covariant2(arg, local_geometry))
+    v₃[i, j, vt] = weighted(Geometry.covariant3(arg, local_geometry))
 end
diff --git a/ext/cuda/operators_spectral_element.jl b/ext/cuda/operators_spectral_element.jl
index 4be97e1e53..2a4e98130b 100644
--- a/ext/cuda/operators_spectral_element.jl
+++ b/ext/cuda/operators_spectral_element.jl
@@ -5,8 +5,9 @@ using CUDA
 import ClimaCore.Operators: AbstractSpectralStyle, strip_space
 import ClimaCore.Operators: SpectralBroadcasted, set_node!, get_node
 import ClimaCore.Operators: get_local_geometry
+import ClimaCore.Operators: Divergence, SplitDivergence, Gradient, Curl
 import ClimaCore.Operators:
-    Divergence, WeakDivergence, SplitDivergence, Gradient, WeakGradient, Curl, WeakCurl
+    form_deriv_entry, form_jacobian_rescale, form_weight_rescale
 import Base.Broadcast: Broadcasted
 
 """
@@ -181,13 +182,19 @@ Base.@propagate_inbounds function resolve_shmem!(obj, ij, slabidx)
     nothing
 end
 
+# The methods below serve both forms of each operator: form_deriv_entry supplies the
+# derivative matrix entries (transposed and sign-flipped for the weak form), and
+# form_jacobian_rescale or form_weight_rescale divides the form's Jacobian or
+# quadrature-weight factor back out of the result. Every operator keeps one accumulator
+# per dimension and combines them once at the end, which keeps the accumulation loops as
+# independent dependency chains.
 Base.@propagate_inbounds function operator_evaluate(
-    op::Divergence{(1,)},
+    op::Divergence{(1,), F},
     (Jv¹,),
     space,
     ij,
     slabidx,
-)
+) where {F}
     vt = threadIdx().z
     i, _ = ij.I
 
@@ -198,69 +205,19 @@ Base.@propagate_inbounds function operator_evaluate(
 
     local_geometry = get_local_geometry(space, ij, slabidx)
 
-    DJv = D[i, 1] * Jv¹[1, vt]
+    D₁Jv¹ = form_deriv_entry(F(), D, i, 1) * Jv¹[1, vt]
     for k in 2:Nq
-        DJv += D[i, k] * Jv¹[k, vt]
+        D₁Jv¹ += form_deriv_entry(F(), D, i, k) * Jv¹[k, vt]
     end
-    return DJv * local_geometry.invJ
+    return form_jacobian_rescale(F(), local_geometry, D₁Jv¹)
 end
 Base.@propagate_inbounds function operator_evaluate(
-    op::Divergence{(1, 2)},
+    op::Divergence{(1, 2), F},
     (Jv¹, Jv²),
     space,
     ij,
     slabidx,
-)
-    vt = threadIdx().z
-    i, j = ij.I
-
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    D = Quadratures.differentiation_matrix(FT, QS)
-
-    local_geometry = get_local_geometry(space, ij, slabidx)
-
-    DJv = D[i, 1] * Jv¹[1, j, vt]
-    for k in 2:Nq
-        DJv += D[i, k] * Jv¹[k, j, vt]
-    end
-    for k in 1:Nq
-        DJv += D[j, k] * Jv²[i, k, vt]
-    end
-    return DJv * local_geometry.invJ
-end
-
-Base.@propagate_inbounds function operator_evaluate(
-    op::WeakDivergence{(1,)},
-    (WJv¹,),
-    space,
-    ij,
-    slabidx,
-)
-    vt = CUDA.threadIdx().z
-    i, _ = ij.I
-
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    D = Quadratures.differentiation_matrix(FT, QS)
-
-    local_geometry = get_local_geometry(space, ij, slabidx)
-
-    Dᵀ₁WJv¹ = D[1, i] * WJv¹[1, vt]
-    for k in 2:Nq
-        Dᵀ₁WJv¹ += D[k, i] * WJv¹[k, vt]
-    end
-    return -Dᵀ₁WJv¹ / local_geometry.WJ
-end
-Base.@propagate_inbounds function operator_evaluate(
-    op::WeakDivergence{(1, 2)},
-    (WJv¹, WJv²),
-    space,
-    ij,
-    slabidx,
-)
+) where {F}
     vt = threadIdx().z
     i, j = ij.I
 
@@ -271,13 +228,13 @@ Base.@propagate_inbounds function operator_evaluate(
 
     local_geometry = get_local_geometry(space, ij, slabidx)
 
-    Dᵀ₁WJv¹ = D[1, i] * WJv¹[1, j, vt]
-    Dᵀ₂WJv² = D[1, j] * WJv²[i, 1, vt]
+    D₁Jv¹ = form_deriv_entry(F(), D, i, 1) * Jv¹[1, j, vt]
+    D₂Jv² = form_deriv_entry(F(), D, j, 1) * Jv²[i, 1, vt]
     for k in 2:Nq
-        Dᵀ₁WJv¹ += D[k, i] * WJv¹[k, j, vt]
-        Dᵀ₂WJv² += D[k, j] * WJv²[i, k, vt]
+        D₁Jv¹ += form_deriv_entry(F(), D, i, k) * Jv¹[k, j, vt]
+        D₂Jv² += form_deriv_entry(F(), D, j, k) * Jv²[i, k, vt]
     end
-    return -(Dᵀ₁WJv¹ + Dᵀ₂WJv²) / local_geometry.WJ
+    return form_jacobian_rescale(F(), local_geometry, D₁Jv¹ + D₂Jv²)
 end
 
 Base.@propagate_inbounds function operator_evaluate(
@@ -341,12 +298,12 @@ Base.@propagate_inbounds function operator_evaluate(
 end
 
 Base.@propagate_inbounds function operator_evaluate(
-    op::Gradient{(1,)},
-    input,
+    op::Gradient{(1,), F},
+    (input,),
     space,
     ij,
     slabidx,
-)
+) where {F}
     vt = threadIdx().z
     i, _ = ij.I
 
@@ -355,29 +312,32 @@ Base.@propagate_inbounds function operator_evaluate(
     Nq = Quadratures.degrees_of_freedom(QS)
     D = Quadratures.differentiation_matrix(FT, QS)
 
+    local_geometry = get_local_geometry(space, ij, slabidx)
+
     @inbounds begin
-        ∂f∂ξ₁ = D[i, 1] * input[1, vt]
+        ∂f∂ξ₁ = form_deriv_entry(F(), D, i, 1) * input[1, vt]
         for k in 2:Nq
-            ∂f∂ξ₁ += D[i, k] * input[k, vt]
+            ∂f∂ξ₁ += form_deriv_entry(F(), D, i, k) * input[k, vt]
         end
     end
-    if eltype(input) <: Number
-        return Geometry.Covariant1Vector(∂f∂ξ₁)
+    result = if eltype(input) <: Number
+        Geometry.Covariant1Vector(∂f∂ξ₁)
     elseif eltype(input) <: Geometry.AbstractTensor{1}
         tensor_axes = (Geometry.Covariant1Axis(), Geometry.tensor_axes(eltype(input))[1])
         tensor_components = hcat(parent(∂f∂ξ₁))'
-        return Geometry.Tensor(tensor_components, tensor_axes)
+        Geometry.Tensor(tensor_components, tensor_axes)
     else
         error("Unsupported input type for gradient operator: $(eltype(input))")
     end
+    return form_weight_rescale(F(), local_geometry, result)
 end
 Base.@propagate_inbounds function operator_evaluate(
-    op::Gradient{(1, 2)},
-    input,
+    op::Gradient{(1, 2), F},
+    (input,),
     space,
     ij,
     slabidx,
-)
+) where {F}
     vt = threadIdx().z
     i, j = ij.I
 
@@ -386,84 +346,36 @@ Base.@propagate_inbounds function operator_evaluate(
     Nq = Quadratures.degrees_of_freedom(QS)
     D = Quadratures.differentiation_matrix(FT, QS)
 
+    local_geometry = get_local_geometry(space, ij, slabidx)
+
     @inbounds begin
-        ∂f∂ξ₁ = D[i, 1] * input[1, j, vt]
-        ∂f∂ξ₂ = D[j, 1] * input[i, 1, vt]
+        ∂f∂ξ₁ = form_deriv_entry(F(), D, i, 1) * input[1, j, vt]
+        ∂f∂ξ₂ = form_deriv_entry(F(), D, j, 1) * input[i, 1, vt]
         for k in 2:Nq
-            ∂f∂ξ₁ += D[i, k] * input[k, j, vt]
-            ∂f∂ξ₂ += D[j, k] * input[i, k, vt]
+            ∂f∂ξ₁ += form_deriv_entry(F(), D, i, k) * input[k, j, vt]
+            ∂f∂ξ₂ += form_deriv_entry(F(), D, j, k) * input[i, k, vt]
         end
     end
-    if eltype(input) <: Number
-        return Geometry.Covariant12Vector(∂f∂ξ₁, ∂f∂ξ₂)
+    result = if eltype(input) <: Number
+        Geometry.Covariant12Vector(∂f∂ξ₁, ∂f∂ξ₂)
     elseif eltype(input) <: Geometry.AbstractTensor{1}
         tensor_axes = (Geometry.Covariant12Axis(), Geometry.tensor_axes(eltype(input))[1])
         tensor_components =
             hcat(parent(∂f∂ξ₁), parent(∂f∂ξ₂))'
-        return Geometry.Tensor(tensor_components, tensor_axes)
+        Geometry.Tensor(tensor_components, tensor_axes)
     else
         error("Unsupported input type for gradient operator: $(eltype(input))")
     end
+    return form_weight_rescale(F(), local_geometry, result)
 end
 
 Base.@propagate_inbounds function operator_evaluate(
-    op::WeakGradient{(1,)},
-    (Wf,),
-    space,
-    ij,
-    slabidx,
-)
-    vt = threadIdx().z
-    i, _ = ij.I
-
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    D = Quadratures.differentiation_matrix(FT, QS)
-
-    local_geometry = get_local_geometry(space, ij, slabidx)
-    W = local_geometry.WJ * local_geometry.invJ
-
-    Dᵀ₁Wf = D[1, i] * Wf[1, vt]
-    for k in 2:Nq
-        Dᵀ₁Wf += D[k, i] * Wf[k, vt]
-    end
-    return Geometry.Covariant1Vector(-Dᵀ₁Wf) / W
-end
-Base.@propagate_inbounds function operator_evaluate(
-    op::WeakGradient{(1, 2)},
-    (Wf,),
-    space,
-    ij,
-    slabidx,
-)
-    vt = threadIdx().z
-    i, j = ij.I
-
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    D = Quadratures.differentiation_matrix(FT, QS)
-
-    local_geometry = get_local_geometry(space, ij, slabidx)
-    W = local_geometry.WJ * local_geometry.invJ
-
-    Dᵀ₁Wf = D[1, i] * Wf[1, j, vt]
-    Dᵀ₂Wf = D[1, j] * Wf[i, 1, vt]
-    for k in 2:Nq
-        Dᵀ₁Wf += D[k, i] * Wf[k, j, vt]
-        Dᵀ₂Wf += D[k, j] * Wf[i, k, vt]
-    end
-    return Geometry.Covariant12Vector(-Dᵀ₁Wf, -Dᵀ₂Wf) / W
-end
-
-Base.@propagate_inbounds function operator_evaluate(
-    op::Curl{(1,)},
+    op::Curl{(1,), F},
     work,
     space,
     ij,
     slabidx,
-)
+) where {F}
     vt = threadIdx().z
     i, _ = ij.I
 
@@ -474,22 +386,22 @@ Base.@propagate_inbounds function operator_evaluate(
     local_geometry = get_local_geometry(space, ij, slabidx)
 
     _, v₂, v₃ = work
-    D₁v₂ = D[i, 1] * v₂[1, vt]
-    D₁v₃ = D[i, 1] * v₃[1, vt]
+    D₁v₂ = form_deriv_entry(F(), D, i, 1) * v₂[1, vt]
+    D₁v₃ = form_deriv_entry(F(), D, i, 1) * v₃[1, vt]
     @simd for k in 2:Nq
-        D₁v₂ += D[i, k] * v₂[k, vt]
-        D₁v₃ += D[i, k] * v₃[k, vt]
+        D₁v₂ += form_deriv_entry(F(), D, i, k) * v₂[k, vt]
+        D₁v₃ += form_deriv_entry(F(), D, i, k) * v₃[k, vt]
     end
     result = Geometry.Contravariant123Vector(zero(FT), -D₁v₃, D₁v₂)
-    return result * local_geometry.invJ
+    return form_jacobian_rescale(F(), local_geometry, result)
 end
 Base.@propagate_inbounds function operator_evaluate(
-    op::Curl{(1, 2)},
+    op::Curl{(1, 2), F},
     work,
     space,
     ij,
     slabidx,
-)
+) where {F}
     vt = threadIdx().z
     i, j = ij.I
 
@@ -500,73 +412,16 @@ Base.@propagate_inbounds function operator_evaluate(
     local_geometry = get_local_geometry(space, ij, slabidx)
 
     v₁, v₂, v₃ = work
-    D₁v₂ = D[i, 1] * v₂[1, j, vt]
-    D₂v₁ = D[j, 1] * v₁[i, 1, vt]
-    D₁v₃ = D[i, 1] * v₃[1, j, vt]
-    D₂v₃ = D[j, 1] * v₃[i, 1, vt]
+    D₁v₂ = form_deriv_entry(F(), D, i, 1) * v₂[1, j, vt]
+    D₂v₁ = form_deriv_entry(F(), D, j, 1) * v₁[i, 1, vt]
+    D₁v₃ = form_deriv_entry(F(), D, i, 1) * v₃[1, j, vt]
+    D₂v₃ = form_deriv_entry(F(), D, j, 1) * v₃[i, 1, vt]
     @simd for k in 2:Nq
-        D₁v₂ += D[i, k] * v₂[k, j, vt]
-        D₂v₁ += D[j, k] * v₁[i, k, vt]
-        D₁v₃ += D[i, k] * v₃[k, j, vt]
-        D₂v₃ += D[j, k] * v₃[i, k, vt]
+        D₁v₂ += form_deriv_entry(F(), D, i, k) * v₂[k, j, vt]
+        D₂v₁ += form_deriv_entry(F(), D, j, k) * v₁[i, k, vt]
+        D₁v₃ += form_deriv_entry(F(), D, i, k) * v₃[k, j, vt]
+        D₂v₃ += form_deriv_entry(F(), D, j, k) * v₃[i, k, vt]
     end
     result = Geometry.Contravariant123Vector(D₂v₃, -D₁v₃, D₁v₂ - D₂v₁)
-    return result * local_geometry.invJ
-end
-
-Base.@propagate_inbounds function operator_evaluate(
-    op::WeakCurl{(1,)},
-    work,
-    space,
-    ij,
-    slabidx,
-)
-    vt = threadIdx().z
-    i, _ = ij.I
-
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    D = Quadratures.differentiation_matrix(FT, QS)
-    local_geometry = get_local_geometry(space, ij, slabidx)
-
-    _, Wv₂, Wv₃ = work
-    Dᵀ₁Wv₂ = D[1, i] * Wv₂[1, vt]
-    Dᵀ₁Wv₃ = D[1, i] * Wv₃[1, vt]
-    @simd for k in 2:Nq
-        Dᵀ₁Wv₂ += D[k, i] * Wv₂[k, vt]
-        Dᵀ₁Wv₃ += D[k, i] * Wv₃[k, vt]
-    end
-    result = Geometry.Contravariant123Vector(zero(FT), Dᵀ₁Wv₃, -Dᵀ₁Wv₂)
-    return result / local_geometry.WJ
-end
-Base.@propagate_inbounds function operator_evaluate(
-    op::WeakCurl{(1, 2)},
-    work,
-    space,
-    ij,
-    slabidx,
-)
-    vt = threadIdx().z
-    i, j = ij.I
-
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    D = Quadratures.differentiation_matrix(FT, QS)
-    local_geometry = get_local_geometry(space, ij, slabidx)
-
-    Wv₁, Wv₂, Wv₃ = work
-    Dᵀ₁Wv₂ = D[1, i] * Wv₂[1, j, vt]
-    Dᵀ₂Wv₁ = D[1, j] * Wv₁[i, 1, vt]
-    Dᵀ₁Wv₃ = D[1, i] * Wv₃[1, j, vt]
-    Dᵀ₂Wv₃ = D[1, j] * Wv₃[i, 1, vt]
-    @simd for k in 2:Nq
-        Dᵀ₁Wv₂ += D[k, i] * Wv₂[k, j, vt]
-        Dᵀ₂Wv₁ += D[k, j] * Wv₁[i, k, vt]
-        Dᵀ₁Wv₃ += D[k, i] * Wv₃[k, j, vt]
-        Dᵀ₂Wv₃ += D[k, j] * Wv₃[i, k, vt]
-    end
-    result = Geometry.Contravariant123Vector(-Dᵀ₂Wv₃, Dᵀ₁Wv₃, Dᵀ₂Wv₁ - Dᵀ₁Wv₂)
-    return result / local_geometry.WJ
+    return form_jacobian_rescale(F(), local_geometry, result)
 end
diff --git a/src/Operators/Operators.jl b/src/Operators/Operators.jl
index 39adc43e54..033123158d 100644
--- a/src/Operators/Operators.jl
+++ b/src/Operators/Operators.jl
@@ -26,5 +26,6 @@ include("finitedifference.jl")
 include("remapping.jl")
 include("integrals.jl")
 include("columnwise.jl")
+include("deprecated.jl")
 
 end # module
diff --git a/src/Operators/deprecated.jl b/src/Operators/deprecated.jl
new file mode 100644
index 0000000000..5011cd9e62
--- /dev/null
+++ b/src/Operators/deprecated.jl
@@ -0,0 +1,124 @@
+# Backwards-compatibility aliases for the pre-FormType spectral element
+# operators. The weak-form operators are now the WeakForm variants of their
+# strong-form counterparts (e.g., Divergence{I, WeakForm} instead of a separate
+# WeakDivergence type); the aliases below keep the old names working. As with
+# src/DataLayouts/deprecated.jl, these are plain aliases without deprecation
+# warnings. Remove this file when all downstream consumers have migrated to the
+# parameterized names (type-alias cleanup, Phase 4 of
+# https://github.com/CliMA/ClimaCore.jl/issues/2554).
+
+"""
+    wdiv = WeakDivergence()
+    wdiv.(u)
+
+Computes the "weak divergence" of a vector field `u`.
+
+This is defined as the scalar field ``\\theta \\in \\mathcal{V}_0`` such that
+for all ``\\phi\\in \\mathcal{V}_0``
+```math
+\\int_\\Omega \\phi \\theta \\, d \\Omega
+=
+- \\int_\\Omega (\\nabla \\phi) \\cdot u \\,d \\Omega
+```
+where ``\\mathcal{V}_0`` is the space of ``u``.
+
+This arises as the contribution of the volume integral after applying
+integration by parts to the weak form expression of the divergence
+```math
+\\int_\\Omega \\phi (\\nabla \\cdot u) \\, d \\Omega
+=
+- \\int_\\Omega (\\nabla \\phi) \\cdot u \\,d \\Omega
++ \\oint_{\\partial \\Omega} \\phi (u \\cdot n) \\,d \\sigma
+```
+
+It can be written in matrix form as
+```math
+ϕ^\\top WJ θ = - \\sum_i (D_i ϕ)^\\top WJ u^i
+```
+which reduces to
+```math
+θ = -(WJ)^{-1} \\sum_i D_i^\\top WJ u^i
+```
+where
+ - ``J`` is the diagonal Jacobian matrix
+ - ``W`` is the diagonal matrix of quadrature weights
+ - ``D_i`` is the derivative matrix along the ``i``th dimension
+"""
+const WeakDivergence{I} = Divergence{I, WeakForm}
+WeakDivergence() = WeakDivergence{()}()
+
+"""
+    wgrad = WeakGradient()
+    wgrad.(f)
+
+Compute the "weak gradient" of `f` on each element.
+
+This is defined as the the vector field ``\\theta \\in \\mathcal{V}_0`` such
+that for all ``\\phi \\in \\mathcal{V}_0``
+```math
+\\int_\\Omega \\phi \\cdot \\theta \\, d \\Omega
+=
+- \\int_\\Omega (\\nabla \\cdot \\phi) f \\, d\\Omega
+```
+where ``\\mathcal{V}_0`` is the space of ``f``.
+
+This arises from the contribution of the volume integral after by applying
+integration by parts to the weak form expression of the gradient
+```math
+\\int_\\Omega \\phi \\cdot (\\nabla f) \\, d \\Omega
+=
+- \\int_\\Omega f (\\nabla \\cdot \\phi) \\, d\\Omega
++ \\oint_{\\partial \\Omega} f (\\phi \\cdot n) \\, d \\sigma
+```
+
+In matrix form, this becomes
+```math
+{\\phi^i}^\\top W J \\theta_i = - ( J^{-1} D_i J \\phi^i )^\\top W J f
+```
+which reduces to
+```math
+\\theta_i = -W^{-1} D_i^\\top W f
+```
+where ``D_i`` is the derivative matrix along the ``i``th dimension.
+"""
+const WeakGradient{I} = Gradient{I, WeakForm}
+WeakGradient() = WeakGradient{()}()
+
+"""
+    wcurl = WeakCurl()
+    wcurl.(u)
+
+Computes the "weak curl" on each element of a covariant vector field `u`.
+
+Note: The vector field ``u`` needs to be excliclty converted to a `CovaraintVector`,
+as then the `WeakCurl` is independent of the local metric tensor.
+
+This is defined as the vector field ``\\theta \\in \\mathcal{V}_0`` such that
+for all ``\\phi \\in \\mathcal{V}_0``
+```math
+\\int_\\Omega \\phi \\cdot \\theta \\, d \\Omega
+=
+\\int_\\Omega (\\nabla \\times \\phi) \\cdot u \\,d \\Omega
+```
+where ``\\mathcal{V}_0`` is the space of ``f``.
+
+This arises from the contribution of the volume integral after by applying
+integration by parts to the weak form expression of the curl
+```math
+\\int_\\Omega \\phi \\cdot (\\nabla \\times u) \\,d\\Omega
+=
+\\int_\\Omega (\\nabla \\times \\phi) \\cdot u \\,d \\Omega
+- \\oint_{\\partial \\Omega} (\\phi \\times u) \\cdot n \\,d\\sigma
+```
+
+In matrix form, this becomes
+```math
+{\\phi_i}^\\top W J \\theta^i = (J^{-1} \\epsilon^{kji} D_j \\phi_i)^\\top W J u_k
+```
+which, by using the anti-symmetry of the Levi-Civita symbol, reduces to
+```math
+\\theta^i = - \\epsilon^{ijk} (WJ)^{-1} D_j^\\top W u_k
+```
+"""
+const WeakCurl{I} = Curl{I, WeakForm}
+WeakCurl() = WeakCurl{()}()
diff --git a/src/Operators/spectralelement.jl b/src/Operators/spectralelement.jl
index 361cec1cd1..271fea3bd1 100644
--- a/src/Operators/spectralelement.jl
+++ b/src/Operators/spectralelement.jl
@@ -99,6 +99,120 @@ Adapt.adapt_structure(to, sbc::SpectralBroadcasted{Style}) where {Style} =
 
 return_space(::SpectralElementOperator, space, args...) = space
 
+"""
+    FormType
+
+Supertype of the singleton types [`StrongForm`](@ref) and [`WeakForm`](@ref),
+which distinguish the variational form of a spectral element operator.
+
+The strong and weak variants of an operator share the same interior
+computation; they differ only in three form-dependent factors:
+ - whether the derivative matrix is applied directly or transposed with a sign
+   flip (from integration by parts); see `form_deriv_entry`,
+ - whether the argument is weighted by the quadrature weights `W` or by the
+   Jacobian factor; see `form_weighted_arg` and `form_jacobian`,
+ - whether the result is rescaled by `J` or by `WJ`; see
+   `form_jacobian_rescale`, and by `W` or not at all; see
+   `form_weight_rescale`.
+
+Operators with strong/weak variants carry a `FormType` as their second type
+parameter, e.g. `Divergence{I, StrongForm}`, with the weak variant available
+under an alias, e.g. `WeakDivergence{I} = Divergence{I, WeakForm}`.
+"""
+abstract type FormType end
+
+"""
+    StrongForm()
+
+The [`FormType`](@ref) of an operator that discretizes a derivative directly at
+the quadrature points (e.g. [`Divergence`](@ref), [`Gradient`](@ref),
+[`Curl`](@ref)).
+"""
+struct StrongForm <: FormType end
+
+"""
+    WeakForm()
+
+The [`FormType`](@ref) of an operator that discretizes the volume-integral
+contribution of the corresponding weak-form expression, obtained after
+integration by parts (e.g. [`WeakDivergence`](@ref), [`WeakGradient`](@ref),
+[`WeakCurl`](@ref)).
+"""
+struct WeakForm <: FormType end
+
+"""
+    form_deriv_entry(form, D, ii, i)
+
+Entry of the derivative matrix `D` applied by an operator of the given
+[`FormType`](@ref) when accumulating the contribution of quadrature point `i`
+to quadrature point `ii`: `D[ii, i]` for the strong form, and `-D[i, ii]` (the
+transpose, with the sign flip from integration by parts) for the weak form.
+"""
+@inline form_deriv_entry(::StrongForm, D, ii, i) = D[ii, i]
+@inline form_deriv_entry(::WeakForm, D, ii, i) = -D[i, ii]
+
+"""
+    form_weighted_arg(form, local_geometry, x)
+
+The argument value `x`, weighted as required by an operator of the given
+[`FormType`](@ref) whose weak variant integrates against test functions: `x`
+itself for the strong form, and `W * x` (with the quadrature weights
+`W = WJ * J⁻¹`) for the weak form. Used by [`Gradient`](@ref) and
+[`Curl`](@ref).
+"""
+@inline form_weighted_arg(::StrongForm, local_geometry, x) = x
+@inline form_weighted_arg(::WeakForm, local_geometry, x) =
+    (local_geometry.WJ * local_geometry.invJ) * x
+
+"""
+    form_jacobian(form, local_geometry)
+
+The Jacobian factor that scales the contravariant components summed by an
+operator of the given [`FormType`](@ref): `J` for the strong form, and `WJ`
+for the weak form. Used by [`Divergence`](@ref).
+"""
+@inline form_jacobian(::StrongForm, local_geometry) = local_geometry.J
+@inline form_jacobian(::WeakForm, local_geometry) = local_geometry.WJ
+
+"""
+    form_jacobian_rescale(form, local_geometry, x)
+
+The result value `x`, divided by the `form_jacobian` of the given
+[`FormType`](@ref): `x * J⁻¹` for the strong form (using the precomputed
+inverse), and `x / WJ` for the weak form. Used by [`Divergence`](@ref) and
+[`Curl`](@ref).
+"""
+@inline form_jacobian_rescale(::StrongForm, local_geometry, x) =
+    x * local_geometry.invJ
+@inline form_jacobian_rescale(::WeakForm, local_geometry, x) =
+    x / local_geometry.WJ
+
+"""
+    form_weight_rescale(form, local_geometry, x)
+
+The result value `x`, divided by the quadrature weights `W = WJ * J⁻¹` if the
+given [`FormType`](@ref) requires it: `x` itself for the strong form, and
+`x / W` for the weak form. Used by [`Gradient`](@ref), whose weak variant
+weights its argument by `W` without a Jacobian factor to divide out.
+
+The CPU `apply_operator` methods for [`Gradient`](@ref) inline this rescale
+behind an `F === WeakForm` branch instead of calling it, so that the strong form
+skips the loop over quadrature points entirely; on GPUs each thread rescales
+only its own point, so there is no loop to skip.
+"""
+@inline form_weight_rescale(::StrongForm, local_geometry, x) = x
+@inline form_weight_rescale(::WeakForm, local_geometry, x) =
+    x / (local_geometry.WJ * local_geometry.invJ)
+
+"""
+    rebuild_operator(op, space)
+
+Reconstruct a `SpectralElementOperator` with its `operator_axes` reset to those
+of `space`, preserving all other type parameters (in particular the
+[`FormType`](@ref) of operators with strong/weak variants).
+"""
+rebuild_operator(op::SpectralElementOperator, space) =
+    unionall_type(typeof(op)){()}(space)
 
 function Base.Broadcast.broadcasted(op::SpectralElementOperator, args...)
     args′ = map(Base.Broadcast.broadcastable, args)
@@ -143,9 +257,8 @@ function Base.Broadcast.instantiate(sbc::SpectralBroadcasted)
         RT = operator_return_eltype(op, map(eltype, args)...)
         return Broadcast.broadcasted(Returns(zero(RT)), Fields.coordinate_field(axes))
     end
-    # If we've already instantiated, then we need to strip the type parameters,
-    # for example, `Divergence{()}(axes)`.
-    op = unionall_type(typeof(op)){()}(axes)
+    # If we've already instantiated, then we need to reset the operator axes.
+    op = rebuild_operator(op, axes)
     Style = AbstractSpectralStyle(ClimaComms.device(axes))
     return SpectralBroadcasted{Style}(op, args, axes)
 end
@@ -164,13 +277,13 @@ function Base.Broadcast.instantiate(
     end
     # For FiniteDifferenceSpace with operators, return zeros for horizontal operators
     if axes isa Spaces.FiniteDifferenceSpace && bc.f isa SpectralElementOperator
-        op = unionall_type(typeof(bc.f)){()}(axes)
+        op = rebuild_operator(bc.f, axes)
         RT = operator_return_eltype(op, map(eltype, args)...)
         return Broadcast.broadcasted(Returns(zero(RT)), Fields.coordinate_field(axes))
     end
 
     if bc.f isa SpectralElementOperator
-        op = unionall_type(typeof(bc.f)){()}(axes)
+        op = rebuild_operator(bc.f, axes)
         Style = AbstractSpectralStyle(ClimaComms.device(axes))
         return Base.Broadcast.Broadcasted{Style}(op, args, axes)
     else
@@ -578,14 +691,21 @@ where ``D_i`` is the derivative matrix along the ``i``th dimension
 ## References
 - [Taylor2010](@cite), equation 15
 """
-struct Divergence{I} <: SpectralElementOperator{I} end
-Divergence() = Divergence{()}()
-Divergence{()}(space) = Divergence{operator_axes(space)}()
+struct Divergence{I, F <: FormType} <: SpectralElementOperator{I} end
+Divergence() = Divergence{(), StrongForm}()
+Divergence{I}() where {I} = Divergence{I, StrongForm}()
+rebuild_operator(::Divergence{I, F}, space) where {I, F} =
+    Divergence{operator_axes(space), F}()
 
 operator_return_eltype(op::Divergence{I}, ::Type{S}) where {I, S} =
     Geometry.divergence_result_type(S)
 
-function apply_operator(op::Divergence{(1,)}, space, slabidx, arg)
+# The strong divergence is J⁻¹ ∑ᵢ Dᵢ (J uⁱ), while the weak divergence is
+# -(WJ)⁻¹ ∑ᵢ Dᵢᵀ (WJ uⁱ); see form_deriv_entry, form_jacobian, and
+# form_jacobian_rescale for the form-dependent factors.
+
+function apply_operator(op::Divergence{(1,), F}, space, slabidx, arg) where {F}
+    form = F()
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -598,25 +718,28 @@ function apply_operator(op::Divergence{(1,)}, space, slabidx, arg)
         ij = CartesianIndex((i,))
         local_geometry = get_local_geometry(space, ij, slabidx)
         v = get_node(space, arg, ij, slabidx)
-        Jv¹ = local_geometry.J * Geometry.contravariant1(v, local_geometry)
+        Jv¹ =
+            form_jacobian(form, local_geometry) *
+            Geometry.contravariant1(v, local_geometry)
         for ii in 1:Nq
-            out[ii] += D[ii, i] * Jv¹
+            out[ii] += form_deriv_entry(form, D, ii, i) * Jv¹
         end
     end
     @inbounds for i in 1:Nq
         ij = CartesianIndex((i,))
         local_geometry = get_local_geometry(space, ij, slabidx)
-        out[i] *= local_geometry.invJ
+        out[i] = form_jacobian_rescale(form, local_geometry, out[i])
     end
     return Field(immutable_slab_data(out), space)
 end
 
 Base.@propagate_inbounds function apply_operator(
-    op::Divergence{(1, 2)},
+    op::Divergence{(1, 2), F},
     space,
     slabidx,
     arg,
-)
+) where {F}
+    form = F()
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -629,19 +752,23 @@ Base.@propagate_inbounds function apply_operator(
         ij = CartesianIndex((i, j))
         local_geometry = get_local_geometry(space, ij, slabidx)
         v = get_node(space, arg, ij, slabidx)
-        Jv¹ = local_geometry.J * Geometry.contravariant1(v, local_geometry)
+        Jv¹ =
+            form_jacobian(form, local_geometry) *
+            Geometry.contravariant1(v, local_geometry)
         for ii in 1:Nq
-            out[1, ii, j, 1] += D[ii, i] * Jv¹
+            out[1, ii, j, 1] += form_deriv_entry(form, D, ii, i) * Jv¹
         end
-        Jv² = local_geometry.J * Geometry.contravariant2(v, local_geometry)
+        Jv² =
+            form_jacobian(form, local_geometry) *
+            Geometry.contravariant2(v, local_geometry)
         for jj in 1:Nq
-            out[1, i, jj, 1] += D[jj, j] * Jv²
+            out[1, i, jj, 1] += form_deriv_entry(form, D, jj, j) * Jv²
         end
     end
     @inbounds for j in 1:Nq, i in 1:Nq
         ij = CartesianIndex((i, j))
         local_geometry = get_local_geometry(space, ij, slabidx)
-        out[1, i, j, 1] *= local_geometry.invJ
+        out[1, i, j, 1] = form_jacobian_rescale(form, local_geometry, out[1, i, j, 1])
     end
     return Field(immutable_slab_data(out), space)
 end
@@ -831,106 +958,6 @@ function apply_operator(op::SplitDivergence{(1, 2)}, space, slabidx, arg1, arg2)
     return Field(immutable_slab_data(out), space)
 end
 
-"""
-    wdiv = WeakDivergence()
-    wdiv.(u)
-
-Computes the "weak divergence" of a vector field `u`.
-
-This is defined as the scalar field ``\\theta \\in \\mathcal{V}_0`` such that
-for all ``\\phi\\in \\mathcal{V}_0``
-```math
-\\int_\\Omega \\phi \\theta \\, d \\Omega
-=
-- \\int_\\Omega (\\nabla \\phi) \\cdot u \\,d \\Omega
-```
-where ``\\mathcal{V}_0`` is the space of ``u``.
-
-This arises as the contribution of the volume integral after applying
-integration by parts to the weak form expression of the divergence
-```math
-\\int_\\Omega \\phi (\\nabla \\cdot u) \\, d \\Omega
-=
-- \\int_\\Omega (\\nabla \\phi) \\cdot u \\,d \\Omega
-+ \\oint_{\\partial \\Omega} \\phi (u \\cdot n) \\,d \\sigma
-```
-
-It can be written in matrix form as
-```math
-ϕ^\\top WJ θ = - \\sum_i (D_i ϕ)^\\top WJ u^i
-```
-which reduces to
-```math
-θ = -(WJ)^{-1} \\sum_i D_i^\\top WJ u^i
-```
-where
- - ``J`` is the diagonal Jacobian matrix
- - ``W`` is the diagonal matrix of quadrature weights
- - ``D_i`` is the derivative matrix along the ``i``th dimension
-"""
-struct WeakDivergence{I} <: SpectralElementOperator{I} end
-WeakDivergence() = WeakDivergence{()}()
-WeakDivergence{()}(space) = WeakDivergence{operator_axes(space)}()
-
-operator_return_eltype(::WeakDivergence{I}, ::Type{S}) where {I, S} =
-    Geometry.divergence_result_type(S)
-
-function apply_operator(op::WeakDivergence{(1,)}, space, slabidx, arg)
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    D = Quadratures.differentiation_matrix(FT, QS)
-    # allocate temp output
-    RT = operator_return_eltype(op, eltype(arg))
-    out = slab_data(RT, FT, Nq)
-    fill!(parent(out), zero(FT))
-    @inbounds for i in 1:Nq
-        ij = CartesianIndex((i,))
-        local_geometry = get_local_geometry(space, ij, slabidx)
-        v = get_node(space, arg, ij, slabidx)
-        WJv¹ = local_geometry.WJ * Geometry.contravariant1(v, local_geometry)
-        for ii in 1:Nq
-            out[ii] += D[i, ii] * WJv¹
-        end
-    end
-    @inbounds for i in 1:Nq
-        ij = CartesianIndex((i,))
-        local_geometry = get_local_geometry(space, ij, slabidx)
-        out[i] /= -local_geometry.WJ
-    end
-    return Field(immutable_slab_data(out), space)
-end
-
-function apply_operator(op::WeakDivergence{(1, 2)}, space, slabidx, arg)
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    D = Quadratures.differentiation_matrix(FT, QS)
-    RT = operator_return_eltype(op, eltype(arg))
-    out = slab_data(RT, FT, Nq, Nq)
-    fill!(parent(out), zero(FT))
-
-    @inbounds for j in 1:Nq, i in 1:Nq
-        ij = CartesianIndex((i, j))
-        local_geometry = get_local_geometry(space, ij, slabidx)
-        v = get_node(space, arg, ij, slabidx)
-        WJv¹ = local_geometry.WJ * Geometry.contravariant1(v, local_geometry)
-        for ii in 1:Nq
-            out[1, ii, j, 1] += D[i, ii] * WJv¹
-        end
-        WJv² = local_geometry.WJ * Geometry.contravariant2(v, local_geometry)
-        for jj in 1:Nq
-            out[1, i, jj, 1] += D[j, jj] * WJv²
-        end
-    end
-    @inbounds for j in 1:Nq, i in 1:Nq
-        ij = CartesianIndex((i, j))
-        local_geometry = get_local_geometry(space, ij, slabidx)
-        out[1, i, j, 1] /= -local_geometry.WJ
-    end
-    return Field(immutable_slab_data(out), space)
-end
-
 """
     grad = Gradient()
     grad.(f)
@@ -953,14 +980,22 @@ where ``D_i`` is the derivative matrix along the ``i``th dimension.
 ## References
 - [Taylor2010](@cite), equation 16
 """
-struct Gradient{I} <: SpectralElementOperator{I} end
-Gradient() = Gradient{()}()
-Gradient{()}(space) = Gradient{operator_axes(space)}()
+struct Gradient{I, F <: FormType} <: SpectralElementOperator{I} end
+Gradient() = Gradient{(), StrongForm}()
+Gradient{I}() where {I} = Gradient{I, StrongForm}()
+rebuild_operator(::Gradient{I, F}, space) where {I, F} =
+    Gradient{operator_axes(space), F}()
 
 operator_return_eltype(::Gradient{I}, ::Type{S}) where {I, S} =
     Geometry.gradient_result_type(Val(I), S)
 
-function apply_operator(op::Gradient{(1,)}, space, slabidx, arg)
+# The strong gradient is Dᵢ f, while the weak gradient is -W⁻¹ Dᵢᵀ (W f); see
+# form_deriv_entry and form_weighted_arg for the form-dependent factors. Only
+# the weak form needs the final W⁻¹ rescale, which stays inlined in each
+# apply_operator so that the strong form can skip that loop entirely.
+
+function apply_operator(op::Gradient{(1,), F}, space, slabidx, arg) where {F}
+    form = F()
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -971,21 +1006,32 @@ function apply_operator(op::Gradient{(1,)}, space, slabidx, arg)
     fill!(parent(out), zero(FT))
     @inbounds for i in 1:Nq
         ij = CartesianIndex((i,))
-        x = get_node(space, arg, ij, slabidx)
+        local_geometry = get_local_geometry(space, ij, slabidx)
+        x = form_weighted_arg(form, local_geometry, get_node(space, arg, ij, slabidx))
         for ii in 1:Nq
-            ∂f∂ξ = Geometry.Covariant1Vector(D[ii, i]) ⊗ x
+            ∂f∂ξ =
+                Geometry.Covariant1Vector(form_deriv_entry(form, D, ii, i)) ⊗ x
             out[ii] += ∂f∂ξ
         end
     end
+    if F === WeakForm
+        @inbounds for i in 1:Nq
+            ij = CartesianIndex((i,))
+            local_geometry = get_local_geometry(space, ij, slabidx)
+            W = local_geometry.WJ * local_geometry.invJ
+            out[i] /= W
+        end
+    end
     return Field(immutable_slab_data(out), space)
 end
 
 Base.@propagate_inbounds function apply_operator(
-    op::Gradient{(1, 2)},
+    op::Gradient{(1, 2), F},
     space,
     slabidx,
     arg,
-)
+) where {F}
+    form = F()
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -997,118 +1043,33 @@ Base.@propagate_inbounds function apply_operator(
 
     @inbounds for j in 1:Nq, i in 1:Nq
         ij = CartesianIndex((i, j))
-        x = get_node(space, arg, ij, slabidx)
+        local_geometry = get_local_geometry(space, ij, slabidx)
+        x = form_weighted_arg(form, local_geometry, get_node(space, arg, ij, slabidx))
         for ii in 1:Nq
-            ∂f∂ξ₁ = Geometry.Covariant12Vector(D[ii, i], zero(eltype(D))) ⊗ x
+            ∂f∂ξ₁ =
+                Geometry.Covariant12Vector(
+                    form_deriv_entry(form, D, ii, i),
+                    zero(eltype(D)),
+                ) ⊗ x
             out[1, ii, j, 1] += ∂f∂ξ₁
         end
         for jj in 1:Nq
-            ∂f∂ξ₂ = Geometry.Covariant12Vector(zero(eltype(D)), D[jj, j]) ⊗ x
+            ∂f∂ξ₂ =
+                Geometry.Covariant12Vector(
+                    zero(eltype(D)),
+                    form_deriv_entry(form, D, jj, j),
+                ) ⊗ x
             out[1, i, jj, 1] += ∂f∂ξ₂
         end
     end
-    return Field(immutable_slab_data(out), space)
-end
-
-"""
-    wgrad = WeakGradient()
-    wgrad.(f)
-
-Compute the "weak gradient" of `f` on each element.
-
-This is defined as the the vector field ``\\theta \\in \\mathcal{V}_0`` such
-that for all ``\\phi \\in \\mathcal{V}_0``
-```math
-\\int_\\Omega \\phi \\cdot \\theta \\, d \\Omega
-=
-- \\int_\\Omega (\\nabla \\cdot \\phi) f \\, d\\Omega
-```
-where ``\\mathcal{V}_0`` is the space of ``f``.
-
-This arises from the contribution of the volume integral after by applying
-integration by parts to the weak form expression of the gradient
-```math
-\\int_\\Omega \\phi \\cdot (\\nabla f) \\, d \\Omega
-=
-- \\int_\\Omega f (\\nabla \\cdot \\phi) \\, d\\Omega
-+ \\oint_{\\partial \\Omega} f (\\phi \\cdot n) \\, d \\sigma
-```
-
-In matrix form, this becomes
-```math
-{\\phi^i}^\\top W J \\theta_i = - ( J^{-1} D_i J \\phi^i )^\\top W J f
-```
-which reduces to
-```math
-\\theta_i = -W^{-1} D_i^\\top W f
-```
-where ``D_i`` is the derivative matrix along the ``i``th dimension.
-"""
-struct WeakGradient{I} <: SpectralElementOperator{I} end
-WeakGradient() = WeakGradient{()}()
-WeakGradient{()}(space) = WeakGradient{operator_axes(space)}()
-
-operator_return_eltype(::WeakGradient{I}, ::Type{S}) where {I, S} =
-    Geometry.gradient_result_type(Val(I), S)
-
-function apply_operator(op::WeakGradient{(1,)}, space, slabidx, arg)
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    D = Quadratures.differentiation_matrix(FT, QS)
-    # allocate temp output
-    RT = operator_return_eltype(op, eltype(arg))
-    out = slab_data(RT, FT, Nq)
-    fill!(parent(out), zero(FT))
-    @inbounds for i in 1:Nq
-        ij = CartesianIndex((i,))
-        local_geometry = get_local_geometry(space, ij, slabidx)
-        W = local_geometry.WJ * local_geometry.invJ
-        Wx = W * get_node(space, arg, ij, slabidx)
-        for ii in 1:Nq
-            Dᵀ₁Wf = Geometry.Covariant1Vector(D[i, ii]) ⊗ Wx
-            out[ii] -= Dᵀ₁Wf
-        end
-    end
-    @inbounds for i in 1:Nq
-        ij = CartesianIndex((i,))
-        local_geometry = get_local_geometry(space, ij, slabidx)
-        W = local_geometry.WJ * local_geometry.invJ
-        out[i] /= W
-    end
-    return Field(immutable_slab_data(out), space)
-end
-
-function apply_operator(op::WeakGradient{(1, 2)}, space, slabidx, arg)
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    D = Quadratures.differentiation_matrix(FT, QS)
-    # allocate temp output
-    RT = operator_return_eltype(op, eltype(arg))
-    out = slab_data(RT, FT, Nq, Nq)
-    fill!(parent(out), zero(FT))
-
-    @inbounds for j in 1:Nq, i in 1:Nq
-        ij = CartesianIndex((i, j))
-        local_geometry = get_local_geometry(space, ij, slabidx)
-        W = local_geometry.WJ * local_geometry.invJ
-        Wx = W * get_node(space, arg, ij, slabidx)
-        for ii in 1:Nq
-            Dᵀ₁Wf = Geometry.Covariant12Vector(D[i, ii], zero(eltype(D))) ⊗ Wx
-            out[1, ii, j, 1] -= Dᵀ₁Wf
-        end
-        for jj in 1:Nq
-            Dᵀ₂Wf = Geometry.Covariant12Vector(zero(eltype(D)), D[j, jj]) ⊗ Wx
-            out[1, i, jj, 1] -= Dᵀ₂Wf
+    if F === WeakForm
+        @inbounds for j in 1:Nq, i in 1:Nq
+            ij = CartesianIndex((i, j))
+            local_geometry = get_local_geometry(space, ij, slabidx)
+            W = local_geometry.WJ * local_geometry.invJ
+            out[1, i, j, 1] /= W
         end
     end
-    @inbounds for j in 1:Nq, i in 1:Nq
-        ij = CartesianIndex((i, j))
-        local_geometry = get_local_geometry(space, ij, slabidx)
-        W = local_geometry.WJ * local_geometry.invJ
-        out[1, i, j, 1] /= W
-    end
     return Field(immutable_slab_data(out), space)
 end
 
@@ -1155,14 +1116,21 @@ Note that unused dimensions will be dropped: e.g. the 2D curl of a
 ## References
 - [Taylor2010](@cite), equation 17
 """
-struct Curl{I} <: CurlSpectralElementOperator{I} end
-Curl() = Curl{()}()
-Curl{()}(space) = Curl{operator_axes(space)}()
+struct Curl{I, F <: FormType} <: CurlSpectralElementOperator{I} end
+Curl() = Curl{(), StrongForm}()
+Curl{I}() where {I} = Curl{I, StrongForm}()
+rebuild_operator(::Curl{I, F}, space) where {I, F} =
+    Curl{operator_axes(space), F}()
 
 operator_return_eltype(::Curl{I}, ::Type{S}) where {I, S} =
     Geometry.curl_result_type(Val(I), S)
 
-function apply_operator(op::Curl{(1,)}, space, slabidx, arg)
+# The strong curl is εⁱʲᵏ J⁻¹ Dⱼ uₖ, while the weak curl is
+# -εⁱʲᵏ (WJ)⁻¹ Dⱼᵀ (W uₖ); see form_deriv_entry, form_weighted_arg, and
+# form_jacobian_rescale for the form-dependent factors.
+
+function apply_operator(op::Curl{(1,), F}, space, slabidx, arg) where {F}
+    form = F()
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -1175,11 +1143,19 @@ function apply_operator(op::Curl{(1,)}, space, slabidx, arg)
         ij = CartesianIndex((i,))
         local_geometry = get_local_geometry(space, ij, slabidx)
         v = get_node(space, arg, ij, slabidx)
-        v₂ = Geometry.covariant2(v, local_geometry)
-        v₃ = Geometry.covariant3(v, local_geometry)
+        v₂ = form_weighted_arg(
+            form,
+            local_geometry,
+            Geometry.covariant2(v, local_geometry),
+        )
+        v₃ = form_weighted_arg(
+            form,
+            local_geometry,
+            Geometry.covariant3(v, local_geometry),
+        )
         for ii in 1:Nq
-            D₁v₂ = D[ii, i] * v₂
-            D₁v₃ = D[ii, i] * v₃
+            D₁v₂ = form_deriv_entry(form, D, ii, i) * v₂
+            D₁v₃ = form_deriv_entry(form, D, ii, i) * v₃
             out[ii] +=
                 Geometry.Contravariant123Vector(zero(FT), -D₁v₃, D₁v₂)
         end
@@ -1187,12 +1163,13 @@ function apply_operator(op::Curl{(1,)}, space, slabidx, arg)
     @inbounds for i in 1:Nq
         ij = CartesianIndex((i,))
         local_geometry = get_local_geometry(space, ij, slabidx)
-        out[i] *= local_geometry.invJ
+        out[i] = form_jacobian_rescale(form, local_geometry, out[i])
     end
     return Field(immutable_slab_data(out), space)
 end
 
-function apply_operator(op::Curl{(1, 2)}, space, slabidx, arg)
+function apply_operator(op::Curl{(1, 2), F}, space, slabidx, arg) where {F}
+    form = F()
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -1205,18 +1182,30 @@ function apply_operator(op::Curl{(1, 2)}, space, slabidx, arg)
         ij = CartesianIndex((i, j))
         local_geometry = get_local_geometry(space, ij, slabidx)
         v = get_node(space, arg, ij, slabidx)
-        v₁ = Geometry.covariant1(v, local_geometry)
-        v₂ = Geometry.covariant2(v, local_geometry)
-        v₃ = Geometry.covariant3(v, local_geometry)
+        v₁ = form_weighted_arg(
+            form,
+            local_geometry,
+            Geometry.covariant1(v, local_geometry),
+        )
+        v₂ = form_weighted_arg(
+            form,
+            local_geometry,
+            Geometry.covariant2(v, local_geometry),
+        )
+        v₃ = form_weighted_arg(
+            form,
+            local_geometry,
+            Geometry.covariant3(v, local_geometry),
+        )
         for ii in 1:Nq
-            D₁v₃ = D[ii, i] * v₃
-            D₁v₂ = D[ii, i] * v₂
+            D₁v₃ = form_deriv_entry(form, D, ii, i) * v₃
+            D₁v₂ = form_deriv_entry(form, D, ii, i) * v₂
             out[1, ii, j, 1] +=
                 Geometry.Contravariant123Vector(zero(FT), -D₁v₃, D₁v₂)
         end
         for jj in 1:Nq
-            D₂v₃ = D[jj, j] * v₃
-            D₂v₁ = D[jj, j] * v₁
+            D₂v₃ = form_deriv_entry(form, D, jj, j) * v₃
+            D₂v₁ = form_deriv_entry(form, D, jj, j) * v₁
             out[1, i, jj, 1] +=
                 Geometry.Contravariant123Vector(D₂v₃, zero(FT), -D₂v₁)
         end
@@ -1224,119 +1213,7 @@ function apply_operator(op::Curl{(1, 2)}, space, slabidx, arg)
     @inbounds for j in 1:Nq, i in 1:Nq
         ij = CartesianIndex((i, j))
         local_geometry = get_local_geometry(space, ij, slabidx)
-        out[1, i, j, 1] *= local_geometry.invJ
-    end
-    return Field(immutable_slab_data(out), space)
-end
-
-"""
-    wcurl = WeakCurl()
-    wcurl.(u)
-
-Computes the "weak curl" on each element of a covariant vector field `u`.
-
-Note: The vector field ``u`` needs to be excliclty converted to a `CovaraintVector`,
-as then the `WeakCurl` is independent of the local metric tensor.
-
-This is defined as the vector field ``\\theta \\in \\mathcal{V}_0`` such that
-for all ``\\phi \\in \\mathcal{V}_0``
-```math
-\\int_\\Omega \\phi \\cdot \\theta \\, d \\Omega
-=
-\\int_\\Omega (\\nabla \\times \\phi) \\cdot u \\,d \\Omega
-```
-where ``\\mathcal{V}_0`` is the space of ``f``.
-
-This arises from the contribution of the volume integral after by applying
-integration by parts to the weak form expression of the curl
-```math
-\\int_\\Omega \\phi \\cdot (\\nabla \\times u) \\,d\\Omega
-=
-\\int_\\Omega (\\nabla \\times \\phi) \\cdot u \\,d \\Omega
-- \\oint_{\\partial \\Omega} (\\phi \\times u) \\cdot n \\,d\\sigma
-```
-
-In matrix form, this becomes
-```math
-{\\phi_i}^\\top W J \\theta^i = (J^{-1} \\epsilon^{kji} D_j \\phi_i)^\\top W J u_k
-```
-which, by using the anti-symmetry of the Levi-Civita symbol, reduces to
-```math
-\\theta^i = - \\epsilon^{ijk} (WJ)^{-1} D_j^\\top W u_k
-```
-"""
-struct WeakCurl{I} <: CurlSpectralElementOperator{I} end
-WeakCurl() = WeakCurl{()}()
-WeakCurl{()}(space) = WeakCurl{operator_axes(space)}()
-
-operator_return_eltype(::WeakCurl{I}, ::Type{S}) where {I, S} =
-    Geometry.curl_result_type(Val(I), S)
-
-function apply_operator(op::WeakCurl{(1,)}, space, slabidx, arg)
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    D = Quadratures.differentiation_matrix(FT, QS)
-    # allocate temp output
-    RT = operator_return_eltype(op, eltype(arg))
-    out = slab_data(RT, FT, Nq)
-    fill!(parent(out), zero(FT))
-    @inbounds for i in 1:Nq
-        ij = CartesianIndex((i,))
-        local_geometry = get_local_geometry(space, ij, slabidx)
-        v = get_node(space, arg, ij, slabidx)
-        W = local_geometry.WJ * local_geometry.invJ
-        Wv₂ = W * Geometry.covariant2(v, local_geometry)
-        Wv₃ = W * Geometry.covariant3(v, local_geometry)
-        for ii in 1:Nq
-            Dᵀ₁Wv₂ = D[i, ii] * Wv₂
-            Dᵀ₁Wv₃ = D[i, ii] * Wv₃
-            out[ii] +=
-                Geometry.Contravariant123Vector(zero(FT), Dᵀ₁Wv₃, -Dᵀ₁Wv₂)
-        end
-    end
-    @inbounds for i in 1:Nq
-        ij = CartesianIndex((i,))
-        local_geometry = get_local_geometry(space, ij, slabidx)
-        out[i] /= local_geometry.WJ
-    end
-    return Field(immutable_slab_data(out), space)
-end
-
-function apply_operator(op::WeakCurl{(1, 2)}, space, slabidx, arg)
-    FT = Spaces.undertype(space)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
-    D = Quadratures.differentiation_matrix(FT, QS)
-    # allocate temp output
-    RT = operator_return_eltype(op, eltype(arg))
-    out = slab_data(RT, FT, Nq, Nq)
-    fill!(parent(out), zero(FT))
-    @inbounds for j in 1:Nq, i in 1:Nq
-        ij = CartesianIndex((i, j))
-        local_geometry = get_local_geometry(space, ij, slabidx)
-        v = get_node(space, arg, ij, slabidx)
-        W = local_geometry.WJ * local_geometry.invJ
-        Wv₁ = W * Geometry.covariant1(v, local_geometry)
-        Wv₂ = W * Geometry.covariant2(v, local_geometry)
-        Wv₃ = W * Geometry.covariant3(v, local_geometry)
-        for ii in 1:Nq
-            Dᵀ₁Wv₃ = D[i, ii] * Wv₃
-            Dᵀ₁Wv₂ = D[i, ii] * Wv₂
-            out[1, ii, j, 1] +=
-                Geometry.Contravariant123Vector(zero(FT), Dᵀ₁Wv₃, -Dᵀ₁Wv₂)
-        end
-        for jj in 1:Nq
-            Dᵀ₂Wv₃ = D[j, jj] * Wv₃
-            Dᵀ₂Wv₁ = D[j, jj] * Wv₁
-            out[1, i, jj, 1] +=
-                Geometry.Contravariant123Vector(-Dᵀ₂Wv₃, zero(FT), Dᵀ₂Wv₁)
-        end
-    end
-    @inbounds for j in 1:Nq, i in 1:Nq
-        ij = CartesianIndex((i, j))
-        local_geometry = get_local_geometry(space, ij, slabidx)
-        out[1, i, j, 1] /= local_geometry.WJ
+        out[1, i, j, 1] = form_jacobian_rescale(form, local_geometry, out[1, i, j, 1])
     end
     return Field(immutable_slab_data(out), space)
 end
diff --git a/test/Operators/spectralelement/unit_form_types.jl b/test/Operators/spectralelement/unit_form_types.jl
new file mode 100644
index 0000000000..92a7233ff8
--- /dev/null
+++ b/test/Operators/spectralelement/unit_form_types.jl
@@ -0,0 +1,99 @@
+#=
+julia --project
+using Revise; include(joinpath("test", "Operators", "spectralelement", "unit_form_types.jl"))
+=#
+using Test
+using ClimaComms
+import ClimaCore:
+    Domains, Meshes, Topologies, Spaces, Quadratures, Operators as O
+
+function test_space()
+    domain = Domains.SphereDomain(1.0)
+    mesh = Meshes.EquiangularCubedSphere(domain, 2)
+    topology = Topologies.Topology2D(
+        ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()),
+        mesh,
+    )
+    return Spaces.SpectralElementSpace2D(topology, Quadratures.GLL{3}())
+end
+
+# The strong and weak variants of an operator differ only in a FormType type
+# parameter, with the weak variants available under their original names.
+@testset "FormType parameter and weak-form aliases" begin
+    for (strong, weak) in
+        (
+        (O.Divergence, O.WeakDivergence),
+        (O.Gradient, O.WeakGradient),
+        (O.Curl, O.WeakCurl),
+    )
+        # The alias is the WeakForm variant, not a distinct type.
+        @test weak{()} === strong{(), O.WeakForm}
+        @test weak{(1, 2)} === strong{(1, 2), O.WeakForm}
+
+        # Constructors default to the strong form, and every spelling of a form
+        # produces the same singleton.
+        @test strong() === strong{(), O.StrongForm}()
+        @test strong{(1, 2)}() === strong{(1, 2), O.StrongForm}()
+        @test weak() === strong{(), O.WeakForm}()
+        @test weak{(1, 2)}() === strong{(1, 2), O.WeakForm}()
+
+        # Unifying the types makes the weak variants subtypes of the strong
+        # names, so code that must distinguish forms has to name the FormType.
+        @test weak() isa strong
+        @test !(weak() isa strong{(), O.StrongForm})
+        @test !(strong() isa strong{(), O.WeakForm})
+
+        # rebuild_operator resets the operator axes to those of the space while
+        # preserving the form, which is what Base.Broadcast.instantiate relies on.
+        space = test_space()
+        for (form, op) in ((O.StrongForm, strong()), (O.WeakForm, weak()))
+            @test O.rebuild_operator(op, space) === strong{(1, 2), form}()
+        end
+    end
+end
+
+# The form-dependent factors that let a single operator body serve both forms.
+@testset "form-dependent factors" begin
+    for FT in (Float32, Float64)
+        D = FT[1 2 3; 4 5 6; 7 8 9]
+        local_geometry = (; J = FT(2), invJ = FT(0.5), WJ = FT(3))
+        W = local_geometry.WJ * local_geometry.invJ
+        x = FT(7)
+
+        # The weak form transposes D and flips its sign, which is the only
+        # change integration by parts makes to the accumulation loop.
+        @test O.form_deriv_entry(O.StrongForm(), D, 1, 3) == D[1, 3]
+        @test O.form_deriv_entry(O.WeakForm(), D, 1, 3) == -D[3, 1]
+
+        # Only the weak form weights its argument, by W = WJ J⁻¹.
+        @test O.form_weighted_arg(O.StrongForm(), local_geometry, x) === x
+        @test O.form_weighted_arg(O.WeakForm(), local_geometry, x) == W * x
+
+        @test O.form_jacobian(O.StrongForm(), local_geometry) == local_geometry.J
+        @test O.form_jacobian(O.WeakForm(), local_geometry) == local_geometry.WJ
+
+        @test O.form_jacobian_rescale(O.StrongForm(), local_geometry, x) ==
+              x / local_geometry.J
+        @test O.form_jacobian_rescale(O.WeakForm(), local_geometry, x) ==
+              x / local_geometry.WJ
+
+        # Only the weak gradient divides out W; the strong form is the identity,
+        # so it must return its argument untouched rather than rescale by one.
+        @test O.form_weight_rescale(O.StrongForm(), local_geometry, x) === x
+        @test O.form_weight_rescale(O.WeakForm(), local_geometry, x) == x / W
+
+        # Unified operator bodies accumulate with form_deriv_entry and then
+        # apply the strong form's sign pattern to the result. Folding the weak
+        # form's sign flip into each term must give bitwise the same value as
+        # negating the completed transposed sum, or the unified bodies would
+        # change the weak form's rounding.
+        values = FT.([0.1, 0.2, 0.3])
+        weak_sum = O.form_deriv_entry(O.WeakForm(), D, 1, 1) * values[1]
+        transposed_sum = D[1, 1] * values[1]
+        for k in 2:3
+            weak_sum += O.form_deriv_entry(O.WeakForm(), D, 1, k) * values[k]
+            transposed_sum += D[k, 1] * values[k]
+        end
+        @test weak_sum === -transposed_sum
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index bad5d60f93..719361320e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -37,6 +37,7 @@ UnitTest("Sphere spaces"                           ,"Spaces/sphere.jl"),
 UnitTest("Fields"                                  ,"Fields/unit_field.jl"), # has benchmarks
 UnitTest("Reinstantiate broadcasted"               ,"Operators/unit_reinstantiate_bc.jl"),
 UnitTest("Placeholder Fields"                      ,"Operators/unit_common.jl"),
+UnitTest("Spectral elem - form types"              ,"Operators/spectralelement/unit_form_types.jl"),
 UnitTest("Spectral elem - rectilinear"             ,"Operators/spectralelement/rectilinear.jl"),
 UnitTest("Spectral elem - opt"                     ,"Operators/spectralelement/opt.jl"),
 UnitTest("Spectral elem - gradient tensor"         ,"Operators/spectralelement/covar_deriv_ops.jl"),