pure subroutine wavespeed_cg1_column(nz, h_rak, rho_rak, rho0, cg1)
!! First-baroclinic wave speed for ONE column. `h_rak`/`rho_rak`
!! are in Roundabout ordering (k=1 bed, k=nz surface), fixed-size
!! NZ_STACK_MAX arrays; only 1..nz are read. Returns `cg1` (m/s),
!! 0 for land / homogeneous / `kc<2` / sub-floor columns.
!! Gathers+flips surface-down, backtracking convective merge,
!! symmetric tridiag, fixed-budget Sturm-count bisection.
!$acc routine seq
integer, intent(in) :: nz
real(wp), intent(in) :: h_rak(NZ_STACK_MAX), rho_rak(NZ_STACK_MAX)
real(wp), intent(in) :: rho0
real(wp), intent(out) :: cg1
! column arrays (surface-down): 5 total (spec §7)
real(wp) :: hc(NZ_STACK_MAX), rc(NZ_STACK_MAX)
real(wp) :: gprime(NZ_STACK_MAX + 1)
real(wp) :: igu(NZ_STACK_MAX + 1), igl(NZ_STACK_MAX + 1)
! scalar carries (names avoid intrinsic shadows: no count/sum/sign/mod/...)
integer :: k, kg, kc, ki, it, n_chg, slo, smid
real(wp) :: g_rho0, drxh, hnew, hnew2, rnew2, gp
real(wp) :: s2tot, lam_lo, lam_hi, lam_mid, lam_probe
logical :: do_weld, do_weld2
cg1 = 0.0_wp
if (nz < 2) return
g_rho0 = GRAVITY/rho0
! ---- gather (FLIP): rdb k=1(bed)..nz(surf) -> surface-down ----
do k = 1, nz
kg = nz + 1 - k ! k=1 -> surface layer
hc(k) = h_rak(kg)
rc(k) = rho_rak(kg)
end do
! ---- drxh_sum on the RAW (unmerged) surface-down profile ----
drxh = 0.0_wp
do k = 2, nz
drxh = drxh + 0.5_wp*(hc(k - 1) + hc(k))*max(0.0_wp, rc(k) - rc(k - 1))
end do
! ---- backtracking convective merge (last-committed compare) ----
! kc counts committed layers; the committed layer kc lives at
! hc(kc)/rc(kc). Walk the remaining raw layers (still in hc/rc at
! index k) and weld or commit. We compact in place: the committed
! stack occupies hc(1..kc); the next raw layer is hc(k) (k>kc).
kc = 1
do k = 2, nz
do_weld = ((rc(k) - rc(kc))*(hc(kc) + hc(k)) < 2.0_wp*TOL_MERGE*drxh)
if (do_weld) then
hnew = hc(kc) + hc(k)
rc(kc) = (hc(kc)*rc(kc) + hc(k)*rc(k))/hnew
hc(kc) = hnew
! backtrack: undo any inversion the weld created (looser tol)
do
if (kc < 2) exit
do_weld2 = ((rc(kc) - rc(kc - 1))*(hc(kc) + hc(kc - 1)) &
< TOL_MERGE*drxh)
if (.not. do_weld2) exit
hnew2 = hc(kc) + hc(kc - 1)
rnew2 = (hc(kc)*rc(kc) + hc(kc - 1)*rc(kc - 1))/hnew2
rc(kc - 1) = rnew2
hc(kc - 1) = hnew2
kc = kc - 1
end do
else
kc = kc + 1
hc(kc) = hc(k)
rc(kc) = rc(k)
end if
end do
if (kc < 2) return
! ---- gprime, Igu, Igl at interfaces K=2..kc ----
s2tot = 0.0_wp
do ki = 2, kc
gp = g_rho0*(rc(ki) - rc(ki - 1))
gprime(ki) = gp
if (gp > 0.0_wp) then
igu(ki) = 1.0_wp/(gp*hc(ki - 1))
igl(ki) = 1.0_wp/(gp*hc(ki))
s2tot = s2tot + gp*(hc(ki - 1) + hc(ki))
else
igu(ki) = 0.0_wp
igl(ki) = 0.0_wp
end if
end do
if (s2tot <= MIN_SPEED2) return
! ---- bracket via Sturm-count doubling (robust for all nz) ----
lam_probe = LAM_SEED
do it = 1, MAX_DBL
if (lam_probe >= 1.0_wp/MIN_SPEED2) then
lam_probe = 1.0_wp/MIN_SPEED2
exit
end if
n_chg = sturm_count(igu, igl, kc, lam_probe)
if (n_chg >= 1) exit
lam_probe = lam_probe*2.0_wp
end do
! safety: no mode-1 found below the cap
if (sturm_count(igu, igl, kc, lam_probe) < 1) return
lam_lo = lam_probe*0.5_wp
lam_hi = lam_probe
! ---- fixed-budget bisection on sign(det) (NO early exit) ----
slo = det_sign(igu, igl, kc, lam_lo)
do it = 1, MAX_ITT
lam_mid = 0.5_wp*(lam_lo + lam_hi)
smid = det_sign(igu, igl, kc, lam_mid)
if (smid == slo) then
lam_lo = lam_mid
else
lam_hi = lam_mid
end if
end do
lam_mid = 0.5_wp*(lam_lo + lam_hi)
if (lam_mid > 0.0_wp) cg1 = 1.0_wp/sqrt(lam_mid)
end subroutine wavespeed_cg1_column