From b7ba4ca0b1fa2762f3ce9a623942ee859d857c95 Mon Sep 17 00:00:00 2001
From: Robert Hallberg <Robert.Hallberg@noaa.gov>
Date: Thu, 17 Apr 2025 18:14:59 -0400
Subject: [PATCH] Zero out bt_pgf_u at OBC points

  Add code to zero out bt_pgf_u and bt_pgf_v at open boundary condition points,
thereby altering the pressure diagnostics that are derived from these arrays.
At OBCs, the pressure force is not well defined, and in fact the u_accel_bt and
v_accel_bt arrays upon which recently added calculation of bt_pgf_u and bt_pgf_v
are based is immediately replaced by a different OBC-specific expression for
the accelerations at OBC points.  Moreover, at OBC points bt_pgf_u and bt_pgf_v
were being calculated using values of eta_anom that were projected outward
across OBCs, but where this means that two different OBCs could be trying to set
these projected values of eta_anom to two different values (e.g, at a convex
corner in OBC segments or at back-to-back OBC segments separated by a single
land point). In such cases the diagnostic becomes indeterminate depending on
which OBC segment is applied first.  Zeroing out this diagnostic removes this
ambiguity.  This change does change (and correct) some diagnostics in energy
budgets at OBC points, but all solutions are bitwise identical and public
interfaces are unaltered.
---
 src/core/MOM_barotropic.F90 | 36 ++++++++++++++++++++++++++----------
 1 file changed, 26 insertions(+), 10 deletions(-)

diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 6d8370062a..b8824abeae 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -1966,16 +1966,32 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
     endif
 
     ! Copy decomposed barotropic accelerations to ADp
-    if (associated(ADp%bt_pgf_u)) then ; do k=1,nz ; do j=js,je ; do I=is-1,ie
-      ADp%bt_pgf_u(I,j,k) = PFu_avg(I,j) - &
-        (((pbce(i+1,j,k) - gtot_W(i+1,j)) * e_anom(i+1,j)) - &
-         ((pbce(i,j,k) - gtot_E(i,j)) * e_anom(i,j))) * CS%IdxCu(I,j)
-    enddo ; enddo ; enddo ; endif
-    if (associated(ADp%bt_pgf_v)) then ; do k=1,nz ; do J=js-1,je ; do i=is,ie
-      ADp%bt_pgf_v(i,J,k) = PFv_avg(i,J) - &
-        (((pbce(i,j+1,k) - gtot_S(i,j+1)) * e_anom(i,j+1)) - &
-         ((pbce(i,j,k) - gtot_N(i,j)) * e_anom(i,j))) * CS%IdyCv(i,J)
-    enddo ; enddo ; enddo ; endif
+    if (associated(ADp%bt_pgf_u)) then
+      do k=1,nz ; do j=js,je ; do I=is-1,ie
+        ADp%bt_pgf_u(I,j,k) = PFu_avg(I,j) - &
+          (((pbce(i+1,j,k) - gtot_W(i+1,j)) * e_anom(i+1,j)) - &
+           ((pbce(i,j,k) - gtot_E(i,j)) * e_anom(i,j))) * CS%IdxCu(I,j)
+      enddo ; enddo ; enddo
+      ! The pressure gradient at OBC points is not meaningful and needs to be zeroed out.
+      if (CS%BT_OBC%apply_u_OBCs) then ; do j=js,je ; do I=is-1,ie
+        if (OBC%segnum_u(I,j) /= OBC_NONE) then
+          do k=1,nz ; ADp%bt_pgf_u(I,j,k) = 0.0 ; enddo
+        endif
+      enddo ; enddo ; endif
+    endif
+    if (associated(ADp%bt_pgf_v)) then
+      do k=1,nz ; do J=js-1,je ; do i=is,ie
+        ADp%bt_pgf_v(i,J,k) = PFv_avg(i,J) - &
+          (((pbce(i,j+1,k) - gtot_S(i,j+1)) * e_anom(i,j+1)) - &
+           ((pbce(i,j,k) - gtot_N(i,j)) * e_anom(i,j))) * CS%IdyCv(i,J)
+      enddo ; enddo ; enddo
+      ! The pressure gradient at OBC points is not meaningful and needs to be zeroed out.
+      if (CS%BT_OBC%apply_v_OBCs) then ; do J=js-1,je ; do i=is,ie
+        if (OBC%segnum_v(i,J) /= OBC_NONE) then
+          do k=1,nz ; ADp%bt_pgf_v(i,J,k) = 0.0 ; enddo
+        endif
+      enddo ; enddo ; endif
+    endif
 
     if (associated(ADp%bt_cor_u)) then ; do j=js,je ; do I=is-1,ie
       ADp%bt_cor_u(I,j) = Coru_avg(I,j)