pure subroutine set_local_BT_cont_types(grid, metrics, bt_work, ms, dt_outer)
!! Populate `bt_work%BTCL_u/v` — the per-face flux-closure coefficients
!! consumed by `find_uhbt` — from the current `h_layer`. Upstream-h-sum
!! approach: FA_u_W0=FA_u_E0=Σ_k h_face (centred); FA_u_WW=Σ_k h_layer(west)
!! and FA_u_EE=Σ_k h_layer(east) (saturated-regime upstream draw); uBT_WW/EE
!! = ±VOL_CFL·dx/dt_outer pin the saturation velocity; uh_crv*/uh_** are the
!! C¹-matching coefficients (Hallberg & Adcroft 2009). No-op when
!! `use_bt_cont_type = .false.` (BTCL_u/v unallocated).
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(barotropic_workstate_t), intent(inout) :: bt_work
type(multilayer_state_t), intent(in) :: ms
real(wp), intent(in) :: dt_outer
real(wp), parameter :: C1_3 = 1.0_wp/3.0_wp
integer :: i, j, k, nx, ny, nz
real(wp) :: fa_centre, fa_up_W, fa_up_E, fa_up_N, fa_up_S
real(wp) :: h_face, u_cfl_x, u_cfl_y, inv_ucfl_x2, inv_ucfl_y2
real(wp) :: inv_dt
! u_cfl_x/y + inv_ucfl_x2/y2 are written per-iteration as DC locals
! (declared here so the `local()` clause can name them).
! Boundary u/v-faces are not written — they keep their type-default zero;
! find_uhbt(0, zero_BTC) = 0, matching the wall-zero ubt they carry.
if (.not. bt_work%use_bt_cont_type) return
nx = grid%nx_total
ny = grid%ny_total
nz = ms%nz_ml
inv_dt = 1.0_wp/dt_outer
! ---- u-faces (i in 2..nx) ----
! Per-face CFL velocity BTC_VOL_CFL·dxCu/dt_outer.
do concurrent(j=1:ny, i=2:nx) &
local(k, fa_centre, fa_up_W, fa_up_E, h_face, u_cfl_x, inv_ucfl_x2)
u_cfl_x = BTC_VOL_CFL*metrics%dxCu(i, j)*inv_dt
inv_ucfl_x2 = 0.0_wp
if (u_cfl_x > 0.0_wp) inv_ucfl_x2 = 1.0_wp/(u_cfl_x*u_cfl_x)
fa_centre = 0.0_wp
fa_up_W = 0.0_wp
fa_up_E = 0.0_wp
do k = 1, nz
h_face = 0.5_wp*(ms%h_layer(i - 1, j, k) + ms%h_layer(i, j, k))
if (h_face > BTC_H_NEGLECT) fa_centre = fa_centre + h_face
if (ms%h_layer(i - 1, j, k) > BTC_H_NEGLECT) then
fa_up_W = fa_up_W + ms%h_layer(i - 1, j, k)
end if
if (ms%h_layer(i, j, k) > BTC_H_NEGLECT) then
fa_up_E = fa_up_E + ms%h_layer(i, j, k)
end if
end do
bt_work%BTCL_u(i, j)%FA_u_W0 = fa_centre
bt_work%BTCL_u(i, j)%FA_u_E0 = fa_centre
bt_work%BTCL_u(i, j)%FA_u_WW = fa_up_W
bt_work%BTCL_u(i, j)%FA_u_EE = fa_up_E
bt_work%BTCL_u(i, j)%uBT_WW = u_cfl_x
bt_work%BTCL_u(i, j)%uBT_EE = -u_cfl_x
bt_work%BTCL_u(i, j)%uh_crvW = C1_3*(fa_up_W - fa_centre)*inv_ucfl_x2
bt_work%BTCL_u(i, j)%uh_crvE = C1_3*(fa_up_E - fa_centre)*inv_ucfl_x2
bt_work%BTCL_u(i, j)%uh_WW = u_cfl_x*C1_3*(2.0_wp*fa_centre + fa_up_W)
bt_work%BTCL_u(i, j)%uh_EE = -u_cfl_x*C1_3*(2.0_wp*fa_centre + fa_up_E)
end do
! ---- v-faces (j in 2..ny) ----
do concurrent(j=2:ny, i=1:nx) &
local(k, fa_centre, fa_up_N, fa_up_S, h_face, u_cfl_y, inv_ucfl_y2)
u_cfl_y = BTC_VOL_CFL*metrics%dyCv(i, j)*inv_dt
inv_ucfl_y2 = 0.0_wp
if (u_cfl_y > 0.0_wp) inv_ucfl_y2 = 1.0_wp/(u_cfl_y*u_cfl_y)
fa_centre = 0.0_wp
fa_up_N = 0.0_wp
fa_up_S = 0.0_wp
do k = 1, nz
h_face = 0.5_wp*(ms%h_layer(i, j - 1, k) + ms%h_layer(i, j, k))
if (h_face > BTC_H_NEGLECT) fa_centre = fa_centre + h_face
if (ms%h_layer(i, j - 1, k) > BTC_H_NEGLECT) then
fa_up_S = fa_up_S + ms%h_layer(i, j - 1, k)
end if
if (ms%h_layer(i, j, k) > BTC_H_NEGLECT) then
fa_up_N = fa_up_N + ms%h_layer(i, j, k)
end if
end do
! Convention: vBT_SS > 0 → southward-draw saturation (v > 0 flow,
! upstream is cell j-1, the "S" column). vBT_NN < 0 → northward-
! draw saturation. Mirrors MOM6's local_BT_cont_v_type.
bt_work%BTCL_v(i, j)%FA_v_S0 = fa_centre
bt_work%BTCL_v(i, j)%FA_v_N0 = fa_centre
bt_work%BTCL_v(i, j)%FA_v_SS = fa_up_S
bt_work%BTCL_v(i, j)%FA_v_NN = fa_up_N
bt_work%BTCL_v(i, j)%vBT_SS = u_cfl_y
bt_work%BTCL_v(i, j)%vBT_NN = -u_cfl_y
bt_work%BTCL_v(i, j)%vh_crvS = C1_3*(fa_up_S - fa_centre)*inv_ucfl_y2
bt_work%BTCL_v(i, j)%vh_crvN = C1_3*(fa_up_N - fa_centre)*inv_ucfl_y2
bt_work%BTCL_v(i, j)%vh_SS = u_cfl_y*C1_3*(2.0_wp*fa_centre + fa_up_S)
bt_work%BTCL_v(i, j)%vh_NN = -u_cfl_y*C1_3*(2.0_wp*fa_centre + fa_up_N)
end do
! Boundary u-faces (i=1, i=nx_face) and v-faces (j=1, j=ny_face)
! stay at their type-default zero — wall faces in our setup carry
! ubt=0 and find_uhbt(0, anything)=0, so no flux through them.
end subroutine set_local_BT_cont_types