diff --git a/Makefile b/Makefile
index e808a162333e833a24663374a9e3f106f39f03f8..d8f9cd131c06c5fa675444dbd4904947c55b313c 100644
--- a/Makefile
+++ b/Makefile
@@ -147,7 +147,8 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
  $(OBJDIR)/auxval.o : src/auxval.F90 \
  	 $(OBJDIR)/fourier_mod.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o \
 	 $(OBJDIR)/geometry_mod.o  $(OBJDIR)/grid_mod.o $(OBJDIR)/numerics_mod.o \
-	 $(OBJDIR)/parallel_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/CLA_mod.o
+	 $(OBJDIR)/parallel_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/CLA_mod.o \
+	 $(OBJDIR)/ExB_shear_flow_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/auxval.F90 -o $@
 
  $(OBJDIR)/basic_mod.o : src/basic_mod.F90 \
@@ -227,7 +228,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/geometry_mod.F90 -o $@
 
  $(OBJDIR)/ghosts_mod.o : src/ghosts_mod.F90 \
- 	 $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o \
+ 	 $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o\
 	 $(OBJDIR)/geometry_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ghosts_mod.F90 -o $@
 
@@ -256,7 +257,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
 
  $(OBJDIR)/memory.o : src/memory.F90 $ \
  	 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/collision_mod.o\
-   $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o \
+     $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o \
 	 $(OBJDIR)/grid_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/memory.F90 -o $@
 
@@ -327,7 +328,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
 	 $(OBJDIR)/array_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/fields_mod.o \
 	 $(OBJDIR)/ghosts_mod.o $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/solve_EM_fields.o\
 	 $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o \
-	 $(OBJDIR)/CLA_mod.o
+	 $(OBJDIR)/CLA_mod.o $(OBJDIR)/ExB_shear_flow_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/stepon.F90 -o $@
 
  $(OBJDIR)/tesend.o : src/tesend.F90 \
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index b6c70efdb48a65b914f9ee8424cee76b1cb6dd4c..b9aa6fe0b3d9c109569a0881eb77c47897d42823 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -1,4 +1,4 @@
-function [FIGURE] = mode_growth_meter(DATA,OPTIONS)
+function [FIGURE, kykx, wkykx, ekykx] = mode_growth_meter(DATA,OPTIONS)
 
 NORMALIZED = OPTIONS.NORMALIZED;
 Nma   = OPTIONS.NMA; %Number moving average
@@ -30,7 +30,9 @@ TW = [OPTIONS.KY_TW; OPTIONS.KX_TW];
 
 [~,ikzf] = max(mean(squeeze(abs(field(1,:,FRAMES))),2));
 
+% Amplitude ratio method to measure growth rates and freq.
 [wkykx, ekykx] = compute_growth_rates(DATA.PHI(:,:,:,FRAMES),DATA.Ts3D(FRAMES));
+kykx = meshgrid(DATA.grids.ky,DATA.grids.kx)';
 
 FIGURE.fig = figure; %set(gcf, 'Position',  [100 100 1200 700])
 FIGURE.FIGNAME = 'mode_growth_meter';
@@ -95,35 +97,13 @@ for i = 1:2
 
     gamma = MODES;
     amp   = MODES;
-    % w_ky = zeros(numel(MODES),numel(FRAMES)-1);
-    % ce   = zeros(numel(MODES),numel(FRAMES));
     i_=1;
+    % Alternative method to measure growth
     for ik = MODES
-
         to_measure = log(plt(field,ik));
         gr = polyfit(t(it1:it2),to_measure(it1:it2),1);
         gamma(i_) = gr(1);
         amp(i_)   = gr(2);
-
-        % % Second method for measuring growth rate (measures frequ. too)
-        % if MODES_SELECTOR == 1
-        %     to_measure = reshape(DATA.PHI(ik,1,:,FRAMES),DATA.grids.Nz,numel(FRAMES));
-        % else
-        %     to_measure = reshape(DATA.PHI(1,ik,:,FRAMES),DATA.grids.Nz,numel(FRAMES));
-        % end
-        % for it = 2:numel(FRAMES)
-        %     phi_n   = to_measure(:,it); 
-        %     phi_nm1 = to_measure(:,it-1);
-        %     dt      = t(it)-t(it-1);
-        %     ZS      = sum(phi_nm1,1); %sum over z
-        % 
-        %     wl          = log(phi_n./phi_nm1)/dt;
-        %     w_ky(i_,it) = squeeze(sum(wl.*phi_nm1,1)./ZS);
-        % 
-        %     % for iky = 1:numel(w_ky(:,it))
-        %     %     ce(iky,it)   = abs(sum(abs(w_ky(iky,it)-wl(iky,:)).^2.*phi_nm1(iky,:),2)./ZS(iky,:));
-        %     % end
-        % end
         i_=i_+1;
     end
     mod2plot = 2:min(Nmax,OPTIONS.NMODES+1);
@@ -135,7 +115,7 @@ for i = 1:2
         pclr = pcolor(XX,YY,abs(plt(fftshift(field,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none');  hold on;
         set(gca,'YDir','normal')
     %     xlim([t(1) t(end)]); %ylim([1e-5 1])
-        xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /\rho_s$');
+        xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /R_0$');
         title([MODESTR,', ',zstrcomp])  
 
 %     subplot(2+d,3,2+3*(i-1))
@@ -151,13 +131,24 @@ for i = 1:2
         plot(t(it1)*[1 1],[1e-10 1],'--k')
         plot(t(it2)*[1 1],[1e-10 1],'--k')
         xlim([t(1) t(end)]); %ylim([1e-5 1])
-        xlabel('$t c_s /\rho_s$'); ylabel(['$|\phi_{',kstr,'}|$']);
+        xlabel('$t c_s /R_0$'); ylabel(['$|\phi_{',kstr,'}|$']);
         title('Measure time window')
 
 %     subplot(2+d,3,3+3*(i-1))
     FIGURE.axes(3+3*(i-1)) = subplot(2+d,3,3+3*(i-1),'parent',FIGURE.fig);
     % yyaxis("left")
-        errorbar(k(MODES),real(omega(ik)),real(err(ik)),'-k',...
+    if OPTIONS.GOK2
+        x_ = k(MODES);
+        y_ = real(omega(ik))./k(MODES).^2;
+        e_ = real(err(ik))./k(MODES).^2;
+        lb_= '\gamma/k^2';
+    else
+        x_ = k(MODES);
+        y_ = real(omega(ik));
+        e_ = real(err(ik));
+        lb_= '\gamma';
+    end
+        errorbar(x_,y_,e_,'-ok',...
                 'LineWidth',1.5,...
                 'DisplayName',...
                 ['$\gamma$, (',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']); 
@@ -170,12 +161,12 @@ for i = 1:2
         end
         % ylabel('$\gamma$');
     % yyaxis("right")
-        errorbar(k(MODES),imag(omega(ik)),imag(err(ik)),'--k',...
+        errorbar(k(MODES),imag(omega(ik)),imag(err(ik)),'--ok',...
                 'LineWidth',1.5,...
                 'DisplayName',...
                 ['$\omega$, (',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']); 
         hold on;
-        ylabel('$\gamma,\omega$');
+        ylabel(['$',lb_,',\omega$']);
         xlabel(['$',kstr,'\rho_s$']);
         title('Growth rates')
 end
diff --git a/matlab/plot/quick_lin_gene_plot.m b/matlab/plot/quick_lin_gene_plot.m
new file mode 100644
index 0000000000000000000000000000000000000000..2a5a42195f8328bb04daa29d092fe160a22f605f
--- /dev/null
+++ b/matlab/plot/quick_lin_gene_plot.m
@@ -0,0 +1,112 @@
+%% gyroLES plot
+% genepath = '/misc/gene_results';
+% % simpath  = '/NL_Zpinch_Kn_1.7_eta_0.25_nuSG_1e-1_gyroLES/';
+% simpath  = '/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_gyroLES/';
+% filename = 'GyroLES_ions.dat';
+% toload = [genepath simpath filename];
+% data = readtable(toload);
+% t = data.Var1;
+% mu_x = data.Var2;
+% mu_y = data.Var3;
+% figure;
+% plot(t,mu_x); hold on
+% plot(t,mu_y);
+% legend('mu_x','mu_y'); xlabel('time'); ylabel('Hyper-diff');
+% title('gyroLES results for H&P fig. 2c')
+
+%% Plot linear results
+% path = '/home/ahoffman/gene/linear_zpinch_results/';
+% fname ='tmp.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSG_0.047_32x16.txt';
+% fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.047_32x16.txt';
+% fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.0235_32x16.txt';
+% fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.0047_32x16.txt';
+% fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.0047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.4_KT_0.35_nuSG_0.0047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.6_KT_0.4_nuSG_0.0047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.6_KT_0.4_nuSGDK_0.0047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.6_KT_0.4_nuLDDK_0.0047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.4_nuLDDK_0.0047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.45_nu_0_32x16.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSG_0.0235_64x32.txt';
+% fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.235_64x32.txt';
+% fname ='GENE_LIN_Kn_1.7_KT_0.425_nuSG_0.235_64x32.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSGDK_0.047_32x16.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.475_nu_0_mu_5e-2.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSGDK_0.0235_64x32.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSGDK_0.0235_32x16.txt';
+% fname ='GENE_LIN_Kn_2.0_KT_0.5_nu_0_32x16.txt';
+% fname ='GENE_LIN_Kn_2.0_KT_0.5_nuSGDK_0.0235_32x16.txt';
+% fname ='GENE_LIN_Kn_1.6_KT_0.4_nu_0_32x16.txt';
+% fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt';
+path = '/home/ahoffman/gene/linear_CBC_results/';
+% fname = 'CBC_100_20x1x32x30x14_Lv_3_Lw_12_circ.txt';
+% fname = 'CBC_100_20x1x32x32x12_Lv_3_Lw_12.txt';
+% fname = 'CBC_KT_4_20x1x32x32x12_Lv_3_Lw_12.txt';
+% fname = 'CBC_KT_4_20x1x32x64x24_Lv_6_Lw_24.txt';
+% fname = 'CBC_KT_5.3_20x1x32x32x12_Lv_3_Lw_12.txt';
+% fname = 'CBC_KT_5.3_32x1x48x40x16_Lv_3_Lw_12.txt';
+% fname = 'CBC_ky_0.3_20x1x32x32x12_Lv_3_Lw_12.txt';
+% fname = 'CBC_ky_0.3_20x1x32x32x12_Lv_3_Lw_12_nuv_1e-3.txt';
+% fname = 'CBC_KT_11_20x1x16x24x10_Lv_3_Lw_12.txt';
+% fname = 'CBC_KT_11_20x1x32x30x14_Lv_3_Lw_12.txt';
+% fname = 'CBC_ky_0.3_20x1x16x24x10_Lv_3_Lw_12_nuv_1e-3.txt';
+%----------Shearless CBC
+% fname = 'CBC_salpha_s0_nz_24_nv_48_nw_16_adiabe.txt';
+% fname = 'CBC_salpha_s0_nz_24_nv_48_nw_16_kine.txt';
+% fname = 'CBC_salpha_s0_nz_24_nv_48_nw_16_kine_beta_1e-4.txt';
+% fname = 'CBC_miller_s0_nz_24_nv_48_nw_16_adiabe.txt';
+% fname = 'CBC_miller_s0_nz_24_nv_48_nw_16_kine.txt';
+%----------Shearless pITG
+% fname = 'pITG_salpha_s0_nz_24_nv_48_nw_16_adiabe.txt';
+% fname = 'pITG_miller_s0_nz_24_nv_48_nw_16_adiabe.txt';
+% fname = 'pITG_salpha_s0_nz_24_nv_48_nw_16_kine.txt';
+% fname = 'pITG_miller_s0_nz_24_nv_48_nw_16_kine.txt';
+%----------Convergence nvpar shearless pITG
+% fname = 'pITG_salpha_s0_nz_24_nv_scan_nw_16_adiabe.txt';
+% fname = 'pITG_miller_s0_nz_24_nv_scan_nw_16_adiabe.txt';
+% fname = 'pITG_salpha_s0_nz_24_nv_scan_nw_16_kine.txt';
+% fname = 'pITG_miller_s0_nz_24_nv_scan_nw_16_kine.txt';
+% fname = 'pITG_salpha_s0_nz_24_nv_scan_nw_24_adiabe.txt';
+% fname = 'pITG_miller_s0_nz_24_nv_scan_nw_24_adiabe.txt';
+%----------Convergence nvpar shearless CBC
+% fname = 'CBC_salpha_nz_24_nv_scan_nw_16_adiabe.txt';
+% fname = 'CBC_miller_nz_24_nv_scan_nw_16_adiabe.txt';
+% fname = 'CBC_salpha_nz_24_nv_scan_nw_16_kine.txt';
+% fname = 'CBC_miller_nz_24_nv_scan_nw_16_kine.txt';
+%---------- CBC
+% fname = 'CBC_salpha_nx_8_nz_18_nv_12_nw_8_adiabe.txt';
+% fname = 'CBC_salpha_nx_8_nz_18_nv_18_nw_8_adiabe.txt';
+
+% fname = 'CBC_salpha_nx_8_nz_24_nv_6_nw_4_adiabe.txt';
+% fname = 'CBC_salpha_nx_8_nz_24_nv_6_nw_4_adiabe_0.5L.txt';
+% fname = 'CBC_salpha_nx_8_nz_24_nv_8_nw_4_adiabe.txt';
+% fname = 'CBC_salpha_nx_8_nz_24_nv_8_nw_4_adiabe_0.5L.txt';
+fname = 'CBC_salpha_nx_8_nz_24_nv_8_nw_4_adiabe_0.2L.txt';
+% fname = 'CBC_salpha_nx_8_nz_24_nv_16_nw_8_adiabe.txt';
+% fname = 'CBC_salpha_nx_8_nz_24_nv_16_nw_8_adiabe_0.5L.txt';
+% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
+% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_adiabe_0.5L.txt';
+% fname = 'CBC_salpha_nx_8_nz_24_nv_48_nw_16_adiabe.txt';
+
+% fname = 'kT_5.3_salpha_nx_8_nz_24_nv_60_nw_30_adiabe.txt';
+
+% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_kine.txt';
+% fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_adiabe.txt';
+% fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
+% fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_kine.txt';
+% fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_kine.txt';
+%----------Convergence nv CBC
+% fname = 'CBC_ky_0.3_nv_scan_8x1x24_nw_8_Lv_3_Lw_6.txt';
+% fname = 'CBC_ky_0.3_nv_scan_8x1x24_nw_16_Lv_3_Lw_6.txt';
+% fname = 'CBC_ky_0.3_nv_scan_8x1x24_nw_24_Lv_3_Lw_6.txt';
+
+data_ = load([path,fname]);
+
+figure
+% plot(data_(:,2),data_(:,3),'-dk','DisplayName',fname); hold on;
+% plot(data_(:,2),data_(:,4),'--*k','DisplayName',fname);
+
+plot(data_(:,2),data_(:,3)./data_(:,2).^2,'-dk','DisplayName',fname); hold on;
+legend show
diff --git a/src/ExB_shear_flow_mod.F90 b/src/ExB_shear_flow_mod.F90
index 508449fe46fab7c9e8bc2e69671db40b40b16e80..649f99664d72a5f4fee7e762ebd8f9b4eddef064 100644
--- a/src/ExB_shear_flow_mod.F90
+++ b/src/ExB_shear_flow_mod.F90
@@ -39,7 +39,7 @@ CONTAINS
     SUBROUTINE Apply_ExB_shear_flow
         USE basic,      ONLY: chrono_ExBs, start_chrono, stop_chrono
         USE grid,       ONLY: local_nky, kxarray, update_grids, &
-            local_nkx, deltakx, kx_min, kx_max
+            total_nkx, deltakx, kx_min, kx_max
         USE prec_const, ONLY: PI
         USE geometry,   ONLY: gxx,gxy,inv_hatB2, evaluate_magn_curv
         USE numerics,   ONLY: evaluate_EM_op, evaluate_kernels
@@ -55,13 +55,17 @@ CONTAINS
             IF(shiftnow_ExB(iky)) THEN
                 print*, "SHIFT ARRAYS"
                 ! shift all fields
-                DO ikx = 1,local_nkx
+                DO ikx = 1,total_nkx
                     ikx_s = ikx + jump_ExB(iky)
-                    IF( (kxarray(iky,ikx) .GE. kx_min) .AND. (kxarray(iky,ikx) .LE. kx_max) ) THEN
+                    ! We test if the shifted modes are still in contained in our resolution
+                    ! IF( (kxarray(iky,ikx)-sky_ExB(iky) .GE. kx_min) .AND. (kxarray(iky,ikx)-sky_ExB(iky) .LE. kx_max) ) THEN
+                    IF ( ((ikx_s .GT. 0 ) .AND. (ikx_s .LE. total_nkx )) .AND. &
+                         (((ikx   .LE. (total_nkx/2+1)) .AND. (ikx_s .LE. (total_nkx/2+1))) .OR. &
+                          ((ikx   .GT. (total_nkx/2+1)) .AND. (ikx_s .GT. (total_nkx/2+1)))) ) THEN
                         moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:)
                         phi(iky,ikx,:)             = phi(iky,ikx_s,:)
                         psi(iky,ikx,:)             = psi(iky,ikx_s,:)
-                    ELSE
+                    ELSE ! if it is not, it is lost (~dissipation for high modes)
                         moments(:,:,:,iky,ikx,:,:) = 0._xp
                         phi(iky,ikx,:)             = 0._xp
                         psi(iky,ikx,:)             = 0._xp
@@ -76,9 +80,9 @@ CONTAINS
         CALL update_grids(sky_ExB,gxx,gxy,inv_hatB2)
 
         !   update the EM op., the kernels and the curvature op.
-        CALL evaluate_kernels
-        CALL evaluate_EM_op
-        CALL evaluate_magn_curv
+        ! CALL evaluate_kernels
+        ! CALL evaluate_EM_op
+        ! CALL evaluate_magn_curv
 
         CALL stop_chrono(chrono_ExBs)
     END SUBROUTINE Apply_ExB_shear_flow
diff --git a/wk/lin_run_script.m b/wk/lin_run_script.m
index aad2bfd6be684b33e912e0f3eac6868196955e1e..b4832e7da8a2dc06bf8c925f311142352b0c82c2 100644
--- a/wk/lin_run_script.m
+++ b/wk/lin_run_script.m
@@ -15,7 +15,7 @@ addpath(genpath([gyacomodir,'matlab/load']))    % Add load folder
 addpath(genpath([gyacomodir,'wk/parameters']))  % Add parameters folder
 
 %% Setup run or load an executable
-RUN = 1; % To run or just to load
+RUN = 0; % To run or just to load
 default_plots_options
 % EXECNAME = 'gyacomo23_sp_save'; % single precision
 EXECNAME = 'gyacomo23_sp'; % single precision
@@ -32,15 +32,15 @@ EXECNAME = 'gyacomo23_sp'; % single precision
 run lin_ITG
 % run lin_KBM
 %% Change parameters
-% EXBRATE = 0.001;              % Background ExB shear flow
-% NY   = 2;
-% NX   = 4;
-% PMAX = 2;
-% JMAX = PMAX/2;
-% ky   = 0.5;
-% LY   = 2*pi/ky;
-% DT   = 1e-3;
-% TMAX = 10;
+EXBRATE = 0.0;              % Background ExB shear flow
+NY   = 40;
+NX   = 8;
+PMAX = 16;
+JMAX = PMAX/2;
+ky   = 0.05;
+LY   = 2*pi/ky;
+DT   = 5e-3;
+TMAX = 50;
 % % SIGMA_E = 0.04;
 % TMAX     = 10;
 % DTSAVE0D = 200*DT;
@@ -52,9 +52,9 @@ setup
 % Run linear simulation
 if RUN
     MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
-    RUN  =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
+    % RUN  =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
    % RUN  =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
-     % RUN  =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;'];
+     RUN  =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;'];
     % RUN  =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
       % RUN = ['./../../../bin/gyacomo23_sp 0;'];
     MVOUT='cd ../../../wk;';
@@ -86,8 +86,9 @@ options.NMA    = 1; % Set NMA option to 1
 options.NMODES = 999; % Set how much modes we study
 options.iz     = 'avg'; % Compressing z
 options.ik     = 1; %
+options.GOK2   = 0; % plot gamma/k^2
 options.fftz.flag = 0; % Set fftz.flag option to 0
-fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
+[fig, kykx, wkykx, ekykx] = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
 end
 
 if (1 && NZ>4)
diff --git a/wk/load_and_replot.m b/wk/load_and_replot.m
new file mode 100644
index 0000000000000000000000000000000000000000..3d0e7ee51fd0f3156026196434d16f4f95057fb9
--- /dev/null
+++ b/wk/load_and_replot.m
@@ -0,0 +1,23 @@
+% Load the .fig file
+figFile = '/home/ahoffman/paper2_fig/31_CBC_linear_convergence.fig'; % Replace with the actual path to your .fig file
+% figFile = 'tmp.fig'; % Replace with the actual path to your .fig file
+figHandle = openfig(figFile);
+
+%% Extract data from the figure
+% Here, we assume that the plotted data is stored in a single line plot
+lineHandle = findobj(figHandle, 'Type', 'line');
+xData = get(lineHandle, 'XData');
+yData = get(lineHandle, 'YData');
+% Extract legend labels
+legendHandle = findobj(figHandle, 'Type', 'legend');
+legendEntries = flip(get(legendHandle, 'String'));
+% legendEntries = legendEntries{end:-1:1};
+% Close the loaded figure
+close(figHandle);
+%%
+    figure
+for i = numel(xData):-1:1
+    % plot(xData{i},yData{i},'DisplayName',legendEntries{i}); hold on
+    plot(xData{i},yData{i}./xData{i}.^2,'DisplayName',legendEntries{i}); hold on
+end
+legend show
\ No newline at end of file
diff --git a/wk/load_metadata_scan.m b/wk/load_metadata_scan.m
index f107d4afec12cbd35cd8844c282d73a5cc056d90..7ef57a678dfc60073e788d78e736c638f4d0a4be 100644
--- a/wk/load_metadata_scan.m
+++ b/wk/load_metadata_scan.m
@@ -15,7 +15,11 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
 % datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_1_P_2_8_DGGK_0.05_be_0.0039.mat';
 % datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_0.9_P_2_8_kN_1.33_DGGK_0.56901_be_0.0039.mat';
 % datafname = 'lin_DTT_AB_rho85_PT_scan/8x24_ky_0.05_1.5_P_2_8_kN_1.33_DGGK_0.1_be_0.0039.mat';
-datafname = 'lin_JET_rho97_scan/8x32_ky_0.031623_31.6228_P_1_1_kN_31.4286_DGGK_0.1_be_0.0031.mat';
+% datafname = 'lin_JET_rho97_scan/4x32_ky_0.01_10_P_2_8_kN_10_DGGK_0.1_be_0.0031.mat';
+% datafname = 'lin_JET_rho97_scan/4x32_ky_0.01_10_P_2_8_kN_10_DGGK_0.05_be_0.0031.mat';
+% datafname = 'lin_JET_rho97_scan/8x32_ky_0.01_10_P_2_8_kN_10_DGGK_0.05_be_0.0031.mat';
+% datafname = 'lin_JET_rho97_scan/8x32_ky_0.01_10_P_2_2_kN_10_DGGK_0.05_be_0.0031.mat';
+datafname = 'lin_JET_rho97_scan/8x32_ky_0.031623_31.6228_P_2_8_kN_10_DGGK_0.1_be_0.0031.mat';
 %% Chose if we filter gamma>0.05
 FILTERGAMMA = 1;