diff --git a/wk/AUG_ky_PJ_scan.m b/wk/AUG_ky_PJ_scan.m
deleted file mode 100644
index 2b562a895e1d8c90af3472024071db00a88bbe74..0000000000000000000000000000000000000000
--- a/wk/AUG_ky_PJ_scan.m
+++ /dev/null
@@ -1,224 +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_AUG_scan'; % Name of the simulation
-RERUN = 1; % rerun if the data does not exist
-RUN = 1;
-K_Ne = 34; % ele Density '''
-K_Te = 56; % ele Temperature '''
-K_Ni = 34; % ion Density gradient drive
-K_Ti = 21; % 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 = 1.52e-4; % electron plasma beta
-MHD_PD = 0; % MHD pressure drift
-% Geometry
-GEOMETRY= 'miller';
-% GEOMETRY= 's-alpha';
-% 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.38; % 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= 'miller';
- 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;
- MHD_PD = 0; % MHD pressure drift
- %% 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;
- ADIAB_I = 0; % adiabatic ion model
- %% 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,' 1 2 2 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/parameters/lin_AUG.m b/wk/parameters/lin_AUG.m
new file mode 100644
index 0000000000000000000000000000000000000000..2d3e19bc32aaabeb1981c5568b7de5cb9d66851f
--- /dev/null
+++ b/wk/parameters/lin_AUG.m
@@ -0,0 +1,98 @@
+%% Set simulation parameters
+SIMID = 'lin_AUG'; % Name of the simulation
+
+%% 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
\ No newline at end of file