From f07ca9ac521c72188fd9bf65378c73608d01a02b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Fri, 4 Nov 2016 18:01:50 +0100
Subject: [PATCH] Fortran interface/fftw: - set precision kind as a input
 parameter; fix f2py bug - update fftw interface (fortran) and diffusion op

---
 CMakeLists.txt                                |  27 +++-
 cmake/fortran_utils.cmake                     |   6 +-
 docs/sphinx/users_guide/diffusion.rst         |  43 ++++++
 docs/sphinx/users_guide/operators.rst         |   1 +
 hysop/common_f/template.f95                   |  15 +-
 hysop/common_f/template.pyf                   |   3 +-
 hysop/f2py/fftw2py.pyf                        | 130 ------------------
 hysop/fortran/template.f95                    |   6 +-
 hysop/fortran/template.pyf                    |   5 +-
 {src/fftw => hysop/numerics/fftw_f}/fft2d.f90 |  58 ++++----
 {src/fftw => hysop/numerics/fftw_f}/fft3d.f90 | 129 +++++++++--------
 hysop/{f2py => numerics/fftw_f}/fftw2py.f90   | 104 +++++++-------
 hysop/numerics/fftw_f/fftw2py.pyf             | 130 ++++++++++++++++++
 hysop/operator/computational.py               |   4 +-
 hysop/operator/diffusion.py                   |  11 +-
 hysop/operator/discrete/diffusion_fft.py      |  66 +++++----
 hysop/operator/tests/test_diffusion.py        |  22 +--
 hysop/operators.py                            |   2 +
 precision.conf.in                             |  37 +++++
 setup.py.in                                   |  23 ++--
 20 files changed, 468 insertions(+), 354 deletions(-)
 create mode 100644 docs/sphinx/users_guide/diffusion.rst
 delete mode 100644 hysop/f2py/fftw2py.pyf
 rename {src/fftw => hysop/numerics/fftw_f}/fft2d.f90 (93%)
 rename {src/fftw => hysop/numerics/fftw_f}/fft3d.f90 (93%)
 rename hysop/{f2py => numerics/fftw_f}/fftw2py.f90 (77%)
 create mode 100644 hysop/numerics/fftw_f/fftw2py.pyf
 create mode 100644 precision.conf.in

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3589391f0..ccbee338c 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -28,10 +28,10 @@ include(OutOfSourceBuild)
 include(MyTools)
 
 # User defined options
-option(VERBOSE_MODE "enable verbose mode for cmake exec. Default = on" ON)
-option(DOUBLEPREC "set precision for real numbers to double precision when this mode is enable. Default = on." ON)
-option(USE_MPI "compile and link HySoP with mpi when this mode is enable. Default = on." ON)
-option(WITH_TESTS "Enable testing. Default = off" ON)
+option(VERBOSE_MODE "Print cmake-hysop state summary. Default = on" ON)
+option(DOUBLEPREC "switch between single and double precision for real numbers. Default = on, ie double precision." ON)
+option(USE_MPI "compile and link HySoP with mpi when this mode is enable. Default = on. WARNING: off mode is bugged." ON)
+option(WITH_TESTS "Enable testing. Default = OFF" ON)
 option(BUILD_SHARED_LIBS "Enable dynamic library build, default = ON." ON)
 option(USE_CXX "Expand hysop with some new functions from a generated c++ to python interface, wrapped into hysop.cpp2hysop module. Default = ON." ON)
 option(WITH_SCALES "compile/create scales lib and link it with HySoP. Default = ON." ON)
@@ -74,7 +74,9 @@ endif()
 if(NOT USE_MPI)
   message(FATAL_ERROR "No-mpi version of hysop is broken, please enable mpi.")
 endif()
-set(USE_FORTRAN "ON")
+
+# Fortran interface is not used by default.
+set(USE_FORTRAN "OFF")
 if(WITH_FFTW OR WITH_SCALES OR WITH_EXTRAS)
   set(USE_FORTRAN "ON")
 endif()
@@ -273,10 +275,13 @@ include(HySoPInstallSetup)
 # set precision for real numbers
 # used to generate precision.f95 file/module and constant.py
 if(DOUBLEPREC)
+  # -- Variables required in fortran files --
   set(WORKING_PRECISION dp)
   set(MPI_PRECISION MPI_DOUBLE_PRECISION)
+  # -- Variables required in python files --
   set(PYTHON_PREC np.float64)
   set(MPI_PYTHON_PREC MPI.DOUBLE)
+  # -- f2py stuff --
   set(f2pymap_for_real double)
 else()
   set(WORKING_PRECISION sp)
@@ -286,7 +291,7 @@ else()
   set(f2pymap_for_real float)
 endif()
 
-# set data layout ('fortran' order or 'C' order)
+# set data layout ('fortran' order or 'C' order) in python (numpy).
 if(FORTRAN_LAYOUT)
   set(DATA_LAYOUT 'F')
   set(MPI_DATA_LAYOUT MPI.ORDER_F)
@@ -325,6 +330,12 @@ if(USE_FORTRAN)
   # add a line 'include fftw2py.py'
   include(fortran_utils)
   write_main_pyf_file(f2hysop)
+  # f2py failed with:
+  # real(kind=wp), dimension(:) :: tab
+  # whit wp = real64 or real32 from iso_fortran_env
+  # a const parameter defined in another module
+  # AND tab is an assumed shaped array, even when using a f2py_f2cmap file.
+  # The solution is to set directly real(kind=8) in pyf file.
 
   # --- fortran libraries setup ---
 
@@ -346,11 +357,13 @@ if(USE_FORTRAN)
 
   set(Fortran_FLAGS ${CMAKE_Fortran_FLAGS})
   append_flags(Fortran_FLAGS ${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}})
-  # --- Generate precision.f95 file ---
+  # --- Generate precision.f95 and precision.pyf files ---
   if(EXISTS ${CMAKE_SOURCE_DIR}/precision.conf.in)
     message(STATUS "Generate precision.f95 file ...")
     configure_file(${CMAKE_SOURCE_DIR}/precision.conf.in
       ${GENERATED_FORTRAN_FILES_DIR}/precision.f95)
+    configure_file(${CMAKE_SOURCE_DIR}/precision.conf.in
+      ${GENERATED_FORTRAN_FILES_DIR}/precision.pyf)
   endif()
 endif()
 
diff --git a/cmake/fortran_utils.cmake b/cmake/fortran_utils.cmake
index 96c2f7ddb..b98468b62 100644
--- a/cmake/fortran_utils.cmake
+++ b/cmake/fortran_utils.cmake
@@ -25,6 +25,9 @@ function(write_main_pyf_file filename)
 ! Note: the context of this file is case sensitive.\n
 python module f2hysop ! in\n
   interface\n")
+ file(APPEND ${_file} 
+      "      ! Example
+      include '@GENERATED_FORTRAN_FILES_DIR@/precision.pyf'\n")
  file(APPEND ${_file} 
       "      ! Example
       include '@CMAKE_SOURCE_DIR@/hysop/common_f/template.pyf'\n")
@@ -34,7 +37,7 @@ python module f2hysop ! in\n
     if(WITH_FFTW)
       file(APPEND ${_file}
       "      ! fftw
-      include '@CMAKE_SOURCE_DIR@/hysop/f2py/fftw2py.pyf'\n")
+      include '@CMAKE_SOURCE_DIR@/hysop/numerics/fftw_f/fftw2py.pyf'\n")
     endif()
     if(WITH_SCALES)
       file(APPEND ${_file}
@@ -53,3 +56,4 @@ message(STATUS "Generate f2hysop.pyf file ...")
 configure_file(${_file} ${GENERATED_FORTRAN_FILES_DIR}/${filename}.pyf @ONLY)
 
 endfunction()
+
diff --git a/docs/sphinx/users_guide/diffusion.rst b/docs/sphinx/users_guide/diffusion.rst
new file mode 100644
index 000000000..2356fbe56
--- /dev/null
+++ b/docs/sphinx/users_guide/diffusion.rst
@@ -0,0 +1,43 @@
+.. _diffusion:
+.. currentmodule:: hysop.operator.diffusion
+
+Diffusion Operator
+------------------
+
+
+The class :class:`Diffusion` is used to model and solve the following problem :
+
+.. math::
+   
+   \frac{\partial \omega}{\partial t} &=& \nu\Delta\omega
+
+for the unknown vector field :math:`\omega` and a given constant viscosity value :math:`\nu`,
+
+with
+
+.. math::
+   \omega = \nabla \times velocity \\
+   \nabla . velocity = 0
+
+
+
+Examples
+^^^^^^^^
+
+
+::
+
+  from hysop.operator.diffusion import Diffusion
+  dom = Box()
+  d3d = Discretization([33, 33, 33])
+  topo = dom.create_topology(d3d)
+  om = Field(name='vorticity', is_vector=True, domain=dom)
+  om.randomize(topo)
+  op = Diffusion(vorticity=om, viscosity=0.3, discretization=topo)
+  op.discretize()
+  op.setup()
+  sim = Simulation(nb_iter=10)
+  sim.initialize()
+  while not sim.is_over:
+      op.apply(sim)
+      sim.advance()
diff --git a/docs/sphinx/users_guide/operators.rst b/docs/sphinx/users_guide/operators.rst
index 861af66d7..6da590e09 100644
--- a/docs/sphinx/users_guide/operators.rst
+++ b/docs/sphinx/users_guide/operators.rst
@@ -110,6 +110,7 @@ Review
    
    absorption_bc
    poisson
+   diffusion
    forces
    penalisation
    adaptive_time_step
diff --git a/hysop/common_f/template.f95 b/hysop/common_f/template.f95
index 1a178b268..cbe57da08 100644
--- a/hysop/common_f/template.f95
+++ b/hysop/common_f/template.f95
@@ -4,8 +4,6 @@
 
 module template_f2py
 
-  use precision
- 
   implicit none
 
   !> a global var
@@ -25,5 +23,18 @@ contains
     output(:,:) = 2 * input(:,:)
     print *, 'aha hah', output(1,1)
   end subroutine check_f2py
+
+  subroutine check_f2py_not(input, output)
+    use precision
+    !> input array
+    real(kind=wp), dimension(:,:), intent(in) :: input
+    !> output array
+    real(kind=wp), dimension(:,:), intent(inout) :: output
+
+    print *, 'template f2py for tab'
+    output(:,:) = 2 * input(:,:)
+    print *, 'aha hah', output(1,1)
+  end subroutine check_f2py_not
+
   
 end module template_f2py
diff --git a/hysop/common_f/template.pyf b/hysop/common_f/template.pyf
index 6245f4750..f637d0be3 100644
--- a/hysop/common_f/template.pyf
+++ b/hysop/common_f/template.pyf
@@ -1,8 +1,7 @@
 !    -*- f90 -*-
 ! Note: the context of this file is case sensitive.
 module template_f2py ! in :template_f2py:template.f95
-  use precision
-  use ISO_FORTRAN_ENV
+
   integer, parameter,optional :: var1=12
   subroutine check_f2py(input,output) ! in :template_f2py:template.f95:template_f2py
     use precision
diff --git a/hysop/f2py/fftw2py.pyf b/hysop/f2py/fftw2py.pyf
deleted file mode 100644
index 860933a2d..000000000
--- a/hysop/f2py/fftw2py.pyf
+++ /dev/null
@@ -1,130 +0,0 @@
-!    -*- f90 -*-
-! Note: the context of this file is case sensitive.
-
-module fftw2py ! in fftw2py.f90
-    use fft3d
-    !use mpi
-    use client_data
-    use hysopparam
-    use fft2d
-    subroutine init_fftw_solver(resolution,lengths,comm,datashape,offset,dim,fftw_type_real) ! in fftw2py.f90:fftw2py
-        integer dimension(dim),intent(in) :: resolution
-        real(kind=pk) dimension(dim),intent(in),depend(dim) :: lengths
-        integer intent(in) :: comm
-        integer(kind=ik) dimension(dim),intent(out),depend(dim) :: datashape
-        integer(kind=ik) dimension(dim),intent(out),depend(dim) :: offset
-        integer, optional,intent(hide),check(len(resolution)>=dim),depend(resolution) :: dim=len(resolution)
-        logical, optional,intent(in) :: fftw_type_real=1
-    end subroutine init_fftw_solver
-    subroutine init_fftw_solver_scalar(resolution,lengths,comm,datashape,offset,dim,fftw_type_real) ! in fftw2py.f90:fftw2py
-        integer dimension(dim),intent(in) :: resolution
-        real(kind=pk) dimension(dim),intent(in),depend(dim) :: lengths
-        integer intent(in) :: comm
-        integer(kind=ik) dimension(dim),intent(out),depend(dim) :: datashape
-        integer(kind=ik) dimension(dim),intent(out),depend(dim) :: offset
-        integer, optional,intent(hide),check(len(resolution)>=dim),depend(resolution) :: dim=len(resolution)
-        logical, optional,intent(in) :: fftw_type_real=1
-    end subroutine init_fftw_solver_scalar
-    subroutine clean_fftw_solver(dim) ! in fftw2py.f90:fftw2py
-        integer intent(in) :: dim
-    end subroutine clean_fftw_solver
-    subroutine solve_poisson_2d(omega,velocity_x,velocity_y,ghosts_vort,ghosts_velo) ! in fftw2py.f90:fftw2py
-        real(kind=pk) dimension(:,:),intent(in) :: omega
-        real(kind=pk) dimension(size(omega,1),size(omega,2)),intent(in,out),depend(omega,omega) :: velocity_x
-        real(kind=pk) dimension(size(omega,1),size(omega,2)),intent(in,out),depend(omega,omega) :: velocity_y
-        integer dimension(2),intent(in) :: ghosts_vort
-        integer dimension(2),intent(in) :: ghosts_velo
-    end subroutine solve_poisson_2d
-    subroutine solve_diffusion_2d(nudt,omega,ghosts_vort) ! in fftw2py.f90:fftw2py
-        real(kind=pk) intent(in) :: nudt
-        real(kind=pk) dimension(:,:),intent(in,out) :: omega
-        integer dimension(2),intent(in) :: ghosts_vort
-    end subroutine solve_diffusion_2d
-    subroutine solve_poisson_3d(omega_x,omega_y,omega_z,velocity_x,velocity_y,velocity_z,ghosts_vort,ghosts_velo) ! in fftw2py.f90:fftw2py
-        real(kind=pk) dimension(:,:,:),intent(in) :: omega_x
-        real(kind=pk) dimension(:,:,:),intent(in) :: omega_y
-        real(kind=pk) dimension(:,:,:),intent(in) :: omega_z
-        real(kind=pk) dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(in,out),depend(omega_x,omega_y,omega_z) :: velocity_x
-        real(kind=pk) dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(in,out),depend(omega_x,omega_y,omega_z) :: velocity_y
-        real(kind=pk) dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(in,out),depend(omega_x,omega_y,omega_z) :: velocity_z
-        integer dimension(3),intent(in) :: ghosts_vort
-        integer dimension(3),intent(in) :: ghosts_velo
-    end subroutine solve_poisson_3d
-    subroutine solve_poisson_2d_c(omega,velocity_x,velocity_y) ! in fftw2py.f90:fftw2py
-        complex(kind=pk) dimension(:,:),intent(in) :: omega
-        complex(kind=pk) dimension(size(omega,1),size(omega,2)),intent(in,out),depend(omega,omega) :: velocity_x
-        complex(kind=pk) dimension(size(omega,1),size(omega,2)),intent(in,out),depend(omega,omega) :: velocity_y
-    end subroutine solve_poisson_2d_c
-    subroutine solve_poisson_3d_c(omega_x,omega_y,omega_z,velocity_x,velocity_y,velocity_z) ! in fftw2py.f90:fftw2py
-        complex(kind=pk) dimension(:,:,:),intent(in) :: omega_x
-        complex(kind=pk) dimension(:,:,:),intent(in) :: omega_y
-        complex(kind=pk) dimension(:,:,:),intent(in) :: omega_z
-        complex(kind=pk) dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(in,out),depend(omega_x,omega_y,omega_z) :: velocity_x
-        complex(kind=pk) dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(in,out),depend(omega_x,omega_y,omega_z) :: velocity_y
-        complex(kind=pk) dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(in,out),depend(omega_x,omega_y,omega_z) :: velocity_z
-    end subroutine solve_poisson_3d_c
-    subroutine solve_curl_diffusion_3d(nudt,velocity_x,velocity_y,velocity_z,omega_x,omega_y,omega_z,ghosts_velo,ghosts_vort) ! in fftw2py.f90:fftw2py
-        real(kind=pk) intent(in) :: nudt
-        real(kind=pk) dimension(:,:,:),intent(in) :: velocity_x
-        real(kind=pk) dimension(:,:,:),intent(in) :: velocity_y
-        real(kind=pk) dimension(:,:,:),intent(in) :: velocity_z
-        real(kind=pk) dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(in,out),depend(velocity_x,velocity_y,velocity_z) :: omega_x
-        real(kind=pk) dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(in,out),depend(velocity_x,velocity_y,velocity_z) :: omega_y
-        real(kind=pk) dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(in,out),depend(velocity_x,velocity_y,velocity_z) :: omega_z
-        integer dimension(3),intent(in) :: ghosts_velo
-        integer dimension(3),intent(in) :: ghosts_vort
-    end subroutine solve_curl_diffusion_3d
-    subroutine solve_diffusion_3d(nudt,omega_x,omega_y,omega_z,ghosts) ! in fftw2py.f90:fftw2py
-        real(kind=pk) intent(in) :: nudt
-        real(kind=pk) dimension(:,:,:),intent(in,out) :: omega_x
-        real(kind=pk) dimension(:,:,:),intent(in,out) :: omega_y
-        real(kind=pk) dimension(:,:,:),intent(in,out) :: omega_z
-        integer dimension(3),intent(in) :: ghosts
-    end subroutine solve_diffusion_3d
-    subroutine projection_om_3d(omega_x,omega_y,omega_z,ghosts) ! in fftw2py.f90:fftw2py
-        real(kind=pk) dimension(:,:,:),intent(in,out) :: omega_x
-        real(kind=pk) dimension(:,:,:),intent(in,out) :: omega_y
-        real(kind=pk) dimension(:,:,:),intent(in,out) :: omega_z
-        integer dimension(3),intent(in) :: ghosts
-    end subroutine projection_om_3d
-    subroutine multires_om_3d(dxf,dyf,dzf,omega_x,omega_y,omega_z,ghosts) ! in fftw2py.f90:fftw2py
-        real(kind=pk) intent(in) :: dxf
-        real(kind=pk) intent(in) :: dyf
-        real(kind=pk) intent(in) :: dzf
-        real(kind=pk) dimension(:,:,:),intent(in,out) :: omega_x
-        real(kind=pk) dimension(:,:,:),intent(in,out) :: omega_y
-        real(kind=pk) dimension(:,:,:),intent(in,out) :: omega_z
-        integer dimension(3),intent(in) :: ghosts
-    end subroutine multires_om_3d
-    subroutine pressure_3d(pressure,ghosts) ! in fftw2py.f90:fftw2py
-        real(kind=pk) dimension(:,:,:),intent(in,out) :: pressure
-        integer dimension(3),intent(in) :: ghosts
-    end subroutine pressure_3d
-    subroutine solve_curl_3d(velocity_x,velocity_y,velocity_z,omega_x,omega_y,omega_z,ghosts_velo,ghosts_vort) ! in fftw2py.f90:fftw2py
-        real(kind=pk) dimension(:,:,:),intent(in) :: velocity_x
-        real(kind=pk) dimension(:,:,:),intent(in) :: velocity_y
-        real(kind=pk) dimension(:,:,:),intent(in) :: velocity_z
-        real(kind=pk) dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(in,out),depend(velocity_x,velocity_y,velocity_z) :: omega_x
-        real(kind=pk) dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(in,out),depend(velocity_x,velocity_y,velocity_z) :: omega_y
-        real(kind=pk) dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(in,out),depend(velocity_x,velocity_y,velocity_z) :: omega_z
-        integer dimension(3),intent(in) :: ghosts_velo
-        integer dimension(3),intent(in) :: ghosts_vort
-    end subroutine solve_curl_3d
-    subroutine solve_curl_2d(velocity_x,velocity_y,omega_z,ghosts_velo,ghosts_vort) ! in fftw2py.f90:fftw2py
-        real(kind=pk) dimension(:,:),intent(in) :: velocity_x
-        real(kind=pk) dimension(:,:),intent(in) :: velocity_y
-        real(kind=pk) dimension(size(velocity_x,1),size(velocity_x,2)),intent(in,out),depend(velocity_x,velocity_x) :: omega_z
-        integer dimension(2),intent(in) :: ghosts_velo
-        integer dimension(2),intent(in) :: ghosts_vort
-    end subroutine solve_curl_2d
-    subroutine spectrum_3d(field,spectrum,wavelengths,ghosts,length) ! in fftw2py.f90:fftw2py
-        real(kind=pk) dimension(:,:,:),intent(in) :: field
-        real(kind=pk) dimension(:),intent(in,out) :: spectrum
-        real(kind=pk) dimension(:),intent(in,out) :: wavelengths
-        integer dimension(3),intent(in) :: ghosts
-        real(kind=pk) intent(in) :: length
-    end subroutine spectrum_3d
-end module fftw2py
-
-! This file was auto-generated with f2py (version:2).
-! See http://cens.ioc.ee/projects/f2py2e/
diff --git a/hysop/fortran/template.f95 b/hysop/fortran/template.f95
index 77685fb29..1e5eadf57 100644
--- a/hysop/fortran/template.f95
+++ b/hysop/fortran/template.f95
@@ -13,11 +13,11 @@ contains
 
   !> do something ...
   subroutine check_f2py(input, output)
-    
+    use precision
     !> input array
-    real(kind=8), dimension(:,:), intent(in) :: input
+    real(kind=wp), dimension(:,:), intent(in) :: input
     !> output array
-    real(kind=8), dimension(:,:), intent(inout) :: output
+    real(kind=wp), dimension(:,:), intent(inout) :: output
 
     print *, 'template f2py for tab'
     output(:,:) = 2 * input(:,:)
diff --git a/hysop/fortran/template.pyf b/hysop/fortran/template.pyf
index 1564dabe0..301c7fa72 100644
--- a/hysop/fortran/template.pyf
+++ b/hysop/fortran/template.pyf
@@ -4,8 +4,9 @@
 module template_f2py ! in template.f95
     integer, parameter,optional :: var1=12
     subroutine check_f2py(input,output) ! in template.f95:template_f2py
-        real(kind=8), dimension(:, :), intent(in) :: input
-        real(kind=8), dimension(size(input,1), size(input,2)), intent(in,out), depend(input,input) :: output
+      use precision
+        real(kind=wp), dimension(:, :), intent(in) :: input
+        real(kind=wp), dimension(size(input,1), size(input,2)), intent(in,out), depend(input,input) :: output
     end subroutine check_f2py
 end module template_f2py
 
diff --git a/src/fftw/fft2d.f90 b/hysop/numerics/fftw_f/fft2d.f90
similarity index 93%
rename from src/fftw/fft2d.f90
rename to hysop/numerics/fftw_f/fft2d.f90
index ab7277265..e3703c922 100755
--- a/src/fftw/fft2d.f90
+++ b/hysop/numerics/fftw_f/fft2d.f90
@@ -19,8 +19,8 @@
 !! dedicated to tests and validation.
 module fft2d
 
-  use client_data
-
+  use, intrinsic :: iso_c_binding
+  use precision
   implicit none
   include 'fftw3-mpi.f03'
 
@@ -79,8 +79,8 @@ contains
   subroutine init_c2c_2d(resolution,lengths)
 
     !! global domain resolution
-    integer, dimension(2), intent(in) :: resolution
-    real(mk),dimension(2), intent(in) :: lengths
+    integer(kind=ip), dimension(2), intent(in) :: resolution
+    real(wp),dimension(2), intent(in) :: lengths
 
     !! Size of the local memory required for fftw (cbuffer)
     integer(C_INTPTR_T) :: alloc_local
@@ -130,7 +130,7 @@ contains
 
     call computeKxC(lengths(c_X))
     call computeKy(lengths(c_Y))
-    normFFT =  1./(fft_resolution(c_X)*fft_resolution(c_Y))
+    normFFT =  ONE/(fft_resolution(c_X)*fft_resolution(c_Y))
     !! call fft2d_diagnostics(alloc_local)
 
     is2DUpToDate = .true.
@@ -142,8 +142,8 @@ contains
 
   !> Execute fftw forward transform, according to the pre-defined plans.
   subroutine c2c_2d(inputData,velocity_x,velocity_y)
-    complex(mk),dimension(:,:) :: velocity_x,velocity_y
-    complex(mk),dimension(:,:),intent(in) :: inputData
+    complex(wp),dimension(:,:) :: velocity_x,velocity_y
+    complex(wp),dimension(:,:),intent(in) :: inputData
 
     integer(C_INTPTR_T) :: i, j
 
@@ -190,8 +190,8 @@ contains
   !! @param[in] resolution global domain resolution
   subroutine init_r2c_2d(resolution,lengths)
 
-    integer, dimension(2), intent(in) :: resolution
-    real(mk),dimension(2), intent(in) :: lengths
+    integer(kind=ip), dimension(2), intent(in) :: resolution
+    real(wp),dimension(2), intent(in) :: lengths
     !! Size of the local memory required for fftw (cbuffer)
     integer(C_INTPTR_T) :: alloc_local,halfLength
 
@@ -233,7 +233,7 @@ contains
 
     call computeKx(lengths(c_X))
     call computeKy(lengths(c_Y))
-    normFFT = 1./(fft_resolution(c_X)*fft_resolution(c_Y))
+    normFFT = ONE/(fft_resolution(c_X)*fft_resolution(c_Y))
     !! call fft2d_diagnostics(alloc_local)
 !!$
 !!$    write(*,'(a,i5,a,16f10.4)') 'kx[',rank,'] ', kx
@@ -248,8 +248,8 @@ contains
   !! @param input data
   subroutine r2c_scalar_2d(inputData, ghosts)
 
-    real(mk),dimension(:,:), intent(in) :: inputData
-    integer, dimension(2), intent(in) :: ghosts
+    real(wp),dimension(:,:), intent(in) :: inputData
+    integer(kind=ip), dimension(2), intent(in) :: ghosts
 
     integer(C_INTPTR_T) :: i, j, ig, jg
     ! init
@@ -278,8 +278,8 @@ contains
   !! @param[in] input data
   subroutine r2c_2d(input_x, input_y, ghosts)
 
-    real(mk),dimension(:,:), intent(in) :: input_x, input_y
-    integer, dimension(2), intent(in) :: ghosts
+    real(wp),dimension(:,:), intent(in) :: input_x, input_y
+    integer(kind=ip), dimension(2), intent(in) :: ghosts
 
     integer(C_INTPTR_T) :: i, j, ig, jg
     ! init
@@ -300,8 +300,8 @@ contains
 
   !> Backward transform
   subroutine c2r_2d(velocity_x,velocity_y, ghosts)
-    real(mk),dimension(:,:),intent(inout) :: velocity_x,velocity_y
-    integer, dimension(2), intent(in) :: ghosts
+    real(wp),dimension(:,:),intent(inout) :: velocity_x,velocity_y
+    integer(kind=ip), dimension(2), intent(in) :: ghosts
     integer(C_INTPTR_T) :: i, j, ig, jg
 
     call fftw_mpi_execute_dft_c2r(plan_backward1,dataout1,rdatain1)
@@ -330,8 +330,8 @@ contains
 
   !> Backward transform for scalar field
   subroutine c2r_scalar_2d(omega, ghosts)
-    real(mk),dimension(:,:),intent(inout) :: omega
-    integer, dimension(2), intent(in) :: ghosts
+    real(wp),dimension(:,:),intent(inout) :: omega
+    integer(kind=ip), dimension(2), intent(in) :: ghosts
     integer(C_INTPTR_T) :: i, j, ig, jg
 
     call fftw_mpi_execute_dft_c2r(plan_backward1,dataout1,rdatain1)
@@ -365,7 +365,7 @@ contains
   !> @param lengths size of the domain
   subroutine computeKx(length)
 
-    real(mk),intent(in) :: length
+    real(wp),intent(in) :: length
 
     !! Local loops indices
     integer(C_INTPTR_T) :: i
@@ -383,7 +383,7 @@ contains
   !> @param lengths size of the domain
   subroutine computeKxC(length)
 
-    real(mk),intent(in) :: length
+    real(wp),intent(in) :: length
 
     !! Local loops indices
     integer(C_INTPTR_T) :: i
@@ -420,7 +420,7 @@ contains
   !> Computation of frequencies coeff, over non-distributed direction(s)
   !> @param lengths size of the domain
   subroutine computeKy(length)
-    real(mk), intent(in) :: length
+    real(wp), intent(in) :: length
 
     !! Local loops indices
     integer(C_INTPTR_T) :: i
@@ -491,7 +491,7 @@ contains
 
     do i = 1,local_resolution(c_X)
        do j = 1, fft_resolution(c_Y)
-          coeff = 1./(1. + nudt * (kx(i)**2+ky(j)**2))
+          coeff = ONE/(ONE + nudt * (kx(i)**2+ky(j)**2))
           dataout1(j,i) = coeff*dataout1(j,i)
        end do
     end do
@@ -533,13 +533,13 @@ contains
   subroutine fft2d_diagnostics(nbelem)
     integer(C_INTPTR_T), intent(in) :: nbelem
     complex(C_DOUBLE_COMPLEX) :: memoryAllocated
-    memoryAllocated = real(nbelem*sizeof(memoryAllocated),mk)*1e-6
+    memoryAllocated = real(nbelem*sizeof(memoryAllocated),wp)*1e-6
     write(*,'(a,i5,a,i12,f10.2)') '[',rank,'] size of each buffer (elements / memory in MB):', &
          nbelem, memoryAllocated
     write(*,'(a,i5,a,2i12)') '[',rank,'] size of kx,y,z vectors (number of elements):', &
          size(kx), size(ky)
     write(*,'(a,i5,a,4i5)') '[',rank,'] local resolution and offset :', local_resolution, local_offset
-    memoryAllocated = 2*memoryAllocated + real(sizeof(kx) + sizeof(ky), mk)*1e-6
+    memoryAllocated = 2*memoryAllocated + real(sizeof(kx) + sizeof(ky), wp)*1e-6
     write(*,'(a,i5,a,f10.2)') '[',rank,'] Total memory used for fftw buffers (MB):', memoryAllocated
 
   end subroutine fft2d_diagnostics
@@ -548,8 +548,8 @@ contains
   !! @param datashape local shape of the input field for the fftw process
   !! @param offset index of the first component of the local field (in each dir) in the global set of indices
   subroutine getParamatersTopologyFFTW2d(datashape,offset)
-    integer(C_INTPTR_T), intent(out),dimension(2) :: datashape
-    integer(C_INTPTR_T), intent(out),dimension(2) :: offset
+    integer(kind=ip), intent(out),dimension(2) :: datashape
+    integer(kind=ip), intent(out),dimension(2) :: offset
     integer(C_INTPTR_T) :: offsetx = 0
     datashape = (/fft_resolution(c_X), local_resolution(c_Y)/)
     offset = (/ offsetx, local_offset(c_Y)/)
@@ -559,8 +559,8 @@ contains
   !! @param[in] resolution global domain resolution
   subroutine init_r2c_2dBIS(resolution,lengths)
 
-    integer, dimension(2), intent(in) :: resolution
-    real(mk),dimension(2), intent(in) :: lengths
+    integer(kind=ip), dimension(2), intent(in) :: resolution
+    real(wp),dimension(2), intent(in) :: lengths
     !! Size of the local memory required for fftw (cbuffer)
     integer(C_INTPTR_T) :: alloc_local,halfLength,howmany
     integer(C_INTPTR_T), dimension(2) :: n
@@ -607,7 +607,7 @@ contains
 
     call computeKx(lengths(c_X))
     call computeKy(lengths(c_Y))
-    normFFT = 1./(fft_resolution(c_X)*fft_resolution(c_Y))
+    normFFT = ONE/(fft_resolution(c_X)*fft_resolution(c_Y))
     !! call fft2d_diagnostics(alloc_local)
 
   end subroutine init_r2c_2dBIS
diff --git a/src/fftw/fft3d.f90 b/hysop/numerics/fftw_f/fft3d.f90
similarity index 93%
rename from src/fftw/fft3d.f90
rename to hysop/numerics/fftw_f/fft3d.f90
index b9b5d922b..b7afd6c97 100755
--- a/src/fftw/fft3d.f90
+++ b/hysop/numerics/fftw_f/fft3d.f90
@@ -15,8 +15,9 @@
 !! dedicated to tests and validation.
 module fft3d
 
-  use client_data
+  use, intrinsic :: iso_c_binding
   use mpi
+  use precision
   implicit none
   include 'fftw3-mpi.f03'
 
@@ -80,8 +81,8 @@ contains
   !! @param[in] lengths width of each side of the domain
   subroutine init_c2c_3d(resolution,lengths)
 
-    integer, dimension(3), intent(in) :: resolution
-    real(mk),dimension(3), intent(in) :: lengths
+    integer(kind=ip), dimension(3), intent(in) :: resolution
+    real(wp),dimension(3), intent(in) :: lengths
 
     !! Size of the local memory required for fftw (cbuffer)
     integer(C_INTPTR_T) :: alloc_local
@@ -143,7 +144,7 @@ contains
 
     !! call fft3d_diagnostics(alloc_local)
 
-    normFFT = 1./(fft_resolution(c_X)*fft_resolution(c_Y)*fft_resolution(c_Z))
+    normFFT = one/(fft_resolution(c_X)*fft_resolution(c_Y)*fft_resolution(c_Z))
     manycase = .false.
 
     is3DUpToDate = .true.
@@ -158,8 +159,8 @@ contains
   !!  @param[in] velocity_y 3d scalar field, y-component of the output vector field
   !!  @param[in] velocity_z 3d scalar field, z-component of the output vector field
   subroutine c2c_3d(omega_x,omega_y,omega_z,velocity_x,velocity_y,velocity_z)
-    complex(mk),dimension(:,:,:) :: velocity_x,velocity_y,velocity_z
-    complex(mk),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
+    complex(wp),dimension(:,:,:) :: velocity_x,velocity_y,velocity_z
+    complex(wp),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
 
     integer(C_INTPTR_T) :: i,j,k
     ! Copy input data into the fftw buffer
@@ -205,8 +206,8 @@ contains
   !! @param[in] lengths width of each side of the domain
   subroutine init_r2c_3d(resolution,lengths)
 
-    integer, dimension(3), intent(in) :: resolution
-    real(mk),dimension(3), intent(in) :: lengths
+    integer(kind=ip), dimension(3), intent(in) :: resolution
+    real(wp),dimension(3), intent(in) :: lengths
     !! Size of the local memory required for fftw (cbuffer)
     integer(C_INTPTR_T) :: alloc_local,halfLength
 
@@ -263,7 +264,7 @@ contains
     call computeKy(lengths(c_Y))
     call computeKz(lengths(c_Z))
 
-    normFFT = 1./(fft_resolution(c_X)*fft_resolution(c_Y)*fft_resolution(c_Z))
+    normFFT = one/(fft_resolution(c_X)*fft_resolution(c_Y)*fft_resolution(c_Z))
     !! call fft3d_diagnostics(alloc_local)
     manycase = .false.
     is3DUpToDate = .true.
@@ -275,8 +276,8 @@ contains
   !! @param[in] lengths width of each side of the domain
   subroutine init_r2c_scalar_3d(resolution,lengths)
 
-    integer, dimension(3), intent(in) :: resolution
-    real(mk),dimension(3), intent(in) :: lengths
+    integer(kind=ip), dimension(3), intent(in) :: resolution
+    real(wp),dimension(3), intent(in) :: lengths
     !! Size of the local memory required for fftw (cbuffer)
     integer(C_INTPTR_T) :: alloc_local,halfLength
 
@@ -317,7 +318,7 @@ contains
     call computeKy(lengths(c_Y))
     call computeKz(lengths(c_Z))
 
-    normFFT = 1./(fft_resolution(c_X)*fft_resolution(c_Y)*fft_resolution(c_Z))
+    normFFT = one/(fft_resolution(c_X)*fft_resolution(c_Y)*fft_resolution(c_Z))
     !! call fft3d_diagnostics(alloc_local)
     manycase = .false.
     is3DUpToDate = .true.
@@ -331,9 +332,9 @@ contains
   !!  @param[in] ghosts, number of points in the ghost layer of input fields.
   subroutine r2c_3d(omega_x,omega_y,omega_z, ghosts)
 
-    real(mk),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
-    integer, dimension(3), intent(in) :: ghosts
-    !real(mk) :: start
+    real(wp),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
+    integer(kind=ip), dimension(3), intent(in) :: ghosts
+    !real(wp) :: start
     integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
 
     ! ig, jg, kg are used to take into account
@@ -368,9 +369,9 @@ contains
   !!  @param[in,out] velocity_z 3d scalar field, z-component of the output vector field
   !!  @param[in] ghosts, number of points in the ghost layer of in/out velocity field.
   subroutine c2r_3d(velocity_x,velocity_y,velocity_z, ghosts)
-    real(mk),dimension(:,:,:),intent(inout) :: velocity_x,velocity_y,velocity_z
-    integer, dimension(3), intent(in) :: ghosts
-    real(mk) :: start
+    real(wp),dimension(:,:,:),intent(inout) :: velocity_x,velocity_y,velocity_z
+    integer(kind=ip), dimension(3), intent(in) :: ghosts
+    real(wp) :: start
     integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
     start = MPI_WTIME()
     call fftw_mpi_execute_dft_c2r(plan_backward1,dataout1,rdatain1)
@@ -398,9 +399,9 @@ contains
   !!  @param[in] ghosts, number of points in the ghost layer of input field.
   subroutine r2c_scalar_3d(scalar, ghosts)
 
-    real(mk),dimension(:,:,:),intent(in) :: scalar
-    integer, dimension(3), intent(in) :: ghosts
-    real(mk) :: start
+    real(wp),dimension(:,:,:),intent(in) :: scalar
+    integer(kind=ip), dimension(3), intent(in) :: ghosts
+    real(wp) :: start
     integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
 
     ! ig, jg, kg are used to take into account
@@ -429,9 +430,9 @@ contains
   !!  @param[in,out] scalar 3d scalar field
   !!  @param[in] ghosts, number of points in the ghost layer of in/out scalar field.
   subroutine c2r_scalar_3d(scalar, ghosts)
-    real(mk),dimension(:,:,:),intent(inout) :: scalar
-    integer, dimension(3), intent(in) :: ghosts
-    real(mk) :: start
+    real(wp),dimension(:,:,:),intent(inout) :: scalar
+    integer(kind=ip), dimension(3), intent(in) :: ghosts
+    real(wp) :: start
     integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
     start = MPI_WTIME()
     call fftw_mpi_execute_dft_c2r(plan_backward1,dataout1,rdatain1)
@@ -459,8 +460,8 @@ contains
   !! @param[in] lengths width of each side of the domain
   subroutine init_r2c_3d_many(resolution,lengths)
 
-    integer, dimension(3), intent(in) :: resolution
-    real(mk),dimension(3), intent(in) :: lengths
+    integer(kind=ip), dimension(3), intent(in) :: resolution
+    real(wp),dimension(3), intent(in) :: lengths
     !! Size of the local memory required for fftw (cbuffer)
     integer(C_INTPTR_T) :: alloc_local,halfLength,howmany, blocksize
     integer(C_INTPTR_T),dimension(3) :: n
@@ -502,7 +503,7 @@ contains
     call computeKy(lengths(c_Y))
     call computeKz(lengths(c_Z))
 
-    normFFT = 1./(fft_resolution(c_X)*fft_resolution(c_Y)*fft_resolution(c_Z))
+    normFFT = one/(fft_resolution(c_X)*fft_resolution(c_Y)*fft_resolution(c_Z))
     !! call fft3d_diagnostics(alloc_local,1)
     manycase = .true.
 
@@ -517,8 +518,8 @@ contains
   !! @param input data
   subroutine r2c_3d_many(omega_x,omega_y,omega_z)
 
-    real(mk),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
-    real(mk) :: start
+    real(wp),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
+    real(wp) :: start
     integer(C_INTPTR_T) :: i,j,k
 
     ! init
@@ -544,8 +545,8 @@ contains
   !!  @param[in,out] velocity_y 3d scalar field, y-component of the output vector field
   !!  @param[in,out] velocity_z 3d scalar field, z-component of the output vector field
   subroutine c2r_3d_many(velocity_x,velocity_y,velocity_z)
-    real(mk),dimension(:,:,:),intent(inout) :: velocity_x,velocity_y,velocity_z
-    real(mk) :: start
+    real(wp),dimension(:,:,:),intent(inout) :: velocity_x,velocity_y,velocity_z
+    real(wp) :: start
     integer(C_INTPTR_T) :: i,j,k
 
     start = MPI_WTIME()
@@ -571,7 +572,7 @@ contains
   !> @param lengths size of the domain
   subroutine computeKx(length)
 
-    real(mk),intent(in) :: length
+    real(wp),intent(in) :: length
 
     !! Local loops indices
     integer(C_INTPTR_T) :: i
@@ -624,7 +625,7 @@ contains
   !> Computation of frequencies coeff, over non-distributed direction(s)
   !> @param length size of the domain
   subroutine computeKz(length)
-    real(mk),intent(in) :: length
+    real(wp),intent(in) :: length
 
     !! Local loops indices
     integer(C_INTPTR_T) :: i
@@ -711,7 +712,7 @@ contains
     do j = 1,local_resolution(c_Y)
        do k = 1, fft_resolution(c_Z)
           do i = 1, local_resolution(c_X)
-             coeff = Icmplx/(1. + nudt * (kx(i)**2+ky(j)**2+kz(k)**2))
+             coeff = Icmplx/(one + nudt * (kx(i)**2+ky(j)**2+kz(k)**2))
              buffer1 = dataout1(i,k,j)
              buffer2 = dataout2(i,k,j)
              dataout1(i,k,j) = coeff*(ky(j)*dataout3(i,k,j)-kz(k)*dataout2(i,k,j))
@@ -736,7 +737,7 @@ contains
     do j = 1,local_resolution(c_Y)
        do k = 1, fft_resolution(c_Z)
           do i = 1, local_resolution(c_X)
-             coeff = 1./(1. + nudt * (kx(i)**2+ky(j)**2+kz(k)**2))
+             coeff = one/(one + nudt * (kx(i)**2+ky(j)**2+kz(k)**2))
              dataout1(i,k,j) = coeff*dataout1(i,k,j)
              dataout2(i,k,j) = coeff*dataout2(i,k,j)
              dataout3(i,k,j) = coeff*dataout3(i,k,j)
@@ -784,7 +785,7 @@ contains
        dataout2(1,1,1) = dataout2(1,1,1)
        dataout3(1,1,1) = dataout3(1,1,1)
     else
-       coeff = 1./(ky(1)**2)
+       coeff = one/(ky(1)**2)
        buffer2 = dataout2(1,1,1)
        dataout1(1,1,1) = dataout1(1,1,1)
        dataout2(1,1,1) = dataout2(1,1,1) - coeff*ky(1)*(ky(1)*buffer2)
@@ -793,7 +794,7 @@ contains
 
     ! !! mind the transpose -> index inversion between y and z
     do i = 2, local_resolution(c_X)
-       coeff = 1./(kx(i)**2+ky(1)**2+kz(1)**2)
+       coeff = one/(kx(i)**2+ky(1)**2+kz(1)**2)
        buffer1 = dataout1(i,1,1)
        buffer2 = dataout2(i,1,1)
        buffer3 = dataout3(i,1,1)
@@ -805,7 +806,7 @@ contains
     ! !! mind the transpose -> index inversion between y and z
     do k = 2, fft_resolution(c_Z)
        do i = 1, local_resolution(c_X)
-          coeff = 1./(kx(i)**2+ky(1)**2+kz(k)**2)
+          coeff = one/(kx(i)**2+ky(1)**2+kz(k)**2)
           buffer1 = dataout1(i,k,1)
           buffer2 = dataout2(i,k,1)
           buffer3 = dataout3(i,k,1)
@@ -819,7 +820,7 @@ contains
     do j = 2,local_resolution(c_Y)
        do k = 1, fft_resolution(c_Z)
           do i = 1, local_resolution(c_X)
-             coeff = 1./(kx(i)**2+ky(j)**2+kz(k)**2)
+             coeff = one/(kx(i)**2+ky(j)**2+kz(k)**2)
              buffer1 = dataout1(i,k,j)
              buffer2 = dataout2(i,k,j)
              buffer3 = dataout3(i,k,j)
@@ -870,20 +871,20 @@ contains
     if(local_offset(c_Y) == 0) then
        dataout1(1,1,1) = 0.0
     else
-       coeff = -1./(ky(1)**2)
+       coeff = -one/(ky(1)**2)
        dataout1(1,1,1) = coeff*dataout1(1,1,1)
     endif
 
     ! !! mind the transpose -> index inversion between y and z
     do i = 2, local_resolution(c_X)
-       coeff = -1./(kx(i)**2+ky(1)**2)
+       coeff = -one/(kx(i)**2+ky(1)**2)
        dataout1(i,1,1) = coeff*dataout1(i,1,1)
     end do
 
     ! !! mind the transpose -> index inversion between y and z
     do k = 2, fft_resolution(c_Z)
        do i = 1, local_resolution(c_X)
-          coeff = -1./(kx(i)**2+ky(1)**2+kz(k)**2)
+          coeff = -one/(kx(i)**2+ky(1)**2+kz(k)**2)
           dataout1(i,k,1) = coeff*dataout1(i,k,1)
        end do
     end do
@@ -892,7 +893,7 @@ contains
     do j = 2,local_resolution(c_Y)
        do k = 1, fft_resolution(c_Z)
           do i = 1, local_resolution(c_X)
-             coeff = -1./(kx(i)**2+ky(j)**2+kz(k)**2)
+             coeff = -one/(kx(i)**2+ky(j)**2+kz(k)**2)
              dataout1(i,k,j) = coeff*dataout1(i,k,j)
           end do
        end do
@@ -970,7 +971,7 @@ contains
     do j = 1,local_resolution(c_Y)
        do k = 1, fft_resolution(c_Z)
           do i = 1, local_resolution(c_X)
-             coeff = Icmplx/(1. + nudt * kx(i)**2+ky(j)**2+kz(k)**2)
+             coeff = Icmplx/(one + nudt * kx(i)**2+ky(j)**2+kz(k)**2)
              buffer1 = dataout_many(1,i,k,j)
              buffer2 = dataout_many(2,i,k,j)
              dataout_many(1,i,k,j) = coeff*(ky(j)*dataout_many(3,i,k,j)-kz(k)*dataout_many(2,i,k,j))
@@ -1014,13 +1015,13 @@ contains
     else
        nbFields = 3
     end if
-    memoryAllocated = real(nbelem*sizeof(memoryAllocated),mk)*1e-6
+    memoryAllocated = real(nbelem*sizeof(memoryAllocated),wp)*1e-6
     write(*,'(a,i5,a,i12,f10.2)') '[',rank,'] size of each buffer (elements / memory in MB):', &
          nbelem, memoryAllocated
     write(*,'(a,i5,a,3i12)') '[',rank,'] size of kx,y,z vectors (number of elements):', &
          size(kx), size(ky),size(kz)
     write(*,'(a,i5,a,6i5)') '[',rank,'] local resolution and offset :', local_resolution, local_offset
-    memoryAllocated = nbFields*memoryAllocated + real(sizeof(kx) + sizeof(ky) + sizeof(kz), mk)*1e-6
+    memoryAllocated = nbFields*memoryAllocated + real(sizeof(kx) + sizeof(ky) + sizeof(kz), wp)*1e-6
     write(*,'(a,i5,a,f10.2)') '[',rank,'] Total memory used for fftw buffers (MB):', memoryAllocated
 
   end subroutine fft3d_diagnostics
@@ -1029,22 +1030,20 @@ contains
   !! @param datashape local shape of the input field for the fftw process
   !! @param offset index of the first component of the local field (in each dir) in the global set of indices
   subroutine getParamatersTopologyFFTW3d(datashape,offset)
-    integer(C_INTPTR_T), intent(out),dimension(3) :: datashape
-    integer(C_INTPTR_T), intent(out),dimension(3) :: offset
-    integer(C_INTPTR_T) :: zero = 0
+    integer(kind=ip), intent(out),dimension(3) :: datashape
+    integer(kind=ip), intent(out),dimension(3) :: offset
+    integer(C_INTPTR_T) :: i_zero = 0
     datashape = (/fft_resolution(c_X), fft_resolution(c_Y), local_resolution(c_Z)/)
-    offset = (/zero, zero, local_offset(c_Z)/)
-
+    offset = (/i_zero, i_zero, local_offset(c_Z)/)
   end subroutine getParamatersTopologyFFTW3d
 
-
   !> forward transform - The result is saved in local buffers
   !!  @param[in] field 3d scalar field, x-component of the input vector field
   !!  @param[in] ghosts, number of points in the ghost layer of input fields.
   subroutine r2c_3d_scal(field, ghosts)
 
-    real(mk),dimension(:,:,:),intent(in) :: field
-    integer, dimension(3), intent(in) :: ghosts
+    real(wp),dimension(:,:,:),intent(in) :: field
+    integer(kind=ip), dimension(3), intent(in) :: ghosts
     !real(8) :: start
     integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
 
@@ -1072,28 +1071,28 @@ contains
   !> Compute spectrum of the given data
   subroutine filter_spectrum_3d(spectrum,wavelengths,length)
 
-    real(mk),dimension(:),intent(inout) :: spectrum
-    real(mk),dimension(:),intent(inout) :: wavelengths
-    real(mk),intent(in) :: length
+    real(wp),dimension(:),intent(inout) :: spectrum
+    real(wp),dimension(:),intent(inout) :: wavelengths
+    real(wp),intent(in) :: length
     integer(C_INTPTR_T) :: i,j,k
-    real(mk) :: c
-    real(mk) :: kk,dk,kc,eps
-    integer(C_INTPTR_T) :: ik
+    real(wp) :: c
+    real(wp) :: kk,dk,kc,eps
+    integer(C_INTPTR_T) :: indk
     spectrum = 0
-    dk = 2.0_mk*pi/length
+    dk = 2.0_wp*pi/length
     kc = pi*fft_resolution(c_X)/length
-    eps=kc/1000000.0_mk
+    eps=kc/1000000.0_wp
 
     !! mind the transpose -> index inversion between y and z
-    c = 1._mk/real(fft_resolution(c_Z)*fft_resolution(c_Y)*fft_resolution(c_X),mk)
+    c = one/real(fft_resolution(c_Z)*fft_resolution(c_Y)*fft_resolution(c_X),wp)
     c = c * c
     do j = 1,local_resolution(c_Y)
        do k = 1, fft_resolution(c_Z)
           do i = 1, local_resolution(c_X)
              kk=sqrt(kx(i)**2+ky(j)**2+kz(k)**2)
              if ((kk.gt.eps).and.(kk.le.kc)) then
-                ik=1+int(kk/dk+0.5_mk)
-                spectrum(ik) = spectrum(ik) + real(dataout1(i,j,k)*conjg(dataout1(i,j,k)), mk) * c
+                indk=1+int(kk/dk+0.5_wp)
+                spectrum(indk) = spectrum(indk) + real(dataout1(i,j,k)*conjg(dataout1(i,j,k)), wp) * c
              end if
           end do
        end do
diff --git a/hysop/f2py/fftw2py.f90 b/hysop/numerics/fftw_f/fftw2py.f90
similarity index 77%
rename from hysop/f2py/fftw2py.f90
rename to hysop/numerics/fftw_f/fftw2py.f90
index 040262cf5..d8cd5b731 100755
--- a/hysop/f2py/fftw2py.f90
+++ b/hysop/numerics/fftw_f/fftw2py.f90
@@ -4,17 +4,17 @@
 !> Interface to mpi-fftw (fortran) utilities
 module fftw2py
 
-  use client_data
-  use hysopparam
+  use precision, only: wp, main_comm, ip
   !> 2d case
-   use fft2d
+  use fft2d
   !> 3d case
   use fft3d
-  use mpi
+  use, intrinsic :: iso_c_binding
+  use mpi, only : mpi_comm_dup, mpi_wtime
   implicit none
 
 contains
-
+  
   !> Initialisation of fftw context : create plans and memory buffers
   !! @param[in] resolution global resolution of the discrete domain
   !! @param[in] lengths width of each side of the domain
@@ -24,10 +24,10 @@ contains
   subroutine init_fftw_solver(resolution,lengths,comm,datashape,offset,dim,fftw_type_real)
 
     integer, intent(in) :: dim
-    integer, dimension(dim),intent(in) :: resolution
-    real(pk),dimension(dim), intent(in) :: lengths
-    integer(ik), dimension(dim), intent(out) :: datashape
-    integer(ik), dimension(dim), intent(out) :: offset
+    integer(kind=ip), dimension(dim),intent(in) :: resolution
+    real(wp),dimension(dim), intent(in) :: lengths
+    integer(kind=ip), dimension(dim), intent(out) :: datashape
+    integer(kind=ip), dimension(dim), intent(out) :: offset
     integer, intent(in)                 :: comm
     logical, optional :: fftw_type_real
     !f2py optional :: dim=len(resolution)
@@ -74,10 +74,10 @@ contains
   subroutine init_fftw_solver_scalar(resolution,lengths,comm,datashape,offset,dim,fftw_type_real)
 
     integer, intent(in) :: dim
-    integer, dimension(dim),intent(in) :: resolution
-    real(pk),dimension(dim), intent(in) :: lengths
-    integer(ik), dimension(dim), intent(out) :: datashape
-    integer(ik), dimension(dim), intent(out) :: offset
+    integer(kind=ip), dimension(dim),intent(in) :: resolution
+    real(wp),dimension(dim), intent(in) :: lengths
+    integer(kind=ip), dimension(dim), intent(out) :: datashape
+    integer(kind=ip), dimension(dim), intent(out) :: offset
     integer, intent(in)                 :: comm
     logical, optional :: fftw_type_real
     !f2py optional :: dim=len(resolution)
@@ -111,9 +111,9 @@ contains
   !! \f[ \nabla (\nabla \times velocity) = - \omega \f]
   !! velocity being a 2D vector field and omega a 2D scalar field.
   subroutine solve_poisson_2d(omega,velocity_x,velocity_y, ghosts_vort, ghosts_velo)
-    real(pk),dimension(:,:),intent(in):: omega
-    real(pk),dimension(size(omega,1),size(omega,2)),intent(out) :: velocity_x,velocity_y
-    integer, dimension(2), intent(in) :: ghosts_vort, ghosts_velo
+    real(wp),dimension(:,:),intent(in):: omega
+    real(wp),dimension(size(omega,1),size(omega,2)),intent(out) :: velocity_x,velocity_y
+    integer(kind=ip), dimension(2), intent(in) :: ghosts_vort, ghosts_velo
     !f2py intent(in,out) :: velocity_x,velocity_y
     
     call r2c_scalar_2d(omega, ghosts_vort)
@@ -129,9 +129,9 @@ contains
   !! \f{eqnarray*} \frac{\partial \omega}{\partial t} &=& \nu \Delta \omega \f}
   !! omega being a 2D scalar field.
   subroutine solve_diffusion_2d(nudt, omega, ghosts_vort)
-    real(pk), intent(in) :: nudt
-    real(pk),dimension(:,:),intent(inout):: omega
-    integer, dimension(2), intent(in) :: ghosts_vort
+    real(wp), intent(in) :: nudt
+    real(wp),dimension(:,:),intent(inout):: omega
+    integer(kind=ip), dimension(2), intent(in) :: ghosts_vort
     !f2py intent(in,out) :: omega
 
     call r2c_scalar_2d(omega, ghosts_vort)
@@ -146,10 +146,10 @@ contains
   !! \f{eqnarray*} \Delta \psi &=& - \omega \\ velocity = \nabla\times\psi \f}
   !! velocity and omega being 3D vector fields.
   subroutine solve_poisson_3d(omega_x,omega_y,omega_z,velocity_x,velocity_y,velocity_z, ghosts_vort, ghosts_velo)
-    real(pk),dimension(:,:,:),intent(in):: omega_x,omega_y,omega_z
-    real(pk),dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(out) :: velocity_x,velocity_y,velocity_z
-    integer, dimension(3), intent(in) :: ghosts_vort, ghosts_velo
-    real(pk) :: start
+    real(wp),dimension(:,:,:),intent(in):: omega_x,omega_y,omega_z
+    real(wp),dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(out) :: velocity_x,velocity_y,velocity_z
+    integer(kind=ip), dimension(3), intent(in) :: ghosts_vort, ghosts_velo
+    real(wp) :: start
     !f2py intent(in,out) :: velocity_x,velocity_y,velocity_z
     start = MPI_WTime()
     call r2c_3d(omega_x,omega_y,omega_z, ghosts_vort)
@@ -165,8 +165,8 @@ contains
   !! \f{eqnarray*} \Delta \psi &=& - \omega \\ velocity = \nabla\times\psi \f}
   !! velocity being a 2D complex vector field and omega a 2D complex scalar field.
   subroutine solve_poisson_2d_c(omega,velocity_x,velocity_y)
-    complex(pk),dimension(:,:),intent(in):: omega
-    complex(pk),dimension(size(omega,1),size(omega,2)),intent(out) :: velocity_x,velocity_y
+    complex(wp),dimension(:,:),intent(in):: omega
+    complex(wp),dimension(size(omega,1),size(omega,2)),intent(out) :: velocity_x,velocity_y
     !f2py intent(in,out) :: velocity_x,velocity_y
 
     call c2c_2d(omega,velocity_x,velocity_y)
@@ -177,8 +177,8 @@ contains
   !!  \f{eqnarray*} \Delta \psi &=& - \omega \\ velocity = \nabla\times\psi \f}
   !! velocity and omega being 3D complex vector fields.
   subroutine solve_poisson_3d_c(omega_x,omega_y,omega_z,velocity_x,velocity_y,velocity_Z)
-    complex(pk),dimension(:,:,:),intent(in):: omega_x,omega_y,omega_z
-    complex(pk),dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(out) :: velocity_x,velocity_y,velocity_z
+    complex(wp),dimension(:,:,:),intent(in):: omega_x,omega_y,omega_z
+    complex(wp),dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(out) :: velocity_x,velocity_y,velocity_z
     !f2py intent(in,out) :: velocity_x,velocity_y,velocity_z
 
     call c2c_3d(omega_x,omega_y,omega_z,velocity_x,velocity_y,velocity_z)
@@ -189,10 +189,10 @@ contains
   !! \f{eqnarray*} \omega &=& \nabla \times v \\ \frac{\partial \omega}{\partial t} &=& \nu \Delta \omega \f}
   !! velocity and omega being 3D vector fields.
   subroutine solve_curl_diffusion_3d(nudt,velocity_x,velocity_y,velocity_z,omega_x,omega_y,omega_z, ghosts_velo, ghosts_vort)
-    real(pk), intent(in) :: nudt
-    real(pk),dimension(:,:,:),intent(in):: velocity_x,velocity_y,velocity_z
-    real(pk),dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(out) :: omega_x,omega_y,omega_z
-    integer, dimension(3), intent(in) :: ghosts_vort, ghosts_velo
+    real(wp), intent(in) :: nudt
+    real(wp),dimension(:,:,:),intent(in):: velocity_x,velocity_y,velocity_z
+    real(wp),dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(out) :: omega_x,omega_y,omega_z
+    integer(kind=ip), dimension(3), intent(in) :: ghosts_vort, ghosts_velo
     !f2py intent(in,out) :: omega_x,omega_y,omega_z
 
     call r2c_3d(velocity_x,velocity_y,velocity_z, ghosts_velo)
@@ -207,9 +207,9 @@ contains
   !! \f{eqnarray*} \frac{\partial \omega}{\partial t} &=& \nu \Delta \omega \f}
   !! omega being 3D vector field.
   subroutine solve_diffusion_3d(nudt,omega_x,omega_y,omega_z, ghosts)
-    real(pk), intent(in) :: nudt
-    real(pk),dimension(:,:,:),intent(inout):: omega_x,omega_y,omega_z
-    integer, dimension(3), intent(in) :: ghosts
+    real(wp), intent(in) :: nudt
+    real(wp),dimension(:,:,:),intent(inout):: omega_x,omega_y,omega_z
+    integer(kind=ip), dimension(3), intent(in) :: ghosts
     !f2py intent(in,out) :: omega_x,omega_y,omega_z
 
     call r2c_3d(omega_x,omega_y,omega_z, ghosts)
@@ -224,8 +224,8 @@ contains
   !! \f{eqnarray*} \omega ' &=& \omega - \nabla\pi \f}
   !! omega being a 3D vector field.
   subroutine projection_om_3d(omega_x,omega_y,omega_z, ghosts)
-    real(pk),dimension(:,:,:),intent(inout):: omega_x,omega_y,omega_z
-   integer, dimension(3), intent(in) :: ghosts
+    real(wp),dimension(:,:,:),intent(inout):: omega_x,omega_y,omega_z
+   integer(kind=ip), dimension(3), intent(in) :: ghosts
     !f2py intent(in,out) :: omega_x,omega_y,omega_z
 
     call r2c_3d(omega_x,omega_y,omega_z, ghosts)
@@ -240,9 +240,9 @@ contains
   !! @param[in] dxf, dyf, dzf: grid filter size = domainLength/(CoarseRes-1)
   !! in the following, omega is the 3D vorticity vector field.
   subroutine multires_om_3d(dxf, dyf, dzf, omega_x,omega_y,omega_z, ghosts)
-    real(pk), intent(in) :: dxf, dyf, dzf
-    real(pk),dimension(:,:,:),intent(inout):: omega_x,omega_y,omega_z
-    integer, dimension(3), intent(in) :: ghosts
+    real(wp), intent(in) :: dxf, dyf, dzf
+    real(wp),dimension(:,:,:),intent(inout):: omega_x,omega_y,omega_z
+    integer(kind=ip), dimension(3), intent(in) :: ghosts
 
     !f2py intent(in,out) :: omega_x,omega_y,omega_z
 
@@ -260,8 +260,8 @@ contains
   !! @param[in, out] pressure
   !! in the following, pressure is used as inout parameter. It must contains the rhs of poisson equation.
   subroutine pressure_3d(pressure, ghosts)
-    integer, dimension(3), intent(in) :: ghosts
-    real(pk),dimension(:,:,:),intent(inout):: pressure
+    integer(kind=ip), dimension(3), intent(in) :: ghosts
+    real(wp),dimension(:,:,:),intent(inout):: pressure
     !f2py intent(in,out) :: pressure
 
     call r2c_scalar_3d(pressure, ghosts)
@@ -276,9 +276,9 @@ contains
   !! \f{eqnarray*} \omega &=& \nabla \times v
   !! velocity and omega being 3D vector fields.
   subroutine solve_curl_3d(velocity_x,velocity_y,velocity_z,omega_x,omega_y,omega_z, ghosts_velo, ghosts_vort)
-    real(pk),dimension(:,:,:),intent(in):: velocity_x,velocity_y,velocity_z
-    real(pk),dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(out) :: omega_x,omega_y,omega_z
-    integer, dimension(3), intent(in) :: ghosts_velo, ghosts_vort
+    real(wp),dimension(:,:,:),intent(in):: velocity_x,velocity_y,velocity_z
+    real(wp),dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(out) :: omega_x,omega_y,omega_z
+    integer(kind=ip), dimension(3), intent(in) :: ghosts_velo, ghosts_vort
     !f2py intent(in,out) :: omega_x,omega_y,omega_z
 
     call r2c_3d(velocity_x,velocity_y,velocity_z, ghosts_velo)
@@ -294,9 +294,9 @@ contains
   !! \f{eqnarray*} \omega &=& \nabla \times v
   !! velocity and omega being 2D vector and scalar fields.
   subroutine solve_curl_2d(velocity_x,velocity_y, omega_z, ghosts_velo, ghosts_vort)
-    real(pk), dimension(:,:), intent(in):: velocity_x,velocity_y
-    real(pk), dimension(size(velocity_x,1), size(velocity_x,2)), intent(out) :: omega_z
-    integer, dimension(2), intent(in) :: ghosts_velo, ghosts_vort
+    real(wp), dimension(:,:), intent(in):: velocity_x,velocity_y
+    real(wp), dimension(size(velocity_x,1), size(velocity_x,2)), intent(out) :: omega_z
+    integer(kind=ip), dimension(2), intent(in) :: ghosts_velo, ghosts_vort
     !f2py intent(in,out) :: omega_z
 
     call r2c_2d(velocity_x,velocity_y, ghosts_velo)
@@ -311,11 +311,11 @@ contains
   !! @param[in] field
   !! @param[out] spectrum
   subroutine spectrum_3d(field, spectrum, wavelengths, ghosts, length)
-    real(pk),dimension(:,:,:),intent(in):: field
-    integer, dimension(3), intent(in) :: ghosts
-    real(pk),dimension(:), intent(inout) :: spectrum
-    real(pk),dimension(:), intent(inout) :: wavelengths
-    real(pk),intent(in) :: length
+    real(wp),dimension(:,:,:),intent(in):: field
+    integer(kind=ip), dimension(3), intent(in) :: ghosts
+    real(wp),dimension(:), intent(inout) :: spectrum
+    real(wp),dimension(:), intent(inout) :: wavelengths
+    real(wp),intent(in) :: length
     !f2py intent(in) :: field
     !f2py intent(inout) :: spectrum
     !f2py intent(inout) :: wavelengths
diff --git a/hysop/numerics/fftw_f/fftw2py.pyf b/hysop/numerics/fftw_f/fftw2py.pyf
new file mode 100644
index 000000000..da6d1d4db
--- /dev/null
+++ b/hysop/numerics/fftw_f/fftw2py.pyf
@@ -0,0 +1,130 @@
+!    -*- f90 -*-
+! Note: the context of this file is case sensitive.
+
+module fftw2py ! in fftw2py.f90
+  use fft2d
+  use fft3d
+  use precision, only: wp, main_comm, ip
+  use mpi
+  subroutine init_fftw_solver(resolution,lengths,comm,datashape,offset,dim,fftw_type_real) ! in fftw2py.f90:fftw2py
+    integer(kind=ip), dimension(dim),intent(in) :: resolution
+    real(kind=wp) dimension(dim),intent(in),depend(dim) :: lengths
+    integer intent(in) :: comm
+    integer(kind=ip) dimension(dim),intent(out),depend(dim) :: datashape
+    integer(kind=ip) dimension(dim),intent(out),depend(dim) :: offset
+    integer, optional,intent(hide), depend(resolution) :: dim=len(resolution)
+    logical, optional,intent(in) :: fftw_type_real=1
+  end subroutine init_fftw_solver
+  subroutine init_fftw_solver_scalar(resolution,lengths,comm,datashape,offset,dim,fftw_type_real) ! in fftw2py.f90:fftw2py
+    integer(kind=ip) dimension(dim),intent(in) :: resolution
+    real(kind=wp) dimension(dim),intent(in),depend(dim) :: lengths
+    integer intent(in) :: comm
+    integer(ip) dimension(dim),intent(out),depend(dim) :: datashape
+    integer(ip) dimension(dim),intent(out),depend(dim) :: offset
+    integer, optional,intent(hide), ,depend(resolution) :: dim=len(resolution)
+    logical, optional,intent(in) :: fftw_type_real=1
+  end subroutine init_fftw_solver_scalar
+  subroutine clean_fftw_solver(dim) ! in fftw2py.f90:fftw2py
+    integer intent(in) :: dim
+  end subroutine clean_fftw_solver
+  subroutine solve_poisson_2d(omega,velocity_x,velocity_y,ghosts_vort,ghosts_velo) ! in fftw2py.f90:fftw2py
+    real(kind=wp) dimension(:,:),intent(in) :: omega
+    real(kind=wp) dimension(size(omega,1),size(omega,2)),intent(in,out),depend(omega,omega) :: velocity_x
+    real(kind=wp) dimension(size(omega,1),size(omega,2)),intent(in,out),depend(omega,omega) :: velocity_y
+    integer(kind=ip) dimension(2),intent(in) :: ghosts_vort
+    integer(kind=ip) dimension(2),intent(in) :: ghosts_velo
+  end subroutine solve_poisson_2d
+  subroutine solve_diffusion_2d(nudt,omega,ghosts_vort) ! in fftw2py.f90:fftw2py
+    real(kind=wp) intent(in) :: nudt
+    real(kind=wp) dimension(:,:),intent(in,out) :: omega
+    integer(kind=ip) dimension(2),intent(in) :: ghosts_vort
+  end subroutine solve_diffusion_2d
+  subroutine solve_poisson_3d(omega_x,omega_y,omega_z,velocity_x,velocity_y,velocity_z,ghosts_vort,ghosts_velo) ! in fftw2py.f90:fftw2py
+    real(kind=wp) dimension(:,:,:),intent(in) :: omega_x
+    real(kind=wp) dimension(:,:,:),intent(in) :: omega_y
+    real(kind=wp) dimension(:,:,:),intent(in) :: omega_z
+    real(kind=wp) dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(in,out),depend(omega_x,omega_y,omega_z) :: velocity_x
+    real(kind=wp) dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(in,out),depend(omega_x,omega_y,omega_z) :: velocity_y
+    real(kind=wp) dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(in,out),depend(omega_x,omega_y,omega_z) :: velocity_z
+    integer(kind=ip) dimension(3),intent(in) :: ghosts_vort
+    integer(kind=ip) dimension(3),intent(in) :: ghosts_velo
+  end subroutine solve_poisson_3d
+  subroutine solve_poisson_2d_c(omega,velocity_x,velocity_y) ! in fftw2py.f90:fftw2py
+    complex(kind=wp) dimension(:,:),intent(in) :: omega
+    complex(kind=wp) dimension(size(omega,1),size(omega,2)),intent(in,out),depend(omega,omega) :: velocity_x
+    complex(kind=wp) dimension(size(omega,1),size(omega,2)),intent(in,out),depend(omega,omega) :: velocity_y
+  end subroutine solve_poisson_2d_c
+  subroutine solve_poisson_3d_c(omega_x,omega_y,omega_z,velocity_x,velocity_y,velocity_z) ! in fftw2py.f90:fftw2py
+    complex(kind=wp) dimension(:,:,:),intent(in) :: omega_x
+    complex(kind=wp) dimension(:,:,:),intent(in) :: omega_y
+    complex(kind=wp) dimension(:,:,:),intent(in) :: omega_z
+    complex(kind=wp) dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(in,out),depend(omega_x,omega_y,omega_z) :: velocity_x
+    complex(kind=wp) dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(in,out),depend(omega_x,omega_y,omega_z) :: velocity_y
+    complex(kind=wp) dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(in,out),depend(omega_x,omega_y,omega_z) :: velocity_z
+  end subroutine solve_poisson_3d_c
+  subroutine solve_curl_diffusion_3d(nudt,velocity_x,velocity_y,velocity_z,omega_x,omega_y,omega_z,ghosts_velo,ghosts_vort) ! in fftw2py.f90:fftw2py
+    real(kind=wp) intent(in) :: nudt
+    real(kind=wp) dimension(:,:,:),intent(in) :: velocity_x
+    real(kind=wp) dimension(:,:,:),intent(in) :: velocity_y
+    real(kind=wp) dimension(:,:,:),intent(in) :: velocity_z
+    real(kind=wp) dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(in,out),depend(velocity_x,velocity_y,velocity_z) :: omega_x
+    real(kind=wp) dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(in,out),depend(velocity_x,velocity_y,velocity_z) :: omega_y
+    real(kind=wp) dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(in,out),depend(velocity_x,velocity_y,velocity_z) :: omega_z
+    integer(kind=ip) dimension(3),intent(in) :: ghosts_velo
+    integer(kind=ip) dimension(3),intent(in) :: ghosts_vort
+  end subroutine solve_curl_diffusion_3d
+  subroutine solve_diffusion_3d(nudt,omega_x,omega_y,omega_z,ghosts) ! in fftw2py.f90:fftw2py
+    real(kind=wp) intent(in) :: nudt
+    real(kind=wp) dimension(:,:,:),intent(in,out) :: omega_x
+    real(kind=wp) dimension(:,:,:),intent(in,out) :: omega_y
+    real(kind=wp) dimension(:,:,:),intent(in,out) :: omega_z
+    integer(kind=ip) dimension(3),intent(in) :: ghosts
+  end subroutine solve_diffusion_3d
+  subroutine projection_om_3d(omega_x,omega_y,omega_z,ghosts) ! in fftw2py.f90:fftw2py
+    real(kind=wp) dimension(:,:,:),intent(in,out) :: omega_x
+    real(kind=wp) dimension(:,:,:),intent(in,out) :: omega_y
+    real(kind=wp) dimension(:,:,:),intent(in,out) :: omega_z
+    integer(kind=ip) dimension(3),intent(in) :: ghosts
+  end subroutine projection_om_3d
+  subroutine multires_om_3d(dxf,dyf,dzf,omega_x,omega_y,omega_z,ghosts) ! in fftw2py.f90:fftw2py
+    real(kind=wp) intent(in) :: dxf
+    real(kind=wp) intent(in) :: dyf
+    real(kind=wp) intent(in) :: dzf
+    real(kind=wp) dimension(:,:,:),intent(in,out) :: omega_x
+    real(kind=wp) dimension(:,:,:),intent(in,out) :: omega_y
+    real(kind=wp) dimension(:,:,:),intent(in,out) :: omega_z
+    integer(kind=ip) dimension(3),intent(in) :: ghosts
+  end subroutine multires_om_3d
+  subroutine pressure_3d(pressure,ghosts) ! in fftw2py.f90:fftw2py
+    real(kind=wp) dimension(:,:,:),intent(in,out) :: pressure
+    integer(kind=ip) dimension(3),intent(in) :: ghosts
+  end subroutine pressure_3d
+  subroutine solve_curl_3d(velocity_x,velocity_y,velocity_z,omega_x,omega_y,omega_z,ghosts_velo,ghosts_vort) ! in fftw2py.f90:fftw2py
+    real(kind=wp) dimension(:,:,:),intent(in) :: velocity_x
+    real(kind=wp) dimension(:,:,:),intent(in) :: velocity_y
+    real(kind=wp) dimension(:,:,:),intent(in) :: velocity_z
+    real(kind=wp) dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(in,out),depend(velocity_x,velocity_y,velocity_z) :: omega_x
+    real(kind=wp) dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(in,out),depend(velocity_x,velocity_y,velocity_z) :: omega_y
+    real(kind=wp) dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(in,out),depend(velocity_x,velocity_y,velocity_z) :: omega_z
+    integer(kind=ip) dimension(3),intent(in) :: ghosts_velo
+    integer(kind=ip) dimension(3),intent(in) :: ghosts_vort
+  end subroutine solve_curl_3d
+  subroutine solve_curl_2d(velocity_x,velocity_y,omega_z,ghosts_velo,ghosts_vort) ! in fftw2py.f90:fftw2py
+    real(kind=wp) dimension(:,:),intent(in) :: velocity_x
+    real(kind=wp) dimension(:,:),intent(in) :: velocity_y
+    real(kind=wp) dimension(size(velocity_x,1),size(velocity_x,2)),intent(in,out),depend(velocity_x,velocity_x) :: omega_z
+    integer(kind=ip) dimension(2),intent(in) :: ghosts_velo
+    integer(kind=ip) dimension(2),intent(in) :: ghosts_vort
+  end subroutine solve_curl_2d
+  subroutine spectrum_3d(field,spectrum,wavelengths,ghosts,length) ! in fftw2py.f90:fftw2py
+    real(kind=wp) dimension(:,:,:),intent(in) :: field
+    real(kind=wp) dimension(:),intent(in,out) :: spectrum
+    real(kind=wp) dimension(:),intent(in,out) :: wavelengths
+    integer(kind=ip) dimension(3),intent(in) :: ghosts
+    real(kind=wp) intent(in) :: length
+  end subroutine spectrum_3d
+
+end module fftw2py
+
+! This file was auto-generated with f2py (version:2).
+! See http://cens.ioc.ee/projects/f2py2e/
diff --git a/hysop/operator/computational.py b/hysop/operator/computational.py
index 626436545..d02936b5a 100755
--- a/hysop/operator/computational.py
+++ b/hysop/operator/computational.py
@@ -8,6 +8,9 @@ from hysop.mpi.topology import Cartesian
 from hysop.tools.parameters import Discretization
 from hysop.tools.profiler import profile
 import hysop.tools.numpywrappers as npw
+from hysop import __FFTW_ENABLED__
+if __FFTW_ENABLED__:
+    from hysop.f2hysop import fftw2py
 
 
 class Computational(Operator):
@@ -259,7 +262,6 @@ class Computational(Operator):
         assert self._single_topo, 'All fields must use the same topology.'
         # Get local mesh parameters from fftw
         comm = self._mpis.comm
-        from hysop.f2hysop import fftw2py
         if build_topos:
             # In that case, self._discretization must be
             # a Discretization object, used for all fields.
diff --git a/hysop/operator/diffusion.py b/hysop/operator/diffusion.py
index c72f69564..bbdf1239f 100644
--- a/hysop/operator/diffusion.py
+++ b/hysop/operator/diffusion.py
@@ -1,8 +1,8 @@
 # -*- coding: utf-8 -*-
-"""
-@file diffusion.py
+"""Operator for diffusion problem.
+
+See :ref:`diffusion` in HySoP user guide.
 
-Operator for diffusion problem.
 
 """
 from hysop.operator.computational import Computational
@@ -15,8 +15,9 @@ from hysop import __FFTW_ENABLED__
 
 
 class Diffusion(Computational):
-    """
-    Diffusion operator
+    """ Diffusion of a field.
+
+
     \f{eqnarray*}
     \omega = Op(\omega)
     \f} with :
diff --git a/hysop/operator/discrete/diffusion_fft.py b/hysop/operator/discrete/diffusion_fft.py
index daa409eeb..f93ce20d3 100644
--- a/hysop/operator/discrete/diffusion_fft.py
+++ b/hysop/operator/discrete/diffusion_fft.py
@@ -1,37 +1,45 @@
 # -*- coding: utf-8 -*-
-"""
-@file diffusion_fft.py
-Discrete Diffusion operator using FFTW (fortran)
+"""Discrete Diffusion operator using FFTW (fortran)
 """
 try:
     from hysop.f2hysop import fftw2py
 except ImportError:
-    from hysop.fakef2py import fftw2py
+    msg = 'fftw package not available for your hysop install.'
+    msg += 'Try to recompile with WITH_FFTW=ON'
+    raise ImportError(msg)
+
 from hysop.operator.discrete.discrete import DiscreteOperator
 from hysop.constants import debug
 from hysop.tools.profiler import profile
 
 
 class DiffusionFFT(DiscreteOperator):
-    """
-    Discretized Poisson operator based on FFTW.
+    """Discretized Poisson operator based on FFTW.
     See details in hysop.operator.diffusion.
 
     """
     @debug
     def __init__(self, vorticity, viscosity, **kwds):
+        """Discrete diffusion operator, based on fftw solver.
+
+        Parameters
+        ----------
+        vorticity : :class:`~hysop.fields.discrete.DiscreteField`
+             vorticity field, in/out parameter
+        viscosity : double
+        kwds : base class arguments
         """
-        Constructor.
-        @param[in,out] vorticity :  discretisation of the field \f$ \omega \f$.
-        @param[in] viscosity : \f$\nu\f$, viscosity of the considered medium.
-        """
-        ## Discretisation of the solution field
+        # Discretisation of the solution field
         self.vorticity = vorticity
-        ## Viscosity.
+        # Viscosity.
         self.viscosity = viscosity
-
-        if self.vorticity.dimension == 1:
-            raise AttributeError("1D case not yet implemented.")
+        dim = self.vorticity.dimension
+        if dim == 3:
+            self._apply = self.apply_3d
+        elif dim == 2:
+            self._apply = self.apply_2d
+        else:
+            raise AttributeError(dim + "D case not yet implemented.")
         # Base class initialisation
         assert 'variables' not in kwds, 'variables parameter is useless.'
         super(DiffusionFFT, self).__init__(variables=[vorticity],
@@ -46,22 +54,22 @@ class DiffusionFFT(DiscreteOperator):
             "Missing dt value for diffusion computation."
         dt = simulation.time_step
         ghosts = self.vorticity.topology.ghosts()
+        self._apply(dt, ghosts)
 
-        if self.vorticity.dimension == 2:
-            self.vorticity.data = fftw2py.solve_diffusion_2d(
-                self.viscosity * dt, self.vorticity.data, ghosts)
+    def apply_2d(self, dt, ghosts):
+        """2d implementation of apply function"""
+        self.vorticity.data = fftw2py.solve_diffusion_2d(
+            self.viscosity * dt, self.vorticity.data, ghosts)
 
-        elif self.vorticity.dimension == 3:
-            self.vorticity.data[0], self.vorticity.data[1],\
-                self.vorticity.data[2] = \
-                fftw2py.solve_diffusion_3d(self.viscosity * dt,
-                                           self.vorticity.data[0],
-                                           self.vorticity.data[1],
-                                           self.vorticity.data[2],
-                                           ghosts)
-
-        else:
-            raise ValueError("invalid problem dimension")
+    def apply_3d(self, dt, ghosts):
+        """3d implementation of apply function"""
+        self.vorticity.data[0], self.vorticity.data[1],\
+            self.vorticity.data[2] = \
+            fftw2py.solve_diffusion_3d(self.viscosity * dt,
+                                       self.vorticity.data[0],
+                                       self.vorticity.data[1],
+                                       self.vorticity.data[2],
+                                       ghosts)
 
     def finalize(self):
         """
diff --git a/hysop/operator/tests/test_diffusion.py b/hysop/operator/tests/test_diffusion.py
index cc3e80f86..5f87c4b33 100755
--- a/hysop/operator/tests/test_diffusion.py
+++ b/hysop/operator/tests/test_diffusion.py
@@ -34,14 +34,11 @@ def computeVort2D(res, x, y, t):
 
 
 @testsenv.fftw_failed
-def test_Diffusion3D():
+def test_diffusion_3d():
+    """Vector field diffusion, 3d domain"""
     dom = pp.Box(length=LL)
-
-    # Fields
     vorticity = pp.Field(domain=dom, formula=computeVort,
                          name='Vorticity', is_vector=True)
-
-    # Definition of the Poisson operator
     diff = Diffusion(viscosity=0.3, vorticity=vorticity, discretization=d3D)
     diff.discretize()
     diff.setup()
@@ -53,14 +50,11 @@ def test_Diffusion3D():
 
 
 @testsenv.fftw_failed
-def test_Diffusion3D_2():
+def test_diffusion_3d_2():
+    """Vector field diffusion, 3d domain, """
     dom = pp.Box(length=LL)
-
-    # Fields
     vorticity = pp.Field(domain=dom, formula=computeVort,
                          name='Vorticity', is_vector=True)
-
-    # Definition of the Poisson operator
     diff = Diffusion(viscosity=0.3, variables={vorticity: d3D})
     diff.discretize()
     diff.setup()
@@ -72,13 +66,10 @@ def test_Diffusion3D_2():
 
 
 @testsenv.fftw_failed
-def test_Diffusion2D():
+def test_diffusion_2d():
+    """Vector field diffusion 2d domain"""
     dom = pp.Box(length=LL[:2])
-
-    # Fields
     vorticity = pp.Field(domain=dom, formula=computeVort2D, name='Vorticity')
-
-    # Definition of the Poisson operator
     diff = Diffusion(viscosity=0.3, vorticity=vorticity, discretization=d2D)
     diff.discretize()
     diff.setup()
@@ -87,4 +78,3 @@ def test_Diffusion2D():
     vorticity.initialize(topo=topo)
     diff.apply(simu)
     diff.finalize()
-
diff --git a/hysop/operators.py b/hysop/operators.py
index 9dfb7a2eb..c4cb5241d 100644
--- a/hysop/operators.py
+++ b/hysop/operators.py
@@ -13,6 +13,7 @@ import hysop.operator.hdf_io
 import hysop.operator.poisson
 import hysop.operator.advection
 import hysop.operator.energy_enstrophy
+import hysop.operator.diffusion
 Stretching = hysop.operator.stretching.Stretching
 Curl = hysop.operator.differential.Curl
 Grad = hysop.operator.differential.Grad
@@ -27,3 +28,4 @@ HDF_Writer = hysop.operator.hdf_io.HDF_Writer
 Poisson = hysop.operator.poisson.Poisson
 Advection = hysop.operator.advection.Advection
 EnergyEnstrophy = hysop.operator.energy_enstrophy.EnergyEnstrophy
+Diffusion = hysop.operator.diffusion.Diffusion
diff --git a/precision.conf.in b/precision.conf.in
new file mode 100644
index 000000000..7255565e4
--- /dev/null
+++ b/precision.conf.in
@@ -0,0 +1,37 @@
+!> Select float precision for the whole code.
+!! This is a generated file, do not edit.
+!! Usage :
+!! cmake -DDOUBLEPREC=ON ...
+!! with value = simple or value = double
+module precision
+
+  use mpi, only: MPI_DOUBLE_PRECISION, MPI_FLOAT
+  use ISO_FORTRAN_ENV
+  use iso_c_binding
+
+  implicit none
+  !> Floats precision
+  integer, parameter :: sp = REAL32
+  integer, parameter :: dp = 8!REAL64
+  !> Chosen precision (set during config. using -DDOUBLEPREC=... with cmake)
+  integer, parameter :: wp = @WORKING_PRECISION@
+  !> integer precision
+  integer, parameter :: ip = 8!int64
+  !> MPI type for float
+  integer, parameter :: MPI_REAL_WP = @MPI_PRECISION@
+  !> alias to 1. in wp
+  real(kind=wp), parameter :: ONE = 1.
+  !> alias to 0. in wp
+  real(kind=wp), parameter :: ZERO = 0.
+  !>  trick to identify coordinates in a more user-friendly way
+  integer, parameter :: c_X=1,c_Y=2,c_Z=3
+  !> MPI main communicator
+  integer :: main_comm
+  !> Rank of the mpi current process
+  integer :: rank ! current mpi-processus rank
+  !> i (sqrt(-1) ...)
+  complex(C_DOUBLE_COMPLEX), parameter :: Icmplx = cmplx(zero, one, kind=wp)
+  !> Pi constant
+  real(wp), parameter :: pi = 4.0*atan(one)
+
+end module precision
diff --git a/setup.py.in b/setup.py.in
index 1676f83c0..10ee92f14 100644
--- a/setup.py.in
+++ b/setup.py.in
@@ -150,6 +150,7 @@ def create_fortran_extension(name, pyf_file=None, src_dirs=None, sources=None,
     # we trust cmake for external libraries and
     # add them to linker, without using libraries option
     extra_link_args = hysop_link_libraries
+    print 'IN CREATE', sources
     ext_fort = Extension(name=name,
                          sources=sources,
                          f2py_options=f2py_options,
@@ -270,12 +271,17 @@ if enable_fortran is "ON":
     #####hysoplib = ['@HYSOP_LIBRARY_NAME@']
     f2py_options = ['--no-lower']
     fortran_src = set([])
+    # -- list of directories which contains fortran sources --
+    # path must be relative to hysop package directory
+    subdirs = ['common_f', ]
 
     ####fortran_src.add('f2py/parameters.f90')
     # -- fftw fortran sources --
     withfftw = "@WITH_FFTW@"
     if withfftw is "ON":
-        fortran_src.add('f2py/fftw2py.f90')
+        #        fortran_src.add('f2py/fftw2py.f90')
+        subdirs.append(os.path.join('numerics/fftw_f'))
+        #fortran_src.add('f2py/fftw2py.f90')
         fftwdir = '@FFTWLIB@'
         #hysoplib.append('fftw3')
         #hysoplib.append('fftw3_mpi')
@@ -291,12 +297,15 @@ if enable_fortran is "ON":
     for i in xrange(len(fortran_src)):
         fortran_src[i] = os.path.join(fortran_root, fortran_src[i])
 
+    precision_file = os.path.join('@GENERATED_FORTRAN_FILES_DIR@',
+                                  'precision.f95')
+    fortran_src.insert(0, precision_file)
+
+    print "OKOKOKOK", fortran_src
     # === Draft for future implementation of fortran interface ===
     # -- f2py signature file --
     pyf_file = os.path.join('@GENERATED_FORTRAN_FILES_DIR@', 'f2hysop.pyf')
-    # -- list of directories which contains fortran sources --
-    # path must be relative to hysop package directory
-    subdirs = ['common_f', ]
+    # change from relative path in subdirs to absolute path in num_dirs
     num_dirs = []
     for sd in subdirs:
         num_dirs.append(os.path.join(fortran_root, sd))
@@ -317,12 +326,6 @@ if enable_fortran is "ON":
     for ex in ext:
         ext_modules.append(ext[ex])
 
-else:
-    packages.append('hysop.fakef2py')
-    packages.append('hysop.fakef2py.scales2py')
-    packages.append('hysop.fakef2py.fftw2py')
-
-
 # --- C++ files and swig interface --
 
 if enable_cpp is "ON":
-- 
GitLab