+(*)EPBL_BBL_EFFIC_BUG & DRAG_DIFFUSIVITY_ANSWER_DATE

  Previously, visc%BBL_meanKE_loss (which was called visc%TKE_BBL before it was
renamed) was missing a factor of the square root of the drag coefficient
compared with the net loss of kinetic energy.  This was compensated for by
multiplication by factors of the square root of the drag coefficient in
add_drag_diffusivity and add_LOTW_BBL_diffusivity, but when an equivalent
expression was added in the ePBL BBL code this was erroneously omitted.
Moreover, Alistair has had a comment questioning this added factor in
add_LOTW_BBL_diffusivity for a decade without adequate resolution.

  To correct this confusing situation, visc%BBL_meanKE_loss has been changed to
include the missing factor of the square root of the drag coefficient, while the
new variable visc%BBL_meanKE_loss_sqrtCd was added to allow for the previous
answers to be recovered when the new runtime parameter EPBL_BBL_EFFIC_BUG is set
to true and DRAG_DIFFUSIVITY_ANSWER_DATE is set below 20250302.  Because the
ePBL bottom boundary layer was only added to dev/gfdl a month ago and has not
yet been merged into main, we can be confident that it has only received very
limited use as yet, so the default for EPBL_BBL_EFFIC_BUG is false but this
default will change answers when EPBL_BBL_EFFIC > 0.  The default for
DRAG_DIFFUSIVITY_ANSWER_DATE is 20250101, which will preserve the previous
answers, but the default should later be taken from DEFAULT_ANSWER_DATE.  By
default, answers are unchanged in any cases that are more than a month old, but
answers can change by default in a few very recent experiments.  There are two
new runtime parameters in some MOM_parameter_doc files.

src/core/MOM_variables.F90

@@ -239,6 +239,9 @@ module MOM_variables
                    !! h points [H T-1 ~> m s-1 or kg m-2 s-1].
     BBL_meanKE_loss, & !< The viscous loss of mean kinetic energy in the bottom boundary layer
                    !! [H L2 T-3 ~> m3 s-3 or W m-2].
+    BBL_meanKE_loss_sqrtCd, & !< The viscous loss of mean kinetic energy in the bottom boundary layer
+                   !! divided by the square root of the drag coefficient [H L2 T-3 ~> m3 s-3 or W m-2].
+                   !! This is being set only to retain old answers, and should be phased out.
     taux_shelf, &  !< The zonal stresses on the ocean under shelves [R Z L T-2 ~> Pa].
     tauy_shelf     !< The meridional stresses on the ocean under shelves [R Z L T-2 ~> Pa].
   real, allocatable, dimension(:,:) :: tbl_thick_shelf_u

src/parameterizations/vertical/MOM_energetic_PBL.F90

@@ -198,6 +198,9 @@ module MOM_energetic_PBL
                              !! energetics that accounts for an exponential decay of TKE from a
                              !! near-bottom source and an assumed piecewise linear linear profile
                              !! of the buoyancy flux response to a change in a diffusivity.
+  logical :: BBL_effic_bug   !< If true, overestimate the efficiency of the non-tidal ePBL bottom boundary
+                             !! layer diffusivity by a factor of 1/sqrt(CDRAG), which is often a factor of
+                             !! about 18.3.
 
   !/ Options for documenting differences from parameter choices
   integer :: options_diff    !< If positive, this is a coded integer indicating a pair of
@@ -662,7 +665,11 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV,
       if (BBL_mixing) then
         if (CS%MLD_iteration_guess .and. (CS%BBL_depth(i,j) > 0.0)) BBLD_io = CS%BBL_depth(i,j)
         BBLD_in = BBLD_io
-        BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%BBL_meanKE_loss(i,j)
+        if (CS%BBL_effic_bug) then
+          BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%BBL_meanKE_loss_sqrtCd(i,j)
+        else
+          BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%BBL_meanKE_loss(i,j)
+        endif
         u_star_BBL = max(visc%ustar_BBL(i,j), CS%ustar_min*GV%Z_to_H)
 
         ! Add in tidal dissipation energy at the bottom, noting that fluxes%BBL_tidal_dis is
@@ -3752,6 +3759,10 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS)
                  "length scale by rotation in the bottom boundary layer.  Making this larger "//&
                  "decreases the bottom boundary layer diffusivity.", &
                  units="nondim", default=CS%Ekman_scale_coef, do_not_log=no_BBL)
+  call get_param(param_file, mdl, "EPBL_BBL_EFFIC_BUG", CS%BBL_effic_bug, &
+                 "If true, overestimate the efficiency of the non-tidal ePBL bottom boundary "//&
+                 "layer diffusivity by a factor of 1/sqrt(CDRAG), which is often a factor of "//&
+                 "about 18.3.", default=.false., do_not_log=(CS%ePBL_BBL_effic<=0.0))
 
   call get_param(param_file, mdl, "DECAY_ADJUSTED_BBL_TKE", CS%decay_adjusted_BBL_TKE, &
                  "If true, include an adjustment factor in the bottom boundary layer energetics "//&

src/parameterizations/vertical/MOM_set_diffusivity.F90

@@ -163,11 +163,17 @@ module MOM_set_diffusivity
                               !! 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.  Values above 20240630 use more accurate
-                              !! expressions for cases where USE_LOTW_BBL_DIFFUSIVITY is true.
+                              !! expressions for cases where USE_LOTW_BBL_DIFFUSIVITY is true.  Values
+                              !! above 20250301 use less confusing expressions to set the bottom-drag
+                              !! generated diffusivity when USE_LOTW_BBL_DIFFUSIVITY is false.
   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.
+  integer :: drag_diff_answer_date !< The vintage of the order of arithmetic in the drag diffusivity
+                              !! calculations.  Values above 20250301 use less confusing expressions
+                              !! to set the bottom-drag generated diffusivity when
+                              !! USE_LOTW_BBL_DIFFUSIVITY is false.
 
   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
@@ -1414,7 +1420,11 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE,
       ! If ustar_h = 0, this is land so this value doesn't matter.
       I2decay(i) = 0.5*CS%IMax_decay
     endif
-    TKE(i) = ((CS%BBL_effic * cdrag_sqrt) * exp(-I2decay(i)*h(i,j,nz)) ) * visc%BBL_meanKE_loss(i,j)
+    if (CS%drag_diff_answer_date <= 20250301) then
+      TKE(i) = ((CS%BBL_effic * cdrag_sqrt) * exp(-I2decay(i)*h(i,j,nz)) ) * visc%BBL_meanKE_loss_sqrtCd(i,j)
+    else
+      TKE(i) = (CS%BBL_effic * exp(-I2decay(i)*h(i,j,nz)) ) * visc%BBL_meanKE_loss(i,j)
+    endif
 
     if (associated(fluxes%BBL_tidal_dis)) &
       TKE(i) = TKE(i) + fluxes%BBL_tidal_dis(i,j) * GV%RZ_to_H * &
@@ -1644,8 +1654,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bo
 
     ! Energy input at the bottom [H Z2 T-3 ~> m3 s-3 or W m-2].
     ! (Note that visc%BBL_meanKE_loss is in [H L2 T-3 ~> m3 s-3 or W m-2], set in set_BBL_TKE().)
-    !### I am still unsure about sqrt(cdrag) in this expressions - AJA
-    BBL_meanKE_dis = cdrag_sqrt * visc%BBL_meanKE_loss(i,j)
+    BBL_meanKE_dis = visc%BBL_meanKE_loss(i,j)
     ! Add in tidal dissipation energy at the bottom [H Z2 T-3 ~> m3 s-3 or W m-2].
     ! Note that BBL_tidal_dis is in [R Z L2 T-3 ~> W m-2].
     if (associated(fluxes%BBL_tidal_dis)) &
@@ -1956,6 +1965,9 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC)
     if (allocated(visc%BBL_meanKE_loss)) then
       do j=js,je ; do i=is,ie ; visc%BBL_meanKE_loss(i,j) = 0.0 ; enddo ; enddo
     endif
+    if (allocated(visc%BBL_meanKE_loss_sqrtCd)) then
+      do j=js,je ; do i=is,ie ; visc%BBL_meanKE_loss_sqrtCd(i,j) = 0.0 ; enddo ; enddo
+    endif
     return
   endif
 
@@ -2079,7 +2091,13 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC)
                   (G%areaCu(I,j)*(ustar(I)*ustar(I)))) + &
                  ((G%areaCv(i,J-1)*(vstar(i,J-1)*vstar(i,J-1))) + &
                   (G%areaCv(i,J)*(vstar(i,J)*vstar(i,J)))) ) )
-      visc%BBL_meanKE_loss(i,j) = &
+      visc%BBL_meanKE_loss(i,j) = cdrag_sqrt * &
+                 ((((G%areaCu(I-1,j)*(ustar(I-1)*u2_bbl(I-1))) + &
+                    (G%areaCu(I,j) * (ustar(I)*u2_bbl(I)))) + &
+                   ((G%areaCv(i,J-1)*(vstar(i,J-1)*v2_bbl(i,J-1))) + &
+                    (G%areaCv(i,J) * (vstar(i,J)*v2_bbl(i,J)))) )*G%IareaT(i,j))
+      ! The following line could be omitted if SET_DIFF_ANSWER_DATE > 20250301 and EPBL_BBL_EFFIC_BUG is false.
+      visc%BBL_meanKE_loss_sqrtCd(i,j) = &
                  ((((G%areaCu(I-1,j)*(ustar(I-1)*u2_bbl(I-1))) + &
                     (G%areaCu(I,j) * (ustar(I)*u2_bbl(I)))) + &
                    ((G%areaCv(i,J-1)*(vstar(i,J-1)*v2_bbl(i,J-1))) + &
@@ -2289,7 +2307,9 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_
   call get_param(param_file, mdl, "SET_DIFF_ANSWER_DATE", CS%answer_date, &
                "The vintage of the order of arithmetic and expressions in the set diffusivity "//&
                "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.", &
+               "while higher values use updated and more robust forms of the same expressions.  "//&
+               "Values above 20250301 also use less confusing expressions to set the bottom-drag "//&
+               "generated diffusivity when USE_LOTW_BBL_DIFFUSIVITY is false.", &
                default=default_answer_date, do_not_log=.not.GV%Boussinesq)
   if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701)
 
@@ -2405,6 +2425,12 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_
                "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.
+  call get_param(param_file, mdl, "DRAG_DIFFUSIVITY_ANSWER_DATE", CS%drag_diff_answer_date, &
+               "The vintage of the order of arithmetic in the drag diffusivity calculations.  "//&
+               "Values above 20250301 use less confusing expressions to set the bottom-drag "//&
+               "generated diffusivity when USE_LOTW_BBL_DIFFUSIVITY is false. ", &
+               default=20250101, do_not_log=CS%use_LOTW_BBL_diffusivity.or.(CS%BBL_effic<=0.0))
+               !### 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)

src/parameterizations/vertical/MOM_set_viscosity.F90

@@ -3140,6 +3140,7 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS
     allocate(visc%kv_bbl_v(isd:ied,JsdB:JedB), source=0.0)
     allocate(visc%ustar_bbl(isd:ied,jsd:jed), source=0.0)
     allocate(visc%BBL_meanKE_loss(isd:ied,jsd:jed), source=0.0)
+    allocate(visc%BBL_meanKE_loss_sqrtCd(isd:ied,jsd:jed), source=0.0)
 
     CS%id_bbl_thick_u = register_diag_field('ocean_model', 'bbl_thick_u', &
        diag%axesCu1, Time, 'BBL thickness at u points', 'm', conversion=US%Z_to_m)
@@ -3215,6 +3216,7 @@ subroutine set_visc_end(visc, CS)
   if (associated(visc%Kv_shear_Bu)) deallocate(visc%Kv_shear_Bu)
   if (allocated(visc%ustar_bbl)) deallocate(visc%ustar_bbl)
   if (allocated(visc%BBL_meanKE_loss)) deallocate(visc%BBL_meanKE_loss)
+  if (allocated(visc%BBL_meanKE_loss_sqrtCd)) deallocate(visc%BBL_meanKE_loss_sqrtCd)
   if (allocated(visc%taux_shelf)) deallocate(visc%taux_shelf)
   if (allocated(visc%tauy_shelf)) deallocate(visc%tauy_shelf)
   if (allocated(visc%tbl_thick_shelf_u)) deallocate(visc%tbl_thick_shelf_u)