diff --git a/matlab/plots/plot_radial_transport_and_shear.m b/matlab/plots/plot_radial_transport_and_shear.m
index dc94c2cf94022a1767c5325a2746382368f29c3d..30a3128e1ad8cc04b1a850a75bbf9ef45e380e25 100644
--- a/matlab/plots/plot_radial_transport_and_shear.m
+++ b/matlab/plots/plot_radial_transport_and_shear.m
@@ -15,10 +15,10 @@ shear_infty_std = std (mean(shear_maxx_avgy(:,its2D:ite2D),1));
 fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position',  [100, 100, 1200, 600])
     subplot(311)
 %     yyaxis left
-        plot(Ts0D,PGAMMA_RI*SCALE,'DisplayName','$\langle n_i d\phi/dz \rangle_z$'); hold on;
+        plot(Ts0D,PGAMMA_RI*SCALE,'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
         plot(Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',...
             'DisplayName',['$\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]);
-        grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_r$')
+        grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_x$')
         ylim([0,5*abs(gamma_infty_avg)]); xlim([Ts0D(1),Ts0D(end)]);
         title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
         ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
@@ -29,10 +29,10 @@ fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position',
         lstyle   = line_styles(1);
 %         plt = @(x_) mean(x_,1);
         plt = @(x_) x_(1,:);
-        plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{x,y}(s_\phi)$'); hold on;
-        plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{x}\langle s_\phi\rangle_y$'); hold on;
-        plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{y}\langle s_\phi\rangle_x$'); hold on;
-        plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s_\phi\rangle_{x,y}$'); hold on;
+        plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{x,y}(\partial^2_x\phi)$'); hold on;
+        plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{x}\langle \partial^2_x\phi\rangle_y$'); hold on;
+        plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{y}\langle \partial^2_x\phi\rangle_x$'); hold on;
+        plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle \partial^2_x\phi\rangle_{x,y}$'); hold on;
         plot(Ts3D(its2D:ite2D),ones(ite2D-its2D+1,1)*shear_infty_avg, '-k',...
         'DisplayName',['$s^{\infty} = $',num2str(shear_infty_avg),'$\pm$',num2str(shear_infty_std)]);
         ylim([0,shear_infty_avg*5.0]); xlim([Ts0D(1),Ts0D(end)]);
diff --git a/matlab/plots/plot_space_time_diagrams.m b/matlab/plots/plot_space_time_diagrams.m
index 582f52200a3641454409d0a89f45425c3b04f5f7..7f101ff0d2809e9b12bf4704981a679e0e9b6868 100644
--- a/matlab/plots/plot_space_time_diagrams.m
+++ b/matlab/plots/plot_space_time_diagrams.m
@@ -4,9 +4,9 @@ trange = itstart:itend;
 [TY,TX] = meshgrid(x,Ts3D(trange));
 fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position',  [100, 100, 1200, 600])
     subplot(211)
-%         pclr = pcolor(TX,TY,squeeze(mean(dens_i(:,:,trange).*dyphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
-        pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,1,trange).*dyphi(:,:,1,trange),2))'); set(pclr, 'edgecolor','none',...
-            'DisplayName','$\langle n_i\partial_z\phi\rangle_z$'); colorbar;
+        pclr = pcolor(TX,TY,squeeze(mean(dens_i(:,:,1,trange).*dyphi(:,:,1,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
+%         pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,1,trange).*dyphi(:,:,1,trange),2))'); set(pclr, 'edgecolor','none',...
+%             'DisplayName','$\langle n_i\partial_z\phi\rangle_z$'); colorbar;
         shading interp
         colormap hot;
 %         caxis([0.0,0.05*max(max(mean(ni00(:,:,its2D:ite2D).*dyphi(:,:,1,its2D:ite2D),2)))]);
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index c0bd159881a84245c615436410ac2f52df5f053d..3061a6179193704c5fc4883b26d8606571110667 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -6,7 +6,13 @@ outfile ='';
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='simulation_A/cw_SGGK_mu_1e-2';
+outfile ='';
+outfile ='';
+outfile ='HD_study/200x100_L_200_P_2_J_1_eta_0.6_nu_1e+00_DGGK_CLOS_0_mu_0e+00';
+% outfile ='simulation_B/cw_DGGK'; % to analyse
+% outfile ='simulation_B/cw_SGGK_mu_1e-1';
+% outfile ='simulation_B/cw_DGGK';
+% outfile ='simulation_A/cw_DGGK_eta_0.5';
 % outfile ='simulation_A/LBDK_damping_150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00';
     BASIC.RESDIR      = ['../results/',outfile,'/'];
     BASIC.MISCDIR     = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
@@ -17,15 +23,14 @@ else% Marconi results
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='';
-outfile ='';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt';
 BASIC.RESDIR      = ['../',outfile(46:end-8),'/'];
 BASIC.MISCDIR     = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
 end
 
 %% Load the results
 % Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 06; JOBNUMMAX = 20; 
+JOBNUMMIN = 00; JOBNUMMAX = 20; 
 % JOBNUMMIN = 07; JOBNUMMAX = 20; % For CO damping sim A 
 compile_results %Compile the results from first output found to JOBNUMMAX if existing
 
@@ -44,13 +49,13 @@ end
 
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 1.4e4; TAVG_1 = 1.5e4; % Averaging times duration
+TAVG_0 = 1.2e4; TAVG_1 = 1.3e4; % Averaging times duration
 plot_radial_transport_and_shear
 end
 
 if 0
 %% Space time diagramms
-cmax = 0.0001 % max of the colorbar for transport
+cmax = 0.00001 % max of the colorbar for transport
 tstart = 0; tend = Ts3D(end); % time window
 plot_space_time_diagrams
 end
@@ -72,7 +77,7 @@ end
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
-t0    =10950; iz = 1; ix = 1; iy = 1;
+t0    =0; iz = 1; ix = 1; iy = 1;
 skip_ =1; DELAY = 2e-3*skip_;
 [~, it03D] = min(abs(Ts3D-t0)); FRAMES_3D = it03D:skip_:numel(Ts3D);
 [~, it05D] = min(abs(Ts5D-t0)); FRAMES_5D = it05D:skip_:numel(Ts5D);
@@ -88,7 +93,7 @@ INTERP = 0; T = Ts3D; FRAMES = FRAMES_3D;
 % FIELD = temp_i-Z_T_i; NAME = 'Ti_NZ';FIELDNAME = 'T_i^{NZ}';
 % FIELD = ne00;         NAME = 'ne00'; FIELDNAME = 'n_e^{00}';
 % FIELD = ni00;   NAME = 'ni00'; FIELDNAME = 'n_i^{00}';
-% FIELD = phi;    NAME = 'phi'; FIELDNAME = '\phi';
+FIELD = phi;    NAME = 'phi'; FIELDNAME = '\phi';
 % FIELD = Gamma_x;  NAME = 'Gamma_x'; FIELDNAME = '\Gamma_x';
 
 % Sliced
@@ -129,7 +134,7 @@ FIELD = phi;           FNAME = 'phi';     FIELDLTX = '\phi';
 % FIELD = dens_e-Z_n_e-(Z_phi-phi);       FNAME = 'Non_adiab_part'; FIELDLTX = 'n_e^{NZ}-\phi^{NZ}';
 
 % Chose when to plot it
-tf = 11000:50:11300;
+tf = 128:1:133;
 
 % Sliced
 ix = 1; iy = 1; iz = 1;
diff --git a/wk/continue_multiple_runs_marconi.m b/wk/continue_multiple_runs_marconi.m
index 97346616ee7aabc408182b9c50dbf63fe1b7bc31..fc4fa8b74a005be18aa32f19d6465cb8d6ab8e1f 100644
--- a/wk/continue_multiple_runs_marconi.m
+++ b/wk/continue_multiple_runs_marconi.m
@@ -1,6 +1,5 @@
 %% Paste the list of continue_run calls
-
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_7e-02_SGGK_mu_0e+00/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt')
 
 %% Functions to modify preexisting fort.90 input file and launch on marconi
 function [] = continue_run(outfilename)
@@ -51,7 +50,7 @@ function [] = continue_run(outfilename)
     A{39} = ['  job2load      = ',num2str(J2L)];
     disp(A{39})
     % Change time step
-    A{3} = ['  dt     = 0.005'];
+    A{3} = ['  dt     = 0.01'];
     % Increase endtime
     A{4} = ['  tmax      = 20000'];
     % Change collision operator
@@ -63,7 +62,7 @@ function [] = continue_run(outfilename)
     % change HD
     line_= A{47};
     mu_old = str2num(line_(13:end));
-    A{47} = ['  mu      = ',num2str(mu_old)];
+    A{47} = ['  mu      = ',num2str(0)];
     % change L
     line_= A{14};
     L_old = str2num(line_(12:end));
diff --git a/wk/linear_study.m b/wk/linear_study.m
index 3f174a9689652fcba4772db5ebeb84d855ea70a3..58b08996205a97d8b824f2147f260fe00a917692 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -1,4 +1,4 @@
-for CO = [1]
+for CO = [2]
     RUN = 1; % To run or just to load
 addpath(genpath('../matlab')) % ... add
 default_plots_options
@@ -7,12 +7,12 @@ default_plots_options
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1.0;   % Collision frequency
+NU      = 0.1;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
 ETAB    = 1.0;
 ETAN    = 1/0.6;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
-NU_HYP  = 0.0;   % Hyperdiffusivity coefficient
+NU_HYP  = 3.0;   % Hyperdiffusivity coefficient
 LAMBDAD = 0.0;
 NOISE0  = 1.0e-5; % Init noise amplitude
 BCKGD0  = 0.0;    % Init background
@@ -111,8 +111,8 @@ for i = 1:Nparam
         [~,itend]   = min(abs(Ts3D-tend));
         gamma_Ni00(i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(Ni00(ikx,1,1,itstart:itend))))));
         gamma_phi (i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(PHI (ikx,1,1,itstart:itend))))));
-        Ni00_ST(i,ikx,:) = squeeze((Ni00(ikx,1,1,1:TMAX/SPS3D)));
-         PHI_ST(i,ikx,:) = squeeze((PHI (ikx,1,1,1:TMAX/SPS3D)));
+%         Ni00_ST(i,ikx,:) = squeeze((Ni00(ikx,1,1,1:TMAX/SPS3D)));
+%          PHI_ST(i,ikx,:) = squeeze((PHI (ikx,1,1,1:TMAX/SPS3D)));
     end
     tend   = Ts5D(end); tstart   = 0.4*tend;
     [~,itstart] = min(abs(Ts5D-tstart));
diff --git a/wk/load_multiple_outputs_marconi.m b/wk/load_multiple_outputs_marconi.m
index 4d2acbcb73bf877973a22e069db26cebde95a53a..87fe34405d3dddf761422a0cc24658d362ea2b49 100644
--- a/wk/load_multiple_outputs_marconi.m
+++ b/wk/load_multiple_outputs_marconi.m
@@ -1,10 +1,6 @@
 addpath(genpath('../matlab')) % ... add
 %% Paste the list of simulation results to load
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_3e-02/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_3e-02/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_2e-01_SGGK_mu_0e+00/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_7e-02_SGGK_mu_0e+00/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_SGGK_mu_3e-02/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e+00_SGGK_mu_0e+00/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_DGGK_mu_0e+00/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_DGGK_mu_0e+00/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt')
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 06111db21cc42d0a955e3f3474cca6b8e35c5ae9..1620081fad564ff2af6cd3a96dee5554c77e738f 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -18,28 +18,28 @@ NP_P          = 2;          % MPI processes along p
 NP_KX         = 24;         % MPI processes along kx
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1.0;   % Collision frequency
+NU      = 0.5;   % Collision frequency
 ETAN    = 1.0/0.6;    % Density gradient drive (R/Ln)
 NU_HYP  = 0.0;
 %% GRID PARAMETERS
-N       = 150;     % Frequency gridpoints (Nkx = N/2)
-L       = 100;     % Size of the squared frequency domain
+N       = 300;     % Frequency gridpoints (Nkx = N/2)
+L       = 120;     % Size of the squared frequency domain
 Nz      = 1;      % number of perpendicular planes (parallel grid)
 q0      = 1.0;    % q factor ()
 shear   = 0.0;    % magnetic shear
 eps     = 0.0;    % inverse aspect ratio
-P       = 4;
-J       = 2;
+P       = 8;
+J       = 4;
 %% TIME PARAMETERS
 TMAX    = 10000;  % Maximal time unit
 DT      = 5e-3;   % Time step
 SPS0D   = 1;      % Sampling per time unit for profiler
 SPS2D   = 1;      % Sampling per time unit for 2D arrays
 SPS3D   = 1/2;      % Sampling per time unit for 3D arrays
-SPS5D   = 1/500;  % Sampling per time unit for 5D arrays
+SPS5D   = 1/100;  % Sampling per time unit for 5D arrays
 SPSCP   = 0;    % Sampling per time unit for checkpoints/10
 RESTART = 1;      % To restart from last checkpoint
-JOB2LOAD= 3;
+JOB2LOAD= 0;
 %% OPTIONS AND NAMING
 % Collision operator
 % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK)
@@ -47,7 +47,7 @@ CO      = 2;
 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   = 'test_3D_marconi';  % Name of the simulation
-SIMID   = 'HD_study';  % Name of the simulation
+SIMID   = 'simulation_B';  % Name of the simulation
 % SIMID   = ['v3.0_P_',num2str(P),'_J_',num2str(J)];  % Name of the simulation
 NON_LIN = 1;   % activate non-linearity (is cancelled if KXEQ0 = 1)
 % INIT options
@@ -79,6 +79,7 @@ HD_CO   = 0.5;    % Hyper diffusivity cutoff ratio
 % kmaxcut = 2.5;
 MU      = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
 NOISE0  = 1.0e-5;
+BCKGD0  = 0.0;    % Init background
 TAU     = 1.0;    % e/i temperature ratio
 ETAT    = 0.0;    % Temperature gradient
 ETAB    = 1.0;    % Magnetic gradient (1.0 to set R=LB)
diff --git a/wk/open_figure_script.m b/wk/open_figure_script.m
index 48416378d5fdc9a6bcdb36305a2431a96e0bcbea..51c651f3b1412bfd754522c2b3cfc3910729a3b5 100644
--- a/wk/open_figure_script.m
+++ b/wk/open_figure_script.m
@@ -15,7 +15,7 @@ simname_ = fname_(54:end-8);
 % simname_ = '';
 % simname_ = '';
 % simname_ = '';
-simname_ = 'simulation_A/DGGK_damping_150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00';
+simname_ = 'simulation_A/cw_SGGK_mu_1e-2';
 
 
 
diff --git a/wk/plot_cosol_mat.m b/wk/plot_cosol_mat.m
index 82a9306680170474ab662c7b277779d3671282a7..60c8e489fa1f1050f0bd5949b6bedbe7836d045c 100644
--- a/wk/plot_cosol_mat.m
+++ b/wk/plot_cosol_mat.m
@@ -92,12 +92,8 @@ subplot(224)
     
 %% Single eigenvalue analysis
 
-% mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_2_J_1_N_1_kpm_4.0_NFLR_5.h5';
-% mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_8.h5';
-% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb.NFLR/k.4/self.NFLR.8.0.h5';
-% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_6_3_NFLR_8_kp_4_npiflr_8_sigma_3e-3/scanfiles_00000/self.0.h5';
-% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_6_3_NFLR_8_kp_4_npiflr_0_sigma_5e-4/scanfiles_00000/self.0.h5';
-mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_6_3_NFLR_8_kp_4_bjf/scanfiles_00000/self.0.h5';
+% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_6_3_NFLR_8_kp_4_bjf/scanfiles_00000/self.0.h5';
+mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_6_3_NFLR_8_kp_4_new_Tljpmf_fort_bjf/scanfiles_00000/self.0.h5';
 matidx = 74;
 
 matidx = sprintf('%5.5i',matidx);disp(matidx);