diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 383f85a537a733c189bfd46440d45f8969586440..b0594ebab94efdd2710e8feb9a90a3a58f2e6a59 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -8,6 +8,10 @@ MODULE array
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_rhs_e ! (ip,ij,ikr,ikz,updatetlevel)
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_rhs_i ! (ip,ij,ikr,ikz,updatetlevel)
+ ! To load collision matrix
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: iCe ! (ip,ij)
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: iCi ! (ip,ij)
+
! Intermediate steps in rhs of equations
!COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE:: moments_Apl, moments_Bpl, moments_Cpl, moments_Dpl,&
! moments_Epl, moments_Fpl, moments_Gpl, moments_Hpl
diff --git a/src/inital.F90 b/src/inital.F90
index 77755599494cd97f14356c9dd649426f8628bda5..c60072f615117f5bd3958b5c56c3aa26a3b185b2 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -1,14 +1,21 @@
SUBROUTINE inital
- use basic
-
- use prec_const
+ USE basic
+ USE model, ONLY : CO
+ USE prec_const
implicit none
WRITE(*,'(a/)') '=== Set initial conditions ==='
CALL init_profiles
+ IF (CO .EQ. -1) THEN
+ WRITE(*,'(a/)') '=== Load Full Coulomb matrix ==='
+
+ CALL load_FC_mat
+ ENDIF
+ !
+
END SUBROUTINE inital
@@ -18,10 +25,10 @@ SUBROUTINE init_profiles
USE basic
USE fourier_grid
USE fields
- use initial_par
+ USE initial_par
USE time_integration
- use prec_const
+ USE prec_const
IMPLICIT NONE
INTEGER :: ip,ij, ikr,ikz
@@ -62,3 +69,9 @@ SUBROUTINE init_profiles
END SUBROUTINE init_profiles
+SUBROUTINE load_FC_mat
+ USE basic
+ USE array
+ IMPLICIT NONE
+
+END SUBROUTINE load_FC_mat
\ No newline at end of file
diff --git a/src/memory.F90 b/src/memory.F90
index af5933d0be418d4640d608491677f041311ef3fd..410abfdc6004c8904666c20901f89ef490627180 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -6,8 +6,9 @@ SUBROUTINE memory
USE fields
USE fourier_grid
USE time_integration
+ USE model, ONLY: CO
- use prec_const
+ USE prec_const
IMPLICIT NONE
! Moments and moments rhs
@@ -19,4 +20,9 @@ SUBROUTINE memory
! Electrostatic potential
CALL allocate_array(phi, ikrs,ikre, ikzs,ikze)
+ ! Collision matrix
+ IF (CO .EQ. -1) THEN
+ CALL allocate_array( iCe, ips_i,ipe_i, ijs_i,ije_i )
+ CALL allocate_array( iCi, ips_i,ipe_i, ijs_i,ije_i )
+ ENDIF
END SUBROUTINE memory
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index fe657e64cd89065f8819fb074be4feb42369dbd0..f30447b0b3eb49652bf995e4cb34f1ee148a3b58 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -5,17 +5,18 @@ MODULE model
IMPLICIT NONE
PRIVATE
- REAL(dp), PUBLIC, PROTECTED :: nu = 1._dp ! Collision frequency
- REAL(dp), PUBLIC, PROTECTED :: tau_e = 1._dp ! Collision frequency
- REAL(dp), PUBLIC, PROTECTED :: tau_i = 1._dp ! Collision frequency
- REAL(dp), PUBLIC, PROTECTED :: sigma_e = 1._dp ! Collision frequency
- REAL(dp), PUBLIC, PROTECTED :: sigma_i = 1._dp ! Collision frequency
- REAL(dp), PUBLIC, PROTECTED :: q_e = 1._dp ! Collision frequency
- REAL(dp), PUBLIC, PROTECTED :: q_i = 1._dp ! Collision frequency
- REAL(dp), PUBLIC, PROTECTED :: eta_n = 1._dp ! Collision frequency
- REAL(dp), PUBLIC, PROTECTED :: eta_T = 1._dp ! Collision frequency
- REAL(dp), PUBLIC, PROTECTED :: eta_B = 1._dp ! Collision frequency
- REAL(dp), PUBLIC, PROTECTED :: lambdaD = 1._dp ! Collision frequency
+ INTEGER(dp), PUBLIC, PROTECTED :: CO = 0 ! Collision Operator
+ REAL(dp), PUBLIC, PROTECTED :: nu = 1._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_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 :: lambdaD = 1._dp ! Debye length
PUBLIC :: model_readinputs, model_outputinputs
@@ -28,7 +29,7 @@ CONTAINS
USE prec_const
IMPLICIT NONE
- NAMELIST /MODEL_PAR/ nu, tau_e, tau_i, sigma_e, sigma_i, &
+ NAMELIST /MODEL_PAR/ CO, nu, tau_e, tau_i, sigma_e, sigma_i, &
q_e, q_i, eta_n, eta_T, eta_B, lambdaD
READ(lu_in,model_par)
@@ -49,6 +50,7 @@ CONTAINS
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
+ CALL attach(fidres, TRIM(str), "CO", CO)
CALL attach(fidres, TRIM(str), "nu", nu)
CALL attach(fidres, TRIM(str), "tau_e", tau_e)
CALL attach(fidres, TRIM(str), "tau_i", tau_i)
diff --git a/wk/kz_linear.m b/wk/kz_linear.m
index 6ed8cc1280a58433fe933f470d1242b137eaa9d1..0c0036d220d2b89df001ead22f62943cc8e6c38f 100644
--- a/wk/kz_linear.m
+++ b/wk/kz_linear.m
@@ -13,18 +13,19 @@ OUTPUTS.write_phi = '.true.';
OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results''';
%% Grid parameters
-GRID.pmaxe = 81;
-GRID.jmaxe = 20;
-GRID.pmaxi = 81;
-GRID.jmaxi = 20;
+GRID.pmaxe = 25;
+GRID.jmaxe = 9;
+GRID.pmaxi = 25;
+GRID.jmaxi = 9;
GRID.nkr = 1;
GRID.krmin = 0.;
GRID.krmax = 0.;
-GRID.nkz = 10;
+GRID.nkz = 20;
GRID.kzmin = 0.1/sqrt(2.0);
-GRID.kzmax = 2.8/sqrt(2.0);
+GRID.kzmax = 5.0/sqrt(2.0);
%% Model parameters
-MODEL.nu = 0.001;
+MODEL.CO = -1; % Collision operator (-1 = Full Coulomb, 0 = Dougherty)
+MODEL.nu = 0.001; % collisionality nu*L_perp/Cs0
MODEL.tau_e = 1.0;
MODEL.tau_i = 1.0;
MODEL.sigma_e = 0.0233380;