(*)Parenthesize MOM_set_diffusivity for FMAs

  Added parentheses to 4 expressions in add_drag_diffusivity, set_BBL_TKE and
add_LOTW_BBL_diffusivity setting the bottom-drag contributions to TKE and
friction velocity so that they will exhibit rotationally consistent solutions
when fused-multiply-adds are enabled.  All answers are bitwise identical in
cases without FMAs, but answers could change when FMAs are enabled.

src/parameterizations/vertical/MOM_set_diffusivity.F90

@@ -1324,10 +1324,10 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE,
 
       ! TKE_Ray has been initialized to 0 above.
       if (Rayleigh_drag) TKE_Ray = 0.5*CS%BBL_effic * US%L_to_Z**2 * G%IareaT(i,j) * &
-            ((G%areaCu(I-1,j) * visc%Ray_u(I-1,j,k) * u(I-1,j,k)**2 + &
-              G%areaCu(I,j)   * visc%Ray_u(I,j,k)   * u(I,j,k)**2) + &
-             (G%areaCv(i,J-1) * visc%Ray_v(i,J-1,k) * v(i,J-1,k)**2 + &
-              G%areaCv(i,J)   * visc%Ray_v(i,J,k)   * v(i,J,k)**2))
+            (((G%areaCu(I-1,j) * visc%Ray_u(I-1,j,k) * u(I-1,j,k)**2) + &
+              (G%areaCu(I,j)   * visc%Ray_u(I,j,k)   * u(I,j,k)**2)) + &
+             ((G%areaCv(i,J-1) * visc%Ray_v(i,J-1,k) * v(i,J-1,k)**2) + &
+              (G%areaCv(i,J)   * visc%Ray_v(i,J,k)   * v(i,J,k)**2)))
 
       if (TKE_to_layer + TKE_Ray > 0.0) then
         if (CS%BBL_mixing_as_max) then
@@ -1514,10 +1514,10 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bo
       ! Add in additional energy input from bottom-drag against slopes (sides)
       if (Rayleigh_drag) TKE_remaining = TKE_remaining + &
             0.5*CS%BBL_effic * US%L_to_Z**2 * G%IareaT(i,j) * &
-            ((G%areaCu(I-1,j) * visc%Ray_u(I-1,j,k) * u(I-1,j,k)**2 + &
-              G%areaCu(I,j)   * visc%Ray_u(I,j,k)   * u(I,j,k)**2) + &
-             (G%areaCv(i,J-1) * visc%Ray_v(i,J-1,k) * v(i,J-1,k)**2 + &
-              G%areaCv(i,J)   * visc%Ray_v(i,J,k)   * v(i,J,k)**2))
+            (((G%areaCu(I-1,j) * visc%Ray_u(I-1,j,k) * u(I-1,j,k)**2) + &
+              (G%areaCu(I,j)   * visc%Ray_u(I,j,k)   * u(I,j,k)**2)) + &
+             ((G%areaCv(i,J-1) * visc%Ray_v(i,J-1,k) * v(i,J-1,k)**2) + &
+              (G%areaCv(i,J)   * visc%Ray_v(i,J,k)   * v(i,J,k)**2)))
 
       ! Exponentially decay TKE across the thickness of the layer.
       ! This is energy loss in addition to work done as mixing, apparently to Joule heating.
@@ -1910,15 +1910,15 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC)
 
     do i=is,ie
       visc%ustar_BBL(i,j) = sqrt(0.5*G%IareaT(i,j) * &
-                ((G%areaCu(I-1,j)*(ustar(I-1)*ustar(I-1)) + &
-                  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))) ) )
+                (((G%areaCu(I-1,j)*(ustar(I-1)*ustar(I-1))) + &
+                  (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%TKE_BBL(i,j) = US%L_to_Z**2 * &
-                 (((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))
+                 ((((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))
     enddo
   enddo
   !$OMP end parallel