pure subroutine ks_find_kappa_tke(nz, tke_min, f2_val, &
ri_crit, shearmix_rate, fri_curvature, &
c_n2, c_s2, ilambda2, kappa_0, &
kappa_trunc, tke_bg, tol_err, max_inner_it, &
n2_in, s2_in, kappa_seed, k_q_io, &
idz_s, hint_s, il2_s, e1_s, &
tke_o, kappa_o, &
ksrc_sc, tkedec_sc, aq_sc, dq_sc, cq_sc, &
dk_sc, ck_sc, ild2_sc)
!! The inner Picard solve (design doc section 5.4): alternate a TKE
!! tridiagonal sweep (Dirichlet surface, e1 tail below the deepest
!! active interface) with a kappa tridiagonal sweep (smooth
!! truncation ramp + active-range tracking) until the Picard
!! increment converges. Scratch arrays are supplied by the caller
!! to avoid double-allocating per-thread stack.
!$acc routine seq
integer, intent(in) :: nz, max_inner_it
real(wp), intent(in) :: tke_min, f2_val
real(wp), intent(in) :: ri_crit, shearmix_rate, fri_curvature
real(wp), intent(in) :: c_n2, c_s2, ilambda2, kappa_0
real(wp), intent(in) :: kappa_trunc, tke_bg, tol_err
real(wp), intent(in) :: n2_in(NZLI), s2_in(NZLI)
real(wp), intent(in) :: kappa_seed(NZLI)
real(wp), intent(inout) :: k_q_io(NZLI)
real(wp), intent(in) :: idz_s(NZL), hint_s(NZLI), il2_s(NZLI), e1_s(NZLI)
real(wp), intent(out) :: tke_o(NZLI), kappa_o(NZLI)
real(wp), intent(inout) :: ksrc_sc(NZLI), tkedec_sc(NZLI)
real(wp), intent(inout) :: aq_sc(NZL), dq_sc(NZLI), cq_sc(NZLI)
real(wp), intent(inout) :: dk_sc(NZLI), ck_sc(NZLI), ild2_sc(NZLI)
integer :: kk, k, k2, it
integer :: ks_src, ke_src, ks_kap, ke_kap, ks_kp, ke_kp, ke_tke
integer :: ks_new, ke_new, k_lo, k_hi
real(wp) :: cqc_l, ckc_l, bqd1_l, bq_l, tsrc_l, raw_l, dnom_l
real(wp) :: bkd1_l, bk_l, trv, tr2v, lhs_l, rhs_l
logical :: conv
ks_src = nz + 2
ke_src = 0
ksrc_sc(1) = 0.0_wp
ksrc_sc(nz + 1) = 0.0_wp
do kk = 2, nz
ksrc_sc(kk) = ks_src_func(ri_crit, shearmix_rate, fri_curvature, &
n2_in(kk), s2_in(kk))
if (ksrc_sc(kk) > 0.0_wp) then
if (ks_src > kk) ks_src = kk
ke_src = kk
end if
end do
if (ks_src > ke_src) then
do kk = 1, nz + 1
tke_o(kk) = tke_min
kappa_o(kk) = 0.0_wp
k_q_io(kk) = 0.0_wp
end do
return
end if
do kk = 2, nz
tkedec_sc(kk) = sqrt(c_n2*n2_in(kk) + c_s2*s2_in(kk))
end do
tke_o(1) = tke_bg
do kk = 2, nz
if (kappa_seed(kk) > 0.0_wp .and. k_q_io(kk) > 0.0_wp) then
tke_o(kk) = kappa_seed(kk)/k_q_io(kk)
else
tke_o(kk) = tke_min
end if
end do
tke_o(nz + 1) = tke_min
do kk = 1, nz + 1
kappa_o(kk) = kappa_seed(kk)
end do
kappa_o(1) = 0.0_wp
kappa_o(nz + 1) = 0.0_wp
ks_kap = 2
ke_kap = nz
ks_kp = 2
ke_kp = nz
do it = 1, max_inner_it
! (a) TKE tridiagonal sweep. ke_tke is an INTERFACE index
! (1..nz+1; bed = nz+1), so the clamp is nz+1 and the bed
! Dirichlet branch fires at ke_tke == nz+1 (design doc 5.4a,
! ambiguity #6 — "0-based nz" = the bed interface = local nz+1).
ke_tke = min(max(ke_kap, ke_kp) + 1, nz + 1)
do k = 1, min(ke_tke, nz)
aq_sc(k) = (0.5_wp*(kappa_o(k) + kappa_o(k + 1)) + kappa_0)*idz_s(k)
end do
dq_sc(1) = -tke_o(1)
tke_o(1) = tke_bg
cq_sc(2) = 0.0_wp
cqc_l = 1.0_wp
do kk = 2, ke_tke - 1
dq_sc(kk) = -tke_o(kk)
tsrc_l = (kappa_o(kk) + kappa_0)*s2_in(kk) + tke_bg*tkedec_sc(kk)
bqd1_l = hint_s(kk)*(tkedec_sc(kk) + n2_in(kk)*k_q_io(kk)) + &
cqc_l*aq_sc(kk - 1)
bq_l = 1.0_wp/(bqd1_l + aq_sc(kk))
tke_o(kk) = bq_l*(hint_s(kk)*tsrc_l + aq_sc(kk - 1)*tke_o(kk - 1))
cq_sc(kk + 1) = aq_sc(kk)*bq_l
cqc_l = bqd1_l*bq_l
end do
if (ke_tke == nz + 1) then
tke_o(nz + 1) = tke_min
dq_sc(nz + 1) = 0.0_wp
else
kk = ke_tke
tsrc_l = kappa_0*s2_in(kk) + tke_bg*tkedec_sc(kk)
bq_l = 1.0_wp/(hint_s(kk)*tkedec_sc(kk) + cqc_l*aq_sc(kk - 1) + &
aq_sc(kk))
cq_sc(kk + 1) = aq_sc(kk)*bq_l
dq_sc(kk) = -tke_o(kk)
raw_l = bq_l*(hint_s(kk)*tsrc_l + aq_sc(kk - 1)*tke_o(kk - 1))
dnom_l = 1.0_wp - cq_sc(kk + 1)*e1_s(kk + 1)
if (abs(dnom_l) > 1.0e-30_wp) then
tke_o(kk) = max((raw_l + cq_sc(kk + 1)* &
(tke_o(kk + 1) - e1_s(kk + 1)*tke_o(kk)))/dnom_l, &
tke_min)
else
tke_o(kk) = max(raw_l, tke_min)
end if
dq_sc(kk) = tke_o(kk) + dq_sc(kk)
do k2 = ke_tke + 1, nz + 1
dq_sc(k2) = e1_s(k2)*dq_sc(k2 - 1)
tke_o(k2) = max(tke_o(k2) + dq_sc(k2), tke_min)
if (abs(dq_sc(k2)) < 1.0e-16_wp*tke_o(k2)) exit
end do
end if
do kk = ke_tke - 1, 1, -1
tke_o(kk) = max(tke_o(kk) + cq_sc(kk + 1)*tke_o(kk + 1), tke_min)
dq_sc(kk) = tke_o(kk) + dq_sc(kk)
end do
! (b) kappa tridiagonal sweep with truncation ramp + range track.
ks_kp = ks_kap
ke_kp = ke_kap
do kk = 2, nz
if (tke_o(kk) > 0.0_wp) then
ild2_sc(kk) = (n2_in(kk)*ilambda2 + f2_val)/tke_o(kk) + il2_s(kk)
else
ild2_sc(kk) = 1.0e30_wp
end if
end do
dk_sc(1) = 0.0_wp
ck_sc(2) = 0.0_wp
ckc_l = 1.0_wp
ke_new = 0
ks_new = nz
do kk = 2, nz
dk_sc(kk) = -kappa_o(kk)
bkd1_l = hint_s(kk)*ild2_sc(kk) + ckc_l*idz_s(kk - 1)
bk_l = 1.0_wp/(bkd1_l + idz_s(kk))
kappa_o(kk) = bk_l*(idz_s(kk - 1)*kappa_o(kk - 1) + &
hint_s(kk)*ksrc_sc(kk))
ck_sc(kk + 1) = idz_s(kk)*bk_l
ckc_l = bkd1_l*bk_l
trv = ckc_l*kappa_trunc
tr2v = 2.0_wp*trv
if (kappa_o(kk) < trv) then
kappa_o(kk) = 0.0_wp
if (kk > ke_src) then
ke_kap = kk - 1
k_q_io(kk) = 0.0_wp
exit
end if
else if (kappa_o(kk) < tr2v) then
kappa_o(kk) = 2.0_wp*(kappa_o(kk) - trv)
end if
ke_new = kk
end do
if (ke_new > 0) ke_kap = ke_new
if (ke_kap >= 1 .and. tke_o(ke_kap) > 0.0_wp) then
k_q_io(ke_kap) = kappa_o(ke_kap)/tke_o(ke_kap)
end if
dk_sc(ke_kap) = dk_sc(ke_kap) + kappa_o(ke_kap)
do kk = ke_kap + 2, ke_kp + 1
if (kk > nz) exit
dk_sc(kk) = -kappa_o(kk)
kappa_o(kk) = 0.0_wp
k_q_io(kk) = 0.0_wp
end do
ks_new = 2
do kk = ke_kap - 1, 2, -1
kappa_o(kk) = kappa_o(kk) + ck_sc(kk + 1)*kappa_o(kk + 1)
if (kappa_o(kk) <= kappa_trunc) then
kappa_o(kk) = 0.0_wp
if (kk < ks_src) then
ks_kap = kk + 1
k_q_io(kk) = 0.0_wp
exit
end if
else if (kappa_o(kk) < 2.0_wp*kappa_trunc) then
kappa_o(kk) = 2.0_wp*(kappa_o(kk) - kappa_trunc)
end if
dk_sc(kk) = dk_sc(kk) + kappa_o(kk)
if (tke_o(kk) > 0.0_wp) then
k_q_io(kk) = kappa_o(kk)/tke_o(kk)
else
k_q_io(kk) = 0.0_wp
end if
ks_new = kk
end do
ks_kap = max(ks_new, 2)
do kk = ks_kp, ks_kap - 2
if (kk < 2) cycle
kappa_o(kk) = 0.0_wp
k_q_io(kk) = 0.0_wp
end do
! (c) Picard convergence test.
k_lo = min(ks_kap, ks_kp)
k_hi = max(ke_kap, ke_kp)
conv = .true.
do kk = k_lo, k_hi
lhs_l = abs(dk_sc(kk))
rhs_l = tol_err*(kappa_0 + kappa_o(kk) - 0.5_wp*dk_sc(kk))
if (lhs_l > rhs_l) then
conv = .false.
exit
end if
end do
if (conv) exit
end do
kappa_o(1) = 0.0_wp
kappa_o(nz + 1) = 0.0_wp
end subroutine ks_find_kappa_tke