diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 29e39214be9b1c7f300ec3862fc2581b14127b3a..e4b95996ac3575f8670fc8a47f5c34a64640f3ee 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -6,153 +6,154 @@ SUBROUTINE moments_eq_rhs
USE fields
USE fourier_grid
USE model
-
- use prec_const
+ USE prec_const
IMPLICIT NONE
INTEGER :: ip,ij, ikr,ikz ! loops indices
REAL(dp) :: ip_dp, ij_dp
- REAL(dp) :: kr, kz, kperp2
- REAL(dp) :: taue_qe_etaB, taui_qi_etaB
- REAL(dp) :: kernelj, kerneljp1, kerneljm1, b_e2, b_i2 ! Kernel functions and variable
- REAL(dp) :: factj, xb_e2, xb_i2 ! Auxiliary variables
- REAL(dp) :: xNapj, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! factors depending on the pj loop
- REAL(dp) :: xphij, xphijp1, xphijm1, xColl
+ REAL(dp) :: kr, kz, kperp2
+ REAL(dp) :: taue_qe_etaB, taui_qi_etaB
+ REAL(dp) :: kernelj, kerneljp1, kerneljm1, b_e2, b_i2 ! Kernel functions and variable
+ REAL(dp) :: factj, sigmae2_taue_o2, sigmai2_taui_o2 ! Auxiliary variables
+ REAL(dp) :: xNapj, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! factors depending on the pj loop
+ REAL(dp) :: xphij, xphijp1, xphijm1, xColl
COMPLEX(dp) :: TNapj, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi
COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs
!Precompute species dependant factors
- taue_qe_etaB = tau_e/q_e * eta_B
- xb_e2 = sigma_e**2 * tau_e!/2.0 ! species dependant factor of the Kernel, squared
- taui_qi_etaB = tau_i/q_i * eta_B
- xb_i2 = sigma_i**2 * tau_i!/2.0 ! species dependant factor of the Kernel, squared
+ taue_qe_etaB = tau_e/q_e * eta_B ! factor of the magnetic moment coupling
+ taui_qi_etaB = tau_i/q_i * eta_B
+ sigmae2_taue_o2 = sigma_e**2 * tau_e/2.0 ! factor of the Kernel argument
+ sigmai2_taui_o2 = sigma_i**2 * tau_i/2.0
!!!!!!!!! Electrons moments RHS !!!!!!!!!
- Tphi = 0 ! electrostatic potential term
ploope : DO ip = ips_e, ipe_e ! This loop is from 1 to pmaxe+1
- ip_dp = REAL(ip-1.,dp) ! REAL index is one minus the loop index (0 to pmax)
+ ip_dp = REAL(ip-1,dp) ! REAL index is one minus the loop index (0 to pmax)
! N_e^{p+2,j} multiplicator
- IF (ip+2 .LE. ipe_e) THEN
+ IF (ip+2 .LE. pmaxe+1) THEN
xNapp2j = -taue_qe_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.))
ELSE
xNapp2j = 0.
ENDIF
! N_e^{p-2,j} multiplicator
- IF (ip-2 .GE. ips_e) THEN
+ IF (ip-2 .GE. 1) THEN
xNapm2j = -taue_qe_etaB * SQRT(ip_dp*(ip_dp - 1.))
ELSE
xNapm2j = 0.
ENDIF
+ factj = 1.0 ! Start of the recursive factorial
+
jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxe+1
- ij_dp = REAL(ij-1.,dp) ! REAL index is one minus the loop index (0 to jmax)
+ ij_dp = REAL(ij-1,dp) ! REAL index is one minus the loop index (0 to jmax)
+
+ IF (ij_dp .GT. 0) THEN
+ factj = factj * ij_dp; ! Recursive factorial
+ ENDIF
! N_e^{p,j+1} multiplicator
- IF (ij+1 .LE. ije_e) THEN
+ IF (ij+1 .LE. jmaxe+1) THEN
xNapjp1 = taue_qe_etaB * (ij_dp + 1.)
ELSE
xNapjp1 = 0.
ENDIF
! N_e^{p,j-1} multiplicator
- IF (ij-1 .GE. ijs_e) THEN
+ IF (ij-1 .GE. 1) THEN
xNapjm1 = taue_qe_etaB * ij_dp
ELSE
xNapjm1 = 0.
ENDIF
! N_e^{pj} multiplicator
- xNapj = -taue_qe_etaB * 2.*(ip_dp + ij_dp + 1.)
+ xNapj = -taue_qe_etaB * 2.*(ip_dp + ij_dp + 1.)
! Collision operator (DK Lenard-Bernstein basis)
xColl = ip_dp + 2.*ij_dp
- ! ... adding Dougherty terms
- IF ( (ip .EQ. 1) .AND. (ij .EQ. 2) ) THEN ! kronecker p0 * j1
- TColl01 = 2.0/3.0*(SQRT2*moments_e(3,1,ikr,ikz,updatetlevel)&
- - 2.0*moments_e(1,2,ikr,ikz,updatetlevel))
- TColl20 = 0.0; TColl10 = 0.0;
- ELSEIF ( (ip .EQ. 3) .AND. (ij .EQ. 1) ) THEN ! kronecker p2 * j0
- TColl20 = -SQRT2/3.0*(SQRT2*moments_e(3,1,ikr,ikz,updatetlevel)&
- - 2.0*moments_e(1,2,ikr,ikz,updatetlevel))
- TColl10 = 0.0; TColl01 = 0.0;
- ELSEIF ( (ip .EQ. 2) .AND. (ij .EQ. 1) ) THEN ! kronecker p1 * j0
- TColl10 = moments_e(2,1,ikr,ikz,updatetlevel)
- TColl20 = 0.0; TColl01 = 0.0;
- ELSE
- TColl10 = 0.0; TColl20 = 0.0; TColl01 = 0.0;
- ENDIF
! phi multiplicator for different kernel numbers
IF (ip .EQ. 1) THEN !(kronecker delta_p^0)
xphij = (eta_n + 2.*ij_dp*eta_T - 2.*eta_B*(ij_dp+1.) )
xphijp1 = -(eta_T - eta_B)*(ij_dp+1.)
xphijm1 = -(eta_T - eta_B)* ij_dp
- factj = REAL(Factorial(ij-1),dp)
ELSE IF (ip .EQ. 3) THEN !(kronecker delta_p^2)
xphij = (eta_T/SQRT2 - SQRT2*eta_B)
- factj = REAL(Factorial(ij-1),dp)
+ xphijp1 = 0.; xphijm1 = 0.
ELSE
xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
- factj = 1
ENDIF
- !write(*,*) '(ip,ij) = (', ip,',', ij,')'
-
! Loop on kspace
krloope : DO ikr = ikrs,ikre
kzloope : DO ikz = ikzs,ikze
kr = krarray(ikr) ! Poloidal wavevector
kz = kzarray(ikz) ! Toroidal wavevector
kperp2 = kr**2 + kz**2 ! perpendicular wavevector
- b_e2 = kperp2 * xb_e2 ! Bessel argument
+ b_e2 = kperp2 * sigmae2_taue_o2 ! Bessel argument
!! Compute moments and mixing terms
! term propto N_e^{p,j}
TNapj = moments_e(ip,ij,ikr,ikz,updatetlevel) * xNapj
! term propto N_e^{p+2,j}
- IF (ip+2 .LE. ipe_e) THEN
+ IF (ip+2 .LE. pmaxe+1) THEN
TNapp2j = moments_e(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j
ELSE
TNapp2j = 0.
ENDIF
! term propto N_e^{p-2,j}
- IF (ip-2 .GE. ips_e) THEN
+ IF (ip-2 .GE. 1) THEN
TNapm2j = moments_e(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j
ELSE
TNapm2j = 0.
ENDIF
! xterm propto N_e^{p,j+1}
- IF (ij+1 .LE. ije_e) THEN
+ IF (ij+1 .LE. jmaxe+1) THEN
TNapjp1 = moments_e(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1
ELSE
TNapjp1 = 0.
ENDIF
! term propto N_e^{p,j-1}
- IF (ij-1 .GE. ijs_e) THEN
+ IF (ij-1 .GE. 1) THEN
TNapjm1 = moments_e(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
ELSE
TNapjm1 = 0.
ENDIF
- ! Collision term completed (DK Dougherty)
+ ! Dougherty Collision terms
+ IF ( (ip .EQ. 1) .AND. (ij .EQ. 2) .AND. (pmaxe .GE. 2) ) THEN ! kronecker p0 * j1
+ TColl01 = 2.0/3.0*(SQRT2*moments_e(3,1,ikr,ikz,updatetlevel)&
+ - 2.0*moments_e(1,2,ikr,ikz,updatetlevel))
+ TColl20 = 0.0; TColl10 = 0.0;
+ ELSEIF ( (ip .EQ. 3) .AND. (ij .EQ. 1) .AND. (jmaxe .GE. 1)) THEN ! kronecker p2 * j0
+ TColl20 = -SQRT2/3.0*(SQRT2*moments_e(3,1,ikr,ikz,updatetlevel)&
+ - 2.0*moments_e(1,2,ikr,ikz,updatetlevel))
+ TColl10 = 0.0; TColl01 = 0.0;
+ ELSEIF ( (ip .EQ. 2) .AND. (ij .EQ. 1) ) THEN ! kronecker p1 * j0
+ TColl10 = moments_e(2,1,ikr,ikz,updatetlevel)
+ TColl20 = 0.0; TColl01 = 0.0;
+ ELSE
+ TColl10 = 0.0; TColl20 = 0.0; TColl01 = 0.0;
+ ENDIF
+ ! Total collisional term
TColl = -nu * (xColl * moments_e(ip,ij,ikr,ikz,updatetlevel) &
- + TColl01 + TColl10 + TColl20)
+ + TColl01 + TColl10 + TColl20)
!! Electrical potential term
- Tphi = 0
- IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! 0 otherwise (krokecker delta_p^0)
+ IF ( (ip .EQ. 1) .OR. (ip .EQ. 3) ) THEN ! kronecker delta_p^0, delta_p^2
kernelj = b_e2**(ij-1) * exp(-b_e2)/factj
kerneljp1 = kernelj * b_e2 /(ij_dp + 1.)
kerneljm1 = kernelj * ij_dp / b_e2
- Tphi = (xphij*Kernelj + xphijp1*Kerneljp1 + xphijm1*Kerneljm1) * phi(ikr,ikz)
+ Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
+ ELSE
+ Tphi = 0
ENDIF
! Sum of all precomputed terms
moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl
-
+
END DO kzloope
END DO krloope
@@ -160,37 +161,41 @@ SUBROUTINE moments_eq_rhs
END DO ploope
!!!!!!!!! Ions moments RHS !!!!!!!!!
- Tphi = 0 ! electrostatic potential term
-
ploopi : DO ip = ips_i, ipe_i
- ip_dp = REAL(ip-1.,dp)
+ ip_dp = REAL(ip-1,dp)
! x N_i^{p+2,j}
- IF (ip+2 .LE. ipe_i) THEN
+ IF (ip+2 .LE. pmaxi+1) THEN
xNapp2j = -taui_qi_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.))
ELSE
xNapp2j = 0.
ENDIF
! x N_i^{p-2,j}
- IF (ip-2 .GE. ips_i) THEN
+ IF (ip-2 .GE. 1) THEN
xNapm2j = -taui_qi_etaB * SQRT(ip_dp * (ip_dp - 1.))
ELSE
xNapm2j = 0.
ENDIF
+ factj = 1.0 ! Start of the recursive factorial
+
jloopi : DO ij = ijs_i, ije_i
- ij_dp = REAL(ij-1.,dp)
+ ij_dp = REAL(ij-1,dp)
+
+ IF (ij_dp .GT. 0) THEN
+ factj = factj * ij_dp; ! Recursive factorial
+ ENDIF
! x N_i^{p,j+1}
- IF (ij+1 .LE. ije_i) THEN
+ IF (ij+1 .LE. jmaxi+1) THEN
xNapjp1 = taui_qi_etaB * (ij_dp + 1.)
ELSE
xNapjp1 = 0.
ENDIF
! x N_i^{p,j-1}
- IF (ij-1 .GE. ijs_i) THEN
+ IF (ij-1 .GE. 1) THEN
xNapjm1 = taui_qi_etaB * ij_dp
ELSE
xNapjm1 = 0.
@@ -199,38 +204,19 @@ SUBROUTINE moments_eq_rhs
! x N_i^{pj}
xNapj = -taui_qi_etaB * 2.*(ip_dp + ij_dp + 1.)
- ! Collision term completed (DK Dougherty)
- TColl = -nu * (xColl * moments_i(ip,ij,ikr,ikz,updatetlevel) &
- + TColl01 + TColl10 + TColl20)
-
- ! ... adding Dougherty terms
- IF ( (ip .EQ. 1) .AND. (ij .EQ. 2) ) THEN ! kronecker p0 * j1
- TColl01 = 2.0/3.0*(SQRT2*moments_i(3,1,ikr,ikz,updatetlevel)&
- - 2.0*moments_i(1,2,ikr,ikz,updatetlevel))
- TColl20 = 0.0; TColl10 = 0.0;
- ELSEIF ( (ip .EQ. 3) .AND. (ij .EQ. 1) ) THEN ! kronecker p2 * j0
- TColl20 = -SQRT2/3.0*(SQRT2*moments_i(3,1,ikr,ikz,updatetlevel)&
- - 2.0*moments_i(1,2,ikr,ikz,updatetlevel))
- TColl10 = 0.0; TColl01 = 0.0;
- ELSEIF ( (ip .EQ. 2) .AND. (ij .EQ. 1) ) THEN ! kronecker p1 * j0
- TColl10 = moments_i(2,1,ikr,ikz,updatetlevel)
- TColl20 = 0.0; TColl01 = 0.0;
- ELSE
- TColl10 = 0.0; TColl20 = 0.0; TColl01 = 0.0;
- ENDIF
+ ! Collision operator (DK Lenard-Bernstein basis)
+ xColl = ip_dp + 2.*ij_dp
! x phi
IF (ip .EQ. 1) THEN !(krokecker delta_p^0)
xphij = (eta_n + 2.*ij_dp*eta_T - 2.*eta_B*(ij_dp+1.) )
xphijp1 = -(eta_T - eta_B)*(ij_dp+1.)
xphijm1 = -(eta_T - eta_B)* ij_dp
- factj = REAL(Factorial(ij-1),dp)
ELSE IF (ip .EQ. 3) THEN !(krokecker delta_p^2)
xphij = (eta_T/SQRT2 - SQRT2*eta_B)
- factj = REAL(Factorial(ij-1),dp)
+ xphijp1 = 0.; xphijm1 = 0.
ELSE
xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
- factj = 1.
ENDIF
! Loop on kspace
@@ -239,52 +225,68 @@ SUBROUTINE moments_eq_rhs
kr = krarray(ikr) ! Poloidal wavevector
kz = kzarray(ikz) ! Toroidal wavevector
kperp2 = kr**2 + kz**2 ! perpendicular wavevector
- b_i2 = kperp2 * xb_i2 ! Bessel argument
+ b_i2 = kperp2 * sigmai2_taui_o2 ! Bessel argument
!! Compute moments and mixing terms
! term propto N_i^{p,j}
TNapj = moments_i(ip,ij,ikr,ikz,updatetlevel) * xNapj
! term propto N_i^{p+2,j}
- IF (ip+2 .LE. ipe_i) THEN
+ IF (ip+2 .LE. pmaxi+1) THEN
TNapp2j = moments_i(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j
ELSE
TNapp2j = 0.
ENDIF
! term propto N_i^{p-2,j}
- IF (ip-2 .GE. ips_i) THEN
+ IF (ip-2 .GE. 1) THEN
TNapm2j = moments_i(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j
ELSE
TNapm2j = 0.
ENDIF
! xterm propto N_i^{p,j+1}
- IF (ij+1 .LE. ije_i) THEN
+ IF (ij+1 .LE. jmaxi+1) THEN
TNapjp1 = moments_i(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1
ELSE
TNapjp1 = 0.
ENDIF
! term propto N_i^{p,j-1}
- IF (ij-1 .GE. ijs_i) THEN
+ IF (ij-1 .GE. 1) THEN
TNapjm1 = moments_i(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
ELSE
TNapjm1 = 0.
ENDIF
- ! Collision term completed (Dougherty)
- TColl = -nu* (xColl * moments_i(ip,ij,ikr,ikz,updatetlevel))
+ ! Dougherty Collision terms
+ IF ( (ip .EQ. 1) .AND. (ij .EQ. 2) .AND. (pmaxi .GE. 2)) THEN ! kronecker p0 * j1
+ TColl01 = 2.0/3.0*(SQRT2*moments_i(3,1,ikr,ikz,updatetlevel)&
+ - 2.0*moments_i(1,2,ikr,ikz,updatetlevel))
+ TColl20 = 0.0; TColl10 = 0.0;
+ ELSEIF ( (ip .EQ. 3) .AND. (ij .EQ. 1) .AND. (jmaxi .GE. 1)) THEN ! kronecker p2 * j0
+ TColl20 = -SQRT2/3.0*(SQRT2*moments_i(3,1,ikr,ikz,updatetlevel)&
+ - 2.0*moments_i(1,2,ikr,ikz,updatetlevel))
+ TColl10 = 0.0; TColl01 = 0.0;
+ ELSEIF ( (ip .EQ. 2) .AND. (ij .EQ. 1) ) THEN ! kronecker p1 * j0
+ TColl10 = moments_i(2,1,ikr,ikz,updatetlevel)
+ TColl20 = 0.0; TColl01 = 0.0;
+ ELSE
+ TColl10 = 0.0; TColl20 = 0.0; TColl01 = 0.0;
+ ENDIF
+ ! Total collisional term
+ TColl = -nu * (xColl * moments_e(ip,ij,ikr,ikz,updatetlevel) &
+ + TColl01 + TColl10 + TColl20)
!! Electrical potential term
- Tphi = 0
- IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! 0 otherwise (krokecker delta_p^0, delta_p^2)
+ IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! kronecker delta_p^0, delta_p^2
kernelj = b_i2**(ij-1) * exp(-b_i2)/factj
kerneljp1 = kernelj * b_i2 /(ij_dp + 1.)
kerneljm1 = kernelj * ij_dp / b_i2
- Tphi = (xphij*Kernelj + xphijp1*Kerneljp1 + xphijm1*Kerneljm1) * phi(ikr,ikz)
+ Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
+ ELSE
+ Tphi = 0
ENDIF
-
! Sum of all precomputed terms
moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl
-
+
END DO kzloopi
END DO krloopi
diff --git a/src/poisson.F90 b/src/poisson.F90
index 1819af6fe1338202ced63129a2a8c8a2f12897b5..b3b5d68422bdf06be70d62eeacfc3442a2e7b767 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -13,82 +13,74 @@ SUBROUTINE poisson
INTEGER :: ikr,ikz, ini,ine, i0j
REAL(dp) :: ini_dp, ine_dp
REAL(dp) :: kr, kz, kperp2
- REAL(dp) :: Kn, Kn2 ! sub kernel factor for recursive build
+ REAL(dp) :: Kne, Kni ! sub kernel factor for recursive build
REAL(dp) :: b_e2, b_i2, alphaD
REAL(dp) :: sigmae2_taue_o2, sigmai2_taui_o2, qe2_taue, qi2_taui ! To avoid redondant computation
- COMPLEX(dp) :: gammaD, gammaphi
- COMPLEX(dp) :: sum_kernel2_e, sum_kernel2_i ! Store to compute sum Kn^2
- COMPLEX(dp) :: sum_kernel_mom_e, sum_kernel_mom_i ! Store to compute sum Kn*Napn
-
+ 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
+ COMPLEX(dp) :: gammaD_phi
+
!Precompute species dependant factors
- sigmae2_taue_o2 = sigma_e**2 * tau_e/2.
+ sigmae2_taue_o2 = sigma_e**2 * tau_e/2. ! factor of the Kernel argument
sigmai2_taui_o2 = sigma_i**2 * tau_i/2.
- qe2_taue = (q_e**2)/tau_e
+ 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
- kr = krarray(ikr)
- kz = kzarray(ikz)
+ kr = krarray(ikr)
+ kz = kzarray(ikz)
kperp2 = kr**2 + kz**2
!! Electrons sum(Kernel * Ne0n) (skm) and sum(Kernel**2) (sk2)
b_e2 = kperp2 * sigmae2_taue_o2 ! non dim kernel argument (kperp2 sigma_a sqrt(2 tau_a))
! Initialization for n = 0 (ine = 1)
- Kn = 1.
- Kn2 = 1.
- sum_kernel_mom_e = moments_e(1,1,ikr,ikz,updatetlevel)
- sum_kernel2_e = Kn2
+ Kne = exp(-b_e2)
+ sum_kernel_mom_e = Kne*moments_e(1, 1, ikr, ikz, updatetlevel)
+ sum_kernel2_e = Kne**2
! loop over n only if the max polynomial degree is 1 or more
if (jmaxe .GT. 0) then
- DO ine=1,jmaxe+1 ! ine = n+1
+ DO ine=2,jmaxe+1 ! ine = n+1
ine_dp = REAL(ine-1,dp)
- Kn = Kn * b_e2/(ine_dp+1) ! We update iteratively the kernel functions (to spare factorial computations)
- Kn2 = Kn2 *(b_e2/(ine_dp+1))**2 ! the exp term is still missing but does not depend on ine ...
+ Kne = Kne * b_e2/ine_dp ! We update iteratively the kernel functions (to spare factorial computations)
- sum_kernel_mom_e = sum_kernel_mom_e + Kn * moments_e(1,ine,ikr,ikz,updatetlevel)
- sum_kernel2_e = sum_kernel2_e + Kn2 ! ... sum recursively ...
+ sum_kernel_mom_e = sum_kernel_mom_e + Kne * moments_e(1, ine, ikr, ikz, updatetlevel)
+ sum_kernel2_e = sum_kernel2_e + Kne**2 ! ... sum recursively ...
END DO
endif
- sum_kernel_mom_e = sum_kernel_mom_e * exp(-b_e2) ! ... multiply the final sum with the missing exponential term
- sum_kernel2_e = sum_kernel2_e * exp(-2.*b_e2) ! and its squared using exp(x)^2 = exp(2x).
!! Ions sum(Kernel * Ni0n) (skm) and sum(Kernel**2) (sk2)
b_i2 = kperp2 * sigmai2_taui_o2
! Initialization for n = 0 (ini = 1)
- Kn = 1.
- Kn2 = 1.
- sum_kernel_mom_i = moments_i(1, 1, ikr, ikz, updatetlevel)
- sum_kernel2_i = Kn2
+ Kni = exp(-b_i2)
+ sum_kernel_mom_i = Kni*moments_i(1, 1, ikr, ikz, updatetlevel)
+ sum_kernel2_i = Kni**2
! loop over n only if the max polynomial degree is 1 or more
if (jmaxi .GT. 0) then
- DO ini=1,jmaxi+1
+ DO ini=2,jmaxi+1
ini_dp = REAL(ini-1,dp) ! Real index (0 to jmax)
- Kn = Kn * b_i2/(ini_dp+1) ! We update iteratively the kernel functions this time for ions ...
- Kn2 = Kn2 *(b_i2/(ini_dp+1.))**2
+ Kni = Kni * b_i2/ini_dp ! We update iteratively the kernel functions this time for ions ...
- sum_kernel_mom_i = sum_kernel_mom_i + Kn * moments_i(1,ini,ikr,ikz,updatetlevel)
- sum_kernel2_i = sum_kernel2_i + Kn2 ! ... sum recursively ...
+ sum_kernel_mom_i = sum_kernel_mom_i + Kni * moments_i(1, ini, ikr, ikz, updatetlevel)
+ sum_kernel2_i = sum_kernel2_i + Kni**2 ! ... sum recursively ...
END DO
endif
- sum_kernel_mom_i = sum_kernel_mom_i * exp(-b_i2) ! ... multiply the final sum with the missing exponential term.
- sum_kernel2_i = sum_kernel2_i * exp(-2.*b_i2)
!!! Assembling the poisson equation
- alphaD = kperp2 * lambdaD**2.
+ alphaD = kperp2 * lambdaD**2
gammaD = alphaD + qe2_taue * (1. - sum_kernel2_e) &
+ qi2_taui * (1. - sum_kernel2_i)
-
- gammaphi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i ! gamma_D * phi term
+ gammaD_phi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i ! gamma_D * phi term
- phi(ikr, ikz) = gammaphi/gammaD
-
+ phi(ikr, ikz) = gammaD_phi/gammaD
+
END DO
END DO