fill_uv_layer_ghosts Subroutine

   pure subroutine fill_uv_layer_ghosts(u_layer, v_layer, &
                                        nx_total, ny_total, nz, nghost, &
                                        nx_phys, ny_phys, bc_w, bc_e, bc_s, bc_n)
      !! Zero-gradient fill of per-layer face velocities into the OPEN-edge
      !! ghost region, corners included.  u_layer (nx_total+1, ny_total, nz);
      !! v_layer (nx_total, ny_total+1, nz).  x-pass over the full cross extent;
      !! y-pass clipped to the physical span unless the adjacent x-edge is open.
      !! Gated on open-ish tags ⇒ WALL/PERIODIC ⇒ no DC ⇒ bit-identical.
      integer, intent(in) :: nx_total, ny_total, nz, nghost, nx_phys, ny_phys
      real(wp), intent(inout) :: u_layer(nx_total + 1, ny_total, nz)
      real(wp), intent(inout) :: v_layer(nx_total, ny_total + 1, nz)
      integer, intent(in) :: bc_w, bc_e, bc_s, bc_n
      integer :: i, j, k, g, ng
      integer :: iu_lo, iu_hi, iv_lo, iv_hi
      logical :: fw, fe, fs, fn

      ng = nghost
      ! Fill set = SAME as the barotropic-substep η-ghost fill (OPEN/TIDAL/
      ! CHAPMAN/CLAMPED), deliberately NOT `is_open_ish` (excludes NESTED, which
      ! the C-grid barotropic treats as a wall) so velocity and η fills stay
      ! consistent and don't re-trip the corner NaN.
      fw = open_ghost_fill_edge(bc_w)
      fe = open_ghost_fill_edge(bc_e)
      fs = open_ghost_fill_edge(bc_s)
      fn = open_ghost_fill_edge(bc_n)

      ! ---- x-pass: west / east, over the full j extent (covers corners) ----
      if (fw) then
         do concurrent(k=1:nz, j=1:ny_total, g=1:ng)
            u_layer(g, j, k) = u_layer(ng + 1, j, k)             ! u faces 1..ng <- wall face ng+1
         end do
         do concurrent(k=1:nz, j=1:ny_total + 1, g=1:ng)
            v_layer(g, j, k) = v_layer(ng + 1, j, k)             ! v cols 1..ng <- first phys col
         end do
      end if
      if (fe) then
         do concurrent(k=1:nz, j=1:ny_total, g=1:ng)
            u_layer(ng + nx_phys + 1 + g, j, k) = u_layer(ng + nx_phys + 1, j, k)  ! east ghost faces
         end do
         do concurrent(k=1:nz, j=1:ny_total + 1, g=1:ng)
            v_layer(ng + nx_phys + g, j, k) = v_layer(ng + nx_phys, j, k)          ! east ghost cols
         end do
      end if

      ! ---- y-pass: south / north, over the corner-safe i extent ----
      ! Extend the y-pass into an x-ghost column for any NON-PERIODIC x-edge so a
      ! south/north-open edge fills its ghost CORNER even against a WALL side.  A
      ! PERIODIC x-edge is excluded — it owns its ghost via the wrap (preserves
      ! periodic+open mixed configs, e.g. Eady).
      iu_lo = ng + 1                 ! u face i-range (physical span by default)
      iu_hi = ng + nx_phys + 1
      iv_lo = ng + 1                 ! v cell i-range (physical span by default)
      iv_hi = ng + nx_phys
      if (bc_w /= OBC_PERIODIC) then
         iu_lo = 1
         iv_lo = 1
      end if
      if (bc_e /= OBC_PERIODIC) then
         iu_hi = nx_total + 1
         iv_hi = nx_total
      end if
      if (fs) then
         do concurrent(k=1:nz, g=1:ng, i=iu_lo:iu_hi)
            u_layer(i, g, k) = u_layer(i, ng + 1, k)            ! u south ghost rows <- first phys row
         end do
         do concurrent(k=1:nz, g=1:ng, i=iv_lo:iv_hi)
            v_layer(i, g, k) = v_layer(i, ng + 1, k)            ! v south ghost faces <- wall face
         end do
      end if
      if (fn) then
         do concurrent(k=1:nz, g=1:ng, i=iu_lo:iu_hi)
            u_layer(i, ng + ny_phys + g, k) = u_layer(i, ng + ny_phys, k)          ! u north ghost rows
         end do
         do concurrent(k=1:nz, g=1:ng, i=iv_lo:iv_hi)
            v_layer(i, ng + ny_phys + 1 + g, k) = v_layer(i, ng + ny_phys + 1, k)  ! v north ghost faces
         end do
      end if
   end subroutine fill_uv_layer_ghosts