+(*)Add and use G%Coriolis2Bu

  Added the new element Coriolis2Bu to the ocean_grid_type and the
dyn_horgrid_type to hold the square of the Coriolis parameter, and use this
array in 10 routines (including btstep, set_dtbt, calculate_diagnostic_fields,
VarMix_init, propagate_int_tide, Calculate_kappa_shear, Calc_kappa_shear_vertex
and add_MLrad_diffusivity) that had been calculating and averaging the square of
the Coriolis parameter.  This could change some answers with FMAs enabled
because the compilers were previously free to split up some of the squares
when averaging the squared Coriolis parameter, but without FMAs all answers are
bitwise identical.  This commit does add a new element to two transparent
types.

src/core/MOM_barotropic.F90

@@ -1649,8 +1649,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
               gtot_W(i,j) * (Datu(I-1,j)*G%IdxCu(I-1,j))) + &
              (gtot_N(i,j) * (Datv(i,J)*G%IdyCv(i,J)) + &
               gtot_S(i,j) * (Datv(i,J-1)*G%IdyCv(i,J-1)))) + &
-            ((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
-             (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2)) * CS%BT_Coriolis_scale**2 )
+            ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
+             (G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1))) * CS%BT_Coriolis_scale**2 )
       H_eff_dx2 = max(H_min_dyn * ((G%IdxT(i,j))**2 + (G%IdyT(i,j))**2), &
                       G%IareaT(i,j) * &
                         ((Datu(I,j)*G%IdxCu(I,j) + Datu(I-1,j)*G%IdxCu(I-1,j)) + &
@@ -2906,8 +2906,8 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
     Idt_max2 = 0.5 * (1.0 + 2.0*CS%bebt) * (G%IareaT(i,j) * &
       ((gtot_E(i,j)*Datu(I,j)*G%IdxCu(I,j) + gtot_W(i,j)*Datu(I-1,j)*G%IdxCu(I-1,j)) + &
        (gtot_N(i,j)*Datv(i,J)*G%IdyCv(i,J) + gtot_S(i,j)*Datv(i,J-1)*G%IdyCv(i,J-1))) + &
-      ((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
-       (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2)) * CS%BT_Coriolis_scale**2 )
+      ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
+       (G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1))) * CS%BT_Coriolis_scale**2 )
     if (Idt_max2 * min_max_dt2 > 1.0) min_max_dt2 = 1.0 / Idt_max2
   enddo ; enddo
   dtbt_max = sqrt(min_max_dt2 / dgeo_de)

src/core/MOM_grid.F90

@@ -171,7 +171,8 @@ module MOM_grid
     Dblock_v, &   !< Topographic depths at v-points at which the flow is blocked [Z ~> m].
     Dopen_v       !< Topographic depths at v-points at which the flow is open at width dx_Cv [Z ~> m].
   real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
-    CoriolisBu    !< The Coriolis parameter at corner points [T-1 ~> s-1].
+    CoriolisBu, & !< The Coriolis parameter at corner points [T-1 ~> s-1].
+    Coriolis2Bu   !< The square of the Coriolis parameter at corner points [T-2 ~> s-2].
   real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: &
     df_dx, &      !< Derivative d/dx f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1].
     df_dy         !< Derivative d/dy f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1].
@@ -581,6 +582,7 @@ subroutine allocate_metrics(G)
 
   ALLOC_(G%bathyT(isd:ied, jsd:jed)) ; G%bathyT(:,:) = -G%Z_ref
   ALLOC_(G%CoriolisBu(IsdB:IedB, JsdB:JedB)) ; G%CoriolisBu(:,:) = 0.0
+  ALLOC_(G%Coriolis2Bu(IsdB:IedB, JsdB:JedB)) ; G%Coriolis2Bu(:,:) = 0.0
   ALLOC_(G%dF_dx(isd:ied, jsd:jed)) ; G%dF_dx(:,:) = 0.0
   ALLOC_(G%dF_dy(isd:ied, jsd:jed)) ; G%dF_dy(:,:) = 0.0
 
@@ -626,8 +628,8 @@ subroutine MOM_grid_end(G)
 
   DEALLOC_(G%dx_Cv) ; DEALLOC_(G%dy_Cu)
 
-  DEALLOC_(G%bathyT)  ; DEALLOC_(G%CoriolisBu)
-  DEALLOC_(G%dF_dx)  ; DEALLOC_(G%dF_dy)
+  DEALLOC_(G%bathyT)  ; DEALLOC_(G%CoriolisBu) ; DEALLOC_(G%Coriolis2Bu)
+  DEALLOC_(G%dF_dx)   ; DEALLOC_(G%dF_dy)
   DEALLOC_(G%sin_rot) ; DEALLOC_(G%cos_rot)
 
   DEALLOC_(G%porous_DminU) ; DEALLOC_(G%porous_DmaxU) ; DEALLOC_(G%porous_DavgU)

src/core/MOM_transcribe_grid.F90

@@ -105,6 +105,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US)
     oG%dyBu(I,J) = dG%dyBu(I+ido,J+jdo)
     oG%areaBu(I,J) = dG%areaBu(I+ido,J+jdo)
     oG%CoriolisBu(I,J) = dG%CoriolisBu(I+ido,J+jdo)
+    oG%Coriolis2Bu(I,J) = dG%Coriolis2Bu(I+ido,J+jdo)
     oG%mask2dBu(I,J) = dG%mask2dBu(I+ido,J+jdo)
   enddo ; enddo
 
@@ -165,6 +166,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US)
   call pass_var(oG%geoLatBu, oG%Domain, position=CORNER)
   call pass_vector(oG%dxBu, oG%dyBu, oG%Domain, To_All+Scalar_Pair, BGRID_NE)
   call pass_var(oG%CoriolisBu, oG%Domain, position=CORNER)
+  call pass_var(oG%Coriolis2Bu, oG%Domain, position=CORNER)
   call pass_var(oG%mask2dBu, oG%Domain, position=CORNER)
 
   if (oG%bathymetry_at_vel) then
@@ -263,6 +265,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US)
     dG%dyBu(I,J) = oG%dyBu(I+ido,J+jdo)
     dG%areaBu(I,J) = oG%areaBu(I+ido,J+jdo)
     dG%CoriolisBu(I,J) = oG%CoriolisBu(I+ido,J+jdo)
+    dG%Coriolis2Bu(I,J) = oG%Coriolis2Bu(I+ido,J+jdo)
     dG%mask2dBu(I,J) = oG%mask2dBu(I+ido,J+jdo)
   enddo ; enddo
 
@@ -324,6 +327,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US)
   call pass_var(dG%geoLatBu, dG%Domain, position=CORNER)
   call pass_vector(dG%dxBu, dG%dyBu, dG%Domain, To_All+Scalar_Pair, BGRID_NE)
   call pass_var(dG%CoriolisBu, dG%Domain, position=CORNER)
+  call pass_var(dG%Coriolis2Bu, dG%Domain, position=CORNER)
   call pass_var(dG%mask2dBu, dG%Domain, position=CORNER)
 
   if (dG%bathymetry_at_vel) then

src/diagnostics/MOM_diagnostics.F90

@@ -679,8 +679,8 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
       do j=js,je ; do i=is,ie
         ! Blend the equatorial deformation radius with the standard one.
         f2_h = absurdly_small_freq2 + 0.25 * &
-            ((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
-             (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2))
+            ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
+             (G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1)))
         mag_beta = sqrt(0.5 * ( &
             (((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2 + &
              ((G%CoriolisBu(I,J-1)-G%CoriolisBu(I-1,J-1)) * G%IdxCv(i,J-1))**2) + &
@@ -729,8 +729,8 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
       do j=js,je ; do i=is,ie
         ! Blend the equatorial deformation radius with the standard one.
         f2_h = absurdly_small_freq2 + 0.25 * &
-            ((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
-             (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2))
+            ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
+             (G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1)))
         mag_beta = sqrt(0.5 * ( &
             (((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2 + &
              ((G%CoriolisBu(I,J-1)-G%CoriolisBu(I-1,J-1)) * G%IdxCv(i,J-1))**2) + &

src/framework/MOM_dyn_horgrid.F90

@@ -169,7 +169,8 @@ module MOM_dyn_horgrid
     Dblock_v, &   !< Topographic depths at v-points at which the flow is blocked [Z ~> m].
     Dopen_v       !< Topographic depths at v-points at which the flow is open at width dx_Cv [Z ~> m].
   real, allocatable, dimension(:,:) :: &
-    CoriolisBu    !< The Coriolis parameter at corner points [T-1 ~> s-1].
+    CoriolisBu, & !< The Coriolis parameter at corner points [T-1 ~> s-1].
+    Coriolis2Bu   !< The square of the Coriolis parameter at corner points [T-2 ~> s-2].
   real, allocatable, dimension(:,:) :: &
     df_dx, &      !< Derivative d/dx f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1].
     df_dy         !< Derivative d/dy f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1].
@@ -289,6 +290,7 @@ subroutine create_dyn_horgrid(G, HI, bathymetry_at_vel)
 
   allocate(G%bathyT(isd:ied, jsd:jed), source=0.0)
   allocate(G%CoriolisBu(IsdB:IedB, JsdB:JedB), source=0.0)
+  allocate(G%Coriolis2Bu(IsdB:IedB, JsdB:JedB), source=0.0)
   allocate(G%dF_dx(isd:ied, jsd:jed), source=0.0)
   allocate(G%dF_dy(isd:ied, jsd:jed), source=0.0)
 
@@ -360,6 +362,7 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns)
   call rotate_array_pair(G_in%dxBu, G_in%dyBu, turns, G%dxBu, G%dyBu)
   call rotate_array(G_in%areaBu, turns, G%areaBu)
   call rotate_array(G_in%CoriolisBu, turns, G%CoriolisBu)
+  call rotate_array(G_in%Coriolis2Bu, turns, G%Coriolis2Bu)
   call rotate_array(G_in%mask2dBu, turns, G%mask2dBu)
 
   ! Topography at the cell faces
@@ -528,8 +531,8 @@ subroutine destroy_dyn_horgrid(G)
   deallocate(G%porous_DminU) ; deallocate(G%porous_DmaxU) ; deallocate(G%porous_DavgU)
   deallocate(G%porous_DminV) ; deallocate(G%porous_DmaxV) ; deallocate(G%porous_DavgV)
 
-  deallocate(G%bathyT)  ; deallocate(G%CoriolisBu)
-  deallocate(G%dF_dx)  ; deallocate(G%dF_dy)
+  deallocate(G%bathyT)  ; deallocate(G%CoriolisBu) ; deallocate(G%Coriolis2Bu)
+  deallocate(G%dF_dx)   ; deallocate(G%dF_dy)
   deallocate(G%sin_rot) ; deallocate(G%cos_rot)
 
   if (allocated(G%Dblock_u)) deallocate(G%Dblock_u)

src/initialization/MOM_fixed_initialization.F90

@@ -60,14 +60,15 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir)
   logical,                 intent(in)    :: write_geom !< If true, write grid geometry files.
   character(len=*),        intent(in)    :: output_dir !< The directory into which to write files.
 
-  ! Local
+  ! Local variables
   character(len=200) :: inputdir   ! The directory where NetCDF input files are.
   character(len=200) :: config
   logical            :: read_porous_file
   character(len=40)  :: mdl = "MOM_fixed_initialization" ! This module's name.
+  integer :: I, J
   logical :: debug
-! This include declares and sets the variable "version".
-#include "version_variable.h"
+  ! This include declares and sets the variable "version".
+# include "version_variable.h"
 
   call callTree_enter("MOM_initialize_fixed(), MOM_fixed_initialization.F90")
   call log_version(PF, mdl, version, "")
@@ -156,8 +157,14 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir)
   call MOM_initialize_rotation(G%CoriolisBu, G, PF, US=US)
 !   Calculate the components of grad f (beta)
   call MOM_calculate_grad_Coriolis(G%dF_dx, G%dF_dy, G, US=US)
+!   Calculate the square of the Coriolis parameter
+  do I=G%IsdB,G%IedB ; do J=G%JsdB,G%JedB
+    G%Coriolis2Bu(I,J) = G%CoriolisBu(I,J)**2
+  enddo ; enddo
+
   if (debug) then
     call qchksum(G%CoriolisBu, "MOM_initialize_fixed: f ", G%HI, scale=US%s_to_T)
+    call qchksum(G%Coriolis2Bu, "MOM_initialize_fixed: f2 ", G%HI, scale=US%s_to_T**2)
     call hchksum(G%dF_dx, "MOM_initialize_fixed: dF_dx ", G%HI, scale=US%m_to_L*US%s_to_T)
     call hchksum(G%dF_dy, "MOM_initialize_fixed: dF_dy ", G%HI, scale=US%m_to_L*US%s_to_T)
   endif

src/parameterizations/lateral/MOM_internal_tides.F90

@@ -361,8 +361,8 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
   if (CS%energized_angle <= 0) then
     frac_per_sector = 1.0 / real(CS%nAngle)
     do m=1,CS%nMode ; do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie
-      f2 = 0.25*((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
-                 (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2))
+      f2 = 0.25*((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
+                 (G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1)))
       if (CS%frequency(fr)**2 > f2) &
         CS%En(i,j,a,fr,m) = CS%En(i,j,a,fr,m) + dt*frac_per_sector*(1.0-CS%q_itides) * &
                             CS%fraction_tidal_input(fr,m) * TKE_itidal_input(i,j,fr)
@@ -371,8 +371,8 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
     frac_per_sector = 1.0
     a = CS%energized_angle
     do m=1,CS%nMode ; do fr=1,CS%nFreq ; do j=js,je ; do i=is,ie
-      f2 = 0.25*((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
-                 (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2))
+      f2 = 0.25*((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
+                 (G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1)))
       if (CS%frequency(fr)**2 > f2) &
         CS%En(i,j,a,fr,m) = CS%En(i,j,a,fr,m) + dt*frac_per_sector*(1.0-CS%q_itides) * &
                             CS%fraction_tidal_input(fr,m) * TKE_itidal_input(i,j,fr)
@@ -630,8 +630,8 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
       do j=js,je ; do i=is,ie
         id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging
         ! Calculate horizontal phase velocity magnitudes
-        f2 = 0.25*((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
-                   (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2))
+        f2 = 0.25*((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
+                   (G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1)))
         Kmag2 = (freq2 - f2) / (cn(i,j,m)**2 + cn_subRO**2)
         c_phase = 0.0
         if (Kmag2 > 0.0) then
@@ -1134,8 +1134,8 @@ subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang)
 
   ! Do the refraction.
     do i=is,ie
-      f2 = 0.25* ((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
-                 (G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2))
+      f2 = 0.25* ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
+                 (G%Coriolis2Bu(I,J-1) + G%Coriolis2Bu(I-1,J)))
       favg = 0.25*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + &
                    (G%CoriolisBu(I,J-1) + G%CoriolisBu(I-1,J)))
       df_dx = 0.5*((G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1)) - &
@@ -1355,7 +1355,7 @@ subroutine propagate(En, cn, freq, dt, G, US, CS, NAngle, residual_loss)
     ! Fix indexing here later
     speed(:,:) = 0.0
     do J=jsh-1,jeh ; do I=ish-1,ieh
-      f2 = G%CoriolisBu(I,J)**2
+      f2 = G%Coriolis2Bu(I,J)
       speed(I,J) = 0.25*((cn(i,j) + cn(i+1,j+1)) + (cn(i+1,j) + cn(i,j+1))) * &
                      sqrt(max(freq2 - f2, 0.0)) * Ifreq
     enddo ; enddo
@@ -1385,12 +1385,12 @@ subroutine propagate(En, cn, freq, dt, G, US, CS, NAngle, residual_loss)
     enddo
 
     do j=jsh,jeh ; do I=ish-1,ieh
-      f2 = 0.5 * (G%CoriolisBu(I,J)**2 + G%CoriolisBu(I,J-1)**2)
+      f2 = 0.5 * (G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I,J-1))
       speed_x(I,j) = 0.5*(cn(i,j) + cn(i+1,j)) * G%mask2dCu(I,j) * &
                      sqrt(max(freq2 - f2, 0.0)) * Ifreq
     enddo ; enddo
     do J=jsh-1,jeh ; do i=ish,ieh
-      f2 = 0.5 * (G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J)**2)
+      f2 = 0.5 * (G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J))
       speed_y(i,J) = 0.5*(cn(i,j) + cn(i,j+1)) * G%mask2dCv(i,J) * &
                      sqrt(max(freq2 - f2, 0.0)) * Ifreq
     enddo ; enddo

src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90

@@ -1516,7 +1516,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
 
     do J=js-1,Jeq ; do I=is-1,Ieq
       CS%f2_dx2_q(I,J) = (G%dxBu(I,J)**2 + G%dyBu(I,J)**2) * &
-                         max(G%CoriolisBu(I,J)**2, absurdly_small_freq**2)
+                         max(G%Coriolis2Bu(I,J), absurdly_small_freq**2)
       CS%beta_dx2_q(I,J) = oneOrTwo * ((G%dxBu(I,J))**2 + (G%dyBu(I,J))**2) * (sqrt(0.5 * &
           ( (((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2 + &
              ((G%CoriolisBu(I+1,J)-G%CoriolisBu(I,J)) * G%IdxCv(i+1,J))**2) + &
@@ -1526,7 +1526,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
 
     do j=js,je ; do I=is-1,Ieq
       CS%f2_dx2_u(I,j) = (G%dxCu(I,j)**2 + G%dyCu(I,j)**2) * &
-          max(0.5* (G%CoriolisBu(I,J)**2+G%CoriolisBu(I,J-1)**2), absurdly_small_freq**2)
+          max(0.5* (G%Coriolis2Bu(I,J)+G%Coriolis2Bu(I,J-1)), absurdly_small_freq**2)
       CS%beta_dx2_u(I,j) = oneOrTwo * ((G%dxCu(I,j))**2 + (G%dyCu(I,j))**2) * (sqrt( &
           0.25*( (((G%CoriolisBu(I,J-1)-G%CoriolisBu(I-1,J-1)) * G%IdxCv(i,J-1))**2 + &
                   ((G%CoriolisBu(I+1,J)-G%CoriolisBu(I,J)) * G%IdxCv(i+1,J))**2) + &
@@ -1537,7 +1537,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
 
     do J=js-1,Jeq ; do i=is,ie
       CS%f2_dx2_v(i,J) = ((G%dxCv(i,J))**2 + (G%dyCv(i,J))**2) * &
-          max(0.5*(G%CoriolisBu(I,J)**2+G%CoriolisBu(I-1,J)**2), absurdly_small_freq**2)
+          max(0.5*(G%Coriolis2Bu(I,J)+G%Coriolis2Bu(I-1,J)), absurdly_small_freq**2)
       CS%beta_dx2_v(i,J) = oneOrTwo * ((G%dxCv(i,J))**2 + (G%dyCv(i,J))**2) * (sqrt( &
           ((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2 + &
           0.25*( (((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) * G%IdyCu(I,j))**2 + &
@@ -1572,8 +1572,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
     allocate(CS%f2_dx2_h(isd:ied,jsd:jed), source=0.0)
     do j=js-1,je+1 ; do i=is-1,ie+1
       CS%f2_dx2_h(i,j) = (G%dxT(i,j)**2 + G%dyT(i,j)**2) * &
-          max(0.25 * ((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
-                      (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2)), &
+          max(0.25 * ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
+                      (G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1))), &
               absurdly_small_freq**2)
       CS%beta_dx2_h(i,j) = oneOrTwo * ((G%dxT(i,j))**2 + (G%dyT(i,j))**2) * (sqrt(0.5 * &
           ( (((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2 + &

src/parameterizations/vertical/MOM_kappa_shear.F90

@@ -279,8 +279,8 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, &
         do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ; enddo
       endif
 
-      f2 = 0.25 * ((G%CoriolisBu(I,j)**2 + G%CoriolisBu(I-1,J-1)**2) + &
-                   (G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2))
+      f2 = 0.25 * ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
+                   (G%Coriolis2Bu(I,J-1) + G%Coriolis2Bu(I-1,J)))
       surface_pres = 0.0 ; if (associated(p_surf)) surface_pres = p_surf(i,j)
 
     ! ----------------------------------------------------    I_Ld2_1d, dz_Int_1d
@@ -551,7 +551,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_
         do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ; enddo
       endif
 
-      f2 = G%CoriolisBu(I,J)**2
+      f2 = G%Coriolis2Bu(I,J)
       surface_pres = 0.0
       if (associated(p_surf)) then
         if (CS%psurf_bug) then

src/parameterizations/vertical/MOM_set_diffusivity.F90

@@ -1631,8 +1631,8 @@ subroutine add_MLrad_diffusivity(dz, fluxes, tv, j, Kd_int, G, GV, US, CS, TKE_t
     if (CS%ML_omega_frac >= 1.0) then
       f_sq = 4.0 * Omega2
     else
-      f_sq = 0.25 * ((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
-                     (G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2))
+      f_sq = 0.25 * ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
+                     (G%Coriolis2Bu(I,J-1) + G%Coriolis2Bu(I-1,J)))
       if (CS%ML_omega_frac > 0.0) &
         f_sq = CS%ML_omega_frac * 4.0 * Omega2 + (1.0 - CS%ML_omega_frac) * f_sq
     endif