diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index d37f265030..8f0f0e7462 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -1297,7 +1297,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
   if ((CS%thickness_diffuse .or. CS%interface_filter) .and. &
       .not.CS%thickness_diffuse_first) then
 
-    if (CS%debug) call hchksum(h,"Pre-thickness_diffuse h", G%HI, haloshift=0, scale=GV%H_to_MKS)
+    if (CS%debug) call hchksum(h,"Pre-thickness_diffuse h", G%HI, haloshift=0, unscale=GV%H_to_MKS)
 
     if (CS%thickness_diffuse) then
       call cpu_clock_begin(id_clock_thick_diff)
@@ -1306,7 +1306,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
       call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt, G, GV, US, &
                              CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp)
 
-      if (CS%debug) call hchksum(h,"Post-thickness_diffuse h", G%HI, haloshift=1, scale=GV%H_to_MKS)
+      if (CS%debug) call hchksum(h,"Post-thickness_diffuse h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
       call cpu_clock_end(id_clock_thick_diff)
       call pass_var(h, G%Domain, clock=id_clock_pass, halo=max(2,CS%cont_stencil))
       if (showCallTree) call callTree_waypoint("finished thickness_diffuse (step_MOM)")
@@ -1327,9 +1327,9 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
   ! apply the submesoscale mixed layer restratification parameterization
   if (CS%mixedlayer_restrat) then
     if (CS%debug) then
-      call hchksum(h,"Pre-mixedlayer_restrat h", G%HI, haloshift=1, scale=GV%H_to_MKS)
+      call hchksum(h,"Pre-mixedlayer_restrat h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
       call uvchksum("Pre-mixedlayer_restrat uhtr", &
-                    CS%uhtr, CS%vhtr, G%HI, haloshift=0, scale=GV%H_to_MKS*US%L_to_m**2)
+                    CS%uhtr, CS%vhtr, G%HI, haloshift=0, unscale=GV%H_to_MKS*US%L_to_m**2)
     endif
     call cpu_clock_begin(id_clock_ml_restrat)
     call mixedlayer_restrat(h, CS%uhtr, CS%vhtr, CS%tv, forces, dt, CS%visc%MLD, CS%visc%h_ML, &
@@ -1337,9 +1337,9 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
     call cpu_clock_end(id_clock_ml_restrat)
     call pass_var(h, G%Domain, clock=id_clock_pass, halo=max(2,CS%cont_stencil))
     if (CS%debug) then
-      call hchksum(h,"Post-mixedlayer_restrat h", G%HI, haloshift=1, scale=GV%H_to_MKS)
+      call hchksum(h,"Post-mixedlayer_restrat h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
       call uvchksum("Post-mixedlayer_restrat [uv]htr", &
-                    CS%uhtr, CS%vhtr, G%HI, haloshift=0, scale=GV%H_to_MKS*US%L_to_m**2)
+                    CS%uhtr, CS%vhtr, G%HI, haloshift=0, unscale=GV%H_to_MKS*US%L_to_m**2)
     endif
   endif
 
@@ -1410,15 +1410,15 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local)
 
   if (CS%debug) then
     call cpu_clock_begin(id_clock_other)
-    call hchksum(h,"Pre-advection h", G%HI, haloshift=1, scale=GV%H_to_MKS)
+    call hchksum(h,"Pre-advection h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
     call uvchksum("Pre-advection uhtr", CS%uhtr, CS%vhtr, G%HI, &
-                  haloshift=0, scale=GV%H_to_MKS*US%L_to_m**2)
-    if (associated(CS%tv%T)) call hchksum(CS%tv%T, "Pre-advection T", G%HI, haloshift=1, scale=US%C_to_degC)
-    if (associated(CS%tv%S)) call hchksum(CS%tv%S, "Pre-advection S", G%HI, haloshift=1, scale=US%S_to_ppt)
+                  haloshift=0, unscale=GV%H_to_MKS*US%L_to_m**2)
+    if (associated(CS%tv%T)) call hchksum(CS%tv%T, "Pre-advection T", G%HI, haloshift=1, unscale=US%C_to_degC)
+    if (associated(CS%tv%S)) call hchksum(CS%tv%S, "Pre-advection S", G%HI, haloshift=1, unscale=US%S_to_ppt)
     if (associated(CS%tv%frazil)) call hchksum(CS%tv%frazil, "Pre-advection frazil", G%HI, haloshift=0, &
-                                               scale=US%Q_to_J_kg*US%RZ_to_kg_m2)
+                                               unscale=US%Q_to_J_kg*US%RZ_to_kg_m2)
     if (associated(CS%tv%salt_deficit)) call hchksum(CS%tv%salt_deficit, &
-                   "Pre-advection salt deficit", G%HI, haloshift=0, scale=US%S_to_ppt*US%RZ_to_kg_m2)
+                   "Pre-advection salt deficit", G%HI, haloshift=0, unscale=US%S_to_ppt*US%RZ_to_kg_m2)
   ! call MOM_thermo_chksum("Pre-advection ", CS%tv, G, US)
     call cpu_clock_end(id_clock_other)
   endif
@@ -1593,10 +1593,10 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
   call cpu_clock_begin(id_clock_thermo)
   if (.not.CS%adiabatic) then
     if (CS%debug) then
-      call uvchksum("Pre-diabatic [uv]", u, v, G%HI, haloshift=2, scale=US%L_T_to_m_s)
-      call hchksum(h,"Pre-diabatic h", G%HI, haloshift=1, scale=GV%H_to_MKS)
+      call uvchksum("Pre-diabatic [uv]", u, v, G%HI, haloshift=2, unscale=US%L_T_to_m_s)
+      call hchksum(h,"Pre-diabatic h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
       call uvchksum("Pre-diabatic [uv]h", CS%uhtr, CS%vhtr, G%HI, &
-                    haloshift=0, scale=GV%H_to_MKS*US%L_to_m**2)
+                    haloshift=0, unscale=GV%H_to_MKS*US%L_to_m**2)
     ! call MOM_state_chksum("Pre-diabatic ", u, v, h, CS%uhtr, CS%vhtr, G, GV, vel_scale=1.0)
       call MOM_thermo_chksum("Pre-diabatic ", tv, G, US, haloshift=0)
       call check_redundant("Pre-diabatic ", u, v, G, unscale=US%L_T_to_m_s)
@@ -1634,8 +1634,8 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
 
       if (CS%debug) then
         call MOM_state_chksum("Pre-ALE ", u, v, h, CS%uh, CS%vh, G, GV, US, omit_corners=.true.)
-        call hchksum(tv%T,"Pre-ALE T", G%HI, haloshift=1, omit_corners=.true., scale=US%C_to_degC)
-        call hchksum(tv%S,"Pre-ALE S", G%HI, haloshift=1, omit_corners=.true., scale=US%S_to_ppt)
+        call hchksum(tv%T,"Pre-ALE T", G%HI, haloshift=1, omit_corners=.true., unscale=US%C_to_degC)
+        call hchksum(tv%S,"Pre-ALE S", G%HI, haloshift=1, omit_corners=.true., unscale=US%S_to_ppt)
         call check_redundant("Pre-ALE ", u, v, G, unscale=US%L_T_to_m_s)
       endif
       call cpu_clock_begin(id_clock_ALE)
@@ -1720,8 +1720,8 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
 
     if (CS%debug .and. CS%use_ALE_algorithm) then
       call MOM_state_chksum("Post-ALE ", u, v, h, CS%uh, CS%vh, G, GV, US)
-      call hchksum(tv%T, "Post-ALE T", G%HI, haloshift=1, scale=US%C_to_degC)
-      call hchksum(tv%S, "Post-ALE S", G%HI, haloshift=1, scale=US%S_to_ppt)
+      call hchksum(tv%T, "Post-ALE T", G%HI, haloshift=1, unscale=US%C_to_degC)
+      call hchksum(tv%S, "Post-ALE S", G%HI, haloshift=1, unscale=US%S_to_ppt)
       call check_redundant("Post-ALE ", u, v, G, unscale=US%L_T_to_m_s)
     endif
 
@@ -1734,18 +1734,18 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
     call postALE_tracer_diagnostics(CS%tracer_Reg, G, GV, CS%diag, dtdia)
 
     if (CS%debug) then
-      call uvchksum("Post-diabatic u", u, v, G%HI, haloshift=2, scale=US%L_T_to_m_s)
-      call hchksum(h, "Post-diabatic h", G%HI, haloshift=1, scale=GV%H_to_MKS)
+      call uvchksum("Post-diabatic u", u, v, G%HI, haloshift=2, unscale=US%L_T_to_m_s)
+      call hchksum(h, "Post-diabatic h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
       call uvchksum("Post-diabatic [uv]h", CS%uhtr, CS%vhtr, G%HI, &
-                    haloshift=0, scale=GV%H_to_MKS*US%L_to_m**2)
+                    haloshift=0, unscale=GV%H_to_MKS*US%L_to_m**2)
     ! call MOM_state_chksum("Post-diabatic ", u, v, &
     !                       h, CS%uhtr, CS%vhtr, G, GV, haloshift=1)
-      if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", G%HI, haloshift=1, scale=US%C_to_degC)
-      if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", G%HI, haloshift=1, scale=US%S_to_ppt)
+      if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", G%HI, haloshift=1, unscale=US%C_to_degC)
+      if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", G%HI, haloshift=1, unscale=US%S_to_ppt)
       if (associated(tv%frazil)) call hchksum(tv%frazil, "Post-diabatic frazil", G%HI, haloshift=0, &
-                                              scale=US%Q_to_J_kg*US%RZ_to_kg_m2)
+                                              unscale=US%Q_to_J_kg*US%RZ_to_kg_m2)
       if (associated(tv%salt_deficit)) call hchksum(tv%salt_deficit, &
-                               "Post-diabatic salt deficit", G%HI, haloshift=0, scale=US%RZ_to_kg_m2)
+                               "Post-diabatic salt deficit", G%HI, haloshift=0, unscale=US%RZ_to_kg_m2)
     ! call MOM_thermo_chksum("Post-diabatic ", tv, G, US)
       call check_redundant("Post-diabatic ", u, v, G, unscale=US%L_T_to_m_s)
     endif
@@ -1765,8 +1765,8 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
       call create_group_pass(pass_T_S, tv%S, G%Domain, To_All+Omit_Corners, halo=dynamics_stencil)
       call do_group_pass(pass_T_S, G%Domain, clock=id_clock_pass)
       if (CS%debug) then
-        if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", G%HI, haloshift=1, scale=US%C_to_degC)
-        if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", G%HI, haloshift=1, scale=US%S_to_ppt)
+        if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", G%HI, haloshift=1, unscale=US%C_to_degC)
+        if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", G%HI, haloshift=1, unscale=US%S_to_ppt)
       endif
 
       ! Update derived thermodynamic quantities.
@@ -3031,7 +3031,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
 
   if (use_ice_shelf .and. CS%debug) then
     call hchksum(CS%frac_shelf_h, "MOM:frac_shelf_h", G%HI, haloshift=0)
-    call hchksum(CS%mass_shelf, "MOM:mass_shelf", G%HI, haloshift=0,scale=US%RZ_to_kg_m2)
+    call hchksum(CS%mass_shelf, "MOM:mass_shelf", G%HI, haloshift=0, unscale=US%RZ_to_kg_m2)
   endif
 
   call cpu_clock_end(id_clock_MOM_init)
@@ -3070,8 +3070,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
     ! \todo This block exists for legacy reasons and we should phase it out of
     ! all examples. !###
     if (CS%debug) then
-      call uvchksum("Pre ALE adjust init cond [uv]", CS%u, CS%v, G%HI, haloshift=1)
-      call hchksum(CS%h,"Pre ALE adjust init cond h", G%HI, haloshift=1, scale=GV%H_to_MKS)
+      call uvchksum("Pre ALE adjust init cond [uv]", CS%u, CS%v, G%HI, haloshift=1, unscale=US%L_T_to_m_s)
+      call hchksum(CS%h,"Pre ALE adjust init cond h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
     endif
     call callTree_waypoint("Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)")
     call adjustGridForIntegrity(CS%ALE_CSp, G, GV, CS%h )
@@ -3128,11 +3128,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
     call cpu_clock_end(id_clock_pass_init)
 
     if (CS%debug) then
-      call uvchksum("Post ALE adjust init cond [uv]", CS%u, CS%v, G%HI, haloshift=1)
-      call hchksum(CS%h, "Post ALE adjust init cond h", G%HI, haloshift=2, scale=GV%H_to_MKS)
+      call uvchksum("Post ALE adjust init cond [uv]", CS%u, CS%v, G%HI, haloshift=1, unscale=US%L_T_to_m_s)
+      call hchksum(CS%h, "Post ALE adjust init cond h", G%HI, haloshift=2, unscale=GV%H_to_MKS)
       if (use_temperature) then
-        call hchksum(CS%tv%T, "Post ALE adjust init cond T", G%HI, haloshift=2, scale=US%C_to_degC)
-        call hchksum(CS%tv%S, "Post ALE adjust init cond S", G%HI, haloshift=2, scale=US%S_to_ppt)
+        call hchksum(CS%tv%T, "Post ALE adjust init cond T", G%HI, haloshift=2, unscale=US%C_to_degC)
+        call hchksum(CS%tv%S, "Post ALE adjust init cond S", G%HI, haloshift=2, unscale=US%S_to_ppt)
       endif
     endif
   endif
@@ -4130,7 +4130,7 @@ subroutine get_ocean_stocks(CS, mass, heat, salt, on_PE_only)
     heat = CS%US%Q_to_J_kg*CS%tv%C_p * &
            global_mass_integral(CS%h, CS%G, CS%GV, CS%tv%T, on_PE_only=on_PE_only, tmp_scale=CS%US%C_to_degC)
   if (present(salt)) &
-    salt = 1.0e-3 * global_mass_integral(CS%h, CS%G, CS%GV, CS%tv%S, on_PE_only=on_PE_only, scale=CS%US%S_to_ppt)
+    salt = 1.0e-3 * global_mass_integral(CS%h, CS%G, CS%GV, CS%tv%S, on_PE_only=on_PE_only, unscale=CS%US%S_to_ppt)
 
 end subroutine get_ocean_stocks
 
diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 7d112de337..a132200759 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -1725,23 +1725,23 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
 
   if (CS%debug) then
     call uvchksum("BT [uv]hbt", uhbt, vhbt, CS%debug_BT_HI, haloshift=0, &
-                  scale=US%s_to_T*US%L_to_m**2*GV%H_to_m)
-    call uvchksum("BT Initial [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=0, scale=US%L_T_to_m_s)
-    call hchksum(eta, "BT Initial eta", CS%debug_BT_HI, haloshift=0, scale=GV%H_to_MKS)
+                  unscale=US%s_to_T*US%L_to_m**2*GV%H_to_m)
+    call uvchksum("BT Initial [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=0, unscale=US%L_T_to_m_s)
+    call hchksum(eta, "BT Initial eta", CS%debug_BT_HI, haloshift=0, unscale=GV%H_to_MKS)
     call uvchksum("BT BT_force_[uv]", BT_force_u, BT_force_v, &
-                  CS%debug_BT_HI, haloshift=0, scale=US%L_T2_to_m_s2)
+                  CS%debug_BT_HI, haloshift=0, unscale=US%L_T2_to_m_s2)
     if (interp_eta_PF) then
-      call hchksum(eta_PF_1, "BT eta_PF_1",CS%debug_BT_HI,haloshift=0, scale=GV%H_to_MKS)
-      call hchksum(d_eta_PF, "BT d_eta_PF",CS%debug_BT_HI,haloshift=0, scale=GV%H_to_MKS)
+      call hchksum(eta_PF_1, "BT eta_PF_1",CS%debug_BT_HI,haloshift=0, unscale=GV%H_to_MKS)
+      call hchksum(d_eta_PF, "BT d_eta_PF",CS%debug_BT_HI,haloshift=0, unscale=GV%H_to_MKS)
     else
-      call hchksum(eta_PF, "BT eta_PF",CS%debug_BT_HI,haloshift=0, scale=GV%H_to_MKS)
-      call hchksum(eta_PF_in, "BT eta_PF_in",G%HI,haloshift=0, scale=GV%H_to_MKS)
+      call hchksum(eta_PF, "BT eta_PF",CS%debug_BT_HI,haloshift=0, unscale=GV%H_to_MKS)
+      call hchksum(eta_PF_in, "BT eta_PF_in",G%HI,haloshift=0, unscale=GV%H_to_MKS)
     endif
-    call uvchksum("BT Cor_ref_[uv]", Cor_ref_u, Cor_ref_v, CS%debug_BT_HI, haloshift=0, scale=US%L_T2_to_m_s2)
+    call uvchksum("BT Cor_ref_[uv]", Cor_ref_u, Cor_ref_v, CS%debug_BT_HI, haloshift=0, unscale=US%L_T2_to_m_s2)
     call uvchksum("BT [uv]hbt0", uhbt0, vhbt0, CS%debug_BT_HI, haloshift=0, &
-                  scale=US%L_to_m**2*US%s_to_T*GV%H_to_m)
+                  unscale=US%L_to_m**2*US%s_to_T*GV%H_to_m)
     if (.not. use_BT_cont) then
-      call uvchksum("BT Dat[uv]", Datu, Datv, CS%debug_BT_HI, haloshift=1, scale=US%L_to_m*GV%H_to_m)
+      call uvchksum("BT Dat[uv]", Datu, Datv, CS%debug_BT_HI, haloshift=1, unscale=US%L_to_m*GV%H_to_m)
     endif
     call uvchksum("BT wt_[uv]", wt_u, wt_v, G%HI, haloshift=0, &
                   symmetric=.true., omit_corners=.true., scalar_pair=.true.)
@@ -1749,9 +1749,9 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
                   symmetric=.true., omit_corners=.true., scalar_pair=.true.)
     call uvchksum("BT visc_rem_[uv]", visc_rem_u, visc_rem_v, G%HI, haloshift=0, &
                   symmetric=.true., omit_corners=.true., scalar_pair=.true.)
-    call uvchksum("BT bc_accel_[uv]", bc_accel_u, bc_accel_v, G%HI, haloshift=0, scale=US%L_T2_to_m_s2)
+    call uvchksum("BT bc_accel_[uv]", bc_accel_u, bc_accel_v, G%HI, haloshift=0, unscale=US%L_T2_to_m_s2)
     call uvchksum("BT IDat[uv]", CS%IDatu, CS%IDatv, G%HI, haloshift=0, &
-                  scale=GV%m_to_H, scalar_pair=.true.)
+                  unscale=GV%m_to_H, scalar_pair=.true.)
     call uvchksum("BT visc_rem_[uv]", visc_rem_u, visc_rem_v, G%HI, &
                   haloshift=1, scalar_pair=.true.)
   endif
@@ -2324,22 +2324,22 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
     if (CS%debug_bt) then
       write(mesg,'("BT vel update ",I4)') n
       call uvchksum(trim(mesg)//" PF[uv]", PFu, PFv, CS%debug_BT_HI, haloshift=iev-ie, &
-                    scale=US%L_T_to_m_s*US%s_to_T)
+                    unscale=US%L_T_to_m_s*US%s_to_T)
       call uvchksum(trim(mesg)//" Cor_[uv]", Cor_u, Cor_v, CS%debug_BT_HI, haloshift=iev-ie, &
-                    scale=US%L_T_to_m_s*US%s_to_T)
+                    unscale=US%L_T_to_m_s*US%s_to_T)
       call uvchksum(trim(mesg)//" BT_force_[uv]", BT_force_u, BT_force_v, CS%debug_BT_HI, haloshift=iev-ie, &
-                    scale=US%L_T_to_m_s*US%s_to_T)
+                    unscale=US%L_T_to_m_s*US%s_to_T)
       call uvchksum(trim(mesg)//" BT_rem_[uv]", BT_rem_u, BT_rem_v, CS%debug_BT_HI, &
                     haloshift=iev-ie, scalar_pair=.true.)
       call uvchksum(trim(mesg)//" [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=iev-ie, &
-                    scale=US%L_T_to_m_s)
+                    unscale=US%L_T_to_m_s)
       call uvchksum(trim(mesg)//" [uv]bt_trans", ubt_trans, vbt_trans, CS%debug_BT_HI, haloshift=iev-ie, &
-                    scale=US%L_T_to_m_s)
+                    unscale=US%L_T_to_m_s)
       call uvchksum(trim(mesg)//" [uv]hbt", uhbt, vhbt, CS%debug_BT_HI, haloshift=iev-ie, &
-                    scale=US%s_to_T*US%L_to_m**2*GV%H_to_m)
+                    unscale=US%s_to_T*US%L_to_m**2*GV%H_to_m)
       if (integral_BT_cont) &
         call uvchksum(trim(mesg)//" [uv]hbt_int", uhbt_int, vhbt_int, CS%debug_BT_HI, haloshift=iev-ie, &
-                      scale=US%L_to_m**2*GV%H_to_m)
+                      unscale=US%L_to_m**2*GV%H_to_m)
     endif
 
     if (find_PF) then
@@ -2426,10 +2426,10 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
 
     if (CS%debug_bt) then
       call uvchksum("BT [uv]hbt just after OBC", uhbt, vhbt, CS%debug_BT_HI, haloshift=iev-ie, &
-                    scale=US%s_to_T*US%L_to_m**2*GV%H_to_m)
+                    unscale=US%s_to_T*US%L_to_m**2*GV%H_to_m)
       if (integral_BT_cont) &
         call uvchksum("BT [uv]hbt_int just after OBC", uhbt_int, vhbt_int, CS%debug_BT_HI, &
-                      haloshift=iev-ie, scale=US%L_to_m**2*GV%H_to_m)
+                      haloshift=iev-ie, unscale=US%L_to_m**2*GV%H_to_m)
     endif
 
     if (integral_BT_cont) then
@@ -2463,8 +2463,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
     if (CS%debug_bt) then
       write(mesg,'("BT step ",I4)') n
       call uvchksum(trim(mesg)//" [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=iev-ie, &
-                    scale=US%L_T_to_m_s)
-      call hchksum(eta, trim(mesg)//" eta", CS%debug_BT_HI, haloshift=iev-ie, scale=GV%H_to_MKS)
+                    unscale=US%L_T_to_m_s)
+      call hchksum(eta, trim(mesg)//" eta", CS%debug_BT_HI, haloshift=iev-ie, unscale=GV%H_to_MKS)
     endif
 
     if (GV%Boussinesq) then
@@ -2947,8 +2947,8 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
   CS%dtbt_max = dtbt_max
 
   if (CS%debug) then
-    call chksum0(CS%dtbt, "End set_dtbt dtbt", scale=US%T_to_s)
-    call chksum0(CS%dtbt_max, "End set_dtbt dtbt_max", scale=US%T_to_s)
+    call chksum0(CS%dtbt, "End set_dtbt dtbt", unscale=US%T_to_s)
+    call chksum0(CS%dtbt_max, "End set_dtbt dtbt_max", unscale=US%T_to_s)
   endif
 
 end subroutine set_dtbt
@@ -3698,9 +3698,9 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
                   scalar_pair=.true.)
     if (present(h_u) .and. present(h_v)) &
       call uvchksum("btcalc h_[uv]", h_u, h_v, G%HI, haloshift=0, &
-                    symmetric=.true., omit_corners=.true., scale=GV%H_to_MKS, &
+                    symmetric=.true., omit_corners=.true., unscale=GV%H_to_MKS, &
                     scalar_pair=.true.)
-    call hchksum(h, "btcalc h",G%HI, haloshift=1, scale=GV%H_to_MKS)
+    call hchksum(h, "btcalc h",G%HI, haloshift=1, unscale=GV%H_to_MKS)
   endif
 
 end subroutine btcalc
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index 671f707583..2296496144 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -584,7 +584,7 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
   call disable_averaging(CS%diag)
 
   if (CS%debug) then
-    call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
+    call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
   endif
   call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
   call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
@@ -668,10 +668,10 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
   call cpu_clock_end(id_clock_mom_update)
 
   if (CS%debug) then
-    call uvchksum("Predictor 1 [uv]", up, vp, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
-    call hchksum(h, "Predictor 1 h", G%HI, haloshift=1, scale=GV%H_to_MKS)
+    call uvchksum("Predictor 1 [uv]", up, vp, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
+    call hchksum(h, "Predictor 1 h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
     call uvchksum("Predictor 1 [uv]h", uh, vh, G%HI,haloshift=2, &
-                  symmetric=sym, scale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
+                  symmetric=sym, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
 !   call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1)
     call MOM_accel_chksum("Predictor accel", CS%CAu_pred, CS%CAv_pred, CS%PFu, CS%PFv, &
              CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, symmetric=sym)
@@ -685,7 +685,7 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
 ! u_av  <- u_av  + dt_pred d/dz visc d/dz u_av
   call cpu_clock_begin(id_clock_vertvisc)
   if (CS%debug) then
-    call uvchksum("0 before vertvisc: [uv]p", up, vp, G%HI,haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
+    call uvchksum("0 before vertvisc: [uv]p", up, vp, G%HI,haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
   endif
 
   if (CS%fpmix) then
@@ -754,12 +754,12 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
   if (associated(CS%OBC)) then
 
     if (CS%debug) &
-      call uvchksum("Pre OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
+      call uvchksum("Pre OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
 
     call radiation_open_bdry_conds(CS%OBC, u_av, u_old_rad_OBC, v_av, v_old_rad_OBC, G, GV, US, dt_pred)
 
     if (CS%debug) &
-      call uvchksum("Post OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
+      call uvchksum("Post OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
 
     ! These should be done with a pass that excludes uh & vh.
 !   call do_group_pass(CS%pass_hp_uv, G%Domain, clock=id_clock_pass)
@@ -837,8 +837,8 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
 
   if (CS%debug) then
     call MOM_state_chksum("Predictor ", up, vp, hp, uh, vh, G, GV, US, symmetric=sym)
-    call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
-    call hchksum(h_av, "Predictor avg h", G%HI, haloshift=2, scale=GV%H_to_MKS)
+    call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
+    call hchksum(h_av, "Predictor avg h", G%HI, haloshift=2, unscale=GV%H_to_MKS)
   ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
     call check_redundant("Predictor up ", up, vp, G, unscale=US%L_T_to_m_s)
     call check_redundant("Predictor uh ", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
@@ -929,10 +929,10 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
   call cpu_clock_end(id_clock_mom_update)
 
   if (CS%debug) then
-    call uvchksum("Corrector 1 [uv]", u_inst, v_inst, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
-    call hchksum(h, "Corrector 1 h", G%HI, haloshift=1, scale=GV%H_to_MKS)
+    call uvchksum("Corrector 1 [uv]", u_inst, v_inst, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
+    call hchksum(h, "Corrector 1 h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
     call uvchksum("Corrector 1 [uv]h", uh, vh, G%HI, haloshift=2, &
-                  symmetric=sym, scale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
+                  symmetric=sym, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
   ! call MOM_state_chksum("Corrector 1", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1)
     call MOM_accel_chksum("Corrector accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv, &
                           CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, &
@@ -1158,8 +1158,8 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
 
   if (CS%debug) then
     call MOM_state_chksum("Corrector ", u_inst, v_inst, h, uh, vh, G, GV, US, symmetric=sym)
-    call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
-    call hchksum(h_av, "Corrector avg h", G%HI, haloshift=1, scale=GV%H_to_MKS)
+    call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
+    call hchksum(h_av, "Corrector avg h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
  !  call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
   endif
 
diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90
index 9c5c248c3c..021b3ceb0e 100644
--- a/src/core/MOM_dynamics_split_RK2b.F90
+++ b/src/core/MOM_dynamics_split_RK2b.F90
@@ -600,7 +600,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
   call disable_averaging(CS%diag)
 
   if (CS%debug) then
-    call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
+    call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
   endif
   call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
   call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
@@ -691,10 +691,10 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
   call cpu_clock_end(id_clock_mom_update)
 
   if (CS%debug) then
-    call uvchksum("Predictor 1 [uv]", up, vp, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
-    call hchksum(h, "Predictor 1 h", G%HI, haloshift=1, scale=GV%H_to_MKS)
+    call uvchksum("Predictor 1 [uv]", up, vp, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
+    call hchksum(h, "Predictor 1 h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
     call uvchksum("Predictor 1 [uv]h", uh, vh, G%HI,haloshift=2, &
-                  symmetric=sym, scale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
+                  symmetric=sym, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
 !   call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1)
     call MOM_accel_chksum("Predictor accel", CS%CAu_pred, CS%CAv_pred, CS%PFu, CS%PFv, &
              CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, symmetric=sym)
@@ -708,7 +708,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
 ! u_av  <- u_av  + dt_pred d/dz visc d/dz u_av
   call cpu_clock_begin(id_clock_vertvisc)
   if (CS%debug) then
-    call uvchksum("0 before vertvisc: [uv]p", up, vp, G%HI,haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
+    call uvchksum("0 before vertvisc: [uv]p", up, vp, G%HI,haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
   endif
 
   !  if (CS%fpmix) then
@@ -770,12 +770,12 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
 
   if (associated(CS%OBC)) then
     if (CS%debug) &
-      call uvchksum("Pre OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
+      call uvchksum("Pre OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
 
     call radiation_open_bdry_conds(CS%OBC, u_av, u_old_rad_OBC, v_av, v_old_rad_OBC, G, GV, US, dt_pred)
 
     if (CS%debug) &
-      call uvchksum("Post OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
+      call uvchksum("Post OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
   endif
 
   ! h_av = (h + hp)/2
@@ -834,8 +834,8 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
 
   if (CS%debug) then
     call MOM_state_chksum("Predictor ", up, vp, hp, uh, vh, G, GV, US, symmetric=sym)
-    call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
-    call hchksum(h_av, "Predictor avg h", G%HI, haloshift=2, scale=GV%H_to_MKS)
+    call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
+    call hchksum(h_av, "Predictor avg h", G%HI, haloshift=2, unscale=GV%H_to_MKS)
   ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
     call check_redundant("Predictor up ", up, vp, G, unscale=US%L_T_to_m_s)
     call check_redundant("Predictor uh ", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
@@ -923,10 +923,10 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
   call cpu_clock_end(id_clock_mom_update)
 
   if (CS%debug) then
-    call uvchksum("Corrector 1 [uv]", u_inst, v_inst, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
-    call hchksum(h, "Corrector 1 h", G%HI, haloshift=1, scale=GV%H_to_MKS)
+    call uvchksum("Corrector 1 [uv]", u_inst, v_inst, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
+    call hchksum(h, "Corrector 1 h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
     call uvchksum("Corrector 1 [uv]h", uh, vh, G%HI, haloshift=2, &
-                  symmetric=sym, scale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
+                  symmetric=sym, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
   ! call MOM_state_chksum("Corrector 1", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1)
     call MOM_accel_chksum("Corrector accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv, &
                           CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, &
@@ -1119,7 +1119,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
 
   if (CS%debug) then
     call MOM_state_chksum("Corrector ", u_av, v_av, h, uh, vh, G, GV, US, symmetric=sym)
-    ! call uvchksum("Corrector inst [uv]", u_inst, v_inst, G%HI, symmetric=sym, scale=US%L_T_to_m_s)
+    ! call uvchksum("Corrector inst [uv]", u_inst, v_inst, G%HI, symmetric=sym, unscale=US%L_T_to_m_s)
   endif
 
   if (showCallTree) call callTree_leave("step_MOM_dyn_split_RK2b()")
diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90
index 3891c86e3a..6e272f7b41 100644
--- a/src/core/MOM_interface_heights.F90
+++ b/src/core/MOM_interface_heights.F90
@@ -307,11 +307,11 @@ subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug)
     tv%valid_SpV_halo = halos
 
     if (do_debug) then
-      call hchksum(h, "derived_thermo h", G%HI, haloshift=halos, scale=GV%H_to_MKS)
+      call hchksum(h, "derived_thermo h", G%HI, haloshift=halos, unscale=GV%H_to_MKS)
       if (associated(tv%p_surf)) call hchksum(tv%p_surf, "derived_thermo p_surf", G%HI, &
-                                              haloshift=halos, scale=US%RL2_T2_to_Pa)
-      call hchksum(tv%T, "derived_thermo T", G%HI, haloshift=halos, scale=US%C_to_degC)
-      call hchksum(tv%S, "derived_thermo S", G%HI, haloshift=halos, scale=US%S_to_ppt)
+                                              haloshift=halos, unscale=US%RL2_T2_to_Pa)
+      call hchksum(tv%T, "derived_thermo T", G%HI, haloshift=halos, unscale=US%C_to_degC)
+      call hchksum(tv%S, "derived_thermo S", G%HI, haloshift=halos, unscale=US%S_to_ppt)
     endif
   elseif (allocated(tv%Spv_avg)) then
     do k=1,nz ; SpV_lay(k) = 1.0 / GV%Rlay(k) ; enddo
diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90
index 3674e63c31..9b8d26cb09 100644
--- a/src/core/MOM_open_boundary.F90
+++ b/src/core/MOM_open_boundary.F90
@@ -3318,22 +3318,22 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US,
     sym = G%Domain%symmetric
     if (OBC%radiation_BCs_exist_globally) then
       call uvchksum("radiation_OBCs: OBC%r[xy]_normal", OBC%rx_normal, OBC%ry_normal, G%HI, &
-                  haloshift=0, symmetric=sym, scale=1.0)
+                  haloshift=0, symmetric=sym, unscale=1.0)
     endif
     if (OBC%oblique_BCs_exist_globally) then
       call uvchksum("radiation_OBCs: OBC%r[xy]_oblique_[uv]", OBC%rx_oblique_u, OBC%ry_oblique_v, G%HI, &
-                  haloshift=0, symmetric=sym, scale=1.0/US%L_T_to_m_s**2)
+                  haloshift=0, symmetric=sym, unscale=1.0/US%L_T_to_m_s**2)
       call uvchksum("radiation_OBCs: OBC%r[yx]_oblique_[uv]", OBC%ry_oblique_u, OBC%rx_oblique_v, G%HI, &
-                  haloshift=0, symmetric=sym, scale=1.0/US%L_T_to_m_s**2)
+                  haloshift=0, symmetric=sym, unscale=1.0/US%L_T_to_m_s**2)
       call uvchksum("radiation_OBCs: OBC%cff_normal_[uv]", OBC%cff_normal_u, OBC%cff_normal_v, G%HI, &
-                  haloshift=0, symmetric=sym, scale=1.0/US%L_T_to_m_s**2)
+                  haloshift=0, symmetric=sym, unscale=1.0/US%L_T_to_m_s**2)
     endif
     if (OBC%ntr == 0) return
     if (.not. allocated (OBC%tres_x) .or. .not. allocated (OBC%tres_y)) return
     do m=1,OBC%ntr
       write(var_num,'(I3.3)') m
       call uvchksum("radiation_OBCs: OBC%tres_[xy]_"//var_num, OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%HI, &
-                    haloshift=0, symmetric=sym, scale=1.0)
+                    haloshift=0, symmetric=sym, unscale=1.0)
     enddo
   endif
 
diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index 37c1607826..5b7740230a 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -597,15 +597,15 @@ subroutine MOM_thermovar_chksum(mesg, tv, G, US)
   ! counts, there must be no redundant points, so all variables use is..ie
   ! and js...je as their extent.
   if (associated(tv%T)) &
-    call hchksum(tv%T, mesg//" tv%T", G%HI, scale=US%C_to_degC)
+    call hchksum(tv%T, mesg//" tv%T", G%HI, unscale=US%C_to_degC)
   if (associated(tv%S)) &
-    call hchksum(tv%S, mesg//" tv%S", G%HI, scale=US%S_to_ppt)
+    call hchksum(tv%S, mesg//" tv%S", G%HI, unscale=US%S_to_ppt)
   if (associated(tv%frazil)) &
-    call hchksum(tv%frazil, mesg//" tv%frazil", G%HI, scale=US%Q_to_J_kg*US%RZ_to_kg_m2)
+    call hchksum(tv%frazil, mesg//" tv%frazil", G%HI, unscale=US%Q_to_J_kg*US%RZ_to_kg_m2)
   if (associated(tv%salt_deficit)) &
-    call hchksum(tv%salt_deficit, mesg//" tv%salt_deficit", G%HI, scale=US%RZ_to_kg_m2*US%S_to_ppt)
+    call hchksum(tv%salt_deficit, mesg//" tv%salt_deficit", G%HI, unscale=US%RZ_to_kg_m2*US%S_to_ppt)
   if (associated(tv%TempxPmE)) &
-    call hchksum(tv%TempxPmE, mesg//" tv%TempxPmE", G%HI, scale=US%RZ_to_kg_m2*US%C_to_degC)
+    call hchksum(tv%TempxPmE, mesg//" tv%TempxPmE", G%HI, unscale=US%RZ_to_kg_m2*US%C_to_degC)
 end subroutine MOM_thermovar_chksum
 
 end module MOM_variables
diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index d5ba50e99a..b819c39ef1 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -1375,7 +1375,7 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv
     do j=js,je ; do i=is,ie
       work_2d(i,j) = G%mask2dT(i,j) * (ssh(i,j) + G%bathyT(i,j))
     enddo ; enddo
-    volo = global_area_integral(work_2d, G, scale=US%Z_to_m)
+    volo = global_area_integral(work_2d, G, unscale=US%Z_to_m)
     call post_data(IDs%id_volo, volo, diag)
   endif
 
diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90
index 719eb7d9e4..324808e374 100644
--- a/src/framework/MOM_horizontal_regridding.F90
+++ b/src/framework/MOM_horizontal_regridding.F90
@@ -612,7 +612,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr
     tr_prev(:,:) = tr_z(:,:,k)
 
     if (debug) then
-      call hchksum(tr_prev, 'field after fill ', G%HI, scale=I_scale)
+      call hchksum(tr_prev, 'field after fill ', G%HI, unscale=I_scale)
     endif
 
   enddo ! kd
@@ -921,7 +921,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, &
       call fill_miss_2d(tr_outf, good2, fill2, tr_prev, G, dtr_iter_stop, answer_date=ans_date)
 
 !     if (debug) then
-!       call hchksum(tr_outf, 'field from fill_miss_2d ', G%HI, scale=I_scale)
+!       call hchksum(tr_outf, 'field from fill_miss_2d ', G%HI, unscale=I_scale)
 !       call myStats(tr_outf, missing_value, G, k, 'field from fill_miss_2d()', unscale=I_scale)
 !     endif
 
@@ -930,7 +930,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, &
       tr_prev(:,:) = tr_z(:,:,k)
 
       if (debug) then
-        call hchksum(tr_prev, 'field after fill ', G%HI, scale=I_scale)
+        call hchksum(tr_prev, 'field after fill ', G%HI, unscale=I_scale)
       endif
 
     enddo ! kd
diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90
index a4b06b6150..eb021d8ffb 100644
--- a/src/ice_shelf/MOM_ice_shelf.F90
+++ b/src/ice_shelf/MOM_ice_shelf.F90
@@ -421,12 +421,12 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
 
   if (CS%debug) then
     call hchksum(fluxes_in%frac_shelf_h, "frac_shelf_h before apply melting", CS%Grid_in%HI, haloshift=0)
-    call hchksum(sfc_state_in%sst, "sst before apply melting", CS%Grid_in%HI, haloshift=0, scale=US%C_to_degC)
-    call hchksum(sfc_state_in%sss, "sss before apply melting", CS%Grid_in%HI, haloshift=0, scale=US%S_to_ppt)
+    call hchksum(sfc_state_in%sst, "sst before apply melting", CS%Grid_in%HI, haloshift=0, unscale=US%C_to_degC)
+    call hchksum(sfc_state_in%sss, "sss before apply melting", CS%Grid_in%HI, haloshift=0, unscale=US%S_to_ppt)
     call uvchksum("[uv]_ml before apply melting", sfc_state_in%u, sfc_state_in%v, &
-                  CS%Grid_in%HI, haloshift=0, scale=US%L_T_to_m_s)
+                  CS%Grid_in%HI, haloshift=0, unscale=US%L_T_to_m_s)
     call hchksum(sfc_state_in%ocean_mass, "ocean_mass before apply melting", CS%Grid_in%HI, haloshift=0, &
-                 scale=US%RZ_to_kg_m2)
+                 unscale=US%RZ_to_kg_m2)
   endif
 
   ! Calculate the friction velocity under ice shelves, using taux_shelf and tauy_shelf if possible.
@@ -774,9 +774,9 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
     if (bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_bdott(is:ie,js:je)
 
     if (CS%debug) then
-      call hchksum(ISS%h_shelf, "h_shelf after change thickness using melt", G%HI, haloshift=0, scale=US%Z_to_m)
+      call hchksum(ISS%h_shelf, "h_shelf after change thickness using melt", G%HI, haloshift=0, unscale=US%Z_to_m)
       call hchksum(ISS%mass_shelf, "mass_shelf after change thickness using melt", G%HI, haloshift=0, &
-                   scale=US%RZ_to_kg_m2)
+                   unscale=US%RZ_to_kg_m2)
     endif
   endif
 
@@ -790,9 +790,9 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
     if (bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_bdott(is:ie,js:je)
 
     if (CS%debug) then
-      call hchksum(ISS%h_shelf, "h_shelf after change thickness using melt", G%HI, haloshift=0, scale=US%Z_to_m)
+      call hchksum(ISS%h_shelf, "h_shelf after change thickness using melt", G%HI, haloshift=0, unscale=US%Z_to_m)
       call hchksum(ISS%mass_shelf, "mass_shelf after change thickness using melt", G%HI, haloshift=0, &
-                   scale=US%RZ_to_kg_m2)
+                   unscale=US%RZ_to_kg_m2)
     endif
 
     if (smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je)
@@ -800,9 +800,9 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
     if (smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_adott(is:ie,js:je)
 
     if (CS%debug) then
-      call hchksum(ISS%h_shelf, "h_shelf after change thickness using surf acc", G%HI, haloshift=0, scale=US%Z_to_m)
+      call hchksum(ISS%h_shelf, "h_shelf after change thickness using surf acc", G%HI, haloshift=0, unscale=US%Z_to_m)
       call hchksum(ISS%mass_shelf, "mass_shelf after change thickness using surf acc", G%HI, haloshift=0, &
-                   scale=US%RZ_to_kg_m2)
+                   unscale=US%RZ_to_kg_m2)
     endif
 
     update_ice_vel = .false.
@@ -1119,7 +1119,7 @@ subroutine add_shelf_forces(Ocn_grid, US, CS, forces, do_shelf_area, external_ca
   if (CS%debug) then
     call uvchksum("rigidity_ice_[uv]", forces%rigidity_ice_u, &
         forces%rigidity_ice_v, CS%Grid%HI, symmetric=.true., &
-        scale=US%L_to_m**3*US%L_to_Z*US%s_to_T, scalar_pair=.true.)
+        unscale=US%L_to_m**3*US%L_to_Z*US%s_to_T, scalar_pair=.true.)
     call uvchksum("frac_shelf_[uv]", forces%frac_shelf_u, &
         forces%frac_shelf_v, CS%Grid%HI, symmetric=.true., &
         scalar_pair=.true.)
@@ -1218,7 +1218,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes, time_step)
   if (CS%debug) then
     if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then
       call uvchksum("tau[xy]_shelf", sfc_state%taux_shelf, sfc_state%tauy_shelf, &
-                    G%HI, haloshift=0, scale=US%RZ_T_to_kg_m2s*US%L_T_to_m_s)
+                    G%HI, haloshift=0, unscale=US%RZ_T_to_kg_m2s*US%L_T_to_m_s)
     endif
   endif
 
@@ -1263,8 +1263,8 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes, time_step)
   endif ; enddo ; enddo
 
   if (CS%debug) then
-    call hchksum(ISS%water_flux, "water_flux add shelf fluxes", G%HI, haloshift=0, scale=US%RZ_T_to_kg_m2s)
-    call hchksum(ISS%tflux_ocn, "tflux_ocn add shelf fluxes", G%HI, haloshift=0, scale=US%QRZ_T_to_W_m2)
+    call hchksum(ISS%water_flux, "water_flux add shelf fluxes", G%HI, haloshift=0, unscale=US%RZ_T_to_kg_m2s)
+    call hchksum(ISS%tflux_ocn, "tflux_ocn add shelf fluxes", G%HI, haloshift=0, unscale=US%QRZ_T_to_W_m2)
     call MOM_forcing_chksum("After adding shelf fluxes", fluxes, G, CS%US, haloshift=0)
   endif
 
@@ -1920,8 +1920,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init,
       endif
     enddo ; enddo
     if (CS%debug) then
-      call hchksum(ISS%mass_shelf, "IS init: mass_shelf", G%HI, haloshift=0, scale=US%RZ_to_kg_m2)
-      call hchksum(ISS%area_shelf_h, "IS init: area_shelf", G%HI, haloshift=0, scale=US%L_to_m*US%L_to_m)
+      call hchksum(ISS%mass_shelf, "IS init: mass_shelf", G%HI, haloshift=0, unscale=US%RZ_to_kg_m2)
+      call hchksum(ISS%area_shelf_h, "IS init: area_shelf", G%HI, haloshift=0, unscale=US%L_to_m*US%L_to_m)
       call hchksum(ISS%hmask, "IS init: hmask", G%HI, haloshift=0)
     endif
 
@@ -1956,7 +1956,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init,
   enddo ; enddo
 
   if (CS%debug) then
-    call hchksum(ISS%area_shelf_h, "IS init: area_shelf_h", G%HI, haloshift=0, scale=US%L_to_m*US%L_to_m)
+    call hchksum(ISS%area_shelf_h, "IS init: area_shelf_h", G%HI, haloshift=0, unscale=US%L_to_m*US%L_to_m)
   endif
 
   CS%Time = Time
diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90
index 837789716b..7816df32de 100644
--- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90
+++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90
@@ -1441,8 +1441,8 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i
                                ISS%hmask, conv_flag, iters, time, CS%Phi, CS%Phisub)
 
     if (CS%debug) then
-      call qchksum(u_shlf, "u shelf", G%HI, haloshift=2, scale=US%L_T_to_m_s)
-      call qchksum(v_shlf, "v shelf", G%HI, haloshift=2, scale=US%L_T_to_m_s)
+      call qchksum(u_shlf, "u shelf", G%HI, haloshift=2, unscale=US%L_T_to_m_s)
+      call qchksum(v_shlf, "v shelf", G%HI, haloshift=2, unscale=US%L_T_to_m_s)
     endif
 
     write(mesg,*) "ice_shelf_solve_outer: linear solve done in ",iters," iterations"
@@ -3790,7 +3790,7 @@ subroutine ice_shelf_temp(CS, ISS, G, US, time_step, melt_rate, Time)
   call pass_var(CS%tmask, G%domain, complete=.true.)
 
   if (CS%debug) then
-    call hchksum(CS%t_shelf, "temp after front", G%HI, haloshift=3, scale=US%C_to_degC)
+    call hchksum(CS%t_shelf, "temp after front", G%HI, haloshift=3, unscale=US%C_to_degC)
   endif
 
 end subroutine ice_shelf_temp
diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90
index 322abc6d5e..10236994e0 100644
--- a/src/initialization/MOM_fixed_initialization.F90
+++ b/src/initialization/MOM_fixed_initialization.F90
@@ -101,7 +101,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir)
   call open_boundary_impose_land_mask(OBC, G, G%areaCu, G%areaCv, US)
 
   if (debug) then
-    call hchksum(G%bathyT, 'MOM_initialize_fixed: depth ', G%HI, haloshift=1, scale=US%Z_to_m)
+    call hchksum(G%bathyT, 'MOM_initialize_fixed: depth ', G%HI, haloshift=1, unscale=US%Z_to_m)
     call hchksum(G%mask2dT, 'MOM_initialize_fixed: mask2dT ', G%HI)
     call uvchksum('MOM_initialize_fixed: mask2dC[uv]', G%mask2dCu, &
                   G%mask2dCv, G%HI)
@@ -157,9 +157,9 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir)
 !   Calculate the components of grad f (beta)
   call MOM_calculate_grad_Coriolis(G%dF_dx, G%dF_dy, G, US=US)
   if (debug) then
-    call qchksum(G%CoriolisBu, "MOM_initialize_fixed: f ", G%HI, scale=US%s_to_T)
-    call hchksum(G%dF_dx, "MOM_initialize_fixed: dF_dx ", G%HI, scale=US%m_to_L*US%s_to_T)
-    call hchksum(G%dF_dy, "MOM_initialize_fixed: dF_dy ", G%HI, scale=US%m_to_L*US%s_to_T)
+    call qchksum(G%CoriolisBu, "MOM_initialize_fixed: f ", G%HI, unscale=US%s_to_T)
+    call hchksum(G%dF_dx, "MOM_initialize_fixed: dF_dx ", G%HI, unscale=US%m_to_L*US%s_to_T)
+    call hchksum(G%dF_dy, "MOM_initialize_fixed: dF_dy ", G%HI, unscale=US%m_to_L*US%s_to_T)
   endif
 
   call initialize_grid_rotation_angle(G, PF)
diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90
index 8bea8fe6e9..ef78a896c3 100644
--- a/src/initialization/MOM_grid_initialize.F90
+++ b/src/initialization/MOM_grid_initialize.F90
@@ -128,28 +128,28 @@ subroutine grid_metrics_chksum(parent, G, US)
   halo = min(G%ied-G%iec, G%jed-G%jec, 1)
 
   call hchksum_pair(trim(parent)//': d[xy]T', G%dxT, G%dyT, G%HI, &
-                    haloshift=halo, scale=US%L_to_m, scalar_pair=.true.)
+                    haloshift=halo, unscale=US%L_to_m, scalar_pair=.true.)
 
-  call uvchksum(trim(parent)//': dxC[uv]', G%dxCu, G%dyCv, G%HI, haloshift=halo, scale=US%L_to_m)
+  call uvchksum(trim(parent)//': dxC[uv]', G%dxCu, G%dyCv, G%HI, haloshift=halo, unscale=US%L_to_m)
 
-  call uvchksum(trim(parent)//': dxC[uv]', G%dyCu, G%dxCv, G%HI, haloshift=halo, scale=US%L_to_m)
+  call uvchksum(trim(parent)//': dxC[uv]', G%dyCu, G%dxCv, G%HI, haloshift=halo, unscale=US%L_to_m)
 
-  call Bchksum_pair(trim(parent)//': dxB[uv]', G%dxBu, G%dyBu, G%HI, haloshift=halo, scale=US%L_to_m)
+  call Bchksum_pair(trim(parent)//': dxB[uv]', G%dxBu, G%dyBu, G%HI, haloshift=halo, unscale=US%L_to_m)
 
   call hchksum_pair(trim(parent)//': Id[xy]T', G%IdxT, G%IdyT, G%HI, &
-                    haloshift=halo, scale=US%m_to_L, scalar_pair=.true.)
+                    haloshift=halo, unscale=US%m_to_L, scalar_pair=.true.)
 
-  call uvchksum(trim(parent)//': Id[xy]C[uv]', G%IdxCu, G%IdyCv, G%HI, haloshift=halo, scale=US%m_to_L)
+  call uvchksum(trim(parent)//': Id[xy]C[uv]', G%IdxCu, G%IdyCv, G%HI, haloshift=halo, unscale=US%m_to_L)
 
-  call uvchksum(trim(parent)//': Id[xy]C[uv]', G%IdyCu, G%IdxCv, G%HI, haloshift=halo, scale=US%m_to_L)
+  call uvchksum(trim(parent)//': Id[xy]C[uv]', G%IdyCu, G%IdxCv, G%HI, haloshift=halo, unscale=US%m_to_L)
 
-  call Bchksum_pair(trim(parent)//': Id[xy]B[uv]', G%IdxBu, G%IdyBu, G%HI, haloshift=halo, scale=US%m_to_L)
+  call Bchksum_pair(trim(parent)//': Id[xy]B[uv]', G%IdxBu, G%IdyBu, G%HI, haloshift=halo, unscale=US%m_to_L)
 
-  call hchksum(G%areaT, trim(parent)//': areaT',G%HI, haloshift=halo, scale=US%L_to_m**2)
-  call Bchksum(G%areaBu, trim(parent)//': areaBu',G%HI, haloshift=halo, scale=US%L_to_m**2)
+  call hchksum(G%areaT, trim(parent)//': areaT',G%HI, haloshift=halo, unscale=US%L_to_m**2)
+  call Bchksum(G%areaBu, trim(parent)//': areaBu',G%HI, haloshift=halo, unscale=US%L_to_m**2)
 
-  call hchksum(G%IareaT, trim(parent)//': IareaT',G%HI, haloshift=halo, scale=US%m_to_L**2)
-  call Bchksum(G%IareaBu, trim(parent)//': IareaBu',G%HI, haloshift=halo, scale=US%m_to_L**2)
+  call hchksum(G%IareaT, trim(parent)//': IareaT',G%HI, haloshift=halo, unscale=US%m_to_L**2)
+  call Bchksum(G%IareaBu, trim(parent)//': IareaBu',G%HI, haloshift=halo, unscale=US%m_to_L**2)
 
   call hchksum(G%geoLonT,trim(parent)//': geoLonT',G%HI, haloshift=halo)
   call hchksum(G%geoLatT,trim(parent)//': geoLatT',G%HI, haloshift=halo)
diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90
index 857cad578d..1467cdaaad 100644
--- a/src/initialization/MOM_state_initialization.F90
+++ b/src/initialization/MOM_state_initialization.F90
@@ -450,7 +450,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
            "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True")
 
   if (new_sim .and. debug .and. (depress_sfc .or. trim_ic_for_p_surf)) &
-    call hchksum(h, "Pre-depress: h ", G%HI, haloshift=1, scale=GV%H_to_MKS)
+    call hchksum(h, "Pre-depress: h ", G%HI, haloshift=1, unscale=GV%H_to_MKS)
 
   ! Remove the mass that would be displaced by an ice shelf or inverse barometer.
   if (depress_sfc) then
@@ -479,7 +479,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
                      units="s", scale=US%s_to_T, fail_if_missing=.true.)
 
       if (new_sim .and. debug) &
-        call hchksum(h, "Pre-ALE_regrid: h ", G%HI, haloshift=1, scale=GV%H_to_MKS)
+        call hchksum(h, "Pre-ALE_regrid: h ", G%HI, haloshift=1, unscale=GV%H_to_MKS)
       call ALE_regrid_accelerated(ALE_CSp, G, GV, US, h, tv, regrid_iterations, u, v, OBC, tracer_Reg, &
                                   dt=dt, initial=.true.)
     endif
@@ -517,7 +517,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
 
   if (new_sim) call pass_vector(u, v, G%Domain)
   if (debug .and. new_sim) then
-    call uvchksum("MOM_initialize_state [uv]", u, v, G%HI, haloshift=1, scale=US%L_T_to_m_s)
+    call uvchksum("MOM_initialize_state [uv]", u, v, G%HI, haloshift=1, unscale=US%L_T_to_m_s)
   endif
 
   ! This is the end of the block of code that might have initialized fields
@@ -552,14 +552,14 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
   call pass_var(h, G%Domain)
 
   if (debug) then
-    call hchksum(h, "MOM_initialize_state: h ", G%HI, haloshift=1, scale=GV%H_to_MKS)
-    if ( use_temperature ) call hchksum(tv%T, "MOM_initialize_state: T ", G%HI, haloshift=1, scale=US%C_to_degC)
-    if ( use_temperature ) call hchksum(tv%S, "MOM_initialize_state: S ", G%HI, haloshift=1, scale=US%S_to_ppt)
+    call hchksum(h, "MOM_initialize_state: h ", G%HI, haloshift=1, unscale=GV%H_to_MKS)
+    if ( use_temperature ) call hchksum(tv%T, "MOM_initialize_state: T ", G%HI, haloshift=1, unscale=US%C_to_degC)
+    if ( use_temperature ) call hchksum(tv%S, "MOM_initialize_state: S ", G%HI, haloshift=1, unscale=US%S_to_ppt)
     if ( use_temperature .and. debug_layers) then ; do k=1,nz
       write(mesg,'("MOM_IS: T[",I2,"]")') k
-      call hchksum(tv%T(:,:,k), mesg, G%HI, haloshift=1, scale=US%C_to_degC)
+      call hchksum(tv%T(:,:,k), mesg, G%HI, haloshift=1, unscale=US%C_to_degC)
       write(mesg,'("MOM_IS: S[",I2,"]")') k
-      call hchksum(tv%S(:,:,k), mesg, G%HI, haloshift=1, scale=US%S_to_ppt)
+      call hchksum(tv%S(:,:,k), mesg, G%HI, haloshift=1, unscale=US%S_to_ppt)
     enddo ; endif
   endif
 
@@ -2310,7 +2310,7 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p
             call MOM_error(FATAL, " initialize_oda_incupd_uv: Unable to open "//trim(filename))
     allocate(tmp_u(G%IsdB:G%IedB,jsd:jed,nz_data), source=0.0)
     allocate(tmp_v(isd:ied,G%JsdB:G%JedB,nz_data), source=0.0)
-    call MOM_read_vector(filename, uinc_var, vinc_var, tmp_u, tmp_v, G%Domain,scale=US%m_s_to_L_T)
+    call MOM_read_vector(filename, uinc_var, vinc_var, tmp_u, tmp_v, G%Domain, scale=US%m_s_to_L_T)
     call set_up_oda_incupd_vel_field(tmp_u, tmp_v, G, GV, oda_incupd_CSp)
     deallocate(tmp_u, tmp_v)
   endif
diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90
index bfa79ba053..f3bd9e8934 100644
--- a/src/parameterizations/lateral/MOM_MEKE.F90
+++ b/src/parameterizations/lateral/MOM_MEKE.F90
@@ -253,17 +253,17 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
   case(EKE_PROG)
     if (CS%debug) then
       if (allocated(MEKE%mom_src)) &
-        call hchksum(MEKE%mom_src, 'MEKE mom_src', G%HI, scale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
+        call hchksum(MEKE%mom_src, 'MEKE mom_src', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
       if (allocated(MEKE%GME_snk)) &
-        call hchksum(MEKE%GME_snk, 'MEKE GME_snk', G%HI, scale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
+        call hchksum(MEKE%GME_snk, 'MEKE GME_snk', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
       if (allocated(MEKE%GM_src)) &
-        call hchksum(MEKE%GM_src, 'MEKE GM_src', G%HI, scale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
+        call hchksum(MEKE%GM_src, 'MEKE GM_src', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
       if (allocated(MEKE%MEKE)) &
-        call hchksum(MEKE%MEKE, 'MEKE MEKE', G%HI, scale=US%L_T_to_m_s**2)
-      call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI, scale=US%s_to_T, &
+        call hchksum(MEKE%MEKE, 'MEKE MEKE', G%HI, unscale=US%L_T_to_m_s**2)
+      call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI, unscale=US%s_to_T, &
                     scalar_pair=.true.)
       call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=0, symmetric=.true., &
-                    scale=GV%H_to_m*(US%L_to_m**2))
+                    unscale=GV%H_to_m*US%L_to_m**2)
     endif
 
     sdt = dt*CS%MEKE_dtScale ! Scaled dt to use for time-stepping
@@ -375,12 +375,12 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
     if (CS%debug) then
       if (CS%visc_drag) &
         call uvchksum("MEKE drag_vel_[uv]", drag_vel_u, drag_vel_v, G%HI, &
-                      scale=GV%H_to_mks*US%s_to_T, scalar_pair=.true.)
-      call hchksum(mass, 'MEKE mass',G%HI,haloshift=1, scale=US%RZ_to_kg_m2)
-      call hchksum(drag_rate_visc, 'MEKE drag_rate_visc', G%HI, scale=GV%H_to_mks*US%s_to_T)
+                      unscale=GV%H_to_mks*US%s_to_T, scalar_pair=.true.)
+      call hchksum(mass, 'MEKE mass',G%HI,haloshift=1, unscale=US%RZ_to_kg_m2)
+      call hchksum(drag_rate_visc, 'MEKE drag_rate_visc', G%HI, unscale=GV%H_to_mks*US%s_to_T)
       call hchksum(bottomFac2, 'MEKE bottomFac2', G%HI)
       call hchksum(barotrFac2, 'MEKE barotrFac2', G%HI)
-      call hchksum(LmixScale, 'MEKE LmixScale', G%HI,scale=US%L_to_m)
+      call hchksum(LmixScale, 'MEKE LmixScale', G%HI, unscale=US%L_to_m)
     endif
 
     ! Aggregate sources of MEKE (background, frictional and GM)
@@ -427,7 +427,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
     endif
 
     if (CS%debug) then
-      call hchksum(src, "MEKE src", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T**3)
+      call hchksum(src, "MEKE src", G%HI, haloshift=0, unscale=US%L_to_m**2*US%s_to_T**3)
     endif
 
     ! Increase EKE by a full time-steps worth of source
@@ -630,7 +630,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
     endif ! MEKE_KH>=0
 
     if (CS%debug) then
-      call hchksum(MEKE%MEKE, "MEKE post-update MEKE", G%HI, haloshift=0, scale=US%L_T_to_m_s**2)
+      call hchksum(MEKE%MEKE, "MEKE post-update MEKE", G%HI, haloshift=0, unscale=US%L_T_to_m_s**2)
     endif
 
   case(EKE_FILE)
diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index b9e50f6b80..2a4181804b 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -594,7 +594,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
     call pass_vector(KH_u_GME, KH_v_GME, G%domain, To_All+Scalar_Pair)
 
     if (CS%debug) &
-    call uvchksum("GME KH[u,v]_GME", KH_u_GME, KH_v_GME, G%HI, haloshift=2, scale=US%L_to_m**2*US%s_to_T)
+      call uvchksum("GME KH[u,v]_GME", KH_u_GME, KH_v_GME, G%HI, haloshift=2, unscale=US%L_to_m**2*US%s_to_T)
 
   endif ! use_GME
 
@@ -1996,11 +1996,11 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
 
   if (CS%debug) then
     if (CS%Laplacian) then
-      call hchksum(Kh_h, "Kh_h", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T)
-      call Bchksum(Kh_q, "Kh_q", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T)
+      call hchksum(Kh_h, "Kh_h", G%HI, haloshift=0, unscale=US%L_to_m**2*US%s_to_T)
+      call Bchksum(Kh_q, "Kh_q", G%HI, haloshift=0, unscale=US%L_to_m**2*US%s_to_T)
     endif
-    if (CS%biharmonic) call hchksum(Ah_h, "Ah_h", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T)
-    if (CS%biharmonic) call Bchksum(Ah_q, "Ah_q", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T)
+    if (CS%biharmonic) call hchksum(Ah_h, "Ah_h", G%HI, haloshift=0, unscale=US%L_to_m**4*US%s_to_T)
+    if (CS%biharmonic) call Bchksum(Ah_q, "Ah_q", G%HI, haloshift=0, unscale=US%L_to_m**4*US%s_to_T)
   endif
 
   if (CS%id_FrictWorkIntz > 0) then
@@ -2722,8 +2722,8 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
         CS%Kh_Max_xy(I,J) = CS%bound_coef * 0.25 * Idt / denom
     enddo ; enddo
     if (CS%debug) then
-      call hchksum(CS%Kh_Max_xx, "Kh_Max_xx", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T)
-      call Bchksum(CS%Kh_Max_xy, "Kh_Max_xy", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T)
+      call hchksum(CS%Kh_Max_xx, "Kh_Max_xx", G%HI, haloshift=0, unscale=US%L_to_m**2*US%s_to_T)
+      call Bchksum(CS%Kh_Max_xy, "Kh_Max_xy", G%HI, haloshift=0, unscale=US%L_to_m**2*US%s_to_T)
     endif
   endif
   ! The biharmonic bounds should avoid overshoots when CS%bound_coef < 0.5, but
@@ -2778,8 +2778,8 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
         CS%Ah_Max_xy(I,J) = CS%bound_coef * 0.5 * Idt / denom
     enddo ; enddo
     if (CS%debug) then
-      call hchksum(CS%Ah_Max_xx, "Ah_Max_xx", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T)
-      call Bchksum(CS%Ah_Max_xy, "Ah_Max_xy", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T)
+      call hchksum(CS%Ah_Max_xx, "Ah_Max_xx", G%HI, haloshift=0, unscale=US%L_to_m**4*US%s_to_T)
+      call Bchksum(CS%Ah_Max_xy, "Ah_Max_xy", G%HI, haloshift=0, unscale=US%L_to_m**4*US%s_to_T)
     endif
   endif
   ! Register fields for output from this module.
diff --git a/src/parameterizations/lateral/MOM_interface_filter.F90 b/src/parameterizations/lateral/MOM_interface_filter.F90
index 07b698e294..ea86b80594 100644
--- a/src/parameterizations/lateral/MOM_interface_filter.F90
+++ b/src/parameterizations/lateral/MOM_interface_filter.F90
@@ -142,9 +142,9 @@ subroutine interface_filter(h, uhtr, vhtr, tv, dt, G, GV, US, CDp, CS)
 
   if (CS%debug) then
     call uvchksum("Kh_[uv]", Lsm2_u, Lsm2_v, G%HI, haloshift=hs, &
-                  scale=US%L_to_m**2, scalar_pair=.true.)
-    call hchksum(h, "interface_filter_1 h", G%HI, haloshift=hs+1, scale=GV%H_to_m)
-    call hchksum(e, "interface_filter_1 e", G%HI, haloshift=hs+1, scale=US%Z_to_m)
+                  unscale=US%L_to_m**2, scalar_pair=.true.)
+    call hchksum(h, "interface_filter_1 h", G%HI, haloshift=hs+1, unscale=GV%H_to_m)
+    call hchksum(e, "interface_filter_1 e", G%HI, haloshift=hs+1, unscale=US%Z_to_m)
   endif
 
   ! Calculate uhD, vhD from h, e, Lsm2_u, Lsm2_v
@@ -225,10 +225,10 @@ subroutine interface_filter(h, uhtr, vhtr, tv, dt, G, GV, US, CDp, CS)
 
   if (CS%debug) then
     call uvchksum("interface_filter [uv]hD", uhD, vhD, &
-                  G%HI, haloshift=0, scale=GV%H_to_m*US%L_to_m**2)
+                  G%HI, haloshift=0, unscale=GV%H_to_m*US%L_to_m**2)
     call uvchksum("interface_filter [uv]htr", uhtr, vhtr, &
-                  G%HI, haloshift=0, scale=US%L_to_m**2*GV%H_to_m)
-    call hchksum(h, "interface_filter h", G%HI, haloshift=0, scale=GV%H_to_m)
+                  G%HI, haloshift=0, unscale=US%L_to_m**2*GV%H_to_m)
+    call hchksum(h, "interface_filter h", G%HI, haloshift=0, unscale=GV%H_to_m)
   endif
 
 end subroutine interface_filter
diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90
index 3f67e0e2b5..8b3a2c2dbe 100644
--- a/src/parameterizations/lateral/MOM_internal_tides.F90
+++ b/src/parameterizations/lateral/MOM_internal_tides.F90
@@ -906,7 +906,7 @@ subroutine sum_En(G, US, CS, En, label)
 
   En_sum = 0.0
   do a=1,CS%nAngle
-    En_sum = En_sum + global_area_integral(En(:,:,a), G, scale=US%RZ3_T3_to_W_m2*US%T_to_s)
+    En_sum = En_sum + global_area_integral(En(:,:,a), G, unscale=US%RZ3_T3_to_W_m2*US%T_to_s)
   enddo
   CS%En_sum = En_sum
   !En_sum_diff = En_sum - CS%En_sum
diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
index 7db0dc4e90..7315c82601 100644
--- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
+++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
@@ -455,9 +455,9 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS)
   endif
 
   if (CS%debug) then
-    call hchksum(CS%cg1, "calc_resoln_fn cg1", G%HI, haloshift=1, scale=US%L_T_to_m_s)
+    call hchksum(CS%cg1, "calc_resoln_fn cg1", G%HI, haloshift=1, unscale=US%L_T_to_m_s)
     call uvchksum("Res_fn_[uv]", CS%Res_fn_u, CS%Res_fn_v, G%HI, haloshift=0, &
-                  scale=1.0, scalar_pair=.true.)
+                  unscale=1.0, scalar_pair=.true.)
   endif
 
 end subroutine calc_resoln_function
@@ -656,11 +656,11 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C
 
   if (CS%debug) then
     call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, G%HI, &
-                  scale=US%Z_to_L, haloshift=1)
+                  unscale=US%Z_to_L, haloshift=1)
     call uvchksum("calc_Visbeck_coeffs_old N2_u, N2_v", N2_u, N2_v, G%HI, &
-                  scale=US%L_to_Z**2*US%s_to_T**2, scalar_pair=.true.)
+                  unscale=US%L_to_Z**2*US%s_to_T**2, scalar_pair=.true.)
     call uvchksum("calc_Visbeck_coeffs_old SN_[uv]", CS%SN_u, CS%SN_v, G%HI, &
-                  scale=US%s_to_T, scalar_pair=.true.)
+                  unscale=US%s_to_T, scalar_pair=.true.)
   endif
 
 end subroutine calc_Visbeck_coeffs_old
@@ -700,9 +700,9 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN,
   crop = CS%cropping_distance>=0. ! Only filter out in-/out-cropped interface is parameter if non-negative
 
   if (CS%debug) then
-    call uvchksum("calc_Eady_growth_rate_2D dz[uv]", dzu, dzv, G%HI, scale=US%Z_to_m, scalar_pair=.true.)
+    call uvchksum("calc_Eady_growth_rate_2D dz[uv]", dzu, dzv, G%HI, unscale=US%Z_to_m, scalar_pair=.true.)
     call uvchksum("calc_Eady_growth_rate_2D dzS2N2[uv]", dzSxN, dzSyN, G%HI, &
-                  scale=US%Z_to_m*US%s_to_T, scalar_pair=.true.)
+                  unscale=US%Z_to_m*US%s_to_T, scalar_pair=.true.)
   endif
 
   !$OMP parallel do default(shared)
@@ -813,7 +813,7 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN,
 
   if (CS%debug) then
     call uvchksum("calc_Eady_growth_rate_2D SN_[uv]", CS%SN_u, CS%SN_v, G%HI, &
-                  scale=US%s_to_T, scalar_pair=.true.)
+                  unscale=US%s_to_T, scalar_pair=.true.)
   endif
 
 end subroutine calc_Eady_growth_rate_2D
diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
index 4361712dde..d1e0d15293 100644
--- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
+++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
@@ -345,8 +345,8 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix,
   ! Apply time filter (to remove diurnal cycle)
   if (CS%MLE_MLD_decay_time>0.) then
     if (CS%debug) then
-      call hchksum(CS%MLD_filtered, 'mixed_layer_restrat: MLD_filtered', G%HI, haloshift=1, scale=GV%H_to_mks)
-      call hchksum(h_MLD, 'mixed_layer_restrat: MLD in', G%HI, haloshift=1, scale=GV%H_to_mks)
+      call hchksum(CS%MLD_filtered, 'mixed_layer_restrat: MLD_filtered', G%HI, haloshift=1, unscale=GV%H_to_mks)
+      call hchksum(h_MLD, 'mixed_layer_restrat: MLD in', G%HI, haloshift=1, unscale=GV%H_to_mks)
     endif
     aFac = CS%MLE_MLD_decay_time / ( dt + CS%MLE_MLD_decay_time )
     bFac = dt / ( dt + CS%MLE_MLD_decay_time )
@@ -362,8 +362,9 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix,
   ! Apply slower time filter (to remove seasonal cycle) on already filtered MLD_fast
   if (CS%MLE_MLD_decay_time2>0.) then
     if (CS%debug) then
-      call hchksum(CS%MLD_filtered_slow, 'mixed_layer_restrat: MLD_filtered_slow', G%HI, haloshift=1, scale=GV%H_to_mks)
-      call hchksum(MLD_fast, 'mixed_layer_restrat: MLD fast', G%HI, haloshift=1, scale=GV%H_to_mks)
+      call hchksum(CS%MLD_filtered_slow, 'mixed_layer_restrat: MLD_filtered_slow', G%HI, &
+                   haloshift=1, unscale=GV%H_to_mks)
+      call hchksum(MLD_fast, 'mixed_layer_restrat: MLD fast', G%HI, haloshift=1, unscale=GV%H_to_mks)
     endif
     aFac = CS%MLE_MLD_decay_time2 / ( dt + CS%MLE_MLD_decay_time2 )
     bFac = dt / ( dt + CS%MLE_MLD_decay_time2 )
@@ -490,11 +491,11 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix,
   endif
 
   if (CS%debug) then
-    call hchksum(h, 'mixed_layer_restrat: h', G%HI, haloshift=1, scale=GV%H_to_mks)
-    call hchksum(U_star_2d, 'mixed_layer_restrat: u*', G%HI, haloshift=1, scale=GV%H_to_m*US%s_to_T)
-    call hchksum(MLD_fast, 'mixed_layer_restrat: MLD', G%HI, haloshift=1, scale=GV%H_to_mks)
+    call hchksum(h, 'mixed_layer_restrat: h', G%HI, haloshift=1, unscale=GV%H_to_mks)
+    call hchksum(U_star_2d, 'mixed_layer_restrat: u*', G%HI, haloshift=1, unscale=GV%H_to_m*US%s_to_T)
+    call hchksum(MLD_fast, 'mixed_layer_restrat: MLD', G%HI, haloshift=1, unscale=GV%H_to_mks)
     call hchksum(Rml_av_fast, 'mixed_layer_restrat: rml', G%HI, haloshift=1, &
-                 scale=GV%m_to_H*US%L_T_to_m_s**2)
+                 unscale=GV%m_to_H*US%L_T_to_m_s**2)
   endif
 
 ! TO DO:
@@ -871,16 +872,16 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
   call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1)
 
   if (CS%debug) then
-    call hchksum(h,'mixed_Bodner: h', G%HI, haloshift=1, scale=GV%H_to_mks)
-    call hchksum(BLD, 'mle_Bodner: BLD', G%HI, haloshift=1, scale=US%Z_to_m)
-    call hchksum(h_MLD, 'mle_Bodner: h_MLD', G%HI, haloshift=1, scale=GV%H_to_mks)
+    call hchksum(h,'mixed_Bodner: h', G%HI, haloshift=1, unscale=GV%H_to_mks)
+    call hchksum(BLD, 'mle_Bodner: BLD', G%HI, haloshift=1, unscale=US%Z_to_m)
+    call hchksum(h_MLD, 'mle_Bodner: h_MLD', G%HI, haloshift=1, unscale=GV%H_to_mks)
     if (associated(bflux)) &
-      call hchksum(bflux, 'mle_Bodner: bflux', G%HI, haloshift=1, scale=US%Z_to_m**2*US%s_to_T**3)
-    call hchksum(U_star_2d, 'mle_Bodner: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T)
+      call hchksum(bflux, 'mle_Bodner: bflux', G%HI, haloshift=1, unscale=US%Z_to_m**2*US%s_to_T**3)
+    call hchksum(U_star_2d, 'mle_Bodner: u*', G%HI, haloshift=1, unscale=US%Z_to_m*US%s_to_T)
     call hchksum(CS%MLD_filtered, 'mle_Bodner: MLD_filtered 1', &
-                 G%HI, haloshift=1, scale=GV%H_to_mks)
+                 G%HI, haloshift=1, unscale=GV%H_to_mks)
     call hchksum(CS%MLD_filtered_slow,'mle_Bodner: MLD_filtered_slow 1', &
-                 G%HI, haloshift=1, scale=GV%H_to_mks)
+                 G%HI, haloshift=1, unscale=GV%H_to_mks)
   endif
 
   ! Apply time filter to h_MLD (to remove diurnal cycle) to obtain "little h".
@@ -977,13 +978,13 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
   endif
 
   if (CS%debug) then
-    call hchksum(little_h,'mle_Bodner: little_h', G%HI, haloshift=1, scale=GV%H_to_mks)
-    call hchksum(big_H,'mle_Bodner: big_H', G%HI, haloshift=1, scale=GV%H_to_mks)
+    call hchksum(little_h,'mle_Bodner: little_h', G%HI, haloshift=1, unscale=GV%H_to_mks)
+    call hchksum(big_H,'mle_Bodner: big_H', G%HI, haloshift=1, unscale=GV%H_to_mks)
     call hchksum(CS%MLD_filtered,'mle_Bodner: MLD_filtered 2', &
-                 G%HI, haloshift=1, scale=GV%H_to_mks)
+                 G%HI, haloshift=1, unscale=GV%H_to_mks)
     call hchksum(CS%MLD_filtered_slow,'mle_Bodner: MLD_filtered_slow 2', &
-                 G%HI, haloshift=1, scale=GV%H_to_mks)
-    call hchksum(wpup,'mle_Bodner: wpup', G%HI, haloshift=1, scale=US%L_to_m*GV%H_to_mks*US%s_to_T**2)
+                 G%HI, haloshift=1, unscale=GV%H_to_mks)
+    call hchksum(wpup,'mle_Bodner: wpup', G%HI, haloshift=1, unscale=US%L_to_m*GV%H_to_mks*US%s_to_T**2)
   endif
 
   ! Calculate the average density in the "mixed layer".
@@ -1047,10 +1048,10 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
   enddo
 
   if (CS%debug) then
-    call hchksum(htot,'mle_Bodner: htot', G%HI, haloshift=1, scale=GV%H_to_mks)
+    call hchksum(htot,'mle_Bodner: htot', G%HI, haloshift=1, unscale=GV%H_to_mks)
     call hchksum(vol_dt_avail,'mle_Bodner: vol_dt_avail', G%HI, haloshift=1, &
-                 scale=US%L_to_m**2*GV%H_to_mks*US%s_to_T)
-    call hchksum(buoy_av,'mle_Bodner: buoy_av', G%HI, haloshift=1, scale=GV%m_to_H*US%L_T_to_m_s**2)
+                 unscale=US%L_to_m**2*GV%H_to_mks*US%s_to_T)
+    call hchksum(buoy_av,'mle_Bodner: buoy_av', G%HI, haloshift=1, unscale=GV%m_to_H*US%L_T_to_m_s**2)
   endif
 
   ! U - Component
diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90
index e9d9b7f50d..d735248417 100644
--- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90
+++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90
@@ -463,23 +463,23 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
 
   if (CS%debug) then
     call uvchksum("Kh_[uv]", Kh_u, Kh_v, G%HI, haloshift=0, &
-                  scale=(US%L_to_m**2)*US%s_to_T, scalar_pair=.true.)
+                  unscale=(US%L_to_m**2)*US%s_to_T, scalar_pair=.true.)
     call uvchksum("Kh_[uv]_CFL", Kh_u_CFL, Kh_v_CFL, G%HI, haloshift=0, &
-                  scale=(US%L_to_m**2)*US%s_to_T, scalar_pair=.true.)
+                  unscale=(US%L_to_m**2)*US%s_to_T, scalar_pair=.true.)
     if (Resoln_scaled) then
       call uvchksum("Res_fn_[uv]", VarMix%Res_fn_u, VarMix%Res_fn_v, G%HI, haloshift=0, &
-                    scale=1.0, scalar_pair=.true.)
+                    unscale=1.0, scalar_pair=.true.)
     endif
     call uvchksum("int_slope_[uv]", int_slope_u, int_slope_v, G%HI, haloshift=0)
-    call hchksum(h, "thickness_diffuse_1 h", G%HI, haloshift=1, scale=GV%H_to_m)
-    call hchksum(e, "thickness_diffuse_1 e", G%HI, haloshift=1, scale=US%Z_to_m)
+    call hchksum(h, "thickness_diffuse_1 h", G%HI, haloshift=1, unscale=GV%H_to_m)
+    call hchksum(e, "thickness_diffuse_1 e", G%HI, haloshift=1, unscale=US%Z_to_m)
     if (use_stored_slopes) then
       call uvchksum("VarMix%slope_[xy]", VarMix%slope_x, VarMix%slope_y, &
-                    G%HI, haloshift=0, scale=US%Z_to_L)
+                    G%HI, haloshift=0, unscale=US%Z_to_L)
     endif
     if (associated(tv%eqn_of_state)) then
-      call hchksum(tv%T, "thickness_diffuse T", G%HI, haloshift=1, scale=US%C_to_degC)
-      call hchksum(tv%S, "thickness_diffuse S", G%HI, haloshift=1, scale=US%S_to_ppt)
+      call hchksum(tv%T, "thickness_diffuse T", G%HI, haloshift=1, unscale=US%C_to_degC)
+      call hchksum(tv%S, "thickness_diffuse S", G%HI, haloshift=1, unscale=US%S_to_ppt)
     endif
   endif
 
@@ -588,10 +588,10 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
 
   if (CS%debug) then
     call uvchksum("thickness_diffuse [uv]hD", uhD, vhD, &
-                  G%HI, haloshift=0, scale=GV%H_to_m*US%L_to_m**2*US%s_to_T)
+                  G%HI, haloshift=0, unscale=GV%H_to_m*US%L_to_m**2*US%s_to_T)
     call uvchksum("thickness_diffuse [uv]htr", uhtr, vhtr, &
-                  G%HI, haloshift=0, scale=US%L_to_m**2*GV%H_to_m)
-    call hchksum(h, "thickness_diffuse h", G%HI, haloshift=0, scale=GV%H_to_m)
+                  G%HI, haloshift=0, unscale=US%L_to_m**2*GV%H_to_m)
+    call hchksum(h, "thickness_diffuse h", G%HI, haloshift=0, unscale=GV%H_to_m)
   endif
 
 end subroutine thickness_diffuse
@@ -1589,11 +1589,11 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
       enddo ; enddo ; enddo
     endif
     if (CS%debug) then
-      call hchksum(MEKE%GM_src, 'MEKE%GM_src', G%HI, scale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
-      call uvchksum("KH_[uv]", Kh_u, Kh_v, G%HI, scale=US%L_to_m**2*US%s_to_T, &
+      call hchksum(MEKE%GM_src, 'MEKE%GM_src', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
+      call uvchksum("KH_[uv]", Kh_u, Kh_v, G%HI, unscale=US%L_to_m**2*US%s_to_T, &
                     scalar_pair=.true.)
-      call uvchksum("Slope_[xy]_PE", Slope_x_PE, Slope_y_PE, G%HI, scale=US%Z_to_L)
-      call uvchksum("hN2_[xy]_PE", hN2_x_PE, hN2_y_PE, G%HI, scale=GV%H_to_mks*US%L_to_Z**2*US%s_to_T**2, &
+      call uvchksum("Slope_[xy]_PE", Slope_x_PE, Slope_y_PE, G%HI, unscale=US%Z_to_L)
+      call uvchksum("hN2_[xy]_PE", hN2_x_PE, hN2_y_PE, G%HI, unscale=GV%H_to_mks*US%L_to_Z**2*US%s_to_T**2, &
                     scalar_pair=.true.)
     endif
   endif ; endif
diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90
index f480c655d7..200cb02443 100644
--- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90
+++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90
@@ -664,11 +664,11 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, &
       "KPP_calculate: The Waves control structure must be associated if STOKES_MIXING is True.")
 
   if (CS%debug) then
-    call hchksum(h, "KPP in: h",G%HI,haloshift=0, scale=GV%H_to_m)
-    call hchksum(uStar, "KPP in: uStar",G%HI,haloshift=0, scale=US%Z_to_m*US%s_to_T)
-    call hchksum(buoyFlux, "KPP in: buoyFlux",G%HI,haloshift=0, scale=US%L_to_m**2*US%s_to_T**3)
-    call hchksum(Kt, "KPP in: Kt",G%HI,haloshift=0, scale=GV%HZ_T_to_m2_s)
-    call hchksum(Ks, "KPP in: Ks",G%HI,haloshift=0, scale=GV%HZ_T_to_m2_s)
+    call hchksum(h, "KPP in: h", G%HI, haloshift=0, unscale=GV%H_to_m)
+    call hchksum(uStar, "KPP in: uStar", G%HI, haloshift=0, unscale=US%Z_to_m*US%s_to_T)
+    call hchksum(buoyFlux, "KPP in: buoyFlux", G%HI, haloshift=0, unscale=US%L_to_m**2*US%s_to_T**3)
+    call hchksum(Kt, "KPP in: Kt", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
+    call hchksum(Ks, "KPP in: Ks", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
   endif
 
   nonLocalTrans(:,:) = 0.0
@@ -920,8 +920,8 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, &
   call cpu_clock_end(id_clock_KPP_calc)
 
   if (CS%debug) then
-    call hchksum(Kt, "KPP out: Kt", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
-    call hchksum(Ks, "KPP out: Ks", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+    call hchksum(Kt, "KPP out: Kt", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
+    call hchksum(Ks, "KPP out: Ks", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
   endif
 
   ! send diagnostics to post_data
@@ -1022,10 +1022,10 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
       "KPP_compute_BLD: The Waves control structure must be associated if STOKES_MIXING is True.")
 
   if (CS%debug) then
-    call hchksum(Salt, "KPP in: S", G%HI, haloshift=0, scale=US%S_to_ppt)
-    call hchksum(Temp, "KPP in: T", G%HI, haloshift=0, scale=US%C_to_degC)
-    call hchksum(u, "KPP in: u",G%HI,haloshift=0,scale=US%L_T_to_m_s)
-    call hchksum(v, "KPP in: v",G%HI,haloshift=0,scale=US%L_T_to_m_s)
+    call hchksum(Salt, "KPP in: S", G%HI, haloshift=0, unscale=US%S_to_ppt)
+    call hchksum(Temp, "KPP in: T", G%HI, haloshift=0, unscale=US%C_to_degC)
+    call hchksum(u, "KPP in: u", G%HI, haloshift=0, unscale=US%L_T_to_m_s)
+    call hchksum(v, "KPP in: v", G%HI, haloshift=0, unscale=US%L_T_to_m_s)
   endif
 
   call cpu_clock_begin(id_clock_KPP_compute_BLD)
diff --git a/src/parameterizations/vertical/MOM_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90
index 19744cb6c5..ce592a2d9c 100644
--- a/src/parameterizations/vertical/MOM_CVMix_conv.F90
+++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90
@@ -281,13 +281,13 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux)
   enddo
 
   if (CS%debug) then
-    ! if (CS%id_N2 > 0) call hchksum(N2_3d, "MOM_CVMix_conv: N2",G%HI,haloshift=0)
+    ! if (CS%id_N2 > 0) call hchksum(N2_3d, "MOM_CVMix_conv: N2", G%HI, haloshift=0, unscale=US%s_to_T**2)
     ! if (CS%id_kd_conv > 0) &
-    !   call hchksum(Kd_conv, "MOM_CVMix_conv: Kd_conv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
+    !   call hchksum(Kd_conv, "MOM_CVMix_conv: Kd_conv", G%HI, haloshift=0, unscale=US%Z2_T_to_m2_s)
     ! if (CS%id_kv_conv > 0) &
-    !   call hchksum(Kv_conv, "MOM_CVMix_conv: Kv_conv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
-    call hchksum(Kd, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
-    call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+    !   call hchksum(Kv_conv, "MOM_CVMix_conv: Kv_conv", G%HI, haloshift=0, unscale=US%Z2_T_to_m2_s)
+    call hchksum(Kd, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
+    call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
   endif
 
   ! send diagnostics to post_data
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index 89ed3bada9..6d8afa1493 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -648,7 +648,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
     call MOM_state_chksum("after set_diffusivity ", u, v, h, G, GV, US, haloshift=0)
     call MOM_forcing_chksum("after set_diffusivity ", fluxes, G, US, haloshift=0)
     call MOM_thermovar_chksum("after set_diffusivity ", tv, G, US)
-    call hchksum(Kd_Int, "after set_diffusivity Kd_Int", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+    call hchksum(Kd_Int, "after set_diffusivity Kd_Int", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
   endif
 
   ! Set diffusivities for heat and salt separately
@@ -669,8 +669,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
     endif
 
     if (CS%debug) then
-      call hchksum(Kd_heat, "after set_diffusivity Kd_heat", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
-      call hchksum(Kd_salt, "after set_diffusivity Kd_salt", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+      call hchksum(Kd_heat, "after set_diffusivity Kd_heat", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
+      call hchksum(Kd_salt, "after set_diffusivity Kd_salt", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
     endif
 
     call cpu_clock_begin(id_clock_kpp)
@@ -727,12 +727,12 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
       call MOM_state_chksum("after KPP", u, v, h, G, GV, US, haloshift=0)
       call MOM_forcing_chksum("after KPP", fluxes, G, US, haloshift=0)
       call MOM_thermovar_chksum("after KPP", tv, G, US)
-      call hchksum(Kd_heat, "after KPP Kd_heat", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
-      call hchksum(Kd_salt, "after KPP Kd_salt", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+      call hchksum(Kd_heat, "after KPP Kd_heat", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
+      call hchksum(Kd_salt, "after KPP Kd_salt", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
       call hchksum(KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, &
-                   scale=US%C_to_degC*GV%H_to_m*US%s_to_T)
+                   unscale=US%C_to_degC*GV%H_to_m*US%s_to_T)
       call hchksum(KPP_salt_flux, "before KPP_applyNLT netSalt", G%HI, haloshift=0, &
-                   scale=US%S_to_ppt*GV%H_to_m*US%s_to_T)
+                   unscale=US%S_to_ppt*GV%H_to_m*US%s_to_T)
       call hchksum(KPP_NLTheat, "before KPP_applyNLT NLTheat", G%HI, haloshift=0)
       call hchksum(KPP_NLTscalar, "before KPP_applyNLT NLTscalar", G%HI, haloshift=0)
     endif
@@ -802,7 +802,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
     call MOM_forcing_chksum("after calc_entrain ", fluxes, G, US, haloshift=0)
     call MOM_thermovar_chksum("after calc_entrain ", tv, G, US)
     call MOM_state_chksum("after calc_entrain ", u, v, h, G, GV, US, haloshift=0)
-    call hchksum(ent_s, "after calc_entrain ent_s", G%HI, haloshift=0, scale=GV%H_to_MKS)
+    call hchksum(ent_s, "after calc_entrain ent_s", G%HI, haloshift=0, unscale=GV%H_to_MKS)
   endif
 
   ! Save fields before boundary forcing is applied for tendency diagnostics
@@ -828,19 +828,19 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
             CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD_h=visc%h_ML)
 
     if (CS%debug) then
-      call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, scale=GV%H_to_mks)
-      call hchksum(ent_s, "after applyBoundaryFluxes ent_s", G%HI, haloshift=0, scale=GV%H_to_mks)
+      call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, unscale=GV%H_to_mks)
+      call hchksum(ent_s, "after applyBoundaryFluxes ent_s", G%HI, haloshift=0, unscale=GV%H_to_mks)
       call hchksum(cTKE, "after applyBoundaryFluxes cTKE", G%HI, haloshift=0, &
-                   scale=US%RZ3_T3_to_W_m2*US%T_to_s)
+                   unscale=US%RZ3_T3_to_W_m2*US%T_to_s)
       call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT", G%HI, haloshift=0, &
-                   scale=US%kg_m3_to_R*US%degC_to_C)
+                   unscale=US%kg_m3_to_R*US%degC_to_C)
       call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS", G%HI, haloshift=0, &
-                   scale=US%kg_m3_to_R*US%ppt_to_S)
-      call hchksum(h, "after applyBoundaryFluxes h", G%HI, haloshift=0, scale=GV%H_to_mks)
-      call hchksum(tv%T, "after applyBoundaryFluxes tv%T", G%HI, haloshift=0, scale=US%C_to_degC)
-      call hchksum(tv%S, "after applyBoundaryFluxes tv%S", G%HI, haloshift=0, scale=US%S_to_ppt)
+                   unscale=US%kg_m3_to_R*US%ppt_to_S)
+      call hchksum(h, "after applyBoundaryFluxes h", G%HI, haloshift=0, unscale=GV%H_to_mks)
+      call hchksum(tv%T, "after applyBoundaryFluxes tv%T", G%HI, haloshift=0, unscale=US%C_to_degC)
+      call hchksum(tv%S, "after applyBoundaryFluxes tv%S", G%HI, haloshift=0, unscale=US%S_to_ppt)
       call hchksum(SkinBuoyFlux, "after applyBdryFlux SkinBuoyFlux", G%HI, haloshift=0, &
-                   scale=US%Z_to_m**2*US%s_to_T**3)
+                   unscale=US%Z_to_m**2*US%s_to_T**3)
     endif
 
     call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
@@ -877,9 +877,9 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
     enddo ; enddo ; enddo
 
     if (CS%debug) then
-      call hchksum(ent_t, "after ePBL ent_t", G%HI, haloshift=0, scale=GV%H_to_MKS)
-      call hchksum(ent_s, "after ePBL ent_s", G%HI, haloshift=0, scale=GV%H_to_MKS)
-      call hchksum(Kd_ePBL, "after ePBL Kd_ePBL", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+      call hchksum(ent_t, "after ePBL ent_t", G%HI, haloshift=0, unscale=GV%H_to_MKS)
+      call hchksum(ent_s, "after ePBL ent_s", G%HI, haloshift=0, unscale=GV%H_to_MKS)
+      call hchksum(Kd_ePBL, "after ePBL Kd_ePBL", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
     endif
 
   else
@@ -922,8 +922,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
   if (associated(tv%T)) then
 
     if (CS%debug) then
-      call hchksum(ent_t, "before triDiagTS ent_t ", G%HI, haloshift=0, scale=GV%H_to_MKS)
-      call hchksum(ent_s, "before triDiagTS ent_s ", G%HI, haloshift=0, scale=GV%H_to_MKS)
+      call hchksum(ent_t, "before triDiagTS ent_t ", G%HI, haloshift=0, unscale=GV%H_to_MKS)
+      call hchksum(ent_s, "before triDiagTS ent_s ", G%HI, haloshift=0, unscale=GV%H_to_MKS)
     endif
 
     call cpu_clock_begin(id_clock_tridiag)
@@ -1255,7 +1255,7 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
     call MOM_state_chksum("after set_diffusivity ", u, v, h, G, GV, US, haloshift=0)
     call MOM_forcing_chksum("after set_diffusivity ", fluxes, G, US, haloshift=0)
     call MOM_thermovar_chksum("after set_diffusivity ", tv, G, US)
-    call hchksum(Kd_heat, "after set_diffusivity Kd_heat", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+    call hchksum(Kd_heat, "after set_diffusivity Kd_heat", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
   endif
 
   ! Set diffusivities for heat and salt separately, and possibly change the meaning of Kd_heat.
@@ -1274,8 +1274,8 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
   endif
 
   if (CS%debug) then
-    call hchksum(Kd_heat, "after double diffuse Kd_heat", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
-    call hchksum(Kd_salt, "after double diffuse Kd_salt", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+    call hchksum(Kd_heat, "after double diffuse Kd_heat", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
+    call hchksum(Kd_salt, "after double diffuse Kd_salt", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
   endif
 
   if (CS%useKPP) then
@@ -1322,12 +1322,12 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
       call MOM_state_chksum("after KPP", u, v, h, G, GV, US, haloshift=0)
       call MOM_forcing_chksum("after KPP", fluxes, G, US, haloshift=0)
       call MOM_thermovar_chksum("after KPP", tv, G, US)
-      call hchksum(Kd_heat, "after KPP Kd_heat", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
-      call hchksum(Kd_salt, "after KPP Kd_salt", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+      call hchksum(Kd_heat, "after KPP Kd_heat", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
+      call hchksum(Kd_salt, "after KPP Kd_salt", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
       call hchksum(KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, &
-                   scale=US%C_to_degC*GV%H_to_m*US%s_to_T)
+                   unscale=US%C_to_degC*GV%H_to_m*US%s_to_T)
       call hchksum(KPP_salt_flux, "before KPP_applyNLT netSalt", G%HI, haloshift=0, &
-                   scale=US%S_to_ppt*GV%H_to_m*US%s_to_T)
+                   unscale=US%S_to_ppt*GV%H_to_m*US%s_to_T)
       call hchksum(KPP_NLTheat, "before KPP_applyNLT NLTheat", G%HI, haloshift=0)
       call hchksum(KPP_NLTscalar, "before KPP_applyNLT NLTscalar", G%HI, haloshift=0)
     endif
@@ -1377,14 +1377,14 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
             CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD_h=visc%h_ML)
 
     if (CS%debug) then
-      call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, scale=GV%H_to_MKS)
-      call hchksum(ent_s, "after applyBoundaryFluxes ent_s", G%HI, haloshift=0, scale=GV%H_to_MKS)
+      call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, unscale=GV%H_to_MKS)
+      call hchksum(ent_s, "after applyBoundaryFluxes ent_s", G%HI, haloshift=0, unscale=GV%H_to_MKS)
       call hchksum(cTKE, "after applyBoundaryFluxes cTKE", G%HI, haloshift=0, &
-                   scale=US%RZ3_T3_to_W_m2*US%T_to_s)
+                   unscale=US%RZ3_T3_to_W_m2*US%T_to_s)
       call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT", G%HI, haloshift=0, &
-                   scale=US%kg_m3_to_R*US%degC_to_C)
+                   unscale=US%kg_m3_to_R*US%degC_to_C)
       call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS", G%HI, haloshift=0, &
-                   scale=US%kg_m3_to_R*US%ppt_to_S)
+                   unscale=US%kg_m3_to_R*US%ppt_to_S)
     endif
 
     call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
@@ -1412,9 +1412,9 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
     enddo ; enddo ; enddo
 
     if (CS%debug) then
-      call hchksum(ent_t, "after ePBL ent_t", G%HI, haloshift=0, scale=GV%H_to_MKS)
-      call hchksum(ent_s, "after ePBL ent_s", G%HI, haloshift=0, scale=GV%H_to_MKS)
-      call hchksum(Kd_ePBL, "after ePBL Kd_ePBL", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+      call hchksum(ent_t, "after ePBL ent_t", G%HI, haloshift=0, unscale=GV%H_to_MKS)
+      call hchksum(ent_s, "after ePBL ent_s", G%HI, haloshift=0, unscale=GV%H_to_MKS)
+      call hchksum(Kd_ePBL, "after ePBL Kd_ePBL", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
     endif
 
   else
@@ -1449,8 +1449,8 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
   if (associated(tv%T)) then
 
     if (CS%debug) then
-      call hchksum(ent_t, "before triDiagTS ent_t ", G%HI, haloshift=0, scale=GV%H_to_MKS)
-      call hchksum(ent_s, "before triDiagTS ent_s ", G%HI, haloshift=0, scale=GV%H_to_MKS)
+      call hchksum(ent_t, "before triDiagTS ent_t ", G%HI, haloshift=0, unscale=GV%H_to_MKS)
+      call hchksum(ent_s, "before triDiagTS ent_s ", G%HI, haloshift=0, unscale=GV%H_to_MKS)
     endif
 
     call cpu_clock_begin(id_clock_tridiag)
@@ -1844,8 +1844,8 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
     if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then
       call find_uv_at_h(u, v, h_orig, u_h, v_h, G, GV, US, eaml, ebml)
       if (CS%debug) then
-        call hchksum(eaml, "after find_uv_at_h eaml", G%HI, scale=GV%H_to_MKS)
-        call hchksum(ebml, "after find_uv_at_h ebml", G%HI, scale=GV%H_to_MKS)
+        call hchksum(eaml, "after find_uv_at_h eaml", G%HI, unscale=GV%H_to_MKS)
+        call hchksum(ebml, "after find_uv_at_h ebml", G%HI, unscale=GV%H_to_MKS)
       endif
     else
       call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
@@ -1883,8 +1883,8 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
     call MOM_state_chksum("after set_diffusivity ", u, v, h, G, GV, US, haloshift=0)
     call MOM_forcing_chksum("after set_diffusivity ", fluxes, G, US, haloshift=0)
     call MOM_thermovar_chksum("after set_diffusivity ", tv, G, US)
-    call hchksum(Kd_lay, "after set_diffusivity Kd_lay", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
-    call hchksum(Kd_Int, "after set_diffusivity Kd_Int", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+    call hchksum(Kd_lay, "after set_diffusivity Kd_lay", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
+    call hchksum(Kd_Int, "after set_diffusivity Kd_Int", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
   endif
 
 
@@ -1959,8 +1959,8 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
       call MOM_state_chksum("after KPP", u, v, h, G, GV, US, haloshift=0)
       call MOM_forcing_chksum("after KPP", fluxes, G, US, haloshift=0)
       call MOM_thermovar_chksum("after KPP", tv, G, US)
-      call hchksum(Kd_lay, "after KPP Kd_lay", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
-      call hchksum(Kd_Int, "after KPP Kd_Int", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+      call hchksum(Kd_lay, "after KPP Kd_lay", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
+      call hchksum(Kd_Int, "after KPP Kd_Int", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
     endif
   endif  ! endif for KPP
 
@@ -1973,9 +1973,9 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
     call cpu_clock_begin(id_clock_kpp)
     if (CS%debug) then
       call hchksum(KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, &
-                   scale=US%C_to_degC*GV%H_to_m*US%s_to_T)
+                   unscale=US%C_to_degC*GV%H_to_m*US%s_to_T)
       call hchksum(KPP_salt_flux, "before KPP_applyNLT netSalt", G%HI, haloshift=0, &
-                   scale=US%S_to_ppt*GV%H_to_m*US%s_to_T)
+                   unscale=US%S_to_ppt*GV%H_to_m*US%s_to_T)
       call hchksum(KPP_NLTheat, "before KPP_applyNLT NLTheat", G%HI, haloshift=0)
       call hchksum(KPP_NLTscalar, "before KPP_applyNLT NLTscalar", G%HI, haloshift=0)
     endif
@@ -2032,8 +2032,8 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
     call MOM_forcing_chksum("after calc_entrain ", fluxes, G, US, haloshift=0)
     call MOM_thermovar_chksum("after calc_entrain ", tv, G, US)
     call MOM_state_chksum("after calc_entrain ", u, v, h, G, GV, US, haloshift=0)
-    call hchksum(ea, "after calc_entrain ea", G%HI, haloshift=0, scale=GV%H_to_MKS)
-    call hchksum(eb, "after calc_entrain eb", G%HI, haloshift=0, scale=GV%H_to_MKS)
+    call hchksum(ea, "after calc_entrain ea", G%HI, haloshift=0, unscale=GV%H_to_MKS)
+    call hchksum(eb, "after calc_entrain eb", G%HI, haloshift=0, unscale=GV%H_to_MKS)
   endif
 
   ! Save fields before boundary forcing is applied for tendency diagnostics
@@ -2188,8 +2188,8 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
         eb(i,j,k) = eb(i,j,k) + ebml(i,j,k)
       enddo ; enddo ; enddo
       if (CS%debug) then
-        call hchksum(ea, "after ea = ea + eaml", G%HI, haloshift=0, scale=GV%H_to_MKS)
-        call hchksum(eb, "after eb = eb + ebml", G%HI, haloshift=0, scale=GV%H_to_MKS)
+        call hchksum(ea, "after ea = ea + eaml", G%HI, haloshift=0, unscale=GV%H_to_MKS)
+        call hchksum(eb, "after eb = eb + ebml", G%HI, haloshift=0, unscale=GV%H_to_MKS)
       endif
     endif
 
@@ -2234,8 +2234,8 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
     if (associated(tv%T)) then
 
       if (CS%debug) then
-        call hchksum(ea, "before triDiagTS ea ", G%HI, haloshift=0, scale=GV%H_to_MKS)
-        call hchksum(eb, "before triDiagTS eb ", G%HI, haloshift=0, scale=GV%H_to_MKS)
+        call hchksum(ea, "before triDiagTS ea ", G%HI, haloshift=0, unscale=GV%H_to_MKS)
+        call hchksum(eb, "before triDiagTS eb ", G%HI, haloshift=0, unscale=GV%H_to_MKS)
       endif
       call cpu_clock_begin(id_clock_tridiag)
 
@@ -2282,8 +2282,8 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
   if (CS%debug) then
     call MOM_state_chksum("after mixed layer ", u, v, h, G, GV, US, haloshift=0)
     call MOM_thermovar_chksum("after mixed layer ", tv, G, US)
-    call hchksum(ea, "after mixed layer ea", G%HI, scale=GV%H_to_MKS)
-    call hchksum(eb, "after mixed layer eb", G%HI, scale=GV%H_to_MKS)
+    call hchksum(ea, "after mixed layer ea", G%HI, unscale=GV%H_to_MKS)
+    call hchksum(eb, "after mixed layer eb", G%HI, unscale=GV%H_to_MKS)
   endif
 
   call cpu_clock_begin(id_clock_remap)
@@ -2473,8 +2473,8 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
 ! mixed layer turbulence is applied elsewhere.
   if (CS%use_bulkmixedlayer) then
     if (CS%debug) then
-      call hchksum(ea, "before net flux rearrangement ea", G%HI, scale=GV%H_to_MKS)
-      call hchksum(eb, "before net flux rearrangement eb", G%HI, scale=GV%H_to_MKS)
+      call hchksum(ea, "before net flux rearrangement ea", G%HI, unscale=GV%H_to_MKS)
+      call hchksum(eb, "before net flux rearrangement eb", G%HI, unscale=GV%H_to_MKS)
     endif
     !$OMP parallel do default(shared) private(net_ent)
     do j=js,je
@@ -2485,8 +2485,8 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
       enddo ; enddo
     enddo
     if (CS%debug) then
-      call hchksum(ea, "after net flux rearrangement ea", G%HI, scale=GV%H_to_MKS)
-      call hchksum(eb, "after net flux rearrangement eb", G%HI, scale=GV%H_to_MKS)
+      call hchksum(ea, "after net flux rearrangement ea", G%HI, unscale=GV%H_to_MKS)
+      call hchksum(eb, "after net flux rearrangement eb", G%HI, unscale=GV%H_to_MKS)
     endif
   endif
 
@@ -2517,9 +2517,9 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
   !  or enters the ocean with the surface velocity.
   if (CS%debug) then
     call MOM_state_chksum("before u/v tridiag ", u, v, h, G, GV, US, haloshift=0)
-    call hchksum(ea, "before u/v tridiag ea", G%HI, scale=GV%H_to_MKS)
-    call hchksum(eb, "before u/v tridiag eb", G%HI, scale=GV%H_to_MKS)
-    call hchksum(hold, "before u/v tridiag hold", G%HI, scale=GV%H_to_MKS)
+    call hchksum(ea, "before u/v tridiag ea", G%HI, unscale=GV%H_to_MKS)
+    call hchksum(eb, "before u/v tridiag eb", G%HI, unscale=GV%H_to_MKS)
+    call hchksum(hold, "before u/v tridiag hold", G%HI, unscale=GV%H_to_MKS)
   endif
   call cpu_clock_begin(id_clock_tridiag)
 
diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90
index 7280106125..aa02a42ea9 100644
--- a/src/parameterizations/vertical/MOM_internal_tide_input.F90
+++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90
@@ -179,9 +179,9 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS)
   endif
 
   if (CS%debug) then
-    call hchksum(N2_bot,"N2_bot",G%HI,haloshift=0, scale=US%s_to_T**2)
-    call hchksum(CS%TKE_itidal_input,"TKE_itidal_input",G%HI,haloshift=0, &
-                 scale=US%RZ3_T3_to_W_m2)
+    call hchksum(N2_bot, "N2_bot", G%HI, haloshift=0, unscale=US%s_to_T**2)
+    call hchksum(CS%TKE_itidal_input,"TKE_itidal_input", G%HI, haloshift=0, &
+                 unscale=US%RZ3_T3_to_W_m2)
   endif
 
   call enable_averages(dt, time_end, CS%diag)
diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90
index dac8508b4d..f553108ac0 100644
--- a/src/parameterizations/vertical/MOM_kappa_shear.F90
+++ b/src/parameterizations/vertical/MOM_kappa_shear.F90
@@ -337,8 +337,8 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, &
   enddo ! end of j-loop
 
   if (CS%debug) then
-    call hchksum(kappa_io, "kappa", G%HI, scale=GV%HZ_T_to_m2_s)
-    call hchksum(tke_io, "tke", G%HI, scale=US%Z_to_m**2*US%s_to_T**2)
+    call hchksum(kappa_io, "kappa", G%HI, unscale=GV%HZ_T_to_m2_s)
+    call hchksum(tke_io, "tke", G%HI, unscale=US%Z_to_m**2*US%s_to_T**2)
   endif
 
   if (CS%id_Kd_shear > 0) call post_data(CS%id_Kd_shear, kappa_io, CS%diag)
@@ -621,8 +621,8 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_
   enddo ! end of J-loop
 
   if (CS%debug) then
-    call hchksum(kappa_io, "kappa", G%HI, scale=GV%HZ_T_to_m2_s)
-    call Bchksum(tke_io, "tke", G%HI, scale=US%Z_to_m**2*US%s_to_T**2)
+    call hchksum(kappa_io, "kappa", G%HI, unscale=GV%HZ_T_to_m2_s)
+    call Bchksum(tke_io, "tke", G%HI, unscale=US%Z_to_m**2*US%s_to_T**2)
   endif
 
   if (CS%id_Kd_shear > 0) call post_data(CS%id_Kd_shear, kappa_io, CS%diag)
diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90
index d78b675c1a..754d2c4b8e 100644
--- a/src/parameterizations/vertical/MOM_set_diffusivity.F90
+++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90
@@ -348,7 +348,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
 
   if (CS%useKappaShear) then
     if (CS%debug) then
-      call hchksum_pair("before calc_KS [uv]_h", u_h, v_h, G%HI, scale=US%L_T_to_m_s)
+      call hchksum_pair("before calc_KS [uv]_h", u_h, v_h, G%HI, unscale=US%L_T_to_m_s)
     endif
     call cpu_clock_begin(id_clock_kappaShear)
     if (CS%Vertex_shear) then
@@ -359,18 +359,18 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
                                    visc%TKE_turb, visc%Kv_shear_Bu, dt, G, GV, US, CS%kappaShear_CSp)
       if (associated(visc%Kv_shear)) visc%Kv_shear(:,:,:) = 0.0 ! needed for other parameterizations
       if (CS%debug) then
-        call hchksum(visc%Kd_shear, "after calc_KS_vert visc%Kd_shear", G%HI, scale=GV%HZ_T_to_m2_s)
-        call Bchksum(visc%Kv_shear_Bu, "after calc_KS_vert visc%Kv_shear_Bu", G%HI, scale=GV%HZ_T_to_m2_s)
-        call Bchksum(visc%TKE_turb, "after calc_KS_vert visc%TKE_turb", G%HI, scale=US%Z_to_m**2*US%s_to_T**2)
+        call hchksum(visc%Kd_shear, "after calc_KS_vert visc%Kd_shear", G%HI, unscale=GV%HZ_T_to_m2_s)
+        call Bchksum(visc%Kv_shear_Bu, "after calc_KS_vert visc%Kv_shear_Bu", G%HI, unscale=GV%HZ_T_to_m2_s)
+        call Bchksum(visc%TKE_turb, "after calc_KS_vert visc%TKE_turb", G%HI, unscale=US%Z_to_m**2*US%s_to_T**2)
       endif
     else
       ! Changes: visc%Kd_shear ;  Sets: visc%Kv_shear and visc%TKE_turb
       call calculate_kappa_shear(u_h, v_h, h, tv, fluxes%p_surf, visc%Kd_shear, visc%TKE_turb, &
                                  visc%Kv_shear, dt, G, GV, US, CS%kappaShear_CSp)
       if (CS%debug) then
-        call hchksum(visc%Kd_shear, "after calc_KS visc%Kd_shear", G%HI, scale=GV%HZ_T_to_m2_s)
-        call hchksum(visc%Kv_shear, "after calc_KS visc%Kv_shear", G%HI, scale=GV%HZ_T_to_m2_s)
-        call hchksum(visc%TKE_turb, "after calc_KS visc%TKE_turb", G%HI, scale=US%Z_to_m**2*US%s_to_T**2)
+        call hchksum(visc%Kd_shear, "after calc_KS visc%Kd_shear", G%HI, unscale=GV%HZ_T_to_m2_s)
+        call hchksum(visc%Kv_shear, "after calc_KS visc%Kv_shear", G%HI, unscale=GV%HZ_T_to_m2_s)
+        call hchksum(visc%TKE_turb, "after calc_KS visc%TKE_turb", G%HI, unscale=US%Z_to_m**2*US%s_to_T**2)
       endif
     endif
     call cpu_clock_end(id_clock_kappaShear)
@@ -379,8 +379,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
     !NOTE{BGR}: this needs to be cleaned up.  It works in 1D case, but has not been tested outside.
     call calculate_CVMix_shear(u_h, v_h, h, tv, visc%Kd_shear, visc%Kv_shear, G, GV, US, CS%CVMix_shear_CSp)
     if (CS%debug) then
-      call hchksum(visc%Kd_shear, "after CVMix_shear visc%Kd_shear", G%HI, scale=GV%HZ_T_to_m2_s)
-      call hchksum(visc%Kv_shear, "after CVMix_shear visc%Kv_shear", G%HI, scale=GV%HZ_T_to_m2_s)
+      call hchksum(visc%Kd_shear, "after CVMix_shear visc%Kd_shear", G%HI, unscale=GV%HZ_T_to_m2_s)
+      call hchksum(visc%Kv_shear, "after CVMix_shear visc%Kv_shear", G%HI, unscale=GV%HZ_T_to_m2_s)
     endif
   elseif (associated(visc%Kv_shear)) then
     visc%Kv_shear(:,:,:) = 0.0 ! needed if calculate_kappa_shear is not enabled
@@ -389,15 +389,15 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
   ! Smooth the properties through massless layers.
   if (use_EOS) then
     if (CS%debug) then
-      call hchksum(tv%T, "before vert_fill_TS tv%T", G%HI, scale=US%C_to_degC)
-      call hchksum(tv%S, "before vert_fill_TS tv%S", G%HI, scale=US%S_to_ppt)
-      call hchksum(h, "before vert_fill_TS h",G%HI, scale=GV%H_to_m)
+      call hchksum(tv%T, "before vert_fill_TS tv%T", G%HI, unscale=US%C_to_degC)
+      call hchksum(tv%S, "before vert_fill_TS tv%S", G%HI, unscale=US%S_to_ppt)
+      call hchksum(h, "before vert_fill_TS h",G%HI, unscale=GV%H_to_m)
     endif
     call vert_fill_TS(h, tv%T, tv%S, kappa_dt_fill, T_f, S_f, G, GV, US, larger_h_denom=.true.)
     if (CS%debug) then
-      call hchksum(tv%T, "after vert_fill_TS tv%T", G%HI, scale=US%C_to_degC)
-      call hchksum(tv%S, "after vert_fill_TS tv%S", G%HI, scale=US%S_to_ppt)
-      call hchksum(h, "after vert_fill_TS h",G%HI, scale=GV%H_to_m)
+      call hchksum(tv%T, "after vert_fill_TS tv%T", G%HI, unscale=US%C_to_degC)
+      call hchksum(tv%S, "after vert_fill_TS tv%S", G%HI, unscale=US%S_to_ppt)
+      call hchksum(h, "after vert_fill_TS h",G%HI, unscale=GV%H_to_m)
     endif
   endif
 
@@ -596,30 +596,30 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
   endif
 
   if (CS%debug) then
-    if (present(Kd_lay)) call hchksum(Kd_lay, "Kd_lay", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+    if (present(Kd_lay)) call hchksum(Kd_lay, "Kd_lay", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
 
-    if (CS%useKappaShear) call hchksum(visc%Kd_shear, "Turbulent Kd", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+    if (CS%useKappaShear) call hchksum(visc%Kd_shear, "Turbulent Kd", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
 
     if (CS%use_CVMix_ddiff) then
-      call hchksum(Kd_extra_T, "MOM_set_diffusivity: Kd_extra_T", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
-      call hchksum(Kd_extra_S, "MOM_set_diffusivity: Kd_extra_S", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s)
+      call hchksum(Kd_extra_T, "MOM_set_diffusivity: Kd_extra_T", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
+      call hchksum(Kd_extra_S, "MOM_set_diffusivity: Kd_extra_S", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
     endif
 
     if (allocated(visc%kv_bbl_u) .and. allocated(visc%kv_bbl_v)) then
       call uvchksum("BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, G%HI, &
-                    haloshift=0, symmetric=.true., scale=GV%HZ_T_to_m2_s, &
+                    haloshift=0, symmetric=.true., unscale=GV%HZ_T_to_m2_s, &
                     scalar_pair=.true.)
     endif
 
     if (allocated(visc%bbl_thick_u) .and. allocated(visc%bbl_thick_v)) then
       call uvchksum("BBL bbl_thick_[uv]", visc%bbl_thick_u, visc%bbl_thick_v, &
-                    G%HI, haloshift=0, symmetric=.true., scale=US%Z_to_m, &
+                    G%HI, haloshift=0, symmetric=.true., unscale=US%Z_to_m, &
                     scalar_pair=.true.)
     endif
 
     if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) then
       call uvchksum("Ray_[uv]", visc%Ray_u, visc%Ray_v, G%HI, 0, &
-          symmetric=.true., scale=GV%H_to_m*US%s_to_T, scalar_pair=.true.)
+          symmetric=.true., unscale=GV%H_to_m*US%s_to_T, scalar_pair=.true.)
     endif
 
   endif
diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90
index e59e5d99aa..88845a390f 100644
--- a/src/parameterizations/vertical/MOM_set_viscosity.F90
+++ b/src/parameterizations/vertical/MOM_set_viscosity.F90
@@ -321,16 +321,18 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
   if (.not.CS%bottomdraglaw) return
 
   if (CS%debug) then
-    call uvchksum("Start set_viscous_BBL [uv]", u, v, G%HI, haloshift=1, scale=US%L_T_to_m_s)
-    call hchksum(h,"Start set_viscous_BBL h", G%HI, haloshift=1, scale=GV%H_to_m)
-    if (associated(tv%T)) call hchksum(tv%T, "Start set_viscous_BBL T", G%HI, haloshift=1, scale=US%C_to_degC)
-    if (associated(tv%S)) call hchksum(tv%S, "Start set_viscous_BBL S", G%HI, haloshift=1, scale=US%S_to_ppt)
+    call uvchksum("Start set_viscous_BBL [uv]", u, v, G%HI, haloshift=1, unscale=US%L_T_to_m_s)
+    call hchksum(h,"Start set_viscous_BBL h", G%HI, haloshift=1, unscale=GV%H_to_m)
+    if (associated(tv%T)) call hchksum(tv%T, "Start set_viscous_BBL T", G%HI, haloshift=1, unscale=US%C_to_degC)
+    if (associated(tv%S)) call hchksum(tv%S, "Start set_viscous_BBL S", G%HI, haloshift=1, unscale=US%S_to_ppt)
     if (allocated(tv%SpV_avg)) &
-      call hchksum(tv%SpV_avg, "Start set_viscous_BBL SpV_avg", G%HI, haloshift=1, scale=US%kg_m3_to_R)
+      call hchksum(tv%SpV_avg, "Start set_viscous_BBL SpV_avg", G%HI, haloshift=1, unscale=US%kg_m3_to_R)
     if (allocated(tv%SpV_avg)) call hchksum(tv%SpV_avg, "Cornerless SpV_avg", G%HI, &
-                                            haloshift=1, omit_corners=.true., scale=US%kg_m3_to_R)
-    if (associated(tv%T)) call hchksum(tv%T, "Cornerless T", G%HI, haloshift=1, omit_corners=.true., scale=US%C_to_degC)
-    if (associated(tv%S)) call hchksum(tv%S, "Cornerless S", G%HI, haloshift=1, omit_corners=.true., scale=US%S_to_ppt)
+                                            haloshift=1, omit_corners=.true., unscale=US%kg_m3_to_R)
+    if (associated(tv%T)) call hchksum(tv%T, "Cornerless T", G%HI, haloshift=1, &
+                                       omit_corners=.true., unscale=US%C_to_degC)
+    if (associated(tv%S)) call hchksum(tv%S, "Cornerless S", G%HI, haloshift=1, &
+                                       omit_corners=.true., unscale=US%S_to_ppt)
   endif
 
   use_BBL_EOS = associated(tv%eqn_of_state) .and. CS%BBL_use_EOS
@@ -1088,13 +1090,13 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
   if (CS%debug) then
     if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) &
         call uvchksum("Ray [uv]", visc%Ray_u, visc%Ray_v, G%HI, haloshift=0, &
-                      scale=GV%H_to_m*US%s_to_T, scalar_pair=.true.)
+                      unscale=GV%H_to_m*US%s_to_T, scalar_pair=.true.)
     if (allocated(visc%kv_bbl_u) .and. allocated(visc%kv_bbl_v)) &
         call uvchksum("kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, G%HI, &
-                      haloshift=0, scale=GV%HZ_T_to_m2_s, scalar_pair=.true.)
+                      haloshift=0, unscale=GV%HZ_T_to_m2_s, scalar_pair=.true.)
     if (allocated(visc%bbl_thick_u) .and. allocated(visc%bbl_thick_v)) &
         call uvchksum("bbl_thick_[uv]", visc%bbl_thick_u, visc%bbl_thick_v, &
-                      G%HI, haloshift=0, scale=US%Z_to_m, scalar_pair=.true.)
+                      G%HI, haloshift=0, unscale=US%Z_to_m, scalar_pair=.true.)
   endif
 
 end subroutine set_viscous_BBL
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index c26ee4ac75..95856f013e 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -1886,12 +1886,12 @@ subroutine vertvisc_coef(u, v, h, dz, forces, visc, tv, dt, G, GV, US, CS, OBC,
 
   if (CS%debug) then
     call uvchksum("vertvisc_coef h_[uv]", CS%h_u, CS%h_v, G%HI, haloshift=0, &
-                  scale=GV%H_to_m, scalar_pair=.true.)
+                  unscale=GV%H_to_m, scalar_pair=.true.)
     call uvchksum("vertvisc_coef a_[uv]", CS%a_u, CS%a_v, G%HI, haloshift=0, &
-                  scale=GV%H_to_m*US%s_to_T, scalar_pair=.true.)
+                  unscale=GV%H_to_m*US%s_to_T, scalar_pair=.true.)
     if (allocated(hML_u) .and. allocated(hML_v)) &
       call uvchksum("vertvisc_coef hML_[uv]", hML_u, hML_v, G%HI, &
-                    haloshift=0, scale=US%Z_to_m, scalar_pair=.true.)
+                    haloshift=0, unscale=US%Z_to_m, scalar_pair=.true.)
   endif
 
 ! Offer diagnostic fields for averaging.
diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90
index ef11739f33..fe20daaefd 100644
--- a/src/tracer/MOM_CFC_cap.F90
+++ b/src/tracer/MOM_CFC_cap.F90
@@ -565,7 +565,7 @@ subroutine CFC_cap_set_forcing(sfc_state, fluxes, day_start, day_interval, G, US
   if (CS%debug) then
     do m=1,NTR
       call hchksum(CS%CFC_data(m)%sfc_flux, trim(CS%CFC_data(m)%name)//" sfc_flux", G%HI, &
-                   scale=US%RZ_T_to_kg_m2s)
+                   unscale=US%RZ_T_to_kg_m2s)
     enddo
   endif
 
diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90
index 6168ec1d70..5591e74d97 100644
--- a/src/tracer/MOM_generic_tracer.F90
+++ b/src/tracer/MOM_generic_tracer.F90
@@ -570,7 +570,7 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml,
       surface_field(i,j) = tv%S(i,j,1)
       dz_ml(i,j) = US%Z_to_m * Hml(i,j)
     enddo ; enddo
-    sosga = global_area_mean(surface_field, G, scale=US%S_to_ppt)
+    sosga = global_area_mean(surface_field, G, unscale=US%S_to_ppt)
 
     !
     !Calculate tendencies (i.e., field changes at dt) from the sources / sinks
@@ -1003,7 +1003,7 @@ subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, GV, CS)
 
     dzt(:,:,:) = GV%H_to_m * h(:,:,:)
 
-    sosga = global_area_mean(sfc_state%SSS, G, scale=G%US%S_to_ppt)
+    sosga = global_area_mean(sfc_state%SSS, G, unscale=G%US%S_to_ppt)
 
     if ((G%US%C_to_degC == 1.0) .and. (G%US%S_to_ppt == 1.0)) then
       call generic_tracer_coupler_set(sfc_state%tr_fields, &
diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90
index b31ebba7c8..f772e0bc8a 100644
--- a/src/tracer/MOM_offline_main.F90
+++ b/src/tracer/MOM_offline_main.F90
@@ -306,8 +306,8 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, C
   enddo ; enddo ; enddo
 
   if (CS%debug) then
-    call hchksum(h_pre, "h_pre before transport", G%HI, scale=GV%H_to_MKS)
-    call uvchksum("[uv]htr_sub before transport", uhtr_sub, vhtr_sub, G%HI, scale=HL2_to_kg_scale)
+    call hchksum(h_pre, "h_pre before transport", G%HI, unscale=GV%H_to_MKS)
+    call uvchksum("[uv]htr_sub before transport", uhtr_sub, vhtr_sub, G%HI, unscale=HL2_to_kg_scale)
   endif
   tot_residual = remaining_transport_sum(G, GV, US, uhtr, vhtr, h_new)
   if (CS%print_adv_offline) then
@@ -325,8 +325,8 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, C
     enddo ; enddo ; enddo
 
     if (CS%debug) then
-      call hchksum(h_vol, "h_vol before advect", G%HI, scale=HL2_to_kg_scale)
-      call uvchksum("[uv]htr_sub before advect", uhtr_sub, vhtr_sub, G%HI, scale=HL2_to_kg_scale)
+      call hchksum(h_vol, "h_vol before advect", G%HI, unscale=HL2_to_kg_scale)
+      call uvchksum("[uv]htr_sub before advect", uhtr_sub, vhtr_sub, G%HI, unscale=HL2_to_kg_scale)
       write(debug_msg, '(A,I4.4)') 'Before advect ', iter
       call MOM_tracer_chkinv(debug_msg, G, GV, h_pre, CS%tracer_reg)
     endif
@@ -347,7 +347,7 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, C
       ! Do ALE remapping/regridding to allow for more advection to occur in the next iteration
       call pass_var(h_new,G%Domain)
       if (CS%debug) then
-        call hchksum(h_new,"h_new before ALE", G%HI, scale=GV%H_to_MKS)
+        call hchksum(h_new,"h_new before ALE", G%HI, unscale=GV%H_to_MKS)
         write(debug_msg, '(A,I4.4)') 'Before ALE ', iter
         call MOM_tracer_chkinv(debug_msg, G, GV, h_new, CS%tracer_reg)
       endif
@@ -373,7 +373,7 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, C
       call cpu_clock_end(id_clock_ALE)
 
       if (CS%debug) then
-        call hchksum(h_new, "h_new after ALE", G%HI, scale=GV%H_to_MKS)
+        call hchksum(h_new, "h_new after ALE", G%HI, unscale=GV%H_to_MKS)
         write(debug_msg, '(A,I4.4)') 'After ALE ', iter
         call MOM_tracer_chkinv(debug_msg, G, GV, h_new, CS%tracer_reg)
       endif
@@ -415,8 +415,8 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, C
   call pass_vector(uhtr, vhtr, G%Domain)
 
   if (CS%debug) then
-    call hchksum(h_pre, "h after offline_advection_ale", G%HI, scale=GV%H_to_MKS)
-    call uvchksum("[uv]htr after offline_advection_ale", uhtr, vhtr, G%HI, scale=HL2_to_kg_scale)
+    call hchksum(h_pre, "h after offline_advection_ale", G%HI, unscale=GV%H_to_MKS)
+    call uvchksum("[uv]htr after offline_advection_ale", uhtr, vhtr, G%HI, unscale=HL2_to_kg_scale)
     call MOM_tracer_chkinv("After offline_advection_ale", G, GV, h_pre, CS%tracer_reg)
   endif
 
@@ -501,7 +501,7 @@ subroutine offline_redistribute_residual(CS, G, GV, US, h_pre, uhtr, vhtr, conve
 
         if (CS%debug) then
           call MOM_tracer_chksum("Before upwards redistribute ", CS%tracer_Reg, G)
-          call uvchksum("[uv]tr before upwards redistribute", uhtr, vhtr, G%HI, scale=HL2_to_kg_scale)
+          call uvchksum("[uv]tr before upwards redistribute", uhtr, vhtr, G%HI, unscale=HL2_to_kg_scale)
         endif
 
         if (x_before_y) then
@@ -542,7 +542,7 @@ subroutine offline_redistribute_residual(CS, G, GV, US, h_pre, uhtr, vhtr, conve
 
         if (CS%debug) then
           call MOM_tracer_chksum("Before barotropic redistribute ", CS%tracer_Reg, G)
-          call uvchksum("[uv]tr before upwards redistribute", uhtr, vhtr, G%HI, scale=HL2_to_kg_scale)
+          call uvchksum("[uv]tr before upwards redistribute", uhtr, vhtr, G%HI, unscale=HL2_to_kg_scale)
         endif
 
         if (x_before_y) then
@@ -602,8 +602,8 @@ subroutine offline_redistribute_residual(CS, G, GV, US, h_pre, uhtr, vhtr, conve
   if (CS%id_vhr>0) call post_data(CS%id_vhr, vhtr, CS%diag)
 
   if (CS%debug) then
-    call hchksum(h_pre, "h_pre after redistribute", G%HI, scale=GV%H_to_MKS)
-    call uvchksum("uhtr after redistribute", uhtr, vhtr, G%HI, scale=HL2_to_kg_scale)
+    call hchksum(h_pre, "h_pre after redistribute", G%HI, unscale=GV%H_to_MKS)
+    call uvchksum("uhtr after redistribute", uhtr, vhtr, G%HI, unscale=HL2_to_kg_scale)
     call MOM_tracer_chkinv("after redistribute ", G, GV, h_new, CS%tracer_Reg)
   endif
 
@@ -685,9 +685,9 @@ subroutine offline_diabatic_ale(fluxes, Time_start, Time_end, G, GV, US, CS, h_p
   call MOM_mesg("Applying tracer source, sinks, and vertical mixing")
 
   if (CS%debug) then
-    call hchksum(h_pre, "h_pre before offline_diabatic_ale", G%HI, scale=GV%H_to_MKS)
-    call hchksum(eatr, "eatr before offline_diabatic_ale", G%HI, scale=GV%H_to_MKS)
-    call hchksum(ebtr, "ebtr before offline_diabatic_ale", G%HI, scale=GV%H_to_MKS)
+    call hchksum(h_pre, "h_pre before offline_diabatic_ale", G%HI, unscale=GV%H_to_MKS)
+    call hchksum(eatr, "eatr before offline_diabatic_ale", G%HI, unscale=GV%H_to_MKS)
+    call hchksum(ebtr, "ebtr before offline_diabatic_ale", G%HI, unscale=GV%H_to_MKS)
     call MOM_tracer_chkinv("Before offline_diabatic_ale", G, GV, h_pre, CS%tracer_reg)
   endif
 
@@ -751,9 +751,9 @@ subroutine offline_diabatic_ale(fluxes, Time_start, Time_end, G, GV, US, CS, h_p
   endif
 
   if (CS%debug) then
-    call hchksum(h_pre, "h_pre after offline_diabatic_ale", G%HI, scale=GV%H_to_MKS)
-    call hchksum(eatr, "eatr after offline_diabatic_ale", G%HI, scale=GV%H_to_MKS)
-    call hchksum(ebtr, "ebtr after offline_diabatic_ale", G%HI, scale=GV%H_to_MKS)
+    call hchksum(h_pre, "h_pre after offline_diabatic_ale", G%HI, unscale=GV%H_to_MKS)
+    call hchksum(eatr, "eatr after offline_diabatic_ale", G%HI, unscale=GV%H_to_MKS)
+    call hchksum(ebtr, "ebtr after offline_diabatic_ale", G%HI, unscale=GV%H_to_MKS)
     call MOM_tracer_chkinv("After offline_diabatic_ale", G, GV, h_pre, CS%tracer_reg)
   endif
 
@@ -794,7 +794,7 @@ subroutine offline_fw_fluxes_into_ocean(G, GV, CS, fluxes, h, in_flux_optional)
   enddo ; enddo
 
   if (CS%debug) then
-    call hchksum(h, "h before fluxes into ocean", G%HI, scale=GV%H_to_MKS)
+    call hchksum(h, "h before fluxes into ocean", G%HI, unscale=GV%H_to_MKS)
     call MOM_tracer_chkinv("Before fluxes into ocean", G, GV, h, CS%tracer_reg)
   endif
   do m = 1,CS%tracer_reg%ntr
@@ -804,7 +804,7 @@ subroutine offline_fw_fluxes_into_ocean(G, GV, CS, fluxes, h, in_flux_optional)
                                         CS%evap_CFL_limit, CS%minimum_forcing_depth, update_h_opt=update_h)
   enddo
   if (CS%debug) then
-    call hchksum(h, "h after fluxes into ocean", G%HI, scale=GV%H_to_MKS)
+    call hchksum(h, "h after fluxes into ocean", G%HI, unscale=GV%H_to_MKS)
     call MOM_tracer_chkinv("After fluxes into ocean", G, GV, h, CS%tracer_reg)
   endif
 
@@ -833,7 +833,7 @@ subroutine offline_fw_fluxes_out_ocean(G, GV, CS, fluxes, h, out_flux_optional)
     call MOM_error(WARNING, "Negative freshwater fluxes with non-zero tracer concentration not supported yet")
 
   if (CS%debug) then
-    call hchksum(h, "h before fluxes out of ocean", G%HI, scale=GV%H_to_MKS)
+    call hchksum(h, "h before fluxes out of ocean", G%HI, unscale=GV%H_to_MKS)
     call MOM_tracer_chkinv("Before fluxes out of ocean", G, GV, h, CS%tracer_reg)
   endif
   do m = 1, CS%tracer_reg%ntr
@@ -843,7 +843,7 @@ subroutine offline_fw_fluxes_out_ocean(G, GV, CS, fluxes, h, out_flux_optional)
                                         CS%evap_CFL_limit, CS%minimum_forcing_depth, update_h_opt = update_h)
   enddo
   if (CS%debug) then
-    call hchksum(h, "h after fluxes out of ocean", G%HI, scale=GV%H_to_MKS)
+    call hchksum(h, "h after fluxes out of ocean", G%HI, unscale=GV%H_to_MKS)
     call MOM_tracer_chkinv("Before fluxes out of ocean", G, GV, h, CS%tracer_reg)
   endif
 
@@ -1041,10 +1041,10 @@ subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale)
 
   if (CS%debug) then
     call uvchksum("[uv]htr before update_offline_fields", CS%uhtr, CS%vhtr, G%HI, &
-                  scale=US%L_to_m**2*GV%H_to_kg_m2)
-    call hchksum(CS%h_end, "h_end before update_offline_fields", G%HI, scale=GV%H_to_MKS)
-    call hchksum(CS%tv%T, "Temp before update_offline_fields", G%HI, scale=US%C_to_degC)
-    call hchksum(CS%tv%S, "Salt before update_offline_fields", G%HI, scale=US%S_to_ppt)
+                  unscale=US%L_to_m**2*GV%H_to_kg_m2)
+    call hchksum(CS%h_end, "h_end before update_offline_fields", G%HI, unscale=GV%H_to_MKS)
+    call hchksum(CS%tv%T, "Temp before update_offline_fields", G%HI, unscale=US%C_to_degC)
+    call hchksum(CS%tv%S, "Salt before update_offline_fields", G%HI, unscale=US%S_to_ppt)
   endif
 
   ! Store a copy of the layer thicknesses before ALE regrid/remap
@@ -1063,9 +1063,9 @@ subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale)
   endif
   if (CS%debug) then
     call uvchksum("[uv]h after update offline from files and arrays", CS%uhtr, CS%vhtr, G%HI, &
-                  scale=US%L_to_m**2*GV%H_to_kg_m2)
-    call hchksum(CS%tv%T, "Temp after update offline from files and arrays", G%HI, scale=US%C_to_degC)
-    call hchksum(CS%tv%S, "Salt after update offline from files and arrays", G%HI, scale=US%S_to_ppt)
+                  unscale=US%L_to_m**2*GV%H_to_kg_m2)
+    call hchksum(CS%tv%T, "Temp after update offline from files and arrays", G%HI, unscale=US%C_to_degC)
+    call hchksum(CS%tv%S, "Salt after update offline from files and arrays", G%HI, unscale=US%S_to_ppt)
   endif
 
   ! If using an ALE-dependent vertical coordinate, fields will need to be remapped
@@ -1083,8 +1083,8 @@ subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale)
     if (CS%id_h_regrid>0) call post_data(CS%id_h_regrid, h, CS%diag)
     if (CS%debug) then
       call uvchksum("[uv]htr after ALE regridding/remapping of inputs", CS%uhtr, CS%vhtr, G%HI, &
-                    scale=US%L_to_m**2*GV%H_to_kg_m2)
-      call hchksum(h_start,"h_start after ALE regridding/remapping of inputs", G%HI, scale=GV%H_to_MKS)
+                    unscale=US%L_to_m**2*GV%H_to_kg_m2)
+      call hchksum(h_start,"h_start after ALE regridding/remapping of inputs", G%HI, unscale=GV%H_to_MKS)
     endif
   endif
 
@@ -1131,10 +1131,10 @@ subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale)
 
   if (CS%debug) then
     call uvchksum("[uv]htr after update_offline_fields", CS%uhtr, CS%vhtr, G%HI, &
-                  scale=US%L_to_m**2*GV%H_to_kg_m2)
-    call hchksum(CS%h_end, "h_end after update_offline_fields", G%HI, scale=GV%H_to_MKS)
-    call hchksum(CS%tv%T, "Temp after update_offline_fields", G%HI, scale=US%C_to_degC)
-    call hchksum(CS%tv%S, "Salt after update_offline_fields", G%HI, scale=US%S_to_ppt)
+                  unscale=US%L_to_m**2*GV%H_to_kg_m2)
+    call hchksum(CS%h_end, "h_end after update_offline_fields", G%HI, unscale=GV%H_to_MKS)
+    call hchksum(CS%tv%T, "Temp after update_offline_fields", G%HI, unscale=US%C_to_degC)
+    call hchksum(CS%tv%S, "Salt after update_offline_fields", G%HI, unscale=US%S_to_ppt)
   endif
 
   call callTree_leave("update_offline_fields")
diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90
index 2b1530e94d..5538e210da 100644
--- a/src/tracer/MOM_tracer_hor_diff.F90
+++ b/src/tracer/MOM_tracer_hor_diff.F90
@@ -660,7 +660,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_
 
   if (CS%debug) then
     call uvchksum("After tracer diffusion khdt_[xy]", khdt_x, khdt_y, &
-                  G%HI, haloshift=0, symmetric=.true., scale=US%L_to_m**2, &
+                  G%HI, haloshift=0, symmetric=.true., unscale=US%L_to_m**2, &
                   scalar_pair=.true.)
   endif
 
diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90
index 0a5a8d4efd..835b93bb82 100644
--- a/src/tracer/MOM_tracer_registry.F90
+++ b/src/tracer/MOM_tracer_registry.F90
@@ -760,7 +760,7 @@ subroutine tracer_array_chksum(mesg, Tr, ntr, G)
   integer :: m
 
   do m=1,ntr
-    call hchksum(Tr(m)%t, mesg//trim(Tr(m)%name), G%HI, scale=Tr(m)%conc_scale)
+    call hchksum(Tr(m)%t, mesg//trim(Tr(m)%name), G%HI, unscale=Tr(m)%conc_scale)
   enddo
 
 end subroutine tracer_array_chksum
@@ -776,7 +776,7 @@ subroutine tracer_Reg_chksum(mesg, Reg, G)
   if (.not.associated(Reg)) return
 
   do m=1,Reg%ntr
-    call hchksum(Reg%Tr(m)%t, mesg//trim(Reg%Tr(m)%name), G%HI, scale=Reg%Tr(m)%conc_scale)
+    call hchksum(Reg%Tr(m)%t, mesg//trim(Reg%Tr(m)%name), G%HI, unscale=Reg%Tr(m)%conc_scale)
   enddo
 
 end subroutine tracer_Reg_chksum
diff --git a/src/tracer/pseudo_salt_tracer.F90 b/src/tracer/pseudo_salt_tracer.F90
index ade36bad19..b737dba5b7 100644
--- a/src/tracer/pseudo_salt_tracer.F90
+++ b/src/tracer/pseudo_salt_tracer.F90
@@ -214,7 +214,7 @@ subroutine pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G
   if (.not.associated(CS%ps)) return
 
   if (debug) then
-    call hchksum(tv%S,"salt pre pseudo-salt vertdiff", G%HI, scale=US%S_to_ppt)
+    call hchksum(tv%S,"salt pre pseudo-salt vertdiff", G%HI, unscale=US%S_to_ppt)
     call hchksum(CS%ps,"pseudo_salt pre pseudo-salt vertdiff", G%HI)
   endif
 
@@ -267,7 +267,7 @@ subroutine pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G
   endif
 
   if (debug) then
-    call hchksum(tv%S, "salt post pseudo-salt vertdiff", G%HI, scale=US%S_to_ppt)
+    call hchksum(tv%S, "salt post pseudo-salt vertdiff", G%HI, unscale=US%S_to_ppt)
     call hchksum(CS%ps, "pseudo_salt post pseudo-salt vertdiff", G%HI)
   endif