diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90
index 443188fb74..f3b72f9415 100644
--- a/config_src/solo_driver/MOM_driver.F90
+++ b/config_src/solo_driver/MOM_driver.F90
@@ -52,7 +52,7 @@ program MOM_main
   use MOM_error_handler,   only : callTree_enter, callTree_leave, callTree_waypoint
   use MOM_file_parser,     only : read_param, get_param, log_param, log_version, param_file_type
   use MOM_file_parser,     only : close_param_file
-  use MOM_forcing_type,    only : forcing, forcing_diagnostics, mech_forcing_diags
+  use MOM_forcing_type,    only : forcing, forcing_diagnostics, mech_forcing_diags, MOM_forcing_chksum
   use MOM_get_input,       only : directories
   use MOM_grid,            only : ocean_grid_type
   use MOM_io,              only : file_exists, open_file, close_file
@@ -385,6 +385,10 @@ program MOM_main
     ! Set the forcing for the next steps.
     call set_forcing(state, fluxes, Time, Time_step_ocean, grid, &
                      surface_forcing_CSp)
+    if (MOM_CSp%debug) then
+      call MOM_forcing_chksum("After set forcing", fluxes, grid, haloshift=0)
+    endif
+
     if (use_ice_shelf) then
       call shelf_calc_flux(state, fluxes, Time, time_step, ice_shelf_CSp)
 !###IS     call add_shelf_flux_forcing(fluxes, ice_shelf_CSp)
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index caf9e11041..c6017c7465 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -352,6 +352,7 @@ module MOM
 use MOM_error_handler,        only : callTree_enter, callTree_leave, callTree_waypoint
 use MOM_file_parser,          only : read_param, get_param, log_version, param_file_type
 use MOM_fixed_initialization, only : MOM_initialize_fixed
+use MOM_forcing_type,         only : MOM_forcing_chksum
 use MOM_get_input,            only : Get_MOM_Input, directories
 use MOM_io,                   only : MOM_io_init, vardesc
 use MOM_obsolete_params,      only : find_obsolete_params
@@ -850,6 +851,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
 
   if (CS%debug) then
     call MOM_state_chksum("Before steps ", u, v, h, CS%uh, CS%vh, G)
+    call MOM_forcing_chksum("Before steps", fluxes, G, haloshift=0)
     call check_redundant("Before steps ", u, v, G)
   endif
 
@@ -1172,9 +1174,9 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
       if (CS%debug) then
         call uchksum(u,"Pre-advection u", G, haloshift=2)
         call vchksum(v,"Pre-advection v", G, haloshift=2)
-        call hchksum(h,"Pre-advection h", G, haloshift=1)
-        call uchksum(CS%uhtr,"Pre-advection uhtr", G, haloshift=0)
-        call vchksum(CS%vhtr,"Pre-advection vhtr", G, haloshift=0)
+        call hchksum(h*G%H_to_m,"Pre-advection h", G, haloshift=1)
+        call uchksum(CS%uhtr*G%H_to_m,"Pre-advection uhtr", G, haloshift=0)
+        call vchksum(CS%vhtr*G%H_to_m,"Pre-advection vhtr", G, haloshift=0)
       ! call MOM_state_chksum("Pre-advection ", u, v, &
       !                       h, CS%uhtr, CS%vhtr, G, haloshift=1)
           if (associated(CS%tv%T)) call hchksum(CS%tv%T, "Pre-advection T", G, haloshift=1)
@@ -1224,12 +1226,13 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
         if (CS%debug) then
           call uchksum(u,"Pre-diabatic u", G, haloshift=2)
           call vchksum(v,"Pre-diabatic v", G, haloshift=2)
-          call hchksum(h,"Pre-diabatic h", G, haloshift=1)
-          call uchksum(CS%uhtr,"Pre-diabatic uh", G, haloshift=0)
-          call vchksum(CS%vhtr,"Pre-diabatic vh", G, haloshift=0)
+          call hchksum(h*G%H_to_m,"Pre-diabatic h", G, haloshift=1)
+          call uchksum(CS%uhtr*G%H_to_m,"Pre-diabatic uh", G, haloshift=0)
+          call vchksum(CS%vhtr*G%H_to_m,"Pre-diabatic vh", G, haloshift=0)
         ! call MOM_state_chksum("Pre-diabatic ",u, v, h, CS%uhtr, CS%vhtr, G)
           call MOM_thermo_chksum("Pre-diabatic ", CS%tv, G,haloshift=0)
           call check_redundant("Pre-diabatic ", u, v, G)
+          call MOM_forcing_chksum("Pre-diabatic", fluxes, G, haloshift=0)
         endif
 
         if (CS%split .and. CS%legacy_split) then  ! The dt here is correct. -RWH
@@ -1271,9 +1274,9 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
         if (CS%debug) then
           call uchksum(u,"Post-diabatic u", G, haloshift=2)
           call vchksum(v,"Post-diabatic v", G, haloshift=2)
-          call hchksum(h,"Post-diabatic h", G, haloshift=1)
-          call uchksum(CS%uhtr,"Post-diabatic uh", G, haloshift=0)
-          call vchksum(CS%vhtr,"Post-diabatic vh", G, haloshift=0)
+          call hchksum(h*G%H_to_m,"Post-diabatic h", G, haloshift=1)
+          call uchksum(CS%uhtr*G%H_to_m,"Post-diabatic uh", G, haloshift=0)
+          call vchksum(CS%vhtr*G%H_to_m,"Post-diabatic vh", G, haloshift=0)
         ! call MOM_state_chksum("Post-diabatic ", u, v, &
         !                       h, CS%uhtr, CS%vhtr, G, haloshift=1)
           if (associated(CS%tv%T)) call hchksum(CS%tv%T, "Post-diabatic T", G, haloshift=1)
diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index a7c1281469..08a7d62cb6 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -1460,8 +1460,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, &
   endif
 
   if (CS%debug) then
-    call uchksum(uhbt, "BT uhbt",CS%debug_BT_G,haloshift=0)
-    call vchksum(vhbt, "BT vhbt",CS%debug_BT_G,haloshift=0)
+    call uchksum(uhbt*G%H_to_m, "BT uhbt",CS%debug_BT_G,haloshift=0)
+    call vchksum(vhbt*G%H_to_m, "BT vhbt",CS%debug_BT_G,haloshift=0)
     call uchksum(ubt, "BT Initial ubt",CS%debug_BT_G,haloshift=0)
     call vchksum(vbt, "BT Initial vbt",CS%debug_BT_G,haloshift=0)
     call hchksum(G%H_to_kg_m2*eta, "BT Initial eta",CS%debug_BT_G,haloshift=0)
@@ -1476,8 +1476,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, &
     endif
     call uchksum(Cor_ref_u, "BT Cor_ref_u",CS%debug_BT_G,haloshift=0)
     call vchksum(Cor_ref_v, "BT Cor_ref_v",CS%debug_BT_G,haloshift=0)
-    call uchksum(uhbt0, "BT uhbt0",CS%debug_BT_G,haloshift=0)
-    call vchksum(vhbt0, "BT vhbt0",CS%debug_BT_G,haloshift=0)
+    call uchksum(uhbt0*G%H_to_m, "BT uhbt0",CS%debug_BT_G,haloshift=0)
+    call vchksum(vhbt0*G%H_to_m, "BT vhbt0",CS%debug_BT_G,haloshift=0)
     if (.not. use_BT_cont) then
       call uchksum(G%H_to_m*Datu, "BT Datu",CS%debug_BT_G,haloshift=1)
       call vchksum(G%H_to_m*Datv, "BT Datv",CS%debug_BT_G,haloshift=1)
@@ -1969,8 +1969,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, &
     endif
 
     if (CS%debug_bt) then
-      call uchksum(uhbt, "BT uhbt just after OBC",CS%debug_BT_G,haloshift=iev-ie)
-      call vchksum(vhbt, "BT vhbt just after OBC",CS%debug_BT_G,haloshift=iev-ie)
+      call uchksum(uhbt*G%H_to_m, "BT uhbt just after OBC",CS%debug_BT_G,haloshift=iev-ie)
+      call vchksum(vhbt*G%H_to_m, "BT vhbt just after OBC",CS%debug_BT_G,haloshift=iev-ie)
     endif
 
 !$OMP parallel do default(none) shared(isv,iev,jsv,jev,n,eta,eta_src,dtbt,CS,uhbt,vhbt,eta_wtd,wt_eta)
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index 3c602a5612..570424e17f 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -568,11 +568,16 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
       vp(i,j,k) = G%mask2dCv(i,j) * (v(i,j,k) + dt * v_bc_accel(i,J,k))
     enddo ; enddo 
   enddo
+
   call enable_averaging(dt, Time_local, CS%diag)
   call set_viscous_ML(u, v, h, tv, fluxes, visc, dt, G, &
                       CS%set_visc_CSp)
   call disable_averaging(CS%diag)
 
+  if (CS%debug) then
+    call uchksum(up,"before vertvisc: up",G,haloshift=0)
+    call vchksum(vp,"before vertvisc: vp",G,haloshift=0)
+  endif
   call vertvisc_coef(up, vp, h, fluxes, visc, dt, G, CS%vertvisc_CSp)
   call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, CS%vertvisc_CSp)
   call cpu_clock_end(id_clock_vertvisc)
@@ -671,6 +676,10 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
 ! up <- up + dt_pred d/dz visc d/dz up
 ! 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 uchksum(up,"0 before vertvisc: up",G,haloshift=0)
+    call vchksum(vp,"0 before vertvisc: vp",G,haloshift=0)
+  endif
   call vertvisc_coef(up, vp, h, fluxes, visc, dt_pred, G, CS%vertvisc_CSp)
   call vertvisc(up, vp, h, fluxes, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, G, &
                 CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot)
diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90
index 95675b7433..450da3ebdf 100644
--- a/src/initialization/MOM_state_initialization.F90
+++ b/src/initialization/MOM_state_initialization.F90
@@ -291,12 +291,6 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, PF, dirs, &
       endif
     endif  ! not from_Z_file.
 
-    if (debug) then
-      call hchksum(h, "MOM_initialize_state: h ", G, haloshift=1)
-      if ( use_temperature ) call hchksum(tv%T, "MOM_initialize_state: T ", G, haloshift=1)
-      if ( use_temperature ) call hchksum(tv%S, "MOM_initialize_state: S ", G, haloshift=1)
-    endif
-
 !   Initialize velocity components, u and v
     call get_param(PF, mod, "VELOCITY_CONFIG", config, &
          "A string that determines how the initial velocities \n"//&
@@ -324,7 +318,6 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, PF, dirs, &
     if (debug) call uchksum(u, "MOM_initialize_state: u ", G, haloshift=1)
     if (debug) call vchksum(v, "MOM_initialize_state: v ", G, haloshift=1)
 
-
 !   Optionally convert the thicknesses from m to kg m-2.  This is particularly
 ! useful in a non-Boussinesq model.
     call get_param(PF, mod, "CONVERT_THICKNESS_UNITS", convert, &
diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90
index e06fa2aff4..ab936f2a5e 100644
--- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90
+++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90
@@ -310,7 +310,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS)
     call vchksum(Kh_v(:,:,:),"Kh_v",G,haloshift=0)
     call uchksum(int_slope_u(:,:,:),"int_slope_u",G,haloshift=0)
     call vchksum(int_slope_v(:,:,:),"int_slope_v",G,haloshift=0)
-    call hchksum(h(:,:,:),"thickness_diffuse_1 h",G,haloshift=0)
+    call hchksum(h(:,:,:)*G%H_to_m,"thickness_diffuse_1 h",G,haloshift=0)
     call hchksum(e(:,:,:),"thickness_diffuse_1 e",G,haloshift=0)
   endif
 
@@ -349,11 +349,11 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS)
   endif
 
   if (CS%debug) then
-    call uchksum(uhD(:,:,:),"thickness_diffuse uhD",G,haloshift=0)
-    call vchksum(vhD(:,:,:),"thickness_diffuse vhD",G,haloshift=0)
-    call uchksum(uhtr(:,:,:),"thickness_diffuse uhtr",G,haloshift=0)
-    call vchksum(vhtr(:,:,:),"thickness_diffuse vhtr",G,haloshift=0)
-    call hchksum(h(:,:,:),"thickness_diffuse h",G,haloshift=0)
+    call uchksum(uhD(:,:,:)*G%H_to_m,"thickness_diffuse uhD",G,haloshift=0)
+    call vchksum(vhD(:,:,:)*G%H_to_m,"thickness_diffuse vhD",G,haloshift=0)
+    call uchksum(uhtr(:,:,:)*G%H_to_m,"thickness_diffuse uhtr",G,haloshift=0)
+    call vchksum(vhtr(:,:,:)*G%H_to_m,"thickness_diffuse vhtr",G,haloshift=0)
+    call hchksum(h(:,:,:)*G%H_to_m,"thickness_diffuse h",G,haloshift=0)
   endif
 
 ! Offer diagnostic fields for averaging.
diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90
index 25ee559b02..be40cf8151 100644
--- a/src/parameterizations/vertical/MOM_KPP.F90
+++ b/src/parameterizations/vertical/MOM_KPP.F90
@@ -433,7 +433,7 @@ subroutine KPP_calculate(CS, G, h, Temp, Salt, u, v, EOS, uStar, &
 
 #ifdef __DO_SAFETY_CHECKS__
   if (CS%debug) then
-    call hchksum(h, "KPP in: h",G,haloshift=0)
+    call hchksum(h*G%H_to_m, "KPP in: h",G,haloshift=0)
     call hchksum(Temp, "KPP in: T",G,haloshift=0)
     call hchksum(Salt, "KPP in: S",G,haloshift=0)
     call hchksum(u, "KPP in: u",G,haloshift=0)
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index 0eedb54441..9a81cf5131 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -374,6 +374,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS)
 
   if (CS%debug) then
     call MOM_state_chksum("Start of diabatic ", u(:,:,:), v(:,:,:), h(:,:,:), G)
+    call MOM_forcing_chksum("Start of diabatic", fluxes, G, haloshift=0)
   endif
   if (CS%debugConservation) call MOM_state_stats('Start of diabatic', u, v, h, tv%T, tv%S, G)
 
@@ -416,6 +417,10 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS)
     call set_opacity(CS%optics, fluxes, G, CS%opacity_CSp)
 
   if (CS%bulkmixedlayer) then
+    if (CS%debug) then
+      call MOM_forcing_chksum("Before mixedlayer", fluxes, G, haloshift=0)
+    endif
+
     if (CS%ML_mix_first > 0.0) then
 !  This subroutine (1)  Cools the mixed layer.
 !    (2) Performs convective adjustment by mixed layer entrainment.
@@ -445,7 +450,10 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS)
       if (ASSOCIATED(tv%S) .and. ASSOCIATED(tv%salt_deficit)) &
         call adjust_salt(h, tv, G, CS)
       call cpu_clock_end(id_clock_mixedlayer)
-      if (CS%debug) call MOM_state_chksum("After mixedlayer ", u, v, h, G)
+      if (CS%debug) then
+        call MOM_state_chksum("After mixedlayer ", u, v, h, G)
+        call MOM_forcing_chksum("After mixedlayer", fluxes, G, haloshift=0)
+      endif
       if (showCallTree) call callTree_waypoint("done with 1st bulkmixedlayer (diabatic)")
       if (CS%debugConservation) call MOM_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, G)
     endif
@@ -597,10 +605,10 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS)
   if (CS%useKPP) then
     call cpu_clock_begin(id_clock_kpp)
     if (CS%debug) then
+      call hchksum(CS%netHeatMinusSW*G%H_to_m, "before KPP_applyNLT netHeat",G,haloshift=0)
+      call hchksum(CS%netSalt*G%H_to_m, "before KPP_applyNLT netSalt",G,haloshift=0)
       call hchksum(CS%KPP_NLTheat, "before KPP_applyNLT NLTheat",G,haloshift=0)
-      call hchksum(CS%netHeatMinusSW, "before KPP_applyNLT netHeat",G,haloshift=0)
       call hchksum(CS%KPP_NLTscalar, "before KPP_applyNLT NLTscalar",G,haloshift=0)
-      call hchksum(CS%netSalt, "before KPP_applyNLT netSalt",G,haloshift=0)
     endif
     ! Apply non-local transport of heat and salt
     ! Changes: tv%T, tv%S
diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90
index c2c576b9b5..44ea4cf499 100644
--- a/src/parameterizations/vertical/MOM_set_diffusivity.F90
+++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90
@@ -519,13 +519,13 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, G, C
     if (CS%debug) then
       call hchksum(tv%T, "before vert_fill_TS tv%T",G)
       call hchksum(tv%S, "before vert_fill_TS tv%S",G)
-      call hchksum(h, "before vert_fill_TS h",G)
+      call hchksum(h*G%H_to_m, "before vert_fill_TS h",G)
     endif
     call vert_fill_TS(h, tv%T, tv%S, kappa_fill, dt_fill, T_f, S_f, G)
     if (CS%debug) then
       call hchksum(tv%T, "after vert_fill_TS tv%T",G)
       call hchksum(tv%S, "after vert_fill_TS tv%S",G)
-      call hchksum(h, "after vert_fill_TS h",G)
+      call hchksum(h*G%H_to_m, "after vert_fill_TS h",G)
     endif
   endif
 
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index 8fdcf34772..fbab853cbc 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -76,7 +76,7 @@ module MOM_vert_friction
 
 use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
 use MOM_diag_mediator, only : diag_ctrl
-use MOM_checksums, only : uchksum, vchksum
+use MOM_checksums, only : uchksum, vchksum, hchksum
 use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE
 use MOM_file_parser, only : get_param, log_version, param_file_type
 use MOM_forcing_type, only : forcing
@@ -893,12 +893,12 @@ subroutine vertvisc_coef(u, v, h, fluxes, visc, dt, G, CS)
 
 
   if (CS%debug) then
-    call uchksum(CS%h_u,"vertvisc_coef h_u",G,haloshift=0)
-    call vchksum(CS%h_v,"vertvisc_coef h_v",G,haloshift=0)
+    call uchksum(CS%h_u*G%H_to_m,"vertvisc_coef h_u",G,haloshift=0)
+    call vchksum(CS%h_v*G%H_to_m,"vertvisc_coef h_v",G,haloshift=0)
     call uchksum(CS%a_u,"vertvisc_coef a_u",G,haloshift=0)
     call vchksum(CS%a_v,"vertvisc_coef a_v",G,haloshift=0)
-    if (allocated(hML_u)) call uchksum(hML_u,"vertvisc_coef hML_u",G,haloshift=0)
-    if (allocated(hML_v)) call vchksum(hML_v,"vertvisc_coef hML_v",G,haloshift=0)
+    if (allocated(hML_u)) call uchksum(hML_u*G%H_to_m,"vertvisc_coef hML_u",G,haloshift=0)
+    if (allocated(hML_v)) call vchksum(hML_v*G%H_to_m,"vertvisc_coef hML_v",G,haloshift=0)
   endif
 
 ! Offer diagnostic fields for averaging.