pure subroutine ks_gather_corner(nx, ny, nz, ic, jc, h_layer, u_face, &
v_face, hT, hS, wet_t, wet_u, wet_v, &
h_sd, u_sd, v_sd, t_sd, s_sd)
!! Assemble the surface-down column at corner (ic,jc) from the
!! 2x2 cell patch + 4 adjacent faces (JHL08 vertex form; the
!! interpolation recipes of the reference implementation):
!! u,v : 2-point THICKNESS-weighted average across the corner,
!! with the face thickness itself a mask-weighted 2-cell
!! average (recomputed inline — deterministic, so the
!! repeated evaluation is bitwise identical to MOM6's
!! precomputed h_at_u/h_at_v Pass A, non-OBC-bug form).
!! T,S : 4-cell mask-AND-thickness-weighted average. The
!! registry stores hTr = h*T, which is exactly the
!! weighted quantity, so we sum wet*hTr directly.
!! h : 4-cell mask-weighted average (no thickness weight —
!! it IS the thickness).
!! Returns RAW h (no floor) — the caller decides floor vs
!! massless-merge exactly as the column path does. The deliberate
!! (SW+NE)+(SE+NW) bracketing is reproducible-sum ordering — do
!! not reassociate.
!$acc routine seq
integer, intent(in) :: nx, ny, nz, ic, jc
real(wp), intent(in) :: h_layer(nx, ny, nz)
real(wp), intent(in) :: u_face(nx + 1, ny, nz)
real(wp), intent(in) :: v_face(nx, ny + 1, nz)
real(wp), intent(in) :: hT(nx, ny, nz)
real(wp), intent(in) :: hS(nx, ny, nz)
real(wp), intent(in) :: wet_t(nx, ny)
real(wp), intent(in) :: wet_u(nx + 1, ny)
real(wp), intent(in) :: wet_v(nx, ny + 1)
real(wp), intent(out) :: h_sd(NZL), u_sd(NZL), v_sd(NZL)
real(wp), intent(out) :: t_sd(NZL), s_sd(NZL)
integer :: k, kg
real(wp) :: w_sw, w_se, w_nw, w_ne
real(wp) :: h_sw, h_se, h_nw, h_ne
real(wp) :: hu_s, hu_n, hv_w, hv_e
real(wp) :: hwt, i_hwt
w_sw = wet_t(ic - 1, jc - 1)
w_se = wet_t(ic, jc - 1)
w_nw = wet_t(ic - 1, jc)
w_ne = wet_t(ic, jc)
do k = 1, nz
kg = nz + 1 - k ! global bottom-up layer -> local surface-down
h_sw = h_layer(ic - 1, jc - 1, kg)
h_se = h_layer(ic, jc - 1, kg)
h_nw = h_layer(ic - 1, jc, kg)
h_ne = h_layer(ic, jc, kg)
! Face thicknesses (mask-weighted 2-cell averages, inline).
hu_s = wet_u(ic, jc - 1)*(w_sw*h_sw + w_se*h_se)/ &
(w_sw + w_se + MASK_SUM_EPS)
hu_n = wet_u(ic, jc)*(w_nw*h_nw + w_ne*h_ne)/ &
(w_nw + w_ne + MASK_SUM_EPS)
hv_w = wet_v(ic - 1, jc)*(w_sw*h_sw + w_nw*h_nw)/ &
(w_sw + w_nw + MASK_SUM_EPS)
hv_e = wet_v(ic, jc)*(w_se*h_se + w_ne*h_ne)/ &
(w_se + w_ne + MASK_SUM_EPS)
! Thickness-weighted 2-point transverse velocity averages.
u_sd(k) = ((u_face(ic, jc - 1, kg)*hu_s) + &
(u_face(ic, jc, kg)*hu_n))/ &
((hu_s + hu_n) + H_TINY_CORNER)
v_sd(k) = ((v_face(ic - 1, jc, kg)*hv_w) + &
(v_face(ic, jc, kg)*hv_e))/ &
((hv_w + hv_e) + H_TINY_CORNER)
! 4-cell mask*thickness weight (diagonal pairs first).
hwt = ((w_sw*h_sw + w_ne*h_ne) + (w_se*h_se + w_nw*h_nw))
i_hwt = 1.0_wp/(hwt + H_DIV_EPS)
t_sd(k) = ((w_sw*hT(ic - 1, jc - 1, kg) + w_ne*hT(ic, jc, kg)) + &
(w_se*hT(ic, jc - 1, kg) + w_nw*hT(ic - 1, jc, kg)))*i_hwt
s_sd(k) = ((w_sw*hS(ic - 1, jc - 1, kg) + w_ne*hS(ic, jc, kg)) + &
(w_se*hS(ic, jc - 1, kg) + w_nw*hS(ic - 1, jc, kg)))*i_hwt
! 4-cell mask-weighted thickness (mean of the WET cells only —
! a partially-wet corner is not diluted by land zeros).
h_sd(k) = hwt/((w_sw + w_ne) + (w_se + w_nw) + MASK_SUM_EPS)
end do
end subroutine ks_gather_corner