diff --git a/test/runtests.jl b/test/runtests.jl
index 544a71ad..f331ecda 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -101,7 +101,7 @@ for dev in devices
   
   @testset "SurfaceQG" begin
     include("test_surfaceqg.jl")
-    
+      
     @test test_sqg_kineticenergy_buoyancyvariance(dev)
     @test test_sqg_advection(0.0005, "ForwardEuler", dev)
     @test test_sqg_deterministicforcing_buoyancy_variance_budget(dev)
@@ -113,6 +113,26 @@ for dev in devices
     @test SurfaceQG.nothingfunction() == nothing
   end
   
+  if dev == CPU()
+    @testset "MultilayerQG" begin
+      include("test_multilayerqg.jl")
+      
+      @test test_pvtofromstreamfunction_2layer(dev)
+      @test test_pvtofromstreamfunction_3layer(dev)
+      @test test_mqg_rossbywave("RK4", 1e-2, 20, dev)
+      @test test_mqg_nonlinearadvection(0.005, "ForwardEuler", dev)
+      @test test_mqg_linearadvection(0.005, "ForwardEuler", dev)
+      @test test_mqg_energies(dev)
+      @test test_mqg_energysinglelayer(dev)
+      @test test_mqg_fluxes(dev)
+      @test test_mqg_fluxessinglelayer(dev)
+      @test test_mqg_setqsetψ(dev)
+      @test test_mqg_paramsconstructor(dev)
+      @test test_mqg_stochasticforcedproblemconstructor(dev)
+      @test test_mqg_problemtype(dev, Float32)
+      @test MultilayerQG.nothingfunction() == nothing
+    end
+  end
 end
 
 dev = CPU()
diff --git a/test/test_multilayerqg.jl b/test/test_multilayerqg.jl
index a351315c..9125778a 100644
--- a/test/test_multilayerqg.jl
+++ b/test/test_multilayerqg.jl
@@ -224,11 +224,11 @@ function test_mqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); n=128, L=2
   @views ψf[:, :, 2] = ψ2
 
   MultilayerQG.set_q!(prob, qf)
-
+  
   stepforward!(prob, nt)
   MultilayerQG.updatevars!(prob)
-
-  return isapprox(vs.q, qf, rtol=rtol_multilayerqg) && isapprox(vs.ψ, ψf, rtol=rtol_multilayerqg)
+  
+  return isapprox(vs.q, qf, rtol=rtol_multilayerqg) && isapprox(vs.ψ, ψf, rtol=rtol_multilayerqg)		
 end
 
 """
diff --git a/test/test_surfaceqg.jl b/test/test_surfaceqg.jl
index 54ae154f..4ee776be 100644
--- a/test/test_surfaceqg.jl
+++ b/test/test_surfaceqg.jl
@@ -113,9 +113,9 @@ function test_sqg_noforcing(dev::Device=CPU())
   prob_unforced = SurfaceQG.Problem(dev; nx=n, Lx=L, stepper="ForwardEuler")
 
   SurfaceQG.addforcing!(prob_unforced.timestepper.N, prob_unforced.sol, prob_unforced.clock.t, prob_unforced.clock, prob_unforced.vars, prob_unforced.params, prob_unforced.grid)
-  
+    
   function calcF!(Fh, sol, t, clock, vars, params, grid)
-    Fh .= 2*ones(size(sol))
+    Fh .= 2 * ArrayType(dev)(ones(size(sol)))
     return nothing
   end
   
@@ -123,7 +123,7 @@ function test_sqg_noforcing(dev::Device=CPU())
 
   SurfaceQG.addforcing!(prob_forced.timestepper.N, prob_forced.sol, prob_forced.clock.t, prob_forced.clock, prob_forced.vars, prob_forced.params, prob_forced.grid)
   
-  return prob_unforced.timestepper.N == Complex.(zeros(size(prob_unforced.sol))) && prob_forced.timestepper.N == Complex.(2*ones(size(prob_unforced.sol)))
+  return prob_unforced.timestepper.N == Complex.(ArrayType(dev)(zeros(size(prob_unforced.sol)))) && prob_forced.timestepper.N == Complex.(2*ArrayType(dev)(ones(size(prob_unforced.sol))))
 end
 
 function test_sqg_deterministicforcing_buoyancy_variance_budget(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, tf=10.0)