diff --git a/wk/KBM_ky_PJ_scan.m b/wk/KBM_ky_PJ_scan.m
deleted file mode 100644
index 45aff676d44927d3799d9689154b55756a104b68..0000000000000000000000000000000000000000
--- a/wk/KBM_ky_PJ_scan.m
+++ /dev/null
@@ -1,218 +0,0 @@
-gyacomodir  = pwd;
-gyacomodir = gyacomodir(1:end-2);
-addpath(genpath([gyacomodir,'matlab'])) % ... add
-addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
-addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
-addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
-EXECNAME = 'gyacomo23_sp';
-% EXECNAME = 'gyacomo23_dp';
-% EXECNAME = 'gyacomo23_debug';
-CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
-%%
-SIMID = 'lin_KBM';  % Name of the simulation
-RERUN   = 0; % rerun if the  data does not exist
-RUN     = 1;
-K_Ne    = 3;                % ele Density '''
-K_Te    = 4.5;              % ele Temperature '''
-K_Ni    = 3;                % ion Density gradient drive
-K_Ti    = 8;                % ion Temperature '''
-P_a = [2 4 8 16];
-% P_a = 10;
-ky_a= 0.05:0.1:0.75;
-% ky_a = 0.05;
-% collision setting
-CO        = 'DG';
-GKCO      = 1; % gyrokinetic operator
-COLL_KCUT = 1.75;
-NU        = 1e-2;
-% model
-NA      = 2;
-BETA    = 0.03;     % electron plasma beta
-% Geometry
-GEOMETRY= 'miller';
-% GEOMETRY= 's-alpha';
-SHEAR   = 0.8;    % magnetic shear
-% time and numerical grid
-DT0    = 5e-3;
-TMAX   = 15;
-% arrays for the result
-g_ky = zeros(numel(ky_a),numel(P_a),2);
-g_avg= g_ky*0;
-g_std= g_ky*0;
-% Naming of the collision operator
-if GKCO
-    CONAME = [CO,'GK'];
-else
-    CONAME = [CO,'DK'];
-end
-j = 1;
-for P = P_a
-    J = P/2;
-i = 1;
-for ky = ky_a
-    %% PHYSICAL PARAMETERS
-    TAU     = 1.0;            % e/i temperature ratio
-    % SIGMA_E = 0.05196152422706632;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-    SIGMA_E = 0.0233380;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-    %% GRID PARAMETERS
-    DT = DT0/sqrt(J);
-    PMAX   = P;     % Hermite basis size
-    JMAX   = P/2;     % Laguerre "
-    NX      = 12;    % real space x-gridpoints
-    NY      = 2;
-    LX      = 2*pi/0.8;   % Size of the squared frequency domain
-    LY      = 2*pi/ky;     % Size of the squared frequency domain
-    NZ      = 24;    % number of perpendicular planes (parallel grid)
-    NPOL    = 1;
-    SG      = 0;     % Staggered z grids option
-    NEXC    = 1;     % To extend Lx if needed (Lx = Nexc/(kymin*shear))
-    %% GEOMETRY
-    % GEOMETRY= 's-alpha';
-    EPS     = 0.18;   % inverse aspect ratio
-    Q0      = 1.4;    % safety factor
-    KAPPA   = 1.0;    % elongation
-    DELTA   = 0.0;    % triangularity
-    ZETA    = 0.0;    % squareness
-    PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
-%     PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
-    SHIFT_Y = 0.0;
-    %% TIME PARMETERS
-    DTSAVE0D = 1;      % Sampling per time unit for 0D arrays
-    DTSAVE2D = -1;     % Sampling per time unit for 2D arrays
-    DTSAVE3D = 0.5;      % Sampling per time unit for 3D arrays
-    DTSAVE5D = 100;     % Sampling per time unit for 5D arrays
-    JOB2LOAD = -1;     % Start a new simulation serie
-    %% OPTIONS
-    LINEARITY = 'linear';   % activate non-linearity (is cancelled if KXEQ0 = 1)
-    % Collision operator
-    ABCO    = 1; % INTERSPECIES collisions
-    INIT_ZF = 0; ZF_AMP = 0.0;
-    CLOS    = 0;   % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
-    NL_CLOS = 0;   % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
-    KERN    = 0;   % Kernel model (0 : GK)
-    INIT_OPT= 'phi';   % Start simulation with a noisy mom00/phi/allmom
-    NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
-    %% OUTPUTS
-    W_DOUBLE = 0;
-    W_GAMMA  = 1; W_HF     = 1;
-    W_PHI    = 1; W_NA00   = 1;
-    W_DENS   = 0; W_TEMP   = 1;
-    W_NAPJ   = 0; W_SAPJ   = 0;
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    % unused
-    ADIAB_E = (NA==1);
-    HD_CO   = 0.0;    % Hyper diffusivity cutoff ratio
-    MU      = 0.0; % Hyperdiffusivity coefficient
-    INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
-    MU_X    = MU;     %
-    MU_Y    = MU;     %
-    N_HD    = 4;
-    HYP_V   = 'none';
-    HRCY_CLOS = 'truncation';   % Closure model for higher order moments
-    DMAX      = -1;
-    NLIN_CLOS = 'truncation';   % Nonlinear closure model for higher order moments
-    NMAX      = 0;
-    MU_Z    = 1.0;     %
-    MU_P    = 0.0;     %
-    MU_J    = 0.0;     %
-    LAMBDAD = 0.0;
-    NOISE0  = 1.0e-5; % Init noise amplitude
-    BCKGD0  = 0.0;    % Init background
-    k_gB   = 1.0;
-    k_cB   = 1.0;
-    %% RUN
-    setup
-    % naming
-    filename = [SIMID,'/',PARAMS,'/'];
-    LOCALDIR  = [gyacomodir,'results/',filename,'/'];
-    % check if data exist to run if no data
-    data_ = {};
-    try
-        data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
-        Ntime = numel(data_.Ts0D);
-    catch
-        data_.outfilenames = [];
-    end
-    if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10)
-        % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
-        system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 2 2 1 0; cd ../../../wk'])
-        % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
-    end
-    data_    = compile_results_low_mem(data_,LOCALDIR,00,00);
-    [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
-    if numel(data_.Ts3D)>10
-        if numel(data_.Ts3D)>5
-        % Load results after trying to run
-        filename = [SIMID,'/',PARAMS,'/'];
-        LOCALDIR  = [gyacomodir,'results/',filename,'/'];
-
-        data_    = compile_results_low_mem(data_,LOCALDIR,00,00);
-        [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
-
-        % linear growth rate (adapted for 2D zpinch and fluxtube)
-        options.TRANGE = [0.5 1]*data_.Ts3D(end);
-        options.NPLOTS = 0; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
-        options.GOK    = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
-
-        [~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
-        [~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
-        field   = 0;
-        field_t = 0;
-        for ik = 2:NY/2+1
-            field   = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
-            field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
-            to_measure  = log(field_t(it1:it2));
-            tw = double(data_.Ts3D(it1:it2));
-    %         gr = polyfit(tw,to_measure,1);
-            gr = fit(tw,to_measure,'poly1');
-            err= confint(gr);
-            g_ky(i,j,ik)  = gr.p1;
-            g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
-        end
-        [gmax, ikmax] = max(g_ky(i,j,:));
-
-        msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
-        end
-    end
-    
-    i = i + 1;
-end
-j = j + 1;
-end
-
-if 0
-%% Check time evolution
-figure;
-to_measure  = log(field_t);
-plot(data_.Ts3D,to_measure); hold on
-plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
-end
-
-%% take max growth rate among z coordinate
-y_ = g_ky(:,:,2); 
-e_ = g_std(:,:,2);
-
-%%
-if(numel(ky_a)>1 && numel(P_a)>1)
-%% Save metadata
- pmin  = num2str(min(P_a));   pmax = num2str(max(P_a));
- kymin = num2str(min(ky_a));  kymax= num2str(max(ky_a));
-filename = [num2str(NX),'x',num2str(NZ),'_ky_',kymin,'_',kymax,...
-            '_P_',pmin,'_',pmax,'_',CONAME,'_',num2str(NU),'_be_',num2str(BETA),'.mat'];
-metadata.name   = filename;
-metadata.kymin  = ky;
-metadata.title  = ['$\nu_{',CONAME,'}=$',num2str(NU),'$\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N)];
-metadata.par    = [num2str(NX),'x1x',num2str(NZ)];
-metadata.nscan  = 2;
-metadata.s2name = '$P$';
-metadata.s2     = P_a;
-metadata.s1name = '$ky$';
-metadata.s1     = ky_a;
-metadata.dname  = '$\gamma c_s/R$';
-metadata.data   = y_;
-metadata.err    = e_;
-save([SIMDIR,filename],'-struct','metadata');
-disp(['saved in ',SIMDIR,filename]);
-clear metadata tosave
-end
diff --git a/wk/KBM_load_metadata_scan.m b/wk/KBM_load_metadata_scan.m
deleted file mode 100644
index b86161b30dce5a2e4937409ecc781439b1afbecc..0000000000000000000000000000000000000000
--- a/wk/KBM_load_metadata_scan.m
+++ /dev/null
@@ -1,172 +0,0 @@
-% Metadata path
-gyacomodir  = pwd;
-gyacomodir = gyacomodir(1:end-2);
-addpath(genpath([gyacomodir,'matlab'])) % ... add
-addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
-addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
-addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
-
-% datafname = 'lin_KBM/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.03.mat';
-% datafname = 'lin_KBM/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.03.mat';
-datafname ='lin_AUG_scan/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.000152.mat';
-%% Chose if we filter gamma>0.05
-FILTERGAMMA = 0;
-
-%% Load data
-fname = ['../results/',datafname];
-d = load(fname);
-if FILTERGAMMA
-    d.data = d.data.*(d.data>0.025);
-    d.err  = d.err.*(d.data>0.025);
-end
-if 0
-%% Pcolor of the peak
-figure;
-% [XX_,YY_] = meshgrid(d.s1,d.s2);
-[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
-pclr=imagesc_custom(XX_,YY_,d.data'.*(d.data>0)');
-% pclr=contourf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
-% pclr=surf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
-title(d.title);
-xlabel(d.s1name); ylabel(d.s2name);
-set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
-set(gca,'YTick',1:numel(d.s2),'YTicklabel',d.s2)
-colormap(jet)
-colormap(bluewhitered)
-clb=colorbar; 
-clb.Label.String = '$\gamma c_s/R$';
-clb.Label.Interpreter = 'latex';
-clb.Label.FontSize= 18;
-end
-if 1
-%% Scan along first dimension
-figure
-colors_ = jet(numel(d.s2));
-for i = 1:numel(d.s2)
-    % plot(d.s1,d.data(:,i),'s-',...
-    plot(d.s1(d.data(:,i)>0),d.data((d.data(:,i)>0),i),'s-',...
-        'LineWidth',2.0,...
-        'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
-        'color',colors_(i,:)); 
-    % errorbar(d.s1,d.data(:,i),d.err(:,i),'s-',...
-    % 'LineWidth',2.0,...
-    % 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
-    % 'color',colors_(i,:)); 
-    hold on;
-end
-xlabel(d.s1name); ylabel(d.dname);title(d.title);
-xlim([d.s1(1) d.s1(end)]);
-colormap(colors_);
-clb = colorbar;
-% caxis([d.s2(1)-0.5,d.s2(end)+0.5]);
-clim([1 numel(d.s2)+1]);
-clb.Ticks=linspace(d.s2(1),d.s2(end),numel(d.s2));
-clb.Ticks    =1.5:numel(d.s2)+1.5;
-clb.TickLabels=d.s2;
-clb.Label.String = d.s2name;
-clb.Label.Interpreter = 'latex';
-clb.Label.FontSize= 18;
-end
-if 0
-%% Scan along second dimension
-figure
-colors_ = jet(numel(d.s1));
-for i = 1:numel(d.s1)
-    plot(d.s2,d.data(i,:),'s-',...
-        'LineWidth',2.0,...
-        'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
-        'color',colors_(i,:)); 
-%     errorbar(d.s2,d.data(i,:),d.err(i,:),'s-',...
-%         'LineWidth',2.0,...
-%         'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
-%         'color',colors_(i,:)); 
-    hold on;
-end
-xlabel(d.s2name); ylabel(d.dname);title(d.title);
-xlim([d.s2(1) d.s2(end)]);
-colormap(jet(numel(d.s1)));
-clb = colorbar;
-caxis([d.s1(1)-0.5,d.s1(end)+0.5]);
-clb.Ticks=linspace(d.s1(1),d.s1(end),numel(d.s1));
-clb.YTick=d.s1;
-clb.Label.String = d.s1name;
-clb.TickLabelInterpreter = 'latex';
-clb.Label.Interpreter = 'latex';
-clb.Label.FontSize= 18;
-end
-
-if 0
-%% Convergence analysis
-figure
-% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
-colors_ = jet(numel(d.s1));
-for i = 1:numel(d.s1)
-%     target_ = d.data(i,end);
-    target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
-    if target_ > 0
-    eps_    = abs(target_ - d.data(i,1:end-1))/abs(target_);
-    semilogy(d.s2(1:end-1),eps_,'-s',...
-        'LineWidth',2.0,...
-        'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
-        'color',colors_(i,:));
-    hold on;
-    end
-end
-xlabel(d.s2name); ylabel('$\epsilon_r$');title(d.title);
-xlim([d.s2(1) d.s2(end)]);
-colormap(colors_);
-clb = colorbar;
-caxis([d.s1(1)-0.5,d.s1(end)+0.5]);
-clb.Ticks=linspace(d.s1(1),d.s1(end),numel(d.s1));
-clb.YTick=d.s1;
-clb.Label.String = d.s1name;
-clb.TickLabelInterpreter = 'latex';
-clb.Label.Interpreter = 'latex';
-clb.Label.FontSize= 18;
-grid on;
-end
-if 0
-%% Pcolor of the error
-figure;  i_ = 0;
-% target_ = 2.72724991618068013377e-01; % Value for nuDGDK = 1.0, kT=6.96, (40,20), Nkx=8
-% target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
-% target_ = 2.71455e-01; % Value for nuDGDK = 0.001, kT=6.96, (40,20), Nkx=8
-% target_ = 1.39427e-01; % Value for nuDGDK = 0.001, kT=5.3, (30,16), Nkx=8
-target_ = 2.72510405826983714839e-01  % Value for nuDGDK = 0.001, kT=6.96, (50,25), Nkx=8
-% target_ = 2.73048910051283844069e-01; % Value for nuDGDK = 0.0, kT=6.96, (40,20), Nkx=8
-% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
-% eps_    = log(abs(target_ - d.data)/abs(target_));
-eps_    = max(-10,log(abs(target_ - d.data)/abs(target_)));
-sign_   = 1;%sign(d.data - target_);
-eps_ = d.data;
-for i = 1:numel(d.s1)
-    for j = 1:numel(d.s2)
-        % target_ = d.data(i,end);
-        % target_ = d.data(i,end);
-        % target_ = d.data(end,j);
-        eps_(i,j) = log(abs(target_ - d.data(i,j))/target);
-        if target_ > 0
-    %     eps_(i,:) = max(-12,log(abs(target_ - d.data(i,1:end))/abs(target_)));
-    %     eps_(i,:) = log(abs(target_ - d.data(i,1:end))/abs(target_));
-    %     eps_(i,:) = min(100,100*abs(target_ - d.data(i,1:end))/abs(target_));
-        else
-        end
-    end
-end
-[XX_,YY_] = meshgrid(d.s1,d.s2);
-[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
-pclr=imagesc_custom(XX_,YY_,eps_'.*(d.data>0)'.*sign_');
-% pclr=contourf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_',5);
-% pclr=surf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_');
-title(d.title);
-xlabel(d.s1name); ylabel(d.s2name);
-set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
-set(gca,'YTick',1:numel(d.s2),'YTicklabel',d.s2)
-% colormap(jet)
-colormap(bluewhitered)
-% caxis([-10, 0]);
-clb=colorbar; 
-clb.Label.String = '$\log(\epsilon_r)$';
-clb.Label.Interpreter = 'latex';
-clb.Label.FontSize= 18;
-end
\ No newline at end of file
diff --git a/wk/lin_AUG.m b/wk/lin_AUG.m
deleted file mode 100644
index 8cfded1e89c2b8bb6f06fbc5abfb222e00f48ff3..0000000000000000000000000000000000000000
--- a/wk/lin_AUG.m
+++ /dev/null
@@ -1,183 +0,0 @@
-%% QUICK RUN SCRIPT
-% This script creates a directory in /results and runs a simulation directly
-% from the Matlab framework. It is meant to run only small problems in linear
-% for benchmarking and debugging purposes since it makes Matlab "busy".
-
-%% Set up the paths for the necessary Matlab modules
-gyacomodir = pwd;
-gyacomodir = gyacomodir(1:end-2);
-addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
-addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
-addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
-addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
-
-%% Set simulation parameters
-SIMID   = 'lin_KBM';  % Name of the simulation
-RUN = 1; % To run or just to load
-default_plots_options
-% EXECNAME = 'gyacomo23_sp'; % single precision
-EXECNAME = 'gyacomo23_dp'; % double precision
-% EXECNAME = 'gyacomo23_debug'; % single precision
-
-%% Set up physical parameters
-CLUSTER.TIME = '99:00:00';  % Allocation time hh:mm:ss
-NU = 0.1;                   % Collision frequency
-TAU = 1.38;                  % e/i temperature ratio
-K_Ne    = 34;              % ele Density '''
-K_Te    = 56;               % ele Temperature '''
-K_Ni    = 34;              % ion Density gradient drive
-K_Ti    = 21;               % ion Temperature '''
-SIGMA_E = 0.0233380;        % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-NA = 2;                     % number of kinetic species
-ADIAB_E = (NA==1);          % adiabatic electron model
-BETA    = 1.52e-4;             % electron plasma beta
-MHD_PD  = 0;                % MHD pressure drift
-NPOL    = 1;                % Number of poloidal turns
-%% Set up grid parameters
-P = 4;
-J = P/2;%P/2;
-PMAX = P;                   % Hermite basis size
-JMAX = J;                   % Laguerre basis size
-NX = 6;                     % real space x-gridpoints
-NY = 8;                    % real space y-gridpoints
-LX = 2*pi/0.3;              % Size of the squared frequency domain in x direction
-LY = 2*pi/0.4;              % Size of the squared frequency domain in y direction
-NZ = 32*NPOL;                    % number of perpendicular planes (parallel grid)
-SG = 0;                     % Staggered z grids option
-NEXC = 1;                   % To extend Lx if needed (Lx = Nexc/(kymin*shear))
-%% GEOMETRY
-% GEOMETRY= 's-alpha';
-GEOMETRY= 'miller';
-% GEOMETRY= 'circular';
-EPS     = 0.3;   % inverse aspect ratio
-Q0      = 5.7;    % safety factor
-SHEAR   = 4.6;%*EPS;    % magnetic shear
-KAPPA   = 1.8;   % elongation
-S_KAPPA = 0.0;
-DELTA   = 0.4;  % triangularity
-S_DELTA = 0.0;
-ZETA    = 0.0; % squareness
-S_ZETA  = 0.0;
-PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
-SHIFT_Y = 0.0;    % Shift in the periodic BC in z
-PB_PHASE= 0;
-%% TIME PARAMETERS
-TMAX     = 5;  % Maximal time unit
-DT       = 2e-3;   % Time step
-DTSAVE0D = 1;      % Sampling per time unit for 0D arrays
-DTSAVE2D = -1;     % Sampling per time unit for 2D arrays
-DTSAVE3D = 3;      % Sampling per time unit for 3D arrays
-DTSAVE5D = 100;     % Sampling per time unit for 5D arrays
-JOB2LOAD = -1;     % Start a new simulation serie
-
-%% OPTIONS
-LINEARITY = 'linear';   % activate non-linearity (is cancelled if KXEQ0 = 1)
-CO        = 'DG';       % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-GKCO      = 1;          % Gyrokinetic operator
-ABCO      = 1;          % INTERSPECIES collisions
-INIT_ZF   = 0;          % Initialize zero-field quantities
-HRCY_CLOS = 'truncation';   % Closure model for higher order moments
-DMAX      = -1;
-NLIN_CLOS = 'truncation';   % Nonlinear closure model for higher order moments
-NMAX      = 0;
-KERN      = 0;   % Kernel model (0 : GK)
-INIT_OPT  = 'mom00';   % Start simulation with a noisy mom00/phi/allmom
-NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
-
-%% OUTPUTS
-W_DOUBLE = 1;     % Output flag for double moments
-W_GAMMA  = 1;     % Output flag for gamma (Gyrokinetic Energy)
-W_HF     = 1;     % Output flag for high-frequency potential energy
-W_PHI    = 1;     % Output flag for potential
-W_NA00   = 1;     % Output flag for nalpha00 (density of species alpha)
-W_DENS   = 1;     % Output flag for total density
-W_TEMP   = 1;     % Output flag for temperature
-W_NAPJ   = 1;     % Output flag for nalphaparallel (parallel momentum of species alpha)
-W_SAPJ   = 0;     % Output flag for saparallel (parallel current of species alpha)
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% UNUSED PARAMETERS
-% These parameters are usually not to play with in linear runs
-MU      = 0.1;    % Hyperdiffusivity coefficient
-MU_X    = MU;     % Hyperdiffusivity coefficient in x direction
-MU_Y    = MU;     % Hyperdiffusivity coefficient in y direction
-N_HD    = 4;      % Degree of spatial-hyperdiffusivity
-MU_Z    = 10.0;    % Hyperdiffusivity coefficient in z direction
-HYP_V   = 'hypcoll'; % Kinetic-hyperdiffusivity model
-MU_P    = 0.0;    % Hyperdiffusivity coefficient for Hermite
-MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
-LAMBDAD = 0.0;    % Lambda Debye
-NOISE0  = 0.0e-5; % Initial noise amplitude
-BCKGD0  = 1.0e-5;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
-COLL_KCUT = 1; % Cutoff for collision operator
-ADIAB_I = 0;          % adiabatic ion model
-
-%%-------------------------------------------------------------------------
-%% RUN
-setup
-% system(['rm fort*.90']);
-% 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 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
-%     RUN  =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
-    % RUN  =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
-    MVOUT='cd ../../../wk;';
-    system([MVIN,RUN,MVOUT]);
-end
-
-%% Analysis
-% load
-filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
-LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
-FIGDIR   = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
-% Load outputs from jobnummin up to jobnummax
-J0 = 0; J1 = 0;
-data = {}; % Initialize data as an empty cell array
-% load grids, inputs, and time traces
-data = compile_results_low_mem(data,LOCALDIR,J0,J1); 
-
-
-if 0 % Activate or not
-%% plot mode evolution and growth rates
-% Load phi
-[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
-options.NORMALIZED = 0; 
-options.TIME   = data.Ts3D;
- % Time window to measure the growth of kx/ky modes
-options.KX_TW  = [0.2 1]*data.Ts3D(end);
-options.KY_TW  = [0.2 1]*data.Ts3D(end);
-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.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
-end
-
-if 1
-%% Ballooning plot
-[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
-if data.inputs.BETA > 0
-[data.PSI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'psi');
-end
-options.time_2_plot = [120];
-options.kymodes     = [0.1:0.1:1];
-options.normalized  = 1;
-options.PLOT_KP     = 0;
-% options.field       = 'phi';
-fig = plot_ballooning(data,options);
-title(GEOMETRY)
-end
-
-if 0
-%% plot geom and metric coefficients
-options.SHOW_FLUXSURF = 1;
-options.SHOW_METRICS  = 1;
-[fig, geo_arrays] = plot_metric(data,options);
-title(GEOMETRY)
-end
\ No newline at end of file
diff --git a/wk/lin_DTT_AB_rho85_PT.m b/wk/parameters/lin_DTT_AB_rho85_PT.m
similarity index 57%
rename from wk/lin_DTT_AB_rho85_PT.m
rename to wk/parameters/lin_DTT_AB_rho85_PT.m
index d78648aa85ab16518dd607609f21ab23bccdc820..4fa60133f5ffb1811c686a3f74ed2fc6900cee96 100644
--- a/wk/lin_DTT_AB_rho85_PT.m
+++ b/wk/parameters/lin_DTT_AB_rho85_PT.m
@@ -1,25 +1,5 @@
-%% QUICK RUN SCRIPT
-% This script creates a directory in /results and runs a simulation directly
-% from the Matlab framework. It is meant to run only small problems in linear
-% for benchmarking and debugging purposes since it makes Matlab "busy".
-
-%% Set up the paths for the necessary Matlab modules
-gyacomodir = pwd;
-gyacomodir = gyacomodir(1:end-2);
-% mpirun     = 'mpirun';
-mpirun     = '/opt/homebrew/bin/mpirun';
-addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
-addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
-addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
-addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
-
 %% Set simulation parameters
-SIMID   = 'lin_DTT';  % Name of the simulation
-RUN = 1; % To run or just to load
-default_plots_options
-EXECNAME = 'gyacomo23_sp'; % single precision
-% EXECNAME = 'gyacomo23_dp'; % double precision
-
+SIMID   = 'lin_DTT_AB_rho85_PT';  % Name of the simulation
 %% Set up physical parameters
 CLUSTER.TIME = '99:00:00';  % Allocation time hh:mm:ss
 NU = 0.05;                  % Collision frequency
@@ -38,11 +18,11 @@ P = 4;
 J = P/2;%P/2;
 PMAX = P;                   % Hermite basis size
 JMAX = J;                   % Laguerre basis size
-NX = 32;                    % real space x-gridpoints
+NX = 16;                    % real space x-gridpoints
 NY = 16;                     % real space y-gridpoints
 LX = 2*pi/0.1;              % Size of the squared frequency domain in x direction
 LY = 2*pi/0.2;             % Size of the squared frequency domain in y direction
-NZ = 32;                    % number of perpendicular planes (parallel grid)
+NZ = 24;                    % number of perpendicular planes (parallel grid)
 SG = 0;                     % Staggered z grids option
 NEXC = 1;                   % To extend Lx if needed (Lx = Nexc/(kymin*shear))
 
@@ -114,64 +94,4 @@ BCKGD0  = 0.0e-5;    % Initial background
 k_gB   = 1.0;     % Magnetic gradient strength
 k_cB   = 1.0;     % Magnetic curvature strength
 COLL_KCUT = 1; % Cutoff for collision operator
-ADIAB_I = 0;          % adiabatic ion model
-
-%%-------------------------------------------------------------------------
-%% RUN
-setup
-% system(['rm fort*.90']);
-% 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 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
-     RUN  =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 1 4 2 0;'];
-%     RUN  =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
-      % RUN = ['./../../../bin/gyacomo23_sp 0;'];
-    MVOUT='cd ../../../wk;';
-    system([MVIN,RUN,MVOUT]);
-end
-
-%% Analysis
-% load
-filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
-LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
-FIGDIR   = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
-% Load outputs from jobnummin up to jobnummax
-J0 = 0; J1 = 0;
-data = {}; % Initialize data as an empty cell array
-% load grids, inputs, and time traces
-data = compile_results_low_mem(data,LOCALDIR,J0,J1); 
-
-
-if 0 % Activate or not
-%% plot mode evolution and growth rates
-% Load phi
-[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
-options.NORMALIZED = 0; 
-options.TIME   = data.Ts3D;
- % Time window to measure the growth of kx/ky modes
-options.KX_TW  = [0.2 1]*data.Ts3D(end);
-options.KY_TW  = [0.2 1]*data.Ts3D(end);
-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.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
-end
-
-if 1
-%% Ballooning plot
-[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
-if data.inputs.BETA > 0
-[data.PSI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'psi');
-end
-options.time_2_plot = [120];
-options.kymodes     = [0.25];
-options.normalized  = 1;
-options.PLOT_KP     = 0;
-% options.field       = 'phi';
-fig = plot_ballooning(data,options);
-end
-
+ADIAB_I = 0;          % adiabatic ion model
\ No newline at end of file
diff --git a/wk/lin_ETG_adiab_i.m b/wk/parameters/lin_ETG_adiab_i.m
similarity index 63%
rename from wk/lin_ETG_adiab_i.m
rename to wk/parameters/lin_ETG_adiab_i.m
index 1c3d907e7082d74f0629cdd4f99507184b42f221..b3b3d7ed66d0eca85839704488aea817953df80b 100644
--- a/wk/lin_ETG_adiab_i.m
+++ b/wk/parameters/lin_ETG_adiab_i.m
@@ -1,23 +1,5 @@
-%% QUICK RUN SCRIPT
-% This script creates a directory in /results and runs a simulation directly
-% from the Matlab framework. It is meant to run only small problems in linear
-% for benchmarking and debugging purposes since it makes Matlab "busy".
-
-%% Set up the paths for the necessary Matlab modules
-gyacomodir = pwd;
-gyacomodir = gyacomodir(1:end-2);
-addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
-addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
-addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
-addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
-
 %% Set simulation parameters
-SIMID   = 'ETG_adiab_i';  % Name of the simulation
-RUN = 1; % To run or just to load
-default_plots_options
-EXECNAME = 'gyacomo23_sp'; % single precision
-% EXECNAME = 'gyacomo23_dp'; % double precision
-% EXECNAME = 'gyacomo23_debug'; % single precision
+SIMID   = 'lin_ETG_adiab_i';  % Name of the simulation
 
 %% Set up physical parameters
 CLUSTER.TIME = '99:00:00';  % Allocation time hh:mm:ss
@@ -113,50 +95,4 @@ NOISE0  = 1.0e-8; % Initial noise amplitude
 BCKGD0  = 0.0e-8;    % Initial background
 k_gB   = 1.0;     % Magnetic gradient strength
 k_cB   = 1.0;     % Magnetic curvature strength
-COLL_KCUT = 1; % Cutoff for collision operator
-
-%%-------------------------------------------------------------------------
-%% RUN
-setup
-% system(['rm fort*.90']);
-% 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 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
-    RUN  =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 3 2 0;'];
-    % RUN  =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
-    MVOUT='cd ../../../wk;';
-    system([MVIN,RUN,MVOUT]);
-end
-
-%% Analysis
-% load
-filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
-LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
-FIGDIR   = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
-% Load outputs from jobnummin up to jobnummax
-J0 = 0; J1 = 0;
-data = {}; % Initialize data as an empty cell array
-% load grids, inputs, and time traces
-data = compile_results_low_mem(data,LOCALDIR,J0,J1); 
-
-
-if 1 % Activate or not
-%% plot mode evolution and growth rates
-% Load phi
-[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
-options.NORMALIZED = 0; 
-options.TIME   = data.Ts3D;
- % Time window to measure the growth of kx/ky modes
-options.KX_TW  = [0.2 1]*data.Ts3D(end);
-options.KY_TW  = [0.2 1]*data.Ts3D(end);
-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.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
-end
-
-
+COLL_KCUT = 1; % Cutoff for collision operator
\ No newline at end of file
diff --git a/wk/lin_Entropy.m b/wk/parameters/lin_Entropy.m
similarity index 63%
rename from wk/lin_Entropy.m
rename to wk/parameters/lin_Entropy.m
index 6197a97f2f5691770d7a97e6eaf2453d2071c865..18a3c2ad714091e40d46d2afcf811db471ab4952 100644
--- a/wk/lin_Entropy.m
+++ b/wk/parameters/lin_Entropy.m
@@ -1,23 +1,5 @@
-%% QUICK RUN SCRIPT
-% This script creates a directory in /results and runs a simulation directly
-% from the Matlab framework. It is meant to run only small problems in linear
-% for benchmarking and debugging purposes since it makes Matlab "busy".
-
-%% Set up the paths for the necessary Matlab modules
-gyacomodir = pwd;
-gyacomodir = gyacomodir(1:end-2);
-addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
-addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
-addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
-addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
-
 %% Set simulation parameters
-SIMID   = 'ETG_adiab_i';  % Name of the simulation
-RUN = 1; % To run or just to load
-default_plots_options
-EXECNAME = 'gyacomo23_sp'; % single precision
-% EXECNAME = 'gyacomo23_dp'; % double precision
-% EXECNAME = 'gyacomo23_debug'; % single precision
+SIMID   = 'lin_Entropy';  % Name of the simulation
 
 %% Set up physical parameters
 CLUSTER.TIME = '99:00:00';  % Allocation time hh:mm:ss
@@ -114,49 +96,3 @@ BCKGD0  = 0.0e-8;    % Initial background
 k_gB   = 1.0;     % Magnetic gradient strength
 k_cB   = 1.0;     % Magnetic curvature strength
 COLL_KCUT = 1; % Cutoff for collision operator
-
-%%-------------------------------------------------------------------------
-%% RUN
-setup
-% system(['rm fort*.90']);
-% 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 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0;'];
-    % RUN  =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 3 2 0;'];
-    % RUN  =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
-    MVOUT='cd ../../../wk;';
-    system([MVIN,RUN,MVOUT]);
-end
-
-%% Analysis
-% load
-filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
-LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
-FIGDIR   = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
-% Load outputs from jobnummin up to jobnummax
-J0 = 0; J1 = 0;
-data = {}; % Initialize data as an empty cell array
-% load grids, inputs, and time traces
-data = compile_results_low_mem(data,LOCALDIR,J0,J1); 
-
-
-if 1 % Activate or not
-%% plot mode evolution and growth rates
-% Load phi
-[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
-options.NORMALIZED = 0; 
-options.TIME   = data.Ts3D;
- % Time window to measure the growth of kx/ky modes
-options.KX_TW  = [0.2 1]*data.Ts3D(end);
-options.KY_TW  = [0.2 1]*data.Ts3D(end);
-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.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
-end
-
-
diff --git a/wk/lin_ITG_salpha.m b/wk/parameters/lin_ITG.m
similarity index 62%
rename from wk/lin_ITG_salpha.m
rename to wk/parameters/lin_ITG.m
index 67c91f8fa728f839ccfd890da7ffe36c2d7b5d7e..6519a58f31a5d685e45ac6330efffde09de31414 100644
--- a/wk/lin_ITG_salpha.m
+++ b/wk/parameters/lin_ITG.m
@@ -1,23 +1,5 @@
-%% QUICK RUN SCRIPT
-% This script creates a directory in /results and runs a simulation directly
-% from the Matlab framework. It is meant to run only small problems in linear
-% for benchmarking and debugging purposes since it makes Matlab "busy".
-
-%% Set up the paths for the necessary Matlab modules
-gyacomodir = pwd;
-gyacomodir = gyacomodir(1:end-2);
-addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
-addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
-addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
-addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
-
 %% Set simulation parameters
 SIMID = 'lin_ITG'; % Name of the simulation
-RUN = 1; % To run or just to load
-default_plots_options
-% EXECNAME = 'gyacomo23_dp'; % single precision
-EXECNAME = 'gyacomo23_sp'; % double precision
-% EXECNAME = 'gyacomo23_debug'; % double precision
 
 %% Set up physical parameters
 CLUSTER.TIME = '99:00:00';  % Allocation time hh:mm:ss
@@ -114,55 +96,4 @@ S_DELTA = 0.0;
 S_ZETA  = 0.0;
 PB_PHASE= 0;
 ADIAB_I = 0;
-MHD_PD  = 0;
-%%-------------------------------------------------------------------------
-%% RUN
-setup
-% system(['rm fort*.90']);
-% Run linear simulation
-if RUN
-    MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
-%     RUN  =['time mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
-    % RUN  =['mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
-%     RUN  =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
-    RUN  =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
-    MVOUT='cd ../../../wk;';
-    system([MVIN,RUN,MVOUT]);
-end
-
-%% Analysis
-% load
-filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
-LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
-FIGDIR   = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
-% Load outputs from jobnummin up to jobnummax
-J0 = 0; J1 = 0;
-data = {}; % Initialize data as an empty cell array
-% load grids, inputs, and time traces
-data = compile_results_low_mem(data,LOCALDIR,J0,J1); 
-
-if 0
-%% Plot heat flux evolution
-figure
-semilogy(data.Ts0D,data.HFLUX_X);
-xlabel('$tc_s/R$'); ylabel('$Q_x$');
-end
-if 1 % Activate or not
-%% plot mode evolution and growth rates
-% Load phi
-[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
-options.NORMALIZED = 0; 
-options.TIME   = data.Ts3D;
- % Time window to measure the growth of kx/ky modes
-options.KX_TW  = [0.2 1]*data.Ts3D(end);
-options.KY_TW  = [0.2 1]*data.Ts3D(end);
-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.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
-end
-
-
-
+MHD_PD  = 0;
\ No newline at end of file
diff --git a/wk/lin_Ivanov.m b/wk/parameters/lin_Ivanov.m
similarity index 57%
rename from wk/lin_Ivanov.m
rename to wk/parameters/lin_Ivanov.m
index 5ecf60b0002a993319ed9fe98ae466292fc7544e..b1b7e7eb78be86b697d965abddffb1946594e2be 100644
--- a/wk/lin_Ivanov.m
+++ b/wk/parameters/lin_Ivanov.m
@@ -1,24 +1,6 @@
-%% QUICK RUN SCRIPT
-% This script creates a directory in /results and runs a simulation directly
-% from the Matlab framework. It is meant to run only small problems in linear
-% for benchmarking and debugging purposes since it makes Matlab "busy".
-
-%% Set up the paths for the necessary Matlab modules
-gyacomodir = pwd;
-gyacomodir = gyacomodir(1:end-2);
-addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
-addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
-addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
-addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
-
 %% Set simulation parameters
-SIMID = 'dbg'; % Name of the simulation
-RUN = 1; % To run or just to load
-default_plots_options
-% EXECNAME = 'gyacomo23_dp'; % double precision
-% To compile single precision gyacomo do 'make clean; make sp' in the /gyacomo folder
-EXECNAME = 'gyacomo23_sp'; % single precision
-% EXECNAME = 'gyacomo23_debug'; % single precision
+SIMID = 'lin_Ivanov'; % Name of the simulation
+
 %% Set up physical parameters
 CLUSTER.TIME = '99:00:00';  % Allocation time hh:mm:ss
 TAU = 1e-3;                  % e/i temperature ratio
@@ -111,70 +93,3 @@ BCKGD0  = 0.0;    % Initial background
 k_gB   = 1.0;     % Magnetic gradient strength
 k_cB   = 1.0;     % Magnetic curvature strength
 COLL_KCUT = 1.0; % Cutoff for collision operator
-
-%%-------------------------------------------------------------------------
-%% RUN
-setup
-% system(['rm fort*.90']);
-% 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 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0;'];
-%     RUN  =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
-%     RUN  =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
-    MVOUT='cd ../../../wk;';
-    system([MVIN,RUN,MVOUT]);
-end
-
-%% Analysis
-% load
-filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
-LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
-FIGDIR   = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
-% Load outputs from jobnummin up to jobnummax
-J0 = 0; J1 = 0;
-data = {}; % Initialize data as an empty cell array
-% load grids, inputs, and time traces
-data = compile_results_low_mem(data,LOCALDIR,J0,J1); 
-
-if 0
-%% Plot heat flux evolution
-figure
-semilogy(data.Ts0D,data.HFLUX_X);
-xlabel('$tc_s/R$'); ylabel('$Q_x$');
-end
-if 1 % Activate or not
-%% plot mode evolution and growth rates
-% Load phi
-[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
-options.NORMALIZED = 0; 
-options.TIME   = data.Ts3D;
- % Time window to measure the growth of kx/ky modes
-options.KX_TW  = [0.5 1]*data.Ts3D(end);
-options.KY_TW  = [0.5 1]*data.Ts3D(end);
-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.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
-end
-
-if 0
-%% Hermite-Laguerre spectrum
-[data.Napjz, data.Ts3D] = compile_results_3Da(LOCALDIR,J0,J1,'Napjz');
-options.ST         = 1;
-options.NORMALIZED = 0;
-options.LOGSCALE   = 1;
-options.FILTER     = 0; %filter the 50% time-average of the spectrum from
-fig = show_moments_spectrum(data,options);
-end
-
-% filename = '/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/Ivanov_2020_fig2_kT_0.26_chi_0.1.txt';
-% % filename = '/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/Ivanov_2020_fig2_kT_1_chi_0.1.txt';
-% dIV = load(filename);
-% 
-% figure 
-% plot(dIV(:,1),2*dIV(:,2))
-
diff --git a/wk/lin_KBM.m b/wk/parameters/lin_KBM.m
similarity index 61%
rename from wk/lin_KBM.m
rename to wk/parameters/lin_KBM.m
index 40b3db95a90a41493ddc721ff55cc747b698e9df..dab2fe08308a13edcecfe55bd7658366c7aec5c2 100644
--- a/wk/lin_KBM.m
+++ b/wk/parameters/lin_KBM.m
@@ -1,22 +1,5 @@
-%% QUICK RUN SCRIPT
-% This script creates a directory in /results and runs a simulation directly
-% from the Matlab framework. It is meant to run only small problems in linear
-% for benchmarking and debugging purposes since it makes Matlab "busy".
-
-%% Set up the paths for the necessary Matlab modules
-gyacomodir = pwd;
-gyacomodir = gyacomodir(1:end-2);
-addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
-addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
-addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
-addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
-
 %% Set simulation parameters
 SIMID   = 'lin_KBM';  % Name of the simulation
-RUN = 1; % To run or just to load
-default_plots_options
-EXECNAME = 'gyacomo23_sp'; % single precision
-% EXECNAME = 'gyacomo23_dp'; % double precision
 
 %% Set up physical parameters
 CLUSTER.TIME = '99:00:00';  % Allocation time hh:mm:ss
@@ -108,61 +91,3 @@ BCKGD0  = 1.0e-5;    % Initial background
 k_gB   = 1.0;     % Magnetic gradient strength
 k_cB   = 1.0;     % Magnetic curvature strength
 COLL_KCUT = 1; % Cutoff for collision operator
-
-%%-------------------------------------------------------------------------
-%% RUN
-setup
-% system(['rm fort*.90']);
-% 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 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
-%     RUN  =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
-%     RUN  =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
-    MVOUT='cd ../../../wk;';
-    system([MVIN,RUN,MVOUT]);
-end
-
-%% Analysis
-% load
-filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
-LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
-FIGDIR   = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
-% Load outputs from jobnummin up to jobnummax
-J0 = 0; J1 = 0;
-data = {}; % Initialize data as an empty cell array
-% load grids, inputs, and time traces
-data = compile_results_low_mem(data,LOCALDIR,J0,J1); 
-
-
-if 0 % Activate or not
-%% plot mode evolution and growth rates
-% Load phi
-[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
-options.NORMALIZED = 0; 
-options.TIME   = data.Ts3D;
- % Time window to measure the growth of kx/ky modes
-options.KX_TW  = [0.2 1]*data.Ts3D(end);
-options.KY_TW  = [0.2 1]*data.Ts3D(end);
-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.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
-end
-
-if 1
-%% Ballooning plot
-[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
-if data.inputs.BETA > 0
-[data.PSI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'psi');
-end
-options.time_2_plot = [120];
-options.kymodes     = [0.25];
-options.normalized  = 1;
-% options.field       = 'phi';
-fig = plot_ballooning(data,options);
-end
-