+(*)Use tv%SpV in MOM_sponge code

  Use tv%SpV convert thicknesses to vertical distances in apply_sponge when it
is allocated to replace multiplication by GV%H_to_Z, thereby eliminating the
dependence on GV%RHo_0 in this modue in non-Boussinesq mode.  The new internal
array dz_to_h is used to reduce the code duplication as a result of these
changes.  All answers in Boussinesq test cases are bitwise identical, but
answers change in fully non-Boussinesq mode.  In semi-Boussinesq mode answers
are mathematically equivalent but change at roundoff unless RHO_0 is an integer
power of 2.

src/parameterizations/vertical/MOM_sponge.F90

@@ -347,10 +347,14 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
                    ! give 0 at the surface [nondim].
 
   real :: e(SZK_(GV)+1)  ! The interface heights [Z ~> m], usually negative.
+  real :: dz_to_h(SZK_(GV)+1)  ! Factors used to convert interface height movement
+                   ! to thickness fluxes [H Z-1 ~> nondim or kg m-3]
   real :: e0       ! The height of the free surface [Z ~> m].
   real :: e_str    ! A nondimensional amount by which the reference
                    ! profile must be stretched for the free surfaces
                    ! heights in the two profiles to agree [nondim].
+  real :: w_mean   ! The vertical displacement of water moving upward through an
+                   ! interface within 1 timestep [Z ~> m].
   real :: w        ! The thickness of water moving upward through an
                    ! interface within 1 timestep [H ~> m or kg m-2].
   real :: wm       ! wm is w if w is negative and 0 otherwise [H ~> m or kg m-2].
@@ -384,9 +388,15 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
                   "work properly with i-mean sponges and a bulk mixed layer.")
 
     do j=js,je ; do i=is,ie ; e_D(i,j,nz+1) = -G%bathyT(i,j) ; enddo ; enddo
-    do k=nz,1,-1 ; do j=js,je ; do i=is,ie
-      e_D(i,j,K) = e_D(i,j,K+1) + h(i,j,k)*GV%H_to_Z
-    enddo ; enddo ; enddo
+    if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then
+      do k=nz,1,-1 ; do j=js,je ; do i=is,ie
+        e_D(i,j,K) = e_D(i,j,K+1) + GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)
+      enddo ; enddo ; enddo
+    else
+      do k=nz,1,-1 ; do j=js,je ; do i=is,ie
+        e_D(i,j,K) = e_D(i,j,K+1) + h(i,j,k)*GV%H_to_Z
+      enddo ; enddo ; enddo
+    endif
     do j=js,je
       do i=is,ie
         dilate(i) = (G%bathyT(i,j) + G%Z_ref) / (e_D(i,j,1) + G%bathyT(i,j))
@@ -424,20 +434,39 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
       do K=2,nz+1 ; do i=is,ie
         h_above(i,K) = h_above(i,K-1) + max(h(i,j,k-1)-GV%Angstrom_H, 0.0)
       enddo ; enddo
-      do K=2,nz
-        ! w is positive for an upward (lightward) flux of mass, resulting
-        ! in the downward movement of an interface.
-        w = damp_1pdamp * eta_mean_anom(j,K) * GV%Z_to_H
-        do i=is,ie
+
+      ! In both blocks below, w is positive for an upward (lightward) flux of mass,
+      ! resulting in the downward movement of an interface.
+      if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then
+        do K=2,nz
+          w_mean = damp_1pdamp * eta_mean_anom(j,K)
+          do i=is,ie
+            w = w_mean * 2.0*GV%RZ_to_H / (tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k))
+            if (w > 0.0) then
+              w_int(i,j,K) = min(w, h_below(i,K))
+              eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,K)
+            else
+              w_int(i,j,K) = max(w, -h_above(i,K))
+              ea(i,j,k) = ea(i,j,k) - w_int(i,j,K)
+            endif
+          enddo
+        enddo
+      else
+        do K=2,nz
+          w = damp_1pdamp * eta_mean_anom(j,K) * GV%Z_to_H
           if (w > 0.0) then
-            w_int(i,j,K) = min(w, h_below(i,K))
-            eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,K)
+            do i=is,ie
+              w_int(i,j,K) = min(w, h_below(i,K))
+              eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,K)
+            enddo
           else
-            w_int(i,j,K) = max(w, -h_above(i,K))
-            ea(i,j,k) = ea(i,j,k) - w_int(i,j,K)
+            do i=is,ie
+              w_int(i,j,K) = max(w, -h_above(i,K))
+              ea(i,j,k) = ea(i,j,k) - w_int(i,j,K)
+            enddo
           endif
         enddo
-      enddo
+      endif
       do k=1,nz ; do i=is,ie
         ea_k = max(0.0, -w_int(i,j,K))
         eb_k = max(0.0, w_int(i,j,K+1))
@@ -462,9 +491,20 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
     damp = dt * CS%Iresttime_col(c)
 
     e(1) = 0.0 ; e0 = 0.0
-    do K=1,nz
-      e(K+1) = e(K) - h(i,j,k)*GV%H_to_Z
-    enddo
+    if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then
+      do K=1,nz
+        e(K+1) = e(K) - GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)
+      enddo
+      dz_to_h(1) = GV%RZ_to_H / tv%SpV_avg(i,j,1)
+      do K=2,nz
+        dz_to_h(K) = 2.0*GV%RZ_to_H / (tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k))
+      enddo
+    else
+      do K=1,nz
+        e(K+1) = e(K) - h(i,j,k)*GV%H_to_Z
+        dz_to_h(K) = GV%Z_to_H
+      enddo
+    endif
     e_str = e(nz+1) / CS%Ref_eta(nz+1,c)
 
     if ( CS%bulkmixedlayer ) then
@@ -481,7 +521,7 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
       wpb = 0.0; wb = 0.0
       do k=nz,nkmb+1,-1
         if (GV%Rlay(k) > Rcv_ml(i,j)) then
-          w = MIN((((e(K)-e0) - e_str*CS%Ref_eta(K,c)) * damp)*GV%Z_to_H, &
+          w = MIN((((e(K)-e0) - e_str*CS%Ref_eta(K,c)) * damp)*dz_to_h(K), &
                     ((wb + h(i,j,k)) - GV%Angstrom_H))
           wm = 0.5*(w-ABS(w))
           do m=1,CS%fldno
@@ -537,7 +577,7 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
       wpb = 0.0
       wb = 0.0
       do k=nz,1,-1
-        w = MIN((((e(K)-e0) - e_str*CS%Ref_eta(K,c)) * damp)*GV%Z_to_H, &
+        w = MIN((((e(K)-e0) - e_str*CS%Ref_eta(K,c)) * damp)*dz_to_h(K), &
                   ((wb + h(i,j,k)) - GV%Angstrom_H))
         wm = 0.5*(w - ABS(w))
         do m=1,CS%fldno