diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index c8055420c488c07179728f05806b1b710f6027c0..d3c48200e77ac6c39be794af4a903c5127e1c698 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -1,7 +1,7 @@
module closure
! Contains the routines to define closures
USE basic
-USE model, ONLY: CLOS, tau_e, tau_i, q_e, q_i, eta_B, nu
+USE model, ONLY: CLOS, tau_e, tau_i, q_e, q_i, nu
USE grid
USE array, ONLY: kernel_e, kernel_i
USE fields, ONLY: moments_e, moments_i
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index 967e1039db9e1f93187521bd8857382a73b38297..ca6120cd4deecf850baf675b1860824977e3f02d 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -43,7 +43,9 @@ zloop: DO iz = izs,ize
real_data_c = 0._dp ! initialize sum over real nonlinear term
! Set non linear sum truncation
- IF (NL_CLOS .EQ. -2) THEN
+ IF (NL_CLOS .EQ. -3) THEN
+ return
+ ELSEIF (NL_CLOS .EQ. -2) THEN
nmax = Jmaxe
ELSEIF (NL_CLOS .EQ. -1) THEN
nmax = Jmaxe-(ij-1)
@@ -139,7 +141,9 @@ zloop: DO iz = izs,ize
real_data_c = 0._dp ! initialize sum over real nonlinear term
! Set non linear sum truncation
- IF (NL_CLOS .EQ. -2) THEN
+ IF (NL_CLOS .EQ. -3) THEN
+ return
+ ELSEIF (NL_CLOS .EQ. -2) THEN
nmax = Jmaxi
ELSEIF (NL_CLOS .EQ. -1) THEN
nmax = Jmaxi-(ij-1)
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index 053a1097e8ddac72e09392e34665b474883b9200..c3a4e82d92f404b8b1980a4e18ab23b1a9f56c50 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -32,11 +32,11 @@ END SUBROUTINE update_ghosts
! [a b|C D|e f] : proc n has moments a to f where a,b,e,f are ghosts
!
!proc 0: [0 1 2 3 4|5 6]
-! V V ^ ^
+! V V ^ ^
!proc 1: [3 4|5 6 7 8|9 10]
! V V ^ ^
!proc 2: [7 8|9 10 11 12|13 14]
-! V V ^ ^
+! V V ^ ^
!proc 3: [11 12|13 14 15 16|17 18]
! ^ ^
!Closure by zero truncation : 0 0
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 312b52d0dee2367599db28a50fc8b16a5c7ab135..1471d9757f22503616158b0c7d4fd34f02afe9b5 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -342,7 +342,6 @@ CONTAINS
DO ikp = ikps,ikpe
kparray(ikp) = REAL(ikp-1,dp) * deltakx
ENDDO
- write(*,*) rank_kx, ': ikps = ', ikps, 'ikpe = ',ikpe
two_third_kpmax = SQRT(two_third_kxmax**2+two_third_kymax**2)
kp_max = 3._dp/2._dp * two_third_kpmax
END SUBROUTINE
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index b89be50bc62e4b990bfae1ca82b4fbbbdf44e244..570bbd04ebda33997e2d8b6dbf6d16f7f73fed66 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -20,9 +20,11 @@ MODULE model
REAL(dp), PUBLIC, PROTECTED :: sigma_i = 1._dp !
REAL(dp), PUBLIC, PROTECTED :: q_e = -1._dp ! Charge
REAL(dp), PUBLIC, PROTECTED :: q_i = 1._dp !
- REAL(dp), PUBLIC, PROTECTED :: eta_n = 1._dp ! Density gradient
- REAL(dp), PUBLIC, PROTECTED :: eta_T = 1._dp ! Temperature gradient
- REAL(dp), PUBLIC, PROTECTED :: eta_B = 1._dp ! Magnetic gradient
+ REAL(dp), PUBLIC, PROTECTED :: K_n = 1._dp ! Density drive
+ REAL(dp), PUBLIC, PROTECTED :: K_T = 0._dp ! Temperature drive
+ REAL(dp), PUBLIC, PROTECTED :: K_E = 0._dp ! Backg. electric field drive
+ REAL(dp), PUBLIC, PROTECTED :: GradB = 1._dp ! Magnetic gradient
+ REAL(dp), PUBLIC, PROTECTED :: CurvB = 1._dp ! Magnetic curvature
REAL(dp), PUBLIC, PROTECTED :: lambdaD = 1._dp ! Debye length
REAL(dp), PUBLIC, PROTECTED :: taue_qe ! factor of the magnetic moment coupling
@@ -50,7 +52,7 @@ CONTAINS
IMPLICIT NONE
NAMELIST /MODEL_PAR/ CO, CLOS, NL_CLOS, KERN, NON_LIN, mu, mu_p, mu_j, nu, tau_e, tau_i, sigma_e, sigma_i, &
- q_e, q_i, eta_n, eta_T, eta_B, lambdaD
+ q_e, q_i, K_n, K_T, K_E, GradB, CurvB, lambdaD
READ(lu_in,model_par)
@@ -92,22 +94,22 @@ CONTAINS
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
- CALL attach(fidres, TRIM(str), "CO", CO)
- CALL attach(fidres, TRIM(str), "CLOS", CLOS)
- CALL attach(fidres, TRIM(str), "KERN", KERN)
- CALL attach(fidres, TRIM(str), "NON_LIN", NON_LIN)
- CALL attach(fidres, TRIM(str), "nu", nu)
- CALL attach(fidres, TRIM(str), "mu", mu)
- CALL attach(fidres, TRIM(str), "tau_e", tau_e)
- CALL attach(fidres, TRIM(str), "tau_i", tau_i)
- CALL attach(fidres, TRIM(str), "sigma_e", sigma_e)
- CALL attach(fidres, TRIM(str), "sigma_i", sigma_i)
- CALL attach(fidres, TRIM(str), "q_e", q_e)
- CALL attach(fidres, TRIM(str), "q_i", q_i)
- CALL attach(fidres, TRIM(str), "eta_n", eta_n)
- CALL attach(fidres, TRIM(str), "eta_T", eta_T)
- CALL attach(fidres, TRIM(str), "eta_B", eta_B)
- CALL attach(fidres, TRIM(str), "lambdaD", lambdaD)
+ CALL attach(fidres, TRIM(str), "CO", CO)
+ CALL attach(fidres, TRIM(str), "CLOS", CLOS)
+ CALL attach(fidres, TRIM(str), "KERN", KERN)
+ CALL attach(fidres, TRIM(str), "NON_LIN", NON_LIN)
+ CALL attach(fidres, TRIM(str), "nu", nu)
+ CALL attach(fidres, TRIM(str), "mu", mu)
+ CALL attach(fidres, TRIM(str), "tau_e", tau_e)
+ CALL attach(fidres, TRIM(str), "tau_i", tau_i)
+ CALL attach(fidres, TRIM(str), "sigma_e", sigma_e)
+ CALL attach(fidres, TRIM(str), "sigma_i", sigma_i)
+ CALL attach(fidres, TRIM(str), "q_e", q_e)
+ CALL attach(fidres, TRIM(str), "q_i", q_i)
+ CALL attach(fidres, TRIM(str), "K_n", K_n)
+ CALL attach(fidres, TRIM(str), "K_T", K_T)
+ CALL attach(fidres, TRIM(str), "K_E", K_E)
+ CALL attach(fidres, TRIM(str), "lambdaD", lambdaD)
END SUBROUTINE model_outputinputs
END MODULE model
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 606994486f825b4f26e2bb7b1bd88b3e2504d7b3..0dd19cb0898e95a3e0f6cfbaa690f6627f40d710 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -44,7 +44,7 @@ SUBROUTINE moments_eq_rhs_e
ELSE
! Loop on kspace
zloope : DO iz = izs,ize
- ! Periodic BC for first order derivative
+ ! 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
@@ -116,15 +116,23 @@ 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)&
+ ! Parallel coupling (Landau Damping)
- Tpare/2._dp/deltaz*gradz_coeff(iz) &
+ ! Mirror term (parallel magnetic gradient)
- gradzB(iz)* Tmir *gradz_coeff(iz) &
- - i_ky * Tphi &
+ ! Drives (density + temperature gradients)
+ - i_ky * Tphi &
+ ! Electrostatic background gradients
+ - i_ky * K_E * moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
+ ! Numerical hyperdiffusion (totally artificial, for stability purpose)
- mu*kperp2**2 * moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
+ ! Collision term
+ TColl
!! Adding non linearity
- IF ( NON_LIN ) THEN
+ IF ( NON_LIN .AND. (NL_CLOS .GT. -3)) 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
@@ -258,17 +266,25 @@ SUBROUTINE moments_eq_rhs_i
TColl = TColl_i(ip,ij,ikx,iky,iz)
ENDIF
- !! Sum of linear terms
- moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = &
+ !! 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) * (Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1)&
+ ! Parallel coupling (Landau Damping)
- Tpari/2._dp/deltaz*gradz_coeff(iz) &
+ ! Mirror term (parallel magnetic gradient)
- gradzB(iz)* Tmir *gradz_coeff(iz) &
- - i_ky * Tphi &
+ ! Drives (density + temperature gradients)
+ - i_ky * Tphi &
+ ! Electrostatic background gradients
+ - i_ky * K_E * moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
+ ! Numerical hyperdiffusion (totally artificial, for stability purpose)
- mu*kperp2**2 * moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
+ ! Collision term
+ TColl
!! Adding non linearity
- IF ( NON_LIN ) THEN
+ IF ( NON_LIN .AND. (NL_CLOS .GT. -3)) 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
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index fbbaf93a7226b1a3dc0ea181d00acb9be91c4a28..fc02b56551efb6d81609ec2d6e8d6ee06e4325d5 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -158,7 +158,8 @@ 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 model, ONLY: taue_qe, taui_qi, sqrtTaue_qe, sqrtTaui_qi, &
+ K_T, K_n, CurvB, GradB
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
@@ -173,7 +174,8 @@ SUBROUTINE compute_lin_coeff
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)
+ xnepj(ip,ij) = taue_qe*(CurvB*(2._dp*p_dp + 1._dp) &
+ +GradB*(2._dp*j_dp + 1._dp))
ynepp1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*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)
@@ -205,7 +207,8 @@ SUBROUTINE compute_lin_coeff
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)
+ xnepj(ip,ij) = taui_qi*(CurvB*(2._dp*p_dp + 1._dp) &
+ +GradB*(2._dp*j_dp + 1._dp))
ynipp1j (ip,ij) = -SQRT(tau_i)/sigma_i * (j_dp+1)*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)
@@ -238,11 +241,11 @@ SUBROUTINE compute_lin_coeff
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
- xphijp1(ip,ij) =-eta_T*(j_dp+1._dp)
- xphijm1(ip,ij) =-eta_T* j_dp
+ xphij(ip,ij) =+K_n + 2.*j_dp*K_T
+ xphijp1(ip,ij) =-K_T*(j_dp+1._dp)
+ xphijm1(ip,ij) =-K_T* j_dp
ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
- xphij(ip,ij) =+eta_T/SQRT2
+ xphij(ip,ij) =+K_T/SQRT2
xphijp1(ip,ij) = 0._dp; xphijm1(ip,ij) = 0._dp;
ELSE
xphij(ip,ij) = 0._dp; xphijp1(ip,ij) = 0._dp
diff --git a/src/srcinfo.h b/src/srcinfo.h
new file mode 100644
index 0000000000000000000000000000000000000000..04893679d4eb2fb6fcfd964e27c172e34f789636
--- /dev/null
+++ b/src/srcinfo.h
@@ -0,0 +1,10 @@
+character(len=40) VERSION
+character(len=40) BRANCH
+character(len=20) AUTHOR
+character(len=40) EXECDATE
+character(len=40) HOST
+parameter (VERSION='d118c2b-dirty')
+parameter (BRANCH='master')
+parameter (AUTHOR='ahoffman')
+parameter (EXECDATE='Thu Oct 7 17:18:16 CEST 2021')
+parameter (HOST ='spcpc606')
diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h
new file mode 100644
index 0000000000000000000000000000000000000000..04893679d4eb2fb6fcfd964e27c172e34f789636
--- /dev/null
+++ b/src/srcinfo/srcinfo.h
@@ -0,0 +1,10 @@
+character(len=40) VERSION
+character(len=40) BRANCH
+character(len=20) AUTHOR
+character(len=40) EXECDATE
+character(len=40) HOST
+parameter (VERSION='d118c2b-dirty')
+parameter (BRANCH='master')
+parameter (AUTHOR='ahoffman')
+parameter (EXECDATE='Thu Oct 7 17:18:16 CEST 2021')
+parameter (HOST ='spcpc606')