diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl
index 9043ed4cec..d56ece72f8 100644
--- a/src/prognostic_equations/remaining_tendency.jl
+++ b/src/prognostic_equations/remaining_tendency.jl
@@ -24,7 +24,6 @@ end
 
 NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
     viscous_sponge_tendency!(Yₜ, Y, p, t, p.atmos.viscous_sponge)
-    surface_temp_tendency!(Yₜ, Y, p, t, p.atmos.surface_model)
 
     # Vertical tendencies
     rayleigh_sponge_tendency!(Yₜ, Y, p, t, p.atmos.rayleigh_sponge)
@@ -80,6 +79,10 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
         p.atmos.turbconv_model,
     )
 
+    # NOTE: Precipitation tendencies should be applied before calling this function,
+    # because precipitation cache is used in this function
+    surface_temp_tendency!(Yₜ, Y, p, t, p.atmos.surface_model)
+
     # NOTE: All ρa tendencies should be applied before calling this function
     pressure_work_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
 
diff --git a/src/solver/solve.jl b/src/solver/solve.jl
index 7cfe75e917..d044b4b79b 100644
--- a/src/solver/solve.jl
+++ b/src/solver/solve.jl
@@ -218,14 +218,14 @@ Return:
 """
 function check_conservation(sol)
     # energy
-    energy_total = sum(sol.u[end].c.ρe_tot)
+    energy_total = sum(sol.u[1].c.ρe_tot)
     energy_atmos_change = sum(sol.u[end].c.ρe_tot) - sum(sol.u[1].c.ρe_tot)
     p = sol.prob.p
     sfc = p.atmos.surface_model
     if sfc isa PrognosticSurfaceTemperature
         sfc_cρh = sfc.ρ_ocean * sfc.cp_ocean * sfc.depth_ocean
         energy_total +=
-            horizontal_integral_at_boundary(sol.u[end].sfc.T .* sfc_cρh)
+            horizontal_integral_at_boundary(sol.u[1].sfc.T .* sfc_cρh)
         energy_surface_change =
             horizontal_integral_at_boundary(
                 sol.u[end].sfc.T .- sol.u[1].sfc.T,
@@ -254,6 +254,10 @@ function check_conservation(sol)
             sol.u[end].sfc.water .- sol.u[1].sfc.water,
         )
 
+        mass_conservation =
+            (sum(sol.u[end].c.ρ) - sum(sol.u[1].c.ρ) + water_surface_change) /
+            sum(sol.u[1].c.ρ)
+
         water_conservation =
             abs(water_atmos_change + water_surface_change) / water_total
     end