diff --git a/src/diagnostics/MOM_diagnose_KdWork.F90 b/src/diagnostics/MOM_diagnose_KdWork.F90
new file mode 100644
index 0000000000..12f8191619
--- /dev/null
+++ b/src/diagnostics/MOM_diagnose_KdWork.F90
@@ -0,0 +1,1175 @@
+!> Provides diagnostics of work due to a given diffusivity
+module MOM_diagnose_kdwork
+
+! This file is part of MOM6. See LICENSE.md for the license.
+
+use MOM_diag_mediator, only : diag_ctrl, time_type, post_data, register_diag_field
+use MOM_diag_mediator, only : register_scalar_field
+use MOM_error_handler, only : MOM_error, FATAL, WARNING
+use MOM_grid,          only : ocean_grid_type
+use MOM_unit_scaling,  only : unit_scale_type
+use MOM_variables,     only : thermo_var_ptrs
+use MOM_verticalGrid,  only : verticalGrid_type
+use MOM_spatial_means, only : global_area_integral
+
+implicit none ; private
+
+#include <MOM_memory.h>
+
+public vbf_CS
+public kdwork_diagnostics
+public Allocate_VBF_CS
+public Deallocate_VBF_CS
+public KdWork_init
+public KdWork_end
+
+!> This structure has memory for used in calculating diagnostics of diffusivity
+!! many of the diffusivity diagnostics are copies of other 3d arrays.  It could
+!! be written more efficiently, but it is less intrusive to copy into this structure
+!! and do all calculations in this module.  These diagnostics may be expensive for
+!! routine use.
+type vbf_CS
+  ! 3d varying Kd contributions
+  real, pointer, dimension(:,:,:) :: &
+    Bflx_salt => NULL(), & !< Salinity contribution to buoyancy flux at interfaces
+                           !! [H Z T-3 ~> m2 s-3 or kg m-1 s-3 = W m-3]
+    Bflx_temp => NULL(), & !< Temperature contribution to buoyancy flux at interfaces
+                           !! [H Z T-3 ~> m2 s-3 or kg m-1 s-3 = W m-3]
+    Bflx_salt_dz => NULL(), & !< Salinity contribution to integral of buoyancy flux over layer
+                              !! [H Z2 T-3 ~> m3 s-3 or kg m-1 s-3 = W m-2]
+    Bflx_temp_dz => NULL(), & !< Temperature contribution to integral of buoyancy flux over layer
+                              !! [H Z2 T-3 ~> m3 s-3 or kg m-1 s-3 = W m-2]
+    ! The following are all allocatable arrays that store copies of process driven Kd, so that
+    ! the process driven buoyancy flux and work can be derived at the end of the time step.
+    Kd_salt => NULL(), &   !< total diapycnal diffusivity of salt at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_temp => NULL(), &   !< total diapycnal diffusivity of heat at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_BBL => NULL(), &    !< diapycnal diffusivity due to BBL at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_ePBL => NULL(), &   !< diapycnal diffusivity due to ePBL at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_KS => NULL(), &     !< diapycnal diffusivity due to Kappa Shear at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_bkgnd => NULL(), &  !< diapycnal diffusivity due to Kd_bkgnd at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_ddiff_S => NULL(), &!< diapycnal diffusivity due to double diffusion of salt at interfaces
+                           !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_ddiff_T => NULL(), &!< diapycnal diffusivity due to double diffusion of heat at interfaces
+                           !![H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_leak => NULL(), &   !< diapycnal diffusivity due to Kd_leak at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_quad => NULL(), &   !< diapycnal diffusivity due to Kd_quad at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_itidal => NULL(), & !< diapycnal diffusivity due to Kd_itidal at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_Froude => NULL(), & !< diapycnal diffusivity due to Kd_Froude at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_slope => NULL(), &  !< diapycnal diffusivity due to Kd_slope at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_lowmode => NULL(), &!< diapycnal diffusivity due to Kd_lowmode at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_Niku => NULL(), &   !< diapycnal diffusivity due to Kd_Niku at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+    Kd_itides => NULL()    !< diapycnal diffusivity due to Kd_itides at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+
+  ! Constant Kd contributions
+  real :: Kd_add !< spatially uniform additional diapycnal diffusivity at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+                 !! a diagnostic for this diffusivity is not yet included, but this makes it straightforward to add
+
+  !>@{ Diagnostic IDs
+  integer :: id_Bdif     = -1, id_Bdif_salt  = -1, id_Bdif_temp  = -1
+  integer :: id_Bdif_dz  = -1, id_Bdif_salt_dz  = -1, id_Bdif_temp_dz  = -1
+  integer :: id_Bdif_idz  = -1, id_Bdif_salt_idz  = -1, id_Bdif_temp_idz  = -1
+  integer :: id_Bdif_idV  = -1, id_Bdif_salt_idV  = -1, id_Bdif_temp_idV  = -1
+  integer :: id_Bdif_ePBL  = -1, id_Bdif_dz_ePBL  = -1, id_Bdif_idz_ePBL  = -1, id_Bdif_idV_ePBL  = -1
+  integer :: id_Bdif_BBL  = -1, id_Bdif_dz_BBL  = -1, id_Bdif_idz_BBL  = -1, id_Bdif_idV_BBL  = -1
+  integer :: id_Bdif_KS  = -1, id_Bdif_dz_KS  = -1, id_Bdif_idz_KS  = -1, id_Bdif_idV_KS  = -1
+  integer :: id_Bdif_bkgnd  = -1, id_Bdif_dz_bkgnd  = -1, id_Bdif_idz_bkgnd  = -1, id_Bdif_idV_bkgnd  = -1
+  integer :: id_Bdif_ddiff_temp  = -1, id_Bdif_ddiff_salt  = -1
+  integer :: id_Bdif_dz_ddiff_temp  = -1, id_Bdif_dz_ddiff_salt  = -1
+  integer :: id_Bdif_idz_ddiff_temp  = -1, id_Bdif_idz_ddiff_salt  = -1
+  integer :: id_Bdif_idV_ddiff_temp  = -1, id_Bdif_idV_ddiff_salt  = -1
+  integer :: id_Bdif_leak  = -1, id_Bdif_dz_leak  = -1, id_Bdif_idz_leak  = -1, id_Bdif_idV_leak  = -1
+  integer :: id_Bdif_quad  = -1, id_Bdif_dz_quad  = -1, id_Bdif_idz_quad  = -1, id_Bdif_idV_quad  = -1
+  integer :: id_Bdif_itidal  = -1, id_Bdif_dz_itidal  = -1, id_Bdif_idz_itidal  = -1, id_Bdif_idV_itidal  = -1
+  integer :: id_Bdif_Froude  = -1, id_Bdif_dz_Froude  = -1, id_Bdif_idz_Froude  = -1, id_Bdif_idV_Froude  = -1
+  integer :: id_Bdif_slope  = -1, id_Bdif_dz_slope  = -1, id_Bdif_idz_slope  = -1, id_Bdif_idV_slope  = -1
+  integer :: id_Bdif_lowmode  = -1, id_Bdif_dz_lowmode  = -1, id_Bdif_idz_lowmode  = -1, id_Bdif_idV_lowmode  = -1
+  integer :: id_Bdif_Niku  = -1, id_Bdif_dz_Niku  = -1, id_Bdif_idz_Niku  = -1, id_Bdif_idV_Niku  = -1
+  integer :: id_Bdif_itides  = -1, id_Bdif_dz_itides  = -1, id_Bdif_idz_itides  = -1, id_Bdif_idV_itides  = -1
+  !>@}
+
+  logical :: do_bflx_salt = .false.  !< Logical flag to indicate if N2_salt should be computed
+  logical :: do_bflx_temp = .false.  !< Logical flag to indicate if N2_temp should be computed
+  logical :: do_bflx_salt_dz = .false.  !< Logical flag to indicate if N2_salt should be computed
+  logical :: do_bflx_temp_dz = .false.  !< Logical flag to indicate if N2_temp should be computed
+
+end type vbf_CS
+
+! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
+! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
+! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
+! vary with the Boussinesq approximation, the Boussinesq variant is given first.
+
+contains
+
+!> Loop over all implemented diffusivities to diagnose and output Kd Work/buoyancy fluxes
+subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz)
+  type(ocean_grid_type),      intent(in)    :: G       !< Grid type
+  type(verticalGrid_type),    intent(in)    :: GV      !< ocean vertical grid structure
+  type(unit_scale_type),      intent(in)    :: US      !< A dimensional unit scaling type
+  type(diag_ctrl), target,    intent(inout) :: diag    !< regulates diagnostic output
+  type (vbf_CS),              intent(inout) :: VBF     !< Vertical buoyancy flux structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
+                              intent(in)    :: N2_Salt !< Buoyancy frequency [T-2 ~> s-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
+                              intent(in)    :: N2_Temp !< Buoyancy frequency [T-2 ~> s-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                              intent(in)    :: dz      !< Grid spacing [Z ~> m]
+
+  ! Work arrays for computing buoyancy flux integrals
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: work3d_i
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work3d_l
+  real, dimension(SZI_(G),SZJ_(G)) :: work2d, work2d_salt, work2d_temp
+  real :: work, work_salt, work_temp
+
+  integer :: i, j, k, nz, isc, iec, jsc, jec
+
+  isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ;
+
+  nz = GV%ke
+
+  ! Compute total fluxes
+  if (VBF%id_Bdif_dz>0 .or. VBF%id_Bdif_salt_dz>0 .or. VBF%id_Bdif_temp_dz>0 .or. &
+      VBF%id_Bdif_idz>0 .or. VBF%id_Bdif_salt_idz>0 .or. VBF%id_Bdif_temp_idz>0 .or. &
+      VBF%id_Bdif_idV>0 .or. VBF%id_Bdif_salt_idV>0 .or. VBF%id_Bdif_temp_idV>0 ) then ! Doing vertical integrals
+    ! Do Salt
+    if (VBF%id_Bdif_salt_dz>0 .or. VBF%id_Bdif_dz>0 .or. VBF%id_Bdif_salt>0 .or. VBF%id_Bdif>0 .or. &
+        VBF%id_Bdif_idz>0 .or. VBF%id_Bdif_salt_idz>0 .or. VBF%id_Bdif_idV>0 .or. VBF%id_Bdif_salt_idV>0) &
+      call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_salt, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    ! Do Temp
+    if (VBF%id_Bdif_temp_dz>0 .or. VBF%id_Bdif_dz>0 .or. VBF%id_Bdif_temp>0 .or. VBF%id_Bdif>0 .or. &
+        VBF%id_Bdif_idz>0 .or. VBF%id_Bdif_temp_idz>0 .or. VBF%id_Bdif_idV>0 .or. VBF%id_Bdif_temp_idV>0) &
+       call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_temp, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_temp_idz>0 .or. VBF%id_Bdif_idz>0) then
+      work2d_temp(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d_temp(i,j) = work2d_temp(i,j) + VBF%Bflx_temp_dz(i,j,k)
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_temp_idV>0 .or. VBF%id_Bdif_idV>0) then
+      work_temp = 0.0
+      do k = 1,nz
+        work_temp = work_temp + global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, &
+                    tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+      enddo
+    endif
+    if (VBF%id_Bdif_salt_idz>0 .or. VBF%id_Bdif_idz>0) then
+      work2d_salt(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d_salt(i,j) = work2d_salt(i,j) + VBF%Bflx_salt_dz(i,j,k)
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_salt_idV>0 .or. VBF%id_Bdif_idV>0) then
+      work_salt = 0.0
+      do k = 1,nz
+        work_salt = work_salt + global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, &
+                    tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+      enddo
+    endif
+    work = work_temp + work_salt
+    do j = jsc,jec ; do i = isc,iec
+      work2d(i,j) = work2d_temp(i,j) + work2d_salt(i,j)
+    enddo ; enddo
+  elseif (VBF%id_Bdif>0 .or. VBF%id_Bdif_salt>0 .or. VBF%id_Bdif_temp>0) then ! Not doing vertical integrals
+    ! Do Salt
+    if (VBF%id_Bdif_salt>0 .or. VBF%id_Bdif>0) &
+      call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_salt, VBF%Bflx_salt)
+    if (VBF%id_Bdif_temp>0 .or. VBF%id_Bdif>0) &
+      call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_temp, VBF%Bflx_temp)
+  endif
+  ! Post total fluxes
+  if (VBF%id_Bdif_salt>0) call post_data(VBF%id_Bdif_salt, VBF%Bflx_salt, diag)
+  if (VBF%id_Bdif_temp>0) call post_data(VBF%id_Bdif_temp, VBF%Bflx_temp, diag)
+  if (VBF%id_Bdif>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_salt_dz>0) call post_data(VBF%id_Bdif_salt_dz, VBF%Bflx_salt_dz, diag)
+  if (VBF%id_Bdif_temp_dz>0) call post_data(VBF%id_Bdif_temp_dz, VBF%Bflx_temp_dz, diag)
+  if (VBF%id_Bdif_dz>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_salt_idz>0) call post_data(VBF%id_Bdif_salt_idz, work2d_salt, diag)
+  if (VBF%id_Bdif_temp_idz>0) call post_data(VBF%id_Bdif_temp_idz, work2d_temp, diag)
+  if (VBF%id_Bdif_idz>0) call post_data(VBF%id_Bdif_idz, work2d, diag)
+  if (VBF%id_Bdif_salt_idV>0) call post_data(VBF%id_Bdif_salt_idV, work_salt, diag)
+  if (VBF%id_Bdif_temp_idV>0) call post_data(VBF%id_Bdif_temp_idV, work_temp, diag)
+  if (VBF%id_Bdif_idV>0) call post_data(VBF%id_Bdif_idV, work, diag)
+
+  ! Compute ePBL fluxes
+  if (VBF%id_Bdif_dz_ePBL>0.or.VBF%id_Bdif_idz_ePBL>0.or.VBF%id_Bdif_idV_ePBL>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_ePBL, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_ePBL, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_ePBL>0) then
+      work2d(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k))
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_ePBL>0) then
+      work = 0.0
+      do k = 1,nz
+        work = work + &
+               (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + &
+                global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3))
+      enddo
+    endif
+  elseif (VBF%id_Bdif_ePBL>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_ePBL, VBF%Bflx_salt)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_ePBL, VBF%Bflx_temp)
+  endif
+  ! Post ePBL fluxes
+  if (VBF%id_Bdif_ePBL>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_ePBL, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_dz_ePBL>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz_ePBL, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_idz_ePBL>0) call post_data(VBF%id_Bdif_idz_ePBL, work2d, diag)
+  if (VBF%id_Bdif_idV_ePBL>0) call post_data(VBF%id_Bdif_idV_ePBL, work, diag)
+
+  ! Compute BBL fluxes
+  if (VBF%id_Bdif_dz_BBL>0.or.VBF%id_Bdif_idz_BBL>0.or.VBF%id_Bdif_idV_BBL>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_BBL, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_BBL, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_BBL>0) then
+      work2d(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k))
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_BBL>0) then
+      work = 0.0
+      do k = 1,nz
+        work = work + &
+               (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + &
+                global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3))
+      enddo
+    endif
+  elseif (VBF%id_Bdif_BBL>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_BBL, VBF%Bflx_salt)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_BBL, VBF%Bflx_temp)
+  endif
+  ! Post BBL fluxes
+  if (VBF%id_Bdif_BBL>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_BBL, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_dz_BBL>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz_BBL, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_idz_BBL>0) call post_data(VBF%id_Bdif_idz_BBL, work2d, diag)
+  if (VBF%id_Bdif_idV_BBL>0) call post_data(VBF%id_Bdif_idV_BBL, work, diag)
+
+  ! Compute Kappa Shear fluxes
+  if (VBF%id_Bdif_dz_KS>0.or.VBF%id_Bdif_idz_KS>0.or.VBF%id_Bdif_idV_KS>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_KS, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_KS, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_KS>0) then
+      work2d(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k))
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_KS>0) then
+      work = 0.0
+      do k = 1,nz
+        work = work + &
+               (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + &
+                global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3))
+      enddo
+    endif
+  elseif (VBF%id_Bdif_KS>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_KS, VBF%Bflx_salt)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_KS, VBF%Bflx_temp)
+  endif
+  ! Post Kappa Shear fluxes
+  if (VBF%id_Bdif_KS>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_KS, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_dz_KS>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz_KS, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_idz_KS>0) call post_data(VBF%id_Bdif_idz_KS, work2d, diag)
+  if (VBF%id_Bdif_idV_KS>0) call post_data(VBF%id_Bdif_idV_KS, work, diag)
+
+  ! Compute bkgnd fluxes
+  if (VBF%id_Bdif_dz_bkgnd>0.or.VBF%id_Bdif_idz_bkgnd>0.or.VBF%id_Bdif_idV_bkgnd>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_bkgnd, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_bkgnd, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_bkgnd>0) then
+      work2d(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k))
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_bkgnd>0) then
+      work = 0.0
+      do k = 1,nz
+        work = work + &
+               (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + &
+                global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3))
+      enddo
+    endif
+  elseif (VBF%id_Bdif_ePBL>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_bkgnd, VBF%Bflx_salt)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_bkgnd, VBF%Bflx_temp)
+  endif
+  ! Post bkgnd fluxes
+  if (VBF%id_Bdif_bkgnd>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_bkgnd, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_dz_bkgnd>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz_bkgnd, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_idz_bkgnd>0) call post_data(VBF%id_Bdif_idz_bkgnd, work2d, diag)
+  if (VBF%id_Bdif_idV_bkgnd>0) call post_data(VBF%id_Bdif_idV_bkgnd, work, diag)
+
+  ! Compute double diffusion fluxes
+  if (VBF%id_Bdif_dz_ddiff_temp>0.or.VBF%id_Bdif_idz_ddiff_temp>0.or.VBF%id_Bdif_idV_ddiff_temp>0) then
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_ddiff_T, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_ddiff_temp>0) then
+      work2d_temp(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d_temp(i,j) = work2d_temp(i,j) + VBF%Bflx_temp_dz(i,j,k)
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_ddiff_temp>0) then
+      work_temp = 0.0
+      do k = 1,nz
+        work_temp = work_temp + global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, &
+                    tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+      enddo
+    endif
+  elseif (VBF%id_Bdif_ddiff_temp>0) then
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_ddiff_T, VBF%Bflx_temp)
+  endif
+  if (VBF%id_Bdif_dz_ddiff_salt>0.or.VBF%id_Bdif_idz_ddiff_salt>0.or.VBF%id_Bdif_idV_ddiff_salt>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_ddiff_S, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    if (VBF%id_Bdif_idz_ddiff_salt>0) then
+      work2d_salt(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d_salt(i,j) = work2d_salt(i,j) + VBF%Bflx_salt_dz(i,j,k)
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_ddiff_salt>0) then
+      work_salt = 0.0
+      do k = 1,nz
+        work_salt = work_salt + global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, &
+                    tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+      enddo
+    endif
+  elseif (VBF%id_Bdif_ddiff_salt>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_ddiff_S, VBF%Bflx_salt)
+  endif
+  ! Post double diffusion fluxes
+  if (VBF%id_Bdif_ddiff_temp>0) call post_data(VBF%id_Bdif_ddiff_temp, VBF%Bflx_temp, diag)
+  if (VBF%id_Bdif_dz_ddiff_temp>0) call post_data(VBF%id_Bdif_dz_ddiff_temp, VBF%Bflx_temp_dz, diag)
+  if (VBF%id_Bdif_idz_ddiff_temp>0) call post_data(VBF%id_Bdif_idz_ddiff_temp, work2d_temp, diag)
+  if (VBF%id_Bdif_idV_ddiff_temp>0) call post_data(VBF%id_Bdif_idV_ddiff_temp, work_temp, diag)
+  if (VBF%id_Bdif_ddiff_salt>0) call post_data(VBF%id_Bdif_ddiff_salt, VBF%Bflx_salt, diag)
+  if (VBF%id_Bdif_dz_ddiff_salt>0) call post_data(VBF%id_Bdif_dz_ddiff_salt, VBF%Bflx_salt_dz, diag)
+  if (VBF%id_Bdif_idz_ddiff_salt>0) call post_data(VBF%id_Bdif_idz_ddiff_salt, work2d_salt, diag)
+  if (VBF%id_Bdif_idV_ddiff_salt>0) call post_data(VBF%id_Bdif_idV_ddiff_salt, work_salt, diag)
+
+  ! Compute Kd_leak fluxes
+  if (VBF%id_Bdif_dz_leak>0.or.VBF%id_Bdif_idz_leak>0.or.VBF%id_Bdif_idV_leak>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_leak, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_leak, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_leak>0) then
+      work2d(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k))
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_leak>0) then
+      work = 0.0
+      do k = 1,nz
+        work = work + &
+               (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + &
+                global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3))
+      enddo
+    endif
+  elseif (VBF%id_Bdif_leak>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_leak, VBF%Bflx_salt)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_leak, VBF%Bflx_temp)
+  endif
+  ! Post Kd_leak fluxes
+  if (VBF%id_Bdif_leak>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_leak, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_dz_leak>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz_leak, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_idz_leak>0) call post_data(VBF%id_Bdif_idz_leak, work2d, diag)
+  if (VBF%id_Bdif_idV_leak>0) call post_data(VBF%id_Bdif_idV_leak, work, diag)
+
+  ! Compute Kd_quad fluxes
+  if (VBF%id_Bdif_dz_quad>0.or.VBF%id_Bdif_idz_quad>0.or.VBF%id_Bdif_idV_quad>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_quad, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_quad, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_quad>0) then
+      work2d(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k))
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_quad>0) then
+      work = 0.0
+      do k = 1,nz
+        work = work + &
+               (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + &
+                global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3))
+      enddo
+    endif
+  elseif (VBF%id_Bdif_quad>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_quad, VBF%Bflx_salt)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_quad, VBF%Bflx_temp)
+  endif
+  ! Post Kd_quad fluxes
+  if (VBF%id_Bdif_quad>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_quad, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_dz_quad>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz_quad, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_idz_quad>0) call post_data(VBF%id_Bdif_idz_quad, work2d, diag)
+  if (VBF%id_Bdif_idV_quad>0) call post_data(VBF%id_Bdif_idV_quad, work, diag)
+
+  ! Compute Kd_itidal fluxes
+  if (VBF%id_Bdif_dz_itidal>0.or.VBF%id_Bdif_idz_itidal>0.or.VBF%id_Bdif_idV_itidal>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_itidal, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_itidal, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_itidal>0) then
+      work2d(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k))
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_itidal>0) then
+      work = 0.0
+      do k = 1,nz
+        work = work + &
+               (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + &
+                global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3))
+      enddo
+    endif
+  elseif (VBF%id_Bdif_itidal>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_itidal, VBF%Bflx_salt)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_itidal, VBF%Bflx_temp)
+  endif
+  ! Post Kd_itidal fluxes
+  if (VBF%id_Bdif_itidal>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_itidal, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_dz_itidal>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k)+VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz_itidal, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_idz_itidal>0) call post_data(VBF%id_Bdif_idz_itidal, work2d, diag)
+  if (VBF%id_Bdif_idV_itidal>0) call post_data(VBF%id_Bdif_idV_itidal, work, diag)
+
+  ! Compute Kd_Froude fluxes
+  if (VBF%id_Bdif_dz_Froude>0.or.VBF%id_Bdif_idz_Froude>0.or.VBF%id_Bdif_idV_Froude>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_Froude, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_Froude, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_Froude>0) then
+      work2d(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k))
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_Froude>0) then
+      work = 0.0
+      do k = 1,nz
+        work = work + &
+               (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + &
+                global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3))
+      enddo
+    endif
+  elseif (VBF%id_Bdif_Froude>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_Froude, VBF%Bflx_salt)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_Froude, VBF%Bflx_temp)
+  endif
+  ! Post Kd_Froude fluxes
+  if (VBF%id_Bdif_Froude>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_Froude, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_dz_Froude>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz_Froude, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_idz_Froude>0) call post_data(VBF%id_Bdif_idz_Froude, work2d, diag)
+  if (VBF%id_Bdif_idV_Froude>0) call post_data(VBF%id_Bdif_idV_Froude, work, diag)
+
+  ! Compute Kd_slope fluxes
+  if (VBF%id_Bdif_dz_slope>0.or.VBF%id_Bdif_idz_slope>0.or.VBF%id_Bdif_idV_slope>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_slope, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_slope, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_slope>0) then
+      work2d(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k))
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_slope>0) then
+      work = 0.0
+      do k = 1,nz
+        work = work + &
+               (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + &
+                global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3))
+      enddo
+    endif
+  elseif (VBF%id_Bdif_slope>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_slope, VBF%Bflx_salt)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_slope, VBF%Bflx_temp)
+  endif
+  ! Post Kd_slope fluxes
+  if (VBF%id_Bdif_slope>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_slope, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_dz_slope>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz_slope, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_idz_slope>0) call post_data(VBF%id_Bdif_idz_slope, work2d, diag)
+  if (VBF%id_Bdif_idV_slope>0) call post_data(VBF%id_Bdif_idV_slope, work, diag)
+
+  ! Compute Kd_lowmode fluxes
+  if (VBF%id_Bdif_dz_lowmode>0.or.VBF%id_Bdif_idz_lowmode>0.or.VBF%id_Bdif_idV_lowmode>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_lowmode, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_lowmode, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_lowmode>0) then
+      work2d(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k))
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_lowmode>0) then
+      work = 0.0
+      do k = 1,nz
+        work = work + &
+               (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + &
+                global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3))
+      enddo
+    endif
+  elseif (VBF%id_Bdif_lowmode>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_lowmode, VBF%Bflx_salt)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_lowmode, VBF%Bflx_temp)
+  endif
+  ! Post Kd_lowmode fluxes
+  if (VBF%id_Bdif_lowmode>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_lowmode, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_dz_lowmode>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz_lowmode, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_idz_lowmode>0) call post_data(VBF%id_Bdif_idz_lowmode, work2d, diag)
+  if (VBF%id_Bdif_idV_lowmode>0) call post_data(VBF%id_Bdif_idV_lowmode, work, diag)
+
+  ! Compute Kd_Niku fluxes
+  if (VBF%id_Bdif_dz_Niku>0 .or. VBF%id_Bdif_idz_Niku>0 .or. VBF%id_Bdif_idV_Niku>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_Niku, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_Niku, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_Niku>0) then
+      work2d(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k))
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_Niku>0) then
+      work = 0.0
+      do k = 1,nz
+        work = work + &
+               (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + &
+                global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3))
+      enddo
+    endif
+  elseif (VBF%id_Bdif_Niku>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_Niku, VBF%Bflx_salt)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_Niku, VBF%Bflx_temp)
+  endif
+  ! Post Kd_Niku fluxes
+  if (VBF%id_Bdif_Niku>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_lowmode, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_dz_Niku>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz_Niku, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_idz_Niku>0) call post_data(VBF%id_Bdif_idz_Niku, work2d, diag)
+  if (VBF%id_Bdif_idV_Niku>0) call post_data(VBF%id_Bdif_idV_Niku, work, diag)
+
+  ! Compute Kd_itides fluxes
+  if (VBF%id_Bdif_dz_itides>0 .or. VBF%id_Bdif_idz_itides>0 .or. VBF%id_Bdif_idV_itides>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_itides, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_itides, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz)
+    if (VBF%id_Bdif_idz_itides>0) then
+      work2d(:,:) = 0.0
+      do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+        work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k))
+      enddo ; enddo ; enddo
+    endif
+    if (VBF%id_Bdif_idV_itides>0) then
+      work = 0.0
+      do k = 1,nz
+        work = work + &
+               (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + &
+                global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3))
+      enddo
+    endif
+  elseif (VBF%id_Bdif_itides>0) then
+    call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_itides, VBF%Bflx_salt)
+    call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_itides, VBF%Bflx_temp)
+  endif
+  ! Post Kd_itides fluxes
+  if (VBF%id_Bdif_itides>0) then
+    work3d_i(:,:,:) = 0.0
+    do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
+      work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_itides, work3d_i, diag)
+  endif
+  if (VBF%id_Bdif_dz_itides>0) then
+    work3d_l(:,:,:) = 0.0
+    do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
+      work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k)
+    enddo ; enddo ; enddo
+    call post_data(VBF%id_Bdif_dz_itides, work3d_l, diag)
+  endif
+  if (VBF%id_Bdif_idz_itides>0) call post_data(VBF%id_Bdif_idz_itides, work2d, diag)
+  if (VBF%id_Bdif_idV_itides>0) call post_data(VBF%id_Bdif_idV_itides, work, diag)
+
+end subroutine KdWork_Diagnostics
+
+!> Diagnose the implied "work", or buoyancy forcing & its integral, due to a given diffusivity and column state.
+subroutine diagnoseKdWork(G, GV, N2, Kd, Bdif_flx, dz, Bdif_flx_dz)
+  type(ocean_grid_type),   intent(in)  :: G    !< Grid type
+  type(verticalGrid_type), intent(in)  :: GV   !< ocean vertical grid structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
+                           intent(in)  :: N2   !< Buoyancy frequency [T-2 ~> s-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
+                           intent(in)  :: Kd   !< Diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
+                           intent(out) :: Bdif_flx !< Buoyancy flux [H Z T-3 ~> m2 s-3 or kg m-1 s-3 = W m-3]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                 intent(in), optional :: dz    !< Grid spacing [Z ~> m]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                intent(out), optional :: Bdif_flx_dz !< Buoyancy flux over layer [H Z2 T-3 ~> m3 s-3 or kg s-3 = W m-2]
+
+  integer :: i, j, k
+
+  Bdif_flx(:,:,1) = 0.0
+  Bdif_flx(:,:,GV%ke+1) = 0.0
+  !$OMP parallel do default(shared)
+  do K=2,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
+    Bdif_flx(i,j,K) = - N2(i,j,K) * Kd(i,j,K)
+  enddo ; enddo; enddo
+
+  if (present(Bdif_flx_dz) .and. present(dz)) then
+    !$OMP parallel do default(shared)
+    do K=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
+      Bdif_flx_dz(i,j,k) = 0.5*(Bdif_flx(i,j,K)+Bdif_flx(i,j,K+1))*dz(i,j,k)
+    enddo ; enddo; enddo
+  endif
+
+end subroutine diagnoseKdWork
+
+!> Allocates arrays only when needed
+subroutine Allocate_VBF_CS(G, GV, VBF)
+  type(ocean_grid_type),   intent(in) :: G   !< ocean grid structure
+  type(verticalGrid_type), intent(in) :: GV  !< ocean vertical grid structure
+  type (vbf_CS),        intent(inout) :: VBF !< Vertical buoyancy flux structure
+
+  integer :: isd, ied, jsd, jed, nz
+
+  isd  = G%isd  ; ied = G%ied  ; jsd = G%jsd  ; jed = G%jed ; nz = GV%ke
+
+  if (VBF%do_bflx_salt) &
+    allocate(VBF%Bflx_salt(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%do_bflx_salt_dz) &
+     allocate(VBF%Bflx_salt_dz(isd:ied,jsd:jed,nz), source=0.0)
+  if (VBF%do_bflx_temp) &
+     allocate(VBF%Bflx_temp(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%do_bflx_temp_dz) &
+     allocate(VBF%Bflx_temp_dz(isd:ied,jsd:jed,nz), source=0.0)
+
+  if (VBF%do_bflx_salt .or. VBF%do_bflx_salt_dz ) &
+    allocate(VBF%Kd_salt(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%do_bflx_temp .or. VBF%do_bflx_temp_dz ) &
+    allocate(VBF%Kd_temp(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_BBL>0 .or. VBF%id_Bdif_dz_BBL>0 .or. VBF%id_Bdif_idV_BBL>0) &
+    allocate(VBF%Kd_BBL(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_ePBL>0 .or. VBF%id_Bdif_dz_ePBL>0 .or. VBF%id_Bdif_idV_ePBL>0) &
+    allocate(VBF%Kd_ePBL(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_KS>0 .or. VBF%id_Bdif_dz_KS>0 .or. VBF%id_Bdif_idV_KS>0) &
+    allocate(VBF%Kd_KS(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_bkgnd>0 .or. VBF%id_Bdif_dz_bkgnd>0 .or. VBF%id_Bdif_idV_bkgnd>0) &
+    allocate(VBF%Kd_bkgnd(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_ddiff_temp>0 .or. VBF%id_Bdif_dz_ddiff_temp>0 .or. VBF%id_Bdif_idV_ddiff_temp>0) &
+    allocate(VBF%Kd_ddiff_T(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_ddiff_salt>0 .or. VBF%id_Bdif_dz_ddiff_salt>0 .or. VBF%id_Bdif_idV_ddiff_salt>0) &
+    allocate(VBF%Kd_ddiff_S(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_leak>0 .or. VBF%id_Bdif_dz_leak>0 .or. VBF%id_Bdif_idV_leak>0) &
+    allocate(VBF%Kd_leak(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_quad>0 .or. VBF%id_Bdif_dz_quad>0 .or. VBF%id_Bdif_idV_quad>0) &
+    allocate(VBF%Kd_quad(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_itidal>0 .or. VBF%id_Bdif_dz_itidal>0 .or. VBF%id_Bdif_idV_itidal>0) &
+    allocate(VBF%Kd_itidal(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_Froude>0 .or. VBF%id_Bdif_dz_Froude>0 .or. VBF%id_Bdif_idV_Froude>0) &
+    allocate(VBF%Kd_Froude(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_slope>0 .or. VBF%id_Bdif_dz_slope>0 .or. VBF%id_Bdif_idV_slope>0) &
+    allocate(VBF%Kd_slope(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_lowmode>0 .or. VBF%id_Bdif_dz_lowmode>0 .or. VBF%id_Bdif_idV_lowmode>0) &
+    allocate(VBF%Kd_lowmode(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_Niku>0 .or. VBF%id_Bdif_dz_Niku>0 .or. VBF%id_Bdif_idV_Niku>0) &
+    allocate(VBF%Kd_Niku(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (VBF%id_Bdif_itides>0 .or. VBF%id_Bdif_dz_itides>0 .or. VBF%id_Bdif_idV_itides>0) &
+    allocate(VBF%Kd_itides(isd:ied,jsd:jed,nz+1), source=0.0)
+
+end subroutine Allocate_VBF_CS
+
+!> Deallocate any arrays that were allocated
+subroutine Deallocate_VBF_CS(VBF)
+  type (vbf_CS), intent(inout) :: VBF !< Vertical buoyancy flux structure
+
+  if (associated(VBF%Bflx_salt)) &
+    deallocate(VBF%Bflx_salt)
+  if (associated(VBF%Bflx_temp)) &
+    deallocate(VBF%Bflx_temp)
+  if (associated(VBF%Bflx_salt_dz)) &
+    deallocate(VBF%Bflx_salt_dz)
+  if (associated(VBF%Bflx_temp_dz)) &
+    deallocate(VBF%Bflx_temp_dz)
+  if (associated(VBF%Kd_salt)) &
+    deallocate(VBF%Kd_salt)
+  if (associated(VBF%Kd_temp)) &
+    deallocate(VBF%Kd_temp)
+  if (associated(VBF%Kd_BBL)) &
+    deallocate(VBF%Kd_BBL)
+  if (associated(VBF%Kd_ePBL)) &
+    deallocate(VBF%Kd_ePBL)
+  if (associated(VBF%Kd_KS)) &
+    deallocate(VBF%Kd_KS)
+  if (associated(VBF%Kd_bkgnd)) &
+    deallocate(VBF%Kd_bkgnd)
+  if (associated(VBF%Kd_ddiff_T)) &
+    deallocate(VBF%Kd_ddiff_T)
+  if (associated(VBF%Kd_ddiff_S)) &
+    deallocate(VBF%Kd_ddiff_S)
+  if (associated(VBF%Kd_leak)) &
+    deallocate(VBF%Kd_leak)
+  if (associated(VBF%Kd_quad)) &
+    deallocate(VBF%Kd_quad)
+  if (associated(VBF%Kd_itidal)) &
+    deallocate(VBF%Kd_itidal)
+  if (associated(VBF%Kd_Froude)) &
+    deallocate(VBF%Kd_Froude)
+  if (associated(VBF%Kd_slope)) &
+    deallocate(VBF%Kd_slope)
+  if (associated(VBF%Kd_lowmode)) &
+    deallocate(VBF%Kd_lowmode)
+  if (associated(VBF%Kd_Niku)) &
+    deallocate(VBF%Kd_Niku)
+  if (associated(VBF%Kd_itides)) &
+    deallocate(VBF%Kd_itides)
+
+end subroutine Deallocate_VBF_CS
+
+!> Handles all KdWork diagnostics and flags which calculations should be done.
+subroutine KdWork_init(Time, G,GV,US,diag,VBF,Use_KdWork_diag)
+  type(time_type), target                :: Time             !< model time
+  type(ocean_grid_type),   intent(in)    :: G        !< ocean grid structure
+  type(verticalGrid_type), intent(in)    :: GV       !< ocean vertical grid structure
+  type(unit_scale_type),   intent(in)    :: US       !< A dimensional unit scaling type
+  type(diag_ctrl), target, intent(inout) :: diag     !< regulates diagnostic output
+  type (vbf_CS), pointer,  intent(inout) :: VBF      !< Vertical buoyancy flux structure
+  logical,                 intent(out)   :: Use_KdWork_diag !< Flag if any output was turned on
+
+  allocate(VBF)
+
+  VBF%do_bflx_salt = .false.
+  VBF%do_bflx_salt_dz = .false.
+  VBF%do_bflx_temp = .false.
+  VBF%do_bflx_temp_dz = .false.
+
+  VBF%id_Bdif = register_diag_field('ocean_model',"Bflx_dia_diff", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz = register_diag_field('ocean_model',"Bflx_dia_diff_dz", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz = register_diag_field('ocean_model',"Bflx_dia_diff_idz", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV = register_scalar_field('ocean_model',"Bflx_dia_diff_idV", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_salt = register_diag_field('ocean_model',"Bflx_salt_dia_diff", diag%axesTi, &
+        Time, "Salinity contribution to diffusive diapycnal buoyancy flux across interfaces", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_salt_dz = register_diag_field('ocean_model',"Bflx_salt_dia_diff_dz", diag%axesTl, &
+        Time, "Salinity contribution to layer integral of diffusive diapycnal buoyancy flux.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_salt_idz = register_diag_field('ocean_model',"Bflx_salt_dia_diff_idz", diag%axesT1, &
+        Time, "Salinity contribution to layer integrated diffusive diapycnal buoyancy flux.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_salt_idV = register_scalar_field('ocean_model',"Bflx_salt_dia_diff_idV", Time, diag, &
+        "Salinity contribution to global integrated diffusive diapycnal buoyancy flux.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_temp = register_diag_field('ocean_model',"Bflx_temp_dia_diff", diag%axesTi, &
+        Time, "Temperature contribution to diffusive diapycnal buoyancy flux across interfaces", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_temp_dz = register_diag_field('ocean_model',"Bflx_temp_dia_diff_dz", diag%axesTl, &
+        Time, "Temperature contribution to layer integral of diffusive diapycnal buoyancy flux.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_temp_idz = register_diag_field('ocean_model',"Bflx_temp_dia_diff_idz", diag%axesT1, &
+        Time, "Temperature contribution to layer integrated diffusive diapycnal buoyancy flux.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_temp_idV = register_scalar_field('ocean_model',"Bflx_temp_dia_diff_idV", Time, diag, &
+        "Temperature contribution to global integrated diffusive diapycnal buoyancy flux.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_BBL = register_diag_field('ocean_model',"Bflx_dia_diff_BBL", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to the BBL parameterization.", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_BBL = register_diag_field('ocean_model',"Bflx_dia_diff_dz_BBL", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to the BBL parameterization.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_BBL = register_diag_field('ocean_model',"Bflx_dia_diff_idz_BBL", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to the BBL parameterization.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_BBL = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_BBL", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to BBL.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_ePBL = register_diag_field('ocean_model',"Bflx_dia_diff_ePBL", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to ePBL", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_ePBL = register_diag_field('ocean_model',"Bflx_dia_diff_dz_ePBL", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to ePBL.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_ePBL = register_diag_field('ocean_model',"Bflx_dia_diff_idz_ePBL", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to ePBL.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_ePBL = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_ePBL", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to ePBL.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_KS = register_diag_field('ocean_model',"Bflx_dia_diff_KS", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kappa Shear", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_KS = register_diag_field('ocean_model',"Bflx_dia_diff_dz_KS", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to Kappa Shear.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_KS = register_diag_field('ocean_model',"Bflx_dia_diff_idz_KS", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to Kappa Shear.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_KS = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_KS", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to Kappa Shear.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_bkgnd = register_diag_field('ocean_model',"Bflx_dia_diff_bkgnd", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to bkgnd mixing", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_bkgnd = register_diag_field('ocean_model',"Bflx_dia_diff_dz_bkgnd", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_bkgnd = register_diag_field('ocean_model',"Bflx_dia_diff_idz_bkgnd", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_bkgnd = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_bkgnd", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to Kd_bkgnd.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_ddiff_temp = register_diag_field('ocean_model',"Bflx_dia_diff_ddiff_heat", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to double diffusion of heat", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_ddiff_temp = register_diag_field('ocean_model',"Bflx_dia_diff_dz_ddiff_heat", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to double diffusion of heat.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_ddiff_temp = register_diag_field('ocean_model',"Bflx_dia_diff_idz_ddiff_heat", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to double diffusion of heat.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_ddiff_temp = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_ddiff_heat", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to double diffusion of heat.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_ddiff_salt = register_diag_field('ocean_model',"Bflx_dia_diff_ddiff_salt", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to double diffusion of salt", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_ddiff_salt = register_diag_field('ocean_model',"Bflx_dia_diff_dz_ddiff_salt", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to double diffusion of salt.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_ddiff_salt = register_diag_field('ocean_model',"Bflx_dia_diff_idz_ddiff_salt", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to double diffusion of salt.", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_ddiff_salt = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_ddiff_salt", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to double diffusion of salt.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_leak = register_diag_field('ocean_model',"Bflx_dia_diff_leak", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_leak mixing", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_leak = register_diag_field('ocean_model',"Bflx_dia_diff_dz_leak", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_leak = register_diag_field('ocean_model',"Bflx_dia_diff_idz_leak", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_leak = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_leak", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to Kd_leak.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_quad = register_diag_field('ocean_model',"Bflx_dia_diff_quad", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_quad mixing", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_quad = register_diag_field('ocean_model',"Bflx_dia_diff_dz_quad", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_quad = register_diag_field('ocean_model',"Bflx_dia_diff_idz_quad", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_quad = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_quad", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to Kd_quad.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_itidal = register_diag_field('ocean_model',"Bflx_dia_diff_itidal", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_itidal mixing", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_itidal = register_diag_field('ocean_model',"Bflx_dia_diff_dz_itidal", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_itidal = register_diag_field('ocean_model',"Bflx_dia_diff_idz_itidal", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_itidal = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_itidal", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to Kd_itidal.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_Froude = register_diag_field('ocean_model',"Bflx_dia_diff_Froude", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_Froude mixing", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_Froude = register_diag_field('ocean_model',"Bflx_dia_diff_dz_Froude", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_Froude = register_diag_field('ocean_model',"Bflx_dia_diff_idz_Froude", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_Froude = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_Froude", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to Kd_Froude.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_slope = register_diag_field('ocean_model',"Bflx_dia_diff_slope", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_slope mixing", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_slope = register_diag_field('ocean_model',"Bflx_dia_diff_dz_slope", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_slope = register_diag_field('ocean_model',"Bflx_dia_diff_idz_slope", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_slope = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_slope", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to Kd_slope.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_lowmode = register_diag_field('ocean_model',"Bflx_dia_diff_lowmode", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_lowmode mixing", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_lowmode = register_diag_field('ocean_model',"Bflx_dia_diff_dz_lowmode", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_lowmode = register_diag_field('ocean_model',"Bflx_dia_diff_idz_lowmode", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_lowmode = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_lowmode", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to Kd_lowmode.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_Niku = register_diag_field('ocean_model',"Bflx_dia_diff_Niku", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_Niku mixing", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_Niku = register_diag_field('ocean_model',"Bflx_dia_diff_dz_Niku", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_Niku = register_diag_field('ocean_model',"Bflx_dia_diff_idz_Niku", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_Niku = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_Niku", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to Kd_Niku.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  VBF%id_Bdif_itides = register_diag_field('ocean_model',"Bflx_dia_diff_itides", diag%axesTi, &
+        Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_itides mixing", &
+        "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3)
+  VBF%id_Bdif_dz_itides = register_diag_field('ocean_model',"Bflx_dia_diff_dz_itides", diag%axesTl, &
+        Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idz_itides = register_diag_field('ocean_model',"Bflx_dia_diff_idz_itides", diag%axesT1, &
+        Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
+        "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)
+  VBF%id_Bdif_idV_itides = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_itides", Time, diag, &
+        "Global integrated diffusive diapycnal buoyancy flux due to Kd_itides.", &
+        units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2)
+
+  if (VBF%id_Bdif_dz>0 .or. VBF%id_Bdif_salt_dz>0 .or. VBF%id_Bdif_dz_BBL>0 .or. &
+      VBF%id_Bdif_dz_ePBL>0 .or. VBF%id_Bdif_dz_KS>0 .or. VBF%id_Bdif_dz_bkgnd>0 .or. &
+      VBF%id_Bdif_dz_ddiff_salt>0 .or. VBF%id_Bdif_dz_leak>0 .or. VBF%id_Bdif_dz_quad>0 .or. &
+      VBF%id_Bdif_dz_itidal>0 .or. VBF%id_Bdif_dz_Froude>0 .or. VBF%id_Bdif_dz_slope>0 .or. &
+      VBF%id_Bdif_dz_lowmode>0 .or. VBF%id_Bdif_dz_Niku>0 .or. VBF%id_Bdif_dz_itides>0 .or. &
+      VBF%id_Bdif_idV>0 .or. VBF%id_Bdif_salt_idV>0 .or. VBF%id_Bdif_idV_BBL>0 .or. &
+      VBF%id_Bdif_idV_ePBL>0 .or. VBF%id_Bdif_idV_KS>0 .or. VBF%id_Bdif_idV_bkgnd>0 .or. &
+      VBF%id_Bdif_idV_ddiff_salt>0 .or. VBF%id_Bdif_idV_leak>0 .or. VBF%id_Bdif_idV_quad>0 .or. &
+      VBF%id_Bdif_idV_itidal>0 .or. VBF%id_Bdif_idV_Froude>0 .or. VBF%id_Bdif_idV_slope>0 .or. &
+      VBF%id_Bdif_idV_lowmode>0 .or. VBF%id_Bdif_idV_Niku>0 .or. VBF%id_Bdif_idV_itides>0 .or. &
+      VBF%id_Bdif_idz>0 .or. VBF%id_Bdif_salt_idz>0 .or. VBF%id_Bdif_idz_BBL>0 .or. &
+      VBF%id_Bdif_idz_ePBL>0 .or. VBF%id_Bdif_idz_KS>0 .or. VBF%id_Bdif_idz_bkgnd>0 .or. &
+      VBF%id_Bdif_idz_ddiff_salt>0 .or. VBF%id_Bdif_idz_leak>0 .or. VBF%id_Bdif_idz_quad>0 .or. &
+      VBF%id_Bdif_idz_itidal>0 .or. VBF%id_Bdif_idz_Froude>0 .or. VBF%id_Bdif_idz_slope>0 .or. &
+      VBF%id_Bdif_idz_lowmode>0 .or. VBF%id_Bdif_idz_Niku>0 .or. VBF%id_Bdif_idz_itides>0 ) then
+    VBF%do_bflx_salt_dz = .true.
+  endif
+  if (VBF%id_Bdif_dz>0 .or. VBF%id_Bdif_temp_dz>0 .or. VBF%id_Bdif_dz_BBL>0 .or. &
+      VBF%id_Bdif_dz_ePBL>0 .or. VBF%id_Bdif_dz_KS>0 .or. VBF%id_Bdif_dz_bkgnd>0 .or. &
+      VBF%id_Bdif_dz_ddiff_temp>0 .or. VBF%id_Bdif_dz_leak>0 .or. VBF%id_Bdif_dz_quad>0 .or. &
+      VBF%id_Bdif_dz_itidal>0 .or. VBF%id_Bdif_dz_Froude>0 .or. VBF%id_Bdif_dz_slope>0 .or. &
+      VBF%id_Bdif_dz_lowmode>0 .or. VBF%id_Bdif_dz_Niku>0 .or. VBF%id_Bdif_dz_itides>0 .or. &
+      VBF%id_Bdif_idV>0 .or. VBF%id_Bdif_temp_idV>0 .or. VBF%id_Bdif_idV_BBL>0 .or. &
+      VBF%id_Bdif_idV_ePBL>0 .or. VBF%id_Bdif_idV_KS>0 .or. VBF%id_Bdif_idV_bkgnd>0 .or. &
+      VBF%id_Bdif_idV_ddiff_temp>0 .or. VBF%id_Bdif_idV_leak>0 .or. VBF%id_Bdif_idV_quad>0 .or. &
+      VBF%id_Bdif_idV_itidal>0 .or. VBF%id_Bdif_idV_Froude>0 .or. VBF%id_Bdif_idV_slope>0 .or. &
+      VBF%id_Bdif_idV_lowmode>0 .or. VBF%id_Bdif_idV_Niku>0 .or. VBF%id_Bdif_idV_itides>0 .or. &
+      VBF%id_Bdif_idz>0 .or. VBF%id_Bdif_temp_idz>0 .or. VBF%id_Bdif_idz_BBL>0 .or. &
+      VBF%id_Bdif_idz_ePBL>0 .or. VBF%id_Bdif_idz_KS>0 .or. VBF%id_Bdif_idz_bkgnd>0 .or. &
+      VBF%id_Bdif_idz_ddiff_temp>0 .or. VBF%id_Bdif_idz_leak>0 .or. VBF%id_Bdif_idz_quad>0 .or. &
+      VBF%id_Bdif_idz_itidal>0 .or. VBF%id_Bdif_idz_Froude>0 .or. VBF%id_Bdif_idz_slope>0 .or. &
+      VBF%id_Bdif_idz_lowmode>0 .or. VBF%id_Bdif_idz_Niku>0 .or. VBF%id_Bdif_idz_itides>0 ) then
+    VBF%do_bflx_temp_dz = .true.
+  endif
+  if (VBF%id_Bdif>0 .or. VBF%id_Bdif_salt>0 .or. VBF%id_Bdif_BBL>0 .or. &
+      VBF%id_Bdif_ePBL>0 .or. VBF%id_Bdif_KS>0 .or. VBF%id_Bdif_bkgnd>0 .or. &
+      VBF%id_Bdif_ddiff_salt>0 .or. VBF%id_Bdif_leak>0 .or. VBF%id_Bdif_quad>0 .or. &
+      VBF%id_Bdif_itidal>0 .or. VBF%id_Bdif_Froude>0 .or. VBF%id_Bdif_slope>0 .or. &
+      VBF%id_Bdif_lowmode>0 .or. VBF%id_Bdif_Niku>0 .or. VBF%id_Bdif_itides>0 .or. &
+      VBF%do_bflx_salt_dz) then
+    VBF%do_bflx_salt = .true.
+  endif
+  if (VBF%id_Bdif>0 .or. VBF%id_Bdif_temp>0 .or. VBF%id_Bdif_BBL>0 .or. &
+      VBF%id_Bdif_ePBL>0 .or. VBF%id_Bdif_KS>0 .or. VBF%id_Bdif_bkgnd>0 .or. &
+      VBF%id_Bdif_ddiff_temp>0 .or. VBF%id_Bdif_leak>0 .or. VBF%id_Bdif_quad>0 .or. &
+      VBF%id_Bdif_itidal>0 .or. VBF%id_Bdif_Froude>0 .or. VBF%id_Bdif_slope>0 .or. &
+      VBF%id_Bdif_lowmode>0 .or. VBF%id_Bdif_Niku>0 .or. VBF%id_Bdif_itides>0 .or. &
+      VBF%do_bflx_temp_dz) then
+    VBF%do_bflx_temp = .true.
+  endif
+
+  Use_KdWork_diag = (VBF%do_bflx_salt .or. VBF%do_bflx_temp .or. VBF%do_bflx_salt_dz .or. VBF%do_bflx_temp_dz)
+
+end subroutine KdWork_init
+
+!> Deallocates control structrue
+subroutine KdWork_end(VBF)
+  type (vbf_CS), pointer,  intent(inout) :: VBF      !< Vertical buoyancy flux structure
+
+  if (associated(VBF)) deallocate(VBF)
+
+end subroutine KdWork_end
+
+!> \namespace mom_diagnose_kdwork
+!!
+!!    The subroutine diagnoseKdWork diagnoses the energetics associated with various vertical diffusivities
+!!    inside MOM6 diabatic routines.
+!!
+
+end module MOM_diagnose_kdwork
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index cb17b54e6d..7b2fcd5b5f 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -22,6 +22,8 @@ module MOM_diabatic_driver
 use MOM_diag_mediator,       only : diag_copy_diag_to_storage, diag_copy_storage_to_diag
 use MOM_diag_mediator,       only : diag_save_grids, diag_restore_grids
 use MOM_diagnose_mld,        only : diagnoseMLDbyDensityDifference, diagnoseMLDbyEnergy
+use MOM_diagnose_kdwork,     only : vbf_CS, KdWork_init, KdWork_end, KdWork_diagnostics
+use MOM_diagnose_kdwork,     only : Allocate_VBF_CS, Deallocate_VBF_CS
 use MOM_diapyc_energy_req,   only : diapyc_energy_req_init, diapyc_energy_req_end
 use MOM_diapyc_energy_req,   only : diapyc_energy_req_calc, diapyc_energy_req_test, diapyc_energy_req_CS
 use MOM_CVMix_conv,          only : CVMix_conv_init, CVMix_conv_cs
@@ -33,7 +35,8 @@ module MOM_diabatic_driver
 use MOM_energetic_PBL,       only : energetic_PBL_get_MLD
 use MOM_entrain_diffusive,   only : entrainment_diffusive, entrain_diffusive_init
 use MOM_entrain_diffusive,   only : entrain_diffusive_CS
-use MOM_EOS,                 only : calculate_density, calculate_TFreeze, EOS_domain
+use MOM_EOS,                 only : calculate_density, calculate_density_derivs, calculate_TFreeze
+use MOM_EOS,                 only : calculate_specific_vol_derivs, EOS_domain
 use MOM_error_handler,       only : MOM_error, FATAL, WARNING, callTree_showQuery,MOM_mesg
 use MOM_error_handler,       only : callTree_enter, callTree_leave, callTree_waypoint
 use MOM_file_parser,         only : get_param, log_version, param_file_type, read_param
@@ -177,12 +180,15 @@ module MOM_diabatic_driver
   real    :: MLD_En_vals(3)          !< Energy values for energy mixed layer diagnostics [R Z3 T-2 ~> J m-2]
   real    :: ref_h_mld = 0.0         !< The depth of the "surface"  density used in a difference mixed based
                                      !! MLD calculation [Z ~> m].
+  logical :: Use_KdWork_diag = .false.  !< Logical flag to indicate if any Kd_work diagnostics are on.
+  logical :: Use_N2_diag = .false.   !< Logical flag to indicate if any N2 diagnostics are on.
 
   !>@{ Diagnostic IDs
   integer :: id_ea       = -1, id_eb       = -1 ! used by layer diabatic
   integer :: id_ea_t     = -1, id_eb_t     = -1, id_ea_s   = -1, id_eb_s     = -1
   integer :: id_Kd_heat  = -1, id_Kd_salt  = -1, id_Kd_int = -1, id_Kd_ePBL  = -1
   integer :: id_Tdif     = -1, id_Sdif     = -1, id_Tadv   = -1, id_Sadv     = -1
+  integer :: id_N2_dd    = -1, id_N2_salt_dd = -1, id_N2_temp_dd
   ! These are handles to diagnostics related to the mixed layer properties.
   integer :: id_MLD_003 = -1, id_MLD_0125 = -1, id_MLD_user = -1, id_mlotstsq = -1
   integer :: id_MLD_003_zr = -1, id_MLD_003_rr = -1
@@ -236,6 +242,8 @@ module MOM_diabatic_driver
   type(diapyc_energy_req_CS),   pointer :: diapyc_en_rec_CSp     => NULL() !< Control structure for a child module
   type(oda_incupd_CS),          pointer :: oda_incupd_CSp        => NULL() !< Control structure for a child module
   type(int_tide_CS),            pointer :: int_tide_CSp          => NULL() !< Control structure for a child module
+  type(vbf_CS),                 pointer :: VBF                   => NULL() !< Control structure for a child module
+
 
   type(bulkmixedlayer_CS) :: bulkmixedlayer         !< Bulk mixed layer control structure
   type(CVMix_conv_CS) :: CVMix_conv                 !< CVMix convection control structure
@@ -564,7 +572,9 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
     KPP_NLTscalar, & ! KPP non-local transport for scalars [nondim]
     KPP_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3]
     Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
-    Sdif_flx    ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
+    Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
+    N2_salt, &  !< Salinity contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2]
+    N2_temp     !< Temperature contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2]
 
   real, dimension(SZI_(G),SZJ_(G)) :: &
     U_star, &    ! The friction velocity [Z T-1 ~> m s-1].
@@ -572,16 +582,29 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
     KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
     SkinBuoyFlux ! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
 
+  real, dimension(SZI_(G)) :: &
+    p_i ,&      ! Pressure at the interface [R L2 T-2 ~> Pa]
+    d_pres, &   ! pressure change across a layer [R L2 T-2 ~> Pa]
+    T_i, &      ! Temperature at the interface [C ~> degC]
+    S_i, &      ! Salinity at the interface [S ~> ppt]
+    drhodS, &   ! Local change in density w.r.t. salinity using model EOS & state [R C-1 ~> kg m-3 ppt-1]
+    drhodT, &   ! Local change in density w.r.t. temperature using model EOS & state [R C-1 ~> kg m-3 degC-1]
+    dSpV_dT, &  ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
+    dSpV_dS     ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
+
   logical, dimension(SZI_(G)) :: &
     in_boundary  ! True if there are no massive layers below, where massive is defined as
                  ! sufficiently thick that the no-flux boundary conditions have not restricted
                  ! the entrainment - usually sqrt(Kd*dt).
 
+  real :: h_neglect    ! A thickness that is so small it is usually lost
+                       ! in roundoff and can be neglected [H ~> m or kg m-2]
   real :: dz_neglect   ! A vertical distance that is so small it is usually lost
                        ! in roundoff and can be neglected [Z ~> m]
   real :: dz_neglect2  ! dz_neglect^2 [Z2 ~> m2]
   real :: add_ent      ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2]
   real :: I_dzval      ! The inverse of the thicknesses averaged to interfaces [Z-1 ~> m-1]
+  real :: I_h          ! The inverse of the thicknesses averaged to interfaces [H-1 ~> m-1 or m2 kg-1]
   real :: Tr_ea_BBL    ! The diffusive tracer thickness in the BBL that is
                        ! coupled to the bottom within a timestep [H ~> m or kg m-2]
   real :: Kd_add_here    ! An added diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
@@ -589,6 +612,12 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
 
   real :: Ent_int ! The diffusive entrainment rate at an interface [H ~> m or kg m-2]
   real :: Idt     ! The inverse time step [T-1 ~> s-1]
+  real :: g_Rho0       ! G_Earth/Rho0 [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2]
+  real :: H_to_pres    ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
+  real :: alt_H_to_pres! A conversion factor from thicknesses to pressure w/ alternative scaling [R Z T-2 ~> Pa m-1]
+  logical :: nonBous   ! True if not using the Boussinesq approximation
+
+  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
 
   logical :: showCallTree ! If true, show the call tree
   integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
@@ -596,6 +625,12 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
   is   = G%isc  ; ie  = G%iec  ; js  = G%jsc  ; je  = G%jec ; nz = GV%ke
   Isq  = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
   dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect
+  h_neglect = GV%H_subroundoff
+
+  nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq)
+  g_Rho0  = GV%g_Earth_Z_T2 / GV%H_to_RZ
+  H_to_pres = GV%H_to_RZ * GV%g_Earth
+  alt_H_to_pres = H_to_pres * US%L_to_Z**2 * GV%Z_to_H
 
   Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0
 
@@ -605,6 +640,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
   ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
   call enable_averages(dt, Time_end, CS%diag)
 
+  if (CS%Use_KdWork_diag) call Allocate_VBF_CS(G, GV, CS%VBF)
+
   if (CS%use_geothermal) then
     call cpu_clock_begin(id_clock_geothermal)
     call geothermal_in_place(h, tv, dt, G, GV, US, CS%geothermal, halo=CS%halo_TS_diff)
@@ -641,10 +678,10 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
     call MOM_state_chksum("before set_diffusivity", u, v, h, G, GV, US, haloshift=CS%halo_TS_diff)
   if (CS%double_diffuse) then
     call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_int, G, GV, US, &
-                         CS%set_diff_CSp, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S)
+                         CS%set_diff_CSp, CS%VBF, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S)
   else
     call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_int, G, GV, US, &
-                         CS%set_diff_CSp)
+                         CS%set_diff_CSp, CS%VBF)
   endif
   call cpu_clock_end(id_clock_set_diffusivity)
   if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)")
@@ -991,6 +1028,12 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
   ! target grids for vertical remapping may need to be regenerated.
   call diag_update_remap_grids(CS%diag)
 
+  ! Set diffusivities for VBF diagnostics if enabled
+  if (CS%use_energetic_PBL .and. associated(CS%VBF%Kd_ePBL)) CS%VBF%Kd_ePBL(:,:,:) = Kd_ePBL(:,:,:)
+  if (associated(CS%VBF%Kd_salt)) CS%VBF%Kd_temp(:,:,:) = Kd_heat(:,:,:)
+  if (associated(CS%VBF%Kd_temp)) CS%VBF%Kd_salt(:,:,:) = Kd_salt(:,:,:)
+
+
   ! Diagnose the diapycnal diffusivities and other related quantities.
   if (CS%id_Kd_int  > 0) call post_data(CS%id_Kd_int,  Kd_int,  CS%diag)
   if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag)
@@ -1025,6 +1068,68 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
     if (CS%id_Sdif > 0) call post_data(CS%id_Sdif, Sdif_flx, CS%diag)
   endif
 
+  if (CS%Use_KdWork_diag .or. CS%Use_N2_diag) then
+    N2_salt(:,:,:) = 0.0
+    N2_temp(:,:,:) = 0.0
+    !Compute N2 and don't mask negatives here
+    EOSdom(:) = EOS_domain(G%HI)
+    if (nonBous) then
+      !$OMP parallel do default(shared)
+      do j=js,je
+        if (associated(tv%p_surf)) then
+          do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo
+        else
+          do i=is,ie ; p_i(i) = 0.0 ; enddo
+        endif
+        do K=2,nz
+          do i=is,ie
+            p_i(i) = p_i(i) + H_to_pres * h(i,j,k-1)
+          enddo
+          T_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))
+          S_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))
+          call calculate_specific_vol_derivs(T_i, S_i, p_i, dSpV_dT, dSpV_dS, tv%eqn_of_state, EOSdom)
+          do i=is,ie
+            I_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k)))
+            N2_salt(i,j,K) = (tv%S(i,j,k-1) - tv%S(i,j,k)) * (dSpv_dS(i) * (alt_H_to_pres * I_dzval))
+            N2_temp(i,j,K) = (tv%T(i,j,k-1) - tv%T(i,j,k)) * (dSpV_dT(i) * (alt_H_to_pres * I_dzval))
+          enddo
+        enddo
+      enddo
+    else
+      !$OMP parallel do default(shared)
+      do j=js,je
+        if (associated(tv%p_surf)) then
+          do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo
+        else
+          do i=is,ie ; p_i(i) = 0.0 ; enddo
+        endif
+        do K=2,nz
+          do i=is,ie
+            p_i(i) = p_i(i) + H_to_pres* h(i,j,k-1)
+          enddo
+          T_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))
+          S_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))
+          call calculate_density_derivs(T_i, S_i, p_i, dRhodT, dRhodS, tv%eqn_of_state, EOSdom)
+          do i=is,ie
+            I_h = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))
+            N2_salt(i,j,K) = -(tv%S(i,j,k-1) - tv%S(i,j,k)) * (dRhodS(i) * (g_rho0 * I_h))
+            N2_temp(i,j,K) = -(tv%T(i,j,k-1) - tv%T(i,j,k)) * (dRhodT(i) * (g_rho0 * I_h))
+          enddo
+        enddo
+      enddo
+    endif
+    if (CS%id_N2_dd>0)      call post_data(CS%id_N2_dd, N2_salt(:,:,:)+N2_temp(:,:,:), CS%diag)
+    if (CS%id_N2_salt_dd>0) call post_data(CS%id_N2_salt_dd, N2_salt, CS%diag)
+    if (CS%id_N2_temp_dd>0) call post_data(CS%id_N2_temp_dd, N2_temp, CS%diag)
+
+    if (CS%Use_KdWork_diag) then
+       call KdWork_diagnostics(G,GV,US,CS%diag,CS%VBF,N2_salt,N2_temp,dz)
+    endif
+
+    call deallocate_VBF_CS(CS%VBF)
+
+  endif
+
   ! mixing of passive tracers from massless boundary layers to interior
   call cpu_clock_begin(id_clock_tracers)
 
@@ -1179,7 +1284,9 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
     KPP_NLTscalar, & ! KPP non-local transport for scalars [nondim]
     KPP_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3]
     Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
-    Sdif_flx    ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
+    Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
+    N2_salt, &  !< Salinity contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2]
+    N2_temp     !< Temperature contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2]
 
   real, dimension(SZI_(G),SZJ_(G)) :: &
     U_star, &    ! The friction velocity [Z T-1 ~> m s-1].
@@ -1191,23 +1298,50 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
     in_boundary  ! True if there are no massive layers below, where massive is defined as
                  ! sufficiently thick that the no-flux boundary conditions have not restricted
                  ! the entrainment - usually sqrt(Kd*dt).
+
+  real, dimension(SZI_(G)) :: &
+    p_i ,&      ! Pressure at the interface [R L2 T-2 ~> Pa]
+    d_pres, &   ! pressure change across a layer [R L2 T-2 ~> Pa]
+    T_i, &      ! Temperature at the interface [C ~> degC]
+    S_i, &      ! Salinity at the interface [S ~> ppt]
+    drhodS, &   ! Local change in density w.r.t. salinity using model EOS & state [R C-1 ~> kg m-3 ppt-1]
+    drhodT, &   ! Local change in density w.r.t. temperature using model EOS & state [R C-1 ~> kg m-3 degC-1]
+    dSpV_dT, &  ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
+    dSpV_dS     ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
+
+  real :: h_neglect    ! A thickness that is so small it is usually lost
+                       ! in roundoff and can be neglected [H ~> m or kg m-2]
   real :: dz_neglect   ! A vertical distance that is so small it is usually lost
                        ! in roundoff and can be neglected [Z ~> m]
   real :: dz_neglect2  ! dz_neglect^2 [Z2 ~> m2]
   real :: add_ent      ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2]
   real :: I_dzval      ! The inverse of the thicknesses averaged to interfaces [Z-1 ~> m-1]
+  real :: I_h          ! The inverse of the thicknesses averaged to interfaces [H-1 ~> m-1 or m2 kg-1]
   real :: Tr_ea_BBL    ! The diffusive tracer thickness in the BBL that is
                        ! coupled to the bottom within a timestep [H ~> m or kg m-2]
   real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2].
   real :: Kd_add_here    ! An added diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
-  real :: Idt     ! The inverse time step [T-1 ~> s-1]
+  real :: Idt          ! The inverse time step [T-1 ~> s-1]
+  real :: g_Rho0       ! G_Earth/Rho0 [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2]
+  real :: H_to_pres    ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
+  real :: alt_H_to_pres! A conversion factor from thicknesses to pressure w/ alternative scaling [R Z T-2 ~> Pa m-1]
+  logical :: nonBous   ! True if not using the Boussinesq approximation
+
+  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
 
   logical :: showCallTree ! If true, show the call tree
-  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
+  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, Isq, Ieq, Jsq, Jeq, nz
 
   is   = G%isc  ; ie  = G%iec  ; js  = G%jsc  ; je  = G%jec ; nz = GV%ke
+  isd  = G%isd  ; ied = G%ied  ; jsd = G%jsd  ; jed = G%jed
   Isq  = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
   dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect
+  h_neglect = GV%H_subroundoff
+
+  nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq)
+  g_Rho0  = GV%g_Earth_Z_T2 / GV%H_to_RZ
+  H_to_pres = GV%H_to_RZ * GV%g_Earth
+  alt_H_to_pres = H_to_pres * US%L_to_Z**2 * GV%Z_to_H
 
   Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0
   ent_s(:,:,:) = 0.0 ; ent_t(:,:,:) = 0.0
@@ -1221,6 +1355,8 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
   ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
   call enable_averages(dt, Time_end, CS%diag)
 
+  if (CS%Use_KdWork_diag) call Allocate_VBF_CS(G, GV, CS%VBF)
+
   if (CS%use_geothermal) then
     call cpu_clock_begin(id_clock_geothermal)
     call geothermal_in_place(h, tv, dt, G, GV, US, CS%geothermal, halo=CS%halo_TS_diff)
@@ -1257,10 +1393,10 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
     call MOM_state_chksum("before set_diffusivity", u, v, h, G, GV, US, haloshift=CS%halo_TS_diff)
   if (CS%double_diffuse) then
     call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_heat, G, GV, US, &
-                         CS%set_diff_CSp, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S)
+                         CS%set_diff_CSp, CS%VBF, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S)
   else
     call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_heat, G, GV, US, &
-                         CS%set_diff_CSp)
+                         CS%set_diff_CSp, CS%VBF)
   endif
   call cpu_clock_end(id_clock_set_diffusivity)
   if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)")
@@ -1535,12 +1671,18 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
   ! target grids for vertical remapping may need to be regenerated.
   call diag_update_remap_grids(CS%diag)
 
+  ! Set diffusivities for VBF diagnostics if enabled
+  if (CS%use_energetic_PBL .and. associated(CS%VBF%Kd_ePBL)) CS%VBF%Kd_ePBL(:,:,:) = Kd_ePBL(:,:,:)
+  if (associated(CS%VBF%Kd_salt)) CS%VBF%Kd_temp(:,:,:) = Kd_heat(:,:,:)
+  if (associated(CS%VBF%Kd_temp)) CS%VBF%Kd_salt(:,:,:) = Kd_salt(:,:,:)
+
   ! Diagnose the diapycnal diffusivities and other related quantities.
   if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag)
   if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag)
   if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, Kd_ePBL, CS%diag)
   if (CS%id_Kd_int > 0) then
     if (CS%double_diffuse .or. CS%useKPP) then
+      ! Using this as a work array might cause confusion.
       do K=1,nz ; do j=js,je ; do i=is,ie
         Kd_heat(i,j,k) = min(Kd_heat(i,j,k), Kd_salt(i,j,k))
       enddo ; enddo ; enddo
@@ -1575,6 +1717,68 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
     if (CS%id_Sdif > 0) call post_data(CS%id_Sdif, Sdif_flx, CS%diag)
   endif
 
+  if (CS%Use_KdWork_diag .or. CS%Use_N2_diag) then
+    N2_salt(:,:,:) = 0.0
+    N2_temp(:,:,:) = 0.0
+    !Compute N2 and don't mask negatives here
+    EOSdom(:) = EOS_domain(G%HI)
+    if (nonBous) then
+      !$OMP parallel do default(shared)
+      do j=js,je
+        if (associated(tv%p_surf)) then
+          do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo
+        else
+          do i=is,ie ; p_i(i) = 0.0 ; enddo
+        endif
+        do K=2,nz
+          do i=is,ie
+            p_i(i) = p_i(i) + H_to_pres * h(i,j,k-1)
+          enddo
+          T_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))
+          S_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))
+          call calculate_specific_vol_derivs(T_i, S_i, p_i, dSpV_dT, dSpV_dS, tv%eqn_of_state, EOSdom)
+          do i=is,ie
+            I_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k)))
+            N2_salt(i,j,K) = (tv%S(i,j,k-1) - tv%S(i,j,k)) * (dSpv_dS(i) * (alt_H_to_pres * I_dzval))
+            N2_temp(i,j,K) = (tv%T(i,j,k-1) - tv%T(i,j,k)) * (dSpV_dT(i) * (alt_H_to_pres * I_dzval))
+          enddo
+        enddo
+      enddo
+    else
+      !$OMP parallel do default(shared)
+      do j=js,je
+        if (associated(tv%p_surf)) then
+          do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo
+        else
+          do i=is,ie ; p_i(i) = 0.0 ; enddo
+        endif
+        do K=2,nz
+          do i=is,ie
+            p_i(i) = p_i(i) + H_to_pres* h(i,j,k-1)
+          enddo
+          T_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k))
+          S_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k))
+          call calculate_density_derivs(T_i, S_i, p_i, dRhodT, dRhodS, tv%eqn_of_state, EOSdom)
+          do i=is,ie
+            I_h = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))
+            N2_salt(i,j,K) = -(tv%S(i,j,k-1) - tv%S(i,j,k)) * (dRhodS(i) * (g_rho0 * I_h))
+            N2_temp(i,j,K) = -(tv%T(i,j,k-1) - tv%T(i,j,k)) * (dRhodT(i) * (g_rho0 * I_h))
+          enddo
+        enddo
+      enddo
+    endif
+    if (CS%id_N2_dd>0)      call post_data(CS%id_N2_dd, N2_salt(:,:,:)+N2_temp(:,:,:), CS%diag)
+    if (CS%id_N2_salt_dd>0) call post_data(CS%id_N2_salt_dd, N2_salt, CS%diag)
+    if (CS%id_N2_temp_dd>0) call post_data(CS%id_N2_temp_dd, N2_temp, CS%diag)
+
+    if (CS%Use_KdWork_diag) then
+       call KdWork_diagnostics(G,GV,US,CS%diag,CS%VBF,N2_salt,N2_temp,dz)
+    endif
+
+    call deallocate_VBF_CS(CS%VBF)
+
+  endif
+
   ! mixing of passive tracers from massless boundary layers to interior
   call cpu_clock_begin(id_clock_tracers)
 
@@ -1893,10 +2097,10 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
     call MOM_state_chksum("before set_diffusivity", u, v, h, G, GV, US, haloshift=CS%halo_TS_diff)
   if (CS%double_diffuse) then
     call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_int, G, GV, US, &
-                         CS%set_diff_CSp, Kd_lay=Kd_lay, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S)
+                         CS%set_diff_CSp, CS%VBF, Kd_lay=Kd_lay, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S)
   else
     call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_int, G, GV, US, &
-                         CS%set_diff_CSp, Kd_lay=Kd_lay)
+                         CS%set_diff_CSp, CS%VBF, Kd_lay=Kd_lay)
   endif
   call cpu_clock_end(id_clock_set_diffusivity)
   if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)")
@@ -3244,6 +3448,18 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
           Time, "Advective diapycnal salnity flux across interfaces", &
           units="psu m s-1", conversion=US%S_to_ppt*GV%H_to_m*US%s_to_T)
     endif
+
+    CS%id_N2_dd = register_diag_field('ocean_model',"N2_diabatic", diag%axesTi, &
+        Time, "Squared buoyancy frequency diagnosed after diffusion applied in thermodynamic timestep.", &
+        "s-2", conversion=US%s_to_T**2)
+    CS%id_N2_temp_dd = register_diag_field('ocean_model',"N2_temp_diabatic", diag%axesTi, &
+         Time, "Squared buoyancy frequency due to temperature stratification diagnosed after diffusion applied "//&
+         "in thermodynamic timestep.", "s-2", conversion=US%s_to_T**2)
+    CS%id_N2_salt_dd = register_diag_field('ocean_model',"N2_salt_diabatic", diag%axesTi, &
+         Time, "Squared buoyancy frequency due to salinity stratification diagnosed after diffusion applied "//&
+         "in thermodynamic timestep.", "s-2", conversion=US%s_to_T**2)
+    if (CS%id_N2_dd>0 .or. CS%id_N2_temp_dd>0 .or. CS%id_N2_salt_dd>0) CS%Use_N2_diag = .true.
+
     CS%id_MLD_003 = register_diag_field('ocean_model', 'MLD_003', diag%axesT1, Time, &
         'Mixed layer depth (delta rho = 0.03)', units='m', conversion=US%Z_to_m, &
         cmor_field_name='mlotst', cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', &
@@ -3284,6 +3500,13 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
            'Reference density for MLD (delta rho = 0.03)', units='kg/m3', conversion=US%R_to_kg_m3)
     endif
   endif
+
+  call KdWork_init(Time, G,GV,US,diag,CS%VBF,CS%Use_KdWork_diag)
+  if (CS%Use_KdWork_diag.and.(.not.useALEalgorithm)) &
+    call MOM_error(WARNING,"The KdWork diagnostics are not fully implemented for use in layer mode.")
+  if (CS%Use_KdWork_diag.and.(CS%use_legacy_diabatic)) &
+    call MOM_error(WARNING,"The KdWork diagnostics are only approximate with the legacy diabatic driver.")
+
   call get_param(param_file, mdl, "DIAG_MLD_DENSITY_DIFF", CS%MLDdensityDifference, &
                  "The density difference used to determine a diagnostic mixed "//&
                  "layer depth, MLD_user, following the definition of Levitus 1982. "//&
@@ -3607,6 +3830,10 @@ end subroutine register_diabatic_restarts
 subroutine diabatic_driver_end(CS)
   type(diabatic_CS), intent(inout) :: CS  !< module control structure
 
+  if (associated(CS%VBF)) then
+    call KdWork_end(CS%VBF)
+  endif
+
   if (associated(CS%optics)) then
     call opacity_end(CS%opacity, CS%optics)
     deallocate(CS%optics)
diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90
index 138a87754d..67e57c7cdf 100644
--- a/src/parameterizations/vertical/MOM_set_diffusivity.F90
+++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90
@@ -13,6 +13,7 @@ module MOM_set_diffusivity
 use MOM_CVMix_shear,         only : CVMix_shear_end
 use MOM_diag_mediator,       only : diag_ctrl, time_type
 use MOM_diag_mediator,       only : post_data, register_diag_field
+use MOM_diagnose_kdwork,     only : vbf_CS
 use MOM_debugging,           only : hchksum, uvchksum, Bchksum, hchksum_pair
 use MOM_EOS,                 only : calculate_density, calculate_density_derivs, EOS_domain
 use MOM_error_handler,       only : MOM_error, is_root_pe, FATAL, WARNING, NOTE
@@ -238,7 +239,7 @@ module MOM_set_diffusivity
 contains
 
 subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_int, &
-                           G, GV, US, CS, Kd_lay, Kd_extra_T, Kd_extra_S)
+                           G, GV, US, CS, VBF, Kd_lay, Kd_extra_T, Kd_extra_S)
   type(ocean_grid_type),     intent(in)    :: G    !< The ocean's grid structure.
   type(verticalGrid_type),   intent(in)    :: GV   !< The ocean's vertical grid structure.
   type(unit_scale_type),     intent(in)    :: US   !< A dimensional unit scaling type
@@ -264,9 +265,11 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
                              intent(out)   :: Kd_int !< Diapycnal diffusivity at each interface
                                                    !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
   type(set_diffusivity_CS),  pointer       :: CS   !< Module control structure.
+  type(vbf_CS),           pointer          :: VBF  !< A diagnostic control structure for vertical buoyancy fluxes
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                    optional, intent(out)   :: Kd_lay !< Diapycnal diffusivity of each layer
                                                    !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
+
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
                    optional, intent(out)   :: Kd_extra_T !< The extra diffusivity at interfaces of
                                                    !! temperature due to double diffusion relative
@@ -286,6 +289,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
 
   type(diffusivity_diags)  :: dd ! structure with arrays of available diags
 
+
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
     T_f, S_f      ! Temperature and salinity [C ~> degC] and [S ~> ppt] with properties in massless layers
                   ! filled vertically by diffusion or the properties after full convective adjustment.
@@ -380,7 +384,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
   if ((CS%double_diffusion) .and. (CS%id_KS_extra > 0)) &
     allocate(dd%KS_extra(isd:ied,jsd:jed,nz+1), source=0.0)
   if (CS%id_R_rho > 0) allocate(dd%drho_rat(isd:ied,jsd:jed,nz+1), source=0.0)
-  if (CS%id_Kd_BBL > 0) allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1), source=0.0)
+  if (CS%id_Kd_BBL > 0 .or. associated(VBF%Kd_BBL)) &
+    allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1), source=0.0)
 
   if (CS%id_Kd_bkgnd > 0) allocate(dd%Kd_bkgnd(isd:ied,jsd:jed,nz+1), source=0.)
   if (CS%id_Kv_bkgnd > 0) allocate(dd%Kv_bkgnd(isd:ied,jsd:jed,nz+1), source=0.)
@@ -432,6 +437,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
     endif
     call cpu_clock_end(id_clock_kappaShear)
     if (showCallTree) call callTree_waypoint("done with calculate_kappa_shear (set_diffusivity)")
+    if (associated(VBF%Kd_KS)) then ; do K=1,nz+1 ; do i=is,ie ; do j=js,je
+      VBF%Kd_KS(i,j,K) = visc%Kd_shear(i,j,K)
+    enddo ; enddo ; enddo ; endif
   elseif (CS%use_CVMix_shear) then
     !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)
@@ -486,6 +494,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
     if (CS%id_Kd_bkgnd > 0) then ; do K=1,nz+1 ; do i=is,ie
       dd%Kd_bkgnd(i,j,K) = Kd_int_2d(i,K)
     enddo ; enddo ; endif
+    if (associated(VBF%Kd_bkgnd)) then ; do K=1,nz+1 ; do i=is,ie
+      VBF%Kd_bkgnd(i,j,K) = Kd_int_2d(i,K)
+    enddo ; enddo ; endif
 
     ! Double-diffusion (old method)
     if (CS%double_diffusion) then
@@ -515,6 +526,13 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
       if (associated(dd%KS_extra)) then ; do K=1,nz+1 ; do i=is,ie
         dd%KS_extra(i,j,K) = KS_extra(i,K)
       enddo ; enddo ; endif
+
+      if (associated(VBF%Kd_ddiff_T)) then ; do K=1,nz+1 ; do i=is,ie
+        VBF%Kd_ddiff_T(i,j,K) = KT_extra(i,K)
+      enddo ; enddo ; endif
+      if (associated(VBF%Kd_ddiff_S)) then ; do K=1,nz+1 ; do i=is,ie
+        VBF%Kd_ddiff_S(i,j,K) = KS_extra(i,K)
+      enddo ; enddo ; endif ;
     endif
 
     ! Apply double diffusion via CVMix
@@ -527,6 +545,12 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
       else
         call compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_extra_T, Kd_extra_S, CS%CVMix_ddiff_csp)
       endif
+      if (associated(VBF%Kd_ddiff_T)) then ; do K=1,nz+1 ; do i=is,ie
+        VBF%Kd_ddiff_T(i,j,K) = KT_extra(i,K)
+      enddo ; enddo ; endif
+      if (associated(VBF%Kd_ddiff_S)) then ; do K=1,nz+1 ; do i=is,ie
+        VBF%Kd_ddiff_S(i,j,K) = KS_extra(i,K)
+      enddo ; enddo ; endif ;
       call cpu_clock_end(id_clock_CVMix_ddiff)
     endif
 
@@ -575,7 +599,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
     if (CS%use_tidal_mixing) &
       call calculate_tidal_mixing(dz, j, N2_bot, rho_bot, N2_lay, N2_int, TKE_to_Kd, &
                                   maxTKE, G, GV, US, CS%tidal_mixing, &
-                                  CS%Kd_max, visc%Kv_slow, Kd_lay_2d, Kd_int_2d)
+                                  CS%Kd_max, visc%Kv_slow, Kd_lay_2d, Kd_int_2d, VBF)
 
     ! Add diffusivity from internal tides ray tracing
     if (CS%use_int_tides) then
@@ -606,6 +630,21 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
       if (CS%id_Kd_slope > 0) then ; do K=1,nz+1 ; do i=is,ie
         dd%Kd_slope(i,j,K) = Kd_slope_2d(i,K)
       enddo ; enddo ; endif
+      if (associated (VBF%Kd_leak)) then ; do K=1,nz+1 ; do i=is,ie
+        VBF%Kd_leak(i,j,K) = Kd_leak_2d(i,K)
+      enddo ; enddo ; endif
+      if (associated (VBF%Kd_quad)) then ; do K=1,nz+1 ; do i=is,ie
+        VBF%Kd_quad(i,j,K) = Kd_quad_2d(i,K)
+      enddo ; enddo ; endif
+      if (associated (VBF%Kd_itidal)) then ; do K=1,nz+1 ; do i=is,ie
+        VBF%Kd_itidal(i,j,K) = Kd_itidal_2d(i,K)
+      enddo ; enddo ; endif
+      if (associated (VBF%Kd_Froude)) then ; do K=1,nz+1 ; do i=is,ie
+        VBF%Kd_Froude(i,j,K) = Kd_Froude_2d(i,K)
+      enddo ; enddo ; endif
+      if (associated (VBF%Kd_slope)) then ; do K=1,nz+1 ; do i=is,ie
+        VBF%Kd_slope(i,j,K) = Kd_slope_2d(i,K)
+      enddo ; enddo ; endif
 
       if (CS%id_prof_leak > 0) then ; do k=1,nz; do i=is,ie
         dd%prof_leak(i,j,k) = prof_leak_2d(i,k)
@@ -633,6 +672,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
         call add_drag_diffusivity(h, u, v,  tv, fluxes, visc, j, TKE_to_Kd, &
                                   maxTKE, kb, rho_bot, G, GV, US, CS, Kd_lay_2d, Kd_int_2d, dd%Kd_BBL)
       endif
+      if (associated(VBF%Kd_BBL)) then ; do K=1,nz+1 ; do i=is,ie
+        VBF%Kd_BBL(i,j,K) = dd%Kd_BBL(i,j,K)
+      enddo ; enddo ; endif
     endif
 
     if (CS%limit_dissipation) then
@@ -651,9 +693,12 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
     endif
 
     ! Optionally add a uniform diffusivity at the interfaces.
-    if (CS%Kd_add > 0.0) then ; do K=1,nz+1 ; do i=is,ie
-      Kd_int_2d(i,K) = Kd_int_2d(i,K) + CS%Kd_add
-    enddo ; enddo ; endif
+    if (CS%Kd_add > 0.0) then
+      do K=1,nz+1 ; do i=is,ie
+        Kd_int_2d(i,K) = Kd_int_2d(i,K) + CS%Kd_add
+      enddo; enddo
+      VBF%Kd_add = CS%Kd_add
+    endif
 
     ! Copy the 2-d slices into the 3-d array that is exported.
     do K=1,nz+1 ; do i=is,ie
diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90
index 5b57103078..b045fe6bde 100644
--- a/src/parameterizations/vertical/MOM_tidal_mixing.F90
+++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90
@@ -5,6 +5,7 @@ module MOM_tidal_mixing
 
 use MOM_diag_mediator,      only : diag_ctrl, time_type, register_diag_field
 use MOM_diag_mediator,      only : safe_alloc_ptr, post_data
+use MOM_diagnose_Kdwork,    only : vbf_CS
 use MOM_debugging,          only : hchksum
 use MOM_error_handler,      only : MOM_error, is_root_pe, FATAL, WARNING, NOTE
 use MOM_file_parser,        only : openParameterBlock, closeParameterBlock
@@ -695,7 +696,7 @@ end function tidal_mixing_init
 !! tidal dissipation and to add the effect of internal-tide-driven mixing to the layer or interface
 !! diffusivities.
 subroutine calculate_tidal_mixing(dz, j, N2_bot, Rho_bot, N2_lay, N2_int, TKE_to_Kd, max_TKE, &
-                                  G, GV, US, CS, Kd_max, Kv, Kd_lay, Kd_int)
+                                  G, GV, US, CS, Kd_max, Kv, Kd_lay, Kd_int, VBF)
   type(ocean_grid_type),            intent(in)    :: G      !< The ocean's grid structure
   type(verticalGrid_type),          intent(in)    :: GV     !< The ocean's vertical grid structure
   type(unit_scale_type),            intent(in)    :: US     !< A dimensional unit scaling type
@@ -729,13 +730,14 @@ subroutine calculate_tidal_mixing(dz, j, N2_bot, Rho_bot, N2_lay, N2_int, TKE_to
   real, dimension(SZI_(G),SZK_(GV)+1), &
                           optional, intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces
                                                             !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
+  type(vbf_CS), pointer                           :: VBF    !< A diagnostic structure for vertical buoyancy fluxes
 
   if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. CS%Lowmode_itidal_dissipation) then
     if (CS%use_CVMix_tidal) then
       call calculate_CVMix_tidal(dz, j, N2_int, G, GV, US, CS, Kv, Kd_lay, Kd_int)
     else
       call add_int_tide_diffusivity(dz, j, N2_bot, Rho_bot, N2_lay, TKE_to_Kd, max_TKE, &
-                                    G, GV, US, CS, Kd_max, Kd_lay, Kd_int)
+                                    G, GV, US, CS, Kd_max, Kd_lay, Kd_int, VBF)
     endif
   endif
 end subroutine calculate_tidal_mixing
@@ -992,7 +994,7 @@ end subroutine calculate_CVMix_tidal
 !! Will eventually need to add diffusivity due to other wave-breaking processes (e.g. Bottom friction,
 !! Froude-number-depending breaking, PSI, etc.).
 subroutine add_int_tide_diffusivity(dz, j, N2_bot, Rho_bot, N2_lay, TKE_to_Kd, max_TKE, &
-                                    G, GV, US, CS, Kd_max, Kd_lay, Kd_int)
+                                    G, GV, US, CS, Kd_max, Kd_lay, Kd_int, VBF)
   type(ocean_grid_type),             intent(in)    :: G      !< The ocean's grid structure
   type(verticalGrid_type),           intent(in)    :: GV     !< The ocean's vertical grid structure
   type(unit_scale_type),             intent(in)    :: US     !< A dimensional unit scaling type
@@ -1022,6 +1024,7 @@ subroutine add_int_tide_diffusivity(dz, j, N2_bot, Rho_bot, N2_lay, TKE_to_Kd, m
   real, dimension(SZI_(G),SZK_(GV)+1), &
                            optional, intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces
                                                              !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
+  type(vbf_CS), pointer                            :: VBF    !< A diagnostics structure for vertical buoyancy fluxes
 
   ! local
 
@@ -1302,38 +1305,57 @@ subroutine add_int_tide_diffusivity(dz, j, N2_bot, Rho_bot, N2_lay, TKE_to_Kd, m
       endif
 
       ! diagnostics
-      if (allocated(CS%dd%Kd_itidal)) then
+      if (allocated(CS%dd%Kd_itidal).or.(associated(VBF%Kd_itidal))) then
         ! If at layers, CS%dd%Kd_itidal is just TKE_to_Kd(i,k) * TKE_itide_lay
         ! The following sets the interface diagnostics.
         Kd_add = TKE_to_Kd(i,k) * TKE_itide_lay
         if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max)
-        if (k>1)  CS%dd%Kd_itidal(i,j,K)   = CS%dd%Kd_itidal(i,j,K)   + 0.5*Kd_add
-        if (k<nz) CS%dd%Kd_itidal(i,j,K+1) = CS%dd%Kd_itidal(i,j,K+1) + 0.5*Kd_add
+        if (allocated(CS%dd%Kd_itidal)) then
+          if (k>1)  CS%dd%Kd_itidal(i,j,K)   = CS%dd%Kd_itidal(i,j,K)   + 0.5*Kd_add
+          if (k<nz) CS%dd%Kd_itidal(i,j,K+1) = CS%dd%Kd_itidal(i,j,K+1) + 0.5*Kd_add
+        endif
+        if (associated(VBF%Kd_itides)) then
+          !Not to be confused w/ Kd_itidal (this is to be consistent w/ output parameter names)
+          if (k>1)  VBF%Kd_itides(i,j,K)   = VBF%Kd_itides(i,j,K)   + 0.5*Kd_add
+          if (k<nz) VBF%Kd_itides(i,j,K+1) = VBF%Kd_itides(i,j,K+1) + 0.5*Kd_add
+        endif
       endif
       if (allocated(CS%dd%Kd_Itidal_work)) &
         CS%dd%Kd_itidal_work(i,j,k) = GV%H_to_RZ * TKE_itide_lay
       if (allocated(CS%dd%Fl_itidal)) &
         CS%dd%Fl_itidal(i,j,k) = TKE_itidal_rem(i)
 
-      if (allocated(CS%dd%Kd_Niku)) then
+      if (allocated(CS%dd%Kd_Niku).or.(associated(VBF%Kd_Niku))) then
         ! If at layers, CS%dd%Kd_Niku(i,j,K) is just TKE_to_Kd(i,k) * TKE_Niku_lay
         ! The following sets the interface diagnostics.
         Kd_add = TKE_to_Kd(i,k) * TKE_Niku_lay
         if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max)
-        if (k>1) CS%dd%Kd_Niku(i,j,K)    = CS%dd%Kd_Niku(i,j,K)   + 0.5*Kd_add
-        if (k<nz) CS%dd%Kd_Niku(i,j,K+1) = CS%dd%Kd_Niku(i,j,K+1) + 0.5*Kd_add
+        if (allocated(CS%dd%Kd_Niku)) then
+          if (k>1) CS%dd%Kd_Niku(i,j,K)    = CS%dd%Kd_Niku(i,j,K)   + 0.5*Kd_add
+          if (k<nz) CS%dd%Kd_Niku(i,j,K+1) = CS%dd%Kd_Niku(i,j,K+1) + 0.5*Kd_add
+        endif
+        if (associated(VBF%Kd_Niku)) then
+          if (k>1) VBF%Kd_Niku(i,j,K)    = VBF%Kd_Niku(i,j,K)   + 0.5*Kd_add
+          if (k<nz) VBF%Kd_Niku(i,j,K+1) = VBF%Kd_Niku(i,j,K+1) + 0.5*Kd_add
+        endif
       endif
 !     if (associated(CS%dd%Kd_Niku)) CS%dd%Kd_Niku(i,j,K) = TKE_to_Kd(i,k) * TKE_Niku_lay
       if (allocated(CS%dd%Kd_Niku_work)) &
         CS%dd%Kd_Niku_work(i,j,k) = GV%H_to_RZ * TKE_Niku_lay
 
-      if (allocated(CS%dd%Kd_lowmode)) then
+      if (allocated(CS%dd%Kd_lowmode).or.(associated(VBF%Kd_lowmode))) then
         ! If at layers, CS%dd%Kd_lowmode is just TKE_to_Kd(i,k) * TKE_lowmode_lay
         ! The following sets the interface diagnostics.
         Kd_add = TKE_to_Kd(i,k) * TKE_lowmode_lay
         if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max)
-        if (k>1)  CS%dd%Kd_lowmode(i,j,K)   = CS%dd%Kd_lowmode(i,j,K)   + 0.5*Kd_add
-        if (k<nz) CS%dd%Kd_lowmode(i,j,K+1) = CS%dd%Kd_lowmode(i,j,K+1) + 0.5*Kd_add
+        if (allocated(CS%dd%Kd_lowmode)) then
+          if (k>1)  CS%dd%Kd_lowmode(i,j,K)   = CS%dd%Kd_lowmode(i,j,K)   + 0.5*Kd_add
+          if (k<nz) CS%dd%Kd_lowmode(i,j,K+1) = CS%dd%Kd_lowmode(i,j,K+1) + 0.5*Kd_add
+        endif
+        if (associated(VBF%Kd_lowmode)) then
+          if (k>1)  VBF%Kd_lowmode(i,j,K)   = VBF%Kd_lowmode(i,j,K)   + 0.5*Kd_add
+          if (k<nz) VBF%Kd_lowmode(i,j,K+1) = VBF%Kd_lowmode(i,j,K+1) + 0.5*Kd_add
+        endif
       endif
       if (allocated(CS%dd%Kd_lowmode_work)) &
         CS%dd%Kd_lowmode_work(i,j,k) = GV%H_to_RZ * TKE_lowmode_lay