diff --git a/fluidsim/util/util.py b/fluidsim/util/util.py
index 7160f8322242561b247f1268b4c7b5713975e436_Zmx1aWRzaW0vdXRpbC91dGlsLnB5..559f4f76d333535d23881481bd6280f9d6231225_Zmx1aWRzaW0vdXRpbC91dGlsLnB5 100644
--- a/fluidsim/util/util.py
+++ b/fluidsim/util/util.py
@@ -292,6 +292,9 @@
     if merge_missing_params:
         merge_params(params, solver.Simul.create_default_params())
 
+    if hasattr(params, "ONLY_COARSE_OPER") and params.ONLY_COARSE_OPER:
+        params.ONLY_COARSE_OPER = False
+
     params.path_run = path_dir
     params.NEW_DIR_RESULTS = False
     if modif_save_params:
diff --git a/scripts/bugs/init_from_file_not_running.py b/scripts/bugs/init_from_file_not_running.py
deleted file mode 100644
index 7160f8322242561b247f1268b4c7b5713975e436_c2NyaXB0cy9idWdzL2luaXRfZnJvbV9maWxlX25vdF9ydW5uaW5nLnB5..0000000000000000000000000000000000000000
--- a/scripts/bugs/init_from_file_not_running.py
+++ /dev/null
@@ -1,72 +0,0 @@
-"""
-init_from_file_not_running.py
-=============================
-
-# Environment
---------------
-Python 3.6.6 (Conda)
-
-# Explanation
--------------
-It does not compute when params.init_fields.type = "from_file" for
-both sequential and MPI.
-
-Bug in my terminal and cluster.
-
-Both fluidfft and fluidsim are the last versions (tests OK).
-
-# Traceback
-------------
-No traceback.
-
-# To run the bug
-----------------
-sequential
-python init_from_file_not_running.py
-
-parallel (4 proc.)
-mpirun -np 4 init_from_file_not_running.py
-
-"""
-from math import pi
-
-import os
-from glob import glob
-from fluidsim.solvers.ns2d.strat.solver import Simul
-
-# Create parameters
-params = Simul.create_default_params()
-
-# Operator parameters
-params.oper.nx = nx = 960
-params.oper.ny = nx // 4
-params.oper.Lx = Lx = 2 * pi
-params.oper.Ly = Lx / 4
-params.oper.NO_SHEAR_MODES = True
-params.oper.coef_dealiasing = 0.6666
-params.oper.type_fft = None
-
-# Without forcing
-params.forcing.enable = False
-
-# Parameters time stepping
-params.time_stepping.USE_CFL = True
-params.time_stepping.t_end = 2.
-params.time_stepping.cfl_coef_group = None
-
-# Parameters initialization
-path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim"
-paths_sims = sorted(glob(os.path.join(path_root, f"sim{nx}_no_shear_modes",
-                                      "NS2D*")))
-path_file = glob(paths_sims[0] + "/state_phys*")[-10]
-params.init_fields.type = "from_file"
-params.init_fields.from_file.path = path_file
-
-# Parameters output
-params.output.HAS_TO_SAVE = False
-params.output.periods_print.print_stdout = 1e-2
-
-# Launch...
-sim = Simul(params)
-sim.time_stepping.start()
-
diff --git a/scripts/ns2d.strat/compute_length_scales.py b/scripts/ns2d.strat/compute_length_scales.py
new file mode 100644
index 0000000000000000000000000000000000000000..559f4f76d333535d23881481bd6280f9d6231225_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfbGVuZ3RoX3NjYWxlcy5weQ==
--- /dev/null
+++ b/scripts/ns2d.strat/compute_length_scales.py
@@ -0,0 +1,84 @@
+"""
+compute_length_scales.py
+==========================
+
+Computes vertical length from Appendix B. Brethouwer 2007.
+"""
+
+import h5py
+import numpy as np
+
+from fluidsim import load_params_simul
+
+def compute_length_scales(path_simulation, tmin=None):
+    """
+    Compute length scales.
+    """
+    # Print out
+    res_out = float(path_simulation.split("NS2D.strat_")[1].split("x")[0])
+    gamma_str = path_simulation.split("_gamma")[1].split("_")[0]
+    if gamma_str.startswith("0"):
+        gamma_out = float(gamma_str[0] + "." + gamma_str[1])
+    else:
+        gamma_out = float(gamma_str)
+
+    print("Compute dissipation nx = {} and gamma {}..".format(res_out, gamma_out))
+
+    # Load object parameters
+    params = load_params_simul(path_simulation)
+
+    # Load data energy spectra
+    with h5py.File(path_simulation + "/spectra1D.h5", "r") as f:
+        times = f["times"].value
+        kx = f["kxE"].value
+        ky = f["kyE"].value
+        spectrum1Dkx_EK_ux = f["spectrum1Dkx_EK_ux"].value
+        spectrum1Dky_EK_uy = f["spectrum1Dky_EK_uy"].value
+
+    # Temporal average spectra
+    dt = np.median(np.diff(times))
+    if not tmin:
+        nb_files = 100
+        tmin = np.max(times) - (nb_files * dt)
+    itmin = np.argmin(abs(times - tmin))
+
+    # Compute time average
+    spectrum1Dkx_EK_ux = np.mean(spectrum1Dkx_EK_ux[itmin:, :], axis=0)
+    spectrum1Dky_EK_uy = np.mean(spectrum1Dky_EK_uy[itmin:, :], axis=0)
+
+    # Compute spectrum dealiased
+    kxmax_dealiasing = params.oper.coef_dealiasing * np.max(abs(kx))
+    kymax_dealiasing = params.oper.coef_dealiasing * np.max(abs(ky))
+
+    ikx_dealiasing = np.argmin(abs(kx - kxmax_dealiasing))
+    iky_dealiasing = np.argmin(abs(ky - kymax_dealiasing))
+
+    kx = kx[0:ikx_dealiasing]
+    ky = ky[0:iky_dealiasing]
+    spectrum1Dkx_EK_ux = spectrum1Dkx_EK_ux[:ikx_dealiasing]
+    spectrum1Dky_EK_uy = spectrum1Dky_EK_uy[:iky_dealiasing]
+
+    # Compute delta_k
+    delta_kx = np.median(np.diff(abs(kx)))
+    delta_ky = np.median(np.diff(abs(ky)))
+
+    # Compute vertical length scale
+    lx = (np.sum(spectrum1Dkx_EK_ux * delta_kx) /
+          np.sum(kx * spectrum1Dkx_EK_ux * delta_kx))
+
+    lz = (np.sum(spectrum1Dky_EK_uy * delta_ky) /
+          np.sum(ky * spectrum1Dky_EK_uy * delta_ky))
+
+    return lx, lz
+
+
+
+if __name__ == "__main__":
+    # path_simulation = "/fsnet/project/meige/2015/15DELDUCA/DataSim/sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma02_2018-08-14_09-59-55"
+
+    path_simulation = "/fsnet/project/meige/2015/15DELDUCA/DataSim/sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma1_2018-08-14_10-01-22"
+
+    lx, lz = compute_length_scales(path_simulation)
+
+    print("lx = ", lx)
+    print("lz = ", lz)
diff --git a/scripts/ns2d.strat/compute_ratio_dissipation.py b/scripts/ns2d.strat/compute_ratio_dissipation.py
index 7160f8322242561b247f1268b4c7b5713975e436_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfcmF0aW9fZGlzc2lwYXRpb24ucHk=..559f4f76d333535d23881481bd6280f9d6231225_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfcmF0aW9fZGlzc2lwYXRpb24ucHk= 100644
--- a/scripts/ns2d.strat/compute_ratio_dissipation.py
+++ b/scripts/ns2d.strat/compute_ratio_dissipation.py
@@ -54,4 +54,17 @@
             params.forcing.nkmax_forcing * max(delta_kx, delta_kz))
     ik_fx = np.argmin(abs(kx - k_fx ))
 
-    return D_kx[ik_fx] / D_kx[-1]
+    ratio_dissipation_old = D_kx[ik_fx] / D_kx[-1]
+
+    # Compute wave-number for 1/2 dissipation
+    ik_half = np.argmin(abs(D_kx - (D_kx[-1] / 2)))
+
+    return D_kx[ik_half] / D_kx[ik_fx]
+
+if __name__ == "__main__":
+    path_simulation = "/fsnet/project/meige/2015/15DELDUCA/DataSim/" + \
+                      "sim960_no_shear_modes/" + \
+                      "NS2D.strat_960x240_S2pix1.571_F07_gamma02_2018-08-10_14-45-06"
+
+    ratio = compute_ratio_dissipation(path_simulation, tmin=None)
+    print(f"Ratio between ik_half / ik_fx = {ratio}")
diff --git a/scripts/ns2d.strat/compute_vertical_length.py b/scripts/ns2d.strat/compute_vertical_length.py
new file mode 100644
index 0000000000000000000000000000000000000000..559f4f76d333535d23881481bd6280f9d6231225_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfdmVydGljYWxfbGVuZ3RoLnB5
--- /dev/null
+++ b/scripts/ns2d.strat/compute_vertical_length.py
@@ -0,0 +1,85 @@
+"""
+compute_length_scales.py
+==========================
+
+Computes vertical length from Appendix B. Brethouwer 2007.
+"""
+
+import h5py
+import numpy as np
+
+from fluidsim import load_params_simul
+
+def compute_length_scales(path_simulation, tmin=None):
+    """
+    Compute length scales.
+    """
+    # Print out
+    res_out = float(path_simulation.split("NS2D.strat_")[1].split("x")[0])
+    gamma_str = path_simulation.split("_gamma")[1].split("_")[0]
+    if gamma_str.startswith("0"):
+        gamma_out = float(gamma_str[0] + "." + gamma_str[1])
+    else:
+        gamma_out = float(gamma_str)
+
+    print("Compute dissipation nx = {} and gamma {}..".format(res_out, gamma_out))
+
+    # Load object parameters
+    params = load_params_simul(path_simulation)
+
+    # Load data energy spectra
+    with h5py.File(path_simulation + "/spectra1D.h5", "r") as f:
+        times = f["times"].value
+        kx = f["kxE"].value
+        ky = f["kyE"].value
+        spectrum1Dkx_EK_ux = f["spectrum1Dkx_EK_ux"].value
+        spectrum1Dky_EK_uy = f["spectrum1Dky_EK_uy"].value
+
+    # Temporal average spectra
+    dt = np.median(np.diff(times))
+    if not tmin:
+        nb_files = 100
+        tmin = np.max(times) - (nb_files * dt)
+    itmin = np.argmin(abs(times - tmin))
+
+    # Compute time average
+    spectrum1Dkx_EK_ux = np.mean(spectrum1Dkx_EK_ux[itmin:, :], axis=0)
+    spectrum1Dky_EK_uy = np.mean(spectrum1Dky_EK_uy[itmin:, :], axis=0)
+
+    # Compute spectrum dealiased
+    kxmax_dealiasing = params.oper.coef_dealiasing * np.max(abs(kx))
+    kymax_dealiasing = params.oper.coef_dealiasing * np.max(abs(ky))
+
+    ikx_dealiasing = np.argmin(abs(kx - kxmax_dealiasing))
+    iky_dealiasing = np.argmin(abs(ky - kymax_dealiasing))
+
+    kx = kx[0:ikx_dealiasing]
+    ky = ky[0:iky_dealiasing]
+    spectrum1Dkx_EK_ux = spectrum1Dkx_EK_ux[:ikx_dealiasing]
+    spectrum1Dky_EK_uy = spectrum1Dky_EK_uy[:iky_dealiasing]
+
+    # Compute delta_k
+    delta_kx = np.median(np.diff(abs(kx)))
+    delta_ky = np.median(np.diff(abs(ky)))
+
+    # Compute vertical length scale
+    lx = (np.sum(spectrum1Dkx_EK_ux * delta_kx) /
+          np.sum(kx * spectrum1Dkx_EK_ux * delta_kx))
+
+    lz = (np.sum(spectrum1Dky_EK_uy * delta_ky) /
+          np.sum(ky * spectrum1Dky_EK_uy * delta_ky))
+
+    return lx, lz
+
+
+
+if __name__ == "__main__":
+    # path_simulation = "/fsnet/project/meige/2015/15DELDUCA/DataSim/sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma02_2018-08-14_09-59-55"
+
+    path_simulation = "/fsnet/project/meige/2015/15DELDUCA/DataSim/sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma1_2018-08-14_10-01-22"
+
+    print(path_simulation)
+    lx, lz = compute_length_scales(path_simulation)
+
+    print("lx = ", lx)
+    print("lz = ", lz)
diff --git a/scripts/ns2d.strat/make_table_parameters.py b/scripts/ns2d.strat/make_table_parameters.py
index 7160f8322242561b247f1268b4c7b5713975e436_c2NyaXB0cy9uczJkLnN0cmF0L21ha2VfdGFibGVfcGFyYW1ldGVycy5weQ==..559f4f76d333535d23881481bd6280f9d6231225_c2NyaXB0cy9uczJkLnN0cmF0L21ha2VfdGFibGVfcGFyYW1ldGVycy5weQ== 100644
--- a/scripts/ns2d.strat/make_table_parameters.py
+++ b/scripts/ns2d.strat/make_table_parameters.py
@@ -13,8 +13,8 @@
 from fluidsim import load_sim_for_plot
 
 # Argparse arguments
-nx = 1920
-MAKE_TABLE = False
+nx = 3840
+MAKE_TABLE = True
 
 # Parameters script
 n_files_tmean = 100
@@ -24,6 +24,8 @@
 
 if nx == 1920:
     directory = "sim1920_no_shear_modes"
+elif nx == 3840:
+    directory = "sim3840_modif_res_no_shear_modes"
 
 path_simulations = sorted(glob(os.path.join(path_root, directory, "NS2D*")))
 
@@ -62,8 +64,8 @@
 
     with h5py.File(path_spectra, "r") as f:
         times_spectra = f["times"].value
-        kx = f["kyE"].value
+        kx = f["kxE"].value
         spectrum1Dkx_EK_ux = f["spectrum1Dkx_EK_ux"].value
 
     spectrum1Dkx_EK_ux = np.mean(spectrum1Dkx_EK_ux[-100:, :], axis=0)
     ## Remove modes with dealiasing
@@ -66,8 +68,9 @@
         spectrum1Dkx_EK_ux = f["spectrum1Dkx_EK_ux"].value
 
     spectrum1Dkx_EK_ux = np.mean(spectrum1Dkx_EK_ux[-100:, :], axis=0)
     ## Remove modes with dealiasing
-    ikxmax = np.argmin(abs(kx - sim.oper.kxmax_dealiasing))
+    kxmax_dealiasing = sim.params.oper.coef_dealiasing * np.max(abs(kx))
+    ikxmax = np.argmin(abs(kx - kxmax_dealiasing))
     kx = kx[:ikxmax]
     spectrum1Dkx_EK_ux = spectrum1Dkx_EK_ux[:ikxmax]
     delta_kx = sim.oper.deltakx
diff --git a/scripts/ns2d.strat/plot_length_scales.py b/scripts/ns2d.strat/plot_length_scales.py
new file mode 100644
index 0000000000000000000000000000000000000000..559f4f76d333535d23881481bd6280f9d6231225_c2NyaXB0cy9uczJkLnN0cmF0L3Bsb3RfbGVuZ3RoX3NjYWxlcy5weQ==
--- /dev/null
+++ b/scripts/ns2d.strat/plot_length_scales.py
@@ -0,0 +1,151 @@
+"""
+plot_length_scales.py
+=======================
+
+"""
+import os
+import numpy as np
+import h5py
+import matplotlib.pyplot as plt
+
+from glob import glob
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+from compute_anisotropy import compute_anisotropy
+from compute_ratio_dissipation import compute_ratio_dissipation
+from compute_reynolds_froude import compute_buoyancy_reynolds
+from compute_length_scales import compute_length_scales
+from fluiddyn.output.rcparams import set_rcparams
+from fluidsim import load_params_simul
+def _get_resolution_from_dir(path_simulation):
+    return path_simulation.split("NS2D.strat_")[1].split("x")[0]
+def _get_gamma_str_from_path(path_simulation):
+    return path_simulation.split("_gamma")[1].split("_")[0]
+
+SAVE = False
+
+# Create path simulations
+path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim"
+directories = ["sim960_no_shear_modes",
+               "sim960_no_shear_modes_transitory",
+               "sim1920_no_shear_modes",
+               "sim1920_modif_res_no_shear_modes",
+               "sim3840_modif_res_no_shear_modes"]
+
+# directories = ["sim960_no_shear_modes",
+#                "sim1920_no_shear_modes",
+#                "sim1920_modif_res_no_shear_modes"]
+
+# directories = ["sim960_no_shear_modes"]
+paths_simulations = []
+for directory in directories:
+    paths_simulations += sorted(glob(os.path.join(path_root, directory, "NS2D*")))
+
+froudes = []
+reynoldsb = []
+reynolds8 = []
+anisotropies = []
+dissipations = []
+markers = []
+lzs_billant = []
+
+set_rcparams(fontsize=14, for_article=True)
+
+fig, ax = plt.subplots()
+ax.set_ylabel(r"$l_zN/U$", fontsize=18)
+ax.set_xlabel(r"$\mathcal{R}_8$", fontsize=18)
+ax.set_xscale("log")
+ax.set_yscale("log")
+# ax.text(0.6, 2e-12, r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", fontsize=16)
+ax.tick_params(axis="x", labelsize=18)
+ax.tick_params(axis="y", labelsize=18)
+# ax.set_xlim([0.003, 1])
+# ax.set_ylim([1e-9, 10])
+
+
+for path in paths_simulations:
+    gamma_str = _get_gamma_str_from_path(path)
+    if gamma_str.startswith("2"):
+        continue
+    else:
+        F_h, Re_8, R_b = compute_buoyancy_reynolds(path)
+        anisotropy = compute_anisotropy(path)
+        dissipation = compute_ratio_dissipation(path)
+        res = _get_resolution_from_dir(path)
+        lx, lz = compute_length_scales(path)
+
+        # Compute
+        ux_rms = []
+        path_phys_files = glob(path + "/state_phys_t*")
+        for path_file in path_phys_files[-10:]:
+            with h5py.File(path_file, "r") as f:
+                ux = f["state_phys"]["ux"].value
+            ux_rms.append(np.sqrt(np.mean(ux**2)))
+
+        # Load parameters
+        params = load_params_simul(path)
+
+        lzs_billant.append(lz * params.N / np.mean(ux_rms))
+
+        froudes.append(F_h)
+        reynoldsb.append(R_b)
+        reynolds8.append(Re_8)
+        anisotropies.append(np.log10(2 * anisotropy))
+        dissipations.append(np.log10(dissipation))
+
+        if res == "960":
+            markers.append("o")
+        elif res == "1920":
+            markers.append("s")
+        elif res == "3840":
+            markers.append("^")
+
+        print("F_h", F_h)
+        print("Re_8", Re_8)
+        print("R_b", R_b)
+
+        for _f, _r, _a, _d, _m, _l in zip(froudes, reynoldsb, anisotropies, dissipations, markers, lzs_billant):
+            scatter = ax.scatter(_r, _l, s=100 * (1**1), c="k", vmin=0, vmax=1, marker=_m, alpha=0.7)
+
+        # for _f, _r, _a, _d, _m, _l in zip(froudes, reynoldsb, anisotropies, dissipations, markers, lzs_billant):
+        #     scatter = ax.scatter(_r, _l, s=2000 * (_a**1) + 20, c=_d, vmin=0, vmax=1, marker=_m, alpha=0.7)
+
+
+# ax.scatter(0.1, 1e-5, marker="o", s=2000 * (np.log10(2*0.5)**1) + 20, color="red")
+# ax.scatter(0.1, 1e-6, marker="o", s=2000 * (np.log10(2*0.75)**1) + 20, color="red")
+# ax.scatter(0.1, 2e-8, marker="o", s=2000 * (np.log10(2)**1) + 20, color="red")
+
+# ax.text(0.16, 8e-6, r"$anisotropy = 0$", color="red", fontsize=14)
+# ax.text(0.16, 5e-7, r"$anisotropy = 0.5$", color="red", fontsize=14)
+# ax.text(0.16, 2e-8, r"$anisotropy = 1$", color="red", fontsize=14)
+
+# divider = make_axes_locatable(ax)
+# cax = divider.append_axes("right", size="5%", pad=0.05)
+# cax.tick_params(labelsize=14)
+# fig.colorbar(scatter, cax=cax)
+
+# Legend
+import matplotlib.lines as mlines
+import matplotlib.pyplot as plt
+
+blue_star = mlines.Line2D([], [], color='k', marker='o', linestyle='None',
+                          markersize=8, label=r'$n_x = 960$')
+red_square = mlines.Line2D([], [], color='k', marker='s', linestyle='None',
+                           markersize=8, label=r'$n_x = 1920$')
+purple_triangle = mlines.Line2D([], [], color='k', marker='^', linestyle='None',
+                                markersize=8, label=r'$n_x = 3840$')
+
+ax.legend(handles=[blue_star, red_square, purple_triangle],
+          loc="upper center",
+          bbox_to_anchor=(0.5,1.1),
+          borderaxespad=0.,
+          ncol=len(markers),
+          handletextpad=0.1,
+          fontsize=14)
+
+fig.tight_layout(pad=0.4)
+
+if SAVE:
+    path_save = "/home/users/calpelin7m/Phd/docs/EFMC18/figures"
+    fig.savefig(path_save + "/vertical_length_billant.eps", format="eps",
+                bbox_inches="tight")
+plt.show()
diff --git a/scripts/ns2d.strat/plot_reynolds_froude.py b/scripts/ns2d.strat/plot_reynolds_froude.py
index 7160f8322242561b247f1268b4c7b5713975e436_c2NyaXB0cy9uczJkLnN0cmF0L3Bsb3RfcmV5bm9sZHNfZnJvdWRlLnB5..559f4f76d333535d23881481bd6280f9d6231225_c2NyaXB0cy9uczJkLnN0cmF0L3Bsb3RfcmV5bm9sZHNfZnJvdWRlLnB5 100644
--- a/scripts/ns2d.strat/plot_reynolds_froude.py
+++ b/scripts/ns2d.strat/plot_reynolds_froude.py
@@ -11,7 +11,7 @@
 import matplotlib.pyplot as plt
 
 from glob import glob
-
+from mpl_toolkits.axes_grid1 import make_axes_locatable
 from compute_anisotropy import compute_anisotropy
 from compute_ratio_dissipation import compute_ratio_dissipation
 from compute_reynolds_froude import compute_buoyancy_reynolds
@@ -19,7 +19,11 @@
 
 def _get_resolution_from_dir(path_simulation):
     return path_simulation.split("NS2D.strat_")[1].split("x")[0]
+def _get_gamma_str_from_path(path_simulation):
+    return path_simulation.split("_gamma")[1].split("_")[0]
+
+SAVE = False
 
 # Create path simulations
 path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim"
 directories = ["sim960_no_shear_modes",
@@ -22,8 +26,9 @@
 
 # Create path simulations
 path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim"
 directories = ["sim960_no_shear_modes",
+               "sim960_no_shear_modes_transitory",
                "sim1920_no_shear_modes",
                "sim1920_modif_res_no_shear_modes",
                "sim3840_modif_res_no_shear_modes"]
 
@@ -26,9 +31,9 @@
                "sim1920_no_shear_modes",
                "sim1920_modif_res_no_shear_modes",
                "sim3840_modif_res_no_shear_modes"]
 
-# directories = ["sim960_no_shear_modes",
-#                "sim1920_no_shear_modes"]
+# directories = ["sim960_no_shear_modes"]
+
 
 paths_simulations = []
 for directory in directories:
@@ -36,6 +41,7 @@
 
 froudes = []
 reynoldsb = []
+reynolds8 = []
 anisotropies = []
 dissipations = []
 markers = []
@@ -43,7 +49,7 @@
 set_rcparams(fontsize=14, for_article=True)
 
 fig, ax = plt.subplots()
-ax.set_xlabel(r"$F_h$")
-ax.set_ylabel(r"$\mathcal{R}$")
+ax.set_xlabel(r"$F_h$", fontsize=18)
+ax.set_ylabel(r"$\mathcal{R}_8$", fontsize=18)
 ax.set_xscale("log")
 ax.set_yscale("log")
@@ -48,6 +54,10 @@
 ax.set_xscale("log")
 ax.set_yscale("log")
-fig.text(0.8, 4e-7, r"$\frac{D(k_{fx})}{D(k_x)}$", fontsize=16)
+ax.text(0.6, 2e-12, r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", fontsize=16)
+ax.tick_params(axis="x", labelsize=18)
+ax.tick_params(axis="y", labelsize=18)
+ax.set_xlim([0.003, 1])
+ax.set_ylim([1e-9, 1e7])
 
 
 for path in paths_simulations:
@@ -51,13 +61,12 @@
 
 
 for path in paths_simulations:
-    F_h, Re_8, R_b = compute_buoyancy_reynolds(path)
-    anisotropy = compute_anisotropy(path)
-    dissipation = compute_ratio_dissipation(path)
-    res = _get_resolution_from_dir(path)
-
-    froudes.append(F_h)
-    reynoldsb.append(R_b)
-    anisotropies.append(anisotropy)
-    dissipations.append(dissipation)
+    gamma_str = _get_gamma_str_from_path(path)
+    if gamma_str.startswith("2"):
+        continue
+    else:
+        F_h, Re_8, R_b = compute_buoyancy_reynolds(path)
+        anisotropy = compute_anisotropy(path)
+        dissipation = compute_ratio_dissipation(path)
+        res = _get_resolution_from_dir(path)
 
@@ -63,8 +72,7 @@
 
-    if res == "960":
-        markers.append("o")
-    elif res == "1920":
-        markers.append("s")
-    elif res == "3840":
-        markers.append("^")
+        froudes.append(F_h)
+        reynoldsb.append(R_b)
+        reynolds8.append(Re_8)
+        anisotropies.append(np.log10(2 * anisotropy))
+        dissipations.append(np.log10(dissipation))
 
@@ -70,6 +78,16 @@
 
-    print("F_h", F_h)
-    print("Re_8", Re_8)
-    print("R_b", R_b)
+        if res == "960":
+            markers.append("o")
+        elif res == "1920":
+            markers.append("s")
+        elif res == "3840":
+            markers.append("^")
+
+        print("F_h", F_h)
+        print("Re_8", Re_8)
+        print("R_b", R_b)
+
+        # for _f, _r, _a, _d, _m in zip(froudes, reynoldsb, anisotropies, dissipations, markers):
+        #     scatter = ax.scatter(_f, _r, s=2000 * (_a**1) + 20, c=_d, vmin=0, vmax=1, marker=_m, alpha=0.7)
 
 for _f, _r, _a, _d, _m in zip(froudes, reynoldsb, anisotropies, dissipations, markers):
@@ -74,5 +92,4 @@
 
 for _f, _r, _a, _d, _m in zip(froudes, reynoldsb, anisotropies, dissipations, markers):
-    scatter = ax.scatter(_f, _r, s=250 * (_a**2), c=_d, vmin=0, vmax=0.3, marker=_m)
-# plt.show()
+            scatter = ax.scatter(_f, _r, s=2000 * (_a**1) + 20, c=_d, vmin=0, vmax=1, marker=_m, alpha=0.7)
 
@@ -78,11 +95,14 @@
 
-# areas = 250 * np.asarray(anisotropies)**2
-# scatter = ax.scatter(froudes, reynoldsb, s=areas, c=dissipations, alpha=0.7, vmin=0, vmax=0.3)
-ax.scatter(max(froudes), 1e2 * min(reynoldsb), s=250 * np.asarray(0.5)**2, c="red")
-ax.scatter(max(froudes), min(reynoldsb), s=250 * np.asarray(1.0)**2, c="red")
-ax.text(0.12, 90 * min(reynoldsb), "anisotropy=0", fontsize=12, color="r")
-ax.text(0.12, min(reynoldsb), "anisotropy=1", fontsize=12, color="r")
-fig.colorbar(scatter)
+ax.scatter(0.1, 1e-5, marker="o", s=2000 * (np.log10(2*0.5)**1) + 20, color="red")
+ax.scatter(0.1, 1e-6, marker="o", s=2000 * (np.log10(2*0.75)**1) + 20, color="red")
+ax.scatter(0.1, 2e-8, marker="o", s=2000 * (np.log10(2)**1) + 20, color="red")
+ax.text(0.16, 8e-6, r"$anisotropy = 0$", color="red", fontsize=14)
+ax.text(0.16, 5e-7, r"$anisotropy = 0.5$", color="red", fontsize=14)
+ax.text(0.16, 2e-8, r"$anisotropy = 1$", color="red", fontsize=14)
+divider = make_axes_locatable(ax)
+cax = divider.append_axes("right", size="5%", pad=0.05)
+cax.tick_params(labelsize=14)
+fig.colorbar(scatter, cax=cax)
 
 # Legend
 import matplotlib.lines as mlines
@@ -91,5 +111,5 @@
 blue_star = mlines.Line2D([], [], color='red', marker='o', linestyle='None',
                           markersize=8, label=r'$n_x = 960$')
 red_square = mlines.Line2D([], [], color='red', marker='s', linestyle='None',
-                          markersize=8, label=r'$n_x = 1920$')
+                           markersize=8, label=r'$n_x = 1920$')
 purple_triangle = mlines.Line2D([], [], color='red', marker='^', linestyle='None',
@@ -95,5 +115,5 @@
 purple_triangle = mlines.Line2D([], [], color='red', marker='^', linestyle='None',
-                          markersize=8, label=r'$n_x = 3840$')
+                                markersize=8, label=r'$n_x = 3840$')
 
 ax.legend(handles=[blue_star, red_square, purple_triangle],
           loc="upper center",
@@ -101,5 +121,12 @@
           borderaxespad=0.,
           ncol=len(markers),
           handletextpad=0.1,
-          fontsize=12)
+          fontsize=14)
+
+fig.tight_layout(pad=0.4)
+
+if SAVE:
+    path_save = "/home/users/calpelin7m/Phd/docs/EFMC18/figures"
+    fig.savefig(path_save + "/reynoldsb_froude.eps", format="eps",
+                bbox_inches="tight")
 plt.show()