diff --git a/matlab/MOLI_kperp_scan.m b/matlab/MOLI_kperp_scan.m
index c3c8d4e..9920550 100644
--- a/matlab/MOLI_kperp_scan.m
+++ b/matlab/MOLI_kperp_scan.m
@@ -1,166 +1,166 @@
 %% Run MOLI for a scan in kperp
 %% Move to MOLI workspace
 cd ../../MoliSolver/MOLI/workspace/
 %% Add paths
 my_paths
 
 %% Directories
 ROOT = '/home/ahoffman/Documents/MoliSolver';
 options.dirs.COSOlverdir = fullfile(ROOT,'COSOlver');
 options.dirs.MOLIdir = fullfile(ROOT,'MOLI');
 
 %% MOLI Physical and Main Parameters
 
 % Solve DK AND/OR GK Linear Moment Hierarchy.
 options.DKI  = 0;   % First-order DK
 options.DKII = 0;   % Second-order DK
 options.GK   = 1;   % Gyrokinetic GK
 options.EM   = 0;   % Include Electromagnetic effects (only for GK=1)
 options.GD   = 0;   % Gyro-Drift
 options.GDI  = 0;   % First/second order
 
 % Solve MOLI
 options.MOLI = 1;    % 1 -> Solve MOLI, 0 -> off
 
 % MOLI Solver
 options.solver.solver = 3;
 
 % MOLI Linear Fit Solver
 options.LinFitSolver = 2;
 
 %% Main parameter scan
 
 % Closure by truncation
 params.Pmaxi = GRID.pmaxi;           % parallel ion Hermite moments
 params.Jmaxi = GRID.jmaxi;           % perpendicular ion Laguerre moments
 params.Pmaxe = GRID.pmaxe;           % parallel electron Hermite moments
 params.Jmaxe = GRID.jmaxe;           % perpendicular electron Laguerre moments
 
 % w/wo ions
 options.ions = 1;           % if ions are present -> 1, 0 otherwise
 
 % Adiabatic electrons
 options.electrons = 1;      % 0 ->  adiabatic electrons, 1 no adiabatic electrons
 
 % w/wo soundwaves
 options.sw = 1;             % 1 -> sound waves on, 0 -> put ion parallel velocity row/column to 0
 
 %% Collision Operator Models and COSOlver Input Parameters
 options.collI  = MODEL.CO;         % collI=-2 -> Dougherty, -1 -> COSOlver, 0 -> Lenard-Bernstein, other -> hyperviscosity exponent
 options.collGK = 0;         % collDKGK =1 -> GK collision operator, else DK collision operator
 options.COSOlver.GKE = 0;
 options.COSOlver.GKI = 0;
 
 % COSOlver Input Parameters (if collI = -1 only)
 options.COSOlver.eecolls = 1;	   % 1 -> electron-electron collisions, 0 -> off
 options.COSOlver.iicolls = 1;      % 1 -> ion-ion collisions, 0 -> off
 options.COSOlver.eicolls = 1;      % 1 -> electron-ion collisions (e-i) on, 0 -> off
 options.COSOlver.iecolls = 1;      % 1 -> ion-electron collisions (i-e) on, 0 -> off
 
 % Collisional Coulomb sum bounds (only if collI = -1, i.e. Coulomb)
 options.COSOlver.lmaxx = 10;                        % upper bound collision operator first sum first species
 options.COSOlver.kmaxx = 10;                        % upper bound collision operator second sum first species
 options.COSOlver.nmaxx = options.COSOlver.lmaxx;    % upper bound collision operator first sum second species
 options.COSOlver.qmaxx = options.COSOlver.kmaxx;    % upper bound collision operator second sum second species
 
 % Collsion FLR sum bounds
 options.COSOlver.nemaxxFLR = 0;         % upper bound FLR electron collision
 options.COSOlver.nimaxxFLR = 0;         % upper bound FLR ion collision
 
 % Collision Operator Model
 % Set electron/ion test and back-reaction model operator
 %
 %  0 => Coulomb Collisions
 options.COSOlver.ETEST = 1; % 0 --> Buffer Operator, 1 --> Coulomb, 2 --> Lorentz
 options.COSOlver.EBACK = 1;
 options.COSOlver.ITEST = 1;
 options.COSOlver.IBACK = 1;
 options.COSOlver.ESELF = 1;
 options.COSOlver.ISELF = 1;
 
 options.COSOlver.OVERWRITE = 0;    % overwrite collisional matrices even if exists
 
 options.COSOlver.cmd   = 'mpirun -np 6 ./bin/CO 2 2 2';
 %% Physical Parameters
 
 % Toroidal effects
 options.magnetic        = 1;           % 1-> Add toroidal magnetic gradient drift resonance effects
 
 % Physical Parameters
 params.tau              = MODEL.tau_i; % Ti/Te
 params.nu               = MODEL.nu;    % electron/ion collision frequency ... only for nu/ omega_pe < nuoveromegapemax (electron plasma frequency) [See Banks et al. (2017)]
 params.nuoveromegapemax = inf;         % Maximum ratio between electron/ion collision frequency and electron plasma frequency [See Banks et al. (2017)]. Set to inf if not desired !!!
 params.mu               = MODEL.sigma_e;   % sqrt(m_e/m_i)
 params.kpar             = 0.0;         % normalized parallel wave number to the major radius
 params.kperp            = GRID.kzmin;  % normalized perpendicular toroidal wave number to the soundLarmor radius. Note: If ions ==0 (e.g. EPW), kperp --> b
 params.kr               = GRID.krmin;  % Radial component of perpendicular vector
 params.alphaD           = 0.0;         % (k*Debye length)^2
 params.Rn               = MODEL.eta_n; % Major Radius / Background density gradient length
 params.RTe              = MODEL.eta_T; % Major Radius * normalized kperp / Background electron temperature gradient length
 params.RTi              = MODEL.eta_T; % Major Radius * normalized kperp / Background ion temperature gradient length
 params.Rphi             = 0.0;         % Major Radius * normalized kperp / Background potentiel gradient length [presence of shear] - only for GK
 params.betae            = 1e-6;        % Electron Beta plasma.
 
 params.rhostar          = 1e-5;        % sound Larmor Radius/Major Radius ~ sqrt(Te)/(R_0*B).
-params.n0               = INITIAL.initback_moments;        % initial density perturbation
+params.n0               = INITIAL.init_background;        % initial density perturbation
 
 params.gradB            = MODEL.eta_B;         % Magnetic field gradient
 params.curvB            = MODEL.eta_B;         % Curvature of B
 params.trappB           = 0.0;         % Trapping term
 
 %% MOLI Options
 
 % Save data in results dir
 options.save    = 0;
 options.verbose = 0;
 options.dbg     = 0;
 
 options.DR        = 0;      % 1 -> Solve kinetic dispersion relation,
 options.KineticDR = 0;      % Solve kinetic dispersion relation (Landau integral) for the given theory
 
 % Compute the kinetic susceptibility for EPW only
 options.SPTBLTY = 0;
 
 options.nharm   = 1;        % Number of harmonics in disp. rel. 1 and 4
 
 wlim = 5.0;
 options.DRsolver.wr_min = -wlim;     % Minimum of real part.
 options.DRsolver.wr_max =  wlim;     % Maximum of real part.
 options.DRsolver.wi_min = -wlim;     % Minimum of imag part.
 options.DRsolver.wi_max =  wlim;     % Maximum of imag part.
 options.DRsolver.nw     =  300;      % Grid resolution
 
 % Disp. Rel. Options
 options.FLRmodel    = 0;   % 1 -> Truncated Laguerre, 0 -> Exact representation
 options.FluidLandau = 0;   % 1 -> Add Landau Fluid Closure to Fluid Dispersion Relation, 0 -> off
 options.deltaLandau = 0;   % 1 -> Hammet-Perkins closure on, 0 -> off
 
 % Fluid dispersion relation
 options.FluidDR = 0;	     % Solve annamaria's fluid equations
 options.Fluid.sITGmm = 0;
 
 % Define scan parameters
 options.fscan = 1;         % 1 -> peform scan over scan.list, 0-> off
 
 options.scan.list = {'kperp'};% List of scan parameters. If empty, solve MOLI with params
 
 options.scan.kperp.array  = linspace(GRID.kzmin, GRID.kzmax, GRID.nkz);
 
 % Time-Evolution Problem [Solver==3] ...
 options.solver.TimeSolver.dt         = BASIC.dt;		% timestep of time evolution (R/c_s or 1/(k v/the) units)
 options.solver.TimeSolver.tmax       = BASIC.tmax;
 options.solver.TimeSolver.Trun       = BASIC.tmax;		  % total time to run time evolution
 options.solver.TimeSolver.t_fit_min  = 0.05;    % Phase-Mixing fit Lower time limit
 options.solver.TimeSolver.t_fit_max  = 8;       % Phase-Mixing fit Upper time limit
 options.solver.TimeSolver.en_fit_min = 0.15;    % Entropy Mode fit Lower time limit
 options.solver.TimeSolver.en_fit_max = 0.3;     % Entropy Mode fit Upper time limit
 options.solver.TimeSolver.movie      = 0;       % Display movie if 1, last frame otherwise
 options.solver.TimeSolver.save       = 0;       % 1 --> save during fscan, Warning: memory storage
 
 
 %% Run MOLI
 % Solve the MOLI
 [results,params,options] = MOLI_Control(params,options);
 
 %% Return to HeLaZ workspace
 cd ../../../HeLaZ/wk
diff --git a/matlab/MOLI_time_solver.m b/matlab/MOLI_time_solver.m
index d9b4d80..c6c6ed3 100644
--- a/matlab/MOLI_time_solver.m
+++ b/matlab/MOLI_time_solver.m
@@ -1,164 +1,164 @@
 %% Run MOLI for a time evolution of the moments at a given kperp
 %% Move to MOLI workspace
 cd ../../MoliSolver/MOLI/workspace/
 %% Add paths
 my_paths
 
 %% Directories
 ROOT = '/home/ahoffman/Documents/MoliSolver';
 options.dirs.COSOlverdir = fullfile(ROOT,'COSOlver');
 options.dirs.MOLIdir = fullfile(ROOT,'MOLI');
 
 %% MOLI Physical and Main Parameters
 
 % Solve DK AND/OR GK Linear Moment Hierarchy.
 options.DKI  = 0;   % First-order DK
 options.DKII = 0;   % Second-order DK
 options.GK   = 1;   % Gyrokinetic GK
 options.EM   = 0;   % Include Electromagnetic effects (only for GK=1)
 options.GD   = 0;   % Gyro-Drift
 options.GDI  = 0;   % First/second order
 
 % Solve MOLI
 options.MOLI = 1;    % 1 -> Solve MOLI, 0 -> off
 
 % MOLI Solver
 options.solver.solver = 3;
 
 % MOLI Linear Fit Solver
 options.LinFitSolver = 0;
 
 %% Main parameter scan
 
 % Closure by truncation
 params.Pmaxi = GRID.pmaxi;           % parallel ion Hermite moments
 params.Jmaxi = GRID.jmaxi;           % perpendicular ion Laguerre moments
 params.Pmaxe = GRID.pmaxe;           % parallel electron Hermite moments
 params.Jmaxe = GRID.jmaxe;           % perpendicular electron Laguerre moments
 
 % w/wo ions
 options.ions = 1;           % if ions are present -> 1, 0 otherwise
 
 % Adiabatic electrons
 options.electrons = 1;      % 0 ->  adiabatic electrons, 1 no adiabatic electrons
 
 % w/wo soundwaves
 options.sw = 1;             % 1 -> sound waves on, 0 -> put ion parallel velocity row/column to 0
 
 %% Collision Operator Models and COSOlver Input Parameters
 options.collI  = MODEL.CO;         % collI=-2 -> Dougherty, -1 -> COSOlver, 0 -> Lenard-Bernstein, other -> hyperviscosity exponent
 options.collGK = 0;         % collDKGK =1 -> GK collision operator, else DK collision operator
 options.COSOlver.GKE = 0;
 options.COSOlver.GKI = 0;
 
 % COSOlver Input Parameters (if collI = -1 only)
 options.COSOlver.eecolls = 1;	   % 1 -> electron-electron collisions, 0 -> off
 options.COSOlver.iicolls = 1;      % 1 -> ion-ion collisions, 0 -> off
 options.COSOlver.eicolls = 1;      % 1 -> electron-ion collisions (e-i) on, 0 -> off
 options.COSOlver.iecolls = 1;      % 1 -> ion-electron collisions (i-e) on, 0 -> off
 
 % Collisional Coulomb sum bounds (only if collI = -1, i.e. Coulomb)
 options.COSOlver.lmaxx = 10;                        % upper bound collision operator first sum first species
 options.COSOlver.kmaxx = 10;                        % upper bound collision operator second sum first species
 options.COSOlver.nmaxx = options.COSOlver.lmaxx;    % upper bound collision operator first sum second species
 options.COSOlver.qmaxx = options.COSOlver.kmaxx;    % upper bound collision operator second sum second species
 
 % Collsion FLR sum bounds
 options.COSOlver.nemaxxFLR = 0;         % upper bound FLR electron collision
 options.COSOlver.nimaxxFLR = 0;         % upper bound FLR ion collision
 
 % Collision Operator Model
 % Set electron/ion test and back-reaction model operator
 %
 %  0 => Coulomb Collisions
 options.COSOlver.ETEST = 1; % 0 --> Buffer Operator, 1 --> Coulomb, 2 --> Lorentz
 options.COSOlver.EBACK = 1;
 options.COSOlver.ITEST = 1;
 options.COSOlver.IBACK = 1;
 options.COSOlver.ESELF = 1;
 options.COSOlver.ISELF = 1;
 
 options.COSOlver.OVERWRITE = 0;    % overwrite collisional matrices even if exists
 
 options.COSOlver.cmd   = 'mpirun -np 6 ./bin/CO 2 2 2';
 %% Physical Parameters
 
 % Toroidal effects
 options.magnetic        = 1;           % 1-> Add toroidal magnetic gradient drift resonance effects
 
 % Physical Parameters
 params.tau              = MODEL.tau_i; % Ti/Te
 params.nu               = MODEL.nu;    % electron/ion collision frequency ... only for nu/ omega_pe < nuoveromegapemax (electron plasma frequency) [See Banks et al. (2017)]
 params.nuoveromegapemax = inf;         % Maximum ratio between electron/ion collision frequency and electron plasma frequency [See Banks et al. (2017)]. Set to inf if not desired !!!
 params.mu               = MODEL.sigma_e;   % sqrt(m_e/m_i)
 params.kpar             = 0.0;         % normalized parallel wave number to the major radius
 params.kperp            = kz_MOLI;  % normalized perpendicular toroidal wave number to the soundLarmor radius. Note: If ions ==0 (e.g. EPW), kperp --> b
 params.kr               = kr_MOLI;  % Radial component of perpendicular vector
 params.alphaD           = 0.0;         % (k*Debye length)^2
 params.Rn               = MODEL.eta_n; % Major Radius / Background density gradient length
 params.RTe              = MODEL.eta_T; % Major Radius * normalized kperp / Background electron temperature gradient length
 params.RTi              = MODEL.eta_T; % Major Radius * normalized kperp / Background ion temperature gradient length
 params.Rphi             = 0.0;         % Major Radius * normalized kperp / Background potentiel gradient length [presence of shear] - only for GK
 params.betae            = 1e-6;        % Electron Beta plasma.
 
 params.rhostar          = 1e-5;        % sound Larmor Radius/Major Radius ~ sqrt(Te)/(R_0*B).
-params.n0               = INITIAL.initback_moments;        % initial density perturbation
+params.n0               = INITIAL.init_background;        % initial density perturbation
 
 params.gradB            = MODEL.eta_B;         % Magnetic field gradient
 params.curvB            = MODEL.eta_B;         % Curvature of B
 params.trappB           = 0.0;         % Trapping term
 
 %% MOLI Options
 
 % Save data in results dir
 options.save    = 0;
 options.verbose = 0;
 options.dbg     = 0;
 
 options.DR        = 0;      % 1 -> Solve kinetic dispersion relation,
 options.KineticDR = 0;      % Solve kinetic dispersion relation (Landau integral) for the given theory
 
 % Compute the kinetic susceptibility for EPW only
 options.SPTBLTY = 0;
 
 options.nharm   = 1;        % Number of harmonics in disp. rel. 1 and 4
 
 wlim = 5.0;
 options.DRsolver.wr_min = -wlim;     % Minimum of real part.
 options.DRsolver.wr_max =  wlim;     % Maximum of real part.
 options.DRsolver.wi_min = -wlim;     % Minimum of imag part.
 options.DRsolver.wi_max =  wlim;     % Maximum of imag part.
 options.DRsolver.nw     =  300;      % Grid resolution
 
 % Disp. Rel. Options
 options.FLRmodel    = 0;   % 1 -> Truncated Laguerre, 0 -> Exact representation
 options.FluidLandau = 0;   % 1 -> Add Landau Fluid Closure to Fluid Dispersion Relation, 0 -> off
 options.deltaLandau = 0;   % 1 -> Hammet-Perkins closure on, 0 -> off
 
 % Fluid dispersion relation
 options.FluidDR = 0;	     % Solve annamaria's fluid equations
 options.Fluid.sITGmm = 0;
 
 % Define scan parameters
 options.fscan = 0;         % 1 -> peform scan over scan.list, 0-> off
 
 options.scan.list = {};% List of scan parameters. If empty, solve MOLI with params
 
 % Time-Evolution Problem [Solver==3] ...
 options.solver.TimeSolver.dt         = BASIC.dt;		% timestep of time evolution (R/c_s or 1/(k v/the) units)
 options.solver.TimeSolver.tmax       = BASIC.tmax;
 options.solver.TimeSolver.Trun       = BASIC.tmax;		  % total time to run time evolution
 options.solver.TimeSolver.t_fit_min  = 0.05;    % Phase-Mixing fit Lower time limit
 options.solver.TimeSolver.t_fit_max  = 8;       % Phase-Mixing fit Upper time limit
 options.solver.TimeSolver.en_fit_min = 0.15;    % Entropy Mode fit Lower time limit
 options.solver.TimeSolver.en_fit_max = 0.3;     % Entropy Mode fit Upper time limit
 options.solver.TimeSolver.movie      = 0;       % Display movie if 1, last frame otherwise
 options.solver.TimeSolver.save       = 0;       % 1 --> save during fscan, Warning: memory storage
 
 
 %% Run MOLI
 % Solve the MOLI
 [results,params,options] = MOLI_Control(params,options);
 
 %% Return to HeLaZ workspace
 cd ../../../HeLaZ/wk
diff --git a/matlab/MOLI_time_solver_2D.m b/matlab/MOLI_time_solver_2D.m
index 02f8133..a356c5e 100644
--- a/matlab/MOLI_time_solver_2D.m
+++ b/matlab/MOLI_time_solver_2D.m
@@ -1,164 +1,164 @@
 %% Run MOLI for a time evolution of the moments at a given kperp
 %% Move to MOLI workspace
 cd ../../MoliSolver/MOLI/workspace/
 %% Add paths
 my_paths
 
 %% Directories
 ROOT = '/home/ahoffman/Documents/MoliSolver';
 options.dirs.COSOlverdir = fullfile(ROOT,'COSOlver');
 options.dirs.MOLIdir = fullfile(ROOT,'MOLI');
 
 %% MOLI Physical and Main Parameters
 
 % Solve DK AND/OR GK Linear Moment Hierarchy.
 options.DKI  = 0;   % First-order DK
 options.DKII = 0;   % Second-order DK
 options.GK   = 1;   % Gyrokinetic GK
 options.EM   = 0;   % Include Electromagnetic effects (only for GK=1)
 options.GD   = 0;   % Gyro-Drift
 options.GDI  = 0;   % First/second order
 
 % Solve MOLI
 options.MOLI = 1;    % 1 -> Solve MOLI, 0 -> off
 
 % MOLI Solver
 options.solver.solver = 3;
 
 % MOLI Linear Fit Solver
 options.LinFitSolver = 0;
 
 %% Main parameter scan
 
 % Closure by truncation
 params.Pmaxi = GRID.pmaxi;           % parallel ion Hermite moments
 params.Jmaxi = GRID.jmaxi;           % perpendicular ion Laguerre moments
 params.Pmaxe = GRID.pmaxe;           % parallel electron Hermite moments
 params.Jmaxe = GRID.jmaxe;           % perpendicular electron Laguerre moments
 
 % w/wo ions
 options.ions = 1;           % if ions are present -> 1, 0 otherwise
 
 % Adiabatic electrons
 options.electrons = 1;      % 0 ->  adiabatic electrons, 1 no adiabatic electrons
 
 % w/wo soundwaves
 options.sw = 1;             % 1 -> sound waves on, 0 -> put ion parallel velocity row/column to 0
 
 %% Collision Operator Models and COSOlver Input Parameters
 options.collI  = MODEL.CO;         % collI=-2 -> Dougherty, -1 -> COSOlver, 0 -> Lenard-Bernstein, other -> hyperviscosity exponent
 options.collGK = 0;         % collDKGK =1 -> GK collision operator, else DK collision operator
 options.COSOlver.GKE = 0;
 options.COSOlver.GKI = 0;
 
 % COSOlver Input Parameters (if collI = -1 only)
 options.COSOlver.eecolls = 1;	   % 1 -> electron-electron collisions, 0 -> off
 options.COSOlver.iicolls = 1;      % 1 -> ion-ion collisions, 0 -> off
 options.COSOlver.eicolls = 1;      % 1 -> electron-ion collisions (e-i) on, 0 -> off
 options.COSOlver.iecolls = 1;      % 1 -> ion-electron collisions (i-e) on, 0 -> off
 
 % Collisional Coulomb sum bounds (only if collI = -1, i.e. Coulomb)
 options.COSOlver.lmaxx = 10;                        % upper bound collision operator first sum first species
 options.COSOlver.kmaxx = 10;                        % upper bound collision operator second sum first species
 options.COSOlver.nmaxx = options.COSOlver.lmaxx;    % upper bound collision operator first sum second species
 options.COSOlver.qmaxx = options.COSOlver.kmaxx;    % upper bound collision operator second sum second species
 
 % Collsion FLR sum bounds
 options.COSOlver.nemaxxFLR = 0;         % upper bound FLR electron collision
 options.COSOlver.nimaxxFLR = 0;         % upper bound FLR ion collision
 
 % Collision Operator Model
 % Set electron/ion test and back-reaction model operator
 %
 %  0 => Coulomb Collisions
 options.COSOlver.ETEST = 1; % 0 --> Buffer Operator, 1 --> Coulomb, 2 --> Lorentz
 options.COSOlver.EBACK = 1;
 options.COSOlver.ITEST = 1;
 options.COSOlver.IBACK = 1;
 options.COSOlver.ESELF = 1;
 options.COSOlver.ISELF = 1;
 
 options.COSOlver.OVERWRITE = 0;    % overwrite collisional matrices even if exists
 
 options.COSOlver.cmd   = 'mpirun -np 6 ./bin/CO 2 2 2';
 %% Physical Parameters
 
 % Toroidal effects
 options.magnetic        = 1;           % 1-> Add toroidal magnetic gradient drift resonance effects
 
 % Physical Parameters
 params.tau              = MODEL.tau_i; % Ti/Te
 params.nu               = MODEL.nu;    % electron/ion collision frequency ... only for nu/ omega_pe < nuoveromegapemax (electron plasma frequency) [See Banks et al. (2017)]
 params.nuoveromegapemax = inf;         % Maximum ratio between electron/ion collision frequency and electron plasma frequency [See Banks et al. (2017)]. Set to inf if not desired !!!
 params.mu               = MODEL.sigma_e;   % sqrt(m_e/m_i)
 params.kpar             = 0.0;         % normalized parallel wave number to the major radius
 params.kperp            = kz;  % normalized perpendicular toroidal wave number to the soundLarmor radius. Note: If ions ==0 (e.g. EPW), kperp --> b
 params.kr               = kr;  % Radial component of perpendicular vector
 params.alphaD           = 0.0;         % (k*Debye length)^2
 params.Rn               = MODEL.eta_n; % Major Radius / Background density gradient length
 params.RTe              = MODEL.eta_T; % Major Radius * normalized kperp / Background electron temperature gradient length
 params.RTi              = MODEL.eta_T; % Major Radius * normalized kperp / Background ion temperature gradient length
 params.Rphi             = 0.0;         % Major Radius * normalized kperp / Background potentiel gradient length [presence of shear] - only for GK
 params.betae            = 1e-6;        % Electron Beta plasma.
 
 params.rhostar          = 1e-5;        % sound Larmor Radius/Major Radius ~ sqrt(Te)/(R_0*B).
-params.n0               = INITIAL.initback_moments;        % initial density perturbation
+params.n0               = INITIAL.init_background;        % initial density perturbation
 
 params.gradB            = MODEL.eta_B;         % Magnetic field gradient
 params.curvB            = MODEL.eta_B;         % Curvature of B
 params.trappB           = 0.0;         % Trapping term
 
 %% MOLI Options
 
 % Save data in results dir
 options.save    = 0;
 options.verbose = 0;
 options.dbg     = 0;
 
 options.DR        = 0;      % 1 -> Solve kinetic dispersion relation,
 options.KineticDR = 0;      % Solve kinetic dispersion relation (Landau integral) for the given theory
 
 % Compute the kinetic susceptibility for EPW only
 options.SPTBLTY = 0;
 
 options.nharm   = 1;        % Number of harmonics in disp. rel. 1 and 4
 
 wlim = 5.0;
 options.DRsolver.wr_min = -wlim;     % Minimum of real part.
 options.DRsolver.wr_max =  wlim;     % Maximum of real part.
 options.DRsolver.wi_min = -wlim;     % Minimum of imag part.
 options.DRsolver.wi_max =  wlim;     % Maximum of imag part.
 options.DRsolver.nw     =  300;      % Grid resolution
 
 % Disp. Rel. Options
 options.FLRmodel    = 0;   % 1 -> Truncated Laguerre, 0 -> Exact representation
 options.FluidLandau = 0;   % 1 -> Add Landau Fluid Closure to Fluid Dispersion Relation, 0 -> off
 options.deltaLandau = 0;   % 1 -> Hammet-Perkins closure on, 0 -> off
 
 % Fluid dispersion relation
 options.FluidDR = 0;	     % Solve annamaria's fluid equations
 options.Fluid.sITGmm = 0;
 
 % Define scan parameters
 options.fscan = 0;         % 1 -> peform scan over scan.list, 0-> off
 
 options.scan.list = {};% List of scan parameters. If empty, solve MOLI with params
 
 % Time-Evolution Problem [Solver==3] ...
 options.solver.TimeSolver.dt         = BASIC.dt;		% timestep of time evolution (R/c_s or 1/(k v/the) units)
 options.solver.TimeSolver.tmax       = BASIC.tmax;
 options.solver.TimeSolver.Trun       = BASIC.tmax;		  % total time to run time evolution
 options.solver.TimeSolver.t_fit_min  = 0.05;    % Phase-Mixing fit Lower time limit
 options.solver.TimeSolver.t_fit_max  = 8;       % Phase-Mixing fit Upper time limit
 options.solver.TimeSolver.en_fit_min = 0.15;    % Entropy Mode fit Lower time limit
 options.solver.TimeSolver.en_fit_max = 0.3;     % Entropy Mode fit Upper time limit
 options.solver.TimeSolver.movie      = 0;       % Display movie if 1, last frame otherwise
 options.solver.TimeSolver.save       = 0;       % 1 --> save during fscan, Warning: memory storage
 
 
 %% Run MOLI
 % Solve the MOLI
 [results,params,options] = MOLI_Control(params,options);
 
 %% Return to HeLaZ workspace
 cd ../../../HeLaZ/wk
diff --git a/matlab/setup.m b/matlab/setup.m
index e5ac181..d8c53dd 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -1,177 +1,177 @@
 %% ________________________________________________________________________
 SIMDIR = ['../results/',SIMID,'/'];
 % Grid parameters
 GRID.pmaxe = PMAXE;  % Electron Hermite moments
 GRID.jmaxe = JMAXE;  % Electron Laguerre moments
 GRID.pmaxi = PMAXI;  % Ion Hermite moments
 GRID.jmaxi = JMAXI;  % Ion Laguerre moments
 GRID.Nr    = N; % r grid resolution
 GRID.Lr    = L; % r length
 GRID.Nz    = N * (1-KREQ0) + KREQ0; % z ''
 GRID.Lz    = L * (1-KREQ0); % z ''
 GRID.kpar  = KPAR;
 
 % Model parameters
 MODEL.CO      = CO;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
 MODEL.CLOS    = CLOS;
 MODEL.NL_CLOS = NL_CLOS;
 if NON_LIN; MODEL.NON_LIN = '.true.'; else; MODEL.NON_LIN = '.false.';end;
 MODEL.mu      = MU;
 MODEL.mu_p    = MU_P;
 MODEL.mu_j    = MU_J;
 MODEL.nu      = NU; % hyper diffusive coefficient nu for HW
 % temperature ratio T_a/T_e
 MODEL.tau_e   = TAU;
 MODEL.tau_i   = TAU;
 % mass ratio sqrt(m_a/m_i)
 MODEL.sigma_e = 0.0233380;
 MODEL.sigma_i = 1.0;
 % charge q_a/e
 MODEL.q_e     =-1.0;
 MODEL.q_i     = 1.0;
 if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end;
 % gradients L_perp/L_x
 MODEL.eta_n   = ETAN;        % source term kappa for HW
 MODEL.eta_T   = ETAT;        % Temperature
 MODEL.eta_B   = ETAB;        % Magnetic
 MODEL.lambdaD = LAMBDAD;
 % if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KR0KH),'_A_',num2str(A0KH)]; end;
 % Time integration and intialization parameters
 TIME_INTEGRATION.numerical_scheme  = '''RK4''';
 if INIT_PHI; INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end;
-INITIAL.initback_moments  = 0.0e-5;
-INITIAL.initnoise_moments = NOISE0;
+INITIAL.init_background  = 0.0e-5;
+INITIAL.init_noiselvl = NOISE0;
 INITIAL.iseed             = 42;
 INITIAL.selfmat_file = '''null''';
 INITIAL.eimat_file = '''null''';
 INITIAL.iemat_file = '''null''';
 if (CO == -3) % Write matrice filename for Full Coulomb DK
     cmat_pmaxe = 25;
     cmat_jmaxe = 18;
     cmat_pmaxi = 25;
     cmat_jmaxi = 18;
     INITIAL.selfmat_file = ...
         ['''../../../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
         num2str(cmat_jmaxi),'_pamaxx_10.h5'''];
     INITIAL.eimat_file = ...
         ['''../../../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
         num2str(cmat_jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
     INITIAL.iemat_file = ...
         ['''../../../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
         num2str(cmat_jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
 elseif (CO == -2) % Write matrice filename for DK Sugama
     cmat_pmaxe = 10;
     cmat_jmaxe = 5;
     cmat_pmaxi = 10;
     cmat_jmaxi = 5;
     INITIAL.selfmat_file = ...
         ['''../../../iCa/self_Coll_GKE_0_GKI_0_ESELF_3_ISELF_3_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
         num2str(cmat_jmaxi),'_JE_12.h5'''];
     INITIAL.eimat_file = ...
         ['''../../../iCa/ei_Coll_GKE_0_GKI_0_ETEST_3_EBACK_3_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
         num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233.h5'''];
     INITIAL.iemat_file = ...
         ['''../../../iCa/ie_Coll_GKE_0_GKI_0_ITEST_3_IBACK_3_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
         num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233.h5'''];
 elseif (CO == 2) % Write matrice filename for Sugama GK
     cmat_pmaxe = 10;
     cmat_jmaxe = 5;
     cmat_pmaxi = 10;
     cmat_jmaxi = 5;
     INITIAL.selfmat_file = ...
         ['''../../../iCa/self_Coll_GKE_1_GKI_1_ESELF_3_ISELF_3_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
         num2str(cmat_jmaxi),'_JE_12_'''];
     INITIAL.eimat_file = ...
         ['''../../../iCa/ei_Coll_GKE_1_GKI_1_ETEST_3_EBACK_3_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
         num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233_'''];
     INITIAL.iemat_file = ...
         ['''../../../iCa/ie_Coll_GKE_1_GKI_1_ITEST_3_IBACK_3_Pmaxe_',num2str(cmat_pmaxe),...
         '_Jmaxe_',num2str(cmat_jmaxe),'_Pmaxi_',num2str(cmat_pmaxi),'_Jmaxi_',...
         num2str(cmat_jmaxi),'_JE_12_tau_1.0000_mu_0.0233_'''];
 elseif (CO == 3) % Full Coulomb GK
     disp('Warning, FCGK not implemented yet')
 elseif (CO == -1) % DGDK
     disp('Warning, DGDK not debugged')
 end
 
 % Naming and creating input file
 if    (CO == -3); CONAME = 'FCDK';
 elseif(CO == -2); CONAME = 'SGDK';
 elseif(CO == -1); CONAME = 'DGDK';
 elseif(CO ==  0); CONAME = 'LB';
 elseif(CO ==  1); CONAME = 'DGGK';
 elseif(CO ==  2); CONAME = 'SGGK';
 elseif(CO ==  3); CONAME = 'FCGK';
 end
 if    (CLOS == 0); CLOSNAME = 'Trunc.';
 elseif(CLOS == 1); CLOSNAME = 'Clos. 1';
 elseif(CLOS == 2); CLOSNAME = 'Clos. 2';
 end
 if (PMAXE == PMAXI) && (JMAXE == JMAXI)
     degngrad   = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE)];
 else
     degngrad   = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),...
         '_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)];
 end
 degngrad = [degngrad,'_eta_',num2str(ETAB/ETAN),'_nu_%0.0e_',...
         CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e'];
     
 degngrad   = sprintf(degngrad,[NU,MU]);
 if ~NON_LIN; degngrad = ['lin_',degngrad]; end
 resolution = [num2str(GRID.Nr),'x',num2str(GRID.Nz/2),'_'];
 gridname   = ['L_',num2str(L),'_'];
 if (exist('PREFIX','var') == 0); PREFIX = []; end;
 if (exist('SUFFIX','var') == 0); SUFFIX = []; end;
 PARAMS = [PREFIX,resolution,gridname,degngrad,SUFFIX];
 BASIC.RESDIR = [SIMDIR,PARAMS,'/'];
 BASIC.PARAMS = PARAMS;
 BASIC.SIMID  = SIMID;
 BASIC.nrun       = 1e8;
 BASIC.dt         = DT;
 BASIC.tmax       = TMAX;    %time normalized to 1/omega_pe
 BASIC.maxruntime = str2num(CLUSTER.TIME(1:2))*3600 ...
                    + str2num(CLUSTER.TIME(4:5))*60 ...
                    + str2num(CLUSTER.TIME(7:8));
 % Outputs parameters
 if RESTART; BASIC.RESTART = '.true.'; else; BASIC.RESTART = '.false.';end;
 OUTPUTS.nsave_0d = floor(1.0/SPS0D/DT);
 OUTPUTS.nsave_1d = -1;
 OUTPUTS.nsave_2d = floor(1.0/SPS2D/DT);
 OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT);
 OUTPUTS.nsave_cp = floor(1.0/SPSCP/DT);
 if W_DOUBLE; OUTPUTS.write_doubleprecision = '.true.'; else; OUTPUTS.write_doubleprecision = '.false.';end;
 if W_GAMMA;  OUTPUTS.write_gamma = '.true.'; else; OUTPUTS.write_gamma = '.false.';end;
 if W_PHI;    OUTPUTS.write_phi   = '.true.'; else; OUTPUTS.write_phi   = '.false.';end;
 if W_NA00;   OUTPUTS.write_Na00  = '.true.'; else; OUTPUTS.write_Na00  = '.false.';end;
 if W_NAPJ;   OUTPUTS.write_Napj  = '.true.'; else; OUTPUTS.write_Napj  = '.false.';end;
 if W_SAPJ;   OUTPUTS.write_Sapj  = '.true.'; else; OUTPUTS.write_Sapj  = '.false.';end;
 OUTPUTS.resfile0    = '''outputs''';
 OUTPUTS.rstfile0    = '''checkpoint''';
 OUTPUTS.job2load    = JOB2LOAD;
 %% Create directories
 if ~exist(SIMDIR, 'dir')
    mkdir(SIMDIR)
 end
 if ~exist(BASIC.RESDIR, 'dir')
 mkdir(BASIC.RESDIR)
 end
 %% Compile and WRITE input file
 INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
 nproc = 1;
 MAKE  = 'cd ..; make; cd wk';
 system(MAKE);
 %%
 disp(['Set up ',SIMID]);
 disp([resolution,gridname,degngrad]);
 if RESTART
 	disp(['- restarting from JOBNUM = ',num2str(JOB2LOAD)]); else
 	disp(['- starting from T = 0']);
 end
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 0bd29c1..aef016f 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -1,81 +1,83 @@
 function [INPUT] = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC)
 % Write the input script "fort.90" with desired parameters
 INPUT = 'fort.90';
 fid = fopen(INPUT,'wt');
 
 fprintf(fid,'&BASIC\n');
 fprintf(fid,['  nrun   = ', num2str(BASIC.nrun),'\n']);
 fprintf(fid,['  dt     = ', num2str(BASIC.dt),'\n']);
 fprintf(fid,['  tmax   = ', num2str(BASIC.tmax),'\n']);
 fprintf(fid,['  RESTART   = ', BASIC.RESTART,'\n']);
 fprintf(fid,['  maxruntime = ', num2str(BASIC.maxruntime),'\n']);
 fprintf(fid,'/\n');
 
 fprintf(fid,'&GRID\n');
 fprintf(fid,['  pmaxe =', num2str(GRID.pmaxe),'\n']);
 fprintf(fid,['  jmaxe = ', num2str(GRID.jmaxe),'\n']);
 fprintf(fid,['  pmaxi = ', num2str(GRID.pmaxi),'\n']);
 fprintf(fid,['  jmaxi = ', num2str(GRID.jmaxi),'\n']);
 fprintf(fid,['  Nr   = ', num2str(GRID.Nr),'\n']);
 fprintf(fid,['  Lr = ', num2str(GRID.Lr),'\n']);
 fprintf(fid,['  Nz   = ', num2str(GRID.Nz),'\n']);
 fprintf(fid,['  Lz = ', num2str(GRID.Lz),'\n']);
 fprintf(fid,['  kpar = ', num2str(GRID.kpar),'\n']);
 fprintf(fid,'/\n');
 
 fprintf(fid,'&OUTPUT_PAR\n');
 fprintf(fid,['  nsave_0d = ', num2str(OUTPUTS.nsave_0d),'\n']);
 fprintf(fid,['  nsave_1d = ', num2str(OUTPUTS.nsave_1d),'\n']);
 fprintf(fid,['  nsave_2d = ', num2str(OUTPUTS.nsave_2d),'\n']);
 fprintf(fid,['  nsave_5d = ', num2str(OUTPUTS.nsave_5d),'\n']);
 fprintf(fid,['  nsave_cp = ', num2str(OUTPUTS.nsave_cp),'\n']);
 fprintf(fid,['  write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']);
 fprintf(fid,['  write_gamma = ', OUTPUTS.write_gamma,'\n']);
 fprintf(fid,['  write_phi   = ', OUTPUTS.write_phi,'\n']);
 fprintf(fid,['  write_Na00 = ', OUTPUTS.write_Na00,'\n']);
 fprintf(fid,['  write_Napj = ', OUTPUTS.write_Napj,'\n']);
 fprintf(fid,['  write_Sapj = ', OUTPUTS.write_Sapj,'\n']);
 fprintf(fid,['  resfile0      = ', OUTPUTS.resfile0,'\n']);
 fprintf(fid,['  rstfile0      = ', OUTPUTS.rstfile0,'\n']);
 fprintf(fid,['  job2load      = ', num2str(OUTPUTS.job2load),'\n']);
 fprintf(fid,'/\n');
 
 fprintf(fid,'&MODEL_PAR\n');
 fprintf(fid,'  ! Collisionality\n');
 fprintf(fid,['  CO      = ', num2str(MODEL.CO),'\n']);
 fprintf(fid,['  CLOS    = ', num2str(MODEL.CLOS),'\n']);
 fprintf(fid,['  NL_CLOS = ', num2str(MODEL.NL_CLOS),'\n']);
 fprintf(fid,['  NON_LIN = ', MODEL.NON_LIN,'\n']);
+fprintf(fid,['  ZF_DAMP = ', MODEL.ZF_DAMP,'\n']);
 fprintf(fid,['  mu      = ', num2str(MODEL.mu),'\n']);
 fprintf(fid,['  mu_p    = ', num2str(MODEL.mu_p),'\n']);
 fprintf(fid,['  mu_j    = ', num2str(MODEL.mu_j),'\n']);
 fprintf(fid,['  nu      = ', num2str(MODEL.nu),'\n']);
 fprintf(fid,['  tau_e   = ', num2str(MODEL.tau_e),'\n']);
 fprintf(fid,['  tau_i   = ', num2str(MODEL.tau_i),'\n']);
 fprintf(fid,['  sigma_e = ', num2str(MODEL.sigma_e),'\n']);
 fprintf(fid,['  sigma_i = ', num2str(MODEL.sigma_i),'\n']);
 fprintf(fid,['  q_e     = ', num2str(MODEL.q_e),'\n']);
 fprintf(fid,['  q_i     = ', num2str(MODEL.q_i),'\n']);
 fprintf(fid,['  eta_n   = ', num2str(MODEL.eta_n),'\n']);
 fprintf(fid,['  eta_T   = ', num2str(MODEL.eta_T),'\n']);
 fprintf(fid,['  eta_B   = ', num2str(MODEL.eta_B),'\n']);
 fprintf(fid,['  lambdaD = ', num2str(MODEL.lambdaD),'\n']);
 fprintf(fid,'/\n');
 
 fprintf(fid,'&INITIAL_CON\n');
 fprintf(fid,['  INIT_NOISY_PHI    =', INITIAL.init_noisy_phi,'\n']);
-fprintf(fid,['  initback_moments  =', num2str(INITIAL.initback_moments),'\n']);
-fprintf(fid,['  initnoise_moments =', num2str(INITIAL.initnoise_moments),'\n']);
+fprintf(fid,['  INIT_ZF_PHI       =', INITIAL.init_noisy_phi,'\n']);
+fprintf(fid,['  init_background  =', num2str(INITIAL.init_background),'\n']);
+fprintf(fid,['  init_noiselvl =', num2str(INITIAL.init_noiselvl),'\n']);
 fprintf(fid,['  iseed             =', num2str(INITIAL.iseed),'\n']);
 fprintf(fid,['  selfmat_file      =', INITIAL.selfmat_file,'\n']);
 fprintf(fid,['  eimat_file        =', INITIAL.eimat_file,'\n']);
 fprintf(fid,['  iemat_file        =', INITIAL.iemat_file,'\n']);
 fprintf(fid,'/\n');
 
 fprintf(fid,'&TIME_INTEGRATION_PAR\n');
 fprintf(fid,['  numerical_scheme=', TIME_INTEGRATION.numerical_scheme,'\n']);
 fprintf(fid,'/');
 
 fclose(fid);
 system(['cp fort.90 ',BASIC.RESDIR,'/.']);
 end
diff --git a/src/inital.F90 b/src/inital.F90
index a4502cb..680f4f8 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -1,347 +1,347 @@
 !******************************************************************************!
 !!!!!! initialize the moments and load/build coeff tables
 !******************************************************************************!
 SUBROUTINE inital
 
   USE basic
   USE model, ONLY : CO, NON_LIN
   USE initial_par
   USE prec_const
   USE coeff
   USE time_integration
   USE array, ONLY : Sepj,Sipj
   USE collision
   USE closure
   USE ghosts
   USE restarts
 
   implicit none
 
   CALL set_updatetlevel(1)
 
   IF (my_id .EQ. 0) WRITE(*,*) 'Evaluate kernels'
   CALL evaluate_kernels
 
   !!!!!! Set the moments arrays Nepj, Nipj and phi!!!!!!
   IF ( RESTART ) THEN
     IF (my_id .EQ. 0) WRITE(*,*) 'Load moments'
     CALL load_moments ! get N_0
 
     IF (my_id .EQ. 0) WRITE(*,*) 'Init phi with Poisson'
     CALL poisson ! compute phi_0=phi(N_0)
 
   ELSE
     IF (INIT_NOISY_PHI) THEN
       IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
       CALL init_phi ! init noisy phi_0, N_0 = 0
     ELSE
       IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy moments and ghosts'
       CALL init_moments ! init noisy N_0
       IF (my_id .EQ. 0) WRITE(*,*) 'Init phi with Poisson'
       CALL poisson ! get phi_0 = phi(N_0)
     ENDIF
 
   ENDIF
 
   IF (my_id .EQ. 0) WRITE(*,*) 'Apply closure'
   CALL apply_closure_model
 
   IF (my_id .EQ. 0) WRITE(*,*) 'Ghosts communication'
   CALL update_ghosts
 
   !!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!!
   IF ( NON_LIN ) THEN;
     IF (my_id .EQ. 0) WRITE(*,*) 'Building Dnjs table'
     CALL build_dnjs_table
 
     IF (my_id .EQ. 0) WRITE(*,*) 'Init Sapj'
     CALL compute_Sapj ! compute S_0 = S(phi_0,N_0)
   ENDIF
 
   !!!!!! Load the COSOlver collision operator coefficients !!!!!!
   IF (ABS(CO) .GT. 1) THEN
     CALL load_COSOlver_mat
     ! Compute collision
     CALL compute_TColl ! compute C_0 = C(N_0)
   ENDIF
 
 END SUBROUTINE inital
 !******************************************************************************!
 
 !******************************************************************************!
 !!!!!!! Initialize the moments randomly
 !******************************************************************************!
 SUBROUTINE init_moments
   USE basic
   USE grid
   USE fields
   USE prec_const
   USE utility, ONLY: checkfield
   USE initial_par
   USE model, ONLY : NON_LIN
   IMPLICIT NONE
 
   REAL(dp) :: noise
   REAL(dp) :: kr, kz, sigma, gain, kz_shift
   INTEGER, DIMENSION(12) :: iseedarr
 
   ! Seed random number generator
   iseedarr(:)=iseed
   CALL RANDOM_SEED(PUT=iseedarr+my_id)
 
     !**** Broad noise initialization *******************************************
     DO ip=ips_e,ipe_e
       DO ij=ijs_e,ije_e
 
         DO ikr=ikrs,ikre
           DO ikz=ikzs,ikze
             CALL RANDOM_NUMBER(noise)
-            moments_e( ip,ij, ikr,ikz, :) = (initback_moments + initnoise_moments*(noise-0.5_dp))
+            moments_e( ip,ij, ikr,ikz, :) = (init_background + init_noiselvl*(noise-0.5_dp))
           END DO
         END DO
 
         IF ( contains_kr0 ) THEN
           DO ikz=2,Nkz/2 !symmetry at kr = 0
             moments_e( ip,ij,ikr_0,ikz, :) = moments_e( ip,ij,ikr_0,Nkz+2-ikz, :)
           END DO
         ENDIF
 
       END DO
     END DO
 
     DO ip=ips_i,ipe_i
       DO ij=ijs_i,ije_i
 
         DO ikr=ikrs,ikre
           DO ikz=ikzs,ikze
             CALL RANDOM_NUMBER(noise)
-            moments_i( ip,ij, ikr,ikz, :) = (initback_moments + initnoise_moments*(noise-0.5_dp))
+            moments_i( ip,ij, ikr,ikz, :) = (init_background + init_noiselvl*(noise-0.5_dp))
           END DO
         END DO
 
         IF ( contains_kr0 ) THEN
           DO ikz=2,Nkz/2 !symmetry at kr = 0
             moments_i( ip,ij,ikr_0,ikz, :) = moments_i( ip,ij,ikr_0,Nkz+2-ikz, :)
           END DO
         ENDIF
 
       END DO
     END DO
 
     ! Putting to zero modes that are not in the 2/3 Orszag rule
     IF (NON_LIN) THEN
       DO ip=ips_e,ipe_e
       DO ij=ijs_e,ije_e
       DO ikr=ikrs,ikre
       DO ikz=ikzs,ikze
         moments_e( ip,ij,ikr,ikz, :) = moments_e( ip,ij,ikr,ikz, :)*AA_r(ikr)*AA_z(ikz)
       ENDDO
       ENDDO
       ENDDO
       ENDDO
       DO ip=ips_i,ipe_i
       DO ij=ijs_i,ije_i
       DO ikr=ikrs,ikre
       DO ikz=ikzs,ikze
         moments_i( ip,ij,ikr,ikz, :) = moments_i( ip,ij,ikr,ikz, :)*AA_r(ikr)*AA_z(ikz)
       ENDDO
       ENDDO
       ENDDO
       ENDDO
     ENDIF
 END SUBROUTINE init_moments
 !******************************************************************************!
 
 
 !******************************************************************************!
 !!!!!!! Initialize a noisy ES potential and cancel the moments
 !******************************************************************************!
 SUBROUTINE init_phi
   USE basic
   USE grid
   USE fields
   USE prec_const
   USE initial_par
   IMPLICIT NONE
 
   REAL(dp) :: noise
   REAL(dp) :: kr, kz, sigma, gain, kz_shift
   INTEGER, DIMENSION(12) :: iseedarr
 
   IF (INIT_NOISY_PHI) THEN
 
     ! Seed random number generator
     iseedarr(:)=iseed
     CALL RANDOM_SEED(PUT=iseedarr+my_id)
 
       !**** noise initialization *******************************************
 
       DO ikr=ikrs,ikre
         DO ikz=ikzs,ikze
           CALL RANDOM_NUMBER(noise)
-          phi(ikr,ikz) = (initback_moments + initnoise_moments*(noise-0.5_dp))*AA_r(ikr)*AA_z(ikz)
+          phi(ikr,ikz) = (init_background + init_noiselvl*(noise-0.5_dp))*AA_r(ikr)*AA_z(ikz)
         END DO
       END DO
 
       !symmetry at kr = 0 to keep real inverse transform
       IF ( contains_kr0 ) THEN
         DO ikz=2,Nkz/2
           phi(ikr_0,ikz) = phi(ikr_0,Nkz+2-ikz)
         END DO
         phi(ikr_0,Nz/2) = REAL(phi(ikr_0,Nz/2)) !origin must be real
       ENDIF
 
       !**** Cancel previous moments initialization
       DO ip=ips_e,ipe_e
         DO ij=ijs_e,ije_e
           DO ikr=ikrs,ikre
             DO ikz=ikzs,ikze
               moments_e( ip,ij,ikr,ikz, :) = 0._dp
             ENDDO
           ENDDO
         ENDDO
       ENDDO
       DO ip=ips_i,ipe_i
         DO ij=ijs_i,ije_i
           DO ikr=ikrs,ikre
             DO ikz=ikzs,ikze
               moments_i( ip,ij,ikr,ikz, :) = 0._dp
             ENDDO
           ENDDO
         ENDDO
       ENDDO
 
     ELSE ! we compute phi from noisy moments and poisson
 
       CALL poisson
     ENDIF
 
 END SUBROUTINE init_phi
 !******************************************************************************!
 
 !******************************************************************************!
 !!!!!!! Build the Laguerre-Laguerre coupling coefficient table for nonlin
 !******************************************************************************!
 SUBROUTINE build_dnjs_table
   USE basic
   USE array, Only : dnjs
   USE grid, Only : jmaxe, jmaxi
   USE coeff
   IMPLICIT NONE
 
   INTEGER :: in, ij, is, J
   INTEGER :: n_, j_, s_
 
   J = max(jmaxe,jmaxi)
 
   DO in = 1,J+1 ! Nested dependent loops to make benefit from dnjs symmetry
     n_ = in - 1
     DO ij = in,J+1
       j_ = ij - 1
       DO is = ij,J+1
         s_ = is - 1
 
         dnjs(in,ij,is) = TO_DP(ALL2L(n_,j_,s_,0))
         ! By symmetry
         dnjs(in,is,ij) = dnjs(in,ij,is)
         dnjs(ij,in,is) = dnjs(in,ij,is)
         dnjs(ij,is,in) = dnjs(in,ij,is)
         dnjs(is,ij,in) = dnjs(in,ij,is)
         dnjs(is,in,ij) = dnjs(in,ij,is)
       ENDDO
     ENDDO
   ENDDO
 END SUBROUTINE build_dnjs_table
 !******************************************************************************!
 
 !******************************************************************************!
 !!!!!!! Evaluate the kernels once for all
 !******************************************************************************!
 SUBROUTINE evaluate_kernels
   USE basic
   USE array, Only : kernel_e, kernel_i
   USE grid
   use model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, CLOS
   IMPLICIT NONE
 
   REAL(dp)    :: factj, j_dp, j_int
   REAL(dp)    :: sigmae2_taue_o2, sigmai2_taui_o2
   REAL(dp)    :: be_2, bi_2, alphaD
   REAL(dp)    :: kr, kz, kperp2
 
   !!!!! Electron kernels !!!!!
   !Precompute species dependant factors
   sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument
 
   factj = 1.0 ! Start of the recursive factorial
   DO ij = 1, jmaxe+1
     j_int = jarray_e(ij)
     j_dp = REAL(j_int,dp) ! REAL of degree
 
     ! Recursive factorial
     IF (j_dp .GT. 0) THEN
       factj = factj * j_dp
     ELSE
       factj = 1._dp
     ENDIF
 
     DO ikr = ikrs,ikre
       kr     = krarray(ikr)
       DO ikz = ikzs,ikze
         kz    = kzarray(ikz)
 
         be_2  =  (kr**2 + kz**2) * sigmae2_taue_o2
 
         kernel_e(ij, ikr, ikz) = be_2**j_int * exp(-be_2)/factj
 
       ENDDO
     ENDDO
   ENDDO
   ! Kernels closure
   DO ikr = ikrs,ikre
     kr     = krarray(ikr)
     DO ikz = ikzs,ikze
       kz    = kzarray(ikz)
       be_2  =  (kr**2 + kz**2) * sigmae2_taue_o2
       kernel_e(ijeg_e,ikr,ikz) = be_2/(real(ijeg_e,dp))*kernel_e(ije_e,ikr,ikz)
     ENDDO
   ENDDO
 
   !!!!! Ion kernels !!!!!
   sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp ! (b_a/2)^2 = (kperp sqrt(2 tau_a) sigma_a/2)^2
 
   factj = 1.0 ! Start of the recursive factorial
   DO ij = 1, jmaxi+1
     j_int = jarray_e(ij)
     j_dp = REAL(j_int,dp) ! REAL of degree
 
     ! Recursive factorial
     IF (j_dp .GT. 0) THEN
       factj = factj * j_dp
     ELSE
       factj = 1._dp
     ENDIF
 
     DO ikr = ikrs,ikre
       kr     = krarray(ikr)
       DO ikz = ikzs,ikze
         kz    = kzarray(ikz)
 
         bi_2  =  (kr**2 + kz**2) * sigmai2_taui_o2
 
         kernel_i(ij, ikr, ikz) = bi_2**j_int * exp(-bi_2)/factj
 
       ENDDO
     ENDDO
   ENDDO
   ! Kernels closure
   DO ikr = ikrs,ikre
     kr     = krarray(ikr)
     DO ikz = ikzs,ikze
       kz    = kzarray(ikz)
       bi_2  =  (kr**2 + kz**2) * sigmai2_taui_o2
       kernel_i(ijeg_i,ikr,ikz) = bi_2/(real(ijeg_i,dp))*kernel_e(ije_i,ikr,ikz)
     ENDDO
   ENDDO
 END SUBROUTINE evaluate_kernels
 !******************************************************************************!
diff --git a/src/initial_par_mod.F90 b/src/initial_par_mod.F90
index fea75a0..56fa6a2 100644
--- a/src/initial_par_mod.F90
+++ b/src/initial_par_mod.F90
@@ -1,72 +1,75 @@
 MODULE initial_par
   !   Module for initial parameters
 
   USE prec_const
   IMPLICIT NONE
   PRIVATE
 
   ! Initial background level
-  REAL(dp), PUBLIC, PROTECTED :: initback_moments=0._dp
-  ! Initial background level
+  REAL(dp), PUBLIC, PROTECTED :: init_background=0._dp
+  ! Initialization through a noisy phi
   LOGICAL,  PUBLIC, PROTECTED :: INIT_NOISY_PHI = .false.
+  ! Initialization through a zonal flow phi
+  LOGICAL,  PUBLIC, PROTECTED :: INIT_ZF_PHI    = .false.
   ! Initial background noise amplitude
-  REAL(dp), PUBLIC, PROTECTED :: initnoise_moments=1E-6_dp
+  REAL(dp), PUBLIC, PROTECTED :: init_noiselvl=1E-6_dp
   ! Initialization for random number generator
   INTEGER,  PUBLIC, PROTECTED :: iseed=42
 
   ! Parameters of initial smooth sine profiles
   REAL(dp), PUBLIC, PROTECTED ::  init_nb_oscil_density=2._dp ! Number of oscillations
   REAL(dp), PUBLIC, PROTECTED ::  init_nb_oscil_temp=2._dp
   REAL(dp), PUBLIC, PROTECTED ::  init_ampli_density=0.1_dp ! Oscillation amplitude
   REAL(dp), PUBLIC, PROTECTED ::  init_ampli_temp=0.1_dp
 
   CHARACTER(len=128), PUBLIC :: selfmat_file  ! COSOlver matrix file names
   CHARACTER(len=128), PUBLIC :: iemat_file  ! COSOlver matrix file names
   CHARACTER(len=128), PUBLIC :: eimat_file  ! COSOlver matrix file names
 
   PUBLIC :: initial_outputinputs, initial_readinputs
 
 CONTAINS
 
 
   SUBROUTINE initial_readinputs
     ! Read the input parameters
 
     USE basic, ONLY : lu_in, RESTART
     USE prec_const
     IMPLICIT NONE
 
     NAMELIST /INITIAL_CON/ INIT_NOISY_PHI
-    NAMELIST /INITIAL_CON/ initback_moments
-    NAMELIST /INITIAL_CON/ initnoise_moments
+    NAMELIST /INITIAL_CON/ INIT_ZF_PHI
+    NAMELIST /INITIAL_CON/ init_background
+    NAMELIST /INITIAL_CON/ init_noiselvl
     NAMELIST /INITIAL_CON/ iseed
     NAMELIST /INITIAL_CON/ selfmat_file
     NAMELIST /INITIAL_CON/ iemat_file
     NAMELIST /INITIAL_CON/ eimat_file
 
     READ(lu_in,initial_con)
     !WRITE(*,initial_con)
 
   END SUBROUTINE initial_readinputs
 
 
   SUBROUTINE initial_outputinputs(fidres, str)
     ! Write the input parameters to the results_xx.h5 file
 
     USE futils, ONLY: attach
     USE prec_const
     IMPLICIT NONE
     INTEGER, INTENT(in) :: fidres
     CHARACTER(len=256), INTENT(in) :: str
 
     CALL attach(fidres, TRIM(str), "INIT_NOISY_PHI", INIT_NOISY_PHI)
 
-    CALL attach(fidres, TRIM(str), "initback_moments", initback_moments)
+    CALL attach(fidres, TRIM(str), "init_background", init_background)
 
-    CALL attach(fidres, TRIM(str), "initnoise_moments", initnoise_moments)
+    CALL attach(fidres, TRIM(str), "init_noiselvl", init_noiselvl)
 
     CALL attach(fidres, TRIM(str), "iseed", iseed)
 
   END SUBROUTINE initial_outputinputs
 
 END MODULE initial_par