From 5ed2781da4fe08dc72d3537cefa22902d8f75454 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Mon, 15 Mar 2021 16:19:05 +0100
Subject: [PATCH] Corrections of Dougherty operator
---
src/collision_mod.F90 | 267 ++++++++++++++++++++++++------------------
1 file changed, 155 insertions(+), 112 deletions(-)
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index b236f536..406ba4c8 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -22,82 +22,103 @@ CONTAINS
USE grid, ONLY: jmaxe, parray_e, jarray_e, krarray, kzarray
USE prec_const
USE time_integration, ONLY: updatetlevel
- USE model, ONLY: sigmae2_taue_o2, tau_e, nu_e
+ USE model, ONLY: sigmae2_taue_o2, tau_e, nu_e, q_e
IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip,ij,ikr,ikz
+ INTEGER, INTENT(IN) :: ip,ij,ikr,ikz
COMPLEX(dp), INTENT(INOUT) :: TColl
COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_
COMPLEX(dp) :: Dpj, Ppj, T_, ibe_
COMPLEX(dp) :: nadiab_moment_0j, nadiab_moment_0jp1, nadiab_moment_0jm1
INTEGER :: in_
- REAL :: n_dp, j_dp, p_dp, be_2
+ REAL(dp) :: n_dp, j_dp, p_dp, be_2, q_e_tau_e
!** Auxiliary variables **
- p_dp = REAL(parray_e(ij),dp)
- j_dp = REAL(jarray_e(ij),dp)
- be_2 = (krarray(ikr)**2 + kzarray(ikz)**2) * sigmae2_taue_o2
- ibe_ = imagu*2._dp*SQRT(be_2)
-
- !** build fluid moments used for collison operator **
- n_ = 0._dp
- upar_ = 0._dp; uperp_ = 0._dp
- Tpar_ = 0._dp; Tperp_ = 0._dp
- DO in_ = 1,jmaxe+1
- n_dp = REAL(in_-1,dp)
-
- ! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = moments_e(1,in_,ikr,ikz,updatetlevel) + kernel_e(in_,ikr,ikz)*phi(ikr,ikz)/tau_e
- nadiab_moment_0jp1 = 0._dp; nadiab_moment_0jm1 = 0._dp
- IF (n_dp+1 .LE. jmaxe)& ! Truncation
- nadiab_moment_0jp1 = moments_e(1,in_+1,ikr,ikz,updatetlevel) + kernel_e(in_+1,ikr,ikz)*phi(ikr,ikz)/tau_e
- IF (n_dp-1 .GE. 0)& ! Truncation
- nadiab_moment_0jm1 = moments_e(1,in_-1,ikr,ikz,updatetlevel) + kernel_e(in_-1,ikr,ikz)*phi(ikr,ikz)/tau_e
-
- ! Density
- n_ = n_ + Kernel_e(in_,ikr,ikz) * nadiab_moment_0j
- ! Parallel velocity
- upar_ = upar_ + Kernel_e(in_,ikr,ikz) * moments_e(2,in_,ikr,ikz,updatetlevel)
- ! Perpendicular velocity
- uperp_ = uperp_ + ibe_*0.5_dp*Kernel_e(in_,ikr,ikz) * (nadiab_moment_0j -nadiab_moment_0jp1)
- ! Parallel temperature
- Tpar_ = Tpar_ + Kernel_e(in_,ikr,ikz) * (SQRT2*moments_e(3,in_,ikr,ikz,updatetlevel) + nadiab_moment_0j)
- ! Perpendicular temperature
- Tperp_ = Tperp_ + Kernel_e(in_,ikr,ikz) * ((2._dp*n_dp+1._dp)*nadiab_moment_0j - n_dp*nadiab_moment_0jm1 - (n_dp+1)*nadiab_moment_0jp1)
-
- ENDDO
- T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
+ p_dp = REAL(parray_e(ip),dp)
+ j_dp = REAL(jarray_e(ij),dp)
+ be_2 = (krarray(ikr)**2 + kzarray(ikz)**2) * sigmae2_taue_o2
+ ibe_ = imagu*2._dp*SQRT(be_2)
+ q_e_tau_e = q_e/tau_e
!** Assembling collison operator **
- ! Velocity-space diffusion
- TColl = -(2._dp*p_dp + j_dp + be_2)*moments_e(ip,ij,ikr,ikz,updatetlevel)
- IF( p_dp .EQ. 0 ) & ! Get adiabatic moment
- TColl = TColl - (2._dp*p_dp + j_dp + be_2) * Kernel_e(ij,ikr,ikz)*phi(ikr,ikz)/tau_e
-
- ! Add energy restoring term
- IF( p_dp .eq. 0 ) THEN ! kronecker p0
- TColl = TColl + T_* 4._dp*j_dp * Kernel_e(ij,ikr,ikz)
- IF( j_dp+1 .LE. jmaxe ) & ! Truncation
- TColl = TColl - T_*2._dp*(j_dp + 1._dp)*Kernel_e(ij+1,ikr,ikz)
- IF( j_dp-1 .GE. 0 ) & ! Truncation
- TColl = TColl - T_*2._dp*j_dp*Kernel_e(ij-1,ikr,ikz)
- ENDIF
- IF( p_dp .eq. 2 ) & ! kronecker p2
- TColl = TColl + T_*SQRT2*Kernel_e(ij,ikr,ikz)
-
- !Add momentum restoring term
- IF( p_dp .eq. 1 ) & ! kronecker p1
+ ! Velocity-space diffusion (similar to Lenhard Bernstein)
+ TColl = -(p_dp + 2._dp*j_dp + be_2)*moments_e(ip,ij,ikr,ikz,updatetlevel)
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
+ ! Get adiabatic moment
+ TColl = TColl - (p_dp + 2._dp*j_dp + be_2) * q_e_tau_e * Kernel_e(ij,ikr,ikz)*phi(ikr,ikz)
+ !** build required fluid moments **
+ n_ = 0._dp
+ upar_ = 0._dp; uperp_ = 0._dp
+ Tpar_ = 0._dp; Tperp_ = 0._dp
+ DO in_ = 1,jmaxe+1
+ n_dp = REAL(in_-1,dp)
+ ! Nonadiabatic moments (only different from moments when p=0)
+ nadiab_moment_0j = moments_e(1,in_ ,ikr,ikz,updatetlevel) + q_e_tau_e * kernel_e(in_ ,ikr,ikz)*phi(ikr,ikz)
+ nadiab_moment_0jp1 = moments_e(1,in_+1,ikr,ikz,updatetlevel) + q_e_tau_e * kernel_e(in_+1,ikr,ikz)*phi(ikr,ikz)
+ nadiab_moment_0jm1 = moments_e(1,in_-1,ikr,ikz,updatetlevel) + q_e_tau_e * kernel_e(in_-1,ikr,ikz)*phi(ikr,ikz)
+ ! Density
+ n_ = n_ + Kernel_e(in_,ikr,ikz) * nadiab_moment_0j
+ ! Perpendicular velocity
+ uperp_ = uperp_ + ibe_*0.5_dp*Kernel_e(in_,ikr,ikz) * (nadiab_moment_0j - nadiab_moment_0jp1)
+ ! Parallel temperature
+ Tpar_ = Tpar_ + Kernel_e(in_,ikr,ikz) * (SQRT2*moments_e(3,in_,ikr,ikz,updatetlevel) + nadiab_moment_0j)
+ ! Perpendicular temperature
+ Tperp_ = Tperp_ + Kernel_e(in_,ikr,ikz) * ((2._dp*n_dp+1._dp)* nadiab_moment_0j &
+ - n_dp * nadiab_moment_0jm1 &
+ - (n_dp+1)* nadiab_moment_0jp1)
+ ENDDO
+ T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
+ ! Add energy restoring term
+ TColl = TColl + T_* 4._dp * j_dp * Kernel_e(ij ,ikr,ikz)
+ TColl = TColl - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij+1,ikr,ikz)
+ TColl = TColl - T_* 2._dp * j_dp * Kernel_e(ij-1,ikr,ikz)
+ TColl = TColl + uperp_*ibe_*((j_dp + 1._dp)* Kernel_e(ij ,ikr,ikz) &
+ -j_dp * Kernel_e(ij-1,ikr,ikz))
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ELSEIF( p_dp .eq. 1 ) THEN ! kronecker p1
+ !** build required fluid moments **
+ upar_ = 0._dp
+ DO in_ = 1,jmaxe+1
+ ! Parallel velocity
+ upar_ = upar_ + Kernel_e(in_,ikr,ikz) * moments_e(2,in_,ikr,ikz,updatetlevel)
+ ENDDO
TColl = TColl + upar_*Kernel_e(ij,ikr,ikz)
- IF( p_dp .eq. 0 ) & ! kronecker p0
- TColl = TColl + uperp_*ibe_*( (j_dp + 1._dp)*Kernel_e(ij,ikr,ikz) - j_dp*Kernel_e(ij-1,ikr,ikz))
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ELSEIF( p_dp .eq. 2 ) THEN ! kronecker p2
+ !** build required fluid moments **
+ n_ = 0._dp
+ Tpar_ = 0._dp; Tperp_ = 0._dp
+ DO in_ = 1,jmaxe+1
+ n_dp = REAL(in_-1,dp)
+ ! Nonadiabatic moments (only different from moments when p=0)
+ nadiab_moment_0j = moments_e(1,in_ ,ikr,ikz,updatetlevel) + q_e_tau_e * kernel_e(in_ ,ikr,ikz)*phi(ikr,ikz)/tau_e
+ nadiab_moment_0jp1 = moments_e(1,in_+1,ikr,ikz,updatetlevel) + q_e_tau_e * kernel_e(in_+1,ikr,ikz)*phi(ikr,ikz)/tau_e
+ nadiab_moment_0jm1 = moments_e(1,in_-1,ikr,ikz,updatetlevel) + q_e_tau_e * kernel_e(in_-1,ikr,ikz)*phi(ikr,ikz)/tau_e
+ ! Density
+ n_ = n_ + Kernel_e(in_,ikr,ikz) * nadiab_moment_0j
+ ! Parallel temperature
+ Tpar_ = Tpar_ + Kernel_e(in_,ikr,ikz) * (SQRT2*moments_e(3,in_,ikr,ikz,updatetlevel) + nadiab_moment_0j)
+ ! Perpendicular temperature
+ Tperp_ = Tperp_ + Kernel_e(in_,ikr,ikz) * ((2._dp*n_dp+1._dp)*nadiab_moment_0j - n_dp*nadiab_moment_0jm1 - (n_dp+1)*nadiab_moment_0jp1)
+ ENDDO
+ T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
+ TColl = TColl + T_*SQRT2*Kernel_e(ij,ikr,ikz)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ENDIF
! Multiply by electron-ion collision coefficient
TColl = nu_e * TColl
END SUBROUTINE DoughertyGK_e
!******************************************************************************!
- !! Doughtery gyrokinetic collision operator for ions
+ !! Doughtery gyrokinetic collision operator for electrons
SUBROUTINE DoughertyGK_i(ip,ij,ikr,ikz,TColl)
USE fields, ONLY: moments_i, phi
USE array, ONLY: Kernel_i
@@ -105,74 +126,96 @@ CONTAINS
USE grid, ONLY: jmaxi, parray_i, jarray_i, krarray, kzarray
USE prec_const
USE time_integration, ONLY: updatetlevel
- USE model, ONLY: sigmai2_taui_o2, tau_i, nu_i
+ USE model, ONLY: sigmai2_taui_o2, tau_i, nu_i, q_i
IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip,ij,ikr,ikz
+ INTEGER, INTENT(IN) :: ip,ij,ikr,ikz
COMPLEX(dp), INTENT(INOUT) :: TColl
COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_
COMPLEX(dp) :: Dpj, Ppj, T_, ibi_
COMPLEX(dp) :: nadiab_moment_0j, nadiab_moment_0jp1, nadiab_moment_0jm1
INTEGER :: in_
- REAL :: n_dp, j_dp, p_dp, bi_2
+ REAL(dp) :: n_dp, j_dp, p_dp, bi_2, q_i_tau_i
!** Auxiliary variables **
- p_dp = REAL(parray_i(ij),dp)
- j_dp = REAL(jarray_i(ij),dp)
- bi_2 = (krarray(ikr)**2 + kzarray(ikz)**2) * sigmai2_taui_o2
- ibi_ = imagu*2._dp*SQRT(bi_2)
-
- !** build fluid moments used for collison operator **
- n_ = 0._dp
- upar_ = 0._dp; uperp_ = 0._dp
- Tpar_ = 0._dp; Tperp_ = 0._dp
- DO in_ = 1,jmaxi+1
- n_dp = REAL(in_-1,dp)
-
- ! Adiabatic moments (only different from moments when j=0)
- nadiab_moment_0j = moments_i(1,in_,ikr,ikz,updatetlevel) + Kernel_i(in_ ,ikr,ikz)*phi(ikr,ikz)/tau_i
- nadiab_moment_0jp1 = 0._dp; nadiab_moment_0jm1 = 0._dp
- IF (n_dp+1 .LE. jmaxi)&
- nadiab_moment_0jp1 = moments_i(1,in_+1,ikr,ikz,updatetlevel) + Kernel_i(in_+1,ikr,ikz)*phi(ikr,ikz)/tau_i
- IF (n_dp-1 .GE. 0)&
- nadiab_moment_0jm1 = moments_i(1,in_-1,ikr,ikz,updatetlevel) + Kernel_i(in_-1,ikr,ikz)*phi(ikr,ikz)/tau_i
-
- ! Density
- n_ = n_ + Kernel_i(in_,ikr,ikz) * nadiab_moment_0j
- ! Parallel velocity
- upar_ = upar_ + Kernel_i(in_,ikr,ikz) * moments_i(2,in_,ikr,ikz,updatetlevel)
- ! Perpendicular velocity
- uperp_ = uperp_ + ibi_*0.5_dp*Kernel_i(in_,ikr,ikz) * (nadiab_moment_0j -nadiab_moment_0jp1)
- ! Parallel temperature
- Tpar_ = Tpar_ + Kernel_i(in_,ikr,ikz) * (SQRT2*moments_i(3,in_,ikr,ikz,updatetlevel) + nadiab_moment_0j)
- ! Perpendicular temperature
- Tperp_ = Tperp_ + Kernel_i(in_,ikr,ikz) * ((2._dp*n_dp+1._dp)*nadiab_moment_0j - n_dp*nadiab_moment_0jm1 - (n_dp+1)*nadiab_moment_0jp1)
-
- ENDDO
- T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
+ p_dp = REAL(parray_i(ip),dp)
+ j_dp = REAL(jarray_i(ij),dp)
+ bi_2 = (krarray(ikr)**2 + kzarray(ikz)**2) * sigmai2_taui_o2
+ ibi_ = imagu*2._dp*SQRT(bi_2)
+ q_i_tau_i = q_i/tau_i
!** Assembling collison operator **
- ! Velocity-space diffusion
- TColl = -(2._dp*p_dp + j_dp + bi_2)*moments_i(ip,ij,ikr,ikz,updatetlevel)
- IF( p_dp .EQ. 0 ) & ! Get adiabatic moment
- TColl = TColl - (2._dp*p_dp + j_dp + bi_2) * Kernel_i(ij,ikr,ikz)*phi(ikr,ikz)/tau_i
-
- ! Add energy restoring term
- IF( p_dp .EQ. 0 ) THEN
- TColl = TColl + T_* 4._dp*j_dp * Kernel_i(ij,ikr,ikz)
- IF( j_dp+1 .LE. jmaxi ) &
- TColl = TColl - T_*2._dp*(j_dp + 1._dp)*Kernel_i(ij+1,ikr,ikz)
- IF( j_dp-1 .GE. 0 ) &
- TColl = TColl - T_*2._dp*j_dp*Kernel_i(ij-1,ikr,ikz)
- ENDIF
- IF( p_dp .eq. 2 ) &
+ ! Velocity-space diffusion (similar to Lenhard Bernstein)
+ TColl = -(p_dp + 2._dp*j_dp + bi_2)*moments_i(ip,ij,ikr,ikz,updatetlevel)
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
+ ! Get adiabatic moment
+ TColl = TColl - (p_dp + 2._dp*j_dp + bi_2) * q_i_tau_i * Kernel_i(ij,ikr,ikz)*phi(ikr,ikz)
+ !** build required fluid moments **
+ n_ = 0._dp
+ upar_ = 0._dp; uperp_ = 0._dp
+ Tpar_ = 0._dp; Tperp_ = 0._dp
+ DO in_ = 1,jmaxi+1
+ n_dp = REAL(in_-1,dp)
+ ! Nonadiabatic moments (only different from moments when p=0)
+ nadiab_moment_0j = moments_i(1,in_ ,ikr,ikz,updatetlevel) + q_i_tau_i * kernel_i(in_ ,ikr,ikz)*phi(ikr,ikz)
+ nadiab_moment_0jp1 = moments_i(1,in_+1,ikr,ikz,updatetlevel) + q_i_tau_i * kernel_i(in_+1,ikr,ikz)*phi(ikr,ikz)
+ nadiab_moment_0jm1 = moments_i(1,in_-1,ikr,ikz,updatetlevel) + q_i_tau_i * kernel_i(in_-1,ikr,ikz)*phi(ikr,ikz)
+ ! Density
+ n_ = n_ + Kernel_i(in_,ikr,ikz) * nadiab_moment_0j
+ ! Perpendicular velocity
+ uperp_ = uperp_ + ibi_*0.5_dp*Kernel_i(in_,ikr,ikz) * (nadiab_moment_0j -nadiab_moment_0jp1)
+ ! Parallel temperature
+ Tpar_ = Tpar_ + Kernel_i(in_,ikr,ikz) * (SQRT2*moments_i(3,in_,ikr,ikz,updatetlevel) + nadiab_moment_0j)
+ ! Perpendicular temperature
+ Tperp_ = Tperp_ + Kernel_i(in_,ikr,ikz) * ((2._dp*n_dp+1._dp)* nadiab_moment_0j &
+ - n_dp * nadiab_moment_0jm1 &
+ - (n_dp+1)* nadiab_moment_0jp1)
+ ENDDO
+ T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
+ ! Add energy restoring term
+ TColl = TColl + T_* 4._dp * j_dp * Kernel_i(ij ,ikr,ikz)
+ TColl = TColl - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij+1,ikr,ikz)
+ TColl = TColl - T_* 2._dp * j_dp * Kernel_i(ij-1,ikr,ikz)
+ TColl = TColl + uperp_*ibi_*((j_dp + 1._dp)* Kernel_i(ij ,ikr,ikz) &
+ -j_dp * Kernel_i(ij-1,ikr,ikz))
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ELSEIF( p_dp .eq. 1 ) THEN ! kronecker p1
+ !** build required fluid moments **
+ upar_ = 0._dp
+ DO in_ = 1,jmaxi+1
+ ! Parallel velocity
+ upar_ = upar_ + Kernel_i(in_,ikr,ikz) * moments_i(2,in_,ikr,ikz,updatetlevel)
+ ENDDO
+ TColl = TColl + upar_*Kernel_i(ij,ikr,ikz)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ELSEIF( p_dp .eq. 2 ) THEN ! kronecker p2
+ !** build required fluid moments **
+ n_ = 0._dp
+ Tpar_ = 0._dp; Tperp_ = 0._dp
+ DO in_ = 1,jmaxi+1
+ n_dp = REAL(in_-1,dp)
+ ! Nonadiabatic moments (only different from moments when p=0)
+ nadiab_moment_0j = moments_i(1,in_ ,ikr,ikz,updatetlevel) + q_i_tau_i * kernel_i(in_ ,ikr,ikz)*phi(ikr,ikz)
+ nadiab_moment_0jp1 = moments_i(1,in_+1,ikr,ikz,updatetlevel) + q_i_tau_i * kernel_i(in_+1,ikr,ikz)*phi(ikr,ikz)
+ nadiab_moment_0jm1 = moments_i(1,in_-1,ikr,ikz,updatetlevel) + q_i_tau_i * kernel_i(in_-1,ikr,ikz)*phi(ikr,ikz)
+ ! Density
+ n_ = n_ + Kernel_i(in_,ikr,ikz) * nadiab_moment_0j
+ ! Parallel temperature
+ Tpar_ = Tpar_ + Kernel_i(in_,ikr,ikz) * (SQRT2*moments_i(3,in_,ikr,ikz,updatetlevel) + nadiab_moment_0j)
+ ! Perpendicular temperature
+ Tperp_ = Tperp_ + Kernel_i(in_,ikr,ikz) * ((2._dp*n_dp+1._dp)*nadiab_moment_0j - n_dp*nadiab_moment_0jm1 - (n_dp+1)*nadiab_moment_0jp1)
+ ENDDO
+ T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
TColl = TColl + T_*SQRT2*Kernel_i(ij,ikr,ikz)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Add momentum restoring term
- IF( p_dp .eq. 1 ) &
- TColl = TColl + upar_*Kernel_i(ij,ikr,ikz)
- IF( p_dp .eq. 0 ) &
- TColl = TColl + uperp_*ibi_*( (j_dp + 1._dp)*Kernel_i(ij,ikr,ikz) - j_dp*Kernel_i(ij-1,ikr,ikz))
+ ENDIF
! Multiply by ion-ion collision coefficient
TColl = nu_i * TColl
--
GitLab