diff --git a/examples/BigCosmoVolume/makeIC.py b/examples/BigCosmoVolume/makeIC.py
index 0994e1c95e053defe7766122c52bc405c7239776..3020feaf753f817f039d2fd09c4fa4f7fb69b896 100644
--- a/examples/BigCosmoVolume/makeIC.py
+++ b/examples/BigCosmoVolume/makeIC.py
@@ -77,7 +77,7 @@ indices = indices < numPart
 coords = coords[indices,:]
 v = v[indices,:]
 m = m[indices]
-h = h[indices]
+h = h[indices] / 1.825742 # Correct from Gadget defintion of h to physical definition
 u = u[indices]
 ids = ids[indices]
 
diff --git a/examples/GreshoVortex/makeIC.py b/examples/GreshoVortex/makeIC.py
index 6aceeed559324f97b0b1e388ff0c3524498b52e4..12edcb6e8154ec6f865d28a6daeb02d385d14bbf 100644
--- a/examples/GreshoVortex/makeIC.py
+++ b/examples/GreshoVortex/makeIC.py
@@ -30,6 +30,7 @@ factor = 3
 boxSize = [ 1.0 , 1.0, 1.0/factor ]
 L = 120           # Number of particles along one axis
 gamma = 5./3.     # Gas adiabatic index
+eta = 1.2349      # 48 ngbs with cubic spline kernel
 rho = 1           # Gas density
 P0 = 0.           # Constant additional pressure (should have no impact on the dynamics)
 fileName = "greshoVortex.hdf5" 
@@ -73,7 +74,7 @@ for i in range(L):
             v[index,1] =  v_phi * (x - boxSize[0] / 2) / r
             v[index,2] = 0.
             m[index] = mass
-            h[index] = 2.251 * boxSize[0] / L
+            h[index] = eta * boxSize[0] / L
             P = P0
             if r < 0.2:
                 P = P + 5. + 12.5*r2
@@ -105,6 +106,14 @@ grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 #Particle group
 grp = file.create_group("/PartType0")
 ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
diff --git a/examples/MultiTypes/makeIC.py b/examples/MultiTypes/makeIC.py
index 3a41910c22c260086b5384b248a5c86ab6340a5e..cf889f9b6eab502f692cd6c8b4506c31664ecdcb 100644
--- a/examples/MultiTypes/makeIC.py
+++ b/examples/MultiTypes/makeIC.py
@@ -32,6 +32,7 @@ Lgas = int(sys.argv[1])  # Number of particles along one axis
 rhoGas = 2.              # Density
 P = 1.                   # Pressure
 gamma = 5./3.            # Gas adiabatic index
+eta = 1.2349             # 48 ngbs with cubic spline kernel
 rhoDM = 1.
 Ldm = int(sys.argv[2])  # Number of particles along one axis
 
@@ -61,11 +62,18 @@ grp.attrs["NumFilesPerSnapshot"] = 1
 grp.attrs["MassTable"] = [0.0, massDM, 0.0, 0.0, 0.0, 0.0]
 grp.attrs["Flag_Entropy_ICs"] = 0
 
-
 #Runtime parameters
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 
 # Gas Particle group
 grp = file.create_group("/PartType0")
@@ -80,7 +88,7 @@ ds = grp.create_dataset('Masses', (numGas,1), 'f')
 ds[()] = m
 m = zeros(1)
 
-h = full((numGas, 1), 1.1255 * boxSize / Lgas)
+h = full((numGas, 1), eta * boxSize / Lgas)
 ds = grp.create_dataset('SmoothingLength', (numGas,1), 'f')
 ds[()] = h
 h = zeros(1)
diff --git a/examples/PerturbedBox/makeIC.py b/examples/PerturbedBox/makeIC.py
index 69c1a69199c9a5262f5ae6c4e95ca14699300fd4..ee1d845fc2149892909a54bf588046b0b1691b03 100644
--- a/examples/PerturbedBox/makeIC.py
+++ b/examples/PerturbedBox/makeIC.py
@@ -90,6 +90,14 @@ grp.attrs["NumPart_Total"] = numPart
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 #Particle group
 grp = file.create_group("/PartType0")
 ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
diff --git a/examples/SedovBlast/makeIC.py b/examples/SedovBlast/makeIC.py
index 75ff81165df51780848e3d8ac679a6dbeb17a039..e64942e8e92ee6fe67142f841f566019b1a668be 100644
--- a/examples/SedovBlast/makeIC.py
+++ b/examples/SedovBlast/makeIC.py
@@ -33,6 +33,7 @@ P = 1.e-5         # Pressure
 E0= 1.e2          # Energy of the explosion
 pert = 0.1
 gamma = 5./3.     # Gas adiabatic index
+eta = 1.2349      # 48 ngbs with cubic spline kernel
 fileName = "sedov.hdf5" 
 
 
@@ -67,7 +68,7 @@ for i in range(L):
             v[index,1] = 0.
             v[index,2] = 0.
             m[index] = mass
-            h[index] = 1.1255 * boxSize / L
+            h[index] = eta * boxSize / L
             u[index] = internalEnergy
             ids[index] = index + 1
             if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 2.01 * boxSize/L:
@@ -98,6 +99,14 @@ grp.attrs["Flag_Entropy_ICs"] = 0
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 #Particle group
 grp = file.create_group("/PartType0")
 grp.create_dataset('Coordinates', data=coords, dtype='d')
diff --git a/examples/SedovBlast/makeIC_fcc.py b/examples/SedovBlast/makeIC_fcc.py
index 17f07440909cb5478d09a5b7a1444c72af2f3a47..0d3a017a9b7f3b30b61e723e3d1646d7797b40a4 100644
--- a/examples/SedovBlast/makeIC_fcc.py
+++ b/examples/SedovBlast/makeIC_fcc.py
@@ -33,6 +33,7 @@ P = 1.e-5         # Pressure
 E0= 1.e2          # Energy of the explosion
 pert = 0.025
 gamma = 5./3.     # Gas adiabatic index
+eta = 1.2349          # 48 ngbs with cubic spline kernel
 fileName = "sedov.hdf5" 
 
 
@@ -70,7 +71,7 @@ for i in range(L):
                 v[index,1] = 0.
                 v[index,2] = 0.
                 m[index] = mass
-                h[index] = 1.1255 * hbox
+                h[index] = eta * hbox
                 u[index] = internalEnergy
                 ids[index] = index + 1
                 if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 1.2 * hbox:
@@ -101,6 +102,14 @@ grp.attrs["Flag_Entropy_ICs"] = 0
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 #Particle group
 grp = file.create_group("/PartType0")
 grp.create_dataset('Coordinates', data=coords, dtype='d')
diff --git a/examples/SodShock/makeIC.py b/examples/SodShock/makeIC.py
index 0ac0564116f8a6ceb57b4f41d23eb9907df0440d..8ae19050c11c0712579b44646c8870d7574d113b 100644
--- a/examples/SodShock/makeIC.py
+++ b/examples/SodShock/makeIC.py
@@ -43,14 +43,14 @@ vol = boxSize[0] * boxSize[1] * boxSize[2]
 glass1 = h5py.File("glass_001.hdf5")
 pos1 = glass1["/PartType0/Coordinates"][:,:]
 pos1 = pos1 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
-
+glass_h1 = glass1["/PartType0/SmoothingLength"][:] / factor
 
 #Read in high density glass
 # glass2 = h5py.File("../Glass/glass_50000.hdf5")
 glass2 = h5py.File("glass_002.hdf5")
 pos2 = glass2["/PartType0/Coordinates"][:,:]
 pos2 = pos2 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
-
+glass_h2 = glass2["/PartType0/SmoothingLength"][:] / factor
 
 #Generate high density region
 rho1 = 1.
@@ -61,9 +61,10 @@ coord1 = append(coord1, coord1 + [0.25, 0, 0], 0)
 # coord1 = append(coord1, pos1 + [0, 0.5, 0.5], 0)
 N1 = size(coord1)/3
 v1 = zeros((N1, 3))
-h1 = ones(N1) * 2.251 * 0.5 * vol / (size(pos1)/3)**(1./3.)
 u1 = ones(N1) * P1 / ((gamma - 1.) * rho1)
 m1 = ones(N1) * vol * 0.5 * rho1 / N1
+h1 = append(glass_h1, glass_h1, 0)
+h1 = append(h1, h1, 0)
 
 #Generate low density region
 rho2 = 0.25
@@ -74,9 +75,10 @@ coord2 = append(coord2, coord2 + [0.25, 0, 0], 0)
 # coord2 = append(coord2, pos2 + [0, 0.5, 0.5], 0)
 N2 = size(coord2)/3
 v2 = zeros((N2, 3))
-h2 = ones(N2) * 2.251 * 0.5 * vol / (size(pos2)/3)**(1./3.)
 u2 = ones(N2) * P2 / ((gamma - 1.) * rho2)
 m2 = ones(N2) * vol * 0.5 * rho2 / N2
+h2 = append(glass_h2, glass_h2, 0)
+h2 = append(h2, h2, 0)
 
 #Merge arrays
 numPart = N1 + N2
@@ -89,8 +91,8 @@ ids = zeros(numPart, dtype='L')
 for i in range(1, numPart+1):
     ids[i-1] = i
 
-#Final operations
-h /= 2
+#Final operation since we come from Gadget-2 cubic spline ICs
+h /= 1.825752
 
 #File
 file = h5py.File(fileName, 'w')
@@ -110,6 +112,14 @@ grp.attrs["Flag_Entropy_ICs"] = 0
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 #Particle group
 grp = file.create_group("/PartType0")
 grp.create_dataset('Coordinates', data=coords, dtype='d')
diff --git a/examples/UniformBox/makeIC.py b/examples/UniformBox/makeIC.py
index c175349e658799cbcb30dfe2619a1594bafc18b9..1484f60596e68734f0f98685ab2ab845f2e0b407 100644
--- a/examples/UniformBox/makeIC.py
+++ b/examples/UniformBox/makeIC.py
@@ -32,6 +32,7 @@ L = int(sys.argv[1])  # Number of particles along one axis
 rho = 2.              # Density
 P = 1.                # Pressure
 gamma = 5./3.         # Gas adiabatic index
+eta = 1.2349          # 48 ngbs with cubic spline kernel
 fileName = "uniformBox.hdf5" 
 
 #---------------------------------------------------
@@ -55,11 +56,18 @@ grp.attrs["NumFilesPerSnapshot"] = 1
 grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 grp.attrs["Flag_Entropy_ICs"] = 0
 
-
 #Runtime parameters
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 #Particle group
 grp = file.create_group("/PartType0")
 
@@ -73,7 +81,7 @@ ds = grp.create_dataset('Masses', (numPart,1), 'f')
 ds[()] = m
 m = zeros(1)
 
-h = full((numPart, 1), 1.1255 * boxSize / L)
+h = full((numPart, 1), eta * boxSize / L)
 ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
 ds[()] = h
 h = zeros(1)
diff --git a/examples/UniformBox/makeICbig.py b/examples/UniformBox/makeICbig.py
index e475fdcbd9f3c4811e3dcfdf20bbd321be3d8b29..bd5cf627fb535595b3abb224cbc8de50589f12cf 100644
--- a/examples/UniformBox/makeICbig.py
+++ b/examples/UniformBox/makeICbig.py
@@ -32,6 +32,7 @@ N = int(sys.argv[2])  # Write N particles at a time to avoid requiring a lot of
 rho = 2.              # Density
 P = 1.                # Pressure
 gamma = 5./3.         # Gas adiabatic index
+eta = 1.2349      # 48 ngbs with cubic spline kernel
 fileName = "uniformBox_%d.hdf5"%L
 
 #---------------------------------------------------
@@ -62,11 +63,19 @@ grp.attrs["NumFilesPerSnapshot"] = 1
 grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 grp.attrs["Flag_Entropy_ICs"] = 0
 
-
 #Runtime parameters
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+
 #Particle group
 grp = file.create_group("/PartType0")
 
@@ -89,7 +98,7 @@ for n in range(n_iterations):
     ds_m[offset:offset+N] = m
     m = zeros(1)
 
-    h = full((N, 1), 1.1255 * boxSize / L)
+    h = full((N, 1), eta * boxSize / L)
     ds_h[offset:offset+N] = h
     h = zeros(1)
 
@@ -122,7 +131,7 @@ m = full((remainder, 1), mass)
 ds_m[offset:offset+remainder] = m
 m = zeros(1)
 
-h = full((remainder, 1), 1.1255 * boxSize / L)
+h = full((remainder, 1), eta * boxSize / L)
 ds_h[offset:offset+remainder] = h
 h = zeros(1)
 
@@ -139,7 +148,7 @@ coords = zeros((remainder, 3))
 coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
 coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
 coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
-ds_x[offset:offset+remainder,:] = coords
+ods_x[offset:offset+remainder,:] = coords
 
 print "Done", offset+remainder,"/", numPart
 
diff --git a/examples/UniformDMBox/makeIC.py b/examples/UniformDMBox/makeIC.py
index 061b4d0ad1959d9e25356aff80e78adb9c1c4faa..449d780fb31bc23dd194f772be45d35e6b0bbe3f 100644
--- a/examples/UniformDMBox/makeIC.py
+++ b/examples/UniformDMBox/makeIC.py
@@ -52,11 +52,19 @@ grp.attrs["NumFilesPerSnapshot"] = 1
 grp.attrs["MassTable"] = [0.0, mass, 0.0, 0.0, 0.0, 0.0]
 grp.attrs["Flag_Entropy_ICs"] = 0
 
-
 #Runtime parameters
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+
 #Particle group
 grp = file.create_group("/PartType1")
 
diff --git a/src/debug.c b/src/debug.c
index 4c1434118c98aab7def28d3a53493767d249d774..53a03d66aee2c169a555ed00a2efa2d5b984066a 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -60,7 +60,7 @@ void printParticle(struct part *parts, struct xpart *xparts, long long int id,
   /* Look for the particle. */
   for (size_t i = 0; i < N; i++)
     if (parts[i].id == id) {
-      printf("## Particle[%zd]:\n id=%lld", i, parts[i].id);
+      printf("## Particle[%zd]:\n id=%lld ", i, parts[i].id);
       hydro_debug_particle(&parts[i], &xparts[i]);
       found = 1;
       break;
@@ -76,12 +76,12 @@ void printgParticle(struct gpart *gparts, long long int id, size_t N) {
   /* Look for the particle. */
   for (size_t i = 0; i < N; i++)
     if (gparts[i].id == -id) {
-      printf("## gParticle[%zd] (DM) :\n id=%lld", i, -gparts[i].id);
+      printf("## gParticle[%zd] (DM) :\n id=%lld ", i, -gparts[i].id);
       gravity_debug_particle(&gparts[i]);
       found = 1;
       break;
     } else if (gparts[i].id > 0 && gparts[i].part->id == id) {
-      printf("## gParticle[%zd] (hydro) :\n id=%lld", i, gparts[i].id);
+      printf("## gParticle[%zd] (hydro) :\n id=%lld ", i, gparts[i].id);
       gravity_debug_particle(&gparts[i]);
       found = 1;
       break;
diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h
index 46e156bb99015069f9958aeea05954e2be6db5e0..a4d1f7dd4397ebfc850b582e1ca81fc0d4edb76a 100644
--- a/src/hydro/Gadget2/hydro_debug.h
+++ b/src/hydro/Gadget2/hydro_debug.h
@@ -23,13 +23,13 @@ __attribute__((always_inline))
       "x=[%.3e,%.3e,%.3e], "
       "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
       "h=%.3e, "
-      "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, "
+      "wcount=%d, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, "
       "dS/dt=%.3e, c=%.3e\n"
       "divV=%.3e, curlV=%.3e, rotV=[%.3e,%.3e,%.3e]  \n "
       "v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
       xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
-      p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho,
+      p->h, (int)p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh, p->rho,
       p->force.pressure, p->entropy, p->entropy_dt, p->force.soundspeed,
       p->div_v, p->force.curl_v, p->density.rot_v[0], p->density.rot_v[1],
       p->density.rot_v[2], p->force.v_sig, p->h_dt, p->ti_begin, p->ti_end);