diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 3b3c5fdee72e727b36bd684b1ad9a30ee66b764a..a93026b14c35d54e5af32ba3f553c5c143fdf308 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -1,4 +1,10 @@
-SUBROUTINE moments_eq_rhs
+!_____________________________________________________________________________!
+!_____________________________________________________________________________!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!! Electrons moments RHS !!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!_____________________________________________________________________________!
+SUBROUTINE moments_eq_rhs_e
USE basic
USE time_integration
@@ -20,14 +26,11 @@ SUBROUTINE moments_eq_rhs
REAL(dp) :: xCapj, xCa20, xCa01, xCa10 ! Coll. factors depending on the pj loop
COMPLEX(dp) :: TNapj, TNapp1j, TNapm1j, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi
COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs
- COMPLEX(dp) :: test_nan, i_kz
+ COMPLEX(dp) :: i_kz
! Measuring execution time
CALL cpu_time(t0_rhs)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !!!!!!!!! Electrons moments RHS !!!!!!!!!
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ploope : DO ip = ips_e, ipe_e ! loop over Hermite degree indices
p_int= parray_e(ip) ! Hermite polynom. degree
p_dp = REAL(p_int,dp) ! REAL of the Hermite degree
@@ -42,8 +45,6 @@ SUBROUTINE moments_eq_rhs
! N_e^{p-2,j} coeff
xNapm2j = taue_qe_etaB * SQRT(p_dp * (p_dp - 1._dp))
- ! factj = 1.0 ! Start of the recursive factorial
-
jloope : DO ij = ijs_e, ije_e ! loop over Laguerre degree indices
j_int= jarray_e(ij) ! Laguerre polynom. degree
j_dp = REAL(j_int,dp) ! REAL of degree
@@ -92,13 +93,6 @@ SUBROUTINE moments_eq_rhs
xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar= 0._dp
ENDIF
- ! Recursive factorial
- ! IF (j_dp .GT. 0) THEN
- ! factj = factj * j_dp
- ! ELSE
- ! factj = 1._dp
- ! ENDIF
-
! Loop on kspace
krloope : DO ikr = ikrs,ikre
kzloope : DO ikz = ikzs,ikze
@@ -113,69 +107,36 @@ SUBROUTINE moments_eq_rhs
!! Compute moments and mixing terms
! term propto N_e^{p,j}
TNapj = xNapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
+ TNapp1j = 0._dp; TNapm1j = 0._dp
+ TNapp2j = 0._dp; TNapm2j = 0._dp
+ TNapjp1 = 0._dp; TNapjm1 = 0._dp
! term propto N_e^{p+1,j} and kparallel
- IF ( (p_int+1 .LE. pmaxe) ) THEN ! OoB check
- TNapp1j = xNapp1j * moments_e(ip+1,ij,ikr,ikz,updatetlevel)
- ELSE
- TNapp1j = 0._dp
- ENDIF
+ IF ( p_int+1 .LE. pmaxe ) TNapp1j = xNapp1j * moments_e(ip+1,ij,ikr,ikz,updatetlevel)
! term propto N_e^{p-1,j} and kparallel
- IF ( (p_int-1 .GE. 0) ) THEN ! OoB check
- TNapm1j = xNapm1j * moments_e(ip-1,ij,ikr,ikz,updatetlevel)
- ELSE
- TNapm1j = 0._dp
- ENDIF
+ IF ( p_int-1 .GE. 0 ) TNapm1j = xNapm1j * moments_e(ip-1,ij,ikr,ikz,updatetlevel)
! term propto N_e^{p+2,j}
- IF (p_int+2 .LE. pmaxe) THEN ! OoB check
- TNapp2j = xNapp2j * moments_e(ip+2,ij,ikr,ikz,updatetlevel)
- ELSE
- TNapp2j = 0._dp
- ENDIF
+ IF ( p_int+2 .LE. pmaxe ) TNapp2j = xNapp2j * moments_e(ip+2,ij,ikr,ikz,updatetlevel)
! term propto N_e^{p-2,j}
- IF (p_int-2 .GE. 0) THEN ! OoB check
- TNapm2j = xNapm2j * moments_e(ip-2,ij,ikr,ikz,updatetlevel)
- ELSE
- TNapm2j = 0._dp
- ENDIF
+ IF ( p_int-2 .GE. 0 ) TNapm2j = xNapm2j * moments_e(ip-2,ij,ikr,ikz,updatetlevel)
! xterm propto N_e^{p,j+1}
- IF (j_int+1 .LE. jmaxe) THEN ! OoB check
- TNapjp1 = xNapjp1 * moments_e(ip,ij+1,ikr,ikz,updatetlevel)
- ELSE
- TNapjp1 = 0._dp
- ENDIF
+ IF ( j_int+1 .LE. jmaxe ) TNapjp1 = xNapjp1 * moments_e(ip,ij+1,ikr,ikz,updatetlevel)
! term propto N_e^{p,j-1}
- IF (j_int-1 .GE. 0) THEN ! OoB check
- TNapjm1 = xNapjm1 * moments_e(ip,ij-1,ikr,ikz,updatetlevel)
- ELSE
- TNapjm1 = 0._dp
- ENDIF
+ IF ( j_int-1 .GE. 0 ) TNapjm1 = xNapjm1 * moments_e(ip,ij-1,ikr,ikz,updatetlevel)
!! Collision
IF (CO .EQ. -3) THEN ! GK Dougherty
CALL DoughertyGK_e(ip,ij,ikr,ikz,TColl)
ELSEIF (CO .EQ. -2) THEN ! DK Dougherty
- IF ( (pmaxe .GE. 2) ) THEN ! OoB check
- TColl20 = xCa20 * moments_e(3,1,ikr,ikz,updatetlevel)
- ELSE
- TColl20 = 0._dp
- ENDIF
- IF ( (jmaxe .GE. 1) ) THEN ! OoB check
- TColl01 = xCa01 * moments_e(1,2,ikr,ikz,updatetlevel)
- ELSE
- TColl01 = 0._dp
- ENDIF
- IF ( (pmaxe .GE. 1) ) THEN ! OoB + odd number for Hermite degrees check
- TColl10 = xCa10 * moments_e(2,1,ikr,ikz,updatetlevel)
- ELSE
- TColl10 = 0._dp
- ENDIF
+ TColl20 = 0._dp; TColl01 = 0._dp; TColl10 = 0._dp
+ IF ( (pmaxe .GE. 2) ) TColl20 = xCa20 * moments_e(3,1,ikr,ikz,updatetlevel)
+ IF ( (jmaxe .GE. 1) ) TColl01 = xCa01 * moments_e(1,2,ikr,ikz,updatetlevel)
+ IF ( (pmaxe .GE. 1) ) TColl10 = xCa10 * moments_e(2,1,ikr,ikz,updatetlevel)
! Total collisional term
TColl = xCapj* moments_e(ip,ij,ikr,ikz,updatetlevel)&
+ TColl20 + TColl01 + TColl10
ELSEIF (CO .EQ. -1) THEN ! Full Coulomb for electrons (COSOlver matrix)
-
TColl = 0._dp ! Initialization
ploopee: DO ip2 = 1,pmaxe+1 ! sum the electron-self and electron-ion test terms
@@ -198,23 +159,15 @@ SUBROUTINE moments_eq_rhs
ELSEIF (CO .EQ. 0) THEN ! Lenard Bernstein
TColl = xCapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
+
ENDIF
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
kernelj = kernel_e(ij, ikr, ikz)
- IF ( j_int+1 .LE. jmaxe ) THEN
- kerneljp1 = kernel_e(ij+1, ikr, ikz)
- ELSE
- kerneljp1 = 0._dp
- ENDIF
-
- IF ( j_int-1 .GE. 0 ) THEN
- kerneljm1 = kernel_e(ij-1, ikr, ikz)
- ELSE
- kerneljm1 = 0._dp
- ENDIF
-
+ kerneljp1 = 0._dp; kerneljm1 = 0._dp
+ IF ( j_int+1 .LE. jmaxe ) kerneljp1 = kernel_e(ij+1, ikr, ikz)
+ IF ( j_int-1 .GE. 0 ) kerneljm1 = kernel_e(ij-1, ikr, ikz)
Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
ELSE
Tphi = 0._dp
@@ -239,12 +192,44 @@ SUBROUTINE moments_eq_rhs
END DO jloope
END DO ploope
- !_____________________________________________________________________________!
- !_____________________________________________________________________________!
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !!!!!!!!! Ions moments RHS !!!!!!!!!
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !_____________________________________________________________________________!
+ ! Execution time end
+ CALL cpu_time(t1_rhs)
+ tc_rhs = tc_rhs + (t1_rhs-t0_rhs)
+
+END SUBROUTINE moments_eq_rhs_e
+!_____________________________________________________________________________!
+!_____________________________________________________________________________!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!! Ions moments RHS !!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!_____________________________________________________________________________!
+SUBROUTINE moments_eq_rhs_i
+
+ USE basic
+ USE time_integration, ONLY: updatetlevel
+ USE array
+ USE fields
+ USE grid
+ USE model
+ USE prec_const
+ USE utility, ONLY : is_nan
+ USE collision
+ IMPLICIT NONE
+
+ INTEGER :: ip2, ij2, p_int, j_int, p2_int, j2_int ! loops indices and polynom. degrees
+ REAL(dp) :: p_dp, j_dp
+ REAL(dp) :: kr, kz, kperp2
+ REAL(dp) :: kernelj, kerneljp1, kerneljm1 ! Kernel functions and variable
+ REAL(dp) :: xNapj, xNapp1j, xNapm1j, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! Mom. factors depending on the pj loop
+ REAL(dp) :: xphij, xphijp1, xphijm1, xphijpar ! ESpot. factors depending on the pj loop
+ REAL(dp) :: xCapj, xCa20, xCa01, xCa10 ! Coll. factors depending on the pj loop
+ COMPLEX(dp) :: TNapj, TNapp1j, TNapm1j, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi
+ COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs
+ COMPLEX(dp) :: i_kz
+
+ ! Measuring execution time
+ CALL cpu_time(t0_rhs)
+
ploopi : DO ip = ips_i, ipe_i ! Hermite loop
p_int= parray_i(ip) ! Hermite degree
p_dp = REAL(p_int,dp) ! REAL of Hermite degree
@@ -253,14 +238,11 @@ SUBROUTINE moments_eq_rhs
xNapp1j = sqrtTaui_qi * SQRT(p_dp + 1)
! N_i^{p-1,j} coeff
xNapm1j = sqrtTaui_qi * SQRT(p_dp)
-
! x N_i^{p+2,j} coeff
xNapp2j = taui_qi_etaB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
! x N_i^{p-2,j} coeff
xNapm2j = taui_qi_etaB * SQRT(p_dp * (p_dp - 1._dp))
- ! factj = 1._dp ! Start of the recursive factorial
-
jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1
j_int= jarray_i(ij) ! REALof Laguerre degree
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
@@ -309,13 +291,6 @@ SUBROUTINE moments_eq_rhs
xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar = 0._dp
ENDIF
- ! ! Recursive factorial
- ! IF (j_dp .GT. 0) THEN
- ! factj = factj * j_dp
- ! ELSE
- ! factj = 1._dp
- ! ENDIF
-
! Loop on kspace
krloopi : DO ikr = ikrs,ikre
kzloopi : DO ikz = ikzs,ikze
@@ -328,63 +303,31 @@ SUBROUTINE moments_eq_rhs
!! Compute moments and mixing terms
! term propto N_i^{p,j}
TNapj = xNapj * moments_i(ip,ij,ikr,ikz,updatetlevel)
+ TNapp1j = 0._dp; TNapm1j = 0._dp
+ TNapp2j = 0._dp; TNapm2j = 0._dp
+ TNapjp1 = 0._dp; TNapjm1 = 0._dp
! term propto N_i^{p+1,j}
- IF ( (p_int+1 .LE. pmaxi) ) THEN ! OoB check
- TNapp1j = xNapp1j * moments_i(ip+1,ij,ikr,ikz,updatetlevel)
- ELSE
- TNapp1j = 0._dp
- ENDIF
+ IF ( p_int+1 .LE. pmaxi ) TNapp1j = xNapp1j * moments_i(ip+1,ij,ikr,ikz,updatetlevel)
! term propto N_i^{p-1,j}
- IF ( (p_int-1 .GE. 0) ) THEN ! OoB check
- TNapm1j = xNapm1j * moments_i(ip-1,ij,ikr,ikz,updatetlevel)
- ELSE
- TNapm1j = 0._dp
- ENDIF
+ IF ( p_int-1 .GE. 0 ) TNapm1j = xNapm1j * moments_i(ip-1,ij,ikr,ikz,updatetlevel)
! term propto N_i^{p+2,j}
- IF (p_int+2 .LE. pmaxi) THEN ! OoB check
- TNapp2j = xNapp2j * moments_i(ip+2,ij,ikr,ikz,updatetlevel)
- ELSE
- TNapp2j = 0._dp
- ENDIF
+ IF ( p_int+2 .LE. pmaxi ) TNapp2j = xNapp2j * moments_i(ip+2,ij,ikr,ikz,updatetlevel)
! term propto N_i^{p-2,j}
- IF (p_int-2 .GE. 0) THEN ! OoB check
- TNapm2j = xNapm2j * moments_i(ip-2,ij,ikr,ikz,updatetlevel)
- ELSE
- TNapm2j = 0._dp
- ENDIF
+ IF ( p_int-2 .GE. 0 ) TNapm2j = xNapm2j * moments_i(ip-2,ij,ikr,ikz,updatetlevel)
! xterm propto N_i^{p,j+1}
- IF (j_int+1 .LE. jmaxi) THEN ! OoB check
- TNapjp1 = xNapjp1 * moments_i(ip,ij+1,ikr,ikz,updatetlevel)
- ELSE
- TNapjp1 = 0._dp
- ENDIF
+ IF ( j_int+1 .LE. jmaxi ) TNapjp1 = xNapjp1 * moments_i(ip,ij+1,ikr,ikz,updatetlevel)
! term propto N_i^{p,j-1}
- IF (j_int-1 .GE. 0) THEN ! OoB check
- TNapjm1 = xNapjm1 * moments_i(ip,ij-1,ikr,ikz,updatetlevel)
- ELSE
- TNapjm1 = 0._dp
- ENDIF
+ IF ( j_int-1 .GE. 0 ) TNapjm1 = xNapjm1 * moments_i(ip,ij-1,ikr,ikz,updatetlevel)
!! Collision
IF (CO .EQ. -3) THEN ! Gyrokin. Dougherty Collision terms
CALL DoughertyGK_i(ip,ij,ikr,ikz,TColl)
ELSEIF (CO .EQ. -2) THEN ! Dougherty Collision terms
- IF ( (pmaxi .GE. 2) ) THEN ! OoB check
- TColl20 = xCa20 * moments_i(3,1,ikr,ikz,updatetlevel)
- ELSE
- TColl20 = 0._dp
- ENDIF
- IF ( (jmaxi .GE. 1) ) THEN ! OoB check
- TColl01 = xCa01 * moments_i(1,2,ikr,ikz,updatetlevel)
- ELSE
- TColl01 = 0._dp
- ENDIF
- IF ( (pmaxi .GE. 1) ) THEN ! OoB check
- TColl10 = xCa10 * moments_i(2,1,ikr,ikz,updatetlevel)
- ELSE
- TColl10 = 0._dp
- ENDIF
+ TColl20 = 0._dp; TColl01 = 0._dp; TColl10 = 0._dp
+ IF ( (pmaxe .GE. 2) ) TColl20 = xCa20 * moments_i(3,1,ikr,ikz,updatetlevel)
+ IF ( (jmaxe .GE. 1) ) TColl01 = xCa01 * moments_i(1,2,ikr,ikz,updatetlevel)
+ IF ( (pmaxe .GE. 1) ) TColl10 = xCa10 * moments_i(2,1,ikr,ikz,updatetlevel)
! Total collisional term
TColl = xCapj* moments_i(ip,ij,ikr,ikz,updatetlevel)&
+ TColl20 + TColl01 + TColl10
@@ -417,16 +360,9 @@ SUBROUTINE moments_eq_rhs
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
kernelj = kernel_i(ij, ikr, ikz)
- IF (j_int+1 .LE. jmaxi) THEN
- kerneljp1 = kernel_i(ij+1, ikr, ikz)
- ELSE
- kerneljp1 = 0.0_dp
- ENDIF
- IF ( j_int-1 .GE. 0 ) THEN
- kerneljm1 = kernel_i(ij-1, ikr, ikz)
- ELSE
- kerneljm1 = 0._dp
- ENDIF
+ kerneljp1 = 0._dp; kerneljm1 = 0._dp
+ IF ( j_int+1 .LE. jmaxe ) kerneljp1 = kernel_i(ij+1, ikr, ikz)
+ IF ( j_int-1 .GE. 0 ) kerneljm1 = kernel_i(ij-1, ikr, ikz)
Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
ELSE
Tphi = 0._dp
@@ -455,4 +391,4 @@ SUBROUTINE moments_eq_rhs
CALL cpu_time(t1_rhs)
tc_rhs = tc_rhs + (t1_rhs-t0_rhs)
-END SUBROUTINE moments_eq_rhs
+END SUBROUTINE moments_eq_rhs_i
diff --git a/src/stepon.F90 b/src/stepon.F90
index 07906739aca661ae73584ac61a81dedaa4981d2e..56522ca160112d9e8fd9a27cff3fc7a9b8bd77ec 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -17,7 +17,8 @@ SUBROUTINE stepon
DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
! Compute right hand side of moments hierarchy equation
- CALL moments_eq_rhs
+ CALL moments_eq_rhs_e
+ CALL moments_eq_rhs_i
! Advance from updatetlevel to updatetlevel+1 (according to num. scheme)
CALL advance_time_level