diff --git a/HySoP/CMakeLists.txt b/HySoP/CMakeLists.txt
index ce87f558d8f86f634647cac65cbca5cbb43ccf61..d73ebe57e75ac1fe4a7de390d4af087f8693aac6 100644
--- a/HySoP/CMakeLists.txt
+++ b/HySoP/CMakeLists.txt
@@ -38,6 +38,7 @@ option(WITH_SCALES "compile/create parmesscales lib and link it with Parmes. Def
 option(WITH_FFTW "Link with fftw library (required for some Parmes solvers), default = ON" ON)
 option(WITH_GPU "Use of GPU (required for some Parmes solvers), default = ON" ON)
 option(WITH_MAIN_FORTRAN "Create an executable (test purpose) from fortran sources in src/main, linked with libparmes, default = ON" ON)
+option(DEBUG "Enable debug mode for Parmes (0:disabled, 1:verbose, 2:trace, 3:verbose+trace). Default = 0" 0)
 
 if(NOT WITH_LIB_FORTRAN)
   message(WARNING "You deactivate libparmes (fortran) generation. This will disable fftw and scales fonctionnalities.")
@@ -211,6 +212,7 @@ if(VERBOSE_MODE)
   message(STATUS " Project uses FFTW : ${WITH_FFTW}")
   message(STATUS " Project uses GPU : ${WITH_GPU}")
   message(STATUS " Python packages will be installed in ${${PROJECT_NAME}_INSTALL_DIR}")
+  message(STATUS " ${PROJECT_NAME} debug mode : ${DEBUG}")
   message(STATUS "====================== ======= ======================")
 endif()
 
diff --git a/HySoP/hysop/__init__.py.in b/HySoP/hysop/__init__.py.in
index 0308bed2ed1a09ab218fd7a4f18023195cdf3c4f..918cc3472a5a83f1ee43cd8a210feec0d4029662 100755
--- a/HySoP/hysop/__init__.py.in
+++ b/HySoP/hysop/__init__.py.in
@@ -8,8 +8,16 @@ on hybrid architectures (MPI-GPU)
 __version__ = 1.00
 __all__ = ['Box', 'CartesianTopology']#, 'ScalarField']
 
+# Compilation flags
+__MPI_ENABLED__ = "@USE_MPI@" is "ON"
+__GPU_ENABLED__ = "@WITH_GPU@" is "ON"
+__FFTW_ENABLED__ = "@WITH_FFTW@" is "ON"
+__SCALES_ENABLED__ = "@WITH_SCALES@" is "ON"
+__VERBOSE__ = "@DEBUG@" in ["1", "3"]
+__DEBUG__ = "@DEBUG@" in ["2", "3"]
+
 # MPI
-if("@USE_MPI@" is "ON"):
+if __MPI_ENABLED__:
     #import parmepy.mpi as mpi
     from mpi import topology
     CartesianTopology = topology.CartesianTopology
diff --git a/HySoP/hysop/constants.py b/HySoP/hysop/constants.py
index ae69f46eaa4970a31148eb2b19356de71576b6f4..75e0b3b263842f5c97a151594c203dd4a118a6d0 100644
--- a/HySoP/hysop/constants.py
+++ b/HySoP/hysop/constants.py
@@ -4,7 +4,7 @@
 Constant parameters required for the parmepy package (internal use).
 
 """
-
+from parmepy import __DEBUG__
 import numpy as np
 import math
 import mpi4py.MPI as MPI
@@ -32,3 +32,29 @@ PERIODIC = 0
 ## Directions string
 S_DIR = ["_X", "_Y", "_Z"]
 
+
+if __DEBUG__:
+#define debug decorator:
+    def debugdecorator(f):
+        ## function f is being decorated
+        def decorator(*args, **kw):
+            # Print informations on decorated function
+            if f.__name__ is '__new__':
+                fullclassname = args[0].__mro__[0].__module__ + '.'
+                fullclassname += args[0].__mro__[0].__name__
+                print 'Instanciate :', fullclassname,
+                print ' (Inherits from : ',
+                print [c.__name__ for c in args[0].__mro__[1:]], ')'
+            else:
+                print 'Call : ', f.__name__, 'of ', f.__module__
+            ## Calling f
+            r = f(*args, **kw)
+            if f.__name__ is '__new__':
+                print '       |-> : ', repr(r)
+            return r
+        return decorator
+else:
+#define empty debug decorator:
+    def debugdecorator(f):
+        return f
+debug = debugdecorator
diff --git a/HySoP/hysop/domain/box.py b/HySoP/hysop/domain/box.py
index cb85d61f5b1fba9609a8629a8897f982219d3b33..9975780556568a573f081add1ee29dca8c96a861 100644
--- a/HySoP/hysop/domain/box.py
+++ b/HySoP/hysop/domain/box.py
@@ -3,7 +3,8 @@
 Classes for box domains description.
 """
 from domain import Domain
-from parmepy.constants import np, PARMES_REAL, PARMES_INTEGER, PERIODIC
+from parmepy.constants import np, PARMES_REAL, PARMES_INTEGER, PERIODIC, \
+    debug
 
 
 class Box(Domain):
@@ -13,6 +14,7 @@ class Box(Domain):
     @todo Have different BC
     """
 
+    @debug
     def __init__(self, dimension=3, length=[1.0, 1.0, 1.0],
                  origin=[0., 0., 0.]):
         """
diff --git a/HySoP/hysop/domain/domain.py b/HySoP/hysop/domain/domain.py
index 759452b8dddbbc439edb11affa304385273ce2b9..7cfd811082e1fd82b1773fed3e1449204d6bebff 100644
--- a/HySoP/hysop/domain/domain.py
+++ b/HySoP/hysop/domain/domain.py
@@ -3,13 +3,19 @@
 Classes for physical domains description.
 """
 from abc import ABCMeta, abstractmethod
+from parmepy.constants import debug
 
 
-class Domain:
+class Domain(object):
     """Abstract description of physical domain. """
 
     __metaclass__ = ABCMeta
 
+    @debug
+    def __new__(cls, *args, **kw):
+        return object.__new__(cls, *args, **kw)
+
+    @debug
     @abstractmethod
     def __init__(self, dimension):
         """ Constructor
@@ -20,6 +26,7 @@ class Domain:
         ## Domain topologies
         self.topologies = []
 
+    @debug
     def addTopology(self, topo):
         try:
             idTopo = self.topologies.index(topo)
diff --git a/HySoP/hysop/fields/analytical.py b/HySoP/hysop/fields/analytical.py
index 5d6780c4c4c4753a456fafb83da1622cab22b7a9..c47729ff3f26680c511f710f3f2a50c2e7297f24 100644
--- a/HySoP/hysop/fields/analytical.py
+++ b/HySoP/hysop/fields/analytical.py
@@ -3,6 +3,7 @@
 
 Continuous variable description defined with an analytic formula.
 """
+from parmepy.constants import debug
 from parmepy.fields.continuous import Field
 
 
@@ -12,6 +13,7 @@ class AnalyticalField(Field):
 
     """
 
+    @debug
     def __init__(self, domain, formula, name="", vector=False):
         """
         Create an AnalyticalField from a formula.
@@ -36,6 +38,7 @@ class AnalyticalField(Field):
         """
         return self.formula(*pos)
 
+    @debug
     def initialize(self):
         if self._fieldId == -1:
             raise ValueError("Cannot initialise analytical field non discretized.")
diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
index 3993e32143be8558186b2aa5315467be59fb840d..4c11fd723d112036d3ca03ee5fd40b8de51c2b89 100644
--- a/HySoP/hysop/fields/continuous.py
+++ b/HySoP/hysop/fields/continuous.py
@@ -3,7 +3,8 @@
 
 Continuous variable description
 """
-from parmepy.constants import np, PARMES_INTEGER
+from parmepy import __VERBOSE__
+from parmepy.constants import np, PARMES_INTEGER, debug
 from parmepy.fields.scalar import ScalarField
 from parmepy.fields.vector import VectorField
 
@@ -11,6 +12,11 @@ from parmepy.fields.vector import VectorField
 class Field(object):
     """ A variable defined on a physical domain """
 
+    @debug
+    def __new__(cls, *args, **kw):
+        return object.__new__(cls, *args, **kw)
+
+    @debug
     def __init__(self, domain, name="?", vector=False):
         """
         Create a Field.
@@ -36,6 +42,7 @@ class Field(object):
         ## Is field is a vector field
         self.vector = vector
 
+    @debug
     def discretize(self, topo):
         """
         Discretization of the field on a topology.
@@ -43,11 +50,10 @@ class Field(object):
         @return discrete field and its index.
 
         """
-        print self.name,
         # Find if topology already use to discretize self
         try:
             fieldId = self.topologies.index(topo)
-            print "Already discretized with topology:"
+            #print "Already discretized with topology:"
             # topo already used
         except ValueError:
             # Create a new discrete field with topo
@@ -64,10 +70,11 @@ class Field(object):
                                          name=self.name + "_D_" + str(fieldId),
                                          idFromParent=fieldId))
             self.topologies.insert(fieldId, topo)
-            print "Discretized with topology:"
-        print self.topologies[fieldId]
+            #print "Discretized with topology:"
+        #print self.topologies[fieldId]
         return self.discreteField[fieldId], fieldId
 
+    @debug
     def initialize(self):
         """
         Initialize a field.
@@ -77,10 +84,12 @@ class Field(object):
         if self._fieldId == -1:
             raise ValueError("Cannot initialise analytical" +
                              " field non discretized.")
-        print self.name,
+        if __VERBOSE__:
+            print self.name,
         for dF in self.discreteField:
             dF.initialize()
-        print " .Done"
+        if __VERBOSE__:
+            print " .Done"
 
     def __str__(self):
         """ Field info display """
diff --git a/HySoP/hysop/fields/discrete.py b/HySoP/hysop/fields/discrete.py
index 5dcfe1393aa1281096ecd32ee4f05aa884c7b076..79d77b0998375315ae01533daed309bbd04bdaf1 100644
--- a/HySoP/hysop/fields/discrete.py
+++ b/HySoP/hysop/fields/discrete.py
@@ -4,6 +4,7 @@
 Discretes variables (scalar and vectors) descriptions.
 """
 from abc import ABCMeta, abstractmethod
+from parmepy.constants import debug
 
 
 class DiscreteField(object):
@@ -11,6 +12,11 @@ class DiscreteField(object):
 
     __metaclass__ = ABCMeta
 
+    @debug
+    def __new__(cls, *args, **kw):
+        return object.__new__(cls, *args, **kw)
+
+    @debug
     @abstractmethod
     def __init__(self, parent, topology=None, name="?", idFromParent=None):
         """ Creates a ScalaField.
@@ -54,6 +60,7 @@ class DiscreteField(object):
         """ Access to the content of the field."""
         raise ValueError("Abstract method")
 
+    @debug
     def initialize(self, formula=None):
         """Initialize values with given formula."""
         raise ValueError("Abstract method")
diff --git a/HySoP/hysop/fields/gpu_discrete.py b/HySoP/hysop/fields/gpu_discrete.py
index 7911790ab6dca823c164ad08778801f0a3959433..b1c912315775a32fb864a0a9dfec9be90a663368 100644
--- a/HySoP/hysop/fields/gpu_discrete.py
+++ b/HySoP/hysop/fields/gpu_discrete.py
@@ -1,6 +1,7 @@
 """
 @package parmepy.fields.gpu_discrete
 """
+from parmepy import __VERBOSE__
 from parmepy.constants import ORDER, np
 from parmepy.gpu import cl, PARMES_REAL_GPU
 from parmepy.fields.vector import VectorField
@@ -16,7 +17,6 @@ class GPUDiscreteField(object):
     def __init__(self, queue):
         self.queue = queue
         self.gpu_data = None
-        print self.name, "buffer allocation "
         self.mem_size = 0
 
 
@@ -38,6 +38,11 @@ class GPUVectorField(VectorField, GPUDiscreteField):
                                          cl.mem_flags.READ_WRITE,
                                          size=self.data[d].nbytes)
             self.mem_size += self.gpu_data[d].size
+        if __VERBOSE__:
+            print self.name, self.mem_size, "Bytes (",
+            print self.mem_size / (1024 ** 2), "MB)"
+        else:
+            print ".",
 
     @classmethod
     def fromVectorField(cls, queue, vfield, precision=PARMES_REAL_GPU):
@@ -51,6 +56,11 @@ class GPUVectorField(VectorField, GPUDiscreteField):
                                            cl.mem_flags.READ_WRITE,
                                            size=vfield.data[d].nbytes)
             vfield.mem_size += vfield.gpu_data[d].size
+        if __VERBOSE__:
+            print vfield.name, vfield.mem_size, "Bytes (",
+            print vfield.mem_size / (1024 ** 2), "MB)"
+        else:
+            print ".",
 
     def toDevice(self):
         """
@@ -89,9 +99,12 @@ class GPUVectorField(VectorField, GPUDiscreteField):
             if e is not None:
                 time += (e.profile.end - e.profile.start) * 1e-9
         self.data_on_device = True
-        print self.mem_size, "Bytes transfered at ",
-        print "{0:.3f} GBytes/sec".format(
-            self.mem_size / (time * 1024 ** 3))
+        if __VERBOSE__:
+            print self.mem_size, "Bytes transfered at ",
+            print "{0:.3f} GBytes/sec".format(
+                self.mem_size / (time * 1024 ** 3))
+        else:
+            print "Transfer host->device"
         return self.mem_size, time
 
     def toHost(self):
@@ -131,10 +144,13 @@ class GPUVectorField(VectorField, GPUDiscreteField):
                 if e is not None:
                     time += (e.profile.end - e.profile.start) * 1e-9
             self.data_on_device = False
-            print self.mem_size, "Bytes transfered at ",
-            print "{0:.3f} GBytes/sec".format(
-                self.mem_size / (time * 1024 ** 3))
-            return self.mem_size, time
+            if __VERBOSE__:
+                print self.mem_size, "Bytes transfered at ",
+                print "{0:.3f} GBytes/sec".format(
+                    self.mem_size / (time * 1024 ** 3))
+                return self.mem_size, time
+            else:
+                print 'Transfer device->host'
 
 
 class GPUScalarField(ScalarField, GPUDiscreteField):
@@ -149,7 +165,6 @@ class GPUScalarField(ScalarField, GPUDiscreteField):
                                              name, idFromParent)
         self.queue = queue
         self.gpu_data = None
-        print self.name, "buffer allocation "
         self.mem_size = 0
         self.data = np.asarray(self.data,
                                dtype=precision, order=ORDER)
@@ -157,6 +172,11 @@ class GPUScalarField(ScalarField, GPUDiscreteField):
                                   cl.mem_flags.READ_WRITE,
                                   size=self.data.nbytes)
         self.mem_size = self.gpu_data.size
+        if __VERBOSE__:
+            print self.name, self.mem_size, "Bytes (",
+            self.mem_size / (1024 ** 2), "MB)"
+        else:
+            print '.',
 
     @classmethod
     def fromScalarField(cls, queue, sfield, precision=PARMES_REAL_GPU):
@@ -168,6 +188,11 @@ class GPUScalarField(ScalarField, GPUDiscreteField):
                                     cl.mem_flags.READ_WRITE,
                                     size=sfield.data.nbytes)
         sfield.mem_size = sfield.gpu_data.size
+        if __VERBOSE__:
+            print sfield.name, sfield.mem_size, "Bytes (",
+            print sfield.mem_size / (1024 ** 2), "MB)"
+        else:
+            print '.',
 
     def toDevice(self):
         """
@@ -198,9 +223,12 @@ class GPUScalarField(ScalarField, GPUDiscreteField):
             if e is not None:
                 time += (e.profile.end - e.profile.start) * 1e-9
         self.data_on_device = True
-        print self.mem_size, "Bytes transfered at ",
-        print "{0:.3f} GBytes/sec".format(
-            self.mem_size / (time * 1024 ** 3))
+        if __VERBOSE__:
+            print self.mem_size, "Bytes transfered at ",
+            print "{0:.3f} GBytes/sec".format(
+                self.mem_size / (time * 1024 ** 3))
+        else:
+            print 'Transfer host->device'
         return self.mem_size, time
 
     def toHost(self):
@@ -227,7 +255,10 @@ class GPUScalarField(ScalarField, GPUDiscreteField):
                 if e is not None:
                     time += (e.profile.end - e.profile.start) * 1e-9
             self.data_on_device = False
-            print self.mem_size, "Bytes transfered at ",
-            print "{0:.3f} GBytes/sec".format(
-                self.mem_size / (time * 1024 ** 3))
+            if __VERBOSE__:
+                print self.mem_size, "Bytes transfered at ",
+                print "{0:.3f} GBytes/sec".format(
+                    self.mem_size / (time * 1024 ** 3))
+            else:
+                print 'Transfer device->host'
             return self.mem_size, time
diff --git a/HySoP/hysop/fields/scalar.py b/HySoP/hysop/fields/scalar.py
index 18ce6931d7708c2b235d2dc2a5a094960d3183fe..b30d040ff8954759cb8dcb5c45ca7db527e51405 100644
--- a/HySoP/hysop/fields/scalar.py
+++ b/HySoP/hysop/fields/scalar.py
@@ -3,7 +3,8 @@
 
 Discrete variable scalar description.
 """
-from parmepy.constants import np, ORDER, PARMES_REAL
+from parmepy import __VERBOSE__
+from parmepy.constants import np, ORDER, PARMES_REAL, debug
 from parmepy.fields.discrete import DiscreteField
 
 
@@ -14,6 +15,7 @@ class ScalarField(DiscreteField):
 
     """
 
+    @debug
     def __init__(self, parent, topology=None, name="?", idFromParent=None):
         """Constructor for a scalar field"""
         super(ScalarField, self).__init__(parent, topology,
@@ -42,6 +44,7 @@ class ScalarField(DiscreteField):
         """
         self.data.__setitem__(i, value)
 
+    @debug
     def initialize(self, formula=None):
         """
         Initialize values with given formula.
@@ -55,7 +58,8 @@ class ScalarField(DiscreteField):
         where x,y,z are point coordinates.
         """
         if formula is not None:
-            print "...",
+            if __VERBOSE__:
+                print "...",
             v_formula = np.vectorize(formula)
             if self.dimension == 3:
                 self.data = v_formula(*self.topology.mesh.coords)
@@ -65,7 +69,8 @@ class ScalarField(DiscreteField):
                 self.data = v_formula(self.topology.mesh.coords)
             self.contains_data = True
         else:
-            print "No formula",
+            if __VERBOSE__:
+                print "No formula",
 
     def __str__(self):
         """ Display class information. """
diff --git a/HySoP/hysop/fields/vector.py b/HySoP/hysop/fields/vector.py
index 02ffb0c7c54f44bd19821a7ae1a7c2b4ed642749..b61bb5a85f4c320f21bf95d9b46198ffcdc291bd 100644
--- a/HySoP/hysop/fields/vector.py
+++ b/HySoP/hysop/fields/vector.py
@@ -3,7 +3,8 @@
 
 Discrete variable scalar description.
 """
-from parmepy.constants import np, ORDER, PARMES_REAL
+from parmepy import __VERBOSE__
+from parmepy.constants import np, ORDER, PARMES_REAL, debug
 from parmepy.fields.discrete import DiscreteField
 
 
@@ -14,6 +15,7 @@ class VectorField(DiscreteField):
     Vector field is composed of several ScalarFields
     """
 
+    @debug
     def __init__(self, parent, topology=None, name="?", idFromParent=None):
         """Constructor for a vector field"""
         super(VectorField, self).__init__(parent, topology,
@@ -71,6 +73,7 @@ class VectorField(DiscreteField):
         elif len(i) > len(self.data):
             self.data[i[0]].__setitem__(tuple(i[1:]), value)
 
+    @debug
     def initialize(self, formula=None):
         """
         Initialize values with given formula.
@@ -84,7 +87,8 @@ class VectorField(DiscreteField):
         where x,y,z are point coordinates.
         """
         if formula is not None:
-            print "...",
+            if __VERBOSE__:
+                print "...",
             v_formula = np.vectorize(formula)
             if self.dimension == 3:
                 self.data[0], self.data[1], self.data[2] = \
@@ -105,7 +109,8 @@ class VectorField(DiscreteField):
 
             self.contains_data = True
         else:
-            print "No formula",
+            if __VERBOSE__:
+                print "No formula",
 
     def __str__(self):
         """ Display class information """
diff --git a/HySoP/hysop/gpu/tools.py b/HySoP/hysop/gpu/tools.py
index 5b5d8c115be72129b0829820337ba60b3ac3b270..fbec5f3e8242583f6e2d90fe7a2fd60adee12f50 100644
--- a/HySoP/hysop/gpu/tools.py
+++ b/HySoP/hysop/gpu/tools.py
@@ -3,11 +3,13 @@
 
 Tools for gpu management.
 """
+from parmepy import __VERBOSE__
 from parmepy.gpu import cl, PARMES_REAL_GPU
 
 
 def get_opencl_environment(platform_id, device_id, device_type):
-    print "=== OpenCL environment ==="
+    if __VERBOSE__:
+        print "=== OpenCL environment ==="
     #Get platform.
     try:
         ## OpenCL platform
@@ -17,9 +19,10 @@ def get_opencl_environment(platform_id, device_id, device_type):
         print " Only ", len(cl.get_platforms()), " available.",
         print " Getting defalut platform. "
         platform = cl.get_platforms()[0]
-    print "  Platform   "
-    print "  - Name       :", platform.name
-    print "  - Version    :", platform.version
+    if __VERBOSE__:
+        print "  Platform   "
+        print "  - Name       :", platform.name
+        print "  - Version    :", platform.version
     #Get device.
     try:
         ## OpenCL device
@@ -31,16 +34,17 @@ def get_opencl_environment(platform_id, device_id, device_type):
     except AttributeError as e:
         print "AttributeError:", e
         device = cl.create_some_context().devices[0]
-    print "  Device"
-    print "  - Name                :",
-    print device.name
-    print "  - Type                :",
-    print cl.device_type.to_string(device.type)
-    print "  - C Version           :",
-    print device.opencl_c_version
-    print "  - Global mem size     :",
-    print device.global_mem_size / (1024 ** 3), "GB"
-    print "===\n"
+    if __VERBOSE__:
+        print "  Device"
+        print "  - Name                :",
+        print device.name
+        print "  - Type                :",
+        print cl.device_type.to_string(device.type)
+        print "  - C Version           :",
+        print device.opencl_c_version
+        print "  - Global mem size     :",
+        print device.global_mem_size / (1024 ** 3), "GB"
+        print "===\n"
     #Creates GPU Context
     ## OpenCL context
     ctx = cl.Context([device])
@@ -76,10 +80,12 @@ def explore():
 
 def _check_double_precision_capability(device, prec):
     ## Precision supported
-    print " Precision capability  ",
+    if __VERBOSE__:
+        print " Precision capability  ",
     capable = True
     if prec is PARMES_REAL_GPU:
-        print "for Single Precision: "
+        if __VERBOSE__:
+            print "for Single Precision: "
         for v in ['DENORM', 'INF_NAN',
                   'ROUND_TO_NEAREST', 'ROUND_TO_ZERO', 'ROUND_TO_INF',
                   'FMA', 'CORRECTLY_ROUNDED_DIVIDE_SQRT', 'SOFT_FLOAT']:
@@ -87,14 +93,17 @@ def _check_double_precision_capability(device, prec):
                 if eval('(device.single_fp_config &' +
                         ' cl.device_fp_config.' +
                         v + ') == cl.device_fp_config.' + v):
-                    print v
+                    if __VERBOSE__:
+                        print v
             except AttributeError as ae:
                 if v is 'CORRECTLY_ROUNDED_DIVIDE_SQRT':
                     capable = False
-                print v, 'is not supported in OpenCL C 1.2.\n',
-                print '   Exception catched : ', ae
+                if __VERBOSE__:
+                    print v, 'is not supported in OpenCL C 1.2.\n',
+                    print '   Exception catched : ', ae
     else:
-        print "for Double Precision: "
+        if __VERBOSE__:
+            print "for Double Precision: "
         if device.double_fp_config > 0:
             for v in ['DENORM', 'INF_NAN', 'FMA',
                       'ROUND_TO_NEAREST', 'ROUND_TO_ZERO', 'ROUND_TO_INF',
@@ -103,12 +112,14 @@ def _check_double_precision_capability(device, prec):
                     if eval('(device.double_fp_config &' +
                             ' cl.device_fp_config.' +
                             v + ') == cl.device_fp_config.' + v):
-                        print v
+                        if __VERBOSE__:
+                            print v
                 except AttributeError as ae:
                     if v is 'CORRECTLY_ROUNDED_DIVIDE_SQRT':
                         capable = False
-                    print  v, 'is supported in OpenCL C 1.2.\n',
-                    print '   Exception catched : ', ae
+                    if __VERBOSE__:
+                        print  v, 'is supported in OpenCL C 1.2.\n',
+                        print '   Exception catched : ', ae
         else:
             raise ValueError("Double Precision is not supported by device")
     return capable
diff --git a/HySoP/hysop/mpi/mesh.py b/HySoP/hysop/mpi/mesh.py
index 5ce63966df82ffd9c3cc2fc8b1943a9a1959158c..02ff8aba5a50554147190926ff0859807370efcf 100644
--- a/HySoP/hysop/mpi/mesh.py
+++ b/HySoP/hysop/mpi/mesh.py
@@ -5,7 +5,8 @@ Cartesian mesh class for local mesh.
 @todo rename to CartesianMesh
 @todo write tests in parmepy/domain/tests/test_mesh.py
 """
-from parmepy.constants import np, PARMES_REAL, PARMES_INTEGER
+from parmepy.constants import np, PARMES_REAL, PARMES_INTEGER, \
+    debug
 
 
 class LocalMesh(object):
@@ -60,6 +61,12 @@ class SubMesh(object):
     Description of a local mesh, on one process of
     the topology.
     """
+
+    @debug
+    def __new__(cls, *args, **kw):
+        return object.__new__(cls, *args, **kw)
+
+    @debug
     def __init__(self, topo, g_start, resolution):
 
         ## Topology that creates (and owns) this mesh
diff --git a/HySoP/hysop/mpi/topology.py b/HySoP/hysop/mpi/topology.py
index 6523dcd3a1aaa2aea7d3ec34cfdb18738f633b4e..7cd6ee7f0172c00a2262df69f1b1b9db1dc07b23 100644
--- a/HySoP/hysop/mpi/topology.py
+++ b/HySoP/hysop/mpi/topology.py
@@ -4,7 +4,7 @@ MPI Topologies
 @todo merge CartesianTopology and Topology
 """
 
-from parmepy.constants import ORDER, np, PARMES_INTEGER, XDIR, YDIR, ZDIR, PARMES_INDEX
+from parmepy.constants import ORDER, np, PARMES_INTEGER, XDIR, YDIR, ZDIR, PARMES_INDEX, debug
 from parmepy.mpi.mesh import LocalMesh, SubMesh
 from main_var import main_comm, MPI, main_size
 from itertools import count
@@ -137,9 +137,14 @@ class Cartesian(object):
 
     """
 
+    @debug
+    def __new__(cls, *args, **kw):
+        return object.__new__(cls, *args, **kw)
+
     ## Counter of topology.Cartesian instances to set idTopo.
     __ids = count(0)
 
+    @debug
     def __init__(self, domain, dim, globalMeshResolution,
                  comm=main_comm, periods=None, ghosts=None, cutdir=None):
         ## Associated domain
@@ -270,6 +275,7 @@ class Cartesian(object):
         ##L_end = self.localGridResolution.copy() - self.ghosts[:] - 1
         self.mesh = None
 
+    @debug
     def setUp(self):
         """Topology set up.
 
diff --git a/HySoP/hysop/numerics/gpu_kernel.py b/HySoP/hysop/numerics/gpu_kernel.py
index 16459e4f3363894cc524912ae620ce2e60b1ef53..062a6428a87fa332b6b8e463698337876dc7ec77 100644
--- a/HySoP/hysop/numerics/gpu_kernel.py
+++ b/HySoP/hysop/numerics/gpu_kernel.py
@@ -1,6 +1,8 @@
 """
 @package parmepy.numerics.gpu_kernel
 """
+from parmepy.constants import debug
+from parmepy import __VERBOSE__
 from parmepy.numerics.method import NumMethod
 from parmepy.gpu import cl
 
@@ -11,6 +13,7 @@ class KernelListLauncher:
 
     Manage launching of OpenCL kernels as a list.
     """
+    @debug
     def __init__(self, kernel, queue, gsize, lsize=None):
         """
         Create a kernel list launcher.
@@ -21,7 +24,7 @@ class KernelListLauncher:
         """
         ## OpenCL Kernel list
         self.kernel = kernel
-        print [k.function_name for k in self.kernel]
+        #print [k.function_name for k in self.kernel]
         ## OpenCL command queue
         self.queue = queue
         ## OpenCL global size index.
@@ -30,7 +33,8 @@ class KernelListLauncher:
         self.local_size = lsize
         self.call_number = [0 for k in self.kernel]
 
-    def launch(self, d, *args):
+    @debug
+    def __call__(self, d, *args):
         """
         Launch a kernel.
 
@@ -44,6 +48,7 @@ class KernelListLauncher:
         return KernelListLauncher.launch_sizes_in_args(
             self, d, self.global_size[d], self.local_size[d], *args)
 
+    @debug
     def launch_sizes_in_args(self, d, *args):
         """
         Launch a kernel.
@@ -54,7 +59,8 @@ class KernelListLauncher:
 
         Opencl global and local sizes are given in args.
         """
-        #print self.kernel[d].function_name, args[0], args[1]
+        if __VERBOSE__:
+            print self.kernel[d].function_name, args[0], args[1]
         self.call_number[d] += 1
         evt = self.kernel[d](self.queue, *args)
         return evt
@@ -78,6 +84,7 @@ class KernelLauncher(KernelListLauncher):
     Manage launching of one OpenCL kernel as a KernelListLauncher
     with a list of one kernel.
     """
+    @debug
     def __init__(self, kernel, queue, gsize, lsize):
         """
         Create a KernelLauncher.
@@ -91,6 +98,7 @@ class KernelLauncher(KernelListLauncher):
         """
         KernelListLauncher.__init__(self, [kernel], queue, [gsize], [lsize])
 
+    @debug
     def launch_sizes_in_args(self, *args):
         """
         Launch the kernel.
@@ -106,7 +114,8 @@ class KernelLauncher(KernelListLauncher):
         """Wrapping OpenCL queue.finish() method."""
         self.queue.finish()
 
-    def launch(self, *args):
+    @debug
+    def __call__(self, *args):
         """
         Launch the kernel.
 
@@ -116,7 +125,7 @@ class KernelLauncher(KernelListLauncher):
         OpenCL global size and local sizes are not given in args.
         Class member are used.
         """
-        return KernelListLauncher.launch(self, 0, *args)
+        return KernelListLauncher.__call__(self, 0, *args)
 
     def function_name(self):
         """Prints OpenCL Kernel function name informations"""
diff --git a/HySoP/hysop/numerics/method.py b/HySoP/hysop/numerics/method.py
index bfe58adc62d97ed035d8d925afa19757dd7ab671..3c6eac9a40daaee424326707fa585c51184da018 100644
--- a/HySoP/hysop/numerics/method.py
+++ b/HySoP/hysop/numerics/method.py
@@ -4,6 +4,7 @@
 Numerical method interface.
 """
 from abc import ABCMeta, abstractmethod
+from parmepy.constants import debug
 
 
 class NumMethod(object):
@@ -11,11 +12,17 @@ class NumMethod(object):
 
     __metaclass__ = ABCMeta
 
+    @debug
+    def __new__(cls, *args, **kw):
+        return object.__new__(cls, *args, **kw)
+
+    @debug
     @abstractmethod
     def __init__(self):
         """ Creates a numerical method."""
         self.name = ''
 
+    @debug
     @abstractmethod
     def __call__(self):
         """Call the method"""
diff --git a/HySoP/hysop/operator/advection.py b/HySoP/hysop/operator/advection.py
index e3a282aa4a2dc0bc4c759634d9017832f7f48527..cb8930c747a47ae3bae8bd6191d6af654078e85f 100644
--- a/HySoP/hysop/operator/advection.py
+++ b/HySoP/hysop/operator/advection.py
@@ -3,9 +3,9 @@
 
 Advection operator representation.
 """
+from parmepy.constants import debug
 from parmepy.operator.continuous import Operator
 from parmepy.mpi.topology import Cartesian
-from parmepy.fields.scalar import ScalarField
 
 
 class Advection(Operator):
@@ -13,6 +13,7 @@ class Advection(Operator):
     Advection operator representation.
     """
 
+    @debug
     def __init__(self, velocity, advectedFields, resolutions=None, method='',
                  **other_config):
         """
@@ -29,12 +30,13 @@ class Advection(Operator):
         Operator.__init__(self, v)
         ## Transport velocity
         self.velocity = velocity
-        ## Transported scalar
+        ## Transported fields
         self.advectedFields = advectedFields
         self.resolutions = resolutions
         self.method = method
         self.config = other_config
 
+    @debug
     def setUp(self):
         """
         Transport operator discretization method.
@@ -42,7 +44,6 @@ class Advection(Operator):
 
         """
         Operator.setUp(self)
-        print self.method
         #variables discretization
         for v in self.variables:
             topoId, topo = self.domain.addTopology(
@@ -52,9 +53,10 @@ class Advection(Operator):
             self.discreteFieldId[v] = dfId
         if self.discreteOperator is None:
             if self.method.find('gpu_2k') >= 0:
-                if not isinstance(self.advectedFields.discreteField[
-                        self.discreteFieldId[self.advectedFields]],
-                        ScalarField):
+                if self.advectedFields.vector:
+                #if not isinstance(self.advectedFields.discreteField[
+                #        self.discreteFieldId[self.advectedFields]],
+                #        ScalarField):
                     raise ValueError("Not implemented yet. GPU advection" +
                                      " is only for scalar fields")
                 from parmepy.operator.gpu_particle_advection_2k \
@@ -63,9 +65,10 @@ class Advection(Operator):
                     GPUParticleAdvection2k(self, method=self.method,
                                            **self.config)
             elif self.method.find('gpu_1k') >= 0:
-                if not isinstance(self.advectedFields.discreteField[
-                        self.discreteFieldId[self.advectedFields]],
-                        ScalarField):
+                if self.advectedFields.vector:
+                #if not isinstance(self.advectedFields.discreteField[
+                #        self.discreteFieldId[self.advectedFields]],
+                #        ScalarField):
                     raise ValueError("Not implemented yet. GPU advection" +
                                      " is only for scalar fields")
                 from parmepy.operator.gpu_particle_advection_1k \
@@ -74,6 +77,8 @@ class Advection(Operator):
                     GPUParticleAdvection1k(self, method=self.method,
                                            **self.config)
             elif self.method.find('scales') >= 0:
+                if not self.advectedFields.dimension == 3:
+                    raise ValueError("Scales Advection not implemented in 2D.")
                 from parmepy.operator.scales_advection \
                     import ScalesAdvection
                 self.discreteOperator = \
@@ -81,7 +86,7 @@ class Advection(Operator):
                                     **self.config)
             else:
                 print "Using default advection operator"
-                if not isinstance(self.advectedFields, ScalarField):
+                if self.advectedFields.vector:
                     raise ValueError("Not implemented yet. default advection" +
                                      " is only for scalar fields")
                 from parmepy.operator.particle_advection \
diff --git a/HySoP/hysop/operator/continuous.py b/HySoP/hysop/operator/continuous.py
index bbb873794a63d2c58f8ab78f21da79e4e8c794b1..11b96a3dbc482e050bbf5b4ef0e950156b9ea1e4 100644
--- a/HySoP/hysop/operator/continuous.py
+++ b/HySoP/hysop/operator/continuous.py
@@ -4,6 +4,7 @@
 Operator representation.
 """
 from abc import ABCMeta, abstractmethod
+from parmepy.constants import debug
 
 
 class Operator(object):
@@ -13,6 +14,11 @@ class Operator(object):
 
     __metaclass__ = ABCMeta
 
+    @debug
+    def __new__(cls, *args, **kw):
+        return object.__new__(cls, *args, **kw)
+
+    @debug
     @abstractmethod
     def __init__(self, variables):
         """
@@ -39,10 +45,12 @@ class Operator(object):
             if self.variables.count(v) == 0:
                 self.variables.append(v)
 
+    @debug
     def setUp(self):
         if self.resolutions is None:
             raise ValueError("Advection setUp: no resolution given.")
 
+    @debug
     def finalize(self):
         """
         Method called on simulation finalization.
@@ -50,6 +58,7 @@ class Operator(object):
         """
         self.discreteOperator.finalize()
 
+    @debug
     def apply(self, *args):
         """
         Apply operator.
diff --git a/HySoP/hysop/operator/diffusion.py b/HySoP/hysop/operator/diffusion.py
index 2cd398b161228539c87ab8b23997dd1c83adf3f3..1bd91d4c216531d49471aef99671954517acd4df 100644
--- a/HySoP/hysop/operator/diffusion.py
+++ b/HySoP/hysop/operator/diffusion.py
@@ -12,6 +12,7 @@ class Diffusion(Operator):
     Diffusion operator representation using FFT
     """
 
+    @debug
     def __init__(self, velocity=None, vorticity=None, viscosity=0.1,
                  resolutions=None, method='',
                  **other_config):
@@ -31,6 +32,7 @@ class Diffusion(Operator):
         self.method = method
         self.config = other_config
 
+    @debug
     def setUp(self):
         """
         Diffusion operator discretization method.
diff --git a/HySoP/hysop/operator/diffusion_fft.py b/HySoP/hysop/operator/diffusion_fft.py
index e696f5dc300cd7d78fe7be1bb1e123890e79c488..d121a58e8a4f782dcbf1619ed99f6ff75db2ba8f 100644
--- a/HySoP/hysop/operator/diffusion_fft.py
+++ b/HySoP/hysop/operator/diffusion_fft.py
@@ -14,6 +14,7 @@ class DiffusionFFT(DiscreteOperator):
     Operator representation.
     """
 
+    @debug
     def __init__(self, diff, method='',
                  viscosity=0.1):
         """
@@ -30,6 +31,7 @@ class DiffusionFFT(DiscreteOperator):
         self.viscosity = viscosity
         self.compute_time = 0.
 
+    @debug
     def setUp(self, method='', topology=None):
         ## La topologie se recupere depuis les variables discretes
         self.topology = self.velocity.topology
@@ -38,6 +40,7 @@ class DiffusionFFT(DiscreteOperator):
         self.resolution = self.topology.mesh.resolution
         self.dim = self.topology.dim
 
+    @debug
     def apply(self, t, dt, ite):
         """
         Apply operator.
diff --git a/HySoP/hysop/operator/discrete.py b/HySoP/hysop/operator/discrete.py
index cff87465103625444484b37c9b3d46f66a48e64f..f3767d91ce6e10553a72a76440887594963ae8aa 100644
--- a/HySoP/hysop/operator/discrete.py
+++ b/HySoP/hysop/operator/discrete.py
@@ -4,6 +4,7 @@
 Discrete operator representation.
 """
 from abc import ABCMeta, abstractmethod
+from parmepy.constants import debug
 
 
 class DiscreteOperator(object):
@@ -13,6 +14,11 @@ class DiscreteOperator(object):
 
     __metaclass__ = ABCMeta
 
+    @debug
+    def __new__(cls, *args, **kw):
+        return object.__new__(cls, *args, **kw)
+
+    @debug
     @abstractmethod
     def __init__(self, variables, method=None):
         """
@@ -40,6 +46,7 @@ class DiscreteOperator(object):
         raise NotImplementedError("Need to override method in " +
                                   "a subclass of DiscreteOperator")
 
+    @debug
     @abstractmethod
     def apply(self, t, dt, ite):
         """
@@ -49,6 +56,7 @@ class DiscreteOperator(object):
         raise NotImplementedError("Need to override method in " +
                                   "a subclass of DiscreteOperator")
 
+    @debug
     def finalize(self):
         """
         Method called on simulation finalization.
diff --git a/HySoP/hysop/operator/gpu_particle_advection.py b/HySoP/hysop/operator/gpu_particle_advection.py
index f26f427525fa61aa29a414aaee1d9498d5d417f5..f82db3717bb56a7ea91615b7496377a9049d7179 100644
--- a/HySoP/hysop/operator/gpu_particle_advection.py
+++ b/HySoP/hysop/operator/gpu_particle_advection.py
@@ -4,7 +4,8 @@
 Discrete advection representation
 """
 from abc import ABCMeta, abstractmethod
-from parmepy.constants import np
+from parmepy import __VERBOSE__
+from parmepy.constants import np, debug
 from parmepy.operator.particle_advection import ParticleAdvection
 from parmepy.fields.continuous import Field
 from parmepy.fields.gpu_discrete import GPUVectorField
@@ -25,6 +26,7 @@ class GPUParticleAdvection(ParticleAdvection):
     """
     __metaclass__ = ABCMeta
 
+    @debug
     @abstractmethod
     def __init__(self, advec,
                  platform_id=0, device_id=0,
@@ -46,7 +48,7 @@ class GPUParticleAdvection(ParticleAdvection):
         self.user_gpu_src = src
         self.gpu_precision = precision
         self.num_method = None
-        self.resolution = self.scalar.topology.mesh.resolution
+        self.resolution = self.advectedFields.topology.mesh.resolution
         self.workItemNumber = None
         self.platform, self.device, self.ctx, self.queue = \
             get_opencl_environment(platform_id, device_id, device_type)
@@ -88,19 +90,22 @@ class GPUParticleAdvection(ParticleAdvection):
             lwi = (int(workItemNumber), 1)
         return workItemNumber, gwi, lwi
 
+    @debug
     def setUp(self):
         self.workItemNumber, self.gwi, self.lwi = self._set_WorkItems()
-        print "Work-items basic config : ",
-        print self.workItemNumber, self.gwi, self.lwi
-        dim = self.scalar.topology.dim
+        if __VERBOSE__:
+            print "Work-items basic config : ",
+            print self.workItemNumber, self.gwi, self.lwi
+        dim = self.advectedFields.topology.dim
         self.coord_min = np.ones(4, dtype=self.gpu_precision)
         self.mesh_size = np.ones(4, dtype=self.gpu_precision)
-        self.coord_min[0:dim] = np.asarray(self.scalar.topology.mesh.origin,
+        self.coord_min[0:dim] = np.asarray(self.advectedFields.topology.mesh.origin,
                                            dtype=self.gpu_precision)
         self.mesh_size[0:dim] = np.asarray(
-            self.scalar.topology.mesh.space_step,
+            self.advectedFields.topology.mesh.space_step,
             dtype=self.gpu_precision)
-        print "=== Kernel sources ==="
+        if __VERBOSE__:
+            print "=== Kernel sources ==="
         self.build_options = ""
         if _check_double_precision_capability(self.device,
                                               self.gpu_precision):
@@ -111,36 +116,37 @@ class GPUParticleAdvection(ParticleAdvection):
         self.build_options += " -D WGN=" + str(self.workItemNumber)
         self.gpu_src = self._collect_usr_cl_src(self.user_gpu_src)
         self.gpu_src += self._collect_kernels_cl_src()
-        print "===\n"
-        print "=== Kernel sources compiling ==="
+        if __VERBOSE__:
+            print "===\n"
+            print "=== Kernel sources compiling ==="
         ## Build code
         if self.gpu_precision is PARMES_REAL_GPU:
             prg = cl.Program(self.ctx, self.gpu_src.replace('.0', '.0f'))
         else:
             prg = cl.Program(self.ctx, self.gpu_src.replace('float', 'double'))
         ## OpenCL program
-        print self.build_options
         self.prg = prg.build(self.build_options)
-        print "Compiler options : ",
-        print self.prg.get_build_info(
-            self.device, cl.program_build_info.OPTIONS)
-        print "Compiler status : ",
-        print self.prg.get_build_info(
-            self.device, cl.program_build_info.STATUS)
-        print "Compiler log : ",
-        print self.prg.get_build_info(self.device, cl.program_build_info.LOG)
-        print "===\n"
+        if __VERBOSE__:
+            print self.build_options
+            print "Compiler options : ",
+            print self.prg.get_build_info(
+                self.device, cl.program_build_info.OPTIONS)
+            print "Compiler status : ",
+            print self.prg.get_build_info(
+                self.device, cl.program_build_info.STATUS)
+            print "Compiler log : ",
+            print self.prg.get_build_info(self.device,
+                                          cl.program_build_info.LOG)
+            print "===\n"
         print "=== OpenCL Buffer allocations ==="
         self._buffer_allocations()
         print "===\n"
         print "=== OpenCL Buffer initialisation ==="
         self._buffer_initialisations()
         print "===\n"
-        print "=== OpenCL Kernels affectation ==="
         self._kernel_affectations()
-        print "===\n"
-        self.numMethod = Splitting(self._apply_in_dir, self.scalar.dimension)
-        print "Method used:", self.name
+        self.numMethod = Splitting(self._apply_in_dir,
+                                   self.advectedFields.dimension)
 
     def _buffer_initialisations(self):
         data, transfer_time, compute_time = 0, 0., 0.
@@ -160,15 +166,18 @@ class GPUParticleAdvection(ParticleAdvection):
                             for d in xrange(len(self.resolution))]
                     args.append(self.coord_min)
                     args.append(self.mesh_size)
-                    evt = kl.launch(*tuple(args))
+                    evt = kl(*tuple(args))
                 else:
                     args = [gpudf.gpu_data]
                     args.append(self.coord_min)
                     args.append(self.mesh_size)
-                    evt = kl.launch(*tuple(args))
+                    evt = kl(*tuple(args))
                 self.queue.finish()
                 temp_time = (evt.profile.end - evt.profile.start) * 1e-9
-                print "Done in ", temp_time, "sec"
+                if __VERBOSE__:
+                    print "Done in ", temp_time, "sec"
+                else:
+                    print ".",
                 compute_time += temp_time
                 gpudf.contains_data = False
                 gpudf.data_on_device = True
@@ -178,6 +187,7 @@ class GPUParticleAdvection(ParticleAdvection):
                     temp_data, temp_time = gpudf.toDevice()
                     data += temp_data
                     transfer_time += temp_time
+        print ""
         if data > 0:
             print "Total Transfers : ", data, "Bytes transfered at ",
             print "{0:.3f} GBytes/sec".format((data * 1e-9) / transfer_time)
@@ -220,10 +230,11 @@ class GPUParticleAdvection(ParticleAdvection):
                                  int(self.resolution[1] / 16))
                 self.copy_lwi = (16, 2)
         f = open(GPU_SRC + src_copy)
-        print  "   - copy:", f.name
         gpu_src += "".join(f.readlines())
         f.close()
-        print "       - config :", self.copy_gwi, self.copy_lwi
+        if __VERBOSE__:
+            print  "   - copy:", GPU_SRC + src_copy
+            print "       - config :", self.copy_gwi, self.copy_lwi
         return gpu_src
 
     def _collect_kernels_cl_src_transpositions(self):
@@ -276,45 +287,50 @@ class GPUParticleAdvection(ParticleAdvection):
                                          int(self.resolution[1]) / 16)
                 self.transpose_xy_lwi = (8, 2)
         f = open(GPU_SRC + src_transpose_xy)
-        print  "   - transpose_xy:", f.name
         gpu_src += "".join(f.readlines())
         f.close()
-        print "       - config :", self.transpose_xy_gwi, self.transpose_xy_lwi
         if len(self.resolution) == 3:
-            print  "   - transpose_xz:", f.name
             f = open(GPU_SRC + src_transpose_xz)
             gpu_src += "".join(f.readlines())
             f.close()
-            print "       - config :",
-            print self.transpose_xz_gwi, self.transpose_xz_lwi
+        if __VERBOSE__:
+            print  "   - transpose_xy:", GPU_SRC + src_transpose_xy
+            print "       - config :", self.transpose_xy_gwi,
+            print self.transpose_xy_lwi
+            if len(self.resolution) == 3:
+                print  "   - transpose_xz:", GPU_SRC + src_transpose_xz
+                print "       - config :",
+                print self.transpose_xz_gwi, self.transpose_xz_lwi
         return gpu_src
 
     def _collect_usr_cl_src(self, usr_src):
-        print "usr src", usr_src
         gpu_src = ""
+        if __VERBOSE__:
+            print "Sources files (user): ", usr_src
         if usr_src is not None:
             if isinstance(usr_src, list):
                 for src in usr_src:
-                    print "Sources files (user): ", src
                     f = open(src, 'r')
                     gpu_src += "".join(f.readlines())
                     f.close()
             else:
-                print "Sources files (user): ", usr_src
                 f = open(usr_src, 'r')
                 gpu_src += "".join(f.readlines())
                 f.close()
         return gpu_src
 
+    @debug
     def apply(self, t, dt, ite):
         self.numMethod(t, dt)
 
+    @debug
     def finalize(self):
         """
         Free OpenCL memory.
         """
         for gpudf in self.variables:
-            print "deallocate :", gpudf.name
+            if __VERBOSE__:
+                print "deallocate :", gpudf.name
             if gpudf.vector:
                 [gpudf.gpu_data[d].release()
                  for d in xrange(len(self.resolution))]
diff --git a/HySoP/hysop/operator/gpu_particle_advection_1k.py b/HySoP/hysop/operator/gpu_particle_advection_1k.py
index b07179e3a467ed99e0356d03e97b6c2f3605cdf4..8f18a53857335c1e191a34c9842c8196d69f5b87 100644
--- a/HySoP/hysop/operator/gpu_particle_advection_1k.py
+++ b/HySoP/hysop/operator/gpu_particle_advection_1k.py
@@ -3,7 +3,8 @@
 
 Discrete advection representation
 """
-from parmepy.constants import np
+from parmepy.constants import np, debug
+from parmepy import __VERBOSE__
 from parmepy.operator.gpu_particle_advection import GPUParticleAdvection
 from parmepy.fields.continuous import Field
 from parmepy.fields.gpu_discrete import GPUVectorField
@@ -23,6 +24,7 @@ class GPUParticleAdvection1k(GPUParticleAdvection):
 
     """
 
+    @debug
     def __init__(self, advec,
                  platform_id=0, device_id=0,
                  device_type='gpu',
@@ -43,6 +45,7 @@ class GPUParticleAdvection1k(GPUParticleAdvection):
                                       method, src,
                                       precision)
 
+    @debug
     def setUp(self):
         GPUParticleAdvection.setUp(self)
 
@@ -77,22 +80,23 @@ class GPUParticleAdvection1k(GPUParticleAdvection):
         GPUVectorField.fromVectorField(self.queue, self.velocity,
                                        self.gpu_precision)
         ## Transported scalar.
-        GPUScalarField.fromScalarField(self.queue, self.scalar,
+        GPUScalarField.fromScalarField(self.queue, self.advectedFields,
                                        self.gpu_precision)
         ## Result scalar
-        particle_scalar = Field(self.scalar.topology.domain,
+        particle_scalar = Field(self.advectedFields.topology.domain,
                                 "Particle_Scalar",
                                 vector=False)
         self.part_scalar, dfId = particle_scalar.discretize(
-            self.scalar.topology)
+            self.advectedFields.topology)
         GPUScalarField.fromScalarField(
             self.queue, self.part_scalar,
             self.gpu_precision)
-        self.variables = [self.scalar, self.velocity,
+        self.variables = [self.advectedFields, self.velocity,
                           self.part_scalar]
         total_mem_used = 0
         for gpuv in self.variables:
             total_mem_used += gpuv.mem_size
+        print ""
         print "Total Global Memory used : ", total_mem_used,
         print "Bytes (", total_mem_used / (1024 ** 2), "MB)",
         print "({0:.3f} %)".format(
@@ -173,16 +177,18 @@ class GPUParticleAdvection1k(GPUParticleAdvection):
         self.build_options += str(self.WINb_advec_and_remesh)
         f = open(GPU_SRC + src_remesh_weights)
         gpu_src += "".join(f.readlines())
-        print  "   - remeshing weights:", f.name
         f.close()
         f = open(GPU_SRC + src_advec_and_remesh)
-        print  "   - advection and remeshing:", f.name
         gpu_src += "".join(f.readlines())
         f.close()
-        print "       - config :",
-        print self.gwi_advec_and_remesh, self.lwi_advec_and_remesh
         gpu_src += self._collect_kernels_cl_src_copy()
         gpu_src += self._collect_kernels_cl_src_transpositions()
+        if __VERBOSE__:
+            print  "   - remeshing weights:", GPU_SRC + src_remesh_weights
+            print  "   - advection and remeshing:",
+            print GPU_SRC + src_advec_and_remesh
+            print "       - config :",
+            print self.gwi_advec_and_remesh, self.lwi_advec_and_remesh
         return gpu_src
 
     def _apply_in_dir(self, t, dt, d, split_id):
@@ -205,29 +211,29 @@ class GPUParticleAdvection1k(GPUParticleAdvection):
         """
         # Particle init
         if split_id == 0:
-            evt_init = self.init_copy.launch(self.scalar.gpu_data,
-                                             self.part_scalar.gpu_data)
+            evt_init = self.init_copy(self.advectedFields.gpu_data,
+                                      self.part_scalar.gpu_data)
         else:
-            if self.scalar.topology.dim == 2:
+            if self.advectedFields.topology.dim == 2:
                 evt_init = \
-                    self.init_transpose.launch(self.scalar.gpu_data,
-                                               self.part_scalar.gpu_data)
+                    self.init_transpose(self.advectedFields.gpu_data,
+                                        self.part_scalar.gpu_data)
             else:
                 if min(split_id, d) == 0:
                     evt_init = \
-                        self.init_transpose.launch(0, self.scalar.gpu_data,
-                                                   self.part_scalar.gpu_data)
+                        self.init_transpose(0, self.advectedFields.gpu_data,
+                                            self.part_scalar.gpu_data)
                 else:
                     evt_init = \
-                        self.init_transpose.launch(1, self.scalar.gpu_data,
-                                                   self.part_scalar.gpu_data)
+                        self.init_transpose(1, self.advectedFields.gpu_data,
+                                            self.part_scalar.gpu_data)
         # Advection
-        evt_advec = self.num_advec_and_remesh.launch(self.velocity.gpu_data[d],
-                                                     self.part_scalar.gpu_data,
-                                                     self.scalar.gpu_data,
-                                                     self.gpu_precision(dt),
-                                                     self.coord_min[d],
-                                                     self.mesh_size[d])
+        evt_advec = self.num_advec_and_remesh(self.velocity.gpu_data[d],
+                                              self.part_scalar.gpu_data,
+                                              self.advectedFields.gpu_data,
+                                              self.gpu_precision(dt),
+                                              self.coord_min[d],
+                                              self.mesh_size[d])
         for df in self.output:
             df.data_on_device = True
             df.contains_data = False
@@ -239,7 +245,7 @@ class GPUParticleAdvection1k(GPUParticleAdvection):
         if split_id == 0:
             self.compute_time_copy[d] += c_time_init
         else:
-            if self.scalar.topology.dim == 2:
+            if self.advectedFields.topology.dim == 2:
                 self.compute_time_swap[0] += c_time_init
             else:
                 if min(split_id, d) == 0:
diff --git a/HySoP/hysop/operator/gpu_particle_advection_2k.py b/HySoP/hysop/operator/gpu_particle_advection_2k.py
index 698f0e468a6fe580fec7a75fa3bd636a28f2b95f..bd878cae7ad0f41c62f0d71db617e551aa0e6344 100644
--- a/HySoP/hysop/operator/gpu_particle_advection_2k.py
+++ b/HySoP/hysop/operator/gpu_particle_advection_2k.py
@@ -3,7 +3,8 @@
 
 Discrete advection representation
 """
-from parmepy.constants import np
+from parmepy.constants import np, debug
+from parmepy import __VERBOSE__
 from parmepy.operator.gpu_particle_advection import GPUParticleAdvection
 from parmepy.fields.continuous import Field
 from parmepy.fields.gpu_discrete import GPUVectorField
@@ -23,6 +24,7 @@ class GPUParticleAdvection2k(GPUParticleAdvection):
 
     """
 
+    @debug
     def __init__(self, advec,
                  platform_id=0, device_id=0,
                  device_type='gpu',
@@ -44,6 +46,7 @@ class GPUParticleAdvection2k(GPUParticleAdvection):
                                       precision)
         self.compute_time_remesh = [0., 0., 0.]
 
+    @debug
     def setUp(self):
         GPUParticleAdvection.setUp(self)
 
@@ -81,31 +84,32 @@ class GPUParticleAdvection2k(GPUParticleAdvection):
         GPUVectorField.fromVectorField(self.queue, self.velocity,
                                        self.gpu_precision)
         ## Transported scalar.
-        GPUScalarField.fromScalarField(self.queue, self.scalar,
+        GPUScalarField.fromScalarField(self.queue, self.advectedFields,
                                        self.gpu_precision)
         ## Particle position
-        particle_position = Field(self.scalar.topology.domain,
+        particle_position = Field(self.advectedFields.topology.domain,
                                   "Particle_Position",
                                   vector=False)
         self.part_position, dfID = particle_position.discretize(
-            self.scalar.topology)
+            self.advectedFields.topology)
         GPUScalarField.fromScalarField(
             self.queue, self.part_position,
             self.gpu_precision)
         ## Result scalar
-        particle_scalar = Field(self.scalar.topology.domain,
+        particle_scalar = Field(self.advectedFields.topology.domain,
                                 "Particle_Scalar",
                                 vector=False)
         self.part_scalar, dfId = particle_scalar.discretize(
-            self.scalar.topology)
+            self.advectedFields.topology)
         GPUScalarField.fromScalarField(
             self.queue, self.part_scalar,
             self.gpu_precision)
-        self.variables = [self.scalar, self.velocity,
+        self.variables = [self.advectedFields, self.velocity,
                           self.part_position, self.part_scalar]
         total_mem_used = 0
         for gpuv in self.variables:
             total_mem_used += gpuv.mem_size
+        print ""
         print "Total Global Memory used : ", total_mem_used,
         print "Bytes (", total_mem_used / (1024 ** 2), "MB)",
         print "({0:.3f} %)".format(
@@ -113,7 +117,6 @@ class GPUParticleAdvection2k(GPUParticleAdvection):
 
     def _collect_kernels_cl_src(self):
         gpu_src = ""
-        print self.method
         if len(self.resolution) == 3:
             if self.gpu_precision is PARMES_REAL_GPU:
                 src_advec = 'advection_3D_opt4_builtin.cl'
@@ -132,11 +135,8 @@ class GPUParticleAdvection2k(GPUParticleAdvection):
             self.lwi_advec = self._set_WorkItems(advec_vector_width)
         self.build_options += " -D WGNADVEC=" + str(self.WINb_advec)
         f = open(GPU_SRC + src_advec)
-        print  "   - advection:", f.name
         gpu_src += "".join(f.readlines())
         f.close()
-        print "       - config :",
-        print self.gwi_advec, self.lwi_advec
         ## remeshing
         remesh_vector_width = 4
         if len(self.resolution) == 3:
@@ -172,16 +172,20 @@ class GPUParticleAdvection2k(GPUParticleAdvection):
         self.build_options += " -D WGNREMESH=" + str(self.WINb_remesh)
         f = open(GPU_SRC + src_remesh_weights)
         gpu_src += "".join(f.readlines())
-        print  "   - remeshing weights:", f.name
         f.close()
         f = open(GPU_SRC + src_remesh)
-        print  "   - remeshing:", f.name
         gpu_src += "".join(f.readlines())
         f.close()
-        print "       - config :",
-        print self.gwi_remesh, self.lwi_remesh
         gpu_src += self._collect_kernels_cl_src_copy()
         gpu_src += self._collect_kernels_cl_src_transpositions()
+        if __VERBOSE__:
+            print  "   - advection:", GPU_SRC + src_advec
+            print "       - config :",
+            print self.gwi_advec, self.lwi_advec
+            print  "   - remeshing weights:", GPU_SRC + src_remesh_weights
+            print  "   - remeshing:", GPU_SRC + src_remesh
+            print "       - config :",
+            print self.gwi_remesh, self.lwi_remesh
         return gpu_src
 
 
@@ -205,35 +209,35 @@ class GPUParticleAdvection2k(GPUParticleAdvection):
         """
         # Particle init
         if split_id == 0:
-            evt_init = self.init_copy.launch(self.scalar.gpu_data,
-                                             self.part_scalar.gpu_data)
+            evt_init = self.init_copy(self.advectedFields.gpu_data,
+                                      self.part_scalar.gpu_data)
         else:
-            if self.scalar.topology.dim == 2:
+            if self.advectedFields.topology.dim == 2:
                 evt_init = \
-                    self.init_transpose.launch(self.scalar.gpu_data,
-                                               self.part_scalar.gpu_data)
+                    self.init_transpose(self.advectedFields.gpu_data,
+                                        self.part_scalar.gpu_data)
             else:
                 if min(split_id, d) == 0:
                     evt_init = \
-                        self.init_transpose.launch(0, self.scalar.gpu_data,
-                                                   self.part_scalar.gpu_data)
+                        self.init_transpose(0, self.advectedFields.gpu_data,
+                                            self.part_scalar.gpu_data)
                 else:
                     evt_init = \
-                        self.init_transpose.launch(1, self.scalar.gpu_data,
-                                                   self.part_scalar.gpu_data)
+                        self.init_transpose(1, self.advectedFields.gpu_data,
+                                            self.part_scalar.gpu_data)
         # Advection
-        evt_advec = self.num_advec.launch(self.velocity.gpu_data[d],
-                                          self.part_position.gpu_data,
-                                          self.gpu_precision(dt),
-                                          self.coord_min[d],
-                                          self.mesh_size[d])
+        evt_advec = self.num_advec(self.velocity.gpu_data[d],
+                                   self.part_position.gpu_data,
+                                   self.gpu_precision(dt),
+                                   self.coord_min[d],
+                                   self.mesh_size[d])
         self.num_advec.finish()
         # remeshing
-        evt_remesh = self.num_remesh.launch(self.part_position.gpu_data,
-                                            self.part_scalar.gpu_data,
-                                            self.scalar.gpu_data,
-                                            self.coord_min[d],
-                                            self.mesh_size[d])
+        evt_remesh = self.num_remesh(self.part_position.gpu_data,
+                                     self.part_scalar.gpu_data,
+                                     self.advectedFields.gpu_data,
+                                     self.coord_min[d],
+                                     self.mesh_size[d])
         for df in self.output:
             df.data_on_device = True
             df.contains_data = False
@@ -248,7 +252,7 @@ class GPUParticleAdvection2k(GPUParticleAdvection):
         if split_id == 0:
             self.compute_time_copy[d] += c_time_init
         else:
-            if self.scalar.topology.dim == 2:
+            if self.advectedFields.topology.dim == 2:
                 self.compute_time_swap[0] += c_time_init
             else:
                 if min(split_id, d) == 0:
diff --git a/HySoP/hysop/operator/multiphase.py b/HySoP/hysop/operator/multiphase.py
index 6b40bc1e7edf55389d62fcd58bdffb4b2824d3a1..0fbc9a078b7c4f7efc051dbaa1d39b9570862f77 100644
--- a/HySoP/hysop/operator/multiphase.py
+++ b/HySoP/hysop/operator/multiphase.py
@@ -14,6 +14,7 @@ class Pressure(Operator):
     Pressure operator representation
     """
 
+    @debug
     def __init__(self, velocity, vorticity, viscosity, density,
                  resolutions=None, method='',
                  **other_config):
@@ -33,6 +34,7 @@ class Pressure(Operator):
         self.method = method
         self.config = other_config
 
+    @debug
     def setUp(self):
         """
         Pressure operator discretization method.
diff --git a/HySoP/hysop/operator/multiphase_d.py b/HySoP/hysop/operator/multiphase_d.py
index 3abee35a49341785eb8992b64b9c39216f592490..9ef3b51f9ba8729b26ae2a3fbe7b3ae62aebb6bf 100644
--- a/HySoP/hysop/operator/multiphase_d.py
+++ b/HySoP/hysop/operator/multiphase_d.py
@@ -16,6 +16,7 @@ class Pressure_d(DiscreteOperator):
     Operator representation.
     """
 
+    @debug
     def __init__(self,pressureOp, method=''):
         """
         Constructor.
@@ -36,6 +37,7 @@ class Pressure_d(DiscreteOperator):
             pressureOp.resolutions[pressureOp.density])
         self.compute_time = 0.
 
+    @debug
     def setUp(self):
         ## La topologie se recupere depuis les variables discretes
         self.topology = self.velocity_old.topology
@@ -46,6 +48,7 @@ class Pressure_d(DiscreteOperator):
         self.resolution = self.topology.mesh.resolution
         self.meshSize = self.topology.mesh.size
 
+    @debug
     def apply(self, t, dt, ite):
         """
         Apply operator.
diff --git a/HySoP/hysop/operator/particle_advection.py b/HySoP/hysop/operator/particle_advection.py
index bb0da9e2f706e315bf6ea225fdebb0b79423e394..7a2366f95ad1004b2e4e4471b52a14fbf6c7bbc7 100644
--- a/HySoP/hysop/operator/particle_advection.py
+++ b/HySoP/hysop/operator/particle_advection.py
@@ -3,7 +3,7 @@
 
 Discrete advection representation
 """
-from parmepy.constants import np
+from parmepy.constants import np, debug
 from parmepy.operator.discrete import DiscreteOperator
 from parmepy.fields.continuous import Field
 from parmepy.gpu import PARMES_REAL_GPU
@@ -18,30 +18,29 @@ class ParticleAdvection(DiscreteOperator):
 
     """
 
+    @debug
     def __init__(self, advec, method=''):
         """
         Create a Advection operator.
-        Work on a given scalar at a given velocity to produce scalar
+        Work on a given advectedFields at a given velocity
+        to produce advectedFields
         distribution at new positions.
 
         @param advec : Transport operator
-        @param idVelocityD : Index of velocity discretisation to use.
-        @param idScalarD : Index of scalar discretisation to use.
-        @param result_position : result position.
-        @param result_scalar : result scalar.
         @param method : the method to use.
         """
         self.advec = advec
-        self.method = method
         ## Velocity.
         self.velocity = self.advec.velocity.discreteField[
             self.advec.discreteFieldId[self.advec.velocity]]
-        ## Transported scalar.
-        self.scalar = self.advec.advectedFields.discreteField[
+        ## Advected Fields.
+        self.advectedFields = self.advec.advectedFields.discreteField[
             self.advec.discreteFieldId[self.advec.advectedFields]]
-        DiscreteOperator.__init__(self, [self.velocity, self.scalar], method)
-        self.input = [self.velocity, self.scalar]
-        self.output = [self.scalar]
+        DiscreteOperator.__init__(self,
+                                  [self.velocity, self.advectedFields],
+                                  method)
+        self.input = [self.velocity, self.advectedFields]
+        self.output = [self.advectedFields]
         ## Compute time detailed per directions
         self.compute_time = [0., 0., 0.]
         ## Compute time for copy detailed per directions
@@ -52,19 +51,22 @@ class ParticleAdvection(DiscreteOperator):
         self.name = "advection_default"
         self.num_method = None
 
+    @debug
     def setUp(self):
         ## Result position
         particle_position = Field(self.advec.advectedFields.domain,
                                   "Particle_Position",
                                   vector=False)
-        particle_position.discretize(self.advec.resolutions[self.advec.advectedFields])
+        particle_position.discretize(
+            self.advec.resolutions[self.advec.advectedFields])
         self.part_position = particle_position.discreteField[0]
-        ## Result scalar
-        particle_scalar = Field(self.advec.advectedFields.domain,
-                                "Particle_Scalar",
-                                vector=False)
-        particle_scalar.discretize(self.advec.resolutions[self.advec.advectedFields])
-        self.part_scalar = particle_scalar.discreteField[0]
+        ## Result advectedFields
+        particle_advectedFields = Field(self.advec.advectedFields.domain,
+                                        "Particle_AdvectedFields",
+                                        vector=False)
+        particle_advectedFields.discretize(
+            self.advec.resolutions[self.advec.advectedFields])
+        self.part_advectedFields = particle_advectedFields.discreteField[0]
         if self.method.find('rk2') >= 0:
             self.name += '_rk2'
             from parmepy.numerics.integrators.runge_kutta2 import RK2
@@ -76,14 +78,15 @@ class ParticleAdvection(DiscreteOperator):
         if self.method.find('m6prime') >= 0:
             self.name += '_m6prime'
             from parmepy.numerics.remeshing.m6prime import M6Prime
-            self.num_remesh = M6Prime(self.scalar.dimension)
+            self.num_remesh = M6Prime(self.advectedFields.dimension)
         else:
             self.name += '_m6prime'
             from parmepy.numerics.remeshing.m6prime import M6Prime
-            self.num_remesh = M6Prime(self.scalar.dimension)
+            self.num_remesh = M6Prime(self.advectedFields.dimension)
         self.interpolation = [Linear(self.velocity[d])
-                              for d in xrange(self.scalar.dimension)]
-        self.numMethod = Splitting(self._apply_in_dir, self.scalar.dimension)
+                              for d in xrange(self.advectedFields.dimension)]
+        self.numMethod = Splitting(
+            self._apply_in_dir, self.advectedFields.dimension)
         print "Method used:", self.name
 
     def _apply_in_dir(self, t, dt, d, split_id):
@@ -96,23 +99,25 @@ class ParticleAdvection(DiscreteOperator):
 
         Advection algorithm:
         @li 1. Particle initialization : \n
-                 - by copy scalar from grid to particles if previous splitting
-                 direction equals current splitting direction.\n
-                 - by transposition of scalar from grid to particle.
+                 - by copy advectedFields from grid to particles
+                 if previous splitting direction equals current
+                 splitting direction.\n
+                 - by transposition of advectedFields from grid to particle.
         @li 2. Particle advection :\n
                  - compute particle position in splitting direction as a
-                 scalar. Performs a RK2 resolution of dx_p/dt = a_p.
+                 advectedFields. Performs a RK2 resolution of dx_p/dt = a_p.
         @li 3. Profile timings of OpenCL kernels.
         """
-        self.part_scalar[...] = self.scalar[...]
-        self.part_position[...] = self.scalar.topology.mesh.coords[d]
+        self.part_advectedFields[...] = self.advectedFields[...]
+        self.part_position[...] = self.advectedFields.topology.mesh.coords[d]
         self.part_position[...] = self.num_advec(self.interpolation[d],
                                                  t, dt,
                                                  self.part_position.data)
-        self.scalar[...] = self.num_remesh(self.part_position.data,
-                                           self.part_scalar.data,
-                                           d)
+        self.advectedFields[...] = \
+            self.num_remesh(self.part_position.data,
+                            self.part_advectedFields.data, d)
 
+    @debug
     def apply(self, t, dt, ite):
         self.numMethod(t, dt)
 
diff --git a/HySoP/hysop/operator/penalization.py b/HySoP/hysop/operator/penalization.py
index bb8c171b179a4437f368cb046ae7da28b2d4f75b..0e817874778b09a65690538e1ad9a28159830fed 100644
--- a/HySoP/hysop/operator/penalization.py
+++ b/HySoP/hysop/operator/penalization.py
@@ -14,6 +14,7 @@ class Penalization(Operator):
     Penalization operator representation
     """
 
+    @debug
     def __init__(self, velocity, vorticity, obstacle, lambd,
                  resolutions=None, method='',
                  **other_config):
@@ -36,6 +37,7 @@ class Penalization(Operator):
         self.method = method
         self.config = other_config
 
+    @debug
     def setUp(self):
         """
         Penalization operator discretization method.
diff --git a/HySoP/hysop/operator/penalization_d.py b/HySoP/hysop/operator/penalization_d.py
index e0b636de643aab12f39f228699c7ee6270814a27..60c162de50143d7863be37b5d843fa4e7bfb1f56 100644
--- a/HySoP/hysop/operator/penalization_d.py
+++ b/HySoP/hysop/operator/penalization_d.py
@@ -15,6 +15,7 @@ class Penalization_d(DiscreteOperator):
     DiscreteOperator.DiscreteOperator specialization.
     """
 
+    @debug
     def __init__(self, penal, method=''):
         """
         Constructor.
@@ -38,9 +39,11 @@ class Penalization_d(DiscreteOperator):
         ## Method
         self.method = method
 
+    @debug
     def setUp(self):
         self.obstacle.definedObst.setUp(self.velocity.topology)
 
+    @debug
     def apply(self, t, dt, ite):
 
         self.compute_time = time.time()
diff --git a/HySoP/hysop/operator/poisson.py b/HySoP/hysop/operator/poisson.py
index e96183ba85cc1b46ea4ff35bd6d5cc83d62bf803..e70f8d0755508568c6f0d63dc95cced1ec7be4f9 100644
--- a/HySoP/hysop/operator/poisson.py
+++ b/HySoP/hysop/operator/poisson.py
@@ -12,6 +12,7 @@ class Poisson(Operator):
     Poisson solver operator representation using FFT
     """
 
+    @debug
     def __init__(self, velocity=None, vorticity=None,
                  resolutions=None, method='',
                  **other_config):
@@ -29,6 +30,7 @@ class Poisson(Operator):
         self.method = method
         self.config = other_config
 
+    @debug
     def setUp(self):
         """
         Poisson solver operator discretization method.
diff --git a/HySoP/hysop/operator/poisson_fft.py b/HySoP/hysop/operator/poisson_fft.py
index f6d8353976f16f92c93da630db2075ef8e1f7072..866c3eb72a58a1b2a1c9843a26a181b651605af1 100644
--- a/HySoP/hysop/operator/poisson_fft.py
+++ b/HySoP/hysop/operator/poisson_fft.py
@@ -14,6 +14,7 @@ class PoissonFFT(DiscreteOperator):
     Operator representation.
     """
 
+    @debug
     def __init__(self, poisson, method=''):
         """
         Constructor.
@@ -29,6 +30,7 @@ class PoissonFFT(DiscreteOperator):
         self.compute_time = 0.
         self.method = method
 
+    @debug
     def setUp(self):
         #on recupere la topologie par les variables discretes
         self.topology = self.velocity.topology
@@ -38,6 +40,7 @@ class PoissonFFT(DiscreteOperator):
         self.dim = self.topology.dim
         self.coords = self.topology.mesh.coords
 
+    @debug
     def apply(self, t, dt, ite):
         """
         Apply operator.
diff --git a/HySoP/hysop/operator/scales_advection.py b/HySoP/hysop/operator/scales_advection.py
index cd78c0ce535ca9e5361048609389101a4b31c29a..1e02728467b9f0b670e41b131142debd3b5d9694 100644
--- a/HySoP/hysop/operator/scales_advection.py
+++ b/HySoP/hysop/operator/scales_advection.py
@@ -6,7 +6,7 @@ Discrete Advection operator using scales code (Jean-Baptiste)
 from parmepy.f2py import scales2py
 from parmepy.operator.discrete import DiscreteOperator
 from parmepy.fields.vector import VectorField
-from parmepy.constants import np, ORDER, PARMES_REAL
+from parmepy.constants import np, ORDER, PARMES_REAL, debug
 import time
 
 
@@ -15,6 +15,7 @@ class ScalesAdvection(DiscreteOperator):
     Operator representation.
     """
 
+    @debug
     def __init__(self, advec, method=''):
         """
         Constructor.
@@ -22,14 +23,14 @@ class ScalesAdvection(DiscreteOperator):
         self.advec = advec
         self.method = method
         self.velocity = self.advec.velocity.discreteField[
-            self.advec.dicreteFieldId[self.advec.velocity]]
+            self.advec.discreteFieldId[self.advec.velocity]]
         variables = [self.velocity]
         if not isinstance(self.advec.advectedFields, list):
             self.advectedFields = [self.advec.advectedFields.discreteField[
-                self.advec.dicreteFieldId[self.advec.advectedFields]]]
+                self.advec.discreteFieldId[self.advec.advectedFields]]]
         else:
             self.advectedFields = [
-                advecF.discreteField[self.advec.dicreteFieldId[advecF]]
+                advecF.discreteField[self.advec.discreteFieldId[advecF]]
                 for advecF in self.advec.advectedFields]
         [variables.append(advecF) for advecF in self.advectedFields]
         DiscreteOperator.__init__(self, variables,
@@ -38,13 +39,17 @@ class ScalesAdvection(DiscreteOperator):
         self.output = list(variables).remove(self.velocity)
         self.compute_time = 0.
 
+    @debug
     def setUp(self):
         #self.stab_coeff = stab_coeff
         self.topology = self.velocity.topology
         self.ghosts = self.topology.ghosts
         self.resolution = self.topology.mesh.resolution
         self.dim = self.topology.dim
+        # TODO scales setup
+        #scales2py.setup('p_O2', True, 'strang')
 
+    @debug
     def apply(self, t, dt, ite):
         """
         Apply operator.
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index 6337a152b2861f8498e98845c5d7e70e5ad12fda..928f11249075eeef24c0d12def682b702801d96a 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -16,6 +16,7 @@ class Stretching(Operator):
     Stretching operator representation
     """
 
+    @debug
     def __init__(self, velocity, vorticity,
                  resolutions=None, method='FD_order4 RK3',
                  **other_config):
@@ -36,6 +37,7 @@ class Stretching(Operator):
         self.method = method
         self.config = other_config
 
+    @debug
     def setUp(self):
         """
         Stretching operator discretization method.
@@ -43,10 +45,10 @@ class Stretching(Operator):
 
         @param idVelocityD : Index of velocity discretisation to use.
         @param idVorticityD : Index of vorticity discretisation to use.
-        @param propertyOp : property garantied by the stretching Op 
+        @param propertyOp : property garantied by the stretching Op
             (divConservation, massConservation) default : divConservation
         @param method : the method to use :
-            space : (FD_order2, FD_order4, ..) default : FD_order4 
+            space : (FD_order2, FD_order4, ..) default : FD_order4
             and time : (Euler, RK, ..) default : RK3
         """
         Operator.setUp(self)
diff --git a/HySoP/hysop/operator/stretching_d.py b/HySoP/hysop/operator/stretching_d.py
index ccdd63286a1509729c8b8edb3b780cbb6954b449..c243d1b5d311d6abdccdbca6a5551fe7b046d2c4 100755
--- a/HySoP/hysop/operator/stretching_d.py
+++ b/HySoP/hysop/operator/stretching_d.py
@@ -21,6 +21,7 @@ class Stretching_d(DiscreteOperator):
     DiscreteOperator.DiscreteOperator specialization.
     """
 
+    @debug
     def __init__(self, stretch, method='FD_order4 RK3', propertyOp='divConservation'):
         """
         Constructor.
@@ -44,6 +45,7 @@ class Stretching_d(DiscreteOperator):
         self.propertyOp = propertyOp
         self.method = method
 
+    @debug
     def setUp(self):
         self.cststretch = 0.
         if self.method.find('Euler') >= 0:
@@ -62,6 +64,7 @@ class Stretching_d(DiscreteOperator):
             self.num_method = RK3
             self.cststretch = 2.5127
 
+    @debug
     def apply(self, t, dt, ite):
         """
         Apply Stretching operator.
@@ -82,7 +85,7 @@ class Stretching_d(DiscreteOperator):
 
         if self.num_method is not None:
 
-            if (self.propertyOp.find('divConservation') < 0 
+            if (self.propertyOp.find('divConservation') < 0
                 and self.propertyOp.find('massConservation') < 0):
                 print 'No or wrong stretching operator property specified, use default divConservation'
                 self.propertyOp = 'divConservation'
@@ -117,7 +120,7 @@ class Stretching_d(DiscreteOperator):
 
                 ## fct2Op
                 self.function = Fct2Op(self.vorticity.data,
-                                       gradResult, 
+                                       gradResult,
                                        method='FD_order4',
                                        choice='gradV',
                                        topology=self.vorticity.topology)
diff --git a/HySoP/hysop/operator/velocity.py b/HySoP/hysop/operator/velocity.py
index 9fab39833f33328537a2261fabbbe208db3ca8f6..d1018a9f68a74c2da7af310e702ae5f220e145d6 100644
--- a/HySoP/hysop/operator/velocity.py
+++ b/HySoP/hysop/operator/velocity.py
@@ -13,6 +13,7 @@ class Velocity(Operator):
     Velocity operator representation.
     """
 
+    @debug
     def __init__(self, velocity, resolutions=None, method='',
                  **other_config):
         """
@@ -27,6 +28,7 @@ class Velocity(Operator):
         self.method = method
         self.config = other_config
 
+    @debug
     def setUp(self):
         """
         Advection operator discretization method.
diff --git a/HySoP/hysop/operator/velocity_d.py b/HySoP/hysop/operator/velocity_d.py
index 2d08222ffe827b2bcbe0e7023ee0833cf1667337..6e0943e68c7996150601880d0c08c1d3fd165c95 100644
--- a/HySoP/hysop/operator/velocity_d.py
+++ b/HySoP/hysop/operator/velocity_d.py
@@ -14,6 +14,7 @@ class Velocity_P(DiscreteOperator):
     Velocity computation from analytical expression.x
     """
 
+    @debug
     def __init__(self, velocity_Op, method=''):
         """
         Create a Velocity operator.
@@ -33,9 +34,11 @@ class Velocity_P(DiscreteOperator):
         self.compute_time = 0.
         self.name = "velocity"
 
+    @debug
     def setUp(self, method=''):
         raise NotImplemented("Not implemented")
 
+    @debug
     def apply(self, t, dt):
         c_time = 0.
         if self.numMethod is not None:
diff --git a/HySoP/hysop/problem/problem.py b/HySoP/hysop/problem/problem.py
index 69f0ac13b09c1d6f5989ff0571cb8b84ae66627c..961d298cc796e1c97c34ded2a54e937fe6180540 100644
--- a/HySoP/hysop/problem/problem.py
+++ b/HySoP/hysop/problem/problem.py
@@ -4,6 +4,8 @@
 Complete problem description.
 """
 from abc import ABCMeta, abstractmethod
+from parmepy.constants import debug
+from parmepy import __VERBOSE__
 from parmepy.operator.monitors.monitoring import Monitoring
 
 
@@ -21,6 +23,11 @@ class Problem(object):
 
     __metaclass__ = ABCMeta
 
+    @debug
+    def __new__(cls, *args, **kw):
+        return object.__new__(cls, *args, **kw)
+
+    @debug
     @abstractmethod
     def __init__(self, operators, monitors=[]):
         """
@@ -43,6 +50,7 @@ class Problem(object):
                     raise ValueError("Problem must have only one " +
                                      "domain for variables.")
 
+    @debug
     def setUp(self, t_end, dt):
         """
         Set solver.
@@ -51,10 +59,15 @@ class Problem(object):
         """
         [op.setUp() for op in self.operators]
         self.timer = Timer(t_end, dt)
+        if  __VERBOSE__:
+            print "==== Variables initialization ===="
         for op in self.operators:
             for v in op.variables:
                 v.initialize()
+        if __VERBOSE__:
+            print "===="
 
+    @debug
     def solve(self):
         """
         Solve problem.
diff --git a/HySoP/setup.py.in b/HySoP/setup.py.in
index c0de913d7206152db0dcdb7dc4c07278b1aa54fe..539f32ebc15edb3c8657808516d6436dc6ac3721 100644
--- a/HySoP/setup.py.in
+++ b/HySoP/setup.py.in
@@ -83,6 +83,14 @@ else:
     packages.append('parmepy.f2py.scales2py')
     packages.append('parmepy.f2py.fftw2py')
 
+data_files = []
+
+if("@WITH_GPU@" is "ON"):
+    data_files.append(('./parmepy/gpu/cl_src',
+                       ['@CMAKE_SOURCE_DIR@/parmepy/gpu/cl_src/' + cl_file
+                        for cl_file in os.listdir('@CMAKE_SOURCE_DIR@/parmepy/gpu/cl_src/')
+                        if cl_file[0]!='.' and cl_file[-3:]=='.cl']))
+
 config = Configuration(name=name,
       version='@PYPACKAGE_VERSION@',
       description='Particular Methods implementation in Python',
@@ -93,22 +101,7 @@ config = Configuration(name=name,
       package_dir={'': '@CMAKE_SOURCE_DIR@'},
       ext_modules=ext_modules,
       packages=packages,
-      data_files=[('./parmepy/gpu/cl_src',
-                   ['@CMAKE_SOURCE_DIR@/parmepy/gpu/cl_src/' + cl_file
-                    for cl_file in os.listdir('@CMAKE_SOURCE_DIR@/parmepy/gpu/cl_src/')
-                    if cl_file[0]!='.' and cl_file[-3:]=='.cl'])]
+      data_files=data_files,
 )
 
 setup(**config.todict())
-
-# In case of cmake WITH_TESTS=ON
-## if('@WITH_TESTS@' is "ON"):
-##     cTestFile = open('@CMAKE_CURRENT_BINARY_DIR@/CTestTestfile.cmake','r')
-##     cTestStr = ""
-##     for line in cTestFile:
-##         cTestStr += line.replace('BUILD_MAIN_PYTEST', config.get_build_temp_dir().replace('temp','lib')+'/parmepy/test/main_unit_tests.py')
-##     cTestFile.close()
-##     cTestFile = open('@CMAKE_CURRENT_BINARY_DIR@/CTestTestfile.cmake','w')
-##     cTestFile.write(cTestStr)
-##     cTestFile.close()
-