SUBROUTINE moments_eq_rhs
USE basic
USE time_integration
USE array
USE fields
USE fourier_grid
USE model
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, 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 ! 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 !!!!!!!!!
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)
! N_e^{p+2,j} multiplicator
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. 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)
IF (ij_dp .GT. 0) THEN
factj = factj * ij_dp; ! Recursive factorial
ENDIF
! N_e^{p,j+1} multiplicator
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. 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.)
! Collision operator (DK Lenard-Bernstein basis)
xColl = ip_dp + 2.*ij_dp
! 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
ELSE IF (ip .EQ. 3) THEN !(kronecker delta_p^2)
xphij = (eta_T/SQRT2 - SQRT2*eta_B)
xphijp1 = 0.; xphijm1 = 0.
ELSE
xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
ENDIF
! 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 * 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. 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. 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. 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. 1) THEN
TNapjm1 = moments_e(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
ELSE
TNapjm1 = 0.
ENDIF
! 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)
!! Electrical potential term
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)
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
END DO jloope
END DO ploope
!!!!!!!!! Ions moments RHS !!!!!!!!!
ploopi : DO ip = ips_i, ipe_i
ip_dp = REAL(ip-1,dp)
! x N_i^{p+2,j}
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. 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)
IF (ij_dp .GT. 0) THEN
factj = factj * ij_dp; ! Recursive factorial
ENDIF
! x N_i^{p,j+1}
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. 1) THEN
xNapjm1 = taui_qi_etaB * ij_dp
ELSE
xNapjm1 = 0.
ENDIF
! x N_i^{pj}
xNapj = -taui_qi_etaB * 2.*(ip_dp + ij_dp + 1.)
! 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
ELSE IF (ip .EQ. 3) THEN !(krokecker delta_p^2)
xphij = (eta_T/SQRT2 - SQRT2*eta_B)
xphijp1 = 0.; xphijm1 = 0.
ELSE
xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
ENDIF
! Loop on kspace
krloopi : DO ikr = ikrs,ikre
kzloopi : DO ikz = ikzs,ikze
kr = krarray(ikr) ! Poloidal wavevector
kz = kzarray(ikz) ! Toroidal wavevector
kperp2 = kr**2 + kz**2 ! perpendicular wavevector
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. 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. 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. 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. 1) THEN
TNapjm1 = moments_i(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
ELSE
TNapjm1 = 0.
ENDIF
! 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
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)
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
END DO jloopi
END DO ploopi
END SUBROUTINE moments_eq_rhs