diff --git a/Makefile b/Makefile
index 0e18cad629a8ff74aff5a45eaf969cd55038198f..456be934ecabdd52c3034b6992e5629457b36c76 100644
--- a/Makefile
+++ b/Makefile
@@ -42,7 +42,7 @@ cleanbin:
$(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 \
@@ -68,6 +68,9 @@ $(OBJDIR)/utility_mod.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 $@
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..ff27b697a2ab8c926a13c71f82f0bb068014ed2b
--- /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 a6cefadbf6e3e51dd64361ddb53986d5738af87f..7d1dfb6ab64666b7cc8f99fa066408790b4cb8bb 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -62,8 +62,8 @@ SUBROUTINE compute_Sapj
! 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
@@ -75,8 +75,8 @@ SUBROUTINE compute_Sapj
! 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
@@ -140,8 +140,8 @@ SUBROUTINE compute_Sapj
! 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
@@ -153,8 +153,8 @@ SUBROUTINE compute_Sapj
! 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
diff --git a/src/inital.F90 b/src/inital.F90
index d3e3b8c827c42f335fe1837f18237ba6035085b6..3909e64cbaae590424a318f5044fcb77dce6a7dc 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -9,12 +9,16 @@ SUBROUTINE inital
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
@@ -34,6 +38,7 @@ SUBROUTINE inital
! 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 ==='
@@ -71,7 +76,7 @@ SUBROUTINE init_moments
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
@@ -90,7 +95,7 @@ SUBROUTINE init_moments
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
@@ -164,7 +169,7 @@ 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
@@ -195,48 +200,77 @@ SUBROUTINE build_dnjs_table
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
!******************************************************************************!