diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index ce5cedd1159dcae0a9e6b31e351ab085fd740881..edd297b22578b08c921d1cbd39c9f77efcd1b944 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -54,105 +54,34 @@ END SUBROUTINE build_dnjs_table
!******************************************************************************!
SUBROUTINE evaluate_kernels
USE basic
- USE array, Only : kernel_e, kernel_i, gxy, gyy, gxx
+ USE array, Only : kernel_e, kernel_i, kparray
USE grid
USE model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, CLOS, sigmae2_taue_o2, sigmai2_taui_o2
IMPLICIT NONE
+ INTEGER :: j_int
+ REAL(dp) :: j_dp, y_, kp2_, kx_, ky_
- REAL(dp) :: factj, j_dp, j_int
- REAL(dp) :: be_2, bi_2, alphaD
- REAL(dp) :: kx, ky, kperp2
-
+DO ikx = ikxs,ikxe
+DO iky = ikys,ikye
+DO iz = izs,ize
!!!!! Electron kernels !!!!!
- !Precompute species dependant factors
- factj = 1.0 ! Start of the recursive factorial
- DO ij = 1, jmaxe+1
+ DO ij = ijsg_e, ijeg_e
j_int = jarray_e(ij)
- j_dp = REAL(j_int,dp) ! REAL of degree
-
- ! Recursive factorial
- IF (j_dp .GT. 0) THEN
- factj = factj * j_dp
- ELSE
- factj = 1._dp
- ENDIF
-
- DO ikx = ikxs,ikxe
- kx = kxarray(ikx)
- DO iky = ikys,ikye
- ky = kyarray(iky)
- DO iz = izs,ize
- kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
- be_2 = kperp2 * sigmae2_taue_o2
- kernel_e(ij, ikx, iky,iz) = be_2**j_int * exp(-be_2)/factj
- ENDDO
- ENDDO
- ENDDO
+ j_dp = REAL(j_int,dp)
+ y_ = sigmae2_taue_o2 * kparray(ikx,iky,iz)**2
+ kernel_e(ij,ikx,iky,iz) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
ENDDO
- ! Kernels closure
- DO ikx = ikxs,ikxe
- kx = kxarray(ikx)
- DO iky = ikys,ikye
- ky = kyarray(iky)
- DO iz = izs,ize
- kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
- be_2 = kperp2 * sigmae2_taue_o2
- ! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj (/!\ ij = j+1)
- kernel_e(ijeg_e,ikx,iky,iz) = be_2/(real(ijeg_e-1,dp))*kernel_e(ije_e,ikx,iky,iz)
- ! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0)
- IF ( be_2 .NE. 0 ) THEN
- kernel_e(ijsg_e,ikx,iky,iz) = (real(ijsg_e,dp))/be_2*kernel_e(ijs_e,ikx,iky,iz)
- ELSE
- kernel_e(ijsg_e,ikx,iky,iz) = 0._dp
- ENDIF
- ENDDO
- ENDDO
- ENDDO
-
!!!!! Ion kernels !!!!!
- factj = 1.0 ! Start of the recursive factorial
- DO ij = 1, jmaxi+1
- j_int = jarray_e(ij)
- j_dp = REAL(j_int,dp) ! REAL of degree
-
- ! Recursive factorial
- IF (j_dp .GT. 0) THEN
- factj = factj * j_dp
- ELSE
- factj = 1._dp
- ENDIF
-
- DO ikx = ikxs,ikxe
- kx = kxarray(ikx)
- DO iky = ikys,ikye
- ky = kyarray(iky)
- DO iz = izs,ize
- kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
- bi_2 = kperp2 * sigmai2_taui_o2
- kernel_i(ij, ikx, iky,iz) = bi_2**j_int * exp(-bi_2)/factj
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ! Kernels closure
- DO ikx = ikxs,ikxe
- kx = kxarray(ikx)
- DO iky = ikys,ikye
- ky = kyarray(iky)
- DO iz = izs,ize
- kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
- bi_2 = kperp2 * sigmai2_taui_o2
- ! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj
- kernel_i(ijeg_i,ikx,iky,iz) = bi_2/(real(ijeg_i-1,dp))*kernel_i(ije_i,ikx,iky,iz)
- ! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0)
- IF ( bi_2 .NE. 0 ) THEN
- kernel_i(ijsg_i,ikx,iky,iz) = (real(ijsg_i,dp))/bi_2*kernel_i(ijs_i,ikx,iky,iz)
- ELSE
- kernel_i(ijsg_i,ikx,iky,iz) = 0._dp
- ENDIF
- ENDDO
- ENDDO
+ DO ij = ijsg_i, ijeg_i
+ j_int = jarray_i(ij)
+ j_dp = REAL(j_int,dp)
+ y_ = sigmai2_taui_o2 * kparray(ikx,iky,iz)**2
+ kernel_i(ij,ikx,iky,iz) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
ENDDO
+ENDDO
+ENDDO
+ENDDO
+
END SUBROUTINE evaluate_kernels
!******************************************************************************!
@@ -174,30 +103,35 @@ SUBROUTINE compute_lin_coeff
DO ij = ijs_e, ije_e
j_int= jarray_e(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
+ ! All Napj terms
xnepj(ip,ij) = taue_qe*(CurvB*(2._dp*p_dp + 1._dp) &
+GradB*(2._dp*j_dp + 1._dp))
- ynepp1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp+1._dp)
- ynepm1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp)
- ynepp1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp+1._dp)
- ynepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp)
- zNepm1j (ip,ij) = -SQRT(tau_e)/sigma_e * (2._dp*j_dp+1_dp)*SQRT(p_dp)
- zNepm1jp1(ip,ij) = +SQRT(tau_e)/sigma_e * (j_dp+1_dp)*SQRT(p_dp)
- zNepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp)
+ ! Mirror force terms
+ ynepp1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp+1._dp)
+ ynepm1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp)
+ ynepp1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp+1._dp)
+ ynepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp)
+ zNepm1j (ip,ij) = +SQRT(tau_e)/sigma_e * (2._dp*j_dp+1_dp)*SQRT(p_dp)
+ zNepm1jp1(ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1_dp)*SQRT(p_dp)
+ zNepm1jm1(ip,ij) = -SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp)
ENDDO
ENDDO
DO ip = ips_e, ipe_e
p_int= parray_e(ip) ! Hermite degree
p_dp = REAL(p_int,dp) ! REAL of Hermite degree
- xnepp1j(ip) = sqrtTaue_qe * SQRT(p_dp + 1_dp)
- xnepm1j(ip) = sqrtTaue_qe * SQRT(p_dp)
- xnepp2j(ip) = taue_qe * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
- xnepm2j(ip) = taue_qe * SQRT(p_dp * (p_dp - 1._dp))
+ ! Landau damping coefficients (ddz napj term)
+ xnepp1j(ip) = SQRT(tau_e)/sigma_e * SQRT(p_dp + 1_dp)
+ xnepm1j(ip) = SQRT(tau_e)/sigma_e * SQRT(p_dp)
+ ! Magnetic curvature term
+ xnepp2j(ip) = taue_qe * CurvB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
+ xnepm2j(ip) = taue_qe * CurvB * SQRT(p_dp * (p_dp - 1._dp))
ENDDO
DO ij = ijs_e, ije_e
j_int= jarray_e(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
- xnepjp1(ij) = -taue_qe * (j_dp + 1._dp)
- xnepjm1(ij) = -taue_qe * j_dp
+ ! Magnetic gradient term
+ xnepjp1(ij) = -taue_qe * GradB * (j_dp + 1._dp)
+ xnepjm1(ij) = -taue_qe * GradB * j_dp
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Ions linear coefficients for moment RHS !!!!!!!!!!
@@ -207,30 +141,36 @@ SUBROUTINE compute_lin_coeff
DO ij = ijs_i, ije_i
j_int= jarray_i(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
- xnepj(ip,ij) = taui_qi*(CurvB*(2._dp*p_dp + 1._dp) &
+ ! All Napj terms
+ xnipj(ip,ij) = taui_qi*(CurvB*(2._dp*p_dp + 1._dp) &
+GradB*(2._dp*j_dp + 1._dp))
- ynipp1j (ip,ij) = -SQRT(tau_i)/sigma_i * (j_dp+1)*SQRT(p_dp+1._dp)
- ynipm1j (ip,ij) = -SQRT(tau_i)/sigma_i * (j_dp+1)*SQRT(p_dp)
- ynipp1jm1(ip,ij) = +SQRT(tau_i)/sigma_i * j_dp*SQRT(p_dp+1._dp)
- ynipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i * j_dp*SQRT(p_dp)
- zNipm1j (ip,ij) = -SQRT(tau_i)/sigma_i * (2._dp*j_dp+1_dp)*SQRT(p_dp)
- zNipm1jp1(ip,ij) = +SQRT(tau_i)/sigma_i * (j_dp+1_dp)*SQRT(p_dp)
- zNipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i * j_dp*SQRT(p_dp)
+ ! Mirror force terms
+ ynipp1j (ip,ij) = -SQRT(tau_i)/sigma_i* (j_dp+1)*SQRT(p_dp+1._dp)
+ ynipm1j (ip,ij) = -SQRT(tau_i)/sigma_i* (j_dp+1)*SQRT(p_dp)
+ ynipp1jm1(ip,ij) = +SQRT(tau_i)/sigma_i* j_dp*SQRT(p_dp+1._dp)
+ ynipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i* j_dp*SQRT(p_dp)
+ ! Trapping terms
+ zNipm1j (ip,ij) = +SQRT(tau_i)/sigma_i* (2._dp*j_dp+1_dp)*SQRT(p_dp)
+ zNipm1jp1(ip,ij) = -SQRT(tau_i)/sigma_i* (j_dp+1_dp)*SQRT(p_dp)
+ zNipm1jm1(ip,ij) = -SQRT(tau_i)/sigma_i* j_dp*SQRT(p_dp)
ENDDO
ENDDO
DO ip = ips_i, ipe_i
p_int= parray_i(ip) ! Hermite degree
p_dp = REAL(p_int,dp) ! REAL of Hermite degree
- xnipp1j(ip) = sqrtTaui_qi * SQRT(p_dp + 1._dp)
- xnipm1j(ip) = sqrtTaui_qi * SQRT(p_dp)
- xnipp2j(ip) = taui_qi * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
- xnipm2j(ip) = taui_qi * SQRT(p_dp * (p_dp - 1._dp))
+ ! Landau damping coefficients (ddz napj term)
+ xnipp1j(ip) = SQRT(tau_i)/sigma_i * SQRT(p_dp + 1._dp)
+ xnipm1j(ip) = SQRT(tau_i)/sigma_i * SQRT(p_dp)
+ ! Magnetic curvature term
+ xnipp2j(ip) = taui_qi * CurvB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
+ xnipm2j(ip) = taui_qi * CurvB * SQRT(p_dp * (p_dp - 1._dp))
ENDDO
DO ij = ijs_i, ije_i
j_int= jarray_i(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
- xnipjp1(ij) = -taui_qi * (j_dp + 1._dp)
- xnipjm1(ij) = -taui_qi * j_dp
+ ! Magnetic gradient term
+ xnipjp1(ij) = -taui_qi * GradB * (j_dp + 1._dp)
+ xnipjm1(ij) = -taui_qi * GradB * j_dp
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! ES linear coefficients for moment RHS !!!!!!!!!!