Added comments highlighting bugs

  Added comments denoted with "###" indicating bugs or conceptual errors in
update_OBC_segment_data, horizontal_viscosity, thickness_diffuse and wave_speed.
Only comments and the case of some indices are changed in this commit, and all
answers are bitwise identical.  Actually correcting these bugs would probably
change answers in some cases.

src/core/MOM_open_boundary.F90

@@ -3802,29 +3802,31 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
     ishift=0;jshift=0
     if (segment%is_E_or_W) then
       allocate(normal_trans_bt(segment%HI%IsdB:segment%HI%IedB,segment%HI%jsd:segment%HI%jed))
-      normal_trans_bt(:,:)=0.0
+      normal_trans_bt(:,:) = 0.0
       if (segment%direction == OBC_DIRECTION_W) ishift=1
       I=segment%HI%IsdB
       do j=segment%HI%jsd,segment%HI%jed
-        segment%Cg(I,j) = sqrt(GV%g_prime(1)*G%bathyT(i+ishift,j))
-        segment%Htot(I,j)=0.0
+        segment%Htot(I,j) = 0.0
         do k=1,GV%ke
           segment%h(I,j,k) = h(i+ishift,j,k)
-          segment%Htot(I,j)=segment%Htot(I,j)+segment%h(I,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)
       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))
-      normal_trans_bt(:,:)=0.0
+      normal_trans_bt(:,:) = 0.0
       if (segment%direction == OBC_DIRECTION_S) jshift=1
       J=segment%HI%JsdB
       do i=segment%HI%isd,segment%HI%ied
-        segment%Cg(i,J) = sqrt(GV%g_prime(1)*G%bathyT(i,j+jshift))
-        segment%Htot(i,J)=0.0
+        segment%Htot(i,J) = 0.0
         do k=1,GV%ke
           segment%h(i,J,k) = h(i,j+jshift,k)
-          segment%Htot(i,J)=segment%Htot(i,J)+segment%h(i,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,j+jshift))
+        !### This should be: segment%Cg(i,J) = sqrt(GV%g_prime(1)*segment%Htot(i,J)*GV%H_to_Z)
       enddo
     endif
 
@@ -4715,7 +4717,7 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC)
   integer :: i, j
   integer :: l_seg
   logical :: fatal_error = .False.
-  real    :: min_depth
+  real    :: min_depth ! The minimum depth for ocean points [Z ~> m]
   integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2
   character(len=256) :: mesg    ! Message for error messages.
   type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list
@@ -4730,7 +4732,6 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC)
   allocate(color(G%isd:G%ied, G%jsd:G%jed)) ; color = 0
   allocate(color2(G%isd:G%ied, G%jsd:G%jed)) ; color2 = 0
 
-
   ! Paint a frame around the outside.
   do j=G%jsd,G%jed
     color(G%isd,j) = cedge

src/diagnostics/MOM_wave_speed.F90

@@ -444,7 +444,8 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_
               hw = 0.5*(Hc(k-1)+Hc(k))
               gp = gprime(K)
               if (l_mono_N2_column_fraction>0. .or. l_mono_N2_depth>=0.) then
-                if ( ((G%bathyT(i,j)-sum_hc < l_mono_N2_column_fraction*G%bathyT(i,j)) .or. &
+                !### Change to: if ( ((htot(i) - sum_hc < l_mono_N2_column_fraction*htot(i)) .or. & ) )
+                if ( ((G%bathyT(i,j) - sum_hc < l_mono_N2_column_fraction*G%bathyT(i,j)) .or. &
                       ((l_mono_N2_depth >= 0.) .and. (sum_hc > l_mono_N2_depth))) .and. &
                      (L2_to_Z2*gp > N2min*hw) ) then
                   ! Filters out regions where N2 increases with depth but only in a lower fraction

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -496,7 +496,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     enddo ; enddo
 
     do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
-      grad_vel_mag_bt_q(I,J) = boundary_mask_q(I,J) * (dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
+      grad_vel_mag_bt_q(I,J) = boundary_mask_q(I,J) * (dvdx_bt(I,J)**2 + dudy_bt(I,J)**2 + &
             (0.25*((dudx_bt(i,j)+dudx_bt(i+1,j+1))+(dudx_bt(i,j+1)+dudx_bt(i+1,j))))**2 + &
             (0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1))+(dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2)
     enddo ; enddo
@@ -1389,7 +1389,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
         if (grad_vel_mag_bt_h(i,j)>0) then
-          GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * &
+          GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j) / CS%GME_h0, 1.0)**2) * &
             (0.25*(KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)+KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k)))
         else
           GME_coeff = 0.0
@@ -1405,8 +1405,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       enddo ; enddo
 
       do J=js-1,Jeq ; do I=is-1,Ieq
-        if (grad_vel_mag_bt_q(i,j)>0) then
-          GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * &
+        if (grad_vel_mag_bt_q(I,J)>0) then
+          !### This expression is not rotationally invariant - bathyT is to the SW of q points,
+          !    and it needs parentheses in the sum of the 4 diffusivities.
+          GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j) / CS%GME_h0, 1.0)**2) * &
               (0.25*(KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)+KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k)))
         else
           GME_coeff = 0.0

src/parameterizations/lateral/MOM_thickness_diffuse.F90

@@ -481,14 +481,19 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
         do j=js,je ; do I=is-1,ie
           hu(I,j)       = 2.0*h(i,j,k)*h(i+1,j,k)/(h(i,j,k)+h(i+1,j,k)+h_neglect)
           if (hu(I,j) /= 0.0) hu(I,j) = 1.0
+          !### The same result would be accomplished with the following without a division:
+          ! hu(I,j) = 0.0 ; if (h(i,j,k)*h(i+1,j,k) /= 0.0) hu(I,j) = 1.0
           KH_u_lay(I,j) = 0.5*(KH_u(I,j,k)+KH_u(I,j,k+1))
         enddo ; enddo
         do J=js-1,je ; do i=is,ie
           hv(i,J)       = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect)
           if (hv(i,J) /= 0.0) hv(i,J) = 1.0
+          !### The same result would be accomplished with the following without a division:
+          ! hv(i,J) = 0.0 ; if (h(i,j,k)*h(i,j+1,k) /= 0.0) hv(i,J) = 1.0
           KH_v_lay(i,J) = 0.5*(KH_v(i,J,k)+KH_v(i,J,k+1))
         enddo ; enddo
         ! diagnose diffusivity at T-point
+        !### Because hu and hv are nondimensional here, the denominator is dimensionally inconsistent.
         do j=js,je ; do i=is,ie
           Kh_t(i,j,k) = ((hu(I-1,j)*KH_u_lay(i-1,j)+hu(I,j)*KH_u_lay(I,j))  &
                         +(hv(i,J-1)*KH_v_lay(i,J-1)+hv(i,J)*KH_v_lay(i,J))) &
@@ -505,7 +510,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
         enddo
 
         do j=js,je ; do i=is,ie
-          MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(1.0,G%bathyT(i,j))
+          MEKE%Kh_diff(i,j) = GV%H_to_Z * MEKE%Kh_diff(i,j) / MAX(1.0*US%m_to_Z, G%bathyT(i,j))
         enddo ; enddo
       endif