diff --git a/.travis.yml b/.travis.yml
index bba66ee6..058ed2f0 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -12,9 +12,9 @@ matrix:
   allow_failures:
     - julia: nightly
 
-after_success:
-  - julia --project -e 'using Pkg; Pkg.add("Coverage"); import GeophysicalFlows; joinpath(dirname(pathof(GeophysicalFlows)), ".."); using Coverage; Coveralls.submit(Coveralls.process_folder());'
-  - julia --project -e 'using Pkg; Pkg.add("Coverage"); import GeophysicalFlows; joinpath(dirname(pathof(GeophysicalFlows)), ".."); using Coverage; Codecov.submit(Codecov.process_folder())'
+# after_success:
+#   - julia --project -e 'using Pkg; Pkg.add("Coverage"); import GeophysicalFlows; joinpath(dirname(pathof(GeophysicalFlows)), ".."); using Coverage; Coveralls.submit(Coveralls.process_folder());'
+#   - julia --project -e 'using Pkg; Pkg.add("Coverage"); import GeophysicalFlows; joinpath(dirname(pathof(GeophysicalFlows)), ".."); using Coverage; Codecov.submit(Codecov.process_folder())'
 
 jobs:
   include:
diff --git a/Project.toml b/Project.toml
index 9db76377..f8726929 100644
--- a/Project.toml
+++ b/Project.toml
@@ -22,8 +22,8 @@ JLD2 = "≥ 0.1.2"
 julia = "≥ 1.0.0"
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Coverage = "a2441757-f6aa-5fb2-8edb-039e3f45d037"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test", "Coverage"]
diff --git a/src/barotropicqg.jl b/src/barotropicqg.jl
index 0acfa323..b46dd400 100644
--- a/src/barotropicqg.jl
+++ b/src/barotropicqg.jl
@@ -22,13 +22,6 @@ using LinearAlgebra: mul!, ldiv!
 using FourierFlows: getfieldspecs, structvarsexpr, parsevalsum, parsevalsum2
 
 abstract type BarotropicQGVars <: AbstractVars end
-abstract type BarotropicQGForcedVars <: BarotropicQGVars end
-
-const realscalarvars = [:U]
-const physicalvars = [:q, :zeta, :psi, :u, :v]
-const transformvars = [ Symbol(var, :h) for var in physicalvars ]
-const forcedvars = [:Fqh]
-const stochforcedvars = [:prevsol]
 
 nothingfunction(args...) = nothing
 
@@ -40,57 +33,44 @@ Construct a BarotropicQG turbulence problem.
 
 function Problem(;
     # Numerical parameters
-            nx = 256,
-            Lx = 2π,
-            ny = nx,
-            Ly = Lx,
-            dt = 0.01,
-    # Physical parameters
-            f0 = 1.0,
-          beta = 0.0,
-           eta = nothing,
-    # Drag and/or hyper-/hypo-viscosity
-            nu = 0.0,
-           nnu = 1,
-            mu = 0.0,
-   # Timestepper and eqn options
-       stepper = "RK4",
-        calcFU = nothingfunction,
-        calcFq = nothingfunction,
-    stochastic = false,
-             T = Float64)
+          nx = 256,
+          Lx = 2π,
+          ny = nx,
+          Ly = Lx,
+          dt = 0.01,
+  # Physical parameters
+          f0 = 1.0,
+           β = 0.0,
+         eta = nothing,
+  # Drag and/or hyper-/hypo-viscosity
+           ν = 0.0,
+          nν = 1,
+           μ = 0.0,
+  # Timestepper and eqn options
+     stepper = "RK4",
+      calcFU = nothingfunction,
+      calcFq = nothingfunction,
+  stochastic = false,
+           T = Float64,
+         dev = CPU())
 
   # the grid
-  gr  = TwoDGrid(nx, Lx, ny, Ly; T=T)
+    gr = TwoDGrid(dev, nx, Lx, ny, Ly; T=T)
   x, y = gridpoints(gr)
 
   # topographic PV
   if eta==nothing
-    eta = 0*x
-    etah = rfft(eta)
+     eta = zeros(dev, T, (nx, ny))
   end
 
-  if typeof(eta)!=Array{T, 2} #this is true if eta was passes in Problem as a function
-    pr = Params(gr, f0, beta, eta, mu, nu, nnu, calcFU, calcFq)
-  else
-    pr = Params{T}(f0, beta, eta, rfft(eta), mu, nu, nnu, calcFU, calcFq)
-  end
+  pr = !(typeof(eta)<:ArrayType(dev)) ? Params(gr, f0, β, eta, μ, ν, nν, calcFU, calcFq) : Params(f0, β, eta, rfft(eta), μ, ν, nν, calcFU, calcFq)
 
-  if calcFq == nothingfunction && calcFU == nothingfunction
-      vs = Vars(gr)
-  else
-      if stochastic ==  false
-          vs = ForcedVars(gr)
-      elseif stochastic == true
-          vs = StochasticForcedVars(gr)
-      end
-  end
+  vs = (calcFq == nothingfunction && calcFU == nothingfunction) ? Vars(dev, gr) : (stochastic ? StochasticForcedVars(dev, gr) : ForcedVars(dev, gr))
+  
   eq = Equation(pr, gr)
-  FourierFlows.Problem(eq, stepper, dt, gr, vs, pr)
+  FourierFlows.Problem(eq, stepper, dt, gr, vs, pr, dev)
 end
 
-InitialValueProblem(; kwargs...) = Problem(; kwargs...)
-ForcedProblem(; kwargs...) = Problem(; kwargs...)
 
 
 # ----------
@@ -98,33 +78,32 @@ ForcedProblem(; kwargs...) = Problem(; kwargs...)
 # ----------
 
 """
-    Params(g::TwoDGrid, f0, beta, FU, eta, mu, nu, nnu, calcFU, calcFq)
+    Params(g::TwoDGrid, f0, β, FU, eta, μ, ν, nν, calcFU, calcFq)
 
 Returns the params for an unforced two-dimensional barotropic QG problem.
 """
-struct Params{T} <: AbstractParams
-  f0::T                      # Constant planetary vorticity
-  beta::T                    # Planetary vorticity y-gradient
-  eta::Array{T,2}            # Topographic PV
-  etah::Array{Complex{T},2}  # FFT of Topographic PV
-  mu::T                      # Linear drag
-  nu::T                      # Viscosity coefficient
-  nnu::Int                   # Hyperviscous order (nnu=1 is plain old viscosity)
-  calcFU::Function    # Function that calculates the forcing F(t) on
-                      # domain-averaged zonal flow U(t)
-  calcFq!::Function   # Function that calculates the forcing on QGPV q
+struct Params{T, Aphys, Atrans} <: AbstractParams
+       f0 :: T            # Constant planetary vorticity
+        β :: T            # Planetary vorticity y-gradient
+      eta :: Aphys        # Topographic PV
+     etah :: Atrans       # FFT of Topographic PV
+        μ :: T            # Linear drag
+        ν :: T            # Viscosity coefficient
+       nν :: Int          # Hyperviscous order (nν=1 is plain old viscosity)
+   calcFU :: Function     # Function that calculates the forcing F(t) on
+                          # domain-averaged zonal flow U(t)
+  calcFq! :: Function     # Function that calculates the forcing on QGPV q
 end
 
 """
-    Params(g::TwoDGrid, f0, beta, eta::Function, mu, nu, nnu, calcFU, calcFq)
+    Params(g::TwoDGrid, f0, β, eta::Function, μ, ν, nν, calcFU, calcFq)
 
 Constructor for Params that accepts a generating function for the topographic PV.
 """
-function Params(g::AbstractGrid{T}, f0, beta, eta::Function, mu, nu, nnu, calcFU, calcFq)  where T
-  x, y = gridpoints(g)
-  etagrid = eta(x, y)
-  etah = rfft(etagrid)
-  Params{T}(f0, beta, etagrid, etah, mu, nu, nnu, calcFU, calcFq)
+function Params(g::AbstractGrid{T, A}, f0, β, eta::Function, μ, ν, nν::Int, calcFU, calcFq) where {T, A}
+  etagrid = A([eta(g.x[i], g.y[j]) for i=1:g.nx, j=1:g.ny])
+     etah = rfft(etagrid)
+  Params(f0, β, etagrid, etah, μ, ν, nν, calcFU, calcFq)
 end
 
 
@@ -137,8 +116,8 @@ end
 
 Returns the equation for two-dimensional barotropic QG problem with params p and grid g.
 """
-function Equation(p::Params, g::AbstractGrid{T}) where T
-  L = @. -p.mu - p.nu*g.Krsq^p.nnu + im*p.beta*g.kr*g.invKrsq
+function Equation(p::Params, g::AbstractGrid)
+  L = @. - p.μ - p.ν*g.Krsq^p.nν + im*p.β*g.kr*g.invKrsq
   L[1, 1] = 0
   FourierFlows.Equation(L, calcN!, g)
 end
@@ -148,52 +127,59 @@ end
 # Vars
 # ----
 
-varspecs = cat(
-  getfieldspecs(realscalarvars, :(Array{T,0})),
-  getfieldspecs(physicalvars, :(Array{T,2})),
-  getfieldspecs(transformvars, :(Array{Complex{T},2})),
-  dims=1)
-
-forcedvarspecs = cat(varspecs, getfieldspecs(forcedvars, :(Array{Complex{T},2})), dims=1)
-stochforcedvarspecs = cat(forcedvarspecs, getfieldspecs(stochforcedvars, :(Array{Complex{T},2})), dims=1)
+struct Vars{Ascalar, Aphys, Atrans, F, P} <: BarotropicQGVars
+        U :: Ascalar
+        q :: Aphys
+     zeta :: Aphys
+      psi :: Aphys
+        u :: Aphys
+        v :: Aphys
+       qh :: Atrans
+    zetah :: Atrans
+     psih :: Atrans
+       uh :: Atrans
+       vh :: Atrans
+      Fqh :: F
+  prevsol :: P
+end
 
-# Construct Vars types
-eval(structvarsexpr(:Vars, varspecs; parent=:BarotropicQGVars))
-eval(structvarsexpr(:ForcedVars, forcedvarspecs; parent=:BarotropicQGForcedVars))
-eval(structvarsexpr(:StochasticForcedVars, stochforcedvarspecs; parent=:BarotropicQGForcedVars))
+const ForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, <:AbstractArray, Nothing}
+const StochasticForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, <:AbstractArray, <:AbstractArray}
 
 """
-    Vars(g)
+    Vars(dev, g)
 
-Returns the vars for unforced two-dimensional barotropic QG problem with grid g.
+Returns the vars for unforced two-dimensional barotropic QG problem on device dev and with grid g.
 """
-function Vars(g::AbstractGrid{T}) where T
-  U = Array{T,0}(undef, ); U[] = 0
-  @createarrays T (g.nx, g.ny) q u v psi zeta
-  @createarrays Complex{T} (g.nkr, g.nl) qh uh vh psih zetah
-  Vars(U, q, zeta, psi, u, v, qh, zetah, psih, uh, vh)
+function Vars(dev::Dev, g::AbstractGrid{T}) where {Dev, T}
+  U = ArrayType(dev, T, 0)(undef, ); U[] = 0
+  @devzeros Dev T (g.nx, g.ny) q u v psi zeta
+  @devzeros Dev Complex{T} (g.nkr, g.nl) qh uh vh psih zetah
+  Vars(U, q, zeta, psi, u, v, qh, zetah, psih, uh, vh, nothing, nothing)
 end
 
 """
-    ForcedVars(g)
+    ForcedVars(dev, g)
 
-Returns the vars for forced two-dimensional barotropic QG problem with grid g.
+Returns the vars for forced two-dimensional barotropic QG problem on device dev and with grid g.
 """
-function ForcedVars(g::AbstractGrid{T}) where T
-  v = Vars(g)
-  Fqh = zeros(Complex{T}, (g.nkr, g.nl))
-  ForcedVars(getfield.(Ref(v), fieldnames(typeof(v)))..., Fqh)
+function ForcedVars(dev::Dev, g::AbstractGrid{T}) where {Dev, T}
+  U = ArrayType(dev, T, 0)(undef, ); U[] = 0
+  @devzeros Dev T (g.nx, g.ny) q u v psi zeta
+  @devzeros Dev Complex{T} (g.nkr, g.nl) qh uh vh psih zetah Fqh
+  Vars(U, q, zeta, psi, u, v, qh, zetah, psih, uh, vh, Fqh, nothing)
 end
 
 """
-    StochasticForcedVars(g)
+    StochasticForcedVars(dev, g)
 
-Returns the vars for stochastically forced two-dimensional barotropic QG problem with grid g.
+Returns the vars for stochastically forced two-dimensional barotropic QG problem on device dev and with grid g.
 """
-function StochasticForcedVars(g::AbstractGrid{T}) where T
-  v = ForcedVars(g)
-  prevsol = zeros(Complex{T}, (g.nkr, g.nl))
-  StochasticForcedVars(getfield.(Ref(v), fieldnames(typeof(v)))..., prevsol)
+function StochasticForcedVars(dev::Dev, g::AbstractGrid{T}) where {Dev, T}
+  U = ArrayType(dev, T, 0)(undef, ); U[] = 0
+  @devzeros Dev T (g.nx, g.ny) q u v psi zeta
+  @devzeros Dev Complex{T} (g.nkr, g.nl) qh uh vh psih zetah Fqh prevsol
+  Vars(U, q, zeta, psi, u, v, qh, zetah, psih, uh, vh, Fqh, prevsol)
 end
 
 
@@ -301,7 +287,7 @@ function set_zeta!(sol, v::Vars, p, g, zeta)
   nothing
 end
 
-function set_zeta!(sol, v::BarotropicQGForcedVars, p, g, zeta)
+function set_zeta!(sol, v::Union{ForcedVars, StochasticForcedVars}, p, g, zeta)
   v.U[] = deepcopy(sol[1, 1])
   mul!(v.zetah, g.rfftplan, zeta)
   v.zetah[1, 1] = 0.0
@@ -357,18 +343,18 @@ meanenergy(prob) = real(0.5*prob.sol[1, 1].^2)
 """
 Returns the enstrophy of the domain-averaged U.
 """
-meanenstrophy(prob) = real(prob.params.beta*prob.sol[1, 1])
+meanenstrophy(prob) = real(prob.params.β*prob.sol[1, 1])
 
 """
     dissipation(prob)
     dissipation(s, v, p, g)
 
-Returns the domain-averaged dissipation rate. nnu must be >= 1.
+Returns the domain-averaged dissipation rate. nν must be >= 1.
 """
 @inline function dissipation(sol, v, p, g)
-  @. v.uh = g.Krsq^(p.nnu-1) * abs2(sol)
+  @. v.uh = g.Krsq^(p.nν-1) * abs2(sol)
   v.uh[1, 1] = 0
-  p.nu/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
+  p.ν/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
 end
 
 @inline dissipation(prob) = dissipation(prob.sol, prob.vars, prob.params, prob.grid)
@@ -396,12 +382,12 @@ end
     drag(prob)
     drag(s, v, p, g)
 
-Returns the extraction of domain-averaged energy by drag mu.
+Returns the extraction of domain-averaged energy by drag μ.
 """
 @inline function drag(sol, v, p, g)
   @. v.uh = g.Krsq^(-1) * abs2(sol)
   v.uh[1, 1] = 0
-  p.mu/(g.Lx*g.Ly)*FourierFlows.parsevalsum(v.uh, g)
+  p.μ/(g.Lx*g.Ly)*FourierFlows.parsevalsum(v.uh, g)
 end
 
 @inline drag(prob) = drag(prob.sol, prob.vars, prob.params, prob.grid)
diff --git a/test/runtests.jl b/test/runtests.jl
index 8754aa85..bca2ce33 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -22,6 +22,7 @@ devices = (CPU(),)
 const rtol_lambdipole = 1e-2 # tolerance for lamb dipole tests
 const rtol_multilayerqg = 1e-13 # tolerance for multilayerqg forcing tests
 const rtol_twodturb = 1e-13 # tolerance for twodturb forcing tests
+const rtol_barotropicQG = 1e-13 # tolerance for twodturb forcing tests
 
 "Get the CFL number, assuming a uniform grid with `dx=dy`."
 cfl(u, v, dt, dx) = maximum([maximum(abs.(u)), maximum(abs.(v))]*dt/dx)
@@ -54,32 +55,31 @@ for dev in devices
     @test TwoDTurb.nothingfunction() == nothing
   end
 
+  @testset "BarotropicQG" begin
+    include("test_barotropicqg.jl")
+
+    @test test_bqg_rossbywave("ETDRK4", 1e-2, 20, dev)
+    @test test_bqg_rossbywave("FilteredETDRK4", 1e-2, 20, dev)
+    @test test_bqg_rossbywave("RK4", 1e-2, 20, dev)
+    @test test_bqg_rossbywave("FilteredRK4", 1e-2, 20, dev)
+    @test test_bqg_rossbywave("AB3", 1e-3, 200, dev)
+    @test test_bqg_rossbywave("FilteredAB3", 1e-3, 200, dev)
+    @test test_bqg_rossbywave("ForwardEuler", 1e-4, 2000, dev)
+    @test test_bqg_rossbywave("FilteredForwardEuler", 1e-4, 2000, dev)
+    @test test_bqg_stochasticforcingbudgets(dev)
+    @test test_bqg_deterministicforcingbudgets(dev)
+    @test test_bqg_advection(0.0005, "ForwardEuler", dev)
+    @test test_bqg_formstress(0.01, "ForwardEuler", dev)
+    @test test_bqg_energyenstrophy(dev)
+    @test test_bqg_meanenergyenstrophy(dev)
+    @test BarotropicQG.nothingfunction() == nothing
+  end
+
 end
 
 
 println("rest of tests only on CPU")
 
-@testset "BarotropicQG" begin
-  include("test_barotropicqg.jl")
-
-  @test test_bqg_rossbywave("ETDRK4", 1e-2, 20)
-  @test test_bqg_rossbywave("FilteredETDRK4", 1e-2, 20)
-  @test test_bqg_rossbywave("RK4", 1e-2, 20)
-  @test test_bqg_rossbywave("FilteredRK4", 1e-2, 20)
-  @test test_bqg_rossbywave("AB3", 1e-3, 200)
-  @test test_bqg_rossbywave("FilteredAB3", 1e-3, 200)
-  @test test_bqg_rossbywave("ForwardEuler", 1e-4, 2000)
-  @test test_bqg_rossbywave("FilteredForwardEuler", 1e-4, 2000)
-  @test test_bqg_stochasticforcingbudgets()
-  @test test_bqg_deterministicforcingbudgets()
-  @test test_bqg_advection(0.0005, "ForwardEuler")
-  @test test_bqg_formstress(0.01, "ForwardEuler")
-  @test test_bqg_energyenstrophy()
-  @test test_bqg_meanenergyenstrophy()
-  @test test_bqg_problemtype(Float32)
-  @test BarotropicQG.nothingfunction() == nothing
-end
-
 @testset "BarotropicQGQL" begin
   include("test_barotropicqgql.jl")
 
diff --git a/test/test_barotropicqg.jl b/test/test_barotropicqg.jl
index ce488ce5..d3ec7a39 100644
--- a/test/test_barotropicqg.jl
+++ b/test/test_barotropicqg.jl
@@ -3,25 +3,23 @@
 
 Evolves a Rossby wave and compares with the analytic solution.
 """
-function test_bqg_rossbywave(stepper, dt, nsteps)
+function test_bqg_rossbywave(stepper, dt, nsteps, dev::Device=CPU())
     nx = 64
-  beta = 2.0
     Lx = 2π
-    mu = 0.0
-    nu = 0.0
-
-  gr  = TwoDGrid(nx, Lx)
-  x, y = gridpoints(gr)
+     β = 2.0
+     μ = 0.0
+     ν = 0.0
+     T = Float64
 
   # the following if statement is called so that all the cases of
   # Problem() fuction are tested
   if stepper=="ForwardEuler"
-    eta = zeros(nx, nx)
+    eta = zeros(dev, T, (nx, nx))
   else
     eta(x, y) = 0*x
   end
 
-  prob = BarotropicQG.InitialValueProblem(nx=nx, Lx=Lx, eta=eta, beta=beta, mu=mu, nu=nu, stepper=stepper, dt=dt)
+  prob = BarotropicQG.Problem(nx=nx, Lx=Lx, eta=eta, β=β, μ=μ, ν=ν, stepper=stepper, dt=dt, dev=dev)
   sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
 
   x, y = gridpoints(g)
@@ -30,7 +28,7 @@ function test_bqg_rossbywave(stepper, dt, nsteps)
    ampl = 1e-2
   kwave = 3.0*2π/g.Lx
   lwave = 2.0*2π/g.Ly
-      ω = -p.beta*kwave/(kwave^2.0 + lwave^2.0)
+      ω = -p.β*kwave/(kwave^2.0 + lwave^2.0)
      ζ0 = @. ampl*cos(kwave*x)*cos(lwave*y)
     ζ0h = rfft(ζ0)
 
@@ -50,24 +48,23 @@ end
 
 Tests if the energy budgets are closed for BarotropicQG with stochastic forcing.
 """
-function test_bqg_stochasticforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7, nnu=2, mu=1e-1, message=false)
+function test_bqg_stochasticforcingbudgets(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1, T=Float64)
   n, L  = 256, 2π
-  nu, nnu = 1e-7, 2
-  mu = 1e-1
-  dt, tf = 0.005, 0.1/mu
+  ν, nν = 1e-7, 2
+  μ = 1e-1
+  dt, tf = 0.005, 0.1/μ
   nt = round(Int, tf/dt)
-  ns = 1
-
+   
   # Forcing
   kf, dkf = 12.0, 2.0
   ε = 0.1
 
-  gr  = TwoDGrid(n, L)
+  gr = TwoDGrid(dev, n, L)
   x, y = gridpoints(gr)
 
-  Kr = [ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]
+  Kr = ArrayType(dev)([ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl])
 
-  force2k = zero(gr.Krsq)
+  force2k = zeros(dev, T, (gr.nkr, gr.nl))
   @. force2k = exp.(-(sqrt(gr.Krsq)-kf)^2/(2*dkf^2))
   @. force2k[gr.Krsq .< 2.0^2 ] = 0
   @. force2k[gr.Krsq .> 20.0^2 ] = 0
@@ -78,14 +75,14 @@ function test_bqg_stochasticforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7, nnu
   Random.seed!(1234)
 
   function calcFq!(Fqh, sol, t, cl, v, p, g)
-    eta = exp.(2π*im*rand(Float64, size(sol)))/sqrt(cl.dt)
+    eta = ArrayType(dev)(exp.(2π*im*rand(T, size(sol)))/sqrt(cl.dt))
     eta[1, 1] = 0
     @. Fqh = eta * sqrt(force2k)
     nothing
   end
 
-  prob = BarotropicQG.ForcedProblem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, dt=dt,
-   stepper="RK4", calcFq=calcFq!, stochastic=true)
+  prob = BarotropicQG.Problem(nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt,
+   stepper="RK4", calcFq=calcFq!, stochastic=true, dev=dev)
 
   sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
 
@@ -102,11 +99,9 @@ function test_bqg_stochasticforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7, nnu
 
   BarotropicQG.updatevars!(prob)
 
-  cfl = cl.dt*maximum([maximum(v.v)/g.dx, maximum(v.u)/g.dy])
-
   E, D, W, R = diags
 
-  t = round(mu*cl.t, digits=2)
+  t = round(μ*cl.t, digits=2)
 
   i₀ = 1
   dEdt = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt
@@ -121,10 +116,6 @@ function test_bqg_stochasticforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7, nnu
 
   residual = dEdt - total
 
-  if message
-    println("step: %04d, t: %.1f, cfl: %.3f, time: %.2f s\n", cl.step, cl.t, cfl, t)
-  end
-  # println(mean(abs.(residual)))
   isapprox(mean(abs.(residual)), 0, atol=1e-4)
 end
 
@@ -133,29 +124,27 @@ end
 
 Tests if the energy budgets are closed for BarotropicQG with stochastic forcing.
 """
-function test_bqg_deterministicforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7, nnu=2, mu=1e-1, message=false)
+function test_bqg_deterministicforcingbudgets(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1)
   n, L  = 256, 2π
-  nu, nnu = 1e-7, 2
-  mu = 1e-1
-  dt, tf = 0.005, 0.1/mu
+  ν, nν = 1e-7, 2
+  μ = 1e-1
+  dt, tf = 0.005, 0.1/μ
   nt = round(Int, tf/dt)
-  ns = 1
-
 
-  # Forcing = 0.01cos(4x)cos(5y)cos(2t)
-  gr  = TwoDGrid(n, L)
+  gr  = TwoDGrid(dev, n, L)
   x, y = gridpoints(gr)
   k0, l0 = gr.kr[2], gr.l[2]
 
-  f = @. 0.01*cos(4*k0*x)*cos(5*l0*y)
+  # Forcing = 0.01cos(4x)cos(5y)cos(2t)
+  f = @. 0.01*cos(4k0*x)*cos(5l0*y)
   fh = rfft(f)
   function calcFq!(Fqh, sol, t, cl, v, p, g)
-    @. Fqh = fh*cos(2*t)
+    Fqh = fh*cos(2t)
     nothing
   end
 
-  prob = BarotropicQG.ForcedProblem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, dt=dt,
-   stepper="RK4", calcFq=calcFq!, stochastic=false)
+  prob = BarotropicQG.Problem(nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt,
+   stepper="RK4", calcFq=calcFq!, stochastic=false, dev=dev)
 
   sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
 
@@ -176,7 +165,7 @@ function test_bqg_deterministicforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7,
 
   E, D, W, R = diags
 
-  t = round(mu*cl.t, digits=2)
+  t = round(μ*cl.t, digits=2)
 
   i₀ = 1
   dEdt = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt
@@ -188,10 +177,6 @@ function test_bqg_deterministicforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7,
 
   residual = dEdt - total
 
-  if message
-    println("step: %04d, t: %.1f, cfl: %.3f, time: %.2f s\n", prob.step, cl.t, cfl, tc)
-  end
-  # println(mean(abs.(residual)))
   isapprox(mean(abs.(residual)), 0, atol=1e-8)
 end
 
@@ -202,25 +187,25 @@ Tests the advection term in the twodturb module by timestepping a
 test problem with timestep dt and timestepper identified by the string stepper.
 The test problem is derived by picking a solution ζf (with associated
 streamfunction ψf) for which the advection term J(ψf, ζf) is non-zero. Next, a
-forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - nuΔζf. One solution
+forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - ν∇²ζf. One solution
 to the vorticity equation forced by this Ff is then ζf. (This solution may not
 be realized, at least at long times, if it is unstable.)
 """
-function test_bqg_advection(dt, stepper; n=128, L=2π, nu=1e-2, nnu=1, mu=0.0, message=false)
+function test_bqg_advection(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0)
   n, L  = 128, 2π
-  nu, nnu = 1e-2, 1
-  mu = 0.0
+  ν, nν = 1e-2, 1
+  μ = 0.0
   tf = 1.0
   nt = round(Int, tf/dt)
 
-  gr  = TwoDGrid(n, L)
+  gr  = TwoDGrid(dev, n, L)
   x, y = gridpoints(gr)
 
   psif = @. sin(2x)*cos(2y) + 2sin(x)*cos(3y)
-  qf = @. -8sin(2x)*cos(2y) - 20sin(x)*cos(3y)
+    qf = @. -8sin(2x)*cos(2y) - 20sin(x)*cos(3y)
 
   Ff = @. -(
-    nu*( 64sin(2x)*cos(2y) + 200sin(x)*cos(3y) )
+    ν*( 64sin(2x)*cos(2y) + 200sin(x)*cos(3y) )
     + 8*( cos(x)*cos(3y)*sin(2x)*sin(2y) - 3cos(2x)*cos(2y)*sin(x)*sin(3y) )
   )
 
@@ -232,14 +217,15 @@ function test_bqg_advection(dt, stepper; n=128, L=2π, nu=1e-2, nnu=1, mu=0.0, m
     nothing
   end
 
-  prob = BarotropicQG.ForcedProblem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, dt=dt, stepper=stepper, calcFq=calcFq!)
+  prob = BarotropicQG.Problem(nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper, calcFq=calcFq!, dev=dev)
   sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
   BarotropicQG.set_zeta!(prob, qf)
 
   # Step forward
   stepforward!(prob, round(Int, nt))
   BarotropicQG.updatevars!(prob)
-  isapprox(v.q, qf, rtol=1e-13)
+  
+  isapprox(v.q, qf, rtol=rtol_barotropicQG)
 end
 
 """
@@ -247,35 +233,35 @@ end
 
 Tests the form stress term that forces the domain-averaged zonal flow U(t).
 """
-function test_bqg_formstress(dt, stepper; n=128, L=2π, nu=0.0, nnu=1, mu=0.0, message=false)
+function test_bqg_formstress(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=0.0, nν=1, μ=0.0)
   n, L  = 128, 2π
-  nu, nnu = 1e-2, 1
-  mu = 0.0
+  ν, nν = 1e-2, 1
+  μ = 0.0
   tf = 1
   nt = 1
 
-  gr  = TwoDGrid(n, L)
+  gr  = TwoDGrid(dev, n, L)
   x, y = gridpoints(gr)
 
-  zetai = @. -20*sin(10*x)*cos(10*y)
+  zetai = @. -20*sin(10x)*cos(10y)
   topoPV(x, y) = @. cos(10x)*cos(10y)
   F(t) = 0 #no forcing
 
   answer = 0.25 # this is what <v*eta> should be
 
-  prob = BarotropicQG.ForcedProblem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, dt=dt, stepper=stepper, eta=topoPV, calcFU = F)
+  prob = BarotropicQG.Problem(nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper, eta=topoPV, calcFU=F, dev=dev)
   BarotropicQG.set_zeta!(prob, zetai)
   BarotropicQG.updatevars!(prob)
 
   # Step forward
   stepforward!(prob, nt)
-  isapprox(prob.timestepper.N[1, 1], answer, rtol=1e-13)
+  isapprox(prob.timestepper.N[1, 1], answer, rtol=rtol_barotropicQG)
 end
 
-function test_bqg_energyenstrophy()
+function test_bqg_energyenstrophy(dev::Device=CPU())
   nx, Lx  = 64, 2π
   ny, Ly  = 64, 3π
-  g  = TwoDGrid(nx, Lx, ny, Ly)
+  g  = TwoDGrid(dev, nx, Lx, ny, Ly)
   k0, l0 = g.k[2], g.l[2] # fundamental wavenumbers
   x, y = gridpoints(g)
 
@@ -286,7 +272,7 @@ function test_bqg_energyenstrophy()
    psi0 = @. sin(2k0*x)*cos(2l0*y) + 2sin(k0*x)*cos(3l0*y)
   zeta0 = @. -((2k0)^2+(2l0)^2)*sin(2k0*x)*cos(2l0*y) - (k0^2+(3l0)^2)*2sin(k0*x)*cos(3l0*y)
 
-  prob = BarotropicQG.InitialValueProblem(nx=nx, Lx=Lx, ny=ny, Ly=Ly, eta = eta, stepper="ForwardEuler")
+  prob = BarotropicQG.Problem(nx=nx, Lx=Lx, ny=ny, Ly=Ly, eta = eta, stepper="ForwardEuler", dev=dev)
   sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
   BarotropicQG.set_zeta!(prob, zeta0)
   BarotropicQG.updatevars!(prob)
@@ -294,14 +280,14 @@ function test_bqg_energyenstrophy()
   energyzeta0 = BarotropicQG.energy(prob)
   enstrophyzeta0 = BarotropicQG.enstrophy(prob)
 
-  isapprox(energyzeta0, energy_calc, rtol=1e-13) && isapprox(enstrophyzeta0, enstrophy_calc, rtol=1e-13) &&
+  isapprox(energyzeta0, energy_calc, rtol=rtol_barotropicQG) && isapprox(enstrophyzeta0, enstrophy_calc, rtol=rtol_barotropicQG) &&
   BarotropicQG.addforcing!(prob.timestepper.N, sol, cl.t, cl, v, p, g)==nothing
 end
 
-function test_bqg_meanenergyenstrophy()
+function test_bqg_meanenergyenstrophy(dev::Device=CPU())
   nx, Lx  = 64, 2π
   ny, Ly  = 96, 3π
-  g  = TwoDGrid(nx, Lx, ny, Ly)
+  g = TwoDGrid(dev, nx, Lx, ny, Ly)
   k0, l0 = g.k[2], g.l[2] # fundamental wavenumbers
   x, y = gridpoints(g)
 
@@ -309,14 +295,14 @@ function test_bqg_meanenergyenstrophy()
   eta(x, y) = @. cos(10x)*cos(10y)
   psi0 = @. sin(2k0*x)*cos(2l0*y) + 2sin(k0*x)*cos(3l0*y)
  zeta0 = @. -((2k0)^2+(2l0)^2)*sin(2k0*x)*cos(2l0*y) - (k0^2+(3l0)^2)*2sin(k0*x)*cos(3l0*y)
-  beta = 10.0
+  β = 10.0
   U = 1.2
 
   energy_calc = 29/9
   enstrophy_calc = 2701/162
 
-  prob = BarotropicQG.ForcedProblem(nx=nx, Lx=Lx, ny=ny, Ly=Ly, beta=beta, eta=eta, calcFU = calcFU,
-                                    stepper="ForwardEuler")
+  prob = BarotropicQG.Problem(nx=nx, Lx=Lx, ny=ny, Ly=Ly, β=β, eta=eta, calcFU = calcFU,
+                                    stepper="ForwardEuler", dev=dev)
 
   BarotropicQG.set_zeta!(prob, zeta0)
   BarotropicQG.set_U!(prob, U)
@@ -328,15 +314,9 @@ function test_bqg_meanenergyenstrophy()
   energyzeta0 = BarotropicQG.energy(prob)
   enstrophyzeta0 = BarotropicQG.enstrophy(prob)
 
-  (isapprox(energyU, 0.5*U^2, rtol=1e-13) &&
-    isapprox(enstrophyU, beta*U, rtol=1e-13) &&
-    isapprox(energyzeta0, energy_calc, rtol=1e-13) &&
-    isapprox(enstrophyzeta0, enstrophy_calc, rtol=1e-13)
+  (isapprox(energyU, 0.5*U^2, rtol=rtol_barotropicQG) &&
+    isapprox(enstrophyU, β*U, rtol=rtol_barotropicQG) &&
+    isapprox(energyzeta0, energy_calc, rtol=rtol_barotropicQG) &&
+    isapprox(enstrophyzeta0, enstrophy_calc, rtol=rtol_barotropicQG)
   )
 end
-
-function test_bqg_problemtype(T=Float32)
-  prob = BarotropicQG.Problem(T=T)
-
-  (typeof(prob.sol)==Array{Complex{T},2} && typeof(prob.grid.Lx)==T && typeof(prob.grid.x)==Array{T,2} && typeof(prob.vars.u)==Array{T,2})
-end
diff --git a/test/test_twodturb.jl b/test/test_twodturb.jl
index 5be8e293..56d440ea 100644
--- a/test/test_twodturb.jl
+++ b/test/test_twodturb.jl
@@ -140,7 +140,7 @@ Tests the advection term in the twodturb module by timestepping a
 test problem with timestep dt and timestepper identified by the string stepper.
 The test problem is derived by picking a solution ζf (with associated
 streamfunction ψf) for which the advection term J(ψf, ζf) is non-zero. Next, a
-forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - νΔζf. One solution
+forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - ν∇²ζf. One solution
 to the vorticity equation forced by this Ff is then ζf. (This solution may not
 be realized, at least at long times, if it is unstable.)
 """
@@ -151,7 +151,7 @@ function test_twodturb_advection(dt, stepper, dev::Device=CPU(); n=128, L=2π, 
   tf = 1.0
   nt = round(Int, tf/dt)
 
-  gr   = TwoDGrid(dev, n, L)
+  gr = TwoDGrid(dev, n, L)
   x, y = gridpoints(gr)
 
    psif = @.   sin(2x)*cos(2y) +  2sin(x)*cos(3y)