diff --git a/benchmarks/3d/se_kernels.jl b/benchmarks/3d/se_kernels.jl
index e618df6dac..c384be4efe 100644
--- a/benchmarks/3d/se_kernels.jl
+++ b/benchmarks/3d/se_kernels.jl
@@ -11,6 +11,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     DataLayouts,
     RecursiveApply
@@ -73,7 +74,7 @@ function initialize_mwe(device, ::Type{FT}) where {FT}
         horizontal_mesh,
         Topologies.spacefillingcurve(horizontal_mesh),
     )
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     h_space = Spaces.SpectralElementSpace2D(horizontal_topology, quad)
 
     # vertical space
diff --git a/benchmarks/3d/vector_laplacian.jl b/benchmarks/3d/vector_laplacian.jl
index 1a7179979c..9e1c09b484 100644
--- a/benchmarks/3d/vector_laplacian.jl
+++ b/benchmarks/3d/vector_laplacian.jl
@@ -1,13 +1,13 @@
 
 using ClimaComms
 using ClimaCore:
-    Geometry, Domains, Meshes, Topologies, Spaces, Fields, Operators
+    Geometry, Domains, Meshes, Topologies, Spaces, Fields, Operators, Quadratures
 using CUDA, BenchmarkTools
 
 hdomain = Domains.SphereDomain(6.37122e6)
 hmesh = Meshes.EquiangularCubedSphere(hdomain, 30)
 htopology = Topologies.Topology2D(hmesh)
-hspace = Spaces.SpectralElementSpace2D(htopology, Spaces.Quadratures.GLL{4}())
+hspace = Spaces.SpectralElementSpace2D(htopology, Quadratures.GLL{4}())
 
 vdomain = Domains.IntervalDomain(
     Geometry.ZPoint(0.0),
diff --git a/benchmarks/bickleyjet/bickleyjet_dg_reference.jl b/benchmarks/bickleyjet/bickleyjet_dg_reference.jl
index 50571a8b85..85eec66537 100644
--- a/benchmarks/bickleyjet/bickleyjet_dg_reference.jl
+++ b/benchmarks/bickleyjet/bickleyjet_dg_reference.jl
@@ -1,11 +1,12 @@
 using Base.Threads
 import ClimaCore.Spaces
+import ClimaCore.Quadratures
 using CUDA
 
 function spaceconfig(::Val{Nq}, ::Type{DA} = Array) where {Nq, DA}
-    quad = Spaces.Quadratures.GLL{Nq}()
-    ξ, W = Spaces.Quadratures.quadrature_points(Float64, quad)
-    D = Spaces.Quadratures.differentiation_matrix(Float64, quad)
+    quad = Quadratures.GLL{Nq}()
+    ξ, W = Quadratures.quadrature_points(Float64, quad)
+    D = Quadratures.differentiation_matrix(Float64, quad)
     return (DA(ξ), DA(W), DA(D))
 end
 
diff --git a/benchmarks/bickleyjet/core_vs_ref.jl b/benchmarks/bickleyjet/core_vs_ref.jl
index 934f783ed5..b6e8bb4745 100644
--- a/benchmarks/bickleyjet/core_vs_ref.jl
+++ b/benchmarks/bickleyjet/core_vs_ref.jl
@@ -22,7 +22,7 @@ for Nq in Nqs
     # setup core
     mesh = Meshes.RectilinearMesh(domain, n1, n2)
     grid_topology = Topologies.Topology2D(mesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
     y0 = init_state.(Fields.coordinate_field(space), Ref(parameters))
diff --git a/benchmarks/bickleyjet/run_ref_thread.jl b/benchmarks/bickleyjet/run_ref_thread.jl
index a9897b108b..c791fb96bf 100644
--- a/benchmarks/bickleyjet/run_ref_thread.jl
+++ b/benchmarks/bickleyjet/run_ref_thread.jl
@@ -10,7 +10,7 @@ Nq = 4
 
 mesh = Meshes.RectilinearMesh(domain, n1, n2)
 grid_topology = Topologies.Topology2D(mesh)
-quad = Spaces.Quadratures.GLL{Nq}()
+quad = Quadratures.GLL{Nq}()
 space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
 X = coordinates(Val(Nq), n1, n2)
diff --git a/docs/src/lib/ClimaCoreTempestRemap.md b/docs/src/lib/ClimaCoreTempestRemap.md
index ad73102b59..88399fc6bb 100644
--- a/docs/src/lib/ClimaCoreTempestRemap.md
+++ b/docs/src/lib/ClimaCoreTempestRemap.md
@@ -49,7 +49,7 @@ apply_remap
 The following example converts an OrdinaryDiffEq solution object `sol` to a netcdf file, and remaps it to an regular latitude-longitude (RLL) grid.
 
 ```julia
-using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces
+using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Quadratures
 using NCDatasets, ClimaCoreTempestRemap
 
 # sol is the integrator solution
@@ -58,7 +58,7 @@ using NCDatasets, ClimaCoreTempestRemap
 
 # the issue is that the Space types changed since this changed
 # we can reconstruct it by digging around a bit
-Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(cspace))
+Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(cspace))
 
 datafile_cc = "test.nc"
 NCDataset(datafile_cc, "c") do nc
@@ -115,7 +115,7 @@ remap_weights(
     meshfile_rll,
     meshfile_overlap;
     in_type = "cgll",
-    in_np = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(cspace)),
+    in_np = Quadratures.degrees_of_freedom(Spaces.quadrature_style(cspace)),
 )
 
 # apply remap
diff --git a/docs/tutorials/introduction.jl b/docs/tutorials/introduction.jl
index 460dc6fa3b..411166908d 100644
--- a/docs/tutorials/introduction.jl
+++ b/docs/tutorials/introduction.jl
@@ -84,7 +84,7 @@ column_face_space =
 # These nodes are chosen by a particular *quadrature rule*, which allows us to integrate functions over the domain. The only supported choice for now is a Gauss-Legendre-Lobatto rule.
 
 ## Gauss-Legendre-Lobatto quadrature with 4 nodes in each direction, so 16 in each element
-quad = ClimaCore.Spaces.Quadratures.GLL{4}()
+quad = ClimaCore.Quadratures.GLL{4}()
 rectangle_space =
     ClimaCore.Spaces.SpectralElementSpace2D(rectangle_topology, quad)
 #----------------------------------------------------------------------------
diff --git a/examples/bickleyjet/bickleyjet_cg.jl b/examples/bickleyjet/bickleyjet_cg.jl
index 6c4a40d016..622529a481 100644
--- a/examples/bickleyjet/bickleyjet_cg.jl
+++ b/examples/bickleyjet/bickleyjet_cg.jl
@@ -2,7 +2,7 @@ using ClimaComms
 using LinearAlgebra
 
 import ClimaCore:
-    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies
+    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Quadratures
 import ClimaCore.Geometry: ⊗
 
 using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
@@ -40,10 +40,10 @@ Nq = 4
 Nqh = 7
 mesh = Meshes.RectilinearMesh(domain, n1, n2)
 grid_topology = Topologies.Topology2D(context, mesh)
-quad = Spaces.Quadratures.GLL{Nq}()
+quad = Quadratures.GLL{Nq}()
 space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
-Iquad = Spaces.Quadratures.GLL{Nqh}()
+Iquad = Quadratures.GLL{Nqh}()
 Ispace = Spaces.SpectralElementSpace2D(grid_topology, Iquad)
 
 function init_state(coord, p)
diff --git a/examples/bickleyjet/bickleyjet_cg_invariant_hypervisc.jl b/examples/bickleyjet/bickleyjet_cg_invariant_hypervisc.jl
index 3fbaccc140..1f92e04aef 100644
--- a/examples/bickleyjet/bickleyjet_cg_invariant_hypervisc.jl
+++ b/examples/bickleyjet/bickleyjet_cg_invariant_hypervisc.jl
@@ -9,6 +9,7 @@ import ClimaCore:
     Operators,
     Spaces,
     Topologies,
+    Quadratures,
     DataLayouts
 using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
 
@@ -55,7 +56,7 @@ domain = Domains.RectangleDomain(
 )
 n1, n2 = 16, 16
 Nq = 4
-quad = Spaces.Quadratures.GLL{Nq}()
+quad = Quadratures.GLL{Nq}()
 mesh = Meshes.RectilinearMesh(domain, n1, n2)
 grid_topology = Topologies.Topology2D(context, mesh)
 if usempi
diff --git a/examples/bickleyjet/bickleyjet_cg_unsmesh.jl b/examples/bickleyjet/bickleyjet_cg_unsmesh.jl
index 783c99d7e7..4b3e197b1b 100644
--- a/examples/bickleyjet/bickleyjet_cg_unsmesh.jl
+++ b/examples/bickleyjet/bickleyjet_cg_unsmesh.jl
@@ -2,7 +2,7 @@ using ClimaComms
 using LinearAlgebra
 
 import ClimaCore:
-    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies
+    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Quadratures
 import ClimaCore.Geometry: ⊗
 
 using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
@@ -41,10 +41,10 @@ Nqh = 7
 
 mesh = Meshes.RectilinearMesh(domain, n1, n2)
 grid_topology = Topologies.Topology2D(context, mesh)
-quad = Spaces.Quadratures.GLL{Nq}()
+quad = Quadratures.GLL{Nq}()
 space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
-Iquad = Spaces.Quadratures.GLL{Nqh}()
+Iquad = Quadratures.GLL{Nqh}()
 Ispace = Spaces.SpectralElementSpace2D(grid_topology, Iquad)
 
 function init_state(coord, p)
diff --git a/examples/bickleyjet/bickleyjet_dg.jl b/examples/bickleyjet/bickleyjet_dg.jl
index 8d659aa161..e5e80436e5 100644
--- a/examples/bickleyjet/bickleyjet_dg.jl
+++ b/examples/bickleyjet/bickleyjet_dg.jl
@@ -9,6 +9,7 @@ import ClimaCore:
     Operators,
     RecursiveApply,
     Spaces,
+    Quadratures,
     Topologies
 import ClimaCore.Geometry: ⊗
 import ClimaCore.RecursiveApply: ⊞, rdiv, rmap
@@ -51,10 +52,10 @@ Nq = 4
 Nqh = 7
 mesh = Meshes.RectilinearMesh(domain, n1, n2)
 grid_topology = Topologies.Topology2D(context, mesh)
-quad = Spaces.Quadratures.GLL{Nq}()
+quad = Quadratures.GLL{Nq}()
 space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
-Iquad = Spaces.Quadratures.GLL{Nqh}()
+Iquad = Quadratures.GLL{Nqh}()
 Ispace = Spaces.SpectralElementSpace2D(grid_topology, Iquad)
 
 function init_state(coord, p)
@@ -206,7 +207,7 @@ function rhs!(dydt, y, (parameters, numflux), t)
     # 6. Solve for final result
     dydt_data = Fields.field_values(dydt)
     dydt_data .= RecursiveApply.rdiv.(dydt_data, space.local_geometry.WJ)
-    M = Spaces.Quadratures.cutoff_filter_matrix(
+    M = Quadratures.cutoff_filter_matrix(
         Float64,
         Spaces.quadrature_style(space),
         3,
diff --git a/examples/common_spaces.jl b/examples/common_spaces.jl
index 1148b7ccf3..81a93ca4ef 100644
--- a/examples/common_spaces.jl
+++ b/examples/common_spaces.jl
@@ -1,5 +1,5 @@
 using ClimaComms
-using ClimaCore: Geometry, Domains, Meshes, Topologies, Spaces
+using ClimaCore: Geometry, Domains, Meshes, Topologies, Spaces, Quadratures
 
 function periodic_line_mesh(; x_max, x_elem)
     domain = Domains.IntervalDomain(
@@ -36,7 +36,7 @@ function make_horizontal_space(
     npoly,
     context::ClimaComms.SingletonCommsContext,
 )
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     if mesh isa Meshes.AbstractMesh1D
         topology = Topologies.IntervalTopology(mesh)
         space = Spaces.SpectralElementSpace1D(topology, quad)
@@ -52,7 +52,7 @@ function make_horizontal_space(
     npoly,
     comms_ctx::ClimaComms.MPICommsContext,
 )
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     if mesh isa Meshes.AbstractMesh1D
         error("Distributed mode does not work with 1D horizontal spaces.")
     elseif mesh isa Meshes.AbstractMesh2D
diff --git a/examples/hybrid/box/bubble_3d_invariant_rhoe.jl b/examples/hybrid/box/bubble_3d_invariant_rhoe.jl
index 0af040705c..8559239df1 100644
--- a/examples/hybrid/box/bubble_3d_invariant_rhoe.jl
+++ b/examples/hybrid/box/bubble_3d_invariant_rhoe.jl
@@ -24,6 +24,7 @@ import ClimaCore:
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators
 using ClimaCorePlots, Plots
@@ -149,7 +150,7 @@ function hvspace_3D(
     )
     Nv = Meshes.nelements(vertmesh)
     Nf_center, Nf_face = 2, 1 #1 + 3 + 1
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzmesh = Meshes.RectilinearMesh(horzdomain, xyelem, xyelem)
     horztopology = Topologies.Topology2D(comms_ctx, horzmesh)
     horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
@@ -191,7 +192,7 @@ end
 function compute_κ₄(sim_params::SimulationParameters{FT}) where {FT}
     (; lxy, xyelem, npoly) = sim_params
     quad_points, _ =
-        Spaces.Quadratures.quadrature_points(FT, Quadratures.GLL{npoly + 1}())
+        Quadratures.quadrature_points(FT, Quadratures.GLL{npoly + 1}())
     Δx = (lxy / xyelem) * diff(quad_points)[1] / 2
     κ₄ = 1.0e6 * (Δx / lxy)^3.2
     return κ₄
diff --git a/examples/hybrid/box/bubble_3d_invariant_rhotheta.jl b/examples/hybrid/box/bubble_3d_invariant_rhotheta.jl
index 0abd187775..a5b3ce5d11 100644
--- a/examples/hybrid/box/bubble_3d_invariant_rhotheta.jl
+++ b/examples/hybrid/box/bubble_3d_invariant_rhotheta.jl
@@ -5,12 +5,12 @@ using LinearAlgebra
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators
 using ClimaCore.Geometry
@@ -56,7 +56,7 @@ function hvspace_3D(
     vertmesh = Meshes.IntervalMesh(zdomain, nelems = zelem)
     vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)
 
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
     hv_center_space =
diff --git a/examples/hybrid/box/bubble_3d_rhotheta.jl b/examples/hybrid/box/bubble_3d_rhotheta.jl
index c2cde998d7..df20085eee 100644
--- a/examples/hybrid/box/bubble_3d_rhotheta.jl
+++ b/examples/hybrid/box/bubble_3d_rhotheta.jl
@@ -5,12 +5,12 @@ using LinearAlgebra, StaticArrays
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators
 import ClimaCore.Geometry: ⊗
@@ -55,7 +55,7 @@ function hvspace_3D(
     vertmesh = Meshes.IntervalMesh(zdomain, nelems = zelem)
     vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)
 
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
     hv_center_space =
diff --git a/examples/hybrid/box/limiters_advection.jl b/examples/hybrid/box/limiters_advection.jl
index 16c71e2778..f9c38aa4de 100644
--- a/examples/hybrid/box/limiters_advection.jl
+++ b/examples/hybrid/box/limiters_advection.jl
@@ -16,6 +16,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     Limiters,
     slab
@@ -80,7 +81,7 @@ function hvspace_3D(
     z_topology = Topologies.IntervalTopology(context, z_mesh)
     z_space = Spaces.CenterFiniteDifferenceSpace(z_topology)
 
-    quad = Spaces.Quadratures.GLL{Nij}()
+    quad = Quadratures.GLL{Nij}()
     hspace = Spaces.SpectralElementSpace2D(horztopology, quad)
     cspace = Spaces.ExtrudedFiniteDifferenceSpace(hspace, z_space)
     fspace = Spaces.FaceExtrudedFiniteDifferenceSpace(cspace)
diff --git a/examples/hybrid/hybrid3dcs_dss.jl b/examples/hybrid/hybrid3dcs_dss.jl
index ecffa77ada..cc5d2a1245 100644
--- a/examples/hybrid/hybrid3dcs_dss.jl
+++ b/examples/hybrid/hybrid3dcs_dss.jl
@@ -12,6 +12,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     DataLayouts
 
@@ -69,7 +70,7 @@ function hybrid3dcubedsphere_dss_profiler(
     z_stretch_string = "uniform"
     horizontal_mesh = cubed_sphere_mesh(; radius = R, h_elem = h_elem)
 
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     h_topology = Topologies.Topology2D(
         comms_ctx,
         horizontal_mesh,
diff --git a/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl b/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl
index fcb6110105..030e6e3028 100644
--- a/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl
+++ b/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl
@@ -12,6 +12,7 @@ import ClimaCore:
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators
 using ClimaCore.Geometry
@@ -44,7 +45,7 @@ function hvspace_2D(
     horzmesh = Meshes.IntervalMesh(horzdomain, nelems = xelem)
     horztopology = Topologies.IntervalTopology(horzmesh)
 
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)
 
     hv_center_space =
diff --git a/examples/hybrid/plane/density_current_2d_flux_form.jl b/examples/hybrid/plane/density_current_2d_flux_form.jl
index d326095a64..800a37821a 100644
--- a/examples/hybrid/plane/density_current_2d_flux_form.jl
+++ b/examples/hybrid/plane/density_current_2d_flux_form.jl
@@ -4,12 +4,12 @@ using LinearAlgebra, StaticArrays
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators
 import ClimaCore.Geometry: ⊗
@@ -43,7 +43,7 @@ function hvspace_2D(
     horzmesh = Meshes.IntervalMesh(horzdomain; nelems = helem)
     horztopology = Topologies.IntervalTopology(horzmesh)
 
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)
 
     hv_center_space =
diff --git a/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl b/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl
index d9516d152d..65c7533ac9 100644
--- a/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl
+++ b/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl
@@ -6,12 +6,12 @@ using StaticArrays, IntervalSets, LinearAlgebra
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators
 using ClimaCore.Geometry
@@ -44,7 +44,7 @@ function hvspace_2D(
     horzmesh = Meshes.IntervalMesh(horzdomain, nelems = xelem)
     horztopology = Topologies.IntervalTopology(horzmesh)
 
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)
 
     hv_center_space =
diff --git a/examples/hybrid/plane/isothermal_channel.jl b/examples/hybrid/plane/isothermal_channel.jl
index 9171950281..5c987c37cf 100644
--- a/examples/hybrid/plane/isothermal_channel.jl
+++ b/examples/hybrid/plane/isothermal_channel.jl
@@ -4,12 +4,12 @@ using StaticArrays, IntervalSets, LinearAlgebra
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators,
     Hypsography
@@ -53,7 +53,7 @@ function hvspace_2D(
     )
     horzmesh = Meshes.IntervalMesh(horzdomain, nelems = xelem)
     horztopology = Topologies.IntervalTopology(horzmesh)
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)
     z_surface = warp_fn.(Fields.coordinate_field(horzspace))
     hv_face_space = Spaces.ExtrudedFiniteDifferenceSpace(
diff --git a/examples/hybrid/plane/topo_agnesi_nh.jl b/examples/hybrid/plane/topo_agnesi_nh.jl
index 7f1e029c04..8ece11136e 100644
--- a/examples/hybrid/plane/topo_agnesi_nh.jl
+++ b/examples/hybrid/plane/topo_agnesi_nh.jl
@@ -4,12 +4,11 @@ using StaticArrays, IntervalSets, LinearAlgebra
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
-    Spaces,
+    Quadratures,
     Fields,
     Operators,
     Hypsography
@@ -66,7 +65,7 @@ function hvspace_2D(
     )
     horzmesh = Meshes.IntervalMesh(horzdomain, nelems = xelem)
     horztopology = Topologies.IntervalTopology(horzmesh)
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)
 
     z_surface = warp_fn.(Fields.coordinate_field(horzspace))
diff --git a/examples/hybrid/plane/topo_schar_nh.jl b/examples/hybrid/plane/topo_schar_nh.jl
index 64c5368896..bbd9bd17d2 100644
--- a/examples/hybrid/plane/topo_schar_nh.jl
+++ b/examples/hybrid/plane/topo_schar_nh.jl
@@ -4,12 +4,12 @@ using StaticArrays, IntervalSets, LinearAlgebra
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators,
     Hypsography
@@ -89,7 +89,7 @@ function hvspace_2D(
     )
     horzmesh = Meshes.IntervalMesh(horzdomain, nelems = xelem)
     horztopology = Topologies.IntervalTopology(horzmesh)
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)
 
     z_surface = warp_fn.(Fields.coordinate_field(horzspace))
diff --git a/examples/hybrid/sphere/deformation_flow.jl b/examples/hybrid/sphere/deformation_flow.jl
index f56071e2ce..f9253e5995 100644
--- a/examples/hybrid/sphere/deformation_flow.jl
+++ b/examples/hybrid/sphere/deformation_flow.jl
@@ -4,7 +4,7 @@ using Test
 using Statistics: mean
 
 using ClimaCore:
-    Geometry, Domains, Meshes, Topologies, Spaces, Fields, Operators, Limiters
+    Geometry, Domains, Meshes, Topologies, Spaces, Fields, Operators, Limiters, Quadratures
 using ClimaTimeSteppers
 
 using Logging
@@ -200,7 +200,7 @@ function run_deformation_flow(use_limiter, fct_op)
     horz_domain = Domains.SphereDomain(R)
     horz_mesh = Meshes.EquiangularCubedSphere(horz_domain, helem)
     horz_topology = Topologies.Topology2D(context, horz_mesh)
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horz_space = Spaces.SpectralElementSpace2D(horz_topology, quad)
 
     cent_space =
diff --git a/examples/hybrid/sphere/hadley_circulation.jl b/examples/hybrid/sphere/hadley_circulation.jl
index 619038387b..3eebc63124 100644
--- a/examples/hybrid/sphere/hadley_circulation.jl
+++ b/examples/hybrid/sphere/hadley_circulation.jl
@@ -5,12 +5,12 @@ using LinearAlgebra
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators
 import ClimaCore.Utilities: half
@@ -72,7 +72,7 @@ function sphere_3D(
     horzdomain = Domains.SphereDomain(R)
     horzmesh = Meshes.EquiangularCubedSphere(horzdomain, helem)
     horztopology = Topologies.Topology2D(context, horzmesh)
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
     hv_center_space =
@@ -270,7 +270,7 @@ remap_weights(
     meshfile_rll,
     meshfile_overlap;
     in_type = "cgll",
-    in_np = Spaces.Quadratures.degrees_of_freedom(
+    in_np = Quadratures.degrees_of_freedom(
         Spaces.quadrature_style(hv_center_space),
     ),
 )
diff --git a/examples/hybrid/sphere/nonhydrostatic_gravity_wave.jl b/examples/hybrid/sphere/nonhydrostatic_gravity_wave.jl
index 8eb24ae410..02bf96e567 100644
--- a/examples/hybrid/sphere/nonhydrostatic_gravity_wave.jl
+++ b/examples/hybrid/sphere/nonhydrostatic_gravity_wave.jl
@@ -5,12 +5,12 @@ using LinearAlgebra
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators
 
@@ -63,7 +63,7 @@ function sphere_3D(
     horzdomain = Domains.SphereDomain(R)
     horzmesh = Meshes.EquiangularCubedSphere(horzdomain, helem)
     horztopology = Topologies.Topology2D(context, horzmesh)
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
     hv_center_space =
diff --git a/examples/hybrid/sphere/remap_pipeline.jl b/examples/hybrid/sphere/remap_pipeline.jl
index 1fffecc833..1fbc752f44 100644
--- a/examples/hybrid/sphere/remap_pipeline.jl
+++ b/examples/hybrid/sphere/remap_pipeline.jl
@@ -1,5 +1,5 @@
 import ClimaCore
-using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Operators
+using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Operators, Quadratures
 using NCDatasets
 using ClimaCoreTempestRemap
 
@@ -53,7 +53,7 @@ function remap2latlon(filein, nc_dir, nlat, nlon)
     # reconstruct space, obtain Nq from space
     cspace = axes(Y.c)
     hspace = Spaces.horizontal_space(cspace)
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(hspace))
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(hspace))
 
     # create a temporary dir for intermediate data
     remap_tmpdir = nc_dir * "remaptmp/"
diff --git a/examples/hybrid/sphere/solid_body_rotation_3d.jl b/examples/hybrid/sphere/solid_body_rotation_3d.jl
index 1aa6b72616..5089373b68 100644
--- a/examples/hybrid/sphere/solid_body_rotation_3d.jl
+++ b/examples/hybrid/sphere/solid_body_rotation_3d.jl
@@ -5,12 +5,12 @@ using LinearAlgebra
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators
 import ClimaCore.Utilities: half
@@ -62,7 +62,7 @@ function sphere_3D(
     horzdomain = Domains.SphereDomain(R)
     horzmesh = Meshes.EquiangularCubedSphere(horzdomain, helem)
     horztopology = Topologies.Topology2D(context, horzmesh)
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
     hv_center_space =
diff --git a/examples/plane/limiters_advection.jl b/examples/plane/limiters_advection.jl
index 0bf3a490f0..5fc432fb82 100644
--- a/examples/plane/limiters_advection.jl
+++ b/examples/plane/limiters_advection.jl
@@ -2,7 +2,7 @@ using ClimaComms
 using LinearAlgebra
 
 import ClimaCore:
-    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Limiters
+    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Limiters, Quadratures
 
 using OrdinaryDiffEq: ODEProblem, solve
 using ClimaTimeSteppers
@@ -140,7 +140,7 @@ Nq = 4
 for (k, ne) in enumerate(ne_seq)
     mesh = Meshes.RectilinearMesh(domain, ne, ne)
     grid_topology = Topologies.Topology2D(context, mesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
     # Initialize state
diff --git a/examples/sphere/limiters_advection.jl b/examples/sphere/limiters_advection.jl
index d5bc54e0d0..51c0c59eef 100644
--- a/examples/sphere/limiters_advection.jl
+++ b/examples/sphere/limiters_advection.jl
@@ -2,7 +2,7 @@ using ClimaComms
 using LinearAlgebra
 
 import ClimaCore:
-    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Limiters
+    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Limiters, Quadratures
 
 using OrdinaryDiffEq, Test
 
@@ -102,7 +102,7 @@ for (k, ne) in enumerate(ne_seq)
     domain = Domains.SphereDomain(R)
     mesh = Meshes.EquiangularCubedSphere(domain, ne)
     grid_topology = Topologies.Topology2D(context, mesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space =
         Spaces.SpectralElementSpace2D(grid_topology, quad; enable_bubble = true)
 
diff --git a/examples/sphere/shallow_water.jl b/examples/sphere/shallow_water.jl
index c6a8a030c3..99b8ecca07 100644
--- a/examples/sphere/shallow_water.jl
+++ b/examples/sphere/shallow_water.jl
@@ -12,6 +12,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     DataLayouts
 
@@ -539,7 +540,7 @@ function shallow_water_driver(ARGS, ::Type{FT}) where {FT}
 
     domain = Domains.SphereDomain(test.params.R)
     mesh = Meshes.EquiangularCubedSphere(domain, ne)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     grid_topology = Topologies.Topology2D(context, mesh)
     if usempi
         global_grid_topology =
diff --git a/examples/sphere/shallow_water_cuda.jl b/examples/sphere/shallow_water_cuda.jl
index 407402e6d7..36b3f6c259 100644
--- a/examples/sphere/shallow_water_cuda.jl
+++ b/examples/sphere/shallow_water_cuda.jl
@@ -18,6 +18,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     DataLayouts,
     InputOutput
@@ -523,7 +524,7 @@ function shallow_water_driver_cuda(ARGS, ::Type{FT}) where {FT}
 
     domain = Domains.SphereDomain(test.params.R)
     mesh = Meshes.EquiangularCubedSphere(domain, ne)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     grid_topology = Topologies.Topology2D(context, mesh)
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
     f = set_coriolis_parameter(space, test)
diff --git a/examples/sphere/shallow_water_dss.jl b/examples/sphere/shallow_water_dss.jl
index 0a6ef86843..7ca6a690f4 100644
--- a/examples/sphere/shallow_water_dss.jl
+++ b/examples/sphere/shallow_water_dss.jl
@@ -13,6 +13,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     DataLayouts
 
@@ -66,7 +67,7 @@ function shallow_water_dss_profiler(::Type{FT}, npoly) where {FT}
     R = FT(6.37122e6) # radius of earth
     domain = Domains.SphereDomain(R)
     mesh = Meshes.EquiangularCubedSphere(domain, ne)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     grid_topology =
         Topologies.Topology2D(context, mesh, Topologies.spacefillingcurve(mesh))
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
diff --git a/examples/sphere/solidbody.jl b/examples/sphere/solidbody.jl
index a7edc694d4..c3c255ce32 100644
--- a/examples/sphere/solidbody.jl
+++ b/examples/sphere/solidbody.jl
@@ -2,7 +2,7 @@ using ClimaComms
 using LinearAlgebra
 
 import ClimaCore:
-    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies
+    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Quadratures
 
 using OrdinaryDiffEq: ODEProblem, solve, SSPRK33
 
@@ -80,7 +80,7 @@ for (k, ne) in enumerate(ne_seq)
     domain = Domains.SphereDomain(R)
     mesh = Meshes.EquiangularCubedSphere(domain, ne)
     grid_topology = Topologies.Topology2D(context, mesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
     coords = Fields.coordinate_field(space)
diff --git a/lib/ClimaCoreMakie/examples/sphere.jl b/lib/ClimaCoreMakie/examples/sphere.jl
index 7683a20c21..18630b5c2f 100644
--- a/lib/ClimaCoreMakie/examples/sphere.jl
+++ b/lib/ClimaCoreMakie/examples/sphere.jl
@@ -7,7 +7,7 @@ R = 6.37122e6
 domain = ClimaCore.Domains.SphereDomain(R)
 mesh = ClimaCore.Meshes.EquiangularCubedSphere(domain, 32)
 grid_topology = ClimaCore.Topologies.Topology2D(mesh)
-quad = ClimaCore.Spaces.Quadratures.GLL{5}()
+quad = ClimaCore.Quadratures.GLL{5}()
 space = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)
 coords = ClimaCore.Fields.coordinate_field(space)
 
diff --git a/lib/ClimaCoreMakie/test/runtests.jl b/lib/ClimaCoreMakie/test/runtests.jl
index baa8601856..6301b879d8 100644
--- a/lib/ClimaCoreMakie/test/runtests.jl
+++ b/lib/ClimaCoreMakie/test/runtests.jl
@@ -18,7 +18,7 @@ OUTPUT_DIR = mkpath(get(ENV, "CI_OUTPUT_DIR", tempname()))
         ClimaComms.SingletonCommsContext(),
         mesh,
     )
-    quad = ClimaCore.Spaces.Quadratures.GLL{5}()
+    quad = ClimaCore.Quadratures.GLL{5}()
     space = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)
     coords = ClimaCore.Fields.coordinate_field(space)
 
@@ -56,8 +56,8 @@ end
         ClimaComms.SingletonCommsContext(),
         mesh,
     )
-    #quad = ClimaCore.Spaces.Quadratures.GLL{Nq}()
-    quad = ClimaCore.Spaces.Quadratures.ClosedUniform{Nq + 1}()
+    #quad = ClimaCore.Quadratures.GLL{Nq}()
+    quad = ClimaCore.Quadratures.ClosedUniform{Nq + 1}()
     space = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)
     coords = ClimaCore.Fields.coordinate_field(space)
 
@@ -82,7 +82,7 @@ end
     )
     mesh = ClimaCore.Meshes.IntervalMesh(domain, nelems = 32)
     topology = ClimaCore.Topologies.IntervalTopology(mesh)
-    quad = ClimaCore.Spaces.Quadratures.GLL{5}()
+    quad = ClimaCore.Quadratures.GLL{5}()
     horz_space = ClimaCore.Spaces.SpectralElementSpace1D(topology, quad)
     horz_coords = ClimaCore.Fields.coordinate_field(horz_space)
 
diff --git a/lib/ClimaCorePlots/src/ClimaCorePlots.jl b/lib/ClimaCorePlots/src/ClimaCorePlots.jl
index 989bfab9e6..2e1e228b51 100644
--- a/lib/ClimaCorePlots/src/ClimaCorePlots.jl
+++ b/lib/ClimaCorePlots/src/ClimaCorePlots.jl
@@ -10,6 +10,7 @@ import ClimaCore:
     Geometry,
     Meshes,
     Operators,
+    Quadratures,
     Topologies,
     Spaces
 
@@ -26,11 +27,11 @@ RecipesBase.@recipe function f(
     nelems = Topologies.nlocalelems(space)
     QS = Spaces.quadrature_style(space)
     quad_name = Base.typename(typeof(QS)).name
-    Nq = Spaces.Quadratures.degrees_of_freedom(QS)
+    Nq = Quadratures.degrees_of_freedom(QS)
     Nu = max(interpolate, Nq)
     coord_field = Fields.coordinate_field(space)
 
-    lagrange_quad = Spaces.Quadratures.ClosedUniform{Nu}()
+    lagrange_quad = Quadratures.ClosedUniform{Nu}()
     dataspace =
         Spaces.SpectralElementSpace1D(Spaces.topology(space), lagrange_quad)
     interp = Operators.Interpolate(dataspace)
@@ -64,7 +65,7 @@ end
 RecipesBase.@recipe function f(space::Spaces.RectilinearSpectralElementSpace2D;)
     quad = Spaces.quadrature_style(space)
     quad_name = Base.typename(typeof(quad)).name
-    dof = Spaces.Quadratures.degrees_of_freedom(quad)
+    dof = Quadratures.degrees_of_freedom(quad)
 
     topology = Spaces.topology(space)
     mesh = topology.mesh
@@ -103,7 +104,7 @@ RecipesBase.@recipe function f(space::Spaces.ExtrudedFiniteDifferenceSpace)
 
     quad = Spaces.quadrature_style(hspace)
     quad_name = Base.typename(typeof(quad)).name
-    dof = Spaces.Quadratures.degrees_of_freedom(quad)
+    dof = Quadratures.degrees_of_freedom(quad)
 
     coord_symbols = propertynames(coord_field)
     hcoord = vec(parent(coord_field)[:, :, 1, :])
@@ -207,7 +208,7 @@ function _slice_triplot(field, hinterpolate, ncolors)
 
     if hinterpolate ≥ 1
         Nu = hinterpolate
-        uquad = Spaces.Quadratures.ClosedUniform{Nu}()
+        uquad = Quadratures.ClosedUniform{Nu}()
         M_hcoord = Operators.matrix_interpolate(hcoord_field, uquad)
         M_vcoord = Operators.matrix_interpolate(vcoord_field, uquad)
         M_data = Operators.matrix_interpolate(field, uquad)
@@ -402,8 +403,8 @@ function _unfolded_pannel_matrix(field, interpolate)
     panel_size = mesh.ne
 
     quad_from = Spaces.quadrature_style(space)
-    quad_to = Spaces.Quadratures.Uniform{interpolate}()
-    Imat = Spaces.Quadratures.interpolation_matrix(FT, quad_to, quad_from)
+    quad_to = Quadratures.Uniform{interpolate}()
+    Imat = Quadratures.interpolation_matrix(FT, quad_to, quad_from)
 
     dof = interpolate
 
@@ -481,7 +482,7 @@ RecipesBase.@recipe function f(
     space = axes(field)
     nelem = Topologies.nlocalelems(Spaces.topology(space))
     quad_from = Spaces.quadrature_style(space)
-    dof_in = Spaces.Quadratures.degrees_of_freedom(quad_from)
+    dof_in = Quadratures.degrees_of_freedom(quad_from)
     quad_from_name = Base.typename(typeof(quad_from)).name
 
     # set the plot attributes
@@ -516,7 +517,7 @@ RecipesBase.@recipe function f(
     nlevel = Spaces.nlevels(space)
     nelem = Topologies.nlocalelems(Spaces.topology(space))
     quad_from = Spaces.quadrature_style(space)
-    dof_in = Spaces.Quadratures.degrees_of_freedom(quad_from)
+    dof_in = Quadratures.degrees_of_freedom(quad_from)
     quad_from_name = Base.typename(typeof(quad_from)).name
 
     # set the plot attributes
diff --git a/lib/ClimaCorePlots/test/runtests.jl b/lib/ClimaCorePlots/test/runtests.jl
index 8127a0f09f..53eed69487 100644
--- a/lib/ClimaCorePlots/test/runtests.jl
+++ b/lib/ClimaCorePlots/test/runtests.jl
@@ -17,7 +17,7 @@ OUTPUT_DIR = mkpath(get(ENV, "CI_OUTPUT_DIR", tempname()))
     )
     mesh = ClimaCore.Meshes.IntervalMesh(domain; nelems = 5)
     grid_topology = ClimaCore.Topologies.IntervalTopology(mesh)
-    quad = ClimaCore.Spaces.Quadratures.GLL{5}()
+    quad = ClimaCore.Quadratures.GLL{5}()
     space = ClimaCore.Spaces.SpectralElementSpace1D(grid_topology, quad)
     coords = ClimaCore.Fields.coordinate_field(space)
 
@@ -40,7 +40,7 @@ end
         ClimaComms.SingletonCommsContext(),
         mesh,
     )
-    quad = ClimaCore.Spaces.Quadratures.GLL{5}()
+    quad = ClimaCore.Quadratures.GLL{5}()
     space = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)
     coords = ClimaCore.Fields.coordinate_field(space)
 
@@ -68,7 +68,7 @@ end
         mesh,
         ClimaCore.Topologies.spacefillingcurve(mesh),
     )
-    quad = ClimaCore.Spaces.Quadratures.GLL{5}()
+    quad = ClimaCore.Quadratures.GLL{5}()
     space = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)
     coords = ClimaCore.Fields.coordinate_field(space)
 
@@ -102,7 +102,7 @@ end
         ClimaComms.SingletonCommsContext(),
         horz_mesh,
     )
-    quad = ClimaCore.Spaces.Quadratures.GLL{5}()
+    quad = ClimaCore.Quadratures.GLL{5}()
     horz_space =
         ClimaCore.Spaces.SpectralElementSpace2D(horz_grid_topology, quad)
 
@@ -166,7 +166,7 @@ end
         mesh,
     )
 
-    quad = ClimaCore.Spaces.Quadratures.ClosedUniform{Nq + 1}()
+    quad = ClimaCore.Quadratures.ClosedUniform{Nq + 1}()
     space = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)
     coords = ClimaCore.Fields.coordinate_field(space)
 
@@ -210,7 +210,7 @@ end
     horzmesh = ClimaCore.Meshes.IntervalMesh(horzdomain; nelems = helem)
     horztopology = ClimaCore.Topologies.IntervalTopology(horzmesh)
 
-    quad = ClimaCore.Spaces.Quadratures.GLL{npoly + 1}()
+    quad = ClimaCore.Quadratures.GLL{npoly + 1}()
     horzspace = ClimaCore.Spaces.SpectralElementSpace1D(horztopology, quad)
 
     hv_center_space = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(
@@ -269,7 +269,7 @@ end
         horzmesh,
     )
 
-    quad = ClimaCore.Spaces.Quadratures.GLL{npoly + 1}()
+    quad = ClimaCore.Quadratures.GLL{npoly + 1}()
     horzspace = ClimaCore.Spaces.SpectralElementSpace2D(horztopology, quad)
 
     hv_center_space = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(
diff --git a/lib/ClimaCoreSpectra/examples/gcm_remap_and_spectra.jl b/lib/ClimaCoreSpectra/examples/gcm_remap_and_spectra.jl
index 26c736df46..84de1efc5c 100644
--- a/lib/ClimaCoreSpectra/examples/gcm_remap_and_spectra.jl
+++ b/lib/ClimaCoreSpectra/examples/gcm_remap_and_spectra.jl
@@ -1,12 +1,12 @@
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators
 
@@ -26,7 +26,7 @@ Nq = 5
 domain = Domains.SphereDomain(R)
 mesh = Meshes.EquiangularCubedSphere(domain, ne)
 topology = Topologies.Topology2D(mesh)
-quad = Spaces.Quadratures.GLL{Nq}()
+quad = Quadratures.GLL{Nq}()
 space = Spaces.SpectralElementSpace2D(topology, quad)
 coords = Fields.coordinate_field(space)
 
diff --git a/lib/ClimaCoreTempestRemap/src/ClimaCoreTempestRemap.jl b/lib/ClimaCoreTempestRemap/src/ClimaCoreTempestRemap.jl
index aed902e5dd..dc1683ce3a 100644
--- a/lib/ClimaCoreTempestRemap/src/ClimaCoreTempestRemap.jl
+++ b/lib/ClimaCoreTempestRemap/src/ClimaCoreTempestRemap.jl
@@ -6,7 +6,7 @@ export def_time_coord, def_space_coord
 
 using ClimaComms
 import ClimaCore
-using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Fields
+using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Fields, Quadratures
 
 using NCDatasets, Dates, PkgVersion, LinearAlgebra, TempestRemap_jll
 
diff --git a/lib/ClimaCoreTempestRemap/src/onlineremap.jl b/lib/ClimaCoreTempestRemap/src/onlineremap.jl
index 14cb831b00..7414b68a24 100644
--- a/lib/ClimaCoreTempestRemap/src/onlineremap.jl
+++ b/lib/ClimaCoreTempestRemap/src/onlineremap.jl
@@ -178,11 +178,11 @@ function generate_map(
             meshfile_target,
             meshfile_overlap;
             in_type = in_type,
-            in_np = Spaces.Quadratures.degrees_of_freedom(
+            in_np = Quadratures.degrees_of_freedom(
                 Spaces.quadrature_style(source_space),
             ),
             out_type = out_type,
-            out_np = Spaces.Quadratures.degrees_of_freedom(
+            out_np = Quadratures.degrees_of_freedom(
                 Spaces.quadrature_style(target_space),
             ),
         )
diff --git a/lib/ClimaCoreTempestRemap/test/exodus.jl b/lib/ClimaCoreTempestRemap/test/exodus.jl
index 440eabeeca..a18d6d1126 100644
--- a/lib/ClimaCoreTempestRemap/test/exodus.jl
+++ b/lib/ClimaCoreTempestRemap/test/exodus.jl
@@ -1,6 +1,6 @@
 import ClimaCore
 using ClimaComms
-using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces
+using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Quadratures
 using NCDatasets
 using TempestRemap_jll
 using Test
diff --git a/lib/ClimaCoreTempestRemap/test/mpi_tests/exodus.jl b/lib/ClimaCoreTempestRemap/test/mpi_tests/exodus.jl
index bf14603538..3b16753d03 100644
--- a/lib/ClimaCoreTempestRemap/test/mpi_tests/exodus.jl
+++ b/lib/ClimaCoreTempestRemap/test/mpi_tests/exodus.jl
@@ -1,6 +1,6 @@
 import ClimaCore
 using ClimaComms
-using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces
+using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Quadratures
 using NCDatasets
 using TempestRemap_jll
 using Test
diff --git a/lib/ClimaCoreTempestRemap/test/mpi_tests/online_remap.jl b/lib/ClimaCoreTempestRemap/test/mpi_tests/online_remap.jl
index 73c6e23c1f..38d50cff8d 100644
--- a/lib/ClimaCoreTempestRemap/test/mpi_tests/online_remap.jl
+++ b/lib/ClimaCoreTempestRemap/test/mpi_tests/online_remap.jl
@@ -1,7 +1,7 @@
 import ClimaCore
 import ClimaCoreTempestRemap as CCTR
 using ClimaComms
-using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Fields
+using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Fields, Quadratures
 using Test
 # use these packages for manual inspection of solutions
 # using Plots
@@ -28,7 +28,7 @@ using Test
     topology_i_distr = Topologies.Topology2D(comms_ctx, mesh_i)
 
     nq_i = 3 # polynomial order for SE discretization
-    quad_i = Spaces.Quadratures.GLL{nq_i}()
+    quad_i = Quadratures.GLL{nq_i}()
     space_i_distr = Spaces.SpectralElementSpace2D(topology_i_distr, quad_i)
 
     topology_i_singleton = Topologies.Topology2D(comms_ctx_singleton, mesh_i)
@@ -41,7 +41,7 @@ using Test
     topology_o_distr = Topologies.Topology2D(comms_ctx, mesh_o)
 
     nq_o = 4
-    quad_o = Spaces.Quadratures.GLL{nq_o}()
+    quad_o = Quadratures.GLL{nq_o}()
     space_o_distr = Spaces.SpectralElementSpace2D(topology_o_distr, quad_o)
 
     topology_o_singleton = Topologies.Topology2D(comms_ctx_singleton, mesh_o)
diff --git a/lib/ClimaCoreTempestRemap/test/netcdf.jl b/lib/ClimaCoreTempestRemap/test/netcdf.jl
index 7385b91806..a71b18b781 100644
--- a/lib/ClimaCoreTempestRemap/test/netcdf.jl
+++ b/lib/ClimaCoreTempestRemap/test/netcdf.jl
@@ -1,7 +1,7 @@
 import ClimaCore
 using ClimaComms
 using ClimaCore:
-    Geometry, Meshes, Domains, Topologies, Spaces, Fields, Hypsography
+    Geometry, Meshes, Domains, Topologies, Spaces, Fields, Hypsography, Quadratures
 using NCDatasets
 using TempestRemap_jll
 using Test
@@ -17,7 +17,7 @@ OUTPUT_DIR = mkpath(get(ENV, "CI_OUTPUT_DIR", tempname()))
     domain = Domains.SphereDomain(R)
     mesh = Meshes.EquiangularCubedSphere(domain, ne)
     topology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), mesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(topology, quad)
     coords = Fields.coordinate_field(space)
 
@@ -85,7 +85,7 @@ end
     hdomain = Domains.SphereDomain(R)
     hmesh = Meshes.EquiangularCubedSphere(hdomain, ne)
     htopology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), hmesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     hspace = Spaces.SpectralElementSpace2D(htopology, quad)
 
     vdomain = Domains.IntervalDomain(
@@ -168,7 +168,7 @@ end
     hdomain = Domains.SphereDomain(R)
     hmesh = Meshes.EquiangularCubedSphere(hdomain, ne)
     htopology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), hmesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     hspace = Spaces.SpectralElementSpace2D(htopology, quad)
 
     vdomain = Domains.IntervalDomain(
@@ -303,7 +303,7 @@ end
     hdomain = Domains.SphereDomain(R)
     hmesh = Meshes.EquiangularCubedSphere(hdomain, ne)
     htopology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), hmesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     hspace = Spaces.SpectralElementSpace2D(htopology, quad)
 
     vdomain = Domains.IntervalDomain(
@@ -411,7 +411,7 @@ end
     domain = Domains.RectangleDomain(x_domain, y_domain)
     hmesh = Meshes.RectilinearMesh(domain, x_elem, y_elem)
 
-    quad = Spaces.Quadratures.GL{1}()
+    quad = Quadratures.GL{1}()
     htopology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), hmesh)
     hspace = Spaces.SpectralElementSpace2D(htopology, quad)
 
diff --git a/lib/ClimaCoreTempestRemap/test/online_remap.jl b/lib/ClimaCoreTempestRemap/test/online_remap.jl
index 193a4618e6..e9a1a99ff9 100644
--- a/lib/ClimaCoreTempestRemap/test/online_remap.jl
+++ b/lib/ClimaCoreTempestRemap/test/online_remap.jl
@@ -1,6 +1,6 @@
 import ClimaCore
 using ClimaComms
-using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Fields
+using ClimaCore: Geometry, Meshes, Domains, Topologies, Spaces, Fields, Quadratures
 using NCDatasets
 using Test
 using ClimaCoreTempestRemap
@@ -61,7 +61,7 @@ end
     )
     space_i = Spaces.SpectralElementSpace2D(
         topology_i,
-        Spaces.Quadratures.GLL{nq_i}(),
+        Quadratures.GLL{nq_i}(),
     )
     coords_i = Fields.coordinate_field(space_i)
 
@@ -73,7 +73,7 @@ end
     )
     space_o = Spaces.SpectralElementSpace2D(
         topology_o,
-        Spaces.Quadratures.GLL{nq_o}(),
+        Quadratures.GLL{nq_o}(),
     )
     coords_o = Fields.coordinate_field(space_o)
 
diff --git a/lib/ClimaCoreVTK/src/ClimaCoreVTK.jl b/lib/ClimaCoreVTK/src/ClimaCoreVTK.jl
index e03b90e7a0..8cee98d072 100644
--- a/lib/ClimaCoreVTK/src/ClimaCoreVTK.jl
+++ b/lib/ClimaCoreVTK/src/ClimaCoreVTK.jl
@@ -6,7 +6,7 @@ module ClimaCoreVTK
 export writevtk, writepvd
 
 using WriteVTK
-import ClimaCore: Fields, Geometry, Spaces, Topologies, Operators
+import ClimaCore: Fields, Geometry, Spaces, Topologies, Operators, Quadratures
 
 include("space.jl")
 include("addfield.jl")
diff --git a/lib/ClimaCoreVTK/src/space.jl b/lib/ClimaCoreVTK/src/space.jl
index e5ca6f25cb..5547bd953f 100644
--- a/lib/ClimaCoreVTK/src/space.jl
+++ b/lib/ClimaCoreVTK/src/space.jl
@@ -104,7 +104,7 @@ an unstuctured mesh of linear cells, suitable for passing to
 `vtk_grid`.
 """
 function vtk_cells_linear(gridspace::Spaces.SpectralElementSpace2D)
-    Nq = Spaces.Quadratures.degrees_of_freedom(
+    Nq = Quadratures.degrees_of_freedom(
         Spaces.quadrature_style(gridspace),
     )
     Nh = Topologies.nlocalelems(gridspace)
@@ -125,7 +125,7 @@ end
 
 function vtk_cells_linear(gridspace::Spaces.FaceExtrudedFiniteDifferenceSpace2D)
     hspace = Spaces.horizontal_space(gridspace)
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(hspace))
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(hspace))
     Nh = Topologies.nlocalelems(hspace)
     Nv = Spaces.nlevels(gridspace)
     ind = LinearIndices((1:Nv, 1:Nq, 1:Nh)) # assume VIFH
@@ -144,7 +144,7 @@ function vtk_cells_linear(gridspace::Spaces.FaceExtrudedFiniteDifferenceSpace2D)
 end
 function vtk_cells_linear(gridspace::Spaces.FaceExtrudedFiniteDifferenceSpace3D)
     hspace = Spaces.horizontal_space(gridspace)
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(hspace))
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(hspace))
     Nh = Topologies.nlocalelems(hspace)
     Nv = Spaces.nlevels(gridspace)
     ind = LinearIndices((1:Nv, 1:Nq, 1:Nq, 1:Nh)) # assumes VIJFH
@@ -179,19 +179,19 @@ This generally does two things:
  - Modifies the vertical space to be on the faces.
 """
 function vtk_grid_space(space::Spaces.SpectralElementSpace1D)
-    if Spaces.quadrature_style(space) isa Spaces.Quadratures.ClosedUniform
+    if Spaces.quadrature_style(space) isa Quadratures.ClosedUniform
         return space
     end
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
-    lagrange_quad = Spaces.Quadratures.ClosedUniform{Nq}()
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
+    lagrange_quad = Quadratures.ClosedUniform{Nq}()
     return Spaces.SpectralElementSpace1D(Spaces.topology(space), lagrange_quad)
 end
 function vtk_grid_space(space::Spaces.SpectralElementSpace2D)
-    if Spaces.quadrature_style(space) isa Spaces.Quadratures.ClosedUniform
+    if Spaces.quadrature_style(space) isa Quadratures.ClosedUniform
         return space
     end
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
-    lagrange_quad = Spaces.Quadratures.ClosedUniform{Nq}()
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
+    lagrange_quad = Quadratures.ClosedUniform{Nq}()
     return Spaces.SpectralElementSpace2D(Spaces.topology(space), lagrange_quad)
 end
 function vtk_grid_space(space::Spaces.FaceExtrudedFiniteDifferenceSpace)
@@ -227,20 +227,20 @@ This generally does two things:
 """
 function vtk_cell_space(gridspace::Spaces.SpectralElementSpace1D)
     @assert Spaces.quadrature_style(gridspace) isa
-            Spaces.Quadratures.ClosedUniform
-    Nq = Spaces.Quadratures.degrees_of_freedom(
+            Quadratures.ClosedUniform
+    Nq = Quadratures.degrees_of_freedom(
         Spaces.quadrature_style(gridspace),
     )
-    quad = Spaces.Quadratures.Uniform{Nq - 1}()
+    quad = Quadratures.Uniform{Nq - 1}()
     return Spaces.SpectralElementSpace1D(Spaces.topology(gridspace), quad)
 end
 function vtk_cell_space(gridspace::Spaces.SpectralElementSpace2D)
     @assert Spaces.quadrature_style(gridspace) isa
-            Spaces.Quadratures.ClosedUniform
-    Nq = Spaces.Quadratures.degrees_of_freedom(
+            Quadratures.ClosedUniform
+    Nq = Quadratures.degrees_of_freedom(
         Spaces.quadrature_style(gridspace),
     )
-    quad = Spaces.Quadratures.Uniform{Nq - 1}()
+    quad = Quadratures.Uniform{Nq - 1}()
     return Spaces.SpectralElementSpace2D(Spaces.topology(gridspace), quad)
 end
 function vtk_cell_space(gridspace::Spaces.FaceExtrudedFiniteDifferenceSpace)
diff --git a/lib/ClimaCoreVTK/test/runtests.jl b/lib/ClimaCoreVTK/test/runtests.jl
index f4dc39c2f5..cc05ab0c67 100644
--- a/lib/ClimaCoreVTK/test/runtests.jl
+++ b/lib/ClimaCoreVTK/test/runtests.jl
@@ -3,7 +3,7 @@ using ClimaComms
 using ClimaCoreVTK
 using IntervalSets
 import ClimaCore:
-    Geometry, Domains, Meshes, Topologies, Spaces, Fields, Operators
+    Geometry, Domains, Meshes, Topologies, Spaces, Fields, Operators, Quadratures
 
 OUTPUT_DIR = mkpath(get(ENV, "CI_OUTPUT_DIR", tempname()))
 mkpath(joinpath(OUTPUT_DIR, "series"))
@@ -15,7 +15,7 @@ mkpath(joinpath(OUTPUT_DIR, "series"))
     mesh = Meshes.EquiangularCubedSphere(domain, 4)
     grid_topology =
         Topologies.Topology2D(ClimaComms.SingletonCommsContext(), mesh)
-    quad = Spaces.Quadratures.GLL{5}()
+    quad = Quadratures.GLL{5}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
     coords = Fields.coordinate_field(space)
 
@@ -140,7 +140,7 @@ end
     mesh = Meshes.RectilinearMesh(domain, n1, n2)
     grid_topology =
         Topologies.Topology2D(ClimaComms.SingletonCommsContext(), mesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
     coords = Fields.coordinate_field(space)
@@ -164,7 +164,7 @@ end
     )
     hmesh = Meshes.IntervalMesh(hdomain, nelems = 10)
     htopology = Topologies.IntervalTopology(hmesh)
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     hspace = Spaces.SpectralElementSpace1D(htopology, quad)
 
     vdomain = Domains.IntervalDomain(
@@ -202,7 +202,7 @@ end
 
     hmesh = Meshes.RectilinearMesh(hdomain, 4, 4)
     htopology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), hmesh)
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     hspace = Spaces.SpectralElementSpace2D(htopology, quad)
 
     vdomain = Domains.IntervalDomain(
@@ -235,7 +235,7 @@ end
     hdomain = Domains.SphereDomain(R)
     hmesh = Meshes.EquiangularCubedSphere(hdomain, 4)
     htopology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), hmesh)
-    quad = Spaces.Quadratures.GLL{5}()
+    quad = Quadratures.GLL{5}()
     hspace = Spaces.SpectralElementSpace2D(htopology, quad)
 
 
diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl
index 0af45218b8..2ec4bc808c 100644
--- a/src/Fields/Fields.jl
+++ b/src/Fields/Fields.jl
@@ -5,6 +5,7 @@ import ..slab, ..slab_args, ..column, ..column_args, ..level
 import ..DataLayouts: DataLayouts, AbstractData, DataStyle
 import ..Domains
 import ..Topologies
+import ..Quadratures
 import ..Grids: ColumnIndex
 import ..Spaces: Spaces, AbstractSpace, AbstractPointSpace
 import ..Geometry: Geometry, Cartesian12Vector
diff --git a/src/Fields/indices.jl b/src/Fields/indices.jl
index 89cac213b7..72d4bc11eb 100644
--- a/src/Fields/indices.jl
+++ b/src/Fields/indices.jl
@@ -67,7 +67,7 @@ function bycolumn(
     ::ClimaComms.CPUSingleThreaded,
 )
     Nh = Topologies.nlocalelems(space)
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
     @inbounds begin
         for h in 1:Nh
             for i in 1:Nq
@@ -83,7 +83,7 @@ function bycolumn(
     ::ClimaComms.CPUMultiThreaded,
 )
     Nh = Topologies.nlocalelems(space)
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
     @inbounds begin
         Threads.@threads for h in 1:Nh
             for i in 1:Nq
@@ -99,7 +99,7 @@ function bycolumn(
     ::ClimaComms.CPUSingleThreaded,
 )
     Nh = Topologies.nlocalelems(space)
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
     @inbounds begin
         for h in 1:Nh
             for j in 1:Nq, i in 1:Nq
@@ -115,7 +115,7 @@ function bycolumn(
     ::ClimaComms.CPUMultiThreaded,
 )
     Nh = Topologies.nlocalelems(space)
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
     @inbounds begin
         Threads.@threads for h in 1:Nh
             for j in 1:Nq, i in 1:Nq
@@ -152,12 +152,12 @@ ncolumns(space::Spaces.ExtrudedFiniteDifferenceSpace) =
 
 function ncolumns(space::Spaces.SpectralElementSpace1D)
     Nh = Topologies.nlocalelems(space)
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
     return Nh * Nq
 end
 function ncolumns(space::Spaces.SpectralElementSpace2D)
     Nh = Topologies.nlocalelems(space)
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
     return Nh * Nq * Nq
 end
 
diff --git a/src/Operators/numericalflux.jl b/src/Operators/numericalflux.jl
index d767c4cedc..7d965d8103 100644
--- a/src/Operators/numericalflux.jl
+++ b/src/Operators/numericalflux.jl
@@ -24,7 +24,7 @@ See also:
 """
 function add_numerical_flux_internal!(fn, dydt, args...)
     space = axes(dydt)
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
     topology = Spaces.topology(space)
 
     for (iface, (elem⁻, face⁻, elem⁺, face⁺, reversed)) in
@@ -100,7 +100,7 @@ end
 
 function add_numerical_flux_boundary!(fn, dydt, args...)
     space = axes(dydt)
-    Nq = Spaces.Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
+    Nq = Quadratures.degrees_of_freedom(Spaces.quadrature_style(space))
     topology = Spaces.topology(space)
 
     for (iboundary, boundarytag) in
diff --git a/src/Remapping/Remapping.jl b/src/Remapping/Remapping.jl
index 45807ab19b..3432c97bb3 100644
--- a/src/Remapping/Remapping.jl
+++ b/src/Remapping/Remapping.jl
@@ -10,6 +10,7 @@ import ..DataLayouts,
     ..Topologies,
     ..Meshes,
     ..Operators,
+    ..Quadratures,
     ..Fields,
     ..Hypsography
 
diff --git a/src/Remapping/distributed_remapping.jl b/src/Remapping/distributed_remapping.jl
index 7617730fd7..76a5cb867b 100644
--- a/src/Remapping/distributed_remapping.jl
+++ b/src/Remapping/distributed_remapping.jl
@@ -151,7 +151,7 @@ function Remapper(target_hcoords, target_zcoords, space)
 
     interpolation_coeffs = map(local_target_hcoords) do hcoord
         quad = Spaces.quadrature_style(space)
-        quad_points, _ = Spaces.Quadratures.quadrature_points(FT, quad)
+        quad_points, _ = Quadratures.quadrature_points(FT, quad)
         return interpolation_weights(horz_mesh, hcoord, quad_points)
     end
 
diff --git a/src/Remapping/interpolate_array.jl b/src/Remapping/interpolate_array.jl
index f1a8531f80..725de231ab 100644
--- a/src/Remapping/interpolate_array.jl
+++ b/src/Remapping/interpolate_array.jl
@@ -458,7 +458,7 @@ function interpolate_array(
         hcoord = xcoord
         helem = Meshes.containing_element(horz_mesh, hcoord)
         quad = Spaces.quadrature_style(space)
-        quad_points, _ = Spaces.Quadratures.quadrature_points(FT, quad)
+        quad_points, _ = Quadratures.quadrature_points(FT, quad)
         weights = interpolation_weights(horz_mesh, hcoord, quad_points)
         h = helem
 
@@ -497,7 +497,7 @@ function interpolate_array(
         hcoord = Geometry.product_coordinates(xcoord, ycoord)
         helem = Meshes.containing_element(horz_mesh, hcoord)
         quad = Spaces.quadrature_style(space)
-        quad_points, _ = Spaces.Quadratures.quadrature_points(FT, quad)
+        quad_points, _ = Quadratures.quadrature_points(FT, quad)
         weights = interpolation_weights(horz_mesh, hcoord, quad_points)
         gidx = horz_topology.orderindex[helem]
         h = gidx
@@ -528,8 +528,8 @@ function interpolation_weights(
 )
     helem = Meshes.containing_element(horz_mesh, hcoord)
     ξ1, ξ2 = Meshes.reference_coordinates(horz_mesh, helem, hcoord)
-    WI1 = Spaces.Quadratures.interpolation_matrix(SVector(ξ1), quad_points)
-    WI2 = Spaces.Quadratures.interpolation_matrix(SVector(ξ2), quad_points)
+    WI1 = Quadratures.interpolation_matrix(SVector(ξ1), quad_points)
+    WI2 = Quadratures.interpolation_matrix(SVector(ξ2), quad_points)
     return (WI1, WI2)
 end
 
@@ -540,7 +540,7 @@ function interpolation_weights(
 )
     helem = Meshes.containing_element(horz_mesh, hcoord)
     ξ1, = Meshes.reference_coordinates(horz_mesh, helem, hcoord)
-    WI1 = Spaces.Quadratures.interpolation_matrix(SVector(ξ1), quad_points)
+    WI1 = Quadratures.interpolation_matrix(SVector(ξ1), quad_points)
     return (WI1,)
 end
 
@@ -550,7 +550,7 @@ end
 Interpolate the given `field` on the given points assuming the given interpolation_matrix
 and global index in the topology.
 
-The coefficients `weights` are computed with `Spaces.Quadratures.interpolation_matrix`.
+The coefficients `weights` are computed with `Quadratures.interpolation_matrix`.
 See also `interpolate_array`.
 
 Keyword arguments
diff --git a/src/Spaces/quadrature.jl b/src/Spaces/quadrature.jl
deleted file mode 100644
index 5e9b892904..0000000000
--- a/src/Spaces/quadrature.jl
+++ /dev/null
@@ -1,286 +0,0 @@
-
-module Quadratures
-
-import GaussQuadrature
-import StaticArrays: SVector, SMatrix, MMatrix
-import LinearAlgebra: Diagonal
-
-export QuadratureStyle,
-    GLL, GL, polynomial_degree, degrees_of_freedom, quadrature_points
-
-"""
-   QuadratureStyle
-
-Quadrature style supertype. See sub-types:
- - [`GLL`](@ref)
- - [`GL`](@ref)
- - [`Uniform`](@ref)
-"""
-abstract type QuadratureStyle end
-
-"""
-    polynomial_degree(QuadratureStyle) -> Int
-
-Returns the polynomial degree of the `QuadratureStyle` concrete type
-"""
-function polynomial_degree end
-
-
-"""
-    degrees_of_freedom(QuadratureStyle) -> Int
-
-Returns the degrees_of_freedom of the `QuadratureStyle` concrete type
-"""
-function degrees_of_freedom end
-
-"""
-    points, weights = quadrature_points(::Type{FT}, quadrature_style)
-
-The points and weights of the quadrature rule in floating point type `FT`.
-"""
-function quadrature_points end
-
-
-"""
-    GLL{Nq}()
-
-Gauss-Legendre-Lobatto quadrature using `Nq` quadrature points.
-"""
-struct GLL{Nq} <: QuadratureStyle end
-
-Base.show(io::IO, ::GLL{Nq}) where {Nq} =
-    print(io, Nq, "-point Gauss-Legendre-Lobatto quadrature")
-
-@inline polynomial_degree(::GLL{Nq}) where {Nq} = Int(Nq - 1)
-@inline degrees_of_freedom(::GLL{Nq}) where {Nq} = Int(Nq)
-
-@generated function quadrature_points(::Type{FT}, ::GLL{Nq}) where {FT, Nq}
-    points, weights = GaussQuadrature.legendre(FT, Nq, GaussQuadrature.both)
-    :($(SVector{Nq}(points)), $(SVector{Nq}(weights)))
-end
-
-"""
-    GL{Nq}()
-
-Gauss-Legendre quadrature using `Nq` quadrature points.
-"""
-struct GL{Nq} <: QuadratureStyle end
-Base.show(io::IO, ::GL{Nq}) where {Nq} =
-    print(io, Nq, "-point Gauss-Legendre quadrature")
-
-@inline polynomial_degree(::GL{Nq}) where {Nq} = Int(Nq - 1)
-@inline degrees_of_freedom(::GL{Nq}) where {Nq} = Int(Nq)
-
-@generated function quadrature_points(::Type{FT}, ::GL{Nq}) where {FT, Nq}
-    points, weights = GaussQuadrature.legendre(FT, Nq, GaussQuadrature.neither)
-    :($(SVector{Nq}(points)), $(SVector{Nq}(weights)))
-end
-
-"""
-    Uniform{Nq}()
-
-Uniformly-spaced quadrature.
-"""
-struct Uniform{Nq} <: QuadratureStyle end
-
-@inline polynomial_degree(::Uniform{Nq}) where {Nq} = Int(Nq - 1)
-@inline degrees_of_freedom(::Uniform{Nq}) where {Nq} = Int(Nq)
-
-@generated function quadrature_points(::Type{FT}, ::Uniform{Nq}) where {FT, Nq}
-    points = SVector{Nq}(range(-1 + 1 / Nq, step = 2 / Nq, length = Nq))
-    weights = SVector{Nq}(ntuple(i -> 2 / Nq, Nq))
-    :($points, $weights)
-end
-
-"""
-    ClosedUniform{Nq}()
-
-Uniformly-spaced quadrature including boundary.
-"""
-struct ClosedUniform{Nq} <: QuadratureStyle end
-
-@inline polynomial_degree(::ClosedUniform{Nq}) where {Nq} = Int(Nq - 1)
-@inline degrees_of_freedom(::ClosedUniform{Nq}) where {Nq} = Int(Nq)
-
-@generated function quadrature_points(
-    ::Type{FT},
-    ::ClosedUniform{Nq},
-) where {FT, Nq}
-    points = SVector{Nq}(range(FT(-1), FT(1), length = Nq))
-    weights = SVector{Nq}(
-        1 / (Nq - 1),
-        ntuple(i -> 2 / (Nq - 1), Nq - 2)...,
-        1 / (Nq - 1),
-    )
-    :($points, $weights)
-end
-
-
-"""
-    barycentric_weights(x::SVector{Nq}) where {Nq}
-
-The barycentric weights associated with the array of point locations `x`:
-
-```math
-w_j = \\frac{1}{\\prod_{k \\ne j} (x_i - x_j)}
-```
-
-See [Berrut2004](@cite), equation 3.2.
-"""
-function barycentric_weights(r::SVector{Nq, T}) where {Nq, T}
-    SVector{Nq}(ntuple(Nq) do i
-        w = one(T)
-        for j in 1:Nq
-            if j != i
-                w *= (r[j] - r[i])
-            end
-        end
-        inv(w)
-    end)
-end
-@generated function barycentric_weights(
-    ::Type{FT},
-    quadstyle::QuadratureStyle,
-) where {FT}
-    barycentric_weights(quadrature_points(FT, quadstyle())[1])
-end
-
-"""
-    interpolation_matrix(x::SVector, r::SVector{Nq})
-
-The matrix which interpolates the Lagrange polynomial of degree `Nq-1` through
-the points `r`, to points `x`. The matrix coefficients are computed using the
-Barycentric formula of [Berrut2004](@cite), section 4:
-```math
-I_{ij} = \\begin{cases}
-1 & \\text{if } x_i = r_j, \\\\
-0 & \\text{if } x_i = r_k \\text{ for } k \\ne j, \\\\
-\\frac{\\displaystyle \\frac{w_j}{x_i - r_j}}{\\displaystyle \\sum_k \\frac{w_k}{x_i - r_k}} & \\text{otherwise,}
-\\end{cases}
-```
-where ``w_j`` are the barycentric weights, see [`barycentric_weights`](@ref).
-"""
-function interpolation_matrix(
-    points_to::SVector{Nto},
-    points_from::SVector{Nfrom},
-) where {Nto, Nfrom}
-    T = eltype(points_to)
-    bw = barycentric_weights(points_from)
-    M = zeros(MMatrix{Nto, Nfrom, T, Nto * Nfrom})
-    for i in 1:Nto
-        x_to = points_to[i]
-        skip_row = false
-        for j in 1:Nfrom
-            if x_to == points_from[j]
-                # assign to one to avoid singularity condition
-                M[i, j] = one(T)
-                # skip over the equal boundry condition
-                skip_row = true
-            end
-            skip_row && break
-        end
-        skip_row && continue
-        w = bw ./ (x_to .- points_from)
-        M[i, :] .= w ./ sum(w)
-    end
-    return SMatrix(M)
-end
-
-@generated function interpolation_matrix(
-    ::Type{FT},
-    quadto::QuadratureStyle,
-    quadfrom::QuadratureStyle,
-) where {FT}
-    interpolation_matrix(
-        quadrature_points(FT, quadto())[1],
-        quadrature_points(FT, quadfrom())[1],
-    )
-end
-
-"""
-    V = orthonormal_poly(points, quad)
-
-`V_{ij}` contains the `j-1`th Legendre polynomial evaluated at `points[i]`.
-i.e. it is the mapping from the modal to the nodal representation.
-"""
-function orthonormal_poly(
-    points::SVector{Np, FT},
-    quad::GLL{Nq},
-) where {FT, Np, Nq}
-    N = Nq - 1
-    a, b = GaussQuadrature.legendre_coefs(FT, N)
-    if N == 0
-        return SMatrix{Np, 1}(ntuple(x -> b[1], Np))
-    end
-    return SMatrix{Np, Nq}(GaussQuadrature.orthonormal_poly(points, a, b))
-end
-
-function spectral_filter_matrix(
-    quad::GLL{Nq},
-    Σ::SVector{Nq, FT},
-) where {Nq, FT}
-    points, _ = quadrature_points(FT, quad)
-    V = orthonormal_poly(points, quad)
-    return V * Diagonal(Σ) / V
-end
-
-function cutoff_filter_matrix(
-    ::Type{FT},
-    quad::GLL{Nq},
-    Nc::Integer,
-) where {FT, Nq}
-    Σ = SVector(ntuple(i -> i <= Nc ? FT(1) : FT(0), Nq))
-    return spectral_filter_matrix(quad, Σ)
-end
-
-"""
-    differentiation_matrix(r::SVector{Nq, T}) where {Nq, T}
-
-The spectral differentiation matrix for the Lagrange polynomial of degree `Nq-1`
-interpolating at points `r`.
-
-The matrix coefficients are computed using the [Berrut2004](@cite), section 9.3:
-```math
-D_{ij} = \\begin{cases}
-    \\displaystyle
-    \\frac{w_j}{w_i (x_i - x_j)} &\\text{ if } i \\ne j \\\\
-    -\\sum_{k \\ne j} D_{kj} &\\text{ if } i = j
-\\end{cases}
-```
-where ``w_j`` are the barycentric weights, see [`barycentric_weights`](@ref).
-"""
-function differentiation_matrix(r::SVector{Nq, T}) where {Nq, T}
-    wb = barycentric_weights(r)
-    SMatrix{Nq, Nq, T, Nq * Nq}(
-        begin
-            if i == j
-                D = zero(T)
-                for l in 1:Nq
-                    if l != i
-                        D += one(T) / (r[i] - r[l])
-                    end
-                end
-                D
-            else
-                (wb[i] / wb[j]) / (r[j] - r[i])
-            end
-        end for j in 1:Nq, i in 1:Nq
-    )
-end
-
-"""
-    differentiation_matrix(FT, quadstyle::QuadratureStyle)
-
-The spectral differentiation matrix at the quadrature points of `quadstyle`,
-using floating point types `FT`.
-"""
-@generated function differentiation_matrix(
-    ::Type{FT},
-    quadstyle::QuadratureStyle,
-) where {FT}
-    differentiation_matrix(quadrature_points(FT, quadstyle())[1])
-end
-
-
-
-end # module
diff --git a/test/Fields/field.jl b/test/Fields/field.jl
index c1218873aa..3e10244cd7 100644
--- a/test/Fields/field.jl
+++ b/test/Fields/field.jl
@@ -35,7 +35,7 @@ function spectral_space_2D(; n1 = 1, n2 = 1, Nij = 4)
     grid_topology =
         Topologies.Topology2D(ClimaComms.SingletonCommsContext(device), mesh)
 
-    quad = Spaces.Quadratures.GLL{Nij}()
+    quad = Quadratures.GLL{Nij}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
     return space
 end
@@ -341,7 +341,7 @@ end
     mesh_xy = Meshes.RectilinearMesh(domain_xy, 10, 10)
     topology_xy = Topologies.Topology2D(context, mesh_xy)
 
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
 
     space_vf = Spaces.CenterFiniteDifferenceSpace(topology_z)
     space_ifh = Spaces.SpectralElementSpace1D(topology_x, quad)
@@ -680,7 +680,7 @@ end
     mesh_xy = Meshes.RectilinearMesh(domain_xy, 10, 10)
     topology_xy = Topologies.Topology2D(context, mesh_xy)
 
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
 
     space_vf = Spaces.CenterFiniteDifferenceSpace(topology_z)
     space_ifh = Spaces.SpectralElementSpace1D(topology_x, quad)
@@ -828,7 +828,7 @@ end
         TU.bycolumnable(space) || continue
         hspace = Spaces.horizontal_space(space)
         Nh = Topologies.nlocalelems(hspace)
-        Nq = Spaces.Quadratures.degrees_of_freedom(
+        Nq = Quadratures.degrees_of_freedom(
             Spaces.quadrature_style(hspace),
         )
         if nameof(typeof(space)) == :SpectralElementSpace1D
diff --git a/test/Fields/field_opt.jl b/test/Fields/field_opt.jl
index c0a9b1c5d6..295d7dffb6 100644
--- a/test/Fields/field_opt.jl
+++ b/test/Fields/field_opt.jl
@@ -5,7 +5,7 @@ import ClimaCore
 import ClimaCore.Utilities: PlusHalf, half
 import ClimaCore.DataLayouts: IJFH
 import ClimaCore:
-    Fields, slab, Domains, Topologies, Meshes, Operators, Spaces, Geometry
+    Fields, slab, Domains, Topologies, Meshes, Operators, Spaces, Geometry, Quadratures
 
 using FastBroadcast
 using LinearAlgebra: norm
@@ -265,7 +265,7 @@ end
         hdomain = Domains.SphereDomain(FT(1e7))
         hmesh = Meshes.EquiangularCubedSphere(hdomain, helem)
         htopology = Topologies.Topology2D(hmesh)
-        quad = Spaces.Quadratures.GLL{npoly + 1}()
+        quad = Quadratures.GLL{npoly + 1}()
         hspace = Spaces.SpectralElementSpace2D(htopology, quad)
         vdomain = Domains.IntervalDomain(
             Geometry.ZPoint{FT}(zero(FT)),
diff --git a/test/Fields/reduction_cuda.jl b/test/Fields/reduction_cuda.jl
index cc26a79911..9c2333d332 100644
--- a/test/Fields/reduction_cuda.jl
+++ b/test/Fields/reduction_cuda.jl
@@ -11,6 +11,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     DataLayouts
 
@@ -32,7 +33,7 @@ include("reduction_cuda_utils.jl")
     R = FT(6.37122e6) # radius of earth
     domain = Domains.SphereDomain(R)
     mesh = Meshes.EquiangularCubedSphere(domain, ne)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     grid_topology = Topologies.Topology2D(context, mesh)
     grid_topology_cpu = Topologies.Topology2D(context_cpu, mesh)
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
@@ -112,7 +113,7 @@ end
         horizontal_mesh,
         Topologies.spacefillingcurve(horizontal_mesh),
     )
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     h_space = Spaces.SpectralElementSpace2D(horizontal_topology, quad)
     h_space_cpu = Spaces.SpectralElementSpace2D(horizontal_topology_cpu, quad)
 
@@ -276,7 +277,7 @@ end
     mesh = Meshes.RectilinearMesh(domain, 3, 3)
     topology = Topologies.Topology2D(comms_ctx, mesh)
     topology_cpu = Topologies.Topology2D(comms_ctx_cpu, mesh)
-    quad = Spaces.Quadratures.GLL{5}()
+    quad = Quadratures.GLL{5}()
     space = Spaces.SpectralElementSpace2D(topology, quad)
     space_cpu = Spaces.SpectralElementSpace2D(topology_cpu, quad)
     coords = Fields.coordinate_field(space)
diff --git a/test/Fields/reduction_cuda_distributed.jl b/test/Fields/reduction_cuda_distributed.jl
index 9d5ae378c5..475c39cba2 100644
--- a/test/Fields/reduction_cuda_distributed.jl
+++ b/test/Fields/reduction_cuda_distributed.jl
@@ -12,6 +12,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     DataLayouts
 
@@ -35,7 +36,7 @@ include("reduction_cuda_utils.jl")
     R = FT(6.37122e6) # radius of earth
     domain = Domains.SphereDomain(R)
     mesh = Meshes.EquiangularCubedSphere(domain, ne)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     grid_topology = Topologies.Topology2D(context, mesh)
     grid_topology_cpu = Topologies.Topology2D(context_cpu, mesh)
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
@@ -122,7 +123,7 @@ end
         horizontal_mesh,
         Topologies.spacefillingcurve(horizontal_mesh),
     )
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     h_space = Spaces.SpectralElementSpace2D(horizontal_topology, quad)
     h_space_cpu = Spaces.SpectralElementSpace2D(horizontal_topology_cpu, quad)
 
diff --git a/test/Hypsography/2d.jl b/test/Hypsography/2d.jl
index e061dc24a4..30ee941228 100644
--- a/test/Hypsography/2d.jl
+++ b/test/Hypsography/2d.jl
@@ -3,12 +3,12 @@ import ClimaCore
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators,
     InputOutput,
@@ -33,7 +33,7 @@ horzdomain = Domains.IntervalDomain(
 )
 horzmesh = Meshes.IntervalMesh(horzdomain, nelems = 20)
 horztopology = Topologies.IntervalTopology(horzmesh)
-quad = Spaces.Quadratures.GLL{4 + 1}()
+quad = Quadratures.GLL{4 + 1}()
 horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)
 
 z_surface = sin.(Fields.coordinate_field(horzspace).x) .+ 1
diff --git a/test/Hypsography/3dsphere.jl b/test/Hypsography/3dsphere.jl
index d7d0edce90..00889f2083 100644
--- a/test/Hypsography/3dsphere.jl
+++ b/test/Hypsography/3dsphere.jl
@@ -4,12 +4,12 @@ import ClimaCore
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators,
     InputOutput,
@@ -33,7 +33,7 @@ horztopology = Topologies.Topology2D(
     ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()),
     horzmesh,
 )
-quad = Spaces.Quadratures.GLL{4 + 1}()
+quad = Quadratures.GLL{4 + 1}()
 horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
 z_surface =
diff --git a/test/InputOutput/hybrid2dbox.jl b/test/InputOutput/hybrid2dbox.jl
index 81f7e2fe4e..63263d36d5 100644
--- a/test/InputOutput/hybrid2dbox.jl
+++ b/test/InputOutput/hybrid2dbox.jl
@@ -5,12 +5,12 @@ import ClimaComms
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators,
     InputOutput
@@ -41,7 +41,7 @@ function hvspace_2D(
     horzmesh = Meshes.IntervalMesh(horzdomain, nelems = xelem)
     horztopology = Topologies.IntervalTopology(context, horzmesh)
 
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)
 
     hv_center_space =
diff --git a/test/InputOutput/hybrid2dbox_stretched.jl b/test/InputOutput/hybrid2dbox_stretched.jl
index a5793d1a23..670b8e26b6 100644
--- a/test/InputOutput/hybrid2dbox_stretched.jl
+++ b/test/InputOutput/hybrid2dbox_stretched.jl
@@ -6,12 +6,12 @@ import ClimaComms
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators,
     InputOutput
@@ -43,7 +43,7 @@ function hvspace_2D(
     horzmesh = Meshes.IntervalMesh(horzdomain, nelems = xelem)
     horztopology = Topologies.IntervalTopology(context, horzmesh)
 
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)
 
     hv_center_space =
diff --git a/test/InputOutput/hybrid2dbox_topography.jl b/test/InputOutput/hybrid2dbox_topography.jl
index 7225a8bd8b..e19d9834da 100644
--- a/test/InputOutput/hybrid2dbox_topography.jl
+++ b/test/InputOutput/hybrid2dbox_topography.jl
@@ -5,12 +5,12 @@ import ClimaComms
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators,
     Hypsography,
@@ -42,7 +42,7 @@ function hvspace_2D(
     horzmesh = Meshes.IntervalMesh(horzdomain, nelems = xelem)
     horztopology = Topologies.IntervalTopology(context, horzmesh)
 
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)
 
     z_surface =
diff --git a/test/InputOutput/hybrid3dbox.jl b/test/InputOutput/hybrid3dbox.jl
index aa4ba23324..ffd2a319c0 100644
--- a/test/InputOutput/hybrid3dbox.jl
+++ b/test/InputOutput/hybrid3dbox.jl
@@ -5,13 +5,13 @@ import ClimaCore
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     DataLayouts,
     Spaces,
+    Quadratures,
     Fields,
     Operators,
     InputOutput
@@ -43,7 +43,7 @@ function hvspace_3D(
     )
     Nv = Meshes.nelements(vertmesh)
     Nf_center, Nf_face = 2, 1 #1 + 3 + 1
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzmesh = Meshes.RectilinearMesh(horzdomain, xelem, yelem)
     context = ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded())
     horztopology = Topologies.Topology2D(context, horzmesh)
diff --git a/test/InputOutput/hybrid3dcubedsphere.jl b/test/InputOutput/hybrid3dcubedsphere.jl
index 1a3b4edf37..45cd0ce185 100644
--- a/test/InputOutput/hybrid3dcubedsphere.jl
+++ b/test/InputOutput/hybrid3dcubedsphere.jl
@@ -6,6 +6,7 @@ using ClimaCore:
     Meshes,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     DataLayouts,
     InputOutput
@@ -34,7 +35,7 @@ end
         horizontal_mesh,
         Topologies.spacefillingcurve(horizontal_mesh),
     )
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     h_space = Spaces.SpectralElementSpace2D(topology, quad)
     # vertical space
     z_domain = Domains.IntervalDomain(
diff --git a/test/InputOutput/hybrid3dcubedsphere_topography.jl b/test/InputOutput/hybrid3dcubedsphere_topography.jl
index a03eba871a..981c0fab29 100644
--- a/test/InputOutput/hybrid3dcubedsphere_topography.jl
+++ b/test/InputOutput/hybrid3dcubedsphere_topography.jl
@@ -6,6 +6,7 @@ using ClimaCore:
     Meshes,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     DataLayouts,
     Hypsography,
@@ -35,7 +36,7 @@ end
         horizontal_mesh,
         Topologies.spacefillingcurve(horizontal_mesh),
     )
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     h_space = Spaces.SpectralElementSpace2D(topology, quad)
     # vertical space
     z_domain = Domains.IntervalDomain(
diff --git a/test/InputOutput/spectralelement2d.jl b/test/InputOutput/spectralelement2d.jl
index 961245aeb1..4ab617ba21 100644
--- a/test/InputOutput/spectralelement2d.jl
+++ b/test/InputOutput/spectralelement2d.jl
@@ -9,6 +9,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     DataLayouts,
     InputOutput
@@ -66,7 +67,7 @@ end
     )
     n1, n2 = 16, 16
     Nq = 4
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     mesh = Meshes.RectilinearMesh(domain, n1, n2)
     device = ClimaComms.device()
     @info "Using device" device
diff --git a/test/Limiters/distributed/dlimiter.jl b/test/Limiters/distributed/dlimiter.jl
index b8df70c1f9..17d93bb0d3 100644
--- a/test/Limiters/distributed/dlimiter.jl
+++ b/test/Limiters/distributed/dlimiter.jl
@@ -1,5 +1,5 @@
 using ClimaCore:
-    DataLayouts, Fields, Domains, Geometry, Topologies, Meshes, Spaces, Limiters
+    DataLayouts, Fields, Domains, Geometry, Topologies, Meshes, Spaces, Limiters, Quadratures
 using ClimaCore.RecursiveApply
 using ClimaCore: slab
 using Test
@@ -41,7 +41,7 @@ zdomain = Domains.IntervalDomain(
 vertmesh = Meshes.IntervalMesh(zdomain, nelems = zelems)
 vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)
 
-quad = Spaces.Quadratures.GLL{Nij}()
+quad = Quadratures.GLL{Nij}()
 horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
 hv_center_space =
diff --git a/test/Limiters/limiter.jl b/test/Limiters/limiter.jl
index 7895791cae..02eeca984e 100644
--- a/test/Limiters/limiter.jl
+++ b/test/Limiters/limiter.jl
@@ -6,7 +6,7 @@ import CUDA
 CUDA.allowscalar(false)
 using ClimaComms
 using ClimaCore:
-    DataLayouts, Fields, Domains, Geometry, Topologies, Meshes, Spaces, Limiters
+    DataLayouts, Fields, Domains, Geometry, Topologies, Meshes, Spaces, Limiters, Quadratures
 using ClimaCore.RecursiveApply
 using ClimaCore: slab
 using Test
@@ -44,7 +44,7 @@ function rectangular_mesh_space(
     )
     mesh = Meshes.RectilinearMesh(domain, n1, n2)
     topology = Topologies.Topology2D(comms_ctx, mesh)
-    quad = Spaces.Quadratures.GLL{Nij}()
+    quad = Quadratures.GLL{Nij}()
     return Spaces.SpectralElementSpace2D(topology, quad)
 end
 
@@ -87,7 +87,7 @@ function hvspace_3D(
 
     vert_center_space = Spaces.CenterFiniteDifferenceSpace(z_topology)
 
-    quad = Spaces.Quadratures.GLL{Nij}()
+    quad = Quadratures.GLL{Nij}()
     horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
     hv_center_space =
diff --git a/test/MatrixFields/matrix_field_test_utils.jl b/test/MatrixFields/matrix_field_test_utils.jl
index e71ed5b999..cc3790da9f 100644
--- a/test/MatrixFields/matrix_field_test_utils.jl
+++ b/test/MatrixFields/matrix_field_test_utils.jl
@@ -5,7 +5,7 @@ import Random: seed!
 
 import ClimaComms
 import ClimaCore:
-    Geometry, Domains, Meshes, Topologies, Hypsography, Spaces, Fields
+    Geometry, Domains, Meshes, Topologies, Hypsography, Spaces, Fields, Quadratures
 using ClimaCore.MatrixFields
 
 # Test that an expression is true and that it is also type-stable.
@@ -195,7 +195,7 @@ function test_spaces(::Type{FT}) where {FT}
     hdomain = Domains.SphereDomain(FT(10))
     hmesh = Meshes.EquiangularCubedSphere(hdomain, helem)
     htopology = Topologies.Topology2D(comms_ctx, hmesh)
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     hspace = Spaces.SpectralElementSpace2D(htopology, quad)
     vdomain = Domains.IntervalDomain(
         Geometry.ZPoint(FT(0)),
diff --git a/test/Operators/finitedifference/implicit_stencils_utils.jl b/test/Operators/finitedifference/implicit_stencils_utils.jl
index fb65fcee14..4142e712ed 100644
--- a/test/Operators/finitedifference/implicit_stencils_utils.jl
+++ b/test/Operators/finitedifference/implicit_stencils_utils.jl
@@ -3,7 +3,7 @@ using ClimaComms
 using Random: seed!
 seed!(1) # ensures reproducibility
 
-using ClimaCore: Geometry, Domains, Meshes, Topologies, Spaces, Fields
+using ClimaCore: Geometry, Domains, Meshes, Topologies, Spaces, Fields, Quadratures
 using ClimaCore: Operators
 import ClimaCore.Operators as OP
 
@@ -62,7 +62,7 @@ function get_space(::Type{FT}) where {FT}
     hdomain = Domains.SphereDomain(radius)
     hmesh = Meshes.EquiangularCubedSphere(hdomain, helem)
     htopology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), hmesh)
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     hspace = Spaces.SpectralElementSpace2D(htopology, quad)
 
     vdomain = Domains.IntervalDomain(
diff --git a/test/Operators/finitedifference/linsolve.jl b/test/Operators/finitedifference/linsolve.jl
index 5ab3e6549d..f87deaa149 100644
--- a/test/Operators/finitedifference/linsolve.jl
+++ b/test/Operators/finitedifference/linsolve.jl
@@ -5,7 +5,7 @@ import ClimaCore
 # To avoid JET failures in the error message
 ClimaCore.Operators.allow_mismatched_fd_spaces() = true
 
-using ClimaCore: Geometry, Domains, Meshes, Topologies, Spaces, Fields
+using ClimaCore: Geometry, Domains, Meshes, Topologies, Spaces, Fields, Quadratures
 
 FT = Float32
 radius = FT(1e7)
@@ -16,7 +16,7 @@ velem = 4
 hdomain = Domains.SphereDomain(radius)
 hmesh = Meshes.EquiangularCubedSphere(hdomain, helem)
 htopology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), hmesh)
-quad = Spaces.Quadratures.GLL{npoly + 1}()
+quad = Quadratures.GLL{npoly + 1}()
 hspace = Spaces.SpectralElementSpace2D(htopology, quad)
 
 vdomain = Domains.IntervalDomain(
diff --git a/test/Operators/finitedifference/opt_implicit_stencils.jl b/test/Operators/finitedifference/opt_implicit_stencils.jl
index 40b2cf0522..e74708500a 100644
--- a/test/Operators/finitedifference/opt_implicit_stencils.jl
+++ b/test/Operators/finitedifference/opt_implicit_stencils.jl
@@ -2,7 +2,7 @@ using Test
 using ClimaComms
 using Random: seed!
 
-using ClimaCore: Geometry, Domains, Meshes, Topologies, Spaces, Fields
+using ClimaCore: Geometry, Domains, Meshes, Topologies, Spaces, Fields, Quadratures
 using ClimaCore: Operators
 
 struct CurriedTwoArgOperator{O, A}
@@ -29,7 +29,7 @@ Operators.Operator2Stencil(op::CurriedTwoArgOperator) =
     hdomain = Domains.SphereDomain(radius)
     hmesh = Meshes.EquiangularCubedSphere(hdomain, helem)
     htopology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), hmesh)
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     hspace = Spaces.SpectralElementSpace2D(htopology, quad)
 
     vdomain = Domains.IntervalDomain(
diff --git a/test/Operators/finitedifference/tensor.jl b/test/Operators/finitedifference/tensor.jl
index 56345c00a1..38e9ab1d45 100644
--- a/test/Operators/finitedifference/tensor.jl
+++ b/test/Operators/finitedifference/tensor.jl
@@ -2,7 +2,7 @@ using Test
 
 using ClimaComms
 using ClimaCore:
-    Geometry, Domains, Meshes, Topologies, Spaces, Fields, Operators
+    Geometry, Domains, Meshes, Topologies, Spaces, Fields, Operators, Quadratures
 using LinearAlgebra
 
 for FT in (Float32, Float64)
@@ -10,7 +10,7 @@ for FT in (Float32, Float64)
     hmesh = Meshes.EquiangularCubedSphere(hdomain, 30)
     htopology = Topologies.Topology2D(hmesh)
     hspace =
-        Spaces.SpectralElementSpace2D(htopology, Spaces.Quadratures.GLL{4}())
+        Spaces.SpectralElementSpace2D(htopology, Quadratures.GLL{4}())
 
     vdomain = Domains.IntervalDomain(
         Geometry.ZPoint{FT}(0.0),
diff --git a/test/Operators/finitedifference/wfact.jl b/test/Operators/finitedifference/wfact.jl
index 592c5f8e14..9f8afca9fb 100644
--- a/test/Operators/finitedifference/wfact.jl
+++ b/test/Operators/finitedifference/wfact.jl
@@ -6,7 +6,7 @@ import ClimaCore
 ClimaCore.Operators.allow_mismatched_fd_spaces() = true
 
 using ClimaCore:
-    Geometry, Domains, Meshes, Topologies, Spaces, Fields, Operators
+    Geometry, Domains, Meshes, Topologies, Spaces, Fields, Operators, Quadratures
 
 import ClimaCore.Utilities: half
 import LinearAlgebra: norm_sqr
@@ -20,7 +20,7 @@ velem = 4
 hdomain = Domains.SphereDomain(radius)
 hmesh = Meshes.EquiangularCubedSphere(hdomain, helem)
 htopology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), hmesh)
-quad = Spaces.Quadratures.GLL{npoly + 1}()
+quad = Quadratures.GLL{npoly + 1}()
 hspace = Spaces.SpectralElementSpace2D(htopology, quad)
 
 vdomain = Domains.IntervalDomain(
diff --git a/test/Operators/hybrid/2d.jl b/test/Operators/hybrid/2d.jl
index 6931fb1a5d..21e3a401de 100644
--- a/test/Operators/hybrid/2d.jl
+++ b/test/Operators/hybrid/2d.jl
@@ -4,12 +4,12 @@ using StaticArrays, IntervalSets, LinearAlgebra
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators,
     level
@@ -44,7 +44,7 @@ function hvspace_2D(;
     horzmesh = Meshes.IntervalMesh(horzdomain, nelems = helem)
     horztopology = Topologies.IntervalTopology(horzmesh)
 
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)
 
     hv_center_space =
diff --git a/test/Operators/hybrid/3d.jl b/test/Operators/hybrid/3d.jl
index 232c8b63b9..d802316951 100644
--- a/test/Operators/hybrid/3d.jl
+++ b/test/Operators/hybrid/3d.jl
@@ -5,12 +5,12 @@ using StaticArrays, IntervalSets, LinearAlgebra
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators
 import ClimaCore.Geometry: WVector
@@ -40,7 +40,7 @@ device = ClimaComms.device()
         ClimaComms.SingletonCommsContext(device),
         horzmesh,
     )
-    quad = Spaces.Quadratures.GLL{3 + 1}()
+    quad = Quadratures.GLL{3 + 1}()
     horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
     hv_center_space =
@@ -87,7 +87,7 @@ function hvspace_3D(
         horzmesh,
     )
 
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
     hv_center_space =
diff --git a/test/Operators/hybrid/cuda.jl b/test/Operators/hybrid/cuda.jl
index 9f97f98ab7..4adb11bc09 100644
--- a/test/Operators/hybrid/cuda.jl
+++ b/test/Operators/hybrid/cuda.jl
@@ -2,7 +2,7 @@ using Test
 using StaticArrays
 using ClimaComms, ClimaCore
 import ClimaCore:
-    Geometry, Fields, Domains, Topologies, Meshes, Spaces, Operators
+    Geometry, Fields, Domains, Topologies, Meshes, Spaces, Operators, Quadratures
 using LinearAlgebra, IntervalSets
 using CUDA
 using OrdinaryDiffEq
@@ -39,7 +39,7 @@ function hvspace_3D_box(
     )
     horzmesh = Meshes.RectilinearMesh(horzdomain, xelem, yelem)
 
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
 
     # Define horz topology and space
     horztopology = Topologies.Topology2D(context, horzmesh)
@@ -70,7 +70,7 @@ function hvspace_3D_sphere(context)
     horzdomain = Domains.SphereDomain(FT(30.0))
     horzmesh = Meshes.EquiangularCubedSphere(horzdomain, 4)
 
-    quad = Spaces.Quadratures.GLL{3 + 1}()
+    quad = Quadratures.GLL{3 + 1}()
 
     # Define horz topology and space
     horztopology = Topologies.Topology2D(context, horzmesh)
diff --git a/test/Operators/hybrid/dss_opt.jl b/test/Operators/hybrid/dss_opt.jl
index c183bae59d..cb859c3b3e 100644
--- a/test/Operators/hybrid/dss_opt.jl
+++ b/test/Operators/hybrid/dss_opt.jl
@@ -7,12 +7,12 @@ import ClimaCore
 import ClimaCore:
     ClimaCore,
     slab,
-    Spaces,
     Domains,
     Meshes,
     Geometry,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     Operators
 
@@ -35,7 +35,7 @@ horztopology = Topologies.Topology2D(
     ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()),
     horzmesh,
 )
-quad = Spaces.Quadratures.GLL{5}()
+quad = Quadratures.GLL{5}()
 horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
 hv_center_space =
diff --git a/test/Operators/hybrid/extruded_3dbox_cuda.jl b/test/Operators/hybrid/extruded_3dbox_cuda.jl
index 1d241bb36c..baa3baa90b 100644
--- a/test/Operators/hybrid/extruded_3dbox_cuda.jl
+++ b/test/Operators/hybrid/extruded_3dbox_cuda.jl
@@ -6,7 +6,7 @@ using LinearAlgebra, IntervalSets
 using CUDA
 using ClimaComms, ClimaCore
 import ClimaCore:
-    Domains, Topologies, Meshes, Spaces, Geometry, column, Fields, Operators
+    Domains, Topologies, Meshes, Spaces, Geometry, column, Fields, Operators, Quadratures
 using Test
 
 function get_space(context)
@@ -32,7 +32,7 @@ function get_space(context)
         x2periodic = true,
     )
     horzmesh = Meshes.RectilinearMesh(horzdomain, 17, 16)
-    quad = Spaces.Quadratures.GLL{3 + 1}()
+    quad = Quadratures.GLL{3 + 1}()
 
     # Define horz topology and space
     horztopology = Topologies.Topology2D(context, horzmesh)
diff --git a/test/Operators/hybrid/extruded_sphere_cuda.jl b/test/Operators/hybrid/extruded_sphere_cuda.jl
index d1afd47a2c..cc8b78d426 100644
--- a/test/Operators/hybrid/extruded_sphere_cuda.jl
+++ b/test/Operators/hybrid/extruded_sphere_cuda.jl
@@ -6,7 +6,7 @@ using LinearAlgebra, IntervalSets
 using CUDA
 using ClimaComms, ClimaCore
 import ClimaCore:
-    Domains, Topologies, Meshes, Spaces, Geometry, column, Fields, Operators
+    Domains, Topologies, Meshes, Spaces, Geometry, column, Fields, Operators, Quadratures
 using Test
 
 function get_space(context)
@@ -28,7 +28,7 @@ function get_space(context)
     horzdomain = Domains.SphereDomain(FT(30.0))
     horzmesh = Meshes.EquiangularCubedSphere(horzdomain, 4)
 
-    quad = Spaces.Quadratures.GLL{3 + 1}()
+    quad = Quadratures.GLL{3 + 1}()
 
     # Define horz topology and space
     horztopology = Topologies.Topology2D(context, horzmesh)
diff --git a/test/Operators/hybrid/opt.jl b/test/Operators/hybrid/opt.jl
index 6fcebff200..86e21681ad 100644
--- a/test/Operators/hybrid/opt.jl
+++ b/test/Operators/hybrid/opt.jl
@@ -6,7 +6,7 @@ using IntervalSets
 
 import ClimaCore
 
-import ClimaCore: Domains, Meshes, Topologies, Spaces, Fields, Operators
+import ClimaCore: Domains, Meshes, Topologies, Spaces, Fields, Operators, Quadratures
 import ClimaCore.Domains: Geometry
 
 # We need to pull these broadcasted expressions out as
@@ -212,7 +212,7 @@ function hspace1d(FT)
     hmesh = Meshes.IntervalMesh(hdomain, nelems = 3)
     htopology = Topologies.IntervalTopology(hmesh)
     Nq = 3
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     return Spaces.SpectralElementSpace1D(htopology, quad)
 end
 
@@ -224,7 +224,7 @@ function hspace2d(FT)
         x2periodic = true,
     )
     Nq = 3
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     hmesh = Meshes.RectilinearMesh(hdomain, 3, 3)
     htopology = Topologies.Topology2D(
         ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()),
diff --git a/test/Operators/remapping.jl b/test/Operators/remapping.jl
index ac7b29ee61..24250578ae 100644
--- a/test/Operators/remapping.jl
+++ b/test/Operators/remapping.jl
@@ -1,10 +1,10 @@
 using Test
 using ClimaComms
 using ClimaCore:
-    Domains, Meshes, Topologies, Geometry, Operators, Spaces, Fields
+    Domains, Meshes, Topologies, Geometry, Operators, Spaces, Fields, Quadratures
 using ClimaCore.Operators: local_weights, LinearRemap, remap, remap!
 using ClimaCore.Topologies: Topology2D
-using ClimaCore.Spaces: AbstractSpace, Quadratures
+using ClimaCore.Spaces: AbstractSpace
 using ClimaCore.DataLayouts: IJFH
 using IntervalSets, LinearAlgebra, SparseArrays
 
@@ -408,7 +408,7 @@ end
             for nq in [2, 5, 10]
                 space2 = make_space(domain, nq)
 
-                _, w = Spaces.Quadratures.quadrature_points(
+                _, w = Quadratures.quadrature_points(
                     FT,
                     Quadratures.GLL{nq}(),
                 )
@@ -515,7 +515,7 @@ end
             for nq1 in [3, 5, 9]
                 space1 = make_space(domain, nq1)
 
-                _, w = Spaces.Quadratures.quadrature_points(
+                _, w = Quadratures.quadrature_points(
                     FT,
                     Quadratures.GLL{nq1}(),
                 )
@@ -524,7 +524,7 @@ end
                 for nq2 in [2, 4, 10]
                     space2 = make_space(domain, nq2)
 
-                    _, w = Spaces.Quadratures.quadrature_points(
+                    _, w = Quadratures.quadrature_points(
                         FT,
                         Quadratures.GLL{nq2}(),
                     )
diff --git a/test/Operators/spectralelement/benchmark_utils.jl b/test/Operators/spectralelement/benchmark_utils.jl
index ae522f9bd3..bb222eeb49 100644
--- a/test/Operators/spectralelement/benchmark_utils.jl
+++ b/test/Operators/spectralelement/benchmark_utils.jl
@@ -12,6 +12,7 @@ import ClimaCore.Meshes as Meshes
 import ClimaCore.Spaces as Spaces
 import ClimaCore.Topologies as Topologies
 import ClimaCore.Geometry as Geometry
+import ClimaCore.Quadratures as Quadratures
 import ClimaCore.Fields as Fields
 import ClimaCore as CC
 import ClimaComms
@@ -118,7 +119,7 @@ function create_space(
     hdomain = Domains.SphereDomain(earth_radius)
     hmesh = Meshes.EquiangularCubedSphere(hdomain, panel_size)
     htopology = Topologies.Topology2D(context, hmesh)
-    quad = Spaces.Quadratures.GLL{poly_nodes}()
+    quad = Quadratures.GLL{poly_nodes}()
     space = if space_type == "SpectralElementSpace2D"
         Spaces.SpectralElementSpace2D(htopology, quad)
     elseif space_type == "ExtrudedFiniteDifferenceSpace"
diff --git a/test/Operators/spectralelement/diffusion2d.jl b/test/Operators/spectralelement/diffusion2d.jl
index 601cd76b31..2fdb49fa93 100644
--- a/test/Operators/spectralelement/diffusion2d.jl
+++ b/test/Operators/spectralelement/diffusion2d.jl
@@ -3,7 +3,7 @@ using ClimaComms
 using StaticArrays, IntervalSets
 import ClimaCore.DataLayouts: IJFH
 import ClimaCore:
-    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
+    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry, Quadratures
 using StaticArrays, IntervalSets, LinearAlgebra
 
 using OrdinaryDiffEq
@@ -31,8 +31,8 @@ using OrdinaryDiffEq
     )
 
     Nq = 6
-    quad = Spaces.Quadratures.GLL{Nq}()
-    points, weights = Spaces.Quadratures.quadrature_points(Float64, quad)
+    quad = Quadratures.GLL{Nq}()
+    points, weights = Quadratures.quadrature_points(Float64, quad)
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
     # Define eigensolution
diff --git a/test/Operators/spectralelement/opt.jl b/test/Operators/spectralelement/opt.jl
index 90fec6dc23..6ebc0d7aec 100644
--- a/test/Operators/spectralelement/opt.jl
+++ b/test/Operators/spectralelement/opt.jl
@@ -4,7 +4,7 @@ using ClimaComms
 using LinearAlgebra, IntervalSets
 
 import ClimaCore:
-    Geometry, Fields, Domains, Topologies, Meshes, Spaces, Operators
+    Geometry, Fields, Domains, Topologies, Meshes, Spaces, Operators, Quadratures
 
 
 # We need to pull these broadcasted expressions out as
@@ -137,7 +137,7 @@ end
             )
 
             Nq = 3
-            quad = Spaces.Quadratures.GLL{Nq}()
+            quad = Quadratures.GLL{Nq}()
             mesh = Meshes.RectilinearMesh(domain, 3, 3)
 
             topology = Topologies.Topology2D(
@@ -158,7 +158,7 @@ end
                 )
 
             INq = 4
-            Iquad = Spaces.Quadratures.GLL{INq}()
+            Iquad = Quadratures.GLL{INq}()
             Ispace =
                 Spaces.SpectralElementSpace2D(Spaces.topology(space), Iquad)
 
@@ -206,7 +206,7 @@ end
                 horzmesh,
             )
 
-            quad = Spaces.Quadratures.GLL{npoly + 1}()
+            quad = Quadratures.GLL{npoly + 1}()
             horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
             hv_center_space = Spaces.ExtrudedFiniteDifferenceSpace(
diff --git a/test/Operators/spectralelement/plane.jl b/test/Operators/spectralelement/plane.jl
index bb65f029d5..b12d05d1e9 100644
--- a/test/Operators/spectralelement/plane.jl
+++ b/test/Operators/spectralelement/plane.jl
@@ -3,7 +3,7 @@ using StaticArrays
 using ClimaComms
 import ClimaCore.DataLayouts: IJFH, VF
 import ClimaCore:
-    Geometry, Fields, Domains, Topologies, Meshes, Spaces, Operators
+    Geometry, Fields, Domains, Topologies, Meshes, Spaces, Operators, Quadratures
 using LinearAlgebra, IntervalSets
 
 FT = Float64
@@ -13,7 +13,7 @@ hdomain = Domains.IntervalDomain(
 )
 
 Nq = 5
-quad = Spaces.Quadratures.GLL{Nq}()
+quad = Quadratures.GLL{Nq}()
 device = ClimaComms.CPUSingleThreaded()
 hmesh = Meshes.IntervalMesh(hdomain, nelems = 16)
 htopology =
diff --git a/test/Operators/spectralelement/rectilinear.jl b/test/Operators/spectralelement/rectilinear.jl
index 935412c208..679ead9cae 100644
--- a/test/Operators/spectralelement/rectilinear.jl
+++ b/test/Operators/spectralelement/rectilinear.jl
@@ -3,7 +3,7 @@ using StaticArrays
 using ClimaComms
 import ClimaCore.DataLayouts: IJFH, VF
 import ClimaCore:
-    Geometry, Fields, Domains, Topologies, Meshes, Spaces, Operators
+    Geometry, Fields, Domains, Topologies, Meshes, Spaces, Operators, Quadratures
 using LinearAlgebra, IntervalSets
 
 FT = Float64
@@ -15,7 +15,7 @@ domain = Domains.RectangleDomain(
 )
 
 Nq = 5
-quad = Spaces.Quadratures.GLL{Nq}()
+quad = Quadratures.GLL{Nq}()
 device = ClimaComms.CPUSingleThreaded()
 grid_mesh = Meshes.RectilinearMesh(domain, 17, 16)
 grid_topology =
@@ -36,7 +36,7 @@ ts_test_setup = (ts_topology, ts_space, ts_coords)
 
     for (topology, space, coords) in (grid_test_setup, ts_test_setup)
         INq = 9
-        Iquad = Spaces.Quadratures.GLL{INq}()
+        Iquad = Quadratures.GLL{INq}()
         Ispace = Spaces.SpectralElementSpace2D(topology, Iquad)
 
         I = Operators.Interpolate(Ispace)
diff --git a/test/Operators/spectralelement/rectilinear_cuda.jl b/test/Operators/spectralelement/rectilinear_cuda.jl
index 1de296a85c..116499504b 100644
--- a/test/Operators/spectralelement/rectilinear_cuda.jl
+++ b/test/Operators/spectralelement/rectilinear_cuda.jl
@@ -2,7 +2,7 @@ using Test
 using StaticArrays
 using ClimaComms, ClimaCore
 import ClimaCore:
-    Geometry, Fields, Domains, Topologies, Meshes, Spaces, Operators
+    Geometry, Fields, Domains, Topologies, Meshes, Spaces, Operators, Quadratures
 using LinearAlgebra, IntervalSets
 using CUDA
 
@@ -15,7 +15,7 @@ domain = Domains.RectangleDomain(
 )
 
 Nq = 5
-quad = Spaces.Quadratures.GLL{Nq}()
+quad = Quadratures.GLL{Nq}()
 
 grid_mesh = Meshes.RectilinearMesh(domain, 17, 16)
 
diff --git a/test/Operators/spectralelement/sphere_curl.jl b/test/Operators/spectralelement/sphere_curl.jl
index a18a9b5542..7d4a191f59 100644
--- a/test/Operators/spectralelement/sphere_curl.jl
+++ b/test/Operators/spectralelement/sphere_curl.jl
@@ -3,7 +3,7 @@ using ClimaComms
 using StaticArrays, IntervalSets
 import ClimaCore.DataLayouts: IJFH
 import ClimaCore:
-    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
+    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry, Quadratures
 using StaticArrays, IntervalSets, LinearAlgebra
 
 FT = Float64
@@ -24,7 +24,7 @@ wcurl = Operators.WeakCurl()
         mesh,
     )
 
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
     coords = Fields.coordinate_field(space)
 
@@ -82,7 +82,7 @@ convergence_rate(err, Δh) =
                 mesh,
             )
 
-            quad = Spaces.Quadratures.GLL{Nq}()
+            quad = Quadratures.GLL{Nq}()
             space = Spaces.SpectralElementSpace2D(grid_topology, quad)
             coords = Fields.coordinate_field(space)
 
diff --git a/test/Operators/spectralelement/sphere_diffusion.jl b/test/Operators/spectralelement/sphere_diffusion.jl
index 304d4aae29..ea2b40c097 100644
--- a/test/Operators/spectralelement/sphere_diffusion.jl
+++ b/test/Operators/spectralelement/sphere_diffusion.jl
@@ -3,7 +3,7 @@ using StaticArrays, IntervalSets
 using ClimaComms
 import ClimaCore.DataLayouts: IJFH
 import ClimaCore:
-    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
+    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry, Quadratures
 using StaticArrays, IntervalSets, LinearAlgebra
 
 @testset "Scalar Poisson problem - ∇⋅∇ = f on the cubed-sphere" begin
@@ -27,7 +27,7 @@ using StaticArrays, IntervalSets, LinearAlgebra
         ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()),
         mesh,
     )
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
     # Define eigensolution
diff --git a/test/Operators/spectralelement/sphere_diffusion_vec.jl b/test/Operators/spectralelement/sphere_diffusion_vec.jl
index dcaeae5306..4fe66248d3 100644
--- a/test/Operators/spectralelement/sphere_diffusion_vec.jl
+++ b/test/Operators/spectralelement/sphere_diffusion_vec.jl
@@ -3,7 +3,7 @@ using ClimaComms
 using StaticArrays, IntervalSets
 import ClimaCore.DataLayouts: IJFH
 import ClimaCore:
-    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
+    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry, Quadratures
 using StaticArrays, IntervalSets, LinearAlgebra
 
 include("sphere_sphericalharmonics.jl")
@@ -40,7 +40,7 @@ end
         mesh,
     )
 
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
     coords = Fields.coordinate_field(space)
 
@@ -88,7 +88,7 @@ convergence_rate(err, Δh) =
                 mesh,
             )
 
-            quad = Spaces.Quadratures.GLL{Nq}()
+            quad = Quadratures.GLL{Nq}()
             space = Spaces.SpectralElementSpace2D(grid_topology, quad)
             coords = Fields.coordinate_field(space)
 
diff --git a/test/Operators/spectralelement/sphere_divergence.jl b/test/Operators/spectralelement/sphere_divergence.jl
index 352fe5c840..f506193389 100644
--- a/test/Operators/spectralelement/sphere_divergence.jl
+++ b/test/Operators/spectralelement/sphere_divergence.jl
@@ -3,7 +3,7 @@ using ClimaComms
 using StaticArrays, IntervalSets
 import ClimaCore.DataLayouts: IJFH
 import ClimaCore:
-    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
+    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry, Quadratures
 using StaticArrays, IntervalSets, LinearAlgebra
 
 FT = Float64
@@ -24,7 +24,7 @@ wdiv = Operators.WeakDivergence()
         mesh,
     )
 
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
     coords = Fields.coordinate_field(space)
 
@@ -78,7 +78,7 @@ convergence_rate(err, Δh) =
                 mesh,
             )
 
-            quad = Spaces.Quadratures.GLL{Nq}()
+            quad = Quadratures.GLL{Nq}()
             space = Spaces.SpectralElementSpace2D(grid_topology, quad)
             coords = Fields.coordinate_field(space)
 
diff --git a/test/Operators/spectralelement/sphere_geometry.jl b/test/Operators/spectralelement/sphere_geometry.jl
index 98b2404911..fab597ea2f 100644
--- a/test/Operators/spectralelement/sphere_geometry.jl
+++ b/test/Operators/spectralelement/sphere_geometry.jl
@@ -1,7 +1,7 @@
 using LinearAlgebra, IntervalSets
 using ClimaComms
 import ClimaCore:
-    Domains, Topologies, Meshes, Spaces, Geometry, Operators, Fields
+    Domains, Topologies, Meshes, Spaces, Geometry, Operators, Fields, Quadratures
 
 using Test
 
@@ -31,7 +31,7 @@ end
             ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()),
             mesh,
         )
-        quad = Spaces.Quadratures.GLL{Nq}()
+        quad = Quadratures.GLL{Nq}()
         space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
         @test sum(ones(space)) ≈ 4pi * radius^2 rtol = 1e-3
diff --git a/test/Operators/spectralelement/sphere_geometry_distributed.jl b/test/Operators/spectralelement/sphere_geometry_distributed.jl
index 82193a5c50..e9c57c3f2d 100644
--- a/test/Operators/spectralelement/sphere_geometry_distributed.jl
+++ b/test/Operators/spectralelement/sphere_geometry_distributed.jl
@@ -1,6 +1,6 @@
 using LinearAlgebra, IntervalSets
 import ClimaCore:
-    Domains, Topologies, Meshes, Spaces, Geometry, Operators, Fields
+    Domains, Topologies, Meshes, Spaces, Geometry, Operators, Fields, Quadratures
 
 using Test
 
@@ -46,7 +46,7 @@ end
             mesh,
             Topologies.spacefillingcurve(mesh),
         )
-        quad = Spaces.Quadratures.GLL{Nq}()
+        quad = Quadratures.GLL{Nq}()
         space = Spaces.SpectralElementSpace2D(grid_topology, quad)
         # for comparison with serial results
         grid_topology_serial =
diff --git a/test/Operators/spectralelement/sphere_gradient.jl b/test/Operators/spectralelement/sphere_gradient.jl
index af9db9052c..675fe0b8a6 100644
--- a/test/Operators/spectralelement/sphere_gradient.jl
+++ b/test/Operators/spectralelement/sphere_gradient.jl
@@ -3,7 +3,7 @@ using ClimaComms
 using StaticArrays, IntervalSets
 import ClimaCore.DataLayouts: IJFH
 import ClimaCore:
-    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
+    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry, Quadratures
 using StaticArrays, IntervalSets, LinearAlgebra
 
 FT = Float64
@@ -24,7 +24,7 @@ wgrad = Operators.WeakGradient()
         mesh,
     )
 
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
     coords = Fields.coordinate_field(space)
 
@@ -77,7 +77,7 @@ convergence_rate(err, Δh) =
                 mesh,
             )
 
-            quad = Spaces.Quadratures.GLL{Nq}()
+            quad = Quadratures.GLL{Nq}()
             space = Spaces.SpectralElementSpace2D(grid_topology, quad)
             coords = Fields.coordinate_field(space)
 
diff --git a/test/Operators/spectralelement/sphere_hyperdiffusion.jl b/test/Operators/spectralelement/sphere_hyperdiffusion.jl
index 1985977f88..7fb365d16f 100644
--- a/test/Operators/spectralelement/sphere_hyperdiffusion.jl
+++ b/test/Operators/spectralelement/sphere_hyperdiffusion.jl
@@ -3,7 +3,7 @@ using ClimaComms
 using StaticArrays, IntervalSets
 import ClimaCore.DataLayouts: IJFH
 import ClimaCore:
-    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
+    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry, Quadratures
 using StaticArrays, IntervalSets, LinearAlgebra
 
 using OrdinaryDiffEq
@@ -30,7 +30,7 @@ include("sphere_sphericalharmonics.jl")
     )
 
     Nq = 6
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
     coords = Fields.coordinate_field(space)
 
diff --git a/test/Operators/spectralelement/sphere_hyperdiffusion_vec.jl b/test/Operators/spectralelement/sphere_hyperdiffusion_vec.jl
index ea91a173ca..c012f64ee2 100644
--- a/test/Operators/spectralelement/sphere_hyperdiffusion_vec.jl
+++ b/test/Operators/spectralelement/sphere_hyperdiffusion_vec.jl
@@ -3,7 +3,7 @@ using ClimaComms
 using StaticArrays, IntervalSets
 import ClimaCore.DataLayouts: IJFH
 import ClimaCore:
-    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
+    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry, Quadratures
 using StaticArrays, IntervalSets, LinearAlgebra
 
 include("sphere_sphericalharmonics.jl")
@@ -46,7 +46,7 @@ end
         mesh,
     )
 
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
     coords = Fields.coordinate_field(space)
 
@@ -94,7 +94,7 @@ convergence_rate(err, Δh) =
                 mesh,
             )
 
-            quad = Spaces.Quadratures.GLL{Nq}()
+            quad = Quadratures.GLL{Nq}()
             space = Spaces.SpectralElementSpace2D(grid_topology, quad)
             coords = Fields.coordinate_field(space)
 
diff --git a/test/Remapping/distributed_remapping.jl b/test/Remapping/distributed_remapping.jl
index 92dff5db08..a49563ffbb 100644
--- a/test/Remapping/distributed_remapping.jl
+++ b/test/Remapping/distributed_remapping.jl
@@ -9,6 +9,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     Remapping,
     Hypsography
@@ -44,7 +45,7 @@ if !(device isa ClimaComms.CUDADevice)
             periodic = true,
         )
 
-        quad = Spaces.Quadratures.GLL{4}()
+        quad = Quadratures.GLL{4}()
         horzmesh = Meshes.IntervalMesh(horzdomain, nelems = 10)
         horztopology = Topologies.IntervalTopology(
             ClimaComms.SingletonCommsContext(device),
@@ -140,7 +141,7 @@ end
         x2periodic = true,
     )
 
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     horzmesh = Meshes.RectilinearMesh(horzdomain, 10, 10)
     horztopology = Topologies.Topology2D(
         ClimaComms.SingletonCommsContext(device),
@@ -276,7 +277,7 @@ end
         x2periodic = true,
     )
 
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     horzmesh = Meshes.RectilinearMesh(horzdomain, 10, 10)
     horztopology = Topologies.Topology2D(
         ClimaComms.SingletonCommsContext(device),
@@ -408,7 +409,7 @@ end
 
     horzdomain = Domains.SphereDomain(1e6)
 
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     horzmesh = Meshes.EquiangularCubedSphere(horzdomain, 6)
     horztopology = Topologies.Topology2D(context, horzmesh)
     horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
diff --git a/test/Remapping/interpolate_array.jl b/test/Remapping/interpolate_array.jl
index 0dd8398eaa..303c9bb11c 100644
--- a/test/Remapping/interpolate_array.jl
+++ b/test/Remapping/interpolate_array.jl
@@ -1,7 +1,7 @@
 
 using ClimaComms
 using ClimaCore:
-    Geometry, Domains, Meshes, Topologies, Spaces, Fields, Remapping
+    Geometry, Domains, Meshes, Topologies, Spaces, Fields, Remapping, Quadratures
 using IntervalSets
 using Test
 
@@ -26,7 +26,7 @@ device = ClimaComms.CPUSingleThreaded()
         periodic = true,
     )
 
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     horzmesh = Meshes.IntervalMesh(horzdomain, nelems = 10)
     horztopology = Topologies.IntervalTopology(
         ClimaComms.SingletonCommsContext(device),
@@ -76,7 +76,7 @@ end
         x2periodic = true,
     )
 
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     horzmesh = Meshes.RectilinearMesh(horzdomain, 10, 10)
     horztopology = Topologies.Topology2D(
         ClimaComms.SingletonCommsContext(device),
@@ -131,7 +131,7 @@ end
         x2periodic = true,
     )
 
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     horzmesh = Meshes.RectilinearMesh(horzdomain, 10, 10)
     horztopology = Topologies.Topology2D(
         ClimaComms.SingletonCommsContext(device),
@@ -183,7 +183,7 @@ end
 
     horzdomain = Domains.SphereDomain(1e6)
 
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     horzmesh = Meshes.EquiangularCubedSphere(horzdomain, 6)
     horztopology = Topologies.Topology2D(
         ClimaComms.SingletonCommsContext(device),
diff --git a/test/Spaces/ddss1.jl b/test/Spaces/ddss1.jl
index f6d61e81d3..2a25994dba 100644
--- a/test/Spaces/ddss1.jl
+++ b/test/Spaces/ddss1.jl
@@ -10,6 +10,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     DataLayouts
 
@@ -42,7 +43,7 @@ function distributed_space(
     )
     mesh = Meshes.RectilinearMesh(domain, n1, n2)
     topology = Topologies.Topology2D(context, mesh, Meshes.elements(mesh))
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(topology, quad)
 
     return (space, context)
diff --git a/test/Spaces/ddss1_cs.jl b/test/Spaces/ddss1_cs.jl
index dbf68d69df..ea90819fcd 100644
--- a/test/Spaces/ddss1_cs.jl
+++ b/test/Spaces/ddss1_cs.jl
@@ -7,6 +7,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     DataLayouts
 
@@ -22,7 +23,7 @@ import ClimaCore:
     mesh = Meshes.EquiangularCubedSphere(domain, 3)
     topology = Topologies.Topology2D(context, mesh)
     topology_cpu = Topologies.Topology2D(context_cpu, mesh)
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     space = Spaces.SpectralElementSpace2D(topology, quad)
     space_cpu = Spaces.SpectralElementSpace2D(topology_cpu, quad)
 
@@ -82,7 +83,7 @@ end
         horizontal_mesh,
         Topologies.spacefillingcurve(horizontal_mesh),
     )
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     h_space = Spaces.SpectralElementSpace2D(horizontal_topology, quad)
     h_space_cpu = Spaces.SpectralElementSpace2D(horizontal_topology_cpu, quad)
 
diff --git a/test/Spaces/distributed/ddss_setup.jl b/test/Spaces/distributed/ddss_setup.jl
index 590f7fe434..48cf087687 100644
--- a/test/Spaces/distributed/ddss_setup.jl
+++ b/test/Spaces/distributed/ddss_setup.jl
@@ -2,7 +2,7 @@ using Logging
 using Test
 
 import ClimaCore:
-    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies
+    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Quadratures
 
 using ClimaComms
 const context = ClimaComms.MPICommsContext()
@@ -40,7 +40,7 @@ function distributed_space(
     )
     mesh = Meshes.RectilinearMesh(domain, n1, n2)
     topology = Topologies.Topology2D(context, mesh, Meshes.elements(mesh))
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(topology, quad)
 
     return (space, context)
diff --git a/test/Spaces/distributed/gather4.jl b/test/Spaces/distributed/gather4.jl
index 09b99791ce..9e0cb445ab 100644
--- a/test/Spaces/distributed/gather4.jl
+++ b/test/Spaces/distributed/gather4.jl
@@ -8,6 +8,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     DataLayouts
 
@@ -39,7 +40,7 @@ domain = Domains.RectangleDomain(
 )
 n1, n2 = 16, 16
 Nq = 4
-quad = Spaces.Quadratures.GLL{Nq}()
+quad = Quadratures.GLL{Nq}()
 mesh = Meshes.RectilinearMesh(domain, n1, n2)
 
 grid_topology = Topologies.Topology2D(comms_ctx, mesh, Meshes.elements(mesh))
diff --git a/test/Spaces/distributed_cuda/ddss2.jl b/test/Spaces/distributed_cuda/ddss2.jl
index b41ec33830..b7ac329a03 100644
--- a/test/Spaces/distributed_cuda/ddss2.jl
+++ b/test/Spaces/distributed_cuda/ddss2.jl
@@ -4,7 +4,7 @@ using Logging
 using Test
 
 import ClimaCore:
-    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies
+    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Quadratures
 
 using ClimaComms
 using CUDA
@@ -55,7 +55,7 @@ pid, nprocs = ClimaComms.init(context)
     )
     mesh = Meshes.RectilinearMesh(domain, n1, n2)
     topology = Topologies.Topology2D(context, mesh, Meshes.elements(mesh))
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(topology, quad)
 
 
diff --git a/test/Spaces/distributed_cuda/ddss3.jl b/test/Spaces/distributed_cuda/ddss3.jl
index 67a90c87be..e32cc63ac3 100644
--- a/test/Spaces/distributed_cuda/ddss3.jl
+++ b/test/Spaces/distributed_cuda/ddss3.jl
@@ -4,7 +4,7 @@ using Logging
 using Test
 
 import ClimaCore:
-    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies
+    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Quadratures
 
 using ClimaComms
 using CUDA
@@ -76,7 +76,7 @@ partition numbers
     )
     mesh = Meshes.RectilinearMesh(domain, n1, n2)
     topology = Topologies.Topology2D(context, mesh, Meshes.elements(mesh))
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(topology, quad)
 
     @test Topologies.nlocalelems(Spaces.topology(space)) == (pid == 1 ? 6 : 5)
diff --git a/test/Spaces/distributed_cuda/ddss4.jl b/test/Spaces/distributed_cuda/ddss4.jl
index 2ff2c3dcf9..ba3c21c027 100644
--- a/test/Spaces/distributed_cuda/ddss4.jl
+++ b/test/Spaces/distributed_cuda/ddss4.jl
@@ -2,7 +2,7 @@ using Logging
 using Test
 
 import ClimaCore:
-    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies
+    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Quadratures
 
 using ClimaComms
 using CUDA
@@ -49,7 +49,7 @@ pid, nprocs = ClimaComms.init(context)
     )
     mesh = Meshes.RectilinearMesh(domain, n1, n2)
     topology = Topologies.Topology2D(context, mesh, Meshes.elements(mesh))
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     space = Spaces.SpectralElementSpace2D(topology, quad)
     init_state(local_geometry, p) = (ρ = 1.0)
     y0 = init_state.(Fields.local_geometry_field(space), Ref(nothing))
diff --git a/test/Spaces/distributed_cuda/ddss_ne32_cs.jl b/test/Spaces/distributed_cuda/ddss_ne32_cs.jl
index 57f9037445..20287a2dbf 100644
--- a/test/Spaces/distributed_cuda/ddss_ne32_cs.jl
+++ b/test/Spaces/distributed_cuda/ddss_ne32_cs.jl
@@ -8,6 +8,7 @@ import ClimaCore:
     Meshes,
     Operators,
     Spaces,
+    Quadratures,
     Topologies,
     DataLayouts
 
@@ -30,7 +31,7 @@ import ClimaCore:
     mesh = Meshes.EquiangularCubedSphere(domain, 32)
     topology_cuda = Topologies.Topology2D(context_cuda, mesh)
     topology_cpu = Topologies.Topology2D(context_cpu, mesh)
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     space_cuda = Spaces.SpectralElementSpace2D(topology_cuda, quad)
     space_cpu = Spaces.SpectralElementSpace2D(topology_cpu, quad)
     x_cuda = ones(space_cuda)
diff --git a/test/Spaces/distributed_cuda/space_construction.jl b/test/Spaces/distributed_cuda/space_construction.jl
index b421f5556d..3d5a2b58fa 100644
--- a/test/Spaces/distributed_cuda/space_construction.jl
+++ b/test/Spaces/distributed_cuda/space_construction.jl
@@ -4,7 +4,7 @@ using ClimaComms
 using CUDA
 
 import ClimaCore:
-    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies
+    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Quadratures
 
 @testset "Distributed extruded CUDA space" begin
     # initializing MPI
@@ -35,7 +35,7 @@ import ClimaCore:
     hdomain = Domains.SphereDomain(radius)
     hmesh = Meshes.EquiangularCubedSphere(hdomain, helem)
     htopology = Topologies.Topology2D(context, hmesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     hspace = Spaces.SpectralElementSpace2D(htopology, quad)
     space = Spaces.ExtrudedFiniteDifferenceSpace(hspace, vspace)
 
diff --git a/test/Spaces/spaces.jl b/test/Spaces/spaces.jl
index 910c7478ee..eec06779fb 100644
--- a/test/Spaces/spaces.jl
+++ b/test/Spaces/spaces.jl
@@ -9,6 +9,7 @@ import ClimaCore:
     Meshes,
     Topologies,
     Spaces,
+    Quadratures,
     Fields,
     DataLayouts,
     Geometry,
@@ -28,8 +29,8 @@ on_gpu = ClimaComms.device() isa ClimaComms.CUDADevice
     mesh = Meshes.IntervalMesh(domain; nelems = 1)
     topology = Topologies.IntervalTopology(mesh)
 
-    quad = Spaces.Quadratures.GLL{4}()
-    points, weights = Spaces.Quadratures.quadrature_points(FT, quad)
+    quad = Quadratures.GLL{4}()
+    points, weights = Quadratures.quadrature_points(FT, quad)
 
     space = Spaces.SpectralElementSpace1D(topology, quad)
 
@@ -95,7 +96,7 @@ on_gpu || @testset "extruded (2d 1×3) finite difference space" begin
     )
     horzmesh = Meshes.IntervalMesh(horzdomain; nelems = 5)
     horztopology = Topologies.IntervalTopology(horzmesh)
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
 
     hspace = Spaces.SpectralElementSpace1D(horztopology, quad)
     # Extrusion
@@ -154,7 +155,7 @@ end
     domain = Domains.RectangleDomain(x_domain, y_domain)
     hmesh = Meshes.RectilinearMesh(domain, x_elem, y_elem)
 
-    quad = Spaces.Quadratures.GL{1}()
+    quad = Quadratures.GL{1}()
     htopology = Topologies.Topology2D(context, hmesh)
     hspace = Spaces.SpectralElementSpace2D(htopology, quad)
 
@@ -188,8 +189,8 @@ end
     mesh = Meshes.RectilinearMesh(domain, 1, 1)
     grid_topology = Topologies.Topology2D(context, mesh)
 
-    quad = Spaces.Quadratures.GLL{4}()
-    points, weights = Spaces.Quadratures.quadrature_points(FT, quad)
+    quad = Quadratures.GLL{4}()
+    points, weights = Quadratures.quadrature_points(FT, quad)
 
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
     @test repr(space) == """
@@ -259,7 +260,7 @@ end
     )
     n1, n2 = 2, 2
     Nq = 5
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     mesh = Meshes.RectilinearMesh(domain, n1, n2)
     grid_topology = Topologies.Topology2D(context, mesh)
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
@@ -304,8 +305,8 @@ end
     mesh = Meshes.RectilinearMesh(domain, n1, n2)
     grid_topology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), mesh)
 
-    quad = Spaces.Quadratures.GLL{4}()
-    points, weights = Spaces.Quadratures.quadrature_points(FT, quad)
+    quad = Quadratures.GLL{4}()
+    points, weights = Quadratures.quadrature_points(FT, quad)
 
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
@@ -384,8 +385,8 @@ end
     mesh = Meshes.RectilinearMesh(domain, n1, n2)
     grid_topology = Topologies.Topology2D(ClimaComms.SingletonCommsContext(), mesh)
 
-    quad = Spaces.Quadratures.GLL{Nij}()
-    points, weights = Spaces.Quadratures.quadrature_points(FT, quad)
+    quad = Quadratures.GLL{Nij}()
+    points, weights = Quadratures.quadrature_points(FT, quad)
 
     space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
diff --git a/test/Spaces/sphere.jl b/test/Spaces/sphere.jl
index 02632830be..0716088dfc 100644
--- a/test/Spaces/sphere.jl
+++ b/test/Spaces/sphere.jl
@@ -1,6 +1,6 @@
 using LinearAlgebra, IntervalSets
 using ClimaComms
-import ClimaCore: Domains, Topologies, Meshes, Spaces, Geometry, column
+import ClimaCore: Domains, Topologies, Meshes, Spaces, Geometry, column, Quadratures
 
 using Test
 
@@ -13,7 +13,7 @@ using Test
         domain = Domains.SphereDomain(radius)
         mesh = Meshes.EquiangularCubedSphere(domain, ne)
         topology = Topologies.Topology2D(context, mesh)
-        quad = Spaces.Quadratures.GLL{Nq}()
+        quad = Quadratures.GLL{Nq}()
         space = Spaces.SpectralElementSpace2D(topology, quad)
 
         # surface area
@@ -66,7 +66,7 @@ end
             domain = Domains.SphereDomain(radius)
             mesh = Meshes.EquiangularCubedSphere(domain, ne)
             topology = Topologies.Topology2D(context, mesh)
-            quad = Spaces.Quadratures.GLL{Nq}()
+            quad = Quadratures.GLL{Nq}()
             no_bubble_space = Spaces.SpectralElementSpace2D(topology, quad)
             # check surface area
             @test sum(ones(no_bubble_space)) ≈ FT(4pi * radius^2) rtol =
@@ -103,7 +103,7 @@ end
     horzdomain = Domains.SphereDomain(radius)
     horzmesh = Meshes.EquiangularCubedSphere(horzdomain, helem)
     horztopology = Topologies.Topology2D(context, horzmesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
     hv_center_space =
diff --git a/test/Spaces/terrain_warp.jl b/test/Spaces/terrain_warp.jl
index 7706245e83..21e46205d4 100644
--- a/test/Spaces/terrain_warp.jl
+++ b/test/Spaces/terrain_warp.jl
@@ -10,6 +10,7 @@ import ClimaCore:
     Operators,
     Meshes,
     Spaces,
+    Quadratures,
     Topologies,
     Hypsography
 
@@ -58,7 +59,7 @@ function generate_base_spaces(
     vert_face_space = Spaces.FaceFiniteDifferenceSpace(vertmesh)
 
     # Generate Horizontal Space
-    quad = Spaces.Quadratures.GLL{npoly + 1}()
+    quad = Quadratures.GLL{npoly + 1}()
     if ndims == 2
         horzdomain = Domains.IntervalDomain(
             Geometry.XPoint{FT}(xlim[1]),
@@ -437,7 +438,7 @@ end
             horzmesh = Meshes.IntervalMesh(horzdomain, nelems = nh)
             horztopology = Topologies.IntervalTopology(horzmesh)
 
-            quad = Spaces.Quadratures.GLL{np + 1}()
+            quad = Quadratures.GLL{np + 1}()
             hspace = Spaces.SpectralElementSpace1D(horztopology, quad)
 
             # Generate surface elevation profile
@@ -491,7 +492,7 @@ end
             horzmesh = Meshes.IntervalMesh(horzdomain, nelems = nh)
             horztopology = Topologies.IntervalTopology(horzmesh)
 
-            quad = Spaces.Quadratures.GLL{np + 1}()
+            quad = Quadratures.GLL{np + 1}()
             hspace = Spaces.SpectralElementSpace1D(horztopology, quad)
 
             # Generate surface elevation profile
diff --git a/test/TestUtilities/TestUtilities.jl b/test/TestUtilities/TestUtilities.jl
index e80e12a650..720393c3f0 100644
--- a/test/TestUtilities/TestUtilities.jl
+++ b/test/TestUtilities/TestUtilities.jl
@@ -13,6 +13,7 @@ using IntervalSets
 import ClimaComms
 import ClimaCore.Fields as Fields
 import ClimaCore.Utilities as Utilities
+import ClimaCore.Quadratures
 import ClimaCore.Geometry as Geometry
 import ClimaCore.Meshes as Meshes
 import ClimaCore.Spaces as Spaces
@@ -39,7 +40,7 @@ function SpectralElementSpace1D(
     )
     mesh = Meshes.IntervalMesh(domain; nelems = 1)
     topology = Topologies.IntervalTopology(context, mesh)
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     return Spaces.SpectralElementSpace1D(topology, quad)
 end
 
@@ -57,7 +58,7 @@ function SpectralElementSpace2D(
     )
     mesh = Meshes.RectilinearMesh(domain, 1, 1)
     topology = Topologies.Topology2D(context, mesh)
-    quad = Spaces.Quadratures.GLL{4}()
+    quad = Quadratures.GLL{4}()
     return Spaces.SpectralElementSpace2D(topology, quad)
 end
 
@@ -98,7 +99,7 @@ function SphereSpectralElementSpace(
     domain = Domains.SphereDomain(radius)
     mesh = Meshes.EquiangularCubedSphere(domain, ne)
     topology = Topologies.Topology2D(context, mesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     return Spaces.SpectralElementSpace2D(topology, quad)
 end
 
@@ -123,7 +124,7 @@ function CenterExtrudedFiniteDifferenceSpace(
     hdomain = Domains.SphereDomain(radius)
     hmesh = Meshes.EquiangularCubedSphere(hdomain, helem)
     htopology = Topologies.Topology2D(context, hmesh)
-    quad = Spaces.Quadratures.GLL{Nq}()
+    quad = Quadratures.GLL{Nq}()
     hspace = Spaces.SpectralElementSpace2D(htopology, quad)
     return Spaces.ExtrudedFiniteDifferenceSpace(hspace, vspace)
 end