diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90
index ef2e4ed5f6..d78b675c1a 100644
--- a/src/parameterizations/vertical/MOM_set_diffusivity.F90
+++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90
@@ -157,7 +157,12 @@ module MOM_set_diffusivity
   integer :: answer_date      !< The vintage of the order of arithmetic and expressions in this module's
                               !! calculations.  Values below 20190101 recover the answers from the
                               !! end of 2018, while higher values use updated and more robust forms
-                              !! of the same expressions.
+                              !! of the same expressions.  Values above 20240630 use more accurate
+                              !! expressions for cases where USE_LOTW_BBL_DIFFUSIVITY is true.
+  integer :: LOTW_BBL_answer_date !< The vintage of the order of arithmetic and expressions
+                              !! in the LOTW_BBL calculations.  Values below 20240630 recover the
+                              !! original answers, while higher values use more accurate expressions.
+                              !! This only applies when USE_LOTW_BBL_DIFFUSIVITY is true.
 
   character(len=200) :: inputdir !< The directory in which input files are found
   type(user_change_diff_CS), pointer :: user_change_diff_CSp => NULL() !< Control structure for a child module
@@ -1426,6 +1431,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bo
 
   ! Local variables
   real :: dz(SZI_(G),SZK_(GV)) ! Height change across layers [Z ~> m]
+  real :: dz_above(SZK_(GV)+1) ! Distance from each interface to the surface [Z ~> m]
   real :: TKE_column       ! net TKE input into the column [H Z2 T-3 ~> m3 s-3 or W m-2]
   real :: TKE_remaining    ! remaining TKE available for mixing in this layer and above [H Z2 T-3 ~> m3 s-3 or W m-2]
   real :: TKE_consumed     ! TKE used for mixing in this layer [H Z2 T-3 ~> m3 s-3 or W m-2]
@@ -1500,7 +1506,15 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bo
     TKE_column = CS%BBL_effic * TKE_column ! Only use a fraction of the mechanical dissipation for mixing.
 
     TKE_remaining = TKE_column
-    total_depth = ( sum(dz(i,:)) + GV%dz_subroundoff ) ! Total column thickness [Z ~> m].
+    if (CS%LOTW_BBL_answer_date > 20240630) then
+      dz_above(1) = GV%dz_subroundoff  ! This could perhaps be 0 instead.
+      do K=2,GV%ke+1
+        dz_above(K) = dz_above(K-1) + dz(i,k-1)
+      enddo
+      total_depth = dz_above(GV%ke+1)
+    else
+      total_depth = ( sum(dz(i,:)) + GV%dz_subroundoff ) ! Total column thickness [Z ~> m].
+    endif
     ustar_D = ustar * total_depth
     h_bot = 0.
     z_bot = 0.
@@ -1525,7 +1539,11 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bo
 
       z_bot = z_bot + dz(i,k)  ! Distance between upper interface of layer and the bottom [Z ~> m].
       h_bot = h_bot + h(i,j,k) ! Thickness between upper interface of layer and the bottom [H ~> m or kg m-2].
-      D_minus_z = max(total_depth - z_bot, 0.) ! Thickness above layer [H ~> m or kg m-2].
+      if (CS%LOTW_BBL_answer_date > 20240630) then
+        D_minus_z = dz_above(K)
+      else
+        D_minus_z = max(total_depth - z_bot, 0.) ! Distance from the interface to the surface [Z ~> m].
+      endif
 
       ! Diffusivity using law of the wall, limited by rotation, at height z [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
       ! This calculation is at the upper interface of the layer
@@ -2191,7 +2209,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_
                  "law of the form c_drag*|u|*u. The velocity magnitude "//&
                  "may be an assumed value or it may be based on the actual "//&
                  "velocity in the bottommost HBBL, depending on LINEAR_DRAG.", default=.true.)
-  if  (CS%bottomdraglaw) then
+  if (CS%bottomdraglaw) then
     call get_param(param_file, mdl, "CDRAG", CS%cdrag, &
                  "The drag coefficient relating the magnitude of the "//&
                  "velocity field to the bottom stress. CDRAG is only used "//&
@@ -2231,6 +2249,14 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_
   else
     CS%use_LOTW_BBL_diffusivity = .false. ! This parameterization depends on a u* from viscous BBL
   endif
+  call get_param(param_file, mdl, "LOTW_BBL_ANSWER_DATE", CS%LOTW_BBL_answer_date, &
+               "The vintage of the order of arithmetic and expressions in the LOTW_BBL "//&
+               "calculations.  Values below 20240630 recover the original answers, while "//&
+               "higher values use more accurate expressions.  This only applies when "//&
+               "USE_LOTW_BBL_DIFFUSIVITY is true.", &
+               default=20190101, do_not_log=.not.CS%use_LOTW_BBL_diffusivity)
+               !### Set default as default=default_answer_date, or use SET_DIFF_ANSWER_DATE.
+
   CS%id_Kd_BBL = register_diag_field('ocean_model', 'Kd_BBL', diag%axesTi, Time, &
                  'Bottom Boundary Layer Diffusivity', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)