diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl
index f576c8b..7951873 100644
--- a/src/LatticeGPU.jl
+++ b/src/LatticeGPU.jl
@@ -31,6 +31,6 @@ include("YM/YM.jl")
 using .YM
 export YMworkspace, GaugeParm, force0_wilson!, field, field_pln, randomize!, zero!, norm2
 export gauge_action, hamiltonian, plaquette, HMC!, OMF4!
-export wfl_euler, wfl_rk3
+export wfl_euler, wfl_rk3, zfl_euler, zfl_rk3
 
 end # module
diff --git a/src/YM/YM.jl b/src/YM/YM.jl
index b5672e5..dcc4aee 100644
--- a/src/YM/YM.jl
+++ b/src/YM/YM.jl
@@ -94,6 +94,6 @@ include("YMhmc.jl")
 export gauge_action, hamiltonian, plaquette, HMC!, OMF4!
 
 include("YMflow.jl")
-export wfl_euler, wfl_rk3
+export wfl_euler, wfl_rk3, zfl_euler, zfl_rk3
 
 end
diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl
index 75105a2..83a7e03 100644
--- a/src/YM/YMact.jl
+++ b/src/YM/YMact.jl
@@ -70,6 +70,7 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, lp::SpaceParm{N,M,D}) wher
                 gc = U[b2,id2,r2]
             end
             h3 = (ga*gb)/gc
+            # END staples
             
             if ru2 == r
                 gb = Ush[bu2,1]
diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl
index ba4e517..932f676 100644
--- a/src/YM/YMflow.jl
+++ b/src/YM/YMflow.jl
@@ -42,3 +42,36 @@ function wfl_rk3(U, ns, eps, lp::SpaceParm, ymws::YMworkspace)
     return nothing
 end
 
+function zfl_euler(U, ns, eps, lp::SpaceParm, ymws::YMworkspace)
+
+    for i in 1:ns
+        force_gauge(ymws, U, 5.0/3.0, lp)
+        U .= expm.(U, ymws.frc1, 2*eps)
+    end
+    
+    return nothing
+end
+
+function zfl_rk3(U, ns, eps, lp::SpaceParm, ymws::YMworkspace)
+
+    for i in 1:ns
+        c0 = eps/2
+        force_gauge(ymws, U, 5.0/3.0, lp)
+        ymws.mom .= ymws.frc1
+        U .= expm.(U, ymws.mom, c0)
+
+        c0 = -34*eps/36
+        c1 = 16*eps/9
+        force_gauge(ymws, U, 5.0/3.0, lp)
+        ymws.mom .= c0.*ymws.mom .+ c1.*ymws.frc1
+        U .= expm.(U, ymws.mom)
+
+        c1 = 6*eps/4
+        force_gauge(ymws, U, 5.0/3.0, lp)
+        ymws.mom .= c1.*ymws.frc1 .- ymws.mom 
+        U .= expm.(U, ymws.mom)
+    end
+
+    return nothing
+end
+
diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl
index 0a3f4ca..55bfbbc 100644
--- a/src/YM/YMhmc.jl
+++ b/src/YM/YMhmc.jl
@@ -13,7 +13,7 @@
     
     function gauge_action(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace)
 
-Returns the value of the Wilson plaquette action for the configuration U. 
+Returns the value of the gauge plaquette action for the configuration U. The parameters `\beta` and `c0` are taken from the `gp` structure. 
 """
 function gauge_action(U, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where T <: AbstractFloat
 
diff --git a/src/main/test.jl b/src/main/test.jl
index 1268753..9ac1b9f 100644
--- a/src/main/test.jl
+++ b/src/main/test.jl
@@ -12,7 +12,7 @@ println(CUDA.device())
 GRP  = SU3
 ALG  = SU3alg
 PREC = Float64
-lp = SpaceParm{4}((32,32,32,32), (4,4,4,4))
+lp = SpaceParm{4}((64,64,64,64), (4,4,4,4))
 gp = GaugeParm{PREC}(6.0, 5.0/3.0, (0.0,0.0), 3)
 
 println("Space  Parameters: ", lp)