diff --git a/Project.toml b/Project.toml
index 6eb55c7b..48fe837a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,7 +2,7 @@ name = "GeophysicalFlows"
 uuid = "44ee3b1c-bc02-53fa-8355-8e347616e15e"
 license = "MIT"
 authors = ["Navid C. Constantinou <navidcy@gmail.com>", "Gregory L. Wagner <wagner.greg@gmail.com>", "and co-contributors"]
-version = "0.14.0"
+version = "0.14.1"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
diff --git a/src/multilayerqg.jl b/src/multilayerqg.jl
index 3f42c310..940777a5 100644
--- a/src/multilayerqg.jl
+++ b/src/multilayerqg.jl
@@ -916,7 +916,7 @@ The kinetic energy at the ``j``-th fluid layer is
 while the potential energy that corresponds to the interface ``j+1/2`` (i.e., the interface 
 between the ``j``-th and ``(j+1)``-th fluid layer) is
 ```math
-𝖯𝖤_{j+1/2} = \\int \\frac1{2} \\frac{f₀^2}{g'_{j+1/2}} (ψ_j - ψ_{j+1})^2 \\frac{𝖽x 𝖽y}{L_x L_y} = \\frac1{2} \\frac{f₀^2}{g'_{j+1/2}} \\sum_{𝐤} |ψ_j - ψ_{j+1}|², \\ j = 1, ..., n-1 .
+𝖯𝖤_{j+1/2} = \\int \\frac1{2} \\frac{f₀^2}{g'_{j+1/2} H} (ψ_j - ψ_{j+1})^2 \\frac{𝖽x 𝖽y}{L_x L_y} = \\frac1{2} \\frac{f₀^2}{g'_{j+1/2} H} \\sum_{𝐤} |ψ̂_j - ψ̂_{j+1}|², \\ j = 1, ..., n-1 .
 ```
 """
 function energies(vars, params, grid, sol)
@@ -934,7 +934,7 @@ function energies(vars, params, grid, sol)
   end
 
   for j = 1:nlayers-1
-    CUDA.@allowscalar PE[j] = 1 / (2 * grid.Lx * grid.Ly) * params.f₀^2 / params.g′[j] * parsevalsum(abs2.(vars.ψh[:, :, j+1] .- vars.ψh[:, :, j]), grid)
+    CUDA.@allowscalar PE[j] = 1 / (2 * grid.Lx * grid.Ly * sum(params.H)) * params.f₀^2 / params.g′[j] * parsevalsum(abs2.(vars.ψh[:, :, j+1] .- vars.ψh[:, :, j]), grid)
   end
 
   return KE, PE
@@ -956,7 +956,7 @@ function energies(vars, params::TwoLayerParams, grid, sol)
     CUDA.@allowscalar KE[j] = @views 1 / (2 * grid.Lx * grid.Ly) * parsevalsum(abs²∇𝐮h[:, :, j], grid) * params.H[j] / sum(params.H)
   end
 
-  PE = @views 1 / (2 * grid.Lx * grid.Ly) * params.f₀^2 / params.g′ * parsevalsum(abs2.(ψ2h .- ψ1h), grid)
+  PE = @views 1 / (2 * grid.Lx * grid.Ly * sum(params.H)) * params.f₀^2 / params.g′ * parsevalsum(abs2.(ψ2h .- ψ1h), grid)
   
   return KE, PE
 end
diff --git a/test/test_multilayerqg.jl b/test/test_multilayerqg.jl
index 849a0876..5f973ddd 100644
--- a/test/test_multilayerqg.jl
+++ b/test/test_multilayerqg.jl
@@ -347,29 +347,31 @@ function test_mqg_energies(dev::Device=CPU();
 
   x, y = gridpoints(gr)
   k₀, l₀ = 2π/gr.Lx, 2π/gr.Ly # fundamental wavenumbers
-
-  nlayers = 2       # these choice of parameters give the
-  f₀, g = 1, 1      # desired PV-streamfunction relations
-  H = [0.2, 0.8]    # q1 = Δψ1 + 25*(ψ2-ψ1), and
-  ρ = [4.0, 5.0]    # q2 = Δψ2 + 25/4*(ψ1-ψ2).
-
-  prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, g, H, ρ)
+  nlayers = 2
   
+  f₀, g = 1, 1
+  H = [2.5, 7.5] # sum(params.H) = 10
+  ρ = [1.0, 2.0] # Make g′ = 1/2
+  
+  prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, g, H, ρ)
   sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid
 
-  ψ1, ψ2, q1, q2, ψ1x, ψ2x, q1x, q2x, Δψ2, Δq1, Δq2 = constructtestfields_2layer(gr)
+  ψ = zeros(dev, eltype(gr), (gr.nx, gr.ny, nlayers))
+  
+  CUDA.@allowscalar @views @. ψ[:, :, 1] =       cos(2k₀ * x) * cos(2l₀ * y)
+  CUDA.@allowscalar @views @. ψ[:, :, 2] = 1/2 * cos(2k₀ * x) * cos(2l₀ * y)
 
-  qf = zeros(dev, eltype(gr), (gr.nx, gr.ny, nlayers))
-  CUDA.@allowscalar @views qf[:, :, 1] .= q1
-  CUDA.@allowscalar @views qf[:, :, 2] .= q2
+  MultiLayerQG.set_ψ!(prob, ψ)
 
-  MultiLayerQG.set_q!(prob, qf)
+  KE1_calc = 1/4      # = H1/H
+  KE2_calc = 3/4/4    # = H2/H/4
+  PE_calc  = 1/16/10  # = 1/16/H
 
   KE, PE = MultiLayerQG.energies(prob)
 
-  return isapprox(KE[1], 61/640*1e-6, rtol=rtol_multilayerqg) &&
-         isapprox(KE[2], 3*1e-6, rtol=rtol_multilayerqg) &&
-         isapprox(PE[1], 1025/1152*1e-6, rtol=rtol_multilayerqg) &&
+  return isapprox(KE[1], KE1_calc, rtol=rtol_multilayerqg) &&
+         isapprox(KE[2], KE2_calc, rtol=rtol_multilayerqg) &&
+         isapprox(PE[1], PE_calc, rtol=rtol_multilayerqg) &&
          MultiLayerQG.addforcing!(prob.timestepper.RHS₁, sol, cl.t, cl, vs, pr, gr)==nothing
 end