diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 0445ae1bee607b7306c2752f9cf0499495be1f3d..ed45fad1e2b78eec06e42b199050e7c32547daf4 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -280,7 +280,7 @@ CONTAINS
SUBROUTINE set_kygrid
USE prec_const
- USE model, ONLY: LINEARITY
+ USE model, ONLY: LINEARITY, N_HD
IMPLICIT NONE
INTEGER :: i_, in, istart, iend
Nky = Ny/2+1 ! Defined only on positive kx since fields are real
@@ -293,7 +293,7 @@ CONTAINS
deltaky = 2._dp*PI/Ly
ky_max = Nky*deltaky
ky_min = deltaky
- diff_ky_coeff= (1._dp/ky_max)**4
+ diff_ky_coeff= (1._dp/ky_max)**N_HD
ENDIF
! Build the full grids on process 0 to diagnose it without comm
ALLOCATE(kyarray_full(1:Nky))
@@ -345,7 +345,7 @@ CONTAINS
SUBROUTINE set_kxgrid(shear)
USE prec_const
- USE model, ONLY: LINEARITY
+ USE model, ONLY: LINEARITY, N_HD
IMPLICIT NONE
REAL(dp), INTENT(IN) :: shear
INTEGER :: i_, counter
@@ -375,7 +375,7 @@ CONTAINS
deltakx = 2._dp*PI/Lx
kx_max = (Nkx/2)*deltakx
kx_min = deltakx
- diff_kx_coeff= (1._dp/kx_max)**4
+ diff_kx_coeff= (1._dp/kx_max)**N_HD
! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
local_kxmax = 0._dp
DO ikx = ikxs,ikxe
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 2e1f4d9d055da973e9dd5476fb7bb10b3f5b4e8a..043c5b02da8ce832b9e45381678865884a8d453f 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -16,10 +16,11 @@ MODULE model
REAL(dp), PUBLIC, PROTECTED :: mu_z = 0._dp ! spatial z-Hyperdiffusivity coefficient (for num. stability)
REAL(dp), PUBLIC, PROTECTED :: mu_p = 0._dp ! kinetic para hyperdiffusivity coefficient (for num. stability)
REAL(dp), PUBLIC, PROTECTED :: mu_j = 0._dp ! kinetic perp hyperdiffusivity coefficient (for num. stability)
- REAL(dp), PUBLIC, PROTECTED :: nu = 1._dp ! Collision frequency
+ INTEGER, PUBLIC, PROTECTED :: N_HD = 4 ! order of numerical spatial diffusion
+ REAL(dp), PUBLIC, PROTECTED :: nu = 0._dp ! Collision frequency
REAL(dp), PUBLIC, PROTECTED :: tau_e = 1._dp ! Temperature
REAL(dp), PUBLIC, PROTECTED :: tau_i = 1._dp !
- REAL(dp), PUBLIC, PROTECTED :: sigma_e = 1._dp ! Mass
+ REAL(dp), PUBLIC, PROTECTED :: sigma_e = 0.023338_dp! sqrt of electron ion mass ratio
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 !
@@ -57,7 +58,7 @@ CONTAINS
IMPLICIT NONE
NAMELIST /MODEL_PAR/ CLOS, NL_CLOS, KERN, LINEARITY, KIN_E, &
- mu_x, mu_y, mu_z, mu_p, mu_j, nu,&
+ mu_x, mu_y, N_HD, mu_z, mu_p, mu_j, nu,&
tau_e, tau_i, sigma_e, sigma_i, q_e, q_i,&
K_n, K_T, K_E, GradB, CurvB, lambdaD
@@ -109,9 +110,11 @@ CONTAINS
CALL attach(fidres, TRIM(str), "KIN_E", KIN_E)
CALL attach(fidres, TRIM(str), "nu", nu)
CALL attach(fidres, TRIM(str), "mu", 0)
- CALL attach(fidres, TRIM(str), "mu_x", mu_x)
- CALL attach(fidres, TRIM(str), "mu_y", mu_y)
- CALL attach(fidres, TRIM(str), "mu_z", mu_z)
+ CALL attach(fidres, TRIM(str), "mu_x", mu_x)
+ CALL attach(fidres, TRIM(str), "mu_y", mu_y)
+ CALL attach(fidres, TRIM(str), "mu_z", mu_z)
+ CALL attach(fidres, TRIM(str), "mu_p", mu_p)
+ CALL attach(fidres, TRIM(str), "mu_j", mu_j)
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)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 3a493643875aa862251991ed6eed83655f2d8625..c30af83562b13fee14fa59c0fd7725bc7ebaa321 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -37,7 +37,7 @@ SUBROUTINE moments_eq_rhs_e
COMPLEX(dp) :: Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments
COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi
COMPLEX(dp) :: Unepm1j, Unepm1jp1, Unepm1jm1 ! Terms from mirror force with adiab moments
- COMPLEX(dp) :: i_ky,phikykxz
+ COMPLEX(dp) :: i_kx,i_ky,phikykxz
! Measuring execution time
CALL cpu_time(t0_rhs)
@@ -46,6 +46,7 @@ SUBROUTINE moments_eq_rhs_e
zloope : DO iz = izs,ize
kxloope : DO ikx = ikxs,ikxe
kx = kxarray(ikx) ! radial wavevector
+ i_kx = imagu * kx ! toroidal derivative
kyloope : DO iky = ikys,ikye
ky = kyarray(iky) ! toroidal wavevector
@@ -110,11 +111,10 @@ SUBROUTINE moments_eq_rhs_e
! Drives (density + temperature gradients)
- i_ky * Tphi &
! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose)
- - mu_x*diff_kx_coeff*kx**4*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
- - mu_y*diff_ky_coeff*ky**4*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
- ! - (mu_x*kx**2 + mu_y*ky**2)*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
+ - mu_x*diff_kx_coeff*kx**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
+ - mu_y*diff_ky_coeff*ky**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" see Pueschel 2010 (eq 25)
- + mu_z * diff_dz_coeff * ddz4_Nepj(ip,ij,iky,ikx,iz) &
+ - mu_z * diff_dz_coeff * ddz4_Nepj(ip,ij,iky,ikx,iz) &
! Collision term
+ TColl_e(ip,ij,iky,ikx,iz) &
! Nonlinear term
@@ -166,7 +166,7 @@ SUBROUTINE moments_eq_rhs_i
COMPLEX(dp) :: 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) :: Tperp, Tpar, Tmir, Tphi
- COMPLEX(dp) :: i_ky, phikykxz
+ COMPLEX(dp) :: i_kx, i_ky, phikykxz
! Measuring execution time
CALL cpu_time(t0_rhs)
@@ -176,7 +176,7 @@ SUBROUTINE moments_eq_rhs_i
zloopi : DO iz = izs,ize
kxloopi : DO ikx = ikxs,ikxe
kx = kxarray(ikx) ! radial wavevector
-
+ i_kx = imagu * kx ! toroidal derivative
kyloopi : DO iky = ikys,ikye
ky = kyarray(iky) ! toroidal wavevector
i_ky = imagu * ky ! toroidal derivative
@@ -242,8 +242,8 @@ SUBROUTINE moments_eq_rhs_i
! Drives (density + temperature gradients)
- i_ky * Tphi &
! Numerical hyperdiffusion (totally artificial, for stability purpose)
- - (mu_x*kx**4 + mu_y*ky**4)*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
- ! - (mu_x*kx**2 + mu_y*ky**2)*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
+ + mu_x*diff_kx_coeff*i_kx**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
+ + mu_y*diff_ky_coeff*i_ky**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)"
+ mu_z * diff_dz_coeff * ddz4_Nipj(ip,ij,iky,ikx,iz) &
! Collision term