diff --git a/ci/scripts/test.sh b/ci/scripts/test.sh
index 4b8827696792e96cae03812f304ddee3cf19b776..a1581e72bc67219aa96c4f54538046c0f4d48757 100755
--- a/ci/scripts/test.sh
+++ b/ci/scripts/test.sh
@@ -68,8 +68,8 @@ python "$HYSOP_DIR/operator/tests/test_diffusion.py"
 python "$HYSOP_DIR/operator/tests/test_directional_stretching.py"
 
 # If scales (fortran advection library) is installed
-python -c "from hysop.f2hysop import scales2py as scales" && python "$HYSOP_DIR/operator/tests/test_bilevel_advection.py"
 python -c "from hysop.f2hysop import scales2py as scales" && python "$HYSOP_DIR/operator/tests/test_scales_advection.py"
+python -c "from hysop.f2hysop import scales2py as scales" && python "$HYSOP_DIR/operator/tests/test_bilevel_advection.py"
 
 EXAMPLE_DIR="$HYSOP_DIR/../examples"
 EXAMPLE_OPTIONS='-maxit 1'
diff --git a/examples/taylor_green/taylor_green.py b/examples/taylor_green/taylor_green.py
index f00d48c77dafd2f87bdbc70cc4979dceee9a6260..0adc5ae8683473842f684e3acb9cd5c6e7ed5531 100644
--- a/examples/taylor_green/taylor_green.py
+++ b/examples/taylor_green/taylor_green.py
@@ -89,7 +89,7 @@ def compute(args):
     
     ### Build the directional operators
     if (impl is Implementation.FORTRAN):
-        advec = Advection(implementation=Implementation.FORTRAN,
+        advec = Advection(implementation=impl,
                 name='advec',
                 velocity = velo,       
                 advected_fields = (vorti,),
diff --git a/hysop/operator/base/advection_dir.py b/hysop/operator/base/advection_dir.py
index 71480c78ffa06656628749752463b9232c0c836c..f56f585fc817201fab25765084837fc142c3415a 100644
--- a/hysop/operator/base/advection_dir.py
+++ b/hysop/operator/base/advection_dir.py
@@ -185,20 +185,25 @@ class DirectionalAdvectionBase(object):
         velocity_cfl           = self.velocity_cfl
         v_topo, v_requirements = requirements.get_input_requirement(velocity)
         v_dx                   = v_topo.space_step
-
+        
         dim  = velocity.domain.dim
         ndir = dim-1-direction
         
-        advection_ghosts = self.velocity_advection_ghosts(velocity_cfl)
-        min_velocity_ghosts = npw.integer_zeros(v_dx.shape)
-        min_velocity_ghosts[ndir] = advection_ghosts
-        v_requirements.min_ghosts = npw.maximum(v_requirements.min_ghosts, min_velocity_ghosts)
-        
         scalar = self.first_scalar
         s_topo, _ = requirements.get_input_requirement(scalar)
         s_dx = s_topo.space_step
         scalar_cfl = velocity_cfl * (v_dx[ndir] / s_dx[ndir])
+
+        advection_ghosts = self.velocity_advection_ghosts(velocity_cfl)
+        min_velocity_ghosts = npw.integer_zeros(v_dx.shape)
+        has_bilevel = any(v_dx != s_dx)
+        if has_bilevel:
+            # add one ghost in each direction for nd interpolation
+            min_velocity_ghosts[...] = 1
+        min_velocity_ghosts[ndir] = advection_ghosts
+        v_requirements.min_ghosts = npw.maximum(v_requirements.min_ghosts, min_velocity_ghosts)
         
+
         scalar_advection_ghosts = self.velocity_advection_ghosts(scalar_cfl)
         remesh_ghosts = self.scalars_out_cache_ghosts(scalar_cfl, self.remesh_kernel)
         assert remesh_ghosts > scalar_advection_ghosts
@@ -256,7 +261,7 @@ class DirectionalAdvectionBase(object):
         self.dvelocity            = dvelocity
         self.dadvected_fields_in  = dadvected_fields_in
         self.dadvected_fields_out = dadvected_fields_out
-
+        
         self.is_bilevel = None
         if any(self.dvelocity.compute_resolution != \
                self.dadvected_fields_out.values()[0].compute_resolution):