diff --git a/src/core/MOM_unit_tests.F90 b/src/core/MOM_unit_tests.F90
index bd449d0b39..b817493b44 100644
--- a/src/core/MOM_unit_tests.F90
+++ b/src/core/MOM_unit_tests.F90
@@ -3,6 +3,7 @@ module MOM_unit_tests
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+use MOM_diag_buffers,               only : diag_buffer_unit_tests_2d, diag_buffer_unit_tests_3d
 use MOM_error_handler,              only : MOM_error, FATAL, is_root_pe
 use MOM_hor_bnd_diffusion,          only : near_boundary_unit_tests
 use MOM_intrinsic_functions,        only : intrinsic_functions_unit_tests
@@ -49,6 +50,10 @@ subroutine unit_tests(verbosity)
        "MOM_unit_tests: CFC_cap_unit_tests FAILED")
     if (mixedlayer_restrat_unit_tests(verbose)) call MOM_error(FATAL, &
        "MOM_unit_tests: mixedlayer_restrat_unit_tests FAILED")
+    if (diag_buffer_unit_tests_2d(verbose)) call MOM_error(FATAL, &
+       "MOM_unit_tests: diag_buffer_unit_tests_2d FAILED")
+    if (diag_buffer_unit_tests_3d(verbose)) call MOM_error(FATAL, &
+       "MOM_unit_tests: diag_buffer_unit_tests_3d FAILED")
   endif
 
 end subroutine unit_tests
diff --git a/src/framework/MOM_diag_buffers.F90 b/src/framework/MOM_diag_buffers.F90
new file mode 100644
index 0000000000..246d1c0437
--- /dev/null
+++ b/src/framework/MOM_diag_buffers.F90
@@ -0,0 +1,526 @@
+!> Provides buffers that can dynamically grow as needed. These are primarily intended for the
+!! diagnostics which need to store intermediate or partial states of state variables
+module MOM_diag_buffers
+
+use MOM_io, only : stdout, stderr
+
+! This file is part of MOM6. See LICENSE.md for the license.
+
+implicit none ; private
+
+public :: diag_buffer_unit_tests_2d, diag_buffer_unit_tests_3d
+
+type, abstract :: buffer_base
+end type buffer_base
+
+!> Holds a 2d field
+type, extends(buffer_base) :: buffer_2d
+  real, dimension(:,:), allocatable :: field !< The actual 2d field to be stored [arbitrary]
+end type buffer_2d
+
+!> Holds a 3d field
+type, extends(buffer_base) :: buffer_3d
+  real, dimension(:,:,:), allocatable :: field !< The actual 3d field to be stored [arbitrary]
+end type buffer_3d
+
+!> The base class for the diagnostic buffers in this module
+type, abstract :: diag_buffer_base ; private
+  integer :: is !< The start slot of the array i-direction
+  integer :: js !< The start slot of the array j-direction
+  integer :: ie !< The end slot of the array i-direction
+  integer :: je !< The end slot of the array j-direction
+  real :: fill_value = 0. !< Set the fill value to use when growing the buffer [arbitrary]
+
+  integer, allocatable, dimension(:) :: ids  !< List of diagnostic ids whose slot corresponds to the row in the buffer
+  integer :: length = 0 !< The number of slots in the buffer
+
+  contains
+
+  procedure(a_grow), deferred :: grow !< Increase the size of the buffer
+  procedure, public :: set_fill_value !< Set the fill value to use when growing the buffer
+  procedure, public :: check_capacity_by_id !< Check the size size of the buffer and increase if necessary
+  procedure, public :: set_horizontal_extents !< Define the horizontal extents of the arrays
+  procedure, public :: mark_available !< Mark that a slot in the buffer can be reused
+  procedure, public :: grow_ids !< Increase the size of the vector storing diagnostic ids
+  procedure, public :: find_buffer_slot !< Find the slot corresponding to a specific diagnostic id
+end type diag_buffer_base
+
+!> Dynamically growing buffer for 2D arrays.
+type, extends(diag_buffer_base), public :: diag_buffer_2d; private
+  type(buffer_2d), public, dimension(:), allocatable :: buffer !< The actual 2D buffer which will dynamically grow
+
+  contains
+
+  procedure, public :: grow => grow_2d !< Increase the size of the buffer
+  procedure, public :: store => store_2d !< Store a field in the buffer, increasing as necessary
+end type diag_buffer_2d
+
+!> Dynamically growing buffer for 3D arrays.
+type, extends(diag_buffer_base), public :: diag_buffer_3d ; private
+  type(buffer_3d), public, dimension(:), allocatable :: buffer !< The actual 2D buffer which will dynamically grow
+  integer :: ks !< The start slot in the k-dimension
+  integer :: ke !< The last slot in the k-dimension
+
+  contains
+
+  procedure, public :: set_vertical_extent !< Set the vertical extents of the buffer
+  procedure, public :: grow => grow_3d !< Increase the size of the buffer
+  procedure, public :: store => store_3d !< Store a field in the buffer, increasing as necessary
+end type diag_buffer_3d
+
+contains
+
+!> Signature for the grow methods on n-dimension diagnostic buffer types
+subroutine a_grow(this)
+  class(diag_buffer_base), intent(inout) :: this !< The diagnostic buffer
+end subroutine
+
+!> Set the fill value to use when growing the buffer
+subroutine set_fill_value(this, fill_value)
+  class(diag_buffer_base), intent(inout) :: this !< The diagnostic buffer
+  real,                    intent(in)    :: fill_value !< The fill value to use when growing the buffer [arbitrary]
+
+  this%fill_value = fill_value
+end subroutine set_fill_value
+
+!> Mark a slot in the buffer as unused based on a diagnostic id. For example,
+!! the data in that slot has already been consumed and can thus be overwritten
+subroutine mark_available(this, id)
+  class(diag_buffer_base), intent(inout) :: this !< The diagnostic buffer
+  integer,                 intent(in)    :: id   !< The diagnostic id
+  integer :: slot
+
+  slot = this%find_buffer_slot(id)
+  this%ids(slot) = 0
+end subroutine mark_available
+
+!> Return the slot of the buffer corresponding to the diagnostic id
+pure function find_buffer_slot(this, id) result(slot)
+  class(diag_buffer_base), intent(in) :: this !< The diagnostic buffer
+  integer, intent(in) :: id !< The diagnostic id
+
+  integer, dimension(1) :: temp
+  integer :: slot !< The slot in the buffer corresponding to the diagnostic id
+
+  if (allocated(this%ids)) then
+    !NOTE: Alternatively could do slot = SUM(findloc(...))
+    temp = findloc(this%ids(:), id)
+    slot = temp(1)
+  else
+    slot = 0
+  endif
+
+end function find_buffer_slot
+
+!> Grow the ids array by one
+subroutine grow_ids(this)
+  class(diag_buffer_base), intent(inout) :: this !< This buffer
+
+  integer, allocatable, dimension(:) :: temp
+  integer :: n
+
+  n = this%length
+
+  allocate(temp(n+1))
+  if (n>0) temp(1:n) = this%ids(:)
+  call move_alloc(temp, this%ids)
+end subroutine grow_ids
+
+!> Check whether the id already has a slot reserved. If not, find a new empty slot and if
+!! need be, grow the buffer.
+impure function check_capacity_by_id(this, id) result(slot)
+  class(diag_buffer_base), intent(inout) :: this !< This 2d buffer
+  integer,                 intent(in)    :: id   !< The diagnostic id
+  integer :: slot
+
+  slot = this%find_buffer_slot(id)
+  if (slot==0) then
+    ! Check to see if there is an open slot
+    if (allocated(this%ids)) slot = this%find_buffer_slot(0)
+    ! If slot is still 0, then the buffer must grow
+    if (slot==0) then
+      call this%grow()
+      slot = this%length
+    endif
+    this%ids(slot) = id
+  endif
+end function check_capacity_by_id
+
+!> Set the horizontal extents of the buffer
+subroutine set_horizontal_extents(this, is, ie, js, je)
+  class(diag_buffer_base), intent(inout) :: this !< The diagnostic buffer
+  integer,               intent(in)    :: is !< The start slot of the array i-direction
+  integer,               intent(in)    :: ie !< The end slot of the array i-direction
+  integer,               intent(in)    :: js !< The start slot of the array j-direction
+  integer,               intent(in)    :: je !< The end slot of the array j-direction
+
+  this%is = is ; this%ie = ie ; this%js = js ; this%je = je
+end subroutine set_horizontal_extents
+
+!> Set the vertical extent of the buffer
+subroutine set_vertical_extent(this, ks, ke)
+  class(diag_buffer_3d), intent(inout) :: this !< The diagnostic buffer
+  integer,               intent(in)    :: ks !< The start slot of the array i-direction
+  integer,               intent(in)    :: ke !< The end slot of the array i-direction
+
+  this%ks = ks; this%ke = ke
+end subroutine set_vertical_extent
+
+!> Grow a 2d diagnostic buffer
+subroutine grow_2d(this)
+  class(diag_buffer_2d), intent(inout) :: this
+
+  integer :: i, n
+  integer :: is, ie, js, je
+  type(buffer_2d), dimension(:), allocatable :: new_buffer
+
+  ! Grow the ID array
+  call this%grow_ids()
+
+  is = this%is; ie=this%ie; js=this%js; je=this%je
+  n = this%length
+
+  allocate(new_buffer(n+1))
+  do i=1,n
+    allocate(new_buffer(i)%field(is:ie,js:je))
+    new_buffer(i)%field(:,:) = this%buffer(i)%field(:,:)
+  enddo
+  allocate(new_buffer(n+1)%field(is:ie,js:je), source=this%fill_value)
+  call move_alloc(new_buffer, this%buffer)
+  this%length = n+1
+
+end subroutine grow_2d
+
+!> Store a 2D array into this buffer
+subroutine store_2d(this, data, id)
+  class(diag_buffer_2d), intent(inout) :: this !< This 2d buffer
+  real, dimension(:,:),  intent(in)    :: data !< The data to be stored in the buffer [arbitrary]
+  integer,               intent(in)    :: id !< The diagnostic id
+
+  integer :: slot
+
+  slot = this%check_capacity_by_id(id)
+  this%buffer(slot)%field(:,:) = data(:,:)
+end subroutine store_2d
+
+!> Grow a 2d diagnostic buffer
+subroutine grow_3d(this)
+  class(diag_buffer_3d), intent(inout) :: this
+
+  integer :: i, n
+  integer :: is, ie, js, je, ks, ke
+  type(buffer_3d), dimension(:), allocatable :: new_buffer
+
+  ! Grow the ID array
+  call this%grow_ids()
+
+  is = this%is; ie=this%ie; js=this%js; je=this%je; ks=this%ks; ke=this%ke
+  n = this%length
+
+  allocate(new_buffer(n+1))
+  do i=1,n
+    allocate(new_buffer(i)%field(is:ie,js:je,ks:ke))
+    new_buffer(i)%field(:,:,:) = this%buffer(i)%field(:,:,:)
+  enddo
+  allocate(new_buffer(n+1)%field(is:ie,js:je,ks:ke), source=this%fill_value)
+  call move_alloc(new_buffer, this%buffer)
+  this%length = n+1
+
+end subroutine grow_3d
+
+!> Store a 3d array into this buffer
+subroutine store_3d(this, data, id)
+  class(diag_buffer_3d),  intent(inout) :: this !< This 3d buffer
+  real, dimension(:,:,:), intent(in)    :: data !< The data to be stored in the buffer [arbitrary]
+  integer,                intent(in)    :: id !< The diagnostic id
+
+  integer :: slot
+
+  ! Find the first slot in the ids array that is 0, i.e. this is a portion of the buffer that can be reused
+  slot = this%check_capacity_by_id(id)
+  this%buffer(slot)%field(:,:,:) = data(:,:,:)
+end subroutine store_3d
+
+!> Unit tests for the 2d version of the diag buffer
+function diag_buffer_unit_tests_2d(verbose) result(fail)
+  logical, intent(in) :: verbose !< If true, write results to stdout
+  logical :: fail !< True if any of the unit tests fail
+
+  fail = .false.
+  write(stdout,*) '==== MOM_diag_buffers: diag_buffers_unit_tests_2d ==='
+  fail = fail .or. new_buffer_2d()
+  fail = fail .or. grow_buffer_2d()
+  fail = fail .or. fill_value_2d()
+  fail = fail .or. store_buffer_2d()
+  fail = fail .or. reuse_buffer_2d()
+
+  contains
+
+  !> Ensure properties of a newly initialized buffer
+  function new_buffer_2d() result(local_fail)
+    type(diag_buffer_2d) :: buffer
+    logical :: local_fail !< True if any of the unit tests fail
+    local_fail = .false.
+    local_fail = local_fail .or. allocated(buffer%buffer)
+    if (verbose) write(stdout,*) "new_buffer_2d: ", local_fail
+    local_fail = local_fail .or. allocated(buffer%ids)
+    if (verbose) write(stdout,*) "new_buffer_2d: ", local_fail
+    local_fail = local_fail .or. buffer%length /= 0
+    if (verbose) write(stdout,*) "new_buffer_2d: ", local_fail
+  end function new_buffer_2d
+
+  !> Test the growing of a buffer
+  function grow_buffer_2d() result(local_fail)
+    type(diag_buffer_2d) :: buffer
+    logical :: local_fail !< True if any of the unit tests fail
+    integer, parameter :: is=1, ie=2, js=3, je=6
+    integer :: i
+
+    local_fail = .false.
+
+    call buffer%set_horizontal_extents(is=is, ie=ie, js=js, je=je)
+    ! Grow the buffer 3 times
+    do i=1,3
+      call buffer%grow()
+      local_fail = local_fail .or. (buffer%length /= i)
+      local_fail = local_fail .or. (lbound(buffer%buffer(i)%field, 1) /= is)
+      local_fail = local_fail .or. (ubound(buffer%buffer(i)%field, 1) /= ie)
+      local_fail = local_fail .or. (lbound(buffer%buffer(i)%field, 2) /= js)
+      local_fail = local_fail .or. (ubound(buffer%buffer(i)%field, 2) /= je)
+    enddo
+    if (verbose) write(stdout,*) "grow_buffer_2d: ", local_fail
+  end function grow_buffer_2d
+
+  !> Test that growing new buffer fills the array with a set fill value
+  function fill_value_2d() result(local_fail)
+    type(diag_buffer_2d) :: buffer
+    logical :: local_fail !< True if any of the unit tests fail
+    integer, parameter :: is=1, ie=2, js=3, je=6
+    real, parameter :: fill_value = -123.456
+    integer :: i
+
+
+    local_fail = .false.
+
+    call buffer%set_horizontal_extents(is=is, ie=ie, js=js, je=je)
+    call buffer%set_fill_value(fill_value)
+    ! Grow the buffer 3 times
+    call buffer%grow()
+    if (any(buffer%buffer(1)%field(:,:) /= fill_value)) local_fail = .true.
+    if (verbose) write(stdout,*) "fill_value_2d: ", local_fail
+  end function fill_value_2d
+
+  !> Test storing a buffer based on a unique id
+  function store_buffer_2d() result(local_fail)
+    type(diag_buffer_2d) :: buffer
+    logical :: local_fail !< True if any of the unit tests fail
+
+    integer, parameter :: is=1, ie=2, js=3, je=6, nlen=3
+    integer :: i, slot
+    real, allocatable, dimension(:,:,:) :: test_2d
+
+    local_fail = .false.
+
+    allocate(test_2d(nlen, is:ie, js:je))
+    call random_number(test_2d)
+    buffer%is = is
+    buffer%ie = ie
+    buffer%js = js
+    buffer%je = je
+
+    do i=1,nlen
+      call buffer%store(test_2d(i,:,:), i*3)
+      slot = buffer%find_buffer_slot(i*3)
+      local_fail = local_fail .or. ANY(buffer%buffer(slot)%field(:,:) /= test_2d(i,:,:))
+    enddo
+
+    if (verbose) write(stdout,*) "store_buffer_2d: ", local_fail
+  end function store_buffer_2d
+
+  !> Test the reuse of a buffer. Fill it first like store_buffer_2d. Then,
+  !! loop through again, but use the slots of the buffer in the following
+  !! order: 2, 1, 3
+  function reuse_buffer_2d() result(local_fail)
+    type(diag_buffer_2d) :: buffer
+    logical :: local_fail !< True if any of the unit tests fail
+
+    integer, parameter :: is=1, ie=2, js=3, je=6, nlen=3
+    integer :: i, new_i, id, new_id
+    real, dimension(nlen, is:ie, js:je) :: test_2d_first, test_2d_second
+    integer, dimension(nlen) :: reorder = [2,1,3]
+
+    local_fail = .false.
+    call random_number(test_2d_first)
+    call random_number(test_2d_second)
+
+    call buffer%set_horizontal_extents(is=is, ie=ie, js=js, je=je)
+
+    do i=1,nlen
+      call buffer%store(test_2d_first(i,:,:), id=i*3)
+    enddo
+
+    do i=1,nlen
+      new_i = reorder(i)
+      ! id and new_id are multiplied by primes to make sure they are unique
+      id = reorder(i)*3
+      new_id = i*7
+      call buffer%mark_available(id=reorder(i)*3)
+      call buffer%store(test_2d_second(i,:,:), id=new_id)
+      local_fail = local_fail .or. buffer%find_buffer_slot(new_id) /= new_i
+      test_2d_first(new_i,:,:) = test_2d_second(i,:,:)
+    enddo
+    local_fail = local_fail .or. any(buffer%ids /= [14, 7, 21])
+    do i=1,nlen
+      local_fail = local_fail .or. any(buffer%buffer(i)%field(:,:) /= test_2d_first(i,:,:))
+    enddo
+    if (verbose) write(stdout,*) "reuse_buffer_2d: ", local_fail
+  end function reuse_buffer_2d
+
+end function diag_buffer_unit_tests_2d
+
+!> Test the 3d version of the buffer
+function diag_buffer_unit_tests_3d(verbose) result(fail)
+  logical, intent(in) :: verbose !< If true, write results to stdout
+  logical :: fail !< True if any of the unit tests fail
+
+  fail = .false.
+  write(stdout,*) '==== MOM_diag_buffers: diag_buffers_unit_tests_3d ==='
+  fail = fail .or. new_buffer_3d()
+  fail = fail .or. grow_buffer_3d()
+  fail = fail .or. fill_value_3d()
+  fail = fail .or. store_buffer_3d()
+  fail = fail .or. reuse_buffer_3d()
+
+  contains
+
+  !> Ensure properties of a newly initialized buffer
+  function new_buffer_3d() result(local_fail)
+    type(diag_buffer_3d) :: buffer
+    logical :: local_fail !< True if any of the unit tests fail
+    local_fail = .false.
+    local_fail = local_fail .or. allocated(buffer%buffer)
+    local_fail = local_fail .or. allocated(buffer%ids)
+    local_fail = local_fail .or. buffer%length /= 0
+    if (verbose) write(stdout,*) "new_buffer_3d: ", local_fail
+  end function new_buffer_3d
+
+  !> Test the growing of a buffer
+  function grow_buffer_3d() result(local_fail)
+    type(diag_buffer_3d) :: buffer
+    logical :: local_fail !< True if any of the unit tests fail
+    integer, parameter :: is=1, ie=2, js=3, je=6, ks=1, ke=10
+    integer :: i
+
+    local_fail = .false.
+
+    call buffer%set_horizontal_extents(is=is, ie=ie, js=js, je=je)
+    call buffer%set_vertical_extent(ks=ks, ke=ke)
+    ! Grow the buffer 3 times
+    do i=1,3
+      call buffer%grow()
+      local_fail = local_fail .or. (buffer%length /= i)
+      local_fail = local_fail .or. (lbound(buffer%buffer(i)%field, 1) /= is)
+      local_fail = local_fail .or. (ubound(buffer%buffer(i)%field, 1) /= ie)
+      local_fail = local_fail .or. (lbound(buffer%buffer(i)%field, 2) /= js)
+      local_fail = local_fail .or. (ubound(buffer%buffer(i)%field, 2) /= je)
+      local_fail = local_fail .or. (lbound(buffer%buffer(i)%field, 3) /= ks)
+      local_fail = local_fail .or. (ubound(buffer%buffer(i)%field, 3) /= ke)
+    if (verbose) write(stdout,*) "grow_buffer_3d: ", local_fail
+    enddo
+    if (verbose) write(stdout,*) "grow_buffer_3d: ", local_fail
+  end function grow_buffer_3d
+
+  !> Test that growing new buffer fills the array with a set fill value
+  function fill_value_3d() result(local_fail)
+    type(diag_buffer_3d) :: buffer
+    logical :: local_fail !< True if any of the unit tests fail
+    integer, parameter :: is=1, ie=2, js=3, je=6
+    real, parameter :: fill_value = -123.456
+    integer :: i
+
+
+    local_fail = .false.
+
+    call buffer%set_horizontal_extents(is=is, ie=ie, js=js, je=je)
+    call buffer%set_fill_value(fill_value)
+    ! Grow the buffer 3 times
+    call buffer%grow()
+    if (any(buffer%buffer(1)%field(:,:,:) /= fill_value)) local_fail = .true.
+    if (verbose) write(stdout,*) "fill_value_3d: ", local_fail
+  end function fill_value_3d
+
+  !> Test storing a buffer based on a unique id
+  function store_buffer_3d() result(local_fail)
+    type(diag_buffer_3d) :: buffer
+    logical :: local_fail !< True if any of the unit tests fail
+
+    integer, parameter :: is=1, ie=2, js=3, je=6, ks=1, ke=10, nlen=3
+    integer :: i, slot
+    real, dimension(nlen,is:ie,js:je,ks:ke) :: test_3d
+
+    local_fail = .false.
+    call random_number(test_3d)
+    buffer%is = is
+    buffer%ie = ie
+    buffer%js = js
+    buffer%je = je
+    buffer%ks = ks
+    buffer%ke = ke
+
+    do i=1,nlen
+      call buffer%store(test_3d(i,:,:,:), i*3)
+      slot = buffer%find_buffer_slot(i*3)
+      local_fail = local_fail .or. ANY(buffer%buffer(slot)%field(:,:,:) /= test_3d(i,:,:,:))
+    enddo
+
+    if (verbose) write(stdout,*) "store_buffer_3d: ", local_fail
+  end function store_buffer_3d
+
+  !> Test the reuse of a buffer. Fill it first like store_buffer_3d. Then,
+  !! loop through again, but use the slots of the buffer in the following
+  !! order: 2, 1, 3
+  function reuse_buffer_3d() result(local_fail)
+    type(diag_buffer_3d) :: buffer
+    logical :: local_fail !< True if any of the unit tests fail
+
+    integer, parameter :: is=1, ie=2, js=3, je=6, ks=1, ke=10, nlen=3
+    integer :: i, new_i, id, new_id
+    real, dimension(nlen, is:ie, js:je, ks:ke) :: test_3d_first, test_3d_second
+    integer, dimension(nlen) :: reorder = [2,1,3]
+
+    local_fail = .false.
+    call random_number(test_3d_first)
+    call random_number(test_3d_second)
+
+    buffer%is = is
+    buffer%ie = ie
+    buffer%js = js
+    buffer%je = je
+    buffer%ks = ks
+    buffer%ke = ke
+
+    do i=1,nlen
+      call buffer%store(test_3d_first(i,:,:,:), id=i*3)
+    enddo
+
+    do i=1,nlen
+      new_i = reorder(i)
+      ! id and new_id are multiplied by primes to make sure they are unique
+      id = reorder(i)*3
+      new_id = i*7
+      call buffer%mark_available(id=reorder(i)*3)
+      call buffer%store(test_3d_second(i,:,:,:), id=new_id)
+      local_fail = local_fail .or. buffer%find_buffer_slot(new_id) /= new_i
+      test_3d_first(new_i,:,:,:) = test_3d_second(i,:,:,:)
+    enddo
+    local_fail = local_fail .or. any(buffer%ids /= [14, 7, 21])
+    do i=1,nlen
+      local_fail = local_fail .or. any(buffer%buffer(i)%field(:,:,:) /= test_3d_first(i,:,:,:))
+    enddo
+    if (verbose) write(stdout,*) "reuse_buffer_3d: ", local_fail
+  end function reuse_buffer_3d
+
+end function diag_buffer_unit_tests_3d
+
+end module MOM_diag_buffers
+
diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90
index 1a43739147..cdb1e2dc42 100644
--- a/src/framework/MOM_diag_mediator.F90
+++ b/src/framework/MOM_diag_mediator.F90
@@ -9,6 +9,7 @@ module MOM_diag_mediator
 use MOM_coms,             only : PE_here
 use MOM_cpu_clock,        only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
 use MOM_cpu_clock,        only : CLOCK_MODULE, CLOCK_ROUTINE
+use MOM_diag_buffers,     only : diag_buffer_2d, diag_buffer_3d
 use MOM_diag_manager_infra, only : MOM_diag_manager_init, MOM_diag_manager_end
 use MOM_diag_manager_infra, only : diag_axis_init=>MOM_diag_axis_init, get_MOM_diag_axis_name
 use MOM_diag_manager_infra, only : send_data_infra, MOM_diag_field_add_attribute, EAST, NORTH
@@ -43,6 +44,7 @@ module MOM_diag_mediator
 #define MAX_DSAMP_LEV 2
 
 public set_axes_info, post_data, register_diag_field, time_type
+public post_data_3d_by_column, post_data_3d_final
 public post_product_u, post_product_sum_u, post_product_v, post_product_sum_v
 public set_masks_for_axes
 ! post_data_1d_k is a deprecated interface that can be replaced by a call to post_data, but
@@ -134,6 +136,10 @@ module MOM_diag_mediator
   real, pointer, dimension(:,:)   :: mask2d => null() !< Mask for 2d (x-y) axes [nondim]
   real, pointer, dimension(:,:,:) :: mask3d => null() !< Mask for 3d axes [nondim]
   type(diag_dsamp), dimension(2:MAX_DSAMP_LEV) :: dsamp !< Downsample container
+
+  ! For diagnostics posted piecemeal
+  type(diag_buffer_2d) :: piecemeal_2d !< A dynamically reallocated buffer for 2d piecemeal diagnostics
+  type(diag_buffer_3d) :: piecemeal_3d !< A dynamically reallocated buffer for 3d piecemeal diagnostics
 end type axes_grp
 
 !> Contains an array to store a diagnostic target grid
@@ -1102,6 +1108,9 @@ subroutine define_axes_group(diag_cs, handles, axes, nz, vertical_coordinate_num
     if (axes%is_u_point) axes%mask2d => diag_cs%mask2dCu
     if (axes%is_v_point) axes%mask2d => diag_cs%mask2dCv
     if (axes%is_q_point) axes%mask2d => diag_cs%mask2dBu
+    call axes%piecemeal_2d%set_horizontal_extents(lbound(axes%mask2d,1), ubound(axes%mask2d,1), &
+                                               lbound(axes%mask2d,2), ubound(axes%mask2d,2))
+    call axes%piecemeal_2d%set_fill_value(diag_cs%missing_value)
   endif
   ! A static 3d mask for non-native coordinates can only be setup when a grid is available
   axes%mask3d => null()
@@ -1118,8 +1127,13 @@ subroutine define_axes_group(diag_cs, handles, axes, nz, vertical_coordinate_num
       if (axes%is_v_point) axes%mask3d => diag_cs%mask3dCvi
       if (axes%is_q_point) axes%mask3d => diag_cs%mask3dBi
     endif
+    call axes%piecemeal_3d%set_horizontal_extents(is=lbound(axes%mask3d,1), ie=ubound(axes%mask3d,1), &
+                                               js=lbound(axes%mask3d,2), je=ubound(axes%mask3d,2))
+    call axes%piecemeal_3d%set_vertical_extent(ks=lbound(axes%mask3d,3), ke=ubound(axes%mask3d,3))
+    call axes%piecemeal_3d%set_fill_value(diag_cs%missing_value)
   endif
 
+
 end subroutine define_axes_group
 
 !> Defines a group of downsampled "axes" from list of handles
@@ -1922,6 +1936,58 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask)
 
 end subroutine post_data_3d_low
 
+!> Put data into the buffer for a diagnostic one column at a time
+subroutine post_data_3d_by_column(diag_field_id, field, diag_cs, i, j)
+  integer,            intent(in) :: diag_field_id !< The id for an output variable returned by a
+                                                  !! previous call to register_diag_field.
+  real, dimension(:), intent(in) :: field         !< 3-d array being offered for output or averaging
+                                                  !! in internally scaled arbitrary units [A ~> a]
+  type(diag_ctrl), target, intent(in) :: diag_CS  !< Structure used to regulate diagnostic output
+  integer,            intent(in) :: i             !< The i-index to post the data in the buffer
+  integer,            intent(in) :: j             !< The j-index to post the data in the buffer
+
+  type(diag_type), pointer :: diag => null()
+  integer :: buffer_slot
+
+  diag => diag_cs%diags(diag_field_id)
+  buffer_slot = diag%axes%piecemeal_3d%check_capacity_by_id(diag_field_id)
+  diag%axes%piecemeal_3d%buffer(buffer_slot)%field(i,j,:) = field(:)
+end subroutine post_data_3d_by_column
+
+!> Put data into the buffer for a diagnostic one point at a time
+subroutine post_data_3d_by_point(diag_field_id, field, diag_cs, i, j, k)
+  integer,           intent(in) :: diag_field_id !< The id for an output variable returned by a
+                                                 !! previous call to register_diag_field.
+  real,              intent(in) :: field         !< 3-d array being offered for output or averaging
+                                                 !! in internally scaled arbitrary units [A ~> a]
+  type(diag_ctrl), target, intent(in) :: diag_CS !< Structure used to regulate diagnostic output
+  integer,           intent(in) :: i             !< The i-index to post the data in the buffer
+  integer,           intent(in) :: j             !< The j-index to post the data in the buffer
+  integer,           intent(in) :: k             !< The k-index to post the data in the buffer
+
+  type(diag_type), pointer :: diag => null()
+  integer :: buffer_slot
+
+  diag => diag_cs%diags(diag_field_id)
+  buffer_slot = diag%axes%piecemeal_3d%find_buffer_slot(diag_field_id)
+  diag%axes%piecemeal_3d%buffer(buffer_slot)%field(i,j,k) = field
+end subroutine post_data_3d_by_point
+
+!> Post the final buffer using the standard post_data interface
+subroutine post_data_3d_final(diag_field_id, diag_cs)
+  integer,           intent(in) :: diag_field_id !< The id for an output variable returned by a
+                                                 !! previous call to register_diag_field.
+  type(diag_ctrl), target, intent(in) :: diag_CS !< Structure used to regulate diagnostic output
+
+  type(diag_type), pointer :: diag => null()
+  integer :: buffer_slot
+
+  diag => diag_cs%diags(diag_field_id)
+  buffer_slot = diag%axes%piecemeal_3d%find_buffer_slot(diag_field_id)
+  call post_data(diag_field_id, diag%axes%piecemeal_3d%buffer(buffer_slot)%field(:,:,:), diag_CS)
+  call diag%axes%piecemeal_3d%mark_available(diag_field_id)
+end subroutine post_data_3d_final
+
 !> Calculate and write out diagnostics that are the product of two 3-d arrays at u-points
 subroutine post_product_u(id, u_a, u_b, G, nz, diag, mask, alt_h)
   integer,                  intent(in) :: id   !< The ID for this diagnostic
diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90
index 7a67cbb5a5..3bd4dfc0de 100644
--- a/src/parameterizations/vertical/MOM_energetic_PBL.F90
+++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90
@@ -7,6 +7,7 @@ module MOM_energetic_PBL
 use MOM_coms,           only : EFP_type, real_to_EFP, EFP_to_real, operator(+), assignment(=), EFP_sum_across_PEs
 use MOM_debugging,      only : hchksum
 use MOM_diag_mediator,  only : post_data, register_diag_field, safe_alloc_alloc
+use MOM_diag_mediator,  only : post_data_3d_by_column, post_data_3d_final
 use MOM_diag_mediator,  only : time_type, diag_ctrl
 use MOM_domains,        only : create_group_pass, do_group_pass, group_pass_type
 use MOM_error_handler,  only : MOM_error, FATAL, WARNING, MOM_mesg
@@ -261,6 +262,7 @@ module MOM_energetic_PBL
   type(EFP_type), dimension(2) :: sum_its_BBL !< The total number of iterations and columns worked on
 
   !>@{ Diagnostic IDs
+  integer :: id_Kd_ePBL_col_by_col = -1
   integer :: id_ML_depth = -1, id_hML_depth = -1, id_TKE_wind = -1, id_TKE_mixing = -1
   integer :: id_ustar_ePBL = -1, id_bflx_ePBL = -1
   integer :: id_TKE_MKE = -1, id_TKE_conv = -1, id_TKE_forcing = -1
@@ -691,6 +693,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV,
                          u_star, u_star_mean, mech_TKE, dt, MLD_io, Kd, mixvel, mixlen, GV, &
                          US, CS, eCD, Waves, G, i, j)
       endif
+      call post_data_3d_by_column(CS%id_Kd_ePBL_col_by_col, Kd, CS%diag, i, j)
 
       ! Add the diffusivity due to bottom boundary layer mixing, if there is energy to drive this mixing.
       if (BBL_mixing) then
@@ -826,6 +829,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV,
     do K=1,nz+1 ; do i=is,ie ; Kd_int(i,j,K) = Kd_2d(i,K) ; enddo ; enddo
 
   enddo ! j-loop
+  call post_data_3d_final(CS%id_Kd_ePBL_col_by_col, CS%diag)
 
   if (CS%debug .and. BBL_mixing) then
     call hchksum(visc%BBL_meanKE_loss, "ePBL visc%BBL_meanKE_loss", G%HI, &
@@ -4334,6 +4338,8 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS)
 
 
 !/ Checking output flags
+  CS%id_Kd_ePBL_col_by_col = register_diag_field('ocean_model', 'Kd_ePBL_col_by_col', diag%axesTi, Time, &
+      'ePBL diapycnal diffusivity at interfaces posted column by column', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
   CS%id_ML_depth = register_diag_field('ocean_model', 'ePBL_h_ML', diag%axesT1, &
       Time, 'Surface boundary layer depth', units='m', conversion=US%Z_to_m, &
       cmor_long_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')