diff --git a/cpp/full_code/ornstein_uhlenbeck_process.cpp b/cpp/full_code/ornstein_uhlenbeck_process.cpp
index b6af42d324a27cddedd56aefbf7c57632ce186c6..93edcadb6ee96f3945998461c498c6deedf33403 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.cpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.cpp
@@ -16,6 +16,7 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
     double ou_kmin,
     double ou_kmax,
     double ou_energy_amplitude,
+    double ou_gamma,
     double DKX,
     double DKY,
     double DKZ,
@@ -44,6 +45,7 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
     this->ou_kmin_squ = pow(ou_kmin,2);
     this->ou_kmax_squ = pow(ou_kmax,2);
     this->ou_energy_amplitude = ou_energy_amplitude;
+    this->ou_gamma = ou_gamma;
     this->epsilon = pow((this->ou_energy_amplitude/this->kolmogorov_constant), 3./2.);
 
     assert(this->kk->kM2 >= this->ou_kmax_squ); 
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.hpp b/cpp/full_code/ornstein_uhlenbeck_process.hpp
index a1aeaa98dd948ca39b2c96177341890b345c9953..6c6bfb33e4a8b560a55d36e170d8de276e500ba2 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.hpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.hpp
@@ -19,7 +19,8 @@ class ornstein_uhlenbeck_process{
         double ou_kmin_squ;
         double ou_kmax_squ;
         double ou_energy_amplitude;
-        double kolmogorov_constant = 2;
+        double ou_gamma;
+        const double kolmogorov_constant = 2;
         double epsilon;
 
         field<rnumber,be,THREE> *ou_field;
@@ -38,6 +39,7 @@ class ornstein_uhlenbeck_process{
             double ou_kmin,
             double ou_kmax,
             double ou_energy_amplitude,
+            double ou_gamma,
             double DKX = 1.0,
             double DKY = 1.0,
             double DKZ = 1.0,
@@ -52,8 +54,13 @@ class ornstein_uhlenbeck_process{
 
         inline double gamma(double kabs)
         {
-            return this->kolmogorov_constant*sqrt(this->kolmogorov_constant/
-                            this->ou_energy_amplitude);
+            //return this->kolmogorov_constant*sqrt(this->kolmogorov_constant/
+            //                this->ou_energy_amplitude);
+            //return pow(this->ou_energy_amplitude/this->kolmogorov_constant,-0.5)
+            //        *pow(kabs,-2.0/3.0);
+            
+            return this->ou_gamma*pow(kabs,-2/3);
+                  
         }
         
 
diff --git a/cpp/full_code/ou_vorticity_equation.hpp b/cpp/full_code/ou_vorticity_equation.hpp
index e36292c2271616794815d43d084810b9daeb3fc9..9e4d329c6adad32ad5eaf04d23663a21dd3a3426 100644
--- a/cpp/full_code/ou_vorticity_equation.hpp
+++ b/cpp/full_code/ou_vorticity_equation.hpp
@@ -24,6 +24,7 @@ class ou_vorticity_equation : public vorticity_equation<rnumber, be>
         double ou_kmin,
         double ou_kmax,
         double ou_energy_amplitude,
+        double ou_gamma,
         double DKX = 1.0,
         double DKY = 1.0,
         double DKZ = 1.0,
@@ -34,7 +35,7 @@ class ou_vorticity_equation : public vorticity_equation<rnumber, be>
               this->ou = new ornstein_uhlenbeck_process<rnumber,be>(
                 NAME,
                 nx, ny, nz,
-                ou_kmin, ou_kmax, ou_energy_amplitude,
+                ou_kmin, ou_kmax, ou_energy_amplitude, ou_gamma,
                 DKX, DKY, DKZ, FFTW_PLAN_RIGOR);
               this->ou_forcing_type = "replace";
             }