(*)Fix northern Flather OBCs in southern halos

  Added additional logic to avoid segmentation faults from applying Flather open
boundary conditions near the opposite side of the computational domain (e.g.,
applying a Northern Flather boundary condition on the southern edge of the
domain).  This situation does not arise when open boundaries are only found on
the edges of a domain, but in parallel cases with sloped OBCs (such as the
ESMG-configs rotated seamount test case) this atypical situation can arise.  The
key point to the specific fix is to note that when such a case arises, the OBC
velocity in the outer halo point does not subsequently impact the solution on
that PE either because its effects are masked out by land or OBC points at the
neighboring orthogonal velocities, and then they are overwritten by a halo
update (perhaps several steps later).  All answers are bitwise identical in all
cases that have been tested and were working before, and the
ESMG/rotated_seamount test case is now running correctly (and giving identical
answers) for a variety of PE counts.

Author: Hallberg-NOAA

src/core/MOM_barotropic.F90

@@ -3686,7 +3686,11 @@ subroutine apply_u_velocity_OBCs(ubt, uhbt, ubt_trans, eta, SpV_avg, ubt_old, BT
     elseif (BT_OBC%u_OBC_type(I,j) == FLATHER_OBC) then  ! Eastern Flather OBC
       cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL
       u_inlet = cfl*ubt_old(I-1,j) + (1.0-cfl)*ubt_old(I,j)  ! Valid for cfl<1
-      if (GV%Boussinesq) then
+      if (I <= MS%isdw) then
+        ! Do not apply an Eastern Flather OBC at the western halo points on a PE, as doing so would
+        ! create a segmentation fault and this velocity will not propagate through to the next iteration.
+        ssh_in = BT_OBC%SSH_outer_u(I,j)
+      elseif (GV%Boussinesq) then
         ssh_in = GV%H_to_Z*(eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i-1,j)))      ! internal
       else
         ssh_1 = GV%H_to_RZ * eta(i,j) * SpV_avg(i,j) - (CS%bathyT(i,j) + G%Z_ref)
@@ -3739,7 +3743,11 @@ subroutine apply_u_velocity_OBCs(ubt, uhbt, ubt_trans, eta, SpV_avg, ubt_old, BT
     elseif (BT_OBC%u_OBC_type(I,j) == -FLATHER_OBC) then  ! Western Flather OBC
       cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL
       u_inlet = cfl*ubt_old(I+1,j) + (1.0-cfl)*ubt_old(I,j)  ! Valid for cfl<1
-      if (GV%Boussinesq) then
+      if (I >= MS%iedw-1) then
+        ! Do not apply a Western Flather OBC at the eastern halo points on a PE, as doing so would
+        ! create a segmentation fault and this velocity will not propagate through to the next iteration.
+        ssh_in = BT_OBC%SSH_outer_u(I,j)
+      elseif (GV%Boussinesq) then
         ssh_in = GV%H_to_Z*(eta(i+1,j) + (0.5-cfl)*(eta(i+1,j)-eta(i+2,j)))  ! internal
       else
         ssh_1 = GV%H_to_RZ * eta(i+1,j) * SpV_avg(i+1,j) - (CS%bathyT(i+1,j) + G%Z_ref)
@@ -3872,7 +3880,11 @@ subroutine apply_v_velocity_OBCs(vbt, vhbt, vbt_trans, eta, SpV_avg, vbt_old, BT
     elseif (BT_OBC%v_OBC_type(i,J) == FLATHER_OBC) then  ! Northern Flather OBC
       cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(i,J) ! CFL
       v_inlet = cfl*vbt_old(i,J-1) + (1.0-cfl)*vbt_old(i,J)  ! Valid for cfl<1
-      if (GV%Boussinesq) then
+      if (J <= MS%jsdw) then
+        ! Do not apply a Northern Flather OBC at the southern halo points on a PE, as doing so would
+        ! create a segmentation fault and this velocity will not propagate through to the next iteration.
+        ssh_in = BT_OBC%SSH_outer_v(i,J)
+      elseif (GV%Boussinesq) then
         ssh_in = GV%H_to_Z*(eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i,j-1)))      ! internal
       else
         ssh_1 = GV%H_to_RZ * eta(i,j) * SpV_avg(i,j) - (CS%bathyT(i,j) + G%Z_ref)
@@ -3926,7 +3938,11 @@ subroutine apply_v_velocity_OBCs(vbt, vhbt, vbt_trans, eta, SpV_avg, vbt_old, BT
     elseif (BT_OBC%v_OBC_type(i,J) == -FLATHER_OBC) then  ! Southern Flather OBC
       cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(i,J) ! CFL
       v_inlet = cfl*vbt_old(i,J+1) + (1.0-cfl)*vbt_old(i,J)  ! Valid for cfl <1
-      if (GV%Boussinesq) then
+      if (J >= MS%jedw-1) then
+        ! Do not apply a Southern Flather OBC at the northern halo points on a PE, as doing so would
+        ! create a segmentation fault and this velocity will not propagate through to the next iteration.
+        ssh_in = BT_OBC%SSH_outer_v(i,J)
+      elseif (GV%Boussinesq) then
         ssh_in = GV%H_to_Z*(eta(i,j+1) + (0.5-cfl)*(eta(i,j+1)-eta(i,j+2)))  ! internal
       else
         ssh_1 = GV%H_to_RZ * eta(i,j+1) * SpV_avg(i,j+1) - (CS%bathyT(i,j+1) + G%Z_ref)
@@ -3987,7 +4003,7 @@ subroutine initialize_BT_OBC(OBC, BT_OBC, G, CS)
   real :: OBC_sign  ! A sign encoding the direction of the OBC being used at a point [nondim]
   real :: OBC_type  ! A real copy of the integer encoding the type of OBC being used at a point [nondim]
   integer :: i, j, isdw, iedw, jsdw, jedw
-  integer :: l_seg
+  integer :: l_seg, Flather_OBC_in_halo
 
   isdw = CS%isdw ; iedw = CS%iedw ; jsdw = CS%jsdw ; jedw = CS%jedw
 
@@ -4037,27 +4053,48 @@ subroutine initialize_BT_OBC(OBC, BT_OBC, G, CS)
   BT_OBC%is_v_N_obc = iedw + 1 ; BT_OBC%ie_v_N_obc = isdw - 1
   BT_OBC%Js_v_N_obc = jedw + 1 ; BT_OBC%Je_v_N_obc = jsdw - 2
 
+  Flather_OBC_in_halo = 0
   do j=jsdw,jedw ; do I=isdw-1,iedw
     BT_OBC%u_OBC_type(I,j) = nint(u_OBC(I,j))
     if (BT_OBC%u_OBC_type(I,j) < 0) then ! This point has OBC_DIRECTION_W.
-      BT_OBC%Is_u_W_obc = min(I, BT_OBC%Is_u_W_obc) ; BT_OBC%Ie_u_W_obc = max(I, BT_OBC%Ie_u_W_obc)
-      BT_OBC%js_u_W_obc = min(j, BT_OBC%js_u_W_obc) ; BT_OBC%je_u_W_obc = max(j, BT_OBC%je_u_W_obc)
+      if ((BT_OBC%u_OBC_type(I,j) == -FLATHER_OBC) .and. (I >= iedw-1)) then
+        ! There is no need to specify the OBC at this point, but the stencil might need to be increased.
+        Flather_OBC_in_halo = 1
+      else
+        BT_OBC%Is_u_W_obc = min(I, BT_OBC%Is_u_W_obc) ; BT_OBC%Ie_u_W_obc = max(I, BT_OBC%Ie_u_W_obc)
+        BT_OBC%js_u_W_obc = min(j, BT_OBC%js_u_W_obc) ; BT_OBC%je_u_W_obc = max(j, BT_OBC%je_u_W_obc)
+      endif
     endif
     if (BT_OBC%u_OBC_type(I,j) > 0) then ! This point has OBC_DIRECTION_E.
-      BT_OBC%Is_u_E_obc = min(I, BT_OBC%Is_u_E_obc) ; BT_OBC%Ie_u_E_obc = max(I, BT_OBC%Ie_u_E_obc)
-      BT_OBC%js_u_E_obc = min(j, BT_OBC%js_u_E_obc) ; BT_OBC%je_u_E_obc = max(j, BT_OBC%je_u_E_obc)
+      if ((BT_OBC%u_OBC_type(I,j) == FLATHER_OBC) .and. (I <= isdw)) then
+        ! There is no need to specify the OBC at this point, but the stencil might need to be increased.
+        Flather_OBC_in_halo = 1
+      else
+        BT_OBC%Is_u_E_obc = min(I, BT_OBC%Is_u_E_obc) ; BT_OBC%Ie_u_E_obc = max(I, BT_OBC%Ie_u_E_obc)
+        BT_OBC%js_u_E_obc = min(j, BT_OBC%js_u_E_obc) ; BT_OBC%je_u_E_obc = max(j, BT_OBC%je_u_E_obc)
+      endif
     endif
   enddo ; enddo
 
   do J=jsdw-1,jedw ; do i=isdw,iedw
     BT_OBC%v_OBC_type(i,J) = nint(v_OBC(i,J))
-    if (BT_OBC%v_OBC_type(i,J) < 0) then ! This point has OBC_DIRECTION_S.
-      BT_OBC%is_v_S_obc = min(i, BT_OBC%is_v_S_obc) ; BT_OBC%ie_v_S_obc = max(i, BT_OBC%ie_v_S_obc)
-      BT_OBC%Js_v_S_obc = min(J, BT_OBC%Js_v_S_obc) ; BT_OBC%Je_v_S_obc = max(J, BT_OBC%Je_v_S_obc)
+    if (BT_OBC%v_OBC_type(i,J) < 0)  then ! This point has OBC_DIRECTION_S.
+      if ((BT_OBC%v_OBC_type(i,J) == -FLATHER_OBC) .and. (J >= jedw-1)) then
+        ! There is no need to specify the OBC at this point, but the stencil might need to be increased.
+        Flather_OBC_in_halo = 1
+      else
+        BT_OBC%is_v_S_obc = min(i, BT_OBC%is_v_S_obc) ; BT_OBC%ie_v_S_obc = max(i, BT_OBC%ie_v_S_obc)
+        BT_OBC%Js_v_S_obc = min(J, BT_OBC%Js_v_S_obc) ; BT_OBC%Je_v_S_obc = max(J, BT_OBC%Je_v_S_obc)
+      endif
     endif
     if (BT_OBC%v_OBC_type(i,J) > 0) then ! This point has OBC_DIRECTION_N.
-      BT_OBC%is_v_N_obc = min(i, BT_OBC%is_v_N_obc) ; BT_OBC%ie_v_N_obc = max(i, BT_OBC%ie_v_N_obc)
-      BT_OBC%Js_v_N_obc = min(J, BT_OBC%Js_v_N_obc) ; BT_OBC%Je_v_N_obc = max(J, BT_OBC%Je_v_N_obc)
+      if ((BT_OBC%v_OBC_type(i,J) == FLATHER_OBC) .and. (J <= jsdw)) then
+        ! There is no need to specify the OBC at this point, but the stencil might need to be increased.
+        Flather_OBC_in_halo = 1
+      else
+        BT_OBC%is_v_N_obc = min(i, BT_OBC%is_v_N_obc) ; BT_OBC%ie_v_N_obc = max(i, BT_OBC%ie_v_N_obc)
+        BT_OBC%Js_v_N_obc = min(J, BT_OBC%Js_v_N_obc) ; BT_OBC%Je_v_N_obc = max(J, BT_OBC%Je_v_N_obc)
+      endif
     endif
   enddo ; enddo