diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90 index 58ccf7f..1914e38 100644 --- a/src/moments_eq_rhs.F90 +++ b/src/moments_eq_rhs.F90 @@ -1,279 +1,279 @@ !_____________________________________________________________________________! !_____________________________________________________________________________! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! Electrons moments RHS !!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !_____________________________________________________________________________! SUBROUTINE moments_eq_rhs_e USE basic USE time_integration USE array USE fields USE grid USE model USE prec_const USE collision use geometry IMPLICIT NONE INTEGER :: p_int, j_int ! loops indices and polynom. degrees REAL(dp) :: kx, ky, kperp2, dzlnB_o_J COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1, Tpare, Tphi ! Terms from b x gradB and drives COMPLEX(dp) :: Tmir, Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments COMPLEX(dp) :: UNepm1j, UNepm1jp1, UNepm1jm1 ! Terms from mirror force with adiab moments COMPLEX(dp) :: TColl ! terms of the rhs COMPLEX(dp) :: i_ky REAL(dp) :: delta_p0, delta_p1, delta_p2 INTEGER :: izprev,iznext ! indices of previous and next z slice ! Measuring execution time CALL cpu_time(t0_rhs) ploope : DO ip = ips_e, ipe_e ! loop over Hermite degree indices p_int= parray_e(ip) ! Hermite polynom. degree delta_p0 = 0._dp; delta_p1 = 0._dp; delta_p2 = 0._dp IF(p_int .EQ. 0) delta_p0 = 1._dp IF(p_int .EQ. 1) delta_p1 = 1._dp IF(p_int .EQ. 2) delta_p2 = 1._dp jloope : DO ij = ijs_e, ije_e ! loop over Laguerre degree indices ! Loop on kspace zloope : DO iz = izs,ize ! Periodic BC for first order derivative iznext = iz+1; izprev = iz-1; IF(iz .EQ. 1) izprev = Nz IF(iz .EQ. Nz) iznext = 1 kxloope : DO ikx = ikxs,ikxe kx = kxarray(ikx) ! radial wavevector kyloope : DO iky = ikys,ikye ky = kyarray(iky) ! toroidal wavevector i_ky = imagu * ky ! toroidal derivative IF (Nky .EQ. 1) i_ky = imagu * kxarray(ikx) ! If 1D simulation we put kx as ky kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2 !! Compute moments mixing terms ! Perpendicular dynamic ! term propto n_e^{p,j} Tnepj = xnepj(ip,ij) * (moments_e(ip,ij,ikx,iky,iz,updatetlevel) & +kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p0) ! term propto n_e^{p+2,j} Tnepp2j = xnepp2j(ip) * moments_e(ip+2,ij,ikx,iky,iz,updatetlevel) ! term propto n_e^{p-2,j} Tnepm2j = xnepm2j(ip) * (moments_e(ip-2,ij,ikx,iky,iz,updatetlevel) & +kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p2) ! term propto n_e^{p,j+1} Tnepjp1 = xnepjp1(ij) * (moments_e(ip,ij+1,ikx,iky,iz,updatetlevel) & +kernel_e(ij+1,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p0) ! term propto n_e^{p,j-1} Tnepjm1 = xnepjm1(ij) * (moments_e(ip,ij-1,ikx,iky,iz,updatetlevel) & +kernel_e(ij-1,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p0) ! Parallel dynamic ! term propto centered FDF dz(n_a) Tpare = xnepp1j(ip) * & (moments_e(ip+1,ij,ikx,iky,iznext,updatetlevel)& -moments_e(ip+1,ij,ikx,iky,izprev,updatetlevel)) & +xnepm1j(ip) * & (moments_e(ip-1,ij,ikx,iky,iznext,updatetlevel)+kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,iznext)*delta_p1& -moments_e(ip-1,ij,ikx,iky,izprev,updatetlevel)-kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,izprev)*delta_p1) ! Mirror terms Tnepp1j = ynepp1j(ip,ij) * moments_e(ip+1, ij,ikx,iky,iz,updatetlevel) Tnepp1jm1 = ynepp1jm1(ip,ij) * moments_e(ip+1,ij-1,ikx,iky,iz,updatetlevel) Tnepm1j = ynepm1j(ip,ij) * (moments_e(ip-1, ij,ikx,iky,iz,updatetlevel) & +kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p1) Tnepm1jm1 = ynepm1jm1(ip,ij) * (moments_e(ip-1,ij-1,ikx,iky,iz,updatetlevel) & +kernel_e(ij-1,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p1) UNepm1j = zNepm1j(ip,ij) * moments_e(ip+1, ij,ikx,iky,iz,updatetlevel) UNepm1jp1 = zNepm1jp1(ip,ij) * moments_e(ip-1,ij+1,ikx,iky,iz,updatetlevel) UNepm1jm1 = zNepm1jm1(ip,ij) * moments_e(ip-1,ij-1,ikx,iky,iz,updatetlevel) Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + UNepm1j + UNepm1jp1 + UNepm1jm1 !! Electrical potential term IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 Tphi = phi(ikx,iky,iz) * (xphij(ip,ij)*kernel_e(ij,ikx,iky,iz) & + xphijp1(ip,ij)*kernel_e(ij+1,ikx,iky,iz) & + xphijm1(ip,ij)*kernel_e(ij-1,ikx,iky,iz) ) ELSE Tphi = 0._dp ENDIF !! Collision IF (CO .EQ. 0) THEN ! Lenhard Bernstein TColl = -nu_ee*(ip+2*ij-3)*moments_e(ip,ij,ikx,iky,iz,updatetlevel) ELSEIF (CO .EQ. 1) THEN ! GK Dougherty CALL DoughertyGK_e(ip,ij,ikx,iky,iz,TColl) ELSE ! COSOLver matrix TColl = TColl_e(ip,ij,ikx,iky,iz) ENDIF !! Sum of all linear terms (the sign is inverted to match RHS) moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = & - imagu*Ckxky(ikx,iky,iz) * (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)& - Tpare/2._dp/deltaz*gradz_coeff(iz) & - gradzB(iz)* Tmir *gradz_coeff(iz) & - + i_ky * Tphi & + - i_ky * Tphi & - mu*kperp2**2 * moments_e(ip,ij,ikx,iky,iz,updatetlevel) & + TColl !! Adding non linearity IF ( NON_LIN ) THEN moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = & moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) + Sepj(ip,ij,ikx,iky,iz) ENDIF END DO kyloope END DO kxloope END DO zloope END DO jloope END DO ploope ! 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 collision IMPLICIT NONE INTEGER :: p_int, j_int ! loops indices and polynom. degrees REAL(dp) :: kx, ky, kperp2 COMPLEX(dp) :: Tnipj, Tnipp2j, Tnipm2j, Tnipjp1, Tnipjm1, Tpari, Tphi COMPLEX(dp) :: Tmir, Tnipp1j, Tnipm1j, Tnipp1jm1, Tnipm1jm1 ! Terms from mirror force with non adiab moments COMPLEX(dp) :: UNipm1j, UNipm1jp1, UNipm1jm1 ! Terms from mirror force with adiab moments COMPLEX(dp) :: TColl ! terms of the rhs COMPLEX(dp) :: i_ky REAL(dp) :: delta_p0, delta_p1, delta_p2 INTEGER :: izprev,iznext ! indices of previous and next z slice ! Measuring execution time CALL cpu_time(t0_rhs) ploopi : DO ip = ips_i, ipe_i ! Hermite loop p_int= parray_i(ip) ! Hermite degree delta_p0 = 0._dp; delta_p1 = 0._dp; delta_p2 = 0._dp IF(p_int .EQ. 0) delta_p0 = 1._dp IF(p_int .EQ. 1) delta_p1 = 1._dp IF(p_int .EQ. 2) delta_p2 = 1._dp jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1 ! Loop on kspace zloopi : DO iz = izs,ize ! Periodic BC for first order derivative iznext = iz+1; izprev = iz-1; IF(iz .EQ. 1) izprev = Nz IF(iz .EQ. Nz) iznext = 1 kxloopi : DO ikx = ikxs,ikxe kx = kxarray(ikx) ! radial wavevector kyloopi : DO iky = ikys,ikye ky = kyarray(iky) ! toroidal wavevector i_ky = imagu * ky ! toroidal derivative IF (Nky .EQ. 1) i_ky = imagu * kxarray(ikx) ! If 1D simulation we put kx as ky kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2 !! Compute moments mixing terms ! Perpendicular dynamic ! term propto n_i^{p,j} Tnipj = xnipj(ip,ij) * (moments_i(ip,ij,ikx,iky,iz,updatetlevel) & +kernel_i(ij,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p0) ! term propto n_i^{p+2,j} Tnipp2j = xnipp2j(ip) * moments_i(ip+2,ij,ikx,iky,iz,updatetlevel) ! term propto n_i^{p-2,j} Tnipm2j = xnipm2j(ip) * (moments_i(ip-2,ij,ikx,iky,iz,updatetlevel) & +kernel_i(ij,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p2) ! term propto n_e^{p,j+1} Tnipjp1 = xnipjp1(ij) * (moments_i(ip,ij+1,ikx,iky,iz,updatetlevel) & +kernel_i(ij+1,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p0) ! term propto n_e^{p,j-1} Tnipjm1 = xnipjm1(ij) * (moments_i(ip,ij-1,ikx,iky,iz,updatetlevel) & +kernel_i(ij-1,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p0) ! Parallel dynamic ! term propto N_i^{p,j+1}, centered FDF Tpari = xnipp1j(ip) * & (moments_i(ip+1,ij,ikx,iky,iznext,updatetlevel)& -moments_i(ip+1,ij,ikx,iky,izprev,updatetlevel)) & +xnipm1j(ip) * & (moments_i(ip-1,ij,ikx,iky,iznext,updatetlevel)+qi_taui*kernel_i(ij,ikx,iky,iz)*phi(ikx,iky,iznext)*delta_p1& -moments_i(ip-1,ij,ikx,iky,izprev,updatetlevel)-qi_taui*kernel_i(ij,ikx,iky,iz)*phi(ikx,iky,izprev)*delta_p1) ! Mirror terms Tnipp1j = ynipp1j(ip,ij) * moments_i(ip+1, ij,ikx,iky,iz,updatetlevel) Tnipp1jm1 = ynipp1jm1(ip,ij) * moments_i(ip+1,ij-1,ikx,iky,iz,updatetlevel) Tnipm1j = ynipm1j(ip,ij) * (moments_i(ip-1, ij,ikx,iky,iz,updatetlevel) & +kernel_i(ij,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p1) Tnipm1jm1 = ynipm1jm1(ip,ij) * (moments_i(ip-1,ij-1,ikx,iky,iz,updatetlevel) & +kernel_i(ij-1,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p1) Unipm1j = znipm1j(ip,ij) * moments_i(ip+1,ij,ikx,iky,iz,updatetlevel) Unipm1jp1 = znipm1jp1(ip,ij) * moments_i(ip-1,ij+1,ikx,iky,iz,updatetlevel) Unipm1jm1 = znipm1jm1(ip,ij) * moments_i(ip-1,ij-1,ikx,iky,iz,updatetlevel) Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + UNipm1j + UNipm1jp1 + UNipm1jm1 !! Electrical potential term IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 Tphi = phi(ikx,iky,iz) * (xphij(ip,ij)*kernel_i(ij,ikx,iky,iz) & + xphijp1(ip,ij)*kernel_i(ij+1,ikx,iky,iz) & + xphijm1(ip,ij)*kernel_i(ij-1,ikx,iky,iz) ) ELSE Tphi = 0._dp ENDIF !! Collision IF (CO .EQ. 0) THEN ! Lenhard Bernstein TColl = -nu_i*(ip+2._dp*ij-3)*moments_i(ip,ij,ikx,iky,iz,updatetlevel) ELSEIF (CO .EQ. 1) THEN ! GK Dougherty CALL DoughertyGK_i(ip,ij,ikx,iky,iz,TColl) ELSE! COSOLver matrix (Sugama, Coulomb) TColl = TColl_i(ip,ij,ikx,iky,iz) ENDIF !! Sum of linear terms moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = & - imagu*Ckxky(ikx,iky,iz) * (Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1)& - Tpari/2._dp/deltaz*gradz_coeff(iz) & - gradzB(iz)* Tmir *gradz_coeff(iz) & - + i_ky * Tphi & + - i_ky * Tphi & - mu*kperp2**2 * moments_i(ip,ij,ikx,iky,iz,updatetlevel) & + TColl !! Adding non linearity IF ( NON_LIN ) THEN moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = & moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) + Sipj(ip,ij,ikx,iky,iz) ENDIF END DO kyloopi END DO kxloopi END DO zloopi END DO jloopi END DO ploopi ! Execution time end CALL cpu_time(t1_rhs) tc_rhs = tc_rhs + (t1_rhs-t0_rhs) END SUBROUTINE moments_eq_rhs_i diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90 index b215be7..1289550 100644 --- a/src/numerics_mod.F90 +++ b/src/numerics_mod.F90 @@ -1,256 +1,256 @@ MODULE numerics USE basic USE prec_const USE grid USE utility USE coeff implicit none PUBLIC :: compute_derivatives, build_dnjs_table, evaluate_kernels, compute_lin_coeff CONTAINS ! Compute the 2D particle temperature for electron and ions (sum over Laguerre) SUBROUTINE compute_derivatives END SUBROUTINE compute_derivatives !******************************************************************************! !!!!!!! Build the Laguerre-Laguerre coupling coefficient table for nonlin !******************************************************************************! SUBROUTINE build_dnjs_table USE array, Only : dnjs USE coeff IMPLICIT NONE INTEGER :: in, ij, is, J INTEGER :: n_, j_, s_ J = max(jmaxe,jmaxi) DO in = 1,J+1 ! Nested dependent loops to make benefit from dnjs symmetry n_ = in - 1 DO ij = in,J+1 j_ = ij - 1 DO is = ij,J+1 s_ = is - 1 dnjs(in,ij,is) = TO_DP(ALL2L(n_,j_,s_,0)) ! By symmetry dnjs(in,is,ij) = dnjs(in,ij,is) dnjs(ij,in,is) = dnjs(in,ij,is) dnjs(ij,is,in) = dnjs(in,ij,is) dnjs(is,ij,in) = dnjs(in,ij,is) dnjs(is,in,ij) = dnjs(in,ij,is) ENDDO ENDDO ENDDO END SUBROUTINE build_dnjs_table !******************************************************************************! !******************************************************************************! !!!!!!! Evaluate the kernels once for all !******************************************************************************! SUBROUTINE evaluate_kernels USE basic USE array, Only : kernel_e, kernel_i, gxy, gyy 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 REAL(dp) :: factj, j_dp, j_int REAL(dp) :: be_2, bi_2, alphaD REAL(dp) :: kx, ky, kperp2 !!!!! Electron kernels !!!!! !Precompute species dependant factors factj = 1.0 ! Start of the recursive factorial DO ij = 1, jmaxe+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= 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 ENDDO ! Kernels closure DO ikx = ikxs,ikxe kx = kxarray(ikx) DO iky = ikys,ikye ky = kyarray(iky) DO iz = izs,ize kperp2= 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= 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= 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 ENDDO END SUBROUTINE evaluate_kernels !******************************************************************************! SUBROUTINE compute_lin_coeff USE array USE model, ONLY: taue_qe, taui_qi, sqrtTaue_qe, sqrtTaui_qi, eta_T, eta_n USE prec_const USE grid, ONLY: parray_e, parray_i, jarray_e, jarray_i, & ip,ij, ips_e,ipe_e, ips_i,ipe_i, ijs_e,ije_e, ijs_i,ije_i IMPLICIT NONE INTEGER :: p_int, j_int ! polynom. degrees REAL(dp) :: p_dp, j_dp REAL(dp) :: kx, ky, z !! Electrons linear coefficients for moment RHS !!!!!!!!!! DO ip = ips_e, ipe_e p_int= parray_e(ip) ! Hermite degree p_dp = REAL(p_int,dp) ! REAL of Hermite degree DO ij = ijs_e, ije_e j_int= jarray_e(ij) ! Laguerre degree j_dp = REAL(j_int,dp) ! REAL of Laguerre degree xnepj(ip,ij) = taue_qe * 2._dp * (p_dp + j_dp + 1._dp) ynepp1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp+1._dp) - ynepp1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*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(2._dp*p_dp) - zNepm1jp1(ip,ij) = +SQRT(tau_e)/sigma_e * (j_dp+1_dp)*SQRT(2._dp*p_dp) - zNepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(2._dp*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)) 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 ENDDO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Ions linear coefficients for moment RHS !!!!!!!!!! DO ip = ips_i, ipe_i p_int= parray_i(ip) ! Hermite degree p_dp = REAL(p_int,dp) ! REAL of Hermite degree DO ij = ijs_i, ije_i j_int= jarray_i(ij) ! Laguerre degree j_dp = REAL(j_int,dp) ! REAL of Laguerre degree xnipj(ip,ij) = taui_qi * 2._dp * (p_dp + j_dp + 1._dp) ynipp1j (ip,ij) = -SQRT(tau_i)/sigma_i * (j_dp+1)*SQRT(p_dp+1._dp) - ynipp1jm1(ip,ij) = +SQRT(tau_i)/sigma_i * j_dp*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(2._dp*p_dp) - zNipm1jp1(ip,ij) = +SQRT(tau_i)/sigma_i * (j_dp+1_dp)*SQRT(2._dp*p_dp) - zNipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i * j_dp*SQRT(2._dp*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) 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)) 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 ENDDO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! ES linear coefficients for moment RHS !!!!!!!!!! DO ip = ips_i, ipe_i p_int= parray_i(ip) ! Hermite degree DO ij = ijs_i, ije_i j_int= jarray_i(ij) ! REALof Laguerre degree j_dp = REAL(j_int,dp) ! REALof Laguerre degree !! Electrostatic potential pj terms IF (p_int .EQ. 0) THEN ! kronecker p0 - xphij(ip,ij) = eta_n + 2.*j_dp*eta_T + xphij(ip,ij) =+eta_n + 2.*j_dp*eta_T xphijp1(ip,ij) =-eta_T*(j_dp+1._dp) xphijm1(ip,ij) =-eta_T* j_dp ELSE IF (p_int .EQ. 2) THEN ! kronecker p2 - xphij(ip,ij) = eta_T/SQRT2 + xphij(ip,ij) =+eta_T/SQRT2 xphijp1(ip,ij) = 0._dp; xphijm1(ip,ij) = 0._dp; ELSE xphij(ip,ij) = 0._dp; xphijp1(ip,ij) = 0._dp xphijm1(ip,ij) = 0._dp; ENDIF ENDDO ENDDO END SUBROUTINE compute_lin_coeff END MODULE numerics