diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index 07a5a4dbe01b4b6bcc9a12de8deefba22c74acec..185fa13b3afb92de78d3d5274a26692e02e8049c 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -21,14 +21,11 @@ SUBROUTINE compute_Sapj
INTEGER :: in, is
REAL(dp):: kr, kz, kerneln
- REAL(dp):: sigmae2_taue_o2, sigmai2_taui_o2
! Execution time start
CALL cpu_time(t0_Sapj)
!!!!!!!!!!!!!!!!!!!! ELECTRON non linear term computation (Sepj)!!!!!!!!!!
- sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the kerneln argument
-
ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
real_data_c = 0._dp ! initialize sum over real nonlinear term
@@ -101,8 +98,6 @@ SUBROUTINE compute_Sapj
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!! ION non linear term computation (Sipj)!!!!!!!!!!
- sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp ! factor of the kerneln argument
-
ploopi: DO ip = ips_e,ipe_e ! Loop over Hermite moments
jloopi: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
real_data_c = 0._dp ! initialize sum over real nonlinear term
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 58604cf58b28a9922094500c838d6c54fa8732a0..1b0967a12d91e485e6820cf248b1d63a150f4887 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -21,6 +21,18 @@ MODULE model
REAL(dp), PUBLIC, PROTECTED :: eta_B = 1._dp ! Magnetic gradient
REAL(dp), PUBLIC, PROTECTED :: lambdaD = 1._dp ! Debye length
+ REAL(dp), PUBLIC, PROTECTED :: taue_qe_etaB ! factor of the magnetic moment coupling
+ REAL(dp), PUBLIC, PROTECTED :: taui_qi_etaB !
+ REAL(dp), PUBLIC, PROTECTED :: sqrtTaue_qe ! factor of parallel moment term
+ REAL(dp), PUBLIC, PROTECTED :: sqrtTaui_qi !
+ REAL(dp), PUBLIC, PROTECTED :: qe_sigmae_sqrtTaue ! factor of parallel phi term
+ REAL(dp), PUBLIC, PROTECTED :: qi_sigmai_sqrtTaui !
+ REAL(dp), PUBLIC, PROTECTED :: sigmae2_taue_o2 ! factor of the Kernel argument
+ REAL(dp), PUBLIC, PROTECTED :: sigmai2_taui_o2 !
+ REAL(dp), PUBLIC, PROTECTED :: nu_e, nu_i ! electron-ion, ion-ion collision frequency
+ REAL(dp), PUBLIC, PROTECTED :: nu_ee, nu_ie ! e-e, i-e coll. frequ.
+ REAL(dp), PUBLIC, PROTECTED :: qe2_taue, qi2_taui ! factor of the gammaD sum
+
PUBLIC :: model_readinputs, model_outputinputs
CONTAINS
@@ -41,6 +53,28 @@ CONTAINS
! Collision Frequency Normalization ... to match fluid limit
nu = nu*0.532_dp
+ !Precompute species dependant factors
+ IF( q_e .NE. 0._dp ) THEN
+ taue_qe_etaB = tau_e/q_e * eta_B ! factor of the magnetic moment coupling
+ sqrtTaue_qe = sqrt(tau_e)/q_e ! factor of parallel moment term
+ ELSE
+ taue_qe_etaB = 0._dp
+ sqrtTaue_qe = 0._dp
+ ENDIF
+
+ taui_qi_etaB = tau_i/q_i * eta_B ! factor of the magnetic moment coupling
+ sqrtTaui_qi = sqrt(tau_i)/q_i ! factor of parallel moment term
+ qe_sigmae_sqrtTaue = q_e/sigma_e/SQRT(tau_e) ! factor of parallel phi term
+ qi_sigmai_sqrtTaui = q_i/sigma_i/SQRT(tau_i)
+ qe2_taue = (q_e**2)/tau_e ! factor of the gammaD sum
+ qi2_taui = (q_i**2)/tau_i
+ sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument
+ sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp
+ nu_e = nu ! electron-ion collision frequency (where already multiplied by 0.532)
+ nu_i = nu * sigma_e * (tau_i)**(-3._dp/2._dp)/SQRT2 ! ion-ion collision frequ.
+ nu_ee = nu_e/SQRT2 ! e-e coll. frequ.
+ nu_ie = nu*sigma_e**2 ! i-e coll. frequ.
+
END SUBROUTINE model_readinputs
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 13b955ac321f6da385496cc0969b5b42c1b2cb7b..3b3c5fdee72e727b36bd684b1ad9a30ee66b764a 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -8,46 +8,23 @@ SUBROUTINE moments_eq_rhs
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) :: taue_qe_etaB, taui_qi_etaB
- REAL(dp) :: sqrtTaue_qe, sqrtTaui_qi, qe_sigmae_sqrtTaue, qi_sigmai_sqrtTaui
REAL(dp) :: kernelj, kerneljp1, kerneljm1 ! Kernel functions and variable
- REAL(dp) :: factj, sigmae2_taue_o2, sigmai2_taui_o2 ! Auxiliary variables
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) :: test_nan, i_kz
- REAL(dp) :: nu_e, nu_i, nu_ee, nu_ie ! Species collisional frequency
! Measuring execution time
CALL cpu_time(t0_rhs)
- !Precompute species dependant factors
- IF( q_e .NE. 0._dp ) THEN
- taue_qe_etaB = tau_e/q_e * eta_B ! factor of the magnetic moment coupling
- sqrtTaue_qe = sqrt(tau_e)/q_e ! factor of parallel moment term
- ELSE
- taue_qe_etaB = 0._dp
- sqrtTaue_qe = 0._dp
- ENDIF
-
- taui_qi_etaB = tau_i/q_i * eta_B ! factor of the magnetic moment coupling
- sqrtTaui_qi = sqrt(tau_i)/q_i ! factor of parallel moment term
- qe_sigmae_sqrtTaue = q_e/sigma_e/SQRT(tau_e) ! factor of parallel phi term
- qi_sigmai_sqrtTaui = q_i/sigma_i/SQRT(tau_i)
- sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument
- sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp
- nu_e = nu ! electron-ion collision frequency (where already multiplied by 0.532)
- nu_i = nu * sigma_e * (tau_i)**(-3._dp/2._dp)/SQRT2 ! ion-ion collision frequ.
- nu_ee = nu_e/SQRT2 ! e-e coll. frequ.
- nu_ie = nu*sigma_e**2 ! i-e coll. frequ.
-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!! Electrons moments RHS !!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -65,7 +42,7 @@ 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
+ ! 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
@@ -116,11 +93,11 @@ SUBROUTINE moments_eq_rhs
ENDIF
! Recursive factorial
- IF (j_dp .GT. 0) THEN
- factj = factj * j_dp
- ELSE
- factj = 1._dp
- ENDIF
+ ! IF (j_dp .GT. 0) THEN
+ ! factj = factj * j_dp
+ ! ELSE
+ ! factj = 1._dp
+ ! ENDIF
! Loop on kspace
krloope : DO ikr = ikrs,ikre
@@ -137,46 +114,49 @@ SUBROUTINE moments_eq_rhs
! term propto N_e^{p,j}
TNapj = xNapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
! term propto N_e^{p+1,j} and kparallel
- IF ( (ip+1 .LE. pmaxe+1) ) THEN ! OoB check
+ IF ( (p_int+1 .LE. pmaxe) ) THEN ! OoB check
TNapp1j = xNapp1j * moments_e(ip+1,ij,ikr,ikz,updatetlevel)
ELSE
TNapp1j = 0._dp
ENDIF
! term propto N_e^{p-1,j} and kparallel
- IF ( (ip-1 .GE. 1) ) THEN ! OoB check
+ IF ( (p_int-1 .GE. 0) ) THEN ! OoB check
TNapm1j = xNapm1j * moments_e(ip-1,ij,ikr,ikz,updatetlevel)
ELSE
TNapm1j = 0._dp
ENDIF
! term propto N_e^{p+2,j}
- IF (ip+(2-pskip) .LE. pmaxe+1) THEN ! OoB check
- TNapp2j = xNapp2j * moments_e(ip+(2-pskip),ij,ikr,ikz,updatetlevel)
+ IF (p_int+2 .LE. pmaxe) THEN ! OoB check
+ TNapp2j = xNapp2j * moments_e(ip+2,ij,ikr,ikz,updatetlevel)
ELSE
TNapp2j = 0._dp
ENDIF
! term propto N_e^{p-2,j}
- IF (ip-(2-pskip) .GE. 1) THEN ! OoB check
- TNapm2j = xNapm2j * moments_e(ip-(2-pskip),ij,ikr,ikz,updatetlevel)
+ IF (p_int-2 .GE. 0) THEN ! OoB check
+ TNapm2j = xNapm2j * moments_e(ip-2,ij,ikr,ikz,updatetlevel)
ELSE
TNapm2j = 0._dp
ENDIF
! xterm propto N_e^{p,j+1}
- IF (ij+1 .LE. jmaxe+1) THEN ! OoB check
+ IF (j_int+1 .LE. jmaxe) THEN ! OoB check
TNapjp1 = xNapjp1 * moments_e(ip,ij+1,ikr,ikz,updatetlevel)
ELSE
TNapjp1 = 0._dp
ENDIF
! term propto N_e^{p,j-1}
- IF (ij-1 .GE. 1) THEN ! OoB check
+ IF (j_int-1 .GE. 0) THEN ! OoB check
TNapjm1 = xNapjm1 * moments_e(ip,ij-1,ikr,ikz,updatetlevel)
ELSE
TNapjm1 = 0._dp
ENDIF
!! Collision
- IF (CO .EQ. -2) THEN ! Dougherty Collision terms
- IF ( (pmaxe .GE. 2-pskip) ) THEN ! OoB check
- TColl20 = xCa20 * moments_e(3-pskip,1,ikr,ikz,updatetlevel)
+ 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
@@ -190,7 +170,6 @@ SUBROUTINE moments_eq_rhs
ELSE
TColl10 = 0._dp
ENDIF
-
! Total collisional term
TColl = xCapj* moments_e(ip,ij,ikr,ikz,updatetlevel)&
+ TColl20 + TColl01 + TColl10
@@ -203,22 +182,17 @@ SUBROUTINE moments_eq_rhs
p2_int = parray_e(ip2)
jloopee: DO ij2 = 1,jmaxe+1
j2_int = jarray_e(ij2)
-
TColl = TColl + moments_e(ip2,ij2,ikr,ikz,updatetlevel) &
*( nu_e * CeipjT(bare(p_int,j_int), bare(p2_int,j2_int)) &
+nu_ee * Ceepj (bare(p_int,j_int), bare(p2_int,j2_int)))
-
ENDDO jloopee
ENDDO ploopee
-
ploopei: DO ip2 = 1,pmaxi+1 ! sum the electron-ion field terms
p2_int = parray_i(ip2)
jloopei: DO ij2 = 1,jmaxi+1
j2_int = jarray_i(ij2)
-
TColl = TColl + moments_i(ip2,ij2,ikr,ikz,updatetlevel) &
*(nu_e * CeipjF(bare(p_int,j_int), bari(p2_int,j2_int)))
-
END DO jloopei
ENDDO ploopei
@@ -285,7 +259,7 @@ SUBROUTINE moments_eq_rhs
! 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
+ ! 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
@@ -335,12 +309,12 @@ 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
+ ! ! 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
@@ -355,46 +329,49 @@ SUBROUTINE moments_eq_rhs
! term propto N_i^{p,j}
TNapj = xNapj * moments_i(ip,ij,ikr,ikz,updatetlevel)
! term propto N_i^{p+1,j}
- IF ( (ip+1 .LE. pmaxi+1) ) THEN ! OoB check
+ IF ( (p_int+1 .LE. pmaxi) ) THEN ! OoB check
TNapp1j = xNapp1j * moments_i(ip+1,ij,ikr,ikz,updatetlevel)
ELSE
TNapp1j = 0._dp
ENDIF
! term propto N_i^{p-1,j}
- IF ( (ip-1 .GE. 1) ) THEN ! OoB check
+ IF ( (p_int-1 .GE. 0) ) THEN ! OoB check
TNapm1j = xNapm1j * moments_i(ip-1,ij,ikr,ikz,updatetlevel)
ELSE
TNapm1j = 0._dp
ENDIF
! term propto N_i^{p+2,j}
- IF (ip+(2-pskip) .LE. pmaxi+1) THEN ! OoB check
- TNapp2j = xNapp2j * moments_i(ip+(2-pskip),ij,ikr,ikz,updatetlevel)
+ IF (p_int+2 .LE. pmaxi) THEN ! OoB check
+ TNapp2j = xNapp2j * moments_i(ip+2,ij,ikr,ikz,updatetlevel)
ELSE
TNapp2j = 0._dp
ENDIF
! term propto N_i^{p-2,j}
- IF (ip-(2-pskip) .GE. 1) THEN ! OoB check
- TNapm2j = xNapm2j * moments_i(ip-(2-pskip),ij,ikr,ikz,updatetlevel)
+ IF (p_int-2 .GE. 0) THEN ! OoB check
+ TNapm2j = xNapm2j * moments_i(ip-2,ij,ikr,ikz,updatetlevel)
ELSE
TNapm2j = 0._dp
ENDIF
! xterm propto N_i^{p,j+1}
- IF (ij+1 .LE. jmaxi+1) THEN ! OoB check
+ IF (j_int+1 .LE. jmaxi) THEN ! OoB check
TNapjp1 = xNapjp1 * moments_i(ip,ij+1,ikr,ikz,updatetlevel)
ELSE
TNapjp1 = 0._dp
ENDIF
! term propto N_i^{p,j-1}
- IF (ij-1 .GE. 1) THEN ! OoB check
+ IF (j_int-1 .GE. 0) THEN ! OoB check
TNapjm1 = xNapjm1 * moments_i(ip,ij-1,ikr,ikz,updatetlevel)
ELSE
TNapjm1 = 0._dp
ENDIF
!! Collision
- IF (CO .EQ. -2) THEN ! Dougherty Collision terms
- IF ( (pmaxi .GE. 2-pskip) ) THEN ! OoB check
- TColl20 = xCa20 * moments_i(3-pskip,1,ikr,ikz,updatetlevel)
+ 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
@@ -408,40 +385,33 @@ SUBROUTINE moments_eq_rhs
ELSE
TColl10 = 0._dp
ENDIF
-
! Total collisional term
TColl = xCapj* moments_i(ip,ij,ikr,ikz,updatetlevel)&
+ TColl20 + TColl01 + TColl10
ELSEIF (CO .EQ. -1) THEN !!! Full Coulomb for ions (COSOlver matrix) !!!
-
TColl = 0._dp ! Initialization
-
ploopii: DO ip2 = 1,pmaxi+1 ! sum the ion-self and ion-electron test terms
p2_int = parray_i(ip2)
jloopii: DO ij2 = 1,jmaxi+1
j2_int = jarray_i(ij2)
-
TColl = TColl + moments_i(ip2,ij2,ikr,ikz,updatetlevel) &
*( nu_ie * CiepjT(bari(p_int,j_int), bari(p2_int,j2_int)) &
+nu_i * Ciipj (bari(p_int,j_int), bari(p2_int,j2_int)))
-
ENDDO jloopii
ENDDO ploopii
-
ploopie: DO ip2 = 1,pmaxe+1 ! sum the ion-electron field terms
p2_int = parray_e(ip2)
jloopie: DO ij2 = 1,jmaxe+1
j2_int = jarray_e(ij2)
-
TColl = TColl + moments_e(ip2,ij2,ikr,ikz,updatetlevel) &
*(nu_ie * CiepjF(bari(p_int,j_int), bare(p2_int,j2_int)))
-
ENDDO jloopie
ENDDO ploopie
ELSEIF (CO .EQ. 0) THEN! Lenhard Bernstein
TColl = xCapj * moments_i(ip,ij,ikr,ikz,updatetlevel)
+
ENDIF
!! Electrical potential term
diff --git a/src/poisson.F90 b/src/poisson.F90
index f126d60ba58f8e662d7ca241a8471ff008cb7e4b..dc5c93a4558a319efdbef0bd51aae93e19b00d6d 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -6,7 +6,7 @@ SUBROUTINE poisson
USE array
USE fields
USE grid
- use model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, DK
+ use model, ONLY : qe2_taue, qi2_taui, q_e, q_i, lambdaD
USE prec_const
IMPLICIT NONE
@@ -14,7 +14,6 @@ SUBROUTINE poisson
INTEGER :: ini,ine
REAL(dp) :: Kne, Kni ! sub kernel factor for recursive build
REAL(dp) :: alphaD
- REAL(dp) :: qe2_taue, qi2_taui ! To avoid redondant computation
REAL(dp) :: sum_kernel2_e, sum_kernel2_i ! Store sum Kn^2
COMPLEX(dp) :: sum_kernel_mom_e, sum_kernel_mom_i ! Store sum Kn*Napn
REAL(dp) :: gammaD
@@ -23,10 +22,6 @@ SUBROUTINE poisson
! Execution time start
CALL cpu_time(t0_poisson)
- !Precompute species dependant factors
- qe2_taue = (q_e**2)/tau_e ! factor of the gammaD sum
- qi2_taui = (q_i**2)/tau_i
-
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze