diff --git a/Examples/NS_planeJet.py b/Examples/NS_planeJet.py
new file mode 100644
index 0000000000000000000000000000000000000000..a050b048cf593fa2d4ceafd2ba49ebb6c7a26fff
--- /dev/null
+++ b/Examples/NS_planeJet.py
@@ -0,0 +1,264 @@
+#!/usr/bin/python
+
+"""
+Taylor Green 3D : see paper van Rees 2011.
+
+All parameters are set and defined in python module dataTG.
+
+"""
+
+import parmepy as pp
+from parmepy.constants import ORDER, np, PARMES_REAL
+from parmepy.f2py import fftw2py
+import numpy as np
+from parmepy.fields.continuous import Field
+from parmepy.variables.variables import Variables
+from parmepy.mpi.topology import Cartesian
+from parmepy.operator.advection import Advection
+from parmepy.operator.stretching import Stretching
+from parmepy.operator.poisson import Poisson
+from parmepy.operator.diffusion import Diffusion
+from parmepy.operator.adapt_timestep import AdaptTimeStep
+from parmepy.operator.redistribute import Redistribute
+from parmepy.problem.navier_stokes import NSProblem
+from parmepy.operator.monitors.printer import Printer
+from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
+from parmepy.problem.simulation import Simulation
+from parmepy.operator.adapt_timestep import AdaptTimeStep
+from parmepy.methods_keys import Scales
+from parmepy.operator.differential import Curl
+from parmepy.numerics.updateGhosts import UpdateGhosts
+from parmepy.mpi.main_var import main_size
+from parmepy.methods_keys import Scales, TimeIntegrator, Interpolation,\
+ Remesh, Support, Splitting, dtAdvecCrit, SpaceDiscretisation, GhostUpdate
+from parmepy.numerics.integrators.runge_kutta2 import RK2 as RK2
+from parmepy.numerics.integrators.runge_kutta3 import RK3 as RK3
+from parmepy.numerics.integrators.runge_kutta4 import RK4 as RK4
+from parmepy.numerics.finite_differences import FD_C_4, FD_C_2
+from parmepy.numerics.interpolation import Linear
+from parmepy.numerics.remeshing import L4_2 as rmsh
+from parmepy.f2py import fftw2py
+
+# problem dimension
+dim = 3
+# resolution
+nb = 129
+# number of ghosts in usual cartesian topo
+NBGHOSTS = 2
+ADVECTION_METHOD = {Scales: 'p_66'}
+# ADVECTION_METHOD = {TimeIntegrator: RK2,
+# Interpolation: Linear,
+# Remesh: rmsh,
+# Support: '',
+# Splitting: 'o2_FullHalf'}
+VISCOSITY = 5e-6
+
+
+## ----------- A 3d problem -----------
+print " ========= Start Navier-Stokes 3D (Taylor Green benchmark) ========="
+## pi constant
+pi = np.pi
+cos = np.cos
+sin = np.sin
+exp = np.exp
+abs = np.abs
+tanh = np.tanh
+## Domain
+box = pp.Box(dim, length=[1., ]*3)
+
+## Global resolution
+nbElem = [nb] * dim
+randX = np.asarray(np.random.random([nb - 1, nb - 1,(nb - 1)/main_size]),
+ dtype=PARMES_REAL, order=ORDER) - 0.5
+randY = np.asarray(np.random.random([nb - 1, nb - 1,(nb - 1)/main_size]),
+ dtype=PARMES_REAL, order=ORDER) - 0.5
+randZ = np.asarray(np.random.random([nb - 1, nb - 1,(nb - 1)/main_size]),
+ dtype=PARMES_REAL, order=ORDER) - 0.5
+
+# ##initjet.f
+# width = 0.01
+# ampl3 = 0.3
+# ampl2 = 0.
+# ampl = 0.05
+# def computeVel(res, x, y, z, t):
+# yy = abs(y - 0.5)
+# aux=(0.1-2.*yy)/(4.*width)
+# res[0][...] =0.5*(1.+tanh(aux))*(1.+ampl3*sin(8.*pi*(x)))*(1.+ampl*exp(-abs((0.1-2.*yy)/(4.*width)**2))*randX)
+# res[1][...] = ampl*exp(-abs((0.1-2.*yy)/(4.*width)**2))*randY
+# res[2][...] = ampl*exp(-abs((0.1-2.*yy)/(4.*width)**2))*randZ
+# return res
+
+
+# def initScal(res, x, y, z, t):
+# yy = abs(y - 0.5)
+# aux=(0.1-2.*yy)/(4.*width)
+# res[0][...] = 0.5*(1.+tanh(aux))*(1.+ampl2*exp(-abs((0.1-2.*yy)/(4.*width)**2))*randZ)
+# return res
+
+
+# ## JCP initial condition
+def computeVel(res, x, y, z, t):
+ res[0][...] = 0.5*(1.-tanh((abs(y-0.5)-0.1/2.)/0.02))*(1.+0.3*sin(8*pi*x))
+ res[1][...] = 0.
+ res[2][...] = 0.
+ return res
+
+
+def initScal(res, x, y, z, t):
+ res[0][...] = 0.5*(1.-tanh((abs(y-0.5)-0.1/2.)/0.02))*(1.+0.3*sin(8*pi*x))
+ return res
+
+
+## Fields
+velo = Field(domain=box, formula=computeVel,
+ name='Velocity', isVector=True)
+vorti = Field(domain=box,
+ name='Vorticity', isVector=True)
+scal = Field(domain=box, formula=initScal,
+ name='PassiveScalar', isVector=False)
+
+## Variables
+dt_adapt = Variables(domain=box, name='adaptative time step',
+ data=[0.01])
+
+## Variables
+## Usual Cartesian topology definition
+# At the moment we use two (or three?) topologies :
+# - "topo" for Stretching and all operators based on finite differences.
+# --> ghost layer = 2
+# - topo from Advection operator for all the other operators.
+# --> no ghost layer
+# - topo from fftw for Poisson and Diffusion.
+# Todo : check compat between scales and fft operators topologies.
+ghosts = [NBGHOSTS] * box.dimension
+topo = Cartesian(box, box.dimension, nbElem,
+ ghosts=ghosts)
+
+## Navier Stokes Operators
+advec = Advection(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ method=ADVECTION_METHOD
+ )
+advecScal = Advection(velo, scal,
+ resolutions={velo: nbElem,
+ scal: nbElem},
+ method=ADVECTION_METHOD
+ )
+
+stretch = Stretching(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ topo=topo
+ )
+
+diffusion = Diffusion(vorti,
+ resolution=nbElem,
+ viscosity=VISCOSITY
+ )
+
+poisson = Poisson(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ projection=[True, 1, False],
+ )
+
+c = Curl(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ method={SpaceDiscretisation: fftw2py, GhostUpdate: False},
+ )
+
+## Tools Operators
+dtAdapt = AdaptTimeStep(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ dt_adapt=dt_adapt,
+ method={TimeIntegrator: RK3,
+ SpaceDiscretisation: FD_C_4,
+ dtAdvecCrit: 'deform'},
+ topo=topo,
+ lcfl=0.125,
+ cfl=None, # 0.5,
+ prefix='./res_PS/dt_2p')
+
+## Simulation
+simu = Simulation(tinit=0.0,
+ tend=4.5, timeStep=dt_adapt,
+ iterMax=10000)
+
+
+
+# Bridges between the different topologies in order to
+# redistribute data.
+# 1 -Advection to stretching
+toGhost_vorti = Redistribute([vorti], advec, stretch)
+toGhost_velo = Redistribute([velo], poisson, stretch)
+# 2 - Stretching to Poisson/Diffusion
+fromGhost_vorti = Redistribute([vorti], stretch, diffusion)
+# 3 - Poisson to TimeStep
+distrPoissTimeStep = Redistribute([velo, vorti], poisson, dtAdapt)
+
+
+# Define the problem to solve
+pb = NSProblem(operators=[toGhost_velo, advecScal, advec,
+ toGhost_vorti,
+ stretch,
+ fromGhost_vorti,
+ diffusion,
+ poisson,
+ distrPoissTimeStep, dtAdapt],
+ simulation=simu, dumpFreq=-1)
+
+## Setting solver to Problem (only operators for computational tasks)
+pb.pre_setUp()
+## Diagnostics related to the problem
+topofft = poisson.discreteFields[poisson.vorticity].topology
+# energy = Energy_enstrophy(velo, vorti, isNormalized=True,
+# topo=topofft,
+# viscosity=VISCOSITY,
+# frequency=1,
+# prefix='./res_PS/energy_S128_2p.dat')
+
+# initialisation
+vorti.setTopoInit(topofft)
+velo.setTopoInit(topofft)
+scal.setTopoInit(topofft)
+#pb.addMonitors([energy])
+pb.setUp()
+
+c.discretize()
+c.setUp()
+c.apply(simu)
+
+# p = Printer(variables=[scal,],
+# topo=topofft,
+# frequency=1,
+# prefix='./res_PS/scal_S128_2p',
+# ext=".vtk")
+# pb.addMonitors([p])
+# p._printStep()
+
+# pf = Printer(variables=[velo, vorti],
+# topo=topofft,
+# frequency=1,
+# prefix='./res_PS/fields_S128_2p',
+# ext=".vtk")
+# pb.addMonitors([pf])
+# pf._printStep()
+
+
+def run():
+ pb.solve()
+
+## Solve problem
+from parmepy.mpi import MPI
+print "Start computation ..."
+time = MPI.Wtime()
+run()
+print 'total time (rank):', MPI.Wtime() - time, '(', topo.rank, ')'
+
+
+pb.finalize()
+## Clean memory buffers
+#fftw2py.clean_fftw_solver(box.dimension)