diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index c82024512f1ef1a9e5741e006f6a3725870b8b60..2d2c5186e0422856c90814b2cdd076979439a40a 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -1,24 +1,121 @@
module collision
! contains the Hermite-Laguerre collision operators. Solved using COSOlver.
+USE array
+USE basic
+USE fields
+USE futils
+USE grid
+USE model
+USE prec_const
+USE time_integration
+USE utility
IMPLICIT NONE
+PRIVATE
+! Set the collision model to use
+! (Lenard-Bernstein: 'LB', Dougherty: 'DG', Sugama: 'SG', Lorentz: 'LR', Landau: 'LD')
+CHARACTER(len=2),PUBLIC,PROTECTED :: collision_model
+LOGICAL, PUBLIC, PROTECTED :: gyrokin_CO =.true. ! activates GK effects on CO
+LOGICAL, PUBLIC, PROTECTED :: interspecies =.true. ! activates interpecies collision
+CHARACTER(len=128), PUBLIC :: mat_file ! COSOlver matrix file names
+LOGICAL, PUBLIC, PROTECTED :: cosolver_coll ! check if cosolver matrices are used
+
+PUBLIC :: collision_readinputs, coll_outputinputs
PUBLIC :: compute_TColl
+PUBLIC :: compute_lenard_bernstein, compute_dougherty
PUBLIC :: LenardBernstein_e, LenardBernstein_i!, LenardBernstein GK
PUBLIC :: DoughertyGK_e, DoughertyGK_i!, Dougherty GK
-PUBLIC :: load_COSOlver_mat
+PUBLIC :: load_COSOlver_mat, compute_cosolver_coll
PUBLIC :: apply_COSOlver_mat_e, apply_COSOlver_mat_i
CONTAINS
+ SUBROUTINE collision_readinputs
+ ! Read the input parameters
+ IMPLICIT NONE
+ NAMELIST /COLLISION_PAR/ collision_model, gyrokin_CO, interspecies, mat_file
+ READ(lu_in,collision_par)
+ SELECT CASE(collision_model)
+ CASE ('LB') ! Lenhard bernstein
+ cosolver_coll = .false.
+ interspecies = .false.
+ CASE ('DG') ! Dougherty
+ cosolver_coll = .false.
+ interspecies = .false.
+ CASE ('SG') ! Sugama
+ cosolver_coll = .true.
+ CASE ('LR') ! Lorentz (Pitch angle)
+ cosolver_coll = .true.
+ interspecies = .false.
+ CASE ('LD') ! Landau
+ cosolver_coll = .true.
+ END SELECT
+
+ END SUBROUTINE collision_readinputs
+
+ SUBROUTINE coll_outputinputs(fidres, str)
+ ! Write the input parameters to the results_xx.h5 file
+ IMPLICIT NONE
+ INTEGER, INTENT(in) :: fidres
+ CHARACTER(len=256), INTENT(in) :: str
+ CHARACTER(len=2) :: gkco = 'dk'
+ CHARACTER(len=2) :: abco = 'aa'
+ CHARACTER(len=6) :: coname
+
+ IF (gyrokin_CO) gkco = 'GK'
+ IF (interspecies) abco = 'ab'
+ WRITE(coname,'(a2,a2,a2)') collision_model,gkco,abco
+
+ CALL attach(fidres, TRIM(str), "CO", coname)
+ CALL attach(fidres, TRIM(str), "matfilename", mat_file)
+ END SUBROUTINE coll_outputinputs
+
+ SUBROUTINE compute_TColl
+ USE basic
+ IMPLICIT NONE
+ ! Execution time start
+ CALL cpu_time(t0_coll)
+
+ SELECT CASE(collision_model)
+ CASE ('LB')
+ CALL compute_lenard_bernstein
+ CASE ('DG')
+ CALL compute_dougherty
+ CASE DEFAULT
+ CALL compute_cosolver_coll
+ END SELECT
+
+ ! Execution time end
+ CALL cpu_time(t1_coll)
+ tc_coll = tc_coll + (t1_coll - t0_coll)
+ END SUBROUTINE compute_TColl
+
+ !******************************************************************************!
+ !! Lenard Bernstein collision operator
+ !******************************************************************************!
+ SUBROUTINE compute_lenard_bernstein
+ IMPLICIT NONE
+ COMPLEX(dp) :: TColl_
+ IF (KIN_E) THEN
+ DO ip = ips_e,ipe_e;DO ij = ijs_e,ije_e
+ DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize
+ CALL LenardBernstein_e(ip,ij,ikx,iky,iz,TColl_)
+ TColl_e(ip,ij,ikx,iky,iz) = TColl_
+ ENDDO;ENDDO;ENDDO
+ ENDDO;ENDDO
+ ENDIF
+ DO ip = ips_i,ipe_i;DO ij = ijs_i,ije_i
+ DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize
+ CALL LenardBernstein_i(ip,ij,ikx,iky,iz,TColl_)
+ TColl_i(ip,ij,ikx,iky,iz) = TColl_
+ ENDDO;ENDDO;ENDDO
+ ENDDO;ENDDO
+ END SUBROUTINE compute_lenard_bernstein
+
!******************************************************************************!
- !! Lenard Bernstein collision operator for electrons
+ !! for electrons
!******************************************************************************!
SUBROUTINE LenardBernstein_e(ip_,ij_,ikx_,iky_,iz_,TColl_)
- USE fields, ONLY: moments_e
- USE grid, ONLY: parray_e, jarray_e, kxarray, kyarray
- USE basic
- USE model, ONLY: sigmae2_taue_o2, nu_ee
- USE time_integration, ONLY : updatetlevel
IMPLICIT NONE
INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
@@ -37,7 +134,7 @@ CONTAINS
END SUBROUTINE LenardBernstein_e
!******************************************************************************!
- !! Lenard Bernstein collision operator for electrons
+ !! for ions
!******************************************************************************!
SUBROUTINE LenardBernstein_i(ip_,ij_,ikx_,iky_,iz_,TColl_)
USE fields, ONLY: moments_i
@@ -63,8 +160,27 @@ CONTAINS
END SUBROUTINE LenardBernstein_i
!******************************************************************************!
- !! Doughtery gyrokinetic collision operator for electrons
+ !! Doughtery collision operator
!******************************************************************************!
+ SUBROUTINE compute_dougherty
+ IMPLICIT NONE
+ COMPLEX(dp) :: TColl_
+ IF (KIN_E) THEN
+ DO ip = ips_e,ipe_e;DO ij = ijs_e,ije_e
+ DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize
+ CALL DoughertyGK_e(ip,ij,ikx,iky,iz,TColl_)
+ TColl_e(ip,ij,ikx,iky,iz) = TColl_
+ ENDDO;ENDDO;ENDDO
+ ENDDO;ENDDO
+ ENDIF
+ DO ip = ips_i,ipe_i;DO ij = ijs_i,ije_i
+ DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize
+ CALL DoughertyGK_i(ip,ij,ikx,iky,iz,TColl_)
+ TColl_i(ip,ij,ikx,iky,iz) = TColl_
+ ENDDO;ENDDO;ENDDO
+ ENDDO;ENDDO
+ END SUBROUTINE compute_dougherty
+
SUBROUTINE DoughertyGK_e(ip_,ij_,ikx_,iky_,iz_,TColl_)
USE fields, ONLY: moments_e, phi
USE grid, ONLY: parray_e, jarray_e, kxarray, kyarray, kparray, Jmaxe, ip0_e, ip1_e, ip2_e
@@ -173,12 +289,6 @@ CONTAINS
!! Doughtery gyrokinetic collision operator for ions
!******************************************************************************!
SUBROUTINE DoughertyGK_i(ip_,ij_,ikx_,iky_,iz_,TColl_)
- USE fields, ONLY: moments_i, phi
- USE grid, ONLY: parray_i, jarray_i, kxarray, kyarray, kparray, Jmaxi, ip0_i, ip1_i, ip2_i
- USE array, ONLY: kernel_i
- USE basic
- USE model, ONLY: sigmai2_taui_o2, qi_taui, nu_i, sqrt_sigmai2_taui_o2
- USE time_integration, ONLY : updatetlevel
IMPLICIT NONE
INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
@@ -281,15 +391,7 @@ CONTAINS
!******************************************************************************!
!! compute the collision terms in a (Np x Nj x Nkx x Nky) matrix all at once
!******************************************************************************!
- SUBROUTINE compute_TColl
- USE fields
- USE grid
- USE array
- USE basic
- USE prec_const
- USE time_integration
- USE model
- USE utility
+ SUBROUTINE compute_cosolver_coll
IMPLICIT NONE
COMPLEX(dp), DIMENSION(1:total_np_e) :: local_sum_e, buffer_e, total_sum_e
COMPLEX(dp), DIMENSION(ips_e:ipe_e) :: TColl_distr_e
@@ -297,81 +399,64 @@ CONTAINS
COMPLEX(dp), DIMENSION(ips_i:ipe_i) :: TColl_distr_i
COMPLEX(dp) :: TColl
INTEGER :: ikxs_C, ikxe_C, ikys_C, ikye_C
-
- ! Execution time start
- CALL cpu_time(t0_coll)
-
- IF (ABS(CO) .GE. 2) THEN !compute only if COSOlver matrices are used
- DO ikx = ikxs,ikxe
- DO iky = ikys,ikye
- DO iz = izs,ize
- IF(KIN_E) THEN
- ! Electrons
- DO ij = 1,Jmaxe+1
- ! Loop over all p to compute sub collision term
- DO ip = 1,total_np_e
- CALL apply_COSOlver_mat_e(ip,ij,ikx,iky,iz,TColl)
- local_sum_e(ip) = TColl
- ENDDO
- IF (num_procs_p .GT. 1) THEN
- ! Sum up all the sub collision terms on root 0
- CALL MPI_REDUCE(local_sum_e, buffer_e, total_np_e, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
- ! distribute the sum over the process among p
- CALL MPI_SCATTERV(buffer_e, counts_np_e, displs_np_e, MPI_DOUBLE_COMPLEX,&
- TColl_distr_e, local_np_e, MPI_DOUBLE_COMPLEX,&
- 0, comm_p, ierr)
- ELSE
- TColl_distr_e = local_sum_e
- ENDIF
- ! Write in output variable
- DO ip = ips_e,ipe_e
- TColl_e(ip,ij,ikx,iky,iz) = TColl_distr_e(ip)
- ENDDO
+ DO ikx = ikxs,ikxe
+ DO iky = ikys,ikye
+ DO iz = izs,ize
+ IF(KIN_E) THEN
+ ! Electrons
+ DO ij = 1,Jmaxe+1
+ ! Loop over all p to compute sub collision term
+ DO ip = 1,total_np_e
+ CALL apply_COSOlver_mat_e(ip,ij,ikx,iky,iz,TColl)
+ local_sum_e(ip) = TColl
ENDDO
+ IF (num_procs_p .GT. 1) THEN
+ ! Sum up all the sub collision terms on root 0
+ CALL MPI_REDUCE(local_sum_e, buffer_e, total_np_e, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
+ ! distribute the sum over the process among p
+ CALL MPI_SCATTERV(buffer_e, counts_np_e, displs_np_e, MPI_DOUBLE_COMPLEX,&
+ TColl_distr_e, local_np_e, MPI_DOUBLE_COMPLEX,&
+ 0, comm_p, ierr)
+ ELSE
+ TColl_distr_e = local_sum_e
ENDIF
- ! Ions
- DO ij = 1,Jmaxi+1
- DO ip = 1,total_np_i
- CALL apply_COSOlver_mat_i(ip,ij,ikx,iky,iz,TColl)
- local_sum_i(ip) = TColl
- ENDDO
- IF (num_procs_p .GT. 1) THEN
- ! Reduce the local_sums to root = 0
- CALL MPI_REDUCE(local_sum_i, buffer_i, total_np_i, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
- ! buffer contains the entire collision term along p, we scatter it between
- ! the other processes (use of scatterv since Pmax/Np is not an integer)
- CALL MPI_SCATTERV(buffer_i, counts_np_i, displs_np_i, MPI_DOUBLE_COMPLEX,&
- TColl_distr_i, local_np_i, MPI_DOUBLE_COMPLEX, &
- 0, comm_p, ierr)
- ELSE
- TColl_distr_i = local_sum_i
- ENDIF
- ! Write in output variable
- DO ip = ips_i,ipe_i
- TColl_i(ip,ij,ikx,iky,iz) = TColl_distr_i(ip)
- ENDDO
+ ! Write in output variable
+ DO ip = ips_e,ipe_e
+ TColl_e(ip,ij,ikx,iky,iz) = TColl_distr_e(ip)
+ ENDDO
+ ENDDO
+ ENDIF
+ ! Ions
+ DO ij = 1,Jmaxi+1
+ DO ip = 1,total_np_i
+ CALL apply_COSOlver_mat_i(ip,ij,ikx,iky,iz,TColl)
+ local_sum_i(ip) = TColl
+ ENDDO
+ IF (num_procs_p .GT. 1) THEN
+ ! Reduce the local_sums to root = 0
+ CALL MPI_REDUCE(local_sum_i, buffer_i, total_np_i, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
+ ! buffer contains the entire collision term along p, we scatter it between
+ ! the other processes (use of scatterv since Pmax/Np is not an integer)
+ CALL MPI_SCATTERV(buffer_i, counts_np_i, displs_np_i, MPI_DOUBLE_COMPLEX,&
+ TColl_distr_i, local_np_i, MPI_DOUBLE_COMPLEX, &
+ 0, comm_p, ierr)
+ ELSE
+ TColl_distr_i = local_sum_i
+ ENDIF
+ ! Write in output variable
+ DO ip = ips_i,ipe_i
+ TColl_i(ip,ij,ikx,iky,iz) = TColl_distr_i(ip)
ENDDO
ENDDO
ENDDO
ENDDO
- ENDIF
-
- ! Execution time end
- CALL cpu_time(t1_coll)
- tc_coll = tc_coll + (t1_coll - t0_coll)
- END SUBROUTINE compute_TColl
+ ENDDO
+ END SUBROUTINE compute_cosolver_coll
!******************************************************************************!
!!!!!!! Compute ion collision term
!******************************************************************************!
SUBROUTINE apply_COSOlver_mat_e(ip_,ij_,ikx_,iky_,iz_,TColl_)
- USE fields, ONLY: moments_e, moments_i
- USE grid
- USE array
- USE basic
- USE time_integration, ONLY: updatetlevel
- USE utility
- USE model, ONLY: CO, nu_e, nu_ee, CLOS
IMPLICIT NONE
INTEGER, INTENT(IN) :: ip_, ij_ ,ikx_, iky_, iz_
@@ -380,9 +465,9 @@ CONTAINS
INTEGER :: ip2,ij2, p_int,j_int, p2_int,j2_int, ikx_C, iky_C
p_int = parray_e_full(ip_); j_int = jarray_e_full(ij_);
- IF (CO .GT. 0) THEN ! GK operator (k-dependant)
+ IF (gyrokin_CO) THEN ! GK operator (k-dependant)
ikx_C = ikx_; iky_C = iky_
- ELSEIF (CO .LT. 0) THEN ! DK operator (only one mat for every k)
+ ELSE ! DK operator (only one mat for every k)
ikx_C = 1; iky_C = 1
ENDIF
@@ -394,7 +479,7 @@ CONTAINS
jloopee: DO ij2 = ijs_e,ije_e
j2_int = jarray_e(ij2)
IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxe))&
- TColl_ = TColl_ + moments_e(ip2,ij2,ikx_,iky_,iz_,updatetlevel) &
+ TColl_ = TColl_ + nadiab_moments_e(ip2,ij2,ikx_,iky_,iz_) &
*( nu_e * CeipjT(bare(p_int,j_int), bare(p2_int,j2_int),ikx_C, iky_C) &
+nu_ee * Ceepj (bare(p_int,j_int), bare(p2_int,j2_int),ikx_C, iky_C))
ENDDO jloopee
@@ -406,7 +491,7 @@ CONTAINS
jloopei: DO ij2 = ijs_i,ije_i
j2_int = jarray_i(ij2)
IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxi))&
- TColl_ = TColl_ + moments_i(ip2,ij2,ikx_,iky_,iz_,updatetlevel) &
+ TColl_ = TColl_ + nadiab_moments_i(ip2,ij2,ikx_,iky_,iz_) &
*(nu_e * CeipjF(bare(p_int,j_int), bari(p2_int,j2_int),ikx_C, iky_C))
END DO jloopei
ENDDO ploopei
@@ -417,13 +502,6 @@ CONTAINS
!!!!!!! Compute ion collision term
!******************************************************************************!
SUBROUTINE apply_COSOlver_mat_i(ip_,ij_,ikx_,iky_,iz_,TColl_)
- USE fields, ONLY : moments_e, moments_i
- USE grid
- USE array
- USE basic
- USE time_integration, ONLY : updatetlevel
- USE utility
- USE model, ONLY: CO, nu_i, nu_ie, CLOS, KIN_E
IMPLICIT NONE
INTEGER, INTENT(IN) :: ip_, ij_ ,ikx_, iky_, iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
@@ -431,9 +509,9 @@ CONTAINS
INTEGER :: ip2,ij2, p_int,j_int, p2_int,j2_int, ikx_C, iky_C
p_int = parray_i_full(ip_); j_int = jarray_i_full(ij_);
- IF (CO .GT. 0) THEN ! GK operator (k-dependant)
+ IF (gyrokin_CO) THEN ! GK operator (k-dependant)
ikx_C = ikx_; iky_C = iky_
- ELSEIF (CO .LT. 0) THEN ! DK operator (only one mat for every k)
+ ELSE ! DK operator (only one mat for every k)
ikx_C = 1; iky_C = 1
ENDIF
@@ -445,10 +523,10 @@ CONTAINS
j2_int = jarray_i(ij2)
IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxi))&
! Ion-ion collision
- TColl_ = TColl_ + moments_i(ip2,ij2,ikx_,iky_,iz_,updatetlevel) &
+ TColl_ = TColl_ + nadiab_moments_i(ip2,ij2,ikx_,iky_,iz_) &
* nu_i * Ciipj (bari(p_int,j_int), bari(p2_int,j2_int), ikx_C, iky_C)
IF(KIN_E) & ! Ion-electron collision test
- TColl_ = TColl_ + moments_i(ip2,ij2,ikx_,iky_,iz_,updatetlevel) &
+ TColl_ = TColl_ + nadiab_moments_i(ip2,ij2,ikx_,iky_,iz_) &
* nu_ie * CiepjT(bari(p_int,j_int), bari(p2_int,j2_int), ikx_C, iky_C)
ENDDO jloopii
ENDDO ploopii
@@ -459,26 +537,17 @@ CONTAINS
jloopie: DO ij2 = ijs_e,ije_e
j2_int = jarray_e(ij2)
IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxe))&
- TColl_ = TColl_ + moments_e(ip2,ij2,ikx_,iky_,iz_,updatetlevel) &
+ TColl_ = TColl_ + nadiab_moments_e(ip2,ij2,ikx_,iky_,iz_) &
*(nu_ie * CiepjF(bari(p_int,j_int), bare(p2_int,j2_int), ikx_C, iky_C))
ENDDO jloopie
ENDDO ploopie
ENDIF
-
END SUBROUTINE apply_COSOlver_mat_i
!******************************************************************************!
!!!!!!! Load the collision matrix coefficient table from COSOlver results
!******************************************************************************!
SUBROUTINE load_COSOlver_mat ! Load a sub matrix from iCa files (works for pmaxa,jmaxa<=P_full,J_full)
- use futils
- use initial_par
- USE grid
- USE array, ONLY: Ceepj, Ciipj, CeipjF, CeipjT, CiepjF, CiepjT
- USE basic
- USE time_integration, ONLY : updatetlevel
- USE utility
- USE model, ONLY: CO, NON_LIN, sigmae2_taue_o2, sigmai2_taui_o2
IMPLICIT NONE
! Indices for row and columns of the COSOlver matrix (4D compressed 2D matrices)
INTEGER :: irow_sub, irow_full, icol_sub, icol_full
@@ -507,22 +576,21 @@ CONTAINS
CHARACTER(len=128) :: var_name, kperp_string, ikp_string
- LOGICAL :: CO_AA_ONLY = .false. ! Flag to remove ei ie collision
-
!! Some terminal info
- IF (CO .EQ. 2) THEN
- IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK Sugama matrix ==='
- ELSEIF(CO .EQ. 3) THEN
- IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK pitch angle matrix ==='
- ELSEIF(CO .EQ. 4) THEN
- IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK Coulomb matrix ==='
- ELSEIF(CO .EQ. -2) THEN
- IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK Sugama matrix ==='
- ELSEIF(CO .EQ. -3) THEN
- IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK pitch angle matrix ==='
- ELSEIF(CO .EQ. -4) THEN
- IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK Coulomb matrix ==='
- ENDIF
+ SELECT CASE (collision_model)
+ CASE ('SG')
+ WRITE(*,*) '=== Load Sugama matrix ==='
+ CASE ('LR')
+ IF (my_id .EQ. 0) WRITE(*,*) '=== Load Lorentz matrix ==='
+ CASE ('LD')
+ IF (my_id .EQ. 0) WRITE(*,*) '=== Load Landau matrix ==='
+ END SELECT
+ SELECT CASE (gyrokin_CO)
+ CASE (.true.)
+ IF (my_id .EQ. 0) WRITE(*,*) '..gyrokinetic model..'
+ CASE (.false.)
+ IF (my_id .EQ. 0) WRITE(*,*) '..driftkinetic model..'
+ END SELECT
! Opening the compiled cosolver matrices results
if(my_id.EQ.0)write(*,*) mat_file
@@ -541,16 +609,15 @@ CONTAINS
CALL getarr(fid, '/coordkperp', kp_grid_mat)
! check that we have enough kperps mat
- WRITE(*,*) kp_max
- IF (NON_LIN .GT. 0) THEN
+ IF (LINEARITY .NE. 'linear') THEN
IF ( (kp_grid_mat(nkp_mat) .LT. 2./3.*kp_max) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Matrix kperp grid too small !!'
ELSE
IF ( (kp_grid_mat(nkp_mat) .LT. kp_max) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Matrix kperp grid too small !!'
ENDIF
- IF (CO .GT. 0) THEN ! GK operator (k-dependant)
+ IF (gyrokin_CO) THEN ! GK operator (k-dependant)
ikps_C = 1; ikpe_C = nkp_mat
- ELSEIF (CO .LT. 0) THEN ! DK operator (only one mat for all k)
+ ELSE ! DK operator (only one mat for all k)
ikps_C = 1; ikpe_C = 1
ENDIF
@@ -563,7 +630,7 @@ CONTAINS
DO ikp = ikps_C,ikpe_C ! Loop over everz kperp values
! we put zeros if kp>2/3kpmax because thoses frequencies are filtered through AA
- IF( (kp_grid_mat(ikp) .GT. two_third_kpmax) .AND. (NON_LIN .GT. 0)) THEN
+ IF( (kp_grid_mat(ikp) .GT. two_third_kpmax) .AND. (LINEARITY .NE. 'linear')) THEN
CiepjT_kp(:,:,ikp) = 0._dp
CiepjF_kp(:,:,ikp) = 0._dp
CeipjT_kp(:,:,ikp) = 0._dp
@@ -572,7 +639,7 @@ CONTAINS
Ciipj__kp(:,:,ikp) = 0._dp
ELSE
! Kperp value in string format
- IF (CO .GT. 0) THEN
+ IF (gyrokin_CO) THEN
write(ikp_string,'(i5.5)') ikp-1
ELSE
write(ikp_string,'(i5.5)') 0
@@ -621,7 +688,7 @@ CONTAINS
ENDDO
DEALLOCATE(Ciipj_full)
- IF(abs(CO) .NE. 3) THEN ! Pitch angle is only applied on like-species
+ IF(interspecies) THEN ! Pitch angle is only applied on like-species
!!!!!!!!!!!!!!! E-I matrices !!!!!!!!!!!!!!
! Get test and field e-i collision matrices
CALL allocate_array( CeipjT_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1))
@@ -692,13 +759,13 @@ CONTAINS
ENDDO
CALL closef(fid)
- IF (CO .GT. 0) THEN ! Interpolation of the kperp matrix values on kx ky grid
+ IF (gyrokin_CO) THEN ! Interpolation of the kperp matrix values on kx ky grid
IF (my_id .EQ. 0 ) WRITE(*,*) '...Interpolation from matrices kperp to simulation kx,ky...'
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
! Check for nonlinear case if we are in the anti aliased domain or the filtered one
kperp_sim = kparray(ikx,iky,1,0) ! current simulation kperp
- IF( (NON_LIN .EQ. 0) .OR. (kperp_sim .LT. two_third_kpmax) ) THEN
+ IF( (LINEARITY .EQ. 'linear') .OR. (kperp_sim .LT. two_third_kpmax) ) THEN
! Find the interval in kp grid mat where kperp_sim is contained
! Loop over the whole kp mat grid to find the smallest kperp that is
! larger than the current kperp_sim (brute force...)
@@ -711,20 +778,26 @@ CONTAINS
! interval boundaries
ikp_next = ikp_mat !index of k1 (larger than kperp_sim thanks to previous loop)
ikp_prev = ikp_mat - 1 !index of k0 (smaller neighbour to interpolate inbetween)
- ! write(*,*) kp_grid_mat(ikp_prev), '<', kperp_sim, '<', kp_grid_mat(ikp_next)
if ( (kp_grid_mat(ikp_prev) .GT. kperp_sim) .OR. (kp_grid_mat(ikp_next) .LT. kperp_sim) ) THEN
- write(*,*) 'Warning, linear interp of collision matrix failed!! '
+ ! write(*,*) 'Warning, linear interp of collision matrix failed!! '
+ ! write(*,*) kp_grid_mat(ikp_prev), '<', kperp_sim, '<', kp_grid_mat(ikp_next)
ENDIF
! 0->1 variable for linear interp, i.e. zero2one = (k-k0)/(k1-k0)
zerotoone = (kperp_sim - kp_grid_mat(ikp_prev))/(kp_grid_mat(ikp_next) - kp_grid_mat(ikp_prev))
-
! Linear interpolation between previous and next kperp matrix values
Ceepj (:,:,ikx,iky) = (Ceepj__kp(:,:,ikp_next) - Ceepj__kp(:,:,ikp_prev))*zerotoone + Ceepj__kp(:,:,ikp_prev)
- CeipjT(:,:,ikx,iky) = (CeipjT_kp(:,:,ikp_next) - CeipjT_kp(:,:,ikp_prev))*zerotoone + CeipjT_kp(:,:,ikp_prev)
- CeipjF(:,:,ikx,iky) = (CeipjF_kp(:,:,ikp_next) - CeipjF_kp(:,:,ikp_prev))*zerotoone + CeipjF_kp(:,:,ikp_prev)
Ciipj (:,:,ikx,iky) = (Ciipj__kp(:,:,ikp_next) - Ciipj__kp(:,:,ikp_prev))*zerotoone + Ciipj__kp(:,:,ikp_prev)
- CiepjT(:,:,ikx,iky) = (CiepjT_kp(:,:,ikp_next) - CiepjT_kp(:,:,ikp_prev))*zerotoone + CiepjT_kp(:,:,ikp_prev)
- CiepjF(:,:,ikx,iky) = (CiepjF_kp(:,:,ikp_next) - CiepjF_kp(:,:,ikp_prev))*zerotoone + CiepjF_kp(:,:,ikp_prev)
+ IF(interspecies) THEN
+ CeipjT(:,:,ikx,iky) = (CeipjT_kp(:,:,ikp_next) - CeipjT_kp(:,:,ikp_prev))*zerotoone + CeipjT_kp(:,:,ikp_prev)
+ CeipjF(:,:,ikx,iky) = (CeipjF_kp(:,:,ikp_next) - CeipjF_kp(:,:,ikp_prev))*zerotoone + CeipjF_kp(:,:,ikp_prev)
+ CiepjT(:,:,ikx,iky) = (CiepjT_kp(:,:,ikp_next) - CiepjT_kp(:,:,ikp_prev))*zerotoone + CiepjT_kp(:,:,ikp_prev)
+ CiepjF(:,:,ikx,iky) = (CiepjF_kp(:,:,ikp_next) - CiepjF_kp(:,:,ikp_prev))*zerotoone + CiepjF_kp(:,:,ikp_prev)
+ ELSE
+ CeipjT(:,:,ikx,iky) = 0._dp
+ CeipjF(:,:,ikx,iky) = 0._dp
+ CiepjT(:,:,ikx,iky) = 0._dp
+ CiepjF(:,:,ikx,iky) = 0._dp
+ ENDIF
ELSE ! Out of anti anti aliased domain
Ceepj (:,:,ikx,iky) = 0._dp
CeipjT(:,:,ikx,iky) = 0._dp
@@ -747,7 +820,7 @@ CONTAINS
DEALLOCATE (Ceepj__kp); DEALLOCATE (CeipjT_kp); DEALLOCATE (CeipjF_kp)
DEALLOCATE (Ciipj__kp); DEALLOCATE (CiepjT_kp); DEALLOCATE (CiepjF_kp)
- IF( CO_AA_ONLY ) THEN
+ IF( .NOT. interspecies ) THEN
CeipjF = 0._dp;
CeipjT = 0._dp;
CiepjF = 0._dp;