+Add a better variant of add_LOTW_BBL_diffusivity

  Added a refactored version of add_LOTW_BBL_diffusivity that is more accurate
by avoiding a subtraction of two large numbers to determine the distance from
the surface and more reproducible by avoiding the use of the Fortran sum
function.  This new version of the code is mathematically equivalent to the
original, and it is selected by setting the new runtime parameter
LOTW_BBL_ANSWER_DATE to be greater than 20240630.  This probably could have been
done reusing the existing parameter SET_DIFF_ANSWER_DATE, but this would have
meant answer changes for some cases that set this to a high value.  By default
the original code is used, but the default should be changed later to follow the
value of DEFAULT_ANSWER_DATE.  By default, all answers are bitwise identical but
there is a new runtime parameter in some cases.

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)