diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index d0a324f96f..d3ad0a0a92 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -14,6 +14,8 @@ module MOM_dynamics_split_RK2
 use MOM_cpu_clock,         only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
 use MOM_diag_mediator,     only : diag_mediator_init, enable_averages
 use MOM_diag_mediator,     only : disable_averaging, post_data, safe_alloc_ptr
+use MOM_diag_mediator,     only : post_product_u, post_product_sum_u
+use MOM_diag_mediator,     only : post_product_v, post_product_sum_v
 use MOM_diag_mediator,     only : register_diag_field, register_static_field
 use MOM_diag_mediator,     only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids
 use MOM_domains,           only : To_South, To_West, To_All, CGRID_NE, SCALAR_PAIR
@@ -344,36 +346,6 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     v_av, & ! The meridional velocity time-averaged over a time step [L T-1 ~> m s-1].
     h_av    ! The layer thickness time-averaged over a time step [H ~> m or kg m-2].
 
-  ! real, allocatable, dimension(:,:,:) :: &
-    ! hf_PFu, hf_PFv, & ! Pressure force accel. x fract. thickness [L T-2 ~> m s-2].
-    ! hf_CAu, hf_CAv, & ! Coriolis force accel. x fract. thickness [L T-2 ~> m s-2].
-    ! hf_u_BT_accel, hf_v_BT_accel ! barotropic correction accel. x fract. thickness [L T-2 ~> m s-2].
-    ! 3D diagnostics hf_PFu etc. are commented because there is no clarity on proper remapping grid option.
-    ! The code is retained for debugging purposes in the future.
-
-  real, allocatable, dimension(:,:) :: &
-    hf_PFu_2d, hf_PFv_2d, & ! Depth integral of hf_PFu, hf_PFv [L T-2 ~> m s-2].
-    hf_CAu_2d, hf_CAv_2d, & ! Depth integral of hf_CAu, hf_CAv [L T-2 ~> m s-2].
-    hf_u_BT_accel_2d, hf_v_BT_accel_2d ! Depth integral of hf_u_BT_accel, hf_v_BT_accel
-
-  ! Diagnostics for thickness x momentum budget terms
-  real, allocatable, dimension(:,:,:) :: &
-    h_PFu, h_PFv, & ! Pressure force accel. x thickness [H L T-2 ~> m2 s-2].
-    h_CAu, h_CAv, & ! Coriolis force accel. x thickness [H L T-2 ~> m2 s-2].
-    h_u_BT_accel, h_v_BT_accel ! barotropic correction accel. x thickness [H L T-2 ~> m2 s-2].
-
-  ! Diagnostics for layer-sum of thickness x momentum budget terms
-  real, dimension(SZIB_(G),SZJ_(G)) :: &
-    intz_PFu_2d, intz_CAu_2d, intz_u_BT_accel_2d ! [H L T-2 ~> m2 s-2].
-  real, dimension(SZI_(G),SZJB_(G)) :: &
-    intz_PFv_2d, intz_CAv_2d, intz_v_BT_accel_2d ! [H L T-2 ~> m2 s-2].
-
-  ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv],
-  real, allocatable, dimension(:,:,:) :: &
-    PFu_visc_rem, PFv_visc_rem, & ! Pressure force accel. x visc_rem_[uv] [L T-2 ~> m s-2].
-    CAu_visc_rem, CAv_visc_rem, & ! Coriolis force accel. x visc_rem_[uv] [L T-2 ~> m s-2].
-    u_BT_accel_visc_rem, v_BT_accel_visc_rem ! barotropic correction accel. x visc_rem_[uv] [L T-2 ~> m s-2].
-
   real :: dt_pred   ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
 
   logical :: dyn_p_surf
@@ -400,8 +372,6 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     do j=G%jsdB,G%jedB ; do i=G%isd,G%ied   ;  vp(i,j,k) = 0.0 ; enddo ; enddo
     do j=G%jsd,G%jed   ; do i=G%isd,G%ied   ;  hp(i,j,k) = h(i,j,k) ; enddo ; enddo
   enddo
-  if (CS%id_ueffA > 0) ueffA(:,:,:) = 0
-  if (CS%id_veffA > 0) veffA(:,:,:) = 0
 
   ! Update CFL truncation value as function of time
   call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp)
@@ -918,261 +888,76 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
 
   ! Calculate effective areas and post data
   if (CS%id_ueffA > 0) then
-     do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-        if (abs(up(I,j,k)) > 0.) ueffA(I,j,k) = uh(I,j,k)/up(I,j,k)
-     enddo ; enddo ; enddo
-     call post_data(CS%id_ueffA, ueffA, CS%diag)
-  endif
-
-  if (CS%id_veffA > 0) then
-     do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-        if (abs(vp(i,J,k)) > 0.) veffA(i,J,k) = vh(i,J,k)/vp(i,J,k)
-     enddo ; enddo ; enddo
-     call post_data(CS%id_veffA, veffA, CS%diag)
-  endif
-
-
-  ! Diagnostics for terms multiplied by fractional thicknesses
-
-  ! 3D diagnostics hf_PFu etc. are commented because there is no clarity on proper remapping grid option.
-  ! The code is retained for debugging purposes in the future.
-  !if (CS%id_hf_PFu > 0) then
-  !  allocate(hf_PFu(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
-  !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-  !    hf_PFu(I,j,k) = CS%PFu(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k)
-  !  enddo ; enddo ; enddo
-  !  call post_data(CS%id_hf_PFu, hf_PFu, CS%diag)
-  !endif
-  !if (CS%id_hf_PFv > 0) then
-  !  allocate(hf_PFv(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
-  !  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-  !    hf_PFv(i,J,k) = CS%PFv(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k)
-  !  enddo ; enddo ; enddo
-  !  call post_data(CS%id_hf_PFv, hf_PFv, CS%diag)
-  !endif
-  if (CS%id_intz_PFu_2d > 0) then
-    intz_PFu_2d(:,:) = 0.0
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      intz_PFu_2d(I,j) = intz_PFu_2d(I,j) + CS%PFu(I,j,k) * CS%ADp%diag_hu(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_intz_PFu_2d, intz_PFu_2d, CS%diag)
-  endif
-  if (CS%id_intz_PFv_2d > 0) then
-    intz_PFv_2d(:,:) = 0.0
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      intz_PFv_2d(i,J) = intz_PFv_2d(i,J) + CS%PFv(i,J,k) * CS%ADp%diag_hv(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_intz_PFv_2d, intz_PFv_2d, CS%diag)
-  endif
-
-  if (CS%id_hf_PFu_2d > 0) then
-    allocate(hf_PFu_2d(G%IsdB:G%IedB,G%jsd:G%jed), source=0.0)
+    ueffA(:,:,:) = 0
     do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      hf_PFu_2d(I,j) = hf_PFu_2d(I,j) + CS%PFu(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k)
+      if (abs(up(I,j,k)) > 0.) ueffA(I,j,k) = uh(I,j,k) / up(I,j,k)
     enddo ; enddo ; enddo
-    call post_data(CS%id_hf_PFu_2d, hf_PFu_2d, CS%diag)
-    deallocate(hf_PFu_2d)
-  endif
-  if (CS%id_hf_PFv_2d > 0) then
-    allocate(hf_PFv_2d(G%isd:G%ied,G%JsdB:G%JedB), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      hf_PFv_2d(i,J) = hf_PFv_2d(i,J) + CS%PFv(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_hf_PFv_2d, hf_PFv_2d, CS%diag)
-    deallocate(hf_PFv_2d)
+    call post_data(CS%id_ueffA, ueffA, CS%diag)
   endif
 
-  if (CS%id_h_PFu > 0) then
-    allocate(h_PFu(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      h_PFu(I,j,k) = CS%PFu(I,j,k) * CS%ADp%diag_hu(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_h_PFu, h_PFu, CS%diag)
-    deallocate(h_PFu)
-  endif
-  if (CS%id_h_PFv > 0) then
-    allocate(h_PFv(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      h_PFv(i,J,k) = CS%PFv(i,J,k) * CS%ADp%diag_hv(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_h_PFv, h_PFv, CS%diag)
-    deallocate(h_PFv)
-  endif
-
-  !if (CS%id_hf_CAu > 0) then
-  !  allocate(hf_CAu(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
-  !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-  !    hf_CAu(I,j,k) = CS%CAu(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k)
-  !  enddo ; enddo ; enddo
-  !  call post_data(CS%id_hf_CAu, hf_CAu, CS%diag)
-  !endif
-  !if (CS%id_hf_CAv > 0) then
-  !  allocate(hf_CAv(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
-  !  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-  !    hf_CAv(i,J,k) = CS%CAv(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k)
-  !  enddo ; enddo ; enddo
-  !  call post_data(CS%id_hf_CAv, hf_CAv, CS%diag)
-  !endif
-  if (CS%id_intz_CAu_2d > 0) then
-    intz_CAu_2d(:,:) = 0.0
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      intz_CAu_2d(I,j) = intz_CAu_2d(I,j) + CS%CAu(I,j,k) * CS%ADp%diag_hu(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_intz_CAu_2d, intz_CAu_2d, CS%diag)
-  endif
-  if (CS%id_intz_CAv_2d > 0) then
-    intz_CAv_2d(:,:) = 0.0
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      intz_CAv_2d(i,J) = intz_CAv_2d(i,J) + CS%CAv(i,J,k) * CS%ADp%diag_hv(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_intz_CAv_2d, intz_CAv_2d, CS%diag)
-  endif
-
-  if (CS%id_hf_CAu_2d > 0) then
-    allocate(hf_CAu_2d(G%IsdB:G%IedB,G%jsd:G%jed), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      hf_CAu_2d(I,j) = hf_CAu_2d(I,j) + CS%CAu(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_hf_CAu_2d, hf_CAu_2d, CS%diag)
-    deallocate(hf_CAu_2d)
-  endif
-  if (CS%id_hf_CAv_2d > 0) then
-    allocate(hf_CAv_2d(G%isd:G%ied,G%JsdB:G%JedB), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      hf_CAv_2d(i,J) = hf_CAv_2d(i,J) + CS%CAv(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_hf_CAv_2d, hf_CAv_2d, CS%diag)
-    deallocate(hf_CAv_2d)
-  endif
-
-  if (CS%id_h_CAu > 0) then
-    allocate(h_CAu(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      h_CAu(I,j,k) = CS%CAu(I,j,k) * CS%ADp%diag_hu(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_h_CAu, h_CAu, CS%diag)
-    deallocate(h_CAu)
-  endif
-  if (CS%id_h_CAv > 0) then
-    allocate(h_CAv(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      h_CAv(i,J,k) = CS%CAv(i,J,k) * CS%ADp%diag_hv(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_h_CAv, h_CAv, CS%diag)
-    deallocate(h_CAv)
-  endif
-
-  !if (CS%id_hf_u_BT_accel > 0) then
-  !  allocate(hf_u_BT_accel(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
-  !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-  !    hf_u_BT_accel(I,j,k) = CS%u_accel_bt(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k)
-  !  enddo ; enddo ; enddo
-  !  call post_data(CS%id_hf_u_BT_accel, hf_u_BT_accel, CS%diag)
-  !endif
-  !if (CS%id_hf_v_BT_accel > 0) then
-  !  allocate(hf_v_BT_accel(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
-  !  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-  !    hf_v_BT_accel(i,J,k) = CS%v_accel_bt(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k)
-  !  enddo ; enddo ; enddo
-  !  call post_data(CS%id_hf_v_BT_accel, hf_v_BT_accel, CS%diag)
-  !endif
-  if (CS%id_intz_u_BT_accel_2d > 0) then
-    intz_u_BT_accel_2d(:,:) = 0.0
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      intz_u_BT_accel_2d(I,j) = intz_u_BT_accel_2d(I,j) + CS%u_accel_bt(I,j,k) * CS%ADp%diag_hu(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_intz_u_BT_accel_2d, intz_u_BT_accel_2d, CS%diag)
-  endif
-  if (CS%id_intz_v_BT_accel_2d > 0) then
-    intz_v_BT_accel_2d(:,:) = 0.0
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      intz_v_BT_accel_2d(i,J) = intz_v_BT_accel_2d(i,J) + CS%v_accel_bt(i,J,k) * CS%ADp%diag_hv(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_intz_v_BT_accel_2d, intz_v_BT_accel_2d, CS%diag)
-  endif
-
-  if (CS%id_hf_u_BT_accel_2d > 0) then
-    allocate(hf_u_BT_accel_2d(G%IsdB:G%IedB,G%jsd:G%jed), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      hf_u_BT_accel_2d(I,j) = hf_u_BT_accel_2d(I,j) + CS%u_accel_bt(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_hf_u_BT_accel_2d, hf_u_BT_accel_2d, CS%diag)
-    deallocate(hf_u_BT_accel_2d)
-  endif
-  if (CS%id_hf_v_BT_accel_2d > 0) then
-    allocate(hf_v_BT_accel_2d(G%isd:G%ied,G%JsdB:G%JedB), source=0.0)
+  if (CS%id_veffA > 0) then
+    veffA(:,:,:) = 0
     do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      hf_v_BT_accel_2d(i,J) = hf_v_BT_accel_2d(i,J) + CS%v_accel_bt(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k)
+      if (abs(vp(i,J,k)) > 0.) veffA(i,J,k) = vh(i,J,k) / vp(i,J,k)
     enddo ; enddo ; enddo
-    call post_data(CS%id_hf_v_BT_accel_2d, hf_v_BT_accel_2d, CS%diag)
-    deallocate(hf_v_BT_accel_2d)
+    call post_data(CS%id_veffA, veffA, CS%diag)
   endif
 
-  if (CS%id_h_u_BT_accel > 0) then
-    allocate(h_u_BT_accel(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      h_u_BT_accel(I,j,k) = CS%u_accel_bt(I,j,k) * CS%ADp%diag_hu(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_h_u_BT_accel, h_u_BT_accel, CS%diag)
-    deallocate(h_u_BT_accel)
-  endif
-  if (CS%id_h_v_BT_accel > 0) then
-    allocate(h_v_BT_accel(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      h_v_BT_accel(i,J,k) = CS%v_accel_bt(i,J,k) * CS%ADp%diag_hv(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_h_v_BT_accel, h_v_BT_accel, CS%diag)
-    deallocate(h_v_BT_accel)
-  endif
+  ! Diagnostics of the fractional thicknesses times momentum budget terms
+  ! 3D diagnostics hf_PFu etc. are commented because there is no clarity on proper remapping grid option.
+  ! The code is retained for debugging purposes in the future.
+  !if (CS%id_hf_PFu > 0) call post_product_u(CS%id_hf_PFu, CS%PFu, CS%ADp%diag_hfrac_u, G, nz, CS%diag)
+  !if (CS%id_hf_PFv > 0) call post_product_v(CS%id_hf_PFv, CS%PFv, CS%ADp%diag_hfrac_v, G, nz, CS%diag)
+  !if (CS%id_hf_CAu > 0) call post_product_u(CS%id_hf_CAu, CS%CAu, CS%ADp%diag_hfrac_u, G, nz, CS%diag)
+  !if (CS%id_hf_CAv > 0) call post_product_v(CS%id_hf_CAv, CS%CAv, CS%ADp%diag_hfrac_v, G, nz, CS%diag)
+  !if (CS%id_hf_u_BT_accel > 0) &
+  !  call post_product_u(CS%id_hf_u_BT_accel, CS%u_accel_bt, CS%ADp%diag_hfrac_u, G, nz, CS%diag)
+  !if (CS%id_hf_v_BT_accel > 0) &
+  !  call post_product_v(CS%id_hf_v_BT_accel, CS%v_accel_bt, CS%ADp%diag_hfrac_v, G, nz, CS%diag)
+
+  ! Diagnostics for the vertical sum of layer thickness x prssure force accelerations
+  if (CS%id_intz_PFu_2d > 0) call post_product_sum_u(CS%id_intz_PFu_2d, CS%PFu, CS%ADp%diag_hu, G, nz, CS%diag)
+  if (CS%id_intz_PFv_2d > 0) call post_product_sum_v(CS%id_intz_PFv_2d, CS%PFv, CS%ADp%diag_hv, G, nz, CS%diag)
+
+  ! Diagnostics for thickness-weighted vertically averaged prssure force accelerations
+  if (CS%id_hf_PFu_2d > 0) call post_product_sum_u(CS%id_hf_PFu_2d, CS%PFu, CS%ADp%diag_hfrac_u, G, nz, CS%diag)
+  if (CS%id_hf_PFv_2d > 0) call post_product_sum_v(CS%id_hf_PFv_2d, CS%PFv, CS%ADp%diag_hfrac_v, G, nz, CS%diag)
+
+  ! Diagnostics for thickness x prssure force accelerations
+  if (CS%id_h_PFu > 0) call post_product_u(CS%id_h_PFu, CS%PFu, CS%ADp%diag_hu, G, nz, CS%diag)
+  if (CS%id_h_PFv > 0) call post_product_v(CS%id_h_PFv, CS%PFv, CS%ADp%diag_hv, G, nz, CS%diag)
+
+  ! Diagnostics of Coriolis acceleratations
+  if (CS%id_intz_CAu_2d > 0) call post_product_sum_u(CS%id_intz_CAu_2d, CS%CAu, CS%ADp%diag_hu, G, nz, CS%diag)
+  if (CS%id_intz_CAv_2d > 0) call post_product_sum_v(CS%id_intz_CAv_2d, CS%CAv, CS%ADp%diag_hv, G, nz, CS%diag)
+  if (CS%id_hf_CAu_2d > 0) call post_product_sum_u(CS%id_hf_CAu_2d, CS%CAu, CS%ADp%diag_hfrac_u, G, nz, CS%diag)
+  if (CS%id_hf_CAv_2d > 0) call post_product_sum_v(CS%id_hf_CAv_2d, CS%CAv, CS%ADp%diag_hfrac_v, G, nz, CS%diag)
+  if (CS%id_h_CAu > 0) call post_product_u(CS%id_h_CAu, CS%CAu, CS%ADp%diag_hu, G, nz, CS%diag)
+  if (CS%id_h_CAv > 0) call post_product_v(CS%id_h_CAv, CS%CAv, CS%ADp%diag_hv, G, nz, CS%diag)
+
+  ! Diagnostics of barotropic solver acceleratations
+  if (CS%id_intz_u_BT_accel_2d > 0) &
+    call post_product_sum_u(CS%id_intz_u_BT_accel_2d, CS%u_accel_bt, CS%ADp%diag_hu, G, nz, CS%diag)
+  if (CS%id_intz_v_BT_accel_2d > 0) &
+    call post_product_sum_v(CS%id_intz_v_BT_accel_2d, CS%v_accel_bt, CS%ADp%diag_hv, G, nz, CS%diag)
+  if (CS%id_hf_u_BT_accel_2d > 0) &
+    call post_product_sum_u(CS%id_hf_u_BT_accel_2d, CS%u_accel_bt, CS%ADp%diag_hfrac_u, G, nz, CS%diag)
+  if (CS%id_hf_v_BT_accel_2d > 0) &
+    call post_product_sum_v(CS%id_hf_v_BT_accel_2d, CS%v_accel_bt, CS%ADp%diag_hfrac_v, G, nz, CS%diag)
+  if (CS%id_h_u_BT_accel > 0) &
+    call post_product_u(CS%id_h_u_BT_accel, CS%u_accel_bt, CS%ADp%diag_hu, G, nz, CS%diag)
+  if (CS%id_h_v_BT_accel > 0) &
+    call post_product_v(CS%id_h_v_BT_accel, CS%v_accel_bt, CS%ADp%diag_hv, G, nz, CS%diag)
 
-  if (CS%id_PFu_visc_rem > 0) then
-    allocate(PFu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      PFu_visc_rem(I,j,k) = CS%PFu(I,j,k) * CS%ADp%visc_rem_u(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_PFu_visc_rem, PFu_visc_rem, CS%diag)
-    deallocate(PFu_visc_rem)
-  endif
-  if (CS%id_PFv_visc_rem > 0) then
-    allocate(PFv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      PFv_visc_rem(i,J,k) = CS%PFv(i,J,k) * CS%ADp%visc_rem_v(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_PFv_visc_rem, PFv_visc_rem, CS%diag)
-    deallocate(PFv_visc_rem)
-  endif
-  if (CS%id_CAu_visc_rem > 0) then
-    allocate(CAu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      CAu_visc_rem(I,j,k) = CS%CAu(I,j,k) * CS%ADp%visc_rem_u(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_CAu_visc_rem, CAu_visc_rem, CS%diag)
-    deallocate(CAu_visc_rem)
-  endif
-  if (CS%id_CAv_visc_rem > 0) then
-    allocate(CAv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      CAv_visc_rem(i,J,k) = CS%CAv(i,J,k) * CS%ADp%visc_rem_v(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_CAv_visc_rem, CAv_visc_rem, CS%diag)
-    deallocate(CAv_visc_rem)
-  endif
-  if (CS%id_u_BT_accel_visc_rem > 0) then
-    allocate(u_BT_accel_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      u_BT_accel_visc_rem(I,j,k) = CS%u_accel_bt(I,j,k) * CS%ADp%visc_rem_u(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_u_BT_accel_visc_rem, u_BT_accel_visc_rem, CS%diag)
-    deallocate(u_BT_accel_visc_rem)
-  endif
-  if (CS%id_v_BT_accel_visc_rem > 0) then
-    allocate(v_BT_accel_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      v_BT_accel_visc_rem(i,J,k) = CS%v_accel_bt(i,J,k) * CS%ADp%visc_rem_v(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_v_BT_accel_visc_rem, v_BT_accel_visc_rem, CS%diag)
-    deallocate(v_BT_accel_visc_rem)
-  endif
+  ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv],
+  if (CS%id_PFu_visc_rem > 0) call post_product_u(CS%id_PFu_visc_rem, CS%PFu, CS%ADp%visc_rem_u, G, nz, CS%diag)
+  if (CS%id_PFv_visc_rem > 0) call post_product_v(CS%id_PFv_visc_rem, CS%PFv, CS%ADp%visc_rem_v, G, nz, CS%diag)
+  if (CS%id_CAu_visc_rem > 0) call post_product_u(CS%id_CAu_visc_rem, CS%CAu, CS%ADp%visc_rem_u, G, nz, CS%diag)
+  if (CS%id_CAv_visc_rem > 0) call post_product_v(CS%id_CAv_visc_rem, CS%CAv, CS%ADp%visc_rem_v, G, nz, CS%diag)
+  if (CS%id_u_BT_accel_visc_rem > 0) &
+    call post_product_u(CS%id_u_BT_accel_visc_rem, CS%u_accel_bt, CS%ADp%visc_rem_u, G, nz, CS%diag)
+  if (CS%id_v_BT_accel_visc_rem > 0) &
+    call post_product_v(CS%id_v_BT_accel_visc_rem, CS%v_accel_bt, CS%ADp%visc_rem_v, G, nz, CS%diag)
 
   if (CS%debug) then
     call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym)
@@ -1551,7 +1336,6 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
        'Effective V-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, &
        x_cell_method='sum', v_extensive = .true.)
 
-
   !CS%id_hf_PFu = register_diag_field('ocean_model', 'hf_PFu', diag%axesCuL, Time, &
   !    'Fractional Thickness-weighted Zonal Pressure Force Acceleration', &
   !    'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
@@ -1583,13 +1367,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   if (CS%id_hf_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_h_PFu = register_diag_field('ocean_model', 'h_PFu', diag%axesCuL, Time, &
-      'Thickness Multiplied Zonal Pressure Force Acceleration', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      'Thickness Multiplied Zonal Pressure Force Acceleration', &
+      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if(CS%id_h_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_h_PFv = register_diag_field('ocean_model', 'h_PFv', diag%axesCvL, Time, &
-      'Thickness Multiplied Meridional Pressure Force Acceleration', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      'Thickness Multiplied Meridional Pressure Force Acceleration', &
+      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if(CS%id_h_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
   CS%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, Time, &
@@ -1613,13 +1397,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   if (CS%id_hf_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_h_CAu = register_diag_field('ocean_model', 'h_CAu', diag%axesCuL, Time, &
-      'Thickness Multiplied Zonal Coriolis and Advective Acceleration', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      'Thickness Multiplied Zonal Coriolis and Advective Acceleration', &
+      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if(CS%id_h_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_h_CAv = register_diag_field('ocean_model', 'h_CAv', diag%axesCvL, Time, &
-      'Thickness Multiplied Meridional Coriolis and Advective Acceleration', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      'Thickness Multiplied Meridional Coriolis and Advective Acceleration', &
+      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if(CS%id_h_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
   CS%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, Time, &
@@ -1663,13 +1447,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   if (CS%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_h_u_BT_accel = register_diag_field('ocean_model', 'h_u_BT_accel', diag%axesCuL, Time, &
-      'Thickness Multiplied Barotropic Anomaly Zonal Acceleration', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      'Thickness Multiplied Barotropic Anomaly Zonal Acceleration', &
+      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if(CS%id_h_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_h_v_BT_accel = register_diag_field('ocean_model', 'h_v_BT_accel', diag%axesCvL, Time, &
-      'Thickness Multiplied Barotropic Anomaly Meridional Acceleration', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      'Thickness Multiplied Barotropic Anomaly Meridional Acceleration', &
+      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if(CS%id_h_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
   CS%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, Time, &
@@ -1683,30 +1467,30 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   if (CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
   CS%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, Time, &
-      'Zonal Pressure Force Acceleration multiplied by the viscous remnant', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Zonal Pressure Force Acceleration multiplied by the viscous remnant', &
+      'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_PFu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
   CS%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, Time, &
-      'Meridional Pressure Force Acceleration multiplied by the viscous remnant', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Meridional Pressure Force Acceleration multiplied by the viscous remnant', &
+      'm s-2', conversion=US%L_T2_to_m_s2)
   if(CS%id_PFv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, Time, &
-      'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', &
+      'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_CAu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
   CS%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, Time, &
-      'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', &
+      'm s-2', conversion=US%L_T2_to_m_s2)
   if(CS%id_CAv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, Time, &
-      'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', &
+      'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_u_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
   CS%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, Time, &
-      'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', &
+      'm s-2', conversion=US%L_T2_to_m_s2)
   if(CS%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
 
   id_clock_Cor        = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index 439ed242b8..9979ecb5b1 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -9,6 +9,8 @@ module MOM_diagnostics
 use MOM_coupler_types,     only : coupler_type_send_data
 use MOM_density_integrals, only : int_density_dz
 use MOM_diag_mediator,     only : post_data, get_diag_time_end
+use MOM_diag_mediator,     only : post_product_u, post_product_sum_u
+use MOM_diag_mediator,     only : post_product_v, post_product_sum_v
 use MOM_diag_mediator,     only : register_diag_field, register_scalar_field
 use MOM_diag_mediator,     only : register_static_field, diag_register_area_ids
 use MOM_diag_mediator,     only : diag_ctrl, time_type, safe_alloc_ptr
@@ -226,8 +228,6 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
                                                  !! previous call to diagnostics_init.
 
   ! Local variables
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: usq ! squared eastward velocity  [L2 T-2 ~> m2 s-2]
-  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vsq ! squared northward velocity [L2 T-2 ~> m2 s-2]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G))  :: uv  ! u x v at h-points          [L2 T-2 ~> m2 s-2]
 
   integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
@@ -238,12 +238,6 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
   real :: work_2d(SZI_(G),SZJ_(G))         ! A 2-d temporary work array.
   real :: rho_in_situ(SZI_(G))             ! In situ density [R ~> kg m-3]
 
-  real, allocatable, dimension(:,:) :: &
-    hf_du_dt_2d, hf_dv_dt_2d ! z integeral of hf_du_dt, hf_dv_dt [L T-2 ~> m s-2].
-
-  real, allocatable, dimension(:,:,:) :: h_du_dt ! h x dudt [H L T-2 ~> m2 s-2]
-  real, allocatable, dimension(:,:,:) :: h_dv_dt ! h x dvdt [H L T-2 ~> m2 s-2]
-
   ! tmp array for surface properties
   real :: surface_field(SZI_(G),SZJ_(G))
   real :: pressure_1d(SZI_(G)) ! Temporary array for pressure when calling EOS [R L2 T-2 ~> Pa]
@@ -278,70 +272,32 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
     call diag_copy_storage_to_diag(CS%diag, diag_pre_sync)
 
     if (CS%id_h_pre_sync > 0) &
-        call post_data(CS%id_h_pre_sync, diag_pre_sync%h_state, CS%diag, alt_h = diag_pre_sync%h_state)
+        call post_data(CS%id_h_pre_sync, diag_pre_sync%h_state, CS%diag, alt_h=diag_pre_sync%h_state)
 
-    if (CS%id_du_dt>0) call post_data(CS%id_du_dt, CS%du_dt, CS%diag, alt_h = diag_pre_sync%h_state)
+    if (CS%id_du_dt>0) call post_data(CS%id_du_dt, CS%du_dt, CS%diag, alt_h=diag_pre_sync%h_state)
 
-    if (CS%id_dv_dt>0) call post_data(CS%id_dv_dt, CS%dv_dt, CS%diag, alt_h = diag_pre_sync%h_state)
+    if (CS%id_dv_dt>0) call post_data(CS%id_dv_dt, CS%dv_dt, CS%diag, alt_h=diag_pre_sync%h_state)
 
-    if (CS%id_dh_dt>0) call post_data(CS%id_dh_dt, CS%dh_dt, CS%diag, alt_h = diag_pre_sync%h_state)
+    if (CS%id_dh_dt>0) call post_data(CS%id_dh_dt, CS%dh_dt, CS%diag, alt_h=diag_pre_sync%h_state)
 
     !! Diagnostics for terms multiplied by fractional thicknesses
 
     ! 3D diagnostics hf_du(dv)_dt are commented because there is no clarity on proper remapping grid option.
-    ! The code is retained for degugging purposes in the future.
+    ! The code is retained for debugging purposes in the future.
     !if (CS%id_hf_du_dt > 0) then
-    !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-    !    CS%hf_du_dt(I,j,k) = CS%du_dt(I,j,k) * ADp%diag_hfrac_u(I,j,k)
-    !  enddo ; enddo ; enddo
-    !  call post_data(CS%id_hf_du_dt, CS%hf_du_dt, CS%diag, alt_h = diag_pre_sync%h_state)
-    !endif
-
-    !if (CS%id_hf_dv_dt > 0) then
-    !  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-    !    CS%hf_dv_dt(i,J,k) = CS%dv_dt(i,J,k) * ADp%diag_hfrac_v(i,J,k)
-    !  enddo ; enddo ; enddo
-    !  call post_data(CS%id_hf_dv_dt, CS%hf_dv_dt, CS%diag, alt_h = diag_pre_sync%h_state)
-    !endif
-
-    if (CS%id_hf_du_dt_2d > 0) then
-      allocate(hf_du_dt_2d(G%IsdB:G%IedB,G%jsd:G%jed))
-      hf_du_dt_2d(:,:) = 0.0
-      do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-        hf_du_dt_2d(I,j) = hf_du_dt_2d(I,j) + CS%du_dt(I,j,k) * ADp%diag_hfrac_u(I,j,k)
-      enddo ; enddo ; enddo
-      call post_data(CS%id_hf_du_dt_2d, hf_du_dt_2d, CS%diag)
-      deallocate(hf_du_dt_2d)
-    endif
+    !  call post_product_u(CS%id_hf_du_dt, CS%du_dt, ADp%diag_hfrac_u, G, nz, CS%diag, alt_h=diag_pre_sync%h_state)
+    !if (CS%id_hf_dv_dt > 0) &
+    !  call post_product_v(CS%id_hf_dv_dt, CS%dv_dt, ADp%diag_hfrac_v, G, nz, CS%diag, alt_h=diag_pre_sync%h_state)
 
-    if (CS%id_hf_dv_dt_2d > 0) then
-      allocate(hf_dv_dt_2d(G%isd:G%ied,G%JsdB:G%JedB))
-      hf_dv_dt_2d(:,:) = 0.0
-      do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-        hf_dv_dt_2d(i,J) = hf_dv_dt_2d(i,J) + CS%dv_dt(i,J,k) * ADp%diag_hfrac_v(i,J,k)
-      enddo ; enddo ; enddo
-      call post_data(CS%id_hf_dv_dt_2d, hf_dv_dt_2d, CS%diag)
-      deallocate(hf_dv_dt_2d)
-    endif
+    if (CS%id_hf_du_dt_2d > 0) &
+      call post_product_sum_u(CS%id_hf_du_dt_2d, CS%du_dt, ADp%diag_hfrac_u, G, nz, CS%diag)
+    if (CS%id_hf_dv_dt_2d > 0) &
+      call post_product_sum_v(CS%id_hf_dv_dt_2d, CS%dv_dt, ADp%diag_hfrac_v, G, nz, CS%diag)
 
-    if (CS%id_h_du_dt > 0) then
-      allocate(h_du_dt(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
-      h_du_dt(:,:,:) = 0.0
-      do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-        h_du_dt(I,j,k) = CS%du_dt(I,j,k) * ADp%diag_hu(I,j,k)
-      enddo ; enddo ; enddo
-      call post_data(CS%id_h_du_dt, h_du_dt, CS%diag)
-      deallocate(h_du_dt)
-    endif
-    if (CS%id_h_dv_dt > 0) then
-      allocate(h_dv_dt(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
-      h_dv_dt(:,:,:) = 0.0
-      do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-        h_dv_dt(i,J,k) = CS%dv_dt(i,J,k) * ADp%diag_hv(i,J,k)
-      enddo ; enddo ; enddo
-      call post_data(CS%id_h_dv_dt, h_dv_dt, CS%diag)
-      deallocate(h_dv_dt)
-    endif
+    if (CS%id_h_du_dt > 0) &
+      call post_product_u(CS%id_h_du_dt, CS%du_dt, ADp%diag_hu, G, nz, CS%diag)
+    if (CS%id_h_dv_dt > 0) &
+      call post_product_v(CS%id_h_dv_dt, CS%dv_dt, ADp%diag_hv, G, nz, CS%diag)
 
     call diag_restore_grids(CS%diag)
 
@@ -362,24 +318,14 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
 
   if (CS%id_h > 0) call post_data(CS%id_h, h, CS%diag)
 
-  if (CS%id_usq > 0) then
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      usq(I,j,k) = u(I,j,k) * u(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_usq, usq, CS%diag)
-  endif
+  if (CS%id_usq > 0) call post_product_u(CS%id_usq, u, u, G, nz, CS%diag)
 
-  if (CS%id_vsq > 0) then
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      vsq(i,J,k) = v(i,J,k) * v(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_vsq, vsq, CS%diag)
-  endif
+  if (CS%id_vsq > 0) call post_product_v(CS%id_vsq, v, v, G, nz, CS%diag)
 
   if (CS%id_uv > 0) then
     do k=1,nz ; do j=js,je ; do i=is,ie
       uv(i,j,k) = (0.5*(u(I-1,j,k) + u(I,j,k))) * &
-                (0.5*(v(i,J-1,k) + v(i,J,k)))
+                  (0.5*(v(i,J-1,k) + v(i,J,k)))
     enddo ; enddo ; enddo
     call post_data(CS%id_uv, uv, CS%diag)
   endif
diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90
index 374f54548e..eb24c994f8 100644
--- a/src/framework/MOM_diag_mediator.F90
+++ b/src/framework/MOM_diag_mediator.F90
@@ -39,6 +39,7 @@ module MOM_diag_mediator
 #define MAX_DSAMP_LEV 2
 
 public set_axes_info, post_data, register_diag_field, time_type
+public post_product_u, post_product_sum_u, post_product_v, post_product_sum_v
 public set_masks_for_axes
 public post_data_1d_k
 public safe_alloc_ptr, safe_alloc_alloc
@@ -1802,6 +1803,108 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask)
 
 end subroutine post_data_3d_low
 
+!> Calculate and write out diagnostics that are the product of two 3-d arrays at u-points
+subroutine post_product_u(id, u_a, u_b, G, nz, diag, mask, alt_h)
+  integer,                  intent(in) :: id   !< The ID for this diagnostic
+  type(ocean_grid_type),    intent(in) :: G    !< ocean grid structure
+  integer,                  intent(in) :: nz   !< The size of the arrays in the vertical
+  real, dimension(G%IsdB:G%IedB, G%jsd:G%jed, nz), &
+                            intent(in) :: u_a  !< The first u-point array in arbitrary units [A]
+  real, dimension(G%IsdB:G%IedB, G%jsd:G%jed, nz), &
+                            intent(in) :: u_b  !< The second u-point array in arbitrary units [B]
+  type(diag_ctrl),          intent(in) :: diag !< regulates diagnostic output
+  real,           optional, intent(in) :: mask(:,:,:)  !< If present, use this real array as the data mask [nondim]
+  real,   target, optional, intent(in) :: alt_h(:,:,:) !< An alternate thickness to use for vertically
+                                               !! remapping this diagnostic [H ~> m or kg m-2]
+
+  ! Local variables
+  real, dimension(G%IsdB:G%IedB, G%jsd:G%jed, nz) :: u_prod ! The product of u_a and u_b [A B]
+  integer :: i, j, k
+
+  if (id <= 0) return
+
+  do k=1,nz ; do j=G%jsc,G%jec ; do I=G%IscB,G%IecB
+    u_prod(I,j,k) = u_a(I,j,k) * u_b(I,j,k)
+  enddo ; enddo ; enddo
+  call post_data(id, u_prod, diag, mask=mask, alt_h=alt_h)
+
+end subroutine post_product_u
+
+!> Calculate and write out diagnostics that are the vertical sum of the product of two 3-d arrays at u-points
+subroutine post_product_sum_u(id, u_a, u_b, G, nz, diag)
+  integer,                  intent(in) :: id   !< The ID for this diagnostic
+  type(ocean_grid_type),    intent(in) :: G    !< ocean grid structure
+  integer,                  intent(in) :: nz   !< The size of the arrays in the vertical
+  real, dimension(G%IsdB:G%IedB, G%jsd:G%jed, nz), &
+                            intent(in) :: u_a  !< The first u-point array in arbitrary units [A]
+  real, dimension(G%IsdB:G%IedB, G%jsd:G%jed, nz), &
+                            intent(in) :: u_b  !< The second u-point array in arbitrary units [B]
+  type(diag_ctrl),          intent(in) :: diag !< regulates diagnostic output
+
+  real, dimension(G%IsdB:G%IedB, G%jsd:G%jed) :: u_sum  ! The vertical sum of the product of u_a and u_b [A B]
+  integer :: i, j, k
+
+  if (id <= 0) return
+
+  u_sum(:,:) = 0.0
+  do k=1,nz ; do j=G%jsc,G%jec ; do I=G%IscB,G%IecB
+    u_sum(I,j) = u_sum(I,j) + u_a(I,j,k) * u_b(I,j,k)
+  enddo ; enddo ; enddo
+  call post_data(id, u_sum, diag)
+
+end subroutine post_product_sum_u
+
+!> Calculate and write out diagnostics that are the product of two 3-d arrays at v-points
+subroutine post_product_v(id, v_a, v_b, G, nz, diag, mask, alt_h)
+  integer,                  intent(in) :: id   !< The ID for this diagnostic
+  type(ocean_grid_type),    intent(in) :: G    !< ocean grid structure
+  integer,                  intent(in) :: nz   !< The size of the arrays in the vertical
+  real, dimension(G%isd:G%ied, G%JsdB:G%JedB, nz), &
+                            intent(in) :: v_a  !< The first v-point array in arbitrary units [A]
+  real, dimension(G%isd:G%ied, G%JsdB:G%JedB, nz), &
+                            intent(in) :: v_b  !< The second v-point array in arbitrary units [B]
+  type(diag_ctrl),          intent(in) :: diag !< regulates diagnostic output
+  real,           optional, intent(in) :: mask(:,:,:)  !< If present, use this real array as the data mask [nondim]
+  real,   target, optional, intent(in) :: alt_h(:,:,:) !< An alternate thickness to use for vertically
+                                               !! remapping this diagnostic [H ~> m or kg m-2]
+
+  ! Local variables
+  real, dimension(G%isd:G%ied, G%JsdB:G%JedB, nz) :: v_prod ! The product of v_a and v_b [A B]
+  integer :: i, j, k
+
+  if (id <= 0) return
+
+  do k=1,nz ; do J=G%JscB,G%JecB ; do i=G%isc,G%iec
+    v_prod(i,J,k) = v_a(i,J,k) * v_b(i,J,k)
+  enddo ; enddo ; enddo
+  call post_data(id, v_prod, diag, mask=mask, alt_h=alt_h)
+
+end subroutine post_product_v
+
+!> Calculate and write out diagnostics that are the vertical sum of the product of two 3-d arrays at v-points
+subroutine post_product_sum_v(id, v_a, v_b, G, nz, diag)
+  integer,                  intent(in) :: id   !< The ID for this diagnostic
+  type(ocean_grid_type),    intent(in) :: G    !< ocean grid structure
+  integer,                  intent(in) :: nz   !< The size of the arrays in the vertical
+  real, dimension(G%isd:G%ied, G%JsdB:G%JedB, nz), &
+                            intent(in) :: v_a  !< The first v-point array in arbitrary units [A]
+  real, dimension(G%isd:G%ied, G%JsdB:G%JedB, nz), &
+                            intent(in) :: v_b  !< The second v-point array in arbitrary units [B]
+  type(diag_ctrl),          intent(in) :: diag !< regulates diagnostic output
+
+  real, dimension(G%isd:G%ied, G%JsdB:G%JedB) :: v_sum ! The vertical sum of the product of v_a and v_b [A B]
+  integer :: i, j, k
+
+  if (id <= 0) return
+
+  v_sum(:,:) = 0.0
+  do k=1,nz ; do J=G%JscB,G%JecB ; do i=G%isc,G%iec
+    v_sum(i,J) = v_sum(i,J) + v_a(i,J,k) * v_b(i,J,k)
+  enddo ; enddo ; enddo
+  call post_data(id, v_sum, diag)
+
+end subroutine post_product_sum_v
+
 !> Post the horizontally area-averaged diagnostic
 subroutine post_xy_average(diag_cs, diag, field)
   type(diag_type),   intent(in) :: diag !< This diagnostic
diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index 4381af9d84..6a9b49683c 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -6,6 +6,8 @@ module MOM_hor_visc
 use MOM_checksums,             only : hchksum, Bchksum
 use MOM_coms,                  only : min_across_PEs
 use MOM_diag_mediator,         only : post_data, register_diag_field, safe_alloc_ptr
+use MOM_diag_mediator,         only : post_product_u, post_product_sum_u
+use MOM_diag_mediator,         only : post_product_v, post_product_sum_v
 use MOM_diag_mediator,         only : diag_ctrl, time_type
 use MOM_domains,               only : pass_var, CORNER, pass_vector, AGRID, BGRID_NE
 use MOM_error_handler,         only : MOM_error, FATAL, WARNING, is_root_pe
@@ -1662,90 +1664,30 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     call post_data(CS%id_FrictWorkIntz, FrictWorkIntz, CS%diag)
   endif
 
-  ! Diagnostics for terms multiplied by fractional thicknesses
-
-  ! 3D diagnostics hf_diffu(diffv) are commented because there is no clarity on proper remapping grid option.
-  ! The code is retained for degugging purposes in the future.
-  !if (present(ADp) .and. (CS%id_hf_diffu > 0)) then
-  !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-  !    CS%hf_diffu(I,j,k) = diffu(I,j,k) * ADp%diag_hfrac_u(I,j,k)
-  !  enddo ; enddo ; enddo
-  !  call post_data(CS%id_hf_diffu, CS%hf_diffu, CS%diag)
-  !endif
-  !if (present(ADp) .and. (CS%id_hf_diffv > 0)) then
-  !  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-  !    CS%hf_diffv(i,J,k) = diffv(i,J,k) * ADp%diag_hfrac_v(i,J,k)
-  !  enddo ; enddo ; enddo
-  !  call post_data(CS%id_hf_diffv, CS%hf_diffv, CS%diag)
-  !endif
-
   if (present(ADp)) then
-    if (CS%id_hf_diffu_2d > 0) then
-      hf_diffu_2d(:,:) = 0.0
-      do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-        hf_diffu_2d(I,j) = hf_diffu_2d(I,j) + diffu(I,j,k) * ADp%diag_hfrac_u(I,j,k)
-      enddo ; enddo ; enddo
-      call post_data(CS%id_hf_diffu_2d, hf_diffu_2d, CS%diag)
-    endif
-
-    if (CS%id_hf_diffv_2d > 0) then
-      hf_diffv_2d(:,:) = 0.0
-      do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-        hf_diffv_2d(i,J) = hf_diffv_2d(i,J) + diffv(i,J,k) * ADp%diag_hfrac_v(i,J,k)
-      enddo ; enddo ; enddo
-      call post_data(CS%id_hf_diffv_2d, hf_diffv_2d, CS%diag)
-    endif
-
-    if (CS%id_intz_diffu_2d > 0) then
-      intz_diffu_2d(:,:) = 0.0
-      do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-        intz_diffu_2d(I,j) = intz_diffu_2d(I,j) + diffu(I,j,k) * ADp%diag_hu(I,j,k)
-      enddo ; enddo ; enddo
-      call post_data(CS%id_intz_diffu_2d, intz_diffu_2d, CS%diag)
-    endif
-
-    if (CS%id_intz_diffv_2d > 0) then
-      intz_diffv_2d(:,:) = 0.0
-      do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-        intz_diffv_2d(i,J) = intz_diffv_2d(i,J) + diffv(i,J,k) * ADp%diag_hv(i,J,k)
-      enddo ; enddo ; enddo
-      call post_data(CS%id_intz_diffv_2d, intz_diffv_2d, CS%diag)
-    endif
+    ! Diagnostics of the fractional thicknesses times momentum budget terms
+    ! 3D diagnostics of hf_diffu(diffv) are commented because there is no clarity on proper remapping grid option.
+    ! The code is retained for debugging purposes in the future.
+    !if (CS%id_hf_diffu > 0) call post_product_u(CS%id_hf_diffu, diffu, ADp%diag_hfrac_u, G, nz, CS%diag)
+    !if (CS%id_hf_diffv > 0) call post_product_v(CS%id_hf_diffv, diffv, ADp%diag_hfrac_v, G, nz, CS%diag)
+
+    ! Diagnostics for thickness-weighted vertically averaged momentum budget terms
+    if (CS%id_hf_diffu_2d > 0) call post_product_sum_u(CS%id_hf_diffu_2d, diffu, ADp%diag_hfrac_u, G, nz, CS%diag)
+    if (CS%id_hf_diffv_2d > 0) call post_product_sum_v(CS%id_hf_diffv_2d, diffv, ADp%diag_hfrac_v, G, nz, CS%diag)
+
+    ! Diagnostics for the vertical sum of layer thickness x momentum budget terms
+    if (CS%id_intz_diffu_2d > 0) call post_product_sum_u(CS%id_intz_diffu_2d, diffu, ADp%diag_hu, G, nz, CS%diag)
+    if (CS%id_intz_diffv_2d > 0) call post_product_sum_v(CS%id_intz_diffv_2d, diffv, ADp%diag_hv, G, nz, CS%diag)
+
+    ! Diagnostics for thickness x momentum budget terms
+    if (CS%id_h_diffu > 0) call post_product_u(CS%id_h_diffu, diffu, ADp%diag_hu, G, nz, CS%diag)
+    if (CS%id_h_diffv > 0) call post_product_v(CS%id_h_diffv, diffv, ADp%diag_hv, G, nz, CS%diag)
+
+    ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv],
+    if (CS%id_diffu_visc_rem > 0) call post_product_u(CS%id_diffu_visc_rem, diffu, ADp%visc_rem_u, G, nz, CS%diag)
+    if (CS%id_diffv_visc_rem > 0) call post_product_v(CS%id_diffv_visc_rem, diffv, ADp%visc_rem_v, G, nz, CS%diag)
   endif
 
-  if (present(ADp) .and. (CS%id_h_diffu > 0)) then
-    allocate(h_diffu(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      h_diffu(I,j,k) = diffu(I,j,k) * ADp%diag_hu(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_h_diffu, h_diffu, CS%diag)
-    deallocate(h_diffu)
-  endif
-  if (present(ADp) .and. (CS%id_h_diffv > 0)) then
-    allocate(h_diffv(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      h_diffv(i,J,k) = diffv(i,J,k) * ADp%diag_hv(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_h_diffv, h_diffv, CS%diag)
-    deallocate(h_diffv)
-  endif
-
-  if (present(ADp) .and. (CS%id_diffu_visc_rem > 0)) then
-    allocate(diffu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      diffu_visc_rem(I,j,k) = diffu(I,j,k) * ADp%visc_rem_u(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_diffu_visc_rem, diffu_visc_rem, CS%diag)
-    deallocate(diffu_visc_rem)
-  endif
-  if (present(ADp) .and. (CS%id_diffv_visc_rem > 0)) then
-    allocate(diffv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      diffv_visc_rem(i,J,k) = diffv(i,J,k) * ADp%visc_rem_v(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_diffv_visc_rem, diffv_visc_rem, CS%diag)
-    deallocate(diffv_visc_rem)
-  endif
 end subroutine horizontal_viscosity
 
 !> Allocates space for and calculates static variables used by horizontal_viscosity().
@@ -2467,15 +2409,15 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
   endif
 
   CS%id_h_diffu = register_diag_field('ocean_model', 'h_diffu', diag%axesCuL, Time, &
-      'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', &
+      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if ((CS%id_h_diffu > 0) .and. (present(ADp))) then
     call safe_alloc_ptr(ADp%diag_hu,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
   endif
 
   CS%id_h_diffv = register_diag_field('ocean_model', 'h_diffv', diag%axesCvL, Time, &
-      'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', &
+      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if ((CS%id_h_diffv > 0) .and. (present(ADp))) then
     call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
   endif
@@ -2495,15 +2437,15 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
   endif
 
   CS%id_diffu_visc_rem = register_diag_field('ocean_model', 'diffu_visc_rem', diag%axesCuL, Time, &
-      'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant', &
+      'm s-2', conversion=US%L_T2_to_m_s2)
   if ((CS%id_diffu_visc_rem > 0) .and. (present(ADp))) then
     call safe_alloc_ptr(ADp%visc_rem_u,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
   endif
 
   CS%id_diffv_visc_rem = register_diag_field('ocean_model', 'diffv_visc_rem', diag%axesCvL, Time, &
-      'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant', &
+      'm s-2', conversion=US%L_T2_to_m_s2)
   if ((CS%id_diffv_visc_rem > 0) .and. (present(ADp))) then
     call safe_alloc_ptr(ADp%visc_rem_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
   endif
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index c0ea3aff53..77ec87b230 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -15,6 +15,7 @@ module MOM_diabatic_driver
 use MOM_diabatic_aux,        only : triDiagTS_Eulerian, find_uv_at_h, diagnoseMLDbyDensityDifference
 use MOM_diabatic_aux,        only : applyBoundaryFluxesInOut, diagnoseMLDbyEnergy, set_pen_shortwave
 use MOM_diag_mediator,       only : post_data, register_diag_field, safe_alloc_ptr
+use MOM_diag_mediator,       only : post_product_sum_u, post_product_sum_v
 use MOM_diag_mediator,       only : diag_ctrl, time_type, diag_update_remap_grids
 use MOM_diag_mediator,       only : diag_ctrl, query_averaging_enabled, enable_averages, disable_averaging
 use MOM_diag_mediator,       only : diag_grid_storage, diag_grid_storage_init, diag_grid_storage_end
@@ -2499,24 +2500,11 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
   if (CS%id_Sdif > 0) call post_data(CS%id_Sdif, Sdif_flx, CS%diag)
   if (CS%id_Sadv > 0) call post_data(CS%id_Sadv, Sadv_flx, CS%diag)
 
-  !! Diagnostics for terms multiplied by fractional thicknesses
-  if (CS%id_hf_dudt_dia_2d > 0) then
-    allocate(hf_dudt_dia_2d(G%IsdB:G%IedB,G%jsd:G%jed), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      hf_dudt_dia_2d(I,j) = hf_dudt_dia_2d(I,j) + ADp%du_dt_dia(I,j,k) * ADp%diag_hfrac_u(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_hf_dudt_dia_2d, hf_dudt_dia_2d, CS%diag)
-    deallocate(hf_dudt_dia_2d)
-  endif
-
-  if (CS%id_hf_dvdt_dia_2d > 0) then
-    allocate(hf_dvdt_dia_2d(G%isd:G%ied,G%JsdB:G%JedB), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      hf_dvdt_dia_2d(i,J) = hf_dvdt_dia_2d(i,J) + ADp%dv_dt_dia(i,J,k) * ADp%diag_hfrac_v(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_hf_dvdt_dia_2d, hf_dvdt_dia_2d, CS%diag)
-    deallocate(hf_dvdt_dia_2d)
-  endif
+  ! Diagnostics for thickness-weighted vertically averaged diapycnal accelerations
+  if (CS%id_hf_dudt_dia_2d > 0) &
+    call post_product_sum_u(CS%id_hf_dudt_dia_2d, ADp%du_dt_dia, ADp%diag_hfrac_u, G, nz, CS%diag)
+  if (CS%id_hf_dvdt_dia_2d > 0) &
+    call post_product_sum_v(CS%id_hf_dvdt_dia_2d, ADp%dv_dt_dia, ADp%diag_hfrac_v, G, nz, CS%diag)
 
   call disable_averaging(CS%diag)
 
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index b332346c6c..adac9e83f4 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -4,23 +4,25 @@ module MOM_vert_friction
 ! This file is part of MOM6. See LICENSE.md for the license.
 use MOM_domains,       only : pass_var, To_All, Omit_corners
 use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
+use MOM_diag_mediator, only : post_product_u, post_product_sum_u
+use MOM_diag_mediator, only : post_product_v, post_product_sum_v
 use MOM_diag_mediator, only : diag_ctrl
-use MOM_debugging, only : uvchksum, hchksum
+use MOM_debugging,     only : uvchksum, 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 : mech_forcing
-use MOM_get_input, only : directories
-use MOM_grid, only : ocean_grid_type
+use MOM_file_parser,   only : get_param, log_version, param_file_type
+use MOM_forcing_type,  only : mech_forcing
+use MOM_get_input,     only : directories
+use MOM_grid,          only : ocean_grid_type
 use MOM_open_boundary, only : ocean_OBC_type, OBC_SIMPLE, OBC_NONE, OBC_DIRECTION_E
 use MOM_open_boundary, only : OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S
-use MOM_PointAccel, only : write_u_accel, write_v_accel, PointAccel_init
-use MOM_PointAccel, only : PointAccel_CS
-use MOM_time_manager, only : time_type, time_type_to_real, operator(-)
-use MOM_unit_scaling, only : unit_scale_type
-use MOM_variables, only : thermo_var_ptrs, vertvisc_type
-use MOM_variables, only : cont_diag_ptrs, accel_diag_ptrs
-use MOM_variables, only : ocean_internal_state
-use MOM_verticalGrid, only : verticalGrid_type
+use MOM_PointAccel,    only : write_u_accel, write_v_accel, PointAccel_init
+use MOM_PointAccel,    only : PointAccel_CS
+use MOM_time_manager,  only : time_type, time_type_to_real, operator(-)
+use MOM_unit_scaling,  only : unit_scale_type
+use MOM_variables,     only : thermo_var_ptrs, vertvisc_type
+use MOM_variables,     only : cont_diag_ptrs, accel_diag_ptrs
+use MOM_variables,     only : ocean_internal_state
+use MOM_verticalGrid,  only : verticalGrid_type
 use MOM_wave_interface, only : wave_parameters_CS
 implicit none ; private
 
@@ -136,11 +138,6 @@ module MOM_vert_friction
   type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() !< A pointer to the control structure
                               !! for recording accelerations leading to velocity truncations
 
-  ! real, pointer :: hf_du_dt_visc(:,:,:)  => NULL() ! Zonal friction accel. x fract. thickness [L T-2 ~> m s-2].
-  ! real, pointer :: hf_dv_dt_visc(:,:,:)  => NULL() ! Merdional friction accel. x fract. thickness [L T-2 ~> m s-2].
-  ! 3D diagnostics hf_du(dv)_dt_visc are commented because there is no clarity on proper remapping grid option.
-  ! The code is retained for degugging purposes in the future.
-
 end type vertvisc_CS
 
 contains
@@ -217,16 +214,6 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
   real :: surface_stress(SZIB_(G))! The same as stress, unless the wind stress
                            ! stress is applied as a body force [H L T-1 ~> m2 s-1 or kg m-1 s-1].
 
-  real, allocatable, dimension(:,:) :: hf_du_dt_visc_2d ! Depth sum of hf_du_dt_visc [L T-2 ~> m s-2]
-  real, allocatable, dimension(:,:) :: hf_dv_dt_visc_2d ! Depth sum of hf_dv_dt_visc [L T-2 ~> m s-2]
-
-  real, allocatable, dimension(:,:,:) :: h_du_dt_visc ! h x du_dt_visc [H L T-2 ~> m2 s-2]
-  real, allocatable, dimension(:,:,:) :: h_dv_dt_visc ! h x dv_dt_visc [H L T-2 ~> m2 s-2]
-  real, allocatable, dimension(:,:,:) :: h_du_dt_str ! h x du_dt_str [H L T-2 ~> m2 s-2]
-  real, allocatable, dimension(:,:,:) :: h_dv_dt_str ! h x dv_dt_str [H L T-2 ~> m2 s-2]
-  real, allocatable, dimension(:,:,:) :: du_dt_str_visc_rem ! du_dt_str x visc_rem_u [L T-2 ~> m s-2]
-  real, allocatable, dimension(:,:,:) :: dv_dt_str_visc_rem ! dv_dt_str x visc_rem_v [L T-2 ~> m s-2]
-
   logical :: do_i(SZIB_(G))
   logical :: DoStokesMixing
 
@@ -525,92 +512,36 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
   if (CS%id_dv_dt_str > 0) &
     call post_data(CS%id_dv_dt_str, ADp%dv_dt_str, CS%diag)
 
-  ! Diagnostics for terms multiplied by fractional thicknesses
-
-  ! 3D diagnostics hf_du(dv)_dt_visc are commented because there is no clarity on proper remapping grid option.
-  ! The code is retained for degugging purposes in the future.
-  !if (CS%id_hf_du_dt_visc > 0) then
-  !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-  !    CS%hf_du_dt_visc(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%diag_hfrac_u(I,j,k)
-  !  enddo ; enddo ; enddo
-  !  call post_data(CS%id_hf_du_dt_visc, CS%hf_du_dt_visc, CS%diag)
-  !endif
-  !if (CS%id_hf_dv_dt_visc > 0) then
-  !  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-  !    CS%hf_dv_dt_visc(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%diag_hfrac_v(i,J,k)
-  !  enddo ; enddo ; enddo
-  !  call post_data(CS%id_hf_dv_dt_visc, CS%hf_dv_dt_visc, CS%diag)
-  !endif
-  if (CS%id_hf_du_dt_visc_2d > 0) then
-    allocate(hf_du_dt_visc_2d(G%IsdB:G%IedB,G%jsd:G%jed), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      hf_du_dt_visc_2d(I,j) = hf_du_dt_visc_2d(I,j) + ADp%du_dt_visc(I,j,k) * ADp%diag_hfrac_u(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_hf_du_dt_visc_2d, hf_du_dt_visc_2d, CS%diag)
-    deallocate(hf_du_dt_visc_2d)
-  endif
-  if (CS%id_hf_dv_dt_visc_2d > 0) then
-    allocate(hf_dv_dt_visc_2d(G%isd:G%ied,G%JsdB:G%JedB), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      hf_dv_dt_visc_2d(i,J) = hf_dv_dt_visc_2d(i,J) + ADp%dv_dt_visc(i,J,k) * ADp%diag_hfrac_v(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_hf_dv_dt_visc_2d, hf_dv_dt_visc_2d, CS%diag)
-    deallocate(hf_dv_dt_visc_2d)
+  if (associated(ADp%du_dt_visc) .and. associated(ADp%du_dt_visc)) then
+    ! Diagnostics of the fractional thicknesses times momentum budget terms
+    ! 3D diagnostics of hf_du(dv)_dt_visc are commented because there is no clarity on proper remapping grid option.
+    ! The code is retained for debugging purposes in the future.
+    !if (CS%id_hf_du_dt_visc > 0) &
+    !  call post_product_u(CS%id_hf_du_dt_visc, ADp%du_dt_visc, ADp%diag_hfrac_u, G, nz, CS%diag)
+    !if (CS%id_hf_dv_dt_visc > 0) &
+    !  call post_product_v(CS%id_hf_dv_dt_visc, ADp%dv_dt_visc, ADp%diag_hfrac_v, G, nz, CS%diag)
+
+    ! Diagnostics for thickness-weighted vertically averaged viscous accelerations
+    if (CS%id_hf_du_dt_visc_2d > 0) &
+      call post_product_sum_u(CS%id_hf_du_dt_visc_2d, ADp%du_dt_visc, ADp%diag_hfrac_u, G, nz, CS%diag)
+    if (CS%id_hf_dv_dt_visc_2d > 0) &
+      call post_product_sum_v(CS%id_hf_dv_dt_visc_2d, ADp%dv_dt_visc, ADp%diag_hfrac_v, G, nz, CS%diag)
+
+    ! Diagnostics for thickness x viscous accelerations
+    if (CS%id_h_du_dt_visc > 0) call post_product_u(CS%id_h_du_dt_visc, ADp%du_dt_visc, ADp%diag_hu, G, nz, CS%diag)
+    if (CS%id_h_dv_dt_visc > 0) call post_product_v(CS%id_h_dv_dt_visc, ADp%dv_dt_visc, ADp%diag_hv, G, nz, CS%diag)
   endif
 
-  if (CS%id_h_du_dt_visc > 0) then
-    allocate(h_du_dt_visc(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      h_du_dt_visc(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%diag_hu(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_h_du_dt_visc, h_du_dt_visc, CS%diag)
-    deallocate(h_du_dt_visc)
-  endif
-  if (CS%id_h_dv_dt_visc > 0) then
-    allocate(h_dv_dt_visc(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      h_dv_dt_visc(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%diag_hv(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_h_dv_dt_visc, h_dv_dt_visc, CS%diag)
-    deallocate(h_dv_dt_visc)
-  endif
-
-  if (CS%id_h_du_dt_str > 0) then
-    allocate(h_du_dt_str(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
-    h_du_dt_str(:,:,:) = 0.0
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      h_du_dt_str(I,j,k) = ADp%du_dt_str(I,j,k) * ADp%diag_hu(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_h_du_dt_str, h_du_dt_str, CS%diag)
-    deallocate(h_du_dt_str)
-  endif
-  if (CS%id_h_dv_dt_str > 0) then
-    allocate(h_dv_dt_str(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
-    h_dv_dt_str(:,:,:) = 0.0
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      h_dv_dt_str(i,J,k) = ADp%dv_dt_str(i,J,k) * ADp%diag_hv(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_h_dv_dt_str, h_dv_dt_str, CS%diag)
-    deallocate(h_dv_dt_str)
-  endif
+  if (associated(ADp%du_dt_str) .and.  associated(ADp%dv_dt_str)) then
+    ! Diagnostics for thickness x wind stress acclerations
+    if (CS%id_h_du_dt_str > 0) call post_product_u(CS%id_h_du_dt_str, ADp%du_dt_str, ADp%diag_hu, G, nz, CS%diag)
+    if (CS%id_h_dv_dt_str > 0) call post_product_v(CS%id_h_dv_dt_str, ADp%dv_dt_str, ADp%diag_hv, G, nz, CS%diag)
 
-  if (CS%id_du_dt_str_visc_rem > 0) then
-    allocate(du_dt_str_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
-    du_dt_str_visc_rem(:,:,:) = 0.0
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      du_dt_str_visc_rem(I,j,k) = ADp%du_dt_str(I,j,k) * ADp%visc_rem_u(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_du_dt_str_visc_rem, du_dt_str_visc_rem, CS%diag)
-    deallocate(du_dt_str_visc_rem)
-  endif
-  if (CS%id_dv_dt_str_visc_rem > 0) then
-    allocate(dv_dt_str_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
-    dv_dt_str_visc_rem(:,:,:) = 0.0
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      dv_dt_str_visc_rem(i,J,k) = ADp%dv_dt_str(i,J,k) * ADp%visc_rem_v(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(CS%id_dv_dt_str_visc_rem, dv_dt_str_visc_rem, CS%diag)
-    deallocate(dv_dt_str_visc_rem)
+    ! Diagnostics for wind stress accelerations multiplied by visc_rem_[uv],
+    if (CS%id_du_dt_str_visc_rem > 0) &
+      call post_product_u(CS%id_du_dt_str_visc_rem, ADp%du_dt_str, ADp%visc_rem_u, G, nz, CS%diag)
+    if (CS%id_dv_dt_str_visc_rem > 0) &
+      call post_product_v(CS%id_dv_dt_str_visc_rem, ADp%dv_dt_str, ADp%visc_rem_v, G, nz, CS%diag)
   endif
 
 end subroutine vertvisc
@@ -1922,7 +1853,6 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
   !    'Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity', &
   !    'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
   !if (CS%id_hf_du_dt_visc > 0) then
-  !  call safe_alloc_ptr(CS%hf_du_dt_visc,IsdB,IedB,jsd,jed,nz)
   !  call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
   !  call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
   !endif
@@ -1931,7 +1861,6 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
   !    'Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity', &
   !    'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
   !if (CS%id_hf_dv_dt_visc > 0) then
-  !  call safe_alloc_ptr(CS%hf_dv_dt_visc,isd,ied,JsdB,JedB,nz)
   !  call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
   !  call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   !endif
@@ -1953,48 +1882,48 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
   endif
 
   CS%id_h_du_dt_visc = register_diag_field('ocean_model', 'h_du_dt_visc', diag%axesCuL, Time, &
-      'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', &
+      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_h_du_dt_visc > 0) then
     call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
   endif
 
   CS%id_h_dv_dt_visc = register_diag_field('ocean_model', 'h_dv_dt_visc', diag%axesCvL, Time, &
-      'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', &
+      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_h_dv_dt_visc > 0) then
     call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz)
   endif
 
   CS%id_h_du_dt_str = register_diag_field('ocean_model', 'h_du_dt_str', diag%axesCuL, Time, &
-      'Thickness Multiplied Zonal Acceleration from Surface Wind Stresses', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      'Thickness Multiplied Zonal Acceleration from Surface Wind Stresses', &
+      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_h_du_dt_str > 0) then
     call safe_alloc_ptr(ADp%du_dt_str,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
   endif
 
   CS%id_h_dv_dt_str = register_diag_field('ocean_model', 'h_dv_dt_str', diag%axesCvL, Time, &
-      'Thickness Multiplied Meridional Acceleration from Surface Wind Stresses', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      'Thickness Multiplied Meridional Acceleration from Surface Wind Stresses', &
+      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_h_dv_dt_str > 0) then
     call safe_alloc_ptr(ADp%dv_dt_str,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz)
   endif
 
   CS%id_du_dt_str_visc_rem = register_diag_field('ocean_model', 'du_dt_str_visc_rem', diag%axesCuL, Time, &
-      'Zonal Acceleration from Surface Wind Stresses multiplied by viscous remnant', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Zonal Acceleration from Surface Wind Stresses multiplied by viscous remnant', &
+      'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_du_dt_str_visc_rem > 0) then
     call safe_alloc_ptr(ADp%du_dt_str,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
   endif
 
   CS%id_dv_dt_str_visc_rem = register_diag_field('ocean_model', 'dv_dt_str_visc_rem', diag%axesCvL, Time, &
-      'Meridional Acceleration from Surface Wind Stresses multiplied by viscous remnant', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Meridional Acceleration from Surface Wind Stresses multiplied by viscous remnant', &
+      'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_dv_dt_str_visc_rem > 0) then
     call safe_alloc_ptr(ADp%dv_dt_str,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)