diff --git a/src/array_mod.F90 b/src/array_mod.F90
index c22e586a2c638ae40d180b4621d1f5461aaf9c11..c401b0430630e3e2e9f7004e4d9fdcef903d52df 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -17,8 +17,8 @@ MODULE array
   COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: dnjs
 
   ! Kernel function evaluation
-  COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: kernel_e
-  COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: kernel_i
+  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: kernel_e
+  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: kernel_i
 
   ! Non linear term array (ip,ij,ikr,ikz)
   COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Sepj ! electron
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index 8019271b11001ffdaa50e5160a54101ab22d8b6d..7cbf0affb3c295c8c6e4d20256796bc1c9086f71 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -8,13 +8,30 @@ use grid
 use basic
 use futils
 use initial_par
+use model
 
 implicit none
 
-PUBLIC :: load_FC_mat, load_FC_mat_old
+PUBLIC :: load_FC_mat
 
 CONTAINS
 
+!******************************************************************************!
+SUBROUTINE LenardBernsteinDK
+
+END SUBROUTINE LenardBernsteinDK
+
+!******************************************************************************!
+SUBROUTINE DoughertyDK
+
+END SUBROUTINE DoughertyDK
+
+!******************************************************************************!
+SUBROUTINE FullCoulombDK
+
+END SUBROUTINE FullCoulombDK
+
+
 !******************************************************************************!
 !!!!!!! Load the Full coulomb coefficient table from COSOlver results
 !******************************************************************************!
@@ -63,7 +80,6 @@ SUBROUTINE load_FC_mat ! Load a sub matrix from iCa files (works for pmaxa,jmaxa
 
   CALL closef(fid1)
   DEALLOCATE(Ceepj_full)
-  IF (my_id .EQ. 0) WRITE(*,*) '..done'
 
   IF (my_id .EQ. 0) WRITE(*,*) 'Load ei FC mat...'
   ! get the Test and Back field electron ion collision matrix
@@ -104,7 +120,6 @@ SUBROUTINE load_FC_mat ! Load a sub matrix from iCa files (works for pmaxa,jmaxa
   CALL closef(fid2)
   DEALLOCATE(CeipjF_full)
   DEALLOCATE(CeipjT_full)
-  IF (my_id .EQ. 0) WRITE(*,*) '..done'
 
   !!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!!
   IF (my_id .EQ. 0) WRITE(*,*) 'Load ii FC mat...'
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 6c6ab63fe803a3eb60d4f60ba76b9ff73db19149..4d5757cdb852aa8ae51de1b416cbca21aa7b419e 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -210,13 +210,13 @@ end
 
 if 0
 %% Photomaton : real space
-FIELD = ni00; FNAME = 'ni';
+% FIELD = ni00; FNAME = 'ni';
 % FIELD = ne00; FNAME = 'ne';
-% FIELD = phi; FNAME = 'phi';
-tf = 19;  [~,it1] = min(abs(Ts2D-tf));
-tf = 20;  [~,it2] = min(abs(Ts2D-tf)); 
-tf = 21; [~,it3] = min(abs(Ts2D-tf));
-tf = 22; [~,it4] = min(abs(Ts2D-tf));
+FIELD = phi; FNAME = 'phi';
+tf = 60;  [~,it1] = min(abs(Ts2D-tf));
+tf = 120;  [~,it2] = min(abs(Ts2D-tf)); 
+tf = 250; [~,it3] = min(abs(Ts2D-tf));
+tf = 400; [~,it4] = min(abs(Ts2D-tf));
 fig = figure; FIGNAME = [FNAME,'_snaps']; set(gcf, 'Position',  [100, 100, 1500, 400])
 plt = @(x) x;%./max(max(x));
     subplot(141)
diff --git a/wk/debug_script.m b/wk/debug_script.m
index 5aa5d036a7b0db39db3724f8273e24c298db211d..ea5d3cae7cb76ca2b6ddb639be87e03274860077 100644
--- a/wk/debug_script.m
+++ b/wk/debug_script.m
@@ -5,9 +5,9 @@ addpath(genpath('../matlab')) % ... add
 CLUSTER.TIME  = '24:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0e2;   % Collision frequency
+NU      = 1e1;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-ETAB    = 0.4;    % Magnetic gradient
+ETAB    = 0.5;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
 HD_CO   = 0.5;    % Hyper diffusivity cutoff ratio
@@ -20,8 +20,8 @@ JMAXE   = 1;     % Highest ''       Laguerre ''
 PMAXI   = 2;     % Highest ion      Hermite polynomial degree
 JMAXI   = 1;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS
-TMAX    = 100;  % Maximal time unit
-DT      = 1e-2;   % Time step
+TMAX    = 400;  % Maximal time unit
+DT      = 3e-2;   % Time step
 SPS0D   = 1/DT;    % Sampling per time unit for profiler
 SPS2D   = 2;      % Sampling per time unit for 2D arrays
 SPS5D   = 1/5;    % Sampling per time unit for 5D arrays
@@ -29,7 +29,7 @@ SPSCP   = 1/10;    % Sampling per time unit for checkpoints
 RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 0;
 %% OPTIONS
-SIMID   = 'test_load_FC';  % Name of the simulation
+SIMID   = 'test_precomp_kernel';  % Name of the simulation
 CO      = -1;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/wk/linear_study.m b/wk/linear_study.m
index b3d0d5caa70195fb4c6077b00b077ce7e0e45cf9..2fd3194361088b16d34211fb2795e073cef605b5 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -5,9 +5,9 @@ default_plots_options
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1e-1;   % Collision frequency
+NU      = 1e+0;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-ETAB    = 0.66;
+ETAB    = 0.5;
 ETAN    = 1.0;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
 MU      = 1e-3;   % Hyper diffusivity coefficient
@@ -17,9 +17,9 @@ NOISE0  = 1.0e-5;
 N       = 100;     % Frequency gridpoints (Nkr = N/2)
 L       = 100;     % Size of the squared frequency domain
 KREQ0   = 1;      % put kr = 0
-%% TIME PARAMETERS
+%% TIME PARMETERS
 TMAX    = 100;  % Maximal time unit
-DT      = 1e-2;   % Time step
+DT      = 5e-2;   % Time step
 SPS0D   = 0.5;      % Sampling per time unit for 2D arrays
 SPS2D   = 1;      % Sampling per time unit for 2D arrays
 SPS5D   = 0.1;    % Sampling per time unit for 5D arrays
@@ -41,12 +41,12 @@ KPAR    = 0.0;    % Parellel wave vector component
 %% PARAMETER SCANS
 if 1
 %% Parameter scan over PJ
-PA = [2, 6, 8, 12];
-JA = [1, 3, 4,  6];
-% PA = [2];
-% JA = [1];
+% PA = [2, 6, 8, 12];
+% JA = [1, 3, 4,  6];
+PA = [4];
+JA = [2];
 Nparam = numel(PA);
-
+param_name = 'PJ';
 gamma_Ni = zeros(Nparam,N/2+1);
 Ni00_ST  = zeros(Nparam,N/2+1,TMAX);
 for i = 1:Nparam
@@ -56,7 +56,7 @@ for i = 1:Nparam
     setup
     % Run linear simulation
     system(...
-        ['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 8 ./../../../bin/helaz; cd ../../../wk']...
+        ['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ./../../../bin/helaz; cd ../../../wk']...
     )
     % Load and process results
     load_results