From 122a8cf5893866aa72ebefcf6e2b13e4bf2a305f Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Thu, 4 Apr 2024 09:18:26 +0200
Subject: [PATCH] Save matlab scripts add the new magnetic tuners have a
 consistent writing with capital k

---
 matlab/setup.m                         |  12 +-
 matlab/write_fort90.m                  |   7 +-
 wk/lin_run_script.m                    |  35 +----
 wk/lin_scan_script.m                   |  17 ++-
 wk/lin_scan_script_P_J.m               | 203 +++++++++++++++++++++++++
 wk/lin_scan_script_P_ky.m              | 201 ++++++++++++++++++++++++
 wk/load_metadata_scan.m                |  40 ++---
 wk/parameters/lin_AUG.m                |   7 +-
 wk/parameters/lin_DIIID_AUSTIN.m       |   7 +-
 wk/parameters/lin_DIIID_LM_rho90.m     |   7 +-
 wk/parameters/lin_DIIID_LM_rho95.m     |   7 +-
 wk/parameters/lin_DIIID_LM_rho95_HEL.m |   7 +-
 wk/parameters/lin_DIIID_data.m         |   7 +-
 wk/parameters/lin_DTT_HM_rho85.m       |   7 +-
 wk/parameters/lin_DTT_HM_rho98.m       |   7 +-
 wk/parameters/lin_ETG_adiab_i.m        |   7 +-
 wk/parameters/lin_Entropy.m            |   7 +-
 wk/parameters/lin_GASTD.m              |   7 +-
 wk/parameters/lin_ITG.m                |   7 +-
 wk/parameters/lin_Ivanov.m             |   7 +-
 wk/parameters/lin_JET_rho97.m          |   7 +-
 wk/parameters/lin_KBM.m                |   7 +-
 wk/parameters/lin_RHT.m                |   7 +-
 wk/parameters/lin_STEP_EC_HD_psi49.m   |   7 +-
 wk/parameters/lin_STEP_EC_HD_psi71.m   |   7 +-
 25 files changed, 548 insertions(+), 93 deletions(-)
 create mode 100644 wk/lin_scan_script_P_J.m
 create mode 100644 wk/lin_scan_script_P_ky.m

diff --git a/matlab/setup.m b/matlab/setup.m
index 7fd5e06d..ba45da52 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -61,8 +61,16 @@ MODEL.K_Ni    = K_Ni;
 MODEL.K_Ne    = K_Ne;
 MODEL.K_Ti    = K_Ti;    
 MODEL.K_Te    = K_Te;    
-MODEL.k_gB   = k_gB;      % Magnetic gradient
-MODEL.k_cB   = k_cB;      % Magnetic curvature
+MODEL.K_gB   = K_gB;      % artificial magnetic gradient  tuner
+MODEL.K_cB   = K_cB;      % artificial magnetic curvature tuner
+try
+    K_mB; K_tB; K_ldB;
+catch
+    K_mB=1; K_tB=1; K_ldB=1;
+end
+MODEL.K_mB   = K_mB;      % artificial mirror force   tuner
+MODEL.K_tB   = K_tB;      % artificial trapping term  tuner
+MODEL.K_ldB  = K_ldB;     % artificial Landau damping tuner
 MODEL.lambdaD = LAMBDAD;
 % CLOSURE parameters
 CLOSURE.hierarchy_closure = ['''',HRCY_CLOS,''''];
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index f0bda4cd..0d81ae12 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -68,8 +68,11 @@ fprintf(fid,['  HYP_V   = ', MODEL.HYP_V,'\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,['  k_gB    = ', num2str(MODEL.k_gB),'\n']);
-fprintf(fid,['  k_cB    = ', num2str(MODEL.k_cB),'\n']);
+fprintf(fid,['  k_gB    = ', num2str(MODEL.K_gB),'\n']);
+fprintf(fid,['  k_cB    = ', num2str(MODEL.K_cB),'\n']);
+fprintf(fid,['  k_mB    = ', num2str(MODEL.K_mB),'\n']);
+fprintf(fid,['  k_tB    = ', num2str(MODEL.K_tB),'\n']);
+fprintf(fid,['  k_ldB   = ', num2str(MODEL.K_ldB),'\n']);
 fprintf(fid,['  lambdaD = ', num2str(MODEL.lambdaD),'\n']);
 fprintf(fid,['  beta    = ', num2str(MODEL.beta),'\n']);
 fprintf(fid,['  ExBrate = ', num2str(MODEL.ExBrate),'\n']);
diff --git a/wk/lin_run_script.m b/wk/lin_run_script.m
index 09ace133..79127a73 100644
--- a/wk/lin_run_script.m
+++ b/wk/lin_run_script.m
@@ -48,16 +48,18 @@ end
 % SIMID = ['rho_scan_DIIID_AUSTIN_2019/3x2x192x96x32/rho',num2str(rho)];
 %% Change parameters
 % GEOMETRY = 's-alpha';
-PMAX  = 2; JMAX = 1;
-DELTA =0.0; 
+PMAX  = 4; JMAX = 2;
+DELTA = 0.2; 
+% K_tB = 0; K_mB = 0; K_ldB = 0;
 % K_Ni = 0; K_Ne = 0;
-% DELTA = 0.0; 
+% K_Ti = 0; TAU = 1/3;
+DELTA = 0.0; 
 % DELTA = 0.2; 
 S_DELTA = DELTA/2;
-LY   = 2*pi/0.05;
+LY   = 2*pi/0.75;
 TMAX = 40;
 NY   = 2;
-DT   = 0.01;
+DT   = 0.0025;
 % TAU  = 1; NU = 0.05;
 % TAU = 1e-3; K_Ti = K_Ti/2/TAU; NU = 3*NU/8/TAU; ADIAB_E = 1; NA = 1;
 % MU_X = 1; MU_Y = 1;
@@ -114,29 +116,6 @@ gkxky = real(wkykx(2:end,1:data.grids.Nx/2))';
 gkxky(isnan(gkxky)) =0;
 gkxky(isinf(gkxky)) =0;
 figure; plot(ky,gkxky(1,:));
-% gkxky(gkxky<0)      =0;
-% % gkxky = imgaussfilt(gkxky,1);
-% %
-% wkxky = imag(wkykx(2:end,1:data.grids.Nx/2))';
-% wkxky(isnan(wkxky)) =0;
-% wkxky(isinf(wkxky)) =0;
-% % wkxky(wkxky<0)      =0;
-% % wkxky = imgaussfilt(wkxky,1.5);
-% %
-% figure; 
-% subplot(121)
-%     contourf(kx,ky,gkxky',10)
-%     % clim(0.5*[0 1]); 
-%     % colormap(bluewhitered); colorbar;
-%     xlim([0.025 1]);
-%     xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$')
-% subplot(122)
-%     contourf(kx,ky,wkxky',10)
-%     % clim(1*[0 1]); 
-%     % colormap(bluewhitered); colorbar 
-%     xlim([0.025 1]);
-%     xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$')
-% % save_figure(data,fig,'.png')
 end
 
 if (0 && NZ>4)
diff --git a/wk/lin_scan_script.m b/wk/lin_scan_script.m
index f95c86bb..cad3ccb8 100644
--- a/wk/lin_scan_script.m
+++ b/wk/lin_scan_script.m
@@ -39,20 +39,22 @@ run lin_DIIID_LM_rho95
 % NU   = 1;
 % TAU  = 1;
 NY   = 2;
+DELTA = 0.0; TRIANG = '';
+S_DELTA = DELTA/2;
 % EXBRATE = 0;
 % S_DELTA = min(2.0,S_DELTA);
 % SIGMA_E  = 0.023;
 % NEXC = 0;
 LX   = 120;
 %% Scan parameters
-SIMID = [SIMID,'_scan'];
-P_a   = [2 4]; J_a = [1 1];
+SIMID = [SIMID,TRIANG,'_scan'];
+P_a   = [2]; J_a = [1];
 % P_a   = 2;
 % ky_a  = [0.01 0.02 0.05 0.1  0.2  0.5  1.0  2.0  5.0  10.0];
-ky_a  = [0.025:0.025:1.1];
+ky_a  = linspace(0.1,1.1,25);
 % ky_a  = 4.0;
 % dt_a  = logspace(-2,-3,numel(ky_a));
-dt_a  = linspace(0.01,0.01,numel(ky_a));
+dt_a  = linspace(0.001,0.001,numel(ky_a));
 CO    = 'DG';
 %% Scan loop
 % arrays for the result
@@ -68,8 +70,8 @@ for PMAX = P_a
         LY   = 2*pi/ky;
         DT   = dt_a(i);%1e-3;%/(1+log(ky/0.05));%min(1e-2,1e-3/ky);
         TMAX = DT*10000;%2;%min(10,1.5/ky);
-        DTSAVE0D = 100*DT;
-        DTSAVE3D =  10*DT;
+        DTSAVE0D =  1.0;
+        DTSAVE3D =  0.5;
         %% RUN
         setup
         % naming
@@ -142,7 +144,8 @@ if(numel(ky_a)>1 || numel(P_a)>1)
                 '_ky_',kymin,'_',kymax,...
                 '_P_',pmin,'_',pmax,...
                 '_kN_',num2str(K_Ni),...
-                '_',CONAME,'_',num2str(NU),'_be_',num2str(BETA),'.mat'];
+                '_',CONAME,'_',num2str(NU),'_be_',num2str(BETA),...,
+                '_d_',num2str(DELTA),'.mat'];
     metadata.name   = filename;
     metadata.kymin  = ky;
     metadata.title  = ['$\nu_{',CONAME,'}=$',num2str(NU),'$\kappa_T=$',num2str(K_Ti),', $\kappa_N=$',num2str(K_Ni)];
diff --git a/wk/lin_scan_script_P_J.m b/wk/lin_scan_script_P_J.m
new file mode 100644
index 00000000..a94fd90f
--- /dev/null
+++ b/wk/lin_scan_script_P_J.m
@@ -0,0 +1,203 @@
+%% 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
+wkdir = pwd;
+gyacomodir = wkdir(1:end-2);
+mpirun     = 'mpirun';
+% mpirun     = '/opt/homebrew/bin/mpirun'; % for macos
+addpath(genpath([gyacomodir,'matlab']))         % Add matlab folder
+addpath(genpath([gyacomodir,'matlab/plot']))    % Add plot folder
+addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute folder
+addpath(genpath([gyacomodir,'matlab/load']))    % Add load folder
+addpath(genpath([gyacomodir,'wk/parameters']))  % Add parameters folder
+
+%% Setup run or load an executable
+RUN     = 1; % To run or just to load
+RERUN   = 0; % rerun if the  data does not exist
+default_plots_options
+% EXECNAME = 'gyacomo23_sp'; % single precision
+EXECNAME = 'gyacomo23_dp'; % double precision
+
+%% Setup parameters
+% run lin_DTT_AB_rho85
+% run lin_DTT_AB_rho98
+% run lin_JET_rho97
+% run lin_Entropy
+% run lin_ITG
+% run lin_RHT
+% rho  = 0.95; TRIANG = 'PT'; READPROF = 1; 
+% prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/'];
+% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson/',TRIANG,'/'];
+% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/'];
+% run lin_DIIID_data
+run lin_DIIID_LM_rho95
+
+%% Change parameters
+% NU   = 1;
+% TAU  = 1;
+NY   = 2;
+DELTA = 0.0; TRIANG = '';
+S_DELTA = DELTA/2;
+% EXBRATE = 0;
+% S_DELTA = min(2.0,S_DELTA);
+% SIGMA_E  = 0.023;
+% NEXC = 0;
+LX   = 120;
+%% Scan parameters
+SIMID = [SIMID,TRIANG,'_scan'];
+P_a   = [2 4 6 8]; 
+J_a   = [1 2 3 4];
+ky    = 0.3;
+DT    = 0.001; TMAX = 40;
+CO    = 'DG';
+%% Scan loop
+% arrays for the result
+g_ky = zeros(numel(P_a),numel(J_a));
+g_std= g_ky*0;
+w_ky = g_ky*0;
+w_std= g_ky*0;
+j = 1;
+for PMAX = P_a
+    i = 1;
+    for JMAX = J_a
+        LY   = 2*pi/ky;
+        DTSAVE0D =  1.0;
+        DTSAVE3D =  0.5;
+        %% 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))
+            MVIN =['cd ',LOCALDIR,';'];
+            % RUNG  =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
+            RUNG  =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
+            % RUNG  =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;'];
+            % RUNG  =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
+            % RUNG = ['./../../../bin/gyacomo23_sp 0;'];
+            MVOUT=['cd ',wkdir,';'];
+            system([MVIN,RUNG,MVOUT]);
+        end
+        data_    = compile_results_low_mem(data_,LOCALDIR,00,00);
+        [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+        if numel(data_.Ts0D)>10
+            % 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');
+            options.NORMALIZED = 0; 
+            options.TIME   = data_.Ts3D;
+             % Time window to measure the growth of kx/ky modes
+            options.KY_TW  = [0.7 1.0]*data_.Ts3D(end);
+            options.KX_TW  = [0.7 1.0]*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.GOK2   = 0; % plot gamma/k^2
+            options.fftz.flag = 0; % Set fftz.flag option to 0
+            options.FIELD = 'phi';
+            options.SHOWFIG = 0;
+            [fig, wkykx, ekykx] = mode_growth_meter(data_,options);
+            % [wkykx,ekykx] = compute_growth_rates(data_.PHI(:,:,:,it1:it2),data_.Ts3D(it1:it2));
+            g_ky (i,j) = real(wkykx(2,1));
+            g_std(i,j) = real(ekykx(2,1));
+            w_ky (i,j) = imag(wkykx(2,1));
+            w_std(i,j) = imag(ekykx(2,1));
+            [gmax, ikmax] = max(g_ky(i,j));
+
+            msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
+        end
+        i = i + 1;
+    end
+    j = j + 1;
+end
+
+%% take max growth rate among z coordinate
+y_ = g_ky + 1i*w_ky; 
+e_ = g_std+ 1i*w_std;
+
+%% Save scan results (gamma)
+metadata = {};
+if(numel(ky_a)>1 || numel(P_a)>1)
+    pmin  = num2str(min(P_a));   pmax = num2str(max(P_a));
+    jmin  = num2str(min(J_a));   jmax = num2str(max(J_a));
+    kymin = num2str(min(ky_a));  kymax= num2str(max(ky_a));
+    filename = [num2str(NX),'x',num2str(NZ),...
+                '_P_',pmin,'_',pmax,...
+                '_J_',jmin,'_',jmax,...
+                '_ky_',ky,...
+                '_kN_',num2str(K_Ni),...
+                '_',CONAME,'_',num2str(NU),'_be_',num2str(BETA),...,
+                '_d_',num2str(DELTA),'.mat'];
+    metadata.name   = filename;
+    metadata.kymin  = ky;
+    metadata.title  = [...
+        '$k_y=$',num2str(ky),...
+        ', $\nu_{',CONAME,'}=$',num2str(NU),...
+        ', $\kappa_T=$',num2str(K_Ti),...
+        ', $\kappa_N=$',num2str(K_Ni)];
+    metadata.par    = [num2str(NX),'x1x',num2str(NZ)];
+    metadata.nscan  = 2;
+    metadata.s1name = '$P$';
+    metadata.s1     = P_a;
+    metadata.s2name = '$J$';
+    metadata.s2     = J_a;
+    metadata.dname  = '$\gamma c_s/R$';
+    metadata.data   = y_;
+    metadata.err    = e_;
+    save([SIMDIR,filename],'-struct','metadata');
+    disp(['saved in ',SIMDIR,filename]);
+if 1
+    gamma = real(metadata.data); g_err = real(metadata.err);
+    omega = imag(metadata.data); w_err = imag(metadata.err);
+    gamma = gamma.*(gamma>0.025);
+    figure
+    colors_ = jet(numel(metadata.s2));
+    subplot(121)
+    for i = 1:numel(metadata.s2)
+        errorbar(metadata.s1,gamma(:,i),0*g_err(:,i),'s-',...
+        'LineWidth',2.0,...
+        'DisplayName',[metadata.s2name,'=',num2str(metadata.s2(i))],...
+        'color',colors_(i,:)); 
+        hold on;
+    end
+    xlabel(metadata.s1name); ylabel(metadata.dname);title(metadata.title);
+    xlim([metadata.s1(1) metadata.s1(end)]);
+    
+    subplot(122)
+    for i = 1:numel(metadata.s2)
+        errorbar(metadata.s1,omega(:,i),w_err(:,i),'s-',...
+        'LineWidth',2.0,...
+        'DisplayName',[metadata.s2name,'=',num2str(metadata.s2(i))],...
+        'color',colors_(i,:)); 
+        hold on;
+    end
+    xlabel(metadata.s1name); ylabel('$\omega R/c_s$');title(metadata.title);
+    xlim([metadata.s1(1) metadata.s1(end)]);
+    
+    colormap(colors_);
+    clb = colorbar;
+    clim([1 numel(metadata.s2)+1]);
+    clb.Ticks=linspace(metadata.s2(1),metadata.s2(end),numel(metadata.s2));
+    clb.Ticks    =1.5:numel(metadata.s2)+1.5;
+    clb.TickLabels=metadata.s2;
+    clb.Label.String = metadata.s2name;
+    clb.Label.Interpreter = 'latex';
+    clb.Label.FontSize= 18;
+end
+    % clear metadata tosave
+end
diff --git a/wk/lin_scan_script_P_ky.m b/wk/lin_scan_script_P_ky.m
new file mode 100644
index 00000000..b0780091
--- /dev/null
+++ b/wk/lin_scan_script_P_ky.m
@@ -0,0 +1,201 @@
+%% 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
+wkdir = pwd;
+gyacomodir = wkdir(1:end-2);
+mpirun     = 'mpirun';
+% mpirun     = '/opt/homebrew/bin/mpirun'; % for macos
+addpath(genpath([gyacomodir,'matlab']))         % Add matlab folder
+addpath(genpath([gyacomodir,'matlab/plot']))    % Add plot folder
+addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute folder
+addpath(genpath([gyacomodir,'matlab/load']))    % Add load folder
+addpath(genpath([gyacomodir,'wk/parameters']))  % Add parameters folder
+
+%% Setup run or load an executable
+RUN     = 1; % To run or just to load
+RERUN   = 0; % rerun if the  data does not exist
+default_plots_options
+% EXECNAME = 'gyacomo23_sp'; % single precision
+EXECNAME = 'gyacomo23_dp'; % double precision
+
+%% Setup parameters
+% run lin_DTT_AB_rho85
+% run lin_DTT_AB_rho98
+% run lin_JET_rho97
+% run lin_Entropy
+% run lin_ITG
+% run lin_RHT
+% rho  = 0.95; TRIANG = 'PT'; READPROF = 1; 
+% prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/'];
+% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson/',TRIANG,'/'];
+% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/'];
+% run lin_DIIID_data
+run lin_DIIID_LM_rho95
+
+%% Change parameters
+% NU   = 1;
+% TAU  = 1;
+NY   = 2;
+DELTA =-0.2; TRIANG = 'NT';
+S_DELTA = DELTA/2;
+% EXBRATE = 0;
+% S_DELTA = min(2.0,S_DELTA);
+% SIGMA_E  = 0.023;
+% NEXC = 0;
+LX   = 120;
+%% Scan parameters
+SIMID = [SIMID,TRIANG,'_scan'];
+P_a   = [2 4 6 8 16]; J_a = [1 2 3 4 8];
+% P_a   = 2;
+% ky_a  = [0.01 0.02 0.05 0.1  0.2  0.5  1.0  2.0  5.0  10.0];
+ky_a  = [0.05 linspace(0.1,1.1,16)]; ky_a = ky_a(1:end-2);
+% ky_a  = 4.0;
+% dt_a  = logspace(-2,-3,numel(ky_a));
+DT    = 0.0005;
+CO    = 'DG';
+DTSAVE3D = 0.002; TMAX = 40;
+%% Scan loop
+% arrays for the result
+g_ky = zeros(numel(ky_a),numel(P_a));
+g_std= g_ky*0;
+w_ky = g_ky*0;
+w_std= g_ky*0;
+j = 1;
+for PMAX = P_a
+    JMAX = J_a(j);
+    i = 1;
+    for ky = ky_a
+        LY   = 2*pi/ky;
+        DTSAVE0D =  1.0;
+        DTSAVE3D =  0.5;
+        %% 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))
+            MVIN =['cd ',LOCALDIR,';'];
+            % RUNG  =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
+            RUNG  =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
+            % RUNG  =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;'];
+            % RUNG  =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
+            % RUNG = ['./../../../bin/gyacomo23_sp 0;'];
+            MVOUT=['cd ',wkdir,';'];
+            system([MVIN,RUNG,MVOUT]);
+        end
+        data_    = compile_results_low_mem(data_,LOCALDIR,00,00);
+        [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+        if numel(data_.Ts0D)>10
+            % 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');
+            options.NORMALIZED = 0; 
+            options.TIME   = data_.Ts3D;
+             % Time window to measure the growth of kx/ky modes
+            options.KY_TW  = [0.7 1.0]*data_.Ts3D(end);
+            options.KX_TW  = [0.7 1.0]*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.GOK2   = 0; % plot gamma/k^2
+            options.fftz.flag = 0; % Set fftz.flag option to 0
+            options.FIELD = 'phi';
+            options.SHOWFIG = 0;
+            [fig, wkykx, ekykx] = mode_growth_meter(data_,options);
+            % [wkykx,ekykx] = compute_growth_rates(data_.PHI(:,:,:,it1:it2),data_.Ts3D(it1:it2));
+            g_ky (i,j) = real(wkykx(2,1));
+            g_std(i,j) = real(ekykx(2,1));
+            w_ky (i,j) = imag(wkykx(2,1));
+            w_std(i,j) = imag(ekykx(2,1));
+            [gmax, ikmax] = max(g_ky(i,j));
+
+            msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
+        end
+        i = i + 1;
+    end
+    j = j + 1;
+end
+
+%% take max growth rate among z coordinate
+y_ = g_ky + 1i*w_ky; 
+e_ = g_std+ 1i*w_std;
+
+%% Save scan results (gamma)
+if(numel(ky_a)>1 || numel(P_a)>1)
+    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,...
+                '_kN_',num2str(K_Ni),...
+                '_',CONAME,'_',num2str(NU),'_be_',num2str(BETA),...,
+                '_d_',num2str(DELTA),'.mat'];
+    metadata.name   = filename;
+    metadata.kymin  = ky;
+    metadata.title  = ['$\nu_{',CONAME,'}=$',num2str(NU),'$\kappa_T=$',num2str(K_Ti),', $\kappa_N=$',num2str(K_Ni)];
+    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]);
+    % plot
+if 1
+    gamma = real(metadata.data); g_err = real(metadata.err);
+    omega = imag(metadata.data); w_err = imag(metadata.err);
+    gamma = gamma.*(gamma>0.025);
+    figure
+    colors_ = jet(numel(metadata.s2));
+    subplot(121)
+    for i = 1:numel(metadata.s2)
+        errorbar(metadata.s1,gamma(:,i),0*g_err(:,i),'s-',...
+        'LineWidth',2.0,...
+        'DisplayName',[metadata.s2name,'=',num2str(metadata.s2(i))],...
+        'color',colors_(i,:)); 
+        hold on;
+    end
+    xlabel(metadata.s1name); ylabel(metadata.dname);title(metadata.title);
+    xlim([metadata.s1(1) metadata.s1(end)]);
+    
+    subplot(122)
+    for i = 1:numel(metadata.s2)
+        errorbar(metadata.s1,omega(:,i),w_err(:,i),'s-',...
+        'LineWidth',2.0,...
+        'DisplayName',[metadata.s2name,'=',num2str(metadata.s2(i))],...
+        'color',colors_(i,:)); 
+        hold on;
+    end
+    xlabel(metadata.s1name); ylabel('$\omega R/c_s$');title(metadata.title);
+    xlim([metadata.s1(1) metadata.s1(end)]);
+    
+    colormap(colors_);
+    clb = colorbar;
+    clim([1 numel(metadata.s2)+1]);
+    clb.Ticks=linspace(metadata.s2(1),metadata.s2(end),numel(metadata.s2));
+    clb.Ticks    =1.5:numel(metadata.s2)+1.5;
+    clb.TickLabels=metadata.s2;
+    clb.Label.String = metadata.s2name;
+    clb.Label.Interpreter = 'latex';
+    clb.Label.FontSize= 18;
+end
+end
diff --git a/wk/load_metadata_scan.m b/wk/load_metadata_scan.m
index 89db4acb..823f506c 100644
--- a/wk/load_metadata_scan.m
+++ b/wk/load_metadata_scan.m
@@ -11,7 +11,12 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
 % datafname = 'lin_DIIID_Oak_Nelson_high_density_NT_scan/6x32_ky_0.01_10_P_2_8_kN_1.0883_LDGK_0.0080915_be_0.0015991.mat';
 % rho = 0.95
 % datagname = 'lin_DIIID_Oak_Nelson_high_density_PT_scan/6x32_ky_0.01_10_P_2_8_kN_0.62888_LDGK_0.0046858_be_0.0048708.mat';
-datafname = 'lin_DIIID_Oak_Nelson_high_density_PT_scan/6x32_ky_0.01_10_P_2_4_kN_0.62888_DGGK_0.0046858_be_0.0048708.mat';
+% datafname = 'lin_DIIID_Oak_Nelson_high_density_PT_scan/6x32_ky_0.01_10_P_2_4_kN_0.62888_DGGK_0.0046858_be_0.0048708.mat';
+% datafname = 'lin_DIIID_LM_rho95_scan/6x32_ky_0.025_1.1_P_2_4_kN_1.7_DGGK_0.02_be_0.000759_d_0.0.mat';
+% datafname = 'lin_DIIID_LM_rho95_scan/6x32_ky_0.025_1.1_P_2_6_kN_1.7_DGGK_0.02_be_0.000759_d_0.mat';
+% datafname = 'lin_DIIID_LM_rho95_scan/6x32_ky_0.1_1.1_P_2_16_kN_1.7_DGGK_0.02_be_0.000759_d_0.mat';
+% datafname = 'lin_DIIID_LM_rho95_scan/6x32_ky_0.05_0.96667_P_2_16_kN_1.7_DGGK_0.02_be_0.000759_d_0.mat';
+datafname = 'lin_DIIID_LM_rho95PT_scan/6x32_ky_0.05_0.96667_P_2_16_kN_1.7_DGGK_0.02_be_0.000759_d_0.2.mat';
 %% Chose if we filter gamma>0.05
 FILTERGAMMA = 1;
 
@@ -49,14 +54,14 @@ colors_ = jet(numel(d.s2));
 subplot(121)
 for i = 1:numel(d.s2)
     % plot(d.s1,gamma(:,i),'s-',...
-    % plot(d.s1(gamma(:,i)>0),gamma((gamma(:,i)>0),i),'s-',...
-    %     'LineWidth',2.0,...
-    %     'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
-    %     'color',colors_(i,:)); 
-    errorbar(d.s1,gamma(:,i),g_err(:,i),'s-',...
-    'LineWidth',2.0,...
-    'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
-    'color',colors_(i,:)); 
+    plot(d.s1(gamma(:,i)>0),gamma((gamma(:,i)>0),i),'s-',...
+        'LineWidth',2.0,...
+        'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
+        'color',colors_(i,:)); 
+    % errorbar(d.s1,gamma(:,i),g_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);
@@ -64,15 +69,14 @@ xlim([d.s1(1) d.s1(end)]);
 
 subplot(122)
 for i = 1:numel(d.s2)
-    % plot(d.s1,gamma(:,i),'s-',...
-    % plot(d.s1(gamma(:,i)>0),gamma((gamma(:,i)>0),i),'s-',...
-    %     'LineWidth',2.0,...
-    %     'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
-    %     'color',colors_(i,:)); 
-    errorbar(d.s1,omega(:,i),w_err(:,i),'s-',...
-    'LineWidth',2.0,...
-    'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
-    'color',colors_(i,:)); 
+    plot(d.s1,omega(:,i),'s-',...
+        'LineWidth',2.0,...
+        'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
+        'color',colors_(i,:)); 
+    % errorbar(d.s1,omega(:,i),w_err(:,i),'s-',...
+    % 'LineWidth',2.0,...
+    % 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
+    % 'color',colors_(i,:)); 
     hold on;
 end
 xlabel(d.s1name); ylabel('$\omega R/c_s$');title(d.title);
diff --git a/wk/parameters/lin_AUG.m b/wk/parameters/lin_AUG.m
index bfa0eddb..c3714ec1 100644
--- a/wk/parameters/lin_AUG.m
+++ b/wk/parameters/lin_AUG.m
@@ -91,7 +91,10 @@ 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
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 ADIAB_I = 0;          % adiabatic ion model
\ No newline at end of file
diff --git a/wk/parameters/lin_DIIID_AUSTIN.m b/wk/parameters/lin_DIIID_AUSTIN.m
index 1ba1c97d..d5bc43ca 100644
--- a/wk/parameters/lin_DIIID_AUSTIN.m
+++ b/wk/parameters/lin_DIIID_AUSTIN.m
@@ -123,7 +123,10 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 1.0e-5; % Initial noise amplitude
 BCKGD0  = 0.0e-5;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 ADIAB_I = 0;          % adiabatic ion model
\ No newline at end of file
diff --git a/wk/parameters/lin_DIIID_LM_rho90.m b/wk/parameters/lin_DIIID_LM_rho90.m
index c699e5e9..85c52caf 100644
--- a/wk/parameters/lin_DIIID_LM_rho90.m
+++ b/wk/parameters/lin_DIIID_LM_rho90.m
@@ -93,7 +93,10 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 1.0e-5; % Initial noise amplitude
 BCKGD0  = 0.0e-5;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 ADIAB_I = 0;          % adiabatic ion model
\ No newline at end of file
diff --git a/wk/parameters/lin_DIIID_LM_rho95.m b/wk/parameters/lin_DIIID_LM_rho95.m
index 3dfccf17..ba750005 100644
--- a/wk/parameters/lin_DIIID_LM_rho95.m
+++ b/wk/parameters/lin_DIIID_LM_rho95.m
@@ -91,8 +91,11 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 1.0e-5; % Initial noise amplitude
 BCKGD0  = 0.0e-5;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 ADIAB_I = 0;          % adiabatic ion model
 EXBRATE = 0;
\ No newline at end of file
diff --git a/wk/parameters/lin_DIIID_LM_rho95_HEL.m b/wk/parameters/lin_DIIID_LM_rho95_HEL.m
index cbd94090..0a6515a6 100644
--- a/wk/parameters/lin_DIIID_LM_rho95_HEL.m
+++ b/wk/parameters/lin_DIIID_LM_rho95_HEL.m
@@ -93,8 +93,11 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 1.0e-5; % Initial noise amplitude
 BCKGD0  = 0.0e-5;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 ADIAB_I = 0;          % adiabatic ion model
 EXBRATE = 0;
\ No newline at end of file
diff --git a/wk/parameters/lin_DIIID_data.m b/wk/parameters/lin_DIIID_data.m
index 043de849..50ac9c72 100644
--- a/wk/parameters/lin_DIIID_data.m
+++ b/wk/parameters/lin_DIIID_data.m
@@ -128,7 +128,10 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 1.0e-5; % Initial noise amplitude
 BCKGD0  = 0.0e-5;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 ADIAB_I = 0;          % adiabatic ion model
\ No newline at end of file
diff --git a/wk/parameters/lin_DTT_HM_rho85.m b/wk/parameters/lin_DTT_HM_rho85.m
index dfdf7fa1..8fb86b3a 100644
--- a/wk/parameters/lin_DTT_HM_rho85.m
+++ b/wk/parameters/lin_DTT_HM_rho85.m
@@ -102,7 +102,10 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 1.0e-5; % Initial noise amplitude
 BCKGD0  = 0.0e-5;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 ADIAB_I = 0;          % adiabatic ion model
\ No newline at end of file
diff --git a/wk/parameters/lin_DTT_HM_rho98.m b/wk/parameters/lin_DTT_HM_rho98.m
index 12c49a63..b0e204c0 100644
--- a/wk/parameters/lin_DTT_HM_rho98.m
+++ b/wk/parameters/lin_DTT_HM_rho98.m
@@ -103,7 +103,10 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 1.0e-5; % Initial noise amplitude
 BCKGD0  = 0.0e-5;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 ADIAB_I = 0;          % adiabatic ion model
\ No newline at end of file
diff --git a/wk/parameters/lin_ETG_adiab_i.m b/wk/parameters/lin_ETG_adiab_i.m
index 219bcce7..ba65082d 100644
--- a/wk/parameters/lin_ETG_adiab_i.m
+++ b/wk/parameters/lin_ETG_adiab_i.m
@@ -93,6 +93,9 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 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
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
\ No newline at end of file
diff --git a/wk/parameters/lin_Entropy.m b/wk/parameters/lin_Entropy.m
index 31c57f71..ef242a93 100644
--- a/wk/parameters/lin_Entropy.m
+++ b/wk/parameters/lin_Entropy.m
@@ -93,6 +93,9 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 1.0e-4; % Initial noise amplitude
 BCKGD0  = 0.0e-8;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
diff --git a/wk/parameters/lin_GASTD.m b/wk/parameters/lin_GASTD.m
index 14d13da7..0fcb56a9 100644
--- a/wk/parameters/lin_GASTD.m
+++ b/wk/parameters/lin_GASTD.m
@@ -93,7 +93,10 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 1.0e-5; % Initial noise amplitude
 BCKGD0  = 0.0e-5;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 ADIAB_I = 0;          % adiabatic ion model
\ No newline at end of file
diff --git a/wk/parameters/lin_ITG.m b/wk/parameters/lin_ITG.m
index d1c06b66..e93816e2 100644
--- a/wk/parameters/lin_ITG.m
+++ b/wk/parameters/lin_ITG.m
@@ -87,8 +87,11 @@ 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
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 S_KAPPA = 0.0;
 S_DELTA = 0.0;
diff --git a/wk/parameters/lin_Ivanov.m b/wk/parameters/lin_Ivanov.m
index a578f0bc..ef267364 100644
--- a/wk/parameters/lin_Ivanov.m
+++ b/wk/parameters/lin_Ivanov.m
@@ -93,8 +93,11 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 1.0e-5; % Initial noise amplitude
 BCKGD0  = 0.0;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1.0; % Cutoff for collision operator
 S_KAPPA = 0.0;
 S_DELTA = 0.0;
diff --git a/wk/parameters/lin_JET_rho97.m b/wk/parameters/lin_JET_rho97.m
index 72ae55cd..9fbd5884 100644
--- a/wk/parameters/lin_JET_rho97.m
+++ b/wk/parameters/lin_JET_rho97.m
@@ -127,7 +127,10 @@ 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
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 ADIAB_I = 0;          % adiabatic ion model
 ADIAB_E = (NA==1);          % adiabatic electron model
\ No newline at end of file
diff --git a/wk/parameters/lin_KBM.m b/wk/parameters/lin_KBM.m
index fc2000de..b848aa41 100644
--- a/wk/parameters/lin_KBM.m
+++ b/wk/parameters/lin_KBM.m
@@ -88,6 +88,9 @@ 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
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
diff --git a/wk/parameters/lin_RHT.m b/wk/parameters/lin_RHT.m
index 0f39498d..a7e2e9a3 100644
--- a/wk/parameters/lin_RHT.m
+++ b/wk/parameters/lin_RHT.m
@@ -89,8 +89,11 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 0.0e-5; % Initial noise amplitude
 BCKGD0  = 1.0;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 S_KAPPA = 0.0;
 S_DELTA = 0.0;
diff --git a/wk/parameters/lin_STEP_EC_HD_psi49.m b/wk/parameters/lin_STEP_EC_HD_psi49.m
index d1fa57ef..8ca4cb51 100644
--- a/wk/parameters/lin_STEP_EC_HD_psi49.m
+++ b/wk/parameters/lin_STEP_EC_HD_psi49.m
@@ -110,7 +110,10 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 1.0e-5; % Initial noise amplitude
 BCKGD0  = 0.0e-5;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 ADIAB_I = 0;          % adiabatic ion model
\ No newline at end of file
diff --git a/wk/parameters/lin_STEP_EC_HD_psi71.m b/wk/parameters/lin_STEP_EC_HD_psi71.m
index 9881468e..41c60816 100644
--- a/wk/parameters/lin_STEP_EC_HD_psi71.m
+++ b/wk/parameters/lin_STEP_EC_HD_psi71.m
@@ -110,7 +110,10 @@ MU_J    = 0.0;    % Hyperdiffusivity coefficient for Laguerre
 LAMBDAD = 0.0;    % Lambda Debye
 NOISE0  = 1.0e-5; % Initial noise amplitude
 BCKGD0  = 0.0e-5;    % Initial background
-k_gB   = 1.0;     % Magnetic gradient strength
-k_cB   = 1.0;     % Magnetic curvature strength
+K_gB   = 1.0;     % Magnetic gradient tuner  (1 is the real value)
+K_cB   = 1.0;     % Magnetic curvature tuner (1 is the real value)
+K_mB   = 1.0;     % mirror force tuner       (1 is the real value)
+K_tB   = 1.0;     % trapping term tuner      (1 is the real value)
+K_ldB  = 1.0;     % Landau damping tuner     (1 is the real value)
 COLL_KCUT = 1; % Cutoff for collision operator
 ADIAB_I = 0;          % adiabatic ion model
\ No newline at end of file
-- 
GitLab