pure subroutine ks_projected_state(nz, dt_now, ks, ke, vel_underflow, &
dbuoy_t, dbuoy_s, h_sd, idz_int_s, &
u0, v0, t0, s0, kappa_ps, &
u_o, v_o, t_o, s_o, c1_o, n2_o, s2_o)
!! Mix (u0,v0,T0,S0) implicitly with `kappa_ps` over `dt_now`,
!! restricted to the layer band [ks,ke], and recompute N^2/S^2 at
!! interfaces (band edges blend mixed inside / original outside).
!! Backward-Euler tridiagonal; no-slip for u,v iff the band
!! reaches the bed (ke==nz), insulating T,S. N^2 floored at 0
!! (design doc section 5.5).
!$acc routine seq
integer, intent(in) :: nz, ks, ke
real(wp), intent(in) :: dt_now, vel_underflow
real(wp), intent(in) :: dbuoy_t(NZLI), dbuoy_s(NZLI)
real(wp), intent(in) :: h_sd(NZL), idz_int_s(NZLI)
real(wp), intent(in) :: u0(NZL), v0(NZL), t0(NZL), s0(NZL)
real(wp), intent(in) :: kappa_ps(NZLI)
real(wp), intent(out) :: u_o(NZL), v_o(NZL), t_o(NZL), s_o(NZL)
real(wp), intent(out) :: c1_o(NZLI), n2_o(NZLI), s2_o(NZLI)
integer :: k, kk
real(wp) :: a_b, a_a, b1, b1nz, d1, bd1
real(wp) :: ua, ub, va, vb, ta, tb, sa, sb, n2v
do k = 1, nz
u_o(k) = u0(k)
v_o(k) = v0(k)
t_o(k) = t0(k)
s_o(k) = s0(k)
end do
do k = 1, nz + 1
c1_o(k) = 0.0_wp
end do
if (ks <= ke .and. dt_now > 0.0_wp) then
! Forward sweep (top layer ks).
a_b = dt_now*kappa_ps(ks + 1)*idz_int_s(ks + 1)
b1 = 1.0_wp/(h_sd(ks) + a_b)
c1_o(ks + 1) = a_b*b1
d1 = h_sd(ks)*b1
u_o(ks) = b1*h_sd(ks)*u0(ks)
v_o(ks) = b1*h_sd(ks)*v0(ks)
t_o(ks) = b1*h_sd(ks)*t0(ks)
s_o(ks) = b1*h_sd(ks)*s0(ks)
do k = ks + 1, ke - 1
a_a = a_b
a_b = dt_now*kappa_ps(k + 1)*idz_int_s(k + 1)
bd1 = h_sd(k) + d1*a_a
b1 = 1.0_wp/(bd1 + a_b)
c1_o(k + 1) = a_b*b1
d1 = bd1*b1
u_o(k) = b1*(h_sd(k)*u0(k) + a_a*u_o(k - 1))
v_o(k) = b1*(h_sd(k)*v0(k) + a_a*v_o(k - 1))
t_o(k) = b1*(h_sd(k)*t0(k) + a_a*t_o(k - 1))
s_o(k) = b1*(h_sd(k)*s0(k) + a_a*s_o(k - 1))
end do
! Bottom layer of the band (ke).
a_a = a_b
if (ke > ks) then
b1 = 1.0_wp/(h_sd(ke) + d1*a_a)
t_o(ke) = b1*(h_sd(ke)*t0(ke) + a_a*t_o(ke - 1))
s_o(ke) = b1*(h_sd(ke)*s0(ke) + a_a*s_o(ke - 1))
if (ke == nz) then
b1nz = 1.0_wp/((h_sd(ke) + d1*a_a) + &
dt_now*kappa_ps(nz + 1)*idz_int_s(nz + 1))
else
b1nz = b1
end if
u_o(ke) = b1nz*(h_sd(ke)*u0(ke) + a_a*u_o(ke - 1))
v_o(ke) = b1nz*(h_sd(ke)*v0(ke) + a_a*v_o(ke - 1))
else
b1 = 1.0_wp/(h_sd(ke) + a_a)
t_o(ke) = b1*h_sd(ke)*t0(ke)
s_o(ke) = b1*h_sd(ke)*s0(ke)
if (ke == nz) then
b1nz = 1.0_wp/(h_sd(ke) + a_a + &
dt_now*kappa_ps(nz + 1)*idz_int_s(nz + 1))
else
b1nz = b1
end if
u_o(ke) = b1nz*h_sd(ke)*u0(ke)
v_o(ke) = b1nz*h_sd(ke)*v0(ke)
end if
if (abs(u_o(ke)) < vel_underflow) u_o(ke) = 0.0_wp
if (abs(v_o(ke)) < vel_underflow) v_o(ke) = 0.0_wp
! Back-substitution upward.
do k = ke - 1, ks, -1
u_o(k) = u_o(k) + c1_o(k + 1)*u_o(k + 1)
v_o(k) = v_o(k) + c1_o(k + 1)*v_o(k + 1)
t_o(k) = t_o(k) + c1_o(k + 1)*t_o(k + 1)
s_o(k) = s_o(k) + c1_o(k + 1)*s_o(k + 1)
if (abs(u_o(k)) < vel_underflow) u_o(k) = 0.0_wp
if (abs(v_o(k)) < vel_underflow) v_o(k) = 0.0_wp
end do
else
do k = 1, nz
if (abs(u_o(k)) < vel_underflow) u_o(k) = 0.0_wp
if (abs(v_o(k)) < vel_underflow) v_o(k) = 0.0_wp
end do
end if
! N^2, S^2 — mixed inside the band, original values outside.
n2_o(1) = 0.0_wp
n2_o(nz + 1) = 0.0_wp
s2_o(1) = 0.0_wp
s2_o(nz + 1) = 0.0_wp
do kk = 2, nz
if (kk - 1 >= ks .and. kk - 1 <= ke) then
ua = u_o(kk - 1)
va = v_o(kk - 1)
ta = t_o(kk - 1)
sa = s_o(kk - 1)
else
ua = u0(kk - 1)
va = v0(kk - 1)
ta = t0(kk - 1)
sa = s0(kk - 1)
end if
if (kk >= ks .and. kk <= ke) then
ub = u_o(kk)
vb = v_o(kk)
tb = t_o(kk)
sb = s_o(kk)
else
ub = u0(kk)
vb = v0(kk)
tb = t0(kk)
sb = s0(kk)
end if
n2v = idz_int_s(kk)*(dbuoy_t(kk)*(ta - tb) + dbuoy_s(kk)*(sa - sb))
if (n2v < 0.0_wp) n2v = 0.0_wp
n2_o(kk) = n2v
s2_o(kk) = ((ua - ub)**2 + (va - vb)**2)*idz_int_s(kk)**2
end do
end subroutine ks_projected_state