diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 8e367c4ef9d3073c134761ba5f24c4263e9cfeee..2200bbb05eec22902a58a33da4e809c1c157b027 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -11,6 +11,7 @@ DT_EVOL  = []; %
 % FIELDS
 Nipj_    = []; Nepj_    = [];
 Ni00_    = []; Ne00_    = [];
+HFLUX_   = [];
 GGAMMA_  = [];
 PGAMMA_  = [];
 PHI_     = [];
@@ -74,10 +75,17 @@ while(CONTINUE)
             end            
         end
         
+        if W_GAMMA || W_HF
+            Ts0D_   = cat(1,Ts0D_,Ts0D);
+        end      
+        
         if W_GAMMA
             GGAMMA_ = cat(1,GGAMMA_,GGAMMA_RI);
             PGAMMA_ = cat(1,PGAMMA_,PGAMMA_RI);
-            Ts0D_   = cat(1,Ts0D_,Ts0D);
+        end
+        
+        if W_HF
+            HFLUX_ = cat(1,HFLUX_,HFLUX_X);
         end
         
         if W_PHI || W_NA00
@@ -131,11 +139,11 @@ while(CONTINUE)
     end
     JOBNUM   = JOBNUM + 1;
 end
-GGAMMA_RI = GGAMMA_; PGAMMA_RI = PGAMMA_; Ts0D = Ts0D_;
+GGAMMA_RI = GGAMMA_; PGAMMA_RI = PGAMMA_; HFLUX_X = HFLUX_; Ts0D = Ts0D_;
 Nipj = Nipj_; Nepj = Nepj_; Ts5D = Ts5D_;
 Ni00 = Ni00_; Ne00 = Ne00_; PHI = PHI_; Ts3D = Ts3D_;
 DENS_E = DENS_E_; DENS_I = DENS_I_; TEMP_E = TEMP_E_; TEMP_I = TEMP_I_;
-clear Nipj_ Nepj_ Ni00_ Ne00_ PHI_ Ts2D_ Ts5D_ GGAMMA_ PGAMMA_ Ts0D_
+clear Nipj_ Nepj_ Ni00_ Ne00_ PHI_ Ts2D_ Ts5D_ GGAMMA_ PGAMMA_ Ts0D_ HFLUX_
 
 Sipj = Sipj_; Sepj = Sepj_;
 clear Sipj_ Sepj_
diff --git a/matlab/create_gif.m b/matlab/create_gif.m
index 0f49d65cc61f09e2d28959743164b02ddad1bcfb..1d0d5717940267f0b3f3343c10e3fc2632250a72 100644
--- a/matlab/create_gif.m
+++ b/matlab/create_gif.m
@@ -17,7 +17,9 @@ else
 % Setup figure frame
 fig  = figure('Color','white','Position', [100, 100, 400, 400]);
     pcolor(X,Y,FIELD(:,:,1)); % to set up
+    if BWR
     colormap(bluewhitered)
+    end
     axis tight manual % this ensures that getframe() returns a consistent size
     if INTERP
         shading interp;
@@ -45,7 +47,9 @@ fig  = figure('Color','white','Position', [100, 100, 400, 400]);
             shading interp; 
         end
         set(pclr, 'edgecolor','none'); axis square;
+        if BWR
         colormap(bluewhitered)
+        end
         xlabel(XNAME); ylabel(YNAME); %colorbar;
         drawnow 
         % Capture the plot as an image 
diff --git a/matlab/load_results.m b/matlab/load_results.m
index fd7e85d694289c9dd6d6fc4eb36e31ff47c6101f..4ccc0671856229edefdc121e6cab796f27b83904 100644
--- a/matlab/load_results.m
+++ b/matlab/load_results.m
@@ -7,12 +7,13 @@ DT_SIM    = h5readatt(filename,'/data/input','dt');
 [Pe, Je, Pi, Ji, kx, ky, z] = load_grid_data(filename);
 
 W_GAMMA   = strcmp(h5readatt(filename,'/data/input','write_gamma'),'y');
-W_PHI     = strcmp(h5readatt(filename,'/data/input','write_phi')  ,'y');
-W_NA00    = strcmp(h5readatt(filename,'/data/input','write_Na00') ,'y');
-W_NAPJ    = strcmp(h5readatt(filename,'/data/input','write_Napj') ,'y');
-W_SAPJ    = strcmp(h5readatt(filename,'/data/input','write_Sapj') ,'y');
-W_DENS    = strcmp(h5readatt(filename,'/data/input','write_dens') ,'y');
-W_TEMP    = strcmp(h5readatt(filename,'/data/input','write_temp') ,'y');
+W_HF      = strcmp(h5readatt(filename,'/data/input','write_hf'   ),'y');
+W_PHI     = strcmp(h5readatt(filename,'/data/input','write_phi'  ),'y');
+W_NA00    = strcmp(h5readatt(filename,'/data/input','write_Na00' ),'y');
+W_NAPJ    = strcmp(h5readatt(filename,'/data/input','write_Napj' ),'y');
+W_SAPJ    = strcmp(h5readatt(filename,'/data/input','write_Sapj' ),'y');
+W_DENS    = strcmp(h5readatt(filename,'/data/input','write_dens' ),'y');
+W_TEMP    = strcmp(h5readatt(filename,'/data/input','write_temp' ),'y');
 
 
 if W_GAMMA
@@ -20,6 +21,10 @@ if W_GAMMA
       PGAMMA_RI              = load_0D_data(filename, 'pflux_ri');
 end
 
+if W_HF
+    [ HFLUX_X, Ts0D, dt0D] = load_0D_data(filename, 'hflux_x');
+end
+
 if W_PHI
     [ PHI, Ts3D, dt3D] = load_3D_data(filename, 'phi');
 end
diff --git a/matlab/plots/plot_time_evolution_and_gr.m b/matlab/plots/plot_time_evolution_and_gr.m
index 4427e3f250c2d32bfeb0d892169a9a9ccf1b98ad..17dbbb666db7a7d84171f9019ada51459bb86c64 100644
--- a/matlab/plots/plot_time_evolution_and_gr.m
+++ b/matlab/plots/plot_time_evolution_and_gr.m
@@ -29,11 +29,12 @@ subplot(111);
     end
     grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); xlabel('$t c_s/R$')
     subplot(222)
-%         plot(Ts0D,GGAMMA_RI*(2*pi/Nx/Ny)^2); hold on;
         yyaxis left; ylabel('$\Gamma_x$');
-        plot(Ts3D,squeeze(sum(sum(sum(Gamma_x,1),2),3)));
+%         plot(Ts3D,squeeze(sum(sum(sum(Gamma_x,1),2),3)));
+        plot(Ts0D,PGAMMA_RI*SCALE); hold on;
         yyaxis right; ylabel('$Q_x$');
-        plot(Ts3D,squeeze(sum(sum(sum(Q_x,1),2),3)));
+%         plot(Ts3D,squeeze(sum(sum(sum(Q_x,1),2),3)));
+        plot(Ts0D,HFLUX_X*SCALE); hold on;
         grid on; xlabel('$t c_s/R$');
     if(~isnan(max(max(g_I(1,:,:)))))
     subplot(223)
diff --git a/matlab/profiler.m b/matlab/profiler.m
index 879cd72eb6c0e693844d5997b650fa792bcb5424..2988e22c285184edcc1afd027e79118ffe67e0d0 100644
--- a/matlab/profiler.m
+++ b/matlab/profiler.m
@@ -1,6 +1,7 @@
 %% load profiling
 % filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00);
-filename_ = ['/misc/HeLaZ_outputs',filename(3:end)];
+% filename_ = ['/misc/HeLaZ_outputs',filename(3:end)];
+filename_ = ['/home/ahoffman/HeLaZ',filename(3:end)];
 CPUTIME   = double(h5readatt(filename_,'/data/input','cpu_time'));
 DT_SIM    = h5readatt(filename_,'/data/input','dt');
 
diff --git a/matlab/setup.m b/matlab/setup.m
index 6452cceaa4f74e73961ba41ebd24c2c06b65c478..2587142b70d8088d9d2ed4c290e4981f8f6c61d8 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -126,8 +126,6 @@ if W_NAPJ;   OUTPUTS.write_Napj  = '.true.'; else; OUTPUTS.write_Napj  = '.false
 if W_SAPJ;   OUTPUTS.write_Sapj  = '.true.'; else; OUTPUTS.write_Sapj  = '.false.';end;
 if W_DENS;   OUTPUTS.write_dens  = '.true.'; else; OUTPUTS.write_dens  = '.false.';end;
 if W_TEMP;   OUTPUTS.write_temp  = '.true.'; else; OUTPUTS.write_temp  = '.false.';end;
-OUTPUTS.resfile0    = '''outputs''';
-OUTPUTS.rstfile0    = '''checkpoint''';
 OUTPUTS.job2load    = JOB2LOAD;
 %% Create directories
 if ~exist(SIMDIR, 'dir')
@@ -143,7 +141,7 @@ end
 INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
 nproc = 1;
 MAKE  = 'cd ..; make; cd wk';
-system(MAKE);
+% system(MAKE);
 %%
 disp(['Set up ',SIMID]);
 disp([resolution,gridname,degngrad]);
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index b4bf4f0b1ab0de8c124e12fe3497dcb644c8b5ec..19b918c376ac2ca883b08cc613ca58e06814af43 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -40,8 +40,6 @@ fprintf(fid,['  write_Napj  = ', OUTPUTS.write_Napj,'\n']);
 fprintf(fid,['  write_Sapj  = ', OUTPUTS.write_Sapj,'\n']);
 fprintf(fid,['  write_dens  = ', OUTPUTS.write_dens,'\n']);
 fprintf(fid,['  write_temp  = ', OUTPUTS.write_temp,'\n']);
-fprintf(fid,['  resfile0    = ', OUTPUTS.resfile0,'\n']);
-fprintf(fid,['  rstfile0    = ', OUTPUTS.rstfile0,'\n']);
 fprintf(fid,['  job2load    = ', num2str(OUTPUTS.job2load),'\n']);
 fprintf(fid,'/\n');
 
diff --git a/matlab/write_sbash_marconi.m b/matlab/write_sbash_marconi.m
index 461b0c53e820cacf100cd8bff565e7b35c7c47ff..481b72548d9ad8c175cf22fa3b78372606813185 100644
--- a/matlab/write_sbash_marconi.m
+++ b/matlab/write_sbash_marconi.m
@@ -12,6 +12,7 @@ fprintf(fid,[...
 'cd ',BASIC.RESDIR,'\n',...
 'cp $HOME/HeLaZ/wk/fort*.90 .\n',...
 'cp $HOME/HeLaZ/wk/batch_script.sh .\n',...
+'rm $HOME/HeLaZ/wk/fort*.90\n',...
 ...
 SBATCH_CMD,...
 'echo tail -f $CINECA_SCRATCH/HeLaZ',BASIC.RESDIR(3:end),'out']);
@@ -31,8 +32,8 @@ fprintf(fid,[...
 '#SBATCH --cpus-per-task=', CLUSTER.CPUPT,'\n',...
 '#SBATCH --ntasks-per-node=', CLUSTER.NTPN,'\n',...
 '#SBATCH --mem=', CLUSTER.MEM,'\n',...
-'#SBATCH --error=err',num2str(JOB2LOAD+1),'.txt\n',...
-'#SBATCH --output=out_',num2str(JOB2LOAD+1),'.txt\n',...
+'#SBATCH --error=err_',sprintf('%2.2d',JOB2LOAD+1),'.txt\n',...
+'#SBATCH --output=out_',sprintf('%2.2d',JOB2LOAD+1),'.txt\n',...
 '#SBATCH --account=FUA35_TSVVT421\n',...
 '#SBATCH --partition=skl_fua_',CLUSTER.PART,'\n',...
 'module load autoload hdf5 fftw\n',...
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index dc2ba46b60a9e81c3b2e4a6b82ceaa884d46ab85..e661b401c309ccd1600d8f4f3a28a0a737f82172 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -5,24 +5,21 @@ outfile ='';
 if 1% Local results
 outfile ='';
 outfile ='';
-% outfile ='artificial_ZF_freeze/sim_A';
-% outfile ='simulation_B/cw_FCGK_kp_3.0';
-outfile ='nonlin_FCGK/150x75_L_200_P_4_J_2_eta_0.6_nu_1e-01_FCGK_mu_0e+00';
-% outfile ='nonlin_PAGK/100x50_L_200_P_4_J_2_eta_0.6_nu_1e-01_PAGK_mu_0e+00';
-% outfile ='nonlin_FCGK/100x50_L_200_P_4_J_2_eta_0.6_nu_1e-01_FCGK_mu_0e+00';
-% outfile ='simulation_A';
-% outfile ='simulation_B/cw_SGGK_like_species';
-% outfile ='simulation_A/CO_damping_SGGK';
-% outfile ='simulation_A/cw_DGGK_eta_0.5';
-% outfile ='simulation_B/cw_DGGK';
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='test_even_p/100x50_L_100_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_0e+00';
     BASIC.RESDIR      = ['../results/',outfile,'/'];
     BASIC.MISCDIR     = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
     system(['mkdir -p ',BASIC.MISCDIR]);
     CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
     system(CMD);
 else% Marconi results
-outfile ='/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';
-% 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';
+% outfile ='/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';
+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
@@ -77,35 +74,36 @@ if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
 t0    =000; iz = 1; ix = 1; iy = 1;
-skip_ =4; DELAY = 2e-3*skip_;
+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);
 T = Ts3D; FRAMES = FRAMES_3D;
-INTERP = 0; NORMALIZED = 0; CONST_CMAP = 0;% Gif options
+INTERP = 0; NORMALIZED = 0; CONST_CMAP = 0; BWR =1;% Gif options
 % Field to plot
-% FIELD = dens_e;       NAME = 'ne';   FIELDNAME = 'n_e';
 % FIELD = dens_i;       NAME = 'ni';   FIELDNAME = 'n_i';
-% FIELD = dens_e-Z_n_e; NAME = 'ne_NZ';FIELDNAME = 'n_e^{NZ}';
-FIELD = dens_i-Z_n_i; NAME = 'ni_NZ';FIELDNAME = 'n_i^{NZ}';
-% FIELD = temp_e;       NAME = 'Te';   FIELDNAME = 'n_i';
+% FIELD = dens_i-Z_n_i; NAME = 'ni_NZ';FIELDNAME = 'n_i^{NZ}';
 % FIELD = temp_i;       NAME = 'Ti';   FIELDNAME = 'n_i';
-% FIELD = temp_e-Z_T_e; NAME = 'Te_NZ';FIELDNAME = 'T_e^{NZ}';
 % 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 = Z_phi;    NAME = 'Zphi'; FIELDNAME = '\phi_Z';
-% FIELD = Gamma_x;  NAME = 'Gamma_x'; FIELDNAME = '\Gamma_x';
+% FIELD = ni00;         NAME = 'ni00'; FIELDNAME = 'n_i^{00}';
+FIELD = phi;          NAME = 'phi'; FIELDNAME = '\phi';
+% FIELD = Z_phi;        NAME = 'Zphi'; FIELDNAME = '\phi_Z';
+% FIELD = Gamma_x;      NAME = 'Gamma_x'; FIELDNAME = '\Gamma_x';
 
 % Sliced
 % plt = @(x) real(x(ix, :, :,:)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; 
 % plt = @(x) real(x( :,iy, :,:)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; 
-% plt = @(x) real(x( :, :,iz,:)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; 
+plt = @(x) real(x( :, :,iz,:)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; 
+
+% % K-space
+% % FIELD = PHI;          NAME = 'PHI'; FIELDNAME = '\tilde \phi';
+% FIELD = Ni00;         NAME = 'Ni00'; FIELDNAME = 'N_i^{00}';
+% plt = @(x) fftshift((abs(x( :, :,1,:))),2); X = fftshift(KX,2); Y = fftshift(KY,2); XNAME = 'k_x'; YNAME = 'k_y'; 
 
 % Averaged
 % plt = @(x) mean(x,1); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; 
 % plt = @(x) mean(x,2); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; 
-plt = @(x) mean(x,3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; 
+% plt = @(x) mean(x,3); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; 
 
 FIELD = squeeze(plt(FIELD));
 
@@ -233,7 +231,7 @@ end
 
 if 1
 %% zonal vs nonzonal energies for phi(t)
-trange = 200:Ns3D;
+trange = 1:Ns3D;
 Ephi_Z           = zeros(1,Ns3D);
 Ephi_NZ_kgt0      = zeros(1,Ns3D);    
 Ephi_NZ_kgt1      = zeros(1,Ns3D);    
diff --git a/wk/continue_multiple_runs_marconi.m b/wk/continue_multiple_runs_marconi.m
index 944063f38c24f180a09d8b755907689bda0dff84..ea1fa75c00f0a242920a366cc19587bda10be206 100644
--- a/wk/continue_multiple_runs_marconi.m
+++ b/wk/continue_multiple_runs_marconi.m
@@ -4,7 +4,8 @@ continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/
 
 %% Functions to modify preexisting fort.90 input file and launch on marconi
 function [] = continue_run(outfilename)
-    EXECNAME = 'helaz_3.2';
+    EXECNAME = 'helaz_3.81';
+    SBATCH_CMD = 'sbatch batch_script.sh\n';
     %% CLUSTER PARAMETERS
     CLUSTER.PART  = 'prod';     % dbg or prod
     CLUSTER.TIME  = '24:00:00'; % allocation time hh:mm:ss
@@ -45,7 +46,7 @@ function [] = continue_run(outfilename)
         line = A{39};
         line = line(end-2:end);
         if(line(1) == '='); line = line(end); end;
-        J2L = str2num(line)+1;
+        J2L = str2num(line)+1; JOB2LOAD=J2L;
     end
     % Change job 2 load in fort.90
     A{39} = ['  job2load      = ',num2str(J2L)];
diff --git a/wk/linear_study.m b/wk/linear_study.m
index 547bb0c58d65809259c17223bae679ffbe415790..eed525c5b6e8ed10e735cf5aafa17f3598a32ddd 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -1,4 +1,4 @@
-for CO = [1]
+for CO = [3]
     RUN = 1; % To run or just to load
 addpath(genpath('../matlab')) % ... add
 default_plots_options
@@ -10,7 +10,7 @@ CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 NU      = 0.1;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
 ETAB    = 1.0;
-ETAN    = 1/0.6;    % Density gradient
+ETAN    = 1./0.6;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
 NU_HYP  = 0.0;   % Hyperdiffusivity coefficient
 LAMBDAD = 0.0;
@@ -23,45 +23,38 @@ KXEQ0   = 1;      % put kx = 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    = 400;  % Maximal time unit
+TMAX    = 200;  % Maximal time unit
 DT      = 1e-2;   % Time step
 SPS0D   = 1;      % Sampling per time unit for 2D arrays
 SPS2D   = 0;      % Sampling per time unit for 2D arrays
 SPS3D   = 2;      % Sampling per time unit for 2D arrays
 SPS5D   = 2;    % Sampling per time unit for 5D arrays
 SPSCP   = 0;    % Sampling per time unit for checkpoints
-RESTART = 0;      % To restart from last checkpoint
-JOB2LOAD= 00;
+JOB2LOAD= -1;
 %% OPTIONS
 % SIMID   = 'v3.6_kobayashi_lin';  % Name of the simulation
 % SIMID   = 'v3.2_CO_damping';  % Name of the simulation
 % SIMID   = 'CO_Patchwork_damping';  % Name of the simulation
-SIMID   = 'test_GF_closure';  % Name of the simulation
+SIMID   = 'test_even_p';  % Name of the simulation
 % SIMID   = 'v3.2_entropy_mode_linear';  % Name of the simulation
 NON_LIN = 0 *(1-KXEQ0);   % activate non-linearity (is cancelled if KXEQ0 = 1)
 % Collision operator
 % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle, 4 : Full Couloumb ; +/- for GK/DK)
 % CO      = 1;
 INIT_ZF = 0; ZF_AMP = 0.0;
-CLOS    = 1;   % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
+CLOS    = 0;   % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
 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
 %% OUTPUTS
-W_DOUBLE = 0;
-W_GAMMA  = 0;
-W_PHI    = 1;
-W_NA00   = 1;
-W_NAPJ   = 1;
-W_SAPJ   = 0;
-W_DENS   = 0;
-W_TEMP   = 0;
+W_DOUBLE = 1;
+W_GAMMA  = 1; W_HF     = 1;
+W_PHI    = 1; W_NA00   = 1;
+W_DENS   = 1; W_TEMP   = 1;
+W_NAPJ   = 1; W_SAPJ   = 0;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % unused
-% DK    = 0;  % Drift kinetic model (put every kernel_n to 0 except n=0 to 1)
-JOBNUM  = 00;
-KPAR    = 0.0;    % Parellel wave vector component
 HD_CO   = 0.5;    % Hyper diffusivity cutoff ratio
 kmax    = N*pi/L;% Highest fourier mode
 MU      = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
@@ -76,8 +69,8 @@ if 1
 %% Parameter scan over PJ
 % PA = [2 4];
 % JA = [1 2];
-PA = [5];
-JA = [2];
+PA = [8];
+JA = [4];
 DTA= DT*ones(size(JA));%./sqrt(JA);
 % DTA= DT;
 mup_ = MU_P;
@@ -95,10 +88,12 @@ for i = 1:Nparam
     JMAXE = JA(i); JMAXI = JA(i);
     DT = DTA(i);
     setup
+    system(['rm fort*.90']);
     % Run linear simulation
     if RUN
-        system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6; cd ../../../wk'])
-%         system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz; cd ../../../wk'])
+        system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6 0; cd ../../../wk'])
+%         system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3.82 1 6 0; cd ../../../wk'])
+%         system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
     end
 %     Load and process results
     %%
diff --git a/wk/local_run.m b/wk/local_run.m
index 06bcd4039f5c60008d0db15551403c539f7957ef..3a3472529a42e1555a946eadd641983c3c36a8e1 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -8,8 +8,8 @@ NU      = 0.1;   % Collision frequency
 ETAN    = 1.0/0.6;    % Density gradient drive (R/Ln)
 NU_HYP  = 0.0;
 %% GRID AND GEOMETRY PARAMETERS
-N       = 150;     % Frequency gridpoints (Nkx = N/2)
-L       = 200;     % Size of the squared frequency domain
+N       = 100;     % Frequency gridpoints (Nkx = N/2)
+L       = 100;     % Size of the squared frequency domain
 Nz      = 1;      % number of perpendicular planes (parallel grid)
 q0      = 1.0;    % safety factor
 shear   = 0.0;    % magnetic shear
@@ -32,7 +32,7 @@ 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   = 'nonlin_FCGK';  % Name of the simulation
-SIMID   = 'test';  % Name of the simulation
+SIMID   = 'test_even_p';  % 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
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 9c2b2b67b9a80bfbc8167e44cf5a29d8a9abec47..64959f91fad0e6143610d831fa34d532c8d16e10 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -15,7 +15,7 @@ CLUSTER.TIME  = '20:00:00'; % allocation time hh:mm:ss
 if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME  = '00:30:00'; end;
 CLUSTER.MEM   = '128GB';     % Memory
 CLUSTER.JNAME = 'HeLaZ';% Job name
-NP_P          = 2;          % MPI processes along p
+NP_P          = 1;          % MPI processes along p
 NP_KX         = 48;         % MPI processes along kx
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
@@ -23,27 +23,27 @@ NU      = 0.1;   % Collision frequency
 ETAN    = 1.0/0.6;    % Density gradient drive (R/Ln)
 NU_HYP  = 0.0;
 %% GRID PARAMETERS
-N       = 300;     % Frequency gridpoints (Nkx = N/2)
+N       = 500;     % 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       = 8;
-J       = 4;
+P       = 4;
+J       = 2;
 %% TIME PARAMETERS
 TMAX    = 10000;  % Maximal time unit
-DT      = 8e-3;   % Time step
+DT      = 1e-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/100;  % Sampling per time unit for 5D arrays
 SPSCP   = 0;    % Sampling per time unit for checkpoints/10
-JOB2LOAD= -1; % start from t=0 if <0, else restart from outputs_$job2load
+JOB2LOAD= 0; % start from t=0 if <0, else restart from outputs_$job2load
 %% OPTIONS AND NAMING
 % Collision operator
 % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK)
-CO      = 3;
+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   = 'test_chained_job';  % Name of the simulation
@@ -94,7 +94,6 @@ CLUSTER.CPUPT = '1';        % CPU per task
 setup
 SBATCH_CMD = 'sbatch batch_script.sh\n';
 write_sbash_marconi
-system('rm fort*.90 setup_and_run.sh batch_script.sh');
 if(mod(NP_P*NP_KX,48)~= 0)
     disp('WARNING : unused cores (ntot cores must be a 48 multiple)');
 end
@@ -114,5 +113,6 @@ if(SUBMIT)
         end
     end
 end
+system('rm fort*.90');
 disp('done');
 end