pure function ks_adaptive_dt(nz, dt_rem, itt_outer, max_substep_it, &
ri_crit, shearmix_rate, fri_curvature, &
src_max_chg, tol_err, vel_underflow, &
dbuoy_t, dbuoy_s, &
h_s, u_cur, v_cur, t_cur, s_cur, &
kappa_out_s, kappa_src_s, local_src_s, &
local_src_avg_s, ks_kap, ke_kap, &
idz_int_s) result(dt_now_r)
!! Largest dt_test <= dt_rem such that, after mixing for
!! 0.5*dt_test with `kappa_out_s`, the regenerated source stays
!! within the tolerance bands of the accepted-state source
!! (design doc section 5.3): a halving pass followed by a 5-step
!! refinement pass.
!$acc routine seq
integer, intent(in) :: nz, itt_outer, max_substep_it, ks_kap, ke_kap
real(wp), intent(in) :: dt_rem, ri_crit, shearmix_rate, fri_curvature
real(wp), intent(in) :: src_max_chg, tol_err, vel_underflow
real(wp), intent(in) :: dbuoy_t(NZLI), dbuoy_s(NZLI)
real(wp), intent(in) :: h_s(NZL), u_cur(NZL), v_cur(NZL)
real(wp), intent(in) :: t_cur(NZL), s_cur(NZL)
real(wp), intent(in) :: kappa_out_s(NZLI), kappa_src_s(NZLI)
real(wp), intent(in) :: local_src_s(NZLI), local_src_avg_s(NZLI)
real(wp), intent(in) :: idz_int_s(NZLI)
real(wp) :: dt_now_r
real(wp) :: tol_max(NZLI), tol_min_a(NZLI), tol_chg_a(NZLI)
real(wp) :: u_pr(NZL), v_pr(NZL), t_pr(NZL), s_pr(NZL)
real(wp) :: c1_pr(NZLI), n2_pr(NZLI), s2_pr(NZLI)
integer :: kk, k_lo, k_hi, ih, ir, max_halvings, ks_lyr, ke_lyr
real(wp) :: dt_tst, dt_inc, dt_try, idtt, tol_dksrc_low
real(wp) :: ksrc_tst, upper_t, lower_t, tol2_l
logical :: valid_dt, valid_try
if (src_max_chg == 10.0_wp) then
tol_dksrc_low = 0.95_wp
else
tol_dksrc_low = (src_max_chg - 0.5_wp)/src_max_chg
end if
tol2_l = 2.0_wp*tol_err
do kk = 1, nz + 1
tol_max(kk) = kappa_src_s(kk) + src_max_chg*local_src_s(kk)
tol_min_a(kk) = kappa_src_s(kk) - tol_dksrc_low*local_src_s(kk)
tol_chg_a(kk) = tol2_l*local_src_avg_s(kk)
end do
k_lo = max(ks_kap - 1, 2)
k_hi = min(ke_kap + 1, nz)
ks_lyr = max(ks_kap - 1, 1)
ke_lyr = min(ke_kap, nz)
dt_tst = dt_rem
valid_dt = .false.
max_halvings = (max_substep_it + 1 - itt_outer)/2
! Halving pass.
do ih = 1, max(max_halvings, 1)
call ks_projected_state(nz, 0.5_wp*dt_tst, ks_lyr, ke_lyr, &
vel_underflow, dbuoy_t, dbuoy_s, h_s, &
idz_int_s, u_cur, v_cur, t_cur, s_cur, &
kappa_out_s, u_pr, v_pr, t_pr, s_pr, &
c1_pr, n2_pr, s2_pr)
idtt = 0.0_wp
if (dt_tst > 0.0_wp) idtt = 1.0_wp/dt_tst
valid_dt = .true.
do kk = k_lo, k_hi
if (n2_pr(kk) < ri_crit*s2_pr(kk)) then
ksrc_tst = ks_src_func(ri_crit, shearmix_rate, fri_curvature, &
n2_pr(kk), s2_pr(kk))
upper_t = max(tol_max(kk), kappa_src_s(kk) + idtt*tol_chg_a(kk))
lower_t = min(tol_min_a(kk), kappa_src_s(kk) - idtt*tol_chg_a(kk))
if (ksrc_tst > upper_t .or. ksrc_tst < lower_t) then
valid_dt = .false.
exit
end if
else
lower_t = min(tol_min_a(kk), kappa_src_s(kk) - idtt*tol_chg_a(kk))
if (0.0_wp < lower_t) then
valid_dt = .false.
exit
end if
end if
end do
if (valid_dt) exit
dt_tst = 0.5_wp*dt_tst
end do
! Refinement pass.
dt_inc = 0.0_wp
if (dt_tst < dt_rem .and. valid_dt) then
dt_inc = 0.5_wp*dt_tst
do ir = 1, 5
dt_try = dt_tst + dt_inc
call ks_projected_state(nz, 0.5_wp*dt_try, ks_lyr, ke_lyr, &
vel_underflow, dbuoy_t, dbuoy_s, h_s, &
idz_int_s, u_cur, v_cur, t_cur, s_cur, &
kappa_out_s, u_pr, v_pr, t_pr, s_pr, &
c1_pr, n2_pr, s2_pr)
idtt = 0.0_wp
if (dt_try > 0.0_wp) idtt = 1.0_wp/dt_try
valid_try = .true.
do kk = k_lo, k_hi
if (n2_pr(kk) < ri_crit*s2_pr(kk)) then
ksrc_tst = ks_src_func(ri_crit, shearmix_rate, &
fri_curvature, n2_pr(kk), s2_pr(kk))
upper_t = max(tol_max(kk), kappa_src_s(kk) + idtt*tol_chg_a(kk))
lower_t = min(tol_min_a(kk), &
kappa_src_s(kk) - idtt*tol_chg_a(kk))
if (ksrc_tst > upper_t .or. ksrc_tst < lower_t) then
valid_try = .false.
exit
end if
else
lower_t = min(tol_min_a(kk), &
kappa_src_s(kk) - idtt*tol_chg_a(kk))
if (0.0_wp < lower_t) then
valid_try = .false.
exit
end if
end if
end do
if (valid_try) dt_tst = dt_try
dt_inc = 0.5_wp*dt_inc
end do
end if
dt_now_r = min(dt_tst*(1.0_wp + tol_err) + dt_inc, dt_rem)
end function ks_adaptive_dt