diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index a249b2d194182008eefb877953a1debaf63278f7..81a433f205d51325c74014c67ade7904426fe26a 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -39,6 +39,9 @@ class NavierStokes(bfps.code):
         self.parameters['nu'] = 0.1
         self.parameters['famplitude'] = 1.0
         self.parameters['fmode'] = 1
+        self.parameters['fk0'] = 0.0
+        self.parameters['fk1'] = 3.0
+        self.parameters['forcing_type'] = 'linear'
         self.parameters['nparticles'] = 0
         self.fluid_includes = '#include "fluid_solver.hpp"\n'
         self.fluid_variables = ''
@@ -112,6 +115,9 @@ class NavierStokes(bfps.code):
                 fs->nu = nu;
                 fs->fmode = fmode;
                 fs->famplitude = famplitude;
+                fs->fk0 = fk0;
+                fs->fk1 = fk1;
+                strncpy(fs->forcing_type, forcing_type, 128);
                 fs->iteration = iter0;
                 fs->read('v', 'c');
                 fs->low_pass_Fourier(fs->cvorticity, 3, fs->kM);
@@ -309,7 +315,7 @@ def test(opt):
     c.parameters['nu'] = 1e-1
     c.parameters['dt'] = 2e-3
     c.parameters['niter_todo'] = opt.nsteps
-    c.parameters['famplitude'] = 0.0
+    c.parameters['famplitude'] = 0.1
     c.parameters['nparticles'] = 32
     c.add_particles()
     c.add_particles(kcut = 'fs->kM/2')
@@ -324,7 +330,6 @@ def test(opt):
               simname = 'test')
     stats = c.read_stats()
     k, espec = c.read_spec()
-    print k
 
     # plot spectra
     fig = plt.figure(figsize=(6,6))
diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp
index 2cccd8431e0737099722dd2f11027d866f09ca11..522a7a6968c226b9459f7df955b7529a5a222a94 100644
--- a/src/fluid_solver.cpp
+++ b/src/fluid_solver.cpp
@@ -153,9 +153,13 @@ fluid_solver<R>::fluid_solver( \
             this->rv[2], this->cv[2], \
             MPI_COMM_WORLD, FFTW_ESTIMATE | FFTW_MPI_TRANSPOSED_OUT); \
  \
-    /* ``physical'' parameters etc */ \
+    /* ``physical'' parameters etc, initialized here just in case */ \
  \
     this->nu = 0.1; \
+    this->fmode = 1; \
+    this->famplitude = 1.0; \
+    this->fk0  = 0; \
+    this->fk1 = 3.0; \
 } \
  \
 template<> \
@@ -249,16 +253,39 @@ template<> \
 void fluid_solver<R>::add_forcing(\
         C *field, R factor) \
 { \
-    ptrdiff_t cindex; \
-    if (this->cd->myrank == this->cd->rank[this->fmode]) \
+    if (strcmp(this->forcing_type, "none") == 0) \
+        return; \
+    if (strcmp(this->forcing_type, "Kolmogorov") == 0) \
     { \
-        cindex = ((this->fmode - this->cd->starts[0]) * this->cd->sizes[1])*this->cd->sizes[2]*3; \
-        field[cindex+2][0] -= this->famplitude*factor/2; \
+        ptrdiff_t cindex; \
+        if (this->cd->myrank == this->cd->rank[this->fmode]) \
+        { \
+            cindex = ((this->fmode - this->cd->starts[0]) * this->cd->sizes[1])*this->cd->sizes[2]*3; \
+            field[cindex+2][0] -= this->famplitude*factor/2; \
+        } \
+        if (this->cd->myrank == this->cd->rank[this->cd->sizes[0] - this->fmode]) \
+        { \
+            cindex = ((this->cd->sizes[0] - this->fmode - this->cd->starts[0]) * this->cd->sizes[1])*this->cd->sizes[2]*3; \
+            field[cindex+2][0] -= this->famplitude*factor/2; \
+        } \
+        return; \
     } \
-    if (this->cd->myrank == this->cd->rank[this->cd->sizes[0] - this->fmode]) \
+    if (strcmp(this->forcing_type, "linear") == 0) \
     { \
-        cindex = ((this->cd->sizes[0] - this->fmode - this->cd->starts[0]) * this->cd->sizes[1])*this->cd->sizes[2]*3; \
-        field[cindex+2][0] -= this->famplitude*factor/2; \
+        double knorm; \
+        CLOOP( \
+                knorm = sqrt(this->kx[xindex]*this->kx[xindex] + \
+                             this->ky[yindex]*this->ky[yindex] + \
+                             this->kz[zindex]*this->kz[zindex]); \
+                if ((this->fk0  <= knorm) && \
+                    (this->fk1 >= knorm)) \
+                    for (int c=0; c<3; c++) \
+                    { \
+                        field[cindex*3+c][0] += this->famplitude*this->cvorticity[cindex*3+c][0]*factor; \
+                        field[cindex*3+c][1] += this->famplitude*this->cvorticity[cindex*3+c][1]*factor; \
+                    } \
+             ); \
+        return; \
     } \
 } \
  \
diff --git a/src/fluid_solver.hpp b/src/fluid_solver.hpp
index bb7843a57468e3a41f93a701d4cc330a608bb0fd..0e0fdb14bcdf81b14211b3d718c8581caf4939fe 100644
--- a/src/fluid_solver.hpp
+++ b/src/fluid_solver.hpp
@@ -60,8 +60,10 @@ class fluid_solver:public fluid_solver_base<rnumber>
 
         /* physical parameters */
         double nu;
-        int fmode;
-        double famplitude;
+        int fmode;         // for Kolmogorov flow
+        double famplitude; // both for Kflow and band forcing
+        double fk0, fk1;   // for band forcing
+        char forcing_type[128];
 
         /* methods */
         fluid_solver(
diff --git a/test.py b/test.py
index b26e05e5355c11a584413cc9c3c29992fac9d843..60538494d84f2124351110f0d6cb66faf2fbed21 100755
--- a/test.py
+++ b/test.py
@@ -270,7 +270,7 @@ if __name__ == '__main__':
     parser.add_argument('--run', dest = 'run', action = 'store_true')
     parser.add_argument('--ncpu', type = int, dest = 'ncpu', default = 2)
     parser.add_argument('--nsteps', type = int, dest = 'nsteps', default = 16)
-    parser.add_argument('-n', type = int, dest = 'n', default = 32)
+    parser.add_argument('-n', type = int, dest = 'n', default = 64)
     opt = parser.parse_args()
     NStest(opt)