(*)Use summed thicknesses in OBC wave speed

  Changed the calculation of the external mode gravity wave speed used by the
radiation open boundary conditions and in the Kelvin-wave test case to use the
summed layer thicknesses rather than bathyT.  This change will prevent taking
the square root of negative thicknesses, even if there is wetting and drying so
some of the bathymetry is above the mean sea level, and it is more physically
accurate.  This PR will at least partly address MOM6 issue mom-ocean#1447.  Some
solutions will exhibit changing answers, including the barotropic Kelvin wave
test case, but after a broad discussion it was decided not to reproduce the
previous solutions with a runtime parameter.

src/core/MOM_open_boundary.F90

@@ -3811,8 +3811,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
           segment%h(I,j,k) = h(i+ishift,j,k)
           segment%Htot(I,j) = segment%Htot(I,j) + segment%h(I,j,k)
         enddo
-        segment%Cg(I,j) = sqrt(GV%g_prime(1)*G%bathyT(i+ishift,j))
-        !### This should be: segment%Cg(I,j) = sqrt(GV%g_prime(1)*segment%Htot(I,j)*GV%H_to_Z)
+        segment%Cg(I,j) = sqrt(GV%g_prime(1)*segment%Htot(I,j)*GV%H_to_Z)
       enddo
     else! (segment%direction == OBC_DIRECTION_N .or. segment%direction == OBC_DIRECTION_S)
       allocate(normal_trans_bt(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB))
@@ -3825,8 +3824,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
           segment%h(i,J,k) = h(i,j+jshift,k)
           segment%Htot(i,J) = segment%Htot(i,J) + segment%h(i,J,k)
         enddo
-        segment%Cg(i,J) = sqrt(GV%g_prime(1)*G%bathyT(i,j+jshift))
-        !### This should be: segment%Cg(i,J) = sqrt(GV%g_prime(1)*segment%Htot(i,J)*GV%H_to_Z)
+        segment%Cg(i,J) = sqrt(GV%g_prime(1)*segment%Htot(i,J)*GV%H_to_Z)
       enddo
     endif
 

src/user/Kelvin_initialization.F90

@@ -211,15 +211,11 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
   km_to_L_scale = 1000.0*US%m_to_L
 
   do j=jsd,jed ; do i=isd,ied
-    depth_tot(i,j) = G%bathyT(i,j)
+    depth_tot(i,j) = 0.0
   enddo ; enddo
-  !### Instead this should be:
-  ! do j=jsd,jed ; do i=isd,ied
-  !   depth_tot(i,j) = 0.0
-  ! enddo ; enddo
-  ! do j=jsd,jed ; do i=isd,ied
-  !   depth_tot(i,j) = depth_tot(i,j) + GV%H_to_Z * h(i,j,k)
-  ! enddo ; enddo
+  do k=1,nz ; do j=jsd,jed ; do i=isd,ied
+    depth_tot(i,j) = depth_tot(i,j) + GV%H_to_Z * h(i,j,k)
+  enddo ; enddo ; enddo
 
   if (CS%mode == 0) then
     mag_SSH = 1.0*US%m_to_Z