diff --git a/Makefile b/Makefile
index 0e18cad..456be93 100644
--- a/Makefile
+++ b/Makefile
@@ -1,138 +1,141 @@
 include local/dirs.inc
 include local/make.inc
 
 EXEC = $(BINDIR)/helaz
 
 F90 = mpiifort
 
 # Add Multiple-Precision Library
 EXTLIBS += -L$(FMDIR)/lib
 EXTINC += -I$(FMDIR)/mod
 
 EXTLIBS += -L$(FFTWDIR)/lib
 EXTINC += -I$(FFTWDIR)/include
 
 all: dirs src/srcinfo.h $(EXEC)
 
 run: all
 	(cd wk; $(EXEC);)
 
 dirs:
 	mkdir -p $(BINDIR)
 	mkdir -p $(OBJDIR)
 	mkdir -p $(MODDIR)
 
 src/srcinfo.h:
 	( cd src/srcinfo; $(MAKE);)
 
 clean: cleanobj cleanmod
 	@rm -f src/srcinfo.h
 	@rm -f src/srcinfo/srcinfo.h
 
 cleanobj:
 	@rm -f $(OBJDIR)/*o
 
 cleanmod:
 	@rm -f $(MODDIR)/*mod
 	@rm -f *.mod
 
 cleanbin:
 	@rm -f $(EXEC)
 
 $(OBJDIR)/diagnose.o : src/srcinfo.h
 
 FOBJ=$(OBJDIR)/advance_field.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o $(OBJDIR)/basic_mod.o \
-$(OBJDIR)/coeff_mod.o $(OBJDIR)/compute_Sapj.o $(OBJDIR)/control.o $(OBJDIR)/fourier_mod.o \
+$(OBJDIR)/coeff_mod.o $(OBJDIR)/collision_mod.o $(OBJDIR)/compute_Sapj.o $(OBJDIR)/control.o $(OBJDIR)/fourier_mod.o \
 $(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o $(OBJDIR)/fields_mod.o \
 $(OBJDIR)/inital.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/main.o $(OBJDIR)/memory.o \
 $(OBJDIR)/model_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o \
 $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/readinputs.o \
 $(OBJDIR)/grid_mod.o $(OBJDIR)/stepon.o $(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o \
 $(OBJDIR)/utility_mod.o
 
  $(EXEC): $(FOBJ)
 	$(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@
 
  $(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/advance_field.F90 -o $@
 
  $(OBJDIR)/array_mod.o : src/array_mod.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/array_mod.F90 -o $@
 
  $(OBJDIR)/auxval.o : src/auxval.F90 $(OBJDIR)/fourier_mod.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o  $(OBJDIR)/grid_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/auxval.F90 -o $@
 
  $(OBJDIR)/basic_mod.o : src/basic_mod.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/basic_mod.F90 -o $@
 
  $(OBJDIR)/coeff_mod.o : src/coeff_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o
 		$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/coeff_mod.F90 -o $@
 
+ $(OBJDIR)/collision_mod.o : src/collision_mod.F90 $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o
+	 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/collision_mod.F90 -o $@
+
  $(OBJDIR)/compute_Sapj.o : src/compute_Sapj.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
 		$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/compute_Sapj.F90 -o $@
 
  $(OBJDIR)/control.o : src/control.F90 $(OBJDIR)/auxval.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/readinputs.o $(OBJDIR)/tesend.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/control.F90 -o $@
 
  $(OBJDIR)/fourier_mod.o : src/fourier_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_mod.F90 -o $@
 
  $(OBJDIR)/diagnose.o : src/diagnose.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/diagnose.F90 -o $@
 
  $(OBJDIR)/diagnostics_par_mod.o : src/diagnostics_par_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/diagnostics_par_mod.F90 -o $@
 
  $(OBJDIR)/endrun.o : src/endrun.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/endrun.F90 -o $@
 
  $(OBJDIR)/fields_mod.o : src/fields_mod.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fields_mod.F90 -o $@
 
  $(OBJDIR)/grid_mod.o : src/grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o
 		$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/grid_mod.F90 -o $@
 
  $(OBJDIR)/inital.o : src/inital.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/inital.F90 -o $@
 
  $(OBJDIR)/initial_par_mod.o : src/initial_par_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/initial_par_mod.F90 -o $@
 
  $(OBJDIR)/main.o : src/main.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/main.F90 -o $@
 
  $(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/grid_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/memory.F90 -o $@
 
  $(OBJDIR)/model_mod.o : src/model_mod.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/model_mod.F90 -o $@
 
  $(OBJDIR)/moments_eq_rhs.o : src/moments_eq_rhs.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs.F90 -o $@
 
  $(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/poisson.F90 -o $@
 
  $(OBJDIR)/ppexit.o : src/ppexit.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppexit.F90 -o $@
 
  $(OBJDIR)/ppinit.o : src/ppinit.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppinit.F90 -o $@
 
  $(OBJDIR)/prec_const_mod.o : src/prec_const_mod.F90
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/prec_const_mod.F90 -o $@
 
  $(OBJDIR)/readinputs.o : src/readinputs.F90  $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/readinputs.F90 -o $@
 
  $(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/stepon.F90 -o $@
 
  $(OBJDIR)/tesend.o : src/tesend.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/tesend.F90 -o $@
 
  $(OBJDIR)/time_integration_mod.o : src/time_integration_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/time_integration_mod.F90 -o $@
 
  $(OBJDIR)/utility_mod.o : src/utility_mod.F90  $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
new file mode 100644
index 0000000..ff27b69
--- /dev/null
+++ b/src/collision_mod.F90
@@ -0,0 +1,219 @@
+module collision
+! contains the Hermite-Laguerre collision operators. Solved using COSOlver.
+
+use prec_const
+use fields
+use array
+use grid
+use basic
+use futils
+use initial_par
+
+implicit none
+
+PUBLIC :: load_FC_mat, load_FC_mat_old
+
+CONTAINS
+
+!******************************************************************************!
+!!!!!!! Load the Full coulomb coefficient table from COSOlver results
+!******************************************************************************!
+SUBROUTINE load_FC_mat ! Works only for a partiular file for now with P,J <= 15,6
+  IMPLICIT NONE
+
+  INTEGER :: irow_sub, irow_full, icol_sub, icol_full
+  INTEGER :: fid1, fid2, fid3, fid4
+
+  INTEGER :: ip_e, ij_e, il_e, ik_e
+  INTEGER :: pdime, jdime
+  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ceepj_full, CeipjT_full
+  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CeipjF_full
+  INTEGER :: ip_i, ij_i, il_i, ik_i
+  INTEGER :: pdimi, jdimi
+  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ciipj_full, CiepjT_full
+  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CiepjF_full
+
+  !!!!!!!!!!!! Electron matrices !!!!!!!!!!!!
+  IF (my_id .EQ. 0) WRITE(*,*) 'Load ee FC mat...'
+  ! get the self electron colision matrix
+  CALL openf(selfmat_file,fid1, 'r', 'D',mpicomm=MPI_COMM_WORLD)
+
+  CALL getatt(fid1,'/Caapj/Ceepj/','Pmaxe',pdime)
+  CALL getatt(fid1,'/Caapj/Ceepj/','Jmaxe',jdime)
+
+  IF ( ((pdime .LT. pmaxe) .OR. (jdime .LT. jmaxe)) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! FC Matrix too small !!'
+
+  CALL allocate_array(  Ceepj_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1))
+  CALL getarr(fid1, '/Caapj/Ceepj', Ceepj_full) ! get array (moli format)
+
+  ! Fill sub array with only usefull polynmials degree
+  DO ip_e = 0,pmaxe ! Loop over rows
+  DO ij_e = 0,jmaxe
+        irow_sub  = (jmaxe +1)*ip_e + ij_e +1
+        irow_full = (jdime +1)*ip_e + ij_e +1
+        DO il_e = 0,pmaxe ! Loop over columns
+        DO ik_e = 0,jmaxe
+              icol_sub  = (jmaxe +1)*il_e + ik_e +1
+              icol_full = (jdime +1)*il_e + ik_e +1
+              Ceepj (irow_sub,icol_sub) = Ceepj_full (irow_full,icol_full)
+        ENDDO
+        ENDDO
+  ENDDO
+  ENDDO
+
+  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
+  CALL openf(eimat_file,fid2, 'r', 'D');
+
+  CALL getatt(fid2,'/Ceipj/CeipjT/','Pmaxi',pdimi)
+  CALL getatt(fid2,'/Ceipj/CeipjT/','Jmaxi',jdimi)
+  IF ( (pdimi .LT. pmaxi) .OR. (jdimi .LT. jmaxi) ) WRITE(*,*) '!! FC Matrix too small !!'
+
+  CALL allocate_array( CeipjT_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1))
+  CALL allocate_array( CeipjF_full, 1,(pdime+1)*(jdime+1), 1,(pdimi+1)*(jdimi+1))
+
+  CALL getarr(fid2, '/Ceipj/CeipjT', CeipjT_full)
+  CALL getarr(fid2, '/Ceipj/CeipjF', CeipjF_full)
+
+  ! Fill sub array with only usefull polynmials degree
+  DO ip_e = 0,pmaxe ! Loop over rows
+  DO ij_e = 0,jmaxe
+        irow_sub  = (jmaxe +1)*ip_e + ij_e +1
+        irow_full = (jdime +1)*ip_e + ij_e +1
+        DO il_e = 0,pmaxe ! Loop over columns
+        DO ik_e = 0,jmaxe
+              icol_sub  = (jmaxe +1)*il_e + ik_e +1
+              icol_full = (jdime +1)*il_e + ik_e +1
+              CeipjT(irow_sub,icol_sub) = CeipjT_full(irow_full,icol_full)
+        ENDDO
+        ENDDO
+        DO il_i = 0,pmaxi ! Loop over columns
+        DO ik_i = 0,jmaxi
+              icol_sub  = (Jmaxi +1)*il_i + ik_i +1
+              icol_full = (jdimi +1)*il_i + ik_i +1
+              CeipjF(irow_sub,icol_sub) = CeipjF_full(irow_full,icol_full)
+        ENDDO
+        ENDDO
+  ENDDO
+  ENDDO
+
+  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...'
+  ! get the self electron colision matrix
+  CALL openf(selfmat_file, fid3, 'r', 'D',mpicomm=MPI_COMM_WORLD);
+
+  CALL allocate_array(  Ciipj_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1))
+
+  IF ( (pmaxe .EQ. pmaxi) .AND. (jmaxe .EQ. jmaxi) ) THEN ! if same degrees ion and electron matrices are alike so load Ceepj
+    CALL getarr(fid3, '/Caapj/Ceepj', Ciipj_full) ! get array (moli format)
+  ELSE
+    CALL getarr(fid3, '/Caapj/Ciipj', Ciipj_full) ! get array (moli format)
+  ENDIF
+
+  ! Fill sub array with only usefull polynmials degree
+  DO ip_i = 0,Pmaxi ! Loop over rows
+  DO ij_i = 0,Jmaxi
+        irow_sub  = (Jmaxi +1)*ip_i + ij_i +1
+        irow_full = (jdimi +1)*ip_i + ij_i +1
+        DO il_i = 0,Pmaxi ! Loop over columns
+        DO ik_i = 0,Jmaxi
+              icol_sub  = (Jmaxi +1)*il_i + ik_i +1
+              icol_full = (jdimi +1)*il_i + ik_i +1
+              Ciipj (irow_sub,icol_sub) = Ciipj_full (irow_full,icol_full)
+        ENDDO
+        ENDDO
+  ENDDO
+  ENDDO
+
+  CALL closef(fid3)
+  DEALLOCATE(Ciipj_full)
+
+  ! get the Test and Back field electron ion collision matrix
+  CALL openf(iemat_file,fid4, 'r', 'D');
+
+  CALL allocate_array( CiepjT_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1))
+  CALL allocate_array( CiepjF_full, 1,(pdimi+1)*(jdimi+1), 1,(pdime+1)*(jdime+1))
+
+  IF (my_id .EQ. 0) WRITE(*,*) 'Load ie FC mat...'
+  CALL getarr(fid4, '/Ciepj/CiepjT', CiepjT_full)
+  CALL getarr(fid4, '/Ciepj/CiepjF', CiepjF_full)
+
+  ! Fill sub array with only usefull polynmials degree
+  DO ip_i = 0,Pmaxi ! Loop over rows
+  DO ij_i = 0,Jmaxi
+        irow_sub  = (Jmaxi +1)*ip_i + ij_i +1
+        irow_full = (jdimi +1)*ip_i + ij_i +1
+        DO il_i = 0,Pmaxi ! Loop over columns
+        DO ik_i = 0,Jmaxi
+              icol_sub  = (Jmaxi +1)*il_i + ik_i +1
+              icol_full = (jdimi +1)*il_i + ik_i +1
+              CiepjT(irow_sub,icol_sub) = CiepjT_full(irow_full,icol_full)
+        ENDDO
+        ENDDO
+        DO il_e = 0,pmaxe ! Loop over columns
+        DO ik_e = 0,jmaxe
+              icol_sub  = (jmaxe +1)*il_e + ik_e +1
+              icol_full = (jdime +1)*il_e + ik_e +1
+              CiepjF(irow_sub,icol_sub) = CiepjF_full(irow_full,icol_full)
+        ENDDO
+        ENDDO
+  ENDDO
+  ENDDO
+
+  CALL closef(fid4)
+  DEALLOCATE(CiepjF_full)
+  DEALLOCATE(CiepjT_full)
+
+END SUBROUTINE load_FC_mat
+!******************************************************************************!
+
+SUBROUTINE load_FC_mat_old
+	  USE basic
+	  USE grid
+	  USE initial_par
+	  USE array
+	  use model
+	  USE futils, ONLY : openf, getarr, closef
+	  IMPLICIT NONE
+
+	  INTEGER :: ip2,ij2
+	  INTEGER :: fid1, fid2, fid3, fid4
+
+	  !!!!!!!!!!!! Electron matrices !!!!!!!!!!!!
+	  ! get the self electron colision matrix
+	  CALL openf(selfmat_file,fid1, 'r', 'D',mpicomm=MPI_COMM_WORLD);
+	  CALL getarr(fid1, '/Caapj/Ceepj', Ceepj) ! get array (moli format)
+	  CALL closef(fid1)
+	  ! get the Test and Back field electron ion collision matrix
+	  CALL openf(eimat_file,fid2, 'r', 'D');
+	  CALL getarr(fid2, '/Ceipj/CeipjT', CeipjT)
+	  CALL getarr(fid2, '/Ceipj/CeipjF', CeipjF)
+	  CALL closef(fid2)
+
+	  !!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!!
+	  ! get the self electron colision matrix
+	  CALL openf(selfmat_file, fid3, 'r', 'D',mpicomm=MPI_COMM_WORLD);
+	  IF ( (pmaxe .EQ. pmaxi) .AND. (jmaxe .EQ. jmaxi) ) THEN ! if same degrees ion and electron matrices are alike so load Ceepj
+	    CALL getarr(fid3, '/Caapj/Ceepj', Ciipj) ! get array (moli format)
+	  ELSE
+	    CALL getarr(fid3, '/Caapj/Ciipj', Ciipj) ! get array (moli format)
+	  ENDIF
+	  CALL closef(fid3)
+	  ! get the Test and Back field ion electron collision matrix
+	  CALL openf(iemat_file, fid4, 'r', 'D');
+	  CALL getarr(fid4, '/Ciepj/CiepjT', CiepjT)
+	  CALL getarr(fid4, '/Ciepj/CiepjF', CiepjF)
+	  CALL closef(fid4)
+	END SUBROUTINE load_FC_mat_old
+	!******************************************************************************!
+
+end module collision
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index a6cefad..7d1dfb6 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -1,189 +1,189 @@
 SUBROUTINE compute_Sapj
   ! This routine is meant to compute the non linear term for each specie and degree
   !! In real space Sapj ~ b*(grad(phi) x grad(g)) which in moments in fourier becomes
   !! Sapj = Sum_n (ikr Kn phi)#(ikz Sum_s d_njs Naps) - (ikz Kn phi)#(ikr Sum_s d_njs Naps)
   !! where # denotes the convolution.
   USE array, ONLY : dnjs, Sepj, Sipj
   USE basic
   USE fourier
   USE fields!, ONLY : phi, moments_e, moments_i
   USE grid
   USE model
   USE prec_const
   USE time_integration!, ONLY : updatetlevel
   IMPLICIT NONE
   INCLUDE 'fftw3-mpi.f03'
 
   COMPLEX(dp), DIMENSION(ikrs:ikre,ikzs:ikze) :: Fr_cmpx, Gz_cmpx
   COMPLEX(dp), DIMENSION(ikrs:ikre,ikzs:ikze) :: Fz_cmpx, Gr_cmpx, F_conv_G
   REAL(dp),    DIMENSION(irs:ire,izs:ize)     :: fr_real, gz_real
   REAL(dp),    DIMENSION(irs:ire,izs:ize)     :: fz_real, gr_real, f_times_g
 
   INTEGER :: in, is
   REAL(dp):: kr, kz, kerneln, be_2, bi_2, factn
   REAL(dp):: sigmae2_taue_o2, sigmai2_taui_o2
 
   ! Execution time start
   CALL cpu_time(t0_Sapj)
 
   !!!!!!!!!!!!!!!!!!!! ELECTRON non linear term computation (Sepj)!!!!!!!!!!
   sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the kerneln argument
 
   ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
     jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
       real_data_c = 0._dp ! initialize sum over real nonlinear term
       factn = 1
 
       nloope: DO in = 1,jmaxe+1 ! Loop over laguerre for the sum
 
         krloope: DO ikr = ikrs,ikre ! Loop over kr
           kzloope: DO ikz = ikzs,ikze ! Loop over kz
             kr     = krarray(ikr)
             kz     = kzarray(ikz)
             be_2    = sigmae2_taue_o2 * (kr**2 + kz**2)
             kerneln = be_2**(in-1)/factn * EXP(-be_2)
             ! First convolution terms
             Fr_cmpx(ikr,ikz) = imagu*kr* phi(ikr,ikz) * kerneln
             Fz_cmpx(ikr,ikz) = imagu*kz* phi(ikr,ikz) * kerneln
             ! Second convolution terms
             Gz_cmpx(ikr,ikz) = 0._dp ! initialization of the sum
             Gr_cmpx(ikr,ikz) = 0._dp ! initialization of the sum
             DO is = 1, MIN( in+ij-1, jmaxe+1 ) ! sum truncation on number of moments
               Gz_cmpx(ikr,ikz) = Gz_cmpx(ikr,ikz) + &
                dnjs(in,ij,is) * moments_e(ip,is,ikr,ikz,updatetlevel)
               Gr_cmpx(ikr,ikz) = Gr_cmpx(ikr,ikz) + &
                dnjs(in,ij,is) * moments_e(ip,is,ikr,ikz,updatetlevel)
             ENDDO
             Gz_cmpx(ikr,ikz) = imagu*kz*Gz_cmpx(ikr,ikz)
             Gr_cmpx(ikr,ikz) = imagu*kr*Gr_cmpx(ikr,ikz)
           ENDDO kzloope
         ENDDO krloope
 
         ! First term drphi x dzf
         DO ikr = ikrs, ikre
           DO ikz = ikzs, ikze
-            cmpx_data_f(ikz,ikr-local_nkr_offset) = Fr_cmpx(ikr,ikz)
-            cmpx_data_g(ikz,ikr-local_nkr_offset) = Gz_cmpx(ikr,ikz)
+            cmpx_data_f(ikz,ikr-local_nkr_offset) = Fr_cmpx(ikr,ikz)*AA_r(ikr)*AA_z(ikz)
+            cmpx_data_g(ikz,ikr-local_nkr_offset) = Gz_cmpx(ikr,ikz)*AA_r(ikr)*AA_z(ikz)
           ENDDO
         ENDDO
 
         call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
         call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g)
 
         real_data_c = real_data_c + real_data_f/Nz/Nr  * real_data_g/Nz/Nr
 
         ! Second term -dzphi x drf
         DO ikr = ikrs, ikre
           DO ikz = ikzs, ikze
-            cmpx_data_f(ikz,ikr-local_nkr_offset) = Fz_cmpx(ikr,ikz)
-            cmpx_data_g(ikz,ikr-local_nkr_offset) = Gr_cmpx(ikr,ikz)
+            cmpx_data_f(ikz,ikr-local_nkr_offset) = Fz_cmpx(ikr,ikz)*AA_r(ikr)*AA_z(ikz)
+            cmpx_data_g(ikz,ikr-local_nkr_offset) = Gr_cmpx(ikr,ikz)*AA_r(ikr)*AA_z(ikz)
           ENDDO
         ENDDO
 
         call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
         call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g)
 
         real_data_c = real_data_c - real_data_f/Nz/Nr  * real_data_g/Nz/Nr
 
         IF ( in + 1 .LE. jmaxe+1 ) THEN
           factn = real(in,dp) * factn ! compute (n+1)!
         ENDIF
       ENDDO nloope
 
       ! Put the real nonlinear product into k-space
       call fftw_mpi_execute_dft_r2c(planf, real_data_c, cmpx_data_c)
 
       ! Retrieve convolution in input format
       DO ikr = ikrs, ikre
         DO ikz = ikzs, ikze
           Sepj(ip,ij,ikr,ikz) = cmpx_data_c(ikz,ikr-local_nkr_offset)*AA_r(ikr)*AA_z(ikz)
         ENDDO
       ENDDO
 
     ENDDO jloope
   ENDDO ploope
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   !!!!!!!!!!!!!!!!!!!! ION non linear term computation (Sipj)!!!!!!!!!!
   sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp ! factor of the kerneln argument
 
   ploopi: DO ip = ips_e,ipe_e ! Loop over Hermite moments
     jloopi: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
       real_data_c = 0._dp ! initialize sum over real nonlinear term
       factn = 1
 
       nloopi: DO in = 1,jmaxi+1 ! Loop over laguerre for the sum
 
         krloopi: DO ikr = ikrs,ikre ! Loop over kr
           kzloopi: DO ikz = ikzs,ikze ! Loop over kz
             kr      = krarray(ikr)
             kz      = kzarray(ikz)
             bi_2    = sigmai2_taui_o2 * (kr**2 + kz**2)
             kerneln = bi_2**(in-1)/factn * EXP(-bi_2)
             ! First convolution terms
             Fr_cmpx(ikr,ikz) = imagu*kr* phi(ikr,ikz) * kerneln
             Fz_cmpx(ikr,ikz) = imagu*kz* phi(ikr,ikz) * kerneln
             ! Second convolution terms
             Gz_cmpx(ikr,ikz) = 0._dp ! initialization of the sum
             Gr_cmpx(ikr,ikz) = 0._dp ! initialization of the sum
             DO is = 1, MIN( in+ij-1, jmaxi+1 ) ! sum truncation on number of moments
               Gz_cmpx(ikr,ikz) = Gz_cmpx(ikr,ikz) + &
                dnjs(in,ij,is) * moments_i(ip,is,ikr,ikz,updatetlevel)
               Gr_cmpx(ikr,ikz) = Gr_cmpx(ikr,ikz) + &
                dnjs(in,ij,is) * moments_i(ip,is,ikr,ikz,updatetlevel)
             ENDDO
             Gz_cmpx(ikr,ikz) = imagu*kz*Gz_cmpx(ikr,ikz)
             Gr_cmpx(ikr,ikz) = imagu*kr*Gr_cmpx(ikr,ikz)
           ENDDO kzloopi
         ENDDO krloopi
 
         ! First term drphi x dzf
         DO ikr = ikrs, ikre
           DO ikz = ikzs, ikze
-            cmpx_data_f(ikz,ikr-local_nkr_offset) = Fr_cmpx(ikr,ikz)
-            cmpx_data_g(ikz,ikr-local_nkr_offset) = Gz_cmpx(ikr,ikz)
+            cmpx_data_f(ikz,ikr-local_nkr_offset) = Fr_cmpx(ikr,ikz)*AA_r(ikr)*AA_z(ikz)
+            cmpx_data_g(ikz,ikr-local_nkr_offset) = Gz_cmpx(ikr,ikz)*AA_r(ikr)*AA_z(ikz)
           ENDDO
         ENDDO
 
         call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
         call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g)
 
         real_data_c = real_data_c + real_data_f/Nz/Nr  * real_data_g/Nz/Nr
 
         ! Second term -dzphi x drf
         DO ikr = ikrs, ikre
           DO ikz = ikzs, ikze
-            cmpx_data_f(ikz,ikr-local_nkr_offset) = Fz_cmpx(ikr,ikz)
-            cmpx_data_g(ikz,ikr-local_nkr_offset) = Gr_cmpx(ikr,ikz)
+            cmpx_data_f(ikz,ikr-local_nkr_offset) = Fz_cmpx(ikr,ikz)*AA_r(ikr)*AA_z(ikz)
+            cmpx_data_g(ikz,ikr-local_nkr_offset) = Gr_cmpx(ikr,ikz)*AA_r(ikr)*AA_z(ikz)
           ENDDO
         ENDDO
 
         call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
         call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g)
 
         real_data_c = real_data_c - real_data_f/Nz/Nr  * real_data_g/Nz/Nr
 
         IF ( in + 1 .LE. jmaxe+1 ) THEN
           factn = real(in,dp) * factn ! compute (n+1)!
         ENDIF
       ENDDO nloopi
 
       ! Put the real nonlinear product into k-space
       call fftw_mpi_execute_dft_r2c(planf, real_data_c, cmpx_data_c)
 
       ! Retrieve convolution in input format
       DO ikr = ikrs, ikre
         DO ikz = ikzs, ikze
           Sipj(ip,ij,ikr,ikz) = cmpx_data_c(ikz,ikr-local_nkr_offset)*AA_r(ikr)*AA_z(ikz)
         ENDDO
       ENDDO
 
     ENDDO jloopi
   ENDDO ploopi
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   ! Execution time END
   CALL cpu_time(t1_Sapj)
   tc_Sapj = tc_Sapj + (t1_Sapj - t0_Sapj)
 
 END SUBROUTINE compute_Sapj
diff --git a/src/inital.F90 b/src/inital.F90
index d3e3b8c..3909e64 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -1,242 +1,276 @@
 !******************************************************************************!
 !!!!!! initialize the moments and load/build coeff tables
 !******************************************************************************!
 SUBROUTINE inital
 
   USE basic
   USE model, ONLY : CO, NON_LIN
   USE prec_const
   USE coeff
   USE time_integration
   USE array, ONLY : Sepj,Sipj
+  USE collision
 
   implicit none
 
   real :: start_init, end_init, time_estimation
 
   CALL set_updatetlevel(1)
+
+  CALL evaluate_kernels
+
   !!!!!! Set the moments arrays Nepj, Nipj !!!!!!
   ! WRITE(*,*) 'Init moments'
   IF ( RESTART ) THEN
     CALL load_cp
   ELSE
     CALL init_moments
     !!!!!! Set phi !!!!!!
     IF (my_id .EQ. 0) WRITE(*,*) 'Init phi'
     CALL poisson
   ENDIF
 
 
   !!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!!
   IF ( NON_LIN ) THEN;
     IF (my_id .EQ. 0) WRITE(*,*) 'Init Sapj'
     CALL compute_Sapj
     ! WRITE(*,*) 'Building Dnjs table'
     CALL build_dnjs_table
   ENDIF
+
   !!!!!! Load the full coulomb collision operator coefficients !!!!!!
   IF (CO .EQ. -1) THEN
     IF (my_id .EQ. 0) WRITE(*,*) '=== Load Full Coulomb matrix ==='
     CALL load_FC_mat
     IF (my_id .EQ. 0) WRITE(*,*) '..done'
   ENDIF
 
 END SUBROUTINE inital
 !******************************************************************************!
 
 !******************************************************************************!
 !!!!!!! Initialize the moments randomly
 !******************************************************************************!
 SUBROUTINE init_moments
   USE basic
   USE grid
   USE fields
   USE prec_const
   USE utility, ONLY: checkfield
   USE initial_par
   IMPLICIT NONE
 
   REAL(dp) :: noise
   REAL(dp) :: kr, kz, sigma, gain, kz_shift
   INTEGER, DIMENSION(12) :: iseedarr
 
   ! Seed random number generator
   iseedarr(:)=iseed
   CALL RANDOM_SEED(PUT=iseedarr+my_id)
 
     !**** Broad noise initialization *******************************************
     DO ip=ips_e,ipe_e
       DO ij=ijs_e,ije_e
 
         DO ikr=ikrs,ikre
           DO ikz=ikzs,ikze
             CALL RANDOM_NUMBER(noise)
-            moments_e( ip,ij,     ikr,    ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+            moments_e( ip,ij,     ikr,    ikz, :) = (initback_moments + initnoise_moments*(noise-0.5_dp))*AA_r(ikr)*AA_z(ikz)
           END DO
         END DO
 
         IF ( ikrs .EQ. 1 ) THEN
           DO ikz=2,Nkz/2 !symmetry at kr = 0
             moments_e( ip,ij,1,ikz, :) = moments_e( ip,ij,1,Nkz+2-ikz, :)
           END DO
         ENDIF
 
       END DO
     END DO
 
     DO ip=ips_i,ipe_i
       DO ij=ijs_i,ije_i
 
         DO ikr=ikrs,ikre
           DO ikz=ikzs,ikze
             CALL RANDOM_NUMBER(noise)
-            moments_i( ip,ij,ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+            moments_i( ip,ij,ikr,ikz, :) = (initback_moments + initnoise_moments*(noise-0.5_dp))*AA_r(ikr)*AA_z(ikz)
           END DO
         END DO
 
         IF ( ikrs .EQ. 1 ) THEN
           DO ikz=2,Nkz/2 !symmetry at kr = 0
             moments_i( ip,ij,1,ikz, :) = moments_i( ip,ij,1,Nkz+2-ikz, :)
           END DO
         ENDIF
 
       END DO
     END DO
 
   ! ENDIF
 END SUBROUTINE init_moments
 !******************************************************************************!
 
 !******************************************************************************!
 !!!!!!! Load moments from a previous save
 !******************************************************************************!
 SUBROUTINE load_cp
   USE basic
   USE futils,          ONLY: openf, closef, getarr, getatt, isgroup, isdataset
   USE grid
   USE fields
   USE diagnostics_par
   USE time_integration
   IMPLICIT NONE
 
   INTEGER :: rank, sz_, n_
   INTEGER ::  dims(1) = (/0/)
   CHARACTER(LEN=50) :: dset_name
 
   WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(rstfile0),'_',job2load,'.h5'
 
   IF (my_id .EQ. 0) WRITE(*,'(3x,a)') "Resume from previous run"
 
   CALL openf(rstfile, fidrst,mpicomm=MPI_COMM_WORLD)
 
     n_ = 0
     WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_
     DO WHILE (isdataset(fidrst, dset_name))
       n_ = n_ + 1
       WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_
     ENDDO
     n_ = n_ - 1
     WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_
     WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_i", n_
 
     ! Read state of system from restart file
     WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_i", n_
     CALL getarr(fidrst, dset_name, moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),pardim=3)
     WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_
     CALL getarr(fidrst, dset_name, moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),pardim=3)
     WRITE(dset_name, "(A, '/', i6.6)") "/Basic/phi", n_
     CALL getarr(fidrst, dset_name, phi(ikrs:ikre,ikzs:ikze),pardim=1)
 
     ! Read time dependent attributes
     CALL getatt(fidrst, dset_name, 'cstep', cstep)
     CALL getatt(fidrst, dset_name, 'time', time)
     CALL getatt(fidrst, dset_name, 'jobnum', jobnum)
     jobnum = jobnum+1
     CALL getatt(fidrst, dset_name, 'iframe2d',iframe2d)
     CALL getatt(fidrst, dset_name, 'iframe5d',iframe5d)
     iframe2d = iframe2d-1; iframe5d = iframe5d-1
 
   CALL closef(fidrst)
 
   IF (my_id .EQ. 0) WRITE(*,'(3x,a)') "Reading from restart file "//TRIM(rstfile)//" completed!"
 
 END SUBROUTINE load_cp
 !******************************************************************************!
 
 !******************************************************************************!
-!!!!!!! Build the Laguerre-Laguerre coupling coefficient table
+!!!!!!! Build the Laguerre-Laguerre coupling coefficient table for nonlin
 !******************************************************************************!
 SUBROUTINE build_dnjs_table
   USE basic
   USE array, Only : dnjs
   USE grid, Only : jmaxe, jmaxi
   USE coeff
   IMPLICIT NONE
 
   INTEGER :: in, ij, is, J
   INTEGER :: n_, j_, s_
 
   J = max(jmaxe,jmaxi)
 
   DO in = 1,J+1 ! Nested dependent loops to make benefit from dnjs symmetry
     n_ = in - 1
     DO ij = in,J+1
       j_ = ij - 1
       DO is = ij,J+1
         s_ = is - 1
 
         dnjs(in,ij,is) = TO_DP(ALL2L(n_,j_,s_,0))
         ! By symmetry
         dnjs(in,is,ij) = dnjs(in,ij,is)
         dnjs(ij,in,is) = dnjs(in,ij,is)
         dnjs(ij,is,in) = dnjs(in,ij,is)
         dnjs(is,ij,in) = dnjs(in,ij,is)
         dnjs(is,in,ij) = dnjs(in,ij,is)
       ENDDO
     ENDDO
   ENDDO
-END SUBROUTINE
+END SUBROUTINE build_dnjs_table
 !******************************************************************************!
 
 !******************************************************************************!
-!!!!!!! Load the Full coulomb coefficient table from COSOlver results
+!!!!!!! Evaluate the kernels once for all
 !******************************************************************************!
-SUBROUTINE load_FC_mat ! Works only for a partiular file for now with P,J <= 15,6
+SUBROUTINE evaluate_kernels
   USE basic
+  USE array, Only : kernel_e, kernel_i
   USE grid
-  USE initial_par
-  USE array
-  use model
-  USE futils, ONLY : openf, getarr, closef
+  use model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, DK
   IMPLICIT NONE
 
-  INTEGER :: ip2,ij2
-  INTEGER :: fid1, fid2, fid3, fid4
-
-  !!!!!!!!!!!! Electron matrices !!!!!!!!!!!!
-  ! get the self electron colision matrix
-  CALL openf(selfmat_file,fid1, 'r', 'D',mpicomm=MPI_COMM_WORLD);
-  CALL getarr(fid1, '/Caapj/Ceepj', Ceepj) ! get array (moli format)
-  CALL closef(fid1)
-  ! get the Test and Back field electron ion collision matrix
-  CALL openf(eimat_file,fid2, 'r', 'D');
-  CALL getarr(fid2, '/Ceipj/CeipjT', CeipjT)
-  CALL getarr(fid2, '/Ceipj/CeipjF', CeipjF)
-  CALL closef(fid2)
-
-  !!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!!
-  ! get the self electron colision matrix
-  CALL openf(selfmat_file, fid3, 'r', 'D',mpicomm=MPI_COMM_WORLD);
-  IF ( (pmaxe .EQ. pmaxi) .AND. (jmaxe .EQ. jmaxi) ) THEN ! if same degrees ion and electron matrices are alike so load Ceepj
-    CALL getarr(fid3, '/Caapj/Ceepj', Ciipj) ! get array (moli format)
-  ELSE
-    CALL getarr(fid3, '/Caapj/Ciipj', Ciipj) ! get array (moli format)
-  ENDIF
-  CALL closef(fid3)
-  ! get the Test and Back field ion electron collision matrix
-  CALL openf(iemat_file, fid4, 'r', 'D');
-  CALL getarr(fid4, '/Ciepj/CiepjT', CiepjT)
-  CALL getarr(fid4, '/Ciepj/CiepjF', CiepjF)
-  CALL closef(fid4)
-END SUBROUTINE load_FC_mat
+  REAL(dp)    :: factj, j_dp, j_int
+  REAL(dp)    :: sigmae2_taue_o2, sigmai2_taui_o2
+  REAL(dp)    :: be_2, bi_2, alphaD
+  REAL(dp)    :: kr, kz, kperp2
+
+  !!!!! Electron kernels !!!!!
+  !Precompute species dependant factors
+  sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument
+
+  factj = 1.0 ! Start of the recursive factorial
+  DO ij = ijs_e, ije_e
+    j_int = jarray_e(ij)
+    j_dp = REAL(j_int,dp) ! REAL of degree
+
+    ! Recursive factorial
+    IF (j_dp .GT. 0) THEN
+      factj = factj * j_dp
+    ELSE
+      factj = 1._dp
+    ENDIF
+
+    DO ikr = ikrs,ikre
+      kr     = krarray(ikr)
+      DO ikz = ikzs,ikze
+        kz    = kzarray(ikz)
+        be_2  =  (kr**2 + kz**2) * sigmae2_taue_o2
+
+        kernel_e(ij, ikr, ikz) = be_2**(ij)/factj * exp(-be_2)
+
+      ENDDO
+    ENDDO
+  ENDDO
+
+  !!!!! Ion kernels !!!!!
+  sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp ! (b_a/2)^2 = (kperp sqrt(2 tau_a) sigma_a/2)^2
+  factj = 1.0 ! Start of the recursive factorial
+  DO ij = ijs_i, ije_i
+    j_int = jarray_e(ij)
+    j_dp = REAL(j_int,dp) ! REAL of degree
+
+    ! Recursive factorial
+    IF (j_dp .GT. 0) THEN
+      factj = factj * j_dp
+    ELSE
+      factj = 1._dp
+    ENDIF
+
+    DO ikr = ikrs,ikre
+      kr     = krarray(ikr)
+      DO ikz = ikzs,ikze
+        kz    = kzarray(ikz)
+        bi_2  =  (kr**2 + kz**2) * sigmai2_taui_o2
+
+        kernel_i(ij, ikr, ikz) = bi_2**(ij)/factj * exp(-bi_2)
+
+      ENDDO
+    ENDDO
+  ENDDO
+
+END SUBROUTINE evaluate_kernels
 !******************************************************************************!