diff --git a/CMakeLists.txt b/CMakeLists.txt
index 349201e..37225d5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -49,7 +49,7 @@ message(STATUS "CMAKE_CXX_COMPILER: " ${CMAKE_CXX_COMPILER_ID})
 
 if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
     if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 6.0) # Because of GCC bug 48891
-        message(FATAL_ERROR "GCC version 6.0 or higher required")    
+        message(FATAL_ERROR "GCC version 6.0 or higher required")
     endif()
 endif()
 
@@ -128,19 +128,14 @@ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}\
 set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG}\
                              -O0 -g")
 
-set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror -Wno-error=unused-parameter\
-                        -Wno-error=unused-function -Wno-error=unknown-pragmas")
-
-# Also warn about implicit conversions if the compiler supports it
-if (${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU")
-    set (CXX_FLAGS_WARNING "${CXX_FLAGS_WARNING} -Wdouble-promotion -Wfloat-conversion")
+if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
+    set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror")# -Wdouble-promotion -Wfloat-conversion")
+elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
+    set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror-all -Wsign-conversion")
+else()
+    message(WARNING "Using an unknown compiler. Compilation warning flags were not set.")
 endif()
 
-# Other flags. -D_FORCE_INLINES is a workaround to some CUDA/C++ "feature"
-# which botches the compilation ("memcpy was not declared in this scope")
-# (Not required with cc >= 3.0)
-#set(CXX_FLAGS_ETC "-D_FORCE_INLINES")
-
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}\
                     ${CXX_FLAGS_WARNING}\
                     ${CXX_FLAGS_ETC}\
diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps
index c1a2214..eb0fc0c 100644
--- a/acc/mhd_solver/stencil_process.sps
+++ b/acc/mhd_solver/stencil_process.sps
@@ -12,6 +12,9 @@ uniform Scalar dsx;
 uniform Scalar dsy;
 uniform Scalar dsz;
 
+uniform Scalar lnT0;
+uniform Scalar lnrho0;
+
 uniform int nx_min;
 uniform int ny_min;
 uniform int nz_min;
@@ -58,9 +61,10 @@ continuity(in Vector uu, in Scalar lnrho) {
 Vector
 momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa) {
     const Matrix S = stress_tensor(uu);
-    const Scalar cs2 = cs2_sound * exp(gamma * value(ss) / cp_sound + (gamma - 1) * (value(lnrho) - LNRHO0));
+    const Scalar cs2 = cs2_sound * exp(gamma * value(ss) / cp_sound + (gamma - 1) * (value(lnrho) - lnrho0));
     const Vector  j = (Scalar(1.) / mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
     const Vector B = curl(aa);
+    //TODO: DOES INTHERMAL VERSTION INCLUDE THE MAGNETIC FIELD?
     const Scalar inv_rho = Scalar(1.) / exp(value(lnrho));
 
     // Regex replace CPU constants with get\(AC_([a-zA-Z_0-9]*)\)
@@ -100,21 +104,23 @@ momentum(in Vector uu, in Scalar lnrho, in Scalar tt) {
 #else
 Vector
 momentum(in Vector uu, in Scalar lnrho) {
-    // !!!!!!!!!!!!!!!!%JP: NOTE TODO IMPORTANT!!!!!!!!!!!!!!!!!!!!!!!!
-    // NOT CHECKED FOR CORRECTNESS: USE AT YOUR OWN RISK
-    const Matrix S = stress_tensor(uu);
-    const Scalar cs2 = cs2_sound * exp((gamma - 1) * (value(lnrho) - LNRHO0));
-    // Regex replace CPU constants with get\(AC_([a-zA-Z_0-9]*)\)
-    // \1
-    const Vector mom = - mul(gradients(uu), value(uu))
-                                                       - cs2 * gradient(lnrho)
-                                                       + nu_visc * (
-                                                            laplace_vec(uu)
-                                                        + Scalar(1. / 3.) * gradient_of_divergence(uu)
-                                                        + Scalar(2.) * mul(S, gradient(lnrho))
-                                                        )
-                                                        + zeta * gradient_of_divergence(uu);
-    return mom;
+  Vector mom;
+
+  const Matrix S = stress_tensor(uu);
+
+    // Isothermal: we have constant speed of sound
+
+  mom = -mul(gradients(uu), value(uu)) -
+    cs2_sound * gradient(lnrho) +
+    nu_visc *
+    (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
+      Scalar(2.) * mul(S, gradient(lnrho))) + zeta * gradient_of_divergence(uu);
+
+  #if LGRAVITY
+  mom = mom - (Vector){0, 0, -10.0};
+  #endif
+
+  return mom;
 }
 #endif
 
@@ -139,8 +145,8 @@ induction(in Vector uu, in Vector aa) {
 #if LENTROPY
 Scalar
 lnT( in Scalar ss, in Scalar lnrho) {
-  const Scalar lnT = LNT0 + gamma * value(ss) / cp_sound +
-    (gamma - Scalar(1.)) * (value(lnrho) - LNRHO0);
+  const Scalar lnT = lnT0 + gamma * value(ss) / cp_sound +
+    (gamma - Scalar(1.)) * (value(lnrho) - lnrho0);
   return lnT;
 }
 
@@ -193,8 +199,12 @@ heat_transfer(in Vector uu, in Scalar lnrho, in Scalar tt)
 }
 #endif
 
+
+
 #if LFORCING
 
+
+
 Vector
 simple_vortex_forcing(Vector a, Vector b, Scalar magnitude)
 {
@@ -207,24 +217,68 @@ simple_outward_flow_forcing(Vector a, Vector b, Scalar magnitude)
     return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
 }
 
+
+// The Pencil Code forcing_hel_noshear(), manual Eq. 222, inspired forcing function with adjustable helicity
 Vector
-forcing(int3 globalVertexIdx)
+helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi)
+{
+
+    xx.x = xx.x*(2.0*M_PI/(dsx*(DCONST_INT(AC_ny_max) - DCONST_INT(AC_ny_min))));
+    xx.y = xx.y*(2.0*M_PI/(dsy*(DCONST_INT(AC_ny_max) - DCONST_INT(AC_ny_min))));
+    xx.z = xx.z*(2.0*M_PI/(dsz*(DCONST_INT(AC_ny_max) - DCONST_INT(AC_ny_min))));
+
+    Scalar cos_phi = cos(phi);
+    Scalar sin_phi = sin(phi);
+    Scalar cos_k_dox_x     = cos(dot(k_force, xx));
+    Scalar sin_k_dox_x     = sin(dot(k_force, xx));
+    // Phase affect only the x-component
+    //Scalar real_comp       = cos_k_dox_x;
+    //Scalar imag_comp       = sin_k_dox_x;
+    Scalar real_comp_phase = cos_k_dox_x*cos_phi - sin_k_dox_x*sin_phi;
+    Scalar imag_comp_phase = cos_k_dox_x*sin_phi + sin_k_dox_x*cos_phi;
+
+
+    Vector force = (Vector){ ff_re.x*real_comp_phase - ff_im.x*imag_comp_phase,
+                             ff_re.y*real_comp_phase - ff_im.y*imag_comp_phase,
+                             ff_re.z*real_comp_phase - ff_im.z*imag_comp_phase};
+
+    return force;
+}
+
+Vector
+forcing(int3 globalVertexIdx, Scalar dt)
 {
     Vector a = Scalar(.5) * (Vector){globalGrid.n.x * dsx,
                                      globalGrid.n.y * dsy,
                                      globalGrid.n.z * dsz}; // source (origin)
-    Vector b = (Vector){(globalVertexIdx.x - nx_min) * dsx,
+    Vector xx = (Vector){(globalVertexIdx.x - nx_min) * dsx,
                         (globalVertexIdx.y - ny_min) * dsy,
                         (globalVertexIdx.z - nz_min) * dsz}; // sink (current index)
+    const Scalar cs2 = cs2_sound;
+    const Scalar cs = sqrt(cs2);
+
+    //Placeholders until determined properly
+    Scalar magnitude = DCONST_REAL(AC_forcing_magnitude);
+    Scalar phase     = DCONST_REAL(AC_forcing_phase);
+    Vector k_force   = (Vector){  DCONST_REAL(AC_k_forcex),   DCONST_REAL(AC_k_forcey),   DCONST_REAL(AC_k_forcez)};
+    Vector ff_re     = (Vector){DCONST_REAL(AC_ff_hel_rex), DCONST_REAL(AC_ff_hel_rey), DCONST_REAL(AC_ff_hel_rez)};
+    Vector ff_im     = (Vector){DCONST_REAL(AC_ff_hel_imx), DCONST_REAL(AC_ff_hel_imy), DCONST_REAL(AC_ff_hel_imz)};
 
-    Scalar magnitude = 0.05;
 
     //Determine that forcing funtion type at this point.
-    Vector c = simple_vortex_forcing(a, b, magnitude);
-    //Vector c = simple_outward_flow_forcing(a, b, magnitude);
+    //Vector force = simple_vortex_forcing(a, xx, magnitude);
+    //Vector force = simple_outward_flow_forcing(a, xx, magnitude);
+    Vector force   = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase);
 
-    if (is_valid(c)) { return c; }
-    else             { return (Vector){0, 0, 0}; }
+    //Scaling N = magnitude*cs*sqrt(k*cs/dt)  * dt
+    const Scalar NN = cs*sqrt(DCONST_REAL(AC_kaver)*cs);
+    //MV: Like in the Pencil Code. I don't understandf the logic here.
+    force.x = sqrt(dt)*NN*force.x;
+    force.y = sqrt(dt)*NN*force.y;
+    force.z = sqrt(dt)*NN*force.z;
+
+    if (is_valid(force)) { return force; }
+    else                 { return (Vector){0, 0, 0}; }
 }
 #endif // LFORCING
 
@@ -272,7 +326,7 @@ solve(Scalar dt) {
 
     #if LFORCING
     if (step_number == 2) {
-        out_uu = out_uu + dt * forcing(globalVertexIdx);
+        out_uu = out_uu + forcing(globalVertexIdx, dt);
     }
     #endif
 }
diff --git a/analysis/python/jupyter/notebook_example.ipynb b/analysis/python/jupyter/notebook_example.ipynb
index e85aef1..f7e4956 100644
--- a/analysis/python/jupyter/notebook_example.ipynb
+++ b/analysis/python/jupyter/notebook_example.ipynb
@@ -24,7 +24,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "meshdir  = \"/scratch/data/mvaisala/forcingtest/\""
+    "meshdir  = \"/scratch/data/mvaisala/forcingtest/\"\n",
+    "#meshdir  = \"/scratch/data/mvaisala/normaltest/\""
    ]
   },
   {
@@ -34,7 +35,7 @@
    "outputs": [],
    "source": [
     "#imesh = 30000\n",
-    "imesh = 30000\n",
+    "imesh = 100\n",
     "mesh = ad.read.Mesh(imesh, fdir=meshdir)"
    ]
   },
@@ -101,15 +102,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "vis.lineplot.plot_ts(ts, lnrho=1, uutot=1, uux=1, uuy=1, ss=1)"
+    "vis.lineplot.plot_ts(ts, lnrho=1, uutot=1, uux=1, uuy=1, uuz=1, ss=1)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/config/astaroth.conf b/config/astaroth.conf
index 6d83764..41b7e51 100644
--- a/config/astaroth.conf
+++ b/config/astaroth.conf
@@ -38,6 +38,10 @@ AC_chi = 0.0001
 
 // Forcing
 AC_relhel = 0.0
+AC_forcing_magnitude = 1e-5
+AC_kmin              = 0.8
+AC_kmax              = 1.2 
+
 
 // Entropy
 AC_cp_sound = 1.0
diff --git a/include/astaroth.h b/include/astaroth.h
index dbe9b17..c354f1e 100644
--- a/include/astaroth.h
+++ b/include/astaroth.h
@@ -136,7 +136,6 @@ extern "C" {
         FUNC(AC_cs_sound), \
         FUNC(AC_eta), \
         FUNC(AC_mu0), \
-        FUNC(AC_relhel), \
         FUNC(AC_cp_sound), \
         FUNC(AC_gamma), \
         FUNC(AC_cv_sound), \
@@ -152,6 +151,23 @@ extern "C" {
         FUNC(AC_angl_uu), \
         FUNC(AC_lnrho_edge),\
         FUNC(AC_lnrho_out),\
+	/* Forcing parameters. User configured. */\
+        FUNC(AC_forcing_magnitude),\
+        FUNC(AC_relhel), \
+        FUNC(AC_kmin), \
+        FUNC(AC_kmax), \
+	/* Forcing parameters. Set by the generator. */\
+        FUNC(AC_forcing_phase),\
+        FUNC(AC_k_forcex),\
+        FUNC(AC_k_forcey),\
+        FUNC(AC_k_forcez),\
+        FUNC(AC_kaver),\
+        FUNC(AC_ff_hel_rex),\
+        FUNC(AC_ff_hel_rey),\
+        FUNC(AC_ff_hel_rez),\
+        FUNC(AC_ff_hel_imx),\
+        FUNC(AC_ff_hel_imy),\
+        FUNC(AC_ff_hel_imz),\
         /* Additional helper params */\
         /* (deduced from other params do not set these directly!) */\
         FUNC(AC_G_CONST),\
@@ -386,6 +402,15 @@ AcResult acQuit(void);
     unless otherwise stated. */
 AcResult acSynchronize(void);
 
+/** */
+AcResult acLoadDeviceConstant(const AcRealParam param, const AcReal value);
+
+/** Tool for loading forcing vector information into the device memory
+ */
+AcResult acForcingVec(const AcReal forcing_magnitude, const AcReal3 k_force,
+                      const AcReal3 ff_hel_re, const AcReal3 ff_hel_im, const AcReal forcing_phase,
+                      const AcReal kaver);
+
 /* End extern "C" */
 #ifdef __cplusplus
 }
diff --git a/src/core/astaroth.cu b/src/core/astaroth.cu
index 9c37117..0c49477 100644
--- a/src/core/astaroth.cu
+++ b/src/core/astaroth.cu
@@ -505,3 +505,44 @@ acSynchronize(void)
 
     return AC_SUCCESS;
 }
+
+AcResult
+acLoadDeviceConstant(const AcRealParam param, const AcReal value)
+{
+    for (int i = 0; i < num_devices; ++i) {
+        loadDeviceConstant(devices[i], param, value);
+    }
+    return AC_SUCCESS;
+}
+
+// Tool for loading forcing vector information into the device memory
+// %JP: Added a generic function for loading device constants (acLoadDeviceConstant).
+// This acForcingVec should go outside the core library since it references user-defined
+// parameters such as AC_forcing_magnitude which may not be defined in all projects.
+// host_forcing.cc is probably a good place for this.
+AcResult
+acForcingVec(const AcReal forcing_magnitude, const AcReal3 k_force, const AcReal3 ff_hel_re,
+             const AcReal3 ff_hel_im, const AcReal forcing_phase, const AcReal kaver)
+{
+
+    for (int i = 0; i < num_devices; ++i) {
+        loadDeviceConstant(devices[i], AC_forcing_magnitude, forcing_magnitude);
+        loadDeviceConstant(devices[i], AC_forcing_phase, forcing_phase);
+
+        loadDeviceConstant(devices[i], AC_k_forcex, k_force.x);
+        loadDeviceConstant(devices[i], AC_k_forcey, k_force.y);
+        loadDeviceConstant(devices[i], AC_k_forcez, k_force.z);
+
+        loadDeviceConstant(devices[i], AC_ff_hel_rex, ff_hel_re.x);
+        loadDeviceConstant(devices[i], AC_ff_hel_rey, ff_hel_re.y);
+        loadDeviceConstant(devices[i], AC_ff_hel_rez, ff_hel_re.z);
+
+        loadDeviceConstant(devices[i], AC_ff_hel_imx, ff_hel_im.x);
+        loadDeviceConstant(devices[i], AC_ff_hel_imy, ff_hel_im.y);
+        loadDeviceConstant(devices[i], AC_ff_hel_imz, ff_hel_im.z);
+
+        loadDeviceConstant(devices[i], AC_kaver, kaver);
+    }
+
+    return AC_SUCCESS;
+}
diff --git a/src/core/kernels/kernels.cuh b/src/core/kernels/kernels.cuh
index 1ff9051..a8d4d50 100644
--- a/src/core/kernels/kernels.cuh
+++ b/src/core/kernels/kernels.cuh
@@ -543,39 +543,6 @@ normalized(const AcReal3& vec)
     return inv_len * vec;
 }
 
-// Sinusoidal forcing
-// https://arxiv.org/pdf/1704.04676.pdf
-// NOTE: This method of forcing is depracated. However, it will remain in here
-// until a corresponding scheme exists in the new code.
-__constant__ AcReal3 forcing_vec;
-__constant__ AcReal forcing_phi;
-static __device__ __forceinline__ AcReal3
-DEPRECATED_forcing(const int i, const int j, const int k)
-{
-#define DOMAIN_SIZE_X (DCONST_INT(AC_nx) * DCONST_REAL(AC_dsx))
-#define DOMAIN_SIZE_Y (DCONST_INT(AC_ny) * DCONST_REAL(AC_dsy))
-#define DOMAIN_SIZE_Z (DCONST_INT(AC_nz) * DCONST_REAL(AC_dsz))
-    const AcReal3 k_vec = (AcReal3){
-        (i - DCONST_INT(AC_nx_min)) * DCONST_REAL(AC_dsx) - AcReal(.5) * DOMAIN_SIZE_X,
-        (j - DCONST_INT(AC_ny_min)) * DCONST_REAL(AC_dsy) - AcReal(.5) * DOMAIN_SIZE_Y,
-        (k - DCONST_INT(AC_nz_min)) * DCONST_REAL(AC_dsz) - AcReal(.5) * DOMAIN_SIZE_Z};
-    AcReal inv_len = reciprocal_len(k_vec);
-    if (isnan(inv_len) || isinf(inv_len))
-        inv_len = 0;
-    if (inv_len > 2) // hack to make it cool
-        inv_len = 2;
-    const AcReal k_dot_x = dot(k_vec, forcing_vec);
-
-    const AcReal waves = cos(k_dot_x) * cos(forcing_phi) - sin(k_dot_x) * sin(forcing_phi);
-
-    return inv_len * inv_len * waves * forcing_vec;
-}
-
-// Note: LNT0 and LNRHO0 must be set very carefully: if the magnitude is different that other values
-// in the mesh, then we will inherently lose precision
-#define LNT0 (AcReal(0.0))
-#define LNRHO0 (AcReal(0.0))
-
 #define H_CONST (AcReal(0.0))
 #define C_CONST (AcReal(0.0))
 
@@ -739,31 +706,11 @@ read_out(const int idx, AcReal* __restrict__ field[], const int3 handle)
  * =============================================================================
  */
 
-static AcReal
-randf(void)
-{
-    return AcReal(rand()) / AcReal(RAND_MAX);
-}
-
 AcResult
 rk3_step_async(const cudaStream_t stream, const int& step_number, const int3& start,
                const int3& end, const AcReal dt, VertexBufferArray* buffer)
 {
     const dim3 tpb(32, 1, 4);
-/////////////////// Forcing
-#if LFORCING
-    const AcReal ff_scale = AcReal(.2);
-    static AcReal3 ff     = ff_scale * (AcReal3){1, 0, 0};
-    const AcReal radians  = randf() * AcReal(2 * M_PI) / 360 / 8;
-    const AcMatrix rotz   = create_rotz(radians);
-    ff                    = mul(rotz, ff);
-    cudaMemcpyToSymbolAsync(forcing_vec, &ff, sizeof(ff), 0, cudaMemcpyHostToDevice, stream);
-
-    const AcReal ff_phi = AcReal(M_PI); // AcReal(2 * M_PI) * randf();
-    cudaMemcpyToSymbolAsync(forcing_phi, &ff_phi, sizeof(ff_phi), 0, cudaMemcpyHostToDevice,
-                            stream);
-#endif // LFORCING
-    //////////////////////////
 
     const int nx = end.x - start.x;
     const int ny = end.y - start.y;
diff --git a/src/standalone/autotest.cc b/src/standalone/autotest.cc
index 107b950..0f961fa 100644
--- a/src/standalone/autotest.cc
+++ b/src/standalone/autotest.cc
@@ -1,5 +1,5 @@
 /*
-   Copyright (C) 2014-2018, Johannes Pekkilae, Miikka Vaeisalae.
+   Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
 
    This file is part of Astaroth.
 
@@ -32,6 +32,7 @@
 #include "core/math_utils.h"
 #include "model/host_memory.h"
 #include "model/host_timestep.h"
+#include "model/host_forcing.h"
 #include "model/model_boundconds.h"
 #include "model/model_reduce.h"
 #include "model/model_rk3.h"
@@ -699,6 +700,13 @@ run_autotest(void)
                                                              VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
                 const AcReal dt   = host_timestep(umax, config);
 
+#if LFORCING
+  
+                //CURRENTLY AUTOTEST NOT SUPPORTED WITH FORCING!!!
+
+#endif
+
+
                 // Host integration step
                 model_rk3(dt, model_mesh);
                 boundconds(config, model_mesh);
diff --git a/src/standalone/model/host_forcing.cc b/src/standalone/model/host_forcing.cc
new file mode 100644
index 0000000..4ce7892
--- /dev/null
+++ b/src/standalone/model/host_forcing.cc
@@ -0,0 +1,177 @@
+/*
+    Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
+
+    This file is part of Astaroth.
+
+    Astaroth is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    Astaroth is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Astaroth.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+/**
+ * @file
+ * \brief Brief info.
+ *
+ * Detailed info.
+ *
+ */
+#include "host_forcing.h"
+
+#include "core/math_utils.h"
+
+// The is a wrapper for genering random numbers with a chosen system.
+AcReal
+get_random_number_01()
+{
+    // TODO: Implement better randon number generator http://www.cplusplus.com/reference/random/
+    return AcReal(rand()) / AcReal(RAND_MAX);
+}
+
+AcReal3
+cross(const AcReal3& a, const AcReal3& b)
+{
+    AcReal3 c;
+
+    c.x = a.y * b.z - a.z * b.y;
+    c.y = a.z * b.x - a.x * b.z;
+    c.z = a.x * b.y - a.y * b.x;
+
+    return c;
+}
+
+AcReal
+dot(const AcReal3& a, const AcReal3& b)
+{
+    return a.x * b.x + a.y * b.y + a.z * b.z;
+}
+
+AcReal3
+vec_norm(const AcReal3& a)
+{
+    AcReal3 c;
+    AcReal norm = dot(a, a);
+
+    c.x = a.x / sqrt(norm);
+    c.y = a.y / sqrt(norm);
+    c.z = a.z / sqrt(norm);
+
+    return c;
+}
+
+AcReal3
+vec_multi_scal(const AcReal scal, const AcReal3& a)
+{
+    AcReal3 c;
+
+    c.x = a.x * scal;
+    c.y = a.y * scal;
+    c.z = a.z * scal;
+
+    return c;
+}
+
+// Generate forcing wave vector k_force
+AcReal3
+helical_forcing_k_generator(const AcReal kmax, const AcReal kmin)
+{
+    AcReal phi, theta, kk; // Spherical direction coordinates
+    AcReal3 k_force;       // forcing wave vector
+
+    AcReal delta_k = kmax - kmin;
+
+    // Generate vector in spherical coordinates
+    phi   = AcReal(2.0) * AcReal(M_PI) * get_random_number_01();
+    theta = AcReal(M_PI) * get_random_number_01();
+    kk    = delta_k * get_random_number_01() + kmin;
+
+    // Cast into Cartesian form
+    k_force = (AcReal3){AcReal(kk * sin(theta) * cos(phi)), //
+                        AcReal(kk * sin(theta) * sin(phi)), //
+                        AcReal(kk * cos(theta))};
+
+    // printf("k_force.x %f, k_force.y %f, k_force.z %f \n", k_force.x, k_force.y, k_force.z);
+
+    // Round the numbers. In that way k(x/y/z) will get complete waves.
+    k_force.x = round(k_force.x);
+    k_force.y = round(k_force.y);
+    k_force.z = round(k_force.z);
+
+    // printf("After rounding --> k_force.x %f, k_force.y %f, k_force.z %f \n", k_force.x,
+    // k_force.y, k_force.z);
+
+    return k_force;
+}
+
+// Generate the unit perpendicular unit vector e required for helical forcing
+// Addapted from Pencil code forcing.f90 hel_vec() subroutine.
+void
+helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force)
+{
+
+    AcReal3 k_cross_e         = cross(k_force, *e_force);
+    k_cross_e                 = vec_norm(k_cross_e);
+    AcReal3 k_cross_k_cross_e = cross(k_force, k_cross_e);
+    k_cross_k_cross_e         = vec_norm(k_cross_k_cross_e);
+    AcReal phi                = AcReal(2.0) * AcReal(M_PI) * get_random_number_01();
+    AcReal3 ee_tmp1           = vec_multi_scal(cos(phi), k_cross_e);
+    AcReal3 ee_tmp2           = vec_multi_scal(sin(phi), k_cross_k_cross_e);
+
+    *e_force = (AcReal3){ee_tmp1.x + ee_tmp2.x, ee_tmp1.y + ee_tmp2.y, ee_tmp1.z + ee_tmp2.z};
+}
+
+// PC Manual Eq. 223
+void
+helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcReal3 k_force,
+                               const AcReal3 e_force, const AcReal relhel)
+{
+
+    // k dot e
+    AcReal3 kdote;
+    kdote.x = k_force.x * e_force.x;
+    kdote.y = k_force.y * e_force.y;
+    kdote.z = k_force.z * e_force.z;
+
+    // k cross e
+    AcReal3 k_cross_e;
+    k_cross_e.x = k_force.y * e_force.z - k_force.z * e_force.y;
+    k_cross_e.y = k_force.z * e_force.x - k_force.x * e_force.z;
+    k_cross_e.z = k_force.x * e_force.y - k_force.y * e_force.x;
+
+    // k cross k cross e
+    AcReal3 k_cross_k_cross_e;
+    k_cross_k_cross_e.x = k_force.y * k_cross_e.z - k_force.z * k_cross_e.y;
+    k_cross_k_cross_e.y = k_force.z * k_cross_e.x - k_force.x * k_cross_e.z;
+    k_cross_k_cross_e.z = k_force.x * k_cross_e.y - k_force.y * k_cross_e.x;
+
+    // abs(k)
+    AcReal kabs = sqrt(k_force.x * k_force.x + k_force.y * k_force.y + k_force.z * k_force.z);
+
+    AcReal denominator = sqrt(AcReal(1.0) + relhel * relhel) * kabs *
+                         sqrt(kabs * kabs -
+                              (kdote.x * kdote.x + kdote.y * kdote.y + kdote.z * kdote.z));
+
+    // MV: I suspect there is a typo in the Pencil Code manual!
+    //*ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator,
+    //                       -relhel*kabs*k_cross_e.y/denominator,
+    //                       -relhel*kabs*k_cross_e.z/denominator};
+
+    //*ff_hel_im = (AcReal3){k_cross_k_cross_e.x/denominator,
+    //                       k_cross_k_cross_e.y/denominator,
+    //                       k_cross_k_cross_e.z/denominator};
+
+    // See PC forcing.f90 forcing_hel_both()
+    *ff_hel_re = (AcReal3){kabs * k_cross_e.x / denominator, kabs * k_cross_e.y,
+                           kabs * k_cross_e.z};
+
+    *ff_hel_im = (AcReal3){relhel * k_cross_k_cross_e.x / denominator, relhel * k_cross_k_cross_e.y,
+                           relhel * k_cross_k_cross_e.z};
+}
diff --git a/src/standalone/model/host_forcing.h b/src/standalone/model/host_forcing.h
new file mode 100644
index 0000000..569485e
--- /dev/null
+++ b/src/standalone/model/host_forcing.h
@@ -0,0 +1,45 @@
+
+/*
+    Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
+
+    This file is part of Astaroth.
+
+    Astaroth is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    Astaroth is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Astaroth.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+/**
+ * @file
+ * \brief Brief info.
+ *
+ * Detailed info.
+ *
+ */
+#pragma once
+#include "astaroth.h"
+
+AcReal get_random_number_01();
+
+AcReal3 cross(const AcReal3& a, const AcReal3& b);
+
+AcReal dot(const AcReal3& a, const AcReal3& b);
+
+AcReal3 vec_norm(const AcReal3& a);
+
+AcReal3 vec_multi_scal(const AcReal scal, const AcReal3& a);
+
+AcReal3 helical_forcing_k_generator(const AcReal kmax, const AcReal kmin);
+
+void helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force);
+
+void helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcReal3 k_force, const AcReal3 e_force, const AcReal relhel);
diff --git a/src/standalone/model/host_timestep.cc b/src/standalone/model/host_timestep.cc
index 9c51d83..54f47cc 100644
--- a/src/standalone/model/host_timestep.cc
+++ b/src/standalone/model/host_timestep.cc
@@ -50,7 +50,8 @@ host_timestep(const AcReal& umax, const AcMeshInfo& mesh_info)
     // New, closer to the actual Courant timestep
     // See Pencil Code user manual p. 38 (timestep section)
     const long double uu_dt   = cdt * dsmin / (fabsl(umax) + sqrtl(cs2_sound + 0.0l));
-    const long double visc_dt = cdtv * dsmin * dsmin / max(max(nu_visc, eta), max(gamma, chi)) + 1; // TODO NOTE: comment the +1 out to get scientifically accurate results
+    const long double visc_dt = cdtv * dsmin * dsmin / max(max(nu_visc, eta), max(gamma, chi));// + 1; // TODO NOTE: comment the +1 out to get scientifically accurate results
+    //MV: White the +1? It was messing up my computations! 
 
     const long double dt = min(uu_dt, visc_dt);
     return AcReal(timescale) * AcReal(dt);
diff --git a/src/standalone/model/model_rk3.cc b/src/standalone/model/model_rk3.cc
index a18f1b8..dbb4413 100644
--- a/src/standalone/model/model_rk3.cc
+++ b/src/standalone/model/model_rk3.cc
@@ -1,5 +1,5 @@
 /*
-    Copyright (C) 2014-2018, Johannes Pekkilae, Miikka Vaeisalae.
+    Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
 
     This file is part of Astaroth.
 
@@ -553,11 +553,6 @@ normalized(const ModelVector& vec)
     return inv_len * vec;
 }
 
-// Note: LNT0 and LNRHO0 must be set very carefully: if the magnitude is different that other values
-// in the mesh, then we will inherently lose precision
-#define LNT0 (ModelScalar(0.0))
-#define LNRHO0 (ModelScalar(0.0))
-
 #define H_CONST (ModelScalar(0.0))
 #define C_CONST (ModelScalar(0.0))
 
@@ -571,9 +566,10 @@ momentum(const ModelVectorData& uu, const ModelScalarData& lnrho
 {
 #if LENTROPY
     const ModelMatrix S   = stress_tensor(uu);
-    const ModelScalar cs2 = get(AC_cs2_sound) * expl(get(AC_gamma) * value(ss) / get(AC_cp_sound) +
-                                                     (get(AC_gamma) - 1) * (value(lnrho) - LNRHO0));
-    const ModelVector j   = (ModelScalar(1.) / get(AC_mu0)) *
+    const ModelScalar cs2 = get(AC_cs2_sound) *
+                            expl(get(AC_gamma) * value(ss) / get(AC_cp_sound) +
+                                 (get(AC_gamma) - 1) * (value(lnrho) - get(AC_lnrho0)));
+    const ModelVector j = (ModelScalar(1.) / get(AC_mu0)) *
                           (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
     const ModelVector B       = curl(aa);
     const ModelScalar inv_rho = ModelScalar(1.) / expl(value(lnrho));
@@ -591,7 +587,8 @@ momentum(const ModelVectorData& uu, const ModelScalarData& lnrho
     // !!!!!!!!!!!!!!!!%JP: NOTE TODO IMPORTANT!!!!!!!!!!!!!!!!!!!!!!!!
     // NOT CHECKED FOR CORRECTNESS: USE AT YOUR OWN RISK
     const ModelMatrix S   = stress_tensor(uu);
-    const ModelScalar cs2 = get(AC_cs2_sound) * expl((get(AC_gamma) - 1) * (value(lnrho) - LNRHO0));
+    const ModelScalar cs2 = get(AC_cs2_sound) *
+                            expl((get(AC_gamma) - 1) * (value(lnrho) - get(AC_lnrho0)));
 
     const ModelVector mom = -mul(gradients(uu), value(uu)) - cs2 * gradient(lnrho) +
                             get(AC_nu_visc) * (laplace_vec(uu) +
@@ -623,8 +620,8 @@ induction(const ModelVectorData& uu, const ModelVectorData& aa)
 static inline ModelScalar
 lnT(const ModelScalarData& ss, const ModelScalarData& lnrho)
 {
-    const ModelScalar lnT = LNT0 + get(AC_gamma) * value(ss) / get(AC_cp_sound) +
-                            (get(AC_gamma) - ModelScalar(1.)) * (value(lnrho) - LNRHO0);
+    const ModelScalar lnT = get(AC_lnT0) + get(AC_gamma) * value(ss) / get(AC_cp_sound) +
+                            (get(AC_gamma) - ModelScalar(1.)) * (value(lnrho) - get(AC_lnrho0));
     return lnT;
 }
 
@@ -678,41 +675,96 @@ entropy(const ModelScalarData& ss, const ModelVectorData& uu, const ModelScalarD
     */
 }
 
-static bool
+static inline bool
 is_valid(const ModelScalar a)
 {
     return !isnan(a) && !isinf(a);
 }
 
-static bool
+static inline bool
 is_valid(const ModelVector& a)
 {
     return is_valid(a.x) && is_valid(a.y) && is_valid(a.z);
 }
 
-static inline ModelVector
-forcing(int3 globalVertexIdx)
-{
-    // source (origin)
-    ModelVector a = ModelScalar(.5) * (ModelVector){get(AC_nx) * get(AC_dsx),
-                                                    get(AC_ny) * get(AC_dsy),
-                                                    get(AC_nz) * get(AC_dsz)};
-    // sink (current index)
-    ModelVector b = (ModelVector){(globalVertexIdx.x - get(AC_nx_min)) * get(AC_dsx),
-                                  (globalVertexIdx.y - get(AC_ny_min)) * get(AC_dsy),
-                                  (globalVertexIdx.z - get(AC_nz_min)) * get(AC_dsz)};
+#if 0
+//FORCING NOT SUPPORTED FOR AUTOTEST
 
-    ModelScalar magnitude = 0.05;
-    // Vector c = magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
-    ModelVector c = magnitude * cross(normalized(b - a), (ModelVector){0, 0, 1}); // Vortex
-    if (is_valid(c)) {
-        return c;
-    }
-    else {
-        return (ModelVector){0, 0, 0};
-    }
+static inline ModelVector
+simple_vortex_forcing(ModelVector a, ModelVector b, ModelScalar magnitude)
+{
+    return magnitude * cross(normalized(b - a), (ModelVector){0, 0, 1}); // Vortex
 }
 
+static inline ModelVector
+simple_outward_flow_forcing(ModelVector a, ModelVector b, ModelScalar magnitude)
+{
+    return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
+}
+
+
+// The Pencil Code forcing_hel_noshear(), manual Eq. 222, inspired forcing function with adjustable helicity
+static inline ModelVector
+helical_forcing(ModelScalar magnitude, ModelVector k_force, ModelVector xx, ModelVector ff_re, ModelVector ff_im, ModelScalar phi)
+{
+
+    xx.x = xx.x*(2.0*M_PI/(get(AC_dsx)*(get(AC_ny_max) - get(AC_ny_min))));
+    xx.y = xx.y*(2.0*M_PI/(get(AC_dsy)*(get(AC_ny_max) - get(AC_ny_min))));
+    xx.z = xx.z*(2.0*M_PI/(get(AC_dsz)*(get(AC_ny_max) - get(AC_ny_min))));
+
+    ModelScalar cos_phi = cos(phi);
+    ModelScalar sin_phi = sin(phi);
+    ModelScalar cos_k_dox_x     = cos(dot(k_force, xx));
+    ModelScalar sin_k_dox_x     = sin(dot(k_force, xx));
+    // Phase affect only the x-component
+    ModelScalar real_comp_phase = cos_k_dox_x*cos_phi - sin_k_dox_x*sin_phi;
+    ModelScalar imag_comp_phase = cos_k_dox_x*sin_phi + sin_k_dox_x*cos_phi;
+
+
+    ModelVector force = (ModelVector){ ff_re.x*real_comp_phase - ff_im.x*imag_comp_phase,
+                             ff_re.y*real_comp_phase - ff_im.y*imag_comp_phase,
+                             ff_re.z*real_comp_phase - ff_im.z*imag_comp_phase};
+
+    return force;
+}
+
+static inline ModelVector
+forcing(int3 globalVertexIdx, ModelScalar dt)
+{
+    ModelVector a = ModelScalar(.5) * (ModelVector){ get(AC_nx) * get(AC_dsx),
+                                     get(AC_ny) * get(AC_dsy),
+                                     get(AC_nz) * get(AC_dsz)}; // source (origin)
+    ModelVector xx = (ModelVector){(globalVertexIdx.x - get(AC_nx_min)) * get(AC_dsx),
+                        (globalVertexIdx.y - get(AC_ny_min) * get(AC_dsy)),
+                        (globalVertexIdx.z - get(AC_nz_min) * get(AC_dsz))}; // sink (current index)
+    const ModelScalar cs2 = get(AC_cs2_sound);
+    const ModelScalar cs = sqrt(cs2);
+
+    //Placeholders until determined properly
+    ModelScalar magnitude = get(AC_forcing_magnitude);
+    ModelScalar phase     = get(AC_forcing_phase);
+    ModelVector k_force   = (ModelVector){  get(AC_k_forcex),   get(AC_k_forcey),   get(AC_k_forcez)};
+    ModelVector ff_re     = (ModelVector){get(AC_ff_hel_rex), get(AC_ff_hel_rey), get(AC_ff_hel_rez)};
+    ModelVector ff_im     = (ModelVector){get(AC_ff_hel_imx), get(AC_ff_hel_imy), get(AC_ff_hel_imz)};
+
+
+    //Determine that forcing funtion type at this point.
+    //ModelVector force = simple_vortex_forcing(a, xx, magnitude);
+    //ModelVector force = simple_outward_flow_forcing(a, xx, magnitude);
+    ModelVector force   = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase);
+
+    //Scaling N = magnitude*cs*sqrt(k*cs/dt)  * dt
+    const ModelScalar NN = cs*sqrt(get(AC_kaver)*cs);
+    //MV: Like in the Pencil Code. I don't understandf the logic here.
+    force.x = sqrt(dt)*NN*force.x;
+    force.y = sqrt(dt)*NN*force.y;
+    force.z = sqrt(dt)*NN*force.z;
+
+    if (is_valid(force)) { return force; }
+    else                 { return (ModelVector){0, 0, 0}; }
+}
+#endif
+
 static void
 solve_alpha_step(const int step_number, const ModelScalar dt, const int i, const int j, const int k,
                  ModelMesh& in, ModelMesh* out)
@@ -774,12 +826,13 @@ solve_beta_step(const int step_number, const ModelScalar dt, const int i, const
     for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w)
         out->vertex_buffer[w][idx] += beta[step_number] * in.vertex_buffer[w][idx];
 
+    (void)dt; // Suppress unused variable warning if forcing not used
 #if LFORCING
     if (step_number == 2) {
-        ModelVector force = forcing((int3){i, j, k});
-        out->vertex_buffer[VTXBUF_UUX][idx] += force.x * dt;
-        out->vertex_buffer[VTXBUF_UUY][idx] += force.y * dt;
-        out->vertex_buffer[VTXBUF_UUZ][idx] += force.z * dt;
+        ModelVector force = forcing((int3){i, j, k}, dt);
+        out->vertex_buffer[VTXBUF_UUX][idx] += force.x;
+        out->vertex_buffer[VTXBUF_UUY][idx] += force.y;
+        out->vertex_buffer[VTXBUF_UUZ][idx] += force.z;
     }
 #endif
 }
diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc
index 57cd0fe..156e23a 100644
--- a/src/standalone/simulation.cc
+++ b/src/standalone/simulation.cc
@@ -29,6 +29,7 @@
 #include "config_loader.h"
 #include "core/errchk.h"
 #include "core/math_utils.h"
+#include "model/host_forcing.h"
 #include "model/host_memory.h"
 #include "model/host_timestep.h"
 #include "model/model_reduce.h"
@@ -202,7 +203,7 @@ run_simulation(void)
     load_config(&mesh_info);
 
     AcMesh* mesh = acmesh_create(mesh_info);
-    //TODO: This need to be possible to define in astaroth.conf 
+    // TODO: This need to be possible to define in astaroth.conf
     acmesh_init_to(INIT_TYPE_GAUSSIAN_RADIAL_EXPL, mesh);
 
     acInit(mesh_info);
@@ -229,17 +230,54 @@ run_simulation(void)
     acStore(mesh);
     save_mesh(*mesh, 0, t_step);
 
-    const int max_steps  = mesh_info.int_params[AC_max_steps];
-    const int save_steps = mesh_info.int_params[AC_save_steps];
+    const int max_steps      = mesh_info.int_params[AC_max_steps];
+    const int save_steps     = mesh_info.int_params[AC_save_steps];
     const int bin_save_steps = mesh_info.int_params[AC_bin_steps]; // TODO Get from mesh_info
 
     AcReal bin_save_t = mesh_info.real_params[AC_bin_save_t];
     AcReal bin_crit_t = bin_save_t;
 
+    // Forcing properties
+    // AcReal relhel    = mesh_info.real_params[AC_relhel]; // JP: gcc warning: unused
+    // AcReal magnitude = mesh_info.real_params[AC_forcing_magnitude]; // JP: gcc warning: unused
+    // AcReal kmin = mesh_info.real_params[AC_kmin]; // JP: gcc warning: unused
+    // AcReal kmax = mesh_info.real_params[AC_kmax]; // JP: gcc warning: unused
+
+    // AcReal kaver = (kmax - kmin) / AcReal(2.0); // JP: gcc warning: unused
+
+    /* initialize random seed: */
+    srand(312256655);
+
     /* Step the simulation */
     for (int i = 1; i < max_steps; ++i) {
         const AcReal umax = acReduceVec(RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
         const AcReal dt   = host_timestep(umax, mesh_info);
+
+#if LFORCING
+        // Generate a forcing vector before canculating an integration step.
+
+        // Generate forcing wave vector k_force
+        AcReal3 k_force;
+        k_force = helical_forcing_k_generator(kmax, kmin);
+
+        // Randomize the phase
+        AcReal phase = AcReal(2.0) * AcReal(M_PI) * get_random_number_01();
+
+        // Generate e for k. Needed for the sake of isotrophy.
+        AcReal3 e_force;
+        if ((k_force.y == 0.0) && (k_force.z == 0.0)) {
+            e_force = (AcReal3){0.0, 1.0, 0.0};
+        }
+        else {
+            e_force = (AcReal3){1.0, 0.0, 0.0};
+        }
+        helical_forcing_e_generator(&e_force, k_force);
+
+        AcReal3 ff_hel_re, ff_hel_im;
+        helical_forcing_special_vector(&ff_hel_re, &ff_hel_im, k_force, e_force, relhel);
+        acForcingVec(magnitude, k_force, ff_hel_re, ff_hel_im, phase, kaver);
+#endif
+
         acIntegrate(dt);
 
         t_step += dt;