diff --git a/src/memory.F90 b/src/memory.F90
index 2d1dee0a5895d10d1b39c2e3d7c745b25799fa81..e6dc6a2abb935e9721364d5c2e9247b384d6722e 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -41,8 +41,8 @@ SUBROUTINE memory
CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
! Non linear terms and dnjs table
- CALL allocate_array( nadiab_moments_e, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( nadiab_moments_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( nadiab_moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( nadiab_moments_i, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize)
! Collision term
CALL allocate_array( TColl_e, ips_e,ipe_e, ijs_e,ije_e , ikxs,ikxe, ikys,ikye, izs,ize)
@@ -90,7 +90,8 @@ SUBROUTINE memory
CALL allocate_array( xphijm1, ips_i,ipe_i, ijs_i,ije_i)
! Curvature and geometry
- CALL allocate_array( Ckxky, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Ckxky, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( kparray, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array(Jacobian,izs,ize)
CALL allocate_array(gxx, izs,ize)
CALL allocate_array(gxy, izs,ize)
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 0dd19cb0898e95a3e0f6cfbaa690f6627f40d710..1660a18484201e08ca799f119ba0456d3f7b17b8 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -25,7 +25,7 @@ SUBROUTINE moments_eq_rhs_e
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
+ INTEGER :: izm2, izm1, izp1, izp2 ! indices for centered FDF ddz
! Measuring execution time
CALL cpu_time(t0_rhs)
@@ -39,60 +39,55 @@ SUBROUTINE moments_eq_rhs_e
jloope : DO ij = ijs_e, ije_e ! loop over Laguerre degree indices
j_int = jarray_e(ij)
- GF_CLOSURE_e : IF( (CLOS .EQ. 1) .AND. (p_int+2*j_int .GT. Dmaxe) ) THEN
- !skip
- ELSE
! Loop on kspace
zloope : DO iz = izs,ize
- ! Periodic BC for parallel centered finite differences
- iznext = iz+1; izprev = iz-1;
- IF(iz .EQ. 1) izprev = Nz
- IF(iz .EQ. Nz) iznext = 1
-
+ ! Obtain the index with an array that accounts for boundary conditions
+ ! e.g. : 4 stencil with periodic BC, izarray(Nz+2) = 2, izarray(-1) = Nz-1
+ izp1 = izarray(iz+1); izp2 = izarray(iz+2);
+ izm1 = izarray(iz-1); izm2 = izarray(iz-2);
+ !
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= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ ! kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ kperp2= kparray(ikx,iky,iz)**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)
+ Tnepj = xnepj(ip,ij)* nadiab_moments_e(ip,ij,ikx,iky,iz)
! term propto n_e^{p+2,j}
- Tnepp2j = xnepp2j(ip) * moments_e(ip+(2/deltape),ij,ikx,iky,iz,updatetlevel)
+ Tnepp2j = xnepp2j(ip) * nadiab_moments_e(ip+pp2,ij,ikx,iky,iz)
! term propto n_e^{p-2,j}
- Tnepm2j = xnepm2j(ip) * (moments_e(ip-(2/deltape),ij,ikx,iky,iz,updatetlevel) &
- +kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p2)
+ Tnepm2j = xnepm2j(ip) * nadiab_moments_e(ip-pp2,ij,ikx,iky,iz)
! 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)
+ Tnepjp1 = xnepjp1(ij) * nadiab_moments_e(ip,ij+1,ikx,iky,iz)
! 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)
+ Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,ikx,iky,iz)
! Parallel dynamic
- ! term propto centered FDF dz(n_a)
+ ! ddz derivative for Landau damping term
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)
-
+ ( onetwelfth*nadiab_moments_e(ip+1,ij,ikx,iky,izm2)&
+ - twothird*nadiab_moments_e(ip+1,ij,ikx,iky,izm1)&
+ + twothird*nadiab_moments_e(ip+1,ij,ikx,iky,izp1)&
+ -onetwelfth*nadiab_moments_e(ip+1,ij,ikx,iky,izp2))&
+ +xnepm1j(ip) * &
+ ( onetwelfth*nadiab_moments_e(ip-1,ij,ikx,iky,izm2)&
+ - twothird*nadiab_moments_e(ip-1,ij,ikx,iky,izm1)&
+ + twothird*nadiab_moments_e(ip-1,ij,ikx,iky,izp1)&
+ -onetwelfth*nadiab_moments_e(ip-1,ij,ikx,iky,izp2))
! 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)
+ Tnepp1j = ynepp1j(ip,ij) * nadiab_moments_e(ip+1,ij ,ikx,iky,iz)
+ Tnepp1jm1 = ynepp1jm1(ip,ij) * nadiab_moments_e(ip+1,ij-1,ikx,iky,iz)
+ Tnepm1j = ynepm1j(ip,ij) * nadiab_moments_e(ip-1,ij ,ikx,iky,iz)
+ Tnepm1jm1 = ynepm1jm1(ip,ij) * nadiab_moments_e(ip-1,ij-1,ikx,iky,iz)
+ ! Trapping terms
+ UNepm1j = zNepm1j(ip,ij) * nadiab_moments_e(ip-1,ij ,ikx,iky,iz)
+ UNepm1jp1 = zNepm1jp1(ip,ij) * nadiab_moments_e(ip-1,ij+1,ikx,iky,iz)
+ UNepm1jm1 = zNepm1jm1(ip,ij) * nadiab_moments_e(ip-1,ij-1,ikx,iky,iz)
Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + UNepm1j + UNepm1jp1 + UNepm1jm1
@@ -117,9 +112,9 @@ SUBROUTINE moments_eq_rhs_e
!! Sum of all linear terms (the sign is inverted to match RHS)
moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
- - imagu*Ckxky(ikx,iky,iz) * (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)&
+ - imagu*Ckxky(ikx,iky,iz)*hatR(iz)* (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)&
! Parallel coupling (Landau Damping)
- - Tpare/2._dp/deltaz*gradz_coeff(iz) &
+ - Tpare*inv_deltaz*gradz_coeff(iz) &
! Mirror term (parallel magnetic gradient)
- gradzB(iz)* Tmir *gradz_coeff(iz) &
! Drives (density + temperature gradients)
@@ -132,7 +127,7 @@ SUBROUTINE moments_eq_rhs_e
+ TColl
!! Adding non linearity
- IF ( NON_LIN .AND. (NL_CLOS .GT. -3)) THEN
+ 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
@@ -140,8 +135,6 @@ SUBROUTINE moments_eq_rhs_e
END DO kyloope
END DO kxloope
END DO zloope
-
- ENDIF GF_CLOSURE_e
END DO jloope
END DO ploope
@@ -176,7 +169,7 @@ SUBROUTINE moments_eq_rhs_i
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
+ INTEGER :: izm2, izm1, izp1, izp2 ! indices for centered FDF ddz
! Measuring execution time
CALL cpu_time(t0_rhs)
@@ -191,60 +184,55 @@ SUBROUTINE moments_eq_rhs_i
jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1
j_int = jarray_i(ij)
- GF_CLOSURE_i : IF( (CLOS .EQ. 1) .AND. (p_int+2*j_int .GT. Dmaxe) ) THEN
- !skip
- ELSE
! 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
-
+ ! Obtain the index with an array that accounts for boundary conditions
+ ! e.g. : 4 stencil with periodic BC, izarray(Nz+2) = 2, izarray(-1) = Nz-1
+ izp1 = izarray(iz+1); izp2 = izarray(iz+2);
+ izm1 = izarray(iz-1); izm2 = izarray(iz-2);
+ !
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= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ ! kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ kperp2= kparray(ikx,iky,iz)**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)
+ Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip,ij,ikx,iky,iz)
! term propto n_i^{p+2,j}
- Tnipp2j = xnipp2j(ip) * moments_i(ip+(2/deltapi),ij,ikx,iky,iz,updatetlevel)
+ Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij,ikx,iky,iz)
! term propto n_i^{p-2,j}
- Tnipm2j = xnipm2j(ip) * (moments_i(ip-(2/deltapi),ij,ikx,iky,iz,updatetlevel) &
- +kernel_i(ij,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p2)
+ Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij,ikx,iky,iz)
! 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)
+ Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip,ij+1,ikx,iky,iz)
! 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)
+ Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip,ij-1,ikx,iky,iz)
! 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)
-
+ ( onetwelfth*nadiab_moments_i(ip+1,ij,ikx,iky,izm2)&
+ - twothird*nadiab_moments_i(ip+1,ij,ikx,iky,izm1)&
+ + twothird*nadiab_moments_i(ip+1,ij,ikx,iky,izp1)&
+ -onetwelfth*nadiab_moments_i(ip+1,ij,ikx,iky,izp2))&
+ +xnipm1j(ip) * &
+ ( onetwelfth*nadiab_moments_i(ip-1,ij,ikx,iky,izm2)&
+ - twothird*nadiab_moments_i(ip-1,ij,ikx,iky,izm1)&
+ + twothird*nadiab_moments_i(ip-1,ij,ikx,iky,izp1)&
+ -onetwelfth*nadiab_moments_i(ip-1,ij,ikx,iky,izp2))
! 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)
+ Tnipp1j = ynipp1j(ip,ij) * nadiab_moments_i(ip+1,ij ,ikx,iky,iz)
+ Tnipp1jm1 = ynipp1jm1(ip,ij) * nadiab_moments_i(ip+1,ij-1,ikx,iky,iz)
+ Tnipm1j = ynipm1j(ip,ij) * nadiab_moments_i(ip-1,ij ,ikx,iky,iz)
+ Tnipm1jm1 = ynipm1jm1(ip,ij) * nadiab_moments_i(ip-1,ij-1,ikx,iky,iz)
+ ! Trapping terms
+ Unipm1j = znipm1j(ip,ij) * nadiab_moments_i(ip-1,ij ,ikx,iky,iz)
+ Unipm1jp1 = znipm1jp1(ip,ij) * nadiab_moments_i(ip-1,ij+1,ikx,iky,iz)
+ Unipm1jm1 = znipm1jm1(ip,ij) * nadiab_moments_i(ip-1,ij-1,ikx,iky,iz)
Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + UNipm1j + UNipm1jp1 + UNipm1jm1
@@ -267,13 +255,13 @@ SUBROUTINE moments_eq_rhs_i
ENDIF
!! Sum of all linear terms (the sign is inverted to match RHS)
- moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = &
+ moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
- - imagu*Ckxky(ikx,iky,iz) * (Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1)&
+ - imagu*Ckxky(ikx,iky,iz)*hatR(iz)*(Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1)&
! Parallel coupling (Landau Damping)
- - Tpari/2._dp/deltaz*gradz_coeff(iz) &
+ - Tpari*inv_deltaz*gradz_coeff(iz) &
! Mirror term (parallel magnetic gradient)
- - gradzB(iz)* Tmir *gradz_coeff(iz) &
+ - gradzB(iz)*Tmir*gradz_coeff(iz) &
! Drives (density + temperature gradients)
- i_ky * Tphi &
! Electrostatic background gradients
@@ -284,7 +272,7 @@ SUBROUTINE moments_eq_rhs_i
+ TColl
!! Adding non linearity
- IF ( NON_LIN .AND. (NL_CLOS .GT. -3)) THEN
+ 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
@@ -292,8 +280,6 @@ SUBROUTINE moments_eq_rhs_i
END DO kyloopi
END DO kxloopi
END DO zloopi
-
- ENDIF GF_CLOSURE_i
END DO jloopi
END DO ploopi
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 73583aa4c7dbb9f32652cb22ddd3d60116b9c789..768cc94315754a43bec2948f330672168099dcfe 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -143,7 +143,7 @@ SUBROUTINE compute_nadiab_moments
DO ij=ijsg_e,ijeg_e
nadiab_moments_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize)&
= moments_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) &
- + qe_taue*kernel_e(ij,ikxs:ikxe,ikys:ikye,izs:ize)*phi(ikx,iky,iz)
+ + qe_taue*kernel_e(ij,ikxs:ikxe,ikys:ikye,izs:ize)*phi(ikxs:ikxe,ikys:ikye,izs:ize)
ENDDO
ELSE
nadiab_moments_e(ip,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize) &
@@ -156,7 +156,7 @@ SUBROUTINE compute_nadiab_moments
DO ij=ijsg_i,ijeg_i
nadiab_moments_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize)&
= moments_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) &
- + qi_taui*kernel_i(ij,ikxs:ikxe,ikys:ikye,izs:ize)*phi(ikx,iky,iz)
+ + qi_taui*kernel_i(ij,ikxs:ikxe,ikys:ikye,izs:ize)*phi(ikxs:ikxe,ikys:ikye,izs:ize)
ENDDO
ELSE
nadiab_moments_i(ip,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize) &