subroutine ocean_obc_fill_ghosts(grid, bc, ms)
!! Fill h_layer and tracer hTr ghosts at open-ish edges.
!! h_layer: zero-gradient (copy adjacent interior column).
!! Tracer hTr per (j,k): outflow ⇒ ghost := interior (zero-gradient);
!! inflow ⇒ ghost hTr := clamped_tracer(it) * h_ghost. Inflow criterion
!! uses the per-layer wall-face velocity sign (outward-normal convention):
!! west inflow u_wall>0, east u_wall<0, south v_wall>0, north v_wall<0.
!! No-op when no edge is open-ish. Per-tracer loop outside the DCs.
type(hgrid_t), intent(in) :: grid
type(ocean_bc_state_t), intent(in) :: bc
type(multilayer_state_t), intent(inout) :: ms
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
! MPI-seam neutralisation (O0): a seam edge carries no physical BC.
! WALL is a no-op throughout this routine (verified: all dispatch is
! guarded by is_open_ish, which excludes OBC_WALL), so remapping the
! cached tag makes every per-edge block skip the seam.
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
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 ----
call fill_h_layer_ghosts(ms%h_layer, nxt, nyt, nz, &
i_w, i_e, j_s, j_n, grid%nghost, &
bc_w, bc_e, bc_s, bc_n)
! ---- Tracer hTr ghosts: upwind-aware or reservoir-based ----
! Reservoir active (tres_* allocated): unconditional hTr_ghost = tres*h_ghost
! (Marchesiello et al. 2001). Absent (default): sign-switch path,
! bit-identical to prior behaviour. Dispatch is per-edge-per-direction.
if (allocated(ms%tracers)) then
do it = 1, size(ms%tracers)
if (.not. allocated(ms%tracers(it)%hTr)) cycle
! Zero-gradient CONCENTRATION pre-fill over the FULL open-edge ghost
! region (corners included): the per-edge upwind fills below only
! touch the physical cross-extent, so unfilled ghost×ghost corners
! would accumulate garbage under inflow and blow up T/S. Upwind
! fills then overwrite the physical edges, leaving corners zero-grad.
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)
! ---- West ghost fill ----
if (is_open_ish(bc_w)) then
if (allocated(bc%tres_west)) then
call fill_ghost_west_res( &
ms%tracers(it)%hTr, ms%h_layer, bc%tres_west, &
nxt, nyt, nz, i_w, j_s, j_n, grid%nghost, it, size(bc%tres_west, 3))
else
call fill_ghost_west_sign(ms%tracers(it)%hTr, ms%h_layer, &
ms%u_face_x_layer, &
nxt, nyt, nz, i_w, j_s, j_n, grid%nghost, &
get_clamped_tracer(bc%west%clamped_tracer, it))
end if
end if
! ---- East ghost fill ----
if (is_open_ish(bc_e)) then
if (allocated(bc%tres_east)) then
call fill_ghost_east_res( &
ms%tracers(it)%hTr, ms%h_layer, bc%tres_east, &
nxt, nyt, nz, i_e, j_s, j_n, grid%nghost, it, size(bc%tres_east, 3))
else
call fill_ghost_east_sign(ms%tracers(it)%hTr, ms%h_layer, &
ms%u_face_x_layer, &
nxt, nyt, nz, i_e, j_s, j_n, grid%nghost, &
get_clamped_tracer(bc%east%clamped_tracer, it))
end if
end if
! ---- South ghost fill ----
if (is_open_ish(bc_s)) then
if (allocated(bc%tres_south)) then
call fill_ghost_south_res( &
ms%tracers(it)%hTr, ms%h_layer, bc%tres_south, &
nxt, nyt, nz, j_s, i_w, i_e, grid%nghost, it, size(bc%tres_south, 3))
else
call fill_ghost_south_sign(ms%tracers(it)%hTr, ms%h_layer, &
ms%v_face_y_layer, &
nxt, nyt, nz, j_s, i_w, i_e, grid%nghost, &
get_clamped_tracer(bc%south%clamped_tracer, it))
end if
end if
! ---- North ghost fill ----
if (is_open_ish(bc_n)) then
if (allocated(bc%tres_north)) then
call fill_ghost_north_res( &
ms%tracers(it)%hTr, ms%h_layer, bc%tres_north, &
nxt, nyt, nz, j_n, i_w, i_e, grid%nghost, it, size(bc%tres_north, 3))
else
call fill_ghost_north_sign(ms%tracers(it)%hTr, ms%h_layer, &
ms%v_face_y_layer, &
nxt, nyt, nz, j_n, i_w, i_e, grid%nghost, &
get_clamped_tracer(bc%north%clamped_tracer, it))
end if
end if
end do
end if
end subroutine ocean_obc_fill_ghosts