MEKE: No-division implementation of src_btm_drag

This patch modifies the calculation of `src_btm_drag` to avoid
additional division operations.

We first note the following variable renames:

* `Isfac` is now `damping`, since it describes the net damping or
  reduction of any damped fields.  `sfac` has been removed since it no
  longer appears in the expressions.

* `ldamping` is now `damp_rate`.  This is primarily to avoid confusion
  with `damping`, but also to more accurately describe its role.

* `ldamping_Strang1` is replaced to `damp_rate_s1`.  It is used for a
  similar purpose (cacheing of values from the first stage) but now
  holds a slightly different value.

The following modifications were made.

1. The conditional split of `sdt` into `sdt_damping` substeps is
   quantified by a new term, `damp_step` which is equal to 0.5 or 1.
   `sdt_damp` is computed as `sdt * damp_step`.

   The presence of `sdt_damp / sdt` in our expressions is replaced with
   `damp_step`, which avoids an extra division.

2. The first expression (using new notation) was originally

    sfac = 1 + sdt_damp * damp_rate
    D = M * (1 - sfac) / (sdt * sfac)

   where `D` is `src_btm_drag` and `M` is `MEKE_current(i,j)`

   This has been transformed to

    D = - M * (sdt_damp * damp_rate)
            / (sdt * (1 + sdt_damp * damp_rate))

      = - M * (sdt_damp / sft) * damp_rate
            * (1 / (1 + sdt_damp * damp_rate))

      = - M * damp_step * damp_rate * damping

   This new expression for `D` no longer requires a division.

3. In the second stage expression for `src_btm_drag`, we again use
   `damp_rate` to replace `ldamping`.  We also use `damp_rate1` to
   denote the original value of `ldaming_Strang1`.

    sfac = (1 + sdt_damp * damp_rate1) * (1 + sdt_damp * damp_rate)
    D = M * (1 - sfac) / (sdt * sfac)

   Using `damping1` to denote the damping from the first stage, `D` can
   be transformed as follows.

    D = -M * (sdt_damp * damp_rate + sdt_damp * damp_rate1
                + sdt_damp**2 * damp_rate * damp_rate1)
      / (sdt * (1 + sdt_damp * damp_rate) * (1 + sdt_damp * damp_rate1))

      = -M * (sdt_damp / sdt) * (1 / (1 + sdt_damp * damp_rate))
        * (damp_rate + damp_rate1 + sdt_damp * damp_rate * damp_rate1)
        / (1 + sdt_damp * damp_rate1)

      = -M * damp_step * damping
        * (damp_rate * (1 + sdt_damp * damp_rate1) + damp_rate1)
        / (1 + sdt_damp * damp_rate1)

      = -M * damp_step * damping * (damp_rate + damp_rate1 * damping1)

   We now store `damp_rate1 * damping1` in `damp_rate_s1(:,:)` rather
   than just `damping1`, and can reduce this expression to

    D = -M * damp_step * damping * (damp_rate + damp_rate_s1(i,j))

   which requires no division.  It also requires no additional
   read or writes.

Assuming there are no mistakes here, this should only cause negligible
bitwise differences in the results.

Also, comments have been added which explain that we cannot modify the
existing expressions for `MEKE%MEKE`, since they would alter the bitwise
results of existing runs.

src/parameterizations/lateral/MOM_MEKE.F90

@@ -201,7 +201,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
     src_GM, &       ! The MEKE source/tendency from the thickness mixing (GM) [L2 T-3 ~> W kg-1] (= m2 s-3).
     src_mom_lp, &   ! The MEKE source/tendency from the Laplacian of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3).
     src_mom_bh, &   ! The MEKE source/tendency from the biharmonic of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3).
-    ldamping_Strang1, &  ! The MEKE damping rate computed at the 1st Strang splitting stage [T-1 ~> s-1].
+    damp_rate_s1, & ! The MEKE damping rate computed at the 1st Strang splitting stage [T-1 ~> s-1].
     MEKE_current, & ! A copy of MEKE for use in computing the MEKE damping [L2 T-2 ~> m2 s-2].
     drag_rate_visc, & ! Near-bottom velocity contribution to bottom drag [H T-1 ~> m s-1 or kg m-2 s-1]
     drag_rate, &    ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1].
@@ -233,11 +233,11 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
   real :: cdrag2    ! The square of the drag coefficient times unit conversion factors [H2 L-2 ~> nondim or kg2 m-6]
   real :: advFac    ! The product of the advection scaling factor and 1/dt [T-1 ~> s-1]
   real :: mass_neglect ! A negligible mass [R Z ~> kg m-2].
-  real :: ldamping  ! The MEKE damping rate [T-1 ~> s-1].
   real :: sdt       ! dt to use locally [T ~> s] (could be scaled to accelerate)
   real :: sdt_damp  ! dt for damping [T ~> s] (sdt could be split).
-  real :: sfac      ! A factor needed to compute damping due to Strang splitting [nondim]
-  real :: Isfac     ! Inverse of sfac [nondim]
+  real :: damp_step ! Size of damping timestep relative to sdt [nondim]
+  real :: damp_rate ! The MEKE damping rate [T-1 ~> s-1].
+  real :: damping   ! The net damping of a field after sdt_damp [nondim]
   logical :: use_drag_rate ! Flag to indicate drag_rate is finite
   integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
   real(kind=real32), dimension(size(MEKE%MEKE),NUM_FEATURES) :: features_array ! The array of features
@@ -287,7 +287,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
 
     ! With a depth-dependent (and possibly strong) damping, it seems
     ! advisable to use Strang splitting between the damping and diffusion.
-    sdt_damp = sdt ; if (CS%MEKE_KH >= 0.0 .or. CS%MEKE_K4 >= 0.) sdt_damp = 0.5*sdt
+    damp_step = 1.
+    if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) damp_step = 0.5
+    sdt_damp = sdt * damp_step
 
     ! Calculate depth integrated mass exchange if doing advection [R Z L2 ~> kg]
     if (CS%MEKE_advection_factor>0.) then
@@ -479,19 +481,29 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
     ! First stage of Strang splitting
     !$OMP parallel do default(shared)
     do j=js,je ; do i=is,ie
-      ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
-      if (MEKE%MEKE(i,j) < 0.) ldamping = 0.
+      damp_rate = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
+      if (MEKE%MEKE(i,j) < 0.) damp_rate = 0.
       ! notice that the above line ensures a damping only if MEKE is positive,
       ! while leaving MEKE unchanged if it is negative
-      Isfac = 1. / (1. + sdt_damp * ldamping)
-      MEKE%MEKE(i,j) =  MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
-      MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
-      ldamping_Strang1(i,j) = ldamping
-      src_GM(i,j) = src_GM(i,j) * Isfac
-      src_mom_lp(i,j) = src_mom_lp(i,j) * Isfac
-      src_mom_bh(i,j) = src_mom_bh(i,j) * Isfac
-      sfac = ( 1.0 + sdt_damp*ldamping )
-      src_btm_drag(i,j) = MEKE_current(i,j) * ( (1.0 - sfac) / ( sdt * sfac ) )
+
+      damping = 1. / (1. + sdt_damp * damp_rate)
+
+      ! NOTE: MEKE%MEKE should use `damping` but we must preserve the existing
+      ! expression for bit reproducibility
+      MEKE%MEKE(i,j) =  MEKE%MEKE(i,j) / (1. + sdt_damp * damp_rate)
+      MEKE_decay(i,j) = damp_rate * G%mask2dT(i,j)
+
+      src_GM(i,j) = src_GM(i,j) * damping
+      src_mom_lp(i,j) = src_mom_lp(i,j) * damping
+      src_mom_bh(i,j) = src_mom_bh(i,j) * damping
+
+      src_btm_drag(i,j) = - MEKE_current(i,j) * ( &
+          damp_step * (damp_rate * damping) &
+      )
+
+      ! Store the effective damping rate if sdt is split
+      if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) &
+        damp_rate_s1(i,j) = damp_rate * damping
     enddo ; enddo
 
     if (CS%kh_flux_enabled .or. CS%MEKE_K4 >= 0.0) then
@@ -658,20 +670,27 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
         endif
         !$OMP parallel do default(shared)
         do j=js,je ; do i=is,ie
-          ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
-          if (MEKE%MEKE(i,j) < 0.) ldamping = 0.
+          damp_rate = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
+          if (MEKE%MEKE(i,j) < 0.) damp_rate = 0.
           ! notice that the above line ensures a damping only if MEKE is positive,
           ! while leaving MEKE unchanged if it is negative
-          Isfac = 1. / (1. + sdt_damp*ldamping)
-          MEKE%MEKE(i,j) =  MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
-          MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
-          src_GM(i,j) = src_GM(i,j) * Isfac
-          src_mom_lp(i,j) = src_mom_lp(i,j) * Isfac
-          src_mom_bh(i,j) = src_mom_bh(i,j) * Isfac
-          src_adv(i,j) = src_adv(i,j) * Isfac
-          src_mom_K4(i,j) = src_mom_K4(i,j) * Isfac
-          sfac = ( 1.0 + sdt_damp*ldamping_Strang1(i,j) ) * ( 1.0 + sdt_damp*ldamping )
-          src_btm_drag(i,j) = MEKE_current(i,j) * ( (1.0 - sfac) / ( sdt * sfac ) )
+
+          damping = 1. / (1. + sdt_damp * damp_rate)
+
+          ! NOTE: As above, MEKE%MEKE should use `damping` but we must preserve
+          !   the existing expression for bit reproducibility.
+          MEKE%MEKE(i,j) =  MEKE%MEKE(i,j) / (1.0 + sdt_damp*damp_rate)
+          MEKE_decay(i,j) = damp_rate*G%mask2dT(i,j)
+
+          src_GM(i,j) = src_GM(i,j) * damping
+          src_mom_lp(i,j) = src_mom_lp(i,j) * damping
+          src_mom_bh(i,j) = src_mom_bh(i,j) * damping
+          src_adv(i,j) = src_adv(i,j) * damping
+          src_mom_K4(i,j) = src_mom_K4(i,j) * damping
+
+          src_btm_drag(i,j) = -MEKE_current(i,j) * ( &
+             damp_step * damping * (damp_rate + damp_rate_s1(i,j)) &
+          )
         enddo ; enddo
       endif
     endif ! MEKE_KH>=0