subroutine addsource_compact(qdt, dtfactor, qtC, wCTprim1, wCTprim2, wCTprim3, qt, wnew, x, dx, &
qsourcesplit)
use mod_global_parameters, only: dt, ndir, smalldouble, courantpar
real(dp), intent(in) :: qdt, dtfactor, qtC, qt
real(dp), intent(in) :: wCTprim1(nw_phys,3),wCTprim2(nw_phys,3),wCTprim3(nw_phys,3)
real(dp), intent(in) :: x(1:ndim), dx(1:ndim)
real(dp), intent(inout) :: wnew(nw_phys)
logical, intent(in) :: qsourcesplit
! .. local ..
real(dp) :: laplb_cd2
real(dp) :: Jdir1,Jdir2,Jdir3
integer :: idir
real(dp) :: Te_c, rho_c, pth_c
real(dp) :: gradT(3), Bmag2, Bmag, Bmag2_safe, bgradT
real(dp) :: bhat(3), gradT_perp(3), gradTperp_mag
real(dp) :: chi, sig_par, sig_perp
real(dp) :: q_sat, f_sat_par, f_sat_perp, tau_par, tau_perp
real(dp) :: cf2_tc, cmax2_tc
real(dp) :: Qpar_proj, Qperp_proj_k
integer :: k_tc
if (.not. qsourcesplit) then
!---------------------------------
! unsplit sources
!---------------------------------
! using the compact stencil formulation and adopting constant eta
! > eta*(laplacian of B)_idir added to B_idir
! > B_idir*[eta*(laplacian of B)_idir] added to e_
do idir = 1,ndim
laplb_cd2 = &
( &
wCTprim1(iw_mag(1)-1+idir,3) &
- 2*wCTprim1(iw_mag(1)-1+idir,2) &
+ wCTprim1(iw_mag(1)-1+idir,1) &
) &
/ dx(1)**2 + &
( &
wCTprim2(iw_mag(1)-1+idir,3) &
- 2*wCTprim2(iw_mag(1)-1+idir,2) &
+ wCTprim2(iw_mag(1)-1+idir,1) &
) &
/ dx(2)**2 + &
( &
wCTprim3(iw_mag(1)-1+idir,3) &
- 2*wCTprim3(iw_mag(1)-1+idir,2) &
+ wCTprim3(iw_mag(1)-1+idir,1) &
) &
/ dx(3)**2
wnew(iw_mag(1)-1+idir) = wnew(iw_mag(1)-1+idir) + qdt*mhd_eta*laplb_cd2
wnew(iw_e) = wnew(iw_e) + qdt * mhd_eta * laplb_cd2 * wCTprim1(iw_mag(1)-1+idir,2)
enddo
! > eta* J**2 added to e
Jdir1 = (wCTprim2(iw_mag(3),3) - wCTprim2(iw_mag(3),1)) &
/ 2.0_dp/dx(2) &
- (wCTprim3(iw_mag(2),3) - wCTprim3(iw_mag(2),1)) &
/ 2.0_dp/dx(3)
Jdir2 = (wCTprim3(iw_mag(1),3) - wCTprim3(iw_mag(1),1)) &
/2.0_dp/dx(3) &
- (wCTprim1(iw_mag(3),3) - wCTprim1(iw_mag(3),1)) &
/ 2.0_dp/dx(1)
Jdir3 = (wCTprim1(iw_mag(2),3) - wCTprim1(iw_mag(2),1)) &
/2.0_dp/dx(1) &
- (wCTprim2(iw_mag(1),3) - wCTprim2(iw_mag(1),1)) &
/ 2.0_dp/dx(2)
wnew(iw_e) = wnew(iw_e) + qdt*mhd_eta*(Jdir1**2+Jdir2**2+Jdir3**2)
! ---- 1. Local thermodynamic state ----
Te_c = wCTprim1(iw_e,2) / wCTprim1(iw_rho,2)
rho_c = wCTprim1(iw_rho,2)
pth_c = wCTprim1(iw_e,2)
! ---- 2. Temperature gradient ----
gradT(1) = (wCTprim1(iw_e,3)/wCTprim1(iw_rho,3) - wCTprim1(iw_e,1)/wCTprim1(iw_rho,1)) / (2.d0*dx(1))
gradT(2) = (wCTprim2(iw_e,3)/wCTprim2(iw_rho,3) - wCTprim2(iw_e,1)/wCTprim2(iw_rho,1)) / (2.d0*dx(2))
gradT(3) = (wCTprim3(iw_e,3)/wCTprim3(iw_rho,3) - wCTprim3(iw_e,1)/wCTprim3(iw_rho,1)) / (2.d0*dx(3))
!---------------------------------
! Anisotropic TC code here
!---------------------------------
! ---- 3. Field direction ----
Bmag2 = 0.0d0
do k_tc = 1, ndir
Bmag2 = Bmag2 + wCTprim1(iw_mag(k_tc),2)**2
end do
Bmag = sqrt(Bmag2)
Bmag2_safe = max(Bmag2, smalldouble**2)
do k_tc = 1, ndir
bhat(k_tc) = wCTprim1(iw_mag(k_tc),2) * Bmag / Bmag2_safe
end do
! ---- 4. Decompose gradT into parallel scalar + perpendicular vector ----
bgradT = bhat(1)*gradT(1) + bhat(2)*gradT(2) + bhat(3)*gradT(3)
do k_tc = 1, ndir
gradT_perp(k_tc) = gradT(k_tc) - bgradT * bhat(k_tc)
end do
gradTperp_mag = sqrt(gradT_perp(1)**2 + gradT_perp(2)**2 + gradT_perp(3)**2)
! ---- 5. Conductivities ----
if (tc_kappa_par > 0.0d0) then
sig_par = tc_kappa_par
else
sig_par = tc_kappa0_par * Te_c**2.5d0
end if
if (tc_kappa_perp > 0.0d0) then
sig_perp = tc_kappa_perp
else
chi = htc_Cchi * Bmag * Te_c**1.5d0 / rho_c
sig_perp = sig_par / (1.0d0 + chi**2)
end if
! ---- 6. Saturation and relaxation times (per channel) ----
cf2_tc = (Bmag2 + mhd_gamma*pth_c) / rho_c
cmax2_tc = 0.0d0
do k_tc = 1, ndim
cmax2_tc = max(cmax2_tc, 0.5d0*(cf2_tc + sqrt(abs( &
cf2_tc**2 - 4.0d0*mhd_gamma*pth_c*wCTprim1(iw_mag(k_tc),2)**2/rho_c**2))))
end do
q_sat = 1.5d0 * rho_c * (pth_c / rho_c)**1.5d0
f_sat_par = 1.0d0 / (1.0d0 + abs(sig_par * bgradT ) / q_sat)
f_sat_perp = 1.0d0 / (1.0d0 + abs(sig_perp * gradTperp_mag) / q_sat)
tau_par = max(4.d0*dt, f_sat_par *sig_par *Te_c*courantpar**2*(mhd_gamma-1.0d0) / (pth_c*cmax2_tc))
tau_perp = max(4.d0*dt, f_sat_perp*sig_perp*Te_c*courantpar**2*(mhd_gamma-1.0d0) / (pth_c*cmax2_tc))
! ---- 7. Decompose current q onto current field direction ----
Qpar_proj = wCTprim1(q_(1),2)*bhat(1) + wCTprim1(q_(2),2)*bhat(2) + wCTprim1(q_(3),2)*bhat(3)
! ---- 8. Relax parallel and perpendicular parts, each at its own rate ----
do k_tc = 1, ndir
Qperp_proj_k = wCTprim1(q_(k_tc),2) - Qpar_proj * bhat(k_tc)
wnew(q_(k_tc)) = wnew(q_(k_tc)) &
- qdt*(f_sat_par *sig_par *bgradT + Qpar_proj)*bhat(k_tc)/tau_par &
- qdt*(f_sat_perp*sig_perp*gradT_perp(k_tc) + Qperp_proj_k)/tau_perp
end do
!---------------------------------
! Field-aligned TC code here
!---------------------------------
else
!---------------------------------
! split sources
!---------------------------------
! Not yet implemented
end if
end subroutine addsource_compact