diff --git a/Makefile b/Makefile
index 7c7294a126136b4f385c1e688d13398dbac7ec44_TWFrZWZpbGU=..204e59ec7e8081cc658c07dd10441e7771db33c4_TWFrZWZpbGU= 100644
--- a/Makefile
+++ b/Makefile
@@ -11,10 +11,10 @@
 	find fluidsim -name "__pycache__" -type d | xargs rm -rf
 
 cleanpythran:
-	find fluidsim -type d -name _pythran | xargs rm -rf
+	find fluidsim -type d -name __pythran__ | xargs rm -rf
 	find fluidsim -name "*pythran*.cpp" -delete
 	find fluidsim -name "*pythran*.so" -delete
 
 clean:
 	rm -rf build
 
@@ -15,10 +15,10 @@
 	find fluidsim -name "*pythran*.cpp" -delete
 	find fluidsim -name "*pythran*.so" -delete
 
 clean:
 	rm -rf build
 
-cleanall: clean clean_so
+cleanall: clean clean_so cleanpythran
 
 black:
 	black -l 82 fluidsim
diff --git a/README.rst b/README.rst
index 7c7294a126136b4f385c1e688d13398dbac7ec44_UkVBRE1FLnJzdA==..204e59ec7e8081cc658c07dd10441e7771db33c4_UkVBRE1FLnJzdA== 100644
--- a/README.rst
+++ b/README.rst
@@ -45,7 +45,7 @@
 
 For a **basic installation** it should be sufficient to run::
 
-  pip install fluiddyn fluidpythran [--user]
+  pip install fluiddyn [--user]
   pip install fluidsim [--user]
 
 Add ``--user`` flag if you are installing without setting up a virtual
@@ -66,7 +66,7 @@
 4. `A C++11 compiler and BLAS
    libraries <https://github.com/serge-sans-paille/pythran#installation>`__ and
    `configure
-   <https://fluidsim.readthedocs.io/en/latest/install.html#dependencies>`__ 
+   <https://fluidsim.readthedocs.io/en/latest/install.html#dependencies>`__
    ``~/.pythranrc`` to customize compilation of Pythran extensions
 5. `Configure
    <https://fluidfft.readthedocs.io/en/latest/install.html#basic-installation-with-pip>`__
@@ -130,4 +130,3 @@
 french license.
 
 .. _CeCILL: http://www.cecill.info/index.en.html
-
diff --git a/config.py b/config.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Y29uZmlnLnB5..204e59ec7e8081cc658c07dd10441e7771db33c4_Y29uZmlnLnB5 100644
--- a/config.py
+++ b/config.py
@@ -11,8 +11,6 @@
 
 """
 
-from __future__ import print_function
-
 import os
 import sys
 import subprocess
@@ -28,10 +26,9 @@
 except ImportError:
     from distutils.command.build_ext import build_ext
     from setuptools import Extension
-try:  # python 3
-    from configparser import ConfigParser
-except ImportError:  # python 2.7
-    from ConfigParser import ConfigParser
+
+from configparser import ConfigParser
+
 
 try:
     import colorlog as logging
@@ -222,19 +219,6 @@
 def get_default_config():
     """Generate default configuration."""
     config = ConfigParser()
-    section = "exclude_pythran"
-    config.add_section(section)
-    exclude_pythran = (
-        "fluidsim.solvers.plate2d",
-        "fluidsim.solvers.plate2d.output",
-        "fluidsim.solvers.sw1l",
-        "fluidsim.solvers.sw1l.output",
-        "fluidsim.solvers.ns2d.strat",
-        "fluidsim.solvers.ns2d.bouss",
-    )
-    for excluded in exclude_pythran:
-        config.set(section, excluded, "True")
-
     config.add_section("environ")
 
     return config
diff --git a/fluidsim/__init__.py b/fluidsim/__init__.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vX19pbml0X18ucHk=..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vX19pbml0X18ucHk= 100644
--- a/fluidsim/__init__.py
+++ b/fluidsim/__init__.py
@@ -25,6 +25,39 @@
 """
 
 from pathlib import Path
+import os
+import sys
+
+if "FLUIDSIM_PATH" in os.environ:
+    os.environ["FLUIDPYTHRAN_DIR"] = str(
+        Path(os.environ["FLUIDSIM_PATH"] / ".fluidpythran")
+    )
+
+if any(
+    any(test_tool in arg for arg in sys.argv)
+    for test_tool in ("pytest", "unittest", "fluidsim-test")
+):
+    from fluiddyn.util import mpi
+
+    mpi.printby0(
+        "Fluidsim guesses that it is tested so it "
+        "loads the Agg Matplotlib backend."
+    )
+    import matplotlib
+
+    matplotlib.use("Agg")
+    import matplotlib.pyplot as plt
+
+    def _show(*args, **kwargs):
+        pass
+
+    plt.show = _show
+
+    if "FLUID_COMPILE_CACHEDJIT" not in os.environ:
+        mpi.printby0("Compilation of cachedjit functions disabled.")
+        from fluidpythran import set_compile_cachedjit
+
+        set_compile_cachedjit(False)
 
 from ._version import __version__, get_local_version
 
diff --git a/fluidsim/base/output/increments.py b/fluidsim/base/output/increments.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvaW5jcmVtZW50cy5weQ==..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvaW5jcmVtZW50cy5weQ== 100644
--- a/fluidsim/base/output/increments.py
+++ b/fluidsim/base/output/increments.py
@@ -1,7 +1,4 @@
-from __future__ import division, print_function
-
-from builtins import range
 import h5py
 import os
 import numpy as np
 
@@ -4,7 +1,10 @@
 import h5py
 import os
 import numpy as np
 
+# pythran import numpy as np
+
+from fluidpythran import pythran_def, Array
 from fluiddyn.util import mpi
 
 from .base import SpecificOutput
@@ -8,7 +8,28 @@
 from fluiddyn.util import mpi
 
 from .base import SpecificOutput
-from .util_pythran import strfunc_from_pdf
+
+A = Array[float, "2d"]
+# pythran def strfunc_from_pdf(
+#     int32[], float64[][], float64[][], float, bool)
+
+
+@pythran_def
+def strfunc_from_pdf(
+    rxs: A, pdf: A, values: A, order: float, absolute: bool = False
+):
+    """Compute structure function of specified order from pdf for increments
+    module.
+
+    """
+    S_order = np.empty(rxs.shape)
+    if absolute:
+        values = abs(values)
+    for irx in range(rxs.size):
+        deltainc = abs(values[irx, 1] - values[irx, 0])
+        S_order[irx] = deltainc * np.sum(pdf[irx] * values[irx] ** order)
+
+    return S_order
 
 
 class Increments(SpecificOutput):
diff --git a/fluidsim/base/output/phys_fields.py b/fluidsim/base/output/phys_fields.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMucHk=..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMucHk= 100644
--- a/fluidsim/base/output/phys_fields.py
+++ b/fluidsim/base/output/phys_fields.py
@@ -12,8 +12,6 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
 from builtins import str
 
diff --git a/fluidsim/base/output/util_pythran.py b/fluidsim/base/output/util_pythran.py
deleted file mode 100644
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvdXRpbF9weXRocmFuLnB5..0000000000000000000000000000000000000000
--- a/fluidsim/base/output/util_pythran.py
+++ /dev/null
@@ -1,23 +0,0 @@
-# pythran export strfunc_from_pdf(
-#     float64[][], float64[][], float64[][], float, bool)
-
-# pythran export strfunc_from_pdf(
-#     int32[], float64[][], float64[][], float, bool)
-
-
-import numpy as np
-
-
-def strfunc_from_pdf(rxs, pdf, values, order, absolute=False):
-    """Compute structure function of specified order from pdf for increments
-    module.
-
-    """
-    S_order = np.empty(rxs.shape)
-    if absolute:
-        values = abs(values)
-    for irx in range(rxs.size):
-        deltainc = abs(values[irx, 1] - values[irx, 0])
-        S_order[irx] = deltainc * np.sum(pdf[irx] * values[irx] ** order)
-
-    return S_order
diff --git a/fluidsim/base/time_stepping/pseudo_spect.py b/fluidsim/base/time_stepping/pseudo_spect.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vYmFzZS90aW1lX3N0ZXBwaW5nL3BzZXVkb19zcGVjdC5weQ==..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vYmFzZS90aW1lX3N0ZXBwaW5nL3BzZXVkb19zcGVjdC5weQ== 100644
--- a/fluidsim/base/time_stepping/pseudo_spect.py
+++ b/fluidsim/base/time_stepping/pseudo_spect.py
@@ -29,8 +29,6 @@
 
 from .base import TimeSteppingBase
 
-from .rk_pythran import step0_RK2_pythran, step1_RK2_pythran
-
 fp = FluidPythran()
 
 T = Type(np.complex128)
@@ -242,23 +240,6 @@
             state_spect * diss + dt * diss2 * tendencies_n12
         )
 
-    def _time_step_RK2_pythran(self):
-        dt = self.deltat
-        diss, diss2 = self.exact_linear_coefs.get_updated_coefs()
-
-        compute_tendencies = self.sim.tendencies_nonlin
-        state_spect = self.sim.state.state_spect
-
-        tendencies_n = compute_tendencies()
-
-        state_spect_n12 = self._state_spect_tmp
-        # state_spect_n12[:] = (state_spect + dt / 2 * tendencies_n) * diss2
-        step0_RK2_pythran(state_spect_n12, state_spect, tendencies_n, diss2, dt)
-
-        tendencies_n12 = compute_tendencies(state_spect_n12, old=tendencies_n)
-        # state_spect[:] = state_spect * diss + dt * diss2 * tendencies_n12
-        step1_RK2_pythran(state_spect, tendencies_n12, diss, diss2, dt)
-
     def _time_step_RK2_fluidpythran(self):
         dt = self.deltat
         diss, diss2 = self.exact_linear_coefs.get_updated_coefs()
diff --git a/fluidsim/base/time_stepping/rk_pythran.py b/fluidsim/base/time_stepping/rk_pythran.py
deleted file mode 100644
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vYmFzZS90aW1lX3N0ZXBwaW5nL3JrX3B5dGhyYW4ucHk=..0000000000000000000000000000000000000000
--- a/fluidsim/base/time_stepping/rk_pythran.py
+++ /dev/null
@@ -1,40 +0,0 @@
-# pythran export step0_RK2_pythran(
-#     complex128[][][],
-#     complex128[][][],
-#     complex128[][][],
-#     float64[][] or complex128[][],
-#     float
-# )
-
-# pythran export step0_RK2_pythran(
-#     complex128[][][][],
-#     complex128[][][][],
-#     complex128[][][][],
-#     float64[][][] or complex128[][][],
-#     float
-# )
-
-
-def step0_RK2_pythran(state_spect_n12, state_spect, tendencies_n, diss2, dt):
-    state_spect_n12[:] = (state_spect + dt / 2 * tendencies_n) * diss2
-
-
-# pythran export step1_RK2_pythran(
-# complex128[][][],
-# complex128[][][],
-# float64[][] or complex128[][],
-# float64[][] or complex128[][],
-# float
-# )
-
-# pythran export step1_RK2_pythran(
-# complex128[][][][],
-# complex128[][][][],
-# float64[][][] or complex128[][][],
-# float64[][][] or complex128[][][],
-# float
-# )
-
-
-def step1_RK2_pythran(state_spect, tendencies_n12, diss, diss2, dt):
-    state_spect[:] = state_spect * diss + dt * diss2 * tendencies_n12
diff --git a/fluidsim/operators/operators2d.py b/fluidsim/operators/operators2d.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczJkLnB5..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczJkLnB5 100644
--- a/fluidsim/operators/operators2d.py
+++ b/fluidsim/operators/operators2d.py
@@ -14,7 +14,9 @@
 
 import numpy as np
 
+from fluidpythran import pythran_def, Array, FluidPythran
+
 from fluiddyn.util import mpi
 
 from fluidfft.fft2d.operators import OperatorsPseudoSpectral2D as _Operators
 
@@ -17,13 +19,6 @@
 from fluiddyn.util import mpi
 
 from fluidfft.fft2d.operators import OperatorsPseudoSpectral2D as _Operators
 
-from . import util2d_pythran
-from .util2d_pythran import (
-    dealiasing_setofvar,
-    laplacian_fft,
-    invlaplacian_fft,
-    compute_increments_dim1,
-)
 from ..base.setofvariables import SetOfVariables
 
@@ -28,4 +23,6 @@
 from ..base.setofvariables import SetOfVariables
 
-if not hasattr(util2d_pythran, "__pythran__"):
+fp = FluidPythran()
+
+if not fp.is_transpiling and not fp.is_compiled:
     warn(
@@ -31,5 +28,5 @@
     warn(
-        "util2d_pythran has to be pythranized to be efficient! "
+        "operators2d.py has to be pythranized to be efficient! "
         "Install pythran and recompile."
     )
 
@@ -33,6 +30,47 @@
         "Install pythran and recompile."
     )
 
+
+Af = Array[np.float64, "2d"]
+Ac = Array[np.complex128, "2d"]
+
+
+@pythran_def
+def dealiasing_setofvar(setofvar_fft: "complex128[][][]", where: "uint8[][]", n0: int, n1: int):
+    """Dealiasing of a setofvar arrays."""
+    nk = setofvar_fft.shape[0]
+
+    for i0 in range(n0):
+        for i1 in range(n1):
+            if where[i0, i1]:
+                for ik in range(nk):
+                    setofvar_fft[ik, i0, i1] = 0.0
+
+
+@pythran_def
+def laplacian_fft(a_fft: Ac, Kn: Af):
+    """Compute the n-th order Laplacian."""
+    return a_fft * Kn
+
+
+@pythran_def
+def invlaplacian_fft(a_fft: Ac, Kn_not0: Af, rank: int):
+    """Compute the n-th order inverse Laplacian."""
+    invlap_afft = a_fft / Kn_not0
+    if rank == 0:
+        invlap_afft[0, 0] = 0.0
+    return invlap_afft
+
+
+@pythran_def
+def compute_increments_dim1(var: Af, irx: int):
+    """Compute the increments of var over the dim 1."""
+    n1 = var.shape[1]
+    n1new = n1 - irx
+    inc_var = var[:, irx:n1] - var[:, 0:n1new]
+    return inc_var
+
+
 nb_proc = mpi.nb_proc
 rank = mpi.rank
 if nb_proc > 1:
diff --git a/fluidsim/operators/operators3d.py b/fluidsim/operators/operators3d.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczNkLnB5..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczNkLnB5 100644
--- a/fluidsim/operators/operators3d.py
+++ b/fluidsim/operators/operators3d.py
@@ -18,6 +18,7 @@
 
 import numpy as np
 
+from fluidpythran import pythran_def, Array, FluidPythran
 from fluiddyn.util.mpi import nb_proc, rank
 
 from fluidsim.base.setofvariables import SetOfVariables
@@ -26,6 +27,34 @@
 
 from .operators2d import OperatorsPseudoSpectral2D as OpPseudoSpectral2D
 
-from . import util3d_pythran
+fp = FluidPythran()
+
+Asov = Array[np.complex128, "4d"]
+Aui8 = Array[np.uint8, "3d"]
+Ac = Array[np.complex128, "3d"]
+Af = Array[np.float64, "3d"]
+
+@pythran_def
+def dealiasing_setofvar(sov: Asov, where_dealiased: Aui8):
+    """Dealiasing 3d setofvar object.
+
+    Parameters
+    ----------
+
+    sov : 4d ndarray
+        A set of variables array.
+
+    where_dealiased : 3d ndarray
+        A 3d array of "booleans" (actually uint8).
+
+    """
+    nk, n0, n1, n2 = sov.shape
+
+    for i0 in range(n0):
+        for i1 in range(n1):
+            for i2 in range(n2):
+                if where_dealiased[i0, i1, i2]:
+                    for ik in range(nk):
+                        sov[ik, i0, i1, i2] = 0.0
 
 
@@ -30,7 +59,19 @@
 
 
-def dealiasing_setofvar_numpy(sov, where_dealiased):
+@pythran_def
+def dealiasing_variable(ff_fft: Ac, where_dealiased: Aui8):
+    """Dealiasing 3d array"""
+    n0, n1, n2 = ff_fft.shape
+
+    for i0 in range(n0):
+        for i1 in range(n1):
+            for i2 in range(n2):
+                if where_dealiased[i0, i1, i2]:
+                    ff_fft[i0, i1, i2] = 0.0
+
+
+def dealiasing_setofvar_numpy(sov: Asov, where_dealiased: Aui8):
     for i in range(sov.shape[0]):
         sov[i][np.nonzero(where_dealiased)] = 0.0
 
 
@@ -33,8 +74,8 @@
     for i in range(sov.shape[0]):
         sov[i][np.nonzero(where_dealiased)] = 0.0
 
 
-def dealiasing_variable_numpy(ff_fft, where_dealiased):
+def dealiasing_variable_numpy(ff_fft: Ac, where_dealiased: Aui8):
     ff_fft[np.nonzero(where_dealiased)] = 0.0
 
 
@@ -38,11 +79,8 @@
     ff_fft[np.nonzero(where_dealiased)] = 0.0
 
 
-if hasattr(util3d_pythran, "__pythran__"):
-    dealiasing_variable = util3d_pythran.dealiasing_variable
-    dealiasing_setofvar = util3d_pythran.dealiasing_setofvar
-else:
+if not fp.is_transpiling and not fp.is_compiled:
     dealiasing_variable = dealiasing_variable_numpy
     dealiasing_setofvar = dealiasing_setofvar_numpy
 
 
@@ -45,7 +83,22 @@
     dealiasing_variable = dealiasing_variable_numpy
     dealiasing_setofvar = dealiasing_setofvar_numpy
 
 
+@pythran_def
+def urudfft_from_vxvyfft(vx_fft: Ac, vy_fft: Ac, kx: Af, ky: Af, rank: int):
+    k2 = kx ** 2 + ky ** 2
+    k2[k2 == 0.0] = 1e-10
+
+    divh_fft = 1j * (kx * vx_fft + ky * vy_fft)
+    urx_fft = vx_fft - divh_fft * kx / k2
+    ury_fft = vy_fft - divh_fft * ky / k2
+
+    udx_fft = vx_fft - urx_fft
+    udy_fft = vy_fft - ury_fft
+
+    return urx_fft, ury_fft, udx_fft, udy_fft
+
+
 class OperatorsPseudoSpectral3D(_Operators):
     """Provides fast Fourier transform functions and 3D operators.
 
@@ -241,7 +294,7 @@
         """Compute toroidal and poloidal horizontal velocities
 
         """
-        return util3d_pythran.urudfft_from_vxvyfft(
+        return urudfft_from_vxvyfft(
             vx_fft, vy_fft, self.Kx, self.Ky, rank
         )
 
diff --git a/fluidsim/operators/util2d_pythran.py b/fluidsim/operators/util2d_pythran.py
deleted file mode 100644
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vb3BlcmF0b3JzL3V0aWwyZF9weXRocmFuLnB5..0000000000000000000000000000000000000000
--- a/fluidsim/operators/util2d_pythran.py
+++ /dev/null
@@ -1,58 +0,0 @@
-"""
-Pythran compatible functions: 2d operators (:mod:`fluidsim.operators.util2d_pythran`)
-=====================================================================================
-
-.. autofunction:: dealiasing_setofvar
-
-.. autofunction:: laplacian_fft
-
-.. autofunction:: invlaplacian_fft
-
-.. autofunction:: compute_increments_dim1
-
-"""
-
-# import numpy as np
-
-# pythran export dealiasing_setofvar(complex128[][][], uint8[][], int, int)
-
-
-def dealiasing_setofvar(setofvar_fft, where, n0, n1):
-    """Dealiasing of a setofvar arrays."""
-    nk = setofvar_fft.shape[0]
-
-    for i0 in range(n0):
-        for i1 in range(n1):
-            if where[i0, i1]:
-                for ik in range(nk):
-                    setofvar_fft[ik, i0, i1] = 0.0
-
-
-# pythran export laplacian_fft(complex128[][], float64[][])
-
-
-def laplacian_fft(a_fft, Kn):
-    """Compute the n-th order Laplacian."""
-    return a_fft * Kn
-
-
-# pythran export invlaplacian_fft(complex128[][], float64[][], int)
-
-
-def invlaplacian_fft(a_fft, Kn_not0, rank):
-    """Compute the n-th order inverse Laplacian."""
-    invlap_afft = a_fft / Kn_not0
-    if rank == 0:
-        invlap_afft[0, 0] = 0.0
-    return invlap_afft
-
-
-# pythran export compute_increments_dim1(float64[][], int)
-
-
-def compute_increments_dim1(var, irx):
-    """Compute the increments of var over the dim 1."""
-    n1 = var.shape[1]
-    n1new = n1 - irx
-    inc_var = var[:, irx:n1] - var[:, 0:n1new]
-    return inc_var
diff --git a/fluidsim/operators/util3d_pythran.py b/fluidsim/operators/util3d_pythran.py
deleted file mode 100644
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vb3BlcmF0b3JzL3V0aWwzZF9weXRocmFuLnB5..0000000000000000000000000000000000000000
--- a/fluidsim/operators/util3d_pythran.py
+++ /dev/null
@@ -1,68 +0,0 @@
-"""
-Pythran compatible functions: 3d operators (:mod:`fluidsim.operators.util3d_pythran`)
-=====================================================================================
-
-.. autofunction:: dealiasing_setofvar
-
-.. autofunction:: dealiasing_variable
-
-.. autofunction:: urudfft_from_vxvyfft
-
-"""
-
-# pythran export dealiasing_setofvar(complex128[][][][], uint8[][][])
-
-
-def dealiasing_setofvar(sov, where_dealiased):
-    """Dealiasing 3d setofvar object.
-
-    Parameters
-    ----------
-
-    sov : 4d ndarray
-        A set of variables array.
-
-    where_dealiased : 3d ndarray
-        A 3d array of "booleans" (actually uint8).
-
-    """
-    nk, n0, n1, n2 = sov.shape
-
-    for i0 in range(n0):
-        for i1 in range(n1):
-            for i2 in range(n2):
-                if where_dealiased[i0, i1, i2]:
-                    for ik in range(nk):
-                        sov[ik, i0, i1, i2] = 0.0
-
-
-# pythran export dealiasing_variable(complex128[][][], uint8[][][])
-
-
-def dealiasing_variable(ff_fft, where_dealiased):
-    """Dealiasing 3d array"""
-    n0, n1, n2 = ff_fft.shape
-
-    for i0 in range(n0):
-        for i1 in range(n1):
-            for i2 in range(n2):
-                if where_dealiased[i0, i1, i2]:
-                    ff_fft[i0, i1, i2] = 0.0
-
-
-# pythran export urudfft_from_vxvyfft(
-#     complex128[][][], complex128[][][], float64[][][], float64[][][], int)
-
-
-def urudfft_from_vxvyfft(vx_fft, vy_fft, kx, ky, rank):
-    k2 = kx ** 2 + ky ** 2
-    k2[k2 == 0.0] = 1e-10
-
-    divh_fft = 1j * (kx * vx_fft + ky * vy_fft)
-    urx_fft = vx_fft - divh_fft * kx / k2
-    ury_fft = vy_fft - divh_fft * ky / k2
-
-    udx_fft = vx_fft - urx_fft
-    udy_fft = vy_fft - ury_fft
-
-    return urx_fft, ury_fft, udx_fft, udy_fft
diff --git a/fluidsim/solvers/ns2d/bouss/solver.py b/fluidsim/solvers/ns2d/bouss/solver.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9uczJkL2JvdXNzL3NvbHZlci5weQ==..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vc29sdmVycy9uczJkL2JvdXNzL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ns2d/bouss/solver.py
+++ b/fluidsim/solvers/ns2d/bouss/solver.py
@@ -6,8 +6,5 @@
    :private-members:
 
 """
-
-from __future__ import division
-
 import numpy as np
 
@@ -12,3 +9,5 @@
 import numpy as np
 
+from fluidpythran import cachedjit, Array
+
 from fluidsim.base.setofvariables import SetOfVariables
@@ -14,4 +13,3 @@
 from fluidsim.base.setofvariables import SetOfVariables
-
 from fluidsim.solvers.ns2d.solver import InfoSolverNS2D, Simul as SimulNS2D
 
@@ -16,6 +14,16 @@
 from fluidsim.solvers.ns2d.solver import InfoSolverNS2D, Simul as SimulNS2D
 
-from .util_pythran import tendencies_nonlin_ns2dbouss
+
+AF = Array[np.float64, "2d"]
+
+
+@cachedjit
+def tendencies_nonlin_ns2dbouss(
+    ux: AF, uy: AF, px_rot: AF, py_rot: AF, px_b: AF, py_b: AF
+):
+    Frot = -ux * px_rot - uy * py_rot + px_b
+    Fb = -ux * px_b - uy * py_b
+    return Frot, Fb
 
 
 class InfoSolverNS2DBouss(InfoSolverNS2D):
diff --git a/fluidsim/solvers/ns2d/bouss/util_pythran.py b/fluidsim/solvers/ns2d/bouss/util_pythran.py
deleted file mode 100644
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9uczJkL2JvdXNzL3V0aWxfcHl0aHJhbi5weQ==..0000000000000000000000000000000000000000
--- a/fluidsim/solvers/ns2d/bouss/util_pythran.py
+++ /dev/null
@@ -1,11 +0,0 @@
-# pythran export tendencies_nonlin_ns2dbouss(
-#     float64[][], float64[][], float64[][], float64[][],
-#     float64[][], float64[][])
-
-
-def tendencies_nonlin_ns2dbouss(ux, uy, px_rot, py_rot, px_b, py_b):
-
-    Frot = -ux * px_rot - uy * py_rot + px_b
-    Fb = -ux * px_b - uy * py_b
-
-    return Frot, Fb
diff --git a/fluidsim/solvers/ns2d/solver.py b/fluidsim/solvers/ns2d/solver.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3NvbHZlci5weQ==..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vc29sdmVycy9uczJkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ns2d/solver.py
+++ b/fluidsim/solvers/ns2d/solver.py
@@ -13,7 +13,6 @@
    :private-members:
 
 """
-from __future__ import division
 
 import sys
 
@@ -17,6 +16,10 @@
 
 import sys
 
+import numpy as np
+
+from fluidpythran import pythran_def, Array
+
 from fluidsim.base.setofvariables import SetOfVariables
 
 from fluidsim.base.solvers.pseudo_spect import (
@@ -24,7 +27,16 @@
     InfoSolverPseudoSpectral,
 )
 
-from fluidsim.solvers.ns2d.util_pythran import compute_Frot
+Af = Array[np.float64, "2d"]
+
+
+@pythran_def
+def compute_Frot(ux: Af, uy: Af, px_rot: Af, py_rot: Af, beta: float = 0):
+    if beta == 0:
+        return -ux * px_rot - uy * py_rot
+
+    else:
+        return -ux * px_rot - uy * (py_rot + beta)
 
 
 class InfoSolverNS2D(InfoSolverPseudoSpectral):
diff --git a/fluidsim/solvers/ns2d/strat/solver.py b/fluidsim/solvers/ns2d/strat/solver.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3NvbHZlci5weQ==..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ns2d/strat/solver.py
+++ b/fluidsim/solvers/ns2d/strat/solver.py
@@ -6,7 +6,6 @@
    :private-members:
 
 """
-from __future__ import division
 
 import numpy as np
 
@@ -10,7 +9,9 @@
 
 import numpy as np
 
+from fluidpythran import cachedjit, Array
+
 from fluidsim.base.setofvariables import SetOfVariables
 
 from fluidsim.solvers.ns2d.solver import InfoSolverNS2D, Simul as SimulNS2D
 
@@ -13,8 +14,17 @@
 from fluidsim.base.setofvariables import SetOfVariables
 
 from fluidsim.solvers.ns2d.solver import InfoSolverNS2D, Simul as SimulNS2D
 
-from .util_pythran import tendencies_nonlin_ns2dstrat
+AF = Array[np.float64, "2d"]
+
+
+@cachedjit
+def tendencies_nonlin_ns2dstrat(
+    ux: AF, uy: AF, px_rot: AF, py_rot: AF, px_b: AF, py_b: AF, N: float
+):
+    Frot = -ux * px_rot - uy * py_rot
+    Fb = -ux * px_b - uy * py_b - N ** 2 * uy
+    return Frot, Fb
 
 
 class InfoSolverNS2DStrat(InfoSolverNS2D):
diff --git a/fluidsim/solvers/ns2d/strat/util_pythran.py b/fluidsim/solvers/ns2d/strat/util_pythran.py
deleted file mode 100644
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3V0aWxfcHl0aHJhbi5weQ==..0000000000000000000000000000000000000000
--- a/fluidsim/solvers/ns2d/strat/util_pythran.py
+++ /dev/null
@@ -1,11 +0,0 @@
-# pythran export tendencies_nonlin_ns2dstrat(
-#     float64[][], float64[][], float64[][], float64[][],
-#     float64[][], float64[][], float)
-
-
-def tendencies_nonlin_ns2dstrat(ux, uy, px_rot, py_rot, px_b, py_b, N):
-
-    Frot = -ux * px_rot - uy * py_rot
-    Fb = -ux * px_b - uy * py_b - N ** 2 * uy
-
-    return Frot, Fb
diff --git a/fluidsim/solvers/ns2d/util_pythran.py b/fluidsim/solvers/ns2d/util_pythran.py
deleted file mode 100644
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3V0aWxfcHl0aHJhbi5weQ==..0000000000000000000000000000000000000000
--- a/fluidsim/solvers/ns2d/util_pythran.py
+++ /dev/null
@@ -1,10 +0,0 @@
-# pythran export compute_Frot(
-#     float64[][], float64[][], float64[][], float64[][], float)
-
-
-def compute_Frot(ux, uy, px_rot, py_rot, beta=0):
-    if beta == 0:
-        return -ux * px_rot - uy * py_rot
-
-    else:
-        return -ux * px_rot - uy * (py_rot + beta)
diff --git a/fluidsim/solvers/plate2d/operators.py b/fluidsim/solvers/plate2d/operators.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL29wZXJhdG9ycy5weQ==..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL29wZXJhdG9ycy5weQ== 100644
--- a/fluidsim/solvers/plate2d/operators.py
+++ b/fluidsim/solvers/plate2d/operators.py
@@ -9,4 +9,5 @@
 
 """
 
+import numpy as np
 
@@ -12,5 +13,5 @@
 
-import numpy as np
+from fluidpythran import cachedjit, Array
 
 from ...operators.operators2d import OperatorsPseudoSpectral2D
 
@@ -14,7 +15,26 @@
 
 from ...operators.operators2d import OperatorsPseudoSpectral2D
 
-from .util_oper_pythran import monge_ampere_step0, monge_ampere_step1
+AC2 = Array[np.complex128, "2d"]
+AF2 = Array[np.float64, "2d"]
+
+
+@cachedjit
+def monge_ampere_step0(a_fft: AC2, b_fft: AC2, KX2: AF2, KY2: AF2, KXKZ: AF2):
+    pxx_a_fft = -a_fft * KX2
+    pyy_a_fft = -a_fft * KY2
+    pxy_a_fft = -a_fft * KXKZ
+    pxx_b_fft = -b_fft * KX2
+    pyy_b_fft = -b_fft * KY2
+    pxy_b_fft = -b_fft * KXKZ
+    return pxx_a_fft, pyy_a_fft, pxy_a_fft, pxx_b_fft, pyy_b_fft, pxy_b_fft
+
+
+@cachedjit
+def monge_ampere_step1(
+    pxx_a: AF2, pyy_a: AF2, pxy_a: AF2, pxx_b: AF2, pyy_b: AF2, pxy_b: AF2
+):
+    return pxx_a * pyy_b + pyy_a * pxx_b - 2 * pxy_a * pxy_b
 
 
 class OperatorsPseudoSpectralPlate2D(OperatorsPseudoSpectral2D):
@@ -37,10 +57,6 @@
         self.tmp_pyy_b = np.empty_like(self.XX)
         self.tmp_pxy_b = np.empty_like(self.XX)
 
-    # def monge_ampere_from_fft2(self, a_fft, b_fft):
-    #     return monge_ampere_from_fft(
-    #         a_fft, b_fft, self.KX, self.KY, self.ifft2)
-
     def monge_ampere_from_fft(self, a_fft, b_fft):
         """Compute the Monge-Ampere operator"""
 
@@ -65,7 +81,7 @@
         )
 
 
-def monge_ampere_from_fft2(a_fft, b_fft, KX, KY, ifft2):
+def monge_ampere_from_fft_numpy(a_fft, b_fft, KX, KY, ifft2):
 
     pxx_a = -ifft2(a_fft * KX ** 2)
     pyy_a = -ifft2(a_fft * KY ** 2)
diff --git a/fluidsim/solvers/plate2d/output/correlations_freq.py b/fluidsim/solvers/plate2d/output/correlations_freq.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9jb3JyZWxhdGlvbnNfZnJlcS5weQ==..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9jb3JyZWxhdGlvbnNfZnJlcS5weQ== 100644
--- a/fluidsim/solvers/plate2d/output/correlations_freq.py
+++ b/fluidsim/solvers/plate2d/output/correlations_freq.py
@@ -17,6 +17,8 @@
 from builtins import range
 import os
 import numpy as np
+
+# pythran import numpy as np
 import h5py
 import matplotlib.pyplot as plt
 
@@ -20,8 +22,9 @@
 import h5py
 import matplotlib.pyplot as plt
 
+from fluidpythran import cachedjit
 from fluiddyn.util import mpi
 
 from fluidsim.base.output.base import SpecificOutput
 from fluiddyn.calcul.easypyfft import FFTW1DReal2Complex
 
@@ -23,7 +26,38 @@
 from fluiddyn.util import mpi
 
 from fluidsim.base.output.base import SpecificOutput
 from fluiddyn.calcul.easypyfft import FFTW1DReal2Complex
 
-# from fluidsim.operators.miscellaneous import compute_correl4, compute_correl2
+
+@cachedjit
+def compute_correl4_seq(
+    q_fftt: "complex128[][]", iomegas1: "int32[]", nb_omegas: int, nb_xs_seq: int
+):
+    r"""Compute the correlations 4.
+
+    .. math::
+       C_4(\omega_1, \omega_2, \omega_3, \omega_4) =
+       \langle
+       \tilde w(\omega_1, \mathbf{x})
+       \tilde w(\omega_2, \mathbf{x})
+       \tilde w(\omega_3, \mathbf{x})^*
+       \tilde w(\omega_4, \mathbf{x})^*
+       \rangle_\mathbf{x},
+
+    where
+
+    .. math::
+       \omega_2 = \omega_3 + \omega_4 - \omega_1
+
+    and :math:`\omega_1 > 0`, :math:`\omega_3 > 0` and
+    :math:`\omega_4 > 0`. Thus, this function produces an array
+    :math:`C_4(\omega_1, \omega_3, \omega_4)`.
+
+    """
+    q_fftt_conj = np.conj(q_fftt)
+
+    nx = q_fftt.shape[0]
+    n0 = iomegas1.shape[0]
+
+    corr4 = np.zeros((len(iomegas1), nb_omegas, nb_omegas), dtype=np.complex128)
 
@@ -29,5 +63,71 @@
 
-from .util_pythran import compute_correl2_seq, compute_correl4_seq
+    for i1 in range(n0):
+        io1 = iomegas1[i1]
+        # this loop could be parallelized (OMP)
+        for io3 in range(nb_omegas):
+            # we use the symmetry omega_3 <--> omega_4
+            for io4 in range(io3 + 1):
+                io2 = io3 + io4 - io1
+                if io2 < 0:
+                    io2 = -io2
+                    for ix in range(nx):
+                        corr4[i1, io3, io4] += (
+                            q_fftt[ix, io1]
+                            * q_fftt_conj[ix, io3]
+                            * q_fftt_conj[ix, io4]
+                            * q_fftt_conj[ix, io2]
+                        )
+                elif io2 >= nb_omegas:
+                    io2 = 2 * nb_omegas - 1 - io2
+                    for ix in range(nx):
+                        corr4[i1, io3, io4] += (
+                            q_fftt[ix, io1]
+                            * q_fftt_conj[ix, io3]
+                            * q_fftt_conj[ix, io4]
+                            * q_fftt_conj[ix, io2]
+                        )
+                else:
+                    for ix in range(nx):
+                        corr4[i1, io3, io4] += (
+                            q_fftt[ix, io1]
+                            * q_fftt_conj[ix, io3]
+                            * q_fftt_conj[ix, io4]
+                            * q_fftt[ix, io2]
+                        )
+                # symmetry omega_3 <--> omega_4:
+                corr4[i1, io4, io3] = corr4[i1, io3, io4]
+    return corr4
+
+
+@cachedjit
+def compute_correl2_seq(
+    q_fftt: "complex128[][]", iomegas1: "int32[]", nb_omegas: int, nb_xs_seq: int
+):
+    r"""Compute the correlations 2.
+
+    .. math::
+       C_2(\omega_1, \omega_2) =
+       \langle
+       \tilde w(\omega_1, \mathbf{x})
+       \tilde w(\omega_2, \mathbf{x})^*
+       \rangle_\mathbf{x},
+
+    where :math:`\omega_1 = \omega_2`. Thus, this function
+    produces an array :math:`C_2(\omega)`.
+
+    """
+    q_fftt_conj = np.conj(q_fftt)
+    nx = q_fftt.shape[0]
+
+    corr2 = np.zeros((nb_omegas, nb_omegas), dtype=np.complex128)
+
+    for io3 in range(nb_omegas):
+        for io4 in range(io3 + 1):
+            for ix in range(nx):
+                corr2[io3, io4] += q_fftt[ix, io3] * q_fftt_conj[ix, io4]
+            corr2[io4, io3] = np.conj(corr2[io3, io4])
+
+    return corr2
 
 
 def compute_correl4(q_fftt, iomegas1, nb_omegas, nb_xs_seq):
diff --git a/fluidsim/solvers/plate2d/output/util_pythran.py b/fluidsim/solvers/plate2d/output/util_pythran.py
deleted file mode 100644
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC91dGlsX3B5dGhyYW4ucHk=..0000000000000000000000000000000000000000
--- a/fluidsim/solvers/plate2d/output/util_pythran.py
+++ /dev/null
@@ -1,102 +0,0 @@
-import numpy as np
-
-
-# pythran export compute_correl4_seq(complex128[][], int32[], int, int)
-
-
-def compute_correl4_seq(q_fftt, iomegas1, nb_omegas, nb_xs_seq):
-    r"""Compute the correlations 4.
-
-    .. math::
-       C_4(\omega_1, \omega_2, \omega_3, \omega_4) =
-       \langle
-       \tilde w(\omega_1, \mathbf{x})
-       \tilde w(\omega_2, \mathbf{x})
-       \tilde w(\omega_3, \mathbf{x})^*
-       \tilde w(\omega_4, \mathbf{x})^*
-       \rangle_\mathbf{x},
-
-    where
-
-    .. math::
-       \omega_2 = \omega_3 + \omega_4 - \omega_1
-
-    and :math:`\omega_1 > 0`, :math:`\omega_3 > 0` and
-    :math:`\omega_4 > 0`. Thus, this function produces an array
-    :math:`C_4(\omega_1, \omega_3, \omega_4)`.
-
-    """
-    q_fftt_conj = np.conj(q_fftt)
-
-    nx = q_fftt.shape[0]
-    n0 = iomegas1.shape[0]
-
-    corr4 = np.zeros((len(iomegas1), nb_omegas, nb_omegas), dtype=np.complex128)
-
-    for i1 in range(n0):
-        io1 = iomegas1[i1]
-        # this loop could be parallelized (OMP)
-        for io3 in range(nb_omegas):
-            # we use the symmetry omega_3 <--> omega_4
-            for io4 in range(io3 + 1):
-                io2 = io3 + io4 - io1
-                if io2 < 0:
-                    io2 = -io2
-                    for ix in range(nx):
-                        corr4[i1, io3, io4] += (
-                            q_fftt[ix, io1]
-                            * q_fftt_conj[ix, io3]
-                            * q_fftt_conj[ix, io4]
-                            * q_fftt_conj[ix, io2]
-                        )
-                elif io2 >= nb_omegas:
-                    io2 = 2 * nb_omegas - 1 - io2
-                    for ix in range(nx):
-                        corr4[i1, io3, io4] += (
-                            q_fftt[ix, io1]
-                            * q_fftt_conj[ix, io3]
-                            * q_fftt_conj[ix, io4]
-                            * q_fftt_conj[ix, io2]
-                        )
-                else:
-                    for ix in range(nx):
-                        corr4[i1, io3, io4] += (
-                            q_fftt[ix, io1]
-                            * q_fftt_conj[ix, io3]
-                            * q_fftt_conj[ix, io4]
-                            * q_fftt[ix, io2]
-                        )
-                # symmetry omega_3 <--> omega_4:
-                corr4[i1, io4, io3] = corr4[i1, io3, io4]
-    return corr4
-
-
-# pythran export compute_correl2_seq(complex128[][], int32[], int, int)
-
-
-def compute_correl2_seq(q_fftt, iomegas1, nb_omegas, nb_xs_seq):
-    r"""Compute the correlations 2.
-
-    .. math::
-       C_2(\omega_1, \omega_2) =
-       \langle
-       \tilde w(\omega_1, \mathbf{x})
-       \tilde w(\omega_2, \mathbf{x})^*
-       \rangle_\mathbf{x},
-
-    where :math:`\omega_1 = \omega_2`. Thus, this function
-    produces an array :math:`C_2(\omega)`.
-
-    """
-    q_fftt_conj = np.conj(q_fftt)
-    nx = q_fftt.shape[0]
-
-    corr2 = np.zeros((nb_omegas, nb_omegas), dtype=np.complex128)
-
-    for io3 in range(nb_omegas):
-        for io4 in range(io3 + 1):
-            for ix in range(nx):
-                corr2[io3, io4] += q_fftt[ix, io3] * q_fftt_conj[ix, io4]
-            corr2[io4, io3] = np.conj(corr2[io3, io4])
-
-    return corr2
diff --git a/fluidsim/solvers/plate2d/util_oper_pythran.py b/fluidsim/solvers/plate2d/util_oper_pythran.py
deleted file mode 100644
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL3V0aWxfb3Blcl9weXRocmFuLnB5..0000000000000000000000000000000000000000
--- a/fluidsim/solvers/plate2d/util_oper_pythran.py
+++ /dev/null
@@ -1,22 +0,0 @@
-# pythran export monge_ampere_step0(
-#     complex128[][], complex128[][],
-#     float64[][], float64[][], float64[][])
-
-
-def monge_ampere_step0(a_fft, b_fft, KX2, KY2, KXKZ):
-    pxx_a_fft = -a_fft * KX2
-    pyy_a_fft = -a_fft * KY2
-    pxy_a_fft = -a_fft * KXKZ
-    pxx_b_fft = -b_fft * KX2
-    pyy_b_fft = -b_fft * KY2
-    pxy_b_fft = -b_fft * KXKZ
-    return pxx_a_fft, pyy_a_fft, pxy_a_fft, pxx_b_fft, pyy_b_fft, pxy_b_fft
-
-
-# pythran export monge_ampere_step1(
-#     float64[][], float64[][], float64[][],
-#     float64[][], float64[][], float64[][])
-
-
-def monge_ampere_step1(pxx_a, pyy_a, pxy_a, pxx_b, pyy_b, pxy_b):
-    return pxx_a * pyy_b + pyy_a * pxx_b - 2 * pxy_a * pxy_b
diff --git a/fluidsim/solvers/sw1l/operators.py b/fluidsim/solvers/sw1l/operators.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL29wZXJhdG9ycy5weQ==..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vc29sdmVycy9zdzFsL29wZXJhdG9ycy5weQ== 100644
--- a/fluidsim/solvers/sw1l/operators.py
+++ b/fluidsim/solvers/sw1l/operators.py
@@ -9,5 +9,11 @@
 
 """
 
+import numpy as np
+
+# pythran import numpy as np
+
+from fluidpythran import cachedjit, Array
+
 from fluidsim.operators.operators2d import OperatorsPseudoSpectral2D, rank
 
@@ -12,6 +18,90 @@
 from fluidsim.operators.operators2d import OperatorsPseudoSpectral2D, rank
 
-from . import util_oper_pythran
+AC = Array[np.complex128, "2d"]
+AF = Array[np.float64, "2d"]
+
+
+@cachedjit
+def _qapamfft_from_uxuyetafft(
+    ux_fft: AC,
+    uy_fft: AC,
+    eta_fft: AC,
+    n0: int,
+    n1: int,
+    KX: AF,
+    KY: AF,
+    K2: AF,
+    Kappa_over_ic: AC,
+    f: float,
+    c2: float,
+    rank: int,
+):
+    """Calculate normal modes from primitive variables."""
+    freq_Corio = f
+    f_over_c2 = freq_Corio / c2
+
+    q_fft = np.empty([n0, n1], dtype=np.complex128)
+    ap_fft = np.empty([n0, n1], dtype=np.complex128)
+    am_fft = np.empty([n0, n1], dtype=np.complex128)
+
+    if freq_Corio != 0:
+        for i0 in range(n0):
+            for i1 in range(n1):
+                if i0 == 0 and i1 == 0 and rank == 0:
+                    q_fft[i0, i1] = 0
+                    ap_fft[i0, i1] = ux_fft[0, 0] + 1.0j * uy_fft[0, 0]
+                    am_fft[i0, i1] = ux_fft[0, 0] - 1.0j * uy_fft[0, 0]
+                else:
+
+                    rot_fft = 1j * (
+                        KX[i0, i1] * uy_fft[i0, i1] - KY[i0, i1] * ux_fft[i0, i1]
+                    )
+
+                    q_fft[i0, i1] = rot_fft - freq_Corio * eta_fft[i0, i1]
+
+                    a_over2_fft = 0.5 * (
+                        K2[i0, i1] * eta_fft[i0, i1] + f_over_c2 * rot_fft
+                    )
+
+                    Deltaa_over2_fft = (
+                        0.5j
+                        * Kappa_over_ic[i0, i1]
+                        * (
+                            KX[i0, i1] * ux_fft[i0, i1]
+                            + KY[i0, i1] * uy_fft[i0, i1]
+                        )
+                    )
+
+                    ap_fft[i0, i1] = a_over2_fft + Deltaa_over2_fft
+                    am_fft[i0, i1] = a_over2_fft - Deltaa_over2_fft
+
+    else:  # (freq_Corio == 0.)
+        for i0 in range(n0):
+            for i1 in range(n1):
+                if i0 == 0 and i1 == 0 and rank == 0:
+                    q_fft[i0, i1] = 0
+                    ap_fft[i0, i1] = ux_fft[0, 0] + 1.0j * uy_fft[0, 0]
+                    am_fft[i0, i1] = ux_fft[0, 0] - 1.0j * uy_fft[0, 0]
+                else:
+                    q_fft[i0, i1] = 1j * (
+                        KX[i0, i1] * uy_fft[i0, i1] - KY[i0, i1] * ux_fft[i0, i1]
+                    )
+
+                    a_over2_fft = 0.5 * K2[i0, i1] * eta_fft[i0, i1]
+
+                    Deltaa_over2_fft = (
+                        0.5j
+                        * Kappa_over_ic[i0, i1]
+                        * (
+                            KX[i0, i1] * ux_fft[i0, i1]
+                            + KY[i0, i1] * uy_fft[i0, i1]
+                        )
+                    )
+
+                    ap_fft[i0, i1] = a_over2_fft + Deltaa_over2_fft
+                    am_fft[i0, i1] = a_over2_fft - Deltaa_over2_fft
+
+    return q_fft, ap_fft, am_fft
 
 
 class OperatorsPseudoSpectralSW1L(OperatorsPseudoSpectral2D):
@@ -31,7 +121,7 @@
         f = float(params.f)
         c2 = float(params.c2)
 
-        return util_oper_pythran.qapamfft_from_uxuyetafft(
+        return _qapamfft_from_uxuyetafft(
             ux_fft,
             uy_fft,
             eta_fft,
diff --git a/fluidsim/solvers/sw1l/output/__init__.py b/fluidsim/solvers/sw1l/output/__init__.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9fX2luaXRfXy5weQ==..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9fX2luaXRfXy5weQ== 100644
--- a/fluidsim/solvers/sw1l/output/__init__.py
+++ b/fluidsim/solvers/sw1l/output/__init__.py
@@ -15,7 +15,6 @@
    normal_mode
 
 """
-from __future__ import division
 
 import numpy as np
 
@@ -19,4 +18,6 @@
 
 import numpy as np
 
+from fluidpythran import cachedjit
+
 from fluidsim.base.output import OutputBasePseudoSpectral
@@ -22,5 +23,23 @@
 from fluidsim.base.output import OutputBasePseudoSpectral
-from .util_pythran import linear_eigenmode_from_values_1k
+
+
+@cachedjit
+def linear_eigenmode_from_values_1k(
+    ux_fft: np.complex128,
+    uy_fft: np.complex128,
+    eta_fft: np.complex128,
+    kx: float,
+    ky: float,
+    f: "float or int",
+    c2: "float or int",
+):
+    """Compute q, d, a (fft) for a single wavenumber."""
+    div_fft = 1j * (kx * ux_fft + ky * uy_fft)
+    rot_fft = 1j * (kx * uy_fft - ky * ux_fft)
+    q_fft = rot_fft - f * eta_fft
+    k2 = kx ** 2 + ky ** 2
+    ageo_fft = f * rot_fft / c2 + k2 * eta_fft
+    return q_fft, div_fft, ageo_fft
 
 
 class OutputBaseSW1L(OutputBasePseudoSpectral):
diff --git a/fluidsim/solvers/sw1l/output/normal_mode.py b/fluidsim/solvers/sw1l/output/normal_mode.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9ub3JtYWxfbW9kZS5weQ==..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9ub3JtYWxfbW9kZS5weQ== 100644
--- a/fluidsim/solvers/sw1l/output/normal_mode.py
+++ b/fluidsim/solvers/sw1l/output/normal_mode.py
@@ -20,4 +20,8 @@
 
 import numpy as np
 
+# pythran import numpy as np
+
+from fluidpythran import cachedjit, Array
+
 from fluiddyn.util import mpi
@@ -23,5 +27,41 @@
 from fluiddyn.util import mpi
-from .util_pythran import get_qmat
+
+
+AF = Array[float, "2d"]
+
+
+@cachedjit
+def get_qmat(
+    f: "float or int",
+    c: "float or int",
+    sigma: AF,
+    KX: AF,
+    KY: AF,
+    K: AF,
+    K2: AF,
+    K_not0: AF,
+):
+    """Compute Q matrix to transform q, ap, am (fft) -> b0, b+, b- (fft) with
+    dimensions of velocity.
+
+    """
+    ck = c * K_not0
+    qmat = np.array(
+        [
+            [
+                -1j * 2.0 ** 0.5 * ck * KY,
+                +1j * f * KY + KX * sigma,
+                +1j * f * KY - KX * sigma,
+            ],
+            [
+                +1j * 2.0 ** 0.5 * ck * KX,
+                -1j * f * KX + KY * sigma,
+                -1j * f * KX - KY * sigma,
+            ],
+            [2.0 ** 0.5 * f * K, c * K2, c * K2],
+        ]
+    ) / (2.0 ** 0.5 * sigma * K_not0)
+    return qmat
 
 
 class NormalModeBase(object):
@@ -79,7 +119,7 @@
         return bvec_fft
 
     def bvecfft_from_uxuyetafft(self, ux_fft, uy_fft, eta_fft):
-        r""" Compute normal mode vector, :math:`\mathbf{B}` 
+        r""" Compute normal mode vector, :math:`\mathbf{B}`
         with dimensions of velocity from primitive variables.  """
         q_fft, ap_fft, am_fft = self.oper.qapamfft_from_uxuyetafft(
             ux_fft, uy_fft, eta_fft
diff --git a/fluidsim/solvers/sw1l/output/util_pythran.py b/fluidsim/solvers/sw1l/output/util_pythran.py
deleted file mode 100644
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC91dGlsX3B5dGhyYW4ucHk=..0000000000000000000000000000000000000000
--- a/fluidsim/solvers/sw1l/output/util_pythran.py
+++ /dev/null
@@ -1,45 +0,0 @@
-import numpy as np
-
-
-# pythran export get_qmat(
-#     float or int, float or int, float64[][], float64[][], float64[][],
-#     float64[][], float64[][], float64[][])
-
-
-def get_qmat(f, c, sigma, KX, KY, K, K2, K_not0):
-    """Compute Q matrix to transform q, ap, am (fft) -> b0, b+, b- (fft) with
-    dimensions of velocity.
-
-    """
-    ck = c * K_not0
-    qmat = np.array(
-        [
-            [
-                -1j * 2.0 ** 0.5 * ck * KY,
-                +1j * f * KY + KX * sigma,
-                +1j * f * KY - KX * sigma,
-            ],
-            [
-                +1j * 2.0 ** 0.5 * ck * KX,
-                -1j * f * KX + KY * sigma,
-                -1j * f * KX - KY * sigma,
-            ],
-            [2.0 ** 0.5 * f * K, c * K2, c * K2],
-        ]
-    ) / (2.0 ** 0.5 * sigma * K_not0)
-    return qmat
-
-
-# pythran export linear_eigenmode_from_values_1k(
-#     complex128, complex128, complex128, float, float,
-#     float or int, float or int)
-
-
-def linear_eigenmode_from_values_1k(ux_fft, uy_fft, eta_fft, kx, ky, f, c2):
-    """Compute q, d, a (fft) for a single wavenumber."""
-    div_fft = 1j * (kx * ux_fft + ky * uy_fft)
-    rot_fft = 1j * (kx * uy_fft - ky * ux_fft)
-    q_fft = rot_fft - f * eta_fft
-    k2 = kx ** 2 + ky ** 2
-    ageo_fft = f * rot_fft / c2 + k2 * eta_fft
-    return q_fft, div_fft, ageo_fft
diff --git a/fluidsim/solvers/sw1l/solver.py b/fluidsim/solvers/sw1l/solver.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3NvbHZlci5weQ==..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/sw1l/solver.py
+++ b/fluidsim/solvers/sw1l/solver.py
@@ -16,7 +16,9 @@
 
 
 """
-from __future__ import division
+
+from fluidpythran import cachedjit, Array
+from fluiddyn.util import mpi
 
 from fluidsim.base.setofvariables import SetOfVariables
 from fluidsim.base.solvers.pseudo_spect import (
@@ -24,5 +26,6 @@
     InfoSolverPseudoSpectral,
 )
 
-from fluiddyn.util import mpi
+A = Array[float, "2d"]
+
 
@@ -28,5 +31,16 @@
 
-from .util_pythran import compute_Frot
+@cachedjit
+def compute_Frot(rot: A, ux: A, uy: A, f: "float or int"):
+    """Compute cross-product of absolute potential vorticity with velocity."""
+    if f != 0:
+        rot_abs = rot + f
+    else:
+        rot_abs = rot
+
+    F1x = rot_abs * uy
+    F1y = -rot_abs * ux
+
+    return F1x, F1y
 
 
 class InfoSolverSW1L(InfoSolverPseudoSpectral):
diff --git a/fluidsim/solvers/sw1l/util_oper_pythran.py b/fluidsim/solvers/sw1l/util_oper_pythran.py
deleted file mode 100644
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3V0aWxfb3Blcl9weXRocmFuLnB5..0000000000000000000000000000000000000000
--- a/fluidsim/solvers/sw1l/util_oper_pythran.py
+++ /dev/null
@@ -1,78 +0,0 @@
-import numpy as np
-
-
-# pythran export qapamfft_from_uxuyetafft(
-#     complex128[][], complex128[][], complex128[][],
-#     int, int, float64[][], float64[][], float64[][],
-#     complex128[][], float, float, int)
-
-
-def qapamfft_from_uxuyetafft(
-    ux_fft, uy_fft, eta_fft, n0, n1, KX, KY, K2, Kappa_over_ic, f, c2, rank
-):
-    """Calculate normal modes from primitive variables."""
-    freq_Corio = f
-    f_over_c2 = freq_Corio / c2
-
-    q_fft = np.empty([n0, n1], dtype=np.complex128)
-    ap_fft = np.empty([n0, n1], dtype=np.complex128)
-    am_fft = np.empty([n0, n1], dtype=np.complex128)
-
-    if freq_Corio != 0:
-        for i0 in range(n0):
-            for i1 in range(n1):
-                if i0 == 0 and i1 == 0 and rank == 0:
-                    q_fft[i0, i1] = 0
-                    ap_fft[i0, i1] = ux_fft[0, 0] + 1.0j * uy_fft[0, 0]
-                    am_fft[i0, i1] = ux_fft[0, 0] - 1.0j * uy_fft[0, 0]
-                else:
-
-                    rot_fft = 1j * (
-                        KX[i0, i1] * uy_fft[i0, i1] - KY[i0, i1] * ux_fft[i0, i1]
-                    )
-
-                    q_fft[i0, i1] = rot_fft - freq_Corio * eta_fft[i0, i1]
-
-                    a_over2_fft = 0.5 * (
-                        K2[i0, i1] * eta_fft[i0, i1] + f_over_c2 * rot_fft
-                    )
-
-                    Deltaa_over2_fft = (
-                        0.5j
-                        * Kappa_over_ic[i0, i1]
-                        * (
-                            KX[i0, i1] * ux_fft[i0, i1]
-                            + KY[i0, i1] * uy_fft[i0, i1]
-                        )
-                    )
-
-                    ap_fft[i0, i1] = a_over2_fft + Deltaa_over2_fft
-                    am_fft[i0, i1] = a_over2_fft - Deltaa_over2_fft
-
-    else:  # (freq_Corio == 0.)
-        for i0 in range(n0):
-            for i1 in range(n1):
-                if i0 == 0 and i1 == 0 and rank == 0:
-                    q_fft[i0, i1] = 0
-                    ap_fft[i0, i1] = ux_fft[0, 0] + 1.0j * uy_fft[0, 0]
-                    am_fft[i0, i1] = ux_fft[0, 0] - 1.0j * uy_fft[0, 0]
-                else:
-                    q_fft[i0, i1] = 1j * (
-                        KX[i0, i1] * uy_fft[i0, i1] - KY[i0, i1] * ux_fft[i0, i1]
-                    )
-
-                    a_over2_fft = 0.5 * K2[i0, i1] * eta_fft[i0, i1]
-
-                    Deltaa_over2_fft = (
-                        0.5j
-                        * Kappa_over_ic[i0, i1]
-                        * (
-                            KX[i0, i1] * ux_fft[i0, i1]
-                            + KY[i0, i1] * uy_fft[i0, i1]
-                        )
-                    )
-
-                    ap_fft[i0, i1] = a_over2_fft + Deltaa_over2_fft
-                    am_fft[i0, i1] = a_over2_fft - Deltaa_over2_fft
-
-    return q_fft, ap_fft, am_fft
diff --git a/fluidsim/solvers/sw1l/util_pythran.py b/fluidsim/solvers/sw1l/util_pythran.py
deleted file mode 100644
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3V0aWxfcHl0aHJhbi5weQ==..0000000000000000000000000000000000000000
--- a/fluidsim/solvers/sw1l/util_pythran.py
+++ /dev/null
@@ -1,16 +0,0 @@
-# pythran export compute_Frot(
-#     float64[][], float64[][], float64[][],
-#     float or int)
-
-
-def compute_Frot(rot, ux, uy, f):
-    """Compute cross-product of absolute potential vorticity with velocity."""
-    if f != 0:
-        rot_abs = rot + f
-    else:
-        rot_abs = rot
-
-    F1x = rot_abs * uy
-    F1y = -rot_abs * ux
-
-    return F1x, F1y
diff --git a/fluidsim/util/util.py b/fluidsim/util/util.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_Zmx1aWRzaW0vdXRpbC91dGlsLnB5..204e59ec7e8081cc658c07dd10441e7771db33c4_Zmx1aWRzaW0vdXRpbC91dGlsLnB5 100644
--- a/fluidsim/util/util.py
+++ b/fluidsim/util/util.py
@@ -196,7 +196,7 @@
     for ii, name in enumerate(name_files):
         tmp = ".".join(name[ind_start_time:].split(".")[:2])
         if "_" in tmp:
-            tmp = tmp[:tmp.index("_")]
+            tmp = tmp[: tmp.index("_")]
         times[ii] = float(tmp)
     if t_approx is None:
         t_approx = times.max()
diff --git a/pyproject.toml b/pyproject.toml
index 7c7294a126136b4f385c1e688d13398dbac7ec44_cHlwcm9qZWN0LnRvbWw=..204e59ec7e8081cc658c07dd10441e7771db33c4_cHlwcm9qZWN0LnRvbWw= 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,2 +1,2 @@
 [build-system]
-requires = ["setuptools", "wheel", "numpy", "fluidpythran>=0.1.0.post0"]
\ No newline at end of file
+requires = ["setuptools", "wheel", "numpy", "fluidpythran>=0.1.2"]
\ No newline at end of file
diff --git a/setup.cfg b/setup.cfg
index 7c7294a126136b4f385c1e688d13398dbac7ec44_c2V0dXAuY2Zn..204e59ec7e8081cc658c07dd10441e7771db33c4_c2V0dXAuY2Zn 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -21,5 +21,5 @@
 
 [options]
 setup_requires =
-    fluidpythran>=0.1.1
+    fluidpythran>=0.1.2
     numpy
diff --git a/setup.py b/setup.py
index 7c7294a126136b4f385c1e688d13398dbac7ec44_c2V0dXAucHk=..204e59ec7e8081cc658c07dd10441e7771db33c4_c2V0dXAucHk= 100644
--- a/setup.py
+++ b/setup.py
@@ -53,8 +53,14 @@
 
 from fluidpythran.dist import make_pythran_files
 
-paths = ["fluidsim/base/time_stepping/pseudo_spect.py"]
+paths = [
+    "fluidsim/base/time_stepping/pseudo_spect.py",
+    "fluidsim/base/output/increments.py",
+    "fluidsim/operators/operators2d.py",
+    "fluidsim/operators/operators3d.py",
+    "fluidsim/solvers/ns2d/solver.py",
+]
 make_pythran_files(
     [here / path for path in paths],
     mocked_modules=(
         "psutil",
@@ -57,8 +63,9 @@
 make_pythran_files(
     [here / path for path in paths],
     mocked_modules=(
         "psutil",
+        "h5py",
         "fluiddyn",
         "fluiddyn.io",
         "fluiddyn.util",
         "fluiddyn.util.paramcontainer",
@@ -61,8 +68,8 @@
         "fluiddyn",
         "fluiddyn.io",
         "fluiddyn.util",
         "fluiddyn.util.paramcontainer",
-        "h5py",
+        "fluiddyn.output",
     ),
 )
 
@@ -133,7 +140,7 @@
     "future >= 0.16",
     "h5py",
     "h5netcdf",
-    "fluidpythran >= 0.1.1",
+    "fluidpythran>=0.1.2",
     "setuptools_scm",
     "xarray",
 ]
@@ -148,9 +155,7 @@
 
 
 def make_pythran_extensions(modules):
-    exclude_pythran = tuple(
-        key for key, value in config["exclude_pythran"].items() if value
-    )
+    exclude_pythran = tuple()
     if len(exclude_pythran) > 0:
         logger.info(
             "Pythran files in the packages "
@@ -197,7 +202,7 @@
         for name in files:
             if (
                 name.endswith("_pythran.py")
-                or path_dir.name == "_pythran"
+                or path_dir.name == "__pythran__"
                 and name.endswith(".py")
             ):
                 path = os.path.join(root, name)
diff --git a/site.cfg.default b/site.cfg.default
index 7c7294a126136b4f385c1e688d13398dbac7ec44_c2l0ZS5jZmcuZGVmYXVsdA==..204e59ec7e8081cc658c07dd10441e7771db33c4_c2l0ZS5jZmcuZGVmYXVsdA== 100644
--- a/site.cfg.default
+++ b/site.cfg.default
@@ -1,11 +1,3 @@
-[exclude_pythran]
-fluidsim.solvers.plate2d = True
-fluidsim.solvers.plate2d.output = True
-fluidsim.solvers.sw1l = True
-fluidsim.solvers.sw1l.output = True
-fluidsim.solvers.ns2d.strat = True
-fluidsim.solvers.ns2d.bouss = True
-
 [environ]
 
 
diff --git a/tox.ini b/tox.ini
index 7c7294a126136b4f385c1e688d13398dbac7ec44_dG94LmluaQ==..204e59ec7e8081cc658c07dd10441e7771db33c4_dG94LmluaQ== 100644
--- a/tox.ini
+++ b/tox.ini
@@ -22,7 +22,7 @@
 deps =
     pip==18.0
     coverage
-    fluidpythran>=0.1.1
+    fluidpythran>=0.1.2
     mako
     numpy
     matplotlib