diff --git a/examples/barotropicqg_acc.jl b/examples/barotropicqg_acc.jl
index a5a8818b..6b65f41c 100644
--- a/examples/barotropicqg_acc.jl
+++ b/examples/barotropicqg_acc.jl
@@ -4,13 +4,17 @@
 #md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/barotropicqgtopography.ipynb).
 #
 # An idealized version of the Southern Ocean. We solve the barotropic 
-# quasi-geostrophic eddy dynamics in a flud with variable depth $H-h(x,y)$. We 
+# quasi-geostrophic eddy dynamics in a flud with variable depth ``H-h(x,y)``. We 
 # also include an "Antarctic Circumpolar Current," i.e.,  a domain-average zonal 
-# velocity $U(t)$ which is forced by constant wind  stress $F$ and influenced by 
+# velocity ``U(t)`` which is forced by constant wind  stress ``F`` and influenced by 
 # bottom drag and topographic form stress. The equations solved are: 
-# $\partial_t \nabla^2\psi + \mathsf{J}(\psi-U y, \nabla^2\psi + \beta y + \eta) = -\mu\nabla^2\psi$
-# and 
-# $\partial_t U = F - \mu U -  \langle\psi\partial_x\eta\rangle$.
+# ```math
+# \partial_t \nabla^2 \psi + \mathsf{J}( \psi - U y, \nabla^2 \psi + \beta y + \eta ) = - \mu \nabla^2 \psi \, ,
+# ```
+# and
+# ```math 
+# \partial_t U = F - \mu U -  \langle \psi \partial_x \eta \rangle \, .
+# ```
 
 
 using FourierFlows, Plots, Printf
@@ -28,11 +32,11 @@ nothing # hide
 
 
 # ## Numerical parameters and time-stepping parameters
-     nx = 128            # 2D resolution = nx^2
-stepper = "FilteredRK4"  # timestepper
-    dt  = 0.1            # timestep
- nsteps = 10000          # total number of time-steps
- nsubs  = 25             # number of time-steps for intermediate logging/plotting (nsteps must be multiple of nsubs)
+      n = 128               # 2D resolution = n²
+stepper = "FilteredETDRK4"  # timestepper
+     dt = 0.1               # timestep
+ nsteps = 8000              # total number of time-steps
+ nsubs  = 25                # number of time-steps for intermediate logging/plotting (nsteps must be multiple of nsubs)
 nothing # hide
 
 
@@ -46,31 +50,31 @@ nν = 4         # viscosity order
  F = 0.0012    # normalized wind stress forcing on domain-averaged zonal flow U(t) flow
 nothing # hide
 
-# Define the topographic potential vorticity, $f_0 h(x, y)/H$
-topoPV(x, y) = 2cos(4x)*cos(4y)
+# Define the topographic potential vorticity, ``f_0 h(x, y)/H``,
+topoPV(x, y) = 2 * cos(4x) * cos(4y)
 nothing # hide
 
-# and the forcing function $F$ (here forcing is constant in time) that acts on the domain-averaged $U$ equation.
+# and the forcing function ``F`` (here forcing is constant in time) that acts on the domain-averaged ``U`` equation.
 calcFU(t) = F
 nothing # hide
 
 
 # ## Problem setup
 # We initialize a `Problem` by providing a set of keyword arguments,
-prob = BarotropicQG.Problem(dev; nx=nx, Lx=Lx, β=β, eta=topoPV,
+prob = BarotropicQG.Problem(dev; nx=n, Lx=Lx, β=β, eta=topoPV,
                   calcFU=calcFU, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper)
 nothing # hide
 
 # and define some shortcuts.
-sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
-x, y = gr.x, gr.y
+sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
+x, y = grid.x, grid.y
 nothing # hide
 
 
 # ## Setting initial conditions
 
 # Our initial condition is simply fluid at rest.
-BarotropicQG.set_zeta!(prob, zeros(gr.nx, gr.ny))
+BarotropicQG.set_zeta!(prob, zeros(grid.nx, grid.ny))
 
 
 # ## Diagnostics
@@ -100,7 +104,7 @@ nothing # hide
 
 # and then create Output.
 get_sol(prob) = sol # extracts the Fourier-transformed solution
-get_u(prob) = irfft(im*gr.lr.*gr.invKrsq.*sol, gr.nx)
+get_u(prob) = irfft(im * grid.lr .* grid.invKrsq .* sol, grid.nx)
 out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
 nothing # hide
 
@@ -111,13 +115,15 @@ nothing # hide
 # of energy and enstrophy.
 
 function plot_output(prob)
-
-  pq = heatmap(x, y, vs.q,
+  
+  grid = prob.grid
+  
+  pq = heatmap(grid.x, grid.y, vars.q',
                  c = :balance,
               clim = (-2, 2),
        aspectratio = 1,
-             xlims = (-gr.Lx/2, gr.Lx/2),
-             ylims = (-gr.Ly/2, gr.Ly/2),
+             xlims = (-grid.Lx/2, grid.Lx/2),
+             ylims = (-grid.Ly/2, grid.Ly/2),
             xticks = -3:3,
             yticks = -3:3,
             xlabel = "x",
@@ -129,7 +135,7 @@ function plot_output(prob)
              label = ["eddy energy" "mean energy"],
          linewidth = 2,
              alpha = 0.7,
-             xlims = (-0.1, 10.1),
+             xlims = (-0.1, 1.01 * μ * nsteps * dt),
              ylims = (0, 0.0008),
             xlabel = "μt")
           
@@ -137,11 +143,11 @@ function plot_output(prob)
              label = ["eddy enstrophy" "mean enstrophy"],
          linewidth = 2,
              alpha = 0.7,
-             xlims = (-0.1, 10.1),
+             xlims = (-0.1, 1.01 * μ * nsteps * dt),
              ylims = (-0.02, 0.12),
             xlabel = "μt")
 
-  l = @layout [ a{0.5w} grid(2, 1) ]
+  l = @layout [a{0.5w} Plots.grid(2, 1)]
   p = plot(pq, pE, pQ, layout=l, size = (900, 600))
 
   return p
@@ -157,34 +163,35 @@ p = plot_output(prob)
 
 startwalltime = time()
 
-anim = @animate for j=0:Int(nsteps/nsubs)
+anim = @animate for j = 0:round(Int, nsteps / nsubs)
   
-  cfl = cl.dt*maximum([maximum(vs.U.+vs.u)/gr.dx, maximum(vs.v)/gr.dy])
+  if j % (2000 / nsubs) == 0
+    cfl = clock.dt * maximum([maximum(vars.U .+ vars.u) / grid.dx, maximum(vars.v) / grid.dy])
+    
+    log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min",
+      clock.step, clock.t, cfl, E.data[E.i], Q.data[Q.i], (time()-startwalltime)/60)
+
+    println(log)
+  end
   
-  log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min",
-    cl.step, cl.t, cfl, E.data[E.i], Q.data[Q.i], (time()-startwalltime)/60)
-  
-  if j%(2000/nsubs)==0; println(log) end
-  
-  p[1][1][:z] = Array(vs.q)
-  p[1][:title] = "∇²ψ + η, μt="*@sprintf("%.2f", μ*cl.t)
-  push!(p[2][1], μ*E.t[E.i], E.data[E.i])
-  push!(p[2][2], μ*Emean.t[Emean.i], Emean.data[Emean.i])
-  push!(p[3][1], μ*Q.t[Q.i], Q.data[Q.i])
-  push!(p[3][2], μ*Qmean.t[Qmean.i], Qmean.data[Qmean.i])
+  p[1][1][:z] = vars.q
+  p[1][:title] = "∇²ψ + η, μt = "*@sprintf("%.2f", μ * clock.t)
+  push!(p[2][1], μ * E.t[E.i], E.data[E.i])
+  push!(p[2][2], μ * Emean.t[Emean.i], Emean.data[Emean.i])
+  push!(p[3][1], μ * Q.t[Q.i], Q.data[Q.i])
+  push!(p[3][2], μ * Qmean.t[Qmean.i], Qmean.data[Qmean.i])
 
   stepforward!(prob, diags, nsubs)
   BarotropicQG.updatevars!(prob)
-
 end
 
 mp4(anim, "barotropicqg_acc.mp4", fps=18)
 
-# Note that since mean flow enstrophy is $Q_U = \beta U$ it can attain negative values. 
+# Note that since mean flow enstrophy is ``Q_U = \beta U`` it can attain negative values. 
 
 
 # ## Save
 
 # Finally save the last snapshot.
-savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
+savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), clock.step)
 savefig(savename)
diff --git a/examples/barotropicqg_betadecay.jl b/examples/barotropicqg_betadecay.jl
index 851dc244..e86e0a76 100644
--- a/examples/barotropicqg_betadecay.jl
+++ b/examples/barotropicqg_betadecay.jl
@@ -22,7 +22,7 @@ nothing # hide
 
 # ## Numerical parameters and time-stepping parameters
 
-      n = 128            # 2D resolution = n^2
+      n = 128            # 2D resolution = n²
 stepper = "FilteredRK4"  # timestepper
      dt = 0.05           # timestep
  nsteps = 2000           # total number of time-steps
@@ -39,48 +39,48 @@ nothing # hide
 
 # ## Problem setup
 # We initialize a `Problem` by providing a set of keyword arguments. Not providing
-# a viscosity coefficient ν leads to the module's default value: ν=0. In this
-# example numerical instability due to accumulation of enstrophy in high wavenumbers
+# a viscosity coefficient `ν` leads to the module's default value: `ν=0`. In this
+# example numerical instability due to accumulation of enstrophy at high wavenumbers
 # is taken care with the `FilteredTimestepper` we picked. 
 prob = BarotropicQG.Problem(dev; nx=n, Lx=L, β=β, μ=μ, dt=dt, stepper=stepper)
 nothing # hide
 
 # and define some shortcuts
-sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
-x, y = gr.x, gr.y
+sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
+x, y = grid.x, grid.y
 nothing # hide
 
 
 # ## Setting initial conditions
 
 # Our initial condition consist of a flow that has power only at wavenumbers with
-# $8 < \frac{L}{2\pi} \sqrt{k_x^2+k_y^2} < 10$ and initial energy $E_0$:
+# ``8 < \frac{L}{2\pi} \sqrt{k_x^2 + k_y^2} < 10`` and initial energy ``E_0``:
 
-E0 = 0.1 # energy of initial condition
+E₀ = 0.1 # energy of initial condition
 
-K = @. sqrt(gr.Krsq)                      # a 2D array with the total wavenumber
-k = [gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]  # a 2D array with the zonal wavenumber
+K = @. sqrt(grid.Krsq)                          # a 2D array with the total wavenumber
+k = [grid.kr[i] for i=1:grid.nkr, j=1:grid.nl]  # a 2D array with the zonal wavenumber
 
 Random.seed!(1234)
-qih = randn(Complex{eltype(gr)}, size(sol))
-qih[K .<  8 * 2π/L] .= 0
-qih[K .> 10 * 2π/L] .= 0
-qih[k .== 0] .= 0         # no power at zonal wavenumber k=0 component
-Ein = energy(qih, gr)     # compute energy of qi
-qih *= sqrt(E0 / Ein)     # normalize qi to have energy E0
-qi  = irfft(qih, gr.nx)
+qih = randn(Complex{eltype(grid)}, size(sol))
+@. qih = ifelse(K < 2  * 2π/L, 0, qih)
+@. qih = ifelse(K > 10 * 2π/L, 0, qih)
+@. qih = ifelse(k == 0 * 2π/L, 0, qih)   # no power at zonal wavenumber k=0 component
+Ein = energy(qih, grid)           # compute energy of qi
+qih *= sqrt(E₀ / Ein)             # normalize qi to have energy E₀
+qi = irfft(qih, grid.nx)
 
 BarotropicQG.set_zeta!(prob, qi)
 nothing #hide
 
 # Let's plot the initial vorticity field:
 
-p1 = heatmap(x, y, vs.q,
+p1 = heatmap(x, y, vars.q',
          aspectratio = 1,
               c = :balance,
            clim = (-12, 12),
-          xlims = (-gr.Lx/2, gr.Lx/2),
-          ylims = (-gr.Ly/2, gr.Ly/2),
+          xlims = (-grid.Lx/2, grid.Lx/2),
+          ylims = (-grid.Ly/2, grid.Ly/2),
          xticks = -3:3,
          yticks = -3:3,
          xlabel = "x",
@@ -88,13 +88,13 @@ p1 = heatmap(x, y, vs.q,
           title = "initial vorticity ζ=∂v/∂x-∂u/∂y",
      framestyle = :box)
 
-p2 = contourf(x, y, vs.psi,
+p2 = contourf(x, y, vars.psi',
         aspectratio = 1,
              c = :viridis,
         levels = range(-0.65, stop=0.65, length=10), 
           clim = (-0.65, 0.65),
-         xlims = (-gr.Lx/2, gr.Lx/2),
-         ylims = (-gr.Ly/2, gr.Ly/2),
+         xlims = (-grid.Lx/2, grid.Lx/2),
+         ylims = (-grid.Ly/2, grid.Ly/2),
         xticks = -3:3,
         yticks = -3:3,
         xlabel = "x",
@@ -102,7 +102,7 @@ p2 = contourf(x, y, vs.psi,
          title = "initial streamfunction ψ",
     framestyle = :box)
 
-l = @layout grid(1, 2)
+l = @layout Plots.grid(1, 2)
 p = plot(p1, p2, layout=l, size=(900, 400))
 
 
@@ -131,14 +131,14 @@ nothing # hide
 
 # and then create Output.
 get_sol(prob) = sol # extracts the Fourier-transformed solution
-get_u(prob) = irfft(im*gr.l.*gr.invKrsq.*sol, gr.nx)
+get_u(prob) = irfft(im * grid.l .* grid.invKrsq .* sol, grid.nx)
 out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
 nothing # hide
 
 
 # ## Visualizing the simulation
 
-# We define a function that plots the vorticity and streamfunction fields and 
+# We define a function that plots the vorticity and streamfunction and 
 # their corresponding zonal mean structure.
 
 function plot_output(prob)
@@ -147,13 +147,13 @@ function plot_output(prob)
   ζ̄ = mean(ζ, dims=1)'
   ū = mean(prob.vars.u, dims=1)'
 
-  pζ = heatmap(x, y, ζ,
+  pζ = heatmap(x, y, ζ',
        aspectratio = 1,
             legend = false,
                  c = :balance,
               clim = (-12, 12),
-             xlims = (-gr.Lx/2, gr.Lx/2),
-             ylims = (-gr.Ly/2, gr.Ly/2),
+             xlims = (-grid.Lx/2, grid.Lx/2),
+             ylims = (-grid.Ly/2, grid.Ly/2),
             xticks = -3:3,
             yticks = -3:3,
             xlabel = "x",
@@ -161,14 +161,14 @@ function plot_output(prob)
              title = "vorticity ζ=∂v/∂x-∂u/∂y",
         framestyle = :box)
 
-  pψ = contourf(x, y, ψ,
+  pψ = contourf(x, y, ψ',
        aspectratio = 1,
             legend = false,
                  c = :viridis,
             levels = range(-0.65, stop=0.65, length=10), 
               clim = (-0.65, 0.65),
-             xlims = (-gr.Lx/2, gr.Lx/2),
-             ylims = (-gr.Ly/2, gr.Ly/2),
+             xlims = (-grid.Lx/2, grid.Lx/2),
+             ylims = (-grid.Ly/2, grid.Ly/2),
             xticks = -3:3,
             yticks = -3:3,
             xlabel = "x",
@@ -196,7 +196,7 @@ function plot_output(prob)
             ylabel = "y")
   plot!(pum, 0*y, y, linestyle=:dash, linecolor=:black)
 
-  l = @layout grid(2, 2)
+  l = @layout Plots.grid(2, 2)
   p = plot(pζ, pζm, pψ, pum, layout = l, size = (900, 800))
   
   return p
@@ -212,29 +212,32 @@ startwalltime = time()
 
 p = plot_output(prob)
 
-anim = @animate for j=0:Int(nsteps/nsubs)
+anim = @animate for j = 0:round(Int, nsteps/nsubs)
 
-  log = @sprintf("step: %04d, t: %d, E: %.4f, Q: %.4f, walltime: %.2f min",
-    cl.step, cl.t, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)
+  if j % (1000 / nsubs) == 0
+    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])
 
-  if j%(1000/nsubs)==0; println(log) end  
+    log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min",
+      clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)
 
+    println(log)
+  end  
 
-  p[1][1][:z] = Array(vs.zeta)
-  p[1][:title] = "vorticity, t="*@sprintf("%.2f", cl.t)
-  p[3][1][:z] = Array(vs.psi)
-  p[2][1][:x] = mean(vs.zeta, dims=1)'
-  p[4][1][:x] = mean(vs.u, dims=1)'
+  p[1][1][:z] = vars.zeta
+  p[1][:title] = "vorticity, t="*@sprintf("%.2f", clock.t)
+  p[3][1][:z] = vars.psi
+  p[2][1][:x] = mean(vars.zeta, dims=1)'
+  p[4][1][:x] = mean(vars.u, dims=1)'
 
   stepforward!(prob, diags, nsubs)
   BarotropicQG.updatevars!(prob)
 
 end
 
-mp4(anim, "barotropicqg_betadecay.mp4", fps=14)
+mp4(anim, "barotropicqg_betadecay.mp4", fps=12)
 
 # ## Save
 
 # Finally save the last snapshot.
-savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
+savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), clock.step)
 savefig(savename)
diff --git a/examples/barotropicqg_betaforced.jl b/examples/barotropicqg_betaforced.jl
index 61069419..a09700ee 100644
--- a/examples/barotropicqg_betaforced.jl
+++ b/examples/barotropicqg_betaforced.jl
@@ -45,24 +45,24 @@ nothing # hide
 #
 # We force the vorticity equation with stochastic excitation that is delta-correlated
 # in time and while spatially homogeneously and isotropically correlated. The forcing
-# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and 
-# width $\delta k_f$, and it injects energy per unit area and per unit time equal 
-# to $\varepsilon$.
+# has a spectrum with power in a ring in wavenumber space of radius ``k_f`` and 
+# width ``\delta k_f``, and it injects energy per unit area and per unit time equal 
+# to ``\varepsilon``.
 
-forcing_wavenumber = 14.0    # the central forcing wavenumber for a spectrum that is a ring in wavenumber space
-forcing_bandwidth  = 1.5     # the width of the forcing spectrum 
-ε = 0.001                    # energy input rate by the forcing
+forcing_wavenumber = 14.0 * 2π/L  # the central forcing wavenumber for a spectrum that is a ring in wavenumber space
+forcing_bandwidth  = 1.5  * 2π/L  # the width of the forcing spectrum 
+ε = 0.001                         # energy input rate by the forcing
 
-gr  = TwoDGrid(n, L)
+grid = TwoDGrid(n, L)
 
-K  = @. sqrt(gr.Krsq)                      # a 2D array with the total wavenumber
-k = [gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]   # a 2D array with the zonal wavenumber
+K  = @. sqrt(grid.Krsq)                          # a 2D array with the total wavenumber
+k = [grid.kr[i] for i=1:grid.nkr, j=1:grid.nl]   # a 2D array with the zonal wavenumber
 
 forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2))
-forcing_spectrum[K .<  2 * 2π/L] .= 0    # no power at low wavenumbers
-forcing_spectrum[K .> 20 * 2π/L] .= 0    # no power at high wavenumbers
-forcing_spectrum[k .< 2π/L] .= 0         # make sure forcing does not have power at k=0
-ε0 = parsevalsum(forcing_spectrum .* gr.invKrsq / 2, gr) / (gr.Lx * gr.Ly)
+@. forcing_spectrum = ifelse(K < 2  * 2π/L, 0, forcing_spectrum)      # no power at low wavenumbers
+@. forcing_spectrum = ifelse(K > 20 * 2π/L, 0, forcing_spectrum)      # no power at high wavenumbers
+@. forcing_spectrum = ifelse(k < 2π/L, 0, forcing_spectrum)    # make sure forcing does not have power at k=0
+ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly)
 @. forcing_spectrum *= ε/ε0               # normalize forcing to inject energy at rate ε
 
 seed!(1234) # reset of the random number generator for reproducibility
@@ -71,8 +71,9 @@ nothing # hide
 # Next we construct function `calcF!` that computes a forcing realization every timestep
 function calcFq!(Fh, sol, t, clock, vars, params, grid)
   ξ = ArrayType(dev)(exp.(2π * im * rand(eltype(grid), size(sol))) / sqrt(clock.dt))
-  @. Fh = ξ*sqrt.(forcing_spectrum)
-  Fh[abs.(grid.Krsq) .== 0] .= 0
+  @. Fh = ξ * sqrt.(forcing_spectrum)
+  @. Fh = ifelse(abs(grid.Krsq) == 0, 0, Fh)
+
   return nothing
 end
 nothing # hide
@@ -88,20 +89,20 @@ prob = BarotropicQG.Problem(dev; nx=n, Lx=L, β=β, μ=μ, dt=dt, stepper=steppe
 nothing # hide
 
 # Let's define some shortcuts.
-sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
-x, y = gr.x, gr.y
+sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
+x, y = grid.x, grid.y
 nothing # hide
 
 
 # First let's see how a forcing realization looks like.
-calcFq!(vs.Fqh, sol, 0.0, cl, vs, pr, gr)
+calcFq!(vars.Fqh, sol, 0.0, clock, vars, params, grid)
 
-heatmap(x, y, irfft(vs.Fqh, gr.nx),
+heatmap(x, y, irfft(vars.Fqh, grid.nx)',
      aspectratio = 1,
                c = :balance,
             clim = (-8, 8),
-           xlims = (-gr.Lx/2, gr.Lx/2),
-           ylims = (-gr.Ly/2, gr.Ly/2),
+           xlims = (-grid.Lx/2, grid.Lx/2),
+           ylims = (-grid.Ly/2, grid.Ly/2),
           xticks = -3:3,
           yticks = -3:3,
           xlabel = "x",
@@ -113,7 +114,7 @@ heatmap(x, y, irfft(vs.Fqh, gr.nx),
 # ## Setting initial conditions
 
 # Our initial condition is simply fluid at rest.
-BarotropicQG.set_zeta!(prob, zeros(gr.nx, gr.ny))
+BarotropicQG.set_zeta!(prob, zeros(grid.nx, grid.ny))
 
 
 # ## Diagnostics
@@ -141,7 +142,7 @@ nothing # hide
 
 # and then create Output.
 get_sol(prob) = sol # extracts the Fourier-transformed solution
-get_u(prob) = irfft(im*gr.l.*gr.invKrsq.*sol, gr.nx)
+get_u(prob) = irfft(im * grid.l .* grid.invKrsq .* sol, grid.nx)
 out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
 nothing # hide
 
@@ -157,13 +158,13 @@ function plot_output(prob)
   ζ̄ = mean(ζ, dims=1)'
   ū = mean(prob.vars.u, dims=1)'
   
-  pζ = heatmap(x, y, ζ,
+  pζ = heatmap(x, y, ζ',
        aspectratio = 1,
             legend = false,
                  c = :balance,
               clim = (-8, 8),
-             xlims = (-gr.Lx/2, gr.Lx/2),
-             ylims = (-gr.Ly/2, gr.Ly/2),
+             xlims = (-grid.Lx/2, grid.Lx/2),
+             ylims = (-grid.Ly/2, grid.Ly/2),
             xticks = -3:3,
             yticks = -3:3,
             xlabel = "x",
@@ -171,15 +172,15 @@ function plot_output(prob)
              title = "vorticity ζ=∂v/∂x-∂u/∂y",
         framestyle = :box)
 
-  pψ = contourf(x, y, ψ,
+  pψ = contourf(x, y, ψ',
             levels = -0.32:0.04:0.32,
        aspectratio = 1,
          linewidth = 1,
             legend = false,
               clim = (-0.22, 0.22),
                  c = :viridis,
-             xlims = (-gr.Lx/2, gr.Lx/2),
-             ylims = (-gr.Ly/2, gr.Ly/2),
+             xlims = (-grid.Lx/2, grid.Lx/2),
+             ylims = (-grid.Ly/2, grid.Ly/2),
             xticks = -3:3,
             yticks = -3:3,
             xlabel = "x",
@@ -225,7 +226,7 @@ function plot_output(prob)
              ylims = (0, 2.5),
             xlabel = "μt")
 
-  l = @layout grid(2, 3)
+  l = @layout Plots.grid(2, 3)
   p = plot(pζ, pζm, pE, pψ, pum, pZ, layout=l, size = (1000, 600), dpi=150)
 
   return p
@@ -241,27 +242,27 @@ startwalltime = time()
 
 p = plot_output(prob)
 
-anim = @animate for j=0:Int(nsteps/nsubs)
-        
-  cfl = cl.dt*maximum([maximum(vs.u)/gr.dx, maximum(vs.v)/gr.dy])
+anim = @animate for j = 0:Int(nsteps / nsubs)
   
-  log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min",
-  cl.step, cl.t, cfl, E.data[E.i], Z.data[Z.i], 
-  (time()-startwalltime)/60)
+  if j % (1000 / nsubs) == 0
+    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])
   
-  if j%(1000/nsubs)==0; println(log) end  
+    log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min", 
+    clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)
+
+    println(log)
+  end  
   
-  p[1][1][:z] = Array(vs.zeta)
-  p[1][:title] = "vorticity, μt="*@sprintf("%.2f", μ*cl.t)
-  p[4][1][:z] = Array(vs.psi)
-  p[2][1][:x] = mean(vs.zeta, dims=1)'
-  p[5][1][:x] = mean(vs.u, dims=1)'
-  push!(p[3][1], μ*E.t[E.i], E.data[E.i])
-  push!(p[6][1], μ*Z.t[Z.i], Z.data[Z.i])
+  p[1][1][:z] = vars.zeta
+  p[1][:title] = "vorticity, μt="*@sprintf("%.2f", μ * clock.t)
+  p[4][1][:z] = vars.psi
+  p[2][1][:x] = mean(vars.zeta, dims=1)'
+  p[5][1][:x] = mean(vars.u, dims=1)'
+  push!(p[3][1], μ * E.t[E.i], E.data[E.i])
+  push!(p[6][1], μ * Z.t[Z.i], Z.data[Z.i])
   
   stepforward!(prob, diags, nsubs)
   BarotropicQG.updatevars!(prob)
-  
 end
 
 mp4(anim, "barotropicqg_betaforced.mp4", fps=18)
@@ -270,5 +271,5 @@ mp4(anim, "barotropicqg_betaforced.mp4", fps=18)
 # ## Save
 
 # Finally save the last snapshot.
-savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
+savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), clock.step)
 savefig(savename)
diff --git a/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl
index 265bc92c..ca1353a1 100644
--- a/examples/barotropicqgql_betaforced.jl
+++ b/examples/barotropicqgql_betaforced.jl
@@ -46,24 +46,24 @@ nothing # hide
 
 # We force the vorticity equation with stochastic excitation that is delta-correlated
 # in time and while spatially homogeneously and isotropically correlated. The forcing
-# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and 
-# width $\delta k_f$, and it injects energy per unit area and per unit time equal 
-# to $\varepsilon$.
+# has a spectrum with power in a ring in wavenumber space of radius ``k_f`` and 
+# width ``\delta k_f``, and it injects energy per unit area and per unit time equal 
+# to ``\varepsilon``.
 
 forcing_wavenumber = 14.0    # the central forcing wavenumber for a spectrum that is a ring in wavenumber space
 forcing_bandwidth  = 1.5     # the width of the forcing spectrum 
 ε = 0.001                    # energy input rate by the forcing
 
-gr  = TwoDGrid(n, L)
+grid = TwoDGrid(n, L)
 
-K = @. sqrt(gr.Krsq)                      # a 2D array with the total wavenumber
-k = [gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]  # a 2D array with the zonal wavenumber
+K = @. sqrt(grid.Krsq)                          # a 2D array with the total wavenumber
+k = [grid.kr[i] for i=1:grid.nkr, j=1:grid.nl]  # a 2D array with the zonal wavenumber
 
 forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2))
-forcing_spectrum[K .<  2 * 2π/L] .= 0    # no power at low wavenumbers
-forcing_spectrum[K .> 20 * 2π/L] .= 0    # no power at high wavenumbers
-forcing_spectrum[k .< 2π/L ] .= 0        # make sure forcing does not have power at k=0
-ε0 = parsevalsum(forcing_spectrum .* gr.invKrsq / 2, gr) / (gr.Lx * gr.Ly)
+@. forcing_spectrum = ifelse(K < 2  * 2π/L, 0, forcing_spectrum)      # no power at low wavenumbers
+@. forcing_spectrum = ifelse(K > 20 * 2π/L, 0, forcing_spectrum)      # no power at high wavenumbers
+@. forcing_spectrum = ifelse(k < 2π/L, 0, forcing_spectrum)    # make sure forcing does not have power at k=0
+ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly)
 @. forcing_spectrum *= ε/ε0               # normalize forcing to inject energy at rate ε
 
 seed!(1234) # reset of the random number generator for reproducibility
@@ -73,7 +73,8 @@ nothing # hide
 function calcF!(Fh, sol, t, clock, vars, params, grid)
   ξ = ArrayType(dev)(exp.(2π * im * rand(eltype(grid), size(sol))) / sqrt(clock.dt))
   @. Fh = ξ * sqrt.(forcing_spectrum)
-  Fh[abs.(grid.Krsq) .== 0] .= 0
+  @. Fh = ifelse(abs(grid.Krsq) == 0, 0, Fh)
+
   return nothing
 end
 nothing # hide
@@ -89,20 +90,20 @@ prob = BarotropicQGQL.Problem(dev; nx=n, Lx=L, β=β, μ=μ, dt=dt, stepper=step
 nothing # hide
 
 # and define some shortcuts.
-sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
-x, y = gr.x, gr.y
+sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
+x, y = grid.x, grid.y
 nothing # hide
 
 
 # First let's see how a forcing realization looks like.
-calcF!(vs.Fh, sol, 0.0, cl, vs, pr, gr)
+calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid)
 
-heatmap(x, y, irfft(vs.Fh, gr.nx),
+heatmap(x, y, irfft(vars.Fh, grid.nx)',
      aspectratio = 1,
                c = :balance,
             clim = (-8, 8),
-           xlims = (-gr.Lx/2, gr.Lx/2),
-           ylims = (-gr.Ly/2, gr.Ly/2),
+           xlims = (-grid.Lx/2, grid.Lx/2),
+           ylims = (-grid.Ly/2, grid.Ly/2),
           xticks = -3:3,
           yticks = -3:3,
           xlabel = "x",
@@ -114,7 +115,7 @@ heatmap(x, y, irfft(vs.Fh, gr.nx),
 # ## Setting initial conditions
 
 # Our initial condition is simply fluid at rest.
-BarotropicQGQL.set_zeta!(prob, zeros(gr.nx, gr.ny))
+BarotropicQGQL.set_zeta!(prob, zeros(grid.nx, grid.ny))
 nothing # hide
 
 # ## Diagnostics
@@ -154,7 +155,7 @@ nothing # hide
 
 # and then create Output.
 get_sol(prob) = sol # extracts the Fourier-transformed solution
-get_u(prob) = irfft(im*g.l.*g.invKrsq.*sol, g.nx)
+get_u(prob) = irfft(im * g.l .* g.invKrsq .* sol, g.nx)
 out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
 nothing # hide
 
@@ -173,13 +174,13 @@ function plot_output(prob)
   ζ̄ₘ = mean(ζ̄, dims=1)'
   ūₘ = mean(prob.vars.U, dims=1)'
 
-  pζ = heatmap(x, y, ζ,
+  pζ = heatmap(x, y, ζ',
        aspectratio = 1,
             legend = false,
                  c = :balance,
               clim = (-8, 8),
-             xlims = (-gr.Lx/2, gr.Lx/2),
-             ylims = (-gr.Ly/2, gr.Ly/2),
+             xlims = (-grid.Lx/2, grid.Lx/2),
+             ylims = (-grid.Ly/2, grid.Ly/2),
             xticks = -3:3,
             yticks = -3:3,
             xlabel = "x",
@@ -187,15 +188,15 @@ function plot_output(prob)
              title = "vorticity ζ=∂v/∂x-∂u/∂y",
         framestyle = :box)
 
-  pψ = contourf(x, y, ψ,
+  pψ = contourf(x, y, ψ',
             levels = -0.32:0.04:0.32,
        aspectratio = 1,
          linewidth = 1,
             legend = false,
               clim = (-0.22, 0.22),
                  c = :viridis,
-             xlims = (-gr.Lx/2, gr.Lx/2),
-             ylims = (-gr.Ly/2, gr.Ly/2),
+             xlims = (-grid.Lx/2, grid.Lx/2),
+             ylims = (-grid.Ly/2, grid.Ly/2),
             xticks = -3:3,
             yticks = -3:3,
             xlabel = "x",
@@ -241,7 +242,7 @@ function plot_output(prob)
              ylims = (0, 5),
             xlabel = "μt")
 
-  l = @layout grid(2, 3)
+  l = @layout Plots.grid(2, 3)
   p = plot(pζ, pζm, pE, pψ, pum, pZ, layout=l, size = (1000, 600), dpi=150)
   
   return p
@@ -255,27 +256,28 @@ p = plot_output(prob)
 
 startwalltime = time()
 
-anim = @animate for j=0:Int(nsteps/nsubs)
+anim = @animate for j = 0:round(Int, nsteps / nsubs)
 
-  cfl = cl.dt*maximum([maximum(vs.v)/gr.dy, maximum(vs.u+vs.U)/gr.dx])
+  if j % (1000 / nsubs) == 0
+    cfl = clock.dt * maximum([maximum(vars.v) / grid.dy, maximum(vars.u .+ vars.U) / grid.dx])
 
-  log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min",
-    cl.step, cl.t, cfl, E.data[E.i], Z.data[Z.i],
-    (time()-startwalltime)/60)
+    log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min",
+      clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i],
+      (time()-startwalltime)/60)
 
-  if j%(1000/nsubs)==0; println(log) end
+    println(log)
+  end
 
-  p[1][1][:z] = @. vs.zeta + vs.Zeta
-  p[1][:title] = "vorticity, μt="*@sprintf("%.2f", μ*cl.t)
-  p[4][1][:z] = @. vs.psi + vs.Psi
-  p[2][1][:x] = mean(vs.Zeta, dims=1)'
-  p[5][1][:x] = mean(vs.U, dims=1)'
-  push!(p[3][1], μ*E.t[E.i], E.data[E.i])
-  push!(p[6][1], μ*Z.t[Z.i], Z.data[Z.i])
+  p[1][1][:z] = @. vars.zeta + vars.Zeta
+  p[1][:title] = "vorticity, μt="*@sprintf("%.2f", μ * clock.t)
+  p[4][1][:z] = @. vars.psi + vars.Psi
+  p[2][1][:x] = mean(vars.Zeta, dims=1)'
+  p[5][1][:x] = mean(vars.U, dims=1)'
+  push!(p[3][1], μ * E.t[E.i], E.data[E.i])
+  push!(p[6][1], μ * Z.t[Z.i], Z.data[Z.i])
   
   stepforward!(prob, diags, nsubs)
   BarotropicQGQL.updatevars!(prob)
-
 end
 
 mp4(anim, "barotropicqgql_betaforced.mp4", fps=18)
@@ -284,5 +286,5 @@ mp4(anim, "barotropicqgql_betaforced.mp4", fps=18)
 # ## Save
 
 # Finally save the last snapshot.
-savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
+savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), clock.step)
 savefig(savename)
diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl
index 44172034..ae57cc08 100644
--- a/examples/multilayerqg_2layer.jl
+++ b/examples/multilayerqg_2layer.jl
@@ -4,7 +4,7 @@
 #md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/multilayerqg_2layer.ipynb).
 #
 # A simulation of the growth of barolinic instability in the Phillips 2-layer model
-# when we impose a vertical mean flow shear as a difference $\Delta U$ in the
+# when we impose a vertical mean flow shear as a difference ``\Delta U`` in the
 # imposed, domain-averaged, zonal flow at each layer.
 
 using FourierFlows, Plots, Printf
@@ -54,8 +54,8 @@ prob = MultilayerQG.Problem(nlayers, dev; nx=nx, Lx=Lx, f₀=f₀, g=g, H=H, ρ=
 nothing # hide
 
 # and define some shortcuts.
-sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid
-x, y = gr.x, gr.y
+sol, clock, params, vars, grid = prob.sol, prob.clock, prob.params, prob.vars, prob.grid
+x, y = grid.x, grid.y
 nothing # hide
 
 
@@ -66,7 +66,7 @@ nothing # hide
 
 q_i  = 4e-3randn((nx, ny, nlayers))
 qh_i = prob.timestepper.filter .* rfft(q_i, (1, 2)) # only apply rfft in dims=1, 2
-q_i  = irfft(qh_i, gr.nx, (1, 2)) # only apply irfft in dims=1, 2
+q_i  = irfft(qh_i, grid.nx, (1, 2)) # only apply irfft in dims=1, 2
 
 MultilayerQG.set_q!(prob, q_i)
 nothing # hide
@@ -97,12 +97,16 @@ nothing # hide
 # And then create Output
 get_sol(prob) = sol # extracts the Fourier-transformed solution
 function get_u(prob)
-  @. v.qh = sol
-  streamfunctionfrompv!(v.ψh, v.qh, p, g)
-  @. v.uh = -im * g.l * v.ψh
-  invtransform!(v.u, v.uh, p)
-  return v.u
+  sol, params, vars, grid = prob.sol, prob.params, prob.vars, prob.grid
+  
+  @. vars.qh = sol
+  streamfunctionfrompv!(vars.ψh, vars.qh, params, grid)
+  @. vars.uh = -im * grid.l * vars.ψh
+  invtransform!(vars.u, vars.uh, params)
+  
+  return vars.u
 end
+
 out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
 nothing #hide
 
@@ -115,17 +119,18 @@ nothing #hide
 symlims(data) = maximum(abs.(extrema(data))) |> q -> (-q, q)
 
 function plot_output(prob)
+  Lx, Ly = prob.grid.Lx, prob.grid.Ly
   
-  l = @layout grid(2, 3)
+  l = @layout Plots.grid(2, 3)
   p = plot(layout=l, size = (1000, 600), dpi=150)
   
   for m in 1:nlayers
-    heatmap!(p[(m-1)*3+1], x, y, vs.q[:, :, m],
+    heatmap!(p[(m-1) * 3 + 1], x, y, vars.q[:, :, m]',
          aspectratio = 1,
               legend = false,
                    c = :balance,
-               xlims = (-gr.Lx/2, gr.Lx/2),
-               ylims = (-gr.Ly/2, gr.Ly/2),
+               xlims = (-Lx / 2, Lx / 2),
+               ylims = (-Ly / 2, Ly / 2),
                clims = symlims,
               xticks = -3:3,
               yticks = -3:3,
@@ -134,13 +139,13 @@ function plot_output(prob)
                title = "q_"*string(m),
           framestyle = :box)
 
-    contourf!(p[(m-1)*3+2], x, y, vs.ψ[:, :, m],
+    contourf!(p[(m-1) * 3 + 2], x, y, vars.ψ[:, :, m]',
               levels = 8,
          aspectratio = 1,
               legend = false,
                    c = :viridis,
-               xlims = (-gr.Lx/2, gr.Lx/2),
-               ylims = (-gr.Ly/2, gr.Ly/2),
+               xlims = (-Lx / 2, Lx / 2),
+               ylims = (-Ly / 2, Ly / 2),
                clims = symlims,
               xticks = -3:3,
               yticks = -3:3,
@@ -185,22 +190,23 @@ p = plot_output(prob)
 
 startwalltime = time()
 
-anim = @animate for j=0:Int(nsteps/nsubs)
-  
-  cfl = cl.dt*maximum([maximum(vs.u)/gr.dx, maximum(vs.v)/gr.dy])
-  
-  log = @sprintf("step: %04d, t: %.1f, cfl: %.2f, KE1: %.4f, KE2: %.4f, PE: %.4f, walltime: %.2f min", cl.step, cl.t, cfl, E.data[E.i][1][1], E.data[E.i][1][2], E.data[E.i][2][1], (time()-startwalltime)/60)
+anim = @animate for j = 0:round(Int, nsteps / nsubs)
+  if j % (1000 / nsubs) == 0
+    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])
+    
+    log = @sprintf("step: %04d, t: %.1f, cfl: %.2f, KE1: %.3e, KE2: %.3e, PE: %.3e, walltime: %.2f min", clock.step, clock.t, cfl, E.data[E.i][1][1], E.data[E.i][1][2], E.data[E.i][2][1], (time()-startwalltime)/60)
 
-  if j%(1000/nsubs)==0; println(log) end
+    println(log)
+  end
   
   for m in 1:nlayers
-    p[(m-1)*3+1][1][:z] = @. vs.q[:, :, m]
-    p[(m-1)*3+2][1][:z] = @. vs.ψ[:, :, m]
+    p[(m-1) * 3 + 1][1][:z] = @. vars.q[:, :, m]
+    p[(m-1) * 3 + 2][1][:z] = @. vars.ψ[:, :, m]
   end
   
-  push!(p[3][1], μ*E.t[E.i], E.data[E.i][1][1])
-  push!(p[3][2], μ*E.t[E.i], E.data[E.i][1][2])
-  push!(p[6][1], μ*E.t[E.i], E.data[E.i][2][1])
+  push!(p[3][1], μ * E.t[E.i], E.data[E.i][1][1])
+  push!(p[3][2], μ * E.t[E.i], E.data[E.i][1][2])
+  push!(p[6][1], μ * E.t[E.i], E.data[E.i][2][1])
   
   stepforward!(prob, diags, nsubs)
   MultilayerQG.updatevars!(prob)
@@ -212,5 +218,5 @@ mp4(anim, "multilayerqg_2layer.mp4", fps=18)
 # ## Save
 
 # Finally save the last snapshot.
-savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
+savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), clock.step)
 savefig(savename)
diff --git a/examples/surfaceqg_decaying.jl b/examples/surfaceqg_decaying.jl
index 7c6c60ca..81304dff 100644
--- a/examples/surfaceqg_decaying.jl
+++ b/examples/surfaceqg_decaying.jl
@@ -171,27 +171,28 @@ startwalltime = time()
 
 p = plot_output(prob)
 
-anim = @animate for j=0:Int(nsteps/nsubs)
+anim = @animate for j = 0:round(Int, nsteps/nsubs)
+  if j % (500 / nsubs) == 0
+    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])
 
-  cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])
-
-  if j%(500/nsubs)==0
     log1 = @sprintf("step: %04d, t: %.1f, cfl: %.3f, walltime: %.2f min",
           clock.step, clock.t, cfl, (time()-startwalltime)/60)
+
     log2 = @sprintf("buoyancy variance: %.2e, buoyancy variance dissipation: %.2e",
               B.data[B.i], Dᵇ.data[Dᵇ.i])
+
     println(log1)
+
     println(log2)
   end
 
   p[1][1][:z] = vars.b
-  p[1][:title] = "buoyancy, t="*@sprintf("%.2f", clock.t)
+  p[1][:title] = "buoyancy, t=" * @sprintf("%.2f", clock.t)
   push!(p[2][1], KE.t[KE.i], KE.data[KE.i])
   push!(p[3][1], B.t[B.i], B.data[B.i])
 
   stepforward!(prob, diags, nsubs)
   SurfaceQG.updatevars!(prob)
-
 end
 
 mp4(anim, "sqg_ellipticalvortex.mp4", fps=14)
@@ -208,7 +209,7 @@ pu = heatmap(x, y, vars.u',
           yticks = -3:3,
           xlabel = "x",
           ylabel = "y",
-           title = "uₛ(x, y, t="*@sprintf("%.2f", clock.t)*")",
+           title = "uₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")",
       framestyle = :box)
 
 pv = heatmap(x, y, vars.v',
@@ -221,7 +222,7 @@ pv = heatmap(x, y, vars.v',
           yticks = -3:3,
           xlabel = "x",
           ylabel = "y",
-           title = "vₛ(x, y, t="*@sprintf("%.2f", clock.t)*")",
+           title = "vₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")",
       framestyle = :box)
 
 pb = heatmap(x, y, vars.b',
@@ -234,7 +235,7 @@ pb = heatmap(x, y, vars.b',
           yticks = -3:3,
           xlabel = "x",
           ylabel = "y",
-           title = "bₛ(x, y, t="*@sprintf("%.2f", clock.t)*")",
+           title = "bₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")",
       framestyle = :box)
 
 layout = @layout [a{0.5h}; b{0.5w} c{0.5w}]
diff --git a/examples/twodnavierstokes_decaying.jl b/examples/twodnavierstokes_decaying.jl
index 85ed2759..c6b65f43 100644
--- a/examples/twodnavierstokes_decaying.jl
+++ b/examples/twodnavierstokes_decaying.jl
@@ -42,8 +42,8 @@ prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ny=n, Ly=L, dt=dt, stepper="Fil
 nothing # hide
 
 # Next we define some shortcuts for convenience.
-sol, cl, vs, gr = prob.sol, prob.clock, prob.vars, prob.grid
-x, y = gr.x, gr.y
+sol, clock, vars, grid = prob.sol, prob.clock, prob.vars, prob.grid
+x, y = grid.x, grid.y
 nothing # hide
 
 
@@ -53,12 +53,12 @@ nothing # hide
 # in the paper by McWilliams (_JFM_, 1984)
 seed!(1234)
 k₀, E₀ = 6, 0.5
-ζ₀ = peakedisotropicspectrum(gr, k₀, E₀, mask=prob.timestepper.filter)
+ζ₀ = peakedisotropicspectrum(grid, k₀, E₀, mask=prob.timestepper.filter)
 TwoDNavierStokes.set_zeta!(prob, ζ₀)
 nothing # hide
 
 # Let's plot the initial vorticity field:
-heatmap(x, y, vs.zeta,
+heatmap(x, y, vars.zeta',
          aspectratio = 1,
               c = :balance,
            clim = (-40, 40),
@@ -96,8 +96,8 @@ if !isdir(plotpath); mkdir(plotpath); end
 nothing # hide
 
 # And then create Output
-get_sol(prob) = Array(prob.sol) # extracts the Fourier-transformed solution
-get_u(prob) = Array(irfft(im*gr.l.*gr.invKrsq.*sol, gr.nx))
+get_sol(prob) = prob.sol # extracts the Fourier-transformed solution
+get_u(prob) = irfft(im * grid.l .* grid.invKrsq .* sol, grid.nx)
 out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
 saveproblem(out)
 nothing # hide
@@ -108,7 +108,7 @@ nothing # hide
 # We initialize a plot with the vorticity field and the time-series of
 # energy and enstrophy diagnostics.
 
-p1 = heatmap(x, y, vs.zeta,
+p1 = heatmap(x, y, vars.zeta',
          aspectratio = 1,
                    c = :balance,
                 clim = (-40, 40),
@@ -118,7 +118,7 @@ p1 = heatmap(x, y, vs.zeta,
               yticks = -3:3,
               xlabel = "x",
               ylabel = "y",
-               title = "vorticity, t="*@sprintf("%.2f", cl.t),
+               title = "vorticity, t=" * @sprintf("%.2f", clock.t),
           framestyle = :box)
 
 p2 = plot(2, # this means "a plot with two series"
@@ -130,7 +130,7 @@ p2 = plot(2, # this means "a plot with two series"
                xlims = (0, 41),
                ylims = (0, 1.1))
 
-l = @layout grid(1, 2)
+l = @layout Plots.grid(1, 2)
 p = plot(p1, p2, layout = l, size = (900, 400))
 
 
@@ -141,20 +141,22 @@ p = plot(p1, p2, layout = l, size = (900, 400))
 startwalltime = time()
 
 anim = @animate for j = 0:Int(nsteps/nsubs)
+  if j % (1000 / nsubs) == 0
+    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])
     
-  log = @sprintf("step: %04d, t: %d, ΔE: %.4f, ΔZ: %.4f, walltime: %.2f min",
-      cl.step, cl.t, E.data[E.i]/E.data[1], Z.data[Z.i]/Z.data[1], (time()-startwalltime)/60)
-  
-  if j%(1000/nsubs)==0; println(log) end  
+    log = @sprintf("step: %04d, t: %d, cfl: %.2f, ΔE: %.4f, ΔZ: %.4f, walltime: %.2f min",
+        clock.step, clock.t, cfl, E.data[E.i]/E.data[1], Z.data[Z.i]/Z.data[1], (time()-startwalltime)/60)
 
-  p[1][1][:z] = vs.zeta
-  p[1][:title] = "vorticity, t="*@sprintf("%.2f", cl.t)
+    println(log)
+  end  
+
+  p[1][1][:z] = vars.zeta
+  p[1][:title] = "vorticity, t=" * @sprintf("%.2f", clock.t)
   push!(p[2][1], E.t[E.i], E.data[E.i]/E.data[1])
   push!(p[2][2], Z.t[Z.i], Z.data[Z.i]/Z.data[1])
 
   stepforward!(prob, diags, nsubs)
   TwoDNavierStokes.updatevars!(prob)  
-  
 end
 
 mp4(anim, "twodturb.mp4", fps=18)
@@ -169,9 +171,9 @@ saveoutput(out)
 # After the simulation is done we plot the radial energy spectrum to illustrate
 # how `FourierFlows.radialspectrum` can be used,
 
-E  = @. 0.5*(vs.u^2 + vs.v^2) # energy density
+E  = @. 0.5 * (vars.u^2 + vars.v^2) # energy density
 Eh = rfft(E)                  # Fourier transform of energy density
-kr, Ehr = FourierFlows.radialspectrum(Eh, gr, refinement=1) # compute radial specturm of `Eh`
+kr, Ehr = FourierFlows.radialspectrum(Eh, grid, refinement=1) # compute radial specturm of `Eh`
 nothing # hide
 
 # and we plot it.
@@ -179,8 +181,7 @@ plot(kr, abs.(Ehr),
     linewidth = 2,
         alpha = 0.7,
        xlabel = "kᵣ", ylabel = "∫ |Ê| kᵣ dk_θ",
-        xlims = (5e-1, gr.nx),
+        xlims = (5e-1, grid.nx),
        xscale = :log10, yscale = :log10,
         title = "Radial energy spectrum",
        legend = false)
-       
\ No newline at end of file
diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl
index 1995b431..ae682476 100644
--- a/examples/twodnavierstokes_stochasticforcing.jl
+++ b/examples/twodnavierstokes_stochasticforcing.jl
@@ -40,22 +40,22 @@ nothing # hide
 
 # We force the vorticity equation with stochastic excitation that is delta-correlated
 # in time and while spatially homogeneously and isotropically correlated. The forcing
-# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and
-# width $\delta k_f$, and it injects energy per unit area and per unit time equal
-# to $\varepsilon$.
+# has a spectrum with power in a ring in wavenumber space of radius ``k_f`` and
+# width ``\delta k_f``, and it injects energy per unit area and per unit time equal
+# to ``\varepsilon``.
 
-forcing_wavenumber = 14.0    # the central forcing wavenumber for a spectrum that is a ring in wavenumber space
-forcing_bandwidth  = 1.5     # the width of the forcing spectrum
-ε = 0.1                      # energy input rate by the forcing
+forcing_wavenumber = 14.0 * 2π/L   # the central forcing wavenumber for a spectrum that is a ring in wavenumber space
+forcing_bandwidth  = 1.5  * 2π/L   # the width of the forcing spectrum
+ε = 0.1                            # energy input rate by the forcing
 
 grid = TwoDGrid(dev, n, L)
 
 K = @. sqrt(grid.Krsq)
 forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2))
-forcing_spectrum[K .< ( 2 * 2π/L)] .= 0 # no power at low wavenumbers
-forcing_spectrum[K .> (20 * 2π/L)] .= 0 # no power at high wavenumbers
+@. forcing_spectrum = ifelse(K < 2  * 2π/L, 0, forcing_spectrum)      # no power at low wavenumbers
+@. forcing_spectrum = ifelse(K > 20 * 2π/L, 0, forcing_spectrum)      # no power at high wavenumbers
 ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly)
-@. forcing_spectrum *= ε/ε0             # normalize forcing to inject energy at rate ε
+@. forcing_spectrum *= ε / ε0             # normalize forcing to inject energy at rate ε
 
 seed!(1234)
 nothing # hide
@@ -65,7 +65,8 @@ function calcF!(Fh, sol, t, clock, vars, params, grid)
   ξ = ArrayType(dev)(exp.(2π * im * rand(eltype(grid), size(sol))) / sqrt(clock.dt))
   ξ[1, 1] = 0
   @. Fh = ξ * sqrt(forcing_spectrum)
-  nothing
+  
+  return nothing
 end
 nothing # hide
 
@@ -89,7 +90,7 @@ nothing # hide
 # go back to physical space.
 calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid)
 
-heatmap(x, y, irfft(vars.Fh, grid.nx),
+heatmap(x, y, irfft(vars.Fh, grid.nx)',
      aspectratio = 1,
                c = :balance,
             clim = (-200, 200),
@@ -124,7 +125,7 @@ nothing # hide
 # energy and enstrophy diagnostics. To plot energy and enstrophy on the same
 # axes we scale enstrophy with $k_f^2$.
 
-p1 = heatmap(x, y, vars.zeta,
+p1 = heatmap(x, y, vars.zeta',
          aspectratio = 1,
                    c = :balance,
                 clim = (-40, 40),
@@ -134,7 +135,7 @@ p1 = heatmap(x, y, vars.zeta,
               yticks = -3:3,
               xlabel = "x",
               ylabel = "y",
-               title = "vorticity, t="*@sprintf("%.2f", clock.t),
+               title = "vorticity, t=" * @sprintf("%.2f", clock.t),
           framestyle = :box)
 
 p2 = plot(2, # this means "a plot with two series"
@@ -155,21 +156,23 @@ p = plot(p1, p2, layout = l, size = (900, 420))
 # Finally, we time-step the `Problem` forward in time.
 
 startwalltime = time()
-anim = @animate for j = 0:Int(nsteps/nsubs)
+
+anim = @animate for j = 0:round(Int, nsteps / nsubs)
+  if j % (1000/nsubs) == 0
+    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])
     
-  log = @sprintf("step: %04d, t: %d, E: %.4f, Z: %.4f, walltime: %.2f min",
-      clock.step, clock.t, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)
-  
-  if j%(1000/nsubs)==0; println(log) end  
+    log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Z: %.4f, walltime: %.2f min",
+          clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)
+    println(log)
+  end  
 
   p[1][1][:z] = vars.zeta
-  p[1][:title] = "vorticity, μ t = "*@sprintf("%.2f", μ * clock.t)
+  p[1][:title] = "vorticity, μt = " * @sprintf("%.2f", μ * clock.t)
   push!(p[2][1], μ * E.t[E.i], E.data[E.i])
   push!(p[2][2], μ * Z.t[Z.i], Z.data[Z.i] / forcing_wavenumber^2)
 
   stepforward!(prob, diags, nsubs)
   TwoDNavierStokes.updatevars!(prob)  
-  
 end
 
 mp4(anim, "twodturb_forced.mp4", fps=18)
diff --git a/examples/twodnavierstokes_stochasticforcing_budgets.jl b/examples/twodnavierstokes_stochasticforcing_budgets.jl
index a4a630ff..3904902b 100644
--- a/examples/twodnavierstokes_stochasticforcing_budgets.jl
+++ b/examples/twodnavierstokes_stochasticforcing_budgets.jl
@@ -43,20 +43,21 @@ nothing # hide
 
 # We force the vorticity equation with stochastic excitation that is delta-correlated
 # in time and while spatially homogeneously and isotropically correlated. The forcing
-# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and
-# width $\delta k_f$, and it injects energy per unit area and per unit time equal
-# to $\varepsilon$.
+# has a spectrum with power in a ring in wavenumber space of radius ``k_f`` and
+# width ``\delta k_f``, and it injects energy per unit area and per unit time equal
+# to ``\varepsilon``.
 
-forcing_wavenumber = 14.0    # the central forcing wavenumber for a spectrum that is a ring in wavenumber space
-forcing_bandwidth  = 1.5     # the width of the forcing spectrum
-ε = 0.1                      # energy input rate by the forcing
+forcing_wavenumber = 14.0 * 2π/L   # the central forcing wavenumber for a spectrum that is a ring in wavenumber space
+forcing_bandwidth  = 1.5  * 2π/L   # the width of the forcing spectrum
+ε = 0.1                            # energy input rate by the forcing
 
 grid = TwoDGrid(dev, n, L)
 
 K = @. sqrt(grid.Krsq)
+
 forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2))
-forcing_spectrum[K .< ( 2 * 2π/L)] .= 0 # no power at low wavenumbers
-forcing_spectrum[K .> (20 * 2π/L)] .= 0 # no power at high wavenumbers
+@. forcing_spectrum = ifelse(K < 2  * 2π/L, 0, forcing_spectrum)      # no power at low wavenumbers
+@. forcing_spectrum = ifelse(K > 20 * 2π/L, 0, forcing_spectrum)      # no power at high wavenumbers
 ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly)
 @. forcing_spectrum *= ε/ε0             # normalize forcing to inject energy at rate ε
 
@@ -68,7 +69,8 @@ function calcF!(Fh, sol, t, clock, vars, params, grid)
   ξ = ArrayType(dev)(exp.(2π * im * rand(eltype(grid), size(sol))) / sqrt(clock.dt))
   ξ[1, 1] = 0
   @. Fh = ξ * sqrt(forcing_spectrum)
-  nothing
+  
+  return nothing
 end
 nothing # hide
 
@@ -92,7 +94,7 @@ nothing # hide
 # go back to physical space.
 calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid)
 
-heatmap(x, y, irfft(vars.Fh, grid.nx),
+heatmap(x, y, irfft(vars.Fh, grid.nx)',
      aspectratio = 1,
                c = :balance,
             clim = (-200, 200),
@@ -132,8 +134,8 @@ nothing # hide
 # We define a function that plots the vorticity field and the evolution of
 # the diagnostics: energy, enstrophy, and all terms involved in the energy and 
 # enstrophy budgets. Last, we also check (by plotting) whether the energy and enstrophy
-# budgets are accurately computed, e.g., $\mathrm{d}E/\mathrm{d}t = W^\varepsilon - 
-# R^\varepsilon - D^\varepsilon$.
+# budgets are accurately computed, e.g., ``\mathrm{d}E/\mathrm{d}t = W^\varepsilon - 
+# R^\varepsilon - D^\varepsilon``.
 
 function computetendencies_and_makeplot(prob, diags)
   sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
@@ -155,7 +157,7 @@ function computetendencies_and_makeplot(prob, diags)
 
   εᶻ = parsevalsum(forcing_spectrum / 2, grid) / (grid.Lx * grid.Ly)
 
-  pzeta = heatmap(x, y, vars.zeta,
+  pzeta = heatmap(x, y, vars.zeta',
             aspectratio = 1,
             legend = false,
                  c = :viridis,
@@ -166,7 +168,7 @@ function computetendencies_and_makeplot(prob, diags)
             yticks = -3:3,
             xlabel = "μt",
             ylabel = "y",
-             title = "∇²ψ(x, y, μt="*@sprintf("%.2f", μ*clock.t)*")",
+             title = "∇²ψ(x, y, μt=" * @sprintf("%.2f", μ * clock.t) * ")",
         framestyle = :box)
 
   pζ = plot(pzeta, size = (400, 400))