diff --git a/benchmark/benchmark_static_ocean.jl b/benchmark/benchmark_static_ocean.jl
index 3800e83307..88c920cda7 100644
--- a/benchmark/benchmark_static_ocean.jl
+++ b/benchmark/benchmark_static_ocean.jl
@@ -19,7 +19,7 @@ for arch in archs, float_type in float_types, N in Ns
     Nx, Ny, Nz = N
     Lx, Ly, Lz = 1, 1, 1
 
-    model = Model(architecture=arch, float_type=float_type, grid=RegularCartesianGrid(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz)))
+    model = Model(architecture=arch, float_type=float_type, grid=RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)))
     time_step!(model, Ni, 1)
 
     bname =  benchmark_name(N, "", arch, float_type)
diff --git a/benchmark/benchmark_tracers.jl b/benchmark/benchmark_tracers.jl
index db0e990d75..b28be2a528 100644
--- a/benchmark/benchmark_tracers.jl
+++ b/benchmark/benchmark_tracers.jl
@@ -57,10 +57,10 @@ for arch in archs, test_case in test_cases
 
     bname =  benchmark_name(N, "$na active + $(lpad(np, 2)) passive", arch, FT)
     @printf("Running benchmark: %s...\n", bname)
-    
+
     model = Model(architecture = arch,
                     float_type = FT,
-                          grid = RegularCartesianGrid(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz)),
+                          grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
                       buoyancy = na2buoyancy(na),
                        tracers = tracers)
     time_step!(model, Ni, 1)
@@ -72,4 +72,3 @@ end
 
 print_timer(timer, title="Tracer benchmarks", sortby=:name)
 println("")
-
diff --git a/examples/internal_wave.jl b/examples/internal_wave.jl
index cd39fa1422..40fd03bc43 100644
--- a/examples/internal_wave.jl
+++ b/examples/internal_wave.jl
@@ -56,7 +56,7 @@ a(x, z) = A * exp( -( (x - x₀)^2 + (z - z₀)^2 ) / 2δ^2 )
 u₀(x, y, z) = a(x, z) * U * cos(k*x + m*z)
 v₀(x, y, z) = a(x, z) * V * sin(k*x + m*z)
 w₀(x, y, z) = a(x, z) * W * cos(k*x + m*z)
-b₀(x, y, z) = a(x, z) * B * sin(k*x + m*z) + N^2 * z 
+b₀(x, y, z) = a(x, z) * B * sin(k*x + m*z) + N^2 * z
 
 # We are now ready to instantiate our model on a uniform grid.
 # We give the model a constant rotation rate with background vorticity `f`,
@@ -64,9 +64,9 @@ b₀(x, y, z) = a(x, z) * B * sin(k*x + m*z) + N^2 * z
 # and diffusivity to stabilize the model.
 
 model = Model(
-        grid = RegularCartesianGrid(N=(Nx, 1, Nx), L=(Lx, Lx, Lx)),
+        grid = RegularCartesianGrid(size=(Nx, 1, Nx), length=(Lx, Lx, Lx)),
      closure = ConstantIsotropicDiffusivity(ν=1e-6, κ=1e-6),
-    coriolis = FPlane(f=f), 
+    coriolis = FPlane(f=f),
      tracers = :b,
     buoyancy = BuoyancyTracer()
 )
diff --git a/examples/netcdf_ouput_example.jl b/examples/netcdf_ouput_example.jl
index a9c48259db..af6788333b 100644
--- a/examples/netcdf_ouput_example.jl
+++ b/examples/netcdf_ouput_example.jl
@@ -8,7 +8,7 @@ Nx, Ny, Nz = 16, 16, 16       # No. of grid points in x, y, and z, respectively.
 Lx, Ly, Lz = 100, 100, 100    # Length of the domain in x, y, and z, respectively (m).
 tf = 5000                     # Length of the simulation (s)
 
-model = Model(grid=RegularCartesianGrid(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz)),
+model = Model(grid=RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
               closure=ConstantIsotropicDiffusivity())
 
 # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
diff --git a/examples/ocean_convection_with_plankton.jl b/examples/ocean_convection_with_plankton.jl
index da34d4d150..35e45dc961 100644
--- a/examples/ocean_convection_with_plankton.jl
+++ b/examples/ocean_convection_with_plankton.jl
@@ -6,7 +6,7 @@
 #   * to set boundary conditions;
 #   * to defined and insert a user-defined forcing function into a simulation.
 #   * to use the `TimeStepWizard` to manage and adapt the simulation time-step.
-# 
+#
 # To begin, we load Oceananigans, a plotting package, and a few miscellaneous useful packages.
 
 using Oceananigans, PyPlot, Random, Printf
@@ -14,12 +14,12 @@ using Oceananigans, PyPlot, Random, Printf
 # ## Parameters
 #
 # We choose a modest two-dimensional resolution of 128² in a 64² m² domain ,
-# implying a resolution of 0.5 m. Our fluid is initially stratified with 
+# implying a resolution of 0.5 m. Our fluid is initially stratified with
 # a squared buoyancy frequency
 #
 # $$ N\^2 = 10^{-5} \, \mathrm{s^{-2}} $$
 #
-# and a surface buoyancy flux 
+# and a surface buoyancy flux
 #
 # $$ Q_b2 = 10^{-8} \, \mathrm{m^3 \, s^{-2}} $$
 #
@@ -39,7 +39,7 @@ end_time = 1day
 # Create boundary conditions. Note that temperature is buoyancy in our problem.
 #
 
-buoyancy_bcs = HorizontallyPeriodicBCs(   top = BoundaryCondition(Flux, Qb), 
+buoyancy_bcs = HorizontallyPeriodicBCs(   top = BoundaryCondition(Flux, Qb),
                                        bottom = BoundaryCondition(Gradient, N²))
 
 # ## Define a forcing function
@@ -51,10 +51,10 @@ buoyancy_bcs = HorizontallyPeriodicBCs(   top = BoundaryCondition(Flux, Qb),
 growth_and_decay = SimpleForcing((x, y, z, t) -> exp(z/16) - 1)
 
 ## Instantiate the model
-model = Model(      
-                   grid = RegularCartesianGrid(N = (Nz, 1, Nz), L = (Lz, Lz, Lz)),
+model = Model(
+                   grid = RegularCartesianGrid(size = (Nz, 1, Nz), length = (Lz, Lz, Lz)),
                 closure = ConstantIsotropicDiffusivity(ν=1e-4, κ=1e-4),
-               coriolis = FPlane(f=1e-4), 
+               coriolis = FPlane(f=1e-4),
                 tracers = (:b, :plankton),
                buoyancy = BuoyancyTracer(),
                 forcing = ModelForcing(plankton=growth_and_decay),
diff --git a/examples/ocean_wind_mixing_and_convection.jl b/examples/ocean_wind_mixing_and_convection.jl
index 45bd92a6d0..d7eddb72cd 100644
--- a/examples/ocean_wind_mixing_and_convection.jl
+++ b/examples/ocean_wind_mixing_and_convection.jl
@@ -18,30 +18,30 @@ using Oceananigans, PyPlot, Random, Printf
 # Here we use an isotropic, cubic grid with `Nz` grid points and grid spacing
 # `Δz = 1` meter. We specify fluxes of heat, momentum, and salinity via
 #
-#   1. A temperature flux `Qᵀ` at the top of the domain, which is related to heat flux 
-#       by `Qᵀ = Qʰ / (ρ₀ * cᴾ)`, where `Qʰ` is the heat flux, `ρ₀` is a reference density, 
-#       and `cᴾ` is the heat capacity of seawater. With a reference density 
-#       `ρ₀ = 1026 kg m⁻³`and heat capacity `cᴾ = 3991`, our chosen temperature flux of 
-#       `Qᵀ = 5 × 10⁻⁵ K m⁻¹ s⁻¹` corresponds to a heat flux of `Qʰ = 204.7 W m⁻²`, a 
+#   1. A temperature flux `Qᵀ` at the top of the domain, which is related to heat flux
+#       by `Qᵀ = Qʰ / (ρ₀ * cᴾ)`, where `Qʰ` is the heat flux, `ρ₀` is a reference density,
+#       and `cᴾ` is the heat capacity of seawater. With a reference density
+#       `ρ₀ = 1026 kg m⁻³`and heat capacity `cᴾ = 3991`, our chosen temperature flux of
+#       `Qᵀ = 5 × 10⁻⁵ K m⁻¹ s⁻¹` corresponds to a heat flux of `Qʰ = 204.7 W m⁻²`, a
 #       relatively powerful cooling rate.
 #
 #   2. A velocity flux `Qᵘ` at the top of the domain, which is related
 #       to the `x` momentum flux `τˣ` via `τˣ = ρ₀ * Qᵘ`, where `ρ₀` is a reference density.
 #       Our chosen value of `Qᵘ = -2 × 10⁻⁵ m² s⁻²` roughly corresponds to atmospheric winds
-#       of `uᵃ = 2.9 m s⁻¹` in the positive `x`-direction, using the parameterization 
+#       of `uᵃ = 2.9 m s⁻¹` in the positive `x`-direction, using the parameterization
 #       `τ = 0.0025 * |uᵃ| * uᵃ`.
 #
-#   3. An evaporation rate `evaporation = 10⁻⁷ m s⁻¹`, or approximately 0.1 millimeter per 
+#   3. An evaporation rate `evaporation = 10⁻⁷ m s⁻¹`, or approximately 0.1 millimeter per
 #       hour.
 #
 # Finally, we use an initial temperature gradient of `∂T/∂z = 0.005 K m⁻¹`,
 # which implies an iniital buoyancy frequency `N² = α * g * ∂T/∂z = 9.8 × 10⁻⁶ s⁻²`
-# with a thermal expansion coefficient `α = 2 × 10⁻⁴ K⁻¹` and gravitational acceleration 
-# `g = 9.81 s⁻²`. Note that, by default, the `SeawaterBuoyancy` model uses a gravitational 
+# with a thermal expansion coefficient `α = 2 × 10⁻⁴ K⁻¹` and gravitational acceleration
+# `g = 9.81 s⁻²`. Note that, by default, the `SeawaterBuoyancy` model uses a gravitational
 # acceleration `gᴱᵃʳᵗʰ = 9.80665 s⁻²`.
 
          Nz = 48       # Number of grid points in x, y, z
-         Δz = 1.0      # Grid spacing in x, y, z (meters)         
+         Δz = 1.0      # Grid spacing in x, y, z (meters)
          Qᵀ = 5e-5     # Temperature flux at surface
          Qᵘ = -2e-5    # Velocity flux at surface
        ∂T∂z = 0.005    # Initial vertical temperature gradient
@@ -70,9 +70,9 @@ S_bcs = HorizontallyPeriodicBCs(top = BoundaryCondition(Flux, Qˢ))
 
 # ## Model instantiation
 #
-# We instantiate a horizontally-periodic `Model` on the CPU with on a `RegularCartesianGrid`, 
-# using a `FPlane` model for rotation (constant rotation rate), a linear equation 
-# of state for temperature and salinity, the Anisotropic Minimum Dissipation closure 
+# We instantiate a horizontally-periodic `Model` on the CPU with on a `RegularCartesianGrid`,
+# using a `FPlane` model for rotation (constant rotation rate), a linear equation
+# of state for temperature and salinity, the Anisotropic Minimum Dissipation closure
 # to model the effects of unresolved turbulence, and the previously defined boundary
 # conditions for `u`, `T`, and `S`. We also pass the evaporation rate to the container
 # model.parameters for use in the boundary condition function that calculates the salinity
@@ -80,7 +80,7 @@ S_bcs = HorizontallyPeriodicBCs(top = BoundaryCondition(Flux, Qˢ))
 
 model = Model(
          architecture = CPU(),
-                 grid = RegularCartesianGrid(N=(Nz, Nz, Nz), L=(Δz*Nz, Δz*Nz, Δz*Nz)),
+                 grid = RegularCartesianGrid(size=(Nz, Nz, Nz), length=(Δz*Nz, Δz*Nz, Δz*Nz)),
              coriolis = FPlane(f=f),
              buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=α, β=β)),
               closure = AnisotropicMinimumDissipation(),
@@ -95,14 +95,14 @@ model = Model(
 # ```
 #
 # To change the `architecture` to `GPU`, replace the `architecture` keyword argument with
-# 
+#
 # ```julia
 # architecture = GPU()
 # ```
 #
 # ## Initial conditions
 #
-# Out initial condition for temperature consists of a linear stratification superposed with 
+# Out initial condition for temperature consists of a linear stratification superposed with
 # random noise damped at the walls, while our initial condition for velocity consists
 # only of random noise.
 
@@ -120,24 +120,24 @@ set!(model, u=u₀, w=u₀, T=T₀, S=35)
 # ## Set up output
 #
 # We set up an output writer that saves all velocity fields, tracer fields, and the subgrid
-# turbulent diffusivity associated with `model.closure`. The `prefix` keyword argument 
-# to `JLD2OutputWriter` indicates that output will be saved in 
+# turbulent diffusivity associated with `model.closure`. The `prefix` keyword argument
+# to `JLD2OutputWriter` indicates that output will be saved in
 # `ocean_wind_mixing_and_convection.jld2`.
 
 ## Create a NamedTuple containing all the fields to be outputted.
 fields_to_output = merge(model.velocities, model.tracers, (νₑ=model.diffusivities.νₑ,))
 
 ## Instantiate a JLD2OutputWriter to write fields.
-field_writer = JLD2OutputWriter(model, FieldOutputs(fields_to_output); interval=hour/4, 
+field_writer = JLD2OutputWriter(model, FieldOutputs(fields_to_output); interval=hour/4,
                                 prefix="ocean_wind_mixing_and_convection", force=true)
-                                
+
 ## Add the output writer to the models `output_writers`.
 model.output_writers[:fields] = field_writer
 
 # ## Running the simulation
 #
 # To run the simulation, we instantiate a `TimeStepWizard` to ensure stable time-stepping
-# with a Courant-Freidrichs-Lewy (CFL) number of 0.2. 
+# with a Courant-Freidrichs-Lewy (CFL) number of 0.2.
 
 wizard = TimeStepWizard(cfl=0.2, Δt=1.0, max_change=1.1, max_Δt=5.0)
 
@@ -153,11 +153,11 @@ fig, axs = subplots(ncols=3, figsize=(12, 5))
 """
     makeplot!(axs, model)
 
-Make a triptych of x-z slices of vertical velocity, temperature, and salinity 
+Make a triptych of x-z slices of vertical velocity, temperature, and salinity
 associated with `model` in `axs`.
 """
 function makeplot!(axs, model)
-    jhalf = floor(Int, model.grid.Nz/2) 
+    jhalf = floor(Int, model.grid.Nz/2)
 
     ## Coordinate arrays for plotting
     xC = repeat(model.grid.xC, 1, model.grid.Nz)
@@ -169,7 +169,7 @@ function makeplot!(axs, model)
     pcolormesh(xC, zF, data(model.velocities.w)[:, jhalf, :])
     xlabel("\$ x \$ (m)"); ylabel("\$ z \$ (m)")
 
-    sca(axs[2]); cla() 
+    sca(axs[2]); cla()
     title("Temperature")
     pcolormesh(xC, zC, data(model.tracers.T)[:, jhalf, :])
     xlabel("\$ x \$ (m)")
@@ -196,9 +196,9 @@ while model.clock.time < end_time
     walltime = @elapsed time_step!(model, 10, wizard.Δt)
 
     ## Print a progress message
-    @printf("i: %04d, t: %s, Δt: %s, wmax = %.1e ms⁻¹, wall time: %s\n", 
-            model.clock.iteration, prettytime(model.clock.time), prettytime(wizard.Δt), 
-            wmax(model), prettytime(walltime))           
+    @printf("i: %04d, t: %s, Δt: %s, wmax = %.1e ms⁻¹, wall time: %s\n",
+            model.clock.iteration, prettytime(model.clock.time), prettytime(wizard.Δt),
+            wmax(model), prettytime(walltime))
 
     model.architecture == CPU() && makeplot!(axs, model)
 end
diff --git a/examples/simple_diffusion.jl b/examples/simple_diffusion.jl
index 7f425c7137..2d4ef33fb5 100644
--- a/examples/simple_diffusion.jl
+++ b/examples/simple_diffusion.jl
@@ -26,7 +26,7 @@ using PyPlot, Printf
 # `Model` constructor:
 
 model = Model(
-    grid = RegularCartesianGrid(N = (1, 1, 128), L = (1, 1, 1)),
+    grid = RegularCartesianGrid(size = (1, 1, 128), length = (1, 1, 1)),
     closure = ConstantIsotropicDiffusivity(κ = 1.0)
 )
 
diff --git a/examples/two_dimensional_turbulence.jl b/examples/two_dimensional_turbulence.jl
index a91260d46c..02352b1c96 100644
--- a/examples/two_dimensional_turbulence.jl
+++ b/examples/two_dimensional_turbulence.jl
@@ -22,8 +22,8 @@ using Oceananigans.TurbulenceClosures: ∂x_faa, ∂y_afa
 # We instantiate the model with a simple isotropic diffusivity
 
 model = Model(
-        grid = RegularCartesianGrid(N=(128, 128, 1), L=(2π, 2π, 2π)),
-    buoyancy = nothing, 
+        grid = RegularCartesianGrid(size=(128, 128, 1), length=(2π, 2π, 2π)),
+    buoyancy = nothing,
      tracers = nothing,
      closure = ConstantIsotropicDiffusivity(ν=1e-3, κ=1e-3)
 )
diff --git a/src/fields.jl b/src/fields.jl
index 26ba8115c5..ee3c4842c1 100644
--- a/src/fields.jl
+++ b/src/fields.jl
@@ -254,7 +254,7 @@ a function with arguments `(x, y, z)`, or any data type for which a
 Example
 =======
 ```julia
-model = Model(grid=RegularCartesianGrid(N=(32, 32, 32), L=(1, 1, 1))
+model = Model(grid=RegularCartesianGrid(size=(32, 32, 32), length=(1, 1, 1))
 
 # Set u to a parabolic function of z, v to random numbers damped
 # at top and bottom, and T to some silly array of half zeros,
diff --git a/src/grids.jl b/src/grids.jl
index 291762689f..fa03315752 100644
--- a/src/grids.jl
+++ b/src/grids.jl
@@ -1,10 +1,10 @@
 """
-    RegularCartesianGrid{T<:AbstractFloat, R<:AbstractRange} <: AbstractGrid{T}
+    RegularCartesianGrid{FT<:AbstractFloat, R<:AbstractRange} <: AbstractGrid{FT}
 
 A Cartesian grid with with constant grid spacings `Δx`, `Δy`, and `Δz` between cell centers
 and cell faces.
 """
-struct RegularCartesianGrid{T<:AbstractFloat, R<:AbstractRange} <: AbstractGrid{T}
+struct RegularCartesianGrid{FT<:AbstractFloat, R<:AbstractRange} <: AbstractGrid{FT}
     # Number of grid points in (x,y,z).
     Nx::Int
     Ny::Int
@@ -18,19 +18,13 @@ struct RegularCartesianGrid{T<:AbstractFloat, R<:AbstractRange} <: AbstractGrid{
     Ty::Int
     Tz::Int
     # Domain size [m].
-    Lx::T
-    Ly::T
-    Lz::T
+    Lx::FT
+    Ly::FT
+    Lz::FT
     # Grid spacing [m].
-    Δx::T
-    Δy::T
-    Δz::T
-    # Cell face areas [m²].
-    Ax::T
-    Ay::T
-    Az::T
-    # Volume of a cell [m³].
-    V::T
+    Δx::FT
+    Δy::FT
+    Δz::FT
     # Range of coordinates at the centers of the cells.
     xC::R
     yC::R
@@ -43,13 +37,19 @@ struct RegularCartesianGrid{T<:AbstractFloat, R<:AbstractRange} <: AbstractGrid{
 end
 
 """
-    RegularCartesianGrid([T=Float64]; N, L)
+    RegularCartesianGrid([FT=Float64]; size, length, x, y, z)
 
-Creates a `RegularCartesianGrid` with `N = (Nx, Ny, Nz)` grid points and domain size
-`L = (Lx, Ly, Lz)`, where constants are stored using floating point values of type `T`.
+Creates a `RegularCartesianGrid` with `size = (Nx, Ny, Nz)` grid points.
 
-Additional properties
-=====================
+The physical length of the domain can be specified via `x`, `y`, and `z` keyword arguments
+indicating the left and right endpoints of each dimensions, e.g. `x=(-π, π)` or via
+the `length` argument, e.g. `length=(Lx, Ly, Lz)` which specifies the length of each dimension
+in which case 0 ≤ x ≤ Lx, 0 ≤ y ≤ Ly, and -Lz ≤ z ≤ 0.
+
+Constants are stored using floating point values of type `FT`.
+
+Grid properties
+===============
 - `(xC, yC, zC)::AbstractRange`: (x, y, z) coordinates of cell centers
 - `(xF, yF, zF)::AbstractRange`: (x, y, z) coordinates of cell faces
 - `(Hx, Hy, Hz)::Int`: Halo size in the (x, y, z)-direction
@@ -57,35 +57,70 @@ Additional properties
 
 Examples
 ========
-```
-julia> grid = RegularCartesianGrid(N=(32, 32, 32), L=(1, 1, 1))
+```julia
+julia> grid = RegularCartesianGrid(size=(32, 32, 32), length=(1, 2, 3))
 RegularCartesianGrid{Float64}
+domain: x ∈ [0.0, 1.0], y ∈ [0.0, 2.0], z ∈ [0.0, -3.0]
   resolution (Nx, Ny, Nz) = (32, 32, 32)
    halo size (Hx, Hy, Hz) = (1, 1, 1)
-      domain (Lx, Ly, Lz) = (1.0, 1.0, 1.0)
-grid spacing (Δx, Δy, Δz) = (0.03125, 0.03125, 0.03125)
+grid spacing (Δx, Δy, Δz) = (0.03125, 0.0625, 0.09375)
 ```
-```
-julia> grid = RegularCartesianGrid(Float32; N=(32, 32, 16), L=(8, 8, 2))
+```julia
+julia> grid = RegularCartesianGrid(Float32; size=(32, 32, 16), x=(0, 8), y=(-10, 10), z=(-π, π))
 RegularCartesianGrid{Float32}
+domain: x ∈ [0.0, 8.0], y ∈ [-10.0, 10.0], z ∈ [3.141592653589793, -3.141592653589793]
   resolution (Nx, Ny, Nz) = (32, 32, 16)
    halo size (Hx, Hy, Hz) = (1, 1, 1)
-      domain (Lx, Ly, Lz) = (8.0f0, 8.0f0, 2.0f0)
-grid spacing (Δx, Δy, Δz) = (0.25f0, 0.25f0, 0.125f0)
+grid spacing (Δx, Δy, Δz) = (0.25f0, 0.625f0, 0.3926991f0)
 ```
 """
-function RegularCartesianGrid(T, N, L)
-    length(N) == 3 || throw(ArgumentError("N=$N must be a tuple of length 3."))
-    length(L) == 3 || throw(ArgumentError("L=$L must be a tuple of length 3."))
-
-    all(isa.(N, Integer)) || throw(ArgumentError("N=$N should contain integers."))
-    all(isa.(L, Number))  || throw(ArgumentError("L=$L should contain numbers."))
-
-    all(N .>= 1) || throw(ArgumentError("N=$N must be nonzero and positive!"))
-    all(L .> 0)  || throw(ArgumentError("L=$L must be nonzero and positive!"))
-
-    Nx, Ny, Nz = N
-    Lx, Ly, Lz = L
+function RegularCartesianGrid(FT=Float64; size, length=nothing, x=nothing, y=nothing, z=nothing)
+    # Hack that allows us to use `size` and `length` as keyword arguments but then also
+    # use the `size` and `length` functions.
+    sz, len = size, length
+    length = Base.length
+
+    if isnothing(len) && (isnothing(x) || isnothing(y) || isnothing(z))
+        throw(ArgumentError("Must supply length or x, y, z keyword arguments."))
+    end
+
+    if !isnothing(len) && !isnothing(x) && !isnothing(y) && !isnothing(z)
+        throw(ArgumentError("Cannot specify both length and x, y, z keyword arguments."))
+    end
+
+    length(sz) == 3        || throw(ArgumentError("length($sz) must be 3."))
+    all(isa.(sz, Integer)) || throw(ArgumentError("size=$sz should contain integers."))
+    all(sz .>= 1)          || throw(ArgumentError("size=$sz must be nonzero and positive!"))
+
+    if !isnothing(len)
+        length(len) == 3       || throw(ArgumentError("length($len) must be 3."))
+        all(isa.(len, Number)) || throw(ArgumentError("length=$len should contain numbers."))
+        all(len .>= 0)         || throw(ArgumentError("length=$len must be nonzero and positive!"))
+
+        Lx, Ly, Lz = len
+        x = (0, Lx)
+        y = (0, Ly)
+        z = (-Lz, 0)
+    end
+
+    function coord2xyz(c)
+        c == 1 && return "x"
+        c == 2 && return "y"
+        c == 3 && return "z"
+    end
+
+    for (i, c) in enumerate((x, y, z))
+        name = coord2xyz(i)
+        length(c) == 2       || throw(ArgumentError("$name length($c) must be 2."))
+        all(isa.(c, Number)) || throw(ArgumentError("$name=$c should contain numbers."))
+        c[2] >= c[1]         || throw(ArgumentError("$name=$c should be an increasing interval."))
+    end
+
+    x₁, x₂ = x[1], x[2]
+    y₁, y₂ = y[1], y[2]
+    z₁, z₂ = z[1], z[2]
+    Nx, Ny, Nz = sz
+    Lx, Ly, Lz = x₂-x₁, y₂-y₁, z₂-z₁
 
     # Right now we only support periodic horizontal boundary conditions and
     # usually use second-order advection schemes so halos of size Hx, Hy = 1 are
@@ -106,29 +141,25 @@ function RegularCartesianGrid(T, N, L)
 
     V = Δx*Δy*Δz
 
-    xC = range(Δx/2, Lx-Δx/2; length=Nx)
-    yC = range(Δy/2, Ly-Δy/2; length=Ny)
-    zC = range(-Lz+Δz/2, -Δz/2; length=Nz)
+    xC = range(x₁ + Δx/2, x₂ - Δx/2; length=Nx)
+    yC = range(y₁ + Δy/2, y₂ - Δy/2; length=Ny)
+    zC = range(z₁ + Δz/2, z₂ - Δz/2; length=Nz)
 
-    xF = range(0, Lx; length=Nx+1)
-    yF = range(0, Ly; length=Ny+1)
-    zF = range(-Lz, 0; length=Nz+1)
+    xF = range(x₁, x₂; length=Nx+1)
+    yF = range(y₁, y₂; length=Ny+1)
+    zF = range(z₁, z₂; length=Nz+1)
 
-    RegularCartesianGrid{T, typeof(xC)}(Nx, Ny, Nz, Hx, Hy, Hz, Tx, Ty, Tz,
-                                        Lx, Ly, Lz, Δx, Δy, Δz, Ax, Ay, Az, V,
-                                        xC, yC, zC, xF, yF, zF)
+    RegularCartesianGrid{FT, typeof(xC)}(Nx, Ny, Nz, Hx, Hy, Hz, Tx, Ty, Tz,
+                                         Lx, Ly, Lz, Δx, Δy, Δz, xC, yC, zC, xF, yF, zF)
 end
 
-RegularCartesianGrid(N, L) = RegularCartesianGrid(Float64, N, L)
-
-RegularCartesianGrid(T=Float64; N, L) = RegularCartesianGrid(T, N, L)
-
-size(g::RegularCartesianGrid) = (g.Nx, g.Ny, g.Nz)
-eltype(g::RegularCartesianGrid{T}) where T = T
+size(grid::RegularCartesianGrid)   = (grid.Nx, grid.Ny, grid.Nz)
+length(grid::RegularCartesianGrid) = (grid.Lx, grid.Ly, grid.Lz)
+eltype(grid::RegularCartesianGrid{FT}) where FT = FT
 
 show(io::IO, g::RegularCartesianGrid) =
     print(io, "RegularCartesianGrid{$(eltype(g))}\n",
+              "domain: x ∈ [$(g.xF[1]), $(g.xF[end])], y ∈ [$(g.yF[1]), $(g.yF[end])], z ∈ [$(g.zF[end]), $(g.zF[1])]", '\n',
               "  resolution (Nx, Ny, Nz) = ", (g.Nx, g.Ny, g.Nz), '\n',
               "   halo size (Hx, Hy, Hz) = ", (g.Hx, g.Hy, g.Hz), '\n',
-              "      domain (Lx, Ly, Lz) = ", (g.Lx, g.Ly, g.Lz), '\n',
               "grid spacing (Δx, Δy, Δz) = ", (g.Δx, g.Δy, g.Δz))
diff --git a/src/models.jl b/src/models.jl
index 0026006f9d..949e9e8e4c 100644
--- a/src/models.jl
+++ b/src/models.jl
@@ -98,7 +98,7 @@ ChannelModel(; boundary_conditions=ChannelSolutionBCs(), kwargs...) =
 function BasicChannelModel(; N, L, ν=ν₀, κ=κ₀, float_type=Float64,
                            boundary_conditions=ChannelSolutionBCs(), kwargs...)
 
-    grid = RegularCartesianGrid(float_type, N, L)
+    grid = RegularCartesianGrid(float_type; size=N, length=L)
     closure = ConstantIsotropicDiffusivity(float_type, ν=ν, κ=κ)
 
     return Model(; float_type=float_type, grid=grid, closure=closure,
@@ -114,7 +114,7 @@ precision `float_type`, and constant isotropic viscosity and diffusivity `ν`, a
 Additional `kwargs` are passed to the regular `Model` constructor.
 """
 function BasicModel(; N, L, ν=ν₀, κ=κ₀, float_type=Float64, kwargs...)
-    grid = RegularCartesianGrid(float_type, N, L)
+    grid = RegularCartesianGrid(float_type; size=N, length=L)
     closure = ConstantIsotropicDiffusivity(float_type, ν=ν, κ=κ)
     return Model(; float_type=float_type, grid=grid, closure=closure, kwargs...)
 end
@@ -132,7 +132,7 @@ precision `float_type`, and the four non-dimensional numbers:
 for characteristic velocity scale `U`, length-scale `λ`, viscosity `ν`,
 tracer diffusivity `κ`, and Coriolis parameter `f`. Buoyancy is scaled
 with `λ U²`, so that the Richardson number is `Ri=B`, where `B` is a
-non-dimensional buoyancy scale set by the user via initial conditions or 
+non-dimensional buoyancy scale set by the user via initial conditions or
 forcing.
 
 Note that `N`, `L`, and `Re` are required.
@@ -141,7 +141,7 @@ Additional `kwargs` are passed to the regular `Model` constructor.
 """
 function NonDimensionalModel(; N, L, Re, Pr=0.7, Ro=Inf, float_type=Float64, kwargs...)
 
-         grid = RegularCartesianGrid(float_type, N, L)
+         grid = RegularCartesianGrid(float_type; size=N, length=L)
       closure = ConstantIsotropicDiffusivity(float_type, ν=1/Re, κ=1/(Pr*Re))
      coriolis = VerticalRotationAxis(float_type, f=1/Ro)
      buoyancy = BuoyancyTracer()
diff --git a/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl b/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl
index 27209fc739..1a438fec29 100644
--- a/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl
+++ b/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl
@@ -63,7 +63,7 @@ function run_ocean_large_eddy_simulation_regression_test(arch, closure)
     # Model instantiation
     model = Model(
              architecture = arch,
-                     grid = RegularCartesianGrid(N=(16, 16, 16), L=(16, 16, 16)),
+                     grid = RegularCartesianGrid(size=(16, 16, 16), length=(16, 16, 16)),
                  coriolis = FPlane(f=1e-4),
                  buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=2e-4, β=8e-4)),
                   closure = closure,
diff --git a/test/regression_tests/rayleigh_benard_regression_test.jl b/test/regression_tests/rayleigh_benard_regression_test.jl
index d913c605b3..77577a7408 100644
--- a/test/regression_tests/rayleigh_benard_regression_test.jl
+++ b/test/regression_tests/rayleigh_benard_regression_test.jl
@@ -32,7 +32,7 @@ function run_rayleigh_benard_regression_test(arch)
 
     model = Model(
                architecture = arch,
-                       grid = RegularCartesianGrid(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz)),
+                       grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
                     closure = ConstantIsotropicDiffusivity(ν=ν, κ=κ),
                     tracers = (:b, :c),
                    buoyancy = BuoyancyTracer(),
diff --git a/test/test_buoyancy.jl b/test/test_buoyancy.jl
index e9bcb112aa..e33d5401fd 100644
--- a/test/test_buoyancy.jl
+++ b/test/test_buoyancy.jl
@@ -9,34 +9,34 @@ function instantiate_roquet_equations_of_state(FT, flavor; coeffs=nothing)
     return typeof(eos.polynomial_coeffs.R₁₀₀) == FT
 end
 
-function instantiate_seawater_buoyancy(FT, EquationOfState) 
+function instantiate_seawater_buoyancy(FT, EquationOfState)
     buoyancy = SeawaterBuoyancy(FT, equation_of_state=EquationOfState(FT))
     return typeof(buoyancy.gravitational_acceleration) == FT
 end
 
 function density_perturbation_works(arch, FT, eos)
-    grid = RegularCartesianGrid(FT, N=(3, 3, 3), L=(1, 1, 1))
+    grid = RegularCartesianGrid(FT; size=(3, 3, 3), length=(1, 1, 1))
     C = datatuple(TracerFields(arch, grid, (:T, :S)))
     density_anomaly = ρ′(2, 2, 2, grid, eos, C)
     return true
 end
 
 function buoyancy_frequency_squared_works(arch, FT, buoyancy)
-    grid = RegularCartesianGrid(N=(3, 3, 3), L=(1, 1, 1))
+    grid = RegularCartesianGrid(FT; size=(3, 3, 3), length=(1, 1, 1))
     C = datatuple(TracerFields(arch, grid, required_tracers(buoyancy)))
     N² = buoyancy_frequency_squared(2, 2, 2, grid, buoyancy, C)
     return true
 end
 
 function thermal_expansion_works(arch, FT, eos)
-    grid = RegularCartesianGrid(FT, N=(3, 3, 3), L=(1, 1, 1))
+    grid = RegularCartesianGrid(FT; size=(3, 3, 3), length=(1, 1, 1))
     C = datatuple(TracerFields(arch, grid, (:T, :S)))
     α = thermal_expansion(2, 2, 2, grid, eos, C)
     return true
 end
 
 function haline_contraction_works(arch, FT, eos)
-    grid = RegularCartesianGrid(N=(3, 3, 3), L=(1, 1, 1))
+    grid = RegularCartesianGrid(FT; size=(3, 3, 3), length=(1, 1, 1))
     C = datatuple(TracerFields(arch, grid, (:T, :S)))
     β = haline_contraction(2, 2, 2, grid, eos, C)
     return true
@@ -44,7 +44,7 @@ end
 
 @testset "Buoyancy" begin
     println("Testing buoyancy...")
-    
+
     @testset "Equations of State" begin
         for FT in float_types
             @test instantiate_linear_equation_of_state(FT, 0.1, 0.3)
@@ -75,7 +75,7 @@ end
             end
 
             for arch in archs
-                for EOS in EquationsOfState 
+                for EOS in EquationsOfState
                     @test thermal_expansion_works(arch, FT, EOS())
                     @test haline_contraction_works(arch, FT, EOS())
                 end
diff --git a/test/test_dynamics.jl b/test/test_dynamics.jl
index 97d4603c53..9e63be435b 100644
--- a/test/test_dynamics.jl
+++ b/test/test_dynamics.jl
@@ -165,7 +165,7 @@ function taylor_green_vortex_test(arch; FT=Float64, N=64, Nt=10)
     @inline v(x, y, z, t) =  sin(2π*x) * exp(-4π^2 * ν * t)
 
     model = Model(architecture = arch,
-                          grid = RegularCartesianGrid(FT; N=(Nx, Ny, Nz), L=(Lx, Ly, Lz)),
+                          grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
                        closure = ConstantIsotropicDiffusivity(FT; ν=1, κ=0),  # Turn off diffusivity.
                        tracers = nothing,
                       buoyancy = nothing) # turn off buoyancy
diff --git a/test/test_fields.jl b/test/test_fields.jl
index 09eb67f117..b3d1aa8ea0 100644
--- a/test/test_fields.jl
+++ b/test/test_fields.jl
@@ -46,7 +46,7 @@ end
     @testset "Field initialization" begin
         println("  Testing field initialization...")
         for arch in archs, FT in float_types
-            grid = RegularCartesianGrid(FT, N, L)
+            grid = RegularCartesianGrid(FT; size=N, length=L)
 
             for fieldtype in fieldtypes
                 @test correct_field_size(arch, grid, fieldtype)
@@ -65,7 +65,7 @@ end
         println("  Testing field setting...")
 
         for arch in archs, FT in float_types
-            grid = RegularCartesianGrid(FT, N, L)
+            grid = RegularCartesianGrid(FT; size=N, length=L)
 
             for fieldtype in fieldtypes, val in vals
                 @test correct_field_value_was_set(arch, grid, fieldtype, val)
diff --git a/test/test_grids.jl b/test/test_grids.jl
index e3f135ae8b..95f8788759 100644
--- a/test/test_grids.jl
+++ b/test/test_grids.jl
@@ -1,28 +1,18 @@
 function correct_grid_size(FT)
-    grid = RegularCartesianGrid(FT, (4, 6, 8), (2π, 4π, 9π))
+    grid = RegularCartesianGrid(FT; size=(4, 6, 8), length=(2π, 4π, 9π))
 
     # Checking ≈ as the grid could be storing Float32 values.
     return (grid.Nx ≈ 4  && grid.Ny ≈ 6  && grid.Nz ≈ 8 &&
             grid.Lx ≈ 2π && grid.Ly ≈ 4π && grid.Lz ≈ 9π)
 end
 
-function correct_cell_volume(FT)
-    Nx, Ny, Nz = 19, 13, 7
-    Δx, Δy, Δz = 0.1, 0.2, 0.3
-    Lx, Ly, Lz = Nx*Δx, Ny*Δy, Nz*Δz
-    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (Lx, Ly, Lz))
-
-    # Checking ≈ as the grid could be storing Float32 values.
-    return grid.V ≈ Δx*Δy*Δz
-end
-
 function faces_start_at_zero(FT)
-    grid = RegularCartesianGrid(FT, (10, 10, 10), (2π, 2π, 2π))
+    grid = RegularCartesianGrid(FT; size=(10, 10, 10), length=(2π, 2π, 2π))
     return grid.xF[1] == 0 && grid.yF[1] == 0 && grid.zF[end] == 0
 end
 
 function end_faces_match_grid_length(FT)
-    grid = RegularCartesianGrid(FT, (12, 13, 14), (π, π^2, π^3))
+    grid = RegularCartesianGrid(FT; size=(12, 13, 14), length=(π, π^2, π^3))
     return (grid.xF[end] - grid.xF[1] ≈ π   &&
             grid.yF[end] - grid.yF[1] ≈ π^2 &&
             grid.zF[end] - grid.zF[1] ≈ π^3)
@@ -30,7 +20,7 @@ end
 
 function ranges_have_correct_length(FT)
     Nx, Ny, Nz = 8, 9, 10
-    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (1, 1, 1))
+    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(1, 1, 1))
     return (length(grid.xC) == Nx && length(grid.xF) == Nx+1 &&
             length(grid.yC) == Ny && length(grid.yF) == Ny+1 &&
             length(grid.zC) == Nz && length(grid.zF) == Nz+1)
@@ -39,7 +29,7 @@ end
 # See: https://github.com/climate-machine/Oceananigans.jl/issues/480
 function no_roundoff_error_in_ranges(FT)
     Nx, Ny, Nz = 1, 1, 64
-    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (1, 1, π/2))
+    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(1, 1, π/2))
     return length(grid.zC) == Nz
 end
 
@@ -50,7 +40,6 @@ end
         println("  Testing grid initialization...")
         for FT in float_types
             @test correct_grid_size(FT)
-            @test correct_cell_volume(FT)
             @test faces_start_at_zero(FT)
             @test end_faces_match_grid_length(FT)
             @test ranges_have_correct_length(FT)
@@ -59,25 +48,34 @@ end
     end
 
     @testset "Grid dimensions" begin
-        println("  Testing grid dimensions...")
+        println("  Testing grid constructor errors...")
         L = (100, 100, 100)
         for FT in float_types
-            @test isbitstype(typeof(RegularCartesianGrid(FT, (16, 16, 16), (1, 1, 1))))
+            @test isbitstype(typeof(RegularCartesianGrid(FT; size=(16, 16, 16), length=(1, 1, 1))))
+
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32,), length=L)
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 64), length=L)
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32, 16), length=L)
 
-            @test_throws ArgumentError RegularCartesianGrid(FT, (32,), L)
-            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 64), L)
-            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32, 16), L)
-            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32), (100,))
-            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32), (100, 100))
-            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32), (100, 100, 1, 1))
-            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32), (100, 100, -100))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32.0), length=(1, 1, 1))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(20.1, 32, 32), length=(1, 1, 1))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, nothing, 32), length=(1, 1, 1))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, "32", 32), length=(1, 1, 1))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32), length=(1, nothing, 1))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32), length=(1, "1", 1))
 
-            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32.0), (1, 1, 1))
-            @test_throws ArgumentError RegularCartesianGrid(FT, (20.1, 32, 32), (1, 1, 1))
-            @test_throws ArgumentError RegularCartesianGrid(FT, (32, nothing, 32), (1, 1, 1))
-            @test_throws ArgumentError RegularCartesianGrid(FT, (32, "32", 32), (1, 1, 1))
-            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32), (1, nothing, 1))
-            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32), (1, "1", 1))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=2)
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), y=[1, 2])
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), z=(-π, π))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=1, y=2, z=3)
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(0, 1), y=(0, 2), z=4)
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(-1//2, 1), y=(1//7, 5//7), z=("0", "1"))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(-1//2, 1), y=(1//7, 5//7), z=(1, 2, 3))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(1, 0), y=(1//7, 5//7), z=(1, 2))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(0, 1), y=(1, 5), z=(π, -π))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(0, 1), y=(1, 5), z=(π, -π))
+            @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), length=(1, 2, 3), x=(0, 1), y=(1, 5), z=(-π, π))
         end
     end
 end
diff --git a/test/test_halo_regions.jl b/test/test_halo_regions.jl
index b6fd7706cf..68a5b39865 100644
--- a/test/test_halo_regions.jl
+++ b/test/test_halo_regions.jl
@@ -2,7 +2,7 @@ function halo_regions_initalized_correctly(arch, FT, Nx, Ny, Nz)
     # Just choose something anisotropic to catch Δx/Δy type errors.
     Lx, Ly, Lz = 10, 20, 30
 
-    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (Lx, Ly, Lz))
+    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))
     field = CellField(FT, arch, grid)
 
     # Fill the interior with random numbers.
@@ -23,7 +23,7 @@ function halo_regions_correctly_filled(arch, FT, Nx, Ny, Nz)
     # Just choose something anisotropic to catch Δx/Δy type errors.
     Lx, Ly, Lz = 100, 200, 300
 
-    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (Lx, Ly, Lz))
+    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))
     field = CellField(FT, arch, grid)
     fbcs = HorizontallyPeriodicBCs()
 
diff --git a/test/test_models.jl b/test/test_models.jl
index dec37d5a20..88c4b9ae00 100644
--- a/test/test_models.jl
+++ b/test/test_models.jl
@@ -74,7 +74,7 @@ end
             N = (16, 16, 8)
             L = (2π, 3π, 5π)
 
-            grid = RegularCartesianGrid(FT, N, L)
+            grid = RegularCartesianGrid(FT; size=N, length=L)
             xF = reshape(grid.xF[1:end-1], N[1], 1, 1)
             yC = reshape(grid.yC, 1, N[2], 1)
             zC = reshape(grid.zC, 1, 1, N[3])
diff --git a/test/test_operators.jl b/test/test_operators.jl
index 8387f0d5be..93aa1c3751 100644
--- a/test/test_operators.jl
+++ b/test/test_operators.jl
@@ -7,7 +7,7 @@
         Lx, Ly, Lz = 100, 100, 100
 
         arch = CPU()
-        grid = RegularCartesianGrid((Nx, Ny, Nz), (Lx, Ly, Lz))
+        grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))
 
         Hx, Hy, Hz = grid.Hx, grid.Hy, grid.Hz
         Tx, Ty, Tz = grid.Tx, grid.Ty, grid.Tz
@@ -18,7 +18,7 @@
 
         # A yz-slice with Nx==1.
         A2yz = OffsetArray(zeros(1+2Hx, Ty, Tz), 1-Hx:1+Hx, 1-Hy:Ny+Hy, 1-Hz:Nz+Hz)
-        grid_yz = RegularCartesianGrid((1, Ny, Nz), (Lx, Ly, Lz))
+        grid_yz = RegularCartesianGrid(size=(1, Ny, Nz), length=(Lx, Ly, Lz))
 
         # Manually fill in halos for the slice.
         A2yz[0:2, 0:Ny+1, 1:Nz] .= A3[1:1, 0:Ny+1, 1:Nz]
@@ -27,7 +27,7 @@
 
         # An xz-slice with Ny==1.
         A2xz = OffsetArray(zeros(Tx, 1+2Hy, Tz), 1-Hx:Nx+Hx, 1-Hy:1+Hy, 1-Hz:Nz+Hz)
-        grid_xz = RegularCartesianGrid((Nx, 1, Nz), (Lx, Ly, Lz))
+        grid_xz = RegularCartesianGrid(size=(Nx, 1, Nz), length=(Lx, Ly, Lz))
 
         # Manually fill in halos for the slice.
         A2xz[0:Nx+1, 0:2, 1:Nz] .= A3[0:Nx+1, 1:1, 1:Nz]
diff --git a/test/test_poisson_solvers.jl b/test/test_poisson_solvers.jl
index 0427a4b393..d0865db0b0 100644
--- a/test/test_poisson_solvers.jl
+++ b/test/test_poisson_solvers.jl
@@ -1,4 +1,4 @@
-function ∇²!(grid::RegularCartesianGrid, f, ∇²f)
+function ∇²!(grid, f, ∇²f)
     @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
         @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
             @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
@@ -9,14 +9,14 @@ function ∇²!(grid::RegularCartesianGrid, f, ∇²f)
 end
 
 function fftw_planner_works(FT, Nx, Ny, Nz, planner_flag)
-    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (100, 100, 100))
+    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(100, 100, 100))
     solver = PoissonSolver(CPU(), PPN(), grid)
     true  # Just making sure the PoissonSolver does not error/crash.
 end
 
 function poisson_ppn_planned_div_free_cpu(FT, Nx, Ny, Nz, planner_flag)
     arch = CPU()
-    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (1.0, 2.5, 3.6))
+    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(1.0, 2.5, 3.6))
     solver = PoissonSolver(arch, PPN(), grid)
     fbcs = HorizontallyPeriodicBCs()
 
@@ -43,7 +43,7 @@ end
 
 function poisson_pnn_planned_div_free_cpu(FT, Nx, Ny, Nz, planner_flag)
     arch = CPU()
-    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (1.0, 2.5, 3.6))
+    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(1.0, 2.5, 3.6))
     solver = PoissonSolver(arch, PNN(), grid)
     fbcs = ChannelBCs()
 
@@ -72,7 +72,7 @@ end
 
 function poisson_ppn_planned_div_free_gpu(FT, Nx, Ny, Nz)
     arch = GPU()
-    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (1.0, 2.5, 3.6))
+    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(1.0, 2.5, 3.6))
     solver = PoissonSolver(arch, PPN(), grid)
     fbcs = HorizontallyPeriodicBCs()
 
@@ -109,7 +109,7 @@ end
 
 function poisson_pnn_planned_div_free_gpu(FT, Nx, Ny, Nz)
     arch = GPU()
-    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (1.0, 2.5, 3.6))
+    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(1.0, 2.5, 3.6))
     solver = PoissonSolver(arch, PNN(), grid)
     fbcs = ChannelBCs()
 
@@ -154,7 +154,7 @@ by giving it the source term or right hand side (RHS), which is
     -((\\pi m_z / L_z)^2 + (2\\pi m_y / L_y)^2 + (2\\pi m_x/L_x)^2) \\Psi(x, y, z)``.
 """
 function poisson_ppn_recover_sine_cosine_solution(FT, Nx, Ny, Nz, Lx, Ly, Lz, mx, my, mz)
-    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (Lx, Ly, Lz))
+    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))
     solver = PoissonSolver(CPU(), PPN(), grid)
 
     xC, yC, zC = grid.xC, grid.yC, grid.zC
diff --git a/test/test_time_stepping.jl b/test/test_time_stepping.jl
index 4da6789ab9..c66350b66e 100644
--- a/test/test_time_stepping.jl
+++ b/test/test_time_stepping.jl
@@ -29,7 +29,7 @@ function compute_w_from_continuity(arch, FT)
     Nx, Ny, Nz = 16, 16, 16
     Lx, Ly, Lz = 16, 16, 16
 
-    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (Lx, Ly, Lz))
+    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))
     bcs = HorizontallyPeriodicSolutionBCs()
 
     u = FaceFieldX(FT, arch, grid)
@@ -114,7 +114,7 @@ function tracer_conserved_in_channel(arch, FT, Nt)
     νv, κv = α*νh, α*κh
 
     model = ChannelModel(architecture = arch, float_type = FT,
-                         grid = RegularCartesianGrid(N = (Nx, Ny, Nz), L = (Lx, Ly, Lz)),
+                         grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
                          closure = ConstantAnisotropicDiffusivity(νh=νh, νv=νv, κh=κh, κv=κv))
 
     Ty = 1e-4  # Meridional temperature gradient [K/m].
diff --git a/test/test_turbulence_closures.jl b/test/test_turbulence_closures.jl
index 9aa39dee54..787a9468b4 100644
--- a/test/test_turbulence_closures.jl
+++ b/test/test_turbulence_closures.jl
@@ -15,7 +15,7 @@ function test_calculate_diffusivities(arch, closurename, FT=Float64; kwargs...)
       tracernames = (:b,)
           closure = getproperty(TurbulenceClosures, closurename)(FT; kwargs...)
           closure = with_tracers(tracernames, closure)
-             grid = RegularCartesianGrid(FT, (3, 3, 3), (3, 3, 3))
+             grid = RegularCartesianGrid(FT; size=(3, 3, 3), length=(3, 3, 3))
     diffusivities = TurbulentDiffusivities(arch, grid, tracernames, closure)
          buoyancy = BuoyancyTracer()
        velocities = Oceananigans.VelocityFields(arch, grid)
@@ -36,7 +36,7 @@ end
 function test_constant_isotropic_diffusivity_fluxdiv(FT=Float64; ν=FT(0.3), κ=FT(0.7))
           arch = CPU()
        closure = ConstantIsotropicDiffusivity(FT, κ=(T=κ, S=κ), ν=ν)
-          grid = RegularCartesianGrid(FT, (3, 1, 4), (3, 1, 4))
+          grid = RegularCartesianGrid(FT; size=(3, 1, 4), length=(3, 1, 4))
            bcs = SolutionBoundaryConditions((:T, :S), HorizontallyPeriodicSolutionBCs())
     velocities = Oceananigans.VelocityFields(arch, grid)
        tracers = Oceananigans.TracerFields(arch, grid, (:T, :S))
@@ -63,7 +63,7 @@ end
 function test_anisotropic_diffusivity_fluxdiv(FT=Float64; νh=FT(0.3), κh=FT(0.7), νv=FT(0.1), κv=FT(0.5))
           arch = CPU()
        closure = ConstantAnisotropicDiffusivity(FT, νh=νh, νv=νv, κh=(T=κh, S=κh), κv=(T=κv, S=κv))
-          grid = RegularCartesianGrid(FT, (3, 1, 4), (3, 1, 4))
+          grid = RegularCartesianGrid(FT; size=(3, 1, 4), length=(3, 1, 4))
            bcs = SolutionBoundaryConditions((:T, :S), HorizontallyPeriodicSolutionBCs())
       buoyancy = SeawaterBuoyancy(FT, gravitational_acceleration=1, equation_of_state=LinearEquationOfState(FT))
     velocities = Oceananigans.VelocityFields(arch, grid)
@@ -97,7 +97,7 @@ function test_anisotropic_diffusivity_fluxdiv(FT=Float64; νh=FT(0.3), κh=FT(0.
 end
 
 function test_function_interpolation(T=Float64)
-    grid = RegularCartesianGrid(T, (3, 3, 3), (3, 3, 3))
+    grid = RegularCartesianGrid(T; size=(3, 3, 3), length=(3, 3, 3))
     ϕ = rand(T, 3, 3, 3)
     ϕ² = ϕ.^2
 
@@ -125,7 +125,7 @@ function test_function_interpolation(T=Float64)
 end
 
 function test_function_differentiation(T=Float64)
-    grid = RegularCartesianGrid(T, (3, 3, 3), (3, 3, 3))
+    grid = RegularCartesianGrid(T; size=(3, 3, 3), length=(3, 3, 3))
     ϕ = rand(T, 3, 3, 3)
     ϕ² = ϕ.^2
 
diff --git a/verification/stratified_couette_flow/stratified_couette_flow.jl b/verification/stratified_couette_flow/stratified_couette_flow.jl
index e856cedc39..fd1822c513 100644
--- a/verification/stratified_couette_flow/stratified_couette_flow.jl
+++ b/verification/stratified_couette_flow/stratified_couette_flow.jl
@@ -108,7 +108,7 @@ function simulate_stratified_couette_flow(; Nxy, Nz, arch=GPU(), h=1, U_wall=1,
 
     model = Model(
        architecture = arch,
-               grid = RegularCartesianGrid(N = (Nxy, Nxy, Nz), L = (4π*h, 2π*h, 2h)),
+               grid = RegularCartesianGrid(size = (Nxy, Nxy, Nz), length = (4π*h, 2π*h, 2h)),
            buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=1.0, β=0.0)),
             closure = AnisotropicMinimumDissipation(ν=ν, κ=κ),
 boundary_conditions = HorizontallyPeriodicSolutionBCs(u=ubcs, v=vbcs, T=Tbcs),
@@ -257,4 +257,3 @@ boundary_conditions = HorizontallyPeriodicSolutionBCs(u=ubcs, v=vbcs, T=Tbcs),
                 wizard.Δt, prettytime(walltime / Ni))
     end
 end
-