diff --git a/samples/standalone/model/host_timestep.cc b/samples/standalone/model/host_timestep.cc
index fb1b52c..897adb2 100644
--- a/samples/standalone/model/host_timestep.cc
+++ b/samples/standalone/model/host_timestep.cc
@@ -31,7 +31,7 @@
 static AcReal timescale = AcReal(1.0);
 
 AcReal
-host_timestep(const AcReal& umax, const AcMeshInfo& mesh_info)
+host_timestep(const AcReal& umax, const AcReal& vAmax, const AcMeshInfo& mesh_info)
 {
     const long double cdt  = mesh_info.real_params[AC_cdt];
     const long double cdtv = mesh_info.real_params[AC_cdtv];
@@ -40,7 +40,7 @@ host_timestep(const AcReal& umax, const AcMeshInfo& mesh_info)
     const long double nu_visc   = mesh_info.real_params[AC_nu_visc];
     const long double eta       = mesh_info.real_params[AC_eta];
     const long double chi       = 0; // mesh_info.real_params[AC_chi]; // TODO not calculated
-    const long double gamma     = mesh_info.real_params[AC_gamma];
+    const long double gamma     = mesh_info.real_params[AC_gamma]; //TODO this does not make sense here at all. 
     const long double dsmin     = mesh_info.real_params[AC_dsmin];
 
     // Old ones from legacy Astaroth
@@ -49,12 +49,11 @@ 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 uu_dt   = cdt * dsmin / (fabsl(umax) + sqrtl(cs2_sound + 0.0l));
+    const long double uu_dt   = cdt * dsmin / (fabsl(umax) + sqrtl(cs2_sound + vAmax*vAmax));
     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!
+                                    max(gamma, chi));
 
     const long double dt = min(uu_dt, visc_dt);
     return AcReal(timescale) * AcReal(dt);
diff --git a/samples/standalone/model/host_timestep.h b/samples/standalone/model/host_timestep.h
index dd5cdf7..9074b26 100644
--- a/samples/standalone/model/host_timestep.h
+++ b/samples/standalone/model/host_timestep.h
@@ -27,6 +27,6 @@
 #pragma once
 #include "astaroth.h"
 
-AcReal host_timestep(const AcReal& umax, const AcMeshInfo& mesh_info);
+AcReal host_timestep(const AcReal& umax, const AcReal& vAmax, const AcMeshInfo& mesh_info);
 
 void set_timescale(const AcReal scale);