diff --git a/matlab/create_mov.m b/matlab/create_mov.m
index 66da1aac2742e96ea68587c0cb89529a526b9d94..93bc14cbc81fac00d4beae3f159ad9c1e2ce725d 100644
--- a/matlab/create_mov.m
+++ b/matlab/create_mov.m
@@ -16,7 +16,7 @@ figure('Color','white','Position', [100, 100, 400, 400]);
     pcolor(X,Y,FIELD(:,:,1)); % to set up
     colormap gray
     axis tight manual % this ensures that getframe() returns a consistent size
-    if INTERP
+    if 1
         shading interp;
     end
     in      = 1;
@@ -24,7 +24,7 @@ figure('Color','white','Position', [100, 100, 400, 400]);
     for n = FRAMES % loop over selected frames
         scale = max(max(abs(FIELD(:,:,n))));
         pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot
-        if INTERP
+        if 1
             shading interp; 
         end
         set(pclr, 'edgecolor','none'); axis square;
diff --git a/matlab/load_marconi.m b/matlab/load_marconi.m
index fa714cbb18ae2f2ac69654c847f1d3dc9143bcf3..05c03182ee796a02a5d90dc6b7e263a8aaf351c4 100644
--- a/matlab/load_marconi.m
+++ b/matlab/load_marconi.m
@@ -3,15 +3,23 @@ function [ RESDIR ] = load_marconi( outfilename )
 %   Detailed explanation goes here
     hostfolder  = outfilename(1:end-8);
     hostfile    = [hostfolder,'/out*'];
-    MISCDIR     = ['/misc/HeLaZ_outputs/results/',outfilename(46:end-8),'/'];
+    MISCDIR     = ['/misc/HeLaZ_outputs/',outfilename(46:end-8),'/'];
     RESDIR      = ['../',outfilename(46:end-8),'/'];
     miscfolder =  [MISCDIR,'.'];
-    localfolder = [RESDIR,'.'];
+    system(['mkdir -p ',miscfolder]);
+    disp(['mkdir -p ',miscfolder]);
+    resultfolder = [RESDIR,'.'];
+    % SCP the output file from marconi to misc folder of SPCPC
     CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',miscfolder];
     disp(CMD);
     system(CMD);
+    % Load the fort.90 as well in misc folder
     CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfolder,'/fort.90',' ',miscfolder];
     disp(CMD);
     system(CMD);
+    % Put it also in the result directory
+    CMD = ['scp -r ahoffman@login.marconi.cineca.it:',hostfolder,'/fort.90',' ',resultfolder];
+    disp(CMD);
+    system(CMD);
 end
 
diff --git a/matlab/photomaton.m b/matlab/photomaton.m
index a73cc5f8c244879ec8d7c78631bcb9a9a3c33d17..6f4343921e1a6ee8b1feb4b4bea2994322674d58 100644
--- a/matlab/photomaton.m
+++ b/matlab/photomaton.m
@@ -1,14 +1,13 @@
 if 0
 %% Photomaton : real space
-% FIELD = ni00; FNAME = 'ni'; XX = RR; YY = ZZ;
-% FIELD = phi; FNAME = 'phi'; XX = RR; YY = ZZ;
-FIELD = dr2phi; FNAME = 'dr2phi'; XX = RR; YY = ZZ;
-% FIELD = fftshift(abs(Ni00),2); FNAME = 'Fni'; XX = fftshift(KR,2); YY = fftshift(KZ,2);
-% FIELD = fftshift(abs(PHI),2);  FNAME = 'Fphi'; XX = fftshift(KR,2); YY = fftshift(KZ,2);
-tf = 150;  [~,it1] = min(abs(Ts2D-tf));
-tf = 300;  [~,it2] = min(abs(Ts2D-tf)); 
-tf = 800; [~,it3] = min(abs(Ts2D-tf));
-tf = 1200; [~,it4] = min(abs(Ts2D-tf));
+% FIELD = ni00; FNAME = 'ni'; FIELDLTX = '$N_i^{00}$'; XX = RR; YY = ZZ;
+% FIELD = ne00; FNAME = 'ne'; FIELDLTX = '$N_e^{00}$'; XX = RR; YY = ZZ;
+FIELD = phi; FNAME = 'phi'; FIELDLTX = '$\phi$'; XX = RR; YY = ZZ;
+% FIELD = dr2phi; FNAME = 'dr2phi'; XX = RR; YY = ZZ;
+tf = 800;  [~,it1] = min(abs(Ts2D-tf));
+tf = 900;  [~,it2] = min(abs(Ts2D-tf)); 
+tf = 1000; [~,it3] = min(abs(Ts2D-tf));
+tf = 1100; [~,it4] = min(abs(Ts2D-tf));
 fig = figure; FIGNAME = [FNAME,'_snaps','_',PARAMS]; set(gcf, 'Position',  [100, 100, 1500, 350])
 plt = @(x) x;%./max(max(x));
     subplot(141)
@@ -35,30 +34,74 @@ plt = @(x) x;%./max(max(x));
         colormap gray
         xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); 
         title(sprintf('$t c_s/R=%.0f$',Ts2D(it4)));
-suptitle(['$\',FNAME,'$, $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
-    ', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']);
+    suptitle([FIELDLTX,', $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
+        ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+        ' $\mu_{hd}=$',num2str(MU)]);
 save_figure
 end
 
+if 0
+%% Photomaton : real space
+% FIELD = Ni00; FNAME = 'fni00'; FIELDLTX = '$\tilde N_i^{00}$'; XX = RR; YY = ZZ;
+% FIELD = Ne00; FNAME = 'fne00'; FIELDLTX = '$\tilde N_e^{00}$'; XX = RR; YY = ZZ;
+FIELD = PHI; FNAME = 'fphi'; FIELDLTX = '$\tilde \phi$'; XX = RR; YY = ZZ;
+tf = 700;  [~,it1] = min(abs(Ts2D-tf));
+tf = 800;  [~,it2] = min(abs(Ts2D-tf)); 
+tf = 900; [~,it3] = min(abs(Ts2D-tf));
+tf = 1500; [~,it4] = min(abs(Ts2D-tf));
+XLIM = [0,1.0]; YLIM = [-.5;.5];
+fig = figure; FIGNAME = [FNAME,'_snaps','_',PARAMS]; set(gcf, 'Position',  [100, 100, 1500, 400])
+plt = @(x) fftshift((abs(x)),2)/max(max(abs(x)));
+    subplot(141)
+        DATA = plt(FIELD(:,:,it1));
+        pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        colormap gray
+        xlabel('$k_r \rho_s$'); ylabel('$k_z \rho_s$'); xlim(XLIM); ylim(YLIM); 
+        title(sprintf('$t c_s/R=%.0f$',Ts2D(it1)));
+    subplot(142)
+        DATA = plt(FIELD(:,:,it2));
+        pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        colormap gray
+        xlabel('$k_r \rho_s$'); ylabel('$k_z \rho_s$'); set(gca,'ytick',[]);  xlim(XLIM); ylim(YLIM); 
+        title(sprintf('$t c_s/R=%.0f$',Ts2D(it2)));
+    subplot(143)
+        DATA = plt(FIELD(:,:,it3));
+        pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        colormap gray
+        xlabel('$k_r \rho_s$'); ylabel('$k_z \rho_s$');set(gca,'ytick',[]);  xlim(XLIM); ylim(YLIM); 
+        title(sprintf('$t c_s/R=%.0f$',Ts2D(it3)));
+    subplot(144)
+        DATA = plt(FIELD(:,:,it4));
+        pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        colormap gray
+        xlabel('$k_r \rho_s$'); ylabel('$k_z \rho_s$'); set(gca,'ytick',[]);  xlim(XLIM); ylim(YLIM); 
+        title(sprintf('$t c_s/R=%.0f$',Ts2D(it4)));
+    suptitle([FIELDLTX,', $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
+        ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+        ' $\mu_{hd}=$',num2str(MU)]);
+save_figure
+end
+
+
 %%
 if 0
 %% Show frame in kspace
-tf = 1000; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
+tf = 800; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
 fig = figure; FIGNAME = ['krkz_',sprintf('t=%.0f',Ts2D(it2)),'_',PARAMS];set(gcf, 'Position',  [100, 100, 700, 600])
-CLIM = [0,1e3]
-    subplot(223); plt = @(x) fftshift((abs(x)),2);
+CLIM = [0,1];
+    subplot(223); plt = @(x) fftshift((abs(x)),2)/max(max(abs(x)));
         pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it2))); set(pclr, 'edgecolor','none');
         caxis(CLIM); colormap hot
         xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); legend('$|\hat\phi|$');
-    subplot(222); plt = @(x) fftshift(abs(x),2);
+    subplot(222);
         pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it2))); set(pclr, 'edgecolor','none');
         caxis(CLIM); colormap hot
         xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{00}|$');
-    subplot(221); plt = @(x) fftshift(abs(x),2);
+    subplot(221);
         pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ne00(:,:,it2))); set(pclr, 'edgecolor','none');
         caxis(CLIM); colormap hot
         xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_e^{00}|$');
-    subplot(224); plt = @(x) fftshift((abs(x)),2);
+    subplot(224);
     colorbar;
         caxis(CLIM); colormap hot
 %         pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar;
diff --git a/matlab/setup.m b/matlab/setup.m
index d8c53ddff95de4e29764017a527da53f503f0fcc..a6b17e7b5f76e012ce3a806de627b2141c2c500b 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -38,18 +38,21 @@ MODEL.lambdaD = LAMBDAD;
 % if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KR0KH),'_A_',num2str(A0KH)]; end;
 % Time integration and intialization parameters
 TIME_INTEGRATION.numerical_scheme  = '''RK4''';
-if INIT_PHI; INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end;
-INITIAL.init_background  = 0.0e-5;
+if (INIT_PHI && INIT_ZF == 0); INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end;
+INITIAL.INIT_ZF = INIT_ZF;
+INITIAL.init_background  = (INIT_ZF>0)*ZF_AMP;
 INITIAL.init_noiselvl = NOISE0;
 INITIAL.iseed             = 42;
 INITIAL.selfmat_file = '''null''';
 INITIAL.eimat_file = '''null''';
 INITIAL.iemat_file = '''null''';
-if (CO == -3) % Write matrice filename for Full Coulomb DK
+INITIAL.mat_file   = '''null''';
+if (CO == -3) % Write matrice filename for Full Coulomb
     cmat_pmaxe = 25;
     cmat_jmaxe = 18;
     cmat_pmaxi = 25;
     cmat_jmaxi = 18;
+    INITIAL.mat_file = ['''../../../iCa/FC_P_25_J_18_N_200_dk_0.05236_MFLR_0.h5'''];
     INITIAL.selfmat_file = ...
         ['''../../../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
@@ -62,11 +65,12 @@ if (CO == -3) % Write matrice filename for Full Coulomb DK
         ['''../../../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
         num2str(cmat_jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
-elseif (CO == -2) % Write matrice filename for DK Sugama
+elseif (CO == -2) % Write matrice filename for Sugama
     cmat_pmaxe = 10;
     cmat_jmaxe = 5;
     cmat_pmaxi = 10;
     cmat_jmaxi = 5;
+    INITIAL.mat_file = ['''../../../iCa/SG_P_10_J_5_N_200_dk_0.05236_MFLR_0.h5'''];
     INITIAL.selfmat_file = ...
         ['''../../../iCa/self_Coll_GKE_0_GKI_0_ESELF_3_ISELF_3_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
@@ -84,6 +88,7 @@ elseif (CO == 2) % Write matrice filename for Sugama GK
     cmat_jmaxe = 5;
     cmat_pmaxi = 10;
     cmat_jmaxi = 5;
+    INITIAL.mat_file = ['''../../../iCa/SG_P_10_J_5_N_200_dk_0.05236_MFLR_0.h5'''];
     INITIAL.selfmat_file = ...
         ['''../../../iCa/self_Coll_GKE_1_GKI_1_ESELF_3_ISELF_3_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 88e496c867f310a2af767384a2fd74bdeb76f26f..ed9b721d09c438e4dc578d75f634ba614beac0dc 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -64,13 +64,14 @@ fprintf(fid,'/\n');
 
 fprintf(fid,'&INITIAL_CON\n');
 fprintf(fid,['  INIT_NOISY_PHI    =', INITIAL.init_noisy_phi,'\n']);
-fprintf(fid,['  INIT_ZF_PHI       =', INITIAL.INIT_ZF_PHI,'\n']);
+fprintf(fid,['  INIT_ZF       =', num2str(INITIAL.INIT_ZF),'\n']);
 fprintf(fid,['  init_background  =', num2str(INITIAL.init_background),'\n']);
 fprintf(fid,['  init_noiselvl =', num2str(INITIAL.init_noiselvl),'\n']);
 fprintf(fid,['  iseed             =', num2str(INITIAL.iseed),'\n']);
 fprintf(fid,['  selfmat_file      =', INITIAL.selfmat_file,'\n']);
 fprintf(fid,['  eimat_file        =', INITIAL.eimat_file,'\n']);
 fprintf(fid,['  iemat_file        =', INITIAL.iemat_file,'\n']);
+fprintf(fid,['  mat_file        =', INITIAL.mat_file,'\n']);
 fprintf(fid,'/\n');
 
 fprintf(fid,'&TIME_INTEGRATION_PAR\n');
diff --git a/matlab/write_sbash_marconi.m b/matlab/write_sbash_marconi.m
index e544d84c83af6f67b157847de6d3e3dcc15bffdf..8caf105d24c4fc6b9a529f47ae167e7b8812e1da 100644
--- a/matlab/write_sbash_marconi.m
+++ b/matlab/write_sbash_marconi.m
@@ -36,7 +36,7 @@ fprintf(fid,[...
 '#SBATCH --account=FUA35_TSVVT421\n',...
 '#SBATCH --partition=skl_fua_',CLUSTER.PART,'\n',...
 'module load autoload hdf5 fftw\n',...
-'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KR)]);
+'srun --cpu-bind=cores ./../../../bin/',EXECNAME,' ',num2str(NP_P),' ',num2str(NP_KR)]);
 
 fclose(fid);
 system(['cp batch_script.sh ',BASIC.RESDIR,'/.']);
diff --git a/wk/ZF_fourier_analysis.m b/wk/ZF_fourier_analysis.m
index c63ea03c1e6fa89cc6f72c0df99d7fa4b5b9d03b..3610ba525141bf10329705abd631a62e9d5adcfa 100644
--- a/wk/ZF_fourier_analysis.m
+++ b/wk/ZF_fourier_analysis.m
@@ -1,42 +1,50 @@
 %% Zonal flow spectral analysis
 fig = figure; FIGNAME = ['zonal_flow_spectral_analysis_',PARAMS];
-tend = Ts0D(end); tstart = tend - TAVG;
+tend = Ts0D(end); tstart = tend-TAVG ;
 [~,its0D] = min(abs(Ts0D-tstart));
 [~,ite0D]   = min(abs(Ts0D-tend));
+[~,its2D] = min(abs(Ts2D-tstart));
+[~,ite2D]   = min(abs(Ts2D-tend));
 set(gcf, 'Position',  [100, 100, 800, 400])
     % Time series analysis (burst period and time frequencies spectrum)
     subplot(121)
     samplerate = Ts0D(2)-Ts0D(1);
     Y = log(PGAMMA_RI(its0D:ite0D)*(2*pi/Nr/Nz)^2);
     [n,~] = size(Y);
-    Yy= fft(y);
+    Yy= fft(Y);
     Pot = Yy .* conj(Yy) / n;
     Freq = (samplerate / n * (1:n))';
     Pot(1) = 0;
-    nmax = min(200,round(n/2));
+    nmax = min(20,round(n/2));
     [amax, itmax] = max(Pot);
     plot((0:nmax-1) , Pot(1:nmax)/amax,'DisplayName','$\Gamma_r(\omega)$');hold on;
-    plot([itmax-1,itmax-1],[0,1],'--k', 'DisplayName',['$T_{per}\approx',num2str(round(2*pi/Freq(itmax))),'\rho_s/c_s$']);
+    plot([itmax-1,itmax-1],[0,1],'--k', 'DisplayName',['$T_{per}\approx',num2str(round(1/Freq(itmax))),'L_\perp/c_s$']);
     legend('show'); grid on; box on; xlabel('Period number'); yticks([]);
-    title('ZF temporal spectrum')
+    title('$\Gamma_r$ temporal spectrum')
     % Space analysis (spatial period of ZF)
     subplot(122)
-    Yy = mean(mean(1i*KR.*PHI(:,:,its2D:ite2D),3),2); hold on;
-    [n,~] = size(Yy);
-    Pot = Yy .* conj(Yy) / n;
-    [amax, ikmax] = max(Pot);
-    nmax = 20;
-    plot(0:nmax,Pot(1:nmax+1)/amax,'DisplayName','$\tilde V_E(k_r)$');
+    nmax = 20; n = numel(r);
+    [TT,NN] = meshgrid(Ts2D(its2D:ite2D),0:n-1);
+    Pot = NN;
+    for it = 1:ite2D-its2D+1
+        Y = mean(real(drphi(:,:,it)),2);
+        Yy = fft(Y);
+        [n,~] = size(Yy);
+        Pot(:,it) = Yy .* conj(Yy) / n;
+    end
+    [amax, ikmax] = max(mean(Pot,2));
+%     pclr = pcolor(NN(1:nmax,:),TT(1:nmax,:),Pot(1:nmax,:)); set(pclr, 'edgecolor','none'); hold on;
+    plot(0:nmax,mean(Pot(1:nmax+1,:),2)/amax,'DisplayName','$\langle\partial_r\phi\rangle_z (k_r)$'); hold on;
     plot([ikmax-1,ikmax-1],[0,1],'--k', 'DisplayName',['$L_z=',num2str(2*pi/kr(ikmax)),'\rho_s$']);
     grid on; box on;
     title('ZF spatial spectrum')
-    xlabel('mode number');  yticks([]); legend('show')
+    xlabel('radial mode number');  yticks([]); legend('show')
 save_figure
 
 %% Shear and phi amp phase space
 fig = figure; FIGNAME = ['phi_shear_phase_space_',PARAMS];
 set(gcf, 'Position',  [100, 100, 700, 500])
-t1 = Ts2D(end); t0 = t1 - TAVG;
+t1 = Ts2D(end); t0 = 0;
 [~,its2D] = min(abs(Ts2D-t0)); [~,ite2D] = min(abs(Ts2D-t1));
 scatter(phi_maxr_avgz(its2D:ite2D),shear_maxr_avgz(its2D:ite2D),35,Ts2D(its2D:ite2D),'.',...
     'DisplayName',PARAMS); cbar = colorbar;ylabel(cbar,'$t c_s/\rho_s$','Interpreter','LaTeX')
@@ -48,6 +56,22 @@ grid on; title('ES pot. vs Shear phase space')
 % plot(phi_avgr_avgz(its2D:ite2D),shear_avgr_avgz(its2D:ite2D),'-')
 save_figure
 
+if 0
+%% density and phi phase space
+fig = figure; FIGNAME = ['phi_ni_phase_space_',PARAMS];
+set(gcf, 'Position',  [100, 100, 700, 500])
+t1 = Ts2D(end); t0 = 0;
+[~,its2D] = min(abs(Ts2D-t0)); [~,ite2D] = min(abs(Ts2D-t1));
+scatter3(max(mean(ni00(:,:,its2D:ite2D),2),[],1),phi_maxr_avgz(its2D:ite2D),shear_maxr_avgz(its2D:ite2D),35,Ts2D(its2D:ite2D),'.',...
+    'DisplayName',PARAMS); cbar = colorbar;ylabel(cbar,'$t c_s/\rho_s$','Interpreter','LaTeX')
+hold on
+xlabel('$\langle n_i^{00} \rangle_z^r$'); ylabel('$\langle \phi \rangle_z^r$'); zlabel('$\langle dV_E/dr \rangle_z^r$')
+grid on; title('ES pot. vs Shear phase space')
+% plot(phi_avgr_maxz(its2D:ite2D),shear_avgr_maxz(its2D:ite2D),'-')
+% plot(phi_maxr_maxz(its2D:ite2D),shear_maxr_maxz(its2D:ite2D),'-')
+% plot(phi_avgr_avgz(its2D:ite2D),shear_avgr_avgz(its2D:ite2D),'-')
+% save_figure
+end
 %% Non zonal quantities
 PHI_NZ = PHI;
 PHI_NZ(ikmax-1:ikmax+1,:,:) = 0;
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index e66ebd2a12931ef09ed498f381e7c503e78c9eb3..29245fc1903201ea874a6a3eecf0b06564befd86 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -1,20 +1,21 @@
 addpath(genpath('../matlab')) % ... add
 %% Load results
 outfile ='';
-if 0 % Local results
-outfile ='';
+if 1 % Local results
 outfile ='';
 outfile ='';
 outfile ='';
+outfile ='v2.6_P_2_J_1/200x100_L_120_P_2_J_1_eta_0.4_nu_1e-01_DGGK_CLOS_0_mu_1e-01';
+
     BASIC.RESDIR      = ['../results/',outfile,'/'];
     BASIC.MISCDIR     = ['../results/',outfile,'/'];
 end
-if 1 % Marconi results
+if 0 % Marconi results
 outfile ='';
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.9_nu_5e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt';
 % load_marconi(outfile);
 BASIC.RESDIR      = ['../',outfile(46:end-8),'/'];
 BASIC.MISCDIR     = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
@@ -182,7 +183,8 @@ fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM),'_',PARAMS];
 set(gcf, 'Position',  [100, 100, 900, 800])
 subplot(111); 
     suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
-        ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$']);
+        ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+        ' $\mu_{hd}=$',num2str(MU)]);
     subplot(421); 
     for ip = 1:Npe
         for ij = 1:Nje
@@ -215,25 +217,38 @@ subplot(111);
         legend(['Gyro. flux';'Part. flux']);
         grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma_{r,i}$')
 %         ylim([0,2.0]);
+    if(1)
     subplot(223)
         plot(kz,g_(1,:),'-','DisplayName','Linear growth rate'); hold on;
         plot([max(kz)*2/3,max(kz)*2/3],[0,10],'--k', 'DisplayName','2/3 Orszag AA');
         grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma R/c_s$'); legend('show');
         ylim([0,max(g_(1,:))]); xlim([0,max(kz)]);
-    subplot(224)
+        shearplot = 426; phiplot = 428;
+    else
+    shearplot = 223; phiplot = 224;      
+    end
+    subplot(shearplot)
+        clr      = line_colors(min(ip,numel(line_colors(:,1))),:);
+        lstyle   = line_styles(min(ij,numel(line_styles)));
+        plot(Ts2D,shear_maxr_maxz,'DisplayName','$\max_{r,z}(s)$'); hold on;
+        plot(Ts2D,shear_maxr_avgz,'DisplayName','$\max_{r}\langle s \rangle_z$'); hold on;
+        plot(Ts2D,shear_avgr_maxz,'DisplayName','$\max_{z}\langle s \rangle_r$'); hold on;
+        plot(Ts2D,shear_avgr_avgz,'DisplayName','$\langle s \rangle_{r,z}$'); hold on;
+    grid on; xlabel('$t c_s/R$'); ylabel('$shear$'); 
+    subplot(phiplot)
         clr      = line_colors(min(ip,numel(line_colors(:,1))),:);
         lstyle   = line_styles(min(ij,numel(line_styles)));
         plot(Ts2D,phi_maxr_maxz,'DisplayName','$\max_{r,z}(\phi)$'); hold on;
         plot(Ts2D,phi_maxr_avgz,'DisplayName','$\max_{r}\langle\phi\rangle_z$'); hold on;
         plot(Ts2D,phi_avgr_maxz,'DisplayName','$\max_{z}\langle\phi\rangle_r$'); hold on;
         plot(Ts2D,phi_avgr_avgz,'DisplayName','$\langle\phi\rangle_{r,z}$'); hold on;
-    grid on; xlabel('$t c_s/R$'); ylabel('$T_e/e$'); legend('show');
+    grid on; xlabel('$t c_s/R$'); ylabel('$E.S. pot$');
 save_figure
 end
 
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
-TAVG = 1500; % Averaging time duration
+TAVG = 3500; % Averaging time duration
 %Compute steady radial transport
 tend = Ts0D(end); tstart = tend - TAVG;
 [~,its0D] = min(abs(Ts0D-tstart));
@@ -283,7 +298,7 @@ fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position',
 save_figure
 end
 
-if 0
+if 1
 %% Space time diagramms
 tstart = 0; tend = Ts2D(end);
 [~,itstart] = min(abs(Ts2D-tstart));
@@ -295,24 +310,41 @@ fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position',  [100, 10
         pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
         shading interp
         colormap hot;
-        caxis([0.0,0.02]);
+        caxis([0.0,max(max(mean(ni00(:,:,its2D:ite2D).*dzphi(:,:,its2D:ite2D),2)))]);
          xticks([]); ylabel('$r/\rho_s$')
-        legend('$\Gamma_r(r)=\langle n_i^{GC}\partial_z\phi\rangle_z$')
+        legend('Radial part. transport $\langle n_i^{00}\partial_z\phi\rangle_z$')
         title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
         ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
         ' $\mu_{hd}=$',num2str(MU)]);
     subplot(212)
         pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
+        fieldmax = max(max(mean(abs(drphi(:,:,its2D:ite2D)),2)));
+        caxis([-fieldmax,fieldmax]);
         xlabel('$t c_s/R$'), ylabel('$r/\rho_s$')
-        legend('$\langle \partial_r\phi\rangle_z$')        
+        legend('Zonal flow $\langle \partial_r\phi\rangle_z$')        
 save_figure
 end
 
+if 0
+%% Averaged phi, ZF and shear profiles
+figure;
+plt = @(x) squeeze(mean(mean(real(x(:,:,its2D:ite2D)),2),3))/max(abs(squeeze(mean(mean(real(x(:,:,its2D:ite2D)),2),3))));
+subplot(311)
+plot(r,plt(phi)); hold on;
+xlim([-60,60]); xticks([]); yticks([]);
+subplot(312)
+plot(r,plt(drphi)); hold on;
+xlim([-60,60]);xticks([]);yticks([]);
+subplot(313)
+plot(r,plt(dr2phi)); hold on;
+xlim([-60,60]); xlabel('$r/\rho_s$');xticks([]);yticks([]);
+end
+
 %%
-t0    = 0000;
+t0    = 500;
 [~, it02D] = min(abs(Ts2D-t0));
 [~, it05D] = min(abs(Ts5D-t0));
-skip_ = 10; 
+skip_ = 5; 
 DELAY = 1e-3*skip_;
 FRAMES_2D = it02D:skip_:numel(Ts2D);
 FRAMES_5D = it05D:skip_:numel(Ts5D);
@@ -362,10 +394,26 @@ FIELDNAME = '$\langle\phi\rangle_{z}$'; XNAME = '$r/\rho_s$';
 create_gif_1D_phi
 end
 if 0
+%% flow averaged on z
+GIFNAME = ['zf_z0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING=1;
+FIELD =(squeeze(mean(real(drphi),2))); linestyle = '-.k'; FRAMES = FRAMES_2D;
+X = (r); T = Ts2D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r);
+FIELDNAME = '$\langle\phi\rangle_{z}$'; XNAME = '$r/\rho_s$';
+create_gif_1D_phi
+end
+if 0
+%% shear averaged on z
+GIFNAME = ['shear_z0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING=1;
+FIELD =(squeeze(mean(real(dr2phi),2))); linestyle = '-.k'; FRAMES = FRAMES_2D;
+X = (r); T = Ts2D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r);
+FIELDNAME = '$\langle\phi\rangle_{z}$'; XNAME = '$r/\rho_s$';
+create_gif_1D_phi
+end
+if 0
 %% phi @ kz = 0
 GIFNAME = ['phi_kz0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING = 0;
-FIELD =squeeze(log10(abs(PHI(no_AA,1,:)))); linestyle = '-.'; FRAMES = FRAMES_2D;
-X = kr(no_AA); T = Ts2D; YMIN = -30; YMAX = 6; XMIN = min(kr); XMAX = max(kr);
+FIELD =squeeze(log10(abs(PHI(:,1,:)))); linestyle = '-.'; FRAMES = FRAMES_2D;
+X = kr; T = Ts2D; YMIN = -0; YMAX = 6; XMIN = min(kr); XMAX = max(kr);
 FIELDNAME = '$\log_{10}|\tilde\phi(k_z=0)|$'; XNAME = '$k_r\rho_s$';
 create_gif_1D
 end
diff --git a/wk/compile_cosolver_mat.m b/wk/compile_cosolver_mat.m
new file mode 100644
index 0000000000000000000000000000000000000000..5a3970d26613508c74ae589f45b93cbd97178b64
--- /dev/null
+++ b/wk/compile_cosolver_mat.m
@@ -0,0 +1,141 @@
+addpath(genpath('../matlab')) % ... add
+%% Grid configuration
+N       = 200;     % Frequency gridpoints (Nkr = N/2)
+L       = 120;     % Size of the squared frequency domain
+dk      = 2.*pi/L;
+kmax    = N/2*dk;
+kr      = dk*(0:N/2);
+kz      = dk*(0:N/2);
+[KZ, KR]= meshgrid(kz,kr);
+KPERP   = sqrt(KR.^2 + KZ.^2);
+kperp   = reshape(KPERP,[1,numel(kr)^2]);
+kperp   = uniquetol(kperp,1e-14);
+Nperp   = numel(kperp);
+%% Model
+    % ! 0 = nothing 
+    % ! 1 = coulomb
+    % ! 2 = pitch-angle (with constant coll.freq.)
+    % ! 3 = sugama
+    % ! 4 = pitch-angle with veloty dependent coll. freq.
+    % ! 5 = improved sugama
+    % ! 6 = Hirschman-Sigmar Clarke
+    % ! 7 = Abel (only for like species)
+    % ! 8 = conserving momentun pitch-angle operator (with veloty dependent coll. freq.)
+    % ! 9 = GK coulomb polarization term
+CO = 3;
+GK = 1;
+P = 10;
+J = 5;
+M_FLR= 0; %to increase the NFLR
+if 0
+    %% plot the kperp distribution
+   figure
+   plot(kperp)
+end
+% %% Check if the differences btw kperp is larger than naming precision
+%%
+%% We compute only on a kperp grid with dk space from 0 to kperpmax
+kperpmax = sqrt(2) * kmax;
+kperp = unique([0:dk:kperpmax,kperpmax]);
+%% Naming
+if CO == 1; CONAME = 'FC'; end;
+if CO == 2; CONAME = 'PA'; end;
+if CO == 3; CONAME = 'SG'; end;
+matfilename = ['../iCa/',CONAME,'_P_',num2str(P),'_J_',num2str(J),...
+    '_N_',num2str(N),'_dk_',num2str(dk),'_MFLR_',num2str(M_FLR),'.h5'];
+
+n_ = 1;
+for k_ = kperp
+    disp(['-Writing matrix for kperp = ',num2str(k_)])
+    %% Script to run COSOlver in order to create needed collision matrices
+    COSOlver_path = '../../Documents/MoliSolver/COSOlver/';
+    COSOLVER.pmaxe = P;
+    COSOLVER.jmaxe = J;
+    COSOLVER.pmaxi = P;
+    COSOLVER.jmaxi = J;
+    COSOLVER.kperp = k_;
+
+    COSOLVER.neFLR    = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2)))+M_FLR; % rule of thumb for sum truncation
+    COSOLVER.niFLR    = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2)))+M_FLR;
+    COSOLVER.idxT4max = 40;
+
+    COSOLVER.neFLRs = 0; %  ... only for GK abel 
+    COSOLVER.npeFLR = 0; %  ... only for GK abel 
+    COSOLVER.niFLRs = 0; %  ... only for GK abel 
+    COSOLVER.npiFLR = 0; %  ... only for GK abel 
+
+    COSOLVER.gk = GK; 
+    COSOLVER.co = CO;
+    
+    k_string      = sprintf('%0.4f',k_);
+    self_mat_file_name = ['../iCa/self_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
+        '_ESELF_',num2str(COSOLVER.co),'_ISELF_',num2str(COSOLVER.co),...
+        '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
+        '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
+        '_JE_12',...
+        '_NFLR_',num2str(COSOLVER.neFLR),...
+        '_kperp_',k_string,'.h5'];
+    ei_mat_file_name = ['../iCa/ei_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
+        '_ETEST_',num2str(COSOLVER.co),'_EBACK_',num2str(COSOLVER.co),...
+        '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
+        '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
+        '_JE_12_tau_1.0000_mu_0.0233',...
+        '_NFLRe_',num2str(COSOLVER.neFLR),'_NFLRi_',num2str(COSOLVER.neFLR)...
+        '_kperp_',k_string,'.h5'];
+    ie_mat_file_name = ['../iCa/ie_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
+        '_ITEST_',num2str(COSOLVER.co),'_IBACK_',num2str(COSOLVER.co),...
+        '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
+        '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
+        '_JE_12_tau_1.0000_mu_0.0233',...
+        '_NFLRe_',num2str(COSOLVER.neFLR),'_NFLRi_',num2str(COSOLVER.neFLR)...
+        '_kperp_',k_string,'.h5'];
+   
+    %% Load self matrix
+    C_self = h5read(self_mat_file_name,'/Caapj/Ceepj');
+    sz_ = size(C_self);
+    % Write it in the compiled file
+    h5create(matfilename,['/Caapj/',k_string],sz_)
+    h5write(matfilename,['/Caapj/',k_string],C_self)
+    
+    %% Load ei matrices
+    % Field
+    C_eiF = h5read(ei_mat_file_name,'/Ceipj/CeipjF');
+    sz_ = size(C_eiF);
+    h5create(matfilename,['/CeipjF/',k_string],sz_)
+    h5write(matfilename,['/CeipjF/',k_string],C_eiF)
+    % Test
+    C_eiT = h5read(ei_mat_file_name,'/Ceipj/CeipjT');
+    sz_ = size(C_eiT);
+    h5create(matfilename,['/CeipjT/',k_string],sz_)
+    h5write(matfilename,['/CeipjT/',k_string],C_eiT)
+    
+    %% Load ie matrices
+    % Field
+    C_ieF = h5read(ie_mat_file_name,'/Ciepj/CiepjF');
+    sz_ = size(C_ieF);
+    h5create(matfilename,['/CiepjF/',k_string],sz_)
+    h5write(matfilename,['/CiepjF/',k_string],C_ieF)
+    % Test
+    C_ieT = h5read(ie_mat_file_name,'/Ciepj/CiepjT');
+    sz_ = size(C_eiT);
+    h5create(matfilename,['/CiepjT/',k_string],sz_)
+    h5write(matfilename,['/CiepjT/',k_string],C_ieT)
+    
+    %% Copy fort.90 input file and put grid params
+    if(k_ == 0)
+        h5create(matfilename,'/dk',1);
+        h5write (matfilename,'/dk',dk);   
+        h5create(matfilename,'/N',1);
+        h5write (matfilename,'/N',N);
+        h5create(matfilename,'/Pmaxe',1);
+        h5write (matfilename,'/Pmaxe',P);   
+        h5create(matfilename,'/Jmaxe',1);
+        h5write (matfilename,'/Jmaxe',J);   
+        h5create(matfilename,'/Pmaxi',1);
+        h5write (matfilename,'/Pmaxi',P);  
+        h5create(matfilename,'/Jmaxi',1);
+        h5write (matfilename,'/Jmaxi',J);   
+    end
+    
+end
+disp(['File saved @ :',matfilename])
\ No newline at end of file
diff --git a/wk/compute_collision_mat.m b/wk/compute_collision_mat.m
index 6a1f5e2842c7b67821fa6a4e4649d0e1468c5f05..a429236ee9d6dd6d6816c4e33a446120b7778bd4 100644
--- a/wk/compute_collision_mat.m
+++ b/wk/compute_collision_mat.m
@@ -11,22 +11,28 @@ KPERP   = sqrt(KR.^2 + KZ.^2);
 kperp   = reshape(KPERP,[1,numel(kr)^2]);
 kperp   = uniquetol(kperp,1e-14);
 Nperp   = numel(kperp);
-if 0
-    %% plot the kperp distribution
-   figure
-   plot(kperp)
-end
-% %% Check if the differences btw kperp is larger than naming precision
-% dkperp  = diff(kperp);
-% warning = sum(dkperp<0.0002);
-% if warning > 0
-%     disp('Warning : dkperp < 0.0002');
-% end
+%% Model
+    % ! 0 = nothing 
+    % ! 1 = coulomb
+    % ! 2 = pitch-angle (with constant coll.freq.)
+    % ! 3 = sugama
+    % ! 4 = pitch-angle with veloty dependent coll. freq.
+    % ! 5 = improved sugama
+    % ! 6 = Hirschman-Sigmar Clarke
+    % ! 7 = Abel (only for like species)
+    % ! 8 = conserving momentun pitch-angle operator (with veloty dependent coll. freq.)
+    % ! 9 = GK coulomb polarization term
+CO = 3;
+GK = 1;
+
 %%
 %% We compute only on a kperp grid with dk space from 0 to kperpmax
 kperpmax = sqrt(2) * kmax;
 kperp = unique([0:dk:kperpmax,kperpmax]);
 %%
+%% If DK collision, compute only kperp = 0
+if (GK == 0); kperp = 0; end;
+%% Kperp scan
 n_ = 1;
 for k_ = kperp
     disp(['----------Computing matrix ',num2str(n_),'/',num2str(numel(kperp)),'----------'])
@@ -38,8 +44,8 @@ for k_ = kperp
     COSOLVER.jmaxi = 05;
     COSOLVER.kperp = k_;
 
-    COSOLVER.neFLR    = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2)))+2; % rule of thumb for sum truncation
-    COSOLVER.niFLR    = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2)))+2;
+    COSOLVER.neFLR    = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2))); % rule of thumb for sum truncation
+    COSOLVER.niFLR    = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2)));
     COSOLVER.idxT4max = 40;
 
     COSOLVER.neFLRs = 0; %  ... only for GK abel 
@@ -47,18 +53,8 @@ for k_ = kperp
     COSOLVER.niFLRs = 0; %  ... only for GK abel 
     COSOLVER.npiFLR = 0; %  ... only for GK abel 
 
-    COSOLVER.gk = 1; 
-    COSOLVER.co = 3;
-    % ! 0 = nothing 
-    % ! 1 = coulomb
-    % ! 2 = pitch-angle (with constant coll.freq.)
-    % ! 3 = sugama
-    % ! 4 = pitch-angle with veloty dependent coll. freq.
-    % ! 5 = improved sugama
-    % ! 6 = Hirschman-Sigmar Clarke
-    % ! 7 = Abel (only for like species)
-    % ! 8 = conserving momentun pitch-angle operator (with veloty dependent coll. freq.)
-    % ! 9 = GK coulomb polarization term
+    COSOLVER.gk = GK; 
+    COSOLVER.co = CO;
 
     write_fort90_COSOlver
     
diff --git a/wk/continue_multiple_runs_marconi.m b/wk/continue_multiple_runs_marconi.m
index 32202e1723f08691133a749766a3467fee9b2196..827b4713f69d41c233ad741df2c61f6f36140e62 100644
--- a/wk/continue_multiple_runs_marconi.m
+++ b/wk/continue_multiple_runs_marconi.m
@@ -1,18 +1,11 @@
 %% Paste the list of continue_run calls
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.5_nu_1e+00_SGGK_CLOS_0_mu_2e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e+00_SGGK_CLOS_0_mu_2e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.7_nu_1e+00_SGGK_CLOS_0_mu_2e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.8_nu_1e+00_SGGK_CLOS_0_mu_2e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_5e-01_DGGK_CLOS_0_mu_2e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.7_nu_5e-01_DGGK_CLOS_0_mu_2e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.8_nu_5e-01_DGGK_CLOS_0_mu_2e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.9_nu_5e-01_DGGK_CLOS_0_mu_2e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.9_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.5_nu_1e-02_SGGK_CLOS_0_mu_2e-02/out.txt')
 
 %% Functions to modify preexisting fort.90 input file and launch on marconi
 function [] = continue_run(outfilename)
+    EXECNAME = 'helaz_2.62';
     %% CLUSTER PARAMETERS
-    CLUSTER.PART  = 'prod';     % dbg or prod
+    CLUSTER.PART  = 'dbg';     % dbg or prod
     CLUSTER.TIME  = '24:00:00'; % allocation time hh:mm:ss
     if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME  = '00:30:00'; end;
     CLUSTER.MEM   = '64GB';     % Memory
@@ -51,7 +44,7 @@ function [] = continue_run(outfilename)
         line = A{33};
         line = line(end-2:end);
         if(line(1) == '='); line = line(end); end;
-        J2L = str2num(line) + 1;
+        J2L = str2num(line)+1;
     end
     % Change job 2 load in fort.90
     A{33} = ['  job2load      = ',num2str(J2L)];
diff --git a/wk/linear_study.m b/wk/linear_study.m
index 7b86263349924b29407c961dd2ad76987062828a..9bf350cf332112ca36160a6f11c980c0dd10f0e8 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -1,6 +1,6 @@
 for NU = [1.0]
 for ETAB = [0.5]
-for CO = [1,2]
+for CO = [-2 2]
 %clear all;
 addpath(genpath('../matlab')) % ... add
 default_plots_options
@@ -14,18 +14,18 @@ TAU     = 1.0;    % e/i temperature ratio
 % ETAB    = 0.5;
 ETAN    = 1.0;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
-NU_HYP  = 1.0;   % Hyperdiffusivity coefficient
+NU_HYP  = 0.0;   % Hyperdiffusivity coefficient
 LAMBDAD = 0.0;
 NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
-N       = 200;     % Frequency gridpoints (Nkr = N/2)
+N       = 20;     % Frequency gridpoints (Nkr = N/2)
 L       = 120;     % Size of the squared frequency domain
 KREQ0   = 1;      % put kr = 0
 MU_P    = 0.0;     % Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
 MU_J    = 0.0;     % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
 %% TIME PARMETERS
-TMAX    = 200;  % Maximal time unit
-DT      = 2e-2;   % Time step
+TMAX    = 100;  % Maximal time unit
+DT      = 1e-2;   % Time step
 SPS0D   = 0.5;      % Sampling per time unit for 2D arrays
 SPS2D   = 1;      % Sampling per time unit for 2D arrays
 SPS5D   = 1/2;    % Sampling per time unit for 5D arrays
@@ -33,7 +33,8 @@ SPSCP   = 0;    % Sampling per time unit for checkpoints
 RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 00;
 %% OPTIONS
-SIMID   = 'linear_study';  % Name of the simulation
+% SIMID   = 'test';  % Name of the simulation
+SIMID   = 'v2.6_lin_analysis';  % Name of the simulation
 NON_LIN = 0 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
 % Collision operator
 % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
@@ -41,7 +42,8 @@ NON_LIN = 0 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
 CLOS    = 0;   % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
 NL_CLOS = 0;   % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS)
 KERN    = 0;   % Kernel model (0 : GK)
-INIT_PHI= 0;   % Start simulation with a noisy phi and moments
+INIT_PHI= 0;   % Start simulation with a noisy phi
+INIT_ZF = 0;   % Start simulation with a noisy phi
 %% OUTPUTS
 W_DOUBLE = 0;
 W_GAMMA  = 1;
@@ -62,11 +64,12 @@ MU      = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
 %% PARAMETER SCANS
 if 1
 %% Parameter scan over PJ
-% PA = [2, 3, 4, 6, 8, 10];
-% JA = [1, 2, 2, 3, 4,  5];
+% PA = [2, 3, 4, 6];%, 8, 10];
+% JA = [1, 2, 2, 3];%, 4,  5];
 PA = [4];
-JA = [2];
-DTA= DT./sqrt(JA);
+JA = [4];
+DTA= DT./sqrt(JA)/4;
+% DTA= DT;
 mup_ = MU_P;
 muj_ = MU_J;
 Nparam = numel(PA);
@@ -80,35 +83,36 @@ for i = 1:Nparam
     PMAXE = PA(i); PMAXI = PA(i);
     JMAXE = JA(i); JMAXI = JA(i);
     DT = DTA(i);
-    MU_P = mup_/PMAXE^2;
-    MU_J = muj_/JMAXE^3;
+%     MU_P = mup_/PMAXE^2;
+%     MU_J = muj_/JMAXE^3;
     setup
     % Run linear simulation
 %     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ./../../../bin/helaz 1 1; cd ../../../wk'])
-    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6; cd ../../../wk'])
-%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 2 3; cd ../../../wk'])
-    % Load and process results
+%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6; cd ../../../wk'])
+    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 2 3; cd ../../../wk'])
+%     Load and process results
     %%
+    filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5'];
     load_results
     tend   = Ts2D(end); tstart   = 0.4*tend;
     [~,itstart] = min(abs(Ts2D-tstart));
     [~,itend]   = min(abs(Ts2D-tend));
     for ikr = 1:N/2+1
-        gamma_Ni00(i,ikr) = LinearFit_s(Ts2D(itstart:itend),squeeze(abs(Ni00(ikr,1,itstart:itend))));
+        gamma_Ni00(i,ikr) = (LinearFit_s(Ts2D(itstart:itend)',(squeeze(abs(Ni00(ikr,1,itstart:itend))))));
         Ni00_ST(i,ikr,1:numel(Ts2D)) = squeeze((Ni00(ikr,1,:)));
     end
     tend   = Ts5D(end); tstart   = 0.4*tend;
     [~,itstart] = min(abs(Ts5D-tstart));
     [~,itend]   = min(abs(Ts5D-tend));
     for ikr = 1:N/2+1
-        gamma_Nipj(i,ikr) = LinearFit_s(Ts5D(itstart:itend),squeeze(max(max(abs(Nipj(:,:,ikr,1,itstart:itend)),[],1),[],2)));
+        gamma_Nipj(i,ikr) = LinearFit_s(Ts5D(itstart:itend)',squeeze(max(max(abs(Nipj(:,:,ikr,1,itstart:itend)),[],1),[],2)));
     end
     gamma_Ni00(i,:) = real(gamma_Ni00(i,:) .* (gamma_Ni00(i,:)>=0.0));
     gamma_Nipj(i,:) = real(gamma_Nipj(i,:) .* (gamma_Nipj(i,:)>=0.0));
 %     kzmax = abs(kr(ikzmax));
 %     Bohm_transport(i) = ETAB/ETAN*gmax/kzmax^2;
     % Clean output
-    system(['rm -r ',BASIC.RESDIR])
+    system(['rm -r ',BASIC.RESDIR]);
 end
 
 if 1
@@ -116,7 +120,7 @@ if 1
 SCALE = 1;%sqrt(2);
 fig = figure; FIGNAME = 'linear_study';
 plt = @(x) x;
-subplot(211)
+% subplot(211)
     for i = 1:Nparam
         clr       = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:);
         linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1);
@@ -126,22 +130,22 @@ subplot(211)
             'DisplayName',['$P=$',num2str(PA(i)),', $J=$',num2str(JA(i))]);
         hold on;
     end
-    grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{00})\rho_s/c_s$'); xlim([0.0,max(kr)]);
-    title(['$\eta_B=',num2str(ETAB),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, ', CLOSNAME])
-    legend('show')
-subplot(212)
-    for i = 1:Nparam
-        clr       = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:);
-        linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1);
-        plot(plt(SCALE*kr),plt(gamma_Nipj(i,:)),...
-            'Color',clr,...
-            'LineStyle',linestyle{1},...
-            'DisplayName',['$P=$',num2str(PA(i)),', $J=$',num2str(JA(i))]);
-        hold on;
-    end
-    grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(\max_{pj}N_i^{pj})\rho_s/c_s$'); xlim([0.0,max(kr)]);
+    grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{00})L_\perp/c_s$'); xlim([0.0,max(kr)]);
     title(['$\eta_B=',num2str(ETAB),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, ', CLOSNAME])
     legend('show')
+% subplot(212)
+%     for i = 1:Nparam
+%         clr       = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:);
+%         linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1);
+%         plot(plt(SCALE*kr),plt(gamma_Nipj(i,:)),...
+%             'Color',clr,...
+%             'LineStyle',linestyle{1},...
+%             'DisplayName',['$P=$',num2str(PA(i)),', $J=$',num2str(JA(i))]);
+%         hold on;
+%     end
+%     grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(\max_{pj}N_i^{pj})\rho_s/c_s$'); xlim([0.0,max(kr)]);
+%     title(['$\eta_B=',num2str(ETAB),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, ', CLOSNAME])
+%     legend('show')
 saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
 saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']);
 end
diff --git a/wk/load_multiple_outputs_marconi.m b/wk/load_multiple_outputs_marconi.m
index 955bbee4e4b9b4ff4da5913c3d49b45c01f86784..9863a0234ebdfb8c6e9d4b2c602eac5e70fbabdb 100644
--- a/wk/load_multiple_outputs_marconi.m
+++ b/wk/load_multiple_outputs_marconi.m
@@ -1 +1,32 @@
+addpath(genpath('../matlab')) % ... add
 %% Paste the list of simulation results to load
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.5_nu_1e-02_SGDK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGDK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.7_nu_1e-02_SGDK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.8_nu_1e-02_SGDK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.9_nu_1e-02_SGDK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.5_nu_1e-03_SGDK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-03_SGDK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.7_nu_1e-03_SGDK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.8_nu_1e-03_SGDK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.9_nu_1e-03_SGDK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.5_nu_1e-02_DGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.7_nu_1e-02_DGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.8_nu_1e-02_DGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.9_nu_1e-02_DGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.5_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.7_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.8_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.9_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.5_nu_1e-02_SGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.7_nu_1e-02_SGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.8_nu_1e-02_SGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.9_nu_1e-02_SGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.5_nu_1e-03_SGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-03_SGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.7_nu_1e-03_SGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.8_nu_1e-03_SGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.9_nu_1e-03_SGGK_CLOS_0_mu_2e-02/out.txt')
diff --git a/wk/local_run.m b/wk/local_run.m
index b8dd706583ae8ff824ee9106691bed332fcdb9bb..f3de8f68e06c90c6149e6b96fffcbf2d98719cbc 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -4,19 +4,19 @@ addpath(genpath('../matlab')) % ... add
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.1;   % Collision frequency
-ETAB    = 0.6;    % Magnetic gradient
+NU      = 1.0e-1;   % Collision frequency
+ETAB    = 0.4;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
-NU_HYP  = 0.2;
+NU_HYP  = 1.0;
 %% GRID PARAMETERS
-N       = 200;     % Frequency gridpoints (Nkr = N/2)
+N       = 20;     % Frequency gridpoints (Nkr = N/2)
 L       = 120;     % Size of the squared frequency domain
-P       = 2;
+P       = 0;
 J       = 1;
-MU_P    = 0.0/P^2;     % Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J    = 0.0/J^2;     % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
+MU_P    = 0.0;     % Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J    = 0.0;     % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
 %% TIME PARAMETERS
-TMAX    = 500;  % Maximal time unit
+TMAX    = 1000;  % Maximal time unit
 DT      = 1e-2;   % Time step
 SPS0D   = 1;    % Sampling per time unit for profiler
 SPS2D   = 1;      % Sampling per time unit for 2D arrays
@@ -27,13 +27,13 @@ JOB2LOAD= 0;
 %% OPTIONS AND NAMING
 % Collision operator
 % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
-CO      = 2;
+CO      = 1;
 CLOS    = 0;   % Closure model (0: =0 truncation)
 NL_CLOS = -1;   % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-% SIMID   = 'local_low_pj_analysis';  % Name of the simulation
-% SIMID   = 'Sugama_no_ei_ie_coll';  % Name of the simulation
-SIMID   = ['v2.5_P_',num2str(P),'_J_',num2str(J)];  % Name of the simulation
+SIMID   = 'test';  % Name of the simulation
+% SIMID   = ['v2.6_P_',num2str(P),'_J_',num2str(J)];  % Name of the simulation
 NON_LIN = 1;   % activate non-linearity (is cancelled if KREQ0 = 1)
+INIT_ZF = 0; ZF_AMP = 0.0;
 %% OUTPUTS
 W_DOUBLE = 0;
 W_GAMMA  = 1;
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 868a20cbe9da870391c7cc6cf64331a5afb2472d..ad34af7e6e4c38cd31bc63b8e071bf6b1a92ddad 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -1,5 +1,9 @@
 clear all;
 addpath(genpath('../matlab')) % ... add
+SUBMIT = 1; % To submit the job automatically
+% EXECNAME = 'helaz_dbg';
+  EXECNAME = 'helaz_2.62';
+for ETAB = [0.6 0.7]
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -7,37 +11,38 @@ addpath(genpath('../matlab')) % ... add
 CLUSTER.PART  = 'prod';     % dbg or prod
 CLUSTER.TIME  = '24:00:00'; % allocation time hh:mm:ss
 if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME  = '00:30:00'; end;
-CLUSTER.MEM   = '64GB';     % Memory
+CLUSTER.MEM   = '128GB';     % Memory
 CLUSTER.JNAME = 'HeLaZ';% Job name
 NP_P          = 2;          % MPI processes along p  
 NP_KR         = 24;         % MPI processes along kr
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1.0;   % Collision frequency
-ETAB    = 0.6;   % Magnetic gradient
+NU      = 1e-1;   % Collision frequency
+% ETAB    = 0.7;   % Magnetic gradient
 NU_HYP  = 1.0;   % Hyperdiffusivity coefficient
 NL_CLOS = -1;   % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
 % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
-CO      = -3;
+CO      = -2;
+INIT_ZF = 0; ZF_AMP = 0.0;
 %% GRID PARAMETERS
 N       = 200;    % Frequency gridpoints (Nkr = N/2)
 L       = 120;    % Size of the squared frequency domain
-P       = 10;     % Electron and Ion highest Hermite polynomial degree
-J       = 05;     % Electron and Ion highest Laguerre polynomial degree
-MU_P    = 0.0/(P/2)^4;% Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J    = 0.0/(J/2)^2;% Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
+P       = 06;     % Electron and Ion highest Hermite polynomial degree
+J       = 03;     % Electron and Ion highest Laguerre polynomial degree
+MU_P    = 0.0;% Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J    = 0.0;% Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
 %% TIME PARAMETERS
 TMAX    = 5000;  % Maximal time unit
 DT      = 1e-2;  % Time step
 SPS0D   = 1;     % Sampling per time unit for profiler
-SPS2D   = 1/2;     % Sampling per time unit for 2D arrays
+SPS2D   = 1/4;     % Sampling per time unit for 2D arrays
 SPS5D   = 1/100;  % Sampling per time unit for 5D arrays
 SPSCP   = 0;     % Sampling per time unit for checkpoints
-RESTART = 1;     % To restart from last checkpoint
+RESTART = 0;     % To restart from last checkpoint
 JOB2LOAD= 0;
 %% Naming
 % SIMID   = 'test';  % Name of the simulation
-SIMID   = ['v2.5_P_',num2str(P),'_J_',num2str(J)];  % Name of the simulation
+SIMID   = ['v2.6_P_',num2str(P),'_J_',num2str(J)];  % Name of the simulation
 PREFIX  =[];
 % PREFIX  = sprintf('%d_%d_',NP_P, NP_KR);
 %% Options
@@ -66,7 +71,7 @@ JMAXI   = J;     % Highest ''       Laguerre ''
 kmax    = N*pi/L;% Highest fourier mode
 % kmax    = 2/3*N*pi/L;% Highest fourier mode with AA
 HD_CO   = 0.5;    % Hyper diffusivity cutoff ratio
-MU      = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
+MU      = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
 NOISE0  = 1.0e-5;
 ETAT    = 0.0;    % Temperature gradient
 ETAN    = 1.0;    % Density gradient
@@ -82,7 +87,11 @@ CLUSTER.CPUPT = '1';        % CPU per task
 setup
 write_sbash_marconi
 system('rm fort.90 setup_and_run.sh batch_script.sh');
-disp('done');
 if(mod(NP_P*NP_KR,48)~= 0)
     disp('WARNING : unused cores (ntot cores must be a 48 multiple)');
+end
+if(SUBMIT)
+    system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh');
+end
+disp('done');
 end
\ No newline at end of file
diff --git a/wk/open_figure_script.m b/wk/open_figure_script.m
index 2b82c3a0df89b2c8f6be6343e5d5c7615955a2a4..0ec2a71fd41f5bb31c0937017701871e5b7370fb 100644
--- a/wk/open_figure_script.m
+++ b/wk/open_figure_script.m
@@ -1,18 +1,25 @@
 simname_ = '';
+fname_='';
+%% Marconi output file
+fname_='';
+fname_='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/out.txt';
+% fname_='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_2e-02/out.txt';
+simname_ = fname_(54:end-8);
 
 %%
 % simname_ = 'v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.7_nu_1e+00_SGGK_CLOS_0_mu_2e-02';
 % simname_ = 'v2.5_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.8_nu_1e+00_SGGK_CLOS_0_mu_2e-02';
 
-% simname_ = 'v2.5_P_2_J_1/200x100_L_120_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02';
+% simname_ = 'v2.6_P_2_J_1/200x100_L_120_P_2_J_1_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-01';
 % simname_ = 'v2.5_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02';
-simname_ = 'v2.5_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_5e-01_DGGK_CLOS_0_mu_2e-02';
+% simname_ = 'v2.5_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.8_nu_1e-01_DGGK_CLOS_0_mu_2e-02';
+% simname_ = 'v2.6_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02';
 
 
 %%
-% figname_ = '/fig/ZF_transport_drphi_';
+figname_ = '/fig/ZF_transport_drphi_';
 % figname_ = '/fig/space_time_';
-figname_ = '/fig/phi_shear_phase_space_';
+% figname_ = '/fig/phi_shear_phase_space_';
 
 path_    = '../results/';
 
diff --git a/wk/transport_results_2_5.m b/wk/transport_results_2_5.m
index c47aca0dedb0a2a246cee4751a70f0eeb5084db8..46d2362717c3c1aa38acfbb4b9a398490fe553ac 100644
--- a/wk/transport_results_2_5.m
+++ b/wk/transport_results_2_5.m
@@ -1,99 +1,151 @@
 %% nuDGGK = 1.0
-figure
+figure; set(gcf, 'Position',  [100, 100, 1200, 350])
 % 6,3 nuDGGK=1.0
 eta   = [0.5,0.6,0.7,0.8];
 gamma = [32.6443 3.6895 0.3744 0.056];
-std_  = [6.1529 0.7986 0.0493 0.0134];
-errorbar(eta,gamma,std_,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=1.0$',...
+std_g = [6.1529 0.7986 0.0493 0.0134];
+subplot(121)
+errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=1.0$',...
     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(1,:)); hold on;
+subplot(122)
 
 % 10,5 nuDGGK = 1.0
-eta = [0.6,0.7,0.8];
-gamma = [3.5546 0.3917 0.0421];
-std_   = [0.5846 0.0486 0.0203];
-errorbar(eta,gamma,std_,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=1.0$',...
+eta = [0.5 0.6,0.7,0.8];
+gamma = [32.6 3.5546 0.3917 0.0500];
+std_g = [7.7 0.5846 0.0486 0.0088];
+shear = [1.8505 0.60866 0.048249];
+std_s = [0.1599 0.00614 0.0061403];
+subplot(121)
+errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=1.0$',...
+    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(2,:)); hold on;
+subplot(122)
+errorbar(eta(2:end),shear,std_s,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=1.0$',...
     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(2,:)); hold on;
 
 % 12,6 nuDGGK = 1.0
 eta = [0.6];
 gamma = [4.064];
-std_   = [0.7964];
-errorbar(eta,gamma,std_,'-','DisplayName','$P,J=12,06$, $\nu_{DGGK}=1.0$',...
+std_g = [0.7964];
+subplot(121)
+errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=12,06$, $\nu_{DGGK}=1.0$',...
     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(3,:)); hold on;
+subplot(122)
 
 
 % Mix len Ricci 2006
+subplot(121)
 semilogy([0.5  1.0],[10  1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
 set(gca, 'YScale', 'log'); grid on; legend('show')
 xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$') 
+subplot(122)
+grate = [0.2131 0.106 0.02021];
+semilogy([0.6 0.7 0.8],grate,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
+set(gca, 'YScale', 'log'); grid on; legend('show')
+xlabel('$L_n/R$'); ylabel('$V_E$''') 
 
 %% nuDGGK = 0.5
-figure
+figure; set(gcf, 'Position',  [100, 100, 1200, 350])
 % 6,3 nuDGGK=0.5
 eta = [0.5,0.6,0.7,0.8];
 gamma = [20.511 2.6292 0.2353 0.057];
-std_   = [3.67 1.13 0.055 0.004];
-errorbar(eta,gamma,std_,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.5$',...
+std_g = [3.67 1.2 0.055 0.004];
+shear = [nan 1.7417  0.57345 0.25155];
+std_s = [nan 0.35991 0.041 0.00913];
+subplot(121)
+errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.5$',...
+    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(4,:)); hold on;
+subplot(122)
+errorbar(eta,shear,std_s,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.5$',...
     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(4,:)); hold on;
 
 % 10,5 nuDGGK=0.5
-eta    = [0.6,0.7,0.8];
-gamma  = [2.4949 0.2578 0.058];
-std_   = [0.8696 0.1191  0.011];
-errorbar(eta,gamma,std_,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.5$',...
+eta    = [0.6,0.7,0.8 0.9];
+gamma  = [2.4949 0.23368 0.057792 0.023572];
+std_g  = [0.8696 0.085267 0.0060463 0.0046137];
+shear = [1.7094 0.55278 0.2054 0.085678];
+std_s = [0.2428 0.068223 0.012734 0.012291];
+subplot(121)
+errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.5$',...
+    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(5,:)); hold on;
+subplot(122)
+errorbar(eta,shear,std_s,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.5$',...
     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(5,:)); hold on;
 
 % Mix len Ricci 2006
+subplot(121)
 semilogy([0.5  1.0],[10  1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
 set(gca, 'YScale', 'log'); grid on; legend('show')
 xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$') 
+subplot(122)
+grate = [0.2194 0.129 0.05084 0.01346];
+semilogy([0.6 0.7 0.8 0.9],grate,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
+set(gca, 'YScale', 'log'); grid on; legend('show')
+xlabel('$L_n/R$'); ylabel('$V_E$''') 
 
 %% nuDGGK = 0.1
-figure
+figure; set(gcf, 'Position',  [100, 100, 1200, 350])
 % 6,3
 eta =   [0.6        0.7     0.8];
 gamma = [0.24321    0.085   0.0367];
-std_  = [0.295      0.05    0.0023];
-bmax  = [1.0        0.21    0.0367];
-bmin  = [0.02       0.04    0.0367];
+std_g  = [0.295      0.05    0.0023];
+gbmax  = [1.0        0.21    0.0367];
+gbmin  = [0.02       0.04    0.0367];
 % errorbar(eta,gamma,std_,'-','DisplayName','$P,J=06,03$, $\nu_{DGGK}=0.1$',...
 %     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
+subplot(121)
 plot(eta,gamma,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(6,:)); hold on;
-plot(eta,bmax,'-.','DisplayName',' ','Color',line_colors(6,:)); hold on;
-plot(eta,bmin,'-.','DisplayName',' ','Color',line_colors(6,:)); hold on;
+plot(eta,gbmax,'-.','DisplayName',' ','Color',line_colors(6,:)); hold on;
+plot(eta,gbmin,'-.','DisplayName',' ','Color',line_colors(6,:)); hold on;
 
 % 10,5
 eta =   [0.5        0.6        0.7      0.8       0.9];
-gamma = [12.2       0.2133     0.088    0.04253   0.02435];
-std_   = [4.7       0.21328    0.065    0.009     0.00386];
-bmax  = [12.2       0.8        0.25     0.04253   0.02435];
-bmin  = [12.2       0.02       0.03     0.04253   0.02435];
+gamma = [12.2       0.19761 0.088 0.04253 0.026037];
+std_g   = [4.7      0.21328 0.065 0.009 0.00118];
+gbmax  = [12.2      0.8 0.25 0.04253 0.026037];
+gbmin  = [12.2      0.02 0.03 0.04253 0.026037];
+shear = [0.65361 0.46548 0.30645 0.21123 ];
+std_s  = [0.21288 0.10233 0.02886 0.0044664 ];
+sbmax  = [1 0.7 0.30645 0.21123 ];
+sbmin  = [0.4 0.35 0.30645 0.21123 ];
 % errorbar(eta,gamma,std_,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
 %     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
+subplot(121)
 plot(eta,gamma,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
-plot(eta,bmax,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
-plot(eta,bmin,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
+plot(eta,gbmax,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
+plot(eta,gbmin,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
+subplot(122)
+plot(eta(2:end),shear,'-','DisplayName','$P,J=10,05$, $\nu_{DGGK}=0.1$',...
+    'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(7,:)); hold on;
+plot(eta(2:end),sbmax,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
+plot(eta(2:end),sbmin,'-.','DisplayName',' ','Color',line_colors(7,:)); hold on;
+
 % Mix len Ricci 2006
+subplot(121)
 semilogy([0.5  1.0],[10  1e-1],'--','color',[0,0,0]+0.6, 'DisplayName','Mix. Length');
 set(gca, 'YScale', 'log'); grid on; legend('show')
 xlabel('$L_n/R$'); ylabel('$\Gamma^\infty_{part}$') 
+subplot(122)
+grate = [0.236 0.174 0.112 0.053];
+semilogy([0.6 0.7 0.8 0.9],grate,'--','color',[0,0,0]+0.6, 'DisplayName','Lin. Max Growth');
+set(gca, 'YScale', 'log'); grid on; legend('show')
+xlabel('$L_n/R$'); ylabel('$V_E$''') 
 
 %% nuSGGK = 1.0
 figure
 % 6,3 nuDGGK=1.0
 eta = [0.5 0.6,0.7,0.8];
 gamma = [2.3 0.2215 0.0133 0.0032];
-std_   = [3.1 0.22 0.0019 0.0006];
-errorbar(eta,gamma,std_,'-','DisplayName','$P,J=06,03$, $\nu_{SGGK}=1.0$',...
+std_g   = [3.1 0.22 0.0019 0.0006];
+errorbar(eta,gamma,std_g,'-','DisplayName','$P,J=06,03$, $\nu_{SGGK}=1.0$',...
     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(8,:)); hold on;
 
 
 % 10,5 nuDGGK = 1.0
 eta = [0.5 0.6,0.7,0.8];
 gamma = [10 0.319 0.009 0.0026];
-std_   = [1.34 0.228 0.001 0.001];
+std_g   = [1.34 0.228 0.001 0.001];
 % errorbar(eta,gamma,err,'-.','DisplayName','$P,J=10,5$, $\nu_{SGGK}=1.0$',...
 %     'LineWidth',2.0,'MarkerSize',5, 'Color', line_colors(1,:)); hold on;