ocean_obc_refill_ghost_ssh Subroutine

   subroutine ocean_obc_refill_ghost_ssh(grid, bc, ms, bt_H_ref)
      !! Re-establish a zero-gradient free surface in the open-edge GHOST
      !! columns, called at the END of the outer step (after the ALE remap)
      !! so the diagnostic manager sees a consistent halo.
      !!
      !! Why this is needed: the ghost cells are updated by the slow
      !! continuity over the full array (`i = 1..nx_total`) using the
      !! zeroed array-edge fluxes, and the conservative ALE remap
      !! (`remap_h_ref = Σh − bt_eta`) PRESERVES that spurious ghost
      !! transport divergence.  Nothing re-imposes the open-boundary
      !! zero-gradient invariant between the remap and the diagnostic
      !! read, so `SSH = Σ_k h_layer − b` in the ghosts drifts to tens of
      !! metres (worst where the boundary bathymetry is steep and at the
      !! corners where two open edges meet), while the physical interior
      !! stays healthy.
      !!
      !! Fix: overwrite each open-edge ghost column with the nearest
      !! interior column SCALED so the ghost column total equals
      !! `bt_H_ref(ghost) + η_interior`, i.e. the free-surface anomaly is
      !! flat across the open boundary.  Since `bt_H_ref == b`, this makes
      !! `SSH_ghost = η_interior` for ANY bathymetry — including formula
      !! topographies whose ghost `b` differs from the first interior cell.
      !! For file bathymetry (ghost `b` constant-extrapolated) the scale
      !! collapses to 1, i.e. a plain zero-gradient re-copy that simply
      !! discards the accumulated ghost drift.
      !!
      !! Physics-neutral: the next outer step's stage-1 `ocean_obc_fill_ghosts`
      !! re-fills the ghost h_layer (thickness copy) BEFORE any physics
      !! kernel reads it, so this pass only affects what the diagnostics
      !! (and other halo consumers) observe — the prognostic trajectory is
      !! bit-identical.  Gated on open-ish edges ⇒ WALL/PERIODIC ⇒ no-op.
      type(hgrid_t), intent(in) :: grid
      type(ocean_bc_state_t), intent(in) :: bc
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: bt_H_ref(grid%nx_total, grid%ny_total)
         !! Mode-split reference column depth (== seeded bathymetry `b`).

      integer :: bc_w, bc_e, bc_s, bc_n
      logical :: any_open
      integer :: it, nz, nxt, nyt, nx, ny
      integer :: i_w, i_e, j_s, j_n

      bc_w = bc%west%bc_type
      bc_e = bc%east%bc_type
      bc_s = bc%south%bc_type
      bc_n = bc%north%bc_type
      ! An MPI seam is not an open edge: its ghost columns are a neighbour's
      ! interior, filled by the halo exchange.  Overwriting them with the
      ! zero-gradient copy made every decomposed OBC run diverge from the
      ! serial one (same gate as `ocean_obc_fill_ghosts`).
      if (.not. bc%has_west) bc_w = OBC_WALL
      if (.not. bc%has_east) bc_e = OBC_WALL
      if (.not. bc%has_south) bc_s = OBC_WALL
      if (.not. bc%has_north) bc_n = OBC_WALL

      any_open = is_open_ish(bc_w) .or. is_open_ish(bc_e) .or. &
                 is_open_ish(bc_s) .or. is_open_ish(bc_n)
      if (.not. any_open) return
      if (.not. allocated(ms%h_layer)) return

      nz = ms%nz_ml
      nxt = grid%nx_total
      nyt = grid%ny_total
      nx = grid%nx_phys
      ny = grid%ny_phys
      i_w = grid%nghost + 1          ! west physical wall-face / first interior cell
      i_e = grid%nghost + nx         ! last interior cell (east)
      j_s = grid%nghost + 1          ! south physical wall-face / first interior cell
      j_n = grid%nghost + ny         ! last interior cell (north)

      ! ---- h_layer ghosts: zero-gradient free surface (scaled to ghost b) ----
      call refill_h_ghost_scaled(ms%h_layer, bt_H_ref, nxt, nyt, nz, &
                                 i_w, i_e, j_s, j_n, grid%nghost, &
                                 bc_w, bc_e, bc_s, bc_n)

      ! ---- Tracer hTr ghosts: zero-gradient concentration onto the new h ----
      if (allocated(ms%tracers)) then
         do it = 1, size(ms%tracers)
            if (.not. allocated(ms%tracers(it)%hTr)) cycle
            call fill_tracer_ghosts_zerograd(ms%tracers(it)%hTr, ms%h_layer, &
                                             nxt, nyt, nz, i_w, i_e, j_s, j_n, &
                                             grid%nghost, bc_w, bc_e, bc_s, bc_n)
         end do
      end if
   end subroutine ocean_obc_refill_ghost_ssh