diff --git a/.gitignore b/.gitignore
index 0b32f91f21a82f38540d5e890279fb4e5a9c0278..41e06c0521f0d41a7b64b76e7ea7edce32cd114f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -31,3 +31,4 @@ wk/fort.90
 .directory
 checkpoint/
 FM/
+*.out
diff --git a/src/inital.F90 b/src/inital.F90
index 5ee70cf33be802eeb8ba7f6898108001192455ee..916e6c878f67c593eb1aed39bd5159ca8863b36d 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -127,9 +127,7 @@ SUBROUTINE load_cp
   WRITE(*,'(3x,a)') "Resume from previous run"
 
   CALL openf(rstfile, fidrst,mpicomm=MPI_COMM_WORLD)
-
-  ! IF (isgroup(fidrst,'/Basic/moments_e')) THEN
-    ! Getting the last saved checkpoint
+  
     n_ = 0
     WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_
     DO WHILE (isdataset(fidrst, dset_name))
@@ -153,20 +151,6 @@ SUBROUTINE load_cp
     CALL getatt(fidrst, dset_name, 'iframe5d',iframe5d)
     iframe2d = iframe2d-1; iframe5d = iframe5d-1
 
-  ! ELSE
-  !   CALL getatt(fidrst, '/Basic', 'cstep', cstep)
-  !   CALL getatt(fidrst, '/Basic', 'time', time)
-  !   CALL getatt(fidrst, '/Basic', 'jobnum', jobnum)
-  !   jobnum = jobnum+1
-  !   CALL getatt(fidrst, '/Basic', 'iframe2d',iframe2d)
-  !   CALL getatt(fidrst, '/Basic', 'iframe5d',iframe5d)
-  !   iframe2d = iframe2d-1; iframe5d = iframe5d-1
-  !
-  !   ! Read state of system from restart file
-  !   CALL getarr(fidrst, '/Basic/moments_e', moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),pardim=3)
-  !   CALL getarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),pardim=3)
-  ! ENDIF
-
   CALL closef(fidrst)
 
   WRITE(*,'(3x,a)') "Reading from restart file "//TRIM(rstfile)//" completed!"
@@ -226,7 +210,7 @@ SUBROUTINE load_FC_mat ! Works only for a partiular file for now with P,J <= 15,
 
   !!!!!!!!!!!! Electron matrices !!!!!!!!!!!!
   ! get the self electron colision matrix
-  CALL openf(selfmat_file,fid1, 'r', 'D');
+  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
@@ -237,7 +221,7 @@ SUBROUTINE load_FC_mat ! Works only for a partiular file for now with P,J <= 15,
 
   !!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!!
   ! get the self electron colision matrix
-  CALL openf(selfmat_file, fid3, 'r', 'D');
+  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
diff --git a/src/srcinfo.h b/src/srcinfo.h
index d9f601e55ad0d1c4261bf844a62604ac06f30ee8..a6564e943784033d91b96ac8cf1eecfa2ddc3fa6 100644
--- a/src/srcinfo.h
+++ b/src/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
 character(len=20) AUTHOR
 character(len=40) EXECDATE
 character(len=40) HOST
-parameter (VERSION='003c315-dirty')
+parameter (VERSION='2c9ad20-dirty')
 parameter (BRANCH='MPI')
 parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Fri Nov 20 14:14:05 CET 2020')
+parameter (EXECDATE='Wed Nov 25 17:26:17 CET 2020')
 parameter (HOST ='spcpc606')
diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h
index d9f601e55ad0d1c4261bf844a62604ac06f30ee8..a6564e943784033d91b96ac8cf1eecfa2ddc3fa6 100644
--- a/src/srcinfo/srcinfo.h
+++ b/src/srcinfo/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
 character(len=20) AUTHOR
 character(len=40) EXECDATE
 character(len=40) HOST
-parameter (VERSION='003c315-dirty')
+parameter (VERSION='2c9ad20-dirty')
 parameter (BRANCH='MPI')
 parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Fri Nov 20 14:14:05 CET 2020')
+parameter (EXECDATE='Wed Nov 25 17:26:17 CET 2020')
 parameter (HOST ='spcpc606')
diff --git a/wk/flux_results.m b/wk/flux_results.m
index d1c37d7fe1421c1473535a6ecdc5aafe2b4f4712..2ff3e2232d0bbe9da3cffaef1452d5d476f5be65 100644
--- a/wk/flux_results.m
+++ b/wk/flux_results.m
@@ -1,5 +1,5 @@
 default_plots_options
-if 0
+if 1
 %% Compute time average and std of the mean flow
 t0 = 180; t1 = 400; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D)); 
 range  = it0:it1;
@@ -8,6 +8,7 @@ stdev  = std(Flux_ri(range))^(.5)
 figure
 hist(Flux_ri(range),20); xlabel('$\Gamma$')
 end
+if 0
 %% Handwritten results for nu = 0.01
 % High definition results (256x128)
 Results_256x128.Gamma = [0.02, 0.03, 0.20, 0.037,  2.7, 2.25,    4, 5e-4, 2e-3, 0.03];
@@ -89,4 +90,5 @@ grid on; title('$\nu = 0.1$')
 legend('Mix. Length, Ricci 2006','$P=2$, $J=1$','$P=3$, $J=2$','$P=4$, $J=2$')
 FIGNAME = [SIMDIR,'flux_study_nu_1e-2.png'];
 saveas(fig,FIGNAME);
-disp(['Figure saved @ : ',FIGNAME])
\ No newline at end of file
+disp(['Figure saved @ : ',FIGNAME])
+end
\ No newline at end of file
diff --git a/wk/fort.90 b/wk/fort.90
index c56c511af13a5bff5fde9c1b0a3b8c2a52b232f8..34f66b4e388173bb8292de9b681f57cb0feaed56 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,42 +1,42 @@
 &BASIC
   nrun   = 100000000
   dt     = 0.01
-  tmax   = 400
+  tmax   = 100
   RESTART = .false.
 /
 &GRID
-  pmaxe =3
-  jmaxe = 2
-  pmaxi = 3
-  jmaxi = 2
-  Nr   = 256
-  Lr = 66
-  Nz   = 256
-  Lz = 66
+  pmaxe =2
+  jmaxe = 1
+  pmaxi = 2
+  jmaxi = 1
+  Nr   = 400
+  Lr = 100
+  Nz   = 400
+  Lz = 100
   kpar = 0
   CANCEL_ODD_P = .false.
 /
 &OUTPUT_PAR
-  nsave_0d = 50
+  nsave_0d = 100
   nsave_1d = 0
-  nsave_2d = 50
-  nsave_5d = 200
+  nsave_2d = 100
+  nsave_5d = 1000
   write_Na00    = .true.
   write_moments = .true.
   write_phi     = .true.
   write_non_lin     = .true.
   write_doubleprecision = .true.
-  resfile0      = '../results/test_parallel/256x128_L_66_Pe_3_Je_2_Pi_3_Ji_2_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-06/outputs'
-  rstfile0      = '../results/test_parallel/256x128_L_66_Pe_3_Je_2_Pi_3_Ji_2_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-06/checkpoint'
+  resfile0      = '../results/test_parallel/400x200_L_100_Pe_2_Je_1_Pi_2_Ji_1_nB_0.6_nN_1_nu_1e-01_FC_mu_5e-04/outputs'
+  rstfile0      = '../results/test_parallel/400x200_L_100_Pe_2_Je_1_Pi_2_Ji_1_nB_0.6_nN_1_nu_1e-01_FC_mu_5e-04/checkpoint'
   job2load      = 0
 /
 &MODEL_PAR
   ! Collisionality
-  CO      = -2
+  CO      = -1
   DK      = .false.
   NON_LIN = .true.
-  mu      = 5e-06
-  nu      = 0.01
+  mu      = 0.0005
+  nu      = 0.1
   tau_e   = 1
   tau_i   = 1
   sigma_e = 0.023338
@@ -45,7 +45,7 @@
   q_i     = 1
   eta_n   = 1
   eta_T   = 0
-  eta_B   = 0.5
+  eta_B   = 0.6
   lambdaD = 0
   kr0KH   = 0
   A0KH    = 0
@@ -55,9 +55,9 @@
   initback_moments  =0
   initnoise_moments =1e-05
   iseed             =42
-  selfmat_file      ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_3_Jmaxe_2_Pmaxi_3_Jmaxi_2_pamaxx_10.h5'
-  eimat_file        ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_3_Jmaxe_2_Pmaxi_3_Jmaxi_2_pamaxx_10_tau_1.0000_mu_0.0233.h5'
-  iemat_file        ='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_3_Jmaxe_2_Pmaxi_3_Jmaxi_2_pamaxx_10_tau_1.0000_mu_0.0233.h5'
+  selfmat_file      ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10.h5'
+  eimat_file        ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5'
+  iemat_file        ='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5'
 /
 &TIME_INTEGRATION_PAR
   numerical_scheme='RK4'
diff --git a/wk/parameters_ZP.m b/wk/parameters_ZP.m
index 5025785f199d7622f3d96476429bd0bc7ad4ce52..e41c7578b3ff70d2f81b23d06fe0903f60208947 100644
--- a/wk/parameters_ZP.m
+++ b/wk/parameters_ZP.m
@@ -14,19 +14,18 @@ NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
 N       = 256;     % Frequency gridpoints (Nkr = N/2)
 L       = 66;     % Size of the squared frequency domain
-PMAXE   = 5;     % Highest electron Hermite polynomial degree
-JMAXE   = 3;     % Highest ''       Laguerre ''
-PMAXI   = 5;     % Highest ion      Hermite polynomial degree
-JMAXI   = 3;     % Highest ''       Laguerre ''
-CANCEL_ODD_P = 0;% Cancels the odd polynomials degree
+PMAXE   = 2;     % Highest electron Hermite polynomial degree
+JMAXE   = 1;     % Highest ''       Laguerre ''
+PMAXI   = 2;     % Highest ion      Hermite polynomial degree
+JMAXI   = 1;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS 
 TMAX    = 400;  % Maximal time unit
-DT      = 4e-3;   % Time step
-SPS0D   = 1.0;      % Sampling per time unit for 2D arrays
+DT      = 1e-2;   % Time step
+SPS0D   = 1;      % 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
 RESTART = 1;      % To restart from last checkpoint
-JOB2LOAD= 0;
+JOB2LOAD= 1;
 %% OPTIONS
 SIMID   = 'ZP';  % Name of the simulation
 CO      = -1;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
@@ -40,5 +39,6 @@ KREQ0   = 0;      % put kr = 0
 KPAR    = 0.0;    % Parellel wave vector component
 LAMBDAD = 0.0; 
 NON_LIN = 1 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
+CANCEL_ODD_P = 0;% Cancels the odd polynomials degree
 
 setup
diff --git a/wk/test_parallel.m b/wk/test_parallel.m
index cf70fd464c7a286883d3829f9f1a2a7fd33a5990..923583ae26a0949489c316ab7dd25abf91d903d1 100644
--- a/wk/test_parallel.m
+++ b/wk/test_parallel.m
@@ -4,31 +4,31 @@ addpath(genpath('../matlab')) % ... add
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1e-2;   % Collision frequency
+NU      = 1e-1;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-ETAB    = 0.5;    % Magnetic gradient
+ETAB    = 0.6;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
-MU      = 5e-6;   % Hyper diffusivity coefficient
+MU      = 5e-4;   % Hyper diffusivity coefficient
 NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
-N       = 256;     % Frequency gridpoints (Nkr = N/2)
-L       = 66;     % Size of the squared frequency domain
-PMAXE   = 3;     % Highest electron Hermite polynomial degree
-JMAXE   = 2;     % Highest ''       Laguerre ''
-PMAXI   = 3;     % Highest ion      Hermite polynomial degree
-JMAXI   = 2;     % Highest ''       Laguerre ''
+N       = 400;     % Frequency gridpoints (Nkr = N/2)
+L       = 100;     % Size of the squared frequency domain
+PMAXE   = 2;     % Highest electron Hermite polynomial degree
+JMAXE   = 1;     % Highest ''       Laguerre ''
+PMAXI   = 2;     % Highest ion      Hermite polynomial degree
+JMAXI   = 1;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS 
-TMAX    = 400;  % Maximal time unit
+TMAX    = 100;  % Maximal time unit
 DT      = 1e-2;   % Time step
-SPS0D   = 2;    % Sampling per time unit for profiler
-SPS2D   = 2;      % Sampling per time unit for 2D arrays
-SPS5D   = 0.5;    % Sampling per time unit for 5D arrays
+SPS0D   = 1;    % Sampling per time unit for profiler
+SPS2D   = 1;      % Sampling per time unit for 2D arrays
+SPS5D   = 0.1;    % Sampling per time unit for 5D arrays
 RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 0;
 %% OPTIONS
 SIMID   = 'test_parallel';  % Name of the simulation
-CO      = -2;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+CO      = -1;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% unused