diff --git a/src/dynamics/diagnostic_variables.jl b/src/dynamics/diagnostic_variables.jl
index 5337820da..ad2df9a92 100644
--- a/src/dynamics/diagnostic_variables.jl
+++ b/src/dynamics/diagnostic_variables.jl
@@ -57,7 +57,7 @@ Base.@kwdef struct DynamicsVariables{NF<:AbstractFloat, Grid<:AbstractGrid{NF}}
     trunc::Int
 
     # MULTI-PURPOSE VECTOR (a, b), work array to be reused in various places, examples:
-    # uω_coslat⁻¹, vω_coslat⁻¹ = a, b                        (all models)
+    # uω_coslat⁻¹, vω_coslat⁻¹ = a, b                        (all models)
     # uω_coslat⁻¹_grid, vω_coslat⁻¹_grid = a_grid, b_grid    (all models)
     # uh_coslat⁻¹, vh_coslat⁻¹ = a, b                        (ShallowWaterModel)
     # uh_coslat⁻¹_grid, vh_coslat⁻¹_grid = a_grid, b_grid   (ShallowWaterModel)
@@ -138,8 +138,8 @@ Base.@kwdef struct SurfaceVariables{NF<:AbstractFloat, Grid<:AbstractGrid{NF}}
     div_mean::LTM{Complex{NF}} = zeros(LTM{Complex{NF}}, trunc+2, trunc+1)    # divergence (in spectral though)
     
     precip_large_scale::Grid = zeros(Grid, nlat_half)    # large scale precipitation (for output)
-    precip_convection::Grid = zeros(Grid, nlat_half)     # convective precipitation (for output)
-    cloud_top::Grid = zeros(Grid, nlat_half)             # cloud top [hPa]
+    precip_convection::Grid = zeros(Grid, nlat_half)     # convective precipitation (for output)
+    cloud_top::Grid = zeros(Grid, nlat_half)             # cloud top [hPa]
     soil_moisture_availability::Grid = zeros(Grid, nlat_half)
 end
 
@@ -217,10 +217,10 @@ function Base.zeros(
     nthreads = Threads.nthreads()
     columns = [ColumnVariables{NF}(; nlev) for _ in 1:nthreads]
 
-    # particle work arrays
+    # particle work arrays
     particles = zeros(ParticleVariables, SG)
 
-    scale = Ref(convert(SG.NF, SG.radius))
+    scale = Ref(one(NF))
 
     return DiagnosticVariables{NF, Grid{NF}, Model}(
         layers, surface, columns, particles,
diff --git a/src/dynamics/scaling.jl b/src/dynamics/scaling.jl
index 716edc12d..801f05508 100644
--- a/src/dynamics/scaling.jl
+++ b/src/dynamics/scaling.jl
@@ -41,11 +41,15 @@ $(TYPEDSIGNATURES)
 Scales the prognostic variables vorticity and divergence with
 the Earth's radius which is used in the dynamical core."""
 function scale!(progn::PrognosticVariables,
+                diagn::DiagnosticVariables,
                 scale::Real)
     new_scale = scale/progn.scale[]     # undo previous scale and new scale in one go
     scale!(progn, :vor, new_scale)
     scale!(progn, :div, new_scale)
     progn.scale[] = scale               # store scaling information
+    diagn.scale[] = scale               # store scaling information
+    # no need to actually scale the diagnostic variables as they will be
+    # overwritten by the transform of the prognostic variables anyway
 end
 
 """
diff --git a/src/dynamics/time_integration.jl b/src/dynamics/time_integration.jl
index 33ae7a4cf..0c836f0bc 100644
--- a/src/dynamics/time_integration.jl
+++ b/src/dynamics/time_integration.jl
@@ -383,7 +383,7 @@ function time_stepping!(
     (; time_stepping) = model
 
     # SCALING: we use vorticity*radius, divergence*radius in the dynamical core
-    scale!(progn, model.spectral_grid.radius)
+    scale!(progn, diagn, model.spectral_grid.radius)
 
     # OUTPUT INITIALISATION AND STORING INITIAL CONDITIONS + FEEDBACK
     # propagate spectral state to grid variables for initial condition output
@@ -410,7 +410,7 @@ function time_stepping!(
         callback!(model.callbacks, progn, diagn, model)
     end
     
-    # UNSCALE, CLOSE, FINISH
+    # UNSCALE, CLOSE, FINISH
     finish!(feedback)                       # finish the progress meter, do first for benchmark accuracy
     unscale!(progn)                         # undo radius-scaling for vor, div from the dynamical core
     unscale!(diagn)                         # undo radius-scaling for vor, div from the dynamical core