(*)Set dSV_dT and dSV_dS with unassociated fluxes

  Set dSV_dT and dSV_dS if present in applyBoundaryFluxesInOut, even if boundary
fluxes are not associated.  With this change, setting BUOY_CONFIG='NONE' and
BUOY_CONFIG='zero' both work and give similar (but not identical) answers in
some test cases with an ePBL boundary layer parameterization, whereas before
answers were tainted with uninitialized values when BUOY_CONFIG='NONE'.  All
answers in the existing MOM6-examples test suite are bitwise identical, but
answers can change in other cases.

src/parameterizations/vertical/MOM_diabatic_aux.F90

@@ -822,19 +822,18 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
 
-  ! Only apply forcing if fluxes%sw is associated.
-  if (.not.associated(fluxes%sw)) return
-
-#define _OLD_ALG_
   Idt = 1.0 / dt
 
   calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS))
   calculate_buoyancy = present(SkinBuoyFlux)
   if (calculate_buoyancy) SkinBuoyFlux(:,:) = 0.0
+  if (present(cTKE)) cTKE(:,:,:) = 0.0
   g_Hconv2 = (US%L_to_Z**2*GV%g_Earth * GV%H_to_RZ) * GV%H_to_RZ
   EOSdom(:) = EOS_domain(G%HI)
 
-  if (present(cTKE)) cTKE(:,:,:) = 0.0
+  ! Only apply forcing if fluxes%sw is associated.
+  if (.not.associated(fluxes%sw) .and. .not.calculate_energetics) return
+
   if (calculate_buoyancy) then
     SurfPressure(:) = 0.0
     GoRho = US%L_to_Z**2*GV%g_Earth / GV%Rho0
@@ -874,7 +873,6 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
       h2d(i,k) = h(i,j,k)
       T2d(i,k) = tv%T(i,j,k)
     enddo ; enddo
-    if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=(1.0/GV%m_to_H))
 
     if (calculate_energetics) then
       ! The partial derivatives of specific volume with temperature and
@@ -898,6 +896,11 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
       pen_TKE_2d(:,:) = 0.0
     endif
 
+    ! Nothing more is done on this j-slice if there is no buoyancy forcing.
+    if (.not.associated(fluxes%sw)) cycle
+
+    if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=(1.0/GV%m_to_H))
+
     ! The surface forcing is contained in the fluxes type.
     ! We aggregate the thermodynamic forcing for a time step into the following:
     ! netMassInOut = surface water fluxes [H ~> m or kg m-2] over time step

src/parameterizations/vertical/MOM_vert_friction.F90

@@ -1143,12 +1143,12 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i,
   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 :: z2      ! A copy of z_i [nondim]
+  real :: botfn   ! A function that is 1 at the bottom and small far from it [nondim]
   real :: topfn   ! A function that is 1 at the top and small far from it [nondim]
   real :: kv_top  ! A viscosity associated with the top boundary layer [Z2 T-1 ~> m2 s-1]
   logical :: do_shelf, do_OBCs
   integer :: i, k, is, ie, max_nk
   integer :: nz
-  real    :: botfn
 
   a_cpl(:,:) = 0.0
   Kv_tot(:,:) = 0.0